High-Density Surface EMG Decomposition by Combining Iterative Convolution Kernel Compensation With an Energy-Specific Peel-off Strategy

Objective- This study aims to develop a novel framework for high-density surface electromyography (HD-sEMG) signal decomposition with superior decomposition yield and accuracy, especially for low-energy MUs. Methods- An iterative convolution kernel compensation-peel off (ICKC-P) framework is proposed, which consists of three steps: decomposition of the motor units (MUs) with relatively large energy by using the iterative convolution kernel compensation (ICKC) method and extraction of low-energy MUs with a Post-Processor and novel ‘peel-off’ strategy. Results- The performance of the proposed framework was evaluated by both simulated and experimental HD-sEMG signals. Our simulation results demonstrated that, with 120 simulated MUs, the proposed framework extracts more MUs compared to K-means convolutional kernel compensation (KmCKC) approach across six noise levels. And the proposed ‘peel-off’ strategy estimates more accurate MUAP waveforms at six noise levels than the ‘peel-off’ strategy proposed in the progressive FastICA peel-off (PFP) framework. For the experimental sEMG signals recorded from biceps brachii, an average of 16.1 ±3.4 MUs were identified from each contraction, while only 10.0 ± 2.8 MUs were acquired by the KmCKC method. Conclusion- The high yield and accuracy of MUs decomposed from simulated and experimental HD-sEMG signals demonstrate the superiority of the proposed framework in decomposing low-energy MUs compared to existing methods for HD-sEMG signal decomposition. Significance- The proposed framework enables us to construct a more representative motor unit pool, consequently enhancing our understanding pertaining to various neuropathological conditions and providing invaluable information for the diagnosis and treatment of neuromuscular disorders and motor neuron diseases.


I. INTRODUCTION
D URING the past decades, surface electromyography (sEMG) has received growing attention due to its noninvasive nature, high yield of motor units (MUs), and application to high contraction levels.By taking advantage of the significant advancements in sEMG signal collection and processing, this technique plays a critical role in understanding the neurophysiology of the neuromuscular system, as well as in the diagnosis and treatment of motor neuron diseases and neuromuscular disorders.The sEMG signal is the summation of motor unit action potentials (MUAPs) from a larger number of active motor units (MUs) within the electrode recording range.Decomposition, a technique that can break down the composite sEMG signal into constituent pulse trains, is indispensable for the study of MU firing rate and MUAP morphology.
Tremendous efforts have been made by researchers and various decomposition algorithms have been proposed.Previously, pattern recognition was one of the most adopted sEMG decomposition techniques.By combining wavelet transform and adaptive resonance theory network classification, Gazzoni et al. proposed an automatic decomposition algorithm for MUAP waveforms detection and recognition in 2004 [1].Nawab et al. extended their knowledge-based artificial intelligence framework [2], which was originally developed for intramuscular EMG decomposition to sEMG but was only limited to low contraction force [3].Later, this approach was refined by Nawab et al. in 2010 to make it applicable for high contraction force levels [4].Despite these efforts, the low-pass filter effect caused by skin and subcutaneous fat separating the muscle from the sEMG electrodes can lead to severe signal superposition, MUAP duration increment, and shape difference reduction.All these factors working together pose significant technical challenges to the pattern-recognition type of sEMG signal decomposition.
In this study, a novel framework for high-density sEMG (HD-sEMG) decomposition was developed by combining the iterative convolution kernel compensation (ICKC) method with a post-processor and a new 'peel-off' strategy.This approach allows for the extraction of MUs with small MUAP energy thus achieving high yield and accurate HD-sEMG decomposition results.

