Real-time modeling and feature extraction method of surface electromyography signal for hand movement classification based on oscillatory theory

A total of eight healthy volunteers (six males and two females, aged between 23 and 33 years old) participated in the study. Their information are listed in table 1. Before the experiment, the leader of the experiment introduced the experiment content to each participant, and they signed an ethical agreement and a non-disclosure agreement with the ethics committee of Hefei University of Technology. The ethics committee approved all the undertaken work prior to data collection. In addition, the whole experiment process was non-invasive and safe. Each participant fully understood the experiment content and participated in the experiment voluntarily.

Before the experiment, an electrostatic elimination course was performed for each participant, and the participants were informed of the whole experiment process and experiment specification. In addition, participants with a lot of body hair were shaved and swabbed with alcohol. The experiment was conducted in an empty room.

A 16-channel analog to digital (AD) acquisition card (USB1252A, Smacq Ltd), five-channel three-pole (positive, negative, and reference poles) dry sEMG sensors (20 mm electrode spacing, Qingdao Zhituo Ltd), a 12 V battery, and a direct current to direct current (DC–DC) module constituted the experimental platform, as shown in figure 2. The maximum allowable voltage of the sEMG sensor was 3.3 V. The signal acquisition software was developed based on LabVIEW development framework, as shown in figure 2. After rubbing the corresponding surface with alcohol, the dry sEMG sensor was placed at the corresponding position and secured with medical pressure-sensitive tape (Heiner Ltd Company). The flexor pollicis longus, flexor digitorum superficialis, flexor carpi ulnaris, extensor digitorum and extensor carpi radialis brevis were chosen as the main muscles according to [21]. The sampling rate was 2000 Hz. The corresponding actions to be identified were M0–M5, as shown in figure 3. The experimental flow chart is shown in figure 4. During the experiment, a metronome was used to regulate the motion duration and resting time of the participants. When the participant was ready, the experiment was started. Each participant was required to conduct the experiment ten times, and each experiment was divided into three sessions. Session 1 and session 2 were separated by 5 h. Session 2 and session 3 were separated by 1 d to avoid muscle fatigue. Each session contained a calibration part and ten trial parts. Each trial consisted of 20 random sets of six actions. The ratio of the calibration and trial parts was 7:3. Each trial consisted of six random actions, each lasting 2 s, then resting for 2 s, as shown in figure 4.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The method developed in this work is shown in figure 5. The sEMG signal was firstly band-pass filtered with 10–500 Hz (fourth-order Butterworth band-pass filter). The RMS signal of the sEMG was used to distinguish between the stationary state and the motion state, which were introduced in our previous work [20]. The sEMG modeling method is partly based on the oscillatory mechanism. Based on the combination of multiple linear Fourier functions [31, 32], attenuation factors were added to form the following basic function. The curve of ${e^{ - v \,{{\left( {k\Delta t} \right)}^2}/2}}$ over time is shown in figure 6. In this work, $v$ is equal to 0.3. When the motion state was detected, the sEMG modeling and feature extraction course could be conducted:

$$\begin{equation}{x_{jk}} = \sin \left( {2\pi \left( {{f_0} + \frac{j}{G}} \right)k\Delta t} \right){e^{\frac{{ - v \,{{\left( {k\Delta t} \right)}^2}}}{2}}}{ }\end{equation} \tag{ 1 }$$

$$\begin{equation}{z_{jk}} = \cos \left( {2\pi \left( {{f_0} + \frac{j}{G}} \right)k\Delta t} \right){e^{\frac{{ - v \,{{\left( {k\Delta t} \right)}^2}}}{2}}}{ }\end{equation} \tag{ 2 }$$

$$\begin{equation}{{\boldsymbol{{a}}}_{\boldsymbol{{k}}}} = \left[ {{x_{1k}},{x_{2k}}, \ldots ,{x_{Jk}},{z_{1k}},{z_{2k}}, \ldots ,{z_{Jk}}} \right]\,\end{equation} \tag{ 3 }$$

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

