Automatic Diagnosis of Epileptic Seizures in EEG Signals Using Fractal Dimension Features and Convolutional Autoencoder Method

Abstract

This paper proposes a new method for epileptic seizure detection in electroencephalography (EEG) signals using nonlinear features based on fractal dimension (FD) and a deep learning (DL) model. Firstly, Bonn and Freiburg datasets were used to perform experiments. The Bonn dataset consists of binary and multi-class classification problems, and the Freiburg dataset consists of two-class EEG classification problems. In the preprocessing step, all datasets were prepossessed using a Butterworth band pass filter with 0.5–60 Hz cut-off frequency. Then, the EEG signals of the datasets were segmented into different time windows. In this section, dual-tree complex wavelet transform (DT-CWT) was used to decompose the EEG signals into the different sub-bands. In the following section, in order to feature extraction, various FD techniques were used, including Higuchi (HFD), Katz (KFD), Petrosian (PFD), Hurst exponent (HE), detrended fluctuation analysis (DFA), Sevcik, box counting (BC), multiresolution box-counting (MBC), Margaos-Sun (MSFD), multifractal DFA (MF-DFA), and recurrence quantification analysis (RQA). In the next step, the minimum redundancy maximum relevance (mRMR) technique was used for feature selection. Finally, the k-nearest neighbors (KNN), support vector machine (SVM), and convolutional autoencoder (CNN-AE) were used for the classification step. In the classification step, the K-fold cross-validation with k = 10 was employed to demonstrate the effectiveness of the classifier methods. The experiment results show that the proposed CNN-AE method achieved an accuracy of 99.736% and 99.176% for the Bonn and Freiburg datasets, respectively.

1. Introduction

Epilepsy is a neurological disorder that has threatened many lives [1,2]. If diagnosed late, the epileptic seizure will cause irreparable damages to the patient [1,2]. Generally, epileptic seizures are divided into three groups of general [3], focal [4], and epilepsy with unknown symptoms [5]. According to the World Health Organization (WHO), more than 60 million individuals suffer from different types of epileptic seizures, and their lives are threatened with severe health issues [6,7].

Neuroimaging modalities are among the essential techniques for epileptic seizure detection, including functional and structural techniques [8,9]. The analysis of the anatomy and anatomical connections of various parts of the brain is performed using structural neuroimaging modalities, such as structural magnetic resonance imaging (sMRI) and diffusion tensor imaging (DTI) [10]. On the other hand, by using functional neuroimaging modalities, specialists can investigate the activities and functional connectivity of the brain accurately [8]. EEG is the primary functional neuroimaging modality with utmost importance for epileptic seizure detection [11].

EEG signals provide specialist physicians with essential information regarding brain performance at the time of epileptic seizures [11,12,13]. Compared to other neuroimaging methods, this modality is relatively cost-effective and has suitable performance. On the other hand, EEG has different channels at the time of recording the brain, making its signals more complicated [5]. EEG signals have different artifacts at the time of recording, such as urban electricity, muscles, etc. These issues create some challenges for specialist physicians in the diagnosis of epileptic seizures using EEG signals [14].

The diagnosis of epileptic seizures using EEG signals is of significant importance to specialist physicians. Nowadays, a wide range of research is being conducted on the implementation of computer-aided diagnosis system (CADS) based on artificial intelligence (AI) in neurological disorders [15,16,17,18], especially epileptic seizures detection [19,20,21,22]. The CADS based on AI consists of preprocessing, feature extraction, feature selection, and classification steps [19,20,21]. Feature extraction is the essential step of CADS in epileptic seizure detection using EEG signals [7]. The feature extraction techniques from EEG signals consists of time [22], frequency [23], time-frequency [24], and nonlinear [25] fields. EEG signals have chaotic features, which has led researchers to employ different feature extraction methods in studies on epileptic seizure detection, such as entropies [26], largest Lyapunov exponent (LLE) [27,28], FDs techniques, etc.

In this paper, nonlinear feature extraction techniques based on FDs are employed in the epileptic seizure detection. Figure 1 shows the block diagram of the proposed method for epileptic seizure detection. The Bonn [29] and Freiburg [30] datasets are employed to perform experiments. In the preprocessing step, all datasets are prepossessed using a Butterworth band pass filter with 0.5–60 Hz cut-off frequency. Following this, the EEG signals of both datasets are segmented into various time windows. Afterward, decomposition of the EEG signals into different sub-bands is carried out using DT-CWT [31,32].

The DT-CWT was first proposed by Kingsbury [33]. The DT-CWT has two parallel structures for obtaining real and imaginary parts of EEG signals. As compared to DWT, DT-CWT has several advantages, such as approximate shift immutability [34,35]. This helps to extract meaningful frequency information from each EEG signal [36]. The DT-CWT conversion is used in this paper with five decomposition levels, and a variety of fractal-based features are extracted from each decomposition level.

In the feature extraction step, various nonlinear features based on FDs are extracted from DT-CWT sub-bands, including HE [37], PFD [38], KFD [38], HFD [38], BC [39], MBC [39], Sevcik [40], DFA [41], RQA [42], MSFD [43], and MF-DFA [44]. FD features contain crucial information about EEG signals, and for this purpose they were used in the present study [45]. This research is the first to combine these feature extraction methods, highlighting the significance of this article.

In the feature selection step, the mRMR method [46] is employed. In this section, the proposed method selected 25 features. Ultimately, the KNN [47], SVM [48,49], and CNN-AE methods are employed for classification. The CNN-AE classifier method with a proposed number of layers is used for the first time in this paper.

Other sections of the article are presented as follows: Section 2 explores previous research carried out on epileptic seizure detection using the FD techniques. The proposed method for epileptic seizure detection is presented in Section 3. The experiment results are presented in Section 4. Section 5 presents statistical parameters for classification performance. Then, the discussion is presented in Section 6. Finally, conclusions and future works are presented in Section 7.

2. Related Works

This section presents a literature review of previous investigations into the diagnosis of epileptic seizures via EEG signals and FD feature extraction techniques. The first study by Wijayanto et al. [50] proposed an automated epileptic seizure detection based on the FD theory. First, the EEG signals were decomposed into five sub-bands. Afterward, the KFD was employed for feature extraction. Ultimately, the SVM method was employed for the classification of features. The results of the research showed that they could obtain favorable results.

El-Kishky et al. [51] employed different features for epileptic seizure detection. In this research, the EEG signals were first preprocessed. Afterward, various features, including HE, HFD, and Shannon entropy, were extracted from the EEG signals.

In another study, Jacob et al. [52] employed discrete wavelet transform (DWT) in the preprocessing step. Then, KFD and HFD features were extracted from DWT sub-bands. Finally, the SVM classifier was employed for classification.

Yuan et al. [51] proposed an intelligent method of epileptic seizure detection using EEG signals. In their proposed method, different techniques were employed in the preprocessing step. Afterward, differential box-counting (DBC) was investigated for feature extraction in EEG signals. They employed the SLFN method in the classification step and obtained favorable results.

In [52], the proposed method of epileptic seizure detection is implemented on the Freiburg dataset. In this study, the noises of the EEG signals were eliminated and then decomposed into time periods. Afterward, the FDs and gradient boosting methods were employed in the feature extraction, and classification steps, respectively.

Yang et al. [53] employed a new feature extraction method based on the FD theory for diagnosis of epileptic seizures using EEG signals. The dataset they employed for carrying out the experiments was of clinical type. The GFD feature extraction method was employed in this study. Following this, the run covering and SVM methods were employed in the feature selection and classification steps, respectively.

Moctezuma et al. [54] employed the CHB-MIT dataset to implement their proposed epileptic seizure diagnosis method. In this research, decomposition and preprocessing of EEG signals into different sub-bands was carried out using the empirical mode decomposition (EMD) method. Next, energy distribution and FD features were extracted from each EMD sub-band. Following this, K-means clustering and SVM methods were employed in the feature selection and classification steps, respectively.

