Wavelet Based Non Linear Thresholding Techniques for Pre Processing ECG Signals (original) (raw)

International Journal of Biomedical And Advance Research

ISSN: 2229-3809 (Online)
Journal DOI:10.7439/ijbar
CODEN:IJBABN

Research Article

Wavelet Based Non Linear Thresholding Techniques for Pre Processing ECG Signals

Nagendra H∗\mathbf{H}^{*}, S Mukherjee and Vinod Kuamr
Department of Electrical Engineering,Indian Institute of Technology,Roorkee,Uttarakhand -247 667 (India)

*Correspondence Info:

Nagendra H
Department of Electrical Engineering,
Indian Institute of Technology,Roorkee,Uttarakhand -247 667 (India)
Email: hnagendral@gmail.com

Abstract

The ECG recording is highly vulnerable to various kind of noises from different sources, such as electrocardiogram (EMG), power supply hum ( 50 Hz or 60 Hz ), measuring devices such as amplifiers, ADC etc. Hence it is very difficult and challenging to interpret and analyze raw ECG data for medical applications. A number of techniques are available to deal with these types of noises efficiently both during recording and pre processing of ECG data. In this paper, four different wavelet threshold denoising techniques are proposed to deal with the issue of noises in ECG recording. Analysis and their performances have been evaluated in terms of SNR and RMSE. The standard MITBIH arrhythmia data from physionet is used for the purpose. The procedure is implemented in MATLAB environment.

Keywords: CWT, DWT, ECG, RMSE, power spectral density SNR. Wavelet transforms.

1. Introduction To Ecg Signal

The human heart consists of two upper and two lower chambers. These upper and lower chambers are called atria and ventricles respectively. The activity of the heart is not a random process but beats in a synchronized manner. Thus the heart produces a regular rhythmic activity, called electrical wave. The depolarization cycle starts at the sinoatrial (SA) node and this wave propagates through through the atrium to atrioventricualr (AV) node and then to the bundle of His and rest of the ventricles [14][15].After depolarization phase the repolarisation state begins. This trace of electrical waves is called an electrocardiogram (ECG).The standard ECG waveform for one cardiac cycle is shown in figure 1.

The heart generates very strong electric field and hence the electrocardiogram (ECG) can be almost measured everywhere on the human body.

Figure1.ECG waveform for one cardiac cycle
img-0.jpeg

The ECG is characterized by many peaks and valley denoted by P-QRS-T and U. The P wave is associated with the activation of the atrium, the QRS complex and the T wave with ventricular repolarization and depolarization respectively. But the U wave is not consistent and is invisible among 70%70 \% of the people. Clinically the U wave is not important. The human ECG signals are in the mV range and the frequency range is 0.05−100 Hz0.05-100 \mathrm{~Hz} and most of the useful information is contained in the range of 0.5−45 Hz0.5-45 \mathrm{~Hz}. Electrocardiography is the starting point to the diagnosis of many heart disorders.

2. Fundamentals of wavelet transform