where j is the fractional factor index, and its value ranges from 1 to J. k represents the data sampling index, ranging from 1 to N. N is the number of signal data points. As the sampling rate was 2000 Hz, the parameter $\Delta t$ is equal to 0.0005. G represents the frequency division factor. The larger its value is, the more intensive the frequency division is, as shown in figure 7. In this figure, the horizontal coordinate represents the signal frequency of concern. ${f_0}\,$ and $f\,$ are the minimum and maximum frequencies of the sEMG of each participant. ${f_0}$ and $f$ of each participant were obtained by measuring the sEMG signal corresponding to each hand movement and analyzing the spectrum preliminarily. Apparently, J = $\left( {\,f - {f_0}} \right)G$. The optimal value of G will be discussed in section 3. Equations (1) and (2) are the basis functions used to construct the sEMG signal. The basic function consists of the vector ${{\boldsymbol{{a}}}_{\boldsymbol{{k}}}}$ as shown in equation (3). Multiple component factors constitute the estimated sEMG signal, as shown in the following formula:

$$\begin{equation}{\hat y_{m,k}} = \mathop \sum \limits_{j = 1}^J {w_{mjk}}{x_{jk}} + {v_{mjk}}{z_{jk}}{ }\end{equation} \tag{ 4 }$$

Zoom In Zoom Out Reset image size

where m represents the number of muscles, m is from 1 to 5. ${w_{mjk}}$ and ${v_{mjk}}$ are weights to be estimated, and the weight vector ${{\boldsymbol{{w}}}_{{\boldsymbol{{m}}},{\boldsymbol{{k}}}}}$ is constituted as follows:

$$\begin{equation}{{\boldsymbol{{w}}}_{{\boldsymbol{{m}}},{\boldsymbol{{k}}}}} = \left[ {{w_{m,1k}},{w_{m,2k}}, \ldots \!,{w_{m,\,Jk}},{v_{m,1k}},{v_{m,2k}}, \ldots \!,{v_{m,Jk}}} \right].\end{equation} \tag{ 5 }$$

The linear Kalman filter [30] was adopted as the weight-updating method, as follows:

$$\begin{equation}{{\boldsymbol{{K}}}_{{\boldsymbol{{m}}},{\boldsymbol{{k}}}}} = \frac{{{{\boldsymbol{{P}}}_{{\boldsymbol{{m}}},{\boldsymbol{{k}}}}}{{\boldsymbol{{a}}}_{\boldsymbol{{k}}}}}}{{{\boldsymbol{{a}}}_{\boldsymbol{{k}}}^{\mathbf{T}}{{\boldsymbol{{P}}}_{{\boldsymbol{{m}}},{\boldsymbol{{k}}}}}{{\boldsymbol{{a}}}_{\boldsymbol{{k}}}} + R}}\end{equation} \tag{ 6 }$$

$$\begin{equation}{{\boldsymbol{{w}}}_{{\boldsymbol{{m}}},{\boldsymbol{{k}}} + 1}} = {{\boldsymbol{{w}}}_{{\boldsymbol{{m}}},{\boldsymbol{{k}}}}} + {{\boldsymbol{{K}}}_{{\boldsymbol{{m}}},{\boldsymbol{{k}}}}}\left( {{y_k} - \mathop \sum \limits_{m = 1}^{\boldsymbol{{M}}} {{\hat y}_{m,k}}} \right)\end{equation} \tag{ 7 }$$

$$\begin{equation}{{\boldsymbol{{P}}}_{{\boldsymbol{{m}}},{\boldsymbol{{k}}} + 1}} = \left[ {{\boldsymbol{{I}}} - {{\boldsymbol{{K}}}_{{\boldsymbol{{m}}},{\boldsymbol{{k}}}}}{\boldsymbol{{a}}}_{\boldsymbol{{k}}}^{\mathbf{T}}} \right]{{\boldsymbol{{P}}}_{{\boldsymbol{{m}}},{\boldsymbol{{k}}}}} + {\boldsymbol{{Q\,\,}}}\end{equation} \tag{ 8 }$$

