IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 42, NO.

1, JANUARY 1995

21

Detection of ECG Characteristic Points Using Wavelet Transforms

Cuiwei Li, Chongxun Zheng, and Changfeng Tai
Abstract-An algorithm based on wavelet transforms (WTs) has been developed for detecting ECG characteristic points. With the multiscale feature of WTs, the QRS complex can be distinguished from high P or T waves, noise, baseline drift, and artifacts. The relation between the characteristic points of ECG signal and those of modulus maximum pairs of its WTs is illustrated. By using this method, the detection rate of QRS complexes is above 99.8% for the MIT/BIH database and the P and T waves can also be detected, even with serious baseline drift and noise.

I. INTRODUCTION

HE AUTOMATIC detection of ECG waves is important to cardiac disease diagnosis. A good performance of an automatic ECG analyzing system depends heavily upon the accurate and reliable detection of the QRS complex, as well as the T and P waves. The detection of the QRS complex is the most important task in automatic ECG signal analysis. Once the QRS complex has been identified, a more detailed examination of ECG signal, including the heart rate, the ST segment, etc., can be performed. The algorithms for QRS detectors can generally 11. THEORY be divided into three categories: 1) nonsyntactic, 2) syntactic, and 3) hybrid. The algorithms based on a syntactic approach are time-consuming, due to the need for grammar inference A. Wavelet Transform for each class of pattems [15]. So, most of the applicable The WT of a signal f(x) is defined as QRS detectors are nonsyntactic [l], [7], [9], [lo]. Generally, these detectors first filter the ECG signal with a bandpass filter (or a matched filter) to suppress the P and T waves and noise. Then, the signal will be passed through a nonlinear where s is scale factor. Q s ( x )= $e(:)is the dilation of a transformation, for example, derivative and square, etc., to basic wavelet ) I. ( by the scale factor s. Let s = 2j (J Z, Z enhance the QRS complexes. Finally, decision rules are used to is the integral set), then the WT is called dyadic WT [3]. The determine whether QRS complexes are present in the signal. dyadic WT of a digital signal f ( n ) can be calculated with These techniques mainly suffer from two problems: 1) the Mallat algorithm [ 111 as follows: signal frequency band of the QRS complex is different for different subjects and even for different beats of the same f(n) = h k S 2 , - 1 f ( n - 2-k) (2) kEZ subject; and 2) the frequency bands of the noise and QRS complex overlap. Although a matched filter [l] can improve (3) w2, = g k S 2 - I - l . f (n - 2-llC 1 f(n) the signal-to-noise ratio, its effect is limited by the variability kEZ of QRS waveforms for different beats of the same subject. The detection of the P wave is difficult because this wave where S2, is a smoothing operator and S 2 0 f ( n ) = cl,. d, is the digital signal to be analyzed [13], which is the ECG is small and sometimes is embedded in noise. One method signal used in this paper. W 2 , f(n)is the WT of digital signal to detect the P wave has been described by Jenkins [4], in f(n). { h k 1 IC E 2) and { g k 1 k E 2) are coefficients of a Manuscript received December 27, 1993; revised September 8, 1994. This lowpass filter H ( w ) and a highpass filter G(w),respectively work was supported by the Foundation of the State Education Commission [ 131; that means of China under grant No. 9269813.
The authors are with the Biomedical Engineering Institute of Xian Jiaotong University, Xian, Shaanxi, 710049, P. R. China. IEEE Log Number 9406725. 0018-9294/95$04.00 0 1995 IEEE
kEZ kEZ

which an esophageal electrode is used to get high-amplitude P waves. This method can not be widely applied due to the uncomfortable sensation caused by the esophageal electrode and lead wire. Other algorithms for P wave detection, including the syntactic method [5] and the hidden Markov methold [6], are complex and time-consuming. Gritzali proposed a simple method to detect P and T waves by length transformation [7], but it is not robust to noise. Wavelet transform is a very promising technique for timefrequency analysis. By decomposing signals into elementary building blocks that are well localized both in time and frequency, the WT can characterize the local regularity of signals [3]. This feature can be used to distinguish ECG waves from serious noise, artifacts and baseline drift. In this paper, an algorithm based on the WT for detecting QRS complex, .P and T waves is proposed. A dyadic wavelet transform is used for extracting ECG characteristic points. The local maxima of the WT modulus at different scales can be used to locate the sharp variation points of ECG signals. The algorithm first detects the QRS complex, then the T wave, and finally the P wave.

22

IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 42, NO. 1, JANUARY 1995

The wavelet we used is a quadratic spline wavelet with compact support and one vanishing moment. It is a first derivative of a smooth function [13]. The Fourier transform of G ( x ) is
4

n.

The symbol A represents the discrete Fourier transform. The filters H ( w ) and G ( w ) are

B. Equivalent Filter of WT
The discrete Fourier transform of WT is [13]
W 2 )f ( w )

=j ( w ) 4 ( 2 j w )
G(W)h)&$