The word wavelet means a “small wave”, which is localized in time and frequency and extends from −∞-\infty to +∞+\infty for a finite duration of time 2,16,17{ }^{2,16,17}. The wavelet transform maps a 1-D function x(t)x(t) into a two dimensional time-scale plane and is denoted by
Wx(s,τ)=1s∫−∞∞x(t)hs(t−τs)dtW_{x(s, \tau)}=\frac{1}{\sqrt{s}} \int_{-\infty}^{\infty} x(t) h^{s}\left(\frac{t-\tau}{s}\right) d t
=∫x(t)hs(s,τ)(t)dt=\int x(t) h^{s}{ }_{(s, \tau)}(t) d t
here h(t)h(t) is called mother wavelet. This mother wavelet is a unscaled and untranslated function. A set of basis functions can be generated by translation and scaling the mother wavelet, called daughter wavelets and is given by
h(s,τ)(t)=1sh(t−τs)h_{(s, \tau)}(t)=\frac{1}{\sqrt{s}} h\left(\frac{t-\tau}{s}\right)
The scaling parameter ss is positive and varies from 0 to ∞\infty for s<1s<1, the transform performs compression of the signal and for s>1s>1,the transform performs the dilation of the signal. The signal x(t)x(t) can be recovered from the wavelet coefficients Wx(s,τ)W_{x}(s, \tau) by the inverse wavelet transform, and is given by the following equation.
x(t)=1c∫−∞0∞∞∫0Wx(s,τ)h(t−τs)dss2dbx(t)=\frac{1}{c} \int_{-\infty 0}^{\infty \infty} \int_{0} W_{x}(s, \tau) h\left(\frac{t-\tau}{s}\right) \frac{d s}{s^{2}} d b
where cc is unknown constant and has to be less than ∞\infty for the wavelet transform to be valid. τ\tau is the time parameter and the parameter ss is called scale, which is inversely proportional to the frequency. The wavelet transform itself varies with both ss and τ\tau. The inclusion of ss and τ\tau allows the function to be scaled and translated for different values of ss and tt. The admissibility condition is given as
c=∫−∞∞∣H(w)∣2∣w∣dw<∞c=\int_{-\infty}^{\infty} \frac{|H(w)|^{2}}{|w|} d w<\infty
The term dss2\frac{d s}{s^{2}} in equation (5) represents differential change in frequency and is obtained from differentiating the relationship between frequency and scale
dw=−dss2d w=\frac{-d s}{s^{2}}
A wavelet has to satisfy at least the following two conditions.
The integral of the function h(t)h(t) over the interval is zero.

∫−∞∞h(t)dt=0\int_{-\infty}^{\infty} h(t) d t=0

i.e. the wavelet function has an equal area above and below zero.
2. The square of h(t)h(t) has an integral equal to unity.

∫−∞∞h2(t)dt=1\int_{-\infty}^{\infty} h^{2}(t) d t=1

This function approaches zero at positive and negative infinite and decays away from origin unlike sinusoidal or other infinite waves.

2.1.The Continuous Wavelet Transforms (CWT)

The continuous wavelet transform and is defined by the equation

Wx(s,τ)=∫−∞∞x(t)ψs,t(t)dt=⟨hs,τ′′,x(t)⟩W_{x}(s, \tau)=\int_{-\infty}^{\infty} x(t) \psi_{s, t}(t) d t=\left\langle h_{s, \tau}^{\prime \prime}, x(t)\right\rangle

where the kernel functions are constructed by following equation

hs,τ(t)=1sψ∗t−τsh_{s, \tau}(t)=\frac{1}{\sqrt{s}} \psi^{*} \frac{t-\tau}{s}

The subscrip * in the above equation refers to the complex conjugate. The square of the norm is given by

∥hs,τ(t)∥2=∫−∞∞∣hs,∞(t)∣2dt=∫−∞∞∣h(t)∣2dt\left\|h_{s, \tau}(t)\right\|^{2}=\int_{-\infty}^{\infty}\left|h_{s, \infty}(t)\right|^{2} d t=\int_{-\infty}^{\infty}|h(t)|^{2} d t

This means that the dilated/compressed version of the signal has the same energy.

2.2. The Discrete wavelet transform (DWT)

In practice, when implemented on a computer, a continuous wavelet transform must be discrete. Therefore, it is necessary to discretize the continuous wavelet transform.

For DWT the scale parameter ss and translation parameter τ\tau can be taken as s=a0m□s=a_{0}^{m} \square where a0a_{0} is constant ( ao≠1a_{o} \neq 1 ) and mm is any positive integer (m∈z),τ=nb0a0m(m \in z), \tau=n b_{0} a_{0}^{m}, where b0b_{0} is any positive constant and nn is any integer n∈zn \in z. The daughter wavelet can be written as
hn,m(t)=a0−m/2h(a0−mt−nb0)h_{n, m}(t)=a_{0}^{-m / 2} h\left(a_{0}^{-m} t-n b_{0}\right)
For dyadic scaling a0=2,b0=1a_{0}=2, b_{0}=1, then above equation can be written as
hn,m=2−m/2h(2−mt−n)h_{n, m}=2^{-m / 2} h\left(2^{-m} t-n\right)
Therefore the signal x(t)x(t) can be written as

