A Robust Automatic Epilepsy Seizure Detection Algorithm Based on Interpretable Features and Machine Learning

Abstract

Epilepsy, as a serious neurological disorder, can be detected by analyzing the brain signals produced by neurons. Electroencephalogram (EEG) signals are the most important data source for monitoring these brain signals. However, these complex, noisy, nonlinear and nonstationary signals make detecting seizures become a challenging task. Feature-based seizure detection algorithms have become a dominant approach for automatic seizure detection. This study presents an algorithm for automatic seizure detection based on novel features with clinical and statistical significance. Our algorithms achieved the best results on two benchmark datasets, outperforming traditional feature-based methods and state-of-the-art deep learning algorithms. Accuracy exceeded 99.99% on both benchmark public datasets, with the 100% correct detection of all seizures on the second one. Due to the interpretability and robustness of our algorithm, combined with its minimal computational resource requirements and time consumption, it exhibited substantial potential value in the realm of clinical application. The coefficients of variation of datasets proposed by us makes the algorithm data-specific and can give theoretical guidance on the selection of appropriate random spectral features for different datasets. This will broaden the applicability scenario of our feature-based approach.

1. Introduction

1.1. Background

Epilepsy, as a common neurological disorder, affects over 70 million patients worldwide [1]. It usually leads to 2–3 times more premature mortality than ordinary, burdening families and society [2]. As an unstable situation in epilepsy patients due to the excessive discharge of brain cells, seizures can cause a range of physiological symptoms, such as partial or entire body reactions that are uncontrollable to the complete loss of consciousness. It can even lead to severe injury or death. Fortunately, the World Health Organization (WHO) reports that 70% of people with epilepsy can avoid seizures if they receive timely treatment [3]. Electroencephalogram (EEG) signals are the most important data source for monitoring these brain signals. Therefore, it is essential to identify seizure events from the EEG and diagnose patients with epilepsy early. EEG signals are obtained by amplifying and recording the spontaneous biological potentials of the brain from the scalp through a precision electronic instrument. They are spontaneous and rhythmic electrical activities of brain cells recorded through electrodes. Epileptic disorders are diagnosed by EEG signals on the basis that they capture and record the abnormal electrical activities in the cerebral cortex. The clinical analysis and interpretation of EEG are time-consuming and put a heavy burden on neurologists [4,5]. To reduce the workload of neurologists by giving more timely diagnostic results, many studies focus on the automatic detection of seizures [6]. Data-driven methods such as machine learning and deep learning classifiers are able to classify EEG data and detect epileptic seizures without compromising performance, while revealing underlying patterns [7]. Therefore, relevant algorithms provide enlightening insights into the automatic detection of seizures.

There are four main stages of epilepsy, including interictal state, preictal state, ictal state, and postictal state [8]. The interictal state is the period between two adjacent seizure events. The ictal state is the period during which a seizure occurs. Preictal and postictal states are the states before and after a seizure event, respectively. Broadly speaking, both the preictal and postictal states can be regarded as part of the interictal state. The goal of seizure detection is to distinguish between ictal and interictal state, to assist doctors in timely diagnosis and treatment, and to achieve the early detection of epilepsy patients. At present, data-driven methods in seizure recognition research fall into two categories: end-to-end deep learning methods and feature extraction combined with classifiers. The end-to-end deep learning method can automatically extract features, reducing human intervention. At the same time, it has a good effect on nonlinear and non-stationary signal processing such as EEG, but the features extracted by the deep learning methods have poor interpretation, so it is often not trusted by neurologists and difficult to be used in practice. However, the effectiveness of methods based on the combination of the feature extraction and classifier mainly depends on the feature extraction process, and therefore, may not produce satisfactory results like deep learning methods. Therefore, a robust automatic seizure detection algorithm based on interpretable features and machine learning was proposed in this paper, combining the interpretable features and the reliability of the results. The main contributions of this paper are as follows:

• Proposed an algorithm based on the combination of feature extraction and machine learning classifier to detect seizure activities, and the results are optimal compared with the state-of-the-art algorithms on two EEG datasets;
• Based on the experience of clinical experts and the theory of random matrix, the characteristics of epileptic wave numbers and random spectral features were proposed to distinguish ictal and interictal state;
• Presented an unsupervised detection algorithm tailored to these epileptiform wave features;
• Defined the coefficient of variation for datasets. It provides a theoretical basis for the adaptive extraction of random spectral features from different datasets.

1.2. Related Works

1.2.1. End-to-End Deep Learning Methods

End-to-end deep learning methods are used to directly and automatically carry out the feature extraction and classification of epilepsy EEG data. According to the framework used, they can be roughly divided into networks based on convolutional neural network (CNN), based on full convolutional network (FCN), based on random convolutional kernel transformation (Rocket), and based on Transformer.

