**1. Introduction**

As is the case for most of the potential technical equipment for rehabilitation training and movement assistance, the exoskeleton system has always attracted attention in related research fields for the elderly and the disabled [1,2]; however, with some key issues not being resolved, the performance of this wearable robot remains relatively limited [3]. Among them, human intention perception belongs to a pretty critical research point and needs to be studied further. Traditional methods try to achieve this function by monitoring limb kinematics data [4] or human-machine interaction information [5], but problems such as response lag greatly restrict its actual effect. In recent years, recognition methods based on biological signals have emerged and gained wide attention [6], which are expected to realize more effective intention understanding.

Biological signals are often generated before the execution of corresponding actions, which show a certain degree of motion predictability and can make intention recognition more timely and accurate. In addition to eye tracking and galvanic skin response (GSR) [7], which are not suitable for combining with exoskeleton control technology, commonly used types in current research mainly contain electroencephalogram (EEG) [8], surface electromyogram (sEMG) [9], and mechanomyography [10]. The EEG signal originates from the potential of the external electrical field that fluctuates around nerve cells in the brain, and is often applied to classify specific movement patterns [11,12]. Due to its instability and susceptibility to interference, EEG-based methods for data collection, signal processing, and intent identification needs to be further improved and optimized. The sEMG signal

**Citation:** Shi, Y.; Dong, W.; Lin, W.; He, L.; Wang, X.; Li, P.; Gao, Y. Human Joint Torque Estimation Based on Mechanomyography for Upper Extremity Exosuit. *Electronics* **2022**, *11*, 1335. https://doi.org/ 10.3390/electronics11091335

Academic Editors: João Paulo Morais Ferreira; Tao Liu and Nicola Francesco Lopomo

Received: 25 February 2022 Accepted: 20 April 2022 Published: 22 April 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

represents the total action potentials of different motor units innervated by certain motor neurons, which can be obtained through electrodes placed on the corresponding muscle tissue [13]. It has the advantage of high sensitivity and low delay, but many factors will diminish its effectiveness for data collection, such as skin surface cleanliness, air humidity, and patch electrode position [14]. MMG reflects the mechanical vibration signal generated when the skeletal muscle contracts, and it contains abundant information, such as the number of muscle fibers participating in the exercise, and the amplitude and frequency of vibration [15]. Compared with the above two signals, it has the following advantages [16].


However, there are still some limitations in the practical application of MMG. For example, it is easily contaminated by low-frequency motion artifacts, and sensitive to sudden step noises; therefore, a suitable signal processing algorithm is needed to extract a relatively pure sequence.

At present, classification and regression algorithms have been widely used in many studies, such as the trajectory control of a redundant robot [17] and hand gesture recognition for teleoperated surgical robot systems [18], which can also complete the analysis of human motion intention. As a relatively common and easy-to-use method, motion pattern recognition [19,20] aims at establishing the mapping relationship from sensor information to finite human motion states through classification models. Then, the control subsystem will generate corresponding commands according to current motion pattern and will deliver them to underlying drivers. Joint angle prediction attempts to calculate the limb position at the next moment based on regression models [21,22], which can effectively avoid the response lag of the exosuit. After that, the motion state of the power-assisted system can be dynamically adjusted through position closed-loop control. Joint torque estimation also belongs to a direct and effective method of obtaining intentions [23,24]. The desired torque can be calculated through the Hill type model [25] or obtained from the biological signals using machine learning algorithms [26], and can act as the input parameter of torque closed-loop control to regulate the motor output. The last method will provide a reference for the direct control of external torque that assists joint motion, which is quite practical for upper extremity exosuits that need to achieve an expected power-assisted efficiency.

If we want to use the abovementioned MMG signal to estimate the joint torque, it is obviously quite difficult through conventional mathematical derivation. The machine learning algorithm can train the mapping model very well based on the existing data, and it shows an excellent fitting ability in many research fields, such as breathing pattern detection [27] and human activity identification [28]; therefore, this method should be able to describe the complex and nonlinear relationship between MMG and joint torques.

