Improvements in Classification of Left and Right Foot Motor Intention Using Modulated Steady-State Somatosensory Evoked Potential Induced by Electrical Stimulation and Motor Imagery

In recent years, motor imagery-based brain-computer interface (MI-BCI) has been applied to motor rehabilitation in patients with motor dysfunction. However, traditional MI-BCI is rarely used for foot motor intention recognition because the motor cortex regions of both feet are anatomically close to each other, and traditional event-related desynchronization (ERD) patterns for MI-BCI have insufficient spatial discrimination. This study introduced steady-state somatosensory evoked potentials (SSSEPs) by synchronous bilateral feet electrical stimulation at different frequencies, which were used as carrier signals modulated by unilateral foot motor intention. Fifteen subjects participated in MI and MI-SSSEP tasks. A Riemannian geometry classifier with a task-related component analysis (TRCA) spatial filter was proposed to demodulate the variation in SSSEP features and discriminate the left and right foot motor intentions. The feature outcomes indicated that the amplitude and phase synchronization of the SSSEPs could be well modulated by unilateral foot MI tasks under the MI-SSSEP paradigm. The classification results revealed that the modulated SSSEP features played an important role in the recognition of left-right foot discrimination. Moreover, the proposed TRCA-based method outperformed the other three methods and improved the foot average classification accuracy to 81.07± 13.29%, with the highest accuracy attained at 97.00%. Compared with the traditional MI paradigm, the foot motor intention recognition accuracy of the MI-SSSEP paradigm was significantly improved, from nearly 60% to more than 80%. This work provides a practical method for left-right foot motor intention recognition and expands the application of MI-BCI in the field of lower-extremity motor function rehabilitation.


Improvements in Classification of Left and Right
Foot Motor Intention Using Modulated Steady-State Somatosensory Evoked Potential Induced by Electrical Stimulation and Motor Imagery Yan Bian , Li Zhao, Senior Member, IEEE, Jiaying Li, Tong Guo , Xing Fu, and Hongzhi Qi , Member, IEEE Abstract -In recent years, motor imagery-based braincomputer interface (MI-BCI) has been applied to motor rehabilitation in patients with motor dysfunction. However, traditional MI-BCI is rarely used for foot motor intention recognition because the motor cortex regions of both feet are anatomically close to each other, and traditional event-related desynchronization (ERD) patterns for MI-BCI have insufficient spatial discrimination. This study introduced steady-state somatosensory evoked potentials (SSSEPs) by synchronous bilateral feet electrical stimulation at different frequencies, which were used as carrier signals modulated by unilateral foot motor intention. Fifteen subjects participated in MI and MI-SSSEP tasks. A Riemannian geometry classifier with a task-related component analysis (TRCA) spatial filter was proposed to demodulate the variation in SSSEP features and discriminate the left and right foot motor intentions. The feature outcomes indicated that the amplitude and phase synchronization of the SSSEPs could be well modulated by unilateral foot MI tasks under the MI-SSSEP paradigm. The classification results revealed that the modulated SSSEP features played an important role in the recognition of left-right foot discrimination. Moreover, the proposed TRCA-based method outperformed the other three methods and improved the foot average classification accuracy to 81.07±13.29%, with the highest accuracy attained at 97.00%. Compared with the traditional MI paradigm, the foot motor intention recognition accuracy of the MI-SSSEP paradigm was significantly improved, from nearly 60% to more than 80%. This work provides a practical method for left-right foot Index Terms-Brain-computer interface, motor imagery, foot motor intention recognition, steady-state somatosensory evoked potentials, task-related component analysis.

I. INTRODUCTION
B RAIN-COMPUTER interface (BCI) can replace, repair, enhance, supplement, or improve the normal output of the central nervous system and change its interaction with the internal/external environment [1]. In recent years, BCIs have been applied to the plasticity rehabilitation of limb motor ability in patients with amyotrophic lateral sclerosis, spinal cord injury, stroke, or other diagnoses, as they can promote the recovery of cerebral function or behavior by manipulating neurophysiological activities and represent a potential rehabilitation method [2], [3]. For example, in a lower-limb gait control system, BCIs can be exploited to effectively manipulate rehabilitation devices, such as foot nerve prostheses or gait orthoses [4], even for chronic stroke patients [5]. For motor recovery, whether the patient's motor intention can be correctly decoded is crucial to promote BCI performance.
Motor imagery (MI) refers to an inner cerebral process of motor behavior but without any actual body movement; therefore, MI-BCI is often used to detect patients' motor intentions and has been widely used in limb motor rehabilitation. In MI-BCI, an increase or decrease in electroencephalogram (EEG) energy occurs in alpha or beta frequency bands, mainly over the primary sensorimotor area during task periods, namely, event-related synchronization (ERS) and event-related desynchronization (ERD) [6], [7], which are vital features to identify target MI and recognize corresponding motor intentions. However, in classical MI-BCI, the types of imagery paradigms are very limited [8], [9], especially for lowerlimb MI, where both feet are usually treated as one imagery paradigm without unilateral discrimination. The main reason is that in primary sensorimotor cortices, the foot area locates near the interhemispheric fissure [10], [11], and thus, the left and right foot areas from the two cerebral hemispheres are This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ anatomically close to each other. Due to low spatial resolution of EEG recordings, the left and right foot MI produce nearly identical EEG patterns [11], [12], making the classification of left and right foot motor intention difficult. Therefore, few reports have focused on lower-limb motor intention discrimination. According to the available documents, the possible EEG features for foot unilateral classification are primarily mu-ERD, beta-ERD, and beta-ERS in traditional MI-BCI. However, these previous results seem controversial; one work reports an ERD/ERS left-right difference for all three features during a foot MI task [12]. In contrast, another similar study believed that mu-ERD and beta-ERD were not prominently lateralized for foot motor intention discrimination [11], and unlike the above two literatures, another study demonstrated both mu-ERD and beta-ERS exhibited no significant difference between left and right foot MI [13]. These conclusions indicate that ERD/ERS left-right foot discrimination may be unreliable and insufficient, which may be associated with MI paradigm, reference method, subject population, and so on. Moreover, the classification results for ERD/ERS left-right foot discrimination remain unsatisfactory, that is, the average accuracies of most studies are approximately 70%, which is unsuitable for practical applications. For example, for left-right foot discrimination, Hashimoto et al. reported an average classification accuracy of 69.3% [11], and Tariq et al. reported an average accuracy of 70.28% [14]. Recently, some studies have attempted to develop classification algorithms, such as common spatial patterns and their extensions, support vector machine (SVM)-based algorithms, and the k-nearest neighbor (KNN) classifier [10], [12], [14], to aid in enhancing leftright foot discrimination, but the small improvements could not facilitate the discrimination to a practical level. Therefore, further investigations are still essential to develop reliable features and methods for higher foot classification accuracy.
To find a reliable and discriminable feature that can lead to high classification accuracy for left-right foot discrimination, this paper introduces a different feature unlike the common ones mentioned in previous studies for foot motor intention recognition, namely, modulated steady-state somatosensory evoked potential (SSSEP), which is induced by electrical stimulation and MI tasks. SSSEP, as a type of steady-state somatosensory characteristic potential, has been introduced into hybrid BCIs in recent years [15]. As it has been proven that the combination of SSSEP features and MI tasks can enhance the spatial resolution of MI recognition and improve the recognition accuracy of upper-limb classification performance [9], [16], [17], it is natural to assume the good performance of introducing SSSEP into MI-BCI for foot motor intention recognition. The electrical stimulation could produce a stable SSSEP of the same stimulation frequency and activate the cerebral sensory cortex. As it only relies on the human sensory system, somatosensory stimulation maintains the advantages of a traditional MI-BCI, as compared with visual or auditory stimulation. Furthermore, SSSEP could be modulated by left or right foot motor intention; the execution of a unilateral foot MI task would suppress the amplitude and phase synchronization of the SSSEP induced on the same site; consequently, the left and right foot motor intentions can be loaded on the modulated SSSEP. Therefore, the dis-crimination problem of spontaneous EEG signals with a low signal-to-noise ratio can be transformed into the demodulation of a reliable SSSEP feature with a high signal-to-noise ratio.
In addition to introducing a new feature for left-right foot discrimination, a suitable demodulation method is also proposed according to the lock-in time and phase properties of the SSSEP, namely, a Riemannian geometry classifier with a task-related component analysis (TRCA) spatial filter, which is sensitive to SSSEP variations. In addition, we also compare the proposed methods with other typical feature-extraction and classification algorithms, such as a filter bank common spatial pattern (FBCSP) algorithm and a SVM classifier, which is prevailing for MI-BCI, to analyze the classification performance of different paradigms, features, and methods. This work not only provides a practical method for left-right foot motor intention recognition, but also enriches the types of limited imagery paradigms and expands application of MI-BCI in the field of lower extremity motor function rehabilitation, such as in the gait control system.

A. Participants
Fifteen healthy subjects (seven male and eight female), with an average age of 23.7±1.9 years, were recruited in this experiment. They were all college students with righthandedness, normal or corrected normal vision and sensation, and no history of trauma or neurological-related diseases. Before the formal experiment, each subject was trained for three days (half an hour each day) to familiarize themselves with electrical stimulation, MI tasks, and specific experimental processes as quickly as possible. Before the experiment, all subjects were informed of the experimental process and provided signed informed. This study was approved by the ethics committee of Tianjin University.

