Exploring spatial-frequency-sequential relationships for motor imagery classification with recurrent neural network

Background Conventional methods of motor imagery brain computer interfaces (MI-BCIs) suffer from the limited number of samples and simplified features, so as to produce poor performances with spatial-frequency features and shallow classifiers. Methods Alternatively, this paper applies a deep recurrent neural network (RNN) with a sliding window cropping strategy (SWCS) to signal classification of MI-BCIs. The spatial-frequency features are first extracted by the filter bank common spatial pattern (FB-CSP) algorithm, and such features are cropped by the SWCS into time slices. By extracting spatial-frequency-sequential relationships, the cropped time slices are then fed into RNN for classification. In order to overcome the memory distractions, the commonly used gated recurrent unit (GRU) and long-short term memory (LSTM) unit are applied to the RNN architecture, and experimental results are used to determine which unit is more suitable for processing EEG signals. Results Experimental results on common BCI benchmark datasets show that the spatial-frequency-sequential relationships outperform all other competing spatial-frequency methods. In particular, the proposed GRU-RNN architecture achieves the lowest misclassification rates on all BCI benchmark datasets. Conclusion By introducing spatial-frequency-sequential relationships with cropping time slice samples, the proposed method gives a novel way to construct and model high accuracy and robustness MI-BCIs based on limited trials of EEG signals.

conventional classification of EEG signals was usually adapted by "shallow" classifiers (linear discriminant analysis (LDA), support vector machine (SVM), and neural network (NN)) [21][22][23][24][25], such classifiers are appropriate for small sample size. Hence, a complete entity of each motor imagery trial's spatial-frequency features was fed into these classifiers for classification. Due to the difficulty of obtaining motor imagery trials, public or private datasets have limited amounts of EEG trials from MI-BCIs [26,27]. Thus, "shallow" classifiers with less data will produce poor performances of classification.
Second, except for spatial-frequency features, EEG signals' sequential relationship is another useful feature for motor imagery classification. By cropping the spatialfrequency features into several time slices, each time slice can be treated as time-series, which contains sequential relationships over time. If the sequential relationships can be modeled by classifiers, the novel spatial-frequencysequential relationships will significantly improve the performances and robustness of motor imagery classification.
To solve the two major challenges, this paper introduces a deep recurrent neural network (RNN) architecture for the classification based on FB-CSP algorithm [28,29]. Also, by modeling EEG signals by RNN, an optimal number of hidden layers are obtained for RNN. Then, a sliding window cropping strategy (SWCS) is used to crop the entity trial into several time slices to increase the number of samples by the optimal number. Since the deep neural networks have dramatically improved the state-ofthe-art methods in signal processing and classification, researches on EEG signals have been developed by using deep learning techniques to extract essential feature representations. The sequential relationships of EEG signals are easy to be extracted by RNN architecture. Therefore, the two contributions of this study are as follow:

A deep RNN architecture is applied to the FB-CSP
features to extract the spatial-frequency-sequential relationships for motor imagery classification. The abundant features will improve the performances of classification. Also, two different memory units, long short-term memory (LSTM) unit [30] and gated recurrent unit (GRU) [31], are included in the RNN architecture. 2. The FB-CSP features extracted from a complete entity motor imagery trial are cropped by the SWCS with an optimal number. The strategy will increase a large a e deep neural networks.

Conventional methods
Manual feature extraction methods and shallow classifiers are developed for conventional motor imagery classification. These features are usually extracted from the spatial-frequency features and sequential relationships of EEG signals. Table 1 illustrates the related work regarding feature extraction methods and the corresponding classifiers in the state of the art methods. From Table 1, we found CSPs algorithm [8,21] is the key algorithm for extracting spatial features in motor imagery classification. Other researchers improve the CSPs algorithm by a probabilistic model [32] or the genetic algorithm (GA) [33]. Except for spatial features, the frequency features of power spectrum density (PSD) and sequential relationships of adaptive auto regression (AAR) are also used in motor imagery classification [22,34]. Besides, the "time-frequency" features combine frequency features and sequential relationships for classification [17]. For the classification, conventional classifiers focus on shallow machine learning models. In some cases, the pre-processing algorithm multivariate empirical mode decomposition (MEMD) has been used to improve signalnoise ratio and classification accuracy [25]. The related works used manual features and shallow classifiers for the following reasons: on the one hand, because public datasets have limited EEG samples, they are more suited for classification by LDA/SVM/Naive Bayes classifiers; on the other hand, the EEG signals are regarded as a complete entity, and the entity is classified by spatial, frequency features or sequential relationships. However, if signals belong to time-series data, sequential relationships over time will provide the discriminant features for motor imagery classification.

