# sEMG-MMG State-Space Model for the Continuous Estimation of Multijoint Angle

Authors: Xugang Xi; Chen Yang; Seyed M. Miran; Yun-Bo Zhao; Shuliang Lin; Zhizeng Luo
Journal: Complexity
 (2020)

---

## Abstract

Continuous joint angle estimation plays an important role in motion intention recognition and rehabilitation training. In this study, a surface electromyography- (sEMG-) mechanomyography (MMG) state-space model is proposed to estimate continuous multijoint movements from sEMG and MMG signals accurately. The model combines forward dynamics with a Hill-based muscle model that estimates joint torque only in a nonfeedback form, making the extended model capable of predicting the multijoint motion directly. The sEMG and MMG features, including the Wilson amplitude and permutation entropy, are then extracted to construct a measurement equation to reduce system error and external disturbances. Using the proposed model, a closed-loop prediction-correction approach, unscented particle filtering, is used to estimate the joint angle from sEMG and MMG signals. Comprehensive experiments are conducted on the human elbow and shoulder joint, and remarkable improvements are demonstrated compared with conventional methods.


---

## Body

## 1. Introduction

In human-machine interaction (HMI), surface electromyography (sEMG) is often used to serve as the input signal source. The sEMG signal is a weak electrical potential generated by muscle cells upon electrical or neurological activation, and it is detected from superficial muscles by using surface electrodes. The signal contains abundant information and has a distinct characteristic. Moreover, its collection is simple and not traumatic. The sEMG signal has become a research hotspot in the field of human-computer interaction technology, especially in the manufacturing of exoskeleton robots [1], intelligent prosthetics [2], and rehabilitation robots [3].sEMG studies usually focus on feature extraction and pattern classification. Identifying a plurality of discrete action categories has an accuracy rate of approximately 90% [4]. Studies on the classification of discrete actions have been done from the laboratory to the market [5]. Estimating continuous human joint movement is a critical focus in sEMG at present. Many achievements have already been made, especially in the field of rehabilitation robots, where forecasting the continuous motion variables of patients is vital to achieve smooth control of the rehabilitation robot [6]. However, the sEMG signal is susceptible to interference from sweat and skin impedance, which consequently deteriorates the control accuracy.Mechanomyography (MMG) is a mechanical signal produced by contracting muscles, and it reflects muscle activity in the form of low-frequency (<50 Hz) vibration. In contrast to sEMG signals, MMG signals are unaffected by sweat and changes in skin impedance. Recent studies have proven the efficacy of MMG, such as in upper-limb prostheses control [7–9] and hand movement analysis [10–13]. However, the MMG signal has considerably lower signal-to-noise ratios than the sEMG signal. External noise caused by movement artifact can also introduce interference to MMG measurement. Thus, combining the advantages of both sEMG and MMG would be of great significance to understanding muscle activities and enhancing HMI.Extracting features such as root mean square and the absolute value of amplitude and then establishing a motion model with neural networks constitute an effective method for estimating continuous joint motion [14, 15]. The estimation can also be performed through the physiological muscle model, for which the Hill-based muscle model (HMM) is often used [16, 17]. However, the HMM presents two challenges: (1) the HMM involves many complex parameters that are difficult to identify and is computationally heavy, and (2) in the HMM, motion states are indirectly calculated from the sEMG and MMG signals, thereby generating accumulative errors and worsening the estimation accuracy.To address the above challenges, this study developed a new model that fuses sEMG and MMG signals and brings together HMM and joint dynamics. This model can directly calculate the joint motion from sEMG and MMG signals. Wilson amplitude (WA) and permutation entropy (PE) are extracted to establish the measurement equations. Then, an sEMG–MMG state-space model is created to estimate the mutijoint angle.


## 2. Methodology

As shown in Figure1, the sEMG-MMG state-space model consists of three parts: (1) the sEMG and MMG signals to the joint forward dynamics, (2) the fusion of sEMG and MMG features, and (3) the prediction algorithm. The model aims to estimate the joint moving states through the muscle model and fusion features.Figure 1
Structure of the sEMG-MMG state-space model. (1) The sEMG and MMG signals to the joint forward dynamics, (2) the fusion of sEMG and MMG features, and (3) the prediction algorithm.


### 2.1. Data Acquisition

