Respiratory curve extraction method

Abstract

To provide a pneumogram extraction method and an apneic event detection method that can accurately extract a pneumogram.SOLUTION: A pneumogram extraction method includes a conversion step for measuring two-dimensional human body oscillatory waveforms and transforming the measured first data to frequency region data, a GMCA step for subjecting the second data, transformed in the conversion step, to morphological component analysis, to obtain second data, and an inverse transformation step for inversely transforming the third data, obtained in the GMCA step, to discrete signal data.SELECTED DRAWING: Figure 1

Description

[Technical field]

[0001] The present invention relates to a method for extracting a respiratory curve and a method for extracting an abnormal respiratory event. The method extracts a respiratory curve, and in particular extracts apnea events and abnormal events in the respiratory system. [Background technology]

[0002] Conventionally, there have been methods for extracting a respiratory curve. In Patent Document 1, the respiratory curve is calculated by the reciprocal of the standard deviation of the respiratory frequency. [Prior art documents] [Patent documents]

[0003] [Patent Document 1] Patent No. 5670071 Summary of the Invention [Problem to be solved by the invention]

[0004] However, the method of Patent Document 1 was unable to extract a respiration curve accurately.

[0005] Therefore, an object of the present application is to provide a respiration curve extraction method and an abnormal respiratory event extraction method that can accurately extract a respiration curve. [Means for solving the problem]

[0006] In order to solve the above problem, a respiration curve extraction method is used, which includes a conversion step of 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 (GMCA) on the second data converted in the conversion step based on the sparsity of the signal representation by morphological component analysis to obtain the second data, and an inverse conversion step of converting the third data obtained in the GMCA step into discrete signal data.

[0007] Also, an abnormal respiratory event detection method is used, which includes a division step of dividing a respiratory curve into a plurality of epochs, a classification step of classifying the plurality of epochs into a plurality of clusters, and an abnormal respiratory interval extraction step of extracting clusters of abnormal respiratory intervals. Effect of the Invention

[0008] According to the respiration curve extraction method and abnormal respiratory event detection method of the present invention, a respiration curve can be accurately extracted, and abnormal respiratory events can also be extracted. [Brief description of the drawings]

[0009] [Figure 1] FIG. 1 is a diagram showing an outline of a respiratory motion extraction method according to an embodiment. [Diagram 2] FIG. 2 is a plan view illustrating pressure measurement according to an embodiment; [Diagram 3] FIG. 13 is a diagram showing average values ​​of correlation coefficients when initial values ​​are set, which are considered in the BSS process according to the embodiment. [Figure 4] FIG. 13 is a diagram showing the distribution of correlation coefficients under each initial value condition of the BSS process according to the embodiment; [Diagram 5] The figure shows the average correlation coefficient between the respiratory curves obtained by the Nishida, ICA, and GMCA methods and RIPsum. [Figure 6] Distribution of correlation coefficients for each method [Figure 7] The figure shows the average correlation coefficient for each subject in the results of respiratory curve extraction using GMCA with the initial condition Initial10. [Figure 8] Diagram showing the outline of the singular spectrum transformation [Figure 9] This 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. [Figure 10] A dendrogram showing the hierarchical structure of the connections of partial respiration curves. [Figure 11] (a) to (d) Graphs showing the results of dividing and clustering the respiratory curves. DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENTS

[0010] (Embodiment 1) Respiration curve extraction using generalized morphological component analysis <Outline of the respiratory curve extraction method> First, we will explain how to analyze a person's pressure signal (observed signal) and obtain a respiratory curve. <Respiration motion extraction method> FIG. 1 shows an outline of a respiratory motion extraction method according to an embodiment.

[0011] 1-1: Measure the pressure signal Bias is removed from the pressure time series at each measurement point by subtracting the epoch-averaged pressure value at the measurement point. The bias-removed pressure time series is used to create an observed signal matrix X. The measured body pressure distribution time series is processed by the following GMCA (described below) for each sleep epoch (30 seconds).

[0012] 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 ​​of the mixing matrix and the original signal matrix, ~A and ~ΘS, are obtained.

[0013] 1-4: Reverse data conversion Then, an inverse DCT is applied to ~ΘS. An estimate ˜S of the source signal matrix 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 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.

[0014] <Details of the method for extracting respiratory 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.