II. METHODS AND MATERIALS
As shown in Fig. 1, the proposed iterative decomposition algorithm consists of three parts: ICKC, Post-Processor (including a strategy for screening IPT with different energies and a strict ICKC method), and Peel-off.The ICKC was firstly applied to the filtered sEMG data to acquire IPTs.Then a screening strategy was adopted to divide the IPTs into highenergy IPTs and low-energy IPTs by considering the pulseto-noise-ratio (PNR), coefficient of variation of inter-spike intervals (COV I S I ), coefficient of variation of the amplitude of spikes (COV A ) and maximum cross-correlation coefficient (ρ).For the high-energy IPTs, the firing instants of MUs can be easily acquired by peak detection; while for the lowenergy IPTs, the strict iterative cycle CKC (SICKC) method was applied to precisely extract the firing instants of the MU.Finally, a new 'peel-off' strategy was adopted to subtract the MUAPs of decomposed MUs from the filtered HD-sEMG Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
data.This entire process will be iteratively performed until no reliable spike trains can be extracted.

A. Data Model
The sEMG can be described as a multiple-input-multipleoutput system: where y(n) is an observation, n is the discrete time, H is the mixing matrix, which consists of all of the channel responses h i j = [h i j (0), . . ., h i j (P − 1)] (the j th MU in surface EMG signals appearing in the i th channel) of length P samples.P is the duration of the action potentials.
. .]T , and w(n) is a noise vector.Then the global pulse train activity index γ (n) reflecting the activity of all active MUs can be calculated by: where C yy stands for the correlation matrix of extended observation y(n)= y 1 (n) , y 1 (n − 1) , . . ., y 1 (n − L) , . . .T .L is the expansion factor and is used to increase the conditionality of the mixing process by adding delayed versions of each observation.This is done to increase the ratio between the number of observations and sources, which can improve the accuracy and quality of the mixing process.Empirically, L was set to 10.Based on the linear minimum mean square error (LMMSE) theory, the j th MU's pulse train can be estimated as [12] and [13]: where Ĉt j y can be estimated by: where Ψ j denotes a set of firing times of j th MU, and car d(Ψ j ) denotes the cardinal number of Ψ j .

B. Iterative CKC Method
In this study, an iterative CKC (ICKC) Method was first adopted to identify the initial IPTs.It includes two main parts: rough iteration and precise iteration.For each rough iteration, two time instants (adjacent or not) corresponding to the two maximum values in the pulse train t j (n) were detected and added to the set Ψ j for the update of t j (n) based on equations (3) and (4).For each precise iteration, two non-adjacent time instants (at least 30 sample points apart) corresponding to the two peak values in the pulse train t j (n) were detected to update the t j (n).Each extraction of IPT requires h1 rough iterations and h2 precise iterations.Empirically, for a 10s signal, h1 and h2 were set to 30-40 and 40-50, respectively.

C. Post-Processor
If only a cursory screening of initial IPTs is performed based on PNR and COV I S I metrics (PNR > 36 dB and COV I S I < 0.3) only high-energy MUs will be retained, while low-energy MUs will usually be discarded.The primary reason lies in that low-energy MUs have many common firing instants with high-energy MUs [31].During the decomposition process, if too many common firing instants are added to the time set of low-energy MUs, the subsequent decomposition will ignore the firing instants of low-energy MUs and extract only the high-energy MU instants.Consequently, the final IPT obtained will be of extremely low quality and will then be discarded.To solve this problem, this study proposes a Post-Processor for specifically identifying the low-energy IPTs.The basic steps of post-processing are as follows.In the first step, the N IPTs obtained using the ICKC are first classified into high-energy IPTs and low-energy IPTs.The classification procedure is performed by considering the PNR, COV I S I , COV A , and maximum cross-correlation coefficient (ρ) of each IPT.According to previous research result [22], the thresholds for PNR, COV I S I , COV A , ρ are set as 36 dB (c 1 ), 0.3 (c 2 ), 0.15 (c 3 ) and 0.3 (c 4 ), respectively.In the second part, the SICKC method is proposed to accurately extract the IPTs with relatively low MUAP energy.Specifically, the SICKC method consists of two parts, de-duplication, and precise iteration.For the de-duplication part, we assume that M spike trains were extracted from M high-energy IPTs by peak detection and grouped them as Ψ h j = n j1 , . . ., n jr , j = 1, . . ., M; while G spike trains were extracted from G low-energy IPTs and grouped them as Ψ li = {n i1 , . . ., n ir } , i = 1, . . ., G.For any low-energy IPTs t i (n), if the number of common time instants between Ψ h j and Ψ li exceeds 30% of the cardinal number of Ψ li , the common time instants will be deleted from Ψ li to make the information pertaining to the low-energy MUs stand out; otherwise, no further procedure is needed.Then the ICKC is applied to refine the LIPTs.The 2 nd part of the pseudocode presented in Fig. 2 clearly shows the flowchart of the postprocessor algorithm.

