Electroencephalography-Based Brain-Computer Interface Motor Imagery Classification.

Background: Advances in the medical applications of brain-computer interface, like the motor imagery systems, are highly contributed to making the disabled live better. One of the challenges with such systems is to achieve high classification accuracy.
Methods: A highly accurate classification algorithm with low computational complexity is proposed here to classify different motor imageries and execution tasks. An experimental study is performed on two electroencephalography datasets (Iranian Brain-Computer Interface competition [iBCIC] dataset and the world BCI Competition IV dataset 2a) to validate the effectiveness of the proposed method. For lower complexity, the common spatial pattern is applied to decrease the 64 channel signal to four components, in addition to increase the class separability. From these components, first, some features are extracted in the time and time-frequency domains, and next, the best linear combination of these is selected by adopting the stepwise linear discriminant analysis (LDA) method, which are then applied in training and testing the classifiers including LDA, random forest, support vector machine, and K nearest neighbors. The classification strategy is of majority voting among the results of the binary classifiers.
Results: The experimental results indicate that the proposed algorithm accuracy is much higher than that of the winner of the first iBCIC. As to dataset 2a of the world BCI competition IV, the obtained results for subjects 6 and 9 outperform their counterparts. Moreover, this algorithm yields a mean kappa value of 0.53, which is higher than that of the second winner of the competition.
Conclusion: The results indicate that this method is able to classify motor imagery and execution tasks in both effective and automatic manners. Copyright:

Introduction

The idea of the interaction of minds with machines is a major concern in the human imagination. Brain–computer interface (BCI) is a system where brain signals control the functionality of a given device.[1] BCI is commonly applied in security, lie detection, alertness monitoring, telepresence, gaming, education, art, human augmentation, and virtual reality apparatus.[2] Applications of BCI in recent advances made in medicine set the stage for restoring and potentially upgrading human physical and mental capabilities. Some of such medical applications consist of deep brain stimulation in Parkinson's disease, bypassing the motor difficulties in amyotrophic lateral sclerosis patients, and cochlear implants.[2]

There exist many methods in extracting the brain signal features such as electrocorticography (ECoG), magnetoencephalography (MEG), functional magnetic resonance imaging (fMRI), functional near-infrared spectroscopy (fNIRS), and electroencephalography (EEG). Since intracortical implants and ECoG are invasive methods, they are not widely employed adopted. The methods MEG and fMRI are not economically feasible. The spatial and temporal resolution of fNIRS is moderate and it is more noise resistant than EEG. Since EEG has appropriate temporal resolution and the device applied therein is of low cost, it has fewer environmental restrictions, making it one of the most applicable methods in this field, although its spatial resolution is low.[3] These advantages make the EEG signal recording the basis of many BCI systems. One of the important branches of the EEG-BCI systems is the motor imagery where the EEG signals’ recording occurs when the subject is imagining a special action (e.g., movement of the hand, foot, or other body organs).

Figure 1 displays a generic pipeline of a motor imagery-based EEG-BCI system. According to it, the operation of a motor imagery system consists of four steps. First, the EEG must be acquired. Next, in the pre-processing step, some filtering and artifact removing techniques are applied to the recorded signal. Then, the system extracts discriminative features applied in training an appropriate classifier, and finally, the system identifies the movement that the subject has imagined.[4]

Figure 1

The schematic of a general Motor Imagery EEG-Brain–Computer Interface system

The recorded EEG signal is blurred due to volume conduction and is likely to be contaminated with various artifacts. Therefore, the EEG signal preprocessing and feature extraction are necessary to represent the input signal in its modified state in reduced feature space which, in turn, improve the motor imagery pattern recognition.[5] The feature extraction methods may be adopted in different spatial, temporal, frequency, and hybrid domains, depending on a variety of criteria.

In the spatial domain, common spatial pattern (CSP) was first proposed in multi-channel EEG classification. CSP is frequently applied in BCI applications to portion a multivariate signal into subcomponents with a maximum difference in variance between two classes. This technique is contributive to classifying the unknown data in the classification stage more effectively.[6789] There exist different CSP-based algorithms which improve accuracy by increasing the complexity of the algorithm.[1011121314] Moreover, temporal features have widely been and are being applied in EEG analysis. In some studies,[1516] the Hjorth parameters of the EEG signals are applied for analyzing its signals. In some studies,[17] to improve classification, time features and signal derivatives are assessed. In the frequency domain, different frequency features such as power spectral density[18] and bispectrum[19] are applied in EEG analysis. For the hybrid method, feature extraction is made by applying different temporal-spatial-frequency representations like wavelet packet decomposition,[2021] time–frequency approach,[222324] and CSP.[25]

The major concern in the BCI system is its performance in terms of accuracy and robustness.[262728293031] The obtained results so far indicate that selecting appropriate features may be more effective than selecting an appropriate classifier in the ultimate performance of a BCI system.[3233]

The objective of this article is to obtain a high-accuracy algorithm while benefiting from low computational complexity. The contributions here consist of applying CSP for the dimensional reduction in the preprocessing step, extracting appropriate temporal and time–frequency features, and selecting a proper classifier to obtain a high-accuracy algorithm. This article is organized as follows: the method is introduced in Section 2; the results are expressed in Section 3 and the article is discussed and concluded in Section 4.

Subjects and Methods

The method applied here is analytic.

Data descriptions

First Iranian brain–computer interface competition dataset

The data are cumulated by the National Brain Mapping Lab for the First Iranian BCI Competition (iBCIC).[3435] The EEG signals of 7 healthy right-handed individuals (3 women and 4 men) with a 31 years age average are recorded. A 64-channel signal recording system with a sampling frequency of 2400 Hz is applied. Channel 33 is the reference electrode connected to the right ear, while electrode 3 (Fpz) is grounded.

To cumulate EEG data from the subjects, first, the EEG signals of the subject in the relaxation and closed eyes state are recorded for 2 min followed by 2 min for the relaxation state with open eyes, after which the experiment protocol begins. The subject sits in a comfortable chair and a monitor within a half-meter distance and is told to react to perform motor execution and imagery tasks according to the signals in the monitor. The limbs engaged are the thumb of the right hand, the right leg, and the right arm. First, the “+” sign is shown on the monitor for 2 s, during which the subject should not think about anything and only wait for the cue. In the following 2 s, the cue is shown, and after another 2 s, it disappears and the word “Go” appears, upon which subject has 3 s to perform the task, Figure 2.

Figure 2

Experiment protocol

Each cue is shown twenty times during signal recording in a random manner. The cue indicating the imagination of a particular task follows the execution process. Both the motor execution and images are run twice for each subject with a 5-min interval.

Data formats

The data provided here consist of 3-s trials, collected from 64 channels. Motor imagery and motor execution data are in separate files with their own training and test data. There is a training data file containing data related to 40 trials for each task of each subject. There exist four classes for motor execution: (1) right arm; (2) right foot; (3) right thumb, and (4) no motor execution. There exist three motor imagery tasks: (1) right arm, (2) right thumb, and (3) right foot. There exist nearly 20 trials for each task in the test data when the participators were required to recognize the true class of each piece of data related to each trial for both motor execution and motor imagery data. The correct test labels are published after the competition.

Brain–computer interface competition IV data sets 2a

To test the generality and effectiveness of this proposed algorithm, it is applied to BCI Competition IV Dataset.[36] Data sets 2a is one of the subsets of the competition dataset including motor imagery data of 9 subjects in 4 classes. The EEG data are extracted from 22 channels with a 250 Hz sampling rate which are first filtered with a notch filter to remove electrical line noise and next followed by a low-pass filter allowing only frequencies between 0.5 and 100 Hz to pass through.

