- Email: [email protected]
Isolating low frequency activity in EEG spectrum analysis
Electroencephalography and Clinical Neurophysiology, 1979, 4 6 : 2 2 4 - - 2 2 6
224
© Elsevier/North-Holland Scientific Publishers, Ltd.
Technical contribution ISOLATING LOW F R E Q U E N C Y ACTIVITY IN EEG SPECTRUM ANALYSIS RICHARD COPPOLA
Biological Psychiatry Branch, National Institute o f Mental Health, Bethesda, Md. 20014 (U.S.A.) ( A c c e p t e d for publication: August 23, 1978)
Digital s p e c t r u m analysis o f EEG began with the indirect m e t h o d (Walter 1963) and currently is usually p e r f o r m e d directly using the Fast Fourier Transform ( F F T ) algorithm ( D u m e r m u t h and Keller 1973). Of concern to m a n y investigators are the large low f r e q u e n c y values often seen in EEG spectra. In a discussion of EEG spectra (Kellaway and Petersen 1973, p. 200) several investigators n o t e d that m a n y published s p e c t r u m patterns show activity below 2 c/sec that is n o t visible in the raw EEG tracing. T h e y offered the possibility that it was due to baseline shifting and one investigator's c o m m e n t was ' t o let the specialists hassle.' The concern over low frequencies is especially a p r o b l e m when dealing with digital equally spaced data since all frequencies d o w n to zero may be present (Blackman and T u k e y 1958). As a result e x t r e m e l y slow drift if present may show up as subharmonic leakage into the low end s p e c t r u m estimates. Whether such activity is due to amplifier drift or to s o m e physiological process is another question. This technical note presents a filter t e c h n i q u e for isolating low frequency c o m p o n e n t s . Hopefully, with its use the proper disposition of such low f r e q u e n c y activity m a y be determined.
and phase characteristic (Otnes and E n o c h s o n 1972): ¢(f) =
a r c t a n [ --a sin(2~fAt) \1Z a ~ t i !
t (3)
If the phase shift introduced by this filter is of concern it may be eliminated by applying the filter twice, first in the normal time direction and second with the t i m e index reversed. In this m a n n e r the phase shift from the first operation would be cancelled by the phase shift of the second. Orr and H o f f m a n (1974) p o i n t out that this was first suggested by Tukey. A n o t h e r advantage of the double application is that an even steeper f r e q u e n c y characteristic results. Fig. 1 shows the transfer p o w e r gain resulting from the double application for various values of the filter parameter and a sampling time of 10 msec. N o t e that this transfer characteristic is that of the trend that is being removed. Table I gives a F O R T R A N routine that realizes this filter.
Data collection and analysis To d e m o n s t r a t e the effects of this filter digitized 1.0-
Methods
Autoregressive trend removal A simple autoregressive filter may be applied to a t i m e series to obtain the low f r e q u e n c y portion. The s p e c t r u m analysis m a y t h e n be p e r f o r m e d on the residuals after subtracting o u t this low frequency activity. A filter of the f o r m Yi = aYi--1 + (1 - - a)x i
Z .5
9
,8
(1)
is a low pass filter with the f r e q u e n c y response :
1
2
3
4 HZ
iH(f)l 2 =
(1 - - a ) 2 1 -- 2a cos(2~'fAt) + a 2
(2)
Fig. 1. Transfer gain of the double time direction autoregressive filter shown for various values o f the filter parameter and for 100 Hz sampling rate.
I S O L A T I N G LOW F R E Q U E N C Y A C T I V I T Y IN E E G
225
TABLE I
A
B
D o u b l e t i m e d i r e c t i o n autoregressive filter r o u t i n e . C C C C C C C
A is t h e filter p a r a m e t e r A r r a y X c o n t a i n s t h e data t o be filtered The data w i t h t h e t r e n d r e m o v e d is left in X A r r a y Y is t e m p o r a r y storage N is t h e n u m b e r of data p o i n t s
100
AA
= 1 .--A
Y(1) DO100I
= X(1) = 2,N
Y(I) YLAG X(N) DO 150 YLAG X(I)
150
= = = I = = =
2~
one sec
Fig. 2. A: t o p l i n e is t h e s p e c t r u m o f t h e raw E E G , b o t t o m lines are t h e s p e c t r a of t h e E E G w i t h t h e t r e n d s r e m o v e d . B: t o p l i n e is raw E E G , b o t t o m lines are the autoreKressive t r e n d s r e m o v e d b y various values o f t h e filter p a r a m e t e r .
raw E E G was collected a n d several analyses were p e r f o r m e d using a PDP-11 l a b o r a t o r y c o m p u t e r s y s t e m ( C o p p o l a 1977). V e r t e x to r i g h t ear recordings were t a k e n f r o m a n o r m a l v o l u n t e e r using a n E E G a m p l i f i e r w i t h a 3 dB b a n d p a s s o f 0 . 5 - - 100 Hz f o l l o w e d b y a n active low pass filter flat to 4 0 Hz a n d d o w n 42 dB at 60 Hz. Several 2 0 sec or longer e p o c h s were digitized to 10 bits w i t h a 10 msec s a m p l e spacing. S p e c t r u m e s t i m a t e s were c o m p u t e d in t h e following m a n n e r . A n artifact-free E E G r e c o r d o f 1 0 2 4 samples was selected a n d t h e autoregressive t r e n d was r e m o v e d as above. T w o d i f f e r e n t m e t h o d s o f s m o o t h ing t h e raw s p e c t r u m e s t i m a t e s were used in o r d e r t o increase reliability. F o r t h e first m e t h o d , t h a t o f t i m e averaging, t h e 1 0 2 4 samples were divided i n t o eight 128 p o i n t c o n t i g u o u s segments. F o r each s e g m e n t a cosine t a p e r was used as t h e t i m e w i n d o w . T h e t a p e r had the following form
W(n) = 1
'°/U ,6 c/~c
A*Y (I-- 1) + AA~X(1) Y(N) X(N) -- YLAG N-- 1,1,--1 A~YLAG + AA~Y(1) X(I)YLAG
= I ( I -- c o s ( 2 ? r n / P N ) )
J
;0 ~ n ~ (PN/2) ; ( P N / 2 ) < n < N -- ( P N / 2 )
= 1(1 -- cos(2n(N-- n)/PN)) ; N - - (PN/2) ~ n ~ N (4) w h e r e P is t h e p r o p o r t i o n t o t a p e r ( h a l f of t h e port i o n is applied t o e i t h e r e n d o f t h e d a t a ) a n d N - - 1 is t h e n u m b e r of samples in t h e record. T h i s w i n d o w gives good s i d e l o b e s u p p r e s s i o n w i t h o u t a p p r e c i a b l e w i d e n i n g of t h e m a i n l o b e a n d since it is c o n t r o l l e d b y t h e p a r a m e t e r P it is easy to s t u d y t h e effects of
m o r e s t r i n g e n t tapering. T h e r e s u l t i n g 128 data p o i n t s were t h e n p r o c e s s e d b y a s t a n d a r d F F T t o give 64 complex Fourier coefficients from which the modified p e r i o d o g r a m was c o m p u t e d . T h e resulting 8 spectra were t h e n averaged b y a d d i n g t h e c o r r e s p o n d ing 64 s p e c t r u m e s t i m a t e s t o g e t h e r a n d dividing by 8. F o r t h e s e c o n d m e t h o d , t h a t of f r e q u e n c y s m o o t h ing, t h e cosine t a p e r was applied t o t h e e n t i r e filtered data r e c o r d a n d a 1 0 2 4 p o i n t F F T p e r f o r m e d . T h e p e r i o d o g r a m was again c o m p u t e d giving s p e c t r u m e s t i m a t e s at 5 1 2 f r e q u e n c y steps. These e s t i m a t e s were c o m p r e s s e d b y s u m m i n g all a d j a c e n t 8 estim a t e s t h u s resulting in a 64 step s p e c t r u m comparable to t h e t i m e averaged m e t h o d above. Results Fig. 2 s h o w s for various values of t h e filter p a r a m eter t h e autoregressive t r e n d t h a t is r e m o v e d f r o m t h e E E G a n d t h e c o r r e s p o n d i n g s p e c t r a for t h e residual EEG. T h e figure d e m o n s t r a t e s t h a t t h e r e is considerable low f r e q u e n c y c o n t e n t in t h e E E G t h a t is n o t very a p p a r e n t in t h e raw record. T h e s p e c t r a s h o w n are t h e t i m e averaged s p e c t r a c o m p u t e d w i t h a t a p e r of 20%. T h e r e were slight d i f f e r e n c e s b e t w e e n t h e t i m e averaged s p e c t r a a n d t h e f r e q u e n c y s m o o t h e d spectra. T h e m o s t n o t i c e a b l e is for t h e case w h e r e n o filter has b e e n applied. T h e side lobe s u p p r e s s i o n gained b y f r e q u e n c y s m o o t h i n g limits t h e s u b h a r m o n i c cont a m i n a t i o n to t h e first h a r m o n i c .
226
R. COPPOLA
Discussion For the example shown, the largest low frequency contribution is at 1 c/sec and below. This large activity causes some elevation of the spectrum estimates at 2--4 c/sec. This elevation is most probably due to leakage. If a filter parameter of 0.95 is used only the activity below 1 c/sec is removed leaving the 2--4 c/sec range unchanged (see Fig. 1). Thus if accurate spectrum estimates in the 2--4 c/sec region were sought, it would seem desirable to filter the lower activity in this case. Whether the activity below 1 c/sec is physiological or electronic is uncertain; however, it would seem more likely to attribute it to the amplifier. Application of this filter technique will give investigators control over the low frequency content of their EEG data and will hopefully lead to a resolution of this issue.
Summary A method for isolating and removing very slow activity from EEG records prior to spectrum computation is presented. This is accomplished by an autoregressive filter whose frequency transfer characteristic can be easily controlled. Examples of spectra computed for various filter parameter values are shown.
R~sum~
Isolement de l'activitd d l'analyse spectrale de I'EEG
basse
frgquence
dans
L'auteur d6crit une m~thode permettant d'isoler et de supprimer des ondes tr~s lentes sur les enregistre-
ments EEG avant analyse spectrale. Ceci est r~alis6 au moyen d'un filtre autor~gressif dont la caract6ristique de transfert de fr~quence peut 6tre facilement control6e. Des exemples de spectres corrig6s pour diverses valeurs de filtre sont donn6s en illustration.
References Btackman, R.B. and Tukey, J.W. The Measurement of Power Spectra. Dover Publications, New York, 1958, p. 48. Coppola, R. A table driven system for stimulus response experiments. Proc. Digital Equipment User's Soc., 1977, 3: 1219--1222. Dumermuth, G. and Keller, E. EEG spectral analysis by means of Fast Fourier Transform. In: P. Kellaway and I. Petersen (Eds.), Automation of Clinical Electroencephalography. Raven Press, New York, 1973: 145--160. Kellaway, P. and Petersen, I. (Eds.) Automation of Clinical Electroencephalography. Raven Press, New York, 1973. Orr, W.C. and Hoffman, H.J. A 90-min cardiac biorhythm: methodology and data analysis using modified periodograms a n d complex demodulation. IEEE Trans. biomed. Engng, 1974, BME-21 : 130--143. Otnes, R.K. and Enochson, I. Digital Time Series Analysis. Wiley, New York, 1972. Walter, D.O. Spectral analysis for electroencephalograms: mathematical determination of neurophysiological relationships from records of limited duration. Exp. Neurol., 1963, 8: 155--181.
224
© Elsevier/North-Holland Scientific Publishers, Ltd.
Technical contribution ISOLATING LOW F R E Q U E N C Y ACTIVITY IN EEG SPECTRUM ANALYSIS RICHARD COPPOLA
Biological Psychiatry Branch, National Institute o f Mental Health, Bethesda, Md. 20014 (U.S.A.) ( A c c e p t e d for publication: August 23, 1978)
Digital s p e c t r u m analysis o f EEG began with the indirect m e t h o d (Walter 1963) and currently is usually p e r f o r m e d directly using the Fast Fourier Transform ( F F T ) algorithm ( D u m e r m u t h and Keller 1973). Of concern to m a n y investigators are the large low f r e q u e n c y values often seen in EEG spectra. In a discussion of EEG spectra (Kellaway and Petersen 1973, p. 200) several investigators n o t e d that m a n y published s p e c t r u m patterns show activity below 2 c/sec that is n o t visible in the raw EEG tracing. T h e y offered the possibility that it was due to baseline shifting and one investigator's c o m m e n t was ' t o let the specialists hassle.' The concern over low frequencies is especially a p r o b l e m when dealing with digital equally spaced data since all frequencies d o w n to zero may be present (Blackman and T u k e y 1958). As a result e x t r e m e l y slow drift if present may show up as subharmonic leakage into the low end s p e c t r u m estimates. Whether such activity is due to amplifier drift or to s o m e physiological process is another question. This technical note presents a filter t e c h n i q u e for isolating low frequency c o m p o n e n t s . Hopefully, with its use the proper disposition of such low f r e q u e n c y activity m a y be determined.
and phase characteristic (Otnes and E n o c h s o n 1972): ¢(f) =
a r c t a n [ --a sin(2~fAt) \1Z a ~ t i !
t (3)
If the phase shift introduced by this filter is of concern it may be eliminated by applying the filter twice, first in the normal time direction and second with the t i m e index reversed. In this m a n n e r the phase shift from the first operation would be cancelled by the phase shift of the second. Orr and H o f f m a n (1974) p o i n t out that this was first suggested by Tukey. A n o t h e r advantage of the double application is that an even steeper f r e q u e n c y characteristic results. Fig. 1 shows the transfer p o w e r gain resulting from the double application for various values of the filter parameter and a sampling time of 10 msec. N o t e that this transfer characteristic is that of the trend that is being removed. Table I gives a F O R T R A N routine that realizes this filter.
Data collection and analysis To d e m o n s t r a t e the effects of this filter digitized 1.0-
Methods
Autoregressive trend removal A simple autoregressive filter may be applied to a t i m e series to obtain the low f r e q u e n c y portion. The s p e c t r u m analysis m a y t h e n be p e r f o r m e d on the residuals after subtracting o u t this low frequency activity. A filter of the f o r m Yi = aYi--1 + (1 - - a)x i
Z .5
9
,8
(1)
is a low pass filter with the f r e q u e n c y response :
1
2
3
4 HZ
iH(f)l 2 =
(1 - - a ) 2 1 -- 2a cos(2~'fAt) + a 2
(2)
Fig. 1. Transfer gain of the double time direction autoregressive filter shown for various values o f the filter parameter and for 100 Hz sampling rate.
I S O L A T I N G LOW F R E Q U E N C Y A C T I V I T Y IN E E G
225
TABLE I
A
B
D o u b l e t i m e d i r e c t i o n autoregressive filter r o u t i n e . C C C C C C C
A is t h e filter p a r a m e t e r A r r a y X c o n t a i n s t h e data t o be filtered The data w i t h t h e t r e n d r e m o v e d is left in X A r r a y Y is t e m p o r a r y storage N is t h e n u m b e r of data p o i n t s
100
AA
= 1 .--A
Y(1) DO100I
= X(1) = 2,N
Y(I) YLAG X(N) DO 150 YLAG X(I)
150
= = = I = = =
2~
one sec
Fig. 2. A: t o p l i n e is t h e s p e c t r u m o f t h e raw E E G , b o t t o m lines are t h e s p e c t r a of t h e E E G w i t h t h e t r e n d s r e m o v e d . B: t o p l i n e is raw E E G , b o t t o m lines are the autoreKressive t r e n d s r e m o v e d b y various values o f t h e filter p a r a m e t e r .
raw E E G was collected a n d several analyses were p e r f o r m e d using a PDP-11 l a b o r a t o r y c o m p u t e r s y s t e m ( C o p p o l a 1977). V e r t e x to r i g h t ear recordings were t a k e n f r o m a n o r m a l v o l u n t e e r using a n E E G a m p l i f i e r w i t h a 3 dB b a n d p a s s o f 0 . 5 - - 100 Hz f o l l o w e d b y a n active low pass filter flat to 4 0 Hz a n d d o w n 42 dB at 60 Hz. Several 2 0 sec or longer e p o c h s were digitized to 10 bits w i t h a 10 msec s a m p l e spacing. S p e c t r u m e s t i m a t e s were c o m p u t e d in t h e following m a n n e r . A n artifact-free E E G r e c o r d o f 1 0 2 4 samples was selected a n d t h e autoregressive t r e n d was r e m o v e d as above. T w o d i f f e r e n t m e t h o d s o f s m o o t h ing t h e raw s p e c t r u m e s t i m a t e s were used in o r d e r t o increase reliability. F o r t h e first m e t h o d , t h a t o f t i m e averaging, t h e 1 0 2 4 samples were divided i n t o eight 128 p o i n t c o n t i g u o u s segments. F o r each s e g m e n t a cosine t a p e r was used as t h e t i m e w i n d o w . T h e t a p e r had the following form
W(n) = 1
'°/U ,6 c/~c
A*Y (I-- 1) + AA~X(1) Y(N) X(N) -- YLAG N-- 1,1,--1 A~YLAG + AA~Y(1) X(I)YLAG
= I ( I -- c o s ( 2 ? r n / P N ) )
J
;0 ~ n ~ (PN/2) ; ( P N / 2 ) < n < N -- ( P N / 2 )
= 1(1 -- cos(2n(N-- n)/PN)) ; N - - (PN/2) ~ n ~ N (4) w h e r e P is t h e p r o p o r t i o n t o t a p e r ( h a l f of t h e port i o n is applied t o e i t h e r e n d o f t h e d a t a ) a n d N - - 1 is t h e n u m b e r of samples in t h e record. T h i s w i n d o w gives good s i d e l o b e s u p p r e s s i o n w i t h o u t a p p r e c i a b l e w i d e n i n g of t h e m a i n l o b e a n d since it is c o n t r o l l e d b y t h e p a r a m e t e r P it is easy to s t u d y t h e effects of
m o r e s t r i n g e n t tapering. T h e r e s u l t i n g 128 data p o i n t s were t h e n p r o c e s s e d b y a s t a n d a r d F F T t o give 64 complex Fourier coefficients from which the modified p e r i o d o g r a m was c o m p u t e d . T h e resulting 8 spectra were t h e n averaged b y a d d i n g t h e c o r r e s p o n d ing 64 s p e c t r u m e s t i m a t e s t o g e t h e r a n d dividing by 8. F o r t h e s e c o n d m e t h o d , t h a t of f r e q u e n c y s m o o t h ing, t h e cosine t a p e r was applied t o t h e e n t i r e filtered data r e c o r d a n d a 1 0 2 4 p o i n t F F T p e r f o r m e d . T h e p e r i o d o g r a m was again c o m p u t e d giving s p e c t r u m e s t i m a t e s at 5 1 2 f r e q u e n c y steps. These e s t i m a t e s were c o m p r e s s e d b y s u m m i n g all a d j a c e n t 8 estim a t e s t h u s resulting in a 64 step s p e c t r u m comparable to t h e t i m e averaged m e t h o d above. Results Fig. 2 s h o w s for various values of t h e filter p a r a m eter t h e autoregressive t r e n d t h a t is r e m o v e d f r o m t h e E E G a n d t h e c o r r e s p o n d i n g s p e c t r a for t h e residual EEG. T h e figure d e m o n s t r a t e s t h a t t h e r e is considerable low f r e q u e n c y c o n t e n t in t h e E E G t h a t is n o t very a p p a r e n t in t h e raw record. T h e s p e c t r a s h o w n are t h e t i m e averaged s p e c t r a c o m p u t e d w i t h a t a p e r of 20%. T h e r e were slight d i f f e r e n c e s b e t w e e n t h e t i m e averaged s p e c t r a a n d t h e f r e q u e n c y s m o o t h e d spectra. T h e m o s t n o t i c e a b l e is for t h e case w h e r e n o filter has b e e n applied. T h e side lobe s u p p r e s s i o n gained b y f r e q u e n c y s m o o t h i n g limits t h e s u b h a r m o n i c cont a m i n a t i o n to t h e first h a r m o n i c .
226
R. COPPOLA
Discussion For the example shown, the largest low frequency contribution is at 1 c/sec and below. This large activity causes some elevation of the spectrum estimates at 2--4 c/sec. This elevation is most probably due to leakage. If a filter parameter of 0.95 is used only the activity below 1 c/sec is removed leaving the 2--4 c/sec range unchanged (see Fig. 1). Thus if accurate spectrum estimates in the 2--4 c/sec region were sought, it would seem desirable to filter the lower activity in this case. Whether the activity below 1 c/sec is physiological or electronic is uncertain; however, it would seem more likely to attribute it to the amplifier. Application of this filter technique will give investigators control over the low frequency content of their EEG data and will hopefully lead to a resolution of this issue.
Summary A method for isolating and removing very slow activity from EEG records prior to spectrum computation is presented. This is accomplished by an autoregressive filter whose frequency transfer characteristic can be easily controlled. Examples of spectra computed for various filter parameter values are shown.
R~sum~
Isolement de l'activitd d l'analyse spectrale de I'EEG
basse
frgquence
dans
L'auteur d6crit une m~thode permettant d'isoler et de supprimer des ondes tr~s lentes sur les enregistre-
ments EEG avant analyse spectrale. Ceci est r~alis6 au moyen d'un filtre autor~gressif dont la caract6ristique de transfert de fr~quence peut 6tre facilement control6e. Des exemples de spectres corrig6s pour diverses valeurs de filtre sont donn6s en illustration.
References Btackman, R.B. and Tukey, J.W. The Measurement of Power Spectra. Dover Publications, New York, 1958, p. 48. Coppola, R. A table driven system for stimulus response experiments. Proc. Digital Equipment User's Soc., 1977, 3: 1219--1222. Dumermuth, G. and Keller, E. EEG spectral analysis by means of Fast Fourier Transform. In: P. Kellaway and I. Petersen (Eds.), Automation of Clinical Electroencephalography. Raven Press, New York, 1973: 145--160. Kellaway, P. and Petersen, I. (Eds.) Automation of Clinical Electroencephalography. Raven Press, New York, 1973. Orr, W.C. and Hoffman, H.J. A 90-min cardiac biorhythm: methodology and data analysis using modified periodograms a n d complex demodulation. IEEE Trans. biomed. Engng, 1974, BME-21 : 130--143. Otnes, R.K. and Enochson, I. Digital Time Series Analysis. Wiley, New York, 1972. Walter, D.O. Spectral analysis for electroencephalograms: mathematical determination of neurophysiological relationships from records of limited duration. Exp. Neurol., 1963, 8: 155--181.