Continuous Decoding of Daily-Life Hand Movements from Forearm Muscle Activity for Enhanced Myoelectric Control of Hand Prostheses*

To validate our method, we tried to learn a stable functional mapping between the forearm EMG activities and hand kinematics recorded during the execution of activities of daily living (ADLs). To this end, we used the KIN-MUS UJI dataset [29]. This dataset contains recordings of muscle activities (7 channels) and hand kinematics (18 DOFs) of 22 participants during the execution of 26 representative activities of daily living. Each activity was performed by the participants only once and lasted several seconds. Further details about the performed activities are reported for convenience in Table I.

The muscle activities were recorded using surface bipolar electrodes, whose locations were chosen to maximize the extraction of information generated by the forearm muscles [30]. The hand kinematics were recorded using an instrumented glove (CyberGlove, CyberGlove Systems LLC), which tracked fingers’ and wrist’s joint angles. The recorded ADLs include the 20 actions of the Sollerman Hand Function Test (SHFT) [31], which are commonly used to assess hand function in clinical settings and involve the interaction with objects of different sizes and weights, such as coins, cutlery, and irons. All the analyses reported in this work are based on the data from participants 1 through 20. Data from participants 21 and 22 presented missing data and were thus discarded.

Importantly, most proposed methods for the online mapping of forearm EMGs onto hand kinematics and for the classification of hand movements are validated using the Ninapro database [32]. However, the datasets included in this database mainly consist of simple and highly controlled movements, such as fingers’ flexion with static wrist posture and different types of grips. This is in stark contrast with the dataset we used in this work, where every task involves a markedly different form of interaction with daily living objects, in addition to a reaching and a release phase. For this reason, mappings learned on this dataset promise a superior ability to generalize to online settings.

The hand joint angles were acquired at a sampling rate of 100Hz. The EMG signals were acquired at a sampling rate of 1000Hz with a passband between 20Hz and 460Hz. The kinematic data were then zero-phase lowpass filtered with a Butterworth filter (cut-off frequency: 5Hz). To extract the EMG envelopes, the EMGs were rectified and zero-phase lowpass filtered with a Butterworth filter (half power frequency: 3Hz). The EMG signals were subsequently subsampled to 100Hz.

The reaching, manipulation, and release phases present in the original dataset were appropriately concatenated. Since previous studies have reported significant correlations between muscle activity and joint accelerations [33], we chose to adopt an acceleration-based kinematic representation. To this aim, we numerically differentiated the hand joint angles twice to compute joint accelerations. All the analyses present in this work were conducted in MATLAB (MATLAB 2020a, The MathWorks, Natick, MA).

An influential study on feature selection for the classification of EMG signals [34] reported that features based on signal energy and those based on frequency content tend to provide superior classification performance. Inspired by this work, we built an EMG feature set composed of EMG envelopes and EMG spectrograms. The spectrograms were computed on 50ms Kaiser windows within the frequency range 20Hz-400Hz; the resulting signals were upsampled in time to 100Hz and binned in frequency (500 bins).
Subsequently, we applied principal component analysis (PCA) to reduce the dimensionality of the EMG feature set to 25; the extracted projections retained about 99% of the original data variance.

Since EMG signals might give rise to different movements depending on the starting kinematic properties, such as posture [15], we extended our feature set by also including information about the previous movements. Specifically, we included the hand joint accelerations observed 30ms in the past. The complete feature set thus includes EMG features and kinematic features. Finally, the resulting features were z-scored.

To avoid over-reliance on the kinematic features, we augmented the training dataset by adding, for each action, 60 additional training examples. In these examples, the kinematic features were corrupted by noise, while the EMG features were left unchanged. Out of the 60 additional examples, 30 were corrupted by white Gaussian noise and 30 by colored noise with a maximum frequency of $10^{-7}{\rm rad}/{\rm sample}$; in both cases the maximum amplitude was $0.1$. Overall, this procedure ensures that the networks learn robust features [35, 36].

To learn the mapping between the complete feature set and hand joint accelerations, we used LSTM-based networks [28]. This choice was motivated by previous work, which showed that the ability of gated recurrent neural networks to learn long-term relationships between signals is instrumental in allowing the prediction of whole-body human motion [37, 38]. Moreover, LSTM-based networks have recently been proved to be effective at predicting wrist position (3 DOFs) [39] and wrist flexion/extension (1 DOFs) [40] by tracking the EMG activity of 2-5 shoulder and arm muscles.

In brief, LSTM networks are recurrent networks of special units. Each unit is endowed with two memory cells (namely, an output memory ${\bf y}$ and a hidden memory ${\bf h}$) and three gates (namely, an input gate ${\bf i}$, a forget gate ${\bf f}$, and an output gate ${\bf o}$). The hidden memory cell ${\bf h}$ allows the network to store and use important input features indefinitely. The input gate ${\bf i}$ modulates the extent to which incoming information is retained; the forget gate ${\bf f}$ modulates the extent to which old information is erased; finally, the output gate $\bf o$ extracts information from the hidden memory useful to generate an appropriate output ${\bf y}$ for the current input ${\bf u}$.

The output, once computed, is stored in the output memory cell.

