CN108836305A - A kind of ECG feature extracting method of fusion Butterworth filtering and wavelet transformation - Google Patents

Abstract

An ECG feature extraction method combining Butterworth filter and wavelet transform belongs to the field of medical computer technology. The core idea is: to denoise the original ECG signal power frequency interference; at the same time, use the curve fitting method to remove the baseline drift; The Q wave, R wave, and S wave are positioned separately, compared with the information marked by experts, and the threshold, scale, and translation parameters are adjusted respectively to obtain the best positioning results. Based on the above results, for Q wave, R wave and S wave, the differential threshold or wavelet transform method is selected for positioning. The invention can perform more accurate denoising processing on the original ECG signal, and can realize the precise positioning of the QRS wave, so as to provide a basis for the digitized analysis of the electrocardiogram.

Description

An ECG Feature Extraction Method Combining Butterworth Filter and Wavelet Transform

technical field

The invention relates to a method for extracting features of electrocardiographic signals, in particular to a method for extracting ECG features combined with Butterworth filter and wavelet transform, and belongs to the technical field of medical computers.

Background technique

For a long time, cardiovascular disease has seriously threatened human health and life because of its high morbidity and high mortality. The World Health Organization (WHO) estimates that 36 million people die each year from non-communicable diseases such as cardiovascular disease, diabetes, respiratory system disease and malignant tumors in the world, accounting for 2/3 of the total global deaths. By 2020, The number will climb to 44 million.

The "China Cardiovascular Disease Report (2015)" shows that in 2014, the mortality rate of cardiovascular disease in China still ranked first in the composition of disease deaths, higher than that of tumors and other diseases. In 2014, the mortality rate of cardiovascular disease in rural areas was 295.63/100,000, of which the mortality rate of heart disease was 143.72/100,000, and the mortality rate of cerebrovascular disease was 151.91/100,000 (cerebral hemorrhage 74.51/100,000, cerebral infarction 45.30/100,000) The mortality rate of cardiovascular diseases in the city is 261.99/100,000, of which the mortality rate of heart disease is 136.21/100,000, and the mortality rate of cerebrovascular disease is 125.78/100,000 (cerebral hemorrhage 52.25/100,000, cerebral infarction 41.99/100,000). With the acceleration of the aging of Chinese society, these data are increasing rapidly, and cardiovascular disease has become one of the major public health problems in our country. Therefore, the prevention, diagnosis and treatment of heart diseases are particularly important, which is a huge challenge and research hotspot for medical personnel today.

Electrocardiogram (ECG) records the current generated by muscle cells during the contraction and relaxation of each motor cycle through electrodes placed on different parts of the human body surface, and also records the change trend of the current with the direction and strength of the cardiac cycle And draw the ECG change curve in real time to judge the electrical activity of the heart. It can determine whether the myocardial cells are dead, whether there is ischemia, whether there is abnormal impulse conduction, and record the characteristic ECG changes with clinical significance. In addition, there are also characteristic changes in the electrocardiogram of certain electrolyte disorders. For example, hyperkalemia will lead to high peak T waves, shortened QT interval, and low blood calcium will cause ST segment prolongation, etc. The feature extraction of ECG signals is of great significance to clinical research.

However, in the detection process of the ECG signal, there are factors such as power frequency interference and environmental noise, which affect the accuracy of the measurement. And at present, because the accuracy and reliability of the intelligent analysis of the electrocardiogram cannot be compared with the analysis of cardiologists, so its proportion in actual clinical application is not large.

Contents of the invention

The purpose of the present invention is to further improve the acquisition accuracy of existing electrocardiographic signals and the precise detection of each QRS band, and proposes an ECG feature extraction method combining Butterworth filter and wavelet transform.

The core idea of the present invention is: to denoise the original ECG signal power frequency interference; at the same time use the way of curve fitting to remove the baseline drift; Methods The Q wave, R wave and S wave were respectively positioned, compared with the information marked by experts, and the threshold, scale and translation parameters were adjusted respectively to obtain the best positioning results.

An ECG feature extraction method that combines Butterworth filter and wavelet transform, comprising the steps of:

Step A, performing spectrum analysis on the original ECG signal to obtain the highest frequency and the lowest frequency, and then using the Butterworth filter to filter the original ECG signal to obtain the filtered ECG ₁ signal; Step A specifically includes the following sub-steps:

SA.1 Use the discrete Fourier transform method to analyze the spectrum of the ECG signal through the formula (1) to obtain the ECG signal spectrum X(e ^jω ), and further obtain the highest frequency f ₁ of the ECG signal, power frequency interference and other The lowest frequency f ₂ of the noise interference frequency band;

Wherein, X(e ^jω ) is the frequency spectrum of the ECG signal, ω is the frequency of the ECG signal, and ECG(n) is the nth amplitude of the ECG signal;

In the obtained X(e ^jω ), the highest frequency of the low frequency band is the highest frequency f ₁ of the ECG signal; the lowest frequency of the high frequency band is the lowest frequency f ₂ of the power frequency interference and other noise interference frequency bands;

SA.2 Filter the ECG signal through the Butterworth filter to obtain the filtered ECG signal ECG ₁ , specifically carry out the Butterworth filter through the formula (2):

