Motion Artifact Suppression for Insulated EMG to Control Myoelectric Prostheses

Myoelectric prostheses help amputees to regain independence and a higher quality of life. These prostheses are controlled by electromyography, which measures an electrical signal at the skin surface during muscle contractions. In this contribution, the electromyography is measured with innovative flexible insulated sensors, which separate the skin and the sensor area by a dielectric layer. Electromyography sensors, and biosignal sensors in general, are striving for higher robustness against motion artifacts, which are a major obstacle in real-world environment. The motion artifact suppression algorithms presented in this article, prevent the activation of the prosthesis drive during artifacts, thereby achieving a substantial performance boost. These algorithms classify the signal into muscle contractions and artifacts. Therefore, new time domain features, such as Mean Crossing Rate are introduced and well-established time domain features (e.g., Zero-Crossing Rate, Slope Sign Change) are modified and implemented. Various artificial intelligence models, which require low calculation resources for an application in a wearable device, were investigated. These models are neural networks, recurrent neural networks, decision trees and logistic regressions. Although these models are designed for a low-power real-time embedded system, high accuracies in discriminating artifacts to contractions of up to 99.9% are achieved. The models were implemented and trained for fast response leading to a high performance in real-world environment. For highest accuracies, recurrent neural networks are suggested and for minimum runtime (0.99–1.15 μs), decision trees are preferred.


Introduction
Electromyography (EMG) sensors are applied, beyond others, for myoelectric prosthesis control, diagnostic purposes and exoskeletons. The state-of-the-art dry sensors, which require a conductive connection to the skin, are typically applied for myoelectric prosthesis control. This research group has already published the developed flexible insulated EMG sensors [1,2] to overcome the disadvantages associated with dry conductive EMG sensors. The insulated sensors avoid pressure marks and are independent of sweat and the skin condition in general.
Various research has presented the control of high-level dexterity prostheses [3][4][5]. Well-established EMG features have been implemented [6][7][8] and different signal processing algorithms have been developed [9][10][11]. Although these algorithms achieve high accuracies for distinguishing hand movements, they often lack robustness in real-world environment [12,13]. However, in terms of practical application, robustness is preferred over technologically complex and unreliable systems [14].
The main reason for the lacking robustness are the motion artifacts [15], which have various origins. They typically occur in low frequency range. However, these artifacts also appear in the frequency range of the contraction EMG, which limits the possibilities for filtering. To achieve high robustness, real-time motion artifact suppression algorithms were implemented on an ultra-low-power Figure 1 shows the signal flow graph of the capacitive EMG signal. The data acquisition, filtering, rectification and smoothing were described in detail by Roland et al. [1,2]. The block diagram of the digital signal processing is plotted in Figure 2. Additional data pre-processing, the feature calculation, the decision algorithm (logistic regressions, decision trees, NN or recurrent neural networks (RNN)) and the data post-processing as well as the signal delay are described in this section.  Figure 1. First, the EMG signal coupled from the human body is amplified and filtered by the analog circuit (ASP). In the µC, the signal is then digitally filtered and evaluated (DSP). The contraction EMG signal is provided at the prosthesis drive via the digital-to-analog converter (DAC). For experimental applications, the signal can be connected to the oscilloscope or Bluetooth (BLE) module (Adapted from Roland et al. [2]).

