Feature extraction of four-class motor imagery EEG signals based on functional brain network

LCD is a signal decomposition method that decomposes any complex signal $x(t)(t>0)$ into the sum of $n$ ISC component ${{c}_{i}}\left(t \right)\left(i=1,2,\ldots ,n \right)$ and a residue ${{u}_{n}}(t)$ (see equation (1)). The ISC component must satisfy two conditions: (1) its local waveform is approximately a sine wave, and (2) the ISC of a single mode will not generate a negative frequency. The signal decomposition expression is as follows:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle x\left(t \right)=\underset{i=1}{\overset{n}{\mathop \sum }}\,{{c}_{i}}\left(t \right)+{{u}_{n}}\left(t \right).\nonumber \end{align} \tag{ 1 }$$

First, the baseline of the original signal is calculated using the cubic spline function, and then the baseline is subtracted from the original signal. If the residual signal satisfies two conditions of the ISC component, the signal is an ISC component, otherwise, the signal is taken as the original signal and the above process is repeated. Meanwhile, after each ISC component is obtained, the standard deviation (SD) is calculated according to formula (2), and the iteration is terminated if the SD is less than 0.05.

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {\rm SD}=\underset{t=0}{\overset{T}{\mathop \sum }}\,\left[ \frac{{{\left| {{h}_{ip}}\left(t \right)-{{h}_{i\left(\,p-1 \right)}}\left(t \right) \right|}^{2}}}{h_{i\left(\,p-1 \right)}^{2}\left(t \right)} \right]\nonumber \end{align} \tag{ 2 }$$

where ${{h}_{ip}}\left(t \right)$ represents the ith ISC component obtained after looping p  times. Owing to the computation demand of the entire LCD calculation, only three channels (C3, C4, and Cz) that contribute most to the classification are selected for LCD decomposition [23]. Then, by conducting a lot of experiments and considering the computation time, the signals of these three channels are decomposed by the LCD and only the first three ISC components are taken, so nine ISC components can be obtained from the C3, C4, and Cz channels in one experiment.

Perform Hilbert transform on each ISC component ${{c}_{i}}\left(t \right)\left(i=1,2,\ldots ,K \right)$, formulated as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{y}_{i}}\left(t \right)=\frac{1}{\pi }\int_{-\infty }^{\infty }{\frac{{{c}_{i}}\left(\tau \right)}{t-\tau }}d\tau \nonumber \end{align} \tag{ 3 }$$

where K is the number of ISC components in each experiment, and K is set to nine in our experiment. The parsing signal ${{z}_{i}}\left(t \right)$ is then constructed as follows:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{z}_{i}}\left(t \right)={{c}_{i}}\left(t \right)+j{{y}_{i}}\left(t \right)={{a}_{i}}\left(t \right){{e}^{\,j{{f}_{i}}\left(t \right)}}\nonumber \end{align} \tag{ 4 }$$

where ${{a}_{i}}\left(t \right)$ and ${{f}_{i}}\left(t \right)$ represent the instantaneous amplitude and frequency of the i-th ISC component, respectively. Then, the instantaneous frequency is sorted, and part of the values are selected from the sorted instantaneous frequency ${{f}^{\prime}}(t)$ at medium intervals as the frequency features ${{F}_{1}}=\left[ {{f}_{11}},{{f}_{12}}\ldots ,{{f}_{1K}} \right]\in {{R}^{1\times KP}}$ of the MI EEG signals, where P is the number of features selected from the ISC component, and P is set to 20 through experimental experience. Each of the eigenvalues ${{f}_{11}},{{f}_{12}}\ldots$ in ${{F}_{1}}$ is a 1  ×  P- dimensional vector.

The CSP algorithm uses the theory of matrix simultaneous diagonalization in algebra to find a set of spatial filters in order to maximize the variance of one class of signals while minimizing the variance of the other class of signals.

Denote the original EEG signal of a trial as ${{E}_{N\times T}}$, where N is the number of electrode leads, and T is the number of data samples. Here, we take an example of a two-class experiment, where data are collected from two types of tasks named the left-hand and the right-hand MI tasks. The subjects are instructed to imagine the movement of their left hands or right hands, but without actual muscle activations in their hands.

The CSP feature is calculated with the following steps:

• (1)  
  Calculate the covariance C of the two-class MI signals in each experiment, formulated as

  $$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle C=\frac{E{{E}^{T}}}{trace(E{{E}^{T}})}\nonumber \end{align} \tag{ 5 }$$

  where $trace\left(X \right)$ is the trace of matrix X, which is the sum of the diagonal elements of matrix X. The average covariance of all experiments is then calculated by summing up covariance matrices, in this case, the left-hand and right-hand MI data:

  $$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{C}_{l}}=\sum\limits_{i=1}^{n}{{{C}_{l,i}}}\quad {{C}_{r}}=\sum\limits_{i=1}^{n}{{{C}_{r,i}}}.\nonumber \end{align} \tag{ 6 }$$

  Then, the sum of the two types of covariance matrices is obtained:

  $$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{C}_{c}}={{C}_{l}}+{{C}_{r}}\nonumber \end{align} \tag{ 7 }$$
• (2)  
  Perform an eigenvalue decomposition of the mixed spatial covariance, formulated below:

  $$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{C}_{c}}={{U}_{c}}{{A}_{c}}U_{c}^{T}\nonumber \end{align} \tag{ 8 }$$

  where ${{A}_{c}}$ is the eigenvalue diagonal matrix, and ${{U}_{c}}$ is the corresponding eigenvector matrix.
• (3)  
  Construct the whitening transformation matrix first. Then, using the features of ${{S}_{l}}$, ${{S}_{r}}$ with the same feature vector, decompose its eigenvalue:

  $$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \begin{array}{@{}cccccccccccccccccccc@{}} & \quad \quad P=A_{c}^{-1/2}U_{c}^{T} \nonumber \\ & {{S}_{l}}=P{{C}_{l}}{{P}^{T}}\ \ {{S}_{r}}=P{{C}_{r}}{{P}^{T}} \nonumber \\ \end{array}\nonumber \end{align} \tag{ 9 }$$

  $$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{S}_{l}}=B{{A}_{t}}{{B}^{T}}\ \ {{S}_{r}}=B{{A}_{r}}{{B}^{T}}.\nonumber \end{align} \tag{ 10 }$$

  The desired space filter $W={{\left({{B}^{T}}P \right)}^{T}}$ is obtained, and the filter is used to obtain ${{Z}_{N\times T}}={{W}_{N\times N}}{{E}_{N\times T}}$.
• (4)  
  Find the eigenvector $f$.

The dimension of f  can be adjusted according to the quality of the EEG signals and the classifier requirements, but should not exceed the number of electrode leads N. Extract the first m rows and the last $m$ rows of $Z\left(2m<N \right)$, take the p th row of Z as ${{Z}_{p}}$. Then, ${{f}_{p}}$ is calculated as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{f}_{p}}=\log \left(\frac{{\rm var}({{Z}_{p}})}{\sum\nolimits_{n=1}^{\,p}{\operatorname{var}(Z{{}}_{l})}} \right)\quad p=1:2m\nonumber \end{align} \tag{ 11 }$$

where $\operatorname{var}\left(X \right)$ represents the variance of the time-series signals.

Since there are four classes of MI tasks that need to be classified, it is necessary to expand the CSP to meet the technical requirements. There are two common methods to expand: one to one and one to the other. The one-to-one method is adopted to expand the CSP in this paper. For each experiment, six projection matrices are generated, and each projection matrix is concatenated to form a complete spatial feature vector ${{F}_{2}}$.

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{F}_{2}}=\left[ {{f}_{21}},{{f}_{22}},{{f}_{23}},{{f}_{24}},{{f}_{25}},{{f}_{26}} \right]\in {{R}^{1\times 6{{N}_{1}}}}\nonumber \end{align} \tag{ 12 }$$

where ${{N}_{1}}$ represents the number of channels in each experiment. The signal of one experiment here includes 22 channel data and nine ISC component data, forming a total ${{N}_{1}}$ of 31. Each of the eigenvalues ${{f}_{21}},{{f}_{22}}\ldots$ in ${{F}_{2}}$ is a 1  ×  ${{N}_{1}}$- dimensional vector.

Researchers have found that biological networks generally have properties of small-world networks. This means the network has a large clustering coefficient and a short characteristic path length [24, 25]. Considering the independence of EEG nodes and the synergistic effect between nodes, we can use the complex network theory to construct an EEG functional brain network. In a functional network based on EEG, each node corresponds to the brain regions detected by different leads, and the collected EEG signals constitute the time series of this node.

The definition of edges in functional networks is based on the functional connections, and the weights of the edges can be determined with various methods, which can be broadly categorized into two main branches: linear and nonlinear. The linear methods include the Pearson correlation, partial correlation, and partial coherence. The nonlinear methods include the synchronization likelihood, canonical correlation analysis (CCA), and mutual information.