In [55], an epileptic seizure diagnosis method via EEG signals based on the GFD technique was proposed. The Bonn dataset was employed to carry out experiments. In this research, the noises of the EEG signals were first eliminated in the preprocessing step and then decomposed into different sub-bands using DWT. The GFD method and the Renyi entropy were extracted from DWT sub-bands in the feature extraction step. Following this, the analysis of variance (ANOVA) test with the box plot method was employed in the feature selection method.

In [56], a method of epileptic seizure diagnosis using different FD methods was proposed. Different methods such as noise elimination, windowing of signals, and finally decomposition of EEG signals into different sub-bands were implemented in the preprocessing step. Next, the HFD, KFD, and Sevcik features were extracted from DWT sub-bands. Finally, the SVM was used for classification.

Sikdar et al. [57] employed the DWT to analyze the EEG signals of the Bonn dataset into different sub-bands. They investigated the MFDFA method for feature extraction. Finally, the statistical analysis and SVM methods were employed in the feature selection and classification steps, respectively.

Dalal et al. [58] used the FAWT method in the preprocessing step of the EEG signals. Then, they used the HFD and Kruskal–Wallis statistical test methods in the feature extraction and feature selection steps. Then, the SVM method was employed in the classification step.

This section introduced several articles on the diagnosis of epileptic seizures via EEG signals using the fractal technique. Other articles are introduced in

Table 1. Related works on diagnosis of epileptic seizures using FD features.

Works	Dataset	Preprocessing	Feature Extraction	Feature Selection	Classifier
[50]	Bonn	Decomposing into 5 Sub-bands	KFD	NA	SVM
[51]	Freiburg	Sliding window, Filtering	DBC	NA	SLFN
[52]	Freiburg	Filtering, Segmentation	FD	NA	Gradient Boosting
[53]	Clinical	Filtering	GFD	Run Covering	SVM
[54]	CHB-MIT	Segmentation, EMD	Energy Distribution Features, FDs	K Means	SVM
[55]	Bonn	Filtering, DWT	GFD, Rayni Entropy	ANOVA	NA
[56]	CHB-MIT	Filtering, Windowing, EMD	HFD, KFD, Sevcik	--	SVM
[57]	Bonn	Filtering, DWT	MF-DFA	Statistical Analysis	SVM
[58]	Bern-Barcelona	FAWT	HFD	Kruskal–Wallis Statistical Test	RELS-TSVM
[59]	Bonn	Filtering, Segmentation, Normalization	Shannon entropy, HFD, HE	Two-sample t-test	NA
[60]	Clinical	Filtering, DWT	HFD, KFD	NA	SVM
[61]	Clinical	Filtering	Average Entropy Feature and Hjorth Parameters	ANOVA	K Means
[62]	Clinical	Segmentation	GFD, Rényi Entropy	ANOVA	--
[63]	Bonn	Filtering, ATFFWT	Shift-invariance Property, Tunable Oscillatory Attribute, FD	Student’s t-test	LS-SVM
[64]	Bonn	Filtering, TQWT	FDs	Student T-Test, SFS	LS-SVM
[65]	Clinical	Filtering, DWT	FDs	GA	SVM
[66]	Bonn	Decomposition	HFD	--	SVM
[67]	Bonn	Filtering, Signal Decomposition Through Naive Method	HFD	--	--
[68]	Bern	Sliding Window Technique, Clustering	Undirected Graph Construction- Fractal Dimension (FD-NG)	--	SCA-SVM
[69]	Bonn	--	HFD, KFD	--	SVM

.

3. Material and Methods

In this section, the steps of the proposed method are presented. First, the EEG dataset for epileptic seizures detection is introduced. In the second step, the preprocessing of EEG signals is conducted according to the DT-CWT method. In the third step, various feature extraction techniques based on FD are introduced. The fourth step is pertinent to the feature selection technique. Finally, various types of classifier algorithms are presented in step five.

3.1. Dataset

3.1.1. The Bonn Dataset

The Bonn dataset consists of five sets of D, C, B, A, and E, and each one includes 100 single-channel signals with a time duration of 23.6 s [29]. All EEG signals from this dataset are recorded in the form of 12-bit with a sampling frequency of 173.61 Hz. In this dataset, the A and B sets of the EEG signals of five individuals are normal, where the A set of EEG signals is pertinent to the case of eyes open and the B set is pertinent to the case of eyes closed. The C and D sets include the EEG signals of five patients at the pre-ictal time [24]. Ultimately, the E set includes the seizure time of five patients with epilepsy with ictal time. In this dataset, the EEG signals of the A and B sets are recorded from the scalp, and the D, C, and E sets are recorded in an invasive way. Invasive electrodes are inserted into the hippocampal formation site symmetrically. All EEG signals are recorded by a 128-channel amplifier using a common, normal reference [29]. Further information regarding the dataset is presented in

Table 2. Thorough explanation of five subsets of dataset.

Sets	Subjects
	Patient Stage	Num. of Cases	Num. of Data	Length of Segments
Set A	Eye Open	5	100	4097
Set B	Eye Closed	5	100	4097
Set C	Seizure Free	5	100	4097
Set D	Seizure Free	5	100	4097
Set E	Seizure Activity	5	100	4097

and

Table 3. Further details about six problem classification.

Row	Classifications	Description
1	A-E	Healthy Control—Ictal
2	B-E	Healthy Control—Ictal
3	C-E	Interictal—Ictal
4	D-E	Interictal—Ictal
5	ABCD-E	Normal—Seizure
6	AB-CD-E	Healthy Control—Interictal—Seizure

. In addition, Figure 2 shows a sample of five sets of the Bonn dataset.

3.1.2. The Freiburg Dataset

In this dataset, the Intracranial EEG (IEEG) signals were recorded in the Freiburg hospital, Germany, from 21 patients with focal epilepsy [30]. All IEEG signals from this dataset are recorded with a sampling frequency of 256 Hz and strip, depth, and grid electrodes [25]. The age range of the patients in this dataset was between 10 and 50 years, and thirteen individuals were female, and eight individuals were male. In this dataset, there are different pre-ictal, ictal, and inter-ictal times of patients with epileptic seizures. Three different types of seizures are reported among the patients, including complex partial (CP), generalized tonic-clonic (GTC), and simple partial (SP) [30]. Each one of the patients has experienced at least two types of these seizures. Further details regarding this dataset are provided in

Table 4. Further details about the Freiburg dataset.

Patient	Age	Gender	Seizure Origin	Seizure Type	Number of Seizures
1	15	Female	Temporal	SP, CP	4
2	38	Male	Frontal	SP, CP, GTC	3
3	14	Male	Temporal	SP, CP	5
4	26	Female	Temporal	SP, CP, GTC	5
5	16	Female	Frontal	SP, CP, GTC	5
6	31	Female	Temporal	CP, GTC	3
7	42	Female	Temporal	SP, CP, GTC	3
8	32	Female	Temporal	SP, CP	2
9	44	Male	Frontal	CP, GTC	5
10	47	Male	Frontal	SP, CP, GTC	5
11	10	Female	Frontal	SP, CP, GTC	4
12	42	Female	Frontal	SP, CP, GTC	4
13	22	Female	Temporal	SP, CP, GTC	2
14	41	Female	Temporal	CP, GTC	4
15	31	Male	Frontal	SP, CP, GTC	4
16	50	Female	Temporal	SP, CP, GTC	5
17	28	Male	Temporal	SP, CP, GTC	5
18	25	Female	Temporal	SP, CP	5
19	28	Female	Frontal	SP, CP, GTC	4
20	33	Male	Temporal	SP, CP, GTC	5
21	13	Male	Temporal	SP, CP	5

[30]. Figure 3 shows a section of the EEG signals from the Freiburg dataset.

3.2. Preprocessing

This DT-CWT was first proposed by Kingsbury [33] and then improved by Ivan [34,35,36], in which two wavelet filter trees are used to obtain the wavelet coefficients corresponding to the real and imaginary parts of the wavelet [33]. Figure 4 shows the DT-CWT structural. As shown in Figure 4, $h_0$ and $h_1$ show low and high-pass filters for real coefficients. Furthermore, $g_0$ and $g_1$ represent low and high-pass filters for imaginary coefficients. In addition, (↓2) indicates a decrease in the sampling rate [33].

