Spatio-temporal warping for myoelectric control: an offline, feasibility study

The DTW method measures the similarity between two sequences of different temporal dynamics. It has been used in speech recognition [47], data mining [48], gesture recognition [49], robotics [50], and myoelectric control [51–53]. Previous research utilising DTW in myoelectric control considered time-series similarity across training and testing EMG sequences from individual channels only. Using DTW, we created a similarity matrix between all pairs of EMG signals. DTW can also be applied across a set of signals together.

Assume two temporal sequences of length l: $\mathbf{a} = \left\{a_1,a_2,\ldots, a_l\right\}$ and $\mathbf{b} = \left\{b_1,b_2,\ldots, b_l\right\}$, representing two EMG signals. We define $\mathbf{D}(\textbf{a},\textbf{b})$ as an l × l distance matrix between a and b, with $\mathbf{D}_{ij} = (a_i - b_j)^2$. In the case of Euclidean distance, we generally set i = j, that is the distance along the same index from the two time series. The path for the Euclidean distance is hence the one along the diagonal of the matrix D, where i = j. Although, due to its simplicity and efficiency, Euclidean distance is widely used in different approaches, it suffers from two main problems: first, it is not practical for situation that the length of the sequences is not equal and second, it is sensitive to time shifting. Finding an optimal path using DTW can address these challenges. In the case of DTW method, the distance is calculated through a different warping path $\mathcal{P}$ generated by traversing the matrix D along ordered pairs of positions: $\mathcal{P} = \lt(e_1,f_1), (e_2,f_2),\cdots, (e_l,f_l)\gt$ where $e_i \in [1:l]$ and $f_i \in [1:l]$ are the positions that make the warping path. A valid warping path must satisfy the conditions $(e_1, f_1) = (1,1)$, $(e_l, f_l) = (l, l)$, $0~\leq e_{i+1} - e_i \leq 1$ and $0~\leq f_{i+1} - f_i \leq 1$ for all i < l.

A warping limit is typically imposed on the distances, that is $\left|e_i -f_i\right| \le w.l$, $\forall(e_i,f_i) \in P^*$. The value of w is the maximum allowed warping path to deviate from the diagonal. According to [54], the distance D to any path $\mathcal{P}$ is given with $D_\mathcal{P}(\textbf{a},\textbf{b}) = \sum_{i = 1}^{s} p_i$, where p_i = $\mathbf{D}_{{e_i},{f_i}}$ is the distance between element at position e_i of a and at position f_i of b for the ith pair of points in a proposed warping path P. Considering a space of all possible paths as $\mathcal{P}$, then the DTW path $P^*$ is the one that has the minimum distance: $P^* = \min_{P \in \mathcal{P}}~ \mathbf{D}_\mathcal{P}(\mathbf{a},\mathbf{b})$.

We utilised the aforementioned distances between all NC EMG signals, i.e. $(NC \times (NC-1))/2$ pairs, as shown in figure 2. Inspired by [19, 25] the DTW features of all pairs of first and second derivatives of EMG signals were also computed. Hence, our method extracts $3~\times (NC \times (NC-1))/2$ features, in total. All calculated DTW distances were then concatenated to form a features vector and passed through logarithmic scaling for improved accuracy, similar to [55].

We first multiplied the features vector extracted from the current windows, DTW_t , with that extracted from the previous window $DTW_{t-n}$. To account for the long-term memory, our algorithm borrowed the cell state (c) concept of the LSTM method. This was the information highway that passes through all the LSTM cells across all time steps to evolve a recursive representation of the features. It was therefore continuously updated within each cell. At the same time, the contents of each cell were updated with a weighted contribution of the cell state through the parameter β. The selection of the parameter β was empirically made with the objective of achieving the best classification performance across several datasets.

$$\begin{equation} f_t = DTW_t \odot DTW_{t-1} + \beta c \end{equation} \tag{ 1 }$$

where f_t representing the extracted features from the current time t, the symbol $\odot$ denotes the element-wise multiplication operation, and c indicates the cell state.

The STW algorithm presented a dynamic normalization mechanism; borrowing from neural machine translation [56]. Attention mechanism provides modeling of dependencies without considering their distance in the input or output sequences. We used this method within each cell to normalise the current windows features and before adding it to the aforementioned features to normalise the long-term memory component. The final features representation is

$$\begin{equation} {f_t}_i = \frac{{f_t}_i}{\sum_i {f_t}_i} + \mathrm{log} \left(1 + \frac{c_i}{\sum_i c_i}\right) \end{equation} \tag{ 2 }$$

