In this work, the low-dimensional dynamics of brain activity associated with acupuncture stimuli was probed. We found that manual acupuncture stimuli can reduce the dimensionality of brain activity, which results from the enhancement of oscillatory activity in the delta and alpha frequency bands induced by acupuncture. Moreover, it was found that large-scale brain activity could be approximated through the dynamics of a relatively simple attractor contained within a low-dimensional neural space, and the attractor's morphology was closely related to the frequency of acupuncture stimulation. These results shed light on the large-scale brain response to manual acupuncture stimuli.

**Author Contributions:** Conceptualization, X.G. and J.W.; methodology, X.G. and J.W.; software, X.G. and J.W.; validation, X.G. and J.W.; formal analysis, X.G. and J.W.; investigation, X.G. and J.W.; resources, X.G. and J.W.; data curation, X.G. and J.W.; writing—original draft preparation, X.G. and J.W.; writing—review and editing, X.G. and J.W.; visualization, X.G. and J.W.; supervision, X.G. and J.W.; project administration, X.G. and J.W.; funding acquisition, X.G. and J.W. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the National Natural Science Foundation of China, grant numbers 62071324, 61871287, and 61471265.

**Institutional Review Board Statement:** The study was conducted and approved by Institutional Review Board of Armed Police Logistics College Affiliated Hospital, P. R. China (protocol code LLKYPJ2010005, date of approval: 5 May 2010).

**Informed Consent Statement:** Informed consent was obtained from all subjects involved in the study.

**Data Availability Statement:** Not applicable.

**Acknowledgments:** This work was supported by the National Natural Science Foundation of China (Grant Numbers: 62071324, 61871287, and 61471265).

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

## *Article* **Evaluation of Feature Selection Methods for Classification of Epileptic Seizure EEG Signals**

**Sergio E. Sánchez-Hernández, Ricardo A. Salido-Ruiz, Sulema Torres-Ramos and Israel Román-Godínez \***

Division of Cyber-Human Interaction Technologies, University of Guadalajara (UdG), Guadalajara 44100, Jalisco, Mexico<EMAIL_ADDRESS>(S.E.S.-H<EMAIL_ADDRESS>(R.A.S.-R<EMAIL_ADDRESS>(S.T.-R.) **\*** Correspondence<EMAIL_ADDRESS>
**Abstract:** Epilepsy is a disease that decreases the quality of life of patients; it is also among the most common neurological diseases. Several studies have approached the classification and prediction of seizures by using electroencephalographic data and machine learning techniques. A large diversity of features has been extracted from electroencephalograms to perform classification tasks; therefore, it is important to use feature selection methods to select those that leverage pattern recognition. In this study, the performance of a set of feature selection methods was compared across different classification models; the classification task consisted of the detection of ictal activity from the CHB-MIT and Siena Scalp EEG databases. The comparison was implemented for different feature sets and the number of features. Furthermore, the similarity between selected feature subsets across classification models was evaluated. The best F1-score (0.90) was reported by the K-nearest neighbor along with the CHB-MIT dataset. Results showed that none of the feature selection methods clearly outperformed the rest of the methods, as the performance was notably affected by the classifier, dataset, and feature set. Two of the combinations (classifier/feature selection method) reporting the best results were K-nearest neighbor/support vector machine and random forest/embedded random forest.

**Keywords:** EEG; epilepsy; seizure detection; machine learning; features; feature selection

#### **1. Introduction**

Epilepsy is one of the most common neurological diseases, affecting around 50 million people of all ages globally [1]. The Center for Surveillance, Epidemiology and Laboratory Services of the United States of America estimated that in 2010, the number of adults with active epilepsy in the United States was 2.3 million; by 2015, the estimate increased to 3 million adults [2]. An additional study calculated that the cumulative incidence of epilepsy, in Norwegian children at the age of ten was around 0.66%, with 0.62% having active epilepsy [3]. The authors in [4] reported six studies of epilepsy prevalence in Mexico, which found prevalence rates of 3.9 to 42.2 per 1000 inhabitants.

An important tool for the diagnosis and management of epilepsy is the electroencephalogram (EEG). As mentioned in [5] (p. ii2), EEG is a "convenient and relatively inexpensive way to demonstrate the physiological manifestations of abnormal cortical excitability that underlie epilepsy." Other diagnostic techniques used in conjunction with EEG include neuroimaging, metabolic tests, and genetic tests. EEG can be processed and classified by using machine learning methods. Several studies have applied machine learning to classify ictal EEG signals [6,7] and predict seizures [8,9]. Concerning signal description, several metrics have been used to characterize such problems, and to build machine learning models, some of those metrics are computed from an EEG. Some researchers have approached the classification problem by calculating statistical, entropy, and univariate linear metrics from EEG [6,8–11]; those metrics can be fed to a model in the form of a vector or a matrix. In addition, metrics can be estimated for the entire EEG

**Citation:** Sánchez-Hernández, S.E.; Salido-Ruiz, R.A.; Torres-Ramos, S.; Román-Godínez, I. Evaluation of Feature Selection Methods for Classification of Epileptic Seizure EEG Signals. *Sensors* **2022**, *22*, 3066. https://doi.org/10.3390/s22083066

Academic Editors: Yvonne Tran and Christian Baumgartner

Received: 14 January 2022 Accepted: 13 April 2022 Published: 16 April 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

bandwidth or for smaller sub-bands [12–14], the latter with the intention of obtaining a more detailed view of the signal. Furthermore, transformation of the EEG signal to different domains has been explored by applying Fourier [15], short-term Fourier [7,16,17], wavelet transform [15,16,18–20], and contourlet transform [21]. It has also been analyzed as a graph [8] and an image [16]. As a result of the increasing interest in this topic, Ref. [22] presented a complete summary of several descriptors for time-domain, frequency-domain, and time–frequency-domain, along with their interpretation, while applying to the epileptic seizure detection on EEG signals.

Regarding the classification problem, a diversity of machine learning models have been tested for epilepsy prediction; among these we have: decision trees [23–25], support vector machine [23,26,27], K-nearest neighbor [23], and random forest [27]. In recent years, diverse deep learning models have been tried for epilepsy-related tasks, e.g., convolutional neural networks [16,17,28] and long-short term memory [8,28].

An important step in classifier modeling is feature selection. When performing feature selection, it is important to take into account several facts: (a) there are numerous attributes that can be calculated from EEG signals—each one describing a particular aspect of EEG; (b) there is a strong relation between features and model accuracy; (c) the curse of dimensionality, which is related to the difficulty of optimizing a solution in high-dimensional spaces; (d) the complexity and interpretability of the resulting models, which are about reducing the time and costs by training simpler learning models or selecting the features that are more relevant and meaningful from the problem perspective [29]. In this sense, several feature selection methods have been used to overcome the aforementioned issues. Some of those techniques are statistical tests, information gain [30], principal component analysis [31], permutation importance [32], and recursive feature elimination (RFE) [33], among others [29].

In this regard, there are several works in the state-of-the-art focused on epileptic EEG signal classification; some of them use feature selection methods (FSMs) to improve their results and reduce the dimension of the feature vectors.

For example, principal component analysis was applied in [34] to obtain less correlated features; however, their main goal was to evaluate the effects of channel selection on epileptic analysis over adults and children, without considering the effect of the feature selection method in the classification's performance, and only one classification method was considered, in this case linear discriminant analysis.

In [24], RFE was used to rank features and a support vector machine (SVM) for classifying epilepsy, autism, and control groups in children. Even though they determined which features and combinations of features contributed the most to the classification accuracy, they did not analyze either several feature selection methods or other classification methods.

The authors of [25] evaluated one feature selection method (recursive feature elimination) and one feature set, in combination with seven classifiers to improve the classification accuracy of automatic seizure diagnosis. From the 11 features that were calculated from the EEG signal, they reduced them to 8 features. All the experiments were performed only on one adult dataset.

In [12], six FSMs along with nine classifiers were used for automatic seizure detection. The FSMs were evaluated to rank and reduce the number of features, ranking the important features using a *t*-test and selecting the top 20 or 25, without testing additional alternatives that may result in different rankings. Furthermore, the authors provided experimental results based on one dataset that belonged to children. Additionally, the hold-out crossvalidation methodology was used, which is commonly used for bigger datasets.

FSMs have also been tested in signals other than EEG such as magnetic resonance images, as shown in [35]. The authors' objective was to compare three different FSMs (i.e., *t*-test filtering, the sparse-constrained dimensionality reduction model, and the support vector machine-recursive feature elimination) to determine which of them performed better when using an SVM as the classifier. However, the authors only tested the performance

of the FSM on one classifier (SVM). Furthermore, they only tested on one dataset without considering if their result may change when either tested on different datasets or FSM sets.

In summary, it can be observed that, even though there is extensive research about the feature selection methods applied in seizure and seizure-free EEG signals' classification, such approximations do not allow having a general perspective of the advantages and disadvantages of using a particular combination of the feature selection method with a classification algorithm (C-FSM) to determine: (1) the effect of the dimensionality reduction on the performance of the classifiers; (2) the best C-FSM combinations; and (3) the amount of coincidence of the best-selected features among FSMs; all of that considering two different feature sets and two different databases (adults and children).

Hence, in this work, the CHB-MIT and Siena Scalp EEG databases were used along with two different feature sets to evaluate the combination of six FSMs and five classification models. The results of this work allow determining: the minimum number of features that can be chosen for each FSM without scarifying the classifiers' performances; the performance of several C-FSM combinations in order to discover if a relationship exists between the FSMs with a particular classification algorithm; and if there is a feature or feature set that remains across different C-FSM combinations.

#### **2. Materials and Methods**

In this section, the features, models, and training procedures are described. A general methodology overview is shown in Figure 1.

**Figure 1.** Flow diagram of the methodology.

#### *2.1. Datasets*

CHB-MIT was one of the two datasets used for this research. The dataset is available for download at Physionet [36] under Open Data Commons Attribution License v1.0.

The data were collected at the Boston Children's Hospital. The database contains scalp electroencephalograms of 23 pediatric patients having epilepsy [37]. The number of recordings per patient varies from 9 to 45; all of them contain a metadata file listing the channels' names and ictal activity intervals; most of the records have a duration of one hour. EEG signals were sampled at 256 Hz; electrodes were placed according to the 10–20 system [37]. Most of the recordings are provided following a bipolar longitudinal montage.

After exploring the data, recordings not containing ictal activity were discarded. Thus, an amount of 137 EEG recordings containing an overall amount of 181 ictal activity intervals from 23 patients were accomplished.

The second dataset used in this research was the Siena Scalp EEG database; it is also available on the Physionet site under a Creative Commons Attribution 4.0 license. The dataset was collected by the Unit of Neurology and Neurophysiology at the University of Siena. It contains EEG recordings of 14 patients, 9 males (ages 36–71) and 5 females (ages 20–58). There are a total of 41 EEG recordings, and these include 47 ictal activity intervals; the recordings' duration is variable, from 1 to 13 h. The start and end of each seizure are detailed in the metadata file provided by the original authors [38]. EEG signals were sampled at 512 Hz. Electrodes were placed according to the 10–20 system [38]. EEG channels provided in the dataset are monopolar.

Recordings were further processed to design a dataset for a bi-class classification problem: seizure or seizure-free.

#### *2.2. Data Pre-Processing*

The number of conserved bipolar channels for the CHB-MIT recordings was 21. A few of the channels do not follow the 10–20 positioning, but channels were considered because these are included in every recording of the dataset. The Siena Scalp recordings were converted to a longitudinal bipolar montage, to have a fair comparison between both datasets; 18 channels were conserved.

A second-order Butterworth high-pass filter was used for removing frequencies below 0.5 Hz; then, a notch filter was applied to remove power line frequency (60 Hz and 50 Hz for CHB-MIT and Siena datasets, respectively).

A window length of 2 s was chosen. Epoch length was based on [13], as they mentioned, when extracting spectral features, it is important to choose small epochs due to the non-stationarity of the EEG.

As the number of ictal segments was significantly lower than the non-ictal ones, a 50% overlap was applied to ictal windows; besides that, some non-ictal windows were removed in order to keep an approximate ratio of 9:1 between both classes.

Ten percent of the CHB-MIT epochs were kept apart for adjusting the parameters of the classification models.

Inspired by [12], windowed EEG signals were separated into the following 5 subbands: complete bandwidth (0.5–30 Hz), delta (0.5–4 Hz), theta (4–8 Hz), alpha (8–12 Hz), and beta (12–25 Hz). Second-order Butterworth band-pass filters were used for band separation.

#### *2.3. Feature Extraction*

Two different features sets were evaluated to observe the effect of varying the metrics and determine which one outperforms (see Table 1). Therefore, half of Feature Set 1 (FS1) was conformed by statistical metrics and was applied in the time-domain or frequencydomain. Specifically, the median frequency (i) was estimated to characterize the power spectrum of the EEG data. On the other hand, the variance (ii), skewness (iii), and kurtosis (iv) were estimated in time-domain to characterize the variability and the distribution of the EEG data. Other features such as the peak frequency (v) were used to describe the frequency of the highest peak in the power spectral density; the root mean square (vi), range (vii), and the number of zero crossings (viii) in the time-domain were used respectively to estimate the effective value of the signal, to measure the maximum wave amplitude, and to count the number of points where the EEG wave cuts the horizontal axis, changing its state from positive to negative and vice versa.

Those features were selected based on previous studies where the authors performed an extensive review on the state-of-the-art on epileptic seizure detection based on EEG [22] along with some other works [6,8,10,12,25,34] that have successfully applied such feature metrics. In particular, Features i–vi were selected based on [12]. The number of zero crossings (ZCs) was estimated instead of the ZC rate because all EEG epochs had the same length.

On the other hand, FS2 was formed by the minimum (ix), which describes the minimum value that a signal can take, complexity (x), used to describe the change in frequency, mobility (xi), which is a measure of the mean frequency, interquartile range (xii), which is a measure of statistical dispersion, the spread of the data or observations, median absolute deviation (xiii), which is a robust measure of how spread out a set of data is, sample entropy (xiv), used for assessing the complexity of physiological time series signals, mean (xv), and standard deviation (xvi), statistical parameters that describe the average value and the amount of variability, or dispersion, from the individual data values to this average value. These features have been utilized across different studies for the classification of ictal EEG data [6,10–12,34]. Parameters used for the sample entropy (SampEn) are r = 0.2, m = 2. When applied to seizure discrimination in EEG, it was previously observed [44] that there is no best parameter combination, but several optimal combinations, one of which is the one used in this work.

The features from FS1 and FS2 were computed per channel and frequency band. For both feature sets (FS1 and FS2), feature vectors of 840 and 720 in length were obtained for the CHB-MIT and Siena datasets, respectively.

After that, correlations between features in all bands and channels were identified in both feature sets. For each pair of features, the Pearson correlation coefficient was computed. If a combination had a coefficient over 0.95, one of the features of the pair would be discarded. As a result, from FS1, 105 and 113 features were removed for the CHB-MIT and Siena datasets, respectively; on the other hand, from FS2, 304 and 257 features were removed for the CHB-MIT and Siena datasets, respectively.

#### *2.4. Classification Methods*

The following five algorithms were selected for classifying ictal EEG signals:

1. Decision tree (DT): This is a hierarchical model composed of decision nodes and terminal leaves, each leaf have an output label. Decision nodes implement a test function *f*(*x*), which is a discriminator dividing the input space into smaller regions. Among all possible splits, the DT looks for the one that minimizes impurity. There are several impurity measures, e.g., the Gini index and entropy. For a two-class problem, the Gini index is defined as [45]:

$$
\phi(p, 1-p) = 2p(1-p) \tag{1}
$$

where *p* is the probability of a sample reaching a node *m*, to belong to a class *C*. The classification and regression trees algorithm (CART) was applied in this research.

2. Support vector machine (SVM): This constructs a hyperplane or set of hyperplanes in a high-dimensional space that can be used for classification. Those hyperplanes have the largest distance to the nearest training data points (also known as the functional margin) [23]. The task of finding the optimal separating hyperplanes can be defined as [45]:

$$\min \frac{1}{2} ||w||^2 \text{ subject to } r^t(w^T x^t + w\_0) \ge 1, \forall t \tag{2}$$

where *w* are the parameters defining the hyperplane, *x<sup>t</sup>* are the instances of the training set, and *r<sup>t</sup>* is the actual label. If the problem is not linearly separable, the problem can be mapped to a new space by using non-linear basis functions [45].

3. K-nearest neighbor (KNN): This is a classifier that learns by analogy. A target unknown instance is compared to all the instances in the training set, locating the *k* closest instances; the algorithm assigns the class that corresponds to the majority. "Closeness" is measured by using a distance metric; in this study, we used the Manhattan distance (selected through a parameter grid search). The Manhattan distance for a *p*-dimension space is defined as [46]:

$$d(i,j) = |\mathbf{x}\_{i1} - \mathbf{x}\_{j1}| + |\mathbf{x}\_{i2} - \mathbf{x}\_{j2}| + \dots + |\mathbf{x}\_{ip} - \mathbf{x}\_{jp}|\tag{3}$$

**Table 2.** Structure of ANN layers.

To select the best-suited parameters for the above machine learning models, a grid search was performed; accuracy was utilized as a comparison metric. For selecting the ANN structure, a different number of neurons and hidden layers (1 and 2 layers) were evaluated, resulting in our best architecture, the one depicted in Table 2. The configuration settings and parameters of the grid search are listed in Table 3 for FS1 and in Table 4 for FS2. Both FSs used the same parameter grid detailed in Table 3. It is worth mentioning that the selected classification algorithms have been previously implemented for epilepsy-related tasks; however, the model parameters were not inspired by any specific work; on the contrary, they were selected by performing a grid search.

**Table 3.** Configuration setting per classification method. It was adjusted from FS1.

**Table 4.** Configuration setting per classification method. It was adjusted from FS2.

#### *2.5. Feature Selection Methods*

In this research, six FSMs were used. The metrics used to assign the importance of each feature are detailed below. The parameters selected for training each algorithm are mentioned in Table 5.

1. Decision tree (DT): The measure used to assign the feature importance is the Gini importance. As described in [47], the importance of feature *Xm* in an RF can be measured by Equation (4)

$$Imp(X\_{\mathfrak{m}}) = \frac{1}{N\_T} \sum\_{T} \sum\_{t \in T: v(s\_t) = X\_{\mathfrak{m}}} p(t) \triangle i(s\_t, t) \tag{4}$$

where *T* is a set of DTs, *v*(*st*) is the feature used to split node *t*, *i*(*st*, *t*) is the impurity decrease in node *t*, and *p*(*t*) is:

$$p(t) = N\_t / N\tag{5}$$

*N* is the number of training samples, and *Nt* is the number of samples reaching the node *t*. As this study uses a single DT, the size of *T* is 1.

$$\mathfrak{g}^{\mathfrak{r}}(\mathfrak{x}) = \operatorname\*{argmin}\_{\mathfrak{g}\in G} L(f, \mathfrak{g}, \pi\_{\mathfrak{x}}) + \Omega(\mathfrak{g}) \tag{6}$$

Function *L* measures the approximation of *g* to *f* in the locality defined by *πx*. Function Ω(*g*) is a measure of the complexity of *g*. As per [49], the exponential kernel was used for *πx*, the weighted square loss for *L*, and the linear model for *G*.

4. Shapley additive explanations (SHAP): This is a unified approach to interpreting model predictions. It assigns each feature an importance value or SHAP value. These are the Shapley values of a conditional expectation function of the original model [50]. The implementation applied in this study, kernel SHAP, optimizes (6), but it uses different forms of *π<sup>x</sup>* , *L*, and Ω:

$$
\Omega(\mathcal{g}) = 0 \tag{7}
$$

$$
\pi\_{\mathbf{x'}}(z') = \frac{M - 1}{(M \, close \, | \, |z'|) |z'| (M - |z'|)} \tag{8}
$$

$$L(f, \mathbf{g}, \pi\_x') = \sum\_{z' \in \mathbb{Z}} [f(h\_x(z')) - \mathbf{g}(z')]^2 \pi\_{x'}(z') \tag{9}$$

where *g* is the explanation model and follows a linear form, *f* is the classification model, *M* is the total number of features, and |*z* | is the number of used features. Parameter *<sup>z</sup>* is the set of features represented as {0, 1}*M*. As explained in [51], the function *hx* maps the 1s of *z* to the value from the instance to be explained *x*, and 0s are replaced by a random feature value of another instance sampled from the data.

$$r(f) = \frac{1}{\Sigma\_{\bar{j}}\frac{1}{r\_{\bar{j}}(f)}}\tag{10}$$

where *rj*(*f*) is the ranking position assigned to a feature by the rest of the FSMs.

Instances of the SHAP background dataset (Table 5) were estimated by applying K-means (*k* = 50) to a subset of the training dataset. This was done in order to use a small, but representative set of instances for the estimation of the SHAP values. For feature evaluation in ERF, first, RF was trained, then permutation importance was estimated in a separate dataset.

The DT and SVM were chosen due to simplicity and training speed. SHAP, ERF, and RR were selected because, in [54], these methods returned good performance and

consistency for a prediction task related to environmental datasets. LIME was considered because it is model-agnostic and allows model interpretation, similar to SHAP.

#### *2.6. Model Training and Evaluation*

Feature importance rankings were obtained for each combination of feature selection and classification method, and as a result, 30 feature rankings were estimated per dataset. For convenience, if any negative importance score was assigned to a feature, the absolute value was calculated. As LIME and SHAP compute feature importance per instance, the average was computed across all instances. When model training was required during ranking computation, 50% of the training dataset was passed to the model, and min–max normalization [46] was used for data scaling.

During the training and evaluation of the classification models, data were min–maxnormalized and a 5 × 2 cross-validation (2-fold, 5 repetitions) was implemented for model evaluation [55]. On each database (i.e., CHB-MIT and Siena), patients' epochs were merged into a single dataset; then, each new dataset was randomized after each iteration of the validation methodology.

The 5 × 2 CV F-test procedure was originally proposed to compare supervised classification algorithms, even though it has been previously implemented for comparison of FSMs (not applied to EEG data) [56,57].

First, the classification models were trained by keeping all the features in the training set and assessed in order to compute their classification performance. Then, features having the smallest ranking criterion were removed, and the model was re-trained and re-evaluated. Feature removal was performed in steps of 50 features at a time. There were 25, 12, 6, and 1 features also considered during the evaluation. Rankings were computed at the beginning of the process.

It should be noted that the above pipeline was repeated 6 times per classification method, as there were 6 FSMs. In addition, models were trained per each dataset separately.

#### *2.7. Computing and Software*

The experiments were run on 2 different computing devices: a computer with Intel Core i7, 12 GB of RAM, and Ubuntu 18.04 and a server with Intel Xeon Gold and NVIDIA Tesla P100. Python 3.7 was used for coding all the experiments. Numpy (1.19.5) [58], pandas (1.2.2) [59], and scipy (1.6.1) [60] were used for data engineering, scikit-learn (0.24.1) [61] and tensorflow (2.4.1) [62] for building machine learning models and feature selection algorithms, and matplotlib (3.3.4) [63] for plotting. Other needed libraries were lime (<IP_ADDRESS>) [49] and shap (0.39.0) [50]. Some processing pipelines were run on a Jupyter Notebook in order to visualize the charts.

#### **3. Results**

#### *3.1. Evaluation of Feature Dimensionality Reduction*

This analysis was performed in order to visualize the effect of the reduction of the feature vector size on the classification models' performance. By doing this, it was possible to have a general overview of the robustness of the classification models regarding the reducing of the feature vectors that may result from feature selection methods.

Figures 2–5 present five plots depicting the average F1-score for the combination classifier (C), feature selection method (FSM), and a number of features (NF) (C-FSM-NF). Every subplot corresponds to a different classification model, namely (a) decision tree (DT), (b) support vector machine (SVM), (c) artificial neural network (ANN), (d) random forest (RF), and (e) K-nearest neighbor (KNN). Colored lines indicate the average performance for different NF values, and each color indicates a different FSM.

When the CHB-MIT dataset and the FS1 were used, the best performances were reported by the ANN model and the worsts by the DT and SVM models. The former reported the best F1-score with 0.86 corresponding to the combination ANN-ERF-250/200 (see Figure 2c). Additionally, it was observed that for every classifier, the RR feature selection method (brown line) decreased from the early stages (see Figure 2).

The second-best-performing classifier was RF (see Figure 2d); most of the F1-score values were between 0.80 and 0.85. Most of the curves showed similar tendencies, but the RR curve (brown line) went down faster, again.

KNN (see Figure 2e) returned few F1-scores that overcame the RF classifier, but it was less stable (i.e., its performance was more dependent on the FSM). A pair of the curves was over 0.80 (yellow and purple lines corresponding to the SVM and ERF feature selection methods, respectively). The rest of them presented a diminishing tendency that started from the early beginning.

The classifiers having the lowest performance were DT and SVM (see Figure 2a,b). The F1-score curves of the DT classifier had values slightly under 0.75. On the other hand, some of the SVM's F1-scores reached 0.75.

**Figure 2.** F1-score of the model was obtained by using a different number of features, using the CHB-MIT dataset and FS1. Classification model used: (**a**) decision tree, (**b**) support vector machine, (**c**) artificial neural network, (**d**) random forest, and (**e**) K-nearest neighbor.

When the Siena dataset and the FS1 were used, several classification models returned lower performance in comparison with the CHB-MIT experiments (see Figure 3a–c,e).

The best classification model was RF (see Figure 3d), and the combination RF-ERF reached an F1-score of 0.85. SVM and RR decreased faster than the rest of the FSMs.

The ANN model reached a performance of 0.80 for some combinations (ANN-LIME, ANN-SHAP). The RR curve decreased faster than the rest of the FSMs. SVM and the DT had a better performance than RR, but not as good as LIME and SHAP (see Figure 3c). The DT (see Figure 3a) showed a steadier behavior compared to all the classifiers, but its performance was around 0.70

Figure 3b depicts an interesting pattern. The F1-score was 0.5 at the beginning of the training, then several curves decreased almost from the beginning of the training (brown, blue, and purple lines). On the other hand, the LIME curve rose markedly as the number of features reduced.

**Figure 3.** F1-score of the model is obtained by using a different number of features, using the Siena dataset and FS1. Classification model used: (**a**) decision tree, (**b**) support vector machine, (**c**) artificial neural network, (**d**) random forest, and (**e**) K-nearest neighbor.

When the CHB-MIT dataset and FS2 were used, the performance was better than the performance obtained during the FS1 experiments. Once again, the RR experiments tended to show poor performance and a faster decay in comparison to the rest of the FSMs.

The best performances were reported by the KNN model (see Figure 4e) and the worst by the DT and SVM models (see Figure 4a,b). KNN reported the best F1-score (0.90), and it corresponded to the experiments that removed a low number of features (e.g., 500 and 450). The combinations KNN-ERF-400/350/300/250/200 also reported an F1-score of 0.90.

The second-best-performing classifier was the ANN (see Figure 4c); most of the F1 score values were between 0.85 and 0.90. Most of the curves depicted a similar tendency, but the RR curve. RF (see Figure 4d) returned some F1-scores around 0.85, and these were slightly lower than the ANN scores.

When the Siena dataset and FS2 were used, the performances were a bit worse than the CHB-MIT experiments. The SVM classifier (see Figure 5b) depicted a similar pattern to the one observed in Figure 3b); when the number of features was reduced, the LIME and SHAP curves showed an increase in performance.

The best classification models were RF and KNN (see Figure 5d,e). These classifiers reported F1-scores around 0.85. For the KNN case, the RR curve (brown line) did not decay as fast as it did for the rest of the classifiers; however, RR was still the FSM with the worst performance.

Finally, to perform the comparison of the C-FSM combination and the feature selected, different cutoffs were selected. It is important to mention that based on the average of the F1-scores across every experiment, the decrease in classification performance between two consecutive cutoffs points was approximately equal, so based on the visual inspection of the Figures 2–5, four cutoffs were defined: 450, 150, 100, and 50 features. Notice that we discarded analyzing cutoff points under 50 features, because several models showed an F1-score lower than 0.6.

**Figure 4.** F1-score of the model is obtained by using a different number of features, using the CHB-MIT dataset and FS2. Classification model used: (**a**) decision tree, (**b**) support vector machine, (**c**) artificial neural network, (**d**) random forest, and (**e**) K-nearest neighbor.

**Figure 5.** F1-score of the model is obtained by using a different number of features, using the Siena dataset and FS2. Classification model used: (**a**) decision tree, (**b**) support vector machine, (**c**) artificial neural network, (**d**) random forest, and (**e**) K-nearest neighbor.

#### *3.2. Comparison of C-FSM Combinations*

This analysis intended to observe, in detail, the performance of several C-FSM combinations. To do so, several cutoff points were chosen, then, to every cutoff point, the F1-score, sensitivity, and accuracy of all the combinations for a model classifier and feature selection method were computed. For the best-performing combinations of C-FSMs, the 5 × 2 CV F-test was applied to find statistically significant differences between the error rates.

Tables 6–9 depict the different C-FSM tuples and their respective F1-score, sensitivity, and accuracy for different sizes of the feature vector. In particular, these tables show the performances when the 450, 150, 100, and 50 best features were kept.

When 450 features were used, it is depicted in Table 6 that there was not an FSM that outperformed the rest. The best performances were obtained by the KNN and ANN models, the former using the FS2 and the combination KNN-SVM/SHAP/ERF (0.90) and the latter using FS1; the best combination was ANN-SVM/SHAP/ERF (0.84). On the other hand, the F1-scores of the Siena dataset experiments were lower, on average. The combinations having the best performance were KNN-SVM/LIME (0.83) and KNN-DT/SVM/LIME/SHAP/ERF/RR (0.86), for FS1 and FS2, respectively.

Observe that the experiments using FS2 tended to report better performances than those in the FS1 experiments, no matter the dataset used. Furthermore, it should be noted that there were large differences between the accuracy and the F1-score values, due to the large class imbalance. The performances of the DT and SVM classifiers were noticeably lower than ANN/RF/KNN, the SVM being the classifier with the worst performance values.

**Table 6.** Performance comparison using 450 features. The largest F1-scores are displayed in **red bold** and **blue bold** for FS1 and FS2, respectively. F1s = F1-score, Se = sensitivity, Acc = accuracy.

Table 7 shows the classification metrics for the 150-feature experiments. The best performances, using the CHB-MIT dataset and FS1, were returned by the combination ANN-LIME/SHAP/ERF (0.85); for the FS2 case, this was KNN-SVM/ERF (0.89). On the other hand, the Siena experiments showed that the best-performing combinations included the ERF as an FSM; RF-ERF returned the largest F1-scores, and these were 0.85 and 0.86 for FS1 and FS2, respectively.

**Table 7.** Performance comparison using 150 features. The largest F1-scores are displayed in **red bold** and **blue bold** for FS1 and FS2, respectively. F1s = F1-score, Se = sensitivity, Acc = accuracy.

In the 100-feature case (see Table 8), the behavior was similar to the 450-feature case, in the sense that the best-suited FSM depended on the feature set and dataset. For the CHB-MIT dataset, the combinations with the best performances were KNN-SVM (0.85) and KNN-DT/SVM/ERF (0.87) for FS1 and FS2, respectively. For the Siena dataset, the best performance was reported by RF-ERF (0.85 and 0.86).

**Table 8.** Performance comparison using 100 features. The largest F1-scores are displayed in **red bold** and **blue bold** for FS1 and FS2, respectively. F1s = F1-score, Se = sensitivity, Acc = accuracy.

**Table 8.** *Cont*.

Table 9 shows that classifiers accounting for the best performances, RF and KNN, 0.84 and 0.85 being the largest values for both of the datasets. It is interesting to note that the best performances were similar for both datasets and the feature sets.

**Table 9.** Performance comparison using 50 features. The largest F1-scores are displayed in **red bold** and **blue bold** for FS1 and FS2, respectively. F1s = F1-score, Se = sensitivity, Acc = accuracy.

It is worth noticing that the RF classifier presented the most steady performance (≈ 0.8) no matter the FSM, dataset, and feature reduction (see Tables 6–9).

In order to identify the significant differences between FSMs, Tables 10 and 11 show the results of the F-test; given a dataset, a number of features, and a feature set (FS1 or FS2), the F-test was applied to the results of the 5 × 2 CV experiments. During the testing, the best-performing C-FSM (as per the F1-score) was compared to the rest of the FSMs. It must

be considered that the 5 × 2 CV F-test evaluates the error rates, not the F1-scores, so the accuracy is computed and depicted in Tables 6–9.

In the case of several combinations having the same F1-score, the sensitivity and accuracy were considered to choose the best-performing combinations.

**Table 10.** Comparison of the best C-FSM when FS1 is used. It is compared against the rest of the FSMs. • denotes *p* < 0.05.

**Table 11.** Comparison of the best C-FSM when FS2 is used. It is compared against the rest of the FSMs. • denotes *p* < 0.05.

The 450 feature section of Tables 10 and 11 shows that most of the test results were not statistically significant (*p* > 0.05). In Table 10, there were three cases where the FSM error rates resulted in being statistically different from the rest of the FSMs, and these were RF-ERF-150 (Siena dataset), RF-ERF-100 (Siena dataset), and KNN-SVM-50 (CHB-MIT dataset). It must be noted that the differences in accuracy may be small (see Table 6–9), even if there are statistical differences.

In Table 11, there are no cases where an FSM was statistically different from the rest of the FSMs. It was observed that the RR experiments tended to show a statistical significance that did not depend on the number of features or the dataset.

Finally, Tables 10 and 11 show that the most common combinations were KNN-SVM and RF-ERF.

#### *3.3. Comparison of Selected Features*

To determine if there were coincidences in the features selected by the FSMs and if the FSMs assigned more importance to the same features, the Jaccard index [64] was used to calculate the similarity by pairs of FSMs. For this analysis, the FS2 experiments were chosen because those experiments produced higher performances in comparison with FS1.

When 450 features were used for training the KNN classifier (see Figure 6a,e), the feature sets practically overlapped; indices had values over 0.85 for the CHB-MIT dataset and the Siena dataset. This behavior was expected, as it is highly probable to select similar features when the number of features in a dataset is large, so the 450 case will not be further discussed.

**Figure 6.** The similarity of the sets of features selected per each feature selection method when KNN was used as a classifier. CHB-MIT dataset: (**a**) best 450 features, (**b**) best 150, (**c**) best 100, and (**d**) best 50. Siena dataset: (**e**) best 450 features, (**f**) best 150, (**g**) best 100, and (**h**) best 50.

When 150 features were used during training (see Figure 6b,f), the Jaccard index for SHAP-LIME (0.43) was the largest of all the combinations when the Siena dataset was used; it was equivalent to 91 out of 150 features. For the CHB-MIT dataset, RR-ERF (0.19) and RR-DT (0.19) obtained the largest values. Figure 6c,d,g,h show a low coincidence for most of the FSMs when 100 or 50 features were used, and this applies to both datasets. The only notorious index was 0.41 belonging to SHAP-LIME (see Figure 6g).

In the case of the SVM classifier, by using 150 features (see Figure 7b,f), the largest values were obtained by SHAP-LIME (0.43 and 0.72); the FSMs coincided in selecting 91 features for the CHB-MIT dataset and 126 for the Siena dataset. The Jaccard indices computed for 100 features (see Figure 7c,g) showed a good similarity for SHAP-LIME (0.41), and the selected feature sets coincided in 59 features for both of the datasets. When 50 features were used for training (see Figure 7d,h), several combinations returned a similarity value of 0.2.

**Figure 7.** The similarity of the sets of features selected per each feature selection method when SVM was used as a classifier. CHB-MIT dataset: (**a**) best 450 features, (**b**) best 150, (**c**) best 100, and (**d**) best 50. Siena dataset: (**e**) best 450 features, (**f**) best 150, (**g**) best 100, and (**h**) best 50.

The ANN classifier followed a similar pattern as SVM; SHAP-LIME obtained larger indices than most of the combinations, and these were 0.72 and 0.52 for 150 features, 0.71 and 0.50 for 100, and 1.0 and 0.33 for 50 (see Figure 8b–d,f–h), respectively. It should be noted that a Jaccard index of 1 means that both sets totally overlapped.

**Figure 8.** The similarity of the sets of features selected per each feature selection method when ANN was used as a classifier. CHB-MIT dataset: (**a**) best 450 features, (**b**) best 150, (**c**) best 100, and (**d**) best 50. Siena dataset: (**e**) best 450 features, (**f**) best 150, (**g**) best 100, and (**h**) best 50.

Interestingly, when the DT classifier was used, the SHAP-LIME similarity was equal to or lower than other combinations. When using 150 features, the largest indices were returned by the combination RR-SHAP (see Figure 9b,f). Figure 9c,g,h show that the largest similarity values were returned by combinations including LIME or SHAP.

**Figure 9.** The similarity of the sets of features selected per each feature selection method when DT was used as a classifier. CHB-MIT dataset: (**a**) best 450 features, (**b**) best 150, (**c**) best 100, and (**d**) best 50. Siena dataset: (**e**) best 450 features, (**f**) best 150, (**g**) best 100, and (**h**) best 50.

Comparable to the DT, the RF case showed the largest values when a combination included LIME or SHAP. A remarkable fact is that there were two combinations having an index value of 0.5 (see Figure 10c,d).

**Figure 10.** The similarity of the sets of features selected per each feature selection method when RF was used as a classifier. CHB-MIT dataset: (**a**) best 450 features, (**b**) best 150, (**c**) best 100, and (**d**) best 50. Siena dataset: (**e**) best 450 features, (**f**) best 150, (**g**) best 100, and (**h**) best 50.

In order to compare the selected features for the best C-FSM combinations in Table 9, Figures 11 and 12 show the top-10 features for KNN-SVM and RF-ERF, respectively. The former figure corresponds to the CHB-MIT dataset, while the latter to the Siena dataset. Each feature is defined as follows: EEG band/bipolar channel/metric. For example, in Figure 11, the feature with the greatest importance value is "alpha\_FP1-F3\_skew", which indicates that the most important feature was the skewness measured in the bipolar channel FP1-F3 on the alpha band. It is important to note that the selected features may vary due to several factors, including seizure type, epileptogenic region, and EEG montage, among others.

**Figure 11.** Features having the largest assigned importance. The results correspond to the combination KNN-SVM and the CHB-MIT dataset. Importance values are not normalized.

**Figure 12.** Features having the largest assigned importance. The results correspond to the combination RF-ERF and the Siena dataset. Importance values are not normalized.

#### **4. Discussion**

In the evaluation of feature dimensionality reduction analysis, it was observed that the best machine learning classifiers were ANN, RF, and KNN, taking into account neither the database nor the FS. This is evident by looking at the performance tendency of every combination of C-FSMs (colored lines). It is important to note that there was a recurrent behavior in all the combinations, that is almost all of them started showing a performance decrease when the feature vectors were around a length of 50. Lower than that, the scores began to be around 0.6 or less. Therefore, a feature vector of a length of 50 is the minimum suggested for having good classification performance using some of the C-FSM combinations, for FS in epilepsy databases. Moreover, if FS2 is used, the scores achieved by most of the C-FSM combinations were improved and steadier. In this sense, the DT and RF classifiers were less affected by the dimensionality reduction no matter the database, the FS, nor the FSM. This could be explained by the nature of those algorithms, that is both algorithms use a ranking metric to determine the importance of each feature vector in a particular classification task.

Regarding to the comparison of C-FSM combinations, it can be noticed that the models that were trained with CHB-MIT data had better F1-scores than the ones trained with Siena data. We considered that the performance difference was due to the CHB-MIT database being larger than the Siena Scalp EEG database. Actually, the number of CHB-MIT epochs used during experimentation was more than twice the number of Siena epochs. Another aspect to consider is the number of channels, 21 and 18 for the CHB-MIT and Siena datasets, respectively.

In addition, we can observe in Tables 6–9 that the greater the number of features, the ANN and KNN showed better F1-scores; on the other hand, the lower the number of features, KNN kept showing the best F1-scores, and RF emerged with better F1-scores than the ANN. Note that, even though KNN presented the highest F1-scores, the difference between KNN and RF reduced as the number of features decreased.

Furthermore, in Tables 10 and 11, it can be seen that the smaller the number of features, the greater the number of significant differences was between the best C-FSM and the other FSMs. In this sense, these significant differences in the classifiers' performance indicate a relationship between the classifier and the feature selection method; in this case, for a lower number of features, KNN was better when using SVM and ERF, while RF was better when using ERF.

The main purpose of this study was not to train a classifier for seizure classification from EEG data; however, it can be considered that the performances obtained, in particular for KNN and RF, were similar to other studies when experiments were conducted under similar conditions, specifically the same classifier and EEG database.

In this sense, comparing our work with the state-of-the-art, in [12], the authors used the CHB-MIT database, seven feature selection methods, nine classifiers, and selecting the top 20 or 25 features, obtaining the best mean classification error of 0.12 by using the KNN classifier. Nonetheless, the authors tested their methodology neither using a wider spectrum of feature vector sizes, nor using different databases. In this work, we found that having the top 50 features, on the same database, the combinations KNN-SVM-50 (FS1), KNN-SVM-50 (FS2), and KNN-ERF-50 (FS2) (see Table 9) returned a mean classification error of 0.03. Moreover, the results in [12] correspond to a balanced dataset, while in this work, the proportion between seizure and no-seizure was 1:9, respectively, our methodology being more appropriate for seizure and no-seizure detection considering that epilepsy datasets are naturally unbalanced.

Additionally, Kathi and Ingle [25] presented an accuracy of 0.97 and an F1-score of 0.97 by using KNN and 11 feature metrics, testing their methodology on the Bonn University database [65] and using an equal proportion of healthy and seizure instances for the training and test sets. Furthermore, they reported a reduction in the feature set used to compute the feature vectors instead of directly reducing the size of the feature vectors. Hence, the authors did not evaluate either different sizes of the feature vectors or the different feature selection methods. On the contrary, in the present work, several experiments reported an accuracy of 0.97, but had a lower F1-score (0.84 and 0.85); however, these lower scores were expected considering that, for unbalanced datasets, the F1-score is more reliable.

On the other hand, RR presented the worst results no matter the classifier used, contrasting with [54], who reported RR as one of the FSMs that performed better and also provided good stability across datasets. However, the nature of the datasets used in [54] was different from EEG, which indicates that RR is not appropriate to be used for seizure and no-seizure detection.

Finally, concerning the comparison of selected features, by using the Jaccard index, it was observed that most of the time, SHAP-LIME returned the largest levels of similarity, meaning that both FMSs usually selected the same features, even though they did not present the best performance on the classification; on the contrary, SVM and ERF, which were the best FSMs for KNN and RF, respectively, did not present a higher Jaccard index.

Large values in SHAP-LIME were partially expected because, as explained in [49,50], SHAP and LIME are designed to be model-agnostic and to explain the classification model. In addition, the SHAP implementation applied in this study (kernel SHAP) followed a similar approach to LIME, but different functions were applied for the estimation of the

locality and the similarity between the classification *f* and interpretable *g* model. In the case of SVM and ERF, the low levels of the similarity of features selected were partially explained because the nature and training procedure of SVM are quite different from tree-based models.

Therefore, this experiment showed that two different feature spaces might result in a good classification performance, such as the features selected for KNN-SVM and RF-ERF (Figures 11 and 12). On the contrary, the RR method presented good similarity values compared to other FSMs; nonetheless, it presented the worst classification performances.

#### **5. Conclusions and Future Work**

In the present study, six FSMs were compared to define, first, the minimum number of features that can be chosen for each FSM without sacrificing the classifiers' performances, second, the performance of several C-FSM combinations to discover if a relationship exists between the FSMs with a particular classification algorithm, and third, if there is a feature or feature set that remains across different C-FSM combinations.

We can conclude that when the number of selected features was drastically reduced (100 and 50 features), the differences between classifiers' performance increased, but none of the FSMs showed a predominance over the rest. Furthermore, it was observed that it was possible to perform a large reduction of the number of features while having a low impact on the model performance until having a 50-feature vector length.

The results indicated that the classifiers' performance might be affected by diverse factors such as the EEG database, the features, and the number of features. However, the combinations KNN-SVM and RF-ERF are advised. Furthermore, the use of RR is not appropriate for seizure EEG data, as it yielded lower performances than the rest of the FSMs, and it was more time-consuming.

Regarding the proposed feature sets, FS2 is suggested to be used on seizure and no-seizure classification problems, given that the performance of several C-FSMs was improved while using it.

In future work, some opportunity areas could be explored: First, a more extended analysis is required to evaluate the FSMs in combination with deep learning models, including a more extensive parameter tuning process and the use of more complex features. Furthermore, considering that our results showed that RF-ERF obtained one of the best performances, it would be interesting to perform an evaluation of tree-based models for feature selection and classification of seizure and no-seizure EEG data.

Second, FSM stability was not evaluated in this study; the evaluation of stability with instance perturbation, as proposed by [48], would help to evaluate the robustness of an FSM against small variations in the EEG dataset, because variations can be caused by inter-subject variability or noise.

Finally, assessing a subject cross-validation methodology would be interesting for future work to test inter-subject performances.

**Author Contributions:** Conceptualization, S.E.S.-H.; methodology, S.E.S.-H., R.A.S.-R., S.T.-R., and I.R.-G.; software, S.E.S.-H.; validation, S.E.S.-H., R.A.S.-R., S.T.-R., and I.R.-G.; formal analysis, S.E.S.- H., R.A.S.-R., S.T.-R., and I.R.-G.; investigation, S.E.S.-H., R.A.S.-R., S.T.-R., and I.R.-G.; data curation, S.E.S.-H.; writing—original draft preparation, S.E.S.-H.; writing—review and editing, R.A.S.-R., S.T.- R., and I.R.-G.; visualization, S.E.S.-H., R.A.S.-R., S.T.-R., and I.R.-G.; supervision, R.A.S.-R., S.T.-R., and I.R.-G. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Publicly available datasets were analyzed in this study. The CHB-MIT Scalp EEG Database is available at https://physionet.org/content/chbmit/1.0.0/ (accessed on

18 October 2021) and the Siena Scalp EEG Database is available at https://physionet.org/content/ siena-scalp-eeg/1.0.0/ (accessed on 18 October 2021).

**Acknowledgments:** We appreciate the facilities provided to the project "Identificación de patrones en registros EEG de crisis epilépticas utilizando métodos de inteligencia artificial explicable" by the Data Analysis and Supercomputing Center (CADS, for its acronym in Spanish) of the University of Guadalajara through the use of the Leo Átrox Supercomputer.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

## *Article* **Wavelet Ridges in EEG Diagnostic Features Extraction: Epilepsy Long-Time Monitoring and Rehabilitation after Traumatic Brain Injury †**

**Yury Vladimirovich Obukhov 1, Ivan Andreevich Kershner 1,\*, Renata Alekseevna Tolmacheva 1, Mikhail Vladimirovich Sinkin 2,3 and Ludmila Alekseevna Zhavoronkova <sup>4</sup>**

**Abstract:** Interchannel EEG synchronization, as well as its violation, is an important diagnostic sign of a number of diseases. In particular, during an epileptic seizure, such synchronization occurs starting from some pairs of channels up to many pairs in a generalized seizure. Additionally, for example, after traumatic brain injury, the destruction of interneuronal connections occurs, which leads to a violation of interchannel synchronization when performing motor or cognitive tests. Within the framework of a unified approach to the analysis of interchannel EEG synchronization using the ridges of wavelet spectra, two problems were solved. First, the segmentation of the initial data of long-term monitoring of scalp EEG with various artifacts into fragments suspicious of epileptic seizures in order to reduce the total duration of the fragments analyzed by the doctor. Second, assessments of recovery after rehabilitation of cognitive functions in patients with moderate traumatic brain injury. In the first task, the initial EEG was segmented into fragments in which at least two channels were synchronized, and by the adaptive threshold method into fragments with a high value of the EEG power spectral density. Overlapping in time synchronized fragments with fragments of high spectral power density was determined. As a result, the total duration of the fragments for analysis by the doctor was reduced by more than 60 times. In the second task, the network of phase-related EEG channels was determined during the cognitive test before and after rehabilitation. Calculation-logical and spatial-pattern cognitive tests were used. The positive dynamics of rehabilitation was determined during the initialization of interhemispheric connections and connections in the frontal cortex of the brain.

**Keywords:** electroencephalogram; wavelet spectrum; ridge; segmentation; phase connectivity; epilepsy; traumatic brain injury

#### **1. Introduction**

Wavelet transform (WT) is widely used in the processing and analysis of non-stationary signals [1–5]. Since the 1990s, in various fields of biology and medicine [6], in neurophysiology [7], discrete and continuous wavelet transforms have been used to extract diagnostic information from signals and images of various types.

**Citation:** Obukhov, Y.V.; Kershner, I.A.; Tolmacheva, R.A.; Sinkin, M.V.; Zhavoronkova, L.A. Wavelet Ridges in EEG Diagnostic Features Extraction: Epilepsy Long-Time Monitoring and Rehabilitation after Traumatic Brain Injury. *Sensors* **2021**, *21*, 5989. https://doi.org/10.3390/ s21185989

Academic Editor: Yvonne Tran

Received: 13 August 2021 Accepted: 3 September 2021 Published: 7 September 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

Considering EEG as a simultaneously amplitude and phase modulated analytical signal, and if the scanning wavelet width is narrower than the changes in signal phase, then it is possible to use the property of the wavelet spectrum ridge, namely those that the amplitude and phase of the signal are equal to the amplitude and phase of the wavelet spectrum ridge [8–11]. Thus, defining the ridge as the absolute maximum of the wavelet spectrum at each moment (reference point) of time, we obtain instantaneous values of the amplitude, frequency and phase of the signal. This very useful property of WT ridges makes it easy to find interchannel synchronized EEG fragments during epileptic seizures (ESs) in long-term clinical monitoring data, restoring the cognitive functions of patients after moderate traumatic brain injury using interchannel phase coupling link analysis, and other tasks of EEG diagnostics.

WT is used for EEG decomposition into time-frequency fragments for the subsequent detection of epileptic seizures (ESs). Currently, there are many publications on the use of various classifiers for the detection and prediction of an epileptic seizure in EEG signals using various classifiers [12–19]. Initial data on epilepsy monitoring should be preliminarily processed, including removal of artifacts and filtering noise to get a clean epilepsy EEG signal for the next step, feature extraction and classification [18,19].

In decision support systems, methods based on the analysis of EEG patterns are most often used and one of them is the "Persyst" system by Persyst Development Corpo ation (https://www.persyst.com, accessed on 13 August 2021). To detect ES in the time domain, discrete-time sequences are analyzed into which the original EEG signal is divided. One of such methods is based on tracking successive extrema in the selected time interval of the signal and evaluating the histogram of the amplitude difference and time separation between the maximum and minimum values of the histogram [20]. The different approaches for detecting ESs which were proposed in the time domain are the calculation of signal energy [21]; the frequency characteristics of the signal were studied: the index of the phase slope of multichannel EEG [22]; frequency-moment signatures [23]; entropy features [24]; Bayesian linear discriminant analyses of lacunarity and fluctuation index [25], four-level Daubechies wavelet transform [26] and five-level wavelet decomposition method [27]. The most promising method of EEG analysis is the study of the parameters of the ridges of wavelet spectrograms. In the EEG signals for the detection of epileptic seizures, the dynamics of synchronization and changes in the phase ratio before, during and after the seizures are monitored.

At present, attempts are being made to improve methods for detecting ESs in EEG. Paper [28] describes a way to improve the support vector machine method by adding an adaptive median feature baseline correction method. A combination of methods is also used to search for ESs, for example, complementary ensemble empirical mode decomposition with extreme gradient boosting [15]. A method has been proposed that combines timedomain feature analysis and entropy calculation [16]. A similar combination was also presented in work [14], but the study of parameters in the time domain was used to segment the signal sections, in order to then carry out analysis using machine learning methods. To differentiate ESs from non-seizure events, neural networks [13] and similar methods are used, such as the method of binarization of frequency and temporal features of signal fragments [29].

It should be noted that the estimation of the accuracy and specificity of the classification was carried out on EEG fragments previously selected and annotated by EEG neurosciensists as ictal and interictal events. The most representative databases are the EPILEPSIAE database [30], the Temple University Hospital EEG Data Corpus [31], Bonn epilepsy dataset etc.

Clinical EEG investigations of epilepsy consist in long-term (several days) monitoring of multichannel EEG using scalp or intracrinial electrodes in the presence of various artifacts: the electrical activity that is not recorded in the cerebral zone, such as that due to the equipment, patient behavior or the environment; eye movement and chewing are common events that can often be confused with a spike; signals instrument fluctuations and artifacts of vital activity [32]. It can be seen from the review that the methods for removing artifacts described in the literature are mainly focused on removing one type of possible artifacts or ocular and muscular ones present in the real initial data of long-term EEG monitoring. In general, we can conclude that the problem of automated removal of artifacts of various types from the initial data of long-term EEG monitoring has not been fully resolved. Additionally, this article describes canonical correlation analysis as a successful method for removing muscle artifacts.

One of the most important characteristics of ES is abnormal inter-channel synchronization or so-called coherency. To assess interchannel EEG synchronization, canonical correlation analysis [33], normalized cross-correlation and imaginary part of coherency or phase synchronization are used [34]. The main disadvantage of estimation coherence is the necessity to average it over time epochs and frequency ranges [35]. The study of short-term frequency synchronization of signals in two EEG channels by comparing their WT ridges frequencies during a previously selected by physicians is presented in ES [36].

We did not find in the literature any information on taking into account one of the most important feature of ESs—EEG interchannel synchronization for detecting ESs. So this article describes a new approach to the segmentation of the initial long-term clinical multichannel EEG monitoring data of patients with epilepsy into temporal fragments suspicious of an ES, to reduce the quantity of EEG fragments. This approach is based at first on the EEG WT ridges segmentation of the into frequency-synchronized fragments, and secondly with a thresholding of the ridge spectral power density.

Another part of this paper is devoted to a new approach to the diagnosis and assessment of rehabilitation of patients after traumatic brain injury (TBI). TBI is an insult to the brain from an external mechanical force, which can lead to permanent or temporary impairment of cognitive, physical, and psychosocial functions. The most used EEG methods of investigation TBI are spectral analysis, absolute and relative amplitude and power, coherence, and symmetry between homologous pairs of electrodes (see review [36]). A multivariate support system has been developed to quantify and classify by Random Forest classifier TBI stage based on analysis of EEG power in various frequency ranges [37]. A study [38] investigated the possibility of detecting moderate TBI according to the Glasgow Coma Scale [39] by EEG amplitude analysis and convolutional neural network classification. Recently, a single channel system was developed for real-time mild TBI detection with Convolutional Neural Networks classifier of EEG power in different frequency ranges [40]. The proposed method can be applied for screening of the moderate TBI and for selection of the patients for further diagnostics and treatment. In [41] the analysis of the EEG data applying the energy, sample entropy, approximate entropy, Lempel–Ziv complexity features demonstrated the increase in sample entropy was related with the functional recovery, i.e., the rehabilitation dynamics of the injured brain region. EEG-based neurofeedback is used for cognitive rehabilitation of patients with TBI [42].

Our approach is based on the analysis of the network of phase sinchronized EEG channels WT ridges in patients with moderate TBI. The interchannel phase difference of the EEG is determined during cognitive tests at the points of the frequency-modulated wavelet spectra ridges. We investigate the neurons connectivity disruption of the brain after TBI and consider the inter-channel phase connectivity between EEG channels during cognitive tests. It does not depend on the EEG signal amplitude. In this paper Section 2 contains the basic formulae and conditions for their application. Section 3 describes a new approach to segmentation of long-term EEGs into temporal fragments suspicious of an ES by interchannel WT ridges frequency sincronization and power spectral density thresholding. Section 4 describes a new approach to determine the evaluation of rehabilitation positive dynamics of patients with moderate TBI.

#### **2. Materials and Methods**

We studied long-term (from several hours to several days) initial EEG records of preoperative patients with epilepsy, obtained in laboratory of invasive neurointerfaces

of the Research Institute TechnoBioMed. A.I. Yevdokimov Moscow State University of Medicine and Dentistry. The segmentation method was used for several days' 19-channel EEG. The records were carried out according to the 10–20 system [43] in reference montage with a sampling rate of 256 Hz. Power supply artifacts were removed from all EEG channels using a notch filter at frequencies multiples of 50 Hz. The use of Morlet WT in the frequency range from 0.5 to 22 Hz, so myographic artifacts were rejected.

The records of 19-channel EEG were considered, therefore the quantity of pairs of channels is 171 for the group of control volunteers (18 subjects) and for the group of patients with moderate TBI (12 subjects), where three patients had repeated EEG records after rehabilitation in two cognitive tests. Cognitive tests were calculation-logical (CT1) and spatial-pattern (CT2). During the CT1 test, the doctor randomly spoke words from the category of "clothing" or "food" to the subject. During the test, the subject counted in their mind the number of items belonging to one of these categories. At the end of the test, the subject announced the result of the number of items. On the CT2 test, the doctor named an arbitrary time. The subject had to represent the position of the hands-on-dial in accordance with the indicated time. If both clock hands were in the same half of the dial, he said "yes", and if they were in different halves, he kept silent. Investigations of control volunteers and patients with moderate TBI were carried out at the National Medical Research Center for Neurosurgery named after Academician N.N. Burdenko. All subjects were right-handed and signed written consent to participate in the research in accordance with the provisions of the Helsinki Agreement. The rehabilitation was performed for 1–2 months. The time of the rehabilitation was 40–45 min two times a week. The criteria for the inclusion of patients in the investigation were the ability to stand on their own and the ability to follow the doctor's instructions, and also the absence of hemiparesis and other neurological disorders. The international 10–20 system of the position of scalp electrodes was used for EEG record. The recording time for every test was 60 s. EEG recording was carried out both during the tests and without them. The sampling rate of the EEG was 250 Hz in the processing of EEG signals. The original signals were recorded with a high-pass filter with a cut-off frequency of 0.5 Hz, a low pass filter with a cut-off frequency of 70 Hz. Then, a notch filter at frequencies multiples of 50 Hz and a Butterworth filter were used. The signals were filtered by a fourth-order Butterworth bandpass filter with a bandwidth from 2 to 10 Hz. The EEG records were analyzed without selecting individual fragments of the signal. However, the removal of outliers in the EEG signals was done with the Huber's X84 method [44].

We considered EEG as an analytical signal with time-varying amplitude and frequency. The analytic signal was first locally represented as a modulated oscillation, demodulated by its own instantaneous frequency, and then Taylor-expanded at each point in time. We represent this signal as the following function:

$$S(t) = A\_S(t) \exp(i\Phi\_S(t)),\tag{1}$$

where *AS*(*t*) is the amplitude and Φ*S*(*t*) is the phase of the signal. Continuous wavelet transform of signal *S*(*t*) is represented as:

$$\mathcal{W}(a,b) = \int\_{-\infty}^{\infty} S(t) \psi\_{a,b}^\*(t) \, dt \tag{2}$$

$$
\psi\_{a,b}(t) = \frac{1}{\sqrt{|a|}} \psi\left(\frac{t-b}{a}\right) \tag{3}
$$

where *a*, *b*, *a* = 0 are the real numbers defining the scale and the shift. We used the following function (Morlet mother wavelet) that was employed in the Matlab software:

$$\psi(t) = \frac{1}{\sqrt{\pi f\_b}} \exp(-t^2/f\_b) \exp(2\pi i f\_c t) \tag{4}$$

where *fb* is a positive and related with the variance of Gaussian function and *fc* is a positive value that corresponds with central frequency. The Morlet wavelet transform can be represented as follows:

$$\mathcal{W}(a,b) = M(a,b) \exp(i\Phi(a,b)),\tag{5}$$

where *M*(*a*, *b*) is the absolute value of wavelet transform and Φ(*a*, *b*) is the phase of wavelet transform (2).

Substituting expressions (3) and (4) in formula (2), we obtain:

$$\mathcal{W}(a,b) = \frac{1}{\sqrt{a\pi f\_b}} \int\_{-\infty}^{\infty} A\_S(t) \exp\left(-\frac{(t-b)^2}{a^2 f\_b}\right) \exp\left(i\left(\Phi\_S(t) - 2\pi f\_c \frac{t-b}{a}\right)\right) dt,\tag{6}$$

usually (in Matlab) *fb* = *fc* = 1.

Integral (6) is approximately calculated by the method of stationary phase [45]. Under certain conditions, the main contribution to the integral is made by the imaginary part of the exponential function, since the contributions of rapidly changing phases cancel each other out, and the contribution is made by the values located at the point of the stationary phase. The stationary phase method is applicable when the amplitude *A*(*t*) of the signal exhibits relatively slow changes compared to fast changes in the total signal associated with fast changes in phase, for example, and asymptotic properties are satisfied concerning the window *ψ*(*t*) under the following assumptions; so that the following conditions are satisfied [11]:

$$\left|\frac{d\Phi\_S(t)}{dt}\right| \gg \left|\frac{1}{A\_S(t)}\frac{dA\_S(t)}{dt}\right|, \left|\frac{1}{A\_S(t)}\frac{dA\_S(t)}{dt}\right| \ll \left|\frac{1}{\psi(t)}\frac{d\psi(t)}{dt}\right|\tag{7}$$

The relationship between the phase from expression (5) and the phase from expression (6) is given as follows:

$$
\Phi(t) = \Phi\_S(t) - 2\pi f\_c \left(\frac{t-b}{a}\right) \tag{8}
$$

For the stationary phase Φ(*t*), we have

$$\frac{d\Phi(t)}{dt} = \Phi\_S'(t) - \frac{2\pi f\_c}{a} = 0\tag{9}$$

Such a condition is satisfied at *t* = *t*(*a*). To estimate the integral from formula (6), we expanded the phase Φ(*t*) in a Taylor series up to a polynomial of the second degree in the neighborhood of point *<sup>t</sup>* = *<sup>t</sup>*(*a*) till the order (*<sup>t</sup>* − *<sup>t</sup>*(*a*))2:

$$\Phi(t) \approx \Phi\_S(t(a)) - 2\pi f\_c \left(\frac{t(a) - b}{a}\right) + \frac{1}{2} \Phi\_S^{\prime\prime}(t(a))(t - t(a))^2 \tag{10}$$

Below, we use the notation Φ*<sup>S</sup>* ≡ Φ*S*(*t*(*a*)) and Φ *<sup>S</sup>* ≡ Φ *<sup>S</sup>*(*t*(*a*)).

After substituting formula (10) into formula (6), we obtained an approximate value for the phase and absolute value of the wavelet transform:

$$\Phi(a,b) \approx \Phi\_S - 2\pi f\_c \left(\frac{t(a) - b}{a}\right) + \frac{1}{2} \arctan\left(\frac{a^2 f\_b}{2} \Phi\_S^{\prime\prime}\right) + \frac{2(t(a) - b)^2 \Phi\_S^{\prime\prime}}{4 + a^4 f\_b^2 (\Phi\_S^{\prime\prime})^2} \tag{11}$$

$$M(a,b) \approx A\_S(t(a)) \left(1 + \frac{a^4 f\_b^2}{4} (\Phi\_S^{\prime\prime} t(a))^2\right)^{-\frac{1}{4}} \exp\left(-\frac{(t(a)-b)^2 a^2 f\_b(\Phi\_S^{\prime\prime})^2}{4 + a^4 f\_b^2 (\Phi\_S^{\prime\prime})^2}\right) \tag{12}$$

Expression (12) shows that the maximum of wavelet transform absolute value was reached at *b* = *t*(*a*). The instantaneous frequency at the ridge point *fr* at time moment *t* = *t*(*a*) was calculated using expression (9):

$$f\_r(t(a)) = 2\pi \frac{f\_c}{a} \tag{13}$$

In this case, the maximal value of the wavelet transform (ridge) of the signal is given by

$$\max\_{a} |\mathcal{W}(a, b)| \approx A\_S(t(a)) \left( 1 + \frac{a^4 f\_b^2}{4} (\Phi\_S^{\prime\prime} t(a))^2 \right)^{-\frac{1}{4}} \tag{14}$$

and the phase is approximated of a ridge point as

$$\Phi(a,b) \approx \Phi\_S + \frac{1}{2} \arctan\left(\frac{a^2 f\_b}{2} \Phi\_S^{\prime\prime}\right) \tag{15}$$

As in [46], the relationship between the frequency of Fourier spectrum of the wavelet transform and the scales *a* of the wavelet transform (2) is given as follows:

$$f = \frac{f\_0}{2a} + \frac{\sqrt{2 + 4(\pi f\_0)^2}}{4\pi a} \cong \frac{f\_0}{a} = \frac{1}{a} \tag{16}$$

where *<sup>f</sup>*<sup>0</sup> is a wavelet central frequency and it is considered that 4(*<sup>π</sup> <sup>f</sup>*0)<sup>2</sup> 2. So, for the ridge points (*fr*, *t*) we have:

$$\begin{aligned} \mathcal{W}\_r(t) &= \max\_f |\mathcal{W}(f,t)|, \; f\_r(t) = \underset{f}{\text{argmax}} |\mathcal{W}(f,t)|,\\ \Phi\_S(t) &\cong \Phi\_r(f,t) = \arctan\left(\frac{\operatorname{Im}(\mathcal{W}(t,f\_r))}{\operatorname{Re}(\mathcal{W}(t,f\_r))}\right), \end{aligned} \tag{17}$$

when the condition

$$\frac{\Phi\_S^{\prime\prime}}{2f\_r^2} = \frac{f\_r^{\prime}}{2f\_r^2} \ll 1,\tag{18}$$

is satisfied.

Summarizing, it should be noted that, in contrast to other works, the obtained simple method for determining the ridge points as the maximum of the modulus of the wavelet spectrum at each time point was undoubtedly easy to calculate.

#### **3. EEG Segmentation**

This chapter describes an EEG segmentation method based on the study of Morlet wavelet transform ridges, which allows finding time intervals of interest in ES detection, which is used to analyze continuous long-term EEG monitoring data during postprocessing. The long-term EEG segmentation method consists of the following stages, as shown in Figure 1: 1. signal filtration at frequencies multiples 50 Hz, 2. wavelet Morlet transform of signals, 3. determination of wavelet spectrogram ridges, 4. marking time intervals with interchannel synchronization, 5. marking time intervals with power spectral density (PSD) values above the threshold, 6. intersections of time intervals, 7. visualization of a segmented signal with marked time intervals.

**Figure 1.** Block diagram of a long-term EEG segmentation.

Let us show how segmentation was carried out using the example of an EEG recording fragment containing an ES. For each EEG channel, we calculated the wavelet spectrogram (2) and the ridges of the wavelet spectrogram (17) in frequency range [0.5; 22] Hz.

Generalized ESs were characterized by changes in power in several EEG channels and the synchronization of different channels pairs. In order to estimate the inter-channel synchronization, the modulus of the frequency difference at the points of the ridges was calculated for each pair of channels. If the modulus of the difference was less than *ε*, then there was synchronization *Synci*,*j*, otherwise, it was not:

$$Sym\_{i,j}(k) = \begin{cases} 1, |f\_{ri}(k) - f\_{rj}(k)| \le \varepsilon \\ 0, |f\_{ri}(k) - f\_{rj}(k)| > \varepsilon \end{cases} \tag{19}$$

where *fri*, *frj* are the frequencies of the ridges of the wavelet spectrograms on the *i* and *j* EEG channels, *k* is the point of the ridge.

Figure 2 shows the projection of the wavelet spectrogram onto the PSD-frequency plane of the sinusoidal signals with a frequency of 2 and 2.5 Hz. At *ε* = 0.5 the peaks were distinguishable. On smaller epsilons, the peaks could merge.

**Figure 2.** PSD-Frequency projection of wavelet spectrums of two sinusoidal signals: blue is 2 Hz, orange is 2.5 Hz.

Nearby points at which condition (19) was satisfied were combined into fragments. Fragments between which the time interval was less than 10 s were combined into one. For each fragment, the beginning and end times of synchronization in pairs of channels were calculated. Table 1 shows a histogram of the number of synchronized fragments depending on the duration in 19 pairs of EEG channels. For neurophysiological considerations, this work considered fragments with a duration of 10 s or more.

**Table 1.** Histogram of the number of synchronized fragments depending on the duration in 19 pairs of EEG derivations.

Figure 3 shows an example of ES fragment from observed EEG recording with visualization of the presence of synchronization in channels pairs.

**Figure 3.** Fragment of daily EEG signal with ES illustrating the frequency of the wavelet spectrogram ridge synchronization in different pairs of EEG channels. Black shows the presence of synchronization. The ordinate shows the labels of channel pairs.

As can be seen, EEG synchronization could be observed in not all channels simultaneously (Figure 3). In this example, synchronization was observed in most channel pairs from about 5970 s to 6000 s, but synchronization began to appear earlier in a smaller number of pairs.

The time intervals in which the inter-channel synchronization in the frequency of the ridges was recorded could correspond to both ES and artifacts of chewing, sleep, and random physical influences on the electrodes, which generated artifacts of a non-epileptic nature.

A characteristic feature of the ES was a sharp change in the amplitude over a short period of time. Therefore, in addition to searching for time intervals in which there was synchronization on several pairs of EEG derivations, the detection of areas with high values of the power spectral density (PSD) was carried out. In order to understand the idea of the method, let us consider the histogram of the ridge points of the wavelet spectrogram, calculated for one of the leads. The peak of the histogram contained about 1.2 × 106 points, the maximum PSD value at which the number of ridge points tended to 0, about 2.8 × <sup>10</sup>7<sup>μ</sup> V2/Hz. Such a histogram gave a large peak in the region of low PSD values and did not allow us to estimate the distribution of the ridge points, therefore, Figure 4 shows the "window" of the PSD histogram. As can be seen from the figure, the number of ridge points with low PSD values was large and could be interpreted as noise. It was necessary to separate the informative points of the ridge from the noise.

In order to separate the points of the ridge of the wavelet spectrogram related to high-amplitude electrical activity from noise, it was required to find the threshold value of PSD *Tr*. The ridge *PSDr* values of the wavelet spectrogram were determined as follows:

$$PSD\_r(t) = \begin{cases} PSD\_r(t), PSD\_r(t) \ge Tr \\ 0, PSD\_r(t) < Tr \end{cases} \tag{20}$$

The points of the ridge *PSDr*(*t*) lying between the nearest points *PSDr*(*t*) = 0 was called the ridge segment. Figure 5 shows a histogram of the number of ridge segments from the PSD threshold *Tr* to *Tr* = <sup>5</sup> × <sup>10</sup>5<sup>μ</sup> <sup>V</sup>2/Hz. At large values of *PSDr*, the number of segments tended to be 0. The *PSDr* values could differ greatly not only from patient to patient but also by channels; therefore, it was required to determine the threshold value adaptively. Figure 5 shows a sharp decrease in the number of fragments with an increase in the threshold value. To identify the threshold value, the second derivative of segments

quantity from the threshold changes was analyzed. After reaching a local maximum (circle mark in Figure 5), it became negligible. This means that the segments quantity linearly decreased with growing threshold. In Figure <sup>5</sup> threshold PSD value *Tr* = 1.5 × <sup>10</sup>5<sup>μ</sup> V2/Hz. With this choice of the threshold, most ESs detected by the expert and a small number of artifacts like ES were observed. Figure 6 shows an example of a segmented ridge of a wavelet spectrum containing an ES.

**Figure 4.** The window of the histogram of the PSD of the wavelet spectrogram ridge of the long-term EEG signal.

**Figure 5.** Histogram of the number of ridge segments from the threshold PSD *Tr*. The circle marks the local maxima at threshold value *Tr* <sup>=</sup> 1.27 <sup>×</sup> 105<sup>μ</sup> <sup>V</sup>2/Hz.

Figure 7 shows a fragment of a daily EEG signal showing an ES. Marks of the expert neurophysiologist are green vertical lines; the blue line, repeating the waveform, marks the fragment on which the synchronization was recorded on several pairs of EEG derivations; the dotted rectangles mark the areas found by the threshold method. Thus, the method of application for the search for ES is shown.

For the 5-hour EEG 2017 the overall synchronized fragments duty more than 10 were found (see Table 1), 112 segments were detected by thresholding, and finally we obtained nine intersected segments. with total duration total duration about 4 min. Earlier [47], we showed that by processing synchronous video of these nine fragments four fragments were recognized as moving artefacts. As a result, there were five segments left with a total duration of 4 min.

**Figure 6.** The segmented wavelet spectrogram ridge of the EEG signal fragment typical for ES. The upper figure shows a segmented ridge *PSDr*, the bottom one shows *fr*.

**Figure 7.** Fragment of an EEG with ES. Green vertical lines are expert marks. The blue line repeating the signal waveform is a mark obtained by searching for synchronized channel pairs. Dotted squares are marks obtained by the threshold method for detecting ES.

The detection of ES on EEG is complicated by the presence of many non-epileptic artifacts in signals received from scalp electrodes: electromyographic, motor, instrumental human actions, etc. An overview of various types of artifacts is given in [48]. Therefore, there is a need to develop methods to differentiate ES from artifacts of a non-epileptic nature.

To solve this problem, an algorithm was proposed, which consisted of studying the broadband peaks of the wavelet spectrograms, which were characteristic of an ES and a chewing artifact. Let us make a comparison using the example of wavelet spectra of an epileptic seizure (Figure 8) and chewing (Figure 9). We analyzed slices of wavelet spectrograms frequency *fcur*(*t*) higher ridge frequency *fr*(*t*), for example, at Figures 8 and 9 *fcur*(*t*) = 4 Hz (green line).

For each slice, we calculated Fourier spectra. Figure 10 Fourier spectra of ES and chewing artifact at *fcur* = 4 Hz. The frequency of the main peak and full width at half maximum (FWHM) of the Fourier spectrum were calculated.

**Figure 8.** Wavelet spectrogram of an EEG with ES. Green line corresponds to slice of the wavelet spectrogram at the frequency *fcur* = 4 Hz.

**Figure 9.** Wavelet spectrogram of an EEG with chewing artifact. Green line corresponds to slice of the wavelet spectrogram at the frequency *fcur* = 4 Hz.

**Figure 10.** Fourier spectra of wavelet spectrograms slices *fcur* = 4 Hz. The red line is the Fourier Spectrum of the ES slice; the blue line is the Fourier spectrum of the chewing artifact slice.

Figure 10 shows differences between ES and chewing artifacts. The main peak frequency of the Fourier spectrum was at 0.71 Hz for the chewing artifact, and at 1.86 Hz for the ES. Two peaks could be observed in ES slice spectrum, this can be interpreted as the presence of spike-wave activity. There was a difference in FWHM of the peaks of the Fourier spectrum of the slices: for a chewing artifact, it was almost 2 times more than for an epileptic seizure. This may mean that the seizure period was more stable than chewing.

#### **4. The Estimation of Inter-Channel EEG Phase Connectivity in Patients with TBI**

Various methods of EEG phase coherence are used to estimate the connectivity of brain regions. Usually, the phase coherency of signals is used for the estimation of the interchannel connectivity of EEG [33,49,50]. Coherency *Cohxy*(*f*) is defined by the normalized complex cross-correlation *Cxy* of signals *x*(*t*) and *y*(*t*):

$$\text{Col}\_{xy}(f) = |<\mathcal{C}\_{xy}>|\mathcal{C}\_{xy} = \frac{S\_{xy}(f)}{\left(|S\_{xx}(f)| |S\_{yy}(f)|\right)^{1/2}},\tag{21}$$

and a phase coherency is defined as | < *exp*(*i*ΔΦ) > |, where | < • > | is an averaging [33].

In coherency analysis of non-stationary EEG, it is necessary to average exp(*i*ΔΦ) over different time intervals (epochs), and it is the first problem. The presence of the peak in the histogram of the phase difference in different epochs determines the presence or the absence of the phase synchronization in the absence of a peak. In addition to this, *Cohxy*(*f*) is averaged in preliminary selected frequency bands that are specified using neurophysiological data. Usually, these bands correspond to the delta (2–4 Hz), theta (4–8 Hz), alpha (8–12 Hz), beta (12–25 Hz) EEG, and other rhythms, and this is the second problem. These disadvantages of the coherency analysis that leads to instability in the definition of the inter-channel EEG connectivity. The validity of the coherent analysis of non-stationary EEG signals is questioned [34].

Another method for the estimation of the phase connectivity is to determinate the analytical signal *x*∗(*t*) = *x*(*t*) + *iH*(*x*(*t*)), where *H*(*x*(*t*)) is the Hilbert transform [51]. Then, the phase of signal *x*∗(*t*) is calculated as the arccosine (arcsine) of the ratio of the real (imaginary) part *x*∗(*t*) to its modulus [52]. The phase synchronization of two signals takes place when:

$$|\Delta \Phi\_{x,y}(t)| \le const,\tag{22}$$

where ΔΦ*x*,*y*(*t*) = *n*Φ*x*(*t*) − *m*Φ*y*(*t*), Φ is a phase of the signal; *n*, *m* are integers. Then, the angular frequency of the signal can be found by the phase differentiating with respect to time. Numerical differentiation in the presence of phase fluctuations is an unstable procedure. Additionally, the disadvantage of the approach associated with the calculation of analytical signals is that it is well applicable for narrowband signals and not good enough for broadband signals [53].

The paper describes the methods and results of determining the phase-connected pairs of EEG channels of patients with moderate TBI before and after rehabilitation, which can be used to estimate the dynamics of treatment and rehabilitation of patients. The method of the estimation of the inter-channel EEG phase synchronization at the points of the ridges *fr*(*ti*) of their wavelet-spectrograms (6) is considered as an inverse task for the task of modeling ridges:

$$f\_r(t\_i) = \arg\left\{ \max\_{f(t\_i)\in\{1:25\text{ Hz}\}} \left( |\mathcal{W}(t\_i, f(t\_i))| \right) \right\},\tag{23}$$

on the condition (18). <sup>Φ</sup>*x*(*t*) <sup>∼</sup><sup>=</sup> <sup>Φ</sup>*r*(*<sup>f</sup>* , *<sup>t</sup>*) = arctan *Im*(*W*(*t*, *fr*)) *Re*(*W*(*t*, *fr*))  according to (17), when the condition (18) is satisfied.

However, the ridge *Wr*(*t*) can be considered as the frequency-modulated signal. It is necessary to take an unmodulated oscillation [54]:

$$\mathbf{x} = A\_0 \sin(\omega\_0 t + \Phi\_0),\tag{24}$$

and enter a variable frequency *ω* = *ω*<sup>0</sup> + Δ*ωξ*(*t*) = *ω*(*t*), where *ξ*(*t*) is some unknown function, and *ω*(*t*) is known.

Then if Φ<sup>0</sup> = 0:

$$\mathbf{x} = A\_0 \sin \left( \omega\_0 t + \Lambda \omega \int\_0^t \xi(t) \, dt \right) = A\_0 \sin(\omega(t)t), \tag{25}$$

Then it is possible to estimate the phase of the ridge as [55]:

$$
\Phi(t, f\_r) = 2\pi f\_r(t)t,\tag{26}
$$

Figure 11 represents two ridge frequencies of the Morlet wavelet transform for two EEG channels. Ridge points are points of the maximum power spectral density. Fp1 EEG channel is indicated by the blue line. Fp2 EEG channel is indicated by the red line. The abscissa is the time in seconds; the ordinate is the frequency in Hz.

**Figure 11.** Ridge frequencies of the Morlet wavelet transform for two EEG channels. Fp1 EEG channel is indicated by the blue line. Fp2 EEG channel is indicated by the red line.

The EEG frequencies coincided in some time fragments. The phase of the ridge could be estimated with the formula (26) if the ridge frequency was known.

The phases of the EEG signals were calculated and compared at the points of the ridges (*ti*, *fr*) of their wavelet spectrograms in EEG records both with cognitive tests and without tests. Then, the phase difference of two signals *x*(*t*) and *y*(*t*) in two EEG channels was calculated. Next, the normalized histogram of portions *ρx*,*<sup>y</sup>* = *nx*,*y*/*N* in different pairs of EEG channels was calculated, where *nx*,*<sup>y</sup>* is the quantity of reference points of ridges with |ΔΦ*x*,*y*(*t*)| < 0.01*π*, *N* is a total quantity of EEG signal reference points in the test.

Figure 12 represents the normalized histograms of portions of the phase difference at the ridge points of the wavelet spectrograms of two EEG channels for the case of a phase coupled pair of EEG channels Fp1-Fp2, which were obtained by two methods. Figure 12a demonstrates the first way based on the calculation of the phase according to the formula (26). Figure 12b demonstrates the second way based on the calculation of the phase according to the formula (17).

**Figure 12.** Histograms of portions *ρx*,*y* for the phase difference at ridge points for two EEG channels. This example shows a phase-connected pair of Fp1-Fp2 EEG channels. (**a**) The histogram is obtained from the phase calculation by (26). (**b**) The histogram is obtained from the phase calculation by (17).

Figure 12a demonstrates that the histogram of portions of the phase difference at the points of the ridge of the wavelet spectrograms calculated by the first way (26) had a higher and sharper peak versus the second way (17) (Figure 12b). Below, we will calculate the phase by formula (26).

Let *A* = max(*ρx*,*y*) be the maximum values of the histogram in the cognitive test and let *B* = max(*ρx*,*y*) are the maximum values of the histogram in the EEG record without a test. It is convenient to consider the difference *D* = *A* − *B*, which was sorted in order to increase max(*ρx*,*y*). Figure 13 demonstrates the dependence of *D* sorted in increasing order versus the numbers of a pair of EEG channels and its derivative for a healthy subject in the CT1 test.

**Figure 13.** The dependence of *D* sorted in increasing order (line1) versus the numbers of a pair of EEG channels and its derivative (line 2) for a control subject in the CT1 test.

Figure 13 shows that the curve of the graph appears at some point *D*. It is advisable to consider pairs of channels with numbers greater than at point sharp of increasing of derivative *D* (black point) as phase-connected pairs. Thus, phase-connected pairs of EEG channels were identified before and after rehabilitation of patients with moderate TBI.

Figure 14 shows a block diagram of the developed algorithm for the determination of phase-connected EEG channels. The developed algorithm for the determination of phase-connected EEG channels consisted of the following stages, as shown in Figure 14: 1. Preprocessing of signals. It was outlier removing; notch filter at frequencies multiples of 50 Hz; filtering of signals with a Butterworth filter; 2. Calculation of wavelet spectra and ridges; 3. Calculation of the ridges phase at each point of the wavelet spectra ridges. Calculation histograms of the phase difference portions (*ρx*,*y*) in two channels for 171 channels pairs. The determination max(*ρx*,*y*) for each channel pairs; 4. Calculation of the difference between max(*ρx*,*y*) with cognitive test and without a test (*D*), sorting in in-creasing order of *D*. Calculation derivative (*D*); 5. The determination of phase connected EEG channels. If the derivative (*D*) sharply increased with the growing pair, the pairs with numbers greater than at the sharp point of the increasing derivative (*D*) were considered as a phase connected pairs. If the derivative (*D*) did not sharply increase, it was impossible to identify phase connected pairs.

**Figure 14.** The block diagram of the developed algorithm for the determination of phase-connected EEG channels.

Figure 15 demonstrates the phase-coupled pairs of EEG channels for seven healthy subjects during the EEG recording in the CT1 test.

**Figure 15.** Phase-connected pairs of EEG channels of control subjects during the EEG recording in the CT1 test.

Figure 15 represents that the frontal regions and interhemispheric connections are activated in cognitive tests (CT1). Interhemispheric connections and connections in the frontal cortex in control subjects are activated during CT1 test in accordance with published work [56]. However, each control subject and patient with TBI are characterized by different phase-connected pairs due to the individuality of each person during the CT1 test. Therefore, we considered phase-connected pairs individually for each subject.

Figure 16 demonstrates the phase-connected pairs of EEG channels for the seven control subjects during the EEG recording in the CT2 test.

**Figure 16.** Phase-connected pairs of EEG channels of the control subjects during the EEG recording in the CT2 test.

Figure 16 represents that the frontal regions and interhemispheric connections were activated in cognitive tests (CT2). Interhemispheric connections and connections in the frontal cortex in control subjects were activated during CT2 test in accordance with published work [56]. However, each control subject and patient with TBI were characterized by different phase-connected pairs due to the individuality of each person during the CT2 test. Therefore, we considered phase-connected pairs individually for each subject.

Figure 17 demonstrates phase-connected pairs of EEG channels for three patients with TBI during the EEG recording in CT1 and CT2 tests.

**Figure 17.** Phase-connected pairs of EEG channels of patients with TBI during the EEG recording in CT1 and CT2 tests.

Figure 17 shows that phase-connected pairs appeared more in the parietal and occipital regions, than in the interhemispheric and frontal cortex in patients with TBI during CT1 and CT2 tests.

Additionally, the dynamics of inter-channel EEG synchronization of three patients with TBI before and after the rehabilitation was also investigated. The phase-connected EEG pairs in patients before and after rehabilitation were compared with the phase-connected pairs of the control group for each test. If interhemispheric connections or connections in the frontal cortex were activated in patients, as in control subjects in cognitive tests (CT1 and CT2), it could be concluded that the cognitive function had positive dynamics.

Figure 18 demonstrates that the positive dynamics could be seen of the rehabilitation of a patient with TBI in the CT1 test. If interhemispheric connections or connections in the frontal cortex in the CT1 test appeared after rehabilitation, as in the control subjects, the positive dynamics of rehabilitation could be concluded.

**Figure 18.** Phase-connected pairs of EEG channels of patients with TBI before (dotted lines) and after the rehabilitation (solid lines) in the CT1 test. (**a**) Patient 1. (**b**) Patient 2.

Figure 19 demonstrates that the positive dynamics could be seen of the rehabilitation of patients with TBI in the CT2 test.

**Figure 19.** Phase-connected pairs of EEG channels of a patients with TBI before (dotted lines) and after the rehabilitation (solid lines) in the CT2 test. (**a**) Patient 1. (**b**) Patient 2.

Figure 19 demonstrates that the positive dynamics could be seen of the rehabilitation of a patient with TBI in the CT2 test because interhemispheric connections or connections in the frontal cortex were activated in patients, as in control subjects. If interhemispheric connections or connections in the frontal cortex in the CT2 test appeared after rehabilitation, as in the control subjects, it the positive dynamics of rehabilitation could be concluded.

Let us consider an example of the lack of progress of the rehabilitation of a patient with TBI. Figure 20 demonstrates the dependence of *D* sorted in increasing order versus the numbers of pairs of EEG channels and its derivative for a patient with TBI.

**Figure 20.** The dependence of *D* sorted in increasing order (line 1) versus the numbers of a pair of EEG channels and its derivative (line 2) for a patient with TBI in the CT1 test. (**a**) before the rehabilitation; (**b**) after the rehabilitation.

Figure 20b shows that there was no sharp increase in the *D* derivative after the rehabilitation during the cognitive calculate-logical test, in contrast to Figure 20a. Thus, it was impossible to clearly determine the quantitative signs and identify phase-connected pairs by the suggested method. It could be concluded that there was no progress in rehabilitation during the cognitive calculation-logical test.

#### **5. Conclusions**

The paper presents an approach for segmenting long-term 19-channel EEG monitoring data. For the signals, the ridges of Morlet wavelet transform were calculated. Interchannel synchronization was used as a new feature of epileptic seizure. We also used the adaptive thresholding of the wavelet spectrogram ridges for signal segmentation. The intersection of the synchronized and the power spectral density intervals were obtained. As a result, the total duration of the fragments for analysis by the doctor was reduced by more than 60 times. It was shown that the frequency of the peak of the Fourier spectrum of the cutoff of the wavelet spectrogram at a frequency higher than the frequency of the ridge during an epileptic discharge was 2.5 times higher than the frequency of the Fourier peak corresponding to chewing. The Fourier peak full width at half maximum of the chewing artifact was 2 times larger than that of ES.

A comparison of the phases of EEG at the points of the Morlet wavelet spectrogram ridges were used for evaluation the EEG interchannel phase synchronization during cognitive tests in control subjects and patients with moderate TBI. Calculation-logical and spatial-pattern cognitive tests were used. Interhemispheric connections and connections in the frontal cortex in control subjects are initiated during the cognitive tests. The possibility of determining the positive dynamics of rehabilitation during the initialization of interhemispheric connections and connections in the frontal cortex of the brain or the absence of progress in rehabilitation has been shown.

**Author Contributions:** Conceptualization, Y.V.O., M.V.S. and L.A.Z.; methodology, Y.V.O.; software, I.A.K. and R.A.T.; validation, Y.V.O., M.V.S., L.A.Z., I.A.K. and R.A.T.; formal analysis, Y.V.O., I.A.K. and R.A.T.; investigation, Y.V.O., I.A.K. and R.A.T.; resources, M.V.S. and L.A.Z.; data acquisition, M.V.S., L.A.Z. and R.A.T.; data curation, M.V.S. and L.A.Z.; writing—original draft preparation, Y.V.O., I.A.K. and R.A.T.; writing—review and editing, Y.V.O., I.A.K. and R.A.T.; visualization, I.A.K. and R.A.T.; supervision, Y.V.O.; project administration, Y.V.O.; funding acquisition, Y.V.O. All authors have read and agreed to the published version of the manuscript.

**Funding:** The work was carried out within the framework of the state task.

**Institutional Review Board Statement:** The study was conducted according to the guidelines of the Declaration of Helsinki, and approved by the Ethics Committee of National Medical Research Center for Neurosurgery named after Academician N.N. Burdenko and by the Interuniversity Ethics Committee.

**Informed Consent Statement:** Informed consent was obtained from all subjects involved in the study.

**Data Availability Statement:** The data presented in this study are available on request from the corresponding author. The clinical data are not publicly available due to the ethical policy of the institute.

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analysis, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results.

#### **References**

## *Article* **Motion Artifacts Correction from Single-Channel EEG and fNIRS Signals Using Novel Wavelet Packet Decomposition in Combination with Canonical Correlation Analysis**

**Md Shafayet Hossain 1, Muhammad E. H. Chowdhury 2,\*, Mamun Bin Ibne Reaz 1,\*, Sawal Hamid Md Ali 1, Ahmad Ashrif A. Bakar 1, Serkan Kiranyaz 2, Amith Khandakar 2, Mohammed Alhatou 3, Rumana Habib <sup>4</sup> and Muhammad Maqsud Hossain <sup>5</sup>**

**Abstract:** The electroencephalogram (EEG) and functional near-infrared spectroscopy (fNIRS) signals, highly non-stationary in nature, greatly suffers from motion artifacts while recorded using wearable sensors. Since successful detection of various neurological and neuromuscular disorders is greatly dependent upon clean EEG and fNIRS signals, it is a matter of utmost importance to remove/reduce motion artifacts from EEG and fNIRS signals using reliable and robust methods. In this regard, this paper proposes two robust methods: (i) Wavelet packet decomposition (WPD) and (ii) WPD in combination with canonical correlation analysis (WPD-CCA), for motion artifact correction from single-channel EEG and fNIRS signals. The efficacy of these proposed techniques is tested using a benchmark dataset and the performance of the proposed methods is measured using two wellestablished performance matrices: (i) difference in the signal to noise ratio (Δ*SNR*) and (ii) percentage reduction in motion artifacts (*η*). The proposed WPD-based single-stage motion artifacts correction technique produces the highest average Δ*SNR* (29.44 dB) when db2 wavelet packet is incorporated whereas the greatest average *η* (53.48%) is obtained using db1 wavelet packet for all the available 23 EEG recordings. Our proposed two-stage motion artifacts correction technique, i.e., the WPD-CCA method utilizing db1 wavelet packet has shown the best denoising performance producing an average Δ*SNR* and *η* values of 30.76 dB and 59.51%, respectively, for all the EEG recordings. On the other hand, for the available 16 fNIRS recordings, the two-stage motion artifacts removal technique, i.e., WPD-CCA has produced the best average Δ*SNR* (16.55 dB, utilizing db1 wavelet packet) and largest average *η* (41.40%, using fk8 wavelet packet). The highest average Δ*SNR* and *η* using single-stage artifacts removal techniques (WPD) are found as 16.11 dB and 26.40%, respectively, for all the fNIRS signals using fk4 wavelet packet. In both EEG and fNIRS modalities, the percentage reduction in motion artifacts increases by 11.28% and 56.82%, respectively when two-stage WPD-CCA techniques are employed in comparison with the single-stage WPD method. In addition, the average Δ*SNR* also increases when WPD-CCA techniques are used instead of single-stage WPD for both EEG and fNIRS signals. The increment in both Δ*SNR* and *η* values is a clear indication that two-stage WPD-CCA performs relatively better compared to single-stage WPD. The results reported using the proposed methods outperform most of the existing state-of-the-art techniques.

**Citation:** Hossain, M.S.; Chowdhury, M.E.H.; Reaz, M.B.I.; Ali, S.H.M.; Bakar, A.A.A.; Kiranyaz, S.; Khandakar, A.; Alhatou, M.; Habib, R.; Hossain, M.M. Motion Artifacts Correction from Single-Channel EEG and fNIRS Signals Using Novel Wavelet Packet Decomposition in Combination with Canonical Correlation Analysis. *Sensors* **2022**, *22*, 3169. https://doi.org/10.3390/ s22093169

Academic Editor: Yvonne Tran

Received: 24 February 2022 Accepted: 25 March 2022 Published: 21 April 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

**Keywords:** motion artifact; electroencephalogram (EEG); functional near-infrared spectroscopy (fNIRS); wavelet packet decomposition (WPD); canonical correlation analysis (CCA)

#### **1. Introduction**

Due to the paradigm shift of hospital-based treatment in the direction of wearable and ubiquitous monitoring, nowadays, the acquisition and processing of vital physiological signals have become prevalent in the ambulatory setting. Since the acquisition of physiological signals is inclined to movement artifacts that happen due to the deliberate and/or voluntary movement of the patient during signal procurement utilizing wearable devices, restricting patients totally from physical movements, intentional and/or unintentional, is exceptionally troublesome. As a result, the physiological signals may get corrupted to some degree by motion artifacts. In some instances, this defilement may end up so conspicuous that the recorded signals may lose their usability unless the movement artifacts are diminished significantly.

Electroencephalogram (EEG) measures the electrical activity of the human brain quantitatively which took place due to the firing of neurons [1] and such brain activity is recorded utilizing a good number of cathodes which are located at different regions of the scalp [2]. EEG is one of the key diagnostic tests for epileptic seizure detection [3,4]. Other decisive utilization of EEG includes the estimation of drowsiness levels [5–8], emotion detection [9], cognitive workload [6,10], and brain-computer interfaces (BCIs) [11–16]. All of which have potential applications in the personal healthcare domain. Lately, the implementation of EEG-based biometric systems utilizing the inborn anti-spoofing capability of EEG signals was studied and appeared to be promising [17].

The functional near-infrared spectroscopy (fNIRS), a non-invasive optical brain imaging technique, measures changes in hemoglobin (Hb) concentrations inside the human brain [18] by employing light of various wavelengths in the infrared band and estimating the difference in the optical absorption [19]. Medical applications of fNIRS mainly focus on the noninvasive measurement of brain functions [20,21], cognitive tasks identification [22,23], and BCI [24–26].

Apart from movement artifacts, physiological signals undergo other types of artifacts as well. Gradient artifacts (GA) and pulse artifacts (PA) are the two most frequent artifacts observed in EEG during the simultaneous EEG-fMRI tests [27–29]. On the other hand, event-related fNIRS signals are regularly sullied by heartbeat, breath, Mayer waves, etc., as well as extra-cortical physiological clamors from the superficial layers [30].

Numerous attempts were made to reduce motion artifacts from EEG previously, which were summarized in [31,32]. In [33], the performance of motion artifacts correction techniques utilizing discrete wavelet transform (DWT) [34], empirical mode decomposition (EMD) [35], ensemble empirical mode decomposition (EEMD) [36], EMD along with canonical correlation analysis (EMD-CCA), EMD with independent component analysis (EMD-ICA), EEMD with ICA (EEMD-ICA), and EEMD with CCA (EEMD-CCA) were reported. Maddirala and Shaik [37] used singular spectrum analysis (SSA) [38], whereas DWT along with the thresholding technique was utilized in [39]. Gajbhiye et al. [40] employed wavelet-based transform along with the total variation (TV) and weighted TV (WTV) denoising techniques, whereas in [41], wavelet domain optimized Savitzky–Golay filter was proposed for the removal of motion artifacts from EEG. Recently, Hossain et al. [42] utilized variational mode decomposition (VMD) [43] for the correction of motion artifacts from EEG signals.

In the last few decades, multiple motion artifacts removal techniques were proposed [44–46] for the removal of motion artifacts from the fNIRS signal. Sweeney et al. [47] used adaptive filter, Kalman Filter, and EEMD-ICA. Scholkmann et al. [48] utilized the moving standard deviation and spline interpolation method, whereas in [49], a waveletbased method was proposed. The authors of [33] used DWT, EMD, EEMD, EMD-ICA, EEMD-ICA, EMD-CCA, and EEMD-CCA. In [50], Barker et al. used an autoregressive model-based algorithm, while kurtosis-based wavelet transform was proposed in [51], and Siddiquee et al. [52] utilized nine-degree of freedom inertia measurement unit (IMU) data to mathematically estimate the movement artifacts in the fNIRS signal using autoregressive exogenous (ARX) input model. A hybrid algorithm was proposed in [53] to filter out the movement artifacts from fNIRS signals where both the spline interpolation method and Savitzky–Golay filtering were employed. Very recently, the two-stage VMD-CCA technique was employed in [42].

The development of robust algorithms that can successfully reduce motion artifacts significantly from EEG and fNIRS data is critical; otherwise, the signals' interpretation could be erroneous by medical doctors and/or machine-learning-based applications. As mentioned earlier, DWT, EMD, EEMD, VMD, DWT-ICA, EMD-ICA, EEMD-ICA, EMD-CCA, EEMD-CCA, VMD-CCA, etc. were the most commonly used methods for the correction of motion artifacts from EEG and fNIRS signals. ICA and CCA cannot be used independently for single-channel EEG/fNIRS motion artifacts correction as the input of ICA/CCA algorithms require at least two (or more) channels data, whereas DWT, EMD, EEMD, VMD, etc. algorithms suffer from several limitations which are discussed in the discussion section of this paper. Additionally, there is still room for improvement for Δ*SNR* and *η* values which can be achieved using other effective novel methods. Therefore, in this paper, two novel motion artifacts removal techniques have been proposed which can eliminate motion artifacts from single-channel EEG and fNIRS signals to a great extent. The first is a single-stage motion artifacts correction technique using the wavelet packet decomposition (WPD), whereas the other novel method is WPD in combination with CCA (WPD-CCA), a two-stage motion artifacts removal technique, as the name suggests.

In this extensive study, for the correction of motion artifact from EEG and fNIRS signals using the WPD method, four different wavelet packet families (Daubechies (dbN), Symlets (symN), Coiflets (coifN), Fejer-Korovkin (fkN)) have been used with three different vanishing moments (for each of the wavelet packets) that resulted in a total of 12 different investigations. The wavelet packets used in the WPD method are db1, db2, db3, sym4, sym5, sym6, coif1, coif2, coif3, fk4, fk6, and fk8. To the best of our knowledge, the WPD algorithm has not been used for the removal of motion artifacts from single-channel EEG and fNIRS signals to date. WPD-CCA method is another novel contribution of this research work where Daubechies and Fejer-Korovkin wavelet packet families are utilized. In the WPD-CCA technique, db1, db2, db3, fk4, fk6, and fk8 have been used separately, resulting in six different investigations to reduce motion artifacts from EEG and fNIRS signals more efficiently.

The rest of this paper is organized as follows: Section 2 discusses the theoretical background of the different algorithms (WPD, CCA, WPD-CCA) investigated here, while Section 3 provides brief information about the EEG and fNIRS benchmark dataset and experimental methodology. Section 4 provides the results of the artifact removal techniques proposed in this work and Section 5 covers the discussion. Finally, the paper is concluded in Section 6.

#### **2. Theoretical Background**

#### *2.1. Wavelet Packet Decomposition (WPD)*

Using the WPD technique, signals can be decomposed into a wavelet packet basis at diverse scales [54,55]. For *j*-level decomposition, a wavelet packet basis is represented by multiple signals *<sup>n</sup>* − <sup>2</sup>*<sup>j</sup> k <sup>k</sup>*∈Z, where *<sup>i</sup>* <sup>∈</sup> <sup>Z</sup>+, 0 <sup>≤</sup> *<sup>i</sup>* <sup>≤</sup> <sup>2</sup>*<sup>j</sup>* <sup>−</sup> 1. The wavelet packet bases *ψi j* (*n*) are produced recursively from the scaling and wavelet functions, *ψ*<sup>0</sup> <sup>1</sup>(*n*) = *φ*(*n*) and *ψ*1 <sup>1</sup>(*n*) = *ψ*(*n*)**,** respectively, as follows:

$$
\psi\_j^{2i}(n) = \sum\_k h(k)\psi\_{j-1}^i \left(n - 2^{j-1}k\right) \tag{1}
$$

$$
\Psi\_j^{2i+1}(n) = \sum\_k g(k)\psi\_{j-1}^i \left(n - 2^{j-1}k\right) \tag{2}
$$

where *h*(*n*) represents lowpass filter and *g*(*n*) is the highpass filter defined as [54,56]:

$$h(k) = \left\langle \psi\_{j}^{2i}(u), \psi\_{j-1}^{i}\left(u - 2^{j-1}k\right) \right\rangle \tag{3}$$

$$\log(k) = \left\langle \Psi\_j^{2i+1}(u), \Psi\_{j-1}^i \left(u - 2^{j-1}k\right) \right\rangle \tag{4}$$

The decomposition of a signal *x*(*n*) onto the wavelet basis *j*(*n*) at level *j* can be expressed as: *<sup>x</sup>*(*n*) <sup>=</sup> <sup>∑</sup>*i*,*<sup>k</sup> <sup>X</sup><sup>i</sup>*

$$\kappa(n) = \sum\_{i,k} X\_j^i \psi\_j^i \left( n - 2^j k \right) \tag{5}$$

where *X<sup>i</sup> j* (*k*) signifies the *k*th wavelet coefficient of the packet *i*, at level *j*. Here, *X<sup>i</sup> j* (*k*) represents the intensity of the localized wavelet *ψ<sup>i</sup> j <sup>n</sup>* − <sup>2</sup>*<sup>j</sup> k* , defined by:

$$X\_j^i(k) = \left\langle \mathbf{x}(n), \,\Psi\_j^i\left(n - 2^j k\right) \right\rangle \tag{6}$$

Let *x*(*n*) represent a recorded EEG/fNIRS signal which can be expressed as the sum of a source signal *s*(*n*) and a motion artifact signal *v*(*n*) as follows:

$$\mathbf{x}(n) = \mathbf{s}(n) + \mathbf{v}(n) \tag{7}$$

In general, the source signal *s*(*n*) is assumed to be normally distributed having a mean value equals to zero, *<sup>s</sup>*(*n*) ∼ *<sup>N</sup>*(0, *<sup>σ</sup>*), where *<sup>σ</sup>*<sup>2</sup> characterizes the variance of *<sup>s</sup>*(*n*) [57]. On the other hand, general assumptions regarding the artifact signal *v*(*n*) includes temporal localization, not normally distributed with high local variance.

According to [58], *X<sup>i</sup> j* (*k*), can be represented as the sum of *S<sup>i</sup> j* (*k*) and *V<sup>i</sup> <sup>j</sup>* (*k*), where *Xi j* (*k*), *S<sup>i</sup> j* (*k*), and *V<sup>i</sup> <sup>j</sup>* (*k*) are the wavelet coefficients of *x*(*n*), *s*(*n*), and *v*(*n*), respectively:

$$X\_j^i(k) = S\_j^i(k) + V\_j^i(k) \tag{8}$$

It is noteworthy to mention that the wavelet coefficients *V<sup>i</sup> <sup>j</sup>* (*k*) will be sparse as well as the non-zero coefficients will have a relatively higher magnitude as the variance of *v*(*n*) is locally high, which would cause an increase in the local variance of the recorded EEG/fNIRS signal *x*(*n*).

#### *2.2. Canonical Correlation Analysis (CCA)*

CCA [59] is one of the most popular blind source separation methods which has the capability of dissociating multiple mixed or noisy signals. Assuming linear mixing, square mixing, and stationary mixing [60], the CCA technique computes an un-mixing matrix **W**,

which helps identify the unknown independent components **^ S** from a matrix **X**, which is a recorded multi-channel signal as follows:

$$
\hat{\mathbf{S}} = \mathbf{W}\mathbf{X} \tag{9}
$$

CCA also estimates the unknown independent components **^ S** using Equation (9) utilizing second-order statistics (SOS). CCA forcefully makes the sources to be autocorrelated maximally as well as makes the sources mutually uncorrelated [61]. Let us assume **y** as a linear combination of neighboring samples for an input signal **x** (i.e., *y*(*t*) = *x*(*t* − 1) + *x*(*t* + 1)) [62]. Consider the linear combinations of the components in **x** and **y**, known as the the canonical variates:

$$\mathbf{x} = \mathbf{w}\_{\mathbf{x}}^{T} \left(\mathbf{x} - \mathbf{\bar{x}}\right) \tag{10}$$

$$\mathbf{y} = \mathbf{w}\_{\mathbf{y}}^{T} \left(\mathbf{y} - \mathbf{\bar{y}}\right) \tag{11}$$

where **wx** and **wy** represents the weight matrices. CCA computes **wx** and **wy** in such a way so that the correlation ρ between *x* and *y* will be maximized [62]:

$$\rho = \frac{\mathbf{w}\_{\mathbf{x}}^T \mathbf{C}\_{xy} \mathbf{w}\_{\mathbf{y}}^T}{\sqrt{\mathbf{w}\_{\mathbf{x}}^T \mathbf{C}\_{xx} \mathbf{w}\_{\mathbf{x}} \mathbf{w}\_{\mathbf{y}}^T \mathbf{C}\_{yy} \mathbf{w}\_{\mathbf{y}}}} \tag{12}$$

where **C***xx* and **C***yy* signify the nonsingular within-set covariance matrices and **C**xy represent the between-sets covariance matrix. The maximized ρ is calculated by setting the derivatives of Equation (12) (with respect to **wx** and **wy)** equal to zero:

$$\begin{aligned} \mathbf{C}\_{xx}^{-1} \mathbf{C}\_{xy} \mathbf{C}\_{yy}^{-1} \mathbf{C}\_{yx}^{T} \hat{\mathbf{w}}\_{x} &= \rho^{2} \hat{\mathbf{w}}\_{\mathbf{x}} \\ \mathbf{C}\_{yy}^{-1} \mathbf{C}\_{yx} \mathbf{C}\_{xx}^{-1} \mathbf{C}\_{xy}^{T} \hat{\mathbf{w}}\_{\mathbf{y}} &= \rho^{2} \hat{\mathbf{w}}\_{\mathbf{y}} \end{aligned} \tag{13}$$

**wx** and **wy** can then be found out as the eigenvectors of the matrices **C**−<sup>1</sup> *xx* **C***xy***C**−<sup>1</sup> *yy* **C***<sup>T</sup> yx* and **C**−<sup>1</sup> *yy* **C***yx***C**−<sup>1</sup> *xx* **C***<sup>T</sup> xy*, respectively, and the corresponding eigenvalues ρ<sup>2</sup> are the squared canonical correlations. It is sufficient to solve only one of the eigenvalue equations to obtain the un-mixing matrix **W** as the solutions are related. Furthermore, the underlying source signals **^ S** can be estimated.

The components that seem to be artifacts can then be discarded by simply setting the corresponding columns of the **^ S** matrix to zero before the signal reconstruction.

#### *2.3. WPD-CCA*

The WPD algorithm can be utilized to decompose a single-channel signal into multichannel signal **X** where each column of matrix **X** represents the detailed and approximated sub-band signals. The total number of generated sub-band signals would be equal to 2*<sup>j</sup>* ,

where *<sup>j</sup>* denotes the level, a priori. To estimate the underlying true sources **^ S** (Equation (9)), these generated sub-band signals can then be used as the multi-channel input signals to

the CCA algorithm. After that, the component/s of **^ S** which seem to be artifacts can be

discarded by making the corresponding columns of the matrix **^ S** equal to zero. Bypassing this newly obtained source matrix through the inverse of the un-mixing matrix **W**<sup>−</sup>1, the multi-channel signals **^ <sup>X</sup>** can be obtained. Finally, the cleaner signal **^ x** can be produced by simply summing all the columns of the matrix **^ X**.

#### **3. Methods**

This section describes the benchmark dataset used, pre-processing, study design, motion component identification, and evaluation metrics.

#### *3.1. Dataset Description*

A publicly available PhysioNet dataset [32,33,63] is used in this study that contains "reference ground truth" and motion corrupted signals for both EEG and fNIRS modalities. The details of the data recording procedure for EEG and fNIRS modalities were mentioned in [47]. During the data acquisition, two channels having the same hardware properties were placed on the test subject's scalp at very close proximity (20 mm for EEG modality and 30 mm for fNIRS modality), where the first channel was impacted with motion artifacts for 10–25 s at regular 2 min interval and the second channel was left untouched and undisturbed for the entire recording period. From the unimpacted channel (2nd channel), the EEG/fNIRS signal was extracted, which was free from motion artifacts and referred to as "reference ground truth" signal, whereas the impacted channel (1st channel) provided EEG/fNIRS signal corrupted with motion artifacts. It is worthwhile to mention that both the motion corrupted and "reference ground truth" signals were extracted simultaneously from channels 1 and 2, respectively, for approximately 9 min for each of the trial/test subjects. Additionally, the same channels were used to extract EEG/fNIRS data from all of the test subjects.

Twenty-three sets of EEG recordings, sampled at 2048 Hz, collected from six patients in four different sessions, are available in the database. Each recording consists of one motion corrupted EEG signal and one reference "ground truth" EEG signal. The average correlation coefficient between the reference "ground truth" and motion corrupted EEG signals is very high over the epochs where the motion artifacts are absent and the average correlation coefficient drops significantly during the epochs of motion artifacts [32]. The superimposed reference "ground truth" and motion corrupted EEG signals are illustrated in Figure 1a.

**Figure 1.** Example of motion-corrupted EEG (**a**) and fNIRS (**b**) signals. Two signals (blue: ground truth and red: motion-corrupted) are highly correlated during the motion artifacts free epochs. Boxed areas show the epochs of motion corrupted signals. A zoomed version is presented underneath each sub-plot.

fNIRS signals were recorded at two different wavelengths: 690 nm and 830 nm wavelengths. There were 16 sets of fNIRS recordings (9 recordings at 830 nm wavelength and 7 recordings at 690 nm wavelength) in total from 10 test subjects at a sampling frequency of 25 Hz [33,63]. Like EEG recordings, each recording of fNIRS consists of one motion corrupted fNIRS signal and one "reference ground truth" fNIRS signal. The overlaid "reference ground truth" fNIRS signal and motion artifact contaminated fNIRS signal is depicted in Figure 1b.

#### *3.2. Signal Preprocessing*

Downsampling: As EEG signals can be partitioned into a few sub-bands, specifically delta (1–4 Hz), theta (4–8 Hz), alpha (8–13 Hz), beta (13–30 Hz), and gamma (30–80 Hz) [64], we downsampled all the 23 sets of EEG recordings from 2048 Hz to 256 Hz, which guarantees data reliability without losing any vital signal information and morphology. The fNIRS signals were not upsampled/downsampled as the original sampling rate was 25 Hz during acquisition.

Power line noise removal: To remove power line noise, a third-order Butterworth notch filter with a center frequency of 50 Hz was utilized to remove 50 Hz and its subsequent harmonics as a pre-processing technique for all the EEG and fNIRS signals.

Baseline Drift Correction: Both the EEG and fNIRS signals were found to have significant baseline drift, which is defined as undesired amplitude shifts in the signal that would result in inaccurate results if not corrected. To remove baseline drift from EEG and fNIRS recordings, a polynomial curve fitting method was used to estimate the baseline, which was then subtracted from the recorded raw signal.

#### *3.3. Study Design*

The simulations of this work were carried out in a PC with Intel(R) Core(TM) i5-8250U CPU at 1.80 GHz which was equipped with 8 GB RAM. In-house-built MATLAB code was written to pre-process the EEG and fNIRS data. The single-stage WPD and two stages WPD-CCA methods were deployed in "MATLAB R2020a, The MathWorks, Inc., Natick, Massachusetts, USA". Figure 2 depicts the motion artifacts elimination framework presented in this study. An automated way for identifying motion corrupted components of the preprocessed signal is also discussed.

**Figure 2.** Methodological framework for the motion artifact correction.