A Wavelet Packets Approach to Electrocardiograph Baseline Drift Cancellation

Baseline wander elimination is considered a classical problem. In electrocardiography (ECG) signals, baseline drift can influence the accurate diagnosis of heart disease such as ischemia and arrhythmia. We present a wavelet-transform- (WT-) based search algorithm using the energy of the signal in different scales to isolate baseline wander from the ECG signal. The algorithm computes wavelet packet coefficients and then in each scale the energy of the signal is calculated. Comparison is made and the branch of the wavelet binary tree corresponding to higher energy wavelet spaces is chosen. This algorithm is tested using the data record from MIT/BIH database and excellent results are obtained.


INTRODUCTION
In recent years computer aided ECG signal analysis has gained momentum and a tremendous amount of work has been carried out. One area of interest has been removing artifacts from data records [1,2]. Artifacts are the noise induced to ECG signals that result from movements of electrodes. This in turn causes deformation and change in the electrical characteristics of the skin under and around the electrodes. These electrical changes appear in the ECG as motion artifacts and baseline drifts. Baseline wanders are considered as an artifact which produces inaccurate data when measuring the ECG parameters. The ST-segment measures are especially strongly affected by this wandering. In most of the ECG recordings the respiration, electrode impedance changes due to perspiration and increased body movements are the main causes of the baseline wandering [3]. Therefore, elimination of the baseline drifts can very much change the clinical information of the ECG signal.
The frequency components of the baseline wander are usually below 0.5 Hz which is higher under stress test condition. Baseline removal has been addressed in many different ways in literature. Ensemble averaging is a classical statistical method for baseline cancellation. This is not an accurate way since characteristics of ECG signal change from beat to beat. Therefore, the consequences of mean-based algorithms are ST displacements.
In [4], baseline estimation method using cubic spline which is a portion of the Maclaurin series (higher than the 4th derivatives are neglected) is proposed. This is a third order approximation where the baseline is estimated by polynomial approximations and then subtracted from the original raw ECG signal. This is a nonlinear method, and its performance is based on estimation of reference points in the PR intervals. The main disadvantage of this method is estimating reference points that may not belong to baseline.
In [5], a linear time-varying filtering approach is undertaken to suppress the baseline drift in the ECG signal. Beat average is subtracted from the signal and then decimated. Low-pass filtering is applied to estimate the baseline wander and is interpolated. Then it is subtracted from the original signal. This is a nonlinear approach so it is complex and highly dependent to beat rate calculations and becomes less accurate in low heart rates.
Another linear filtering method has been suggested in [6]. In this approach a digital narrowband linear-phase FIR filter with cut-off frequency of 0.8 Hz is applied. Although it can be implemented in real time, this approach has many serious problems. Another filtering technique using digital and hybrid linear-phase with cut-off frequency of 0.64 Hz is used in [7]. There are numerous problems associated with filtering. First, when the FIR structure is used, the number of coefficients is too high and therefore results in a long impulse 2 International Journal of Biomedical Imaging Forward wavelet transform Inverse wavelet transform response results; second, there is an overlap on the spectral range of ECG and baseline and any cut-off frequency above 0.05 Hz does not comply with the AHA recommendation for ECG recordings [8], and removing them from ECG will introduce distortion on the ST segment as well as the QRS complex. Adaptive filtering has been shown to be useful in many areas of biomedical signal processing, where baseline wander elimination is not an exception. Application of adaptive filters in baseline drift removal can be considered as notch filtering with lower cut-off frequency being at zero. In [9] a single tap adaptive filter with a constant reference signal (DC signal) is proposed. Although it effectively removes baseline, it produces some distortion in the ST-segment due to attenuation of low-frequency components of the ECG signal. A cascade of adaptive filters is used in [10]. The first stage is the replica of the filter scheme used in [9] with a small improvement in the bandwidth. In the second stage, the reference input is changed to an impulse sequence obtained from the QRS complex. Another drawback in these filters is slow tracking properties; therefore, a smoothed ECG signal will result.
Transform techniques have been also used to remove baseline drifts. In [11], it is suggested to use short time Fourier transform (STFT) to detect the presence of baseline drifts in the ECG signal and then use a time-varying filter similar to ones used in [5] to remove it. This method has two problems. First, STFT cannot have good time and frequency resolution at the same time. Second, it will have the same problems that are associated with filtering that is already mentioned above. Discrete wavelet transform (DWT) has also been used in baseline wander removal. In [12], baseline wander is estimated from the DWT course level coefficients at level j, and then it was subtracted from the original ECG signal. In this method Symlets wavelet of the order 10 was used. In [13], wavelet coefficients are calculated up to the sixth resolution level. Then the energy levels of the first and last levels are compared. Up to two stages of filtering are carried out based on the energy levels of the wavelet coefficients.
In this paper, a wavelet-transform-based method to remove baseline wander in ECG signals is described. This method is robust to noise and is based on estimation of the energy of the baseline and it is removed through inverse discrete wavelet transform. In Section 2 a brief introduction to wavelets, its properties, and multiresolution analysis is given. In Section 3, ECG signal characteristic points related to various activities of heart are explained. In Section 4 the pro-posed method is described and its effectiveness is examined. In Section 5 the MIT-BIH ECG database records are used to evaluate the performance of the proposed algorithm. In Sections 6 and 7 discussions and conclusion are made.