D. New 'Peel-off' Strategy
To prevent FastICA from converging to local solutions, Chen and Zhou [22], [23] peeled off the MUAP waveforms constructed by the least squares method from the sEMG signal.However, due to the estimation errors and alignment errors of the MUAP, some noise might be introduced and accumulates in the residual signal after each peeling-off.With the progress of peeling off, the low-energy MUs' information might be submerged by noise, leading to significant difficulty in the extraction of the low-energy MUs.Inspired by the previous work [31], a new MUAP waveform construction strategy was proposed based on the source deflation theory for the peeloff purpose.We assume that k MUs are activated at the time instant i, the set of indexes corresponding to the activated MU is set to S i = s i1 , . . ., s ik , then the t(n) at the time instant i in Eq. ( 1) can be expressed as: where b j (i) denotes the information of the s i j th MU firing at the time instant i and it has the same dimension as t (i).For example, if the s i1 th MU is active at the time instant i, b We first optimize b j (i) in the source space by establishing the corresponding objective function and constraints.Then, by linking the source space and the signal space through Eq. ( 3), we can subtly transform the objective function and its constraints for optimizing b j (i) into the objective function and its constraints for optimizing C j (i), where C j (i) is the MUAP of the s i j th MU at the time instant i.By considering the orthogonality and sparsity of the firing MUs, b j (i) can be optimized by minimizing ( 6) and (7) with the constraint of (8): where E j (i) denotes the proportion of the s i j th MU account to the signal energy at the time instant i.The objective function (6) as well as the constraint (8) can be transformed from the source space to the signal space based on equations ( 3) and ( 5): where C yy stands for the correlation matrix of the observation y(n).The objective function (7) can be approximately expressed as: where the symbol * represents the Hadamard product and b j (i) * b l (i) 1 can be expressed as the inner product of b l (i) and b j (i): Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
To further enhance the sparsity of the objective function (11), we remove the term b j (i) * b j (i) 1 from k l=1 b j (i) * b l (i) 1 .Thus, the objective function ( 11) can be approximated as: In a typical setup, the initial value of C j (i) is estimated using the spike trigger averaging technique, where the value is set to C t j y , as calculated in Eq. ( 4) using the j and y (n).The E j (i) was set to 0.9 y T (i)C −1 yy C t j y .It should be noted that in Eq. ( 13), the MU used for optimizing C j (i) only includes the MU activated at the time instant i except for the s i j th MU.However, for actual optimization, all identified MUs within the whole observed sEMG signal, excluding the s i j th MU, can be utilized to enhance the accuracy of C j (i)'s optimization.Users can select the MU utilized to optimize C j (i) based on their desired balance between computational efficiency and precision.Additionally, low-energy MUAP waveforms, often considered physiological noise, can cause an amplitude shift in the estimated MUAP waveform.Therefore, amplitude alignment of the MUAP waveform is critical (e.g.subtracting the mean value of the MUAP waveform from the MUAP waveform).This study uses the default settings of the fminimax() function in Matlab to perform the optimization calculations.The 3 r d part of the pseudocode presented in Fig. 2 clearly shows the flowchart of the New 'Peel off' Strategy.

