Multichannel Surface EMG Decomposition Based on Measurement Correlation and LMMSE

Ning, Yong; Zhao, Yuming; Juraboev, Akbarjon; Tan, Ping; Ding, Jin; He, Jinbao

doi:https://doi.org/10.1155/2018/2347589

Journal of Healthcare Engineering

On this page

Abstract Introduction Materials and Methods Results and Discussion Discussion Conclusions Data Availability Conflicts of Interest Acknowledgments References Copyright Related Articles

Special Issue

Smart Wearables in Healthcare: Signal Processing, Device Development, and Clinical Applications

View this Special Issue

Research Article | Open Access

Volume 2018 | Article ID 2347589 | https://doi.org/10.1155/2018/2347589

Multichannel Surface EMG Decomposition Based on Measurement Correlation and LMMSE

Yong Ning,^1,2Yuming Zhao,³Akbarjon Juraboev,¹Ping Tan,¹Jin Ding,¹and Jinbao He⁴

Academic Editor: Chengyu Liu

Received19 Jan 2018

Revised09 Apr 2018

Accepted14 May 2018

Published28 Jun 2018

Abstract

A method based on measurement correlation (MC) and linear minimum mean square error (LMMSE) for multichannel surface electromyography (sEMG) signal decomposition was developed in this study. This MC-LMMSE method gradually and iteratively increases the correlation between an optimized vector and a reconstructed matrix that is correlated with the measurement matrix. The performance of the proposed MC-LMMSE method was evaluated with both simulated and experimental sEMG signals. Simulation results show that the MC-LMMSE method can successfully reconstruct up to 53 innervation pulse trains with a true positive rate greater than 95%. The performance of the MC-LMMSE method was also evaluated using experimental sEMG signals collected with a 64-channel electrode array from the first dorsal interosseous muscles of three subjects at different contraction levels. A maximum of 16 motor units were successfully extracted from these multichannel experimental sEMG signals. The performance of the MC-LMMSE method was further evaluated with multichannel experimental sEMG data by using the “two sources” method. The large population of common MUs extracted from the two independent subgroups of sEMG signals demonstrates the reliability of the MC-LMMSE method in multichannel sEMG decomposition.

1. Introduction

Electromyographic (EMG) signals are comprised of action potentials produced by the muscle fibers contained in different motor units (MUs) [1]. It is of great importance for physiological investigation and clinical diagnosis to decompose EMG signals into their constituent motor unit action potential (MUAP) trains. EMG signal decomposing will lead to a better understanding of the properties of MU control and reveal the MUAP changes due to muscle fiber denervation/reinnervation [2]. It will also aid in the examination of neuromuscular diseases (e.g., amyotrophic lateral sclerosis) and the process of evaluating the degree of dysfunction found in upper motoneuron diseases such as Parkinson’s disease [3], cerebral palsy [4], hemiparetic stroke [5], and other disorders [6, 7]. Furthermore, EMG decomposition can facilitate the study of the interpulse interval (IPI) variability [8], recruitment strategies [9], myoelectrical manifestations of fatigue [10], and short-term MU synchronization [11].

EMG signals can be detected by introducing a fine wire or needle sensor into the muscle tissue or by placing sensors on the surface of the skin. In the course of studying these EMG signals, it has been found that the surface detection of EMG provides several advantages over wire or needle detection. For example, surface electrodes can be used quickly and easily, without causing discomfort for the subject or requiring medical supervision [12], and measurements can be performed with a high degree of repeatability. More importantly, surface EMG (sEMG) is also able to obtain global information about muscle activities and consequently records a vast amount of information [12]. This makes it more convenient for studying neuromuscular control mechanisms than the invasive methods, which offer less information about global muscle activities and are more difficult to utilize.

Over the past few decades, great strides have been made in decomposing indwelling EMG signals [13–15]. However surface EMG decomposition remains a difficult task. There is routinely a high level of action potential overlapping and cancellation within sEMG signals. The volume conduction effect for propagating action potentials is also enhanced in surface recordings due to the relatively large distance between electrodes and sources [12]. In addition, there exists the spatial integrating effect caused by surface electrodes. Hence, the differences in surface action potential shapes from different MUs are not as distinguishable as with intramuscular recordings [12]. Together, all of these factors make sEMG decomposition an extremely difficult task, especially at high force levels.

Various approaches for sEMG decomposition have been proposed over the past years in both sEMG recording and processing [16–21]. In particular, the design of surface electrode arrays comprised of a number of tiny electrode probes with a small interelectrode distance promises to increase the motor unit discrimination capacity by reducing MUAP superimposition while providing spatial information across the muscle. The extraction of a single MU from sEMG has also become a feasible task at very low force levels with appropriate signal processing methods such as two-dimensional template matching [20]. Recent developments in sEMG decomposition have further allowed for the extraction of a number of simultaneously firing MUs at relatively high force levels. Nawab et al. [17] developed a remarkable sEMG decomposition technique using a specially designed 5-pin Laplace electrode array in conjunction with a knowledge-based artificial intelligence framework. Holobar and Zazula proposed the convolution kernel compensation (CKC) method [18] and the gradient CKC approach (GCKC) [19] to decompose multichannel sEMG signals recorded with high density electrode arrays. It has been demonstrated that the GCKC method holds the promise of high efficiency and a strong antinoise performance in sEMG decomposition [19], but it has a strict requirement for the length of the EMG signals for its iterative process to converge. It has become easier to a certain extent for decomposing multichannel SEMG signals which are originally difficult to process since the CKC method was introduced into the field of SEMG decomposition. Other multichannel signal processing methods have also been tested with high density sEMG decomposition, including traditional template matching, independent component analysis, higher order cumulants, and correlation measurement, but most of these methods have been limited to relatively low muscle contraction levels.