x(t)=∑n=−∞∞∑m=−∞∞Wnmhnm(t)x(t)=∑n∈z∑m∈zWnmhnm(t)\begin{aligned} & x(t)=\sum_{n=-\infty}^{\infty} \sum_{m=-\infty}^{\infty} W_{n m} h_{n m}(t) \\ & x(t)=\sum_{n \in z} \sum_{m \in z} W_{n m} h_{n m}(t) \end{aligned}

In DWT 13,17{ }^{13,17}, the filters with different cut off frequency are used for the analysis of the signal at different scales. The low pass filter is denoted by G0\mathrm{G}_{0} and high pass is denoted by H0\mathrm{H}_{0}. At each level of decomposition the high pass filter gives detailed coefficients d[n]\mathrm{d}[\mathrm{n}], which generally represent the noise of the signal and output of the low pass filter gives approximation coefficients a[n]\mathrm{a}[\mathrm{n}], also called coarse approximation. The DWT of the original signal can be obtained by concatenating all the coefficients a[n]\mathrm{a}[\mathrm{n}] and d[n]\mathrm{d}[\mathrm{n}], starting from the last level of decomposition.

The decomposition and reconstruction of the signal using wavelet transform is shown in figure (2) and (3) respectively.

img-1.jpeg

Figure (2). Three level wavelet Decomposition
img-2.jpeg

Figure (3). Three level wavelet reconstruction

2.3. Selection of Wavelet

A number of wavelet functions are generally used in data analysis.e.g.Haar wavelet, Daubechies wavelet, Mortlet, Symlet, Mexican hat, Shannon wavelet, coiflets etc. The choice of the wavelet depends upon the type of the application in hand.

2.4.Daubechies wavelets DpD_{p}

Daubechies wavelets are family of orthogonal discrete wavelet transform and characterized by maximal number of vanishing moments for a given support.Daubechies orthogonal wavelets from D2-d20 are commonly used and D4, D6, D8 are the most common. The index number refers to the number of moments. The number of vanishing moments is equal to the half of the number of coefficients. For example D2 has one vanishing moment and D2 has two vanishing moments. From the Daubeches wavelet family D4 is used because it is similar in shape of the real ECG signal.

DB4 scaling signal and wavelet have support of four points. The four values DS4 (i) for the scaling signal for daub4 are
DS4 (1) =1+342,=\frac{1+\sqrt{3}}{4 \sqrt{2}}, \quad DS4 (2) =3+342=\frac{3+\sqrt{3}}{4 \sqrt{2}},
DS4 (3) =3−342,=\frac{3-\sqrt{3}}{4 \sqrt{2}}, \quad DS4 (4) =1−342=\frac{1-\sqrt{3}}{4 \sqrt{2}}
The associated daub4 wavelets are defined by
DW (1) =1−342,=\frac{1-\sqrt{3}}{4 \sqrt{2}}, \quad DW4 (2) =3−342=\frac{\sqrt{3}-3}{4 \sqrt{2}},
DW4 (3) =3+342,=\frac{3+\sqrt{3}}{4 \sqrt{2}}, \quad DW4 (4) =−1−342=\frac{-1-\sqrt{3}}{4 \sqrt{2}}
Restrictions: p≥1\mathrm{p} \geq 1, Support: [0,2p−1][0,2 \mathrm{p}-1], where p is approximation order

3. Noise Removal Process From The Ecg Signal

Electrocardiography (ECG) is one of the main tools for observing the heart activity [3] [11]. The major challenge in using the ECG is that it has very small signal to noise ratio (SNR) and buried in wide variety of noise sources. Therefore it is essential to minimize these noises from the signal to obtain clinically related data to assess the condition of the heart. There are many methods available to remove (or minimize) these noises from the raw data. Since the traditional methods have their own limitations the advanced filtering technique called wavelet filtering for preprocessing an ECG signal are used.

The heart wave (ECG) is highly non-stationary and for the analysis of this kind of signal a wavelet transform technique is most suitable since it gives the information in both time and frequency, required for clinical IJBAR (2013) 04(08)04(08)