There exist two files: training and evaluation: First, the data are labeled as (1) imagining the left-hand movement, (2) imagining the right-hand movement, (3) imagining the right-foot movement, and (4) imagining the movement of the thumb of the right hand. There are 72 trials for each class, thus 4 × 72 = 288 trials for each subject. The competitors are required to determine the class of each trial in the evaluation dataset. The procedure here corresponds to that of some studies.[36]

Preprocessing

CSP is a practical method frequently applied in BCI applications to proportionate a multivariate signal into subcomponents with maximum differences in variance between two classes. This allows for better separability of two classes due to EEG application on the difference between the power of certain bands and the corresponding variances.

To show how CSP[37] works, assume x(t) is a matrix of signals with the N × T dimension, where N is the count of channels and T is the count of samples and the output of CSP transform will be:

x (t) =W * x (t) (1)

Where, WT is the transpose of W, which is an N × M matrix where every column is a special filter. The objective of CSP is to determine these spatial filters in a manner that the variances of one output to the other notion will be maximized. For this purpose, consider is the data matrix representing the data from class and the trials from i=1 to K.. Each trial is an N × T matrix where N is the count of channels and T is the count of samples of the trial. To obtain the spatial filters, the covariance of each class is estimated through Eq. (2):

W1W = M1 (2)

Consequently, W would be computed through the following equations:

W1W = M1 (3)

W2W = M2 (4)

M1 + M2 = I (5)

In Eq. (5), I is the identity matrix. Therefore, the problem will be converted into an eigenvalue problem expressed through Eq. (6):

R1 W = λR2W (6)

The eigenvalues consist of: and , where wj is the jth eigenvector. The eigenvalue and are the elements of M1 and M2, and , respectively. Thus, when, is high, is low, and the variance of the output for one class will be higher than that of the other, thereby, an increase in the separability.

In this proposed algorithm, by applying CSP, 64 channels are converted into 4 subcomponents in the iBCIC dataset and 22 channels are converted into 4 subcomponents in Data sets 2a to allow the occurrence of separability between classes. Since many features will be extracted from the EEG data of each subcomponent in the following processes, the complexity of this method will increase if all subcomponents are applied, while in studies regarding EEG, only a few eigenvectors are applied. Thus, in both the datasets, only 4 eigenvectors (the first two and the last two) are selected because eigenvectors from the top and the bottom cause the highest separation between classes. In this context, channel selection can be made through CSP, and here, the method introduced by[38] is adopted where a channel score based on L1-norm of WT columns is defined.

Feature extraction

Time domain features

To extract time features, the data are low-pass filtered once with the cutoff frequency of 120 Hz, which yields low results compared to the time when no filtering is done. The features extracted from these 4 sub-components consist of: the mean absolute value of the signal, the signal mean, the signal variance, Higuchi's fractal dimension at kmax = 5, Higuchi's fractal dimension at kmax = 10, permutation entropy at the order 13 and permutation entropy at order 8, and thus, the total number of temporal features for each trial is 4 × 7 = 28 for both datasets.

The fractal dimension of a time series at any given time domain is calculated through Higuchi's fractal. A signal is labeled fractal if any of its sections resemble the whole signal. To express how this fractal dimension is estimated, consider X as a discrete signal having N samples X(1), X(2),…, X(N). By introducing some subsequences with the distance of k, from each other and different beginning points of n = 1, 2,… the expression: is yield.[39]

To calculate the fractal dimension, the length of the signal, L(k), for different values of k is calculated. The length of is estimated as:[39]

Where, J is the biggest integer ≤ (N-n)/k and L(k) is the average value over k sets of L. If L(k) ∝ k−, then the curve is fractal with dimension D. D is located between 1 and 2 where 1 is for simple curves like a line and 2 is for complex curves. For EEG signals, D varies between 1.4 and 1.7.