Here, the first and second filter banks are represented by symbols $\psi \left(t\right)$ and ${\psi }^{\prime }\left(t\right)$, respectively, and the complex wavelet base is defined as (1) [33]:

$${\Psi }_t\left(t\right)=\psi \left(t\right)+{\psi }^{\prime }\left(t\right).$$

(1)

The filter bank $\psi \left(t\right)$ is defined as follows [33]:

$$\psi \left(t\right)=\sqrt{2}{\displaystyle \sum }_n{h}_1\left(n\right)\varnothing \left(2t-n\right).$$

(2)

where

$$\varnothing \left(t\right)=\sqrt{2}{\displaystyle \sum }_n{h}_0\left(n\right)\varnothing \left(2t-n\right).$$

(3)

Equations are ${\psi }^{\prime }\left(t\right)$ also similar to Equations (1) and (2), where parameters $h_1$ and $h_0$ are substituted with $g_1$ and $g_0$, respectively. Ideally, ${\psi }^{\prime }\left(t\right)$ is equal to Hilbert transform $\psi \left(t\right)$ as follows [33]:

$${\psi }^{\prime }\left(t\right)=H\left\{\phi \left(t\right)\right\}.$$

(4)

An important step in establishing Equation (4) is that $g_0$ is equal to half the delayed sample of $h_0\left(n\right)$ [33].

$$g_0\left(n\right)\approx h_0\left(n-0·5\right).$$

(5)

In addition, in the case of orthogonal wavelet, Equation (4) is satisfied if [33]:

$$H_0^{\prime }\left(e^{jw}\right)=e^{-j0.5\omega }{H}_0\left(e^{j\omega }\right), for \left|\omega \right|<\pi .$$

(6)

Therefore, ideally, high-pass filters are defined as follows [33]:

$$H_1^{\prime }\left(e^{jw}\right)=-jsgn\left(\omega \right)e^{j0.5\omega }{H}_1\left(e^{j\omega }\right), for \left|\omega \right|<\pi .$$

(7)

For accurate DT-CWT running, the first step of the filter bank entails particular precision, otherwise frequency responses will not be accurately analyzed. Suppose the $k^{th}$ step of the first filter bank that terminates with $H_1\left(e^{jw}\right)$. If $H^{\left(k\right)}\left(e^{jw}\right)$ is the equivalent response, then we can apply new identities for k > 1, as represented as follows [33]:

$$H^k\left(e^{j\omega }\right)=H_1\left(e^{j{2}^{k-1}\omega }\right){\displaystyle \prod }_{m=0}^{k-2}{H}_0\left(e^{j{2}^m\omega }\right), for \left|\omega \right|<\pi .$$

(8)

If the equivalent response of the second filter bank of the corresponding branch, which is obtained by substituting $H_i^{\prime }$ with $H_i$, is shown with $H^{\prime \left(k\right)}\left(e^{jw}\right)$, then we have [33]:

$$H^1\left(e^{j\omega }\right)=H_1\left(e^{j\omega }\right).$$

(9)

$$H^{\prime \left(1\right)}\left(e^{j\omega }\right)=H_1^{\prime }\left(e^{j\omega }\right).$$

(10)

Assuming pairs $\left\{h_0,h_1\right\},\left\{g_0,g_1\right\}$ satisfy Equations (6) and (7), it follows that for all k ≥ 1, $H^{\prime \left(k\right)}\left(e^{jw}\right)$ and $H^{\prime \left(k\right)}\left(e^{jw}\right)$ estimate the following equation [33]:

$$H^{\prime \left(k\right)}\left(e^{j\omega }\right)=-j sgn\left(\omega \right)e^{j0.5\omega }{H}^k\left(e^{j\omega }\right), for \left|\omega \right|<\pi .$$

(11)

That is equal to [33]:

$$\frac{H^{\prime \left(k\right)}\left(e^{j\omega }\right)}{H^{\left(k\right)}\left(e^{j\omega }\right)}=\frac{H_0^{\prime }\left(e^{j\omega }\right)}{H_0\left(e^{jw}\right)}\frac{H^{\prime \left(k-1\right)}\left(e^{j2\omega }\right)}{H^{\left(k-1\right)}\left(e^{j2\omega }\right)}=-jsgn\left(\omega \right)e^{j0.5\omega }, for \left|\omega \right|<\pi .$$

(12)

Furthermore, if the $k^{th}$ response of the new step of the first and second filter banks are shown with $H_{new}^{\left(k\right)}\left(e^{jw}\right)$ and $H_{new}^{\prime \left(k\right)}\left(e^{jw}\right)$, respectively, then we can conclude:

We have pairs of $\left\{h_0,h_1\right\},\left\{g_0,g_1\right\}$ which satisfy Equations (6) and (7) for k > 1 [33]:

$$H_{new}^{\left(k\right)}\left(e^{jw}\right)=H\left\{H_{new}^{\prime \left(k\right)}\left(e^{jw}\right)\right\}.$$

(13)

If and only if [33]:

$$g_0\left(n\right)=h_0\left(n-1\right).$$

(14)

More information regarding DT-CWT is presented in [33]. In this study, we employed the Kingsbury toolbox, which covers various filters including Antonini (An), Q-shift (QS), LeGall (LeG), and near-symmetric (NS) filters [33]. Reference [70] shows that NS 13/19 and QS 14/14 filters give more successful results than other filters. Therefore, we performed our simulations using these filters. We considered five sub-bands because of the five important cerebral frequencies for epileptic seizures detection, including delta, theta, alpha, beta, and gamma waves. Therefore, each EEG signal in the datasets were segmented into five sub-bands, and then FDs features were extracted from each sub-band. Figure 5 shows the five frequency sub-bands of an EEG signal using DT-CWT with NS 13/19 and QS 14/14 filters.

3.3. Feature Extraction

Feature extraction is accounted for as the essential step of epileptic seizure diagnosis. Various nonlinear features are introduced in this section. These features are mostly based on the FD theories.

3.3.1. Hurst Exponent

The HE is employed to calculate the correlation in time series, including EEG signals [37]. HE is widely used to assess whether there is a long-term correlation in a time series. According to the interesting features of HE, it is employed for feature extraction from EEG signals [37]. The HE relationship can be expressed as follows [37]:

$$HE=\frac{\mathit{log}\left(\frac{R}{S}\right)}{\left(T\right)}.$$

(15)

where R is the difference between the maximum and minimum deviation and is pertinent to the mean. S indicates the standard deviation of the time series and T indicates the time duration of data [37].

3.3.2. Higuchi

The FD refers to the fractional or non-integer dimension of a geometric object. This dimension is a relative criterion of the number of structural blocks forming a pattern [38]. This feature is useful in different fields of the medical sciences [38]. Consider the $x\left(1\right),\left(2\right),\dots ,N$ time series. Consider $x_m^k$ Equation that expresses k number of time series:

$$x_m^k=\left\{x\left(m\right).x\left(m+k\right).x\left(m+2k\right)\dots x(m+\lfloor \frac{N-m}{k}\rfloor k)\right\}, for m=1.2\dots k.$$

(16)

In the above equation, m is the initial time, k is a delay, and $a$ is the integer part of a. Mean duration of $L_m\left(k\right)$ for each time series is defined as follows [38]:

$$L_m\left(k\right)=\frac{{\sum }_{i=1}^{\lfloor \frac{N-m}{k}\rfloor }\left|x\left(m+ik\right)-x\left(m+\left(i-1\right)k\right)\right|\left(n-1\right)}{\lfloor \frac{N-m}{k}\rfloor k}.$$

(17)

In the above equation, N and $\frac{\left(n-1\right)}{\lfloor \frac{N-m}{k}\rfloor k}$ are the duration of time series x and normalization factor, respectively. In Equation (18), the sum of mean durations L (k) for each k display is indicated [38].