Among them, N in the formula (2) represents an integer set, |H(ω)| is the modulus of the Butterworth filter, |H(ω)| ² is the square of the modulus; ω _c (2πf ₁ <ω _c <2πf ₂ ) is the cut-off frequency of the Butterworth filter, and p(p∈N) is the order;

Step B. Carry out curve fitting to the ECG ₁ output in step A, obtain polynomial coefficients and remove the baseline drift, obtain the signal after removing the baseline drift, ECG ₂ ;

Step B includes the following sub-steps:

SB.1 Carry out curve fitting based on formula (3):

Among them, y is the polynomial used for fitting, the order of the polynomial is q(q∈N), a ₀ , a ₁ ,..., a _p are the coefficients of the q-order polynomial, x is the sampling time series, and x ⁱ is x to the ith power;

SB.2 Calculate the polynomial coefficient in formula (3) based on formula (4);

in, is a vector of polynomial coefficients, is the magnitude vector of the ECG signal, C=[1 xx ² ... x ^q ] ^T , C ^T is the transpose of C;

SB.3 Based on the formula (5), obtain the signal after removing the baseline drift;

ECG ₂ =ECG ₁ -y (5)

In formula (5), ECG ₂ is the signal after baseline drift removal, ECG ₁ is the signal after filtering; y is the output of formula (3);

Step C, the ECG ₂ signal output in step B uses the method of differential threshold to first locate the R wave, and then perform the positioning of the Q wave and S wave, and output the positions of the R wave, Q wave and S wave respectively;

Step C includes the following sub-steps:

SC.1 Calculate the first-order difference S(n) and the second-order difference W(n) of the ECG ₂ signal sequentially based on formula (6) and formula (7);

S(n)=ECG ₂ (n)-ECG ₂ (n-1) (6)

W(n)=S(n)-S(n-1) (7)

Wherein, S(n) is the nth first-order difference, ECG ₂ (n) is the nth amplitude of the ECG ₂ signal, ECG ₂ (n-1) is the n-1th amplitude of the ECG ₂ signal, W(n) is the nth second-order difference, S(n-1) is the n-1th first-order difference;

SC.2 use the enumeration method to determine the maximum value max 1 of S(n) and the maximum value max ₂ of W(n) output by step SC.1 _;

SC.3 Based on the formula (8), set the threshold Y ₁ :

Y ₁ ＝m ₁ ·max ₁ (0<m ₁ <1) (8)

In the formula (8), the threshold threshold is Y ₁ , and m ₁ is an artificially set coefficient, and its range is (0,1);

SC.4 Judge whether there are consecutive k ₁ (k ₁ ∈ N) points in S(n) exceeding Y ₁ and ECG ₂ (n) > 0, and perform the following operations according to the judgment result:

SC.4A If there are k ₁ (k ₁ ∈ N) consecutive points exceeding Y ₁ in S(n), and ECG ₂ (n) > 0, this point n is the R peak of this range, and its position is r , skip to SC.5;

Wherein, k ₁ is an integer greater than or equal to 2;

SC.4B If there are no consecutive k ₁ (k ₁ ∈ N) points in S(n) exceeding Y ₁ , then k ₁ =k ₁ -1, skip to step SC.4;

SC.5 Based on formulas (9) and (10), set the starting threshold Y ₂ and the ending threshold Y ₃ :

Y ₂ =m ₂ ·max ₂ (0<m ₂ <1) (9)

Y ₃ =m ₃ ·max ₂ (0<m ₃ <1) (10)

In formulas (9) and (10), the thresholds are Y ₂ and Y ₃ respectively, and m ₂ and m ₃ are artificially set coefficients;

SC.6 Search forward from each R wave, judge whether there are consecutive k ₂ (k ₂ ∈ N) points in W(n) smaller than Y ₂ , and perform the following operations according to the judgment result:

SC.6A If there are k ₂ (k ₂ ∈ N) consecutive points in W(n) smaller than Y ₂ , the point n is the starting point s of the QRS wave in this range, skip to SC.5;

Wherein, k ₂ is an integer greater than or equal to 2;

SC.6B If there are no consecutive k ₂ (k ₂ ∈N) points in W(n) smaller than Y ₂ , then k ₂ =k ₂ -1, skip to step SC.6;

SC.7 Search backward from each R wave, judge whether there are consecutive k ₃ (k ₃ ∈ N) points in W(n) smaller than Y ₃ , and perform the following operations according to the judgment result:

Wherein, k ₃ is an integer greater than or equal to 2;

SC.7A If there are k ₃ (k ₃ ∈ N) consecutive points in W(n) smaller than Y ₃ , this point n is the end point e of the QRS wave in this range, skip to SD.8;

SC.7B If there are no consecutive k ₃ (k ₃ ∈N) points in W(n) smaller than Y ₃ , then k ₃ =k ₃ -1, skip to step SC.7;

SC.8 Based on the formula (11), search for the minimum amplitude point Q within the range of the starting point s of the Q wave, R wave and S wave and the position r of the R wave:

Q=min(ECG ₂ (n)),(s<n<r) (11)

In formula (11), min(ECG ₂ (n)) is a function to search for the minimum value of ECG ₂ (n);

Based on the formula (12), the amplitude minimum point S is searched within the range of the R wave and the terminal point e:

S=min(ECG ₂ (n)),(r<n<e) (12)