The inherent characteristics of the soft exosuit based on Bowden cables greatly increases the difficulty for designing control strategies [29,30]. The gravity compensation algorithm is a simple and common method for motion generation, which will output an active torque to offset the joint load imposed by the limb weight [31,32]; however, due to the large error of compensation model and the ignorance of dynamic characteristics, it may cause the system response to, more or less, mismatch the joint movement. The exosuit named CRUX [33] tries to control the target arm to follow the reference trajectory of the healthy one, consequently completing the active rehabilitation training process [34]. This method will limit the subjective initiative of the wearer to a certain extent, and is not suitable for situations where both arms need assistance. Some scholars from Italy have proposed a threshold method based on sEMG [35]. When the signal amplitude of the wearer's measured muscle exceeds the set value, the wearable system starts to pro-

duce a power-assisted effect for motion immediately. The most obvious disadvantage of this method is that it only outputs a constant driving force, but the required joint torque changes dynamically under different motion states. In general, the control logic for the upper extremity exosuit still has some defects, and remains to be further explored. power−assisted effect for motion immediately. The most obvious disadvantage of this method is that it only outputs a constant driving force, but the required joint torque changes dynamically under different motion states. In general, the control logic for the upper extremity exosuit still has some defects, and remains to be further explored.

proposed a threshold method based on sEMG [35]. When the signal amplitude of the wearer's measured muscle exceeds the set value, the wearable system starts to produce a

*Electronics* **2022**, *11*, x FOR PEER REVIEW 3 of 15

In this paper, we intend to complete data preprocessing and feature extraction using MMG, then establish the mapping relationship from this signal to the joint torque based on the machine learning algorithm, and finally apply it to the rehabilitation training research of upper extremity exosuit. MMG is collected through an inertial measurement unit (IMU) and synthesized by the linear accelerations along three axes. The HHT will filter this original signal to obtain a relatively pure result. We extract three features from the processed data and combine them with the collected joint torque to form data sets for training and testing. A RFR model is designed as the algorithm framework for joint torque estimation, and its parameters are determined through offline training. According to the above research results, a control algorithm based on joint torque estimation will take charge of the motion control for an upper extremity exosuit. Eventually, some experiments are carried out to test the accuracy of torque estimation and the efficiency of power assistance. The main contributions and highlights of this study are summarized as follows. In this paper, we intend to complete data preprocessing and feature extraction using MMG, then establish the mapping relationship from this signal to the joint torque based on the machine learning algorithm, and finally apply it to the rehabilitation training research of upper extremity exosuit. MMG is collected through an inertial measurement unit (IMU) and synthesized by the linear accelerations along three axes. The HHT will filter this original signal to obtain a relatively pure result. We extract three features from the processed data and combine them with the collected joint torque to form data sets for training and testing. A RFR model is designed as the algorithm framework for joint torque estimation, and its parameters are determined through offline training. According to the above research results, a control algorithm based on joint torque estimation will take charge of the motion control for an upper extremity exosuit. Eventually, some experiments are carried out to test the accuracy of torque estimation and the efficiency of power assistance. The main contributions and highlights of this study are summarized as follows.


The remaining research content of this article is organized as follows. Section 2 describes the measurement and calculation methods for MMG, its corresponding joint torque, and how to process the original signal to construct data sets. Section 3 introduces the design details of RFR model and uses a large amount of test data to fit the desired quantitative relationship. In Section 4, we have developed a control strategy for the upper extremity exosuit based on torque estimation. Section 5 proves the feasibility of the above methods through some experiments. Section 6 is the conclusion. The remaining research content of this article is organized as follows. Section 2 describes the measurement and calculation methods for MMG, its corresponding joint torque, and how to process the original signal to construct data sets. Section 3 introduces the design details of RFR model and uses a large amount of test data to fit the desired quantitative relationship. In Section 4, we have developed a control strategy for the upper extremity exosuit based on torque estimation. Section 5 proves the feasibility of the above methods through some experiments. Section 6 is the conclusion.