$$L\left(k\right)=\sum _{m=1}^k{L}_m\left(k\right).$$

(18)

3.3.3. Katz

KFD is one of the FD techniques. In the KFD algorithm, the FD is calculated according to Equation (19) [38].

$$D=\frac{lo{g}_{10}\left(L\right)}{lo{g}_{10}\left(d\right)}.$$

(19)

where $L$ Is the length of the whole curve or the sum of distances between consecutive points in the time series, and $d$ is the estimated diameter that is equal to the distance between the first point of the sequence and a point in the sequence. Mathematically, $L$ and $d$ are calculated as follows [38]:

$$d=\underset{2\le i\le N}{\mathit{max}}\left\{\left|x\left(i+1\right)-x\left(i\right)\right|\right\}.$$

(20)

$$L={\displaystyle \sum }_{i=1}^{N-1}\left|x\left(i+1\right)-x\left(i\right)\right|.$$

(21)

where $i$ will be equal to the point with a maximum distance from $1$. The FD calculated by this method depends on the measurement unit. If the measurement unit changes, the calculated FD will change as well. The KFD method solved this problem using a general unit. $L$ and $d$ are normalized using the mean distance between consecutive points of $\underset{\_}{a}$. Therefore, Equation (19) can be expressed as follows [38]:

$$D=\frac{lo{g}_{10}\left(\raisebox{1ex}{$L$}\!\left/ \!\raisebox{-1ex}{$\underset{\_}{a}$}\right.\right)}{lo{g}_{10}\left(\raisebox{1ex}{$d$}\!\left/ \!\raisebox{-1ex}{$\underset{\_}{a}$}\right.\right)}.$$

(22)

The final Kartz equation can be expressed as follows by considering $n=\raisebox{1ex}{$L$}\!\left/ \!\raisebox{-1ex}{$\underset{\_}{a}$}\right.$ [38]:

$$D=\frac{lo{g}_{10}\left(n\right)}{lo{g}_{10}\left(\frac{d}{L}\right)+lo{g}_{10}\left(n\right)}.$$

(23)

3.3.4. Petrosian

PFD is a quick method of estimating FD [38]. The PFD method is expressed as Equation (24) [38]:

$$D=-\frac{lo{g}_{10}n}{lo{g}_{10}n+lo{g}_{10}\left(\frac{n}{n+0.4N\mathrm{\Delta }}\right)}.$$

(24)

where n is the length of the sequence and $N\mathrm{\Delta }$ is the number of changes in generated binary sequence [38].

3.3.5. Sevcik

The Sevcik algorithm is based on the Hausdorff dimension $D_h$, defined as follows [39,40]:

$$D_h=\underset{ϵ\to 0}{lim}\frac{-\mathit{log}N\left(\epsilon \right)}{\mathit{log}\left(\epsilon \right)}.$$

(25)

In the equation above, $N\left(\epsilon \right)$ is the number of balls with ε radius. By defining $N\left(\epsilon \right)=\frac{1}{2\epsilon }$, the equation is written as follows [39,40]:

$$D_h=\underset{ϵ\to 0}{lim}\frac{-\mathit{log}\left(L\right)-\mathit{log}\left(2\epsilon \right)}{\mathit{log}\left(\epsilon \right)}.$$

(26)

$$D_h=\underset{ϵ\to 0}{lim}1-\frac{\mathit{log}\left(L\right)-\mathit{log}\left(2\right)}{\mathit{log}\left(\epsilon \right)}.$$

(27)

In this method, the signals are normalized in two dimensions, which is 1 square meter. The normalized value of $x^{*}$ and ordinate $y^{*}$ are as follows [39,40]:

$$x_i^{*}=\frac{x_i}{x_{max}}. i=1,2,\dots ,N$$

(28)

$$y_i^{*}=\frac{y_i-y_{min}}{y_{max}-y_{min}}. i=1,2,\dots ,N$$

(29)

The square unit can be imagined as a network of N*N cells. The final equation is indicated in Equation (30) [39,40].

$$D_s=\underset{\acute{N}\to 0}{lim}1+\frac{\mathit{log}\left(L\right)-\mathit{log}\left(2\right)}{\mathit{log}\left(2\left(1-N\right)\right)}.$$

(30)

3.3.6. Box Counting

The BC method is one of the most important fractal methods, playing a crucial part in medical signals analysis [39]. In this method, the EEG signals are covered with a number of boxes, and their number is counted with a specific size to check how many are required for complete coverage. The mathematical definition is provided as follows [39]:

$$D_B=\underset{r\to 0}{lim}\left(\mathit{log}N\left(r\right)/\mathit{log}(\frac{1}{r})\right).$$

(31)

In this equation, $N\left(r\right)$ is the number of total boxes with r size. The counting box algorithm estimates the curve by the number of boxes required to cover the multiple box curve [39]. Thus,

$$\mathit{log}N\left(r\right)=D_B\mathit{log}(\frac{1}{r})+c.$$

(32)

where C is a constant, the lowest square of the best straight line is considered as the DB box-counting dimension estimation of the curve [39].

3.3.7. Multiresolution Box-Counting

Assume a discrete-time signal $s=\left\{s\left(1\right),s\left(2\right),\dots ,s\left(N\right)\right\}$ with the sampling frequency of $f_s$ with N samples. Each sample point s(i) is in the form of i = 1,…, N(X(i), y(i)). The steps of this method are defined in the following [39].

Step 1:

Consider two S(i) and S(i + 1) points indicating EEG signals with the interval of $dt=x\left(i+1\right)-x\left(i\right)=\frac{1}{f_s}$, whose points height is in the form of $h=y\left(i+1\right)-y\left(i\right)$. This method is repeated for all points on the EEG signal. In Equation (33), the total number of required boxes are provided to cover EEG signal with r resolution [39].

$$B\left(r\right)=sum\left(b\left(i\right)\right). i=1,\dots ,N-1$$

(33)

Step 2:

Each intermittent point is placed on the signal to obtain the temporal resolution of $r=\frac{2}{f_s}$. It is repeated in this temporal resolution, and the total number of required boxes is calculated to cover the entire curve [39].

Step 3:

We receive the number of boxes to cover the curve by repeating the steps above for multiple temporal resolutions. For $r=\frac{1}{f_s},\frac{2}{ f_s},\dots ,\frac{R}{f_s}$ that $\frac{R}{f_s}$ the maximum temporal resolution is true where the investigated curve is placed [39].

Step 4:

The linear regression coefficient of curve $\mathit{log}\left(B\left(r\right)\right)$ versus $\mathit{log}\left(\frac{1}{r}\right)$ is considered an EEG signal fractal dimension estimation [39].

3.3.8. DFA

In order to extract DFA, the method proposed by Peng et al. is employed [41]. Consider $\{x\left(i\right),i=1,2,\dots ,N\}$ time series. $\left\{y\left(k\right)\right\}$ Time series can be obtained according to Equation (34) using the difference sum of {x(i)} time series from its mean [41].

$$y\left(k\right)={\displaystyle \sum }_{i=0}^k\left[x\left(i\right)-\langle x\rangle \right].$$

(34)

In this equation, <x> is the mean of {x(i)} time series [41].

In the next step, the obtained $\left\{y\left(k\right)\right\}$ time series is divided into $N_n=int\left(N/n\right)$ identical windows with the length of n, and each window is a line on the point sets and is picked using the least square error method. Afterward, mean squared error function $F\left(n\right)$ for whole time series $\left\{y\left(k\right)\right\}$ and window length n is obtained according to Equation (35) [41].

$$F\left(n\right)=\sqrt{\frac{1}{N}{\displaystyle \sum }_{k=1}^N{\left[y\left(k\right)-y_n\left(k\right)\right]}^2}.$$

(35)

In this equation, $y_n\left(k\right)$ is the time series obtained from the line in each window [41].

3.3.9. Recurrent Quantitative Analysis