SC.9 Compare the positions of R wave, Q wave and S wave determined by step SC.4 and step SC.8 with the positions marked by experts, and obtain R wave, Q wave and S wave respectively based on formula (13) The detection accuracy of p ₁₁ , p ₁₂ , p ₁₃ :

In the formula (13), the value of i is 1, 2, 3, N ₁₁ , N ₁₂ , N ₁₃ are respectively the number of R wave, Q wave and S wave in the ECG signal, N ₁₁ ^* , N ₁₂ ^* , N ₁₃ ^* are respectively the number of correctly detected R wave, Q wave and S wave;

Step D, using the wavelet transform method to obtain wavelet decomposition diagrams of different levels of the ECG ₂ signal processed in step B, and then output the positions of Q, R, and S waves by detecting the position and amplitude of the wavelet transform coefficient modulus maximum; Step D comprises the following steps again:

SD.1 Based on the formulas (14) and (15), the biorthogonal spline wavelet is used to perform wavelet transformation on the ECG signal, and the wavelet decomposition diagrams at different levels are obtained;

Among them, x ⁰ (n)=ECG ₂ (n), which is the processed ECG signal; Z is a natural number set (the same below); x ^(j) (n) is the high-frequency signal of x ^j-1 (n) The binary wavelet transform on the j scale, d ^(j) (n) is the binary wavelet transform of the low-frequency signal of x ^j-1 (n) on the j scale, h _j-1 (k), g _{j- 1} (k) is respectively the digital low-pass filter and the high-pass filter coefficient of the k group orthogonality of setting;

SD.2 Based on the wavelet decomposition diagrams at different scales obtained in SD.1, for waveforms at the scale of c(c∈N), detect the amplitude maximum point M and the minimum value point L; detect the amplitude between two points is the point of p(p<ε), where ε is the value whose amplitude is close to zero, and p is the position of the R wave of the original ECG signal;

SD.3 Based on the R wave position p obtained in SD.2 and the formula (16), search for the modulus minimum at the front position of the R wave for the waveform at the b (b∈N) scale, which is the Q wave.

Q=x ^(b) (n)<x ^(b) (n-1), x ^(b) (n)<x ^(b) (n+1), (n<p) (16)

Wherein, x ^(b) (n) is the magnitude of the nth sequence of the ECG signal at the b scale;

Search for the modulus minimum at the position after the R wave, which is the S wave,

S=x ^(b) (n)<x ^(b) (n-1), x ^(b) (n)<x ^(b) (n+1), (n>p) (17)

Wherein, x ^(b) (n) is the magnitude of the nth sequence of the ECG signal at the b scale;

SD.4 Compare the positions of R wave, Q wave and S wave detected by step SD.2 and step SD.3 with the positions marked by experts, and obtain R wave, Q wave and S wave respectively based on formula (18) The detection accuracy of p ₂₁ , p ₂₂ , p ₂₃ :

In the formula (18), the value of i is 1, 2, 3, N ₂₁ , N ₂₂ , N ₂₃ are respectively the number of R wave, Q wave and S wave in the ECG signal, N ₂₁ ^* , N ₂₂ ^* , N ₂₃ ^* are respectively the number of correctly detected R wave, Q wave and S wave;

Wherein, step D and step C can also exchange order;

Step E, repeatedly adjust the artificially set parameters of the difference threshold method through the R wave, Q wave and S wave marked by experts, and repeatedly adjust the wavelet transform hierarchy until the accuracy converges; step E also includes the following steps :

SE.1 Adjust the values of m ₁ , m ₂ , and m ₃ based on the R wave, Q wave, and S wave obtained in step C and the points marked by experts repeatedly, and decide whether to skip to step C, specifically:

SE.1A Based on the R wave obtained in step C, compare it with the R wave marked by experts repeatedly. If the accuracy p ₁₁ improves, adjust the value of m ₁ and skip to step C; if the accuracy p ₁₁ does not improve, Then skip to step SE.1B;

Step SE.1B is based on the Q wave obtained in step C, compared with the Q wave marked by repeated consultation with experts,

If the accuracy p ₁₂ is improved, then adjust the value of m ₂ and skip to step C; if the accuracy p ₁₂ is not improved, then skip to step SE.1C;

Step SE.1C is based on the S wave obtained in step C, compared with the S wave marked by experts repeatedly,

If the accuracy p ₁₃ is improved, then adjust the value of m ₃ and skip to step C; if the accuracy p ₁₃ is not improved, then skip to step SE.2;

SE.2 Based on the R wave, Q wave and S wave obtained in step D and the points marked by experts repeatedly, adjust the values of b and c to decide whether to skip to step D, specifically:

SE.2A Based on the R wave obtained in step D, compare it with the R wave marked by experts repeatedly. If the accuracy p ₂₁ improves, adjust the value of b and skip to step D; if the accuracy p ₂₁ does not improve, then Skip to step SE.2B;

Step SE.2B is based on the Q wave and S wave obtained in step D, compared with the Q wave and S wave marked by experts repeatedly,

If both the accuracy p ₂₂ and p ₂₃ are improved, then adjust the value of c and skip to step D; if at least one of the accuracy p ₂₂ and p ₂₃ is not improved, then skip to step F.

Step F, comparing the detection accuracy of step C and step D respectively, and selecting a detection method with high accuracy for R wave, Q wave and S wave respectively, step F includes the following steps:

SF.1 Compare the size of p ₁₁ and p ₂₁ , if p ₁₁ ≥ p ₂₁ , then skip to step C to extract R wave, if p ₁₁ <p ₂₁ , then skip to step D to extract R wave;

SF.2 Compare the size of p ₁₂ and p ₂₂ , if p ₁₂ ≥ p ₂₂ , skip to step C to extract Q wave, if p ₁₂ < p ₂₂ , skip to step D to extract Q wave;

SF.3 Compare the size of p ₁₃ and p ₂₃ , if p ₁₃ ≥ p ₂₃ , then skip to step C to extract S wave, if p ₁₃ <p ₂₃ , then skip to step D to extract S wave, This method ends.

Beneficial effect

An ECG feature extraction method fused with Butterworth filter and wavelet transform of the present invention has the following beneficial effects compared with the prior art: the original ECG signal can be denoised more accurately, and the QRS wave can be realized Accurate positioning provides the basis for digital analysis of ECG.

Description of drawings

Fig. 1 is a flow chart of extracting Q wave, R wave and S wave of an ECG feature extraction method that combines Butterworth filter and wavelet transform according to the present invention.

Detailed ways

Hereinafter, the present invention will be described in detail with reference to the accompanying drawings and examples.

Example 1

The present invention proposes an ECG signal feature extraction method that combines Butterworth filter and wavelet transform, as shown in FIG. 1 . Specifically include the following steps:

Step 1. Perform spectrum analysis on the original ECG signal to obtain the highest frequency and the lowest frequency, and then use the Butterworth filter to filter the original ECG signal to obtain the filtered ECG ₁ signal; Step 1 specifically includes the following sub-steps:

S1.1 Use the discrete Fourier transform method to analyze the spectrum of the ECG signal through the formula (19), and obtain the highest frequency and the lowest frequency of the ECG signal spectrum, X(e ^jω ):

Wherein, X(e ^jω ) is the frequency spectrum of the ECG signal, ω is the frequency of the ECG signal, and ECG(n) is the nth amplitude of the ECG signal;

In the obtained X(e ^jω ), the highest frequency of the low frequency band is the highest frequency f ₁ of the ECG signal; the lowest frequency of the high frequency band is the lowest frequency f ₂ of the power frequency interference and other noise interference frequency bands;

S1.2 Filter the ECG signal through a Butterworth filter to obtain a filtered ECG signal ECG ₁ , and perform Butterworth filtering based on formula (20);

Wherein, N in the formula (20) represents an integer set, |H(ω)| is the modulus of the Butterworth filter, |H(ω)| ² is the square of the modulus; ω _c =2π×50Hz is the Butterworth The cutoff frequency of the filter, p=2 is the order;

Step 2. Carry out curve fitting to the ECG ₁ output in step SA.2, obtain polynomial coefficients and remove the baseline drift, and obtain the signal ECG ₂ after the baseline drift is removed;

Step 2 includes the following sub-steps:

S2.1 Carry out curve fitting based on formula (21);

Among them, y is the polynomial used for fitting, the order of the polynomial is q=3, a ₀ , a ₁ , a ₂ , and a ₃ are the coefficients of the q-order polynomial, x is the sampling time series, and x ⁱ is the ith of x power;

S2.2 Calculate the polynomial coefficient in formula (3) based on formula (22);

in, is a vector of polynomial coefficients, is the magnitude vector of the ECG signal, C=[1 xx ² x ³ ] ^T , C ^T is the transpose of C;

S2.3 Based on formula (23), obtain the signal after removing the baseline drift;

ECG ₂ = ECG ₁ -y (23)

In formula (23), ECG ₂ is the signal after baseline drift removal, ECG ₁ is the signal after filtering; y is the output of formula (21);

Step 3. The ECG ₂ signal output in step 2 is used to locate the R wave first, then the Q wave and S wave, and output the positions of the R wave, Q wave, and S wave respectively by using the differential threshold method;

Step 3 includes the following sub-steps:

S3.1 Sequentially calculate the first-order difference S(n) and the second-order difference W(n) of the ECG ₂ signal based on formula (24) and formula (25);

S(n)=ECG ₂ (n)-ECG ₂ (n-1) (24)

W(n)=S(n)-S(n-1) (25)

Wherein, S(n) is the nth first-order difference, ECG ₂ (n) is the nth amplitude of the ECG ₂ signal, ECG ₂ (n-1) is the n-1th amplitude of the ECG ₂ signal, W(n) is the nth second-order difference, S(n-1) is the n-1th first-order difference;

S3.2 use enumeration method to determine the maximum value max ₁ of S(n) output by step S3.1 and the maximum value max ₂ of W(n);

S3.3 Based on the formula (26), set the threshold Y ₁ :

Y ₁ =m ₁ ·max ₁ (0<m ₁ <1) (26)

In formula (26), the threshold is Y ₁ , m ₁ =0.7;

S3.4 Judge whether there are consecutive k ₁ (k ₁ ∈ N) points in S(n) exceeding Y ₁ and ECG ₂ (n) > 0, and perform the following operations according to the judgment result:

S3.4A If there are k ₁ (k ₁ ∈ N) consecutive points in S(n) exceeding Y ₁ , and ECG ₂ (n) > 0, the point n is the R peak of this range, and its position is r , skip to S3.5;

Wherein, k ₁ is an integer greater than or equal to 2;