One of the most well-known deep learning techniques is convolutional neural network (CNN), which can extract features by using various filters in the convolutional layer, pooling layer, normalized layer, and fully connected layer [9]. A 1D convolutional neural network-long short-term memory (1D CNN-LSTM) model proposed by Xu et al. [10] combined CNN with the long-short term memory network (LSTM), which can effectively capture and model the long-term dependence relationship to achieve a recognition accuracy of 99.39% in the seizure recognition task. Ismail et al. [11] used Resnet networks to achieve 50 wins out of 85 time classification tasks for UCR and UEA. Resnet adds residual blocks on the basis of traditional CNN, so it can avoid the problem of gradient disappearance during deep network training. Although this kind of deep neural network has excellent classification accuracy, it often has a huge training cost. Eldeleal et al. [12] introduced a new time series lightweight adaptive network (TSLANet) as a general convolutional baseline for different time series tasks. It has an adaptive spectrum module that can use Fourier analysis to enhance feature representation, and an interactive convolution block that enhances the model’s ability to decode complex time patterns, achieving an average classification accuracy of 85.90% across 116 datasets.

In FCN architecture, after the convolutional blocks, the features are input into a global average pooling layer [13] instead of a fully connected layer, which greatly reduces the number of model parameters and is conducive to improving model training efficiency. FCNs are essentially convolutional networks with no local pooling layers. This assures that a time series’ length remains constant throughout the convolutions, which makes it more conducive to fully extracting the information in the time series compared with CNN. The FCN-based temporal classification framework proposed by Wang et al. [14] had an average error per class (MPCE) of 0.023 on 44 benchmark datasets, which is the lowest among all deep learning methods.

Another popular framework with lower computational costs is the random convolutional kernel transform (Rocket). Rocket uses random convolutional kernel to first transform the time series data, and then trains a linear classifier using the transformed features [15]. As the kernel weights are not learned, the computational cost of computing the convolutions is low. However, the original Rocket is a univariate classification model that has not been generalized to multivariate temporal classification tasks. Another limitation is that, when using fixed random convolution kernels, for very large datasets, learning may become “saturated” at some point, making it impossible to continue improving accuracy. To solve these problems, variants such as MiniRocket [16] and MultiRocket [17] were developed. Recently, Dempster et al. [18] proposed Hydra, an algorithm that combines traditional dictionary methods with Rocket, which significantly improved the classification accuracy of Rocket and its variants. Hydra can be trained and tested on 106 datasets in approximately 31 min using a single CPU core.

In recent years, Transformer has attracted great interest from researchers due to its superior performance in natural language processing and speech processing. In particular, Transformer has shown good modeling capabilities in terms of the random dependencies and interactions in sequence data, which is attractive in the field of time series modeling [19]. Baurzhan [20] obtained the optimal result of Parkinson’s gait prediction and recognition by combining transformer encoder layer and two LSTM layers. However, the efficacy of Transformer in time series has also been questioned. For example, the alignment invariance of self-attention mechanisms in Transformer inevitably leads to the loss of time information [21]. In addition, self-attention mechanisms typically require a lot of computational resources, and they are prone to overfitting when trained on smaller datasets.

In summary, although end-to-end deep learning methods generally show superior performance, training models can be time-consuming and it is easier to overfit in mini data [22]. Furthermore, these models are essentially uninterpretable “black boxes” and are often considered unacceptable by neurologists in clinical healthcare settings [23]. Methods that use a combination of feature extraction and machine learning classifiers tend to have lower computational costs and more interpretive results. We combine the recent deep learning models mentioned above, FCN (Wang et al., 2019 [14]), Resnet (Fawaz et al., 2019 [11]), TSLANet (Eldeleal et al., 2024 [12]), Hydra (Dempster et al., 2024 [18]), and Transformer + LSTM (Ignacio, 2023 [20]), which are used as benchmarks to compare with our approach to evaluate the performance of our model.

1.2.2. Feature-Based Methods

The second category methods extract corresponding features according to expert experience or domain knowledge, and then use machine learning classification algorithm to identify seizures. Compared with end-to-end ones, feature engineering makes the results more interpretable and acceptable to neurologists, and the lower computational cost also makes it more potential for application. Traditional EEG analysis extracted features in the time domain, frequency domain, time–frequency domain, space domain, or nonlinear features as input variables, then used a variety of classifiers to classify the EEG data [24,25,26,27,28,29,30,31,32,33,34,35,36]. These features can be categorized based on their interpretation or the domain from which the features are calculated [37]. Chen et al. [38] proposed an algorithm for seizure detection based on time–frequency domain and nonlinear features in 2023. They extracted mixed features such as approximate entropy (ApEn), fuzzy entropy (FuzzyEn), sample entropy (SampEn), and the standard deviation (STD) of EEG signals in the subbands decomposed by discrete wavelet transform (DWT). The classification with CNN network after feature selection through random forest algorithm achieved 99.9% accuracy on the Bonn dataset and 100% accuracy on India dataset. Mahshidh et al. [39] similarly extracted the entropy features of the sub-frequency bands after DWT processing, then used the ensemble algorithm to classify after feature selection, and the seizure detection accuracy rate reached 99.5%. Although the accuracy of these algorithms is high, the calculation of the above time–frequency domain and nonlinear features are often complicated.