E. Simulated EMG Signals
To obtain the simulated EMG signals, we assume that the simulated muscle tissue comprised 120 normally distributed MUs with a diameter of 16 mm.Each MU contained 50-1000 fibers whose intracellular action potentials were modeled by the 3-Gaussian model reported by Ning and Zhang [32], [33].All semi-fiber lengths were set to 50 mm, and the endplate and tendon positions of the fibers were uniformly distributed in the range of ±4 mm.To simplify the model, MU conduction velocities were set to 4.0 m/s.The inter-spike interval variability followed a Gaussian distribution with a coefficient of variation of 20%.An 8 × 8 electrode array grid with a 4 mm interelectrode distance in both directions was generated and centered above the simulated MUs.A sampling rate of 2,000 Hz was assumed.The maximum input excitation level was set to 10% and the maximum number of recruited motor units was 70.Gaussian zero-mean noises with signalto-noise (SNR) of 10, 15, 20, 25, and 30 dB were added to the synthetic noise-free sEMG data to simulate the measured, noise-contaminated sEMG data.

F. Experimental EMG Signals
The proposed framework's performance was evaluated using 128-channel HD-sEMG signals collected from the biceps brachii muscle (BBM) of 10 healthy subjects (4 male and 6 female) with ages ranging from 21 to 67 years and a mean age of 34.The experiment was approved by the University of Houston and the University of Texas Health Science Center-Houston Institutional Review Board.Subjects were seated with the target arm secured firmly in two adjustable metal plates.The skin overlaying the BBM was gently abraded and cleaned with alcohol pads.As shown in Fig. 3, two 8 × 8 electrode arrays (TMSi, Enschede, the Netherlands) with a diameter of 4.5 mm and an inter-electrode distance of 8.5 mm were placed over the muscle belly along the muscle fiber.The reference electrode was placed on the lateral epicondyle of the humerus, and a ground electrode was attached to the wrist of the contralateral arm using a fully soaked Velcro wrist band (TMSi, Enschede, the Netherlands).The subjects were instructed to perform maximal voluntary contractions (MVCs) three times and the highest force value observed was used as the 100% MVC force.Then three 10-second isometric contractions were performed at 10%, 30%, and 50% MVC.A 2 min rest interval was given between trials to avoid muscle fatigue.HD-sEMG signals were recorded via a Refa-136 amplifier (TMSi, Enschede, the Netherlands) at a sampling rate of 2,048 Hz with a 24-bit resolution.

G. Evaluation
For the synthetic sEMG data, the decomposition results of the proposed framework were compared with those of KmCKC in terms of the number of decomposed MUs and the F1-score [21].The rationale for using KmCKC for comparison is that the first step of the proposed framework in this study is a simplified version of KmCKC; while the other two parts are the novelty of the proposed framework which work together to make sure the accurate extraction of lowenergy MUs.Therefore, comparing to KmCKC can give us the most intuitive comparison pertaining to the performance of the proposed framework.To further verify the accuracy of the new 'peel-off' strategy, we first compare the similarity between the MUAP waveforms respectively constructed by our approach and Chen's peel-off strategy [22] concerning the simulated MUAP waveform.The similarity between the MUAP waveform of the j th MU calculated by the new 'peeloff' strategy or Chen's peel-off strategy and the simulated MUAP waveform of the j th MU is calculated as follows: where MU A P k s j is the MUAP waveform train of the j th MU in the simulated signal at the k th channel, and MU A P k cj stands for the MUAP waveform train of the j th MU calculated by the new 'peel-off' strategy or Chen's peel-off strategy in the k th channel.Cn represents the number of channels.| • | means absolute value operation.mean(•) stands for mean operation.
We used the Shapiro-wilk test to check the normality of the following data: the number of MUs decomposed by KMCKC with the proposed method in the simulated signal, as well as the corresponding F1-score, at six noise levels in the simulated signal; the number of MUs decomposed by KMCKC and the proposed method in the experimental data across all subjects at three MVC levels; the similarity of the MUAP waveforms obtained by the proposed new 'peel-off' strategy and Chen's peel-off strategy to the simulated MUAP waveforms at six noise levels, where the MU is obtained by decomposing the initial simulated signal by ICKC and Post-Processor.The test results showed that all of the above data followed a normal distribution, except for the F1-score obtained by KMCKC decomposition of the simulated signal at six noise levels and the number of MUs decomposed by KMCKC in the experimental data across all subjects at three MVC levels.Therefore, we conducted statistical analysis using different tests.Specifically, we used a paired samples t-test to examine statistical differences between two groups of normally distributed samples with identical sample sizes.We used a one-sample t-test to examine statistical differences between a group of normally distributed samples and a fixed value.We used the Mann-Whitney U test to examine statistical differences between two groups of samples, with at least one group not conforming to a normal distribution.It should be noted that due to the random selection of clustering centers, the decomposition results of the KmCKC algorithm for the same data have some variability.To ensure a more reliable and stable outcome, the decomposition results of the KmCKC algorithm are the average of 10 decompositions of the same data.

