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.

(ProQuest: ... denotes non-US-ASCII text omitted.)

Recommended by Chang-Hwan Im

1, Department of Electronics and Computer Technology, ETSIIT, University of Granada, C/ Daniel Saucedo Aranda, 18071 Granada, Spain

Received 16 March 2013; Accepted 13 May 2013

This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

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 f(t)∈^L2 (R) is defined as [11] [figure omitted; refer to PDF] In (2), * denotes complex conjugate, and a is a scale factor and b a translation factor. The normalization factor ^a-1/2 is included so that ||_ψa,b ||=||ψ|| . Thus, _ψa,b (t) 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 m is a function ψ:R[arrow right]R which satisfies the following four properties.

(1) If m>1 , then ψ is (m-1) times differentiable.

(2) ψ∈^L∞ (R) . If m>1 , for each j∈{1,...,m-1}, ^ψ(j) ∈^L∞ (R) .

(3) ψ and all its derivatives up to order (m-1) decay rapidly: for each r>0 there is a γ>0 such that [figure omitted; refer to PDF]

(4) For each j∈{0,1,...,m} , ∫^tj ψ(t)dt=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 f(n) is defined as [figure omitted; refer to PDF] where [figure omitted; refer to PDF] and a=^2j and b=^k2j , with j,k∈... . Thus, DWT is computed in practice through a set of two FIR-like filters, lowpass and highpass, at each decomposition level i , 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 [figure omitted; refer to PDF] for an N -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].

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

(a) [figure omitted; refer to PDF]

(b) [figure omitted; refer to PDF]

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 [figure omitted; refer to PDF] 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 i-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 (dbxx ), Coiflets (coifx ), Symlets (symx ), and biorthogonal (biorx _x ), where x 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 L in order to produce the wavelet approximation coefficients ^an(L) that captures the baseline wander.

(2) Zeroing approximation coefficients: replace the approximation vector ^an(L) 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 L and the detail coefficients of levels from 1 to L , 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: [figure omitted; refer to PDF] where s(n) represents the noisy signal, f(n) is an unknown, deterministic signal, time n is equally spaced, and σ is a noise level. In the simplest model we suppose that e(n) is a Gaussian white noise N(μ,^σ2 )=N(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 f 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 s(n) to produce the noisy wavelet coefficients to the level N , by which we can properly distinguish the presence of partial discharges.

(2) Threshold detail coefficients: for each level from 1 to N , 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 N and the modified detail coefficients of levels from 1 to N , in order to obtain estimated/smoothed signal f^(n) .

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 L in order to produce the approximation coefficients ^an(L) that capture the BW.

(2) Zeroing approximations: the approximation ^an(L) is replaced by all-zero vector.

(3) Threshold details: the level M (with M<L ) allowing to properly distinguish the presence of partial discharges in the noisy details must be selected. Additionally, for each level from i=1 to M , the appropriate threshold limit and rule (soft or hard) are applied to the detail coefficients ^dn(i) for better removing the noise.

(4) Reconstruction: the wavelet reconstruction, based on the zeroing approximations of level L , the modified details of levels 1 to M , and original details of levels from M+1 to L , 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 L 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 L . 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 L , so that ^an(L) 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, N , which is determined by the length of the sampled ECG signal. Decomposition level must be a positive integer not greater than _log2 (_Ls ) , where _Ls 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 ^an(i) for i=5 to 9. This figure shows that ^an(6) and ^an(7) 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 i, ^an(i) , is decomposed to obtain approximation and detail at level i+1, ^an(i+1) , and ^dn(i+1) . Specifically, approximation ^an(i+1) is the result of low-filtering approximation ^an(i) and a downsampling by 2, so that frequencies in ^an(i) 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 i=1 to 7. Let us consider that the most important frequency bands in baseline wander are below a certain frequency _fc . For example, _fc =1 Hz 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 _fc . The decomposition level for wandering suppression can be calculated as follows: [figure omitted; refer to PDF] where _Fmax is the highest frequency component. For example, the sampling frequency for DaISy dataset is 250 Hz, so _Fmax =125 Hz (Nyquist Theorem) and for _fc =1 Hz from (9) we can obtain _LBW =7 . After applying seven low-pass filters and downsampling processes, ^an(7) captures frequencies from 0 Hz to 0.977 Hz and is a good estimation of baseline wander. The visual inspection method also determined that ^an(7) 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 ^an(7) . 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 ^an(7) , and the ECG signal obtained after zeroing ^an(7) .

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

(a) [figure omitted; refer to PDF]

(b) [figure omitted; refer to PDF]

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 omitted; refer to PDF]

Figure 4: Original ECG signal (DaIsy dataset lead 4), estimated BW from _A7 , and BW corrected signal.

[figure omitted; refer to PDF]

3.1.3. Level of Decomposition for Denoising

The maximum level for detail thresholding, M , 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 _Fs =500 Hz , noise content is significant in detail subbands ^dn(1) , ^dn(2) , and ^dn(3) . Thus, in order to avoid losing clinically important components of the signal, such as PQRST morphologies, only detail ^dn(1) , ^dn(2) , and ^dn(3) 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 ^dn(1) and ^dn(2) . 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 ^dn(1) and ^dn(2) 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.

[figure omitted; refer to PDF]

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, N , and a level of decomposition for wavelet denoising, J , it will be possible to know the coefficient length for each level j , denoted by _nj , 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: [figure omitted; refer to PDF]

(ii) universal threshold level dependent: [figure omitted; refer to PDF]

(iii): universal modified threshold level dependent: [figure omitted; refer to PDF]

(iv) exponential threshold: [figure omitted; refer to PDF]

(v) exponential threshold level dependent: [figure omitted; refer to PDF]

(vi) minimax threshold (RefMatlab): [figure omitted; refer to PDF] 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]: [figure omitted; refer to PDF]

(ii) Kurtosis and ECE-based thresholds [15]: [figure omitted; refer to PDF] where _...j is the energy contribution efficiency of j th detail subband, _...J ^ is detail energy contribution efficiency of j th wavelet subband, and _FjSN is the ratio between the Kurtosis value of the signal at subband j to the Kurtosis value of Gaussian noise: [figure omitted; refer to PDF]

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: [figure omitted; refer to PDF] 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]: [figure omitted; refer to PDF] where C is a constant scale factor, which depends on distribution of noise. For normally distributed data C=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 j th level _cj1 ,_cj2 ,...,_cjn , the MAD is defined as [27] [figure omitted; refer to PDF] Thus, estimated noise variance can be written as [figure omitted; refer to PDF] Some predefined applications or functions, such as the wnoisest.m Matlab function, use a simpler expression for σ : [figure omitted; refer to PDF] 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: [figure omitted; refer to PDF] where ^{d^n(j)} and ^d~n(j) represent the modified values of j th level detail coefficient based on the selected threshold and _thj and are an approximation of the detail coefficients of the de-noised transform. The new detail coefficients, ^{d^n(j)} or ^d~n(j) , 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 L , 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.

F s (Hz)	N	F w (Hz)			SNR for estimated BW
A 4	A 5	A 6	A 7	A 8	A 9	A 10	A 11	A 12
250	5400	0.15	1.866	6.563	12.579	25.986	25.484	7.826	- 7.957	- 8.948	- 9.199
0.19	1.866	6.563	12.570	25.636	22.445	- 0.216	- 8.630	- 9.072	- 9.193
0.23	1.866	6.563	12.555	24.591	16.422	2.536	- 8.981	- 9.003	- 9.035
0.27	1.866	6.563	12.548	23.882	10.589	- 7.621	- 8.874	- 8.943	- 8.977
0.31	1.866	6.563	12.561	24.063	6.442	- 8.311	- 8.886	- 9.212	- 9.339
500	10800	0.15	0.612	4.491	10.172	16.583	30.755	24.674	5.845	- 9.762	- 10.373
0.19	0.612	4.491	10.170	16.541	29.129	21.425	- 2.123	- 10.440	- 10.866
0.23	0.612	4.491	10.168	16.485	26.036	14.864	0.825	- 10.778	- 10.796
0.27	0.612	4.491	10.168	16.469	25.625	8.869	- 9.412	- 10.669	- 10.750
0.31	0.612	4.491	10.170	16.521	25.237	4.688	- 10.16	- 10.726	- 10.064
1000	21600	0.15	0.104	1.040	5.684	11.679	18.344	33.993	23.436	4.831	- 10.805
0.19	0.104	1.040	5.684	11.678	18.319	31.447	20.792	- 3.154	- 11.545
0.23	0.104	1.040	5.684	11.677	18.271	26.531	13.913	- 0.238	- 11.897
0.27	0.104	1.040	5.684	11.677	18.255	24.680	7.822	- 10.512	- 11.771
0.31	0.104	1.040	5.684	11.679	18.300	25.393	3.629	- 11.258	- 11.807

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 M , 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, T_hexp 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.

