The electrocardiogram (ECG) records the electrical activity of the heart，which is a noninvasively recording produced by an electrocardiographic device and collected by skin electrodes placed at designated locations on the body. The ECG signal is characterized by six peaks and valleys, which are traditionally labeled P, Q, R, S, T, and U, shown in figure 1.
It has been used extensively for detection of heart disease. ECG is non-stationary bioelectrical signal including valuable clinical information, but frequently the valuable clinical information is corrupted by various kinds of noise. The main sources of noise are: power-line interference from 50–60 Hz pickup and harmonics from the power mains; baseline wanders caused by variable contact between the electrode and the skin and respiration; muscle contraction form electromyogram (EMG) mixed with the ECG signals; electromagnetic interference from other electronic devices and noise coupled from other electronic devices, usually at high frequencies. The noise degrades the accuracy and precision of an analysis. Obtaining true ECG signal from noisy observations can be formulated as the problem of signal estimation or signal denoising. So denoising is the method of estimating the unknown signal from available noisy data. Generally, excellent ECG denoising algorithms should have the following properties: Ameliorate signal-to-noise ratio (SNR) for obtaining clean and readily observable signals; Preserve the original characteristic waveform and especially the sharp Q, R, and S peaks, without distorting the P and T waves.
A lot of methods have been proposed for ECG denoising. In general both linear and nonlinear filters are presented, such as elliptic filter, median filter, Wiener filter and wavelet transform etc. These methods have some drawbacks. They remove not only noise but also the high frequency components of non-stationary signals. In the worse they can remove the characteristic points of signals that are crucial for successful detection of waveform. In recent years wavelet transform (WT) has become favourable technique in the field of signal processing. Donoho et al proposed the denoising method called “wavelet shrinkage”; it has three steps: forward wavelet transform, wavelet coefficients shrinkage at different levels and the inverse wavelet transform, which work in denoising the signals such as Universal threshold, SureShrink, Minimax. Wavelet shrinkage methods have been successful in denoising ECG signals (Agante, P.M&Marques J.P, 1999; Brij N. Singh & Arvind.K, 2006). A New wavelet shrinkage method for denoising of biological signals is proposed based on a new thresholding filter (Prasad V.V.K.D.V; Siddaiah P; Rao BP,2008).De-noising using traditional DWT has a translation variance problem which results in Pseudo-Gibbs phenomenon in the Q and S waves, so the following algorithms tried to solve this problem: used cyclic shift tree de-noising technique for reducing white Gaussian noise or random noise, EMG noise and power line interference (Kumari, R.S.S. et al,2008).The selected optimal wavelets basis has been investigated with suitable shrinkage method to de-noise ECG signals, not only it obtains higher SNR, but preserves the peaks of R wave in ECG(Suyi Li. et al,2009).Scale-dependent threshold methods are successively proposed. A new thresholding procedure is proposed based on wavelet denoising using subband dependent threshold for ECG signals: The S-median-DM and S-median thresholds (Poornachandra.S, 2008).
In this work, in order to enhance ECG, the new adaptive shrunken denoising method based on Ensemble Empirical Mode Decomposition (EEMD) is presented that has a good influence in enhancing the SNR, and also in terms of preserving the original characteristic waveform. The paper is organized as follows: section 2 introduces Empirical Mode Decomposition (EMD) and EEMD is studied in section 3. EMD is a relatively new, data-driven adaptive technique used to decompose ECG signal into a series of Intrinsic Mode Functions (IMFs). The EEMD overcomes largely the mode mixing problem of the original EMD by adding white noise into the targeted signal repeatedly and provides physically unique decompositions. Wavelet shrinkage is studied in section 4; the wavelet shrinkage denoising method is simply signal extraction from noisy signal via wavelet transform. It has been shown to have asymptotic near-optimality properties over a wide class of functions. The crucial points are the selections of threshold value and thresholding function. The generalized threshold function is build. Computationally exact formulas of bias 、variance and risk of generalized threshold function are derived. Section 5 concentrates on adaptive threshold values based on EEMD.Noisy signal is decomposed into a series of IMFs, and then the threshold values are derived by the noise energies of each IMFs. To evaluate the performance of the algorithm, Test signal and Clinic noisy ECG signals are processed in section 6. The results show that the novel adaptive threshold denoising method can achieve the optimal denoising of the ECG signal. Conclusions are presented in section7.
2. Empirical mode decomposition
EMD has recently been proposed by N.E. Huang in 1998 which is developed as a data-driven tool for nonlinear and non-stationary signal processing. EMD can decompose signal into a series of IMFs subjected to the following two conditions:
In the whole dataset, the number of extrema and the number of zero-crossing must either be equal or differ at most by one.
At any time, the mean value of the envelope of the local maxima and the envelope of the local minima must be zero.
Figure.2 shows a classical IMF. The IMFs represent the oscillatory modes embedded in signal. Each IMF actually is a zero mean monocomponent AM-FM signal with the following form:
with time varying amplitude envelop and phase. The amplitude and phase have both physically and mathematically meaning.
Most signals include more than one oscillatory mode, so they are not IMFs. EMD is a numerical sifting process to disintegrate empirically a signal into a finite number of hidden fundamental intrinsic oscillatory modes, that is, IMFs.The sifting process can be separated into following steps:
Finding all the local extrema, including maxima and minima; then connecting all the maxima and minima of signal x(t) using smooth cubic splines to get its upper envelopeand lower envelope.
Subtracting mean of these two envelopes from the signal to get their difference:.
Regarding theas the new data and repeating steps 1 and 2 until the resulting signal meets the two criteria of an IMF, defined as. The first IMF contains the highest frequency component of the signal. The residual signalis given by.
Regarding as new data and repeating steps (1) (2) (3) until extracting all the IMFs. The sifting procedure is terminated until the Mth residue becomes less than a predetermined small number or becomes monotonic.
The original signal x (t) can thus be expressed as following:
is an IMF where j represents the number of corresponding IMF and is residue.
The EMD decomposes non-stationary signals into narrow-band components with decreasing frequency. The decomposition is complete, almost orthogonal, local and adaptive. All IMFs form a completely and “nearly” orthogonal basis for the original signal. The basis directly comes from the signal which guarantees the inherent characteristic of signal and avoids the diffusion and leakage of signal energy. The sifting process eliminates riding waves, so each IMF is more symmetrical and is actually a zero mean AM-FM component.
The major disadvantage of EMD is the so-called mode mixing effect. For example, the simulated signal is defined as follows:
The signal is composed of sine wave and impulse functions, shown as figure3.It is decomposed into a series of IMFs by EMD, illustrated as figure 4. The decomposition is polluted by mode mixing, which indicates that oscillations of different time scales coexist in a given IMF, or that oscillations with the same time scale have been assigned to different IMFs.
3. Ensemble empirical mode decomposition
Ensemble EMD (EEMD) was introduced to remove the mode-mixing effect. The EEMD overcomes largely the mode mixing problem of the original EMD by adding white noise into the targeted signal repeatedly and provides physically unique decompositions when it is applied to data with mixed and intermittent scales.
The EEMD decomposing process can be separated into following steps:
1. Add a white noise seriesto the targeted data, the noise must be zero mean and variance constant, so.
2. Decompose the data with added white noise into Intrinsic Mode Functions (IMFs) and residue rn
3. Repeat step 1 and step 2 N times, but with different white noise serried wi(t) each time, so
4. Obtain the ensemble means of corresponding IMFs of the decompositions as the final result. Each IMF is obtained by decomposed the targeted signal.
This new approach utilizes the full advantage of the statistical characteristics uniform distribution of frequency of white noise to improve the EMD method. The above signal is decomposed into a series of IMFs by EEMD, which is shown in figure 5. Through adding white noise into the targeted signal makes all scaled continues to avoid mode mixing phenomenon. Comparing the IMF component of the same level, EEMD has more concentrated and band limited components.
We consider the following model of a discrete noisy signal:
The vector represents noisy signal and is an unknown original clean signal. is independent identity distribution Gaussian white noise with mean zero and unit variance. For simplicity, we assume intensity of noise is one. The step of wavelet shrinkage is defined as follows:
Apply discrete wavelet transform to observed noisy signal.
Estimate noise and threshold value, thresholding the wavelet coefficients of observed signal.
Apply the inverse discrete wavelet transform to reconstruct the signal.
The wavelet shrinkage method relies on the basic idea that the energy of signal will often be concentrated in a few coefficients in wavelet domain while the energy of noise is spread among all coefficients in wavelet domain. Therefore, the nonlinear shrinkage function in wavelet domain will tend to keep a few larger coefficients over threshold value that represent signal, while noise coefficients down threshold value will tend to reduce to zero.
In the wavelet shrinkage, how to select the threshold function and how to select the threshold value are most crucial. Donoho introduced two kinds of thresholding functions: hard threshold function and soft threshold function.
Hard threshold function (8) results in larger variance and can be unstable because of discontinuous function. Soft threshold function (9) results in unnecessary bias due to shrinkage the large coefficients to zero. We build the generalized threshold function:
When m is even number:
When m is odd number:
When m=1, it is soft threshold function; when m=, it is hard threshold function. When m=2 it is Non-Negative Garrote threshold function. We show slope signal as an example, Figure.6 graphically shows generalized threshold functions for different m. It can be clearly seen that when the coefficient is small, the smaller m is, the closer the generalized function is to the soft threshold function; when the coefficient is big, the bigger m is, the closer the generalized function is to the hard threshold function. As, generalized threshold function achieves a compromise between hard and soft threshold function. With careful selection of m, we can achieve better denoising performance.
We derived the exact formula of mean, bias, variance andrisk for generalized threshold function.
andare density and probability function of standard Gaussian random variable respectively. Then:
,,,are the mean, bias, variance and risk of generalized threshold function When m is 1, 2, , they are the mean, bias, variance and risk of the risk of soft, Non-Negative Garrote, hard threshold functions, respectively.
Soft threshold function provides smoother results in comparison with the hard threshold function; however, the hard threshold function provides better edge preservation in comparison with the soft threshold function. The hard threshold function is discontinuous and this leads to the oscillation of denoised signal. Soft threshold function tends to have bigger bias because of shrinkage, whereas hard threshold function tends to have bigger variance because of discontinuity. Non-Negative Garrote threshold function is the trade-off between the hard and soft threshold function. Firstly it is continuous; secondly the shrinkage amplitude is smaller than the soft threshold function.
5. Adaptive threshold values based on EEMD
Threshold value is a parameter that controls the bias and vriance tradeoff of the risk. If it is too small, the estimators tend to overfit the data, then result is close to the input and the estimate bias is reduced but the variance is increased. If the threshold value is too large, a lot of wavelet coefficients are set zero and the estimators tend to underfit the data; the estimate variance is reduced but the bias is increased. The optimal threshold value is the best compromise between variance and bias and it should minimize the risk of the results as compared with noise-free data.
Several methods have been proposed for the determinations of threshold values. The universal threshold, proposed by Donoho and Johnstone, uses the fixed form threshold equal to the square root of two times the logarithm of the length of the signal. LDT, the level dependent threshold, proposed by I.M.Johnstone, and B.W.Silverman, uses a different threshold for each of the levels based on a single formula. Stein Unbiased Risk Estimate (SURE) is an adaptive threshold selection rule. It is data driven and the threshold value minimizes an estimate of the risk. Other threshold values include minimaxi threshold etc. In this paper, an adaptive threshold method is proposed based on EEMD. The threshold values directly relate to the energy of noise on each IMFs. Next, the derivation of adaptive threshold values is initiated by the characteristic of Fractional Gaussian noise (fGn).
fGn is a generalization of white noise. The statistical properties of fGn are controlled by a single parameter H, and the autocorrelation sequence
This can also be defined as:
is the variance of fGn. The value of H is in the range of 0 to 1. The Fourier transform of (18) gives the power spectral density of fGn:
In the decomposing of a given fGn, EMD is worked as a dyadic filter. Restricting to the band-pass IMFs, self-similarity would mean that
Given the self-similar relation (6) for PSDs for band-pass IMFs we can deduce how the variance should evolve as a function of k:
is the variance of the IMF index.
According to lots of simulation:
denotes the H-dependent variation of the IMF energy. In practice, can be estimated from (23), which also gives the model energy of the noisy signal
represents the first IMF coefficients.
According to (21)
According to the relationship between energy and variance are given by
For white noise,
The energies of each IMFs can be defined as:
is the noise energy that can be achieved by the first IMF variance, which can be achieved by (23).
The adaptive threshold value of each IMF can be identified as:
N is the length of signal.
Given these results, a possible strategy for de-noising a signal (with a known H) is generalized as follows:
Decompose the noisy signal into IMFs with EEMD.
Assuming that the first IMF captures most of the noise, estimate the noise level in the noisy signal by computingfrom (28).
Discarding the first IMF, for other IMFs, calculate the adaptive threshold value from (29); shrink the coefficients using the Non-Negative Garrote threshold function.
Reconstruct the signal by the shrunken IMFs, obtain the denoised signal.
6. Results and discussions
To evaluate the performance of the algorithm, Test signal and Clinic noisy ECG signals are processed.
6.1. Test signal
We choose time shifted sine signal which shapes similarly to ECG to test above method; Gaussian White Noise is added as noise, which is zero mean and standard deviation change with the SNR., var means standard deviation. The SNR of noisy test signals are 5. Figure8 shows the original clean signal; figure9 shows the noisy signal; figure10 shows the denoised signal by the above algorithm, the SNR of which achieve 14. Furthermore, the original characteristic waveform is preserved.
6.2. Clinical noisy ECG signal
The ECG signal as Figure.11 illustrates comes from clinical patient. Signal is sampled at 360 Hz; signal length is 1500; the ECG signal is corrupted by noise. Figure12 shows its phase space diagram, which is a plot of the time derivative of the ECG signal against the ECG signal itself. The derivative can accentuate the noisy and high frequency content in signal, so it can better show dramatic improvement after denoising. The noisy ECG signal is processed using the method mentioned above. For the generalized threshold function, m is selected as 2, which is Non-Negative Garrote threshold function. The noisy ECG signal is decomposed into a series of IMFs by EEMD. The first seven IMFs are shown in figure13; the latter seven IMFs are shown in figure 14.The First IMF is discarded owing to predominant noise. Obtain the adaptive threshold value of each IMFs by formula (29). The values are 0.0422, 0.0297, 0.0210, 0.0148, 0.0104, 0.0074, 0.0052, 0.0037, 0.0026, 0.0018, 0.0013, 0.0009, 0.0006. Then shrink the coefficients of each IMFs by the adaptive threshold values and Non-Negative Garrote threshold function. The first shrunken six IMFs are shown in figure15; the latter shrunken seven IMFs are shown in figure16. Reconstruct the signal by the shrunken 13 IMFs and obtain the denoised signal. The filtered ECG signal is illustrated as figure17. The phase space diagram of filtered ECG signal is shown as figure 18. From visual inspection, the ECG signal is much cleaner after being denoised; the original characteristic waveform, especially the sharp Q, R, and S peaks is preserved, without distorting the P and T waves.The results indicate that the method we have proposed significantly reduces noise and well preserves the characteristics of ECG signal.
Another ECG signal as Fig.19 illustrates comes from clinical patient. Signal is sampled at 360 Hz. Length is 1500; the ECG signal is corrupted by noise. Figure20 shows its phase space diagram. The filtered ECG signal is illustrated as figure 21 using the above algorithm. The phase space diagram of filtered ECG signal is shown as figure 22. It is obvious that the noise is reduced.
In this paper, the adaptive noise removal scheme based EEMD is studied for ECG signal. EEMD reduces the mode mixing existing in EMD. Better filtering performance for EEMD is achieved. The adaptive threshold values guarantee the better estimation of noise. We have demonstrated that the algorithm is useful for removing noise from clinic ECG signal. It not only decreases the signal noise, but also the ECG waveform is better conserved. Application of EEMD with adaptive threshold value also has potential for other biomedical signals or other fields.
This work was supported in part by the public welfare Program of Zhejiang Province Science and technology department under the Grant: 2010C31022.