Noise-assisted multivariate empirical mode decomposition for multichannel EMG signals

Background Ensemble Empirical Mode Decomposition (EEMD) has been popularised for single-channel Electromyography (EMG) signal processing as it can effectively extract the temporal information of the EMG time series. However, few papers examine the temporal and spatial characteristics across multiple muscle groups in relation to multichannel EMG signals. Experiment The experimental data was obtained from the Center for Machine Learning and Intelligent Systems, University of California Irvine (UCI). The data was donated by the Nueva Granada Military University and the Technopark node Manizales in Colombia. The databases of 11 male subjects from the healthy group were taken into the study. The subjects undergo three exercise programs, leg extension from a sitting position (sitting), flexion of the leg up (standing), and gait (walking), while four electrodes were placed on biceps femoris (BF), vastus medialis (VM), rectus femoris (RF), and semitendinosus (ST). Methods Based on the experimental data, a comparative study is provided by assessing the Empirical Mode Decomposition (EMD)-based approaches, EEMD, Multivariate EMD (MEMD), and Noise-Assisted MEMD (NA-MEMD). The outcomes from these approaches are then quantitatively estimated on the basis of three criterions, the number of Intrinsic Mode Functions (IMFs), mode-alignment and mode-mixing. Results Both MEMD and NA-MEMD methods (except EEMD) can guarantee equal numbers of IMFs. For mode-alignment and mode-mixing, NA-MEMD is optimal compared with MEMD and EEMD, and MEMD is merely better than EEMD. Conclusions This study proposes the NA-MEMD approach for multichannel EMG signal processing. This finding implies that NA-MEMD is effective for simultaneously analysing IMFs based frequency bands. It has a vital clinical implication in exploring the neuromuscular patterns that enable the multiple muscle groups to coordinate while performing the functional activities of daily living.


Background
Electromyography (EMG) is the collective electric manifestation during muscle contraction, and indicates the electrophysiological responses of motor units from a muscle group, which is controlled by the nervous system. The surface EMG signal, originating in motor units and then recorded by measurement tools, was often contaminated by various types of noises or artifacts, e.g., power line interference, baseline wandering, electrocardiographic (ECG) artifacts, capacitive effects of the detection site, and the firing rate of motor units [1][2][3][4][5][6]. Therefore, the identity of an actual EMG still remains difficult [7][8][9].
Recently, several methods have been developed to analyse and de-noise the EMG signal [10][11][12]. The conventional techniques based on Fourier analysis (e.g., IIR filters) are widely used for EMG-based filtering. However, Fourier analysis is purely based on predefined basis functions, which not only reduces the noise but also attenuate the EMG signal. As an alternative to the usual Fourier transform method, wavelet analysis is also popularised due to its advantages in terms of the time-frequency representation [13][14][15][16]. The wavelet-based approaches, however, are also suboptimal because the preselected wavelet function is often not suitable for matching the natural property of an EMG signal. Previous studies have also introduced the Empirical Mode Decomposition (EMD) approach to handle EMG signals [17]. Instead of those reported in literatures [18], the EMD is a fully data-driven adaptive time-frequency analysis method, and offers no prior assumption through the overall data processing procedure [19][20][21][22].
The EMD algorithm was put forward by Huang et al. and provided the most successful results for the decomposition and time-frequency analysis of non-stationary signals. It is given as a sifting process that decomposes a signal into a finite set of intrinsic mode functions (IMFs), amplitude-and/or frequency-modulated components representing its inherent oscillatory modes. Adriano et al. first employed this technique to filter EMG signals in background activity attenuation [17]. However, the first version of EMD was only used for a single-channel EMG, and did not focus on the accuracy of the decomposed subfrequency bands. In order to alleviate this problem, the Ensemble EMD (EEMD), an adaptive dyadic filter bank, was introduced. This method can effectively eliminate the mode-mixing and physically produce more unique frequency levels. The literature shows that several studies have investigated the de-noising performance for EMG signals using the EEMD algorithm [23]. However, such single-channel based EMD algorithms cannot be directly applied into the multiple-channel EMG signal processing [24]. Moreover, the EMD or EEMD algorithms cannot guarantee the equality of the number of decomposed IMFs across multichannels, and may lead to subsequent EMG-based analyses being physically meaningless. Accordingly, the multivariate extension of EMD (MEMD) and its noise-assisted analysis method, Noise-Assisted Multivariate EMD (NA-MEMD) have been developed recently to produce the same number of IMFs across all channels thereby facilitating direct multichannel analyses with the consideration of cross-channel interdependence (mode-alignment) and single-channel independence (mode-mixing) [25][26][27][28][29][30].
EEMD has been extensively applied as an accurate and computationally efficient quantitative analysis for electromyography (EMG) signals. The EEMD algorithm can effectively extract the temporal information of EMG time series. However, few papers examine the temporal and spatial characteristics across multiple muscle groups in relation to multichannel EMG signals. In this study, NA-MEMD is proposed to handle the multichannel EMG signal processing. The performance of the proposed method has been validated by comparing it with EEMD and MEMD. The experimental data was obtained from the Center for Machine Learning and Intelligent Systems, University of California Irvine (UCI). The data was donated by the Nueva Granada Military University and the Technopark node Manizales in Colombia. Three criterions are proposed to assess the decomposition performance, (1) the number of intrinsic mode functions; (2) mode-alignment (common frequency scales in the same indexed IMFs across different channels) for the cross-channel interdependence; (3) mode-mixing (a single IMF containing multiple scales and/or a single scale residing in multiple IMFs) for the singlechannel independence. Results indicate that both MEMD and NA-MEMD methods (except EEMD) can guarantee equal numbers of IMFs. Specifically, for mode-alignment and mode-mixing, NA-MEMD is optimal compared to MEMD and EEMD, and MEMD is merely better than EEMD.