#### **2. Data Sets Acquisition 2. Data Sets Acquisition**

#### *2.1. Raw Information Collection 2.1. Raw Information Collection*

In order to obtain the MMG signal and corresponding joint torque at the same time, we built a measurement platform for joint information collection. Figure 1a shows how to use this device to get relevant data about elbow static flexion and extension. In order to obtain the MMG signal and corresponding joint torque at the same time, we built a measurement platform for joint information collection. Figure 1a shows how to use this device to get relevant data about elbow static flexion and extension.

**Figure 1.** Information collection process during elbow static flexion/extension: (**a**) measurement platform; (**b**) joint torque calculation model for elbow.

During this process, a six-dimension force sensor and a high-performance IMU are responsible for detecting forces along three directions at the end of limb, and gathering MMG at the brachioradialis of forearm, respectively. When subjects hold the force measuring rod and try to perform specific actions, this platform will transmit corresponding data to the host computer through a data line and save them in files. During this process, a six−dimension force sensor and a high−performance IMU are responsible for detecting forces along three directions at the end of limb, and gathering MMG at the brachioradialis of forearm, respectively. When subjects hold the force measuring rod and try to perform specific actions, this platform will transmit corresponding data to the host computer through a data line and save them in files.

**Figure 1.** Information collection process during elbow static flexion/extension: (**a**) measurement

*Electronics* **2022**, *11*, x FOR PEER REVIEW 4 of 15

platform; (**b**) joint torque calculation model for elbow.

With readings from a six-dimension force sensor (*F<sup>x</sup>* and *Fy*) as input, the joint torque ( → *Telbow*) can be calculated by a certain mathematical relation. The mechanical model of the elbow is shown in Figure 1b, where → *A*, → *F*, and *L*, respectively, represent the current posture vector of human arm, the resultant force vector at the end, and the length of forearm. Then, we can dynamically obtain the joint torque values during elbow static flexion and extension through the following formulas. With readings from a six−dimension force sensor (*F<sup>x</sup>* and *Fy*) as input, the joint torque (*T*⃗ *elbow*) can be calculated by a certain mathematical relation. The mechanical model of the elbow is shown in Figure 1b, where , , and *L*, respectively, represent the current posture vector of human arm, the resultant force vector at the end, and the length of forearm. Then, we can dynamically obtain the joint torque values during elbow static flexion and extension through the following formulas.

$$
\stackrel{\rightarrow}{F} = \stackrel{\rightarrow}{F}\_x + \stackrel{\rightarrow}{F}\_y \tag{1}
$$

$$
\stackrel{\rightarrow}{T}\_{\text{elbow}} = L(\frac{\stackrel{\rightarrow}{\dot{A}}}{\stackrel{\rightarrow}{|\dot{A}|}} \times \stackrel{\rightarrow}{F}) \tag{2}
$$

During dynamic flexion/extension without using the platform, it can be seen in Figure 2 that when the measured muscle contracts or relaxes and drives the human joint to rotate, the IMU can perceive some regular acceleration signals in the *x*, *y*, and *z* directions synchronously, and the same is true for the process of static joint output. In order to integrate all the effective information, we calculate the sum of linear acceleration vectors along three axes (*ACCx*, *ACCy*, and *ACCz*), and take its modulus as the original MMG signal (*MMG*(*t*)), which can be expressed as the following equation. During dynamic flexion/extension without using the platform, it can be seen in Figure 2 that when the measured muscle contracts or relaxes and drives the human joint to rotate, the IMU can perceive some regular acceleration signals in the *x*, *y*, and *z* directions synchronously, and the same is true for the process of static joint output. In order to integrate all the effective information, we calculate the sum of linear acceleration vectors along three axes (*ACCx*, *ACCy*, and *ACCz*), and take its modulus as the original MMG signal (*MMG*(*t*)), which can be expressed as the following equation.