assessment. The ECG signal is very weak whose amplitude varies in the range of 110μ V110 \mu \mathrm{~V} to 4 mV .Hence, it is easily affected by different sources of noises. The various kind of noises that are present in the ECG signal are broadly classified into the following types [8] [9].
(1) Power line interference (50 Hz/60 Hz)(50 \mathrm{~Hz} / 60 \mathrm{~Hz})
(2) Base line wander (0.5 Hz)(0.5 \mathrm{~Hz}) muscle noise (8 Hz(8 \mathrm{~Hz} to 14 Hz$)$
(3) Internal amplifiers noise
(4) Electrosurgical noise
(5) Electrode contact noise
(6) Polarization noise etc.

Therefore denoising the ECG signal is an essential part of the analysis process.
Detection and classification of signals in the presence of noise and interference is very critical in many areas of signal processing applications, especially signals which are non stationary in their characteristics such as ECG.The signal component to be detected is represented by y(t),n(t)y(t), n(t) be the noise component, then the observed signal x(t)x(t) can be modeled by the following equation:

x(t)=v(t)+n(t)x(t)=v(t)+n(t)

The noise suppression techniques in signal processing are based on representing the signal x(t)x(t) such that it is possible to separate the noise component from the signal component. To solve these difficult problems, the conventional method of signal filtering is often ineffective. Wavelet transform is a new signal processing tool, which is is suitable for processing non-stationary signals. e.g. separating the of signal from noise.

4. Denoising by wavelet shrinkage techniques

The information is mostly contained in few coefficients. By choosing a proper wavelet the one that correlates with the signal to be detected, large wavelet coefficients values are obtained when there is signal information and much smaller coefficients when there is mostly noise. This is the fundamental principle behind the idea of wavelet thresholding, in which the signal reconstruction is achieved by using wavelet coefficients whose magnitudes are above specified threshold values. In any wavelet analysis one should choose the wavelet to match the characteristics of the signal to be analyzed 1,6,5,12{ }^{1,6,5,12}.

4.1.Selection of wavelet shrinkage function

The two basic forms for selecting the shrinkage functions are soft shrinkage and hard shrinkage and they are defined by the following equations, which are developed by Donoho and Johnston 1,7{ }^{1,7}

Soft shrinkage:
δs={0 if ∣x∣<λx−λ if x>λx+λ if x<λ}\delta^{s}=\left\{\begin{array}{ccc}0 & \text { if } & |x|<\lambda \\ x-\lambda & \text { if } & x>\lambda \\ x+\lambda & \text { if } & x<\lambda\end{array}\right\}
Hard shrinkage:
δH(x)={0 if ∣x∣<λx if ∣x∣>λ}\delta^{H}(x)=\left\{\begin{array}{lll}0 & \text { if } & |x|<\lambda \\ x & \text { if } & |x|>\lambda\end{array}\right\}
Soft threshold results in the reduction of wavelet coefficients by a value λ\lambda and for this reason it is also known as wavelet shrinkage. By applying wavelet transform the noise contribution is estimated and proper threshold is chosen. The remaining wavelet coefficients are used to reconstruct the original information signal y(t)y(t).

The decomposition and reconstruction procedure is shown in figure (2) and (3).By choosing a very large value of λ\lambda,the number of coefficients obtained will be very few, resulting in over smoothing of a signal. Thus there will be a possibility that some of the information may be lost. On the other hand, if the threshold is set to too low,

then part of the noise may be maintained at the reconstructed signal. Thus, it is very crucial to select the optimum threshold value to achieve a very good noise rejection performance without major loss in the spatial resolution.

Soft shrinkage suffers from large bias while hard shrinkage has smaller bias but has bigger variance and very sensitive to small changes in the data. In the hard shrinkage function there is a discontinuity in the shrinkage function as shown in figure (5). The coefficient values xx which are above the threshold are untouched. On the other hand the soft thresholding function is continuous since it shrinks the values above the threshold λ\lambda. The soft shrinkage function is preferred for denoising since the noise affects all coefficients.
img-3.jpeg

Figure 4. Hard and Soft thresholding

4.2. Steps for denoising ECG signal

  1. Transform the noisy ECG signal to wavelet domain for finding DWT coefficients of each level (sub band).
  2. Apply thresholding to obtain the estimated wavelet coefficients for each level. It is possible to use different thresholding functions.
  3. Reconstruct the denoised ECG signal from the estimated wavelet coefficients by inverse DWT.

In these three steps the most critical thing is how to select the threshold and quantify them.

4.3. Principle of wavelet threshold denoising

The method is based on taking the discrete wavelet transform (DWT) of a signal and passing this transform through a threshold, which removes the coefficients below a certain value, then taking the inverse discrete wavelet transform (IDWT).

4.4. Methods to choose threshold value

(1) Global threshold: In this method a same value of threshold λ\lambda is to be applied globally to all wavelet coefficients for different decomposition levels mm of the DWT.
(2) Level dependent threshold: In this method different threshold value λj\lambda_{j} is chosen for each resolution level jj to shrink the wavelet coefficients at that level.

Most of the wavelet coefficients containing the signal information are in the lower resolution range. Thus the wavelet coefficients due to signal have a larger magnitude than the corresponding noise coefficients. Hence the proper choice of threshold λ\lambda takes into account at least the factors, the sample size N and noise level σ2\sigma^{2}. This threshold requires the estimation of noise level σ\sigma.

4.5.Rules for selecting the threshold λi\lambda_{i}

Donoho and Johnstone (1994) [7] proposed a robust estimate method to determine the standard deviation (noise level) σ\sigma based on median absolute deviation given as
σ=median⁡{(∣wk∣:k=1,2,3,…..N/2)}0.6745\sigma=\frac{\operatorname{median}\left\{\left(\left|w_{k}\right|: k=1,2,3, \ldots . . N / 2\right)\right\}}{0.6745}
The shrinkage function determines how the thresholds are applied to the data
(i)Universal threshold method: This method is most popular in all kinds of threshold and can be expressed as follows
λ=2log⁡(N)\lambda=\sqrt{2 \log (N)}
If data is not normalized with respect to noise standard deviation, then the threshold is given by

λ=σ2log⁡(N)\lambda=\sigma \sqrt{2 \log (N)}
where NN is the sampling length of the noisy ECG signal, σ\sigma is the estimated standard deviation of zero mean of white Gaussian noise (WGN) and wkw_{k} are the detail coefficients of wavelet transform at first scale. This method gives the largest threshold value and hence results in a relatively high degree of smoothness.
(ii) Minimaxi threshold (vishu shrink) method: The minimaxi threshold minimizes a theoretical upper bound on the asymptotic risk. They are always smaller than the universal threshold for given sample size and thus results in less smoothing. An empirical method to calculate the threshold based on the standard deviation of the noise σ\sigma, is given by the following equation
λ=2log⁡(N)σEσE\lambda=\sqrt{\frac{2 \log (N)}{\sigma \mathbb{E}} \sigma \mathbb{E}}
(iii) SURE Stein’s Unbiased Risk Estimate): This method depends on the shrinkage function and threshold adapts to the signal at each multiresolution level. The threshold is based on the principle of minimizing the Stein’s unbiased \frac{\text { SURE } \text { Stein's Unbiased Risk Estimate): This method depends on the shrinkage function and threshold adapts }}{\text { to the signal at each multiresolution level. The threshold is based on the principle of minimizing the Stein's unbiased }} risk estimate (SURE) at each resolution level. The threshold on sub band S is given as