Deep learning methods
Statistical, integrated, and deep learning are the common classification methods in machine learning [35,36].
In particular, deep learning classification methods have been used gradually for EEG signal classification [37][38][39]. Table 2 illustrates the related works regarding the state of the art of deep learning classifiers. From Table 2, we conclude that deep learning is widely used in EEG signal classification. Convolution Neural Network (CNN) models [40][41][42][43][44] and Deep Belief Network (DBN) models [32,45,46] are most often used in the analysis of EEG signals. Actually, the CNN and DBN models are used to extract the spatial features from EEG signals. These two deep learning models still treat the complete entity trials for classification, so the performance can't be improved much. However, the deep RNN architecture can extract the sequential relationships from EEG signals [47,48]. By using a sliding window cropping strategy, the complete entity trials will be cropped into several time slices for classification. Several multiples growth number of samples by cropping for classification will obtain a significant performance improvement of motor imagery classification. Therefore, the discriminant features for motor imagery classification are extracted by

Methods
By considering references [47,49], our proposed method regards EEG signals as time-series data, and the extracted spatial-frequency features' sequential relationships are represented by RNN architecture. Due to the fact that conventional FB-CSPs algorithms with shallow classifiers do not contain the sequential relationships, and these algorithms regarded the entity of each trial as a single sample for classification. Therefore, two methods are developed to validate and represent spatial-frequencysequential relationships for classification. First, we test a group of smoothing time windows on the FB-CSP features to validate whether the sequential relationships can improve the classification performance of EEG timeseries. Then, a deep RNN architecture is applied to represent spatial-frequency-sequential relationships on FB-CSP features for classification. It is easy to cause overfitting problems and drop the classification performance if the deep neural networks are presented for classification by the entity of a trial [50]. Therefore, before using the deep RNN architecture, a sliding window cropping strategy is applied to crop the entity of each trial into several time slices. Then, each time slice is fed into the deep RNN architecture for the motor imagery classification. The size of each time slice will be set as the same of the optimal number of hidden layers for RNN to obtain the optimal classification performance. The proposed method is illustrated in Fig. 1. In Fig. 1, our proposed method is comprised of four progressive stages of signal processing and machine learning on EEG signals: (1) a filter bank comprising multiple Butterworth band-pass filters to extract frequency features, (2) a CSP algorithm is used to extract spatial features, (3) a sliding window cropping strategy is applied to crop time slices to model the sequential relationships of spatial-frequency features, (4) classification of the spatial-frequency-sequential relationships on time slices by a deep RNN architecture. In the deep RNN architecture, two different memory units, GRU and LSTM unit, are included to compare classification performance and robustness. The CSP projection matrix for each filter band, the discriminative spatial-frequency features, and the deep RNN architectures are computed and trained from training data labeled with the respective motor imagery action. These parameters computed from the training phase then used to validate each single-trial motor imagery action. By using the same cropping strategy in the validation phase, the classification of single-trial motor imagery action will predict several targets. The final evaluated action will be obtained by averaging all predicted targets. The procedure of our proposed methods. Our proposed method is comprised of four progressive stages of signal processing and machine learning on EEG signals: (1) a filter bank comprising multiple Butterworth band-pass filters to extract frequency features, (2) a CSP algorithm is used to extract spatial features, (3) a sliding window cropping strategy is applied to crop time slices to model the sequential relationships of spatial-frequency features, (4) classification of the spatial-frequency-sequential relationships on time slices by a deep RNN architecture. In the deep RNN architecture, two different memory units, LSTM unit and GRU , are included to compare classification performance and robustness

