CA1256947A - Method and apparatus for analyzing electrocardiographic signals - Google Patents

Abstract

X9a ABSTRACT
METHOD AND APPARATUS FOR ANALYZING
ELECTROCARDIOGRAPHIC SIGNALS
ECG signals are digitized and a template representative of an arrhythmia free interval for a particular patient is formed. Subsequent ECG
signals which match the template are averaged to form averaged ECG signals to reduce noise. A fast Fourier transform (FFT) is performed on the terminal 40 millisecond portion of the QRS complex and the ST segment. The 60 db area and 40 Hz intercept of the resultant spectral output are used as a figure of merit in predicting the liklihood of a praticular patient experiencing ventricular tachycaxdia. Alternatively, the FFT magnitude is squared to form energy spectra of the ECG signals and a measure of the energy in a first preselected portion such as the 20 - 50 Hz region is compared with a second preselected portion such as the entire spectra to form a figure of merit for predicting ventricular tachycardia.

Description

X9a METHOD AND APPARATUS FOR ANALYZING
. ~
ELECTROCARD I OGRAPEII C S I GNALS

This inve~tion relates to electrocardio-graphy and, more particularly, to a method and ap~aratus or predicting risk for development of maligna~t ventricular arrh~thmias.
Many patients who have su~fered heart damage via a myocardial infarction are in danser of sudden death from acute arrhythmia. The ability to reliably and noninvasely predict the risk for development of such arrhy~hmias in patients i5 desirable.
Recently high gain amplification and signal processing techni~ues in the time domaln have detected low amplitude, hi~h frequency potentials in the terminal QRS complex and ST segments of signal averaged electrocardiograms (ECG's) obtained during arrhythmia free intervals from patients and exper1mental animals with sustained ventricular tachycardia. Recent studies of these patients wi~h time domain analysis have used a variety of low (25 to 100 Hz3 and high ~250 to 300 ~z) band pass filt2r. A ma~or limitation of these time domain procedure is the lack of a ~riori knowledge of the frequency distribution of the sig~als of interest 30 and hence the inherent risk that filtering will ~zs~
xga

-2-exclude sig~als of particular interest. Fast Fourier Transfonm (FFT) analysis is a powerful analytic method that i5 complimentary to time domain analysis and avoids some of the limitations of a prlorl filtering. Moreover, FFT analysis facilitates identification and characterization of .frequencies independent of regional amplitude a~d is thus particularly well suited for assessing low amplitude signals. Fi~ally, this system and FFT
a~alysis provide flexibility or analyzing different ECG region~ using the sam~ system hardwar0 and ~oftware.
As another esa~ple of time domain a~alysis, U.S. Patent 4,4 ~ ,459 (Simson~ discloses applying ~he latter portions of the X, Y, and Z digital QRS
signals in reverse time order to a digital high pass filter to eliminate a ringing artifact rom th~ filter output. The resulting filter~d outputs are combined to create ~ filtered QRS which is examined and ~he last 40 milliseconds of which is isolated and measured to obtain an indication of the level of high energy conte~t. ~he initial portion o~ the QRS wave form is also processed in a ~orward direction to obtain an indication of the QRS total duration.

This invention provides a noninvasive electrocardiographic system capable o improved objective identiication and characterization of low amplitude potentials in the surace ECG signal.
The X, Y, and Z ECG signals are amplified, recorded over a broad ban~width, converted from analog to digital, processed to select normal QRS
waveforms, and signal averaged to reduce the noise ~25~

X9a ~3--level. A fast Fourier transform is then performed on the terminal 40 milliseconds of the QRS complex and ST segme~ts of each of the avera~ed X, Y) and Z
ECG signals. The resultant spectral outputs are examined and a flgure of merit created to obtain a characterization of the frequency content of the various portions of the signal averayed ECG
exa~ined.
In one embodiment th~ magnitude of the FFT is first squared to form ~nergy spectra of the ECG
signals. An energy mea~ure o~ a first preselected portion ~uch as the 20 to 50 ~z region is compared with an energy measur~ o~ a second presele~ted portion such as the e~tire ~nergy spectra to determine the presence or absence of a predetermined frequency content in a predetermined portion of the ECG signals such as the combined QRS
terminal 40 milli~eco~ds and ST regio~. The comparison include~ taking ~he magnitudes of ~ha peaks in the first pre~elected portion of a spectrum and comparing them with ~he magnitude of the largest peak of ~he spectrum. Also the area under the spectrum curve in the first preselected portion is compared wi~h the axea of the entire spectrum.
The peak mag~itude and area comparisons after further operation are called the mag~itude ratio and area ratio. Since a plurality of ECG leads are used with each patient when monitorin~ the ECG, a mea~ value formed from all ECG leads for ~he magnitude ratio and area ratio is formed for each patient.
BRIEF DESC~IRTION OF THE DRAWINGS
FIG. 1 is a block diagram of a preferred embodiment of the present invention.

X9a FIG. 2A is a representative signal averaged ECG from one of the leads of FIG. 1 showing various segments of the ECG signal and the R-R interval.
FIG. 2B is a spectral plot of the FFT of the terminal Q~S s~gment of FIG. 2A.
FIG. 3 is a flow chart of the method used for .template formation.
FIG. 4 is a simplified flow chart of ~he program for selecting and averaging ECG signals.
FIG. 5 is a flow chart of ~he program used for chaxacterizing the frequency content of the averaged ECG si~nals generated by FIG. 4.
FIGS. 6A through 6F repres~nt the outputs of the program of FIG. S operating o~ a first s~gment of averaged X, Y and Z ECG signals from a patient with prior myocardial infarction without sustained ventricular tachycaxdia (VT) and a patient Wi~l prior myocardial infarction with sustained VT.
FIGS. 7A throug~ 7F represent the outputs of the program of FIG. 5 operating on a second se~m nt of averaged X, Y and Z ECG signals fxom a patient with prior myocaIdi~l infarction without sustained ventri~ular tachycardia (VT) and a patient with prior myocardial in~arction with sustained VT.
FIG. 8 is a flow chart of an alternate progr~m to ~hat of FIG. 5 used for characterizing the ~requency content of the averaged ECG signals generated by FIG. 4.
FIG. 9 is a representative signal averaged ~CG from one of the leads of FIG. 1 showing the combined terminal QRS and ST segment of the ECG
~ignal.
FIG. 10 is a sample energy spectra plot of the signal of FIG. 9 showing the main peak as a solid line on a first scale as marked on the left ~7 X9a -5- .

side and the secondary peaks as a dashed line on a second scale on the right sid~.
FIGS. llA through llF represent the outputs of the program of FIG. a operating on a first segment of averaged X, Y and Z ECG signals from a patient with prior myocardial infarction ~ithout sustained ventricular tachycardia (VT) and a patient with prior myocardial infarctlon with sustained VT.
FIGS. 12~ through 12F represent ~he outputs of the program of FIG. 8 operating on a seg~ent of averaged X, Y c~nd Z ECG sisnals ~rom a patient with prior myocardial infarction without sustained ventri~ular tachycardia (VT~ and a di~fererlt patient from the patient of FIGS. llA ~hrough llF
with pricr myocardial in~arction wlth sustained VT.
FIG. 13 is a ~raph comparing ~ha mean area ratio~ of patients from Groups I, II and III.
FI~. 14 is a graph comparing ~he mean magnitude ratios vs. frequency for patients from the Groups I, II and III.

Referring now to FIG. 1, a block diagram of the present in~ention designated generally 100 is shown including standard bipolar X, Y, and Z
electrocardiogram (ECG) leads 102, 104, and 106, respectively. Standard Frank X, Y, and Z leads, Model No. 1507-llA, are applied to a patient in accordance with the Model 1507-llA recommended electrode placement. Each of the X, Y, and Z lead signals are c~mplified ~1,000 fold) by Hewl~tt-Packard 1057 Amplifiers 110, 112, and 114, respectively. The outputs of the amplifiers are passed through band pass filters 116, 118, and 120, respectively, having a band~idth of 0.05 Hz to 470 X9a ~z. Then the signal~ are converted from analog to digltal signals by A~D converter 12Z at a 1 KHz rate having 12 bits of accuracy and a lowes~
re~olution o~ 1 micxovolt. An ADAC Inc. Model S No.16SE-C-3P6A-P A~D converter was used.
The input to the A/D converter 122 is switched sequentially among the outputs of filters 116, 118 and 120 with a 300 microsecond delay in betwee~. The digital X, Y a~d Z signals from A/D
con~erter 122 are provided in parallel over line 130 to microcomputer 150. A floppy di~k memory 152 and a Selanar raster graphics display 154 with rastex control are coupled to the microcomputer 150. A D~C VT 103 LSI 11/23 microcomputer with 64K
hytes of memory was used. The microcomputer 150 along with disk 152 and display 154 ca~ be mounted on a portable cart or use at the patient~ bedside.
The amplified X, Y, and Z signals are simultaneously displayed on ~he ~ewlett-Packard oscilliscope 160 for continuous real time visual, monitoring.
R~ferring now to FIG. 2A, a representative signal averaged Z ECG lead designated generally 200 is ~hown. Two QRS complexes designated generally 202 and 204 are shown. Tha large signal 206 is known as the R wave, and its peak 208 is called the fiducial point. The time separation between adjacent R waves 206 and 210 is called the R-R
interval. The peak to pcak amplitude of the QRS
signal is the voltaye m~asured from the positive peak of the R wave at 208 to the negative peak of the QRS waveform at 212.
Each ECG signal can be thought of as including: a terminal QRS signal 214 which comprises the last 40 msec of the QRS complex (end ~25~
xsa of QRS complex and preceding 40 msec); and ~n ST
segm~nt 216 which extends from ~he end of the QRS
segmen~ 214 to the beginning of the T segment 218 of the ECG signal 200.
Referring now to FIG. 3, a simplified block diagram for template formation designated generally 300 is shown for selecting certain ones o~ the digital X, Y, and Z signals which are to be con-sidered as arrhy~hmia-free sign ls and which are to be used in forming the signal averaged X, Y, and Z
ECG ignals for furth~r processing. By visual inspection of the display 160, a set of X, Y, and Z
signals is chosen which can be co~sidered to represent a display of normal sinu~ rhythm 302.
The ignal such as the o~e having the largest R
wave amplitude relative to the P and T points of the ECG signal is ~elected 304. The~ a three second sample of the X, Y and Z leads are digitized and stored 306. The preselected si~nal 304 is 20 displayed on display 154. Se~ 308. The R-R
interval (the distance between adjacent R waves of adjacent QRS complexes1 and the fiducial point (peak point 208 of the R wave) are identified using the ad~ustable cursor accompanying the Selanar display 310. At the same time the QRS ~mplitudes (peak to peak voltayes) for each of the X, Y, and Z
signals are determined 314. The R-R interval, fiducial point 208, the X, Y, and Z QRS amplitudes and 40 points of the R waves (20 points on either side of the fiducial point) are then stored as a template in the mini-computer 150 memory. This is done automatically using the DEC graphics software package in c~mhination with the use of the cursors on di~play 154.
Raferring now to FIG. 4 a process desig~ated ;6~a*7 xsa generally 400 for selecting arrhythmia free X, Yand Z signals and for averaging ECG signals is shown. Once again, a three second sample of the heartbeat is chosen, digitized, and the digitized X, Y, and Z signals for the three second interval are stored in a circular memory buffer. See 402.
The ~iddle beat of these stored three second sample is chosen 404 and the selected signal compared with the template. Specifically, the R=R
interval is compared with ~he template ~-R interval 406. If the R-R i~terval i~ not within *20% of the templat~ valu~ the beat i5 rejected. See 408.
O~herwise, the peak to p~ak amplitudes of the X, Y, and Z signals of the chosen beat are compared wi~h the QRS amplitudes of the template 410. If at least 2 of the 3 amplitudes are the same as the template amplitudes then ~he b~at is selected as a ca~didate for averaging, otherwise ~he b~at is rejected 412. Then a 40 point cross correlation 414 of the R wave is made wi~h the template wavefoxm about the fiducial point 208. The 40 points which are chosen are 20 points lying on either side of the fiducial point. If the correlation coefficient is not greater than 98% the beat is rejected 416. If the beat passes the amplitude and cross correlation comparisons, then for the current beat to be acceptable, the beat which went before it and the beat which follows it must all be found to be acceptable when compared with the R-R template. If th~ center b~at is found to be unacceptable the~ it is shifted left in the circular buffer and compared to the next beat.
This process goes on in real time a~d only if buffer overflow occurs do you go back to 404. The X, Y, and Z ECG signals of each acceptable beat are X9a a~eraged, point by point, with the X, Y, and Z ECG
signals of all other acceptable beats until a 100 beat average is created 418.
The above description of the operation of FIG. 4 is done automatically using the FORTRAN
program No. 1 (pages 20- 30).
Refexring now to FIG. 5, a fast Fourier transform (FFT) analysis 500 i5 performed on the terminal 40 milliseconds of the QR5 complex for each of the averaged X, Y, and Z signals. All of the averaged X, Y and Z QRS sisnals are stored in memo~y 152, e.q., the ~veraged X, Y and Z signals are display~d on display 15g. Se~ 502. Using the display cursors, the terminal 40 milliseconds of the QRS signal is determi~ed 504 by finding the end of the QRS complex and sliding 40 msec into ~he R
wave.
A fast Fourier trangform (FFT) 506 i9 performed on the chosen segment of th~ ECG signal, for example, on th~ ~ermi~al 40 millisecond of the QRS complex as de~in~d above. The FFT can be perform~d using any one of a plurality of acceptable commercially available computer programs suitable ~or use on the particulax minicomputer 150, for example, a FORTRAN program from the standard IEEE Library is acceptable and is hereto employed. For the particular segment o~ the ECG
si~nal chosen a 512 point fast Fourier transform was calculated. The selected sample values (those selected automatically by the cursor) were placed at the begi~ning of the 512 input array of the FFT
program and ~he remaining values set to z~ro. This stap permitted maintenance of the ~ame frequency scale in the output data but allowed an arbitrary mem~er o~ input values up to 512.

.,,.",, ' ;6~
xsa Fourier analysis assumes that the signal contained in the sample window interval is a repetitive function. If the initial sample point and final sample point are not isopotential, a sharp discontinuity will be introduced be~ween the end of one cycle and the beginning of the next that will artifactually add both high and low requencies to the original signal. To eliminate this source of error when perorming ~he FFT on discrete components of the ECG (i.~. the terminal 40 msec of the QR5 complex) the segment of intere~t is multiplied by a four term Blackman~arris window function 508 to s~oo~h the data to zero at th~
boundaries 504. 5ee "On the use of windows for harmonic analysis with the dis~reet Fourier trans-form", Harris, F.J., PROC IE~E 66:51,1978 and the FORTRAN Program #2 (pages 3~w40)~ The particular Blackman~Harris funGtion us~d herein had a 6 dB
bandwidth, a 92 dB sidelobe level, and a sidelobe falloff of 6 dB per octave. Multiplication of ~he Blackma~-~arris window must be performed be~ore the FFT is calculated.
The FFT data are plotted with a high resolution plotter (Versatec, Inc.) 510 and transferred to disk 152 for storage. Referring to FIG. 2B, for each plot, the dB drop at 40 Hz 250 is found and the area 252 under the curve from the fundamental frequency to the fre~uency at which the amplitude of the spectrum plot drops to 60 dB below the peak level (60 dB area) is determined. (See the ~ORTRAN program No. 3, pages 41-46). Thase two values together form a flgur~ of merit (FOM) 512, FIG. 5, for the spectral plot. The 40 Hz intercept . was chosen because most of the energy of a normal QRS is less than 35 Hz and because work by others ~2569~
xsa has shown khat fragmented signals have a peak frequency in the ~5 to 50 Hz range.
The above mentioned steps of Fig. 5 are pe.rfonmed simultaneously on the averaged X, Y and Z
signals and three spectral plots are formed. The mean of the 40 Hz intercept and the 60 dB areas of the FOM's for the averaged X, Y and Z signals is calculated for each patient to generate a single 40 ~z intercept and 60 dB area FOM. The process described i~ connection with FIG. 5 given ~bove was also p~rform~d or ST segm~nt 216 o~ the averaged X, ~ a~d Z ECG signals. It should be remembered that th~ ST segment 216 is de~ined as begi~ing at the end of the QRS complex and pxoceeding to the 15 beginning of ~he T wave. The ST segment is defined using the cursor of the Selanar display.
Using the a~ove-mentioned method, 61 patients were grouped according to cli m cal characteristic~.
None was currently receiving antiarrhythmic medications. Group I comprised 16 patients with prior myocardial infarction with at least one documented episode of sustained VT or cardiac arrest that was not associated wi~h a new infarct.
Each of the 13 of these 16 patients who.was studied in the Clinical Electrophysiology Laboratory at Washington University had inducible sustained VT or VF imilar to that occurri~ clinically. Group II
consisted of 35 patients with prior myocardial infarction without a history or documented episode of sustained VT (> 30 seconds in duration), syncope, or cardiac arrest who were admitted to the Barnes Hospital Telemetry Unit. All patients in this group were monitored for at least 7 days.
Seventeen exhibited absent or simple v@ntricular ectopy (Modified Lown Class 0 to 1); la patients ~L~5~ 7 xsa had complex ventricular ectopy (Class 2 to 5), of which nine patients had nonsustained VT. Group III
included 10 healthy male controls 24 to 40 yeaxs of age with no clinical evidence of organic heart disease or arrhy~hmias.
There were no significant differences between .Groups I a~d II wi~h regard to age, infarct location, presence of left ventricular aneurysms, or Q~S duration. Left ~entxicular ejection 10 fraction was signific~ntly less in patients with in~arction who had manife~ted sustained VT or VF
compared with ejection fraGtion in Group II
patient~ ~34% ~ 16 vs. 45% ~ 15, p ~ 0.02).
The FET analysis of FIG. S of signal averaged X, Y, and Z ECG signals showed signi~icant differences in the mean 60 dB area and the mean ~0 ~z intercept of the terminal 40 msec of the QRS and of the ST segment in ~atients with prior myoc~rdial infarction and a subsequent epi ode of sustained VT
or VF compared with values in patients wi~h prior myocardial infaxction without these arrhythmias and with values i~ normal subjects. There were no siy~ificant diff~rences i~ the 60 ~3 area or 40 Hz intercept of the terminal QRS or ST segment in patients with myocardial infarction without sustained VT compared with normal subjects.
Representative plots of FFT analysis of the terminal QRS complex from a patient with prior myocardial infarction without sustained VT and from a patient with prior myocardial infarction with sustained VT are shown in FIGS. 6 A-C and 6 D-F, respecti~ely. Each depicts power vs freguency plots of the terminal 40 msec of signal averaged QRS complexe~ recorded from bipolar X, Y, and Z
leads along with values for the 60 dB area and the i6~
xsa -13~

dB drop at the ~0 ~z intercept. In each lead, the terminal 40 msec of the QRS complex from the patient who had manifest sustained VT contained relatively more high frequency components ~han ~he S complex from the patient who had not, reflected by a greater value for the 60 dB area and a lesser dB
drop at the 40 ~z intercept. Similar plots were obtained for the FFT ana1ysis of the ST segments ~rom the same two patients. In e~ch lead, the ST
10 segment in t~e signal from the patient who had manifest an episode o~ sustained V~ contained relatively mo~e hiyh frequency compoll~nts than the ST segment from the patient without t7T. Spectral diffexences in the tenninal QP~S and ST seganent were 15 consistently most marked at frequencies less than 120 Hz.
Based on results of the FF~ analysis of FIG.
5 o~ the terminal 40 msec of the QRS complex in o~nal subjects (Group III~, mean 60 clB area values ~ s~reater than 2400 ancl me~n 40 ~z intercept values less than 47 dB were defined as abnormal and i~dicative of an incraase in high frequency compone~ts in the terminal QRS. Abnormal values for both th~ 60 dB area and 40 Hz intercept were fou~d in 88% of patie~ts with prior myocardial infarction having a subseguent episode of sustained V~ (Group I~ and 15% o patients with prior myocardial infarction wi~hout sustained VT (Group II). .
Based on results of the FFT analysis of FIG.
5 of the ST segment in Group III controls, mean 60 dB value~ greater than 2500 and maan 40 Hz intercept valuas les~ than 52 dB were defined as abnormal and indicative of increased high frequency components in the sr segment. Abnonmal values for 3L~5~r7 xsa the 60 dB area o the ST segment wexe found in 81%
and 25% of Group I and Group II patients, respectively. ~bno.rmal values for the 40 Hz intercept were found in 76% and 20~ respectively.
Values for ~oth the 60 dB area and 4~ Hz intercept of the terminal QRS and ST se~men~ were . independent of QRS duration, left ventricular ejection fraction, and complexity of spontaneous ventricular ecopty.
An lternat~ embodiment of the above inve~tion is now de~cribed in connection with FIGS.
~14. The combined terminal 40 millisecond Q~S and ST s~gment of each signal av~raged X, Y and Z lead is adenti~ied wi~h ~he use of the display cursor lS and standard electrocardiograph criteria 802 and eo4. See also FIG. 9 which depicts withi~ curgor lines ~he combined termi~al QRS and ST ~egment of interes~. In a manner similar to FIG. 5, the above region o.f intere~t i~ multiplied by a four-term Blackman-~arris window 806. Next, each of the signal averaged signals is scaled by identifying the maximum ma~nitude of ~he averaged signal and setting it to unity ao8. Then a 512 point fast~Fourier tra~s~orm is calculated 810 in a manner as described before in connection with FIG.
5 only now ~he FFT is performed on th~ combined terminal 40 millisecond QRS and ST se~ment. The magnitude of the Fouriex transform is then squared 812 to obtain the ener~y spectrum of the signal and the result is plotted 814.
A represe~tative plot of an energy spectrum is shos~n in FIG. 10. To detect smaller peaks which might ba obscuxed by the dominant amplitudes of low frequency components, a secoIld plot ( the broken curve in FIG. 10) was generated by dividing the :, . .

X9a amplitude scale by 500 (see 816). Data was again plotted 818 but data exceeding values in the reduced scale were not plotted. The data was plotted on a high resolution plotter (Versatec, Inc.) and transferred to disc storage.
For each of the spectral plots generated (one for each signal avera~ed X, Y and Z signal for each patient) the data .is first analyzed to locate peaks between 20 and 50 ~z 820. Peaks for the digital data could be defined in ~everal dif~ere~t ways but in ~he preferred e~bodiment, a peak is a~ i~crease in mag~itud~ for at lea~t two adjacent data points followed by a decrease i~ magnitude for at least one point. The fre~uency range of 20 to 50 Hz was chosen after analysis over the bandwidth had demonstrated that freque~cies above 70 ~z did not contribute substantially to the terminal QRS and ST
segments in any group. Because of potential 60 ~z interference, freguencies between 50 to 70 Hz are not analyzed.
The ma~nitude of each peak frequency measured on the reduced scale ~sp~ctral plot of 818) is divided by the maximum magnitude of the entire si~nal measured on the initial scale (spectral plot 814) and the~ multiplied by 105 to form the magnitude ratio 322. The area under the magniied curve (spectral plot al8) between 20 to 50 Hz is divided by the area under the initial curve (spectral plot 814) and then multipLied by 105 to fonm the area ratio 822. The area ratio is computed to determine the relative contribution of compone~ts at freque~cies between 20 and 50 Hz to the entire signal.
For patient-to~patient comparisonæ, the X, Y
and Z value~ for peak reguencies are averaged ~i~

~L~5~3~i~
xsa toge~her and expressed as a single number, the mean of X, Y and Z values . The X, Y and Z values for the pre-multiplied magnitude and area ratios are averaged after log transformation and the mean value for both the magnitude and area ratios expressed as antilogs are then multiplied by a .constant tl x 105) to facilitate graphic display 824. Program No. 4 impleme~ting the abov~
description of FIG. ~ is provided on pages 47_55.
Using the above d~scxibed analysi~ of ~he ener~y spectra of the si~al ave~aged X, Y and Z
signals, three Group~ of patients were studied.
~ roup I comprised 23 patients, each of whom had sustained prior myocar~ial infarction a~d at least one documented episode of sustained ventricular tachycardia or cardiac arrest that was not associated with a new infarct. Group II
consis~ed of 53 patie~ts, each of whom had sustained prior myocardial infarctio~ without sustained ven~ricular taGhycardia (> 30 seconds in duration or associated with imm~diate hemodynamic d*compensation), syncope, or cardia arrest. G.roup III comprised 11 normal subjects.
Ther~ were no significant differences in pertinent clinical features between Group I and II
wi~h respact to age, locus of infarct, or the prese~ce or absence of left ventricular aneurysm.
QRS duration was significantly greater in patients in Group I compared with values in those in Group 30 II (106 ~ Z3 msec vs 90 t 15 msec; p <0.01). Left ventricular ejection fraction was significantly less in patients in group I compared with those in Group II (35 ~ 14% vs 45 t 15%; p ~0.01).
The above analysis of enargy spectra of - 35 si~nal averaged X, Y, and Z electrocardiographic X9a recordings showed significant differences in the area ratios and magnitude ratios of peak frequ~ncies of the combined terminal 40 msec of the QRS and ST segments in patients with pxior infarction associated with subsequent sustained ventricular tachycardia (Group I) compaxed with values in patients without sustained ventricular tachycardia (Groups II and III). There were no signific~nt differences in the area ratios in recordings Xrom patients wi~h myocardial infarction without sustained ventricular tachycardia compared with values in normal subjects.
Representative plots of the type described above of sgyared fast-Fourier transformed data of the terminal QRS and ST segments from a patient wi~h prior myocardial infarction without sustained ventricular tachycardia and from a patient with prior myocardial i~farction with su~tained ventriculax tachyc rdia axe shown in FIGS. llA~llC
and llD-llF, respectively. Each panel depicts magnitude vs. freguency plots of the combined terminal 40 msec of the QRS complex and ST segments o~ signal avera~ed electrocardiographic complexes recorded ~rom bipolar ~, Y, and Z leads along with values from the area ratios, peak frequencies, and magnitude ratios. In each lead, the terminal QRS
and ST segment from the patie~t who had mani~est sustained ventxicular tachycardia contained relatively more high frequency components than the compl~x from the patient without sustained ventricular tachycardia. The dif~erences in freguency content were more marked between 20 to 50 Hz. ~he terminal QRS and S~ segments ~rom the patient who had no manifested sustained ventricular tachycaxdia did contain components with frequ~ncies 5~7 xsa -18~

above 20 Hz. However, ~he overall contribution of these components to ~le entire signal was modest as refl~cted by the low area and magnitude xatio values.
FIGS. 12A-12C and 12D-12F illustrates another example of energy spectra of the terminal QRS and ST seyments from patients without a~d with sustained ventricular tachycardia, respectively.
The terminal QRS and ST segments from both patients contained components with frequencies between 20 to 50 Hz but differed markedly with respect to ~he relative contributions of these components to the entire signal. In the X and Z leads, the terminal QRS and ST segments from ~he patient who had manifest sustained ventric~lar tachycardia contained a 10 to 100-fold greater proportion of components in the 20 to 50 ~z range compared with corresponding el~ctrocardiographic segments rom ~he patient without sustained ventricular 20 tachycardia .
Comparisons between valuss for the area ratios for ~he thre~ patient groups are shown in ~I&. 13~ There were significant group differences in area ratio values betwee~ patients with prior myocardial in~arction who had manifest sustained ventricular tachycardia compared with patie~ts with prior myocardial infarction without sustained ventricular tachycardia and with normal subjects.
~ 10-fold difference was evident for patients with compared to those without sustained ventricular tachycaxdia. There was no overlap between values among patients who had manifest sustained ventricular tachycardia and values from normal subjects. However, twelve patients with prior myocardial infarction without clinically documented X9a sustained ventricular tachycardia had area ratios overlapping values from patients who had manifest sustained ventricular tachycardia. It is not ye t clear whether these 12 patients are at increased risk for development of sustained v~ntricular tachycardia.
FI~. 14 illustrates the me~n peak frequencies detected between 20 to 50 ~z and the corresponding magnitude ratios for th~ ~hxe~ patient groups.
There were no significant differences in ~he peak frequencies among patients i~ the three groups (28 4 Hz and 41 ~ 4 ~z, in Group I; 30 ~ 3 ~z and 43 ~ 3 ~z, Group II; ~nd 28 ~ 5 ~z and 43 ~ 3 ~2, Group III). ~owever, ~he relative contribution of the magnitudes of ~hes~ peak freguencies to the ovarall magnitude of th~ spectxal plot of their terminal QRS and ST seg~ents differed signi~icantly. Thus, the terminal QRS and ST
seg~ntC in patients with prior ~yocardial infarction who had manifest sustained ventricular tachycaxdia contained a 10 to 100-fold greater contribution from components in the ~0 to 50 Hz range compared with ~orrespo~ding electrocardiogxaphic segme~ts in patients without sustained ventricular tachycardia. No frequencies above 50 ~z contributed substantially to those electrocardiographic segments in any yroup.

~, .

d ~
X9a PROGRAM NO. 1 ECG Signal averaging package:

ECGAVG takes an incoming digi-tized Electrocardiograph signal and averages normal QRS
complexes. ECGAV& has two sections, a learning period and an averaging period.
During the learning period, the clinician may choose one of the three ECG channels X, Y, or Z
(channels 0, 1, or 2 respectively) for analysis.
The waveform of the desired channel is then displayed at the terminal for the clinician to mark the peaks of two adjacent normal QRS's. This serves to set the peak detector threshold, the R-R
interval threshold, and to provide a template waveform for the signal-averaging routines.
During the averaging period, the clinician enters the number of beats to average. Once -this is done, the peak detector FINDPK looks for signal values above (or below) the QRS threshold for positive-going (or negative) QRS's that are within the given R-R interval of each other. If a peak is found, it is then correlated against the template waveform using ECGCOR. If the correlation coeffi.cient of the detected peak and the template is above correlation threshold (right now at about 98% confidence~ then the beat is flagged as a normal QRS. If two such beats are found within the R-R inverval of each other, the second of these two ~2~
X~a is then averaged into an array. When the desired number of beats have been averaged, -the program then dumps the three channel names (0, 1, and 2) in the order that -they have been averaged (integer values) and then waveform arrays in that same order, all in unformat-ted form.

Files for the Signal Averager are:

ECGAVG Signal averager ECGPK - contains subroutine FINDPK, the peak detector.

ECGCOR - contains subroutines I22COR
and I24COR for calculating the correlation coefficient over a 32 sample point window between an I*2 an an I*2 array or between an 1*2 and an I*4 array.

ECGFOR - contains the Selenar Graphic driver routines.
ECGMAC - contains the macro routines caled by ECGAVG and ECGFOR.

XXLIB - contains the ADAC A/D
subroutines.

Once the ECG signal has been averaged, the waveforms may then be displayed on the terminal for windowing using ECGFFT. The unformatted data output files of ECGAVG are read in-to ECGFFT and X9a displayed a-t the terminal. The clinician is then asked how may plot sets ~windowed sections) are requlred for the displayed waveform. The clinician then sets the crosshair cursor at the beginning and ending points of each beat section. The data between the cursors is then multiplied by the window funtion and is written out -to disk in el3.6 format, preceeded by the number of plot sets (16 format~ and the numher of data points (also 16 format). The next set of data points is then windowecl and written to disk also preceeded by its number of data points, until the number of plot sets has been reached. The windowed data can then be sent to the Perkin-Elmer for FFTing and plottiny.

Files related to ECGFFT are:

ECGFFT - Displays channels X,Y,Z and composite waveforms always in that order and allows sections of the waveforms to be windowed and output to disk in formatted form.
ECGFOR - contains subroutine for Selenar Graphis board.

FFTMAC - contains mere routines called by ECGFFT and ECGFOR.

WINDOW - contains WNDW, the windowing function.

Also FF~NOW does the same thing as ECGFFT

X9a but does not window the c Subroutine findpk steps -through NPTS sample C points at DELD intervals starting a-t -the sample C polnt ISMP and looks for a change in first C derivative from pos. to neg. or neg. to pos. If C no such point is found (ie: a sec-tion of rising or C falling or zero slopes) then a 0 is returned.
C

subroutine findpk (ismp,irtsmp,npts,deld,istat) integer sam~9216),smp,ismp,irtsmp,npts,deld,isvmax,i,j,k, tb~sz,istat,loedge,hledge,ampth,iendpt byte nmbr(6) data tbfsz./9216/,nmbr/6*0/
common/smpl/sam,loedge,hledge,ampth smp = lsmp istat = 0 irismp = 0 iendpt = ismp + npts*deld if (iendpt .gt~ 9216) iendpt = iendpt - 9216 c c... Check curwrt buffer edges to see that data being analyzed c... does not get written over. If data is overrun by sampler c... a status = 1 is sent, otherwise status = 05 c if (deld .lt. 0) go to 56 !if neg.directed peak search go to 57 c c lf (ismp .gt. hledge) !if high edge of written go to 5~ buffer is if (ismp .lt. loedge) !greater than ismp or go to 52 low edge of istat = 1 !written buffer is .ge.
than ismp x9~
- 2~ -call cvidec(ismp,nmbr) !-then data is overrun by sampler call print('data overrun a-t ismp:') call print(nmbr) return c 52 if (iendpt .lt. loedge) !if trailing edge of go to 60 curwrt is if (iendpt .gt. hiedge) !.le.ismp,wait for go to 60 next buffer go to 52 c 56 if (iendpt .gt. hiedge) go to 58 if (iendpt .It. loedge) go to 58 istat = 1 call cvidec (ismp,nmbrx) 'then data is overrun by sampler call print ('iendpt at overrun is:') call print (nmbrx) return 25 c 58 if (ismp .lt. loedge) go to 60 if (ismp .gt. hiedge) go to 60 go to 58 c c... Now at peak detecting logic:
c 60 if (ampth .lt. 0) !go to negdirected qrs go to 100 logic if neg qrs ~5~
- 25 _ X9a do 62 i = l.npts !Check for smp above amp theshold J = 1 if (sam~smp~ .ge. ampth) go to 63 smp = smp ~ deld if (smp . gt. tbfsz~ smp = smp-tbfsz !check for positive wraparound if ~ smp . le. 0) smp = smp + tbfsz ! check or negative wraparound 62 continue go to 280 63 lsvmax = sam (smp) do 75 i = j,npts !Check if positon is k = i on downslope smp = smp + deld if (smp .gt. ibfsz) !check for positive smp = smp-tbfsz wraparound i f ( smp .le. 0) 'check for negative smp = smp + tbfsz wraparound if (sam(smp) 'look for downslope .gt. lsvmax) go to 65 isvmax = sam(smp) 75 continue go to 280 65 do 76 i = k,npts !Now that we're on an if (sam~smp~.lt. isvmax) upslope ne~t downslope go to 275 !look for ~/- trans.
isvmax = sam(smp) smp = smp + deld if (smp .gt.tbfsz) !check for postive smp = smp-tbsz wraparound if (smp .le. 0) !check for negative smp = smp = tbfsz wraparound 76 continue go to 280

3~7 X9a 100 do 162 i = l,npts !Check for sample value below -threshold J = i if (sam(smp) .lc. ampth) go to 163 smp = smp + deld if (smp .gt. tbfsz) !check for positive smp = smp - tbfsz wraparound if (smp .le. 0) 'check for negative smp = smp + tbfsz wraparound 162 continue go to 280 163 isvmax = sam(smp) do 175 i = npts k = i smp = smp + deld if (smp .gt. tbfsz) !positive wraparound smp = smp - tbfsz if (smp .le 0~
smp = smp + tbfsz !negative wraparound if ~sam(smp3 .lt. !look for up slope isvmax) go to 165 isvmax = sam(smp) 175 con-tinue go to 280 165 do 176 i = k,np-ts if (sam(smp) .gt. !look for -/+
isvmax) go to 275 transition isvmax = sam(smp) smp = smp ~ deld if (smp .gt. tbfsz) !positive wraparound smp = smp - tbfsz if (smp .le. 0) !negative wraparound smp = smp + tbfsz 176 continue go to 280 275 irtsmp = smp ~ deld X9a 280 return end C Subroutine 124COR takes t~o arrays (one i*2 C and the second, i*4)and computes their C correlation coeeficeint over the NPTS sample C points at intervals of DELD using SYSLIB i*4 math C routines. The return argumen-t CORR is a r*4 C variable.
C The first argument passed is the index to the C sample array SAM. The second arguement is the C array to be correlated against SAM.
c subroutine 124cor(index,array,npis,deld,corr,istat) integer sam~9216),i,j,index,npis,deld,ibeats, istat,loedge,hiedge,iendpt,nmpth,init(l0) integerk4 array(npts),sumO,suml,sum2,tmpl,tmp2,tmp3, ibeat,templ,temp2 real corr,num,denl,den2 byte nmbr (6) equivalence ~sumO,init~l)),suml,init(3)), (sum2,init~5)),(templ,init(7)), (temp2,init(9)) common/smpl/snm,loedge,hiedge,ampth common/chn/ichn,ibeats istat = 0 iendpt = index -~ npts*deld if (lendpt .gt. 9216) iendpt = iendpt-9216 c c.......... Check for data overrun by sampler or sampler overrun by program logic:
c if (index .gt. hiedge) go to 50 if (index .lt. loedge) go to 50 call cvtdec(index,nmbr) CALL PRINT('data overrun at index:') call print(nmbr) xga lS tat=l !error data o~errun return c 50 if (inendp-t .lt. loedge) gO to 60 if ~'iendpt .gt. hiedge) ~o -to 60 go to 50 c 60 call jicvt (ibeats,ibeat) c c... Initialize variables do 62 i - -1,10 init~i) = O
62 continue c do 100 i = l,npts if (index, .gt. 9216) !Check for pos.wraprnd index = index - 9216 call jicvt(sam(index),templ) j = jdiv(array(i),i~eat,temp2 j = jmul(templ,temp2,tmpl) j = jadd(tmpl,sumO,sumO) jmul(templ,templ,tmp2) j = jadd(temp2,suml,suml) j = jmul(temp2,temp2,tmp3) j = jadd(tmp3,sum2,sum2) Index = index -~ deld 100 continue num = ajflt(sumO) den:L = ajflt(suml) den2 = ajflt(sum2) corr = num/s~rt(denl*den2) return end qr7 xga C Subroutine 122COR performs the same function C as 124COR but on two C sections of sam indexed by INDEXl 8 INDEX2 using 1*4 mnth.
C
subroutine 122cor(indexl,index2,npts, deld,corr,istat) integer sam(9216),i,j,indexl,index2, npts,deld,ibeats, istat,loedge,hiedge,ampth,init(10) integer*4 sumO,suml,sum2,templ, temp2,tmpl,tmp2,tmp3 real corr,num,denl,den2 eguivalence (init~l),sumO),(init(3),suml), (init(5),sum2), (init(7),templ),(init(g),temp2) common/smpl/sam,loedge,hledge,ampth common/chn/ichn,ibeats istat = O
c c... ..Initialize variables c do 50 i - 1,10 init(i) = O
50 continue do 100 i = l,npts if (indexl .gt. 9216) !Check for pos.wraprnd indexl = indexl-9216 if(index 2 .~t. 9216) index 2 = index 2-9216 call jicvt(sam~indexl),-templ) call jicvt(sam(index2),temp2) j = jmul(teml,temp2,tmpl) j = jadd(tmpl,sumO,sumO) j = jmul(templ,templ,tmp2) - j = jadd(tmp2,suml,suml~
j = jmul(temp2,temp2,tmp3) ~s~
xga j = jadd(tmp3,sum2,sum2) indexl = indexl -~ deld index2 = index2 + deld 100 continue num = ajflt(sumO) denl = ajflt(suml) den2 = ajflt(sum2) corr = num/sqr-t(denl*den2) return end ,. .

X9a PROGRAM NO. 2 PDP-ll FORTRAN-77 V4.0-3 15:52:38 5-Dec-83 ECGFFT.FOR:l /F77/TR:BLOCKS/WR

c and FFTing, NPLOT is the number of plot sets to be done.
0006 integer nsamp(10) 0007 integer nixsmp(10 ,niysmp(10) 10 0008 integer nplot c c These vars are used to calculate the composite beat.
0009 real cwave(1024) ,templ,temp2,temp3 c c FFT_hold the windowed data to be output to disk.
0010 real fftO(514),fftl(514),fft2(514), fftc(514) c 0011 byte query ~80) 0012 byte cr 0013 byte nmbrx(7) 0014 byte input(6) 0015 byte infile(10) 0016 byte outfil(11) X9a 0017 byte output(6) c 0018 equivalence (ichn,iobuff(1)), (jchn,iobuff(2)),(kchn,iobuff(3)), .(iwave,iobufE(4)),(jwave,iobuff(1028)), kwave,iobuff(2052)) 0019 equivalence (input,infile),(output,outfil) c 0020 common /zoom/gainx,gainy c 0021 data idnam/3RDY /
0022 data dblk/3RDY1,0,0,0/
0023 data nmbrx/7*0/
0024 data nplot/0/
15 0025 data nsamp/10*0/
0026 data infile/6*40,'D','A','T',0/
0027 data outfil/6*40,'.','F','F','T',0/
c c..... Initialize:
c -0028 gainx = 0 0029 gainy = -4 c 0030 call prin-t('Welcome to ECGFFT version 1.0') c c.. ...set up I/0 channel:
c 0031 iochan = igetc() 0032 if (lochan .lt. 0) stop 'No channel available' 0033 25 call print('Enter name of input file: (6chars~') 0034 read(5,27) input 0035 27 format(6a) 0036 72 call print('Enter the number of ~Z~;~9~7 X9a plots to run: (0-10)') 0037 read(5,74) nplot 0038 74 format(16) 0039 if (nplot .eq. 0) go to 69 5 0040 call print ('Enter name of output file: (6chars)') 0041 read(5,27) output c c... Establish input channel:
c PDP-11 FORTRAN-77 V4.0-3 15:52:38 5-Dec-83 ECGFFT.FOR;1 /F77/TR:BLOCKS/WR

c and FFTing. NPLOT is the number of plot sets to be done.
0006 integer nsamp(10) 0007 integer nixsmp(l0),nlysmp(10) 0008 integer nplot c c These vars are used to calculate the composite beat.
0009 real cwave(1024),templ,temp2,temp3 c c FFT_hold the windowed data to be output to disk.
0010 real fft0(514),fftl(514), ft2(514),fftc(514) c 0011 byte query(80) 0012 byte cr 0013 byte nmbrx(7) 0014 byte input~6) 0015 byte infile(l0) 35 0016 byte outfil(ll) X9a 0017 byte output(6) c 0018 equivalence (ichn,iobuff(l)), (jchn,iobuff(2)),(kchn,iobufE(3)), .(iwave,iobuff(4),(jwave,iobuff(1028)), (kwave,iobuff(2052)) equivalence (input,infile), (output,outfil) c 10 0020 common/zoom/gainx,gainy c 0021 data idnam/3RDY /
0022 data dblk/3RDY1,0,O,0/
0023 data nmbrx/7*0/
0024 data nplot/0/
0025 data nsamp/10*0/~
0026 data infile/6*40,'D','A','T',0/
0027 data outfil/6*40,'.i,'F','F','T',0/
c c..... initialize:
c 0023 gainx = 0 002~ yainy =
c 25 0038 call print('Welcome to ECGFFT
version 1.0') c c... set up I/0 channel:
c 0031 iochan = igetc() 0032if(iochan .lt. 0) stop 'No channel available' 003325 call print('Enter name of input file:
(6chars~') 0034read(5,27) inpu-t "~

.

~5~
xsa ~ 35 -0035 27 format(6a) 0036 72 call print('Enter the number of plots to run: (0-10)') 0037 read(5.74)nplot 0038 74 format(16) 0039 if (nplot .eg. 0) go to 69 0040 call print('Enter name of output file:
(6chars)') 0041 read(5,27) output c c... Establish input channel:
c FFT-ll FORTRAN-77 V4.0-3 15:52:30 5-Dec-83 ECGFT.FOR,1 0042 69 i = irnd50(9,infile,db1k(2)) 0043 if (i .ne. 9) goto 25 0044 if (ifetch(idnam) .ne. 0) stop 'fatal ~rror fetching handler' 0045 if (lookup(iochan,dblk) .lt. 0~ stop 'lookup error' c c..... Read in ichn,jchn,kchn and symbol data:
c 0046 icode - irendw(3075,iobuff,0,iochan) 0047 call cvtdec(icode,nmbrx) 0048 call print(nmbrx) 0049 if (icode .lt. 0) stop 'Error inputting data' 30 0050 if (iclose(iochan) .lt. 0) stop 'Error closing file' 0051 if (ifreec(iochan) .lt. 0) stop 'Channel not allocated' c c..... Compute composite channel:

~2S~ 7 xsa 0052 do 75 i = 1,1024 0053 -templ = float(iwave(i)) 0054 temp2 = float(jwave(i)) 5 0055 temp3 = float(ikwave(i)) 0056 cwave(i) = sqrt(templ*templ -~ temp2*temp2 -~ temp3*temp3) 0057 fwave(i) = int(cwave(1)) 0058 75 continue 0059 go to (81,82,83),ichn~1 c c... If ichn = 0 then reset channels as follows c 006081 ichn = 2 0061jchn = 0 0062kchn = 1 0063go to 84 c............... ........ If ichn = 1 006482 ichn = 0 0065jchn = l 0066kchn = 2 0067go to 84 c............... ... If ichn = 2 then 006883 ichn = 1 0069jchn = 2 0070kchn = 0 007184 continue c c..... Enter Selenar mode and plot out data c 0072 call clrtrm() 0073 call setsel() 0074 call clrscn() 35 0075 call stgrph() 3l~5~'~
xga 0076 call stzoom(gainx,ga1ny) 0077 call swvplt(iwave,1024,(4-ichn)*50,1) 0078 call swvplt(jwave,1024,(4-jchn)*50,1) 0079 call swvplt(kwave,1024,(4-kchn)*50,1) 5 0080 call swvplt(fwave,1024,25,1) 0081 call sreset() 0082 call setirm() 0083 if (nplot .eq. 0) go to 250 10 PDP-11 FORTRAN-77 V4.0-3 15:52:38 5-Dec-83 ECGFFT.FOR;1 /F77/TR:BLOCKS/WR

0084 do 100 k = l,nplot 0085 call print('Position cursor at 1st point:
.(press<cr> to cont.)') 0086 call setsel() 0087 call oncrsr() 0088 80 cr = ittinr() 20 0089 if (cr .ne. 13) go to 80 0090 cr = ittinr() 0091 call getcrs(ixl,iy) 0092 call offers() 0093 call stgrph() 25 0094 call sline~ixl,l000,ixl,50) 0095 call sreset() 0096 call settrm(~
0097 call clrtrm() 0098 call print('Position cursor at 2nd point:

.(press <cr> to cont.)') 0099 call setsel() 0100 call oncrsr() 0101 90 cr = ittinr() 35 0102 if (cr .ne. 13) go to 90 ^~L'Z Sfi~
xsa 0103 cr = i-ttinr() 0104 call getcrs~ix2,iy) 0105 call offers() 0106 call stgrph~) 0107 call sline(ix2,1000,1x2,50) 0108 call sreset() 0109 call settrm() 0110 call clrtrm() 0111 NIXSMP(K) = IX1 0112 NIYSMP(K) = IX2 0113 nsamp(k) = ix2-ixl 0114 100 continue c c... Open channel and output number o~ plots:
c 0115 OPEN(UNIT=l,NAME=outfil,TYPE='NEW', FORM='FORMATTED') 0116 write(1,105)nplot 0117 105 format(i6) c c... Load windowed data into fft arrays:
c 0118 do 160 1 = l,nplot 0119 do 110 1 = nixsmp(l),niysmp(1) 0120 iitmp = i-nixsmp(1)+1 0121 fftO(iitmp~ = float(iobuff (i+4-~ichnlO24)) . *wndw(iitmp,nsamp(1)) 30 0122 fftl(iitmp) = float(iobuff (i+4+jchn*1024)) . *wndw(iitmp,nsamp(1)) 0123 fft2(iitmp) = float(iobuff X9a (i+4-~kchn*1024)) . *wndw(iitmp,nsamp(l~) 0125 llO contlnue c c..... Zero out (decimate) the remainder of the FFT arrays:
c PDP-ll FORTRAN-77 V4.-3 15:52:38 5-Dec-83 CCGFFT.FOR;l /FFT/TR:BLOCKS/WR

c With the actual computation and plotting of the FFT's being done on anokher c computer (the inter-data or perkin-elmer) and since the data transfer from c the LSI-11 to the Inter-Data is so slow, only the beginning of the array c with the windowed data is sent out so the decimation here is not really c necessary.
126 do 120 i = nsamp(l)-1,514 127 fftO(i) = 0.
128 fftl(i) = 0.
129 ff~2(i) = 0.
130 f~tc(i) = 0.
131 120 continue c c....... Output the windowed data to disk, the number of points and the arrays:
c 132 write(l,125)nsamp(1) 133 125 format(16) c c... Output arrays:

,, X9a 134 130 format(cl3.6) 135 do 140 i = l,nsamp(1) 136 write(1,130)ff-tO(i) 137 140 continue 138 do 145 i = l,nsamp(l) 139 write(1,130)fftl(i) 140 145 continue 141 do 150 i = 1,usamp(1) 142 write(l,130)fft2(i) 143 150 continue 144 do 155 i = l,nsamp(1) 145 write(1,130)fftc(i) 146 155 continue 147 160 continue c c... Close Channel:

c 148 CLQSE(UNIT=1) ~0 c 149 250 call print('Bye,bye') 150 stop 151 end , X9a P ROGRAM NO . 3 PLOTBS

C. READ FROM MISC AND COMPUTE FFT.
AND PLOT DATA FOR DA
C. MODIFY TO REMOVE ZERO OFF-SET
C. JOANNE MAPLHAM JAN 1~83 C.
IMPLICIT.INTEGER*2 (l-N) DIMENSION.X(512~.Y(1024).ZX~5).ZY(5) COMPLEX C(512) EQUIVALENCE (C(l).Y(I)) C.
PAT=FLOAT(X'5555') C.
~EAD (3) NSET
ISET~l S CONTINUE
2G IRLOT=l ~0=-0 . 1 Y0=7.2 X~IN=0.
XMA~=520.
YMIN=-40.
YMAX=1~0.
XL=6.5 YL=2.
XE=XU+XL
NOP=257 XINC=40.
XT=.5 YINC=40.
~ YT=0.5 CONTINUE
YB=Y0~YL

J=l DO 20 I=1.8 K=J+63 READ (3) (X(N),M=J,K) J=J+64 1001 FORMAT (lOF12.2) D0 40 I=1,512 C(l)=CMPLX(X(I)Ø0) J=J+64 1001 FORMAT (10F12.2) DO 40 I=1,512 C(l)=CMPLX(X~1)Ø0) CONTINUE
M=9 CALL FFT(C,M) C. COMPUTE Y = 20*LOGl0(F),M(K)=K*1000/512 C. CALCULATE AREA FROM 0 TO DECREASE OF 60 DB
C.
1002 FORMAT(10E12.4) J=l Yr~AX=O .
DO 100 I=l,NOP
X(I)+(I-1)*1000./512.

TEMP=SORT(Y(J)*Y(J)+Y(J+l)*Y(J+l)) IF (TEMP.E0.0) GO TO 80 Y(I)=20.*ALOG10(TEMP) Y(I)=YMIN

IF (YMAX.GT.Y(I)) GO TO 100 YMAX=Y(I) 100 J=J+2 1003 FORMAT (5(F10.1,F10.2)) C.

i7 xsa C. AREA
C.
YJ=Y(1)-60.
A=Y(l) D0 110 I=2,500 IF(Y(I).LT.YJ(GO TO 120 110 A=A~Y(I) 120 A=A*1000./512 C. FIND VALUE AT 40HZ
R=(40.-X(21))/(X(22)-X(21)) B=Y(21)~R*(Y(22)-Y~21)) B=Y ( 1 ~ -B
CALL MODE~4,0.07,0.05,0.) ZX(l)=X0 ZX(2~=XE
ZY(l)=Y0 ZY(2)=Y0 CALL DRAW (ZX,ZY,2,1) XX=XO
FN=0.
ZY(l)=Y0-0.04 ZY~2~=Y0 YA=Y0-.15 DO 220 I=1,20 IF (XX.GT.XE) GO TO 240 ZX ( 1 ) =~X
Z~(2)=XX
CALL DRAW (ZX,ZY.2.1) XA=XX-0.11 IF (FN.LT.10.)XA=XX-.03 CALL NOI'E (XA,YA,FN,1000) FN=FN~XINC
XX=XX-~XT

~S~
X9a C.
C. DRAW DASHED LINES FOR X-GRID
C.
XX=XO
CALL MODE (10,PAT,9999.,1.) ZY(l)=Y0 ZY(2)=YE
DO 250 I=1.20 XX-XX~XT
IF (XX~GT~XE~ GO TO 260 ZX(l)=~X
2X(2)=XX
CALL DRAW (ZX,ZY,2,1) C.
C. DRAW DASHED GRIP FOR Y AXIS
C.
YY=YO
Z~(l)=X0 ZX(2)=XE
DO 200 I=1,20 YY=YY~YT
IF (YY.GT.YE) GO TO 300 ZY( 1 )=YY
ZY(2)=YY
280 CALL DRAW (ZX,ZY,2.1) C
C RESET LINE TO SOLID

C

CALL MODE (10.-1.,9999.,1.) C

xsa C DRAW Y AXIS AND LABEL
C

ZX ( 1 ) =XO
ZX(2)=X0 ZY(l)=Y0 ZYI2)=YE
CALL DRAW (ZX,ZY,2,1) Z~l)=Y0-0.05 ZX~2)=X0 1 0 YY=YO
YN=YMIN
DO 320 I=1,20 IF (YY,GT,YE) GO TO 350 ZY(l)=YY
ZY(2)=YY
CALL DRAW (ZX,ZY,2,1) YJ=ZY(1)-0.05 XJ=X0 .3 IF (ABS(YN) .GE.100.) XJ=X0-.36 IF (YN.E0Ø) XJ=X0-.15 CALL NOTE (XJ,YJ,YN.1000) YN=YN+YINC
320 YY=YY~YT
350 CALL MODE (4,0.1,0.067,0.) YJ=Y0~
CALL NOTE ~5.,YJ,'AREA=',5) CALL NOTE (5.6,YJ,A,1000) YJ=YJ-.2 CALL NOTE (5.0,YJ,'F(40HZ)=',8) CALL NOTE (6.0,YJ,B,1000) XS=XT/XINC
YS=YT/YINC
DO 400 I=l,NOP
X(I)=X(I)*XS+X0 Y(I)=(Y(I)-YMIN~*YS+Y0 ,. . .

X9a - ~6 -CALL DRAW (X,Y,NOP,l) XJ=X0-.38 YJ=Y0+.5 CALL MODE (4,0.1,0.067.90) CALL NOTE (XJ,YJ,'FFT-MAG (DB)',12) CALL MODE (4,0.1,0.067.0) XJ=X0+XL/2,-.5 YJ=Y0-.3 CALL NOTE (XJ,YJ,'FREQ (HZ)',9) IF (IPLOT.E0.4) GO TO 500 IPLOT=IPLOT+l Y0=Y0-2.4 500 CALL DRAW(0.,0.,1,9000) IF (ISET.EQ.NSET) GO TO 600 ISET=ISET~l CALL DRAW (0.,0.,0,9999) STOP.
END.

~25~7 xga PROGRAM NO. 4 C READ FROM DISC AND CONPUTE FFT. AND
PLOT DATA FOR DA

C SCALE INPUT DATA TO 1.0 C ~ODIFY FOR 2 PLOTS ON SINGLE SCALE
C
IMPLICIT INTEGER*2 (I-N) DIMENSION X(512),Y(1024),ZX(5),ZY(5),Z(512) DIMENSION IPEAK(10) COMPLEX C(512) EQUIVALENCE (C(l),Y(1)3 EQUIVALENCE tY(513),Z(l)) C

PAT=FLOAT(X'5555') C

READ (3) NSET
ISET=l S CONTINUE
IPLOT=l XO=O .
~0=7.2 XMIN=0.
XMAX=120.
YMIN=0.
XL=6.
YL=2.
XE=XO~XL
NOP=61 XINC=10.
XT=.5 YT=0.5 CONTINUE
YE=YO+YI.
J=l X9a DO 20 I=1.8 K=J+63 RE~D (3) (X(M),M=J,K) J=J+64 1001 FORMAT (10F12.2) S=O .
YM=O.
DO 30 I=1,512 IF (AB5(X(I)).LE.YM) GO TO 30 YM=ABS(X(I)) S=S~X(I) 0 35 I=1,512 K=513-I
IF (X(K).NE.0) GO TO 37 CONTINUE

S-S/K
1009 FORMAT (' MEAN =',F10.4) DO 40 I=1,512 X(I)=X(I)/YM
C(I)=CMPLX(XtI),O.O) CONTINUE
WRITE (6,1031) (X(I),I=l,K) 1031 FORMAT('lATA AFTER SCALING'/(10F 10.4)) WRITE (6,:1009) S
M=9 CALL FFT(C,M) C. COMPUTE POWER SPECTRUM
C.
1002 FORMAT (10E12.4) J=l FMAX=0.
DO 100 I=1,257 X(I) (I-1)*1000./512.
Y(I~=(Y(J)*Y(J)+Y(J+l)*Y(J+l))-~ .

xsa IF (FMAX.GT.Y(I)) GO TO 100 FMAX=Y(I) JK=I
100 J=J~2 1003 FORMAT (5F10.1.F10.2)) WRITE (6,1007) E'MAX,Y(I).I=1.150) 1007 FORMAT ('OAUTO-CORRELATION ',F10.3/(10E12.4)) WRITE (6,1040) IPLOT,ISET
1040 FORMAT ('1 PLOT ',13,' OF SET ',13/) C. AREA FROM 0-20 HZ & 20-50 HZ
AR=0.
DO 700 I=l,ll 700 AR=AR~Y(I) AR=AR*X(2) ARB=0.
DO 720 I=12,27 720 ARB=ARB=Y~I) ARB=ARB*X(2) R=ARB/AR
WRITE (6,1041) AR,AR8,R
1041 FORMAT ('0 AREA FROM 0-20 HZ = ',E14.5./
1'0 AREA FROM 20-50 HZ = ',E14.5,5X, 'RATIO = ',El4.5) WRITE (6,1042) Y(JK),X(JK) 1042 FORMAT ('0 PEAK = ',F10.2. ' AT ',F8,2,' HZ'/) R=AR/FMAX
WRITE(6,1043)R
1043 FORMAT ('0 AREA 0~20 HZ/PEAK = ',F14.4) C. FIND PEAKS FROM 20-50 HZ ( POINTS NO 8-27) NPK=0 IP=0 ID=l SAVE =Y(7) C. DO 880 I=8,27 SUM=Y(I) 3l25~ 7 xsa DIF =SUM-SAVE
IF ~DIF)810,870,840 ID-ID ~l IF (ID.EQ.l) GO TO 820 IP=0 020 IF (IP.LT.2) GO TO 870 IF (NPK.EQ.10) GO TO 890 NPK=NPK+l IPEAK(NPK~=I-l IP=0 840 IP=IP~l 870 SAVE =SUM

890 WRITE (6,1047) 1047 FORMAT ('MORE THAN 10 PEAKS
FOUND--ONLY lST 10 ACCEPTED'~
895 WRITE '6,1046)NPK
1046 FORMAT ('0',I4,' PEAKS FOUND FROM
FREQ 20- 50 HZ'//6X.
l',FREQ',llX,'PEAK',llX,'RATIO'//) IF NPK,EQ,0) GO TO 898 DO 897 I=l,MPK
J=IPEAK(I) R=Y(J),FMAX
897 WRITE (6.1048)X(J),Y(J),R
1048 FORMAT (F10.2,E15.4,E15.1) DV=0.1 120 K=FMAX/DV
IF (K.EQ.0) GO TO 140 9~
- 51 - X9a DV=DV*10.

K=(FMAX*10/DV)-~l YMAX=(K*DV)/10.
YINC=YMAX/4 .
C.
C DRAW X AXIS
CALL MODE (4,0.07,0.05,0.) ZX(l)=X0 ZX(2)=XE
ZY(l)=Y0 ZY(2)=Y0 CALL DRAW ~ZX,ZY,2,1) XX=XO
FN=G.
ZY(l)=Y0-0.04 ZY(2)=Y0 Y~=Y0-.15 DO 220 I=1,20 IF (XX.GT.XE) GO TO 240 ZX ( 1 ~ -XX
ZX(2~=XX
CALL DRAW (ZX,ZY,2,1) XA=XX-0.11 IF (FN.LT.10.) XA=XX-.03 CALL NOTE (XA,YA,FN,1000) FN=FN~XINC
XX=XX+XT

C.
C. DRAW DASHED LINES FOR X-GRID
C.
XX=X0 ~S~4~
xsa CALL MODE (10,PAT,9999.,1.) ZY ( 1 ) -~0 ZY(2)=YE
DO 250 I=1,20 S IF (XX.GT.XE) GO TO 260 ZX(l)=XX
ZX(2)=~X
CALL DRAW (ZX,ZY,2,1) XX=XX+XT

C.
C. DRAW DASHED GRID FOR T AXIS
C.
YY=Y0 ZX(l)=X0 ZX(2)=XE
DO 200 I=1,20 YY=YY+YT
IF (YY.GT.YE) GO TO 300 ZY(l)=YY
ZY(2)=YY
280 CALL DRAW (ZX,ZY,2,1) C.
C RESET LINE TO SOLID
C.

CALL MODE(10,-1.,9999.,1.) C.
C. DRAW Y AXIS AND LABEL
C.
LN=1000 IF (FMAX.LT.1.0) LN-1002 IF (FMAX.LT.10.) LN=1001 ZX(l)-X0-0.05 X9a ZX(2)=X0 YY=YO
YN=YMIN
DO 320 I=l,20 IF (YY.GT.YE) GO TO 350 ZY( 1 )=YY
ZY(2)=YY
CALL DRAW (ZX,ZY,2,1) YJ=ZY(1)-0.05 XJ=~0-.3 IF (ABS(YN) .GE.100.) XJ=X0-.4 IF (YN.EQØ) XJ=X0-.15 CALL NOTE (XJ,YJ,YN,LN) YN=YN~YINC
320 YY=YY~YT

ZX(l)=XE
ZX(2)=XE
ZY(l)=Y0 Z~(2)=YE
LN=1001 CALL DRAW (ZX,ZY,2,1) YINCB=YINC/500.
YN=YMIN
ZX(2'=X~0.05) YY=Yq DO 360 I =1,20 IF (YY.GT.YE) GO TO 370 ZY( 1 )=YY
ZY(2)=YY

CALL DRAW (ZX,ZY,2,1) YJ=Z~(1)-0~05 XJ=XE~0.1 CALL NOTE (XJ,YJ,YN,LN~
YN=YN+YINCB

*7 xsa YY=Yy-tYT

CALL MODE (4,0.1,0.067,0.) YJ=YO~1.8 XS=XT/XINC
YS =YT/Y I NC
1:)0 400 I=l,NOP
X( I )=X( I ) *XS+XO

4 0 0CONT I MUE
CALL MODE (10, PAT , 9 9 9 9 . , 1 . ) CALL DRAW ( X, Z, NOP, 1 ) CALL MODE (10,-1.,9999. ,1. ) YS=YT/Y I NCB
DO 420 I=l,NOP
Y(I)=(Y(I)-YMIN)*YS+Y0 IF (Y(I).GT.YE) Y(I)=YE

CALL DRAW ( X, Y, NOP, 1 ) XJ=XO- . 53 ~J=Y0+ . 5 CALL MODE ( 4,0.1,0.067,90.) CALL NOTE (XJ,YJ, 'AUTO-COEIR ' ,12) CALL MODE ( 4,0.1,0.067,0.) XJ=XO+XL/2 . - . 5 YJ=YO- . 3 CALL NOTE ( XJ, YJ, ' FE`~EQ ( HZ ) ' . 9 ) IF (IPLOT.EQ.4) GO TO 500 IPLOT=IPLOT-tl YO=YO-2.4 500CALL DRAW (0., 0. ,1, 9000) I F ( I SET . EQ . NSET ) GO TO 6 0 0 I SET- I SET-tl ~L2~
xsa CALL DRAW ( 0 ., 0 . . O, 9999 ) STOP
END

Claims (28)

X9a The embodiments of the invention in which an exclusive property or privilege is claimed are defined as follows:

1. A method for analyzing electrocardiogram (ECG) signals to determine the presence or absense of a predetermined frequency content in a preselected portion of said ECG signals comprising the steps of:
converting analog ECG signals to digital ECG
signals;
performing a fast Fourier transform (FFT) on said preselected portion of said digital ECG
signals; and determining a figure of merit (FFT) associated with the frequency content of said at least a preselected portion of said ECG signals from said output of said Fourier transform step, whereby said determining of the presence or absence of said predetermined frequency content can be determined from the vaule of said FOM.

2. The method of Claim 1 wherein said method further comprises the step of:
comparing said FOM with a predetermined FOM
associated with said predetermined frequency content.

3. The method of Claim 1 wherein said method further comprises the step of:
forming an averaged ECG signal by averaging preselected ones of said digital ECG signals before performing a FFT on at least a portion of said averaged ECG signal.

4. The method of Claim 3 wherein said ECG
signals comprises X, Y, and Z signals and said method further comprises:
forming an averaged ECG signal for each of said X, Y, and Z signals;

X9a forming an ECG template;
screening said digital ECG signals through said template to obtain said preselected ones, said step of forming a template further comprising:
selecting model X, Y, and Z signals;
determining the R R interval of a selected one of said model X, Y, and Z signals;
determining the fiducial point of said selected one;
determining the peak to peak amplitude of each of said model X, Y and Z QRS signals.

5. The method of Claim 4 wherein said step of screening said ECG signals comprises the steps of:
choosing candidate X, Y, and Z ECG signals for comparison with said template;
performing a multi-point cross correlation between said selected one of said model X, Y, and æ
QRS signals and the corresponding one of said candidate X, Y, and Z QRS signals;
selecting said candidate ECG signals as a portion of said preselected ones if said R-R
interval of said selected one of said candidate QRS
signal if within ?20% of said template value; said cross correlation coefficient is greater than or equal to 98%; and at least two out of three of the candidate QRS peak to peak amplitudes are the same as the template values.

6. The method of Claim 5 wherein said screening step further comprises selecting said candidate ECG signals as a portion of said preselected ones if and only if the X, Y, and Z ECG
signals before and after said chosen ECG signal are also acceptable when compared with said template.

X9a

7. The method of Claim 1 wherein said preselected portion comprises the terminal 40 milliseconds of the QRS complex.

8. The method of Claim 1 wherein said preselected portion comprises the ST portion of the ECG signal.

9. The method of Claim 2 wherein said step of determining the figure of merit comprises the steps of:
determining the normalized power level of the FFT output at 40 Hz; and determining the 60 dB area of the FFT output.

10. The method of Claim 9 wherein said method further comprises the step of forming a spectral plot of the FFT output before determining said FOM
therefrom.

11. The method of Claim 9 wherein said ECG
signals comprises X, Y and Z ECG signals and said step o f determining said FOM further comprises:
determining said FOM for each of said X, Y, and Z ECG signals; and forming the mean of said 40 Hz intercepts and said 60 dB areas form said X, Y, and Z FOM's.

12. The method of Claim 9 wherein said preselected portion of said digital ECG signals comprises the terminal portion of said QRS complex and said figure of merit associated with said predetermined frequency content comprises a 60 dB
area greater than 2400 and a 40 Hz drop less than 47 dB.

13. The method of Claim 9 wherein said preselected portion of said digital ECG signals comprises the ST portion and said figure of merit associated with said predetermined frequency content comprises a 60 dB area greater than 2500 X9a and 40 Hz drop less than 52 dB.

14. An apparatus for analyzing electrocardiogram (ECG) signals to determine the presence of a predetermined frequency content in the ECG signal comprising:
A/D means for converting analog ECG signals to digital ECG signals;
FFT means for performing a fast Fourier transform (FFT) on at least a portion of said digital ECG signals; and FOM means for determining a figure of merit (FOM) for said FFT, said FOM associated with the frequency content of said at least a portion of said digital ECG signals.

15. The apparatus of Claim 14 wherein said apparatus further comprises comparison means for comparing said FOM with a predetermined FOM associated with said predetermined frequency content.

16. An apparatus for analyzing X, Y, and Z
electrocardiagram (ECG) signals to determine the presence of a predetermined frequency content thereof comprising:
A/D means for converting analog X, Y, and Z
ECG signals to digital ones;
averaging means for averaging selected ones of each of said X, Y, and Z signals for forming averaged X, Y, and Z signals;
FFT means for performing a fast Fourier transform (FFT) on at least a portion of each of said averaged X, Y, and Z signals;
FOM means for determing a figure of merit (FOM) from the FFT of each of said averaged X, Y, and Z signals; and X9a comparison means for comparing said FOM with a predetermined FOM associated with said predetermined ferquency content.

17. The apparatus of Claim 16 wherein said apparatus further comprises:
means for selecting said selected ones of said X, Y, and Z signals for averaging comprising means for comparing said ECG X, Y, and Z signals with a predetermined template of a model ECG
signal.

18. The apparatus of Claim 15 wherein said at least a portion of said digital ECG signals comprises the terminal 40 milliseconds of a QRS
comples of said ECG signals.

19. The apparatus of Claim 15 wherein said at least a portion of said digital ECG signals comprises the ST segment of said ECG signal.

20. The apparatus of Claim 18 wherein said predetermined FOM comprises a 60 DB area greater than 2400 and a 40 Hz intercept less than 47 dB.

21. The apparatus of Claim 19 wherein said predetermined FOM comprises a 60 dB area greater than 2500 and a 40 Hz intercept less than 52 dB.

22. A method for analyzing electrocardiogram (ECG) signals to determine the presence or absence of a predetermined frequency content in a preselected portion of said ECG signals comprising the steps of:
converting analog ECG signals to digital ECG
signals;
forming energy spectra of said preselected portion of said ECG signals; and comparing a measure of the energy content of a first preselected portion of said each of energy spectra with a measure of the energy of a second X9a preselected portion whereby the presence or absence of a predetermined frequency content is determined.

23. The method of Claim 22 wherein the step of forming the energy spectra comprises the steps of:
performing a fast-Fourier transform (FFT) on said ECG preselected portion; and squaring the magnitude of said FFT of said ECG
preselected portion.

24. The method of Claim 23 wherein said ECG
preselected portion comprises a terminal portion of the QRS portion along with the ST portion of said ECG signals.

25. The method of Claim 24 wherein the step of comparing an energy measure of the energy content comprises the step of:
comparing the maximum amplitude of each frequency peak in said first preselected portion with the maximum amplitude of the largest peak in said energy spectra.

26. The method of Claim 24 wherein the step of comparing an energy measure of the energy content comprises the step of comparing the area of the energy spectra in said first preselected portion with the area of the energy spectra in said second preselected portion.

27. The method of Claim 25 wherein said first preselected portion comprises approximately the region from 20 Hz to between 50 and 70 Hz.

28. The method of Claim 24 wherein said first preselected portion comprises approximately the region from at least 20 Hz to 50 Hz and said second preselected portion comprises approximately the entire spectra.