λU=arg⁡min⁡U,SS[SURE⁡(λ,ws)]∑λ,SS[SURE⁡(λ,ws)]\lambda^{U}=\arg \frac{\min _{U, S S}\left[\operatorname{SURE}\left(\lambda, w_{s}\right)\right]}{\sum_{\lambda, S S}\left[\operatorname{SURE}\left(\lambda, w_{s}\right)\right]}

Where wsw_{s} denotes detail coefficients from sub band S and SURE⁡(λ,ws)\operatorname{SURE}\left(\lambda, w_{s}\right) denotes the corresponding Stein’s unbiased estimate of the risk corresponding to a specific shrinkage function. This method performs poorly if the coefficients are very sparse.i.e., most of the coefficients at a level are zero.
(iv) SURE Shrink: This is the hybrid method combining the universal and SURE threshold, and also the best predict variable threshold selection. If SNR is small, SURE estimation will have lots of noise. If this occurs the fixed threshold can be used. This method is preferred where the wavelet coefficients decomposition is sparse.

1Ns∑n=1Ns{(wnσ)−1}≤log⁡NsNs\frac{1}{N_{s}} \sum_{n=1}^{N_{s}}\left\{\left(\frac{w_{n}}{\sigma}\right)-1\right\} \leq \frac{\log N_{s}}{\sqrt{N_{s}}}

