CN110680308A - Electrocardiosignal denoising method based on fusion of improved EMD and threshold method - Google Patents

Electrocardiosignal denoising method based on fusion of improved EMD and threshold method Download PDF

Info

Publication number
CN110680308A
CN110680308A CN201911070409.XA CN201911070409A CN110680308A CN 110680308 A CN110680308 A CN 110680308A CN 201911070409 A CN201911070409 A CN 201911070409A CN 110680308 A CN110680308 A CN 110680308A
Authority
CN
China
Prior art keywords
imf
signal
point
denoising
emd
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
CN201911070409.XA
Other languages
Chinese (zh)
Other versions
CN110680308B (en
Inventor
郭树理
何昆仑
尹俭芳
韩丽娜
范利
曹丰
刘宏斌
王春喜
李玉龙
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Beijing University of Technology
Chinese PLA General Hospital
Original Assignee
Beijing University of Technology
Chinese PLA General Hospital
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by Beijing University of Technology, Chinese PLA General Hospital filed Critical Beijing University of Technology
Priority to CN201911070409.XA priority Critical patent/CN110680308B/en
Publication of CN110680308A publication Critical patent/CN110680308A/en
Application granted granted Critical
Publication of CN110680308B publication Critical patent/CN110680308B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

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
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B5/00—Measuring for diagnostic purposes; Identification of persons
    • A61B5/24—Detecting, measuring or recording bioelectric or biomagnetic signals of the body or parts thereof
    • A61B5/316—Modalities, i.e. specific diagnostic methods
    • A61B5/318—Heart-related electrical modalities, e.g. electrocardiography [ECG]
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B5/00—Measuring for diagnostic purposes; Identification of persons
    • A61B5/72—Signal processing specially adapted for physiological signals or for diagnostic purposes
    • A61B5/7203—Signal processing specially adapted for physiological signals or for diagnostic purposes for noise prevention, reduction or removal
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B5/00—Measuring for diagnostic purposes; Identification of persons
    • A61B5/72—Signal processing specially adapted for physiological signals or for diagnostic purposes
    • A61B5/7235—Details of waveform analysis

Landscapes

  • Health & Medical Sciences (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Engineering & Computer Science (AREA)
  • Surgery (AREA)
  • Animal Behavior & Ethology (AREA)
  • Veterinary Medicine (AREA)
  • Public Health (AREA)
  • Physics & Mathematics (AREA)
  • General Health & Medical Sciences (AREA)
  • Biophysics (AREA)
  • Pathology (AREA)
  • Biomedical Technology (AREA)
  • Heart & Thoracic Surgery (AREA)
  • Medical Informatics (AREA)
  • Molecular Biology (AREA)
  • Signal Processing (AREA)
  • Artificial Intelligence (AREA)
  • Physiology (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Psychiatry (AREA)
  • Cardiology (AREA)
  • Measurement And Recording Of Electrical Phenomena And Electrical Characteristics Of The Living Body (AREA)

Abstract

The invention provides an electrocardiosignal denoising method based on the fusion of an improved EMD (empirical mode decomposition) method and a threshold method, and belongs to the technical field of signal filtering. The method solves the problem of modal aliasing by superposing white noise with different weight coefficients, solves the problem of end points by a method of a least square support vector machine, constructs upper and lower envelope lines of a signal by a method of conformal spline interpolation, constructs cubic spline interpolation with second-order approximation accuracy, less segmentation and small operand by a conformal segmentation method, can inhibit the overshoot/undershoot problem of envelope fitting, provides a criterion of IMF component screening termination by the orthogonality and energy properties of the decomposed IMF, ensures the orthogonality and completeness of EMD decomposition, judges the number of noise contained in the screened IMF signal by the principle of mutual information to determine whether to filter the IMF signal, and increases the rapidity of an EMD algorithm; the threshold function is improved, and the threshold function combines the advantages of soft and hard thresholds to filter IMFs containing noise.

Description

Electrocardiosignal denoising method based on fusion of improved EMD and threshold method
Technical Field
The invention relates to an electrocardiosignal denoising method, in particular to an electrocardiosignal denoising method based on the fusion of an improved EMD and a threshold method, and belongs to the technical field of signal filtering.
Background
The electrocardiosignal is a typical unsteady and weak bioelectric signal and is widely applied to diagnosis and treatment of various heart diseases. The electrocardiosignals are often accompanied by very serious high-frequency and low-frequency noises, and the noise frequency band is often overlapped with the electrocardiosignal frequency band, so that the filtering pretreatment is difficult.
Conventional biomedical signal processing is primarily based on fourier theory. Fourier signal processing techniques are almost irreplaceable in the field of signal spectrum analysis and its associated signal processing fields such as data compression, signal detection, filtering, etc. However, the integral interval of the fourier transform is from positive infinity to negative infinity, and it cannot obtain the frequency spectrum content of the signal in a certain period of time. Due to the excellent time-frequency analysis characteristic and the capability of processing non-stationary random signals, the wavelet transform becomes an effective method for processing biomedical signals such as electrocardio and the like. Also, EMD has begun to be applied in the biomedical processing field due to its good adaptability in analyzing nonlinear and non-stationary signals. Such as electrocardiogram signal analysis, blood pressure signal de-noising, heartbeat signal analysis, etc., have been successfully applied.
Empirical Mode Decomposition (EMD) decomposes a signal into a sum of finite eigenmode functions (IMFs). The EMD decomposition fully considers the local scale characteristics of the signal, and each IMF component obtained in the way represents one scale characteristic of the original signal and contains the real physical information of the original signal. As a new self-adaptive signal time-frequency processing method, the EMD has wide application in the aspects of mechanical fault diagnosis, feature extraction, geophysical detection, medical analysis and the like, and the EMD method is also expanded to the field of two-dimensional signal processing. The method is well applied to the fields of less image edge detection, texture analysis, image fusion, image compression, image filtering and the like, and the effectiveness of the EMD is demonstrated.
The current EMD method has some defects and needs further research. Such as: the method comprises the following steps of research on an EMD quick algorithm, EMD modal aliasing problems, end point effects, envelope fitting problems of signals and criterion research on IMF component screening termination. In the method, the decomposition result directly influenced by the quality of the fitting of the envelope curve is completed by taking the extreme point as a known point in the original method adopting cubic spline function interpolation, and the fitting of the non-uniform point is caused by the non-uniformity of the distribution of the extreme point, so that the problem of extreme value undershoot or overshoot is caused, and a large error of decomposition is caused. The modal aliasing problem is a problem often encountered during EMD use, which is mainly caused by both intermittent events and signal interactions.
Aiming at the defects, the invention provides an electrocardiosignal denoising method based on the fusion of improved EMD and a threshold value method, and aims to improve the quality of the empirical mode decomposition of the electrocardiosignal and improve the effect of removing noise.
Disclosure of Invention
Aiming at the defects in the prior art, the invention aims to provide an improved EMD and threshold fusion technology to decompose electrocardiosignals and remove noise of the electrocardiosignals. The method solves the problem of modal aliasing by superposing white noise with different weight coefficients, solves the problem of end points by a method of a least square support vector machine, constructs upper and lower envelope lines of a signal by a method of conformal spline interpolation, and constructs cubic spline interpolation with second-order approximation precision, less segmentation and small operand by a conformal segmentation method. By means of orthogonality and energy properties of the decomposed IMF, a criterion of 'screening' termination of IMF components is provided, and orthogonality and completeness of EMD decomposition are guaranteed. The screened IMF signals are judged to contain noise according to the principle of mutual information to decide whether to carry out filtering processing on the IMF signals, and therefore the rapidity of the EMD algorithm is greatly improved. And an improved threshold function is provided, which combines the advantages of soft and hard thresholds to filter the IMF containing noise.
The invention aims to solve the problem of modal aliasing caused by EMD decomposition, and provides an electrocardiosignal denoising method based on the fusion of improved EMD and a threshold value method.
An electrocardiosignal denoising method based on the fusion of an improved EMD and a threshold method comprises the following steps:
step 1: performing improved EMD decomposition on the electrocardiosignal;
step 1.1: the method of adding white noise to a cardiac electrical signal can be expressed as:
Figure BDA0002259653930000021
wherein,
Figure BDA0002259653930000022
representing the signal after the positive coefficient white noise is added to the electrocardiosignal for the a time;
Figure BDA0002259653930000023
representing the signal after negative coefficient white noise corresponding to the positive coefficient is added to the electrocardiosignal for the a time; y (t) represents the cardiac signal; n (t) represents white noise with a mean of zero and unit variance; omegaaIs a weightAnd (4) the coefficient. The following steps 1.2 to 1.5 are processing methods for adding the positive coefficient white noise, and similarly, the processing method for adding the negative coefficient white noise only needs to change the "+" mark to the "-" mark, and the meaning indicates the corresponding situation for adding the negative coefficient white noise.
Step 1.2: in order to avoid the phenomenon of two-end divergence during signal interpolation fitting, namely an end effect. Extending two ends of the signal by using a least square support vector machine; the signal processed in step 1.1 is sampled, first, the signal is sampled
Figure BDA0002259653930000024
Sampling is carried out, and the signal sampling sequence is
Figure BDA0002259653930000025
And N is the number of sampling points. The training sample set is A { (h)1,g1),(h2,g2),···,(hl,gl) }, wherein: h isiIn order to input the vector, the vector is input,
Figure BDA0002259653930000026
giin order to output the vector, the vector is,
Figure BDA0002259653930000027
i is more than or equal to 1 and less than or equal to l; l is the number of training sample sets; test sample set B { (h)N-S+1,gN-S+1),(hN-S+2,gN-S+2),···,(hN,gN) }, wherein: s is more than or equal to 1 and is the number of test samples;
the output of the least squares support vector machine is
Figure BDA0002259653930000028
Wherein,
Figure BDA0002259653930000029
for the first signal after continuation, it will beAs new boundary points of the original data, in the same wayObtaining the extension value of the 2 nd data sequence
Figure BDA0002259653930000031
By analogy, all continuation sequences can be obtained according to the number of data points needing to be prolongedWherein M is the number of signal points extended right, the left extension method for a given data sequence is the same as the right extension method, and the left extension sequence is recorded as:
Figure BDA0002259653930000033
the final extended sequence is noted
Step 1.3: and constructing upper and lower envelope lines of the signal by a conformal spline interpolation method.
Step 1.3.1: sequence Z (t) obtained by comparison of step 1.2i) The relationship between the size of a certain point and the size of the adjacent points to judge whether the point is an extreme point. The specific operation is as follows:
when the endpoint is judged, and
Figure BDA0002259653930000035
the method comprises the following steps:
if it is
Figure BDA0002259653930000036
Then point is reached
Figure BDA0002259653930000037
Is the maximum point.
If it isThen point is reachedIs a minimum point.
If it isThen point is reached
Figure BDA00022596539300000311
Is the maximum point.
If it is
Figure BDA00022596539300000312
Then point is reached
Figure BDA00022596539300000313
Is a minimum point.
When judging other points than the above-mentioned end points:
if it is
Figure BDA00022596539300000314
I is more than or equal to M +1 and less than or equal to N + M-1, then point
Figure BDA00022596539300000315
Is the maximum point.
If it is
Figure BDA00022596539300000316
I is more than or equal to M +1 and less than or equal to N + M-1, then point
Figure BDA00022596539300000317
Is a minimum point.
Sequences obtainable by the above methods
Figure BDA00022596539300000318
The maximum point sequence of (a) is:
{xmax(t0),xmax(t1),...,xmax(tb) B is the number of maximum points;
the minimum point sequence is:
{xmin(t0),xmin(t1),...,xmin(tc) Where c is the number of maxima points.
Step 1.3.2: the upper envelope is constructed by the maximum point.
Two value points are inserted between every two maximum value points, and the definition is as follows:
wherein, tijTo be at time ti-1And tiThe j-th insertion point in between.
The sequences after insertion of the numerical points were:
{xmax(t0),xmax(t11),xmax(t12),xmax(t1),xmax(t21),xmax(t22),xmax(t3),...,xmax(tb)}
taking the above sequence as the control vertex of the T-B spline curve, the upper envelope curve, i.e. the T-B spline curve, can be expressed as follows:
Figure BDA0002259653930000041
wherein j represents the fitted j-th section of the T-B spline curve; u represents an upper envelope; h isij(T) is a T-B spline basis function when fitting T0To t1When a curve is segmented, hijThe expression of (t) is as follows:
Figure BDA0002259653930000042
Figure BDA0002259653930000043
Figure BDA0002259653930000044
Figure BDA0002259653930000045
wherein:
Figure BDA0002259653930000046
Figure BDA0002259653930000047
Figure BDA0002259653930000049
Figure BDA00022596539300000410
Figure BDA00022596539300000411
Figure BDA00022596539300000412
Figure BDA0002259653930000051
Figure BDA0002259653930000052
fitting t1:t2、t2:t3、···、tb-1:tbWhen a curve is segmented, hijThe expression of (t) is similar to that described above.
Step 1.3.3: the lower envelope is constructed by minimum points, the construction method is similar to the method for constructing the upper envelope, and the expression of the lower envelope is as follows:
Figure BDA0002259653930000056
wherein D represents a lower envelope;
step 1.4: averaging the upper envelope and the lower envelope:
Figure BDA0002259653930000057
wherein m isa(t) is the average of the upper and lower envelope lines obtained after the positive coefficient white noise is added for the a-th time;
calculating signals
Figure BDA0002259653930000058
And maDifference in (t):
Figure BDA0002259653930000059
if C (t) does not meet the IMF defined cutoff condition, repeating step 1.3; otherwise, extracting C (t) as
Figure BDA00022596539300000510
Figure BDA00022596539300000511
Adding IMF of positive coefficient white noise decomposed by the ith EMD for the ith; residual amount of
Figure BDA00022596539300000512
Step 1.5: and (3) taking the residual amount r (t) in the step 1.4 as a new signal, namely r (t) ═ y (t), and then screening through the steps 1.1 to 1.4 to obtain the next lower-frequency IMF, and stopping screening until the screened IMF component meets the following rule.
The stopping criteria for the termination of the IMF component "screening" are:
wherein: IMFi(t) is the IMF of the ith EMD decomposition, and n is the IMFi(t) length;
when ZSD < θ, sieving is stopped, preferably θ is 0.03.
Step 1.6: the signal after the above decomposition can be expressed as the following formula:
Figure BDA0002259653930000061
wherein: d is the number of IMFs after EMD decomposition of the electrocardiosignals,
Figure BDA0002259653930000062
adding IMF, r of positive coefficient white noise subjected to I EMD decomposition to the ad(t) is the residual error in the d-th decomposition.
Similar to the step 1.1 to the step 1.5, the EMD decomposition after adding the white noise with the negative coefficient to the electrocardiosignal is calculated, and the decomposed signal can be expressed as the following formula:
Figure BDA0002259653930000063
wherein,
Figure BDA0002259653930000064
the i-th EMD decomposed IMF with negative coefficient white noise was added for the a-th time.
Step 1.7: accumulating and averaging the IMF added with the white noise of the positive coefficient and the negative coefficient obtained in the step 1.6 to obtain:
Figure BDA0002259653930000065
wherein:
Figure BDA0002259653930000066
step 2: judging the linear relation between each IMF and the original signal through arranging mutual information, and judging the effective signal of the IMF obtained in the step 1, wherein the formula is as follows:
QI(IMF(k),y)＝S(IMF(k))+S(y)-S(IMF(k),y) (13)
wherein: s (IMF)(k)) And S (y) are respectively IMF(k)(t) and y (t) entropy of alignment, S (IMF)(k)) The calculation formula is as follows:
Figure BDA0002259653930000067
wherein p (-) represents the joint probability density of the permutation (-), and the calculation formula is as follows:
where i ═ 1,2, ·, n- (d-1) τ, n is the length of the time sequence, d is the given embedding dimension, τ is the time delay; # denotes the number of elements in the set,
Figure BDA0002259653930000069
is IMF(k)Is arranged in (pi) with respect toi) (ii) a Arrangement (n)i) Representing time series IMF(k)(t) mapping to a d-dimensional space and then mapping each state vector to a corresponding permutation sequence;
suppose a time series { IMF(k)(t) maps to d-dimensional phase space, defining its joint state vector as:
the state vector is used for representing the track of the ith time sequence, and a track state matrix mapped to a state space can be obtained through the definition:
obtaining an arrangement sequence by comparing the magnitude relation of the adjacent values of the row vector of each state matrix, namely obtaining an arrangement matrix of the track matrix:
s (y) and S (IMF)(k)) The calculation of (a) is similar.
S(IMF(k)Y) represents IMF(k)(t) and y (t) joint permutation entropy, the calculation formula is as follows:
Figure BDA0002259653930000073
wherein p (pi)i,πj) Represents arrangement (pi)i,πj) The calculation formula is as follows:
Figure BDA0002259653930000074
where the subscript i ═ 1,2,. cndot.n- (d-1) τ,a pair of state space trajectory matrices, the corresponding arrangement of which is (pi)i,πj) (ii) a Arrangement (n)i,πj) Representation of IMF(k)(t) and y (t) an arrangement order corresponding to each state vector after the time series is mapped to the d-dimensional space, which is expressed as follows: suppose a time series { IMF(k)(t) and { y (t) } are mapped to d-dimensional phase space, defining its joint state matrix as:
Figure BDA0002259653930000076
the state vector is used for representing the track of the ith time sequence, and a track-like matrix mapped to the state space can be obtained through the definition:
Figure BDA0002259653930000077
obtaining an arrangement sequence by comparing the magnitude relation of the adjacent values of the row vector of each state matrix, namely obtaining an arrangement matrix of the track matrix:
Figure BDA0002259653930000078
and judging the degree of noise interference of each IMF according to the QI value of the mutual information of each IMF and the original signal, wherein the more the noise quantity is, the smaller the QI value is. And selecting the IMF with a small obvious QI value for denoising.
And step 3: carrying out threshold denoising processing on the IMF which is selected in the step 2 and needs to be denoised, wherein the threshold formula is as follows:
Figure BDA0002259653930000081
wherein: adjustable parameter
Figure BDA0002259653930000082
The adjustable parameter alpha is a positive number, and the adjustable parameter m is a positive number; threshold lambda1The calculation formula of (a) is as follows:
threshold value
Figure BDA0002259653930000083
σ is the standard deviation of noise in each IMF component, σ ═ mean (IMF)i) 0.6745, median (-) is the median;
threshold lambda2The calculation formula of (a) is as follows:
and 4, step 4: and (3) performing signal reconstruction on the IMF left in the step (2) and the signal processed in the step (3), wherein a calculation formula is as follows:
Figure BDA0002259653930000085
where y' (t) represents the denoised signal.
So far, the electrocardiosignal denoising process based on the fusion of the improved EMD and the threshold value method is completed from the step 1 to the step 4.
Advantageous effects
Compared with the prior art, the electrocardiosignal denoising method based on the fusion of the improved EMD and the threshold method has the following beneficial effects:
1. the method effectively solves the problem of mode aliasing caused by the interaction of two signals in the EMD by adopting a method of respectively adding positive white noise and negative white noise with different weight coefficients for multiple times, and improves the frequency resolution of the algorithm (the frequency resolution refers to the capability of distinguishing two adjacent frequency components);
2. the method uses a least square support vector machine method to carry out continuation and windowing on the signal, thereby avoiding the distortion phenomenon caused by the adoption of a cubic spline interpolation method to carry out fitting on the signal during decomposition;
3. the method uses a conformal spline interpolation method to construct the upper envelope line and the lower envelope line of a signal, and utilizes a conformal segmentation method to construct a cubic spline interpolation method with second-order approximation accuracy, less segmentation and small operand to inhibit overshoot/undershoot problems of envelope fitting, thereby avoiding the problem that interpolation errors caused by the traditional interpolation method are continuously accumulated along with the continuous proceeding of a decomposition process to cause serious errors.
4. The method directly adopts a threshold value method to denoise the IMF, and the method effectively reduces the operation time. The threshold function is improved, the details and the edge characteristics of the useful signal are better reserved by the threshold function, an effective signal is better provided on the basis of realizing signal noise reduction, the amplitude of the signal is guaranteed, the denoising effect is improved, and the operation time is reduced.
5. The method has simple flow and easy realization, and can be used for the design work of related software in the field of electrocardiosignal analysis.
Drawings
FIG. 1 is a schematic representation of the steps of the present invention;
FIG. 2 is a flow chart of the present invention;
FIG. 3 is a detailed diagram of the EMD decomposition of the present invention;
Detailed Description
The invention is explained in detail below with reference to the figures and examples, but the specific embodiments of the invention are not limited thereto.
The data in the embodiment are derived from the sampling data of the patient in the cardiovascular medical department of the general hospital of people liberation of China, an electrocardiograph is adopted to acquire the electrocardiosignal of the patient, the signal duration is 5s, the sampling frequency is 360Hz, and the electrocardiosignal data of the patient under the condition of daily free activity (non-violent movement) is acquired.
The embodiment explains a process of applying the electrocardiosignal denoising method based on the fusion of the improved EMD and the threshold value method to an electrocardiosignal denoising scene.
Fig. 1 is a schematic diagram of the steps of the method, and fig. 2 is a specific flowchart, and it can be seen from the diagram that the electrocardiosignal is firstly collected, and then the following steps are carried out:
step A: performing improved EMD on the acquired signals; fig. 3 is a detailed diagram of the EMD decomposition of the present invention, which is specifically as follows:
step A1: adding positive omega to the acquired signal y (t)1White noise n (t) with a zero mean value and unit variance;
in the embodiment, the white noise amplitude is 0.05 times of standard deviation of the electrocardiosignal, i.e. ω10.5, the signal after adding noise is recorded as x1(t)；
Step A2: using least square support vector machine to x obtained in step A.11(t) continuation of the signal;
for x1(t) sampling to obtain a sampling sequence { x (1), x (2), x (3), ·, x (N) }, wherein N is the number of sampling points, and the training sample set is B { (h)1,g1),(h2,g2),···,(hl,gl)}。
Specifically, in this embodiment, the number of sampling points is 2000, the training sample is 100, and the training sample can be obtained as follows: the input samples are 100 × 1900 dimensional data.
The first predictor of continuation to the right is x (N + 1); then, x (N +1) is used as a new boundary point of the original data, and the 2 nd data sequence continuation value x (N +2) is obtained by the same method. By analogy, all continuation sequences x (N +1), x (N +2) …, x (N + M) can be obtained according to the number of data points required to be continued, and the forward continuation method for a given data sequence is the same as the backward continuation method.
100 sampling points are respectively extended from two end points of the signal to two sides, namely, an extended sequence with M being 100 is as follows: { x (1), x (2), x (3),. cndot., x (2200) }.
Step A3: constructing a signal envelope line by utilizing a conformal spline interpolation method;
specifically, in this embodiment, the magnitude relationship between each point in the sequence obtained in step a.2 and its adjacent points is compared to determine whether the point is an extreme point, and according to the comparison rule, the compared maximum sequence is { x }max(t0),xmax(t1),...,xmax(tb) B is the number of maximum points; inserting two values between two maximum values of the formula (2) as control vertexes of the T-B spline, and constructing an upper envelope of the signal by using the formula (3) and recording the upper envelope as the upper envelopeThe method of constructing the lower envelope is similar to that of the upper envelope, and is described
Step A4: calculating the mean value m of the upper and lower envelope curves using equation (6)a(t)；
Step A5: calculating the signal x obtained in step A.2 using equation (7)1' (t) and mean value ma(t) the difference;
judging whether C (t) meets the IMF definition interception condition according to the formula (8), namely meeting the IMF definition interception conditionWhen ZSD is less than 0.3, extracting C (t) as IMF; otherwise using the formula
Figure BDA0002259653930000103
Calculating the residual quantity, making the residual quantity be an original signal, and repeating the steps A1-A5 until the screening is stopped when the residual component is a monotonic function.
Step A6: changing the weight factor omega of the step A.1, and carrying out screening of the step A.1 to the step A.5 for m times;
specifically, in this embodiment, the white noise amplitude may be 0.05 times, 0.06 times and 0.04 times of the standard deviation of the electrocardiographic signal, and the sum of the white noise amplitudes is 100 times;
step A7: all IMF components of the signal subjected to EMD can be obtained by using the formula (12);
and B: performing effective signal judgment on each obtained IMF through arranging mutual information;
step B1: solving the permutation entropy of the time series;
specifically to this embodiment, the embedding dimension is 1001 and the time delay is 1. Calculating time series IMF according to Shannon entropy of permutation probability density function(k)Permutation entropy S (IMF) of (t) and y (t)(k)) And S (y); calculating the joint time series (IMF) using equation (18)(k) (t), y (t) joint permutation entropy S (IMF)(k),y)。
Step B2: calculating permutation mutual information by using a formula (13);
step B3: and judging the degree of noise interference of each IMF according to the QI value of the mutual information of each IMF and the original signal, wherein the more the noise quantity is, the smaller the QI value is. And selecting the IMF with a small obvious QI value for denoising.
And C: b, performing threshold denoising treatment on the IMF which is selected in the step B and needs to be denoised;
specifically, in this embodiment, the selected IMF is subjected to denoising processing using equation (23).
Step D: and C, performing signal reconstruction on the IMF remained in the step B and the signal processed in the step C.
The signal is reconstructed by equation (24).
Therefore, the electrocardiosignal denoising method based on the fusion of the improved EMD and the threshold method is completed from the step A to the step D.

