Computer Programs in Biomedicine 6 (1976) 269-276 © North-Holland Publishing Company

PERIOD A N A L Y S I S OF T H E E L E C T R O E N C E P H A L O G R A M * , * * Bernard Allan COHEN

Emory University Regional Rehabilitation Research and Training Center, Center for Rehabilitation Medicine, Atlanta, Georgia30322, USA A program is described which utilizes the Digital Equipment Corporation LINC-8 computer with 4 K of memory to process the electrical activity of the brain into a set of meaningful statistics. The program provides an efficient method of obtaining spectoral information concerning the EEG without the associated cost or complexity of conventional power spectrum techniques. The combination of utilization of a small laboratory computer with an easily programmable method provides an approach for automated EEG analysis which is within the financial and technical scope of most small hospital EEG laboratories. Neurophysiology

Mini-Computer

Electroencephalogram

Period Analysis

1. Introduction

2. Mathematical description

While the technique of period analysis of the electroencephalogram (EEG) has been widely reported in the literature [ 1 - 1 3 ] , a systematic approach towards examining the mathematical properties and derivation of this method is frequently lacking. Since the purpose of the serial EEG analysis in the present study was to monitor long runs of data (e.g., several hours), it was important that efficient computational procedures be developed for deriving the parameters used to track the transient changes which occur. The parameters used in this study were derived from properties of the signal which are related to the signal's autospectrum, autocorrelation function, and its zero crossing behavior. The derivation o f these parameters as well as the mathematical properties and relationships among these parameters are shown in the following discussion.

It has been shown [14] that the relationship between the average number of zero crossings and the power spectral density (PSD) of a normally distributed zero mean process can be given by

* This work was originally in partial fulfillment of the requirements of the degree of Doctor of Philosophy in Biomedical Engineering at Marquette University, Milwaukee, Wisconsin, USA. ** Supported in part by Veterans Administration Center (Wood, Wisconsin) Project No. 1070-1 and National Institutes of Health Grants 5T01-GM-1051 and 1R01-GM-19680.

2°

f2K+2 p ( f ) d f (1) 0

where N K = average rate of zero crossings of the K-th derivative of the signal s(t), and P(f) -- power spectral density of s(t) integrated over the frequencies (f) o f interest. In period analysis only the values N K for K = 0, 1,2 are o f concern [3]. These correspond, in the signal s(t), to the points of baseline crossings, extremal points, and inflection points, respectively. Thus N O = the average rate at which s(t) passes through its zero mean value. N 1 = the average rate at which ~ t t) passes through its zero mean value.

d2s(t)

N 2 = the average rate at which ~ its zero mean value, dt

passes through

B.A. Cohen,Periodanalysisof the EEG

270

To demonstrate the validity of eq. (1) for EEG, it is necessary to compute the spectral moments defined by the right-hand side of eq. (1). This implies a lengthy computation to estimate the components of the Fourier spectra, and then integrating these spectra to compute their moments. However, moments of the power spectral density can be derived from relatively simple operations on the autocorrelation function. The autocorrelation function, 0(r), and power spectral density, P(f), for a finite epoch, have been defined by Rice [14] as ¢(r)

T/2 1 = T f s(t)s(t + r)dt -T/2

POe) 1 i

=

(2)

\ dr 2 ] .

T~2 s(t) e-i2"f dt-12

f

j

-T/2

(3)

Eqs. (2) and (3) are related by the transformations, as T-~ ~, given by:

(8)

For notational simplicity let ¢, _ d ¢ =~

d2¢ - dr 2

¢,, _

and

From the mean value theorem for differentiation [ 16] we have

+~

f

lag derivatives of the autocorrelation function [10]. Eq. (7) represents a major computational efficiency for, instead of computing the components of P(.f) and then integratingf2np(f) to obtain moments, all one needs to do is compute the values of 0(r) at increments of 2n+l of lag and then perform simple finite differencing operations on these values to obtain the desired spectral moments. A further computation efficiency is gained by observing that odd ordered derivatives are known a priori to be zero. Thus, in computing the second moment of P(f) we need to calculate only two correlation values and from eq. (7),

¢)(r)e-i2"frdr

_oo

(9, = 2 ? ¢(r) cos2rrfrdr

(4)

0

and using eq. (9) to compute ¢"(0) gives

and +~

0,0,=

= f em¢ 2"f7df

Zxr/2

0,0,

(lO)

_oo

= 2 ? P(f) cos 2 n f r d f (5) 0 It is noted that eq. (5) is independent of the amplitude probability distribution of s(t), so long as P(f) exists; thus, it may be used in general as a basis for generating moments of P(f) [ 15]. It is also noted that only the even moments do not vanish and thus if the even ordered derivatives of @(r) are considered, co

e2%(r!=2f dr2

0

(-1)n(2,f)2npOOcos2,frdf

¢ " ( 0 ) = 2 ¢ ( A r ) - (p(0)

Now substituting eq. (11) into eq. (8) we have

ff2p(f)df= 0

1 (2rr) 2

~b(O) - ¢(Ar) (zXr)2

(7)

From this it can be seen that the moments of the power spectral density are proportional to the zero

(12)

However, from eq. (1), subject to the previously stated conditions on s(t),

fP(f)df / (-x)n ~*ld2n~(o~ ~2(21r)2n] ~ dr 2n ]

(11)

(At) 2

(6)

and rearranging and setting r = O,

J f2nP(f)df= 0

but ¢(0) = 0 and ¢ ' 1 ~ } is given by eq. (9)so that eq. \ "¢-" / (10) reduces to

! o or using eq. (12) and eq. (7) (with n = O)

B.A. Cohen, Period analysis of the EEG

(-~12

(1

~Az)~ .

(14)

In the EEG data analyses of numerous investigations [ 1,2,3,10,17 ], it has been demonstrated that the second moment of the normalized power spectral density is accurately approximated by NO, the average zero crossing rate ofs(t). Saltzberg [3] has termed the parameter N O the "major period count" and concluded that it provides an important measure for signal tracking because its computational simplicity makes large scale EEG analysis of spectral moments quite practical. Additionally, Burch et al. [2] have studied the validity of assuming the EEG to be Gaussian. They found the square-wave train generated by the primary wave of the EEG to consist of a homogeneous population of major periods. It has been shown in eq. (12) that the second moment of the power spectral density may be approximated by two values of the autocorrelation function and that this result is independent of the constraints imposed on eq. (13) which relates zero crossings to the second moment. A demonstration of the validity of zero crossing approaches to spectral moment calculations (and the subsequent demonstration that the EEG is an excellent realization of a normally distributed zero mean process as required by Rice [14] for eq. (1) to hold) is found by comparing the values obtained from both methods of computation, i.e., zero crossing and power spectrum. This comparison was shown by Saltzberg et al. [3] in both sleep and drug studies. Saltzberg's experiment was designed to ascertain whether significant EEG alterations were reflected in spectral moments and zero crossing tracking parameters. He defined the square root of the second moment of the PSD as the "gyrating frequency" since it has dimensions of frequency and since the equation representing it has the same form as that for the radius of gyration in mechanics [3]. Thus from eq. (14) 1 1/_ NO (15) fg - 2rr&r V 1 - % ( & r ) - 2 ' where Cn(Ar) is the normalized autocorrelation function at the first increment of delay and Ar is the interval between data points. In addition to showing the efficacy of zero-cross parameters, Saltzberg also demonstrated, with long data samples, that the gyrating frequency of each analysis epoch was consistently

271

equal to half the measured average zero crossing rate for the corresponding epoch. Finally, he concluded that the constraints on EEG signal statistics justify using zero crossing parameters to estimate and monitor spectral variance. In 1971 Saltzberg and Burch summarized this previous finding and stated that '~the EEG is a process for which it is legitimate to use average zero crossing rates to calculate moments of the P S D . . . for the purpose of monitoring long-term changes in the statistical properties of the EEG" [10].

3. Discussion The preceding derivation shows the purely mathematical relationship between zero-cross techniques and moments of the power spectral density of the signal. However, when the encephalographer clinically examines the signal his analysis is not based on some mathematical relationship, but rather (in part) on the synchrony of the wave and on certain well-defined landmarks which are explicit frequency components. This is done based on the fact that a certain frequency will occupy a certain amount of time on the chart. Perhaps the only mathematical relationship which enters the analysis, and then only subconsciously, is f = 1/T

where f i s frequency and T is the time interval between repetitive waves. The encephalographer makes very little attempt to analyze the relatively "infrequent" complex portions of the EEG signal (signals riding on signals). In the EEG, since one does not know the source, the signal can only be approximated by a periodic process. One may see vibrations in the signal but really does not know what they are. Nothing is known about its harmonic nature. When any kind of automated analysis (zero-cross, FFT, PSD, etc.) is done, it is always done as a means of identifying parameters that have been empirically chosen by encephalographers to describe a phenomena. These phenomena have been given names so that they may be categorically studied, "time-wise", within this phenomena; there is little attempt to address the problem of the existence of fundamental or hamonic frequency components. If one calls a component DELTA, it is because the encephalographers do. Thus, the problem of sig-

272

B.A. Cohen, Period analysis of the EEG

nal interpretation is really one of visual pattern recognition; an attempt is made to learn a given array of patterns so that they are recognizable in the future. That being the case, the problem is one of feature selection or feature description of some phenomena in an existential fashion. These various features (really they are a phenomena in a phenomena) are called Delta, Theta, Alpha, etc.; indeed, they have been defined in frequency terms which really "correspond" to time intervals. Automated zero-cross analysis, which really is an elaboration of the visual method used by encephalographers, is merely a means of determining features which can be related to clinical pictures or pre-defined patterns, so as to enable the acquisition of clinical correlates. There is no desire or attempt in zero-cross analysis to desynthesize the waveform with the intent of resynthesizing it in the future (though it has been shown [3,17] that it can be done). Rather, it was the intent and purpose of this research to develop the simplest technique possible to extract the most relevant information from the EEG. This goal is in line with one of the most fundamental precepts of pattern recognition, i.e., features that one uses are not necessarily required to have any bearing or real relationship to the process being studied, as long as they exist in a pattern space to be used with a transformation which allows classification of the patterns. Subsequently, any transformation used (FFT, PSD, zero-crossing, etc.) is for the purpose of scaling the data and maneuvering it into a form that can be quantified, with respect to the solution of the actual problem being studied. The application of period analysis in this research fulfills the afore-mentioned requirements and purposes. Nevertheless, questions frequently arise regarding the validity of using period analytic methods when one considers the existence of complex waveforms such as are found in the EEG. It is emphasized that, for the most part, this type of activity (a signal riding on another signal) is quite negligible, especially when the analyses are done over very long periods. Saltzberg has demonstrated [3] that the bandwidth constraints of the EEG are such that there is a high degree of correlation between the interval measurements of period analysis and relative amplitude within each primary zero crossing interval; the longer the time period between axis crossings, the greater the signal amplitude.

Burch has found [17] that, in general, the superimposed frequency is so low in amplitude in relation to the dominant frequency that the resultant is more a distortion of wave shape than a "mixture"; the superimposed frequency itself is not usually apparent. By virtue of the characteristics of the EEG signal, there will rarely (perhaps 90 percent confidence or better) be erroneous results. Rather, information may be missed (this can be avoided by applying period analysis to the signal's first and second derivatives [1, 2,3,10,11,17] ). Errors of m isclassification will usually not occur. In effect what this means is that period analytic data preserves the necessary information for effectively reconstructing the EEG waveshape [3,17]. Finally, it is appropriate to examine some of the advantages and disadvantages of using period analysis methods for interpreting EEG waveforms. One such problem which has been noted is that period analysis can have the disadvantage of showing higher frequency noise as lower frequency sine waves are digitized and analyzed [2]. To avoid this problem, several precautions were taken. First, the period was always measured from an average of the two sample points about the zero axis, thus avoiding the problem of triggering on noise. Second, a dead-band was employed so that any reproduced signals smaller than 1 percent of the maximum expected peak amplitude (100/iV) would be rejected. Finally, the signals were band-pass ffiltered before being digitized thus reducing the noise in the frequency range outside the band of interest. From the point of view of waveform analysis, period analysis' most significant features are its sensitivity to changes in EEG waveshape and its amenability to automatic data processing techniques. It has been noted that frequency analysis of EEGs based upon the FFT, while it has been used extensively, has the disadvantages of sensitivity to non-stationary trends in the EEG and rather extensive computation time. Alternatively, period analysis, based upon the first three terms of the Gram-Charlier series [1], has the advantage of relative insensitivity to non-stationary effects and extreme simplicity of computation. In general, whenever the topic of pros/cons comes up (in EEG analysis), the answer must be specified in terms of what methods are possible, given the limitations of time, money, instruments and knowledge. In this research, while the FFT or PSD are both possible, the amount of data to be analyzed with the

B.A. Cohen, Period analysis of the EEG

existent equipment would render either of these techniques prohibitive. Finally, for the purpose of examining long-term changes in the EEG, it would be inefficient to use a process if it takes ten minutes to analyze one minute of data. A clinician can glance at the record and specify nearly as much information as may be obtained from the computer using elaborately expensive, mathematically pure and exact methods such as the FFT or PSD.

4. Program structure A general flowchart diagram of the computer program is found in fig. I. The program begins by retrieving several significant parameters dealing with the length of the epoch, the tape drives to be used for data storage, the number of epochs to be analyzed, column headings, etc. This information was originally generated in a monitor program used for input/output and operator interaction. The main-line program begins by clearing a variety of epoch counters for each of the eight EEG bands (CLEAR). Sub-routine CROSS then begins sampling the signal from the A-to-D converter. An appropriate zero crossing is evaluated for rejection from the threshold levels that were specified by the operator. If an appropriate set of data points is found such that a time interval between successive zero crossings can be calculated, then control is directed to sub-routine FILTER. Here the decision is made as to which category the appropriate time interval belongs in. An interrupt routine (INTRPT) is included to eliminate artifact information from the summing registers. Subroutine CHKTIM and ADDTIM are used to calculate the total time remaining in a given epoch prior to program transfer to a storage routine. Sub-routines TOOUMB and FRMUMB are used to transfer data from the storage registers to or from the appropriate upper memory bank. Additionally subroutines DUMPMB, DUMPPT, LOADMB, and LOADPT are utilized to transfer the contents of the upper memory bank storage either to or from the appropriate tape drives. After each memory bank transfer, sub-routine CHKEND is used to evaluate whether or not the memory bank is indeed full and a tape transfer is necessary. Following the entire analysis,

273

the tapes containing the data are rewound to access the header information at the beginning of the tape segment containing that data. This is accomplished with sub-routine RETHED and RETDAT while subroutines HEDOUT and DATOUT are used to print out the appropriate headings and data. An octal-todecimal conversion sub-routine is included (DECTYP) as well as sub-routines to stop and start the analog tape recorder (STREC and GOREC). Calculation of the total amount of time remaining in an epoch utilizes sub-routine DIVIDE to calculate a floating point number. Also under operator control is the option of obtaining the printout in the form of IBM coded paper tape which is accomplished through sub-routine IBMCON. Finally, sub-routines are included to control the spacing for hard copy printout (PAPOUT) and for IBM coded paper tape output (ATOUT). The last sub-routine of the overall program is a routine utilized to evaluate the timing mechanism inherent in the computer (CLOCK). This short subroutine is probably the only part of the program which would have to be changed if the overall program were to be utilized on any other Digital Equipment Corporation devices since this particular part of the program is inherent with the LINC-8.

5. Sample run A sample output from the program may be seen in fig. 2. Here the reader will note the appropriate headings give codes for the hospitalization and the patient's problem as well as counters representing the amount of frequency activity (total number of zero crossings per epoch) in each of eight EEG bands. These EEG bands were determined on the basis of the half-cycle timing windows noted in fig. 3. The sample output additionally shows the clock timing for each of the epochs. Each epoch in this particular case was set at 60 sec (600-1/10th sec). It is noted that because of artifact information, on the average about 58 sec of information was generally obtained.

6. Hardware and software specifications The program is coded in LINC assembly language for use on the LINC-8 computer system. The program

Yes *

®

®

Fig. 1. Flow diagram of LINC-8 period analysis program.

Last Epoch Yet?

el, 11~9i;

Gj -

Set cl,

t

(7:z-3

J

Yes

Yg8

Yes

Yes

Yes

Yes

OW.EGA

(

p

~Yes

No

r

~J

.®

DELTA

&

oa

i-.)

B.A. Cohen, Period analysis of the EEG

1278

0576

0022

As written, the system requires only 4 K of core and at least one magnetic tape drive (LINC or DEC tape drives). The time required to run the program is only a function of the amount of data to be input. As set up in its current form, the clock routine enables the user to process data at ten times real time. This effectively means that analog data which required one hour to record now takes approximately 6 min to process. Once all of the data has been processed the printout time remains strictly a function of the output device being utilized. In the case of this research, the output device was the standard ASR33 teletype.

116I 1099 1254 1524

05 74 ~579 0571 0573 0570 0572 0573 0578

00~ 0024 0025 0026 0027 0028 0029 0050

7. Availability of program

~I01075 ~CH

1425-16-0961

10109175

F0C 59 M NOEM 1 0 1 1 0 1 1 5 0P30 AM ~ 0 OF LEAD8 T3-(VERTEX)

O040/0200 ARTIO DFLTA THETA ALPHA SIGM~ GAMMA KAPPA I.AMDA OMEGA 0000 0038 0243 0232 0168 0205 01~8 0120 0084 0000 0026 02=6 0248 0161 0201 0168 C I 4 1 0097 0000 0038 0234 0248 O205 0203 0 t 6 0 0103 0071 OOO0 0059 0222 0250 0160 0200 0141 0146 0088 0000 0027 0225 0278 0203 0221 010a 0102 0070 0000 0040 0 1 8 ] 0256 0217 0240 0206 0156 0097 0000 OOal 0101 02AI 0182 0246 02O4 0157 0094 0000 0030 0194 0248 0216 0258 0212 0150 0003 000¢ 0029 0198 0257 021A 0241 0193 0154 0088 0000 0028 0186 0259 0202 0255 0229 0169 0098 0000 0041 0161 0244 0200 0240 0212 0157 0110 0000 0041 0171 0262 0180 0243 0197 0164 0080 0000 0068 OIA7 0220 0148 0207 0182 0141 0105 0001 0092 0140 C I 9 2 0144 0174 0132 0088 0079 0000 0O94 0169 0208 OIA2 0166 0120 0 l l 8 0059 0000 0104 0158 0204 0153 0151 0110 0106 0065 0000 0076 0166 0227 0186 0170 011~ 0097 0073 0001 0073 0164 0205 0209 0192 0151 C105 0056 0000 00 77 0166 0231 0142 01~3 0140 0110 0055 0000 0085 0156 021~ 0198 0179 0129 00~9 0054 0000 0077 OIYA 0257 0166 014~ Olla 0094 0~59 0000 0062 0146 0259 0225 02a5 0153 0126 0064 0000 0076 0145 0247 0206 0188 0124 0096 0 0 7 5 0000 0079 0171 0275 0193 01~8 0102 0057 0058 0000 0050 0177 0269 0201 0199 0170 0112 0076 0000 0042 0149 0520 0272 0211 0157 0115 0062 0000 0035 0190 0 2 ~ 3 0221 0221 0164 0 1 1 6 0070 0000 0058 0166 0250 0240 0208 0145 O l l C 0065 0000 0069 0165 0248 0104 0230 0121 0072 0068 0000 0070 0184 0274 O 2 [ I 0147 0086 0086 0044 o000 0060 0158 0 3 0 0 026l 0145 0 1 3 2 0074 0053 0000 0058 0143 050~ 0275 0 t 8 5 0082 0070 0037 0000 0071 0176 0310 02~'~ 0J44 O08? 0O68 O0~R 0000 0068 0167 0310 0214 0154 0088 0059 0039 0000 0089 01~3 0283 0195 0120 0077 0 0 4 3 00~1 0000 ~09~ 0174 0257 0175 0109 0066 OOaY 0035 00o0 0000

0096 0095

0178 0177

0269 0287

0147 Olg4

010~ 0115

0062

0000

0107

0183

0233

0171

0116

0073 0042

0000

0091

0190

02fl5

0140

0125

0050

275

TOTAl_ ¢I.~CR FPOCH 1218 056~ OOPI 1288 0563 00C2 1260 057O 0005 1286 0568 0004 1290 0572 0005 1353 0567 0006 1356 0568 0007 1401 0864 0008 1374 0565 0009 1426 0566 0010 1405 0561 0011 1347 0563 0012 1218 0558 0015 1041 0563 0014 1076 0572 0015 1029 0570 0010 1122 0572 0017 1155 0565 0018 1124 0570 0019 III8 0514 C020 1090 O575 0021

13t0

1242 I1~7 II06 1223

0578

0051

It76 1123 1106 1021 0~1

0592 05~2 0583 0586 0587 0587 o584 0549

0032 0053 0034 0035 0036 0057

0039 00al

0028 0026

0964 0~78

0033 0041

0032 0024

0~37 o~a6

0544

0038

0039 OOaO

A manuscript of the program is available from the author in either ASCII coded paper tape or hard copy printout. Additionally a DEC tape or a LINC tape of the manuscript can be obtained.

Fig. 2. Sample output of period analysis program. will run equally well on systems set up for PAL assembly language programs which include many of the other systems available from the Digital Equipment Corporation. The only changes that may be necessary would be in the timing clock routines of the program.

Acknowledgement The author gratefully acknowledges the technical assistance of Mrs. LaDonna Compton.

References EEG Activity

Frequency (HA)

Half-Cycle Timing Window (m. See.)

DELTA

[ A, ~t ]

0.S to 4.0

THETA

[O,0]

4.0 to 8.0

125 to 62.5

ALPHA

[ A, a]

8.0 to 12.0

62.5 to 41.7

SIGMA

Iv,

a]

12.0

to 16.0

4!.7

GAMMA

IF,

"~]

16.0

to 20.0

01.3 to 25.0

KAPPA

[ K, ~ ]

23.0

to 24.0

25.0 to 20.85

LAMDA

[ A, X ]

24.0

to 28.0

OMEGA

[ ~2, ~ ]

28.0 to 32.0

|030

to 125

to 31.3

[ 1 ] B. Saltzberg and N.R. Burch, IRE Trans. Med. Elec. 8 (1957) 24. [2] N.R. Burch, W.J. Nettleton, J. Sweeney and B.J. Edwards, Ann. N.Y. Acad. Sci. 115 (1964) 827. [3] B. Saltzberg, R.J. Edwards, R.G. Heath and N.R. Burch, Synoptic Analysis of EEG Signals, Rochester Conference of Data Acquisition and Processes in Biology and Medicine 5,267-307, Rochester, NY., 1966. [4] F. Ferillo, C. Rivano, G. Rossandini, G. Rossi and C. TureUa, Electroenceph. clin. Neurophysiol. 27 (1969) 700.

20.85

to 17.85

17.85 to 15.62

Fig. 3. EEG major periods and the corresponding half-cycle timing windows. The LAMDA band is intentionally misspelled to facilitate five-letter computer codes.

[5] T.M. Itil, D.M. Shapiro, M. Fink and D. Kassebaum, Electroenceph. clin. Neurophysiol. 27 (1969) 76. [6] H. Legewie and W. Probst, Electroenceph. clin. Neurophysiol. 27 (1969) 533. [7] B. Hjorth, Electroenceph. clin. Neurophysiol. 34 (1973) 321. [8] T.M. Itil, J. Nerv. Ment. Dis. 150 (1970) 201. [9] R. Roessler, F. Collins and R. Ostman, Electroenceph. clin. Neurophysiol. 29 (1970) 358. [10] B. Saltzberg and N.R. Burch, Electroenceph. clin. Neurophysiol. 30 (1971) 568.

276

B.A. Cohen, Period analysis of the EEG

[ 11 ] A.J. Welch, Aerosp. Med. 42 (1971) 601. [12] T.M. Itil, B. Saletu and S. Davis, Biol. Psychiatry 5 (1972) 1. [13] B. Hjorth, Electroenceph. clin. Neurophysiol. 34 (1973) 321. [14] S.O. Rice, Mathematical Analysis of Random Noise, in: Noise and Stochastic Processes, Edited by N. Wax (Dover Publications Inc., New York 1954) pp. 133-294.

[15] J.S. Bendat and A.G. Piersol, Random Data: Analysis and Measurement Procedures (John Wiley and Sons, New York. 1971). [16] C.R. Wylie, Jr., Advanced Engineering Mathematics (3rd ed.) (McGraw-Hill, New York, 1966). [17] N.R. Burch and H.E. Childers, Physiological Data Acquisition, in: Physiological Aspects of Space Flight, Edited by B.E. Flaherty (Columbia University Press, New York, 1961) pp. 195-211.