Under the above condition the SURE shrink uses the universal threshold; otherwise the SURE threshold is used for the coefficients on sub bands. Where NsN_{s} is number of coefficients and wkw_{k} is {ws}\left\{w_{s}\right\}.

Denoising command under MATLAB [2] [4] is:
Four threshold denoising code:
xd1=wden⁡(y,rigrsure′,s′,one′,4,db4′)\mathrm{xd} 1=\operatorname{wden}\left(\mathrm{y}, \mathrm{rigrsure}^{\prime}, \mathrm{s}^{\prime}, \mathrm{one}^{\prime}, 4, \mathrm{db} 4^{\prime}\right);
xd2=wden⁡(y,heursure′,s′,one′,4,db4′)\mathrm{xd} 2=\operatorname{wden}\left(\mathrm{y}, \mathrm{heursure}^{\prime}, \mathrm{s}^{\prime}, \mathrm{one}^{\prime}, 4, \mathrm{db} 4^{\prime}\right);
xd3=wden⁡(y,sqtwolog′,s′,one′,4,db4′)\mathrm{xd} 3=\operatorname{wden}\left(\mathrm{y}, \mathrm{sqtwolog}^{\prime}, \mathrm{s}^{\prime}, \mathrm{one}^{\prime}, 4, \mathrm{db} 4^{\prime}\right);
xd4=wden⁡(y,minimaxi′,s′,one′,4,db4′)\mathrm{xd} 4=\operatorname{wden}\left(\mathrm{y}, \mathrm{minimaxi}^{\prime}, \mathrm{s}^{\prime}, \mathrm{one}^{\prime}, 4, \mathrm{db} 4^{\prime}\right);
The default threshold denoising method code:
[thr,sorh,keepapp]=ddencmp(‘den’,‘wv’,y);
xd5=wdencmp(′gbl′,y,′db4′,3,thr,sorh,keepapp);⁡\mathrm{xd} 5=\operatorname{wdencmp('gbl',y,'db4',3,thr,sorh,keepapp);}
Where yy is the signal through a low-pass filter.
xd1-xd4 are denoised results. The spectrum of the original signal is distributed in frequency range from 0 to 100 Hz and the main component is distributed in 0.05−45 Hz0.05-45 \mathrm{~Hz}.

The results of the four threshold denoising are compared and shown in figure (10). The four threshold methods namely rigrsure, heursure, universal, and minimax, have removed considerable noise present in the signal. One can very easily visualize and see that the rigrsure and minimaxi threshold method can effectively remove noise while

retain the characteristics of ECG, while the heursure and sqtwolog threshold denoising methods may remove some of the useful component from the signal.

5. Statistical Analysis

To evaluate the performances of various techniques the following equations are used in the literature. The important equations are given as below [2] [14] [10]

  1. Standard deviation ( σ\sigma ), which gives the variation of the noise through the signal and is given by

ParseError: KaTeX parse error: Expected '\right', got 'EOF' at position 85: …xt { ecg })^{2}}̲{N}}

where Secg is the original ECG signal, S′ecg\mathrm{S}^{\prime} \mathrm{ecg} is the mean of all the samples in the signal and N is length of the filtered signal.
2. The root mean square error (RMSE) of original signal and denoised signal is given by the following equation

RMSE =1N∑i=1N( Soriginal − Sdenoised )2=\sqrt{\frac{1}{N} \sum_{i=1}^{N}(\text { Soriginal }- \text { Sdenoised })^{2}}
3. Signal to Noise Ratio (SNR), which gives the information about quality of the signal. Higher the SNR better is the performance of the system and the signal to noise ratio is given by the following equation

