Performance of a Novel Muscle Synergy Approach for Continuous Motion Estimation on Untrained Motion

When applying continuous motion estimation (CME) model based on sEMG to human-robot system, it is inevitable to encounter scenarios in which the motions performed by the user are different from the motions in the training stage of the model. It has been demonstrated that the prediction accuracy of the currently effective approaches on untrained motions will be significantly reduced. Therefore, we proposed a novel CME method by introducing muscle synergy as feature to achieve better prediction accuracy on untrained motion tasks. Specifically, deep non-smooth NMF (Deep-nsNMF) was firstly introduced on synergy extraction to improve the efficiency of synergy decomposition. After obtaining the activation primitives from various training motions, we proposed a redundancy classification algorithm (RC) to identify shared and task-specific synergies, optimizing the original redundancy segmentation algorithm (RS). NARX neural network was set as the regression model for training. Finally, the model was tested on prediction tasks of eight untrained motions. The prediction accuracy of the proposed method was found to perform better than using time-domain feature as input of the network. Through Deep-nsNMF with RS, the highest accuracy reached 99.7%. Deep-nsNMF with RC performed similarly well and its stability was relatively higher among different motions and subjects. Limitation of the approach is that the problem of positive correlation between the prediction error and the absolute value of real angle remains to be further addressed. Generally, this research demonstrates the potential for CME models to perform well in complex scenarios, providing the feasibility of dedicating CME to real-world applications.


