JP7630781B2 - How to extract a respiratory curve - Google Patents

Description

The present invention relates to a method for extracting a respiratory curve and a method for extracting abnormal respiratory events. The method extracts a respiratory curve, and in particular extracts apnea and other abnormal phenomena in the respiratory system.

Previously, there was a method for extracting a respiratory curve. In Patent Document 1, the curve was calculated using the reciprocal of the standard deviation of the respiratory frequency.

Patent No. 5670071

However, the method described in Patent Document 1 was unable to accurately extract the respiratory curve.

Therefore, the objective of this application is to provide a method for extracting a respiratory curve and a method for extracting abnormal respiratory events that can accurately extract a respiratory curve.

In order to solve the above problem, a respiration curve extraction method is used that includes a conversion process of measuring a two-dimensional human body vibration waveform and converting the measured first data into frequency domain data, a GMCA process of performing generalized morphological component analysis (GMCA) on the second data converted in the conversion process based on the sparsity of the signal representation by morphological component analysis to obtain the second data, and an inverse conversion process of converting the third data obtained in the GMCA process into discrete signal data.

The method for detecting abnormal breathing events includes a division step of dividing a breathing curve into a number of epochs, a classification step of classifying the epochs into a number of clusters, and an abnormal breathing section extraction step of extracting clusters of abnormal breathing sections.

The respiratory curve extraction method and abnormal respiratory event detection method of the present invention enable accurate extraction of the respiratory curve and extraction of abnormal respiratory events.

FIG. 1 is a diagram showing an outline of a respiratory motion extraction method according to an embodiment. FIG. 2 is a plan view illustrating pressure measurement according to an embodiment; FIG. 13 is a diagram showing the average value of correlation coefficients when initial values are set, which are considered in the BSS process according to the embodiment. FIG. 13 is a diagram showing the distribution of correlation coefficients under each initial value condition of the BSS process according to the embodiment; A graph showing the average correlation coefficients between the respiratory curves and RIPsum obtained by the Nishida, ICA, and GMCA methods. Distribution of correlation coefficients for each method A graph showing the average correlation coefficient for each subject in the results of respiratory curve extraction using GMCA with the initial condition Initial 10. A schematic diagram of the singular spectrum transform FIG. 1 is an example of the results of calculating RIPsum from the RIP signal measured in a PSG test and dividing the respiratory curve using SST. A dendrogram showing the hierarchical structure of the connections of partial respiration curves. (a) to (d) are diagrams showing the results of dividing and clustering a respiratory curve.

(First embodiment) Respiration curve extraction using generalized morphological component analysis <Outline of respiration curve extraction method>
First, the analysis of a person's pressure signal (observed signal) and obtaining a respiration curve will be described.
<Respiratory movement extraction method>
FIG. 1 shows an outline of a respiratory motion extraction method according to an embodiment.

1-1: Measure pressure signal Bias is removed from the pressure time series at each measurement point by subtracting the epoch average pressure value at the measurement point. An observed signal matrix X is created from the pressure time series with the bias removed. The measured body pressure distribution time series is processed by the following GMCA (explained below) for each sleep epoch (30 seconds).

1-2: Data conversion Perform DCT on X. The coefficient matrix obtained by performing DCT on X is ΘX.
1-3: BSS processing (GMCA)
By applying GMCA to ΘX, the estimated values ~A and ~ΘS of the mixing matrix and the original signal matrix are obtained.

1-4: Inverse data transformation Then, an inverse DCT is applied to ~ΘS. An estimate of the source signal matrix &#732;S is obtained. ~S contains n source signal estimates ~sj. The column vector ~aj of ~A represents how much of the corresponding source signal ~sj is contained in the signal at each observation point. Below, the elements ~aji of ~aj are called weighting coefficients. The source signal that is determined to be closest to respiratory motion among the source signals is the respiratory curve estimated by GMCA. The weighting coefficient vector for this respiratory curve represents the spatial distribution of respiratory motion on the surface pressure sensor.

<Details of the method for extracting respiration curves>
<1-1: Measuring pressure signal>
This will be explained with reference to Fig. 2. In Fig. 2, a person 10 is lying on a sensor 11. A control unit 12 can measure pressure at each point of the grid.

<1: Measuring pressure signals and information contained in the time series of body pressure distribution>
The areal pressure sensor is a device that measures pressure at multiple measurement points and measures two-dimensional pressure distribution due to body weight. The time series of body pressure distribution measured by this device includes periodic pressure changes caused by respiratory effort and heart rate.