where the subscript i indicates the normalization applied across features at time t.

Five EMG datasets were utilised to test the performance of the proposed feature extraction algorithm as described in table 1.

We used database-1 (DB1), database-5 (DB5), and database-7 (DB7) of the Ninapro repository [57–60], comprising respectively, 27 limbed-intact subjects (52 movements), 10 limb-intact subjects (52 movements and rest) and 20 limb-intact and two amputees (40 movements and rest). EMG signals were recorded using 10 Otto Bock electrodes and two wearable MYO armbands on the forearm in DB1 and DB5, respectively. In DB7, 12 Trigno (Delsys, USA) EMG sensors were placed on the limb. The third EMG dataset [41] included data from 22 limb-intact subjects each performing eleven hand/wrist gestures using two armbands. The fourth dataset, that is selective classification (SC-EMG), included data from five traumatic long trans-radial amputees [61], each performing seven classes of motion. Ten electrodes were placed around the forearm at the point of the largest diameter.

The below features were extracted from windows of 150 ms at 50 ms increments:

• HTD: Hudgins' TD feature set [18], comprising: mean absolute value (MAV), MAV slope, zero crossings (ZCs), slope sign changes, and waveform length (WL).
• AR-RMS: The 6th-order AR coefficients and the root-mean-square (RMS).
• LSF9: The lower sampling rate features defined in [22], comprising L-scale, maximum fractal length, the mean value of the square root, Willison amplitude, ZC, RMS, integrated absolute value, difference absolute standard deviation value, and variance.
• ATD: Combined AR with TD as defined in [60], comprising MAV, WL, 4th-order AR coefficients and log-variance (LogVar).
• STFS: Spatio-temporal features from [62], comprising integral square descriptor, normalised root-square coefficient of first and second differential derivatives, mean log-kernel, an estimate of mean derivative of the higher-order moments per sliding window, and a measure of spatial muscle information.
• fTDD: The fusion of TD features from [19] with six features representing the first three even power spectrum moments, with an irregularity factor, a sparsity measure, and the ratio of WL of the first derivative to that of the second derivative. The step size was 15.
• TSD: The temporal spatial-descriptors from [15].

To ensure a fair comparison, the dimensionality of all feature sets was reduced to c − 1, where c is the number of classes, using the spectral regression feature projection method [63].

The following classifiers were chosen: linear discriminant analysis (LDA), extreme learning machine (ELM), k-nearest neighbor (KNN), and support vector machines (SVMs). In the ELM classifier, one hidden layer with 1250 neurons was utilised. The parameters of SVM were optimised for each dataset, as the performance of SVM is susceptible to the kernel function parameter γ and the regularization parameter C, while for KNN, the parameter K was empirically set to K = 5. Additionally, three deep learning models were also implemented; LSTM, CNN, and a combination of CNN and LSTM. These models are shown in figure 3, designed empirically to reflect the best classification accuracy results using the raw EMG data. Raw EMG inputs to models like CNN (or even CNN + LSTM) are more successful than the spectrogram images [31, 41]. The RMS of raw EMG signals were therefore calculated, generating NC scalar values for each analysis window of 150 ms. These values were then turned into pseudo-images by multiplying each generated vector of size ($NC~\times$ 1) by its transpose, resulting in an NC × NC images. The images were scaled logarithmically and then provided as inputs to the CNN and CNN + LSTM models. For the LSTM model, the raw EMG samples for each 150 ms window of the NC channels were used as features.

Zoom In Zoom Out Reset image size

The Wilcoxon signed rank test was applied to verify the statistical significance of the achieved results, with the results being considered significant for a p-value smaller than 0.05. Finally, the size of differences observed was measured using Cohen's effect size d for paired samples defined as the difference between two group means divided by the standard deviation [22]. A set of predefined thresholds of 0.2, 0.5, 0.8, 1.2, and $\gt$2 are usually employed to equate the effect size to small, medium, large, very large, and huge effects respectively. Additionally, when reporting Cohen's d value for testing one method versus another, the direction of the effect was calculated by subtracting the mean of the latter from that of the former. MATLAB 2020a was utilised for all experiments on a laptop with an i7 processor, 16 GB of RAM, and a GPU unit (NVIDIA GeForce RTX 2060).

