A Mathematical Review on EEG Channel Selection Techniques for Motor Imagery Classification

1. Introduction

Brain Computer Interface(BCI) systems measure the activities of the brain in order to perform specific tasks such as interpreting one’s will to speak or move. While there are numerous methods to measure the brain’s activity, such as EEG (electroencephalography), fMRI (functional magnetic resonance imaging), and MEG (Magnetoencephalography), BCI systems utilizing scalp EEG(henceforth called EEG) is trending as it has high temporal resolution, is non-invasive, and extremely portable.

Classification of EEG signals follows the procedure shown in Figure 1. [1,2]

EEG signals are measured by attaching electrodes on the scalp according to the montage, which is a standardized arrangement of electrodes and selection of channel pairs. Montages are divided into two groups: referential and bipolar, by how channels are formed. In referential montage, the voltage difference between each electrode and one reference electrode is recorded, and forms a channel. In Bipolar montage, the voltage difference among consecutive two electrodes are recorded, and forms a channel. The recorded signals then go through artifact and noise removal process to remove undesired signals.[3]

Then, channel selection reduces the number of channels used in feature extraction and classification and ensures that only meaningful channels are included in the process afterwards. After the channel selection, only the necessary channels’ data are needed to be acquired from subjects, so it can help reduce the time consuming setup process and inconvenience for subjects. Also, reducing the number of channels helps reduce computational complexity and increases classification accuracy. [2] Channel selection usually follows the following Figure 2:

Motor Imagery EEG (MI-EEG) is a measurement of EEG signals while subject simulates its movement without actually moving. Therefore MI-EEG based BCIs have been suggested as a substitute method for control, especially for people who are disabled or suffering from locked-in syndrome.[4] It is widely accepted that imagination of a movement stimulates same neural pathway as real movement[5], and can even lead to muscle growth[6]. Traditional Statistical filters used for MI-EEG analysis can be divided into CSP (Common Spatial pattern)-based and non CSP-based filters.[2] Therefore, this paper will discuss about CSP based, filters and non-CSP based statistical filters used in MI-EEG classification in the following chapters.

2. Common Spatial Pattern(CSP) Based Filters

2.1. Original CSP Filter

The Common Spatial Pattern algorithm is known as a spatial filter effective for EEG filtering by maximizing one class’s variance while minimizing other class’s variance in n orthogonal bases as Figure 3, and can be mathematically induced as follows Derivation.

Derivation^[7,8,9]

Let $\mathbf{X}\in {\mathbb{R}}^{C\times N\times T}$and $\mathbf{X}\left(t\right)\in {\mathbb{R}}^{C\times N}$ where C is the number of channels, N is the number of samples per epoch, and T is the number of epochs. $\mathbf{X}\left(t\right)$ represents the t th epoch’s signal. In this paper, $\mathbf{X}$ only consists of two classes: $c{l}_1$ (class 1) and $c{l}_2$ (class 2) for simplicity. Then normalized covariance matrix of each classes: ${\mathbf{C}}_{\mathbf{1}}$ and ${\mathbf{C}}_{\mathbf{2}}$ is defined as

$${\mathbf{C}}_{\mathbf{i}}=\sum _{j\in c{l}_i}{\displaystyle \frac{\mathbf{X}\left(j\right)\mathbf{X}{\left(j\right)}^T}{tr\left(\mathbf{X}\left(j\right)\mathbf{X}{\left(j\right)}^T\right)}}$$

(1)

The covariance matrix is normalized with the sum of the diagonal elements of $\mathbf{X}\left(j\right)\mathbf{X}{\left(j\right)}^T$ to remove variances among trials.

• As $\mathbf{X}\left(j\right)\mathbf{X}{\left(j\right)}^T$ is symmetric, ${\mathbf{C}}_{\mathbf{i}}$ is symmetric and therefore if we define $\mathbf{C}\triangleq ({\mathbf{C}}_{\mathbf{1}}+{\mathbf{C}}_{\mathbf{2}})$, it can be orthogonally diagonalized as

$$\mathbf{C}=({\mathbf{C}}_{\mathbf{1}}+{\mathbf{C}}_{\mathbf{2}})=\mathbf{U}\mathrm{\Lambda }{\mathbf{U}}^T$$

