Discriminating three motor imagery states of the same joint for brain-computer interface

Local mean decomposition

LMD is a new adaptive time-frequency analytical algorithm first proposed by Smith (2005). When compared with EMD, LMD can better reduce the mode mixing phenomenon and has lower computational complexity and higher decomposition speed (Zhang & Chen, 2017). Using the LMD method, the extended sequence x_i(t) is decomposed into a set of PFs and a residual component, and the expression is as follows: (2)${\displaystyle x_i\left(t\right)=\sum _{j=1}^k{h}_{ij}\left(t\right)+u_i\left(t\right)}$

where h_ij(t) is the jth PF obtained from the i th electrode, and u_i(t) represents the residual component.

Selection of product functions based on the cloud model

The cloud model is an uncertainty transformation model based on fuzzy set theory and probability theory, which can achieve the transformation between a qualitative concept and its quantitative data. It expresses the qualitative concept through the three digital characteristics {Ex, En, He}. Expectation (Ex) is the central value of a concept, entropy (En) represents the randomness of a qualitative concept, and super entropy (He) is the dispersion degree of a concept (Wang et al., 2019). (3)${\displaystyle \mathrm{Ex}=\frac{1}{n}\sum _{i=1}^n{x}_i}$ (4)${\displaystyle \mathrm{En}=\sqrt{\frac{\pi }{2}}\times \frac{1}{n}\sum _{i=1}^n\left|x_i-\mathrm{Ex}\right|}$ (5)${\displaystyle \mathrm{He}=\sqrt{S-{\mathrm{En}}^2}}$

where S is the second-order central moment of x_i, x_i(i = 1, 2, …,n) represents the quantitative value of n cloud-droplets.

En and He can express the degree of complexity. The larger the values of En and He, the more complex the signals, and vice versa. Since the structure of the real EEG components is generally complicated, the En and He of the real PFs are larger than those of false PFs. Therefore, the parameters En and He of the cloud model are used to select the effective PFs. In our study, through analyzing a large amount of experimental data, we selected the PF1 component as the effective component for subsequent processing. The details of the effective PF selection can be found in the Supplemental Information. For the recognition of other data sets, we need to employ the entropy (En) and super entropy (He) of the cloud model to select effective PF components again.

In this research, the experimental data of all channels were decomposed by LMD in turn, and the PF1 components of each channel were recombined to construct a new signal matrix X ∈R^M×N, where M is the order of the selected PFs, and N represents the sampling point of PF. In the next step, CSP was used to extract the spatial features of X.

Common spatial pattern

The goal of CSP is to design a spatial filter which can maximize the variance of two kinds of motor imagery EEG data. Because CSP is based on the simultaneous diagonalization of a 2-class covariance matrix, it can only be used in binary problems. Therefore, we applied the one-versus-one (OVO) scheme to LMD-CSP for the multiclass problem, so that a k-class problem was transformed into k (k-1)/2 binary class problems. For the signal matrix X_i of the i th experimental data, the covariance matrix is calculated as follows: (6)${\displaystyle R=\frac{X_i{X}_i^{\mathrm{T}}}{\text{trace}\left(X_i{X}_i^{\mathrm{T}}\right)}}$

where “T” represents the transpose operator, and trace(⋅) means to find the trace of a matrix.

For the binary class MI tasks (i.e., shoulder abduction and extension), we calculated the covariance matrix over the trails of each class, and averaged them to obtain the mean covariance matrix R_A and R_E. Then the R_A and R_E were transformed to obtain the spatial filter W_AE (in this paper, W_AE ∈R^{8 ×14}) (Liu et al., 2012). Combined with the OVO scheme, we obtained three spatial filters (W_AE, W_AF and W_EF). Then the three spatial filters were spliced vertically to gain a global spatial filter. The global spatial filter is as follows: (7)${\displaystyle W=\left[W_{\mathrm{AE}};W_{\mathrm{AF}};W_{\mathrm{EF}}\right]}$

where W_AE represents the spatial filter between abduction and extension, W_EF represents the spatial filter between extension and flexion, and W_AF represents the spatial filter between abduction and flexion.

For a single trial data, we obtained the LMD-CSP projection matrix Z = WX_i. The p th row of Z is denoted as Z_p. The required characteristics can be obtained by: (8)${\displaystyle f_p=log\left(\frac{\mathrm{var}\left(Z_p\right)}{\sum _{k=1}^g\mathrm{var}\left({\mathrm{Z}}_k\right)}\right)}$

where p = 1, 2, …g. The features obtained by LMD-CSP are denoted as F = [f ₁, f ₂, …f_g].

Figure 3 shows the distributions of the most significant two LMD-CSP features from subjects 1. Obviously, the features of different MI tasks extracted by LMD-CSP are highly distinguishable and these are easy to separate.

Twin support vector machine

TWSVM, first proposed by Jayadeva, & Khemchandani & Chandra (2007) is a new machine learning method based on traditional SVM. TWSVM aims to construct a hyperplane for each class. It requires that each hyperplane is close to the corresponding class samples as possible, and far away from the other class samples as possible. TWSVM has higher training speed and generalization ability than SVM because the former solves two small quadratic programming problems (QPPs) to construct the hyperplanes like SVM, and the constraint condition of a QPP is only related to one class of samples.