In view of the existing facts, it is hard to decompose complex superposition sEMG signals. Moreover, the decomposing procedure is also a bit cumbersome. For example, it usually needs multiple steps to build a correlation vector between the IPT and the measurements in the past. While in this article, it only needs an iterative procedure to form the correlation vector. A method based on measurement correlation and linear minimum mean square error (MC-LMMSE) was developed in this study to decompose multichannel sEMG signals. The MC-LMMSE method is firstly used to reconstruct a matrix correlated with the measurement matrix. Then, it gradually and iteratively increases the correlation between an optimized vector and a reconstructed matrix until a satisfactory innervation pulse train (IPT) is obtained. The performance of the MC-LMMSE method was assessed with both simulated and experimental sEMG signals. The results demonstrated that the MC-LMMSE method can successfully extract more MUs and reconstruct IPTs with a higher true positive rate (TPR) than the GCKC method, even from complex superposition signals.

2. Materials and Methods

2.1. Data Model

Multichannel sEMG signals can be modelled as a linear-time-invariant multi-input-multi-output system [22] if the muscle contraction is maintained at a constant force level. This system can be represented by the matrix form as follows [18]:where is the M measurements, is the nth sample of the jth measurement, stands for a vector of zero-mean white noise, denotes a mixing matrix which consists of all of the channel responses (the jth source in sEMG signals appearing in the ith measurement) of samples, and is an extended form of the N sources can be described as .

2.2. Method of LMMSE

Given a vector form whose probability density function (PDF) is unknown, the linear estimator of a variable θ related to the X statistics can be written as follows:

Choose the weighting coefficients a_n’s to minimize the Bayesian mean square error (MSE):where the resultant estimator is termed the linear minimum mean square error (LMMSE) estimator [23]. Substituting (2) in (3), then it becomes

Differentiating and setting this equal to zero,

Then, it produces

Substituting (6) in (4), then it becomes

Let , and it has

Because , it has

Equation (9) can be maximized by taking the gradient:and setting it to zero, which results in

Substituting (6) and (11) into (2) produces

If the means of θ and X are zero, then

For multichannel sEMG signals, is the innervation pulse train (IPT) that needs to be estimated, X is the measured multichannel sEMG signal, and C_Xθ is a parameter that needs to be calculated. It has been pointed out in [18] that all firing times of MU need to be known in advance to calculate C_Xθ, which can be written aswhere set contains all firing times of the same MU and is the extended form of the measured signal X. The LMMSE estimator can be obtained after substituting (14) in (13). In fact, it is very hard to know MU firing times beforehand. In view of this, a method that is able to identify complete or most of firing time of MU was proposed; therefore, we can achieve the results or approach results of LMMSE.

2.3. Measurement Matrix Autocorrelation

Multiplied by a 1 × M vector from both sides, (1) becomes

The ith IPT in can be calculated with (16) if (suppose the value of the element is 1, and all other values are 0), and the noise term is negligible.

In practice, it is difficult and even impossible to find such a vector if the mixing matrix G is unknown, but can still be satisfactorily reconstructed as long as one of the elements in is far greater than others.

The similarity between vectors A and B can be evaluated as follows [24]:where denotes the inner product and denotes the norm. The shapes of the MUAPs generated by the same MU should have a certain degree of similarity when the isometric muscle contraction is held at a constant force. Therefore, the inner product of two vectors which are associated with different time instants fired by the same MU, should be relatively large. This property provides the possibility to estimate the IPTs of MUs with the following equation:where is the value of the estimated innervation pulse train at the sample time , denotes the number of sample times in each channel, and is a 1 × M vector. If has a strong correlation with the measurement vectors associated with the time instants fired by a particular MU, the firing pattern of this MU will be easily observed in . The vector , then, increases the values in at time instants when this MU is firing and decreases other values at time instants when it is not. The following average form [18] can be used as to achieve such a purpose:where denotes the time instants fired by the particular MU, is the series of measurement vectors associated with , and is the number of elements in . An ideal will have a stronger correlation with all the measurement vectors contained in and, in this case, due to the average result, should be larger than other values in and can be easily observed. It is difficult, however, to find a satisfactory vector , as the firing pattern of any MU is unknown in practice. As a result, it is necessary to develop an advanced approach to better estimate in order to successfully reconstruct the IPTs.

2.4. Measurement Correlation Based on LMMSE (MC-LMMSE)

An iterative algorithm based on LMMSE is developed in the proposed MC-LMMSE method to gradually optimize the vector f in order to achieve a better IPT reconstruction. Assuming is a matrix which has a certain column correlation with , then the IPT estimation equation can be rewritten aswhere the vector plays the same role as the aforementioned vector . Replace with in (19) and the vector can be rewritten as

In this article, the matrix in the MC-LMMSE method is reconstructed from unitary matrices obtained from the singular value decomposition (SVD) of the measurement matrix (see Step 1 below). Other matrices can also be selected. The high column correlation of the matrix helps the MC-LMMSE increase the values of in (20) at the time instants fired by the same MU. Hence, the influence of noise on its IPT estimation results is significantly suppressed.