The signal σ measured by a pressure sensor at one point can be modeled as follows:
σ(t)=SH(t)+SR(t)+M(t)+B(t)+N(t)
Here, SH and SR represent signals due to breathing and heart rate, respectively, M represents a signal due to body movements such as turning over in bed, B represents a bias due to weight, and N represents noise.
If it were possible to properly separate and extract information on respiratory and cardiac movements from the measured body pressure distribution time series σ(t) using a surface pressure sensor, this would be useful in estimating the subject's sleep state.

<1-2 Data Conversion>
Discrete Cosine Transform (DCT) is a method for transforming a discrete signal into the frequency domain. For a discrete signal X0,...,XN-1 of N samples, when transformed into N coefficients C0,...,CN-1, the DCT is defined as follows:
...Number 1
...Number 2

The discrete cosine transform is one method for transforming a discrete signal into the frequency domain, although other transform methods are also possible.
Similarly, the inverse discrete cosine transform, which will be described later, is one method of transforming the frequency domain into a discrete signal, and other transform methods may also be used.
The embodiment will be described below, but it is merely an example of the invention and the invention is not limited thereto.
<1-3: BSS Processing (GMCA)>: Generalized Morphological Component Analysis When signals output from n signal sources are measured by m sensors to obtain observed signals, the problem of separating and estimating the unknown original signal, i.e., the source signal, from the observed signals is called blind source separation (BSS).

Let x i (i=1, . . . , m) be an observation signal of length N measured at m observation points.
The matrix in which the observed signals are arranged is denoted as X:=[xT1, . . . , xTm]T.
On the other hand, let n source signals be sj (j=1, . . . , n).
Let S:=[sT1,...,sTn]T be the matrix in which the source signals are arranged. If A is the mxn mixing matrix, the relationship between the observed signal matrix X and the source signal matrix S is expressed as X=AS. Since X is known, A and S are unknown and generally cannot be solved, BSS achieves source separation based on some assumptions. Examples of BSS techniques include principal component analysis, independent component analysis, and nonnegative matrix factorization.

Morphological Component Analysis (MCA) (J.-L Starck, Y.A, J. Bobin, Michael Elad, and D. Donoho. Morphological component analysis. Proceedings of SPIE - The International Society for Optical Engineering, Vol. 5914,, 08 2005.) has been proposed as a method for signal decomposition and noise removal, assuming that signals and images can be represented using some of the bases when a dictionary of redundant signal representation bases is given. MCA assumes that a signal y is sparsely represented (the coefficients for most of the bases are zero) in a specific basis Φk in a dictionary D, and can be expressed as a linear combination of D morphological components. Here, the dictionary D is a combination of orthonormal bases Φk such that D = [ΦT1,...,ΦTD]T. In this case, y is expressed as Equation 3.
...Number 3

where αk is a coefficient vector that associates y with Φk.
Since the coefficient vector α is sparse, the minimization problem of the following equation 4 is obtained.
...Number 4
∥... . . ∥0 is the l0 norm and represents the number of non-zero components. GMCA was proposed as a method to solve this minimization problem (generalized MCA). Here, if there exists a unique sparse decomposition of signal y and y=α*D is the solution, then α*＝ΔD(y).

GMCA assumes that the source signal is sparsely represented in a specific basis Φk in a dictionary D and can be expressed as a linear combination of D morphological components. GMCA attempts BSS by relying on the sparsity of the signal representation.
Assuming that each source signal sj is sparsely represented by a linear combination of D morphological components using a specific basis Φk in the dictionary D, the source signal sj is expressed as follows:
...Number 5

Here, αjk is a coefficient vector that associates sj with Φk.
Since the source signal sj is sparsely represented using the dictionary D, the following optimization problem of Equation 6 can be defined for the decomposition of the observed signal matrix X.
...Number 6

By relaxing the l0 norm (∥ ⋅ ∥ 0) to the l1 norm and introducing the Lagrange undetermined multiplier λ j , the problem of Equation 6 can be approximated by the following Equation 7.
...Number 7

Here, ∥·∥ ² ₂ is the Frobenius norm, and ∥P∥ ² ₂ = Trace(P ^T P).
In GMCA, the optimization problem of Equation 7 is solved using an iterative reduced threshold algorithm, and estimates of the source signal matrix and the mixing matrix, ~S and ~A, are obtained from the observed signal matrix X. In other words, by gradually decreasing λj from a large value, the optimization calculation is repeated while relaxing the constraint on sparsity, and a solution to Equation 7 is obtained.

