# A Fusion Algorithm for EEG Signal Processing Based on Motor Imagery Brain-Computer Interface

Authors: Xiaozhong Geng; Song Xue; Ping Yu; Dezhi Li; Mengzhe Yue; Xi Zhang; Linen Wang
Journal: Wireless Communications and Mobile Computing
 (2022)

---

## Abstract

Electroencephalogram (EEG) signal processing is a very important module in the brain-computer interface system. As an important physiological feature of the human body, EEG signals are closely related to the functional state of the cerebral nervous system. However, the EEG signals collected on the scalp are generally weak and inevitably subject to various noise interferences. In order to remove artifacts from the EEG in brain computer interfaces (BCIs), a fusion algorithm for EEG signal preprocessing is proposed. The fusion method includes the following steps: firstly, the raw EEG signals are separated into a set of statistics independent components (ICs) by the improved FastICA algorithm. Then, each independent component is decomposed into a series of intrinsic mode functions (IMFs) by using the improved empirical mode decomposition method (EMD). Many IMFs with high-frequency noise are deleted. The rest of the IMFs are reconstructed. Furthermore, artifacts are further eliminated by iterative process of the improved FastICA algorithm, and then, the EEG signals are reconstructed again by inverse ICA. Finally, the cleaned EEG signal was obtained. The comparative experiment shows that the EMD-ICA fusion algorithm not only accurately eliminates the artifact components but also better retains the local characteristics of the raw EEG. Continuous wavelet transform was used to extract energy features ofμ rhythm and β rhythmic to represent the characteristics of EEG signals under different motor imageries. These two features are normalized and used as the input data of the convolutional neural network (CNN) designed by the paper, and the two kinds of features are learned by CNN, and then, the two-classification problem of motor imagery EEG signals is completed. The experimental results show that the average classification accuracy and kappa value of the proposed method are higher than those of SVM and SAE for most subjects.


---

## Body

## 1. Introduction

Brain-computer interface (BCI) is a communication control system established between the brain and external devices (computers or other electronic devices) through EEG signals generated during brain activity [1]. This system realizes the direct connection between the brain and the external world. With the help of this system, people can directly manipulate some devices such as computers, wheelchairs, artificial limbs, etc. through motor imagery of the brain [2]. This can effectively enhance the ability of patients with severe physical disabilities to communicate with the outside world or control the external environment, so as to improve the quality of life of the patients. Therefore, brain-computer interfaces (BCIs) can be used for the treatment of patients with severe neuromuscular diseases, as well as for the rehabilitation of patients with stroke, head injury, or other diseases [3].Usually, the BCI system is composed of three modules: signal acquisition module, signal processing module, and equipment control module, as shown in Figure1 [4]. The signal processing module mainly involves preprocessing, feature extraction, and feature classification. The signal acquisition module mainly completes recording of EEG signals and transmits the signal to the signal processing module; the signal processing module mainly completes the preprocessing of the collected brain signals, then extracts the features of the signal that reflect the subject’s intention, and at last classifies the extracted features and converts them into control commands; the device control module realizes the control of external devices according to the obtained control commands. The EEG signal processing module is the most critical part of the brain-computer interface system, and it is also a current research hotspot [5].Figure 1
Structure of brain-computer interface system.Motor imagery (MI) EEG signal (EEG) refers to the EEG signal generated when the imagination of the movement of the limbs or other parts is repeatedly simulated and rehearsed in the brain without the actual participation of the human body in the movement behavior [6]. At present, in the research of brain-computer interface systems that use motor imagery EEG as the signal to be classified, the requirements for classification performance are always high. Therefore, it is necessary to study how to improve the classification accuracy of EEG signals and reduce the classification error.Motor imagery potentials (motor imaginary potentials) refer to the changes in scalp potentials caused by neural electrical activity in the brain motor cortex area related to the action when the brain imagines a certain limb action. Xu et al. first discovered imagery action potentials [7]. In the process of studying electroencephalograph (EEG) signals, they noticed that the imagination of limb movement can cause changes in the activity state of a large number of nerve cells in the cortical motor center. It causes the synchronization of certain frequency components in the EEG to increase or decrease, which is the so-called event-related synchronization (ERS) and event-related desynchronization (ERD) phenomenon. ERD corresponds to the decrease of power spectrum intensity, and ERS corresponds to the increase of power spectrum intensity [8]. The ERD/ERS caused by actions of different parts are different in the frequency band and cortical area [9]. Zheng and Guo confirmed the above phenomenon through experiments, pointed out that the ERD/ERS phenomenon is mainly concentrated in the mu rhythm and the beta rhythm segment in the EEG, and proposed a quantitative theory for the ERD/ERS phenomenon [10]. Figure 2 is the distribution of ERD/ER∗S phenomenon in the cerebral cortex for the left- and right-hand motor imagery.Figure 2
The distribution of ERD/ERS phenomenon in the cerebral cortex for the left- and right-hand motor imagery.However, with the increase of motor imaging tasks, the recognition accuracy rate will be greatly deviated. This is because the scalp potential is very susceptible to the interference of artifact signals. The collected original EEG signals are often mixed with electrooculogram (EOG) artifacts, electromyogram (EMG) artifacts, electrocardiogram (ECG) artifacts, and power frequency interference. The existence of artifacts makes EEG analysis and processing a certain degree of difficulty [11].In the process of collecting EEG signals, the location of the sampling electrode is close to the eye, so the EOG artifact has the greatest interference to the EEG, and it often appears randomly in the EEG signal with higher amplitude. A large number of studies have shown that the effect of using traditional frequency domain filtering methods to remove the artifact spectrum mixed in the EEG signal is not obvious.Independent component analysis (ICA) is a multidimensional signal processing method, which belongs to a new blind source signal separation method [12]. The basic principle of ICA is to separate the observed mixed signal into a combination of statistically independent non-Gaussian signal sources through linear transformation; that is, the independent source signal can be analyzed from the observed multichannel mixed signal. It is generally believed that the EEG signal collected by electrodes covering the entire scalp is a linear combination of signals from different active areas of the brain. Various artifacts are propagated through the scalp and mixed with EEG signals. This mixing is also considered to be a linear mixture [13]. Therefore, the independent component analysis method can separate the independent EEG components and artifact interference from the collected signals.Empirical mode decomposition (EMD) was proposed by Chinese American N.E. Huang et al. in 1998 [14]. Empirical mode decomposition (EMD) is an adaptive signal processing technology. It is a time-frequency analysis method for processing nonlinear and nonstationary signals. This method can adaptively decompose the signal into a limited number of intrinsic mode function (IMF) components and a margin that characterizes the signal trend based on the characteristics of the input signal itself without knowing any prior knowledge. This method has unique advantages in nonstationary signal processing, so it has been widely used in EEG signal processing.Grossman and Morlet proposed the Wavelet Transform (WT) in 1984 [15]. The wavelet transform is a new development of the Fourier transform, and the wavelet transform coefficients can reflect the local information of the signal in the time and frequency domains. The main advantage of the wavelet transform is that it has a variable time-frequency analysis window. A wide window can be used to analyze low-frequency signals, and a narrow window can be used to analyze high-frequency signals. In this way, wavelet transform can provide optimal time-frequency resolution for signal analysis in all frequency ranges. Moreover, since the wavelet transform window range can automatically adapt to the instantaneous events of each scale, it is especially suitable for analyzing nonstationary signals, such as EEG signals.Deep learning (DL) is a new research direction in the field of machine learning (ML) [16]. Deep learning originated from artificial neural networks, which simulate the nerve cells in the human brain by using a large number of nodes (neurons). Each neuron has a specific activation function, which uses different weights to represent the synapses that connect two neurons. A large number of neurons with activation functions are connected in different ways to simulate the local nervous system of the human brain. Deep learning technology has powerful signal processing and recognition capabilities, and it has been gradually applied to the field of EEG signal analysis to meet the needs of automatic EEG signal detection in clinical and scientific research. The preprocessed EEG signals are input into convolutional neural networks, cyclic neural networks, deep confidence networks, etc., through the neural network to complete the automatic extraction of EEG features, combined with timing information to complete some end-to-end operations such as detection, recognition, and behavior prediction of EEG signals.In summary, due to the characteristics of the EEG signal itself, it is difficult to achieve an ideal denoising effect simply by using traditional denoising methods. Therefore, this paper proposed the method of combined EMD and ICA, which can effectively remove the influence of noise artifacts.The rest of the paper is organized as follows. Section2 outlines the data sets. Section 3 states algorithm model and the steps of the proposed method. Section 4 presents the experimental results and analysis. Finally, the paper is concluded in Section 5.


## 2. Data Sets

BCI Competition 2008–Graz data set A was used in this paper. It was provided by the Institute for Knowledge Discovery (Laboratory of Brain-Computer Interfaces), Graz University of Technology, (Clemens Brunner, Robert Leeb, Gernot Müller-Putz, Alois Schlögl, and Gert Pfurtscheller). This data set consists of EEG data from 9 subjects. The cue-based BCI paradigm consisted of four different motor imagery tasks, namely, the imagination of movement of the left hand (class 1), right hand (class 2), both feet (class 3), and tongue (class 4). Two sessions on different days were recorded for each subject. Each session is comprised of 6 runs separated by short breaks. One run consists of 48 trials (12 for each of the four possible classes), yielding a total of 288 trials per session [17].At the beginning of each session, a recording of approximately 5 minutes was performed to estimate the EOG influence. The recording was divided into 3 blocks: (1) two minutes with eyes open (looking at a fixation cross on the screen), (2) one minute with eyes closed, and (3) one minute with eye movements.The subjects were sitting in a comfortable armchair in front of a computer screen. At the beginning of a trial (t=0 s), a fixation cross appeared on the black screen. In addition, a short acoustic warning tone was presented. After two seconds (t=2 s), a cue in the form of an arrow pointing either to the left, right, down, or up (corresponding to one of the four classes, left hand, right hand, foot, or tongue) appeared and stayed on the screen for 1.25 s. This prompted the subjects to perform the desired motor imagery task. No feedback was provided. The subjects were asked to carry out the motor imagery task until the fixation cross disappeared from the screen at t=6 s. A short break followed where the screen was black again. The paradigm is illustrated in Figure 3 [18].Figure 3
Timing scheme of the paradigm.Twenty-two Ag/AgCl electrodes (with interelectrode distances of 3.5 cm) are used to record the EEG; the montage is shown in Figure4 [18]. All signals were recorded monopolarly with the left mastoid serving as reference and the right mastoid as ground. The signals were sampled with 250 Hz and bandpass filtered between 0.5 Hz and 100 Hz. The sensitivity of the amplifier was set to 100 μV. An additional 50 Hz notch filter was enabled to suppress line noise.Figure 4
Electrode montage corresponding to the international 10-20 system.In addition to the 22 EEG channels, 3 monopolar EOG channels are recorded and also sampled with 250 Hz. They were bandpass filtered between 0.5 Hz and 100 Hz (with the 50 Hz notch filter enabled), and the sensitivity of the amplifier was set to 1 mV. The EOG channels are provided for the subsequent application of artifact processing methods [19] and must not be used for classification. A visual inspection of all data sets was carried out by an expert, and trials containing artifacts are marked. Eight out of the total of nine data sets were analyzed in [20, 21].


## 3. Algorithm Model

### 3.1. EMD