Using these recurrent maps provides the possibility of not only reviewing but also quantifying the hidden pattern and finding structural changes and the similarities among various patterns in data [42]. Known as RQA, this ability is a nonlinear method of processing time series, capable of quantifying the repeat number and range in a dynamic system, thus, creating its trajectory in the phase space [42]. This method’s various criteria are provided in the following [42].

(1) Recurrent Rate

Recurrent Rate (RR) specifies the percentage of recurrent points in a map pertinent to the total correlation defined as follows [42]:

$$RR=\frac{1}{N^2}{{\displaystyle \sum }}_{i,j=1}^N{R}_{i,j}.$$

(36)

(2) DET

It defines the percentage of existing recurrent points of diagonal lines to the total recurrent points. The diagonal lines consist of two or multiple points that are neighbors in the diagonal form without any other distance [42].

$$DET=\frac{{\sum }_{l=l_{min}}^{N}lP\left(l\right)}{{\sum }_{i,j=1}^N{R}_{i,j}}.$$

(37)

(3) Lmax Index

It defines the biggest diagonal line length [42].

$$L_{max}=max\left(\left\{l_i;i=1,\dots ,N_l\right\}\right).$$

(38)

(4) ENTR Index

It defines the entropy of diagonal lines [42].

$$ENTR={\displaystyle \sum }_{l=l_{min}}^{N}P\left(l\right)ln\left(P\left(l\right)\right).$$

(39)

(5) LAM Index

It defines the percentage of recurrent points creating the perpendicular (vertical) lines [42].

$$LAM=\frac{{\sum }_{v=v_{min}}^{N}vP\left(v\right)}{{\sum }_{v=1}^{N}vP\left(v\right)}.$$

(40)

(6) Trap Time (TT)

It defines the average length of vertical points [42].

$$TT=\frac{{\sum }_{v=v_{min}}^{N}vP\left(v\right)}{{\sum }_{v=v_{min}}^{N}P\left(v\right)}.$$

(41)

3.3.10. Margaos-Sun Fractal Dimension (MSFD)

The MSFD is one of the other fractal methods used for extracting the features from signals or time series [43]. Similar to the aforementioned fractal methods, this method is of paramount importance in multi-dimensional signals analysis. In this method, the EEG signals of both datasets are segmented with different time length windows, and afterward, MSFD is applied to them. Further explanations regarding this method are provided in [43].

3.3.11. MFDFA

In this section, the MFDFA algorithm for the time series $x_k$ signal of length N is presented. This technique consists of five steps, which are described below [44].

In the first step, the cumulative deviation relation $x_k$ for $y\left(i\right)$ is calculated using the Equation (42) [44].

$$Y\left(i\right)={\displaystyle \sum }_{k=1}^i\left[x_k-x\right]. , i=1,\dots ,N.$$

(42)

where

$$⟨x⟩=\frac{1}{N}{\displaystyle \sum }_{k=1}^N{x}_k.$$

(43)

In Equation (43), the signal $x_k$ converts to the time series $y\left(i\right)$ [44].

In the second step, y(i) is divided into Ns non-overlapping segments whose length is equal to s (Ns ≡ int (N/s)) [44].

In the third step, the least-squares for the polynomial series of order m are used to calculate distances v (v = 1,2,…, Ns) [44].

$$F^2\left(s,v\right)=\frac{1}{s}{\displaystyle \sum }_{i=1}^s{\left\{Y\left[\left(v-1\right)s+i\right]-y_v\left(i\right)\right\}}^2.$$

(44)

For each segment v, v = 1,…, Ns and

$$F^2\left(s,v\right)=\frac{1}{s}{\displaystyle \sum }_{i=1}^s{\left\{Y\left[N-\left(v-N_s\right)s+i\right]-y_v\left(i\right)\right\}}^2.$$

(45)

The effects of noise interference as well as non-static processes are eliminated at this stage. This leads to obtaining multi-fractal properties [44].

In the fourth step, the oscillation function of order qth is calculated for each scale s. This function has a qth RMS order [44].

$$F_q\left(s\right)={\left\{\frac{1}{2N_s}{\displaystyle \sum }_{v=1}^{2N_s}{\left[F^2\left(s,v\right)\right]}^{q/2}\right\}}^{1/q}.$$

(46)

Finally, the comparison powers of the log-log diagrams of Fq (s) are calculated for each q in the fifth step. The relation of this section is defined as follows [44]:

$$F_q\left(s\right)\approx S^h\left(q\right).$$

(47)

In this relation, h(q) is generalized Hurst power. The generalized HE is as follows [44]:

$$\tau \left(q\right)\approx qh\left(q\right)-1.$$

(48)

Note that a multi-fractal signal has multiple HEs. Moreover, the parameter τ(q) has a nonlinear relation with q [44].

In addition, it is proved that the singularity spectrum f (α) is related to h (q), which is shown in Equations (49) and (50).

$$\alpha =h\left(q\right)+q{h}^{\prime }\left(q\right).$$

(49)

$$f\left(\alpha \right)=q\left[\alpha -h\left(q\right)\right]+1.$$

(50)

The spectrum can be specified using a quadratic function with minimum squares. It is important to note that the quadratic function in question must lie in the vicinity of α₀.

$$f\left(\alpha \right)=A{\left(\alpha -{\alpha }_0\right)}^2+B\left(\alpha -{\alpha }_0\right)+C.$$

(51)

In this relation, C is an additive constant, and B is a measure of spectral asymmetry. Finally, the spectrum width is defined as follows [44]:

$$W={\alpha }_1-{\alpha }_2.$$

(52)

$$with f\left({\alpha }_1\right)=f\left({\alpha }_1\right)=0.$$

As indicated by width W, the spectrum in this relation is multi-fractal [44].

3.4. Information on Number of Features

This section describes the output values for each feature extraction technique. As discussed earlier, the DT-CWT technique was employed to decompose EEG signals into five sub-bands alpha, beta, gamma, theta, and delta. Then, a total number of 11 features were extracted from each DT-CWT sub-band, including MFDFA, MSFD, MBC, BC, Sevcik, DFA, HE, PFD, KFD, HFD, and RQA, each with its own output values. Additionally, some features had only one output and others possessed multiple outputs. The MSFD, Sevcik, DFA, HE, PFD, KFD, HFD, BC, and MRBC have a single output and MFDFA and RQA have multi outputs. The RQA method has six outputs and the MFDFA has different outputs. For the Bonn dataset with time windows of 5 s, the MFDFA has 51 outputs. Additionally, the MFDFA has 55 outputs for the Freiburg dataset with time windows of 4 s.

For the Bonn dataset, each EEG signal was segmented into time windows of 5 s. Time windows of 5 s have given the greatest accuracy in epileptic seizure detection according to [71]. For the Bonn dataset, the number of segments is equal to 2000, and 330 features (11 FDs) were extracted from each frame. As reported in [51], each EEG signal in the Freiburg dataset was segmented into time windows of 4 s. The number of features extracted for each frame of the Freiburg dataset was 350. Table 4 shows that the Freiburg dataset includes multiple subjects with various data lengths.

3.5. Feature Selection Using mRMR Method

mRMR is a feature selection method proposed by Peng et al. [46]. This method selects a group of features whose mutual information is minimal and with the output (class label) as the maximum value [46]. The mutual information of two features, x, and y, is defined as Equation (53), where p(x), p(y), and p(x, y) are the probability density functions (PDF) [46].

$$I\left(x;y\right)={{\displaystyle \iint }}^{ }p\left(x,y\right)\mathit{log}\frac{p\left(x,y\right)}{p\left(x\right)p\left(y\right)}{d}_x{d}_y.$$

(53)

The maximum correlation or dependence between the features and class label is calculated via the formula as follows [46]:

$$\mathit{max}D\left(S,c\right),D=\frac{1}{\left|S\right|}{\displaystyle \sum }_{x_iϵS}I\left(x_i;c\right).$$

(54)

where S is the set of features and c is the class label. In addition, the minimum redundancy among the features is calculated via the formula as follows [46]:

$$\mathit{min}R\left(S\right),R=\frac{1}{{\left|S\right|}^2}{\displaystyle \sum }_{x_i,X_{j\in S}}I\left(x_i,x_j\right).$$

(55)

In the end, the mRMR equation is calculated as follows [46]:

$$\mathit{max}\varnothing \left(D,R\right),\varnothing =D-R.$$

(56)

3.6. Classification Methods

The proposed classification algorithms are examined in this section. In this study, the proposed classification methods include KNN [47], SVM [48,49], and the CNN-AE method. The novelty of this section is to employ the CNN-AE method with the number of proposed layers. An explanation of each method is provided in the following section.

3.6.1. Support Vector Machine

The SVM method is one of the well-established, popular, and best classification methods [48,49]. This method is widely used in various 2-class and multiclass classification problems with significantly high performance. In the SVM method, the 2-class or multiclass classification is performed using various kernels such as linear and RBF [48,49]. In this article, the RBF kernel is used to classify binary- and multiclass problems.

3.6.2. KNN

The KNN method, also used as a classification method, operates such that it finds the nearest neighbor for each data (according to the Euclidean distance) [47]. Afterward, it selects the K-nearest neighbor, specifying the data label based on the label of the majority. This method operates significantly slowly. More information on this method is provided in [47]. In this article, K is considered to be three.

3.6.3. Convolutional Autoencoder (CNN-AE)

Nowadays, DL methods have a special role in various medical applications, and numerous investigations are being conducted in this field [72,73,74,75]. Various DL models have been provided so far, the combined CNN-AE architecture is among them. Unsupervised learning-based CNN-AE is used for input data training [76,77]. This architecture has two parts, decoder and encoder, each of which uses convolutional layers [76,77]. In this method, first, data are entered into encoder layers, received in the output data in compressed form. Then, the encoder section data is received in the decoder section, and the data recovery operation is performed. In this study, a CNN-AE model together with the number of proposed layers is employed. The proposed method is shown in detail in

Table 5. Details of proposed CNN-AE method.

Architecture	Layers	Filters	Kernel Size	Stride	Activation	Number of Features (Bonn Dataset)	Number of Features (Freiburg Dataset)
Encoder	Input Data	--	--	--	--	868	1024
	Conv1D	32	7	1	ReLU	868	1024
	Max Pooling	--	4	1	--	217	256
	Dropout	--	--	Rate = 0.25	--	217	256
	Conv1D	32	7	1	ReLU	217	256
	Max Pooling	--	4	1	--	54	64
Decoder	Conv1D	32	7	1	ReLU	54	64
	Up Sampling	--	4	1	--	216	256
	Conv1D	32	7	1	ReLU	216	256
	Up Sampling	--	4	1	--	864	1024
	Zero Pad	--	--	--	--	868	1024
	Conv1D	1	7	1	TanH	868	1024

.

A block diagram of the proposed technique is shown in Figure 6, where the parameters of each layer in the Bonn dataset are illustrated for more simplicity. EEG signals of each dataset are first applied to the network’s inputs. As shown in Figure 6, EEG signals pass through the encoder layers allowing them to extract the features of interest. The input EEG signals are reconstructed in decoder layers. The number of features extracted are presented in Figure 6 and Table 5. As illustrated in Figure 6, handcrafted features are merged with features of the CNN-AE model upon feature reduction via mRMR. Here, 25 handcrafted features are merged with features of the CNN-AE model within the concatenate block. Thereupon, data were classified into three fully connected (FC) layers, with the first two layers having 256 neurons and the third layer incorporating the Softmax activation function. The steps used to extract features and classify data for the Freiburg dataset are similar to those used on the Bonn dataset.

4. Statistical Parameters for Classification Performance

The assessment parameters have a crucial role in demonstrating the performance of classification algorithms. In this paper, accuracy (Acc), precision (Prec), specificity (Spec), sensitivity (Sens), and F1_score statistical metrics are used [71]. The pertinent equations of each have been illustrated below. In addition, the K-fold cross validation with k = 10 is used for evaluation [71].

$$Acc=\frac{TP+TN}{FP+FN+TP+TN}$$

(57)

$$Sens=\frac{TP}{FN+TP}.$$

(58)

$$Spec=\frac{TN}{FP+TN}.$$

(59)

$$FS=\frac{2 TP}{2 TP+FP+FN}.$$

(60)

5. Experiment Results

In this section, we present the results of the proposed method. Simulations were performed using a hardware system with 16 GB RAM, GPU Nvidia 1070, and Core i7 processor. In addition, Matlab 2019a software has been used to implement preprocessing, feature extraction, and feature selection. Classification algorithms including KNN, SVM, and CNN-AE have been implemented using Python language and Scikit-Learn [78] and Keras toolboxes [79].

The proposed method of this article includes four sections: preprocessing, feature extraction, feature selection, and classification. In the preprocessing step, firstly, we used a Butterworth band-pass filter with 0.5–60 Hz cut-off frequency for removal of EEG artifacts. Then, the EEG signals of the datasets were segmented into different time windows. In this step, the EEG window length for the Bonn dataset is 5 s, and for the Freiburg dataset is 4 s. The process of choosing time frames is of utmost importance for CADS efficiency, as it largely influences the precise processing of EEG signals for reaching the best detection accuracy. Similar to [71], signals were segmented into time windows of 5 s in the Bonn dataset. Likewise, the EEG signals in the Freiburg dataset were segmented into time windows of 4 s, similar to the method reported in [51]. Additionally, the DT-CWT is used for EEG signal decomposition into different sub-bands.

In the feature extraction step, various FD features were extracted from DT-CWT sub-bands. In this research, FD features consist of HE [37], PFD [38], KFD [38], HFD [38], BC [39], MBC [39], Sevcik [40], DFA [41], RQA [42], MSFD [43], and MFDFA [44]. The combining of FD features for feature extraction from EEG signals is the first innovation of the paper.

The execution time of feature extraction algorithms is of great importance. In practical applications, the execution time of each feature extraction algorithm plays a crucial role in terms of hardware implementation. We examine the execution time of each feature extraction algorithm in this part of the experiment. Accordingly, in

Table 6. The execution time of each FD method.

Dataset	HE	PFD	KFD	HFD	BC	MBC	Sevick	DFA	RQA	MSFD	MF-DFA
Bonn	0.0246	0.0027	0.0107	0.0203	0.0830	0.2183	0.0047	0.0474	0.1039	0.0072	0.1617
Freiburg	0.0242	0.00281	0.0114	0.02001	0.0827	0.2189	0.0051	0.0450	0.1043	0.0077	0.1622

, the execution time of each feature extraction algorithm is shown in the Bonn and Freiburg datasets. The values in Table 6 are in seconds. As evidenced, MSFD has the lowest execution time, whereas MFDFA has a greater execution time.

The third section is pertinent to feature selection. In this section, the mRMR method is employed for feature selection.

In the fourth section, various classification algorithms are used. These methods include KNN, SVM, and CNN-AE. The SVM classification algorithm was chosen with radial basis function (RBF) kernel and KNN method with k = 3. The CNN-AE method is a further novelty of the present study.

In this paper, the CNN-AE proposed model consists of two inputs. In the first input, the raw EEG signals of the datasets are applied to the proposed CNN-AE model. Then, CNN-AE features are combined with FD features. Results show that combining the handcrafted features with the CNN-AE features increases the CADS efficiency in epileptic seizure detection from EEG signals. Furthermore, in this part of the experiment, the K-Fold cross validation with k = 10 is used. In addition, to present logical results, the classifier algorithms have been run ten times under the same conditions, and their average results have been presented in tables. In this section, the results of the proposed method for the Bonn and Freiburg datasets are presented.

Table 7. Results for the Bonn dataset using different classifier methods without mRMR technique.