Spatial-frequency features
The widely used spatial-frequency features extraction algorithm for classification of motor imagery EEG signals was Filter Bank Common Spatial Patterns (FB-CSP) algorithm [8,9]. There are two steps in the FB-CSP method: (1)a group of band-pass filters are presented to the raw EEG data to obtain the subject-specific frequency band.
(2)The CSP algorithm is provided to every filter result to extract the optimal spatial features. Then, a classifier is used in all of the FB-CSP features for motor imagery classification.
To extract CSP features, let X c ∈ R N * T represent one band-pass filtering result, where c is the number of classes, N is the number of potentials of EEG, T is the number of samples in each trial. Each dataset contains L trials of EEG signals and each signal X c is a zero average signal. The purpose of the CSP algorithm is to find an optimal spatial vector, − → w ∈ R M×N , to project the original EEG signal to a new space to obtain good spatial resolution and discrimination between different classes of EEG signals. To calculate the optimal projection matrix, let the average covariance matrix of class "c" be C c , and average power of class "c" be P c = − → w T C c − → w . For an example of two classes on the minimized projected − → w axis, the maximization of the power ratio is written into the Rayleigh quotient form: The Rayleigh quotient is then re-translated into a constrained optimization problem, which is then solved by applying the Lagrange multiplier method to the problem. The optimization results include both eigen-vectors and eigen-values. The optimal CSP spatial filter vector, − → w * ∈ R M×N , is constructed by taking M = 2m, M ≤ N eigenvectors corresponding to the "m" largest and "m" smallest eigen-values: where − → w λ i is the eigen-vector that corresponds to the eigen-value λ i . Each filter band of EEG signals, X c , is spatially filtered by: where Z c ∈ R M×T is the spatial-frequency features. The EEG signals are composed of rapidly changing voltage values; therefore, band power (variance) is used as a feature for the classifier. For multi-class extension to the FB-CSP algorithm, the one-versus-rest (OVR) strategy is presented to solve the multi-class motor imagery BCI classification.

Spatial-frequency-sequential relationships
Conventional algorithms for motor imagery EEG signals fed spatial-frequency features (FB-CSP) into classifiers to discriminate different motor imagery targets. In this paper, the FB-CSP features are fed into a deep RNN architecture to get spatial-frequency-sequential relationships to improve the classification performance of motor imagery. To validate and represent the spatial-frequencysequential relationships, a group of smoothing time windows are put on the FB-CSP features to validate the effect of sequential relationships, and a RNN model with sliding window cropping strategy is applied to represent spatialfrequency-sequential relationships on EEG time-series. To improve the classification performance and overcome the over-fitting problem, the LSTM unit and GRU are used to construct LSTM-RNN architecture or GRU-RNN architecture for EEG signals classification.

Smoothing time windows on FB-CSP features
Since the FB-CSP features are extracted from EEG timeseries, such features also contain sequential relationships. Before we represent the sequential relationships by the RNN architecture, a group of smoothing time windows are put on the FB-CSP features to smooth the sequential relationships. For the classification by FB-CSP features, we will adjust the smoothing time window size, and find the influence of classification performance by the smoothing time windows. If the influence for the performance is large, the sequential relationships on FB-CSP features will be validated to influence the classification performance. Therefore, the RNN architectures with LSTM and GRU memories will be applied to extract spatial-frequencysequential relationships on FB-CSP features. According to the smoothing process, given smoothing window size, ω, the following smoothing operation is applied to the FB-CSP features, Z c : where Z c (t) is the smoothed FB-CSP features. In the experiments, we adjust the parameter ω to obtain different smoothing levels of FB-CSP features, and get the classification performance by support vector machine (SVM). The classification performance will validate and instruct the sequential relationships for EEG signals classification.

