- Email: [email protected]
Computer analysis of rest and exercise electrocardiograms
COMPUTERS
AND
Computer
BIOMEDICAL
Analysis
RESEARCH
5, 329-346 (1972)
of Rest and
Exercise
Electrocardiograms*
H. K. WOLF, P. J. MACINNIS,
S. STOCK, R. K. HELPPI, AND P. M. RAUTAHARJU
Medical Biophysics rind Bioengineering Research Laboratory, Dalhousie University, Halifax, Nova Scotia, Canada Received November 19, 1971
The Dalhousie Program imposes rigid performance specifications and employs some new techniques. It can accommodate the three Frank orthogonal leads or six bipolar leads, which it computationally combines into three orthogonal leads of electrocardiograms (ECG). Each individual component lead is subjected to quality control, including noise abatement, thus providing pure signals. If the root-mean-square value of noise is less than 60 WV then in 95 % of all measurements the measurement precision is f10 msec for QRS duration, +20 msec for other intervals, and 1-25 yV for amplitudes. Patterns of all QRS complexes found are analyzed, using clustering techniques; and the members of the homogeneous majority cluster are averaged, using cross correlation for alignment of the complexes. The result is considered as the representative complex. The diagnostic program of Pipberger is used for clinical interpretation of the rest ECG. A modified diagnostic program is used in epidemiology studies, based on combined utilization of ECG and clinical information. The diagnostic program for interpretation of the exercise ECG is based on vector analysis of the ST-segment waveform.
Despite encouraging results with initial attempts (Z-3) and subsequent refinements (47), the progress in computer analysis of electrocardiograms (ECG) in the past decade has been rather slow. Only recently have attempts been made to use advanced statistical techniques for ECG interpretation (8). It has been very difficult to develop computer programs that produced results comparable to traditional visual measurement and interpretation. The medical profession has also been very reluctant to accept computer programs or systems if they merely duplicated human performance. Acceptance has been more favourable in exercise electrocardiography where the computer provides improved measurement and the interpretation of parameters which lie beyond the reach of visual determination (9-13). * Supported by the Government of Canada Public Health Research Grant (Project 602-7-I 33), the Nova Scotia Heart Foundation, and the Medical Research Council of Canada. Copyright 0 1972 by Academic Press, Inc. 329 All rights of reproduction in any form reserved.
330
ET AI..
WOLF
The Dalhousie ECG-analysis program for processing rest and exercise ECG employs some new concepts in computer analysis of the statistical properties of the ECG signal and noise, using optimising principles to achieve a satisfactory performance. PROGRAM OBJECTIVES
Formulation of the program provide reliable analysis of both 1. Measurement precision: msec; amplitude measurements,
took into account the following requirements, to rest and exercise ECG : QRS duration , -10 msec; all other intervals, ~20 525 PV or 5 % of the signal amplitude, whichever
DETERMINE NOISE LEVEL IN ALL LEADS OF RESTING ECG OR EXERCISE ECG
LOCATE All
MEASURE 1. TIME
ORS
-w FROM
INTERVAL
COMPLEXES
? EACH COMPLEX R;Rn-,
2. AMPLITUDE AT FIDUCIAL 3. QRS DURATION
ASSEMBLE OR5 COMPLEXES AND SELECT MAJORITY
3
PGlNT
INTO CLUSTERS CLUSTER
4
AVERAGE THE MMBERS OF SELECKG CLUSTER
5
WAVE RECOWlTlON ON AVERAGED COMPLEX
6
DIAGNOSIS BASED ON RESTING OR EXERCISE ECG
i
FIG. 1. Flowchart of the Dalhousie ECG Program, indicating the sequence of the seven major processing steps.
is greater. This measurement precision must be maintained in an environment of up to 60 PV root-mean-square of noise (band width 30-120 Hz), and measurements must be within the specified limits in at least 95% of the records. Measurement precision must be verified initially with specified simulation procedures (24). 2. Noise-abatement procedures to be used only to the extent required to maintain measurement precision. 3. Measurements to be only those required for quality control or the diagnostic programs. 4. The program must accommodate three or six simultaneously recorded ECG leads. 5. The program incorporates only limited arrhythmia analysis; therefore, classification is restricted to diagnostic categories that can be made reliably on the
COMPUTERECGANALYSIS
basis of QRS duration analysis of the P wave.
and R-R interval
331
measurements and only rudimentary
Consideration of these five points resulted in the program sequence depicted in Fig. 1. Perhaps the most obvious departure from the logic followed in other ECG programs is the scheme of selecting the representative complex. Detected QRS complexes are clustered according to three simple measurements made from each complex. The majority cluster is selected, and the complex closest to its center is used for comprehensive wave-form analysis. However, if noise level of the record exceeds a pre-set threshold, the members of the majority cluster are averaged and the resulting averaged complex is used. Thus, only one complex is analyzed for diagnostic interpretation. The advantages are greater precision and consistency of measurements than with most other programs, and substantial saving of computer time. MEASUREMENTPROGRAMDESIGN
Noise and Duration Parameters
The first 1.5 set of each ECG lead is used for noise measurement and characterization of the expected signal amplitudes. This period of each of the 3 or 6 simultaneous leads is scanned to locate a QRS complex. Here, as throughout the program, the search is made not from component ECGs but from a derived function, the filtered spatial velocity (FSV). The FSV function is defined in Appendix A. The existence of a QRS complex is established if the FSV exceeds and remains above threshold for a specified length of time, the threshold level being 3% of maxima1 FSV in the 1.5 set sample. When a QRS complex has been positively identified and its beginning (JON)and end (JOFF)have been estimated within the first 1.5 set of the signal, the following steps are taken : 1. The scalar lead with the largest peak amplitude (AMAX) in the QRS complex is selected. 2. The time from JONto peak amplitude in this selected lead, as a fraction of the estimated QRS duration, is defined as the relative tracking interval (RTI).
RTI = (Jmnx - JONMJOFF - JON)
(1)
The 0.3 set segment immediately after the end of the first QRS complex is used to calculate the power spectrum for each scalar lead. Since this part of the ECG contains only the ST-T segment, and sometimes the P wave of the next complex, signal components are predominantly below 20 Hz (Fig. 2); and components in this segment that are above 30 Hz are therefore considered as noise. Based on this logic, three noise parameters are determined : (1) rms amplitude (2) rms amplitude 12
of the 50 Hz components (N50) of the 60 Hz components (N60)
332
WOLF
(3) Total noise (NT) 120 Hz.
ET AL.
as the rms value of the components
between 30 and
120
NT=
J
Z: i:30
(2)
ci2,
where Ci = amplitude of i Hz component. Comparison of the relative magnitude of the three parameters N50, N60 and N7‘ indicates whether noiseis predominantly random or is 50 or 60 Hz interference, and a noise-reduction scheme is selected accordingly. Since noise parameters are determined separately for each scalar lead, any necessary digital filtering can be 1.0
I 1II1
5Hz
25Hz
FREQUENCY
(Hz,
FIG. 2. Power spectrum of typical P, QRS, and T waves, computed from a relatively noise-free X lead. The scaling factor chosen gives the DC component of the T wave an amplitude of I. The 5- and 25-Hz lines are the lower and upper cut-off frequencies for the band-pass filter used for QRS detection.
limited to the lead in which the noise tolerance limit is exceeded,thereby minimizing filtering and resultant signal distortion. QRS Detection
We scan an auxiliary derived function, FSV, to locate QRS-complexes. Figure 3 showsthis function. It also demonstratesthat P and T waves are almost eliminated and the amount of random noise is significantly lessin the FSV than in the original X, Y and Z signals. The FSV threshold value for the detection of a QRS complex is set to 3 % of the peak value in the first 1.5 set sampled. A wave form is accepted as a valid QRS complex if the FSV stays above the threshold for longer than 50 msecand exceeds 20 % of maximal FSV for longer than 16 msec.These conditions guarantee reliable QRS detection in most cases.
COMPUTER
ECG ANALYSIS
333
3. “Noisy” X, Y and 2 signals of an ECG containing 60 PV RMS random noise. Estimated onset and offset are based on FSV function threshold. The three lines bisecting each QRS complex are: (a) QRS onset; (b) tracking amplitude; (c) estimated QRS offset. Because of digital filter delay, FSV is delayed 15 msec; this delay is corrected in the program. FIG. QRS
When a QRS complex has been located and estimates of onset (JON) and offset (JOFF)have been established, two measurementsare made: 1. The time during which the FSV stays above the threshold (estimated QRS duration); QRS, = (JOFF- JON)x SINT, (3) where SINT = Sample Interval. 2. Tracking amplitude (TA), the difference between the amplitude at the relative tracking point and that at the point where FSV first exceeds threshold (estimated QRS onset). TA = IDATA,,,,,,
- IDATAJON,isEr,
(4)
k = JON+ RTI x (JoFF - JON) ISEL = number of the selectedlead IDATA
= sampled ECG data.
Thesetwo measurements,which are determined in the selectedlead only, together with the interval from the preceding QRScomplex, characterize QRScomplexes and thereby reveal diagnostically important categories. TA especially is a sensitive
334
WOLF
ET AL.
indicator of changes in the QRS waveform partly because it is measured close to peak amplitude, at least in those complexes similar to the sample complex of the first 1.5 set from the record. The greatest changes in amplitude occur both sides of the relative tracking point. A change in QRS pattern, or an error in determining the relative tracking point, gives rise to a TA greatly different from that of other complexes. In Fig. 3, for example, X is the selected lead; therefore, TA is close to peak amplitude. The 15msec delay between ECG and FSV caused by the bandpass filter is corrected during estimation of wave onsets and offsets. The parameters measured on each QRS complex (the points at which FSV Hurst rises above the threshold, the point at which it is last above it, and TA) were selected not for their diagnostic significance but to provide minimal necessary characterization of the types of QRS complexes and indicate unambiguously the location of each. These threshold-crossing points appear to be good parameters for categorizing QRS complexes by their morphologic similarity. Selection qf Representative Complex (a) Cluster formation. It is reasonable to assume that two major types of irregularities due to biological variation may occur in the ECG: the so-called respiratory sinus arrhythmia with associated amplitude and interval changes; and irregular beats, which originate from one or more ectopic pacemaker locations (extrasystoles). The latter may have an abnormal wave configuration and can occur at any phase of respiration. Thus, an ECG may have waveform patterns which occur predominantly during inspiration or expiration, assuming that the transition phase is relatively short. If significant variation occurs, the QRScomplexes in any record can becategorized into distinctly different groups. Since it is logical to choose the representative complex from the members of the majority group, an algorithm is used that clusters the QRS complexes into similarity groups (16). All but the first QRS complex are grouped, each complex being characterized by three scaled coordinates-estimated QRS duration, TA, and the interval to the preceding complex. For scaling, the most common value is set at 1 for the first two parameters and at 0.5 for the third; thus, the relative deviations permitted for the R-R interval are twice those permitted for QRS duration and tracking amplitude. Figure 4 shows an ECG record and the result of clustering the population of 23 QRS complexes of this record. An ellipse with axes of 1 SD is drawn around the center of each cluster, and the location ofeach QRScomplex in the three-dimensional space is indicated by a number defining the cluster to which the complex belongs. The two premature atria1 beats in the source record, identified by arrows, are in cluster 3. The slight respirogenic variation in QRS pattern results in formation of clusters 1 and 2. This example illustrates how clustering can improve the homogeneity of the population available for averaging. It is obvious that the scatter of complexes is reduced, if only cluster 1 is utilized for averaging rather than including the whole
COMPUTER
R-R
INTERVRL
ECG
335
ANALYSIS
TRRCKING
AMPLJTUOE
660
ORS-IJJRRTION
ORS-DURATION
FIG. 4. The 23 QRS complexes (the first complex is not considered) of the ECG record shown in the upper part of this figure are assembled into three clusters. The two complexes marked by arrows are premature atria1 beats and appear in cluster 3. The slight respirogenic pattern variation which exists in this record is sufficient to separate the remaining complexes into clusters 1 and 2. Variation of the three parameters is considerably reduced for each of the clusters as compared with the total population. Ellipses with axes of 1 SD are drawn around the center of each cluster. The units for QRS duration and R-R interval are in milliseconds, the axis for tracking amplitude is scaled in arbitrary units.
,O,P
7 lJF CC~NPLEXC'. ,'01 ilifll ,'>i HtlcIRE
mm
% OF COWLEXES . 468 AESTING RECORDS
FIG. 5. Distribution of the percentage of QRS complexes classified into the majority cluster in analysis of (A) 201 exercise, and (B) 461 rest records, each of 15-set duration.
336
WOLF
ET AL.
population for the calculation of the representative complex. Clustering not only removes excessive physiologic ECG variation; it also improves measurement reliability, by eliminating complexes that cause measurement errors in waverecognition programs. This effect is shown in Fig. 5, which depicts clustered QRS complexes of rest and exercise records. In the rest ECGs there is fairly wide scatter in the percentage of complexes included in the majority cluster. Since the physiologic variation is much smaller in an exercise ECG than in a rest ECG, one would expect the percentage of complexes in the majority cluster to be greater in the former, but this is not so, a finding attributed to the greater noise in them. Figure 6 demonstrates the “onset estimate” based on the threshold-crossing point. It shows the frequency distribution of the differences in this estimate between complexes in one cluster. The complex closest to the center of each cluster was taken
FIG. 6. Distribution of the difference between the QRS onset estimate of the center complex and of all other complexes of the majority cluster, after alignment by cross correlation. The histogram was computed from analysis of 120clusters totalling 1114 QRS complexes.
as the standard for the “true” pattern of QRS onset and offset: and the difference in onset estimate between each intracluster complex and the reference complex was calculated after aligning the complexes to maximize covariance. This difference was within f4 msec in 98 % of cases. The distribution skewness toward negative values is the result of noise in the ECG, which favors premature estimation of onset. (b) Final complex selection. When all of the QRS complexes in an ECG have been assembled into similarity clusters, we select the specific complex from which the diagnostically important measurements will be made. One obvious choice is the majority-cluster member closest to the cluster center and, therefore, closest to the mean of that cluster. This selection is preferred because it requires the least amount of computer time and, unlike most averaging methods, does not distort the ECG signal. We accept the center complex of the majority cluster as the representative complex whenever random noise in the record is -90 PV RMS. When this condition is not met, some or all members of the cluster are used to calculate a representative complex, by time-coherent averaging. Since noise is rarely
COMPUTER
ECG
ANALYSIS
337
synchronous with the ECG signal, this method of averaging reduces random interference proportional to the square root of the number of complexes used for averaging. With this type of noise-abatement procedure (17) the predominant problems are errors in smoothing for QRS onset, and, to a lesser degree, QRS offset. Although these errors potentially distort the signal, this disadvantage is usually compensated by greater measurement precision than is possible when noise has not been reduced. Clustering is done before averaging, thereby obviating errors due to averaging of inhomogeneous populations. The adverse effects of averaging on the signal are minimized if the complexes are aligned carefully. We use a cross-correlation method, with the cluster’s center complex as a reference. The complexes are aligned according to their QRS-onset estimates, and then are shifted left or right until covariance with the reference is
FIG. 7. Distribution of time difference in Q&S-complex alignment, using cross correlation and the QRS peak as a reference point, on 774 complexes which were members of the majority cluster of each record.
maximal; thus, all parts of a complex contribute to the final alignment. However, due to the basically triangular shape of the QRS complex, amplitudes around the peak play a dominant role in the establishment of alignment, so that in most cases alignment according to the location of peak amplitude is almost as effective as the cross-correlation method. It should be noted however that the remarkably accurate result depicted in Fig. 7 (error within &I2 msecin 95 %) was achieved on complexes which belonged to the same cluster. The larger errors occurred when the signal contained 60 Hz noise or the complex had two peaks of almost equal amplitude: in such cases, peak amplitude may occur at different points of the QRS complex, depending on the phase relationship between the 60 Hz signal and the ECG or on physiologic factors such as respiration. Figure 8 depicts the alignment procedure; only the selectedlead is plotted, since only that lead is used for alignment. The averaged complex (center.foot) illustrates the improvement in signal/noise ratio, especially in comparison with complex 10.
338
WOLF
ET AL.
Plots of the cross-covariancefunction are calculated both in coarse stepsand in tiner stepsaround the peak. The average standard deviation (SD) for the entire complex is used by some(18) to judge the quality of alignment, but results obtained with this
ms
so IS o.oYoo
FIG. 8. Alignment by cross correlation. Center: 10 complexes of a majority cluster, above the average of these after alignment. Lefi: plots of cross-correlation as a function of displacement for (fur Zefi) coarse and (inside left) fine shifts. Right: enlarged QAS-complex, between the range of +SD for each sample point. Average SD of the entire complex is 2.8 % of peak amplitude.
method are lessprecise: such factors asnoisecontent of the original signal and slight differences in the QRS patterns have a considerable effect on the calculation of SD. However, if correction is made for noise, the measuremay be useful in determining the acceptability of the averaging procedure: for example, large vahresindicate that
COMPUTER
ECG ANALYSIS
339
complexes with large variations in their basic pattern have been included in the averaged cluster. Thus, the average SD reflects more the quality of the clustering than the accuracy of the averaging procedure. Baseline fluctuations in the original ECG, resulting from poor instrumentation or loose electrodes, can introduce amplitude measurement errors which, due to the nature of the error signal, are impossible to correct. In order to safeguard against these errors the complexes of the selected cluster are investigated for baseline drift. A reference amplitude is calculated for each complex and each lead as the average of a 100 msec segment starting 20 msec after the end of the QRS complex. The range of fluctuation of this reference level is taken as an indicator for baseline drift. No attempt is made to correct any detected drift, but in cases of severe baseline fluctuations the diagnostic statements include a warning about possible erroneous amplitude measurements. Ana1ysi.r of the Selected Complex (a) Determination qfzero-reference line. This is one of the most crucial preliminary steps. The data for zero baseline selection are the amplitudes of the X, Y and 2 leads of all samples within a 100 msec window preceding the estimated QRS onset. The zero baseline is chosen as the point where the activities in X, Y and 2 are a minimum, or in other words, where there is little difference between successive sample points. A very convenient way to establish that point is given in the clustering method described above. The clusters formed by successive sample points in the X, Y and Z space are investigated. Initially, each sample is taken as a cluster center, The clustering radius R is increased in 2 PV steps until significant clusters begin to appear. As soon as one cluster contains all the samples of a 16 msec segment of the record, the process is being terminated and the baseline (X0, Y,.Z,) is defined as the average of the cluster members.
Xi, Yi, Zi are the sample points of the Nj members of the,j-th cluster. If there are several similar clusters within certain limits, the one closest in time to the Q/W-onset estimate is used to define the baseline point. This procedure for identification of the baseline is more complicated than other schemes. During the development of the Dalhousie program we found that approaches using simple criteria were adequate to identify the baseline correctly in only 99.4% of records, and that improvement beyond this level required a different method. (b) QRS complex. At this stage the data base is broadened by adding the spatial
340
WOLF
ET AL.
magnitude function (SMF), which is computed as the square root of the sum of the squared X, Y and 2 amplitudes, and the difference function (DF) computed from the SMF. While X, Y and Z are still unprocessed signals, the SMF is passed through a 60-Hz zero-phase shift low-pass filter, and its time derivative (DF) is processed by a 40-Hz zero-phase shift low-pass filter (Fig. 9).
FIG. 9. X, Yand Zcomponents of the selected complex (average of IO complexes of the majority cluster). The spatial magnitude SMF on the bottom of the graph has been processed by a 60 Hz low-pass filter and the time-derivative of SMF, and DF function has been passed through a 40-E filter. The calibration lines on the left represent 1 mV. The time marks on the graph indicate, from left to right, P onset, P offset, baseline reference point, QRS onset, QRS offset, T maximum and T offset.
The onset and offset points of the QRS complex are established in two steps. First, the initial estimate is calculated on the basis of the difference function, a step that yields acceptable results in most cases.Then the X, Y and Z leads are investigated for possible correction of these initial estimates: their slope and curvature are computed on the basis of a second-order least-square fit around onset and offset points. The conditions for the correct onset and offset points are that in at least one lead
COMPUTERECGANALYSIS
341
the slope has to be a certain preset value and that the slope in all other leads is less. The preset values are : for onset ]slope1= 2 pV/msec for offset [slope1= 5 pV/msec This second step usually improves the initial measurements. The QRS complex amplitudes are measured after a second-order fit to the sampled data, a procedure that permits measurements between sample points without reliance on linear interpolation. (Interpolation of this kind is needed, for instance, in determination of the time-normalized samplepoints like the l/8 vectors used in Pipberger’s program, and becomes very important when the data are sampledwith a low sampling rate.) (c) T UYXXZ. The offset point of the T wave is computed solely on the basisof DF. First, the wave’s maximal spatial magnitude is determined at the point where DF crosseszero in negative direction, the searchfor this point being limited-to that part of DF which starts 30 msec past the QRS offset and ends TX msec past the QRS offset. TX is a heart rate (HR) dependent variable and is calculated as
TX=383-
1.2 x HR,
(6)
where TX is expressedin msecand HR in beats/min. We previously establishedthe correlation between TX and HR on the basisof 1200 rest and exerciseECGs. The T offset point is the point where DF approaches the zero line again, after the zero crossing at maximum of the T-wave. (d) P waoe.The beginning and end of the P wave are calculated, using the SMF for P wave. For the onset, the amplitude and the slope of SMF are considered, to guarantee proper detection even when the P wave starts on the trailing part of the preceding T wave. For the P offset point, only the amplitude of SMF is taken into account. DIAGNOSTIC
PROCEDURES
The diagnostic program for clinical interpretation of rest ECGs is the unmodified multivariate analysis program formulated by Pipberger (8), and the Dalhousie program incorporates also, as an option, a specific diagnostic program designed for epidemiologic studies. This program is intended as a replacement of the so-called Minnesota code (19) for ECG classification. The most important diagnostic categories distinguished in our program for population studiesare : 1. “Hard-core” infarction criteria based on direction and on relative magnitudes of the initial QRS vectors. The ECG counterpart of this classification scheme is the Minnesota code category 1.1 (major Q waves). 2. Infarction criteria derived from QRS measurementswith lower diagnostic specificity than those in category 1, combined with T-wave abnormalities. This
342
WOLF
ET AL.
category replaces combinations of the Minnesota code categories 1.2 and 5. I, or 1.I! and 5.2 (lesser Q waves and major T-wave findings), and is expected to provide supportive evidence when the history indicates the likelihood of possible myocardial infarction. 3. Left ventricular hypertrophy. This diagnosis is based on age-specific amplitude criteria. These are derived from the components of the Frank leads and are thought to reflect the magnitude of the maximal left ventricular excitation forces. The diagnostic program for analysis of exercise ECGs (10) is primarily based on vector analysis of the ST-T segment waveform. The lead with the most ischemic ST segment is selected from all leads within a 45” reference cone in direction 1. I. I (left lower anterior octant), and measurements from this lead are used to classify the ECG response. The term ischemic is used in this context to characterize STsegment depression below the baseline and is measured as the time integral of the negative ST amplitudes in the selected lead. A similar analysis procedure is performed on the ten selected bipolar leads. In addition, the lead with a maximum STsegment elevation is identified and the mean amplitude of the ST segment elevation is measured. To characterize T-wave abnormalities, the integral of the T wave is measured from the lead where this integral is smallest (or most negative). Normal limits for all of these ST-T segment variables have been determined from large groups of male subjects examined in the course of longitudinal follow-up studies. The scope of the rhythm analysis routine has been limited to detection of gross abnormalities in order to keep the required length of ECG record to a minimum. The search for rhythm irregularities is based primarily on the analysis of the interbeat intervals. The P-wave in the representative (averaged) complex is taken as an indicator for atria1 activity. In cases of extrasystolic beats, QRS duration and tracking amplitude are utilized as wave form parameters to distinguish between ventricular and supraventricular extrasystoles. Appendix B describes the decision logic for the rhythm analysis, in the form of a flowchart.
KESULTS
The program has been tested for more than a year in cardiovascular epidemiological applications, and its performance was recently evaluated by visual analysis of 10,473 successive ECGs recorded in field conditions, 6,275 on resting patients and 4,198 during exercise. The program failed in QRS-ST measurements in 259 (2.6 %) and because of excessive noise in the record in 130 (1.3 %) cases; QRS-onset detection error >lO msec was observed in 15 records (0.2%), and QRS-offset detection errors >15 msec in 45 (0.4 YJ, and the end of the T wave was missed (error >30 msec) in 69 (0.7 %). f-wave detection failed in 180 (1.8 %) : the P wave was missed, or was detected with an error >20 msec at onset or offset.
COMPUTER
ECG ANALYSIS
343
One continuing weakness of the program is P-wave detection. P-wave measurement in terms of precision and noise tolerance is excellent, failing in less than 2 % of the records, which indicates that the program performs well when P waves are present. Since nearly all of the failures occurred in cases without P waves or where minor irregularities or noise caused false detection of P waves, it is evident that additional steps must be included in the detection logic for a more reliable arrhythmia analysis. The program has been implemented on an XDS Sigma 5 computer, which permits overlap of data acquisition and signal processing. Thus, since QRS detection is carried out in parallel with data acquisition, the locations of all of the QRScomplexes have been established by the time a 15set period of signals has been sampled. Execution of the rest of the program, including diagnostic interpretation takes a further 5 sec. DISCUSSION
One characteristic of the ECG signal is large inter-individual variability in waveforms, a fact that has to be taken into account during the design of an automatic ECG measurement program. Furthermore, self-adjusting procedures have to be incorporated, to permit adaptation of wave-detection algorithms to the specific wave pattern in each record. The determination of “case-specific parameters” (noise and duration) and the incorporation of clustering procedures have facilitated the program development and have improved its performance. Computer measurement of ECGs is further complicated by varying amounts of interference from noncardiac signal sources, such as electromyograms from skeletal muscle activity, electrode-skin interface artifacts, and 60- or 50-Hz interference from power lines. Although there are several methods by which reasonably reliable measurements can be obtained even when the ECG signal is partially obscured by noise, there is usually a substantial overlap of the spectrum of the noise and the ECG signal; therefore any reduction of noise may adversely affect the signal. Since an appropriate noise-reduction procedure can be chosen only with precise knowledge of the magnitude and type of noise in any given record, we measure noise separately for each component ECG before initiating measurement and analysis procedures. Perhaps the most important aspect in which the Dalhousie program differs from others is the selection of the representative ECG complex for measurement and analysis on the basis of clustering of all complexes. There is substantial beat-to-beat variation in the ECG pattern: part is true biological variation in the signal (for instance, due to respiration) and part results from measurement errors due to noise. In some ECG programs the beat-to-beat variation is ignored and measurements are made from one complex only; in others, measurements are made from several complexes and representative values are computed as arithmetic means or medians. The problems inherent with the method of single-complex measurement are obvious when one considers that QRS amplitude changes of 20 % and more are possible as a
344
WOLF
ET AL.
result of respiratory activity. On the other hand, when the mean value of measurements from more than one complex is used, some complexes are selected from the expiratory and others from the inspiratory phase of the respiration cycle, and the calculated average is biased by the number of complexes included and their phase relationship with the respiratory cycle; the resulting average does not necessarily represent even approximately any true wave form in the signal. These problems are avoided to a large extent when the clustering procedure is used. Extrasystoles in general are different from other complexes in at least one of the three clustering parameters used-QRS duration, preceding R-R interval and tracking amplitude. Furthermore, it is reasonable to assume that the ECG amplitude and respiration-induced interval changes follow a sinusoidal or trapezoidal time course. Consequently, complexes that originate during the inspiratory peak will form a relatively uniform group, as will the complexes originating during maximal expiration, and the remaining complexes will form groups with fewer members. The primary purpose of clustering is to identify one or two prominent groups from the population of QRS complexes; it is not critical to which cluster some of the intermediate complexes are classified. Numerous approaches have been proposed foi clustering. Since we are not interested in the optimal separation of all complexes into different groups, we can accept a clustering algorithm which is noniterative and therefore fast, a criterion satisfied by MacQueen’s procedure (f6) with only minor modifications. All investigators engaged in ECG program development have found it relatively easy to achieve a satisfactory measurement performance with even very primitive algorithms in a large majority (>95 %) of records but that any improvement beyond this level is increasingly more difficult. An example of this was the development of criteria for identification of the baseline: the classical approach, based on the use of various fixed threshold criteria, failed in 49 (0.62 %) of 7852 records, and improvement beyond this level became possible only with a more flexible, sophisticated approach, using clustering and dynamic analysis of the characteristics of each individual record rather than utilization of fixed criteria. No baseline failures were detected in the (admittedly fairly small) test sample when the new approach was implemented. APPENDIX
A : THE FUNCTION OF FSV
The function of FSV is calculated from forward differences of the prefiltered 3 or 6 leads. The prefiltering of the individual leads is performed by a digital bandpass filter, similar to the one described by Weaver (15). The lower cut-off frequency of the filter is set to 5 Hz in order to remove the P-wave and T-wave components as is shown in Fig. 2. The upper cut-off frequency was set to 25 Hz, to reduce line frequency interference and disturbance by other high frequency noise.
COMPUTER
ECG
ANALYSIS
345
Hence FSV, = 2 [IDATA:, - IDATA;+4, j]2 j=l N = number of leads (either 3 or 6) i = sample point number IDATA:, = Data point IDATAi,j after processing by bandpass filter. The FSV function therefore represents primarily the squared mode of the spatial velocity of the QRS complex, the P- and T-wave parts being removed by prefiltering. APPENDIX
B: THE DECISION LOGIC FOR RHYTHM ANALYSIS
FIG. 10. The flowchart depicts the decision tree for the rhythm analysis. The algorithm utilizes as source information: (a) the R-R intervals, (b) the QRS durations, (c) the tracking amplitudes, (d) the grouping of QRS-complexes (number of clusters), and (e) the parameters of the P wave in the averaged complex.
REFERENCES F. W. AND PIPBERGER, H. V. Automatic regconition of electrocardiographic waves by digital computer. Civc. Res. 9, 1138-1143 (1961). CADY, L. D., JR., WOODBURY, M. A., TICK, L. J., AND GERTLER, M. M. A method for electrocardiogram wave-pattern estimation. Example: left ventricular hypertrophy. Circ. Res. 9, 1078-1082 (1961). CACERES, C. A., STEINBERG, C. A., ABRAHAM, S., CARBERY, W. T., MCBRIDE, T. M., TOLLES, W. E., AND RIKLI, A. E. Computer extraction of electrocardiographic parameters. Circulation 25,356-362 (1962).
1. STALLMANN, 2.
3.
346
WOLF
ET AL.
H. V. Computer analysis of electrocardiograms. In “Computers in Biomedical Research” (B. D. Waxman and R. W. Stacy, Eds.), Vol. I. Chap. 16, pp. 377-406, Academic Press, New York, 1965. 5. SMITH, R. E. AND HYDE, C. M. A computer system for electrocardiographic analysis. Proc. &d Ann. Rocky Mountain Bio-Engineering Symp., University of Colorado, Boulder, CO May Z--3. 4. PIPBERGER,
1966. 6. OKAJIMA,
M., STARK, L., WHIPPLE, G. H., AND YASUI, S. Computer pattern recognition techniques: Some results with real electrocardiographic data. IEEE Trans. Bio-med. Electronics 10, 106-114 (1963). 7. BONNER, R. E. AND SCHWETMAN, H. D. Computer diagnosis of electrocardiograms Il. A computer program for ECG measurements. Comput. Biomed. Rex 1,366-386 (1968). 8. KLINGEMAN, J. AND PIPBERGER, H. V. Computer classifications of electrocardiograms. Compnt. Biomed. Res. 1, 1-17 (1967). 9. BLOMQVIST, G. The Frank lead exercise electrocardiogram. A quantitative study based on averaging technic and digital computer analysis. Acra Med. Scund. 178 (Suppl. 440), 9-39 (1965). 10. RAUTAHARJU, P. M. AND WOLF, H. K. Computer interpretation and classification of exercise electrocardiograms, vector cardiographic aspects. In “X&h International Symposium on Vectorcardiography, New York” (I. Hoffman, Ed.). North-Holland, Amsterdam, 1971. If. SHEFFIELD, L. T., HOLT, J. H., LESTER, F. M., CONROY,D. V., AND REEVES,T. J. On-line analysis of the exercise electrocardiogram. Circulation 40, 935-944 (1969). 12. VON DER GROEBEN,T., TOOLE, T. G., WEAVER, C. S., AND FITZGERALD, J. W. Noise reduction in exercise electrocardiograms by digital filter techniques. In “Measurements in Exercise Electrocardiography” (H. Blackburn, Ed.). C. C. Thomas Publishers, Springfield. IL. 1969. 13. RAUTAHARJU, P. M. AND BLACKBURN, H. The exercise electrocardiogram: Experience in analysis of “noisy” cardiograms with a small computer. Amer. Heart J. 69,515-520 (1965). 14. WOLF, H. AND RAUTAHARJU,P. M. Assessment of precision and accuracy of computer programs for analysis of electrocardiograms. In “Measurements in Exercise Electrocardiography” (H. Blackburn, Ed.), pp. 169-180. C. C. Thomas, Springfield, IL. 1969. 1.5. WEAVER, C. S., VON DER GROEBEN, T., MANTEY, P. E., TOOLE, T. G., COLE, C. A., FI-IZGERALU. T. W., AND LAWRENCE,R. W. Digital filtering with application to electrocardiogram processing. IEEE Trans. Audio and Electroacoustic. AV-16 (1968). 16. MACQUEEN, T. Some methods for classification and analysis of multivariate analysis. In “Proceedings of Fifth Berkley Symposium on Mathematical Statistics and Probability.” University of California Press, Berkeley, CA, 1966. 17. RAUTAHARJU, P. M. Deterministic type waveform analysis in electrocardiography. Ann. I%‘.I)‘. Acad. Sri. 128,939-954 (1966). 18. BRANDON, C. W. AND BRODY, D. A. A hardware trigger for temporal indexing of the electrocardiographic signal. Comput. Biomed. Res. 3,45-57 (1970). 19. BLACKBURN, H., KEYES, A., SIMONSON,E., RAUTAHARJU, P. M., ANL) PUNSAR, S. The electrocardiogram in population studies. A classification system. Circuhtion 21, 1160-l 175 (1960).
AND
Computer
BIOMEDICAL
Analysis
RESEARCH
5, 329-346 (1972)
of Rest and
Exercise
Electrocardiograms*
H. K. WOLF, P. J. MACINNIS,
S. STOCK, R. K. HELPPI, AND P. M. RAUTAHARJU
Medical Biophysics rind Bioengineering Research Laboratory, Dalhousie University, Halifax, Nova Scotia, Canada Received November 19, 1971
The Dalhousie Program imposes rigid performance specifications and employs some new techniques. It can accommodate the three Frank orthogonal leads or six bipolar leads, which it computationally combines into three orthogonal leads of electrocardiograms (ECG). Each individual component lead is subjected to quality control, including noise abatement, thus providing pure signals. If the root-mean-square value of noise is less than 60 WV then in 95 % of all measurements the measurement precision is f10 msec for QRS duration, +20 msec for other intervals, and 1-25 yV for amplitudes. Patterns of all QRS complexes found are analyzed, using clustering techniques; and the members of the homogeneous majority cluster are averaged, using cross correlation for alignment of the complexes. The result is considered as the representative complex. The diagnostic program of Pipberger is used for clinical interpretation of the rest ECG. A modified diagnostic program is used in epidemiology studies, based on combined utilization of ECG and clinical information. The diagnostic program for interpretation of the exercise ECG is based on vector analysis of the ST-segment waveform.
Despite encouraging results with initial attempts (Z-3) and subsequent refinements (47), the progress in computer analysis of electrocardiograms (ECG) in the past decade has been rather slow. Only recently have attempts been made to use advanced statistical techniques for ECG interpretation (8). It has been very difficult to develop computer programs that produced results comparable to traditional visual measurement and interpretation. The medical profession has also been very reluctant to accept computer programs or systems if they merely duplicated human performance. Acceptance has been more favourable in exercise electrocardiography where the computer provides improved measurement and the interpretation of parameters which lie beyond the reach of visual determination (9-13). * Supported by the Government of Canada Public Health Research Grant (Project 602-7-I 33), the Nova Scotia Heart Foundation, and the Medical Research Council of Canada. Copyright 0 1972 by Academic Press, Inc. 329 All rights of reproduction in any form reserved.
330
ET AI..
WOLF
The Dalhousie ECG-analysis program for processing rest and exercise ECG employs some new concepts in computer analysis of the statistical properties of the ECG signal and noise, using optimising principles to achieve a satisfactory performance. PROGRAM OBJECTIVES
Formulation of the program provide reliable analysis of both 1. Measurement precision: msec; amplitude measurements,
took into account the following requirements, to rest and exercise ECG : QRS duration , -10 msec; all other intervals, ~20 525 PV or 5 % of the signal amplitude, whichever
DETERMINE NOISE LEVEL IN ALL LEADS OF RESTING ECG OR EXERCISE ECG
LOCATE All
MEASURE 1. TIME
ORS
-w FROM
INTERVAL
COMPLEXES
? EACH COMPLEX R;Rn-,
2. AMPLITUDE AT FIDUCIAL 3. QRS DURATION
ASSEMBLE OR5 COMPLEXES AND SELECT MAJORITY
3
PGlNT
INTO CLUSTERS CLUSTER
4
AVERAGE THE MMBERS OF SELECKG CLUSTER
5
WAVE RECOWlTlON ON AVERAGED COMPLEX
6
DIAGNOSIS BASED ON RESTING OR EXERCISE ECG
i
FIG. 1. Flowchart of the Dalhousie ECG Program, indicating the sequence of the seven major processing steps.
is greater. This measurement precision must be maintained in an environment of up to 60 PV root-mean-square of noise (band width 30-120 Hz), and measurements must be within the specified limits in at least 95% of the records. Measurement precision must be verified initially with specified simulation procedures (24). 2. Noise-abatement procedures to be used only to the extent required to maintain measurement precision. 3. Measurements to be only those required for quality control or the diagnostic programs. 4. The program must accommodate three or six simultaneously recorded ECG leads. 5. The program incorporates only limited arrhythmia analysis; therefore, classification is restricted to diagnostic categories that can be made reliably on the
COMPUTERECGANALYSIS
basis of QRS duration analysis of the P wave.
and R-R interval
331
measurements and only rudimentary
Consideration of these five points resulted in the program sequence depicted in Fig. 1. Perhaps the most obvious departure from the logic followed in other ECG programs is the scheme of selecting the representative complex. Detected QRS complexes are clustered according to three simple measurements made from each complex. The majority cluster is selected, and the complex closest to its center is used for comprehensive wave-form analysis. However, if noise level of the record exceeds a pre-set threshold, the members of the majority cluster are averaged and the resulting averaged complex is used. Thus, only one complex is analyzed for diagnostic interpretation. The advantages are greater precision and consistency of measurements than with most other programs, and substantial saving of computer time. MEASUREMENTPROGRAMDESIGN
Noise and Duration Parameters
The first 1.5 set of each ECG lead is used for noise measurement and characterization of the expected signal amplitudes. This period of each of the 3 or 6 simultaneous leads is scanned to locate a QRS complex. Here, as throughout the program, the search is made not from component ECGs but from a derived function, the filtered spatial velocity (FSV). The FSV function is defined in Appendix A. The existence of a QRS complex is established if the FSV exceeds and remains above threshold for a specified length of time, the threshold level being 3% of maxima1 FSV in the 1.5 set sample. When a QRS complex has been positively identified and its beginning (JON)and end (JOFF)have been estimated within the first 1.5 set of the signal, the following steps are taken : 1. The scalar lead with the largest peak amplitude (AMAX) in the QRS complex is selected. 2. The time from JONto peak amplitude in this selected lead, as a fraction of the estimated QRS duration, is defined as the relative tracking interval (RTI).
RTI = (Jmnx - JONMJOFF - JON)
(1)
The 0.3 set segment immediately after the end of the first QRS complex is used to calculate the power spectrum for each scalar lead. Since this part of the ECG contains only the ST-T segment, and sometimes the P wave of the next complex, signal components are predominantly below 20 Hz (Fig. 2); and components in this segment that are above 30 Hz are therefore considered as noise. Based on this logic, three noise parameters are determined : (1) rms amplitude (2) rms amplitude 12
of the 50 Hz components (N50) of the 60 Hz components (N60)
332
WOLF
(3) Total noise (NT) 120 Hz.
ET AL.
as the rms value of the components
between 30 and
120
NT=
J
Z: i:30
(2)
ci2,
where Ci = amplitude of i Hz component. Comparison of the relative magnitude of the three parameters N50, N60 and N7‘ indicates whether noiseis predominantly random or is 50 or 60 Hz interference, and a noise-reduction scheme is selected accordingly. Since noise parameters are determined separately for each scalar lead, any necessary digital filtering can be 1.0
I 1II1
5Hz
25Hz
FREQUENCY
(Hz,
FIG. 2. Power spectrum of typical P, QRS, and T waves, computed from a relatively noise-free X lead. The scaling factor chosen gives the DC component of the T wave an amplitude of I. The 5- and 25-Hz lines are the lower and upper cut-off frequencies for the band-pass filter used for QRS detection.
limited to the lead in which the noise tolerance limit is exceeded,thereby minimizing filtering and resultant signal distortion. QRS Detection
We scan an auxiliary derived function, FSV, to locate QRS-complexes. Figure 3 showsthis function. It also demonstratesthat P and T waves are almost eliminated and the amount of random noise is significantly lessin the FSV than in the original X, Y and Z signals. The FSV threshold value for the detection of a QRS complex is set to 3 % of the peak value in the first 1.5 set sampled. A wave form is accepted as a valid QRS complex if the FSV stays above the threshold for longer than 50 msecand exceeds 20 % of maximal FSV for longer than 16 msec.These conditions guarantee reliable QRS detection in most cases.
COMPUTER
ECG ANALYSIS
333
3. “Noisy” X, Y and 2 signals of an ECG containing 60 PV RMS random noise. Estimated onset and offset are based on FSV function threshold. The three lines bisecting each QRS complex are: (a) QRS onset; (b) tracking amplitude; (c) estimated QRS offset. Because of digital filter delay, FSV is delayed 15 msec; this delay is corrected in the program. FIG. QRS
When a QRS complex has been located and estimates of onset (JON) and offset (JOFF)have been established, two measurementsare made: 1. The time during which the FSV stays above the threshold (estimated QRS duration); QRS, = (JOFF- JON)x SINT, (3) where SINT = Sample Interval. 2. Tracking amplitude (TA), the difference between the amplitude at the relative tracking point and that at the point where FSV first exceeds threshold (estimated QRS onset). TA = IDATA,,,,,,
- IDATAJON,isEr,
(4)
k = JON+ RTI x (JoFF - JON) ISEL = number of the selectedlead IDATA
= sampled ECG data.
Thesetwo measurements,which are determined in the selectedlead only, together with the interval from the preceding QRScomplex, characterize QRScomplexes and thereby reveal diagnostically important categories. TA especially is a sensitive
334
WOLF
ET AL.
indicator of changes in the QRS waveform partly because it is measured close to peak amplitude, at least in those complexes similar to the sample complex of the first 1.5 set from the record. The greatest changes in amplitude occur both sides of the relative tracking point. A change in QRS pattern, or an error in determining the relative tracking point, gives rise to a TA greatly different from that of other complexes. In Fig. 3, for example, X is the selected lead; therefore, TA is close to peak amplitude. The 15msec delay between ECG and FSV caused by the bandpass filter is corrected during estimation of wave onsets and offsets. The parameters measured on each QRS complex (the points at which FSV Hurst rises above the threshold, the point at which it is last above it, and TA) were selected not for their diagnostic significance but to provide minimal necessary characterization of the types of QRS complexes and indicate unambiguously the location of each. These threshold-crossing points appear to be good parameters for categorizing QRS complexes by their morphologic similarity. Selection qf Representative Complex (a) Cluster formation. It is reasonable to assume that two major types of irregularities due to biological variation may occur in the ECG: the so-called respiratory sinus arrhythmia with associated amplitude and interval changes; and irregular beats, which originate from one or more ectopic pacemaker locations (extrasystoles). The latter may have an abnormal wave configuration and can occur at any phase of respiration. Thus, an ECG may have waveform patterns which occur predominantly during inspiration or expiration, assuming that the transition phase is relatively short. If significant variation occurs, the QRScomplexes in any record can becategorized into distinctly different groups. Since it is logical to choose the representative complex from the members of the majority group, an algorithm is used that clusters the QRS complexes into similarity groups (16). All but the first QRS complex are grouped, each complex being characterized by three scaled coordinates-estimated QRS duration, TA, and the interval to the preceding complex. For scaling, the most common value is set at 1 for the first two parameters and at 0.5 for the third; thus, the relative deviations permitted for the R-R interval are twice those permitted for QRS duration and tracking amplitude. Figure 4 shows an ECG record and the result of clustering the population of 23 QRS complexes of this record. An ellipse with axes of 1 SD is drawn around the center of each cluster, and the location ofeach QRScomplex in the three-dimensional space is indicated by a number defining the cluster to which the complex belongs. The two premature atria1 beats in the source record, identified by arrows, are in cluster 3. The slight respirogenic variation in QRS pattern results in formation of clusters 1 and 2. This example illustrates how clustering can improve the homogeneity of the population available for averaging. It is obvious that the scatter of complexes is reduced, if only cluster 1 is utilized for averaging rather than including the whole
COMPUTER
R-R
INTERVRL
ECG
335
ANALYSIS
TRRCKING
AMPLJTUOE
660
ORS-IJJRRTION
ORS-DURATION
FIG. 4. The 23 QRS complexes (the first complex is not considered) of the ECG record shown in the upper part of this figure are assembled into three clusters. The two complexes marked by arrows are premature atria1 beats and appear in cluster 3. The slight respirogenic pattern variation which exists in this record is sufficient to separate the remaining complexes into clusters 1 and 2. Variation of the three parameters is considerably reduced for each of the clusters as compared with the total population. Ellipses with axes of 1 SD are drawn around the center of each cluster. The units for QRS duration and R-R interval are in milliseconds, the axis for tracking amplitude is scaled in arbitrary units.
,O,P
7 lJF CC~NPLEXC'. ,'01 ilifll ,'>i HtlcIRE
mm
% OF COWLEXES . 468 AESTING RECORDS
FIG. 5. Distribution of the percentage of QRS complexes classified into the majority cluster in analysis of (A) 201 exercise, and (B) 461 rest records, each of 15-set duration.
336
WOLF
ET AL.
population for the calculation of the representative complex. Clustering not only removes excessive physiologic ECG variation; it also improves measurement reliability, by eliminating complexes that cause measurement errors in waverecognition programs. This effect is shown in Fig. 5, which depicts clustered QRS complexes of rest and exercise records. In the rest ECGs there is fairly wide scatter in the percentage of complexes included in the majority cluster. Since the physiologic variation is much smaller in an exercise ECG than in a rest ECG, one would expect the percentage of complexes in the majority cluster to be greater in the former, but this is not so, a finding attributed to the greater noise in them. Figure 6 demonstrates the “onset estimate” based on the threshold-crossing point. It shows the frequency distribution of the differences in this estimate between complexes in one cluster. The complex closest to the center of each cluster was taken
FIG. 6. Distribution of the difference between the QRS onset estimate of the center complex and of all other complexes of the majority cluster, after alignment by cross correlation. The histogram was computed from analysis of 120clusters totalling 1114 QRS complexes.
as the standard for the “true” pattern of QRS onset and offset: and the difference in onset estimate between each intracluster complex and the reference complex was calculated after aligning the complexes to maximize covariance. This difference was within f4 msec in 98 % of cases. The distribution skewness toward negative values is the result of noise in the ECG, which favors premature estimation of onset. (b) Final complex selection. When all of the QRS complexes in an ECG have been assembled into similarity clusters, we select the specific complex from which the diagnostically important measurements will be made. One obvious choice is the majority-cluster member closest to the cluster center and, therefore, closest to the mean of that cluster. This selection is preferred because it requires the least amount of computer time and, unlike most averaging methods, does not distort the ECG signal. We accept the center complex of the majority cluster as the representative complex whenever random noise in the record is -90 PV RMS. When this condition is not met, some or all members of the cluster are used to calculate a representative complex, by time-coherent averaging. Since noise is rarely
COMPUTER
ECG
ANALYSIS
337
synchronous with the ECG signal, this method of averaging reduces random interference proportional to the square root of the number of complexes used for averaging. With this type of noise-abatement procedure (17) the predominant problems are errors in smoothing for QRS onset, and, to a lesser degree, QRS offset. Although these errors potentially distort the signal, this disadvantage is usually compensated by greater measurement precision than is possible when noise has not been reduced. Clustering is done before averaging, thereby obviating errors due to averaging of inhomogeneous populations. The adverse effects of averaging on the signal are minimized if the complexes are aligned carefully. We use a cross-correlation method, with the cluster’s center complex as a reference. The complexes are aligned according to their QRS-onset estimates, and then are shifted left or right until covariance with the reference is
FIG. 7. Distribution of time difference in Q&S-complex alignment, using cross correlation and the QRS peak as a reference point, on 774 complexes which were members of the majority cluster of each record.
maximal; thus, all parts of a complex contribute to the final alignment. However, due to the basically triangular shape of the QRS complex, amplitudes around the peak play a dominant role in the establishment of alignment, so that in most cases alignment according to the location of peak amplitude is almost as effective as the cross-correlation method. It should be noted however that the remarkably accurate result depicted in Fig. 7 (error within &I2 msecin 95 %) was achieved on complexes which belonged to the same cluster. The larger errors occurred when the signal contained 60 Hz noise or the complex had two peaks of almost equal amplitude: in such cases, peak amplitude may occur at different points of the QRS complex, depending on the phase relationship between the 60 Hz signal and the ECG or on physiologic factors such as respiration. Figure 8 depicts the alignment procedure; only the selectedlead is plotted, since only that lead is used for alignment. The averaged complex (center.foot) illustrates the improvement in signal/noise ratio, especially in comparison with complex 10.
338
WOLF
ET AL.
Plots of the cross-covariancefunction are calculated both in coarse stepsand in tiner stepsaround the peak. The average standard deviation (SD) for the entire complex is used by some(18) to judge the quality of alignment, but results obtained with this
ms
so IS o.oYoo
FIG. 8. Alignment by cross correlation. Center: 10 complexes of a majority cluster, above the average of these after alignment. Lefi: plots of cross-correlation as a function of displacement for (fur Zefi) coarse and (inside left) fine shifts. Right: enlarged QAS-complex, between the range of +SD for each sample point. Average SD of the entire complex is 2.8 % of peak amplitude.
method are lessprecise: such factors asnoisecontent of the original signal and slight differences in the QRS patterns have a considerable effect on the calculation of SD. However, if correction is made for noise, the measuremay be useful in determining the acceptability of the averaging procedure: for example, large vahresindicate that
COMPUTER
ECG ANALYSIS
339
complexes with large variations in their basic pattern have been included in the averaged cluster. Thus, the average SD reflects more the quality of the clustering than the accuracy of the averaging procedure. Baseline fluctuations in the original ECG, resulting from poor instrumentation or loose electrodes, can introduce amplitude measurement errors which, due to the nature of the error signal, are impossible to correct. In order to safeguard against these errors the complexes of the selected cluster are investigated for baseline drift. A reference amplitude is calculated for each complex and each lead as the average of a 100 msec segment starting 20 msec after the end of the QRS complex. The range of fluctuation of this reference level is taken as an indicator for baseline drift. No attempt is made to correct any detected drift, but in cases of severe baseline fluctuations the diagnostic statements include a warning about possible erroneous amplitude measurements. Ana1ysi.r of the Selected Complex (a) Determination qfzero-reference line. This is one of the most crucial preliminary steps. The data for zero baseline selection are the amplitudes of the X, Y and 2 leads of all samples within a 100 msec window preceding the estimated QRS onset. The zero baseline is chosen as the point where the activities in X, Y and 2 are a minimum, or in other words, where there is little difference between successive sample points. A very convenient way to establish that point is given in the clustering method described above. The clusters formed by successive sample points in the X, Y and Z space are investigated. Initially, each sample is taken as a cluster center, The clustering radius R is increased in 2 PV steps until significant clusters begin to appear. As soon as one cluster contains all the samples of a 16 msec segment of the record, the process is being terminated and the baseline (X0, Y,.Z,) is defined as the average of the cluster members.
Xi, Yi, Zi are the sample points of the Nj members of the,j-th cluster. If there are several similar clusters within certain limits, the one closest in time to the Q/W-onset estimate is used to define the baseline point. This procedure for identification of the baseline is more complicated than other schemes. During the development of the Dalhousie program we found that approaches using simple criteria were adequate to identify the baseline correctly in only 99.4% of records, and that improvement beyond this level required a different method. (b) QRS complex. At this stage the data base is broadened by adding the spatial
340
WOLF
ET AL.
magnitude function (SMF), which is computed as the square root of the sum of the squared X, Y and 2 amplitudes, and the difference function (DF) computed from the SMF. While X, Y and Z are still unprocessed signals, the SMF is passed through a 60-Hz zero-phase shift low-pass filter, and its time derivative (DF) is processed by a 40-Hz zero-phase shift low-pass filter (Fig. 9).
FIG. 9. X, Yand Zcomponents of the selected complex (average of IO complexes of the majority cluster). The spatial magnitude SMF on the bottom of the graph has been processed by a 60 Hz low-pass filter and the time-derivative of SMF, and DF function has been passed through a 40-E filter. The calibration lines on the left represent 1 mV. The time marks on the graph indicate, from left to right, P onset, P offset, baseline reference point, QRS onset, QRS offset, T maximum and T offset.
The onset and offset points of the QRS complex are established in two steps. First, the initial estimate is calculated on the basis of the difference function, a step that yields acceptable results in most cases.Then the X, Y and Z leads are investigated for possible correction of these initial estimates: their slope and curvature are computed on the basis of a second-order least-square fit around onset and offset points. The conditions for the correct onset and offset points are that in at least one lead
COMPUTERECGANALYSIS
341
the slope has to be a certain preset value and that the slope in all other leads is less. The preset values are : for onset ]slope1= 2 pV/msec for offset [slope1= 5 pV/msec This second step usually improves the initial measurements. The QRS complex amplitudes are measured after a second-order fit to the sampled data, a procedure that permits measurements between sample points without reliance on linear interpolation. (Interpolation of this kind is needed, for instance, in determination of the time-normalized samplepoints like the l/8 vectors used in Pipberger’s program, and becomes very important when the data are sampledwith a low sampling rate.) (c) T UYXXZ. The offset point of the T wave is computed solely on the basisof DF. First, the wave’s maximal spatial magnitude is determined at the point where DF crosseszero in negative direction, the searchfor this point being limited-to that part of DF which starts 30 msec past the QRS offset and ends TX msec past the QRS offset. TX is a heart rate (HR) dependent variable and is calculated as
TX=383-
1.2 x HR,
(6)
where TX is expressedin msecand HR in beats/min. We previously establishedthe correlation between TX and HR on the basisof 1200 rest and exerciseECGs. The T offset point is the point where DF approaches the zero line again, after the zero crossing at maximum of the T-wave. (d) P waoe.The beginning and end of the P wave are calculated, using the SMF for P wave. For the onset, the amplitude and the slope of SMF are considered, to guarantee proper detection even when the P wave starts on the trailing part of the preceding T wave. For the P offset point, only the amplitude of SMF is taken into account. DIAGNOSTIC
PROCEDURES
The diagnostic program for clinical interpretation of rest ECGs is the unmodified multivariate analysis program formulated by Pipberger (8), and the Dalhousie program incorporates also, as an option, a specific diagnostic program designed for epidemiologic studies. This program is intended as a replacement of the so-called Minnesota code (19) for ECG classification. The most important diagnostic categories distinguished in our program for population studiesare : 1. “Hard-core” infarction criteria based on direction and on relative magnitudes of the initial QRS vectors. The ECG counterpart of this classification scheme is the Minnesota code category 1.1 (major Q waves). 2. Infarction criteria derived from QRS measurementswith lower diagnostic specificity than those in category 1, combined with T-wave abnormalities. This
342
WOLF
ET AL.
category replaces combinations of the Minnesota code categories 1.2 and 5. I, or 1.I! and 5.2 (lesser Q waves and major T-wave findings), and is expected to provide supportive evidence when the history indicates the likelihood of possible myocardial infarction. 3. Left ventricular hypertrophy. This diagnosis is based on age-specific amplitude criteria. These are derived from the components of the Frank leads and are thought to reflect the magnitude of the maximal left ventricular excitation forces. The diagnostic program for analysis of exercise ECGs (10) is primarily based on vector analysis of the ST-T segment waveform. The lead with the most ischemic ST segment is selected from all leads within a 45” reference cone in direction 1. I. I (left lower anterior octant), and measurements from this lead are used to classify the ECG response. The term ischemic is used in this context to characterize STsegment depression below the baseline and is measured as the time integral of the negative ST amplitudes in the selected lead. A similar analysis procedure is performed on the ten selected bipolar leads. In addition, the lead with a maximum STsegment elevation is identified and the mean amplitude of the ST segment elevation is measured. To characterize T-wave abnormalities, the integral of the T wave is measured from the lead where this integral is smallest (or most negative). Normal limits for all of these ST-T segment variables have been determined from large groups of male subjects examined in the course of longitudinal follow-up studies. The scope of the rhythm analysis routine has been limited to detection of gross abnormalities in order to keep the required length of ECG record to a minimum. The search for rhythm irregularities is based primarily on the analysis of the interbeat intervals. The P-wave in the representative (averaged) complex is taken as an indicator for atria1 activity. In cases of extrasystolic beats, QRS duration and tracking amplitude are utilized as wave form parameters to distinguish between ventricular and supraventricular extrasystoles. Appendix B describes the decision logic for the rhythm analysis, in the form of a flowchart.
KESULTS
The program has been tested for more than a year in cardiovascular epidemiological applications, and its performance was recently evaluated by visual analysis of 10,473 successive ECGs recorded in field conditions, 6,275 on resting patients and 4,198 during exercise. The program failed in QRS-ST measurements in 259 (2.6 %) and because of excessive noise in the record in 130 (1.3 %) cases; QRS-onset detection error >lO msec was observed in 15 records (0.2%), and QRS-offset detection errors >15 msec in 45 (0.4 YJ, and the end of the T wave was missed (error >30 msec) in 69 (0.7 %). f-wave detection failed in 180 (1.8 %) : the P wave was missed, or was detected with an error >20 msec at onset or offset.
COMPUTER
ECG ANALYSIS
343
One continuing weakness of the program is P-wave detection. P-wave measurement in terms of precision and noise tolerance is excellent, failing in less than 2 % of the records, which indicates that the program performs well when P waves are present. Since nearly all of the failures occurred in cases without P waves or where minor irregularities or noise caused false detection of P waves, it is evident that additional steps must be included in the detection logic for a more reliable arrhythmia analysis. The program has been implemented on an XDS Sigma 5 computer, which permits overlap of data acquisition and signal processing. Thus, since QRS detection is carried out in parallel with data acquisition, the locations of all of the QRScomplexes have been established by the time a 15set period of signals has been sampled. Execution of the rest of the program, including diagnostic interpretation takes a further 5 sec. DISCUSSION
One characteristic of the ECG signal is large inter-individual variability in waveforms, a fact that has to be taken into account during the design of an automatic ECG measurement program. Furthermore, self-adjusting procedures have to be incorporated, to permit adaptation of wave-detection algorithms to the specific wave pattern in each record. The determination of “case-specific parameters” (noise and duration) and the incorporation of clustering procedures have facilitated the program development and have improved its performance. Computer measurement of ECGs is further complicated by varying amounts of interference from noncardiac signal sources, such as electromyograms from skeletal muscle activity, electrode-skin interface artifacts, and 60- or 50-Hz interference from power lines. Although there are several methods by which reasonably reliable measurements can be obtained even when the ECG signal is partially obscured by noise, there is usually a substantial overlap of the spectrum of the noise and the ECG signal; therefore any reduction of noise may adversely affect the signal. Since an appropriate noise-reduction procedure can be chosen only with precise knowledge of the magnitude and type of noise in any given record, we measure noise separately for each component ECG before initiating measurement and analysis procedures. Perhaps the most important aspect in which the Dalhousie program differs from others is the selection of the representative ECG complex for measurement and analysis on the basis of clustering of all complexes. There is substantial beat-to-beat variation in the ECG pattern: part is true biological variation in the signal (for instance, due to respiration) and part results from measurement errors due to noise. In some ECG programs the beat-to-beat variation is ignored and measurements are made from one complex only; in others, measurements are made from several complexes and representative values are computed as arithmetic means or medians. The problems inherent with the method of single-complex measurement are obvious when one considers that QRS amplitude changes of 20 % and more are possible as a
344
WOLF
ET AL.
result of respiratory activity. On the other hand, when the mean value of measurements from more than one complex is used, some complexes are selected from the expiratory and others from the inspiratory phase of the respiration cycle, and the calculated average is biased by the number of complexes included and their phase relationship with the respiratory cycle; the resulting average does not necessarily represent even approximately any true wave form in the signal. These problems are avoided to a large extent when the clustering procedure is used. Extrasystoles in general are different from other complexes in at least one of the three clustering parameters used-QRS duration, preceding R-R interval and tracking amplitude. Furthermore, it is reasonable to assume that the ECG amplitude and respiration-induced interval changes follow a sinusoidal or trapezoidal time course. Consequently, complexes that originate during the inspiratory peak will form a relatively uniform group, as will the complexes originating during maximal expiration, and the remaining complexes will form groups with fewer members. The primary purpose of clustering is to identify one or two prominent groups from the population of QRS complexes; it is not critical to which cluster some of the intermediate complexes are classified. Numerous approaches have been proposed foi clustering. Since we are not interested in the optimal separation of all complexes into different groups, we can accept a clustering algorithm which is noniterative and therefore fast, a criterion satisfied by MacQueen’s procedure (f6) with only minor modifications. All investigators engaged in ECG program development have found it relatively easy to achieve a satisfactory measurement performance with even very primitive algorithms in a large majority (>95 %) of records but that any improvement beyond this level is increasingly more difficult. An example of this was the development of criteria for identification of the baseline: the classical approach, based on the use of various fixed threshold criteria, failed in 49 (0.62 %) of 7852 records, and improvement beyond this level became possible only with a more flexible, sophisticated approach, using clustering and dynamic analysis of the characteristics of each individual record rather than utilization of fixed criteria. No baseline failures were detected in the (admittedly fairly small) test sample when the new approach was implemented. APPENDIX
A : THE FUNCTION OF FSV
The function of FSV is calculated from forward differences of the prefiltered 3 or 6 leads. The prefiltering of the individual leads is performed by a digital bandpass filter, similar to the one described by Weaver (15). The lower cut-off frequency of the filter is set to 5 Hz in order to remove the P-wave and T-wave components as is shown in Fig. 2. The upper cut-off frequency was set to 25 Hz, to reduce line frequency interference and disturbance by other high frequency noise.
COMPUTER
ECG
ANALYSIS
345
Hence FSV, = 2 [IDATA:, - IDATA;+4, j]2 j=l N = number of leads (either 3 or 6) i = sample point number IDATA:, = Data point IDATAi,j after processing by bandpass filter. The FSV function therefore represents primarily the squared mode of the spatial velocity of the QRS complex, the P- and T-wave parts being removed by prefiltering. APPENDIX
B: THE DECISION LOGIC FOR RHYTHM ANALYSIS
FIG. 10. The flowchart depicts the decision tree for the rhythm analysis. The algorithm utilizes as source information: (a) the R-R intervals, (b) the QRS durations, (c) the tracking amplitudes, (d) the grouping of QRS-complexes (number of clusters), and (e) the parameters of the P wave in the averaged complex.
REFERENCES F. W. AND PIPBERGER, H. V. Automatic regconition of electrocardiographic waves by digital computer. Civc. Res. 9, 1138-1143 (1961). CADY, L. D., JR., WOODBURY, M. A., TICK, L. J., AND GERTLER, M. M. A method for electrocardiogram wave-pattern estimation. Example: left ventricular hypertrophy. Circ. Res. 9, 1078-1082 (1961). CACERES, C. A., STEINBERG, C. A., ABRAHAM, S., CARBERY, W. T., MCBRIDE, T. M., TOLLES, W. E., AND RIKLI, A. E. Computer extraction of electrocardiographic parameters. Circulation 25,356-362 (1962).
1. STALLMANN, 2.
3.
346
WOLF
ET AL.
H. V. Computer analysis of electrocardiograms. In “Computers in Biomedical Research” (B. D. Waxman and R. W. Stacy, Eds.), Vol. I. Chap. 16, pp. 377-406, Academic Press, New York, 1965. 5. SMITH, R. E. AND HYDE, C. M. A computer system for electrocardiographic analysis. Proc. &d Ann. Rocky Mountain Bio-Engineering Symp., University of Colorado, Boulder, CO May Z--3. 4. PIPBERGER,
1966. 6. OKAJIMA,
M., STARK, L., WHIPPLE, G. H., AND YASUI, S. Computer pattern recognition techniques: Some results with real electrocardiographic data. IEEE Trans. Bio-med. Electronics 10, 106-114 (1963). 7. BONNER, R. E. AND SCHWETMAN, H. D. Computer diagnosis of electrocardiograms Il. A computer program for ECG measurements. Comput. Biomed. Rex 1,366-386 (1968). 8. KLINGEMAN, J. AND PIPBERGER, H. V. Computer classifications of electrocardiograms. Compnt. Biomed. Res. 1, 1-17 (1967). 9. BLOMQVIST, G. The Frank lead exercise electrocardiogram. A quantitative study based on averaging technic and digital computer analysis. Acra Med. Scund. 178 (Suppl. 440), 9-39 (1965). 10. RAUTAHARJU, P. M. AND WOLF, H. K. Computer interpretation and classification of exercise electrocardiograms, vector cardiographic aspects. In “X&h International Symposium on Vectorcardiography, New York” (I. Hoffman, Ed.). North-Holland, Amsterdam, 1971. If. SHEFFIELD, L. T., HOLT, J. H., LESTER, F. M., CONROY,D. V., AND REEVES,T. J. On-line analysis of the exercise electrocardiogram. Circulation 40, 935-944 (1969). 12. VON DER GROEBEN,T., TOOLE, T. G., WEAVER, C. S., AND FITZGERALD, J. W. Noise reduction in exercise electrocardiograms by digital filter techniques. In “Measurements in Exercise Electrocardiography” (H. Blackburn, Ed.). C. C. Thomas Publishers, Springfield. IL. 1969. 13. RAUTAHARJU, P. M. AND BLACKBURN, H. The exercise electrocardiogram: Experience in analysis of “noisy” cardiograms with a small computer. Amer. Heart J. 69,515-520 (1965). 14. WOLF, H. AND RAUTAHARJU,P. M. Assessment of precision and accuracy of computer programs for analysis of electrocardiograms. In “Measurements in Exercise Electrocardiography” (H. Blackburn, Ed.), pp. 169-180. C. C. Thomas, Springfield, IL. 1969. 1.5. WEAVER, C. S., VON DER GROEBEN, T., MANTEY, P. E., TOOLE, T. G., COLE, C. A., FI-IZGERALU. T. W., AND LAWRENCE,R. W. Digital filtering with application to electrocardiogram processing. IEEE Trans. Audio and Electroacoustic. AV-16 (1968). 16. MACQUEEN, T. Some methods for classification and analysis of multivariate analysis. In “Proceedings of Fifth Berkley Symposium on Mathematical Statistics and Probability.” University of California Press, Berkeley, CA, 1966. 17. RAUTAHARJU, P. M. Deterministic type waveform analysis in electrocardiography. Ann. I%‘.I)‘. Acad. Sri. 128,939-954 (1966). 18. BRANDON, C. W. AND BRODY, D. A. A hardware trigger for temporal indexing of the electrocardiographic signal. Comput. Biomed. Res. 3,45-57 (1970). 19. BLACKBURN, H., KEYES, A., SIMONSON,E., RAUTAHARJU, P. M., ANL) PUNSAR, S. The electrocardiogram in population studies. A classification system. Circuhtion 21, 1160-l 175 (1960).