(2)

where $\mathrm{\Lambda }$ is a diagonal matrix consisting eigenvalues and $\mathbf{U}$ is an orthonormal matrix whose columns are eigenvectors of $\mathbf{C}$ and ${\mathbf{U}}^T\mathbf{U}=I$ as ${\mathbf{U}}^T={\mathbf{U}}^{-1}$. Therefore,

$${\mathbf{U}}^T({\mathbf{C}}_{\mathbf{1}}+{\mathbf{C}}_{\mathbf{2}})\mathbf{U}={\mathbf{U}}^T\mathbf{C}\mathbf{U}=\mathrm{\Lambda }$$

(3)

By applying whitening transform and simultaneously diagonalizing using method proposed in [9], we get

$${\overline{\mathbf{U}}}^T({\mathbf{C}}_{\mathbf{1}}+{\mathbf{C}}_{\mathbf{2}})\overline{\mathbf{U}}=\mathbf{I}$$

(4)

where $\overline{\mathbf{U}}=\mathbf{U}{\mathrm{\Lambda }}^{-{\displaystyle \frac{1}{2}}}$, and letting ${\mathbf{S}}_{\mathbf{1}}\triangleq {\overline{\mathbf{U}}}^T{\mathbf{C}}_{\mathbf{1}}\overline{\mathbf{U}}$ and ${\mathbf{S}}_{\mathbf{2}}\triangleq {\overline{\mathbf{U}}}^T{\mathbf{C}}_{\mathbf{2}}\overline{\mathbf{U}}$, ${\mathbf{S}}_{\mathbf{1}}$ and ${\mathbf{S}}_{\mathbf{2}}$ can be simultaneously diagonalized as

$${\mathbf{S}}_{\mathbf{1}}=\mathbf{V}{\mathrm{\Lambda }}_1{\mathbf{V}}^T,{\mathbf{S}}_{\mathbf{2}}=\mathbf{V}{\mathrm{\Lambda }}_2{\mathbf{V}}^T$$

(5)

sharing same eigenvectors satisfying

$${\mathrm{\Lambda }}_1+{\mathrm{\Lambda }}_2=\mathbf{I}$$

(6)

by properties of simultaneous diagonalization.

• Equation (5), (6) are important properties as it shows that when spanned by $\mathbf{V}$, the bigger the variance of one class is in one axis, the smaller the variance of other class in corresponding axis.
• By transforming ${\mathrm{\Lambda }}_1$ and ${\mathrm{\Lambda }}_2$ back to equation about ${\mathbf{C}}_{\mathbf{1}}$ and ${\mathbf{C}}_{\mathbf{2}}$, we get

$${\mathbf{Q}}^T{\mathbf{C}}_{\mathbf{1}}\mathbf{Q}+{\mathbf{Q}}^T{\mathbf{C}}_{\mathbf{2}}\mathbf{Q}=\mathbf{I}$$

(7)

where $\mathbf{Q}=\mathbf{U}{\mathrm{\Lambda }}^{-{\displaystyle \frac{1}{2}}}\mathbf{V}$. For an epoch $\mathbf{X}\in {\mathbb{R}}^{C\times N}$, we can calculate its CSP transformation as

$${\mathbf{X}}_{\mathbf{CSP}}={\displaystyle \frac{diag\left({\overline{Q}}^T\mathbf{X}{\mathbf{X}}^T\overline{Q}\right)}{sum\left(diag\left({\overline{Q}}^T\mathbf{X}{\mathbf{X}}^T\overline{Q}\right)\right)}}$$

(8)

where $\overline{Q}$ is a $2k\times C$ sized matrix which is composed of k best columns of Q that makes class 1’s variance large, and k best columns that makes class 2’s variance large. □

2.2. Filter Bank CSP^[10]

Typical frequency band analyzed in MI-EEG is 8-30Hz[11,12]. However, as EEG signals have variance among trials and subjects, algorithms’ performance can be enhanced by selecting specific frequency bands per subjects[13]. Therefore, the authors suggested to divide the frequency bands in small bands(4Hz in this paper), and apply CSP filters to each bands. This technique allowed automatic selection of meaningful frequency bands per subjects.