RNN architectures with LSTM and GRU memories
To represent the spatial-frequency-sequential relationships, we introduce the RNN architecture in this study [51,52]. The RNN architecture, containing an input layer, recurrent hidden layers and an output layer, is widely used to represent time-series [53,54]. In recurrent hidden layers, a number of simple computation units with weighted interconnections, including delayed feedback [28]. The feedback will give intrinsic states and learn tasks from memory, which is suitable for modeling EEG signals. With the activation functions, the deep RNN architecture is good at learning sequential patterns from EEG signals. Figure 2 illustrates the standard deep RNN architecture. In the figure, the simplified RNN architecture is shown in the box on the left. The box on the right shows the architecture unfolded in a form of time-series [ · · · , t − 1, t, t + 1, · · · ]. In the form of time-series hidden layers, the input of layer "t" contains the output of layer "t + 1", so do the input of layer "t − 1". The sequential relationships propagate from the end of the time-series to the start of the time-series by neurons, which are connected by horizontal lines in the figure.
Recurrent connections between hidden layers are followed by a feed-forward output layer. Hence, the deep RNN architecture is universal approximators of finite states. Therefore, a deep RNN architecture can approximate any finite states with enough recurrent hidden layers and trained weights. Let Z c ∈ R M * T represent FB-CSP features, where M is the features dimension, T is the number of samples in each trial. The RNN architecture can be defined as: where x t is the vector of input layer, which is one of the time slices of the FB-CSP features Z c ∈ R M * T . h t is the vector of hidden layer. y t is the vector of output layer. W , U and b are the recurrent connected weights. σ is the activation functions. Neural networks are processed by back-propagations (BP) algorithm in common. For the RNN architecture, the sequential relationships propagate all steps back through time, so the feedback of hidden layers will be processed by back-propagation through time (BPTT) algorithm [55]. The training procedure of a deep RNN architecture is performed using a stochastic gradient descent (SGD) algorithm. By using SGD algorithm, we can iteratively update the network's weight values based on BPTT algorithm. However, the BPTT algorithm is too sensitive to recent distractions; thus, the error flow tends to vanish as long as the weights have absolute low variations, especially at the onset of the training phase. Long short-term memory (LSTM) unit [30] and Gated recurrent unit (GRU) [31] are proposed to overcome the vanishing gradient problem. The LSTM and GRU architecture is illustrated in Figs. 3 and 4. The introduction of these two architectures are as follow: 1. LSTM architecture: In a LSTM unit [56], input, output and forget gates are used to retain memory contents; these gates also prevent the irrelevant inputs and outputs from entering the memory. Thus, the unit stores the long term memory features of the time-series data. A peephole method [57] will be included in the LSTM architecture to transfer memories for all gates. 2. GRU architecture: A GRU supports each recurrent unit to adaptively obtain dependencies of different time scales. The GRU has "update" and "reset" gates to prevent the error flow of information in the unit. Similarly to the LSTM unit, the gates prevent irrelevant inputs and outputs.
In such "memory units", because these special units have internal states, multiplicative gates are employed to enforce constant error flow. These two different memory units are used in the deep RNN architecture to classify motor imagery tasks through spatial-frequencysequential relationships. For each hidden layer of the RNN architecture, the original hidden layer will be replaced by LSTM unit or GRU to construct LSTM-RNN architecture or GRU-RNN architecture. Classification results are compared and analyzed to show which memory unit is more suitable for MI-BCI.

Sliding window cropping strategy
The conventional trial-wise EEG signals classification algorithms treat the entity duration of a trial as a single sample and the corresponding label as a single target. Then, a shallow classifier is used to train and validate motor imagery tasks. The conventional algorithms will lead to less samples and high dimensionality of features, which will cause the over-fitting problem and drop the accuracy of classification. In this study, a deep RNN architecture is used for the classification of EEG signals, if the entity duration of a trial is fed into deep RNN architecture, the number of hidden layers will be too large to get long-term patterns for the classification of EEG signals. To avoid the over-fitting problem of classification, a sliding window cropping strategy is applied to each trial to crop the entity duration of the trial into several time slices, and the label of the trial will be repeated to all time slices. This strategy will increase the number of training samples for the RNN architecture, which is widely used in the recognition tasks of image, audio and EEG signals by neural networks [58][59][60]. Fig. 3 The LSTM unit architecture. In a LSTM unit, input, output and forget gates are used to retain memory contents; these gates also prevent the irrelevant inputs and outputs from entering the memory. Thus, the unit stores the long term memory features of the time-series data Fig. 4 The GRU architecture. A GRU supports each recurrent unit to adaptively obtain dependencies of different time scales. The GRU has "update" and "reset" gates to prevent the error flow of information in the unit. Similarly to the LSTM unit, the gates prevent irrelevant inputs and outputs In our study, let Z c ∈ R M * T represents the inputs of RNN, the entity duration of a trial includes T time steps. Assumed τ is the cropping size of the sliding window cropping strategy, the time slices of the trial by cropping can be defined as The number of training samples will be increased T − τ times, and all time slices will get the label y c as the same label from the original trial. Since the deep RNN architecture has the ability to extract signals' sequential relationships for classification, we treat the number of hidden layers as the size of time slices. Therefore, we need to confirm the optimal number of hidden layers of the deep RNN architecture for motor imagery EEG signals classification; then, the optimal cropping size will be obtained from the EEG modeling experiment. If the optimal number of hidden layers is confirmed, the cropping size is confirmed. In common, the trial duration used for motor imagery is two seconds, and we obtain 500 samples for a 250 Hz sample rate. If the optimal number of the hidden layers is 20, the original trial will be crop to 480 time slices. The sliding window procedure for cropping a trial into time slices is shown in Fig. 5.

Fig. 5
The sliding window procedure for cropping a trial into time slices. In common, the trial duration used for motor imagery is two seconds, and we obtain 500 samples for a 250 Hz sample rate. If the optimal number of the hidden layers is 20, the original trial will be crop to 480 time slices

Experimental datasets setup
The performances of the algorithms were evaluated on the BCI Competition IV [27] "Dataset 2a" and "Dataset 2b" 1 . The two datasets are compared in Table 3. Figure 6 illustrates how the single-trial EEG data were extracted on "Dataset 2a" and "Dataset 2b". The two datasets share the same procedure. In the motor imagery classification experiments, each subject sat in a soft chair comfortably facing a computer screen. The BCI Competition IV experiments are composed of the following six steps: (1) Each trial started with a warning tone. (2) Simultaneously, a fixation cross was shown on the computer screen for two seconds.
(3) After two seconds, a cue, in the form of an arrow, was randomly shown in lieu of the fixation cross, and the subjects started the corresponding motor imagery task of the cue. (4) After another 1.25 s, the cue reverted to the fixation cross. (5) The motor imagery task continued until the sixth second, at which time the fixation cross disappeared. (6) Finally, there was a short 1.5 s break. The signals were sampled at 250 Hz and recorded. The preprocessing operations on the signals for notch filtered and band-pass filtered were 50Hz and 0.1-100Hz, respectively. The BCI Competition IV "Dataset 2a" is composed of the following four classes of motor imagery EEG measurements from nine subjects: (1) left hand, (2) right hand, (3) feet, and (4) tongue. Two sessions, one for training and another for evaluation, were recorded from each subject.
"Dataset 2b" is composed of two classes of motor imagery EEG measurements from nine subjects: (1) left hand and (2) right hand. Five sessions, the first three for training and the last two for evaluation, were recorded from each subject. According to the extraction procedure, the time range [4, 6s] was chosen for motor imagery classification because of a strong ERD/ERS phenomenon within that range [12,44].
The spatial-frequency features are extracted by the FB-CSPs algorithm. In the division of the whole band (8-30Hz, covered μ and β rhythms) to obtain universality for all subjects, the optimal band width range is 4Hz overlaps the next by 2Hz [5,25]. The optimal division of band-pass filters is shown in Table 4. After the optimal frequency bands filter the raw EEG signals, the CSP algorithm is applied to the filtered EEG signals to obtain spatial-frequency features. In (2) in the CSP algorithm, parameter m for processing "Dataset 2a" and "Dataset 2b" is set to 2 and 1, respectively.
After extraction of spatial-frequency features, two separate experiments to confirm the parameters and validate the performances of spatial-frequency-sequential relationships and the classification of motor imagery are as follows: 1. EEG modeling experiments: First, a size range of [0, 4] smoothing time windows are put on the FB-CSP features to obtain the performance of classification.  Table 4 Optimal division of band-pass filters classification, the spatial-frequency FB-CSP features are fed into the deep RNN architecture to obtain spatial-frequency-sequential relationships. The spatial-frequency features are cropped by a sliding window sized by the optimal number of hidden layers. In the classification by LSTM-RNN architecture and GRU-RNN architecture, the accuracies, errors and efficiency of classification will be compared between spatial-frequency features and spatial-frequency-sequential relationships.

EEG modeling experiments and results
To obtain the performance of classification influenced by the sequential relationships, a group of smoothing windows with the size range [0 ,4] is presented to FB-CSP features. In our experiments, via smoothed FB-CSP features, the SVM classifier with RBF kernel is used for motor imagery classification. Figure 7 illustrates the smoothing time window experimental results for "Dataset 2a" and "Dataset 2b". Among the results, "SW=0" expresses the FB-CSP features without smoothing. From the results, we find that the performance of EEG signals classification was fully influenced by the smoothing time windows. Thus, the RNN architecture is introduced in this study to extract spatial-frequency-sequential relationships from FB-CSP features for classification. However, we must validate the presentation of spatial-frequency-sequential relationships by a RNN architecture at first. There are three steps to validate the presentation of spatial-frequency-sequential relationships by RNN architecture. First, to validate whether the deep RNN architecture can model EEG signals or not, we train a deep RNN architecture by 200 iterations of SGD algorithm over 22 channels of the first three seconds of EEG signals from "Dataset 2a Subject 3". To test the modeling ability, the previous outputs are fed back into model's inputs to predict the current EEG signals. The results on channel "C3" by 20, 30, 90 hidden layers are drawn in Fig. 8. From the results, we find the deep RNN architecture will predict the same level of signals as the number of hidden layers Second, we evaluate the classification performances of the deep RNN architecture by 200 iterations of SGD algorithm over the training data of "Dataset 2a Subject 3". The loss function for EEG signals by RNN architecture is the logarithmic cross-entropy, which is defined as [61]: where c t is the ground truth result, and b t c is the prediction result by deep classifier. The number of iterations for optimizing the loss function is an experience value of controlling training epochs by limiting the number of hidden layers. Figure 9 gives the training and validation crossentropies as the number of hidden layers increased. From the results, we find the training and validation crossentropies have separations over 20 hidden layers. The cross-entropies will not reduce if the signals are overfitted by the RNN architecture. In fact, the cross-entropies will not increase, so the deep RNN architecture continues to learn components of the signals that are common to all of the EEG sequences. Compared with LSTM-RNN architecture and GRU-RNN architecture, the LSTM-RNN architecture needs more hidden layers to achieve a same level of cross-entropy during the classification of EEG signals.
Third, since a large number of hidden layers requires much computational complexity, and causes the overfitting problem to achieve low validation accuracies, Fig. 10 gives the training and validation accuracies as the number of hidden layers increased. From the results, we find the validation accuracies appear peaks with a 20-15 hidden layers. Compared with LSTM-RNN architecture and GRU-RNN architecture, the LSTM-RNN architecture needs more hidden layers to achieve a same level of accuracy during the classification of EEG signals. When the deep RNN architecture is over-fitting, the accuracy of GRU-RNN has a sharp drop than LSTM-RNN.

Classification experiments and results
Let Z c ∈ R M * T represents the spatial-frequency features, where M is the feature dimension, T is the number of samples in each trial. After EEG modeling experiments, the optimal number of hidden layers for LSTM-RNN and GRU-RNN of all subjects are confirmed. Then, the optimal number τ is used for cropping training set and validation set by sliding window cropping strategy. Hence, the samples of each trial in training set and validation set

Number of Hidden Layers Average Cross-Entropy
Training Validation Fig. 9 The curves of LSTM-RNN and GRU-RNN's training and validation cross-entropies as the number of hidden layers increased. The training and validation cross-entropies have separations over 20 hidden layers. The cross-entropies will not reduce if the signals are over-fitting by RNN. In fact, the cross-entropies will not increase, so the deep RNN architecture continues to learn components of the signals that are common to all of the EEG sequences. a The curves of LSTM-RNN and b The curves of GRU-RNN will be increased T − τ times to satisfy the deep RNN architecture. After training procedure, the validation procedure will produce T − τ classification targets in each trial. Finally, the unique target of the trial will be calculated by averaging all targets of time slices.
To confirm the parameters and weights of RNN architecture, the characteristics of non-linearity and nonstationarity in EEG signals will be considered, since the characteristics will limit the reliability of the conventional activation function in the deep learning architecture. Therefore, there are three different activation strategies, "tanh", "sigmoid" and "ReLu", for constructing activation functions [62]. The activation function, "tanh", is applied to cell input activation function of both the LSTM unit and GRU. The activation function, "sigmoid", is applied to the cell output activation function of the LSTM units.  Fig. 10 The curves of LSTM-RNN and GRU-RNN's training and validation accuracies as the number of hidden layers increased. The validation accuracies appear peaks with a 20-15 hidden layers. The results of accuracies will be paradoxical with the results of cross-entropies, and the classification of EEG signals will be quickly over-fitting. The reason is that the deep RNN architecture continues to learn common components of the EEG sequences, while simultaneously learning signal noise and non-discriminative components. Therefore, there is a trade-off between classification and long-term patterns of modeling errors. a The curves of LSTM-RNN and b The curves of GRU-RNN To prevent the vanishing error flow, the "ReLu" activation function is applied to the gates activation function of both the LSTM unit and GRU. The weights of RNN are initialized by a Gaussian distribution N ∼ (0, 0.2). The BPTT algorithm is used to train RNN by minimizing crossentropy (see (8)) loss function. Also, because the Adam strategy [63] is suitable for time-series in deep classifiers and its momentum improves the robustness of error flow, the strategy is applied to compute the learning rate during BPTT. Finally, a "Dropout" strategy is applied to prevent the over-fitting problem [64]. The key idea of the "Dropout" strategy is to randomly eliminate units (along with their connections) from the neural network during training. By experience, the dropout rate is set at 0.2, and the maximum number of iterations is set at 200.
For both "Dataset 2a" and "Dataset 2b", we train two different RNN architectures, each of which includes LSTM unit and GRU. Figure 11 illustrates the learning curves for different memory units in different datasets. In the case of both datasets, GRU-RNN architecture converges faster than LSTM-RNN architecture. To reach lowest loss, GRU-RNN architecture acquire less number of iterations than LSTM-RNN architecture. For some specific subjects, GRU-RNN architecture obtains lower average cross-entropy loss than LSTM-RNN architecture within 200 iterations. Overall, the subjects' EEG signals from "Dataset 2a" and "Dataset 2b" represent similar average cross-entropy between LSTM-RNN architecture and GRU-RNN architecture. Table 5 gives the average training and validation time complexity per trial comparison between spatialfrequency-sequential relationships and spatial-frequency features. To compare time complexity, "Dataset 2a Subject 7" and "Dataset 2b Subject 8" are used to detect the average training and validation time complexity per trial. In Table 5    All classification experimental results of all subjects in "Dataset 2a" and "Dataset 2b" are listed in Tables 6  and 7, respectively. Spatial-frequency-sequential relationships extracted from LSTM-RNN architecture and GRU-RNN architecture; spatial-frequency features extracted from CNN, SVM with linear, polynomial and RBF kernels are used to classify motor imagery for comparison. The results are presented by error rate forms, and a paired t-test statistical technique is used to detect whether the spatial-frequency-sequential relationships significantly outperform than spatial-frequency features in the classification of MI-BCIs. For each subject, we confirm the optimal number of hidden layers (ONHL) for LSTM-RNN architecture and GRU-RNN architecture.
From the results in Tables 6 and 7 To compare the classification performances of GRU-RNN and LSTM-RNN, Fig. 12 gives the classification accuracies for all subjects with all algorithms.