The muscles relevant to elbow and shoulder joint motion were selected, including the biceps brachii (BB), triceps brachii (TB), brachioradialis (BR), trapezius (TZ), teres minor (TM), anterior deltoid (AD), lateral deltoid (LD), and pectoralis major (PM). The sEMG and MMG signals were measured by placing eight surface electrodes on the muscles of the right upper limb (Figure2).Figure 2
sEMG sensor arrangement (Ch1: BB, Ch2: TB, Ch3: BR, Ch4: TZ, Ch5: TM, Ch6: AD, Ch7: LD, and Ch8: PM).The Trigno™ Wireless (Delsys Inc, Natick, MA, USA) was used to record the sEMG and MMG signals at 1600 Hz. The bandwidth was 20–450 Hz, and the baseline noise was less than 1.25μV. The sensors were positioned following the SENIAM recommendation. The reference electrode was attached to the lateral elbow epicondyle.The Codamotion system (Charnwood Dynamics Ltd., UK) was used to record the real joint angle information. This system can capture and analyze motion information in real time through infrared capture. Figure3 shows the Codamotion system and placement of markers.Figure 3
Codamotion system and placement of markers.Ten healthy subjects (five males and five females, 23 ± 2 years old, weighing 50–70 kg) were tested in the experiments. The study was approved by the Hangzhou Dianzi University’s Institutional Review Board, and all subjects provided a written informed consent prior to participation. Figure4 shows the experimental movements. The subject (1) sits on a chair and keeps the muscles relaxed, (2) straightens the right arm forward with elbow extension, (3) performs elbow flexion, (4) straightens the arm upward with elbow extension, (5) straightens the arm to the right side with elbow extension, (6) performs elbow flexion, and (7) returns the upper arm to the initial position while maintaining elbow flexion. The experimental movements involved the main activities of the upper limb, thereby enabling us to better analyze the movement of the shoulder and elbow joints. Each subject completed all movements with steps 1-2-3-4-5-6-7-1 in 15 s and performed the sequence with or without a load (3 kg). The subjects were given an appropriate rest time after each round to reduce muscle fatigue. We completed the data collection within five days, and the subjects repeated the procedure 10 times in each experiment day. In this way, we were able to collect the required data.Figure 4
Upper limb movements recorded: steps 1-2-3-4-5-6-7-1. Schematic of the (a) right-hand side and (b) back side.


### 2.2. Data Preprocessing

Both sEMG and MMG signals are very weak, nonstationary, and random. Many kinds of noise exist in the acquisition process. The sEMG and MMG signals must be denoised for further application. In this experiment, the Trigno Wireless (Delsys Inc.) with a bandpass filter (sEMG and MMG signals of 20–450 and 5–150 Hz, respectively) was used to eliminate most of the inherent noises and motion artifacts.The sampling frequency of the joint angle signal is 100 Hz, and that of the sEMG and MMG signals is 1600 Hz. To keep the sEMG and MMG signals consistent with the joint angle signal, subsampling of the sEMG signals was implemented using feature extraction. Afterward, the subsampling frequency of the sEMG and MMG signals was changed from 1600 Hz to 100 Hz.


### 2.3. sEMG-MMG State-Space Model

#### 2.3.1. Hill-Based Muscle Model

Figure5 illustrates a simplified HMM. The HMM consists of a contractile element and a parallel elastic element. The contractile element produces the active muscle force FAm, and the parallel elastic element produces the passive force FPm. The model can be described as follows:(1)Fm=FPm+FAm=fPl⋅F0m+fAl⋅fVv⋅ak⋅F0m,l=lml0m,v=vmv0m,where Fmt is the musculotendon force, fPl is the normalized passive force-length relationship, F0m is the maximum isometric force, l is the normalized muscle fiber length, lm is the current muscle fiber length, l0m is the optimal fiber length, fAl is the normalized active force-length relationship, fVv is the normalized force-velocity relationship, v is the normalized muscle fiber velocity, vm is the muscle contraction velocity, v0m is the maximum muscle contraction velocity, and ak is the muscle activation [18].Figure 5
Hill-based muscle model.According to the abovementioned equations, the musculotendon force can be described as follows:(2)Fmt=fAl⋅fVv⋅ak+fPl⋅F0m⋅cosϕ,where Fmt is the musculotendon force and ϕ is the pennation angle.The simplifications offAl, fVv, fPl, F0m, lt, l0m, and ϕ can replace the complex biomechanical parameters [19–22]. The musculotendon length can also be simplified as a first-order polynomial [23, 24]. Thus,(3)lmt=b0+b1⋅θ,where lmt is the musculotendon length, b0 and b1 are constants, and θ is the joint angle. Then, the moment arm r and the joint moment τ can be calculated as follows:(4)r=∂lmtθ∂θ=b1,τ=Fmt⋅r.


#### 2.3.2. sEMG-MMG Joint Motion Model

The angular acceleration of the joint can be calculated as follows:(5)θ¨=1Ie⋅τ−τeg,τeg=τegm⋅sinθ,where Ie is the moment of inertia, τeg is the external and gravity torque, and τegm is the maximum of τeg.By combining the above equations, the angular accelerationθ¨ can be described as follows:(6)θ¨=s0+s1⋅θ+s2⋅θ2⋅ak+s3⋅es4⋅θ−s5⋅sinθ,where si (i = 0, 1, …, 5) are constants computed as follows:(7)s0=k0⋅F0m⋅b1⋅cosϕIe+k1⋅F0m⋅b0−lt⋅b1l0m⋅Ie+k2⋅F0m⋅b0−lt2⋅b1l0m2⋅Ie⋅cosϕ,s1=k1⋅F0m⋅b12l0m⋅Ie+2⋅k2⋅b0−lt⋅b12l0m2⋅Ie⋅cosϕ,s2=k2⋅F0m⋅b13l0m2⋅Ie⋅cosϕ,s3=F0m⋅b1⋅cosϕIe⋅exp10⋅b0−lt⋅b12l0m⋅cosϕ−15,s4=10⋅b1l0m⋅cosϕ,s5=τegmIe.Then, we can obtain the sEMG-MMG joint motion model in discrete time as follows:(8)θ¨k+1=s0+s1⋅θk+s2⋅θk2⋅ak+s3⋅es4⋅θk−s5⋅sinθk,θ˙k+1=θ˙k+θ¨k⋅Ts,θk+1=θk+θ˙k⋅Ts,where θ˙k is the joint angular velocity, θk is the joint position, and Ts is the sampling time.