An initial vector f will first be formed from any time instants fired by an MU. The MC-LMMSE method will then be implemented by following the steps listed below to make f approximate the ideal vector in (21) and to reconstruct future IPTs with high accuracy. The schematic outline of the MC-LMMSE is shown in Figure 1.(1)Decompose the matrix using SVD, where T denotes the transpose, and estimate the matrix .(2)Randomly select sEMG signals from a few channels and denote each channel signal by ; calculate the Teager energy operator [25] of , , and set a threshold (thre); identify all the time instants in which satisfy to form .(3)Choose , and then estimate an IPT from each time instant in according to (20).(4)Identify (the subscript k denotes the kth iteration) time instants, corresponding to the highest peaks for each initial IPT , where (A and C are constants greater than or equal to zero, where in most instances ), and then replace with . A new IPT will be obtained by substituting f into (20). The vector f will be gradually improved by repeating this iterative process until ( is a rough estimate number of firing times in each IPT) at which point the final IPT will be obtained.(5)Classify all the IPTs into groups for each specific MU.

Note that after substituting in Step (1) into (20), it is similar to CKC method [18]. Both of them are correlation method in essence. However, it is helpful for simplifying the decomposition expression by using (20) and understanding the distinguishing feature of these correlation methods to decompose sEMG.

The MC-LMMSE and classic CKC have some similarities, which include that (1) they directly estimate IPTs from measurement matrix without involving calculation of unknown mixing matrix G, (2) they all need to select some vectors of measurement corresponding to discharged time instants. However, there are also differences between them that lead to different results (see the following section of results). In MC-LMMSE method, a new way to reconstruct matrix correlated with the measurement matrix was proposed (20). Then, sEMG signals can be decomposed by using the reconstructed matrix. In this article, the SVD method was used to reconstruct the correlation matrix. The measurement matrix itself can also be directly used as the correlation matrix (see the previous section of measurement matrix autocorrelation). In addition, other effective ones such as the measurement matrix transformed by FastICA [26] can also be used as the correlation matrix. Hence, it may further obtain better results if the correlation matrix is properly selected in future. Another difference comparing with CKC is the utilization of iterative technique in Step (4) which can achieve more precise IPTs, and more number of MU firing time is an improved iteration method of CKC. Because it can find more number of time instants discharged by one MU in the process of gradually and iteratively calculating vector in (21) in terms of the characteristics of SEMG and the algorithm. Therefore, the quality of can be improved a lot when comparing with classic CKC and GCKC methods.

2.5. Simulated Signals

2.5.1. Simulated Signals Generated by Random Mixing Matrices

Ten sources were assumed and the IPTs were randomly generated in the simulation with a mean IPI of 100 samples. The lengths of the IPTs were set to 20,000 samples where , was uniformly distributed over the interval [−10, 10]. The zero-mean mixing matrix G was also randomly generated with a length _ij of 10 samples. The number of measurements was set to 25 and the number of delayed repetitions of each original measurement was set to 9. Therefore, the number of extended IPTs was increased to 190 with 250 measurements. Gaussian zero-mean noises were added to each signal with different signal-to-noise ratios (SNRs) of −10 dB, −5 dB, 0 dB, 5 dB, and 10 dB. The measurement matrix autocorrelation method did not need to reconstruct the matrix in Step 1, while in Step 3 and in Step 4 were replaced by and , respectively. The number of channels and threshold value in Step 1 were set to 10 and 0.45 ( denotes the maximum absolute value of in the Step 2), respectively, when the MC-LMMSE method was implemented. The number of iterations to estimate [19] and the number of main decomposition loops were set to 40 and 500, respectively, when the method of GCKC was implemented. The scalar function was taken in (9) from [19]. An IPT was selected as real when its TPR was greater than 75%.

2.5.2. Simulated Signals Generated by Gaussian Function [27]

The extracellular single fiber action potential (SFAP) was depicted by the sum of three basic Gaussian functions [28].where t is time, is the amplitude factor, is the bandwidth, and is the position of the center of the peak. With this equation, one may approximate a particular triphasic action potential waveform with considerable accuracy by adjusting and . Each fiber is assumed to be parallel to the skin surface, so the shape of the SFAP detected by the electrodes is considered to be a function of the physiological parameters, such as the fiber location within a 3-dimensional Cartesian coordinate system, and the muscle fiber conduction velocity. and in (22) were depicted aswhere stands for the vertical fiber depth below the surface of the skin, represents the center position along the fiber in the z-x plane, x is the fiber center position in the z-x plane perpendicular to the direction, and is the conduction velocity of muscle fiber. The MUAP shapes detected by different electrodes were depicted as the summation of the SFAP shapes contained in the MUs. The MUAP trains were then generated by the convolution of the MUAP shapes with their corresponding firing times. Finally, the composite sEMG signals were modelled as linear summations of the MUAP trains. The characteristics of the MUAP, such as the amplitude distribution, shape, and duration, were determined by the morphological properties of the active muscle fibers contained within corresponding MUs. The SEMG signals can be simulated with considerable similarity by adjusting the parameters of Gaussian functions according to the characteristics of real SEMG signals. In this article, the SEMG signals are just roughly simulated. However, it can still demonstrate the basic characteristics of SEMG.