Methods	Sets	Accuracy (%)	Sensitivity (%)	Specificity (%)	F1_Score (%)
SVM	A-E	98.213 ± 1.23	97.933 ± 1.68	98.028 ± 1.37	98.832 ± 0.82
	B-E	98.131 ± 1.54	97.896 ± 2.11	97.842 ± 2.14	98.127 ± 1.08
	C-E	98.928 ± 1.07	97.871 ± 2.09	97.947 ± 1.11	98.935 ± 1.03
	D-E	98.769 ± 0.94	98.724 ± 1.23	98.022 ± 1.80	98.534 ± 1.48
	ABCD-E	98.692 ± 1.28	98.772 ± 0.87	97.979 ± 1.23	98.631 ± 1.11
	AB-CD-E	98.617 ± 0.76	98.634 ± 1.03	98.081 ± 1.14	98.603 ± 1.23
KNN	A-E	96.152 ± 2.36	97.307 ± 1.94	97.192 ± 2.31	96.286 ± 2.29
	B-E	97.024 ± 2.07	96.943 ± 2.48	96.414 ± 0.88	97.119 ± 1.35
	C-E	96.941 ± 1.71	97.219 ± 1.51	97.128 ± 1.24	96.268 ± 2.28
	D-E	96.803 ± 1.63	96.774 ± 2.62	96.915 ± 0.96	96.523 ± 1.16
	ABCD-E	96.538 ± 3.22	96.922 ± 3.07	96.896 ± 1.77	96.487 ± 1.28
	AB-CD-E	96.057 ± 2.04	96.231 ± 1.35	96.093 ± 2.23	96.342 ± 3.14
CNN-AE	A-E	99.613 ± 0.14	99.523 ± 0.21	99.652 ± 0.32	99.742 ± 0.27
	B-E	99.608 ± 0.37	99.431 ± 0.32	99.563 ± 0.35	99.638 ± 0.33
	C-E	99.482 ± 0.27	99.308 ± 0.28	99.798 ± 0.68	99.431 ± 0.25
	D-E	99.558 ± 0.38	99.378 ± 0.17	99.416 ± 0.39	99.483 ± 0.39
	ABCD-E	99.432 ± 0.24	98.604 ± 0.49	99.584 ± 0.18	99.309 ± 0.47
	AB-CD-E	99.389 ± 0.41	99.281 ± 0.28	99.375 ± 0.34	99.736 ± 0.22

provides the results of the proposed method without considering the mRMR feature selection algorithm for the Bonn dataset.

Table 8. Results for the Bonn dataset using different classifier methods alongside mRMR technique.

Methods	Sets	Accuracy (%)	Sensitivity (%)	Specificity (%)	F1_Score (%)
SVM	A-E	99.542 ± 0.41	99.477 ± 0.32	99.482 ± 0.41	99.389 ± 0.43
	B-E	99.384 ± 0.43	99.431 ± 0.53	99.528 ± 0.37	99.682 ± 0.21
	C-E	99.423 ± 0.26	99.569 ± 0.41	99.323 ± 0.59	98.761 ± 0.89
	D-E	99.348 ± 0.31	99.284 ± 0.28	99.411 ± 0.46	99.313 ± 0.51
	ABCD-E	99.212 ± 0.68	99.363 ± 0.39	99.268 ± 0.24	99.084 ± 0.85
	AB-CD-E	99.138 ± 0.52	99.212 ± 0.71	99.381 ± 0.37	99.026 ± 0.62
KNN	A-E	97.308 ± 2.04	97.493 ± 1.63	98.387 ± 1.43	97.482 ± 2.41
	B-E	97.177 ± 1.87	97.241 ± 2.08	97.715 ± 2.32	97.361 ± 1.63
	C-E	97.592 ± 1.44	97.137 ± 1.14	97.093 ± 0.96	97.452 ± 1.92
	D-E	97.431 ± 1.12	97.765 ± 1.96	98.934 ± 0.49	97.524 ± 1.09
	ABCD-E	97.474 ± 2.21	97.913 ± 1.83	98.887 ± 1.11	97.385 ± 1.31
	AB-CD-E	97.361 ± 1.85	97.137 ± 0.84	97.128 ± 2.71	97.276 ± 2.26
CNN-AE	A-E	99.834 ± 0.13	99.772 ± 0.19	99.823 ± 0.15	99.793 ± 0.33
	B-E	99.785 ± 0.29	99.467 ± 0.51	99.732 ± 0.18	99.809 ± 0.19
	C-E	99.718 ± 0.63	99.715 ± 0.22	99.581± 0.39	99.573 ± 0.28
	D-E	99.679 ± 0.26	99.413 ± 0.51	99.611 ± 0.21	99.644 ± 0.31
	ABCD-E	99.668 ± 0.41	99.649 ± 0.29	99.424 ± 0.33	99.818 ± 0.17
	AB-CD-E	99.594 ± 0.36	99.535 ± 0.41	99.662 ± 0.32	99.676 ± 0.24

presents the the results of using the mRMR feature selection algorithm for the Bonn dataset.

Table 7 shows that the proposed CNN-AE method has achieved satisfactory results. Compared to other classifier algorithms, the proposed method for binary and multi-class problems’ classification of the Bonn datasets has higher efficiency. In Table 8, we present the results of the proposed method with the mRMR method for the Bonn dataset.

As shown in Table 8, the use of the mRMR method alongside classifier algorithms led to higher diagnostic accuracy. In this section, we observe that the performance of KNN, SVM and CNN-AE algorithms has been significantly improved. Figure 7 shows the results of the classifier algorithms with the mRMR feature selection method for the different problem classifications in the Bonn dataset. Figure 7 shows that using the mRMR method has increased the efficiency of classifier algorithms in epileptic seizure detection.

Table 9. Results of proposed method for the Freiburg dataset.

Method	Accuracy (%)	Sensitivity (%)	Specificity (%)	F1_Score (%)
SVM	96.832 ± 1.57	97.036 ± 1.82	96.317 ± 2.12	97.109 ± 1.32
mRMR + SVM	98.718 ± 1.07	98.263 ± 0.86	97.828 ± 1.39	98.134 ± 0.74
KNN	95.111 ± 2.62	94.781 ± 2.84	95.017 ± 3.12	95.942 ± 2.53
mRMR + KNN	96.194 ± 1.94	95.962 ± 1.19	96.419 ± 2.08	96.081 ± 1.39
CNN-AE	99.072 ± 0.40	99.167 ± 0.71	98.431 ± 0.62	99.145 ± 0.44
mRMR + CNN-AE	99.176 ± 0.31	99.237 ± 0.24	99.697 ± 0.28	99.204 ± 0.38

shows the results of the proposed method for the Freiburg dataset. As shown in Table 9, the CNN-AE algorithm has been able to achieve the highest accuracy compared to other classifier algorithms when using mRMR.

Figure 8 shows the accuracy of various classifier algorithms with or without the mRMR method. As demonstrated in Table 7, Table 8 and Table 9, the proposed CNN-AE model has the highest detection accuracy for the Bonn and Freiburg datasets.

6. Discussion

This study provides a CADS for epileptic seizure detection via EEG signals using FD features and the proposed CNN-AE method. The proposed CADS include four sections: preprocessing, feature extraction, feature selection, and DL-based classification. First, the EEG signals of both datasets were decomposed into different time windows. Afterward, the DT-CWT technique was used for EEG signal decomposition into different sub-bands. In addition, the NS 13/19 and QS 14/14 filters were used for simulation of the DT-CWT.

Previous studies have indicated that EEG signals have a chaotic nature. Accordingly, using nonlinear feature extraction techniques, e.g., entropies [26], LLE [27], FDs [50,51,52,53,54,55,56,57,58,59,60,61,62,63,64,65,66,67,68,69], etc., leads to CADS performance enhancement in epileptic seizure detection [69]. Thereby, in the feature extraction step, different FD features are extracted from DT-CWT sub-bands. For the first time, different FD feature extraction methods are combined in this study, and this is considered to be the novelty of this paper.

In the following step, in order to reduce the size of the feature matrix, the mRMR method is used. In this step, the proposed feature selection method extracted 25 features.