A. Simulated Surface EMG
Fig. 4 shows the two representative decomposition results of the proposed ICKC-P framework, with an SNR of 10 dB.Fig. 4 (a1) shows a representative example of the LIPT, where it is impossible for us to accurately extract the firing instants of this MU due to the small and irregular amplitude of the IPT.On the contrary, with the application of SICKC, we can obtain an IPT as shown in Fig. 4 (a2).The firing instants of the MU can be accurately detected from Fig. 4 (a2) with much stronger and regular peaks in IPT, which highly demonstrates the ability of SICKC to extract low-energy MUs.Fig. 4 (b1) shows the representative IPT after the application of the 'peeloff' strategy and ICKC.Fig. 4 (b2) shows the IPT generated through equation ( 3) using the firing instants of the IPT in Fig. 4 (b1) and the initial synthetic surface EMG.Despite the use of accurate firing instants to estimate the IPT, the resulting IPT remains irregular, making it challenging to extract the correct firing instants for this low-energy MU during the initial decomposition.This is because the high-energy MU information in the initial synthesis signal masks the lowenergy MU information.However, as shown in Fig. 4 (b1), the 'peel-off' strategy and ICKC enable us to obtain this regular IPT after removing the high-energy MU information from the initial synthesis signal.This suggests that the 'peeloff' strategy can provide further guarantee for the extraction of the low-energy MUs.
Fig. 5 (a) shows the comparison between originally synthetic sEMG signals, summation of decomposed MUAPs and the residual of a typical sEMG channel.The reconstructed signal matched very well with the original signal and the residual is small.The average number of decomposed MUs and the F1-scores achieved by the proposed ICKC-P framework and KmCKC across ten simulation trials were presented in Fig. 5 (b1) and (b2), respectively.The results show that as the SNR decreases, the number of MUs and F1 scores obtained by both decomposition techniques decrease, which is in line with common sense.However, at each noise level, our framework yields a higher number of MUs compared to KmCKC.For instance, at an SNR of 30 dB, 11 and 18 MUs were respectively decomposed by KmCKC and our framework; and at an SNR of 10 dB, 7 and 16 MUs were decomposed, respectively.Table I shows that the number of decomposed MUs for the two methods was significantly different at each noise level ( p = 0.000 < 0.01 for all six noise levels) based on one-sample t-tests.However, no significant difference in F1-Score at each noise level was found between the two methods after Mann-Whitney U tests ( p > 0.05 for all six noise levels).Nevertheless, all the pulse trains decomposed by our framework present 93.04% accuracy, even at a low SNR of 10 dB, which demonstrates the robustness of the proposed framework in the presence of noise.Fig. 5 (b3) shows that at all considered noise levels, compare to the MUAP reconstructed by Chen's peel-off strategy, the MUAP reconstructed by our framework presents a relatively high similarity with the simulated MUAP.As shown in Table I, after performing paired samples t-tests, the two methods were significantly distinct from each other in terms of their similarity to the simulated MUAP waveform at each noise level ( p = 0.000 < 0.01 for all six noise levels).The above results highly demonstrate that the proposed framework outperforms the KmCKC and Chen's peel-off strategy.

B. Experimental Surface EMG Signals
Table II displays the decomposition results yielded by KmCKC and our framework on the experimental sEMG data, respectively.Across ten subjects, 13.1 ± 1.6, 17.1 ± 2.6 and 17.7 ± 4.1 MUs were identified by the proposed framework at 10%, 30% and 50% MVC, with COV I S I of 0.23 ± 0.01, 0.22 ± 0.02 and 0.22 ± 0.03; on average 7.9 ± 1.9, 10.0 ± 1.9 and 10.7 ± 2.2 MUs were identified by KmCKC with COV I S I of 0.24 ± 0.02, 0.23 ± 0.03 and 0.23 ± 0.2.Table III indicates a significant difference in the number of decomposed MUs between the two methods at the three levels of contraction after the Mann-Whitney U test ( p = 0.000 < 0.01 for all three levels).Across all trials of all subjects, an average of 16.1 ± 3.4 MUs were extracted from each trial with a COV I S I of 0.22 ± 0.02; an average of 9.6 ± 2.3 Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.Fig. 4. A demonstration of reconstructing MUs using the proposed framework.The synthetic surface EMG was simulated in a condition of 10 dB SNR.(a1) A representative example of LIPT.(a2) The precise IPT estimated in the last iteration in SICKC, which matches (a1).(b1) After the 'peel-off' strategy, the precise IPT estimated in the last iteration in ICKC.(b2) A representative IPT of low-energy MU estimated using all correct firing instants and initial synthetic surface EMG, which matches (b1).
MUs were extracted by KmCKC with a COV I S I of 0.23 ± 0.04.The average energy of the residual signal was 33.1% ± 6.8% of the original signal, which was acquired by subtracting the summation of the MUAP decomposed by the ICKC-P framework from the original surface signal.Fig. 6 shows the decomposition results of a representative subject (♯3) with a contraction level of 30% MVC.Fig. 6 (a) shows the measured sEMG signal, the summation of decomposed MUAPs, and the residual of a typical channel.As in the simulation, the summation of the decomposed MUAPs matched the measured signal very well, indicated by a very small residual.Fig. 6 (b) shows that for this specific trial, 16 MUs were extracted by using the proposed ICKC-P framework, while only 12 MUs were acquired by KmCKC.