Velichko et al. [40] proposed a new entropy NNetEn, which reflected the reservoir information of the time series in the neural network, and demonstrated the practical application of the feature in distinguishing EEG signals from healthy people and patients with Alzheimer’s disease. However, the entropy characteristics depend on the parameters, which need to be adjusted for different application scenarios. And this feature also suffers from the limitation of computational complexity.

Although the features above were primarily constructed from data characteristics observed during the interictal and ictal phases, they still lack convincing clinical background knowledge. Epileptiform discharges are essential markers for physicians to diagnose the presence of seizures in patients. Slimen and Seddik [41] used machine learning methods and took the epileptiform wave number as an indicator to classify the EEG signals. However, they only considered spike waves, although slow and sharp waves also constitute significant indicators for epileptic events that cannot be overlooked.

Random spectral feature extraction, a prevalent technique based on random projection in high dimensional space, has recently gained considerable interest. Liao and Couillet [42] employed the random matrix theory to examine the impact of various nonlinear mappings on feature extraction. They suggested a selection criterion for nonlinear mapping functions based on variance and mean differences in datasets. Their recommendations provide valuable theoretical guidance for selecting appropriate functions for datasets characterized by varying variances and mean differences, albeit without quantitative measurements of the differences in mean and variance.

Considering the classification accuracy and computational complexity, we introduced wave number features and random spectral features with high single-feature discriminability among different categories, which achieve a better classification accuracy than deep neural networks by combining with traditional EEG features.

2. Materials and Methods

2.1. Datasets

This work used two EEG datasets. The first one was collected from Neurology & Sleep Centre, Hauz Khas, New Delhi [43]. After being acquired at a sampling rate of 200 Hz using Grass Telefactor Comet AS40 Amplification System, the signals were filtered between 0.5 and 70 Hz and then segmented into pre-ictal, interictal, and ictal stages by several experts. There are three folders named according to the seizure phase, and each folder contains 50 mat files. Each mat file consists of a sample of 1024 non-overlapping EEG time series data, corresponding to a time duration of 5.12 s. In this work, we used all the interictal and ictal stages data for research.

To test the robustness of our proposed method, it has been tested on another single-channel dataset which was collected at the University of Bonn, Germany [44]. This database has also become a benchmark for seizure detection. As shown in Figure 1, it consists of five subsets (Set A~Set E). Each set of EEG recordings consists of 100 segments. And each segment represented by 4097 samples, with a single channel of 23.6 s and a sampling rate of 173.61 Hz. The five subgraphs in Figure 1 visually show a segment of the five subsets, respectively, in which the first 4000 sample points of modified segments were selected. These segments were selected from the continuous EEG recordings after visual inspection for artifacts. EEG data were obtained from five patients who achieved complete seizure control after removal of one tumor. Among them, sets D and E contained EEG segments of interictal and ictal periods, respectively. Therefore, data from sets D and E were used in this study.

2.2. Methods

The process of feature-based automatic seizure detection based on EEG signals can be divided into data pre-processing, feature extraction, and seizure detection with classification algorithms [45]. Figure 2 presents the normal process of feature-based EEG signal seizure detection algorithm.

EEG as biomedical signals were contaminated by noise and artifacts, and therefore require a pre-processing step to clean up the data [46]. Noise and artifacts arise from physiological activities such as eye blinks, eye movements, and muscle contractions. In addition, the detection equipment itself can cause interference [47]. These can usually be removed using filtering, wavelet analysis, canonical correlation analysis, etc. Feature extraction, as a core component of any pattern recognition system, typically employs mathematically driven algorithms to extract relevant information from the raw data to better characterize patterns at any given point in time [48]. In many cases, integrating the feature extraction components into a pattern recognition system tends to accurately better distinguish between various patterns of interest and produce results faster than using raw data directly. Therefore, feature extraction is considered essential for the development of an efficient intelligent system for epileptic seizure detection. After extracting the relevant features, machine learning classifiers are applied to give detection results.

Datasets used in this article have already undergone pre-processing to remove artifacts, so no further actions are necessary. And the study of pre-processing methods is beyond the scope of this paper. We next introduce features that are interpretable and robust in terms of clinical significance and statistical significance, respectively.

2.2.1. Wave Numbers

Epileptiform discharges with different changes are generally seizures periods’ characteristics. The appearance of spike and sharp waves are some especially key characteristics of epileptic periods, and seizures are also often accompanied by the appearance of slow waves [49]. In view of these facts, we proposed three novel features—spike wave numbers, sharp wave numbers, and slow wave numbers. For the sake of the calculation of wave numbers, we also proposed an unsupervised wave detection algorithm, which can recognize spike, sharp, and slow waves without relevant labels.