[0015] <1: Measure the pressure signal 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 caused by body weight. The time series of body pressure distribution measured by this device includes periodic pressure changes caused by respiratory effort and heart rate.

[0016] The signal σ measured by a pressure sensor at a point can be modeled as follows: σ(t)=SH(t)+SR(t)+M(t)+B(t)+N(t) Here, SH and SR are the signals due to breathing and heart rate, respectively, M is the signal due to body movements such as turning over in bed, B is the bias due to weight, and N is 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.

[0017] <1-2 Data conversion> The 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, which is transformed into N coefficients C0,...,CN-1, the DCT is defined as follows:

number

number

[0018] 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. An embodiment will be described below, which 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 using 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).

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

[0020] Morphological Component Analysis (MCA) (J.-L Starck, YA, J. Bobin, Michael Elad, and D. Donoho. Morphological component analysis. Proceedings of SPIE - The International Society for Optical Engineering, Vol. 5914,, 08 2005.) was 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 bases are zero) using a specific base Φ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 the following equation 3.

number

[0021] where αk is the coefficient vector that associates y with Φk. Since the coefficient vector α is sparse, the minimization problem of the following equation 4 is obtained.

number

[0022] 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 to achieve 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

[0023] where αjk is the coefficient vector that associates sj with Φk. Since the source signals sj are sparsely represented using the dictionary D, the following optimization problem (6) can be defined for the decomposition of the observed signal matrix X.

number

[0024] By relaxing the l0 norm (∥·∥0) to the l1 norm and introducing the Lagrange multiplier λj, the problem in Equation 6 can be approximated by Equation 7 below.

number

[0025] However, 2 2 is the Frobenius norm, and ∥P∥ 2 2=Trace(P T P). In GMCA, the optimization problem of Equation 7 is solved using an iterative reduced threshold algorithm to obtain the estimated values ​​~S and ~A of the source signal matrix and the mixing matrix 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 sparsity constraint, and the solution of Equation 7 is obtained.

[0026] In the algorithm for solving Equation 7, when estimating the source signal matrix S by iterative calculation, a conversion by D must be performed at each calculation, which results in high computational costs. The fast algorithm for GMCA assumes that an inherent sparse decomposition exists for 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

[0027] Unless the dictionary D is an orthonormal basis, Equation 7 and Equation 8 are not equivalent. When the dictionary D is redundant, good results are obtained by experiment, although no rigorous mathematical proof is given. If A is fixed in equation 8, ΘS can be found by the order.

number

[0028] where A† is the pseudo-inverse matrix. Also, Sλ is a soft-thresholding operator that operates as follows when y[t] is the input:

number

[0029]

number

[0030] 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 represented by a smaller number of bases, and if the threshold is small, it is represented by a larger number of bases. When the threshold λ is large, the source signal is estimated from the most important coefficients, and as the threshold λ is 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.

[0031] 1-3-A Initial Value Settings for 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. The mixing matrix Aini contains information on the respiratory frequency band as prior knowledge. The following four conditions are considered as initial values ​​to be given to GMCA.

[0032] 1. Random: A ini Give a random value that follows the standard normal distribution. 2. PCA-based: Variance-covariance matrix of coefficient matrix of observed signal matrix Σ=Θ X Θ T X The eigenvectors vj (j=1,2,..,n) of A are calculated by principal component analysis (PCA). ini is given by n eigenvectors vj (j=1,2,..,m).

[0033] 3. Initial 1: 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 set as a1 = ΘXΘ T s1 The remaining weight coefficient vector aj(j=1,2,...,n) is given random values ​​following the standard normal distribution.

[0034] 4. Initial 10: As in Initial 1, 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.

[0035] 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.

[0036] 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.

[0037] <1-4: Reverse data conversion> The inverse transform of DCT is called the inverse discrete cosine transform (IDCT), and is defined as follows:

number

