Epileptic seizure detection based on the kernel extreme learning machine

Abstract

This paper presents a pattern recognition model using multiple features and the kernel extreme learning machine (ELM), improving the accuracy of automatic epilepsy diagnosis. After simple preprocessing, temporal- and wavelet-based features are extracted from epileptic EEG signals. A combined kernel-function-based ELM approach is then proposed for feature classification. To further reduce the computation, Cholesky decomposition is introduced during the process of calculating the output weights. The experimental results show that the proposed method can achieve satisfactory accuracy with less computation time.

1. Introduction

Epilepsy is one of the most common chronic neurological disorders worldwide. The hyper-synchronized causes excessive electrical discharges in a group of neurons. Approximately one in every 100 people is afflicted by it [1]. EEG recorded by electrodes placed on the scalp is an efficient method for checking electrical activity in the brain. Conventional diagnosis of neurological disorders based on EEG signals relies on neurologists to visually check the recordings; it is experience-dependent and time-consuming. Therefore, to reduce the artificial workload and improve the efficiency and accuracy of the diagnosis, the application of signal processing and machine learning methods for epileptic EEG signal analysis is valuable.

The process of automatic epileptic seizure analysis is a pattern recognition system that can distinguish seizures from standard EEG signals [2]. In recent years, many problems have been solved, with relevant methods. Suitable features are essential for correct detection of epileptic patterns. To date, temporal [3], frequency [4], wavelet [5], spatial [6], nonlinear dynamic, and deterministic chaos features, as well as multi-feature fusion strategies [7], have been used to describe EEG signals [8, 9]. In addition, many effective classification algorithms have been used in this area, such as support vector machine (SVM) [10, 11], artificial neural networks (ANN) [12], etc. Some of the methods have achieved satisfactory performance. Especially, work [11] has obtained extraordinary recognition accuracy.

However, owing to the poor signal-to-noise ratio (SNR) of raw EEG signals in practical application, the problems of improving recognition performance and identification efficiency still exist. The limits existing in conventional methods including overlong training time, unsatisfactory accuracy or limited generalization performance for clinical application make the automatic detection of epileptic seizure still challenged. Compared with other algorithms, the ELM requires less training time while maintaining satisfactory classification accuracy. Problems existing in traditional feedforward neural network learning algorithms such as local minima and various training parameters are avoided in ELM. Furthermore, it has a higher generalization performance. Owing to its superior property, ELM has thus been applied to EEG signal feature classification based on diverse feature extraction methods and has achieved impressive performance in literatures [5, 13, 14] for instance.

To better balance the conflict between efficiency and accuracy of pattern recognition, this paper focuses primarily on the classification process for seizure recognition. In the proposed algorithm, several features are combined comprehensively to represent the characteristics of EEG for epilepsy. Then a combined kernel function is introduced into ELM. Moreover, Cholesky decomposition is employed to reduce the calculation burden and the effectiveness of the algorithm, validated through experiments.

2. Kernel elm based on cholesky decompositon

Unlike conventional single hidden layer feed forward neural networks (SLFNs), the parameters of hidden layers in ELM are randomly given at the beginning of the training process and fixed it instead of the complicated iterative calculation process. Then, the least-squares method is used to obtain the output weights [15].

Otherwise, ELM tends to minimize the training error and the norm of production weights. The introduction of kernel function gives the algorithm more stability and better generalization ability [16]. All the above attributes give the kernel ELM higher efficiency and better performance [17].

2.1 Kernel ELM

Suppose there are $N$ arbitrary samples $({𝒙}_i,{𝒕}_i)$ , where ${𝒙}_i={[x_{i1},x_{i2},\mathrm{\cdots },x_{in}]}^T\in {𝑹}^n$ , ${𝒕}_i={[t_{i1},t_{i2},\mathrm{\cdots },t_{im}]}^T\in {𝑹}^m$ . Thus, standard SLFNs with $L$ hidden nodes can be mathematically expressed:

$f_L(𝒙)={\displaystyle \sum _{i=1}^L}{𝜷}_i{g}_i(𝒙)={\displaystyle \sum _{i=1}^L}{𝜷}_{i}G({𝒂}_i,b_i,{𝒙}_j),j=1,\mathrm{\cdots },N$(1)

Figure 1.

Single hidden layer feedforward network.

Where ${𝜷}_i={[{\beta }_{i1},{\beta }_{i2},\mathrm{\cdots },{\beta }_{im}]}^T$ is the weight vector connecting the $i$ th hidden node and the output nodes. ${𝒂}_i={[a_{i1},a_{i2},\mathrm{\cdots },a_{in}]}^T$ Is the weight vector connecting the $i$ th hidden node and the input nodes, and $b_i$ is the threshold of the $i$ th hidden node. $g_i$ denotes the output function $G({𝒂}_i,b_i,𝒙)$ of the $i$ th hidden node (cf. Fig. 1). where

$𝑯={\left[\begin{array}{cccccccccccccccccccc}\hfill h({𝒙}_1)\hfill & & & & & & & & & & & & & & & & & & & \\ \hfill \mathrm{⋮}\hfill & & & & & & & & & & & & & & & & & & & \\ \hfill h({𝒙}_N)\hfill & & & & & & & & & & & & & & & & & & & \end{array}\right]}_{N\times L}={\left[\begin{array}{cccccccccccccccccccc}\hfill G({𝒂}_1,b_1,{𝒙}_1)\hfill & \hfill \mathrm{\cdots }\hfill & \hfill G({𝒂}_L,b_L,{𝒙}_1)\hfill & & & & & & & & & & & & & & & & & \\ \hfill \mathrm{⋮}\hfill & \hfill \mathrm{\cdots }\hfill & \hfill \mathrm{⋮}\hfill & & & & & & & & & & & & & & & & & \\ \hfill G({𝒂}_1,b_1,{𝒙}_N)\hfill & \hfill \mathrm{\cdots }\hfill & \hfill G({𝒂}_L,b_L,{𝒙}_N)\hfill & & & & & & & & & & & & & & & & & \end{array}\right]}_{N\times L},$ $𝜷={\left[\begin{array}{cccccccccccccccccccc}\hfill {𝜷}_1^T\hfill & & & & & & & & & & & & & & & & & & & \\ \hfill \mathrm{⋮}\hfill & & & & & & & & & & & & & & & & & & & \\ \hfill {𝜷}_L^T\hfill & & & & & & & & & & & & & & & & & & & \end{array}\right]}_{L\times m},T={\left[\begin{array}{cccccccccccccccccccc}\hfill t_1^T\hfill & & & & & & & & & & & & & & & & & & & \\ \hfill \mathrm{⋮}\hfill & & & & & & & & & & & & & & & & & & & \\ \hfill t_N^T\hfill & & & & & & & & & & & & & & & & & & & \end{array}\right]}_{N\times m}$

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

$\parallel 𝑯\widehat{𝜷}-𝑻\parallel =\underset{𝜷}{\text{min}}\parallel 𝑯𝜷-𝑻\parallel$(2)

The unique solution is

$\widehat{𝜷}={𝑯}^{+}𝑻$(3)

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

Randomly assign input weights $a_i$ and biases $b_i$ according to some continuous probability density function;

Calculate the hidden layer output matrix $𝑯$ ;

Obtain the output weights according to Eq. (3).

If $H{H}^T$ is nonsingular, to improve the stability of ELM, we have:

$𝜷={𝑯}^T{\left({\displaystyle \frac{1}{C}}+𝑯{𝑯}^T\right)}^{-1}𝑻$(4)

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

$h(𝒙)𝜷=h(𝒙){𝑯}^T{\left({\displaystyle \frac{1}{C}}+𝑯{𝑯}^T\right)}^{-1}𝑻$(5)

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

${𝛀}_{ELM}=𝑯{𝑯}^T=\left[\begin{array}{cccccccccccccccccccc}\hfill K({𝒙}_{1,}{𝒙}_1)\hfill & \hfill \mathrm{\cdots }\hfill & \hfill K({𝒙}_{1,}{𝒙}_N)\hfill & & & & & & & & & & & & & & & & & \\ \hfill \mathrm{⋮}\hfill & \hfill \mathrm{\cdots }\hfill & \hfill \mathrm{⋮}\hfill & & & & & & & & & & & & & & & & & \\ \hfill K({𝒙}_{N,}{𝒙}_1)\hfill & \hfill \mathrm{\cdots }\hfill & \hfill K({𝒙}_{N,}{𝒙}_N)\hfill & & & & & & & & & & & & & & & & & \end{array}\right]$(6)

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

$h(𝒙)𝜷={\left[\begin{array}{cccccccccccccccccccc}\hfill K(𝒙,{𝒙}_1)\hfill & & & & & & & & & & & & & & & & & & & \\ \hfill \mathrm{⋮}\hfill & & & & & & & & & & & & & & & & & & & \\ \hfill K(𝒙,{𝒙}_N)\hfill & & & & & & & & & & & & & & & & & & & \end{array}\right]}^T{\left({\displaystyle \frac{1}{C}}+{𝛀}_{\text{𝐸𝐿𝑀}}\right)}^{-1}𝑻$(7)

Then, we have:

$f(𝒙)={\left[\begin{array}{cccccccccccccccccccc}\hfill K(𝒙,{𝒙}_1)\hfill & & & & & & & & & & & & & & & & & & & \\ \hfill \mathrm{⋮}\hfill & & & & & & & & & & & & & & & & & & & \\ \hfill K(𝒙,{𝒙}_N)\hfill & & & & & & & & & & & & & & & & & & & \end{array}\right]}^T{\left({\displaystyle \frac{1}{C}}+{𝛀}_{\text{𝐸𝐿𝑀}}\right)}^{-1}𝑻$(8)

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.2 ELM 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:

$𝜷={𝑯}^T{\left({\displaystyle \frac{1}{C}}+𝑯{𝑯}^T\right)}^{-1}𝑻$(9)

Assuming that

$E=𝑯{𝑯}^T+{\displaystyle \frac{𝑰}{C}}$(10)

Then we have

$𝑬\cdot 𝜷=𝑻$(11)

Apparently,

${𝑬}^T={\left(𝑯{𝑯}^T+{\displaystyle \frac{𝑰}{C}}\right)}^T=𝑯{𝑯}^T+{\displaystyle \frac{𝑰}{C}}=𝑬$(12)

For any $𝒂=[a_1,a_2,\mathrm{\cdots },a_N]\ne \mathrm{𝟎}$ , it holds that ${𝒂}^T\cdot 𝑬\cdot 𝒂>0$ . Additionally, $𝑬{𝑬}^T⩾0$ . That is, $𝑬$ is a positive definite matrix that can be decomposed in Cholesky form; i.e.,

$𝑬=𝑷\cdot {𝑷}^T$(13)

where $𝑷$ is an upper triangular matrix,

$𝑷=\left[\begin{array}{cccccccccccccccccccc}\hfill p_{11}\hfill & \hfill 0\hfill & \hfill \mathrm{\cdots }\hfill & \hfill 0\hfill & & & & & & & & & & & & & & & & \\ \hfill p_{21}\hfill & \hfill p_{22}\hfill & \hfill \mathrm{\cdots }\hfill & \hfill 0\hfill & & & & & & & & & & & & & & & & \\ \hfill \mathrm{⋮}\hfill & \hfill \mathrm{⋮}\hfill & \hfill \mathrm{\ddots }\hfill & \hfill \mathrm{⋮}\hfill & & & & & & & & & & & & & & & & \\ \hfill p_{N1}\hfill & \hfill p_{N2}\hfill & \hfill \mathrm{\cdots }\hfill & \hfill p_{NN}\hfill & & & & & & & & & & & & & & & & \end{array}\right]$

where

$p_{ij}=\{\begin{array}{cccccccccccccccccccc}\sqrt{e_{ij}},\hfill & i=j=1\hfill & & & & & & & & & & & & & & & & & & \\ \sqrt{e_{ii}}-{\displaystyle \sum _{n=1}^{i-1}}{p}_{in}^2,\hfill & i=j>1\hfill & & & & & & & & & & & & & & & & & & \\ (e_{ii}-{\displaystyle \sum _{n=1}^{i-1}}{p}_{in}{p}_{jn})/p_{jj},\hfill & i>j\hfill & & & & & & & & & & & & & & & & & & \end{array}$(14)

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

$𝑷\cdot {𝑷}^T\cdot 𝜷=𝑻$

That is,

${𝑷}^T\cdot 𝜷={𝑷}^{-1}\cdot 𝑻$(15)

Assuming $𝑴={𝑷}^{-1}\cdot 𝑻$ , we have

$𝑷\cdot 𝑴=𝑻$(16)

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

$m_{ij}=\{\begin{array}{cccccccccccccccccccc}{𝒕}_i/p_{ii},\hfill & i=1\hfill & & & & & & & & & & & & & & & & & & \\ ({𝒕}_i-{\displaystyle \sum _{j=1}^{i-1}}{p}_{ii}{𝒎}_j)/p_{ii},\hfill & i>1\hfill & & & & & & & & & & & & & & & & & & \end{array}$(17)

The output weights, $𝜷$ , can be calculated as:

${\beta }_{ij}=\{\begin{array}{cccccccccccccccccccc}{𝒎}_i/p_{ii},\hfill & i=N\hfill & & & & & & & & & & & & & & & & & & \\ ({𝒎}_i-{\displaystyle \sum _{j=i+1}^N}{p}_{ji}{𝜷}_j)/p_{ii},\hfill & i<N\hfill & & & & & & & & & & & & & & & & & & \end{array}$(18)

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.3 Online 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:

${𝑯}_{\text{𝑛𝑒𝑤}}=\left[\begin{array}{c}\hfill 𝑯\hfill \\ \hfill h({𝒙}_{N+1})\hfill \end{array}\right],{𝑻}_{\text{𝑛𝑒𝑤}}=\left[\begin{array}{c}\hfill 𝑻\hfill \\ \hfill {𝒕}_{N+1}\hfill \end{array}\right]$

Thus,

${𝑬}_{\text{𝑛𝑒𝑤}}=\left({𝑯}_{\text{𝑛𝑒𝑤}}{𝑯}_{\text{𝑛𝑒𝑤}}^T+{\displaystyle \frac{{𝑰}_{N+1}}{C}}\right)=\left[\begin{array}{cccccccccccccccccccc}\hfill 𝑬\hfill & \hfill (h_{N+1}\cdot 𝑯)\hfill & & & & & & & & & & & & & & & & & & \\ \hfill {(h_{N+1}\cdot 𝑯)}^T\hfill & \hfill h_{N+1}\cdot h_{N+1}^T+1/C\hfill & & & & & & & & & & & & & & & & & & \end{array}\right]$(19)

For the introduction of the kernel function, ${𝑬}_{\text{𝑛𝑒𝑤}}$ can be calculated through the following expression.

${𝑬}_{\text{𝑛𝑒𝑤}}=\left[\begin{array}{cccccccccccccccccccc}\hfill 𝑬\hfill & \hfill k({𝒙}_{n+1},𝑿)\hfill & & & & & & & & & & & & & & & & & & \\ \hfill {(k({𝒙}_{n+1},𝑿))}^T\hfill & \hfill k{({𝒙}_{n+1},{𝒙}_{n+1})}^T+1/C\hfill & & & & & & & & & & & & & & & & & & \end{array}\right]$(20)

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

Calculate, ${𝑷}_{\text{𝑛𝑒𝑤}}=[𝑷,p_{N+1}]$ , where $p_{N+1}$ can be obtained by Eq. (14);

Calculate, ${𝑴}_{\text{𝑛𝑒𝑤}}=[𝑴,M_{N+1}]$ , where $M_{N+1}$ can be obtained by Eq. (17);

Recalculate the output weights according to Eq. (18).

3. Method

The flowchart of the recognition system mainly includes data collection, preprocessing, feature extraction, and classification (cf. Fig. 2).

Figure 2.

Flowchart of epileptic EEG recognition.

3.1 Preprocessing

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.2 Feature 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:

Crest Factor

$Cf=\text{max}\left({𝒔}_i\right)/\sqrt{{\displaystyle \frac{s_{i1}^2+s_{i2}^2+\mathrm{\cdots }+s_{in}^2}{n}}}$

Kurtosis

$Ku={\displaystyle \sum _{j=1}^n}{\left(s_{ij}-{\overline{𝒔}}_i\right)}^4/(n-1)\text{std}\left({𝒔}_i\right)$

Impulse Factor

$If=\text{max}\left({𝒔}_i\right)/\overline{\left|{𝒔}_i\right|}$

Signal Factor

$Sf=\sqrt{{\displaystyle \frac{s_{i1}^2+s_{i2}^2+\mathrm{\cdots }+s_{in}^2}{n}}}/\overline{\left|{𝒔}_i\right|}$

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, $V_i^j$ 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

$En\left({𝒔}_i\right)={\displaystyle \sum _i}{s}_{ij}^2$ (21) $i=1,\mathrm{\cdots },N,j=1,\mathrm{\cdots },n$

where $s_{ij}$ represents the projection coefficients of a signal ${𝒔}_i$ in an orthonormal basis. The energy feature of each epoch is:

${\text{𝐄𝐧𝐱}}_i=[E{n}_i(s_{i3}^0),E{n}_i(s_{i3}^1),\mathrm{\cdots },E{n}_i(s_{i3}^7)],i=1,2,\mathrm{\cdots },N$(22)

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

$\text{𝐸𝑛𝑡}({𝒔}_i)=-{\displaystyle \sum _i}{s}_{ij}^2\log (s_{ij}^2)$ (23) $i=1,\mathrm{\cdots },N,j=1,\mathrm{\cdots },n$

Thus, the entropy feature vector of each epoch is:

${\text{𝐄𝐧𝐭𝐱}}_i=[{\text{𝐸𝑛𝑡}}_i(s_{i3}^0),{\text{𝐸𝑛𝑡}}_i(s_{i3}^1),\mathrm{\cdots },{\text{𝐸𝑛𝑡}}_i(s_{i3}^7)],i=1,2,\mathrm{\cdots }N$(24)

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

${𝒙}_i=[C{f}_i,K{u}_i,I{f}_i,S{f}_i,{\text{𝐄𝐧𝐱}}_i,{\text{𝐄𝐧𝐭𝐱}}_i],i=1,2,\mathrm{\cdots },N$(25)

3.3 Classification

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.

RBF kernel function:

$k(𝒙,{𝒙}_i)=\exp (-{\parallel 𝒙-{𝒙}_i\parallel }^2/\alpha )$

Polynomial kernel function:

$k(𝒙,{𝒙}_i)={(𝒙{𝒙}_i+1)}^{\gamma }$

Thus, the combined primary function is:

$k(𝒙,{𝒙}_i)=\eta {(𝒙{𝒙}_i+1)}^{\gamma }+(1-\eta )\exp (-{\parallel 𝒙-{𝒙}_i\parallel }^2/\alpha )$(26)

Where $\eta$ 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 $\alpha$ and weighting factor $\eta$ ) [20]. The mean square error of ELM is utilized as the fitness function.

Table 1

Description of the three datasets analyzed

Datasets	Subject	Electrode type	Subject’s state
Set A	Five healthy subjects	Surface	Normal
Set D	Five patients	Intracranial	Seizure-free
Set E	Five patients	Intracranial	Seizure activity

4. Experiments

4.1 Data 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.2 Results

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.

Table 2

Comparison of different algorithms for binary problem

Method	Average classification accuracy (%)	Average training time (s)
BPNN	92.2	69.82
SVM	90.3	23.15
Original ELM	94.2	0.145
Work [11]	100	59.20
Work [5]	94.8	0.980
Proposed	97.7	0.104

Table 3

Comparison of different algorithms for a three-class problem

Method	Average classification accuracy (%)	Average training time (s)
BPNN	90.9	80.97
SVM	88.4	34.50
Original ELM	93.6	0.238
Work [11]	98.2	75.70
Work [7]	96.0	80.71
Proposed	96.5	0.157

Table 4

Confusion matrix of the proposed method

	Set A	Set D	Set E
Set A	100	4	1
Set D	0	94	3
Set E	0	2	96

4.3 Discussion

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.

5. Conclusions

In this paper, we proposed an ELM kernel algorithm by introducing a combined kernel function to address the problem of seizure recognition. By employing Cholesky decomposition and its calculation process, which involves only arithmetic, the calculation efficiency of the proposed method is further improved.

Among different classifiers, a comparative study was conducted, to illustrate effectiveness in our approach. The results show that our method achieves better recognition accuracy with considerably less training time. The overall implementation of the method is easy to understand, and the computation burden is low.

Conflict of interest

None to report.