As described in Algorithm 1, the unsupervised detective algorithm starts with searching for all maximum local peaks in each sample window. A function named find_peaks in python 3.8 package SciPy [50] was used to search peaks (local maxima) based on the simple value comparison of neighboring samples and returns those peaks whose properties match optionally specified conditions (minimum and/or maximum) for their height, prominence, width, threshold, and distance from each other. According to the definition, spikes are the most basic paroxysmal EEG activity, with a duration of 20~70 ms (14.5~50 Hz), and an amplitude of more than 50 μV. As for sharp wave, the time limit is 70~200 ms (5~14.7 Hz), and wave amplitude is approximately between 100 and 200 μV. And, Slow wave’s frequency is lower than 8 Hz [51]. Then, we choose wave peaks in maximum local peaks by judging whether they meet the definition. We will then delete peaks which are local maxima but not real wave peaks. Finally, we sum the wave numbers of each class of waves in each sample to obtain the features that will be input to the machine learning method.

Algorithm 1. Unsupervised wave detective algorithm

Input:   EEG sample X;   Wave-path: ${\{x_{t,i},y_{t,i}\}}_{i=\mathrm{1,2},3}$;   Amplitude range: :${\{x_{a,i},y_{a,i}\}}_{i=\mathrm{1,2},3}$ Output:   Wave number: $n_i$, i = 1,2,3;   Wave center position: ${\mathit{C}}_i=\{c_{1i},c_{2i},\dots ,c_{n_{i}i}\}$, i = 1,2,3; For i = 1 to 3 do Use the function find-peaks to find all local maximum points in X as a set of candidate    peaks for class $i^{th}$ epileptic wave: $P_i=\left\{p_{1i},p_{2i},\dots ,p_{m_{i}i}\right\}\leftarrow$find-peaks(X,i) Select candidate peaks where wave-path and amplitude fall within the range :$\{x_{t,i},y_{t,i}\}$ and $\{x_{a,i},y_{a,i}\}$:    ${\mathit{C}}_i=\{c_{1i},c_{2i},\dots ,c_{n_{i}i}\}\stackrel{wath-path\in \left\{x_{t,i},y_{t,i}\right\},\mathrm{a}\mathrm{m}\mathrm{p}\mathrm{l}\mathrm{i}\mathrm{t}\mathrm{u}\mathrm{d}\mathrm{e}\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{g}\mathrm{e}\in \{x_{a,i},y_{a,i}\}}{\leftarrow }{P}_i$ End for Return $n_i, {\mathit{C}}_i$

where i = 1, 2, and 3, these correspond to spike, sharp, and slow waves, respectively. The wave-path and amplitude range are determined by their definitions.

2.2.2. Random Feature

Random feature map is a method of extracting data features by random projection in high-dimensional space. More precisely, for the data x, we generate a random matrix W, whose elements are independently and identically distributed (i.i.d.) random variables, and choose a nonlinear activation function $\sigma (·)$; then the random character of the data is given by $\sigma \left(\mathbf{W}\mathbf{x}\right)$, where the nonlinear function acts on each element of the vector. This is similar to a MLP layer of a neural network, where W is the weight of that layer.

The mathematical expression for random feature extraction is as follows:

Let there be T independent data vectors obeying a Gaussian mixture model with K categories ${Class}_1$,…${Class}_K$, and the data vectors have p dimensions, i.e., ${\mathit{x}}_1,\dots ,{\mathit{x}}_T\in {\mathit{R}}^p$. Given a certain category a, there are $T_a$ samples under the category, for a data vector ${\mathit{x}}_i$ from ${Class}_a$ has:

$${\mathit{x}}_i=\frac{{\mathit{\mu }}_{\mathit{a}}}{\sqrt{p}}+{\mathit{\omega }}_i$$

(1)

where ${\mathit{\omega }}_i~N(0,{\mathit{C}}_a/p)$ for some mean ${\mathit{\mu }}_a\in {\mathit{R}}^p$ and covariance ${\mathit{C}}_a\in {\mathit{R}}^{p\times p}$ of associated class ${\mathit{C}}_a$. Here, data are normalized by $\sqrt{p}$ to guarantee that $\left|\left|{\mathit{x}}_i\right|\right|=O\left(1\right)$ with high probability when $\left|\left|{\mathit{C}}_a\right|\right|=O\left(1\right)$.

Denote data matrix as $\mathit{X}=[{\mathit{x}}_1,\cdots ,{\mathit{x}}_T]\in {\mathit{R}}^{p\times T}$. Before extracting random spectral features, multiply X by random matrix $\mathit{W}\in {\mathit{R}}^{n\times p}$, which has independent identically distributed entries. Then, use some nonlinear activation function $\sigma (·)$ to obtain the random feature matrix $\mathit{\Sigma }\equiv \sigma \left(\mathit{W}\mathit{X}\right)\in R^{n\times T}$, whose columns are associated random features of ${\mathit{x}}_i$. And $\sigma \left(·\right)$ acts on each element of the vector. Then, we define the Gram matrix of random features as $\mathit{G}\equiv 1/n{\mathit{\Sigma }}^T\mathit{\Sigma }$, the entry(i,j) of which is given by:

$${\mathit{G}}_{ij}=\frac{1}{n}{\sigma \left(\mathit{W}{\mathit{x}}_i\right)}^{\mathit{T}}\sigma \left(\mathit{W}{\mathit{x}}_j\right)=\frac{1}{n}\sum _{k=1}^n{\sigma \left({\mathit{w}}_k^T{\mathit{x}}_i\right)}^{\mathit{T}}\sigma \left({\mathit{w}}_k^T{\mathit{x}}_j\right)$$

(2)

with ${\mathit{w}}_k^T$ the k-th row of W. Since ${\mathit{w}}_k$ follows the same known distribution, taking expectation over $\mathit{w}={\mathit{w}}_k$ of the above equation one results in the average kernel matrix $\mathsf{\Phi }$ with (i,j) entry of which given by:

$$\mathsf{\Phi }\left({\mathit{x}}_i,{\mathit{x}}_j\right)={\mathit{E}}_{\mathit{w}}{\mathit{G}}_{ij}={\mathit{E}}_{\mathit{w}}\sigma \left(w^T{\mathit{x}}_i\right)\sigma \left(w^T{\mathit{x}}_j\right)$$

(3)

Because we assume that the data obey a Gaussian mixture model, under which the data only contain two meaningful statistics: the mean and covariance. Therefore, we are interested in knowing whether different activation functions extract features with different degrees of bias for these two statistics. In order to give further valuable theoretical results, we start with some notational conventions.

Let ${\mathbf{C}}^0=\sum _{a=1}^K\frac{T_a}{T}{\mathit{C}}_a$, for $a\in \{1,\dots ,K\}$, $C_a^0={\mathit{C}}_a-{\mathbf{C}}^0.$ Then, for random elements,

$$\mathsf{\Omega }=\left[{\mathsf{\omega }}_1\cdots ,{\mathsf{\omega }}_{\mathrm{T}}\right]\in {\mathrm{R}}^{\mathrm{p}\times \mathrm{T}},\mathsf{\varphi }={\left\{{‖{\mathsf{\omega }}_{\mathrm{i}}‖}^2-\mathrm{E}{‖{\mathsf{\omega }}_{\mathrm{i}}‖}^2\right\}}_{\mathrm{i}=1}^{\mathrm{T}}\in {\mathrm{R}}^{\mathrm{T}}$$

and, for deterministic elements:

$$\mathbf{t}=1/\sqrt{\mathrm{p}}\left[\mathrm{t}\mathrm{r}\left({\mathbf{C}}_1^{\mathrm{o}}\right),\cdots ,\mathrm{t}\mathrm{r}\left({\mathbf{C}}_{\mathrm{K}}^{\mathrm{o}}\right)\right],\mathbf{M}=\left[{\mathsf{\mu }}_1,\cdots ,{\mathsf{\mu }}_{\mathrm{K}}\right]$$

$$\mathbf{J}=\left[{\mathbf{j}}_1\cdots ,{\mathbf{j}}_{\mathrm{T}}\right]\in {\mathrm{R}}^{\mathrm{T}\times \mathrm{K}}, \mathbf{S}=1/\mathrm{p}{\left(\mathrm{t}\mathrm{r}\right({\mathbf{C}}_{\mathrm{a}}{\mathbf{C}}_{\mathrm{b}}\left)\right)}_{\mathrm{a},\mathrm{b}=1}^{\mathrm{K}}$$

with ${\mathbf{j}}_{\mathrm{a}}$ denotes the canonical vector of ${Class}_a$ such that ${{(\mathbf{j}}_{\mathrm{a}})}_i={\delta }_{{\mathit{x}}_i\in {Class}_a}$.

In fact, under the assumption of a mixed Gaussian distribution, we can obtain an estimate of $\mathsf{\Phi }$, i.e., $\tilde{\mathsf{\Phi }}$, which is defined as:

$$\tilde{\mathsf{\Phi }}=d_1{\left(\mathit{\Omega }+\mathit{M}\frac{{\mathit{J}}^T}{\sqrt{p}}\right)}^T\left(\mathit{\Omega }+\mathit{M}\frac{{\mathit{J}}^T}{\sqrt{p}}\right)+d_2\mathit{U}\mathit{B}{\mathit{U}}^{\mathit{T}}+d_0{\mathit{I}}^T$$

(4)

$$\mathit{U}=\left[\frac{\mathit{J}}{\sqrt{p}},\mathit{\varphi }\right],\mathit{B}=\left[\begin{array}{cc}{\mathit{t}}^T\mathit{t}& \mathit{t}\\ {\mathit{t}}^T& \mathbf{1}\end{array}\right]$$

