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

Authors: E. Castillo; D. P. Morales; A. García; F. Martínez-Martí; L. Parrilla; A. J. Palma
Journal: Journal of Applied Mathematics
 (2013)

---

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


---

## Body

## 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]
(1)W(a,b)=∫-∞+∞f(t)ψa,b(t)dt,(2)ψa,b(t)=1aψ*(t-ba).
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→R which satisfies the following four properties.(1)
Ifm>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
(3)|ψ(j)(t)|<1t,j∈{1,…,m-1}foreach|t|>γ.(4)
Foreachj∈{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
(4)W(a,b)=C(j,k)=∑nf(n)ψj,k(n),
where
(5)ψj,k(n)=2-j/2ψ(2-jn-k)
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
(6)an(i)=∑k=0N-1gka2n-k(i-1)i=1,2,…,J,dn(i)=∑k=0N-1hka2n-k(i-1)an(0)≡xn
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)
(b)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(7)a^m(i-1)={∑k=0(N/2)-1g-2ka^(m/2)-k(i)+∑k=0(N/2)-1h-2kd^(m/2)-k(i)ddddddddddddddddddddddddddddmeven∑k=0(N/2)-1g-2k+1a^((m-1)/2)-k(i)+∑k=0(N/2)-1h-2k+1d^((m-1)/2)-k(i)ddddddddddddddddddddddddddddmodd,
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 levelL in order to produce the wavelet approximation coefficients an(L) that captures the baseline wander.(2)
Zeroing approximation coefficients: replace the approximation vectoran(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 levelL 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:
(8)s(n)=f(n)+σe(n),
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 signals(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 toN, 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 levelN 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.


## 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]
(1)W(a,b)=∫-∞+∞f(t)ψa,b(t)dt,(2)ψa,b(t)=1aψ*(t-ba).
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→R which satisfies the following four properties.(1)
Ifm>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
(3)|ψ(j)(t)|<1t,j∈{1,…,m-1}foreach|t|>γ.(4)
Foreachj∈{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
(4)W(a,b)=C(j,k)=∑nf(n)ψj,k(n),
where
(5)ψj,k(n)=2-j/2ψ(2-jn-k)
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
(6)an(i)=∑k=0N-1gka2n-k(i-1)i=1,2,…,J,dn(i)=∑k=0N-1hka2n-k(i-1)an(0)≡xn
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)
(b)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(7)a^m(i-1)={∑k=0(N/2)-1g-2ka^(m/2)-k(i)+∑k=0(N/2)-1h-2kd^(m/2)-k(i)ddddddddddddddddddddddddddddmeven∑k=0(N/2)-1g-2k+1a^((m-1)/2)-k(i)+∑k=0(N/2)-1h-2k+1d^((m-1)/2)-k(i)ddddddddddddddddddddddddddddmodd,
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 levelL in order to produce the wavelet approximation coefficients an(L) that captures the baseline wander.(2)
Zeroing approximation coefficients: replace the approximation vectoran(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 levelL 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:
(8)s(n)=f(n)+σe(n),
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 signals(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 toN, 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 levelN 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.


## 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 levelL in order to produce the wavelet approximation coefficients an(L) that captures the baseline wander.(2)
Zeroing approximation coefficients: replace the approximation vectoran(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 levelL 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:
(8)s(n)=f(n)+σe(n),
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 signals(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 toN, 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 levelN 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 levelL in order to produce the approximation coefficients an(L) that capture the BW.(2)
Zeroing approximations: the approximationan(L) is replaced by all-zero vector.(3)
Threshold details: the levelM (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 levelL, 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 thewavelet 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 levelL, 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 Section2, 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=1Hz for wandering coming from respiration (0.15–0.3 Hz) [19]. Other wandering components such as motion of the patients and instruments have higher frequencies components. To remove wandering, it should be necessary to select the resolution level such as the approximation captures the ECG components for frequencies lower than this fc. The decomposition level for wandering suppression can be calculated as follows:
(9)LBW=int[log2(Fmaxfc)],
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)
(b)Figure 3
FFT of ECG DaIsy dataset lead 4 ECG1 and FFTs of approximation sequences, at decomposition levels from 1 to 7 and zoom to FFTs approximation and detail level 7.Figure 4
Original ECG signal (DaIsy dataset lead 4), estimated BW fromA7, and BW corrected signal.


#### 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=500Hz, 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.


#### 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:(10)Thuni=2logN,(ii)
universal threshold level dependent:(11)Thuni,j=2lognj,(iii)
universal modified threshold level dependent:(12)Thuni,mod,j=2lognjnj,(iv)
exponential threshold:(13)Thexp=2((j-J)/2)2logN,(v)
exponential threshold level dependent:(14)Thexp,j=2((j-J)/2)2lognj,(vi)
minimax threshold (RefMatlab):(15)Thminimax=0.3936+0.1829*(log(nj)log(2)).
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]:
(16)Thmaxcoef=2n-J,n=round[log2(max{|cj|})],(ii)
Kurtosis and ECE-based thresholds [15]:
(17)ThDF,j=1ϵj×max(cj)FjSN,ThDFj^=1ϵJ^×max(cj)FjSN,
where ϵj is the energy contribution efficiency of jth detail subband, ϵJ^ is detail energy contribution efficiency of jth wavelet subband, and FjSN is the ratio between the Kurtosis value of the signal at subband j to the Kurtosis value of Gaussian noise:
(18)ϵj=EjEt×100,ϵJ^=EjEt^×100,K′=m4m22=(1/N)∑i=1N(xi-media(x))2((1/N)∑i=1N(xi-media(x))2)2.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:(19)thj,rescaled=σj·thj.
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]:
(20)σ=C·MAD(cj),
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 cj1,cj2,…,cjn, the MAD is defined as [27]
(21)MAD=median(|cj-median(cj)|).
Thus, estimated noise variance can be written as
(22)σ=1.4826·median(|cj-median(cj)|).
Some predefined applications or functions, such as the wnoisest.m Matlab function, use a simpler expression for σ:
(23)σ=1.4826·median(|cj|).
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:
(24)d^n(j)={sign(dn(j))(|dn(j)|-thj)if|dn(i)|≥thj0if|dn(i)|<thj,(25)d~n(j)={dn(j)if|dn(i)|≥thj0if|dn(i)|<thj,
where d^n(j) and d~n(j) represent the modified values of jth 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.


## 3.1. Analysis of Parameters

### 3.1.1. Wavelet Family

The selection of thewavelet 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 levelL, 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 Section2, 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=1Hz for wandering coming from respiration (0.15–0.3 Hz) [19]. Other wandering components such as motion of the patients and instruments have higher frequencies components. To remove wandering, it should be necessary to select the resolution level such as the approximation captures the ECG components for frequencies lower than this fc. The decomposition level for wandering suppression can be calculated as follows:
(9)LBW=int[log2(Fmaxfc)],
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)
(b)Figure 3
FFT of ECG DaIsy dataset lead 4 ECG1 and FFTs of approximation sequences, at decomposition levels from 1 to 7 and zoom to FFTs approximation and detail level 7.Figure 4
Original ECG signal (DaIsy dataset lead 4), estimated BW fromA7, and BW corrected signal.


### 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=500Hz, 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.


### 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:(10)Thuni=2logN,(ii)
universal threshold level dependent:(11)Thuni,j=2lognj,(iii)
universal modified threshold level dependent:(12)Thuni,mod,j=2lognjnj,(iv)
exponential threshold:(13)Thexp=2((j-J)/2)2logN,(v)
exponential threshold level dependent:(14)Thexp,j=2((j-J)/2)2lognj,(vi)
minimax threshold (RefMatlab):(15)Thminimax=0.3936+0.1829*(log(nj)log(2)).
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]:
(16)Thmaxcoef=2n-J,n=round[log2(max{|cj|})],(ii)
Kurtosis and ECE-based thresholds [15]:
(17)ThDF,j=1ϵj×max(cj)FjSN,ThDFj^=1ϵJ^×max(cj)FjSN,
where ϵj is the energy contribution efficiency of jth detail subband, ϵJ^ is detail energy contribution efficiency of jth wavelet subband, and FjSN is the ratio between the Kurtosis value of the signal at subband j to the Kurtosis value of Gaussian noise:
(18)ϵj=EjEt×100,ϵJ^=EjEt^×100,K′=m4m22=(1/N)∑i=1N(xi-media(x))2((1/N)∑i=1N(xi-media(x))2)2.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:(19)thj,rescaled=σj·thj.
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]:
(20)σ=C·MAD(cj),
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 cj1,cj2,…,cjn, the MAD is defined as [27]
(21)MAD=median(|cj-median(cj)|).
Thus, estimated noise variance can be written as
(22)σ=1.4826·median(|cj-median(cj)|).
Some predefined applications or functions, such as the wnoisest.m Matlab function, use a simpler expression for σ:
(23)σ=1.4826·median(|cj|).
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:
(24)d^n(j)={sign(dn(j))(|dn(j)|-thj)if|dn(i)|≥thj0if|dn(i)|<thj,(25)d~n(j)={dn(j)if|dn(i)|≥thj0if|dn(i)|<thj,
where d^n(j) and d~n(j) represent the modified values of jth 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.


## 3.1.1. Wavelet Family

The selection of thewavelet 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 levelL, 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 Section2, 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=1Hz for wandering coming from respiration (0.15–0.3 Hz) [19]. Other wandering components such as motion of the patients and instruments have higher frequencies components. To remove wandering, it should be necessary to select the resolution level such as the approximation captures the ECG components for frequencies lower than this fc. The decomposition level for wandering suppression can be calculated as follows:
(9)LBW=int[log2(Fmaxfc)],
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)
(b)Figure 3
FFT of ECG DaIsy dataset lead 4 ECG1 and FFTs of approximation sequences, at decomposition levels from 1 to 7 and zoom to FFTs approximation and detail level 7.Figure 4
Original ECG signal (DaIsy dataset lead 4), estimated BW fromA7, and BW corrected signal.


## 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=500Hz, 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.


## 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:(10)Thuni=2logN,(ii)
universal threshold level dependent:(11)Thuni,j=2lognj,(iii)
universal modified threshold level dependent:(12)Thuni,mod,j=2lognjnj,(iv)
exponential threshold:(13)Thexp=2((j-J)/2)2logN,(v)
exponential threshold level dependent:(14)Thexp,j=2((j-J)/2)2lognj,(vi)
minimax threshold (RefMatlab):(15)Thminimax=0.3936+0.1829*(log(nj)log(2)).
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]:
(16)Thmaxcoef=2n-J,n=round[log2(max{|cj|})],(ii)
Kurtosis and ECE-based thresholds [15]:
(17)ThDF,j=1ϵj×max(cj)FjSN,ThDFj^=1ϵJ^×max(cj)FjSN,
where ϵj is the energy contribution efficiency of jth detail subband, ϵJ^ is detail energy contribution efficiency of jth wavelet subband, and FjSN is the ratio between the Kurtosis value of the signal at subband j to the Kurtosis value of Gaussian noise:
(18)ϵj=EjEt×100,ϵJ^=EjEt^×100,K′=m4m22=(1/N)∑i=1N(xi-media(x))2((1/N)∑i=1N(xi-media(x))2)2.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:(19)thj,rescaled=σj·thj.
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]:
(20)σ=C·MAD(cj),
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 cj1,cj2,…,cjn, the MAD is defined as [27]
(21)MAD=median(|cj-median(cj)|).
Thus, estimated noise variance can be written as
(22)σ=1.4826·median(|cj-median(cj)|).
Some predefined applications or functions, such as the wnoisest.m Matlab function, use a simpler expression for σ:
(23)σ=1.4826·median(|cj|).
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:
(24)d^n(j)={sign(dn(j))(|dn(j)|-thj)if|dn(i)|≥thj0if|dn(i)|<thj,(25)d~n(j)={dn(j)if|dn(i)|≥thj0if|dn(i)|<thj,
where d^n(j) and d~n(j) represent the modified values of jth 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 theecg.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.807Synthetic 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, Thexp achieves best denoising if it is used along with soft thresholding, as it is the case with the combination Thminimax 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.4360To 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 rescalingfor DaIsy dataset,andmultiple rescalingfor 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)
(b)ecgca 746 abdomen lead 1 Noninvasive Fetal ECG Database (Physionet Dataset): (a) signal; (b) BW and noise corrected signal; (c) signals detail.
(a)
(b)
(c)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)
(b)
(c)Figures9, 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 rescalingandsimple rescaling. It can be observed that the best denoising is obtained using multiple rescaling. On the other hand, Figure 11 shows the results for r7 abdomen 2 signal using soft and hard thresholding; soft thresholding gets better denoised signal since hard thresholding introduces some discontinuities into the denoised signal. Finally, universal and minimax thresholds are used for denoising r10 abdomen 2 signal, as Figure 12 shows. Similar results are obtained for these two thresholds, but observing signal details, minimax threshold provides better results.Figure 9
r1 abdomen 1 Abdominal and Direct Fetal Electrocardiogram Database (Physionet Dataset) signal and BW and noise corrected signals using different levels for denoising.Figure 10
r4 abdomen 1 Abdominal and Direct Fetal Electrocardiogram Database (Physionet Dataset) signal and BW and noise corrected signals usingmultiple and simple rescaling.Figure 11
r7 abdomen 2 Abdominal and Direct Fetal Electrocardiogram Database (Physionet Dataset) signal and BW and noise corrected signals using the two types of thresholding rules.Figure 12
r10 abdomen 2 Abdominal and Direct Fetal Electrocardiogram Database (Physionet Dataset) signal and BW and noise corrected signals using two types of thresholds.


## 5. Conclusion

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



---

*Source: 763903-2013-06-13.xml*