S3.4B If there are no consecutive k ₁ (k ₁ ∈ N) points in S(n) exceeding Y ₁ , then k ₁ =k ₁ -1, skip to step SC.4;

S3.5 Based on formulas (27) and (28), set the starting threshold threshold Y ₂ and the ending threshold threshold Y ₃ :

Y ₂ =m ₂ ·max ₂ (0<m ₂ <1) (27)

Y ₃ =m ₃ ·max ₂ (0<m ₃ <1) (28)

In the formulas (27) and (28), the thresholds are respectively Y ₂ and Y ₃ , m ₂ =0.5, m ₃ =0.6;

S3.6 Search forward from each R wave, judge whether there are consecutive k ₂ (k ₂ ∈ N) points in W(n) smaller than Y ₂ , and perform the following operations according to the judgment result:

S3.6A If there are consecutive k ₂ = 5 points in W(n) smaller than Y ₂ , the point n is the starting point s of the QRS wave in this range, skip to S3.5;

Wherein, k ₂ is an integer greater than or equal to 2;

S3.6B If there are no consecutive k ₂ =5 points in W(n) smaller than Y ₂ , then k ₂ =k ₂ -1, skip to step S3.6;

S3.7 Search backward from each R wave, judge whether there are consecutive k ₃ =5 points in W(n) smaller than Y ₃ , and perform the following operations according to the judgment result:

Wherein, k ₃ is an integer greater than or equal to 2;

S3.7A If there are k ₃ = 5 consecutive points in W(n) smaller than Y ₃ , the point n is the end point e of the QRS wave in this range, skip to S4.8;

S3.7B If there are no consecutive k ₃ =5 points in W(n) smaller than Y ₃ , then k ₃ =k ₃ -1, skip to step S3.7;

S3.8 Based on the formula (29), search for the minimum amplitude point Q within the range of the starting point s of the Q wave, R wave and S wave and the position r of the R wave:

Q=min(ECG ₂ (n)),(s<n<r) (29)

In formula (29), min(ECG ₂ (n)) is a function to search for the minimum value of ECG ₂ (n);

Based on the formula (30), the amplitude minimum point S is searched within the range of the R wave and the terminal point e:

S=min(ECG ₂ (n)),(r<n<e) (30)

S3.9 Compare the positions of R wave, Q wave and S wave determined by step S3.4 and step S3.8 with the positions marked by experts, and obtain R wave, Q wave and S wave respectively based on formula (13) The detection accuracy of p ₁₁ , p ₁₂ , p ₁₃ :

In the formula (31), the value of i is 1, 2, 3, N ₁₁ , N ₁₂ , N ₁₃ are respectively the number of R wave, Q wave and S wave in the ECG signal, N ₁₁ ^* , N ₁₂ ^* , N ₁₃ ^* are respectively the number of correctly detected R wave, Q wave and S wave;

Step 4, use the wavelet transform method to obtain the wavelet decomposition diagrams of different levels of the ECG ₂ signal processed in step 2, and then output the position of the Q, R, and S waves by detecting the position and amplitude of the modulus maximum value of the wavelet transform coefficient; Step 4 further includes the following steps:

S4.1 Based on formulas (32) and (33), perform wavelet transformation on the ECG signal using biorthogonal spline wavelet to obtain wavelet decomposition diagrams at different levels;

Among them, x ⁰ (n)=ECG ₂ (n), which is the processed ECG signal; Z is a natural number set (the same below); x ^(j) (n) is the high-frequency signal of x ^j-1 (n) The binary wavelet transform on the j scale, d ^(j) (n) is the binary wavelet transform of the low-frequency signal of x ^j-1 (n) on the j scale, h _j-1 (k), g _{j- 1} (k) is respectively the digital low-pass filter and the high-pass filter coefficient of the k group orthogonality of setting;

S4.2 Based on the wavelet decomposition diagrams at different scales obtained in S4.1, for waveforms at 2 scales, detect the amplitude maximum point M and the minimum value point L; the amplitude between the two points is detected to be p(p< 0.05), p is the position of the R wave of the original ECG signal;

S4.3 Based on the R wave position p obtained in S4.2 and the formula (34) for the waveform at scale 1, search for the modulus minimum value of the front position of the R wave, which is the Q wave.

Q＝x ^(b) (n)<x ^(b) (n-1),x ^(b) (n)<x ^(b) (n+1),(n<p)(34)

Wherein, x ^(b) (n) is the magnitude of the nth sequence of the ECG signal at a scale of 1;

Search for the modulus minimum at the position after the R wave, which is the S wave,

S=x ^(b) (n)<x ^(b) (n-1), x ^(b) (n)<x ^(b) (n+1), (n>p) (35)

Wherein, x ^(b) (n) is the magnitude of the nth sequence of the ECG signal at a scale of 1;

S4.4 Compare the positions of R wave, Q wave and S wave detected by step S4.2 and step S4.3 with the positions marked by experts, and obtain R wave, Q wave and S wave respectively based on formula (36) The detection accuracy of p ₂₁ , p ₂₂ , p ₂₃ :

In the formula (13), the value of i is 1, 2, 3, N ₂₁ , N ₂₂ , N ₂₃ are respectively the number of R wave, Q wave and S wave in the ECG signal, N ₂₁ ^* , N ₂₂ ^* , N ₂₃ ^* are respectively the number of correctly detected R wave, Q wave and S wave;