(5)

$\tilde{\mathsf{\Phi }}$ depends on the activation function only through the constants ${\mathrm{d}}_1$,${\mathrm{d}}_2$. Since data are viewed as a random variable, $\mathsf{\Phi }$ is therefore a random matrix. It can be shown that $‖\tilde{\mathsf{\Phi }}-\mathsf{\Phi }‖⟶0$, as $\mathrm{p},\mathrm{T}\to \mathrm{\infty }$. Note that $\mathsf{\Omega }$ and t as random matrix and vector, which portray the first- and second-order random fluctuations in the kernel matrix. While M, t, S are non-random matrices and vectors that characterize the statistical features that are naturally extracted after a complex nonlinear random mapping. Specifically for K categories of data, M characterizes the mean value of the data, while t and S characterize the information about the covariance of data. That is to say, activation functions extract features with varying degrees of bias for these two types of statistical information by only two constants, i.e.,${ d}_1$ and $d_2$. Thus, we can classify various activation functions according to their theoretical value.

Table 1 presents definitions and examples of different classes of activation functions. Mean-oriented functions are those whose $d_1\ne 0$, $d_2=0$, which asymptotically track only the difference in means. In Table 1, the mean orientation function is sign(t). In contrast, variance-oriented functions which asymptotically track only the difference in variance, since their $d_2\ne 0$, $d_1=0$, such as $\left|t\right|$, $\mathrm{c}\mathrm{o}\mathrm{s}\left(t\right)$ and $e^{-t^2/2}$. If $d_1\ne 0$ and $d_2\ne 0$, the function is balanced, it can extract both the difference in mean and variance, and the ratio can be adjusted by the parameters. We have max(t,0) as an example of balanced function in Table 1.

Empirical estimates of the differences in the means and covariances of datasets can be calculated from $‖{\mathit{M}}^T\mathit{M}‖$ and $‖{\mathit{t}\mathit{t}}^T+2\mathit{S}‖$, respectively. Therefore, by analyzing the differences in the mean and variance statistics of different data classes, we can select activation functions with varying biases to extract features for seizure diagnosis. This approach can enable us make informed decisions about choosing the most suitable activation functions for different types of datasets, improving the accuracy and reliability of the diagnostic process.

To better reflect the difference between the mean and variance of the dataset, we introduce coefficient of variation of dataset as:

$${\mathit{c}}_{\mathit{v}}=\frac{‖{\mathit{t}\mathit{t}}^{\mathit{T}}+\mathbf{2}\mathit{S}‖}{‖{\mathit{M}}^{\mathit{T}}\mathit{M}‖}$$

(6)

The higher the coefficient of variation of a dataset indicates that it is more desirable to extract random spectral features using variance-oriented activation functions. Conversely, random spectral features extracted using mean-oriented activation functions may have better classification results.

2.2.3. Traditional Features

In this paper, we choose several previous prevalent features that have been proven effective in EEG seizure recognition. They are correlation dimension (CD) [25], detrended fluctuation analysis (DFA) [26,27], Lyapunov exponent (LE) [28,29], Higuchi fractal dimension (HFD) [30,31,32], Hurst exponent (HE) [25], singular value decomposition entropy (SVD) [33], Fisher information (FI) [34], and Petrosian fractal dimension (PFD) [35] to describe the EEG state of patients. The feature extraction is performed by PyEEG in Python [52]. Since a single sample in the Bonn dataset has nearly four times as many sample points as in New Delhi dataset, calculating the correlation dimension and Lyapunov exponent can become very complicated. Therefore, we only use traditional features other than these two on the Bonn dataset. These traditional features, together with the wave numbers and random features proposed by us, are used as the input of SVMs (support vector machines) model to obtain the classification result of a seizure or not. Specifically, the detailed information of each traditional feature used in the algorithm is shown in Table 2:

2.3. Settlement and Evaluation

In our feature-based models, prevalent machine learning classifier SVM was used to judge whether patients are suffering from a seizure or not. In order to verify the two new types of features proposed in our algorithm, we performed ablation experiments on the features. To verify the advanced nature of the algorithm, experiments are performed to compare our method with the state-of-the-art methods.

To ensure that the effectiveness of the model depends only on the quality of the features, not on the intermediate steps, all feature-based methods will be trained using the SVM model with default parameters, i.e., the kernel function chooses a Gaussian kernel function with a penalty coefficient C = 1, gamma = 1/(number of features × sample variance). The 100 samples of length 1024 in the New Delhi dataset and the 200 samples of length 4027 in the Bonn dataset were randomly divided into the training and testing set of the original EEG data in the ratio of 6:4, and the corresponding features were extracted and used as SVM model’s input. For end-to-end methods, the raw EEG samples in the divided training set and testing set were directly used as inputs to them, which would automatically train and infer to output classification results.