The key of EMD is to decompose a complex signal into a finite number of intrinsic mode function (IMF), and each of the decomposed IMF components contains local characteristic signals of different time scales of the raw signal.IMF has two limiting conditions [22]:
(1)
In the entire data sequence, the number of extreme points Ne and the number of zero crossing points Nz are equal or differ by at most 1, namely,(1)Ne‐1≤Nz≤Ne+1(2)
At any time pointt, the mean value of the upper envelope curve Mmaxt formed by the local maximum point and the lower envelope curve Nmint formed by the local minimum point is always zero. This is(2)Mmaxt+Nmint2=0Since empirical mode decomposition produces end effects in the decomposition process, which causes distortion problems near the two end points of the signal, a method to suppress the end effects is proposed. The decomposition steps of improved EMD are as follows:(1)
Take the initial end point of the raw signalXt as a maximum starting point, and mark it as M0, select an adjacent maximum point, and mark it as M1, connect the two points, and calculate its slope S1(3)S1=M1−M0tm1‐t0(2)
Calculate the minimum value obtained by the extension at the starting end point(4)N0=N1−S1•tn1−t0,in which, N1 refers to the minimum value closest to the starting point, and tn1 refers to the abscissa of N1. The abscissa of the minimum value is Tm0=t0(3)
Correspondingly, take the terminal end of the raw signal as the minimum point, and mark it asNQmin, select an adjacent minimum point, and mark it as NQmin−1, and calculate the slope S2 of the connection between the two points(5)S2=NQmin−1−NQmintnQmin−1−tQmin(4)
Calculate the maximum value obtained by terminal extension(6)MQmax=MQmax−1−S2•tmQmax−1−tNIts abscissa isTmQmax=tN.
(5)
Find the maximum value MAX and the minimum value MIN in the maximum value sequence. IfMQmax<MIN, define the maximum extension value as M1Qmax=terminalendpointvalue+max/2; if Ν0>max, define the maximum extension value as N10=startingendpointvalue+min/2. This method avoids the situation that the extended extreme value does not meet the interpolation requirements(6)
Find all the maximum and minimum points in the raw data signalXt, and use cubic spline interpolation function to fit the upper envelope Mmaxt and lower envelope Nmint.(7)
Calculate the local average value of the raw signalXt according to the upper and lower envelopes and record it as(7)m1t=Mmaxt+Nmint2(8)
The difference between the raw signal and the local mean is recorded as(8)h1t=Xt−m1t(9)
Ifh1t does not meet the two conditions of IMF, take h1t as the raw signal and go to (7) until the condition is met. At this time, h1t is the first IMF component. And then, obtain the remaining signal r1t by subtracting h1t from the raw signal Xt(9)r1t=Xt−h1t(10)
Regardingr1t as a new signal and repeating the above steps, n IMF components can be screened out(10)rnt=rn−1t−hntWhenrnt is a monotonic sequence, the screening ends.
(11)
Remove the high-noise components both IMF0 and IMF1(12)
The raw data sequence can be reconstructed as(11)X~t=∑i=2nhit+rntX~t is the reconstructed signal after EMD denoising.Extend the minimum and maximum points at both ends of the original signal, which not only improves the EMD method to better suppress the end effect of empirical mode decomposition but also reduces the intrinsic mode function components and postprocessing time.Figure5 is a schematic diagram of the EMD decomposition of the EEG data of the first subject and the instantaneous frequency of the corresponding IMF. It can be seen that the EMD method has good adaptability to the decomposition of EEG signals. After the raw signal is decomposed by EMD, 9 IMF components are obtained. Each order of IMF represents the characteristic signals of different time scales in the raw EEG signal. And the obtained IMF components are arranged in a descending order of frequency. Obviously, IMF0 and IMF1 are high-frequency noise. Therefore, when the EMD method is used to denoise the EEG signal, the denoising effect can be achieved by deleting the IMF component of the high-frequency part, and finally, the remaining part of the IMF is superimposed to obtain the denoised EEG signal.Figure 5
EMD decomposition results of the first subject’s EEG data.Figure6 is a comparison diagram of EEG signals before and after EMD denoising. It can be seen from Figure 6 that the EMD method removes the high-frequency noise in the raw signal, but it is still an unsmooth signal, which still contains more noise signals. The reason for this phenomenon is that the low-frequency part is not processed during the EMD denoising process, and the interference of the ECG and EMG has not been eliminated. Therefore, the EMD denoised EEG signal and the raw signal in the time domain do not seem to change much. The EMD method alone cannot effectively remove the high-frequency noise contained in the EEG signal.Figure 6
Comparison of EMD denoising results before and after.


### 3.2. ICA (Independent Component Analysis)