Although mutual information is not used very common in MI-EEG classification, MIBIF(Mutual Information-based Best Individual Feature) and MIRSR(Mutual Information-based Rough Set Reduction) was used to select features to be trained in this paper, which calculates mutual information among features and selects k features that best represent the data.

2.3. Sparse Common Spatial Pattern (SCSP)

CSP coefficients calculated above, Q, is dense, meaning that most of the entries in the matrix can be considered meaningful, and only some have negligible effect. However, the matrix $\overline{Q}$ that includes the selected channels’ coefficients select only some of the columns of Q; therefore the original signal can not be transformed into a shape that can be discriminated clearest. Therefore, the authors have sparsified the CSP spatial filters, which can increase performance as the removed features include relatively unnecessary information compared to original CSP filters.

As CSP algorithm’s projection matrix $\mathbf{W}$ can be written as

$$\begin{array}{cc}\hfill \underset{{\mathbf{W}}_i}{min}\sum _{i=1}^{i=m}{\mathbf{W}}_i{\mathbf{C}}_{\mathbf{1}}{{\mathbf{W}}_i}^T+& \sum _{i=m+1}^{i=2m}{\mathbf{W}}_i{\mathbf{C}}_{\mathbf{2}}{{\mathbf{W}}_i}^T\hfill \\ \hfill \mathrm{while}{\mathbf{W}}_i({\mathbf{C}}_{\mathbf{1}}+{\mathbf{C}}_{\mathbf{2}}){\mathbf{W}}_i^T=1,& i=\{1,2,\cdots 2m\}\hfill \\ \hfill \mathrm{and}{\mathbf{W}}_i({\mathbf{C}}_{\mathbf{1}}+{\mathbf{C}}_{\mathbf{2}}){\mathbf{W}}_j^T=0,& i,j=\{1,2,\cdots 2m\}i\ne j\hfill \end{array}$$

(9)

Sparsity of W can be increased with several methods, including $L_0$ norm, $L_p$ norm, and $L_1/L_2$ norm, where $L_p$ norm of $\mathbf{X}$ is defined as

$${\left|\right|\mathbf{X}\left|\right|}_p={\left(\sum _{i=1}^n{\left|\mathbf{X}\right|}^p\right)}^{1/p}$$

(10)

where $|·|$ represents magnitude of vector $(·)$, making $L_0$ norm the number of non-zero entries of $\mathbf{X}$ and $L_1$ norm the sum of magnitudes of vectors, also known as Manhattan distance.

According to [15], the $L_1/L_2$ norm is more suit as a measure of non-sparsity than $L_1$ and $L_p$ norm as it is bounded within $(0,\sqrt{n})$, while $L_1$ and $L_p$ norm does not have an upper bound. Also, $L_1/L_2$ norm does not depend on vectors’ size, while $L_1$ and $L_p$ depends on it.

Therefore, $L_1/L_2$ norm is adopted to (9) with the scaling factor r and $(1-r)$ multiplied to it, forming

$$\begin{array}{cc}\hfill \underset{{\mathbf{W}}_i}{min}(1-r)\sum _{i=1}^{i=m}{\mathbf{W}}_i{\mathbf{C}}_{\mathbf{1}}{{\mathbf{W}}_i}^T+& \sum _{i=m+1}^{i=2m}{\mathbf{W}}_i{\mathbf{C}}_{\mathbf{2}}{{\mathbf{W}}_i}^T+r\sum _{i=1}^{2m}{\displaystyle \frac{\left|\right|{\mathbf{w}}_{\mathbf{i}}{\left|\right|}_1}{\left|\right|{\mathbf{w}}_{\mathbf{i}}{\left|\right|}_2}}\hfill \\ \hfill \mathrm{while}{\mathbf{W}}_i({\mathbf{C}}_{\mathbf{1}}+{\mathbf{C}}_{\mathbf{2}}){\mathbf{W}}_i^T=1,& i=\{1,2,\cdots 2m\}\hfill \\ \hfill \mathrm{and}{\mathbf{W}}_i({\mathbf{C}}_{\mathbf{1}}+{\mathbf{C}}_{\mathbf{2}}){\mathbf{W}}_j^T=0,& i,j=\{1,2,\cdots 2m\}i\ne j\hfill \end{array}$$