$$\text{MMG}(t) = \sqrt{\text{ACC}\_x^2 + \text{ACC}\_y^2 + \text{ACC}\_z^2} \tag{3}$$

**Figure 2.** The IMU data during elbow dynamic flexion/extension with low strength: (**a**) three−axis euler angles; (**b**) three-axis accelerations. **Figure 2.** The IMU data during elbow dynamic flexion/extension with low strength: (**a**) three-axis euler angles; (**b**) three-axis accelerations.

#### *2.2. MMG Signal Processing 2.2. MMG Signal Processing*

However, the method of using linear acceleration values for MMG characterization inevitably mixes low−frequency artifacts into the collected original signal, such as gravitational acceleration, IMU posture changing artifacts caused by muscle deformation, and motion artifacts introduced by upper limb movements. Moreover, high−frequency white noise may also be superimposed on it. In order to improve the quality of data and lay the foundation for subsequent feature extraction, an effective filtering method must be applied to eliminate the abovementioned interferences. Traditional signal processing methods are mostly based on Fourier analysis, but these ones have limited effects in practical However, the method of using linear acceleration values for MMG characterization inevitably mixes low-frequency artifacts into the collected original signal, such as gravitational acceleration, IMU posture changing artifacts caused by muscle deformation, and motion artifacts introduced by upper limb movements. Moreover, high-frequency white noise may also be superimposed on it. In order to improve the quality of data and lay the foundation for subsequent feature extraction, an effective filtering method must be applied to eliminate the abovementioned interferences. Traditional signal processing methods are mostly based on Fourier analysis, but these ones have limited effects in practical applications of processing MMG due to its non-linear and non-stationary characteristics. With reference to related literatures, we decide to use HHT to analyze the original data, which is more suitable for these kinds of signals.

The HHT consists of two parts, namely empirical mode decomposition (EMD) and the Hilbert transform. EMD can decompose a complex signal into a limited number of intrinsic mode functions (IMFs) and a residual based on the local time scale characteristics of itself. The specific implementation steps are as follows. First, we find all the local maximum points and local minimum points of the original MMG signal and fit their respective envelopes through the cubic spline curve. Then, the mean value of the upper and lower envelopes (*U*(*t*) and *L*(*t*)) will be subtracted from the original sequence to get the remaining part with the low frequency information removed (*x*(*t*)). The HHT consists of two parts, namely empirical mode decomposition (EMD) and the Hilbert transform. EMD can decompose a complex signal into a limited number of intrinsic mode functions (IMFs) and a residual based on the local time scale characteristics of itself. The specific implementation steps are as follows. First, we find all the local maximum points and local minimum points of the original MMG signal and fit their respective envelopes through the cubic spline curve. Then, the mean value of the upper and lower envelopes (*U*(*t*) and *L*(*t*)) will be subtracted from the original sequence to get the remaining part with the low frequency information removed (*x*(*t*)).

applications of processing MMG due to its non−linear and non−stationary characteristics. With reference to related literatures, we decide to use HHT to analyze the original data,

*Electronics* **2022**, *11*, x FOR PEER REVIEW 5 of 15

which is more suitable for these kinds of signals.

$$\mathbf{x}(t) = \text{MMG}(t) - \frac{1}{2}(\mathcal{U}(t) + L(t))\tag{4}$$

If the number of extreme values and zero crossings on the entire data set of *x*(*t*) differs by 0 or 1, and the average value of its two envelopes remains zero at any time, then *IMF<sup>i</sup>* (*t*) = *x*(*t*), otherwise it is necessary to let *MMG*(*t*) = *x*(*t*), and repeat the above steps until these two conditions are met. Next, we remove the obtained *IMF<sup>i</sup>* (*t*) from *MMG*(*t*) and repeat all the above steps again with the remaining part (*r<sup>i</sup>* (*t*)) to get other IMF components until *r<sup>i</sup>* (*t*) is a constant or monotonic function. As shown in Figure 3a, the original MMG signal is decomposed into 9 IMFs and 1 residual (*res*(*t*)) according to the signal frequency, which can be expressed as follows. If the number of extreme values and zero crossings on the entire data set of *x*(*t*) differs by 0 or 1, and the average value of its two envelopes remains zero at any time, then *IMFi*(*t*) = *x*(*t*), otherwise it is necessary to let *MMG*(*t*) = *x*(*t*), and repeat the above steps until these two conditions are met. Next, we remove the obtained *IMFi*(*t*) from *MMG*(*t*) and repeat all the above steps again with the remaining part (*ri*(*t*)) to get other IMF components until *ri*(*t*) is a constant or monotonic function. As shown in Figure 3a, the original MMG signal is decomposed into 9 IMFs and 1 residual (*res*(*t*)) according to the signal frequency, which can be expressed as follows.