To perform the reliability and the effectiveness of the obtained results, accuracy, precision, specificity, and sensitivity were calculated to estimate the proficiency of methods used in this work:

$$Accuracy=\frac{TP+TN}{TP+TN+FP+FN}\times 100\mathrm{\%}$$

(7)

$$Precision=\frac{TP}{TP+FP}\times 100\mathrm{\%}$$

(8)

$$Specificity=\frac{TN}{TN+FP}\times 100\mathrm{\%}$$

(9)

$$Sensitivity=\frac{TP}{TP+FN}\times 100\mathrm{\%}$$

(10)

To assess the complexity of the model and the likelihood of clinical application, we also calculated the cost of time required for model training and inference.

Due to the stochastic nature of the algorithm results, taking into account fairness and time cost, we conducted 1000 experiments for each of all feature-based methods and 20 experiments for neural network methods. Each method selected its best five results for calculating the average of the evaluation index scores to estimate the performance.

3. Results

3.1. Feature Extraction

3.1.1. Wave Numbers

As shown in Figure 3, the three types of epileptic waves were completely separated, and local maxima which were not peak points were not recognized as peaks. Additionally, epileptic waves numbers between seizure segments and non-seizure segments were significantly different. Only a small number of slow waves were detected in non-seizure samples during the interictal state, which is reasonable because normal electrical activity in the brain also produces slow waves, which tend to follow spikes and sharp waves in abnormal discharges. In contrast, the algorithm detected a large number of spikes, sharp waves, and slow waves in the seizure sample. On the one hand, this result, which is consistent with a priori knowledge, proves the validity of the algorithm. On the other hand, it also implies the rationality of using a number of waveforms as classification features.

3.1.2. Random Spectral Features

We first obtained the variance and mean difference estimates of each dataset, then the coefficient of variation can be calculated. The details are shown in Table 3 as follows.

Table 3 shows that the estimated mean difference and variance difference of the New Delhi dataset are 27.51 and 1989.05, respectively. And, the mean difference estimate for the Bonn dataset is 53.71 and the variance difference estimate is 6547.83. Obviously, both datasets are variance-oriented. So the variance-oriented activation functions, i.e., $\left|t\right|$, $\mathrm{c}\mathrm{o}\mathrm{s}\left(t\right)$ and $e^{-t^2/2}$ were chosen to extract random spectral features separately. And, the random spectral features were the first three eigenvectors of the Gram matrix given the activation function.

3.2. Ablation Experiments

We classified the traditional features used in the algorithm according to the nonlinear domain, time domain, and frequency domain, and then investigated the evaluation metrices after using the original single class features, adding wave number features, adding random spectral features and adding above two novel features at the same time. At the end, we also examined the enhancement effect of wave numbers, random spectral features, and two new types of features when using all traditional features for classification. On the two datasets, our experimental results are shown in Table 4 and Table 5.

From the tables, it can be seen that, in both datasets, the addition of wave numbers and random features significantly improved the classification accuracy of the model, regardless of using a single type of traditional feature or using all traditional features. For the effect of single-type features to improve accuracy, random features are superior to wave numbers. For instance, in the New Delhi dataset, after adding wave numbers to the traditional features, the accuracy reached 96.77%, while, when random features were added, the accuracy reached 99.88%. Among all kinds of traditional features, nonlinear features had the best single type classification accuracy, reaching 81.42%. After adding two new features, the accuracy reached 99.88%. The accuracy of the seizure detection reached 99.99% by combining the two new features with all the traditional EEG features.

As for the Bonn dataset, the random features and wave numbers still have a significant enhancement effect on seizure detection. And, nonlinear features still have the best single-type classification accuracy, reaching 90.15%, but after adding our two new features, the accuracy increases to 99.86%. The accuracy of seizure detection reached 100% by the fusion of our two new features and all the traditional EEG features. This shows that two kinds of new features are conducive to improving the detection effect.

To further validate the advantages of our proposed features over traditional ones, the features used as inputs to the classifier were projected to a two-dimensional space and visualized using the T-SNE (t-distributed stochastic neighbor embedding) algorithm. If the points of high-dimensional feature vectors from the different classes of samples after dimensionality reduction are dissimilar or distant from each other in the low-dimensional space, it means that the features have a strong discriminability.

In Figure 4, we show the visualized results of the T-SNE algorithm. It can be seen from the figure that, compared with traditional features, wave numbers and random features have a better differentiation degree. This is mainly because wave numbers are related to the gold standard of the clinical interpretation of EEG. And, random features are extracted according to the characteristics of the dataset, which have better data adaptability than the traditional features. This explains why our algorithm achieves optimal results on both benchmark datasets.

3.3. Contrast Experiments