Data Acquisition and Pre-Processing
Insulated EMG data were acquired with the sensor developed by Roland et al. [1,2]. Figure 3 shows flexible coupling electrodes, which are an assembly of textiles or a flexprint. The EMG sensing system filtered the signal with a first-order analog bandpass ( f CL = 11 Hz, f CU = 1064 Hz), a digital comb (50 Hz and harmonics) and a first-order digital lowpass ( f C = 531 Hz). In the feature calculation path, a second order digital high pass ( f C = 60 Hz) was implemented. The optimal choice of these above-mentioned parameters was investigated in Roland et al. [1,2]. The parameter values were selected according to these findings. The flexible sensor was placed at the musculus extensor digitorum at the human forearm of one subject. Four different types of data were measured. Data for strong and weak isometric contractions, random sequences of short contractions and measurements of many artifacts were acquired. Therefore, various different motion artifacts were generated, such as sensor lift-off, external mechanical shocks or vibrations.
The data were measured with a digital oscilloscope (HS3 of TiePie Engineering) at the µC's digital-to-analog converter. The 10 s measurements were sampled by the oscilloscope at 10 kHz to fulfill the Nyquist-Shannon sampling theorem [17]. A measured contraction and a motion artifact are presented in Figure 4. The measured EMG data were downsampled by averaging to 2 kHz. According to Roland et al. [2] downsampling to 2 kHz does not reduce the EMG signal quality. This sampling frequency of 2 kHz was applied for the feature calculation. Time windows with small signals, which would not lead to a prosthesis drive activation, were removed from the data set.  After these pre-processing steps, a total of 381 s of contractions and 489 s of artifacts were the remaining training, test and validation data. These data comprised a wide variety of different contractions and artifacts to achieve high robustness in real-world environment.
The contraction and artifact signals were mixed to 1 s long windows. The contraction and artifact signal windows were alternated, hence a fast reaction to the changing signal class was ensured. The features were calculated over these alternating time windows. Features which were slow at the transients led to errors at the transition from one class to the other. By iterating the classes in the data set, these slow features were consequentially removed by the feature selection. This is essential for real-world performance, where artifacts and contractions are alternating and a fast reaction of the decision is required. For the training and validation set of the RNN, the duration of the windows was varied for each snipped, thereby preventing the RNN to learn the periodic behavior of the training set.
The training (70%) and validation (15%) set amounted to 85% of the data. The validation data were selected out of the 85% of the data when training the models with cross validation, see Section 2.6. The remaining 15% were the test set, which was equal for all models. The training, validation and test set comprised an equal percentage of the strong, short, and weak contractions, and artifacts according to their share on the whole data set. A sample of the features highly correlates to the previous and next sample in time, the adjacent feature values are very similar. Drawing individual samples for the allocation to the training, validation and test data would lead to non-generalizable results. Therefore, continuous time sequences were selected in the allocation in order that the training, validation and test data are independent from each other.

Feature Calculation
In this contribution, new features and modifications of features for biosignal processing, and well-established features were included to the feature calculation. The TD features were hand-crafted as the models should achieve high accuracies with a shallow model architecture. The feature vector f eat was calculated by filtering with a first-order lowpass, an exponentially decaying moving average filter (EMA) [18]: where f (x i ) denotes a function of the current comb and lowpass filtered value x i . The function f (x i ) is defined differently for each feature. The coefficient b was selected in the interval [0, 1] (in floating-point format). This b was implemented as a fixed-point with a multiplication and a shift operation to avoid operations with floating-point variables. This implementation requires less calculation resources than a standard moving average filter [2]. The features were designed with the aim to detect differences of contraction EMG and artifacts. These TD features also attempt to detect frequency information, which is relevant for the distinction of the classes, e.g., by the zero-crossing rate. Some features were implemented several times with varied parameters, which enables the extraction of different signal characteristics.
The features were designed with MATLAB R and implemented on the µC in C. The MATLAB R and C implementation were in the same value range and truncation, overflow and saturation were considered. Hence, the values listed in the following section are valid for both implementations. Some features were saturated at an upper bound (ub) and/or a lower bound (lb) to avoid a drift of the feature values, thus achieving a fast transition between the different classes. The selected parameters of the implemented features are listed in Tables 1-7.

Zero-Crossing Rate (ZCR)
The ZCR (Figure 5a) is a well-established feature for EMG signal analysis, which comprises frequency information. When the EMG signal crosses zero, the feature is increased by 100. This feature update parameter is set to 100 to prevent the feature variable from vanishing right away due to truncation caused by the shift operation in the lowpass filter in Equation (1). Moreover, the value was not selected higher to prevent overflows. The parameter was determined, so that no overflows occurred in the training data. These facts are also valid for the selection of this feature update parameter for the following described features. As also noise with low amplitude contributes to the feature, a hysteresis (hyst) was added to detect the crossings only, when the signal exceeds a certain amplitude. An upper bound ub was included in ZCR2S to avoid a drift of the feature value, thereby achieving fast transitions. The conditions for the ZCR were defined as where s denotes the sign of the previous zero-crossing.  The contraction EMG is usually superimposed to a slight artifact at short contractions, due to the relative movement of the muscle tissue to the skin surface at the beginning and the end of a contraction. To detect this superimposed EMG signal, the crossing of the signal through the smoothed signal, as shown in Figure 5b, increases the feature by 100.