9

$$\text{MMG}(t) = \sum\_{i=1}^{9} \text{IMF}\_i(t) + \text{res}(t) \tag{5}$$

**Figure 3.** MMG signal processing: (**a**) decomposing result through EMD; (**b**) autocorrelation function curves of the first four IMFs. **Figure 3.** MMG signal processing: (**a**) decomposing result through EMD; (**b**) autocorrelation function curves of the first four IMFs.

After completing the above analysis, it is time to select the IMFs dominated by MMG through certain methods and reorganize them to obtain a relatively pure signal. To eliminate the noise−dominated IMFs, we introduce the concept of autocorrelation (*RIMF*(*t*1,*t*2)), which reflects the correlation degree of signal values at different times (*t*<sup>1</sup> and *t*2). Its normalized expression form, (*ρIMF*(*τ*)), can be obtained with the following formula, where *τ* = *t*1−*t*2. After completing the above analysis, it is time to select the IMFs dominated by MMG through certain methods and reorganize them to obtain a relatively pure signal. To eliminate the noise-dominated IMFs, we introduce the concept of autocorrelation (*RIMF*(*t*1,*t*2)), which reflects the correlation degree of signal values at different times (*t*<sup>1</sup> and *t*2). Its normalized expression form, (*ρIMF*(*τ*)), can be obtained with the following formula, where *τ* = *t*<sup>1</sup> − *t*2.

