- Email: [email protected]
System identification for the ECG using CZT
COMPUTERS
AND
BIOMEDICAL
System
25, 407-416 (1992)
RESEARCH
Identification
for the ECG Using CZT
U. C. NIRANJAN AND I. S. N. MURTHY Department
of Electrical
Engineering,
Indian
Institute
of Science,
Bangalore
560 012, India
Received April 11. 1991
A new approach for extraction of clinically useful parameters from the ECG signal is presented using the system identification technique of CZT on the DCT-transformed signal. A one to one relationship between the model singularities and the significant points in the time signal is arrived at. The method allows the determination of the R-R interval needed in rhythm analysis. The complex cepstrum is used for identifying and removing the effect of zeros outside the unit circle. A significant data compression of 1 in 10 is achieved. A large number of continuous strips of ECG data are analyzed and the results are presented. 0 1992 Academx
Press. Inc.
1. INTRODUCTION The electrocardiogram (ECG) signal serves as a reliable tool for both morphological and rhythm analysis. The signal is recorded noninvasively in a simple way and is used extensively. Although it is well known that both parametric and nonparametric mathematical models have played a prominent role in the analysis and synthesis of many physiological systems, very few attempts have been made to study this clinically important signal from a system-theoretic point of view, in order to derive the much needed clinical parameters for diagnosis. Approximation of ECG as the impulse response of a system required models of very high order (I), which hindered the attempts to relate the model pole-zero combinations to the component waves in the signal and extraction of clinical information therefrom. In the present paper, we present a new approach for the extraction of clinically useful parameters from an ECG signal. The DCT of a given ECG signal is modeled as the impulse response of a pole-zero system using the chirp 2 transform (CZT) technique. Exploiting the frequency transformation property of the discrete cosine transform (DCT), a one to one relationship, between the model singularities (viz. the poles and zeros) and the significant points in the time signal, is arrived at. It is shown that this method of analysis accomplishes computation of R-R interval and delineation of component waves in a given signal. These claims are substantiated by the analysis of continuous strips of large amounts of ECG data from MIT and our own databases. The method 407 OOlO-4809/92 $5.00 Copyright 0 1992 by Academx Press. Inc. All rights of reproduction in any form reserved.
provides a data compression of 1 in 10. allows the e\;timation 01 model pole\ and zeros one by one in a natural way, eliminates the neceshit! of :t priori fixation of the unknown
order autoregressive methods
model order, and computation
and moving
average polynomials
of the rooks of the higher as required by other
(2. 3). 2. THEORY
2.1, The Chirp Z Transfbrm Apart from least squares the application (DFT) and CZT. The CZT as (4) I/- I S(z) = C s(n) Zi-‘I. ,i=n is through
Algorithm methods, another approach to system identification of transforms such as discrete Fourier transform of a signal s(n). of length ,W samples. is defined
whereZ,
= AW.-“./\
= O,i,3..
..M
- 1.
and where W, and A, are positive real numbers. The points Z,, at which the CZT is computed. lie along a spiral in the Z plane, where W, controls the rate at which the contour spirals and QO is the angular separation between two successive points. The parameters A, and 8, determine the radius and angle of the first point on the spiral. Rabiner et al. (5) have shown that the evaluation of Z transform along a spiral contour can be expressed as the convolution of the weighted signal with another weighting function. Corinthios (6) has shown that considerable reduction in evaluation of CZT in Eq. [l] can be achieved by making W, and e-j’, equal to unity. This in turn amounts to evaluation of the Z transform along constant damping and constant frequency contours, respectively. 2.2. Determir~ation
of System Function
Singularities
Consider a system which can be described by a pole-zero
model of the type
h=l
where factors inside outside The
Ickl and lekl are less than unity, and ldLl is greater than unity, so that (1 - C~Z-‘) and (1 - e,z- ‘) correspond to the model zeros and poles the unit circle and factors (1 - d,z- ‘1 correspond to the model zeros the unit circle. model singularities cL, dk, and e, in Eq. [2] are identified from S(z), the
ECG
SYSTEM
IDENTIFICATION
VIA
409
CZT
2 transform of the signal s(n), such that h(n), the model impulse response, closely approximates s(n). The magnitude spectrum computed along the semicircle (constant damping contour) shows a peak or valley when the contour passes closer to a pole or zero, while the phase spectrum computed along the radial line (constant frequency contour) closer to this peak or valley shows a phase change close to 27r rad. Thus the circular coordinates of the point at which these abrupt changes in magnitude and phase spectrum occur are taken as the magnitude and phase of the pole or zero. The response due to the identified zero or pole is removed from the signal by inverse filtering. The system identification procedure is then repeated on the reduced order signal until all the poles and zeros of the system are identified. 2.3. Removal
of Zeros Outside the Unit Circle
In our study, we noticed that the DCT of many abnormal ECG signals have zeros outside the unit circle and therefore cannot be canceled by inverse filtering. For this purpose the complex cepstrum is used, which converts the convolution operation into simple addition. Alternately, the outside zeros can also be estimated from the complex cepstrum. It is known that the complex cepstrum of a mixed phase signal will have both causal and anticausal parts whose inverse cepstrum corresponds to the minimum and maximum phase parts of the signal, respectively. In case of DCT of the ECG signal, the maximum phase part corresponds to the impulse response of the outside zeros. Also, the linear phase of the signal DCT (obtained from the unwrapped phase) is proportional to the number of outside zeros. The outside zeros can therefore easily be estimated by solving for the roots of the maximum phase signal. 2.4. Application
to Electrocardiogram
(ECG)
The DCT of a signal s(n) is defined as (7) S(k) = *IN)
Cr :r,: 2 s(n) cos [(y--)q
k=0,1,2,.
. .,N-
1, [31
where CL =
N-2,
k = 0,
1,
k# 0.
The inverse discrete cosine transform
(IDCT)
is defined as
N-l
s(n) = V@IN)
2 Ch S(k) h=O
cos
[(2n.&,l’kT],
n=0,1,2,.
. .,N-
1 [41
The DCT-transformed
ECG signal is amenable
for approximation
as the
410
NIKANJAN
AND
Ml.‘KIfil
impulse response of a pole-zero model, while the time domain a~gnal IS not (8). When the signal s(n) in Eq. 133 is expressed as the sum of N delayed impul$cs, the DCT is the sum of N cosinusoids ofamplitudes proportional to the magnitude of the signal samples, i.e.. V?%N)s(n) and frequencie\ proportional to the sample number, i.e.. (27 + lbri2N rad. Due to this frequency transformation property of the DCT, the magnitude spectrum of the transformed signal and the rectified signal look alike. Thus, in the present study, since the poles and zeros in the model of a given signal DCT represent the peaks and valleys in the DCT spectrum, there exists a direct one to one relationship between model singularities and extremum points in the time signal. Also a pole represents a prominent cosinusoid in the transform of a component wave and therefore corresponds to peak sample number in the time signal. Thus, if a peak occurs at sample number /i in a signal of length N samples. then the approximate position of the pole in Z plane would be ri-rnlil(N ~ I) rad. This property of the DCT enables one to choose appropriate initial contours to be scanned along in the Z plane, as the locations of the signal peaks are approximately known. In case of real ECG. the magnitude spectrum of the transformed signal has sometimes been found to be very noisy when evaluated deep inside the unit circle. This in turn resulted in incorrect estimates of magnitude Y and phase 0 of a pole or zero. The magnitude of the pole corresponding to the broad T wave was actually in the range 0.65 to 0.75. However. when the spectrum is noisy. the CZT identified the magnitude of the T-wave pole in the range 0.9 to 0.98. To overcome this problem, a parabolic fit was made to the T-wave part of the smoothed magnitude spectrum by the least squares technique. The parabola has the following form (9), near the peak sample 01of the smoothed spectrum. y(cu) = m2 + ha + (’
I.Fl
Having determined the parameters N, h. and c. the spectral peak location ctP, center frequency P and bandwidth d of the pole in the S-domain are computed as (9) ap = - bJ2u i- =
(17,
B = -{b?
+ a,)f,lN
[61
- 4a[c - osy(a,)]}“y; aN
where n,, = location of discrete spectral peak f, = sampling frequency in s/set. N = number of discrete samples in the magnitude Using the S-plane to Z-plane transformation phase of the T-wave pole are computed as
spectrum.
z = exp(s/j,), the magnitude
and
ECG ,z
SYSTEM
IDENTIFICATION
VIA
411
CZT
0.40
5
: 015 n; -0.10 .E $ -0.35 3 + g -0.60 a
I
0
100
200
0
I
100
Sample
No.
DCT Coefficient
Sample
No.
DCT Coefficient
I
L
200
No.
No.
Normalized
Frequency
No.
FIG. 1. (a) Premature Ventricular Contraction (PVC) signal. (b) DCT of signal in (a). (c) Magnitude spectrum of DCT in (b). (d) Reconstructed ECG signal from the model. (e) Reconstructed DCT signal from the model. (f) Magnitude spectrum of reconstructed DCT in (e).
The above method of estimating when the spectrum is noisy.
Y and 8 is used only for the T-wave pole,
3. EXAMPLES 3.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 coprocessor. The quality of the reconstructed signal is assessed by the measure normalized mean square error (NMSE), defined in (10). The clinical utility of the method is illustrated by the examples given below followed by the summary on analysis of a large number of ECG data. 3.2. Example
1
A Premature Ventricular Contraction (PVC) beat of length 256 samples (N) and its DCT are shown in Figs. la and lb. The complex cepstrum and 512-point CZT of the DCT indicated a complex zero outside the unit circle at Y = 1.4 and
412
NIKANJAN
AND MLIKTH\I
TAHLt POLES.
ZEROS,
4ND
i l%Ah
F\
Pole<
Zeros Magnitude I
Phase (4 (in degrees)
PVC in Fig. la
0.1 0.91 I .4 0.87 0.94
0 k 20 + 36 k 64.93 I80
0.90 0.92 0.89 0.74
LBBB in Fig. 2b
1.2 1.3 0.78 I .53 1.3
0 k 41.83 f 60.78 t 69.2 180
0.87 0.94 0.87 0.67
Signal
SAMPI
Magnitude I
Phax (4 (in degrees)
Acrual peah sample
t’\tlmattxi peak sample 4x 60 96 IX0
+ ~
31,
-.-
t 64.98 2 78.23 f 122.8
il.7 92 II? 180
36 92 III I74
t3 = ?36”. Using the complex cepstrum technique the response of the outside zero was deconvolved from the signal DCT to obtain the minimum phase component (MPC). The magnitude spectrum of this MPC computed along a semicircle of radius 0.92 showed sharp peaks at angles 1~42”. The phase spectrum evaluated along radial lines between 38” and 46” also showed a sharp change of about 25~rad between dampings 0.92 and 0.93, confirming the presence of a pole. The response due to this identified pole was deconvolved from the MPC by inverse filtering. This procedure was repeated until all the poles and zeros were identified. The reconstructed DCT and signal with an NMSE of 10.2% are shown in Figs. le and Id. The identified poles and zeros are given in Table 1. The pole-zero plot is shown in Fig. 3a. Now converting the pole angles (0) 33.75”, 42”, 67.5”, and 127” to sample numbers by the inverse relation (N - 1)8/180, we obtain 48,60,96, and 180. Correlating these sample numbers with Fig. la, we see that they correspond to the R wave, the S wave, the kink in the QRS after the S wave, and the peak in the T wave, respectively. By a similar analysis we see that the zeros at 0”. 20”, and 180” correspond to the isoelectric regions preceding the QRS and after the T wave, while those at 36 and 64.93” correspond to the points in the QRS at which the signal crosses the baseline. The actual peak samples and those estimated from pole locations are given in Table 1. The maximum deviation between the two in most cases is around ?2 samples. Also the magnitude spectra of the DCT of the original and reconstructed signal, shown in Figs. lc and If, respectively, resemble the rectified and smoothed ECG signal.
ECG SYSTEM 320
IDENTIFICATION
VIA CZT
413
a
160 270
320
d
160
320 e 160 0 -160 -320
k 0
100
200
FIG. 2. Sample No. vs Amplitude. (a) Simulated ECG data obtained by Steiglitz-McBride model for the signal in (b). (b) Left Bundle Branch Block (LBBB) beat. (c) Normal ECG signal. (d) Reconstructed signal by the CZT model for the simulated signal in (a). (e) Reconstructed signal by the model for the actual signal in (b). (f) Reconstructed signal by the model for the normal ECG in (c).
3.3. Example 2
The DCT of a Left Bundle Branch Block (LBBB) beat, shown in Fig. 2b is modeled by the Steiglitz-McBride (SM) method (2). The IDCT of the identified model impulse response shown in Fig. 2a is used as the simulated ECG data for analysis. For this case as well as for other simulated data, the identification of model singularities by the CZT technique was found to be relatively easy. The reconstructed signal with an NMSE of 3.8% is shown in Fig. 2d. The model was then estimated using the DCT of the actual signal shown in Fig. 2b. However, in this case the magnitude spectrum of the CZT was noisy and did not show very narrow peaks or valleys. After the larger magnitude poles and zeros were identified and deconvolved from the signal, the T-wave pole was identified from the parameters of the parabolic fit made to the dominant part of the magnitude spectrum. The reconstructed signal is shown in Fig. 2e. The identified model singularities and the peak sample numbers are given in Table 1, while the pole-zero plot is shown in Fig. 3b. Due to the noisy spectrum, the identified parameters were a little inaccurate, which resulted in a higher NMSE of 14.8%. Next consider a normal ECG signal shown in Fig. 2c. All the model zeros for this signal were lying inside the unit circle and thus the model was identified with relative ease. The reconstructed signal is shown in Fig. 2f.
414
-15
NIRANJAN
-1 0
-0 5
FIG. 3. CZT Fig. 2(b).
0.0
pole-zero
0.5
10
plot [CX) pole.
AND
1.5
:MI.Kl
-.
(0) zero]
L. 3
lik
*
i -(
0
--AL,i---05
for (a) PVC signal
..L. OS
00
of Fig.
I(a).
~... ..-i.-+.-.. ‘0
rhr LBBB
‘5
heat of
3.4. Example 3 In this example we have considered a continuous strip of ECG of length 500 samples containing two beats, as shown in Fig. 4a. The poles and zeros identified from the signal CZT, shown in Fig. 4b, are distributed in the Z plane according to the locations of the component waves in the time signal. The poles with largest magnitude in the clusters corresponding to two QRS complexes are separated by an angle of 85”. This corresponds to a separation of 236 samples between the respective R-wave peaks in the time signal, which is the actual z c3 x 0 2 b 3c a -E g
04
b
a
r-pole.
0.2
o-zero
1 t
0 -0.2
0
500 samole
-0.1
-0.05 difference
no
0
0.05
in amplitudes
0.1
0 -10
-5 difference
0
5
?O
in intervals
FIG. 4. Continuous strip of ECG signal containing two beats. I-) Original signal. t”.) Reconstructed signal. tb) Pole-zero plot [(x) pole, (0) zero] for the DCT ofthe signal in (a). (c) A histogram representing the differences in amplitudes of samples in the 200 actual and reconstructed ECG beats, analyzed. td) A histogram representing the differences in actual and measured R-R interval. for 200 ECG beats.
ECG SYSTEM IDENTIFICATION
VIACZT
415
R-R interval for the signal shown. When signal frames of smaller durations with single R waves are considered, the angular separation between predominant pole clusters can be carried forward to the next frame, to determine the R-R interval. Continuous strips of ECG data from our own and MIT databases containing 200 beats have been analyzed. A histogram representing the difference in amplitudes between the actual and the reconstructed signals is shown in Fig. 4c. In order to make a uniform comparison, the original signal was normalized to have amplitudes between 0 and I, and the reconstructed signal was normalized by the same scale factor. As can be seen from the histogram, the maximum difference in amplitudes is 50.075, while most of the differences are less than kO.016. Another histogram representing the discrepancy in the actual and measured R-R interval for 200 ECG beats is shown in Fig. 4d. The maximum difference in R-R interval is r6 samples; however, we do not see any specific pattern in the distribution of the differences from Fig. 4d. 4. DISCUSSIONS The frequency transformation property of DCT allowed demarcation of the isoelectric region, baseline crossing points and peaks in the component waves from the angles subtended by the zeros and poles in the Z plane. Thus the polar plot can be visualized as a representation of significant points in the signal and thus can be made use of for the delineation of component waves. As the QRS usually requires two or more pairs of model poles for adequate representation, it is identified in a given frame of signal by searching for pole clusters. The angular separation between the QRS pole clusters is converted to time information and thereby heart rate. Identification of T-wave pole, which is usually located deep inside the unit circle, posed some problems when the spectrum was noisy. In such cases the pole was estimated from the parameters of a leastsquares parabolic fit made to the dominant part of the spectrum. On an average a signal of length 250 samples can be represented by about six complex conjugate or its equivalent real pole-zero pairs. This results in a data compression of a little more than 1 in 10. The reconstructed signals from the model parameters are of good quality with NMSEs less than 15%. On average the CZT was computed about 12 to 15 times to identify all the model singularities, depending upon the complexity of the ECG signal. However, it must be noted that an erroneously identified pole or zero will introduce its counterpart in the signal while inverse filtering. 5. CONCLUSIONS
A simple method based on repetitive scanning of the Z plane of the discrete cosine transformed ECG signal by the chirp Z transform technique is proposed for the identification of poles and zeros. We demonstrated the utility of the method by relating the model parameters to the features in the waveform.
3Ih
NIRANJAN
AND
MUKI
tti
Specifically, we obtained the R-R interval and locations ot component and isoelectric regions, in addition to a data compression of 1 in IO.
/. 2. 3. -1. 5. 6. 7. 8. Y. IO.
tiave\
MURTHY. I. S. N., RANG4RU. M. R., UDUPA. K. J.. ANDCJOYAI ,\. K. Homom~)rphlc andlyb~\ and modelling of ECG hignats. lEE.E Trur1.r. Biotnr
AND
BIOMEDICAL
System
25, 407-416 (1992)
RESEARCH
Identification
for the ECG Using CZT
U. C. NIRANJAN AND I. S. N. MURTHY Department
of Electrical
Engineering,
Indian
Institute
of Science,
Bangalore
560 012, India
Received April 11. 1991
A new approach for extraction of clinically useful parameters from the ECG signal is presented using the system identification technique of CZT on the DCT-transformed signal. A one to one relationship between the model singularities and the significant points in the time signal is arrived at. The method allows the determination of the R-R interval needed in rhythm analysis. The complex cepstrum is used for identifying and removing the effect of zeros outside the unit circle. A significant data compression of 1 in 10 is achieved. A large number of continuous strips of ECG data are analyzed and the results are presented. 0 1992 Academx
Press. Inc.
1. INTRODUCTION The electrocardiogram (ECG) signal serves as a reliable tool for both morphological and rhythm analysis. The signal is recorded noninvasively in a simple way and is used extensively. Although it is well known that both parametric and nonparametric mathematical models have played a prominent role in the analysis and synthesis of many physiological systems, very few attempts have been made to study this clinically important signal from a system-theoretic point of view, in order to derive the much needed clinical parameters for diagnosis. Approximation of ECG as the impulse response of a system required models of very high order (I), which hindered the attempts to relate the model pole-zero combinations to the component waves in the signal and extraction of clinical information therefrom. In the present paper, we present a new approach for the extraction of clinically useful parameters from an ECG signal. The DCT of a given ECG signal is modeled as the impulse response of a pole-zero system using the chirp 2 transform (CZT) technique. Exploiting the frequency transformation property of the discrete cosine transform (DCT), a one to one relationship, between the model singularities (viz. the poles and zeros) and the significant points in the time signal, is arrived at. It is shown that this method of analysis accomplishes computation of R-R interval and delineation of component waves in a given signal. These claims are substantiated by the analysis of continuous strips of large amounts of ECG data from MIT and our own databases. The method 407 OOlO-4809/92 $5.00 Copyright 0 1992 by Academx Press. Inc. All rights of reproduction in any form reserved.
provides a data compression of 1 in 10. allows the e\;timation 01 model pole\ and zeros one by one in a natural way, eliminates the neceshit! of :t priori fixation of the unknown
order autoregressive methods
model order, and computation
and moving
average polynomials
of the rooks of the higher as required by other
(2. 3). 2. THEORY
2.1, The Chirp Z Transfbrm Apart from least squares the application (DFT) and CZT. The CZT as (4) I/- I S(z) = C s(n) Zi-‘I. ,i=n is through
Algorithm methods, another approach to system identification of transforms such as discrete Fourier transform of a signal s(n). of length ,W samples. is defined
whereZ,
= AW.-“./\
= O,i,3..
..M
- 1.
and where W, and A, are positive real numbers. The points Z,, at which the CZT is computed. lie along a spiral in the Z plane, where W, controls the rate at which the contour spirals and QO is the angular separation between two successive points. The parameters A, and 8, determine the radius and angle of the first point on the spiral. Rabiner et al. (5) have shown that the evaluation of Z transform along a spiral contour can be expressed as the convolution of the weighted signal with another weighting function. Corinthios (6) has shown that considerable reduction in evaluation of CZT in Eq. [l] can be achieved by making W, and e-j’, equal to unity. This in turn amounts to evaluation of the Z transform along constant damping and constant frequency contours, respectively. 2.2. Determir~ation
of System Function
Singularities
Consider a system which can be described by a pole-zero
model of the type
h=l
where factors inside outside The
Ickl and lekl are less than unity, and ldLl is greater than unity, so that (1 - C~Z-‘) and (1 - e,z- ‘) correspond to the model zeros and poles the unit circle and factors (1 - d,z- ‘1 correspond to the model zeros the unit circle. model singularities cL, dk, and e, in Eq. [2] are identified from S(z), the
ECG
SYSTEM
IDENTIFICATION
VIA
409
CZT
2 transform of the signal s(n), such that h(n), the model impulse response, closely approximates s(n). The magnitude spectrum computed along the semicircle (constant damping contour) shows a peak or valley when the contour passes closer to a pole or zero, while the phase spectrum computed along the radial line (constant frequency contour) closer to this peak or valley shows a phase change close to 27r rad. Thus the circular coordinates of the point at which these abrupt changes in magnitude and phase spectrum occur are taken as the magnitude and phase of the pole or zero. The response due to the identified zero or pole is removed from the signal by inverse filtering. The system identification procedure is then repeated on the reduced order signal until all the poles and zeros of the system are identified. 2.3. Removal
of Zeros Outside the Unit Circle
In our study, we noticed that the DCT of many abnormal ECG signals have zeros outside the unit circle and therefore cannot be canceled by inverse filtering. For this purpose the complex cepstrum is used, which converts the convolution operation into simple addition. Alternately, the outside zeros can also be estimated from the complex cepstrum. It is known that the complex cepstrum of a mixed phase signal will have both causal and anticausal parts whose inverse cepstrum corresponds to the minimum and maximum phase parts of the signal, respectively. In case of DCT of the ECG signal, the maximum phase part corresponds to the impulse response of the outside zeros. Also, the linear phase of the signal DCT (obtained from the unwrapped phase) is proportional to the number of outside zeros. The outside zeros can therefore easily be estimated by solving for the roots of the maximum phase signal. 2.4. Application
to Electrocardiogram
(ECG)
The DCT of a signal s(n) is defined as (7) S(k) = *IN)
Cr :r,: 2 s(n) cos [(y--)q
k=0,1,2,.
. .,N-
1, [31
where CL =
N-2,
k = 0,
1,
k# 0.
The inverse discrete cosine transform
(IDCT)
is defined as
N-l
s(n) = V@IN)
2 Ch S(k) h=O
cos
[(2n.&,l’kT],
n=0,1,2,.
. .,N-
1 [41
The DCT-transformed
ECG signal is amenable
for approximation
as the
410
NIKANJAN
AND
Ml.‘KIfil
impulse response of a pole-zero model, while the time domain a~gnal IS not (8). When the signal s(n) in Eq. 133 is expressed as the sum of N delayed impul$cs, the DCT is the sum of N cosinusoids ofamplitudes proportional to the magnitude of the signal samples, i.e.. V?%N)s(n) and frequencie\ proportional to the sample number, i.e.. (27 + lbri2N rad. Due to this frequency transformation property of the DCT, the magnitude spectrum of the transformed signal and the rectified signal look alike. Thus, in the present study, since the poles and zeros in the model of a given signal DCT represent the peaks and valleys in the DCT spectrum, there exists a direct one to one relationship between model singularities and extremum points in the time signal. Also a pole represents a prominent cosinusoid in the transform of a component wave and therefore corresponds to peak sample number in the time signal. Thus, if a peak occurs at sample number /i in a signal of length N samples. then the approximate position of the pole in Z plane would be ri-rnlil(N ~ I) rad. This property of the DCT enables one to choose appropriate initial contours to be scanned along in the Z plane, as the locations of the signal peaks are approximately known. In case of real ECG. the magnitude spectrum of the transformed signal has sometimes been found to be very noisy when evaluated deep inside the unit circle. This in turn resulted in incorrect estimates of magnitude Y and phase 0 of a pole or zero. The magnitude of the pole corresponding to the broad T wave was actually in the range 0.65 to 0.75. However. when the spectrum is noisy. the CZT identified the magnitude of the T-wave pole in the range 0.9 to 0.98. To overcome this problem, a parabolic fit was made to the T-wave part of the smoothed magnitude spectrum by the least squares technique. The parabola has the following form (9), near the peak sample 01of the smoothed spectrum. y(cu) = m2 + ha + (’
I.Fl
Having determined the parameters N, h. and c. the spectral peak location ctP, center frequency P and bandwidth d of the pole in the S-domain are computed as (9) ap = - bJ2u i- =
(17,
B = -{b?
+ a,)f,lN
[61
- 4a[c - osy(a,)]}“y; aN
where n,, = location of discrete spectral peak f, = sampling frequency in s/set. N = number of discrete samples in the magnitude Using the S-plane to Z-plane transformation phase of the T-wave pole are computed as
spectrum.
z = exp(s/j,), the magnitude
and
ECG ,z
SYSTEM
IDENTIFICATION
VIA
411
CZT
0.40
5
: 015 n; -0.10 .E $ -0.35 3 + g -0.60 a
I
0
100
200
0
I
100
Sample
No.
DCT Coefficient
Sample
No.
DCT Coefficient
I
L
200
No.
No.
Normalized
Frequency
No.
FIG. 1. (a) Premature Ventricular Contraction (PVC) signal. (b) DCT of signal in (a). (c) Magnitude spectrum of DCT in (b). (d) Reconstructed ECG signal from the model. (e) Reconstructed DCT signal from the model. (f) Magnitude spectrum of reconstructed DCT in (e).
The above method of estimating when the spectrum is noisy.
Y and 8 is used only for the T-wave pole,
3. EXAMPLES 3.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 coprocessor. The quality of the reconstructed signal is assessed by the measure normalized mean square error (NMSE), defined in (10). The clinical utility of the method is illustrated by the examples given below followed by the summary on analysis of a large number of ECG data. 3.2. Example
1
A Premature Ventricular Contraction (PVC) beat of length 256 samples (N) and its DCT are shown in Figs. la and lb. The complex cepstrum and 512-point CZT of the DCT indicated a complex zero outside the unit circle at Y = 1.4 and
412
NIKANJAN
AND MLIKTH\I
TAHLt POLES.
ZEROS,
4ND
i l%Ah
F\
Pole<
Zeros Magnitude I
Phase (4 (in degrees)
PVC in Fig. la
0.1 0.91 I .4 0.87 0.94
0 k 20 + 36 k 64.93 I80
0.90 0.92 0.89 0.74
LBBB in Fig. 2b
1.2 1.3 0.78 I .53 1.3
0 k 41.83 f 60.78 t 69.2 180
0.87 0.94 0.87 0.67
Signal
SAMPI
Magnitude I
Phax (4 (in degrees)
Acrual peah sample
t’\tlmattxi peak sample 4x 60 96 IX0
+ ~
31,
-.-
t 64.98 2 78.23 f 122.8
il.7 92 II? 180
36 92 III I74
t3 = ?36”. Using the complex cepstrum technique the response of the outside zero was deconvolved from the signal DCT to obtain the minimum phase component (MPC). The magnitude spectrum of this MPC computed along a semicircle of radius 0.92 showed sharp peaks at angles 1~42”. The phase spectrum evaluated along radial lines between 38” and 46” also showed a sharp change of about 25~rad between dampings 0.92 and 0.93, confirming the presence of a pole. The response due to this identified pole was deconvolved from the MPC by inverse filtering. This procedure was repeated until all the poles and zeros were identified. The reconstructed DCT and signal with an NMSE of 10.2% are shown in Figs. le and Id. The identified poles and zeros are given in Table 1. The pole-zero plot is shown in Fig. 3a. Now converting the pole angles (0) 33.75”, 42”, 67.5”, and 127” to sample numbers by the inverse relation (N - 1)8/180, we obtain 48,60,96, and 180. Correlating these sample numbers with Fig. la, we see that they correspond to the R wave, the S wave, the kink in the QRS after the S wave, and the peak in the T wave, respectively. By a similar analysis we see that the zeros at 0”. 20”, and 180” correspond to the isoelectric regions preceding the QRS and after the T wave, while those at 36 and 64.93” correspond to the points in the QRS at which the signal crosses the baseline. The actual peak samples and those estimated from pole locations are given in Table 1. The maximum deviation between the two in most cases is around ?2 samples. Also the magnitude spectra of the DCT of the original and reconstructed signal, shown in Figs. lc and If, respectively, resemble the rectified and smoothed ECG signal.
ECG SYSTEM 320
IDENTIFICATION
VIA CZT
413
a
160 270
320
d
160
320 e 160 0 -160 -320
k 0
100
200
FIG. 2. Sample No. vs Amplitude. (a) Simulated ECG data obtained by Steiglitz-McBride model for the signal in (b). (b) Left Bundle Branch Block (LBBB) beat. (c) Normal ECG signal. (d) Reconstructed signal by the CZT model for the simulated signal in (a). (e) Reconstructed signal by the model for the actual signal in (b). (f) Reconstructed signal by the model for the normal ECG in (c).
3.3. Example 2
The DCT of a Left Bundle Branch Block (LBBB) beat, shown in Fig. 2b is modeled by the Steiglitz-McBride (SM) method (2). The IDCT of the identified model impulse response shown in Fig. 2a is used as the simulated ECG data for analysis. For this case as well as for other simulated data, the identification of model singularities by the CZT technique was found to be relatively easy. The reconstructed signal with an NMSE of 3.8% is shown in Fig. 2d. The model was then estimated using the DCT of the actual signal shown in Fig. 2b. However, in this case the magnitude spectrum of the CZT was noisy and did not show very narrow peaks or valleys. After the larger magnitude poles and zeros were identified and deconvolved from the signal, the T-wave pole was identified from the parameters of the parabolic fit made to the dominant part of the magnitude spectrum. The reconstructed signal is shown in Fig. 2e. The identified model singularities and the peak sample numbers are given in Table 1, while the pole-zero plot is shown in Fig. 3b. Due to the noisy spectrum, the identified parameters were a little inaccurate, which resulted in a higher NMSE of 14.8%. Next consider a normal ECG signal shown in Fig. 2c. All the model zeros for this signal were lying inside the unit circle and thus the model was identified with relative ease. The reconstructed signal is shown in Fig. 2f.
414
-15
NIRANJAN
-1 0
-0 5
FIG. 3. CZT Fig. 2(b).
0.0
pole-zero
0.5
10
plot [CX) pole.
AND
1.5
:MI.Kl
-.
(0) zero]
L. 3
lik
*
i -(
0
--AL,i---05
for (a) PVC signal
..L. OS
00
of Fig.
I(a).
~... ..-i.-+.-.. ‘0
rhr LBBB
‘5
heat of
3.4. Example 3 In this example we have considered a continuous strip of ECG of length 500 samples containing two beats, as shown in Fig. 4a. The poles and zeros identified from the signal CZT, shown in Fig. 4b, are distributed in the Z plane according to the locations of the component waves in the time signal. The poles with largest magnitude in the clusters corresponding to two QRS complexes are separated by an angle of 85”. This corresponds to a separation of 236 samples between the respective R-wave peaks in the time signal, which is the actual z c3 x 0 2 b 3c a -E g
04
b
a
r-pole.
0.2
o-zero
1 t
0 -0.2
0
500 samole
-0.1
-0.05 difference
no
0
0.05
in amplitudes
0.1
0 -10
-5 difference
0
5
?O
in intervals
FIG. 4. Continuous strip of ECG signal containing two beats. I-) Original signal. t”.) Reconstructed signal. tb) Pole-zero plot [(x) pole, (0) zero] for the DCT ofthe signal in (a). (c) A histogram representing the differences in amplitudes of samples in the 200 actual and reconstructed ECG beats, analyzed. td) A histogram representing the differences in actual and measured R-R interval. for 200 ECG beats.
ECG SYSTEM IDENTIFICATION
VIACZT
415
R-R interval for the signal shown. When signal frames of smaller durations with single R waves are considered, the angular separation between predominant pole clusters can be carried forward to the next frame, to determine the R-R interval. Continuous strips of ECG data from our own and MIT databases containing 200 beats have been analyzed. A histogram representing the difference in amplitudes between the actual and the reconstructed signals is shown in Fig. 4c. In order to make a uniform comparison, the original signal was normalized to have amplitudes between 0 and I, and the reconstructed signal was normalized by the same scale factor. As can be seen from the histogram, the maximum difference in amplitudes is 50.075, while most of the differences are less than kO.016. Another histogram representing the discrepancy in the actual and measured R-R interval for 200 ECG beats is shown in Fig. 4d. The maximum difference in R-R interval is r6 samples; however, we do not see any specific pattern in the distribution of the differences from Fig. 4d. 4. DISCUSSIONS The frequency transformation property of DCT allowed demarcation of the isoelectric region, baseline crossing points and peaks in the component waves from the angles subtended by the zeros and poles in the Z plane. Thus the polar plot can be visualized as a representation of significant points in the signal and thus can be made use of for the delineation of component waves. As the QRS usually requires two or more pairs of model poles for adequate representation, it is identified in a given frame of signal by searching for pole clusters. The angular separation between the QRS pole clusters is converted to time information and thereby heart rate. Identification of T-wave pole, which is usually located deep inside the unit circle, posed some problems when the spectrum was noisy. In such cases the pole was estimated from the parameters of a leastsquares parabolic fit made to the dominant part of the spectrum. On an average a signal of length 250 samples can be represented by about six complex conjugate or its equivalent real pole-zero pairs. This results in a data compression of a little more than 1 in 10. The reconstructed signals from the model parameters are of good quality with NMSEs less than 15%. On average the CZT was computed about 12 to 15 times to identify all the model singularities, depending upon the complexity of the ECG signal. However, it must be noted that an erroneously identified pole or zero will introduce its counterpart in the signal while inverse filtering. 5. CONCLUSIONS
A simple method based on repetitive scanning of the Z plane of the discrete cosine transformed ECG signal by the chirp Z transform technique is proposed for the identification of poles and zeros. We demonstrated the utility of the method by relating the model parameters to the features in the waveform.
3Ih
NIRANJAN
AND
MUKI
tti
Specifically, we obtained the R-R interval and locations ot component and isoelectric regions, in addition to a data compression of 1 in IO.
/. 2. 3. -1. 5. 6. 7. 8. Y. IO.
tiave\
MURTHY. I. S. N., RANG4RU. M. R., UDUPA. K. J.. ANDCJOYAI ,\. K. Homom~)rphlc andlyb~\ and modelling of ECG hignats. lEE.E Trur1.r. Biotnr