Automatic Epileptic Seizure Detection in EEG Signals Using Multi-Domain Feature Extraction and Nonlinear Analysis

Abstract

Epileptic seizure detection is commonly implemented by expert clinicians with visual observation of electroencephalography (EEG) signals, which tends to be time consuming and sensitive to bias. The epileptic detection in most previous research suffers from low power and unsuitability for processing large datasets. Therefore, a computerized epileptic seizure detection method is highly required to eradicate the aforementioned problems, expedite epilepsy research and aid medical professionals. In this work, we propose an automatic epilepsy diagnosis framework based on the combination of multi-domain feature extraction and nonlinear analysis of EEG signals. Firstly, EEG signals are pre-processed by using the wavelet threshold method to remove the artifacts. We then extract representative features in the time domain, frequency domain, time-frequency domain and nonlinear analysis features based on the information theory. These features are further extracted in five frequency sub-bands based on the clinical interest, and the dimension of the original feature space is then reduced by using both a principal component analysis and an analysis of variance. Furthermore, the optimal combination of the extracted features is identified and evaluated via different classifiers for the epileptic seizure detection of EEG signals. Finally, the performance of the proposed method is investigated by using a public EEG database at the University Hospital Bonn, Germany. Experimental results demonstrate that the proposed epileptic seizure detection method can achieve a high average accuracy of 99.25%, indicating a powerful method in the detection and classification of epileptic seizures. The proposed seizure detection scheme is thus hoped to eliminate the burden of expert clinicians when they are processing a large number of data by visual observation and to speed-up the epilepsy diagnosis.

1. Introduction

Epilepsy is one of the most common neurological disorders [1]. An estimated 0.6–0.8% of the world’s population or around 50 million people worldwide suffer from epilepsy [2]. Its sudden and recurrent seizures can result from large numbers of neurons going through an excessive and synchronous electrical discharge. Epileptic seizures can be classified into two types: partial and generalized [3]. Partial epileptic seizures occur when a local region of the brain experiences excessive and synchronous electrical discharge, while in a generalized epileptic seizure, the entire brain experiences excessive or synchronous electrical discharge. Both types of epileptic seizures can occur at all ages and are especially prominent in younger and older demographics [4]. Thus, how to diagnose and predict epileptic seizures effectively is still a challenging problem.

Electroencephalography (EEG) is a useful method to monitor the nonlinear electrical function of the brain’s nerve cells; thus, it is a valuable tool for the epilepsy evaluation and treatment [5]. Observation shows that EEG signals include spikes, sharp waves and spike-and-wave complexes not only during a seizure, but also a short time before and between seizures [6]. Currently, the general detection of epileptic activity requires manual scanning of EEG recordings, which commonly takes several days to complete [7]. This epileptic activity detection method is time consuming and error prone [8]. As a result, it is necessary to develop robust and reliable techniques for epileptic activity detection in EEG signals [3,4,9,10].

Recently, several automated methods for detecting the epileptic activity have been proposed [7,10,11]. Most of them employed Fourier spectral analysis for EEG signal extraction based on the assumption that EEG signals are stationary [12], allowing for the transformation of signals from the time to the frequency domain. These methods allow researchers to study the information unique to the frequency domain [12]. EEG signals are typically nonstationary time series and provide only the information available in the time or the frequency domain. However, the prior literature has manifested that the frequency components in EEG signals would change over time, i.e., the EEG signals should not be regarded as having stationary properties [13,14]. Therefore, the time-frequency analysis methods are required to eliminate such a dilemma. A popular short time Fourier transform (STFT) method has been created where the window will move in time to calculate the spectrum density of EEG signals [8,15]. Additionally, the wavelet transform approaches to time-frequency estimation are commonly attractive for EEG signal analysis. For example, the discrete wavelet transform (DWT) method is a classical time-frequency analysis method similar to the short time Fourier transform and has also been used to extract features from EEG signals [7,16]. Faust et al. employed the DWT-based EEG denoising method and the feature extraction for the seizure detection and epilepsy diagnosis [17], and their results showed that the wavelet technique was an effective method for automatic epilepsy diagnosis using EEG signals. Additionally, Hassan et al. further proposed the tunable factor wavelet transform scheme for the automated epilepsy diagnosis. Their seizure detection method can alleviate the burden of medical professionals and help speed up the epileptic seizure diagnosis [4]. However, these approaches, such as the traditional spectrogram STFT, also suffer from the tradeoff issue between time and frequency resolutions. For example, in those traditional approaches, fast Fourier transform (FFT) is used for short overlapping sequences that are assumed to be stationary [3,8]. Additionally, the extracted features, either time or frequency features by wavelet transform from the EEG signals, cannot independently provide an effective way to identify the epileptic seizure activities, without considering the highly complex and nonlinear EEG signals [18]. To some extent, this may deteriorate the performance of identifying EEG signals.

Apart from the time-frequency feature extraction, nonlinear analysis of EEG signals has also received considerable attention to detect epileptic seizures that can be seen as a transition of the human brain [1,19]. Particularly, the nonlinear methods based on the Lyapunov exponent [17], higher order spectra (HOS) [1], information theory and entropy and intrinsic mode functions (IMF) have generally been adopted to analyze and extract nonlinear features for the epileptic seizure diagnosis from EEG signals [10,20]. For example, Acharya et al. extracted the nonlinear HOS features, approximation entropy and sample entropy from the EEG segments, and different classifiers were then used to evaluate the classification performance among normal, interictal and ictal EEG signals [1,21]. Experimental results demonstrated that the nonlinear feature technique is capable of classifying the EEG segments. Additionally, as EEG signals are highly nonlinear and nonstationary, Qu et al. applied the empirical mode decomposition (EMD) to decompose EEG signals into a collection of intrinsic mode functions (IMFs), which is capable of representing nonlinear and nonstationary processes [10]. Furthermore, the IMF features generated by the empirical mode decomposition method were extracted and fed into the support vector machines (SVMs) for epileptic seizure classification of EEG signals. The classification results indicated that the IMFs can achieve high classification accuracy. However, few studies have adopted multimodal features, such as the combination of multi-domain features and nonlinear features, for EEG epileptic activity classification. Therefore, we may integrate these different feature models in a unified framework to include their advantages and try to exclude their limitations for early epilepsy diagnosis.

This paper presents an automatic technique to detect the epileptic activity in EEG signals using multi-domain features and nonlinear analysis to improve the performance of EEG epileptic seizure detection. Specifically, the wavelet threshold method is first applied to remove noise components from the EEG signals prior to EEG signal feature extraction. Then, multiple features are extracted in the time, frequency and time-frequency domains using plural discrete wavelet transforms, as well as the nonlinear analysis, such as IMF, spectrum entropy and approximate entropy (ApEn), from five frequency sub-bands of clinical interest in order to increase the classification robustness and accuracy of the epileptic seizure detection. Furthermore, a dimensionality reduction algorithm of the principal component analysis (PCA), along with the feature ranking method of the analysis of variance (ANOVA), is applied to eliminate irrelevant or redundant features, which may deteriorate the classification performance in the original high-dimensional feature space. The novelty of this paper lies in the combination of multi-domain features and nonlinear features, which increase the classification accuracy to 99.25%. Such a high accuracy indicates that our proposed classification framework could efficiently assist physicians in the detection and classification of brain abnormalities, including epileptic seizures.

The remainder of this paper is organized as follows. Section 2 presents the description of the experimental data used in this work, elucidates a brief description of the feature extraction strategy and explicates the classification model used in the proposed method. Experimental results and the comparison with related literature are presented in Section 3. Finally, the conclusion of this work is given in Section 4.

2. Methodology

EEG signals with weak and low amplitude suffer from complex low frequency noise such as system interference. Thus, the EEG signal pre-processing with noise removal is especially regarded as a significant step for the epileptic seizure analysis and detection [22,23]. Hence, in this work, a wavelet threshold denoising method is employed, which has a superior performance compared to the Fourier transform denoising method [19]. Specifically, the fourth-order Daubechies (db4) wavelet is selected due to its good local approximated performance for nonstationary signals [19,24]. Five frequency sub-bands of clinical interest are then obtained by using the wavelet decomposition and reconstruction: delta (0–4 Hz), theta (4–8 Hz), alpha (8–16 Hz), beta (16–32 Hz) and gamma (32–64 Hz). Herein, wavelet features of its good localizing properties are extracted from each sub-band in the time-frequency domain, followed by a well-known PCA algorithm of the dimensionality reduction in order to remove the irrelevant or spurious features. Finally, different classical classifiers including K-nearest-neighbor (KNN), linear discriminant analysis (LDA), naive Bayesian (NB), logistic regression (LR) and SVM, are used to identify the epileptic seizures from EEG signals. A detailed flowchart of our proposed classification framework is shown in Figure 1.

2.1. Materials

EEG signals used in this work are recorded at the University Hospital Bonn, Germany [18]. The dataset consists of five subsets (denoted as S, F, N, O and Z) recorded with the same 128-channel amplifier system and 12-bit analog-to-digital convertor. Each of the subsets includes 100 segments with a sampling frequency of 173.61 Hz and a duration of 23.6 s, i.e., 4096 sample points, the corresponding Nyquist frequency bandwidth is 86.8 Hz. EEG samples in Datasets O and Z are obtained from five healthy volunteers through external surface electrodes for open and closed eye conditions, while Datasets S, F and N consists of EEG segments recorded from epileptic patients using intracranial electrodes to monitor interictal and ictal epileptic activity. The datasets F and N are acquired during seizure free intervals, while the dataset S only contains the seizure activity. All five Datasets S, F, N, O and Z are tested and categorized into two separate groups using our proposed method. The epileptic seizure (S) class is composed of subset S, and the non-seizure (FNOZ) class includes Subsets F, N, O and Z, respectively.

2.2. Signal Pre-Processing: Wavelet Threshold De-Noising

Generally, physiological signals are contaminated or even distorted by artifacts. Removing noise or artifacts (denoising) is an important step in medical signal processing and analysis. Recently, most automated techniques for signal noise removal are based on the statistical distribution and spectrum characteristics. Particularly, the analysis in the time-frequency domain produces superior signal processing results because it contains both the frequency and time components [25]. For example, the discrete wavelet transform (DWT) allows the signal magnitude to be amplified as larger DWT coefficients, restricting the noise in the total time-frequency domain. In this way, the DWT coefficients of the noise are smaller than the desired signals. Figure 2 presents the block diagrams to depict the decomposition of the input signal using DWT up to the fourth level. Specifically, the original EEG signal (0–64 Hz) is firstly decomposed into its higher frequency part (32–64 Hz) and lower frequency part (0–32 Hz), i.e., the detail and the approximation of the signal at the first level. Then, the approximation of the first decomposition level is additionally divided into its higher (16–32 Hz) and lower (0–16 Hz) frequency part, i.e., the detail and approximation at the second level. Thus, the wavelet threshold method can perform well in denoising nonstationary EEG signals, which is defined by [26]:

$$\lambda =\sigma \sqrt{2\log N}$$

(1)

where $\lambda$ is the wavelet threshold, $\sigma$ is the standard deviation of the noise and N is the length of the sample signals, respectively. Since the noise buried in the EEG signals is practically white noise and has sharp waves [27], it is effective to adopt the wavelet threshold method to remove the noise from EEG signals.

2.3. Feature Extraction

2.3.1. Feature Extraction in the Time Domain, Frequency Domain and Time-Frequency Domain

Features extracted by analyzing the EEG signals solely in the time domain may omit the important frequency information and vice versa. However, the signal processing technique using wavelets can eliminate such deficiency and extract effective features in the time-frequency domain well. Specifically, a signal can be represented as a linear or nonlinear combination of basis functions. In the time-frequency domain, a function $\psi$ formed with limited power, i.e., finite duration, and zero mean can be used as the basis function of the wavelet [28], which is satisfied with:

$$\sum _{N=-\infty }^{\infty }{\left|\psi \left[N\right]\right|}^2<\infty ,\sum _{N=-\infty }^{\infty }\psi \left[N\right]=0$$

(2)

where N is the length of the sample data. The wavelet can also be moved over time by the parameter b and scaled by the dilation parameter a, which is defined as follows:

$${\psi }_{a,b}\left[N\right]=\frac{1}{\sqrt{a}}\psi \left[\frac{N-b}{a}\right]$$

(3)

It should be noted that the wider wavelets (with larger dilation parameter a) are more suitable for extracting slower changes, while the narrower wavelets (with smaller dilation parameter a) can be used to extract faster changes well. Features can be extracted by analyzing EEG signals through the basis function of wavelets, which contain information in both the time and frequency domain. Modifying the parameters a and b allows for an easy calculation of the wavelet transform coefficients in the following formula:

$$w_{a,b}\left[n\right]=\sum _{\tau =1}^{N}x\left[\tau \right]{\psi }_{a,b}\left[n-\tau \right],1\le n\le N$$

(4)

where $x\left[\tau \right]$ is the sample signal with the length of N samples. Note that it is significant to determine the wavelet function type and level before the application of DWT. In this work, the db4 wavelet is employed for its perfect local characteristics in the time-frequency domain [28], because the db4 wavelet’s shape and smoothening feature have the superior approximated performance in the field of nonstationary signal processing techniques. In addition, the EEG signals are decomposed into four levels, resulting in five clinically-significant sub-bands. The standard deviation and the relative power of the DWT transform coefficients from each of five clinical interest sub-bands are extracted as features in the time-frequency domain.

In the time domain, statistic features, including the mean, the variance, the coefficient of variation and the total variation, are extracted. Additionally, the maximum, minimum, mean and the total variation measures of the DWT transform coefficients are also estimated in order to describe the nonstationary signals in detail in the time domain, where the total variation measure is defined as:

$$V_x=\frac{1}{N-1}\frac{{\sum }_{n=2}^N\left|x\left[n\right]-x\left[n-1\right]\right|}{ma{x}_x-mi{n}_x}$$

(5)

where $ma{x}_x$ and $mi{n}_x$ are the maximum and minimum of the signal x, respectively. Note that the value of the total variation $V_x$ is in the range between $\frac{1}{N-1}$ for slow signals and 1 for fast changing signals.

As for the frequency domain, the fast Fourier transform (FFT) is an effective common practice for signal analysis with different frequencies, which cannot be identified in the time domain. Herein, the relative power spectral density estimated by the coefficients of the FFT is extracted, representing the signals’ features in the frequency domain.

Such features of the time-frequency, time and frequency domain are tested for the discriminative power in the epileptic signal classification using the following machine learning procedure.

2.3.2. Feature Extraction via Nonlinear Analysis

The EEG signals are extremely nonstationary and nonlinear. Thus, features extracted via only linear analysis techniques can provide an insufficient approximation of the nonlinear signals. Therefore, nonlinear analysis of EEG signals is more appropriate and reliable.

The empirical mode decomposition (EMD) is an effective method to decompose nonstationary signals into a sum of intrinsic mode functions (IMFs) and has been successfully applied in various applications [24,29]. Each IMF should satisfy two conditions. Firstly, the number of extrema and the number of zero crossings differ at most by one. Secondly, the envelope defined by the local maxima and the envelope defined by the local minima should have a zero mean at any point. IMFs can be extracted from EEG signals $x\left[n\right]$ by the following steps [29]:

• Procedure 1, IMF extraction procedure:
  • Extract local max and local min magnitudes from signal $x\left[n\right]$;
  • Obtain the envelope ${\epsilon }_{max}\left[n\right]$ by connecting all of the maximums with cubic spline interpolation, and similarly obtain the envelope ${\epsilon }_{min}\left[n\right]$ by connecting all of the minimums with cubic spline interpolation;
  • Compute the average of ${\epsilon }_{max}\left[n\right]$ and ${\epsilon }_{min}\left[n\right]$, and denote as $\alpha \left[n\right]$:

    $$\alpha \left[n\right]=\frac{{\epsilon }_{max}\left[n\right]+{\epsilon }_{min}\left[n\right]}{2}$$
  • Extract the detail $d\left[n\right]$ from $x\left[n\right]$ as:

    $$d\left[n\right]=x\left[n\right]-\alpha \left[n\right]$$
  • Check whether the detail $d\left[n\right]$ satisfy the above conditions mentioned for IMF or not;
  • Repeat Steps 1–5, until $d\left[n\right]$ satisfies the conditions for IMF.

Once the first IMF is derived through the procedure described above, which is defined by $v_1\left[n\right]=d\left[n\right]$, another IMF can be extracted similarly in terms of the residue $r\left[n\right]=x\left[n\right]-c_1\left[n\right]$. This procedure is repeated until the residue reaches the tolerance or no more IMF can be extracted. Thus, the original signal $x\left[n\right]$ can be represented by:

$$x\left[n\right]=\sum _{m=1}^M{v}_m\left[n\right]+r_M\left[n\right]$$

(6)

where M is the number of IMFs, $v_m\left[n\right]$ is the m-th IMF and $r_M\left[n\right]$ is the final residue. In this work, the EEG signal is found to be represented by 8 IMFs, i.e., $M=8$.

In order to better characterize the dynamic IMFs, a phase space reconstruction (PSR) approach is applied here to search the feature patterns in a higher dimension. The PSR is used to reconstruct the so-called phase space, which is a very useful tool for characterizing any low-dimensional or high-dimensional dynamic system. Alternatively, a point in the phase space represents the state of the system at any given time [11]. For an IMF $v\left[n\right]$, where $n=1,2,\cdots ,N$, the phase space can be reconstructed by [11]:

$$S_j=\left(v\left[j\right],v\left[j+\tau \right],v\left[j+2\tau \right],\cdots ,v[j+\left(m-1\right)\tau ]\right)$$

(7)

where $j=1,2,\cdots ,N-\left(m-1\right)\tau$, m is the embedding dimension and $\tau$ is the time delay, which is chosen to be 1. Practically, the embedding dimension is confined to two or three to preserve visual simplicity. Thus, the 2D PSR of IMF is a plot of $v\left[j\right]$ and $v\left[j+1\right]$, which is an elliptical area [30]. The major axis radius $r_a$ and the minor axis radius $r_b$ can then be calculated by [30]:

$$\left\{\begin{array}{c}{M}_x=\sqrt{\frac{1}{N-1}{\displaystyle \sum _{k=1}^{N-1}}{v}^2\left[k\right]}\hfill \\ M_y=\sqrt{\frac{1}{N-1}{\displaystyle \sum _{k=1}^{N-1}}{v}^2\left[k+1\right]}\hfill \\ M_{xy}=\frac{1}{N-1}{\displaystyle \sum _{k=1}^{N-1}}v\left[k\right]v\left[k+1\right]\hfill \\ L=\sqrt{\left(M_x^2+M_y^2\right)-4\left(M_x^2{M}_y^2-M_{xy}^2\right)}\hfill \\ r_a=\sqrt{2\left(M_x^2+M_y^2+L\right)}\hfill \\ r_b=\sqrt{2\left(M_x^2+M_y^2-L\right)}\hfill \end{array}\right.$$

(8)

where $M_x$, $M_y$ and $M_{xy}$ are the standard deviation of the IMF on the major and minor axes (denoted as x- and y-axis) of the two-dimensional phase space and the covariance, respectively.

The 3D PSR is a plot of $v\left[j\right]$, $v\left[j+1\right]$ and $v\left[j+2\right]$, which also forms an elliptical ball. The distance of the point ($v\left[j\right]$, $v\left[j+1\right]$, $v\left[j+2\right]$) from the origin (0, 0, 0) is defined by [11]:

$$E_j=\sqrt{{v\left[j\right]}^2+v^2\left[j+1\right]+v^2\left[j+2\right]}$$

(9)

For each value in the 3D PSR space, the distance $E_j$ is calculated. Additionally, the first and third quarter of the distance set $\left\{E_j\right|j=1,2,\cdots ,N-2\}$ are extracted as features (denoted as $l_c$ and $l_d$, respectively), representing the data variability in 3D PSR space [31]. Hence, for each IMF, a feature vector of $(r_a,r_b,l_c,l_d)$ is extracted to represent the 2D and 3D PSR. With $M=8$, a total number of $4\times 8=32$ features is obtained from all IMFs in this work.

Additionally, since Shannon first developed the concept of information entropy as part of information theory [32], various entropies have been developed and widely applied to quantify the dynamics of the nonstationary signal. Such entropies can be broadly classified into two categories as embedding entropies (i.e., approximate entropy) and spectral entropies (i.e., wavelet entropy).

Approximate entropy: The approximate entropy (ApEn) is a measure of data regularity, measuring the predetermination of the current amplitude based on the previous amplitude [30,31]. The ApEn is the likelihood that proximal data trends will remain close together in the next pattern. A high regularity produces smaller ApEn values, while a low regularity produces higher ApEn values. Thus, it can be utilized to represent the overall complexity and predetermination of the EEG signals. Note that the ApEn is scale invariant and model independent. Herein, the approximate entropy of each sub-band EEG signal can be calculated by the following steps:

• Procedure 2, approximate entropy calculation:
  • Let the values containing N samples in each sub-band be $X=\left[x\left(1\right),x\left(2\right),x\left(3\right),\cdots ,x\left(N\right)\right]$;
  • Let $X_m\left(\mathrm{i}\right)$ be a sub-sequence of X such that $X_m\left(\mathrm{i}\right)=[x\left(i\right),x\left(i+1\right),x\left(i+2\right),\cdots ,$$x\left(i+m-1\right)]$ for $1\le i\le N-m+1$, where m is the length of the sub-sequence;
  • Let r represent the noise filter level that is defined as [33]:

    $$r=k\times SD,k=0,0.1,\cdots ,0.9$$

    (10)

    where SD is the standard deviation of the data sequence X.
  • Let $X_m\left(\mathrm{i}\right)$ represent a set of sub-sequences obtained from $X_m\left(\mathrm{j}\right)$ by varying j from 1–($N-m$). Each sequence $X_m\left(\mathrm{j}\right)$ in the set of $X_m\left(\mathrm{j}\right)$ is compared with $X_m\left(\mathrm{i}\right)$, and in this process, two parameters, namely $C_i^m\left(r\right)$ and $C_i^{m+1}\left(r\right)$, are defined as follows:

    $$C_i^m\left(r\right)=\frac{{\sum }_{j=1}^{N-m}{k}_j}{N-m+1},$$

    (11)

    $$\mathrm{where}k=\left\{\begin{array}{c}1,if\left|x\left(i+t\right)-x\left(j+t\right)\right|<rfor0\le t<m\hfill \\ 0,otherwise\hfill \end{array}\right.$$
  • The approximate entropy is calculated by using $C_i^m\left(r\right)$ and $C_i^{m+1}\left(r\right)$ as follows:

    $$ApEn=\frac{1}{N-m+1}{\displaystyle \sum _{i=1}^{N-m+1}}\ln C_i^m\left(r\right)-\frac{1}{N-m}{\displaystyle \sum _{i=1}^{N-m}}\ln C_i^{m+1}\left(r\right)$$

    (12)

Spectral entropies: The spectral entropies use the components of the power spectrum of the signal as the probabilities to describe the variability within the signal. In this study, the entropy of the wavelet coefficients (denoted as wavelet entropy) with different scales is calculated to reveal more information carried by the DWT coefficients. Wavelet entropy measures the regularity degree of signals and indicates the dynamical characteristics of nonlinear systems [34]. Detailed information and approximations can be obtained from the given EEG time series X by applying multilevel wavelet decomposition using the db4 wavelet. The power of the EEG signal X can be expressed as the sum of all component (0–64 Hz) powers, which is shown in the following equation:

$$power=\sum _{i=1}^{P}powe{r}_i$$

(13)

where P is the total number of components and $powe{r}_i$ is calculated by the different component coefficients. If the relative power is defined by $p_i=powe{r}_i/power$, the wavelet entropy ($\mathrm{W}EN$) can be calculated as follows:

$$WEN=-\sum _{i=1}^P{p}_i\log p_i$$

(14)

Obviously, the spectral entropy ($SEN$) can be easily calculated by [35]:

$$SEN=-{\displaystyle \sum _{f=0}^{f=64}}{p}_f\log p_f$$

(15)

where $p_f$ is the relative power of the component with frequency f, i.e., $p_f=powe{r}_f/\sum powe{r}_f$.

2.4. Classification and Performance Evaluation

Redundant features are removed by the principle component analysis (PCA) and analysis of variance (ANOVA) statistical test, where the principle components p are selected to be the optimal value ${10}^{-3}$ (FDR adjusted). The calculated features are fed into different classifiers for the epileptic seizure classification. Particularly, the support vector machine (SVM) classifier is a well-known supervised learning method based on the finite sample theory [36]. In small sample cases, traditional classifiers based on the empirical error minimization are prone to generating the overfitting problem, while the SVM is based on the structural risk minimization principle and can ensure a good generalization ability [37]. Therefore, the effectiveness of the selected optimal feature subset is firstly evaluated by a radial basis function (RBF) kernel-based SVM (RBF-SVM). In addition, the other four classical classifiers, i.e., KNN, LDA, NB and LR, are also used to illustrate the effectiveness of the proposed classification framework.

Generally, the classification performance of the proposed framework can be evaluated by statistical measurements of the sensitivity (SEN), specificity (SPE) and accuracy (ACC), which are defined as follows:

$$\begin{array}{cc}\hfill SEN& =\frac{TP}{TP+FN}\times 100\%\hfill \end{array}$$

(16)

$$\begin{array}{cc}\hfill SPE& =\frac{TN}{TN+FP}\times 100\%\hfill \end{array}$$

(17)

$$\begin{array}{cc}\hfill ACC& =\frac{TP+TN}{TP+TN+FP+FN}\times 100\%\hfill \end{array}$$

(18)

where $TP$ and $TN$ represent the total number of correctly-detected true normal events and true seizure events, the $FP$ and $FN$ indicate the total number of erroneously normal events and erroneously seizure events, respectively.

In order to obtain an unbiased estimation of the classification performance, N-fold cross-validation is applied as the training-test strategy. Specifically, input datasets are randomly divided into N equal parts, where $N-1$ parts are used for the feature selection and classifier training and the remaining part for testing the classification performance. Such a procedure is repeated N times, and each time, a different part is left out and tested, then the average result across the total testing parts is calculated as the performance measurement. In this work, taking different training-testing partitioning into consideration, both 5- and 10-fold cross-validation are performed, each with 100 repetitions, i.e., 100-times N-fold cross-validation can develop robust classifiers and an efficient diagnostic system.

3. Results

3.1. Wavelet Threshold Denoising and Feature Extraction

A signal $x\left[t\right]$ for denoising is randomly selected from the datasets S, F, N, O and Z, respectively. For comparison and clarity, the original signals (Figure 3a) and the processed results (Figure 3b) with a length of 10 s are shown below. Additionally, the wavelet decompositions of the original EEG signal and the denoised signal are shown in Figure 4 and Figure 5, respectively. As for the wavelet threshold denoising, all DWT coefficients smaller than the threshold are declined to zero. Particularly, as shown in Figure 5, the detail at the first decomposition level is nearly zero, indicating that the frequency component of 32–64 Hz of the original EEG signal is totally removed and further reflecting that lower frequencies are more common in normal brain states. It has been reported that there is a shift of the EEG signal energy from lower to higher frequency bands before and during a seizure [38].

After the DWT wavelet threshold denoising, signal features are analyzed in the time, frequency, time-frequency domain, as well as the nonlinear analysis using EMD-PSR and entropies. In total, 83 features are extracted for five datasets including a total of 500 EEG segments. Therefore, each original segment of the EEG signals from the time domain can now be presented in the feature space as $X={\left[x_1,x_2,x_3,\cdots ,x_{83}\right]}^T$. Figure 6 shows the p-values (FDR adjusted) for the feature set via ANOVA. Apart from eight features, all features show a statistically-significant difference under the level of ${10}^{-3}$. Specifically, eight features over the level of ${10}^{-3}$ are involved in the total variation of the beta sub-band, the max DWT coefficient of the beta sub-band, $r_b$ of the fifth IMF, $l_c$ of the third and sixth IMFs and $l_d$ of the first, fourth and seventh IMFs, respectively.

3.2. Dimension Reduction in Feature Space

The original features extracted are commonly irrelevant or redundant. In such a case, it is a trade-off to select an optimal feature subset from the original feature set, and some kinds of feature dimensionality reduction, such as the well-known PCA, are employed to improve the classification performance. The PCA algorithm is implemented to obtain a relatively low dimensional, but significantly discriminative feature set. The result of the PCA dimensionality reduction is shown in

Table 1. Dimension reduction of the feature space based on the PCA.

Feature Analyzed	Dimension
	Before	After
Standard Deviation	5	2
Total Variation	5	2
Relative Power FFT	5	2
Standard Deviation & Relative power DWT	10	2
EMD-PSR	32	2
Entropy	11	2
Min, Max, Mean DWT coefficients	15	2

. Clearly, each original feature space is mapped into a reduced two-dimensional feature space by using the PCA method.

3.3. Experiment Classification Results

How to accurately classify EEG signals recorded from the epileptic subjects, healthy controls or epileptic ones during the interictal time is a challenging problem. Thus, it is highly required to achieve a good classification performance of automatically classifying the epileptic seizure activity. In this study, we use the aforementioned EEG datasets (S, F, N, O, Z) to evaluate the classification performance of the proposed epileptic detection framework. Each EEG segment is firstly denoised via the DWT wavelet threshold method and further decomposed into five frequency sub-bands based on the clinical interest. Then, features from the time domain, frequency domain, time-frequency domain and nonlinear analysis are extracted, and the optimal feature set is further selected. Five classifiers based on the single domain features and multi-domain features are used to detect the epileptic seizure EEGs from the non-seizure EEGs. The results of the classification performance measured by the SEN, SPE and ACC, along with the standard deviation, over the total 100 rounds of both 5- and 10-fold cross-validation are shown in

Table 2. The classification performance of 5-fold and 10-fold CV with wavelet denoising method.

Feature Set	Classifiers	5-Fold CV			10-Fold CV
		SEN	SPE	ACC	SEN	SPE	ACC
Time domain	KNN	97.72 ± 0.83	99.74 ± 0.44	98.73 ± 0.43	91.98 ± 1.10	99.95 ± 0.11	98.36 ± 0.24
	LDA	85.88 ± 1.92	100.00 ± 0.00	92.94 ± 0.96	83.56 ± 0.88	99.75 ± 0.04	96.52 ± 0.18
	NB	92.92 ± 1.31	96.56 ± 0.57	94.74 ± 0.74	92.94 ± 0.70	98.47 ± 0.16	97.37 ± 0.18
	LR	95.84 ± 0.86	98.28 ± 0.69	97.06 ± 0.53	95.76 ± 0.62	98.12 ± 0.65	96.94 ± 0.45
	SVM	97.04 ± 0.82	99.46 ± 0.90	98.25 ± 0.61	93.46 ± 1.08	99.95 ± 0.11	98.65 ± 0.23
Frequency domain	KNN	91.86 ± 1.85	89.36 ± 1.66	90.61 ± 1.15	92.68 ± 1.75	91.00 ± 1.23	91.84 ± 1.03
	LDA	74.50 ± 1.53	95.96 ± 0.87	85.23 ± 0.87	75.10 ± 0.92	95.32 ± 0.55	85.21 ± 0.50
	NB	84.04 ± 0.77	91.50 ± 1.14	87.77 ± 0.71	83.42 ± 0.67	90.70 ± 0.88	87.06 ± 0.64
	LR	88.72 ± 1.59	91.44 ± 0.92	90.08 ± 0.80	90.10 ± 0.94	90.50 ± 0.57	90.30 ± 0.54
	SVM	90.02 ± 1.91	90.10 ± 1.06	90.06 ± 1.10	91.02 ± 1.10	89.70 ± 0.90	90.36 ± 0.69
Time-frequencydomain	KNN	93.96 ± 1.98	98.69 ± 0.28	97.75 ± 0.48	94.48 ± 1.42	98.79 ± 0.16	97.92 ± 0.31
	LDA	72.28 ± 1.59	99.65 ± 0.16	94.18 ± 0.34	72.60 ± 1.40	99.66 ± 0.12	94.25 ± 0.29
	NB	86.12 ± 1.44	95.28 ± 0.63	93.45 ± 0.59	86.80 ± 1.04	95.16 ± 0.43	93.48 ± 0.41
	LR	93.08 ± 1.59	98.34 ± 0.74	97.29 ± 0.66	91.48 ± 1.20	93.82 ± 0.86	92.65 ± 0.68
	SVM	94.48 ± 1.53	98.76 ± 0.31	97.90 ± 0.42	94.26 ± 1.34	98.77 ± 0.27	97.87 ± 0.33
Nonlinear analysis	KNN	96.72 ± 1.37	98.96 ± 0.28	98.51 ± 0.33	97.62 ± 0.82	99.03 ± 0.18	98.75 ± 0.24
	LDA	66.38 ± 2.54	99.60 ± 0.22	92.96 ± 0.54	65.10 ± 2.00	99.70 ± 0.17	92.78 ± 0.42
	NB	76.26 ± 1.87	96.46 ± 0.46	92.42 ± 0.54	77.32 ± 1.22	96.67 ± 0.39	92.80 ± 0.43
	LR	91.52 ± 1.55	98.93 ± 0.31	97.45 ± 0.44	94.24 ± 0.86	96.74 ± 1.09	95.49 ± 0.63
	SVM	96.02 ± 1.83	98.98 ± 0.36	98.39 ± 0.44	96.56 ± 0.85	99.10 ± 0.30	98.59 ± 0.30
Multi-domain and nonlinear analysis	KNN	96.12 ± 1.11	99.15 ± 0.17	98.54 ± 0.25	96.58 ± 0.94	99.16 ± 0.16	98.65 ± 0.22
	LDA	90.00 ± 1.30	99.67 ± 0.14	97.74 ± 0.27	91.16 ± 0.95	99.71 ± 0.10	98.00 ± 0.22
	NB	91.58 ± 1.34	95.26 ± 0.61	94.52 ± 0.58	91.48 ± 0.98	95.30 ± 0.46	94.53 ± 0.42
	LR	95.64 ± 1.32	99.44 ± 0.16	98.68 ± 0.30	95.86 ± 1.02	99.51 ± 0.17	98.78 ± 0.26
	SVM	97.04 ± 1.52	99.54 ± 0.22	99.04 ± 0.34	97.98 ± 1.07	99.56 ± 0.20	99.25 ± 0.28

Where bold values indicate the best classification results.

. It can be observed that the classification accuracy by using the multi-domain features and nonlinear analysis method is higher than that of the single domain features or nonlinear analysis. Particularly, for the comparison of the effectiveness between with and without the wavelet denoising method, the classification results of the EEG signals without using the wavelet denoising method are also presented in

Table 3. The classification performance of 5-fold and 10-fold CV without wavelet denoising method.

Feature Set	Classifiers	5-Fold CV			10-Fold CV
		SEN	SPE	ACC	SEN	SPE	ACC
Time domain	KNN	94.12 ± 1.32	98.62 ± 0.60	96.37 ± 0.73	94.64 ± 1.00	98.66 ± 0.65	96.65 ± 0.58
	LDA	80.36 ± 1.05	97.60 ± 0.92	88.98 ± 0.73	79.92 ± 0.77	97.40 ± 0.69	88.66 ± 0.51
	NB	89.90 ± 0.98	96.60 ± 0.57	93.25 ± 0.53	89.82 ± 0.95	96.84 ± 0.37	93.33 ± 0.52
	LR	95.22 ± 0.97	96.84 ± 0.76	96.03 ± 0.64	95.20 ± 0.63	96.82 ± 0.82	96.01 ± 0.60
	SVM	96.14 ± 0.80	98.34 ± 0.84	97.24 ± 0.63	96.00 ± 0.72	98.56 ± 0.57	97.28 ± 0.46
Frequency domain	KNN	83.84 ± 1.68	83.10 ± 1.93	83.47 ± 1.37	84.40 ± 1.11	83.60 ± 1.60	84.00 ± 1.00
	LDA	76.86 ± 0.96	90.96 ± 0.92	83.91 ± 0.66	76.70 ± 0.92	90.38 ± 0.52	83.54 ± 0.54
	NB	79.80 ± 1.22	89.46 ± 0.96	84.63 ± 0.73	79.66 ± 0.79	89.68 ± 0.68	84.67 ± 0.52
	LR	82.16 ± 1.42	87.24 ± 1.16	84.70 ± 0.78	82.14 ± 1.02	87.00 ± 0.96	84.57 ± 0.66
	SVM	85.04 ± 1.56	84.12 ± 1.61	84.58 ± 0.83	85.30 ± 1.27	83.86 ± 1.17	84.58 ± 0.75
Time-frequencydomain	KNN	95.22 ± 0.92	93.78 ± 1.01	94.50 ± 0.79	95.24 ± 0.84	93.54 ± 1.06	94.39 ± 0.65
	LDA	78.72 ± 1.02	98.04 ± 0.60	88.38 ± 0.62	79.00 ± 0.49	98.18 ± 0.52	88.59 ± 0.40
	NB	88.02 ± 1.10	95.52 ± 1.06	91.77 ± 0.76	88.38 ± 0.85	95.44 ± 1.08	91.91 ± 0.75
	LR	91.22 ± 1.01	94.18 ± 1.52	92.70 ± 0.87	91.38 ± 0.60	94.18 ± 1.16	92.78 ± 0.65
	SVM	94.70 ± 1.20	94.72 ± 1.44	94.71 ± 0.93	94.76 ± 0.97	95.44 ± 1.49	95.10 ± 0.90
Nonlinear analysis	KNN	95.60 ± 1.22	94.04 ± 1.26	94.82 ± 0.82	95.78 ± 0.90	95.00 ± 0.85	95.39 ± 0.61
	LDA	91.38 ± 1.44	94.24 ± 1.48	92.81 ± 1.12	91.48 ± 1.02	94.84 ± 1.03	93.16 ± 0.82
	NB	87.48 ± 1.24	94.12 ± 1.21	90.80 ± 0.84	86.74 ± 1.13	93.76 ± 0.86	90.25 ± 0.76
	LR	94.60 ± 1.34	94.82 ± 1.03	94.71 ± 0.83	95.52 ± 1.06	95.62 ± 0.60	95.57 ± 0.66
	SVM	95.04 ± 1.37	94.74 ± 1.11	94.89 ± 0.92	95.22 ± 1.38	94.72 ± 0.94	94.97 ± 0.80
Multi-domain andnonlinear analysis	KNN	96.80 ± 0.72	95.50 ± 1.12	96.15 ± 0.69	96.86 ± 0.63	95.82 ± 0.55	96.34 ± 0.39
	LDA	91.36 ± 1.52	96.28 ± 1.11	93.82 ± 1.11	91.56 ± 0.83	97.02 ± 0.68	94.29 ± 0.51
	NB	92.52 ± 1.22	95.44 ± 1.04	93.98 ± 0.83	92.42 ± 1.02	95.64 ± 0.89	94.03 ± 0.70
	LR	94.90 ± 1.20	96.98 ± 0.99	95.94 ± 0.87	95.16 ± 0.97	97.04 ± 0.89	96.10 ± 0.73
	SVM	95.98 ± 1.12	96.86 ± 1.10	96.42 ± 0.82	96.04 ± 0.85	97.12 ± 0.74	96.58 ± 0.60

. It is obvious that the DWT wavelet threshold method is used to remove the artifacts, and the classification accuracy is obviously improved. The proposed classification framework based on the multi-domain (time domain, frequency domain and time-frequency domain) feature extraction and nonlinear analysis can obtain a robust and high classification accuracy of epileptic seizure EEGs (Class S) between the seizure-free and normal EEGs (Class FNOZ), with an average accuracy of 99.04% (standard deviation 0.34) and 99.25% (standard deviation 0.28) for both 5-fold and 10-fold cross-validation, respectively. Additionally, the result of most classifiers based on features from only the single domain or nonlinear analysis fades next to that of the proposed classification framework based on all features. Furthermore, our results show that the SVM classifier performs better than the other classifiers by registering 99.25% for the classification accuracy and, for both, the sensitivity of 97.98% and the specificity of 99.95% from the 10-fold cross-validation with the wavelet denoised EEGs. It is evident that using the combination of the multi-domain features and nonlinear analysis, the classification results of the average accuracy of 99.25% is reached, indicating that the proposed method in this paper is an accurate tool in classifying the segments accurately.

4. Discussion and Conclusions

The main objective of this work is to efficiently perform an automated classification framework to detect seizure and non-seizure EEG signals. As shown in Table 2, the proposed classification framework based on multi-domain feature and nonlinear analysis is capable of distinguishing the seizure signals, with the overall average classification accuracy of 100-times 10-fold cross-validation reaching 99.25%, where the 10-fold cross-validation is used to achieve the best classification accuracy compared to the 5-fold cross-validation method. In addition, in order to demonstrate the effect of the wavelet denoising method, the epileptic seizure classification task from the non-seizure EEG segments based on the multimodal feature method without and with the wavelet denoising method is implemented, and the comparison of the experimental results is shown in Table 2 and Table 3. From Table 2, we can conclude that the classification results of different feature methods without using the wavelet denoising method can also perform well, which indicates the good classification performance of the proposed method. Furthermore, the classification results can be improved by using the wavelet denoising method, the effectiveness of which is tested by the wavelet threshold denoising method.

Table 4. The comparison of classification performance from different methods for the same data set.

Problem	Authors	Methods	Accuracy
S-FNOZ	Tzalla et al. [3]	Time-frequency analysis, artificial neural network	97.73%
	Guo et al. [39]	Multiwavelet transform, MLPNN	98.27%
	Rivero et al. [42]	Time frequency analysis, KNN	98.40%
	Kaleem et al. [40]	Variation of empirical mode decomposition	98.20%
	Kai Fu et al. [10]	HMS analysis, SVM	98.80%
	Niknazar M et al. [43]	Wavelet transform, RQA, ECOC	98.67%
	Musa Peker et al. [41]	Dual-tree complex wavelet transform, complex-valued neural networks	99.15%
	Jaiswal et al. [44]	Local neighbor Descriptive pattern, artificial neural network	98.72%
	This work	DWT, multi-domain feature extraction and nonlinear analysis	99.25%

Where bold values indicate the best classification results.

displays the comparison of the classification results between the proposed multimodal feature method and the reported state-of-the-art approaches in the previous literature, where the classification performance is evaluated by the identical EEG signals, and the comparison of the classification results is more feasible and persuasive. It is notable that the majority of existing techniques is focused on feature extraction only in a single domain, such as the time domain, frequency domain, time-frequency domain, or nonlinear analysis. For example, Tzalla et al. [3] proposed a classification framework of the seizure based on the time-frequency analysis of the whole EEG signals and neural network and obtained an average accuracy of 97.72%. Guo et al. [39] applied the multi-wavelet transform to better represent the EEG signals in the time-frequency domain; their accuracy of seizure classification was only 98.27% based on approximate entropy feature and neural network. Additionally, Kaleem et al. [40] and Fu et al. [10] implemented the nonlinear analysis of the empirical mode decomposition on the EEG signals and achieved an accuracy of 98.20% and 98.80% on the seizure classification, respectively. Their work indicated that features from different domains are capable of discriminating the seizure EEG signals. Furthermore, Peker et al. [41] proposed a method for diagnosing epilepsy based on complex-valued classifiers. The classification accuracy of 99.15% was achieved in their proposed hybrid method with 10-fold cross-validation. However, multi-domain-based features are more comprehensive and more prone to achieve a superior performance of seizure detection, with the classification accuracy in our proposed method reaching 99.25%. In Table 4, it is obvious that our proposed framework can obtain better classification results than the rest of the existing reported results.

The evaluated method proposed in this paper does achieve a high accuracy and can thus be used to detect the epileptic seizure activity in the EEG signals. The high classification accuracy is attributed to the highly discriminative features. The proposed multimodal feature method better captured the nonstationary and nonlinear EEG signals. The SVM classifier has superior classification capabilities compared to other classifiers, thus, it leads to the highest accuracy in our work. In addition, the proposed classification framework is fully automated and can be easily implemented as a software application for the clinical diagnosis. The results are objective. However, the proposed classification scheme may lead to higher computational complexity than the single domain feature or nonlinear analysis methods because a large number of redundant features is fed into the classifiers for epileptic seizure activities. Additionally, this effective detection framework may need to be tested by large-scale clinical trials with EEG signals. In order to validate the effectiveness of our proposed technique, the proposed scheme should be valuated by some sparse learning algorithms and a larger multi-ethnic database, so that more bench marking and diversification of data collection will be obtained.

In conclusion, this paper proposes an effective classification framework based on wavelet threshold denoising, multi-domain feature extraction, nonlinear analysis and PCA dimensionality reduction for the automatic detection of epileptic seizures. Experimental results of the proposed method on a dataset of five categories (S and FNOZ) outperform the existing reported results in the literature with an average accuracy of 99.25%. The classification performance is obtained largely due to the combination of features from various domains of time, frequency, time-frequency and nonlinear features, outperforming the methods based on the single-domain feature extraction and classification. We are planning to study larger databases to evaluate our proposed technique and further to establish the significance of the discriminative features discovered in our work.