IV. DISCUSSION
The firing pattern and MUAP morphology offered by the EMG decomposition technique can not only enhance our understanding pertaining to the neurophysiology of the neuromuscular system but also provide vital information for the diagnosis and progress tracking of motor neural diseases and neuromuscular disorders [34], [35], [36], [37].Compared to the intramuscular EMG, the low-pass filter effect caused by the skin and transcutaneous fat separating the surface electrodes from muscle post significant technical challenges to the sEMG signal decomposition.
Currently, blind source separation methods are commonly adopted, and their performance has been widely validated under normal and neuropathological conditions [5], [6], [7], [8], [9], [10], [11], [12], [13], [14], [15], [16], [17], [18], [19], [20], [21], [22], [23], [24], [25], [26], [27], [28], [29], [30], [31], [32], [33].Despite the efforts made, these types of approaches, including the classical CKC, GCKC, and KmCKC present a bias toward the high-energy MUs [18], [21], [39] and no special attention has been given to the low-energy MUs which are critically important in motor unit number estimation and botulinum neural toxin (BoNT) injection in the clinical treatment of spasticity.So far, only the 'peel-off' strategy in the PFP framework [22] and the source deflation procedure proposed by Negro [30] can decompose part of the low-energy MUs.However, neither of the aforementioned two approaches targeted the extraction of the low-energy MUs.In this study, a novel ICKC-P framework targeting both high and low-energy MUs was developed by combing the iterative CKC, Post-Processor and the new 'peel-off' strategy.The Post-Processor was employed to extract the low-energy MUs that are not totally swamped by high-energy MUs and noise by maintaining the sparsity of low-energy MU firing instants.This is a significant improvement over the traditional sEMG decomposition framework.Furtherly, the new 'peel-off' strategy, was dedicated to decomposing the low-energy MUs that are masked by the high-energy MUs by optimizing the MUAP waveform based on the orthogonality and sparsity of the source signal.
To evaluate the performance of our proposed ICKC-P Framework, both computer-simulated and experimental sEMG   11 and 18 IPTs were respectively decomposed by KmCKC and the ICKC-P framework; while with an SNR of 10 dB, the ICKC-P framework totally outperformed the KmCKC with 9 more decomposed MUs.Therefore, the proposed framework presents significant superiority especially in a noisy environment, due to the post-processor and the new 'peel-off' strategy.The similarity between the MUAP waveform calculated by the proposed new 'peel-off' strategy and the simulated MUAP Although the proposed ICKC-P framework boasts superior performance compared to typical methods it has some limitations: 1) Model defects.As with all other blind source separation type of decomposition methods (including the classical CKC [18], GCKC [38], KmCKC [21], the PFP framework [22] and the framework proposed by Negro et al. [30]), the ICKC-P framework is based on the assumption that the sources are independent to each other, i.e., only very limited number of MUs discharge simultaneously at a specific time instant.However, under some circumstances, this condition does not hold, such as when the sampling frequency is too low or too many MUs are activated at high contraction level.Though the great effort was made in the Post-Processor to address the superimposition issue, this problem cannot be fully eliminated, but just alleviated.2) MUAP shape alteration.In the 'peel off' strategy proposed by Chen et al., they assumed that MUAPs from the same motor unit remain stationary over a short period of time (about 10s).Nonetheless different external environments, including temperature [39] or muscle fatigue [40], may cause the deformation of the MUAP morphology, even during constant force isometric contraction.In this paper, though the MUAP morphology alteration was taken into consideration, it is still an approximate estimation.An adaptive MUAP construction algorithm, together with strategies to reduce interference from the external environment remain our future efforts.3) Time Consuming.In the ICKC-P framework, many parameters need to be set and adjusted manually and empirically, such as the h1 in ICKC and the ρ in Post-Processor.Therefore, this algorithm is time consuming and highly depends on the researcher's experience.4) Performance validation.In this study, the performance of the proposed decomposition framework was compared with those acquired by KmCKC, in terms of the number of identified MUs, C O V I S I and the ratio between the residual signal and the source signal.For experimental surface EMG data, '2-source' validation is the most commonly accepted method to quantify the decomposition performance [41].However, the '2-source' validation in the current study was not carried out due to the lack of concurrent intramuscular EMG recordings.
V. CONCLUSION In this study, by combining ICKC, Post-Processor and a new 'peel-off' strategy, a novel framework for HD-sEMG decomposition was developed.The ICKC method was adopted to extract the IPT of large-energy MUs, while the Post-Processor and 'peel-off' strategy were specifically developed to identify the small MUs with low-energy.Both simulation and experimental results highly suggest the proposed ICKC-P framework outperforms the KmCKC approach in terms the number of identified MUs while remaining comparable decomposition accuracy.

Fig. 1 .
Fig. 1.The flowchart of the proposed sEMG decomposition framework.The proposed framework consists of three parts: ICKC, Post-Processor (including a strategy for screening IPT with different energies and a strict ICKC method), and Peel-off.