In the algorithm for solving Equation 7, when estimating the source signal matrix S by iterative calculation, it is necessary to perform a conversion by D at each calculation, which results in a high calculation cost. The fast algorithm of GMCA assumes that there is a sparse decomposition unique to the observed signal xi, and performs source separation in the coefficient space. The coefficient matrix of the sparse decomposition of the observed signal matrix xi is defined as ΘX:=[ΔD(x1)T,...,ΔD(xm)T]T, and the coefficient matrix of the sparse decomposition of the source signal matrix si is defined as ΘS:=[ΔD(s1)T,...,ΔD(sn)T]T. The fast algorithm approximates Equation 7 with the following optimization problem.
^{...Number 8}

Unless the dictionary D is an orthonormal basis, the formulas 7 and 8 are not equivalent. When the dictionary D is redundant, no rigorous mathematical proof is given, but good results are obtained by experiment.
When A is fixed in Equation 8, .THETA.S can be obtained by the order.
...Number 9

Here, A† is a pseudo-inverse matrix. Also, Sλ is a soft-thresholding operator that acts as follows when y[t] is input.
...Number 10
When ΘS is fixed in Equation 8, A becomes a least squares estimate and can be found by degree 11.

...Number 11
The iterative reduced threshold algorithm for solving Equation 8 is as follows:
1. Apply the discrete cosine transform to the observed signal x i to calculate ΘX.
2. Set the initial threshold λ0, the minimum threshold λmin, and the initial value Aini of the mixing matrix.
3.1 Estimate ΘS using Equation 9.
3.2 Estimate A using Equation 11.
3.3 Decrease the threshold λ.
3.4 If λ>λmin, go back to 3.1.

4. When the threshold value λ=λmin is reached, the calculation is terminated.
The smallest possible threshold λmin determines the sparsity of the final source signal. If the threshold is large, the source signal estimate is expressed in a smaller number of bases, and if the threshold is small, it is expressed in a larger number of bases. When the threshold λ is large, the source signal is estimated from the most important coefficients, and as the threshold λ becomes smaller, the error between the product of the mixing matrix and the source signal matrix and the observed signal matrix becomes smaller. This iterative reduced threshold algorithm ensures robustness by first estimating the most important features of the data and gradually estimating the details.

1-3-A Initial Value Setting of Generalized Morphological Component Analysis As shown in Equation 1, the time series of body pressure distribution mainly includes pressure changes in the respiratory frequency band, pressure changes in the heart rate frequency band, and pressure changes due to body movements such as turning over in bed. By providing this prior knowledge to GMCA, the accuracy of extracting respiratory motion is improved.
As prior knowledge to be given, information on the respiratory frequency band is obtained in the mixing matrix Aini.
The following four conditions are considered as initial values to be given to GMCA.

1. Random: A _ini is assigned a random value that follows the standard normal distribution.
2. PCA-based: Eigenvectors vj (j=1, 2, . . . , n) _of the variance-covariance matrix Σ ⁼ _ΘXΘTX of the coefficient matrix of the observed signal matrix are obtained by principal component analysis (PCA). _Aini is given by n eigenvectors vj (j=1, 2, . . . , m).

3. Initial1: For the coefficient vector Θs1 of the first source signal, the coefficient corresponding to the respiratory frequency band is set to 1, and the other coefficients are set to 0. For the weighting coefficient vector, the weighting coefficient vector a1 corresponding to the coefficient vector Θs1 of the first source signal is given as a1 = ΘXΘ ^T _s1 . The remaining weighting coefficient vectors aj (j = 1, 2, ..., n) are given random values following the standard normal distribution.

4. Initial10: As in Initial1, for the coefficient vector Θs1 of the first source signal, the coefficient corresponding to the respiratory frequency band is set to 1, and the other coefficients are set to 0.
Furthermore, for the coefficient vector Θs2 of the second source signal, the coefficient corresponding to the heart rate frequency band is set to 1, and the other coefficients are set to 0.

For the coefficient vector Θs3 of the third source signal, the coefficients corresponding to the frequency band lower than the respiratory frequency band are set to 1, and the other coefficients are set to 0. This corresponds to a body movement component such as turning over in bed.
For the coefficient vector Θs4 of the fourth source signal, the coefficient corresponding to the frequency band between the respiratory frequency band and the heart rate frequency band is set to 1, and the other coefficients are set to 0.

The coefficient vector of the remaining source signal is divided into frequency bands higher than the heart rate frequency band, and the coefficients corresponding to the divided frequency bands are set to 1 and the other coefficients are set to 0.
The initial value of the mixing matrix is given as Aini=ΘXΘ†S using the pseudo-inverse matrix Θ†s of the coefficient matrix ΘS of the source signal matrix defined by the above settings and the coefficient vector ΘX of the observed signal matrix.

