Open Collections

UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Single trial EEG signal analysis using outlier information Birch, Gary Edward 1988

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata

Download

Media
831-UBC_1988_A1 B57.pdf [ 4.79MB ]
Metadata
JSON: 831-1.0097955.json
JSON-LD: 831-1.0097955-ld.json
RDF/XML (Pretty): 831-1.0097955-rdf.xml
RDF/JSON: 831-1.0097955-rdf.json
Turtle: 831-1.0097955-turtle.txt
N-Triples: 831-1.0097955-rdf-ntriples.txt
Original Record: 831-1.0097955-source.json
Full Text
831-1.0097955-fulltext.txt
Citation
831-1.0097955.ris

Full Text

SINGLE TRIAL EEG SIGNAL ANALYSIS USING OUTLIER INFORMATION By GARY EDWARD BIRCH .A.Sc., The U n i v e r s i t y of B r i t i s h Columbia, 1983 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR*OF PHILOSOPHY i n THE FACULTY OF GRADUATE STUDIES (Department of E l e c t r i c a l Engineering) We accept t h i s t hesis as conforming to the required standard THE UNIVERSITY OF BRITISH COLUMBIA June, 1988 © Gary Edward Birch, 1988 In presenting this thesis in partial fulfilment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Department The University of British Columbia 1956 Main Mall Vancouver, Canada V6T 1Y3 DE-6(3/81) ABSTRACT The goal of t h i s thesis work was to study the c h a r a c t e r i s t i c s of the EEG s i g n a l and then, based on the i n s i g h t s gained from these studies, pursue an i n i t i a l i n v e s t i g a t i o n i n t o a processing method that would extract useful event r e l a t e d information from s i n g l e t r i a l EEG. The fundamental t o o l used to study the EEG s i g n a l c h a r a c t e r i s t i c s was autoregressive modeling. E a r l y i n v e s t i g a t i o n s pointed to the need to employ robust techniques i n both model parameter estimation and s i g n a l estimation a p p l i c a t i o n s . Pursuing robust techniques u l t i m a t e l y led to the development of a s i n g l e t r i a l processing method which was based on a simple neurological model that assumed an ad d i t i v e o u t l i e r nature of event r e l a t e d p o t e n t i a l s to the ongoing EEG process. When event r e l a t e d p o t e n t i a l s , such as motor r e l a t e d p o t e n t i a l s , are generated by a unique a d d i t i o n a l process they are "added" in t o the ongoing process and hence, w i l l appear as a d d i t i v e o u t l i e r content when considered from the point of view of the ongoing process. By modeling the EEG with AR models with robustly estimated (GM-estimates) parameters and by using those models i n a robust s i g n a l estimator, a "cleaned" EEG s i g n a l i s obtained. The o u t l i e r content, data that i s extracted from the EEG during cleaning, i s then processed to y i e l d event r e l a t e d information. The EEG from four subjects formed the basis of the i n i t i a l i n v e s t i g a t i o n i n t o the v i a b i l i t y of t h i s s i n g l e t r i a l processing scheme. The EEG was c o l l e c t e d under two conditions: an aictive task i n which subjects performed a i i s k i l l e d thumb movement and an i d l e task i n which subjects remained a l e r t but di d not carry out any motor a c t i v i t y . The o u t l i e r content was processed which provided si n g l e t r i a l o u t l i e r waveforms. In the a c t i v e case these waveforms possessed consistent features which were found to be r e l a t e d to events i n the i n d i v i d u a l thumb movements. In the i d l e case the waveforms did not contain consistent features. Bayesian c l a s s i f i c a t i o n of a c t i v e t r i a l s versus i d l e t r i a l s was c a r r i e d out using a cost s t a t i s t i c r e s u l t i n g from the a p p l i c a t i o n of dynamic time warping to the o u t l i e r waveforms. Across the four subjects, when the d e c i s i o n boundary was set with the cost of m i s c l a s s i f i c a t i o n equal, 93% of the a c t i v e t r i a l s were c l a s s i f i e d c o r r e c t l y and 18% of the i d l e t r i a l s were i n c o r r e c t l y c l a s s i f i e d as a c t i v e . When the cost of m i s c l a s s i f y i n g an i d l e t r i a l was set to be f i v e times greater, 80% of the a c t i v e t r i a l s were c l a s s i f i e d c o r r e c t l y and only 1.7% of the i d l e t r i a l s were i n c o r r e c t l y c l a s s i f i e d as a c t i v e . TABLE OF CONTENTS Page ABSTRACT i i TABLE OF CONTENTS i v LIST OF TABLES v i LIST OF FIGURES . . v i i ACKNOWLEDGEMENTS i x CHAPTER 1: INTRODUCTION 1 1.1 Problem Statement 1 1.2 Motor P o t e n t i a l s 3 CHAPTER 2: MODELING THE EEG SIGNAL 7 2.1 Need for S t a t i s t i c a l Understanding of EEG 7 2.2 Previous Stochastic Studies on EEG 8 2.3 Neural Basis for the Gaussian Nature of EEG 10 2.A Applying AR Modeling to EEG 13 CHAPTER 3: THEORY OF AR MODELING 17 3.1 Conventional AR Parameter Estimates 17 3.2 LSQ Parameter Estimation on Simulated Gaussian Data 24 3.3 Deviations from Gaussianity 35 3.4 GM Estimation on Simulated Contaminated Gaussian Data . . . . 36 i v CHAPTER A: OUTLIER PROCESSING OF SINGLE TRIAL EEG Al A . l Experimental Design and EEG Data A c q u i s i t i o n Al A.2 Neurological Premise . A6 A.3 Signal Cleaning Process A7 A.A Extracting and Processing O u t l i e r Information 51 A.5 AR Spectral Analysis 60 A.6 Applying O u t l i e r Processing to Single T r i a l EEG 67 A.7 S t a t i s t i c a l Analysis of Features i n the O u t l i e r Patterns . . . 78 A.8 A p p l i c a t i o n of Dynamic Time Warping to O u t l i e r Patterns . . . 87 CHAPTER 5: CONCLUSION 100 5.1 Summary of Major Results and Related Conclusions 100 5.2 Areas for Future Investigation 10A 5.3 S i g n i f i c a n t Contributions 106 REFERENCES 108 v LIST OF TABLES Page 4.1 Feature S t a t i s t i c s from Single T r i a l O u t l i e r Patterns 86 4.2 Difference Between Idle and Active Warping Costs 93 4.3 Group S t a t i s t i c s 95 4.4 Bayesian C l a s s i f i c a t i o n 96 4.5 Summary of Bayesian Cross-Matched C l a s s i f i c a t i o n 99 v i LIST OF FIGURES Page 1.1 Primary Motor and Somatosensory Areas of the Cerebral Cortex . . . 4 3.1 Least Squares AR Parameter Estimation 26 3.2 AR Parameter Estimation on Gaussian Data 38 3.3 AR Parameter Estimation on Gaussian Data with 10% AO Contamination 39 3.4 AR Parameter Estimation on Gaussian Data with 20% AO Contamination 40 4.1 Experimental Apparatus and Conventional Grand Averages from Previous Motor P o t e n t i a l Study . . . . . 42 4.2 Robust AR Parameter Estimation 50 4.3 Process to Extract O u t l i e r Points 52 4.4 Correlated O u t l i e r Detection 54 4.5 AR Spectral Estimates of One Second EEG Segments Consecutively Offset by a Third of a Second 61 4.6 Progression of AR Spectral Estimates with Increasing Order using Idle EEG Data 63 4.7 Progression of AR Spectral Estimates with Increasing Order using Active EEG Data 64 4.8 AR Spectral Estimation of Idle Task EEG 66 4.9 Average of Segmented Orginal and Cleaned EEG 68 4.10 Example Waveforms 71 v i i A.11 Average Waveforms 75 A.12 Comparison of N=15 Averaged O u t l i e r Patterns 79 A.13 Single T r i a l GM2 O u t l i e r Patterns with Corresponding Thumb Movements 80 A.lA Single T r i a l Feature D e f i n i t i o n 8A A. 15 Standard O u t l i e r Patterns . . . .* 89 v i i i ACKNOWLEDGEMENTS I would l i k e to thank everyone i n the Psychophysiology Laboratory at UBC who have helped to f a c i l i t a t e t h i s work. In p a r t i c u l a r , I would l i k e to thank Adelle Forth and Timothy Harpur f o r t h e i r invaluable assistance i n gathering the EEG data. I would l i k e to express my deep gratitude to Dr. J . C. Lind, who introduced me to robust s t a t i s t i c s and who reminded me that i n some cases o u t l i e r s are l i k e diamonds i n a sea of semi-precious stones. I am g r a t e f u l to my committee members, Drs. M. P. Beddoes, G. Dumont, and R. D. Hare, for t h e i r support and expert advice. I would l i k e to give s p e c i a l a d d i t i o n a l thanks to Dr. Hare for making h i s laboratory f a c i l i t i e s a v a i l a b l e to me and for h i s strong support of ray work since the beginning. I extend my deepest appreciation to my supervisor Dr. P.D. Lawrence for his guidance, encouragement and dedicated support throughout my graduate career. To my family, I am thankful f o r t h e i r encouragement throughout my education. I would l i k e to e s p e c i a l l y thank Brenda for her patient and nurturing support. i x 1 CHAPTER 1 INTRODUCTION 1.1 Problem Statement Most often the greatest f a i l i n g of t e chnical aids for severely d i s -abled persons i s the inadequacy of the man/machine i n t e r f a c e . With a univer-s a l e f f e c t i v e and e f f i c i e n t i n t e r f a c e , current technology has the c a p a b i l i t y of providing substantial independence and hence, a g r e a t l y improved q u a l i t y of l i f e for even the most severely disabled persons. In pursuit of such an i d e a l i n t e r f a c e , researchers have been studying the f e a s i b i l i t y of u t i l i z i n g e l e c t r i c a l b r a i n p o t e n t i a l s to d i r e c t l y communicate to peripheral devices. Two important example applications would be the supervisory control of a robotic arm and the method of input to a personal computer system. Such an i n t e r f a c i n g c a p a b i l i t y would also prove to be very useful i n many man/machine int e r f a c e problems i n the able-bodied population. The ultimate goal of t h i s researcher i s to provide a d i r e c t communication system between man and external devices using e l e c t r i c a l b r a i n a c t i v i t y . Brain p o t e n t i a l s are comprised of continuous random e l e c t r i c a l a c t i v i t y which when recorded are r e f e r r e d to as electroencephalographic (EEG) s i g n a l s . Embedded within the EEG <are event r e l a t e d p o t e n t i a l s (ERP) and t y p i c a l l y , there i s a -6 db to a -9 db signal-to-noise r a t i o between the ERP and the EEG. This a c t i v i t y i s u s u a l l y recorded using surface electrodes on the scalp, since the current r i s k versus benefit s i t u a t i o n does not j u s t i f y electrode implantation. The EEG signals as measured from the scalp surface are i n the order of 5-50 micro-volts and are e a s i l y contaminated by other 2 b i o - e l e c t r i c a l signals such as electroocular (EOG) p o t e n t i a l s and muscle (EMG) p o t e n t i a l s , e s p e c i a l l y those of the scalp and face. To date, the study of ERP's has mostly been confined to the averaging of EEG recorded during a s p e c i f i c event, such as a f l a s h of l i g h t , over many t r i a l s . With averaging, the random s i g n a l not r e l a t e d to the event tends to average out to zero at a rate generally proportional to the square root of the number of t r i a l s averaged. The s i g n a l r e l a t e d to the event i s assumed to remain constant and hence becomes more recognizable as the background random s i g n a l decreases. This approach has many drawbacks. In the context of t h i s work, two of the greatest drawbacks are f i r s t l y , that i t i s not amenable to real-time process-ing of event r e l a t e d p o t e n t i a l s for use i n closed loop control applications and, secondly, there i s a s i g n i f i c a n t loss of unique si n g l e t r i a l informa-t i o n . However, i t i s a useful t o o l for obtaining a general idea of the underlying waveform of an event r e l a t e d response. One of the most s i g n i f i c a n t obstacles that must be overcome i n pursuing the ultimate goal i s the establishment of a s i g n a l processing method that can extract event r e l a t e d information from si n g l e t r i a l EEG. There have been some sing l e t r i a l processing schemes, proposed by various researchers (see for example [1,2,3]), that were designed to detect features i n the EEG that are r e l a t e d to s p e c i f i c external events. The usefulness of t h e i r r e s u l t s has been generally l i m i t e d because t h e i r schemes have often been p a r t l y dependent on fundamental assumptions about the s t a t i s t i c a l character-i s t i c s of the EEG which at the present time are not well understood and, more importantly, have been c r i t i c a l l y dependent on complex external v i s u a l s t i m u l i . The goal of t h i s t hesis work was to study the c h a r a c t e r i s t i c s of the 3 EEG s i g n a l and then, based on the i n s i g h t s gained from these studies, pursue an i n i t i a l i n v e s t i g a t i o n into a processing method that would extract useful event r e l a t e d information from si n g l e t r i a l EEG. The fundamental t o o l used to study the EEG s i g n a l c h a r a c t e r i s t i c s was autoregressive modeling. E a r l y investigations pointed to the need 1so employ robust techniques i n both model parameter estimation and s i g n a l estimation a p p l i c a t i o n s . Pursuing robust techniques u l t i m a t e l y lead to the development of a s i n g l e t r i a l processing method which was based on a simple neurological model that assumed an addi-t i v e o u t l i e r nature of event r e l a t e d p o t e n t i a l s to the ongoing EEG process. The EEG from four subjects formed the basis of the i n i t i a l i n v e s t i g a t i o n into the v i a b i l i t y of t h i s s i n g l e t r i a l processing scheme. The EEG was c o l l e c t e d under two conditions: an active task i n which subjects performed a s k i l l e d thumb movement and an i d l e task i n which subjects remained a l e r t but d i d not carry out any motor a c t i v i t y . 1.2 Motor Potentials Brain signals that are r e l a t e d to movement are a type of event r e l a t e d p o t e n t i a l and were f i r s t reported by Kornhuber and Deecke [4]. These poten-t i a l s are produced by the sensory-motor cortex p r i o r to and during voluntary movements of the body. Figure 1.1 [5] shows a cross-section of the motor cortex which i l l u s t r a t e s how various parts of the body are f u n c t i o n a l l y mapped onto the motor cortex. The use of ERP's r e l a t e d to movement have many advantages over other ERP's i n the context of an i n t e r f a c e system. Motor p o t e n t i a l s are produced as a r e s u l t of a s e l f - i n i t i a t e d cognitive process. This i s unlike the 4 Primary Motor and Somatosensory Areas of the Cerebral Cortex Primary motor P r i m a r y somatosensory Figure 1.1 5 perhaps better known subset of ERP's, c a l l e d evoked p o t e n t i a l s , where an external stimulus i s required to e l i c i t a response, (for example: v i s u a l , somatosensory or auditory evoked potentials.) Motor p o t e n t i a l s that are associated with parts of the body over which a disabled person does not have cognitive c o n t r o l , could be u t i l i z e d to obtain unique control of a given peripheral device and thereby minimizing the l i k e l i h o o d of unintentional a c t i v a t i o n of that device. As suggested by Figure 1.1, there i s i n fact some c a p a b i l i t y to discern what part of the body i s being moved by u t i l i z i n g s p a t i a l information from across the motor cortex with the use of surface electrodes. Some of the recent work i n t h i s area has been c a r r i e d out by Brunia and Van den Bosch [6] i n which they e x p l o i t the i p s i l a t e r a l and c o n t r a l a t e r a l properties of the motor p o t e n t i a l s to demonstrate an a b i l i t y to discriminate between hand and foot movements as well as r i g h t versus l e f t body movements. The extent to which t h i s d i s c r i m i n a t i o n c a p a b i l i t y can be r e f i n e d needs to be pursed further but at the very l e a s t i t can provide some d i v e r s i t y i n the control functions that could be derived from the motor p o t e n t i a l s . It has also been demonstrated [7] that, i n the averaged motor p o t e n t i a l s , there i s a sustained response ( i n the order of one second) throughout a prolonged task requiring substantial cognitive involvement. This property could be p a r t i c u l a r l y exploited i n an i n t e r f a c e system that requires a continuous control function such as i n the task of guiding a robotic arm. A s p e c i f i c issue that w i l l need to be addressed i n future work, i n terms of an i n t e r f a c e systems for the disabled, i s the r o l e of peripheral aff e r e n t feedback i n the generation of motor responses. Due to the e a r l y work of Vaughan [8] i t was thought that the p o s i t i v e components a f t e r onset 6 of the movement were due to afferent feedback. However, Vaughan et a l . [9] four years l a t e r unexpectedly found a similar p o s i t i v i t y i n deafferented monkeys while they where carrying out a s e l f paced task. Papakostopoulos [10] postulated that as long as there are elements of s k i l l required i n the task, the p o s i t i v i t y w i l l be developed despite the absence of afferent feed-back. The uncertainty surrounding t h i s issue remains i n r e l a t i v e l y recent work (see Grunewald et a l . [7]) and to resolve i t completely w i l l require the a b i l i t y to analyze the motor potentials on a single t r i a l basis because a disabled person can not provide a movement trigger for conventional averaging techniques. The single t r i a l analysis method described i n th i s thesis provides t h i s required a b i l i t y and therefore, w i l l ultimately f a c i l i t a t e the study of motor potentials from persons who lack peripheral afferent feedback. 7 CHAPTER 2 MODELING THE EEG SIGNAL 2.1 Need for S t a t i s t i c a l Understanding of EEG A pr e r e q u i s i t e to the mathematical modeling of a given s i g n a l i s an adequate understanding of i t s fundamental c h a r a c t e r i s t i c s . Previous i n v e s t i -gations (see Persson [11] or McEwen and Anderson [12]) have noted that the random character of EEG makes the theory of random processes applicable to EEG s i g n a l s . Therefore, i f the approach to the si g n a l analysis of EEG i s to be based on random process s i g n a l theory then basic s t a t i s t i c a l character-i s t i c s of the EEG si g n a l should be well understood. S t a t i s t i c a l properties, p a r t i c u l a r l y assumptions about Gaussianity, are often key factors i n the r e s u l t i n g performance of many of the si g n a l processing methods that have been conventionally applied to EEG. For instance, applications of Wiener f i l t e r i n g , Kalman f i l t e r i n g and AR parameter estimation to EEG have met with mixed and inconsistent r e s u l t s (see McGillem et a l . [13]). This could very well be explained i f i n fac t the si g n a l c h a r a c t e r i s t i c s were at d i f f e r e n t times ranging approximations to the required Gaussian assumptions. On occasions when the si g n a l processing was c a r r i e d out while the EEG was r e l a t i v e l y close to Gaussian, the r e s u l t would have been close to optimal and r e l a t i v e l y good performance would be expected. However, as w i l l be demonstrated i n Section 3.4, i f estimation had been c a r r i e d out on EEG that was r e l a t i v e l y non-Gaussian the performance would have l i k e l y been very poor. Hence, to u t i l i z e s i g n a l processing methods which make c e r t a i n s t a t i s -8 t i c a l assumptions, requires both a s t a t i s t i c a l knowledge about the target s i g n a l as well as knowledge about the ramifications of using a given method when those assumptions are not to some extent met. In the type of s t a t i s -t i c a l modeling employed i n t h i s study, d e t a i l s of which are provided i n subsequent sections, a s a t i s f a c t o r y understanding of these issues i s very important. The following section begins t h i s process of understanding by reviewing previous investigations into the s t a t i s t i c a l character of EEG. 2.2 Previous Stochastic Studies on EEG There have been r e l a t i v e l y few investigations into the s t a t i s t i c a l c h a r a c t e r i s t i c s of spontaneous EEG a c t i v i t y . The r e s u l t i n g conclusions from these few investigations have been l a r g e l y contradictory and i n d e c i s i v e . In general, most investigators were attempting to measure the degree of wide-sense s t a t i o n a r i t y and to estimate the amplitude p r o b a b i l i t y d i s t r i b u t i o n . McEwen and Anderson [12] d i d some e a r l y extensive work i n t h i s area. To t e s t for wide-sense s t a t i o n a r i t y they divided a given EEG segment into two halves and then c a r r i e d out a two-sample Kolmogorov-Smirnov (K-S) t e s t on both the sample amplitude and spe c t r a l d i s t r i b u t i o n functions. This t e s t required that both the amplitude and spe c t r a l d i s t r i b u t i o n s from each h a l f could not be s i g n i f i c a n t l y d i f f e r e n t for the whole EEG segment to be con-sidered as wide-sense stationary. They tested for the Gaussianity of a given EEG segment by using i t s amplitude d i s t r i b u t i o n i n a K-S goodness of f i t t e s t with unknown mean and variance using a 0.05 l e v e l of s i g n i f i c a n c e . They rejected t h e i r n u l l hypothesis that EEG from awake r e s t i n g subjects with eyes closed was Gaussian and wide-sense stationary approximately 15% of the time 9 over two second epochs and approximately 60% of the time over 8 second epochs. Persson [14] i n commenting on t h e i r r e s u l t s , pointed out that the s t a t i s t i c a l t e s t s that they used assume independent samples (observations) but the d i g i t i z a t i o n rates used res u l t e d i n samples that were hig h l y correlated. In f a c t , McEwen and Anderson noted that too high a sample rate would cause the e f f i c a c y of the s t a t i s t i c a l t e s t s to be adversely affected and they consequently recommended sampling at a rate as l i t t l e above the Nyquist rate as possible. Persson [14] went on to argue that the maximum t o l e r a b l e c o r r e l a t i o n c o e f f i c i e n t between adjacent samples i s about 0.5 and i n previous work he showed, based on an estimated autocorrelation function from r e a l EEG, that to meet t h i s requirement the sample rates should not be much greater than 10 Hz. The obvious r e s u l t i n g conundrum i s that i f only approximately two second epochs can be considered stationary and a sample rate i n the order of 10 Hz i s used then the r e s u l t i n g number of samples would be so small that a reason-able inference cannot be made about the amplitude d i s t r i b u t i o n . Weiss [15] approached t h i s problem by developing a correction f a c t o r , based on the second and fourth s p e c t r a l moments, for the Kolmogorov-Smirnov goodness of f i t t e s t which i s designed to compensate for c o r r e l a t i o n i n the data two sample points back i n time. He tested t h i s method on simulated EEG data which was generated by a second order autoregressive processes. Although he reports good r e s u l t s on t h i s simulated data, i t s usefulness i s s t i l l generally l i m i t e d by i t s a b i l i t y to compensate for the c o r r e l a t i o n over only two sample points. In addition, i t s effectiveness, i f applied to actual EEG, w i l l be further l i m i t e d by the accuracy of the estimated s p e c t r a l moments. In Section 2.4 a d i f f e r e n t approach i s discussed which would be 10 p o t e n t i a l l y more f l e x i b l e and u l t i m a t e l y provide more information about the EEG signal c h a r a c t e r i s t i c s . 2.3 Neural Basis for the Gaussian Nature of EEG R. E l u l was responsible for the i n i t i a l work devoted to the stochastic aspects of EEG based on neuronal a c t i v i t y . He f i r s t suggested [16,17,18] that each i n d i v i d u a l neuron generator was independent of the summed contribu-tions from a l l the neuron generators. Therefore, t h i s r e s u l t a n t sum, the EEG, could be thought of as the sum of s t a t i s t i c a l l y independent or nearly independent neuronal contributions and since the contribution from each neuron i s very small r e l a t i v e to the r e s u l t i n g EEG there must be a very large number of neurons contributing at any given time. Based on these arguments the a p p l i c a t i o n of the Central Limit Theorem (CLT) i s j u s t i f i e d ; that i s , the sum of neuronal a c t i v i t y w i l l tend toward Gaussianity. However, E l u l [19] l a t e r c a r r i e d out an experiment i n which he administered tetrodotoxin (TTX) into the b r a i n of cats. The amount of TTX that was given to the cats should have caused about a 10% drop i n neural a c t i v i t y . The r e s u l t i n g EEG a c t i v i t y was reduced way below the l e v e l that could be accounted for based on hi s concept of independent or nearly independent neural a c t i v i t y . A. Siegel [20] followed up on E l u l ' s work and he proposed the idea that a substantial proportion of the neurons belong to synchronized groups. These groups, however, would be n e c e s s a r i l y r e s t r i c t e d i n s i z e due to the existence of many competing inputs to a given neuron which would r e s u l t i n attenuation of the synchronizing e f f e c t as one moves along a chain of i n t e r -acting neurons. He further postulated that because of t h i s r e s t r i c t e d s i z e 11 there would s t i l l be a very large number of " i n t e r n a l l y synchronized but mutually unsynchronized groups of neurons". This explanation of neuronal a c t i v i t y was able to predict the dramatic reduction of EEG a c t i v i t y which occurred i n the TTX experiment. I t also allows for E l u l ' s basic idea of the summed a c t i v i t y being independent of the a c t i v i t y of the i n d i v i d u a l generators. E l u l [18] also suggested that d i f f e r e n t degrees of independence between neurons, as was alluded to e a r l i e r when the terms "independent" and "nearly independent" were used, would be the major influence on the degree of Gaussianity: as the dependence becomes greater the r e s u l t i n g d i s t r i b u t i o n becomes less Gaussian. Siegel [20] elaborates on t h i s concept and i n so doing suggests a mechanism that would produce t h i s r e s u l t . He u t i l i z e s Bernstein's [21] theory of applying the CLT to dependent v a r i a b l e s . Roughly speaking, i t states that as long as the dependence between va r i a b l e s decreases with separation then the CLT can be applied. Therefore, Siegel argues that as long as the dependency of two neuronal generators decreases quickly enough with increased separation, the a p p l i c a t i o n of the CLT can s t i l l be j u s t i f i e d . It also follows that, at periods of time when dependency i s l e s s , the e f f e c -t i v e number of independent contributors increases and the CLT i s more c l o s e l y approximated. E l u l [18] suggested the a p p l i c a t i o n of t h i s concept to various l e v e l s of mental a c t i v a t i o n : performing an a c t i v e mental task would require a greater degree of interneuronal coupling than would a mental i d l e state. Hence, the degree of Gaussianity would decrease during performance of mental tasks. He c a r r i e d out some empirical work with EEG, which showed some support for t h i s idea, but the s t a t i s t i c a l analysis suffered from the sample dependency problem that was described i n the previous section as well as 12 other methodological problems. As a t t r a c t i v e as the above concept at f i r s t seems, i t i s however, contradictory to a common assumption about EEG. As E l u l states: "low-vol-tage, f a s t a c t i v i t y implies 'desynchronized 1 (active) EEG, and high-voltage slow a c t i v i t y i s i n d i c a t i v e of 'synchronization'" ( i d l e EEG). Siegel [20] resolves t h i s apparent paradox with the following argument. In the i d l e state there are, as stated previously, groups of neurons which are i n t e r n a l l y synchronized but mutually independent, which r e s u l t s i n summation of a c t i v i t y that i s r e l a t i v e l y high-voltage and appears to be i n r e l a t i v e synchrony. During a mental task the r e l a t i o n s h i p between the "within group" neurons becomes more complicated than simple "in-step" synchrony. As Siegel [20] states: " E s s e n t i a l l y , t h i s i s because neurons must be r e l a t e d i n configura-tions which correspond i n complexity to that of the task i t s e l f . " So the e l e c t r i c a l neuronal a c t i v i t y , although more interdependent does not have the same appearance of synchronicity and hence there w i l l be a greater amount of neuronal a c t i v i t y canceling each other out r e s u l t i n g i n lower voltage EEG. In some recent work by Anninos, Zenone and E l u l [22], they studied neuronal a c t i v i t y and the r e s u l t i n g EEG based on a rigorous a r t i f i c i a l neural net model. Their p r i n c i p l e conclusion was that the main factor i n causing the summed neural a c t i v i t y to deviate from Gaussianity was i n f a c t the l e v e l of interneuronal connectivity: greater connectivity caused greater devia-tions from Gaussianity. This r e s u l t occurred independent of the p r o b a b i l i t y d i s t r i b u t i o n of the membrane p o t e n t i a l i n the i n d i v i d u a l elements. In addition, they discovered that for a given l e v e l of connectivity i n t h e i r model, when external input was applied, as would be the case when afferent signals were applied to the neural net, the r e s u l t i n g d i s t r i b u t i o n became 13 less Gaussian. This finding may add<some i n t e r e s t i n g i n s i g h t to the possible differences between the motor r e l a t e d p o t e n t i a l s i n normal and disabled persons. To date, i t appears that the idea of r e l a t i n g various l e v e l s of mental a c t i v a t i o n to various l e v e l s of Gaussianity has not been c a r e f u l l y confirmed or rejected by empirical measurements. This i s probably due to the f a c t , as noted i n Section 2.2, that a s a t i s f a c t o r y t o o l to measure the EEG s t a t i s t i c s does not seem to be a v a i l a b l e . 2.4 Applying AR Modeling to EEG A very general l i n e a r model for the modeling of stochastic d i s c r e t e -time processes i s the autoregressive moving average (ARMA) model. It i s given by P I X i = E a k X i - k + E b j e i - j 2 ' 1 k=l j=0 where x_^  i s the di s c r e t e s i g n a l sequence of length n, i = l,2...n, e_^  i s the r e s i d u a l error sequence and a, , k=l,2...p and b., j=0,l,2....1 are weighting parameters on past values of the si g n a l and residuals r e s p e c t i v e l y . The autoregressive (AR) model i s the " a l l pole" version of the ARMA model and i t has the following form x. = a.x. , + a„x. „ + ... + a x. + e. 2.2 i 1 l - l 2 i-2 p l - p i where again x^ i s the si g n a l sequence, a^ are weighting parameters, e^ i s the white re s i d u a l error sequence and p i s the order of the AR model. The s i g n a l x^ at a g i v e n time i i s assumed to be a l i n e a r l y weighted sum of p 14 past v a l u e s of x. plus a random (white) error term e.. Often t h i s l a s t term 1 1 i s r e f e r r e d to as the p r e d i c t i v e error since i t i s the differe n c e between the measured value and the predicted value. The AR model has some s i g n i f i c a n t p r a c t i c a l advantages over the more general ARMA model. F i r s t l y , a closed form s o l u t i o n to the minimization problem for the estimation of the ARMA model parameters does not e x i s t and hence i t e r a t i v e numerical optimization approaches must be u t i l i z e d . Whereas, i n the case of the AR model, the closed form s o l u t i o n of the minimization problem does e x i s t and computationally e f f i c i e n t methods have been developed to estimate the AR model parameters. Secondly, as noted by Kay and Marple [23] , the Wold decomposition theorem demonstrates that any stationary ARMA process (in f a c t , any MA processes as well) of f i n i t e variance can be repre-sented by a unique AR model which may be of i n f i n i t e order. The imp l i c a t i o n i s , even i f i t i s argued that for a given s i g n a l an ARMA model i s the most appropriate model, a reasonable approximation can s t i l l be achieved by u t i l i z i n g an AR model with an appropriately chosen model order. S i m i l a r l y , i n a study by Beamish and P r i e s t l y [24] they note that the time s e r i e s does not have to exactly conform to a f i n i t e AR model but rather assumes i t can be modeled by an i n f i n i t e AR model. Then by choosing an appropriate order which w i l l provide i n some sense the optimal f i t with a f i n i t e model, the si g n a l can be well represented. Selection of t h i s appropriate model order i s a very important issue and i t i s dealt with i n d e t a i l i n Section 4.4. Previous work has indic a t e d that AR modeling would prove to be a useful t o o l i n the i n v e s t i g a t i o n of the EEG s i g n a l . Jansen et a l . [25] note the s t a t i s t i c a l d e f i n i t i o n of regression i s : "a functional r e l a t i o n s h i p between two or more correlated v a r i a b l e s used to pre d i c t values of one v a r i -15 able when given values of the others." I f these variables are time r e l a t e d , as would be the case with EEG, then an autoregressive model can be applied [26]. In general, i t i s not clear how to d e f i n i t i v e l y assess whether an AR model at a given model order i s adequately representing a segment of EEG. For purposes of t h i s study, the appropriateness of AR modeling w i l l be discussed i n terms of i t s r e l a t i v e performance when i t i s applied to the spec t r a l estimation of EEG (see Section 4.4). The AR model as applied to EEG can be u t i l i z e d i n several ways. The estimated parameters could p o s s i b l y be used as features i n a s i g n a l detec-t i o n / d i s c r i m i n a t i o n problem. As indicated above, the parameters can be used i n spectrum estimation. Many researchers [23,25,27] have demonstrated that there are some d i s t i n c t advantages of t h i s approach over the conventional FFT methods of spectrum estimation. Another benefit of applying AR modeling to EEG i s the fact that r e s i d u a l s , i d e a l l y , have the c o r r e l a t i o n of the process removed (whitened) and since the AR process i s l i n e a r , the s t a t i s t i c a l char-a c t e r i s t i c s of the o r i g i n a l process are s t i l l contained i n the r e s i d u a l s . Although Andrews [28] demonstrates that residuals from a non-Gaussian process w i l l tend to mask the evidence of non-Gaussianity, Chambers and Heathcote [29] have developed a method of characterizing the Gaussianity of a process based on a scale factor which i s determined by the c h a r a c t e r i s t i c function of the r e s i d u a l error d i s t r i b u t i o n . The main benefit of t h i s approach i s that i t overcomes the problems of corr e l a t e d data samples, as was discussed i n Section 2.2. The greatest benefit of applying the AR model to EEG i n terms of t h i s t h e s i s work l i e s i n the f a c t that i t has a very convenient state-space representation, which allows for the s t r a i g h t forward use of state-space 16 techniques such as Kalman type f i l t e r i n g [30]. These techniques play a major ro l e i n the EEG sin g l e t r i a l processing scheme that i s discussed i n d e t a i l i n Chapter A. 17 CHAPTER 3 THEORY OF AR MODELING 3.1 Conventional AR Parameter Estimates EEG s i g n a l c h a r a c t e r i s t i c s are changing over time and hence, a sing l e time-invariant model can not be applied. This r e s u l t s i n the need to estimate model parameters from the EEG si g n a l i n a manner that w i l l attempt to account for time varying c h a r a c t e r i s t i c s . The estimation of the AR model parameters can generally be c a r r i e d out eit h e r by block mode estimation or by recursive estimation. Recursive methods sequentially update the parameters data point by data point. They have the p o t e n t i a l advantage of being set up such that the estimation of parameters adapts to time varying c h a r a c t e r i s t i c s of a s i g n a l [31]. This i s e s s e n t i a l l y accomplished by assigning greater weight to newer information than to older information. I t has been demon-strated for EEG s i g n a l s , which are slowly time-varying over long epochs, that adaptive recursive estimation schemes can be e f f e c t i v e [32,33]. However, Jansen [25,32] provides evidence which indicates that the adaptation process i s not rapid enough for short segment analysis of EEG signals that are expected to be time-varying r e l a t i v e l y quickly. Short segments i n the order of 1 to 2 seconds are t y p i c a l for the work c a r r i e d out i n t h i s t h e s i s . Hence, a block mode approach, where the AR model parameters are estimated on the basis of a short data segment, i s employed throughout t h i s work. Various block mode methods have been used to estimate a set of AR model parameters from a sample s i g n a l segment. The following discussion describes the most common conventional methods for AR parameter estimation. 18 In the past, the standard method to estimate AR parameters was often based on the Yule-Walker equations. These equations provide a r e l a t i o n s h i p between the autocorrelation function and the AR model parameters (see Ulrych and Bishop [34]. The d e r i v a t i o n of these equations i s reviewed by Kay and Marple [23] and the f i n a l r e s u l t expressed i n matrix form i s 3.1 R (0) R (-1) X X X X . . R U-p) X X R (1) X X R (1) R (0) X X X X . . R (2-p) X X a2 R (2) X X R (p-1) R (p-2) . X X r X X r . . R (0) X X ap R* (p) X X where R ^ t k ) i s the a u t o c o r r e l a t i o n f u n c t i o n for lag k, p i s the AR model o r d e r , and a^ k=l,2...p are the AR model parameters. Therefore, by obtaining estimates of the autocorrelation function, estimates of the AR parameters can be obtained by solving the system of Equations 3.1. Kay and Marple recommend, to achieve low mean-squared error, estimation of the autocorrela-t i o n function at s p e c i f i c lags with the following expression . n-m-1 R (m) = - I x.^ x. xx . l+m I n i=0 3.2 Note also that for a stationary process the conjugate symmetry property of the a u t o c o r r e l a t i o n function R (m)=R (-m) can also be u t i l i z e d i n solving X X X X 6 3.1. In addition, the Yule-Walker equations y i e l d an expression that allows the variance of the residuals to be c a l c u l a t e d [23] R (0) = I a. R (k) + 6* xx . * k xx e k=l 3.3 where 6* i s the v a r i a n c e of the r e s i d u a l sequence. E q u a t i o n 3.1 can be augmented with 3.3 to y i e l d the following a l t e r n a t i v e form of the Yule-Walker equations [23] 19 R (0) xx R (1) xx R (p) xx R (-1) xx R (0) xx R (p-1) xx R x x ( - ( p - l ) ) R (Oj X X a LPJ 6* e 0 3.A Least square (LSQ) estimation of the AR parameters i s another very common method. In f a c t , as w i l l be shown l a t e r , a l l the conventional methods discussed i n t h i s s ection can be shown to be based on le a s t squares minimiza-t i o n c r i t e r i a . The method that i s most commonly r e f e r r e d to as the "LSQ" method can be derived i n the following manner [23]. From Equation 2.2 the p r e d i c t i o n error can be written as e i = x i " a k x i - k k=l 3.5 The sum of squared p r e d i c t i o n errors i s then SSE n i = l 1= X [ x i - X a k x i - k ] z 3.6 i = l * k=l To f i n d the AR parameters that minimize 3.6 the p a r t i a l d e r i v a t i v e s with respect to each a^ are taken and set equal to zero. That i s 3(SSE) 3 a = 0 q=l,2. 3.7 The r e s u l t of applying 3.7 to 3.6 i s P n y a. y x. ,x. , k . , l-k l-c n I x,- x k=l i = l i = l l l - q q= 1,2... p 3.8 Then by s u b s t i t u t i n g 3.8 into 3.6, the minimum SSE can be shown to be [23] 3.9 Now by expanding 3.8 in t o matrix form for model order p r e s u l t s i n n p n SSE . = Y. x ? + E a i E x - x - i mm . . I i . k . , I l - k i = l k=l i = l 20 £ x i - i x i - i £ x i-2 x i - r Ixi_1xi_2 ^ xi-2 Xi-2 Ex. X . . ^ 1-p 1-J Z x i - p z x i - p j a LPJ £ X i X i - l EXiXi-2 Zx.x. _ 1 l-p_ 3.10 _ i - i i -p By c a l c u l a t i n g the summations i n 3.10, from a data sequence of n points, t h i s system of equations can be solved to determine the LSQ estimate of the AR parameters. Note that i f the autocorrelation estimate given i n 3.2 i s used, except for the s c a l i n g factor, 1/n, which does not e f f e c t the s o l u t i o n for the AR parameters, the above equations reduce to the Yule-Walker method of parameter estimation. Therefore, the Yule-Walker estimates are equivalent to the LSQ estimates for s u f f i c i e n t l y large n. An i m p l i c i t assumption, which i s p a r t i c u l a r l y evident i n the develop-ment of the Yule-Walker s o l u t i o n , i s that the autocorrelation function i s assumed to be zero outside the data segment of i n t e r e s t . In p r a c t i c e , when r e l a t i v e l y short segment lengths are used, t h i s truncation of the auto-c o r r e l a t i o n function can r e s u l t i n r e l a t i v e l y poor parameter estimates [25,23], Burg [35] addressed t h i s problem by using extrapolation of the autocorrelation function based on concepts of maximum entropy and he formu-lat e d a method of AR parameter estimation, known as the Burg method or the maximum entropy method (MEM). Many authors have noted the superior perform-ance of t h i s method over the Yule-Walker type methods when applied to r e l a t i v e l y short data segments (for example see [23,25,36,37]). Therefore, i n t h i s thesis work, because t y p i c a l l y r e l a t i v e l y short epochs of EEG are u t i l i z e d , MEM was selected as the conventional AR parameter estimation method. 21 The fundamental idea behind Burg's method i s to provide a non-zero e x t r a p o l a t i o n of t h e ^ a u t o c o r r e l a t i o n function beyond the known lags (up to and including the p lag: see Equation 3.1) as opposed to the implied zero extrapolation as i n the Yule-Walker equations. Burg argued that the extrapo-l a t i o n should impose the fewest possible constraints on the extrapolated autocorrelation function without compromising any information about the known lags. To achieve t h i s , he required that the hypothetical time s e r i e s , which would be represented by the extrapolated autocorrelation function, should have maximum entropy. This requirement maximizes the randomness of that time s e r i e s , given the constraints of the estimate of the function, and hence produces a minimum bias s o l u t i o n . From a spectral point view, the maximum entropy estimate i s based on choosing a spectral estimate such that the entropy (E) per sample l / 2 f E = J" S In F (f) df 3.11 - l / 2 f X s where F (f) i s the spe c t r a l estimate of the data segment and f i s the sample X s frequency, i s maximized subject to l/2f J F (f) exp-(j2TTk—) df = R (k) k=l , 2 , . . .p 3.12 - l / 2 f X f X X s s It can be shown [38] that the spe c t r a l estimate which maximizes entropy subject to the constraint that i t s f i r s t p Fourier c o e f f i c i e n t s correspond exactly to the sample autocorrelation function evaluated at the f i r s t p lags i s the estimate of the spe c t r a l density function of an AR model of order p. In addition, i t i s shown that the estimates of the parameters and the power 22 i n the r e s i d u a l sequence, P , can be written as R (0) . . . R (p) X X xx r R (p) X X . R (0) X X .1 p a 1 p 1 — 0 a 0 Lp_ _ _ 3.13 Note that these equations are of the same form as the augmented Yule-Walker Equations 3.4 where the variance of the residuals for the zero mean re s i d u a l sequence i s equal to P . The algorithm that was u t i l i z e d i n t h i s thesis work to c a l c u l a t e MEM parameter estimates i s based on a procedure outlined by Anderson [36]. The system of Equations 3.13 i s solved i n a sequential manner. Beginning with p=0, P , i s estimated by 1 n P = - X x . 2 < 3.14 o . , 1 n i = l Then the model for p=l i s determined as that which minimizes the power i n the forward p r e d i c t i o n error sequence averaged together with the power i n the backward p r e d i c t i o n error sequence. This average power i s given by 1 1 n I [(x.-a 1(l)x.. 1) 2 + ( x i _ 1 - a 1 ( l ) x i ) 2 ] 3.15 2 n-1 i=l+l For the general case of progressing from order p-1 to order p, i t i s shown that 1 1 2 n-p i=p+l k=l + ( x i - p " J ^ k ^ i - p + k ^ 3.16 23 and i t can be minimized w i t h r e s p e c t to the si n g l e parameter a^(p). Note that the arguments on the parameters i n d i c a t e the current maximum model order at a given stage i n t h i s sequential procedure. The dependence of the other model parameters at the pth stage are given by the Levinson recursion a k(p) = a k ( p - l ) - a p ( p ) a p _ k ( p - l ) 3.17 3TT A p p l y i n g E q u a t i o n 3.17 to 3.16, t a k i n g ^—j = 0 and solving f or a (p) ap P P r e s u l t s i n 2 _ I b ( i - l ) e (i) a (p) = 1 = P + 1 3.18 p r n I ( |b . ( i - D I ' + le , ( i ) l 2 ) • i i P ~ l P ~ l i=p+l P P where b ( i ) = Y a (i)x(i-p+k) and e (i)= Y a ( i ) x ( i - k ) are the backward and P i = 0 P P i = 0 P forward p r e d i c t i o n e r r o r s . In the above notation i t i s assumed that a (0)=1.0. Hence, the value of ap(p) i s calculated v i a Equation 3.18, with the remaining parameters being calculated v i a Equation 3.17 and then the backward and forward p r e d i c t i o n error sequences are updated. The order i s then increased to p+1 and the procedure i s repeated u n t i l the desired model order i s reached. van den Bos [38] noted that MEM i s equivalent to a l e a s t squares f i t t i n g of an AR model. More s p e c i f i c a l l y , Kay and Marple [23] note that MEM can be viewed as a constrained l e a s t squares minimization problem because, as i s i n d i c a t e d by Equation 3.15, the sum of both the forward and backward error energies (squared error terms), as given below 24 n n E = I le (i) 12 + I lb Ci) 12 3.19 i=p+l P i=p+l p i s the minimization c r i t e r i a for MEM, 3.2 LSQ Parameter Estimation on Simulated Gaussian Data LSQ methods, as i s shown i n Section 3.3, are optimal when applied to a Gaussian random process i n terms of maximum l i k e l i h o o d estimation which provides the best l i n e a r unbiased estimates. Simulation studies using LSQ parameter estimation were c a r r i e d out on computer generated Gaussian AR data. The simulated data were generated using a Gaussian sequence, e^ with variance 1.0 and mean 0.0, d r i v i n g an 8th order AR model with the following parameters: a, = 0.838 a. = -0.471 a. = 0.638 a, = -0.429 1 2 3 4 a c = 0.518 a, = -0.304 a.n = 0.182 a_ = -0.243 This set of parameters was c a l c u l a t e d from EEG data and the set was selected as being t y p i c a l for an 8th order AR model of actual EEG. The Gaussian sequence e^ was generated from a uniform white sequence v i a the procedure given by Box and Muller [39]. The uniform pseudo-random number generator was based on a "generalized feedback s h i f t r e g i s t e r algorithm" which i s given by T.G. Lewis and W.H. Payne [40] . They demonstrate that t h i s generator has excellent random properties and a sequence period of 2 r a i s e d to the j t h power, where j i s the integer word length which i n t h i s case was 31. F i f t y sets of estimations using d i f f e r e n t random data segments of lengths n=1000, n=500, n=100, n=64, and n=32 were c a r r i e d out. The mean squared error (MSE) between the actual and estimated parameter values were 25 calculated for each set. The r e s u l t s are summarized i n Figure 3.1. These r e s u l t s show that the MSE began to increase s i g n i f i c a n t l y when the segment length was decreased to n=100 and they became very large when the segment length was reduced to n=32. This points to, i n terms of parameter estimation e f f i c a c y , a p r a c t i c a l lower bound on segment s i z e . I d e a l l y , segment sizes i n the order of n=500 would be preferable, but i n the case of actual EEG, there i s a need (see Section A.5) to make the segment s i z e as small as possible. The EEG data was sampled at 64 Hz (see Section 4.1), and therefore, a one second segment contains only 64 data points. Since the MSE and SE became so large at n=32 and since i t was s i g n i f i c a n t l y less at n=64, i t was decided that one second (n=64) should be the absolute lower bound on segment length for EEG data. Ultimately, a segment s i z e of 1.5 seconds (n=96), moderately above t h i s lower bound, was u t i l i z e d < i n the s i n g l e t r i a l analysis method (see Section 4.4). Further discussion on segment length r e l a t i n g s p e c i f i c a l l y to the EEG data used i n t h i s t hesis work i s provided i n Section 4.5. 3.3 Deviations from Gaussianity Several researchers, such as McGillem et a l . [13], Jansen et a l . [25] and Smith and Lager [41], employed AR parameter estimation i n various a p p l i -cations involving EEG. In t h i s work, the required parameter estimation was based on estimation methods involving least squares, which are generally optimal i n Gaussian processes. However, i t has been shown [42] that the performance of these methods can, under c e r t a i n circumstances, s i g n i f i c a n t l y break down under even s l i g h t deviations from Gaussianity. Although i t was discussed e a r l i e r i n Section 2.2 and 2.3 that there was considerable evidence 0.023 0.021 0.019 0.017 0.015 0.013 0.012 0.010 0 .008 o Z3 Z> 0 .006 c 0 .004 o OJ 0 .002 ^ 0 .000 LEAST SQUARES AR PARAMETER ESTIMATION Gaussian N=50 3 4 5 Parameter Number 7 8 . . LSQ NPTS=100 LSQ NPTS=500 Figure 3.1a LSQ NPTS=1000 LEAST SQUARES AR PARAMETER ESTIMATION Gaussian N=50 0.0904 0.0829 0.0753 — (- 0.0678 o c 0.0603 < _ L x J 0.0527 - T J 0.0452 C U c_ 0.0377 a a 0.0301 co 0.0226 — c 0.0151 o O J 0.0075 0.0000 - / / / - / / / / \ 3 4 5 Parameter Number 7 8 LSQ NPTS=32 LSQ NPTS=64 Figure 3.1b .SQ NPTS=100 tSJ 28 to expect EEG to be generally Gaussian, there i s c e r t a i n l y no guarantee that sample segments of EEG w i l l always be s t r i c t l y Gaussian. This prompts questions about both the Gaussian nature of sample segments of EEG, i n p a r t i c u l a r r e l a t i v e l y short segments, and the r e s u l t i n g performance of AR model parameter estimation techniques under conditions of varying degrees of deviation from Gaussianity. A considerable amount of work has been c a r r i e d out on robust estima-t i o n of l o c a t i o n parameters and l i n e a r regression model parameters i n the case of independent and i d e n t i c a l l y d i s t r i b u t e d ( i . i . d . ) observations. The basic goal of a robust procedure, for the purposes of t h i s t h e s i s , i s to provide good estimates when the data has a small number of o u t l i e r s ( i n the order of 5 to 20 percent) causing the assumed, Gaussian i n the case of EEG, d i s t r i b u t i o n function to be contaminated. In addition, the robust procedure should provide estimation r e s u l t s which are not s i g n i f i c a n t l y d i f f e r e n t from the conventional LSQ methods when the data i s not contaminated with o u t l i e r s . Location Case There are a number of robust methods but the most s a t i s f y i n g to date appear to be those given by modifications to maximum l i k e l i h o o d . Hogg [43] provides a good background t u t o r i a l on robust methods and the following b r i e f d e s c r i p t i o n of modified maximum l i k e l i h o o d methods i s based on that t u t o r i a l . 11 x ^  t ^ 2 * **••»* x^ aire a random sample from a p r o b a b i l i t y density function f(x-O) where 0 i s a l o c a t i o n parameter, then the logarithm of the l i k e l i h o o d function i s given by n n JnL(0) = I inf(x.-e) = - I p(x.-0) 3.20 i = l 1 i = l 1 The maximum l i k e l i h o o d method maximizes In L(0) or i n terms of the p function 29 minimizes n I p(x.-G) = K(0) . 3.21 i= l 1 Assume that the minimization can be achieved by d i f f e r e n t i a t i n g and solving K'(9) = 0. In other words, f i n d the value of G that s a t i s f i e s n I V> (x.-9) = 0 3.22 i = l 1 where \|>(t) = P 1 (t) = - f 1 ( t ) / f (t) 3.23 The value of 0 that minimizes K(0) i s termed the maximum l i k e l i h o o d e s t i m a t e of 0 and i s denoted as 0. Robust M-estimates are generated v i a Equation 3.22 except that d i f f e r e n t p s i functions are used than that described i n Equation 3.23. Each d i f f e r e n t p s i function describes a s p e c i f i c type of M-estimate. The basis of robust M-estimates i s to f i n d p s i (\Ji) functions that w i l l protect against o u t l i e r data points that cause undue influence on the estimation result.* An example of such a p s i function i s that due to Huber [44]. It i s designed to deal with data that i s d i s t r i b u t e d normally i n the middle with double exponential t a i l s ("heavy-tailed" d i s t r i b u t i o n ) . The p s i function i s given by -c t < -c ij)(t) = { t | t | <. c 3.24 c t > c where c i s a tuning constant. A scale invariant version of the M-estimator i s given as n x.-0 I = 0 3.25 i = l s where s i s a robust estimate of the process scale. In t h i s case, c = 1.5 30 i s often selected to allow data from a t r u l y Gaussian d i s t r i b u t i o n to be uninfluenced by t h i s p s i function while s t i l l providing the desired protection from outlying points. The so l u t i o n to Eqn. 3.25 i s t y p i c a l l y found by a block mode i t e r a t i o n scheme, such as the iterated-weighted least squares (IWLS) procedure [45]. Other commonly used p s i functions are Hampel's three part redescending and Tukey's Bi-weight. It should be noted that, the c a l c u l a t i o n of the p s i function for a Gaussian d i s t r i b u t i o n demonstrates that the LSQ estimate i s the maximum l i k e l i h o o d estimator. The Gaussian d i s t r i b u t i o n i s given by [46] r - ( x - 0 ) * n  e x p [ 2o» ] f(x-O) = — — 3.26 V2TTOJ and hence p(x-G) = \ fin(2Tro*) + ( * ~ f } 2 3.27 and \b(x-0) = (x-O) 3.28 Applying Eqn. 3.22 i t follows that the maximum l i k e l i h o o d estimator i s the value of 0 that s a t i s f i e s 1 n ±7 I (x.-6) = 0 s 3.29 a2 . , I i = l y i e l d i n g the well known r e s u l t 0 = x 3.30 which i s exactly the LSQ estimate of 0 [43]. Regression Case The above methods can be extended to the block mode l i n e a r regression case (see Hogg [43]). Given the l i n e a r model 31 Y = Xa + e 3.31 where: X i s n x p data matrix a i s a parameter vector of order p e i s a r e s i d u a l vector of order n p i s the order of the model n i s the number of data points It then follows that the expression to be minimized i s now T n y. - x. a I P( 1 a 1 ) 3.32 i = l e where: y_^  i s the i t h data point x. i s the i t h row of the matrix X - i Considering the p f i r s t p a r t i a l d e r i v a t i v e s , the following set of p equations must be solved T n y -x a _ I x. * ( 1 1 ) ~ 0 3.33 i = l e where s i s a r o b u s t estimate of the standard deviation of the r e s i d u a l s , e Again, t h i s set of equations can be solved using a IWLS procedure and Hogg [43] recommends the following robust estimate median |e. - median(e.)I s = — 3 34 e 0.6745 * Autoregression Case The a p p l i c a t i o n of robust methods to time-series data has lagged behind the a p p l i c a t i o n to the i . i . d . case, probably, as suggested by Martin [47], because of the considerable d i v e r s i t y i n q u a l i t a t i v e features of time-serie s data sets as well as the possible dependency that may a r i s e i n the residuals due to data o u t l i e r s that occur i n patches (correlated). Martin 32 [47] applies robust estimation to time-series data and shows that many of the concepts from the i . i . d . observations case can be applied d i r e c t l y . The greatest dif f e r e n c e seems to occur i n the d e f i n i t i o n s of the robust q u a l i t i e s of the given methods due mostly to the added d i f f i c u l t i e s mentioned above. Martin [47] also points out that for data to q u a l i f y as a time-series o u t l i e r i t only has to be " d i f f e r e n t " on the innovations (residual) scale not the process scale. Since the innovations scale i s t y p i c a l l y 10 to 10,000 times smaller than the process scale the o u t l i e r s w i l l often be impossible to v i s u a l l y detect i n a p l o t of the raw data. Martin [47] i d e n t i f i e s three types of o u t l i e r s that may occur i n time-seri e s data: 1) independent i s o l a t e d gross-error o u t l i e r s which may be caused by various recording (measurement) errors 2) patchy type o u t l i e r s whose behavior seems to be uncorrelated with the behavior of the rest of the data - t h i s may be caused by b r i e f malfunc-tions i n the data c o l l e c t i o n system, inherent behavior of the process or maybe other unaccountable e f f e c t s 3) patchy o u t l i e r s whose behavior does appear to be r e l a t e d to the r e s t of the data with the possible exception of an i n i t i a l jump - t h i s type of o u t l i e r may be caused by unusual events within the process He also suggests that two types of o u t l i e r models can reasonably simulate the above types of o u t l i e r a c t i v i t y . The additive o u t l i e r model would apply for types 1 and 2 while the innovations o u t l i e r model would apply for type 3. The Innovations O u t l i e r (10) model i s described i n [47] as P x. = £ a, x. . + e. 3.35 l i . k l-k l k=l where the i n n o v a t i o n sequence e^ i s i . i . d . with a symmetric d i s t r i b u t i o n G and the observations are given by 33 y. = x. . 3.36 J 1 l Innovation o u t l i e r s occur when G i s heavy-tailed. Martin and Thompson [45] have found that t h i s type of deviation from Gaussianity does not cause, except perhaps i n extreme cases, serious problems for the conventional estimation of parameters from a time-series. In the case of the Additive O u t l i e r (AO) model, the observations are given by y. = x. + v. 3.37 Ji l l where x^ i s defined as i n 3.35 with G Gaussian, v^ i s independent of x^, and v^ has a symmetric d i s t r i b u t i o n . A s u i t a b l e d i s t r i b u t i o n f o r v^for the i . i . d . case i s the contaminated-normal with degenerate c e n t r a l component which has the following form [45] CND(r,o 2) = (l-r)N(0,0) + rN(0,o 2) 3.38 where y i s the proportion of contamination and the notation N(u,o 2) represents a normal (Gaussian) d i s t r i b u t i o n with mean u and variance a 2 . With t h i s d i s t r i b u t i o n , the p r o b a b i l i t y that v^ = 0 i s the p r o b a b i l i t y that y^ = x^ which equals 1-/ and i n t y p i c a l a p plications 0.01 < y < 0.25. In contrast to the 10 case, i t has been shown [45,47,48] that conventional time-serie s parameter estimation i s h i g h l y non-robust under t h i s additive type of contamination. Although, v^ has been r e s t r i c t e d above to the i . i . d . case i t has been found [45] that schemes which deal well with t h i s type of o u t l i e r also deal best with the patchy type of o u t l i e r s . Furthermore, Martin and Thompson [45] point out that, i n p r a c t i c e , the d e t a i l s of the o u t l i e r d i s t r i -b u t i o n f o r v^ are l a r g e l y i r r e l e v a n t because i t would be a poor robust estimator that depended s i g n i f i c a n t l y on a given d i s t r i b u t i o n f or v^. A generalized robust M-estimate (GM-estimate) [47] , a member of the 34 class of bounded-influence autoregression (BIFAR) estimates, was shown to have robust q u a l i t i e s for the AO case and i s given by [42] T I 5 W(x > « s" i" P") = 0 3.39 i=p+l e Note t h a t except for the weighting factors W(x^), these equations are of the same form as the equations given for the M-estimate. As described by Martin [47] , the r o l e of the a d d i t i o n a l weighting factors i s to down-weight the T summands of E q u a t i o n 3.39 f o r which x^a i s a poor predictor because one or more of the values i n x^ are too large. He shows that an appropriate calcu-l a t i o n of these weights can be achieved by l e t t i n g WCx^) = wtcL) 3.40 where d. i s defined as l T -1 3.41 x. C x. d»(x.) = — -- i P and C i s the pxp covariance matrix for the p ^ order AR model of the process and w(.) i s a non-negative decreasing weighting function, t y p i c a l l y of the form w(t) = c \|)(t/c)/t 3.42 where c i s a tuning constant. Hence, there are at least two tuning constants required: one for the p s i function and one for the w function. The value d 2 T -1 i n 3.41 i s p r o p o r t i o n a l t o x^C x^ ; t h i s e x p r e s s i o n , known as the Mahalanobis distance, provides, as noted by Martin [47], a natural metric by which the r e l a t i v e "largeness" of x^ can be determined. The GM-estimates can be solved using an IWLS procedure as given below [45] 35 where I x. W(x. ) B^  (y. - xT a J + 1 ) =0 3.43 i=p+l r r r j = 1,2, ... NIT T ' j y. - x. a J *( 1 - 1 ~ P " ) e * = s J . 3.44 l e T " j y. - x. a J  J l — i - p — The process to obtain GM estimates requires s t a r t i n g with model order p=l and sequentially increasing p u n t i l the desired order i s reached. At each model order i n t h i s sequence, a set of p equations r e s u l t i n g from Equation 3.43 are solved for each of the NIT i t e r a t i o n s . MEM estimates for a and the c o r r e s p o n d i n g r o b u s t e s t i m a t e of s g are used as s t a r t i n g values for the i t e r a t i o n s of Equation 3.43. The reason that the model parameters must be estimated sequentially i s due to the manner i n which a robust estimate of C" 1 i s c alculated. It has been found [47] that a successful approach to obtain a robust estimate of C" 1 i s based on the f a c t o r i z a t i o n C" 1 = A TA 3.45 where A i s upper t r i a n g u l a r and i s given by -a(p-k) X K % > k s e(p-k) s g(p-k) 0 I < k where k = 1,2 p-1 and a(p-k) , s e(p-k) are the parameter estimates and r e s i d u a l standard deviation for model order p-k and the required s t a r t i n g v a l u e of s e 2 ( 0 ) i s set equal to the variance of the o r i g i n a l data sequence 36 y\ . Hence, the pxp matrix C" 1 i s represented i n terms of the AR parameters and the corresponding r e s i d u a l standard deviations derived for model orders up to p-1. Therefore, by f i t t i n g AR models i n succession, the p r i o r GM parameter estimates at model order p-1 w i l l enable the construction of A P, the A matrix at the p ^ i t e r a t i o n , which provides the robust estimate of C" 1 v i a Equation 3.45 to be used for the current GM parameter estimate at order P-3.A GM Estimation on Simulated Contaminated Gaussian Data Simulation studies were c a r r i e d out to determine the r e l a t i v e performance of robust GM estimation methods as compared to the conventional LSQ (MEM) estimation method when applied to eighth order AR Gaussian data with 0%, 10% and 20% l e v e l s of AO contamination. Additive o u t l i e r s are studied because conventional parameter estimation i s not robust under t h i s type of contamination (see Section 3.3) and also because, the sing l e t r i a l a nalysis method (see Section 4.2) s i s l a r g e l y based on additive o u t l i e r concepts. The contamination for these studies was produced based on the AO model given i n Equation 3.37 with the d i s t r i b u t i o n G being Gaussian with v a r i a n c e 1.0 and mean 0.0. The d i s t r i b u t i o n for v. was of the form given i n l & Equation 3.38 with o2=2.0 and T=0, 0.1, and 0.2. The simulations were c a r r i e d out on data segments with a length equal to 100 points because t h i s r e f l e c t s the segment length used on the actual EEG signals as described i n Chapter 4. Figures 3.2, 3.3, and 3.4, show the MSE performance of the LSQ, GM, GM1, and GM2 estimation methods on f i f t y (N=50) random simulated 8th order AR 37 Gaussian processes with 0%, 10%, and 20% l e v e l s of AO contamination respec-t i v e l y . The GM1 and GM2 methods are extensions to the GM estimate which w i l l be described i n Section A.3. In these studies, as suggested by Martin and Thompson [A5] , the Huber p s i function, as described i n Equation 3.2A, was used for the f i r s t two i t e r a t i o n s (j = l,2) of Equation 3.A3 and then Tukey's bi-weight M t ) m r t ( l - (t/c)»)» | t | < c V{t> 1 0 It| * c . J , 4 / was used for the l a s t i t e r a t i o n (j=3=NIT). Through t r i a l and error, the various tuning constants were selected to provide the best performance i n term of MSE and are summarized below c Huber p s i 1.0 Tukey p s i 3.0 w based on Huber p s i 1.3 These r e s u l t s show that the robust methods perform, i n terms of MSE, almost as well as the LSQ method on the uncontaminated Gaussian data. With a 10% l e v e l of AO contamination the LSQ performance f a l l s o f f dramatically. The robust methods perform much better with the GM2 method performing almost as well as i n the uncontaminated case. At the 20% l e v e l of contamination the performance across a l l the methods has dropped s i g n i f i c a n t l y . However, the GM2 method i s s t i l l c l e a r l y the best performer with a MSE generally less than the LSQ method i n the 10% contamination case. The r e s u l t s support the expectations of robust estimation i n that the robust methods perform nearly as well as the LSQ method on uncontaminated data but perform s i g n i f i c a n t l y better than the LSQ method on contaminated data. 0.1200 0.1100 0.1000 0.0900 0.0800 0.0700 0.0600 0.0500 0.0400 0.0300 0.0200 0.0100 0.0000 L_ RR P A R A M E T E R E S T I M A T I O N G A U S S I A N N=50 3 4 5 Parameter Number LSQ GM GM1 Figure 3.2 AR P A R A M E T E R E S T I M A T I O N A D D I T I V E O U T L I E R 10'/. VAR.=2 ,0 N=50 0.122 0.112 — 0.102 0.092 o c _ 0.082 c _ UJ 0.071 -a 0.061 a) <_ 0.051 o r> 0.041 — cr co 0.031 — tr 0.020 o 0.010 . 0.000 3 4 5 Parameter Number 8 ._. LSQ GM . GM1 GM2 F i g u r e 3 . 3 RR P A R A M E T E R E S T I M R T I O N A D D I T I V E O U T L I E R 20% VAR.=2 .0 N=50 . . . LSQ GM GM1 GM2 Figure 3.4 41 CHAPTER 4 OUTLIER PROCESSING OF SINGLE TRIAL EEG 4.1 Experimental Design and EEG Data A c q u i s i t i o n The objective of the experimental data a c q u i s i t i o n was to obtain EEG signals from subjects during a c o n t r o l l e d voluntary s k i l l e d motor a c t i v i t y (active task) and during a c o n t r o l l e d state i n which the subjects were a l e r t but not involved i n any motor a c t i v i t y ( i d l e task). For the active task, subjects placed t h e i r r i g h t hand i n an apparatus which oriented t h e i r hand i n a standard p o s i t i o n . The task required them to aim for a "target" p o s i t i o n by performing a slow smooth (ramped) extension with t h e i r r i g h t thumb. During the movement the t i p of t h e i r thumb pressed against a lever that provided a small opposing force to the thumb movement. In the s t a r t i n g p o s i t i o n the lever rested on a support so that there was no i n i t i a l load and the thumb was i n a relaxed state. A potentiometer attached to the lever provided p o s i t i o n information which was used to derive an encoded thumb movement s i g n a l . A sketch of the apparatus i s given i n F i g . 4.1a. The duration of the ramped < extension was approximately one second long. After subjects completed the extension and had returned t h e i r thumb to the s t a r t i n g p o s i t i o n , v i s u a l feedback, v i a d i f f e r e n t colored l i g h t s , was presented to them i n d i c a t i n g whether they had h i t , overshot, or undershot the target p o s i t i o n . During thumb movements subjects were asked to f i x a t e on the v i s u a l feedback area i n an attempt to minimize eye movements and to prevent subjects from looking at t h e i r thumb. Aft e r subjects were given some pr a c t i c e t r i a l s , they were asked to carry out f i f t y self-paced r e p e t i t i o n s Thumb br a c k e t Lever, i g h t L i n e a r p o t e n t iometer Lever s u p p o r t f o r s t a r t i n g p o s i t i o n l a R i g h t hand g r i p s dowel wit h thumb p l a c e d under the thumb b r a c k e t and p e r p e n d i c u l a r t o the l e v e r 43 with t h e i r r i g h t thumb. On each t r i a l , a c q u i s i t i o n of the EEG started 1 second before thumb movement onset (determined by monitoring the encoded thumb movement signal) and then continued for 4.5 more seconds at which time the feedback was presented. The a c q u i s i t i o n was then halted one second a f t e r presenting the feedback which res u l t e d i n a t o t a l epoch length of 6.5 seconds. The above active task was based on a previous study into motor po t e n t i a l s c a r r i e d out by Grunewald and Grunewald-Zuberbier [7] . In t h e i r analysis they u t i l i z e d conventional averaging techniques and t h e i r grand averages across seven subjects, 35 t r i a l s each, are given i n Figure 4.1b. For the i d l e task the subjects were kept i n the same phys i c a l s i t u a t i o n as i n the active task but i n t h i s case they were not performing any thumb movements. Twenty epochs of EEG, 6.5 seconds i n length, were c o l l e c t e d from the subjects under t h i s condition. A f t e r each epoch was c o l l e c t e d , the feedback l i g h t s flashed to ind i c a t e to the subject that they had 10 seconds to relax before the onset of the next epoch. Af t e r taking a short pause the onus was on the subject to f i x a t e t h e i r eyes on the feedback l i g h t s i n preparation for the onset of the next epoch. The EEG signals were recorded from the scalp using s i l v e r / s i l v e r chloride electrodes. This type of electrode possess the most appropriate c h a r a c t e r i s t i c s for EEG recording i n terms of low p o t e n t i a l d i f f e r e n c e s , long time constants and low resistance between the e l e c t r o l y t e and the metal surface of the electrode. The electrodes were "cupped" shaped with a small hole i n the top and they were f i r m l y attached to the scalp around the rim with the use of c o l l o d i o n . E l e c t r o l y t e j e l l y was i n j e c t e d into the a i r space under the "cup" of the electrode which provided a good e l e c t r i c a l connection between the scalp and the electrode. T y p i c a l l y the electrode impedance LEFT INDEX FIN6ER RIGHT INDEX FINGER A B CD E F A B CD E F Ramp p o s i t i o n i n g movements of the l e f t and r i g h t index f i n g e r . Grand averages across seven r i g h t handed subjects. Recordings of force, r e c t i f i e d EMG, po s i t i o n , v e r t i c a l EOG, and slow p o t e n t i a l s h i f t s i n l e f t and r i g h t precentral (C3' f C4') and postcentral (C3'', C4**) EEG (1cm a n t e r i o r and 2 cm po s t e r i o r to C3, C4 p o s i t i o n s ) . Figure 4.1b 45 between the scalp and the reference electrode was approximately 2500 ohms. The sign a l was i n i t i a l l y recorded from three standard i n t e r n a t i o n a l 10/20 system electrode s i t e s Cz, C3 and C4 (See Jasper [49]). The s i g n a l from each of these electrodes was referenced to linked ear lobes. A bi - p o l a r EOG signal and the corresponding encoded thumb movements were also recorded. The EOG electrodes were placed on the supra o r b i t a l ridge and the external canthi of the r i g h t eye. The EOG s i g n a l was u t i l i z e d i n a very conservative a r t i f a c t r e j e c t i o n c r i t e r i o n which rejected any EEG epoch that had a corresponding EOG s i g n a l that at anytime during the epoch fluctuated above or below baseline by more than a given threshold value. This value was nominally set at 17 microvolts. As w e l l , any EEG epoch that was not rejected by the above c r i t e r i o n but contained peak values that exceeded baseline by 43 microvolts, t y p i c a l l y due to f a c i a l EMG, were also rejected as artifact-contaminated. In addition, for the active case, t r i a l s not containing a reasonable thumb movement were also rejected. Ultimately, based on the above s e l e c t i o n c r i t e r i a , 15 t r i a l s from both the active and i d l e cases were u t i l i z e d from each of the four subjects. The EEG signal that was used for a l l the subsequent s i g n a l analysis was taken from the C3 electrode s i t e which was c o n t r a l a t e r a l to the thumb movement. In a l l cases the EEG and EOG signals were i n i t i a l l y amplified by a Beckman 711 polygraph using an analogue lowpass f i l t e r with a -3dB point at 100Hz (20dB per decade r o l l - o f f ) and a highpass analogue f i l t e r with a time constant of 14.7 seconds. The signals for each epoch were d i g i t i z e d i n r e a l time at a rate of 1024 samples per second and were stored on a hard disk. Before any s i g n a l analysis schemes were applied, the EEG signals were preprocessed by a 201-point phaseless d i g i t a l lowpass f i l t e r which had a cutoff frequency of 29Hz (-3dB p o i n t ) , a t r a n s i t i o n width of 3Hz (-24dB at 46 32Hz), and a minimum stopband attenuation of -27dB. It i s generally agreed [50] that almost a l l the power i n the normal EEG i s between zero and t h i r t y Hertz. Therefore, with the above d i g i t a l f i l t e r the data was resampled at the r e l a t i v e l y low rate of 64 Hz which i s desirable because as the sample rate increases, there i s a corresponding need to increase the AR model order since the time dependency i s spread over a greater number of sample points. Generally, one wants to make the best possible trade-off between using a s u f f i c i e n t l y high sample rate that w i l l allow for the accurate representation of the highest frequency of i n t e r e s t and yet within that framework keep i t as low as possible so that the required model order i s minimized. 4.2 Neurological Premise The concept that event r e l a t e d information i s contained i n EEG time series o u t l i e r s i s based on the following model of summation of e l e c t r i c a l b r ain a c t i v i t y at a given point on the scalp. Under i d l e conditions the ongoing e l e c t r i c a l a c t i v i t y that sums, s p a t i a l l y and temporally, at a given point on the scalp can be modeled as an o v e r a l l ongoing process as "observed" from that point on the scalp during* a p a r t i c u l a r time i n t e r v a l . When event r e l a t e d p o t e n t i a l s , such as motor r e l a t e d p o t e n t i a l s , are generated by a unique a d d i t i o n a l process they are "added" into the p r e - e x i s t i n g ongoing process and would appear as additive o u t l i e r content when considered from the point of view of the ongoing pre - e x i s t i n g process. Therefore, i f one could d i s t i n g u i s h o u t l i e r points from pre - e x i s t i n g process points i n the s i n g l e t r i a l a ctive EEG time seri e s then these o u t l i e r points could be used to provide information about event r e l a t e d p o t e n t i a l s on a single t r i a l b a s i s . 47 The underlying p r i n c i p l e of the s i n g l e t r i a l processing scheme i s to generate an AR model of the active EEG s i g n a l using a robust parameter estimation method that w i l l represent the ongoing, underlying process by down-weighting unusual data points (see Section 3.3). This estimated model i s then used i n a robust s i g n a l estimator (see Section 4.3) which produces an estimated s i g n a l of the ongoing, underlying EEG process. The difference between the o r i g i n a l measured s i g n a l and the estimated s i g n a l i s considered to be additive o u t l i e r content (see Section 4.4.1). The o u t l i e r content i s then processed (see Section 4.4.2) to produce waveform patterns that provide si n g l e t r i a l event r e l a t e d information. 4.3 Signal Cleaning Process An approach to detect o u t l i e r points i n a time ser i e s was proposed by Martin and Thomson [45], They used t h i s system not to study the character, or information contained therein, of the o u t l i e r s but rather to produce a "cleaned" time series which was used as part of a process that produced robust s p e c t r a l estimates. This cleaning process i s based on a " r o b u s t i f i e d " Kalman s i g n a l estimator. The objective of the cleaning process i s to provide an estimate of the o r i g i n a l s i g n a l without the AO content. It r e l i e s on a e s t i m a t e d pth order AR model of the process x^ as g i v e n i n the AO model (3.37) which for convenience i s repeated below = x. + v. l l 4.1 This AR process written i n state v a r i a b l e form i s given by x. — I = $ X. - i - 1 + u. — 1 4.2 48 T T where x. = (x. ,x. x. u. = (e.,0,0 0) , and - i I i - i i-p+1 - i i ' ' ' ' • a l a2 1 0 0 1 0 0 0 0 . 1 0 4.3 Note t h a t , g i v e n the above d e f i n i t i o n , the state x^ i s equal to the current value of x. and past values of x. up<to x. ,,. l I r i-p+1 "Robust" estimates of the state x^ are calculated r e c u r s i v e l y with the following expression [45] m. y. - x x. = <J>x. , + — s. \|) ( - i - i - l s? l y s. l i-1- ) 4.4 where i s the f i r s t column of the pxp m a t r i x which i s r e c u r s i v e l y c a l culated as follows M l + l l x 4.5 and y. - x . , a m.m. r> w / i ~ i - l \ - i - i P. = M. - w ( ) s. 4.6 where Q i s a pxp matrix with a l l zero entries except the f i r s t element which i s equal to a robust estimate of the variance of the re s i d u a l sequence, i . e . Q( l , l ) = S g 2 and s^ i s a time varying scale defined by s±* = M i ( l , l ) 4.7 49 The cleaned data at time i , w i l l then be the f i r s t element of the estimated state x., which i s - i y. = x. 4.8 J 1 C 1 In other words, y^ i s an estimate of the process without the influence of the a d d i t i v e o u t l i e r s , v^. Note that with the s c a l i n g given 4.7 and given that there i s no influence from the p s i function, as would be the desired r e s u l t when there i s no o u t l i e r content, x. = y.. Hampel's three part redescending p s i function [43] was used i n Equation 4.4 and i t i s given as follows: Itl 0 <. | t | < a a <. | t | < b ii(t) = sign(t) { 4.9 c-b b * Itl < c 0 c £ | t | where a, b and c are tuning parameters. In a s i m i l a r fashion as i n the case of GM-estimation the w function i s of the form w(t) = f ( t ) / t 4.10 The cleaning process described above was u t i l i z e d as part of the procedure to obtain GM1 and GM2 AR parameter estimates. A GM1 estimate i s based on the estimated cleaned time ser i e s and a GM2 estimate i s a further i t e r a t i o n where the parameters from GM1 are used i n the cleaner to provide a t h e o r e t i c a l l y improved estimate of the cleaned time s e r i e s . Figure 4.2 provides a block diagram of the procedure to obtain GM, GM1, and GM2 parameter estimates. Through the simulation studies described i n Section 3.4 i t was found that the best performance i n terms of minimizing the MSE was obtained by s e t t i n g the tuning parameters i n Equation 4.9 as given below 50 ROBUST AR PARAMETER ESTIMATION y(k) 3 / GM ESTIMATION y(k) a CLEANER GM ESTIMATE \1 / 9c(k) GM ESTIMATION y(k) .A a CLEANER GM1 ESTIMATE GM ESTIMATION Figure 4.2 GM2 ESTIMATE 51 tuning parameter a 1.8 b 2.2 c 3.0 4.4 Extracting and Processing O u t l i e r Information 4.4.1 Extracting O u t l i e r Information The o u t l i e r e xtraction process for the EEG data i s accomplished by taking the epoch, 6.5 seconds long, and d i v i d i n g i t into 1.5 second segments with each segment overlapped by .75 seconds. Each segment i s modeled at the order expected for i d l e task EEG (12-14) and hence, reducing the a b i l i t y of the model to account for active task information i n the EEG (see Section 4.5). The s i g n a l from each segment i s then cleaned using the estimated model parameters i n the cleaner described i n Section 4.3. The o u t l i e r s are then calculated by taking the diff e r e n c e between the o r i g i n a l and cleaned signals (see Figure 4.3). * The o u t l i e r e xtraction process was i n i t i a l l y tested by applying i t to Gaussian simulated 12th order AR data which contained 10% "patchy" (correlated) additive o u t l i e r contamination. This t e s t confirmed that the extraction process had some d i s t i n c t a b i l i t y to recover the o u t l i e r content from the simulated s i g n a l . Patchy contamination was used since i t was expected that i n the case of r e a l EEG the additive event r e l a t e d p o t e n t i a l s would be correlated. As suggested by Martin and Zeh [51] , the correlated v a l u e s f o r v. were generated by d i v i d i n g the segment i n t o equal halves. P R O C E S S TD E X T R A C T O U T L I E R P O I N T S yCk) 7> P R R R M E T E R E S T I M R T I O N A a yCk) C L E R N E R PtCk) yCkJ-y^k) yCk) O U T L I E R S F i g u r e 4.3 53 Immediately f o l l o w i n g the f i r s t non-zero v^ i n each h a l f the r e s t of the non-zero v^'s were grouped together. These grouped v^'s were used to produce corr e l a t e d v_^'s v i a the following expression v? = 0v? . + (1 - Q 2 ) 1 ' 2 v . 4.11 l i - l l where a value of 0 = 0.6 was used i n these simulation t e s t s . This procedure r e s u l t s i n a correlated o u t l i e r series which has roughly the same variance as i n the independent case [51]. It was found i n t h i s a p p l i c a t i o n that the tuning parameters for the p s i function given i n Equation 4.9 needed to be set such that i t provided a stronger influence than i n the parameter estimation a p p l i c a t i o n discussed i n Section 4.3. By t r i a l and error the best performance of the extraction process was obtained by s e t t i n g the tuning parameters as given below tuning parameter * a 1.0 b 1.2 c 1.8 Some example r e s u l t s from these t e s t s , which q u a l i t a t i v e l y demonstrate the p o t e n t i a l performance of the ex t r a c t i o n process, are shown i n Figure 4.4. c Each example contains the same randomly generated v_^  contamination which was used to contaminate d i f f e r e n t Gaussian AR sequences x^. The broken l i n e i n each p l o t r e p r e s e n t s the actual v^ values and the s o l i d l i n e represents the o u t l i e r content that was extracted v i a the the o u t l i e r e x t r a c t i o n process. Figures 4.4a through 4.4c show three d i f f e r e n t examples of using GM2 and LSQ parameters i n the cleaning process. Figures 4.4d and 4.4e are the second and C O R R E L A T E D O U T L I E R D E T E C T I O N : flR-12 A D D I T I V E O U T L I E R 10% V O R . = 2 . 0 N P T S = 1 0 0 O F F S E T - I 5.892 ••.983 4.073 3.163 2.234 1.344 0.433 -0 .475 — 1.384 -2.294 -3 .203 — 4.113 — S.022 h i • v i o OM2 Parameter Estimates 2 0 3 0 4 -0 Data Points 5 0 7 0 e o 9 0 4.831 j 4.0Z9 j 3.206 ! . 2.383 j 1.560 j 0.737 I — 0.08S j —0.908 j — 1.731 —2.354 ' - 3 . 377 ! — 4.199 — 3.022 b -• i V i o 2 0 3 0 Data Points LSQ Parameter Estimates 5 0 6 . 0 7 0 8 0 ~ 9 0 Figure 4.4a £ x t r a o t « d Outltor* C O R R E L A T E D O U T L I E R D E T E C T I O N : A R - 1 2 A D D I T I V E O U T L I E R 107. V A R . = 2 . 0 N P T S = 1 0 0 O F F S E T - 4 3.271 4,244 i 3.2X7 2.190 1.163 0.133 - 0 . 8 9 2 — 1,919 — 2.946 - 3 . 9 7 3 — 3.OOO -6 .Q28 l O 20 30 Data Points A-O GH2 Parameter Estimates S O I— 6.0 7"b 8 0 9 f J 4.851 4.029 3.20& 2.383 1.3&0 0.737 -o.oes - 0 . 9 0 B — 1.731 — 2.334 — 3.377 — 4.199 — 3.022 l O I 20 SO Data Points A LSQ Parameter Estimates _i_ 6,o > ° " " 8 0 9 0 _ i - J _ - -E x t r a o t a d O u t l i « n Aotual OutlMr* Figure 4.4b C O R R E L A T E D OUTLIER DETECTION: AR-12 ADDIT IVE OUTLIER 107. VAR.=2.0 NPTS=100 O F F S E T = 8 7 , 0 0 0 6 . 0 0 0 5 . 0 0 0 4 . 0 0 0 3 . 0 0 0 2 . 0 0 0 v 1 . 0 0 0 • ° 0 . 0 0 0 ^ — 1 . 0 0 0 = — 2 . 0 0 0 " - 3 . 0 0 0 = « - 4 , 0 0 0 — 5 . 0 0 0 — 6 . 0 0 0 — 7 . 0 0 0 1.— I o 1 0 .... - U . 2 0 3 0 Data Points 4-0 t 5 0 GM2 Parameter Estimates 7 0 s o s>o 7 . 0 0 0 6 . 0 0 0 s . 0 0 0 4 . 0 0 0 3 . 0 0 0 2 . 0 0 0 1 . 0 0 0 0 . 0 0 0 - 1 . 0 0 0 - 2 . O O O - 3 . 0 0 0 - 4 . 0 0 0 - s . 0 0 0 - 6 . 0 0 0 - 7 . 0 0 0 - v r - a . LSQ Parameter Estimates l O : o 3 0 Data Points 4-0 S O 6 . 0 70 8 0 9 0 E x t r a c t e d O u t l i a n ftotuol O u t l i e r s Figure 4.4c C O R R E L A T E D O U T L I E R D E T E C T I O N : A R - 1 2 A D D I T I V E O U T L I E R 107. V A R . = 2 . 0 N P T S = 1 0 0 O F F S B T = 4 6.134 5.132 4.10? 3.087 2.064 1.042 0.019 -1.003 -2.026 •3.048 -4.071 •3.093 -6.116 _L GM1 Parameter Estimates I O 2 0 3 0 Data Points 4-0 5 0 7 0 S O 9 0 4.286 3.431 2.573 1.713 0.858 O.OOO -0.838 -1.713 -2.573 -3.431 -4.288 -5.146 -6.003 ... .1 I O 3 c F " 4 o Data Points A, A s o GM Parameter Estimates 6,0 J 7 0 8 0 L. E x t r a o t « d O u t l i s r s A o t u a l O u t l i e r s — i Figure 4.4d C O R R E L A T E D OUTLIER DETECTION: ADDIT IVE OUTLIER 107. VAR.=2.0 NPTS-100 O F F S E T - - 8 7 . O O O 6 , O O O 5 . O O O 4 , 0 0 0 3 . O O O 2 . O O O 1 , 0 0 0 O . O O O - l . O O O - 2 . O O O - 3 . O O O - 4 . O O O - 5 . O O O - 6 . O O O - 7 . O O O tz L . A . l O 2 0 3 0 4-0 Data Points S O GM1 Parameter Estimates 6.0 7 0 8 0 9 0 7 . O O O 6 . O O O 5 . O O O 4 . O O O 3 . O O O 2 . O O O l . O O O O . O O O - l . O O O - 2 . O O O - 3 . O O O - 4 . O O O — 3 . O O O - 6 . O O O - 7 . O O O l O 2 b 3 0 4-0 Data Points GM Parameter Estimates 5 0 6 0 7 0 8 0 9»0 E x t r a c t e d Outlier* Figure 4.4e Aotual Outliers 59 t h i r d examples repeated using GM and GM1 parameter estimation. It i s clear from these examples that the process performs better with GM2 estimates than with GM estimates and much better than with LSQ estimates. The performance of GM2 and GM1 are quite s i m i l a r with perhaps some subtle improvements i n the case of GM2. Since these t e s t s revealed three c l e a r l y discernable jumps i n performance i n using LSQ, GM, and GM2 parameter estimates, i t was decided that subsequent studies using o u t l i e r detection i n t h i s t h e s i s work would be r e s t r i c t e d to those three estimation methods. 4.4.2 Processing O u t l i e r Information Processing the extracted o u t l i e r information i s accomplished by taking the o u t l i e r content from each EEG segment and averaging i t together with the o u t l i e r content from the corresponding overlapping segment. This r e s u l t s i n an o u t l i e r pattern spanning the whole 6.5 second epoch. The o u t l i e r pattern i s then smoothed by convolving i t with a 16 point tapered smoothing window which i s based on a minimum-bias sp e c t r a l window suggested by Papoulis [52]. It i s given by ,,,,, 1 C , 1 + COS(2TT k/16) W(k) = 16TT2 4.12 [16(2ir k/16) 2 - I T 2 ] 2 where k = 0, +/-1, +/-2 +/-16 The r e s u l t i n g smoothed pattern constitutes the output waveform of the sin g l e t r i a l processing method. 60 A. 5 AR Spectral Analysis Preliminary studies involving AR spe c t r a l analysis were useful i n providing some measure of the a b i l i t y of the AR model to represent the EEG s i g n a l . As w e l l , these studies were instrumental i n e s t a b l i s h i n g appropriate EEG segment lengths and a procedure for the s e l e c t i o n of the AR model order. I t can be shown [23] t h a t t h e AR s p e c t r a l e s t i m a t e i s S (f) S(f) = A.13 P II - I a k e x p - ( ^ ) | 2 k=l s where S (f) i s the power spectrum of the re s i d u a l sequence e. and f i s the 6 I S sample frequency. Since the term S e ( f ) applies to the residuals which are i n theory white, the r e s u l t i n g power density function of the residuals should be f l a t and t h e r e f o r e S e ( f ) w i l l be a constant independent of frequency. Id e a l l y , the value of t h i s constant (noting that the mean of the residuals i s zero) w i l l be proportional to the variance of the residuals [A6]. Hence, the f i n a l expression for the conventional AR spe c t r a l estimate i s obtained by r e p l a c i n g S (f) i n A.13 with s 2 / f , where s 2 i s an estimate of the variance r ° e e s e of the r e s i d u a l s and the 1/f term i s included i n the numerator so that the s true power of the corresponding analogue s i g n a l w i l l be represented [23]. The EEG si g n a l c h a r a c t e r i s t i c s from subjects, p a r t i c u l a r l y during h i g h l y active mental states, are changing r e l a t i v e l y quickly. Single t r i a l AR sp e c t r a l estimates from adjacent one second segments demonstrated that considerable change i n si g n a l c h a r a c t e r i s t i c s could occur over t h i s span of two seconds. An example of t h i s i s provided i n Figure A.5. I t contains four consecutive AR spe c t r a l p l o t s , each derived from a one second segment of 61 Figure 4.5 AR spectral estimates of one second segments of active trial EEG consecutively offset by a third of a second. 62 active EEG and each o f f s e t by 0.333 seconds ( t o t a l span of * o seconds). As was discussed i n Section 3.2, i t was found that a p r a c t i c a l lower bound on segment length, from a parameter estimation point of view, was approximately one second. These findings u l t i m a t e l y lead to the u t i l i z a t i o n , as noted i n Section 4.4.1, of a 1.5 second segment length with an o f f s e t of 0.75 seconds i n the si n g l e t r i a l processing method. This was an attempt to trade o f f the need for short segments because of the r e l a t i v e l y r a pid changing signal c h a r a c t e r i s t i c s with the desire to r a i s e the segment length above the lower bound for purposes of improving the parameter estimation e f f i c a c y . It was found that s e l e c t i n g the model order v i a conventional methods such as Akaike's Information C r i t e r i a (AIC) does not work well with these short segments [53]. Conclusions were s i m i l a r to Jansen [32] i n that the s e l e c t i o n of an appropriate model order requires some t r i a l and error and, i f possi b l e , some a - p r i o r i knowledge of expected r e s u l t s . I t was found useful to t r y a number of orders within a reasonable range (for a sample rate of 64Hz, somewhere between 8 to 25), following the trend of the estimate as the model order was increased. Features were i d e n t i f i e d that seemed reasonable based on both the a - p r i o r i knowledge of the condition under which the EEG was c o l l e c t e d and a conventional FFT based estimate. The order was sequentially increased, expecting the features to become better defined, u n t i l spurious peaks began to occur. The appropriate model order was then selected to be two or three below that value. T y p i c a l l y , model orders were selected i n the range of 12 to 14 from subjects during the i d l e task and i n the range of 18 to 22 from subjects during the active task. Figures 4.6 and 4.7 provide example AR s p e c t r a l p l o t s to demonstrate t h i s model order s e l e c t i o n procedure for the i d l e and active task r e s p e c t i v e l y . 63 Figure 4.6 Progression of AR spectral estimates with increasing order using BEG data from an example idle trial, a) conventional FFT b) model order 8 c) model order 10 d) model order 12 e) model order 14 f) model order 16 64 LO • as •e a *. a* e 0)04 0 ao " 1 2 '' i«'" 'to'' '2C F r a q u a n o y (Hz) Co> 8 s i r F r « q u « n o y (Hz) Cd) 24 LO — 12 I * 2 0 24 F r « q u « n o y (Hs) (•) LO • as _ 4 . a« >a4U a2 ao 12 16 20 r r a q u w i o y CHx) t f ) 24 Figure 4.7 Progression of AR spectral estimates with Increasing order using BEG data from an example active task, a) conventional FFT b) model order 8 c) model order 10 d) model order 12 e) model order 14 f) model order 16 65 In the de r i v a t i o n of the autoregressive s p e c t r a l estimate an al t e r n a t i v e to assuming that the residuals w i l l be p e r f e c t l y white i n the c a l c u l a t i o n of S g ( f ) would be to estimate that quantity with a conventional FFT based estimate. The r e s i d u a l s i g n a l can be thought of as a whitened si g n a l because the information that can be represented by an AR model has been subtracted r e s u l t i n g i n a si g n a l with a much f l a t t e r spectrum. When the FFT i s applied to t h i s prewhitened signal the inherent drawback of leakage i s gre a t l y reduced. A p p l i c a t i o n of conventional leakage c o n t r o l , such as Blackman windowing, serves to further reduce t h i s problem. The prewhitened AR estimation method, therefore, combines the spe c t r a l information from both the AR model and the re s i d u a l FFT spe c t r a l estimate. I t i s given by [45] S(f) = — r 4.14 1 - I a, exp - ( l ^ f ) k=l K f s where S.T (f) i s a s p e c t r a l estimate of the r e s i d u a l sequence e. using a Ne l conventional FFT method. Some i n s i g h t into the a b i l i t y of the AR model to represent short segments of EEG was gained by pursuing studies using prewhitened AR spe c t r a l estimates. These studies demonstrated that when an appropriate model order was u t i l i z e d the conventional AR spe c t r a l estimates were reasonably good compared to the prewhitened AR estimate which makes use of information retained i n the residuals (see Birch et a l . [53]). This indicates that the AR model, although not perfect, does represent much of the information contained i n a short segment of EEG. An example of both a conventional and a prewhitened 12th order AR spe c t r a l estimate of i d l e task EEG i s given i n Figure 4.8. 66 RR SPECTRAL ESTIMATION OF IDLE TASK EEG Frequency Figure 4.6 67 A.6 Applying O u t l i e r Processing to Single T r i a l EEG A.6.1 Comparison of Segmented Cleaned Active, O r i g i n a l Active and O r i g i n a l Idle Signals It would be predicted, given the above neurological premise, that the o r i g i n a l i d l e s i g n a l and the cleaned active s i g n a l should have l i t t l e or no evidence of motor p o t e n t i a l a c t i v i t y whereas the o r i g i n a l a c t i v e s i g n a l should contain evidence of motor p o t e n t i a l a c t i v i t y . Figure A.9 provides two sets of p l o t s : one with N=6 t r i a l s the other with N=15 t r i a l s . Each set contains p l o t s of conventionally averaged cleaned a c t i v e , o r i g i n a l a c t i v e , and o r i g i n a l i d l e s i g n a l s . Motor p o t e n t i a l a c t i v i t y i n the active case should occur, approximately, during the f i r s t three seconds of the epoch, noting that the actual thumb movement began one second into the epoch. These pl o t s demonstrate that the conventional averaging technique reveals some d i s t i n c t motor a c t i v i t y i n the o r i g i n a l active case (raised l e v e l of p o s i t i v i t y i n the averaged s i g n a l during the f i r s t three seconds with a peak at about two seconds). However, i n the cleaned active and o r i g i n a l i d l e cases the averaging does not reveal any d i s t i n c t motor a c t i v i t y . Hence, the above p r e d i c t i o n i s s u b s t a n t i a l l y borne out. The strong negative peak i n the N=6 active task p l o t at about 6 seconds i s the v i s u a l evoked response to the feedback l i g h t and i s not due to motor p o t e n t i a l a c t i v i t y . AVERAGE OF SEGMENTED ORIGINAL AND CLEANED SUBJECT*! N=6 ORIGINAL. A C T I V E Tim* <8*oonds) lO 0 * 2 O —2 — * - lO C L E A N E D A C T I V E — C3M2 O.O 0.5 l.O 1.5 Z.O ~2.5 3.0 3.5 4..'o 4.3" S.O 5.5 ~6.C Tim* ( 8 « e o n d s ) DRIGINflL I D L E —I 1 I I i j i i I i i . , O.O O.S l.O 1.5 2.0 Z.5 3.0 5.5 4.0 4.5 S.O "5.S fe.O Tim* (Saoonds) Figure 4.9a 69 AVERAGE OF SEGMENTED ORIGINAL AND CLEANED EEG SUBJECTS! N=15 O R I G I N A L A C T I V E O . O O . S l . O l . S 2 . 0 2 . S T i m * ( S a o o n d a ) S . O 3 . S 4 . 0 C L E A N E D A C T I V E - C3M2 O . O O . S l . O 1 . 5 2 . 0 2 . S T u n a ( 8 « a o n d « ) 3 . 0 3 . 5 4 . 0 ORIGINAL. I D L E F i g u r e 4.9b 70 4.6.2 Examples of Single T r i a l O u t l i e r Patterns and Single T r i a l Raw EEG To q u a l i t a t i v e l y demonstrate the r e s u l t s of the s i n g l e t r i a l processing method, four example p l o t s of the single t r i a l o u t l i e r patterns using GM2 model parameters paired with the corresponding raw EEG are provided i n Figure 4.10. Note the s i g n i f i c a n t amount of information that i s i n the o u t l i e r patterns which can not be e a s i l y seen i n the raw EEG s i g n a l s . Results provided i n the following sections demonstrate that the information i n these o u t l i e r patterns i s r e l a t e d to the thumb movements. However, at t h i s point i t i s i n t e r e s t i n g to note the many s i m i l a r i t i e s of these s i n g l e t r i a l patterns with the grand average waveforms from the Grunewald study c i t e d i n Section 4.1. 4.6.3 Comparison of Averaged Active O u t l i e r Patterns, Averaged Idle O u t l i e r Patterns and the Conventional Average of Active EEG To demonstrate that there i s some strong consistency i n the active case o u t l i e r patterns and very l i t t l e consistency i n the i d l e case o u t l i e r patterns, the p l o t s i n Figure 4.11 have been provided. These pl o t s contain the averaged o u t l i e r patterns for N=6 and N=15 using GM2 parameter estimates. As w e l l , for comparison purposes, a p l o t of the conventional average for the active case i s also included i n t h i s f i g u r e . The fact that the average active case patterns maintain a general shape s i m i l a r to the s i n g l e t r i a l patterns, strongly indicates that there i s information r e l a t e d to the thumb movement that i s consistent from t r i a l to t r i a l . On the other hand, the f a c t that the average reduces i n magnitude and EXAMPLE VRVEEORMS Subject #1 - Trial #3 O U T L I E R P A T T E R N 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 5.5 6.0 Time (Seconds) R A W E E G 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 5.5 6.0 Time (Seconds) Figure 4.10a EXAMPLE WflVEEORMS Subject #1 - Trial #6 O U T L I E R P A T T E R N 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 5.5 6.0 Time (Seconds) R A W E E G 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 5.5 6.0 Time (Seconds) Figure 4.10b EXAMPLE WAVEFORMS Subject #1 - Trial #15 O U T L I E R P A T T E R N 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 5.5 6.0 Time (Seconds) R A W E E G 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 5.5 6.0 Time (Seconds) Figure 4.10c EXAMPLE VAVEEORMS Subject #1 - Trial #35 O U T L I E R P A T T E R N 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 5.5 6.0 Time (Seconds) 172 111 50 -11 -72 -133 -193 -254 -315 -376 R A W E E G _L J_ _L 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 5.5 Time (Seconds) Figure 4.10d 6.0 75 AVERAGE WAVEFORMS N=6 C O N V E N T I O N A L A V E R A G E — i I 1 1 i i I i i I J i i O . O O . S l . O 1 . 5 2 . 0 2 . 5 3 . 0 3 . S 4 . 0 4 . 3 S . O 5 . 5 6 . 0 T i m * ( S a a o n d a ) A C T I V E O U T L I E R — GM2 — i i 1 l i i 1 i i i i i i _ O . O O . S l . O 1 . 5 2 . 0 2 . 5 3 . 0 3 . 5 4 . 0 4 . 5 S . O S . S e>.0 T B T M I CSakoondart I D L E O U T L I E R — GM2 — i 1 1 1 1 i i I I i i . I I O . O O . S l . O l . S 2 . 0 2 . 5 3 . 0 3 . 5 4 . 0 4 . 5 5 . 0 5 . S 6 . 0 T t m a C8«)Oondaf) Figure 4.11a 76 AVERAGE WAVEEORMS N=15 CONVENTIONAL. A V E R A G E Tim* ( S a o o n d s ) O^ O 5^ 3LO l £ zfe z T g 37o sts 4^3 Tim* W«nond«) 77 subtle features become less pronounced as N i s increased, indicates that there i s a s i g n i f i c a n t amount of uniqueness i n i n d i v i d u a l t r i a l s that i s l o s t as many t r i a l s are averaged together. This uniqueness i s c e r t a i n l y i n part due to the variance i n thumb movements from t r i a l to t r i a l which can be seen c l e a r l y i n Section A.7. I t i s also expected that a d d i t i o n a l t r i a l by t r i a l uniqueness i s due to cognitive factors such as the mental i n t e n s i t y with which the subject c a r r i e d out the task. The average of the i d l e task o u t l i e r patterns c l e a r l y demonstrates that there i s no s i g n i f i c a n t information that i s being r e i n f o r c e d across i d l e task t r i a l s . The conventional average of active t r i a l s shows that with N's of 6 and 15 the motor p o t e n t i a l information i s quite l i m i t e d and the "smearing" e f f e c t of event r e l a t e d information that i s discussed above for the active case o u t l i e r patterns would also be occurring i n these conventional averages. Hence, with the conventional averaging method, even with much greater N's as i n the case of the Grunewald study (see Section A . l ) , the information obtained w i l l be l i m i t e d to that which has remained r e l a t i v e l y consistent across the t r i a l s . A.6.A LSQ Active O u t l i e r Patterns Degrading with Higher Model Orders It would also be expected, based on the neurological premise, that the single t r i a l processing method would perform best when the AR model order was selected to best f i t the i d l e case. As the model order i s increased the AR model would be expected to gain some improved a b i l i t y to represent the motor r e l a t e d a c t i v i t y i n the active task EEG. Hence, the performance of the sin g l e t r i a l method should begin to degrade since the cleaning process, which u t i l i z e s the higher order AR model, would lose some of i t s effectiveness i n 78 detecting motor r e l a t e d o u t l i e r s . A p a i r of averaged active o u t l i e r pattern p l o t s using LSQ parameters for model order 12 (generally appropriate for the i d l e case) and model order 22 (generally appropriate for the active case) are shown i n Figure 4.12. These p l o t s demonstrate that the performance does degrade, i n terms of both the amplitude and the d e t a i l of features i n the averaged o u t l i e r pattern, when the model order i s better matched to the active case. 4.7 S t a t i s t i c a l Analysis of Features i n the O u t l i e r Patterns The set of 15 active t r i a l o u t l i e r patterns superimposed with the corresponding encoded thumb movements from Subject #1 are given i n Figure 4.13. The patterns, although unique from t r i a l to t r i a l , do seem to posses a generally consistent waveform which contains features that appear r e l a t e d to events i n the thumb movements. S t a t i s t i c a l analysis was c a r r i e d out to determine i f features i n the i n d i v i d u a l thumb movements are r e l a t e d to features i n the corresponding o u t l i e r patterns. Two features i n the thumb movement and three features i n the o u t l i e r pattern were u t i l i z e d i n the s t a t i s t i c a l a n a l y s i s . The features are described below and are shown i n Figure 4.14. Feature 1: Time from epoch onset to the point when the thumb movement f i r s t reaches the "on target" p o s i t i o n . Feature 2: Time from epoch onset to the point when the thumb movement f i r s t leaves the "on target" p o s i t i o n . COMPARISON OF N=15 AVERAGED OUTLIER PATTERNS L S Q O R D E R 12 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 Time (Seconds) L S Q O R D E R 22 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 Time (Seconds) Figure 4.12 S I N G L E T R I A L GM2 O U T L I E R P A T T E R N S W I T H C O R R E S P O N D I N G T H U M B M O V E M E N T TRIAL #3 0.0 0.5 1.0 145 2.0 245 3.0 345 4.0 4.5 Tins (Stoonds) 0.0 045 LO 145 2.0 245 3.0 34J 4.0 445 TMM Oaoonda) TRIAL #6 0J0 0J5 1.0 1.5 2.0 245 3.0 345 4.0 445 TIM OMMMS) TRIAL #8 OX) 045 145 2.0 24) 3X1 345 AX) 445 oo o Figure 4.13a R . l o t l v . d a o n l t u d . R . l o t l v . N a g n 11 u d • T8 R a l a t i v * M a g n i t u d e B « l o t l v « N a g n l t u d a Z8 \?Tr£9L£«IRIftL G M 2 OUTLIER PATTERNS WITH CORRESPONDING THUMB MOVEMENT TRIAL #41 TRIAL #44 0.0 0.5 1.0 1.3 2.0 245 3.0 345 44) 4J TfeM tttootxfc) 00) 041 1.0 14$ 24) 245 34) 34S 4.0 445 TIM WtooBdrt 84 SINGLE TRIAL FEATURE DEFINITION SAMPLE THUMB MOVEMENT WITH STANDARD GM2 OUTLIER PATTERN FOR SUBJECT #1 0,5 1,0 1,5 2.0 2.5 3.0 3,5 4.0 Tftne (Seconds) Figure 4.14 85 Feature 3: Time from epoch onset to the f i r s t dominant (greatest amplitude) p o s i t i v e peakn the o u t l i e r pattern. Feature A: Time from epoch onset to the f i r s t negative peak i n the o u t l i e r pattern a f t e r feature 3, that has a minimum of 5 units magnitude peak-to-trough d i f f e r e n c e . Feature 5: Time from epoch onset to the next p o s i t i v e peak i n the o u t l i e r pattern a f t e r feature A, that has a minimum of 20 units magnitude peak-to-trough diffe r e n c e on both sides of the peak. There was an expectation r e s u l t i n g from the e a r l i e r conventional study by Grunewald and Grunewald-Zuberbier [7] and from observations taken from Figure A.13 that feature 1 would be p a r t i c u l a r l y r e l a t e d to features 3 and A whereas feature 2 would be p a r t i c u l a r l y r e l a t e d to feature 5. The sample c o r r e l a t i o n c o e f f i c i e n t s between a l l of the features from Subject #1 were calculated and are summarized i n Table A . l . These r e s u l t s show that the c o r r e l a t i o n c o e f f i c i e n t s between the features that were expected to be p a r t i c u l a r l y r e l a t e d are the strongest, with c o e f f i c i e n t values a l l greater than 0.77. Hence, t h i s demonstrates that there i s a strong consistent r e l a t i o n s h i p between features i n the thumb movement and features i n the sing l e t r i a l o u t l i e r pattern. In p a r t i c u l a r , the r e l a t i o n s h i p between features i n the o u t l i e r pattern and i n the thumb movement was examined using the z-test for the difference between c o r r e l a t i o n s calculated on dependent samples (see Steiger [5A]). The r e s u l t s from these t e s t s are also summarized i n Table A . l . They show that features 3 and A correlated with feature 1 s i g n i f i c a n t l y 86 TABLE 4.1 SINGLE TRIAL FEATURE STATISTICS F e a t u r e C o r r e l a t i o n M a t r i x f e a t u r e 1 f e a t u r e 2 f e a t u r e 3 f e a t u r e 4 f e a t u r e f e a t u r e 1 1 .0 f e a t u r e 2 0.76 1.0 f e a t u r e 3 0.78 0.51 1.0 f e a t u r e 4 0.88' 0.69 0.71 1.0 f e a t u r e 5 0.60 0.80 0.41 0-.57 1.0 z - t e s t on the D i f f e r e n c e Between C o r r e l a t i o n s C a l c u l a t e d on Dependent Samples C o r r e l a t i o n Coef f i c i e n t s z p (one - s i d e d ) f e a t u r e 4 f e a t u r e 5 f e a t u r e 1 0.78 0.51 1.91 0 .029 f e a t u r e 2 0.88 0.69 1 .67 0 .048 f e a t u r e 3 0.60 0.80 1 . 55 0 .060 87 more strongly than with feature 2 (p < 0.05). As expected, the c o r r e l a t i o n between feature 5 and feature 2 was larger than that between feature 5 and feature 1, but t h i s difference achieved only a marginal l e v e l of s i g n i f i c a n c e . 4.8 A p p l i c a t i o n of Dynamic Time Warping to O u t l i e r Patterns A l l the i n i t i a l work with actual EEG was c a r r i e d out on the data from one subject, r e f e r r e d to as Subject #1. These i n i t i a l i n v e s t i g a t i o n s , revealed that the use of dynamic time warping (DTW) provided the best quantitative measure of performance for the s i n g l e t r i a l processing method compared to the other previous a n a l y s i s . Hence, DTW analysis was u l t i m a t e l y applied to a l l four of the subjects included i n t h i s study. S p e c i f i c r e s u l t s are summarized i n the following subsections. 4.8.1 Standard O u t l i e r Patterns using Dynamic Time Warping DTW as described by Roberts et a l . [55] was used to obtain standard (template) representative s i n g l e t r i a l active o u t l i e r patterns for each subject. The time warping procedure attempts to best match waveform A to waveform B by s h i f t i n g , expanding or contracting the time scale of waveform A i n such a manner that minimizes the "cost" of warping. The cost of warping i s based on the following cost function [55] C = Q(A,B,W) - X P(w) 4 * 1 5 where w i s the warped time function (warped time axis) used to warp waveform A, Q i s the c o r r e l a t i o n between warped waveform A and waveform B, P 88 i s a penalty function and X ishe penalty c o e f f i c i e n t . The penalty function i s nonlinear such that the penalty for large expansions or contractions i s proportionately much higher than the penalty for small expansions or contractions. The X c o e f f i c i e n t i s a tuning parameter that d i r e c t l y e f f e c t s how expensive i t i s to warp. In a l l the DTW applications used i n t h i s study a X=75.0 was u t i l i z e d because i t was found by t r i a l and error that t h i s value produced reasonable warped waveforms; smaller values of produced X extreme warpings whereas larger values of X produced warpings that were only s h i f t e d i n time and contained very l i m i t e d expansions or contractions. A standard pattern for each of the four subjects was achieved by using the procedure recommended by Roberts et a l . [55]. The set of active s i n g l e t r i a l patterns was warped against each pattern i n that set. The pattern that produced the lowest mean cost and variance across the set was then selected as the best representative pattern of that set. The standard pattern was then constructed by averaging together a l l the patterns i n the set a f t e r being warped to the above selected pattern. Plots of the standard patterns for a l l four subjects using LSQ, GM and GM2 parameter estimation are provided i n Figure A.15. A.8.2 DTW Cost S t a t i s t i c s on Individual Subjects Once a standard active case pattern was obtained for each subject, i t was warped against the 15 t r i a l s of active o u t l i e r patterns and the 15 t r i a l s of i d l e o u t l i e r patterns for each subject. Each time a pattern i s warped to the standard pattern a cost value i s produced. This cost r e f l e c t s how well the s i n g l e t r i a l pattern " f i t " the standard pattern. The lower the cost the S T A N D A R D O U T L I E R P A T T E R N S GM2 M O D E L P A R A M E T E R S SUBJECT #1 SUBJECT #2 1.0 1.5 2.0 2,5 3.0 3.5 Tn* (Saoonds) 4.0 4.5 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 T I M (Saoonds) 50 • 40 •*» 30 C CS 20 o to • 0 > -10 — -20 • cc -30 -40 1_ SUBJECT #3 0.0 0.5 1.0 _ l _ 1.5 2.0 2.5 TBM (Saoonds) _1_ _ l _ 3.0 3.5 4.0 4.5 SUBJECT #4 Figure 4.15a 1.5 2.0 2.5 TIM (Saoonds) oo 06 S T A N D A R D O U T L I E R P A T T E R N S L S Q M O D E L P A R A M E T E R S SUBJECT #1 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 +.0 4.5 TIM Oaoonda> SO • 40 9. 4» 30 C a 20 o 10 • 0 > 4* -10 o — -20 • OS -30 -40 SUBJECT #2 -L _I_ &6 o5 Co 15 £0 2s 3JO 53 <ui its T I M CSaoonds) _ l _ SUBJECT #3 1.5 2.0 2.5 TBM O c o o n d s ) SUBJECT #4 0.5 1.0 1.5 2.0 2.5 T I M (Seconds) 3.0 3.5 4.0 4.5 Figure 4.15c 92 better the f i t . Therefore, i t would be expected that the costs i n the active cases would be smaller than the costs i n the i d l e case. A t - t e s t designed to te s t the d i f f e r e n c e between two means, given by [56] t = A.16 / + s 2 2 / N df = 2N - 2 = 28 where x1 and x 2 are the sample means, S j 2 and s 2 2 are the sample variances, and N i s the number of active and i d l e cases, was applied to the mean costs for the active and i d l e cases. The r e s u l t s of t h i s t e s t for GM2, GM, and LSQ model parameters for a l l four subjects are summarized i n Table A.2. This te s t shows that i n the GM2 case the difference between the means i s h i g h l y s t a t i s t i c a l l y s i g n i f i c a n t (p < .001). The mean differences i n the GM and LSQ case are also s t a t i s t i c a l l y s i g n i f i c a n t (ranging from p < .001 to p < .01) except for Subject #4 with LSQ parameters where the difference was not s i g n i f i c a n t . These r e s u l t s strongly support the expectation that the costs from the active case would be smaller than the costs from the i d l e case. In turn, t h i s implies that on average the active case patterns f i t the standard patterns much better than the i d l e case patterns. 4.8.3 Grouped DTW Cost S t a t i s t i c s The average active and i d l e cost values from the four subjects were considered together to provide i n f e r e n t i a l s t a t i s t i c s about the actual popu-l a t i o n (for a l l possible s u b j e c t s ^ of mean differences between i d l e and active cases. For each subject, the difference between the average i d l e cost (over the 15 i d l e t r i a l s ) and the average active cost (over the 15 active t r i a l s ) was evaluated. The differences appear i n Table 4.3 and form 93 TABLE 4.2 DIFFERENCE BETWEEN IDLE AND ACTIVE WARPING COSTS t - T e s t R e s u l t s f o r t h e D i f f e r e n c e Between Means: df=28 SUBJECT GM2 MODEL 1 2 3 4 DIFFERENCE BETWEEN THE MEANS 13.3 16.4 1 1.6 11.5 TWO-SIDED t-VALUE 6.05 5.04 4.66 6.25 P < 0.001 0.001 0.001 0.001 GM MODEL 1 2 3 4 8.7 12.7 6.1 4.9 4.18 3.74 2.93 3.35 0.001 0.001 0.01 0.01 LSQ MODEL 1 2 3 4 12.7 17.2 12.9 3.5 4.08 3.65 2.96 0.65 0.001 0.002 0.01 X 9A the basis of the inference to be made for a l l possible subjects. This was c a r r i e d out by considering the n u l l hypothesis that the mean difference of the actual population i s zero and then applying a s t a t i s t i c a l t e s t to deter-mine whether that hypothesis should be rejected. The s t a t i s t i c a l t e s t was a t - t e s t designed to te s t the difference between two means with correlated (paired) samples. I t i s given by [50] MD SEMD df = N - 1 = 3 4.17 where MD i s the sample mean difference (mean of the paired d i f f e r e n c e s ) , S E ^ i s the standard error of t h i s sample mean difference and N equals the number of mean differences (number of subjects). The standard error of the sample mean difference i s given by V N where s^ i s the sample standard deviation of the mean dif f e r e n c e . With three degrees of freedom a 95% confidence i n t e r v a l for the mean difference i s given by [50] C I g 5 = MD ± 3.18(SE M D) A.19 The group s t a t i s t i c s for GM2, GM, and LSQ parameter estimation are summarized i n Table A.3. In the GM2 case the s t a t i s t i c s imply that, even given t h i s small sample of four subjects, the hypothesis that the mean differen c e of the actual population of mean differences i s zero i s strongly rejected (p < .002). A l t e r n a t i v e l y , i n terms of confidence i n t e r v a l s , i t was found that with 95% confidence the i n t e r v a l of 9.7 to 16.8 contains the actual mean difference value. The implications are s i m i l a r but less s i g n i f i -cant i n the GM and LSQ cases. These r e s u l t s imply that the mean cost d i f f e r -ences between active and i d l e t r i a l s based on a l l possible subjects i s highly u n l i k e l y to be zero. 95 TABLE 4.3 GROUP STATISTICS Mean C o s t D i f f e r e n c e Between A c t i v e and I d l e Cases GM2 MODELING GM MODELING LSQ MODELING S u b j e c t # 1 13.3 8.7 12.7 2 16.4 12.7 17.2 3 11.6 6.1 12.9 4 11.5 4.9 3.5 Mean C o s t 13.2 8.1 11.6 S t d . Dev. 2.28 3.45 .5.77 S t d . E r r o r 1.14 1 .73 2.88 Summary f o r GM2 M o d e l i n g Mean D i f f e r e n c e = 13.2 two s i d e d t - t e s t t = 11.58 : H 0 r e j e c t e d p < 0.002 95% C o n f i d e n c e I n t e r v a l : 9.6 t o 16.8 Summary f o r GM M o d e l i n g Mean D i f f e r e n c e = 8.1 two s i d e d t - t e s t t = 4.68 : H Q r e j e c t e d p < 0.05 95% C o n f i d e n c e I n t e r v a l : 2.6 t o 13.6 Summary f o r LSQ M o d e l i n g Mean D i f f e r e n c e = 11.6 two s i d e d t - t e s t . t = 4.03 : Hjj r e j e c t e d p < 0.05 95% C o n f i d e n c e I n t e r v a l : 2.4 t o 20.8 96 TABLE 4.4 BAYESIAN CLASSIFICATION Assuming P r o b a b i l i t i e s of E i t h e r Case O c c u r r i n g are Equal a) Cost of M i s c l a s s i f i c a t i o n set Equal SUBJECT BOUNDARY VALUE ACTIVE CLASSIFIED CORRECTLY FALSE POSITIVE GM2 MODEL 1 2 3 4 GM MODEL 1 2 3 4 LSQ MODEL 1 2 3 4 1 1 .35 19.88 19.10 17.05 15. 19. 16. 13. 28 10 93 26 22.58 41 .06 34.08 20.52 14/15 15/15 13/15 14/15 56/60=93% 12/15 15/15 12/15 13/15 52/60=87% 13/15 14/15 13/15 5/15 45/60=75% 2/15 1/1 5 5/15 3/15 11/60=18% 4/15 4/1 5 8/15 5/15 21/60=35% 4/15 5/15 7/15 2/15 18/60=30% b) Cost of M i s c l a s s i f y i n g an I d l e Case as A c t i v e Set to be 5 Times Greater GM2 MODEL 1 9.18 14/15 0/1 5 2 19.88 12/15 1/1 5 3 19.10 10/15 0/15 4 17.05 12/15 0/15 48/60=80% 1/60=1.7% GM MODEL 1 8.57 9/15 0/15 2 15.58 11/15 2/1 5 3 7.53 2/1 5 0/15 4 6.38 0/15 0/15 22/60=37% 2/60=3.3% LSQ MODEL 1 12.22 4/15 ' 1/15 2 28.68 5/15 2/15 3 21 .35 7/15 3/1 5 4 2.92 0/15 0/1 5 16/60=26% 6/60=10% 97 A.8.4 Bayesian C l a s s i f i c a t i o n of Active Cases versus Idle Cases A Bayesian c l a s s i f i e r was applied to the cost values to c l a s s i f y a ctive cases verses i d l e cases. It was assumed that the cost values had a Gaussian d i s t r i b u t i o n . This c l a s s i f i c a t i o n was c a r r i e d out under two d i f f e r e n t conditions: cost of m i s c l a s s i f i c a t i o n set equal and cost of m i s c l a s s i f y i n g an i d l e case as an active case set to be f i v e times greater. In both conditions the p r o b a b i l i t i e s of e i t h e r case occurring were set equal. The r e s u l t s of the c l a s s i f i c a t i o n under both conditions across the four subjects for GM2, GM and LSQ are provided i n table 4.4. Under the f i r s t condition 93% of the GM2 active cases were c o r r e c t l y c l a s s i f i e d with 18% of the i d l e cases being c l a s s i f i e d as a c t i v e . The number of f a l s e p o s i t i v e s i s f a i r l y high i f the active cases are going to be used i n a control a p p l i c a -t i o n . Hence, the second c l a s s i f y i n g condition was c a r r i e d out were the decision boundary was moved closer to the active case mean by increasing the cost of a f a l s e p o s i t i v e . In t h i s case the percentage of GM2 active cases c o r r e c t l y c l a s s i f i e d was reduced to a s t i l l very respectable 80% but i n so doing reduced the f a l s e p o s i t i v e s to a very low 1.7% (only one i d l e case out of s i x t y was c l a s s i f i e d i n c o r r e c t l y ) . The r e s u l t s using GM and LSQ o u t l i e r patterns under the f i r s t condi-t i o n were f a i r l y good i n c l a s s i f y i n g active cases c o r r e c t l y but both had a very high percentage of f a l s e p o s i t i v e s . Under the second condition, how-ever, the c l a s s i f i c a t i o n performance using GM and LSQ o u t l i e r patterns f e l l o f f dramatically. This i s perhaps the best demonstration of the s u p e r i o r i t y of u t i l i z i n g GM2 parameter estimates i n the si n g l e t r i a l processing method. 98 4.8.5 Cross V a l i d a t i o n and Intersubject R e l i a b i l i t y Using DTW A combined measure of the cross v a l i d a t i o n of a standard pattern on a unrelated set of o u t l i e r patterns and the intersubject r e l i a b i l i t y of the standard patterns was obtained by applying DTW with the standard pattern from one subject to the o u t l i e r patterns from a d i f f e r e n t subject. The GM2 standard patterns from each subject were applied to the GM2 o u t l i e r patterns of each of the other subjects. This resulted i n twelve ad d i t i o n sets of active and i d l e costs. The t - t e s t for the difference between two sample means was applied, i n a s i m i l a r manner as i n Subsection 4.8.2, and the diff e r e n c e was h i g h l y s i g n i f i c a n t (p < 0.001) i n every case. In addition, i n a s i m i l a r manner as i n subsection 4.8.4, Bayesian c l a s s i f i c a t i o n with the cost of m i s c l a s s i f y i n g an i d l e case as an active case set to be f i v e times greater was c a r r i e d out and the r e s u l t s are given i n Table 4.5. On average the percent correct using cross-matched standard patterns f e l l o f f by about 16% when compared to the percent correct using matched standard patterns for subjects 1,3 and 4 while i t stayed exactly the same for subject 2. The percentage of f a l s e p o s i t i v e s was almost exactly same for both the matched and cross-matched cases. The o v e r a l l average of active cases c o r r e c t l y c l a s s i f i e d using standard patterns from one subject on the o u t l i e r patterns from the other subjects was 67% with only a 3% o v e r a l l average of f a l s e p o s i t i v e s . These r e s u l t s strongly i n d i c a t e that the standard patterns do cross v a l i d a t e on data that was not used i n the construction of these p a t t e r n s and t h a t these p a t t e r n s p r o v i d e s u b s t a n t i a l i n t e r s u b j e c t r e l i a b i l i t y . 99 TABLE 4.5 SUMMARY OF BAYESIAN CROSS-MATCHED CLASSIFICATION C o s t of C l a s s i f y i n g an I d l e Case as an A c t i v e Case Set t o be F i v e Times G r e a t e r S t a n d a r d Data from S u b j e c t # P a t t e r n # 1 2 C F C F C F C F 1 14/15 0/15 12/15 1/15 9/15 0/15 8/15 0/15 2 11/15 0/15 12/15 1/15 8/15 0/15 1 1/15 0/15 3 12/15 0/15 12/15 1/15 10/15 0/15 9/15 0/15 4 1 1/15 1/1 5 12/15 1/15 6/1 5 0/15 12/1 5 0/1 5 A c t i v e T r i a l C l a s s i f i e d C o r r e c t l y F = F a l s e P o s i t i v e AVERAGED RESULTS W i t h Matched S t a n d a r d P a t t e r n W i t h C r o s s - M a t c h e d S t a n d a r d P a t t e r n S u b j e c t Number 1 2 3 4 OVERALL AVERAGE A c t i v e C o r r e c t 93% 80% 66% 80% 80% F a l s e P o s i t i v e 0% 7% 0% 0% 2% A c t i v e C o r r e c t 75% 80% 51% 62% 67%" F a l s e P o s i t i v e 2% 7% 0% 0% 3% 100 CHAPTER 5 CONCLUSION 5.1 Summary of Major Results and Related Conclusions AR Parameter Estimation Simulation studies on AR parameter estimation demonstrated that the robust general maximum l i k e l i h o o d (GM) methods performed almost as well as the l e a s t squares (LSQ) method on Gaussian processes and s i g n i f i c a n t l y better on additive o u t l i e r (AO) contaminated Gaussian processes. Amongst the GM methods the GM2 method provided the best performance. However, i t should be noted that the GM2 estimates are also the most computationally expensive to c a l c u l a t e . In terms of the sing l e t r i a l o u t l i e r processing method, the most important fi n d i n g from these simulation studies i s that given the AO model (see equation 3.37) the robust estimation methods, i n p a r t i c u l a r the GM2 method, demonstrated a strong a b i l i t y to model the process x^ without being unduly influenced by the additive o u t l i e r s v^. Neurological Premise and O u t l i e r E x t r a c t i o n The basic neurological premise for the single t r i a l processing method i s that event r e l a t e d p o t e n t i a l s have an additive o u t l i e r e f f e c t on the ongoing EEG process. I f the o u t l i e r content could be extracted from the r e s u l t i n g o v e r a l l combined process, then si n g l e t r i a l event r e l a t e d informa-t i o n could be obtained. Simulation studies demonstrated that a robust s i g n a l estimator, which u t i l i z e s robustly estimated AR model parameters, has a d i s t i n c t a b i l i t y to extract a s i g n i f i c a n t amount of the additive o u t l i e r 101 content from a contaminated process. This extraction process was found to be most e f f e c t i v e when GM2 parameter estimates were u t i l i z e d . This r e s u l t should be expected since the GM2 parameter estimation i s based on an estimated cleaned signa l x^ and hence, i t has the best opportunity to provide a good e s t i m a t e d model r e p r e s e n t i n g the actual process x^. The better the model of x^ the b e t t e r the expected performance of the o u t l i e r e x traction process. Spectral Analysis of EEG AR sp e c t r a l analysis provided a great deal of i n s i g h t i n t o both the EEG s i g n a l i t s e l f and into the a p p l i c a t i o n of AR modeling to the EEG s i g n a l . EEG data used i n t h i s t hesis work was c o l l e c t e d during an ac t i v e task involving motor a c t i v i t y and an i d l e task not involving motor a c t i v i t y . Spectral analysis demonstrated that the si g n a l c h a r a c t e r i s t i c s of these EEG signals were t y p i c a l l y changing at a r e l a t i v e l y r a p id rate. Hence, an EEG segment s i z e as small as p r a c t i c a l parameter estimation considerations would allow was u t i l i z e d . Ultimately, i n the single t r i a l processing method, the EEG epochs were broken down into 1.5 second segments o f f s e t by 0.75 seconds. AR sp e c t r a l analysis was also the key t o o l i n determining appropriate AR model orders. It was found that orders i n the range of 12 to 14 were s u i t -able for i d l e task EEG while orders i n the range of 20 to 24 were su i t a b l e for active task EEG. F i n a l l y , the study of prewhitened AR spe c t r a l analysis was useful i n providing some i n s i g h t into how well the AR model represented the information i n the EEG s i g n a l . This study indicated that, given the s e l e c t i o n of an appropriate model order, the AR model does represent much of the information contained i n a short segment of EEG. 102 Single T r i a l Outlier Processing I n i t i a l investigation into the single t r i a l o u t l i e r processing method was carried out on the EEG data from one subject. I t was shown, through conventional averaging analysis, that the cleaned active task EEG did not contain any s i g n i f i c a n t motor related potentials. This r e s u l t indicated that much of the motor related a c t i v i t y had been extracted from the active task EEG signal by the application of the cleaning process. By using the o u t l i e r information extracted from the active EEG, single t r i a l o u t l i e r patterns were produced. These patterns had strong s i m i l a r i t i e s to previous results using conventional averaging techniques over many t r i a l s of active EEG. By averaging active t r i a l o u t l i e r patterns together i t was demonstrated that much of the information was consistent across active t r i a l s whereas, i n contrast, the average of i d l e case patterns showed that there was no s i g n i f i -cant information that was consistent across i d l e t r i a l s . In addition, t h i s averaging also demonstrated that there was a s i g n i f i c a n t amount of informa-t i o n i n the active task patterns that was unique to individual t r i a l s which was l o s t when the patterns were averaged together. I t i s , therefore, expected that t h i s same loss of information i s occurring i n the conventional averaging method of EEG analysis. I t was shown that consistent features i n the active o u t l i e r patterns were strongly correlated with features i n the thumb movement. This analysis demonstrated that there was a strong relationship between the information i n the o u t l i e r patterns and events i n the thumb movements on a t r i a l by t r i a l basis. I t was found i n the i n i t i a l investigations that dynamic time warping 103 (DTW) analysis provided the best quantitative information on the performance of the sing l e t r i a l processing method. Hence, DTW analysis was applied to a l l four of the subjects used i n t h i s i n i t i a l i n v e s t i g a t i o n into the sing l e t r i a l processing method. Through the a p p l i c a t i o n of DTW, standard repre-sentative active s i n g l e t r i a l o u t l i e r patterns for each subject were obtained. The o u t l i e r patterns from both the active t r i a l s and the i d l e t r i a l s were warped against the standard patterns. With each warping an associated cost value was obtained which r e f l e c t e d how well the o u t l i e r patterns f i t the standard pattern. These cost values revealed that there was a hi g h l y s t a t i s t i c a l l y s i g n i f i c a n t (p < .001) differ e n c e between the i d l e and active mean costs across a l l four subjects with o u t l i e r patterns derived with GM2 parameter estimates. The cost values from the four subjects pooled together i n a group, demonstrated that the mean of the actual population of mean differences between active and i d l e cases was highly (p < .002) u n l i k e l y to be equal zero. Bayesian c l a s s i f i c a t i o n was applied to the warping cost values to c l a s s i f y a ctive patterns versus i d l e patterns. I t was found that with the cost of m i s c l a s s i f y i n g an i d l e case as an active case set to be f i v e times greater, 80% of the GM2 active patterns were c l a s s i f i e d c o r r e c t l y while only 1.7% of the i d l e cases were i n c o r r e c t l y c l a s s i f i e d as a c t i v e . This analysis also demonstrated the s u p e r i o r i t y of u t i l i z i n g GM2 parameter estimates because with the u t i l i z a t i o n of GM and LSQ parameter estimates the c l a s s i f i -c ation performance f e l l o f f dramatically. Bayesian c l a s s i f i c a t i o n was also applied to the cost values obtained from using standard patterns from one subject on the o u t l i e r patterns from the other subjects. These r e s u l t s strongly indicated that the standard 104 patterns do cross v a l i d a t e on data that was not used i n the construction of these patterns and that these patterns provide substantial intersubject r e l i a b i l i t y . In conclusion, the pursuance of robust methods to deal with the ranging Gaussian properties of the EEG si g n a l l e d to the development of a singl e t r i a l processing method based on u t i l i z i n g o u t l i e r information. The v a l i d i t y of t h i s processing method to extract event r e l a t e d information from active task EEG has been established through the r e s u l t s obtained from the investigations undertaken i n t h i s t hesis work. 5.2 Areas for Future Investigation There are many areas involved with the possible further improvement of the s i n g l e t r i a l processing method that should be pursued i n future i n v e s t i -gations. Some of the most important recommended areas are: 1) Pursue models that w i l l better represent the underlying s i g n a l since the processing method i s fundamentally based on these models. This may involve the a p p l i c a t i o n of d i f f e r e n t types of models such as those based on orthonormal functions. Regardless of the type of model employed, the accuracy of the estimated parameters i s an important issue. In terms of using GM estimation on AR models, further work could be c a r r i e d out to improve the parameter estimates. One p a r t i c u l a r aspect to consider i s the u t i l i z a t i o n of methods that are more robust than MEM estimation, such as based methods, to provide the s t a r t i n g estimates i n the GM i t e r a t i o n procedure. 105 2) Study the e f f e c t s of varying the segment s i z e to determine the length that i s best suited for the representation of the underlying EEG process. For instance, i t may be found that the s i g n a l c h a r a c t e r i s t i c s of the under-l y i n g process are changing slowly enough such that the modeling of longer segments would allow for an improved representation of the underlying process. 3) Further investigate the o u t l i e r detection process to determine ways i n which the performance could be improved. One s p e c i f i c area to consider, as suggested by Martin [45], i s to u t i l i z e a cleaning process that makes use of both forward and backward p r e d i c t i o n i n the estimation of the s i g n a l x^. 4) Peruse a l t e r n a t i v e s to the current method of processing the o u t l i e r information. The current method i s r e l a t i v e l y unsophisticated i n that i t i s simply a c a r e f u l smoothing of the extracted o u t l i e r information. Investiga-t i o n into other approaches of processing t h i s information may prove to be b e n e f i c i a l i n revealing a d d i t i o n a l information that may be contained i n the o u t l i e r data. Future empirical EEG experimentation should be c a r r i e d out on the si n g l e t r i a l processing method. The i n i t i a l goals of these i n v e s t i g a t i o n s should be to further v a l i d a t e the method on a new set of EEG data. A recom-mended paradigm would be to c o l l e c t a t r a i n i n g set of a c t i v e t r i a l s for the construction of a standard pattern. * Then c o l l e c t a set of intermixed active and i d l e t r i a l s , perhaps at the d i s c r e t i o n of the subject, on which the c l a s s i f i c a t i o n performance of the processing method could be further evaluated. Later goals of these investigations should be oriented to using the processing method to learn more about motor p o t e n t i a l s , p a r t i c u l a r l y the 106 motor p o t e n t i a l s rom disabled persons, so that when u t i l i z e d i n a control a p p l i c a t i o n these p o t e n t i a l s can be taken advantage of i n the most appro-p r i a t e manner. F i n a l l y , i n vestigations on making the sing l e t r i a l processing method work i n r e a l time must be undertaken. As i t stands, the method i s very computationally i n t e n s i v e . Although some improvements could undoubtedly be made i n the e f f i c i e n c y of these computations, the most s i g n i f i c a n t advances towards t h i s goal would l i k e l y be i n the implementation of some or a l l of the component processes i n s p e c i a l i z e d hardware. 5.3 S i g n i f i c a n t Contributions The s i g n a l processing method developed i n t h i s t hesis work i s a s i g n i f i c a n t contribution. Modeling the underlying s i g n a l and then extracting the o u t l i e r content from the underlying s i g n a l i s a unique approach to deriving very low l e v e l and r e l a t i v e l y short event type information from an ongoing process. The path that l e d to the development of t h i s processing method, also led to the understanding that ranging l e v e l s of Gaussianity i n the EEG si g n a l requires that serious consideration be given to the use of robust methods i n the future a p p l i c a t i o n of various types of EEG si g n a l processing. Also, a successful approach to the s e l e c t i o n of appropriate AR model orders, the s e l e c t i o n of appropriate EEG segment lengths and the assessment of the r e l a t i v e a b i l i t y of AR models to represent the EEG si g n a l were established through the studies on AR sp e c t r a l a n a l y s i s . The a b i l i t y to c o n s i s t e n t l y acquire event r e l a t e d information from a singl e t r i a l i s an important contribution to the f i e l d of EEG si g n a l 107 a n a l y s i s . In addition, the neurological premise inv o l v i n g the way i n which event r e l a t e d information i s contained i n the o v e r a l l EEG s i g n a l i s estab-l i s h e d as a v i a b l e model. This model should be considered when attempting to understand event r e l a t e d p o t e n t i a l s and t h e i r r e l a t i o n s h i p to ongoing EEG processes. F i n a l l y , the work i n t h i s t h e s i s , taken as a whole, represents an important contribution towards the ultimate goal of harnessing EEG signals for control a p p l i c a t i o n s . It overcomes perhaps one of the greatest obstacles by providing the framework for the extraction of u s e f u l information from sing l e t r i a l EEG. 108 REFERENCES 1. Westerkamp,J.J. and Aunon.J.I., "Optimum Multielectrode A P o s t e r i o r i estimates of Single-Response Evoked Po t e n t i a l s " , IEEE Trans. Biomedical Engineering, Vol. BME-34, No. 1, Jan., 1987, pp 13-22. 2. Heinze.H.J., Kunkel.H., "ARMA - F i l t e r i n g of Evoked P o t e n t i a l s " , Methods of Information i n Medicine, 1984, Vol. 23, pp 29-36. 3. V i d a l . J . J . , "Real-Time Detection of Brain Events i n EEG", Proceedings of the IEEE, V o l . 65, No. 5, 1977, pp 633-641. 4. Kornhuber,H.H. und Deecke,L., "Hirnpotentialanderungen beim Menschen vor und nach Willkurbewegungen, D a r g e s t e l l t mit Magnetbandspeicherung und Ruckwartsanalyse", Pflugers Arch. ges. Physiol., 1964, 281:52. 5. Landau,B., E s s e n t i a l Human Anatomy and Physiology, Scott Foresman Company, 1976, pp 227. 6. Brunia.C.H.M. and Van den Bosch,W.E.J., "Movement-related Slow P o t e n t i a l s . I. A Contrast Between Finger and Foot Movements i n Right-handed Subjects", Electroencephalography and C l i n i c a l Neurophysiology, 1984, 57: 515-527. 7. Grunewald.G. and Grunewald-Zuberbier,E. , "Cerebral Potentials During Voluntary Ramp Movements i n Aiming Tasks", T u t o r i a l s i n ERP Research: Endogenous Components, A.W.K. G a i l l a r d and W. R i t t e r (eds.), North-Holland Publishing Company, 1983. 8. Vaughan, J r . , H.G., Costa, L.D. and R i t t e r , W., "Topography of the Human Motor P o t e n t i a l " , Electroencephalography and C l i n i c a l Neurphysiology, v o l . 25, pp. 1-10, 1968. 9. Vaughan, Jr.,H.G., Bossom.J. and Gross,E.G., " C o r t i c a l Motor Potentials i n Monkeys Before and Af t e r Upper Limb Deafferentation", Exp. Neurol., 1970, 26: 253-262. 10. Papakostopoulos,D., "A No-stimulas, No-response, Event-related P o t e n t i a l of the Human Cortex", E l e c t r o e n c e p h a l o g r a p h y and C l i n i c a l Neurophysiology, 1980, 48: 622-638. 11. Persson.J., "Comments on Estimations and Tests of EEG Amplitude D i s t r i b u t i o n s " , Electroencephalography and C l i n i c a l Neurophysiology, 1974, 37: 309-313. 12. McEwen,J.A. and Anderson,G.B., "Modeling the S t a t i o n a r i t y and Gaussianity of Spontaneous E l e c t r o e n c e p h a l o g r a p h i c A c t i v i t y " , IEEE Trans. Biomedical Engineering, Vol. BME-22 No.5, 1975, pp 361-369. 109 13. McGillem.C.D., Aunon.J.I. and Childers,D.G., "Signal Processing i n Evoked Poten t i a l Research: Applications of F i l t e r i n g and Pattern Recognition", CRC C r i t i c a l Reviews i n Bioengineering, 1981, pp 225-265. 14. Persson,J., "Comments on 'Modeling the S t a t i o n a r i t y and Gaussianity of Spontaneous Electroencephalographic A c t i v i t y ' " , IEEE Trans. Biomedical Engineering, Vol. BME-24, 1977, p 302. 15. Weiss,M.C, "A New Test for EEG Gaussian Amplitude D i s t r i b u t i o n " , Proc. IEEE Frontiers of Engineering j.n Health Care, 1979, pp 309-313. 16. E l u l , R . , i n Progress i n B i o m e d i c a l Research, Fogeland L . J . and George,F.W., eds., Vol. 1, p 131, Spartan Books, 1967. 17. Elul.R., i n Data A c q u i s i t i o n and Processing i n Biology and Medicine, Enslein.K., ed., Vol. 5, p 93, Pergaman Press, 1968. 18. Elul.R., "Gaussian Behavior of the Electroencephalogram: Changes during Performance of Mental Task", Science, Vol. 164, 1969, pp 328-331. 19. Elul.R., Internat. Rev. Neurobiol., 15, p 227, 1972. 20. Seigel.A., "Stochastic Aspects of the Generation of The Electroencephal-gram", J. Theor. Biol-92, 1981, pp 317-339. 21. Bernstein, S., Math. Ann., v o l . 97, p. 1, 1927. 22. Anninos.P., Zenone.S. and Elul.R., " A r t i f i c a l Neural Nets: Dependence of the EEG Amplitude's P r o b a b i l i t y D i s t r i b u t i o n on S t a t i s t i c a l Para-meters", J. Theor. Biol-103, 1983, pp 339-348. 23. Kay,S.M. and Marple,S.L., "Spectrum Analysis - A Modern Perspective", Proceedings of the IEEE, Vol.69, No.11, 1981. 24. Beamish,N. and Priestly,M.B., "A study of autoregressive and window sp e c t r a l estimation," Appl. S t a t i s t . , vol.30, n o . l , pp.41-58, 1981. 25. Jansen.B.H., Bourne,J.R. and Ward.J.W., "Autoregressive Estimation of Short Segment Spectra for Computerized EEG Analysis", IEEE Trans. Biomedical Engineering, Vol. BME-28 No. 9, 1981, pp 630-641. 26. Box.G.E.P. and Jenkins,G.M., Time Series Analysis: Forcasting and Control, Holden-Day, 1970. 27. Kavech.M. and Cooper,G.R., "An Empirical Investigation of the Properties of the Autoregressive Spectral Estimator", IEEE Trans. Inform. Theory, Vol. IT-22, 1976, pp 313-323. 28. Andrews,D.F., "The Robustness of Residual Displays", i n Robustness i n S t a t i s t i c s , Launer.R.L. and Wilkinson,G.N. eds., Academic Press, 1979. 110 29. Chambers,R.L. and Heathcote,C,R., "On the Estimation of Slope and the I d e n t i f i c a t i o n of O u t l i e r s i n Linear Regression", Biometrika-68, 1, 1981, pp 21-33. 30. Brown,R.G., Introduction to Random Signal Analysis and Kalman F i l t e r i n g , John Wiley & Sons, 1983. 31. Goodwin,G.C. and Sin.K.S., Adaptive F i l t e r i n g , P r e d i c t i o n , and Control, P r e n t i c e - H a l l , 1984, pp.224-288. 32. Jansen.B.H., "Analysis of biomedical signals by means of l i n e a r modeling," CRC C r i t i c a l Reviews i n Biomedical Engineering, v o l . 12(4), pp.343-392, 1985. 33. Isaksson,A. and Wennberg.A., "Spectral properties of nonstationary EEG signals evaluated by means of Kalman f i l t e r i n g . A p p l i c a t i o n to a v i g i l a n c e t e s t , " i n Quantitative A n a l y t i c Studies i n Epilepsy, P. Kellaway and I. Peterson, Eds., Raven Press, New York, 1976. 34. Ulrych, T.J. and Bishop, T.N., "Maximum Entropy Spectral Analysis and Autoregressive Decomposition", Rev. Geophysics Space Phys., v o l . 13, pp. 183-200, 1975. 35. Burg,J.P., "Maximum Entropy Spectral Analysis", Proceedings of 37th Meeting of the Society of Exploration Geophysicists, Oct. 31, 1967. 36. Anderson,N., "On the C a l c u l a t i o n of F i l t e r C o e f f i c i e n t s for Maximum Entropy Spectral Analysis", Geophysics, Vol. 39, No.l, Feb., 1974, pp. 69-72. 37. Hurvich.C.M., "Data-Dependent Spectral Windows: Generalizing the C l a s s i c a l Framework to Include Maximum Entropy E s t i m a t e s " , Technometrics, Vol. 28, No.3, Aug., 1986, pp 259-268. 38. van den Bos,A., "Alternative i n t e r p r e t a t i o n of maximum entropy s p e c t r a l a n a l y s i s , " IEEE Trans. Inform. Theory, vol.IT-17, pp.493-494, 1971. 39. Box.G.E.P. and Muller.M.E., "A Note on Generation of Normal Deviates", AMS, Vol. 28, 1958, pp 610-611": 40. Lewis,T.G. and Payne,W.H., "Generalized Feedback S h i f t Registor Pseudo-random Number Algorithm", Journal of the Association for Computing Machinery, Vol. 20, No.3, 1973, pp 456-468. 41. Smith,W.D. and Lager,D.L., "Evaluation of simple algorithms for s p e c t r a l parameter analysis of the electroencephalogram," IEEE Trans. Biomed. Eng., vol.BME-33, no.3, pp.352-358, 1986. 42. Kleiner,B., Martin,R.D. and Thomson,D.J., "Robust Estimation of Power Spectra", J. Royal S t a t i s t . S o c , Ser.B, Vol. 41 No. 3, 1979, pp 313-351. I l l 43. Hogg, R.V., "An Introduction to Robust Estimation", i n Robustness i n S t a t i s t i c s , Launer.R.L. and Wilkinson,G.N. eds., Academic Press, 1979. 44. Huber.P.J, "Robust Estimation of Location Parameters", Ann. Math. S t a t i s t . , 35, 1964, pp 73-101. 45. Martin,R.D. and Thomson,D.J., "Robust-Resistant Spectrum Estimation", Proceedings of the IEEE, Vol. 70, No.9, 1982, pp 1097-1115. 46. Ziemer.R.E., and Tranter,W.H., P r i n c i p l e s of Communications; Systems, Modulation, and Noise, Houghton M i f f l i n Company, Boston, 1976, p.57 and p.224. 47. Martin,R.D., "Robust Estimation for Time Series Autoregression", i n Robustness i n S t a t i s t i c s , Launer.R.L. and Wilkinson,G.N. eds., Academic Press, 1979. 48. Martin,R.D., "Robust Estimation of Autoregressive Models", i n Institue of Mathematical S t a t i s t i c s Directions i n Time Series, Brillinger,D.R. and Tiao.G.C, 1978. 49. Ja s p e r , H.H., The Ten-Twenty E l e c t r o d e Placement System of the I n t e r n a t i o n a l F e d e r a t i o n " , E l e c t r o e n c e p h a l o g r a p h y and C l i n i c a l Neurophysiology, v o l . 10, pp. 371-375, 1958. 50. Isaksson.A., Wennberg.A. and Zetterberg,L.."Computer analysis of EEG signals with parametric models," Proc. IEEE, vol.69, no.4, pp.451-461, 1981. 51. Martin,R.D. and Zeh.J.E., "Determining the Character of Time Series O u t l i e r s " , Proceedings of the Amer. S t a t i s t . Assoc., Business and Economics Section, 1977, pp. 818-823. 52. P a p o u l i s . A . , "Minimum-bias windows f o r h i g h - r e s o l u t i o n s p e c t r a l estimates', IEEE Trans. Inform. Theory, vol.IT-19, n o . l , pp.9-12, 1973. 53. Birch.G.E., Lawrence,P.D., Lind.J.C. and Hare.R.D., "Application of Prewhitening to AR Spectral Estimation of EEG", IEEE Trans. Biomedical Engineering, In press, 1988. 54. Steiger, J.H., "Tests for Comparing Elements of a C o r r e l a t i o n Matrix", Psychological B u l l e t i n , v o l . 87, no. 2, pp. 245-251, 1980. 55. Roberts,K., Lawrence.P., Eisen.A. and Hoirch.M., "Enhancement and Dynamic Time Warping of Somatosensory Evoked P o t e n t i a l Components Applied to Patients with Multiple S c l e r o s i s " , IEEE Trans. Biomedical Engineering, Vol. BME-34, No. 6, June 1987, pp. 397-405. 56. Glass,G.V. and Hopkins,K.D., S t a t i s t i c a l Methods i n Education and Psychology, Second E d i t i o n , P r e n t i c e - H a l l , 1984. 

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

Usage Statistics

Share

Embed

Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                        
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            src="{[{embed.src}]}"
                            data-item="{[{embed.item}]}"
                            data-collection="{[{embed.collection}]}"
                            data-metadata="{[{embed.showMetadata}]}"
                            data-width="{[{embed.width}]}"
                            async >
                            </script>
                            </div>
                        
                    
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.831.1-0097955/manifest

Comment

Related Items