#### 2.3.3. sEMG-MMG Features and State-Space Model

The simplified sEMG–MMG joint motion model involves two problems. First, the model omits certain physiological parameters. Second, the biomechanical parameters vary because of different body conditions. These problems cause accumulative errors in the recursive calculation of the joint angle. To eliminate these errors, we constructed a measurement equation providing the joint angle measurements as feedback, which then leads to a sEMG-MMG state-space model. We built the measure equation by Wilson amplitude (WA) and permutation entropy (PE) in this study.(1)
WA:The WA describes the different characteristics among adjacent amplitudes [25]. It can be described as follows:(9)WA=∑i=1Nuxi+1−xi−T,wherexi and xi+1 are adjacent samples, and the threshold T is 0.05 V in this study.(2)
PE:The PE describes the complexity of the data and has the advantages of good robustness and high anti-interference [26, 27]. It can be described as follows:(10)Hn=−∑π=1n!pπlnpπ,wherepπ is the probability statistics in the entire time series, n is an n-dimensional vector making, and π is the different permutation ways.To relate the sEMG and MMG features with the joint motion model, the following second-order polynomial was used as the fitting function:(11)yku=c0u+c1u⋅θ˙k+c2u⋅θk+c3u⋅θ˙k2+c4u⋅θk2+c5u⋅θ˙k⋅θk,where u=1,2, ciu(i = 0, 1, …, 5) are constant parameters, yk1 is the WA of the sEMG and MMG signals, and yk2 is the PE of sEMG and MMG signals.By combining the equations above, the nonlinear expression of the model was obtained:(12)xk+1=fxk,ak+ωk,yk+1=hxk+1+υk+1,fxk,ak=s0+s1⋅θk+s2⋅θk2⋅ak+s3⋅es4⋅θk−s5⋅sinθkθ˙k+θ¨k⋅Tθk+θ˙k⋅T,hxk=c01+c11⋅θ˙k+c21⋅θk+c31⋅θ˙k2+c41⋅θk2+c51⋅θ˙k⋅θkc02+c12⋅θ˙k+c22⋅θk+c32⋅θ˙k2+c42⋅θk2+c52⋅θ˙k⋅θk,where ωk is the process noise, υk is the measurement noise, ak represents the muscle activation ak, xk=θ¨kθ˙kθkT, and yk=yk1yk2T.


#### 2.3.4. Unscented Particle Filtering (UPF)

The UPF consists of the generic particle filter (PF) and the unscented Kalman filter (UKF). The UPF first uses the UKF to generate the proposal distribution for the PF, and the rest of the steps are the same as PF [28]. The steps to generate the proposal distribution using UKF are as follows:Particles are updated with the UKF using the following equations:(13)X0=x¯,Xi=x¯+nx+λPxii=1,…,nx,Xi=x¯−nx+λPxii=nx,…,2nx,W0m=λnx+λ,W0c=W0m+1−α2+β,Wim=12⋅n+λ,λ=α2nx+κ−nx,where nx is the dimension of x; x¯ and Px are the mean and covariance of x, respectively; κ is the scaling parameter; and α, β are the control parameters.Sigma points:(14)Xt−1a=x¯t−1ax¯t−1a±na+λPt−1a.Time update:(15)Xt∣t−1x=fXt−1x,Xt−1v,x¯t∣t−1=∑i=02naWimXi,t∣t−1x,Yt∣t−1=hXt−1x,Xt−1n,y¯t∣t−1=∑i=02naWimYi,t∣t−1x,Pt∣t−1=∑i=02naWicXi,t∣t−1x−x¯t∣t−1Xi,tt−1x−x¯t∣t−1T,whereYi=gXi and y=gx is a nonlinear transformation.Measurement update:(16)Pytyt=∑i=02naWicYi,t∣t−1−y¯t∣t−1Yi,t∣t−1−y¯t∣t−1T,Pxtyt=∑i=02naWicXi,t∣t−1x−x¯t∣t−1Yi,t∣t−1x−y¯t∣t−1T,Kt=PxtytPytyt−1,x¯t=x¯t∣t−1+Ktyt−y¯t∣t−1,Pt=Pt∣t−1−KtPytytKtT.Then, we substituted the variables of the nonlinear expression with the update equation through calculation and simulation.