Wherein, step 4 and step 3 can also exchange order;

Step 5. Repeatedly adjust the artificially set parameters of the differential threshold method through the R wave, Q wave and S wave marked by experts, and repeatedly adjust the wavelet transform hierarchy until the accuracy converges; Step 5 also includes the following steps :

S5.1 Based on the R wave, Q wave and S wave obtained in step 3 and the points marked by experts repeatedly, adjust the values of m ₁ , m ₂ , and m ₃ to decide whether to skip to step 3, specifically:

S5.1A Based on the R wave obtained in step 3, compare it with the R wave marked by experts repeatedly. If the accuracy p ₁₁ improves, adjust the value of m ₁ and skip to step 3; if the accuracy p ₁₁ does not improve, Then skip to step S5.1B;

Step S5.1B is based on the Q wave obtained in step 3, compared with the Q wave marked by experts repeatedly,

If the accuracy _p12 is improved, then adjust the value of _m2 , and skip to step 3; if the accuracy _p12 is not improved, then skip to step S5.1C;

Step S5.1C is based on the S wave obtained in step 3, compared with the S wave marked by experts repeatedly,

If the accuracy _p13 is improved, then adjust the value of _m3 and skip to step 3; if the accuracy _p13 is not improved, then skip to step S5.2;

S5.2 Based on the R wave, Q wave and S wave obtained in step 4 and the points marked by experts repeatedly, adjust the values of b and c to decide whether to skip to step 4, specifically:

S5.2A Based on the R wave obtained in step 4, compare it with the R wave marked by experts repeatedly. If the accuracy p ₂₁ is improved, adjust the value of b and skip to step 4; if the accuracy p ₂₁ does not increase, then Skip to step S5.2B;

Step S5.2B is based on the Q wave and S wave obtained in step 4, compared with the Q wave and S wave marked by experts repeatedly,

If both the accuracy p ₂₂ and p ₂₃ are improved, then adjust the value of c and skip to step 4; if at least one of the accuracy p ₂₂ and p ₂₃ is not improved, then skip to step 6.

Step 6, compare the detection accuracy of step 3 and step 4 respectively, respectively for R wave, Q wave and S wave, select the detection method with high accuracy, step 5 comprises the following steps again:

S6.1 Compare the size of p ₁₁ and p ₂₁ , if p ₁₁ ≥ p ₂₁ , skip to step 3 to extract the R wave, and if p ₁₁ < p ₂₁ , skip to step 4 to extract the R wave;

S6.2 Compare the size of p ₁₂ and p ₂₂ , if p ₁₂ ≥ p ₂₂ , skip to step 3 to extract Q wave, if p ₁₂ < p ₂₂ , skip to step 4 to extract Q wave;

S6.3 Compare the size of p ₁₃ and p ₂₃ , if p ₁₃ ≥ p ₂₃ , skip to step 3 to extract S wave, if p ₁₃ < p ₂₃ , skip to step 4 to extract S wave, This method ends.

It should be noted that what is described in this specification is only a preferred embodiment of the present invention, and the above embodiments are only used to illustrate the technical solution of the present invention rather than limit the present invention. All technical solutions obtained by those skilled in the art through logical analysis, reasoning or limited experiments according to the concept of the present invention shall fall within the scope of the present invention.

Claims (2)

1. An ECG feature extraction method fusing Butterworth filtering and wavelet transformation is characterized in that: carrying out denoising processing on the power frequency interference of the original ECG signal; simultaneously, baseline drift is removed by using a curve fitting mode; after relatively accurate ECG signals are obtained, respectively positioning Q waves, R waves and S waves by using a differential threshold value method and a wavelet transformation method, comparing the positioning Q waves, the R waves and the S waves with information labeled by experts, and respectively adjusting parameters of a threshold value, a scale and a translation quantity to obtain an optimal positioning result;

the method comprises the following steps:

step A,Performing spectral analysis on the original ECG signal to obtain the highest frequency and the lowest frequency, and filtering the original ECG signal by using a Butterworth filter to obtain a filtered ECG₁A signal; step a further comprises the following substeps:

SA.1 analyzing the spectrum of the ECG signal by using the discrete Fourier transform method according to formula (1) to obtain the spectrum X (e) of the ECG signal^jω) Further, the maximum frequency f of the ECG signal is obtained₁And the lowest frequency f of the frequency bands of power frequency interference and other noise interference₂；

Wherein, X (e)^jω) Is the frequency spectrum of the ECG signal, ω is the frequency of the ECG signal, and ECG (n) is the nth amplitude of the ECG signal;

in the obtained X (e)^jω) The highest frequency of the middle and low frequency bands, i.e. the highest frequency f of the ECG signal₁(ii) a The lowest frequency of the high frequency band is the lowest frequency f of the frequency band of power frequency interference and other noise interference₂；

SA.2 filtering the ECG signal with a Butterworth filter to obtain a filtered ECG signal₁Specifically, butterworth filtering is performed by equation (2):

wherein N in formula (2) represents an integer set, | H (ω) | is a mode of the Butterworth filter, | H (ω) | does not include²Is the square of the mode; omega_c(2πf₁＜ω_c＜2πf₂) P (p ∈ N) is the order, which is the cut-off frequency of the Butterworth filter;