The average classification error rates across all folds and subjects are shown in figure 4(a). Furthermore, the corresponding bar plots with standard deviation values are shown in figure 1 of supplementary materials (available online at stacks.iop.org/JNE/18/066028/mmedia). In terms of the deep models, CNN + LSTM (noted with CNNL, for brevity) showed a significantly better performance than the individual models of CNN and LSTM (p < 0.001 for both tests, d = 3.15 for CNN vs. CNN + LSTM and d = 3.55 for LSTM vs. CNN + LSTM). CNN also significantly outperformed LSTM ($p \lt 0.001, d = 2.21$). This finding indicates that the temporal-spatial information captured by CNN + LSTM extracts more information than the individual CNN and LSTM, while the spatial information captured by CNN appears more important on this dataset than the temporal information captured by LSTM.

Zoom In Zoom Out Reset image size

In terms of the traditional hand-crafted feature extraction algorithms, as it is clear in figure 4(a), HTD and ARRMS were the worst performers on the 53 class problem in DB5, with all other methods significantly outperforming them, except for LSTM. When comparing HTD and LSTM, no difference was found ($p = 0.517,\ d = 0.10$), which is in line with the findings in [21].

We observed significant differences between LSF9 and ATD ($p = 0.030,\ d = 0.26$), and large differences between LSF9 and STFS ($p \lt 0.001,\ d = 0.80$) and LSF9 and TSD ($p \lt 0.001,\ d = 0.80$). These results also show that fTDD achieved significantly lower error rates than all other traditional and deep learning feature extraction methods, including CNN + LSTM, except the proposed STW features that significantly outperformed fTDD ($p\lt0.001, d = 2.64$). In contrast, STW significantly outperformed all other feature extraction and/or learning methods, including CNN + LSTM, with all tests having d > 2. Overall, the proposed approach yielded an average decrease in classification error of 19.81% ± 4.36% in comparison to all other methods considered in this work. Additionally, it showed a significantly lower standard deviation.

In addition to the classification accuracy metric, we evaluated the proposed STW method in comparison with other hand-crafted feature extraction algorithms using three other measurements including Fscore, Recall, and Precision [64]. Results are illustrated in figure 3 of supplementary materials. Average values for Fscore, Recall, and Precision, as well as statistical tests, prove that STW can outperform all conventional hand-crafted algorithms. Furthermore, to perform a quantitative comparison; in terms of separability of clusters, between STW and other conventional feature extraction methods, Davies-Bouldin index (DBI) was used [65]. Averaged DBI values calculated across test folds and all subjects are shown in figure 4 of supplementary materials. As it is visible clearly, STW has achieved the lowest DBI values which show its highest capability in feature discrimination criteria in comparison with other feature extraction methods.

The average classification error results on the SC-EMG datasets are shown in figure 4(b) for all traditional hand-crafted and deep learning models (Bar plots of the results are provided in figure 1 of supplementary materials). In terms of the deep models, it is interesting to see that, LSTM performed better than CNN in terms of the average error rates across all amputees. This could be in part attributed to the fact that the amputation affects muscles morphology and synergy and hence the performance of CNN is impacted by its capabilities in capturing spatial information between signals which may vary (i.e. warp) in timing as a result of amputation. In this direction, a recent work by Tallec and Ollivier [66] shows that LSTM networks have the capability to learn to warp input sequences. Despite the differences in the average classification errors of deep models, the Wilcoxon signed rank test though revealed no statistically significant difference between LSTM and CNN, LSTM and CNN + LSTM, and neither between CNN and CNN + LSTM ($p \gt$ 0.05 for all tests). In terms of the hand-crafted methods, fTDD performed significantly worse than TSD (p < 0.001) while the proposed STW significantly outperformed all other hand-crafted and deep models ($d = 2.54, 2.09, 1.87, 2.08, 1.87, 0.89$, and 2.05 for STW vs. HTD, AR-RMS, STFS, LSF9, ATD, TSD, and fTDD, respectively, and p < 0.001 for all comparisons). The performance of the STW method can be further justified by its ability to approximate the warping path that best aligns the two signals of interest in the time domain, supported by the DTW stage. Such a path is a more representative of the distance between the signals than what Euclidean distance can offer; despite nonlinear variations in speed. The interested reader in referred to [47] where the details of the DTW method can be found.