Assume thatSt=s1t,⋯,sntT is an N-dimensional vector composed of the raw signal source [23], where sit is its component, i=1,2,⋯,n. sit is linearly combined by the hybrid system A into an M-dimensional observation vector Xt=x1t,⋯,xmtT. That is,
(12)xit=ai1s1t+ai2s2t+⋯+ains1t=∑j=1naijsj,i=1,2,⋯,m.Formula (12) may be expressed as follows:
(13)Xt=A×St,in which A is an m×n mixing coefficient matrix. The task of ICA is to find an unmixing coefficient matrix W when the coefficient matrix A and the source matrix S are unknown and then separate the source signals St from the observed signal Xt, that is,
(14)Yt=WXt.in which, Yt=y1t,⋯,ymtT is the output after unmixing, W is a n∗m matrix. It is required that Yt be as close to St as possible. The mathematical model of ICA is shown in Figure 7.Figure 7
The mathematical model of ICA.Because each independent source signal and the mixing matrix are unknown, the output of the ICA has uncertainty, including the uncertainty of the amplitude (or variance) and the uncertainty of the output sequence.The independence criteria of ICA mainly include SOBI (Second-Order Blind Identification), FastICA (Hyvärinen’s Fixed Point Algorithm), information maximization (infomax), and JADE (Joint Approximation Diagonalization of Eigenmatrices) [24].This paper uses an improved FastICA algorithm. Specific steps are as follows.(1)
Suppose the raw EEG observation signalXt=x1t,⋯,xmtT, where xit represents the ith component of the raw EEG observation signal; m represents the number of components of the raw EEG observation signal; and t represents the sampling time of the raw EEG signal(2)
Demeaning for the raw EEG observation signalXt to obtain signal X’t:(15)X′t=Xt−EXt(3)
Use the following formula (16) to decompose the processed EEG observation signal X′t into signals Zt=z1t,⋯,zmtT that are not correlated with each other:(16)Zt=B~X′t,(17)B~=Ds−σ2I‐1/2UsT,(18)σ≅EX′t∗X′tTI,in which B~ represents the projection factor; Ds is a diagonal matrix whose diagonal elements are eigenvalues of the covariance matrix CX=EX′t∗X′tT of X’t; and Us is a matrix with the unit norm eigenvector of CX as a column
(4)
Set the initial value of the unmixing matrix toW(5)
The optimal valueW∗ of the unmixing matrix is as follows:(19)W∗=EZtGWTZt−EG′WTZtW(6)
The optimized valueW∗ of the unmixing matrix is subjected to decorrelation and normalization processing, and the processed unmixing matrix W∗¯ is obtained:(20)W∗¯=W∗W∗(7)
Determine whether the processed unmixing matrixW∗¯ is convergent. If it is convergent, the matrix W∗¯ is the final unmixing matrix. If not, W=W∗¯, go to (5)(8)
Use formula (21) to process the EEG signal Zt=z1t,⋯,zmtT obtained in step (3):(21)Yt=W∗¯∗Zt,in which Yt=y1t,⋯,ymtT and y1t,⋯,ymt are not related to each otherAn object is selected from the Graz data set A as the input signal. Figure8 shows the raw EEG signal before FastICA containing the artifacts such as EOG and ECG. It can be clearly seen that the artifacts seriously interfere with the EEG signal. Figure 9 shows the clean EEG signal by FastICA. Comparing the EEG signal of Figure 8 with Figure 9, it is easy to draw the conclusion that the artifacts such as EOG and ECG are effectively removed by FastICA.Figure 8
Raw EEGs before FastICA.Figure 9
EEGs after FastICA.


### 3.3. EMD-ICA Method

In this paper, the EMD-ICA joint algorithm is used to denoise the EEG signal. The algorithm steps are as follows:(1)
The raw EEG signal is decomposed inton independent components X1,X2,⋯,Xnusing the improved FastICA algorithm in Section 3.2(2)
Decompose each independent componentXi into multiple intrinsic mode functions (IMFs) and margins using the improved empirical mode decomposition method (EMD) in Section 3.1(3)
Delete IMF0 and IMF1 components with higher noise(4)
Superimpose and reconstruct the decomposed IMF component and margin to obtainY(5)
Use the inverse improved FastICA algorithm in Section3.2 to reconstruct signal Y to obtain cleaned EEG signal ZThe EMD-ICA joint algorithm is used to preprocess the EEG signal. The joint algorithm eliminates the interference of noise to the maximum on a finer scale, while retaining most of the useful EEG signals, so as to better complete the purpose of denoising for EEG signals. The denoising process of EEG based on the above method is shown in Figure10.Figure 10
EMD-ICA fusion denoising process.The raw data of the first subject is selected for EMD-ICA decomposition, and the experimental results are shown in Figure11. It can be seen that the raw EEG signal is not smooth and is mixed with noise signals with higher amplitude. After EMD-ICA denoising the raw data, the obtained signal is more in line with the natural trend of the original data. This indicates that the EMD-ICA fusion algorithm not only accurately eliminates the artifact components but also better retains the local characteristics of the raw EEG.Figure 11
Experimental results of EMD-ICA fusion method.


### 3.4. Continue Wavelet Transform

Continuous wavelet transform is to project a one-dimensional signal onto a two-dimensional time-scale plane, which is established by translation and scaling of a mother wavelet. The definition of the wavelet function is as follows [25, 26]:
(22)∫‐∞+∞ϕtdt=0, in which, ϕt is the mother wavelet function, ϕt=1. The wavelet function generated by ϕt can be expressed as
(23)ϕa,bt=1bϕt−ab,in which a and b are translation and scale parameters, respectively.For any signal functionXt, decompose Xt on this wavelet function, and its wavelet transform can be obtained as
(24)WXa,b=X,ϕa,b=∫−∞+∞Xt1bϕt−ab¯dt,in which W is the result of wavelet transform; X,ϕa,b represents the inner product of X and ϕa,b.The formula for the reconstruction of the signalXt is
(25)Xt=1λ∫0+∞∫−∞+∞WXa,b1bϕt−abdadbb2,where
(26)λ=∫0+∞ϕ∧ω2+ωdω,in which ω represents the analog angular frequency. ϕ^ω is the Fourier transform of the mother wavelet function ϕx.The diversification of the mother wavelet function is one of the main advantages of the wavelet transform. Commonly used mother wavelet functions include Haar wavelet, Morlet wavelet, Mexican Hat (Mexh) wavelet, and Daubechies (DbN) wavelet [27].Left-hand and right-hand motor imagery tasks led to ERD and ERS phenomena on the left and right sides of the cerebral motor cortex, respectively, which affected the EEG signals of the C4 and C3 electrodes. Figures12 and 13 are wavelet time-frequency maps of the C3 and C4 channels during the left-hand motor imagery task.Figure 12
Wavelet time-frequency map of C3.Figure 13
Wavelet time-frequency map of C4.It can be seen from Figures12 and 13 that the energy of the C4 channel dropped significantly at the 3rd second after the start of the experiment and recovered after a period of time; that is, the ERD phenomenon occurred. However, the energy of the C3 channel remains high rather than decreasing, which is known as the ERS phenomenon. The maps of C3 and C4 are combined into the final input image, as shown in Figure 14.Figure 14
The combined input map of C3 and C4.