where ${{\boldsymbol{{K}}}_{{\boldsymbol{{m}}},{\boldsymbol{{k}}}}}$ represents the Kalman gain of the mth component at time k, ${{\boldsymbol{{P}}}_{{\boldsymbol{{m}}},{\boldsymbol{{k}}}}}$ represents the error matrix of the mth component at time k, ${\boldsymbol{{Q}}}$ represents the prediction noise covariance matrix, $R$ represents the measurement noise covariance matrix, and ${\boldsymbol{{I}}}$ is the identity matrix.

The linear Kalman filter is easy to calculate and can be used to update the signal weight, which has been widely used. Therefore, the linear Kalman filter was adopted in this study as shown in equations (6)–(8). Compared with the traditional least mean square (LMS) algorithm [33], its calculation accuracy is better. We found that calculating the Kalman filter was one of the most time-consuming parts of the algorithm. When the processor is weak, the LMS method can also be considered to replace the Kalman filter to update the weight.

Based on the sEMG modeling method mentioned above, real-time data of each component signal can be obtained. The weights ${w_{mjk}}$ and ${v_{mjk}}$ of each component signal were used to construct the following features of the sEMG, as shown in equation (9):

$$\begin{equation}{F_{m,j}} = \sum\limits_{k = 1}^N w_{m,jk}^2 + v_{m,\,jk}^2.\end{equation} \tag{ 9 }$$

As the weights ${w_{mjk}}$ and ${v_{mjk}}$ represent the variation trend of multiple linear Fourier functions, as shown in equations (1) and (2), the feature can reflect the time-frequency characteristics of the sEMG. Combined with common RMS features, the feature set required for subsequent action recognition can be formed. The number of features of each sEMG channel is J. The obtained features were sorted from large to small, and the ReliefF [34] algorithm was used for feature selection and to eliminate unnecessary features. The ReliefF algorithm assigns different weights to features based on the correlation of each feature and category. Features with weights less than a certain threshold value will be removed.

The Adaboost method is a common supervised classification algorithm [35], and is an important one in ensemble learning. By integrating multiple sub-classifiers (such as the decision tree), the optimal classification results can be obtained according to the results of multiple sub-classifiers combined with decision methods. The algorithm was proposed based on the boosting algorithm by Freund and Schapire [35]. By selecting the weak classifier with the lowest weight coefficient from the trained weak classifier and adjusting the sample weight and weight of the weak classifier, a final strong classifier can be combined from them. In this work, the decision tree was used as the weak classifiers. The pseudo code given below is taken/adapted from [35].

Adaboost [35]

Given $\left( {{x_1},{y_1}} \right),\, \ldots ,\,\left( {{x_m},{y_m}} \right)$ where ${x_i} \in \aleph ,\,{y_i} \in \left\{ { - 1, + 1} \right\}$

Initialize: ${D_1}\left( i \right) = 1/m{\text{ for }}i = 1, \ldots ,m$

For t = 1,...,T:

Train weak learner using distribution ${D_t}$.

Get weak hypothesis ${h_t}:\,\aleph \to \left\{ { - 1, + 1} \right\}$.

Aim: select ${h_t}$ with low weighted error:

$$\begin{equation}{\varepsilon _t} = {\text{P}}{{\text{r}}_{i\sim {D_t}}}\left[ {{h_t}\left( {{x_i}} \right) \ne {y_i}} \right].\end{equation} \tag{ 10 }$$

Choose ${\alpha _t} = \frac{1}{2}\ln \left( {\frac{{1 - {\varepsilon _t}}}{{{\varepsilon _t}}}} \right)$

Update, for i = 1,...,m:

$$\begin{equation}{D_{t + 1}}\left( i \right) = \frac{{{D_t}\left( i \right)\exp \left( { - {\alpha _t}{y_i}{h_t}\left( {{x_i}} \right)} \right)}}{{{Z_t}}}\end{equation} \tag{ 11 }$$

where ${Z_t}$ is a normalization factor (chose so that ${D_{t + 1}}$ will be distribution).

Output the final hypothesis:

$$\begin{equation}H\left( x \right) = {\text{sign}}\left( {\sum\limits_{t = 1}^T {\alpha _t}{h_t}\left( x \right)} \right).\end{equation} \tag{ 12 }$$

The root mean square difference (RMSD) was used to evaluate the performance of the sEMG signal modeling. The smaller the RMSD value was, the better the modeling performance was:

$$\begin{equation}{\text{RMSD}} = \sqrt {\frac{{{{\mathop \sum \nolimits }}_{k = 1}^N{{\left( {{y_k} - {{\hat y}_k}} \right)}^2}}}{{{{\mathop \sum \nolimits }}_{k = 1}^Ny_k^2}}} .\end{equation} \tag{ 13 }$$

The accuracy, precision, recall and F1score were used as evaluation indices of hand movement recognition, and the corresponding calculation formulas are as follows:

$$\begin{equation}{\text{Accuracy}} = \frac{{{\text{TP}} + {\text{TN}}}}{{{\text{TP}} + {\text{FP}} + {\text{TN}} + {\text{FN}}}} \times 100\% { }\end{equation} \tag{ 14 }$$

$$\begin{equation}{\text{Precisio}}{{\text{n}}_\zeta } = \frac{{{\text{T}}{{\text{P}}_i}}}{{{\text{T}}{{\text{P}}_i} + {\text{F}}{{\text{P}}_i}}} \times 100\% { }\end{equation} \tag{ 15 }$$

$$\begin{equation}{\text{Recal}}{{\text{l}}_\zeta } = \frac{{{\text{T}}{{\text{P}}_i}}}{{{\text{T}}{{\text{P}}_i} + {\text{F}}{{\text{N}}_i}}} \times 100\% { }\end{equation} \tag{ 16 }$$

$$\begin{equation}{\text{F}}1{\text{scor}}{{\text{e}}_\zeta } = 2 \times \frac{{{\text{Precisio}}{{\text{n}}_i} \times {\text{Recal}}{{\text{l}}_i}}}{{{\text{Precisio}}{{\text{n}}_i} + {\text{Recal}}{{\text{l}}_i}}} \times 100\% \end{equation} \tag{ 17 }$$

where the subscript $\zeta$ indicates indices of different actions, and $\zeta$ = 1,2...6 correspond to M0–M5. TP means true positive, FP is false positive, TN is true negative, and FN is false negative.

When f₀ and f are determined, the value of the parameter G will affect the frequency band subdivision degree, modeling estimation accuracy and real-time performance of the algorithm. The relationship between the RMSD values and G value is shown in figure 8(a), and (b) shows the relationship between the execution time and G value. Figure 8(a) shows that the RMSD value will decrease with the increase of G; figure 8(b) shows that the calculation time will rapidly increase with the increase of G. After comprehensive consideration, the optimal G value of this study was 0.5.

Zoom In Zoom Out Reset image size

