A wavelet-based adaptive filter for removing ECG interference in EMGdi signals

Journal of Electromyography and Kinesiology 20 (2010) 542–549

Contents lists available at ScienceDirect

Journal of Electromyography and Kinesiology journal homepage: www.elsevier.com/locate/jelekin

A wavelet-based adaptive filter for removing ECG interference in EMGdi signals Choujun Zhan a, Lam Fat Yeung a,*, Zhi Yang b a b

Department of Electronic Engineering, City University of Hong Kong, Hong Kong, China School of Information Science and Technology, Sun Yat-Sen University, Guang Zhou, China

a r t i c l e

i n f o

Article history: Received 23 March 2009 Received in revised form 27 May 2009 Accepted 22 July 2009

Keywords: EMGdi ECG Wavelet Adaptive filter

a b s t r a c t Diaphragmatic electromyogram (EMGdi) signals convey important information on respiratory diseases. In this paper, an adaptive filter for removing the electrocardiographic (ECG) interference in EMGdi signals based on wavelet theory is proposed. Power spectrum analysis was performed to evaluate the proposed method. Simulation results show that the power spectral density (PSD) of the extracted EMGdi signal from an ECG corrupted signal is within 1.92% average error relative to the original EMGdi signal. Testing on clinical EMGdi data confirm that this method is also efficient in removing ECG artifacts from the corrupted clinical EMGdi signal. Crown Copyright Ó 2009 Published by Elsevier Ltd. All rights reserved.

1. Introduction Previous research shows that diaphragmatic electromyographic (EMGdi) signals convey important information about the respiratory control mechanism (Gross et al., 1979; Roberts et al., 1996). It can be used to predict electromyographic signs of diaphragmatic fatigue (Levine and Gillen, 1987). EMGdi signal is an effective indicator for breath monitoring and respiratory muscle fatigue disease detection (Deng et al., 1998); such as Chronic obstructive pulmonary disease (COPD) syndrome (Akkiraju and Reddy, 1992), which is related to fatigue of the respiratory muscles and can cause a severe disorder of the respiratory system with a high death rate (Duiverman et al., 2004). A major problem encountered in the analysis of EMGdi signals is the signal quality (Deng et al., 2000) and the contamination of the electrocardiographic (ECG) signal. Spectrum analysis indicates that the ECG spectrum extends from 1 to 75 Hz with most of its energy lying within the 1–50 Hz range, while the EMGdi spectrum extends between 0.5 and 250 Hz with most of the energy lying between 0.5 and 150 Hz (Levine et al., 1986; Schweitzer et al., 1979). Hence, EMGdi signals and ECG signals have a significant spectrum overlap between 1 and 75 Hz. Conventional digital filters may not be able to filter the ECG without introducing significant distortion to the original EMGdi signal (Akkiraju and Reddy, 1992; Liang et al., 2005). Thus, the reduction of ECG interference from an EMGdi signal is a challenging task. Various methods for removing ECG inference from EMGdi signals have been proposed (Akkiraju and Reddy, 1992; Deng et al., * Corresponding author. E-mail addresses: [email protected] (C. Zhan), eelyeung@cityu. edu.hk (L.F. Yeung).

1998, 2000; Drake and Callaghan, 2006; Liang et al., 2005; Levine et al., 1986; Marque et al., 2005; Roberts et al., 1996). The traditional ECG artifacts removing methods are high-pass filter (Panagiotacopulos et al., 1998; Schweitzer et al., 1979), gating technique and subtraction technique (Levine et al., 1986; Roberts et al., 1996). It has been demonstrated that the high-pass filter technique is not efficient in removing the ECG (Akkiraju and Reddy, 1992). The gating technique is simple to implement; however, there are two major limitations: first, segments of an EMGdi signal contaminated by ECG artifacts will be discarded, no matter how important they are, which will lead to the loss of EMGdi information; secondly, if the patient’s heart rate is high, a large proportion of the respiratory EMGdi will be gated off, which can introduce significant distortion to the original EMGdi (Levine et al., 1986). Subtraction techniques can remove the ECG artifacts without sacrificing segments of EMGdi signal corrupted by ECG. However, it needs manual selection of the ECG templates (Levine et al., 1986). Recent advances in applying adaptive filtering techniques can help to remove noise (Akkiraju and Reddy, 1992; Deng et al., 1998, 2000; Liang et al., 2005; Marque et al., 2005). However, most adaptive filter techniques require a reference ECG signal. As a result, extra channels for recording pure ECG reference signal are always needed. In addition, the computation burden of these methods is heavy, which call for high speed hardware in the online ECG cancellation process. Consequently, it is useful to develop a filtering approach which does not need extra channels for ECG reference. Furthermore, it is also efficient in computation, and is capable of removing ECG artifacts from the corrupted clinical EMGdi signals while high fidelity of EMGdi is preserved. Wavelets are widely used in signal processing and are useful in signal denoising (Carre et al., 1998; Donoho, 1994, 1995; Mallat,

1050-6411/$ - see front matter Crown Copyright Ó 2009 Published by Elsevier Ltd. All rights reserved. doi:10.1016/j.jelekin.2009.07.007

543

C. Zhan et al. / Journal of Electromyography and Kinesiology 20 (2010) 542–549

1998). It will be shown that traditional wavelet threshold denoising approaches (Carre et al., 1998; Johnstone and Silverman, 1997; Krim et al., 1999; Pan et al., 1999), extensions of the hard thresholding (Donoho, 1994) and the soft thresholding (Donoho, 1995), are inappropriate in removing ECG from EMGdi. However, wavelet theory is powerful in multiresolution analysis (MRA). In this paper, combining with MRA, a wavelet-based adaptive filter is proposed, which does not require an extra reference ECG signal and is efficient in computation. Experiments show that it needs no more than 35 min to process a 1-h-length clinical EMGdi signal with 7,200,000 data, which means that this technique can be used in on-line ECG elimination. Simulations were initially used to evaluate the performance of this method. To quantify the difference between the simulated pure EMGdi signal and the filtered/extracted EMGdi signal, power spectrum analysis was applied. Experimental results show that total power of EMGdi signals were reduced from 2050 mv2 =Hz, before processing, to 224 mv2 =Hz, after processing. The average relative error between the power spectral density (PSD) of the simulated pure EMGdi and the extracted EMGdi with ECG contamination were within 1.92%, which is suitable for clinical application. Then clinical EMGdi data were processed using this wavelet-based adaptive filter and the results agreed with the simulated ones. The paper is organized as follows: in Section 2, the limitation of traditional threshold denoising approaches in removing the ECG artifacts from EMGdi signal will be reviewed, and the new wavelet-based adaptive filter will be introduced; in Section 3, experiment results will be presented; finally, conclusions will be made in Section 4. 2. ECG elimination