### 3.5. Convolutional Neural Network

As a multilayer neural network, convolutional neural network (CNN) has been successfully applied in the fields of computer vision and image processing due to its powerful feature extraction ability. CNN usually consist of multiple convolutional and pooling layers. A typical CNN has the following layers: input layer, convolutional layer, pooling layer, fully connected layer, and output layer. In the convolution layer, the convolution filter is used to perform convolution filtering on the input data, which has the characteristics of local linking and weight sharing. Compared with the fully connected layer, the convolutional layer can extract the local information of the input data and reduce the number of parameters. In the pooling layer, the method of maximum pooling or average pooling is usually used to reduce the dimension of features, which reduces the complexity of network computation [28].Compared with traditional algorithms, the ability of deep learning algorithm to extract features is greatly improved, and usually, the more complex the network, the more sufficient the features extracted, and the better the obtained classifier effect. However, the advantages of deep learning algorithm classification accuracy are usually only reflected when the number of sample sets is large enough, and the more complex the network, the more parameters to be trained, and the more training set samples are required. The amount of data in the BCI data set is not sufficient for the above classical models, as they are difficult to train and more likely to overfit. Therefore, this paper designs a smaller 2-layer convolutional neural network model as shown in Figure15, which contains 2 convolutional layers, 2 pooling layers, 2 fully connected layers, and 1 softmax output layer. The 2-layer CNN model reduces the size of the feature map and increases the number of channels by alternately computing convolutional layers and pooling layers. The output layer uses the softmax classification function to select an item with a higher probability as the classification result to perform a two-class judgment for the left hand and right hand on the input.Figure 15
Convolutional neural network structural model.The size of the input image in this paper is uniformly adjusted to64∗64, the kernel size is 4∗4 for the first convolutional layer, 2∗2 for the second convolutional layer, the number of kernels in the first convolutional layer is 8, and the number of kernels in the second convolutional layer is 8. The number of cores is 16 [29]. Considering various factors such as server memory size, input sample size, network model complexity, etc., the batch size is 8 and the epoch is 300 times in this experiment. In the training of neural networks, gradient descent algorithm is used iteratively. The adaptive learning rate for each parameter is obtained using the ADAM algorithm [30]. In addition, this paper introduces batch normalization (BN) to alleviate overfitting and uses a rectified linear unit (ReLU) function as an activation function in the convolutional layer to speed up network training. The comparison of classification accuracy of CNN, SA3, and SVM algorithms is shown in Figure 16.Figure 16
Comparison of classification accuracy of CNN, SA3, and SVM algorithms.


## 3.1. EMD

The key of EMD is to decompose a complex signal into a finite number of intrinsic mode function (IMF), and each of the decomposed IMF components contains local characteristic signals of different time scales of the raw signal.IMF has two limiting conditions [22]:
(1)
In the entire data sequence, the number of extreme points Ne and the number of zero crossing points Nz are equal or differ by at most 1, namely,(1)Ne‐1≤Nz≤Ne+1(2)
At any time pointt, the mean value of the upper envelope curve Mmaxt formed by the local maximum point and the lower envelope curve Nmint formed by the local minimum point is always zero. This is(2)Mmaxt+Nmint2=0Since empirical mode decomposition produces end effects in the decomposition process, which causes distortion problems near the two end points of the signal, a method to suppress the end effects is proposed. The decomposition steps of improved EMD are as follows:(1)
Take the initial end point of the raw signalXt as a maximum starting point, and mark it as M0, select an adjacent maximum point, and mark it as M1, connect the two points, and calculate its slope S1(3)S1=M1−M0tm1‐t0(2)
Calculate the minimum value obtained by the extension at the starting end point(4)N0=N1−S1•tn1−t0,in which, N1 refers to the minimum value closest to the starting point, and tn1 refers to the abscissa of N1. The abscissa of the minimum value is Tm0=t0(3)
Correspondingly, take the terminal end of the raw signal as the minimum point, and mark it asNQmin, select an adjacent minimum point, and mark it as NQmin−1, and calculate the slope S2 of the connection between the two points(5)S2=NQmin−1−NQmintnQmin−1−tQmin(4)
Calculate the maximum value obtained by terminal extension(6)MQmax=MQmax−1−S2•tmQmax−1−tNIts abscissa isTmQmax=tN.
(5)
Find the maximum value MAX and the minimum value MIN in the maximum value sequence. IfMQmax<MIN, define the maximum extension value as M1Qmax=terminalendpointvalue+max/2; if Ν0>max, define the maximum extension value as N10=startingendpointvalue+min/2. This method avoids the situation that the extended extreme value does not meet the interpolation requirements(6)
Find all the maximum and minimum points in the raw data signalXt, and use cubic spline interpolation function to fit the upper envelope Mmaxt and lower envelope Nmint.(7)
Calculate the local average value of the raw signalXt according to the upper and lower envelopes and record it as(7)m1t=Mmaxt+Nmint2(8)
The difference between the raw signal and the local mean is recorded as(8)h1t=Xt−m1t(9)
Ifh1t does not meet the two conditions of IMF, take h1t as the raw signal and go to (7) until the condition is met. At this time, h1t is the first IMF component. And then, obtain the remaining signal r1t by subtracting h1t from the raw signal Xt(9)r1t=Xt−h1t(10)
Regardingr1t as a new signal and repeating the above steps, n IMF components can be screened out(10)rnt=rn−1t−hntWhenrnt is a monotonic sequence, the screening ends.
(11)
Remove the high-noise components both IMF0 and IMF1(12)
The raw data sequence can be reconstructed as(11)X~t=∑i=2nhit+rntX~t is the reconstructed signal after EMD denoising.Extend the minimum and maximum points at both ends of the original signal, which not only improves the EMD method to better suppress the end effect of empirical mode decomposition but also reduces the intrinsic mode function components and postprocessing time.Figure5 is a schematic diagram of the EMD decomposition of the EEG data of the first subject and the instantaneous frequency of the corresponding IMF. It can be seen that the EMD method has good adaptability to the decomposition of EEG signals. After the raw signal is decomposed by EMD, 9 IMF components are obtained. Each order of IMF represents the characteristic signals of different time scales in the raw EEG signal. And the obtained IMF components are arranged in a descending order of frequency. Obviously, IMF0 and IMF1 are high-frequency noise. Therefore, when the EMD method is used to denoise the EEG signal, the denoising effect can be achieved by deleting the IMF component of the high-frequency part, and finally, the remaining part of the IMF is superimposed to obtain the denoised EEG signal.Figure 5
EMD decomposition results of the first subject’s EEG data.Figure6 is a comparison diagram of EEG signals before and after EMD denoising. It can be seen from Figure 6 that the EMD method removes the high-frequency noise in the raw signal, but it is still an unsmooth signal, which still contains more noise signals. The reason for this phenomenon is that the low-frequency part is not processed during the EMD denoising process, and the interference of the ECG and EMG has not been eliminated. Therefore, the EMD denoised EEG signal and the raw signal in the time domain do not seem to change much. The EMD method alone cannot effectively remove the high-frequency noise contained in the EEG signal.Figure 6
Comparison of EMD denoising results before and after.


