Summary of the invention:
Purpose of the present invention, be to have in order to overcome existing polygraph that price is high, the shortcoming of complicated operation, solve the problem of traditional analysis technology such as Fourier analysis, auto-regressive analysis method complete failure, a kind of measuring method of the ECG of utilization signal detection sleep apnea is provided.
Purpose of the present invention can reach by taking following measure:
Utilize the measuring method of ECG signal detection sleep apnea, it is characterized in that:
Utilize portable sleep apnea detection record analyser that the patient is carried out whole night monitoring, record electrocardio, mouth and nose air-flow, oxygen saturation signal utilize the one section electrocardiosignal that contains sleep apnea to carry out following steps then and handle,
1) utilizes the slope threshold method to carry out the QRS ripple and detect, calculate R-R interval time series, and utilize the least square line match that R-R interval time series is carried out local detection to obtain heart rate variability signals HRV;
2) HRV with the variation of heart beating number of times as a kind of stochastic signal, select successive 500~800 (the best is 600) data points as sliding window, HRV signal in the window being carried out empirical modal decompose EMD, is one group of inherent mode function IMF with the HRV signal decomposition;
3) each IMF is carried out the Hilbert conversion and obtain the HH spectrum, both the time of the amplitude of HRV signal and frequency distributed;
4) ratio that calculates the energy of standard deviation, each IMF component of meansigma methods, the instantaneous amplitude of the instantaneous frequency of each IMF component and gross energy according to the HH of IMF spectrum goes out the time and the number of times of sleep apnea as eigenvalue according to the change-detection of eigenvalue.
Purpose of the present invention can also reach by taking following measure:
One embodiment of the present invention are: carry out the aforementioned the 1st) go on foot when gathering, in order to get rid of the abnormal R-R interval that causes because of the omission of R ripple in the QRS ripple testing process, adopt these impulsive noises of rolling average window filter filtering, its method is, get successive 41 R-R interval data, calculate the meansigma methods of other 40 data except that the 21st data point, compare with the 21st data point and this local mean values then, if this data point greater than local mean values 20% or be lower than 20% of local mean values, then it is removed; Change the window's position then, repeat aforementioned calculation all data points are carried out Filtering Processing.
One embodiment of the present invention are: the aforementioned the 1st) utilize the least square line match to carry out the method that local detection obtains heart rate variability signals HRV in the step to be, get successive 80 data points and form sliding window, carry out the least square line match, computing formula as shown in the formula:
Si=yi-(axi+b)
Wherein Si is the HRV signal after the detection, and yi is a heart rate signal, and xi is a time data, and the computing formula of a, b is as follows:
One embodiment of the present invention are: carry out the aforementioned the 2nd) when going on foot operation, select successive 600 data points as sliding window, the data in the window are carried out empirical modal decompose EMD, be one group of inherent mode function IMF with the HRV signal decomposition; Then each IMF is carried out the Hilbert conversion and obtain the HH spectrum, both the time of the amplitude of HRV signal and frequency distributed.
One embodiment of the present invention are:
1) the IMF component of the meansigma methods of instantaneous frequency between 0.02Hz-0.055Hz as responsive frequency range, corresponding instantaneous amplitude standard deviation and energy likened to be eigenvalue, if the instantaneous amplitude standard deviation is between 0.02 to 0.6, energy can be judged to sleep apnea than greater than 22.5%;
2) step-length is 100 data points, progressively changes position of window, and the variation by the eigenvalue in the responsive frequency range realizes apneic detection.
The present invention has following beneficial effect:
The present invention presents the obvious periodic fluctuation according to heart rate signal when asphyxia occurs, the spatial distribution of some inband energy of signal is compared during with eupnea respective change can be taken place, calculate the performance number of each frequency range respectively, and each frequency band power accounts for the percentage ratio of general power as eigenvalue, this unconspicuous signal characteristic form with significant energy variation in some frequency band is showed, adopt common cardiac monitoring technology to realize the detection of sleep apnea, compare with the sleep monitor of existing costliness, has simplicity of design, make detection reach minimum, will be extremely important the family oriented of sleep monitoring and the early diagnosis of sleep disordered breathing to the influence of patient's sleep.
The specific embodiment
Present embodiment utilizes portable sleep apnea detection record analyser that the patient is carried out whole night monitoring, record electrocardio, mouth and nose air-flow, oxygen saturation signal, utilize the one section electrocardiosignal that contains sleep apnea to carry out the processing of above-mentioned steps then, Fig. 1 is the R-R interval time series that obtains through the identification of QRS ripple, and Fig. 2 is that the selection window width is the HRV data that 80 detections obtain.
Fig. 3 is the HRV signal of each 600 data point of same different periods of individuality with Fig. 4, obtain six rank IMF by empirical mode decomposition, x among the figure (t) expression be primary HRV signal, imf1-imf6 represents six rank IMF, ref represents is residual error after primary signal deducts each rank IMF.From Fig. 3,4 as can be seen, the first rank IMF comprises the local high-frequency composition of signal, and along with the IMF exponent number increases, the IMF radio-frequency component progressively reduces.Compare with the global frequencies composition decomposition of fourier method, this is olation can disclose the time-varying characteristics of HRV signal intermediate frequency rate well by the decision of local frequencies composition.By the IMF of Fig. 3 and Fig. 4, can calculate the HH marginal spectrum of two sections HRV respectively shown in Fig. 5,6 by (14) formula.From HH limit spectrogram as can be seen, in low-frequency range, the total range value when pairing total range value was apparently higher than eupnea when sleep apnea took place that is to say the fluctuation of the heart rate of following asphyxia and producing, and has changed the HRV energy distributions.But because Fig. 5,6 characterizes is the distribution situation of signal energy accumulation on each Frequency point in the whole time span, changes frequency range greatly in order further to determine energy value, need carry out power spectrumanalysis such as Fig. 7, shown in Figure 8 to each IMF component.
Utilize Hilbert-Huang transfer pair heart rate variability signals to analyze and extract, utilize the variation of eigenvalue to realize that sleep apnea detects with eigenvalue.
Theoretical basis of the present invention is as follows:
1, empirical modal decomposes
Empirical modal decomposes supposes that at first the signal data of being gathered is to be formed by stacking by many basic inherent mode functions (IMF), becomes signal decomposition several eigen mode state functions then.And each inherent mode function has been represented a simple mode of oscillation, and they are linear or non-linear, but must satisfy following two conditions: 1. on whole signal length, the number of extreme point and zero crossing must equate or differ one at the most; 2. at any time, be zero by the coenvelope line of maximum point definition with by the meansigma methods of the lower envelope line of minimum point definition, that is to say that the envelope up and down of signal is symmetrical in time shaft.
In the processing procedure of actual signal, it is unpractical satisfying second condition fully, so as long as the meansigma methods of the two less than a predetermined a small amount of.According to definition, can adopt following method analytic function:
If time series is X (t), then:
(1) finds out all maximum points and the minimum point of X (t), it is fitted to the upper and lower envelope of former data sequence respectively with cubic spline function.The average of upper and lower envelope is average packet winding thread m1, former sequence is deducted m1 just can obtain anew sequences h 1 of removing low frequency, promptly
X(t)-m1=h1 (1)
General h1 is a stationary sequence not necessarily, needs it is repeated said process for this reason.If the average packet winding thread of h1 is m11, the sequence of then removing behind the low-frequency component of this envelope representative is h11, promptly
h1-m11=h11 (2)。
Repeat above process,,, finish to decompose, make h1k become first IMF item, promptly up to satisfying above two qualificationss through k time
h1(k-1)-m1k=h1k (3)
Order
C1=h1k (4)
C1 first IMF component for going out from this data separating, the minimum time yardstick that it comprises signal i.e. the most short-period mode, and the difference that makes primary signal and C1 is residual signal r1
Be X (t)-C1=r1(5)
Because r1 has comprised the component of longer cycle, it can be considered as initial data and repeat above process obtaining C2, its new difference is:
r1-C2=r2 (6)
......
rn-1-Cn=rn (7)
When rn is a monotonic sequence or relatively original signal amplitude is minimum when ignoring, think the process of extracting the inherent mode of signal of finishing.
Comprehensive above equation (5) (6) and (7) obtain at last:
From equation (8) as can be known, primary signal can be expressed as a n inherent mode function Cj and a residual signal rn sum, rn or constant amount or one group of dull data.All IMF components can revert back initial data X (t) at last through reverse stack.
Hilbert (Hilbert) conversion
Hilbert transform is different with other conversion, and it still transforms to time domain to signal from time domain.If X (t) is a time sequence, Y (t) is its Hilbert transform, promptly
X (t) and Y (t) are formed a plural number, can obtain the corresponding analytic signal Z (t) of X (t)
Instantaneous amplitude a (t), the instantaneous phase θ (t) of signal are expressed as follows
Its instantaneous frequency definition is
For actual signal, utilize the HT conversion to obtain its conjugation quadrature component, then real signal is carried out analytic representation, can obtain three instantaneous characteristic parameters of this signal, be instantaneous amplitude, instantaneous phase and instantaneous frequency, thereby realize the truly extraction of instantaneous parameters.But, by the physical significance of instantaneous frequency as can be known, be not that signal can both be discussed with instantaneous frequency arbitrarily.Only contain a kind of mode of oscillation when signal satisfies, do not have the situation of complicated stack ripple just feasible.That is to say that for the signal sequence that contains complicated wave system directly the HT conversion has lost the effectiveness on the method to a certain extent.For this reason, must carry out empirical modal to the signal sequence that contains complicated wave system and decompose (EMD), obtain a series of frequency contents inherent mode function (IMFs) component from high to low, again each component be carried out the HT conversion and obtain the HH spectrum, and then obtain marginal spectrum.
2, HH spectrum and marginal spectrum
Each IMF of formula (12) done to add up after the HT conversion
Here omitted residual error function rn, Re represents treating excess syndrome portion.Express (13) and be called the HH spectrum.
By formula (11) (13), signal amplitude a (t) and instantaneous frequency ω (t) are the functions of time, therefore can be presented at amplitude on frequency-time plane, promptly constitute the HH amplitude spectrum, the HH spectrum has accurately been described the rule that the amplitude of signal changes with frequency and time on whole section.Because energy can be with square the describing of amplitude, thus H (ω t) has also reflected the regularity of distribution of signal energy on the various yardsticks of frequency (or time) to a certain extent.HH spectrum H (ω, t) determine after, just can utilize following formula that time integral is obtained HH marginal spectrum (marginal Hilbert spectrum):
The HH marginal spectrum provides each frequency values pairing total range value, has characterized the distribution situation of signal energy accumulation on each Frequency point in the whole time span on statistical significance.
The solution of HHT method boundary effect
EMD decomposes IMF one by one by calculating envelope repeatedly.In computational process each time, come the local mean values of signal calculated according to the upper and lower envelope of signal; Upper and lower envelope is to be provided by cubic spline interpolation (Cubic Spline) algorithm by the local maximum of signal and minimum, and its structure affects the overall process of EMD, is determining the result that EMD decomposes, and is the key issue of HHT conversion.Because the signal two ends can not be in maximum and minimum simultaneously, thereby the average of boundary needs to estimate, when envelope is carried out spline interpolation, to outwards carry out continuation to signal or its extreme value, guaranteeing that envelope arrives at end points, otherwise or can be because of producing the integrity that bigger swing has a strong impact on data at the end points place, or understand because of the diffusion at end points place spectrum energy, signal is polluted to the center of data by two ends and destroys whole segment data, finally makes the EMD decomposition failure.This paper selects mirror image border continuation method for use, be exactly according to end points with apart from the symmetrical centre of the magnitude relationship between two nearest extreme values of end points decision mirror image continuation, when the value of end points between two nearest extreme values of end points the time, being that point of symmetry carry out the mirror image continuation from the nearest extreme point of end points.Otherwise, be that point of symmetry carry out the mirror image continuation with the end points.
Performing step of the present invention is as follows:
(1) earlier electrocardiogram (ECG) signal that collects is carried out filtering,, utilize the slope threshold method to carry out the QRS ripple then and detect, calculate R-R interval time series with the filtering baseline drift.
(2) in order to get rid of the abnormal R-R interval that causes because of the omission of R ripple in the QRS ripple testing process, we adopt these impulsive noises of rolling average window filter filtering.Its method is to get successive 41 R-R interval data, calculate the meansigma methods of other 40 data except that the 21st data point, compare with the 21st data point and this local mean values then, if this data point greater than local mean values 20% or be lower than 20% of local mean values, then it is removed.Change the window's position then, repeat aforementioned calculation all data points are carried out Filtering Processing.
(3) utilize the least square line match to carry out local detection and obtain heart rate variability (HRV) signal.Its method is to get successive 80 data points to form sliding window, carries out the least square line match, and computing formula is as (16) formula:
Si=yi-(axi+b) (16)
Wherein Si is the HRV signal after the detection, and yi is a heart rate signal, and xi is a time data, and the computing formula of a, b is shown in (17) (18):
(4) HRV with the variation of heart beating number of times as a kind of stochastic signal, select successive 600 data points as sliding window, the data in the window are carried out empirical modal decompose (EMD), be one group of inherent mode function (IMF) with the HRV signal decomposition.Then each IMF is carried out the Hilbert conversion and obtain the HH spectrum, both the time of the amplitude of HRV signal and frequency distributed.
(5) calculate meansigma methods, the standard deviation of instantaneous amplitude, the energy of each IMF component and the ratio of gross energy of the instantaneous frequency of each IMF component according to the HH of IMF spectrum.
(6) the IMF component of the meansigma methods of instantaneous frequency between 0.02Hz-0.055Hz as responsive frequency range, its corresponding instantaneous amplitude standard deviation and energy liken to and are eigenvalue, step-length is 100 data points, progressively change position of window, if when having asphyxia to occur in the window, bigger variation can take place in the eigenvalue in the responsive frequency range, and the variation by eigenvalue realizes apneic detection.
Table 1: each IMF component instantaneous frequency meansigma methods
| IMF1 | IMF2 | IMF3 | IMF4 | IMF5 | IMF6 |
| Eupnea (Hz) | 0.219091 | 0.080622 | 0.039806 | 0.017478 | 0.0084146 | 0.0040864 |
| Asphyxia (Hz) | 0.2139 | 0.079383 | 0.041595 | 0.028977 | 0.062419 | 0.0050818 |
The power spectrum of each IMF component has reflected the regularity of distribution of signal energy on the various yardsticks of frequency, but will accurately analyze the difference between each different objects, also must measure from quantitative angle.Therefore, be necessary on the basis of trying to achieve the instantaneous amplitude regularity of distribution, to extract the index of reflection sleep apnea feature.Back each layer energy decomposed in definition now
Wherein N is the number of data points of each IMF.
Decomposing back each layer gross energy is
Wherein Ej is the energy of j layer, the number of plies of n for decomposing.The energy ratio that defines each layer is:
What standard deviation reflected is overall degree of divergence, has characterized instantaneous amplitude a (t) random fluctuation in time or the intensity of vibration.The energy of table 2 each IMF during for the sleep apnea that calculates and eupnea than and the standard deviation of instantaneous amplitude, data have reflected that substantially sleep apnea causes the fluctuation situation of heart rate in the table.
The standard deviation of the energy of each IMF ratio and instantaneous amplitude when table 2 sleep apnea and eupnea
| IMF1 | IMF2 | IMF3 | IMF4 | IMF5 | IMF6 |
| The asphyxia energy is than Mj (%) | 7.19 | 20.9 | 30.62 | 19.98 | 15.34 | 5.98 |
| Asphyxia instantaneous amplitude standard deviation | 0.01131 | 0.01985 | 0.02684 | 0.01063 | 0.00605 | 0.004619 |
| The eupnea energy is than Mj (%) | 5.17 | 20.97 | 13.01 | 20.98 | 27.74 | 12.14 |
| Eupnea instantaneous amplitude standard deviation | 0.00639 | 0.01714 | 0.01158 | 0.01498 | 0.0118 | 0.002012 |
By repeatedly experiment contrast discovery, the average instantaneous frequency of IMF3 is generally between 0.02Hz to 0.05Hz, the standard deviation of energy ratio and instantaneous amplitude produces bigger variation when sleep apnea occurs, illustrate that these eigenvalues can reflect the change procedure of sleep-respiratory delicately.By changing the position of sliding window, can carry out the detection and the location of sleep apnea quickly and accurately.
Need to prove, Hilbert-Huang transform (Hilbert-Huang Transform, be called for short HHT) be after wavelet transformation, another that was proposed in 1998 by people such as U.S. scientist N.E.Huang is mainly used in the new method that non-stationary signal is analyzed, this method not only goes for the analysis of linear process, and is applicable to non-linear and analysis nonstationary time series.The core concept of HHT is that signal is decomposed (empirical mode decomposition by empirical modal, EMD), resolve into several intrinsic mode functions (Intrinsic Mode Function, IMF) utilize Hilbert transition structure analytic signal then, draw the instantaneous frequency and the amplitude of signal, and then obtain the Hilbert spectrum[4]Because the EMD method is the Time Domain Decomposition that the time-domain information of basis signal itself carries out, the common number of the IMF that obtains be finite sum stably, and be narrow band signal with practical significance, its result of Hilbert conversion who carries out based on these IFM components has reflected real physical message, therefore, its Hilbert spectrum also can accurately reflect signal energy, the distribution of frequency on space or time scale.It is not subjected to Fourier analysis requirement signal must be linear, steady and periodic limitation, has whole advantages of wavelet analysis again, and eliminated the fuzzy and unintelligible of wavelet analysis on resolution, has and composes structure more accurately.Therefore this Hilbert frequency spectrum analysis method based on EMD has very high using value in the analysis of non-linear and non-stationary process.