(11)

The original CSP algorithm is a complementary case when $r=0$.

2.4. Robust Sparsd Common Spatial Pattern(RSCSP)

The RSCSP algorithm is a derivative algorithm of the SCSP algorithm developed to reduce the effect of outliers in the EEG signal. The non-robust SCSP algorithm calculates covariance matrix of class i: $C_i$ as (1). As EEG signals are regularized, $C_i$ can be written in the form of

$${\mathbf{C}}_w={\displaystyle \frac{1}{(t\times n_w)-1}}{\mathbf{E}}_w{\mathbf{E}}_w^{\mathbf{T}}$$

(12)

, where $w=1,2$ and ${\mathbf{E}}_w\in {\mathbb{R}}^{C\times (t·n_w)}$ denotes EEG signals of class w concatenated by channels. As outliers can distort covariance matrices, MCD(Minimum Covariance Determinant) is used. MCD can be described as a process to find a subset H of original matrix $X=[c_1,c_2,\cdots c_n]$, such that its determinant $det\left({\Sigma }_H\right)$ is reduced, and can be mathematically represented as an optimization problem of ${\displaystyle \underset{H\subset X,\left|H\right|=h}{min}det\left({\Sigma }_H\right)}$ where

$${\mathbf{C}}_w={\displaystyle \frac{1}{(\alpha \times t\times n_w)-1}}{\mathbf{E}}_w{\mathbf{E}}_w^{\mathbf{T}}$$

(13)

when $1-\alpha$ is the breakdown point. FASTMCD algorithm can also be adopted to decrease the complexity of the calculation. FASTMCD algorithm first selects random h points and updates the points by ranking the data by their Mahalanobis distances and select h points with the smallest distances to form a new subset. Then if the new subset has a smaller determinant of covariance matrix, it functions as a new initial set and continues this process until the covariance matrix converges.

3. CSP Based Channel Selection Techniques

3.1. Bhattacharyya Bound of CSP Features

In [7], the authors reduced the number of channels by selecting a set of channels using forward search. Bhattacharyya bound, which is the upper bound of bayes error probability was used as evaluation, as the calculation of Bayes error probability is challenging. The statement: Bhattacharyya bound is the upper bound of Bayes error probability can be proved by Proof.

To summarize, the channels are selected by Algorithm 1

• Proof.

The Bayes error probability $P_e$ can be written as

$$P_e\triangleq {\int }_{\mathbf{X}}min(P\left(C_1\right)p\left(x\right|C_1),P\left(C_2\right)p\left(x\right|C_2))dx$$

(14)

where $P\left(C_1\right)$ and $P\left(C_2\right)$ are the prior probabilities of class 1 and 2, respectively. $p\left(x\right|C_1)$ and $p\left(x\right|C_2)$ are are class conditional probability density functions, and x represents an observed data point.

• Sinces $min(a,b)\le \sqrt{ab}$,

$$min(P\left(C_1\right)p\left(x\right|C_1),P\left(C_2\right)p\left(x\right|C_2))\le \sqrt{P\left(C_1\right)P\left(C_2\right)p\left(x\right|C_1)p\left(x\right|C_2)}$$

(15)

and integrating both sides over $\mathbf{X}$ results

$$\begin{array}{cc}\hfill P_e& \le {\int }_{\mathbf{X}}\sqrt{P\left(C_1\right)P\left(C_2\right)p\left(x\right|C_1)p\left(x\right|C_2)}dx\hfill \end{array}$$

(16)

$$\begin{array}{cc}\hfill & =\sqrt{P\left(C_1\right)P\left(C_2\right)}{\int }_{\mathbf{X}}\sqrt{p\left(x\right|C_1)p\left(x\right|C_2)}dx\hfill \end{array}$$

(17)

By [17], (17) is a Bhattacharyya bound, which is a complementary case of Chernoff bound where $\beta =0.5$, yielding

$$P_e\le \sqrt{P\left(C_1\right)P\left(C_2\right)}{e}^{-k(1/2)}$$

(18)

where