Figure 9 shows an example from the sEMG signal modeling results, where G is 0.5. Figure 9(a) represents the estimation of the EMG signal (not normalized). In order to better show the sEMG signal estimation effect of the method developed in this work, part of the actual sEMG signal and the estimated sEMG signal are zoomed. Figure 9(b) shows some sub-signals of the sEMG signal. Figure 9(c) is the frequency spectrum of the sub-signals of the sEMG signal. It can be seen from figure 9(c) that each sub-signal has different frequency characteristics. Combined with figure 9(a), it can be seen that the developed method can achieve a real-time signal separation effect, and the corresponding RMSD value is 0.2318. The sEMG signal estimation results of all experimental objects are shown in figure 10. According to the experimental results, the average RMSD result is 0.3838 ± 0.0591.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Based on the sEMG signal modeling method, relevant features are extracted using equation (9), and feature sets are formed by combining RMS features. Before estimating the hand movements, the feature sets were normalized as shown in figure 11(a), so the ordinate has no unit. Four sample feature signals are listed in the figure as an example. The actual number of feature signals is much larger than 4, and the ReliefF [34] algorithm was used for feature selection and to eliminate unnecessary features. It can be seen from the figure that the extracted feature signal can better reflect the trend of the signal. When the signal is strong, the feature signal is relatively strong; when the muscle is at rest, the feature signal also tends to zero. Figures 11(b)–(f) represent the value of the feature signal in the single movement duration of the six actions of channel 1–channel 5 corresponding to M0–M5. It can be seen that the feature values of the five channels corresponding to different actions are obviously different and easy to distinguish. Combined with the feature selection method, some redundant features can be removed to ensure the classification accuracy and reduce the computational burden.

Zoom In Zoom Out Reset image size

Figure 12 shows the confusion matrix corresponding to different hand movements of all participants. In figure 12, the main diagonal represents the accuracy of the actual action and the predicted action, shown by grayscale proportional to the accuracy. The horizontal axis represents the accuracy of M0, M1, M2, M3, M4 and M5 for the actual action. It can be seen from the figure that the accuracy is higher than 90% except for M3 of S1 and M3 of S7. The average recognition accuracy of each action (M0–M5) is shown in figure 13. It can be seen that M5 has the highest recognition accuracy of 98.36 ± 1.99%, while M3 has the worst classification effect of 92.22 ± 4.91%. Relatively, the M3 action has the largest standard deviation, followed by the M2 action with 95.70 ± 2.83%. Figure 14 shows the average accuracy of all participants. The statistical results show that S5 has the highest average accuracy of 98.17 ± 1.60%, while S1 has the lowest average accuracy of 94.03 ± 1.29%. The average experimental results of all participants are 96.03 ± 1.74%.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The precision, recall and F1score values corresponding to all experimental results are listed in tables 2–4, respectively. The average precision results of all participants are 98.4748 (M0), 97.9882 (M1), 95.9361 (M2), 91.8710 (M3), 94.0819 (M4), 97.7501 (M5), respectively. The average recall values of all participants are 98.0091 (M0), 97.9449 (M1), 95.7000 (M2), 92.2210 (M3), 93.7771 (M4), 98.3608 (M5). The average F1scores of all participants are 98.2402 (M0), 97.9510 (M1), 95.8013 (M2), 92.0095 (M3), 93.8902 (M4), 98.0482 (M5). According to the experimental results, the method developed in this study has good recognition performance for different movements.

In order to compare the accuracy and execution time aspects of the developed method and the existing hand motion intention recognition method, table 5 shows the relevant comparison results. Among them, the WNN with the WT method comes from [11], and the WT method is used to decompose sEMG signals. In the tunable Q wavelet (TQWT) with Kraskov entropy (KRE) and k-Nearest Neighbor (K-NN) method [34], the TQWT method is used to decompose the sEMG signals and extract the KRE features of each sub-signal. EMD with RMS and the bagging method is partly from the literature [36]. The EMD method is used to decompose the sEMG signals. RMS features were extracted from the decomposed sEMG signals. The VMD method combined with the composite permutation entropy index (CPEI) and the bagging method comes from [17]. The particle swarm optimization method [11] was adopted to search for the optimal parameters of each method. All these methods used the same dataset and classification task in this work. In order to compare the computational complexity of each method, the time required for actual 2 s data (2000 Hz) was calculated. The computing platform was a computer equipped with Intel (R) Core (TM) i5-8250U CPU. The running environment is Matlab 2016Rb. The comparison results in table 5 show that the proposed method has the highest recognition accuracy and good real-time performance.