Experiments
The experimental data was obtained from the Center for Machine Learning and Intelligent Systems, University of California Irvine (UCI). The data was donated by the Nueva Granada Military University and the Technopark node Manizales in Colombia. UCI consented to cite these datasets in publications [31]. This work is also approved by the Institution Research Ethics Board of University of Electronic Science and Technology of China (UESTC). The databases of 11 male subjects from the healthy group are taken into the study. The subjects undergo three exercise programs associated with the knee joint, leg extension from a sitting position (sitting), flexion of the leg up (standing), and gait (walking), while four electrodes are placed on biceps femoris (BF), vastus medialis (VM), rectus femoris (RF), and semitendinosus (ST). The goniometer is also used to record the angle of the knee joint during the exercise programs. Each subject is asked to perform these exercise programs once, and each exercise program contains approximately five motion repetitions. The period time of motion is about 4 s, 2 s for motion and 2 s for rest.

EEMD
The Empirical Mode Decomposition (EMD) is a fully data-driven and adaptive timefrequency analysis method. It describes a signal as a linear combination of a finite set of intrinsic mode functions (IMFs) and a residual signal. The mathematic representation of EMD can be depicted as where x(t) is an original signal, c m (t) and r(t) represent the m th IMF and the residual assumed as the (M + 1) th IMF, respectively. These resultant IMFs, c m (t) M m=1 , are sequentially extracted from the original signal by an iteration algorithm called the  sifting processing [19]. For the EMD-based sifting process, the local maxima and minima of x(t) are first identified, and the upper (lower) envelope is constructed by fitting the local maxima (minima) into a cubic-spline curve. The averaged curve of upper and lower envelopes is then intended to update x(t) by subtracting it from x(t). This sifting process will be iteratively executed until each IMF can be determined, while two stoppage criterions should be satisfied, i.e., (a) the number of zero crossings and the number of extrema (inclusive of the total number of the local maxima and minima) should not differ by more than one; (b) the average value of the upper envelope and the lower envelope through the overall data should be zero. After repeating the above sifting process, all IMFs c m (t) M m=1 are obtained when the residual r(t) becomes a monotonic function.
For EEMD, the extra white Gaussian noises (WGNs) w(t) are added with the original signal x(t) to obtain an ensemble of noise-assisted signal s(t), i.e., s(t) = x(t) + w(t), and the ensemble signal is decomposed by using the EMD algorithm. This single noise-added procedure is then repeatedly executed, and for each iteration the different realization of white noise w n (t) is given where n = 1, 2, . . . N representing the number of iterations that is set to 50 in this study. The final IMFs can be calculated by averaging the same indexed IMFs of the decomposition. The EEMD algorithm is provided as follows [29]: Identify all local extrema of x n (t); (4) Find lower and upper envelopes, e n l (t) and e n u (t), which interpolate all local minima and maxima, respectively; (5) Calculate the local mean, m n (t) = 1 2 (e n l (t) + e n u (t)); (6) Subtract the local mean from x n (t), c n m (t) =x n (t) −m n (t) (n is the index number of IMF); (7) Let x n (t) = c n m (t) and go to step 3); repeat until c n m (t) becomes IMFs; (8) Average the corresponding IMFs from the whole ensemble to obtain the averaged IMFs; for instance, the m th IMF can be obtained by using c m (t) = 1 N ( N n=1 c n m (t)).

MEMD
Rehman and Mandic developed multivariate empirical mode decomposition (MEMD), which is a natural extension of the original EMD/EEMD. In MEMD, the multiplechannel EMGs should be first projected into n-dimensional spaces based on low discrepancy Hammersley sequence. The projections along different directions in multidimensional spaces represent the amplitudes of EMGs across four channels. The extrema are interpolated via cubic-spline interpolation in order to obtain the subenvelopes e θ v (t) V v=1 . Those sub-envelopes are then averaged to obtain a local mean of a multiple-channel EMG signal, M(t). Then, the first IMF can be extracted by subtracting the local mean from the input channels. The outline of MEMD algorithm is presented as follows [25]: (1) Choose a suitable point set for sampling a (p − 1) sphere; (2) Calculate a projection, w θ v (t), of N-channel input signals x N (t) (N = 4) along the direction vector d θ v , for all v (the whole set of direction vectors), giving w θ v (t) V v=1 as the set of projections; If c N (t) fulfills the stopping criterion for a multivariate IMFs, apply the above procedure to x N (t) − c N (t); otherwise, apply it to c N (t).
Different with EEMD, the sifting process is followed by , and the threshold value γ was set to 0.2 based on the EMG signals during lower-limb exercises. The sifting process will be continued if Eq. (2) is satisfied.

NA-MEMD
The Noise-Assisted multivariate empirical mode decomposition (NA-MEMD) was also introduced by Rehman and Mandic. This signal processing work exploits the dyadic filter properties of EMD and MEMD. Additionally, it also applies the noise assisted analysis method into MEMD, a dyadic filter bank on each channel while adding certain multidimensional WGNs together with the original signals which are decomposed by using MEMD. More specifically, K-channel (K ≥ 1) uncorrelated WGNs time series of the same length with that of the M-channel EMGs (M = 4) are randomly separately created. Then, a new input multichannel signal is constructed by adding the original EMGs with the noise channel, the resulting (M + K ) − channel multivariate signal. Considering the decomposition of the constructed signal, the remaining procedures are strictly followed by those of MEMD [32]. Figure 3 outlines the processing procedure of NA-MEMD. The effects of the number of noise channels and noise power in NA-MEMD are discussed in [33]. In this study, the average STD based on all EMG channels is selected as the residual noise power, and the number of noise channel is set to four. The schematic diagram for methods of EEMD, MEMD, and NA-MEMD is presented in Fig. 1.

Data preprocessing
The raw EMG data measured from each subject were first segmented. The period of the exercise motion was reserved and labeled. The approaches of EEMD, MEMD, and NA-MEMD are then used to decompose these segmented EMG data, by which the decomposed IMFs are obtained via three methods, and then normalized by using standard deviation. Based on each single normalized IMF data, the alignment of IMF based frequency bands in cases of EEMD, MEMD, and NA-MEMD is estimated by the spectra analysis. The last step is to evaluate three criterions, the number of IMFs (indicating cross-channel interdependence), mode-alignment (estimating the alignment of the frequency bands of the same-index IMFs across channels), and mode-mixing (estimating the similarity of the frequency bands of IMFs within a single channel). Figure 2 depicts the schematics for data flow and evaluation criterions, the number of IMFs, mode-alignment, and mode-mixing.

Mode-alignment
The mode-alignment effect indicates the common frequency scales in the same indexed IMFs across different channels. This effect would numerically analyse the correlations of frequency scales for each component of the EMG channels, and take advantages of comparatively analysing the frequency similarity of the same-indexed IMFs across channels.
In order to obtain this performance, the power spectral density (PSD) of the normalized IMF is first calculated. The PSD correlations between two IMFs are then obtained by c i,j , where i stands for the ith indexed IMF, and j is for the number of the channel. The correlation matrix for all IMFs across channels could be expressed as The elements in the ith row of the correlation matrix are averaged to represent the mode-alignment value in the ith indexed IMF.

Mode-mixing
The mode-mixing effect describes the overlap of frequency information among the decomposed IMFs within one EMG channel, which would reflect whether or not a single IMF contains multiple scales and/or a single scale resides in multiple IMFs [34]. In this study, we used the following equation to quantitatively describe the mode-mixing effects, MM i,j , where f 2i , f 8i , and D i are the PSD of the ith indexed IMF. f 2 , f 8 are the frequencies at which 20 and 80% of the energy of an IMF are reached, respectively. D i is the difference between f 2i and f 8i . Based on Eq. (4), the mode-mixing effect for a single EMG channel could be calculated as c(2, j) . . . . . . . . . . . .   where I is the total number of IMFs. Figure 3 shows an example of the decomposition result in the vastus medialis muscle group for three exercise programs (sitting, standing, and walking) for EEMD, MEMD and NA-MEMD. Since the most predominant energy for an EMG signal is approximately between 20 and 500 Hz [5], the decomposed components that have lower subfrequency bands than 20 Hz are synthesized together (from the 5th to 11th IMFs, the 6th to 16th IMFs, and the 7th to 16th IMFs for EEMD, MEMD, and NA-MEMD, respectively).

Spectra analysis
In order to analyse IMF based frequency components produced by EEMD, MEMD, and NA-MEMD, the decomposed IMFs, specifically representing one of exercise motions, are first normalized. These IMFs are then utilized for the analysis of spectra. Figure 4 indicates the spectra results of IMFs for three exercise motions decomposed by EEMD, MEMD, and NA-MEMD. In this study, we only focus on the shape of individual spectra in considerations of mode-alignment and mode-mixing. The alignment of frequency bands of the same-index IMFs in muscles BF, VM, RF, and ST is closer, the mode-alignment performance is more prominent. The spectra figures only can qualitatively analyze and demonstrate the differences of decomposition by EEMD, MEMD, and NA-MEMD, in which the stabilization of the shape of individual spectra from BF, VM, RF and ST can be observed. Based on these spectra information, the statistical analyses are used to quantitatively estimate the performance of mode-alignment and mode-mixing.

The number of IMFs
A statistical survey is also taken by investigating the EMG signals of muscle groups RF, BR, VM, and ST to sitting, standing, and walking exercises for all subjects. The averaged number of IMFs for each muscle is shown in Table 1. It has been clearly shown that MEMD and NA-MEMD could guarantee the equal number of IMFs across EMG different channels. In addition, the number of IMFs via MEMD and NA-MEMD have a larger amount compared to those of EEMD, indicating that more details of EMG frequency components can be obtains based on MEMD and NA-MEMD results.

Mode-alignment
In order to statistically analyse the mode-alignment performance for multiple-channel EMGs, the correlation matrixes based on the motion segmentations of four-channel EMG signals from all subjects in three exercise programs are calculated. The IMFs with the subfrequency energy less than 20Hz are removed as it contains much noise and has a low signal-to-noise ratio. The mode alignment effects of decomposed IMFs of four-channel EMG data obtained from the health group are identified in Table 2. Based on these results, twoway analysis of variance (ANOVA) is used to examine the influence of exercise programs (i.e., sitting, standing, and walking) and methods (i.e., EEMD, MEMD, and NA-MEMD)  on the performance index of mode-alignment ( Table 3). The assessment results show that the methods have a significant main effect (p < 0.01), and no interaction between exercise programs and methods (p > 0.05). The statistic analysis with no interaction effect confirms that the three types of exercise programs equally represent the characteristic of functional activities of daily living. In order to further evaluate the difference among methods, the mode-alignment values in three exercise programs for each subject are averaged, and then the One-way Repeated Measures ANOVA is used to compare the mode alignment of IMFs by EEMD, MEMD, and NA-MEMD. It is clear that there is a significant difference among three methods (F = 32.022, p = 0.000). By using the Least Significant Difference   (LSD), the mode-alignment effect of NA-MEMD is the best among three methods, and the effect of MEMD is merely better than that of EEMD (Fig. 5).

Mode-mixing
In this study, we also investigate the mode-mixing effect based on each muscle channel by using Eqs. (4) and (5). In order to avoid the effects of inter-subject variability, the mode-mixing effects from all subjects are investigated. For each single subject, the decomposed IMFs for each muscle channel with a central frequency of the spectrum less than 20 Hz are also removed. The mode-alignment MM i,i+1 (i = 1, 2, . . . , I)from the remaining IMFs for each muscle channel are calculated. The mode-alignment effect for each muscle channel M is then obtained by averaging the set of MM i,i+1 (i = 1, 2, . . . , I).
Following this procedure, the mode-mixing effects of four EMG channels for three exercise programs are provided in Table 4. In a similar way, the influence of exercise programs and methods on mode-mixing is first quantitatively analyzed through two-way ANOVA. Table 5 indicates the main effect of methods (p < 0.01) as well as no interaction effect between two factors (p = 0.706). The mode-mixing values in the three exercise programs of each subject are averaged, and the averaged values are applied to test the influence of methods on the performance of mode mixing by using One-way Repeated Measures ANOVA and LSD. It can be seen from Fig. 5 that there are significant differences between two methods (EEMD vs. MEMD, EEMD vs. NA-MEMD, and MEMD vs. NA-MEMD) for the performance of mode mixing of decomposed IMFs. The NA-MEMD achieves the best mode mixing performance compared to EEMD and MEMD, while MEMD far outperforms that of EEMD.