B. Electrical Stimulation
In our study, two-channel current pulses with a pulse width of 0.2 ms were used to induce the SSSEP. Two-channel electrodes were placed at the posterior tibial nerve of the medial malleolus of both the left and right foot. The stimulation position was approximately 2-3 cm behind the medial malleolus of the feet. Before the beginning of the formal experiment, the exact stimulation position and current stimulation intensity of each subject was adjusted to achieve the effect of initial toe contraction without any pain. The lower-limb electrical stimulation intensity of different subjects is generally between 10 and 25 mA, and there are individual differences in the optimal stimulation frequency among different subject [4].
To determine the general optimal stimulation frequency and avoid overlapping with the frequency bands of MI (i.e., alpha and beta bands), the stimulation frequency of the left foot usually ranges from 26 to 30 Hz, and that of the right foot ranges from 31 to 35 Hz. The preliminary experimental results showed that all stimulation frequencies mentioned above could induce a stable SSSEP, but the amplitude of the SSSEP was the most obvious at the 28/33 Hz combination. Therefore, we selected 28 Hz as the left-foot stimulation frequency and 33 Hz as the right-foot stimulation frequency for all subjects. In addition, the two stimulation channels of the left and right foot were executed simultaneously. Figure 1 shows the device, the stimulation positions of both feet and the stimulation results of the SSSEP signal induced by electrical stimulation in channel Cz of a typical subject. The peak frequencies of somatosensory characteristics were consistent with the electrical stimulation frequencies used in our study, which indicated that the SSSEP in the primary sensorimotor area could be effectively induced by using the above electrical stimulation parameters.

C. Experiment Description
The two paradigms in this study, MI and MI-SSSEP, are explained as follows: a) MI is a classical motor imagery paradigm. Under this experimental condition, the subjects needed to imagine the unilateral kinesthetic MI task of the left foot or right foot hook action according to a text prompt of "left foot" or "right foot" on a screen. The imagery task lasted for 5 s. b) MI-SSSEP is a hybrid experimental paradigm, and the subjects needed to carry out the same kinesthetic MI task as in the MI paradigm under the condition of electrical stimulation. The duration of the task was 5 s. One second before the task, electrical stimulation started until the end of the MI task. When electrical stimulation was applied to the left and right foot, the prompts of " --Prepare" or "Prepare --" appeared on the screen, and the subjects needed to prepare for the MI-SSSEP task. One second later, the prompt of "left foot" or "right foot" was displayed with a red "+" indicating that the unilateral MI task should be done under somatosensory stimulation during this period.
The time course of a single trial is shown in Figure 2. Data from 80 trials were collected under each paradigm, including 40 trials of the left foot and 40 trials of the right foot, which appeared in a random order. Each session of the experiment contained 40 trials, and after each session, there was a rest of 3-5 minutes to prevent fatigue and keep the subjects in a good experimental state. A SynAmps2 electrophysiological amplifier (Compumedics Corporation, Australia) and a 64-channel electrode cap (using the international standard 10-20 system) were used for EEG signal acquisition. Reference and ground electrodes were set to the left and right mastoids, respectively. Scan 4.5 software was used for EEG signal acquisition, display, and storage. The data preprocessing included removing useless electrodes (M1, M2, CB1, and CB2) and performing a common average reference. Subsequently, the data were filtered using a 5-100 Hz bandpass filter to eliminate irrelevant frequency bands. Afterwards, blind source separation technology was used to remove artifacts, such as electrooculogram and electromyography. Finally, the data were downsampled to 200 Hz for subsequent analysis.

D. ERD and SSSEP Feature Analysis
To further study the neural activities in the sensorimotor cortex under different paradigms, the event-related spectral perturbation (ERSP) was calculated based on a short-term Fourier transform. The disturbance of the event-related spectrum of the task periods relative to the baseline time period was calculated [7], [18]. Typically, F h ( f, t) is the spectral estimate of trial h at frequency f and time t for n trials, the ERSP can be calculated as follows: In this study, the Hanning window was used for the shortterm Fourier transform. For the experimental paradigms, the average ERSP results of the Cz channel were calculated.
To further explore the modulated SSSEP features (27-29 Hz and 32-34 Hz) of the MI-SSSEP paradigm, a brain topographic map of the left/right foot task was constructed [19]. A grand average brain topographic map of 60 leads was used to compare the differences in the spatial distribution characteristics of the SSSEP spectra.
This study further analyzed whether there were significant differences between the ERD or SSSEP features of left-and right-lateral tasks under the two paradigms. The spectrum values of each subject during the left/right foot tasks in the alpha, beta, or modulated SSSEP frequency bands were calculated on the Cz channel, and then the comparisons of these values were implemented between the left and right foot tasks. As a modulation signal influenced by MI, the phase synchronization of the SSSEP may also change during the MI task periods. Here, we employed interstimulus phase coherence (ISPC) to measure phase synchronization variations under somatosensory stimulation [17]. Assuming that f is the SSSEP frequency, the original EEG S(t) is first filtered to the [ f − 0.5, f + 0.5] SSSEP band, and then the filtered signal S f (t) is segmented by a sliding window, which moves 1/ f seconds between two adjacent sliding windows. The segmented signals are defined as S f k (t), where k denotes the k-th segment. IfŜ f k (t) represents the Hilbert transform of S f k (t), then the k-th segment analysis signal X f k (t) is expressed as Formula (2): ISPC can be defined as follows: where t denotes the time interval 1/ f seconds. The ISPC is the result of averaging N s stimuli, where N s is the number of stimuli contained in one second. The ISPC values are between 0 and 1. A value of 0 indicates that there is absolutely no phase synchronization, and a value of 1 indicates perfect phase synchronization.
To further study the separability of foot motor intention in the SSSEP, we calculated the grand average separability r 2 coefficient during the execution of the left/right foot imagery task under the MI-SSSEP paradigm. The calculation process is shown in Formula (4) [17]. The sample signal,x pj _1 , of the SSSEP band was composed of several trials in a 28-or 33-Hz frequency band during the left foot imagery task on the j -th specific channel, where N pj _1 is the sample length. x pj _2 is a sample signal composed of multiple trials during the right foot imagery task on the j -th channel of the same subject, N pj _2 is the sample length, and P is the total number of subjects.
E. Classification 1) Riemannian Geometry Classifier With TRCA Spatial Filter: a) TRCA spatial filter: The TRCA algorithm could exhibit maximal similarity among task blocks and extract task-related components by properly weighing the observed signals [20]. Assuming task-related signal s(t) ∈ R and task-unrelated signal n(t) ∈ R, the observed multichannel signals x(t) ∈ R N c are composed of the following (5): x j (t) = a 1, j s(t) + a 2, j n(t), j = 1, 2, . . . , N c (5) where j denotes the channel index with N c is the number of channels, and a 1, j , a 2, j are the mixing coefficients of the j -th channel. Then, the task-related components, s(t), can be recovered from a linear sum of the observed signal, x(t). To ensure that the task-related components exhibit maximal temporal similarity among task trials, inter-trial covariance maximization should be converted. If x (h) (t)(h = 1, 2, · · · N t ) is the h-th trial of the observed signal, and y(t) represents the estimated task-related component, then the covariance of all possible combinations of trials is summed as Formula (6): where the symmetric matrix S is defined as follows: To bound the coefficients, the variance of y(t) is restricted as follows: The restricted optimization can be solved using Formula (9): The optimal coefficient vector is obtained as the eigenvector of the matrix Q −1 S. The eigenvalues of matrix Q −1 S, which can be arranged in descending order, imply task consistency among the experimental trials [21]. In our study, we used the TRCA method to design spatial filters during task interval (0-3 s, where 0 s represented the start of the MI task) to obtain the SSSEP task-related components. In addition, this method can reduce the feature dimensions for Riemannian geometry measurement classification. b) Data partitioning, multi-trial reference, and feature extraction: Considering there should be no overlap between the training set and the test set for constructing the spatial filter and task classification, the MI-SSSEP task interval signals were first divided into a training set and a test set according to the 10 × 10-fold cross validation method. Subsequently, the training set was filtered by 27-29 and 32-34 Hz subbands, and then the filtered multiple training trial data X 28 train , X 33 train were used to build TRCA spatial filters and attain w 28 , w 33 . The first two eigenvectors were employed to obtain the averaged task-related component, y re f , which was used as a reference after averaging for multiple training trials. Afterwards, the training and test trials were filtered by TRCA filters separately, and the task-related components y train and y test were obtained. The covariance matrices of y re f y train and y re f y test were estimated, which were used as feature vectors for the Riemannian geometry measurement. The above processes of data partitioning, TRCA filters constructing and feature extracting were visualized in Figure 3. c) Riemannian geometry measurement: The Riemannian distance can be used to classify left-right foot motor intentions. Here, given two samples, x i and x j , the Riemannian distance between them can be computed, as shown in Formula (10) [22]: where i represents the covariance matrix of x i , j denotes the covariance matrix of x j , λ c (c = 1, . . . , N c ) are the real eigenvalues of , and N c is he number of channels. The Riemannian geometric mean G(.) of M covariance matrices is defined by the matrix that minimizes the sum of the squared Riemannian distances, as shown in Formula (11): d) Classification: Given a training set and the corresponding labels, r i {1, 2} (1: left-foot motor intention; 2: right-foot motor intention), the training of the classifier estimates the Riemannian geometric mean for each class as follows (12): where K {1, 2} denotes the class label. The classification of a test trial can be obtained by computing the minimum Riemannian distance to each class, as shown in Formula (13): wherer denotes the predictive class. The classification process was also displayed in Figure 3. Figure 3 showed overall diagrams of the proposed method.
2) Baseline Classification Algorithm: An FBCSP based on filter bank strategy is selected as a baseline feature extraction algorithm under both MI and MI-SSSEP paradigms, as it is a prevailing method in BCI applications [14], [23]. According to the FBCSP algorithm process, the training-set data were filtered by multiple sub-bands, which mainly included alpha and beta bands, as well as the first-and second-harmonic bands of the SSSEP stimulation frequency. The specific frequency bands were occupied at 8-13, 13-26, 27-29, 32-34, 55-57, and 65-67 Hz. These six sub-bands were then selected to construct a spatial filter. After bandpass filtering of each subband in the training process, the FBCSP algorithm was used to build a spatial filter model, and the first two maximum and minimum eigenvalues were selected to form a spatial filter in our study. Afterwards, the spatial filter model was used in the bandpass filtered data of the test set to obtain the feature vectors of each sub-band. After combining the features of all sub-bands, a SVM classifier was used to obtain the accuracies for the two paradigms as a baseline classification method, which is also a classical method in BCI pattern recognition [24], [25]. To ensure the reliability of the classification results, a 10 × 10-fold cross validation was performed, and only training data were involved in forming a spatial filter or SVM classifier model during the process.

F. Statistical Analysis
The Wilcoxon signed-rank test was employed to analyze the modulation effect in the ERD or SSSEP features by unilateral foot tasks. As a kind of non-parametric test, the Wilcoxon signed-rank test is not sensitive to data distribution since not all ERD or SSSEP features satisfied the normal distribution when using Shapiro-Wilk test before the comparisons. Moreover, the results of decoding accuracy recognition among the two experimental paradigms were also analyzed by the Wilcoxon signed-rank test based on MATLAB software. In addition, the classification performances of the MI-SSSEP were evaluated using three different features (i.e., ERD/SSSEP/hybrid features) and the one-way repeated measurement ANOVA with Bonferroni post-hoc test was conducted to obtain the multi-dataset statistical results. The Shapiro-Wilk test was also fulfilled before ANOVA tests to verify normality. Additionally, to explore the optimal classification methods for the modulated SSSEP features, statistical comparisons of four methods were further investigated. The ANOVA tests were performed with IBM SPSS 20 in our study.

A. ERD and SSSEP Features
The characteristic variations in the grand average ERSP with time and frequency under the two experimental paradigms are shown in Figure 4(a). This figure shows that the energies of the alpha and beta bands decreased during the imagery task periods under the MI and MI-SSSEP paradigms, which indicated that ERD phenomena occurred in both paradigms. In addition, under the MI-SSSEP paradigm, SSSEP features (c) r 2 coefficient mean topographical maps of 28/33 Hz in the MI-SSSEP. " * " means P < 0.05, and " * * " means P < 0.01. occurred at 28 and 33 Hz after the beginning of electrical stimulation and lasted until the end of electrical stimulation, and the features induced on the Cz channel were extremely obvious. This phenomenon was consistent with the conclusion that the left and right foot cerebral mapping areas were located near the interhemispheric fissure. The brain topographic map distribution results of grand average SSSEP features at all 60 electrodes, are also shown in Figure 4(a). The spatial distribution results under the MI-SSSEP paradigm revealed that the areas of somatosensory characteristics at 28 and 33 Hz were mainly distributed near the midline of the central area. Moreover, when performing different unilateral MI tasks, the spatial distributions of SSSEP features at different stimulation frequencies under the MI-SSSEP paradigm also appeared to vary, indicating that the unilateral MI tasks could modulate the SSSEP signals well. From the temporal-spatial-frequency analysis results of EEG under the different experimental paradigms, the MI-SSSEP paradigm exhibited an EEG pattern that included both ERD and SSSEP features. However, the features that play a key role in the improvement of BCI classification performance in left/right foot motor intention discrimination still needed to be clarified. To select separable features from the hybrid paradigm, we further analyzed the ERD and SSSEP features on the Cz channel for the two experimental tasks. Figure 4(b) shows the grand average alpha and beta ERD features evaluated from 0.5-4.5 s across all subjects (0 s represented the beginning of the task) in the MI and MI-SSSEP paradigms. The results of the Wilcoxon signedrank test showed that there was no significant difference in ERD characteristics between the left and right foot tasks in either the MI or MI-SSSEP paradigm (MI paradigm alpha: P = 0.5000 >0.05, beta: P = 0.4102 > 0.05; MI-SSSEP paradigm alpha: P = 0.1250 > 0.05, beta: P = 0.4551 > 0.05). Therefore, ERD features alone did not provide a sufficient separation effect on foot movement intention recognition. The Wilcoxon signed-rank test was also used to analyze the grand averaged SSSEP features on the Cz channel under the MI-SSSEP paradigm. The results showed that there were statistically significant differences between the average power spectrum values of the left and right foot MI tasks at the two sensory stimulation frequencies under the MI-SSSEP paradigm (SSSEP 28 Hz: P = 0.0075 < 0.01, SSSEP 33 Hz: P = 0.0416 < 0.05), as shown in Figure 4(b). These results implied that the SSSEP could be well modulated by the MI tasks. Moreover, the SSSEP variations were well separable as motor intention carrier signals under the MI-SSSEP paradigm. Subsequently, we further computed the grand average r 2 coefficient topographic map of 60 channels during the foot imagery tasks under the MI-SSSEP paradigm, as shown in Figure 4(c). The results demonstrated that the separability features mainly distributed in the sensorimotor area of the central area, which agreed with the modulated SSSEP topographic maps. Therefore, it was inferred that the SSSEP features of the sensorimotor cortex could provide good separability for discriminating foot motor intention.
The MI tasks not only had a significant impact on SSSEP amplitudes, but also changed the phase synchronization in the MI-SSSEP, as shown in Figure 5(a). The grand average ISPC curves on the Cz channel demonstrated that the phase synchronization could be improved rapidly after the beginning of electrical stimulation, and that ISPC values decreased approximately 0.5 s ahead of the end of somatosensory stimulation. Therefore, the electrical simulation caused fine phase synchronization at the same stimulation frequencies. In particular, the unilateral MI task obstructed the phase synchronization of the ipsilateral stimulation frequency; for example, when implementing the left foot MI task, the ISPC value of 28 Hz decreased dramatically from the start of the MI task because the left-foot stimulation frequency was 28 Hz, but while fulfilling the right-foot MI task, the phase synchronization of 28 Hz was not influenced during the MI task period, and vice versa. To further study the ISPC feature discrepancy between two unilateral foot MI tasks, the r 2 coefficient and corresponding topographical maps at different times were compared, as seen in Figure 5(b). As shown in the lower panel of Figure 5, the ISPC revealed a high discrepancy during the MI task period, which agreed with the conclusion of upper-limb [17], but the results of lower-limb were much smaller. The topographical maps of the r 2 coefficient spatial distributions showed that the ISPC discrepancies between left-right foot motor intentions mainly focused on the sensorimotor area near the central area, which was consistent with the previous results shown in Figure 4(c). These findings verified that both the amplitude and phase synchronization of the modulated SSSEP under the MI-SSSEP paradigm were closely related to unilateral foot MI tasks and different motor intentions.