SNR=10log⁡10[∑i=1N( filteredsignal )2∑i=1N( originalsignal − filteredsignal )2]S N R=10 \log _{10}\left[\frac{\sum_{i=1}^{N}(\text { filteredsignal })^{2}}{\sum_{i=1}^{N}(\text { originalsignal }- \text { filteredsignal })^{2}}\right]

6. Results and Discussions

The four wavelet denoising algorithms on eight different ECG signals from MIT-BIH arrhythmia data base which is freely available on physionet were tested. The sampling frequency of the sampled signal is 360 Hz with ADC resolution of 11 bits and bit rate of 3960 bps and number of samples used is 650000 . The initial value is 995 and base line wander is 1024.The MLII lead data is used in all the recordings. The insignificant data has been suppressed and the DC component is removed by subtracting the mean from the signal. This signal is treated as the original recorded signal. To this signal power line frequency of 60 Hz and Gaussian white noise is added as shown in figure (5). This signal is passed through an FIR filter of order 64 and cut off frequency of 70 Hz . The filtered response and the filtered signal and spectrum are shown in figure (6). An FIR filter is used since it is an all pole system and always stable. Most of the clinical information is contained in the frequency range of 0.5 Hz to 45 Hz , as shown in figure (15) and the results proved in table 3 from the graph. A band pass filter is designed to pass this frequency range. The command [ba]=[b \mathrm{a}]= butter (64, [0.5/ (Fs/2) 35/ (Fs/2]) provides the desired filter coefficients for 64th 64^{\text {th }} order filter. These coefficients are used in the filter command to band pass the signal between 0.5 Hz to 45 Hz , sampled at 360 Hz . An adaptive filter has been used for attenuating the 60 Hz signal and moving average filter for removing base line wander. The db4 wavelet which resembles the standard ECG signal characteristics is used in the analysis process. The performances of the methods have been evaluated in terms of SNR and RMSE. From the table 1 , it can be observed that the SNR of the rigrsure is slightly higher than the heursure and universal and almost equal to the minimaxi threshold.SNR values of universal threshold are high due to over smoothing effect.

Table1. SNR of different types of wavelet thresholding

ECG data Rigrsure Heursure Univerasal Minimaxi
100 32.5704 32.5704 32.6002 32.5704
102 36.4383 36.4357 36.4895 36.4383
104 35.1712 35.1712 35.2036 35.1712
106 33.5221 33.5267 33.5391 33.5221
108 43.2205 43.2200 43.2049 43.2205
112 34.8709 34.8709 34.8659 34.8709
114 33.5663 33.5713 33.5832 33.5663
200 32.8587 32.8600 32.8574 32.8587

Table2. RMSE of different types of wavelet thresholding

ECG data Rigrsure Heursure Univerasal Minimaxi
100 0.0015 0.0319 0.0015 0.0789
102 0.0075 0.0075 0.0075 0.0098
104 0.0044 0.0044 0.0044 0.0156
106 0.0042 0.0068 0.0042 0.0565
108 0.0174 0.0048 0.0174 0.0507
112 0.0014 0.0014 0.0014 0.0507
114 0.0067 0.0013 0.0067 0.0393
200 0.0179 0.0236 0.0179 0.0837

Table 3.Information contained in major frequency components

Frequency in Hz
10.3687
38.018438.0184
58.7558
86.4055
114.0553
127.8802
img-4.jpeg

Figure 5.Original ECG, Gaussian noise, power frequency, and noisy ECG
img-5.jpeg

Figure 6. Filter response, spectrum of filtered, noisy and smoothed ECG
img-6.jpeg

Figure 7. Signal after smoothing filter is applied
img-7.jpeg

Figure 8. Histogram of original signal
img-8.jpeg

Figure 9. Detailed coefficients at level 1,2,3,4
www.ssjournals.com

img-9.jpeg

Figure 10. Results obtained from four thresholding methods
img-10.jpeg

Figure 11.filtered signal without removed dc component
img-11.jpeg

Figure12. Welch PSD, normal distribution, scatter, quintile plots
img-12.jpeg

Figure13. PSD, Histogram, spectrum and image of filtered ECG
img-13.jpeg

Figure 14. power spectrum
img-14.jpeg

Figure15. Amplitude spectrum

7. Conclusions

In this work we have analyzed a very important signal, the electrocardiography by applying an advanced filtering tool called discrete wavelet transform. A new threshold and shrinkage functions are used to denoise the noisy ECG signal efficiently to keep it distortion free and smooth. In signal denoising applications, soft thresholding method gives better results than hard thresholding method. From simulation results we can observe that the wavelet transform can remove the noise effectively and improve the SNR .By using rigrsure threshold and minimaxi threshold we can effectively reduce the noise and retain the useful information of the ECG signal. Use of heursure and sqtwlog threshold techniques will remove a useful component of ECG signal .These methods are more harm than good. Use of default threshold methods cannot effectively remove the noise and unable to obtain the signal characteristics of the ECG signal.

References

  1. Poornachandra S.,“Wavelet based Denoising using sub band dependent threshold for ECG signal”, Digital Signal Processing, vol.18, pp. 49 - 55, 2008.
  2. Zhang xizheng, Rui yuaquing, Wangweixiong, “An New Filtering Methods in the Wavelet Domain for Bowel Sounds”, (IJACSA) International Journal of Advanced Computer Science and Applications, vol.1, No.5, pp. 26 −31,2010-31,2010,
  3. M.kania, M.Ferenice, R.maniwski, “Wavelet denoising for Multilead High Resolution ECG Signals”, Measurement Science Review, vol.7, Sec.2, No.2, 2007.
  4. Wu Wei, Cai Peisheng, “Wavelet Denoising Simulation based on MATLAB”, Signal and Electronic Engineering, vol.6, No.3, pp. 220 - 222, 2008.
  5. G Umamaheshwara Reddy,M.Murlidhar, S.Vardarajan, “ECG Denoising using Improved Thresholding based on Wavelet Transforms”,IJCSNS International Journal of Computer Science and Network Security, vol.9, pp. 221−225,2009221-225,2009.
  6. Yang Ying,Wei Yusen, “New Threshold and Shrinkage Function for ECG Signal Denoising based on wavelet Transform”, IEEE Trans.78-1-4244-2902,2009.
  7. Dono,D.L., Denoising by soft Threshholding", IEEE Trans.,Inf Theory,1995,vol.41(3),pp. 613 - 627.
  8. Banerjiee,Swati, “ECG Signal Denoising and QRS Complex detection by wavelet transform based Thresholding”, sensors and Transducers,2010.
  9. Pramodkumar, Dewanjali Agnihotri, “Biosignal Denoising via Wavelet Thresholds”,IETE journal of Research, vol 56, issue 3,May-June 2010.
  10. Jaffery Z.A.,Ahmad K.,Afroz,“Performance comparision of Wavelet Threshold Estimators for ECG Signal denoising”, Computing society,pp.248-251, IEEE ,2010 .
  11. S A Chouakri,F Bereksi-Reguig,S Ahmaidi,O Fokapu, “Wavelet Denoising of the Electrocardiogram Signal based on the Corrupted noise estimation”, computers in cardiology,vol.32,pp 1021-1024 ,IEEE 2005
  12. S.Poornachandra,N.Kumaravel, “A novel method for the elimination of power line frequency in ECG signal using hyper shrinkage function”, Digital Signal Processing,vol.18,pp.116-126, Elsevier 2008
  13. Yang Ying,Wei Yaseen, “New Threshold and Shrinkage function fot ECG Signal Denoising based on wavelet Transform”, 971-1-4244, IEEE 2009.
  14. Nagendra H,S.Mukherjee,Vinod Kumar," Application of Wavelet Techniques in ECG signal Processing:An Overview",IJEST,vol.3,No.10,pp.7432-7443,October 2011.
  15. Win Von Drongelen “Signal Processing for neuroscintists: An introduction to the analysis of physiological signals”,An Academic Press.
  16. Pascal Wallish et al,“Matlab for neuroscintists:An introduction to scintific computing in Matlab”,Academic Press. 2009
  17. Agostino Abbate,pankaj K Das et al ,“Wavelets and Sub bands :Fundamnetals and Applications”, 2002.