## 2.1. Data Acquisition

The muscles relevant to elbow and shoulder joint motion were selected, including the biceps brachii (BB), triceps brachii (TB), brachioradialis (BR), trapezius (TZ), teres minor (TM), anterior deltoid (AD), lateral deltoid (LD), and pectoralis major (PM). The sEMG and MMG signals were measured by placing eight surface electrodes on the muscles of the right upper limb (Figure2).Figure 2
sEMG sensor arrangement (Ch1: BB, Ch2: TB, Ch3: BR, Ch4: TZ, Ch5: TM, Ch6: AD, Ch7: LD, and Ch8: PM).The Trigno™ Wireless (Delsys Inc, Natick, MA, USA) was used to record the sEMG and MMG signals at 1600 Hz. The bandwidth was 20–450 Hz, and the baseline noise was less than 1.25μV. The sensors were positioned following the SENIAM recommendation. The reference electrode was attached to the lateral elbow epicondyle.The Codamotion system (Charnwood Dynamics Ltd., UK) was used to record the real joint angle information. This system can capture and analyze motion information in real time through infrared capture. Figure3 shows the Codamotion system and placement of markers.Figure 3
Codamotion system and placement of markers.Ten healthy subjects (five males and five females, 23 ± 2 years old, weighing 50–70 kg) were tested in the experiments. The study was approved by the Hangzhou Dianzi University’s Institutional Review Board, and all subjects provided a written informed consent prior to participation. Figure4 shows the experimental movements. The subject (1) sits on a chair and keeps the muscles relaxed, (2) straightens the right arm forward with elbow extension, (3) performs elbow flexion, (4) straightens the arm upward with elbow extension, (5) straightens the arm to the right side with elbow extension, (6) performs elbow flexion, and (7) returns the upper arm to the initial position while maintaining elbow flexion. The experimental movements involved the main activities of the upper limb, thereby enabling us to better analyze the movement of the shoulder and elbow joints. Each subject completed all movements with steps 1-2-3-4-5-6-7-1 in 15 s and performed the sequence with or without a load (3 kg). The subjects were given an appropriate rest time after each round to reduce muscle fatigue. We completed the data collection within five days, and the subjects repeated the procedure 10 times in each experiment day. In this way, we were able to collect the required data.Figure 4
Upper limb movements recorded: steps 1-2-3-4-5-6-7-1. Schematic of the (a) right-hand side and (b) back side.


## 2.2. Data Preprocessing

Both sEMG and MMG signals are very weak, nonstationary, and random. Many kinds of noise exist in the acquisition process. The sEMG and MMG signals must be denoised for further application. In this experiment, the Trigno Wireless (Delsys Inc.) with a bandpass filter (sEMG and MMG signals of 20–450 and 5–150 Hz, respectively) was used to eliminate most of the inherent noises and motion artifacts.The sampling frequency of the joint angle signal is 100 Hz, and that of the sEMG and MMG signals is 1600 Hz. To keep the sEMG and MMG signals consistent with the joint angle signal, subsampling of the sEMG signals was implemented using feature extraction. Afterward, the subsampling frequency of the sEMG and MMG signals was changed from 1600 Hz to 100 Hz.


## 2.3. sEMG-MMG State-Space Model

### 2.3.1. Hill-Based Muscle Model

Figure5 illustrates a simplified HMM. The HMM consists of a contractile element and a parallel elastic element. The contractile element produces the active muscle force FAm, and the parallel elastic element produces the passive force FPm. The model can be described as follows:(1)Fm=FPm+FAm=fPl⋅F0m+fAl⋅fVv⋅ak⋅F0m,l=lml0m,v=vmv0m,where Fmt is the musculotendon force, fPl is the normalized passive force-length relationship, F0m is the maximum isometric force, l is the normalized muscle fiber length, lm is the current muscle fiber length, l0m is the optimal fiber length, fAl is the normalized active force-length relationship, fVv is the normalized force-velocity relationship, v is the normalized muscle fiber velocity, vm is the muscle contraction velocity, v0m is the maximum muscle contraction velocity, and ak is the muscle activation [18].Figure 5
Hill-based muscle model.According to the abovementioned equations, the musculotendon force can be described as follows:(2)Fmt=fAl⋅fVv⋅ak+fPl⋅F0m⋅cosϕ,where Fmt is the musculotendon force and ϕ is the pennation angle.The simplifications offAl, fVv, fPl, F0m, lt, l0m, and ϕ can replace the complex biomechanical parameters [19–22]. The musculotendon length can also be simplified as a first-order polynomial [23, 24]. Thus,(3)lmt=b0+b1⋅θ,where lmt is the musculotendon length, b0 and b1 are constants, and θ is the joint angle. Then, the moment arm r and the joint moment τ can be calculated as follows:(4)r=∂lmtθ∂θ=b1,τ=Fmt⋅r.


### 2.3.2. sEMG-MMG Joint Motion Model

