Printer Friendly

Noise suppression in ECG signals through efficient one-step wavelet processing techniques.

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) [member of] [L.sup.2](R) is defined as [11]

W(a, b) = [[integral].sup.+[infinity].sub.-[infinity]] f(t) [[psi].sub.a, b](t)dt, (1)

[[psi].sub.a, b](t) = 1/[square root of [alpha]] [[psi].sup.*](t - b/a). (2)

In (2), * denotes complex conjugate, and a is a scale factor and b a translation factor. The normalization factor [a.sup.-1/2] is included so that [parallel][[psi].sub.a, b][parallel] = [parallel][psi][parallel]. Thus, [[psi].sub.a, b](t) represents a shifted and scaled version of the so-called mother wavelet [psi], which is a window function that defines the basis for the wavelet transformation. A mother wavelet [psi] of order m is a function [psi]: R [right arrow] R which satisfies the following four properties.

(1) If m > 1, then [psi] is (m - 1) times differentiable.

(2) [psi] [member of] [L.sup.[infinity]](R). If m > 1, for each j [member of] {1,..., m - 1}, [[psi].sup.(j)] [member of] [L.sup.[infinity]](R).

(3) [psi] and all its derivatives up to order (m - 1) decay rapidly: for each r > 0 there is a [gamma] > 0 such that

[absolute value of [[psi].sup.(j)](t)] < 1/t, j [member of] {1,..., m - 1} for each [absolute value of t] > [gamma]. (3)