(a)

(b)

(C)

(d)

(e)

j=1

j=2 G(24fWf(w)!&) G(2j-'w)H(2j-'w). . . H ( w ) f ( w ) $ ( w ) j > 2 ={ (6)

where 4 is a smooth function, and f^(w)&w) is the discrete Fourier transform of the ECG signal d, ( n E 2 ) [13]. From (6), the WT of f ( n ) at scale 2j is equal to a filtered signal of d, that passed through a digital bandpass filter (or high pass filter for scale 2'). Let Q j ( w ) be the transform function of the equivalent filter. Then, we have
Qj(w) =

Fig. 1. The relation between the characteristic points of simulation waves and those of their WT's at different scales. The uniphase wave (a) is symmetric and T, = 24-1 - 1 = 7 points. The uniphase wave (b) is not symmetric and Tb - T, = 2 points. Complex (c) is symmetric and T, = T, = 7 points. h waves (d) and (e),Tdl < re,,d 2 < Te2.Tel = 24-1 - 1 = 7 points and T Te2 = 24 24-' - 2 = 22 points.

j=1

j =2 z/;:)H(w) G(2j-'w)H(2j-'w) . . . H ( w ) j > 2.

From (4), and (7), we can deduce (see Appendix) (3,

2Jf2J-1-2
&'(U)

83-1

k=1-2J-1

qieikw

(131-2J-lfk

= -q;2+23-1-2-k IC E [l- 23-1,2j

'3

+ 2j-1 - 21

maximum-negative minimum pair with a delay of about 2j-l1 points. The larger the scale, the bigger is the error between the delay and 2j-l- 1. At small scale 2 l , the error is zero. At large scale, the error is still tolerant. For example, in Fig. l(b), the slope of rising edge is three times larger than that of the dropping edge, but the error is only two points at ~ c a l e . 2 ~ . From (8), we can see that the beginning point of q: is at (7) 1 -2jP1, and the ending point of qh is at 2j 2j-l- 2. So, the onset of a uniphase wave corresponds to the onset of a positive maximum-negative minimum pair with a delay of 1-23-l points, and the offset corresponds to an offset with a delay of 2j 2j-l- 2 points, as shown in Fig. l(e). (8) Corresponding to 250/s ECG sampling rate, the amplitudefrequency responses of the equivalent digital filters Q j ( w ) (j= 1, 2, . . . , 5 ) are illustrated in Fig. 2, and their 3-dB bandwidths are listed in Table I.

where Q j ( w ) is an FIR digital filter with generalized linear phase [2]. The filter is antisymmetric and the time delay of its central point is (the delay is considered as 2j-I- 1 points in the algorithm). According to the form of the quadratic spline wavelet [ l 11, we know that every uniphase wave, like the (a) or (b) waves in Fig. 1, corresponds to a positive maximumnegative minimum pair of WZ3f(n) different scales. The at wave rising edge corresponds to a negative minimum, and the dropping edge corresponds to a positive maximum at different scales. The moduli of these maxima or minima corresponding to the same edge are named as the modulus maximum line. If the uniphase wave is symmetric to its peak, as shown in Fig. l(a), then its peak corresponds to the zero-crossing point of the positive maximum-negative minimum pair with a delay of exactly 2j-I- 1 points. If the uniphase wave is not symmetric to its peak, as shown in Fig. l(b), then the peak corresponds to the zero-crossing point of the positive

C. Relation Between the Signal Singularity and its W S

Signal singularities often carry the most important information. It is important to find the location of singularity and characterize the singular degree in signal processing. The singular degree is often described with the Lipschitz exponent [ l l ] . On the pioneer work of Grossman [12], Mallat founded the relation between the WT's and singularities of signals [ 111. If wavelet $(x) is the first derivative of a smooth function, then 1 1) When scale s is small enough, the maxima of 1 W,f (x) indicate the locations of sharp variation signal points.

2 ) The function f(z)is Lipschitz cy (0 < Q < 1) over [a, b ] , if and only if there exists a constant A such that for all x E [ a , b ] ,
I W 2 J f ( Z ) I5

(9)

LI et al.: DETECTION OF ECG CHARACTERISTIC POINTS USING WAVELET TRANSFORMS

23

3, the energy of the QRS complex is decreased further and,

at the same time, the energies of motion artifact and noise are increased. Furthermore, the selection of more scales requires more cdculations. So, we only select scales from 2l to 24. These scales are called characteristic scales in this paper. 2 ) Determination of Modulus Maximum Lines o R Waves: f From Fig. l(a), we can see that the wave corresponds to two modulus maximum lines with opposite signs of W,, f ( n ) . The modulus maximum lines corresponding to R waves at characteristic scales must be determined. A method to select these modulus maxima from large to small scale is as follows: Step 1: Find all of the modulus maxima larger than a threshold 4 at scale 24 to obtain the location set of modulus maxima {n; I k = 1 . . . N } ; Step 2 : Find a modulus maximum larger than the threshold ~3 on the neighborhood of n: at scale 23 and define its location as n:. If several modulus maxima exist, then the largest one is selected. But, the modulus maximum with its " L f (Hz) 25 50 75 100 125 location nearest nf will be selected if the largest one is not larger than 1.2 times the others. If no modulus maximum Fig. 2. The amplitude-frequency responses of equivalent filter Q J( U ) at exists, then set n:, n:, and n: to zero. So, the location set different scales Corresponding to 250/s sampling rate. f = 125w/n. {n: 1 k = 1 . . . N } can be found. Step 3: Similar to Step 2 , the locations of the modulus TABLE I maxima at scales 22 and 2l are found. Then, the location 3-dB BANDWIDTHS EQUIVALENT OF FILTER &'[U) AT D~FFERENT SCALES sets {n: 1 k = 1 . . . N } and {nk 1 k = l . . . N } will be obtained. scale 3 db bandwidth (Hz) By searching modulus maxima at characteristic scales, the s=2I 62.5 w 125.0 location set of modulus maximum lines {nf , n:, n:, nk I k = 1. . . N } is obtained. By eliminating those modulus maximum s=22 18 w 58.5 lines with n; = 0, the location set of remaining modulus s=23 a 27 maximum lines is {n;,ni,n:,nk I k = l...Nl}. Searching the modulus maxima from large to small scale s=24 4 13.5 has two reasons. First, the number of modulus maxima at large scale are much fewer than those at small scale, so searching the s=25 2 6.5 modulus maxima from large to small scale can save calculating time. Second, high frequency noise decays greatly at large scale so that it creates no modulus maximum or a very small 111. DETECTION METHODS one. Therefore, by using an appropriate threshold to detect modulus maxima at large scale, the faulty detections caused A . R Peak Detection by high frequency noise can be reduced greatly. The modulus We analyzed 600 data points of an ECG signal each time. maxima detected at large scale may also include those induced The WT's of ECG signals are calculated with (2) and (3), then by P and T waves, as well as the baseline drift. But their the decision rules are applied for detection. A more detailed effects at small scale are slight because their high frequency description of the detection procedure is given as follows: energies are trifling. Using a threshold for these modulus I ) Selection o Characteristic Scales: From Fig. 2 , it can maxima at small scale can also eliminate their effects. f be seen that the W,,f(n> at small scales reflects the high 3 ) Calculation of Singular Degree: Let a j ( n k ) = IW,, frequency components of the signal and, at large scales, reflects f ( n k ) l and suppose a: is the upper bound of the Lipschitz the low frequency components of the signal. According to the exponent. a: is called the regularity exponent [3]. From (9), power spectra of the ECG signal, noise and artifact [8], and we can deduce the regularity exponent [14] the 3-dB bandwidths of Q j ( w ) shown in Table I, it is clear (10) ( R log, aj+l(nk) - log, aj(nk). 2 5 that most energies of the QRS complex are at the scales of Z3 and 24, and the energy at scale 23 is the largest. From scale 23 to smaller or larger scales, the energy of the QRS complex Let decreases gradually. According to our experiments, for QRS " j = log,, a j + l ( n ; + l ) - log, a j ( n i ) complex with more high frequency components, the energy ' at scale 2 is larger than that at scale Z3, and for the QRS where a j represents the decay of IW,, f ( n k ) l over the scales complex with more low frequency components, the energy at and is an approximation of the Lipschitz exponent at singularscale 24 is larger than that at scale 23. For larger scales 2j (j2 ity nk. If characteristic scales are set from 2 to 2 J f N ,then '
N

24

IEEE TRAlUSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 42, NO. 1, JANUARY 1995

a j can be calculated from j = J to j = J N - 1. When the scale is large enough, a j ( n k ) will decrease to zero, no matter the regularity of the signal at n k [3]. In fact, when j + +CO, then [3]

Next, eliminate redundant modulus maximum lines. Usually, the R wave corresponds to only two modulus maximum lines. But, for some dual R waves or noise in the neighborhood (120 ms) of a modulus maximum line, there are two or more modulus maximum lines, only one of which is useful; the others are redundant. The redundant modulus maximum lines are eliminated by the rules described below. Because the main energy of the QRS complex is at scale 23, Hence, from (1 1), we can deduce aj 5 - 1 if j is big enough. the modulux maxima at this scale is analyzed. For example, So, small scales should be selected [3], in order to avoid in the neighborhood (120 ms) of a positive maximum, there aj 5 - 1. Therefore, a = will be a better are two negative minimum. Let the two minima be Minl and Min2, their absolute values be AI and A2, respectively, approximation of the regularity exponent a. At the characteristic scales from 2l to Z4, we can calculate and the intervals between the minima and maximum beL1 the a1, az, and a3 of a singularity point. The R wave always and L2, respectively. The rules of judging redundant modulus corresponds to a1 > 0, and mostly a2 > 0. Although maximum lines are as follows: some modulus maxima of R waves with higher frequency Rule 1: If A1/L1 > 1.2Az/L2, Min2 is redundant. components decay fast at large scale to make a2 5 0, a1 a2 Rule 2 : If A2/Lz > 1.2A1/L1, Minl is redundant. is still greater than zero. For most R waves, their energies at Rule 3: Otherwise, if Minl and Min2 are on the same side of scale 23 are larger than those at scale 24, and the decay of the positive maximum, then the minimum farther from the I W 2 3 f ( n k ) l from Z3 to 24 is large enough to make not only maximum is redundant. If Minl and Min2 are on different sides of the maximum, then the minimum following the a3 < 0, but also a1 a2 a3 5 0. For high frequency noise and interference with sharp irregularities, there are also maximum is redundant. With these judging rules, the R wave edge with the larger a1 5 0, a2 5 0, and a3 5 0, hence al+az+a3 5 0. So, from the value of a1 a2 a3, the R wave, high frequency noise, slope and amplitude, or the positive R wave (corresponding and interference can not be distinguished. Therefore, only a1 to positive maximum after negative minimum in a modulus and a2 are selected and we let a = 1:2 in order to make maximum pair) is selected. a > 0 for most of the normal R waves. For distorted R waves, A similar procedure is implemented for the case of one negthe increasing of the low frequency components can only make ative minimum and two positive maxima, and the redundant a much larger. So, if the a suddenly decreases greatly or even modulus maximum lines can also be eliminated. becomes negative, the corresponding singularity point must be 5 ) Detection o the R Peak: According to the relation bef noise or interference, which will be eliminated. tween signal characteristic points and those of WT depicted 4 ) Elimination o Isolation and Redundant Modulus Max- in Section 11, R peak can be located at a zero-crossing point f imum Lines: The frequency bands of motion artifact and of a positive maximum-negative minimum pair at scale 2. muscle noise often overlap with that of the QRS complex. After we eliminate the isolation and redundant lines from the Ordinary bandpass filters can not eliminate or even decrease location set of modulus maximum lines, n:(k = 1 , 2 . . . N2) their effects. In the location set of modulus maximum lines, in the remaining set is only composed of the locations of the there may be those created by motion artifact and muscle noise. positive maximum-negative minimum pairs at scale 2. Thus, By eliminating isolation and redundant modulus maximum the zero-crossing points of these positive maximum-negative lines, we can greatly reduce the effects of motion artifact and minimum pairs are found to obtain the locations of R peaks. muscle noise. In R peak detection, four thresholds { 1, 2, 3, 4) are First, eliminate isolation modulus maximum lines. The R used for the modulus maxima detection at four different wave corresponds to a positive maximum-negative minimum scales. These thresholds are updated by the modulus maximum pair at each characteristic scale and the interval of its two of WTs of the most recently recognized QRS complex as modulus maxima at scale 2l is shorter than its width, as follows: shown in Fig. l(a). Suppose ni is the location of a positive If IW23f(ni)I 2A7 then A,+1 = A?; 2 f(n)at scale 2l, and maximum (negative minimum) of W2, else A,+ = 0.875AY (IC = 1 . . . N I ,IC # 1) is the location of a negative minimum 0.1251W23f(ni)l. (positive maximum) of W 2 J f ( nat the same scale. If the ) Let el = 0.3A,+ 3 = 1,2.3,4 interval between n: and (IC # 1) is larger than an interval threshold, the maximum (minimum) at n: is considered the is isolation maximum (minimum). The corresponding modulus where IWZJf(ni)l the modulus maximum of WT corre2 maximum line is an isolation line and should be eliminated sponding to the R wave. When IWZlf(n3,)1 2A,, the from the set of modulus maximum lines. The selected interval parameter A,+ maintains its original value so that the sudden threshold should be approximately the same as the interval of increase of the QRS amplitude would not effect the next the two modulus maxima created by the widest possible QRS judging. complex in order that wide QRS complexes are not lost and In most cases, the algorithm correctly detects the R peak. artifacts and noise are mostly eliminated. In this paper, the However, some tactics based on pioneer work [lo] are still interval threshold is empirically defined as 120 ms. used to improve the detection.

aJ+.z+~v-l
+

+ +

+ +

ni

ni

LI et al.: DETECTION OF ECG CHARACXERISTIC POINTS USING WAVELET TRANSFORMS

25

',

'I

(a)

(b)

(C)

Fig. 3. ECG signals and their WT's at scale 2l to Z 5 . The modulus maxima of WT of the QRS complexes at scale 23 are largest. The effects of high frequency noise are mainly at small scales, as shown in (a). The effects of baseline drift are mainly at large scales, as shown in (h). An artifact (c) creates an isolation modulus maximum in a certain time.

Tactic 1: Blanking, where events immediately following a QRS detection are ignored for a period of 200 ms. Tactic 2: Searching back, where previously rejected events are reevaluated when a significant time has passed without finding a QRS complex. If no QRS complex was detected within 150% of the latest average RR interval, then we detect the modulus maxima again at scale 23 with a threshold of 0.5 3. If a positive maximum and negative minimum of WT are detected and the interval between them is less than 120 ms, then their zero-crossing point is detected and the reevaluated R peak can be found with a modification of 3 (= 23-1- 1) points delay. QRS Onset and OfSset Detection
After the detection of R peaks, the onset and offset of the QRS complex are also detected. The onset of the QRS complex is defined as the beginning of the Q wave (or R wave when Q wave is not present), and the offset of the QRS complex is defined as the ending of the S wave (or R wave when the S wave is not present). Ordinarily, the Q and S waves are high frequency and low amplitude waves and their energies are mainly at small scale. So, we detect the onsets and offsets of QRS complexes at scale 2 l . From the (c) wave in Fig. 1, we can see that the onset of QRS corresponds to the beginning of the first modulus maximum before the modulus maximum pair created by the R wave, and the offset of QRS corresponds to the ending of the first modulus maximum after that modulus maximum pair. From the modulus maximum pair of the R wave, the beginning and ending of the first modulus maxima before and after the modulus maximum pair are detected within a time window. The reason for detecting the beginning and ending points at scale 2 l , rather than at the original signal, is to avoid the effect of baseline drift.

the power spectra of ECG signal, noise and artifact [8] and the passbands of filter Q j ( w ) at different scales in Table I, the energies of T and P waves are mainly at scales z4 and 25. But, baseline drift is serious at scale 25, so scale 24 is selected to detect T and P waves. The T wave creates a pair of modulus maxima with a different sign of WzJf(n) at scale 24, within a time window after the detected R peak. Because the T wave is almost symmetric to its peak, the peak of the T wave corresponds to the zero-crossing point of the modulus maximum pair with a 7 (24-1 - 1) point delay. The onset and offset of the T wave correspond to the onset and offset of the modulus maxima pair with delays of -7 ( 1 -24-1) and 22 (24 24-1- 2) points, respectively. In practice, the points near the onset and offset of the modulus maximum pair are also approximately zero because the T wave is smoothed by the equivalent filter Q j ( U ) at large scale, so those actually detected are inside the interval between the onset and offset of the modulus maximum pair, as shown in Fig. l(d). When we detect the onset and offset of the modulus maxima pair, the onset is actually modified with a delay of -2 points and the offset is 17 points. These practice delays are empirical. The peak, onset, and offset of the P wave are detected similarly to those of the T wave within a time window before the detected R wave.

IV. RESULTS
We used the MIT/BIH arrhythmia database to evaluate our algorithm. Only channel 1 of the two-channel ECG signal in the database was used. The WT's of some ECG signals with serious high frequency noise, baseline drift, and artifacts are shown in Fig. 3. In Table 11, we list QRS detection rates for records of the MIT/BIH database. The WT algorithm produces 65 false positive (FP) beats (0.05%) and 112 false negative (FN) beats (0.1%) for a total detection failure of 177 beats (0.15%). Record 105 is more noisy than the others; Record 108 has unusually high and sharp P waves; Record 203 has a great number of QRS complexes with multiform ventricular

C. T and P Wave Detection

After the detection of the QRS complex, the peaks, onsets, and offsets of T and P waves are also detected. According to

26

IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 42, NO. 1, JANUARY 1995

TABLE II
RESULTSOF THE QRS DETEC~ON ALGOR~THM THE h4IT/TjM DATABASE FOR

Type

Total (beats) 2273 1865 2187 2084 2230 2572 2027 2137 1763 2532 2124 2539 1795 1879 1953 2412 1535 2275 1987 1863 2476 1518 1619 2601 1963 2136 2982 2656 1862 2956 3004 2647 2748 3251 2208 2154 2048 2427 2484 2605 2053 2256 1886 1780 3079 2753

FP
(beats) 0 1 0 0 8 15 2 0 13 0 1 2 2 3 0 0 1 1 1 2 0 0 0 0 1

FN
(beats) 0 0 0 0 2 13 3 0 15

..........................................................

Failed detection (beats) 0 1 0

0

Failed detection (2) 0

0 0.11 0

100 101 102 103 104 105 106 107 io8 109
111

10 28 5
0

0.45 1.09 0.25

0

28

1.59

0
1 1 0 0

0
2 3 2 3 0
1

112 113 114 115 116 117 118 119 121 122 123 124 200 201 202 203 205 207 208 209 210 212 213 217 219 220 221 222 223 228 230 231 232 233 234

0 0.09 0.12
0.11 0.16 0 0.04 0.07 0.04 0.05 0.16 0 0 0 0.04 0.66 0.05 0.87 0.04 0.27 0.14 0 0.23 0 0 0.09 0 0 0.29 0.40 0.08 0.49 0 0 0 0 0

0
1

0 0 0
1

0 0
0 1

1 1 1 3 0 0 0
1

0
2 0 2 0 0 3

12 1 24
1

3
4 0 3 0 0 1 0 0 7 9 2 7 0 0 0 0 0 112

0 0 1 0 0 0
1

13 1 26 1 5 4 0 6 0 0 2
0 0

0 3 0 0
0

7 10 2 10 0
0 0

0 0 65

0
0

-----_____---------_____________________------------------

TOTALS:116137

177

0.15

arrhythmia; and Record 222 has some non-QRS waves with highly unusual morphologies. These records often create more QRS detection errors with other algorithms [9], [lo], but significant improvement was achieved with our WT algorithm. Some ECG records with artifact, serious noise, and baseline drift are illustrated in Fig. 4, in which the detected characteristic points are marked.

V. CONCLUSION DISCUSSION AND

In this paper, we present an algorithm based on WT for the detection of QRS, T, and P waves of ECG. There are several distinct features to our method.

1) With multiscale information, it is easy to characterize the ECG waves, and the QRS complex is easy to distinguish from high P and T waves, noise, baseline drift, and interference. 2) A quadratic spline wavelet with compact support is used. It has generalized linear phase so there is a determinate relation between ECG characteristic points and the modulus maxima, or the zero-crossing points, of the WTs. 3) A QRS complex corresponds to a modulus maxima pair of the WT. On the contrary, an artifact having only a rising or dropping edge in a certain time corresponds