The angular acceleration of the joint can be calculated as follows:(5)θ¨=1Ie⋅τ−τeg,τeg=τegm⋅sinθ,where Ie is the moment of inertia, τeg is the external and gravity torque, and τegm is the maximum of τeg.By combining the above equations, the angular accelerationθ¨ can be described as follows:(6)θ¨=s0+s1⋅θ+s2⋅θ2⋅ak+s3⋅es4⋅θ−s5⋅sinθ,where si (i = 0, 1, …, 5) are constants computed as follows:(7)s0=k0⋅F0m⋅b1⋅cosϕIe+k1⋅F0m⋅b0−lt⋅b1l0m⋅Ie+k2⋅F0m⋅b0−lt2⋅b1l0m2⋅Ie⋅cosϕ,s1=k1⋅F0m⋅b12l0m⋅Ie+2⋅k2⋅b0−lt⋅b12l0m2⋅Ie⋅cosϕ,s2=k2⋅F0m⋅b13l0m2⋅Ie⋅cosϕ,s3=F0m⋅b1⋅cosϕIe⋅exp10⋅b0−lt⋅b12l0m⋅cosϕ−15,s4=10⋅b1l0m⋅cosϕ,s5=τegmIe.Then, we can obtain the sEMG-MMG joint motion model in discrete time as follows:(8)θ¨k+1=s0+s1⋅θk+s2⋅θk2⋅ak+s3⋅es4⋅θk−s5⋅sinθk,θ˙k+1=θ˙k+θ¨k⋅Ts,θk+1=θk+θ˙k⋅Ts,where θ˙k is the joint angular velocity, θk is the joint position, and Ts is the sampling time.


### 2.3.3. sEMG-MMG Features and State-Space Model

The simplified sEMG–MMG joint motion model involves two problems. First, the model omits certain physiological parameters. Second, the biomechanical parameters vary because of different body conditions. These problems cause accumulative errors in the recursive calculation of the joint angle. To eliminate these errors, we constructed a measurement equation providing the joint angle measurements as feedback, which then leads to a sEMG-MMG state-space model. We built the measure equation by Wilson amplitude (WA) and permutation entropy (PE) in this study.(1)
WA:The WA describes the different characteristics among adjacent amplitudes [25]. It can be described as follows:(9)WA=∑i=1Nuxi+1−xi−T,wherexi and xi+1 are adjacent samples, and the threshold T is 0.05 V in this study.(2)
PE:The PE describes the complexity of the data and has the advantages of good robustness and high anti-interference [26, 27]. It can be described as follows:(10)Hn=−∑π=1n!pπlnpπ,wherepπ is the probability statistics in the entire time series, n is an n-dimensional vector making, and π is the different permutation ways.To relate the sEMG and MMG features with the joint motion model, the following second-order polynomial was used as the fitting function:(11)yku=c0u+c1u⋅θ˙k+c2u⋅θk+c3u⋅θ˙k2+c4u⋅θk2+c5u⋅θ˙k⋅θk,where u=1,2, ciu(i = 0, 1, …, 5) are constant parameters, yk1 is the WA of the sEMG and MMG signals, and yk2 is the PE of sEMG and MMG signals.By combining the equations above, the nonlinear expression of the model was obtained:(12)xk+1=fxk,ak+ωk,yk+1=hxk+1+υk+1,fxk,ak=s0+s1⋅θk+s2⋅θk2⋅ak+s3⋅es4⋅θk−s5⋅sinθkθ˙k+θ¨k⋅Tθk+θ˙k⋅T,hxk=c01+c11⋅θ˙k+c21⋅θk+c31⋅θ˙k2+c41⋅θk2+c51⋅θ˙k⋅θkc02+c12⋅θ˙k+c22⋅θk+c32⋅θ˙k2+c42⋅θk2+c52⋅θ˙k⋅θk,where ωk is the process noise, υk is the measurement noise, ak represents the muscle activation ak, xk=θ¨kθ˙kθkT, and yk=yk1yk2T.


### 2.3.4. Unscented Particle Filtering (UPF)

The UPF consists of the generic particle filter (PF) and the unscented Kalman filter (UKF). The UPF first uses the UKF to generate the proposal distribution for the PF, and the rest of the steps are the same as PF [28]. The steps to generate the proposal distribution using UKF are as follows:Particles are updated with the UKF using the following equations:(13)X0=x¯,Xi=x¯+nx+λPxii=1,…,nx,Xi=x¯−nx+λPxii=nx,…,2nx,W0m=λnx+λ,W0c=W0m+1−α2+β,Wim=12⋅n+λ,λ=α2nx+κ−nx,where nx is the dimension of x; x¯ and Px are the mean and covariance of x, respectively; κ is the scaling parameter; and α, β are the control parameters.Sigma points:(14)Xt−1a=x¯t−1ax¯t−1a±na+λPt−1a.Time update:(15)Xt∣t−1x=fXt−1x,Xt−1v,x¯t∣t−1=∑i=02naWimXi,t∣t−1x,Yt∣t−1=hXt−1x,Xt−1n,y¯t∣t−1=∑i=02naWimYi,t∣t−1x,Pt∣t−1=∑i=02naWicXi,t∣t−1x−x¯t∣t−1Xi,tt−1x−x¯t∣t−1T,whereYi=gXi and y=gx is a nonlinear transformation.Measurement update:(16)Pytyt=∑i=02naWicYi,t∣t−1−y¯t∣t−1Yi,t∣t−1−y¯t∣t−1T,Pxtyt=∑i=02naWicXi,t∣t−1x−x¯t∣t−1Yi,t∣t−1x−y¯t∣t−1T,Kt=PxtytPytyt−1,x¯t=x¯t∣t−1+Ktyt−y¯t∣t−1,Pt=Pt∣t−1−KtPytytKtT.Then, we substituted the variables of the nonlinear expression with the update equation through calculation and simulation.