Fig. 2 .
Fig. 2. The pseudo-code of the decomposition framework.The markers 1, 2 and 3 stand for ICKC, Post-Processor and Peel-Off Strategy, respectively.

Fig. 5 .
Fig. 5. (a)Top panel represents the original simulated surface EMG signal (black lines) in one typical channel with SNR = 10 dB, and the mid panel represents the summation of the MUAP trains (red lines) decomposed by the proposed Framework, and the bottom panel represents the residual acquired by subtracting the summation of the decomposed MUAPs from the original surface EMG signal.(b1) shows the number of reconstructed IPTs and (b2) shows the F1-Score of the reconstructed IPTs achieved by the ICKC-P framework and KmCKC at different noise levels.(b3) shows the similarity of the MUAP waveforms obtained with the New 'peel-off' strategy and Chen's peel-off strategy to the simulated MUAP waveforms at different noise levels.The results of KmCKC have averaged over ten simulation trials.N in the horizontal axis denotes conditions for which no noise was added to the simulated signals.

Fig. 6 .
Fig. 6.(a) The top panel shows the sEMG (black lines) recordings from BBM (Subject ♯3) in one typical channel.The mid panel shows summation of the MUAPs (red lines) extracted by the proposed framework, and the bottom panel shows the residual (blue lines) acquired by subtracting the summation of the decomposed MUAPs from the original surface EMG signal.(b) The right panel shows MU firing patterns identified by the KmCKC (red lines) and our framework (blue lines) with a contraction level of 30% MVC.Erroneous spikes and missed spikes are indicated by "•".The left panel shows the constructed MUAP waveform of each MU.

TABLE I STATISTICAL
DIFFERENCES (P-VALUE) BETWEEN THE FOLLOWING SAMPLE GROUPS: NUMBER OF MUS DECOMPOSED BY KMCKC AND THE PROPOSED METHOD, AND THE CORRESPONDING F1 SCORES, AT SIX NOISE LEVELS OF THE SIMULATED SIGNAL, AS WELL AS SIMILARITY OF THE MUAP WAVEFORMS OBTAINED BY THE PROPOSED NEW 'PEEL-OFF' STRATEGY AND CHEN'S PEEL-OFF STRATEGY TO THE SIMULATED MUAP WAVEFORMS signals were adopted in this study.For the simulated sEMG, the number of MUs decomposed by KmCKC and the ICKC-P Framework gradually decreased with the decline of SNR.at the same SNR level, the ICKC-P Framework can always decompose a larger number of MUs than the KmCKC method.Concretely, with an SNR of 30 dB, Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.

TABLE II DECOMPOSING
RESULTS OF THE BBM SIGNALS AT 10%, 30% AND 50% MVC USING KMCKC AND ICKC-P FRAMEWORK TABLE III STATISTICAL DIFFERENCES (P-VALUE) BETWEEN THE NUMBER OF MUS DECOMPOSED BY KMCKC AND THE PROPOSED METHOD IN EXPERIMENTAL DATA ACROSS ALL SUBJECTS AT THREE MVC LEVELS waveform almost keeps constant as the SNR increases.At the same SNR level, our newly proposed 'peel-off' strategy can estimate the MUAP waveforms with higher accuracy indicating higher similarity with the simulated MUAP.By directly subtracting the optimized C j ( j ) from the signal corresponding to the moment in j (Fig. 2), the proposed 'peel-off' strategy could avoid the alignment problem.For the experimental data, by taking advantage of the Post-Processor and the 'peel-off' procedure, low-energy MUs that cannot be decomposed by KmCKC were successfully decomposed by our ICKC-P framework with an average yield of 16.1 ± 3.4 MUs, COV I S I of 0.22 ± 0.02 and residual of 33.1% ± 6.8%; while for the KmCKC, only 9.6 ± 2.3 MUs were extracted with an average COV I S I of 0.23 ± 0.04.