B. Classification Performance
Whether the EEG signals in the SSSEP bands exhibited elevated separability after TRCA spatial filtering was vital to determine the final classification performance. Therefore, the separability of the EEG signals before and after the TRCA w-28/w-33 spatial filters was further analyzed. The TRCA-based filtered results are shown in Figure 6. The left panel in Figure 6(a) shows the r 2 time-frequency maps of data filtered by the w-28/w-33 TRCA spatial filter of a typical subject, which revealed that the filtered data still maintained good separability after the TRCA filtering, especially during the 0-3 s after the start of the MI task. The right panel in Figure 6(b) displays the statistical comparisons of the grand average r 2 coefficients before and after the TRCA filter. The two violin plots indicated that TRCA filtering by w-28 and w-33 could both significantly improve the SSSEP discrepancies (P <0.001); thus, as a carrier signal modulated by unilateral MI tasks, such a substantial promotion was determined to improve foot motion intention recognition based on the TRCA spatial filter method. Comparisons of recognition performance for different features and methods Fig. 6. TRCA filter analysis, including (a) r 2 temporal-frequency maps of data filtered by the SSSEP 28/33 Hz TRCA spatial filter for a typical subject (S4), and (b) grand average r 2 coefficient comparisons before and after the TRCA filter. " * * * " means P < 0.001. are exhibited in Figure 7. As shown in Figure 7(a), "HY" and "MI" represent the results of the MI-SSSEP and MI paradigms, respectively; "FC," "FS," and "FE" denote the baseline method classification using both ERD and SSSEP features (six sub-band features), SSSEP features (four subband features with first and second harmonic), and ERD features (alpha and beta features), respectively. These results revealed that there was a significant classification accuracy difference using the baseline method among different features for the MI-SSSEP paradigm (one-way repeated measurement ANOVA: F (1.062,14.867) = 15.297, P = 0.001). The accuracies of applying hybrid and SSSEP features both significantly outperformed only using ERD features, with P values of 0.001 and 0.014, respectively, while there was no significant classification difference between employing hybrid and SSSEP features (P=0.155). There was no obvious accuracy difference with different features for the MI paradigm (hybrid vs. ERD features). In addition, the Wilcoxon signed-rank test results indicated that the classification accuracy of the MI-SSSEP paradigm was significantly higher than that of the MI paradigm using hybrid features (P = 0.0013), with the average accuracy of 15 subjects able to be improved from 59.57% to 71.22%. However, the statistical results show that there was no significant difference in classification performance between the MI-SSSEP and MI paradigms (P = 0.4670 > 0.05) when both applied ERD features. To further explore the optimal classification method for the modulated SSSEP features, the classification outcomes are recalculated with four different methods, including the proposed TRCA-based Riemannian geometry method (TRCARIE), TRCA-based SVM method (TRCASVM), FBCSP-based Riemannian geometry method (FBCSPRIE), and the baseline method of FBCSP-based SVM classifier (FBCSPSVM). Figure 7(b) shows the classification results of these methods using two basebands of SSSEP features under the MI-SSSEP paradigm and ERD features under the MI paradigm. The best average accuracy of 15 subjects further reached 81.07% for the MI-SSSEP paradigm, with the maximum individual accuracy reaching 97.00%, as obtained by the proposed TRCARIE method. Moreover, the classification accuracy of the TRCARIE method significantly outperformed the other methods for the MI-SSSEP paradigm, namely, the TRCARIE versus baseline method, with the average accuracy increased by 16.41% (P=0.002). These results showed that for the MI-SSSEP paradigm, the TRCA filter with the Riemann classifier could effectively recognize the left-right foot motor intention and was thus especially suitable for the identification of SSSEP carrier signals. Furthermore, a two-way repeated measurement ANOVA confirmed the main effects of the feature method on the accuracy were statistically significant (F (1,14) = 13.185, P= 0.003 < 0.01); the same was the case for the classification method (F (1,14) = 20.718, P < 0.001), and there was no significant interaction between them. Unlike the results of the MI-SSSEP, the average accuracies of the MI paradigm indicated that there were no significant differences among the four methods, which may be caused by the poor separability of the two ERD bands.

