CN110464337A - A kind of heart rate variability signals analysis method based on extreme value Energy Decomposition method - Google Patents
A kind of heart rate variability signals analysis method based on extreme value Energy Decomposition method Download PDFInfo
- Publication number
- CN110464337A CN110464337A CN201910841605.6A CN201910841605A CN110464337A CN 110464337 A CN110464337 A CN 110464337A CN 201910841605 A CN201910841605 A CN 201910841605A CN 110464337 A CN110464337 A CN 110464337A
- Authority
- CN
- China
- Prior art keywords
- signal
- extreme
- component
- energy
- rate variability
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Granted
Links
- 238000000034 method Methods 0.000 title claims abstract description 38
- 238000000354 decomposition reaction Methods 0.000 title claims abstract description 28
- 238000004458 analytical method Methods 0.000 title claims abstract description 18
- 230000002159 abnormal effect Effects 0.000 claims abstract description 7
- 239000013598 vector Substances 0.000 claims description 13
- 238000010183 spectrum analysis Methods 0.000 claims description 6
- 238000005070 sampling Methods 0.000 claims description 3
- 238000012216 screening Methods 0.000 claims description 3
- 230000000747 cardiac effect Effects 0.000 abstract description 4
- 206010019280 Heart failures Diseases 0.000 description 19
- 206010007559 Cardiac failure congestive Diseases 0.000 description 18
- 238000010586 diagram Methods 0.000 description 7
- 230000000694 effects Effects 0.000 description 3
- 230000007774 longterm Effects 0.000 description 3
- 230000033764 rhythmic process Effects 0.000 description 3
- 208000024172 Cardiovascular disease Diseases 0.000 description 2
- 230000007423 decrease Effects 0.000 description 2
- 230000003247 decreasing effect Effects 0.000 description 2
- 238000001514 detection method Methods 0.000 description 2
- 208000019622 heart disease Diseases 0.000 description 2
- 238000005259 measurement Methods 0.000 description 2
- 210000000653 nervous system Anatomy 0.000 description 2
- 238000007781 pre-processing Methods 0.000 description 2
- 238000004393 prognosis Methods 0.000 description 2
- 238000011160 research Methods 0.000 description 2
- 230000011218 segmentation Effects 0.000 description 2
- 238000001228 spectrum Methods 0.000 description 2
- 230000002123 temporal effect Effects 0.000 description 2
- 238000012360 testing method Methods 0.000 description 2
- 241000282412 Homo Species 0.000 description 1
- 238000004364 calculation method Methods 0.000 description 1
- 238000007405 data analysis Methods 0.000 description 1
- 201000010099 disease Diseases 0.000 description 1
- 208000037265 diseases, disorders, signs and symptoms Diseases 0.000 description 1
- 230000003993 interaction Effects 0.000 description 1
- 208000010125 myocardial infarction Diseases 0.000 description 1
- 238000011084 recovery Methods 0.000 description 1
- 238000000638 solvent extraction Methods 0.000 description 1
- 230000003595 spectral effect Effects 0.000 description 1
- 230000005654 stationary process Effects 0.000 description 1
Classifications
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
- A61B5/00—Measuring for diagnostic purposes; Identification of persons
- A61B5/02—Detecting, measuring or recording pulse, heart rate, blood pressure or blood flow; Combined pulse/heart-rate/blood pressure determination; Evaluating a cardiovascular condition not otherwise provided for, e.g. using combinations of techniques provided for in this group with electrocardiography or electroauscultation; Heart catheters for measuring blood pressure
- A61B5/0205—Simultaneously evaluating both cardiovascular conditions and different types of body conditions, e.g. heart and respiratory condition
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
- A61B5/00—Measuring for diagnostic purposes; Identification of persons
- A61B5/02—Detecting, measuring or recording pulse, heart rate, blood pressure or blood flow; Combined pulse/heart-rate/blood pressure determination; Evaluating a cardiovascular condition not otherwise provided for, e.g. using combinations of techniques provided for in this group with electrocardiography or electroauscultation; Heart catheters for measuring blood pressure
- A61B5/024—Detecting, measuring or recording pulse rate or heart rate
- A61B5/02405—Determining heart rate variability
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
- A61B5/00—Measuring for diagnostic purposes; Identification of persons
- A61B5/02—Detecting, measuring or recording pulse, heart rate, blood pressure or blood flow; Combined pulse/heart-rate/blood pressure determination; Evaluating a cardiovascular condition not otherwise provided for, e.g. using combinations of techniques provided for in this group with electrocardiography or electroauscultation; Heart catheters for measuring blood pressure
- A61B5/024—Detecting, measuring or recording pulse rate or heart rate
- A61B5/0245—Detecting, measuring or recording pulse rate or heart rate by using sensing means generating electric signals, i.e. ECG signals
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
- A61B5/00—Measuring for diagnostic purposes; Identification of persons
- A61B5/24—Detecting, measuring or recording bioelectric or biomagnetic signals of the body or parts thereof
- A61B5/316—Modalities, i.e. specific diagnostic methods
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
- A61B5/00—Measuring for diagnostic purposes; Identification of persons
- A61B5/24—Detecting, measuring or recording bioelectric or biomagnetic signals of the body or parts thereof
- A61B5/316—Modalities, i.e. specific diagnostic methods
- A61B5/318—Heart-related electrical modalities, e.g. electrocardiography [ECG]
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
- A61B5/00—Measuring for diagnostic purposes; Identification of persons
- A61B5/72—Signal processing specially adapted for physiological signals or for diagnostic purposes
- A61B5/7203—Signal processing specially adapted for physiological signals or for diagnostic purposes for noise prevention, reduction or removal
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
- A61B5/00—Measuring for diagnostic purposes; Identification of persons
- A61B5/72—Signal processing specially adapted for physiological signals or for diagnostic purposes
- A61B5/7235—Details of waveform analysis
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
- A61B5/00—Measuring for diagnostic purposes; Identification of persons
- A61B5/72—Signal processing specially adapted for physiological signals or for diagnostic purposes
- A61B5/7235—Details of waveform analysis
- A61B5/725—Details of waveform analysis using specific filters therefor, e.g. Kalman or adaptive filters
Landscapes
- Health & Medical Sciences (AREA)
- Life Sciences & Earth Sciences (AREA)
- Engineering & Computer Science (AREA)
- Molecular Biology (AREA)
- Animal Behavior & Ethology (AREA)
- Pathology (AREA)
- Physics & Mathematics (AREA)
- Biomedical Technology (AREA)
- Heart & Thoracic Surgery (AREA)
- Medical Informatics (AREA)
- Cardiology (AREA)
- Surgery (AREA)
- Biophysics (AREA)
- General Health & Medical Sciences (AREA)
- Public Health (AREA)
- Veterinary Medicine (AREA)
- Physiology (AREA)
- Signal Processing (AREA)
- Artificial Intelligence (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Psychiatry (AREA)
- Pulmonology (AREA)
- Measurement And Recording Of Electrical Phenomena And Electrical Characteristics Of The Living Body (AREA)
Abstract
The invention discloses a kind of heart rate variability signals analysis methods based on extreme value Energy Decomposition method, and the ECG signal including obtaining the unknown state under given time and given sample frequency obtains RRI signal x (t) after denoising；It regard RRI signal x (t) as original signal, original signal x (t) is decomposed into n extreme value mode function component and a surplus, the n extreme value mode function component that original signal x (t) is decomposed, the component of original signal different frequency range is represented, determines whether the RRI signal is abnormal cardiac rate Variability Signals according to n extreme value mode function component.The present invention analyzes RRI signal using extreme value Energy Decomposition method, and original signal is decomposed into multiple components, that is, extreme value component function, the energy of each component is calculated, obtains its Energy distribution.
Description
Technical Field
The invention relates to electrocardiogram signal analysis, in particular to a heart rate variability signal analysis method based on an extreme value energy decomposition method.
Background
Physiological signals are generated by the interaction of multiple systems of a living body, and the time and the intensity of the action of different systems are different, so that the physiological signals have complexity in time and space. Heart Rate Variability (HRV) refers to the measurement of the temporal variation between successive cardiac cycles, and specifically, the variation in the difference between the successively occurring normal P-P intervals. However, since P-waves are less pronounced than R-waves or the P-wave tips are sometimes blunt, studies on heart rate variability signals are typically replaced with R-R Interval Signals (RRIs) equal to P-P intervals. Research shows that HRV can be used as noninvasive detection index of activity of the vegetative nervous system, and has important significance in the aspect of judging the prognosis of certain cardiovascular diseases.
Long-term modulated rhythms (<1Hz)) of the heart are studied in the prior art using the heart rate variability signal (HRV signal) as the subject of analysis. Numerous studies have shown that human HRV signals have long-term correlations and nonlinear kinetic complexity, and that age and disease lead to reduced kinetic complexity. The Heart Rate Variability (HRV) signal is usually studied by RR interval (RRI) signal, i.e. the time interval between successive R waves of RRI signal.
The most common method of studying energy changes in HRV signals is power spectral analysis (PSD). The PSD converts the power of the HRV signal into a function of frequency by fourier transform, studies the power level of different frequency domain ranges, and generally the HRV spectrum is divided into several frequency bands, such as High Frequency (HF), Low Frequency (LF), and Very Low Frequency (VLF). The LF/HF ratio has important clinical value. Heart disease causes changes in the HRV power spectrum, such as heart failure and myocardial infarction resulting in increased normalized HF, decreased LF and VLF. However, PSD is not a data-driven based method and is coarser for segmentation of the frequency domain, resulting in loss of detail and not flexible enough for segmentation.
Therefore, it is desired to solve the above problems.
Disclosure of Invention
The purpose of the invention is as follows: the invention aims to provide a heart rate variability signal analysis method based on an extreme value energy decomposition method, which can intuitively reflect the real rule of electrocardiogram energy distribution and signal fluctuation by adopting less data.
The technical scheme is as follows: in order to achieve the above purpose, the invention discloses a heart rate variability signal analysis method based on an extremum energy decomposition method, which comprises the following steps:
(1) acquiring an ECG signal in an unknown state at a given time and a given sampling frequency, then carrying out denoising pretreatment on the ECG signal, and extracting an RRI signal from the ECG signal to obtain an RRI signal x (t) in the unknown state;
(2) taking the RRI signal x (t) as an original signal, calculating all local extreme points of the original signal, and then connecting all the extreme points and all the minimum points of the original signal by adopting a spline curve to respectively form an upper envelope line emaxAnd a lower envelope eminObtaining an envelope mean value signal m (t) ═ e of the upper and lower envelope linesmax+emin)/2；
(3) Subtracting the envelope mean signal m (t) from the original signal x (t) to obtain h (t) ═ x (t) -m (t); then judging whether h (t) meets the judgment condition of the extreme value modal function, if not, returning h (t) to the step (3) as the original signal until hk(t) if the condition for determining the extremum mode function is satisfied, c is recorded1(t)＝hk(t) as a first extreme modal function component;
(4) subtracting the first extreme mode function component c from the original signal x (t)1(t) obtaining the residue r1(t)＝x(t)-c1(t) then judging hk(t) whether a stopping criterion is met, and if not, r1(t) as a new original sequence x (t), returning to steps (2) and (3); if h isk(t) when the stopping criterion is met and n is less than 8, returning to the step (1) to obtain the original signal again; if h isk(t) when the stopping criterion is met and n is more than or equal to 8, obtaining 2 nd, 3 rd, … th extreme value modal function components and margin rn(t) the original signal x (t) is then decomposed into n extreme modal function components and a residual, i.e.
(5) Opposite polar mode function component ci(t), i is 1, 2, …, n, and the center frequency of each extreme modal function component is obtained by performing spectrum analysis;
(6) n extreme value modal function components obtained by decomposing the original signal x (t) represent components of different frequency bands of the original signal, and then the energy of each component is calculated
Ei＝∫|ci(t)|2dt，i＝1，2，…，n
Normalizing each energy value to obtain a normalized energy distribution vector
pi＝Ei/E，i＝1，2，...，n
Wherein,first component p1Representing the energy of the highest frequency band, representing the proportion of the energy distribution of the signal in the highest frequency band, the last component pnRepresenting the proportion of the energy distribution of the signal in the lowest frequency band range; drawing a normalized energy distribution graph according to the normalized energy distribution vector, wherein the abscissa represents component levels, the ordinate represents normalized energy distribution vector values, the curve represents an average value, and the error bar represents a standard deviation;
(7) selecting a second component p2To the seventh component p7Calculating the energy difference value EDV (p)2+p3+p4)-(p5+p6+p7) And when the EDV is less than or equal to a first standard value, the RRI signal is judged to be a normal heart rate variability signal, when the first standard value is less than or equal to the EDV and less than a second standard value, the RRI signal is judged to be a suspected abnormal heart rate variability signal, and when the EDV is more than or equal to the second standard value, the RRI signal is judged to be an abnormal heart rate variability signal.
Wherein, the minimum data amount N of the original signal x (t) is 2n+1Wherein n is inThe number of extremal modal function components resolved.
Preferably, the specific method of denoising pre-processing in step (1) is as follows: the ECG signal is filtered through a 40Hz zero phase FIR low pass filter to remove high frequency noise, and then through a median filter to remove baseline wander.
In addition, the determination conditions of the extreme mode function in the step (3) are as follows: (a) in the whole data sequence, the number of the extreme points is equal to or different from the number of the zero-crossing points by one; (b) and the upper and lower envelopes are symmetric with respect to the time axis at any time.
Further, h in the step (4)k(t) satisfying the stopping criteria the formula:
epsilon represents a screening threshold, and is taken to be 0.2-0.3.
Preferably, the first standard value in the step (7) is-0.15, and the second standard value is 0.08.
Has the advantages that: compared with the prior art, the invention has the following remarkable advantages:
the method adopts an Extreme Energy Decomposition (EED) method to analyze the RRI signal, decomposes an original signal into a plurality of components, namely an extreme component function, and calculates the Energy of each component to obtain the Energy distribution; according to the invention, the signal can be decomposed into different time level signals from high frequency to low frequency according to the fluctuation rule of the RRI signal, and the frequency band is finely divided; the data length obtained by the extremum decomposition on all levels is the same, so that the data length cannot be reduced, and the extremum decomposition method can be used for short-time data analysis, namely, an accurate result can be obtained by analysis with less data quantity; the EED is not susceptible to noise for different levels of component energy analysis.
Drawings
FIG. 1 is a diagram of an original signal in the present invention;
FIG. 2 is a diagram illustrating the envelope of the original signal according to the present invention;
FIG. 3 is a schematic diagram of a subtracted envelope mean signal of an original signal in accordance with the present invention;
FIG. 4 is a schematic diagram of obtaining a first extreme modal function component according to the present invention;
FIG. 5 is a schematic flow chart of the extreme energy decomposition method of the present invention;
FIG. 6 is an EED decomposition diagram of the RRI signal in embodiment 1 of the present invention;
fig. 7 is a diagram illustrating normalized energy distribution in an RRI signal in embodiment 1 of the present invention.
Detailed Description
The technical scheme of the invention is further explained by combining the attached drawings.
As shown in fig. 1, fig. 2, fig. 3, fig. 4 and fig. 5, a method for analyzing heart rate variability signals based on an extremum energy decomposition method of the present invention includes the following steps:
(1) acquiring an ECG signal in an unknown state at a given time and a given sampling frequency, then carrying out denoising pretreatment on the ECG signal, and extracting an RRI signal from the ECG signal to obtain an RRI signal x (t) in the unknown state; the specific method of denoising pretreatment comprises the following steps: because the ECG energy is mainly concentrated in 0-40 Hz, the ECG signal is filtered by a 40Hz zero-phase FIR low-pass filter to eliminate high-frequency noise, and then baseline drift is removed by a median filter;
(2) taking RRI signal x (t) as original signal, where N is 2, the minimum data amount required by original signal x (t)n+1Wherein n is the number of resolved extremum modal function components; all local extreme points of the original signal are solved, and then all the extreme points and all the minimum points of the original signal are connected by a spline curve to form an upper envelope line e respectivelymaxAnd a lower envelope eminObtaining an envelope mean value signal m (t) ═ e of the upper and lower envelope linesmax+emin)/2；
(3) Subtracting the envelope mean signal m (t) from the original signal x (t) to obtain h (t) ═ x (t) -m (t); then judging whether h (t) meets the judgment condition of the extreme value modal function, if not, returning h (t) to the step (2) as an original signal until hk(t) if the condition for determining the extremum mode function is satisfied, c is recorded1(t)＝hk(t) as a first extreme modal function component; the determination condition of the extreme value modal function is as follows: (a) in the whole data sequence, the number of the extreme points is equal to or different from the number of the zero-crossing points by one; (b) at any time, the upper envelope line and the lower envelope line are symmetrical relative to a time axis;
(4) subtracting the first extreme mode function component c from the original signal x (t)1(t) obtaining the residue r1(t)＝x(t)-c1(t) then judging hk(t) whether a stopping criterion is met, and if not, r1(t) as a new original sequence x (t), returning to steps (2) and (3); if h isk(t) when the stopping criterion is met and n is less than 8, returning to the step (1) to obtain the original signal again; if h isk(t) when the stopping criterion is met and n is more than or equal to 8, obtaining 2 nd, 3 rd, … th extreme value modal function components and margin rn(t) the original signal x (t) is then decomposed into n extreme modal function components and a residual, i.e.
Wherein h isk(t) satisfying the stopping criteria the formula:
epsilon represents a screening threshold, and is taken to be 0.2-0.3; the extremum modal decomposition that satisfies the stopping criterion then satisfies the following two conditions: (a) finally obtained extreme value modal function component cn(t) or the remainder rn(t) is less than a predetermined threshold; (b) residual signal rn(t) becomes a monotonous signal from which an extreme mode function signal cannot be extracted again;
(5) opposite polar mode function component ci(t), i is 1, 2, …, n, and the center frequency of each extreme modal function component is obtained by performing spectrum analysis;
(6) n extreme value modal function components obtained by decomposing the original signal x (t) represent components of different frequency bands of the original signal, and then the energy of each component is calculated
Ei＝∫|ci(t)|2dt，i＝1，2，…，n
Normalizing each energy value to obtain a normalized energy distribution vector
pi＝Ei/E，i＝1，2，...，n
Wherein,first component p1Representing the energy of the highest frequency band, representing the proportion of the energy distribution of the signal in the highest frequency band, the last component pnRepresenting the proportion of the energy distribution of the signal in the lowest frequency band range; drawing a normalized energy distribution graph according to the normalized energy distribution vector, wherein the abscissa represents component levels, the ordinate represents normalized energy distribution vector values, the curve represents an average value, and the error bar represents a standard deviation;
(7) selecting a second component p2To the seventh component p7Calculating the energy difference value EDV (p)2+p3+p4)-(p5+p6+p7) When the EDV is less than or equal to a first standard value, the RRI signal is judged to be a normal heart rate variability signal, when the first standard value is less than the EDV and less than a second standard value, the RRI signal is judged to be a suspected abnormal heart rate variability signal, and when the EDV is more than or equal to the second standard value, the RRI signal is judged to be an abnormal heart rate variability signal; wherein the first standard value is-0.15 and the second standard value is 0.08.
The Extreme Energy Decomposition (EED) method adopted by the invention is a method based on the concept of an extreme mode function, wherein the extreme mode function is a signal with single frequency which simultaneously meets the following two conditions:
(a) the number of extreme points (including maxima and minima) and the number of zero-crossing points must be equal or differ by at most one throughout the data sequence;
(b) at any moment, the average value of an upper envelope line formed by the local maximum value points and a lower envelope line formed by the local minimum value points is zero, namely the local upper envelope line and the local lower envelope line are locally symmetrical relative to a time axis;
the above two conditions, condition (a) is similar to the requirement of the gaussian normal stationary process for the traditional narrow band, and condition (b) ensures that the instantaneous frequency calculated by the extreme value modal function has physical significance.
The standard selection of the extreme value modal function decomposition termination is moderate, and the conditions are too strict, so that the last extreme value modal function components lose significance; conditions are too loose, which can result in loss of useful components; in practical application, the number of layers of the extremum modal function components to be decomposed can be set according to requirements, and when the number of layers of decomposition is met, the calculation is terminated.
Example 1
The EED analysis method is used to analyze the energy distribution of the ECG at different levels in healthy persons and CHF patients.
A heart rate variability signal analysis method based on an extreme value energy decomposition method for healthy people comprises the following steps:
(1) acquiring ECG signals of healthy people from an RR interval database nsr2db of physionet; the data comprise 54 healthy people (age 28-76, average 61), then the ECG signals are subjected to denoising pretreatment, RRI signals are extracted from the ECG signals, and RRI signals x (t) are obtained; the specific method of denoising pretreatment comprises the following steps: because the ECG energy is mainly concentrated in 0-40 Hz, the ECG signal is filtered by a 40Hz zero-phase FIR low-pass filter to eliminate high-frequency noise, and then baseline drift is removed by a median filter;
(2) taking RRI signal x (t) as original signal, where N is 2, the minimum data amount required by original signal x (t)n+1＝29Wherein n is the number of the decomposed extreme mode function components, and n is 8; all local extreme points of the original signal are solved, and then all the extreme points and all the minimum points of the original signal are connected by a spline curve to form an upper envelope line e respectivelymaxAnd a lower envelope eminObtaining an envelope mean value signal m (t) ═ e of the upper and lower envelope linesmax+emin)/2；
(3) Subtracting the envelope mean from the original signal x (t)Signal m (t) to yield h (t) x (t) -m (t); then judging whether h (t) meets the judgment condition of the extreme value modal function, if not, returning h (t) to the step (2) as an original signal until hk(t) if the condition for determining the extremum mode function is satisfied, c is recorded1(t)＝hk(t) as a first extreme modal function component; the determination condition of the extreme value modal function is as follows: (a) in the whole data sequence, the number of the extreme points is equal to or different from the number of the zero-crossing points by one; (b) at any time, the upper envelope line and the lower envelope line are symmetrical relative to a time axis;
(5) subtracting the first extreme mode function component c from the original signal x (t)1(t) obtaining the residue r1(t)＝x(t)-c1(t) adding r1(t) as a new original sequence x (t), returning to the steps (2) and (3) to obtain the 2 nd, 3 rd, … th and 8 th extreme value modal function components and the margin r8(t) the original signal x (t) is then decomposed into 8 extreme modal function components and a residual, i.e.
Obtaining an extreme modal decomposition diagram of a healthy person as shown in fig. 6, it can be seen that the component 1 has the highest frequency, the signal fluctuates on the shortest time scale, and the frequency gradually decreases as the serial number of the component increases;
(6) opposite polar mode function component ci(t), i is 1, 2, …, 8, and performing spectrum analysis to obtain the center frequency of each extreme value modal function component, and obtaining a frequency domain analysis result graph;
(7) 8 extreme value modal function components obtained by decomposing the original signal x (t) represent components of different frequency bands of the original signal, and then the energy of each component is calculated
Ei＝∫|ci(t)|2dt，i＝1，2，…，8
Normalizing each energy value to obtain a normalized energy distribution vector
pi＝Ei/E，i＝1，2，...，8
Wherein,first component p1Representing the energy of the highest frequency band, representing the proportion of the energy distribution of the signal in the highest frequency band, the last component pnRepresenting the proportion of the energy distribution of the signal in the lowest frequency range.
(8) Selecting a second component p2To the seventh component p7Calculating the energy difference value EDV (p)2+p3+p4)-(p5+p6+p7)＝-0.2092±0.2940。
A heart rate variability signal analysis method based on an extremum energy decomposition method for a CHF patient comprises the following steps:
(1) acquiring ECG signals from an RR interval database CHF2db of physionet, wherein the chfdb database comprises 29 patients suffering from Congestive Heart Failure (CHF) (age 34-79, average 55), then carrying out denoising preprocessing on the ECG signals, and extracting RRI signals from the preprocessed ECG signals to obtain RRI signals x (t); the specific method of denoising pretreatment comprises the following steps: because the ECG energy is mainly concentrated in 0-40 Hz, the ECG signal is filtered by a 40Hz zero-phase FIR low-pass filter to eliminate high-frequency noise, and then baseline drift is removed by a median filter;
(2) taking RRI signal x (t) as original signal, where N is 2, the minimum data amount required by original signal x (t)n+1＝29Wherein n is the number of the decomposed extreme mode function components, and n is 8; all local extreme points of the original signal are solved, and then all the extreme points and all the minimum points of the original signal are connected by a spline curve to form an upper envelope line e respectivelymaxAnd a lower envelope eminObtaining an envelope mean value signal m (t) ═ e of the upper and lower envelope linesmax+emin)/2；
(3) Subtracting the envelope mean signal m (t) from the original signal x (t) to obtain h (t) ═ x (t) -m (t); then judging whether h (t) meets the judgment condition of the extreme value modal function, if not, returning h (t) to the step (2) as the original signal until the step (2) is finishedhk(t) if the condition for determining the extremum mode function is satisfied, c is recorded1(t)＝hk(t) as a first extreme modal function component; the determination condition of the extreme value modal function is as follows: (a) in the whole data sequence, the number of the extreme points is equal to or different from the number of the zero-crossing points by one; (b) at any time, the upper envelope line and the lower envelope line are symmetrical relative to a time axis;
(4) subtracting the first extreme mode function component c from the original signal x (t)1(t) obtaining the residue r1(t)＝x(t)-c1(t) adding r1(t) as a new original sequence x (t), returning to the steps (2) and (3) to obtain the 2 nd, 3 rd, … th and 8 th extreme value modal function components and the margin r8(t) the original signal x (t) is then decomposed into 8 extreme modal function components and a residual, i.e.
(5) Opposite polar mode function component ci(t), i is 1, 2, …, 8, and obtaining the central frequency of each extreme modal function component by performing spectrum analysis;
(6) 8 extreme value modal function components obtained by decomposing the original signal x (t) represent components of different frequency bands of the original signal, and then the energy of each component is calculated
Ei＝∫|ci(t)|2dt，i＝1，2，…，8
Normalizing each energy value to obtain a normalized energy distribution vector
pi＝Ei/E，i＝1，2，...，8
Wherein,first component p1Representing the energy of the highest frequency band, representing the proportion of the energy distribution of the signal in the highest frequency band, the last component pnRepresenting the proportion of the energy distribution of the signal in the lowest frequency band range; plotting from normalized energy distribution vectors of healthy and CHF patientsPreparing a normalized energy distribution graph, wherein the abscissa represents component levels, the ordinate represents normalized energy distribution vector values, the curve represents an average value, and the error bars represent standard deviations;
(7) selecting a second component p2To the seventh component p7Calculating the energy difference value EDV (p)2+p3+p4)-(p5+p6+p7)＝0.2642±0.4070。
The present invention further analyzes the mean center frequencies of different level components of HRV signals of healthy and CHF patients to obtain the results of Table 1.
TABLE 1 mean center frequencies for different component levels of HRV signals in healthy and CHF patients
It can be concluded that the HRV center frequency gradually decreases as the component level increases.
A typical way of partitioning the frequency domain by a conventional Power Spectral Density (PSD) method is: HF (0.15-0.4 Hz), LF (0.04-0.15 Hz), and VLF (0.0033-0.04 Hz). The EED method of the invention has the frequency higher than HF for the component level 1; level 2 is in the frequency range of HF; levels 3, 4 are in the LF range; levels 5-7 are in the frequency range of VLF; level 8 is below VLF. In addition, it can be seen that CHF patients at the same level have a slightly higher frequency than healthy people, reflecting the effect of heart disease on the rhythm of HRV fluctuations, which are faster in the same level.
The invention separates HRV signals of two groups into extreme value modal function component Ci(t), calculating to obtain a normalized energy distribution vector, and drawing an EED graph, as shown in FIG. 7. FIG. 7 is a graph showing EED analysis results of RR interval signals of healthy and CHF patients, wherein the data length is 10000 points, the curve represents the mean value, and the error bars represent the standard deviation. The symbol above the curve represents the two groups of population energies Ttest p<0.01. In the hierarchy selection, hierarchies 1 and above 7 are removed, including margins. Level 1 is susceptible to noise, resulting in large fluctuations in energy,causes the standard deviation of the result to be too large, and the frequency can be as high as more than a few kHz, so that the clear physiological significance is not provided; levels higher than 7 reflect the long-term rhythm of the signal, are very susceptible to the external environment, and have low frequency and unknown physiological significance. In FIG. 7, the normalized energy values for CHF patients are significantly higher at the component low levels (levels 2, 3) than for healthy people, and at the component high levels (levels)>5) In contrast, the opposite occurs, healthy people have higher energy than CHF patients; the energy of healthy people is stable in levels 2-5, when the level is more than 5, the energy slowly rises, and CHF patients quickly fall in levels 2-4 and then tend to be stable; at 4 levels (2, 3, 6, 7), there is a significant difference in energy between the two (p)<0.01). For comparison, the present invention adds the result of the EED analysis of the replacement data (CHF recovery) generated by the method of randomizing the original data as shown in fig. 7, the energy distribution of the replacement data monotonically decreasing with increasing scale, being higher on the small time scale and lower on the long time scale than in CHF patients.
To evaluate the energy distribution difference of the EED curve at the low level and the high level decomposition, the energy difference value EDV of each population is calculated, and the high EDV value represents the higher component low level energy distribution and the lower component high level energy distribution of the RRI signal. EDV values were calculated for healthy people, CHF patients, and their surrogate data, as shown in table 2.
TABLE 2 EDV values for healthy and CHF patients
Represents healthy humans and CHF patients with T-test results p < 0.0001.
And represents the substitute data and the original data T test result p < 0.0001.
As can be seen from table 2, healthy and CHF patients have significant differences in EDV values, which are much higher than those of healthy people. An EDV value of a healthy person less than 0 indicates a higher energy at the high level of the component, indicating a higher intensity of modulation for the high level component. The EDV values of both populations were significantly different from their surrogate data EDV values, with the EDV of human HRV being significantly less than the random sequence.
Heart Rate Variability (HRV) refers to the measurement of the temporal variation between successive cardiac cycles, and specifically, the variation in the difference between the successively occurring normal P-P intervals. However, since P-waves are less pronounced than R-waves or the P-wave tips are sometimes blunt, studies on heart rate variability signals are typically replaced with R-R Interval Signals (RRIs) equal to P-P intervals. Research shows that HRV can be used as noninvasive detection index of activity of the vegetative nervous system, and has important significance in the aspect of judging the prognosis of certain cardiovascular diseases.
In clinical and practical applications, the present invention proposes that 512 (n-8) short-time cardiac signals with consecutive RR intervals can be used for the above described EED analysis, which is effective and has a certain data margin; the above studies show that EED decomposition has good results for 7 component levels (N-7), and the minimum required data size can be N-2n+1＝27+1＝256。
Claims (6)
1. A heart rate variability signal analysis method based on an extreme value energy decomposition method is characterized by comprising the following steps:
(1) acquiring an ECG signal in an unknown state at a given time and a given sampling frequency, then carrying out denoising pretreatment on the ECG signal, and extracting an RRI signal from the ECG signal to obtain an RRI signal x (t) in the unknown state;
(2) taking the RRI signal x (t) as an original signal, calculating all local extreme points of the original signal, and then connecting all the extreme points and all the minimum points of the original signal by adopting a spline curve to respectively form an upper envelope line emaxAnd a lower envelope eminObtaining an envelope mean value signal m (t) ═ e of the upper and lower envelope linesmax+emin)/2；
(3) Subtracting the envelope mean signal m (t) from the original signal x (t) to obtain h (t) ═ x (t) -m (t); then judging whether h (t) meets the judgment condition of the extreme value modal function, if not, returning h (t) to the step (3) as the original signal until hk(t) If the judgment condition of the extreme value modal function is satisfied, recording c1(t)＝hk(t) as a first extreme modal function component;
(4) subtracting the first extreme mode function component c from the original signal x (t)1(t) obtaining the residue r1(t)＝x(t)-c1(t) then judging hk(t) whether a stopping criterion is met, and if not, r1(t) as a new original sequence x (t), returning to steps (2) and (3); if h isk(t) when the stopping criterion is met and n is less than 8, returning to the step (1) to obtain the original signal again; if h isk(t) when the stopping criterion is met and n is more than or equal to 8, obtaining 2 nd, 3 rd, … th extreme value modal function components and margin rn(t) the original signal x (t) is then decomposed into n extreme modal function components and a residual, i.e.
(5) Opposite polar mode function component ci(t), i is 1, 2, …, n, and the center frequency of each extreme modal function component is obtained by performing spectrum analysis;
(6) n extreme value modal function components obtained by decomposing the original signal x (t) represent components of different frequency bands of the original signal, and then the energy of each component is calculated
Ei＝∫|ci(t)|2dt,i＝1，2，…，n
Normalizing each energy value to obtain a normalized energy distribution vector
pi＝Ei/E,i＝1,2,…,n
Wherein,first component p1Representing the energy of the highest frequency band, representing the proportion of the energy distribution of the signal in the highest frequency band, the last component pnRepresenting the proportion of the energy distribution of the signal in the lowest frequency band range; plotting normalized energy scores from normalized energy distribution vectorsLayout, wherein the abscissa represents the component hierarchy, the ordinate represents the normalized energy distribution vector value, the curve represents the mean, and the error bar represents the standard deviation;
(7) selecting a second component p2To the seventh component p7Calculating the energy difference value EDV (p)2+p3+p4)-(p5+p6+p7) And when the EDV is less than or equal to a first standard value, the RRI signal is judged to be a normal heart rate variability signal, when the first standard value is less than or equal to the EDV and less than a second standard value, the RRI signal is judged to be a suspected abnormal heart rate variability signal, and when the EDV is more than or equal to the second standard value, the RRI signal is judged to be an abnormal heart rate variability signal.
2. A method of analyzing a heart rate variability signal based on the extreme energy decomposition method according to claim 1, characterized in that: the minimum data amount N of the original signal x (t) is 2n+1And n is the number of the decomposed extreme mode function components.
3. A method of analyzing a heart rate variability signal based on the extreme energy decomposition method according to claim 1, characterized in that: the specific method of denoising pretreatment in the step (1) is as follows: the ECG signal is filtered through a 40Hz zero phase FIR low pass filter to remove high frequency noise, and then through a median filter to remove baseline wander.
4. A method of analyzing a heart rate variability signal based on the extreme energy decomposition method according to claim 1, characterized in that: the determination conditions of the extreme mode function in the step (3) are as follows: (a) in the whole data sequence, the number of the extreme points is equal to or different from the number of the zero-crossing points by one; (b) and the upper and lower envelopes are symmetric with respect to the time axis at any time.
5. A method of analyzing a heart rate variability signal based on the extreme energy decomposition method according to claim 1, characterized in that: h in the step (4)k(t) satisfying the stopping criteria the formula:
epsilon represents a screening threshold, and is taken to be 0.2-0.3.
6. A method of analyzing a heart rate variability signal based on the extreme energy decomposition method according to claim 1, characterized in that: the first standard value in the step (7) is-0.15, and the second standard value is 0.08.
Priority Applications (3)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201910841605.6A CN110464337B (en) | 2019-09-06 | 2019-09-06 | Heart rate variability signal analysis method based on extreme value energy decomposition method |
| PCT/CN2019/121416 WO2021042591A1 (en) | 2019-09-06 | 2019-11-28 | Heart rate variability signal analysis method based on extremum energy decomposition method |
| AU2019101756A AU2019101756A4 (en) | 2019-09-06 | 2019-11-28 | Method for analyzing heart rate variability signal based on extremum energy decomposition method |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201910841605.6A CN110464337B (en) | 2019-09-06 | 2019-09-06 | Heart rate variability signal analysis method based on extreme value energy decomposition method |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN110464337A true CN110464337A (en) | 2019-11-19 |
| CN110464337B CN110464337B (en) | 2020-09-01 |
Family
ID=68515060
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201910841605.6A Active CN110464337B (en) | 2019-09-06 | 2019-09-06 | Heart rate variability signal analysis method based on extreme value energy decomposition method |
Country Status (3)
| Country | Link |
|---|---|
| CN (1) | CN110464337B (en) |
| AU (1) | AU2019101756A4 (en) |
| WO (1) | WO2021042591A1 (en) |
Cited By (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| WO2021042591A1 (en) * | 2019-09-06 | 2021-03-11 | 江苏华康信息技术有限公司 | Heart rate variability signal analysis method based on extremum energy decomposition method |
| CN115040136A (en) * | 2022-06-10 | 2022-09-13 | 广东粤港澳大湾区国家纳米科技创新研究院 | Electrocardiosignal baseline noise detection method and system |
| CN116548928A (en) * | 2023-07-11 | 2023-08-08 | 西安浩阳志德医疗科技有限公司 | Nursing service system based on internet |
| CN117503153A (en) * | 2024-01-05 | 2024-02-06 | 北华大学 | Patient postoperative rehabilitation evaluation method based on artificial intelligence |
Families Citing this family (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN116738153B (en) * | 2023-08-15 | 2023-10-17 | 中国科学院东北地理与农业生态研究所 | Organic fertilizer utilization effect evaluation method based on spectral analysis |
| CN117390380B (en) * | 2023-12-12 | 2024-02-13 | 泰安金冠宏食品科技有限公司 | Data analysis method in oil-residue separation system |
| CN118430732B (en) * | 2024-07-02 | 2024-09-10 | 深圳爱递医药科技有限公司 | Hemodialysis data intelligent processing system based on smart phone application program |
Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20080269628A1 (en) * | 2007-04-25 | 2008-10-30 | Siemens Medical Solutions Usa, Inc. | Denoising and Artifact Rejection for Cardiac Signal in a Sensis System |
| CN101496716A (en) * | 2009-02-26 | 2009-08-05 | 周洪建 | Measurement method for detecting sleep apnoea with ECG signal |
| CN101642369A (en) * | 2008-08-04 | 2010-02-10 | 南京大学 | Autonomic nervous function biological feedback method and system |
| CN108814579A (en) * | 2018-04-16 | 2018-11-16 | 西安交通大学 | A method of electrocardio, breathing combined calculation heart rate variability based on EMD decomposition |
Family Cites Families (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN110464337B (en) * | 2019-09-06 | 2020-09-01 | 江苏华康信息技术有限公司 | Heart rate variability signal analysis method based on extreme value energy decomposition method |
-
2019
- 2019-09-06 CN CN201910841605.6A patent/CN110464337B/en active Active
- 2019-11-28 AU AU2019101756A patent/AU2019101756A4/en active Active
- 2019-11-28 WO PCT/CN2019/121416 patent/WO2021042591A1/en active Application Filing
Patent Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20080269628A1 (en) * | 2007-04-25 | 2008-10-30 | Siemens Medical Solutions Usa, Inc. | Denoising and Artifact Rejection for Cardiac Signal in a Sensis System |
| CN101642369A (en) * | 2008-08-04 | 2010-02-10 | 南京大学 | Autonomic nervous function biological feedback method and system |
| CN101496716A (en) * | 2009-02-26 | 2009-08-05 | 周洪建 | Measurement method for detecting sleep apnoea with ECG signal |
| CN108814579A (en) * | 2018-04-16 | 2018-11-16 | 西安交通大学 | A method of electrocardio, breathing combined calculation heart rate variability based on EMD decomposition |
Non-Patent Citations (2)
| Title |
|---|
| 庞春颖等: "一种基于改进经验模态分解的癫痫脑电识别新方法", 《中国生物医学工程学报》 * |
| 董红生等: "基于Hilbert谱的心率变异信号时频分析方法", 《仪器仪表学报》 * |
Cited By (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| WO2021042591A1 (en) * | 2019-09-06 | 2021-03-11 | 江苏华康信息技术有限公司 | Heart rate variability signal analysis method based on extremum energy decomposition method |
| CN115040136A (en) * | 2022-06-10 | 2022-09-13 | 广东粤港澳大湾区国家纳米科技创新研究院 | Electrocardiosignal baseline noise detection method and system |
| CN116548928A (en) * | 2023-07-11 | 2023-08-08 | 西安浩阳志德医疗科技有限公司 | Nursing service system based on internet |
| CN116548928B (en) * | 2023-07-11 | 2023-09-08 | 西安浩阳志德医疗科技有限公司 | Nursing service system based on internet |
| CN117503153A (en) * | 2024-01-05 | 2024-02-06 | 北华大学 | Patient postoperative rehabilitation evaluation method based on artificial intelligence |
| CN117503153B (en) * | 2024-01-05 | 2024-03-15 | 北华大学 | Patient postoperative rehabilitation evaluation method based on artificial intelligence |
Also Published As
| Publication number | Publication date |
|---|---|
| CN110464337B (en) | 2020-09-01 |
| AU2019101756A4 (en) | 2020-11-12 |
| WO2021042591A1 (en) | 2021-03-11 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN110464337B (en) | Heart rate variability signal analysis method based on extreme value energy decomposition method | |
| CN110558973B (en) | Computer equipment for executing electrocardiogram signal quantitative analysis method based on extreme value energy decomposition method | |
| Amann et al. | Reliability of old and new ventricular fibrillation detection algorithms for automated external defibrillators | |
| Salo et al. | Ectopic beats in heart rate variability analysis: effects of editing on time and frequency domain measures | |
| Goovaerts et al. | A machine-learning approach for detection and quantification of QRS fragmentation | |
| Mandala et al. | ECG parameters for malignant ventricular arrhythmias: a comprehensive review | |
| EP0801538A1 (en) | Method and apparatus for assessing cardiovascular risk | |
| Miličević | Low to high frequency ratio of heart rate variability spectra fails to describe sympatho-vagal balance in cardiac patients | |
| Costin et al. | Atrial fibrillation onset prediction using variability of ECG signals | |
| CN110558959B (en) | HRV signal analysis method for meditation training based on extreme value energy decomposition method | |
| Rankawat et al. | Robust heart rate estimation from multimodal physiological signals using beat signal quality index based majority voting fusion method | |
| Long et al. | Time-frequency analysis of heart rate variability for sleep and wake classification | |
| CN110558974B (en) | Electrocardiogram signal analysis method based on extreme value energy decomposition method | |
| Riasi et al. | Prediction of ventricular tachycardia using morphological features of ECG signal | |
| Khavas et al. | Robust heartbeat detection using multimodal recordings and ECG quality assessment with signal amplitudes dispersion | |
| Kuntamalla et al. | The effect of aging on nonlinearity and stochastic nature of heart rate variability signal computed using delay vector variance method | |
| Belinchon et al. | How reproducible are heart rate variability indices along the time to predict cardiovascular events in hypertensive patients? | |
| Xia et al. | Influence of recognition errors of computerised analysis of 24-hour electrocardiograms on the measurement of spectral components of heart rate variability | |
| Ruiz et al. | Thorough assessment of a p-wave delineation algorithm through the use of diverse electrocardiographic databases | |
| WO2005006209A1 (en) | Method and apparatus for detecting paroxysmal atria fibrillation | |
| Patil et al. | Discrimination between atrial fibrillation (AF) & normal sinus rhythm (NSR) using linear parameters | |
| CN117137497B (en) | Cardiac rhythm prediction method, defibrillation control method and corresponding devices | |
| Raju et al. | Automatic detection and classification of cardiac arrhythmia using neural network | |
| Ouyang et al. | Data mining analysis of the influences of electrocardiogram P-wave morphology parameters on atrial fibrillation | |
| Saraiva et al. | Classification of cardiovascular signals |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| CB03 | Change of inventor or designer information |
Inventor after: Ning Xinbao Inventor after: Zeng Peng Inventor after: Zhou Zuojian Inventor after: Jiang Xiaodong Inventor after: Liu Guo Inventor after: Wang Binbin Inventor after: Wang Hua Inventor before: Zhou Zuojian Inventor before: Ning Xinbao Inventor before: Wang Binbin Inventor before: Jiang Xiaodong Inventor before: Wang Hua |
|
| CB03 | Change of inventor or designer information | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |