An Approach of One-vs-Rest Filter Bank Common Spatial Pattern and Spiking Neural Networks for Multiple Motor Imagery Decoding

Motor imagery (MI) is a typical BCI paradigm and has been widely applied into many aspects (e.g. brain-driven wheelchair and motor function rehabilitation training). Although signiﬁcant achievements have been achieved, multiple motor imagery decoding is still unsatisfactory. To deal with this challenging issue, ﬁrstly, a segment of electroencephalogram was extracted and preprocessed. Secondly, we applied a ﬁlter bank common spatial pattern (FBCSP) with one-vs-rest (OVR) strategy to extract the spatio-temporal-frequency features of multiple MI. Thirdly, the F-score was employed to optimise and select these features. Finally, the optimized features were fed to the spiking neural networks (SNN) for classiﬁcation. Evaluation was conducted on two public multiple MI datasets (Dataset IIIa of the BCI competition III and Dataset IIa of the BCI competition IV). Experimental results showed that the average accuracy of the proposed framework reached up to 90.09% (kappa: 0.868) and 81.33% (kappa: 0.751) on the two public datasets, respectively. The achieved performance (accuracy and kappa) was comparable to the best one of the compared methods. This study demonstrated that the proposed method can be used as an alternative approach for multiple MI decoding and it provided a potential solution for online multiple MI detection.

output their thoughts by only imagining without any real movements and speaking, as well as realising practical operations [4].
In the establishment of a BCI system, feature extraction is a crucial step. Common spatial pattern (CSP) is a widely used method to extract spacial features as it effectively constructs the best spatial filter for differentiating two classes of motor imagery. As CSP searches the best spatial filter by considering temporal dynamics, it depends on the information in the temporal domain and is sensitive to temporal noise. Common spatial spectrum pattern (CSSP) was proposed to reduce the interference of such noise [5]. However, CSSP cannot offset the increased complexity of optimization problems because it increases the flexibility of the frequency filter with the delay taps. Afterwards, an improved algorithm was proposed, naming Common Sparse Spectral Spatial Pattern (CSSSP) [6].
Before applying CSP, a bandpass filter is usually used to concentrate on a specific frequency band. The specified frequency band largely affects the performance of CSP. In order to have good performance, frequency band is manually specified for each individual subject. Alternatively, a few frequency bands could be specified to extract features in different bands at the same time. Filter bank common spatial pattern (FBCSP) is one of such methods, extracting discriminative spatio-temporal information [7], [8]. Its variants were subsequently proposed, which are discriminative FBCSP (DFBCSP), sparse FBCSP (SFBCSP) and penalized time-frequency band CSP (PTFBCSP) [8]- [10]. These methods have different advantages. For example, DFBCSP allows to set subject-specific frequency band for filtering based on Fisher's ratio and has a better performance compared to FBCSP. SFBCSP automatically selects features in a few frequency bands based on the regression of the least absolute shrinkage and selection operator (LASSO). PTFBCSP overcomes the limitation of the fixed time period and fixed frequency band.
Spiking neural networks (SNN) has obtained a good classification effect as a classifier in the digital recognition, EEG-based motor imagery, and other classification problems [11]- [13]. Virgilio et al. used the original EEG signal, power spectrum density (PSD) and discrete wavelet transform (DWT) combined with SNN classifier to analyze and study the multiple-classes of motor imagery, and achieve a good classification effect [14]. Carino et al. conducted a taxonomic study by evaluating two different neuronal model coding strategies for two classes of motor imagery time-frequency characteristics in patients with brain death and healthy subjects [15]. Niranjani et al. used the spiking neural classifier based on the Online Meta-neuron based Learning Algorithm (OMLA) to classify MI. Due to the simultaneous use of global and local information of the network, the performance of the spiking neural classifier is better than other classifiers [16]. Salazar et al. demonstrated that the classification performance of spiking neural models (SNM) is a potential choice for EEG when training with fewer data samples [17]. Though the SNN has been applied in motor imagery analysis in the above literatures, the main contribution of this paper is to combine FBCSP with the SNN classifier for multiple MI decoding. The performance was evaluated using the dataset IIIa in the BCI Competition III and the dataset IIa in the BCI Competition IV [18], [19].

II. METHODOLOGY AND MATERIALS A. DATA COLLECTION
We have evaluated our method on the dataset IIIa of the BCI Competition III and dataset IIa of the BCI Competition IV [18], [19]. More details of these multiple-class MI datasets are listed in Table 1. The channel selection for two datasets is based on the selection of channel location in [20]- [22]. These channels are selected mainly because they account for a higher proportion of MI information than other channels that have not been selected, and the channel selection range includes the brain activity areas of the four-classes of motor imagery. Choosing a small number of channels but a large proportion of the weight is beneficial to the subsequent processing speed of the entire experiment and improves the accuracy of multiple-classes of MI.
More specifically, the location of electrodes montage for two datasets are shown in Fig. 1 (red dots stand for the selected channels) and the numbers of channels for EEG collection are 60 and 22, respectively. For the dataset IIIa in BCI Competition III, there are three subjects and the paradigm consists of four MI tasks, namely left hand (class 1), right hand (class 2), both feet (class 3), and tongue (class 4). Each category has 90 trials for the subject K3b while 60 trials for the subjects K6b and L1b. For the Dataset IIa in the BCI Competition IV, there are nine subjects and four categories (left hand, right hand, foot and tongue). Each subject underwent two sessions, which comprised 6 runs separating by a short break. Each run had 48 trials with equivalent numbers for each class.
To analyze datasets, the experimental algorithm is run by using an Intel(R) Core(TM) i5 CPU @2.7 GHz host computer (Lenovo).
In the preprocessing, the spatial filtering of MI data adopts common average reference (CAR) [29]. In band pass filtering, we divided the frequency band into 2, 6, 10, 12 frequency bands with different ranges for band pass filtering. The value ranges of two frequency bands are 7-14 Hz and 14-30 Hz, and 6 frequency bands are 7-12 Hz, 12 First, we define the four classes of MI tasks as x 1 , x 2 , x 3 , and x 4 , respectively. Then the covariance model of multiple-class MI space is shown as following [29]: where x i represents the trial of each class in the entire dataset. The dimension of x i is the number of channels × (time window × sampling frequency). By accumulating each class of experiments in the training set by (1), R i obtained is the sum of the covariances of each class of trial in the training set. The compound covariance matrix formula for CSP is as follows: The singular value decomposition of covariance matrix R x can be carried out as: where U 0 is the unitary matrix of principal components, and C is the diagonal matrix of eigenvalues. After calculating singular value decomposition of eigenvector and eigenvalue matrix, the transformation matrix of covariance matrix can be obtained: For multiple-class MI, multiple-class of OVR-CSP modes (such as left hand) can be calculated as: Then we transform R 1 and R 1 into: Furthermore, we do eigenvalue decomposition for G 1 and G 1 : By combining the equation (1)-(7), we can get: We get the common eigenvector matrix U for R 1 and R 1 . Further, CSP projection matrix V 1 = U T 1 P can be obtained. One trial of the EEG data matrix X i projection is obtained. The selected features of left-hand MI can be obtained by the following formula: M2, M3, and M4 can be obtained in the same way. In this paper, the first four pairs and the last four pairs of eigenvalues in the projection matrix were selected, and CSP filters are used for spatial filtering to form eight eigenvalues in each trial of the training dataset [30].
We carried out the aforementioned feature extraction for each frequency band after preprocessing. Finally, N frequency bands for each trial need N × 8 CSP filters for feature extraction. Eight eigenvalues extracted from each frequency band are merged to form N ×8 eigenvalues. After the features obtained by the CSP, we uniformly normalize them to better perform feature selection and classification.

D. FEATURE SELECTION
The F-score algorithm has a simple but effective advantage for evaluating the discriminative power of each feature in the feature set [31], [32]. F-score is used to sort and select the features. With this method, the features with the best classification effect are selected from the classification order. The optimal feature set that can achieve the best accuracy of classification is the optimal feature set selected by F-score.
Given train feature dataset samples X k , k = 1, 2, . . . , n. The number of positive and negative samples is N + and N − , respectively. Then, the F-score of the f th feature of the dataset is defined as [33]: where,

E. CLASSIFICATION
SNN is proposed based on biological principles, which could simulate the connection and communication between neurons to the great extent. The spike neuron model is the mathematical abstraction of the real neuron, which shows good performance and strong robustness in some pattern recognition problems [13], [17], [34], [35]. Due to such advantage, SNN is used in this study to classify the optimal feature set of MI.
Since the Izhikevich (IZ) model has good performance between biological accuracy and calculation cost, it is selected for the realization of the neuron model. The description of the IZ model is illustrated by the following equations [15], [36]: where k represents the rheobase resistance, v r represents resting membrane potential, v t represents the instantaneous threshold potential, C and v respectively represent the membrane capacitance and the membrane potential, while v peak represents the spike cutoff value, u represents the recovery current, I represents a vector of the input current to the neuron, a represents a recovery time constant, b represents the input resistance, the voltage reset value is expressed by c, and d represents the outwards minus inwards currents that affect the after-spike behavior of the model during the spike period [14], [15]. The parameters used in the SNN are shown in Table 2. The detailed structure of the model can be described in [37].
To produce the desired behavior in the output of the spiking neuron, we must adjust the synaptic weights of the model. In this regard, we use the cuckoo search algorithm as a learning strategy to train the spiking neural model, which corresponds to the training stage of the spiking neuron [36]. In this algorithm, each egg in the nest represents a solution, while the cuckoo egg represents a new solution. Its goal is to use new and potentially better solutions, not bad ones in the nest. The detailed description of the algorithm can be found in [38].
When a new solution X (t+1) is generated for a cuckoo i, a Levy flight is performed, The equation as follows: where α > 0 is the step size. Usually, we can use α = 1. The product ⊕ represents entrywise multiplications. The Levy flight essentially provides a random walk, and the random step size is obtained from Levy distribution: which has an infinite variance with an infinite mean. In the long run, the random walk via Levy flight is more effective in exploring the search space because of its longer step length. Because this kind of walking is more effective in exploring search space, it can regulate the synaptic weights of spike neurons [36]. For the problem of MI recognition, let the optimal CSP feature set S = {X i , n} p i=1 selected by feature selection be the p input patterns, where n = 1, . . . , N is the class to of VOLUME 8, 2020 X i ∈ IR k . First, we need convert each input mode to the current input I . It should be noted that the spike neuron model is stimulated by the injection current I calculated from the input pattern, rather than by the input pattern X i ∈ IR k . Since synaptic weights of the model are directly connected to the input pattern X i ∈ IR k , the injection current generated by this input mode can be calculated as [36]: where W i ∈ IR k is the synaptic weights set of the neuron model, and γ is a gaining factor that helps the neuron fire. This kind of conversion can convert multiple input modes from the same class into the same or similar current, and help neurons generate similar firing rates. Then during T ms, the spiking neuron is stimulated with the input current I , and the firing rate is calculated. These steps are used in each input pattern. Calculate the average firing rate of each class AFR ∈ IR N by obtained the firing rate of each input. If the spiking neuron is trained, the class used to determine the unknown input pattern X belongs will be determined by the firing rates generated by each input pattern. As described in the following equation: where the firing rate of the neuron model stimulated by input pattern X is fr. It is very important to find a set of optimal the neuron model ← − w synaptic weights to improve the accuracy of the spiking neuron. This needs to be adjusted by using the cuckoo search algorithm. Therefore, the current behavior needed to generate the model can be realized by finding a set of synaptic weights through the following fitness function: where ← − w represents the synapse of the model, s is the optimal CSP feature set as a set of input patterns, and performance( ← − w , S) is a function of identifying the accuracy of the number of correct classified divided by the total number of tested.