$$R\_{\rm IMF}(t\_{1'}t\_2) = E[\rm IMF(t\_1) \cdot \rm IMF(t\_2)] \tag{6}$$

$$\rho\_{IMF}(\tau) = \frac{\mathcal{R}\_{IMF}(\tau)}{\mathcal{R}\_{\mathcal{X}}(\mathbf{0})} \tag{7}$$

If the normalized autocorrelation curve belongs to an impulse function close to the zero point, it can be ascertained that the corresponding IMF is dominated by noise, because noise has randomness and a weak correlation at every moment. After the calculations in Figure 3b, the first IMF can be classified as such a disturbance. If the normalized autocorrelation curve belongs to an impulse function close to the zero point, it can be ascertained that the corresponding IMF is dominated by noise, because noise has randomness and a weak correlation at every moment. After the calculations in Figure 3b, the first IMF can be classified as such a disturbance.

To exclude the IMFs dominated by motion artifacts, we try to find the discrimination basis from the energy distribution of each order component. After calculating the energy value (*E i IMF*) of each IMF according to the following Formula (8), it is revealed that energy

contained in each IMF has increased significantly, starting from the sixth order. By analyzing the data of different participants, we find that removing these IMFs as motion artifacts can achieve a better result. lyzing the data of different participants, we find that removing these IMFs as motion artifacts can achieve a better result. In the end, we believe that IMF2~IMF5 can effectively characterize the muscle activ-

contained in each IMF has increased significantly, starting from the sixth order. By ana-

To exclude the IMFs dominated by motion artifacts, we try to find the discrimination basis from the energy distribution of each order component. After calculating the energy

) of each IMF according to the following Formula (8), it is revealed that energy

*Electronics* **2022**, *11*, x FOR PEER REVIEW 6 of 15

In the end, we believe that IMF2~IMF5 can effectively characterize the muscle activity, and the filtered signal can be obtained with them being recombined. It can be seen from Figure 4 that there is a pretty clear correspondence between the processed MMG and joint torque. ity, and the filtered signal can be obtained with them being recombined. It can be seen from Figure 4 that there is a pretty clear correspondence between the processed MMG and joint torque. *n*

$$E\_{\rm IMF}^i = \frac{1}{n} \sum\_{j=1}^n \left[ \rm IMF\_i(j) \right]^2 \tag{8}$$

**Figure 4.** Comparison between the filtered MMG of the brachioradialis muscle and the corresponding elbow joint torque: (**a**) the change curve of joint torque; (**b**) the change curve of filtered MMG. **Figure 4.** Comparison between the filtered MMG of the brachioradialis muscle and the corresponding elbow joint torque: (**a**) the change curve of joint torque; (**b**) the change curve of filtered MMG.

#### *2.3. Feature Extraction 2.3. Feature Extraction*

value (*E*

IMF

*i*

Considering that the change in limb strength is often accompanied by the fluctuation of the muscle fiber's vibration amplitude, we select the root mean square (RMS) as the time domain characteristic of MMG. It reflects the effective value of data amplitude and can be calculated with the following formula. Considering that the change in limb strength is often accompanied by the fluctuation of the muscle fiber's vibration amplitude, we select the root mean square (RMS) as the time domain characteristic of MMG. It reflects the effective value of data amplitude and can be calculated with the following formula.

$$\text{RMS}\_{\text{MMG}} = \sqrt{\frac{1}{N} \sum\_{i=1}^{N} X\_i^2} \tag{9}$$

*i* =1 Due to the correlation between muscle activity and its vibration frequency, mean power frequency (MPF) is used to represent the frequency domain characteristic of MMG. It is necessary to perform the Hilbert transform on the filtered MMG to analyze its frequency spectrum, and then integrate it on the time axis to obtain the Hilbert marginal spectrum that characterizes the relationship between signal frequency (*fi*) and energy (*Ei*). Due to the correlation between muscle activity and its vibration frequency, mean power frequency (MPF) is used to represent the frequency domain characteristic of MMG. It is necessary to perform the Hilbert transform on the filtered MMG to analyze its frequency spectrum, and then integrate it on the time axis to obtain the Hilbert marginal spectrum that characterizes the relationship between signal frequency (*f<sup>i</sup>* ) and energy (*E<sup>i</sup>* ). Then, MPF can be calculated through the following equation.

$$MPF\_{MMG} = \frac{\sum\_{i=1}^{N} f\_i E\_i}{\sum\_{i=1}^{N} E\_i} \tag{10}$$

∑ *E<sup>i</sup> N i =*1 The MMG signal may contain information about the number of muscle fibers involved in power output; therefore, we additionally introduce the concept of sample entropy (SampEn) to describe the characteristic from a non−linear perspective, which can The MMG signal may contain information about the number of muscle fibers involved in power output; therefore, we additionally introduce the concept of sample entropy (SampEn) to describe the characteristic from a non-linear perspective, which can quantify the complexity of the time series.

quantify the complexity of the time series. In order to not lose continuity information in the sequence, we apply the sliding window strategy to extract these characteristic values of the filtered MMG signal. The window length and step length are set as 500 ms and 50 ms, respectively. So far, the data set of elbow static flexion/extension is established with RMS, MPF, and SampEn of the MMG signal as features, and joint torque as the label. We can also use similar methods to obtain In order to not lose continuity information in the sequence, we apply the sliding window strategy to extract these characteristic values of the filtered MMG signal. The window length and step length are set as 500 ms and 50 ms, respectively. So far, the data set of elbow static flexion/extension is established with RMS, MPF, and SampEn of the MMG signal as features, and joint torque as the label. We can also use similar methods to obtain relevant information of shoulder static flexion/extension and static adduction/abduction.

#### relevant information of shoulder static flexion/extension and static adduction/abduction. **3. Off-Line Torque Estimation**