The depths of the centers for all measured MUs were uniformly distributed from 1 mm to 6 mm. A random number of fibers (uniformly distributed between 30 mm and 70 mm) were assumed in active MUs. All semifiber lengths were set to 50 mm, and the tendon and endplate positions of the fibers were uniformly distributed in the range of ±5 mm. The conduction velocities of active MUs were set to 4.0 m/s, the firing rates of the MUs were normally distributed with the mean, and standard deviation of 20 ± 5 Hz and 60 active MUs were assumed in total. The starting times of MUs were chosen from 10 ms to 200 ms. A 16 × 16 electrode-array grid with a 3 mm interelectrode distance in both directions was employed for recording the sEMG signals. This grid center was placed at the center of the muscle and the signals were sampled with frequency of 2,000 Hz. The numbers of fibers, position of the active MUs, and discharging patterns were all randomly generated. The signals were also corrupted by additive Gaussian zero-mean noise with SNR of 20 dB as shown in Figure 2. The number of delayed repetitions of each original measurement was set to 9 [18, 19].

2.6. Experimental Signals

The experimental sEMG signals were collected from the first dorsal interosseous (FDI) muscles of three adult subjects. The procedures were approved by the Institutional Review Board of Northwestern University (Chicago, USA), and all three subjects gave their written consent before the experiment. Subjects were seated upright in a mobile Biodex chair (Biodex, Shirley, NY). A standard 6 degrees-of-freedom load cell (ATI Inc, Apex, NC) setup was used along with standard procedures for minimizing spurious force contributions from unrecorded muscles as described in [29] to accurately record the isometric contraction force of the FDI muscle during index finger abduction. SEMG signals were recorded from the FDI muscle using a flexible 2-dimensional 64-channel surface electrode array (8 × 8 array with the electrode probe diameter of 1.2 mm, and the center-to-center probe distance of 4 mm) (TMS International BV, The Netherlands) [30]. The skin of the tested muscle was carefully prepared and the electrode array was attached to the FDI muscle with a double adhesive sticker and further secured with medical tapes [29]. The maximum voluntary contraction (MVC) was first measured. Each subject was then asked to generate an isometric contraction force of the FDI muscle at the different contraction levels of 2 N, 4 N, 6 N, and 8 N. Multiple trials were performed with one force level being recorded for each trial. The subject was asked to maintain the force as stable as possible for up to 15 s. A Refa amplifier (TMS International BV, The Netherlands) was used to record sEMG signals. The signals were sampled at 2 kHz with a bandpass filter set at 10–500 Hz. The number of delayed repetitions of each original measurement was set to 9 [18].

2.7. Validation

For simulated signals, the parameter TPR and MR defined in (25) and (26) are used to further validate the accuracy of sEMG signal decomposition algorithm, and defined in (26):where TP is the number of correctly identified firing times of pulses in the reconstructed IPT, FP is the number of misplaced discharges, and FN stands for the number of unidentified firing times of pulses in the IPT. For the simulated signals generated by the Gaussian function, the firing time tolerance was set to ±1 sample. Therefore, each identified firing time was considered as true if it was detected within ±0.5 ms (sampling frequency of 2,000 samples/s) from its actual position along the signal. The value defined in (25) was averaged over 10 trials for all identified IPTs. For simulated signals generated by random mixing matrices, the time tolerance was set to 0. The value defined in (25) in this case was also averaged over 10 trials for all identified IPTs.

For experimental signals, to validate the accuracy of MC-LMMSE algorithm, the “two-source” technique, in which all 64 channels of the electrode array were divided into two independent groups with equal number of channels, was used as an alternative to using intramuscular EMG together with surface EMG [16, 17]. The coincident rate of the firing times of the MUs, which are decomposed from both channel groups using the MC-LMMSE algorithm, were calculated, and a high coincident rate was taken to suggest a favourable performance of the algorithm.

3. Results and Discussion

3.1. Simulated Signals

3.1.1. Tests on Signals Generated by Random Mixing Matrices

Ten trials were conducted to test the performance of the proposed MC-LMMSE method in decomposing the sEMG signals simulated by random mixing matrices and the results were averaged over the 10 trials. The number of reconstructed IPTs, corresponding TPR and MR achieved by the measurement matrix autocorrelation, GCKC and the MC-LMMSE method at different SNRs are presented in Table 1. Results show that the measurement matrix autocorrelation method could not completely reconstruct the IPTs even with a high SNR of 10 dB. The GCKC method only reconstructed an average of 5 IPTs when SNR was set to −10 dB. The MC-LMMSE method reconstructed all the 10 IPTs successfully with the high TPRs at all tested SNR levels (−10 dB to 10 dB) and the TPR maintained over 92% even in severely noisy environments (SNR = −10 dB). Results demonstrate that the MC-LMMSE method offers superior performance to the measurement matrix autocorrelation and GCKC methods of sEMG decomposition. In addition, a parameter called pulse-to-noise-ratio (PNR) [31] was also utilized to evaluate the performance of MC-LMMSE method. The average PNR was 12.37 dB and infinite, respectively, when SNR was set at −10 dB and greater than 0 dB.

3.1.2. Tests on Signals Generated by Gaussian Function