F. PERFORMANCE MEASURE
The measurement performance used in this paper is classification accuracy and kappa value, respectively. The classification accuracy can be expressed by the following formula [24]: where C is the number of correctly classified samples of the test set. N represents the total number of test set samples. Kappa is often used to measure multi-class problems, which can be expressed as [14]: where p 0 is the overall agreement of the test set, which is equal to the value of accuracy. p e is the chance agreement probability value of the test set, which can be obtained by the following formula: where a i and b i represent the sum of ith class real samples and ith class predicted samples of the confusion matrix, respectively. N is the total number of test set samples.

III. EXPERIMENTAL RESULTS
To assess the impact of the frequency bands on MI-BCI performance, we evaluated the CSP and FBCSP methods for the following cases: 1) CSP: For the time window optimization problem of two datasets, it has been studied in [20], [39], [40]. On this basis, we manually select the optimal time window of each subject. The topographical maps of FBCSP weights for K3b and A01 are shown in Fig. 2. The color bars indicate changes in weights. Red represents a high weight in the topographic map, and blue represents a low weight in the topographic map. When the subject performed the left-hand motor imagery, the weights on the right brain hemisphere were high (in red color). When the subject performed the right-hand motor imagery, the high weights were on the left brain hemisphere. When the subject performed the foot motor imagery, the high weights appeared on the central brain region. When the subject performed the tongue motor imagery, the high weights were shown around the outer brain region. It can be seen that the activation intensity and topographic map distribution are different between different tasks. These topographic map distributions verify the rationality of the subject's CSP projection matrix neurophysiology.
Since motor imagery of each limb part will cause the weight of the non-motor imagery cortex to decrease, CSP has a greater weight in the large area of the sensory motor cortex. Some areas that are not close to the sensory motor cortex also have larger weights, which may be caused by some artificial factor.
To date, most of researches on the MI classification addressed feature extraction and made relatively less attention on the development of classifier. Also, fixed frequency bands, fixed time windows, and all channels were used, which is time-consuming and is not optimal for the MI classification. In this paper, the OVR-FBCSP was used to extract the features of different time frames of different frequency bands of two datasets, and three classifiers (support vector machine (SVM), relevance vector machine (RVM), and SNN) were used for feature classification. Table 3 showed that the accuracy of OVR-FBCSP-SNN was superior to the other methods in the assessment of the  dataset IIIa of BCI competition III for all subjects except the subject K3b. On average, the highest mean accuracy (90.09%) was achieved by OVR-12FBCSP-SNN. In the assessment of the dataset IIa of BCI competition IV, SNN was better than SVM and RVM in the overall classification of nine subjects. The highest accuracy of OVR-SNN was higher than that of OVR-SVM and OVR-RVM except for A01, A03, A05, and A08 subjects. In particular, the highest average accuracy of OVR-12FBCSP-SNN was 81.33%. 86856 VOLUME 8, 2020  In the assessment of the dataset IIa of BCI competition IV, the overall kappa values of SNN in nine subjects were better than that of SVM and RVM. The highest average kappa value of OVR-12FBCSP-SNN was 0.751. Table 5 showed the training and testing time of two datasets in different frequency band selection range. We can see with the increase of frequency bands number, the corresponding training and testing time also increases.

IV. DISCUSSION
In recent years, great progress has been made in the analysis of MI based on EEG signals, and the mechanism of different MI has gradually been revealed. In this paper, we explored the combinations of FBCSP and SVM, RVM, SNN for multiple-class MI classification. First, we used different time windows according to different subjects. And the corresponding band-pass filtering is carried out for each subject according to the same sub-band selection. Then, we extracted the CSP features in each frequency band and selected features. Finally, different classifiers were used to classify MI. Table 3 and Table 4 showed that the OVR-FBCSP-SNN was effective and feasible. These results were discussed in details below.