Raw Signal Smoothed Hysteresis
A hysteresis hyst was included to avoid noise contributing to the feature. A feature variant with a saturation at an upper bound (ub) and a lower bound (lb) was implemented. The conditions for the MCR were defined as where s denotes the sign of the previous mean crossing and v is defined as The smoothed signal smo i is calculated with an EMA with the coefficient c: The input x i is delayed and scaled by When the EMG signal slope is changing (Figure 5c), the feature is increased by 100. To avoid noise contribution, the feature is calculated with the smoothed signal from Equation (5). Variants of the feature require the signal slope to maintain the direction for a minimum number of data points d min and a maximum number of data points d max . This way capturing the desired frequency information. Again, features with a saturation at a lower lb and an upper bound ub were implemented. The conditions were defined as where s denotes the sign of the previous slope. The variables w and p are counters and u and v are the smoothed signal from Equation (5):  The function f (x i ) for the MAV was defined as The MAV2 is calculated by means of the MAV1 by The feature MAV1S is saturated at the upper bound ub1 and the MAV2S at ub2. When the EMG signal exceeds a defined threshold thresh, the feature is increased by 100. The conditions for WAM1 were defined as For WAM2, the smoothed signal smo i is calculated with Equation (5) and the conditions are The VAR is the squared EMG signal, which was also implemented with a saturation at the upper bound ub. The VAR is calculated by

Downsampling
For the training of the NN, decision tree, and logistic regression, the features were downsampled from 2 kHz to 40 Hz to accelerate model training. In the implementation the features were not downsampled to achieve short time periods between the decisions. For the RNN, the features were downsampled from 2 kHz to 250 Hz, which led to a higher calculation effort in the learning process. However, an equal sampling frequency in the training and implementation is indispensable for high accuracies of RNNs.

Correlation of Features
To reduce the number of features before training, the Pearson correlation coefficients [19] of the feature vectors were calculated by pairwise comparison. For an absolute value of the correlation coefficient greater than 0.9, one feature of the pair was eliminated.

Normalization
For the logistic regression, NN and RNN (but not for the decision tree) the min-max normalization was calculated for each entry of the feature vector f eat: with y max = 1 and y min = −1 in floating-point format. The minimum and maximum ( f eat min , f eat max ) were determined for the feature vector calculated with the training data.

Training of Models
All models were trained with a 5-fold cross validation. The validation data were randomly drawn time sequences of the strong, short and weak contractions and the artifacts. The decisions within 100 ms after a transition of the target value were not included in the error calculation of the validation data. Wrong decisions, starting from 100 ms after a transition of the target value, were recorded as errors, thus achieving fast models. For the test data, this tolerance was set to 150 ms (cf. signal delay, Section 2.8.3). The models were trained with functions provided by the MATLAB R toolbox "Statistics and Machine Learning".

Logistic Regression
The features were selected with the lasso algorithm [20] assuming a binomial distribution. 30 different regularization penalties (λ) were applied. The range of the λ was defined such that the resulting feature vectors were ranging from a maximum number of regression coefficients to maximum sparsity. A vector with the desired number of features was selected and the non-zero entries were the selected features.
With these features the logistic regression model was trained by the Ridge Regression algorithm [21]. Finally, the model with the highest accuracy was selected.

Decision Tree
Binary decision trees were trained with different numbers of maximum splits, which avoids overfitting. The split and the respective feature was selected by calculating the Gini Diversity Index [22].

Neural Network (NN)
The features for the shallow NNs with one hidden layer were chosen by sequential backwards selection [23]. The activation function for the hidden unit was SATLINS, and PURELIN for the output layer. These activation functions allow a fast calculation on an embedded system. The NNs were trained by the Levenberg-Marquardt algorithm [24,25], aiming at the minimization of the cross-entropy loss [26]. The model was trained with a maximum number of epochs of 1000 and a maximum of 6 validation failures.

Recurrent Neural Network (RNN)
For the RNN (Figure 6), the maximum number of epochs was set to 100. All other hyperparameters were set as described in Section 2.6.3.