## 2.3.1. Hill-Based Muscle Model

Figure5 illustrates a simplified HMM. The HMM consists of a contractile element and a parallel elastic element. The contractile element produces the active muscle force FAm, and the parallel elastic element produces the passive force FPm. The model can be described as follows:(1)Fm=FPm+FAm=fPl⋅F0m+fAl⋅fVv⋅ak⋅F0m,l=lml0m,v=vmv0m,where Fmt is the musculotendon force, fPl is the normalized passive force-length relationship, F0m is the maximum isometric force, l is the normalized muscle fiber length, lm is the current muscle fiber length, l0m is the optimal fiber length, fAl is the normalized active force-length relationship, fVv is the normalized force-velocity relationship, v is the normalized muscle fiber velocity, vm is the muscle contraction velocity, v0m is the maximum muscle contraction velocity, and ak is the muscle activation [18].Figure 5
Hill-based muscle model.According to the abovementioned equations, the musculotendon force can be described as follows:(2)Fmt=fAl⋅fVv⋅ak+fPl⋅F0m⋅cosϕ,where Fmt is the musculotendon force and ϕ is the pennation angle.The simplifications offAl, fVv, fPl, F0m, lt, l0m, and ϕ can replace the complex biomechanical parameters [19–22]. The musculotendon length can also be simplified as a first-order polynomial [23, 24]. Thus,(3)lmt=b0+b1⋅θ,where lmt is the musculotendon length, b0 and b1 are constants, and θ is the joint angle. Then, the moment arm r and the joint moment τ can be calculated as follows:(4)r=∂lmtθ∂θ=b1,τ=Fmt⋅r.


## 2.3.2. sEMG-MMG Joint Motion Model

The angular acceleration of the joint can be calculated as follows:(5)θ¨=1Ie⋅τ−τeg,τeg=τegm⋅sinθ,where Ie is the moment of inertia, τeg is the external and gravity torque, and τegm is the maximum of τeg.By combining the above equations, the angular accelerationθ¨ can be described as follows:(6)θ¨=s0+s1⋅θ+s2⋅θ2⋅ak+s3⋅es4⋅θ−s5⋅sinθ,where si (i = 0, 1, …, 5) are constants computed as follows:(7)s0=k0⋅F0m⋅b1⋅cosϕIe+k1⋅F0m⋅b0−lt⋅b1l0m⋅Ie+k2⋅F0m⋅b0−lt2⋅b1l0m2⋅Ie⋅cosϕ,s1=k1⋅F0m⋅b12l0m⋅Ie+2⋅k2⋅b0−lt⋅b12l0m2⋅Ie⋅cosϕ,s2=k2⋅F0m⋅b13l0m2⋅Ie⋅cosϕ,s3=F0m⋅b1⋅cosϕIe⋅exp10⋅b0−lt⋅b12l0m⋅cosϕ−15,s4=10⋅b1l0m⋅cosϕ,s5=τegmIe.Then, we can obtain the sEMG-MMG joint motion model in discrete time as follows:(8)θ¨k+1=s0+s1⋅θk+s2⋅θk2⋅ak+s3⋅es4⋅θk−s5⋅sinθk,θ˙k+1=θ˙k+θ¨k⋅Ts,θk+1=θk+θ˙k⋅Ts,where θ˙k is the joint angular velocity, θk is the joint position, and Ts is the sampling time.


## 2.3.3. sEMG-MMG Features and State-Space Model