<1-4: Reverse data conversion>
The inverse transform of DCT is called the inverse discrete cosine transform (Inverse DCT; IDCT), and is defined as follows:
...Number 12
The coefficients obtained by DCT represent the weights of the cosine functions of each frequency when a discrete signal is expressed as the sum of cosine functions of various frequencies.

<Example>
Extraction of respiration curves using generalized morphological component analysis 1: Data set A SmartRubber (registered trademark) sensor 11 manufactured by Sumitomo Riko Co., Ltd. (hereinafter, SR sensor 11) was used to measure the time series of body pressure distribution. The SR sensor 11 used is a mat-shaped pressure sensor with dimensions of 105 cm x 90 cm, and has 32 x 25 = 800 measurement points. The measurement cycle of each measurement point is 20 Hz. Pressure values are output as integers in the range of approximately 100 to 2000. The data set to be analyzed was measured in parallel with PSG examinations at the Ohta Sleep Science Center.

A PSG test is a type of test called polysomnography (PSG) that is used to confirm sleep apnea syndrome (SAS). This is a test that requires an overnight stay at a medical institution, where electrodes are attached with tape to the necessary parts of the head, face, and body, and the patient actually sleeps overnight while brain waves, breathing, eye and muscle movements are recorded to examine the state of sleep. The results are analyzed and compared with the progression of normal sleep to diagnose sleep apnea syndrome (SAS).

For the measurement, the SR sensor 11 was placed on the mattress of the bed on which the subject (person 10) slept in the examination room. In the PSG examination, a respiratory inductance plethysmograph (RIP) belt was attached to the subject's chest and abdomen to measure respiratory movement. The RIPsum signal rsum, which is the sum of the output potentials of the RIP belts on the chest and abdomen, was used as a reference signal to be compared with the processing results of the body pressure distribution time series. The correlation coefficient with RIPsum was used as the evaluation index for the respiratory movement extraction results. Since there are multiple source signal estimates ~sj output by the proposed method, the correlation coefficient was calculated by allowing a time shift of ±0.5 seconds, and the source signal with the maximum absolute value of the correlation coefficient was used as the respiratory curve, and the absolute value of the correlation coefficient at that time was used as the evaluation index.

2: Extraction accuracy of respiratory curves The respiratory curve extraction method using GMCA was applied to the time series of body pressure distribution during sleep for 23 subjects (10 people) to evaluate the extraction accuracy of the respiratory curves. The number of source signals was set to n = 10. A total of 22,158 epochs were analyzed, and the respiratory curve was extracted from the time series of body pressure distribution during sleep in each epoch and the correlation coefficient was calculated.

Figure 3 shows the average correlation coefficients when the initial values were set as discussed in 1-3-A. The average correlation coefficients were highest for Initial10, followed by Initial1, Random, and PCA-based. Since the average correlation coefficients for Initial10 and Initial1 were high, it can be seen that the accuracy of extracting respiratory movement is improved by using initial values based on the prior knowledge that the body pressure distribution time series contains signals in the respiratory frequency band.

Figure 4 shows the distribution of correlation coefficients for each initial value condition. For the initial value conditions of Initial10, Initial1, and Random, there are two peaks near correlation coefficients of 0.1 and 0.85. Compared to the distributions for the other conditions, the distribution for Initial10 has a smaller distribution near the peak with a correlation coefficient of 0.1 and a larger distribution near the peak with a correlation coefficient of 0.85. There is no significant difference in the position of the two peaks for the three initial value conditions. In other words, the improvement in the accuracy of extraction of respiratory movement with Initial10 compared to the other initial value conditions does not indicate that a respiratory curve closer to RIPsum is extracted, but that it makes it possible to extract a respiratory curve in epochs where a respiratory curve could not be extracted with the other initial value conditions.

Next, the accuracy of respiratory curve extraction was evaluated by applying the method of Nishida et al. (Nishida Yoshifumi, Takeda Masasuke, Mori Taketoshi, Mizo Hiroshi, Sato Tomomasa. Non-invasive and unrestricted measurement of breathing and body position during sleep using pressure sensors. Journal of the Robotics Society of Japan, Vol. 16, No. 5, pp. 705-711, July 1998.: Method for calculating the respiratory curve from the phase difference of signals from multiple sensors), independent component analysis (ICA), and GMCA-based respiratory curve extraction methods to the time series of body pressure distribution during sleep of the same 23 subjects. Initial 10 was used as the initial value when extracting the respiratory curve using GMCA. For ICA, the source signal obtained by performing ICA on the data matrix with a basis number of 5 was used as the extraction result.