The GCKC and MC-LMMSE methods were employed to decompose the sEMG signals generated by a Gaussian function. On average, 26 IPTs were reconstructed by the GCKC method with a TPR of 92.67% and MR of 4.26%; while 53 IPTs were reconstructed by the MC-LMMSE method with a TPR of 97.89% and MR of 1.93% (Figures 3 and 4). The average PNR of MC-LMMSE was 27.39 dB. Figure 3(a) shows the MUAP shapes of one MU, detected by the 16 × 16 electrode array, which were estimated using the spike-triggered averaging method [32]. The innervation zone of the MU and the propagation of MUAPs can also be clearly observed.

(a)

(b)

Figure 3(b) shows the 53 IPTs reconstructed from the signals. The firing times of each extracted MU are indicated by an assigned label at top of the signal in Figure 4. Thirty-five MUs can be correctly identified from this channel and the challenge caused by overlapped action potentials appears to be solved by the proposed MC-LMMSE method. The parameters used in the MC-LMMSE and GCKC methods for this test are the same as those used in Test 1, except that the number of main decomposition loops in the GCKC method was set to 5,000.

3.2. Tests on Experimental Signals

Figure 5(a) shows the force profile and the 16 IPTs identified from the sEMG signals of the FDI muscles by using the MC-LMMSE method. These sEMG recordings were taken during an isometric constant force contraction at 10% of the maximum voluntary contraction (MVC).

(a)

(b)

(c)

Figure 5