In the final step, various machine learning (ML) and DL algorithms are used to classify the input data. The classifier methods consist of KNN, SVM, and CNN-AE methods. The CNN-AE method with proposed layers is another novelty of the paper. As previously mentioned, this architecture has two inputs, one of them is for applying raw EEG signals to the model, and the other input is for applying FD features. In the classification section, in order to calculate valid results, K-Fold with k = 10 is used. The CNN-AE architecture of this study is a feature fusion method since the FD-based handcrafted features and DL features are combined in this method to increase the accuracy of epileptic seizure detection. For the Bonn dataset, the proposed method achieved 99.736% accuracy for two-class classification problems and 99.594% accuracy for three-class classification problem. In addition, the proposed method on the Freiburg dataset achieved 99.176% accuracy. The experiment results show that the proposed method for both datasets achieved satisfactory results.

7. Conclusions and Future Works

Epilepsy is among a group of neural brain disorders taking place at different ages with various symptoms. If not diagnosed in the first stages, epileptic seizures will develop further [1,2,3,4,5]. Thereby, diagnosis of epileptic seizures in the first stages is of paramount importance for specialist physicians [6]. There have been various methods introduced to date for the early diagnosis of epileptic seizures. Among these methods, the recording of EEG signals is crucial for specialist physicians and neurologists. The EEG signals include important information regarding brain function at the time of epileptic seizures, which is significantly important for physicians. Recording EEG signals for epileptic seizure diagnosis is carried out in a long-term form, such that physicians cannot precisely examine all signals in detail. In addition, many factors, e.g., different artifacts in EEG signals, contribute to the fact that physicians encounter challenges in epileptic seizure detection.

In order to overcome the aforementioned obstacles, different AI methods are provided to help epileptic seizure diagnosis via EEG signals, and numerous investigations are being conducted. The AI techniques consist of various ML and DL methods. Each of these methods has some drawbacks and benefits.

This paper provides a novel CADS for epileptic seizure detection in EEG signals using FD features and a CNN-AE model with proposed layers. In this study, the steps of CADS for epileptic seizure detection from EEG signals include preprocessing, feature extraction, feature selection, and classification. Firstly, the Bonn and Freiburg datasets were used for experiments. During preprocessing, a band-pass filter (BPF) with 0.5–60 Hz cut-off frequency was employed for removal of artifacts in the EEG signals. In the next step, the EEG signals of the datasets were segmented into different time windows. For this, each EEG signal was segmented into time windows of 5 s and 4 s for the Bonn and Freiburg datasets, respectively. The DT-CWT with QS 14/14 and NS 13/19 filters was then used to decompose the EEG signals into different sub-bands. Feature extraction was performed employing various FD techniques. This step is the first innovation of the paper. Then, the mRMR technique was used for feature selection and the number of features was reduced to 25. The classification was eventually carried out using various ML models and a CNN-AE model, where the CNN-AE model with proposed layers is another novelty of the paper. The experiment results shows that the proposed CADS achieved the highest accuracy.

Table 10. Comparison of proposed method with other related works on the Bonn dataset.

Work	Dataset	Preprocessing	Feature Exraction	Dimension Reduction	Classification	Acc
[80]	Bonn	FSWT, CSoS	FuzzyEn, HFD	t-DSNE	KNN	93.62%
[81]	Bonn	TQWT	ApEn, HFD, KFD	--	SVM, RF	100%
[82]	Bonn	Windowing	Feature scaling	DNN		97.21%
[83]	Bonn	Standard Preprocessing	MFDFA, HE, FD, Combined Index	Statistical Analysis	N/A	N/A
[84]	Bonn	Windowing, DWT	Different Statistical Features	--	ANN	Multi-Class: 97.33%
[85]	Bonn	Windowing, DWT	Different Linear and Non-Linear Features	--	RF	95.60%
[86]	Bonn	Complete Ensemble Empirical Mode Decomposition (CEEMD) Method	Multiscale Dispersion Entropy (MDE)	Different Methods	ANN	Binary-Class:99%
						Multi-Class: 98.97%
[87]	Bonn	Enhanced Curvelet Transformation (ECT)	MGT. NPT. GLCM	PCA	RF	98%
[88]	Bonn	Filtering, Windowing, SVD	Histogram	--	Different Methods	N/A
[89]	Bonn	Windowing, DWT	Different Statistical Featuers	CCP, PCA	LSTM	99%
[90]	Bonn	STFT	Haralick Features	--	DT	92.50%
[91]	Bonn	EMD. MPF	1D-CNN		Softmax	98.60%
[92]	Bonn	TQWT	AVP and SVD Features	--	KNN	98.80%
[93]	Bonn	EMD with TQWT	Information Potential (IP) Features	Different Techniques	LS-SVM	99%
[94]	Bonn	Generalized Stockwell Transform (GSD)	SVM	--	RF	99.12%
Proposed	Bonn	Filtering, Windowing, DT-CWT	11 Fractal Dimension Features	mRMR	CNN-AE	Binary-Class:99.736%
						Multi-Class: 99.594%

and

Table 11. Comparison of proposed method with other related works of the Freiburg dataset.

Work	Dataset	Preprocessing	Feature Exraction	Dimension Reduction	Classification	Acc
[95]	Freiburg	Filtering	Different Entropies (ApEn, SampEn, PE, PFuzzy)	--	SVM	95.30%
[96]	Freiburg	Windowing, FFT	Proposed CNN Architecutre		Softmax	92%
[97]	Freiburg	DWT	DESA, Temporal and Spatial Averaging	Feature Aggregation	Different Methods	95%
[98]	Freiburg	Windowing, WPT	Different Features	--	Weighted ELM	N/A
[99]	Freiburg	EEG Signals Decomposition	CNN Model		N/A	N/A
[100]	Freiburg	Filtering	Different Linear and Non-Linear Features	--	SVM	96.8%
[101]	Freiburg	WT	Different Statistical Features	Bag-of-Words	SVM	N/A
[102]	Freiburg	Filtering	Proposed LSTM Model		N/A	97.75%
[103]	Freiburg	Filtering, FFT	DL Model (Integer-Net)		N/A	93.2%
[104]	Freiburg	Filtering	Different Linear and Non-Linear Features	--	SVM	97.5%
[105]	Freiburg	Filtering, Highly Adaptive Directional Time–Frequency Distribution (HADTFD)	TF-Flux, TF-Entropy, TF-Flatness	Spatial Averaging	Linear Classifier	98.56%
[106]	Freiburg	Windowing, DWT	1D-LBP Features	--	Different Methods	95.33%
[107]	Freiburg	Filtering	Different Linear and Non-Linear Features	Krill Herd Method	Proposed Classifier Method	98.9%
Proposed Method	Freiburg	Filtering, DT-CWT	11 Fractal Dimension Features	mRMR	CNN-AE	99.176%

display the related works on epileptic seizure detection in EEG signals using the Bonn and Freiburg datasets, respectively.

As shown in Table 10 and Table 11, the proposed method has been able to achieve satisfactory results compared to previous work. The proposed method includes a combination of the most important preprocessing, feature extraction, feature selection, and classification methods for epileptic seizure detection. This has improved the efficiency of the proposed method compared to related work.

The proposed method has higher accuracy than other methods used in practical applications in epileptic seizure detection in EEG signals, such that it can be implemented on hardware platforms such as field programmable gate array (FPGA) [108,109,110] and employed in health centers as an auxiliary tool for quick epileptic seizure diagnosis [111].

In future works, we could use various CNN-AE models such as convolutional De-noising AE (CNN-DAE) [112,113] models for epileptic seizure detection. First, the EEG signals were converted to 2D images using the short-time Fourier transform (STFT) technique and then applied to one of the inputs of the CNN-AE model [114]. Then, various 2D features of STFT were extracted and applied to the other input of the CNN-AE model [115]. Another future study is using local binary pattern (LBP) feature extraction methods [116,117], such that the EEG signals of LBP features are extracted and applied to the CNN-AE input. For future works, we can use new deep feature fusion methods for epileptic seizure diagnosis using EEG signals [118,119]. Another future study is using advanced DL models for epileptic seizure detection in EEG signals [120,121,122,123,124,125,126].