Figure 5 shows the average correlation coefficients between the respiration curves and RIPsum obtained by the Nishida et al. method, ICA, and GMCA methods. This graph shows that the respiration curve extraction method using GMCA can extract respiration curves with higher accuracy than other methods.
Figure 6 shows the distribution of correlation coefficients for each method. Compared to ICA, GMCA has a higher correlation coefficient at the peak position, and it can be seen that it is possible to extract a respiration curve closer to the RIPsum.

Figure 7 shows the average correlation coefficient for each subject in the results of extracting the respiratory curve using GMCA with the initial condition Initial 10. The extraction accuracy of the respiratory curve for subject No. 12 was low. Of the 23 subjects, this subject had the shortest sleep time during the test and repeatedly got up from bed. It is possible that the extraction accuracy is low because it includes such unmeasurable time. However, even during times when the body pressure distribution time series can be stably measured, the extraction accuracy is not very high. Investigating the causes of poor extraction accuracy for some subjects and finding solutions is a future task.

(Embodiment 2) Detection of apnea events from a respiration curve In the above, a respiration curve is extracted. Here, an apnea event is detected from the respiration curve. Note that a respiration curve obtained by another method may be used.
<2-1: Division>
In order to detect apnea events, it is necessary to distinguish between normal breathing intervals and apnea intervals from respiration curves measured by various measurement methods. There are differences in respiration curves between subjects due to differences in the subject's physique, breathing method, and breathing rate. In tests on the same subject, the respiration curve changes due to changes in the subject's condition. The relationship between the testing equipment and the subject easily changes due to the subject's movements during sleep, and this change causes the respiration curve to change. The respiration curve includes not only changes between normal breathing and apnea, but also many changes due to such other factors. For this reason, it is difficult to appropriately divide the respiration curve by modeling all changes in the respiration curve. Therefore, in the embodiment, the respiration curve is divided using SST, which is a non-parametric time series analysis method.

(SST)
Singular Spectrum Transformation (SST) (Tsuyoshi Ide, Masaru Sugiyama. Anomaly Detection and Change Detection. Kodansha, 2015.) is a method for extracting feature patterns of time series data by singular value decomposition and determining the degree of change based on the extracted feature patterns.
An outline is shown in Figure 8.
Consider time series data y = {y1, y2,..., yt,...}. Let the partial time series x(t), which is a compilation of w pieces of data just before time t, be x(t) = [y(t-w), y(t-w+1),..., y(t-1)]T. On the past and future sides of time t, n partial time series are arranged, shifted by one time, and two data matrices X(t) and Z(t) with w rows and n columns are defined as X(t) = [x(t-n+1),..., x(t-1), x(t)].
Construct Z(t) = [x(t-n+1+L),...,x(t-1+L),x(t+L)].

X(t) and Z(t) are called the history matrix and test matrix, respectively. L is a non-negative integer that determines the relative positions of the history matrix and test matrix. The history matrix is composed of data from time t-w-n+1 to t, and contains past information with respect to time t. The test matrix is composed of data from time t-w-n+1+L to t+L, and contains future information with respect to time t.

Next, singular value decomposition is performed on the history matrix X(t) and the test matrix Z(t), and the left singular vectors with the largest singular values are α(t) and β(t), respectively. α(t) and β(t) represent characteristic patterns of the past and future, respectively. If the behavior of time series data y does not change before and after time t, α(t) and β(t) match. Conversely, if the behavior of y changes before and after time t, α(t) and β(t) do not match. The change in the behavior of y can be quantified by evaluating the similarity between α(t) and β(t). The degree of change a(t) is defined as a(t) = 1-α(t)Tβ(t). The degree of change a(t) is minimum when α(t) and β(t) match, and maximum when α(t) and β(t) are orthogonal.

The SST is applied to the respiration curve to calculate the degree of change a(t), and the point where the degree of change a(t) is at its maximum value is defined as the change point. The respiration curve between the change points is defined as a partial respiration curve. It is desirable for the partial respiration curve to include at least one respiratory waveform. Since the time it takes to take one breath during sleep is 3 to 6 seconds, in this embodiment, time t is defined as the change point when the degree of change a(t) is maximum in the 6 seconds before and after time t.

Figure 9 shows an example of the results of calculating RIPsum from the RIP signal measured in a PSG test and dividing the respiration curve using SST. The sampling rate of the RIP signal is 16 [Hz], and the SST parameters are w = n = 20, L = 10. The top row is the RIPsum, and the bottom row is the degree of change a(t) calculated using SST. The red line indicates the point of change. As a result of dividing the respiration curve, the section from 86 [s] to around 105 [s], where the amplitude of RIPsum is decreasing, can be extracted as a partial respiration curve. In the embodiment, the same settings of w = n = 20, L = 10 are used for the SST parameters.

