Journal of Healthcare Engineering

Volume 2018, Article ID 2347589, 12 pages

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

## Multichannel Surface EMG Decomposition Based on Measurement Correlation and LMMSE

^{1}School of Automation and Electrical Engineering, Zhejiang University of Science and Technology, Hangzhou 310023, China^{2}School of Computing, University of Portsmouth, Portsmouth PO1 3HE, UK^{3}China Coal Research Institute, Beijing 100013, China^{4}Ningbo University of Technology, Ningbo 315211, China

Correspondence should be addressed to Yong Ning; moc.621@6180gnoygnin

Received 19 January 2018; Revised 9 April 2018; Accepted 14 May 2018; Published 28 June 2018

Academic Editor: Chengyu Liu

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.

#### 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 *n*th sample of the *j*th measurement, stands for a vector of zero-mean white noise, denotes a mixing matrix which consists of all of the channel responses (the *j*th source in sEMG signals appearing in the *i*th 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 *i*th 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 *k*th 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.