Offline Data Post-Processing
The model decision was further processed to increase the stability against short wrong decisions. The system should be prevented from incorrect reactions due to individual false classified samples. To allow a transition of the decision, a minimum number of equal consecutive decisions was required. The transition at the output was only conducted after (excl.) n slope equal decisions. For the offline post-processing n slope was set to 2 for the NNs, decision trees and logistic regressions, and to 3 for the RNNs. This parameter introduces an additional delay, which is depending on the model decision frequency and the value of n slope . The parameter n slope was selected in a way to achieve similar delays for the different implementations. Despite the additional delay, n slope is essential for a high stability of the system. This parameter value was selected, as it is a good trade-off between delay and stability, however, it can be adapted to the individual amputee's preference.

Implementation
The features and the models were implemented on an ultra-low-power µC, the ATSAML21E18B from Microchip Technology Inc. (Chandler, AZ, USA) with a 32 bit ARM R CORTEX R -M0+ processor, 256 kB flash, and 32 kB SRAM main memory. The clock frequency was set to 48 MHz to test the implementations. The clock frequency can be reduced (depending on the runtime of the selected model) to decrease power consumption [2].
The analog digital sampling of the EMG signal has been described by Roland et al. [2]. The digital signal processing sampling frequency was set to 2 kHz. For the RNN (250 Hz), the 48 MHz clock provides 192,000 clock cycles (4 ms) between two decisions. The NN, decision tree, logistic regression, and the feature calculation (2 kHz) provide 24,000 clock cycles (500 µs) between two samples.

Fixed-Point Representation
All algorithms were implemented in fixed-point format, as described by ARM Ltd. (Cambridge, UK) [27], which is faster than a floating-point implementation on this µC. For high precision and to avoid overflows or saturations the optimal value range of the variables was selected to exploit the full operating range of the data type. The quantization effects and possible overflows or saturations had already been considered in the offline design. The high performance of the µC 32 bit hardware multiplier was exploited with this fixed-point implementation. All divisions were replaced by shift operations to require low calculation resources.

Online Data Post-Processing
As mentioned above, the time between two decisions was deviating in the offline and online implementations for the NN, decision tree and logistic regression. Therefore, different post-processing parameters were selected. The parameter n slope was set to 20 in the implementation of the NN, decision tree and logistic regression. For the RNN, n slope was set to 3, equal to the offline design.

Signal Delay
The filtered EMG signal and the decision were multiplied at the output; therefore, they should correspond to each other in time. Hence, the input signal was delayed, so that it corresponds to the current decision. The smoothing had a time constant T of 51 ms [2], which introduced a delay. An additional delay was set to 100 ms in the implementation. Table 8 lists the accuracies of the resulting models. Please note that these accuracies were calculated for test data with many transitions. Additionally, many artifacts were in a similar frequency and amplitude range than the contractions. Furthermore, transitions longer than 150 ms were recorded as errors and small signals were removed from the data set. Consequently, these challenging data led to higher error rates in the evaluation but to a highly robust behavior in practical application.

Results
The runtime of the implemented features is shown in Table 9 and of the models in Table 10. The runtime of the decision tree is listed for the shortest and longest branch.
The time constant T of the EMA for the features MCR2, WFL2S, WAM1 was 63.5 ms, and 127.5 ms for all other features.
For the logistic regression, NN and RNN, the feature calculation was followed by the normalization to ensure that the features were in an equal value range. For the decision tree, no normalization was required. The coefficients were determined with the features calculated for the training data. These quantized coefficients were implemented on the µC. The resulting runtime for each feature normalization was 0.53 µs.

Selected Models
The choice of the final model depends on the requirements. For a highly robust system, the author suggests the RNN with 9 hidden units and 9 features. Figure 7 shows the output of the selected model for a window of the test data. The confusion matrix of the test data for this RNN is presented in Table 11.
For an application with minimum calculation effort but high performance, the decision tree with 3 features (SSC3, ZCR2, VAR1S) and 4 splits is the best choice. In comparison to the RNN, the decision tree has a longer time lag at the transitions (see Figure 8), which leads to lower accuracies. By increasing the signal delay, the accuracy could be improved; however, in this article, long transitions were recorded as an error, because it was aimed at a fast system. The confusion matrix of the selected decision tree is presented in Table 12. Please note that the decision tree has less test samples due to a lower sampling frequency, see Section 2.3. However, the durations of the test data time sequences of the decision tree and RNN are equal. Table 11. Confusion matrix for selected RNN model (9 hidden units, 9 input features).

