Spatial Filtering to Detect Brain Sources from EEG Measurements P.Vergallo, A.Lay-Ekuakille, Senior Member,IEEE Department of Innovation Engineering University of Salento Lecce, Italy
[email protected]
Abstract— In clinical analysis, neural activity in the brain due to groups of neurons located in the head is recorded by means of EEG electrodes positioned on the scalp of the patient. The identification of the sources responsible for this brain activity is of great importance, especially if neurophysiological disorders are detected. The problem is that there is usually an unknown number of signals simultaneously recorded from the electrodes, each one from unknown direction and with unknown amplitudes, and the measurements are always corrupted by noise too. In this work a method of spatial filtering which combines DOA estimation and Beamforming techniques is presented: a DOA estimation method allows to estimate the directions of activation of the sources, while spatial filters of type Beamforming are designed to pass the electrical activity from a given direction and stop that originated from other directions. In this way the brain is explored to detect sources that have determined the recorded activity from EEG. Keywords— EEG, Beamforming, DOA Estimation, Sources Localization, Brain activity.
I. INTRODUCTION EEG is one of the main methods for the extraction of information from the brain for research and clinical purposes. It records the electrical activity of the brain by electrodes on the human scalp. EEG signal is a random process, but amplitudes or frequencies distribution of EEG signals depends on physiologically state of the patient. EEG signals can include the following waveform rhythms: α(813Hz), β(13-30Hz), ϑ (4-8Hz), δ(0.5-4Hz) [1]. EEG signal is related to the level of consciousness of the person. As the activity increases, the EEG shifts to the highest dominant frequency and lower amplitude. For clinical analysis it takes into account the fact that the rhythm α is index of normal tracing. Since its regularity and its amplitude are larger when the subject has closed eyes, the EEGs were recorded by placing the patient at rest. Two main tests are added to this trace, hyperpnoea and photic stimulation. Therefore the normality of the recorded tracks obtained by these tests is investigated. However in addition to these types of normal rhythms, pathological paroxysmal EEG figures can be present. These suddenly interrupt the underlying brain activity, and just suddenly ceased. They are indicative of a seizure. So, advanced techniques of signal processing, has been studied to analyze complex activities in
978-1-4799-2921-4/14/$31.00 ©2014 IEEE
S.Urooj, Member, EEE, V. Bhateja Dept of Electrical Eng./Dept of Electronics & Comm. Gautam Buddha University/SRM University Greater-Noida (U.P.)/Lucknow-Deva (U.P.) India
order to extract important information, especially when a seizure occurs [2]. Some of techniques concern a univariate time-series analysis from single- channel EEG recording, as Fourier Transform, Discrete wavelet transform [3], or multispectral analysis [4-5]. However multivariate timeseries techniques, which consists in taking more than one channel of EEG into account simultaneously, allows a major understanding of the interactions between the different components of the system under observation [6-7-8]. In particular, multivariate time-series analysis allows an estimation of spatial parameters necessary in the sources localization problem, especially when a seizure occurs. Seizures can belong to one of the following classes: Focal epilepsy: seizures begin within a limited space of the cerebral cortex, called the epileptogenic zone; Generalized epilepsy: the discharge of the seizure affects the entire cerebral cortex from the beginning; Multifocal epilepsy: more different epileptogenic zones are responsible for the onset of crisis. In particular the definition of the epileptogenic area is essential to undergo to surgery the patients whose seizures are refractory to drug therapy. In this work the EEG measurements obtained from multielectrodes on the scalp, in according to international measurement system 10/20, are used in order to estimate the location, direction and magnitude of the neural brain sources responsible for the recorded activity. The problem consists of two steps, to locate the brain sources by DOA (Direction of arrival) method, and to estimate the source signal identifying its amplitude and frequency by solving the socalled waveform estimation problem. This second step is carried out by Beamforming algorithm. Beamforming is a technique of spatial filtering, originally used in combination with a sensor [9] or antenna array [10]. It has all advantages of spatial processing, because it spatially samples the signal and uses the collected spatial information, linearly combined through appropriate weights, to preserve the signal originating from a desired direction, attenuating the signals from other ones. At the same way, spatial filters are designed in neurophysiology to pass the electrical activity from a given position and stop that originated from other
location [7] on the base the estimated directions by DOA estimation method [11-12]. The choice of coefficients to combine linearly the spatial samples depends on the type of Beamforming used. We estimate the directions of action of the sources of interest within the considered volume conductor by an improved version of the algorithm MUSIC (MUltiple Signal Classification) presented in [8], which calculates a spatial sparse solution by focusing the signal energy along only the directions of interest. The waveform of the estimated source is obtained by a Beamformer LCMP (Linearly Constrained Minimum Power). The key point of the proposed method is that the power in output from a spatial filter is the estimation of the neural power originated within the passband of the filter designed. So we define the passband of the our spatial filter through accurate evaluation of the neural power estimated with DOA estimation algorithm. The peaks identified in this map are associated with each potential brain source to be reconstruct by Beamforming technique. A validation of the method is reported considering simulated data, and afterwards real EEG data. II. MODEL FOR THE HEAD AND BRAIN SOURCE A simple spherical model for the head is used where the international measurement system 10/20 by electrodes placed on the scalp is mapped on it (Fig.1) [7], although more sophisticated models use a set of nested concentric spheres for brain, skull and scalp [13-14].
Fig. 2. The dipole parameters for a given current source and current sink configuration (left) and the dipole with its component (right)
So the dipole moment is:
d I 0 pa p
(1)
whit ̅ a unit vector in the direction of the dipole. It can be decomposed in three dipoles each oriented along x, y, z,
d dx ax d y a y dz az
(2)
ax, ay, az are the unity vectors along the three cartesian axes, while d=(dx,dy,dz) are the components of the dipole moment. So a potential P at an arbitrary point generated by a dipole at a position r can be defined as:
P( r , d ) d x P ( r , a x ) d y P ( r , a y ) d z P (r , a x )
(3)
More precisely, if the source is located in a point and if are the positions of the electrodes, considering a spherical surface, the potential P is:
P(rf , rd , d )
d ar 4 (rf rd )2
(4)
where σ is the conductivity of the head, and ̅ is the unit vector in the radial direction. So the distance field pointsource is:
r ( x x ')2 ( y y ')2 ( z z ')2
Fig. 1.
3D placement of electrodes in according to 10/20 system (left) and approximated considered head geometry (right) [7].
When analyzing EEG data, it is often desirable a more realistic model. Despite the simplification of the our model, it is possible to achieve good outcomes thanks to a good definition of a function which characterizes the sensitivity/energy distribution, named lead field matrix. The electrical activity in a limited area of the brain is modeled by a current dipole [15], shown in Fig.2 as two monopoles of opposite sign, but equal strength I0, separated by a very small distance p. It is defined in terms of its location and its orientation.
(5)
The voltage measured at the electrode is expressed as a function of the position and orientation of a unit dipole source defined by the lead field matrix L(rd) [15]. In fact, each electrode on the scalp detects the component of activation of the dipole along its sensitivity direction, and this information is contained in L(rd) matrix. As the medium is linear, the potential to the scalp is the superposition of the potential for many active neurons. Therefore the potential P due to L active sources to the positions qi, with i = 1, ..., L, is: L
P L(rd (qi ))d (qi ) n
(7)
i 1
where n is a noise component. The Eq.(7) represents the socalled solution to forward problem. In particular the concept of lead field is suitable for use of Beamforming techniques in our spatial sampling problem of the EEG recorded signal. As a scanning method, Beamforming examines each point in the considered volume in order to determine whether a presumable source at that point contributes to the recorded signal. This problem could be expensive, and a better solution is obtained estimating a priori the directions of activation of the sources of interest. In such way the Beamformer will point only along these directions to obtain
information concerning amplitude and frequency of the desired brain source. III. PROPOSED METHOD As aforementioned, the starting point for the definition of a spatial filter is the estimate of the directions of the active sources by DOA estimation algorithm [8] that offers good results in terms of spatial resolution. The used algorithm is an improved version of the classical method MUSIC [16], because it uses the sparsity of the field signal in order to focus the energy of the signal along the direction of origin of the sources [10]. The forward problem in Eq.(7) is redefined as P L(qi )s n (8) where L(qi) is the lead field matrix which contains the information related to each source dipole located in (xi,yi,zi) point of the considered volume conductor; s is the field signal matrix where the signals related to possible active source are defined, and it represents the sparse solution to be calculated; n is a Gaussian noise, and P is the matrix which contains EEG measurements. The sparsity of the solution is defined by considering that it contains a no null information at the correspondence of the (xi,yi,zi) point of active sources and zero otherwise. So it is imposed a priori in the problem, and two cost functions must be satisfied, one related to data fitting and other to obtain a sparse solution. The solution is given by solving the Tikhonov regularization problem, 2
J (s(k )) P(k ) L(qi ) s(k ) 2 2 s(k ) s0
2 2
s0
(s(k ))
2
1) 40.72, 21.77 2) 79.68, 119.17 3) 63.64, 90.52
(9)
where k is an index relative to each acquired sample from the electrodes on the scalp, and s0 is the a priori information about the sparsity of the solution to impose in the problem in this way: K
where TD is defined in order to consider the only D singular values, TD=[ID 0]T, with a DxD matrix identity and 0 a D (k-D) matrix of zeros. In this way, Tikhonov problem becomes D-dimensional because Pr is a matrix n D. The scalar λ in Eq.(9) is the regularization parameter balancing the tradeoff between the two costs, calculated through a graphical tool known as L-curve [10]. An initial solution for s0 is calculated, starting from least squares solution, and its components less than 50% with respect to maximum value of same s0 vector are set to zero. A new solution is calculated by solving Tikhonov problem (Eq.(9)), and the new output EEG is reconstructed on which MUSIC algorithm is applied. The procedure is repeated until a defined number of source is found. As result a brain power map is obtained, where the location of maximum power are associated with the direction of activation of the sources of interest. In Fig.3 the result obtained considering three sources is shown. The positions of the sources in terms of azimuth and elevation angles are:
Fig. 3. Sparse solution for DOA estimation problem
(10)
k 1
So, for each point of the volume conductor, the solutions calculated by solving the k Tikhonov problem are added together, and the sparsity of the signal is defined setting the smallest terms in s0 to zero. In this way s0 is such that it contains a great number of elements equal to zero and just few components are different from the null value. However, as Eq.(9) must be solved k times, that is for each sample, it is possible to reduce the computational effort by performing the Singular Value Decomposition (SVD) on the matrix P of EEG signals. So, if P , where n is the number of electrodes and k is the number of acquired samples,
P U V T
(11) with U unitary matrix, the matrix of the singular values of P, V unitary matrix. The D singular values, where D is the number of sources to estimate, corresponds to signal plus noise subspace, while the rest corresponds only to noise subspace. So the measurements are expressed as Pr U TD (12)
Fig. 4. MUSIC solution without sparsity
Spatial resolution (Fig.3) is better than the standard MUSIC method (Fig.4). After, these directions are given to Beamforming algorithm in order to reconstruct the waveform of the sources and obtained information about its frequencies and amplitudes. LCMP Beamformer is designed [7] by including linear constraints to preserve the signal in the desired direction and reduces contributions of noise and interferences in the others directions. With reference to Fig.5 the output of the Beamformer is:
Y wT P
(13)
Fig. 5. Beamformer applied on EEG measurement
In Eq.(13), P are the EEG measurements, w is the vector of the calculated coefficient from LCMP Beamformer. The output Y is a linear combination of the measurements through the appropriate weights w in order to obtain a spatial selectivity. For this aim, the following constrain must be satisfied:
Fig. 6. Brain power map
Beamforming algorithm uses these calculated directions to estimate amplitudes and frequencies of the sources (Fig.7).
wT L(q) g T
(14) where L(q) is the lead field matrix which contains the information for each direction of activation of the sources calculated through DOA estimation algorithm, and g is the gain vector such that a unity gain is associated with the estimated direction, while a zero gain corresponds to other directions. In this way, all sources of interest can be reconstructed. The weights are calculated by minimizing the power in output,
E[Y 2 ] wT Rx w
(15) under the constrain in Eq.(14). Rx in Eq.(15) is the covariance matrix of the measurements, estimated from the available data:
Rx
1 N ( xn x)( xn x)T N 1 i 1
(16)
where ̅ is the sample mean and , n=1,..,N, are the measurements from the electrodes. It is a constrained minimization problem that is solved by means of the Lagrange multipliers: wopt min[wT Rx w l (wT L(q) g H ) l* ( L(q) H w g )] (17) w
obtaining,
wopt g T [ L(q)T Rx1L(q)]1 L(q)T Rx1
(18)
IV. RESULTS At first time, the method is applied on the simulated signals obtained by solving the forward problem (Eq.(7)). For this aim, it is assumed that the activity recorded by the electrodes is due to three sources with known positions within of the volume conductor considered, with frequencies different: the first source has a frequency of 4Hz, the second 8Hz and the last is the sum of components with frequencies 4Hz and 8Hz. From the simulated EEG measurements, the direction of activation of the sources are calculated through the DOA estimation algorithm, applying Eq.(8-12). So the brain power map is obtained (Fig.6), where the peaks correspondence to the directions of the sources.
Fig. 7. Reconstructed sources
It is evident as the frequencies of the sources are well estimated, while the amplitudes present an estimate error which increases with the decrease of the Signal to Noise Ratio (SNR). This is due to fact that the noise introduces a non-uniform component in the map that depends upon the measured sites. The method is following applied on real EEG. In Fig.8, the EEG signals and the estimated sources on the base of the neural power map obtained in DOA estimation step are reported. The sum of these calculated components allows to obtained the EEG signals in according to Eq.(7) (Fig.9). So, by calculating all sources related to neural power map, important information about neural activity in terms of frequency and amplitude are obtained. In this way it is possible to associate each sources with the region of the head affected by this activity. In Fig.10 the neural power map in xy coordinates is reported.
Fig. 8. EEG signals recorded from some electrodes, neural power map calculated by DOA estimation algorithm and some calculated sources by Beamforming algorithm
yet necessary by considering discontinuities and inhomogeneities of the head. REFERENCES [1]
J.Delay, G. Verdeaux, J. Gaches, “Elettroencefalografia clinica”, Masson, 1982.
[2]
S.Tong, N.V.Takor,”Quantitative EEG analysis methods and clinical applications”, Artech House, 2009. N.Hazarika, J. Z. Chenb, A.C. Tsoic, A. Sergejewd,”Classification of EEG signals using the wavelet transform”, Signal Processing, Vol. 59, Issue 1, May 1997, pp. 61–72. A.Lay-Ekuakille, P. Vergallo, D.Caratelli, F. Conversano, S. Casciaro, A. Trabacca, “Multispectrum approach in quantitative EEG:Accuracy and physical effort”, IEEE Sensor Journal, Vol.13, Issue 9, pp.3331-3340, 2013. A.Lay-Ekuakille, P. Vergallo, G. Griffo, F. Conversano, S.Casciaro, S.Urooj, V.Bhateja, A.Trabacca, “Mutidimensional Analysis of EEG Features Using Advanced Spectral Estimates for Diagnosis Accuracy”, 2013 IEEE Memea, 4-5 May, 2013, Gatineau - Ottawa, Canada. N. Mammone, D. Labate, A. Lay-Ekuakille, F.C. Morabito, “Analysis of absence seizure generation using EEG spatial-temporal regularity measures”, International Journal of Neural Systems, Vol. 22, Issue 6, December 2012. P. Vergallo, A. Lay-Ekuakille, N.I. Giannoccaro, A. Massaro, S. Urooj, D. Caratelli, A.Trabacca, “Processing EEG Signals through Beamforming Techniques for Seizure Diagnosis”, ICST 2012, December 18 -22, 2012, Kolkata, India. P. Vergallo, A. Lay-Ekuakille, “Brain Source Localization: A New Method Based On Music Algorithm and Spatial Sparsity”, Rev. Sci. Instrum., vol. 84, 085117, 2013. A. Lay-Ekuakille, P. Vergallo, D. Saracino, A. Trotta,”Optimizing and Post Processing of a Smart Beamformer for Obstacle Retrieval”, IEEE Sensors Journal, vol.12,no.5, May 2012. P.Vergallo, A. Lay-Ekuakille, D. Caratelli, “Sparsity of the field signal-based method for improving spatial resolution in antenna sensor array processing”, Progress in Electromagnetics Research, Vol. 142, pp. 369-388, 2013. B.D.VanVeen, W.VanDrongelen, M.Yuchtman, A.Suzuki, ”Localization of Brain Electrical Activity via Linearly Constrained Minimum Variance Spatial Filtering ”, IEEE Transactions on Biomedical Engineering, vol. 44,no.9, pp.867880, September 1997 G.VanHoey, R.Van de Walle, B. Vanrumste, M. D’Havé, I. Lemahieu, P.Boon, “Beamforming Techniques Applied in EEG Source Analysis”, ProRISC Circuits, Systems and Signal Processing location, Mierlo, The Netherlands, 1999. Y. Salu, L. G. Cohen, D. Rose, S. Sato, C. Kufta, M. Hallett, “An improved method for Localizing Electric Brain Dipoles”, IEEE Transactions on Biomedical Enginnering, vol.37, no.7, July 1990. F. Vatta, P. Bruno, P.Inchingolo, “MultiregionBicentric-Spheres Models of the Head for the Simulation of Bioelectric Phenomena”, IEEE Transactions on Biomedical Engineering, vol.52, no.3, March 2005. J.Malmivuo, R. Plonsey, “Bioelectromagnetism: Principles and Applications of Bioelectric and Biomagnetic Fields ”, Oxford UnivPr, 1995. R. O. Schmidt, “Multiple emitter location and signal parameter estimation in Proc. RADC Special Estimation Workshop, Griffiths AFB, New York, 1986, pp. 243–258, Reprinted in IEEE Trans. AP-34, pp. 276–280.
[3] [4]
Fig. 9. Reconstructed EEG signals
[5]
[6] [7]
[8] [9] [10]
Fig. 10. Neural power map in x-y system
V. CONCLUSIONS In this work a method to estimate the location, direction and magnitude of the neural brain sources is defined by using scalp electrodes. This method is a combination of two algorithms: a DOA estimation algorithm allows to locate the brain sources, while a Beamforming algorithm allows to extract information about frequencies and amplitudes associated with the estimated source. The importance of the proposed method consists in the fact that it is possible to determine the sources which contribute to the recorded EEG signals. So, sudden events as epileptic seizures can be identified and to locate the area of interest. However improvements are
[11]
[12]
[13] [14] [15] [16]