system (GT201/50 produced by ADInstruments, Castle Hill, Australia) which has a 0.5–1000 Hz bandwidth. The raw EMGdi signals were digitized at a sampling frequency of 2000 Hz for each channel by a 16-bit ADC (ML870 PowerLab 8/30 produced by ADInstruments). Finally, the data were stored and processed using a Pentium Dual Core computer (2.13 GHz  2). 2.2. The wavelet transform The method of wavelet-threshold denoising is based on the principle of multiresolution analysis (MRA) (Mallat, 1998). We attempt to recover a signal f ½n from the observation data z½n, which is corrupted by noise w½n

z½n ¼ f ½n þ w½n;

n ¼ 1; 2; . . . ; N:

ð1Þ

The discrete wavelet decomposition of z½n can be expressed by

z½n ¼

X

h½n  2ia1i þ

i

X

1

g½n  2idi :

ð2Þ

i

Let us introduce the fast filter bank algorithm that computes the orthogonal wavelet coefficients of the observation signal z½n. At decomposition level 1, one has

8 1 P > < ai ¼ h½n  2iz½n; n P 1 > : di ¼ g½n  2iz½n:

ð3Þ

n

At decomposition level j þ 1, one has the decomposition Eq. (4)

8 ðjþ1Þ P > ¼ h½n  2iajn ; j ¼ 1; 2; . . . ; M  1; < ai n P ðjþ1Þ > : di ¼ g½n  2iajn ; j ¼ 1; 2; . . . ; M  1:

ð4Þ

n

2.1. Measurement of EMGdi signals

This process splits the aj frequency information roughly in half, par-

Raw EMGdi signals were obtained via a multiple-array esophageal electrode consisting of seven copper rings (2 mm wide and 2 mm in diameter) placed 10 mm apart (shown in Fig. 1a). The multiple-array esophageal electrode was placed via the nose into the esophagus and stomach of the patient (shown in Fig. 1b). The raw EMGdi signals were amplified 1000 times by the PowerLab

titioning it into a set of detail coefficients d and a set of approximation coefficients aðjþ1Þ . This procedure continues iteratively until the desired M-level decomposition is obtained. h½n and g½n are the low-pass and high-pass filters, corresponding to the selected wavelet basis, respectively. For example, with db4 wavelet basis,    pffiffiffi pffiffiffi pffiffiffi pffiffiffi pffiffiffi pffiffiffi h½1 ¼ 1 þ 3 4 2, h½2 ¼ 3 þ 3 4 2, h½3 ¼ 3  3 4 2,