Claims (10)

1. An electrocardiosignal denoising method based on the fusion of an improved EMD and a threshold method comprises the following steps:
A. performing an improved EMD decomposition of the acquired signal, in particular:
A1. adding omega to the collected signal y (t)aThe multiplied mean value is zero, and the positive coefficient white noise n (t) of unit variance is obtained to obtain a signal
Figure FDA0002259653920000011
Subscript a indicates the a-th addition of white noise;
A2. signal pair using least squares support vector machineCarrying out continuation at two ends to obtain a continuation sequence;
A3. constructing the upper envelope of the continuation sequence by using a conformal spline interpolation method
Figure FDA0002259653920000013
And lower envelope
Figure FDA0002259653920000014
j represents the j section of T-B spline fitting curve;
A4. calculating the mean value of the upper envelope line and the lower envelope line obtained after white noise is added for the a-th timeAnd a signal
Figure FDA0002259653920000016
And maDifference of (t)
Figure FDA0002259653920000017
If C (t) does not meet the IMF defined cutoff condition, repeating step A3; otherwiseExtracting C (t) asThe IMF of the ith EMD decomposition of the added positive coefficient white noise is shown; calculating the residual amount
Figure FDA0002259653920000019
A5. Taking the residual amount r (t) of the step a4 as a new signal, namely r (t) ═ y (t), and then screening through the steps a 1-a 4 to obtain the next IMF with lower frequency, and stopping screening until the screened IMF component meets ZSD < theta, wherein theta is a set threshold value, and the ZSD calculation formula is as follows:
Figure FDA00022596539200000110
IMFi(t) is the IMF of the ith EMD decomposition, and n is the IMFi(t) length;
A6. signals decomposed from A1 to A5 are shown as
Figure FDA00022596539200000111
Wherein d is the number of IMFs after EMD decomposition of the electrocardiosignals,
Figure FDA00022596539200000112
adding IMF, r of positive coefficient white noise subjected to I EMD decomposition to the ad(t) is the residual error when the d-th decomposition is performed;
similar to the steps A1-A5, EMD decomposition is performed on the signal after negative coefficient white noise is added to the electrocardiosignal, and the decomposed signal is expressed as
Figure FDA00022596539200000113
Wherein,
Figure FDA00022596539200000114
adding IMF of negative coefficient white noise decomposed by the ith EMD for the ith time;
A7. adding the white noise with the positive coefficient and the negative coefficient obtained in the step A6 and then performing IMF accumulation and averaging to obtain
Figure FDA00022596539200000115
Wherein,
Figure FDA00022596539200000116
B. judging the linear relation between each IMF and the original signal through the permutation mutual information, and judging the effective signal of the IMF obtained in the step A;
C. b, performing threshold denoising treatment on the IMF which is selected in the step B and needs to be denoised;
D. and C, performing signal reconstruction on the IMF remained in the step B and the signal processed in the step C.
2. The method for denoising an electrocardiographic signal according to claim 1, wherein a2 includes: to pair
Figure FDA0002259653920000021
Sampling is carried out to obtain a sampling sequence of
Figure FDA0002259653920000022
N is the number of sampling points; the output of the least squares support vector machine is
Figure FDA0002259653920000023
Wherein,
Figure FDA0002259653920000024
is the first signal after continuation; then will be
Figure FDA0002259653920000025
Obtaining the 2 nd data sequence continuation value as a new boundary point of the original data
Figure FDA0002259653920000026
By analogy, all continuation sequences are obtained according to the number of data points needing to be prolonged
Figure FDA0002259653920000027
Figure FDA0002259653920000028
Wherein M is the number of signal points extended to the right; left continuation for a given data sequence
Figure FDA0002259653920000029
The final continuation sequence is
Figure FDA00022596539200000210
3. The method for denoising an electrocardiographic signal according to claim 1 or 2, wherein a3 comprises: comparing the magnitude relation between a certain point in the continuation sequence and the left and right adjacent points to judge whether the point is an extreme point; constructing an upper envelope by means of maximum points
Figure FDA00022596539200000211
Constructing the lower envelope by minimum points
Figure FDA00022596539200000212
4. The method for denoising an electrocardiographic signal according to claim 3, wherein the method for judging whether the point is an extreme point comprises: for end points
Figure FDA00022596539200000213
And
Figure FDA00022596539200000214
if it is
Figure FDA00022596539200000215
Then point is reached
Figure FDA00022596539200000216
Is the maximum point, otherwise, the point
Figure FDA00022596539200000217
Is a minimum value point; if it is
Figure FDA00022596539200000218
Then point is reachedIs the maximum point, otherwise, the point
Figure FDA00022596539200000220
Is a minimum value point;
for points other than the above endpoints: if it is
Figure FDA00022596539200000221
Then point is reached
Figure FDA00022596539200000222
Is a maximum point; if it is
Figure FDA00022596539200000223
Then point is reached
Figure FDA00022596539200000224
Is a minimum value point, wherein, i is more than or equal to-M +1 and less than or equal to N + M-1.
5. The method for denoising an electrocardiographic signal according to claim 1, wherein the step B comprises:
B1. computing time series IMF(k)Permutation entropy S (IMF) of (t) and y (t)(k)) And S (y);
B2. calculating arrangement mutual information;
B3. and judging the degree of noise interference of each IMF according to the mutual information QI value of each IMF and the original signal, and selecting the IMF with the small QI value for denoising.
6. The electrocardiosignal denoising method of claim 1 or 5, wherein the calculation formula of the arrangement mutual information is:
QI(IMF(k),y)＝S(IMF(k))+S(y)-S(IMF(k)y) where S (IMF)(k)) And S (y) are respectively IMF(k)(t) and y (t) entropy of alignment, S (IMF)(k)Y) is IMF(k)(t) and y (t) joint permutation entropy.
7. The method of denoising of cardiac signals of claim 6, wherein S (IMF)(k)) The calculation formula of the permutation entropy of (a) is:
Figure FDA0002259653920000031
where p (-) represents the joint probability density of the permutation (-).
8. The method of denoising an electrocardiographic signal of claim 6, wherein the IMF is(k)Joint permutation entropy S (IMF) of (t) and y (t)(k)Y) is calculated asWherein, is arranged (pi)i,πj) Representation of IMF(k)(t) and y (t) a permutation order corresponding to each state vector after the time series is mapped to the d-dimensional space, and the probability density is combinedThe subscripts i ═ 1,2, ·, n- (d-1) τ,
Figure FDA0002259653920000034
is a pair of state space trajectory matrices, which correspond to an arrangement of (pi)i,πj)。
9. The method for denoising an electrocardiographic signal according to claim 1, wherein the threshold value in step C is calculated by:
Figure FDA0002259653920000035
wherein the adjustable parameter m is a positive number, and the adjustable parameter m is a positive number
Figure FDA0002259653920000036
Alpha is a positive number, and alpha is a positive number,
Figure FDA0002259653920000037
σ is the standard deviation of noise in each IMF component, σ ═ mean (IMF)i) 0.6745, median (-) is the median;
Figure FDA0002259653920000038
10. the method for denoising the cardiac signal according to claim 1, wherein the signal reconstruction calculation formula in step D is
Figure FDA0002259653920000039
Where y' (t) represents the denoised signal.
CN201911070409.XA 2019-11-04 2019-11-04 Electrocardiosignal denoising method based on fusion of improved EMD and threshold method Active CN110680308B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201911070409.XA CN110680308B (en) 2019-11-04 2019-11-04 Electrocardiosignal denoising method based on fusion of improved EMD and threshold method

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201911070409.XA CN110680308B (en) 2019-11-04 2019-11-04 Electrocardiosignal denoising method based on fusion of improved EMD and threshold method