<2-2: Clustering> Classification of partial respiratory curves In this section, we define the features of partial respiratory curves and explain a method for automatically classifying partial respiratory curves by clustering.
(Hierarchical Clustering)
Clustering is a method of automatically classifying data by dividing similar data into clusters without using external criteria. Hierarchical clustering is one of the clustering methods, in which similar clusters are joined in order from the state of each individual cluster, generating a hierarchical structure of clusters. A dendrogram is a representation of the hierarchical structure of clusters.
The process of hierarchical clustering is as follows.

1. Each piece of data is treated as a cluster, and the distance between clusters for all combinations is calculated.
2. Merge the closest cluster pairs.
3. Recalculate the distances between the clusters.
4. If the number of clusters is 2 or more, return to 2.

The Ward method is a method for recalculating the distance between clusters, and successively merges clusters so that the variance within the cluster is minimized. The Euclidean distance is used for the distance between individual data points. Suppose there is a set of N pieces of data X = {x1,...,xN}, and X is divided into C clusters G1,...,GC. The center of gravity M(Gi) of cluster Gi is calculated using the number of elements in cluster Gi, |Gi|, as M(Gi) = 1|Gi|Σxn∈Gixn. For cluster Gi, the sum of squares of the difference in distance between the center of gravity and each individual, E(Gi), is E(Gi) = Σxn∈Gi∥xn-M(Gi)∥2. Here, ΔE is defined as the distance between clusters, ΔE(Gi,Gj) = E(Gi∪Gj)-E(Gi)-E(Gj), and the cluster pair for which ΔE is the smallest is combined. The total sum E of intra-cluster variance is E = ΣCi = 1E(Gi), and the Ward method generates a cluster structure that minimizes E.

The respiratory curve is a measurement of bodily movement due to breathing. The relationship between the measuring equipment and the subject during sleep can easily change if the subject turns over in their sleep. In addition, bodily movement due to breathing varies depending on the subject, breathing state, and sleep state. As such, the respiratory curve changes due to various factors, making it difficult to model accurately. The criteria for determining apnea/hypopnea in PSG tests do not cover all respiratory curve models. Therefore, partial respiratory curves are classified using clustering, a method for grouping data into clusters without a teacher.

In the embodiment, the features of the partial respiratory curves are extracted from the partial respiratory curves using the autoregressive coefficients and the partial respiratory curve amplitude ratios, and the partial respiratory curves are classified by performing hierarchical clustering using the features of the partial respiratory curves as data. The aim is to automatically classify the partial respiratory curves into normal breathing intervals and apnea intervals. It is also preferable to classify them according to subjects, breathing method, and sensing state.

<Example>
Extraction of apnea intervals 1: Data set The sum of the RIPsum output values of the chest and abdominal RIP belts measured in the PSG test conducted at the Ohta Sleep Science Center was calculated, and the RIPsum signal was analyzed as a respiration curve. 50 intervals of 6 epochs were randomly extracted from the data of 23 subjects, and a total of 300 epochs were analyzed.
The PSG test data is divided into sleep epochs of 30 seconds each. In the PSG test data to be analyzed, epochs in which an apnea interval of 10 seconds or more occurs are labeled by a specialized test technician as abnormal breathing, such as obstructive apnea. Epochs that are not labeled as abnormal breathing are called normal breathing epochs, and epochs that are labeled as abnormal breathing are called apnea/hypopnea epochs. In addition, the occurrence of an apnea interval of 10 seconds or more is called an apnea event.

2: Cluster structure of partial respiration curves First, 50 RIPsum signals with a length of 6 epochs were divided into partial respiration curves using the SST method in 2-1. The number of divided partial respiration curves was 625. Next, the partial respiration curves were clustered using the hierarchical clustering method in 2-2. Figure 10 is a dendrogram showing the hierarchical structure of the combination of partial respiration curves. The vertical axis shows the sum of the variance within the clusters when the clusters are combined. In this analysis, the combination of clusters is stopped when the number of clusters reaches 13, and the partial respiration curves are labeled. In other words, the partial respiration curves are classified into 13 clusters. The red horizontal line indicates the stage at which the combination of clusters is stopped. The partial respiration curves combined below the red line belong to the same cluster, and the partial respiration curves combined above the red line belong to different clusters. The clusters located on the left in Figure 10 are cluster 1, cluster 2, ... , cluster 13, in that order.

3: Identifying a cluster corresponding to an apnea interval Identifying a cluster corresponding to an apnea interval using the label attached in the PSG test. A partial respiratory curve belonging to a cluster corresponding to a normal respiratory interval is called a normal partial respiratory curve, and a partial respiratory curve belonging to a cluster corresponding to an apnea interval is called an apnea partial respiratory curve. When a normal partial respiratory curve is followed by an apnea partial respiratory curve and the length of the apnea interval is 10 seconds or longer, an epoch including a change point between the normal partial respiratory curve and the apnea partial respiratory curve is an apnea/hypopnea epoch. The apnea/hypopnea epoch includes the start point of the apnea partial respiratory curve.

The proportion of normal breathing intervals on the respiratory curve is greater than the proportion of apnea intervals on the respiratory curve. For this reason, the apnea partial breathing curve follows the normal breathing interval, and its starting point often becomes an apnea event. In other words, if a certain cluster corresponds to an apnea interval, the starting points of the partial breathing curves belonging to that cluster will be included in normal breathing epochs at a low rate, and in apnea/hypopnea epochs at a high rate. For each cluster, the number of starting points of the partial breathing curves belonging to that cluster that are included in normal breathing epochs and in apnea/hypopnea epochs are tallied to determine which cluster corresponds to the apnea interval.

Table 1 shows the number of starting points of partial respiratory curves in normal breathing and apnea/hypopnea epochs for each cluster. The clusters with the largest ratios of the number of starting points of partial respiratory curves included in apnea/hypopnea epochs to the number of partial respiratory curves in the cluster are clusters 7, 3, and 2, in descending order. These clusters are highly likely to be clusters corresponding to apnea intervals. However, for all clusters, the number of starting points of partial respiratory curves included in normal breathing epochs is greater than the number of starting points of partial respiratory curves included in apnea/hypopnea epochs.

Figure 11 shows the results of dividing the respiratory curves and clustering them. The vertical lines indicate change points, and the respiratory curves separated by the vertical lines are partial respiratory curves. The numbers shown in the upper left of the partial respiratory curves indicate the cluster number to which the partial respiratory curve belongs. The vertical lines at t = 30, 60, ..., 150 [s] indicate the boundaries of the epochs. The labels given by the technician during the PSG test are shown below the epoch intervals. OSA and OSH indicate that the epochs were labeled as obstructive apnea and obstructive hypopnea, respectively.

In Fig. 11(a), the epoch from t = 60 to 90 [s] is labeled as OSH. The partial respiratory curve from t = 85 to 100 [s] belongs to cluster 3, has a small amplitude, and corresponds to the hypopnea section.
In Fig. 11(b), the epochs t = 30 to 60 [s], t = 60 to 90 [s], and t = 150 to 180 [s] are labeled as obstructive apnea. Considering the amplitude of the respiration curve, the section t = 35 to 100 [s] can be considered an apnea section except for the increase in amplitude around 60 seconds. This section corresponds to the partial respiration curves belonging to clusters 3 and 7. The apnea section of the last epoch begins 5 seconds before the end of the final epoch.

In this embodiment, the respiratory curve was separated into six epochs for analysis, so the data for the apnea period used to determine the last epoch was not included in the analysis. By performing the analysis including the next epoch, it should be possible to detect the presence of partial respiratory curves belonging to cluster 3, as with the apnea event from t = 30 to 60 [s].

In Figure 11(c), the epochs of t = 60 to 90 [s] and t = 120 to 150 [s] are labeled as obstructive hypopnea. The change in the amplitude of the respiratory curve is more gradual than in Figures 11(a) and 11(b), but the two partial respiratory curves of t = 70 to 100 [s] belonging to cluster 7 correspond to the apnea interval. As with the difference in the amplitude change, the cluster corresponding to the apnea interval has changed.

In Figure 11(d), the two partial respiratory curves from t = 135 to 155 [s] belong to cluster 2. There is a correspondence between the apnea period labeled as obstructive hypopnea in the epoch from t = 120 to 150 [s] and the two partial respiratory curves belonging to cluster 2. We were also able to confirm that the partial respiratory curves in cluster 2 matched the apnea period.

From Figure 11 and the above discussion, we can see that the partial respiratory curves belonging to clusters 2, 3, and 7 correspond to apnea intervals. However, according to Table 1, more than half of the partial respiratory curves belonging to clusters 2, 3, and 7 have a normal breathing epoch as their starting point. From Figure 11, we can see that even though the partial respiratory curves belonging to clusters 2, 3, and 7 actually correspond to apnea intervals, there are two cases in which their starting point is a normal breathing epoch.

The first is when a partial respiratory curve belonging to clusters 2, 3, or 7 appears following a normal breathing section, and normal breathing returns within a short period of time. The partial respiratory curve from t = 140 to 150 [s] in Figure 11b is one example. This partial respiratory curve has a smaller amplitude than the partial respiratory curves before and after it, and is an apnea section. However, the length of the section is short at 10 seconds, so it is not determined to be the start of an apnea event. As a result, in Table 1 it is tallied as a partial respiratory curve whose starting point is the normal breathing epoch.

