- Email: [email protected]
AR spectral analysis of EEG signals by using maximum likelihood estimation
Computers in Biology and Medicine 31 (2001) 441–450 www.elsevier.com/locate/compbiomed
AR spectral analysis of EEG signals by using maximum likelihood estimation ˙Inan G,ulera;∗ , M. Kemal Kiymikb , Mehmet Akinc , Ahmet Alkanb a
Department of Electronics–Computer Education, Faculty of Technical Education, Gazi University, 06500 Teknikokullar, Ankara, Turkey b Department of Electrical–Electronics Engineering, Kahramanmaras$ S&utc$u& I˙mam University, 46001 Kahramanmaras$, Turkey c Department of Electrical–Electronics Engineering, Dicle University, Diyarbak*r, Turkey Received 3 April 2001; accepted 4 June 2001
Abstract In this study, EEG signals were analyzed using autoregressive (AR) method. Parameters in AR method were realized by using maximum likelihood estimation (MLE). Results were compared with fast Fourier transform (FFT) method. It is observed that AR method gives better results in the analysis of EEG signals. On the other hand, the results have also showed that AR method can also be used for some other researches and diagnosis c 2001 Elsevier Science Ltd. All rights reserved. of diseases. Keywords: EEG; Autoregressive method (AR); Maximum likelihood estimation (MLE); Fast Fourier transform (FFT)
1. Introduction Power spectrum estimation methods have been used for a long time in the analysis of biological signals. From a historical perspective, it starts with the paper by Robinson and the book by Marple [1]. The classical power spectral estimation methods are mainly based on the fast Fourier transform (FFT), originally introduced by Schuster (1898). The modern model-based or parametric methods were originated by Yule (1927). These two methods were subsequently developed and applied by Walker (1931), Barlett (1948), Parzen (1957), Blackman and Tukey (1958), Burg (1967), and others [1]. EEG signals involve a great deal of information about the function of the brain. Since the classiBcation and evaluation of these signals are limited, there is no deBnite criterion evaluated by the experts for the visual analysis of EEG signals. For example, in the case of dominant alpha activity, ∗
Corresponding author. Tel.: +90-312-212-3976; fax: +90-312-212-0059. E-mail address: [email protected] (˙I. G,uler).
c 2001 Elsevier Science Ltd. All rights reserved. 0010-4825/01/$ - see front matter PII: S 0 0 1 0 - 4 8 2 5 ( 0 1 ) 0 0 0 2 2 - 1
442
I.˙ G&uler et al. / Computers in Biology and Medicine 31 (2001) 441–450
delta and teta activities have not been noticed. Since the routine clinical diagnosis needs to analyze EEG signals, some automation and computer techniques have been used for this purpose. To do that, EEG signals have been analyzed by using diGerent techniques [2,3]. Due to the random Iuctuations of EEG signals, the statistical average characteristics of random signals should be obtained. In particular, the autocorrelation function of a random process is the appropriate statistical average that will be used for characterizing random signals in the time domain. The Fourier transform of the autocorrelation function, which yields the power density spectrum, provides the transformation from the time domain to the frequency domain and estimation of signal variance as a function of frequency. Non-parametric power spectrum estimation methods can be calculated using FFT and understood easily by comparing with the parametric methods. Since FFT method needs long duration data records for suitable frequency resolution, it is applied to windowed data sets which assume all data zero except the data in the window. On the other hand, some spectral losses occur because of windowing process and these losses mask weak signals in the data [4]. Parametric power spectrum estimation methods reduce the spectral losses and give better frequency resolution. Also autoregressive (AR) method has an advantage over FFT in the case of shorter duration data records [5]. In order to model a signal, the properties of this signal must be taken into account. For example, AR model is suitable for the signals which have sudden peaks in their frequency spectrums. Inversely, moving average (MA) model is suitable for signals which have no sharp peaks in their frequency spectrums. On the other hand, the autoregressive-moving average (ARMA) model is suitable for signals which have both type of peaks in their frequency spectrums. Since EEG signals contain peaks at some frequencies, AR model can be used by employing Burg or Levinson–Durbin algorithm [6,7]. In this study, AR spectral analysis of EEG signals using maximum likelihood estimation (MLE) is studied. Data recorded from healthy and epileptic patients are processed using FFT and AR method with MLE in order to obtain the best method in the process. 2. Materials and method 2.1. Spectral analysis of EEG signals EEG signals are analyzed by using spectral analysis methods to diagnose some cerebral diseases. The power spectral density of the signal P(f) is found by applying conventional and modern spectral analysis methods such as FFT and AR. The data acquisition system for the processing of EEG signal is shown in Fig. 1. EEG signals used in this study are taken from Medical Faculty Hospital of Dicle University. These signals belong to several healthy and unhealthy (epileptic patients) persons. The signals are collected by a data acquisition system which contains data acquisition card (PCI MIO-16-E+ type), signal processors and a personnel computer. Data can be taken into computer memory quickly by using this card which is connected to PCI data bus of the computer [8]. For this system C based graphical programming language (LabVIEW) is used. The system provides real time data processing. Discrete Fourier transform (DFT) of a discrete signal x(n) is deBned as follows: N −1 −j2k Xk = ; (1) x(n) exp N n=0
I.˙ G&uler et al. / Computers in Biology and Medicine 31 (2001) 441–450
443
Fig. 1. A scheme of the EEG data acquisition system.
where Xk is expressed as the discrete Fourier coePcient, N is the frame size and x(n) is the input signal on the time domain. To obtain the frequency spectrum of this signal, logarithmic values of the squares of absolute values of Xk are found as follows: P(k) = 10 log |X |2 :
(2)
In the AR modeling method, the amplitude of a signal at a given period is obtained by summing up the diGerent amplitudes of previous samples, and adding the estimation error. The AR coePcients, which identify the amplitude rates, can be calculated by using the Burg or Levinson–Durbin algorithms. In the Levinson–Durbin algorithm, AR parameters are found out by solving the Yule– Walker equations. However, in the Burg algorithm, AR coePcients are found out with advanced reversed errors, which depend on samples taken from the signal. The order of the model, namely the Blter, depends on the number of AR coePcients. In the AR method, the model order is identiBed according to diGerent criteria. In this study, Akaike information criteria (AIC) are taken as the base. After inspecting Refs. [7,9], the model order p = 10 was taken because the determined model order was lower. In the AR model, MLE is used for the solution of the Yule–Walker equations to get AR model parameters [9]. To obtain stable and high performance AR model, some factors must be taken into consideration such as: selection of the right algorithm according to the model, selection of the model order, the length of the signal which will be modeled, and the level of stationary of the data. The order of the model is very important. If the selected order is low, there will be no deBnite peaks in the spectrum. So, the frequency details of the signal cannot be identiBed. If the selected order is very high, faulty peaks may be seen in the frequency spectrum which is not related to the original signal. 2.2. Autoregressive parameter estimation and MLE In the AR model, to Bnd out model parameters of Levinson–Durbin algorithm which makes use of the solution of the Yule–Walker equations is used. Autocorrelation estimation is used for the solution of these equations. After those autocorrelation, AR model parameters are estimated. To do that, biased form of the autocorrelation estimation is used which is given as rxx =
N − m− 1 1 ∗ x (n)x(n + m); N n=0
m ¿ 0:
(3)
I.˙ G&uler et al. / Computers in Biology and Medicine 31 (2001) 441–450
444
The purpose now is to estimate the AR model parameters by using MLE in the solution of the Yule–Walker equations from a record of EEG data. If the maximum likelihood estimate of a parameter exists under regular condition, it is consistent, asymptotically unbiased, ePcient, and normally distributed. Unfortunately, the maximum likelihood (ML) estimator is often too cumbersome to obtain [10]. As this is the case for the EEG model, it is proposed to estimate the model parameters by maximizing an approximation of the log-likelihood function, known as Whittle’s approximation, the derived estimator is expected to retain the properties associated with the ML estimator in an asymptotic sense, but with much less complexity [4]. In fact, Whittle’s estimate asymptotically retains the properties of the ML estimate for Gaussian random processes, but this is not generally true for the non-Gaussian case [2]. In many cases, it is diPcult to evaluate the MLE of the parameter its power spectrum density function (PSDF) is Gaussian due to the need to invert a large dimension covariance matrix. For example, if x ∼ N(0; c()), the MLE of is obtained by maximizing 1 t −1 P(x; ) = (4) e[−1=2X c ()X ] : n=2 (2) det(c()) If the covariance matrix cannot be inverted in the closed form, then a search technique will require inversion of the N ×N matrix for each value of to be searched. An alternative approximate method can be applied when x is the data from a zero mean random process, so that covariance matrix is Toeplitz. In such a case, the asymptotic log-likelihood function is given by N 1=2 N I (f) df; (5) ln P(x; ) = ln 2 − ln Pxx (f) + 2 2 −1=2 Pxx (f) where 1 I (f) = N
2 −1 N x(n)e−(j2fn) n=0
is the periodogram of the data and Pxx (f) is the power spectral density (PSD). The dependence of the the log-likelihood function on is through the PSD. DiGerentiation of (5) produces the necessary conditions for MLE @ ln P(x; ) I (f) @Pxx (f) 1 N 1=2 − 2 =− df (6) @i 2 −1=2 Pxx (f) Pxx (f) @i or 1=2 I (f) @Pxx (f) 1 − 2 df = 0: (7) Pxx (f) @i −1=2 Pxx (f) The second derivative allows the Newton–Raphson or scoring method to be implemented using the asymptotic likelihood function. This leads to simpler iterative procedures and is commonly used in practice. In this study, to Bnd MLE, the asymptotic form of the log-likelihood given by (5) is used. Since the PSD is Pxx (f) =
2u |A(f)|2
(8)
I.˙ G&uler et al. / Computers in Biology and Medicine 31 (2001) 441–450
after some calculations and derivations, the estimated autocorrelation function is [3,4,10] N −1−|k | 1 x(n)x(n + |k|); |k| 6 N − 1; Rˆ xx (k) = N n=0 0; |k| ¿ N
445
(9)
and the set of equations to be solved for the approximate MLE of the AR Blter parameters becomes p a(l) ˆ Rˆ xx (k − l) = − Rˆ xx (k); k = 1; 2; : : : ; p (10) l=1
or in matrix form ˆ Rxx (0) Rˆ xx (1) Rˆ (1) ··· xx .. . ··· ˆ ˆ Rxx (p − 1) Rxx (p − 2)
ˆ Rxx (1) · · · Rˆ xx (p − 1) a(1) ˆ ˆ Rˆ (2) · · · Rˆ xx (1) a(2) xx .. .. : . . ··· ··· ˆ a(p) ˆ Rˆ xx (p) · · · Rxx (0)
(11)
These are the so-called estimated Yule–Walker equations and this is the autocorrelation method of linear prediction. Note that the special form of the matrix and the right-hand vector, which thereby allow a recursive solution known as the Levinson recursion [1]. To complete the discussion explicit form for the MLE of 2u should be determined. From (8) p p 2 ˆu = a[k]R ˆ [ − k] = a[k]R ˆ (12) xx xx [k]: k=0
k=0
These parameters are used to compute the AR spectral power density as YW Pxx (f)
2 ˆwp : −j2fk |2 |1 + pk=1 a(k)e ˆ
(13)
3. Results and discussions In this study, analysis of three real and one simulated EEG signals is presented. For this aim, a software is developed by using laboratory virtual instrument engineering workbench (Lab View) Programming Language. These signals are analyzed by using FFT and AR method which have the MLE optimization. Results of the spectrums of these methods are represented below and compared with each other. In Fig. 2 an epileptic EEG signal and FFT of this signal are given. If the frequency spectrum of FFT is examined, it is seen that there are peaks at 1 and 3 Hz. AR spectrum of the same signal is presented in Fig. 3. There are peaks at 3 Hz with higher amplitude, 6, 9.5 and 13:5 Hz. When we compare these two spectrums it is seen AR spectrum has got sharper peaks and less misleading peaks than FFT. Due to this better frequency solution, explanation and determination of the activities in the signal is easier by using AR method. Since the signal is taken from an epileptic patient, the results Bt with the typical characteristics of epilepsy that is delta activity (low frequency range) [2]. Fig. 4 shows an EEG signal taken from a child and FFT of this signal. Spectrum of the AR method is given in Fig. 5. If these two spectrums are examined although FFT spectrum has got
446
I.˙ G&uler et al. / Computers in Biology and Medicine 31 (2001) 441–450
Fig. 2. Epileptic EEG signal and its FFT.
Fig. 3. AR spectrum of epileptic EEG signal.
wide and misleading peaks AR spectrum has got sharp and clear peaks. If these spectrums are examined, the delta activity (1–3 Hz), the alpha activity (6 –9 Hz), and the beta activity (11–13 Hz) can easily be seen. These results are true because it is a normal EEG signal. Higher variations and misleading peaks in FFT spectrum avoid the dominant alpha and delta activities. Fig. 6 shows an EEG signal of a healthy person and its FFT spectrum. FFT spectrum has peaks at 3.5, 8 and 10 Hz. On the other hand, we can see nearly the same peaks in Fig. 7 but more clearly than FFT spectrum. Fig. 8 shows a simulated EEG signal and its FFT spectrum. The signal contains four sinusoidal signals (at frequencies 2, 6, 9 and 15 Hz) and white noise. AR spectrum which uses the MLE is shown in Fig. 9. When these spectrums are compared with each other, it is seen that in the AR spectrum frequency content of the signal can be seen more clearly than FFT spectrum.
I.˙ G&uler et al. / Computers in Biology and Medicine 31 (2001) 441–450
Fig. 4. Normal EEG signal taken from a three-year old child and its FFT.
Fig. 5. AR spectrum of EEG signal.
Fig. 6. EEG signal taken from a 29-year old healthy person and its FFT.
447
448
I.˙ G&uler et al. / Computers in Biology and Medicine 31 (2001) 441–450
Fig. 7. AR spectrum of EEG signal.
Fig. 8. Simulated EEG signal and its FFT.
Fig. 9. AR spectrum of EEG signal.
I.˙ G&uler et al. / Computers in Biology and Medicine 31 (2001) 441–450
449
4. Conclusions In this study, the FFT and AR spectrums which have MLE optimization of real and simulated EEG signals are plotted and compared. To get AR method model parameters, MLE which has wide applications in statistics is used. After calculations and observations, it is seen that AR method spectrums have sharp and clear peaks. That is, AR spectrums show the frequency content of the signals more clearly than FFT spectrums and these results can be used for diagnosis of pathological events.
References [1] J.G. Proakis, D.G. Manolakis, Digital Signal Processing Principles, Algorithms, and Applications, Prentice-Hall, Englewood CliGs, NJ, 1996. [2] M. Zoubir, B. Boashash, Seizure detection of newborn EEG using a model approach, IEEE Trans. Biomed. Eng. 45(6) (1998) 477–489. [3] R. Srinivasan, L. Nunez, R. Silberstein, Spatial Bltering and neocortical dynamics: estimates of EEG coherence, IEEE Trans. Biomed. Eng. 45(7) (1998) 38–45. [4] K.O. Dhaparidze, A.M. Yaglom, Spectrum parameter estimation in time series analysis, in: P.R. Krishnaiah (Ed.), Developments in Statistics, Vol. 4, Academic Press, New York, 1983, pp. 1–96 (Chapter 1). [5] Y. Bard, Nonlinear Parameter Estimation, Academic Press, New York, 1974. [6] N.F. G,uler, M.K. KTymik, I.˙ G,uler, Comparison of FFT and AR-based sonogram outputs of 20 MHz pulsed Doppler data in real time, Comput. Biol. Med. 25(1995) 383–391. [7] D. Burshtein, E. Weinstein, Some relations between the various criteria for AR model order determination, IEEE Trans. Acoust. Speech Signal Process. ASSP-34=4 (1985) 1017–1019. [8] Instrumentation Reference and Catalogue, National Instruments, 1997. [9] B. Friedlander, B. Porat, A general lower bound for parametric spectrum estimation, IEEE Trans. Acoust. Speech Signal Process. ASSP-32 (1984) 728–733. [10] A.P. Dempster, D.B. Rubin, Maximum likelihood from incomplete data via the EM algorithm, Ann. Roy 39 (1977) 1–38. ˙ Inan G uler graduated from Erciyes University in 1981. He took his M.S. degree from Middle East Technical University ˙ in 1985, and his Ph.D. degree from Technical University of Istanbul in 1990, all in Electronic Engineering. Between 1988 and 1989 he worked in the Department of Medical Physics and Clinical Engineering at Leicester Royal InBrmary, where he studied the experimental part of his Ph.D. He is a Professor at Gazi University where he leads the Biomedical Engineering Group and is the Head of the Department. His areas of interest are biomedical systems, biomedical signal processing, biomedical instrumentation, and electronic circuit design. He has written more than 60 articles on biomedical engineering. M. Kemal Kiymik graduated from Erciyes University in 1982, gained his M.S. degree from Middle East Technical University in 1986, and his Ph.D. degree from Erciyes University in 1989, all in Electronic Engineering. Between 1991 and 1992 he worked as a research fellow at Rensselaer Polytechnic Institute where he carried out his post-doctoral studies. He ˙ is an Associate Professor at the Department of Electronic Engineering in KahramanmaraUs S,utUcu, Imam University where he is the head of department. His areas of interest are signal processing, and linear and nonlinear systems. He has numerous articles in the Beld of biomedical signal processing. Mehmet Akin graduated from Yildiz University in 1988, gained his M.S. and Ph.D. degree from Erciyes University in 1991 and 1995, respectively, all in Electronic Engineering. He is an Assistant Professor at the Department of Electronic
450
I.˙ G&uler et al. / Computers in Biology and Medicine 31 (2001) 441–450
Engineering in Dicle University where he is the head of the department. His area of interest is biomedical signal processing. He has numerous articles on the related areas. Ahmet Alkan graduated from Middle East Technical University in 1995, gained his M.S. degree from KahramanmaraUs ˙ S,utUcu, Imam University 1998, all in Electronic Engineering. He is working as an instructor in the Department of Electronic ˙ Engineering in KahramanmaraUs S,utUcu, Imam University. His area of interest is biomedical signal processing.
AR spectral analysis of EEG signals by using maximum likelihood estimation ˙Inan G,ulera;∗ , M. Kemal Kiymikb , Mehmet Akinc , Ahmet Alkanb a
Department of Electronics–Computer Education, Faculty of Technical Education, Gazi University, 06500 Teknikokullar, Ankara, Turkey b Department of Electrical–Electronics Engineering, Kahramanmaras$ S&utc$u& I˙mam University, 46001 Kahramanmaras$, Turkey c Department of Electrical–Electronics Engineering, Dicle University, Diyarbak*r, Turkey Received 3 April 2001; accepted 4 June 2001
Abstract In this study, EEG signals were analyzed using autoregressive (AR) method. Parameters in AR method were realized by using maximum likelihood estimation (MLE). Results were compared with fast Fourier transform (FFT) method. It is observed that AR method gives better results in the analysis of EEG signals. On the other hand, the results have also showed that AR method can also be used for some other researches and diagnosis c 2001 Elsevier Science Ltd. All rights reserved. of diseases. Keywords: EEG; Autoregressive method (AR); Maximum likelihood estimation (MLE); Fast Fourier transform (FFT)
1. Introduction Power spectrum estimation methods have been used for a long time in the analysis of biological signals. From a historical perspective, it starts with the paper by Robinson and the book by Marple [1]. The classical power spectral estimation methods are mainly based on the fast Fourier transform (FFT), originally introduced by Schuster (1898). The modern model-based or parametric methods were originated by Yule (1927). These two methods were subsequently developed and applied by Walker (1931), Barlett (1948), Parzen (1957), Blackman and Tukey (1958), Burg (1967), and others [1]. EEG signals involve a great deal of information about the function of the brain. Since the classiBcation and evaluation of these signals are limited, there is no deBnite criterion evaluated by the experts for the visual analysis of EEG signals. For example, in the case of dominant alpha activity, ∗
Corresponding author. Tel.: +90-312-212-3976; fax: +90-312-212-0059. E-mail address: [email protected] (˙I. G,uler).
c 2001 Elsevier Science Ltd. All rights reserved. 0010-4825/01/$ - see front matter PII: S 0 0 1 0 - 4 8 2 5 ( 0 1 ) 0 0 0 2 2 - 1
442
I.˙ G&uler et al. / Computers in Biology and Medicine 31 (2001) 441–450
delta and teta activities have not been noticed. Since the routine clinical diagnosis needs to analyze EEG signals, some automation and computer techniques have been used for this purpose. To do that, EEG signals have been analyzed by using diGerent techniques [2,3]. Due to the random Iuctuations of EEG signals, the statistical average characteristics of random signals should be obtained. In particular, the autocorrelation function of a random process is the appropriate statistical average that will be used for characterizing random signals in the time domain. The Fourier transform of the autocorrelation function, which yields the power density spectrum, provides the transformation from the time domain to the frequency domain and estimation of signal variance as a function of frequency. Non-parametric power spectrum estimation methods can be calculated using FFT and understood easily by comparing with the parametric methods. Since FFT method needs long duration data records for suitable frequency resolution, it is applied to windowed data sets which assume all data zero except the data in the window. On the other hand, some spectral losses occur because of windowing process and these losses mask weak signals in the data [4]. Parametric power spectrum estimation methods reduce the spectral losses and give better frequency resolution. Also autoregressive (AR) method has an advantage over FFT in the case of shorter duration data records [5]. In order to model a signal, the properties of this signal must be taken into account. For example, AR model is suitable for the signals which have sudden peaks in their frequency spectrums. Inversely, moving average (MA) model is suitable for signals which have no sharp peaks in their frequency spectrums. On the other hand, the autoregressive-moving average (ARMA) model is suitable for signals which have both type of peaks in their frequency spectrums. Since EEG signals contain peaks at some frequencies, AR model can be used by employing Burg or Levinson–Durbin algorithm [6,7]. In this study, AR spectral analysis of EEG signals using maximum likelihood estimation (MLE) is studied. Data recorded from healthy and epileptic patients are processed using FFT and AR method with MLE in order to obtain the best method in the process. 2. Materials and method 2.1. Spectral analysis of EEG signals EEG signals are analyzed by using spectral analysis methods to diagnose some cerebral diseases. The power spectral density of the signal P(f) is found by applying conventional and modern spectral analysis methods such as FFT and AR. The data acquisition system for the processing of EEG signal is shown in Fig. 1. EEG signals used in this study are taken from Medical Faculty Hospital of Dicle University. These signals belong to several healthy and unhealthy (epileptic patients) persons. The signals are collected by a data acquisition system which contains data acquisition card (PCI MIO-16-E+ type), signal processors and a personnel computer. Data can be taken into computer memory quickly by using this card which is connected to PCI data bus of the computer [8]. For this system C based graphical programming language (LabVIEW) is used. The system provides real time data processing. Discrete Fourier transform (DFT) of a discrete signal x(n) is deBned as follows: N −1 −j2k Xk = ; (1) x(n) exp N n=0
I.˙ G&uler et al. / Computers in Biology and Medicine 31 (2001) 441–450
443
Fig. 1. A scheme of the EEG data acquisition system.
where Xk is expressed as the discrete Fourier coePcient, N is the frame size and x(n) is the input signal on the time domain. To obtain the frequency spectrum of this signal, logarithmic values of the squares of absolute values of Xk are found as follows: P(k) = 10 log |X |2 :
(2)
In the AR modeling method, the amplitude of a signal at a given period is obtained by summing up the diGerent amplitudes of previous samples, and adding the estimation error. The AR coePcients, which identify the amplitude rates, can be calculated by using the Burg or Levinson–Durbin algorithms. In the Levinson–Durbin algorithm, AR parameters are found out by solving the Yule– Walker equations. However, in the Burg algorithm, AR coePcients are found out with advanced reversed errors, which depend on samples taken from the signal. The order of the model, namely the Blter, depends on the number of AR coePcients. In the AR method, the model order is identiBed according to diGerent criteria. In this study, Akaike information criteria (AIC) are taken as the base. After inspecting Refs. [7,9], the model order p = 10 was taken because the determined model order was lower. In the AR model, MLE is used for the solution of the Yule–Walker equations to get AR model parameters [9]. To obtain stable and high performance AR model, some factors must be taken into consideration such as: selection of the right algorithm according to the model, selection of the model order, the length of the signal which will be modeled, and the level of stationary of the data. The order of the model is very important. If the selected order is low, there will be no deBnite peaks in the spectrum. So, the frequency details of the signal cannot be identiBed. If the selected order is very high, faulty peaks may be seen in the frequency spectrum which is not related to the original signal. 2.2. Autoregressive parameter estimation and MLE In the AR model, to Bnd out model parameters of Levinson–Durbin algorithm which makes use of the solution of the Yule–Walker equations is used. Autocorrelation estimation is used for the solution of these equations. After those autocorrelation, AR model parameters are estimated. To do that, biased form of the autocorrelation estimation is used which is given as rxx =
N − m− 1 1 ∗ x (n)x(n + m); N n=0
m ¿ 0:
(3)
I.˙ G&uler et al. / Computers in Biology and Medicine 31 (2001) 441–450
444
The purpose now is to estimate the AR model parameters by using MLE in the solution of the Yule–Walker equations from a record of EEG data. If the maximum likelihood estimate of a parameter exists under regular condition, it is consistent, asymptotically unbiased, ePcient, and normally distributed. Unfortunately, the maximum likelihood (ML) estimator is often too cumbersome to obtain [10]. As this is the case for the EEG model, it is proposed to estimate the model parameters by maximizing an approximation of the log-likelihood function, known as Whittle’s approximation, the derived estimator is expected to retain the properties associated with the ML estimator in an asymptotic sense, but with much less complexity [4]. In fact, Whittle’s estimate asymptotically retains the properties of the ML estimate for Gaussian random processes, but this is not generally true for the non-Gaussian case [2]. In many cases, it is diPcult to evaluate the MLE of the parameter its power spectrum density function (PSDF) is Gaussian due to the need to invert a large dimension covariance matrix. For example, if x ∼ N(0; c()), the MLE of is obtained by maximizing 1 t −1 P(x; ) = (4) e[−1=2X c ()X ] : n=2 (2) det(c()) If the covariance matrix cannot be inverted in the closed form, then a search technique will require inversion of the N ×N matrix for each value of to be searched. An alternative approximate method can be applied when x is the data from a zero mean random process, so that covariance matrix is Toeplitz. In such a case, the asymptotic log-likelihood function is given by N 1=2 N I (f) df; (5) ln P(x; ) = ln 2 − ln Pxx (f) + 2 2 −1=2 Pxx (f) where 1 I (f) = N
2 −1 N x(n)e−(j2fn) n=0
is the periodogram of the data and Pxx (f) is the power spectral density (PSD). The dependence of the the log-likelihood function on is through the PSD. DiGerentiation of (5) produces the necessary conditions for MLE @ ln P(x; ) I (f) @Pxx (f) 1 N 1=2 − 2 =− df (6) @i 2 −1=2 Pxx (f) Pxx (f) @i or 1=2 I (f) @Pxx (f) 1 − 2 df = 0: (7) Pxx (f) @i −1=2 Pxx (f) The second derivative allows the Newton–Raphson or scoring method to be implemented using the asymptotic likelihood function. This leads to simpler iterative procedures and is commonly used in practice. In this study, to Bnd MLE, the asymptotic form of the log-likelihood given by (5) is used. Since the PSD is Pxx (f) =
2u |A(f)|2
(8)
I.˙ G&uler et al. / Computers in Biology and Medicine 31 (2001) 441–450
after some calculations and derivations, the estimated autocorrelation function is [3,4,10] N −1−|k | 1 x(n)x(n + |k|); |k| 6 N − 1; Rˆ xx (k) = N n=0 0; |k| ¿ N
445
(9)
and the set of equations to be solved for the approximate MLE of the AR Blter parameters becomes p a(l) ˆ Rˆ xx (k − l) = − Rˆ xx (k); k = 1; 2; : : : ; p (10) l=1
or in matrix form ˆ Rxx (0) Rˆ xx (1) Rˆ (1) ··· xx .. . ··· ˆ ˆ Rxx (p − 1) Rxx (p − 2)
ˆ Rxx (1) · · · Rˆ xx (p − 1) a(1) ˆ ˆ Rˆ (2) · · · Rˆ xx (1) a(2) xx .. .. : . . ··· ··· ˆ a(p) ˆ Rˆ xx (p) · · · Rxx (0)
(11)
These are the so-called estimated Yule–Walker equations and this is the autocorrelation method of linear prediction. Note that the special form of the matrix and the right-hand vector, which thereby allow a recursive solution known as the Levinson recursion [1]. To complete the discussion explicit form for the MLE of 2u should be determined. From (8) p p 2 ˆu = a[k]R ˆ [ − k] = a[k]R ˆ (12) xx xx [k]: k=0
k=0
These parameters are used to compute the AR spectral power density as YW Pxx (f)
2 ˆwp : −j2fk |2 |1 + pk=1 a(k)e ˆ
(13)
3. Results and discussions In this study, analysis of three real and one simulated EEG signals is presented. For this aim, a software is developed by using laboratory virtual instrument engineering workbench (Lab View) Programming Language. These signals are analyzed by using FFT and AR method which have the MLE optimization. Results of the spectrums of these methods are represented below and compared with each other. In Fig. 2 an epileptic EEG signal and FFT of this signal are given. If the frequency spectrum of FFT is examined, it is seen that there are peaks at 1 and 3 Hz. AR spectrum of the same signal is presented in Fig. 3. There are peaks at 3 Hz with higher amplitude, 6, 9.5 and 13:5 Hz. When we compare these two spectrums it is seen AR spectrum has got sharper peaks and less misleading peaks than FFT. Due to this better frequency solution, explanation and determination of the activities in the signal is easier by using AR method. Since the signal is taken from an epileptic patient, the results Bt with the typical characteristics of epilepsy that is delta activity (low frequency range) [2]. Fig. 4 shows an EEG signal taken from a child and FFT of this signal. Spectrum of the AR method is given in Fig. 5. If these two spectrums are examined although FFT spectrum has got
446
I.˙ G&uler et al. / Computers in Biology and Medicine 31 (2001) 441–450
Fig. 2. Epileptic EEG signal and its FFT.
Fig. 3. AR spectrum of epileptic EEG signal.
wide and misleading peaks AR spectrum has got sharp and clear peaks. If these spectrums are examined, the delta activity (1–3 Hz), the alpha activity (6 –9 Hz), and the beta activity (11–13 Hz) can easily be seen. These results are true because it is a normal EEG signal. Higher variations and misleading peaks in FFT spectrum avoid the dominant alpha and delta activities. Fig. 6 shows an EEG signal of a healthy person and its FFT spectrum. FFT spectrum has peaks at 3.5, 8 and 10 Hz. On the other hand, we can see nearly the same peaks in Fig. 7 but more clearly than FFT spectrum. Fig. 8 shows a simulated EEG signal and its FFT spectrum. The signal contains four sinusoidal signals (at frequencies 2, 6, 9 and 15 Hz) and white noise. AR spectrum which uses the MLE is shown in Fig. 9. When these spectrums are compared with each other, it is seen that in the AR spectrum frequency content of the signal can be seen more clearly than FFT spectrum.
I.˙ G&uler et al. / Computers in Biology and Medicine 31 (2001) 441–450
Fig. 4. Normal EEG signal taken from a three-year old child and its FFT.
Fig. 5. AR spectrum of EEG signal.
Fig. 6. EEG signal taken from a 29-year old healthy person and its FFT.
447
448
I.˙ G&uler et al. / Computers in Biology and Medicine 31 (2001) 441–450
Fig. 7. AR spectrum of EEG signal.
Fig. 8. Simulated EEG signal and its FFT.
Fig. 9. AR spectrum of EEG signal.
I.˙ G&uler et al. / Computers in Biology and Medicine 31 (2001) 441–450
449
4. Conclusions In this study, the FFT and AR spectrums which have MLE optimization of real and simulated EEG signals are plotted and compared. To get AR method model parameters, MLE which has wide applications in statistics is used. After calculations and observations, it is seen that AR method spectrums have sharp and clear peaks. That is, AR spectrums show the frequency content of the signals more clearly than FFT spectrums and these results can be used for diagnosis of pathological events.
References [1] J.G. Proakis, D.G. Manolakis, Digital Signal Processing Principles, Algorithms, and Applications, Prentice-Hall, Englewood CliGs, NJ, 1996. [2] M. Zoubir, B. Boashash, Seizure detection of newborn EEG using a model approach, IEEE Trans. Biomed. Eng. 45(6) (1998) 477–489. [3] R. Srinivasan, L. Nunez, R. Silberstein, Spatial Bltering and neocortical dynamics: estimates of EEG coherence, IEEE Trans. Biomed. Eng. 45(7) (1998) 38–45. [4] K.O. Dhaparidze, A.M. Yaglom, Spectrum parameter estimation in time series analysis, in: P.R. Krishnaiah (Ed.), Developments in Statistics, Vol. 4, Academic Press, New York, 1983, pp. 1–96 (Chapter 1). [5] Y. Bard, Nonlinear Parameter Estimation, Academic Press, New York, 1974. [6] N.F. G,uler, M.K. KTymik, I.˙ G,uler, Comparison of FFT and AR-based sonogram outputs of 20 MHz pulsed Doppler data in real time, Comput. Biol. Med. 25(1995) 383–391. [7] D. Burshtein, E. Weinstein, Some relations between the various criteria for AR model order determination, IEEE Trans. Acoust. Speech Signal Process. ASSP-34=4 (1985) 1017–1019. [8] Instrumentation Reference and Catalogue, National Instruments, 1997. [9] B. Friedlander, B. Porat, A general lower bound for parametric spectrum estimation, IEEE Trans. Acoust. Speech Signal Process. ASSP-32 (1984) 728–733. [10] A.P. Dempster, D.B. Rubin, Maximum likelihood from incomplete data via the EM algorithm, Ann. Roy 39 (1977) 1–38. ˙ Inan G uler graduated from Erciyes University in 1981. He took his M.S. degree from Middle East Technical University ˙ in 1985, and his Ph.D. degree from Technical University of Istanbul in 1990, all in Electronic Engineering. Between 1988 and 1989 he worked in the Department of Medical Physics and Clinical Engineering at Leicester Royal InBrmary, where he studied the experimental part of his Ph.D. He is a Professor at Gazi University where he leads the Biomedical Engineering Group and is the Head of the Department. His areas of interest are biomedical systems, biomedical signal processing, biomedical instrumentation, and electronic circuit design. He has written more than 60 articles on biomedical engineering. M. Kemal Kiymik graduated from Erciyes University in 1982, gained his M.S. degree from Middle East Technical University in 1986, and his Ph.D. degree from Erciyes University in 1989, all in Electronic Engineering. Between 1991 and 1992 he worked as a research fellow at Rensselaer Polytechnic Institute where he carried out his post-doctoral studies. He ˙ is an Associate Professor at the Department of Electronic Engineering in KahramanmaraUs S,utUcu, Imam University where he is the head of department. His areas of interest are signal processing, and linear and nonlinear systems. He has numerous articles in the Beld of biomedical signal processing. Mehmet Akin graduated from Yildiz University in 1988, gained his M.S. and Ph.D. degree from Erciyes University in 1991 and 1995, respectively, all in Electronic Engineering. He is an Assistant Professor at the Department of Electronic
450
I.˙ G&uler et al. / Computers in Biology and Medicine 31 (2001) 441–450
Engineering in Dicle University where he is the head of the department. His area of interest is biomedical signal processing. He has numerous articles on the related areas. Ahmet Alkan graduated from Middle East Technical University in 1995, gained his M.S. degree from KahramanmaraUs ˙ S,utUcu, Imam University 1998, all in Electronic Engineering. He is working as an instructor in the Department of Electronic ˙ Engineering in KahramanmaraUs S,utUcu, Imam University. His area of interest is biomedical signal processing.