CN108836305B - A kind of ECG feature extracting method of fusion Butterworth filtering and wavelet transformation - Google Patents
A kind of ECG feature extracting method of fusion Butterworth filtering and wavelet transformation Download PDFInfo
- Publication number
- CN108836305B CN108836305B CN201810429640.2A CN201810429640A CN108836305B CN 108836305 B CN108836305 B CN 108836305B CN 201810429640 A CN201810429640 A CN 201810429640A CN 108836305 B CN108836305 B CN 108836305B
- Authority
- CN
- China
- Prior art keywords
- wave
- ecg
- formula
- signal
- point
- 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.)
- Active
Links
Classifications
-
- 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/7225—Details of analog processing, e.g. isolation amplifier, gain or sensitivity adjustment, filtering, baseline or drift compensation
-
- 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)
- Surgery (AREA)
- Animal Behavior & Ethology (AREA)
- Veterinary Medicine (AREA)
- Signal Processing (AREA)
- Physics & Mathematics (AREA)
- Public Health (AREA)
- Biophysics (AREA)
- Pathology (AREA)
- Biomedical Technology (AREA)
- Heart & Thoracic Surgery (AREA)
- Medical Informatics (AREA)
- Molecular Biology (AREA)
- General Health & Medical Sciences (AREA)
- Physiology (AREA)
- Artificial Intelligence (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Psychiatry (AREA)
- Power Engineering (AREA)
- Cardiology (AREA)
- Measurement And Recording Of Electrical Phenomena And Electrical Characteristics Of The Living Body (AREA)
Abstract
A kind of ECG feature extracting method of fusion Butterworth filtering and wavelet transformation, belongs to Medical computer technology field.Core concept are as follows: denoising has been carried out to original ECG signal Hz noise；The mode for using curve matching simultaneously, goes baseline drift；After getting more accurate ECG signal, Q wave, R wave and S wave are positioned respectively respectively using the method for difference threshold algorithm and wavelet transformation, with the information comparison of expert's mark, threshold value, scale and translational movement parameter are adjusted separately, optimal positioning result is obtained.Based on the above results, for Q wave, R wave and S wave, chooses differential threshold respectively or small wave converting method is positioned.The present invention can carry out more accurate denoising to original ECG signal and provide basis further, it is possible to realize the accurate positioning of QRS wave for the digital assay of electrocardiogram.
Description
Technical field
It is filtered the present invention relates to a kind of method of electrocardiosignal feature extraction more particularly to a kind of fusion Butterworth and small
The ECG feature extracting method of wave conversion, belongs to Medical computer technology field.
Background technique
For a long time, cardiovascular disease seriously threatens the health and lives of the mankind because of its high incidence and high mortality.
The World Health Organization, office (WHO) estimates have 36,000,000 people to die of cardiovascular disease, diabetes, breathing system every year in the world at present
The system non-communicable diseases such as disease and malignant tumour account for the 2/3 of global dead sum, arrive the year two thousand twenty, number will rise to 4400
Ten thousand.
" Chinese cardiovascular disease reports (2015) " display, Chinese cardiovascular death rate still occupies disease death composition within 2014
First place, be higher than tumour and other diseases.Rural area cardiovascular death rate in 2014 is 295.63/10 ten thousand, and Heart disease is dead
Dying rate is 143.72/10 ten thousand, and mortality of cerebrovascular disease rate is 151.91/10 ten thousand (cerebral hemorrhage 74.51/10 ten thousand, cerebral infarctions 45.30/10
Ten thousand)；City cardiovascular death rate is 261.99/10 ten thousand, and wherein deaths from heart disease rate is 136.21/10 ten thousand, and the cerebrovascular is died of illness
Dying rate is that 125.78/10 ten thousand (cerebral hemorrhage 52.25/10 ten thousand, cerebral infarction 41.99/10 is ten thousand).With Chinese society's aging paces
Accelerate, for these data in rapid growth, cardiovascular disease has become one of the great public health problem in China.So the heart
The prevention of dirty class disease, diagnosing and treating are particularly important, are the huge challenge and research hotspot of current medical personnel.
Electrocardiogram (Electrocardiogram, ECG) checks the electrode by being placed on human body surface different parts, note
The record electric current that also myocyte generates during each also dynamic cyclic contraction and diastole, while record current Unlimited period side
To with strong and weak variation tendency and real-time depiction electrocardio change curve, the case where to judge cardiac electrical activity, it can be determined that the heart
Whether myocyte is dead, whether there is ischemic, is abnormal with the presence or absence of impulse conduction, and records the characteristic heart with clinical meaning
Electrograph changes.In addition to this, the disease electrocardiogram of certain electrolyte disturbances also has feature to sexually revise.As potassemia will appear T
Wave height point, QT interval reduction, blood calcium reduction will appear ST elongated segment etc.；The feature extraction for carrying out ECG signal has clinical research
There is important meaning.
But in the detection process of electrocardiosignal, there are the influences of the factors such as Hz noise, ambient noise, affect measurement
When precision.And at present since the accuracy and reliability of the intellectual analysis of electrocardiogram can not be with the analysis of cardiologist
It mentions in the same breath, so its ratio in practical clinical and little.
Summary of the invention
It is an object of the invention to further increase the acquisition accuracy of existing electrocardiosignal and the essence of each QRS wave section
Really detection proposes the ECG feature extracting method of a kind of fusion Butterworth filtering and wavelet transformation.
The core idea of the invention is as follows: denoising has been carried out to original ECG signal Hz noise；It is quasi- using curve simultaneously
The mode of conjunction, goes baseline drift；After getting more accurate ECG signal, the side of difference threshold algorithm and wavelet transformation is used
Method respectively positions Q wave, R wave and S wave respectively, and the information comparison with expert's mark adjusts separately threshold value, scale and translation
Parameter is measured, optimal positioning result is obtained.
A kind of ECG feature extracting method of fusion Butterworth filtering and wavelet transformation, includes the following steps:
Step A, spectrum analysis is carried out to original ECG signal and obtains highest frequency and low-limit frequency, reuse Butterworth
Filter is filtered original ECG signal, obtains filtered ECG1Signal；Step A specifically includes following sub-step again
It is rapid:
The method that SA.1 uses discrete Fourier transform by formula (1), analyzes the frequency spectrum of ECG signal, obtains
ECG signal frequency spectrum X (ejω), further obtain the highest frequency f of ECG signal1And Hz noise and other noise jamming frequency ranges
Low-limit frequency f2；
Wherein, X (ejω) be ECG signal frequency spectrum, ω be ECG signal frequency, ECG (n) be ECG signal n-th of width
Value；
In obtained X (ejω) in, the highest frequency of low-frequency range, as the highest frequency f of ECG signal1；High band it is minimum
Frequency is the low-limit frequency f of Hz noise He other noise jamming frequency ranges2；
SA.2 is filtered ECG signal through Butterworth filter, obtains filtered ECG signal ECG1, specific logical
It crosses formula (2) and carries out Butterworth filtering:
Wherein, the N in formula (2) represents set of integers, | H (ω) | it is the mould of Butterworth filter, | H (ω) |2For mould
Square；ωc(2πf1< ωc2 π f of <2) be Butterworth filter cutoff frequency, p (p ∈ N) be order；
Step B, to the ECG of step A output1It carries out curve fitting, seek multinomial coefficient and go baseline drift, gone
Signal after baseline drift, ECG2；
Step B includes following sub-step again:
SB.1 is based on formula (3) and carries out curve fitting:
Wherein, y is the multinomial for being fitted, polynomial order q (q ∈ N), a0,a1,...,apFor q rank multinomial
Coefficient, x are sampling time sequence, xiFor the i-th side of x；
SB.2 is based on formula (4), seeks the multinomial coefficient in formula (3)；
Wherein,For multinomial coefficient composition vector,For the amplitude vector of ECG signal, C=
[1 x x2 ... xq]T, CTFor the transposition of C；
SB.3 is based on formula (5), the signal after obtaining baseline drift；
ECG2=ECG1-y (5)
In formula (5), ECG2To remove the signal after baseline drift, ECG1For filtered signal；Y is the defeated of formula public (3)
Out；
Step C, the ECG for exporting step B2Signal first carries out the positioning of R wave using the method for differential threshold, then carries out Q
The positioning of wave and S wave exports the position of R wave, Q wave and S wave respectively；
Step C includes following sub-step again:
SC.1 is based on formula (6) and formula (7), and sequence calculates ECG2The first-order difference S (n) and second differnce W of signal
(n)；
S (n)=ECG2(n)-ECG2(n-1) (6)
W (n)=S (n)-S (n-1) (7)
Wherein, S (n) is n-th of first-order difference, ECG2It (n) is ECG2N-th of amplitude of signal, ECG2It (n-1) is ECG2
(n-1)th amplitude of signal, W (n) are n-th of second differnce, and S (n-1) is (n-1)th first-order difference；
SC.2 uses enumerative technique, determines the maximum value max of the S (n) of step SC.1 output1With the maximum value max of W (n)2；
SC.3 is based on formula (8), given threshold thresholding Y1:
Y1=m1·max1(0 < m1< 1) (8)
In formula (8), threshold value thresholding is Y1, m1It is the coefficient being manually set, in the range of (0,1)；
SC.4 judges whether there is continuous k in S (n)1(k1∈ N) a point is more than Y1And ECG2(n) 0 >, and according to judging result
It proceeds as follows:
If there is continuous k in SC.4A S (n)1(k1∈ N) a point is more than Y1, and ECG2(n) 0 >, point n is this range
R wave crest, present position r skips to SC.5；
Wherein, k1It is greater than the integer equal to 2；
If without continuous k in SC.4B S (n)1(k1∈ N) a point is more than Y1, then k1=k1- 1, skip to step SC.4；
SC.5 is based on formula (9) and (10), sets starting point threshold value thresholding Y2With terminal threshold value thresholding Y3:
Y2=m2·max2(0 < m2< 1) (9)
Y3=m3·max2(0 < m3< 1) (10)
In formula (9) (10), threshold value thresholding is respectively Y2,Y3, m2,m3It is the coefficient being manually set；
SC.6 is searched for before each R wave direction, judges whether there is continuous k in W (n)2(k2∈ N) a point is less than Y2, and according to sentencing
Disconnected result proceeds as follows:
If there is continuous k in SC.6A W (n)2(k2∈ N) a point is less than Y2, point n is the QRS wave that the point is this range
Starting point s, skip to SC.5；
Wherein, k2It is greater than the integer equal to 2；
If without continuous k in SC.6B W (n)2(k2∈ N) a point is less than Y2, then k2=k2- 1, skip to step SC.6；
SC.7 is searched for after each R wave direction, judges whether there is continuous k in W (n)3(k3∈ N) a point is less than Y3, and according to sentencing
Disconnected result proceeds as follows:
Wherein, k3It is greater than the integer equal to 2；
If there is continuous k in SC.7A W (n)3(k3∈ N) a point is less than Y3, point n is the QRS wave that the point is this range
Terminal e, skip to SD.8；
If without continuous k in SC.7B W (n)3(k3∈ N) a point is less than Y3, then k3=k3- 1, skip to step SC.7；
SC.8 is based on formula (11), searches for amplitude in the range of the position r where starting point s and the R wave in Q wave, R wave and S wave
Minimum point Q:
Q=min (ECG2(n)), (s < n < r) (11)
In formula (11), min (ECG2It (n)) is search ECG2(n) function of minimum；
Based on formula (12), amplitude minimum point S is searched in the range of R wave and terminal e:
S=min (ECG2(n)), (r < n < e) (12)
The position of R wave, Q wave and S wave that SC.9 is determined by step SC.4 and step SC.8, with the position of expert's mark into
Row compares, and is based on formula (13), respectively obtains the accuracy in detection p of R wave, Q wave and S wave11,p12,p13:
In formula (13), the value of i is 1,2,3, N11,N12,N13The number of R wave, Q wave and S wave respectively in ECG signal,
N11 *,N12 *,N13 *Respectively detect the number of correct R wave, Q wave and S wave；
Step D, by step B treated ECG2Signal obtains the wavelet decomposition of different estate using the method for wavelet transformation
Figure, then the position by detection wavelet conversion coefficient modulus maximum and amplitude, export Q, R, S wave position；Step D is wrapped again
Containing following steps:
SD.1 is based on formula (14), (15), carries out wavelet transformation to ECG signal using Bi-orthogonal Spline Wavelet Transformation, obtains not
With the wavelet decomposition figure under stratum；
Wherein, x0(n)=ECG2(n), ECG signal that as treated；Z is nature manifold (similarly hereinafter)；x(j)It (n) is xj-1
(n) dyadic wavelet transform of the high frequency band signal on j scale, d(j)It (n) is xj-1(n) low-band signal is on j scale
Dyadic wavelet transform, hj-1(k),gj-1(k) the kth group of respectively setting orthogonal wave digital lowpass filter and high-pass filter system
Number；
Wavelet decomposition figure under the different scale that SD.2 is obtained based on SD.1 detects width to the waveform under c (c ∈ N) scale
It is worth maximum point M and minimum point L；Detection point-to-point transmission amplitude be p (p < ε) point, wherein ε be amplitude level off to zero number
Value, p are the R wave position of former ECG signal；
SD.3 to the waveform under b (b ∈ N) scale, searches for R wavefront position based on the obtained R wave position p of SD.2 and formula (16)
The modulus minimum set, as Q wave,
Q=x(b)(n) < x(b)(n-1),x(b)(n) < x(b)(n+1), (n < p) (16)
Wherein, x(b)It (n) is the amplitude of n-th of sequence of ECG signal under b scale；
The modulus minimum of position, as S wave after search R wave,
S=x(b)(n) < x(b)(n-1),x(b)(n) < x(b)(n+1), (n > p) (17)
Wherein, x(b)It (n) is the amplitude of n-th of sequence of ECG signal under b scale；
The position of R wave, Q wave and S wave that SD.4 is detected by step SD.2 and step SD.3, with the position of expert's mark into
Row compares, and is based on formula (18), respectively obtains the accuracy in detection p of R wave, Q wave and S wave21,p22,p23:
In formula (18), the value of i is 1,2,3, N21,N22,N23The number of R wave, Q wave and S wave respectively in ECG signal,
N21 *,N22 *,N23 *Respectively detect the number of correct R wave, Q wave and S wave；
Wherein, step D and step C can also be with exchange sequences；
Step E, the R wave marked by expert, Q wave and S wave carry out the parameter of difference threshold algorithm being manually set anti-
Polyphony is whole, is adjusted repeatedly to the stratum of wavelet transformation, until accuracy restrains；Step E is comprised the following steps again:
SE.1 is based respectively on the R wave of step C acquisition, Q wave and S wave with by compared with the point of expert's mark, adjusting repeatedly
m1,m2,m3Value, decide whether to skip to step C respectively, specifically:
The R wave that SE.1A is obtained based on step C is compared, if accuracy p with the R wave by marking repeatedly with expert11
It improves, then adjusts m1Value, skip to step C；If accuracy p11It does not improve, then skips to step SE.1B；
The Q wave that step SE.1B is obtained based on step C, with by repeatedly with expert mark Q Bobbi compared with,
If accuracy p12It improves, then adjusts m2Value, skip to step C；If accuracy p12It does not improve, then skips to step
SE.1C；
The S wave that step SE.1C is obtained based on step C, with by repeatedly with expert mark S Bobbi compared with,
If accuracy p13It improves, then adjusts m3Value, skip to step C；If accuracy p13It does not improve, then skips to step
SE.2；
SE.2 is based respectively on the R wave of step D acquisition, Q wave and S wave with by compared with the point of expert's mark, adjusting repeatedly
The value of b, c decides whether to skip to step D respectively, specifically:
The R wave that SE.2A is obtained based on step D is compared, if accuracy p with the R wave by marking repeatedly with expert21
It improves, then adjusts the value of b, skip to step D；If accuracy p21It does not improve, then skips to step SE.2B；
Q wave that step SE.2B is obtained based on step D, S wave, with by repeatedly with the Q wave of expert's mark and S Bobbi compared with,
If accuracy p22、p23It improves, then adjusts the value of c, skip to step D；If accuracy p22、p23At least one does not have
It improves, then skips to step F.
Step F, it is respectively compared the accuracy in detection of step C and step D, is directed to R wave, Q wave and S wave respectively, selection is accurate
High detection method is spent, step F is comprised the following steps again:
SF.1 compares p11And p21Size, if p11≥p21, then step C is skipped to, the extraction of R wave is carried out, if p11< p21, then
Step D is skipped to, the extraction of R wave is carried out；
SF.2 compares p12And p22Size, if p12≥p22, then step C is skipped to, the extraction of Q wave is carried out, if p12< p22, then
Step D is skipped to, the extraction of Q wave is carried out；
SF.3 compares p13And p23Size, if p13≥p23, then step C is skipped to, the extraction of S wave is carried out, if p13< p23, then
Step D is skipped to, the extraction of S wave is carried out, this method terminates.
Beneficial effect
A kind of ECG feature extracting method of fusion Butterworth filtering and wavelet transformation of the present invention, compares prior art, tool
It has the advantages that: more accurate denoising can be carried out to original ECG signal, further, it is possible to realize QRS wave
It is accurately positioned, provides basis for the digital assay of electrocardiogram.
Detailed description of the invention
Fig. 1 be the Q wave of ECG feature extracting method of a kind of fusion Butterworth filtering of the present invention and wavelet transformation, R wave with
And the extraction flow chart of S wave.
Specific embodiment
In the following, developing simultaneously embodiment in conjunction with attached drawing, the present invention will be described in detail.
Embodiment 1
The invention proposes the ECG signal feature extracting method of a kind of fusion Butterworth filtering and wavelet transformation, such as Fig. 1
It is shown.Specifically comprise the following steps:
Step 1 obtains highest frequency and low-limit frequency to original ECG signal progress spectrum analysis, reuses Butterworth
Filter is filtered original ECG signal, obtains filtered ECG1Signal；Step 1 specifically includes following sub-step again
It is rapid:
The method that S1.1 uses discrete Fourier transform by formula (19), analyzes the frequency spectrum of ECG signal, obtains
The highest frequency and low-limit frequency of ECG signal frequency spectrum, X (ejω):
Wherein, X (ejω) be ECG signal frequency spectrum, ω be ECG signal frequency, ECG (n) be ECG signal n-th of width
Value；
In obtained X (ejω) in, the highest frequency of low-frequency range, as the highest frequency f of ECG signal1；High band it is minimum
Frequency is the low-limit frequency f of Hz noise He other noise jamming frequency ranges2；
S1.2 is filtered ECG signal through Butterworth filter, obtains filtered ECG signal ECG1, based on public affairs
Formula (20) carries out Butterworth filtering；
Wherein, the N in formula (20) represents set of integers, | H (ω) | it is the mould of Butterworth filter, | H (ω) |2For mould
Square；ωc=2 π × 50Hz is the cutoff frequency of Butterworth filter, and p=2 is order；
Step 2, the ECG that step SA.2 is exported1It carries out curve fitting, seek multinomial coefficient and go baseline drift, obtain
To removing the signal ECG after baseline drift2；
Step 2 includes following sub-step again:
S2.1 is based on formula (21) and carries out curve fitting；
Wherein, y is the multinomial for being fitted, polynomial order q=3, a0,a1,a2,a3What it is for q rank multinomial is
Number, x are sampling time sequence, xiFor the i-th side of x；
S2.2 is based on formula (22), seeks the multinomial coefficient in formula (3)；
Wherein,For multinomial coefficient composition vector,For the amplitude vector of ECG signal, C=
[1 x x2 x3]T, CTFor the transposition of C；
S2.3 is based on formula (23), the signal after obtaining baseline drift；
ECG2=ECG1-y (23)
In formula (23), ECG2To remove the signal after baseline drift, ECG1For filtered signal；Y is the defeated of formula public (21)
Out；
Step 3, the ECG for exporting step 22Signal first carries out the positioning of R wave using the method for differential threshold, then carries out Q
The positioning of wave and S wave exports the position of R wave, Q wave and S wave respectively；
Step 3 includes following sub-step again:
S3.1 is based on formula (24) and formula (25), and sequence calculates ECG2The first-order difference S (n) and second differnce W of signal
(n)；
S (n)=ECG2(n)-ECG2(n-1) (24)
W (n)=S (n)-S (n-1) (25)
Wherein, S (n) is n-th of first-order difference, ECG2It (n) is ECG2N-th of amplitude of signal, ECG2It (n-1) is ECG2
(n-1)th amplitude of signal, W (n) are n-th of second differnce, and S (n-1) is (n-1)th first-order difference；
S3.2 uses enumerative technique, determines the maximum value max of the S (n) of step S3.1 output1With the maximum value max of W (n)2；
S3.3 is based on formula (26), given threshold thresholding Y1:
Y1=m1·max1(0 < m1< 1) (26)
In formula (26), threshold value thresholding is Y1, m1=0.7；
S3.4 judges whether there is continuous k in S (n)1(k1∈ N) a point is more than Y1And ECG2(n) 0 >, and according to judging result
It proceeds as follows:
If there is continuous k in S3.4A S (n)1(k1∈ N) a point is more than Y1, and ECG2(n) 0 >, point n is this range
R wave crest, present position r skips to S3.5；
Wherein, k1It is greater than the integer equal to 2；
If without continuous k in S3.4B S (n)1(k1∈ N) a point is more than Y1, then k1=k1- 1, skip to step SC.4；
S3.5 is based on formula (27) and (28), sets starting point threshold value thresholding Y2With terminal threshold value thresholding Y3:
Y2=m2·max2(0 < m2< 1) (27)
Y3=m3·max2(0 < m3< 1) (28)
In formula (27) (28), threshold value thresholding is respectively Y2,Y3, m2=0.5, m3=0.6；
S3.6 is searched for before each R wave direction, judges whether there is continuous k in W (n)2(k2∈ N) a point is less than Y2, and according to sentencing
Disconnected result proceeds as follows:
If there is continuous k in S3.6A W (n)2=5 points are less than Y2, point n is rising for the QRS wave that the point is this range
Point s, skips to S3.5；
Wherein, k2It is greater than the integer equal to 2；
If without continuous k in S3.6B W (n)2=5 points are less than Y2, then k2=k2- 1, skip to step S3.6；
S3.7 is searched for after each R wave direction, judges whether there is continuous k in W (n)3=5 points are less than Y3, and tied according to judgement
Fruit proceeds as follows:
Wherein, k3It is greater than the integer equal to 2；
If there is continuous k in S3.7A W (n)3=5 points are less than Y3, point n is the end for the QRS wave that the point is this range
Point e, skips to S4.8；
If without continuous k in S3.7B W (n)3=5 points are less than Y3, then k3=k3- 1, skip to step S3.7；
S3.8 is based on formula (29), searches for amplitude in the range of the position r where starting point s and the R wave in Q wave, R wave and S wave
Minimum point Q:
Q=min (ECG2(n)), (s < n < r) (29)
In formula (29), min (ECG2It (n)) is search ECG2(n) function of minimum；
Based on formula (30), amplitude minimum point S is searched in the range of R wave and terminal e:
S=min (ECG2(n)), (r < n < e) (30)
The position of R wave, Q wave and S wave that S3.9 is determined by step S3.4 and step S3.8, with the position of expert's mark into
Row compares, and is based on formula (13), respectively obtains the accuracy in detection p of R wave, Q wave and S wave11,p12,p13:
In formula (31), the value of i is 1,2,3, N11,N12,N13The number of R wave, Q wave and S wave respectively in ECG signal,
N11 *,N12 *,N13 *Respectively detect the number of correct R wave, Q wave and S wave；
Step 4, by step 2 treated ECG2Signal obtains the wavelet decomposition of different estate using the method for wavelet transformation
Figure, then the position by detection wavelet conversion coefficient modulus maximum and amplitude, export Q, R, S wave position；Step 4 is wrapped again
Containing following steps:
S4.1 is based on formula (32), (33), carries out wavelet transformation to ECG signal using Bi-orthogonal Spline Wavelet Transformation, obtains not
With the wavelet decomposition figure under stratum；
Wherein, x0(n)=ECG2(n), ECG signal that as treated；Z is nature manifold (similarly hereinafter)；x(j)It (n) is xj-1
(n) dyadic wavelet transform of the high frequency band signal on j scale, d(j)It (n) is xj-1(n) low-band signal is on j scale
Dyadic wavelet transform, hj-1(k),gj-1(k) the kth group of respectively setting orthogonal wave digital lowpass filter and high-pass filter system
Number；
Wavelet decomposition figure under the different scale that S4.2 is obtained based on S4.1, to the waveform under 2 scales, detected amplitude is very big
It is worth point M and minimum point L；Detection point-to-point transmission amplitude is the point of p (p < 0.05), and p is the R wave position of former ECG signal；
S4.3 to the waveform under 1 scale, searches for the mould of R wavefront-position based on the obtained R wave position p of S4.2 and formula (34)
Minimum, as Q wave,
Q=x(b)(n) < x(b)(n-1),x(b)(n) < x(b)(n+1), (n < p) (34)
Wherein, x(b)It (n) is the amplitude of n-th of sequence of ECG signal under 1 scale；
The modulus minimum of position, as S wave after search R wave,
S=x(b)(n) < x(b)(n-1),x(b)(n) < x(b)(n+1), (n > p) (35)
Wherein, x(b)It (n) is the amplitude of n-th of sequence of ECG signal under 1 scale；
The position of R wave, Q wave and S wave that S4.4 is detected by step S4.2 and step S4.3, with the position of expert's mark into
Row compares, and is based on formula (36), respectively obtains the accuracy in detection p of R wave, Q wave and S wave21,p22,p23:
In formula (13), the value of i is 1,2,3, N21,N22,N23The number of R wave, Q wave and S wave respectively in ECG signal,
N21 *,N22 *,N23 *Respectively detect the number of correct R wave, Q wave and S wave；
Wherein, step 4 and step 3 can also be with exchange sequences；
Step 5, the R wave marked by expert, Q wave and S wave carry out the parameter of difference threshold algorithm being manually set anti-
Polyphony is whole, is adjusted repeatedly to the stratum of wavelet transformation, until accuracy restrains；Step 5 comprises the following steps again:
S5.1 is based respectively on the R wave of step 3 acquisition, Q wave and S wave with by compared with the point of expert's mark, adjusting repeatedly
m1,m2,m3Value, decide whether to skip to step 3 respectively, specifically:
The R wave that S5.1A is obtained based on step 3 is compared, if accuracy p with the R wave by marking repeatedly with expert11
It improves, then adjusts m1Value, skip to step 3；If accuracy p11It does not improve, then skips to step S5.1B；
The Q wave that step S5.1B is obtained based on step 3, with by repeatedly with expert mark Q Bobbi compared with,
If accuracy p12It improves, then adjusts m2Value, skip to step 3；If accuracy p12It does not improve, then skips to step
S5.1C；
The S wave that step S5.1C is obtained based on step 3, with by repeatedly with expert mark S Bobbi compared with,
If accuracy p13It improves, then adjusts m3Value, skip to step 3；If accuracy p13It does not improve, then skips to step
S5.2；
S5.2 is based respectively on the R wave of step 4 acquisition, Q wave and S wave with by compared with the point of expert's mark, adjusting repeatedly
The value of b, c decides whether to skip to step 4 respectively, specifically:
The R wave that S5.2A is obtained based on step 4 is compared, if accuracy p with the R wave by marking repeatedly with expert21
It improves, then adjusts the value of b, skip to step 4；If accuracy p21It does not improve, then skips to step S5.2B；
Q wave that step S5.2B is obtained based on step 4, S wave, with by repeatedly with the Q wave of expert's mark and S Bobbi compared with,
If accuracy p22、p23It improves, then adjusts the value of c, skip to step 4；If accuracy p22、p23At least one does not have
It improves, then skips to step 6.
Step 6, the accuracy in detection for being respectively compared step 3 and step 4, are directed to R wave, Q wave and S wave respectively, and selection is accurate
High detection method is spent, step 5 comprises the following steps again:
S6.1 compares p11And p21Size, if p11≥p21, then step 3 is skipped to, the extraction of R wave is carried out, if p11< p21, then
Step 4 is skipped to, the extraction of R wave is carried out；
S6.2 compares p12And p22Size, if p12≥p22, then step 3 is skipped to, the extraction of Q wave is carried out, if p12< p22, then
Step 4 is skipped to, the extraction of Q wave is carried out；
S6.3 compares p13And p23Size, if p13≥p23, then step 3 is skipped to, the extraction of S wave is carried out, if p13< p23, then
Step 4 is skipped to, the extraction of S wave is carried out, this method terminates.
It should be noted that being only presently preferred embodiments of the present invention described in this specification, above embodiments are only used
In illustrating technical solution of the present invention rather than limitation of the present invention.All those skilled in the art pass through under this invention's idea to patrol
Analysis, reasoning or the limited available technical solution of experiment are collected, it all should be within the scope of the present invention.
Claims (1)
1. a kind of ECG feature extracting method of fusion Butterworth filtering and wavelet transformation, it is characterised in that: believe raw ECG
Number Hz noise has carried out denoising；The mode for using curve matching simultaneously, goes baseline drift；It is more accurate getting
After ECG signal, adjust threshold value and small wave converting method adjustment scale and translational movement parameter using difference threshold algorithm, respectively to Q wave,
R wave and S wave are positioned respectively, and the information comparison with expert's mark obtains optimal positioning result；Include the following steps:
Step A, spectrum analysis is carried out to original ECG signal and obtains highest frequency and low-limit frequency, reuse Butterworth filtering
Device is filtered original ECG signal, obtains filtered ECG1Signal；Step A specifically includes following sub-step again:
The method that SA.1 uses discrete Fourier transform by formula (1), analyzes the frequency spectrum of ECG signal, obtains ECG letter
Number frequency spectrum X (ejω), further obtain the highest frequency f of ECG signal1And the lowest frequency of Hz noise and other noise jamming frequency ranges
Rate f2；
Wherein, X (ejω) be ECG signal frequency spectrum, ω be ECG signal frequency, ECG (n) be ECG signal n-th of amplitude；
In obtained X (ejω) in, the highest frequency of low-frequency range, as the highest frequency f of ECG signal1；The lowest frequency of high band
Rate is the low-limit frequency f of Hz noise He other noise jamming frequency ranges2；
SA.2 is filtered ECG signal through Butterworth filter, obtains filtered ECG signal ECG1, especially by formula
(2) Butterworth filtering is carried out:
Wherein, p (p ∈ N) is order in formula (2), and N represents set of integers；| H (ω) | it is the mould of Butterworth filter, | H
(ω)|2For square of mould；ωc(2πf1< ωc2 π f of <2) be Butterworth filter cutoff frequency；
Step B, to the ECG of step A output1It carries out curve fitting, seek multinomial coefficient and go baseline drift, obtain baseline
Signal after drift, ECG2；
Step B includes following sub-step again:
SB.1 is based on formula (3) and carries out curve fitting:
Wherein, y is the multinomial for being fitted, polynomial order q (q ∈ N), a0,a1,...,apWhat it is for q rank multinomial is
Number, x are sampling time sequence, xiFor the i-th side of x；
SB.2 is based on formula (4), seeks the multinomial coefficient in formula (3)；
Wherein,For multinomial coefficient composition vector,For the amplitude vector of ECG signal, C=[1 x
x2 ... xq]T, CTFor the transposition of C；
SB.3 is based on formula (5), the signal after obtaining baseline drift；
ECG2=ECG1-y (5)
In formula (5), ECG2To remove the signal after baseline drift, ECG1For filtered signal；Y is the output of formula public (3)；
Step C, the ECG for exporting step B2Signal first carries out the positioning of R wave using the method for differential threshold, then carries out Q wave and S
The positioning of wave exports the position of R wave, Q wave and S wave respectively；
Step C includes following sub-step again:
SC.1 is based on formula (6) and formula (7), and sequence calculates ECG2The first-order difference S (n) and second differnce W (n) of signal；
S (n)=ECG2(n)-ECG2(n-1) (6)
W (n)=S (n)-S (n-1) (7)
Wherein, S (n) is n-th of first-order difference, ECG2It (n) is ECG2N-th of amplitude of signal, ECG2It (n-1) is ECG2Signal
(n-1)th amplitude, W (n) are n-th of second differnce, and S (n-1) is (n-1)th first-order difference；
SC.2 uses enumerative technique, determines the maximum value max of the S (n) of step SC.1 output1With the maximum value max of W (n)2；
SC.3 is based on formula (8), given threshold thresholding Y1:
Y1=m1·max1(0 < m1< 1) (8)
In formula (8), threshold value thresholding is Y1, m1It is the coefficient being manually set, in the range of (0,1)；
SC.4 judges whether there is continuous k in S (n)1(k1∈ N) a point is more than Y1And ECG2(n) 0 >, and carried out according to judging result
Following operation:
If there is continuous k in SC.4A S (n)1(k1∈ N) a point is more than Y1, and ECG2(n) 0 >, point n are the R wave of this range
Peak, present position r, skips to SC.5；
Wherein, k1It is greater than the integer equal to 2；
If without continuous k in SC.4B S (n)1(k1∈ N) a point is more than Y1, then k1=k1- 1, skip to step SC.4；
SC.5 is based on formula (9) and (10), sets starting point threshold value thresholding Y2With terminal threshold value thresholding Y3:
Y2=m2·max2(0 < m2< 1) (9)
Y3=m3·max2(0 < m3< 1) (10)
In formula (9) and (10), Y2、Y3Respectively threshold value thresholding；m2, m3It is the coefficient being manually set；
SC.6 is searched for before each R wave direction, judges whether there is continuous k in W (n)2(k2∈ N) a point is less than Y2, and tied according to judgement
Fruit proceeds as follows:
If there is continuous k in SC.6A W (n)2(k2∈ N) a point is less than Y2, which is the starting point s of the QRS wave of this range, is skipped to
SC.5；
Wherein, k2It is greater than the integer equal to 2；
If without continuous k in SC.6B W (n)2(k2∈ N) a point is less than Y2, then k2=k2- 1, skip to step SC.6；
SC.7 is searched for after each R wave direction, judges whether there is continuous k in W (n)3(k3∈ N) a point is less than Y3, and tied according to judgement
Fruit proceeds as follows:
Wherein, k3It is greater than the integer equal to 2；
If there is continuous k in SC.7A W (n)3(k3∈ N) a point is less than Y3, which is the terminal e of the QRS wave of this range, is skipped to
SD.8；
If without continuous k in SC.7B W (n)3(k3∈ N) a point is less than Y3, then k3=k3- 1, skip to step SC.7；
SC.8 is based on formula (11), and search amplitude is minimum in the range of the position r where starting point s and the R wave in Q wave, R wave and S wave
It is worth point Q:
Q=min (ECG2(n)), (s < n < r) (11)
In formula (11), min (ECG2It (n)) is search ECG2(n) function of minimum；
Based on formula (12), amplitude minimum point S is searched in the range of R wave and terminal e:
S=min (ECG2(n)), (r < n < e) (12)
The position of R wave, Q wave and S wave that SC.9 is determined by step SC.4 and step SC.8 is compared with the position of expert's mark
Compared with, be based on formula (13), respectively obtain the accuracy in detection p of R wave, Q wave and S wave11,p12,p13:
In formula (13), the value of i is 1,2,3, N11,N12,N13The number of R wave, Q wave and S wave, N respectively in ECG signal11 *,
N12 *,N13 *Respectively detect the number of correct R wave, Q wave and S wave；
Step D, by step B treated ECG2Signal obtains the wavelet decomposition figure of different estate using the method for wavelet transformation, so
Position by detection wavelet conversion coefficient modulus maximum and amplitude afterwards, export Q, R, S wave position；Step D includes such as again
Lower step:
SD.1 is based on formula (14), (15), carries out wavelet transformation to ECG signal using Bi-orthogonal Spline Wavelet Transformation, obtains not same order
Wavelet decomposition figure under layer；
Wherein, x0(n)=ECG2(n), ECG signal that as treated；Z is nature manifold (similarly hereinafter)；x(j)It (n) is xj-1(n)
Dyadic wavelet transform of the high frequency band signal on j scale, d(j)It (n) is xj-1(n) low-band signal on j scale two into small
Wave conversion, hj-1(k),gj-1(k) the kth group of respectively setting orthogonal wave digital lowpass filter and high-pass filter coefficient；
Wavelet decomposition figure under the different scale that SD.2 is obtained based on SD.1, to the waveform under c (c ∈ N) scale, detected amplitude pole
Big value point M and minimum point L；Detection point-to-point transmission amplitude is the point of p (p < ε), wherein ε be amplitude level off to zero numerical value, p is
For the R wave position of former ECG signal；
SD.3 to the waveform under b (b ∈ N) scale, searches for R wavefront-position based on the obtained R wave position p of SD.2 and formula (16)
Modulus minimum, as Q wave,
Q=x(b)(n) < x(b)(n-1),x(b)(n) < x(b)(n+1), (n < p) (16)
Wherein, x(b)It (n) is the amplitude of n-th of sequence of ECG signal under b scale；
The modulus minimum of position, as S wave after search R wave,
S=x(b)(n) < x(b)(n-1),x(b)(n) < x(b)(n+1), (n > p) (17)
Wherein, x(b)It (n) is the amplitude of n-th of sequence of ECG signal under b scale；
The position of R wave, Q wave and S wave that SD.4 is detected by step SD.2 and step SD.3 is compared with the position of expert's mark
Compared with, be based on formula (18), respectively obtain the accuracy in detection p of R wave, Q wave and S wave21,p22,p23:
In formula (18), the value of i is 1,2,3, N21,N22,N23The number of R wave, Q wave and S wave, N respectively in ECG signal21 *,
N22 *,N23 *Respectively detect the number of correct R wave, Q wave and S wave；
Wherein, step D and step C can also be with exchange sequences；
Step E, the R wave marked by expert, Q wave and S wave adjust the parameter of difference threshold algorithm being manually set repeatedly
It is whole, the stratum of wavelet transformation is adjusted repeatedly, until accuracy restrains；Step E is comprised the following steps again:
SE.1 is based respectively on the R wave of step C acquisition, Q wave and S wave with by repeatedly compared with the point of expert's mark, adjustment m1,m2,
m3Value, decide whether to skip to step C respectively, specifically:
The R wave that SE.1A is obtained based on step C is compared, if accuracy p with the R wave by marking repeatedly with expert11It improves,
Then adjust m1Value, skip to step C；If accuracy p11It does not improve, then skips to step SE.1B；
The Q wave that step SE.1B is obtained based on step C, with by repeatedly with expert mark Q Bobbi compared with,
If accuracy p12It improves, then adjusts m2Value, skip to step C；If accuracy p12It does not improve, then skips to step SE.1C；
The S wave that step SE.1C is obtained based on step C, with by repeatedly with expert mark S Bobbi compared with,
If accuracy p13It improves, then adjusts m3Value, skip to step C；If accuracy p13It does not improve, then skips to step SE.2；
SE.2 is based respectively on the R wave of step D acquisition, Q wave and S wave with by compared with the point of expert's mark, adjusting b, c repeatedly
Value, decides whether to skip to step D, specifically:
The R wave that SE.2A is obtained based on step D is compared, if accuracy p with the R wave by marking repeatedly with expert21It improves,
The value for then adjusting b skips to step D；If accuracy p21It does not improve, then skips to step SE.2B；
Q wave that step SE.2B is obtained based on step D, S wave, with by repeatedly with the Q wave of expert's mark and S Bobbi compared with,
If accuracy p22、p23It improves, then adjusts the value of c, skip to step D；If accuracy p22、p23At least one is not improved,
Then skip to step F；
Step F, it is respectively compared the accuracy in detection of step C and step D, is directed to R wave, Q wave and S wave respectively, accuracy of selection is high
Detection method, step F comprises the following steps again:
SF.1 compares p11And p21Size, if p11≥p21, then step C is skipped to, the extraction of R wave is carried out, if p11< p21, then skip to
Step D carries out the extraction of R wave；
SF.2 compares p12And p22Size, if p12≥p22, then step C is skipped to, the extraction of Q wave is carried out, if p12< p22, then skip to
Step D carries out the extraction of Q wave；
SF.3 compares p13And p23Size, if p13≥p23, then step C is skipped to, the extraction of S wave is carried out, if p13< p23, then skip to
Step D, carries out the extraction of S wave, and this method terminates.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201810429640.2A CN108836305B (en) | 2018-05-08 | 2018-05-08 | A kind of ECG feature extracting method of fusion Butterworth filtering and wavelet transformation |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201810429640.2A CN108836305B (en) | 2018-05-08 | 2018-05-08 | A kind of ECG feature extracting method of fusion Butterworth filtering and wavelet transformation |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN108836305A CN108836305A (en) | 2018-11-20 |
| CN108836305B true CN108836305B (en) | 2019-04-12 |
Family
ID=64212855
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201810429640.2A Active CN108836305B (en) | 2018-05-08 | 2018-05-08 | A kind of ECG feature extracting method of fusion Butterworth filtering and wavelet transformation |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN108836305B (en) |
Families Citing this family (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN109542021A (en) * | 2018-12-24 | 2019-03-29 | 广东理致技术有限公司 | A kind of sensor weak signal data acquisition method and device |
| CN112244862B (en) * | 2020-10-14 | 2024-05-14 | 哈尔滨理工大学 | Electrocardiogram signal denoising algorithm based on RFDA wavelet threshold |
| CN112611444B (en) * | 2020-12-30 | 2023-04-28 | 西安和其光电科技股份有限公司 | Distributed optical fiber vibration monitoring system and method capable of achieving accurate positioning |
Citations (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US5025794A (en) * | 1988-08-30 | 1991-06-25 | Corazonix Corporation | Method for analysis of electrocardiographic signal QRS complex |
| CN101799974A (en) * | 2010-03-12 | 2010-08-11 | 上海交通大学 | Electrocardio signal transmission method based on self-adaptive codebook |
| CN101828916A (en) * | 2010-05-07 | 2010-09-15 | 深圳大学 | Electrocardiosignal processing system |
| CN106037720A (en) * | 2015-12-04 | 2016-10-26 | 贵州大学 | Application method of hybrid continuous information analysis technology in medicine |
| CN107693011A (en) * | 2017-11-13 | 2018-02-16 | 湖北科技学院 | A kind of ECG signal baseline filtering method |
| CN107788969A (en) * | 2017-09-29 | 2018-03-13 | 成都瑞迪康医疗科技有限公司 | The automatic testing method of QRS complex in a kind of electrocardiosignal |
Family Cites Families (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US9408549B2 (en) * | 2009-11-03 | 2016-08-09 | Vivaquant Llc | Detecting fiducial points in physiological signals |
-
2018
- 2018-05-08 CN CN201810429640.2A patent/CN108836305B/en active Active
Patent Citations (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US5025794A (en) * | 1988-08-30 | 1991-06-25 | Corazonix Corporation | Method for analysis of electrocardiographic signal QRS complex |
| CN101799974A (en) * | 2010-03-12 | 2010-08-11 | 上海交通大学 | Electrocardio signal transmission method based on self-adaptive codebook |
| CN101828916A (en) * | 2010-05-07 | 2010-09-15 | 深圳大学 | Electrocardiosignal processing system |
| CN106037720A (en) * | 2015-12-04 | 2016-10-26 | 贵州大学 | Application method of hybrid continuous information analysis technology in medicine |
| CN107788969A (en) * | 2017-09-29 | 2018-03-13 | 成都瑞迪康医疗科技有限公司 | The automatic testing method of QRS complex in a kind of electrocardiosignal |
| CN107693011A (en) * | 2017-11-13 | 2018-02-16 | 湖北科技学院 | A kind of ECG signal baseline filtering method |
Also Published As
| Publication number | Publication date |
|---|---|
| CN108836305A (en) | 2018-11-20 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN110840402B (en) | Atrial fibrillation signal identification method and system based on machine learning | |
| Haddadi et al. | Discrete wavelet transform based algorithm for recognition of QRS complexes | |
| CN109907752B (en) | Electrocardiogram diagnosis and monitoring system for removing motion artifact interference and electrocardio characteristic detection | |
| CN110327036B (en) | Method for extracting respiratory signal and respiratory frequency from wearable electrocardiogram | |
| CN108836305B (en) | A kind of ECG feature extracting method of fusion Butterworth filtering and wavelet transformation | |
| CN103405227B (en) | Double-layer morphological filter based electrocardiosignal preprocessing method | |
| CN106108889A (en) | Electrocardiogram classification method based on degree of depth learning algorithm | |
| CN204931634U (en) | Based on the depression evaluating system of physiologic information | |
| CN106889981B (en) | A kind of intelligent terminal for being used to extract fetal heart frequency | |
| CN102626310A (en) | Electrocardiogram signal feature detection algorithm based on wavelet transformation lifting and approximate envelope improving | |
| CN114532993B (en) | Automatic detection method for electroencephalogram high-frequency oscillation signals of epileptic patients | |
| CN103761424A (en) | Electromyography signal noise reducing and aliasing removing method based on second-generation wavelets and ICA (independent component analysis) | |
| CN106485208A (en) | The automatic removal method of eye electrical interference in single channel EEG signals | |
| CN108420406A (en) | Method based on pulse wave sleep stage | |
| Li et al. | A low-complexity ECG processing algorithm based on the Haar wavelet transform for portable health-care devices | |
| CN110490257A (en) | It is a kind of based on the electro-physiological signals entropy analysis method for removing trend term | |
| CN102783945A (en) | Fetal electrocardiogram signal extracting method based on wavelet threshold denoising | |
| Liu et al. | A novel P-QRS-T wave localization method in ECG signals based on hybrid neural networks | |
| CN111887811B (en) | Brain abnormal discharge detection method and system based on electroencephalogram signal characteristics | |
| CN103876731B (en) | A kind of Fetal ECG signal extracting device and method | |
| CN110755069B (en) | Dynamic electrocardiosignal baseline drift correction method for jump mutation noise | |
| CN110179456B (en) | Electrocardio noise recognition model training and electrocardio noise detection method and device | |
| CN111481193A (en) | Fall risk assessment and early warning method and system | |
| CN108836316B (en) | Electrocardiosignal R wave extraction method based on BP neural network | |
| CN116269422A (en) | Fetal electrocardio separation acquisition method and device with high signal-to-noise ratio |
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 | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |