Epilepsy & Behavior 31 (2014) 339–345
Contents lists available at ScienceDirect
Epilepsy & Behavior journal homepage: www.elsevier.com/locate/yebeh
Automatic seizure detection using diffusion distance and BLDA in intracranial EEG Shasha Yuan, Weidong Zhou ⁎, Qi Yuan, Yanli Zhang, Qingfang Meng School of Information Science and Engineering, Shandong University, PR China
a r t i c l e
i n f o
Article history: Received 4 May 2013 Revised 3 September 2013 Accepted 3 October 2013 Available online 20 November 2013 Keywords: Seizure detection EEG Discrete wavelet transform (DWT) Diffusion distance Bayesian linear discriminant analysis (BLDA)
a b s t r a c t Approximately 1% of the world's population suffers from epilepsy. An automatic seizure detection system is of great significance in the monitoring and diagnosis of epilepsy. In this paper, a novel method is proposed for automatic seizure detection in intracranial EEG recordings. The EEG recordings are divided into 4-s epochs, and then wavelet decomposition with five scales is performed to the EEG epochs. Detail signals at scales 3, 4, and 5 are selected to form a signal distribution. The diffusion distances are extracted as features, and Bayesian linear discriminant analysis (BLDA) is used as the classifier. A total of 193.75 h of intracranial EEG recordings from 21 patients having 87 seizures are employed to evaluate the system, and the average sensitivity of 94.99%, specificity of 98.74%, and false-detection rate of 0.24/h are achieved. The seizure detection system based on diffusion distance yields a high sensitivity as well as a low false-detection rate for long-term EEG recordings. © 2013 Elsevier Inc. All rights reserved.
1. Introduction Approximately 1% of the world's population suffers from epilepsy at some time in their life [1–3]. Unfortunately, the occurrence of an epileptic seizure is unforeseeable. Electroencephalography (EEG) is a noninvasive method to measure the electrical activity of the brain which contains physiological and pathological information [4,5]. Electroencephalography has become an important tool to comprehend the complex activities of the brain, and long-term EEG recordings are widely used in the monitoring for presurgical evaluation of patients. These EEGs are inspected by neurologists visually to detect epileptic seizures for clinical diagnosis, and this is a very time-consuming and errorprone process because of the huge amount of EEG data [6]. Therefore, developing an automatic seizure detection system is of great importance for speeding up the inspection process and relieving the workload of medical staff. The research on epileptic seizure detection can be traced back to the 1970s, and various methods have been proposed to solve this problem with different degrees of success. One of the earliest seizure detectors was proposed by Gotman in 1982 [7]. It was based on the decomposition of EEG into half-waves and the amplitude of waves relative to the background, the duration and rhythmicity, etc. were extracted for seizure detection. Currently, the techniques of seizure detection
⁎ Corresponding author at: School of Information Science and Engineering, Shandong University, 27 Shanda Road, Jinan, PR China. E-mail address: [email protected] (W. Zhou). 1525-5050/$ – see front matter © 2013 Elsevier Inc. All rights reserved. http://dx.doi.org/10.1016/j.yebeh.2013.10.005
can mainly be divided into five categories: time domain-based, frequency domain-based, time–frequency domain-based, artificial neural network-based, and nonlinear method-based [8]. Murro et al. developed a method that used three features of the time domain and frequency domain including relative amplitude, dominant frequency, and rhythmicity of the EEG for the discriminant analysis [9]. Srinivasan et al. developed a seizure detection system based on a neural network using approximate entropy (ApEn) as the nonlinear feature [10]. The wavelet transform (WT) has become a powerful method in the time–frequency analysis. Its application on feature extraction and property expression has been subjected to an extensive evaluation on EEG signals for seizure detection [6,8,10–12]. Compared with the Fourier transform, the WT can localize abrupt changes in both time and frequency domains and give precise frequency information at low frequencies as well as precise time information at high frequencies. This characteristic makes it suitable for the analysis of nonstationary signals such as the EEG signal. The WT can decompose a signal into different scales with different wavelet coefficients. This property provides flexibility in analyzing signals and extracting features for the classification. Choosing suitable features that can best represent the characteristics of the EEG signals is important for seizure detection in EEGs. Many types of features have been investigated, which include linear features and nonlinear features based on wavelet, spectral, or chaos analysis. These features comprise energy distribution [13], amplitude relative to background activity [9], approximate entropy [8,10], etc. In this study, we propose a novel feature called diffusion distance to quantify the distinction between seizure and nonseizure EEGs. The diffusion distance feature as introduced by Ling and Okada in [14] is a nonlinear distance measurement to measure the dissimilarity between two distributions.
340
S. Yuan et al. / Epilepsy & Behavior 31 (2014) 339–345
It is computed by constructing a Gaussian pyramid from the distribution and summing up the L1 norms of the various levels. Wang Yan et al. proposed a topology-preserved diffusion distance for image retrieval and interest point matching [15]. To our knowledge, the diffusion distance has rarely been applied in the area of seizure detection. A suitable classifier can greatly improve the classification accuracy of a detection system. In previous work, many kinds of classifiers have been used in seizure detection, for instance, support vector machine (SVM), boosting, Fisher linear discriminant analysis (FLDA), and neural network (NN). Bayesian linear discriminant analysis (BLDA) can be regarded as an extension of FLDA and has shown high performance in motor imagery EEG classification [16,17]. Compared with conventional FLDA, BLDA algorithm employs regularization to avoid overfitting to high dimensional and noisy datasets. In this work, we propose a novel method using diffusion distance and BLDA for seizure detection. The remainder of this paper is organized as follows: a brief introduction of the intracranial EEG (iEEG) dataset is given in Section 2. Then in Section 3, the method used in our system is described in detail containing wavelet analysis, feature extraction, and classification. The performance of this experiment is presented in Section 4. Section 5 is a discussion of the results, and a conclusion follows in Section 6. 2. EEG database The long-term EEG database used in this study contains intracranial EEG recordings of 21 patients suffering from medically intractable focal epilepsy, which were recorded during presurgical epilepsy monitoring with invasive electrodes at the Epilepsy Center of the University Hospital of Freiburg, Germany [18]. All EEGs were acquired using a Neurofile NT digital video-EEG system with 128 channels, 256-Hz sampling rate, and a 16-bit analog-to-digital converter. The onset and offset of the seizures for 21 patients were determined by experienced epileptologists based on identification of epileptic patterns preceding clinical manifestation of seizures in EEG recordings, and three in-focus and three extrafocus electrodes were selected for each patient. In our study, we choose three in-focus electrodes for seizure detection. The EEG data used in this study include a total of 193.75h of intracranial EEG recordings and 87 seizures which range from less than 0.2 min to more than 15 min in duration. There are 2–10 h of nonseizure data and 2–5 h of EEG data containing seizure events for each patient. In this work, the seizure data chosen are continuous, and there is 1–3 h of continuous EEG data before and after every seizure as the nonseizure data. The details of this EEG database for each patient are summarized in Table 1.
Table 1 Details of the database we used in this study. The acronyms used in the table are SP: simple partial seizure, CP: complex partial seizure, and GTC: generalized tonic–clonic seizure. Patient
Seizure origin
Seizure type
EEG length (h)
Number of used seizures
Mean seizure duration (s)
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 Total
Temporal Frontal Temporal Temporal, occipital Frontal Temporal Temporal Temporal Frontal Temporal Frontal Frontal Temporal, occipital Temporal Frontal, temporal Temporal Temporal Temporal, occipital Frontal Temporal Temporal –
SP, CP SP, CP, GTC SP, CP SP, CP, GTC SP, CP, GTC CP, GTC SP, CP, GTC SP, CP CP, GTC SP, CP, GTC SP, CP, GTC SP, CP, GTC SP, CP, GTC CP, GTC SP, CP, GTC SP, CP, GTC SP, CP, GTC SP, CP SP, CP, GTC SP, CP, GTC SP, CP –
9 5.62 8.19 10 9.90 6.42 6 3.57 10 10.65 8 8 4 7 10 11.69 14.80 13 13 12.91 12 193.75
4 3 5 5 5 3 3 2 5 5 4 4 2 4 4 5 5 5 4 5 5 87
13.1 118.2 92.7 87.4 44.9 66.9 153.5 163.7 114.7 411.0 157.3 55.1 158.3 216.4 145.4 121.0 86.2 13.7 12.5 85.7 83.1 114.3
3. Methods 3.1. Overview The block diagram of this seizure detection system is exhibited in Fig. 1. On the preprocessing stage, the data of each channel was segmented into 4-second epochs separately without overlap to reduce nonstationarity of EEG signals. Then, we apply discrete wavelet transform (DWT) with five scales on each epoch. For each channel, we randomly choose one nonseizure segment as a reference and compute the diffusion distance between the reference and other EEG epochs as features. The set of features extracted from all channels forms a feature vector to capture spatial information between the channels. The feature vectors are fed to a Bayesian linear discriminant analysis (BLDA) classifier for classification. In the last step, we use a smoothing filter and a collar operation for BLDA outputs to obtain the final detection results. A seizure onset is declared only when two consecutive epochs are classified as seizure epochs. Different steps of the system are described in detail in the following sections. 3.2. Preprocessing
2.1. Training data To train the BLDA classifier, a set of the training data is used. For each patient, one or two seizures are selected randomly as well as the same segments of nonseizure data for training. The diversity of seizures can improve detection performance. In total, 716 segments of seizure and nonseizure data with 1024 points per segment of 21 patients are used as the training data.
The EEG is broken down into epochs using a moving window of size 1024 data points (4 s), which is long enough to capture the main stationary characteristics of the EEG as well as short enough to capture the burst of seizures, for the purpose of feature extraction. Then, the discrete wavelet transform (DWT) is applied to each epoch. The WT decomposes a signal into a set of basic functions, which are called wavelets and obtained from a single prototype wavelet ψ(t) called mother wavelet by dilations and shifting.
2.2. Testing data
1 t−b ψa;b ðt Þ ¼ pffiffiffi ψ a a
A testing dataset is used to test the performance of the proposed algorithm. In this system, 168.63-hour EEG recordings (1.54-h seizure data and 167.09-h nonseizure data) containing 61 seizures of 21 patients are selected as the test data. There are 1387 seizure segments and 105,381 nonseizure segments in total, and the length of each segment is 1024 points. For each patient, the training data and testing data are independent of each other and have no overlap.
ð1Þ
where a is the scaling parameter, and b is the shifting parameter. The continuous wavelet transform (CWT) of a signal f(t) is defined as the convolution between the signal f(t) and the wavelet functions ψa,b(t), Z D E 1 t−b dt: Tf ða; bÞ ¼ f ; ψa;b ¼ pffiffiffi f ðt Þψ a a
ð2Þ
S. Yuan et al. / Epilepsy & Behavior 31 (2014) 339–345
341
S1,1
Channel 1
Wavelet Transform
S1,2
Feature Extraction
X1
Feature Extraction
X2
S1,3
S2,1
Channel 2
Wavelet Transform
S2,2 S2,3
Classifier
Post
Output
Processing
SN,1
Channel N
Wavelet Transform
SN,2
Feature Extraction
XN
SN,3 Fig. 1. Block diagram of the epileptic seizure detection system which consists of preprocessing, feature extraction, classification, and postprocessing.
Furthermore, the scaling factor a and shifting factor b are dispersed as a = 2 j and b = 2 jk, so the amount of computation can be reduced largely and wavelet analysis get more efficient. Mallat developed an effective way for the wavelet decomposition of the signal by passing the signal through a series of low-pass (LP) and high-pass (HP) filter pairs [19]. The original signal is split into different subbands with high-frequency information called the detail coefficient and lowfrequency information called the approximate coefficient. The discrete wavelet transform (DWT) is obtained, and Eq. (2) becomes: D
E
1 Tf ð j; kÞ ¼ f ; ψ j;k ¼ pffiffiffiffiffi 2j
Z
f ðt Þψ
j
!
t−2 k dt: 2j
ð3Þ
It is very essential to select a suitable wavelet and the number of decomposition levels in the application of DWT. Daubechies-4 (db4) wavelet shows better adaptability in matching the change of EEG signals, which have been confirmed in many studies [6,12]. The number of decomposition levels depends on the main domain of frequency to analyze. Even though seizures have a much broader spectrum, the seizures frequently evolve into the frequency components below 30 Hz. Besides, the background noise and artifacts such as muscle activities and power line interference may also exist in the high-frequency band. In this work, we use db-4 wavelet and execute wavelet decomposition with five scales. Electroencephalography epochs are decomposed into five detail coefficients (D1–D5) and an approximation coefficient (A5). Since our data are sampled at 256 Hz, the frequency bands are 64–128 Hz, 32–64 Hz, 16–32 Hz, 8–16 Hz, 4–8 Hz, and 0–4 Hz corresponding to decomposition scales, respectively. In order to capture the major seizure activities and eliminate the high-frequency noise in raw EEG signals, three subbands of scales 3, 4, and 5 are selected, and the EEG features are extracted from these subbands.
in an isolated field satisfies the heat diffusion equation, the calculation formula of T(x, t) is: T ðx; t Þ ¼ T 0 ðxÞ ϕðx; t Þ
ð4Þ
here T 0 ðxÞ ¼ T ðx; 0Þ ¼ dðxÞ ϕðx; t Þ ¼
ð5Þ
2 1 −x e 2t2 : 1=2 ð2πÞ t
ð6Þ
T(x, t) can be regarded as a diffusion process to equalize d1 and d2 through exchanging the distribution value. Hence, a distance between d1 and d2 can be proposed to measure the process as the dissimilarity of two distributions. The distance is defined as: ^ ðd ; d Þ ¼ K 1 2
Z
t 0
kðjT ðx; t ÞjÞdt
ð7Þ
where t is a positive constant, k(.) is a norm that measures how T(x, t) differs from 0, and we use the L1 norm here because it requires a less burdensome calculation. For two m-dimensional distributions d1(x) and d2(x), where x is an ^ ðd1 ; d2 Þ is similar to the one m-dimensional matrix, the definition of K in Eq. (7). The Gaussian pyramid can efficiently describe the continuous diffusion process T(x, t) because aliasing would not be introduced when smoothing and subsampling in the calculation of Gaussian pyramid. Therefore, the diffusion distance based on the Gaussian pyramid is proposed, K ðd1 ; d2 Þ ¼
L X
kðjddl ðxÞjÞ
ð8Þ
l¼0
3.3. Feature extraction
where
A feature which is able to differentiate the nonseizure and seizure data is extremely crucial for seizure detection. In this work, diffusion distance is extracted as the feature for classification. Given two 1D distributions d1(x) and d2(x), the difference between them can be defined as d(x) = d1(x) − d2(x). Assume an isolated temperature field T(x, t), T(x, 0) = d(x) when t = 0. Since the temperature
dd0 ðxÞ ¼ d1 ðxÞ−d2 ðxÞ
ð9Þ
ddl ðxÞ ¼ ½ddl−1 ðxÞ ϕðx; σ Þ ↓2 ðl ¼ 1; …; LÞ ϕðx; t Þ ¼
1 e ð2πÞm=2 t
T −x 2x 2t
:
ð10Þ ð11Þ
342
S. Yuan et al. / Epilepsy & Behavior 31 (2014) 339–345
Diffusion Distance
6
x 10
Seizures Nonseizures
6
Amplitude
Given that likelihood function and prior distribution are Gaussian distribution, the posterior distribution is also Gaussian and can be computed using Bayes' rule. Then, the mean m of the posterior satisfies the following equation: −1 T m ¼ β βXX þ I′ðα Þ Xy: ð16Þ
4 2 0
0
10
20
30
40
50
60
70
80
90
100
EEG epochs Fig. 2. The comparison of the diffusion distance features between the seizure and nonseizure EEG epochs.
L is the number of Gaussian pyramid layers, the sign “↓2” means halfsize downsampling, and σ is the constant standard deviation for the Gaussian filter ϕ. When applying the L1 norm, Eq. (8) is simplified as: K ðd1 ; d2 Þ ¼
L X jddl ðxÞj:
ð12Þ
l¼0
As mentioned above, scales 3, 4, and 5 contain the main information of seizures after wavelet decomposition. Hence, we put the three subsignals together to create a distribution of each epoch, select one epoch randomly as reference, and then compute the diffusion distance with other distributions as the feature fed to the BLDA classifier. Fig. 2 shows the values of diffusion distance extracted from seizure and nonseizure EEG signals. There are 100 seizure epochs and 100 nonseizure epochs respectively chosen randomly from different patients. It can be clearly observed from Fig. 2 that the diffusion distance of a seizure epoch is generally much higher than that of a nonseizure epoch and can be used as features for classification. 3.4. Classification 3.4.1. BLDA classifier Bayesian linear discriminant analysis (BLDA) can be viewed as an extension of Fisher's discriminant analysis (FLDA) which means performing regression in a Bayesian framework [16]. It applies regularization during the training process to prevent overfitting to high dimensional and possibly noisy datasets. No user intervention is required to adjust hyperparameters, and all the model parameters are estimated automatically in the classification process, which reduce time consumption in the training and classification procedure. Assuming the weight vector ω used in Bayesian regression satisfies Gaussian distribution, the likelihood function for ω is: pðDjβ; ωÞ ¼
β 2π
N 2
2 β T exp − X ω−y : 2
ð13Þ
X denotes the matrix that is composed of the training feature vectors, y denotes a vector containing the regression targets, D denotes the pair {X; y}, β denotes the inverse variance of the noise, and N denotes the number of examples in the training set. For the weight vector ω, the expression for its prior distribution is: pðωjα Þ ¼
α D ε 1 1 T 2 2 exp − ω I′ðα Þω 2π 2π 2
ð14Þ
where D is the number of features, and I′(α) is a D + 1 dimensional, diagonal matrix, and ε is a very small value, 2
α 0 60 α I′ðα Þ ¼ 6 4 0 0
3 ⋯ 0 ⋯ 07 7: 5 ⋱ ⋯ ε
ð15Þ
^. y ^ denotes a Thus, we get the linear discriminant equation μ ¼ mT y new input vector for classification. Hence, the linear discriminant result depends on the parameters α and β. According to an iterative algorithm proposed by MacKay [20], the two parameters are optimized to minimize the regression error to obtain the optimal classification discriminant surface. 3.4.2. Postprocessing The postprocessing technique is applied to the outputs of the BLDA classifier to improve the classification accuracy. The postprocessing contains smoothing, applying a threshold, and collar technique. To remove the burrs of the classifier outputs, a central linear moving average filter (MAF) is applied to the BLDA outputs. This filter is a convolution of the inputs with a rectangular pulse of unit area and an effective tool for reducing random noise, as well as keeping the sharpest part [21]. Thus, some false decisions lasting a short time can be modified by this smoothing process. It is defined as: xk ¼
N X 1 ^x 2N þ 1 i¼−N kþi
ð17Þ
where ^x is the input signal, x is the filtered signal, and 2N+1 denotes the smoothing length of MAF. The number of N is patient-specific, and we adjust it on the basis of different patients within the range of no greater than 35 in this work. Then, the smoothed outputs are compared to a threshold. After comparison, we obtain the binary decisions: 1 — nonseizure; 0 — seizure. The value of threshold is different for each patient in this study. However, seizure onset is changing gradually, and the smoothing procedure may make the beginning and ending of seizures obscure. The differences of the features between a seizure and nonseizure are not very clear on both sides of seizures. Hence, a collar operation is used for seizure decision in the last step to compensate for the missing part of a seizure. Each binary decision is extended in x epochs from both sides in this process. A seizure was defined as a paroxysmal burst clearly different from the background and lasting at least 6s [22], so at least two consecutive epochs marked 0 are declared as a seizure event. It helps to avoid false detections due to very brief seizurelike EEG morphologies that were actually non-pathologic. 4. Results In this work, all the experiments were carried out in a MATLAB 7.0 environment running in an AMD Sampson processor with 2.40 GHz. We implement wavelet decomposition to EEG epochs with five scales and choose scales 3, 4, and 5 to extract diffusion distance features using the method described in Section 3.3. The training data and testing date are selected randomly from iEEG recordings for each patient. Afterwards, the features are sent into the BLDA classifier for classification. Finally, the postprocessing is performed. We evaluated the performance of this proposed method with two approaches: the segment-based level and the event-based level. For the segment-based level, the EEG segments classified by our algorithm are compared with those marked by experts. The overall performance of our system on segment-based level is seen in Table 2. There are three statistical measures: sensitivity, specificity, and recognition accuracy, which are defined as follows: • Sensitivity: true positive / the total number of seizures marked by the EEG experts. The number of seizures identified both by our detector and by the EEG experts is defined as true positive (TP).
S. Yuan et al. / Epilepsy & Behavior 31 (2014) 339–345 Table 2 The results of our detection method for each patient from the evaluation of the performance on a segment-based level. Patient
Sensitivity (%)
Specificity (%)
Recognition accuracy (%)
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 Mean
62.5 100 100 100 100 100 100 95.46 96 89.89 100 100 100 98.01 88.06 95.77 100 100 70 99 100 94.99
99.98 99.15 99.49 98.50 96.95 99.94 99.91 96.30 97.57 92.10 99.91 99.79 99.77 99.31 99.15 99.36 99.85 99.79 99.98 97.42 99.23 98.74
99.94 99.16 99.50 98.52 96.97 99.94 99.91 96.29 97.56 92.08 99.91 99.79 99.78 99.26 99.06 99.33 99.85 99.79 99.94 97.44 99.23 98.73
• Specificity: true negative / the total number of nonseizures marked by the EEG experts. True negative (TN) represents the nonseizures identified by our method and by the experts. • Recognition accuracy: number of correctly classified EEG epochs / the total number of EEG epochs. All results for each patient according to the three statistical measures are depicted on Table 2. It can be observed that the best sensitivity of 100%, specificity of 99.98%, and recognition accuracy of 99.94% were achieved for different patients. In addition, the mean of sensitivity for all the 21 patients was nearly 95.0%, and the means of the specificity and recognition accuracy are both over 98.5%. Furthermore, twelve patients (patients 2, 3, 4, 5, 6, 7, 11, 12, 13, 17, 18, and 21), over half of all the patients, had a sensitivity of 100%. The specificities of all patients are greater than 92%, and among them, there were fifteen patients having specificities higher than 99%. For patient 1, the sensitivity was low because the duration of seizure activities was very short, so some seizure epochs have not been identified effectively. Moreover, considering the feasibility of the proposed method, the event-based approach is also employed to evaluate our algorithm. To assess the system at this level, two statistical measures are calculated: the number of true detected seizures and the false-detection rate. Any detected event by our system that overlapped a seizure event labeled by the EEG experts is defined as a true seizure detection. A sequence of consecutive false positives, which are not adjacent with the seizures, is defined as a false detection. False positive (FP) means the number of seizure segments identified only by our algorithm but not by experts. The system performance on the event-based level is shown in Table 3. In our study, 61 seizures were used to test our method totally, and 58 seizures were detected correctly, and 100% of seizure events were detected for 18 patients. Fig. 3 illustrates an example of a detected seizure which had obvious changes in the used feature. Our method only missed 3 seizures, which come from patients 1, 15, and 19. For patients 1 and 19, the duration of missed seizures is shorter than 12s. The missed seizure of patient 15 is presented in Fig. 4. As shown in this figure, the seizure is very short, and there is no significant difference in the frequency and amplitude between the seizure and nonseizure EEG parts. This seizure is missed because of the unobvious change of diffusion distance features. In total, the average false-detection rate is 0.24/h for all the 21 patients, there were no false detections in the records of approximately half of all patients.
343
Table 3 The results of our detection method for each patient from the evaluation of the performance on an event-based level. Patient
Number of seizures marked by experts
Number of true seizures marked by our system
False detections/h
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 Total
2 2 4 4 3 2 2 1 4 4 2 3 1 3 3 3 4 4 2 4 4 61
1 2 4 4 3 2 2 1 4 4 2 3 1 3 2 3 4 4 1 4 4 58
0 0.22 0 0.11 1.14 0 0 1.17 0.56 0.69 0 0 0 0 0.11 0.10 0 0.25 0 0.50 0.27 0.24
5. Discussion Seizure detection is a difficult task because of the vagaries of ictal EEGs. A perfect detection tool, which assists the clinical evaluation of patients with epilepsy, should be able to reduce the amount of inspected data and have high sensitivity, low false-detection rate. An efficient feature will play a significant role in the entire system. The diffusion distance applied in our study shows its powerful ability to distinguish seizure and nonseizure signals. The significant advantage of this feature is that its computational complexity is much less when compared with other distances, such as the earth mover's distance [14,23]. Diffusion distance is robust to the noise in comparison to the EMD. Furthermore, this feature describes the distance between two different distributions.
Fig. 3. An example of a detected seizure of patient 6. (a) One-hour raw EEG data. (b) The raw EEG data from 47 min to 50.5 min containing one seizure. (c) The diffusion distance features extracted from the one-hour EEG data. The seizure event marked by the EEG experts is between the two vertical lines.
344
S. Yuan et al. / Epilepsy & Behavior 31 (2014) 339–345 Table 4 Comparison of the performance for different methods.
Fig. 4. An example of a missed seizure of patient 15. (a) One-hour raw EEG data. (b) The raw EEG data from 56 min to 57.5 min containing one seizure. (c) The diffusion distance features extracted from the one-hour EEG data. The seizure event marked by the EEG experts is between the two vertical lines.
For every epoch of EEG, it could be transformed into a distribution, and we could calculate the diffusion distance of nonseizure and nonseizure epochs as well as the distance of seizure and nonseizure epochs. It is noticed from Fig. 2 that the diffusion distance of the different categories is larger than that of the same categories, and the diffusion distance can be used to discriminate seizure and nonseizure epochs. Fig. 5 shows the mean value and standard deviation of the diffusion distance. This statistical analysis is achieved by SPSS software. We randomly choose 100 segments from seizure data and nonseizure data, respectively. The figure indicates that the mean value of this feature of seizure signals is much higher than that of nonseizure signals. Results of the independent two-sample T-test indicate that the distinctions in diffusion distance between nonseizure and seizure EEGs are statistically significant (p-values b 0.0001). In addition, the postprocessing part enables us to deal with imprecise boundaries between nonseizure and seizure EEGs and locate the seizure more accurately. These processing technologies all make a contribution to improving the precision of classification. In view of the application of wavelet decomposition to eliminate high-frequency components, it may miss the seizure onset
Method
Sensitivity (%)
Number of seizures used
Number of channels used
Multistage seizure detection [24] A fuzzy logic algorithm [25] Differential windowed variance (DWV) algorithm [26] Improved patient-specific seizure detection [27] Our proposed algorithm
87.5 95.8 91.525
24 56 59
3 6 3
78
63
–
94.99
61 (test) + 26 (training)
3
and hence affect the delay between seizure onset and detection. A few short seizures with predominantly high-frequency discharges may also be missed for this reason. However, the collar technique in the postprocessing part compensates for the missing part of a seizure to make the seizure detection more accurate. Moreover, although keeping the high-frequency components is helpful for detecting subclinical electrographic events and short seizures with predominantly high-frequency discharges, the existence of high-frequency artifacts, such as muscle activities and the power line interferences, could result in the obvious increase of the false-detection rate. So this work employs wavelet decomposition to suppress those high-frequency artifacts and achieve a low false-detection rate. The EEG database used in this study has been applied in other seizure detection systems [24–27]. Raghunathan proposed a multistage seizure detection which used a discrete wavelet transform-based filtering block and a “feature extraction block” containing coastline and variance energy features. The average sensitivity was 87.5% from 5 patients with 24 seizures [24]. Majumdar and Vardhan proposed a differential windowed variance (DWV) algorithm to automatically detect seizure onsets in continuous ECoG (depth-EEG) signals [26]. The sensitivity of this method was 91.525% with 59h of seizure from 15 patients with epilepsy. In both of their methods, the channels they used were the same as the ones we used in this study. Chua developed an improved patient-specific seizure detection by applying a subject-independent quadratic discriminant classifier which was tested on 63 seizures of 15 patients and obtained a sensitivity of 78% [27]. The number of channels was not clearly mentioned in Chua's study. Table 4 depicts a comparison on the results between their methods and our method. Compared to other automatic seizure-detection systems, our proposed algorithm yields a higher sensitivity. Rabbi and Fazel Rezai presented a multistage fuzzy logic algorithm for epileptic seizure onset detection using features of amplitude, frequency, and entropy. The method was evaluated on the same iEEG database from 20 patients with 56 seizures. A sensitivity of 95.8% and a false-detection rate of 0.26/h were achieved [25]. Out of 56 seizures analyzed, their system detected 54 seizures, and 2 seizures were missed. The data from patient 10 (of the FSPEEG database) were discarded because of the presence of electrode movement artifacts. In addition, they employed all the six channels, but we used only three. Compared to this system, our method which evaluated 61 seizures from all of the 21 patients yields a comparable sensitivity.
6. Conclusion
Fig. 5. The mean and standard deviation of the diffusion distance feature extracted from EEG signals.
In this work, we presented a novel method of seizure detection using diffusion distance and BLDA on intracranial EEGs. The diffusion distance, as a novel and effective feature, was extracted to distinguish seizure and nonseizure EEG epochs. The Bayesian linear discriminant analysis was used to classify seizure and nonseizure epochs. The algorithm has been validated on an intracranial EEG database, and the experimental results show that this proposed method achieves a promising performance with a high sensitivity of 94.99% and a low false-detection rate
S. Yuan et al. / Epilepsy & Behavior 31 (2014) 339–345
of 0.24/h. Our study indicates the potential clinical application of the proposed method. In future work, we will focus on integrating more available features to improve the performance and robustness of this seizure detector. Acknowledgements The authors Shasha Yuan, Weidong Zhou, Qi Yuan, Yanli Zhang, and Qingfang Meng are also with the Suzhou Institute of Shandong University, Suzhou, PR China. The support of the Program of Science and Technology of Suzhou (No. ZXY2013030) and the Independent Innovation Foundation of Shandong University (No. 2012DX008) is gratefully acknowledged. References [1] Hauser WA, Annegers JF, Rocca WA. Descriptive epidemiology of epilepsy: contributions of population-based studies from Rochester, Minnesota. Mayo Clin Proc 1996;71:578–86. [2] Stam CJ, Pijn JP, Suffczynski P, Lopezda Silva FH. Dynamics of the human alpha rhythm: evidence for non-linearity? Clin Neurophysiol 1999;110(10):1801–13. [3] Iasemidis LD, Shiau DS, Chaovalitwongse W, Sackellares JC, Pardalos PM, Principe JC, et al. Adaptive epileptic seizure prediction system. IEEE Trans Biomed Eng 2003;50(5):616–27. [4] Shoeb A, Edwards H, Connolly J, Bourgeois B, Treves ST, Guttag J. Patient-specific seizure onset detection. Epilepsy Behav 2004;5(4):483–98. [5] Aarabi A, Fazel-Rezai R, Aghakhani Y. A fuzzy rule-based system for epileptic seizure detection in intracranial EEG. Clin Neurophysiol 2009;120(9):1648–57. [6] Khan YU, Gotman J. Wavelet based automatic seizure detection in intracerebral electroencephalogram. Clin Neurophysiol 2003;114(5):898–908. [7] Gotman J. Automatic recognition of epileptic seizures in the EEG. Electroencephalogr Clin Neurophysiol 1982;54(5):530–40. [8] Ocak H. Automatic detection of epileptic seizures in EEG using discrete wavelet transform and approximate entropy. Expert Syst Appl 2009;36:2027–36. [9] Murro AM, King DW, Smith JR, Gallagher BB, Flanigin HF, Meador K. Computerized seizure detection of complex partial seizures. Electroencephalogr Clin Neurophysiol 1991;79:330–3.
345
[10] Srinivasan V, Eswaran C, Sriraam NH. Approximate entropy-based epileptic EEG detection using artificial neural networks. IEEE Trans Inf Technol Biomed 2007;11(3): 288–95. [11] Khan YU, Rafiuddin N, Farooq O. Automated seizure detection in scalp EEG using multiple wavelet scales. IEEE International Conference on Signal Processing, Computing and Control (ISPCC); 2012. p. 1–5. [12] Subasi A. Epileptic seizure detection using dynamic wavelet network. Expert Syst Appl 2005;29:343–55. [13] Qu H, Gotman J. A patient-specific algorithm for the detection of seizure onset in long-term EEG monitoring: possible use as a warning device. IEEE Trans Biomed Eng 1997;44(2):115–22. [14] Ling H, Okada K. Diffusion distance for histogram comparison. Proc IEEE Comput Soc Conf Comput Vis Pattern Recognit 2006;1:246–53. [15] Yan W, Wang Q, Liu Q, Lu H, Ma S. Topology-preserved diffusion distance for histogram comparison. Proceedings of the British Machine Conference; 2007. p. 61.1–61.10. [16] Hoffmann U, Vesin JM, Ebrahimi T, Diserens K. An efficient P300-based brain– computer interface for disabled subjects. J Neurosci Methods 2008;167(1): 115–25. [17] Lei X, Yang P, Yao D. An empirical Bayesian framework for brain–computer interfaces. IEEE Trans Neural Syst Rehabil Eng 2009;17(6):521–9. [18] Freiburg seizure prediction project. Freiburg, Germany 2008. [On-line]. Available: http://epilepsy.uni-freburg.de/freiburg-seizure-prediction-project/eeg-database. [19] Mallat S. A theory for multiresolution signal decomposition: the wavelet representation. IEEE Trans Pattern Anal Mach Intell 1989;11(7):674–93. [20] MacKay DJC. Bayesian interpolation. Neural Comput 1992;4(3):415–47. [21] Temko A, Thomas E, Marnane W, Lightbody G, Boylan G. EEG-based neonatal seizure detection with support vector machines. Clin Neurophysiol 2011;122: 464–73. [22] Grewal S, Gotman J. An automatic warning system for epileptic seizures recorded on intracerebral EEGs. Clin Neurophysiol 2005;116(10):2460–72. [23] Rubner Y, Tomasi C, Guibas LJ. The earth mover's distance as a metric for image retrieval. Int J Comput Vis 2000;40(2):99–121. [24] Raghunathan S, Jaitli A, Irazoqui PP. Multistage seizure detection techniques optimized for low-power hardware platforms. Epilepsy Behav 2011;22:61–8. [25] Rabbi AF, Fazel-Rezai R. A fuzzy logic system for seizure onset detection in intracranial EEG. Comput Intell Neurosci 2012;2012:705140. [26] Majumdar K, Vardhan P. Automatic seizure detection in ECoG by differential operator and windowed variance. IEEE Trans Neural Syst Rehabil Eng 2011;19(4): 356–65. [27] Chua EC, Patel K, Fitzsimons M, Bleakley CJ. Improved patient specific seizure detection during pre-surgical evaluation. Clin Neurophysiol 2011;122:672–9.
Contents lists available at ScienceDirect
Epilepsy & Behavior journal homepage: www.elsevier.com/locate/yebeh
Automatic seizure detection using diffusion distance and BLDA in intracranial EEG Shasha Yuan, Weidong Zhou ⁎, Qi Yuan, Yanli Zhang, Qingfang Meng School of Information Science and Engineering, Shandong University, PR China
a r t i c l e
i n f o
Article history: Received 4 May 2013 Revised 3 September 2013 Accepted 3 October 2013 Available online 20 November 2013 Keywords: Seizure detection EEG Discrete wavelet transform (DWT) Diffusion distance Bayesian linear discriminant analysis (BLDA)
a b s t r a c t Approximately 1% of the world's population suffers from epilepsy. An automatic seizure detection system is of great significance in the monitoring and diagnosis of epilepsy. In this paper, a novel method is proposed for automatic seizure detection in intracranial EEG recordings. The EEG recordings are divided into 4-s epochs, and then wavelet decomposition with five scales is performed to the EEG epochs. Detail signals at scales 3, 4, and 5 are selected to form a signal distribution. The diffusion distances are extracted as features, and Bayesian linear discriminant analysis (BLDA) is used as the classifier. A total of 193.75 h of intracranial EEG recordings from 21 patients having 87 seizures are employed to evaluate the system, and the average sensitivity of 94.99%, specificity of 98.74%, and false-detection rate of 0.24/h are achieved. The seizure detection system based on diffusion distance yields a high sensitivity as well as a low false-detection rate for long-term EEG recordings. © 2013 Elsevier Inc. All rights reserved.
1. Introduction Approximately 1% of the world's population suffers from epilepsy at some time in their life [1–3]. Unfortunately, the occurrence of an epileptic seizure is unforeseeable. Electroencephalography (EEG) is a noninvasive method to measure the electrical activity of the brain which contains physiological and pathological information [4,5]. Electroencephalography has become an important tool to comprehend the complex activities of the brain, and long-term EEG recordings are widely used in the monitoring for presurgical evaluation of patients. These EEGs are inspected by neurologists visually to detect epileptic seizures for clinical diagnosis, and this is a very time-consuming and errorprone process because of the huge amount of EEG data [6]. Therefore, developing an automatic seizure detection system is of great importance for speeding up the inspection process and relieving the workload of medical staff. The research on epileptic seizure detection can be traced back to the 1970s, and various methods have been proposed to solve this problem with different degrees of success. One of the earliest seizure detectors was proposed by Gotman in 1982 [7]. It was based on the decomposition of EEG into half-waves and the amplitude of waves relative to the background, the duration and rhythmicity, etc. were extracted for seizure detection. Currently, the techniques of seizure detection
⁎ Corresponding author at: School of Information Science and Engineering, Shandong University, 27 Shanda Road, Jinan, PR China. E-mail address: [email protected] (W. Zhou). 1525-5050/$ – see front matter © 2013 Elsevier Inc. All rights reserved. http://dx.doi.org/10.1016/j.yebeh.2013.10.005
can mainly be divided into five categories: time domain-based, frequency domain-based, time–frequency domain-based, artificial neural network-based, and nonlinear method-based [8]. Murro et al. developed a method that used three features of the time domain and frequency domain including relative amplitude, dominant frequency, and rhythmicity of the EEG for the discriminant analysis [9]. Srinivasan et al. developed a seizure detection system based on a neural network using approximate entropy (ApEn) as the nonlinear feature [10]. The wavelet transform (WT) has become a powerful method in the time–frequency analysis. Its application on feature extraction and property expression has been subjected to an extensive evaluation on EEG signals for seizure detection [6,8,10–12]. Compared with the Fourier transform, the WT can localize abrupt changes in both time and frequency domains and give precise frequency information at low frequencies as well as precise time information at high frequencies. This characteristic makes it suitable for the analysis of nonstationary signals such as the EEG signal. The WT can decompose a signal into different scales with different wavelet coefficients. This property provides flexibility in analyzing signals and extracting features for the classification. Choosing suitable features that can best represent the characteristics of the EEG signals is important for seizure detection in EEGs. Many types of features have been investigated, which include linear features and nonlinear features based on wavelet, spectral, or chaos analysis. These features comprise energy distribution [13], amplitude relative to background activity [9], approximate entropy [8,10], etc. In this study, we propose a novel feature called diffusion distance to quantify the distinction between seizure and nonseizure EEGs. The diffusion distance feature as introduced by Ling and Okada in [14] is a nonlinear distance measurement to measure the dissimilarity between two distributions.
340
S. Yuan et al. / Epilepsy & Behavior 31 (2014) 339–345
It is computed by constructing a Gaussian pyramid from the distribution and summing up the L1 norms of the various levels. Wang Yan et al. proposed a topology-preserved diffusion distance for image retrieval and interest point matching [15]. To our knowledge, the diffusion distance has rarely been applied in the area of seizure detection. A suitable classifier can greatly improve the classification accuracy of a detection system. In previous work, many kinds of classifiers have been used in seizure detection, for instance, support vector machine (SVM), boosting, Fisher linear discriminant analysis (FLDA), and neural network (NN). Bayesian linear discriminant analysis (BLDA) can be regarded as an extension of FLDA and has shown high performance in motor imagery EEG classification [16,17]. Compared with conventional FLDA, BLDA algorithm employs regularization to avoid overfitting to high dimensional and noisy datasets. In this work, we propose a novel method using diffusion distance and BLDA for seizure detection. The remainder of this paper is organized as follows: a brief introduction of the intracranial EEG (iEEG) dataset is given in Section 2. Then in Section 3, the method used in our system is described in detail containing wavelet analysis, feature extraction, and classification. The performance of this experiment is presented in Section 4. Section 5 is a discussion of the results, and a conclusion follows in Section 6. 2. EEG database The long-term EEG database used in this study contains intracranial EEG recordings of 21 patients suffering from medically intractable focal epilepsy, which were recorded during presurgical epilepsy monitoring with invasive electrodes at the Epilepsy Center of the University Hospital of Freiburg, Germany [18]. All EEGs were acquired using a Neurofile NT digital video-EEG system with 128 channels, 256-Hz sampling rate, and a 16-bit analog-to-digital converter. The onset and offset of the seizures for 21 patients were determined by experienced epileptologists based on identification of epileptic patterns preceding clinical manifestation of seizures in EEG recordings, and three in-focus and three extrafocus electrodes were selected for each patient. In our study, we choose three in-focus electrodes for seizure detection. The EEG data used in this study include a total of 193.75h of intracranial EEG recordings and 87 seizures which range from less than 0.2 min to more than 15 min in duration. There are 2–10 h of nonseizure data and 2–5 h of EEG data containing seizure events for each patient. In this work, the seizure data chosen are continuous, and there is 1–3 h of continuous EEG data before and after every seizure as the nonseizure data. The details of this EEG database for each patient are summarized in Table 1.
Table 1 Details of the database we used in this study. The acronyms used in the table are SP: simple partial seizure, CP: complex partial seizure, and GTC: generalized tonic–clonic seizure. Patient
Seizure origin
Seizure type
EEG length (h)
Number of used seizures
Mean seizure duration (s)
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 Total
Temporal Frontal Temporal Temporal, occipital Frontal Temporal Temporal Temporal Frontal Temporal Frontal Frontal Temporal, occipital Temporal Frontal, temporal Temporal Temporal Temporal, occipital Frontal Temporal Temporal –
SP, CP SP, CP, GTC SP, CP SP, CP, GTC SP, CP, GTC CP, GTC SP, CP, GTC SP, CP CP, GTC SP, CP, GTC SP, CP, GTC SP, CP, GTC SP, CP, GTC CP, GTC SP, CP, GTC SP, CP, GTC SP, CP, GTC SP, CP SP, CP, GTC SP, CP, GTC SP, CP –
9 5.62 8.19 10 9.90 6.42 6 3.57 10 10.65 8 8 4 7 10 11.69 14.80 13 13 12.91 12 193.75
4 3 5 5 5 3 3 2 5 5 4 4 2 4 4 5 5 5 4 5 5 87
13.1 118.2 92.7 87.4 44.9 66.9 153.5 163.7 114.7 411.0 157.3 55.1 158.3 216.4 145.4 121.0 86.2 13.7 12.5 85.7 83.1 114.3
3. Methods 3.1. Overview The block diagram of this seizure detection system is exhibited in Fig. 1. On the preprocessing stage, the data of each channel was segmented into 4-second epochs separately without overlap to reduce nonstationarity of EEG signals. Then, we apply discrete wavelet transform (DWT) with five scales on each epoch. For each channel, we randomly choose one nonseizure segment as a reference and compute the diffusion distance between the reference and other EEG epochs as features. The set of features extracted from all channels forms a feature vector to capture spatial information between the channels. The feature vectors are fed to a Bayesian linear discriminant analysis (BLDA) classifier for classification. In the last step, we use a smoothing filter and a collar operation for BLDA outputs to obtain the final detection results. A seizure onset is declared only when two consecutive epochs are classified as seizure epochs. Different steps of the system are described in detail in the following sections. 3.2. Preprocessing
2.1. Training data To train the BLDA classifier, a set of the training data is used. For each patient, one or two seizures are selected randomly as well as the same segments of nonseizure data for training. The diversity of seizures can improve detection performance. In total, 716 segments of seizure and nonseizure data with 1024 points per segment of 21 patients are used as the training data.
The EEG is broken down into epochs using a moving window of size 1024 data points (4 s), which is long enough to capture the main stationary characteristics of the EEG as well as short enough to capture the burst of seizures, for the purpose of feature extraction. Then, the discrete wavelet transform (DWT) is applied to each epoch. The WT decomposes a signal into a set of basic functions, which are called wavelets and obtained from a single prototype wavelet ψ(t) called mother wavelet by dilations and shifting.
2.2. Testing data
1 t−b ψa;b ðt Þ ¼ pffiffiffi ψ a a
A testing dataset is used to test the performance of the proposed algorithm. In this system, 168.63-hour EEG recordings (1.54-h seizure data and 167.09-h nonseizure data) containing 61 seizures of 21 patients are selected as the test data. There are 1387 seizure segments and 105,381 nonseizure segments in total, and the length of each segment is 1024 points. For each patient, the training data and testing data are independent of each other and have no overlap.
ð1Þ
where a is the scaling parameter, and b is the shifting parameter. The continuous wavelet transform (CWT) of a signal f(t) is defined as the convolution between the signal f(t) and the wavelet functions ψa,b(t), Z D E 1 t−b dt: Tf ða; bÞ ¼ f ; ψa;b ¼ pffiffiffi f ðt Þψ a a
ð2Þ
S. Yuan et al. / Epilepsy & Behavior 31 (2014) 339–345
341
S1,1
Channel 1
Wavelet Transform
S1,2
Feature Extraction
X1
Feature Extraction
X2
S1,3
S2,1
Channel 2
Wavelet Transform
S2,2 S2,3
Classifier
Post
Output
Processing
SN,1
Channel N
Wavelet Transform
SN,2
Feature Extraction
XN
SN,3 Fig. 1. Block diagram of the epileptic seizure detection system which consists of preprocessing, feature extraction, classification, and postprocessing.
Furthermore, the scaling factor a and shifting factor b are dispersed as a = 2 j and b = 2 jk, so the amount of computation can be reduced largely and wavelet analysis get more efficient. Mallat developed an effective way for the wavelet decomposition of the signal by passing the signal through a series of low-pass (LP) and high-pass (HP) filter pairs [19]. The original signal is split into different subbands with high-frequency information called the detail coefficient and lowfrequency information called the approximate coefficient. The discrete wavelet transform (DWT) is obtained, and Eq. (2) becomes: D
E
1 Tf ð j; kÞ ¼ f ; ψ j;k ¼ pffiffiffiffiffi 2j
Z
f ðt Þψ
j
!
t−2 k dt: 2j
ð3Þ
It is very essential to select a suitable wavelet and the number of decomposition levels in the application of DWT. Daubechies-4 (db4) wavelet shows better adaptability in matching the change of EEG signals, which have been confirmed in many studies [6,12]. The number of decomposition levels depends on the main domain of frequency to analyze. Even though seizures have a much broader spectrum, the seizures frequently evolve into the frequency components below 30 Hz. Besides, the background noise and artifacts such as muscle activities and power line interference may also exist in the high-frequency band. In this work, we use db-4 wavelet and execute wavelet decomposition with five scales. Electroencephalography epochs are decomposed into five detail coefficients (D1–D5) and an approximation coefficient (A5). Since our data are sampled at 256 Hz, the frequency bands are 64–128 Hz, 32–64 Hz, 16–32 Hz, 8–16 Hz, 4–8 Hz, and 0–4 Hz corresponding to decomposition scales, respectively. In order to capture the major seizure activities and eliminate the high-frequency noise in raw EEG signals, three subbands of scales 3, 4, and 5 are selected, and the EEG features are extracted from these subbands.
in an isolated field satisfies the heat diffusion equation, the calculation formula of T(x, t) is: T ðx; t Þ ¼ T 0 ðxÞ ϕðx; t Þ
ð4Þ
here T 0 ðxÞ ¼ T ðx; 0Þ ¼ dðxÞ ϕðx; t Þ ¼
ð5Þ
2 1 −x e 2t2 : 1=2 ð2πÞ t
ð6Þ
T(x, t) can be regarded as a diffusion process to equalize d1 and d2 through exchanging the distribution value. Hence, a distance between d1 and d2 can be proposed to measure the process as the dissimilarity of two distributions. The distance is defined as: ^ ðd ; d Þ ¼ K 1 2
Z
t 0
kðjT ðx; t ÞjÞdt
ð7Þ
where t is a positive constant, k(.) is a norm that measures how T(x, t) differs from 0, and we use the L1 norm here because it requires a less burdensome calculation. For two m-dimensional distributions d1(x) and d2(x), where x is an ^ ðd1 ; d2 Þ is similar to the one m-dimensional matrix, the definition of K in Eq. (7). The Gaussian pyramid can efficiently describe the continuous diffusion process T(x, t) because aliasing would not be introduced when smoothing and subsampling in the calculation of Gaussian pyramid. Therefore, the diffusion distance based on the Gaussian pyramid is proposed, K ðd1 ; d2 Þ ¼
L X
kðjddl ðxÞjÞ
ð8Þ
l¼0
3.3. Feature extraction
where
A feature which is able to differentiate the nonseizure and seizure data is extremely crucial for seizure detection. In this work, diffusion distance is extracted as the feature for classification. Given two 1D distributions d1(x) and d2(x), the difference between them can be defined as d(x) = d1(x) − d2(x). Assume an isolated temperature field T(x, t), T(x, 0) = d(x) when t = 0. Since the temperature
dd0 ðxÞ ¼ d1 ðxÞ−d2 ðxÞ
ð9Þ
ddl ðxÞ ¼ ½ddl−1 ðxÞ ϕðx; σ Þ ↓2 ðl ¼ 1; …; LÞ ϕðx; t Þ ¼
1 e ð2πÞm=2 t
T −x 2x 2t
:
ð10Þ ð11Þ
342
S. Yuan et al. / Epilepsy & Behavior 31 (2014) 339–345
Diffusion Distance
6
x 10
Seizures Nonseizures
6
Amplitude
Given that likelihood function and prior distribution are Gaussian distribution, the posterior distribution is also Gaussian and can be computed using Bayes' rule. Then, the mean m of the posterior satisfies the following equation: −1 T m ¼ β βXX þ I′ðα Þ Xy: ð16Þ
4 2 0
0
10
20
30
40
50
60
70
80
90
100
EEG epochs Fig. 2. The comparison of the diffusion distance features between the seizure and nonseizure EEG epochs.
L is the number of Gaussian pyramid layers, the sign “↓2” means halfsize downsampling, and σ is the constant standard deviation for the Gaussian filter ϕ. When applying the L1 norm, Eq. (8) is simplified as: K ðd1 ; d2 Þ ¼
L X jddl ðxÞj:
ð12Þ
l¼0
As mentioned above, scales 3, 4, and 5 contain the main information of seizures after wavelet decomposition. Hence, we put the three subsignals together to create a distribution of each epoch, select one epoch randomly as reference, and then compute the diffusion distance with other distributions as the feature fed to the BLDA classifier. Fig. 2 shows the values of diffusion distance extracted from seizure and nonseizure EEG signals. There are 100 seizure epochs and 100 nonseizure epochs respectively chosen randomly from different patients. It can be clearly observed from Fig. 2 that the diffusion distance of a seizure epoch is generally much higher than that of a nonseizure epoch and can be used as features for classification. 3.4. Classification 3.4.1. BLDA classifier Bayesian linear discriminant analysis (BLDA) can be viewed as an extension of Fisher's discriminant analysis (FLDA) which means performing regression in a Bayesian framework [16]. It applies regularization during the training process to prevent overfitting to high dimensional and possibly noisy datasets. No user intervention is required to adjust hyperparameters, and all the model parameters are estimated automatically in the classification process, which reduce time consumption in the training and classification procedure. Assuming the weight vector ω used in Bayesian regression satisfies Gaussian distribution, the likelihood function for ω is: pðDjβ; ωÞ ¼
β 2π
N 2
2 β T exp − X ω−y : 2
ð13Þ
X denotes the matrix that is composed of the training feature vectors, y denotes a vector containing the regression targets, D denotes the pair {X; y}, β denotes the inverse variance of the noise, and N denotes the number of examples in the training set. For the weight vector ω, the expression for its prior distribution is: pðωjα Þ ¼
α D ε 1 1 T 2 2 exp − ω I′ðα Þω 2π 2π 2
ð14Þ
where D is the number of features, and I′(α) is a D + 1 dimensional, diagonal matrix, and ε is a very small value, 2
α 0 60 α I′ðα Þ ¼ 6 4 0 0
3 ⋯ 0 ⋯ 07 7: 5 ⋱ ⋯ ε
ð15Þ
^. y ^ denotes a Thus, we get the linear discriminant equation μ ¼ mT y new input vector for classification. Hence, the linear discriminant result depends on the parameters α and β. According to an iterative algorithm proposed by MacKay [20], the two parameters are optimized to minimize the regression error to obtain the optimal classification discriminant surface. 3.4.2. Postprocessing The postprocessing technique is applied to the outputs of the BLDA classifier to improve the classification accuracy. The postprocessing contains smoothing, applying a threshold, and collar technique. To remove the burrs of the classifier outputs, a central linear moving average filter (MAF) is applied to the BLDA outputs. This filter is a convolution of the inputs with a rectangular pulse of unit area and an effective tool for reducing random noise, as well as keeping the sharpest part [21]. Thus, some false decisions lasting a short time can be modified by this smoothing process. It is defined as: xk ¼
N X 1 ^x 2N þ 1 i¼−N kþi
ð17Þ
where ^x is the input signal, x is the filtered signal, and 2N+1 denotes the smoothing length of MAF. The number of N is patient-specific, and we adjust it on the basis of different patients within the range of no greater than 35 in this work. Then, the smoothed outputs are compared to a threshold. After comparison, we obtain the binary decisions: 1 — nonseizure; 0 — seizure. The value of threshold is different for each patient in this study. However, seizure onset is changing gradually, and the smoothing procedure may make the beginning and ending of seizures obscure. The differences of the features between a seizure and nonseizure are not very clear on both sides of seizures. Hence, a collar operation is used for seizure decision in the last step to compensate for the missing part of a seizure. Each binary decision is extended in x epochs from both sides in this process. A seizure was defined as a paroxysmal burst clearly different from the background and lasting at least 6s [22], so at least two consecutive epochs marked 0 are declared as a seizure event. It helps to avoid false detections due to very brief seizurelike EEG morphologies that were actually non-pathologic. 4. Results In this work, all the experiments were carried out in a MATLAB 7.0 environment running in an AMD Sampson processor with 2.40 GHz. We implement wavelet decomposition to EEG epochs with five scales and choose scales 3, 4, and 5 to extract diffusion distance features using the method described in Section 3.3. The training data and testing date are selected randomly from iEEG recordings for each patient. Afterwards, the features are sent into the BLDA classifier for classification. Finally, the postprocessing is performed. We evaluated the performance of this proposed method with two approaches: the segment-based level and the event-based level. For the segment-based level, the EEG segments classified by our algorithm are compared with those marked by experts. The overall performance of our system on segment-based level is seen in Table 2. There are three statistical measures: sensitivity, specificity, and recognition accuracy, which are defined as follows: • Sensitivity: true positive / the total number of seizures marked by the EEG experts. The number of seizures identified both by our detector and by the EEG experts is defined as true positive (TP).
S. Yuan et al. / Epilepsy & Behavior 31 (2014) 339–345 Table 2 The results of our detection method for each patient from the evaluation of the performance on a segment-based level. Patient
Sensitivity (%)
Specificity (%)
Recognition accuracy (%)
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 Mean
62.5 100 100 100 100 100 100 95.46 96 89.89 100 100 100 98.01 88.06 95.77 100 100 70 99 100 94.99
99.98 99.15 99.49 98.50 96.95 99.94 99.91 96.30 97.57 92.10 99.91 99.79 99.77 99.31 99.15 99.36 99.85 99.79 99.98 97.42 99.23 98.74
99.94 99.16 99.50 98.52 96.97 99.94 99.91 96.29 97.56 92.08 99.91 99.79 99.78 99.26 99.06 99.33 99.85 99.79 99.94 97.44 99.23 98.73
• Specificity: true negative / the total number of nonseizures marked by the EEG experts. True negative (TN) represents the nonseizures identified by our method and by the experts. • Recognition accuracy: number of correctly classified EEG epochs / the total number of EEG epochs. All results for each patient according to the three statistical measures are depicted on Table 2. It can be observed that the best sensitivity of 100%, specificity of 99.98%, and recognition accuracy of 99.94% were achieved for different patients. In addition, the mean of sensitivity for all the 21 patients was nearly 95.0%, and the means of the specificity and recognition accuracy are both over 98.5%. Furthermore, twelve patients (patients 2, 3, 4, 5, 6, 7, 11, 12, 13, 17, 18, and 21), over half of all the patients, had a sensitivity of 100%. The specificities of all patients are greater than 92%, and among them, there were fifteen patients having specificities higher than 99%. For patient 1, the sensitivity was low because the duration of seizure activities was very short, so some seizure epochs have not been identified effectively. Moreover, considering the feasibility of the proposed method, the event-based approach is also employed to evaluate our algorithm. To assess the system at this level, two statistical measures are calculated: the number of true detected seizures and the false-detection rate. Any detected event by our system that overlapped a seizure event labeled by the EEG experts is defined as a true seizure detection. A sequence of consecutive false positives, which are not adjacent with the seizures, is defined as a false detection. False positive (FP) means the number of seizure segments identified only by our algorithm but not by experts. The system performance on the event-based level is shown in Table 3. In our study, 61 seizures were used to test our method totally, and 58 seizures were detected correctly, and 100% of seizure events were detected for 18 patients. Fig. 3 illustrates an example of a detected seizure which had obvious changes in the used feature. Our method only missed 3 seizures, which come from patients 1, 15, and 19. For patients 1 and 19, the duration of missed seizures is shorter than 12s. The missed seizure of patient 15 is presented in Fig. 4. As shown in this figure, the seizure is very short, and there is no significant difference in the frequency and amplitude between the seizure and nonseizure EEG parts. This seizure is missed because of the unobvious change of diffusion distance features. In total, the average false-detection rate is 0.24/h for all the 21 patients, there were no false detections in the records of approximately half of all patients.
343
Table 3 The results of our detection method for each patient from the evaluation of the performance on an event-based level. Patient
Number of seizures marked by experts
Number of true seizures marked by our system
False detections/h
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 Total
2 2 4 4 3 2 2 1 4 4 2 3 1 3 3 3 4 4 2 4 4 61
1 2 4 4 3 2 2 1 4 4 2 3 1 3 2 3 4 4 1 4 4 58
0 0.22 0 0.11 1.14 0 0 1.17 0.56 0.69 0 0 0 0 0.11 0.10 0 0.25 0 0.50 0.27 0.24
5. Discussion Seizure detection is a difficult task because of the vagaries of ictal EEGs. A perfect detection tool, which assists the clinical evaluation of patients with epilepsy, should be able to reduce the amount of inspected data and have high sensitivity, low false-detection rate. An efficient feature will play a significant role in the entire system. The diffusion distance applied in our study shows its powerful ability to distinguish seizure and nonseizure signals. The significant advantage of this feature is that its computational complexity is much less when compared with other distances, such as the earth mover's distance [14,23]. Diffusion distance is robust to the noise in comparison to the EMD. Furthermore, this feature describes the distance between two different distributions.
Fig. 3. An example of a detected seizure of patient 6. (a) One-hour raw EEG data. (b) The raw EEG data from 47 min to 50.5 min containing one seizure. (c) The diffusion distance features extracted from the one-hour EEG data. The seizure event marked by the EEG experts is between the two vertical lines.
344
S. Yuan et al. / Epilepsy & Behavior 31 (2014) 339–345 Table 4 Comparison of the performance for different methods.
Fig. 4. An example of a missed seizure of patient 15. (a) One-hour raw EEG data. (b) The raw EEG data from 56 min to 57.5 min containing one seizure. (c) The diffusion distance features extracted from the one-hour EEG data. The seizure event marked by the EEG experts is between the two vertical lines.
For every epoch of EEG, it could be transformed into a distribution, and we could calculate the diffusion distance of nonseizure and nonseizure epochs as well as the distance of seizure and nonseizure epochs. It is noticed from Fig. 2 that the diffusion distance of the different categories is larger than that of the same categories, and the diffusion distance can be used to discriminate seizure and nonseizure epochs. Fig. 5 shows the mean value and standard deviation of the diffusion distance. This statistical analysis is achieved by SPSS software. We randomly choose 100 segments from seizure data and nonseizure data, respectively. The figure indicates that the mean value of this feature of seizure signals is much higher than that of nonseizure signals. Results of the independent two-sample T-test indicate that the distinctions in diffusion distance between nonseizure and seizure EEGs are statistically significant (p-values b 0.0001). In addition, the postprocessing part enables us to deal with imprecise boundaries between nonseizure and seizure EEGs and locate the seizure more accurately. These processing technologies all make a contribution to improving the precision of classification. In view of the application of wavelet decomposition to eliminate high-frequency components, it may miss the seizure onset
Method
Sensitivity (%)
Number of seizures used
Number of channels used
Multistage seizure detection [24] A fuzzy logic algorithm [25] Differential windowed variance (DWV) algorithm [26] Improved patient-specific seizure detection [27] Our proposed algorithm
87.5 95.8 91.525
24 56 59
3 6 3
78
63
–
94.99
61 (test) + 26 (training)
3
and hence affect the delay between seizure onset and detection. A few short seizures with predominantly high-frequency discharges may also be missed for this reason. However, the collar technique in the postprocessing part compensates for the missing part of a seizure to make the seizure detection more accurate. Moreover, although keeping the high-frequency components is helpful for detecting subclinical electrographic events and short seizures with predominantly high-frequency discharges, the existence of high-frequency artifacts, such as muscle activities and the power line interferences, could result in the obvious increase of the false-detection rate. So this work employs wavelet decomposition to suppress those high-frequency artifacts and achieve a low false-detection rate. The EEG database used in this study has been applied in other seizure detection systems [24–27]. Raghunathan proposed a multistage seizure detection which used a discrete wavelet transform-based filtering block and a “feature extraction block” containing coastline and variance energy features. The average sensitivity was 87.5% from 5 patients with 24 seizures [24]. Majumdar and Vardhan proposed a differential windowed variance (DWV) algorithm to automatically detect seizure onsets in continuous ECoG (depth-EEG) signals [26]. The sensitivity of this method was 91.525% with 59h of seizure from 15 patients with epilepsy. In both of their methods, the channels they used were the same as the ones we used in this study. Chua developed an improved patient-specific seizure detection by applying a subject-independent quadratic discriminant classifier which was tested on 63 seizures of 15 patients and obtained a sensitivity of 78% [27]. The number of channels was not clearly mentioned in Chua's study. Table 4 depicts a comparison on the results between their methods and our method. Compared to other automatic seizure-detection systems, our proposed algorithm yields a higher sensitivity. Rabbi and Fazel Rezai presented a multistage fuzzy logic algorithm for epileptic seizure onset detection using features of amplitude, frequency, and entropy. The method was evaluated on the same iEEG database from 20 patients with 56 seizures. A sensitivity of 95.8% and a false-detection rate of 0.26/h were achieved [25]. Out of 56 seizures analyzed, their system detected 54 seizures, and 2 seizures were missed. The data from patient 10 (of the FSPEEG database) were discarded because of the presence of electrode movement artifacts. In addition, they employed all the six channels, but we used only three. Compared to this system, our method which evaluated 61 seizures from all of the 21 patients yields a comparable sensitivity.
6. Conclusion
Fig. 5. The mean and standard deviation of the diffusion distance feature extracted from EEG signals.
In this work, we presented a novel method of seizure detection using diffusion distance and BLDA on intracranial EEGs. The diffusion distance, as a novel and effective feature, was extracted to distinguish seizure and nonseizure EEG epochs. The Bayesian linear discriminant analysis was used to classify seizure and nonseizure epochs. The algorithm has been validated on an intracranial EEG database, and the experimental results show that this proposed method achieves a promising performance with a high sensitivity of 94.99% and a low false-detection rate
S. Yuan et al. / Epilepsy & Behavior 31 (2014) 339–345
of 0.24/h. Our study indicates the potential clinical application of the proposed method. In future work, we will focus on integrating more available features to improve the performance and robustness of this seizure detector. Acknowledgements The authors Shasha Yuan, Weidong Zhou, Qi Yuan, Yanli Zhang, and Qingfang Meng are also with the Suzhou Institute of Shandong University, Suzhou, PR China. The support of the Program of Science and Technology of Suzhou (No. ZXY2013030) and the Independent Innovation Foundation of Shandong University (No. 2012DX008) is gratefully acknowledged. References [1] Hauser WA, Annegers JF, Rocca WA. Descriptive epidemiology of epilepsy: contributions of population-based studies from Rochester, Minnesota. Mayo Clin Proc 1996;71:578–86. [2] Stam CJ, Pijn JP, Suffczynski P, Lopezda Silva FH. Dynamics of the human alpha rhythm: evidence for non-linearity? Clin Neurophysiol 1999;110(10):1801–13. [3] Iasemidis LD, Shiau DS, Chaovalitwongse W, Sackellares JC, Pardalos PM, Principe JC, et al. Adaptive epileptic seizure prediction system. IEEE Trans Biomed Eng 2003;50(5):616–27. [4] Shoeb A, Edwards H, Connolly J, Bourgeois B, Treves ST, Guttag J. Patient-specific seizure onset detection. Epilepsy Behav 2004;5(4):483–98. [5] Aarabi A, Fazel-Rezai R, Aghakhani Y. A fuzzy rule-based system for epileptic seizure detection in intracranial EEG. Clin Neurophysiol 2009;120(9):1648–57. [6] Khan YU, Gotman J. Wavelet based automatic seizure detection in intracerebral electroencephalogram. Clin Neurophysiol 2003;114(5):898–908. [7] Gotman J. Automatic recognition of epileptic seizures in the EEG. Electroencephalogr Clin Neurophysiol 1982;54(5):530–40. [8] Ocak H. Automatic detection of epileptic seizures in EEG using discrete wavelet transform and approximate entropy. Expert Syst Appl 2009;36:2027–36. [9] Murro AM, King DW, Smith JR, Gallagher BB, Flanigin HF, Meador K. Computerized seizure detection of complex partial seizures. Electroencephalogr Clin Neurophysiol 1991;79:330–3.
345
[10] Srinivasan V, Eswaran C, Sriraam NH. Approximate entropy-based epileptic EEG detection using artificial neural networks. IEEE Trans Inf Technol Biomed 2007;11(3): 288–95. [11] Khan YU, Rafiuddin N, Farooq O. Automated seizure detection in scalp EEG using multiple wavelet scales. IEEE International Conference on Signal Processing, Computing and Control (ISPCC); 2012. p. 1–5. [12] Subasi A. Epileptic seizure detection using dynamic wavelet network. Expert Syst Appl 2005;29:343–55. [13] Qu H, Gotman J. A patient-specific algorithm for the detection of seizure onset in long-term EEG monitoring: possible use as a warning device. IEEE Trans Biomed Eng 1997;44(2):115–22. [14] Ling H, Okada K. Diffusion distance for histogram comparison. Proc IEEE Comput Soc Conf Comput Vis Pattern Recognit 2006;1:246–53. [15] Yan W, Wang Q, Liu Q, Lu H, Ma S. Topology-preserved diffusion distance for histogram comparison. Proceedings of the British Machine Conference; 2007. p. 61.1–61.10. [16] Hoffmann U, Vesin JM, Ebrahimi T, Diserens K. An efficient P300-based brain– computer interface for disabled subjects. J Neurosci Methods 2008;167(1): 115–25. [17] Lei X, Yang P, Yao D. An empirical Bayesian framework for brain–computer interfaces. IEEE Trans Neural Syst Rehabil Eng 2009;17(6):521–9. [18] Freiburg seizure prediction project. Freiburg, Germany 2008. [On-line]. Available: http://epilepsy.uni-freburg.de/freiburg-seizure-prediction-project/eeg-database. [19] Mallat S. A theory for multiresolution signal decomposition: the wavelet representation. IEEE Trans Pattern Anal Mach Intell 1989;11(7):674–93. [20] MacKay DJC. Bayesian interpolation. Neural Comput 1992;4(3):415–47. [21] Temko A, Thomas E, Marnane W, Lightbody G, Boylan G. EEG-based neonatal seizure detection with support vector machines. Clin Neurophysiol 2011;122: 464–73. [22] Grewal S, Gotman J. An automatic warning system for epileptic seizures recorded on intracerebral EEGs. Clin Neurophysiol 2005;116(10):2460–72. [23] Rubner Y, Tomasi C, Guibas LJ. The earth mover's distance as a metric for image retrieval. Int J Comput Vis 2000;40(2):99–121. [24] Raghunathan S, Jaitli A, Irazoqui PP. Multistage seizure detection techniques optimized for low-power hardware platforms. Epilepsy Behav 2011;22:61–8. [25] Rabbi AF, Fazel-Rezai R. A fuzzy logic system for seizure onset detection in intracranial EEG. Comput Intell Neurosci 2012;2012:705140. [26] Majumdar K, Vardhan P. Automatic seizure detection in ECoG by differential operator and windowed variance. IEEE Trans Neural Syst Rehabil Eng 2011;19(4): 356–65. [27] Chua EC, Patel K, Fitzsimons M, Bleakley CJ. Improved patient specific seizure detection during pre-surgical evaluation. Clin Neurophysiol 2011;122:672–9.
