An ECG Signal De-Noising Approach Based on Wavelet Energy and Sub-Band Smoothing Filter

: Electrocardiographic (ECG) signal is essential to diagnose and analyse cardiac disease. However, ECG signals are susceptible to be contaminated with various noises, which affect the application value of ECG signals. In this paper, we propose an ECG signal de-noising method using wavelet energy and a sub-band smoothing filter. Unlike the traditional wavelet threshold de-noising method, which carries out threshold processing for all wavelet coefficients, the wavelet coefficients that require threshold de-noising are selected according to the wavelet energy and other wavelet coefficients remain unchanged in the proposed method. Moreover, The sub-band smoothing filter is adopted to further de-noise the ECG signal and improve the ECG signal quality. The ECG signals of the standard MIT-BIH database are adopted to verify the proposed method using MATLAB software. The performance of the proposed approach is assessed using Signal-To-Noise ratio (SNR), Mean Square Error (MSE) and percent root mean square difference (PRD). The experimental results illustrate that the proposed method can effectively remove noise from the noisy ECG signals in comparison to the existing methods.


Introduction
ECG signals record the electrical activity of the heart, which can reflect its state. Therefore, ECG signals are widely used in heart disease diagnosis [1]. However, an ECG signal is easily degraded by various noises in the acquisition process, which reduces the application value of ECG. There are many reasons for ECG signal degradation, such as electromyographic (EMG) interference, respiratory interference, frequency interference, etc. Due to the different types of ECG signal noise sources, filtering becomes a difficult problem [2]. Therefore, the pre-processing and de-noising of ECG signals is a crucial process before the analysis and diagnosis of ECG .
At present, there are many methods to filter and pre-process ECG signals. Huang et al. [3] introduced a new signal analysis method called empirical mode decomposition (EMD), which is

Discrete Wavelet Transform
Wavelet transform is often used to analyze instantaneous and time-varying signals. Although a classical Fourier transform can reflect the overall connotation of signals, its expression is usually not intuitive enough. As an alternative to Fourier transform, wavelet transform can decompose signals at different scales. Different decomposition scales can be selected for different processing targets. Moreover, the time and the frequency domains can be simultaneously located. For non-stationary signal analysis, it has great advantages.
Discrete wavelet transform is realized by passing the signal through a series of low-pass and high-pass filters [21,22]. Wavelet transform decomposes the signal into detailed components and approximates components of different scales. Figure 1 is a schematic diagram of discrete wavelet transform, where x[n] is the discrete input signal with length n; g[n] is a low-pass filter, which can filter out the high-frequency part of the input signal and output the low-frequency part. h[n] is a high-pass filter. Contrary to the low-pass filter, it filters out the low-frequency part and outputs the high-frequency part. ↓ is a down-sampling filter. At present, many wavelet bases have been developed, such as Haar, Daubechies (Db), Symlet and so on, for the analysis and synthesis of signals. Correct selection of a wavelet basis function plays an important role in de-noising performance [23]. Table 1 shows the comparison of the signal-to-noise ratio (SNR) and the mean square error (MSE) of the de-noised ECG signals which are de-noised by several typical wavelet basis. It can be seen from the table that SNR of db5 is relatively high and MSE is relatively low. Therefore, we select db5 as a wavelet basis in the process of de-noising ECG signal in the paper.

The Proposed Approach
The framework of the proposed approach is shown in Figure 2. The wavelet coefficients that required threshold de-noising are determined by wavelet energy, and then further de-noising through sub-band smoothing filter to obtain the final de-noised smoothed ECG signal. A detailed explanation on every step is provided below.

Noisy ECG
Pre-processing Denoising based on wavelet engery Denoised ECG Signal reconstruction Signal smooth filter Figure 2. the framework of the proposed approach.

Pre-Processing
It is firstly preprocessed by zero averaging before the wavelet decomposition of the signal. The purpose is to remove the DC component in the signal and make the amplitude of the zero-frequency signal very small. In the subsequent power spectrum analysis, the influence of DC energy on the spectrum analysis of low and intermediate frequency signals can be reduced. The zero averaging formula is given as follows: where X is the original input signal, mean(X) represents the mean of X.
The algorithm defaults a tolerance value of 3. Tolerance value is an empericial value, which is adopted to determine the maximum boundary value of the engery change points. For example, if the tolerence value is 3, the position value of the energy change point cannot exceed 3. That is, if the minimum point of energy change occurs at the fourth position, the coefficients of the forward three decompositions are only considered.

Wavelet Decomposition of ECG
Multi-scale and multi-level wavelet decomposition is performed on the signal after the averaging process. As is known to all, the determination of the number of layers of wavelet decomposition affects the subsequent de-noising effect to a certain extent.
If the number of decomposition layers is too small, the de-noising effect will be unsatisfactory. Generally, in order to eliminate high-frequency noise and extract low-frequency components, the number of decomposition layers will be increased to some extent [21], if the number of layers are too high, the error will be considered as large; moreover, it mainly focuses on the reflection of the characteristics of the wavelet basis rather than the signal being analyzed itself, which results in false frequencies and much information loss. For purpose of obtaining more wavelet coefficients, the number of decomposition layers can be appropriately increased, but when the number of decomposition layers is too much, it is easy to cause false frequency.

Calculating the Energy of Each Decomposition Layer
Firstly, the detail coefficients of each layer are calculated according to the detcoef function of Matlab, and then the energy contained in the detail coefficients of each layer is calculated. The calculation formula is as follows: where E j represents the detail coefficient energy of the jth layer, N represents the number of detail coefficient of the jth layer, and co f i represents the ith detail coefficient of the jth layer.

Selecting the Required De-Noising Points
According to the coefficient energy of each layer calculated above, the wavelet energy curve can be obtained, and then the position of the maximum and minimum points of the energy curve can be calculated, that is the decomposition level where they are.
According to the positions of the maximum and minimum points of the calculated wavelet energy curve, the variation points of the number of detail layers that need to be de-noised can be judged. The previous level of the variation points is considered to be de-noised, while the level of the variation point and the subsequent level are not considered to de-noise.
In practical processing, we only locate the initial several detail coefficients that required processing, without considering the low-frequency components behind.Therefore the first position of the energy change point is particularly important. Selecting the location points required de-noising is divided into the following cases: A. When the minimum and maximum values of the energy curve exist simultaneously: 1. When the locationof the first maximum point is less than the location of the first minimum point and the location of the first maximum point is less than the median value of the location of the first maximum point and the location of the first minimum point (If the median is not an integer, round to zero).
When the position of the first maximum point is greater than the tolerance value, the location point is calculated as: When the position of the first maximum point is less than the tolerance value, the location point is calculated as where chg_locs is the location point to be determined and max_locs(1) is the location of the first maximum point. 2. When the location of the first maximum point is greater than the location of the first minimum point, and the position of the first minimum point is less than the median value of the position of the first maximum point and the position of the first minimum point .
When the position of the first minimum point is greater than the tolerance value, the location point is calculated as: When the position of the first minimum point is less than the tolerance value, the location point is calculated as: where min_locs(1) is the location of the first minimum point.
B. When the minimum point of the energy curve exists and the maximum point does not exist.
1. When the position of the first minimum point is greater than the tolerance value, the location point is calculated as: 2. When the position of the first minimum point is less than the tolerance value, the location point is calculated as: C. When the maximum point of the energy curve exists and the minimum point does not exist.
1. when the position of the first maximum point is greater than the tolerance value, the location point is calculated as: 2. When the position of the first maximum point is less than the tolerance value, the location point is calculated as: D. When the minimum and maximum of the energy curve do not exist, we consider that the default noise only appears in the first and second high frequency detail components, and the locating point is set as follows:

Threshold De-Noising
When the location point required de-noising is selected, the detail coefficients of each layer before the location point are de-noised using threshold. Hard threshold and soft threshold are adopted in this paper. In hard threshold processing, when the absolute value of the wavelet coefficients is less than the given threshold, it will be zero; when the wavelet coefficients are larger than the given threshold, it will remain unchanged. That is as follows: where ω is a certain wavelet coefficient and λ is a given threshold.
In soft threshold processing, when the absolute value of the wavelet coefficient is less than the given threshold value, set it as zero; If it is larger than the given threshold value, let the wavelet coefficient subtract the threshold value.That is as follows: Hard threshold is superior to soft threshold in mean square error, however, the signal will produce additional oscillation and jump point, which does not have the smoothness of the original signal. The wavelet coefficients obtained by soft threshold have good continuity and the tested signal will not produce additional oscillation. However, due to the compression of the signal, there will be a certain deviation, which directly affects the approximation level between the reconstructed signal and the real signal. We adopt a fixed threshold to signal de-noising, and the formula is as follows: where λ is the threshold, N is the number of samples, and σ is the noise standard deviation. The noise standard deviation is calculated separately according to each layer, and the values are different in the threshold processing of different layers.

Signal Reconstruction
According to the detailed coefficients of the location point layer and the following layer of the wavelet decomposition, and the detail coefficients from the first layer to the front layer of the locating point after threshold de-noising, and the approximate coefficients, the signal is reconstructed by the wavelet transform, which is the inverse discrete wavelet transform.

signal Smoothing Processing
After threshold de-noising and wavelet reconstruction, some noises still exist in the reconstructed ECG signal. These noises are clearly visible in the low time varying components of the ECG signal i.e., the region between QRS complexes. Therefore, signal smoothing processing is applied for further enhancement of ECG signal quality. In this paper, we determine the subband information required smoothing by the subbands divided by the maximum of the minimum values and the minimum of the maximum values of reconstructed signal. We define lmax as the maximum of the minimum values of reconstructed signal, and upmin as the minimum of the maximum values of reconstructed signal. The reconstructed signal is smoothed point by point from the second point, the processing can be divided into the following situations: A. when the ith point of the signal, i.e., sig(i), is between lmax and upmin, consider its previous point sig(i − 1) and the next point sig(i + 1), and the filtering is as follows: 1. when both sig(i − 1) and sig(i + 1) are between lmax and upmin, the smoothing filter formula is: 2. when sig(i − 1) is between lmax and upmin, and sig(i + 1) is not between lmax and upmin, the smoothing filter formula is: 3. when sig(i + 1) is between lmax and upmin, and sig(i − 1) is not between lmax and upmin, the smoothing filter formula is: 4. when both sig(i − 1) and sig(i + 1) are not between lmax and upmin, the smoothing filter formula is: B. When the ith point of the signal, i.e., sig(i), is not between lmax and upmin, the smoothing filter formula is:

Experimental Environment
To verify the performance of the proposed algorithm, we conduct a number of experiments in our personal computer. The presented approach is implemented in Matlab2012b. We use MIT-BIH arrhythmia database [24] as a reference for real signals. Each sample length is 650,000. The sampling rate is 360 Hz and the resolution is 11 bps. For simulated electrocardiograms, we used the introduced dynamic model to generate synthetic ECG signals [25,26].
In the qualitative evaluation, we select the signal numbered 114 for direct de-noising observation because the record contains more noise. In addition, since there is an ectopic complexes in 119 record, we add noise to record 119 and compared the effect of de-noising with several traditional methods. Finally, In order to observe the de-noising effect of this method under different types of noise, we add power line interference, Gaussian noise and EMG noise to records 106, 113 and 115, respectively.
In the quantitative evaluation, we consider ECG signals such as 100, 103, 105, 113, 115, 117, 119, 122, 200, 215, 230, which contain time-varying QRS morphology,normal and abnormal ECG beats. In each set of experiments, these signals are tested 100 times and averaged. In addition, the record 115 is quantitatively evaluated by adding different types of noise.
We simulated power line interference, EMG noise and synthetic interference. In order to evaluate the performance of the proposed method for power line interference denoising, the additional functions of these interference are as follows: where A = 0.15 mv and f = 50 HZ. The signal was added to the original electrocardiogram signal to simulate power line interference. Random noise is added to ecg signals [27] to simulate EMG interference (µ = 0, σ = 0.15). In the case of synthetic interference, white gaussian noise (WGN) is used to simulate these noises.

Qualitative Evaluation
The performance of the proposed de-noising algorithm is qualitatively analyzed by visual evaluation. Figure 3 shows the results of de-noising the No. 114 record in the MIT-BIH database. As can be seen from the figure, the signal that is de-noised by wavelet obviously reduces a lot of noise, while the signal smoothed further removes more noise.  In Figure 4, we select 119 records from the MIT-BIH arrhythmia database, which contained some ectopic complexs (there exists a complex in the opposite direction of QRS complex before the appearance of QRS complex,called as ectopic complex). We add 20 dB of white Gaussian noise to it to destroy the record. It can be observed from the denoised signals that these ectopic complexes are still preserved in the reconstructed signals, which will not affect the analysis of ECG signals. Compared with the traditional wavelet de-noising, EMD and EMD soft threshold, the ECG signal de-noised by the proposed method is closer to the original ECG signal and smoother.    Figure 5 shows the de-noising results of signal after adding 50 HZ power line interference to No. 106 record of the MIT-BIH database. Figure 6 displays the results of adding 20 dB WGN to No. 113 record of the MIT-BIH database. Figure 7 shows the de-noising results after adding EMG interference to No. 115 record of the MIT-BIH database. As can be seen from these figures, the de-noising effect is obvious and the distortion is smaller than that of the original signal. The results show that the proposed method can effectively de-noise all kinds of noises, and the de-noised signal retains the characteristics of the original signal. Moreover the signal is very similar to the original signal.
The evaluation also extends to the synthetic ECG signal. We add 15 dB Gauss white noise to the synthetic signal and denoise it. As shown in Figure 8, it can be seen that the algorithm also has excellent denoising effect on the synthetic ECG signal.

Quantitative Comparison
In this paper, the performance of de-noising signal is analyzed by SNR, MSE and PRD, and the quality of reconstructed signal is studied. The SNR can be divided into output SNR and improved SNR.
The output SNR formula is defined as follows: The formula of improved SNR is defined as follows: The MSE is used to measure the deviation between the original signal and the de-noised signal. The formula is given as follows: The percent root mean square difference is defined as follows: where x (i) is the original signals, y (i) is the signal after adding noise, s (i) is the final de-noised and smoothed signal, and L is the signal length.
The experimental comparison shows that the de-noising effect is better when the wavelet basis of DB family is used for wavelet transform. In order to compare with the traditional wavelet-based ECG de-noising scheme, we tested the traditional wavelet-based ECG de-noising scheme (WT), multi-adaptive bionic wavelet transform (MABWT) and our proposed method in MIT-BIH arrhythmia database [28]. Table 2 (Soft Threshold) and Table 3 (Hard Threshold) respectively show the improved SNR results of ECG signals de-noising with WGN (5 dB) added. It can be seen that compared with the traditional wavelet-based de-noising scheme and multi-adaptive bionic wavelet transform, the method proposed in this paper has a higher SNR improvement whether it uses soft threshold or hard threshold. Table 2. SNR imp comparison of de-noising results using soft threshold for ECG signals with WGN (5 dB). Table 3. SNR imp comparison of de-noising results using hard threshold for ECG signals with WGN (5 dB). Tables 4-6 shows the results of quantitative comparison with several other methods after adding different noise to the 115th record for denoising, namely,the parallel-type fractional zero-phase filtering (FZP), zero-phase Butterworth filter (BZP) [29], the Riemann Liouville (RL) integrator [30] and the proposed method. We select record 115 because the record is relatively clean and suitable for observing the de-noising effect of the algorithm on different types of noise. Table 4 gives a comparison of power line interference (50 Hz), with an amplitude of 0.15 mv. In Table 5, EMG interference is compared and these noises are simulated using random noise. Table 6 gives the comparison results in the case of synthetic interference. These noises are simulated using WGN with an input signal-to-noise ratio of 15 dB. Tables 4-6 shows that the proposed algorithm provides better de-noising results and has smaller deviation compared with the original ECG signal. The ability of the algorithm in tracking signal and removing noise is proved.  Figure 9 illustrates the comparison between the proposed de-noising method and the traditional wavelet de-noising method (soft threshold and hard threshold) for a wide range of input SNR [28]. It can be observed that for low SNR, the proposed method has a high improvement SNR. With the increase of the input SNR, the improved SNR is still high and stable. Unlike traditional wavelet denoising, the improved SNR is very low or even negative for high input SNR.  Table 7 and Figure 10 show the comparison of the proposed method with EMD soft threshold [4], wavelet soft threshold [12] and EMD-wavelet [13]. Table 7 shows the results of the MSE for all the comparison methods using different ECG signals under consideration at a particular input SNR level of 20 dB. It is vivid from this table that the proposed method yields the smallest MSE and verify its capability to yield enhanced ECG signal with better quality. A bar diagram presenting the percentage PRD(%) results obtained by using different de-noising methods are plotted in Figure 10. We add Gauss white noise with input signal-to-noise ratio of 20 dB to ECG signal. The histogram clearly shows that the proposed method performs best for any ECG signal because it yields the lowest PRD(%).

Conclusions
Wavelet transform is a signal analysis tool. The tight support and orthogonality of wavelet function cause it to have good performance in signal de-noising. A new de-noising method of ECG signals based on wavelet energy and sub-band smoothing is proposed in this paper. The wavelet coefficients that required threshold de-noising are determined by wavelet energy, and then the sub-band smoothing filter is used for further de-noising. The experimental results show that the proposed method can significantly improve the quality of ECG signal, and the de-noising results have higher SNR, lower MSE and lower PRD. In the future, we will attempt to adopt the image/video processing methods [31][32][33][34][35][36][37][38] to process ECG signal. With the wide application of deep learning in image processing, We will adopt network optimization methods [39][40][41][42][43][44][45][46][47][48] to improve the real-time and high efficiency performance of ECG signal processing.