The same pattern of results was observed for DB7 datasets. Average classification error rates are shown in figure 4(c) as well as corresponding bar plots in figure 1 of supplementary materials. To avoid repetition, we used only an LDA classifier for decoding. HTD and AR-RMS were performed worst. In comparison to HTD and AR-RMS, LSTM performed better than both methods on DB7 datasets (p < 0.01 for both tests, d = 2.06 for HTD vs. LSTM, d = 1.15 for AR-RMS vs. LSTM). On the other hand, both CNN + LSTM and CNN significantly outperformed LSTM (p < 0.001 for both tests, d = 2.97 for LSTM vs. CNN + LSTM and d = 2.36 for LSTM vs. CNN), with CNN + LSTM also significantly outperforming CNN ($p \lt 0.001, d = 1.18$). On the other hand, both fTDD and the proposed STW significantly outperformed all other methods, including CNN + LSTM, while STW also significantly outperforming fTDD ($p \lt 0.001, d = 1.32$).

The average classification values achieved for the 3DC dataset are illustrated in figures 4(d), and (e) as well as in figure 1 of supplementary materials. As two different armbands including MYO and 3DC were used in this dataset, separate analyses were conducted. Instead of raw EMG signals, extracted STW features were fed into the deep models. KNN achieved the lowest average classification error, but statistical tests reveal no significant differences. Similar to the SC-EMG dataset, LSTM could achieve the best performance compared to CNN + LSTM, and CNN in terms of average classification error rates. Wilcoxon signed rank test revealed significant differences between LSTM and CNN and also between LSTM and CNN + LSTM (p < 0.001). With respect to the 3DC armband, similar results were obtained as LSTM significantly outperformed CNN and CNN + LSTM (p < 0.01).

Reporting average classification error or accuracy for the off-line analysis is biased and the comparison of results; achieved by a specific method, between off-line and real-time schemes will not be fair. To partially address this challenge, we used class-wise accuracy standard deviations (CWS) [67]. Averaged confusion matrices across all folds and subjects were calculated and the standard deviation of the main diagonal of the matrices was considered as CWS. In this analysis, DB5 and 3DC were considered.

Averaged confusion matrices across ten subjects and six folds using ELM, KNN, LDA, and SVM classifiers for the DB5 dataset are illustrated in figure 5. The proposed STW method could achieve the lowest CWS values among all feature sets. The results of CWS analysis for the 3DC dataset are presented in figure 6.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

For the DB5 and 3DC, we varied the feature extraction window size from 50 to 250 ms. The average classification error rates across all folds and subjects for DB5, 3DC with 3DC armband, and 3DC with MYO armband are shown in figures 7(a)–(c), respectively. Results show the classification errors decrease with larger windows, as expected. However, the results also show that varying the window size had a medium to large effect between consecutive windows, which could potentially indicate the robustness of the method against varying window sizes.

Zoom In Zoom Out Reset image size

The computation time for each of the feature extraction methods, except deep learning models, was calculated on randomly generated data with 150 samples across ten dimensions; equivalent to 150 ms with ten channels sampled at 1000 Hz. The time required to extract the features was registered and the analysis was repeated 1000 times. Average results are displayed in table 2.

Table 2. Time required to extract handcrafted features reported as average ± standard deviation.

Feature set	Time (ms)
LSF9	$0.8053 \,\,\pm\,\, 0.0620$
TSD	$0.4940 \,\,\pm\,\, 0.0411$
STFS	$0.4986 \,\,\pm\,\, 0.0381$
STW	$\boldsymbol{0.4248} \,\,\pm\,\, {0.0632}$
ATD	$0.2115 \,\,\pm\,\, 0.0237$
HTD	$0.1421 \,\,\pm\,\, 0.0180$
AR-RMS	$0.1660 \,\,\pm\,\, 0.0509$
fTDD	$0.0125 \,\,\pm\,\, 0.0354$

To evaluate the proposed STW method against recently published state-of-the-art features and classifiers a comprehensive comparison is performed on DB1 of the Ninapro database. For a fair comparison, a similar train-test approach was selected according to the [58, 68] in such a way that repetitions 1, 3, 4, 6, 8, 9, and 10 were selected to train the classifiers, and repetitions 2, 5, and 7 were used to evaluate methods. Moreover, the length of the window was selected equal to 200 ms. Average classification accuracy achieved by STW + SVM in comparison with previous traditional features as well as deep learning models is reported in table 3. Moreover, averaged DBI values achieved by conventional feature extraction methods and the proposed STW algorithm are illustrated in figure 4 of supplementary materials. Results show that STW has achieved the lowest DBI value in comparison to other methods.