Results obtained from first dorsal interosseous (FDI) muscle. (a) MU discharge patterns with the force profile identified from the FDI muscle during an isometric constant force contraction at 10% MVC (2 N (Subject A). Each vertical line indicates a MU discharge at a given time instant. ((b) Top panel) the sum of identified MUAP trains (grey lines) compared to the raw sEMG signal (black lines) in one selected channel from the first dorsal interosseous (FDI) muscles during an isometric constant force contraction at 10% MVC (Subject A). ((b) Middle panel) an expended view of the top panel. ((b) Bottom panel) the residual (grey lines) compared to the raw sEMG (black lines) after the subtraction of the reconstructed MUAP trains. (c) MU firing patterns identified from Group 1 (black lines) and Group 2 (grey lines). All 64 channel signals were divided into 2 independent groups, with the even numbered columns selected as one group and the odd numbered columns as the other group.

It can be seen that the firing rate of MUs changes with the fluctuation of the contraction force; Figure 5(b) compares the summation of the identified MUAP trains and their residuals respective to the original sEMG signals, where the signal-to-interference ratio (SIR) [33] between the sum of identified MUAP trains and raw sEMG signal was 59.73%. Figure 6 shows the mean and standard deviation of discharge rates of the extracted 16 MUs from FDI muscles. It can been seen from the figure that the average discharge rates of these extracted 16 MUs range from 7.55 ± 5.1 to 18.1 ± 7.1 pulses/second. Different MUs correspond to different average discharge rate patterns, which are monotonically increasing. Considering the individual differences of the physiological characteristics [34], these values may differ slightly from the previous reported results; however, overall they are similar [34, 35].

The results achieved by the GCKC and MC-LMMSE methods are shown in Table 2. It can be seen that the MC-LMMSE method extracted more MUs than the GCKC method, especially in the cases of high force contraction. The parameters used in the MC-LMMSE and GCKC methods in this test are the same as in Test 1. The performance of the MC-LMMSE method with the experimental electrode array sEMG was further investigated by using the “two sources” method. All of the 64 channel signals recorded at different contraction force levels were divided into 2 independent groups, each with 32 channels. SEMG signals recorded from channels with even column numbers were selected to form Group 1, while signals recorded from channels with odd column numbers were selected to form Group 2. The proposed MC-LMMSE method was applied to each of the groups for sEMG decomposition, and the numbers of MUs extracted from all the channel signals, signals in Group 1 and signals in Group 2, were compared (Table 2 and Figure 5(c)). It can be seen that, overall, the number of extracted MUs decreases as the number of EMG channels decreases. This trend becomes more remarkable in cases where a higher force of contraction was applied. It can also be seen that results achieved from the two independent groups share over 84% of the commonly extracted MUs and show over 90% of the same firing times for the common MUs.

4. Discussion

One important concept to decompose high density array signals like SEMG is proposed in this article. There are two important steps for decomposing signals which lead to its superior performance compared to other decomposition methods. One is the appropriate selection of the matrix which is correlated to the measurement matrix; the other one is the estimation of the reconstructed IPTs with the iterative optimization process presented in Step 4. Both steps are critical in achieving favourable decomposition results. In fact, in [18] can also be considered as a matrix correlated with the measurement matrix. In addition to the mentioned correlated matrix in this article, other matrices have also been found that can decompose sEMG signals. The decomposition results are likely to improve in near future. However, like other decomposition methods, the MC-LMMSE method also has some limitations. For example, there ought to be at least thousands of samples in sEMG signals, otherwise if the length of signals is too short, it will be difficult to obtain satisfactory results. It can be seen from the results of the simulation data that this MC-LMMSE method requires a larger number of detected electrodes to get better results. But only 64 electrodes were used to record the real sEMG signals in this article. Hence, if hundreds of electrodes could be employed to record the real signals, there is hope that a larger number of MUs could be extracted, and the allowable force of muscle contraction could also become larger.

The matrix in the MC-LMMSE method is constructed with a high level of column correlation from the unitary matrices obtained using the SVD of the measurement matrix. This high column correlation is able to help the MC-LMMSE suppress the influence of noise, as the correlation between vector and the other vectors from associated with the firing times of the same MU is further enhanced by the iterative optimization procedure in Step 4. (In fact, the results obtained by MC-LMMSE can better approach the LMMSE estimator when compared with CKC which is derived from LMMSE estimator. Please refer to [18] for further understanding why MC-LMMSE method can get such results.) Therefore, both the employment of a SVD of the measurement matrix and the iterative optimization procedure in the MC-LMMSE contribute to the improvement of the decomposition performance when compared to the other methods tested in this paper (Figure 3, Tables 1 and 2). The time instants in each iteration step corresponding to the highest peaks in are usually the firing times of a particular MU, making it possible to employ such an optimizing approach to improve the vector f. Both the decomposition method presented in this study and CKC method are based on high density surface EMG recordings; however, the MC-LMMSE method employs a different approach for IPT estimation. It differs from CKC that (21) was adopted in the proposed MC-LMMSE algorithm to gradually optimize the vector f and give it a stronger correlation with the different vectors from matrix associated with the firing time instants of a particular MU. Instead, Equation (20) is utilized in the MC-LMMSE method to estimate the IPTs, where is reconstructed by the unitary matrices obtained through the SVD of , and the vector f is obtained by an iterative optimization procedure. The final IPTs can then be obtained by substituting f into (20). CKC and GCKC are the two typical sEMG signal decomposition methods. Moreover, GCKC can get better results compared with CKC [19], hence we chose GCKC method as a comparison here. The following relevant published articles have little improvement in performance and many of them are related applications for decomposition. It should be noted that although the results obtained by MC-LMMSE seem to be superior to CKC, it had better be further confirmed by an independent research team.

IPTs can be relatively easily reconstructed from sEMG signals with a low degree of MUAP superposition as long as in Step 4 is similar to (Table 1). However, it will be difficult to satisfactorily reconstruct the IPTs from sEMG signals with a high degree of MUAP superposition in cases where is small. For these scenarios, more iteration steps will be needed to optimize the vector f to adequately reconstruct the IPTs (Figures 3(b) and 4).

In order to evaluate the performance of the proposed MC-LMMSE method in experimental sEMG decomposition, the 64 sEMG channels were divided into two independent groups with equal numbers of channels in each. The “two sources” method was employed to compare MUs which were independently decomposed from the two groups of the sEMG signals. Comparison results in Table 2 confirm the high stability, efficiency, and accuracy of the MC-LMMSE method in experimental sEMG decomposition. As a correlation method, the MC-LMMSE method requires a relatively large number of electrodes to achieve good decomposition results; a reduction in the number of electrodes leads to a reduction in the amount of correlation information, which will affect the number of reconstructed IPTs in decomposition results. Consequently, it is necessary to increase the number of recording electrodes under the premise that the amount of information in the sEMG is fully provided if a large number of extracted MUs are desired, particularly in cases of relatively high muscle contraction levels (Table 2).

The major challenges in sEMG decomposition can be summarized as follows [17]: (1) the occurrence of large amounts of superposition between the action potentials from different MUs; (2) the changes in shapes of the different action potentials contained in every MUAP train; (3) high degree of similarity in action potential shapes between different MUAP trains. Those challenges can be overcome to some extent by using the proposed MC-LMMSE method. IPTs can still be reconstructed with high accuracy even if they have a high degree of MUAP superposition (Figures 3(b) and 4). The shapes of the action potentials in MUAP trains may change during the isometric muscle contractions as a result of the changes in conduction velocity (e.g., caused by muscle fatigue) or movement of the electrode, making the decomposition task more challenging. However, even if a large degree of change in MUAP shape occurs quickly, the MC-LMMSE method can still be used to reconstruct the IPTs efficiently and accurately by increasing the iterations in Step 4 or by dividing the signal recordings into short epochs, which could then be considered stationary an absent of shape changes. It is unlikely that the shapes of the action potentials contained in different MUAP trains are similar across all observed channels and, as a result, the correlation between the measurements vectors associated with different MUs can be neglected. Note that both MC-LMMSE and CKC build on the low probability of different MUs to share the exact firing time [18]. MU synchronization does affect the decomposition performance. How much it affects in detail depends on the level of the synchronization and the complexity of the sEMG signal, such as the degree of MUAP waveform superposition, the amount of noise, and so on. In fact, it is extremely difficult to encounter a very high synchronous rate when decomposing real sEMG signals. The formula of probability of synchronization rate was given in (12) of [18]. It can also be seen from the formula that the probability of synchronization is very low. The literature [36] also shows that in the case of very high synchronization rate, it can still be decomposed well by GCKC method. CKC, GCKC, and MC-LMMSE are all based on LMMSE and correlation methods. If GCKC and CKC can do that, the method in this study can also do that.

5. Conclusions

In summary, a new MC-LMMSE method was developed for multichannel sEMG decomposition based on the principle that the measurement vectors associated with the firing times of a single MU have a certain degree of similarity. The MC-LMMSE method gradually and iteratively increases the correlation between the optimized vectors and the reconstructed matrix to better decompose complex sEMG signals. The superior performance of the MC-LMMSE method was demonstrated with both simulated and experimental electrode array sEMG signals. The results show that, in each case, the MC-LMMSE method can extract a relatively large number of MUs with strong robustness to noise and excellent accuracy.

Data Availability

The data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.

Acknowledgments

The authors would like to thank Mr. Thomas Potter from the University of Houston for editing the manuscript. This work was supported by National Natural Science Foundation of China (51677171), Zhejiang Provincial Natural Science Foundation of China (LY17C100001), Department of Education of Zhejiang Province (Y201533132), China Scholarship Council, and Ningbo Natural Science Foundation of China (2017A610218).

References

R. Merletti and P. A. Parker, Electromyography: Physiology, Engineering, and Non-Invasive Applications, John Wiley & Sons, Hoboken, NJ, USA, 2004.
C. J. De Luca, P. J. Foley, and Z. Erim, “Motor unit control properties in constant-force isometric contractions,” Journal of Neurophysiology, vol. 76, no. 3, pp. 1503–1516, 1996.
View at: Publisher Site | Google Scholar
V. Ruonala, E. Pekkonen, O. Airaksinen, M. Kankaanpää, P. A. Karjalainen, and S. M. Rissanen, “Levodopa induced changes in electromyographic patterns in patients with advanced Parkinson’s disease,” Frontiers in Neurology, vol. 9, no. 35, 2018.
View at: Publisher Site | Google Scholar
J. W. Yoo, D. R. Lee, Y. J. Cha, and S. H. You, “Augmented effects of EMG biofeedback interfaced with virtual reality on neuromuscular control and movement coordination during reaching in children with cerebral palsy,” Neurorehabilitation, vol. 40, no. 2, pp. 175–185, 2017.
View at: Publisher Site | Google Scholar
T. A. Caires, L. F. R. M. Fernandes, L. J. Patrizzi et al., “Immediate effect of mental practice with and without mirror therapy on muscle activation in hemiparetic stroke patients,” Journal of Bodywork and Movement Therapies, vol. 21, no. 4, pp. 1024–1027, 2017.
View at: Publisher Site | Google Scholar
B. Barth, K. Mayer, U. Strehl, A. J. Fallgatter, and A.-C. Ehlis, “EMG biofeedback training in adult attention deficit/hyperactivity disorder: an active (control) training?” Behavioural Brain Research, vol. 329, pp. 58–66, 2017.
View at: Publisher Site | Google Scholar
M. J. Zwarts and D. F. Stegeman, “Multichannel surface EMG: basic aspects and clinical utility,” Muscle & Nerve, vol. 28, no. 1, pp. 1–17, 2003.
View at: Publisher Site | Google Scholar
E. A. Clancy and N. Hogan, “Probability density of the surface electromyogram and its relation to amplitude detectors,” IEEE Transactions on Biomedical Engineering, vol. 46, no. 6, pp. 730–739, 1999.
View at: Publisher Site | Google Scholar
N. Fallentin, K. Jorgensen, and E. B. Simonsen, “Motor unit recruitment during prolonged isometric contractions,” European Journal of Applied Physiology and Occupational Physiology, vol. 67, no. 4, pp. 335–341, 1993.
View at: Publisher Site | Google Scholar
R. Merletti, M. Knaflitz, and C. J. De Luca, “Myoelectric manifestations of fatigue in voluntary and electrically elicited contractions,” Journal of Applied Physiology, vol. 69, no. 5, pp. 1810–1820, 1990.
View at: Publisher Site | Google Scholar
J. L.F. Weytjens and D. Vansteenberghe, “The effects of motor unit synchronization on the power spectrum of the electromyogram,” Biological Cybernetics, vol. 51, no. 2, pp. 71–77, 1984.
View at: Publisher Site | Google Scholar
P. Zhou and W. Z. Rymer, “MUAP number estimates in surface EMG: Template-matching methods and their performance boundaries,” Annals of Biomedical Engineering, vol. 32, no. 7, pp. 1007–1015, 2004.
View at: Publisher Site | Google Scholar
S. Karimimehr, H. R. Marateb, S. Muceli, M. Mansourian, M. A. Mañanas, and D. Farina, “A real-time method for decoding the neural drive to muscles using single-channel intra-muscular EMG recordings,” International Journal of Neural Systems, vol. 27, no. 6, Article ID 1750025, 2017.
View at: Publisher Site | Google Scholar
J. Roussel, P. Ravier, M. Haritopoulos, D. Farina, and O. Buttelli, “Decomposition of multi-channel intramuscular EMG signals by cyclostationary-based blind source separation,” IEEE Transactions on Neural Systems and Rehabilitation Engineering, vol. 25, no. 11, pp. 2035–2045, 2017.
View at: Publisher Site | Google Scholar
X. Ren, C. Zhang, X. Li, G. Yang, T. Potter, and Y. Zhang, “Intramuscular EMG decomposition basing on MUAPs detection and superposition resolution,” Frontiers in Neurology, vol. 9, no. 2, 2018.
View at: Publisher Site | Google Scholar
M. Chen, X. Zhang, X. Chen, and P. Zhou, “Automatic implementation of progressive FastICA peel-off for high density surface EMG decomposition,” IEEE Transactions on Neural Systems and Rehabilitation Engineering, vol. 26, no. 1, pp. 144–152, 2018.
View at: Publisher Site | Google Scholar
S. H. Nawab, S. S. Chang, and C. J. De Luca, “High-yield decomposition of surface EMG signals,” Clinical Neurophysiology, vol. 121, no. 10, pp. 1602–1615, 2010.
View at: Publisher Site | Google Scholar
A. Holobar and D. Zazula, “Multichannel blind source separation using convolution kernel compensation,” IEEE Transactions on Signal Processing, vol. 55, no. 9, pp. 4487–4496, 2007.
View at: Publisher Site | Google Scholar
A. Holobar and D. Zazula, “Gradient convolution kernel compensation applied to surface electromyograms,” in Proceedings of the Independent Component Analysis and Signal Separation, vol. 4666, pp. 617–624, London, UK, September 2007.
View at: Google Scholar
M. Chen and P. Zhou, “A novel framework based on FastICA for high density surface EMG decomposition,” IEEE Transactions on Neural Systems and Rehabilitation Engineering, vol. 24, no. 1, pp. 117–127, 2016.
View at: Publisher Site | Google Scholar
I. Gligorijevic, J. P. Van Dijk, B. Mijovic, S. Van Huffel, J. H. Blok, and M. De Vos, “A new and fast approach towards sEMG decomposition,” Medical & Biological Engineering & Computing, vol. 51, no. 5, pp. 593–605, 2013.
View at: Publisher Site | Google Scholar
D. Zazula and A. Holobar, “An approach to surface EMG decomposition based on higher-order cumulants,” Computer Methods and Programs in Biomedicine, vol. 80, pp. S51–S60, 2005.
View at: Publisher Site | Google Scholar
S. M. Key, Fundamentals of Statistical Signal Processing: Estimation Theory, Prentice-Hall Int., Englewood Cliffs, NJ, USA, 1993.
Y. Ma, S. H. Lao, E. Takikawa, and M. Kawade, “Discriminant analysis in correlation similarity measure space,” in Proceedings of the 24th International Conference on Machine Learning, pp. 577–584, Corvallis, OR, USA, June 2007.
View at: Google Scholar
H. Teager and S. Teager, Evidence for Nonlinear Sound Production Mechanisms in the Vocal Tract: Speech Production and Speech Modelling, Springer, Berlin, Germany, 1990.
Y. Ning, J. Li, and S. Zhu, “Research on decomposition of surface EMG signals based on FastICA and channel correlation,” Space Medicine & Medical Engineering, vol. 30, no. 3, pp. 191–197, 2017, in Chinese.
View at: Google Scholar
Y. Ning and Y. Zhang, “A new approach for multi-channel surface EMG signal simulation,” Biomedical Engineering Letters, vol. 7, no. 1, pp. 45–53, 2017.
View at: Publisher Site | Google Scholar
J. Clark and R. Plonsey, “A mathematical evaluation of the core conductor model,” Biophysical Journal, vol. 6, no. 1, pp. 95–112, 1966.
View at: Publisher Site | Google Scholar
X. Li, A. Suresh, P. Zhou, and W. Z. Rymer, “Alterations in the peak amplitude distribution of the surface electromyogram poststroke,” IEEE Transactions on Biomedical Engineering, vol. 60, no. 3, pp. 845–852, 2013.
View at: Publisher Site | Google Scholar
P. Zhou, N. L. Suresh, and W. Z. Rymer, “Surface electromyogram analysis of the direction of isometric torque generation by the first dorsal interosseous muscle,” Journal of Neural Engineering, vol. 8, no. 3, Article ID 036028, 2011.
View at: Publisher Site | Google Scholar
A. Holobar, M. A. Minetto, and D. Farina, “Accurate identification of motor unit discharge patterns from high-density surface EMG and validation with a novel signal-based performance metric,” Journal of Neural Engineering, vol. 11, no. 1, p. 016008, 2014.
View at: Publisher Site | Google Scholar
K. G. Keenan, D. Farina, R. Merletti, and R. M. Enoka, “Amplitude cancellation reduces the size of motor unit potentials averaged from the surface EMG,” Journal of Applied Physiology, vol. 100, no. 6, pp. 1928–1937, 2006.
View at: Publisher Site | Google Scholar
A. Holobar, M. A. Minetto, A. Botter, F. Negro, and D. Farina, “Experimental analysis of accuracy in the identification of motor unit spike trains from high-density surface EMG,” IEEE Transactions on Neural Systems and Rehabilitation Engineering, vol. 18, no. 3, pp. 221–229, 2010.
View at: Publisher Site | Google Scholar
P. Zhou and W. Z. Rymer, “Can standard surface EMG processing parameters be used to estimate motor unit global firing rate?” Journal of Neural Engineering, vol. 1, no. 2, pp. 99–110, 2004.
View at: Publisher Site | Google Scholar
C. J. De Luca, R. S. LeFever, M. P. McCue, and A. P. Xenakis, “Behaviour of human motor units in different muscles during linearly varying contractions,” Journal of Physiology, vol. 329, no. 1, pp. 113–128, 1982.
View at: Publisher Site | Google Scholar
A. Holobar, V. Glaser, J. A Gallego, J. L Dideriksen, and J. Farina, “Non-invasive characterization of motor unit behaviour in pathological tremor,” Journal of Neural Engineering, vol. 9, no. 5, Article ID 056011, 2012.
View at: Publisher Site | Google Scholar

Copyright

Copyright © 2018 Yong Ning et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

PDF Download Citation

Download other formats

Order printed copies

Views

2086

Downloads

1057

Citations