Noise Reduction in ECG Signal Using an Effective Hybrid Scheme

Electrocardiogram (ECG) is a critical biological signal, which usually carries a great deal of essential information about patients. The high quality ECG signals are always required for a proper diagnosis of cardiac disorders. However, the raw ECG signals are highly noisy in nature. In the paper, we propose a hybrid denoising scheme to enhance ECG signals by combining high-order synchrosqueezing transform (FSSTH) with non-local means (NLM). With this method, a noisy ECG signal is first decomposed into an ensemble of intrinsic mode functions (IMFs) by FSSTH. Then, some noise is removed by eliminating a set of noisy IMFs that are determined by a scaling exponent obtained by the detrended fluctuation analysis (DFA); while the remaining IMFs are filtered by NLM. Finally, the denoised ECG signal is obtained by reconstructing the processed IMFs. Experiments are carried out using the simulated ECG signals and real ones from the MIT-BIH database, and the denoising performances are evaluated in terms of signal to noise ratio (SNR), root mean squared error (RMSE) and percent root mean square difference (PRD). Results show that the hybrid denoising scheme involving both FSSTH and NLM is able to suppress complex noise from ECG signals more effectively while preserving the details well.


I. INTRODUCTION
The electrocardiogram (ECG) signals, which have been widely used to diagnose early-stage cardiovascular disorders, could reflect the electrical activity of heart and provide the clinical information about heart [1]- [4]. However, the ECG signals are highly noisy due to the inevitably noise introduced during signal collection and transmission. It is known that a noise-free ECG signal is much desired for a reliable and efficient down-stream analysis. Therefore, the noise reduction of ECG signal via appropriate algorithms is of great importance in the field of biomedical engineering [5], [6].
In the past, a number of methods have been proposed for ECG signal denoising, which can generally be divided into four categories. The first category is the spatiotemporal domain method, e. g. classical wiener filters [7], median filter [8] and adaptive filtering [9]- [11]. However, such methods are usually inaccurate for a noisy ECG due to the non-stationary nature of the cardiac signal. The second The associate editor coordinating the review of this manuscript and approving it for publication was Nadeem Iqbal. category is the frequency domain algorithm, for example, bandpass filtering [12] and wavelet-based denoising technique [13], [14]. The former is often incapable of suppressing noise overlapping with cardiac components in the frequency domain. Though the latter achieves noise removal by shrinking the wavelet coefficients of ECG signal and performing inverse wavelet transform, the selection of basis function and threshold is a troublesome process. Meanwhile, the hard thresholding algorithm may cause oscillation in the denoised ECG signal while the soft thresholding algorithm may reduce the amplitudes of R wave [15]. The third one is the statistical method such as independent component analysis (ICA) [16] and principle component analysis (PCA) [17], which are powerful in removing in-band noise by eliminating the dimensions associated with noise. However, the mapping model is very sensitive to minor disturbances in either signal or noise. The last one is the mode decomposition based method, for instance, empirical mode decomposition (EMD) [18], [19], ensemble empirical mode decomposition (EEMD) [20], [21], complete ensemble empirical mode decomposition (CEEMD) [22], extreme-point symmetric mode decomposition (ESMD) [23] and variational mode decomposition (VMD) [24], [25]. The major disadvantage of these methods is mode-mixing. Although the variations of EMD have greatly alleviated this phenomenon, the redundant information still exists in the decomposition result, which limits the further application of such methods in many aspects.
More recently, Pham and Meignen (2017) [26] developed a new adaptive signal analysis algorithm, which is termed as high-order synchrosqueezing transform (FSSTH). It is a new generalization of the STFT-based synchrosqueezing transform by computing more accurate estimates of the instantaneous frequencies using higher order approximations for both the amplitude and phase, which results in perfect concentration and reconstruction for a wider variety of signals [27]. Currently, FSSTH has been successfully applied to the analysis of a transient gravitational-wave signal [26], seismic time-frequency analysis [28] and machine fault diagnosis [29]- [31]. However, the applications in the field of biomedical engineering are seldom reported.
In this paper, inspired by the advantage of FSSTH, we propose a hybrid denoising scheme for the ECG signal based on high-order synchrosqueezing transform and non-local means (NLM). The proposed method starts by decomposing a noisy ECG signal into a series of IMFs using the FSSTH. Then, considering the scaling exponent extracted from the detrended fluctuation analysis (DFA) for each IMF, the obtained IMFs are separated into two parts: (a) signal-dominant components and (b) noise-dominant components. The former are filtered by the NLM while the latter are eliminated. Finally, the denoised ECG signal can be recovered by reconstructing the filtered IMFs. Both the simulated signal and experimental ECG signals from the MIT-BIH database have verified the performance of the proposed method on noise reduction for ECG signals. Our contributions are as follows: (1) we first bring the FSSTH to process ECG signal, (2) an empirical equation is introduced to adaptively determine the number of modes from FSSTH, and (3) a scaling exponent obtained by DFA is utilized as a new threshold to distinguish between signal and noise.
The remainder of this paper is organized as follows. In section II, we describe the general principle of high-order synchrosqueezing transform, non-local means, detrended fluctuation analysis, and the proposed hybrid denoising algorithm for the ECG signal. The experimental results using the simulated signal and real-life ECG signals from the MIT-BIH database are shown in section III, where the traditional thresholding algorithm and non-local means algorithm are compared to the proposed method. A discussion of the results is given in section IV and the concluding remarks are contained in section V.

II. THE PROPOSED METHOD A. HIGH-ORDER SYNCHROSQUEEZING TRANSFORM
The high-order synchrosqueezing transform (FSSTH) is a new extension of the conventional STFT-based SST (FSST), which was first proposed by Thakur and Wu [32]. It achieves more accurate estimates of the instantaneous frequencies by using higher order approximations for both the amplitude and phase [26], thus resulting in a perfect concentration and reconstruction for a wider range of AM-FM signal modes.
An AM-FM signal is represented as: where A (t) and φ (t) are the instantaneous amplitude and phase, respectively.
The STFT of signal f can be defined as: where g * is the complex conjugate of a window function g.
The traditional STFT-based SST (FSST) is expressed as follows: where γ is the threshold and δ is the Dirac distribution. ω f (t, η) denotes the instantaneous frequency estimate at time t and frequency η, and is estimated as: where R {Z } represents the real part of complex number Z , and ∂ t is the partial differentiation with respect to t. Despite a rigorous mathematical framework for FSST, its application is restricted to a class of signal with 'slowly varying' instantaneous frequency, which is due to the requirement of the weak frequency modulation hypothesis for the modes comprising the signal.
The new instantaneous frequency estimate from the high-order SST is based on the high order Taylor expansions of the amplitude and phase of a signal. The Taylor expansion of the signal in equation (1) for τ close to t is given as: where Z (k) (t) is the kth derivative of Z at time t. Therefore, the STFT of this signal at time t and frequency η can be rewritten as: The local instantaneous frequency estimate ω f (t, η) can be written as: Now, introducing the frequency modulation operator q [k,N ] η,f : The N th-order local complex instantaneous frequency, ω [N ] η,f at time t and frequency η, can be represented by: Finally, the high-order SST is defined by replacing η,f (t, η) in equation (3), yielding: And, the mode is approximately recovered by: where d denotes the compensation factor and ϕ (t) is an estimate for φ (t).
A synthetic signal ( Figure 1) characterized by the nonlinear cosinoidal frequency, used previously in [33], is employed to evaluate the performance of the high-order SST. Figures 2(a) and (b) show the time-frequency representations based on FSST and FSSTH, respectively. It is clear that FSST leads to a blurred time-frequency map, in which the low time-frequency resolution will greatly decrease the precision of instantaneous frequency estimation. On the other hand, a sharpened time-frequency representation is provided by FSSTH, thus facilitating the instantaneous frequency feature interpretation. Therefore, FSSTH is more suitable for coping with fast frequency-varying signal compared to FSST.

B. NON-LOCAL MEANS
The non-local means (NLM) [34], [35] is one of the most popular denoising algorithms, which can achieve excellent detail preservation. For a given sample i, the estimate u (i) is where Z (i) = j ω (i, j), and the weights are: where λ is a bandwidth parameter; denotes the neighborhood surrounding i and containing L samples; d 2 (·, ·) represents the squared Euclidean distance between a set of selected samples centered on i and j considering the similarity of their neighborhoods.
It should be noted that there are three key parameters in the NLM algorithm, namely, the bandwidth λ; the patch size, specified as a half-width P (L = 2P + 1); and the size of N (i), specified as a half-width Q. The parameter λ determines the amount of smoothing applied. A smaller λ will lead to noise fluctuations while a larger λ will cause dissimilar patches to appear similar. A good overall selection of λ is 0.5σ [36], where σ is the noise standard deviation. The patch half-width P chooses the scale on which patches are compared, a reasonable selection for P is the half-width of the high-amplitude 'R' ECG complex [35]. In principle, increasing the half-width Q will achieve better performance; meanwhile, the computation burden should be taken into account.

C. DETRENDED FLUCTUATION ANALYSIS
Detrended fluctuation analysis (DFA), proposed in [37], is an appropriate method to obtain more reliable scaling exponent when the signal has non-stationary properties. The key idea of the DFA is to determine how the root mean square (RMS) fluctuation of the time series around the local trend in the box size changes as a function of the time scale n [37].
The first step of the algorithm is to find the integrated time series y (k) after removing the average, which is described as follows: where x is the average of the time series within the range [1, N ], and N is the number of samples. Then, y (k) is divided into n boxes. For each box, the local trend y n (k) is estimated using the residue decomposed by EMD. The RMS fluctuation F (n) is defined as: After n iterations, the slope of the plot of log [F (n)] log (n), α, is called the scaling exponent, and can be formulated as: It should be mentioned that the scaling exponent α provides a quantitative score and can be used to characterize the temporal correlation in the time series. Generally speaking, the value α = 0.5 is corresponding to the completely uncorrelated data. If 0 < α < 0.5, the signal is anticorrelated. When 0.5 < α < 1, the long-range temporal correlations exist [38]- [40]. If α > 1, the correlations do not exhibit power-law behavior with respect to time. In addition, the scaling exponent α also indicates the roughness of the time series. The larger the α, the smoother the time series is, and vice versa. DFA has become a widely used method for the determination of scaling properties and detection of long-range correlations in non-stationary time series owing to its ability of detrending the time series [41]. In the paper, we utilize DFA to determine the number of IMFs and which IMF is the signal or the noise. Figure 3 illustrates the proposed denoising workflow. The specific implementation procedure of the hybrid denoising algorithm for ECG signal based on FSSTH and NLM is summarized as follows:

D. PROCEDURE OF THE HYBRID DENOISING FOR ECG SIGNAL
1) Collect the ECG signal.
2) Decompose the noisy ECG signal using the FSSTH and obtain an ensemble of band-limited IMFs.
3) Estimate the scaling exponent with respect to each IMF by the DFA to determine the number of IMFs from FSSTH. 4) Evaluate each IMF via the DFA and determine the threshold. 5) Eliminate the noisy IMFs based on thresholding criteria. 6) Denoising the rest of IMFs by means of NLM. 7) Obtain the denoised ECG signal by reconstructing the processed IMFs.

III. EXPERIMENTS AND RESULTS
Several metrics, such as signal to noise ratio (SNR), root mean squared error (RMSE) and percent root mean square difference (PRD), are used as quantitative performance estimators in this study.
The SNR indicates the noise reduction level. The greater SNR is, the better denoising performance can be obtained. The SNR is described as: The RMSE is employed to determine the variance between the desired output and the actual one. A lower RMSE corresponds to a smaller difference. The RMSE is defined as: The PRD implies the recovery performance by measuring the error between the input signal and the reconstructed signal. A lower PRD represents a better reconstruction. The PRD is obtained by the following formula: where s i (t) is the original ECG signal, ∧ s i (t) is the denoised ECG signal, and N is the length of the ECG signal.

A. EVALUATION ON SIMULATED ECG SIGNAL
For the simulated ECG signal, Gaussian white noise with SNR of 0, 4, 8, 12 and 16dB, baseline wander (BW), muscle artifact (MA), and electrode motion (EM) that are with SNR of 4dB will be added to the noise-free ECG signal. It is worth noting that the clean ECG signal and four kinds of noises are simulated using Open-Source Electrophysiological Toolbox (OSET) [23]. Figure 4 shows the clean ECG signal (Figure 4(a)) and noisy one with 4 dB SNR (Figure 4(b)) that is corrupted by the aforementioned noise. Now, FSSTH is employed to decompose the noisy ECG signal. To make the FSSTH adaptive, we need to determine the optimal number of IMFs (K ) in a data-driven way. In the paper, we utilize an empirical equation for the determination of K [39]: where α is the scaling exponent from DFA of the input signal.  Figure 5 shows the relationship between scaling exponents and the number of IMFs from FSSTH, whose values are listed in Table 1. The first IMFs with larger scaling exponent are mostly comprised of the signal while the last IMFs with lower scaling exponent mainly carry noise. As shown in Figure 5, if the K is lower, for example 2, 3, 4, the IMFs are signal-dominant, and the scaling exponents are larger than 0.1 ( Figure 5(a)). Whereas, when K is large, for example 6, 7, 8, the scaling exponents do not always decrease with the IMF number increasing. It is more obvious when K = 8 ( Figure 5(c)), which indicates that the IMFs are noise. Figure 5(b) displays the scaling exponent when K equals to 5. The IMF includes both signal and noise, and the α value of the last IMF is lower than 0.1. Based on the observations, we select 5 for K as an appropriate number of IMFs from FSSTH. Figure 6 shows the IMFs extracted by using FSSTH, where the decomposition gives 5 individual modes. The first two IMFs with larger scaling exponent are signal-dominant while the remaining ones are mainly noise (Figure 5(b)). Thus, α = 0.2 is selected as a threshold in this example. In other words, IMF3, IMF4 and IMF5 are removed; IMF1 and IMF2 are filtered by NLM, which are shown in Figure 7. As can be obviously seen, most of the noise has been successfully suppressed and the effective waveforms are outstood compared to the noisy IMFs, especially for IMF2.
Then, the filtered IMFs with higher scaling exponent than the threshold are used for reconstruction in order to obtain the denoised ECG signal. Figures 8(a) and (b) depict the denoised results using the traditional thresholding method that directly neglects a number of IMFs containing noise and NLM, respectively. Figure 8(c) shows the denoised ECG signal using the proposed hybrid denoising scheme. It can be clearly seen that the proposed method exhibits the better performance in noise removal. The waveform components such as QRS   complex, P wave and T wave have become clearer compared to the noisy ECG signal. Meanwhile, PQ and ST segments are also smoother, which is very helpful for capturing the feature VOLUME 8, 2020 of ECG signal. By contrast, some noise still remains in the denoised ECG signals for the other two methods, especially for the traditional thresholding approach. In addition, for better visualization, we also enlarge the comparative areas corresponding to Figure 8 (the blue rectangles), which are plotted in Figure 9. One clearly sees that our proposed method is superior to the traditional thresholding and NLM, because the former is closer to the noise free ECG signal (Figure 4(a)). Figure 10 shows the comparisons based on SNR, RMSE and PRD for all methods implemented on the simulated ECG signal corrupted with several noises of different SNR levels. We can observe that the proposed method outperforms the other two algorithms by achieving much higher SNR and much lower RMSE and PRD. Thus, one can conclude that the proposed method performs better in both noise removal and signal preservation.

B. EVALUATION ON REAL ECG SIGNAL
To further assess the validity of our method, we apply the proposed hybrid denoising algorithm to a series of real ECG signals. Here, we have taken into account the 100m.dat and 111m.dat from the MIT-BIH Arrhythmia Database and 12713_04m.dat and 13420_12m.dat from the MIT-BIH ECG Compression Test Database as examples. The sampling frequencies of the two datasets are 360Hz and 250Hz, respectively. It is worth noting that both databases have already been widely used in the related researches [22], [23], [42]- [44]. Figures 11, 13, 15 and 17 show the denoised results for real ECG signals from the two datasets by using the traditional thresholding algorithm (11(b), 13 (b), 15(b) and 17(b)), NLM (11(c), 13 (c), 15(c) and 17(c)) and the proposed method (11(d), 13 (d), 15(d) and 17(d)), respectively. The enlarged denoised results corresponding to Figures 11, 13, 15 and 17 are illustrated in Figures 12, 14, 16 and 18, respectively. It can be observed that the traditional thresholding method presents the poor performance in noise removal, because some amount of noise still resides in the denoised ECG signals, which may have a great effect on the clinical diagnosis of heart disease. In contrast, NLM and the proposed method do a better job in noise reduction. The main difference between them lies in the relative flat segments in the ECG signals. Compared to NLM, the proposed method performs better with less oscillations and smoother trends (Figures 12, 14, 16 and 18), which is of particular importance in characterizing the PQ and ST segments correctly in the ECG signals.

IV. DISCUSSIONS
The electrocardiography (ECG) signal plays a key role in diagnosing diverse kinds of heart diseases. The noise may  contaminate the ECG signal during its acquisition and transmission. Thus, denoising ECG signal is required to avoid wrong interpretation in ECG.  In this paper, the traditional thresholding algorithm and NLM are employed to compare with the proposed hybrid denoising method. The traditional thresholding VOLUME 8, 2020  algorithm discards some IMFs containing the noise directly after the FSSTH decomposition. Though some strong and high-frequency noise is attenuated, the weak and  low-frequency noise fails to suppress effectively. Therefore, the denoised ECG signals are characterized by the lower SNR due to the presence of residual noise. The NLM, originally  proposed for image denoising [34], can achieve noise suppression by replacing a noisy pixel with the weighted average of all the pixels in a pre-defined search window where the weights are robustly determined by the nonlocal comparison of image patches [23], [35]. However, the NLM is unable to deliver sufficient noise reduction, especially for the ECG signal corrupted highly by noise because of the inaccurate weight computation based on noisy signals [23]. In the proposed method, most of the noise is first removed by means of a new thresholding determined by a scaling exponent obtained by DFA. The remaining noise is then filtered by NLM. Thus, it can achieve adequate noise reduction and excellent signal preservation. In addition, we also evaluate the computation costs of the above methods operating on the 100m.dat and 13420_12m.dat, which is listed in Table 2. It can be seen that the proposed approach takes a longer time than the traditional thresholding algorithm and NLM. The reason lies in the fact that much computational time is involved in the ECG signal decomposition based on FSSTH and the NLM denoising for different IMFs. It should be mentioned that all experiments are done on a PC station equipped with an Intel(R) Core(TM) i5-7500 CPU clocked at 3.4 GHz and 8 GB of RAM. Additionally, we also listed the abbreviations of many technical terms used in the paper and their full names, which are shown in Table 3.
The purpose of this paper is removing the noise in ECG signal using a hybrid denoising scheme based on the high-order synchrosqueezing transform (FSSTH) and non-local means (NLM). The higher time-frequency resolution feature of FSSTH in signal decomposition is combined with the excellent detail preservation performance of NLM in signal denoising. An important finding in applying our method is the considerable reduction of noise with the details preservation for ECG signal, which contributes to more accurate diagnosis of the heart diseases. On the other hand, it is worth noting that the hybrid method gains more attention for the denoising of ECG signal, because it combines the advantage of both algorithms and achieves the excellent denoising performance. We believe that it is the potential avenues for ECG signal denoising method.

V. CONCLUSION
We have proposed a novel hybrid scheme for ECG signals denoising. In the paper, we use the FSSTH to decompose the noisy ECG signal into a series of IMFs, in which the ECG signal and the noise have the different property. Then, an empirical equation based on DFA is introduced to adaptively determine the number of IMFs and a scaling exponent obtained by DFA is utilized as a threshold to distinguish the signal-dominant components with the noise-dominant ones. Finally, the denoised ECG signal can be recovered through the reconstruction of the signal-dominant IMFs filtered by the NLM. We have also applied the proposed method on simulated and real ECG signals and compared the results with the results of the traditional thresholding and NLM methods. We observe that the proposed method can obtain better result in both noise removal and signal preservation with the higher SNR as well as the lower RMSE and PRD. Future work will focus on improving the computational efficiency of the proposed method so that it can be better applied in practice. ZHIHUA ZHANG received the B.E. degree in safety science and engineering from the North University of China, in 2018. He is currently pursuing the M.E. degree with the Beijing University of Chemical Technology. His research interests include signal analysis, deep learning, and intelligent fault diagnosis. VOLUME 8, 2020