- Email: [email protected]
Decomposition of ECG by linear filtering
Co,,,put. Biol. Med. Vol. 22. No. I/2. pp Printed in Great Britain
13-22, 1992 0
DECOMPOSITION
OF ECG BY LINEAR
I. S. N. Department
of Electrical
MURTHY*
Engineering,
Indian
and U. C. Institute
(Receiued 24 June 1991; in revised form 23 September November 1991)
0010.4825192 $55.00+ .oO IY92 Pergamon Press plc
FILTERING
NIRANJAN
of Science,
Bangalore
560 012, India
1991; received for publication 25
Abstract-A simple method is developed for the delineation of a given electrocardiogram (ECG) signal into its component waves. The properties of discrete cosine transform (DCT) are exploited for the purpose. The transformed signal is convolved with appropriate filters and the component waves are obtained by computing the inverse transform (IDCT) of the filtered signals. The filters are derived from the time signal itself. Analysis of continuous strips of ECG signals with various arrhythmias showed that the performance of the method is satisfactory both qualitatively and quantitatively. The small amplitude P wave usually had a high percentage rms difference (PRD) compared to the other large component waves. ECG analysis PRD DFTM
Component delineation FIR filter design
1.
Arrhythmia
Linear
filtering
DCT
INTRODUCTION
Delineation of the component waves in an ECG signal is vital for the detection of cardiac arrhythmias of critically ill patients in an intensive coronary care unit (ICCU). This topic has received considerable attention in the literature. Fairly accurate methods [l-4] are available for the detection of the QRS complex. However, the twin factors of small size and absence of spectral variability for temporally differing P waves have inhibited all the efforts for its detection [5,6]. Thus, one of the main obstacles for the development of a comprehensive arrhythmia analysis system has been the lack of simple and reliable methods for P wave detection. In the present paper a new approach based on a simple linear filtering process is proposed for the delineation of component waves in ECG irrespective of their shape, size and location. The linearity and frequency transformation properties of the discrete cosine transform (DCT) are exploited to convert each component wave into cosinusoids of a distinct frequency band in the transform domain. Furthermore the magnitude spectrum of the transformed signal and the rectified time signal look alike. As such the required filter characteristics such as lowpass (LP ), bandpass (BP) and highpass (HP) can be inferred from each one of them. Appropriate FIR filters [7] are then designed using the information available from the time signal itself. In the case of signals with well separated component waves, flexible filters can be obtained making the delineation relatively easy. For signals with arrhythmias where such clear separation between component waves does not exist, compound waves, containing two or more components are delineated. In contrast to this, the spectra of the component waves in the direct time signal overlap with each other and do not permit recovery of the components by such simple linear filtering [S]. Continuous strips of ECG signals sampled at 500 samples per second (sps) containing about 300 beats have been analyzed. The results obtained have been found to be satisfactory both qualitatively and quantitatively. In Section 2, the DCT and its linearity property are explained. The frequency transformation property of the DCT, derivation of the mathematical expressions for discrete Fourier transform magnitude (DFTM) and filter characteristics are considered in Section 3. In the next section, the scheme of analysis, examples and application of the *To
whom
correspondence
should
be addressed. 13
I. S. N. MLJRTHY and U. C. NIRANIAN
14
method to real ECG are dealt with. Finally, discussion and summary are included in Section 5. 2. TRANSFORM
AND
PROPERTIES
The DCT of a discrete time signal x(n) of length (N+ 1) samples is defined [9] as X(k) = A/@%$,
i c&n) n=O
cos (nknlN);
k= 0,1,2,
...,N
(1)
n = 0, 1,2, . . . , N.
(2)
wherecj=l/fifori=OorN;Ci=lfori#OorN. The inverse transform (IDCT) of X(k) is x(n) = t/c2,N>c,
2
ckX(k) cos (nkn/N);
k=O
Consider an ECG signal x(n), with P, QRS and T waves -x,(n), x,(n), x,(n) each of length (N + 1) samples. Each component wave exists over a small set of samples with the remaining ones as zero, i.e. x(n) = c x,(n); n =0, 1,2, . . . , N. “P.9.’
(3)
Now, if DCT 44
=
z
then
x(n) =
c
i=p4,
-c(n) f
5
2 X,(k) =X(k). i=p.
(4)
q, 1
Thus, by the linearity property, the DCT of a complete beat is expressed as the sum of the transforms of the individual components and vice versa. In other words, addition and deletion of component waves to and from the complete beat is possible. 3. THEORY 3.1. Frequency transformation We now show that the DCT converts each sample in a signal x(n) into a cosinusoid of distinct frequency and a group of consecutive samples into cosinusoids in a frequency band. This property enables component delineation by simple filtering of the transformed signal. Consider a signal x(n), of length (N + 1) samples. Let x(n) contain a single impulse at n = no with an amplitude a,,. Then we see from equation (1) that X(k) is a cosine wave with frequency (n,,zlN) and amplitude -a,,, i.e. X(k) = d@%$+,ckanO cos (n&z/N);
k = 0, 1,2, . . . , N.
(5)
When this impulse is shifted from n = no = 1 to n = no = N in steps of one sample at a time, the frequency of oscillation of X(k) changes from (n/N) to n rad in steps of (n/N) rad. Thus, from a given X(k), the location and amplitude of the impulse in x(n) can be determined using equation (2). We now extend this argument to a set of consecutive samples located between time samples n = n, and n = nl - 1. This group of samples may
Decomposition
of ECG by linear filtering
15
take any arbitrary shape of a component wave such as a P, QRS or T. Expressing this group of samples, x,(n), in terms of x(n), q-1
xl(n)=x(n)
x
6(n-m)
fornO~n~nl-l fornl
Substitution
Ocndn,-1.
of equation (6) into equation (1) leads to
X,(k) =
qzqc,[c~(no)
cos (n&lN)
+ c,+1x(no+ 1)
cos ((no + 1)knlN) + * . . + c,,_,x(q k=O,1,2,..
- 1) cos ((n1- 1) kJr/N)]
. ,N.
(7)
We note that X,(k), the DCT of xl(n) is made up of an equal number of (n,-no) cosinusoids each with a distinct frequency (nnlN) rad and amplitude ma,,; , N. It can be shown that the sum of these cosinusoids in X,(k) leads to a n=O,1,2,... damped cosine wave whose frequency of oscillation is determined by the sample number which has the largest amplitude in xl(n). We can extend this argument to an arbitrary waveform such as the ECG containing more than one component wave. From the linearity property in (4), it is not difficult to see that each component wave in x(n) gets transformed into cosinusoids in a distinct frequency band in X(k). The component waves can be separated out from each other by filtering X(k) in the DCT domain. On the other hand such a simple separation of the component waves is not possible from the direct time signal x(n). 3.2. Discrete Fourier trunsform magnitude (DFTM)
of X(k)
Consider X(n), the 2N point discrete Fourier transform sequence X(k) [9], ZN-
X(a) = x
(DFT) of the (N+ 1) point
1
X(k) e-iznkn’2N;S2 = 0, 1,2, . . . ,2N - 1
(8)
k=O
whereX(k)=OforkaN+l. Substituting X,(k) from equation (7) for X(k) in equation (8), (X,(S2)( the 2N point DFTM of X,(k) is obtained as ZN-
IX,(Q)1 =
I
C cfix,(n)(l/V@Vj) n=O
sin (n(N + 1) (n - Sr)/ZN) e+j”((“-R)‘2) sin (x(n - s2)/2N)
+ e_jn((n+Rj,2j sin (n(N+ I) (n + &)/2N) sin (n(n + 52)/2N) + (1- ti){V5
+ l/ti((
- l)‘“_*‘+ (- l)‘“+*‘)} II sZ=O,1,2 ,...,
2N-1.
(9)
Here, we note that the DFT coefficients X,(a), are the weighted sum of the constituent samples in xl(n). The weighting function inside the square brackets on the right-hand-side (r.h.s.) of equation (9) attains its maximum value in magnitude when SL= n and 52= (2N - n). Also it decreases rapidly in the immediate neighborhood on both sides of these samples. Thus, with the peak of x,(n) occurring at n = np, not only the CBM 22:1/2-B
I. S. N. MURTHY and U. C. NIRANJAN
16
peak value in IX,(Q)] occurs at 8 = np but also the rectified time signal xi(n) looks like the one-sided DFTM in equation (9) under no leakage condition. This concept can be extended to any time signal x(n) with M arbitrary number of component waves. As such, now our claim is that, one can obtain filter characteristics namely the bandwidth, cut-off frequencies and transition bands to recover the component waves from the time signal x(n) itself rather than from IX(Q)]. 3.3. Filter design If a single sample at n = no in an arbitrary signal x(n) is to be recovered, then X(k), its DCT need be filtered with a notch filter of center frequency non/N rad and the inverse transform of the resulting filter output must be computed. Similarly, a set of consecutive samples corresponding to a component wave, x1(n), that exists between sample numbers no and n1 - 1 in x(n) can be recovered by filtering X(k) with a filter whose bandwidth is (n, -n,)rclN rad and cut-off frequencies at n&N and n&N rad. If x1(n) is at the beginning, at an intermediate position or at the end of x(n), then the corresponding filter will be of lowpass, bandpass or highpass type, respectively. Now, consider a signal x(n) of length N+ 1 samples, having M component waves. Let the beginning and end sample numbers of these component waves be no and n, - 1, n, , IZ~_~ and n,+,- 1, respectively, i.e. andn,-l,n,andn,-l,...
X(n)
= i:
Xi(n)
I=1
=x(&u(n-n;_,)-u(n-n,)];n=O,1,2,...,N
(10)
i=l
where no = 0 and II~ = N + 1. This is shown in Fig. lb for a normal ECG with P, QRS and T waves, i.e. for M = 3. The constituent samples in each of the M component waves are mapped as cosinusoids in the frequency bands -rig/N to (n, - 1)x/N, n&N to (n2 - l)n/N, . . . , (n,_,)n/N to (nM- 1)nlN rad respectively in the DCT domain. Therefore, to recover a component wavex;(n),i=1,2,. . . , M we convolve X(k) with a filter whose lower and upper cut-off frequencies are (n,_&N) and (n;- 1)nlN rad. The beginning and end sample numbers n,_, and n; - 1 of the ith component wave in the time signal can be determined using a simple rule as explained in Section 4.2, to fix the filter cut-off frequencies. These two samples need not be symmetrically placed with reference to the peak sample np of a component wave. 4. ANALYSIS 4.1. Data acquisition and evaluation ECG waveforms have been recorded using Franks corrected orthogonal lead system from different subjects and were digitized at 500 sps using a 12 bit ADC. The required computer programs have been developed in Matlab and all processing has been carried out on a PC-386 system with a co-processor. The quality of the recovered waveform is assessed by a cardiologist and also by computing the percentage rms difference (PRD), defined in [lo]. The PRD for the illustrative examples are shown in Figs 3 and 4. For a component wave the PRD is computed by considering the original and reconstructed signals over the points of interest which are identified as the beginning and end samples of the component wave. 4.2. Scheme of analysis The proposed scheme of analysis is shown in the block diagram of Fig. la. The signal x(n) in a given frame is rectified and smoothed by 5 point moving average technique. The
Decomposition
of ECG by linear filtering
17
DCT Beat
. . tm(n) component
m
(b)
l ../3 Time
sample
-’
No.
Fig. 1. (a) Scheme of analysis for component delineation-block signal with components marked.
diagram. (b) Digitized ECG
averaged signal is scanned through and the local peak values are determined. A search is then made in the neighborhood of each peak and the two samples with minimum amplitudes are determined as the beginning (xb) and end samples (xc) of the component wave. Then, it is checked if these minimum amplitude values are retained in at least the next 10 consecutive samples within a limit of + 10% of the peak amplitude. If this criterion is satisfied, then the cut-off frequencies of the filter are decided from sample numbers xb and x,. If the criterion is satisfied by more than 10 consecutive samples as in the case of normal ECG signal, then, it means that the two adjacent component waves under consideration are well separated. The rectified time signal is used only for the purpose of determining the beginning and end sample numbers of a component wave as required in the filter design. The actual filtering operation is carried out on X(k), the DCT of the given signal x(n). In the case of ECG signals with arrhythmias no clear separation may exist between certain groups of components. In those cases we filter out all such components together to get a compound wave. The filter output when inverse transformed gives the desired component wave. The algebraic sum of all the delineated component waves gives back the original signal. 4.3. Results Example 1. Consider three signals x,(n), x*(n), x3(n) of length 100 samples, sampled at 100 sps each with a single impulse at n= 1, 49, and 98 as shown in Fig. 2a, g, m.
I
50
(m)
loo
I
0
0.2 ,
s (n)
(h)
100
0.0
' o
a.0 ,
L
.25 (0)
-aJ
(i)
-5
-0.2-
0
!I0 (P!
100
-0.20
so
(9)
100
-cm0
.25
(4
.5
Fig. 2. (a) Signal x,(n) with a single impulse at n =
I.
1. (b) X,(k), the Da of the signal in (a). (c) ]X,(Q)l, the magnitude spectrum of the signal in (b). X-axis, normalized frequency number; Y-axis, arbitrary units. (d) a,(n), the recovered signal obtained as IDCT of X,(k), the filtered signal shown in (e). (e) X,(k), the signal obtained by filtering the signal in (b) with the filter of(f). (f) Log magnitude spectrum of the filter designed from the time signal x1(n) in (a). (g) Signal x2(n) with an impulse at n = 49. (h) X,(k), the DCT of the signal in (8). (i) IX,(Q)l, the DFTM of the signal in (h). X-axis, normalized frequency number; Y-axis, arbitrary units. (j) .f,(n), the recovered signal obtained as IDCT of X,(k), the filtered signal shown in (k). (k) X2(k), the signal obtained by filtering the signal in (h) with the filter of (1). (I) Log magnitude spectrum of the filter obtained from the time signal x2(n) in (8). (m) Signal x,(n) with a single impulse at n = 98. (n) XX/c), the DCT of the signal in (m). (0) IX,(Q)l, the magnitude spectrum of the signal in (n). X-axis, normalized frequency number; Y-axis, arbitrary units. (p) R&r), the recovered signal obtained as the IDCT of the filtered signal in (q). (q) X3(k), the signal obtained by filtering the signal in (n) with the filter of(r). (r) Log magnitude spectrum of the filter obtained from the time signal x3(n) in (m).
0
o.nJ
(9)
Decomposition
of ECG by linear filtering
19
Neglecting
the side lobes from their single sided DFTMs shown in Fig. 2c, i, o, we notice that they resemble the time signals. To recover these impulses one needs to use notch filters centered at frequencies of 0.5,24.5, and 49 Hz such as the ones shown in Fig. 2f, 1, r. The DCTs in Fig. 2b, h, n are convolved with the respective filters and the filter outputs are shown in Fig. 2e, k, q. The time signals recovered by computing the IDCT of the three filtered outputs are shown in Fig. 2d, j, p. A comparison of these recovered impulses with the original ones shows very good fit between the two sets but for a slight negative shift and a small ringing phenomena caused due to the filtering process.
Example 2. The signal in this case is a left bundle branch block (LBBB) type beat with a characteristic triphasic QRS complex and is of length 251 samples with three component waves. The signal x(n), its DCT X(k) and the 2N point DFIM, IX(Q) 1are shown in Fig. 3a, 3e and 3i respectively. The cut-off frequencies of the filters whose log magnitude spectra are shown in Fig. 3j, k and 1 are decided by the beginning and end sample numbers 0, 90; 91, 150; and 151, 250 of the three component waves. The filters’ cut-off frequencies are modified such that their transition bands are 85~r/250 to 95n/250, 145nl 250 to 155x/250 and 155~/200 to 165x/200 rad respectively. The output signals of the three filters are shown in Fig. 3f-h along with their IDCTs in Fig. 3b-d. The PRDs for the components turned out to be 5.3,4.1,14.9 and are marked on the respective figures. Example 3. As the next example, we consider the case of ventricular flutter or fibrillation. Though the signal in this case is characterized by closely packed repetitive peaks, it did not pose any special problems in the design of filters. The time signal x(n), its DCT X(k), is shown in Fig. 4a, f (for clarity only 71 DCT coefficients are shown). In this case the preprocessing algorithm identified four component waves. The log magnitude spectra of the filters designed are shown in Fig. 4g-j. The IDCT of the filter outputs are shown in Fig. 4b-e. In this example it is not possible to separate out the smaller components from the larger ones as no clear separation of minimum 10 to 15 samples, exists between them. 4.4. Application to real ECG
Continuous strips of ECG data containing about 300 beats recorded from different patients have been analyzed in close association with a cardiologist. The results are evaluated by the clinician and by computing the PRD. It was found that the original and delineated signals in most of the cases were visually indistinguishable. The computational complexity of the method depended upon how often the signal changed its nature in its time course. The method worked quite satisfactorily even on signals with changing heart rates in a frame of analysis and was independent of the shape, size and location of the components waves. 5. DISCUSSION
AND
SUMMARY
5.1. Discussion For a signal sampled at 500 sps, we have used filters with transition bands varying from 6.25 Hz to 20 Hz, with 300 coefficients. When the transition bands become very narrow, then a large number of filter coefficients are needed to keep the PRD low. This necessitates higher computational effort in the filtering process. The computational requirements of the method are computation of DCT, linear convolution, IDCT and FIR filter design. The increase in percentage error of the delineated components compared to that of the whole signal could be attributed to the fact that some finite amount of power (of about 10%) of the component waves leaks to the adjacent filter passbands. The reconstructed component waves have small d.c. shifts and the errors due to these shifts get compensated when all the component waves are added together. When a signal is contaminated with noise the method can be applied efficiently on signals obtained after noise removal.
20
NIRANJAN
I. S. N. MUXTHYand U. C.
8 d
:-:
~
0 ?
In
d
P 4i
?
a0 i
Decomposition LEGEND:
ABSCISSA:
-*w 0
lwd”
FIGS. (a)C FIG. (f) D? FIGS. (C,)-(J)
SAMPLE NO. COEFFICIENT NORMALIZED
of ECG by linear filtering NO. FREQUENCY
NO
21
ORDINATE: FIGS. (a)c) AMPLITUDE IN ARBITRARY unit FIG. (f) D & T COEFFICIENT AMPLITUDE FIGS. (Q&(j) LOG MAGNITUDE IN DECIBELS
m
whole beat PRDlB
PRO 0.8
Fig. 4. (a) Signal with onset of ventricular flutter or fibrillation with four components. (b) Z,(n), recovered first component wave. (c) _&(n), recovered second component wave. (d) &(n), recovered third component wave. (e) .&(n), recovered fourth component wave. (f) X(k), the DCT of the signal in (a), with first 71 coefficients. (g) Log magnitude spectrum of the filter obtained from the signal in (a) for the recovery of the first component wave. (h) Log magnitude spectrum of the filter obtained from the signal in (a) for the recovery of the second component wave. (i) Log magnitude spectrum of the filter obtained from the time signal in (a) for the recovery of the third component wave. (j) Log magnitude spectrum of the filter obtained from the signal in (a) for the recovery of the fourth component wave.
5.2.Summary We have shown that different component waves in a given ECG signal are mapped into distinct and non-overlapping frequency bands in the DCT domain. Since the rectified time signal resembled the DFIM of the transformed signal, appropriate filters are designed from the time signal itself. The transformed signal is then convolved with these filters and the component waves are obtained by computing the IDCT of the filtered signals. There was a good visual fit between the original and delineated components. The complete beat and large amplitude components have a PRD of 510%) whereas the small amplitude P wave has a PRD of 10-U%. The method is simple and is illustrated with suitable examples.
REFERENCES 1. J. Pan and W. J. Tompkins, A real-time QRS detection algorithm, IEEE Trans. BME 32,230-236 (1985). 2. A. Lightenberg and M. Kunt, A robust digital QRS detection algorithm for arrhythmia monitoring, Compur. Biomed. Res. 16, 273-286 (1983). 3. 0. Pahlm and L. Sornmo, Software QRS detection in ambulatory monitoring - A review, Med. Biol. Engng Comput. 22, 289-297 (1984). 4. B. C. Yu, et al. A nonlinear digital filter for cardiac QRS complex detection, J. clin. Engng 10, 193-201 (1985). 5. J. J. Hegenveld and J. H. Van Bemmel, Computer detection of P-waves, Comput. Biomed. Res. 9, 125 (1976). 6. R. A. Dufault, P-wave detection in the surface ECG via the LMS algorithm, Proc. 39th ACEMB Conf. 8 (1986). 7. L. R. Rabiner and B. Gold, Theory and Applications of Digital Signa/ Processing. Prentice-Hall, Englewood Cliffs, N.J., U.S.A. (1975). 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. BME 31, 702-706 (1984). 9. D. F. Elliott, Handbook of Digiral Signal Processing. Engineering Applications, Academic Press (1987). 10. J. P. Abenstein and W. J. Tompkins, A new data reduction algorithm for real-time ECG analysis, IEEE Trans. BME 29, 43-48 (1982).
About the Author-I. S. N. MURTHYreceived his Ph.D. degree in 1968 from the Indian Institute of Science, Bangalore and currently holds the position of a Professor in the Department of Electrical Engineering. He has several publications in control theory, pattern recognition and
22
I. S. N. MURTHY and U. C. NIRANJAN
ECG processing. His current research interests are in biomedical signal processing, biorobotics, and biocybernetics. He is the Editor of the Medical and Life Sciences Engineering Journal of the BME Society of India. About the Author-U. CHOLAYYA NIRANJANreceived his B.E. and M.Tech. degree in Electrical Engineering in 1985 and 1988 respectively. For his Master’s dissertation he worked on aircraft system identification at National Aeronautical Laboratory, Bangalore. After working as a lecturer in Manipal Institute of Technology, for two years, he is working towards the Ph.D. degree at the Department of Electrical Engineering, Indian Institute of Science, Bangalore. His research interests are adaptive signal processing, spectrum estimation, modeling and system identification. He is presently a UGC Research Fellow and is a student member of IEEE.
13-22, 1992 0
DECOMPOSITION
OF ECG BY LINEAR
I. S. N. Department
of Electrical
MURTHY*
Engineering,
Indian
and U. C. Institute
(Receiued 24 June 1991; in revised form 23 September November 1991)
0010.4825192 $55.00+ .oO IY92 Pergamon Press plc
FILTERING
NIRANJAN
of Science,
Bangalore
560 012, India
1991; received for publication 25
Abstract-A simple method is developed for the delineation of a given electrocardiogram (ECG) signal into its component waves. The properties of discrete cosine transform (DCT) are exploited for the purpose. The transformed signal is convolved with appropriate filters and the component waves are obtained by computing the inverse transform (IDCT) of the filtered signals. The filters are derived from the time signal itself. Analysis of continuous strips of ECG signals with various arrhythmias showed that the performance of the method is satisfactory both qualitatively and quantitatively. The small amplitude P wave usually had a high percentage rms difference (PRD) compared to the other large component waves. ECG analysis PRD DFTM
Component delineation FIR filter design
1.
Arrhythmia
Linear
filtering
DCT
INTRODUCTION
Delineation of the component waves in an ECG signal is vital for the detection of cardiac arrhythmias of critically ill patients in an intensive coronary care unit (ICCU). This topic has received considerable attention in the literature. Fairly accurate methods [l-4] are available for the detection of the QRS complex. However, the twin factors of small size and absence of spectral variability for temporally differing P waves have inhibited all the efforts for its detection [5,6]. Thus, one of the main obstacles for the development of a comprehensive arrhythmia analysis system has been the lack of simple and reliable methods for P wave detection. In the present paper a new approach based on a simple linear filtering process is proposed for the delineation of component waves in ECG irrespective of their shape, size and location. The linearity and frequency transformation properties of the discrete cosine transform (DCT) are exploited to convert each component wave into cosinusoids of a distinct frequency band in the transform domain. Furthermore the magnitude spectrum of the transformed signal and the rectified time signal look alike. As such the required filter characteristics such as lowpass (LP ), bandpass (BP) and highpass (HP) can be inferred from each one of them. Appropriate FIR filters [7] are then designed using the information available from the time signal itself. In the case of signals with well separated component waves, flexible filters can be obtained making the delineation relatively easy. For signals with arrhythmias where such clear separation between component waves does not exist, compound waves, containing two or more components are delineated. In contrast to this, the spectra of the component waves in the direct time signal overlap with each other and do not permit recovery of the components by such simple linear filtering [S]. Continuous strips of ECG signals sampled at 500 samples per second (sps) containing about 300 beats have been analyzed. The results obtained have been found to be satisfactory both qualitatively and quantitatively. In Section 2, the DCT and its linearity property are explained. The frequency transformation property of the DCT, derivation of the mathematical expressions for discrete Fourier transform magnitude (DFTM) and filter characteristics are considered in Section 3. In the next section, the scheme of analysis, examples and application of the *To
whom
correspondence
should
be addressed. 13
I. S. N. MLJRTHY and U. C. NIRANIAN
14
method to real ECG are dealt with. Finally, discussion and summary are included in Section 5. 2. TRANSFORM
AND
PROPERTIES
The DCT of a discrete time signal x(n) of length (N+ 1) samples is defined [9] as X(k) = A/@%$,
i c&n) n=O
cos (nknlN);
k= 0,1,2,
...,N
(1)
n = 0, 1,2, . . . , N.
(2)
wherecj=l/fifori=OorN;Ci=lfori#OorN. The inverse transform (IDCT) of X(k) is x(n) = t/c2,N>c,
2
ckX(k) cos (nkn/N);
k=O
Consider an ECG signal x(n), with P, QRS and T waves -x,(n), x,(n), x,(n) each of length (N + 1) samples. Each component wave exists over a small set of samples with the remaining ones as zero, i.e. x(n) = c x,(n); n =0, 1,2, . . . , N. “P.9.’
(3)
Now, if DCT 44
=
z
then
x(n) =
c
i=p4,
-c(n) f
5
2 X,(k) =X(k). i=p.
(4)
q, 1
Thus, by the linearity property, the DCT of a complete beat is expressed as the sum of the transforms of the individual components and vice versa. In other words, addition and deletion of component waves to and from the complete beat is possible. 3. THEORY 3.1. Frequency transformation We now show that the DCT converts each sample in a signal x(n) into a cosinusoid of distinct frequency and a group of consecutive samples into cosinusoids in a frequency band. This property enables component delineation by simple filtering of the transformed signal. Consider a signal x(n), of length (N + 1) samples. Let x(n) contain a single impulse at n = no with an amplitude a,,. Then we see from equation (1) that X(k) is a cosine wave with frequency (n,,zlN) and amplitude -a,,, i.e. X(k) = d@%$+,ckanO cos (n&z/N);
k = 0, 1,2, . . . , N.
(5)
When this impulse is shifted from n = no = 1 to n = no = N in steps of one sample at a time, the frequency of oscillation of X(k) changes from (n/N) to n rad in steps of (n/N) rad. Thus, from a given X(k), the location and amplitude of the impulse in x(n) can be determined using equation (2). We now extend this argument to a set of consecutive samples located between time samples n = n, and n = nl - 1. This group of samples may
Decomposition
of ECG by linear filtering
15
take any arbitrary shape of a component wave such as a P, QRS or T. Expressing this group of samples, x,(n), in terms of x(n), q-1
xl(n)=x(n)
x
6(n-m)
fornO~n~nl-l fornl
Substitution
Ocndn,-1.
of equation (6) into equation (1) leads to
X,(k) =
qzqc,[c~(no)
cos (n&lN)
+ c,+1x(no+ 1)
cos ((no + 1)knlN) + * . . + c,,_,x(q k=O,1,2,..
- 1) cos ((n1- 1) kJr/N)]
. ,N.
(7)
We note that X,(k), the DCT of xl(n) is made up of an equal number of (n,-no) cosinusoids each with a distinct frequency (nnlN) rad and amplitude ma,,; , N. It can be shown that the sum of these cosinusoids in X,(k) leads to a n=O,1,2,... damped cosine wave whose frequency of oscillation is determined by the sample number which has the largest amplitude in xl(n). We can extend this argument to an arbitrary waveform such as the ECG containing more than one component wave. From the linearity property in (4), it is not difficult to see that each component wave in x(n) gets transformed into cosinusoids in a distinct frequency band in X(k). The component waves can be separated out from each other by filtering X(k) in the DCT domain. On the other hand such a simple separation of the component waves is not possible from the direct time signal x(n). 3.2. Discrete Fourier trunsform magnitude (DFTM)
of X(k)
Consider X(n), the 2N point discrete Fourier transform sequence X(k) [9], ZN-
X(a) = x
(DFT) of the (N+ 1) point
1
X(k) e-iznkn’2N;S2 = 0, 1,2, . . . ,2N - 1
(8)
k=O
whereX(k)=OforkaN+l. Substituting X,(k) from equation (7) for X(k) in equation (8), (X,(S2)( the 2N point DFTM of X,(k) is obtained as ZN-
IX,(Q)1 =
I
C cfix,(n)(l/V@Vj) n=O
sin (n(N + 1) (n - Sr)/ZN) e+j”((“-R)‘2) sin (x(n - s2)/2N)
+ e_jn((n+Rj,2j sin (n(N+ I) (n + &)/2N) sin (n(n + 52)/2N) + (1- ti){V5
+ l/ti((
- l)‘“_*‘+ (- l)‘“+*‘)} II sZ=O,1,2 ,...,
2N-1.
(9)
Here, we note that the DFT coefficients X,(a), are the weighted sum of the constituent samples in xl(n). The weighting function inside the square brackets on the right-hand-side (r.h.s.) of equation (9) attains its maximum value in magnitude when SL= n and 52= (2N - n). Also it decreases rapidly in the immediate neighborhood on both sides of these samples. Thus, with the peak of x,(n) occurring at n = np, not only the CBM 22:1/2-B
I. S. N. MURTHY and U. C. NIRANJAN
16
peak value in IX,(Q)] occurs at 8 = np but also the rectified time signal xi(n) looks like the one-sided DFTM in equation (9) under no leakage condition. This concept can be extended to any time signal x(n) with M arbitrary number of component waves. As such, now our claim is that, one can obtain filter characteristics namely the bandwidth, cut-off frequencies and transition bands to recover the component waves from the time signal x(n) itself rather than from IX(Q)]. 3.3. Filter design If a single sample at n = no in an arbitrary signal x(n) is to be recovered, then X(k), its DCT need be filtered with a notch filter of center frequency non/N rad and the inverse transform of the resulting filter output must be computed. Similarly, a set of consecutive samples corresponding to a component wave, x1(n), that exists between sample numbers no and n1 - 1 in x(n) can be recovered by filtering X(k) with a filter whose bandwidth is (n, -n,)rclN rad and cut-off frequencies at n&N and n&N rad. If x1(n) is at the beginning, at an intermediate position or at the end of x(n), then the corresponding filter will be of lowpass, bandpass or highpass type, respectively. Now, consider a signal x(n) of length N+ 1 samples, having M component waves. Let the beginning and end sample numbers of these component waves be no and n, - 1, n, , IZ~_~ and n,+,- 1, respectively, i.e. andn,-l,n,andn,-l,...
X(n)
= i:
Xi(n)
I=1
=x(&u(n-n;_,)-u(n-n,)];n=O,1,2,...,N
(10)
i=l
where no = 0 and II~ = N + 1. This is shown in Fig. lb for a normal ECG with P, QRS and T waves, i.e. for M = 3. The constituent samples in each of the M component waves are mapped as cosinusoids in the frequency bands -rig/N to (n, - 1)x/N, n&N to (n2 - l)n/N, . . . , (n,_,)n/N to (nM- 1)nlN rad respectively in the DCT domain. Therefore, to recover a component wavex;(n),i=1,2,. . . , M we convolve X(k) with a filter whose lower and upper cut-off frequencies are (n,_&N) and (n;- 1)nlN rad. The beginning and end sample numbers n,_, and n; - 1 of the ith component wave in the time signal can be determined using a simple rule as explained in Section 4.2, to fix the filter cut-off frequencies. These two samples need not be symmetrically placed with reference to the peak sample np of a component wave. 4. ANALYSIS 4.1. Data acquisition and evaluation ECG waveforms have been recorded using Franks corrected orthogonal lead system from different subjects and were digitized at 500 sps using a 12 bit ADC. The required computer programs have been developed in Matlab and all processing has been carried out on a PC-386 system with a co-processor. The quality of the recovered waveform is assessed by a cardiologist and also by computing the percentage rms difference (PRD), defined in [lo]. The PRD for the illustrative examples are shown in Figs 3 and 4. For a component wave the PRD is computed by considering the original and reconstructed signals over the points of interest which are identified as the beginning and end samples of the component wave. 4.2. Scheme of analysis The proposed scheme of analysis is shown in the block diagram of Fig. la. The signal x(n) in a given frame is rectified and smoothed by 5 point moving average technique. The
Decomposition
of ECG by linear filtering
17
DCT Beat
. . tm(n) component
m
(b)
l ../3 Time
sample
-’
No.
Fig. 1. (a) Scheme of analysis for component delineation-block signal with components marked.
diagram. (b) Digitized ECG
averaged signal is scanned through and the local peak values are determined. A search is then made in the neighborhood of each peak and the two samples with minimum amplitudes are determined as the beginning (xb) and end samples (xc) of the component wave. Then, it is checked if these minimum amplitude values are retained in at least the next 10 consecutive samples within a limit of + 10% of the peak amplitude. If this criterion is satisfied, then the cut-off frequencies of the filter are decided from sample numbers xb and x,. If the criterion is satisfied by more than 10 consecutive samples as in the case of normal ECG signal, then, it means that the two adjacent component waves under consideration are well separated. The rectified time signal is used only for the purpose of determining the beginning and end sample numbers of a component wave as required in the filter design. The actual filtering operation is carried out on X(k), the DCT of the given signal x(n). In the case of ECG signals with arrhythmias no clear separation may exist between certain groups of components. In those cases we filter out all such components together to get a compound wave. The filter output when inverse transformed gives the desired component wave. The algebraic sum of all the delineated component waves gives back the original signal. 4.3. Results Example 1. Consider three signals x,(n), x*(n), x3(n) of length 100 samples, sampled at 100 sps each with a single impulse at n= 1, 49, and 98 as shown in Fig. 2a, g, m.
I
50
(m)
loo
I
0
0.2 ,
s (n)
(h)
100
0.0
' o
a.0 ,
L
.25 (0)
-aJ
(i)
-5
-0.2-
0
!I0 (P!
100
-0.20
so
(9)
100
-cm0
.25
(4
.5
Fig. 2. (a) Signal x,(n) with a single impulse at n =
I.
1. (b) X,(k), the Da of the signal in (a). (c) ]X,(Q)l, the magnitude spectrum of the signal in (b). X-axis, normalized frequency number; Y-axis, arbitrary units. (d) a,(n), the recovered signal obtained as IDCT of X,(k), the filtered signal shown in (e). (e) X,(k), the signal obtained by filtering the signal in (b) with the filter of(f). (f) Log magnitude spectrum of the filter designed from the time signal x1(n) in (a). (g) Signal x2(n) with an impulse at n = 49. (h) X,(k), the DCT of the signal in (8). (i) IX,(Q)l, the DFTM of the signal in (h). X-axis, normalized frequency number; Y-axis, arbitrary units. (j) .f,(n), the recovered signal obtained as IDCT of X,(k), the filtered signal shown in (k). (k) X2(k), the signal obtained by filtering the signal in (h) with the filter of (1). (I) Log magnitude spectrum of the filter obtained from the time signal x2(n) in (8). (m) Signal x,(n) with a single impulse at n = 98. (n) XX/c), the DCT of the signal in (m). (0) IX,(Q)l, the magnitude spectrum of the signal in (n). X-axis, normalized frequency number; Y-axis, arbitrary units. (p) R&r), the recovered signal obtained as the IDCT of the filtered signal in (q). (q) X3(k), the signal obtained by filtering the signal in (n) with the filter of(r). (r) Log magnitude spectrum of the filter obtained from the time signal x3(n) in (m).
0
o.nJ
(9)
Decomposition
of ECG by linear filtering
19
Neglecting
the side lobes from their single sided DFTMs shown in Fig. 2c, i, o, we notice that they resemble the time signals. To recover these impulses one needs to use notch filters centered at frequencies of 0.5,24.5, and 49 Hz such as the ones shown in Fig. 2f, 1, r. The DCTs in Fig. 2b, h, n are convolved with the respective filters and the filter outputs are shown in Fig. 2e, k, q. The time signals recovered by computing the IDCT of the three filtered outputs are shown in Fig. 2d, j, p. A comparison of these recovered impulses with the original ones shows very good fit between the two sets but for a slight negative shift and a small ringing phenomena caused due to the filtering process.
Example 2. The signal in this case is a left bundle branch block (LBBB) type beat with a characteristic triphasic QRS complex and is of length 251 samples with three component waves. The signal x(n), its DCT X(k) and the 2N point DFIM, IX(Q) 1are shown in Fig. 3a, 3e and 3i respectively. The cut-off frequencies of the filters whose log magnitude spectra are shown in Fig. 3j, k and 1 are decided by the beginning and end sample numbers 0, 90; 91, 150; and 151, 250 of the three component waves. The filters’ cut-off frequencies are modified such that their transition bands are 85~r/250 to 95n/250, 145nl 250 to 155x/250 and 155~/200 to 165x/200 rad respectively. The output signals of the three filters are shown in Fig. 3f-h along with their IDCTs in Fig. 3b-d. The PRDs for the components turned out to be 5.3,4.1,14.9 and are marked on the respective figures. Example 3. As the next example, we consider the case of ventricular flutter or fibrillation. Though the signal in this case is characterized by closely packed repetitive peaks, it did not pose any special problems in the design of filters. The time signal x(n), its DCT X(k), is shown in Fig. 4a, f (for clarity only 71 DCT coefficients are shown). In this case the preprocessing algorithm identified four component waves. The log magnitude spectra of the filters designed are shown in Fig. 4g-j. The IDCT of the filter outputs are shown in Fig. 4b-e. In this example it is not possible to separate out the smaller components from the larger ones as no clear separation of minimum 10 to 15 samples, exists between them. 4.4. Application to real ECG
Continuous strips of ECG data containing about 300 beats recorded from different patients have been analyzed in close association with a cardiologist. The results are evaluated by the clinician and by computing the PRD. It was found that the original and delineated signals in most of the cases were visually indistinguishable. The computational complexity of the method depended upon how often the signal changed its nature in its time course. The method worked quite satisfactorily even on signals with changing heart rates in a frame of analysis and was independent of the shape, size and location of the components waves. 5. DISCUSSION
AND
SUMMARY
5.1. Discussion For a signal sampled at 500 sps, we have used filters with transition bands varying from 6.25 Hz to 20 Hz, with 300 coefficients. When the transition bands become very narrow, then a large number of filter coefficients are needed to keep the PRD low. This necessitates higher computational effort in the filtering process. The computational requirements of the method are computation of DCT, linear convolution, IDCT and FIR filter design. The increase in percentage error of the delineated components compared to that of the whole signal could be attributed to the fact that some finite amount of power (of about 10%) of the component waves leaks to the adjacent filter passbands. The reconstructed component waves have small d.c. shifts and the errors due to these shifts get compensated when all the component waves are added together. When a signal is contaminated with noise the method can be applied efficiently on signals obtained after noise removal.
20
NIRANJAN
I. S. N. MUXTHYand U. C.
8 d
:-:
~
0 ?
In
d
P 4i
?
a0 i
Decomposition LEGEND:
ABSCISSA:
-*w 0
lwd”
FIGS. (a)C FIG. (f) D? FIGS. (C,)-(J)
SAMPLE NO. COEFFICIENT NORMALIZED
of ECG by linear filtering NO. FREQUENCY
NO
21
ORDINATE: FIGS. (a)c) AMPLITUDE IN ARBITRARY unit FIG. (f) D & T COEFFICIENT AMPLITUDE FIGS. (Q&(j) LOG MAGNITUDE IN DECIBELS
m
whole beat PRDlB
PRO 0.8
Fig. 4. (a) Signal with onset of ventricular flutter or fibrillation with four components. (b) Z,(n), recovered first component wave. (c) _&(n), recovered second component wave. (d) &(n), recovered third component wave. (e) .&(n), recovered fourth component wave. (f) X(k), the DCT of the signal in (a), with first 71 coefficients. (g) Log magnitude spectrum of the filter obtained from the signal in (a) for the recovery of the first component wave. (h) Log magnitude spectrum of the filter obtained from the signal in (a) for the recovery of the second component wave. (i) Log magnitude spectrum of the filter obtained from the time signal in (a) for the recovery of the third component wave. (j) Log magnitude spectrum of the filter obtained from the signal in (a) for the recovery of the fourth component wave.
5.2.Summary We have shown that different component waves in a given ECG signal are mapped into distinct and non-overlapping frequency bands in the DCT domain. Since the rectified time signal resembled the DFIM of the transformed signal, appropriate filters are designed from the time signal itself. The transformed signal is then convolved with these filters and the component waves are obtained by computing the IDCT of the filtered signals. There was a good visual fit between the original and delineated components. The complete beat and large amplitude components have a PRD of 510%) whereas the small amplitude P wave has a PRD of 10-U%. The method is simple and is illustrated with suitable examples.
REFERENCES 1. J. Pan and W. J. Tompkins, A real-time QRS detection algorithm, IEEE Trans. BME 32,230-236 (1985). 2. A. Lightenberg and M. Kunt, A robust digital QRS detection algorithm for arrhythmia monitoring, Compur. Biomed. Res. 16, 273-286 (1983). 3. 0. Pahlm and L. Sornmo, Software QRS detection in ambulatory monitoring - A review, Med. Biol. Engng Comput. 22, 289-297 (1984). 4. B. C. Yu, et al. A nonlinear digital filter for cardiac QRS complex detection, J. clin. Engng 10, 193-201 (1985). 5. J. J. Hegenveld and J. H. Van Bemmel, Computer detection of P-waves, Comput. Biomed. Res. 9, 125 (1976). 6. R. A. Dufault, P-wave detection in the surface ECG via the LMS algorithm, Proc. 39th ACEMB Conf. 8 (1986). 7. L. R. Rabiner and B. Gold, Theory and Applications of Digital Signa/ Processing. Prentice-Hall, Englewood Cliffs, N.J., U.S.A. (1975). 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. BME 31, 702-706 (1984). 9. D. F. Elliott, Handbook of Digiral Signal Processing. Engineering Applications, Academic Press (1987). 10. J. P. Abenstein and W. J. Tompkins, A new data reduction algorithm for real-time ECG analysis, IEEE Trans. BME 29, 43-48 (1982).
About the Author-I. S. N. MURTHYreceived his Ph.D. degree in 1968 from the Indian Institute of Science, Bangalore and currently holds the position of a Professor in the Department of Electrical Engineering. He has several publications in control theory, pattern recognition and
22
I. S. N. MURTHY and U. C. NIRANJAN
ECG processing. His current research interests are in biomedical signal processing, biorobotics, and biocybernetics. He is the Editor of the Medical and Life Sciences Engineering Journal of the BME Society of India. About the Author-U. CHOLAYYA NIRANJANreceived his B.E. and M.Tech. degree in Electrical Engineering in 1985 and 1988 respectively. For his Master’s dissertation he worked on aircraft system identification at National Aeronautical Laboratory, Bangalore. After working as a lecturer in Manipal Institute of Technology, for two years, he is working towards the Ph.D. degree at the Department of Electrical Engineering, Indian Institute of Science, Bangalore. His research interests are adaptive signal processing, spectrum estimation, modeling and system identification. He is presently a UGC Research Fellow and is a student member of IEEE.