$$k(1/2)={\displaystyle \frac{1}{8}}{({\mu }_2-{\mu }_1)}^T{\left[{\displaystyle \frac{{\Sigma }_1+{\Sigma }_2}{2}}\right]}^{-1}({\mu }_2-{\mu }_1)+{\displaystyle \frac{1}{2}}ln{\displaystyle \frac{|\frac{{\Sigma }_1+{\Sigma }_2}{2}|}{\sqrt{|{\Sigma }_1\left|\right|{\Sigma }_2|}}}$$

(19)

and ${\mu }_i$ and ${\Sigma }_i$ represent mean vectors and covariance matrix of class i’s CSP feature vectors, respectively.

Algorithm 1 Bhattacharyya Bound Based Channel Selection

3.2. CSP-Rank

CSP-Rank is based on the idea that the matrix $\overline{Q}$ in (8) includes the importance of each channel in classification as it acts as a weight to the original signal. The authors of this paper therefore set $k=1$ and selected two vectors: $S{F}_1$ and $S{F}_2$, where $S{F}_i$ maximizes class i’s variance. Therefore, the channels with the largest entries in each of $S{F}_1$, $S{F}_2$, $S{F}_3$, ⋯ are consecutively selected until the model trained with the selected channels reaches the desired accuracy.

The results are statistically significant at $p<0.05$ significance compared to random selection. Also, the number of channels needed to reach 90% accuracy was 8-38 channels, while SVM-RFE(Support-vector machine recursive feature elimination, explained in Section 4.1) required 12-28 channels, therefore being better than SVM-RFE.

3.3. $L_1$ Norm of CSP Filter

Since it is widely accepted that the effectiveness of the CSP algorithm decreases as the number of channels increases due to overfitting, the authors first applied the CSP filter to all channels and scored each channel based on its $L_1$ norm. Then, the channels with high scores were selected and used to train the model. The following score is used for calculating channels’ weight:

$$SC\left(i\right)={\displaystyle \frac{\left|\right|{\overline{Q}}_i{\left|\right|}_1}{\left|\right|\overline{Q}{\left|\right|}_1}}$$

(20)

where $\overline{Q}\in {\mathbb{R}}^{2k\times C}$ is the CSP projection matrix from (8) , and ${\overline{Q}}_i$ is the ith column of $\overline{Q}$, and can be ordered in decreasing order to select wanted number of channels.

3.4. Cohen’s d Effect Size CSP(E-CSP)

The E-CSP algorithm aims to remove channels that contain redundant information. The algorithm is divided into two major parts: part 1, where the noisy channels are removed, and part 2, where the remaining channels are evaluated to determine if they are meaningful.

First, to remove trials with large artifacts, the z-score calculated by equation, z-score, calculated by (21) is used.

$$Z_{li}^j={\displaystyle \frac{|C_{li}^j-\overline{C_{li}}|}{\sigma }}$$

(21)