Through the analysis of various connection methods, we choose the CCA to calculate the nonlinear correlation between each pair of leads. This method can analyze the signal in the entire EEG frequency band as well as in a specific range of frequency spectra. This is ideal for analyzing instantaneous and unstable signals such as EEGs. CCA considers the linear combination of the two sets of variables and studies the correlation coefficient $\rho (u,v)$ between them. In all linear combinations, we find a linear combination with the largest correlation coefficient and use this maximum correlation coefficient to represent the correlation of the pair of variables. The correlation coefficient is expressed as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{\rho }_{uv}}=\frac{{\rm Cov}(U,V)}{\sqrt{{\rm Var}(U)}\sqrt{{\rm Var}(V)}}.\nonumber \end{align} \tag{ 13 }$$

Mathematically, the functional network obtained is a correlation matrix in which each element represents a correlation between two brain regions. After obtaining the correlation matrix, the next step is to binarize it by setting a threshold ${{C}_{thr}}$. If the value of an element in the matrix is greater than this threshold, it is considered that there is a functional connection between the two brain regions, and the value here is set to 1. Otherwise the value is set to 0, thereby establishing a complete binarized function network. The process of constructing functional brain networks based on EEG signals is shown in figure 1.

Zoom In Zoom Out Reset image size

In order to describe the topological structure of the established brain network, some common measures are provided, and the changes in the topological properties of the brain networks can be studied by analyzing these measures. We select five measures: degree, clustering coefficient, average shortest path length, local efficiency, and betweenness centrality. These complex network measures are described in detail below.

2.1.3.1. Degree.

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle k_{i}^{B}=\sum\limits_{j\in G}{{{a}_{ij}}}\quad {\rm or}\quad k_{i}^{W}=\sum\limits_{j\in G}{{{w}_{ij}}}\nonumber \end{align} \tag{ 14 }$$

where ${{a}_{ij}}$ or ${{w}_{ij}}$ is the corresponding element in the binary or weighted network matrix. The degree of a node is the number of edges of one node. A node with a higher degree is considered more important in the network.

2.1.3.2. Clustering coefficient.

The clustering coefficients of the nodes are defined as follows [26]:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \begin{array}{@{}llcccccccccccccccccc@{}} & c_{i}^{B}=\frac{2R}{k_{i}^{B}(k_{i}^{B}-1)} \nonumber \\ & \quad \quad \quad \quad \quad {\rm or} \nonumber \\ & c_{i}^{W}=\frac{2}{k_{i}^{B}(k_{i}^{B}-1)}{{\sum\limits_{j,k}{({{w}_{ij}}{{w}_{jk}}{{w}_{ki}})}}^{1/3}} \nonumber \\ \end{array}\nonumber \end{align} \tag{ 15 }$$

where $R$ represents the number of directly connected neighbors of node i. The clustering coefficient of a node is defined as the ratio between the actual number of edges existing between the neighbor nodes of the node and the maximum possible number of connected edges. The clustering coefficient reflects the local connectivity and measures the cluster characteristics and closeness within the functional brain network.

2.1.3.3. Average shortest path length.

The shortest path length is the smallest number of edges between two nodes. In other words, it is the minimum number of steps to travel through the network from node i to j . The average shortest path length is defined as the mean number of steps along the shortest paths between all possible pairs of network nodes. The definition is as follows:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle L=\frac{1}{N(N-1)}\sum\limits_{i,j\in V,i\ne j}{{{d}_{ij}}}\nonumber \end{align} \tag{ 16 }$$

where N represents the total number of network nodes. (The number of network nodes is the same as the number of electrode leads in our experiment), d_ij represents the distance between nodes i and j  in the network.

2.1.3.4. Local efficiency.

For a network G with N nodes, the global efficiency is calculated as shown in equation (17) [27]:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{E}_{glob}}(G)=\frac{1}{N(N-1)}\sum\limits_{i\ne j\in G}{\frac{1}{{{d}_{ij}}}}.\nonumber \end{align} \tag{ 17 }$$

The formula for calculating the local efficiency is as follows:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{E}_{loc}}(G)=\frac{1}{N}\sum\limits_{i\in G}{{{E}_{glob}}}({{G}_{i}})\nonumber \end{align} \tag{ 18 }$$

where ${{E}_{glob}}\left({{G}_{i}} \right)$ is the global efficiency of G_i, and G_i is a subgraph composed of the neighbors of node $i$. The global efficiency and local efficiency measure the information transmission ability of the network globally and locally.

2.1.3.5. Betweenness centrality.