I. INTRODUCTION
T HE surface electromyography (sEMG) signals have been used for motion intent recognition, which in turn has been applied to human-robot interaction systems of rehabilitation robots, allowing the robots to actively assist the patients with motor function impairments and broadening the range of applications.Feature extraction of sEMG combined with machine learning has enabled the classification of motion patterns [1], [2], [3], etc.Along with discrete motion classification problem, another important application is continuous motion estimation.Continuous motion estimation tasks based on sEMG use time-or frequency-domain features extracted from filtered sEMG signals to build a model of mapping relationship with kinematic parameters [4], [5], [6].The most significant advantage of using sEMG for human kinematics prediction is that sEMG signals are generated earlier than the actual motion, thus enabling the predicted angle signal to be applied in the active-compliance control strategy of Human-robot Interaction system [7].However, this method is only applicable to the ideal situation, such as healthy, non-fatigued muscle with a single trained motion pattern.When predicting any untrained motion, or when fatigue and compensatory phenomena occur within the muscles, the prediction accuracy undergoes a significant degradation, or worse, it might identify a completely wrong changing direction of the angular signal.Among the current methods using time-domain features and neural networks, the approach with the highest performance for untrained action prediction tasks still has not exceeded 90% accuracy, which is difficult to put into practical applications [8].
Muscle synergy analysis is a method that decomposes the sEMG signal to each precise muscle by potential origin.It has been used to investigate the synergy of healthy subjects' hands [9], shoulder [10], etc., as well as the effect of synergy quantity on gait in non-healthy subjects, such as subjects with neurologic injuries [11].These studies clarify the principles of motion from muscular perspective.Analogically, introducing muscle synergy into continuous motion estimation could potentially allow a more intuitive interpretation of the method from the neuromuscular theoretical perspective.Furthermore, it has been demonstrated that synergies are less likely to be affected relative to time-domain characteristics under different motion conditions, between subjects, and even when subjects perform actions in non-ideal muscle states [12], [13], [14].It could be inferred that by including analysis of the potential characteristics of each muscle, the generalization ability of the prediction model is optimized when subjects perform multiple untrained motion patterns.
Correlative research in neurology has demonstrated that the generation of motion is based on the combination of discrete motion primitives [15].The musculoskeletal system has large number of degrees of freedom.Many hypotheses proposed that the process of muscle motion controlled by CNS was driven by multiple combinations of discrete components, and it has been verified in the research of d'Avella et al. [16] that those components are muscle synergies with specific time-varying characteristics, i.e., coordinated activations of muscle groups.Non-negative matrix factorization (NMF) is a part-based algorithm that was firstly proposed by Lee and Seung [17].Due to the restriction that the result of NMF decomposition is positive, it has been widely used in the study of muscle synergy since it is suitable for the fact that muscle can only be activated positively [18], [19], [20], [21].Cheung et al. [22] further proposed a method to minimized the error in EMG reconstruction by introducing an exponentially distributed type of noise term to the loss function.Besides NMF, there are several algorithms that have been used in muscle synergy, such as independent component analysis, principal component analysis (PCA), autoencoder [23], etc.In the study by Rabbi et al. [24], NMF has been proven to perform the best on synergy extraction tasks among these algorithms.However, muscle synergy using the NMF decomposition method still suffers from problems such as redundancy in the decomposition results, and local minima in the algorithm itself causing the decomposition results not to be closest to the original matrix [25].Le Roux et al. [26] proposed an NMF-based deep network architecture to disassemble the iterations and parameters of standard NMF.Yu et al. [27] proposed deep non-smooth NMF (Deep-nsNMF) by introducing sparseness and deep architecture into NMF, and obtained higher accuracy on face clustering compared to NMF and non-smooth nonnegative matrix factorization (ns-NMF) [28], [29], [30].Yet there have not been any research utilizing Deep-nsNMF on sEMG signal for muscle synergy analysis to our knowledge.
Moreover, due to the two types of synergies, namely shared and task-specific muscle synergies [15], which are included in the synergies extracted from specific training motion, the results extracted from all training motions can also contain redundancy if they are directly input into the angle prediction model as a single feature matrix.Therefore, further optimization of the results from NMF is also necessary.In the research on the redundancy of sEMG signals, Carmona-Saez et al. [28] proposed a correlation-based redundancy segmentation (RS) algorithm to extract an irredundant sEMG subvector.
Qin et al. [32] utilized the cosine similarity value to achieve RS of muscle synergies extracted from time-domain features.However, current RS has not considered the redundancy in synergy weights.Inspired by the two mentioned studies, we proposed the idea of redundancy classification on synergy weight based on correlation coefficient.
Muscle synergy was initially researched on animals.To apply synergy decomposition in human-related tasks, d'Avella et al. [33] implemented the capture of complex spatiotemporal features for reaching specific muscle patterns on the human upper limb by combining a small number of components based on the method from previous frog experiments, verifying that synergy decomposition enables the reconstruction of motion patterns through the downscaling of the raw EMG signals.On this basis, researchers have further applied the results of decomposition to discrete intention recognition or continuous motion estimation problems by using synergies as input to regression models such as machine learning [34] and neural networks [32], [35], [36], [37], [38].In previous research, the data used in the analysis and prediction process were generally sourced from a limited and homogeneous set of motions.Although the results obtained were promising, these synergies might only be effective within each motion type and it is unknown whether the models can represent a generalized control strategy for all types of motions.In the study of d'Avella et al. [33], the muscle patterns of multiple movements under different loads and postures were reconstructed by the synergies and the reconstruction accuracy proved the generalizability of synergies.Huang et al. [34] proposed a feature space decomposition method to improve the accuracy of muscle synergies on recognizing movement pattern switching.Torres-Oviedo et al. [39] demonstrated that at least five synergies are required to adequately reconstruct the decoding of muscle forces in walking locomotion of cats.To the best of our knowledge, there are no studies that address the generalization ability of the method by applying muscle synergy to continuous motion estimation.
To address these issues, this paper presented analysis and training based on goal-oriented motion tasks, and completed prediction and performance evaluation on a new set of life-scene motion tasks to ensure that the method has generalization ability and robustness.The Deep-nsNMF algorithm was first introduced into the synergistic extraction of EMG signals to optimize the local minima problem that exists in NMF algorithm, and to obtain a decomposition result that is closer to the original data.An optimized redundancy classification (RC) algorithm was proposed to analyze the shared muscle synergy and task-specific muscle synergy between different motion modes by classifying the synergy weights according to a redundancy threshold.We used a nonlinear autoregressive neural network with external input (NARX) as the multi-joint angle prediction model, and the feasibility of the optimized approach was verified by the prediction tasks on untrained motions.

II. METHODOLOGY
The continuous motion estimation method proposed in this paper consists of three phases, as shown in Fig. 1.In the first Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
stage, raw sEMG is preprocessed into time-domain feature and the Deep-nsNMF method is used to obtain the activation primitive matrix from all trained motions.Then the RC and RS algorithm is used to downscale the matrix, to reduce redundancy and to identify the type of the synergies.With the non-redundant excitation primitives, the matrix is input into the NARX network and trained to fit the 12 channels of angle motion, then the model is tested on untrained motion to verify the generalization ability of the method.

A. Dataset and Signal Preprocessing
The sEMG and angle signals used in this paper were selected from the Ninapro Database2 (DB2, https://ninapro.hevs.ch/instructions/DB2.html) [40], which consists of 3 modules.Hand kinematics were acquired using a 22-sensor CyberGlove II dataglove with a sampling frequency of 2000 Hz, and then the data were transferred to a PC at a frequency of slightly lower than 25 Hz using a Bluetooth tunnel serial port.sEMG was acquired using 12 Trigno Wireless electrodes (Delsys, Inc, www.delsys.com),which were sampled at 2000 Hz and relayed to the PC using a standard USB connection.Data provided by both sensors was correlated to a precise timestamp through Windows performance counters.The first 8 channels of sEMG consist of signals from isometric electrodes placing around the forearm, and channels 9 to 12 were from targeted muscles: the muscle Flexor Digitorum Superficialis, Extensor Digitorum Superficialis, Biceps Brachii, and Triceps Brachii, respectively.The placement of the angle sensor was referenced in Fig 2(a).
The first eight types of motion in the first module, namely Exercise B, was selected for the model training.These movements consist of individual and collaborative stretching and bending motions of the five fingers and they basically cover the daily motion patterns of the hand, as shown in Fig 2(b).The remaining nine motion patterns are mainly wrist movements, which resulted in similar amplitude patterns for the majority of the angular signal channels, therefore this section of data was not used to ensure training effectiveness and low time cost.
The motions in the second module, Exercise C, were chosen for the untrained motion prediction tasks.Exercise C contains grasping and functional motions that mimic interaction with physical objects in daily life, and in this paper, eight motions were chosen as shown in Fig. 2(c).
Through observing the variation of all 22 angle signal channels over the training motion, ten channels with the most distinct pattern changes, and two channels with the smoothest changes were selected to ensure that the method proposed in this paper can adapt to different characteristics of the motion when applied to other situations.The selected angle channels are 1, 2, 3, 5, 8, 11, 12, 13, 16, 19, 20 and 21.The sEMG signal acquired from the electrodes consists of power line interferences, therefore the noise was shielded using a 50 Hz Hampel filter.The power spectral density (PSD) represents the distribution of signal power in the frequencydomain [41].The PSD of the raw sEMG signal was calculated and plotted against frequency, and it was observed to be highest at range [10,160] Hz.A 4th order Butterworth bandpass filter was used for noise reduction of the EMG signal and the cutoff frequency was set to [10,160] Hz.Similarly, PSD was calculated for angular signals and a 4th order Butterworth low-pass filter with a cutoff frequency of 20 Hz was used for noise reduction.Unipolar signals were obtained by full wave rectification of the filtered sEMG signals, and the mean absolute value (MAV) was extracted as time-domain feature.

B. Synergy Extraction Algorithm
The sEMG signal matrix can be factorized into two matrix, muscle synergy weights and synergy excitation primitives, i.e.W and H, using NMF, where x n ] T is a matrix consists of n channels and t sample point in time.If the number of synergy patterns factorized from the original sEMG matrix was set as k, then w i in W = [w 1 , w 2 , . . ., w k ] represents the muscle activation pattern for all channels in the i-th synergy, H = [h 1 , h 2 , . . ., h k ] T symbolizes the projection of the signal's features in the time-domain into the specific synergy pattern.
In order to obtain the WH that is closest to X , the Euclidean distance was used as the lost function.Iterations of H and W are next performed according to the multiplicative update rules [42], where • represents Hadamard product.Due to the possibility of local minima in the loss function E during decomposition [43], different results are obtained for the same set of sEMG signals when undergoing multiple factorizations for the same defined value of k.The NMF can also converge below the given k value.Therefore, the NMF algorithm needs to be further optimized.
The sparse NMF algorithm (s-NMF) is implemented by introducing sparsity weight [26].Based on s-NMF, ns-NMF improved the sparseness degree by adding a smoothing matrix S to Eq. ( 1), The ns-NMF gains a sparser decomposition than NMF with no constraints.To ensure that the matrix error of before and after the factorization is sufficiently low and to select the most appropriate value of k, the degree of conformity is evaluated using variance accounted for (VAF), For the purpose of ensuring that what is applied to the analysis and prediction for different subjects and different motion patterns is the synergy in the extracted decomposition results that is closest to the real activation pattern of muscle groups, so that the decomposition results reach the highest and stable Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.VAF value, Deep-nsNMF was introduced by decomposing X into k layers: The optimized algorithm yields a VAF that is positively correlated with the value of k.To ensure the factorization accuracy while minimizing the dimensionality of H, the value of k at which the VAF was first higher than 0.95 was chosen as the number of synergies.