The gating mechanism is defined by the following equations:

	$\displaystyle{\bf f}^{(t)}$	$\displaystyle=\sigma({\bf A}_{f}{\bf y}^{(t-1)}$	$\displaystyle+{\bf B}_{f}{\bf u}^{(t)}$	$\displaystyle+{\bf b}_{f})$		(1)
	$\displaystyle{\bf i}^{(t)}$	$\displaystyle=\sigma({\bf A}_{i}{\bf y}^{(t-1)}$	$\displaystyle+{\bf B}_{i}{\bf u}^{(t)}$	$\displaystyle+{\bf b}_{i})$		(2)
	$\displaystyle{\bf o}^{(t)}$	$\displaystyle=\sigma({\bf A}_{o}{\bf y}^{(t-1)}$	$\displaystyle+{\bf B}_{o}{\bf u}^{(t)}$	$\displaystyle+{\bf b}_{o})$		(3)

where $t$ is the current time step and $\sigma()$ is the element-wise sigmoid function.

The hidden memory content is modified according to

	$\displaystyle{\bf{\tilde{h}}}^{(t)}$	$\displaystyle=\tanh({\bf A}_{h}{\bf y}^{(t-1)}+{\bf B}_{h}{\bf u}^{(t)}+{\bf b}_{h})$		(4)
	$\displaystyle{\bf h}^{(t)}$	$\displaystyle={\bf f}^{(t)}\odot{\bf h}^{(t-1)}+{\bf i}^{(t)}\odot{\bf\tilde{h}}^{(t)}$		(5)

where $\odot$ denotes the element-wise product operator. Finally, the output is determined by

$${\bf y}^{(t)}={\bf o}^{(t)}\odot\tanh({\bf h}^{(t)})$$

(6)

In this work, we considered architectures composed of $n_{d}=\{1,2,3\}$ LSTM layers of $n_{h}$ units projecting onto a fully connected layer with 36 units, followed by a dropout layer ([41]) with probability $0.5$, and a final fully connected layer of 18 units. The dynamics of each LSTM layer are defined by equations (1-6). We chose not to consider deeper LSTM architectures since evidence suggests that the benefit of adding more than two layers is generally limited [42].

We used Bayesian hyperparameter optimization (bayesopt() — [43]) to tune $n_{h}$, initial learning rate, gradient threshold, and L2 regularization strength. The search ranges were $[30,600]$, $[10^{-4},10^{-1}]$, $[.2,1]$ and $[10^{-10},10^{2}]$, respectively. The available LSTM units were equally distributed among the $n_{d}$ layers. The selected cost function was the fraction of unexplained variance on the validation set.

Hyperparameter optimization automatically stopped when the maximum number of 100 cost function evaluations was reached. To efficiently evaluate a candidate hyperparameter set, we trained the parameters of the corresponding network for only three epochs. To further speed up the training time, we based hyperparameter optimization only on the data recorded from participant 1. Following this procedure, we selected, for each depth $n_{d}$, the three best networks in terms of validation error. Such networks were further trained for 25 epochs with subject-specific data only. Details about the parameter optimization are provided in the following section. Overall, this strategy allowed us to retrieve a good hyperpameter set reasonably fast: for example, for $n_{d}=1$, hypeparameter optimization was completed in less than one hour (i.e., $54^{\prime}47^{\prime\prime}$).

The networks were implemented and trained in MATLAB, using the Neural Network Toolbox. Parameter optimization was performed using adaptive moment estimation (Adam, [44]) with a piecewise learning rate schedule with a drop period of 10 epochs and a drop factor of 10%. To speed up training, we performed all the optimizations on a shared GPU (NVIDIA GeForce RTX 2070) using mini batches of 40 sequences. This configuration allowed us to fully train a network in a limited amount of time: for example, for $n_{d}=1$, the training speed was $8.25$ sec/epoch. Parameter optimization automatically stopped if the validation error did not decrease for three consecutive training epochs, or if the maximum number of 30 epochs was reached. For training, validation, and test we used the EMG signals and corresponding hand kinematics recorded during the execution of the tasks 1-20, 21-23, and 24-26, respectively. We refer to Table I for further details about the selected tasks.

To measure the reconstruction quality of the network’s predictions, we computed the Pearson’s correlation coefficients ($\rho$) between ground-truth and predicted angular accelerations, and the normalized reconstruction errors ($\epsilon$). The normalized reconstruction error is defined as the fraction of unexplained variance (i.e., $\epsilon=1-R^{2}$).

To assess the potential benefit of the recurrence present in LSTM networks, we also considered standard non-linear feedforward (FF) networks as a baseline. Specifically, we considered FF networks of depth $n_{d}=\{1,3,10\}$ with hyperbolic tangent activation function. A dropout layer with probability $0.5$ was added between each couple of layers. For each depth, similarly to what was done for LSTM networks, we performed Bayesian hyper-parameter optimization. In this case, we set the maximum number of units to 3000. All the other parameter ranges are as above.

To measure the robustness of the trained networks, we also assessed their ability to make accurate predictions in the presence of random noise corrupting the angular acceleration measurements. Specifically, we considered noise levels of amplitudes $\xi_{a}=\{0,1,5,10,15,30\}\%$. The amplitudes are defined as percentages of the maximum angular acceleration measured during a task, across all joints.

In the following sections, we report the test performance (measured on the independent activities of daily living 24, 25, and 26) of the resulting six networks: three LSTM and three FF networks.