CN110464337B - Heart rate variability signal analysis method based on extreme value energy decomposition method - Google Patents

Abstract

The invention discloses a heart rate variability signal analysis method based on an extreme value energy decomposition method, which comprises the steps of obtaining an ECG signal of an unknown state at a given time and a given sampling frequency, and obtaining an RRI signal x (t) after denoising; taking an RRI signal x (t) as an original signal, decomposing the original signal x (t) into n extreme value modal function components and a margin, decomposing the original signal x (t) into n extreme value modal function components representing components of different frequency bands of the original signal, and judging whether the RRI signal is an abnormal heart rate variability signal or not according to the n extreme value modal function components. The invention adopts an extreme value energy decomposition method to analyze the RRI signal, decomposes the original signal into a plurality of components, namely an extreme value component function, and calculates the energy of each component to obtain the energy distribution of the components.

Description

Heart rate variability signal analysis method based on extreme value energy decomposition method

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 e_maxAnd a lower envelope e_minObtaining an envelope mean value signal m (t) ═ e of the upper and lower envelope lines_max+e_min)/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 h_k(t) if the condition for determining the extremum mode function is satisfied, c is recorded₁(t)＝h_k(t) as a first extreme modal function component;

(4) subtracting the first extreme mode function component c from the original signal x (t)₁(t) obtaining the residue r₁(t)＝x(t)-c₁(t) then judging h_k(t) whether a stopping criterion is met, and if not, r₁(t) as a new original sequence x (t), returning to steps (2) and (3); if h is_k(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 is_k(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 r_n(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 c_i(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

E_i＝∫|c_i(t)|²dt，i＝1，2，…，n

Normalizing each energy value to obtain a normalized energy distribution vector

p_i＝E_i/E，i＝1，2，...，n

Wherein the content of the first and second substances,

first component p₁Representing 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 p_nRepresenting 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 p₂To the seventh component p₇Calculating the energy difference value EDV (p)₂+p₃+p₄)-(p₅+p₆+p₇) 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 2ⁿ⁺¹And n is the number of the decomposed extreme mode function components.

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:

and representing a screening threshold, and taking the value between 0.2 and 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)ⁿ⁺¹Wherein 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 respectively_maxAnd a lower envelope e_minObtaining an envelope mean value signal m (t) ═ e of the upper and lower envelope lines_max+e_min)/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 h_k(t) if the condition for determining the extremum mode function is satisfied, c is recorded₁(t)＝h_k(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)₁(t) obtaining the residue r₁(t)＝x(t)-c₁(t) then judging h_k(t) whether a stopping criterion is met, and if not, r₁(t) as a new original sequence x (t), returning to steps (2) and (3); if h is_k(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 is_k(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 r_n(t) the original signal x (t) is then decomposed into n extreme modal function components and a residual, i.e.

Wherein h is_k(t) satisfying the stopping criteria the formula:

representing a screening threshold, and taking the value between 0.2 and 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 c_n(t) or the remainder r_n(t) is less than a predetermined threshold; (b) residual signal r_n(t) becomes a monotonous signal from which an extreme mode function signal cannot be extracted again;

(5) opposite polar mode function component c_i(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

E_i＝∫|c_i(t)|²dt，i＝1，2，…，n

Normalizing each energy value to obtain a normalized energy distribution vector

p_i＝E_i/E，i＝1，2，...，n

Wherein the content of the first and second substances,

first component p₁Representing 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 p_nRepresenting 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 p₂To the seventh component p₇Calculating the energy difference value EDV (p)₂+p₃+p₄)-(p₅+p₆+p₇) 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)ⁿ⁺¹＝2⁹Wherein 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 respectively_maxAnd a lower envelope e_minObtaining an envelope mean value signal m (t) ═ e of the upper and lower envelope lines_max+e_min)/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 h_k(t) if the condition for determining the extremum mode function is satisfied, c is recorded₁(t)＝h_k(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, poleThe number of the value 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)₁(t) obtaining the residue r₁(t)＝x(t)-c₁(t) adding r₁(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 r₈(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 c_i(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

E_i＝∫|c_i(t)|²dt，i＝1，2，…，8

Normalizing each energy value to obtain a normalized energy distribution vector

p_i＝E_i/E，i＝1，2，...，8

Wherein the content of the first and second substances,

first component p₁Representing 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 p_nRepresenting the proportion of the energy distribution of the signal in the lowest frequency range.

(8) Selecting a second component p₂To the seventh component p₇Calculating the energy difference value EDV (p)₂+p₃+p₄)-(p₅+p₆+p₇)＝-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)ⁿ⁺¹＝2⁹Wherein 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 respectively_maxAnd a lower envelope e_minObtaining an envelope mean value signal m (t) ═ e of the upper and lower envelope lines_max+e_min)/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 h_k(t) if the condition for determining the extremum mode function is satisfied, c is recorded₁(t)＝h_k(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) will be provided withOriginal signal x (t) minus first extreme mode function component c₁(t) obtaining the residue r₁(t)＝x(t)-c₁(t) adding r₁(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 r₈(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 c_i(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

E_i＝∫|c_i(t)|²dt，i＝1，2，…，8

Normalizing each energy value to obtain a normalized energy distribution vector

p_i＝E_i/E，i＝1，2，...，8

Wherein the content of the first and second substances,

first component p₁Representing 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 p_nRepresenting the proportion of the energy distribution of the signal in the lowest frequency band range; drawing a normalized energy distribution map according to the normalized energy distribution vectors of healthy people and CHF patients, wherein the abscissa represents the component level, the ordinate represents the normalized energy distribution vector value, the curve represents the average value, and the error bar represents the standard deviation;

(7) selecting a second component p₂To the seventh component p₇Calculating the energy difference value EDV (p)₂+p₃+p₄)-(p₅+p₆+p₇)＝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 C_i(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 easily affected by noise, so that large energy fluctuation is caused, the standard deviation of results is overlarge, and the frequency of the level 1 can be as high as more than several kHz, so that the level 1 has no clear physiological significance; 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 the above, the opposite change occursThe energy of healthy people is higher than that of 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-2ⁿ⁺¹＝2⁷⁺¹＝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 e_maxAnd a lower envelope e_minObtaining an envelope mean value signal m (t) ═ e of the upper and lower envelope lines_max+e_min)/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 h_k(t) if the condition for determining the extremum mode function is satisfied, c is recorded₁(t)＝h_k(t) as a first extreme modal function component;

(4) subtracting the first extreme mode function component c from the original signal x (t)₁(t) obtaining the residue r₁(t)＝x(t)-c₁(t) then judging h_k(t) whether a stopping criterion is met, and if not, r₁(t) As a new original sequence x (t), returning to the steps (2) and (3); if h is_k(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 is_k(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 r_n(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 c_i(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

E_i＝∫|c_i(t)|²dt,i＝1，2，…，n

Normalizing each energy value to obtain a normalized energy distribution vector

p_i＝E_i/E,i＝1,2,…,n

Wherein the content of the first and second substances,

first component p₁Representing 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 p_nRepresenting 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 p₂To the seventh component p₇Calculating the energy difference value EDV (p)₂+p₃+p₄)-(p₅+p₆+p₇) When EDV is less than or equal to firstAnd when the EDV is larger 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 2ⁿ⁺¹And 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:

and representing a screening threshold, and taking the value between 0.2 and 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.

Images