M	Threshold	Wavelet	SNR	SNRIMP	M	Threshold	Wavelet	SNR	SNRIMP
Soft thresholding 15 dB	Hard thresholding 15 dB
3	Thexp	coif3	21.2631	6.4015	4	Thminimax	sym7	20.8648	5.9319
3	Thexp	db7	21.4798	6.3944	4	Thminimax	sym7	21.0110	5.9256
3	Thexp	coif3	21.4773	6.3919	4	Thminimax	sym7	20.7633	5.9017
3	Thexp	db7	21.2360	6.3743	5	Thminimax	sym7	20.8288	5.8958
3	Thexp	sym7	21.4286	6.3432	5	Thminimax	sym7	20.7219	5.8603
3	Thexp	db5	21.4274	6.3420	5	Thminimax	sym7	20.9024	5.8170
3	Thexp	db5	21.2008	6.3391	3	Thminimax	sym7	20.8818	5.7964
Soft thresholding 25 dB	Hard thresholding 25 dB
3	Thexp	bior3.9	29.5524	4.4960	3	Thminimax	bior6.8	29.5025	4.5809
3	Thexp	bior3.9	29.2865	4.3649	3	Thminimax	bior6.8	29.6610	4.5534
3	Thexp	sym7	29.2657	4.3440	4	Thminimax	bior6.8	29.5801	4.4725
3	Thexp	sym7	29.3678	4.3114	4	Thminimax	bior6.8	29.3823	4.4607
3	Thexp	bior3.9	29.3712	4.2784	3	Thminimax	bior6.8	29.5049	4.4485
3	Thexp	bior6.8	29.2479	4.1915	3	Thminimax	sym6	29.5035	4.4471
3	Thminimax	bior3.9	29.0937	4.1720	5	Thminimax	bior6.8	29.5436	4.4360

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 R -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, M=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.

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

(a) [figure omitted; refer to PDF]

(b) [figure omitted; refer to PDF]

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

(a) [figure omitted; refer to PDF]

(b) [figure omitted; refer to PDF]

(c) [figure omitted; refer to PDF]

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.

(a) [figure omitted; refer to PDF]

(b) [figure omitted; refer to PDF]

(c) [figure omitted; refer to PDF]

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, M=3, M=4 , and M=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 omitted; refer to PDF]

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 omitted; refer to PDF]

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 omitted; refer to PDF]

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.

[figure omitted; refer to PDF]

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."

References

[1] J. Lee, Y. Chee, I. Kim, "Personal identification based on vectorcardiogram derived from limb leads electrocardiogram," Journal of Applied Mathematics , vol. 2012, 2012.

[2] D. P. Morales, A. García, E. Castillo, M. A. Carvajal, J. Banqueri, A. J. Palma, "Flexible ECG acquisition system based on analog and digital reconfigurable devices," Sensors and Actuators A , vol. 165, no. 2, pp. 261-270, 2011.

[3] D. P. Morales, A. García, E. Castillo, M. A. Carvajal, L. Parrilla, A. J. Palma, "An application of reconfigurable technologies for non-invasive fetal heart rate extraction," Medical Engineering and Physics , 2013.

[4] S. G. Mallat, "Theory for multiresolution signal decomposition: the wavelet representation," IEEE Transactions on Pattern Analysis and Machine Intelligence , vol. 11, no. 7, pp. 674-693, 1989.

[5] E.-B. Lin, P. C. Liu, "A discrete wavelet analysis of freak waves in the ocean," Journal of Applied Mathematics , vol. 2004, no. 5, pp. 379-394, 2004.

[6] T. Abualrub, I. Sadek, M. Abukhaled, "Optimal control systems by time-dependent coefficients using cas wavelets," Journal of Applied Mathematics , vol. 2009, 2009.

[7] R. Sameni, G. D. Clifford, "A review of fetal ECG signal processing, issues and promising directions," The Open Pacing, Electrophysiology & Therapy Journal , vol. 3, no. 1, pp. 4-20, 2010.

[8] C.-T. Ku, K.-C. Hung, H.-S. Wang, Y.-S. Hung, "High efficient ECG compression based on reversible round-off non-recursive 1-D discrete periodized wavelet transform," Medical Engineering and Physics , vol. 29, no. 10, pp. 1149-1166, 2007.