## 3.2. ICA (Independent Component Analysis)

Assume thatSt=s1t,⋯,sntT is an N-dimensional vector composed of the raw signal source [23], where sit is its component, i=1,2,⋯,n. sit is linearly combined by the hybrid system A into an M-dimensional observation vector Xt=x1t,⋯,xmtT. That is,
(12)xit=ai1s1t+ai2s2t+⋯+ains1t=∑j=1naijsj,i=1,2,⋯,m.Formula (12) may be expressed as follows:
(13)Xt=A×St,in which A is an m×n mixing coefficient matrix. The task of ICA is to find an unmixing coefficient matrix W when the coefficient matrix A and the source matrix S are unknown and then separate the source signals St from the observed signal Xt, that is,
(14)Yt=WXt.in which, Yt=y1t,⋯,ymtT is the output after unmixing, W is a n∗m matrix. It is required that Yt be as close to St as possible. The mathematical model of ICA is shown in Figure 7.Figure 7
The mathematical model of ICA.Because each independent source signal and the mixing matrix are unknown, the output of the ICA has uncertainty, including the uncertainty of the amplitude (or variance) and the uncertainty of the output sequence.The independence criteria of ICA mainly include SOBI (Second-Order Blind Identification), FastICA (Hyvärinen’s Fixed Point Algorithm), information maximization (infomax), and JADE (Joint Approximation Diagonalization of Eigenmatrices) [24].This paper uses an improved FastICA algorithm. Specific steps are as follows.(1)
Suppose the raw EEG observation signalXt=x1t,⋯,xmtT, where xit represents the ith component of the raw EEG observation signal; m represents the number of components of the raw EEG observation signal; and t represents the sampling time of the raw EEG signal(2)
Demeaning for the raw EEG observation signalXt to obtain signal X’t:(15)X′t=Xt−EXt(3)
Use the following formula (16) to decompose the processed EEG observation signal X′t into signals Zt=z1t,⋯,zmtT that are not correlated with each other:(16)Zt=B~X′t,(17)B~=Ds−σ2I‐1/2UsT,(18)σ≅EX′t∗X′tTI,in which B~ represents the projection factor; Ds is a diagonal matrix whose diagonal elements are eigenvalues of the covariance matrix CX=EX′t∗X′tT of X’t; and Us is a matrix with the unit norm eigenvector of CX as a column
(4)
Set the initial value of the unmixing matrix toW(5)
The optimal valueW∗ of the unmixing matrix is as follows:(19)W∗=EZtGWTZt−EG′WTZtW(6)
The optimized valueW∗ of the unmixing matrix is subjected to decorrelation and normalization processing, and the processed unmixing matrix W∗¯ is obtained:(20)W∗¯=W∗W∗(7)
Determine whether the processed unmixing matrixW∗¯ is convergent. If it is convergent, the matrix W∗¯ is the final unmixing matrix. If not, W=W∗¯, go to (5)(8)
Use formula (21) to process the EEG signal Zt=z1t,⋯,zmtT obtained in step (3):(21)Yt=W∗¯∗Zt,in which Yt=y1t,⋯,ymtT and y1t,⋯,ymt are not related to each otherAn object is selected from the Graz data set A as the input signal. Figure8 shows the raw EEG signal before FastICA containing the artifacts such as EOG and ECG. It can be clearly seen that the artifacts seriously interfere with the EEG signal. Figure 9 shows the clean EEG signal by FastICA. Comparing the EEG signal of Figure 8 with Figure 9, it is easy to draw the conclusion that the artifacts such as EOG and ECG are effectively removed by FastICA.Figure 8
Raw EEGs before FastICA.Figure 9
EEGs after FastICA.


## 3.3. EMD-ICA Method