where i, j, l represents channel number, trial,and class, respectively, and $\sigma$ is the standard deviation of lth task’s ith channel, is calculated. The z-score was adopted as it measures how far the signal deviates from the average. Therefore, the ith channel’s jth trial can be considered noisy if $z_{li}^j>\alpha$, where $\alpha$ is the criterion of noisy channel. Therefore, let $s_{ij}^l=\left\{\begin{array}{cc}1,\hfill & (z_{li}^j>\alpha )\hfill \\ 0,\hfill & \left(else\right)\hfill \end{array}\right.$ and let the frequency of lth class’ jth trial being noisy as

$${\displaystyle f_l^j=\frac{{\sum }_{i=1}^N{s}_{ij}^l}{N}\times 100}$$

(22)

If the trial’s frequency exceeds the threshold $\beta \%$, it is considered as a noisy trial. Therefore, we can define the set of noisy trials as

$$K_l=\left\{j\right|f_l^j>\beta \};\forall j$$

(23)

Therefore, after removing the noisy channels, the ith channel’s mean for jth task can be calculated as

$$\overline{c_{li}}={\displaystyle \frac{1}{|J_s^l|}}\sum _{j\in J_s^l}{c}_{li}^j$$

(24)

where $J_s^l$ represents non-noisy channels.

Second, the effect size $d_i$ of channel i is evaluated as a meaningful channel if

$$d_i={\displaystyle \frac{|\overline{c_{1i}}-\overline{c_{2i}}|}{\sigma }}>\gamma$$

(25)

where

$$\sigma ={\displaystyle \frac{{\sigma }_{1i}^j+{\sigma }_{2i}^j}{2}};\forall j\in J_s^l$$

(26)

and ${\sigma }_{li}^j$ stands for the standard deviation of lth class’ ith channel, and $\gamma$ represents the criterion for a channel to be thought as meaningless, as the variance among two class is small.

Finally, the leftover channels L, are selected and CSP filter is applied.

4. Non-CSP Based Methods

4.1. Support-Vector Machine Recursive Feature Elimination (SVM-RFE)^[21,22]

SVM(Hard Margin Support Vector Machine) Algorithm

Support vector machine is an algorithm that aims to find a hyperplane that maximizes the margin among the given classes. Let there be N data points: $(x_1,y_1),(x_2,y_2),\cdots (x_N,y_N)$ where d is the dimension of the point, $x_i\in {\mathbb{R}}^d$ and $y_i\in \{1,2\}$ be a label. The hyperplane can be written as $follows\left(x\right)=w^{T}x+b$ where w is the weight and b is the bias. The optimization problem of w and b to maximize the margin can be written as

$${\displaystyle J=\underset{w,b}{min}\frac{1}{2}{\left|\right|w\left|\right|}^2\mathrm{subject}\mathrm{to}{y}_i(w^T{x}_i+b)\ge 1,\forall i}$$

(27)

, and this problem can be solved using the Lagrange multiplier.

• Let Lagrange function

$$L(w,b,\alpha )={\displaystyle \frac{1}{2}}{\left|\right|w\left|\right|}^2-\sum _{i=1}^N{\alpha }_i(y_i(w^T{x}_i+b)-1)$$

(28)

, where $\alpha \ge 0$ is the Lagrange multiplier. By partial differentiating with w, we get

$${\displaystyle \frac{\partial L}{\partial w}}=w-\sum _{i=1}^N{\alpha }_i{x}_i{y}_i=0⟹w=\sum _{i=1}^N{\alpha }_i{x}_i{y}_i$$

(29)

By partial differentiating with b,

$${\displaystyle \frac{\partial L}{\partial b}}=-\sum _{i=1}^N{\alpha }_i{y}_i=0⟹\sum _{i=1}^N{\alpha }_i{y}_i=0$$

(30)

Therefore, we get a dual Lagrangian function

$$L_D\left(\alpha \right)=\sum _{i=1}^N{\alpha }_i-{\displaystyle \frac{1}{2}}\sum _{i=1}^N\sum _{j=1}^N{\alpha }_i{\alpha }_j{y}_i{y}_j{x}_i^T{x}_j,$$

(31)

making (27)

$$max\sum _{i=1}^N{\alpha }_i-{\displaystyle \frac{1}{2}}\sum _{i=1}^N\sum _{j=1}^N{\alpha }_i{\alpha }_j{y}_i{y}_j(x_i·x_j)\mathrm{subject}\mathrm{to}\sum _{i=1}^N{\alpha }_i{y}_i=0$$

(32)

Soft Margin SVM Algorithm

However, as real datasets are not discriminated using linearly, slack variable ${\xi }_i$ is adopted, changing the objective (27) into

$${\displaystyle J=\underset{w,b}{min}\frac{1}{2}{\left|\right|w\left|\right|}^2+C\sum _{i=1}^N{\xi }_i\mathrm{subject}\mathrm{to}{y}_i(w^T{x}_i+b)\ge 1-{\xi }_i,{\xi }_i\ge 0,\forall i}$$

(33)

and can be derived into a same shape as (32) using the same method from above as ${\displaystyle C\sum _{i=1}^N{\xi }_i}$ is a constant term.

SVM-RFE Algorithm

SVM-RFE algorithm aims to find a subset of size r among n variables $(r<n)$ which maximizes the performance of the classifier using backward sequential selection.

From previous research [22,23], it is well known that $\mathrm{\Delta }J$ can be used as a criterion for elimination. The equation (33) can be transformed into

$$\mathrm{\Delta }J\left(i\right)={\displaystyle \frac{\partial J}{\partial w_i}}\mathrm{\Delta }{w}_i+{\displaystyle \frac{{\partial }^{2}J}{\partial w_i^2}}{(\mathrm{\Delta }{w}_i)}^2$$

(34)

using Taylor series to the second order. As the objective function J reaches its optimal point, ${\displaystyle \frac{\partial J}{\partial w_i}}\mathrm{\Delta }{w}_i=0$. Therefore,

$$\mathrm{\Delta }J\left(i\right)={\displaystyle \frac{{\partial }^{2}J}{\partial w_i^2}}{(\mathrm{\Delta }{w}_i)}^2$$

(35)

. When removing ith feature, the weight $w_i$ must be 0, therefore $\mathrm{\Delta }{w}_i=w_i$. Therefore, the criterion for removing the ith feature can be approximated as

$$C_i=\mathrm{\Delta }J\left(i\right)={\left(w_i\right)}^2$$

(36)

4.2. ${\gamma }^2$ Value Based Channel Selection

Similar to how the Bhattacharyya bound-based CSP selection used the F-ratio for selecting the initial subset of channels, the authors of [19] have utilized the ${\gamma }^2$ value as a control group. ${\gamma }^2$ is defined as

$${\gamma }^2\left(i\right)={\left({\displaystyle \frac{\sqrt{n_1{n}_2}}{n_1+n_2}}{\displaystyle \frac{\overline{\left|\right|x_i^{\left(1\right)}{\left|\right|}_2}-\overline{\left|\right|x_i^{\left(2\right)}{\left|\right|}_2}}{stdev\left(\right|\left|x_i^{(1,2)}\right|{|}_2)}}\right)}^2$$

(37)

where $k\in \left\{1(\mathrm{alias}\mathrm{of}c{l}_1\right),2(\mathrm{alias}\mathrm{of}c{l}_2\left)\right\}$, $x_i^{\left(k\right)}\in {\mathbb{R}}^T$ denotes ith channel of $X^{\left(k\right)}$, which means $X^{\left(k\right)}={[x_1^{\left(k\right)},x_2^{\left(k\right)},\cdots ,x_N^{\left(k\right)}]}^{\mathbf{T}}$, and $n_k$ denotes number of epochs in class k.

4.3. Correlation Based Channel Selection^[24,25]

EEG has a limitation that even though one is repeating the same task, due to other activities and thoughts in the brain, different EEG is recorded, and makes it hard to analyze the task performed. The authors assumed that as channels related to MI-EEG while performing same task will contain similar information across trials, thus having high correlation among trials, while channels that are not related will not have alike information, therefore having low correlation.

Pearson’s correlation is used in both papers, which measures linear dependency among random variables, and is calculated by

$$\rho (X,Y)={\displaystyle \frac{1}{N-1}}\sum _{i=1}^n\left({\displaystyle \frac{X_i-\overline{X}}{{\sigma }_X}}\right)\left({\displaystyle \frac{Y_i-\overline{Y}}{{\sigma }_Y}}\right)$$

(38)

where $X,Y$ are random variables, $\overline{X},\overline{Y}$ are means of random variables, and ${\sigma }_X,{\sigma }_Y$ represents standard deviation of X and Y respectively. Correlations are calculated among every channels, forming a correlation matrix of size $C\times C$. By averaging each row, we can find the channel that is most effective. By repeating this procedure among every epochs, we can select $N_s$ number of channels by selecting the channels that are considered to be effective over certain criteria times.

4.4. Genetic Algorithm (GA) Based Selection

Genetic Algorithm(henceforth GA) is a heuristic methodology for finding an optimal solution by mimicking Charles Darwin’s natural evolution theory. Genetic algorithm mimics how genetic information are passed over through generations. Like DNA is formed with genes, a gene is an element that forms a chromosome. A chromosome is a set of genes that forms a potential solution to the problem. The fitness function evaluates how well each chromosome solves the problem. The scores calculated are used for selection, which is a process of selecting chromosomes that will create an offspring. The higher the score is, the higher chance it has to be selected. Like our DNA, an offspring is created by combining the genetic information of two chromosomes. The point where that crossover will be made is randomly selected. For example, aaaa and bbbb can become aaaa, abbb, aabb, abbb,and bbbb. Mutation happens randomly at each position, changing the value from 0 to 1 or 1 to 0.

The genetic algorithm follows the flow chart Figure 4:

Genetic Algorithm For Selecting Optimal Channels for P300 Based Classification ^[26,27]

EEG signals can be roughly divided into four categories: Visual Evoked Potential (VEP), Sensorimoter Rythm (SRM), Slow Cortical Potentials (SCP), and Event-Related Potentials (ERP). VEP is a signal acquired from three occipital electrodes, while forming two channels by setting the middle electrode as a reference, and represents the pathway of visual activity, from optic nerve to primary optic cortex. SRM is a signal recorded from sensorimotor cortex at mu(7–13 Hz), beta(13–30 Hz), and gamma(30–200 Hz) frequency bands, and is strongly modulated with simple movements and certain cognitive tasks. SCP is a measurement of a slowly changing electrical signal with frequency lower than 1Hz. ERP is not a term referencing an EEG measured in specific location, but is a fluctuation generated by neurons due to some event, such as motor or speech imagery, at a specific time. [26,28,29,30,31,32,33,34].

P300 is an unique ERP, which appears as a positive peak near 300ms after stimulation happens. A 64 channel EEG was recorded and four channels were selected in [26], by defining the fitness function ${\displaystyle f\left(X\right)=\sum _{k=1}^n{A}_K}$ where X is a chromosome consisting four genes, which are channels in this case. $A_k$ is defined as the maximum $r^2$ among four channels, where $r^2$ is a value between 0 and 1, and higher the $r^2$ is, the higher the variance among classes are.

5. Result

In this study, we utilized the EEG dataset from BCI Competition IV Dataset 2a, which is specifically designed for benchmarking motor imagery classification algorithms. The dataset comprises EEG recordings from nine healthy subjects, each performing four different motor imagery tasks: imagining movements of the left hand, right hand, both feet, and tongue. EEG signals were recorded using 22 Ag/AgCl electrodes placed according to the international 10-20 system at a sampling rate of 250Hz [35].

For preprocessing, we applied a band-pass filter between 8Hz and 30Hz, a frequency band commonly used in previous research[11,12]. Baseline correction was performed by subtracting the mean of a 2 second pre-cue interval, mentioned as "Fixation cross" (from -2s to 0s relative to the cue onset), from the subsequent motor imagery data. The entire 4seconds of motor imagery data (from 0s to 4s after the cue onset) were used for classification.

Data were split into training and evaluation datasets using 5-fold cross-validation, ensuring that the classes were balanced in terms of the number of trials in each fold. Cross-validation was performed individually for each subject to account for inter-subject variability of EEG signals.

The accuracies of each algorithm were measured by training a support vector machine (SVM) classifier with a linear kernel using the data from the selected channels.

5.1. CSP and SCSP Algorithm

CSP algorithm can be thought of as a special case of the SCSP algorithm where $r=0$. The accuracy of SCSP algorithm depends on its r value, so accuracy of SCSP algorithm where $r=0,0.1,0.2,0.3,0.4,0.5$ is compared as Figure 5.

5.2. Accuracy of CSP Based Channel Selection Methods

The accuracy of CSP based channel selection methods: CSPRank(Section 3.2), L1 Norm applied CSPRank(Section 3.3), SCSPRank(SCSP used instead of CSP in CSPRank), E-CSP(Section 3.4), and Bhattacharyya bounded CSP(Section 3.1) is as following Table 1 and Figure 6.

Table 1. Comparison of Classification Performance for Various CSP Based EEG Channel Selection Methods Across Different Numbers of Channels.

Table 1. Comparison of Classification Performance for Various CSP Based EEG Channel Selection Methods Across Different Numbers of Channels.

5.3. Accuracy of Non-CSP Based Statistical Channel Selection Methods

The accuracy of Non-CSP based channel selection methods: SVM-RFE(Section 4.1), ${\gamma }^2$ value based channel selection(Section 4.2), Correlation based channel selection(Section 4.3), and Genetic Algorithm based channel selection(Section 4.4) is as following Table 2 and Figure 7.

Table 2. Comparison of Classification Performance for Various Non-CSP Based Channel Selection Methods Across Different Numbers of Channels.

Table 2. Comparison of Classification Performance for Various Non-CSP Based Channel Selection Methods Across Different Numbers of Channels.