(4) For each j [member of] {0,..., 1, m}, [integral] [t.sup.j][psi](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

W(a, b) = C (j, k) = [summation over n] f(n) [[psi].sub.j, k](n), (4)


[[psi].sub.j, k](n) = [2.sup.-j/2][psi]([2.sup.-j]n - k) (5)

and a = [2.sup.j] and b = k[2.sup.j], with j, k [member of] Z. 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 isused 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


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

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


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 [a.sup.(L).sub.n] that captures the baseline wander.

(2) Zeroing approximation coefficients: replace the approximation vector [a.sup.(L).sub.n] 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 carryout 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 fort he noisy signal is basically of the following form:

s(n) = f(n) + [sigma]e(n), (8)

where s(n) represents the noisy signal, f(n) is an unknown, deterministic signal, time n is equally spaced, and [sigma] is a noise level. In the simplest model we suppose that e(n) is a Gaussian white noise N([mu], [[sigma].sup.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 [??](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 [a.sup.(L).sub.n] that capture the BW

(2) Zeroing approximations: the approximation [a.sup.(L).sub.n] 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 [d.sup.(i).sub.n] 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 [a.sup.(L).sub.n] 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 [log.sub.2]([L.sub.s]), where [L.sub.s] 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 [a.sup.(i).sub.n] for i = 5 to 9. This figure shows that [a.sup.(6).sub.n] and [a.sup.(7).sub.n] 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, [a.sup.(i).sub.n], is decomposed to obtain approximation and detail at level i + 1, [a.sup.(i + 1).sub.n], and [d.sup.(i + 1).sub.n]. Specifically, approximation [a.sup.(i + 1).sub.n] is the result of low-filtering approximation [a.sup.(i).sub.n] and a downsampling by 2, so that frequencies in [a.sup.(i).sub.n] 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 [f.sub.c]. For example, [f.sub.c] = 1 Hz for wandering coming from respiration (0.15-0.3Hz) [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 [f.sub.c]. The decomposition level for wandering suppression can be calculated as follows:

[L.sub.BW] = int [[log.sub.2] ([F.sub.max]/[f.sub.c])], (9)

where [F.sub.max] is the highest frequency component. For example, the sampling frequency for DaISy dataset is 250 Hz, so [F.sub.max] = 125 Hz (Nyquist Theorem) and for [f.sub.c] =1 Hz from (9) we can obtain [L.sub.BW] = 7. After applying seven low-pass filters and downsampling processes, [a.sup.(7).sub.n] captures frequencies from 0 Hz to 0.977 Hz and is a good estimation of baseline wander. The visual inspection method also determined that [a.sup.(7).sub.n] 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 [a.sup.(7).sub.n]. 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 [a.sup.(7).sub.n], and the ECG signal obtained after zeroing [a.sup.(7).sub.n].

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 o fECG 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 [F.sub.s] = 500 Hz, noise content is significant in detail subbands [d.sup.(1).sub.n], [d.sup.(2).sub.n], and [d.sup.(3).sub.n]. Thus, in order to avoid losing clinically important components of the signal, such as PQRST morphologies, only detail [d.sup.(1).sub.n], [d.sup.(2).sub.n], and [d.sup.(3).sub.n] 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 [d.sup.(1).sub.n] and [d.sup.(2).sub.n]. 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 [d.sup.(1).sub.n] and [d.sup.(2).sub.n] 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.

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 [n.sub.j], 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.sub.uni] = [square root of 2 log N], (10)

(ii) universal threshold level dependent:

[Th.sub.uni, j] = [square root of 2 log [n.sub.j]], (11)

(iii) universal modified threshold level dependent:

[Th.sub.uni, mod, j] = [square root of 2 log [n.sub.j]]/[square root of [n.sub.j]], (12)

(iv) exponential threshold:

[Th.sub.exp] = [2.sup.((j - J)/2)] [square root of 2 log N], (13)

(v) exponential threshold level dependent:

[Th.sub.exp, j] = [2.sup.((j - J)/2)] [square root of 2 loh [n.sub.j]], (14)

(vi) minimax threshold (RefMatlab):

[Th.sub.minimax] = 0.3936 + 0.1829 * (log([n.sub.j])/log(2)). (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.sub.maxcoef] = [2.sup.n - J], n = round [[log.sub.2] (max{[absolute value of [c.sub.j]]})], (16)

(ii) Kurtosis and ECE-based thresholds [15]:


where [[epsilon].sub.j] is the energy contribution efficiency of jth detail subband, [??] is detail energy contribution efficiency of jth wavelet subband, and [F.sub.jSN] is the ratio between the Kurtosis value of the signal at subband j to the Kurtosis value of Gaussian noise:


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.sub.j, rescaled] = [[sigma].sub.j] x [th.sub.j]. (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 a is calculated based on median absolute deviation (MAD) [15, 25, 26]:

[sigma] = C x MAD ([c.sub.j]), (20)

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 jth level [c.sub.j1], [c.sub.j2],..., [c.sub.jn], the MAD is defined as [27]

MAD = median ([absolute value of ([c.sub.j] - median ([c.sub.j]))]). (21)

Thus, estimated noise variance can be written as

[sigma] = 1.4826 x median ([absolute value of ([c.sub.j] - median ([c.sub.j]))]). (22)

Some predefined applications or functions, such as the wnoisest.m Matlab function, use a simpler expression for [sigma]:

[sigma] = 1.4826 x median([absolute value of [c.sub.j]]). (23)

It must be noted that the rescale factor a 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:



where [[??].sup.(j).sub.n] and [[??].sup.(j).sub.n] represent the modified values of jth level detail coefficient based on the selected threshold and [th.sub.j] and are an approximation of the detail coefficients of the denoised transform. The new detail coefficients, [[??].sup.(j).sub.n] or [[??].sup.(j).sub.n] 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.

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, [Th.sub.exp] achieves best denoising if it is used along with soft thresholding, as it is the case with the combination [Th.sub.minimax] and hard thresholding.

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 1ksps, 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.

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.

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. 10.1155/2013/763903


This work has been partially supported by the CEI BIOTIC Granada under project 2013/81, "Compromiso con la investigation y el desarrollo."


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

[2] D. P. Morales, A. Garcia, E. Castillo, M. A. Carvajal, J. Banqueri, and 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. Garcia, E. Castillo, M. A. Carvajal, L. Parrilla, and 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 and 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, and M. Abukhaled, "Optimal control systems by time-dependent coefficients using cas wavelets," Journal of Applied Mathematics, vol. 2009, Article ID 636271, 10 pages, 2009.

[7] R. Sameni and 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, and 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, and 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, and 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 and A. K. Chan, Fundamentals of Wavelets Theory, Algorithms, and Applications, EE.UU: John Wiley & Sons, Hoboken, NJ, USA, 2nd edition, 2011.

[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 and 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, and 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, and 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, and 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 et al., "A wearable smartphone-based platform for real-time cardiovascular disease detection via electrocardiogram processing," IEEE Transactions on Information Technologyin Biomedicine, vol. 14, no. 3, pp. 734-740, 2010.

[18] B. De Moor, Database for the identification of systems (DaISy), 2010,

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

[20] M. S. Manikandan and 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 and 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,

[23] LabView TM, Advanced Signal Processing Toolkit, Wavelet Analysis Tools User Manual, 2013, 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 and 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 and 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 and 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,

[29] H. G. Rodney Tan, A. C. Tan, P. Y. Khong, and 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 et al., "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, and 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.

E. Castillo, D. P. Morales, A. Garcia, F. Martinez-Marti, L. Parrilla, and A. J. Palma

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

Correspondence should be addressed to D. P. Morales;

Received 16 March 2013; Accepted 13 May 2013

Academic Editor: Chang-Hwan Im

TABLE 1: BW suppression analysis.

                    [F.sub.w]          SNR for estimated BW
[F.sub.s]   N       (Hz)        [A.sub.4]   [A.sub.5]   [A.sub.6]

                    0.15          1.866       6.563      12.579
                    0.19          1.866       6.563      12.570
250         5400    0.23          1.866       6.563      12.555
                    0.27          1.866       6.563      12.548
                    0.31          1.866       6.563      12.561
                    0.15          0.612       4.491      10.172
                    0.19          0.612       4.491      10.170
500         10800   0.23          0.612       4.491      10.168
                    0.27          0.612       4.491      10.168
                    0.31          0.612       4.491      10.170
                    0.15          0.104       1.040       5.684
                    0.19          0.104       1.040       5.684
1000        21600   0.23          0.104       1.040       5.684
                    0.27          0.104       1.040       5.684
                    0.31          0.104       1.040       5.684

                    [F.sub.w]            SNR for estimated BW
[F.sub.s]   N       (Hz)        [A.sub.7]   [A.sub.8]   [A.sub.9]

                    0.15         25.986      25.484       7.826
                    0.19         25.636      22.445      -0.216
250         5400    0.23         24.591      16.422       2.536
                    0.27         23.882      10.589      -7.621
                    0.31         24.063       6.442      -8.311
                    0.15         16.583      30.755      24.674
                    0.19         16.541      29.129      21.425
500         10800   0.23         16.485      26.036      14.864
                    0.27         16.469      25.625       8.869
                    0.31         16.521      25.237       4.688
                    0.15         11.679      18.344      33.993
                    0.19         11.678      18.319      31.447
1000        21600   0.23         11.677      18.271      26.531
                    0.27         11.677      18.255      24.680
                    0.31         11.679      18.300      25.393

                    [F.sub.w]            SNR for estimated BW
[F.sub.s]   N       (Hz)        [A.sub.10]   [A.sub.11]   [A.sub.12]

                    0.15          -7.957       -8.948       -9.199
                    0.19          -8.630       -9.072       -9.193
250         5400    0.23          -8.981       -9.003       -9.035
                    0.27          -8.874       -8.943       -8.977
                    0.31          -8.886       -9.212       -9.339
                    0.15          5.845        -9.762      -10.373
                    0.19          -2.123      -10.440      -10.866
500         10800   0.23          0.825       -10.778      -10.796
                    0.27          -9.412      -10.669      -10.750
                    0.31          -10.16      -10.726      -10.064
                    0.15          23.436       4.831       -10.805
                    0.19          20.792       -3.154      -11.545
1000        21600   0.23          13.913       -0.238      -11.897
                    0.27          7.822       -10.512      -11.771
                    0.31          3.629       -11.258      -11.807

TABLE 2: Denoising evaluation using synthetic ECG signals.

M   Threshold          Wavelet   SNR       SNRIMP

               Soft thresholding 15 dB
3   [Th.sub.exp]       coif3     21.2631   6.4015
3   [Th.sub.exp]       db7       21.4798   6.3944
3   [Th.sub.exp]       coif3     21.4773   6.3919
3   [Th.sub.exp]       db7       21.2360   6.3743
3   [Th.sub.exp]       sym7      21.4286   6.3432
3   [Th.sub.exp]       db5       21.4274   6.3420
3   [Th.sub.exp]       db5       21.2008   6.3391
               Soft thresholding 25 dB
3   [Th.sub.exp]       bior3.9   29.5524   4.4960
3   [Th.sub.exp]       bior3.9   29.2865   4.3649
3   [Th.sub.exp]       sym7      29.2657   4.3440
3   [Th.sub.exp]       sym7      29.3678   4.3114
3   [Th.sub.exp]       bior3.9   29.3712   4.2784
3   [Th.sub.exp]       bior6.8   29.2479   4.1915
3   [Th.sub.minimax]   bior3.9   29.0937   4.1720

M   Threshold          Wavelet   SNR       SNRIMP

              Hard thresholding 15 dB
4   [Th.sub.minimax]   sym7      20.8648   5.9319
4   [Th.sub.minimax]   sym7      21.0110   5.9256
4   [Th.sub.minimax]   sym7      20.7633   5.9017
5   [Th.sub.minimax]   sym7      20.8288   5.8958
5   [Th.sub.minimax]   sym7      20.7219   5.8603
5   [Th.sub.minimax]   sym7      20.9024   5.8170
3   [Th.sub.minimax]   sym7      20.8818   5.7964
             Hard thresholding 25 dB
3   [Th.sub.minimax]   bior6.8   29.5025   4.5809
3   [Th.sub.minimax]   bior6.8   29.6610   4.5534
4   [Th.sub.minimax]   bior6.8   29.5801   4.4725
4   [Th.sub.minimax]   bior6.8   29.3823   4.4607
3   [Th.sub.minimax]   bior6.8   29.5049   4.4485
3   [Th.sub.minimax]   sym6      29.5035   4.4471
5   [Th.sub.minimax]   bior6.8   29.5436   4.4360
COPYRIGHT 2013 Hindawi Limited
No portion of this article can be reproduced without the express written permission from the copyright holder.
Copyright 2013 Gale, Cengage Learning. All rights reserved.

Article Details
Printer friendly Cite/link Email Feedback
Title Annotation:Research Article
Author:Castillo, E.; Morales, D.P.; Garcia, A.; Martinez-Marti, F.; Parrilla, L.; Palma, A.J.
Publication:Journal of Applied Mathematics
Article Type:Report
Date:Jan 1, 2013
Previous Article:Thermal diffusivity identification of distributed parameter systems to sea ice.
Next Article:Stability problems for Chua system with one linear control.

Terms of use | Privacy policy | Copyright © 2021 Farlex, Inc. | Feedback | For webmasters