A. TIME WINDOWS FOR SPATIAL PATTERNS
A fixed time window (0.5-2.5 s with respect to the onset of cues) was usually used in the CSP for feature extraction [44], [45]. This is not very reasonable because the temporal course of brain activity responding to motor imagery is not exact the same and varies from one subject to another. In this study, the optimal time window was manually set for different subjects so that the classification accuracy was improved. For the dataset IIIa of the BCI Competition III, the optimal time windows of K3b, K6b, and L1b are 0.5-4.0 s, 0.5-3.6 s, and 0.5-3.6 s, respectively. For the dataset IIa of the BCI Competition IV, the optimal time windows of A01-A09 are 0.5-3.7 s, 0.5-3.7 s, 0.5-3.9 s, 0.5-3.5 s, 0.5-3.9 s, 0.5-3.8 s, 0.5-3.8 s, 0.5-3.4 s, and 0.5-3.5 s, respectively.

B. DIFFERENT FREQUENCY RANGE SELECTION
Most of studies focused on the band selection of CSP filters [54]- [56]. Through the average division of different frequency bands, the CSP feature extraction is carried out for each sub-band, and then the optimal CSP feature combination in each sub-band is selected through feature selection. The VOLUME 8, 2020 optimal CSP feature combination is classified to get a good classification effect.
The experimental results showed that the FBCSP combined with SNN was an effective method to improve the accuracy of MI classification [27]. Ang et al. used the FBCSP combined with various feature selection and classifiers to identify the two-class MI. The results showed that the FBCSP is superior to other feature extraction methods [7]. Mahmood et al. used the FBCSP combined with SVM to classify multiple-class MI. Compared with other methods, this method obtained the best recognition accuracy [20]. Zhang et al. used OVR-FBCSP combined with CNN+LSTM method for MI recognition [53]. Islam et al. used TSM to find the precise frequency band associated with the MI mission for MI recognition [41]. In their work, the optimal MI recognition accuracy based on the sub-band selection method of mutual information is 82.6%±13.06, which is lower than our experimental result 90.09% and proves that the effect of FBCSP is feasible for MI recognition.

C. DIFFERENT CLASSIFIERS IN MOTOR IMAGERY
For the classification of MI features, the classifier has a great influence on the classification effect. In this study, we use SVM, RVM, and SNN as classifiers and apply the feature selection method of F-score to screen the most distinctive CSP features in each subject, to help improve the accuracy of recognition. For the dataset IIIa of the BCI competition III and dataset IIa of BCI competition IV, the average accuracy of the OVR-FBCSP-SNN classifier is 90.09% and 81.33%, respectively. This result proves that SNN is a good classifier based on EEG, and our classification method is feasible and effective.
Referring to other MI literature, Liu et al. used features to define a score use the step-wise linear discriminant analysis (SWLDA) method performing MI classification [57]. Nguyen et al. used CSP for feature extraction and fuzzy logic system (FLS) for classification. CSP is used to extract significant features, and then these features are input into FLS as classification input [58]. G. S. et al. used Bayesian networks and artificial neural networks to classify MI features [59]. Yang et al. designed a classification method called adaptive kernel fisher support vector machine (KF-SVM) is designed and applied to EEG MI classification in BCI [60]. Komijani et al. presents MI classification for BCI systems using a recurrent adaptive neuro-fuzzy interface system (ANFIS), and the classification system is based on time-series prediction [61]. Lahiri et al. proposed to use the whole classifier composed of the k-nearest neighbor (KNN) layer for classification [62].

D. COMPARISON OF DIFFERENT MOTOR IMAGERY RELATED STUDIES
Obviously, compared with other recent methods (Table 6 and  Table 7), our proposed OVR-FBCSP-SNN method has some advantages. The accuracy and kappa value of the two datasets are close to the highest value of the latest literature. The 86858 VOLUME 8, 2020 results show that the framework can be used as an alternative method.