WAVELETS
Wavelets are transform methods that have received great deal of attention in the field of signal processing over the past several years. A wavelet system is a set of building blocks from which one can construct or represent a signal or a function. It is a two-dimensional expansion set. The wavelet transform is a time-scale representation method that decomposes signal f (t) into basis functions of time and scale which are dilated and translated versions of a basis function ψ(t) which is called mother wavelet [14]. Translation is accomplished by considering all possible integer translations of ψ(t) and dilation is obtained by multiplying t by a scaling factor which is usually factors of 2. The following equation shows how wavelets are generated from the mother wavelet: where j indicates the resolution level and k is the translation in time. This is called dyadic scaling, since the scaling factor is taken to be 2. Wavelet decomposition is a linear expansion and it is expressed as where ϕ(t) is called the scaling function or father wavelet and c k and d j,k are the coarse and detail level expansion coefficients, respectively. In the field of signal processing, most of the results of wavelet theory are developed using filter banks. In applications one never has to deal directly with the scaling functions or wavelets, only the coefficients of the filters in the filter banks are needed. In the forward wavelet transform system, the signal is convolved with a pair of maximally decimated quadrature mirror filters (QMF) [14,15]. Let H 0 (z) and H 1 (z) be low-pass and high-pass filters, respectively, that split the signal's bandwidth to half. These filters along with decimators manipulate the forward wavelet transform. The wavelet decomposition system is shown in Figure 1. In this ccc dcc cdc cdd dcc ddc dcd ddd figure, filters F 0 (z) and F 1 (z) and the interpolator perform the reverse process in order to reconstruct the original signal f (t). These filters are related to wavelet and scaling functions as expressed below [16]: In biorthogonal filter banks, the impulse responses of synthesis filters are related to analysis filters as Theoretically, in (2) the expansion coefficients c k and d j,k are calculated from the inner product of f (t) with ϕ(t) and ψ(t), respectively. The power of wavelet transform is based on the fact that these coefficients are computed in a recursive manner. Once c k is known in a starting scale J, all the coefficients for j = J, J − 1, . . ., are found by a simple linear transform.
A wide variety of functions could be chosen as the mother wavelet as long as (5) is satisfied: Wavelets are useful in applications such as signal denoising, wave detection, data compression, feature extraction, and so forth. There are many techniques based on wavelet theory, such as wavelet packets, wavelet approximation and decomposition, discrete and continuous wavelet transform, and so forth.
In Figure 1 only one-sided expansion of the signal using coarse information was shown. This is not completely satisfactory, and more flexibility is obtained in full wavelet packet decomposition. A full binary tree for tree scale wavelet packet transform is shown in Figure 2. The projection of the func-  tion on the scaling functions and the wavelets are recursively decomposed to obtain the binary tree.
As Figure 2 shows, at each level both high-pass and lowpass filters are applied followed by decimation where both coarse and detail level coefficients are obtained. In this figure "c" and "d" stand for course and details of the signal, respectively. This configuration gives possibility of many different paths for decomposition.
The energy of a signal is given in terms of the wavelet coefficients by Parsaval's relation as We will use this equation for computing the energy of the signal in different scales, where it is explained more in later sections.