The betweenness centrality is defined as the number of shortest paths going through a node or edge [28]. The higher the betweenness centrality of a node, the greater the flow of information carried by the node, and the more significant the impact on the function of the functional brain network. The betweenness centrality is formulated in equation (19):

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{C}_{B}}(i)=\frac{2\sum\limits_{h\ne i,j\ne i,h<j\in V}{{{g}_{hj}}(i)}}{(N-1)(N-2){{g}_{hj}}}\nonumber \end{align} \tag{ 19 }$$

where ${{g}_{hj}}$ is the number of all shortest paths between node $h\in V$ and $j\in V$, and V is the set of all nodes in the network.

After quantifying the relationship between nodes, it is necessary to select an appropriate threshold to binarize the adjacency matrix. Two principles need to be followed to establish the network: To ensure the integrity of the network, it should not contain any isolated node or isolated part; and the small-world characteristics of the network should be ensured [24, 25]. According to the random model of Erdos and Renyi [29], if a graph with N nodes is to be fully connected, the connection sparsity should be greater than $2{\rm ln}\,N/N$. In addition, it should be ensured that its small-world attribute value $\sigma$ is much greater than 1. Through a large number of experiments, the threshold was empirically set to 0.84, and the corresponding sparsity of the brain network was 0.35.

The offline experimental results of chapter 4 below show that the classification effect of the measure of degree in the binary network is better than for other measures. Thus, the measure of degree is used as the feature of the brain network for online experiments.

After establishing a functional brain network based on MI EEG signals, the measures described in the previous section are extracted as the features of brain networks and are then fused with the multiscale features extracted from the CSP and LCD algorithms. There are two feature fusion strategies: parallel feature fusion and serial feature fusion. Compared to parallel feature fusion, one advantage of serial feature fusion is its simplicity in that it requires two steps: normalization and concatenation of multiple features. This effectively retains the discriminative information of various features for classification. For the above reason, the serial feature fusion strategy is adopted in this work.

