Hammerstein–Wiener Multimodel Approach for Fast and Efficient Muscle Force Estimation from EMG Signals

This paper develops a novel approach to characterise muscle force from electromyography (EMG) signals, which are the electric activities generated by muscles. Based on the nonlinear Hammerstein–Wiener model, the first part of this study outlines the estimation of different sub-models to mimic diverse force profiles. The second part fixes the appropriate sub-models of a multimodel library and computes the contribution of sub-models to estimate the desired force. Based on a pre-existing dataset, the obtained results show the effectiveness of the proposed approach to estimate muscle force from EMG signals with reasonable accuracy. The coefficient of determination ranges from 0.6568 to 0.9754 using the proposed method compared with a range of 0.5060 to 0.9329 using an artificial neural network (ANN), generating significantly different accuracy (p < 0.03). Results imply that the use of multimodel approach can improve the accuracy in proportional control of prostheses.


Introduction
Since 1952, the relationship between electromyography (EMG) and muscle force has been the focus of research for many applications ranging from prostheses control to active user-driven exoskeletons. Different approaches based on empirical and mechanical models are proposed for estimating muscle forces from the EMG signals. As the name implies, mechanical models are developed from physical laws (mechanics, biologics, electric, etc.). They are to be helpful in applications where individual muscle kinetics and kinematics are of interest. Hill's muscle model is considered the most used physical model in musculoskeletal application to study phenomena in which only mechanical behaviour is considered [1][2][3][4]. However, this model has not shown its complete suitability in the control of upper limb myoelectric prosthesis.
Within our scope for controlling upper limb myoelectric prostheses, accurate estimation of the overall force from the EMG drive proportional control-based schemes. Therefore, many black-box models and statistical identification techniques are proposed [5][6][7][8][9][10][11]. The degree of linearity between the EMG and muscle force was debated for many years. Staudenmann et al. (2010) showed in a review paper that the relationship between EMG muscle activities and muscle force is not necessarily linear from a physiological and biophysical perspective. It depends on the level of recruitment and, therefore, on the composition of the type of muscle fibres [11][12][13][14][15][16][17]. Indeed, muscle force is modulated by the number of motor units recruited and their activation frequency [18,19], explaining the nonlinear relationship between EMG and muscle force [20,21]. In [22], experiments based on a recruitment range controlled by an electrical stimulation protocol showed that the EMG-muscle force relationship could be assimilated as a linear relation for muscles constituted with one type of fibre. However, mixed-fibre-type muscles, composed of fast and slow switching fibres, exhibit a nonlinear relationship between muscle activities and forces. Furthermore, Kamavuako and Rosenvang (2012) provided evidence of hysteresis in the EMG-force relationship to indicate non-linearity [23].
From a mathematical modelling point of view, several estimation methods [24][25][26][27][28][29][30] were proposed in the literature, including linear and nonlinear approaches with different outcomes. In this perspective, Hahne et al. proposed a comparative study of linear and nonlinear regression techniques to estimate muscle force from the EMG signals for myoelectric control movements. This study was based on four estimation algorithms: linear regression (LR), mixture of linear experts (ME), multilayer perceptron (MLPs), and Kernel ridge regression (KRR). They showed that, on the whole, the interest of nonlinear methods, especially the KRR method, is considered a nonparametric and nonlinear statistical learning method [27]. Luo et al. (2019) [30] proposed a three domains fuzzy wavelet neural network algorithm without prior knowledge of the biomechanical model to analyse the force muscles. This approach is based on the mean absolute value of the electromyographic signal, used as the input for the potential model. Recently, Wimalasena et al. adapted a new AutoLFADS approach based on unsupervised deep learning, more precisely, recurrent neural networks (RNNs). This approach is applied to estimate muscle activities from multi EMG signals [31].
The superiority of nonlinear techniques supports the non-linearity of the EMG-force relationship, especially at higher contraction levels when many motor units are recruited. Non-conventional black-box models, such as artificial neural network-based methods, machine learning, or deep learning approaches, are suggested for modelling complex systems. These approaches are proven effective in predicting force, especially for nonlinear systems; nevertheless, they are considered complex techniques in modelling architecture and training methods [32][33][34]. Despite this, many non-conventional techniques only show good performance in specific tasks due to "Catastrophic Forgetting", which is the problem of losing information on a first task, T1, after training a second task, T2. Otherwise, the performance of T1 will significantly decline [35][36][37][38]. Additionally, non-conventional techniques are costly for production use. Indeed, complex problems require an extensive network that can be exceptionally time-consuming to compute at inference times [5][6][7][8]. To address the challenges encountered in modelling the non-linearity of the EMG-force relationship for upper limb prostheses, and thus to overcome the disadvantages of the techniques above, in this study, we propose a multimodel approach. A multimodel approach is suggested for nonlinear system modelling to decompose a whole operation area of a studied process into a defined number of sub-operating regions. In each one, a local model is computed. This approach is successfully applied even when dynamical system characteristics vary considerably over the operating regime. Thus, the multimodel concept is considered an exciting method to improve the model's performance in terms of precision and without overly increasing the computational burden or the number of parameters to adjust [39][40][41][42][43][44].

Related Studies
We propose a novel MIMO multimodel structure to model muscle force from the EMG signals in the present study. Sub-models of the proposed design are based on the nonlinear Hammerstein-Wiener (H-W) model, which was applied in different fields to estimate complex and nonlinear systems [45][46][47][48][49][50][51][52][53][54][55]. Indeed, the H-W model was used in some works to estimate muscle force from EMG signals. In this context, Kumar et al. (2011) proposed a nonlinear modelling approach to characterise the relationship between EMG and muscle force in [47]. Their study developed a hybrid model based on the fusion technique between multi nonlinear autoregressive exogenous (ARX) and Hammerstein-Wiener (H-W) models designed for different force/EMG data. Each force sensor model's output is first fused using an adaptive probability algorithm. Secondly, the outputs of various sensors are combined by the same adaptive algorithm to estimate the desired force.
Despite the acceptable results, the proposed approach is significantly complex because of the high number of nonlinear models and parameters to be calculated for each hybrid model developed for each force sensor.
In 2011, a H-W based model was proposed to estimate one muscle force from one EMG Signal. The proposed approach shows that despite the acceptable error rate shown in the direct validation of the H-W model, a significant error is obtained in cross-validation. This approach delivers good performance only for particular data [48]. Sebastian et al. (2011) presented an analysis of different filtering methods to measure the dynamics of surface EMG signals. This analysis was compared based on the EMG finger force model based on H-W estimation to show that, in this case, the nonlinear spatial filters gave better-fit values [50]. Indeed, regardless of the application, approaches based on the H-W model show significant performance to mimic nonlinear systems [48][49][50][51][52][53][54][55]. However, this technique requires parameter adjustment for each new dataset. Furthermore, the real and the estimated force error is significant, especially in cross-validation [39][40][41]. This uncertainty can affect the quality of monitoring or control. Additionally, most of the proposed models characterise just one force profile/channel from EMG signals. However, it is most interesting for complex movements to simultaneously model two or more muscle forces. In this work, based on a multimodel approach, we propose a MIMO muscle force estimator from EMG signals of the upper limb using diverse force profiles. The novelty of the proposed approach lies in (1) its efficiency in estimating force with a minimal number of parameters; (2) ability to optimise the number of used EMG signals to estimate multiple muscle force profiles; (3) stable performance with direct and cross-validation with a minimal computational burden. To show the efficiency of the proposed modelling structure, we have compared it with the artificial neural network (ANN) modelling technique known for its effectiveness, especially for data sets having nonlinear relations [27,56].
The paper is organised such that, in Section 3, we describe the proposed multimodel structure of the proposed experimental approach and define a sub-section for data analysis. The validation of the developed method and a comparative study with the ANN model are presented and interpreted in Section 4. Section 5 discusses and interprets the results from different scenarios that we developed to validate the multimodel approach. Finally, the conclusion and perspective are presented in Section 5.

Methods
The multimodel technique consists of replacing a unique representation of a complex dynamic system with a nonlinear regime by combining multiple models, called sub-models. Different measurement pairs developed each sub-model to describe the behaviour of the studied system in its whole operating domain at a particular operating point. Figure 1 show the principle of system modelling using a multimodel approach. Sub-models are grouped into a model base, named library. We note that they can have different structures or/and orders; however, no one can describe the studied nonlinear process in its whole operating area. Based on the interaction of the sub-models, the decision unit allows us to estimate the contribution of each sub-model in a weighted manner according to their posterior likelihoods to foster the most pertinent one at each time. The decision block controls the output block to compute the output of the multimodel structure, Y mm, which is obtained by the fusion or the switching of the different sub-model responses.
To explain and simplify the basic principle of the proposed MIMO multimodel structure, we suggest, for example, fixing the number of inputs and outputs to 2, i.e., we have to estimate two muscle forces, F m1 and F m2 , from two EMGs signals, EMGm1 and EMGm2.
The proposed multimodel structure is described by the following steps (the different steps will be detailed in the next section), as shown in Figure 2: To explain and simplify the basic principle of the proposed MIMO multimodel structure, we suggest, for example, fixing the number of inputs and outputs to 2, i.e., we have to estimate two muscle forces, Fm1 and Fm2, from two EMGs signals, EMGm1 and EMGm2.
The proposed multimodel structure is described by the following steps (the different steps will be detailed in the next section), as shown in Figure 2:  Inputs of the multimodel bi-forces estimator. F m1 and F m2 : Outputs of the multimodel bi-forces estimator.
Output of sub-model (i). F i1 and F i2 : Output of sub-model (i) obtained by applying multimodel inputs, EMGm1 and EMGm2, respectively. err i1 and err i2 : Errors of sub-model (i) computed between its real output, F i , and outputs, F i1 and F i2 µ i1 and µ i2 : Validities of sub-model (i) according to F m1 and F m2 , respectively.

•
Step 1_Elaboration of sub-models: Definition of sub-models of the library: Submodel (i): model defined for the measurements couple (EMGi, F i ).

•
Step 2_Normalised residues estimation: Computation of the normalised error of each sub-model (i): err' i1 and err' i2 .

•
Step 3_Validity Computing: Computation of the weight, also named validity, of each sub-model: µ i1 and µ i2.

•
Step 4_Outputs Computing: Computation of outputs F m1 and F m2 .
Below is the explanation of the four steps.
Step 1_Elaboration of sub-models: In terms of duration, form, and even complexity, the considered experimental data present different force profiles. Furthermore, the complex and unpredictable characteristics of EMG signals and the anatomical and physiological properties of muscles justify the choice of the Hammerstein-Wiener (H-W) structure to estimate muscle force. Indeed, this structure considers inputs and outputs as nonlinear information [51][52][53][54][55]. It allows observing input/output measurements to develop a model that approximates the behaviour of the underlying system. Furthermore, a nonlinear H-W model, considered a black-box method, can be presented as a concatenation mapping between previously observed data and a regression space followed by a nonlinear function. In addition, in our case of study, the data analysis stage, detailed in section B, shows the non-linearity of the inputs and outputs, promoting the use of the H-W regressor and justifies our choice to utilise it for this study.
The following equations can formulate the system identification approach: The following general dynamic model can also represent the Output vector: where: v(k): Noise signal, presenting the prediction of error between actual and estimated output.
In the case of an ideal model and good prediction of F, v(kt) approaches zero, and the mapping function is as follows: with: : Observation vector, also named regression vector, contains inputs and outputs data of : Parameter vector θ = [θ 1 , θ 2 , . . . , θ p ] p is the number of parameters to be estimated.
The parameters of the HW models are calculated during the training phase based on the parametric identification algorithm recursive least squares algorithm (RLS). The evolution of these parameters shows that they converge towards a constant value, making it possible to consider fixed parameters that do not need readjustment. The following equations can describe the RLS algorithm: where: ∧ θ (k): vector of estimated parameters, P (k): adaptation matrix, ε (k): estimated error, y (k): actual output of the system to identify, k: is the discrete-time. The sub-model given above was implemented, and the result is summarised in Figure 3 for a representative subject. We can notice a good correlation between the real muscle force and the response of the model. Parameters estimation of the H-W sub-models is based on the recursive least squares algorithm [57,58]. We note that the actual force is a continuous blue line and the estimated force is the dotted red line. We note that the response of the developed H-W sub-models is similar to the experimental data. However, these models are limited if we change input/output data by applying experimental data from which one sub-model was developed to another sub-model referring to the same type of force profile. In this sense, the limit of H-W models will be overcome with the multimodel approach. As shown in Figure 4, Hammerstein-Wiener can be presented by a serial association of nonlinear, linear, and other nonlinear systems. Step 2_Normalised residues estimation: In the literature, different approaches were proposed to estimate validity degrees of sub-models [41,59,60] that can be defined by a switching function or by fusion of sub-model responses. In our study, we use a fusion technique based on computing normalised residues, err' i1 and err' i2 , of each sub-model to evaluate its pertinence to describe the system in an operating area at each time, estimation of errors between the real and the estimated outputs: Step 3_Validity Computing: Validity coefficients must belong to the interval [0,1] and vary contrary to residual values [42][43][44][45]. Otherwise, the most pertinent model must have the minor residue and the most important validity value.
µ 1i and µ 2i are designed validities of sub-model (i) according to F m1 and F m2 , respectively.
Step 4_Outputs Computing: The multimodel outputs F m1 and F m2 are computed as follows:

Experimental Approach
The proposed study makes use of pre-existing data [56]. In brief, 10 able-bodied human subjects (7 w/3 m; age range, 21-37 years, mean 26.9 years) took part in this experiment following the Declaration of Helsinki (approval no.: N-20080045). Each subject received an information sheet about the investigation and gave written consent before participation. Subjects were asked to follow six force profiles randomly assigned. The six force profiles were defined as follows: (a) a step increase in static grip force with five increments of a 10 s duration of 10N (step), (b) two linear ramps (saw), (c) a freely varying force (vol), and (d) a steady-state force at 50N (single-level); a slow increase and decrease for a duration of 10 s (circle). Four trials were recorded for each force profile, and a rest of 1 min followed each practice. In total, 160 trials were collected. The force produced during power grip was measured using a commercially available handgrip dynamometer (Vernier Software & Technology, accuracy ± 0.6 N, operational range 0-600 N, grip size 50 mm × 25 mm). Surface EMG was measured in a bipolar configuration using disposable Ag/AgCl surface electrodes (Ambu Neuroline 720, Denmark) from the M. extensor carpi radialis (ECR), flexor digitorum superficialis (FDS), and flexor carpi radialis (FCR). The sEMG signals were amplified by a factor 2000 with a multichannel surface EMG amplifier (EMG16, LISiN, Torino, Italy), and band-pass filtered (20-500 Hz) and sampled at 1000 Hz. Figure 5 present the experimental setup to measure the force profile and EMG signals from the forearm muscle.

Data Analysis
A data processing step was applied to digitally rectify the force and EMG signals before the estimation stage. We have recovered the envelope of the signals with a low pass filter using a 6th order Butterworth filter to cut off frequency at 1Hz and finally normalise the data accordingly using Min-Max normalisation in the data preprocessing stage. The Min-Max normalisation is defined as follows: with: d n : the normalised value, d r : the real value, d min : the minimum value, d max : the maximum value. Then, data were sent to the proposed approach, as described above, for force estimation and compared with performance based on the artificial neural network as previously suggested. The appropriate number of sub-models depends on the different operating domains that can describe the whole behaviour of the studied nonlinear process. It is fixed according to statistical analysis of the EMG data. Indeed, the multimode approach is generally proposed to optimise the characterisation of a nonlinear system having different operating points. A sub-model represents each operating point. However, in this study, we have input/output type models (black-box models) possessing the same structure but different non-physical parameters that change according to the change of the data, especially the EMG inputs, considered unpredictable signals and challenging to characterise.
Furthermore, EMG is asynchronous and can be assimilated to a succession of irregular and unpredictable wave trains. Therefore, we analysed the recorded data of EMG signals relating to each type of movement to study and determine the number of data classes in terms of distribution and statistical characteristics. The number of sub-models N is then relative to these classes. Otherwise, each sub-model has a well-determined distribution class of EMG inputs. Figure 6 represent the box plots and the distributions of the EMG signals, measured for the same movement "Two linear ramps (saw)." We remark that we have mainly three distinct behaviours of these data, which allows us to propose three sub-models; each presents a different statistical characteristic. • Scenario2: Predict profiles of arbitrary forces from sub-models characterising 02 known forces. For example, using step and circle to estimate free profiles (vol).
• Scenario3: Predict profiles of arbitrary forces from sub-models characterising 03 known forces. This is the same as scenario 2, but with three standard profiles. For the different scenarios, we note that inputs/outputs of the multimodel approach are considered unknown, and any sub-model of the library does not represent them. Three performance measures: coefficient of determination (R2), root mean squared error (RMSE), and computational time (CT), were used to assess the proposed approach. For each performance metric, a two-way repeated measure analysis of variance (ANOVA, IBM SPSS Statistics 26) with factors methods (multimodel vs. ANN) and scenarios (1,2 and 3) was used to assess the performance of the proposed approach. Table 1 present the computational environment used in this paper.  Based on this consideration and to show the ability of the proposed multimodel approach to estimate different kinds of force profiles (FP), we present the following scenarios: • Scenario1: Predict force profiles from sub-models characterising the same kind of desired forces. For example, predicting a step from another instance of the step profile. • Scenario2: Predict profiles of arbitrary forces from sub-models characterising 02 known forces. For example, using step and circle to estimate free profiles (vol). • Scenario3: Predict profiles of arbitrary forces from sub-models characterising 03 known forces. This is the same as scenario 2, but with three standard profiles. For the different scenarios, we note that inputs/outputs of the multimodel approach are considered unknown, and any sub-model of the library does not represent them. Three performance measures: coefficient of determination (R2), root mean squared error (RMSE), and computational time (CT), were used to assess the proposed approach. For each performance metric, a two-way repeated measure analysis of variance (ANOVA, IBM SPSS Statistics 26) with factors methods (multimodel vs. ANN) and scenarios (1, 2 and 3) was used to assess the performance of the proposed approach. Table 1 present the computational environment used in this paper. The architecture of the considered ANN algorithm can be described as follows, Figure 7:  Figure 8 show the ability of the developed approach to predict different kinds of force profiles from surface EMG signals. This figure illustrates the simulation results of the proposed scenarios to show the interest of the proposed approach compared with the artificial neural network in the prediction of forces. CT was 110.1 ± 13.4 ms and the multimodel approach was significantly lower (p < 0.01) than the ANN (825.0 ± 136.0 ms). Indeed, the fast training of the proposed method allows an online implementation to deal with continuous data preprocessing and generate more efficient real-time estimation. The overall results of the scenarios are summarised in Table 2 for R 2 and RMSE, respectively. For R 2 , the proposed approach (0.876 ± 0.045) outperformed (p < 0.001) the ANN (0.738 ± 0.066). There was a significant difference between scenarios, with scenario 3 performing worst (p = 0.24). Furthermore, Table 1 show a significant interaction (p = 0.028) between approaches and scenarios. The overall results of the scenarios are summarised in Tables 1 and 2 for R 2 and RMSE, respectively. For R 2 , the proposed approach (0.876 ± 0.045) outperformed (p < 0.001) the ANN (0.738 ± 0.066). There was a significant difference between scenarios, with scenario 3 performing worst (p = 0.24). Furthermore, Table 2 showed a significant interaction (p = 0.028) between approaches and scenarios. The RMSE of the proposed approach was 0.026 ± 0.006, significantly (p < 0.001) lower than ANN (0.039 ± 0.009). RMSE was different (p = 0.021) between scenarios but there was no interaction (p = 0.068), Table 3.   Table 4 and Figure 9 show a comparative study between the multimode approach and ANN regarding computation time measured for the different subjects according to the proposed three scenarios. This study shows the interest in the multimode approach compared to the ANN. Indeed, the computation time of the mm approach is 10 to 20 times less than that of the ANN method. It is also almost constant for the different scenarios, unlike the ANN approach, where the computation time depends on the scenario. In addition to the computational time, the memory requirements are also important to consider. In this perspective, Table 5 present the memory requirements for both methods using the computational environment described in Table 1. Indeed, the ANN approach requires almost 25% more memory than the multimodel approach.

Discussion
This study proposed a multimodel approach based on nonlinear Hammerstein-Wiener to estimate force from EMG with the MIMO feature that enables multiple force profiles to be learned simultaneously as sub-models. The concept of force estimation from EMG is essential for the proportional control of myoelectric prostheses [2][3][4]. Overall, the proposed approach outperformed the artificial neural network for all metrics used in the study. The architecture of the ANN used in the study is the same that was used by [47] to allow direct comparison. We have deliberately chosen not to compare with other H-W models because our target application is myoelectric control.
Furthermore, ANNs are established as the golden standard for force estimation due to the poor performance of linear regression, especially for simultaneous force estimation [30,47]. By comparing a previous architecture using the same dataset, the actual performance of the proposed approach can be demonstrated. Regardless of the method used for the estimation, we see that they perform very well for simple problems (scenario 1) with a decrease in performance when the problem becomes more complex. Nevertheless, the proposed approach is more robust to changes in complexity as proven using both R 2 and RMSE. Furthermore, the proposed method requires less computational power and time compared with the ANN. Particular attention can be brought to scenario three, where it is shown that the increase in the number of sub-models does not automatically generate an improvement in force estimation. Indeed, the choice of the type of sub-models that constitute the library of the multimodel structure is critical.
Nevertheless, the multimodel approach allows us to not only characterise arbitrary forces from EMG signals but also improve the prediction's statistical performance, the linear resemblance between the real and the estimated force, and provide multi forces from a simple, well-defined model sub-models library. Despite the improved performance, the dataset used in the study deals with sequential force profiles, while the challenging problem with myoelectric prostheses is the ability to control motions simultaneously. This study intended to pave the way for using the approach as we are designing a novel experimental protocol that will include both offline and real-time estimation of simultaneous force profiles. However, the situation with the COVID-19 pandemic is impeding the recruitment of subjects, especially participants with limb amputations.

Conclusions
This paper aims to characterise muscle forces from electromyography signals of the upper limb. A multimodel approach was developed to extend and increase the flexibility to model several kinds of force profiles to overcome this challenge. The proposed method is mainly based on nonlinear Hammerstein-Wiener models used to reconstruct the library of the multimodel structure. The validity of each H-W sub-model is computed to estimate the contribution of each one to estimate the desired outputs. The multimodel approach was validated for different scenarios to propose libraries with varying numbers of sub-models representing the same or other types of the desired force profile. The proposed method shows acceptable metrics values in R 2 , RMSE, and computed total CPU time. However, it would be of great interest to continue the real-time experimental investigations to assess the method, especially for amputee subjects. In addition, the results of the multimodel approach are promising in terms of efficiency, computation time, and memory requirements.