ECG SIGNAL CHARACTERISTICS
In a typical ECG signal, the dominant parts are considered as P-wave, P-R interval, QRS complex, ST-segment, T-wave, TU-segment, U-wave, and UP-segment. In Figure 3 a typical ECG signal is shown. The P-wave represents both of the artria activation or depolarization. The first half of the P-wave is the activation of the right artrium, while the second half is the activation of the artria septum and left artrium. The QRS complex represents the ventricular depolarization. The Q-wave is the first negative deflection right after the P-wave. The R-wave is the first positive deflection after the P-wave. The S-wave is the first negative deflection after the R-wave.
The T-wave is the ventricular recovery or repolarization [17].
ST-segment has clinical significance in diagnosis of myocardial ischemia and infarction. Slow fluctuation of baseline wander can influence the ST-segment and therefore incorrect diagnosis may happen. In addition, accurate baseline detection is needed for the localization of ventricular arrhythmias with body surface potential mappings as the major frequency components of ST-segment and baseline drifts are close.

METHOD
In this section we discuss a new algorithm that we propose for removal of the baseline wander in ECG signals. The morphological features of the ECG signal often change by various sources, such as 50 Hz power line, motion artifacts baseline drifts, and so forth. It is of great importance for diagnosis purposes to remove the baseline wandering, especially in the case of precise measurement of the ST-segment. This segment is a slow varying part of the ECG signal where parts of its spectral components overlap with spectral components of baseline wandering. Therefore, by simple digital filtering it is not possible to remove baseline drifts. Our proposed algorithm is based on the assumption that the baseline drift and ECG signal constitute a mixture of two independent signals, mixed in a linear fashion. Figure 4 shows the block diagram of the procedure. First, the wavelet transform of this signal is computed. The dyadic wavelet packet decomposition of the ECG signal is carried out using Daub-4 mother wavelet. It is well known that high-frequency components are mostly focused on low-level scales. Therefore, it is expected to observe the baseline drift in larger scales.
In each scale using the wavelet coefficients, the energies of the signal for both the coarse and detail levels are calculated. These energies represent the energy of the decomposed signal in assumed scales as In the above equations, E c j,k is the energy of the signal in the coarse level of scale j (low-pass filtering branch), and E d j,k is the energy in the detail level of the signal at scale j (high-pass filtering branch). In Figure 5, the power spectrum density of the ECG signal for the MIT-BIH record "103" is shown. As the figure shows, concentrations of the energy at high frequencies decrease as the scale is increased. In this figure, the wavelet coefficients d 0 to d 8 which are the details of the signal at corresponding levels, and the coefficient c 8 is the coarse or the approximation at the last scale, are plotted. The next step in the algorithm is to compare these energy levels and then choose the branch of the binary tree that has the higher energy. In this step, best basis functions for decomposition are chosen. The path for higher energy branches will be followed until a point is reached where energy differences exceed a preset threshold level, Eth. In this point the binary tree search is completed, and the baseline wander signal is retrieved by taking the inverse wavelet transform of the wavelet packet coefficients of the last scale. In order to suppress the baseline drift, the estimated baseline wander is subtracted from the original data record and a baseline wander free ECG signal is identified.
Selecting the threshold level is accomplished in the following manner. The normalized energy and normalized bandwidth of the signal are calculated. In every scale, the product of bandwidth and energy is obtained. It was experimentally determined that whenever this value is 0.1% of the total energy-bandwidth product of the signal, then the last scale has been reached.