In this paper, the EMD-ICA joint algorithm is used to denoise the EEG signal. The algorithm steps are as follows:(1)
The raw EEG signal is decomposed inton independent components X1,X2,⋯,Xnusing the improved FastICA algorithm in Section 3.2(2)
Decompose each independent componentXi into multiple intrinsic mode functions (IMFs) and margins using the improved empirical mode decomposition method (EMD) in Section 3.1(3)
Delete IMF0 and IMF1 components with higher noise(4)
Superimpose and reconstruct the decomposed IMF component and margin to obtainY(5)
Use the inverse improved FastICA algorithm in Section3.2 to reconstruct signal Y to obtain cleaned EEG signal ZThe EMD-ICA joint algorithm is used to preprocess the EEG signal. The joint algorithm eliminates the interference of noise to the maximum on a finer scale, while retaining most of the useful EEG signals, so as to better complete the purpose of denoising for EEG signals. The denoising process of EEG based on the above method is shown in Figure10.Figure 10
EMD-ICA fusion denoising process.The raw data of the first subject is selected for EMD-ICA decomposition, and the experimental results are shown in Figure11. It can be seen that the raw EEG signal is not smooth and is mixed with noise signals with higher amplitude. After EMD-ICA denoising the raw data, the obtained signal is more in line with the natural trend of the original data. This indicates that the EMD-ICA fusion algorithm not only accurately eliminates the artifact components but also better retains the local characteristics of the raw EEG.Figure 11
Experimental results of EMD-ICA fusion method.


## 3.4. Continue Wavelet Transform

Continuous wavelet transform is to project a one-dimensional signal onto a two-dimensional time-scale plane, which is established by translation and scaling of a mother wavelet. The definition of the wavelet function is as follows [25, 26]:
(22)∫‐∞+∞ϕtdt=0, in which, ϕt is the mother wavelet function, ϕt=1. The wavelet function generated by ϕt can be expressed as
(23)ϕa,bt=1bϕt−ab,in which a and b are translation and scale parameters, respectively.For any signal functionXt, decompose Xt on this wavelet function, and its wavelet transform can be obtained as
(24)WXa,b=X,ϕa,b=∫−∞+∞Xt1bϕt−ab¯dt,in which W is the result of wavelet transform; X,ϕa,b represents the inner product of X and ϕa,b.The formula for the reconstruction of the signalXt is
(25)Xt=1λ∫0+∞∫−∞+∞WXa,b1bϕt−abdadbb2,where
(26)λ=∫0+∞ϕ∧ω2+ωdω,in which ω represents the analog angular frequency. ϕ^ω is the Fourier transform of the mother wavelet function ϕx.The diversification of the mother wavelet function is one of the main advantages of the wavelet transform. Commonly used mother wavelet functions include Haar wavelet, Morlet wavelet, Mexican Hat (Mexh) wavelet, and Daubechies (DbN) wavelet [27].Left-hand and right-hand motor imagery tasks led to ERD and ERS phenomena on the left and right sides of the cerebral motor cortex, respectively, which affected the EEG signals of the C4 and C3 electrodes. Figures12 and 13 are wavelet time-frequency maps of the C3 and C4 channels during the left-hand motor imagery task.Figure 12
Wavelet time-frequency map of C3.Figure 13
Wavelet time-frequency map of C4.It can be seen from Figures12 and 13 that the energy of the C4 channel dropped significantly at the 3rd second after the start of the experiment and recovered after a period of time; that is, the ERD phenomenon occurred. However, the energy of the C3 channel remains high rather than decreasing, which is known as the ERS phenomenon. The maps of C3 and C4 are combined into the final input image, as shown in Figure 14.Figure 14
The combined input map of C3 and C4.


## 3.5. Convolutional Neural Network

As a multilayer neural network, convolutional neural network (CNN) has been successfully applied in the fields of computer vision and image processing due to its powerful feature extraction ability. CNN usually consist of multiple convolutional and pooling layers. A typical CNN has the following layers: input layer, convolutional layer, pooling layer, fully connected layer, and output layer. In the convolution layer, the convolution filter is used to perform convolution filtering on the input data, which has the characteristics of local linking and weight sharing. Compared with the fully connected layer, the convolutional layer can extract the local information of the input data and reduce the number of parameters. In the pooling layer, the method of maximum pooling or average pooling is usually used to reduce the dimension of features, which reduces the complexity of network computation [28].Compared with traditional algorithms, the ability of deep learning algorithm to extract features is greatly improved, and usually, the more complex the network, the more sufficient the features extracted, and the better the obtained classifier effect. However, the advantages of deep learning algorithm classification accuracy are usually only reflected when the number of sample sets is large enough, and the more complex the network, the more parameters to be trained, and the more training set samples are required. The amount of data in the BCI data set is not sufficient for the above classical models, as they are difficult to train and more likely to overfit. Therefore, this paper designs a smaller 2-layer convolutional neural network model as shown in Figure15, which contains 2 convolutional layers, 2 pooling layers, 2 fully connected layers, and 1 softmax output layer. The 2-layer CNN model reduces the size of the feature map and increases the number of channels by alternately computing convolutional layers and pooling layers. The output layer uses the softmax classification function to select an item with a higher probability as the classification result to perform a two-class judgment for the left hand and right hand on the input.Figure 15
Convolutional neural network structural model.The size of the input image in this paper is uniformly adjusted to64∗64, the kernel size is 4∗4 for the first convolutional layer, 2∗2 for the second convolutional layer, the number of kernels in the first convolutional layer is 8, and the number of kernels in the second convolutional layer is 8. The number of cores is 16 [29]. Considering various factors such as server memory size, input sample size, network model complexity, etc., the batch size is 8 and the epoch is 300 times in this experiment. In the training of neural networks, gradient descent algorithm is used iteratively. The adaptive learning rate for each parameter is obtained using the ADAM algorithm [30]. In addition, this paper introduces batch normalization (BN) to alleviate overfitting and uses a rectified linear unit (ReLU) function as an activation function in the convolutional layer to speed up network training. The comparison of classification accuracy of CNN, SA3, and SVM algorithms is shown in Figure 16.Figure 16
Comparison of classification accuracy of CNN, SA3, and SVM algorithms.


## 4. Results and Analysis