Subjects
Spatial-frequency-sequential relationships Spatial-frequency features

Discussion for sequential relationships
Four different sub-experiments have been created to analyze the application of sequential relationships. From smoothing time window experimental results shown in Fig. 7, a small size of smoothing time window leaded to an improvement of classification accuracy for "Dataset 2a" and a decline of classification variance for "Dataset 2b". However, a large size of smoothing time window leaded performance to decline for both "Dataset 2a" and "Dataset 2b". The smoothing time window with different sliding sizes can change the sequential relationships; in addition, the experimental results demonstrated that the relationships significantly changed the classification performance; since the classification results can be changed if different sizes of smoothing time windows were applied to the sequential relationships. The finding gives us a novel enlightenment to smooth the extracted features to improve the classification performance and robustness [8,12]. In addition, the sequential learning in the NLP also suggested to consider the sequential relationships as the key features for solving the natural language processing (NLP) problems [65,66]. Therefore, due to the EEG signals contained the similar characteristics as the sentence structure in the NLP, the finding can assist us to use the sequential relationships to model the EEG signals.
In the experiments of representing EEG signals' sequential relationships by RNN architectures, we found more hidden layers number of RNN architecture represented the EEG signals well, but more hidden layers number also caused memory vanishing problem(see Figs. 8, 9 and 10). To overcome the memory vanishing problem in the conventional RNN architecture, the LSTM-RNN architecture and GRU-RNN architecture have been introduced for the classification. We have validated the training and validation results in Figs. 9 and 10, the results of accuracies were paradoxical with the results of crossentropies, and the classification of EEG signals will quickly over-fitting. The reason is that the deep RNN architecture continues to learn common components of the EEG sequences, while simultaneously learning signal noise and non-discriminative components [47]. Hence, here we must use propriate numbers of hidden layers to retain the classification performance. From the results in Fig. 10, the number of the hidden layers of the LSTM-RNN was about 30, and the number of hidden layers of GRU-RNN was about 35 [67,68] (see Fig. 10). Therefore, the constraint of LSTM-RNN/GRU-RNN architecture leads us to crop the trial of EEG signals to time slices to feed into the classification architectures. Therefore, the SWCS is introduced on the time-series to crop the entity of a trial into several time slices. The time slices keep the same length of the number of hidden layers in order to well trained the LSTM-RNN architecture and GRU-RNN architecture.