ðjþ1Þ

Fig. 1. (a) Electrodes placement and (b) flow chart of the EMGdi signal recording and denoising procedures.

544

C. Zhan et al. / Journal of Electromyography and Kinesiology 20 (2010) 542–549



pffiffiffi pffiffiffi h½4 ¼ 1  3 4 2, h½n ¼ 0, for n – 0; 1; 2; 3, and g½1 ¼ h½4, g½2 ¼ h½3, g½3 ¼ h½2, g½4 ¼ h½1. A simple discrete wavelet transform (DWT) example is provided in the Appendix and more detail on DWT can be found in Chapter VII of Mallat (1998). For orthogonal wavelets, to reconstruct a signal f ½n from a gi~j ðj ¼ 1; 2; . . . ; MÞ, one can ~M and d ven wavelet coefficients series a use the following reconstruction formula:

~ji ¼ a

X

~nðjþ1Þ þ h½i  2na

X

n

~ ðjþ1Þ g½i  2nd n

for j ¼ 1; 2; . . . ; M  1

n

ð5Þ and

f ½i ¼

X

~1n þ h½i  2na

X

n

~1 : g½i  2nd n

ð6Þ

n

Based on the white noise assumption, Donoho (1994) suggests truncating small wavelet coefficients to 0, as they are dominated by noise and carrying only a small amount of signal information (Grossman and Morlet, 1984; Mallat, 1998). Hence, in order to reconstruct f ½n from z½n, for each decomposition level j, a threshold T j is selected. If the wavelet coefficients are smaller than the threshold T j , these coefficients can be ignored or shirked and the attenu~j ðj ¼ 1; 2; . . . ; MÞ are given by (7) or (8). ~M and d ated coefficients a Then, the original signal can be reconstructed from the attenuated wavelet coefficients using (5) and (6). Donoho and Johnstone proposed two types of shrinkage functions:  ‘‘Hard Thresholding” (Deng et al., 1998), defined as,

 hðxÞ ¼

0 if jxj 6 T j ; x if jxj > T j :

ð7Þ

 ‘‘Soft Thresholding” (Carre et al., 1998) defined as,

 hðxÞ ¼

0

if jxj 6 T j ;

ð8Þ

sgnðxÞ  jx  T j if jxj > T j ;

man, 1997; Krim et al., 1999; Pan et al., 1999). These methods share the same basic idea: truncating or shrinking small wavelet coefficients. The major differences between these methods are the selection of wavelet basis (Carre et al., 1998; Krim et al., 1999) or threshold T j (Johnstone and Silverman, 1997; Pan et al., 1999). We could classify these wavelet thresholding methods as ‘‘small-coefficient-truncated” methods. These methods (Carre et al., 1998; Johnstone and Silverman, 1997; Krim et al., 1999; Pan et al., 1999) are successful in extracting the original signal subject to a white noise contamination. However, in this application the contaminant is an ECG signal (which is far from the white noise assumption). Most energy of the ECG contamination is located on the QRS complex which has a much larger amplitude than the EMGdi signal. On the contrary, EMGdi signals are similar to a band-limited noise, as shown in Fig. 2; actually, EMGdi can be simulated by a zero mean Guassian white noise passed through a band-pass filter (Stashuk, 1993). Based on MRA, we find that: the larger the coefficient is, the more likely it is dominated by ECG signal. Hence, the ‘‘small-coefficienttruncated” methods will remove EMGdi signal but the ECG contamination will be preserved, just as show in Fig. 3d. 2.3. Wavelet-based adaptive filter Here, a new wavelet-based adaptive filter is proposed. The concept is mainly based on the MRA results: most ECG energy are contributed to large wavelet coefficients and located on small-scale spaces. Thus, with a proper threshold, we can remove most of ECG artifacts by truncating or shrinking the ‘‘larger” wavelet coefficients (which are the dominant part of the ECG interference). This is the first difference between the proposed method and traditional methods. Here, the proposed method can be classified as ‘‘largecoefficient-truncated” method. Compared with (7) and (8), the new ‘‘large-coefficient-truncated” shrinkage functions are defined as:



where x is a wavelet coefficient.

hðxÞ ¼