The simplified sEMG–MMG joint motion model involves two problems. First, the model omits certain physiological parameters. Second, the biomechanical parameters vary because of different body conditions. These problems cause accumulative errors in the recursive calculation of the joint angle. To eliminate these errors, we constructed a measurement equation providing the joint angle measurements as feedback, which then leads to a sEMG-MMG state-space model. We built the measure equation by Wilson amplitude (WA) and permutation entropy (PE) in this study.(1)
WA:The WA describes the different characteristics among adjacent amplitudes [25]. It can be described as follows:(9)WA=∑i=1Nuxi+1−xi−T,wherexi and xi+1 are adjacent samples, and the threshold T is 0.05 V in this study.(2)
PE:The PE describes the complexity of the data and has the advantages of good robustness and high anti-interference [26, 27]. It can be described as follows:(10)Hn=−∑π=1n!pπlnpπ,wherepπ is the probability statistics in the entire time series, n is an n-dimensional vector making, and π is the different permutation ways.To relate the sEMG and MMG features with the joint motion model, the following second-order polynomial was used as the fitting function:(11)yku=c0u+c1u⋅θ˙k+c2u⋅θk+c3u⋅θ˙k2+c4u⋅θk2+c5u⋅θ˙k⋅θk,where u=1,2, ciu(i = 0, 1, …, 5) are constant parameters, yk1 is the WA of the sEMG and MMG signals, and yk2 is the PE of sEMG and MMG signals.By combining the equations above, the nonlinear expression of the model was obtained:(12)xk+1=fxk,ak+ωk,yk+1=hxk+1+υk+1,fxk,ak=s0+s1⋅θk+s2⋅θk2⋅ak+s3⋅es4⋅θk−s5⋅sinθkθ˙k+θ¨k⋅Tθk+θ˙k⋅T,hxk=c01+c11⋅θ˙k+c21⋅θk+c31⋅θ˙k2+c41⋅θk2+c51⋅θ˙k⋅θkc02+c12⋅θ˙k+c22⋅θk+c32⋅θ˙k2+c42⋅θk2+c52⋅θ˙k⋅θk,where ωk is the process noise, υk is the measurement noise, ak represents the muscle activation ak, xk=θ¨kθ˙kθkT, and yk=yk1yk2T.


## 2.3.4. Unscented Particle Filtering (UPF)

The UPF consists of the generic particle filter (PF) and the unscented Kalman filter (UKF). The UPF first uses the UKF to generate the proposal distribution for the PF, and the rest of the steps are the same as PF [28]. The steps to generate the proposal distribution using UKF are as follows:Particles are updated with the UKF using the following equations:(13)X0=x¯,Xi=x¯+nx+λPxii=1,…,nx,Xi=x¯−nx+λPxii=nx,…,2nx,W0m=λnx+λ,W0c=W0m+1−α2+β,Wim=12⋅n+λ,λ=α2nx+κ−nx,where nx is the dimension of x; x¯ and Px are the mean and covariance of x, respectively; κ is the scaling parameter; and α, β are the control parameters.Sigma points:(14)Xt−1a=x¯t−1ax¯t−1a±na+λPt−1a.Time update:(15)Xt∣t−1x=fXt−1x,Xt−1v,x¯t∣t−1=∑i=02naWimXi,t∣t−1x,Yt∣t−1=hXt−1x,Xt−1n,y¯t∣t−1=∑i=02naWimYi,t∣t−1x,Pt∣t−1=∑i=02naWicXi,t∣t−1x−x¯t∣t−1Xi,tt−1x−x¯t∣t−1T,whereYi=gXi and y=gx is a nonlinear transformation.Measurement update:(16)Pytyt=∑i=02naWicYi,t∣t−1−y¯t∣t−1Yi,t∣t−1−y¯t∣t−1T,Pxtyt=∑i=02naWicXi,t∣t−1x−x¯t∣t−1Yi,t∣t−1x−y¯t∣t−1T,Kt=PxtytPytyt−1,x¯t=x¯t∣t−1+Ktyt−y¯t∣t−1,Pt=Pt∣t−1−KtPytytKtT.Then, we substituted the variables of the nonlinear expression with the update equation through calculation and simulation.


## 3. Experimental Results