Permutation entropy is a measure through which the frequency of the occurrence of symbolic patterns or motifs in time series like EEG signals is estimated. In this context, consider X a signal with N samples, if m samples, or a tuple, with the τ time delay between each consequent sample are selected by considering the beginning time at t, the tuple and the count of tuples will be and k = N-(m-1)τ, respectively. In this study τ is 1. Each tuple with m samples can have m! different ordinal patterns or motifs. Six different motifs for a tuple with 3 samples are shown in Figure 3.

Figure 3

Different ordinal patterns

The probability of the i motif is calculated through Eq. (8):[40]

Where, k, s, and p are the count of tuples, i motif, and its probability, respectively. Consequently, the permutation entropy is calculated through the Shannon Eq. (9):

Where, m and ln m! are the order of permutation entropy and the normalization factor, respectively, and H increases with the irregularity of the signal.

Time–frequency features

The features here are extracted in the wavelet domain. Before extracting the features, unlike the temporal features, the signal is filtered through a 150-point low-pass filter with the cutoff frequency at 120 Hz. In this study, datasets 2a are not subject to filtering because the data are already filtered. The Daubechies wavelet 20db is applied to extract time-frequency features from this signal. The wavelet is decomposed up to level 5. At each level, g(n) is applied to extract the approximation coefficients of x(n) or the approximation coefficients of the preceding level, while h(n) is applied to obtain the detail coefficients, consequently, x(n) is decomposed into 6 signals: five detail coefficient sets and the approximation coefficient set at level 5. For each one of the 6 coefficient sets, the following 4 features are extracted: first, sum of the absolute values of the first one-third of the coefficient set; second, sum of the absolute values of the second one-third of the coefficient set; third, sum of the absolute values of the third one-third of the coefficient set; and fourth, sum of the absolute values of all coefficients of the set. Here, there exist 6 subbands for each subcomponent and 4 features for each subband. In general, there will be 4 × 6 × 4 = 96 features for each trial.

Feature selection

Feature selection is made with respect to Wilk's lambda criterion, by applying the stepwise linear discriminant analysis (SWLDA) in the Statistical Package for the Social Sciences (SPSS), IBM Company, Armonk (N.Y., USA), Released 2015. IBM SPSS Statistics for Windows, Version 23.0, IBM Corp. SWLDA is determined through the best linear combination of the features, which yields in the lowest P value.

Classification

The binary classifiers are applied where the majority voting among the results would determine classification. In motor execution classification in iBCIC, there exist 4 possible classes, thus, 6 possible binary comparisons: (1) classification of Class 1 and 2, (2) classification of Class 1 and 3, (3) classification of Class 1 and 4, (4) classification of Class 2 and 3, (5) classification of 2 and 4, and (6) classification of 3 and 4. As to the motor imagery data with three tasks, there exist 3 binary classifications following the sample pattern as that of motor execution. When the count of votes for any given two classes is equal, the decision is made by the classifier corresponding to both the classes. The classifiers applied here consist of (1) support vector machine (SVM), (2) K-nearest neighbors (KNN), (3) Random Forest, and (4) linear discriminant analysis (LDA). LDA yields the best results in almost all cases except in some cases in datasets 2a. The K-fold cross-validation is applied in the classification stage to consider the effect of all the data in the classifier training. In this method, the data are divided into k equal sections and in each iteration, one section is selected as the validation data, and the other k-1 sections are applied as training data. Therefore, each piece of data will be included in the validation data only once and k-1 times in the training data. Accuracy is considered as the ratio of the count of correctly labeled epochs to the count of total epochs. The overall accuracy is the mean of the total k iterations. Cohen's kappa coefficient is a statistical measure of inter-rater reliability, generally thought as a more robust measure than accuracy, because in Cohen's kappa the agreement occurring by chance is of concern.

Results

SWLDA is applied to select the features which can lead to high classification performance. The count of proper features for every two comparable classes and the accuracy obtained through LDA with k-fold cross-validation (k = 20) are tabulated in Table 1.

Table 1

Proper count of features and the accuracy for binary classifications