C. Redundant Synergy Extraction Algorithm
After the sEMG signals of all motion types were decomposed, the number of synergies k for which the VAF exceeded the given threshold for the first time was selected to obtain the corresponding W and H, respectively.The most fitting k value for each motion types are different, and the synergy excitation primitives extracted for the n-th action type can be expressed as, Combine all sub-vectors h n i into a total synergy excitation primitive H syn = [H 1 , H 2 , . . ., H 8 ] T , i.e.
The presence of shared synergies in H syn leads to redundant primitives for decomposition, which can increase the computational workload of the training process.Redundancy Segmentation Algorithm [28], [32] enables the extraction of non-redundant features for muscle synergy by calculating the correlation coefficients between excitation primitives after synergy factorization.Firstly, the correlation coefficients between the excitation primitive vectors are calculated, where C a,b is the correlation coefficient of the a-th and b-th excitation primitive, h a and h b are the mean value, S a and S b are the standard deviation, t is the total number of samples.
Then, the threshold θ th is set.When the value of C a,b exceeds θ th , it indicates that the ath and bth synergies are relatively linearly correlated, whereas a value of C a,b below θ th represents that they are two mutually independent synergies.Therefore, H syn can be segmented into two parts, Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.syn .If there exists at least one value of the correlation coefficients corresponding to h i that is greater than or equal to θ th , then h i is considered to have redundant synergies and is stored in H 2 syn ; if all of the correlation coefficients corresponding to h i are less than θ th , then h i is considered an independent synergy and is stored in H 1 syn .The redundancy segmentation is completed after iterating through all the excitation primitive vectors in H n .
2) Optimal Redundancy Classification Algorithm: In the two matrixes obtained by processing the surface EMG signals using the factorization algorithm, the synergy weights W represents the degree of activation of each muscle in each synergy, and the synergy excitation primitives H represents the time-domain variations of the degree of activation corresponding to each synergy.In the previous algorithm, the information contained within W was ignored by RS directly on H. Therefore, a novel W -H redundancy classification algorithm is proposed.
Firstly, all w n i are combined into a general synergistic excitation primitive matrix W syn , Next, RC is performed on W syn .By iterating all synergy weights and calculating the correlation coefficients, the w i of which all C a,b are lower than the threshold θ th , are regarded as new synergy types, and those that exist at least one C a,b greater than the threshold θ th are regarded as redundant synergies and categorized by the maximum C a,b value.The optimized algorithm yields not just two fuzzy categories, redundant and non-redundant, but specifically all the classifications and all the items contained in each class.The pseudocode of the W classification process in RC method is shown in Algorithm 1.
After completing the classification of W , redundant segmentation of H is conducted within each class respectively, and all irredundant H are combined into a single feature matrix, which will be further used as an input to the regression model.

