Identification and Analysis of Limb Rehabilitation Signal Based on Wavelet Transform

Figure 1 explains the basic flow of the formation of sEMG signal in a motor unit during limb rehabilitation exercise. The motor unit, stimulated by the neurotransmitter in synaptic cleft, produces an action potential on the relevant muscle fiber. The action potential will be transmitted electrochemically. The transmission will cause muscle fiber to change and the muscle to contract. The repeated muscle contractions continuously induce and activate the action potential of the motor unit. In this way, a sequence of action potentials will appear. These potentials will superimpose each other in time and space, thereby producing EMG signal.

1.png

Figure 1. Formation of sEMG signal

If the limbs of the patient are gradually improving, the theory of muscle synergy indicates that the central nervous system will send a limb control command. After receiving the command, the corresponding limb muscles will be regulated well. These muscles will be linearly combined in order, and a muscle contraction mode will form. In this way, the limb could complete the right motion function. Mathematically, the i-th correct limb motion function u_i equals the linear weighting of the muscle synergy g_j and activation coefficient q_ij of M groups of limb muscles:

$\begin{align}  & {{u}_{i}}\left( l \right)\approx \sum\nolimits_{j=1}^{M}{{{q}_{ij}}{{g}_{j}}\left( l \right)}\text{   } \\ & \left\{ i=1,...,N,j=1,...,M,M<N \right\} \\\end{align}$                  (1)

Let O be the number of synergic elements; N be the dimensionality of the eigenvector; M be the number of samples; u be the N×M matrix of muscle motion samples of limb rehabilitation; Q be the muscle synergy matrix with N×O motion patterns; G be the O×M activation coefficient matrix of muscle synergy; F be the N×M matrix of noises. Then, the muscle motions of limb rehabilitation can be modeled as:

$U=Q\times G+F$               (2)

Formula 2 is a nonnegative matrix. There are two classic multiplicative iteration methods to decompose such a matrix: Kullback-Leibler (KL) multiplicative iteration, and Euclidean distance iteration. The KL multiplicative iteration can be expressed as:

$\underset{Q\text{ }G}{\mathop{min}}\,T\left( U\left\| QG \right. \right)=\sum{_{ij}}\left( {{U}_{ij}}log\frac{{{u}_{ij}}}{{{\left( QG \right)}_{ij}}} \right)-{{U}_{ij}}+{{\left( QG \right)}_{ij}}$ st. Q∈R^N×O, H∈R^O×M Q, G≥0        (3)

The iterative formula to solve formula 3 can be derived as:

${{Q}_{ij}}\leftarrow {{Q}_{ij}}\frac{\sum{_{ij}}{{G}_{jl}}{{U}_{il}}/{{\left( QG \right)}_{il}}}{\sum{_{O}}{{G}_{jO}}};$

${{G}_{jl}}\leftarrow {{G}_{jl}}\frac{\sum{_{i}}{{Q}_{ij}}{{U}_{il}}/{{\left( QG \right)}_{il}}}{\sum{_{\tau }}{{Q}_{\tau j}}};$                (4)

In Euclidean distance iteration, the muscle motions of limb rehabilitation can be modeled as:

$F=U-QG\text{   }s.t.Q,G\ge 0$                     (5)

Formula 5 shows that ||F||² must be minimized to solve matrices Q and G, while preserving the effective information in the original matrix U. Provided that the noise matrix obeys Gaussian distribution, the maximum likelihood formula can be derived from the training samples of the muscle motions of limb rehabilitation:

$\begin{align}  & K\left( Q,G \right)=\prod{_{i,j}}\frac{1}{\sqrt{2\pi }{{\phi }_{ij}}}exp\left( -\frac{F_{ij}^{2}}{2{{\phi }_{ij}}} \right) \\ & =\prod{_{i,j}}\frac{1}{\sqrt{2\pi {{\phi }_{ij}}}}exp\left( -\frac{{{\left( {{U}_{ij}}-{{\left( QG \right)}_{ij}} \right)}^{2}}}{2{{\phi }_{ij}}} \right) \\\end{align}$               (6)

Taking the logarithm of formula 6, the logarithmic likelihood function can be established as:

$K\left( Q,G \right)=\sum{_{ij}ln\frac{1}{\sqrt{2\pi }{{\phi }_{ij}}}-\frac{1}{2{{\phi }_{ij}}}\sum{_{i,j}}}{{\left[ {{U}_{ij}}-{{\left( QG \right)}_{ij}} \right]}^{2}}$          (7)

Ignoring the variance difference between noises, formula (7) is equivalent to:

$ZS\left( Q,G \right)=\frac{1}{2}\sum{_{ij}}{{\left[ {{U}_{ij}}-{{\left( QG \right)}_{ij}} \right]}^{2}}$                 (8)

where, (QG)_ij can be calculated by:

${{\left( QG \right)}_{ij}}=\sum{_{l}{{Q}_{il}}}{{G}_{lj}}\Rightarrow \frac{\partial {{\left( QG \right)}_{ij}}}{\partial {{Q}_{il}}}={{G}_{kj}}$                     (9)

Combining formulas 8 and 9:

$\begin{align}  & \frac{{{\partial }_{J}}{{\left( QG \right)}_{ij}}}{\partial {{Q}_{il}}}=\sum{_{j}}\left[ -{{G}_{hj}}\left( {{U}_{ij}}-{{\left( QG \right)}_{ij}} \right) \right] \\ & ={{\left( QG{{G}^{\psi }} \right)}_{il}}-{{\left( U{{G}^{\psi }} \right)}_{il}} \\\end{align}$                   (10)

Similarly, we have:

$\frac{{{\partial }_{ZS}}{{\left( QG \right)}_{ij}}}{\partial {{Q}_{il}}}={{\left( {{Q}^{\psi }}QG \right)}_{lj}}-{{\left( {{Q}^{\psi }}U \right)}_{lj}}$                  (11)

Through gradient descent, Q_il can be iteratively calculated by:

${{Q}_{il}}={{Q}_{il}}-{{\beta }_{1}}\cdot =\left[ {{\left( QG{{G}^{\psi }} \right)}_{il}}-{{\left( U{{G}^{\psi }} \right)}_{il}} \right]$                  (12)

Similarly, G_ij can be calculated by:

${{G}_{lj}}={{G}_{lj}}-{{\beta }_{2}}\cdot \left[ {{\left( {{Q}^{\psi }}QG \right)}_{lj}}-{{\left( {{Q}^{T}}U \right)}_{lj}} \right]$        (13)

β₁ and β₂ can be calculated by:

${{\beta }_{1}}=\frac{{{Q}_{il}}}{{{\left( QG{{G}^{\psi }} \right)}_{il}}}\text{   }{{\beta }_{2}}=\frac{{{G}_{lj}}}{{{\left( {{Q}^{\psi }}QG \right)}_{lj}}}$                (14)

Q_il and G_ij can be iteratively calculated by:

${{Q}_{il}}\leftarrow {{Q}_{il}}\cdot \frac{{{\left( U{{G}^{\psi }} \right)}_{il}}}{{{\left( QG{{G}^{\psi }} \right)}_{il}}}\text{   }{{G}_{lj}}\leftarrow {{G}_{lj}}\cdot \frac{{{\left( {{Q}^{\psi }}U \right)}_{lj}}}{{{\left( {{Q}^{\psi }}QG \right)}_{lj}}}$                      (15)

Through the above derivation, it is possible to clarify the muscle synergy process in limb rehabilitation by decomposing the nonnegative matrix (Figure 2).

2.png

Figure 2. Muscle synergy extraction through nonnegative matrix decomposition

Through nonnegative matrix decomposition, matrix U was decomposed into matrices Q and G. However, the number of columns in Q and the number of synergic elements O are both uncertain. This is a key factor in solving muscle synergy during limb rehabilitation analysis. Once O is determined, it is easy to determine the Q and G decomposed from U.

The common way to determine the number of synergic elements O is to calculate the data similarity UDA between the reconstructed matrix U'=Q×G of muscle motions in limb rehabilitation and the original matrix U:

$UDA=1-\frac{{{\left( u-u' \right)}^{2}}}{{{u}^{2}}}$                  (16)

The closer the U' is to U, the greater the UDA. The inverse is also true.

4.1 Denoising principle

During the EMG signal collection from patient limbs during rehabilitation, the EMG signal can be divided into upper limb EMG signal, lower limb EMG signal, and hand EMG signal, according to the placement of electrodes on the patient body. In actual clinical detection, the noninvasive sEMG signal collection technique is often adopted. However, the collected signal is usually too weak and prone to noise disturbance. The noisy signal will cause errors to the classification of mode data. Therefore, it is important to denoise the EMG signal before pattern recognition and evaluation of rehabilitation status. To ensure the denoising effect, this paper combines wavelet thresholding with TVD (Figure 4).

4.png

Figure 4. Flow of proposed combinatory denoising strategy

The total variation of the one-dimensional (1D) EMG signal can be defined as:

$QB\left( a \right)=\sum\limits_{M=1}^{m}{\left| a\left( M+1 \right)-a\left( M \right) \right|}$                  (22)

Formula 22 can be rewritten as QB(a)=||ξa|| ₁, which is the l₁-norm of ξa. Ξ can be expressed as:

$\xi ={{\left( \begin{matrix}   -1 & 1 & {} & {} & {} & {}  \\   {} & -1 & 1 & {} & {} & {}  \\   {} & {} & {} & \ddots  & {} & {}  \\   {} & {} & {} & {} & -1 & 1  \\\end{matrix} \right)}_{m\times \left( m+1 \right)}}$                      (23)

Suppose b= b(M) is a signal with additive noise. Then:

$b=a+Q$                 (24)

where, a =a(M)∈ℝ^m is the clean signal; Q = Q(M)∈ℝ^m is noising signal.

The proposed combinatory denoising strategy is as follows: First, wavelet thresholding is carried out in three steps:

Step 1. Carry out wavelet decomposition of EMG signal b(M) to obtain detail coefficients and scale coefficients, and set thresholds for the obtained scale coefficients.

Step 2. Estimate the standard deviation ε of the noise based on the detail coefficients. Suppose the EMG signal length is m_s. Then, define the optional wavelet threshold as ε·(2log_em_s)^0.5. Choose the threshold just greater than the maximum noise amplitude to process the detail coefficients on each scale.

Step 3. Carry out inverse wavelet transform based on ε, and further reconstruct the EMG signal. That is, perform TVD regularization on the denoised signal a^*(M).

The TVD consists of three steps:

Step 1. Compute the allowance AL(M) of b(M) and a^*(M):

$AL\left( M \right)=b\left( M \right)-{{a}^{*}}\left( M \right)$                  (25)

Step 2. Let β be the regularization parameter to adjust weights; ξ be the first-order differential matrix. Denoise AL(M) through TVD regularization:

$\begin{align}  & A{{L}^{*}}\left( M \right)=\arg \underset{AL}{\mathop{\min }}\,\frac{1}{2}\left\| AL\left( M \right)-A{{L}^{*}}\left( M \right) \right\|_{2}^{2} \\ & +\beta {{\left\| \xi A{{L}^{*}}\left( M \right) \right\|}_{1}} \\\end{align}$                  (26)

Formula 26 shows that the denoised allowance contains a fidelity term and a TVD regularization term. If the latter term does not play a punitive role, β=0; if the latter term dominates the process, β approaches ∞. In the latter case, the TVD regularization term must be very small to reconstruct the EMG signal. However, a small TVD regularization term will dampen the denoising effect, and reduce the fidelity. Therefore, a suitable β value should be selected to optimize the denoising effect.

Step 3. Let a(M) be the noise-free clean EMG signal. Reconstruct the original signal from the denoised signal of wavelet thresholding and that of TVD regularization:

$a\left( M \right)={{a}^{*}}\left( M \right)+A{{L}^{*}}\left( M \right)$                     (27)

It is difficult to find the derivative of the TVD term ||ξAL^||₁ in the absolute form. To overcome the difficulty, this paper iteratively solves the derivative alternating direction method of multipliers (ADMM). According to the TVD regularization denoising model, formula 26 can be converted into a constraint minimization problem:

$A{{L}^{*}}\left( M \right)=arg\underset{A{{L}^{*}}}{\mathop{min}}\,\frac{1}{2}\left\| AL\left( M \right)-A{{L}^{*}}\left( M \right) \right\|_{2}^{2}+\beta {{\left\| v \right\|}_{1}}$                       (28)

The objective function of the optimization problem can be described by:

$\text{ }\xi A{{L}^{*}}\left( M \right)-v=0$                    (29)

Let μ=(μ₁; μ₂; …; μ_m^s_-1) be the augmented Lagrange multiplier; σ be the penalty parameter. Then, formula 29 can be converted into an unconstrained optimization problem:

$\begin{align}  & K\left( A{{L}^{*}}\left( M \right),v,\mu  \right)=arg\underset{v,\overset{\wedge }{\mathop{AL}}\,}{\mathop{min}}\,\frac{1}{2}\left\| AL\left( M \right)-A{{L}^{*}}\left( M \right) \right\|_{2}^{2} \\ & +\beta {{\left\| v \right\|}_{1}}+{{\mu }^{\psi }}\left( \xi A{{L}^{*}}\left( M \right)-v \right)+\frac{\sigma }{2}\left\| \xi A{{L}^{*}}\left( M \right)-v \right\|_{2}^{2} \\\end{align}$                 (30)

Because it is impossible to directly find the derivative of absolute value, make v=v_l and find the derivative of AL* in formula 30:

$\begin{align}  & \frac{\partial K}{\partial A{{L}^{*}}\left( M \right)}=-\left( AL\left( M \right)-A{{L}^{*}}\left( M \right) \right)+{{\xi }^{\psi }}{{\mu }_{l}} \\ & -\sigma {{\xi }^{\psi }}\left( \xi A{{L}^{*}}\left( M \right)-{{v}_{l}} \right)=0 \\\end{align}$                     (31)

Solving formula 31:

$A{{L}^{*}}_{l+1}\left( M \right)=\frac{AL\left( M \right)-{{\xi }^{\psi }}{{\mu }_{l}}+\sigma {{\xi }^{\psi }}{{v}_{l}}}{HS+\sigma {{\xi }^{\psi }}\xi }$             (32)

Let Φ(∙) be the 1D contraction operator. Combining formula 32 with formula 30 to solve v_l+1:

$\begin{align}  & {{v}_{l+1}}=\arg \underset{v}{\mathop{\min }}\,\beta {{\left\| v \right\|}_{1}}+\mu _{^{l}}^{\psi }\left( \xi A{{L}^{*}}_{l+1}\left( M \right)-v \right) \\ & +\frac{\sigma }{2}\left\| \xi A{{L}^{*}}_{l}\left( M \right)-v \right\|_{2}^{2}=\arg \underset{v}{\mathop{\min }}\,\beta {{\left\| v \right\|}_{1}} \\ & +\frac{\sigma }{2}\left\| \xi A{{L}^{*}}_{l+1}\left( M \right)-v+\frac{{{\mu }_{l}}}{\sigma } \right\|_{2}^{2} \\ & =\Phi \left( \xi A{{L}^{*}}_{l+1}\left( M \right)+\frac{{{\mu }_{l}}}{\sigma },+\frac{\beta }{\sigma } \right) \\\end{align}$                      (33)

4.2 Evaluation of denoising effect

This paper evaluates the denoising effect of EMG signal with root mean square error E_RMSE, signal-to-noise ratio E_SNR, and Pearson correlation coefficient E_P. Let CX(M) be the clean EMG signal before adding any noise; a(M) be the denoised EMG signal. Then, E_RMSE can be calculated by:

${{E}_{RMSE}}=\sqrt{\frac{1}{m}\sum\limits_{M=1}^{m}{{{\left( a\left( M \right)-CX\left( M \right) \right)}^{2}}}}$                (34)

The smaller the E_RMSE, the better the denoising effect. E_SNR can be calculated by:

${{E}_{SNR}}=10{{\log }_{10}}\left( \sum\limits_{M=1}^{m}{C{{X}^{2}}\left( M \right)}/\sum\limits_{M=1}^{m}{{{\left( a\left( M \right)-CX\left( M \right) \right)}^{2}}} \right)$                 (35)

The greater the E_SNR, the better the denoising effect. The Pearson correlation coefficient E_P measures the correlation between two signals. The value range of the coefficient is [-1, 1]. If E_P=1, then the two signals are positively correlated; If E_P=-1, then the two signals are negatively correlated. Let CX'(M) and a'(M) be the mean of CX(M) and a(M), respectively. Then, E_P can be calculated by:

$\gamma =\frac{\sum\nolimits_{M=1}^{m}{\left( CX\left( M \right)-C{X}'\left( M \right) \right)\left( a\left( M \right)-{a}'\left( M \right) \right)}}{\sqrt{\sum\nolimits_{M=1}^{m}{{{\left( CX\left( M \right)-C{X}'\left( M \right) \right)}^{2}}}}\sqrt{\sum\nolimits_{M=1}^{m}{{{\left( a\left( M \right)-{a}'\left( M \right) \right)}^{2}}}}}$                     (36)

The greater the E_P, the better the denoising effect.

This paper probes deep into the identification and analysis of limb rehabilitation signal based on wavelet transform. Firstly, the authors detailed the basic flow of sEMG signal generation in motor unit during limb rehabilitation exercise, and presented a limb EMG pattern recognition method. Next, SVM was selected to recognize the pattern of the EMG signal extracted from the limb rehabilitation exercise of patients, and to judge the rehabilitation status. Through experiments, the authors plotted the processed muscle synergy matrices and the activation coefficient matrices, disclosed the distribution of the features of the same motion in limb rehabilitation in the feature space, and judged the EMG pattern and the stage of rehabilitation. Finally, the wavelet thresholding was integrated with TVD to improve the noise removal effect of EMG signals. The EMG signal denoising effects of multiple methods were compared through experiments. The comparison shows that our method achieved the best denoising effect, the highest robustness, and the strongest interference resistance.