Discussion for EEG classification
In the EEG classification experiment, he results of RNN architectures outperformed the state-of-the-art methods (see Tables 6 and 7, and Fig. 12). There are two reasons for the results. The first reason is that EEG signals usually have easy distinguishing parts and difficult distinguishing parts. "Easy parts" for classification are well represented by the spatial-frequency features, since these features are statistical features. However, the "difficult  Fig. 12 The classification accuracy performances for all subjects with all algorithms. CNN and SVM (Linear) outperforms RNN in some subjects with high-level (over 60%) accuracies (S3, S7, S8, S9 in "Dataset 2a" and S4, S6, S8, S9 in "Dataset 2b"). However, in low-level (below 60%) accuracies of subjects, RNN outperforms CNN and SVM-Linear (S2, S4, S6 in "Dataset 2a" and S2, S3 in for "Dataset 2b"). In the average-level accuracies, RNN outperforms CNN and SVM (Linear). a Classification accuracy performances for "Dataset 2a" and b Classification accuracy performances for "Dataset 2b" parts" are nonlinear and non-stationary; therefore, the statistical features cannot well model these parts [69,70]. Since the RNN architectures have enough neurons to fit the sequences' nonlinear and non-stationary characteristics, the introduced spatial-frequency-sequential relationships retain the classification performance of "easy parts", while improve the classification performance of "difficult parts". The second reason is that the conventional spatial-frequency features regarded EEG signals as a complete entity, many factors can influence classification; in particular, the subject-specific diversity might be significant one. Instead of conventional features, spatialfrequency-sequential relationships consider the EEG signals as time-series to take spatial, frequency and temporal features of EEG data into consideration. This idea not only reduces the factors that influence classification, but also increases the corresponding robustness of the subjects by reducing the subject-specific diversity. Hence, the classification accuracies of MI-BCI by spatial-frequencysequential relationships are significantly higher than those using the spatial-frequency features.
To solve the limitation of hidden layers number of RNN architecture, the SWCS is introduced on the time-series to crop the entity of a trial into several time slices. In this way, the number of samples for deep learning models is widely increased; therefore, enough samples are required to satisfy the generalization and performance of classification. In common, trials of EEG signals were obtained from complicated devices and the number of samples was too less to train the deep learning models. The application of SWCS solved the problem of sample number, and can also fit the RNN architecture. In addition, the LSTM and GRU are two different memory units for the RNN architecture. In order to test which memory unit is suited for processing EEG signals, the two memory units were applied to the RNN architectures. Our experimental results showed that GRU-RNN architecture was suited for the EEG signals better, and the spatial-frequency-sequential relationships with GRU memory units outperformed both shallow learning and deep learning models (see Table 5).

Advantages of spatial-frequency-sequential relationships
Since the EEG signals are nonlinear and non-stationary signals, the spatial and frequency features are the robust statistical features, which can be well classified by SVM using FB-CSP features. We considered such features as the "easy parts" for classification; in contrast, we considered the temporal features as the "difficult parts" that are difficult to be classified by the conventional machine learning models. Hence, we introduced the spatialfrequency-sequential relationships by using RNN architecture on FB-CSP features. Experimental results showed that the approach involving spatial-frequency-sequential relationships can have the better classification performance of the "difficult parts"; also, the results outperformed the-state-of-the-art methods. Besides classifying the "difficult parts" of EEG signals, another reason for introducing the RNN architecture was that: the limitation number of hidden layers leaded us to crop the entity of a trial into several time slices. Although the different time slices shared the same labels, the slices increased the diversity of EEG signals during the classification. Therefore, the cropping strategy also improved the classification performance. To sum up, the advantages of introducing spatial-frequency-sequential relationships can improve classification performance and increase EEG samples diversity.

Conclusion
In this paper, an FB-CSPs algorithm was used to extract spatial-frequency features, which were cropped by a sliding window cropping strategy into time slices. Then, the time slices were fed into deep RNN architectures, with two different memory units, to extract spatial-frequencysequential relationships for MI-BCI classification. The extracted relationships included spatial, frequency and temporal characteristics. The experiments on MI-BCI demonstrated that the proposed method owned two advantages: (1) The spatial-frequency-sequential relationships extracted by FB-CSPs and RNN architectures can achieve significantly higher performance than spatialfrequency features. Meanwhile, the relationships had the same level of time complexity with the conventional algorithms. (2) A comparison of the accuracy and efficiency of motor imagery classification between GRUs and LSTM units revealed that GRUs can generate the better results.
Our future work will focus on collecting more EEG data to construct deeper RNNs architecture, exploring the error flow rules in BPTT, and constructing a deeper RNNs architecture that is adapted and generalized for EEG signals.