Subject	Compared classes	Selected features count	Accuracy (using LDA classifier and k-fold cross-validation [k=20])
1	1 versus 2	23	95
2	1 versus 3	46	85
3	1 versus 4	25	100
4	2 versus 3	61	98
5	2 versus 4	15	100
6	3 versus 4	29	100
7	5 versus 6	40	76
8	5 versus 7	78	83
9	6 versus 7	31	69

LDA: Linear discriminant analysis

As observed in Table 1, the classes are (1) right arm motor execution, (2) right leg motor execution, (3) right thumb motor execution, (4) no motor execution (the resting state), (5) right arm motor imagery, (6) right leg motor imagery, and (7) right thumb motor imagery. Here, in general, there exists an indirect relation between Class 4 (rest) and the proper feature count; there exists a direct relation between Class 4 (rest) and accuracy. The classification accuracy of motor execution tasks is better than that of the motor imagery tasks, which is due to the difficulty of the motor imagery tasks and the difference between the natures of source regions of the signals generated by the motor execution and motor imagery tasks. Among the selected features for different classifications, the highest frequencies are related to the wavelet coefficients of level 5. When classifying Class 4 versus other classes, the signal mean is highly contributive, indicating this can be done by applying simpler features.

Classification results for Iranian brain–computer interface competition dataset

To exhibit the performance of this algorithm, different classifiers including Random Forest, KNN, SVM, and LDA are applied, and the best results are related to the LDA classifier because it is compatible with the feature selection method, SWLDA. The outcome of LDA classification for all the motor imagery and motor execution data in the iBCIC dataset, obtained through k-fold cross-validation (k = 20) are tabulated in Table 2, where, as expected, the accuracy of motor execution is higher than that of the motor imagery.

Table 2

Results of k-fold cross-validation in Iranian brain-computer interface competition dataset

Data type	Accuracy	κ
Motor execution	90.36	0.86
Motor imagery	63.10	0.48
Total average	76.73	0.67

Classification results for Iranian brain–computer interface competition test dataset

To label the iBCIC test dataset through this proposed algorithm, the binary classifiers are trained by applying the training data and the optimal features selected by SWLDA, Table 3.

Table 3

Results of the method on the test data of the Iranian brain-computer interface competition dataset

Data type	Accuracy	κ
Motor execution	91.18	0.88
Motor imagery	60.71	0.41
Total average	75.94	0.65

As observed in Table 4, there exists a 23.41% improvement in the total average for this proposed method than the winner of the first iBCIC Competition.

Table 4

Results of the top 6 participants of the competition[3435]

Rank	Execution accuracy	Imagery accuracy	Average
1	0.6807	0.3690	0.5249
2	0.7001	0.3476	0.5238
3	0.7180	0.3190	0.5185
4	0.6490	0.3857	0.5173
5	0.5531	0.3452	0.4491
6	0.2500	0.3238	0.2869

Classification results of brain–computer interface competition Iv datasets 2a

The results of datasets 2a obtained through this algorithm are tabulated in Table 5. As observed in Table 5, in most cases, LDA is the best classifier, while for the iBCIC datasets, LDA is the best for all cases. When compared with subjects 6 and 9, this algorithm outperforms the algorithms in some studies.[36] Moreover, the kappa mean here is higher than that of the second winner. As to subjects 4 and 5, the results here are close to those of the winner.

Table 5

Results for the world brain-computer interface competition IV - datasets 2a

Subject	Optimal number of features for each binary classifications						The best classifier for each binary classifications						Accuracy	kappa
	1 versus 2	1 versus 3	1 versus 4	2 versus 3	2 versus 4	3 versus 4	1 versus 2	1 versus 3	1 versus 4	2 versus 3	2 versus 4	3 versus 4