Publications (2)

Publication Number Publication Date
CN110680308A true CN110680308A (en) 2020-01-14
CN110680308B CN110680308B (en) 2021-04-20

Family

ID=69116188

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201911070409.XA Active CN110680308B (en) 2019-11-04 2019-11-04 Electrocardiosignal denoising method based on fusion of improved EMD and threshold method

Country Status (1)

Country Link
CN (1) CN110680308B (en)

Cited By (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111643068A (en) * 2020-05-07 2020-09-11 长沙理工大学 Electrocardiosignal denoising algorithm, electrocardiosignal denoising equipment and storage medium based on EMD and energy thereof
CN112364291A (en) * 2020-11-17 2021-02-12 哈工大机器人（合肥）国际创新研究院 Pre-filtering extreme point optimization set empirical mode decomposition method and device
CN113253300A (en) * 2021-06-18 2021-08-13 湖南国天电子科技有限公司 Optical echo signal denoising method and system for laser cloud measuring radar machine
CN113503915A (en) * 2021-07-01 2021-10-15 北京工商大学 Pollutant data prediction method and device, storage medium and electronic equipment
CN113616213A (en) * 2021-07-29 2021-11-09 山东大学 Electrocardiosignal denoising method, equipment and storage medium based on BP neural network and improved EMD method
CN114469124A (en) * 2022-01-30 2022-05-13 北京理工大学 Method for identifying abnormal electrocardiosignals in motion process
CN115881276A (en) * 2023-02-22 2023-03-31 合肥工业大学 Time/frequency dual barcode feature image generation method for electrocardiosignal and storage medium
CN117420346A (en) * 2023-12-19 2024-01-19 东莞市兴开泰电子科技有限公司 Circuit protection board overcurrent value detection method and system
CN117435907A (en) * 2023-08-14 2024-01-23 北京体育大学 Sports record detection method and device based on data analysis

Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101853242A (en) * 2010-06-23 2010-10-06 哈尔滨工业大学 Equipment or system built-in test signal false-alarm filtering method based on empirical mode decomposition
CN104702244A (en) * 2013-12-05 2015-06-10 中国科学院深圳先进技术研究院 Adaptive filter for filtering power frequency interference in electromyography signal based on EEMD (Ensemble Empirical Mode Decomposition) algorithm
CN105030232A (en) * 2015-06-30 2015-11-11 广东工业大学 Baseline drift correction method for electrocardiosignal
CN107885940A (en) * 2017-11-10 2018-04-06 吉林大学 A kind of signal characteristic extracting methods for distributed optical fiber vibration sensing system
CN108288058A (en) * 2018-04-12 2018-07-17 大连理工大学 A kind of improved wavelet threshold knee joint swinging signal Denoising Algorithm
CN109589114A (en) * 2018-12-26 2019-04-09 杭州电子科技大学 Myoelectricity noise-eliminating method based on CEEMD and interval threshold
CN109785854A (en) * 2019-01-21 2019-05-21 福州大学 The sound enhancement method that a kind of empirical mode decomposition and wavelet threshold denoising combine
CN109998527A (en) * 2019-04-09 2019-07-12 湖北工业大学 A kind of heart disease detection method based on multi-scale entropy

Patent Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101853242A (en) * 2010-06-23 2010-10-06 哈尔滨工业大学 Equipment or system built-in test signal false-alarm filtering method based on empirical mode decomposition
CN104702244A (en) * 2013-12-05 2015-06-10 中国科学院深圳先进技术研究院 Adaptive filter for filtering power frequency interference in electromyography signal based on EEMD (Ensemble Empirical Mode Decomposition) algorithm
CN105030232A (en) * 2015-06-30 2015-11-11 广东工业大学 Baseline drift correction method for electrocardiosignal
CN107885940A (en) * 2017-11-10 2018-04-06 吉林大学 A kind of signal characteristic extracting methods for distributed optical fiber vibration sensing system
CN108288058A (en) * 2018-04-12 2018-07-17 大连理工大学 A kind of improved wavelet threshold knee joint swinging signal Denoising Algorithm
CN109589114A (en) * 2018-12-26 2019-04-09 杭州电子科技大学 Myoelectricity noise-eliminating method based on CEEMD and interval threshold
CN109785854A (en) * 2019-01-21 2019-05-21 福州大学 The sound enhancement method that a kind of empirical mode decomposition and wavelet threshold denoising combine
CN109998527A (en) * 2019-04-09 2019-07-12 湖北工业大学 A kind of heart disease detection method based on multi-scale entropy

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
张红 等: "《一种用CEEMDAN和排列熵去除脉搏信号噪声的方法》", 《中国科技论文》 *
王亚东等: "《基于LS-SVM回归对EMD中端点效应的研究》", 《河南城建学院学报》 *

Cited By (14)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111643068A (en) * 2020-05-07 2020-09-11 长沙理工大学 Electrocardiosignal denoising algorithm, electrocardiosignal denoising equipment and storage medium based on EMD and energy thereof
CN111643068B (en) * 2020-05-07 2023-05-23 长沙理工大学 Electrocardiogram signal denoising algorithm, equipment and storage medium based on EMD and energy thereof
CN112364291A (en) * 2020-11-17 2021-02-12 哈工大机器人（合肥）国际创新研究院 Pre-filtering extreme point optimization set empirical mode decomposition method and device
CN112364291B (en) * 2020-11-17 2024-05-14 哈工大机器人（合肥）国际创新研究院 Empirical mode decomposition method and device for pre-filtering extreme point optimization set
CN113253300A (en) * 2021-06-18 2021-08-13 湖南国天电子科技有限公司 Optical echo signal denoising method and system for laser cloud measuring radar machine
CN113503915A (en) * 2021-07-01 2021-10-15 北京工商大学 Pollutant data prediction method and device, storage medium and electronic equipment
CN113616213B (en) * 2021-07-29 2023-02-28 山东大学 Electrocardiosignal denoising method, electrocardiosignal denoising device and electrocardiosignal denoising storage medium based on BP neural network and improved EMD method
CN113616213A (en) * 2021-07-29 2021-11-09 山东大学 Electrocardiosignal denoising method, equipment and storage medium based on BP neural network and improved EMD method
CN114469124A (en) * 2022-01-30 2022-05-13 北京理工大学 Method for identifying abnormal electrocardiosignals in motion process
CN114469124B (en) * 2022-01-30 2024-04-09 北京理工大学 Method for identifying abnormal electrocardiosignals in movement process
CN115881276A (en) * 2023-02-22 2023-03-31 合肥工业大学 Time/frequency dual barcode feature image generation method for electrocardiosignal and storage medium
CN117435907A (en) * 2023-08-14 2024-01-23 北京体育大学 Sports record detection method and device based on data analysis
CN117420346A (en) * 2023-12-19 2024-01-19 东莞市兴开泰电子科技有限公司 Circuit protection board overcurrent value detection method and system
CN117420346B (en) * 2023-12-19 2024-02-27 东莞市兴开泰电子科技有限公司 Circuit protection board overcurrent value detection method and system

Also Published As

Publication number Publication date
CN110680308B (en) 2021-04-20

Similar Documents

Publication Publication Date Title
CN110680308B (en) Electrocardiosignal denoising method based on fusion of improved EMD and threshold method
Aqil et al. ECG Signal Denoising by Discrete Wavelet Transform.
Alfaouri et al. ECG signal denoising by wavelet transform thresholding
CN109102550B (en) Full-network low-dose CT imaging method and device based on convolution residual error network
CN104367316B (en) Denoising of ECG Signal based on morphologic filtering Yu lifting wavelet transform
Strasser et al. Motion artifact removal in ECG signals using multi-resolution thresholding
CN111616697B (en) Electrocardiosignal denoising algorithm based on new threshold function wavelet transform
Mithun et al. A wavelet based technique for suppression of EMG noise and motion artifact in ambulatory ECG
WO2017148451A1 (en) Smooth wavelet transform-based method and system for filtering out electromyography interference
Kumar et al. Optimal selection of wavelet function and decomposition level for removal of ECG signal artifacts
CN107341769A (en) Electrocardiosignal denoising method and system
Bhateja et al. A novel approach for suppression of powerline interference and impulse noise in ECG signals
Singh et al. ECG signal denoising based on empirical mode decomposition and moving average filter
Abbaspour et al. ECG artifact removal from surface EMG signal using an automated method based on wavelet-ICA
Lu et al. Model-based ECG denoising using empirical mode decomposition
CN116584960A (en) Diaphragmatic electromyographic signal noise reduction method
CN106236075B (en) A kind of noise-reduction method applied to portable electrocardiograph institute thought-read electrograph
CN111493821A (en) PPG signal real-time denoising method based on MODWT and median filtering
Boucheham et al. Piecewise linear correction of ECG baseline wander: a curve simplification approach
Jain et al. An adaptive method for shrinking of wavelet coefficients for phonocardiogram denoising
CN114098656A (en) Signal noise reduction method and system based on empirical mode decomposition and bit plane conversion
Stalin Subbiah et al. Reduction of noises in ECG signal by various filters
Dubey et al. Two-stage nonlocal means denoising of ECG signals
CN110236528B (en) Method and device for acquiring respiratory information
Kaur et al. Adaptive wavelet thresholding for noise reduction in electrocardiogram (ECG) signals

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant