Citation |

- Permanent Link:
- https://ufdc.ufl.edu/UF00082461/00001
## Material Information- Title:
- Automated detection and quantification of petit mal seizures in the electroencephalogram
- Creator:
- Principe, Jose Carlos Santos Carvalho, 1950-
- Publication Date:
- 1979
- Language:
- English
- Physical Description:
- viii, 364 leaves : ill. ; 28 cm.
## Subjects- Subjects / Keywords:
- Amplitude ( jstor )
Digital filters ( jstor ) Electroencephalography ( jstor ) Microcomputers ( jstor ) Seizures ( jstor ) Signal to noise ratios ( jstor ) Signals ( jstor ) Sine function ( jstor ) Statistical discrepancies ( jstor ) Waveforms ( jstor ) Dissertations, Academic -- Electrical Engineering -- UF Electrical Engineering thesis Ph. D Electroencephalography ( lcsh ) Epilepsy ( lcsh ) - Genre:
- bibliography ( marcgt )
non-fiction ( marcgt )
## Notes- Thesis:
- Thesis--University of Florida.
- Bibliography:
- Bibliography: leaves 333-345.
- General Note:
- Typescript.
- General Note:
- Vita.
- Statement of Responsibility:
- by Jose Carlos Santos Carvalho Principe.
## Record Information- Source Institution:
- University of Florida
- Holding Location:
- University of Florida
- Rights Management:
- Copyright Jose Carlos Santos Carvalho Principe. Permission granted to the University of Florida to digitize, archive and distribute this item for non-profit research and educational purposes. Any reuse of this item in excess of fair use or other copyright exemptions requires permission of the copyright holder.
- Resource Identifier:
- 023383227 ( ALEPH )
06808347 ( OCLC ) AAL4060 ( NOTIS )
## UFDC Membership |

Downloads |

## This item has the following downloads: |

Full Text |

AUTOMATED DETECTION AND QUANTIFICATION OF PETIT MAL SEIZURES IN THE ELECTROENCEPHALOGRAM BY JOSE CARLOS SANTOS CARVALHO PRINCIPE A DISSERTATION PRESENTED TO THE GRADUATE COUNCIL OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA 1979 ACKNOWLEDGEMENT S The author would like to thank his supervisory committee for their aid, guidance and criticism during the preparation of this dissertation. Especially the author would like to thank Dr. Jack Smith for his willingness to share facilities, personal contacts and responsibilities. It was an unconventional but effective learning experience, which went beyond the restricted field of EEG studies. The friendly relationship created during these years is deeply appreciatï¿½d. There are many others who contributed their knowledge, skills and time for this dissertation to become a reality. The author would like to thank Sonny, Yvonne, Shiv, Nate, Dan, Alan and Salim for their large, small, direct or indirect contributions. Special thanks go to Fred for his constant help*in the data collection phase. The author also would like to thank the Department of Electrical Engineering of the University of Florida, Comissao Permanente INVOTAN (NATO--Portugal), and the Institute of International Education for the financial support received throughout his graduate studies. The author sincerely thanks his family for their support and confidence, especially his wife for her constant companionship in spite of the stress created by the long working hours devoted to the laboratory. iii TABLE OF CONTENTS Page ACKNOWLEDGEMENTS . . ii ABSTRACT . . . . . . . . . . . . . . . . . . . . . . vi CHAPTER I LITERATURE SURVEY . . . . . 1 Spectral Analysis . 1 Time Domain Approach . . . . 18 Selection of the Method of Seizure Detection . . . . . . . 25 II MICROCOMPUTER BASED DIGITAL FILTER DESIGN . . 30 Preliminary Considerations . . . . . 30 Design Criteria . . . . . . . . . . 35 Filter Transfer Function . . . 38 Lowpass to Bandpass Transformations . . . 41 Transformations to the Z Plane . . . . 43 Finite Length Effects of the Implementation . . . 53 Practical Considerations . . . . . . 84 III A MODEL FOR PETIT MAL SEIZURES AND ITS IMPLEMENTATION . . 94 Detection Problem . . . . . 94 A Petit Mal Seizure Model . . 98 Definition of the Implementation Scheme . . . 106 System Implementation . . . . . . . 136 Testing of the Program . . . . 157 IV SYSTEM EVALUATION AND PRESENTATION OF RESULTS . 165 Description of Data Collection . . . . . 165 System Evaluation . . . . . . . . . . . . . . 170 Statistics of Seizure Data . . . 187 Analysis of the False Detections and Proposed Improvement in the Detector . . . 224 Proposed System Utilization and Concluding Remarks . . . . . 233 APPENDIX I CHEBYSHEV FILTER DESIGN . . . . . . . Chebyshev Polynomials . . . . . . . . Automated Design . . . . . . . . . . . Filter Design Program . . . . . . . . II DIGITAL FILTER IMPLEMENTATION . . . . Preliminary Considerations . . . . . . Implementation of the Filter Algorithm III FLOW CHARTS AND PROGRAM LISTINGS . . . REFERENCES . . . . . . . . . . . . . . . . . . BIOGRAPHICAL SKETCH . . . . . . . . . . . . . . Page 239 239 242 244 250 250 253 265 333 346 . . . ï¿½ï¿½ï¿½ ï¿½ ï¿½ï¿½ï¿½ ï¿½ï¿½r ï¿½ï¿½ï¿½ Abstract of Dissertation Presented to the Graduate Council of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy AUTOMATED DETECTION AND QUANTIFICATION OF PETIT MAL SEIZURES IN THE ELECTROENCEPHALOGRAM By Jose Carlos Santos Carvalho Principe August 1979 Chairman: Jack R. Smith Major Department: Electrical Engineering This dissertation deals with the automated analysis and quantification of Petit Mal (PM) paroxysms in the human electroencephalogram (EEG). A petit-mal detector was built to help neurologists evaluate drug effectiveness in the treatment of PM epilepsy. The decision of basing the detection scheme in a microcomputer provided completely new directions for the design process because of the versatility and great computational capabilities of the machine. The problem of nonrepeatability of characteristics in the analog filters was avoided by replacing them by digital filters. A microcomputer based digital filter design procedure was developed, taking into consideration the computation speed, the filter internal gain and the noise characteristics of different topologic structures. To yield sufficient output signal to noise ratio a 16 bit microcomputer and a 12 bit A/D converter were utilized. With the procedure, a fourth order bandpass filtering function can be accomplished in less than 300 ps. Time domain parameters were selected to arrive at a PM seizure model which would be applicable to the detection of classical PM and PM variant epilepsies. The parameters were translated into electronic quantities and implemented in a microcomputer. For the first time in the automated analysis of EEG, a microcomputer based system was built for the processing of one channel of data, using a completely digital, real time, detection scheme. All the programs which comprise the detector were written in assembly language, and care was taken to reduce the interactions between program modules. The top-down program approach was utilized. The system was tested and evaluated with data from epileptics collected in the Veterans Administration Hospital, Neurology Service, from an ongoing drug study. A telemetry link was used to record the EEG, allowing the patient to move at will in a 3 x 5 meters room. Six patients who showed the largest number of seizures were selected for the evaluation. A total of 70 hours of data was analyzed. The agreement in the detections with the human scorer, for seizures (Sz) greater than three seconds, was 86 percent. If the sorting of seizures in the groups 3 < Sz vii < 10 sec and Sz > 10 sec is taken into consideration, the agreement decreased to 77 percent. Besides outputting the information about number of seizures, their duration and time of occurrence, the detector also presents a detailed analysis of the detection parameters. These results, the first of their kind ever reported, include the mean and variance for the PM recruiting period, the half period of the slow waves, the delay between the slow wave and the spike, and the amplitude of the filtered spike and slow wave. They show the constancy of the period measures inter and intra seizures for the same patient and the higher variability of the amplitude measures. Preliminary results correlating pairs of the detection parameters are also presented for three patients. The automated detection and quantification of PM seizures, using the present system, are critically brought to a focus. viii CHAPTER I LITERATURE SURVEY Two research methods have been applied to automated electroencephalogram (EEG) studies: spectral analysis and time domain analysis. It is the purpose of this chapter to survey the basic different techniques and assumptions involved and to present a better picture of the advantages and limitations of each method. Our attention will be mainly focused on techniques used or related to detection of pathological paroxysmal events in the EEG, and the literature referenced herein denotes this emphasis. Spectral Analysis Before addressing this subject, let us first summarize the general assumptions made about the EEG that validate the application of the technique. The theory behind the frequency domain approach using a nonparametric model, as in conventional power spectral analysis, handles the EEG signal as a stochastic process. The statistical properties of the process essentially influence the analytical results. One major difficulty is exactly the different opinions about the stationarity and Gaussian behavior of the EEG process which are prerequisites to carry out spectral analysis and to use the power spectrum as a sufficient descriptor of the EEG parameters. At the present time this is a very hot topic where one can read reports that show the short time stationarity of the EEG for sequences below 20 seconds for a probability of error of 10 percent (Cohen & Sances, 1977; Kawabata, 1973; Saltzberg, 1972). However, the work of Elul (1967, 1969) and Dumermuth (1968) points out the nonstationary behavior of EEG sequences as short as 4 seconds. Another important characteristic, albeit less addressed, is the Gaussianity of the EEG. The work of Elul (1967, 1969) and McEwen and Andersen (1975)'suggested that the behavior of the EEG followed a Gaussian distribution 32 to 66 percent of the time depending upon the behavioral state of the subject (active mental task to relaxed state, respectively). It also displayed a spatial distribution with closer Gaussian characteristics in the occipital leads. Some of the discrepancies that are reported in the literature can be explained by the different sampling rates used in the A/D conversion, since oversampling may distort the local statistical properties of the digital waveforms (Persson, 1974). Another reason, and maybe more plausible, just demonstrates the variability of the EEG with subject and behavioral state. In practice the stationarity requirement can be loosened to stationarity during the observation period, and if the process departs slightly from Gaussianity, ancillary parameters like skewness and kurtosis could be employed along with spectral analysis. If stationarity holds, but the segments are highly non-Gaussian, higher order statistics as bispectra (Dumermuth et al., 1970) could be needed. Undoubtedly, the early work on EEG quantification of Grass and Gibbs (1938), Walter and coworkers (1963, 1966, 1967) brought improvement to the field of electroencephalography, which at the time was faced with amplitude measurements of randomlike activity. Maybe more important yet, it brought a consistent analysis technique, well established from other areas with an enormous amount of computational power. After the introduction of the Fast Fourier Transform (FFT) algorithm in 1965 by Cooley and Tukey, this power increased manifold. At this point, however, the means began to obscure the basic constraints of the method. Let us review then the methods of spectral analysis. One of the basic results is to obtain the power spectrum, i.e., an estimate of the mean square value or average intensity of the EEG as a function of frequency. It displays the decomposition of the total variance into contribution from the individual frequency bands. To obtain the power spectrum one can apply Fourier transform to correlograms (indirect method, Blackman and Tukey, 1958), or take directly the Fourier transform of the data and square its absolute value--periodograms (direct method, Dumermuth & Keller, 1973; Matousek, 1973), or use autoregressive techniques (Gersch, 1970; Mathieu, 1970). After the introduction of the FFT the direct method is by far the most widely used. To find the periodogram, the FFT of the finite EEG sequence is taken. The finite time series can be thought of as the multiplication of the infinite time series X(t) with a rectangular window function of duration T. In the frequency domain this corresponds to the convolution of the spectrum of X(t) with a sinf/f which introduces finite resolution and leakage in the spectrum of X(t) and therefore affects any spectral quantity extracted from it. The finite resolution and leakage are always coupled together and derive from general relations between the time (duration) and frequency (bandwidth) domain transformations (Thomas, 1969). It is only possible to improve one at expense of' the other. One family of approximation functions that yield the optimum compromise (prolate spheroidal functions) are quite difficult to implement and are seldom used (Temes et al., 1973), so more readily available functions are preferred. The best solution is to deal with each problem separately: use a window (Hammings, Kaiser, Tukey, Barlett, see Childers & Durling, 1975) to smooth the data which will improve side lobe rejection; and increase the record length to improve the resolution. The last requirement can only be accomplished at the expense of longer computation times and is directly constrained by the stationarity of the data. This concept deserves a further explanation because sometimes the increase of the record length, by appending zeros at the end of the data record, is thought to increase the "resolution." As a matter of fact, the procedure displays more points of the discrete Fourier transform (DFT), and so a better approximation to the Z transform is achieved. However, the Z transform remains the same (as long as the record length is constant), so the initial resolution is kept unchanged. Another problem encountered is the estimation of the power spectrum. Let us resume the estimation process. To compute the essential statistical properties of a Gaussian random process the more general way is to evaluate the autocorrelation function by ensemble average. The method is computationally expensive since various sample sequences of the process have to be available. If the random process is ergodic (Lee, 1960; Thomas, 1969), then it is correct to compute the autocorrelation function by a time average only over one sequence. There are fewer memory requirements involved, but as the properties of the random process are to be inferred from the ones calculated in the sample sequence, consistent estimators have to be used (i.e., estimators that are unbiased and whose variance decreases with the number of data points). It is very easy to find consistent estimators for the autocorrelation functions (Oppenheim and .Shafex:, 1975). The problem arises when the goal is to estimate the power spectrum. It is not true in general that the Fourier transform of a consistent estimator of the autocorrelation function (periodogram) will produce a good estimate of the power spectrum (Barlett, 1953; Jenkins & Watts, 1969). When an infinite number of points is used to evaluate the periodogram, very rapid fluctuations are present in the estimate of the spectrum. This fact shows that increasing the number of data points of the observation sequence will not change the statistical properties of the estimator since the individual significance of the data points remains the same. To increase the stability of the spectral estimates, the number of independent observations have to be higher or in some way the statistics of the observation have to be improved. These two points lead to the two main methods of increasing the stability of the estimation--the averaging of periodograms (Barlett, 1953), and the smoothing by means of a window (Welch, 1967). Let us describe briefly the two methods. The Barlett method averages the periodogram over certain frequency bands. To achieve this the data are divided in k segments; the periodogram is calculated for each segment; and the final estimate is the average of the various periodograms. The procedure leads to a bias estimate, but the variance decreases with k (overlapped segments can be used at expenses of lower independence of individual segments and smaller resolution). The number of points on the power spectrum is decreased by k. The Welch method windows segments of the data before the periodogram is calculated. Generally the triangular window is used to ensure positive estimates of the power spectrum. The methodology is the same as for the Barlett method, and Welch was able to show that the variance also decreased with k and the reduction on.resolution was also present. Both methods give similar results, so generally the Barlett method is preferred to save one step (Dumermuth, 1968; Matousek & Petersen, 1973; Hagne et al., 1973). Although the variance of the estimate can be made arbitrarily small, there is a trade-off for resolution. Therefore, if an unknown spectrum is going to be estimated, great care must be exercised not to miss peaks in the power spectrum and at the same time trust the peaks displayed. In EEG the modified periodogram as an estimator of the power spectrum has been widely used. One of the pioneering works was done at UCLA to investigate the EEG activity of astronaut candidates by D. O. Walter and coworkers (Walter, 1963; Walter, Rhodes, Brown & Adey, 1966; Walter, Kado, Rhodes & Adey, 1967; Walter, Rhodes & Adey, 1967). Also, Dumermuth and coworkers in Switzerland (Dumermuth, Huber, Kleiner & Gasser, 1970; Dumermuth, 1968) use extensively this technique. It has been applied in automated sleep scoring (Walter et al., 1967; Caille, 1967; Rosadini et al., 1968), age-dependent EEG changes (Hagne et al., 1973; Matousek & Petersen, 1973), schizophrenia (Giannitrapini & Kayton, 1974; Etevenon et al., 1976), studies of EEG in twins (Dumermuth, 1968), studies of EEG background activity (Walter et al., 1966; Zetterberg, 1969; Gevins et al., 1975), brain lesions (Walter et al., 1967), and also drug evaluation (Matejcek & Devos, 1976; Kï¿½nkel et al., 1976). Sometimes the power spectrum is not the only parameter used. Dumermuth et al. (1972) uses the bispectrum, that is, the Fourier transform of the second order autocorrelation function R(tl, t2, tl+t2) to investigate coupling between EEG frequencies. Incidentally, the bispectrum analysis shows that, unless for the ac and 8 frequency bands, the function is practically zero, as would be the case for a Gaussian random process. The coherence spectrum --ratio of the square of the cross spectrum over the spectra of the individual waveforms--is also used to investigate correlation between two EEG channels (Dumermuth et al., 1972; Walter et al., 1967). The most bothersome question in the application of the periodogram as an estimator is the fact that there are no "best criteria" to determine the combination, window duration-overlap for a particular data sequence. That is one of the reasons why other more consistent methods (in the sense that an error criterion can be defined and leads to a minimization), as the autoregressive (AR) and autoregressive moving average (ARMA), have been recently introduced (Fenwick et al., 1967; Isaksson, 1975; Gersch & Sharp, 1973). The assumption behind their application is that the EEG sequence selected is a sample from a stationary time series. Therefore, it can be modeled as the output of a linear system (if transients are excluded) with a white noise input. The output can then be interpreted as the linear superposition of the natural modes of the system. The eigenvalues extracted from the parametric time series model of the EEG are the characteristic frequencies and their associated damping factors. The damping reflects the extent to which the EEG will exhibit an oscillatory or random appearance. Unlike nonparametric (conventional) spectra analysis, the relative dominance of the frequencies in the EEG can be associated with the ordering of the magnitudes of the eigenvalues. In the AR model, the linear system transfer function is approximated by an all-pole function. The two problems are the determination of the filter poles location and the order of the approximation polynomial. There are three main minimization procedures to fit the AR model to the data: one evaluates the filter coefficients as the maximum likelihood estimate when the input to the linear system is white Gaussian noise and the output is the data sequence available (maximum likelihood method). Another evaluates the filter coefficients in such a way that the mean square error between the next sample and the predicted one, given the data sequence, is minimized (linear prediction). Finally, the other calculates the filter poles in such a manner that the estimated values of the sequence maximize the entropy of the autocorrelation function (maximum entropy). The three methods have been developed under different constraints, but recently van den Bos (1971) and Smylie et al. (1973) proved their equivalence. It is also interesting to note that all the methods arrive at the same matrix equation involving the autocorrelation function (Yule Walker equation (YW)). To estimate the filter coefficients, the YW equation can be solved. Then the autocorrelation functions for various lags have to be computed and a matrix inversion performed. Algorithms due to Levinson (1947) or Durbin (1960) are generally preferred to solve it recursively. A technique due to Burg (1967) suggests another estimation of the coefficients that does not require prior estimate of the autocorrelation function. It fits successively higher order prediction error operators to the input data by convolving the filter in both forward and backward direction. The square of the two error series are added to obtain the error power,and the filter coefficients are determined by a minimization procedure. A recursion relation is also available (Burg, 1967, 1972; Anderson, 1974). Actually the two different techniques (YW method and Burg) use the data in two different ways. The YW assumes that the data is zero beyond the limits of the observation window. The Burg technique extends the estimation of the autocovariance coefficients beyond the initial number of points in such a way that the entropy of the autocovariance function is maximized at each step. As can be expected, the two techniques have quite different properties. The estimator by the YW method has the attractive property that its mean square error is generally smaller than other estimators (Jenkins & Watts, 1969). However, the estimates are very sensitive to rounding (Box & Jenkins, 1970) mainly for random processes that display peaking power spectrum, which means that the resolution obtained is not always enough. This is to be expected since the YW technique effectively windows the data. The Burg technique does not display this shortcoming of lack of resolution. It has been shown (Pusey, 1975) that it can resolve two tones arbitrarily close if the S/N ratio is high enough. However, the variance of the estimator is larger than for the YW and does not decrease to zero monotonically. The correct identification of the order of the AR model that approximates the data is vital in the computation of the power spectrum (Ulrich & Bishop, 1975). The criterion generally used is Akaike's Final Prediction error (FPE) that is defined as the mean square prediction error (Akaike, 1969, 1970). The YW estimate of the order of the process using this criterion tends to be conservative but with a small variance. The Burg estimate of the order displays a large variance, so generally some upper bound must be imposed in the search. Another approach is to use the first minimum in the estimated error (Ulrich & Bishop, 1975). While the AR model is an all-pole approximation to the random process, the ARMA model is a pole-zero approximation. Computationally it is the more expensive (in time and memory dimensions) and more difficult to implement, but the two techniques for AR computations could be modified for this more general case (Treitel et al., 1977). However, as an infinite AR process can approximate any ARMA process, this property of the AR model is preferred. The FPE criteria to determine the order of the AR model can still be used, but there is a tendency to overestimate the order (Ulrich & Bishop, 1975). After briefly reviewing the three approximation models (MA, AR, ARMA) to estimate the power spectrum, a question must still be answered. From the sequence of the random process available how should one decide what model to use? This question (identification problem) has not yet been-answered in general, and only in some cases the physical knowledge of the generation process of the data has helped. It is surprising that we have not seen published any work in this direction by the users of these techniques in EEG. On the contrary, the approximation properties of the AR model have been exclusively used. There are reports in the / literature (Treitel et al., 1977) that clearly show that this can lead to inappropriate approximations, which means that there is no single correct technique to calculate the spectrum in the absence of knowledge about the physics of generation. In EEG the autoregressive model has been used by Fenwick et al. (1969) to predict evoked responses. Pfurtscheller and Haring (1972) used AR models to attempt data compression. Gersch & Sharp (1973) used AR models for multivariate EEG analysis. Wennberg & Zetterberg (1971) used ARMA models to extract alpha, beta, delta, and theta band intensities for automatic EEG classifications. Mathieu (1970, 1976) used AR models in sleep scoring. Jones (1974) compared the window lag estimation with AR estimator to determine power spectral and coherence functions in the neonate EEGs. Gersch and Yonemoto (1977) compared the AR and ARMA models for the power spectrum estimation of EEG. The same authors (1973) also used parametric AR models to apply the Shannon-Gelfand-Yaglom amount of information measure in sleep scoring. In epilepsy AR models have been used to-analyze the ictal event. Gersch and Goddard (1970) used the partial coherence among time series to extract information about driving. Tharp (1972) and Tharp and Gersch (1975) use a similar procedure to determine seizure focus. Herolf (1973) and Lopes da Silva et al. (1975) also use an AR model to perform seizure detections by inverse filtering. Another concern faced in spectral analysis is the presentation of the results. Due to the specific techniques involved, the results need to be further processed to be readily interpreted by the neurologist. The main goals are data reduction and readability. The compressed spectral array of Bickford et al. (1973) has been extensively used and summarizes pretty well the changes of EEG with some external factors (different behavioral states, tracking of alpha waves with light stimuli, hemispheric symmetry, acute slow wave abnormalities). The canonogram method of Gotman and Gloor (1973) consists of a spatial arrangement of polygons of different sizes related to the ratio of (delta + delta)/(alpha + beta) which is thought to be a good descriptor of localized brain abnormalities. It is considered to be a good descriptor for detection of slow wave abnormalities and asymmetries in the EEG if a preprocessing for gross artifacts in the raw data is performed. The results of optimum filtering have also been applied to extract certain features from the EEG. One of the earliest techniques was matched filtering (Smith et al., 1969; Zetterberg, 1973; Saltzberg, 1971). A matched filter is the filter that maximizes the signal-to-noise ratio, since it transforms all the available information (energy) of the input signal in a voltage at a specific time T. The prerequisites for its construction or modeling are the a priori knowledge of the waveform shape and the stationarity of the background noise. A matched filter is in fact a correlation detector, and a relatively simple way of implementing digitally one (for a white Gaussian noise) is to select a wave pattern, store it backwards in memory, and perform a correlation with the input (template matching). If the noise is not white, we have to prewhiten it and separate the implementation into two steps: first, divide the Fourier transform of the desired wave pattern by the power spectrum of the noise and second, inverse transform the resultant spectrum to obtain the template that can be stored in the computer memory. For EEG, prewhitening is a must, since the background activity (noise, if a particular pattern has to be detected) is not white. Matched filtering has been used for spike detection in EEG by Saltzberg (1972), Zetterberg (1973), and more recently by Barlow and Sokolov (1975) and Eftang (1975). Another result borrowed from communication theory is inverse filtering. The inverse filter is generally applied to extract information about the arrival times of the individual components of a composite waveform (Childers & Durling, 1975; Robinson, 1967). It is required to know the shape of the waveform and to make the stationarity assumption on the background noise. The filter is nonrealizable (the output is theoretically a delta function), and great care must be taken in the required inversion of the Fourier transform of the waveform not to divide by zero. Inverse filtering also deteriorates the S/N ratio, and a compromise has to be reached between resolution and signal-to-noise ratio. In EEG inverse filtering is used in a different way. Barlow and Dubinsky (1976) use inverse filtering to generate the impulse responses of bandreject filters (EEG - power spectrum -* square root + individualize band of interest subtract from uniform spectrum (with amplitude equal to maximum) + IDFT (to get impulse response)). This procedure is highly sensitive to truncation effects, and so poor filters result. Lopes da Silva et al. (1975) use a linear prediction scheme. He implements the inverse filter of the all-pole approximation model of the EEG and monitors the error between the predicted value and the real data points. Deviations above a prescribed threshold are correlated with nonstationarities of the input EEG data. Very similar techniques are used in EEG adaptative segmentation procedures (Bodestein & Praetorious, 1977) which could eventually lead to transient event detection. A combination of the nonstationary information of the phasic events and the sharpness of the spikes was combined by Birkemeier et al. (1978) to increase the resolution separation of the clusters of epileptic transients and background activity. The system requires settings for each patient and works better in epochs with few spikes (otherwise they may bias the estimation). It is sensitive to artifacts. Etevenon et al. (1976) use null-phase inverse filtering to deconvolve the EEG and then detect fast activity. His method is equivalent to the one just described (Lopes da Silva), but it is executed in the frequency domain. From the EEG sequence the power spectrum is calculated through the periodogram, and after smoothing and whitening the amplitude spectrum is obtained by evaluating the square root. Then the inverse of the spectrum is computed, multiplied by the Fourier transform of the incoming EEG, and the inverse FFT taken. This signal accentuates sharp transients contained in the EEG, and after threshold detection they can be extracted. To cope with the EEG nonstationarities, recently there have been reports on Kalman-filtering (Isaksson & Wennberg, 1976). If the stationarity constraints are removed from the noise in the ARMA model, it can be shown that the best fit in the least mean squared sense to the time series is achieved by the Kalman filter (Kalman, 1960; Isaksson, 1975). So, this filter is a generalization of the Wiener filter for nonstationary data and is closely related to autoregressive techniques. It consists basically of the ARMA fit to the time series, but due to the inclusion of a feedback loop with variable gain (Kalman gain) the locations of the poles of the system are allowed to vary. The key parameter in the design is the setting of the Kalman gain since it weighs the new value of the residual signal with its past values and changes accordingly the coefficients of the all-pole filter. The gain controls then the adaptability of the filter. If the gain is set at zero, the Kalman filter is just an ARMA fit for the incoming data and is nonadaptable. If the gain is increased, then we have adaptable properties and the tracking is improved with the gain. One important consequence of the adaptability of the filter is to avoid the smoothing effect present in other spectral methods. Incidentally, by comparing the autocorrelation function obtained by any of the stationary spectral methods with the autocorrelation function estimated from the output of the Kalman filter, some knowledge about the limits of the record length to preserve stationarity can be inferred. Isaksson and Wennberg (1976) point out that about 90 percent of the records of 20 s are stationary which agrees with some of the previous observations. At the present time, only Isaksson (1975) and Isaksson and Wennberg (1976) have used this technique. Time Domain Approach If one can relate the degree of quantization to the maturity of a field of science, electroencephalography is a newborn. The lack of quantization has two main consequences: the inexistence of objective criteria to characterize the EEG phenomenon that weakens any purely analytical research method (like spectral analysis); the resort to a multitude of soft criteria (i.e., not unique), each based on empirically derived parameter values that try to translate the highly nonlinear way the electroencephalographer reads the EEG. Visual analysis is a true time domain method of detection of structural features like diffuse and localized change in frequency and voltage pattern, changes in the topographic distribution and in the interhemispheric symmetry, sharp and rhythmic activity, paroxysms and unstable or irregular time course. The basic patterns and clinical EEG features are composed on the following step-by-step procedure: the analysis of the period and amplitude.forms the wave concept that is associated with the gradient to give the simple grapho-element. These can be integrated in a sequence to produce the complex grapho-elements from which, using superposition, the notion of pattern is acquired. The variability in time and hemispheric location gives rise to the time domain structure and topography of the EEG. The job of the electroencephalographer is to individualize these elements and compare them to his acquired set of "normal activity" in order to make the diagnosis. Of course, the boundaries are quite fuzzy, highly subjective, and not always consistent (Woody, 1968; Rose, 1973). The techniques are then dependent upon the specific detection task desired, and no simple overview is possible besides the translation to electronic terms of the decomposition process: analysis of the wave's amplitude in certain frequency ranges, which implies broad bandpass filtering, zero crossing, and threshold analysis. Depending upon the criteria formulated the next step can be the testing of the grapho-elements throughout further processing (sharpness or repetition period, for instance) or testing of a pattern by coincidence logic. What makes the process somewhat erratic is that there is no methodology beyond the extensive comparison with the results obtained by the electroencephalographer and accordingly modify the implementation or the criteria until an agreement is reached. Due to the diversity of techniques, the detection of abnormal brain activity (spikes) will be emphasized. One of the first criteria to detect spikes was the sharpness (Buckley et al., 1968; Ktonas, 1970; Walter et al., 1973; Gevins et al., 1975). This was achieved by monitoring the first and second derivatives of the EEG and comparing them with a fixed threshold set for each patient. The parameter was found generally unsatisfactory since differentiation increases the energy of the incoming signal at high frequencies (i.e., extends the bandwidth) and the detection system becomes very sensitive to muscle artifacts. The same basic idea was further improved by Carrie (1972a, 1972b, 1972). He established a moving threshold set by background activity. Although the system was still sensitive to high frequencies, since they biased the threshold, he reported better results. Gevins et al. (1975, 1976) also used the curvature parameter (second derivative), but to improve the system performance the duration and frequency of occurrence are also introduced in the detection criteria. The threshold is automatically set for each patient by an empirical algorithm. A more complex model of abnormal spikes which included different slopes for the leading and trailing edges of the spike, plus a parameter related to the time it takes the wave to reach maximum slope, was introduced by Ktonas and Smith (1974), and it seems to describe fairly well spike activity. Smith (1974) used some of these parameters to implement a spike detector that gave good agreement without adjustment of detector thresholds. Basically, it consists of a set of gated monostables with "on" times related to the leading edge first derivative (negative excursion), sharpness of the apex, and trailing edge first derivative (positive excursion). The system requires fairly high sampling rates and is somewhat sensitive to the amplitude of the incoming signal and to muscle artifact, although in a smaller degree than the systems that computed the second derivative of the EEG. Another detection scheme uses the amplitude-frequency characteristics of the background activity and of the spikes to define detection boundaries in a plane. It is a modified scheme of the period amplitude-analysis of Leader et al. (1967) and Carrie and Frost (1971). Ma et al. (1977) used it to establish optimum decision boundaries that turn out to be nonlinear to accommodate the dependence of amplitude threshold with duration. It is necessary to use each patient's background activity and an average value taken from a population of subjects to estimate the boundary (assuming the clusters of normal and abnormal activity are jointly Gaussian distribution with equal a priori probabilities). The results seem promising but are dependent on the availability of the sets and are computationally complicated and expensive (in time and memory dimensions). Harner and Ostergen (1976) also used a similar approach called sequential analysis to display the amplitudeduration of the EEG but preserving the time information. His methods show the clustering of paroxysmal spike activity and some patterns as spike and wave. It is possible to choose between the display of half or full waves and have a better qualitative view of the rising and falling phase of spike and wave activity. At this time, however, the system is only qualitative (Harner, 1973) since the definition of boundaries for spikes or spike and wave complexes need medium computers and more resolving power. One important theoretical defect of amplitudeduration techniques, still not solved, concerns the mixing of two frequencies having similar amplitudes. Pattern recognition was also applied to the detection of paroxysmal events in EEG (Serafini, 1973; Viglioni, 1974). The theoretical method for the determination of the relevant parameters by performing data reduction (linear principal component analysis (Larsen, 1969) and nonlinear homeomorphic procedures (Shepard, 1966)) have hardly been applied to the EEG (Gevins et al., 1975). The parameters chosen are once again few in number and related to visual analysis: amplitude and the average value of the waves, number of zero crossings and mean value of zero crossings. The experimental results show that the criteria are satisfactory, but it requires a training set, is critically dependent upon the length of the normalization interval, and requires a fairly large number of decomposition (orthonormal) functions. Following the same line of pattern recognition Matejcek and Schenk (1974), Remond and Renault (1972), and Schenk (1974, 1976) applied a vectorial iteration technique to decompose the EEG. First, the maxima and minima of the envelope of waveform are estimated, and the higher order component is taken as the mean.value; the second step is an iteration procedure and uses the result of the preceding analysis. The process is repeated until a featureless waveform is obtained. With these wave components the half waves analysis is used to extract further information. The appealing property of this method is the relatively small computation involved (compared to the FFT) and the strict resemblance with the visual (pattern) analysis. Gotman and Gloor (1976) describe and evaluate a set of parameters to describe the EEG at the waveform level. The technique can therefore recognize phasic interictal events in the EEG. Its main application is to generate information about epileptic focus and degree of abnormality of a particular record. No mention is made on its use in the detection and quantification of generalized seizures. The other work reviewed fell in one of the previously described methods. Vera and Blume (1978) use the derivative method to analyze on line 16 channels of EEG. Chick et al. (1976) use a scheme similar to Smith (1974). As far as petit mal (PM) activity is concerned, Jestico et al. (1976) use a bandpass filter 2-4 Hz to detect the slow wave component and measure the duration of the paroxysm. Kaiser (1976) used the duration of the spike and of the slow wave monitored at a certain voltage level to detect the PM activity. In a PDP-12 Ehrenberg and Penry (1976) used zero crossing information and a measure of integrated amplitude from a combination of 4 EEG channels to detect PM activity. The system stores the duration and time of occurrence of the seizure. The overall agreement, consensus versus machine, is 85 percent. It is, however, interesting to note that in the paper's extensive discussion, no mention is made on artifacts that may cause potential problems like chewing and body movements, which suggests that the investigators had available rather clean data. The system is reported to be sensitive to slow wave sleep. Carrie and Frost (1977) also described a spike and wave detector. It is composed of a spike detector, using the sharpness criteria (Carrie, 1972a), an EMG detector and an information of the amplitude of the background activity (Carrie, 1972b). One channel is analyzed on line,and the system stores the duration of the paroxism and its time of occurrence. The agreement is very high for seizures greater than 3 sec (85 percent) but drops off to 25 percent for Sz between 1-3 sec. The patients were reported to have well-defined PM paroxisms and were free of medication. From the EEG data available to us it seems that the large amount of high energy artifacts that resemble the PM pattern when the patients are in an unrestrained environment dictates the use of more powerful pattern recognition algorithms. Johnson (1978) implemented in a microcomputer a PM detector based on the repetition properties of the pattern, i.e., spikes followed by slow waves. He used two analog filters (22-45 Hz spike, and 1.5-4 Hz slow wave), followed by threshold logic to input the basic elements to the microcomputer memory, where the pattern recognition takes place. The system had no false detects in selected epochs, but was very sensitive to irregularities of the pattern in the middle of the ictal event. Selection of the Method of Seizure Detection The requirements put on the analysis method related to the specific problem of detection of paroxismal events are the following: 1) The detection shall be performed in real time in one channel, using a microcomputer. 2) The detector shall have high resolution capabilities and be insensitive to high energy artifacts. 3) It shall also be able to quantify in detail the paroxisms. The computation time constraint is probably enough to make the selected choice of analysis obvious, since FFT algorithms that run real time in simple microcomputers, with workable resolutions, are not known to exist. But even in the affirmative case, it seems that spectral analysis is not tailored to event detection in the EEG because the spectrum produces a smoothing (leakage/resolution). To compensate for it, longer sequences and/or special techniques (parametric spectral analysis) need to be used, which increase the computation time. After obtaining the spectrum, some type of pattern recognition must be utilized to judge about the presence of certain frequency components coupled with the event, which incidentally must assume a complete knowledge about the background spectrum and the relations time patternsfrequency patterns. The question arises, why not work with the time patterns to begin with? Another problem is how can the information from frequency analysis be translated to the clinician who sees time patterns and wants information such as time of occurrence (ï¿½l sec), duration, amplitude, etc.? It seems a much more natural choice to use time domain techniques. The other technique reviewed is matched filtering, which is the optimum detector to extract patterns from a noisy background, when the noise is stationary and the patterns known. The problem with its application to EEG is the variability of the patterns (around 10 percent, Smith, 1978) and the drastic change in the power of the noise (background activity), which deteriorates the detector performance. Yeo (1975) was able to show that a zero crossing detector, although not optimum, performed better when the patterns were allowed to vary. The false alarm rate of the detector is also independent of the noise power. The technique which detects nonstationarities in the EEG have some drawbacks besides being computationally expensive. The inverse filtering can not differentiate among types of nonstationarities. Therefore, a preprocessing or other ancillary methodology needs to be developed. Kalman filtering is quite sensitive to the gain setting. For high gains small changes can be exaggerated, and so one is faced with the problem of calibrating the filter, which is usually done with a control data set (artificially generated). The direct dependence of the gain on the type of data, the dependence of the error signal on the power of the input sequence, and some transient behavior of the system are shortcomings. It seems that the popularity of frequency domain techniques in EEG detection stems from the availability of standard software packages and a fairly well developed (but often forgotten) mathematical methodology. Up to the present, one of the drawbacks of the time domain approach is the implementation medium (hardware) which requires a full engineering development, not always accessible to the clinicians. However, the present innovation trend in microcomputer systems may very well bring the possibility of standardization through the software implementation of detectors, at low price. To build an electronic detector where some parameters (amplitude in a prescribed frequency range, zero crossing, sharpness) have to be monitored and decisions made, the omissions created by the poor quantization (definitions) of the EEG process have to be filled. That is the reason why soft (e.g., nonunique) criteria have to be created and the result of the detection always compared to the EEG scoring by the experts. The attractiveness of the procedure resides in the small number of assumptions needed on the data and the simple and fast implementation of most of the parameters. On the other hand, there is no theoretical analysis backing it; therefore, the methodology almost exclusively used is trial and error. The process is long and sensitive to errors in the choice of the key parameters or misjudgements of the criteria. From the time domain techniques known there was some hesitation between employing the rather well-established techniques developed by Smith et al. (1975) to analyze sleep EEG and the ones proposed by Gotman and Gloor (1976) of implementing definitions in the raw EEG. Since the implementation of a completely digital petit mal detector around a microcomputer was, per se, a task which involved quite a few unknown steps, it was decided to use the technique which was best documented. However, it is expected that this work will give quantitative bases for comparing the two different methodologies, with respect to the appropriateness of using prefiltering in the detection schemes of short nonstationarities, like spikes. Chapter II will be devoted to the design of digital filters using microcomputers and study in detail the hardware selection of A/D converter number of bits and computation wordlength to obtain a certain output signal to noise ratio. In Chapter III the model for the petit mal activity utilized in the detection will be presented, along with its implementation on a microcomputer. The system's testing will also be explained there. As this research work is primarily engineering oriented (a new instrumentation system is designed around a new model of the petit mal activity), the system's evaluation is presented in Chapter IV. Preliminary data demonstrating the high resolution capabilities of the detector and its use in quantifying the petit mal seizure data will also be presented. CHAPTER II MICROCOMPUTER BASED DIGITAL FILTER DESIGN The purpose of this chapter is to analyze the tradeoffs in digital filter implementation using microcomputers and arrive at a comprehensive design procedure. The design approach proceeds from the study of the digital filter noise factor to the choice of a microcomputer wordlength and A/D converter precision, so that a specified output signal to noise ratio could be obtained. This approach is quite general and may be valuable even if future technologies will make practical other implementation media (e.g., bit slice microcomputers). Preliminary Considerations For quite some time the EEG research group at the University of Florida has been developing nonlinear methods for the detection of EEG waveforms. Basically the detectors include bandpass filters followed by zero crossing and threshold analysis of the filtered data to extract information about period, amplitude and number of in-band waves in the raw EEG. The superiority of this detection method for EEG activity has been established when the variability of the patterns is fairly high (greater than 10 percent-Smith, -1978). The bandpass filters have up to the present been implemented in analog form. The remainder of the pattern recognition algorithm is implemented digitally. Analog filters have the disadvantage of being very sensitive to changes of the component values produced by environment parameters like temperature, humidity, and aging. Although there are design methods that can minimize the sensitivity to component changes (e.g., leap frog), they are more difficult to design since in such structures there is as-ubstantial:number of feedback loops and the change of one filter parameter implies the modification of a large number of filter components. Another problem to which analog filters are sensitive is the tolerance in component values. To design two filters with identical characteristics some or all of the components have to be hand matched. However, this procedure only applies to that particular point in time, since the matched components can have different aging coefficients and so evolve differently in the long run. The identical filter characteristics are stressed here because, at the present time, a sensitivity analysis which can be applied to the detection of EEG waveforms is not available. However, it is known from experience that very small differences in the filter parameters and/or pattern criteria produce drastic changes in the detection (Smith, 1979). Therefore, the only way similar detection characteristics can be guaranteed between two systems is to match every detection stage. The substitution of the analog filters by digital ones is also a natural extension of the present instrumentation since the other functions of the pattern recognition algorithm are already digitally implemented. The repeatibility and uniformity of characteristics will be therefore ensured. The implementation of digital filters could be accomplished basically in two different ways: building special purpose machines to implement the filter algorithm (hardware) or making use of general purpose computers and writing adequate software to accomplish the filtering function. The two techniques possess different properties that are worth comparing. The main advantage of the hardware realization is speed since the processor's architecture is intentionally adapted to the special type of processing, e.g., multiplications and additions (Gold et al., 1971; de Mori et al., 1975). Generally, the processor is microprogrammed and the arithmetic unit is implemented in fast ECL logic. There are some minor modifications to the basic procedure, the use of Read Only Memories to implement the multipliers being the most interesting (Peled & Liu, 1974). However, the development of a special processor is very expensive and implies the availability of specialized laboratories and the work of a diversified group of researchers. On the other hand, the software implementation of digital filters can be performed by anyone who masters a computer language as long as the recursion relation is available. This fact per se has the potential to expand the field of application of digital filter signal processing. This approach is also less expensive since the machine may be a general purpose computer, and so there is no extra cost involved in the new application. The drawback is the computer's fixed architecture which is not adapted to calculate recursion relations. For a lot of applications the requirements of the implementation can be met with mini-computers and even microcomputers. The principal elements of the requirement set are speed, type of arithmetic, wordlength. Digital algorithms work in digital representations of real world (analog) signals. It is well known (Thomas) that to make a digital representation of an analog signal unique (i.e., the analog signal can be reconstructed again from the digital samples) a maximum frequency must be assumed in the signal spectrum. This is the same as saying that the time domain representation of the signal is restricted to change at less than a prescribed rate, governed by the time/ frequency inverse relationship. The theorem that relates the uniqueness of the digital representation and the sampling frequency is the Nyquist theorem, and it imposes a lower bound on the sampling frequency (sampling frequency equal to two fM). This theorem assumes two impossibilities: first, that the signal spectra are frequency band limited, and second, that the past and future history of the signals are completely available for reconstructing the analog signal. Here we will not discuss the validity of the hypothesis but will acknowledge the fact that the signal's digital representation is an approximation. To describe sampled waveforms in the time domain a sampling frequency higher than the Nyquist rate should be utilized. Therefore, for real time processing a sampling frequency between 200 and 300 Hz seems adequate for most applications (Smith, 1979), which means, in the worst case, that the computation time per sample must be less than 3.34 ms. An analysis of the clock frequencies and the instruction cycles of today's microcomputers shows that hundreds of operations can be performed in this interval. It seems a comfortable margin when one compares this number with the apparent simplicity of the recursion algorithm--a few multiplications and additions. "Apparent" is stressed because we are generally. led fo think in nterm.offloating point arithmetic. It is a good exercise to estimate the enormous number of operations needed to perform a floating point operation when using numbers greater than one as sole representations. It turns out that the computation time becomes prohibitive if floating point is utilized. There are other types of arithmetic, like the block floating point (Oppenheim, 1970), which do not seem to bring any advantage for our application. The most severe limitation in the use of fixed point arithmetic is the small dynamic range, since the maximum representable number is 2b, where b+l is the number of bits in the processor wordlength. The wordlength is also coupled to the precision of the representation because the numbers that represent the constants (and signals) must be quantized to fit the wordlength. For instance, if an 8 bit processor is used, the constants are represented only by two hexadecimal digits. The choice of software implementation of digital filters in microcomputers imposes stringent constraints which must be carefully examined to ensure that the analog input signal to noise ratio is not degraded. Nevertheless, the use of microcomputers is thought a good choice due to their cost, size, availability, and the ease with which the filtering function is actually performed, requiring primarily software knowledge. In this specific application, as rather sophisticated pattern recognition algorithms will be necessary, the microcomputer will also be time shared, to accomplish these functions. Design Criteria The design of digital filters can be broadly divided in two phases. The first is related with the determination of the filter algorithm and the second with the implementation of the recursion relation. Generally they are taken independently since the type of problems encountered are quite different. In the design of the filter algorithm a machine of infinite wordlength is assumed. The goal is to arrive at a recursion relation that better fits the filter's frequency (or time) specifications. In the second phase the objective is to choose filter implementations which optimize, in some sense, the finite wordlength effects of any practical processor. The most important problems of the finite length effects are the finite dynamic range of the computation, which may cause overflows, and the finite precision of the constants, arithmetic, and input data, which produces error that may degrade the signal to noise ratio. It turns out that the two design phases are not completely independent,as will be shown. The microcomputer implementation is expected to be very sensitive to the above-mentioned parameters since the arithmetic must be fixed point (small dynamic range which increases the probability of overflow) and the wordlength is relatively small, giving a heavier weight to the finite effect errors (also called ï¿½roundoff: errors). For this reason it was thought convenient not to separate the two design phases in order to have a better perspective of the interactions between a specific recursion relation and its implementation. As various routes may be taken for digital filter design, a criterion to compare different design procedures (hence different implementations) is necessary. Here design procedure means selection of appropriate transformation methods to arrive at the algorithm for the filter transfer function, with maximum simplicity in terms of number of operations, providing at the same time good overflow and noise properties. To design a filter with a specific output signal to noise ratio using the minimum requirements of computation wordlength and processing time, the following factors need to be analyzed: 1) Filter transfer function 2) Filter internal magnification 3) Noise characteristics of filter topology 4) Processing speed. The problem as enunciated is relatively different from the design procedure commonly found in the literature, since there the optimization is studied with respect to only one parameter (most generally to sensitivity to roundoff), neglecting any other considerations. It turns out that the optimum structures possess a much higher number of multipliers which will mean slower computation time in our application. Examples are as follows: The synthesis of infinite impulse response filters (IIR) with low roundoff noise based on state space formulations has been presented by Hwang (1976) and Mullis and Roberts (1976). Their structures have N2 more multipliers than the canonical structures, where N is the filter order. Another example will be the structure proposed by Barnes et al. (1977), which is free from overflow oscillations, but requires N+1 more multipliers than the canonical structures. Multiplier extraction was also considered (Szczupack & Mitra, 1975), but the structures obtained lack the optimum sensitivity characteristics. Another interesting idea is presented in Chang (1978), where the decrease in the roundoff noise is accomplished by feeding back the discarded bits, properly weighed, to the input of the adder following the place where the product is quantized. This procedure doubles the multiplier number. Szczupack:and Mitra (1978) propose the reduction of the roundoff noise by ensuring the presence of zeros in the noise transfer functions. This procedure leads to zeros that are complex, therefore requiring extra multipliers for its realization. Many more examples could be given, but the picture remains unchanged; decrease of roundoff noise means more multipliers. For a microcomputer implementation of digital filters multiplication is by far the most time-consuming operation, and it is therefore necessary to minimize its use. It seems appropriate to say that for microcomputer implementation a suboptimal solution would be the best, weighing the effect of the increase in computation time and the decrease in the variance of the output noise. This analysis is not known to exist, so the procedure chosen utilized canonical structures as well as slight variations in their topology (adders as a variable). Filter Transfer Function The characteristics of the filters presently used for sleep EEG studies have been obtained through trial and error. Basically the filters are low Q, second order bandpass filters (Smith, 1978). The present design uses a slightly underdamped frequency response. It is hard to qualify the desirable properties of the filters since there is no general theory of EEG detection. The only possible criteria are to understand the filter function in the detection scheme and hopefully arrive at some rules that will serve as guidelines for filter selection. The purpose of using EEG filters in sleep studies is to attenuate out-of-band activity to enable further signal processing (e.g., zero-crossing). On the other hand, the filters shall not mask out-of-band activity so that it will look like in-band activity for the rest of the processing algorithm. Every filtering function produces a certain masking, since the filter output is the convolution of the filter's impulse response with the input. However, the weight of the smoothing is controlled by the Q of the filter, hence the use of low Q filters. At present no other characteristics of the filters, like phase distortion or group delay, have been brought into the picture of filter selection (Smith, 1978). The in-band amplitude characteristics have also been neglected due to the wide variability of the amplitude of the EEG. The above-mentioned analysis refers to EEG sleep studies. For the epilepsy application the knowledge is still less quantitative due to less experience. From the study of the petit mal pattern it can be concluded that the fast attenuation characteristics are of paramount importance for the slow wave filter (activity as high as 5 Hz) since the spike principal component can be as low as 11 Hz, roughly one octave away. Otherwise, the spike energy will degrade the periodic appearance of the slow wave filter output, preventing the use of full cycle detection, which has been considered the best (Smith, 1978). After this quick analysis it can be said that the overwhelming filter characteristics are the good out-ofband attenuation and the low Q. The Q is a factor controlled by the filter parameters and not by the filter type. Therefore, a "good" EEG filter shall display fast out-of-band attenuation. With this constraint in mind, the filter type chosen will be the Chebyshev, since it displays the optimum out-of-band attenuation rate for a certain order, when all the zeros are assumed at infinity (Weinberg, 1962). Elliptic filters have similar attenuation characteristics but are more difficult to design due to the proximity of the poles and zeros. Also for a digital filter implementation this fact means complex zeros and subsequently more multiplications. To design the digital filters and use the present experience with analog filters it was decided to accomplish the design in the S plane and use one of the mapping rules to the Z domain (impulse invariant or bilinear). It was thought important to have an automated facility to design the filters from the frequency domain filter characteristics, i.e., center frequency, bandwidth, ripple in-band and out-of-band attenuation rates. The use of Chebyshev polynomial is suited for this goal as shown in Appendix I. The procedure to determine the Chebyshev lowpass goes as follows:(Principe et al., 1978): 1. From a given ripple factor (E) in the passband, the lowpass filter poles can be obtained from (1). 2. To determine n, the polynomial order from the attenuation in dB/OCT, (3) can be used. 21 2k+1) s = -sin h4 sin( ï¿½f) w k 2 2n + j cos h42 cos(21T) (1) 1 1 1+ - + + 1 )n sin h)2 E E 2 (2) 1 1 1+ + i + 1 + 1 + n cos hE2 2E log( E n = (3) log(2 + /v) Lowpass to Bandpass Transformations After designing the lowpass filter transfer function, the bandpass can be obtained using one of the standard transformations--the narrow-band or the wide-band (Blinchikoff, 1976). In the narrow-band, the lowpass filter poles are transformed by the relation w2-w = 2 1 i ï¿½jw0 SBP 2 LP where w2, w1 are the upper and lower cut-off frequencies and w0 the center frequency. For the wide-band the poles pairs are given by a>0 a<0 w0 [-+bï¿½j(Y8-a) w0 [- a-b+j(l+a)] Y 2 2 A = 1---(a - ) 4 y y w0 [-YZ+bï¿½j (Y+a)] w0 [-a-bj j-a )] o 2 2 y2 B = 2 a =/ A2+B2 +A 22 b A+B -A 2 and the lowpass poles are -aï¿½j. Here the transformation of the poles is stressed since it will avoid polynomial factorization to obtain the bandpass digital filter in a cascade form. The two transformations possess quite different properties, the most important for this application being the different Q of the poles. Looking closely at (4) it can be seen that the bandpass poles lay on a parallel to the jw Pair 1 Pair 2 where axis. As the Q is proportional to the ratio w/a, the higher frequency poles always possess greater Q. For narrow-band filters the difference is small, but it can become quite appreciable for wide-band filters. Hence, the narrow-band does not transform the time characteristics (overshoot, group delay) of the lowpass design (Blinchikoff, 1976). However, it is widely used due to its simplicity. With the wide-band transformation, the poles of the bandpass filter lay in a straight line from the origin, ensuring that the ratio w/a is constant. The impulse response of the bandpass filter is a time scaled version of the lowpass. The disadvantage of the wide-band is a more tedious design procedure. It will be shown later that the different characteristics of the transformations will be valuable for the control of the frequency response of the digital filter and the finite length effects of the implementation. Transformations to the Z Plane Now that the bandpass pole locations are known, it will be seen how the Z plane pole pairs can be obtained. Basically there are two different types of mapping rules from the S to the Z plane. One, the rational transformations, map the entire jw axis onto the unit circle. There are several ways this can be accomplished (a series expansion of eZt), but the lowest order approximation is the bilinear transformation given by 2 z-1 s + (6) T z+l The particular type of isomorphism (an infinite length line mapped into a finite line) implies a nonlinear mapping, in the sense that a Az has different correspondence to a Aw, depending on the value of z. As a matter of fact, the transformation follows a tangent rule as can be easily shown (Oppenheim & Shafer, 1975). Therefore, to get the correct frequency characteristics in the Z domain, the S plane filter critical parameters (center frequency, bandwidth) have to be properly prewarped (Childers & Durling, 1975). Nevertheless, a piecewise linear response is mapped onto a piecewise linear response, preserving the boundaries. It is important to point out that the prewarping does not change the nonlinear properties of the transformation, but compensates for it. So a distortion of the phase response of the filter will always be present. One of the most interesting properties of the bilinear transform is the absence of aliasing. The theoretical explanation is very easily understood. As the entire jw axis is mapped onto the unit circle, the folding frequency (z=-1) corresponds to w=-, and the attenuation of any realizable network shall also be infinitely large. However, from the Z domain point of view the reasoning does not seem explanatory since the Z axis is the same whether the filter is mapped using a rational or a nonrational transformation. In the particular case of a filter, the responses must "peak" up at the same place in frequency to make any designs useful, and the response is repeated indefinitely at fs intervals (no matter what jw frequency corresponds to fs/2). The answer must be found in the elements that constitute the network in the Z domain and which are dependent upon the transformation used. In the case of the bilinear, there is always at least one zero mapped at z=-l, which sets at this point the amplitude of the frequency response at zero. So the zero is responsible for the good attenuation properties of the bilinear transform near the folding frequency. This observation is worth noting since it will be useful in the design using nonrational mappings. The best known nonrational transformation is the impulsive invariant (Oppenheim & Shafer, 1975). There are a few slight modifications, and the direct substitution method will be employed here (Childers & Durling, 1975), where (s-a.) + (z-e-iT). (7) ai are the S plane filter poles. Due to the periodicity of the exponential function, the transformation is not one to one. The net effect is that each strip of width ws in the S plane will be overlayed onto the unit circle; i.e., the points wl+Kws (K=1,2,.) will all be mapped into the same point as wl. Mathematically, the frequency response of the digital filter H(ejw) is related to the frequency response of the analog filter (Oppenheim & Shafer, 1975) as H(ejw) = H (j +j i k) (8) T a T T Ha(s) is the Laplace transform of the filter impulse response and T the sampling interval. It is clear that only if H(w) is bandlimited to ws/2, (8) will be a reasonable description. This fact excludes the use of the direct substitution to design high pass filters and in general any function which will display appreciable energy at half the sampling frequency. Anyway it can be expected that the attenuation characteristics of the high frequency end of the filters designed with the straight direct substitution method will be worse than if the bilinear was employed. For the case of bandpass filters, the zeros of the S plane transfer function, if not at s=-, can fold back in the bandpass, ruining the design. Therefore, it is a good design procedure not to map the zeros of the S plane transfer function using this transformation. It was decided to place the zeros independently in the Z domain. As the filter is a bandpass,at least one zero must be placed at z=1 to block the d.c. gain. Taking into consideration the preceding discussion of the effect of the zero at z=-l, it was decided to place the other zero(s) of the transfer function at z=-l to obtain the good attenuation characteristics of the bilinear transform. The other convenient properties of the zeros at z=ï¿½l are the implementation with simple adders (or subtractors) which are fast and do not contribute any roundoff errors. It may be expected that the Z domain filter characteristics will not display constant gain across the passband, since a linear relation was used to map the poles and now a zero is being placed at w=ws/2. To compensate for this fact, one technique places a different number of zeros at z=l, according to the relation (Childers & Durling, 1975) l+cos WcT k = (9) 1-cos WcT where wc is filter center frequency. Bearing in mind the discussion of the lowpass to bandpass transformations, it can be expected that the narrowband will compensate in part the asymmetry since the high frequency pole has already a higher Q than the low frequency one. It turns out that the narrowband transformation used with the direct substitution requires one zero at each of the locations. The number of zeros required for the other combinations (i.e., direct substitution with wideband transformation and bilinear) is shown in Table I for second order filters. From the results obtained with the EEG filters (Principe et al., 1979) it can also be concluded that the sensitivity of the passband gain to different location of the filter center frequency is much lower than (9) predicts. TABLE I REQUIRED NUMBER OF ZEROS Z = 1 Z = -1 Direct substitution 1 1 narrow-band Direct substitution 2 1 wide-band Bilinear 2 2 Another conclusion of the analysis which is very important regards the filter passband gain. In digital filter design the control of the filter midband gain is generally accomplished by properly scaling the filter transfer function. It is suggested that the correct handling of the Z plane filter design may constitute an alternate way to control the midband gain, at least in some cases. An example given in (Principe et al., 1979) is worth mentioning. For the delta filter (0.1-3 Hz, fs=100 Hz), the midband gain was decreased from 878 for the filter designed with the narrowband direct substitution to 87 just by using the wideband transformation (both frequency responses met the specifications). Let us analyze in general the effect of the zeros at z=ï¿½l in the transfer function. For the zero at z=l the numerator takes the form z-l = ejWT11 = 4 sin2 wT (10) and for the zero at z=-l Iz+ll = 4 cos2 wT (11) It is concluded that the effect of the zero at z=-l has the symmetric effect of the zero at z=l, with respect to fs/2, so whatever conclusion is reached for one will apply for the other in the symmetric region. Fig. la shows the attenuation characteristics of EQ(10),(11). Now suppose 2 wT 4 sin2 w 2 r/2 . r wT Fig. la. Frequency Response of One Zero 4 wT 16 cos wT 2 16 sin4 wT Ir/2 7r wT Fig. lb. Frequency Response of Two Zeros that instead of one zero there are two. The numerator becomes (Fig. lb) 2 4 wT z-i = 16 sin F. (12) The region of the unit circle where the two zeros give smaller magnification is 4 wT .2 wT 16 sin 2r < 4 sin w or wT < arc sin 1 - wT < 600. (13) 2 2 It can be said that the wide-band transformation will yield smaller magnification transfer functions when the filter passband falls in the region 00 The design procedure just described has been implemented in a Fortran IV program, presented in Appendix I. The inputs are the frequency domain specifications of center frequency, bandwidth, ripple in-band and out-ofband attenuation rates. The program chooses the filter order to meet the specifications and gives the values of the S plane bandpass filter poles, the z domain filter poles and the quadratic factors C and D for each second order resonator given by (a cascade implementation is assumed) C = 2RE(z) (14) D =I (zi)+I2(z.). The number of zeros at each of the points (z=ï¿½l) is printed. The frequency response of the filter is also plotted for a quick check on the design. The above procedures design infinite impulse response (IIR) filters. For EEG applications, where the attenuation characteristics seem to be more important than the linear phase characteristics, the IIR seem to be preferable to the finite impulse response filters (FIR). However, it is worth mentioning that the FIR can be designed to present a linear phase across the passband, which may be important for applications where the phase information of the input will be an important parameter. Another application which may call for FIR is the situation where multi-sampling rates will be used. One of the great problems in digital filter design, as mentioned earlier, is the lack of independent midband gain control. Therefore, when a narrow-band filter needs to be designed, the gain shall be expected to be high. To decrease it, scaling can be utilized at expenses of dynamic range, or for some cases an alternate transformation can give smaller gains. The other possibility is to decrease the sampling frequency, making the filter appear wider and therefore lowering the Q. Although for certain areas like speech processing this is not of great practical value, for EEG it is a viable procedure, since the EEG rhythms can be .considered not coupled. The penalty paid is higher computation time as a digital lowpass filtering shall be performed before the reduction of the sampling frequency to avoid aliasing. Peled (1976) has a detailed analysis of the procedure including setting of lowpass corner frequency and new sampling frequency. The rationale for using lowpass FIR filters is the following: If decimation is desired (lower sampling frequency at the output), only the samples which will be used need to be computed. In the case of FIR filters, where the output is only a function of the past values of the input, this can be accomplished easily. For IIR filters the present output is a function of past values of the input plus of the output, which means that every point must be computed. The threshold where the FIR filter is preferable to the IIR is given by the ratio of decimation and the difference in filter order to keep similar attenuation characteristics. An automated procedure to design IIR bandpass Chebyshev filters has been presented. It was shown that it is possible to design bandpass filters using the direct substitution to map the filter poles and a technique derived from the bilinear to map the filter zeros. The relations between the design and the implementation are apparent in the choice of integer zeros (to cut computation time and roundoff noise) and in the choice of the transformation, which may help controlling the midband filter magnification. Finite Length Effects of the Implementation After obtaining the filter transfer function in the Z domain, the filter algorithm is the inverse Z transform of H(z). For a second order function of the special type which will be used here, H(z) can be written as in (15) and y(nT) as in (16). (z+l) H(z) = (z(+ ) (15) (z2-Cz+D) Yn = Cyn-l-Dyn-2+Xn-Xn-1 (16) The filtering will then be accomplished by computing (16) in the microcomputer. If an infinite wordlength machine is assumed, the implementation problem is nonexistent. Practically finite wordlengths must be used leading to the fact that almost every digital filter is nonlinear. For this reason the output of the digital filter deviates from what is actually desired. This leads to the study of various ways of computing the filter recursion relation. As was said in the introduction for EEG applications, fourth order analog bandpass filters were found satisfactory. Following the basic rules explained in the previous section the transfer function for the class of EEG filters can be written as H(z) = (-1) (z+l)q (17) (z -C1z+D1)(z -C2z+D2) where l1p, qZ2. There are three basic ways to implement (17). One is called the direct realization and implements the fourth order polynomial directly. The high sensitivity of this realization to the finite length effects of the implementation is well known (Gold & Rader, 1969). The other two are the cascade and the parallel form. For filters with the zeros on the unit circle, the cascade form can be implemented with fewer multipliers (Childers & Durling, 1975). Therefore, it will be the only one studied. The variables which condition the implementation are the type of arithmetic (fixed or floating point), the number representation (sign magnitude, two's complement and one's complement), and the quantization methods (truncation or rounding). For fixed point arithmetic, the implementation is based on the assumption that the location of the binary point is fixed, so numbers must be aligned before performing the additions. If two fixed point numbers are added, overflow can occur (Gold & Rader, 1969), and as there is a linear relation between dynamic range and wordlength, the probability of overflow is large. For the fixed point multiplication, overflow can never occur. However, quantization of the true value shall be generally made, since the product of two b-bits numbers is two b-bits long and it is general practice to use a constant wordlength throughout the filter calculations. For floating point representation a positive number F is represented by F=2CM, where M, the mantissa, is a fraction, and c, the characteristic, can be either positive or negative. Quantization of the mantissa is necessary for both the addition and multiplication. Overflow is also theoretically possible, but as there is an exponential relation between wordlength and dynamic range, it can be essentially excluded. It is therefore assumed that floating point is overflow free. Truncation and rounding are the two forms of quantization. The quantization error is defined by E = x-Q(x) (18) where Q(x) is the portion of x, the input signal, retained. The quantization error depends upon the type of number representation used as shown in Fig. 2. In a microcomputer implementation, the two's complement is the natural choice because it is employed in the processor's arithmetic unit. ROUNDING Q (x) x PR(e) /A TRUNCATION Q (x) PR(e) 1/A -b A=2-b , where b+l is register length. Fig. 2. Quantization Error I Therefore, the following analysis will be restricted to two's complement truncation and rounding. For rounding the error is bounded by (Oppenheim & -Shafer, 1975) -2-b- )-b while for truncation it is -2-b+l (20) <2-b+lx<0 where b+l is the number of bits of the wordlength. The truncation error is bigger than the roundoff error, but filters realized with truncation need fewer operations, as truncation is readily done in the microcomputer. Filter Internal Magnification Fixed point arithmetic must be chosen for the implementation of digital filters in today's microcomputers due to computation time constraints. Therefore, one of the serious problems in the implementation phase is the possibility of overflow in the additions. Overflow causes gross errors in the filtered output and even the possibility of sustained oscillations (Ebert et al., 1969). An important consideration in the design is to ensure that when an overflow, produced by the input, occurs, the filter will recover in a short time. When the recovery is short, overflow propagation is limited and the filter is said to have a stable forced response (Claasen et al., 1976). It can be proved that the conditions to guarantee stability of the forced response are equivalent to the zero input stability of the system; i.e., the filter poles must lay in a specific region inside the unit circle. For the special case of second order filters there is a result that can be easily determined and will ensure filters with no overflow propagation. Consider the second order filter (16) with zero input Yn = Cyn-1-DYn-2. (21) The roots of the filter must lay inside the unit circle, which implies that the values C and D must lay inside the triangle shown in Fig. 3b. Ignoring roundoff error and assuming two's complement arithmetic, the output of the filter.subject to overflow will be Yn = f(CYn-1-DYn-2) where f( * ) is a sawtooth waveform (the characteristic of two's complement arithmetic shown in Fig. 3a). It is clear that if I CYn-1-DYn_21<1 (22) then yn will be correctly interpreted. So a sufficient condition to have filters with no limit cycles (nor propagation of overflow) is 60 2b-1 b-1 x Fig. 3a. Overflow Characteristics of Two's Complement Fig. 3b. No Overflow Propagation Region ICI+IDI which is shown in Fig. 3b as the square hatched. This condition is very restrictive for most applications; therefore, most of the practical designs will display the undesirable effect of overflow propagation. Since the overflow is a nonlinearity in the system, it can be expected that various implementations (different ways of calculating the recursion relation) will display different overflow properties. The parameters one can choose to control the overflow are special type of adders, scaling, and exploring topologic differences in the structures. When the result of an addition is bigger than the register length in two's complement representation, there is an abrupt change in sign. It has been shown (Fettweis & Meerkitter, 1972) that at least for the class of wave digital filters, saturation arithmetic allows for the absence of parasitic oscillations. Here saturation arithmetic means the type that holds the signal at maximum level when an overflow occurs. This implies that overflow must be monitored. The relations between the type of saturation arithmetic, distortions at the output and stability are a present area of research (Claasen et al., 1976). No general results are available. Scaling the signal level is the most general way to handle overflow. To prevent adder overflow the input, xn must be scaled so that the output, Yn<1, for all n. There are a number of approaches that can be considered to accomplish this (Oppenheim &'.Shafer, 1975). As Yn = hkxn-k (24) k=O where hk is the filter response, it is easy to see that lynl lIhkllxn-kl. Z (maxlxkl) I hkl. (25) k=O k=O Therefore, by making maxlxk = 1/ [ lhkl, overflow of yn can k=O be prevented. This means on the other hand that the zero forming paths of the filter must be multiplied by k = 1/ _ lhk . Experiments have shown that this scaling k=0 policy is too pessimistic and does not allow good use of the full dynamic range of the wordlength. A second approach would be to assume the input sinusoidal, with a frequency equal to the resonant frequency(es) of the filter. This is equivalent to requiring (H( . ) is the filter's frequency response) IH(eJX)I 2 1 - 1 n (26) and is accomplished by multiplying the coefficients by k, the scaling constant, so that a0+ale +a2e max k 1 (27) -jx -2jX l+ble +b2e where a0, al, a2, bl, b2 are the coefficients of the second order filter. This scaling policy is sometimes too optimistic, that is to say, it may lead to overflow. A third approach uses the frequency domain relation Y(z)=H(z)X(z). The energy of the output signal is k yk, which by k=O Parsevals relation means m 2 12x j j Yk - 2 70 Y(e )y(e- )dX. (28) k=0 Using Schwartz inequality 2rr 21r. 2 27r 12 j x k Yk 0 (f H(e )H(e )dl(O X(e3)X(e )dX) k=0 2 2 rr -jx S xk0 H(ej )H(e )dX. (29) k=0 Therefore, scaling constants K may be chosen satisfying 0 k2 H(e j)H(e )dX z 1. (30) Again, in some applications this policy would be too optimistic (Gold & Rader, 1969). There is available a general theory, which enables signal modelling by constraining the input norms (Jackson, 1970). For the EEG analysis a 1/f type of spectrum may be an adequate model. However, this study was not pursued due to two main reasons. First, the modelling is subject to some gross approximations, mainly for abnormal EEG. Second, no matter what methodology is used, scaling involves the basic approach of increasing the number of multipliers. In our case, where implementation of bandpass transfer functions were obtained without multipliers in the zero forming paths, scaling means increasing tremendously the computation time, as will be shown quantitatively later. It was decided to allow for a controlled increase in the computation wordlength instead of scaling the filter transfer function. Hence, the internal magnification analysis serves only as an indication of how many bits shall be allocated above the input wordlength to avoid overflow (each power of two corresponds to one bit). Sinusoidal analysis was employed to obtain the information of the wordlength increase, due to its simplicity. It consisted of calculating, using a computer simulation, the maximum value of the transfer functions from the input to the output of each node in the structure. It has been shown that scaling.uses most efficiently the dynamic range of the wordlength for a particular structure (Jackson, 1970). However, scaling is essentially a method of controlling the internal magnification. If a certain .filter needs to be built, in a specific machine and with particular structure, scaling just enables its implementation with a good use of dynamic range. Sometimes this methodology may lead to impractical filters, mainly when microcomputers are used as the implementation medium (DeWar et al., 1978). This motivated the search for structures which would display reduced magnifications. Let us take as an example the canonical structure shown in Fig. 4a and also the structure called feedforward of Fig. 4b. They are the only ones which will be compared in this work. Their big advantage from the microcomputer implementation point of view is the minimum (canonic) number of multipliers. The internal magnification in node 1, G1, is the one which will constrain the choice of the wordlength. As is also known, the zeros may counteract the magnification effect of the poles. In the canonical structure, as the feedback signals are always taken before the zeros are implemented, G1 is independent of the zero placement and number. On the other hand, in the feedforward structure, one of the zero forming paths is introduced before the feedback signals are taken, allowing for decreased gains when the filter poles are close to the zeros, i.e., in our application when the filters are low frequency (zero at z=l) or close to fs/2 Fig. 4a. Canonic Structure Fig. 4b. Feedforward Structure (zero at z=-1). This fact has been experimentally noticed in (Mick, 1975). It may also give a."rule of thumb" for the pairing of zeros and poles in the feedforward structure: the zero closer to filter passband shall be implemented first for reduced magnification. The goal is to bound the maximum magnification by the input-output filter gain. When this is possible, both structures present the same wordlength requirements. However, they may display quite different noise power at the output. Fig's. 5a and b present the structures for fourth order filters designed with the direct substitution, narrow-band, wide-band, and the feedforward implementation. The implementation with the canonic form is pictured in Fig. 5c. Noise Analysis The effect of truncation (rounding) is apparent in three main areas: the A/D conversion, the multiplications, and the coefficient quantization. As the phenomenon has the same basic characteristics, its modelling will be overviewed first. Quantization is a nonlinear operation on the signals, but can be modelled, as a superposition on the signal of a set of noise samples shown in Fig. 6 (Rabiner & Gold, 1975). The quantized samples will be expressed by x(n) = Q[x(n)] = x(n)+e(n) (31) where x(n) is the exact sample and e(n) the quantization e G1 G2 .5 Fig. 5a. 4th Order Filter. Feedforward Implementaton. a e a G 3 i e3 4I Fig. 5b. 4th Order Filter. Canonic Implementation. G e 4th Order Filter. Feedforward Implementation (Wide Design). e a Fig. 5c. e G5 Fig. 5d. 4th Order Filter. Canonic Form (1D) with Scaling. xn=x (t) x(n) =x(n) +e (n) e(n) Fig. 6. Additive Model for Quantization error. It is common to make the following assumptions on the additive noise process (Gold & Rader, 1969): 1. The sequence of error samples e(n) is a sample sequence of a stationary random process. 2. The error sequence is uncorrelated with the sequence of exact samples x(n). 3. The random variables of the error process are uncorrelated; i.e., the error is a white noise process. 4. The probability distribution of the error process is uniform over the range of quantization error. It is worth pointing out that the above-mentioned hypotheses are approximations that do fail for special types of signals, the most important being the class of constant inputs. But when the sequence x(n) is a complex signal and the quantization step is small, it gives good results. Therefore, the model presented in Fig. 6 will be used for most of the quantization effects in the filter implementation. The assumption that the error is signal independent is a good approximation for rounding, but for truncation it is clearly not true, since the sign of the error is always opposite to the sign of the signal (20). However, it is easily shown (Oppenheim & Shafer, 1975) that this effect can be incorporated in the mean of the noise process. For truncation, the mean value me and the variance 2 are given bye given by -b -2 me (32) -2b a2 = 2 e 12 and for rounding by me =0 (33) -2b 2 2 e 12 The probability density function for the error process is shown in Fig. 2. The autocovariance sequence of the error is assumed to be See(n) = ao6(n) (34) for both rounding and truncation. The ratio of signal power to noise power is a useful measure of the relative strengths of the signal and the noise. For rounding or truncation the signal to noise ratio is a2 a2 2b 2 s = (12.2 )a (35) 2 -2b s ae 2 /12 e or when expressed with a logarithmic scale 2 SNR = 10 log o ) = 6.02b+10.29-10 log10 (a). (36) e When the signals are properly scaled to avoid clipping, (36) can be written as (Oppenheim & Shafer, 1975) SNR = 6b-1.24 dB : 6b (37) which shows clearly the interrelationship between dynamic range and quantization error in any signal processing algorithm. For each bit added to the wordlength, the signal to noise ratio improves approximately 6 dB. It was stated earlier that the class of constant inputs needs a special modelling. One case that deserves studying is the zero input condition, since it may arise often in real applications. Theoretically, if x(n) is zero in (16), y(n) will then converge to zero if the filter is assumed stable. It has been observed, however, that in the actual system the output does not always display the theoretical behavior. All components of the filter can only attain a finite number of values since they are quantized. This makes the real filter a finite state machine. An immediate consequence of this is that if the filter output does not converge to zero with zero input, it must become periodic after some finite time. The periodic oscillation is referred to as limit cycle or zero input limit cycle (Rabiner, 1972) and is the prototype of correlation between noise samples. The assumptions underlying Fig. 6 are then not applicable, and the best method for investigating zeroinput stability is the second method of Lyapunov (Liu, 1971). The general analysis is very difficult to carry out, and here only the results will be presented for the case of two's complement truncation and the canonical second order structure. Parker and Hess (1971) showed that the occurrence of limit cycles can occur in the horizontally hatched region of the triangle shown in Fig. 7. The cross hatched region is the stability region (Claasen, et al., 1976). Simulations on digital computers have shown that no limit cycles occur in the remaining region inside the triangle. Probably more important is to have an idea on the bound of the amplitude of the limit cycle. The absolute bound on the amplitude of a limit cycle of any frequency is given by , 1 1-cl+lIDI y 1 n max 4 (1- ) 2 /D (l-D) I-C-/4D D 1 D>O and D(2//D-1))JICI~2/D 1 D>O and ICI.2D(2/D-1)2 (38) for complex roots (Peled, 1976). This bound is too pessimistic, and so the RMS bound is usually preferred (Sandberg &' Kaiser,: 1972). It is given by 1 1-ICI+D -1 ((- y' 2 (1-D)1-C2/4D DZO or D>O and IC >4D2/(1+D) D>0 and IC1<4D2/(l+D) (39) Fig. 7. Stability Region for Second Order Filter o o o 0 0 0 tI 1 f I I I 1 1 I 10 12 14 16 18 Wordlength Fig. 8. Quantization Error Versus Wordlength Cascade Form (after Rabiner and Gold) where p is the period of the oscillation. A few points about limit cycles are worth noting. Filters implemented with rounding possess much smaller stability regions than when implemented with truncation. Filters using floating point arithmetic may possess large amplitude limit cycles. In fixed point arithmetic the limit cycle amplitude is generally small (within the noise floor of the filter). It is a general rule that when the amplitude is high, the source of the oscillation is caused by overflow nonlinearities (Claasen et al., 1976). Each implementation of the same second order structure will display its own stability region, the one presented being just a guideline. As the study of zero input limit cycles is very easy to accomplish with the real filter, further discussion will be postponed to the implementation phase. Coefficient quantization. Usually the coefficients of a digital filter are obtained by a design procedure which assumes very large precision. For practical realizations the coefficients must be quantized to a fixed, small, number of bits. As a consequence, the frequency response of the actual filter deviates from the ideal frequency response. Actually this effect is similar to the effect of tolerances in the analog components. There are two general approaches to the analysis. The first treats coefficient quantization errors as intrinsically statistical quantities, and the effects of the quantization can be modelled as a stray transfer function in parallel with the corresponding filter (Rabiner & Gold, 1975). The second approach is to study each individual filter separately and optimize the pole locations by comparing the transfer functions of the filter realized with high coefficient precision and the one obtained with the actual values. Rabiner and Gold present experimental results for the error as a function of the wordlength. The table for cascade realization is shown in Fig. 8. For filters realized as a cascade of second order structures, the normalized error variance is below 10-2 for wordlengths of 16 bits. This is considered a very good practical agreement. The degradation is very steep towards shorter wordlengths, the variance being 8 to 10 bit coefficients. Therefore, 16 bit coefficients are thought necessary for practical designs. The extrapolation of their results for any filter realized as a cascade of second order sections is not adequate, since the grid of possible pole locations inside the unit circle is not constant. The grid is defined by the intersections of concentric circles corresponding to the quantization of r2 and straight lines corresponding to quantization of r cos 6 (standard complex notation is used). The density of pole locations is less in regions close to the real Z axis as can be shown in the following way. Consider a filter with pole pair (L,K) where (40) -1 k w = wT = cos1 k R 2/; Assuming the errors small Ar = r AL + Dr Ak (41) S- 2 AL + D Ak DL ak and substituting (40) yields Ar 1 AL 2r (42) AL _ Ak 2r2tan 8 2r sin 8 Since AO=tAwR, the error sensitivity is directly proportional to the sampling rate. Furthermore, the second equation shows that error in angle 6 is greater when 8 is small (poles with small imaginary part). This point shall be carefully noted since there is a tendency to believe that higher sampling frequencies merely cause the digital filter to behave more like the corresponding analog filter. This is true only if infinite precision in the computations is assumed. For digital filters, it is therefore convenient to define the wide-band concept in the Z domain, since a wide-band analog filter can correspond to a narrow-band digital filter, if the sampling frequency fs utilized is higher than the Nyquist rate. Also, the implementation of bandpass filters with low center frequencies are the ones which are most sensitive to quantization effects and so should be the ones closely monitored. Practically this can be accomplished by checking the frequency response of the filters obtained with the truncated coefficient values to see if they meet the.desigfn specifications. A/D conversion noise. The A/D conversion noise is the uncertainty derived from the use of finite precision to represent the analog signals. The assumption of independent signal and noise allows one to proceed with the noise computations while ignoring the signal. Let h(nT) be the system impulse response and H(Z) the corresponding transfer function. Then the output autocorrelation function R (nT) for the noise e(nT) is R(k) = y(n-k)y(n) n=O = [ h(p)e(n-k-p) 1 h(m)e(n-m) n=O p=O m=O = ( h()h(m)Rx(k-m-p). - (43) pm The variance is obtained for k=0; therefore a = R (0) = [ h(p)h(m)Rx(m-p). (44) y y p M Using the property of the noise autocovariance function expressed by (34) a2 = 02 1 h2(m) (45) y em=0 where a2 is given by (33). Using Parseval's relation (44) e becomes a2 =2 j h2 (m) = ea2 H(z)H(z-l)z-ldz (46) Y em=0e 2j c where c is taken as the unit circle. This expression is often easier to compute than the infinite summation and will be extensively used. The transfer function for the signal and A/D conversion noise are the same. Since bandpass filters will be implemented, and the frequency spectrum of the input noise is white, the bandwidth of the output noise is narrower at the output, yielding effectively a2 A/D conversion noise is not the only source of noise introduced in the computations. As the goal is to determine the signal to noise ratio at the output, the total noise at the output must be evaluated including the contribution of the A/D noise. The result of (46) implies that the system transfer function is realized without any error. In practice this is not the case, but assuming the A/D conversion noise uncorrelated with the other noise sources, the total output noise can be evaluated by simple addition. Roundoff noise in the multiplications. It is general practice to quantize the output of the multipliers to keep the same wordlength throughout the filter calculations. Essentially the multiplication roundoff noise is the same process as the A/D noise. However, the placement of the noise sources depends upon the structure used and how the arithmetic is actually performed. This means that, in general, the transfer function for the multiplication roundoff noise does not coincide with the systems transfer function. Fig's 5a and 5c show the placement of the noise sources (el, e2, e3, e4) for the feedforward and canonical structure. If the addition was performed in double precision, there would be only one noise source at the output of the adders, but the computation time would have been increased. Truncation of the result is accomplished by just dropping the 16 least significant bits of the product. To compute equation (46) the transfer function from the noise sources to the output must be determined, along with the weight a2 of the noise sources. When the compuin tation wordlength is the same as the input number of bits, there is no question what the weight shall be. However, in this case the computation wordlength is bigger than the input number of bits; hence an analysis was undertaken to see which is the limiting factor. Actually this was a very misleading point in the noise analysis, and so two explanations based in different arguments will be presented: From the point of view of the fixed point arithmetic, before an addition is performed, the addends must be aligned. If the filter computations are correct, the part of the product retained must be aligned with the input; otherwise, the additions which involve the input will not yield the correct result. Hence, the least significant bit of the part of the produce retained has the same weight of the input, even when the computation wordlength is bigger than the input number of bits. The alternate argument follows the "binary point" across the various shifting operations to yield the correct alignment. This will be explained in Appendix II. To evaluate the total noise at the filter output, the roundoff noise power produced by each source will be added to the input quantization noise, since it is assumed that the noise sources are uncorrelated. There are three quantities that are important noise parameters: the peak inband noise referred to 1 volt, the A/D converter precision, and a dimensionless number a 2 /l , which is a measure of out in the noise gain of the system and may be translated into the number of bits of the output that are "noisy" (one bit cor2 2 responds to powers of a /a = 3.98). From this ratio out in and the square of the filter midband gain G the noise factor K G2 K 2 2 (47) 2 2in out in can be defined and enables direct comparisons between structures and designs, since it is independent of the input noise variance. The value KZ1 sets the boundary for the case where the noise created by the digitalization and the calculations will degrade, in a r.m.s. sense, the input (analog) signal to noise ratio. The output signal to noise ratio (peak signal power to r.m.s. noise) can be obtained from (47) through a multiplication by A2la/ where in A is the input signal amplitude. This definition of S/N ratio seems appropriate for EEG signal processing since the various signals of interest are narrow-band and fall in the filter passband. The filters being of order four, an exhaustive search to find the best cascading of sections and best pole zero pairing was possible. The Dinap program (Bass et al., 1978) was extensively applied to access the noise (and magnification) properties of different structures. Fig. 9 shows a typical output. In (Balakrishnan, 1979) a detailed explanation of the program and its use is presented. Practical Considerations Processing Speed In the internal magnification section it was said that the filters were going to be implemented with a computation wordlength bigger than the input number of bits. This is not the general way of filter implementation, which scales e MAXIMUM MAGNIF:'ICAI ION= 0. 60016E I02 FRE:UtENCY ANt.ALYSIS CHEIDYSIGEV 1rTV ORDER III'ULS.-E IN'J. IDAnIDPASS C2, D2 -I-1; CI, MI, --1 ,Lil-J WAVF FREQUENCY (H) lAti IrUDE(D[I) MAON) t lr. MAX TMUM MiiAGN IFICArTIO NM O. 22S19DE 0; FRE(I I)IrICY AblALYS1S CHEYSHIEV 4TIH 1 OR[F.R IMPUL.!SE: I.I' BANIPASS C2.D2, t1sC, 1,DI.--1 SL.A-OW HAVE FRFGUENCY (H) MACN ITUDE ( IB) MA(Il TUDEI MAXIMUI'I HA0C-IIA C OI= O 230093E 02 FREUEtINCY ANAI.YSI CHEDYSIIEV -11H ORDER IIPULSE INV. DANDPASS Cp, D;, +1 CI. D . I lUit-I , WAVE FREGUENCY (H) MA I TUDE ' D ) MlAGtl I TiUDF MAXIMUM MlAGNIFICATI ON- O 479117F 02 Ni)DE NN!tIRER 4 PHIASIE (D) NOIIVE NUIIDER 7 PHA.E (D) MODE NUMBER 11 PHIAGE (D) Fig. 9a. Filter Magnification Analysis Using DINAP e e REAL PART REAL PART RIEAL PART IMAG. PART IMAO. PART IMAG. PART e NI ISE ,ANALY S I L . ,l t i: i t .j 'AGN ITUDn .gï¿½ -59. 86543 .ï¿½ --64. 66954. -69. 47366 . -74. 27779. -79. 08191. . -33. O1603. . -8. 69016. --93. 4942. --9b 29840. -103. 10251. -10/. 90663. O. 10E-'01 l.L . ,i Fi i ) C, jiPUL, fV. .': 0 . 1 ( 1 . . . . . . . I . . . . ' . . if * li M~ w it 9 ï¿½ R -Si I San n m il IF* . *15* 1* * * * * *t Sl * *9 ********* O) 79E 0 O. 6E 02 o0. 2I4E 02 O. 3E 02 O 12L 02 . E 02 o. c8E 02 o. 40E 02 0. 36E 02 FREQUENCY (HTZ) W r ((:2/W*iH'rINIERAL) = O. 30E-03 VOLTS FREIU ENCY ANII YSI NJDE NU 1li-IER 16 Fig. 9b. Filter Noise Analysis Using DINAP e FREGUENCY ANALYSIS NODE IUN IhiLH MAGNITUDE 60. 01599 . 54. 01489. 48. 01379. 42. 01270. 36. 01160. 30. 01050. 24. 00940. 18. 00830. * ï¿½, I ilï¿½ k 16 CIiHE YVSilII.V 41H OIRï¿½L INPUL.LE INV. 6ANDPA'Li (D. Dl?. +l1 1, Di 1 ].0t WAVE It 12. 00720. 6. 00610 .Ci 0. 00500. o. 1OE-1 0. 40E 0i 10" 1Iti i S i k L i i i il il .k:k* * h t & k . . . . . . . . . . . . . . . . . . . *t**k****ht*it**t **tta************* .Qt~RQ*R ~Q+RQtR~Q 0 7/E . i 0. 12~ 02 0. J 6E 02 0. 20E 02 O. 24E 02 0. 32E 02 0. 40E 02 o0 2i 02 0. 36E 02 FREGUENCY (I-ITZ) Fig. 9c. Filter Frequency Response (Slow Wave) the signals within the filter structure. It has been shown (Jackson, 1970) that scaling uses more efficiently the computation wordlength. The values for the scaling factors Gl, G2, and G3 are obtained from the internal magnification analysis. In the scaling design the zeros can no longer be realized as simple adders (subtracters). To use the full dynamic range of the processor the zeros shall be realized as multipliers, increasing the required number from 4 (Fig. 5c) to 9 (Fig. 5d). In a microcomputer implementation this increase slows down drastically the computation, because the multiplication is by far the most time consuming operation. To save computation time, the multiplications can be approximated by shift operations (Oppenheim & Shafer, 1975), but it defeats the purpose of the design procedure. Even in this case, the network complexity increases appreciably (5 more operations per cycle). It can be concluded that scaling better uses the dynamic range at expenses of longer computation times. In (Principe et al., 1979) the results were presented for the class of EEG filters which show that with the nonscaling design only one extra bit is needed to obtain the same output signal to noise ratio. Therefore, for the microcomputer implementation the increase in computation time was thought more demanding than the nonoptimum use of dynamic range. The nonscaling design coupled with a careful filter design enables the implementation of the class of EEG filters with only 4 multiplications per cycle, allowing sampling frequencies up to 3 KHZ. Signal to Noise Ratio Specification Now that the output noise can be calculated and the filter magnification is known, the important step of choosing the hardware to obtain a given output signal to noise ratio can be addressed. From the specification involving the output signal to noise ratio and the degradation produced in the filtering, the input signal to noise ratio can be calculated, which sets the input wordlength. If the analog dynamic range is smaller, an A/D converter that spans the analog dynamic range shall be chosen and a shift left operation up to the requirement of the input wordlength shall be performed. From the internal magnification analysis the increase in wordlength can be determined and added to the input wordlength to set the computation number of bits. In the present form a 12 bit A/D converter and a 16 bit microcomputer are being utilized, allowing for the implementation of filters with gain of 16 or less. There may be cases where a longer computation register length may be desirable, but on the other hand a 12 bit A/D converter may be excessive for EEG processing, mainly when the data are prerecorded on analog tape. The dynamic range of a Sangamo 3500 is 46 dB at 60"/sec. (manufacturer specifications), which suggests that an 8 bit A/D is sufficient. However, the usable signal level really depends on the detection scheme used. For example, a sine wave of 30 Hz, -46 dB below maximum signal, was recognized with a zero crossing detector at the tape slowest speed (15/16"). It also depends on which digital parameters the algorithm uses. In speech processing figures on 20 dB between the smallest signal that is going to be processed and the A/D noise are frequently used (Gold & Rader, 1969). The deterioration of the input signal to noise ratio due to the filtering must also be taken into consideration. It was shown (Principe et al., 1979) that for the case of sleep EEG filters, implemented with the nonscaling design, the noise factor of the implementation ranged from 0.45 to 2.14, i.e., yielding in the worst case a degradation of 1 bit. For all these reasons an 8 bit A/D converter is not recommended to represent the EEG data. A 10 bit A/D (dynamic range of 54 dB) is probably a good choice for EEG applications. If the input of 12 bits obtained with the present A/D converter is scaled to 10 bits, filters with internal gains up to 64 could be directly implemented in a 16 bit microcomputer. Structures which will introduce the lowest noise and at the same time present the smallest magnification (i.e., for which the input-output gain is the limiting factor) are desired. These specifications are in fact contradictory, and the solution is generally a trade-off between increase in wordlength and large signal to noise ratio. There is no present theory to arrive at the general rule to pair the poles and zeros and to cascade the second order sections to obtain the "best" solution. The following rule is suggested to set the pairing of the poles and zeros and the order of the cascade. From the general analysis (every possible combination of poles and zeros and ordering) the first thing is to check if the greatest magnification is implementable in the choice of the microcomputer and A/D converter. If it is, then just choose the ordering that 2 displays the lowest a ut If. the. maximum increase in . wordlength required is too big, then choose the pairing which is accommodated by the combination and then pick up the ordering which has the lower noise. As the general analysis may not be available, it is thought important to extend the discussion (although heuristically) of the pairing of the poles and zeros and section ordering. For the case of fourth order filters the rules are simple: If the higher Q pole is realized first, then the output signal to noise ratio can be expected higher. To realize the highest possible signal to noise ratio, the higher Q pole shall be paired with the zero further away from its resonant frequency. However, this methodology can lead to very large magnifications (or the possibility of overflow). To obtain a compromise between the magnification and the signal to noise ratio, the higher Q pole shall be paired with the zero closer to its resonant frequency and realized first. There is no advantage in realizing the higher Q pole last, |

Full Text |

28
soft (e.g., nonunique) criteria have to be created and the result of the detection always compared to the EEG scoring by the experts. The attractiveness of the procedure resides in the small number of assumptions needed on the data and the simple and fast implementation of most of the parameters. On the other hand, there is no theoretical analysis backing it; therefore, the methodology almost exclusively used is trial and error. The process is long and sensitive to errors in the choice of the key parameters or misjudgements of the criteria. From the time domain techniques known there was some hesitation between employing the rather well-established techniques developed by Smith et al. (1975) to analyze sleep EEG and the ones proposed by Gotman and Gloor (1976) of implementing definitions in the raw EEG. Since the im plementation of a completely digital petit mal detector around a microcomputer was, per se, a task which involved quite a few unknown steps, it was decided to use the tech nique which was best documented. However, it is expected that this work will give quantitative bases for comparing the two different methodologies, with respect to the appro priateness of using prefiltering in the detection schemes of short nonstationarities, like spikes. Chapter II will be devoted to the design of digital filters using microcomputers and study in detail the hard ware selection of A/D converter number of bits and 305 TITL *MA I NF * *:*Â£ * -ait-^e-aSt-^c * MA I NF RE G ***-*.#**.*:** *#****# **#:*:***** * MODULE WHICH PERFORMS A FREQ. * ANALYSIS IN A SPECIFIC FREQ. * RANGE =SÂ¡ REQUIRES PEAK PROGRAM AMD F ILT 4. T PC = 2 1 5A,WSP=FE 40) * fe IT NEEDS A SPECIAL INHHANDLER- -T* * OUTPUT IS A STAIRCASE WAVE WHI IS PROPCP.T IONAL TO DATA FREQ ( * IT ALSO CALCULATES THE RATE OF * OF FREQ. (D/A III. * HRATE- HIGH LIMIT FOR PERIOD 'fe LRATE- LOW ' * PERIOD STORES PERIOD * DETEC FLAG FOR DETEC. -fe PWAVE OLD PERIOD 23C 0 V AORG >2300 foso MA I N EQU >FQ8 0 l EF2 AOC EQU >1EF2 FC06 MASK EQU > FCO 6 F C 08 DE TEC EQU >FC08 FC0 A FERIO EQU >FCQ A 00 15 HRAT E E QU >15 0 0 0F LRATE EQU >F 23 Â£0 PEAK EQU >23E0 FCOC P WAVE EQU > FCO C 2300 02E0 L WPI MAIM 2302 FD8 3 2304 02 05 LI 5 >7FFF 2306 7FFF 230 8 C305 MOV 5,8MASK 230 A FC 06 230 C 0300 LIMI 3 230E 0003 23 10 0207 LI 7 >FE40 2312 FE40 2214 02 03 LI 9 > 21 5 A 2316 21 5A 2318 0407 8LWP 7 GET DATA 231 A C 0 02 MOV 2.0 231 C 4020 SZC 2MASK > 0 231 E F C 06 2320 02EO L WPI >FF 3 A 2322 FF8A PROGRAM :h amp. >/ A I ) CHANG m* ro to to m ro to w to to to to r\j to ra to ro to to to m to to ru to tu to to i\> iv ro ra to ro to to w to to ro to to io to t\> ra to f\) io to to n*O' O' o o> o> it.O' o> at 0101 tJi ui o<.(>NnmA>in'^foomn>ijDO\f>Mo(iin>mO'Olyomn>aiff'f-'iaomm>oso'4>toonm>fflO')> 288 0A43 SLA 3.4 3D03 D I V 3,4 Cl 81 MOV 1 ,6 C 841 SR A 1 ,4 COSI MOV 1,2 3381 .MPY 1 ,2 GAC2 SLA 2,12 0943 SRL 3,4 E0C2 SOC 2 ,3 6103 S 3. 4 C046 MOV 6,1 0200 LI 0,>C 000c CACO SLA 0 ,>c E0 40 SOC 0,1 0380 RTWP Cl 84 SMALL MOV 4 ,6 04C0 CLR 0 4120 SZC 2MASK4,4 F C 26 1616 JNE CO NT 0200 LI 0 ,4 0004 Cl 06 MOV 6,4 4120 SZC 2MASKS ,4 FC28 16 02 JNE SHIFT 0200 LI 0,8 0008 0A05 SHIFT SLA 5,0 04C4 CLR 4 3D03 DIV 3.4 CA 06 SLA 6,0 El 06 SOC 6,4 Cl 81 MOV 1,6 OA 02 SLA 2,0 04 Cl CLR 1 3C43 DIV 3,1 0 A 06 SLA 6,0 E0 46 SOC 6, 1 C08 1 MOV 1 ,2 0902 SRL 2,0 1002 JMP CONTI Cl 06 CO NT MOV 6 ,4 C081 MOV 1 ,2 3381 CO NT 1 MPY l 2 61 03 3 3 ,4 0 ACO SLA 0,>C E040 sac 0,1 0330 RTWP END NO, INC R3 3Y 1 HEX DIG. PERFORM DIVISION ALIGN MEAN VARIANCE IN P.4 OUTPUT TAG RO WILL BE SHIFT COUNT -SHIFT i_. 4 POSSIBLE? NO, REAL MAG. RO HAS S.L. COUNT SHIFT L. 3 POSSIBLE? YES RO HAS S.L. COUNT SHIFT REMAINDER OF DIVISION GET MORE DIGITS IN QUOCIENT FORM EXPANDED WORD EXPANDED 1/N*SUM(X*X) IN P.4 ALIGN MEAN IN SAME WAY MEAN IN P.2 X*X IN R2,3 VARIANCE IN P.4 OUTPUT TAG 261 0110 0046 2160 0206 LI 6,ADC 0111 2162 1EF0 0112 0047 2164 04E6 CLR Â§>6(6) SET GAIN TO 1:1 0113 2166 0006 0114 0048 2168 04E6 CLR Â§>8(6) DISABLE AUTO INC MODE 0115 216A 0008 0116 0049 216C 0205 LI 5, MASK R5 IS MASK FOR CONVERTIO 0117 216E 7FFF 0118 0050 ft 0119 0051 * SET UP OF IMS 9901 AS A TIMER 0120 0052 ft 0121 0053 2170 0200 LI 12,>100 R12 HAS ADDR CF 9901 0122 2172 0100 0123 0054 2174 1E00 SBZ 0 ENABLE INTERRUPT 0124 CC55 2176 1D03 SBO 3 PRIORITY SET TO 3 0125 0056 2178 0300 LIMI 3 SET INT MASK 0126 217A 0003 0127 0057 217C 0200 LI 0,3 9901 FOR IMMEDIATE INT 0128 217E CC03 0129 0058 2180 33CO LOCK 0,15 0130 0059 ft 0131 0060 2182 1CFF JMP $ 0132 0061 2184 1000 NOP WAIT FOR INTERRUPT 0133 0062 2186 4985 3ACK SZC 5,Â§>C(6) END OF CONVERSION? 0134 2188 0000 0135 0063 218A 13FD JEQ BACK NO, JUMP BACK 0136 0064 218C 02E0 LWPI SECND 0137 218E 2100 0138 0065 2190 C29B MOV *11,10 R10,RECEIVES DATA POINT 0133 0066 * 0140 0067 CALCULATION CF 1ST RESONATOR 0141 0063 t THE RECURSION RELATION IS 0142 0069 X1+=-DX2+E+C(X1+E) 0143 0070 ft X2+=X1+E 0144 0071 Y=X2 0145 0072 REGISTER USED 0146 0073 ft REGISTER 10 FOR INPJT Â£ 0147 0074 REGISTER 7 FOR X1 0148 0075 ft REGISTER 3 FOR X2 0149 0076 0150 0077 2192 0743 ABS 3 GET SIGN CF Y 0151 0078 2194 3802 MPY 2,3 -DX2 IN R3 R4 0152 0079 2196 1101 JLT P0ST1 JUMP IF R3 PCS 0153 0080 2198 0503 NEC 3 IF NEG COM PL R3 0154 0031 219A ACCA posh A 10,3 -DX2+E IN R3 0155 0082 219C 0243 MOV 3,9 SAVE IT IN 39 0156 0083 215E A1CA A 10,7 X1+E IN R7 0157 0084 21A0 0007 MOV 7,3 X1+E IN R3 0158 0085 21A2 0747 ABS 7 0159 0086 21A4 3901 MPY 1,7 0160 0087 21A6 1505 JGT PGST2 0161 CHEBY PACE 3 0162 RECD LOG OBJ SOURCE STATEMENT 0163 0088 21A8 0A47 SLA 7,4 ALLIGN BINARY POINT 0164 0089 21AA 0908 SRL 8,12 XFER 4 BITS FRCM LOWER REG 0165 0090 21AC E1C8 SOC 8,7 TO HICHER ORDER REG 0166 0091 21AE 0507 NEG 7 0167 0092 21B0 1003 JMP P0ST5 0168 0093 21B2 0A47 PCST2 SLA 7,4 319 2000 **** * **** * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * TI TL OAT A* *** ***************** ************** *** OAT A ************************************* MODULE THAT FORMATS THE DATA TO BE PLOTTED BY PLOT* IT TAKES CARE CF ALL INITIALIZATION SINCE PLOT IS MADE RELOCATADLE INPUT PARAMETERS FD63 0 IF SZ<3 SEC ARE CONSIDERED NOT ZERO OTHERWISE FD6A- LAST ADDR OF TAG+2 FCO0 DISPLACEMENT TO TAG OF X DATA FC02- '* Y DATA FC18- PLOT TYPE 1 POINTS (FROM X AXIS) 2 HISTOGRAM 3 LINEAR INTERPOLATION HYSTOGP AM ONLY ALLOWS + INC IN X DATA IS ASSUMED BEG. AT 2720. ACCORDING TO FORMAT EXPLAINED IN P.M. SZ DETECTOR MAX X APPEARS MAX Y APPEARS IN FD32 IN FD86 ARE USED WITH TEKT, FOLLOWING CONSTANTS DISPLAY SCOPE FCO 4 (CONS) 300 FC06 (XO) SFF FC08 (YO) SFF FC14 (SCALE.) 8FF FC16 (MAX X.Y) 1400 R12 OF WSP PLOT NEEDS TO BE LOADED BEGINNING ADDR OF PLOT PLUS 39A WITH AORG >2000 0 080 MONI T SOU >8 0 FFFF KEY EQU >FFFF 2720 STATS EQU >2720 FC 12 MASK EQU >FC 1 2 FD68 WDATA EQU >FD68 FD 80 WPLOT E QU >FD8 0 2200 OUT^C EQU >2200 FCOO AODX EQU >FC0 0 FC 02 ADDY EQU >FC 02 FC04 CONST EQU >FC04 FC 06 CLD X EQU >FC0o F C 08 CLDY EQU >Fcoa 611 10 matrix equation involving the autocorrelation function (Yule Walker equation (YW)). To estimate the filter coef ficients, the YW equation can be solved. Then the autocor relation functions for various lags have to be computed and a matrix inversion performed. Algorithms due to Levinson (1947) or Durbin (1960) are generally preferred to solve it recursively. A technique due to Burg (1967) suggests another esti mation of the coefficients that does not require prior estimate of the autocorrelation function. It fits succes sively higher order prediction error operators to the input data by convolving the filter in both forward and backward direction. The square of the two error series are added to obtain the error power,and the filter coefficients are determined by a minimization procedure. A recursion rela tion is also available (Burg, 1967, 1972; Anderson, 1974). Actually the two different techniques (YW method and Burg) use the data in two different ways. The YW assumes that the data is zero beyond the limits of the observation win dow. The Burg technique extends the estimation of the autocovariance coefficients beyond the initial number of points in such a way that the entropy of the autocovariance function is maximized at each step. As can be expected, the two techniques have quite different properties. The estimator by the YW method has the attractive property that its mean square error is generally smaller than other esti mators (Jenkins & Watts, 1969). However, the estimates are 274 False Dectection y--| in- rmr-T-Ti TT~ -7 -tt-t--g--rr- r r nnrrin Fig. 49a. Mean and Variance of Repetition Period Versus Time (#7-1) 06 sec 1.4 sec 150 pv i 1sec. Fig. 27b. PM Detector Output (Patient #12) 155 307 140 The block diagram of the detector is shown in Fig. 24. It is composed of two parallel processing channels, one for the slow wave and another for the spike. In a real time operation, with a single processor, the computations must be performed serially. It turns out that these facts (real time operation and parallel channels) put considerable con straints on the program articulation since for each input sample, decisions in each channel must be made (e.g., did a zero crossing occur?), which implies that the linking among the program modules may be data dependent. The goal in the program design, which used the top down approach, was modularity. Each program module was associated with a particular function of the detection algorithm, and care was taken to make the modules as inde pendent as possible. However, the ideal, one entry point one output (Tausworthe, 1977) would be achieved with much higher program complexity. It is considered that a good compromise between a reasonable increase in processing time and program modularity was achieved. A functional description of each module is presented next. Initialization of A/D converter, timer (TMS 9901) to set the sampling frequency, filter coefficients, special pointers and external references. It also decides about the periodicity of PM complexes (beginning of seizure) and outputs flags. MAIN 253 addr + EH. The board has also available 2, 12 bit DAC channels (base addr + 0 and + 2). The TMS 9901 Programmable System Interface chip' was used as a timer to set sampling frequency. It is addressed through the CRU (addr 10 0), and it uses the II interrupt priority 3. The 9901 needs to be set in the timing mode with the interrupt enable and the processor interrupt mask set at three. When the value loaded in bits 1-12 of the CRU, which is decremented at the clock frequency over 64 reaches zero, the interrupt is activated. In the TM 990/101M-1 board interrupt 3 is vectored to the beginning of memory, where PC is loaded with FFAA. and WSP 11 with FF88 (TMS 9900 Data Manual). H The flowchart of the filter program is shown in Fig. II-2, and it is configured in an endless loop to be used as a real time filter. Minor modifications were introduced in this program flow to be used as a subroutine in the detector (instead of closing the loop, the program ends with a RTWP instruction). Implementation of the Filter Algorithm To accomplish a digital filtering function only addi tions, multiplications and delays (memory) need to be implemented. Fig. II-3 presents, for convenience, one sec tion of the structure used as in the implementation of the EEG filters. After the topology of the resonator, the next value x* and x* at the input of the storing elements (state Ill For this value of t2 the error is (t2*-T) sin wt2+t2 sin w(T-t2) sin wt2+sin w(T-t2) (53) Now substituting in our example of the 80 Hz sampling fre quency and 16 Hz sinusoid the maximum error is 0.33 ms. As expected, the interpolation decreased the error a lot. This point brings up the question of "availability" of information in the digital representation of the signal. The Nyquist theorem only shows that, if the signal informa tion is used from the infinite past to the infinite future, the digital samples possess the total information for wave form reconstruction. In practice, the assumption is impos sible, and therefore a certain error must be tolerated. The practical question is what is the best compromise, sampling frequency/reconstruction method from the point of view of accuracy and computation speed, when both quantities are considered as variables. It was shown that for the simple case of frequency determination by period measure ments, the interpolation method was much more accurate than the simple zero-crossing. However, interpolation is much more time consuming. For this.particular application where the arithmetic facilities are scarce, it seems appropriate to increase the overall digital representation of the wave by sampling the analog waveforms at a higher rate, and use fast, easily implementahle detection schemes like -.zero crossing analysis. The savings are even more apparent when 231 and with similar parameters to set the middle of the period/amplitude windows and then allow a certain varia tion. These methods have to be tested before any definite solution is taken. Both detectors employed, the slow wave and spike, worked well with EEG activity. However, the latter is sen sitive to artifacts. Any sharp transient in the EEG was detected as a spike, even when it was "step like." This fact does not validate the initial assumptions and calls for the inclusion of a more sophisticated pattern recogni tion in the spike channel. Fig. 50 shows the performance of the spike filter (downward pulse in channel II) in a moving artifact epoch. The modification that is proposed is to employ the spike detector presented in Smith (1974), comprising requirements in the positive and negative slopes as well as sharpness measurements. This improvement coupled with the analysis of two frontal channels to use the bilateral synchrony of the PM activity will consider-- ably reduce the false positives. As mentioned earlier, with minor modifications the system can process two chan nels of information in real time. If the proposed modifi cation in the spike channel makes the processing time too long, a simplification on the spike bandpass filter will be evaluated. It will consist of coupling a second order high pass filter with a number of zeros at fs/2 such that enough attenuation is obtained at the higher frequencies. The rationale is to save computation time (the zeros are easy 64 1/f type of spectrum may be an adequate model. However, this study was not pursued due to two main reasons. First, the modelling is subject to some gross approximations, mainly for abnormal EEG. Second, no matter what methodol ogy is used, scaling involves the basic approach of in creasing the number of multipliers. In our case, where implementation of bandpass transfer functions were obtained without multipliers in the zero forming paths, scaling means increasing tremendously the computation time, as will be shown quantitatively later. It was decided to allow for a controlled increase in the computation wordlength instead of scaling the filter transfer function. Hence, the internal magnification analysis serves only as an indication of how many bits shall be allocated above the input wordlength to avoid overflow (each power of two corresponds to one bit). Sinusoidal analysis was employed to obtain the information of the wordlength increase, due to its simplicity. It consisted of calculating, using a computer simulation, the maximum value of the transfer functions from the input to the output of each node in the structure. MAX I MUM MAGNIFICAUOM 0.6001 ACE 02 FREQUENCY ANALYSIS NOCE NUMRER 4 CHEBYBI1CV 4TH ORDER IMPULSE INV. BANDPASS ca. 1)2 I-11 Cl, 1)1, -1 FI .on WAVE FREQUENCY (Hi MAONII'UDE(DU) MAGM] ft IDE PHASE (D) REAL PART IMAG. PART MAX IMUM MAGM1FICAX1 ON 0.SP3195E 02 FREQUENCY ANALYSIS NODE NUMBER 7 CHEBYSHEV 4TH ORDER U1PUI ca. Lia. HÂ¡ Cl. 1)1. -1 SLUM -BE INV BANDPASS WAVE FREQUENCY MAGNITUDE PHASE (D) HEAL PART IMAQ. PART MAXI Ml JIT MAONIF i C AT I Oil- 0 2S00H3E 02 FREQUENCY ANALYSIS NODE NUMBER 11 CHEBYSHEV 11H ORDER IMPULSE INV. BANDPASS Ca. 1)2, +1 Â¡ C1. D1, 1 5I..UW WAVE FREQUENCY (H) MAGNITUDE'DR) MAGNITUDE PHASE (D) REAL PART IMAG. PART MAXIMUM MAGNIFICATION- O 477117F 02 03 cn Fig. 9a. Filter Magnification Analysis Using DINAP 59 said to have a stable forced response (Claasen et al., 1976). It can be proved that the conditions to guarantee stability of the forced response are equivalent to the zero input stability of the system; i.e., the filter poles must lay in a specific region inside the unit circle. For the special case of second order filters there is a result that can be easily determined and will ensure filters with no overflow propagation. Consider the second order filter (16) with zero input (21) Yn = Cyn_!-Dyn_2. The roots of the filter must lay inside the unit circle, which implies that the values C and D must lay inside the triangle shown in Fig. 3b. Ignoring roundoff error and assuming two's complement arithmetic, the output of the filter subject to overflow will be where f( ) is a sawtooth waveform (the characteristic of two's complement arithmetic shown in Fig. 3a). It is clear that if cyn-rDyn-21
(22)then yn will be correctly interpreted. So a sufficient condition to have filters with no limit cycles (nor propagation of overflow) is 294 2654 0280 Cl 13.SZMIN INC RESPECTIVE SZ COUNT 2696 0200 2698 1 503 JGT CO NT 269 A 0 5A0 I NC 3LES33 SZ<3 269 C FC2A 269E 1032 JMP CUT 26A 0 0 280 CO NT Cl 13,TEN 26A2 0960 26A4 1503 JGT aiG 26A6 05A0 INC 2LES10 3 FC2C 2 6A A 10 02 JMP STA 26 AC 0 5 AC BIG INC 2GRT 10 3Z> 1 0 2 6AE FC2E READY FOR TOTALS 2 63 0 05A0 STA INC FLAG 2662 FC22 - 2634 0420 5LWP 3STAT1 R.RATE 2666 FCOE 2 638 C360 * MOV 2ST ATI. 13 2 63 A FCOE 2 6 SC C.FAD MOV 221 13) *14+ MEAN 263E 0002 26C0 CFAO MOV 38(13)*414+ VARIAMCE 26C2 0008 1/2 PERIOD OF S.W. 26C4 0420 BLWP 3STAT2 26C6 FC 12 26C8 C360 MOV 2STAT2,13 2 6CA FC 12 26CC CFAD MOV Â£21 13) 4 14+ 26CE 0002 2600 CFAO MOV 23(13),414+ 2602 0008 2604 0420 OLWP 3STAT3 DELAY SPIKE S.W. 2606 FC 16 2608 C360 MOV 2STAT3,13 2 60 A FC 1 6 2 6DC CFAO MOV 32( 13) ,414+ 26DE 0 0 02 26F. 0 CFAD MOV 38(13),414+ 26E2 0 0 08 S.W. AMPLITUDE 2cE4 0420 BLWP 3STAT4 2 6E6 FC IA 26E8 C360 MOV 2STAT4, 13 2 6E A FC 1A 26EC CFAO MOV 32( 13) ,4 14 + 2 6EE 0002 2 6F 0 CFAD > a y 38(13),414+ 26F2 0008 SPIKE AMPLITUDE 2 6F4 0420 3L WP aSTAT5 2 6F6 FC IE 26F3 C3 60 MOV 2STAT5,13 26FA FC IE 2 6FC CFAD MOV 22( 13) ,414 + 26FE 0002 2 70 C CFAD MOV 38(13),414+ 2702 0008 2704 02 00 CUT LI 13,>FFFF OUTPUT SEPARATION TAG 2706 FFFF 270 8 CF 80 MOV 13,414+ 2 70 A C30E MOV 14.aSZTIM STORE MEM POINTER 270 C FC24 2 70E 04 EO CLP 3FLAG NO TOTALS FOR STATIS 2710 FC22 2712 C80F MOV 15,3THRE S RESTORE S.W. AMP THRES 2714 F C 04 2716 C460 B 2SZCUT GO TO MAIN 2718 2324 END 80 = /L 0 = wT = cos-"*' R 2/L Assuming the errors small Ar If AL + w Ak A0 = M AL + M Ak 3L 3k and substituting (40) yields Ar = _ 1 2r AL A0 = AL Ak 2r tan 0 2r sin 0 (40) (41) (42) Since A0=tAwR, the error sensitivity is directly propor tional to the sampling rate. Furthermore, the second equa tion shows that error in angle 0 is greater when 0 is small (poles with small imaginary part). This point shall be carefully noted since there is a tendency to believe that higher sampling frequencies merely cause the digital filter to behave more like the corresponding analog filter. This is true only if infinite precision in the computations is assumed. For digital filters, it is therefore convenient to define the wide-band concept in the z domain, since a wide-band analog filter can correspond to a narrow-band digital filter, if the sampling frequency fs utilized is 236 Hi ru. - I IA - CALIBRATION WITH SWEEP SINE 1-7 Hz inn; nfi H r*i 'U* V r U L inr o in iH ^7 ,Uwl h uJ^L_J^ 4 sec Fig. 51a. Change of Frequency within Seizure (Patient #7-1) 167 Subject Profile: NAME: last first AGE: SEX: M F RACE CLINICAL HISTORY: START FOOT.: END FOOT.: TAPE *: SUBJECT *: I B K other: MEDICATIONS: ADDITIONSL INFORMATION/OBSERVATIONS: Fig. 30. Technical Log Used in the Data Collection 259 period successively reduced. When the computation time is longer than the time between successive interrupts, the program is interrupted in the middle of the loop. Due to the existence of the wait loop at the end of the computa tions, the width of the time step of the D/A converter increases approximately to the double of the value pre viously used. The value loaded in the 9901 that produced the jump was IB, which corresponds to 277 ys. It is worth noting that the interrupt handler is also included in the processing time. The filter implemented needed only the alignment of one coefficient. One version of the program is shown next, along with the corresponding interrupt handler. REFERENCES Akaike, H. "Power Spectrum Estimation through Auto regressive Model Fitting." Ann. Inst. Statist. Math., 21:407-419, 1969. Akaike, H. "Statistical Predictor Identification." Ann. Inst. Statist. Math., 22:203-217, 1970. Andersen, N. "On the Calculations of Filter Coefficients for Maximum Entrophy Spectral Analysis." Geophysics, 39:69-72, 1974. Anderson, R. B. Proving Programs Correct. Wiley, New York, 1979. Balakrishnan, S. K. "Digital Filters for EEG Processing," Master s Thesis, University of Florida, 1979. Barlett, M. S. An Introduction to Stochastic Processes with Special Reference to Methods and Applications. Cambridge University Press, Cambridge, 1953. Barlow, J. S., and E. N. Sokolov. "Selective On-Line EEG Filtering by Means of a Minicomputer." Electroenceph. Clin. Neurophysiol., 39:208, 1975. Barlow, J. S., and J. Dubinsky. "Some Computer Approaches to Continuous Automatic Clinical EEG Monitoring." In Quantitative Analytic Studies in Epilepsy. Raven Press, New York, pp. 309-327, 1976. Barnes, C., Casper, W., and A. Fam. "Minimum Norm Recur sive Digital Filters That Are Free from Overflow Limit Cycles." IEEE Trans. Circuits and Systems, CAS-24:569- 574, 1977. Bass, S. C., Grundmann, J. W., and S. E. Belter. "DINAP II: A Digital Filter Analysis Program." Purdue University Technical Report, TR-EE 78-14, March 1978. Bickford, R. G., Brimm, J., Berger, L., and H. Aung. "Application of Compressed Spectral Array in Clinical EEG." In Automation of Clinical Electroencephalog raphy Raven Press, New York, pp. 55-64, 1973. 333 614 sec 1.39 sec II I i II I LI Fig. 41c. Variance R. Period Versus Duration (#5-3) 214 77 Fig. 7. Stability Region for Second Order Filter u o u M W T5 H 03 e h O Z Wordlength Fig. 8. Quantization Error Versus Wordlength Cascade Form (after Rabiner and Gold) 218 #12, sometimes with a perfect linear trend (Figs. 44a and b). For patients #5 and 12 the half period of the slow waves was insensitive to the repetition rate. The varia bility on patient #7 did not allow any conclusion (false detections?). The amplitude of the waves in patients #12, 5-1 and 5-2 decreased with longer repetition periods (Fig. 45), and for the other no conclusions could be made. Next, the half period of the waves is going to be used as X axis. In all the patients the variability increased with longer half periods (Fig. 46). The slow wave ampli tude decreases with the half period for patients #12-3 and 5-2. For the others no conclusions could be made. The spike amplitude and the slow wave amplitude were also plotted. Only in patient #7-1 do they seem to be correlated in the sense that high spike amplitudes were coupled with high slow wave amplitudes (Fig. 47). In the others no apparent correlation could be perceived. In summary the trends in the data set are the following: 1) The half periods of the waves remained approximately constant across seizure duration. 2) The variability of the repetition rate increased with increased repetition rates. 3) Longer seizures tend to have smaller repetition rates. 4) The duration of the seizures did not affect the amplitude of slow waves and spikes. 6 points is used to evaluate the periodogram, very rapid fluctuations are present in the estimate of the spectrum. This fact shows that increasing the number of data points of the observation sequence will not change the statistical properties of the estimator since the individual signifi cance of the data points remains the same. To increase the stability of the spectral estimates, the number of indepen dent observations have to be higher or in some way the sta tistics of the observation have to be improved. These two points lead to the two main methods of increasing the sta bility of the estimationthe averaging of periodograms (Barlett, 1953), and the smoothing by means of a window (Welch, 1967). Let us describe briefly the two methods. The Barlett method averages the periodogram over certain frequency bands. To achieve this the data are divided in k segments; the periodogram is calculated for each segment; and the final estimate is the average of the various periodograms. The procedure leads to a bias estimate, but the variance decreases with k (overlapped segments can be used at expenses of lower independence of individual segments and smaller resolution). The number of points on the power spectrum is decreased by k. The Welch method windows segments of the data before the periodogram is calculated. Generally the triangular window is used to ensure positive estimates of the power spectrum. The methodology is the same as for the Barlett 35 in the processor wordlength. The wordlength is also coupled to the precision of the representation because the numbers that represent the constants (and signals) must be quantized to fit the wordlength. For instance, if an 8 bit processor is used, the constants are represented only by two hexadecimal digits. The choice of software implementation of digital fil ters in microcomputers imposes stringent constraints which must be carefully examined to ensure that the analog input signal to noise ratio is not degraded. Nevertheless, the use of microcomputers is thought a good choice due to their cost, size, availability, and the ease with which the fil tering function is actually performed, requiring primarily software knowledge. In this specific application, as rather sophisticated pattern recognition algorithms will be necessary, the microcomputer will also be time shared, to accomplish these functions. Design Criteria The design of digital filters can be broadly divided in two phases. The first is related with the determination of the filter algorithm and the second with the implementa tion of the recursion relation. Generally they are taken independently since the type of problems encountered are quite different. In the design of the filter algorithm a machine of infinite wordlength is assumed. The goal is to arrive at a recursion relation that better fits the BIOGRAPHICAL SKETCH Jose Carlos Santos Carvalho Principe was born April 13, 1950, in Porto, Portugal. He received his bachelor's degree in Electrical Engineering from the University of Porto in 1973. He was awarded a Fulbright-ITT scholarship to pursue graduate studies at the University of Florida, where he received the Master of Science degree in 1975. He returned to Portugal and joined the staff of the Elec trical Engineering Department of the University of Aveiro as an assistant professor. He was awarded a NATO-Portugal scholarship to pursue studies towards a Doctor of Philosophy degree at the University of Florida. 346 334 Birkemeier, W. P., Fontaine, A. B., Celesia, G. G., and K. M. Ma. "Pattern Recognition Techniques for the Detection of Epileptic Transients in the EEG." IEEE Trans. Biomed. Engng., BME-25:213-217, 1978. Blackman, R. B., and J. S. Tukey. The Measurement of Power Spectra. Dover Publishing Company, New York, 1958. Blinchikoff, H. J. Filtering in the Time and Frequency Domain. John Wiley & Sons, Inc., New York, 1976. Bodestein, G., and H. M. Praetorious. "Feature Extraction for the EEG by Adaptative Segmentation." Proceedings IEEE, 65:642-652, 1977. Bourne, J. R. "Computer Quantification of EEG Data Recorded from Renal Patients." Computers and Bio medical Research, 8:461-473, 1975. Box, G. E., and G. H. Jenkins. Time Series Analysis: Forecasting and Control. Holden-Day, San Francisco, . 1970. Buckley, J. K., Saltzberg, B., and R. G. Health. "Decision Criteria and Detection Circuitry for Multi-Channel EEG Correlation." IEEE Region III Convention Record, New Orleans, 1968. Burg, J. P. "Maximum Entrophy Spectral Analysis." 37th Meeting of Soc. of Expo. Geophys., Oklahoma City, 1967. Burg, J. P. "The Relation between Maximum Entrophy Spec tral Analysis and Maximum Likelihood Spectra." Geo physics 37:375-376, 1972. Caille, E. J. "Apport de 1'analyse Spectral de l'EEG pendant le Sommeil." Societe Psychol. France, 7, 1967. Carrie, J. R. G., and J. D. Frost. "A Small Computer System for EEG Wavelength Amplitude Profile Analysis." Int. J. Bio-Medical Computing, 2:251-263, 1971. Carrie, J. R. G. "A Technique for Analyzing Transient EEG Abnormalities." Electroenceph. Clin. Neurophysiol., 32:199-201, 1972. Carrie, J. R. G. "A Hybrid Computer Technique for Detect ing Sharp EEG Transients." Electroenceph. Clin. Neurophysiol., 33:336-338, 1972a. o o o o 0001 0 i.i ME ' FIL. TRO JDB * 1 UOt. > c'iJ68* ? .-Â¡j U.-1 .1 L- PPINL IRE j L'LH''=M 0003 ." PASSWORD 0004 RDUTE PRINT REMDTF6 ' 0005 /v EXEC F0PT6.CS 0006 F DPT. SYS IN PL 000? /'INCLUDE CHEEY1 0008 ./ 0009 'bO . SYSLIE BE DSN =SYS 1.FORTLIE?DISP=SHP 0 01 0 LB DSN =UF. 3062099.EELIB.PISP=SHR 0011 , ED DSN =6RTDR.SSPLIE.FORT,DISP=SHR 0012 'SO. SYSIN PL 0013 4 0 404800 15 12 01 0 014 3 0 3 04:30 0 15 12 01 0015 50 5 04P 0 0 15 12 01 0 016 0 0 0 0 0 0 0 0 0 0000 00 0 0 0 00 0017 0018 FDJ END OF li.lORK FILE Fig. 1-2. Sample Input for Filter Program 246 256 variables) are expressed as a function of the past values x^ and. x2 and' the: input E by x* = (-Dx2+E)+C(x1+E) x2 = Xl+E y *2. As the fourth order filter is implemented as a cascade, the corresponding workspaces were overlapped one register to pass automatically the output of the first resonator to the second. The values of the filter coefficients, C and D, were converted to hexadecimal as fractional numbers and stored in the initialization. Generally they are less than one, but C can be greater than unity, which means that if the coefficients were stored preserving their relative mag nitude, part of them will be only represented as 12 bits. This is thought undesirable due to the quantization effects. Therefore, it was decided to store all the coefficients with full precision, leaving the alignment for the computa tions. For instance, for one of the EEG filters the coef ficients are C1 = 1.3848 D = 0.7256 C2 are = 0.5159 stored as D2 = 0.7256 C, = 1628 D, = B9C7 1 H 1 H = 8417 n = B9C7 . 2 H 2 H 74 m. a 2 _ e -2b 2 12 and for rounding by me = (32) (33) The probability density function for the error process is shown in Fig. 2. The autocovariance sequence of the error is assumed to be Yee(n) =Og6(n) (34) for both rounding and truncation. The ratio of signal power to noise power is a useful measure of the relative strengths of the signal and the noise. For rounding or truncation the signal to noise ratio is (12.2 2b (35) or when expressed with a logarithmic scale 2 o 0 SNR = 10 log, f-|) = 6.02b+10.29-rlO log, rt(cT). (36) 10v 2J 10 e- 299 2114 390 MPY 1 ,7 2116 1503 JGT P0ST7 2118 0A17 SLA 7.1 2 11A 0507 NEG 7 211C 1001 JMP POSTS 211 E 0 417 P0ST7 SLA 7, 1 2120 A 1C9 POSTS A 9,7 2122 C3C3 MOV 3, 1 5 2124 0 82F L> SR A 15.2 * A 2ND RESONATOR 2 126 02E0 LWPI FORTH 212 8 FD3E 212A 0 743 A8S 3 212C 38C2 MPY 2, 3 2 12E 1101 JLT POSTA 2130 0503 NEG 3 2132 AOCO POSTA A 0,3 2134 C243 MOV 3.9 2136 6 ICO S C 7 2138 C0C7 MOV 7,3 213A C747 A3S 7 213C 390 MPY 1 ,7 2 13E 1503 JGT POSTB 2140 OA 17 SLA 7,1 2 142 0507 NEG 7 2 144 1 001 JMP POSTC 2 146 CA17 P0ST3 SLA 7, 1 2 143 A1C9 POST C A 9, 7 2 1 4 A 0 43 MOV 3,5 2 14C 0335 SR A 5,3 2.14E C805 MOV 5,aOUTl 2150 FD 66 2 152 0380 RTYJP END GAIN i:i.3 0* NO.OF ERRORS IN THIS ASSEMBLY^ 0000 OF RELOCATABLE LOCATIONS USED = 0000 270 1 309 TITL I NTHANF* ****3$:## *** ******* fr.#^**:*#*-**#:***:* * INTHANDFREQ ***4 -***.*:*<: #**#*=(: 44 44 *444*4 * MODULE TO 3E USED WITH MA I NF * BESIDES SETTING FS IT ALSO * CALCULATES THE FREQ. AND ITS 4 RATE OF CHANGE IN A PRESCRI3ED * FREG, BAND OF THE INPUT 2200 AORG >2200 FCOA PERI C EQU >FC0 A FC08 DE TEC EQU >FC08 FCOC PWAV E EQU >FC0C 2200 0201 LI 1, > 1EFO 2 20 2 1 EFO 2204 C 342 MOV 2,0>A< 1 ) 2206 COOA 220 8 02 OC LI 12.>100 220 A C 1 00 220C 02 03 LI 3>187 220 E 0187 2210 33C3 LDCF 3,15 2212 IE 00 SDZ 0 2214 1 D03 S30 *3 221 6 C300 LI MI 3 2213 0 0 03 221 A 0 60 A DEC 10 22 1C 1 603 JNE TIME 2 21E 0569 INC 9 2220 0 2 0 A LI 10,>F0 2222 0 OFO 2224 05CE TIME INCT 14 2226 06 20 DEC 3PERI0 MORE POINTS? 222S FCOA 2 22 A 16 1C JNE CONT 2 22 C C142 MOV 2,5 MIDDLE 3AND PEP IOD TO RS 2 22E C320 MOV SDETEC.3DETEC I NSAND WAVE? 223 0 FC08 223 2 FC08 2234 1306 JEQ CONTI 2236 C820 MOV FCOC 223 A FCOA 223C 04E0 CLP. SDET EC CLEAR FLAG 223E FC OS 2 24 0 1004 JMP C0NT2 2242 C802 CONTI MOV 2.DPER 10 PERIOD GETS MIDDLE VALUE 2244 FC OA 2 246 C444 MOV 4, *1 OUTPUT 2248 10 0D JMP CONT 224A 6 160 CO NT 2 S SPERIC,5 GET PERIOD DIFF 224 C FCOA 224E 0 A 1 5 SLA 5,1 AMP. IT 2250 C4 45 MOV 5,41 2252 61 85 S 5,6 SUBT. FRO* PREVIOUS 2254 0746 A3S 6 2 256 39A0 MPY 2PERI0,6 GET RATE GF CHANGE 2258 FCOA 2 2SA 1501 JGT POST 225C 0507 NEG 7 225E Cl 85 POST MOV 5,6 UPDATE OLO PERIOD 2260 C 307 MOV 7 ,3> 1EF2 OUTPUT RATE GF CHANGE 2262 1 EF2 2264 0380 CONT RTWP END 281 318 235 the development of a portable seizure detector, using the proposed design. Studies are under way toward this direc tion. The constancy of the period parameters in the PM paroxysms is intriguing, when seen from the variability of the EEG normal rhythms. Also the abrupt beginning and end of the paroxysm are intriguing. Is there any seizure parameter which can be coupled with the end of seizure? Again one of the candidates is the PM recruiting rate, since it is known that it slows down generally towards the end of seizure. A slight modification in the PM detector is being evaluated at this time to study the frequency variation in narrow EEG frequency bands (programs MAINF and INTHANF). The system's output is a pulse, modulated in amplitude with the period of the waves, which meets the period criteria. Figs. 51a and b show examples of the out put (channel III) for two different seizures. The rate of change of frequency is also outputted (channel IV). Another area which will be interesting to study is the relation of background activity and seizure occurrence. The EEG laboratory at the University of Florida is now capable of performing such automated analysis using the described PM detector and the Sleep Analyzing Hybrid Com puter (Smith, 1978), which is really an EEG waveform ana lyzer. Correlations between the presence of seizures and predominant EEG backgrounds can be investigated. 277 2 40E CO 02 MOV 2,0 2410 4020 SZC DMASK,0 241 2 FC06 24 14 0280 Cl C,>0000 2416 00 00 241 a 1 3 04 JEQ POST 24 1 A a 0C2 C 2,3 241 C 13EF JEQ CUME 241 E 15EE JGT CUME 2420 1 ODF JMP NEG 2422 0589 POST INC 9 2424 8 0C2 C 2 ,3 2426 13EA JEQ CUME 2428 1 5E9 JGT CUME 242A 0284 Cl 4,TPLUS 242 C 00 3C 2 42E 15BA JGT PEAK 2430 0284 Cl 4,TMINU 2432 00 14 2434 116 7 JLT PEAK 2436 60CC 5 12 ,3 2438 8 8 03 y** V. 3 cDTHRES 243 A PC 04 2 43C 1 1 B3 JLT PEAK 243E CS03 MOV 3,3WSTA4 2 44 0 FE 00 2442 0420 3LWP SSTAT4 2444 FC1A $ * 3RD LC UP 2446 C0C2 POST 2 MO V 2,3 2448 0581 INC 1 2 4 4 A 0589 INC 9 2 44 C 0407 3LWP 7 2 44 E C802 MOV 2,3>IEF0 2450 1EF 0 2452 0405 3LWP 5 2454 CO 02 MOV 2.0 2456 4020 SZC 2MASK,0 2458 FC06 245A 0280 Cl 0,0000 2 45 C 0000 2456 1 3F3 JEG POST2 246 C C 8 04 MO V 4 S> WSTA2 2462 FDCO 2464 0420 3LWP 3STAT2 2466 FC 12 246 8 0289 Cl 9,TPLUS 246A 003C 2 46C 1598 JGT PEAK 2 46 5 0289 C I 9, TM INU 2470 00 14 2472 1 198 JLT PEAK 2 47 4 0204 LI 4, >800 2476 0300 2478 C 3 04 MOV 4 5 IEF2 247 A IEF2 2 47C 04 58 3 *1 1 END DATA POSIT? INCREASING? NCT MONOTONIC.GO TO NEG ZERO CROSSING REACHED PEAK? VAL-PEAK WITHIN LIMITS RESTART SEARCH GET PEAK-VAL AMP > 75 UV? RESTART SEARCH P ASS PEAK TO ST AT I S POS TO NEG ZERO CROS? PASS VAL-PEAK TIME TC WITHIN LIMITS? NO RESTART SEARCH S.W. RECOGNIZED OUTPUT NEG PULSE GO TO MAIM NO.OF ERRORS IN THIS ASSEMBLY^ 0000 RELOCATABLE LOCATIONS USED = 0000 STAT IS 13 used ARMA models to extract alpha, beta, delta, and theta band intensities for automatic EEG classifications. Mathieu (1970, 1976) used AR models in sleep scoring. Jones (1974) compared the window lag estimation with AR estimator to determine power spectral and coherence func tions in the neonate EEGs. Gersch and Yonemoto (1977) com pared the AR and ARMA models for the power spectrum estima tion of EEG. The same authors (1973) also used parametric AR models to apply the Shannon-Gelfand-Yaglom amount of information measure in sleep scoring. In epilepsy AR models have been used toanalyze the ictal event. Gersch and Goddard (1970) used the partial coherence among time series to extract information about driving. Tharp (1972) and Tharp and Gersch (1975) use a similar procedure to determine seizure focus. Herolf (1973) and Lopes da Silva et al. (1975) also use an AR model to perform seizure detections by inverse filtering. Another concern faced in spectral analysis is the presentation of the results. Due to the specific tech niques involved, the results need to be further processed to be readily interpreted by the neurologist. The main goals are data reduction and readability. The compressed spectral array of Bickford et al. (1973) has been exten sively used and summarizes pretty well the changes of EEG with some external factors (different behavioral states, tracking of alpha waves with light stimuli, hemispheric symmetry, acute slow wave abnormalities). The canonogram 54 it is possible to design bandpass filters using the direct substitution to map the filter poles and a technique derived from the bilinear to map the filter zeros. The relations between the design and the implementation are apparent in the choice of integer zeros (to cut computa tion time and roundoff noise) and in the choice of the transformation, which may help controlling the midband filter magnification. Finite Length Effects of the Implementation After obtaining the filter transfer function in the Z domain, the filter algorithm is the inverse Z transform of H(z). For a second order function of the special type which will be used here, H(z) can be written as in (15) and y(nT) as in (16). H (z) (z+1) (z2-Cz+D) (15) yn = Cyn-l-Dyn-2+xn-xn-l (16) The filtering will then be accomplished by computing (16) in the microcomputer. If an infinite wordlength machine is assumed, the implementation problem is nonexistent. Practically finite wordlengths must be used leading to the fact that almost every digital filter is nonlinear. For this reason the output of the digital filter deviates from what is actually desired. This leads to the study of various ways of 166 F3-A1, F4-A2, C3-A1, C4-A2. Besides the paper record obtained on a Grass Model 6, the sessions were tape recorded on an FM tape recorder Sanborn Model 3900, at 1 7/8"/sec. In the beginning and at the end of each ses sion a 50 yV, 10 Hz sinusoidal calibration signal was recorded on each tape (50"). Two types of logging were utilized: a ten-minute log describing the activities of the patient (awake, sleeping, walking, eating) and a tape log giving important information about tape number, beginning footage, recording conditions, EEG montage used, etc. Fig. 30 presents a copy of the technical log utilized. The specifications of the telemetry system, regarding signal fidelity, were judged appropriate (40 dB S/N ratio, 0.5-100 Hz frequency response) although some clipping occurred (200 yV maximum signal). The tape recorder band width at 1 7/8" was 0-256 Hz, with a signal to noise ratio of 35 dB. Some technical problems arose, mainly with the teleme try unit. Also some of the tapes did not have a proper calibration (very noisy signal), or the technician changed gains in the middle of the session to avoid clipping of the EEG. However, it is estimated that 90 percent of the 110 tapes have good quality data. One aspect that deserves mentioning is the lack of a time code recorded on the tapes. For synchronization purposes its information would be very convenient. 234 to help set dose levels of anticonvulsant medication. For instance, the patient is given a certain drug dosage, and the EEG is analyzed for the next day to assess the reduc tion in seizures. The generation of statistics for the detection param eters can be utilized to further quantify the PM patterns in the EEG. Up to now relations between drug administra tion and changes in the PM patterns have not been addressed. For research purposes, to assess how the drug is actuating and ultimately to understand better the PM generation, this aspect is important. Moreover, it may turn out that cer tain parameters (the repetition rate of the wave complexes is a good candidate) will be sensitive indicators of drug levels and could therefore be used to prescribe the dos ages. The above study requires new sets of data because it will be of paramount importance to quantify adequately the occurrence of the PM in humans before any anticonvulsant medication is administered. It will be important, for instance, to study the ultradian characteristics of PM activity. The correlation of seizure duration-interictal intervals with physiological states seems very interesting as the preliminary results of Stevens et al. (1971) and the observations on patient #4 show. The present system output can be used for this pur pose without modifications. To accomplish this study a 24-hour patient monitoring is desirable, which points to 326 22E8 04 C 3 CLR 8 2 SEA 0588 LOOP INC a 22EC C12A MOV 3>A( 10 ) 4 DEL X IN 4 22EE COOA 22FC 3908 MPY 8,4 N* QELX/DELY IN R4 22F2 3D2 A DI V 3>C( 10 ) 4 22F4 OOOC 22F C16A MOV 321 10),5 OLD X IN R5 2 2F a 0002 2 2FA 4 1E A SZC 3>Â£{ 10 ) ,7 DEL Y +? 22FC QCOE 2 2FE 13 03 JEQ POS 23C 0 062A DEC 14(10) NO, DEC OLDY 2302 0004 2304 10 02 JMP COMI 2306 0 5AA POS INC 24{10) 2308 00 04 22 QA 4 1 AA com SZC 3>E(10) .6 DEL X +? 230 C 00 OE 230E 1601 JNE POS 1 221 0 0504 NEG 4 2312 A 144 FOSl A 4,5 ADO DEL X TO X 2314 CA6A MOV 24(10),32(9) OUTPUT NEX X Y 22 1 6 0004 2313 0002 231 4 C 645 MOV 5, *9 221 C 0 69C GL *12 231 E 8A AA C 34( 10),33( 10 ) F INISHED? 2220 GO 04 2222 0 0 08 2324 16E2 JNE LOOP 222 6 CAAA MOV 36(10),32(10) UPDATE NEW X (ROUNDING) 2223 0006 2 32 A 0002 222C 0380 RTViP * SAME THING BUT DEL Y/DEL X 222 E C4C3 YS IG CLR 8 2330 0538 LCP INC 8 2232 C12A MOV 2>C( 10) ,4 2334 cooc 2336 3 908 MPY 8,4 2333 3D2A DI V 3>A( 10 ) .4 233 A 00 0A 233C Cl 6A MGV 24(10),5 233E 0 0 04 2 340 41 AA SZC 3>E( 10). 2242 00 CE 2344 16 03 JNE PO 32 2346 062 A DEC 22{10) 2343 0 0 02 234 A 1002 JMP CON 3 234C 05 AA P0S2 INC 32(10) TABLE VIIcontinued PATIENT #12 *Just one Sz > 3 sec 195 217 29 computation wordlength to obtain a certain output signal to noise ratio. In Chapter III the model for the petit mal activity utilized in the detection will be presented, along with its implementation on a microcomputer. The system's testing will also be explained there. As this research work is primarily engineering oriented (a new instrumentation system is designed around a new model of the petit mal activity), the system's evalua tion is presented in Chapter IV. Preliminary data demon strating the high resolution capabilities of the detector and its use in quantifying the petit mal seizure data will also be presented. 62 must be scaled so that the output, yn plish this (Oppenheim & .Shafer, 1975). As 00 88 I hk*n-k <24> k=0 where h^ is the filter response, it is easy to see that IyJ l IhkNxn-kI (max|xk|) Â£ |hk|. (25) k=0 k=0 00 Therefore, by making max|xv| = 1/ Â£ |h, |, overflow of y can K k=0 K n be prevented. This means on the other hand that the zero forming paths of the filter must be multiplied by CO k = 1/ Â£ |h |. Experiments have shown that this scaling k=0 k policy is too pessimistic and does not allow good use of the full dynamic range of the wordlength. A second approach would be to assume the input sinus oidal, with a frequency equal to the resonant frequency(es) of the filter. This is equivalent to requiring (H( ) is the filter's frequency response) |H(e^X) | < 1 -TT < X < ir (26) and is accomplished by multiplying the coefficients by k, the scaling constant, so that 285 248 132 131 21 133 13 aa a C 42 2* 110 70 1 5 16 7 17 13 C C C 60 61 62 6 63 C 140 141 00 132 1st,MN X1*2*1*1 hB8*wB/W0 u-2*Psecn*PT:irn*wBB**2 a 31 -wbh *? r cps crn **a-(piM cm **2) AAsSQPTC (SGPTCA**2+B**21i-A1/2.1 8R3QPT( C3QiT(A**?+fl**2)-A) / 2.5 xi2*r PPECII3 aWP*CWP5*PPECIT)-9P) PT^Cin*tiO*<--6B*PICTI)*AA) GO TQ 133.. 00 21 1*1. N , PfiCI3apREri)*Ja PIHfT3ariO*PInCT3*WP oo 134 rt.N PSCI3CMP|.XfPPEf II /PIKCT) ) P7(n*cExPP8nr cofJTr'ue wPI7Â£fj,aai FOUH AT 11H11 SITE(6r<0 CPSCI1. FnWMATClHn50X,7HP33,2F17.71 cnNTijunus xanO Pi33 FILTER TOa*O.Q nM*N 00 43 >1*1,100 TC*0>1PI_XC1 .0,0.0) WC9Â£M)CU0-WCU1Wfl.-FLOAT(M-l)*0.02)*2.+0 MCbHaCHPLXCO.OjWCPCm 00 42 J*1,N. rc*Tc*CWCPHnP3(J) WxCSH-CONJGCPSC Jl) 1 IF CNWIO.EG.Ol MNaN+1 WCDfa)*TC/wra(hl**(Nn-11 AC*CaP3 fWC8CM)**(MM 11/TO 1 IFCiC-T0ai4T,3,4t TOQtaC CONTINUE WRITE (6, M3) TOO FQPmATCIUQ.16H GAIN OF FILTER* r F16,4 3 00 70 H* i ,100 TKP in) ai r AM (A T.liC CWCO CM) ) /REAL (wco m 1 3 TXo (Ml SCABS f 1 ./ creu*wco con 3 ^3,1=1,11 wRITE C, 163 FOUMATItHO) FOHH!rUoY.?WPLQ'r OF SAMO PASS FILTER) call plotcwra*tkb.tool 00 .7 1*1.10 ipitec, m FOHKATC1HO) EVAI UiTIOM nF IMS mAGmITUOE A NO PHASE p THF OIGITaL FUTES USING OIHECT mapping of PnLFS A'JQ SEARCHING SYMMETRY ACROSS THE BANG CASCADE S F AI IZATION . Ka(l .+C03 r40*f31 / f 1 .-COS C'40*T3 3 IFCK+1-M16G,61,62 ko*n-k*i M 3K <0*0 Kt*K <0*0 X!N-1 <2*1 IE CNWIO.Ea.33 X2*2 FOHHa^^i isHNiiMBEP OF ZEROS at THE CRIGTn AMO Z*-1) WPITE Co,6*1XQ,K1 F0WMaT(IHO.30X,3HX0s, 2.30X, 3HKla, l'i) GO TO 147 HTGH PASS toamsfophattcn 00 141 t 1,M . PS(IIS40U/P3(13 PZ(i5*cexP{PsCi)*T) TCC-0,0 00 142 1*1.100 rr*CMPLX(i.0*0.0) 279 Yes 254 j^INTHANDLER ^ / / / / / / / / V t Initiate A/D > > Set sampling frequency \ > ^ Return ^ Fig. II-2. Flow Chart of Filter 291 microcomputer and a 12 bit A/D converter were utilized. With the procedure, a fourth order bandpass filtering func tion can be accomplished in less than 300 ys. Time domain parameters were selected to arrive at a PM seizure model which would be applicable to the detection of classical PM and PM variant epilepsies. The parameters were translated into electronic quantities and implemented in a microcomputer. For the first time in the automated analysis of EEG, a microcomputer based system was built for the processing of one channel of data, using a completely digital, real time, detection scheme. All the programs which comprise the detector were written in assembly language, and care was taken to reduce the interactions between program modules. The top-down program approach was utilized. The system was tested and evaluated with data from epileptics collected in the Veterans Administration Hospi tal, Neurology Service, from an ongoing drug study. A telemetry link was used to record the EEG, allowing the patient to move at will in a 3 x 5 meters room. Six patients who showed the largest number of seizures were selected for the evaluation. A total of 70 hours of data was analyzed. The agreement in the detections with the human scorer, for seizures (Sz) greater than three seconds, was 86 per cent. If the sorting of seizures in the groups 3 < Sz vii 342 Rabiner, L., and B. Gold. Theory and Applications of Digital Signal Processing. Prentice-Hall, Englewood Cliffs, New Jersey, 1975. Remond, A., and B. Renault. "La Theory des Objects Elec- trographique." R. Rev. EEG Neurophysiol., 2:241-256, 1972. Remond, A. "Reviews in Epilepsy." Hndbook of Electroen cephalography and Clinical Neurophysiol., 13a, 1974. Robinson, E. A. Statistical Communication and Detection. Hafner, New York, 1967. Rosadini, G., Cavazza, B., and F. Ferrillo. "On the Organ ization of Electroencephalographic Rhythms in the Sleeping Man." Acta Neurol. Lat. Amer., pp. 200-217, 1968. Rose, S. W. "Reliability and Validity of Visual EEG Assess ment in Third-Grade Children." Clin. EEG, 4:197-205, 1973. Saltzberg, B. "Detection of Focal Depth Spiking in the Scalp of EEG Monkeys." Electroenceph. Clin. Neuro- physiol^, 31:327-333, 1971. Saltzberg, B. "Theoretical and Experimental Investigation of the EEG Signals Using Parameter Tracking and Matched Filtering." Doctoral Dissertation, Marquette Univ., Milwaukee, 1972. Sandberg, I. W., and J. F. Kaiser. "A Bound on Limit Cycles in Fixed Point Implementation of Digital Fil ters." IEEE Trans. Audio Electroacoust., AU-20:110- 112, 1972. Schenk, G. K. "The Quantification of EEG by Vectoral Iter ation Techniques: A Simulation Method of Visual EEG Analysis." Electroenceph. Clin. Neurophysiol., 37:106, 1974. Schenk, G. K. "The Pattern Oriented Aspect of EEG Quanti fication." In Quantitative Studies in Epilepsy. Raven Press, New York, pp. 431-462, 1976. Schmidt, R. P., and B. J. Wilder. Epilepsy. Davis, Phila delphia, 1968. Sellden, U. "Psychotechnical Performance Related to Parox ismal Discharges in the EEG." Clin. Electroenceph., 2:18-27, 1971. 175 (computer versus human) were built for each patient. They included the sorting of seizures (Sz) in three categories (l The computer agreement with the detection of seizure (first line), the computer agreement with the sorting (second line), and the total time in agreement (third line) are also presented. The first information answers the question, "did the computer recognize this seizure?" More important, however, is the second parameter, which is related to the aspect of sorting, since it has a higher physiological significance. To assess the reliability of the computer performance as the "leading" output, not only the agreement with human is important but also the false detections. The picture of computer performance is acquired from the three last lines of the tables. Table IV presents the totals for each patient analyzed. The purpose is to analyze any dependence of performance with the pattern (e.g., PM classical versus PM variant). The computer performs well with respect to sorting in the classical PM patterns as expected, since the key aspect of these is regularity. Take, for instance, #7, where the agreement is 100 percent in Sz>10 sec, and #12, where the agreement is very good for 3 total time). Patient #5 also had an excellent agreement in the two first sessions, but in the third (after Depakene administration) there was a greater variability in the PM 156 rate were recognized), up to the last slow wave before the positive going pulse is outputted. This pulse means that no slow waves were recognized in the last second or that the period does not meet the requirements twice, consecu tively. The statistics are stored in memory beginning in RAM address 2720. The format of the seizure data is the fol- ri lowing, for seizure > 3 sec: FFFF (SEPARATION TAG) Time of occurrence (sec) Duration (# of points) (X) repetition period (# of points) (a ) variance of repetition period (# of points) (X) 1/2 period of waves (# of points) (a ) variance of 1/2 period of waves (# of points) (X) delay spike/S.W. (# of points) 2 (a ) variance of delay spike/S.W. (# of points) (X) amplitude of filtered S.W. (level of A/D conversion) 2 (a ) variance of amplitude of filtered S.W. (level of A/D conversion) (X) amplitude of filtered spike (level of A/D conversion) 2 (a ) variance of amplitude of filtered spike (level of A/D conversion) FFFF (SEPARATION TAG) For seizure less than 3 sec only the two first parameters are stored. The system output also includes the number of seizures (Sz) detected in each of the categories Sz<3 sec, .945 sec o 0.062 sec Fig. 44b. Variance Repetition Period Versus Repetition Period (#5-2) 220 49 Another conclusion of the analysis which is very important regards the filter passband gain. In digital filter design the control of the filter midband gain is generally accomplished by properly scaling the filter transfer function. It is suggested that the correct handling of the Z plane filter design may constitute an alternate way to control the midband gain, at least in some cases. An example given in (Principe et al., 1979) is worth mentioning. For the delta filter (0.1-3 Hz, fs=100 Hz), the midband gain was decreased from 878 for the filter designed with the narrowband direct substitution to 87 just by using the wideband transformation (both fre quency responses met the specifications). Let us analyze in general the effect of the zeros at z=l in the transfer function. For the zero at z=l the numerator takes the form (10) and for the zero at z=-l |z+1| = 4 cos2 (11) It is concluded that the effect of the zero at z=-l has the symmetric effect of the zero at z=l, with respect to fs/2, so whatever conclusion is reached for one will apply for the other in the symmetric region. Fig. la shows the attenuation characteristics of EQ(10),(11). Now suppose 5 fact, the procedure displays more points of the discrete Fourier transform (DFT), and so a better approximation to the Z transform is achieved. However, the Z transform remains the same (as long as the record length is con stant) so the initial resolution is kept unchanged. Another problem encountered is the estimation of the power spectrum. Let us resume the estimation process. To compute the essential statistical properties of a Gaussian random process the more general way is to evaluate the autocorrelation function by ensemble average. The method is computationally expensive since various sample sequences of the process have to be available. If the random process is ergodic (Lee, 1960; Thomas, 1969), then it is correct to compute the autocorrelation function by a time average only over one sequence. There are fewer memory requirements involved, but as the properties of the random process are to be inferred from the ones calculated in the sample sequence, consistent estimators have to be used (i.e., estimators that are unbiased and whose variance decreases with the number of data points). It is very easy to find consistent estimators for the autocorrelation functions (Oppenheim and .Shafer:,. 1975). The problem arises when the goal is to estimate the power spectrum. It is not true in general that the Fourier transform of a consistent esti mator of the autocorrelation function (periodogram) will produce a good estimate of the power spectrum (Barlett, 1953; Jenkins & Watts, 1969). When an infinite number of 183 the seizure count due to desynchronization, still detecting slow waves with the right repetition rate, but not calling the event a new seizure due to the absence of spikes. Fig. 32 exemplifies the explanation. Probably more impor tant than the discrepancy in seizure duration is the appearance of a large number of misses in this group of patients (16.7 percent). They occurred mainly in patient #3 and were due to the nonexistence of well-defined spikes. This fact probably calls for a modification in the detec tion algorithm. Table V presents the final agreement averaged over the six patients studied. The discussion will be directed towards the clinical significant seizures, i.e., seizure > 3 sec. The agreement of the computer versus the human with respect to seizure detection is 100 percent for seizure > 10 sec and 88 percent for 3 < Sz < 10 sec. The agreement for Sz > 3 sec is 91 percent. If the sorting of the sei zures in the categories is taken into consideration, the agreement drops to 86, 73 and 77 percent, respectively. The total time in agreement is 83, 94 and 89 percent, respectively, for 3 < Sz < 10, Sz > 10, Sz > 3 sec. How ever, as these numbers are totals, positive and negative errors can compensate each other. Therefore, some caution shall be exercised in their interpretation. The computer counts are also presented in Table V. They are intended to judge the system usefulness as the only scorer, now that the agreement with the human has been 128 energy (spikes), the output of the filter will tend to display the effect of the impulse response decay, shifting negatively the zero (d.c.) level. To detect the slow waves in the four subjects, the half period requirements were set at High period ... 0.083 sec (14) 11 Low period ... 0.200 sec (30H) which corresponds to frequencies of 6 to 2.5 Hz. Further in the testing it was found that this measure (peak-to- valley) alone is very sensitive to artifacts since the response of the slow wave bandpass filter to a step input mimics exactly the period window. Assuming that the band pass filter preserves the lowpass characteristics, but introduces a time scale warp (Blinchikoff, 1976), the rise time was found to be 0.14 sec. Therefore, a maximum period of 0.5 sec (3C) between the zero crossings was also Â£1 required (points B and D). Periodicity of the PM Complexes If the filters were good enough to attenuate the spike component, only this measure would be necessary to detect the PM slow wave component. However, only higher Q, or higher order filters would be necessary, giving heavier weight to the undesirable ringing phenomena. The period of the slow wave component was measured between peaks, since these points are detected in the pre vious measure. The periodicity window was set between 16 prediction scheme. He implements the inverse filter of the all-pole approximation model of the EEG and monitors the error between the predicted value and the real data points. Deviations above a prescribed threshold are correlated with nbnstationarities of the input EEG data. Very similar techniques are used in EEG adaptative segmentation proce dures (Bodestein & Praetorious, 1977) which could eventual ly lead to transient event detection. A combination of the nonstationary information of the phasic events and the sharpness of the spikes was combined by Birkemeier et al. (1978) to increase the resolution separation of the clusters of epileptic transients and background activity. The system requires settings for each patient and works better in epochs with few spikes (otherwise they may bias the estimation). It is sensitive to artifacts. Etevenon et al. (1976) use null-phase inverse filter ing to deconvolve the EEG and then detect fast activity. His method is equivalent to the one just described (Lopes da Silva), but it is executed in the frequency domain. From the EEG sequence the power spectrum is calculated through the periodogram, and after smoothing and whitening the amplitude spectrum is obtained by evaluating the square root. Then the inverse of the spectrum is computed, multi plied by the Fourier transform of the incoming EEG, and the inverse FFT taken. This signal accentuates sharp tran sients contained in the EEG, and after threshold detection they can be extracted. 136 ia = w-* sin w(tw-At) = w-* sin wt. max cos wAt + A cos wt, max sin wAt = \axU"cos wAt) (54) and At = l/2fs. For fs = 240 Hz and w = 2ttx10 AA (1-0.9914) i.e., the error is better than 0.9 percent. The repeata bility is also good, because an average over 64 cycles is performed. Fig. 21 shows a typical output. System Implementation The PM seizure detector was implemented around a TMS 9900 16 bit microprocessor. The system's block diagram is shown in Fig. 22, and it incorporated the TMS 990/101-M microcomputer board, TMS 201/2 memory expansion board and the Analog Devices RTI 1240 I/O board. The I/O board, which plugs directly in the TI bus in a memory mapped con figuration, is a 16 channel multiplexed, 12 bit D/A and two, 12 bit D/A.channels. The memory mapping of the TI micro computer as used is shown in Fig. 23. 135 beginning a program module looked at the peaks of the calibration sinewave and computed the average over 64 peaks. The number obtained corresponds to 25 yV (the calibration sine is 50 yV p.p.). Then, knowing the mid band gain for each of the filters and the desired threshold level, the amplitude at the filter output that corresponds to 75 yV can be immediately obtained by a multiplication (75/25 x gain). This new number was stored in a memory location used in the program as the threshold. The cali bration was then performed without operator's interven tion, reducing one important step of possible errors. The precision in the determination of the peak was very good, since the sampling frequency is 240 Hz and the calibration sine is 10 Hz. The maximum errors can be computed by the following derivation (complemented by Fig. 20): Fig. 20 45 sec t 202 224C 2 24 E 2250 2252 2254 2256 2253 225 A 225C 225E 226 0 2262 2264 2266 2263 2 26 A 226C 226 E 2 27 0 2272 2274 2 27 6 2278 2 27 A 2 27C 227 E 2230 2282 2234 223 6 2233 229 A 2 28C 228 E 2 29 0 2292 229 4 2296 2 29 8 229 A 2 29C 2 2 9 E 22A0 2 2A2 22A4 22A6 2 2A8 2 2AA 22AC 22AE 22 EG 22B2 2284 2 236 04E0 CLR 3FLAG NO TOT FC22 0201 LI 1 ,PC1 ADDR 0 FCOS 02 02 LI 2 >24A 8 JUMP 24A8 CC42 MOV 2 *1 + 0202 LI 2 >24C0 2ND 24C0 CC42 MOV 2 l + 0202 LI 2 >2400 3RD 2 4D0 CC42 MOV 2, *1 + 02 02 LI 2 >FD AO WS P F DAO CC42 MOV 2 *1 + 0202 LI 2 ,>2590 PC 2590 CC42 MCV 2 ,*1 + 0202 LI 2 >FDC0 WSP FOCO CC42 MOV 2**1 + CC61 MCV 3-4( l ) ,* 1 + 1 FFFC 0202 LI 2 >FDE0 WSP FDEC CC42 MOV 2**1 + CC 61 MOV 2>-4( I ) .*1 + FFFC 0202 LI 2 *>FE0 WS 1 OFEC CC42 MCV 2**1 + CC 61 MOV 3-4 0202 LI 2 >FE20 wsp FE 20 CC42 MOV 2, *1 + CC61 MOV a-4(i), i + ! FFFC 0202 LI 2,> OFFF MASK ! OFFF C802 MOV 2.3MASK4 FC26 02 02 LI 2 >OQFF MASK ! OOFF C802 MOV 2,3 M AS K 8 FC 23 0205 LI 5,>7FFF MA SK 7FFF C 3 05 MOV 5 3MA SK FC06 0300 LI MI 3 I NT R MASK ' 0003 02 05 LI 5,WPIKE WSP 2ND JUMP TO SPIKE RD JUMP TO SPIKE SP FCR ST AT IS 1 PC FCR ST AT IS 1 WSP FOR STATS 2 PC FOR STATS 2 FOR STATS 3 PC FOR STATS 3 FOR STATS 40 PC FOR STATS 4 FOR STATS 5 PC FOR STATS 5 MASK FOR CONVERSION 238 As a concluding remark it is appropriate to praise the use of microcomputers for the versatility given to the detection scheme. The great number of changes in the cri teria performed during the detector design were easily accomplished by correcting programs, not cutting wires. With adequate software developing facilities the design time can be shortened, but better yet, it enables updating or modification of features in a final product. For instance, suppose that the criteria of detection are judged too loose for a specific application (e.g., study only classical PM patients), a simple modification in one of the program modules will change the criteria from three slow waves and at least one spike to spike/slow wave three times, decreasing appreciably the false alarm rate. Another example is to customize the PM parameters for the pattern of a certain patient. The power of the pattern recognition algorithm and the on line statistics generation also attest the advantages of the use of microcomputers. It is our belief that the system presented in this dissertation, as well as its future developments, will become a valuable tool in PM epilepsy research and also in medical diagnosis and health care. 90 filters with only 4 multiplications per cycle, allowing sampling frequencies up to 3 KHZ. Signal to Noise Ratio Specification Now that the output noise can be calculated and the filter magnification is known, the important step of choos ing the hardware to obtain a given output signal to noise ratio can be addressed. From the specification involving the output signal to noise ratio and the degradation pro duced in the filtering, the input signal to noise ratio can be calculated, which sets the input wordlength. If the analog dynamic range is smaller, an A/D converter that spans the analog dynamic range shall be chosen and a shift left operation up to the requirement of the input word- length shall be performed. From the internal magnification analysis the increase in wordlength can be determined and added to the input wordlength to set the computation number of bits. In the present form a 12 bit A/D converter and a 16 bit microcomputer are being utilized, allowing for the implementation of filters with gain of 16 or less. There may be cases where a longer computation register length may be desirable, but on the other hand a 12 bit A/D converter may be excessive for EEG processing, mainly when the data are prerecorded on analog tape. The dynamic range of a Sangamo 3500 is 46 dB at 60"/sec. (manufacturer specifica tions) which suggests that an 8 bit A/D is sufficient. 332 234E 0002 2350 4 1 EA C0N3 SZC Â£>E( 10) 7 2352 OOOE 235 4 1301 JEQ P0S3 2356 0504 NEG 4 2358 A144 P0S3 A 4,5 235A C66A MOV 32(10) ,*9 235 C 0002 235E CA45 MOV 5,32(9) 2360 0002 2362 069C BL *12 2364 SAAA C Â£2(10) ,Â£6(10 ) 2366 0002 2368 00 06 236A 1 6E2 JNO LOP 236C CAAA MOV 23(10),Â£4(10) 23E 0008 2370 0004 2372 0380 RTWP 2374 41 AA FLAT SZC Â£>E (10,6 CHECK SIGN OF DEL X 2376 OOOE 2378 1 6C3 JNE PCS4 237A 0 62 A DEC 32(10) 237C 0002 237E 1002 JMP COM4 2330 05AA PC S4 INC Â£2<10) 2382 0002 2384 C66A C0N4 MOV Â£2(10),*9 OUTPUT X,Y 238 6 0002 2333 CA6A MOV 34{10),Â£2(9) 233A C 004 23SC 0002 238E G69C 3L *12 2330 8AAA C 22( 10) ,Â£6( 10) 2 392 0002 2394 00 06 2336 1 1 EE JLT FLAT 239 3 0380 RTWP * INTE NSIFY CCRU BIT 26) 239A CA 8C MOV 12,Â£>16i10 > SAVE VALUE CF R12 239 C 00 16 239E 020C LI 12,>120 2 3A0 0120 23A2 1 D08 S30 8 OUTPUT *5 V TO CRU 2 2A4 C53C WIDTH INC 12 FORM 4 PULSE 2 3A6 028C Cl 12,>140 23A8 0140 22AA l 1FC JLT WIDTH 23AC 020C LI 12,>120 23AE 0120 2330 1 E08 SSZ a OUTPUT ZERO TO CRU 23B2 C32A MOV 2>16(10),12 233 4 0016 2 336 Q4oB E *11 END )* NO.OF ERRORS IN THIS ASSEMBLY* 0000 JF RELOCATABLE LOCATIONS USEC = 0000 337 Gersch, W. "Spectral Analysis of EEG by Autoregressive Decomposition of Time Series." Math. Bioscience, 7:205-222, 1970. Gersch, W., and D. R. Sharpe. "Estimation of Power Spectra with Finite Order Autoregressive Models." IEEE Trans. Automatic Control, AC-18:367-369, 1973. Gersch, W., and J. Yonemoto. "Parametric Time Series Models for Multivariate EEG Analysis." Computers and Biomedical Research, 10:113-125, 1973. Gersch, W., and J. Yonemoto. "Automatic Classification of Multivariate EEG Using an Amount of Information Meas ure and Eigenvalues of Parametric Time Series." Com puters in Biomedical Research, 10:297-318, 1977. Gersch, W., and G. Goddard. "Locating the Site of Epilep tic Focus by Spectral Analysis." Science, 169:701- 702, 1970. Gevins, A. S., and E. L. Yeager. "EEG Spectral Analysis in Real Time." In DECUS Spring Proceedings, pp. 71-80, 1972. Gevins, A. S., Yeager, C. L., Diamond, S. L., Spire, J. P., Zeitlin, G. M., and A. H. Gevins. "Automated Analysis of the Electrical Activity of the Human Brain: A Progress Report." Proceedings IEEE, 63:1382-1399, 1975. Gevins, A. S., Yeager, C. L., Diamond, S. L., Zeitlin, G. M., Spire, J. P., and A. H. Gevins. "Sharp Tran sient Analysis and Threshold Linear Coherence Spectra of Paroxisms." EEG in Quantitative Studies in Epilep sy. Raven Press, New York, pp. 463-482, 1976. Giannitrapani, D., and L. Kayton. "Schizophrenia and EEG Spectral Analysis." EEG Clinical Neurophy., pp. 377- 386, 1974. Gold, B., and C. M. Rader. Theory and Application of Digital Signal Processing. McGraw-Hill, New York, 1969. Gold, B., Lebow, I. L., Hugh, P. G., and C. M. Rader. "The FDP a Fast Programmable Signal Processor." IEEE Trans. Computers, C-20:33-38, 1971. Gotman, J., and P. Gloor. "Clinical Applications of Spec tral Analysis and Extraction of Features from Electro encephalograms with Slow Waves in Adult Patients." Electroenceph. Clin. Neurophysiol., 35:225-235, 1973. 97 The most striking are the following: definition of begin ning and end of seizure, uncertainties related to the fading of the pattern in the middle of an ictal discharge (single seizure?), modifications of the nature of the pat tern to polyspike and wave or only slow waves (should the count be terminated and restarted?). It is interesting that only one paper was found which dealt exclusively with the morphology of the PM wave com plexes. This advanced the idea of trying not only the detection of the PM paroxysms but also their analysis. Weir (1965) analyzed, using a storage scope, the morphology of the spike and wave complex and concluded that "they are composed by a small short duration negative spike followed by a prominent positive transient from which the second, larger amplitude negative spike originates. The negative wave that follows is not a complete dome but approximates the first three quarters of one (Fig. 10)." He points out the masking effect of the EEG recording time constant for the lack of widespread recognition of these elements. He characterized the complexes as follows: Spike one is small amplitude (25-50 yV), negative and of short duration (10 ms). It has a triphasic potential when considered with the preceding activity. The positive transient is a sharp positive deflection with duration of 100-150 ms. Spike two appears in the positive portion of the complex and was gen- erally of high amplitude (100-200 yV) and with duration between 40 and 90 ms. The amplitude of spike two waxed and I certify that I have read this study and that in my opinion it conforms to acceptable standards of scholarly presentation and is fully adequate, in scope and quality, as a dissertation for the degree of Doctor of Philosophy. Graduate Research Professor of Psychology This dissertation was submitted to the Graduate Faculty of the College of Engineering and to the Graduate Council, and was accepted as partial fulfillment of the requirements for the degree of Doctor of Philosophy. August 1979 Dean, College of Engineering Dean, Graduate School 131 The weight attributed to them will be reduced in this work. There was even an attempt of designing the slow-wave detec tion without an amplitude threshold, but surprisingly there was a lot of EEG activity detected as shown in Fig. 18. When an amplitude threshold of 25 yV was utilized, the only activity besides PM detected in the awake records was the eye movements. In the neurophysiology literature the amplitude of the slow waves is said to be at around 100 yV; hence, it was decided to increase the threshold to 75 yV. It is here acknowledged the desirable effect of increasing the threshold for the rejection of high energy artifacts. The amplitude of the spike threshold was also set at 75 yV. The importance of this threshold is higher than for the slow waves, since the only other parameter used to detect the spike transient is the period. This value dis plays a good rejection to most of the muscle activity present in the EEG when the subject is not moving around. The recordings were obtained with the patients awake for the great majority of the time, and therefore the EEG is expected to be contaminated with muscle activity. Fig. 19 shows a typical performance of the detector in a muscle contraction, subject not moving. One other problem faced was the determination of the amplitude levels at the filter output, since these filters possess gains greater than one. The microcomputer and its arithmetic facilities enabled an elegant solution. In the 112 other parameters of the digital waveforms are necessary, like peak amplitude or slope measurements, because for interpolation every parameter utilized needs a reconstruc tion method. For the detection parameters that are sought a sampling frequency of 240 Hz seems a good choice, since the preci sion on the period measurement of the spike is kept around 1 Hz (10 percent), and the peak amplitude can also be determined with a reasonable accuracy. The other variable that controls the choice of the sampling frequency is the frequency content of the data. The concern is not the accuracy of representing the EEG waveforms but the aliasing effect of muscle activity. It is known that electrical potentials from the cranio-facial muscles have frequency components up to 150. Hz (Smithy 19.79)., while for the . EEG the frequency content can be assumed below 40 Hz. For the representation of muscle activity a sampling frequency of 300 Hz would be necessary, but the sampling frequency can be selected below this quantity as long as the muscle components are folded back to a non-EEG frequency range. For the numbers given above, any sampling frequency above 190 Hz will accomplish the folding to a non-EEG range. This fact is shown in Fig. 13. Selection of the PM Filters Characteristics The purpose of using EEG linear filters in the EEG detection is to attenuate out-of-band activity. Care must I 272 25 yV Fig. 18. Slow Wave Performance. No Threshold. 132 328 TITL ,oL0T * 44 4# #*.4.4 44 4444444 4 4444 44 44444:44 444 * PLOT 4444 444 4 44 *44444 4 4444 444 4 4 44 44 444 * MODULE THAT DISPLAYS POSITIVE * DATA INTO THE TEXT 611 DISPLAY * X DATA I S ASSUMED I N RO * M A X X It ft It R 1 * Y DATA ft It W R2 4 MA X Y Vt t n R3 4 X AXIS IS D/A I 4 * REG 9 IS D/A CONVERTER POINTER * REG 10 IS EASE ADDR * REG 12 NEEDS TO BE LOADED WITH * BEGINING ADOR CF PLCT+39A * * SEE DATA PROGRAM FOR CONSTANT DESCRIPTION * 2200. AORG >2200 1 EFO ADCX EQU ' >1EF0 1EF2 ADCY EQU >1 EF2 FCOA CONST EQU >F C04 FC06 OLDX EQU >FC 0 6 FC 03 CLDY EQU >FC0 8 FCOA NEVrX EQU >F COA FCOC NE'rt Y EQU > FCOC FC OE XDIF EQU >FC0 E FC10 YD IF EGU >FC 1 0 FC12 MASK EQU >FC12 FC 1 A SCALE EQU >FC 1 4 FC 16 LIM EGU >FC 1 6 22 0 0 C12A MOV 2>l 0 CIO ) ,4 R4 GETS XO 2202 00 10 220 A CA44 MOV 4,32C9) OUTPUT xo 2206 0 0 02 2203 C 644 BACK MOV 4, #9 OUTPUT YO 220 A 06 9C BL *1 2 INTENSIFY BEAM 220C 0 5 84 INC 4 220 E 8 A34 C a, a>12c i 0) PEACH THE END? 2210 00 12 22 12 1 5F A JNE BACK 221 A Ci 2A MOV 2> 1 0 C 10 ) ,4 PLOT Y AXIS 2216 00 10 2218 C644 MOV 4 ,*9 221 A CA44 EAQ MOV 4,32(9) 22 1C 0 0 02 221 E 0 69C BL 4 12 2220 0534 INC 4 Fig. 31continued 173 TABLE VIIcontinued PATIENT #16 Session (EEG Mon tage) Repetition Period (sec) Half Period Slow Waves (sec) Delay Slow Wave Spike (sec) Filtered Slow Wave Amplitude (yV) Filtered Spike Amplitude (yv) In Sz In^sT Sessiori\^ X a X a X a X a X 0 \ 0.13 \ 0.03 \ 0.11 \ 6.08 V 5.83 1st V1, 0.37 0.7o\ 0.12 0.04 0.19 0.1l\ 141 5.54V 128 4.6S\ \ 0.10 \ 0.04 \ 0.38 \ 5.08 \ 3.86 2nd 0.31 0.13 \ 0.17 122 106 \ 0.04\ 0.2l\ 7.27\ 6.6r\ 3rd* \ \ *Bad tape 196 152 was implemented. Instead of performing the shifting before the division by N, it seems equally time-consuming, but much simpler, to perform an extra division using the shifted left remainder. Proper care was taken to shift equivalently x. 2 It may happen that Ex^ is so big that the division by N will yield a number greater than the maximum represent able number in a 16 bit machine. In cases like this, the number of occurrences (N) is shifted left 4 bits, and the variance is then calculated with a magnitude of one hexa decimal digit smaller. Although this operation degrades the resolution, it is the only way to perform the computa tion in fixed point. FILTER. This is the module that performs the fourth order bandpass filtering in the range 10-25 Hz and 0.8-6 Hz with sampling frequencies of 240 and 80 Hz respectively. Control is passed through a BLWP using register 8, 9 of MAIN. The modules are 4 independent workspaces (one for each resonator). A detailed description of the filter program is given in Appendix II. When control is passed to FILTER, there is a wait loop which synchronizes program flow and the sampling frequency. Control leaves the loop through a level III interrupt. The 80 Hz sampling fre quency is obtained by calculating the filtering for every one out of three input points. INTHAND. This module sets the sampling frequencies through an interrupt. The TMS 9901 is used as a timer, and 226 with automatisms that produce muscle potentials (e.g., eye blinking, teeth grinding). The inclusion of further proc essing channels to use the bilateral synchrony property of the PM pattern also fails, because chewing is also syn chronous in the frontals. Of course, a brute force solu tion is to include a push button on the detector which will disable the seizure detection. However, this may lead to a lot of errors (accidental triggering, forgetting to push the button, etc.). The inclusion of the statistics as'an output and the constancy of the PM parameters suggest another method to deal with the problem. Only by chance does the chewing possess the same characteristics as the PM patterns; there fore, further processing of the statistics data may elimi nate the false detects in chewing (and the other types of artifacts). The most sensitive parameters seem to be the repetition rate and its variance. Figs. 49a and b show examples of false positives in the plot time of occurrence/ repetition rate/variance of repetition rate for patients #7-1 and #7-2 and Fig. 49c for patient #12-1. The problem as enunciated is a classical communications problem: knowing the a priori probability of false detection (from the evaluation) and assuming gaussian probability density functions for the detection, a boundary can be easily eval uated with the seizure data (Helstrom, 1964). Additional information can be supplied taking into consideration the fact that no false detections greater than five seconds 276 23AE 2 3B0 FC00 1 504 JGT CLEAR 2382 0203 LI 3 > 140 YES.LOAD SAFE VALUE 23B4 228 6 2383 238 A 0140 C803 FCOO 0581 CLEAR MOV I NC 3.3DELAY l INC R. RATE COUNTER 2 23 C C 0 41 MOV 1 1 OVERFLOW (NEG)7 23 BE 22C C 1502 0201 JGT LI OKAY 1>140 LOAD SAFE VALUE 22C2 23C4 0 1 40 04 07 CKAY LVKP 7 GET DATA POINT 22C 6 C802 MOV 2 *2> 1 EFO OUTPUT IT 23C8 2 3CA 1 EFO CAEO CLP 3>1EF2 CLEAR D/A II 22CC 23CE 1EF2 0208 L I 3. >208C SKIP INITIALIZATION 2 2D0 2 20 2 20 8C 0405 BL'aP 5 GO TO SPIKE 2304 C 002 MOV 2.0 GET SIGN CF INPUT 220 6 2308 230 A 22DC 230 E 23E0 4020 FC06 0280 8000 16E2 C0C2 4 * * NEG SZC 3MASK.0 c.i o,>aooo JNE PEAK 1ST LOOP MOV 2 ,3 SAVE DATA IN R3 22E2 0581 INC 1 INC R. RATE COUNTER 23 EA 05 AO INC 3DELAY INC SPIKE DELAY 2 2E6 23E3 23E A 23EC 2 3EE 23F0 FCOO 0407 caoa 1 EFO 0405 S0C2 8L WP MOV BLViP C 7 2,3>1EF0 5 2. 3 REACH NEG PEAK? 2 2F2 2 3F4 2 3F 13F6 1 1F5 04 C4 JEQ JLT CL R NEG MEG 4 YES CLEAR VAL-PEAK COUNTER 23F 8 04C9 CLR 9 CLEAR ZERO CROC COUNTER 2 3F A C303 MOV 3.12 SAVE VALUE OF VALLEY IN RI2 2 3FC C0C2 * * CUME 2ND LG MOV CP 2 .3 2 3FE 2400 2402 2404 2406 24CS 240 A 240 C C 584 0581 C5A 0 FCOO 0407 C 8 02 l EFO 0405 I NC INC I NC 3LWP MOV QLWP 4 1 3DELAY 7 2 3>1EFO 5 39 bandpass filters (Smith, 1978). The present design uses a slightly underdamped frequency response. It is hard to qualify the desirable properties of the filters since there is no general theory of EEG detection. The only possible criteria are to understand the filter function in the detection scheme and hopefully arrive at some rules that will serve as guidelines for filter selection. The purpose of using EEG filters in sleep studies is to attenuate out-of-band activity to enable further signal processing (e.g., zero-crossing). On the other hand, the filters shall not mask out-of-band activity so that it will look like in-band activity for the rest of the processing algorithm. Every filtering function produces a certain masking, since the filter output is the convolution of the filter's impulse response with the input. However, the weight of the smoothing is controlled by the Q of the fil ter, hence the use of low Q filters. At present no other characteristics of the filters, like phase distortion or group delay, have been brought into the picture of filter selection (Smith, 1978). The in-band amplitude character istics have also been neglected due to the wide variability of the amplitude of the EEG. The above-mentioned analysis refers to EEG sleep studies. For the epilepsy application the knowledge is still less quantitative due to less expe rience From the study of the petit mal pattern it can be concluded that the fast attenuation characteristics are of paramount importance for the slow wave filter (activity as Fig. 15d. Spike Impulse Response. 256 Points. 124 44 s 2 z-1 T Z+l (6) The particular type of isomorphism (an infinite length line mapped into a finite line) implies a nonlinear mapping, in the sense that a Az has different correspondence to a Aw, depending on the value of z. As a matter of fact, the transformation follows a tangent rule as can be easily shown (Oppenheim & Shafer, 1975). Therefore, to get the correct frequency characteristics in the Z domain, the S plane filter critical parameters (center frequency, band width) have to be properly prewarped (Childers & Durling, 1975). Nevertheless, a piecewise linear response is mapped onto a piecewise linear response, preserving the boundaries. It is important to point out that the prewarping does not change the nonlinear properties of the transformation, but compensates for it. So a distortion of the phase response of the filter will always be present. One of the most interesting properties of the bilinear transform is the absence of aliasing. The theoretical explanation is very easily understood. As the entire jw axis is mapped onto the unit circle, the folding frequency (z=-l) corresponds to w=, and the attenuation of any realizable network shall also be infinitely large. How ever, from the Z domain point of view the reasoning does not seem explanatory since the Z axis is the same whether the filter is mapped using a rational or a nonrational transformation. In the particular case of a filter, the ACKNOWLEDGEMENT S The author would like to thank his supervisory committee for their aid, guidance and criticism during the preparation of this dissertation. Especially the author would like to thank Dr. Jack Smith for his willingness to share facilities, personal contacts and responsibilities. It was an unconventional but effective learning experience, which went beyond the restricted field of EEG studies. The friendly relationship created during these years is deeply appreciated. There are many others who contributed their knowledge, skills and time for this dissertation to become a reality. The author would like to thank Sonny, Yvonne, Shiv, Nate, Dan, Alan and Salim for their large, small, direct or indirect contributions. Special thanks go to Fred for his constant help-in the data collection phase. The author also would like to thank the Department of Electrical Engineering of the University of Florida, Comissao Permanente INVOTAN (NATOPortugal), and the Institute of International Education for the financial support received throughout his graduate studies. The author sincerely thanks his family for their sup port and confidence, especially his wife for her constant ii 158 TAPE RECORDER Procedure: 1 Tape Recorder to Record Mode 2 Extension Port Switch to Reverse 3 Press Y Key in Terminal (after 10" of Leader) 4 When Finished Switch to Neutral . 28. Record from Microcomputer to Tape Recorder Fig 33 available. This fact per se has the potential to expand the field of application of digital filter signal proces sing. This approach is also less expensive since the machine may be a general purpose computer, and so there is no extra cost involved in the new application. The draw back is the computer's fixed architecture which is not adapted to calculate recursion relations. For a lot of applications the requirements of the implementation can be met with mini-computers and even microcomputers. The principal elements of the requirement set are speed, type of arithmetic, wordlength. Digital algorithms work in digital representations of real world (analog) signals. It is well known (Thomas) that to make a digital representation of an analog signal unique (i.e., the analog signal can be reconstructed again from the digital samples) a maximum frequency must be assumed in the signal spectrum. This is the same as saying that the time domain representation of the signal is restricted to change at less than a prescribed rate, governed by the time/ frequency inverse relationship. The theorem that relates the uniqueness of the digital representation and the sampling frequency is the Nyquist theorem, and it imposes a lower bound on the sampling frequency (sampling frequency equal to two f^). This theorem assumes two impossibilities: first, that the signal spectra are frequency band limited, and second, that the past and future history of the signals are completely available for reconstructing the analog 345 Walter, D. 0., Kado, R. T., Rhodes, J. M., and W. R. Adey. "Electroencephalographic Baselines in Astronaut Candi dates Estimated by Computation and Pattern Recognition Techniques." Aerospace Med., 38:371-376, 1967. Walter, D. 0., Muler, H. F., and R. M. Jell. "Semiauto matic Quantification of Sharpness of EEG Phenomena." IEEE Trans. Biomed. Engng., BME-20:53-55, 1973. Weinberg, L. Network Analysis and Synthesis. McGraw-Hill Book Co., New York, 1962. Weir, B. "The Morphology of the Spike and Wave Complexes." Electroenceph. Clin. Neurophysiol., 19:284-290, 1965. Welch, P. D. "The Use of FFT. for the Estimation of Power Spectra. A Method Based on Time Averaging over Short Modified Periodograms." IEEE Trans. Audio and Elec- troacoust., AU-15:70-73, 1967. Wennberg, A., and L. H. Zetterberg. "Application of a Com puter Based Model for EEG Analysis." Electroenceph. Clin. Neurophysiol., 31:451-468, 1971. Wilder, B. J., Willmore, L. J., Villarreal, H. J., Bruni, J., and R. J. Perchalski. "Valporic AcidAcute Study." Published by Epilepsy Research Foundation of Florida, Inc., 1978. Woody, R. H. "Interjudge Reliability in Clinical Electro encephalography." J. Clin. Psychol., 24:251-256, 1968. Yeo, W. C. "A Signal Detection from Noise Utilizing Zero Crossing Information." Ph.D. Dissertation, University of Florida, 1975. Zetterberg, L. H. "Estimation of Parameters for a Linear Difference Equation with Application to EEG Analysis." Math. Biosci., 5:227-275, 1969. Zetterberg, L. H. "Spike Detection by Computer and by Analog Equipment." In Automation of Clinical Electro- encephalography. Raven Press, New York, pp. 277, 234, 1973. 134 130 wave complex was the only parameter, it would be very difficult to implement the desynchronization of the pat tern. The detection being a two-stage process, it is possible to acknowledge the presence of slow waves even when they do not appear with the prescribed repetition rate, as may happen in the middle of the seizure. This feature enabled the recognition of desynchronizations with out breaking the seizure in two separate events. Period Window for the Spike As it was pointed out earlier, one of the spikes (second) present in the PM pattern is relatively slow and of very high energy. The first spike individualized by (Weir, 1965) will not be detected here. To detect most of the second spikes present in the beginning of the PM parox ysm a time window of 0.09 sec (Ajj) to 0.05 sec (6h) corresponding to 11-20 Hz was necessary. It avoids detec tion in alpha bursts. Amplitude Thresholds Previous work in sleep research has shown the incon sistency of amplitude measures for EEG signal detection. 53 decrease it, scaling can be utilized at expenses of dynamic range, or for some cases an alternate transformation can give smaller gains. The other possibility is to decrease the sampling frequency, making the filter appear wider and therefore lowering the Q. Although for certain areas like speech processing this is not of great practical value, for EEG it is a viable procedure, since the EEG rhythms can be considered not coupled. The penalty paid is higher compu tation time as a digital lowpass filtering shall be per formed before the reduction of the sampling frequency to avoid aliasing. Peled (1976) has a detailed analysis of the procedure including setting of lowpass corner frequency and new sampling frequency. The rationale for using lowpass FIR filters is the following: If decimation is desired (lower sampling fre quency at the output), only the samples which will be used need to be computed. In the case of FIR filters, where the output is only a function of the past values of the input, this can be accomplished easily. For HR filters the present output is a function of past values of the input plus of the output, which means that every point must be computed. The threshold where the FIR filter is prefer able to the IIR is given by the ratio of decimation and the difference in filter order to keep similar attenuation characteristics. An automated procedure to design IIR bandpass Chebyshev filters has been presented. It was shown that FREQUENCY ANALYSIS NODE HUNtiEH 16 MAGNITUDE 60.01599.. ** Ht Ir K a -il 1! I: CHEYSMEV 41H ORDER IMPULSE 1NV. BANDPASS C2. 02. 11 Â¡ C 1, D1 1 SLUM NAVE 54. 014B9. 43. 01379. 42. 01270. . 36. 01160. 30.01050. 24.00940. 18. 00330. 12. 00720. 6. 00610 0. 00500. *** it funs*********#*****###***##**#**. O. 10E-0J O 79E OI O. J6E 02 O. 24E 02 O. 32E 02 O.40E 02 O. 40E 01 O. 12E 02 O.20E 02 0 SHE 02 O. 36E 02 FREQUENCY (IITZ) CO 00 Fig. 9c. Filter Frequency Response (Slow Wave) 110 Referring to Fig. 12, the equation of the straight line is 2~Btl + B+A t t^+t2 t^+t (49) where A = -M sin wt^ and B = M sin wt2* For y=0 (zero crossing) E = Btl~At2 (50) B+A To find the maximum of (50) with respect to t2* the follow ing substitutions of variables were made. 0 = t2~T/2 and a = T/2. So (50) becomes e = 0-a cotg(wa)tg(w0). (51) The maximum is -!Â§ = -14 = 0 -* cos^w0 = aw cotg wa O t o Op fc2 = Â§ + Â£ cos"1 ( |w cotg (2^)) . (52) 108 It is also desirable to build a system which will be repeatable, since the advantages of microcomputers may very well broaden the number and areas of application of the system. The choice of the implementation media imposes a strong bias in the detection methods. It is felt, however, that this does not necessarily mean constraining the detec tor performance, as will be exemplified. It just means that certain possible ways of extracting the information of the basic PM components will not be applied here. For instance, if the frequency information of the basic recruiting rate is desirable, time domain techniques (like zero-crossing, or peak detection) will be preferred to the frequency domain analysis, like FFT. When the character istics of the PM paroxisms were summarized, some highlights about the extraction of the parameters were introduced. Here the complete detection scheme will be analyzed. Sampling Frequency The first building block of the detector is the A/D converter, since the entire detection process will be realized in the digital domain. One of the most crucial parameters to be set is the sampling frequency. For EEG signals the frequency content may be set at 40 Hz even for abnormal EEG, since Smith (1974) found no modification in the spike characteristics used for detection when the data were prefiltered at this frequency. However, when the 100 The slow wave detection that will be used here is adapted from the detection of delta waves in the sleep studies. Basically it consists of a bandpass filter 0.8-6 Hz, 12 dB/oct, and amplitude and period detectors. Spike detection is a very difficult problem. The methods employed so far, besides the computationally expen sive amplitude-duration criteria, rely on the sharpness as the main parameter. They are, however, sensitive to the muscle artifact (see literature survey). As the spike is not a very reliable parameter during the PM paroxism, not a lot of effort was directed towards its detection. One approximate model for the spike is a triangular wave. From the Fourier transform of a triangle, most of the energy can be expected at the fundamental frequency (the even harmon ics are zero, and the third harmonic is 1/9 of the funda mental) From the definition of the spike duration, the fundamental frequency will be between 12 and 20 Hz. Hence, a crude spike detector will be a bandpass filter from 12-25 Hz, 12 dB/oct to extract most of its energy, followed by a zero-crossing analysis of the filtered data to separ ate other EEG activities which will have appreciable energy in the frequency range (e.g., alpha wave bursts). The effect of muscle artifact will be minimized since most of its energy is contained above this region (O'Donnell et al., 1974). The next step is to combine the information of the detectors to reach a pattern recognition algorithm. It is Update mem. pointer 2 descriptor of the EEG parameters. At the present time this is a very hot topic where one can read reports that show the short time stationarity of the EEG for sequences below 20 seconds for a probability of error of 10 percent (Cohen & Sanees, 1977; Kawabata, 1973; Saltzberg, 1972). However, the work of Elul (1967, 1969) and Dumermuth (1968) points out the nonstationary behavior of EEG sequences as short as 4 seconds. Another important characteristic, albeit less addressed, is the Gaussianity of the EEG. The work of Elul (1967, 1969) and McEwen and Andersen (1975)' suggested that the behavior of the EEG followed a Gaussian distribution 32 to 66 percent of the time depending upon the behavioral state of the subject (active mental task to relaxed state, respectively). It also displayed a spatial distribution with closer Gaussian characteristics in the occipital leads. Some of the discrepancies that are reported in the literature can be explained by the different sampling rates used in the A/D conversion, since oversampling may distort the local statistical properties of the digital waveforms (Persson, 1974). Another reason, and maybe more plausible, just demonstrates the variability of the EEG with subject and behavioral state. In practice 'the stationarity requirement can be loos ened to stationarity during the observation period, and if the process departs slightly from Gaussianity, ancillary parameters like skewness and kurtosis could be employed along with spectral analysis. If stationarity holds, but 24 integrated amplitude from a combination of 4 EEG channels to detect PM activity. The system stores the duration and time of occurrence of the seizure. The overall agree ment, consensus versus machine, is 85 percent. It is, how ever, interesting to note that in the paper's extensive discussion, no mention is made on artifacts that may cause potential problems like chewing and body movements, which suggests that the investigators had available rather clean data. The system is reported to be sensitive to slow wave sleep. Carrie and Frost (1977) also described a spike and wave detector. It is composed of a spike detector, using the sharpness criteria (Carrie, 1972a), an EMG detector and an information of the amplitude of the background activity (Carrie, 1972b). One channel is analyzed on line, and the system stores the duration of the paroxism and its time of occurrence. The agreement is very high for seizures greater than 3 sec (85 percent) but drops off to 25 percent for Sz between 1-3 sec. The patients were reported to have well-defined PM paroxisms and were free of medication. From the EEG data available to us it seems that the large amount of high energy artifacts that resemble the PM pat tern when the patients are in an unrestrained environment dictates the use of more powerful pattern recognition algo rithms. Johnson (1978) implemented in a microcomputer a PM detector based on the repetition properties of the pattern, i.e., spikes followed by slow waves. He used two analog 1.5 hour i V I- Fig. 38. Spike Amplitude Versus Time (#12-3) 225 Fig. 48a. Chewing No Detections 149 wave reaches the peak, tests on the valley-to-peak ampli tude and period are made. If any of these conditions fails, control is passed to the beginning of the PEAK pro gram. Otherwise, the time to reach the next zero crossing is monitored (R9). If it is within limits (only a positive threshold is tested), a slow wave was recognized and the values of peak amplitude and 1/2 period (valley-to-peak) are transferred to STATIS. A negative pulse is outputted to D/A II to signal the detection of a slow wave, and con trol is passed to MAIN to check the repetition rate. SPIKE. This program has an independent workspace and is linked through a BLWP using R5, R6 of MAIN. Valid data appear in R3. The first loop looks for a zero crossing - to +. As this program has three basic loops, it has to be indicated to PEAK in which of the three loops is the pro gram control, before jumping back. This was accomplished by loading the value of the PC for the next sample in R6 of MAIN. When the zero crossing is detected, Rl is cleared, and the program looks for the peak of the activity. In the next sample control enters the third loop that looks for the next zero crossing. When found, the program tests the period and amplitude requirement. If met, the amplitude and the value on DELAY (number of points between last slow wave peak and spike peak) are transferred to STATIS; the spike recognition location (DELAY) is cleared; and a posi tive peak is outputted to D/A II. 280 93 since the noise created will appear without attenuation at the output. If the poles have approximately the same Q, the pairing is the more determining factor. Nevertheless, the combination which gives the highest magnification shall be realized first, but it may not be apparent from inspec tion. Supposing that none of the combinations can be imple mented with the available hardware (i.e., A/D and microcom puter) before setting for a change in the hardware, the following should be attempted: Design the filter with another transformation (wide-band), and decrease the sampling frequency (decimation). If all the above fail, then the only thing left is to scale the signals between resonators. The second method is the only one that does not sacrifice output signal to noise ratio. Scaling between resonators is probably the one that most deterio rates the signal to noise ratio, since it is done in a place on the structure where the noise is far from being white. Hence, scaling between resonators means reduce the noise power less than where the noise is wide-band (inter nal nodes of structureJackson, 1970). However, it can be accomplished with only one shift operation. As a concluding remark it is felt that the design pro cedure described is of practical value for EEG signal processing. From the topology point of view the most stringent constraint is the number of multipliers. It is felt that even with the practicability of the design 325 170 System Evaluation The group of patients for the drug study was purposely chosen to be medically refractory to other drugs. They were under heavy medication, and for the most part they had a long history of generalized seizures. It was also common to find patients with mixed seizure types (myoclonic or focal). Since the detector has built in specifically a scheme to detect PM seizures, a selection of the patients was judged appropriate to test the PM model implemented. Another reason for the selection can be found in the effec tiveness of the drug, since most of the patients had a total reduction of clinically significant seizures. So no paroxisms will be found in most of .the post-drug records. As stated in Chapter III, the detection criteria were chosen to apply to the classical PM and PM variants. A group of four patients was selected to be used in the set ting of the detection parameters. The group was formed by patients #3 PM variant #7 classical #12 classical #21 PM variant. Examples of their patterns have been presented in Figs. 11a, 27a, 27b, lid. Portions (30 minutes long) of the tapes containing paroxysms were re-recorded and used in the detector calibration. Needless to say, in this group 18 can be inferred. Isaksson and Wennberg (1976) point out that about 90 percent of the records of 20 s are stationary which agrees with some of the previous observations. At the present time, only Isaksson (1975) and Isaksson and Wennberg (1976) have used this technique. Time Domain Approach If one can relate the degree of quantization to the maturity of a field of science, electroencephalography is a newborn. The lack of quantization has two main conse quences: the inexistence of objective criteria to charac terize the EEG phenomenon that weakens any purely analytical research method (like spectral analysis); the resort to a multitude of soft criteria (i.e., not unique), each based on empirically derived parameter values that try to trans late the highly nonlinear way the electroencephalographer reads the EEG. Visual analysis is a true time domain method of detec tion of structural features like diffuse and localized change in frequency and voltage pattern, changes in the topographic distribution and in the interhemispheric sym metry, sharp and rhythmic activity, paroxysms and unstable or irregular time course. The basic patterns and clinical EEG features are composed on the following step-by-step procedure: the analysis of the period and: amplitude .. forms the wave concept that is associated with the gradient to give the simple grapho-element. These can be integrated 198 The data available to us had unfortunately too many variables to permit conclusions considering the comparison of pre- and post-Depakene administration. The problem is the use of different EEG montages. As can be seen from Table VII, different montages were utilized to record each session. There have been reports of parameter dependency with EEG montage in normal activity (Smith, 1978), and in PM it is well known that the amplitude of the slow wave varies with the electrode location. For this reason it is impossible to couple the increase in the repetition period for the third session of patients #3, 4, and 7 with Depakene administration. However, with a more careful data collection this question can be easily answered with the present system, as long as the false detections are excluded (not to bias results as in the case of patient #5-2.) Another bit of information which can be obtained with the new instrumentation system is the correlations between PM parameters to quantify in more detail the PM phenomena. With the set of data available the following problems should be considered: first, the small data set; second, the existence of a number of uncontrolled variables, such as seizures occurring in different physiologic states (sleeping, awake, medication), utilization of different EEG montages, and all the seizure detections considered (including false positives). Therefore, a nonrigorous 275 TITL PEAK* ** 3jt* *** *5* #*##3:*5}:3Â¡!* * PEAK PROGRAM ******* ***.*;** *:*:Â£#*****#** * MEASURES VALLEY TO PEAK * PERIOD ANO AMPLITUDES IN ' * SLOW WAVE FILTERED DATA. * IT ALSO CHECKS ZERO * CROSSINGS (MAX LIMIT) * THESE PARAMETERS ARE PASSED * TO ST AT IS * NEG PULSE IN D/A II WHEN S.W. * IS RECOGNIZED * * CONTROL IS PASSED TO MAIN <3L) * INPUT,SPIKE AND STATTS (SLWP) * *REGI * * * * * * * STER3 USED RO SIGN OF INPUT RI R. RATE COUNTER R2 INPUT DATA R4 VALLEY TO PEAK COUNTER RS,6 SPIKE LINK R7.Q INPUT LINK R9 ZERO CROSSING COUNTER R 12 AMPLITUDE 20A 0 AORG >23AO OOOA ZPLU S EQU >A FOSO MAIN EQU >F08Q FFF6 ZMI N U EQU >FFF6 003C TPLUS EQU >3C 0014 TMI N U EQU >14 OOOC TM IN I EQU >C 0 01 F TM AX I EQU > 1 F FC 04 THRE S EQU >FC04 07FF ONE E QU > 7FF FC06 MASK EQU > F CO 6 FC 00 DELAY EQU >F CO 0 FC 1 A ST AT 4 EQU >FC 1 A FDCO WSTA2 EQU > F C C 0 FEOO WSTA4 EQU >FE 0 0 FC 1 2 ST AT 2 EQU >FC12 2 2 A 0 02E0 LWPI MAIN 23A2 FD30 22A4 Cl 00 PEAK MOV 0 ,4 2 3 AS 05A0 INC 2DELAY 23AS FCOO 2 2A A C 32 0 MOV DDELAY,00ELAY 23 AC FCOO KEEP S IGN IN R4 INC SPIKE DELAY OVERFLOW(NEG)? 025 sec Fig. 42b. Half Period Versus Duration (#12-3) 216 Microcomputer System Fig. II-l. TI 9900 Microcomputer Block Diagram 252 37 number of operations, providing at the same time good overflow and noise properties. To design a filter with a specific output signal to noise ratio using the minimum requirements of computation wordlength and processing time, the following factors need to be analyzed: 1) Filter transfer function 2) Filter internal magnification 3) Noise characteristics of filter topology 4) Processing speed. The problem as enunciated is relatively different from the design procedure commonly found in the literature, since there the optimization is studied with respect to only one parameter (most generally to sensitivity to round off) neglecting any other considerations. It turns out that the optimum structures possess a much higher number of multipliers which will mean slower computation time in our application. Examples are as follows: The synthesis of in finite impulse response filters (HR) with low roundoff noise based on state space formulations has been presented by Hwang (1976) and Mullis and Roberts (1976). Their struc- 2 tures have N more multipliers than the canonical structures, where N is the filter order. Another example will be the structure proposed by Barnes et al. (1977), which is free from overflow oscillations, but requires N+l more multi pliers than the canonical structures. Multiplier extrac tion was also considered (Szczupack & Mitra, 1975), False Detection tan o * Fig. 49c. Mean and Variance of Repetition Period Versus Time (#12-1) 08 sec 0.639 sec 284 2 40 6 3000 2 438 13 08 JEQ NEG1 2 40 A C SAO MOV QPC2 24 EC FCOA 2 43E OOOC 24C 0 C 083 CUME MG V 3,2 NO,REACH PEAK? 24C2 C380 RTWP 24C4 8083 C 3,2 2 4C6 1 5FC JGT CUME 24C8 1 0E9 JMP DACK 2 4C A C9A0 NEG1 MQV 3PC3 ,SC(6 ) 24CC FCOC 24CE OOOC 24D0 0 4CI CLR 1 CLEAR DURATION COUNTER 2402 Cl 43 NEC MCV 3, 5 2404 0380 RTWP 240 6 C581 I NC 1 2408 8 143 C 3, 5 REACH VALLEY? 2 40 A 1 1FB JLT MEG 2 40 C 60 85 S 5,2 24DE 8802 C 2,3SPTHP BIG ENOUGH? 2 4E C FC 02 2 4E 2 1 ICC JLT BACK 24E4 Cl 00 BAQ MQV 0,4 SAVE SIGN CF INPUT 2 4E6 CO 03 MOV 3 ,0 GET NEW SIGN 24E8 4020 32 C 3MASK 0 2 4EA FC06 2 4EC A987 A 7,a>C(6) PC FOR NEXT JUMP 24EE OOOC 2 4F C 0380 RTWP 2 4F2 0581 INC 1 24F4 8100 C 0,4 NEXT ZERO CROSSING? 2 4F6 13F6 JEQ SAG 24 F8 0281 Cl 1,HRATE WITHIN LIMITS? 2 4F A COOA 24FC 15CF JGT BACK 24FE 0281 Cl 1 LRATE 2500 0005 2 50 2 1 ICC JLT BACK 2504 C820 MOV 3DELAY,3WSTA3 PASS DELAY SPIKE S.W. 2 50 6 FC 00 2508 FDEO 2S0 A C302 MOV 2.3WSTA5 PASS SPIKE AMP TO STAT 2 50C FE2C 250 E 04E0 CLR 2DELAY 2510 FC 00 2512 0420 BLV*P 3STAT3 251 4 F Cl 6 2516 C 420 3LWP 3STAT3 2518 FC IE 251 A 02 04 LI 4,>7FF OUTPUT + 3UL SE TO 0/A 251 C 0 7FF 251 E C804 MOV 4,3>1EF2 2520 1 EF2 2522 1 08C JMP SACK 2524 1000 NOP END NO GF ERRORS IN THIS A 3SE VBLY= 0000 RELOCATABLE LOCATIONS USED = 0000 38 but the structures obtained lack the optimum sensitivity characteristics. Another interesting idea is presented in Chang (1978), where the decrease in the roundoff noise is accomplished by feeding back the discarded bits, properly weighed, to the input of the adder following the place where the product is quantized. This procedure doubles the multiplier number. Szczupahk.and Mitra (1978) propose the reduction of the roundoff noise by ensuring the presence of zeros in the noise transfer functions. This procedure leads to zeros that are complex, therefore requiring extra multipliers for its realization. Many more examples could be given, but the picture remains unchanged; decrease of roundoff noise means more multipliers. For a microcomputer implementation of digital filters multiplication is by far the most time-consuming operation, and it is therefore necessary to minimize its use. It seems appropriate to say that for microcomputer implementation a suboptimal solution would be the best, weighing the effect of the increase in computation time and the decrease in the variance of the output noise. This analysis is not known to exist, so the procedure chosen utilized canonical structures as well as slight variations in their topology (adders as a variable). Filter Transfer Function The characteristics of the filters presently used for sleep EEG studies have been obtained through trial and error. Basically the filters are low Q, second order 157 3 FC2A, FC2C, FC2E. Testing of the Program The efforts undertaken in the program modularity were rewarded in the testing phase. A brief description of the facilities available is thought helpful at this point. To write the programs in assembly language the cross assembler resident in the Amdahl 470-V6 was extensively used. The editing facilities of the operating system and the error messages of the cross assembler were thought adequate for the software develop ment. The same feeling is shared with the cross assembler resident in the NOVA 2/10. The program loading in the TI microcomputer from any of the computers was performed auto matically through the resident TIBUG monitor (LOAD command) via a Modem (Pennywhistle 130). To get a permanent copy of the programs after being loaded in RAM memory a REALISTIC cassette recorder was used, along with the DUMP command of the TIBUG monitor, and a Modem. Fig. 29 shows the respec tive connections. Care must be taken in the volume and tone setting in the cassette recorder. There were no relo catable facilities in the monitor; therefore, only absolute origins were used in the assembly programming. The debugging of the program functions was not as easy, first, due to the complexity and dependence of the program flow in the input data, and second, due to the poor 240 w Passband Trans. Attenuation Band Band Fig. I-I. Chebyshev Low Pass Characteristic 96 discharges, that is, paroxysms with duration greater than 3 seconds. In PM epilepsy the paroxysms are spike and wave com plexes. A complex is a group of two or more waves clearly distinguished from background activity and recurring with a well-organized form. The spike has a duration of 1/12 of a second or less, and a wave has a duration of 1/5 to 1/2 second. The spike is usually biphasic and has a preva lently negative polarity. The amplitude is generally greater than that of background activity and averages 50-150 yV. The wave is a sinusoidal-type waveform with medium to high negative voltage that averages 100 yV. The spike can be less apparent or be substituted by polyspikes. In PM epilepsy the wave complexes appear predominantly in the frontal leads and are bilaterally synchronous and sym metrical, recurring rhythmically at about 3 c/sec (faster in the beginning of the seizure, 3-4 Hz, and slowing down toward the end, 3-2 Hz). Those are the definitions avail able in the clinical neurophysiology literature (Remond, 1974). Although not exactly qualitative the definitions contain a lot of empirical notions as "clearly distin guished from background activity," "recognized form," "prevalently negative polarity," "amplitude generally greater," "about 3 c/s." There are aspects where there are no definitions, as the morphology of the spike and the interval between the spike and the wave. There are aspects open to the personal interpretation of the EEG scorer. CHAPTER IV SYSTEM EVALUATION AND PRESENTATION OF RESULTS Description of Data Collection The data for the PM study was collected in the Veterans Administration Hospital, Neurology Service, from an ongoing drug study. At the time, Dr. B. J. Wilder was evaluating the clinical efficacy and safety of valporic acid (Depakene) in patients with uncontrolled absence sei zures . Twenty-five patients (14 men and 11 women) with absence seizures were the subjects of the study. Almost all the patients (23) were receiving other anticonvulsant drugs. The central part of the study spanned a period of 4 months. One EEG tracing was obtained in the pre-entry period; then a two-week single blind placebo period was followed by the second EEG. Ten weeks of therapy with valporic acid then followed. Another EEG tracing was obtained at the end of this period. The EEG's were for the most part obtained with a 4 channel telemetry system (Datel Model 1000, Physiological Telemetry Systems), and the patients were free to move about in a 3x5 meter room. There were different monopolar montages used, but the ones most frequently utilized were 165 0 TABLE IV COMPUTER/HUMAN AGREEMENT FOR INDIVIDUAL PATIENTS PATIENT #3 Human 1 < Sz < 3 sec 3 < Sz < 10 sec Sz > 10 sec Computer Agreement in Seizure Detection (# of Sz) 399 221 74 90 31 ^ 31 Agreement in Seizure Duration (# of Sz) 200 63 24 Total Time in Agreement (sec) 416 . 421 514 576 Computer Counts 245 92 26 Computer Misses 178 16 0 Computer False Detects 37 0 0 176 TABLE VIIcontinued PATIENT #5 193 820 sec Fig. 46. Variance Half Period Versus Half Period (#12-1) 222 289 143 Fig. 25. PM Detector Flow Chart. Dotted Line Mean Data Dependent Path. 144 entry and output. To turn the communications between these programs bimodal, they would have to be broken down further without showing any relevant advantage, except possibly easier testing. On the other hand, this would have increased the computation time. Another characteristic that is shown in the diagram is the dependence of the com munications on the input data (dotted lines). For instance, the closed loop in PEAK means that control is in the module up to the point where a slow wave is detected. A typical "closed loop" description on the algorithm goes as follows: After initialization (MAIN) the program analyzes serially the occurrence of a spike (SPIKE) and of a slow wave (PEAK) in the filtered data (FILTER). The con trol remains in PEAK up to the point where a slow wave is detected. For every input point, the program also tests for the occurrence of a spike. The communication between SPIKE and PEAK depends on the particular phase of waveform recognition. When three slow waves with the prescribed repetition rate and at least one spike are detected, the beginning of a seizure is acknowledge (MAIN); flags are set at the output; and the time of occurrence is stored in memory (SZURE). The program keeps looking at more slow waves (PEAK) and spikes (SPIKE), storing the amplitude and period for each event (STATIS). The periodicity of the slow waves is checked (SZURE), and when no slow waves are detected for the last 1 sec, seizure count is halted; sta tistics of the detection parameters are obtained (STATIS) 72 Fig. 6. Additive Model for Quantization 19 in a sequence to produce the complex grapho-elements from which, using superposition, the notion of pattern is acquired. The variability in time and hemispheric location gives rise to the time domain structure and topography of the EEG. The job of the electroencephalographer is to individualize these elements and compare them to his acquired set of "normal activity" in order to make the diagnosis. Of course, the boundaries are quite fuzzy, highly subjective, and not always consistent (Woody, 1968; Rose, 1973). . The techniques are then dependent upon the specific detection task desired, and no simple overview is possible besides the translation to electronic terms of the decom position process: analysis of the wave's amplitude in cer tain frequency ranges, which implies broad bandpass filter ing, zero crossing, and threshold analysis. Depending upon the criteria formulated the next step can be the testing of the grapho-elements throughout further processing (sharp ness or repetition period, for instance) or testing of a pattern by coincidence logic. What makes the process some what erratic is that there is no methodology beyond the extensive comparison with the results obtained by the electroencephalographer and accordingly modify the imple mentation or the criteria until an agreement is reached. Due to the diversity of techniques, the detection of abnormal brain activity (spikes) will be emphasized. One of the first criteria to detect spikes was the sharpness 43 TABLE I REQUIRED NUMBER OF ZEROS t3 n i* Z = -1 Direct substitution narrow-band 1 1 Direct substitution wide-band 2 1 Bilinear 2 2 27 expensive. The inverse filtering can not differentiate among types of nonstationarities. Therefore, a preproces sing or other ancillary methodology needs to be developed. Kalman filtering is quite sensitive to the gain setting. For high gains small changes can be exaggerated, and so one is faced with the problem of calibrating the filter, which is usually done with a control data set (artificially gen erated) The direct dependence of the gain on the type of data, the dependence of the error signal on the power of the input sequence, and some transient behavior of the system are shortcomings. It seems that the popularity of frequency domain tech niques in EEG detection stems from the availability of standard software packages and a fairly well developed (but often forgotten) mathematical methodology. Up to the pre sent, one of the drawbacks of the time domain approach is the implementation medium (hardware) which requires a full engineering development, not always accessible to the clinicians. However, the present innovation trend in microcomputer systems may very well bring the possibility of standardization through the software implementation of detectors, at low price. To build an electronic detector where some parameters (amplitude in a prescribed frequency range, zero crossing, sharpness) have to be monitored and decisions made, the omissions created by the poor quantization (definitions) of the EEG process have to be filled. That is the reason why 287 TITL ST A TIS' **** *** ******4e*** fc****^**.* * * STATIS $Â£:#* ifc*#:*:^*-*:* i*:* 5**-#*:*:$:***Â£* * COMPUTES RUNNING AVERAGES. * C3TAINS MEAN VALUES AND * VARIANCES. * MEAN R1 * VARIANCE R 4 * DATA IS FORMATTED TO USE * ALWAYS MOST X OF 3 ITS. CODE (REEFERS TO .MOST SIG. DIGIT) * 0 TRUE MAGNITUDE * 4 DATA SHIFT L. 1 HEX DIG. $ 0 H 2 ti * C ONLY VARIANCE SHIFT R. 1 HEX DIG. REGISTERS USED * RO INPUT DATA * R 1 MOST SIG. WORD OF MEAN R2 LEAST * ' It It R3 X OF ENTRIES * P 4 MOST SIG. WORD OF VARIANCE * R5 LEAST P6 SCRATCH 2590 AORG >2590 FC 22 FLAG EQU >FC22 FC26 MASK 4 EQU >FC26 FC2S MASKS EQU >FC2S 2590 CS20 MOV 2FLAG,3FLA G TOTALS? 2 59 2 FC22 2594 FC22 2596 1600 JNE FINI 2598 A080 A 0.2 NO,FORM PUNNING AVE. 2 59 A 1701 JMC SMAL1 OF MEAN 2 59 C 0581 INC 1 DOUBLE PRECISION 2 59E 0 583 SMAL 1 INC 3 X OF ENTRIES 25A 0 C 180 MOV 0.6 GET X*X 25A2 3980 MPY 0,6 2 5A4 A1 06 A 6 .4 RUNNING AVE FOR VARIANCE 25A6 A l 47 A 7,5 25A3 1701 JNC SMAL2 2 5A A 0584 INC 4 25 AC 0380 SMAL2 RTWP * Â£ TOTALS 25AE 3C43 FI N I D I V 3,1 R1 HAS MEAN 2500 3 C 03 D I V 3,4 R4 HAS 1/N*SUM(X*X) 2 582 1910 JNO SMALL DIVISION POSSIBLE? CHAPTER III A MODEL FOR PETIT MAL SEIZURES AND ITS IMPLEMENTATION Detection Problem The automated detection of pathological paroxysms in the EEG is at the present time an unsolved problem. The main reasons are the variability of the wave patterns and the nonexistence of quantitative definitions of pathologi cal paroxysms and of their elements (spikes and slow waves). A paroxysm is a group of waves which appears and disappears abruptly and which is clearly distinguished from background activity. When the paroxysm is not related to normal brain function, it is called pathological (epilepti form) EEGers differentiate between interictal epileptic discharges (events that are not accompanied by clinical seizures) and ictal epileptic discharges (when clinical seizures are evident). The definitions are in general hard to apply, but for the special case of petit mal (PM) epi lepsy the difference is solely based on the duration. It was found that spike and wave bursts with duration shorter than 3 seconds did not affect appreciably the behavior of the patient (Sellden, 1971). Throughout this work the evaluation of the detector will be restricted to the ictal 95 APPENDIX I CHEBYSHEV FILTER DESIGN Chebyshev Polynomials The general formula for the Chebyshev polynomials of order n is (Weinberg, 1962) T (w) = cos (n arcos w) (1-1) n which approximates a low pass characteristic with a magni tude that does not exceed a prescribed maximum deviation in the passband and displays the fastest possible rate of cut off outside the passband (Fig.1-1). Equation (1) can be transformed into a more recognizable polynomial form by letting , -1 (J) = cos w and substituting again cos c|>=w and sin cj>=/ 1-w^, giving T (w) = cos ncf> = R (cos
n PAGE 68 )LJ D 2YHUIORZ &KDUDFWHULVWLFV RI 7ZRnV &RPSOHPHQW )LJ E 1R 2YHUIORZ 3URSDJDWLRQ 5HJLRQ PAGE 69 _&__'_ f ZKLFK LV VKRZQ LQ )LJ E DV WKH VTXDUH KDWFKHG 7KLV FRQn GLWLRQ LV YHU\ UHVWULFWLYH IRU PRVW DSSOLFDWLRQV WKHUHIRUH PRVW RI WKH SUDFWLFDO GHVLJQV ZLOO GLVSOD\ WKH XQGHVLUDEOH HIIHFW RI RYHUIORZ SURSDJDWLRQ 6LQFH WKH RYHUIORZ LV D QRQOLQHDULW\ LQ WKH V\VWHP LW FDQ EH H[SHFWHG WKDW YDULRXV LPSOHPHQWDWLRQV GLIIHUHQW ZD\V RI FDOFXODWLQJ WKH UHFXUVLRQ UHODWLRQf ZLOO GLVSOD\ GLIIHUHQW RYHUIORZ SURSHUWLHV 7KH SDUDPHWHUV RQH FDQ FKRRVH WR FRQWURO WKH RYHUIORZ DUH VSHFLDO W\SH RI DGGHUV VFDOLQJ DQG H[SORULQJ WRSRORJLF GLIIHUHQFHV LQ WKH VWUXFn WXUHV :KHQ WKH UHVXOW RI DQ DGGLWLRQ LV ELJJHU WKDQ WKH UHJLVWHU OHQJWK LQ WZRnV FRPSOHPHQW UHSUHVHQWDWLRQ WKHUH LV DQ DEUXSW FKDQJH LQ VLJQ ,W KDV EHHQ VKRZQ )HWWZHLV t 0HHUNLWWHU f WKDW DW OHDVW IRU WKH FODVV RI ZDYH GLJLWDO ILOWHUV VDWXUDWLRQ DULWKPHWLF DOORZV IRU WKH DEVHQFH RI SDUDVLWLF RVFLOODWLRQV +HUH VDWXUDWLRQ DULWKn PHWLF PHDQV WKH W\SH WKDW KROGV WKH VLJQDO DW PD[LPXP OHYHO ZKHQ DQ RYHUIORZ RFFXUV 7KLV LPSOLHV WKDW RYHUIORZ PXVW EH PRQLWRUHG 7KH UHODWLRQV EHWZHHQ WKH W\SH RI VDWXUDWLRQ DULWKPHWLF GLVWRUWLRQV DW WKH RXWSXW DQG VWDELOLW\ DUH D SUHVHQW DUHD RI UHVHDUFK &ODDVHQ HW DO f 1R JHQHUDO UHVXOWV DUH DYDLODEOH 6FDOLQJ WKH VLJQDO OHYHO LV WKH PRVW JHQHUDO ZD\ WR KDQGOH RYHUIORZ 7R SUHYHQW DGGHU RYHUIORZ WKH LQSXW [Q PAGE 70 PXVW EH VFDOHG VR WKDW WKH RXWSXW \QO IRU DOO Q 7KHUH DUH D QXPEHU RI DSSURDFKHV WKDW FDQ EH FRQVLGHUHG WR DFFRPn SOLVK WKLV 2SSHQKHLP t 6KDIHU f $V , KNrQN ! N ZKHUH KA LV WKH ILOWHU UHVSRQVH LW LV HDV\ WR VHH WKDW ,\r O ,KN1[QN, r PD[_[N_f e _KN_ f N N 7KHUHIRUH E\ PDNLQJ PD[_[Y_ e _K _ RYHUIORZ RI \ FDQ . N . Q EH SUHYHQWHG 7KLV PHDQV RQ WKH RWKHU KDQG WKDW WKH ]HUR IRUPLQJ SDWKV RI WKH ILOWHU PXVW EH PXOWLSOLHG E\ &2 N e _K _ ([SHULPHQWV KDYH VKRZQ WKDW WKLV VFDOLQJ N N SROLF\ LV WRR SHVVLPLVWLF DQG GRHV QRW DOORZ JRRG XVH RI WKH IXOO G\QDPLF UDQJH RI WKH ZRUGOHQJWK $ VHFRQG DSSURDFK ZRXOG EH WR DVVXPH WKH LQSXW VLQXVn RLGDO ZLWK D IUHTXHQF\ HTXDO WR WKH UHVRQDQW IUHTXHQF\HVf RI WKH ILOWHU 7KLV LV HTXLYDOHQW WR UHTXLULQJ + f LV WKH ILOWHUnV IUHTXHQF\ UHVSRQVHf _+HA;f _ 77 ; LU f DQG LV DFFRPSOLVKHG E\ PXOWLSO\LQJ WKH FRHIILFLHQWV E\ N WKH VFDOLQJ FRQVWDQW VR WKDW PAGE 71 PD[ N M; aM ; DJDAH DH M; M ; EAH W!p f ZKHUH DJ DA D EA E DUH WKH FRHIILFLHQWV RI WKH VHFRQG RUGHU ILOWHU 7KLV VFDOLQJ SROLF\ LV VRPHWLPHV WRR RSWLn PLVWLF WKDW LV WR VD\ LW PD\ OHDG WR RYHUIORZ $ WKLUG DSSURDFK XVHV WKH IUHTXHQF\ GRPDLQ UHODWLRQ <]f +]f;]f 7KH HQHUJ\ RI WKH RXWSXW VLJQDO LV e \ ZKLFK E\ N . 3DUVHYDOV UHODWLRQ PHDQV WW Q L[ -R\N W |