[0038] <Example> Extraction of respiration curves using generalized morphological component analysis 1: Dataset A SmartRubber (registered trademark) sensor 11 manufactured by Sumitomo Riko Co., Ltd. (hereinafter referred to as 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 for 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.

[0039] A PSG test is a test called polysomnography (PSG) 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 head, face, and necessary parts of the body, and brain waves, breathing, eyeball, and muscle movements are recorded while the patient sleeps overnight to examine the state of sleep. The results are analyzed and compared with the normal progression of sleep to diagnose sleep apnea syndrome (SAS).

[0040] 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.

[0041] 2: Accuracy of respiratory curve extraction A respiration curve extraction method using GMCA was applied to the time series of body pressure distribution during sleep for 23 subjects (10 subjects) throughout the night, and the accuracy of extracting the respiration curve was evaluated. The number of source signals was set to n=10. A total of 22,158 epochs were analyzed, and the respiration curve was extracted from the time series of body pressure distribution during sleep for each epoch, and the correlation coefficient was calculated.

[0042] Figure 3 shows the average correlation coefficients when setting the initial values ​​examined in 1-3-A. The average correlation coefficients were highest for Initial10, followed by Initial1, Random, and PCA-based. As 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.

[0043] 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 is smaller near the peak with a correlation coefficient of 0.1, and larger 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 extraction accuracy of respiratory movement with Initial10 compared to the other initial value conditions does not indicate that it extracts a respiratory curve closer to RIPsum, but rather 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.

[0044] 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 respiration and body position during sleep using pressure sensors. Journal of the Robotics Society of Japan, Vol. 16, No. 5, pp. 705-711, jul1998.: A method for calculating a respiratory curve from the phase difference of signals from multiple sensors), independent component analysis (ICA), and GMCA to the time series of body pressure distribution during sleep of the same 23 subjects. Initial value Initial10 was used 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.

[0045] 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 RIPsum.

[0046] 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 12 was low. Of the 23 subjects, this subject had a short sleep time during the test and repeatedly got up from bed. It is possible that the extraction accuracy was low because such unmeasurable time was included. However, even during times when the body pressure distribution time series could be stably measured, the extraction accuracy was not very high. An investigation into the causes of poor extraction accuracy for such subjects and countermeasures are a future task.

[0047] (Embodiment 2) Apnea event detection from respiration curve In the above, the respiration curve is extracted. In this step, apnea events are detected from the respiration curve. Note that the respiration curve may be obtained by another method. <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. When testing 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 while sleeping, and this change causes the respiration curve to change. The respiration curve contains 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.

[0048] (SST) Singular Spectrum Transformation (SST) (Ide Tsuyoshi, Sugiyama Masaru. Anomaly Detection and Change Detection. Kodansha, 2015.) is a method that extracts characteristic patterns of time series data using singular value decomposition and calculates the degree of change based on that. An outline is shown in Figure 8. Consider time series data y={y1,y2,...,yt,...}. The partial time series x(t) which is a compilation of w data just before time t is expressed as x(t)=[y(tw),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 expressed 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)].

[0049] 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 tw-n+1 to t, and contains past information with respect to time t. The test matrix is ​​composed of data from time tw-n+1+L to t+L, and contains future information with respect to time t.

[0050] 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 defined as α(t) and β(t), respectively. α(t) and β(t) represent characteristic patterns of the past and future, respectively. If the behavior of the time series data y does not change before and after time t, α(t) and β(t) will match. Conversely, if the behavior of y changes before and after time t, α(t) and β(t) will 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) are consistent, and maximum when α(t) and β(t) are orthogonal.

[0051] 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 that the partial respiration curve includes 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.

[0052] FIG. 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 upper part is RIPsum, and the lower part is the degree of change a(t) calculated by 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.

[0053] <2-2: Clustering> Classification of partial respiration curves In this section, we define the features of partial respiratory curves and explain a method to automatically classify partial respiratory curves using clustering. (Hierarchical Clustering) Clustering is a method of automatically classifying data by dividing it into clusters that group together similar data 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 in a tree diagram. The process of hierarchical clustering is as follows.

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

[0055] Ward's method is one of the techniques for recalculating the distance between clusters, and it is a method of successively merging clusters so that the variance within the cluster is minimized. The distance between individual data points is determined by Euclidean distance. 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 as M(Gi) = 1|Gi|Σxn∈Gixn, where |Gi| is the number of elements in cluster Gi. 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 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.

[0056] 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 due to the subject turning 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 and 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.

[0057] 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 feature quantities of the partial respiratory curves as data. The purpose is to automatically classify the partial respiratory curves into normal breathing sections and apnea sections. In addition, it is preferable to classify them according to subjects, breathing methods, and sensing conditions.