With the features in the above spatial and frequency domains and the features of the brain network fused, the obtained feature vector of the EEG signals, denoted by $F\in {{R}^{1\times \left(KP+6{{N}_{1}}+MN \right)}}$, is defined below:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \left\{\begin{array}{@{}lllllccccccccccccccc@{}} F=\left[ {{F}_{11}},{{F}_{22}},{{F}_{33}} \right] \nonumber \\ {{F}_{11}}=\left[ \frac{{{f}_{11}}}{\left\Vert {{f}_{11}} \right\Vert},\frac{{{f}_{12}}}{\left\Vert {{f}_{12}} \right\Vert},\ldots ,\frac{{{f}_{1K}}}{\left\Vert {{f}_{1K}} \right\Vert} \right]\in {{R}^{1\times KP}} \nonumber \\ {{F}_{22}}=\left[ \frac{{{f}_{21}}}{\left\Vert {{f}_{21}} \right\Vert},\frac{{{f}_{22}}}{\left\Vert {{f}_{22}} \right\Vert},\ldots ,\frac{{{f}_{26}}}{\left\Vert {{f}_{26}} \right\Vert} \right]\in {{R}^{1\times 6{{N}_{1}}}} \nonumber \\ {{F}_{33}}=\left[ \frac{{{f}_{31}}}{\left\Vert {{f}_{31}} \right\Vert},\frac{{{f}_{32}}}{\left\Vert {{f}_{32}} \right\Vert},\ldots ,\frac{{{f}_{3M}}}{\left\Vert {{f}_{3M}} \right\Vert} \right]\in {{R}^{1\times MN}}. \nonumber \\ \end{array} \right.\nonumber \end{align} \tag{ 20 }$$

The three channels C3, C4, and Cz that contribute most to the classification are selected to perform LCD decomposition [23], which produces nine ISC components and then extracts the frequency domain feature F₁. Next, the obtained nine ISC components are added to the 22 channels of original EEG signals, and the CSP algorithm is used to extract the spatial features from the 31-channel data as ${{F}_{2}}$. The feature vector ${{F}_{3}}$ is extracted from the functional brain network, where $\left\Vert \circ \right\Vert$ is the l2-norm in equation (20). A flowchart of the feature extraction algorithm is shown in figure 2.

Zoom In Zoom Out Reset image size

The EEG data were recorded using a UE-16B EEG amplifier at a sampling rate of 1000 Hz. All 16 channels were selected (Fp1, Fp2, F3, F4, C3, C4, P3, P4, O1, O2, F7, F8, T3, T4, T5, and T6). Since the frequency of the MI EEG signal is concentrated at 8–30 Hz, and the cutoff frequency of the low-pass filter was set to 100 Hz (with a 50 Hz notch filter enabled). The left and right ear electrodes, A1 and A2, were used as reference electrodes, and the forehead was used as ground electrode G. The controlled object in the online BCI system was a humanoid robot, the NAO robot, produced by Aldebaran Robotics.

The subjects were eight graduate students aged from 23 to 26, including two females and six males. All of them were right-handed and had no neurological history. To minimize environmental effects such as light stimulations, the experiments were carried out in a quiet environment with dimmed lighting. The subjects were seated in comfortable backrest chairs to reduce muscle strain that could interfere with the experimental results. Before the experiments, all subjects were instructed and trained about the experimental procedure. The study was conducted with approval from the Wuhan University of Technology.

The MI EEG data acquisition included training data collection and real-time data collection. The training data collection session consisted of four runs with 2 min breaks between two consecutive runs. Each run had 25 MI trials. The MI tasks to be performed in the four runs were left hand, right hand, both foot, and tongue movements. Each trial consisted of a 2 s ready period, a 4 s MI period, and a 2 s break period. When the prompt picture was displayed on the screen, the subject needed to perform the corresponding MI task until the prompt image disappeared. The subject then rested for 2 s and waited for the next experiment to begin. The experimental process for training EEG data acquisition is shown in figure 3.

Zoom In Zoom Out Reset image size

We segmented the data in this duration into several epochs to conduct a series of experiments statistically, and found that the data between 2.5 s and 3.5 s achieved the best performance. Therefore, this 1 s data epoch was selected for feature extraction and classification.

The real-time data acquisition session of the robot control continuously collected data, and the data of the 1 s time period was collected and sent to the signal processing module, where the data of the next second was collected at the same time.

This subsection introduces the BCI system framework designed for evaluating the proposed algorithm. The system includes four functional components: signal acquisition, signal processing, human–computer interaction, and robot control. The signal acquisition module is responsible for the data acquisition, filtering, and amplification of EEG signals, which are then sent to the human–computer interaction module in real time. Next, the human–computer interaction module stores the received EEG data, and the signal processing module starts to process the data that are eventually converted into a set of control commands. The commands are sent to the robot control module through socket communication. The NAO robot interprets the commands and starts executing the corresponding movement behaviors according to the received commands. A block diagram of the overall system is shown in figure 4.

Zoom In Zoom Out Reset image size

The mapping between the MI tasks and the corresponding movement control commands of the robot are listed in table 1 together with their labels used for classification. As can be seen, the robot has four basic behaviors moving in four directions: forward, left, right, and backward. Correspondingly, the subjects are asked to execute the mental activities of moving four body parts, both feet, left hand, right hand, and tongue.

In terms of interacting with the robot, two control strategies were adopted for the motion control of the NAO robot: synchronous control and asynchronous control.

• (1)  
  Synchronous controlThe principle of the NAO robot synchronously controlled by the MI signals is as follows. The robot needs to perform a predefined sequence of motion behaviors. A motion behavior corresponds to a control instruction. In this work, we test with a total of 17 instructions. For example, [2 2 0 0 2 2 0 0 1 0 0 1 1 3 3 3 1] is the sequence of labels (as seen in table 1) of the control commands for testing in this work. The subject then performs the corresponding MI tasks with the given instructions displayed in order. During the experiment, the robot will carry out the corresponding motion only if the classification result of the MI signal matches the preset command. The experiment is completed when the robot completes all instructions.
• (2)  
  Asynchronous controlAsynchronous control differs in that it does not require a pre-known sequence of control commands. The NAO robot starts performing behaviors purely based on the classification results. The task is to control the robot to move toward a predefined goal position in the room (figure 9). The decisions of robot motions are not constrained and are only based on the classification results of the subject's mental activities. However, false classifications are unavoidable. The expected motions from the subject's mental activities cannot always be correctly recognized, resulting in incorrect motion behaviors executed by the NAO robot.

To avoid or minimize the number of false alarms, an error control mechanism is deployed in the experiment of the asynchronous control robot. It is assumed that the subject's mental activities will generally remain temporally consistent, meaning the subject will not change his or her task very quickly. Therefore, the false alarm detector used in the work compares the classification results of two consecutive detections by comparing the current classification result with the previous one. If the two classification results are the same, it is determined that the movement is correct, and the corresponding control instruction is sent to the robot, the sent instruction is also compared with the next classification result. However, if the two are inconsistent, a false alarm will be triggered, and the control command is not sent until the sending condition is met again. This control mechanism greatly increases the classification success rate of MI tasks that the subject needs to perform, and it effectively reduces the chances of performing incorrect movements by the NAO robot.