LI et al.: DETECTION OF ECG CHARACTERISTIC POINTS USING WAVELET TRANSFORMS

21

R R

Fig. 4. The detected ECG waves. The symbol R marks the locations of R peaks, where the positive R waves are marked above and the negative ones are marked below. The symbols Q and S mark the beginning and ending locations of QRS complexes. The symbols P and T mark the peak locations and widths of P and T waves, respectively. (a) and @) Different noise. (c) Baseline drift. (d) Artifacts. (e) High P waves. (f) High T waves. (8) and (h) Multiform ventricular arrhythma beats.

to only one modulus maximum of the WT. So, the algorithm can greatly reduce the effect of artifact. 4) The algorithm can detect not only QRS complex, but also the T and P waves. So, other detailed characteristics, for example, ST segment and RR rate, etc., can easily be obtained. Besides the quadratic spline wavelet, we also used higher order spline wavelets, whose Fourier transforms are as follows: $(w) = i,( sinw/4

-)

2n+1

n = 2,3,4.

The equivalent filters of the WTs using these wavelets have bandwidths approximating those of the quadratic spline wavelet, so the results of ECG detection with these wavelets are almost the same as those with quadratic spline wavelet. But, the coefficient series of their equivalent filters is longer than that of the quadratic spline wavelet, so more time is required to calculate their WTs. Daubechies orthonormal wavelets [16] are also used in our experiments. But, the orthonormal wavelet does not have linear phase or generalized linear phase [2], so no determinate