Ten subjects underwent two different situations, namely, loads of 3 kg and no load. Figure6 shows the eight channels of the sEMG signals and joint motion recorded from a typical subject without load. Figure 7 depicts the eight channels of the MMG signals and joint motion recorded from a typical subject with load. Figure 8 illustrates the sEMG and MMG features from a typical subject.Figure 6
(a–h) sEMG signals from BB, TB, BR, TZ, TM, AD, LD, and PM. (i) Angle signal of shoulder joint. (j) Angle signal of elbow joint motion.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(j)Figure 7
(a–h) MMG signals from BB, TB, BR, TZ, TM, AD, LD, and PM. (i) Angle signal of shoulder joint. (j) Angle signal of elbow joint motion.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(j)Figure 8
(a) sEMG signal from TM (A), WA (B), and PE (C). (b) MMG signals from AD (A), WA (B), and PE (C).
(a)
(b)Table1 shows the identified values of the sEMG-MMG state-space model. The parameters are different in load and no-load conditions.Table 1
Identified parameters of the sEMG-MMG state-space model.
s
i
s
0
s
1
s
2
s
3
s
4
s
5
Load
5.26
−8.57
−5.37
0.02
−0.13
0.58
No load
1.88
−1.97
−0.56
0.08
−0.68
0.24
c
i
1
c
0
1
c
1
1
c
2
1
c
3
1
c
4
1
c
5
1
Load
−2.67
−0.07
4.43
−0.24
−1.59
0.02
No load
−0.39
0.01
1.57
0.12
−0.29
0.24
c
i
2
c
0
2
c
1
2
c
2
2
c
3
2
c
4
2
c
5
2
Load
−6.43
0.22
8.35
0.03
−2.54
−0.04
No load
−2.34
0.23
8.69
0.48
−2.41
−0.15The identified parameters were substituted into the sEMG-MMG state-space model above to directly estimate the joint angle through the sEMG-MMG features and closed-loop prediction approaches.The feature dataset of the multijoint movement was inputted into the proposed model. Five-fold cross validation was used to perform the simulation. Figure9 shows the real and estimated joint angles from a typical subject. Table 2 lists the averaged results of the five subdatasets.Figure 9
(a) Real and estimated shoulder joint angles without load. (b) Real and estimated shoulder joint angles with load. (c) Real and estimated elbow joint angles without load. (d) Real and estimated elbow joint angle with load. Red indicates estimated value; blue indicates Codamotion measurement.
(a)
(b)
(c)
(d)Table 2
Evaluation of the performance of the proposed model (RMSE: root mean square error, CC: correlation coefficient).
RMSE
CC
Elbow joint without load
5.86 ± 2.32
0.95 ± 0.02
Elbow joint with load
5.27 ± 2.03
0.96 ± 0.02
Shoulder joint without load
7.89 ± 1.68
0.92 ± 0.03
Shoulder joint with load
7.26 ± 1.73
0.93 ± 0.02As shown in Figure9 and Table 2, the mean RMSE and CC of the elbow joint with load are 5.27 and 0.96, respectively, which indicates that the sEMG-MMG state-space model successfully kept the estimation convergent. Moreover, the estimation errors of the elbow joint are smaller than those of the shoulder joint, and the load successfully improved the accuracy.Five representative regression techniques, namely, the backpropagation neural network (BPNN) [29], generalized regression neural network (GRNN) [30], support vector regression (SVR) [31], Fisher linear discriminant analysis (FDA) [32], and fuzzy min-max neural network (FMMNN) [33] were compared to verify the effectiveness of the proposed model. The same sEMG and MMG signals and features (WA and PE) and training set of neural networks were used in the experiments. The results are shown in Figure 10 and Table 3.Figure 10
Mean RMSE and CC of seven regression models.
(a)
(b)Table 3
Evaluation of the performance of the proposed model (RMSE: root mean square error, CC: correlation coefficient).
RMSE
CC
Proposed model
6.27 ± 0.63
0.95 ± 0.02
HMM
13.42 ± 2.59
0.79 ± 0.09
BPNN
8.89 ± 1.03
0.92 ± 0.03
GRNN
8.58 ± 0.73
0.89 ± 0.04
SVR
8.43 ± 0.82
0.91 ± 0.05
FDA
9.07 ± 1.03
0.93 ± 0.02
FMMNN
7.29 ± 0.78
0.93 ± 0.01As shown in Figure10 and Table 3, the proposed model exhibits the best performance in the experiments. The mean RMSE and CC of the proposed model are 6.27 and 0.95, respectively. The errors of the proposed model are much smaller than those of the other regression methods.


## 4. Discussion

This study aimed to build an sEMG-MMG state-space model to accurately estimate continuous joint movements from sEMG signals. The results showed that the proposed model improved the estimation accuracy of elbow-joint and shoulder-joint motions, with mean RMSEs of less than 5.3 and 7.3, respectively, indicating that the sEMG-MMG state-space model successfully kept the estimation convergent. Moreover, the proposed model performed better in the experiments than the conventional regression methods. In the comparative experiment, all of these methods were limited by data size, causing the neural networks to perform poorer than the sEMG-MMG state-space model because the training of neural networks requires more data sets.Several basic limitations were associated with this study that need further development to be used in rehabilitation robots under real-world conditions. sEMG and MMG signals with specific motions were collected from healthy volunteers in a laboratory environment. Therefore, the performance of this model in real scenario with unscripted free-form activities performed by elderly or actual patients is unknown. Although no considerable difference in the inherent characteristics of the sEMG and MMG signals was found between subjects with and without disabilities, the amplitude and frequency of the signal still influenced the estimation accuracy. In this study, each individual was trained separately and required multiple repetitions of the task to obtain sufficient data for training and testing purposes. However, its difference from joint motions in real life with a larger task set is unclear. Also, the sEMG-MMG state-space model ignores some complex physiological parameters of HMM, which can result in prediction errors. These conditions need to be investigated before using the proposed model for clinical purposes.


## 5. Conclusions

In this study, a sEMG-MMG state-space model for multijoint motion estimation was constructed, where the continuous angular displacement could be estimated using the sEMG and MMG signals. Extensive experiments and comparison between the estimations and Codamotion measurements demonstrated the effectiveness of the proposed model. Compared with traditional methods for angle estimation, the proposed technique improved the estimation accuracy and real-time characteristics. In the future, the proposed method could be utilized for the training of rehabilitation robots and exoskeleton robots.



---

*Source: 4503271-2020-02-11.xml*