Epileptic seizure detection based on the kernel extreme learning machine

2.1Kernel ELM

Suppose there are N arbitrary samples (𝒙i,𝒕i), where 𝒙i=[xi⁢1,xi⁢2,⋯,xi⁢n]T∈𝑹n, 𝒕i=[ti⁢1,ti⁢2,⋯,ti⁢m]T∈𝑹m. Thus, standard SLFNs with L hidden nodes can be mathematically expressed:

Figure 1.

Single hidden layer feedforward network.

Where 𝜷i=[βi⁢1,βi⁢2,⋯,βi⁢m]T is the weight vector connecting the i th hidden node and the output nodes. 𝒂i=[ai⁢1,ai⁢2,⋯,ai⁢n]T Is the weight vector connecting the i th hidden node and the input nodes, and bi is the threshold of the i th hidden node. gi denotes the output function G⁢(𝒂i,bi,𝒙) of the i th hidden node (cf. Fig. 1). where

The least-squares method is utilized to obtain weight 𝜷 by solving the following equations.

The unique solution is

Where 𝑯+ is the Moore-Penrose generalized inverse 𝑯. The training process of the ELM algorithm is outlined below.

If H⁢HT is nonsingular, to improve the stability of ELM, we have:

where 1/C is a positive value, and the corresponding output function of ELM is:

If the hidden layer feature mapping h⁢(𝒙) is unknown, a kernel function can be constructed to replace the calculation of 𝑯⁢𝑯T [13]:

Thus, Eq. (5) can be written as:

Then, we have:

For binary classification, the class label is determined by the two outputs of ELM through a competition mechanism. For the multi-class problem, the number of products should be the same as the number of categories. The output with the maximum value is the class to which the sample belongs.

2.2ELM based on Cholesky decomposition

In this paper, Cholesky decomposition is used to decrease the computation burden in the process of obtaining the output weights. This factorization method factorizes a matrix into the product of a triangular matrix and its conjugate transpose matrix.

The expression of the output weights rewrote as:

Assuming that

Then we have

Apparently,

For any 𝒂=[a1,a2,⋯,aN]≠𝟎, it holds that 𝒂T⋅𝑬⋅𝒂>0. Additionally, 𝑬⁢𝑬T⩾0. That is, 𝑬 is a positive definite matrix that can be decomposed in Cholesky form; i.e.,

where 𝑷 is an upper triangular matrix,

where

Thus, Eq. (11) can be written as:

That is,

Assuming 𝑴=𝑷-1⋅𝑻, we have

𝑴 can be obtained by solving the linear system of Eq. (16).

The output weights, 𝜷, can be calculated as:

Unlike conventional ELM calculations of the inverse matrix, this method directly obtains the output weights by iterative computations on Eqs (14), (17) and (18). These formulas involve only addition, subtraction, and square root operations; thus, the computational complexity is significantly reduced.

2.3Online update scheme

For the characteristic calculation based on Cholesky decomposition, it is convenient to extend the algorithm to the practical application of detecting epileptic seizures in routine clinical EEG recordings online.

When a new sample arrives, the matrix of the hidden layer can be written as:

Thus,

For the introduction of the kernel function, 𝑬𝑛𝑒𝑤 can be calculated through the following expression.

Through Eq. (14), it is evident there no need to recalculate the 1-N rows and the 1-N columns of 𝑬𝑛𝑒𝑤 during the computation based on Cholesky decomposition. Only a new lower triangular matrix must be calculated, to obtain 𝑬𝑛𝑒𝑤 complex

3.1Preprocessing

Original EEG signals pollute various interference signals, such as power line interference and electrooculograms (EOGs). To eliminate the effect of noise and to obtain cleaner EEG signals, preprocessing is necessary. Preprocessing steps commonly include filtering, data normalizing, artifact rejecting, etc. According to actual needs, reasonable measures can then be selected.

Waves with useful information are predominantly distributed in low-frequency regions. Therefore, for our application requirement, a six-order band-pass Chebyshev Type I filter with cutoff frequencies of 0.5 Hz and 40 Hz was designed to filter each extracted signal. This process comprised the first step of the analysis.

3.2Feature extraction

The purpose of the feature extraction process is to find compelling features to characterize the cognitive components. The extracted feature vectors of different tasks are expected to have distinct differences. Multiple elements are removed, from various EEG signals, including time domain features, wavelet packet energy, and entropy. The mathematical expressions of several featured used are:

The above listed time domain features are the most intuitive and straightforward ways to observe and analyze signals. However, for complex EEG signals with characteristics of nonlinearity, nonstationarity, and time variation, a single analysis method often cannot obtain a good effect. Therefore, we chose wavelet packet decomposition (WPD), which can efficiently locate signals in both time and frequency domains to extract the EEG features [18].

As shown in Fig. 3, through wavelet packet transformation, each epoch is decomposed into three levels. Eight sets of coefficients in the following frequency bands are obtained: 0.5–5 Hz, 5–10 Hz, 15–15 Hz, 15–20 Hz, 20–25 Hz, 25–30 Hz, and 35–40 Hz.

Figure 3.

The structure of WPD, Vij indicates the space expanded by the jth node of the ith layer WPD.

Then, the wavelet packet energy and entropy of each node are calculated as features of the EEG signals. They respectively indicate the strength and complexity of signals. The power of the EEG signal of a finite length is given by

where si⁢j represents the projection coefficients of a signal 𝒔i in an orthonormal basis. The energy feature of each epoch is:

The wavelet packet entropy is calculated according to Eq. (3.2), where the Shannon entropy is employed.

Thus, the entropy feature vector of each epoch is:

Consequently, the feature vector of each time is constructed as follows.

3.3Classification

Any function can be used as the kernel function of ELM as long as it aligns with Mercer’s theorem [19]. Several commonly used core functions exist, including the respective Gaussian kernel, polynomial kernel, perceptron kernel, radial basis function (RBF) kernel, wavelet kernel functions, among others. They each offer different advantages. Whether the selection of basic function is reasonable will directly affect the final classification result. A separate service often cannot achieve a satisfactory approximation effect. Accordingly, a combined kernel function, which is expected to obtain a better result, is constructed by adding different weights to different core roles in this study.

Here, the RBF kernel function, which has a stronger learning capability, and the polynomial kernel function, which offers a better generalization ability, are adopted to construct the combined core function. The expressions are as follows.

Thus, the combined primary function is:

Where η is the weight value.

After the initialization, differential evolution (DE) algorithm is adopted to obtain the optimal values of the three parameters (punishment factor C, kernel parameter α and weighting factor η) [20]. The mean square error of ELM is utilized as the fitness function.

4.1Data description

The Department of Epileptology, Bonn University, Germany [18] obtained the experimental data applied in this study; collecting from five healthy subjects and five epileptic patients. The complete dataset includes five sets (A-E), three of which are analyzed in this paper (A, D, and E). Details of the three datasets are listed in Table 1.

In each dataset, 100 single-channel EEGs of 23.6 s durations were recorded. The data sampling rate was 173.61 Hz. Thus, each EEG epoch had 4,096 sampling points.

4.2Results

In this section, the classification performance of the proposed algorithm is evaluated on the epileptic EEG datasets described above. A binary classifier was established to distinguish samples among healthy subjects (dataset A) and patients (dataset E). Additionally, the three-class problem among the three datasets was solved. The ten-fold cross-validation technique was used to reduce the bias of training and testing data. According to this technique, the dataset was divided into ten subsets [21]. To improve the dependability of this technique, the 10-fold cross-validation procedure was performed 10 times. Each time, one of the ten subsets was utilized as the testing dataset and the other 9 subsets were put together to form the training dataset. In particular, the data from test fold is not involved in the optimization procedure. All final results were averaged over the ten repetitions.

Tables 2 to 4 show the results. Specifically, Tables 2 and 3 compare the correctly classified percentage and the time required for training of different algorithms in the binary-class problem between health and seizures, and the three-class problem among health, seizure-free, and seizure activity. Moreover, to better inspect the performance, a confusion matrix is shown in Table 4.

4.3Discussion

Tables 2 and 3 show that the average recognition accuracies of our method in binary and three-class problems are both better than SVM, back-propagation (BP) neural network (BPNN), and the original ELM adopted. Because the Cholesky decomposition was adopted to simplify the calculation, our method was more time-efficient compared to SVM, BPNN and the method proposed in literature [11]. From the results, we can see that the proposed algorithm is suitable for the recognition of epileptic EEG patterns.

To apply the practice, there needs to be, no time-consuming operation in our approach. During the process of feature extraction, each epoch is decomposed into three levels by wavelet packet transform. More decomposition levels have no significant effects on the results, which was confirmed by experiments. Also, the classifier model can be quickly refreshed online if demanded (e.g., routine clinical applications). No need exists to retrain the entire network; only some parameters must be calculated (as described in Section 2.3). Through incremental recursion; the new training function can obtain new samples.

About future research, because there are various unpredictable interferences in EEG data collected in complex application environments, more pre-processing operations should be taken into account, such as automatic artifact rejection, to enhance the signal-to-noise ratio. Also, a more efficient feature extraction method, such as deep learning, can be a feasible way to improve the performance of the classifier. Additional experiments on practical applications are required, to address the remaining areas of improvement owing to the complexity of seizure recognition.