E. EXPERIMENTAL LIMITATIONS AND FUTURE CONSIDERATIONS
First of all, the dataset of MI is limited to the left hand, right hand, foot, and tongue, which cannot reflect other situations of MI. For example, grabbing, lifting, and other forms of MI are not included in the dataset. In the future, we will realize the recognition of MI with more classes of MI. Secondly, we select the optimal time window and frequency band of each subject manually, which is not an effective method in practical application. It is urgent to improve the use of automatic selection of each subject's optimal time window and frequency band, which is efficient and accurate. Thirdly, we only use SNN as classifier without considering whether it can be used as a feature extraction method.

V. CONCLUSION
To enhance the performance of multiple-class MI classification, a new method of combining one-vs-rest filter bank common spatial pattern with spiking neural networks was proposed. Two datasets of four classes of MI (dataset IIIa of BCI competition III and dataset IIa of BCI competition IV) were used to evaluate the proposed method. Experimental results demonstrated that our method was comparable to the best existing method in terms of classification accuracy. It demonstrated that the proposed method could serve as an alternative for multi-class MI classification. He is currently a Full Professor with the Faculty of Intelligent Manufacturing, Wuyi University, and been selected as a Thousand-Hundred-Ten Talent of Universities in Guangdong Province. He is also the Director of the Jiangmen Brain-like Computation and Hybrid Intelligence Research and Development Center. From January 2017 to January 2019, he was a Visiting Research Fellow with the Centre for Life Sciences, National University of Singapore. His current research interests include the fields of brain-like computation, pattern recognition, deep learning, and hybrid intelligence.  Since July 2006, he has been a Lecturer with the School of Information Engineering, Wuyi University of Technology. His research interests include measurement, signal processing, and automatic control.

CONG TANG was born in
JUNHUA LI (Senior Member, IEEE) received the Ph.D. degree in computer science from the Department of Computer Science and Engineering, Shanghai Jiao Tong University, China, in 2013.
He is currently a Lecturer with the School of Computer Science and Electronic Engineering, University of Essex, Colchester, U.K. He was a Senior Research Fellow with the National University of Singapore, Singapore. His research interests include brain-computer interface, computational neuroscience, data analytics, and machine learning.
ANASTASIOS BEZERIANOS (Senior Member, IEEE) received the degree in physics from Patras University, the degree in telecommunications from Athens University, and the Ph.D. degree in bioengineering from the University of Patras. He is currently the Head of the Cognitive Engineering (COGEN) Laboratory, N.1 Health Institute, National University of Singapore, a Professor with the NUS Graduate School for Integrative Sciences and Engineering, and a Visiting Professor with the Computer Science Department, New South Wales University (NSWU), Canberra, ACT, Australia. He has been a Professor of medical physics with the Medical School of Patras University, Patras, Greece, since 2004. His research entails diverse areas spanning from artificial intelligence and robotics to biomedical signal processing and brain imaging as well as mathematical biology and systems medicine and bioinformatics. His work is summarized in 140 journals and 217 conference proceedings publications, one book, and two patents. He has research collaborations with research institutes and universities in Japan, China, and Europe. He is also a Fellow of the European Alliance for Medical and biological Engineering and Science (EAMBES), and the Founder and the Chairman of the Biannual International Summer School on Emerging Technologies in biomedicine. He is also an Associate Editor of the IEEE TRANSACTIONS ON NEURAL SYSTEMS AND REHABILITATION ENGINEERING (TNSRE), the PLOS ONE, and Neuroscience Journal and a Reviewer for several international scientific journals. He is also a Registered Expert of the Horizon 2020 Program of the European Union and a Reviewer of research grant proposals in Greece, Italy, Cyprus, and Canada. VOLUME 8, 2020