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