[9] C. Ghule, D. G. Wakde, G. Virdi, N. R. Khodke, "Design of portable ARM processor based ECG module for 12 lead ECG data acquisition and analysis," in Proceedings of the 2nd International Conference on Biomedical and Pharmaceutical Engineering (ICBPE '09), pp. 1-8, December 2009.

[10] D. Karadaglic, M. Mirkovic, D. Milosevic, R. Stojanovic, "A FPGA system for QRS complex detection based on Integer Wavelet Transform," Measurement Science Review , vol. 11, no. 4, pp. 131-138, 2011.

[11] J. C. Goswami, A. K. Chan Fundamentals of Wavelets Theory, Algorithms, and Applications , EE.UU: John Wiley & Sons, Hoboken, NJ, USA, 2011., 2nd.

[12] P. S. Addison, "Wavelet transforms and the ECG: a review," Physiological Measurement , vol. 26, no. 5, pp. R155-R199, 2005.

[13] B. N. Singh, A. K. Tiwari, "Optimal selection of wavelet basis function applied to ECG signal denoising," Digital Signal Processing , vol. 16, no. 3, pp. 275-287, 2006.

[14] C. B. Mbachu, G. N. Onoh, E. N. Ifeagwu, S. U. Nnebe, "Processing ECG signal with Kaiser Window-Based FIR digital dilters," International Journal of Engineering Science and Technology , vol. 3, no. 8, 2011.

[15] L. N. Sharma, S. Dandapat, A. Mahanta, "ECG signal denoising using higher order statistics in Wavelet subbands," Biomedical Signal Processing and Control , vol. 5, no. 3, pp. 214-222, 2010.

[16] M. I. Ibrahimy, F. Ahmed, M. A. Mohd Ali, E. Zahedi, "Real-time signal processing for fetal heart rate monitoring," IEEE Transactions on Biomedical Engineering , vol. 50, no. 2, pp. 258-262, 2003.

[17] J. J. Oresko, Z. Jin, J. Cheng, S. Huang, Y. Sun, H. Duschl, A. C. Cheng, "A wearable smartphone-based platform for real-time cardiovascular disease detection via electrocardiogram processing," IEEE Transactions on Information Technology in Biomedicine , vol. 14, no. 3, pp. 734-740, 2010.

[18] B. De Moor Database for the identification of systems (DaISy), 2010, http://homes.esat.kuleuven.be/~smc/daisy/

[19] J. G. Webster Medical Instrumentation, Application and Design , John Wiley & Sons, 1995.

[20] M. S. Manikandan, S. Dandapat, "Wavelet energy based diagnostic distortion measure for ECG," Biomedical Signal Processing and Control , vol. 2, no. 2, pp. 80-96, 2007.

[21] D. Donoho, I. Johnstone, "Adapting to unknown smoothness via wavelet shrinkage," Journal of the American Statistical Association , vol. 90, no. 432, pp. 1200-1224, 1995.

[22] Inc. The MathWorks Denoising: Wavelet shrinkage, nonparametric regression, block thresholding, multisignal thresholding, 2013, http://www.mathworks.es/es/help/wavelet/denoising.html

[23] LabView TM Advanced Signal Processing Toolkit, Wavelet Analysis Tools User Manual, 2013, http://www.ni.com/pdf/manuals/371533a.pdf

[24] S. Sardy, "Minimax threshold for denoising complex signals with waveshrink," IEEE Transactions on Signal Processing , vol. 48, no. 4, pp. 1023-1028, 2000.

[25] P. J. Rousseeuw, C. Croux, "Alternatives to the median absolute deviation," Journal of the American Statistical Association , vol. 88, no. 424, pp. 1273-1283, 1993.

[26] I. M. Johnstone, B. W. Silverman, "Wavelet threshold estimators for data with correlated noise," Journal of the Royal Statistical Society B , vol. 59, no. 2, pp. 319-351, 1997.

[27] W. H. Swallow, F. Kianifard, "Using robust scale estimates in detecting multiple outliers in linear regression," Biometrics , vol. 52, no. 2, pp. 545-556, 1996.

[28] Matlab ECG simulation Using atlab, 2013, http://www.mathworks.com/matlabcentral/fileexchange/10858

[29] H. G. Rodney Tan, A. C. Tan, P. Y. Khong, V. H. Mok, "Best wavelet function identification system for ECG signal denoise applications," in Proceedings of the International Conference on Intelligent and Advanced Systems (ICIAS '07), pp. 631-634, November 2007.

[30] A. L. Goldberger, L. A. Amaral, L. Glass, J. M. Hausdorff, P. C. Ivanov, R. G. Mark, J. E. Mietus, G. B. Moody, C. K. Peng, H. E. Stanley, "PhysioBank, PhysioToolkit, and PhysioNet: components of a new research resource for complex physiologic signals," Circulation , vol. 101, no. 23, pp. E215-E220, 2000.

[31] J. Jezewski, A. Matonia, T. Kupka, D. Roj, R. Czabanski, "Determination of fetal heart rate from abdominal signals: evaluation of beat-to-beat accuracy in relation to the direct fetal electrocardiogram," Biomedical Engineering , vol. 57, no. 5, pp. 383-394, 2012.