step B, ECG outputted from step A₁Performing curve fitting, calculating polynomial coefficient and removing baseline drift to obtain signal after removing baseline drift, ECG₂；

Step B again comprises the following substeps:

sb.1 curve fitting is performed based on equation (3):

where y is the polynomial to be fitted, the order of the polynomial q (q e N), a₀,a₁,...,a_pIs the coefficient of a polynomial of order q, x being the sampling time series, xⁱIs the ith power of x;

SB.2, based on the formula (4), calculating polynomial coefficients in the formula (3);

wherein,is a vector composed of polynomial coefficients,for amplitude vectors of ECG signals, C ═ 1 xx²... x^q]^T，C^TIs the transposition of C;

SB.3 obtaining the signal after baseline drift removal based on the formula (5);

ECG₂＝ECG₁-y (5)

in equation (5), ECG₂For removing the post-baseline-shift signal, ECG₁Is a filtered signal; y is the output of equation (3);

step C, ECG output from step B₂The signal uses a differential threshold method to position R waves, then position Q waves and S waves, and respectively output the positions of the R waves, the Q waves and the S waves;

step C in turn comprises the following substeps:

SC.1 sequentially calculates ECG based on equation (6) and equation (7)₂A first order difference S (n) and a second order difference W (n) of the signals;

S(n)＝ECG₂(n)-ECG₂(n-1) (6)

W(n)＝S(n)-S(n-1) (7)

where S (n) is the nth first order difference, ECG₂(n) is ECG₂Nth amplitude of signal, ECG₂(n-1) is ECG₂The nth-1 amplitude of the signal, W (n) is the nth second order difference, and S (n-1) is the nth-1 first order difference;

SC.2 determines the maximum value max of S (n) output in step SC.1 by using an enumeration method₁And maximum value max of W (n)₂；

SC.3 sets a threshold Y based on formula (8)₁：

Y₁＝m₁·max₁(0＜m₁＜1) (8)

In equation (8), the threshold is Y₁，m₁Is an artificially set coefficient, and the range is (0, 1);

SC.4 judges whether there is a continuous k in S (n)₁(k₁E N) points exceeding Y₁And ECG₂(n) > 0, and according to the judgment result, the following operations are carried out:

SC.4A if there are consecutive k in S (n)₁(k₁E N) points exceeding Y₁And ECG₂(n) > 0, the point n is the R wave crest in the range, the position is R, and the jump is carried out to SC.5;

wherein k is₁Is an integer of 2 or more;

SC.4B if there is no consecutive k in S (n)₁(k₁E N) points exceeding Y₁Then k is₁＝k₁-1, jump to step sc.4;

SC.5 sets the threshold value Y of the starting point based on the equations (9) and (10)₂And an endpoint threshold Y₃：

Y₂＝m₂·max₂(0＜m₂＜1) (9)

Y₃＝m₃·max₂(0＜m₃＜1) (10)

In the formulas (9) and (10), the threshold values are Y₂,Y₃，m₂,m₃Is a coefficient set artificially;

SC.6 search and judge forward from each R waveWhether there is a succession of k in the section W (n)₂(k₂E is N) points less than Y₂And according to the judgment result, the following operations are carried out:

SC.6A if W (n) has consecutive k₂(k₂E is N) points less than Y₂The point n is the starting point s of the QRS wave with the range, and the jump is made to sc.5;

wherein k is₂Is an integer of 2 or more;

SC.6B if there is no consecutive k in W (n)₂(k₂E is N) points less than Y₂Then k is₂＝k₂-1, jump to step sc.6;

SC.7 searches backward from each R wave to determine whether there are consecutive k in W (n)₃(k₃E is N) points less than Y₃And according to the judgment result, the following operations are carried out:

wherein k is₃Is an integer of 2 or more;

SC.7A if W (n) has consecutive k₃(k₃E is N) points less than Y₃The point n is the end point e of the QRS wave with the point in the range, and the jump is carried out to SD.8;

SC.7B if there are no consecutive k in W (n)₃(k₃E is N) points less than Y₃Then k is₃＝k₃-1, jump to step sc.7;

sc.8 searches for an amplitude minimum value point Q in the range of the origin S of Q-wave, R-wave, and S-wave and the R-position where the R-wave is located, based on equation (11):

Q＝min(ECG₂(n)),(s＜n＜r) (11)

in formula (11), min (ECG)₂(n)) is a search ECG₂(n) a function of minima;

based on equation (12), an amplitude minimum value point S is searched for in the range of the R wave and the end point e:

S＝min(ECG₂(n)),(r＜n＜e) (12)

SC.9 compares the positions of the R wave, the Q wave and the S wave determined in the steps SC.4 and SC.8 with the positions marked by the experts, and based on the formula (13), the detection accuracy p of the R wave, the Q wave and the S wave is respectively obtained₁₁,p₁₂,p₁₃：

In the formula (13), the value of i is 1, 2, 3, N₁₁,N₁₂,N₁₃The number of R wave, Q wave and S wave in ECG signal, N₁₁ ^*,N₁₂ ^*,N₁₃ ^*Respectively detecting the number of correct R waves, Q waves and S waves;

step D, processing the ECG processed in the step B₂The signal uses wavelet transform method to get wavelet decomposition pictures of different levels, then outputs Q, R, S wave position by detecting position and amplitude of wavelet transform coefficient modulus maximum; step D comprises the following steps:

SD.1 based on formulas (14) and (15), performing wavelet transformation on the ECG signal by using dual-orthogonal spline wavelets to obtain wavelet decomposition views under different levels;

wherein x is⁰(n)＝ECG₂(n), i.e. the processed ECG signal; z is a natural number set (the same below); x is the number of^(j)(n) is x^j-1(n) a dyadic wavelet transform of the high band signal on the j scale, d^(j)(n) is x^j-1(n) a dyadic wavelet transform of the low band signal on the j scale, h_j-1(k),g_j-1(k) A set kth set of orthogonal digital low pass filter and high pass filter coefficients, respectively;

SD.2, based on wavelet decomposition diagrams under different scales obtained by SD.1, detecting a maximum value point M and a minimum value point L of the amplitude for the waveform under the scale of c (c belongs to N); detecting a point with an amplitude value p (p is less than epsilon) between two points, wherein epsilon is a numerical value with the amplitude value approaching zero, and p is the position of an R wave of the original ECG signal;

SD.3 searches the mode minimum value of the R wave front position based on the R wave position p obtained by SD.2 and the waveform of b (b belongs to N) scale of formula (16), namely Q wave,

Q＝x^(b)(n)＜x^(b)(n-1),x^(b)(n)＜x^(b)(n+1),(n＜p) (16)

wherein x is^(b)(n) is the amplitude of the nth sequence of the ECG signal on the b scale;

searching the mode minimum value of the position after the R wave, namely the S wave,

S＝x^(b)(n)＜x^(b)(n-1),x^(b)(n)＜x^(b)(n+1),(n＞p) (17)

wherein x is^(b)(n) is the amplitude of the nth sequence of the ECG signal on the b scale;

SD.4 compares the positions of the R wave, the Q wave and the S wave detected in the steps SD.2 and SD.3 with the positions marked by experts, and respectively obtains the detection accuracy p of the R wave, the Q wave and the S wave based on a formula (18)₂₁,p₂₂,p₂₃：

In the formula (18), the value of i is 1, 2, 3, N₂₁,N₂₂,N₂₃The number of R wave, Q wave and S wave in ECG signal, N₂₁ ^*,N₂₂ ^*,N₂₃ ^*Respectively detecting the number of correct R waves, Q waves and S waves;

wherein, the step D and the step C can also exchange the sequence;

e, repeatedly adjusting artificially set parameters of a difference threshold method through R waves, Q waves and S waves labeled by experts, and repeatedly adjusting a wavelet transform hierarchy until accuracy is converged; step E also comprises the following steps:

SE.1 adjusting m based on R-, Q-and S-waves obtained in step C, respectively, and by repeated comparison with expert-labeled points₁,m₂,m₃Respectively determining whether to jump to the step C, specifically:

SE.1A based on the R wave obtained in step C, and by repeating the process withComparing the R waves marked by the experts, and judging whether the accuracy is p₁₁Increasing, then adjusting m₁Jumping to step C; if accuracy p₁₁If not, jumping to the step SE.1B;

step se.1b, based on the Q-wave obtained in step C, compares with the Q-wave labeled by the expert repeatedly,

if accuracy p₁₂Increasing, then adjusting m₂Jumping to step C; if accuracy p₁₂If not, jumping to the step SE.1C;

step se.1c, based on the S-wave obtained in step C, compares with the S-wave labeled by the expert repeatedly,

if accuracy p₁₃Increasing, then adjusting m₃Jumping to step C; if accuracy p₁₃If not, jumping to the step SE.2;

and SE.2, respectively based on the R wave, the Q wave and the S wave obtained in the step D and the points marked by the experts through repeated comparison, adjusting the values of b and c, and respectively determining whether to jump to the step D, wherein the specific steps are as follows:

SE.2A based on the R wave obtained in step D, compare with the R wave labeled by the expert through repetition, if the accuracy p₂₁If yes, adjusting the value of b, and jumping to the step D; if accuracy p₂₁If not, jumping to step SE.2B;

step se.2b is based on the Q-wave, S-wave obtained in step D, compared with the Q-wave and S-wave labeled by experts by iteration,

if accuracy p₂₂、p₂₃If the values are all increased, adjusting the value of c, and jumping to the step D; if accuracy p₂₂、p₂₃If at least one is not increased, jumping to the step F;

step F, respectively comparing the detection accuracy of the step C and the detection accuracy of the step D, and selecting a detection method with high accuracy aiming at the R wave, the Q wave and the S wave, wherein the step F comprises the following steps:

SF.1 comparison of p₁₁And p₂₁If p is the size of₁₁≥p₂₁Jumping to the step C to extract R wave, if p₁₁＜p₂₁Jumping to the step D, and extracting R waves;

SF.2 comparison of p₁₂And p₂₂If p is the size of₁₂≥p₂₂Jumping to step C, extracting Q wave, if p₁₂＜p₂₂Jumping to the step D, and extracting the Q wave;

SF.3 comparison of p₁₃And p₂₃If p is the size of₁₃≥p₂₃Jumping to step C, extracting S wave, if p₁₃＜p₂₃And D, jumping to the step D to extract the S wave, and ending the method.

2. The method of extracting features of an ECG incorporating butterworth filtering and wavelet transformation as claimed in claim 1, wherein: step D and step C may also be in an alternating order.