A Discrete-Time Algorithm for Stiffness Extraction from sEMG and Its Application in Antidisturbance Teleoperation

We have developed a new discrete-time algorithm of stiffness extraction from muscle surface electromyography (sEMG) collected from human operator’s arms and have applied it for antidisturbance control in robot teleoperation. The variation of arm stiffness is estimated from sEMG signals and transferred to a telerobot under variable impedance control to imitate human motor control behaviours, particularly for disturbance attenuation. In comparison to the estimation of stiffness from sEMG, the proposed algorithm is able to reduce the nonlinear residual error effect and to enhance robustness and to simplify stiffness calibration. In order to extract a smoothing stiffness enveloping from sEMG signals, two enveloping methods are employed in this paper, namely, fast linear enveloping based on low pass filtering and moving average and amplitude monocomponent and frequency modulating (AM-FM) method. Both methods have been incorporated into the proposed stiffness variance estimation algorithm and are extensively tested. The test results show that stiffness variation extraction based on the two methods is sensitive and robust to attenuation disturbance. It could potentially be applied for teleoperation in the presence of hazardous surroundings or human robot physical cooperation scenarios.


Introduction
Physiological experiments have shown that human arm can be stabilized by human motor control, mainly dependant on mechanical impedance adaptation during interaction with external environment [1]. Such behaviours have received increasing research attention in biomimetic controller design [2][3][4] and have been applied in human robot skill transfer, for example, rehabilitation and manual training.
Conventional intelligent control could deal with system uncertainties given the form of dynamics [5][6][7], while "human-like" controllers could be designed without model information. Human inspired control strategies [3] enable robots to have some human motor features in an economical perspective, and may have great potentials in compliant human robot interactions especially in some physical human robot coupling scenarios, for example, rehabilitation or daily physical cooperation tasks. In addition, spring-like robot joint actuators such as Serial Elastic Actuators (SEA) [8] or Variable Impedance Actuators (VIA) [9] allow researchers to adapt robot impedance parameters physically or programmably [10]. However, it is still difficult to select proper impedance regulations for specific tasks due to complexity and difficulty in human motor precise modelling [11,12].
While there exist different kinds of methods and technologies for human robot skill transfer with various body sensors, for example, vision based kinematic skill transfer from human operator to robot [13][14][15]. Such methods are effective for human robot kinematic skill transfer but difficult to transfer human stiffness regulations which represent human learned skills [16,17]. Another interesting method is impedance learning by demonstration [18] which needs a set of examples for the learning procedure while it may be not efficient for teleoperation because of complexity of environmental modelling. Haptic feedback may enable robot adapt its stiffness in a better human-like manner [19] and  Figure 1: Profile of arm control model and stiffness transfer from human to robot for antidisturbance test [34]: ,̇ ,̈ are planned position, velocity, and acceleration, respectively. is feed-forward force. is feedback force. is the external force. ,̇ ,̈ are output position, velocity, and acceleration, respectively. enhance human robot interaction performance, while it may have limitations as well due to varying environment and different user experience. sEMG based skill transfer from human to robot may be comparatively advantageous to bridge the gap between human and robot interaction, as sEMG signals denote human muscles activations and in proper muscle cocontractions, sEMG signals have an approximately linear correlation with human joint motion, force, stiffness, and so forth [20][21][22][23]. Moreover, sEMG signals are easily accessible to build intuitive human robot interaction interface [24][25][26][27]. It is also a feedforward way to decode human motion intention from sEMG signals [28].
Because of these specific characteristics, sEMG signals have already been employed in teleoperation scenarios. In [26,29,30], a robot teleimpedance system was developed, combined with stiffness estimated from sEMG signals and kinematic motion tracked by vision markers adhered on the human arm.
Efficient sEMG signal processing is essential for sEMG based real time robot control in teleoperation, and discretetime based sEMG enveloping seems to fit for this requirement. In [31], there is a general review for sEMG processing, which indicates that discrete-time filters are more efficient in sEMG signal processing such as Discrete Wavelet Transform (DWT). In [30], sEMG signals are processed to extract incremental stiffness in discrete-time domain to reduce nonlinear stiffness estimation residual errors and calibration time. Its application was tested on robot for disturbance attenuation for both posture control and motion control. However, stiffness estimation in a smooth and robust manner is still a challenging work. For example, human limb stiffness measurement is generally based on mechanical perturbation [32], but it is not suitable for real time human robot interaction, as large errors are generated in human arm continuous stiffness estimation, which may eliminate the performance of telerobot. Alternatively, sEMG based on stiffness estimation could be classified into different types according to muscle cocontraction levels, and several experiments have shown that robot tends to be more stable using such discrete way, which produces a better performance in real-world usage scenarios [33].
In this paper, a profile of teleoperation control approach is developed, based on stiffness variation estimation proposed in our previous work [30]. The whole structure is shown in Figure 1. Human arm control model is simplified as a feedforward force-dependant term and feedback force-dependant term for compensating external force [34] or virtual external force by visual feedback from robot environment. Human operator arm stiffness variation is extracted from sEMG signals during arm antidisturbance tasks. Two sEMG signal enveloping methods are employed and their applications on robot posture maintenance control and path following control have been extensively tested on a Baxter Robot [35], which will be described in detail in the following sections.

Human Arm Stiffness Variation Estimation.
In this paper, only muscle cocontraction force effect on the human arm stiffness variation is considered, except from the effect of arm posture [32], because we mainly focus on the variation of stiffness in different sampling intervals. In addition, arm posture does not change much in the experiment. Therefore, the stiffness variation estimation model is comparatively simplified.
According to the approximate linear relationship between joint force and sEMG signals, human arm joint force estimation via sEMG can be depicted by (1).
∈ and anta ∈ denote antagonist and antagonist of muscle cocontractions, respectively, which can be represented by filtered sEMG signals. ( ) is the nonlinear residual error.
= ℎ , ℎ ∈ 6× is the Jacobian of human operator arm and ∈ 6 is the force exerted on endpoint and can be measured by force sensor. can be extracted from sEMG signals. Thus [ ] can be identified empirically by LS. Correspondingly, human arm joint stiffness can be depicted as (2) according to [20], as follows: denotes the joint stiffness generated by muscles involved, and ὔ ( ) is the nonlinear residual error.
In order to compensate possible nonlinear residual error, we employ a differential method to estimate human arm stiffness variation. So (2) can be modified, as follows: where Δ ℎ = ℎ ( + 1) − ℎ ( ) is stiffness variation between sample instant + 1 and , [ is the enveloping variation extracted from sEMG signals between sample instant +1 and . Note that if Δ ℎ is the vector of endpoint of Cartesian space, it needs to be transformed to robot joint space which will be detailed in the following sections. For simplification, in this paper, we assume Δ ὔ ( ) = 0.

Identification of Human Arm Endpoint Stiffness
Variation Coefficients. In this paper, human arm elbow stiffness variation coefficients estimation is based on human endpoint stiffness measurement. Estimation of force coefficients from sEMG signals is prerequisite for stiffness variation coefficients, as the stiffness variation coefficients are the absolute value of force coefficients as indicated in (2).
Generally, LS is a simple and efficient way to estimate coefficients and has been used widely in human limb stiffness measurements [36,37]. Here, we use (4)-(5) to estimate the stiffness coefficients. During the estimation, subject wears MYO Armband [38] shown in Figure 2 (4 channels for agonist muscles and the other 4 channels for antagonist muscle) on the upper arm near elbow and exerts force (±5 , ±10 , ±15 , ±20 , ±25 ) in , , , with 3 repeatable trails, respectively, with his/her wrist attached to coupling mechanism in order to constrain wrist motion.
The coupling mechanism is properly designed embedded with ATI mini45 force sensor [39] to couple human and robot physically with hand free and wrist motion resistance in human arm endpoint stiffness estimation. The force amplitude and direction are monitored by force presentation interface in the presence of dynamic color bar: In total, we can get differentiated force matrix 9 × 3 × 3 in each coordinate axis and corresponding sEMG signals, which would be normalized to 100% Maximal Voluntary Contraction (MVC) and then differentiated synchronously using MATLAB/Simulink [40]. Stiffness estimation proceeding can be depicted in Figure 3: wherê is estimated coefficients, Δ = [ Δ Δ anta ], and Δ = ℎ ( + 1) − ℎ ( ).

Enveloping of sEMG Signals Based on Fast Smoothing
Algorithm. As mentioned above, sEMG signals are noninvasive and generated in muscle activities, representing the muscle tension, joint force stiffness variation, and so forth [21]. sEMG signals are a linear summation of a compound of motor-unit action potentials (MUAPs) trains triggered by motor units [41]. Normally, 20-500 Hz sEMG signals are extracted for postprocessing [42]. While in this paper, it is enough for us to implement the antidisturbance experiment using sEMG sampling rate at 200 Hz. Additionally, stiffness estimated from sEMG should be smoothed in order to make robot motion compliant. Based on [43,44], a modified fast sEMG envelopig algorithm is developed to extract amplitude of sEMG signals as indicated in (8)-(10) assuming that there is no time delay between procedures: where ( ), = 1, 2, 3, 4, denotes the procedure outputs of sEMG signal smooth enveloping at simple constant ; > 0 is coefficient of squaring sEMG signals; is the length of sEMG array at simple constant ; LPF denotes low pass filter. ZOH( 4 ( ), ) is the zero-order holder to generate the same data dimensions with the reference position in robot motion replication, while the experiment in this paper is implemented under real time, so resampling is not necessary. The whole procedure and example of sEMG enveloping is shown in Figure 4. 5 ( ) is supposed to be the envelop of sEMG signals and it is easy to calculate envelop variation using differential method by (15): 2.4. Enveloping of sEMG Signals Based on Multiband AM-FM Model. According to [27], sEMG signals are viewed as monocomponent amplitude and frequency modulating (AM-FM) signals due to their multiple characteristic bands, which can be depicted as in (12), so here the key point is to extract smoothing ( ) features from sEMG signals: where ( ) and Θ( ) are amplitude and phase of sEMG signals, respectively. As mentioned above, sEMG signals are the summation of MUAP, which can be denoted in (13), so the amplitude of sEMG signals can be enveloped by the primary  AM-FM components though the phase of sEMG signals is not considered in this paper.
where ( ) denotes the instantaneous summation of MUAPs, ( ) denotes residual errors caused by disturbance, modelling error, finite summation, and so forth. In detail, envelop extraction of sEMG signals needs 6 steps as specified in what follows: (1) Preprocessing: raw sEMG signals are detected and preprocessed as the same rate mentioned above.
(2) Signal segmentation: to get a smoothed envelop with comparatively low time delay, sEMG signals are experimentally segmented into 200 ms segments without overlap.
(3) FIR filtering: a bank of 48th-order finite impulse response (FIR) filters is utilised to get clean sEMG signals. Functions implemented in MATLAB are as follows: window = hann (order + 1) , where order = 48, fir1 is designed as a bandpass filter, Energy Separation Algorithm (DESA) proposed in [45] which can be depicted in what follows: where 4 is the output from the 4th step of processing; is samples in one segment; ( ) is one sample at sample instant ; is the length of one segment.
(5) Sequence smoothing: 21-pint of median filter is empirically designed and employed to remove spikes in the IE sequence.
The overall procedure and example of sEMG amplitude envelop are shown in Figure 5.

Human Arm Stiffness Mapping to Robot Joint and Variable
Impedance Interface Design. From the stiffness variation estimation, we can obtain stiffness variation regulation of arm, which can be transferred from Cartesian space to joint space via the following: where ℎ ( ) = diag( ℎ ( )) denotes estimated joint stiffness at sample constant and ∈ 6×7 is human arm Jacobian. In this paper, ( ) is estimated Cartesian endpoint stiffness variation of operator's arm. In particular, for single joint stiffness estimation, Δ ℎ ( ) = 2 , where is the length of link connected to the specific joint. Robot arm joint stiffness is calculated by the following: where ( ) and 0 are the current joint stiffness at sample instant and initial joint stiffness of robot arm, respectively, Δ ℎ ( ) is the variation of stiffness at sample point , and is the sum of previous stiffness variations. Furthermore, it is of importance to normalize the stiffness within the specified stable region for stability during human robot stiffness mapping; otherwise it will make the mapping nonconsistent which may cost a lot of time to tune control parameters.
Here, robot arm stable region max and min can be obtained experimentally. Thus, we can employ the following: where ( ) denotes modified stiffness value based on estimation from sEMG, max is the maximum value of estimated  joint stiffness regulations, min is the minimum value of estimated joint stiffness regulations, max and min are the limit of stability of Baxter Robot, and is the estimated joint stiffness in each sample constant. To sum up, the mapping model and impedance interface design are shown in Figure 6 where the robot commanded motion will be the reference trajectory or a fixed point. The mapping stiffness will be the gain of position error in the feedback loop and the gain for velocity is the damping which will be simplified as (20), where is chosen empirically. So the overall experimental system can be seen as a human robot integrated system which is illustrated in Figure 6. The impedance interface is a PD control in each joint as described in what follows:

Experimental System Set-Up
The hardware set-up for stiffness variation coefficients estimation is shown in Figure 8. Human wrist is fixed by coupling mechanism with force sensor assembled on the base. Human subject elbow is adhered to a portable support to compensate arm gravity. MYO is worn on the upper arm near elbow to detect muscle (see Figure 9) activities.
The experimental robot platform in this work is Baxter Robot which consists of a torso, 2 DOF head, and dual 7 DOF arms, torque sensors, joint encoders, and direct programming access via a standard ROS (Robot Operating System) interface, and so forth. Each joint of the Baxter Robot arm is driven by a Serial Elastic Actuator (SEA), which provides passive compliance to minimize the force of any contact or impact. Baxter Robot shoulder joints 0 and 1 are employed in antidisturbance test, respectively. Overall experimental system is shown in Figure 7. Static and dynamic antidisturbance tests are implemented on Baxter Robot, respectively, as follows: (1) Antidisturbance control under posture maintenance: as shown in Figure 10, 1 joint would respond for external disturbance with variable stiffness adaptation while other joints will be kept in high stiffness mode without external disturbance effect. During the test implementation, the human subject ( 1 ) would give a random disturbance on the robot arm, while human subject ( 2 ) who wears MYO Armband on the forearm would strengthen his arm when seeing robot arm under disturbance by visual feedback.
(2) Antidisturbance control under path following: as shown in Figure 11, 0 joint is chosen for antidisturbance dynamically. Baxter Robot 0 joint would follow a reference trajectory from A to B, while human subject ( 1 ) would apply a random disturbance on Baxter arm when the arm is moving from A to B and the other subject ( 2 ) does similar antidisturbance task as the posture antidisturbance control test when robot is exposed to disturbance. Other robot joints are MYO Armband Arm support Coupling kept in high stiffness mode as well in order to clearly identify the antidisturbance performance of Baxter Robot single joint.

Human Arm Stiffness Variation Coefficients Estimation.
As indicated in Figure 12, the relationship between force and sEMG is estimated by linear model using LS method. It can be concluded that the force and sEMG have an approximate relationship according to both of two sEMG enveloping methods depicted in Section 2. Figure 13 shows the human arm stiffness coefficients estimation results using the two methods. It indicates that both of two methods have the similar properties for coefficients estimation. While fast linear enveloping is faster than AM-FM enveloping, but the later seems to be better in sEMG enveloping as shown in Figures 4 and 5. It can be seen that the enveloping spline in Figure 5 is smoother than Figure 4 C H 8 Figure 9: Muscles, number, and positions of electrodes involved in the elbow joint stiffness estimation [46].  Figure 10: Robot arm antidisturbance test I for posture maintenance (the dash yellow line) control: robot arm imitates human arm stiffness adaptive features in perturbation environment. Human arm will get cocontracted to compensate the external disturbance and this feature is transferred to robot arm by extraction of stiffness variations from sEMG signals. enveloping algorithm is employed for antidisturbance control under posture maintenance. Test results are shown in Figure 14. There are 3 steps to test human operated antidisturbance when doing robot posture control. Human subject 2 strengths his elbow when the robot is exposed to disturbance caused by subject 1 , and when the disturbance disappears, human subject arm 2 returns relaxed. It can be seen that human arm stiffness variation is sensitive and successfully eliminates external disturbance. As indicated from the disturbance force enveloping, when antidisturbance strategy is introduced during the test, amplitude of oscillation on Baxter joint 1 decreases dramatically even though the disturbance force increases.

4.3.
Antidisturbance Control under Path following Results. As shown in Figure 11, antidisturbance control is implemented  Figure 11: Robot arm antidisturbance test II: Baxter Robot 0 joint will follow a reference trajectory while a random disturbance imposes on Baxter Robot arm, and the human subject 2 would strengthen his arm to reduce the external disturbance with the same stiffness transfer interface as antidisturbance test I. under robot path following scenario, based on fast sEMG linear enveloping algorithm. Test results are shown in Figure 15. It can be seen that robot can follow the reference trajectory with no external disturbance, but when the external disturbance is employed on the robot arm, robot joint 0 fails to follow the reference path and even oscillation occurs in the test. However, when human operator transfers his stiffness to Baxter Robot joint 0 , Baxter Robot tends to be more stable and robust to the dynamic disturbance with a low position error.

Conclusion
This paper proposes a new discrete-time algorithm for stiffness extraction from surface electromyography (sEMG), the   performance of which has been tested under teleoperation antidisturbance experiment. We have used two methods of sEMG enveloping methods to extract discrete stiffness variations from sEMG, respectively. Experimental studies show that both methods are efficient and robust to external dynamic environment. The proposed method with these features could be potentially applied in various teleoperation and human robot interact tasks. Discrete Dynamics in Nature and Society

Competing Interests
The authors declare that they have no competing interests.  Step I, robot moves smoothly following a sine wave; Step II, rest for motion in disturbance environment; Step III, robot replicates the path regulation as Step I; Step IV, rest when finishing path following; Step V, human subject 2 does antidisturbance task when human subject 1 exerts disturbance on the Baxter robot 0 joint, which is used to follow a sinusoid reference trajectory.