Discussion
The objective of this study is to evaluate a superior solution for the preprocessing of multichannel EMG signals as well as for the analysing of the IMF based frequency components related to multiple muscle groups. The muscle coordination often occurs in human motions, which is not only indicated by multichannel EMG signals, but also conducted by neuromuscular patterns [35]. Generally, the neuromuscular pattern is intrinsic for the specific exercise motions. Therefore, the single-channel-based analyses for the observation of the nervous system and its corresponding muscle contraction are not sufficient.
Additionally, similar with the ECG lead system [36], it is also desirable to develop the EMG lead system in which the behaviors of motor units can be represented as a set of statistically independent sources.
The human exercise is often supported by multiple relative muscles. For example, the muscle groups of BF, VM, RF and ST are the muscles related to the knee movement such as standing, sitting, and walking. Hence, the use of the so-called four-lead system (the leads placed on muscles BF, VM, RF and ST) would well indicate the overall neuromuscular patterns, which are further controlled by the human brain activity. Moreover, it is   a natural way to simultaneously decompose the multichannel EMG signals and analyse the subfrequency bands of multichannel EMG signals. Although previous literatures have reported the successful applications of EMD/ EEMD in the single-channel EMGs [17], these approaches cannot solve the critical problem about the fusion and analysis of multichannel EMG signals [30,34,37]. Therefore, EEMD, MEMD and NA-MEMD have been investigated in this study for the decomposition performance of four knee muscle groups associated with standing, sitting, and walking. Three criterions (the number of IMFs, mode-alignment and mode-mixing) are employed to quantitatively depict the decomposition efficiency.
It has been confirmed that both MEMD and NA-MEMD (exclusive of EEMD) could provide an equal number of IMFs across EMG different channels. If the number of IMFs is unequal, then the decomposed subfrequency signals cannot be directly applied for the subsequent study. This also leads to the similar oscillation modes appearing in multiple IMFs (Fig. 3a).
The mode-alignment effect focuses on the cross-channel dependence. In order to compare the same indexed IMFs among muscle channels, a similar subfrequency band of the same indexed IMFs should also be observed. The statistics show that there is a significant difference among three methods (F = 32.022, p = 0.000). Moreover, the effect of NA-MEMD is the best among three methods. In addition, the effect of MEMD is better than that of EEMD.
For the assessment of the mode-mixing effect, there are significant differences between two methods (EEMD vs. MEMD, EEMD vs. NA-MEMD, and MEMD vs. NA-MEMD). Specifically, NA-MEMD achieves the best mode-mixing performance compared to EEMD and MEMD, and the effect of MEMD outperforms that of EEMD.

Limitation
The experimental data in this study was obtained from the Center for Machine Learning and Intelligent Systems, UCI. The data was donated by the Nueva Granada Military University and the Technopark node Manizales in Colombia. The physical characteristics of the participants were not recorded in the datasets. There also was no information about the prior nutritional intake, physical activity and environment conditions before all participants engaged in the experimental sessions. In addition, as the exercise programs (i.e., sitting, standing, and walking) are only taken from one measurement, the intrasubject variability such as random errors may not be avoided. The experiment description contained in the datasets did not clearly specify the location of electrodes placed on muscles BF, VM, RF, and ST.

Conclusions
This study proposed the noise-assisted multivariate empirical mode decomposition (NA-MEMD) approach for the preprocessing of multiple channel EMG signals, by which the temporal and spatial characteristics across multiple muscle groups can be quantitatively depicted. The four muscle groups of BF, VM, RF, and ST associated with lower limb exercises (sitting, standing, and walking) of 11 healthy subjects were utilised for the assessment of the EMD-based approaches. A comparative study was provided by assessing the NA-MEMD with Ensemble Empirical Mode Decomposition (EEMD), and Multivariate EMD (MEMD). Three criterions were used to assess the comparative outcomes, i.e., the number of intrinsic mode functions (IMFs), mode-alignment and modemixing. The results indicated that the current EMD-based approach of using EEMD was suboptimal for multichannel EMG signals due to its poor performance in relation to the three criterions. When compared with MEMD and NA-MEMD, both approaches with data from lower limb EMG signals would guarantee an equal number of IMFs across channels. In addition, the statistical results showed that both the mode-alignment and mode-mixing effects of NA-MEMD were superior to those of MEMD. This finding implied that NA-MEMD is effective for simultaneously analysing IMFs based frequency bands. It has a vital clinical implication in terms of exploring the neuromuscular patterns that enable the coordination of multiple muscle groups for the purposes of performing daily activities.