1	12	5	15	27	18	22	LDA	LDA	LDA	LDA	LDA	LDA	71.43	0.62
2	12	26	11	12	25	12	LDA	LDA	LDA	LDA	LDA	LDA	43.87	0.27
3	11	28	25	13	14	16	LDA	LDA	LDA	LDA	LDA	LDA	77.35	0.70
4	12	4	10	7	10	6	LDA	KNN	LDA	LDA	LDA	RF	59.44	0.47
5	5	8	5	6	10	1	LDA	LDA	SVM	LDA	LDA	LDA	51.74	0.36
6	18	4	7	5	5	12	LDA	LDA	LDA	LDA	LDA	LDA	48.20	0.31
7	15	20	18	16	12	23	LDA	LDA	LDA	LDA	LDA	SVM	69.10	0.59
8	16	7	15	5	12	13	LDA	LDA	LDA	LDA	LDA	LDA	74.48	0.66
9	13	4	17	14	22	6	LDA	LDA	KNN	LDA	LDA	LDA	83.97	0.79
Mean													64.40	0.53

LDA: Linear discriminant analysis, KNN: K-nearest neighbors, SVM: Support vector machine

Channel selection results

The topographic maps of different classifiers of the iBCIC Dataset are shown in Figure 4, where the darker the red color, the more appropriate the channel for classification. Since binary classifiers are applied here, selecting the channels are classifier-specific oriented. The topographic map of classifier 5 (right leg vs rest) of iBCIC Dataset is magnified in Figure 5. The topographic maps of BCI Competition IV Dataset 2a are shown in Figure 6, and the topographic map of classifier 1 (left vs right hand) of BCI Competition IV dataset 2a is magnified in Figure 7.

Figure 4

Channel selection results on Iranian brain–computer interface competition Dataset

Figure 5

Channel selection results on Iranian brain–computer interface competition Dataset, classifier 5 (right leg vs. rest)

Figure 6

Channel selection results on Brain–Computer Interface Competition IV Dataset 2a

Figure 7

Channel selection results on Brain–Computer Interface Competition IV Dataset 2a, classifier 1 (left vs. right hand)

Discussion

In this article, a robust pattern recognition algorithm is proposed to classify the motor execution and imagery tasks in an automatic manner. For this purpose, CSP is applied as a spatial filter to reduce the dimensionality of the data. Due to the nonstationary nature of EEG signals, a combination of temporal and time–frequency features is extracted from the data. SWLDA, a powerful method in feature selection, is applied here. The multi-class classification is converted into multiple one-versus-one classifications which outperforms one-versus-all classification and multi-class classification. The four classifiers of KNN, SVM, LDA, and Random Forest are applied, among which the LDA yields the best results. The separability of the resting-state class from the other motor execution classes is better because the motor execution tasks activate the closely related regions of the brain. The classification accuracy of motor execution tasks is better than that of motor imagery tasks, and this is due to the difficulty of the motor imagery tasks and the difference between the natures of the source regions of the signals generated by the motor execution and motor imagery tasks.

As illustrated in Figure 4, the selected channels are mostly in the vicinity of Cz, which are related to motor areas of the cerebral cortex.[41] Since the tasks in iBCIC are only related to the right limbs, the selected channels are mostly located on the left side of the brain. The selected channels here almost resemble the available findings, indicating the correctness of the data and the data acquisition procedure. In Figure 5, the selected channels are in the proximity of CP3, located close to the primary foot area, which corresponds to that of some studies.[4142] In Competition IV, the dataset 2a recorded the data of 22 channels where the resolution is lower than that of the 64-channel process and both the left and right sides are in Figure 6 and Figure 7, which correspond as both hands are moved. The effectiveness of the channels in the vicinity of Cz here resembles that of some studies.[414243]

This newly proposed algorithm consists of a feature extraction method and a method for classification developed here together with CSP and SWLD. The results indicate that this proposed method is simpler and more effective than that of the competition participators in most cases, with the capability of effectively classifying motor tasks in an automatic manner.

Financial support and sponsorship

None.

Conflicts of interest

There are no conflicts of interest.