SIMULATION RESULTS
In this section, examples are presented in order to illustrate the effectiveness of the proposed algorithm. Removing the baseline wander using the above concept was tested on ECG data records. The test database was extracted from the MIT-BIH database. Records "b108," "b210," and "b115" which are severely distorted by baseline wandering were used. Each data is 50 seconds in duration with 11 b/sample resolution and 360 Hz sampling rate. Results are shown in Figures 6  to 8. This algorithm was tested on 15 other data records from MIT-BIH, and was successfully able to remove the baseline drifts efficiently from all of the signals. In addition to tests carried out on these real data records, a low-frequency sine-wave (0.01 Hz) was also added to a baseline free ECG   signal with a ratio of ECG signal power to baseline power of −5 dB, and then algorithm was applied. Results are shown in Figure 9, which indicates that exactly the added signal was removed. This is confirmed by calculating the PRD of the  added sine-wave as  Note that in these figures, both the original ECG signal and ECG with the removed baseline drift as well as the estimated baseline are shown. This analysis is carried out for the signal of Figure 9 and results are shown in Figure 13.
We have also compared R-R intervals before and after applying the algorithm. Any change in R-R interval will cause distortion in the ST-segment and as a result clinical information of the ECG signal will change. We first identified QRS peaks of ECG records using a peak detector scheme, and then time differences between the peaks are evaluated. Results are depicted in Figure 14, which shows R-R intervals exactly before and after coinciding. Note that graphs are plotted apart from each other for comparison purposes.
In the implementation of the wavelet packet transform, in the construction of the filter bank, Daubechies db 4 mother wavelet was used.  Figure 12: (a) 20 seconds of record "b115-1" highlighted, (b) power spectrum density of record "b115-1," (c) power spectrum density of ST-segments of record "b115-1."

DISCUSSIONS
There are number of advantages in using wavelet transfer method in the analysis of ECG signals. First, there is no need for a priori information about the signal being processed. The second advantage is its easy implementation. The ability of the proposed wavelet-based algorithm was demonstrated by test signals. This algorithm could be used in multiple ECG recordings, especially in stress tests. An additional advantage of using wavelets would be utilization of WT in highfrequency noise cancellations. Daubechies wavelets are suitable for ECG signal analysis because they bear some morphological resemblance to them. The proposed method removes the components that are not correlated to ECG signal and has such characteristics that somehow are added to it. Since the spectrum of the baseline is below the spectrum of the ECG signal, therefore its energy concentration in corresponding time-scale plane does  not change much as the scale is changed in the binary tree, but the energy of the ECG signal decreases as the scale is increased. Therefore, in the binary tree search, we reach the point that the energy of the ECG signal almost vanishes (no details in that scale) but we still have considerable energy for baseline wander.
In removing the baseline wander from ECG signal, the clinical information of the ECG record must be preserved. The most important part of the ECG signal that could be affected by the baseline is the ST-segment, whose spectrums overlap. R-R interval change could be another source of distortion. Both of these have been examined using PSD calculations.
The main cause of error in our analysis could be in determining the ST-segment's length for PSD calculations. This is not fully solved because of the complicated structure of the signal.
Although the proposed method is applied only for ECG baseline suppression, it could be used in any other biological signal processing cases for unwanted low-frequency cancellation. We used this method in off-line preprocessing of the ECG signal, it could be used in real-time applications as well.
In the current study we have used a fixed threshold level to identify the energy level of the baseline wander in wavelet domain (Eth in Figure 4). Further study can be carried out for an adaptive threshold determination. Another point of further exploration could be the use of other mother wavelets such as biorthagonal wavelets.

CONCLUSION
In this paper we presented an algorithm based on WT for removing baseline drifts in ECG signals. It has been shown that the presented algorithm can eliminate baseline drifts from ECG signals without introducing any deformation to the signal and losing any clinical information.