Noise Suppression in ECG Signals through Efficient One-Step Wavelet Processing Techniques

Abstract

This paper illustrates the application of the discrete wavelet transform (DWT) for wandering and noise suppression in electrocardiographic (ECG) signals. A novel one-step implementation is presented, which allows improving the overall denoising process. In addition an exhaustive study is carried out, defining threshold limits and thresholding rules for optimal wavelet denoising using this presented technique. The system has been tested using synthetic ECG signals, which allow accurately measuring the effect of the proposed processing. Moreover, results from real abdominal ECG signals acquired from pregnant women are presented in order to validate the presented approach.

1. Introduction

Electrocardiogram (ECG) is a valuable technique that has been in use for over a century, not only for clinic applications [1]. ECG acquisition from the human skin involves the use of high gain instrumentation amplifiers. This fact makes the ECG signal to be contaminated by different sources of noise [2]. This circumstance is highlighted when the target is the measurement of fetal ECG signals acquired over the mother’s abdomen [3]. For processing ECG signals, it is necessary to remove contaminants from these signals that make visual inspection and ECG feature extraction difficult. In general, ECG contaminants can be classified into different categories, including power line interference, electrode pop or contact noise, patient-electrode motion artifacts, electromyographic (EMG) noise, and baseline wandering. Among these noises, the power line interference and the baseline wandering (BW) are the most significant and can strongly affect ECG signal analysis. Except for these two noises, other noises may be wideband and usually involve a complex stochastic process, which also distorts the ECG signal. The power line interference is narrow-band noise centered at 50 Hz or 60 Hz with a bandwidth of less than 1 Hz. Usually the ECG signal acquisition analog hardware can remove the power line interference. However, the baseline wandering and other wideband noises are not easy to be suppressed by analog circuits. Instead, the software scheme is more powerful and feasible for portable ECG signal processing. Thus, denoising this type of signals is decisive for further parameter extraction in clinic applications.

Wavelet transform (WT) [4] is a useful tool for a variety of signal processing [5, 6] and compression applications [7, 8]; its primary, and most advantageous, application areas are those that have to generate or process wideband signals. This transform produces a time-frequency decomposition of the signal under analysis, which separates individual signal components more effectively than the traditional Fourier analysis. This fact makes WT one of the most used tools for biosignal processing, with ECG being an obvious candidate for this type of analysis. Discrete wavelet transform (DWT) provides a multiresolution, analysis, which allows representing a signal by a finite sum of components at different resolutions so that each component can be processed adaptively based on the objectives of the application.

This paper proposes an arrangement of discrete wavelet transform structures for ECG signal processing on portable, embedded computing real-time implementations [9], focusing on the suppression of different types of noise including DC levels and wandering. This suppression is carried out with a one-step wavelet processing, which reduces computing cost. Appropriate software models and parameter selection are presented based on mathematical analysis for the mentioned ECG digital processing. Moreover, these software models will also enable very quick tests of the required signal processing. In order to select an appropriate method for the development of the software modeling, noise suppression and feature extraction procedures on ECG signals are studied in the following.

2. Materials and Methods

2.1. DWT Fundamentals

Wavelet transform (WT) is one of the most used tools in multiresolution signal analysis due to its ability to decompose a signal at various resolutions, which allow observing high-frequency events of short duration in nonstationary signals [10]. The continuous WT of a signal 𝑓(𝑡)∈𝐿²(𝑅) is defined as [11] 𝑊(𝑎,𝑏)=^+∞_−∞𝑓(𝑡)𝜓_𝑎,𝑏(𝑡)𝜓𝑑𝑡,_𝑎,𝑏(𝑡)=1√𝑎𝜓^∗𝑡−𝑏𝑎. (1)(2) In (2), * denotes complex conjugate, and 𝑎 is a scale factor and 𝑏 a translation factor. The normalization factor 𝑎^−1/2 is included so that ‖𝜓_𝑎,𝑏‖=‖𝜓‖. Thus, 𝜓_𝑎,𝑏(𝑡) represents a shifted and scaled version of the so-called mother wavelet 𝜓, which is a window function that defines the basis for the wavelet transformation. A mother wavelet 𝜓 of order 𝑚 is a function 𝜓∶𝑅→𝑅 which satisfies the following four properties.(1)If 𝑚>1, then 𝜓 is (𝑚−1) times differentiable.(2)𝜓∈𝐿^∞(𝑅). If 𝑚>1, for each 𝑗∈{1,…,𝑚−1},𝜓^(𝑗)∈𝐿^∞(𝑅). (3)𝜓 and all its derivatives up to order (𝑚−1) decay rapidly: for each 𝑟>0 there is a 𝛾>0 such that 𝜓^(𝑗)(𝑡<1)𝑡,𝑗∈{1,…,𝑚−1}foreach|𝑡|>𝛾. (3)(4)Foreach𝑗∈{0,1,…,𝑚}, ∫𝑡^𝑗𝜓(𝑡)𝑑𝑡=0.Practical applications rely on the discrete wavelet transform (DWT), since it provides enough information while requiring reasonable computation time and resources. The DWT of a discrete function 𝑓(𝑛) is defined as 𝑊(𝑎,𝑏)=𝐶𝑗,𝑘=_𝑛𝑓(𝑛)𝜓_𝑗,𝑘(𝑛), (4) where 𝜓_𝑗,𝑘(𝑛)=2^−𝑗/2𝜓2^−𝑗𝑛−𝑘 (5) and 𝑎=2^𝑗 and 𝑏=𝑘2^𝑗, with 𝑗,𝑘∈ℤ. Thus, DWT is computed in practice through a set of two FIR-like filters, lowpass and highpass, at each decomposition level 𝑖, followed by a downsampling by 2, which implies a reduction in the sampling frequency. The result of the low-pass filter is usually known as the approximation and is used as input to the following decomposition level, while the result of the high-pass filter is called the detail. Thus, the usual DWT corresponds to the scheme in Figure 1(a), with 𝑎_𝑛^(𝑖)=^𝑁−1_𝑘=0𝑔_𝑘𝑎^(𝑖−1)_2𝑛−𝑘𝑑𝑖=1,2,…,𝐽,_𝑛(𝑖)=^𝑁−1_𝑘=0ℎ_𝑘𝑎^(𝑖−1)_2𝑛−𝑘𝑎_𝑛⁽⁰⁾≡𝑥_𝑛 (6) for an 𝑁-sample input sequence. The coefficients of the low-pass and high-pass filters are defined by the family of wavelet functions used as basis for the transformation: Haar, Daubechies, Quadratic Spline, and so forth [11].

(a)

(b)

Figure 1 

Block diagram of the multiresolution decomposition using the wavelet transform, (a) analysis filter bank and (b) reconstruction filter bank.

Since every detail is the result of high-pass filtering of the previous approximation and approximations are the result of low-pass filtering, these details and approximations at every decomposition level contain information at different frequencies and time scales; this reflects the multiresolution analysis of the WT. Moreover, it is possible to reconstruct the original signal from the set of details and approximations, through the inverse DWT, so 𝑎_𝑚^(𝑖−1)=^(𝑁/2)−1_𝑘∑₌₀𝑔_2𝑘𝑎₍₎^(𝑖)_{𝑚/2−𝑘}+∑^(𝑁/2)−1_𝑘=0ℎ_2𝑘𝑑₍₎^(𝑖)_{𝑚/2−𝑘}𝑚even_{(𝑁/2)−1𝑘}∑₌₀𝑔_2𝑘+1𝑎₎₎^(𝑖)_{((𝑚−1/2−𝑘}+∑^(𝑁/2)−1_𝑘=0ℎ_2𝑘+1𝑑₎₎^(𝑖)_{((𝑚−1/2−𝑘}𝑚odd, (7) where, once more, the corresponding details are high-pass filtered and the approximations are low-pass filtered, this time through the appropriate reconstruction filters. Thus, the reconstructed approximation at level 𝑖−1 is obtained by addition of the output of the reconstruction filters and an upsampling by 2, as shown in Figure 1(b).

2.2. DWT-Based ECG Processing

The special characteristics of ECG signals and their frequency spectrum make the use of traditional Fourier analysis for the detection of ECG features difficult [12]. Wavelet transform can be applied in many fields, but its primary, and most advantageous, application areas are those that have, generate, or process wideband signals. This is due to the multiresolution analysis that wavelet provides, which allows representing a signal by a finite sum of components at different resolutions, so that each component can be processed adaptively based on the objectives of the application. In this way, this technique represents signals compactly and in several levels of resolution, which is ideal for decomposition and reconstruction purposes. Thus, the discrete wavelet transform is an effective way to digitally remove noises within specific subbands for ECG signals.

Several wavelet families have been proposed for ECG processing [12]. Selecting the right wavelet for this application is a task requiring at least a brief discussion. The wavelet type to use for the discrete wavelet analysis is an important decision for this processing. Singh and Tiwari [13] studied different wavelet families and analyzed the associated properties. We can distinguish between two types: the orthogonal wavelets as Haar, Daubechies (db𝑥𝑥), Coiflets (coif𝑥), Symlets (sym𝑥), and biorthogonal (bior𝑥_𝑥), where 𝑥 indicates the order of the wavelet and the higher the order, the smoother the wavelet. The orthogonal wavelets are not redundant and are suitable for signal or image denoising and compression. A biorthogonal wavelet is a wavelet where the associated wavelet transform is invertible but not necessarily orthogonal. Designing biorthogonal wavelets allows more degrees of freedom than orthogonal wavelets. One of these additional degrees of freedom is the possibility to construct symmetric wavelet functions. The biorthogonal wavelets usually have the linear phase property and are suitable for signal or image feature extraction.

2.2.1. Wavelet-Based Wandering Suppression

Baseline wandering (BW) usually comes from respiration at frequencies varying between 0.15 and 0.3 Hz. Different methods for wandering suppression exist, will the common objective of them to make the resulting ECG signals contain as little baseline wandering information as possible, while retaining the main characteristics of the original ECG signal. One of the proposed methods consists in high-pass digital filtering; for example, a Kaiser Window FIR high-pass filter [14] could be designed, where appropriate specifications of the high-pass filter should be selected to remove the baseline wandering.

In addition to digital filters, the wavelet transform can also be used to remove the low frequency trend of a signal. This wavelet-based approach is better because it introduces no latency and less distortion than the digital filter-based approach. The required steps for the application of this wavelet-based processing for baseline wander correction are the following.(1)Decomposition: apply wavelet transform to the signal up to a certain level 𝐿 in order to produce the wavelet approximation coefficients 𝑎_𝑛^(𝐿) that captures the baseline wander.(2)Zeroing approximation coefficients: replace the approximation vector 𝑎_𝑛^(𝐿) for an all-zero vector to subtract this part from the raw ECG signal and remove the wandering.(3)Reconstruction: compute wavelet reconstruction, based on these zeroing approximation coefficients of level 𝐿 and the detail coefficients of levels from 1 to 𝐿, in order to obtain the BW corrected signal.The main idea is to remove the low frequency components, which better estimate the baseline wander. This processing is easy to carry out using wavelet decomposition, for which it is necessary to select the proper wavelet function and resolution level.

2.2.2. Wavelet-Based Denoising

The goal in ECG denoising is to try to recover the clean ECG from the undesired artifacts with minimum distortion. The recovered signal is called denoised signal, and it allows further ECG processing, such as in the case of separation of fetal ECG [3], QRS complex detection, and parameter estimation (such as the cardiac rhythm). The underlying model for the noisy signal is basically of the following form: 𝑠(𝑛)=𝑓(𝑛)+𝜎𝑒(𝑛), (8) where 𝑠(𝑛) represents the noisy signal, 𝑓(𝑛) is an unknown, deterministic signal, time 𝑛 is equally spaced, and 𝜎 is a noise level. In the simplest model we suppose that 𝑒(𝑛) is a Gaussian white noise 𝑁(𝜇,𝜎²)=𝑁(0,1). The denoising objective is to suppress the noise part of the noisy ECG signal and to recover the clean ECG. From a statistical viewpoint, the model is a regression model over time, and the method can be viewed as a nonparametric estimation of the function 𝑓 using orthogonal basis.

Wavelet denoising has emerged as an effective method requiring no complex treatment of the noisy signal [15]. It is due to the sparsity, localitym and multiresolution nature of the wavelet transform. The wavelet transform localizes the most important spatial and frequential features of a regular signal in a limited number of wavelet coefficients. Moreover, the orthogonal transform of stationary white noise results in stationary white noise. Thus, in the wavelet domain the random noise is spread fairly uniformly among all detail coefficients. On the other hand, the signal is represented by a small number of nonzero coefficients with relatively larger values. This sparsity property assures that wavelet shrinkage can reduce noise effectively while preserving the sharp features (peaks of QRS complex). The general wavelet denoising procedure basically proceeds in three steps [11].(1)Decomposition: apply wavelet transform to the noisy signal 𝑠(𝑛) to produce the noisy wavelet coefficients to the level 𝑁, by which we can properly distinguish the presence of partial discharges.(2)Threshold detail coefficients: for each level from 1 to 𝑁, select appropriate threshold limit and apply soft or hard thresholding to the detail coefficients at some particular levels to best remove the noise. (3)Reconstruction: compute wavelet reconstruction, based on the original approximation coefficients of level 𝑁 and the modified detail coefficients of levels from 1 to 𝑁, in order to obtain estimated/smoothed signal 𝑓(𝑛).Although the application of this denoising method is not conceptually complex, some issues are carefully studied and addressed in the following for getting satisfactory results.

3. One-Step DWT-Based Baseline Wandering (BW) and Noise Suppression Proposal

Due to the similar wavelet structure for the application of BW and noise suppression, we propose here to apply in only one step both wavelet-based techniques. It will save important resources and/or time, which would facilitate any future portable hardware implementation [16, 17]. The required steps for the application of this approach are as follows.(1)Decomposition: the wavelet decomposition is applied down to a certain level 𝐿 in order to produce the approximation coefficients 𝑎_𝑛^(𝐿) that capture the BW.(2)Zeroing approximations: the approximation 𝑎_𝑛^(𝐿) is replaced by all-zero vector.(3)Threshold details: the level 𝑀 (with 𝑀<𝐿) allowing to properly distinguish the presence of partial discharges in the noisy details must be selected. Additionally, for each level from 𝑖=1 to 𝑀, the appropriate threshold limit and rule (soft or hard) are applied to the detail coefficients 𝑑_𝑛^(𝑖) for better removing the noise. (4)Reconstruction: the wavelet reconstruction, based on the zeroing approximations of level 𝐿, the modified details of levels 1 to 𝑀, and original details of levels from 𝑀+1 to 𝐿, is computed to obtain the BW corrected and denoised signal.Thus, simultaneous BW and noise suppressions are easy to get using this wavelet-based technique, for which it is necessary to select the proper parameters. One of the aims of this work is to develop a mathematical processing model oriented to portable hardware implementation. Thus, a deep analysis of the involved parameters is detailed in the following.

3.1. Analysis of Parameters

3.1.1. Wavelet Family

The selection of the wavelet family has to be based on the similarities between the ECG basic structure and the wavelet functions and the type of processing to apply. For ECG wandering suppression, the selection of the wavelet family is based on the study of the wavelet that best resembles the most significant and characteristic waveform QRS of the ECG signal. Thus, the detail sequences at different levels of decomposition from 1 to 𝐿 can capture and keep the detail features that are of interest for properly reconstructing the ECG without baseline wander. This is achieved through the elimination of the approximation sequence at level 𝐿. The selection of the wavelet family for ECG denoising is made in a similar way to that for wandering removing; it is also based on the different types of wavelets and their correlation to different signals. The order-6 Daubechies wavelet is a function well suited because of its similarity to an actual ECG. Other features include order-10 Symlets.

3.1.2. Resolution Level for Wandering

Another important parameter for wavelet-based wandering suppression is the decomposition level 𝐿, so that 𝑎_𝑛^(𝐿) can capture the baseline wander without oversmoothing. To select this resolution level, it is important to take into account the maximum number of decomposition levels, 𝑁, which is determined by the length of the sampled ECG signal. Decomposition level must be a positive integer not greater than  log₂(𝐿_𝑠), where 𝐿_𝑠 is the length of the signal, two different methods can be employed for selecting the resolution level.(i)Visual inspection: it consists of plotting the approximation sequences for different decomposition levels and to select the resolution level whose approximation better captures the baseline wandering. ECG signal from lead 4 of the DaISy dataset [18] is chosen to illustrate this method. This dataset contains 8 leads of skin potential recordings of a pregnant woman. The lead recordings, three thoracic and five abdominal, were sampled at 250 sps rate and are 10-second long. Figure 2 plots the ECG signal and approximation sequences 𝑎_𝑛^(𝑖) for 𝑖=5 to 9. This figure shows that 𝑎_𝑛⁽⁶⁾ and 𝑎_𝑛⁽⁷⁾ are the approximation sequences better capturing the baseline wander. However, this method is not effective for real-time processing. (ii)Analytical calculation: this method is based on the calculation of the resolution level for wander suppression as follows. Wavelet decomposition uses half band low-pass and high-pass filters. As it was commented in Section 2, multiresolution decomposition allows representing a signal by a finite sum of components at different resolutions. Each decomposition level contains information at different frequency bands and time scales. The approximation sequence at level 𝑖,𝑎_𝑛^(𝑖), is decomposed to obtain approximation and detail at level 𝑖+1,𝑎_𝑛^(𝑖+1), and 𝑑_𝑛^(𝑖+1). Specifically, approximation 𝑎_𝑛^(𝑖+1) is the result of low-filtering approximation 𝑎_𝑛^(𝑖) and a downsampling by 2, so that frequencies in 𝑎_𝑛^(𝑖) that are above half of the highest frequency component are removed. This means that at every decomposition level the frequency band of the approximation is reduced to the half. Figure 3 shows a study of the frequency subbands of wavelet decomposition of DaISy dataset lead 4. This figure displays approximation FFTs for resolution levels for 𝑖=1 to 7. Let us consider that the most important frequency bands in baseline wander are below a certain frequency 𝑓_𝑐. For example, 𝑓_𝑐=1Hz for wandering coming from respiration (0.15–0.3 Hz) [19]. Other wandering components such as motion of the patients and instruments have higher frequencies components. To remove wandering, it should be necessary to select the resolution level such as the approximation captures the ECG components for frequencies lower than this 𝑓_𝑐. The decomposition level for wandering suppression can be calculated as follows: 𝐿_BW=intlog₂𝐹_max𝑓_𝑐, (9) where 𝐹_max is the highest frequency component. For example, the sampling frequency for DaISy dataset is 250 Hz, so 𝐹_max=125 Hz (Nyquist Theorem) and for 𝑓_𝑐=1 Hz from (9) we can obtain 𝐿_BW=7. After applying seven low-pass filters and downsampling processes, 𝑎_𝑛⁽⁷⁾ captures frequencies from 0 Hz to 0.977 Hz and is a good estimation of baseline wander. The visual inspection method also determined that 𝑎_𝑛⁽⁷⁾ captures the baseline wander. The BW corrected signal can be obtained using wavelet reconstruction based on the detail coefficients of levels from 1 to 7 and zeroing 𝑎_𝑛⁽⁷⁾. Figure 4 illustrates the proposed method for wavelet-based wandering suppression. This figure shows the original ECG signal (DaISy dataset lead 4), the estimated baseline wander with 𝑎_𝑛⁽⁷⁾, and the ECG signal obtained after zeroing 𝑎_𝑛⁽⁷⁾.

(a)

(b)

Figure 2 

ECG data with BW (lead 4 DaIsy dataset) and approximation sequences at decomposition levels from 5 to 9.

Figure 3 

FFT of ECG DaIsy dataset lead 4 ECG1 and FFTs of approximation sequences, at decomposition levels from 1 to 7 and zoom to FFTs approximation and detail level 7.

Figure 4 

Original ECG signal (DaIsy dataset lead 4), estimated BW from 𝐴7, and BW corrected signal.

3.1.3. Level of Decomposition for Denoising

The maximum level for detail thresholding, 𝑀, depends on several factors such as the SNR in the original signal or the sample rate. The level of decomposition specifies the number of levels in the discrete wavelet analysis to take into account for detail thresholding. Unlike conventional techniques, wavelet decomposition produces a family of hierarchically organized decompositions. The selection of a suitable level for the hierarchy will depend on the signal and experience. Most often, the level is chosen based on a desired low-pass cut-off frequency. Measured signal having lower SNR usually needs more levels of wavelet transform to remove most of its noise. On the other hand, a factor that also influences the optimal level decomposition for denoising is the sampling rate of ECG signal. Sharma et al. [15] use a study of the spectrum of the detail coefficient at each level for estimating the optimum level decomposition for denoising. They concluded that noise content is significant in high frequency detail subbands, while most of the spectral energy lies in low frequency subbands. For an ECG signal sampled at 𝐹_𝑠=500Hz, noise content is significant in detail subbands 𝑑_𝑛⁽¹⁾, 𝑑_𝑛⁽²⁾, and 𝑑_𝑛⁽³⁾. Thus, in order to avoid losing clinically important components of the signal, such as PQRST morphologies, only detail 𝑑_𝑛⁽¹⁾, 𝑑_𝑛⁽²⁾, and 𝑑_𝑛⁽³⁾ should be treated for denoising. Manikandan and Dandapat [20] present a wavelet energy-based diagnostic distortion measure to assess the reconstructed signal quality for ECG compression algorithms. Their work includes a study showing, for a given sampling frequency, the information of the ECG signal and its energy contribution at each frequency subband. This helps to understand how the noise perturbs the detail coefficients and the effect over the energy spectrum. The authors also show several examples that clarify the effect of applying “zeroing” of detail coefficients at different decomposition levels. As an example to clarify the correct selection of this parameter, Figure 5 shows 4th level wavelet decomposition for ECG signal DaISy dataset lead 4. It can be observed that small subbands reflect the high frequency components of the signal, and large subbands reflect the low frequency components of the signal. The effect of high frequency artifacts can be seen in detail subbands 𝑑_𝑛⁽¹⁾ and 𝑑_𝑛⁽²⁾. These bands are weighted with small values because the energy contribution to the spectrum is low. For this example, it should be advisable to treat only detail subbands 𝑑_𝑛⁽¹⁾ and 𝑑_𝑛⁽²⁾ for wavelet denoising in order to maintain the main features of the ECG signal. It is due, as the figure shows, to the fact that detail subbands higher than level 2 contain most of the significant information for diagnostic. A more detailed study for the determination of the optimum level for wavelet-based denoising is shown in Section 4.

Figure 5 

fourth level wavelet decomposition for DaISy dataset lead 4.

3.1.4. Threshold Limits

Most algorithms are based on the previous threshold definition established in Donogo’s universal theory [21]. Since this work, modified versions of the universal threshold and new thresholds have been proposed. Matlab functions for denoising (wden.m, thselect.m) [22] and Labview blocks (Wavelet Denoise Express VI) [23] establish some of these thresholds as predefined options: rigsure, sqtwlog, heursure, and minimax. Taking this into account, the aim of this work is to develop processing models for portable computing implementations, so a study of the proposed threshold based on the following criteria was made: computational complexity, delays, latency, and clock frequency. The evaluation of these criteria derives the following threshold classification.

Prefixed Thresholds. For a given signal length, 𝑁, and a level of decomposition for wavelet denoising, 𝐽, it will be possible to know the coefficient length for each level 𝑗, denoted by 𝑛_𝑗, and thus, some thresholds proposed in the literature [13] could be calculated using software tools and stored in memory for later usage. Some of them are the following:(i)universal threshold: Th_uni=2log𝑁, (10)(ii)universal threshold level dependent: Th_uni,𝑗=2log𝑛_𝑗, (11)(iii)universal modified threshold level dependent: Th_uni,mod,𝑗=2log𝑛_𝑗𝑛_𝑗, (12)(iv)exponential threshold: Th_exp=2^{((𝑗−𝐽)/2)}2log𝑁, (13)(v)exponential threshold level dependent: Th_exp,𝑗=2^{((𝑗−𝐽)/2)}2log𝑛_𝑗, (14)(vi)minimax threshold (RefMatlab): Th_minimax𝑛=0.3936+0.1829∗log_𝑗(2).log (15) Minimax threshold uses the minimax principle to estimate the threshold [24].

Nonprefixed Thresholds. There are more threshold proposals that have been positively evaluated [15, 22] and, in some cases, get better results than the prefixed thresholds. However, a previous software calculation would not be possible since samples and/or coefficient data which are initially unknown are needed. On the other hand, the estimation of these thresholds could imply complex operations, and thus an important increment of latency. The following are two examples of this type of thresholds:(i)maxcoef threshold [22]: Th_maxcoef=2^𝑛−𝐽,𝑛=roundlog₂max𝑐_𝑗, (16)(ii)Kurtosis and ECE-based thresholds [15]: Th_DF,𝑗=1𝜖_𝑗×max𝑐_𝑗𝐹_jSN,ThDF_𝑗=1𝜖_𝐽×max𝑐_𝑗𝐹_jSN, (17) where 𝜖_𝑗 is the energy contribution efficiency of 𝑗th detail subband, 𝜖_𝐽 is detail energy contribution efficiency of 𝑗th wavelet subband, and 𝐹_jSN is the ratio between the Kurtosis value of the signal at subband 𝑗 to the Kurtosis value of Gaussian noise: 𝜖_𝑗=𝐸_𝑗𝐸_𝑡𝜖×100,_𝐽=𝐸_𝑗𝐸_𝑡𝐾×100,_′=𝑚₄𝑚²₂=(1/𝑁)∑^𝑁_𝑖=1𝑥_𝑖−media(𝑥)²(1/𝑁)∑^𝑁_𝑖=1𝑥_𝑖−media(𝑥)²². (18)

Analyzing this type of thresholds, their data and coefficients dependence makes them different from prefixed thresholds. On the hand, the implementation of the involved operations of these nonprefixed thresholds requires a high number of operations, division between them. Thus increasing implementation complexity of these thresholds. On the other hand, it will consume a big number of clock cycles and at least two accesses to each stored sample. It could have devastating effects over the total latency and real-time processing.

3.1.5. Threshold Rescaling

For signal denoising, once the threshold to be applied is selected, this threshold is rescaled using noise variance: th_{𝑗,rescaled}=𝜎_𝑗⋅th_𝑗. (19) The noise variance is used to rescale the threshold at each level, so other important setting is related to the method for estimating the noise variance at each level. If the noise is white, the standard deviation from the wavelet coefficients at the first level can be used, and the thresholds can be updated using this value. If the noise is not white, best results for denoising are obtained when estimating the noise standard deviation at each level independently and using each one to rescale the associated wavelet coefficients.

The 𝜎 is calculated based on median absolute deviation (MAD) [15, 25, 26]: 𝜎=𝐶⋅MAD𝑐_𝑗, (20) where 𝐶 is a constant scale factor, which depends on distribution of noise. For normally distributed data 𝐶=1.4826 is the 75th percentile of the normal distribution with unity variance. On the other hand, for a univariate dataset of wavelet details at 𝑗th level 𝑐_𝑗1,𝑐_𝑗2,…,𝑐_𝑗𝑛, the MAD is defined as [27] MAD=median𝑐_𝑗−median𝑐_𝑗. (21) Thus, estimated noise variance can be written as 𝜎=1.4826⋅median𝑐_𝑗−median𝑐_𝑗. (22) Some predefined applications or functions, such as the wnoisest.m Matlab function, use a simpler expression for 𝜎: 𝜎=1.4826⋅median𝑐_𝑗. (23) It must be noted that the rescale factor 𝜎 is not a prefixed value since, as (21) and (22) show, it depends on the wavelet coefficient values. However, memory access is reduced if expression (22) is used (expression (21) implies three memory accesses, while (22) requires only one). Depending on the system architecture and data processing, a higher number of memory accesses could increase the system latency and thus jeopardize real-time processing.

3.1.6. Thresholding Rules

Thresholding can be done using soft or hard thresholding. Hard thresholding is the simplest method. It can be described as the usual process of setting to zero the elements whose absolute values are lower than the threshold. Soft thresholding is an extension of hard thresholding, first setting to zero the elements whose absolute values are lower than the threshold and then shrinking the nonzero coefficients towards 0. The hard procedure creates discontinuities, while the soft procedure does not. The soft and hard thresholding is shown in (24) and (25), respectively: 𝑑_𝑛^(𝑗)=sign𝑑_𝑛^(𝑗)𝑑_𝑛^(𝑗)−th_𝑗if𝑑_𝑛^(𝑖)≥th_𝑗0if𝑑_𝑛^(𝑖)<th_𝑗,𝑑_𝑛^(𝑗)=𝑑_𝑛^(𝑗)if𝑑_𝑛^(𝑖)≥th_𝑗0if𝑑_𝑛^(𝑖)<th_𝑗, (24)(25) where 𝑑_𝑛^(𝑗) and 𝑑_𝑛^(𝑗) represent the modified values of 𝑗th level detail coefficient based on the selected threshold and th_𝑗 and are an approximation of the detail coefficients of the de-noised transform. The new detail coefficients, 𝑑_𝑛^(𝑗) or 𝑑_𝑛^(𝑗), have to be calculated for the wavelet transform levels considered for denoising, as it was pointed previously.

4. Results

This section is devoted to analyze the one-step proposed method through quantitative parameters. When working with real noisy ECG signals, it is not trivial to calculate a parameter that provides a quantitative measure of the benefits of the applied technique. In order to better analyze our proposed one-step model, a separate study of BW and noise suppression has been carried out using synthetic ECG signals. For a quantitative evaluation of the BW suppression, we have employed three types of synthetic ECG signals. These synthetic signals were elaborated from the ecg.m Matlab function [28] and are 21.6 second long and contain 5400, 10800, and 21600 samples thus the resulting sample rates are 250, 500, and 1000 sps, respectively. Signals affected by BW are obtained adding a sine wave plus a DC level, using frequencies from 0.15 to 0.31 Hz that fit to the frequency band in real BW. Our study has estimated the BW of the signals as the approximations from level 1 to 12 and has reconstructed the signal removing the estimated BW. Table 1 resumes the main results, showing the SNRs between the synthetic and the BW corrected signal. According to the expression of 𝐿, for signals sampled at 250, 500, and 1000 sps, the adequate decomposition level for BW suppression will be 7, 8, and 9, respectively, which is corroborated by Table 1. This study also reflects that better BW suppression (higher SNR) is achieved if the signal is sampled at higher rate.

Table 1 

BW suppression analysis.

Synthetic ECG signals were also used to evaluate the performance of the noise suppression. These signals were contaminated adding Gaussian white noise. Thus, the noisy signal is processed by the proposed method to obtain the denoised signal. This scenario is used by several authors [29] and allows visual inspection and quantitative evaluation. There are several parameters to measure the quality of the denoised signal [15, 29], as the SNR Improvement Measure SNRIMP [29]. A study using 3, 4, and 5 as maximum levels for wavelet denoising 𝑀, universal, exponential, and minimax as threshold limits, and simple rescaling and soft and hard thresholding was carried out. A total of 12 wavelet functions were used for this evaluation. The study also considers two noise levels, approximately 15 dB and 25 dB and a maximum of three attempts for each case (for each attempt, all the parameters are the same, including noise level, but the noisy signal is different due to the random distribution of it over the original signal). Table 2 shows the best results for denoising. Observing these summarized results and all the generated data, we can conclude that there are no large differences for the SNRIMP values of the different wavelet functions. Comparing soft and hard thresholding, soft gets for both noise levels higher SNRIMP using a less number of levels. Regarding thresholds, Th_exp achieves best denoising if it is used along with soft thresholding, as it is the case with the combination Th_minimax and hard thresholding.

Table 2 

Denoising evaluation using synthetic ECG signals.

To study BW suppression and denoising for our one-step denoising and BW suppression proposal, visual inspection of the obtained signals is also important, and in some cases it is even more conclusive than quantitative measures. DaISy dataset and Physionet Dataset [30] are targeted for evaluating the proposed one-step BW and noise suppression. The recordings from Noninvasive Fetal ECG Database [30] have two thoracic and four abdominal channels sampled at 1 ksps, all 60 seconds long. The signal bandwidth is 0.01 Hz–100 Hz. In addition, recordings from the Abdominal and Direct Fetal Electrocardiogram Database [31] have been used. Each recording comprises four differential signals acquired from maternal abdomen and the reference direct fetal electrocardiogram registered from the fetus’s head. The fetal 𝑅-wave locations were automatically determined in the direct FECG signal by means of online analysis applied in the KOMPOREL system [31]. The recordings, sampled at 1 ksps, are 5-minute long, and the signal bandwidth is 1 Hz–150 Hz.

For these signals the selected parameters were wavelet function db6, 𝑀=3,   universal threshold, soft thresholding, simple rescaling for DaIsy dataset, and multiple rescaling for Noninvasive Fetal ECG Database. Figure 6 includes the obtained result for lead 1 and lead 2 from DaIsy dataset, with BW and noise corrected signals being shown. Figure 7 shows an example of results for ecgca 746 signal of Noninvasive Fetal ECG Database, including the detail of one of the fetal QRS complexes before and after processing. These figures show that the abdominal ECG signals are BW corrected and denoised while retaining their main characteristics, as the fetal QRS complexes, which are very important for future parameter extraction [7]. Figure 8 shows more examples of results for signals of Noninvasive Fetal ECG Database.

(a)

(b)

Figure 6 

DaIsy database signals: (a) lead 1 and noise and BW corrected signal; (b) lead 2 and noise and BW corrected signals.

(a)

(b)

(c)

Figure 7 

ecgca 746 abdomen lead 1 Noninvasive Fetal ECG Database (Physionet Dataset): (a) signal; (b) BW and noise corrected signal; (c) signals detail.

(a)

(b)

(c)

Figure 8 

Noninvasive Fetal ECG Database (Physionet Dataset): (a) ecgca 906 lead 1 signal, BW and noise corrected signal; (b) ecgca 473 lead 1 signal, BW and noise corrected signal; (c) ecgca 840 lead 1 signal, BW and noise corrected signal.

Figures 9, 10, 11, and 12 show processing examples for signals of Abdominal and Direct Fetal Electrocardiogram Database displayed. Figure 9 uses signal r1 abdomen 1 to make a study of the BW and noise corrected signals using three different levels for signal denoising, 𝑀=3,𝑀=4, and 𝑀=5; as this level increases, the denoised signal loses its main characteristics (i.e., fetal QRS complexes). Figure 10 shows the denoising of r4 signals using multiple rescaling and simple rescaling. It can be observed that the best denoising is obtained using multiple rescaling. On the other hand, Figure 11 shows the results for r7 abdomen 2 signal using soft and hard thresholding; soft thresholding gets better denoised signal since hard thresholding introduces some discontinuities into the denoised signal. Finally, universal and minimax thresholds are used for denoising r10 abdomen 2 signal, as Figure 12 shows. Similar results are obtained for these two thresholds, but observing signal details, minimax threshold provides better results.

Figure 9 

r1 abdomen 1 Abdominal and Direct Fetal Electrocardiogram Database (Physionet Dataset) signal and BW and noise corrected signals using different levels for denoising.

Figure 10 

r4 abdomen 1 Abdominal and Direct Fetal Electrocardiogram Database (Physionet Dataset) signal and BW and noise corrected signals using multiple and simple rescaling.

Figure 11 

r7 abdomen 2 Abdominal and Direct Fetal Electrocardiogram Database (Physionet Dataset) signal and BW and noise corrected signals using the two types of thresholding rules.

Figure 12 

r10 abdomen 2 Abdominal and Direct Fetal Electrocardiogram Database (Physionet Dataset) signal and BW and noise corrected signals using two types of thresholds.

5. Conclusion

This paper presents the mathematical bases for electrocardiogram signal denoising by means of discrete wavelet processing. A novel one-step wavelet-based method has been introduced performing both BW and noise suppression, which makes computationally feasible real-time implementations. The presented approach is performed by only one wavelet decomposition and reconstruction step, which is required for eliminating both types of perturbations. This approach has been linked to an exhaustive study of the related parameters, such as number of decomposition levels, threshold edges, rescaling, and rules that allow an optimal signal denoising and meeting specific ECG signal characteristics including signal shape, sample rate, and noise levels. The presented results for synthetic ECG signals validate this method, while applications on real AECG signal from three different databases have led to improved signals that are valid for further analysis and extraction of parameters such as heart rate variability. The defined algorithm will also allow its compact implementation, thus fitting portable ECG applications.

Acknowledgment

This work has been partially supported by the CEI BIOTIC Granada under project 2013/81, “Compromiso con la investigación y el desarrollo.”