Suppose that training samples of class 1 are denoted as A =[x₁⁽¹⁾, x₂⁽¹⁾, …, x_m1⁽¹⁾]^T ∈R^m1×g, where x_j⁽ⁱ⁾ ∈R^g represents the jth sample of class 1, m1 represents the number of samples. Similarly, training samples of class 2 are denoted as B ∈R^m2×g. When training the TWSVM classifier between class 1 and class 2, two nonparallel hyperplanes are obtained as follows: (9)${\displaystyle K\left(x^T,C^T\right)w_1+b_1=0}$ (10)${\displaystyle K\left(x^T,C^T\right)w_2+b_2=0}$

where C^T =[A^TB^T]^T, w₁ and w₂ are two normal vectors of hyperplanes, b₁ and b₂ are the bias vectors.

TWSVM constructs the two hyperplanes by solving the following optimization problems: (11)${\displaystyle \underset{w^{\left(1\right)},b^{\left(1\right)},{\xi }^{\left(2\right)}}{min}\frac{1}{2}{∥K\left(A,C^T\right)w^{\left(1\right)}+e_1{b}^{\left(1\right)}∥}^2+c_1{e}_2^T{\xi }^{\left(2\right)}}$

$-\left(K\left(B,C^T\right)w^{\left(1\right)}+e_2{b}^{\left(1\right)}\right)\ge e_2-{\xi }^{\left(2\right)}$, ${\xi }^{\left(2\right)}\ge 0$ (12)${\displaystyle \underset{w^{\left(2\right)},b^{\left(2\right)},{\xi }^{\left(1\right)}}{min}\frac{1}{2}{∥K\left(B,C^T\right)w^{\left(2\right)}+e_2{b}^{\left(2\right)}∥}^2+c_2{e}_1^T{\xi }^{\left(1\right)}}$

$-\left(K\left(A,C^T\right)w^{\left(2\right)}+e_1{b}^{\left(2\right)}\right)\ge e_1-{\xi }^{\left(1\right)}$, ${\xi }^{\left(1\right)}\ge 0$

where c ₁ and c ₂ are penalty parameters, e ₁ and e ₂ are column vectors of ones.

In this study, we established a classifier for each binary class, so k (k-1)/2 TWSVM classifiers were constructed. In the process of classification, the feature vectors of each MI task were input into the classifiers, and the final result was obtained by voting. The penalty parameters c ₁ and c ₂ of the TWSVM were set by the following MOGWO.

Multi-objective Grey Wolf Optimizer

MOGWO, proposed by Mirjalili et al. (2016) is a new swarm intelligence optimization algorithm based on the conventional grey wolf optimizer algorithm. When compared with other multi-objective optimization algorithms, MOGWO has higher convergence and coverage (Dilip et al., 2018). In MOGWO, the location of each gray wolf is denoted as a solution, and the first three best solutions are denoted as the alpha (α) wolves, beta (β) wolves, and delta (δ) wolves, the other candidate solutions are omega (ω) wolves. In the iterative process, ω wolves are led by α wolves, β wolves and δ wolves to find global optimal solutions.

MOGWO uses an external archive to store and update the non-dominated Pareto optimal solutions. At the same time, it also employs a leader selection mechanism to search for the least crowded segments from the archive, and three non-dominated solutions of the segments are used as α, β, δ wolves by a roulette-wheel method. The location of each search agent is updated as follows: (13)${\displaystyle \stackrel{\rightharpoonup }{D_i}=|\overrightarrow{C_i}\cdot \overrightarrow{X_i}\left(t\right)-\overrightarrow{X}\left(t\right)|,i\in \left(\alpha ,\beta ,\delta \right)}$ (14)${\displaystyle \overrightarrow{X_i}\left(t\right)=\overrightarrow{X_i}\left(\mathrm{t}\right)-\overrightarrow{A_i}\cdot \overrightarrow{D_i},i\in \left(\alpha ,\beta ,\delta \right)}$ (15)${\displaystyle \overrightarrow{X}\left(t+1\right)=\frac{\overrightarrow{X_{\alpha }}\left(t\right)+\overrightarrow{X_{\beta }}\left(t\right)+\overrightarrow{X_{\delta }}\left(t\right)}{3}}$ where the coefficient vector $\overrightarrow{A}$ and $\overrightarrow{C}$ are calculated as follows: (16)${\displaystyle \overrightarrow{A}=\overrightarrow{2a}\times \overrightarrow{r_1}-\overrightarrow{a}}$ (17)${\displaystyle \overrightarrow{C}=2\cdot \overrightarrow{r_2}}$ where the parameter $\overrightarrow{a}$ is decreased from 2 to 0 in the iterative process, and $\overrightarrow{r_1}$, $\overrightarrow{r_2}$ are random vectors in [0, 1].

Multi-objective function

The average recognition accuracy is usually employed to measure classification performance. However, it only quantifies the overall classification performance and ignores evaluating the classification results of each class. Therefore, we applied the mean recognition accuracy and the recognition of each class as objective functions to evaluate the candidate solutions generated by the MOGWO algorithm. The objective function can be expressed as: (18)${\displaystyle Accuracy=\frac{\sum _{i=1}^{3}N{C}_i}{\sum _{i=1}^{3}N{C}_i+\sum _{i=1}^{3}N{E}_i}}$ (19)${\displaystyle C{R}_i=\frac{N{C}_i}{N{C}_i+N{E}_i},i=1,2,3}$ where Accuracy is the mean recognition accuracy, CR_i represents the recognition accuracy of i th class, NC _i is the number of i th class correctly distinguished samples, NE_i is the number of i th class wrongly distinguished samples.