[0058] <Example> Extraction of apnea periods 1: Dataset The sum of the RIPsum output values ​​of the chest and abdomen RIP belts measured in PSG tests performed at the Ohta Sleep Science Center was calculated, and the RIPsum signal was analyzed as a respiratory curve. 50 sections of 6 epochs were randomly extracted from the data of 23 subjects, and a total of 300 epochs were analyzed. 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 occurred are labeled by a specialized technician as having respiratory abnormalities such as obstructive apnea. Epochs that are not labeled as having respiratory abnormalities are called normal breathing epochs, and epochs that are labeled as having respiratory abnormalities are called apnea / hypopnea epochs. In addition, the occurrence of an apnea interval of 10 seconds or more is called an apnea event.

[0059] 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 a cluster 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 cluster combination 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 order.

[0060] 3: Identifying clusters corresponding to apnea periods Identifying clusters corresponding to apnea intervals Clusters corresponding to apnea intervals are identified using the labels assigned in the PSG test. The partial respiratory curves belonging to the cluster corresponding to normal respiratory intervals are called normal partial respiratory curves, and the partial respiratory curves belonging to the cluster corresponding to apnea intervals are called apnea partial respiratory curves. 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, the epoch that includes the 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.

[0061] The proportion of normal breathing sections in the respiratory curve is greater than the proportion of apnea sections in the respiratory curve. For this reason, the apnea partial breathing curve follows the normal breathing section, and its starting point often becomes an apnea event. In other words, if a certain cluster corresponds to an apnea section, the starting points of the partial breathing curves belonging to that cluster will be included in normal breathing epochs less frequently, and in apnea / hypopnea epochs more frequently. 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 the number that are included in apnea / hypopnea epochs are counted to determine the cluster that corresponds to the apnea section.

[0062] 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 correspond 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. [Table 1]

[0063] 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 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.

[0064] In Fig. 11(a), the epoch from t = 60 to 90 [s] is labeled OSH. The partial respiration curve from t = 85 to 100 [s] belongs to cluster 3, and has a small amplitude corresponding to the hypopnea section. In Figure 11(b), the epochs t = 30-60 [s], t = 60-90 [s], and t = 150-180 [s] are labeled as obstructive apnea. Considering the amplitude of the respiration curve, the section t = 35-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 began 5 seconds before the end of the final epoch.

[0065] In the embodiment, the respiratory curve was separated into six epochs for analysis, so the data of the apnea interval 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 a partial respiratory curve belonging to cluster 3, as with the apnea event at t=30-60[s].

[0066] In Figure 11(c), the epochs 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 at 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.

[0067] 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 interval 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 also matched the apnea interval.

[0068] 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.

[0069] The first is when a partial respiratory curve belonging to clusters 2, 3, or 7 appears following a normal breathing period, and then returns to normal breathing 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 period. However, the length of the period is short at 10 seconds, so it is not determined to be the starting point of an apnea event. As a result, in Table 1 it is tabulated as a partial respiratory curve whose starting point is the normal breathing epoch.

[0070] 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 is seen in the section after t = 128 [s] in Figure 11c. This section is composed of three partial respiratory curves belonging to cluster 7, and corresponds to an apnea interval in terms of 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.

[0071] 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 follows a partial respiratory curve 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.

[0072] Table 2 shows the results of the aggregation, showing the correspondence between the partial respiratory curves of clusters 2, 3, and 7 and the normal and apnea intervals. For all of clusters 2, 3, and 7, more than 66% of the partial respiratory curves correspond to apnea intervals. Among the partial respiratory curves considered to be normal respiratory intervals in cluster 2 in Table 2, there were six that began shortly 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 of cluster 2 that correspond to apnea intervals reaches 80%. From the above discussion, we conclude that clusters 2, 3, and 7 are collections of partial respiratory curves that correspond to apnea intervals. [Table 2]

[0073] (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). [Industrial Applicability]

[0074] The alarm method of the present invention and devices using the same can be used in watches, alarm clocks, portable devices, and mobile phones. [Explanation of symbols]

[0075] 10 people 11 Sensors 12 Control section

Claims

1. 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.

2. 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 third coefficient vector of the 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.

3. 3. The method of claim 1, wherein the transform step uses a discrete cosine transform.

Images