relation between the characteristic points of ECG and those of their WTs can be found. Therefore, these wavelets do not give detection results as good as the spline wavelets. From the results, we know that the onsets and offsets of P and T waves may not be detected to an accuracy required for morphological diagnosis when they are influenced seriously by noise or baseline drift or their amplitudes are too small. Also, only the uniphase P and T waves are considered in this paper. This method can be considered for further development and supplemented with other techniques for the biphase P and T waves, and this will be our further research work. The algorithm is implemented on a 486-DX33 personal computer with C Language. Processing 10 min ECG data (250/s sampling rate) requires 1 min, on average. This speed is not good enough to analyze long-time ambulatory ECG data. But, this algorithm is promising with higher speed computers or microprocessors. The algorithm for detecting ECG characteristic points based on WT shows the potential of the WT, especially for processing time-varying biomedical signals. The power of WT lies in its multiscale information analysis which can characterize a signal very well. It is clear that the WT method will lead to a new way of biomedical signal processing.

28

IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 42, NO. 1, JANUARY 1995

APPENDIX ANTISYMMETRIC CHARACTERISTIC

OF THE EQUIVALENT FILTER

From (4), (3, and (7), we have

j = 1, Q 3 ( w ) = 2(eZw- 1) 1 3 = 2, Q J ( w ) = - ( ~ 4 z w + 4
~ ~ 3 z w + ~ ~ 2 2 w - ~ ~e2 2 W ~ . - -w - )

(AI) From (4), we know 1 ~ ( 2 3 w = - ( e 2 3 + 1 z w + 3e2zw + 3 + e - 2 ~ ~ ) 6 4 2 ) )