The second case is when an apnea interval consisting of multiple partial respiratory curves exists across epochs, and the starting point of the partial respiratory curve exists after the epoch boundary. This can be seen in the interval after t = 128 [s] in Figure 11c. This interval is made up of three partial respiratory curves belonging to cluster 7, and corresponds to an apnea interval when viewed from the amplitude of the respiratory curves. When an apnea interval exists beyond the epoch boundary, only the previous epoch becomes an apnea/hypopnea epoch. The starting point of the partial respiratory curve from t = 128 to 158 [s] is included in the apnea/hypopnea epoch. However, the starting points of the partial respiratory curves from t = 158 to 170 [s] and t = 170 to 180 [s] are counted as normal breathing epochs.

Based on Table 1, the partial respiratory curve sections of clusters that are likely to correspond to apnea sections are counted as either normal breathing sections or apnea sections. Among the partial respiratory curves starting in normal breathing epochs in Table 1, the two patterns in which the partial respiratory curves mentioned above correspond to apnea sections are counted as apnea sections. In other words, for partial respiratory curves that belong to clusters 2, 3, and 7 starting from a normal breathing epoch, the analysis results are visually confirmed and counted for cases where the following partial respiratory curve belongs to a cluster other than cluster 2, 3, and 7 and is short in length, and cases where the following partial respiratory curve is included in clusters 2, 3, and 7 starting from an apnea/hypopnea epoch. Partial respiratory curves that start in normal breathing epochs that do not meet this condition are counted as normal breathing sections. In addition, partial respiratory curves that start in apnea/hypopnea epochs are counted as apnea sections.

Table 2 shows the results of the tabulation, showing the correspondence between the partial respiratory curves and the normal and apnea intervals for clusters 2, 3, and 7. For all clusters 2, 3, and 7, 66% or more of the partial respiratory curves correspond to apnea intervals. Among the partial respiratory curves considered to be normal respiratory intervals for cluster 2 in Table 2, there were six that began slightly before the start of the apnea/hypopnea epoch and most of the time interval was included in the apnea/hypopnea epoch. Taking this into consideration, the proportion of the partial respiratory curves for cluster 2 that correspond to apnea intervals reaches 80%. From the above discussion, we conclude that clusters 2, 3, and 7 are sets of partial respiratory curves that correspond to apnea intervals.

(Overall)
The embodiments can be combined with each other.
The embodiments can be combined with each other.
The above extracts apnea phenomena, but it can also extract respiratory system abnormalities such as COVID-19 and SARS (Severe Acute Respiratory Syndrome).

The alarm method of the present invention and devices using it can be used in watches, alarm clocks, portable devices, and mobile phones.

10 Person 11 Sensor 12 Control unit

Claims (3)

measuring a two-dimensional human body vibration waveform and converting the measured first data into frequency domain data;
a GMCA step of performing generalized morphological component analysis based on sparsity of signal representation on the second data converted in the conversion step to obtain second data;
and an inverse transformation step of transforming the third data obtained in the GMCA step into discrete signal data;
The generalized morphological component analysis is
determining the source signals and the weights by iterative calculation so that the sum of the weighted source signals of the vibration waveform coincides with the vibration waveform and each of the source signals can be expressed by a basis that is the smallest possible signal code;
In the iterative calculation of the generalized morphological component analysis, for a coefficient vector of a first source signal, a coefficient corresponding to a respiratory frequency band is set to 1 and other coefficients are set to 0 as initial values;
Measuring the two-dimensional vibration waveform of the human body means measuring two-dimensional pressure distribution at multiple measurement points using a surface pressure sensor, a respiratory curve extraction method. For the coefficient vector of the second source signal, the coefficient corresponding to the heart rate frequency band is set to 1 and the other coefficients are set to 0 as initial values.
For the coefficient vector of the third source signal, as an initial value, the coefficient corresponding to a frequency band lower than the respiratory frequency band is set to 1, and the other coefficients are set to 0;
For the fourth coefficient vector of the source signal, as an initial value, a coefficient corresponding to a frequency band between the respiratory frequency band and the heart rate frequency band is set to 1, and other coefficients are set to 0;
2. The method for extracting a respiratory curve according to claim 1, wherein, for the remaining coefficient vectors of the source signal, a frequency band higher than the heart rate frequency band is divided as an initial value, and the coefficients corresponding to the divided frequency bands are set to 1 and the other coefficients are set to 0. The respiratory curve extraction method according to claim 1 or 2, wherein the conversion step uses a discrete cosine transform.