Based on Donoho’s works, many wavelet-based denoising methods were developed (Carre et al., 1998; Johnstone and Silver-

0 if jxj > T j ðxÞ; x if jxj 6 T j ðxÞ;

ð9Þ

where x is wavelet coefficient, T j ðxÞ is the adaptive threshold of x.

mv

(a) 400 200 0 −200 −400

0

1000

2000

3000

4000

5000

(b)

mv

100 0 −100

0

500

1000

1500

2000 (c)

2500

3000

3500

4000

mv

100 0 −100

0

500

1000

1500

2000

2500

Fig. 2. (a) ECG signal; (b) EMGdi signal; (c) Gaussian white noise signal.

3000

545

mv

C. Zhan et al. / Journal of Electromyography and Kinesiology 20 (2010) 542–549

600 400 200 0 −200 −400

5

10

15

20

(b)

5

10

15

20

15

20

15

20

(f)

0 −100 5

mv

−200 100

0 −100

mv

(e)

0

100 mv

200

(a)

600 400 200 0 −200 −400

10

15

20

(c)

5

10 (g)

100 0 −100

5

600 400 200 0 −200 −400

10

15

20

(d)

5

10 (h)

100 0 −100

5

10 15 Time (sec)

20

5

10 15 Time (sec)

20

Fig. 3. A simulation example: (a) the simulated ECG-interferent signal; (b) the simulated pure EMGdi signal; (c) the corrupted EMGdi signal; (d) the processed EMGdi signal using hard thresholding method; (e) processed result by ‘‘inverse” hard thresholding method; (f) processed result by high-pass filter; (g) processed result by gating technique; and (h) processed result by wavelet-based adaptive filter.

Note that we use the adaptive threshold T j ðxÞ instead of the constant threshold T j at each level. For a constant threshold, if the threshold is low then both the ECG artifacts together with the EMGdi signals will be removed. On the contrary, if the constant threshold is large enough to keep most of the EMGdi signals, a lot of ECG noise will be kept. Hence, a fixed threshold is inadequate for achieving a good performance, as shown in Fig. 3e. Therefore, an adaptive thresholding method is suggested and this is the second difference. A simple adaptive threshold will be use for illustration by defining it as a gain K multiplying an averaged value of 2L neighbor samples, with upper bound U b and lower bound Lb , i.e.

T j ðxÞ ¼ K  v j ðxÞ;

ð10Þ

where

8 PiþU PiLb b jaMn jþ n¼iU jaMn j > > b < n¼iþLb if 2ðU b Lb þ1Þ v j ðxÞ ¼ > PiþUb j PiLb j > : n¼iþLb jdn jþ n¼iUb jdn j if 2ðU L þ1Þ b

b

x ¼ aM i ; j

x ¼ di

for j ¼ 1; 2; . . . ; M: ð11Þ

Finally, we arrive the following wavelet-based adaptive filter algorithm: A1. Algorithm of the wavelet-based adaptive filter INPUT: EMGdi signals contaminated by ECG interference. OUTPUT: Denoised EMGdi signals. Step 1: Compute the scale decomposition using (3) and (4) j and obtain wavelet coefficients di and aji . Step 2: Use 2L sampling points corresponding to two intervals. Calculate the average value v i and adaptive thresholding T j ðxÞ using (10) and (11). Step 3: Compute the attenuated decomposition coefficients ~j using (9). ~ji and d a i Step 4: Reconstruct the EMGdi signal ~f ½n based on the atten~j using reconstruc~ji and d uated decomposition coefficients a i tion formula (5) and (6). End.

3. Experiment result In order to check the performance of our approach, both computer simulation EMGdi and clinical EMGdi signals are processed by the wavelet-based adaptive filter. 3.1. Simulation result Based on the EMGdi signal generation procedure (Deng et al., 2000), a pure EMGdi signal was generated in two steps: first, a zero mean Gaussian white noise was generated and then modulated by a cosine function; second, the cosine-modulated Gaussian noise was filtered by a filter with the appropriate spectral transfer, as shown in Fig. 3b. The simulated ECG signal (Fig. 3a) was generated by Matlab simulation program (Karthik, XXXX). The corrupted EMGdi signal with ECG inference (Fig. 3c) was generated by adding the simulated pure EMGdi and simulated ECG signals. The processed results of the corrupted EMGdi signal using different approaches are shown in Fig. 3d–h. Fig. 3d shows the processed result by the hard thresholding method, which belongs to the ‘‘small-coefficient-truncated” methods. Note that most of ECG artifacts were not removed. However, most of the EMGdi signal was removed, as the EMGdi signal was treated as white noise and was removed (refer to Section 2.2). For fair comparison, we changed the hard thresholding method a little by using the following shrinkage functions:

 hðxÞ ¼

0 if jxj > T j ; if jxj 6 T j :

x

ð12Þ

Compared with shrinkage functions of hard thresholding (7), which shrinks small x reserves large x, note that (12) preserves the small x and shrinks large x, but the threshold is constant for each level. Here, we call this method the ‘‘inverse” hard thresholding method. The processed result by the ‘‘inverse” hard thresholding is shown in Fig. 3e. One can see that most of ECG artifacts were removed but there are still a lot of ‘spikes’, just like the processed result by high-pass filter (Fig. 3f). The gating technique seems to be able to remove the ECG inference; however, the EMGdi part corrupted by

546

C. Zhan et al. / Journal of Electromyography and Kinesiology 20 (2010) 542–549

120

80

0

4

0

PSD( mv /Hz)

0

0

0

100

5 2

2

20

5

50

60

40

0

50

50

100

150

100 150 Frequency (Hz)

200

250

200

(d)

2

(b)

6

10

(c)

100

2

PSD( mv /Hz)

100

PSD( mv /Hz)

(a) Processed EMGdi Simulated pure EMGdi Simulated corrupted EMGdi

250

0

200

4

3

3

2

2

1

1 0

100

5

(e)

4

0

0

0

100 200 Frequency (Hz)

200 (f)

0

100 200 Frequency (Hz)

Fig. 4. PSD of signal processed by different approaches: (a) PSD of signal processed by the wavelet-based adaptive filter; (b) displays same data as (a) on expanded scale; (c) hard thresholding method; (d) ‘‘inverse” hard thresholding method; (e) high-pass filter; and (f) gating technique.

ECG was also removed at the same time, as shown in Fig. 3g. The processed EMGdi signal using the wavelet-based adaptive filter is shown Fig. 3h, which is close to the simulation pure EMGdi signal (Fig. 3b). Power spectral analysis (PSD) is the commonly used quantitative method for describing an EMG signal (Akkiraju and Reddy, 1992; Deng et al., 1998; Levine et al., 1986; Roberts et al., 1996; Schweitzer et al., 1979). To evaluate the performances of different methods, the PSD of the simulated pure and corrupted EMGdi be-

fore and after process by various methods are given. Fig. 4a shows that maximum amplitude of the PSD of the corrupted EMGdi is decreased from almost 120 mv2 =Hz to about 5 mv2 =Hz after processing by the proposed method. In Fig. 4b, the PSD of the processed EMGdi is almost the same as the simulated pure EMGdi. Fig. 4c–f shows the PSD of EMGdi signal processed by other methods: Fig. 4c, the PSD of the signal processed by the hard threshold is almost the same as the PSD of the simulated corrupted signal, which shows that the hard threshold is inappropriate for removing ECG;

Table 1 Statistical results of different methods. ‘‘Inverse” hard thresholding (%)

High-pass filter (%)

Gating technique (%)

Wavelet-based adaptive filter (%)

9085 3.75e + 5

39.3 15.3

91 5.93

7.08 0.508

1.92 0.158

mv

mv

Mean Variance

Hard thresholding (%)

600 400 200 0 −200 −400 600 400 200 0 −200 −400

(a)

0 −100 5

10

15

(b)

5

10

15

(e)

100 0 −100

5

10

15

(c)

100 mv

(d)

100

5

0

−100

−100 5

10 Time (sec)

15

15

(f)

100

0

10

5

10 Time (sec)

15

Fig. 5. A clinical example: (a) the corrupted EMGdi signal; (b) the processed EMGdi signal using the hard thresholding method; (c) the processed EMGdi signal using ‘‘inverse” hard thresholding method; (d) the processed EMGdi signal using high-pass filter; (e) the processed EMGdi signal using gating technique; and (f) the processed EMGdi signal using the wavelet-based adaptive filter.

547

C. Zhan et al. / Journal of Electromyography and Kinesiology 20 (2010) 542–549

posed method is robust with these wavelet bases. However, the db4 wavelet produces the best performance: the average value of d is 1.92% with a variance of 0.158%, which is very small. The statistical results of different methods are shown in Table 1.

Fig. 4d, the PSD of signal processed by ‘‘inverse” hard thresholding is much better, but still some power of ECG is not removed; Fig. 4e, the PSD of signal processed by high-pass filter is almost the same as the PSD of simulated pure EMGdi when frequency larger than 100 Hz, but between 0 and 75 Hz, the power of EMGdi signal is not preserved; Fig. 4f, the PSD of a signal processed by the gating technique is similar to the PSD of simulated pure EMGdi in shape but smaller in amplitude. To evaluate the difference between power spectra densities of the processed EMGdi signal and simulated pure EMGdi signal, we define a relative error:

P d¼

n ðP½n

P

2 e  P½nÞ

n P½n

2

 100%;

3.2. Application to clinical data Now, the performance of the adaptive truncated thresholding approach is evaluated using clinical data. The signals were processed by a Pentium Dual Core computer (2.13 GHz  2). Fig. 5a shows one segment of the raw clinical EMGdi signal with 36,000 data and Fig. 5b–f shows the denoising results using different methods. One can see that the proposed method gives the best performance (Fig. 5f). Clinical EMGdi signals of different patients were processed and the PSD of the clinical EMGdi before and after the process by various methods are presented. Here, the processed EMGdi results of patient A and B are provided in Figs. 6 and 7, respectively. Note that the maximum amplitude of the PSD of patient A decreases from almost 120 mv2 =Hz (before processing, Fig. 6a) to about 2:5 mv2 =Hz (after processing by the proposed method, Fig. 6b), which has

ð13Þ

e where n is the frequency, P½n and P½n are the spectral densities of the simulated pure EMGdi signal and the processed EMGdi respectively. We applied our algorithm on the corrupted simulation EMGdi signal with different wavelet bases, such as Daubechies wavelet, Meyer wavelet, Symlets wavelet, and etc. Our extensive tests show that with different wavelet bases, the performances are almost the same: the average value of d is around 1:9—3:0%, namely, the pro-

120

Processed EMGdi Clinical EMGdi

100 4

(b)

80

0

2

20

0

0

2 1

3

0

4

100

200

0

0

4

(e)

3

3

2

2

1

1

100

200

(f)

2

1

(d)

3

50

60

40

4

(c)

100

PSD( mv /Hz)

PSD( mv 2/Hz)

PSD( mv 2/Hz)

(a)

0

0

50

50

100

150

100 150 Frequency (Hz)

200

250

200

0

250

0

100 200 Frequency (Hz)

0

0

100 200 Frequency (Hz)

Fig. 6. PSD of EMGdi of Patient-A processed by different approaches: (a) PSD of signal processed by the wavelet-based adaptive filter; (b) displays the same data as (a) on the expanded scale; (c) hard thresholding method; (d) ‘‘inverse” hard thresholding method; (e) high-pass filter; and (f) gating technique.

300

Processed EMGdi Clinical EMGdi

250 4

200

(b)

0

PSD( mv /Hz)

1 50

100

200

(e)

0

0

4

3

3

2

2

0

50

100

150

200

100 150 Frequency (Hz)

200

250

1

1

0 0

0

4

2

0

2 1

150

0

(d)

3

50

3

100

4

(c)

100

100

200

(f)

2

PSD( mv 2/Hz)

PSD( mv 2/Hz)

(a)

50

250

0

100 200 Frequency (Hz)

0

0

100 200 Frequency (Hz)

Fig. 7. PSD of EMGdi of Patient-B processed by different approaches: (a) PSD of the signal processed by the wavelet-based adaptive filter; (b) displays the same data as (a) on the expanded scale; (c) hard thresholding method; (d) ‘‘inverse” hard thresholding method; (e) high-pass filter; and (f) gating technique.

548

C. Zhan et al. / Journal of Electromyography and Kinesiology 20 (2010) 542–549

Table 2 Comparisons of PSD feature. Corrupted EMGdi

Hard thresholding

‘‘Inverse” hard thresholding

High-pass filter

Gating technique

Wavelet-based adaptive filter

H/L ratio PWR fc ðHzÞ

0.096 2052 33.9

Patient A 0.0467 1714 26.8

0.398 137.6 79.0

38856 86.6 155.6

1.10 160.8 87.4

1.170 224.8 88.1

H/L ratio PWR fc ðHzÞ

0.0202 3910 23.8

Patient B 0.0113 3680 21.8

0.1232 92.9 68.2

16128 34.4 222

0.170 100 57.1

0.185 137 57.2

shown most of the power of ECG signal is removed. Similar processing results are also achieved for Patient-B (shown in Fig. 7). In practice, as we can not directly record the pure clinical EMGdi signal, using relative error to evaluate the ECG cancelation performance is difficult. Various features are used to characterize the EMGdi PSD (Levine et al., 1986; Roberts et al., 1996; Schweitzer et al., 1979), such as total power, normalized power, median frequency, standard deviation and etc. In this manuscript, we have only utilized three commonly used features to analyze the processed results: high-to-low ratio (H/L ratio), total power, centroid frequency. The H/L is the ratio of power in the frequency band of 125–150 Hz to that in the frequency band of 25–50 Hz. Total power (PWR) is defined as follows:

PWR ¼

X

P½n:

ð14Þ

n

The work is supported by CERG Grant 9041163 (CityU 122906) and 9041041 (CityU 122305). Appendix A. A simple example of discrete wavelet transform For signal series s ¼ ½ 1

2 3 4 4 3 2 1  with an 8 ele pffiffiffi pffiffiffi using db4 wavelet basis, h½1 ¼ 1 þ 3 =4 2,    pffiffiffi pffiffiffi pffiffiffi pffiffiffi pffiffiffi h½2 ¼ 3 þ 3 =4 2, h½3 ¼ 3  3 =4 2, h½4 ¼ 1  3 pffiffiffi =4 2, h½n ¼ 0 for n – 0; 1; 2; 3, and g½1 ¼ h½4, g½2 ¼ h½3, g½3 ¼ h½2, g½4 ¼ h½1, then the one-level db4 forward transform matrix for an 8 element signal is ment,

2

The centroid frequency fc is

P n n  P½n fc ¼ P : n P½n

Acknowledgement

ð15Þ

The statistical results of the corrupted EMGdi signal and the processed EMGdi by different methods are shown in Table 2. For Patient-A, when comparing the H/L ratio of different methods, we can see that high-pass filter, gating technique and waveletbased adaptive filters are able to remove most of the ECG artifacts. However, the centroid frequency of EMGdi processed by high-pass filter is 155 Hz, which shows that the EMGdi power located in low frequency is also removed. Gating technique shows a good performance in H/L ratio and centroid frequency. However, the PWR is 160 mv2 =Hz, which is much smaller than the PWR of the proposed method ð224 mv2 =HzÞ. Therefore, much of the EMGdi information is missing when using the gating technique. Combining the simulation results, it has been demonstrated that the proposed method is efficient in removing ECG artifacts from the corrupted clinical EMGdi signal with little distortion. Similar processing results are also achieved for Patient-B. 4. Conclusion and discussion In this paper, a wavelet-based adaptive filter was proposed for removing ECG artifacts from EMGdi signals. Compared to other adaptive filters (Akkiraju and Reddy, 1992; Deng et al., 1998, 2000; Marque et al., 2005), this method does not require extra channels for recording ECG reference. The experimental results show that the proposed method needs no more than 35 min to process 1-h-length clinical EMGdi signal, which is efficient. Thus, the proposed method is efficient in computation and is feasible in on-line ECG removal applications. Furthermore, a quantitative comparison to other methods is given; power spectrum analysis shows that the average relative error between the PSD of simulated pure EMGdi and the corrupted EMGdi signal processed by the proposed method is within 1.92%, which is small and can be considered negligible in practice.

3

2 h½1 6 a1 ½2 7 6 7 60 6 6 1 7 6 6 a ½3 7 6 0 7 6 6 a1 ½4 7 6 7 60 6 6 1 7¼6 6 d ½1 7 6 7 6 g½1 6 6 1 7 6 6 d ½2 7 6 7 60 6 6 1 7 6 4 d ½3 5 4 0 a1 ½1

1

d ½4

0

h½2 h½3 h½4 0

0

0

0

3 2

s½1

3

7 6 s½2 7 7 7 6 7 7 6 0 0 0 h½1 h½2 h½3 h½4 7 6 s½3 7 7 7 6 7 6 0 0 0 0 0 h½1 h½2 7 7 6 s½4 7 7: 76 7 6 g½2 g½3 g½4 0 0 0 0 7 7 6 s½5 7 7 6 0 g½1 g½2 g½3 g½4 0 0 7 7 6 s½6 7 7 7 6 0 g½1 g½2 g½3 g½4 5 4 s½7 5 0 0 s½8 0 0 0 0 0 g½1 g½2 0

h½1 h½2 h½3 h½4 0

0

ð16Þ References Akkiraju P, Reddy DC. Adaptive cancellation technique in processing myoelectric activity of respiratory muscles. IEEE Trans Biomed Eng 1992;39(6):652–5. Carre P, Leman H, Fernandez C, Marque C. Denoising of the uterine EHG by an undercimated wavelet transform. IEEE Trans Biomed Eng 1998;45(9):1104–13. Deng Y, Wolf W, Bullemer F. Time frequency analysis of diaphragmatic electromyogram for the detection of diaphragm fatigue. In: Proceeding of the 20th Annual International Conference of the IEEE Engineering and Biology Society. vol. 20, No. 6; 1998. p. 3010–27. Deng Y, Wolf W, Schnell R, Appel U. New aspects to event-synchronous cancellation of ECG interference: an application of the method in diaphragmatic EMG signals. IEEE Trans Biomed Eng 2000;47(9):1177–84. Donoho D. Ideal spatial adaptation via wavelet shrinkage. Biometrika 1994;81(3):425–55. Donoho D. De-noising by soft-thresholding. IEEE Trans Inform Theory 1995;41(3): 612–27. Drake JDM, Callaghan JP. Elimination of electrocardiogram contamination from electromyogram signals: an evaluation of currently used removal techniques. J Electromyogr Kinesiol 2006;16(2):175–87. Duiverman ML1, Eykern van LA, Vennik PW, Koeter GH, Maarsingh EJ, Wijkstra PJ. Reproducibility and responsiveness of a non-invasive EMG technique of the respiratory muscles in COPD patients and in healthy subjects. J Appl Physiol 2004;96(5):1723–9. Gross D, Grassino A, Ross WRD, Macklem PT. Electromyogram pattern of diaphragmatic fatigue. J Appl Physiol 1979;46(1):1–7. Grossman A, Morlet J. Decomposition of hardy function into square integrable wavelets of constant shape. SIAM J Math Anal 1984;15(4):723–36. Johnstone I, Silverman B. Wavelet threshold estimators for data with correlated noise. J Royal Statist Soc 1997;59(2):319–51. Karthik R. ECG simulation using Matlab. .

C. Zhan et al. / Journal of Electromyography and Kinesiology 20 (2010) 542–549 Krim H, Tucker D, Mallat S, Donoho D. On denoising and best signal representation. IEEE Trans Inform Theory 1999;45(7):2225–38. Levine S, Gillen M. Diaphragmatic pressure waveform can predict electromyographic signs of diaphragmatic fatigue. J Appl Physiol 1987;64(5):1681–9. Levine S, Gillen J, Weiser P, Gillen M, Kwatny E. Description and validation of an ECG removal procedure for EMGdi power spectrum analysis. J Appl Physiol 1986;60(3):1073–81. Liang HL, Lin ZY, Wang H. Removal of ECG contamination from diaphragmatic EMG by nonlinear filtering. Nonlinear Anal 2005;63:745–53. Mallat S. A wavelet tour of signal processing. New York: Academic Press; 1998. Marque C, Bisch C, Dantas R, Elayoubi S, Brosse V, Perot C. Adaptive filtering for ECG rejection from surface EMG recordings. J Electromyogr Kinesiol 2005;15(3): 310–5. Panagiotacopulos ND, Lee JS, Pope MH, Friesen K. Evaluation of EMG signals from rehabilitation patients with lower back pain using wavelets. J Electromyogr Kinesiol 1998;8(4):269–78. Pan Q, Zhang L, Dai G, Zhang H. Two denoising methods by wavelet transform. IEEE Trans Signal Proc 1999;47(12):3401–6. Roberts C, Bartolo A, Dzwonczyk RR, Goldman E. Analysis of diaphragm EMG signals: comparison of gating vs. subtraction for removal of ECG contamination. J Appl Physiol 1996;80(6):1892–902. Schweitzer TW, Fitzgerald JW, Bowden JA, Lynne-Davies P. Spectral analysis of human inspiratory diaphragmatic electromyograms. J Appl Physiol 1979;46(1): 152–65. Stashuk DW. Simulation of electromyographic signals. J Electromyogr Kinesiol 1993;3:157–73.

Zhan Choujun received his B.S. degree in automatic control engineering from Sun Yat-Sen University, Guangzhou, China in 2007 and now a Ph.D student in the Electric Engineering department of the City University of Hong Kong. His research interests include systems biology and nonlinear dynamical systems.

549

Dr. Yeung received the B.Sc in Electronic and Electrical Engineering from Portsmouth University, UK in 1978; and Ph.D in Control Engineering from Imperial College of Science and Technology and Medicine in 1990. He joined the Department of Electronic Engineering, City University of Hong Kong in 1992; he is now Associate Professor and Director of the Control Systems laboratory, Director of the Underwater System Laboratory, member of the Hoi Ha Wan Marine Park Management Committee and Laboratory. He is now the convener of the IET Pearl River Delta Development Committee. His current research interests are bioinformatics systems, communications, fuzzy set systems, evolutionary computing and its application to system optimization. He is the co-founder the Peal Technology Co. Ltd., under CityU Enterprise (specialized in underwater robotics and systems).

Zhi Yang is a professor at Department of Electronics and Communication Engineering, Sun Yat-Sen University, China. He received his B.S. and M.S degrees in automatic control from Gansu University of Technology, China, in 1983 and 1988, respectively. His research interests include automatic instruments, CIMS, industrial automation, and complex systems modeling.