D. Design of the NARX Multi-Joint Angle Prediction Model
After obtaining the excitation primitive matrix with RC, the results cannot be directly used as inputs to the regression model because the number of synergies extracted is different for each subject.Furthermore, the trained model will be used on the prediction of a single type of motion, so number of input matrix dimension would be significantly different.Therefore, by calculating the mean value of the k values extracted from all trained motions and setting the rounded mean value as the number of neurons in the input layer of the network, H is downscaled using PCA to achieve the requirement of the network structure setting while minimizing the loss of information [44].The NARX neural network was utilized as the regression model, the hidden layer was set to 15 neurons based on experimental experience, and the output layer corresponded to 12 angle channels.The number of delays is set to 2. The Sigmoid function is used as the activation function and Bayesian Regularization is used as the training algorithm.

E. Evaluation Criteria
Normalized Root Mean Square Error (NRMSE), Pearson's correlation coefficient (R), Coefficient of Determination (COD) and Concordance Correlation Coefficient (CCC) were selected as evaluation criteria.After acquiring the predicted angle data for each channel using the four approaches on the untrained dataset, the NRMSE, R, COD and CCC between the true and predicted angles with CI 95% were provided.

A. Performance of Prediction on Untrained Motion
After completing feature extraction, synergy factorization, and redundancy extraction (see Supplementary Material), we used four approaches, namely the time-domain feature MAV, the synergy excitation primitives extracted from NMF, the excitation primitives after Deep-nsNMF with RS, and the excitation primitives after Deep-nsNMF with RC, respectively as input matrix for the NARX network, and the corresponding angle data in Exercise B was used as the output.After training the network, it was used to predict each of the eight sets of motions in Exercise C. The prediction results were collected and analyzed.
For a specific subject, through visualizing the true and predicted angle as shown in Fig. 3, it was found that when using only the time-domain feature and the synergy excitation primitives from NMF as inputs, the two methods were able to obtain correct predictions for the time periods when the motions were more stable, but there were large errors in the predictions for both methods where the angle showed high variability during the cycle.There were even errors in several channels that gave completely opposite predictions of the trend Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.  of angle, such as the predicted angle in channel 11 using NMF in Fig. 3.The prediction results obtained by the two redundant feature extraction algorithms after synergy factorization using the Deep-nsNMF algorithm are overall better than the first two approaches and the difference between the two results by using RS and RC is not significant.The results of SPM analysis with repeated measures ANOVA for all channels are shown in Fig. 4, where it can be seen that significant differences occurred in the intermediate stages of motion cycle on channels 1,5,7,9 and 12.
The indicators NRMSE and R of the four methods on the eight motions of EC were analyzed independently, as shown in Fig. 5.All subjects except subject 4 performed in a similar manner on each motion.The errors in the predictions using the four methods were higher than overall for motion 2 (index finger extension grasp) and motion 8 (power disk grasp), while the correlation coefficients between the predicted and true values for motion 2 and motion 7 (parallel extension grasp) were noticeably lower than average.The prediction error for Mot 3 to Mot 7 were relatively lower and the difference was non-significant.The prediction error NRMSE of Mot 2 and Mot 8 of MAV and NMF methods were significantly different compared to NRMSE of other untrained motions, while the differences between different motions of the two methods using Deep-nsNMF were not significant.It could be considered that Deep-nsNMF method has improved prediction robustness for the types of untrained motions that highly influenced the performance of the previous methods, MAV and NMF.Analysis of the specific motion patterns reveals that the knife-grip index finger extension grasp in motion 2 and the  parallel extension grasp of a book in motion 7 both contain the feature that the main work-doing element of the motion is the extended fingers, The disc grasping in motion 8 also contains the extension of most joints on the fingers and the contraction of the terminal joints.However, all the motions in the training set EB were captured without interaction with external objects, and none of the designed stretching motions were able to perform the motor function of grasping.Therefore, it was hypothesized that the synergies extracted from the training set did not include the task-specific synergies in motion 2,7 and 8 of the predicting motion set, resulted in inferior training performance compared to the other motions.Fig. 6 illustrated the Bland-Altman analysis of the four methods, taking nine values in one motion cycle for each subject, and the signal residuals were calculated for all untrained motion types, with different colored scatters in the figure representing different motion types.The results of the Bland-Altman analysis for all methods were demonstrated in Table I.It can be seen that the MAV method has the largest proportion of outliers and the NMF method has the largest deviation of the mean line from zero axes.While the Deep-nsNMF with RC method has the smallest mean line to zero axes deviation and the smallest percentage of points out of the consistency limit.
Observations of the overall distributions of the two indicators revealed that the range of the distributions of the error and correlation coefficients decreased with the improvement of the method.The prediction accuracy of Deep-nsNMF with RS and RC are similar, but Deep-nsNMF with RC had a generally smaller margin of error than with RS, which indicated that the proposed method showed a more stable performance among different situations.Overall, the results demonstrated the generalization ability of the proposed method to adapt to differences between subjects, and the improvement of the capability to extract valid information from the sEMG signal.Fig. 7 illustrated the overall performance of all subjects' prediction using the four methods.Table II and Table III showed the confidence interval of NRMSE, R, COD and CCC within all subjects and prediction motion types.It could be seen that there was a significant reduction in NRMSE using the Deep-nsNMF method compared to both MAV and NMF.When R was used as the evaluation criteria of prediction accuracy, there was a significant improvement using the DS and DC methods compared to MAV.When COD and CCC were used for evaluation, there was a significant improvement in the DS method compared to both MAV and NMF, and a significant increase in the DC method compared to the MAV method.It could be seen that the two methods using Deep-nsNMF generally performed better than using MAV and NMF.The prediction performance using NMF is better relative to when only time-domain feature was used as input.
The two redundant feature extraction algorithms gave similar prediction results.The prediction accuracy on untrained motions of Deep-nsNMF with RS was slightly higher than that of Deep-nsNMF with RC if compared carefully.But the general range of NRMSE and R of the proposed method was both smaller than Deep-nsNMF with RS, which proved the optimized stable feature of Deep-nsNMF with RC algorithm.

IV. DISCUSSION
Based on the continuous motion estimation problem using time-domain features of sEMG signals, prediction of multi-joint angles on untrained motion dataset was achieved by introducing synergy factorization algorithm.There are studies that optimize NMF, the most used synergy factorization algorithm, from the principle of the algorithm itself, such as previous work [22] that allowed the algorithm to better adapt to different types of noise through model selection criteria.However, noise is one of the many factors affecting the accuracy of the NMF algorithm, and the instability of the decomposition results caused by local minima to the NMF algorithm has a much more serious impact on the subsequent analysis.Comparing to other optimized NMF, The Deep-nsNMF method proposed in study [27] produces more localized and less overlapping features, and also solved the problem on local minima.Therefore, the potential of using Deep-nsNMF on muscle synergy extraction was validated in our research.As comparison, MAV was used as a time-domain feature and input into the NARX neural network to construct the mapping relationship between the feature and the angle.The synergy factorization was performed using the NMF algorithm, and the k values corresponding to the optimal VAFs were obtained by multiple times of factorizations.In order to make the decomposed result more stable and closer to the original matrix, the Deep-nsNMF was introduced as the optimized factorization algorithm.It was experimentally verified that this method could obtain the positive correlation between VAF and k, and the VAF after decomposition was basically higher than 50%, which revealed the improvement of decomposition calculation efficiency.Through analyzing the decomposition results (see Supplementary Material), it was found that the VAF values of Deep-nsNMF were generally higher than those of NMF decomposition, while the k values of the former were selected higher than those of the latter.This is due to the fact that the decomposition performance of Deep-nsNMF was more stable than that of NMF, and the overall threshold of VAF was also higher, therefore the standard set for VAF was also higher for Deep-nsNMF than NMF.In addition, the VAF-k analysis of NMF reveals that it was difficult to achieve stable VAF value higher than 0.95, and the optimal case was to achieve a VAF of about 0.90.Therefore, it can only be selected according to the optimal case of VAF, which resulted in k-value of Deep-nsNMF being generally higher than NMF.
Due to the existence of shared synergies and task-specific synergies among different motions and even the single motion, there was redundancy in the decomposed excitation primitives.A matrix containing all the independent excitation primitives was obtained using redundant feature extraction algorithms, and input to the regression model.In study [32], redundant feature of activation coefficient vectors was extracted through similarity, but the feature of W between different motions was not taken into account.Considering the theory of redundancy segmentation proposed in [28], we proposed the strategy of analyzing the synergy weights before synergy classification.Through comparing the dominant muscle channels of synergistic weights between different motions and the similarity of activation patterns of each primitive, it was found that the same synergy would present different or similar activation patterns in different action modes.Since the RS algorithm only considers the redundancy of synergistic excitation primitives, this paper optimized the algorithm by firstly classifying the synergy weights for redundancy, and then segmenting the redundancy within each synergy weight pattern.Using the proposed RC algorithm enables intuitive analysis of the composition of shared and task-specific synergies in each motion.The irredundant synergy matrix was then fed into the regression model for training on all Exercise B type motion data.
In contrast to the theory of muscle synergy, the uncontrolled manifolds (UCM) is another common theoretical assumption used for analyzing principle of neural motor control [45], [46].Different from muscle synergies, UCM uses the stability of a particular state as a criterion to determine which elemental variables have control over the system, by analyzing the redundancy between the degrees of freedom of joint space relative to the end-effector.On the other hand, muscle synergies have the fixed channel number of the sEMG signal and assume that the motion is controlled by a variety of specific inter-channel synergies.Related study [47] also explored the possibility of using the UCM hypothesis to identify synergies in a postural task.This study demonstrated a significant correlation between muscle activation patterns and the steady state of the center of pressure, but the application of the UCM to sEMG signal is complex and unable to account for all other potential performance variables.Muscle synergy factorization downgrades the sEMG matrix and uses the feature vectors obtained from the factorization algorithm directly for learning the mapping relationship with the angle of motion, which is more suitable for continuous motion prediction problems.
Through conducting Bland-Altman analysis of the residual signal, it was revealed that the real angle of the 6 th channel was larger than the values of all other channels, and that the prediction error for that channel was correspondingly larger.When we excluded the 6 th channel and analyzed 100 samples drawn from the remaining 11 channels, the mean difference ± SD was −0.2575 ± 9.5709.A one-sample t-test on the residuals yielded P-value of 0.4943 (Null hypothesis: the mean of the residual is zero).Therefore, the bias between the real and predicted angle of the proposed method can be considered not statistically significant.The comprehensive performance of the models trained using the above four methods for predicting the multi-joint angles of the untrained motion set Exercise C is shown in Table II.The overall performance of the two methods using Deep-nsNMF was similar, and also significantly better than the results using the MAV and NMF.The average prediction accuracy when using the RS algorithm was improved by 7.7% over the accuracy using only timedomain feature, and the highest prediction accuracy of 96.2% was achieved when performing prediction on Motion 1 of EC.The accuracy when using RC algorithm was improved by 7.0% over using MAV, and again reaching the highest prediction accuracy of 96.3% at motion 1 of EC.
Although the proposed method presented a significant improvement in performance under different testing criteria, there remain limitations that have not been addressed.As shown in Fig. 5 for the NRMSE after prediction using the Deep-nsNMF method, the difference between various untrained motions was not significant, but the prediction errors for Mot 2 and Mot 8 were still the highest and second highest, just like the performance of the MAV and NMF methods.The authors concluded that an increase in the extent of the difference between the training and testing motions would decrease the predictive performance of the proposed method, which would potentially lead to the failure of the method.However, due to the limitation of the dataset, it is currently difficult to quantitatively analyze the differences between the training and testing motions.Therefore, further research is planned to quantify the difference between trained and untrained motions using a set of specially designed motion types.The Bland-Altman plot indicates that the prediction error increased as the absolute value of ground truth angle increased.This issue can be observed intuitively in the fitting plots of predicted and real angle, and is found to be prevalent in several recent studies on continuous motion estimation [32], [48], [49], [50].The essence of this problem is that the prediction error tends to maximize at the peaks and troughs of the real angle, which the method proposed in this paper has not yet solved.Integrating the Bland-Altman analysis to unravel the relationship between prediction error and the real angle presents a promising avenue for future research.Since this paper focused on the optimization and validation of the prediction methodology, only offline analysis has been conducted, yet the proposed method has not been deployed on a real-time prediction system.For systems with the function of real-time prediction, it is generally recommended that the prediction time of the algorithm is less than 200ms [51].Besides, the time cost of the training phase of the algorithm, the effect on the synergy feature by muscle fatigue, and the applicability of the method to subjects with motion disorders also need to be taken into account.In future studies, the construction of a real-time prediction system and the deployment of online training models will be completed to provide more in-depth study of above-mentioned problems and to validate the method proposed in this paper.
Generally, the experiments demonstrated that the introduction of muscle synergy and Deep-nsNMF significantly improved the prediction accuracy of untrained motions.The prediction performance with the proposed RC algorithm was slightly lower than that of RS algorithm, which might be due to the inadequate number of experimental subjects.In future studies, we will increase the number and age range of subjects to enhance the feasibility of the method.

V. CONCLUSION
In this paper, we aimed to address the problem of multi-joint continuous motion estimation using sEMG signals from untrained motions.Through the theory of muscle synergy, Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
we introduced the Deep-nsNMF algorithm as an improvement of the original decomposition method NMF.To reduce the dimension of input matrix and to realize the extraction of task-specific synergy and shared synergy, we used the RS algorithm and proposed a RC algorithm based on it.The feasibility of the algorithm proposed in this paper was verified by conducting prediction experiments on the untrained motion dataset.The introduction of Deep-nsNMF was proved to have significant improvement on the generalization of the model, and the proposed RC algorithm enhanced the rationalization and interpretability of the synergy process.In future research, we will design the training motions in a more reasonable way to ensure that a limited number of motions can be used to train a model with appropriate prediction accuracies for a variety of untrained motions, and to improve the generalization ability of the model while increasing the efficiency of the training process.

Fig. 2 .
Fig. 2. Information on the dataset collection and motion selection.(a) The placement of the angle sensor.(b) Eight motions selected from EB for model training.(c) Eight motions selected from EC for untrained motion prediction tasks.

Fig. 4 .
Fig. 4. The real angle (black), the predicted angle by MAV-NARX (yellow), NMF-NARX (green), Deep-nsNMF-RS (red) and Deep-nsNMF-RC (blue) respectively, and F-values of the SPM analysis with repeated measures ANOVA for all subjects shown in time-domain, in which dashed red lines indicate P = 0.05 level and grey shaded areas indicate regions with statistically significant differences.

Fig. 6 .
Fig. 6.The bland-Altman analysis of all tested untrained motions by 4 discussed methods.Different colors of the dots represent different types of motions, particularly from motion type 1-8 corresponding to the color is: red, green, blue, yellow, magenta, cyan, gray, brown.The black line is the regression line of all samples.
Algorithm 1 Redundancy Classification Algorithm Initialize: W_syn = the weights extracted from all motion pattems; W_seg = the first weight vector in W_syn; W_num = {1}; Calculate Correlation of W and update: for i = 2 to (the number of weight vectors in W_syn) C_tem = [] for j = 1 to (the number of weight vectors in W_seg) Compute the C a,b of W_ seg(j, 1) and W_ syn(i,:); Store C a,b in C_tem; end if all in C_tem < θ th add W_syn(i,:) in series vertically into W_seg; add i vertically in series into W_num; else max _C = the row number of W_seg with the biggest C a,b ; add i horizontally in line(max _C) of W_num; end end h 1 was firstly stored in H 1 syn , by taking the remaining h i in H 1 syn in turn, C a,b values are computed with respect to all the vectors in H 1

TABLE I THE
BLAND-ALTMAN ANALYSIS OF 4 METHODS ON 6912 SAMPLES