Analysis of the EEG Rhythms Based on the Empirical Mode Decomposition During Motor Imagery When Using a Lower-Limb Exoskeleton. A Case Study

The use of brain-machine interfaces in combination with robotic exoskeletons is usually based on the analysis of the changes in power that some brain rhythms experience during a motion event. However, this variation in power is frequently obtained through frequency filtering and power estimation using the Fourier analysis. This paper explores the decomposition of the brain rhythms based on the Empirical Mode Decomposition, as an alternative for the analysis of electroencephalographic (EEG) signals, due to its adaptive capability to the local oscillations of the data, showcasing it as a viable tool for future BMI algorithms based on motor related events.


INTRODUCTION
A brain-machine interface (BMI) (Rao, 2013) is a system that allows controlling a device through the analysis of the electric biosignals that can be acquired from the brain, with the help of scalp electrodes. Electroencephalographic (EEG) signals are thus acquired in a non-invasive way, for its posterior processing and interpretation by the algorithms associated with the BMI. EEG activity is usually categorized by different rhythms which are associated with certain brain activities and specific frequency bands: activity below 4 Hz is related to Delta band [sleep waves (Amzica and Steriade, 1998) and motion related cortical potentials (Shibasaki and Hallett, 2006)]; the Theta band is usually assessed between 4 and 7 Hz and has been related, in the literature, to response repression (Kirmizi-Alsan et al., 2006); the Alpha band, which oscillates between 8 and 15 Hz, varies when the eyes are closed and it has been used with Beta band (16-31 Hz) for decoding the motor activity (Pfurtscheller et al., 2006); finally, the Gamma band (32-100 Hz) has been used in movement and attentive focus (Rao, 2013;Costa et al., 2016).
BMIs have emerged as a promising tool to rehabilitate or assist people when they are commanding robotic exoskeletons (He et al., 2018). In the case of the control of robotic exoskeletons by means of a BMI, the most common paradigms are based on: • Movement Related Cortical Potential (MRCP) are timefrequency changes of the signal associated with movements or decision making (Shibasaki and Hallett, 2006). However, their low potential and frequency makes them difficult to detect in a single trial due to the necessity to mitigate artifacts or other mental processes, averaging several events. • Event related Des/Synchronization (ERD/ERS) is based on the fact that during motor intention there is a previous decrement of the power in alpha and beta bands during the 2 s before the voluntary movement, followed by an increase of the power (Pfurtscheller and Neuper, 1994;Pfurtscheller et al., 1997Pfurtscheller et al., , 2006. For the clear identification of the ERD/ERS, it is also necessary tp average several trials. This paradigm has also been used during the mental task of motor imagery (MI), as it can produce a similar effect to the real movement (Jeon et al., 2011;Del Castillo et al., 2018). • Steady-State Visual Evoked Potentials (SSVEPs) are related to the response offered by the brain due to an external stimulus. Its nature is not directly related to the movement or MI, but this response can be used to control an exoskeleton with a BMI (Kwak et al., 2015;Zhang et al., 2015;Gui et al., 2017). Although the potential is not related to a movement action, it can be used to train the user to make a connection between the visual stimulus and the desired robotic action.
In these paradigms, the extraction of features is usually accomplished by the filtering of the frequency bands, associated with the rhythms related to the brain activity being analyzed. Although the features are generally extracted using traditional signal processing filters, there are other alternatives in the literature based on time vs. frequency analysis techniques, such as the wavelet transform (Xu and Song, 2008;Yang et al., 2010;Kant et al., 2019) and the Stockwell transform (Ortiz et al., 2017). However, they have not been used in conjunction with an exoskeleton (He et al., 2018). In the present study, a different approach will be studied using a decomposition algorithm, the Empirical Mode Decomposition (EMD). In this study, the EMD algorithm will be used to detect the variations of the EEG signals during MI tasks in comparison to the relaxed state, when a user is commanding a lower-limb exoskeleton. Due to the oscillatory nature of EEG and the lack of stationary behavior of EEG rhythms, EMD is an interesting algorithm for to use for the analysis of these signals. EMD has been used to detect epilepsy (Martis et al., 2012;Li et al., 2013) or even as an EOG removal technique (Looney et al., 2008), however, even as its utility for EEG analysis has been demonstrated (Rutkowski et al., 2010), its application for developing a BMI to command an exoskeleton has not been explored yet. The results will be compared to those obtained by a Butterworth filter. In addition, the activation of the exoskeleton will be analyzed in the time domain, thanks to the properties of the Hilbert transform.

Subjects
Three able-bodied adults voluntarily participated in the study. All the information was given to the subjects before the experiment FIGURE 1 | Image of a subject executing one of the experimental trials. and they agreed by signing an informed consent form. All procedures were approved by the Institutional Review Board of the University of Houston (USA).

Equipment
Two non-invasive bundles of 32 wet electrodes were used for the EEG acquisition with the help of an actiCap (Brain Products GmbH, Germany) unit. The distribution followed the 10 − 10 international system, reallocating four scalp electrodes around the eyes to assess the ocular artifacts in bipolar configuration (TP10-TP9 for Vertical-EOG and PO10-PO9 for Horizontal-EOG). Reference and ground were placed in the lower lobes of each ear. Data were wirelessly transmitted by a Move transmitter and amplified by two BrainAmp amplifiers (Brain Products GmbH, Germany). The lower-limb exoskeleton used was the REX (Rex Bionics, New Zealand). Activation commands were sent to the exoskeleton wirelessly, while the status of the exoskeleton was received by the computer through wire serial port communication. The control of the exoskeleton was done by custom firmware developed in Matlab. Figure 1 shows a subject commanding the REX exoskeleton during one of the trials. Although REX can be commanded without external help, two individuals were placed by the sides of the exoskeleton to avoid any potential risk of losing balance.

Protocol
The objective of the paper was to compare the assessment of the EEG variations observed in different rhythms associated with motor imagery tasks when a subject is using a lower-limb exoskeleton. With this aim in mind, the protocol was designed as detailed below.
The whole protocol of a experimental session lasted around 2 h and involved the following parts: (1) adaptation of the exoskeleton to the subject, (2) placing and gelling the 64 electrodes, achieving an impedance value below 20k for each electrode, (3) execution of several runs to make sure the subject is comfortable with the REX movement, and finally (4) the fulfillment of the 10 experimental trials. Figure 2 shows the details of an experimental trial, the different events, and the status of the exoskeleton.
Each trial started with 15 s of free mental task to help the subject concentrate in the experiment and to allow enough time to pass for the ocular artifact removal algorithm convergence (Kilicarslan et al., 2016). This period of time was not considered for the analysis, as it was not associated with any mental task. After that, an acoustic signal warned the subject to relax and to try imagine a rest state of no movement. This provided the reference state (bottom blue points) to the MI event (red points) that occurs after the rest event. After 10 s of rest, a new acoustic cue indicated to the subject to focus on the mental task of the movement of their legs (MI start). At the same time, a command was issued to move the exoskeleton. However, as REX has a variable time for activation, the start of the movement lagged for up to 2 s (start of first step point). After the first step was completed (around 5 s after the cue signal) the exoskeleton kept walking at a steady pace (cycling walk) recording the EEG signal for MI assessment (red top). After ∼20 s of MI, another cue signal indicated to the subject to start making a reverse count (black top line). This part serves as a control event to be sure that the MI features were not caused by motion artifacts. After another 20 s, a final acoustic cue indicated the final rest event. This final event was not used in this study as it had no associated MI event.

Empirical Mode Decomposition
EMD was developed by Norden E. Huang as a decomposition method to improve the extraction of the instantaneous frequency and amplitude of non-stationary signals (Huang et al., 1998). The algorithm decomposes the original signals in several modes, called intrinsic mode functions (IMF), to provide a better timefrequency representation using the Hilbert transform. Each IMF contains the local oscillation information of the signal, extracted from higher to lower frequencies.
Given a discrete signal s(t), its EMD can be obtained following the next steps: 1. Assign the signal to a temporal signal x(t) = s(t). 2. Find the local extrema of x(t). 3. Find the maximum envelope e + (t) of x(t) by passing a natural cubic spline through the local maxima. Similarly, find the minimum envelope, e − (t), with the local minima. 4. Compute an approximation to the local average: m(t) = (e + (t) + e − (t))/2. 5. Find the proto-mode function z i (t) = x(t) − m(t). 6. Check whether z i (t) is an IMF. An IMF is a wave that satisfies two conditions: the number of extrema and the number of zero crossings may differ by no more than one; and its local average is zero. The threshold used to set this last condition is critical, as a high value can lead to overtraining and a low value can leave components unextracted. Moreover, to avoid the extraction of accidental IMFs, the conditions must be accomplished in at least two or three consecutive iterations. 7. If z i (t) is not an IMF, repeat the loop on z i (t). If z i (t) is an IMF then set IMF n (t) = z i (t) and begin the process again considering x(t) = x(t) − IMF n (t), being n = 1 : M. 8. The process stops when the desired limit of IMFs has been achieved (n = M), or a maximum number of sifting operations have been done i = p, without the IMF conditions fulfilled. In this case, the residual is considered as res(t) = s(t) − n=1 n=M IMF n Nevertheless, the algorithm is not exempt of drawbacks, such as the difficulty to separate near tones or border effects (Rilling and Flandrin, 2008). Border effects and the time of computation can be a drawback for the EMD employment in real time EEG signals. As processing windows are usually short during online processing, a great part of the data is susceptible to being spoiled by the border effects of a bad envelope calculation during the sifting process. As this research performs an offline analysis of the data, it was exempt of this issue. Border effects were outside the processing windows, as the whole trial was used as input data to the EMD algorithm and not partially processed in epochs. However, its application in real time would require the use of overlapped processing epochs, neglecting the border content. Nevertheless, this effect was studied for the EEG signals of this research. For 1 s epochs overlapped by 0.5 s, the time of processing was around 0.3 s per epoch on average and 40% of the window was spoiled. As the time of computation depends on the hardware, and the border effect depends on the sampling frequency, it is difficult to recommend an epoch length. However, an epoch should be as long and with as short a shifting as the hardware allows for computational times.

Pre-processing
The data acquisition was done at a 1 kHz sampling frequency. However, in order to obtain a consistent set of IMFs, the data was resampled to 200 Hz. The number of IMFs and the frequency tone associated with them depends on the sampling frequency, so a 1 kHz sampling frequency would not provide relevant tones over the 100 Hz. Resampling to 200 Hz provide tones below 100 Hz. In addition, the ocular artifact removal algorithm parameters are optimized to work at 100-200 Hz sampling frequencies.
The H ∞ algorithm employs the information of the pairs of electrodes TP10-TP9 and PO10-PO9 (Kilicarslan et al., 2016). The optimized parameters used were the following: γ = 1.15, q = 10 −10 and p 0 = 0.5. Finally, a 60Hz notch filter was applied to mitigate the network power supply component.

IMF Extraction
Once the previous filters were applied, the EMD, as explained in section 2.2, was used to obtain the different modes of the nine electrodes located at the motor cortex: FC1, FCz, FC2, C1, Cz, C2, CP1, CPz, and CP2. Figure 3 shows an example of decomposition of one of the trials for the electrode Cz. Border effects can be seen at the tails of the signal, especially clear in the case of the end of IMF5. In this study, only the first five IMFs were considered (M = 5), as they were enough to obtain the EEG rhythms to be analyzed as instantaneous frequency indicators.
Instantaneous frequency (f (t)) and amplitude (A(t)) were computed using the Hilbert transform (Huang et al., 1998) of each IMF. Both can be obtained based on the derivative of the instantaneous phase and the module of the analytical function. Given a signal s(t) = IMF n (t), its analytical complex function z(t) can be obtained by applying the Hilbert transform (H()) as: obtaining the instantaneous frequency as a function of the derivative of the instantaneous phase: (2) Figure 4 shows the instantaneous frequency for the five IMFs considered. As it can be seen that IMF 1 oscillates in the Gamma band (40-60 Hz), IMF2 in the Beta band (20-30 Hz), IMF3 in the Alpha band (5-15 Hz), IMF4 in the Theta band (2-7 Hz), and IMF5 in the Delta band (<4 Hz). This associates each of the IMFs with a significant EEG rhythm. It is necessary to clarify that the order of the IMF and its relationship with each rhythm can differ depending on the sampling frequency and the possible noise content of the signal. This is one of the reasons why the signals were resampled to 200 Hz in order to have a decomposition similar to the EEG rhythms.

Event Assessment
Two different analyses were carried out. The first analyzed the behavior of the different EEG rhythms during the MI event.
In this analysis, the changes of power in the different IMFs were calculated taking the previous rest state as a reference. In order to avoid taking into account evoked potentials, the events considered only the stationary parts of the events. This means that only the power of the IMFs, when the exoskeleton was stopped or during stable walking, was considered (see Figure 2 bottom blue for rest and top red for MI). For a determined start event, the variation of power was assessed for each IMF (i) and electrode (ch) based on the power of the signal during the MI periods: and the rest periods: as the percentage of variation of the power taking the rest period as reference: As t = t 2 − t 1 can be different for rest and MI periods, power was divided by it.
In the second analysis, the rhythms were evaluated by its instantaneous variation in time. The instantaneous power of each IMF was computed based on its analytical function as: 3. RESULTS

Motor Imagery Analysis (Power Variation)
A statistic analysis was carried out to detect the significant differences for the P(%) between the subjects and the electrodes using SPSS (Andy, 2013). Data included the P(%) of the five extracted IMFs for each of the three subjects and the nine 16.9 ± 2.4 9.5 ± 1.3 −2.5 ± 1.6 −9.1 ± 4.7 −5.7 ± 5.7
A Shaphiro-Wilks test of normality (S-W) was applied to check for the normality assumption (Ghasemi and Zahediasl, 2012). Averaging the electrodes, S1 and S2 followed a normal distribution for all the IMFs (S1: D (10)  Using the individual electrode P(%) output for each IMF as a dependent variable, splitting the data for each subject and electrode, the S-W test results indicated that S1 followed a normal distribution for all the IMFs and electrodes (D(10) = Analyzing P(%) of the IMFs as dependent variables for each subject without splitting the data by electrode, the S-W test indicated that S1 followed a normal distribution for all the IMFs with only a small deviation from normality in the Q-Q plot for IMF1 [D(90) = 0.919, p < 0.05] and a small asymmetry in the histogram toward the high tail. S-W test for S2 indicated a deviation from the normal distribution for IMF1-4 (D(90) = [0.947 − 0.972], p < 0.05) with only a small skewness (S-shape) in the Q-Q plots. S3 distribution was not normal with a clear double peak shape for all the IMFs in the histograms (D(90) = [0.560 − 0.734], p < 0.05) and a clear kurtosis deviation in the Q-Q plots, which indicated different performance in some of the trials. As normality could not be assured for all the subjects and electrodes, data were FIGURE 6 | Instantaneous power of IMF3 after the acoustic cue for MI (black line) for trial 1 of subject S1. The value is compared with the instantaneous power of the signal filtered by a Butterworth filter (5-15 Hz). The peak response is higher using the IMF. REX start of movement is represented by a red spot.
analyzed using non-parametric tests to compare the subject and electrode dependency.
For the averaged electrodes, using the Kruskal-Wallis (K-W) test, the IMF output was significantly affected by the subject (H(2) = [9. 757, 11.027, 10.934, 17.559, 8.565], p < 0.05, appearing in the data for each IMF). The Mann-Whitney test was used to analyze the differences between the subjects. A Bonferroni correction was applied and so all effects are reported at a 0.0167 level of significance. S1 results were significantly different in comparison to S2 for 17,12,0,15] Analyzing each subject individually, the K-W test indicated that there were no differences between the electrodes for the five IMFs of all the subjects: S1 (H(8) = [3.721, 3.279, 10.252, 6.226, 7.732], p > 0.05), S2 (H(8) = [3.174, 2.183, 2.153, 6.357, 3.910], p > 0.05), and S3 (H(8) = [1. 496, 0.501, 6.790, 11.344, 4.984], p > 0.05). This means that the electrodes behaved in a similar way for a given IMF and subject. This matches the boxplot representation. The left part of Figure 5 shows the boxplots of the power variation for the average of the nine electrodes during the 10 trials accomplished by a subject, while the right part shows the same data for the Cz electrode. The associated descriptive statistics to the boxplots can be seen in Table 1 for all the electrodes and its average, and particularized for the Cz electrode in Table 2.
As can be seen, if the electrodes were averaged through the 10 trials, there was a similar trend for the three subjects, with a clear increment of the power, especially for high frequency IMFs associated with gamma (IMF1) and beta bands (IMF2). However, IMFs3-5 (alpha, theta, and delta bands) showed a more erratic behavior. Focusing our attention on one of the most representative electrodes, Cz (see right part of Figure 5), its boxplot follows the same trend than the average. However, the deviation, as can be seen in the boxplot and Table 2, is higher, which indicates the difficulty to see this pattern in a single trial without averaging the electrodes.
The comparison of the averaged P(%) calculated using EMD and the components extracted from the signal by a second order Butterworth filter, show that the results were similar ( Table 3). The table shows, using bold text, the methods with a higher increment for P(%) in the first two IMFs.

Gait Intention Analysis (Instantaneous Power)
In order to study the transient when a subject focuses on the motion intention, an additional study was carried out. The ERD/ERS phenomenon indicates that for certain subjects, the motor intention involves a desynchronization of the mu band (8-12 Hz) followed by a posterior synchronization (Pfurtscheller et al., 2006). Although the phenomenon can be observed in a single trial, it usually requires averaging several trials. Moreover, the effect is also more difficult to notice in foot imagery than in hands. In this study, EMD is proposed as a technique to evaluate the desynchronization and posterior synchronization of the signal during the motion intention. As literature focuses its attention on the mu band, the IMF under study was IMF3, as its instantaneous frequency oscillated around 10 Hz (see Figure 4 IMF3).
With this purpose in mind, the instantaneous power of IMF3 was assessed as indicated by equation 6. Figure 6 shows the instantaneous power during the 5 s after the acoustic cue was issued (t = 0) for one of the trials (black line). A red spot indicates the real start of the REX movement, which indicates that the previous time to the red spot is a MI period with no movement, exempt for this reason of any possible motion artifacts. The instantaneous power was compared to the one obtained by only using a second order Butterworth bandpass filter to the signal between 8 and 15 Hz (blue line). Upon inspection of the figure, it can be seen that there are several peaks during the motor imagery period previous to the activation of the exoskeleton. For several electrodes there is a first peak around 1 s before the movement and another one just before the movement. As can be seen, the peak that appears just before the start of the movement is especially clear, using the IMF power (black line), compared to that of the filtered signal (blue line). This indicates that there is an increment of power around 1 s after the acoustic cue, followed by a decrement and a posterior increment just before the real movement starts, which could be related to the ERD/ERS phenomenon. However, the peak magnitudes are variable, which makes it difficult to quantify the phenomenon in a single trial. Table 4 shows a comparison of both methods for the time at which the maximum peak of the instantaneous power is computed in the five s after the acoustic cue for the MI event is issued. The mean time shows the averaged time for the 10 trials accomplished by each subject per electrode. A negative time indicates that the peak achieved was average for a channel before the REX activation. According to K-W test, the subjects showed no significant differences when using EMD [H(2) = 5.42, p > 0.05], while they did for the frequency filter band [H(2) = 27.13, p < 0.001]. As Table 4 shows, a negative time was achieved for the average of the 10 trials in all the subjects for the electrodes C1, CPz, Cz, and FC2 when using the EMD method. The filtered signal showed a higher dispersion for the times, detecting the intention later for the majority of the electrodes. Furthermore, Figure 6 shows that there are also several peaks for certain electrodes after the exoskeleton started moving. To check the evolution and magnitude of these peaks and its possible relationship with the statistical results of section 3.1, the instantaneous power is represented for the full trial in Figure 7. This was limited to IMF1 and IMF2, as these IMFs were the ones that showed a clearer increment of power during MI in Table 2. To see the behavior in a clearer way, the instantaneous power was averaged with a moving mean of 1 s. Figure 7 shows this averaged value for the MI (red), rest (blue), and the control period of the reverse count (green) events. Upon inspection of Figure 7, it is evident that the MI periods show a clear increment in rest periods. In addition, they also show an increment in reverse count operations, which indicates that the increment is not related to any possible motion or exoskeleton artifact, since during the reverse count periods the exoskeleton was moving in the same way than during the MI events. In fact, reverse count periods show a similar power to rest periods.

DISCUSSION
By analyzing the results, it has been demonstrated that EMD can be a useful tool in characterizing the EEG rhythms during MI events associated with a BMI operating an exoskeleton.
Regarding motor imagery, an increment of power was noticed, not only in average for the MI vs. rest periods, but also in comparison to an alternative mental task during movement (reverse count) in the case of the high frequency IMFs. The power for the IMFs associated with these IMFs showed a different pattern only for the MI events, while showing a similar one for the rest periods. This is consistent with the literature on gamma power increment during motion (Seeber et al., 2014). This behavior could be explained by the attention focus to gait during these periods (Costa-García et al., 2019), which could be useful for the future development of a BMI based on high frequency rhythms assessing the attention to gait. On the other hand, lower frequency IMFs showed a more inconsistent behavior during MI tasks. This can be explained by the nature of EMD, as low frequency IMFs are more susceptible to over-training problems. In addition, in the case of a real-time application, shorter windows of computation could make it more difficult to obtain a consistent result using these IMFs due to possible border effects. This advise against the use of high order IMFs (low frequency ones), at least for BMI control, even as the residual captures the low frequency oscillations of the EEG channel in a precise way (see Figure 3). Furthermore, the statistical analysis revealed that the results were dependant on the subject. However, given a determined subject and IMF, they have a similar behavior for the nine electrodes considered. Nevertheless, the high deviation of the P(%) per individual electrode makes its use in a single trial difficult, without averaging the nine electrodes through the trials.
Regarding motor intention, the analysis revealed important facts. As Figure 6 shows, motor intention analyzed by the IMF associated with alpha band followed a pattern similar to the one explained by the ERD/ERS (Pfurtscheller et al., 2006;Jeon et al., 2011) in a single trial. The effect was more noticeable than when using a frequency filter. However, it requires a processing window of at least 4 s, which could compromise the BMI response in real time, as processing would delay the command and increase the computational times. In addition, the magnitude of the peaks was variable, which makes computing them in a significant stable base, or to estimate a statistical trend to be used in a BMI algorithm, difficult, even as statistically the maximum peak is usually achieved during the first 5 s before the actual start of movement.
The results obtained for the P(%) were similar to the ones obtained for the filtered signal, but clearer for the motor intention despite its higher computational load. Nevertheless, EMD could be an interesting alternative for the assessment of ERD/ERS, but would require further investigation and an experimental protocol specifically designed for its study. Future research will explore the viability of using this paradigm for the detection  The maximum peak is computed for each electrode within the 5 s after the MI acoustic cue is issued. Zero reference is marked by the instant the exoskeleton starts moving. A negative time indicates that the detection is accomplished before the exoskeleton starts moving. The table shows the averaged time for the ten trials of each subject and its 95% confidence interval. Electrodes with negative times are in bold and methods with negative times for all the subjects appear between brackets. Frontiers in Neurorobotics | www.frontiersin.org of motor intention through the analysis of moving epochs and the usability of high frequency features, extracted with EMD or other time-frequency algorithms, for a BMI that commands an exoskeleton in real time.

DATA AVAILABILITY STATEMENT
The datasets presented in this article are not readily available because of University of Houston IRB restrictions.
Requests to access the datasets should be directed to JC-V (jlcontreras-vidal@uh.edu).

ETHICS STATEMENT
The studies involving human participants were reviewed and approved by Institutional Review Board of the University of Houston (USA). The patients/participants provided their written informed consent to participate in this study. Written informed consent was obtained from the individuals for the publication of any potentially identifiable images or data included in this article.

AUTHOR CONTRIBUTIONS
MO and JA: conceptualization. MO: methodology, validation, formal analysis, investigation, writing-original draft preparation, and visualization. MO and EI: software and data curation. JC-V: resources. MO, EI, JC-V, and JA: writing-review and editing. JC-V and JA: supervision and project administration. MO, JC-V, and JA: funding acquisition. All authors contributed to the article and approved the submitted version.

FUNDING
This research was funded by the Spanish Ministry of Education and Professional Formation through grant CAS18/00048 José Castillejo; by the Spanish Ministry of Science and Innovation, the Spanish State Agency of Research, and the European Union through the European Regional Development Fund in the framework of the project Walk-Controlling lower-limb exoskeletons by means of brain-machine interfaces to assist people with walking disabilities (RTI2018-096677-B-I00); and by the Consellería de Innovación, Universidades, Ciencia y Sociedad Digital (Generalitat Valenciana) and the European Social Fund in the framework of the project Desarrollo de nuevas interfaces cerebro-máquina para la rehabilitación de miembro inferior (GV/2019/009).