Our algorithm was compared with the state-of-the-art methods mentioned in the related work section. Among the deep learning methods that serve as benchmarks for comparison are FCN (Wang et al., 2019 [14]), Resnet (Fawaz et al., 2019 [11]), TSLANet (Eldeleal et al., 2024 [12]), Hydra (Dempster et al., 2024 [18]), and Transformer + LSTM (Ignacio, 2023 [20]), which use raw EEG data as input to the model. Benchmarking methods based on features are entropy features (Mahshid and Zahra, 2023 [39]), time–frequency domain and nonlinear features (Chen et al., 2023 [38]), and neural network entropy features (Velichko et al., 2023 [40]). The comparative results of the model evaluation metrics in both datasets are presented in Table 6 and Table 7.

On both datasets, our method achieved the highest accuracy and specificity. Specifically, our method on the Bonn dataset achieved 100% in all metrics, which indicated that our algorithm did not have any leakage and misdetection on this dataset, and outperformed all the state-of-the-art algorithms. And, the accuracy was improved by 1.25% compared to the sub-optimal method.

In the New Delhi dataset, the Transformer + LSTM framework had better precision and sensitivity performance, but our algorithm also achieved sub-optimal 99.99%. In addition, the accuracy and specificity of Transformer + LSTM were 95.5% and 91.3%, while the accuracy and specificity of our algorithm were 99.99% and 100%. This showed that our judgment of the non-seizure segments was more accurate. In terms of seizure recognition, the gap between us and that deep learning algorithm is minimal. In summary, compared with deep learning methods, our algorithm outperformed or came close to them and was more interpretable.

Among all the most advanced feature-based methods, the method using time–frequency and nonlinear features had the best performance, achieving sub-optimal results on both datasets, with an accuracy of 99.34% on the New Delhi dataset and 96.2% on the Bonn dataset. Our algorithm was significantly superior to this method, which further reflected that the features we extract were more applicable. Additionally, compared with feature-based methods, our algorithm combined timing information and clinical experience, made it more conducive to the interpretation of features.

In addition, considering the complexity of the model affects the training cost and inference speed. Moreover, the inference speed further affects the physicians’ experience in using the intelligent algorithms. Therefore, we also compare the training time and inferring time of each model and present the results in Table 8.

As can be seen in Table 8, our algorithm on both datasets took significantly less time than deep learning models such as FCN, Resnet, and Hydra. It was also less time consuming than feature-based methods. Since the computation of traditional features was more complex, the inference of Transformer + LSTM and TSLANet on the New Delhi dataset took less time than ours, but our method for training was far less expensive than either. Moreover, because the complexity of computing wave numbers and random spectral features was very low, our method also took much less time than the other time–frequency domain, nonlinear features, and entropy features for feature extraction.

Our method took 0.49 s to infer on the Bonn dataset and 2.45 s to infer on the New Delhi dataset, which was also acceptable for clinical application. Therefore, taking into account the model accuracy and the time consumption, our model has the most potential to be applied in the clinic.

4. Discussion

From the experimental results of two datasets, both our proposed wave number features and random features significantly improved the performance of the feature-based seizure detection algorithm. This indicated that these two types of features were highly robust and effective. Additionally, it also outperformed all neural network approaches with significantly less time and computational cost.

Furthermore, the significance of our proposed coefficient of variation can be explained by comparing the coefficient of variation of the two datasets with their evaluation scores. Theoretically, the larger the coefficient of variation, the larger the difference in variance compared to the mean. The better seizure diagnosis was based on random spectral features in this case. It can be seen in Table 9 that the New Delhi dataset had a larger coefficient of variation than the Bonn one, and random features did work better on it. A Mann–Whitney test was also conducted for the accuracy of two datasets, and took “α = 0.05”. We obtained “p-value = 0.00 < α”, i.e., random spectral features were statistically more effective on the New Delhi dataset than on the other one.

More importantly, the wave number feature comes from clinical experience and the random spectral feature has a statistical mechanism. This makes our proposed features interpretable, and clearly our approach is more credible compared to those black box models.

5. Conclusions

In this paper, we proposed wave numbers and random spectral features and combined them with a machine learning classification algorithm to form an effective automatic seizure detection algorithm. The coefficient of variation we proposed to provide theoretical guidance for selecting suitable random spectral features for different datasets was also verified to be statistically valid. Our method achieved 99.99% and 100% accuracy on the New Delihi and Bonn epilepsy public datasets, respectively, demonstrating the generalizability of this method.

In ablation experiments, we verified that the proposed features can significantly improve the classification accuracy. In comparison experiments, our seizure detection algorithm achieved the best results on both benchmark datasets. Moreover, our method spent much less time than state-of-the-art deep learning methods in model training and met the requirements of clinical applications in terms of inference time consumption.

The datasets used in this study were all single-channel, and our approach will be extended to multi-channel clinical EEG data in the future. And the validity of the coefficients of variation of the datasets also needs to be validated in more EEG datasets. It is also an interesting topic to further study the theory and threshold of random feature selection based on the size of coefficient of variation.