IV. DISCUSSION
The recognition of left and right foot motor intention is an onerous problem because the motor cortex regions of both feet are anatomically close to each other. Hence, the conventional MI paradigm has difficulty in solving left-right foot discrimination problems, which is mainly due to the lack of reliable and discriminable features. As shown in Figure 7(a), the classification accuracy using only ERD features was less than 70%, which was consistent with previous studies [11]. Figure 7(a) also indicated that the classification accuracy using modulated SSSEP features significantly outperformed that using ERD features, obtaining a similar classification efficiency as using both features under the MI-SSSEP task. These results implied that using SSSEP features was better than ERD features in 2-class foot motor intention recognition in our study. Figures 4 and 5 proved our hypothesis that SSSEP features could be modulated by the ipsilateral MI task, which influenced both the amplitude and phase synchronization, resulting in a new reliable feature for foot leftright discrimination. The modulation effect may be due to the execution of the MI task sharing some sensory neuronal clusters with the sensory resources used to generate SSSEP. Therefore, the amplitude and phase synchronization of SSSEP were suppressed during the ipsilateral MI process as a result of resources competition, while the MI activity had little effect on the SSSEP produced at the contralateral site. Obviously, this top-down neural modulation process was related to the action-perception coupling and interaction of the neuronal populations underlying the primary sensorimotor area and will bring new ideas to the development of BCI paradigms. More importantly, the highest average classification accuracy (15 subjects), obtained by using two bands of modulated SSSEP features, reached 81.07±13.29% under the MI-SSSEP paradigm, an improvement of nearly 20% compared with the average accuracy only using ERD features, which indicated more notable and more likely improvements for foot unilateral classification. To our knowledge, this study is the first to introduce modulated SSSEP as a new feature to discriminate left-right foot motor intention and significantly improves the average classification accuracy to a practical level; thus, the current work provides a novel way of lower-extremity motor function rehabilitation.
Our study used two paradigms: the proposed MI-SSSEP paradigm and the traditional MI paradigm. The experimental results suggested that the introduction of modulated SSSEP features could significantly improve the classification performance of the MI-SSSEP paradigm compared with the traditional MI paradigm when performing the identical MI task. Furthermore, based on the individual classification of 15 subjects under the MI-SSSEP paradigm, 73.3% of subjects achieved more than 70% accuracy using the proposed method, and 26.7% of subjects achieved more than 90%. Due to the obvious individual differences and low accuracy in MI-BCI, the MI-SSSEP paradigm is expected to not only provide feasible foot motor intention discrimination but also supply a practical way to improve the recognition efficiency for most subjects. In our previous research, we tried to apply electrical stimulation to an upper-limb MI task, which was used to recognize the motor intention between the left and right hands or to distinguish the motor intention between the right hand and ipsilateral arms. Compared with the single MI task, the classification accuracy was improved by more than 10% [9], [16]. The current work verified the effectiveness of using SSSEP features in lower-limb motor intention recognition, and it can be supposed that the SSSEP features generated during the MI task could help to enhance the spatial resolution of MI recognition and extend the numbers of limb imagery actions in MI-BCI, such as left/ right foot or right foot/right knee MI tasks. In addition, since the production of SSSEP depends on the human body's own sensory system, the coactivation mode of the sensory cortex and motor cortex is closer to the oscillation mode of traditional MI-BCI, as compared with other types of hybrid paradigms.
To obtain optimal results for foot motor intention recognition, the demodulation method is worthy of intensive study. Figure 7(b) displays the classification performance of the four methods used in this study, which demonstrate that the proposed Riemannian geometry with the TRCA spatial filter method was best for left-right foot discrimination under the MI-SSSEP paradigm. Although there is no report involving the TRCA method used for SSSEP feature extraction, several works reported that the TRCA method showed good performance in extracting task-related components, such as steady-state visual evoked potential (SSVEP) [20], [21]. In terms of the time-locked characteristics of the SSSEP signal compared to that of SSVEP, it is natural to assume that the TRCA spatial filter has good performance in SSSEP feature extraction. As shown in Figure 6(b), the SSSEP separability after TRCA spatial filtering was further significantly improved compared with that using pre-filtered signals, and the feature extraction methods had a significant effect on improving the classification accuracy (TRCA vs. FBCSP methods). Moreover, the statistical results demonstrated the Riemannian geometry classifier could also obtain higher classification performance for the MI-SSSEP paradigm than that of the SVM classifier. The Riemannian distance measurement can take the advantages of both temporal information and spatial information simultaneously, rather than using single information as traditional measurements; thus, this method could facilitate the classification of the modulated SSSEP, which is characterized in both time and frequency domains. Therefore, the proposed method was especially adapted for the MI-SSSEP paradigm as an ideal demodulation algorithm for the modulated SSSEP carrier signals. However, for the MI paradigm, the proposed method did not exhibit good performance. Possible explanations are that the separability of ERD feature was not strong in our study, and the ERD features did not have as good of time-locked and phase-locked characteristics compared to SSSEP features; consequently, the TRCA method may be unsuitable for processing ERD features. Recently, some new methods were proposed to extract features in MI-BCI [26], [27], and the effects for extracting ERD features using these methods are worth further investigation.
Considering the important applications of MI-BCI in the field of motor rehabilitation, the method proposed in this paper not only successfully elicited modulated SSSEP features but also significantly improved left-right foot classification accuracy by using the modulated features. Of course, the introduction of SSSEP into MI-BCI also brings some negative aspects, such as increasing system complexity and being unsuitable for somatosensory pathway impaired patients. However, considering there are many patients whose somatosensory system is intact, this study could provide a new feature, method, and experimental data support to control rehabilitation devices of lower-limb motor function for such a population.

V. CONCLUSION
In this study, based on somatosensory stimulation and MI-BCI, a new feature and method for foot unilateral motor intention recognition was proposed. SSSEP signals induced by electrical stimulation were used as carrier channels, and the unilateral foot motor intention could change the amplitude and phase synchronization of the ipsilateral SSSEP. The Riemannian geometry with a TRCA spatial filter method was introduced in this study to demodulate the variation in SSSEP carrier signals. The average recognition accuracy was improved by 16.41% using the proposed method compared with the baseline algorithms under the MI-SSSEP paradigm, reaching 81.07%, which far exceeded the 70% BCI application level. Thus, the modulated SSSEP features and proposed algorithms could effectively distinguish the foot motor intention of anatomically adjacent areas and help to expand the application of MI-BCI in the field of lower-extremity motor function rehabilitation.