In order to quantitatively evaluate the MD-ICA method algorithm proposed in Section3.3, the raw signal was subjected to denoising preprocessing by three different methods: ICA, ICA-WT, and EMD-ICA. The whole experiment process includes 120 groups of experiments, and 120 groups of vector matrices are obtained, which are divided into 40 groups of test vector matrices and 80 groups of training vector matrices. The classification results are shown in Table 1.Table 1
Comparison of classification accuracy of ICA, ICA-WT, and EMD-ICA denoising algorithms.
MethodsClassification accuracy (%)WorstBestMeanICA63.5576.3071.21ICA-WT78.7083.4980.26EMD-ICA80.5385.1682.45It can be seen from Table1 that the denoising effect of only using the ICA method is the worst. Although the denoising effect of ICA-WT is better than that of only using ICA, the final classification accuracy is still not as good as that of the EMD-ICA algorithm. The experimental results show that the EMD-ICA algorithm proposed in this paper can effectively improve the signal-to-noise ratio, obtain an ideal denoising effect, and at the same time can preserve the effective components of the signals; thus, the EMD-ICA algorithm can effectively improve the classification accuracy.In Section3.4, time-frequency analysis of EEG signals was performed using continuous wavelet transform. The choice of mother wavelet is flexible, which is the advantage of wavelet transform and is an important part of wavelet transform. At present, there is no clear theoretical guidance for the selection of mother wavelet. Generally, these aspects including orthogonality, support length, symmetry, vanishing moment, regularity, and similarity of the wavelet function must be considered. In addition, the calculation amount of different mother wavelet functions often varies greatly, so for some algorithms that have higher requirements on calculation speed, the calculation amount is also a factor to be considered. The classification results using different mother wavelet functions are shown in Table 2. It can be clearly seen that the results obtained by different mother wavelets are very different. The classification results based on Morlet wavelet transform are relatively good, and the calculation speed is lower than that of Mexh and db4 wavelet functions. Finally, this paper decides to take Morlet wavelet as the mother wavelet function.Table 2
Comparison of different mother wavelet functions.
MethodsClassification accuracy (%)Average times (s)WorstBestMeanHaar73.5378.2775.1976.37Morlet81.1285.4983.2679.53Mexh80.4483.1681.4087.09db476.3982.3679.7180.46EMD-ICA denoising was performed on the EEG data of 9 subjects, and the EEG signals were converted into two-dimensional time-frequency images using continuous wavelet transform with Morlet as the mother wavelet. Then, the time-frequency images are classified using the neural network described in Section 3.5 of this paper. At the same time, support vector machines (SVM) and stacked autoencoders (SAE) are also applied to the same input image for classification. As can be seen from Table3 and Figure 16, under the same input data, the average classification accuracy of the CNN method proposed in this paper is higher than that of SVM and SAE.Table 3
Comparison of classification accuracy of different classification algorithms.
SubjectsClassification accuracy (%)CNNSAESVM185.4475.9281.73283.6568.8183.16386.3275.668.65492.3193.4595.38583.879.6481.7685.6773.5983.74775.5367.2976.04882.5475.673.91980.4673.8182.35Average83.9775.9780.74Table4 shows the comparison of kappa values calculated by CNN, SAE, and SVM methods; the kappa value is a classification performance measure that removes the influence of random classification accuracy [31]. The kappa calculation formula is
(27)Kappa=p0−pc1−pc,in which p0 is classification accuracy, and pc is the result of random classification; the value of pc is 0.5 for two kinds of classification.Table 4
Kappa value calculated by different methods.
SubjectsBest kappa valueCNNSAESVM10.6120.3300.56320.2590.3570.18030.5780.3560.26740.8880.8800.88550.7150.5100.69760.4820.5680.52270.3880.3300.39680.5530.5760.42690.3110.3650.573Average0.5320.4750.501It can be seen from Table4 that compared with the SAE and SVM algorithms, the CNN algorithm has a higher kappa value, which further proves that the results obtained by the convolutional neural network classification algorithm are better than those obtained by the traditional algorithm to a certain extent.In summary, on most subjects, the algorithm in this paper can effectively improve the classification and recognition accuracy of motor imagery EEG signals. However, due to the individual differences of EEG data, the algorithm in this paper may not be ideal for individual subjects.


## 5. Conclusion

Brain-computer interface technology has innovatively reformed the way of information interaction between people and the outside world. It repairs and expands the physiological and cognitive functions of the human body. Preprocessing of EEG signals is a hot issue in brain-computer interface research [32]. EEG signal is a bioelectric signal with strong background noise, including 60 Hz power frequency interference and physiological artifacts produced by its own physiological activities. These artifacts seriously affect the feature extraction and classification of motor imagery EEG signals. This paper focuses on how to effectively reduce the physiological artifacts in EEG signals so as to obtain cleaning EEG. Through taking full advantage of EMD and ICA, a new method is proposed to eliminate physiological artifacts. Experimental results show that the fusion method not only effectively removes the artifacts but also preserves the weak EEG activity signals missed in the artifacts. Then, continuous wavelet transform was utilized to extract energy features of μ rhythm and β rhythmic to represent the characteristics of EEG signals under different motor imageries. These two features are normalized and used as the input of the 2-layer convolutional neural network designed by the paper, and the two kinds of features are learned by using the CNN, and then, the two-classification problem of motor imagery EEG signals is completed. Experiments are carried out with the Graz data set A as the data set. The experimental results show that, comparing SA3 and SVM algorithms, the average classification accuracy and kappa value of this paper proposed method are higher than those of SVM and SAE for most subjects. It is further demonstrated that the proposed algorithm can effectively remove physiological artifacts while retaining useful EEG components. The research in this paper has certain significance for the research and application of BCI system. It also provides an idea and method for clinical medicine and brain science research. Deep learning has the advantage of automatically learning complex features from massive data. However, it is difficult for a user to generate a motor imagery EEG data set that can train a deep learning model, so how to train an effective network model with fewer samples is a problem to be solved.



---

*Source: 8935543-2022-03-24.xml*