8

so H(23w) (j 2 ) is a symmetric filter. Let E P ( w ) = H(23-2w) . . H ( w )

j

> 2.

Then, from (A2), we can deduce that P ( w ) is also a symmetric filter, where

[6] D. A. Coast, R. M. Stem, G. G. Cano, and S. A. Briller, An approach to cardiac arrhythmia analysis using hidden Markov model, IEEE Trans. Biomed. Eng., vol. 37, pp. 8 2 6 8 3 5 , 1990. [7] F. Gritzali, G. Frangakis, and G. Papakonstantinou, Detection of the P and T waves in an ECG, Compuf.Biomed. Res., vol. 22, pp. 83-91, 1989. [8] N. V. Thakor, J. G. Webster, and W. J. Tompkins, Estimation of QRS complex power spectra for design of a QRS filter, IEEE Trans. Biomed. Eng., vol. BME-31, pp. 702-705, 1984. [9] J. Pan and W. J. Tompkins, A real-time QRS detection algorithm, IEEE Trans. Biomed. Eng., vol. BME-32, pp. 230-236, 1985. [lo] P. S. Hamilton and W. J. Tompkins, Quantitative investigation of QRS I BH detection rules using the MT I arrhythmia database, IEEE Trans. Biomed. Eng., vol. BME-33, pp. 1157-1187, 1986. [ l 11 S. Mallat, Zero-crossings of a wavelet transform, IEEE Trans. Inform. Theory, vol. 37, pp. 1019-1033, 1991. [12] A. Grossman, Wavelet transform and edge detection, in Stochastic Processes in Physics and Engineering, M. Hazewinkel, Ed. Dodrecht: Reidel, 1986. [I31 S. Mallat, Characterization of signals from multiscale edges, IEEE Trans. Pattern Anal. Machine Infell., vol. 14, pp. 71G732, 1992. [I41 C. Z. Qiang, A study of fault diagnosis and cutting chatter with wavelets, Chinese J.Huazhong Univ. Sei. and Tech., vol. 21, pp. 88-94, 1993. [I51 E. Pietka, Feature extraction in computerized approach to the ECG analysis, Pattern Recog., vol. 24, pp. 139-146, 1991. [ 161 I. Daubechies, Orthonormal bases of compactly supported wavelets, Commun. Pure Appl. Math., vol. 41, pp. 909-996, 1988.

therefore
Qj(w) = G(2j-)wP(w) j

>2

Cuiwei Li received the M.S. degree in automatic control in 1991, and the Ph.D. degree in biomedical engineering in 1994 from Xian Jiaotong University, China. She is now with the Biomedical Engineering Institute of Xian Jiaotong University, where she is working on biomedical signal processing.

From (A3) and (A4), we get

23+23--1-2

> 2,

&(U)

=
k=1-23-1

qhezkw

qi-ZJ-+k

= -&+23-1-2-k

IC

E [l - 23--1,23

+ 23-1

21.

From (Al) and (A5), we know the equivalent filter Q J ( w )is antisymmetric whose central point is at

Chongxun Zheng received the diploma in electrical engineering from Xian Jiaotong University, China, in 1962. He worked in the Department of Electrical Engineering, Xian Jiaotong University from 1962 to 1981. He was as a Visiting Scholar in Kings College of London and Brunel University, U.K., from 1981 to 1983. He joined the Department of Information and Control Engineering, Xian Jiaotong University from 1984 to 1990, where he was an Associate Professor, and then a Professor (since 1988). He is currently a Professor and the Director of the Research Institute of Biomedical Engineering at the same university. His research interests include biomedical signal processing and applications, cardiac electrophysiology, speech recognition, and diagnosis and reliable design of digital systems.

y.

REFERENCES
[I] Q. Z. Xie, Y. H. Hu, and W. J. Tompkins, Neural-network based adaptive matched filtering of QRS detection, IEEE Trans. Biomed. Eng., vol. 39, pp. 317-329, 1992. [2] C. K. Chui, Wavelet Analysis and ifs Applications. New York: Academic Press, 1992. [3] S. Mallat and W. L. Hwang, Singularity detection and processing with wavelets, IEEE Trans. Inform. Theory, vol. 38, pp. 617-643, 1992. [4] J. M. Jenkins, D. Wu, and R. C. Arzbaecher, Computer diagnosis of supraventricular and ventricular arrhythmias, Circ., vol. 60,pp. 977-985, 1979. [ 5 ] E. Skordalakis, Syntactic ECG processing: A review, Pattern Recog., vol. 19, pp. 305-313, 1986.

Changfeng Tai was bom in Shaanxi, China, on May 27, 1966. He received the B.S., M.S., and Ph.D. degrees in biomedical engineenng from Xian Jiaotong University, China, in 1986, 1989, and 1992, respectively. He has been a faculty member at the Biomedical Engineering Institute, Xian Jiaotong University, since 1986, where he teaches electncal neurophysiology and biomedical signal processing. His research interests include functional neuromuscular stimulation, nerve modeling, and biomedical signal processing. Dr. Tai was named the excellent graduate student of Xian Jiaotong Umversity in 1986.