Contraction Artifact
Actual

Conclusions
Prosthesis drive activation due to artifacts is highly unpleasant for amputees. Avoiding incorrect activation leads to higher performance and thereby increases the acceptance of the myoelectric prostheses. With the presented models, the artifacts and contraction EMG signals were successfully distinguished with high accuracies. These algorithms substantially improved the robustness of the flexible insulated EMG sensors in real-world environment. Activation of the prosthesis drive due to artifacts is now effectively prevented by the implemented models. This high performance was achieved also by the algorithms' fast reaction to signal changes. All the models were designed to require low calculation effort for implementation on low-power wearable systems. With the shallow models, high accuracies were accomplished, although many transitions of the classes were included in the test data. The algorithms led to a short time between decisions as no Fourier transform had to be calculated. The presented models use hand-crafted time domain features only. Due to the wisely designed features, high accuracies were achieved with shallow models.
The author suggests an RNN or a regular NN for highest accuracies (99.91%/99.06%). For an implementation, which requires very little calculation resources, the author suggests a decision tree (98.76% accuracy).

Conclusions
Prosthesis drive activation due to artifacts is highly unpleasant for amputees. Avoiding incorrect activation leads to higher performance and thereby increases the acceptance of the myoelectric prostheses. With the presented models, the artifacts and contraction EMG signals were successfully distinguished with high accuracies. These algorithms substantially improved the robustness of the flexible insulated EMG sensors in real-world environment. Activation of the prosthesis drive due to artifacts is now effectively prevented by the implemented models. This high performance was achieved also by the algorithms' fast reaction to signal changes. All the models were designed to require low calculation effort for implementation on low-power wearable systems. With the shallow models, high accuracies were accomplished, although many transitions of the classes were included in the test data. The algorithms led to a short time between decisions as no Fourier transform had to be calculated. The presented models use hand-crafted time domain features only. Due to the wisely designed features, high accuracies were achieved with shallow models.
The author suggests an RNN or a regular NN for highest accuracies (99.91%/99.06%). For an implementation, which requires very little calculation resources, the author suggests a decision tree (98.76% accuracy).

Conclusions
Prosthesis drive activation due to artifacts is highly unpleasant for amputees. Avoiding incorrect activation leads to higher performance and thereby increases the acceptance of the myoelectric prostheses. With the presented models, the artifacts and contraction EMG signals were successfully distinguished with high accuracies. These algorithms substantially improved the robustness of the flexible insulated EMG sensors in real-world environment. Activation of the prosthesis drive due to artifacts is now effectively prevented by the implemented models. This high performance was achieved also by the algorithms' fast reaction to signal changes. All the models were designed to require low calculation effort for implementation on low-power wearable systems. With the shallow models, high accuracies were accomplished, although many transitions of the classes were included in the test data. The algorithms led to a short time between decisions as no Fourier transform had to be calculated. The presented models use hand-crafted time domain features only. Due to the wisely designed features, high accuracies were achieved with shallow models.
The author suggests an RNN or a regular NN for highest accuracies (99.91%/99.06%). For an implementation, which requires very little calculation resources, the author suggests a decision tree (98.76% accuracy).
As mentioned above, the linear separation algorithm presented by Roland et al. [16] has difficulties in detecting short contractions. This algorithm is now outperformed by the new models. Additionally, the need for the calculation of the short-time Fourier transform is obviated.
In the next step, models will be trained on amputees and the EMG sensor system with the presented motion artifact suppression will be applied to real-world prosthesis control. Furthermore, the potential of improving the performance by a random forest or a convolutional NN will be evaluated. The accuracies of ensemble models will be investigated. The runtime of the features and the models will be included to the cost term and a global optimization problem, which also comprises the feature parameters, will be solved. Additionally, the application of these algorithms to the state-of-the-art conductive EMG sensors and to other biosignal sensors (e.g., EEG or ECG) will be investigated. Furthermore, a multiple subject dataset will be acquired, and the presented models will be trained with the new data. The effect of inter-subject variability in the EMG signal will be evaluated. It will be examined whether a model trained on multiple subjects or a model trained on the individual EMG data leads to the best performance.