1. Introduction
The neuromuscular system consists of the nervous system and the muscular system. The main function of the human muscular system is to provide the energy needed by the human body to perform various actions. Exercise-induced muscle fatigue is a physiological phenomenon in which the maximum voluntary contraction (MVC) capacity and output power of the muscle decreases. The cause of exercise-induced muscle fatigue is the accumulation of metabolites (lactic acid, hydrogen ion, inorganic phosphate) in the human blood during exercise [
1]. The risk of sports injuries increases with muscle fatigue. Therefore, the precise positioning of fatigued muscles is the basis for relieving and curing muscle fatigue and has important sports medicine significance. Electromyography (EMG) is an electrical signal generated by skeletal muscles when they contract spontaneously [
2]. EMG signal analysis can provide continuous measurement of the muscle’s state during the continuous fatigue contraction process, which is different from the subjective typical assessment that usually indicates that the subject can no longer perform the test [
3]. The surface electromyography signal (sEMG) is a comprehensive effect of superficial muscle EMG and nerve trunk electrical activity on the skin surface. sEMG can reflect neuromuscular activity to a certain extent, and it has the advantages of non-invasiveness and simple operation [
4]. The analysis of the sEMG signal is performed by a trained neurologist. When there are not enough experts to meet the needs of fatigue warning, it is very important to use deep learning technology to realize neuromuscular fatigue detection and classification on the basis of EMG signal processing.
However, the sEMG signals are affected by many confounding factors, including the following: the shape of the volume conductor; thickness of the subcutaneous tissue layers; distribution and size of motor unit areas in muscles; changes in transmembrane action potential; conductivities of the tissues; electrode position; detection system; etc. [
5]. Wavelet and wavelet packet transform are a multiresolution decomposition method that can be used to analyze signals and images denoising. Both wavelet transform and wavelet packet transform are multiresolution decomposition methods, which can characterize local features in the time domain and frequency domain, and they are usually used to analyze signals and image denoising [
6]. Li et al. analyzed and compared several different wavelet threshold denoising algorithms, and they proposed a new wavelet threshold denoising algorithm to effectively eliminate noise interference under strong noise background and preserve signal details [
7]. Traditional threshold functions include hard threshold and soft threshold functions, which have some deficiencies in signal de-noising. To solve these shortcomings, Zhang et al. proposed an improved wavelet threshold method to denoise MRI images [
8]. The results show that the improved wavelet threshold method in denoising MRI images has better performance than the traditional wavelet threshold function.
When the local muscle is fatigued, the sEMG signal of the muscle will change in the time domain and the frequency domain [
9]. Many scholars have studied muscle fatigue recognition models that can identify muscle fatigue states by detecting changes in sEMG signals. Latasa et al. detected the EMG signals of cyclists in an incremental continuous cycling test and used a multi-segment linear regression algorithm to study the aerobic-anaerobic threshold transition [
3]. Martinez et al. recorded the subjects’ sEMG signals in the fatigue state of the vastus lateralis in the incremental cycling test and analyzed the correlation between muscle fiber conduction velocity (MFCV), instantaneous mean frequency (iMNF), normalized root mean square (RMS), and muscle fatigue state [
10]. Subasi and Kiymik used independent component analysis (ICA) to process the surface EMG signal of the biceps brachii, and they used artificial neural network (ANN) to classify muscle fatigue [
11]. M-wave normalization of the sEMG signal can describe the activity amplitudes in the state of muscle fatigue during repeated sprinting [
12]. Wu et al. proposed a BFA–Gaussian support vector machine (GSVCM) model to improve the accuracy of muscle fatigue recognition [
13]. Hussain and Mamun utilized different wavelet functions (WFs) to analyze the EMG signal of the right rectus femoris muscle to identify the fatigue state of the right rectus femoris muscle during walking [
14]. Consequently, they used these functions for accurate automated muscle fatigue classification and systematically processing the EMG signal based on deep learning and machine learning methods.
This paper researches a new muscle fatigue state recognition model. First of all, we present a new wavelet packet threshold function denoising method to denoise sEMG signals. The signal-to-noise ratio (SNR) and root mean square error (RMSE) parameters of this method, hard threshold, and soft threshold denoising are compared. Secondly, we describe the time-domain and frequency-domain feature extraction applied in the sEMG signals classification process. We divided the sEMG signals into two groups: non-fatigue and fatigue, based on the anaerobic threshold (AT). Finally, the long and short-term memory (LSTM) network was used to identify muscle fatigue and compared with other classification algorithms in terms of classification performance.
2. Materials and Methods
2.1. Participants
Experimental participants: Twenty healthy males (n = 20; Age range: 22–33 years; Height: 1.79 ± 0.09 m; Weight: 65.5 ± 6.6 kg) were recruited to participate in this experiment. All participants were screened and free from cardiovascular, neuromuscular, and metabolic diseases. Before the test, all participants voluntarily participated in the test, fully understood the purpose, details, methods, and potential test risks of the experiment, and signed an informed consent form. They did not engage in vigorous exercise and did not consume caffeine, nicotine, and alcohol within 48 h before testing. The study was conducted in accordance with the Declaration of Helsinki and the National law of China.
2.2. Instrumentation
To ensure that the power is accurately regulated, the exercise test wase performed in the Lode Corival cpet cycle ergometer. The Ag–AgCl surface electrodes were arranged in bipolar configuration (20 mm center-to-center, 1 cm in diameter), and a Noraxon Ultium sEMG sensor were used to record the sEMG signal during the experiment. The sample frequency of the Noraxon Ultium sEMG sensors is 2000 Hz. We use Noraxon MR3 software to analyze the collected sEMG signals offline. The Ultima GX system was used to measure the participants’ ventilation (VE), oxygen uptake (VO2), and carbon dioxide production (VCO2) during the exercise test. We calibrated all the sensors and experimental equipment before the experiment.
The experimental data analysis was carried out on a workstation with an Intel CoreTM i7-9700 and 16 GB memory card. The preparation of the simulation programs was carried out on MATLAB 2020a.
2.3. Procedure
Before the experiment, in order to reduce the impedance, it is necessary to shave the excess hair of the skin of the participant’s legs and wipe it with medical alcohol. Then, we placed sEMG sensors on the appropriate position of the vastus rectus femoris (RF), vastus lateralis (VL), vastus medialis (VM), and gastrocnemius (GA) muscles according to the SENIAM (Surface EMG for Non-Invasive Assessment of Muscles) guidelines [
15]. In order for participants not to be affected during the execution of the procedure, the sports bandages were used to fix the sEMG sensors. The position of the sEMG sensors on the left leg is shown in
Figure 1.
Before the test, participants could use equipment and procedures proficiently. After a 3 min warm-up on a cycle ergometer (Lode Corival cpet), each of the participants performed an incremental protocol in a temperate environment (25–28 ℃), starting at a 100 W initial workload with increases of 25 W every 1 min. The participants were instructed to maintain a pedaling rate within the range of 70–75 r/min throughout the test [
16,
17]. During the test, we strongly encouraged each participant to provide a maximal effort. The test was terminated when the participant was unable to maintain a pedaling rate above 70 r/min due to volitional exhaustion. The participant performs the test as shown in
Figure 2.
After the test, we use the V-slope method to calculate the AT based on analyzing the slopes of VO
2 and VCO
2 volume curves [
18]. We divide the sEMG signals into fatigue and non-fatigued status at the time corresponding to the AT.
2.4. Wavelet Packet Threshold Denoising
The wavelet packet threshold denoising algorithm uses the multi-scale characteristics of wavelet packet analysis to decompose the signal by wavelet packet. The threshold function determines whether the wavelet signals of different layers are noise signals according to the threshold. The wavelet coefficients related to noise become corresponding appropriate values according to a different number of layers. In the threshold denoising method based on wavelet transform, the selection of wavelet basis, the number of wavelet decomposition layers, the threshold, and the threshold function have a great influence on the effect of wavelet threshold denoising [
19]. According to previous research, the Daubechies (Db) wavelet family has the most suitable wavelet functions for sEMG signals denoising analysis [
20,
21]. Therefore, this paper uses the db45 wavelet function for the wavelet packet decomposition. This chapter will mainly describe the optimal wavelet packet decomposition layer algorithm and the improved wavelet threshold function.
2.4.1. Best Tree Wavelet Packet Analysis
The wavelet packet decomposes the sEMG signals to obtain high-frequency coefficients and low-frequency coefficients. Both the high-frequency coefficient and the low-frequency coefficient are decomposed by analogy as the input signal of the next stage until the decomposition reaches the set number of layers. The complete binary tree is produced as shown in
Figure 3. Before the decomposition of the wavelet packet, it is necessary to find the best tree of the wavelet packet. Starting with the root node, the best tree is calculated using Shannon entropy. The Shannon entropy calculation method is as follows.
where
is the wavelet packet coefficient sequence. The following method is used to calculate the best tree. A node N1 is split into two nodes N2 and N3 if and only if the sum of the entropy of N2 and N3 is lower than the entropy of N1; otherwise, node N1 will not be decomposed [
22]. This is a local criterion based on the information available at node N1, as shown in
Figure 4. It is calculated by the Shannon entropy criterion in which the decomposition of four layers is most suitable in the sEMG signals’ wavelet packet decomposition.
2.4.2. Threshold Function
The threshold function is essential for sEMG signal denoising analysis. The hard threshold and soft threshold functions proposed by Donoho [
23] are extensively used. The hard threshold and soft threshold functions are expressed as Equations (2) and (3), respectively.
The hard threshold function is calculated as follows:
The soft threshold function is expressed as follows:
Although the above-mentioned traditional signal denoising methods have some effects in practical applications, it is undeniable that they still have shortcomings. The disadvantage of the hard threshold function itself is discontinuous at threshold
, and the wavelet coefficients (
) larger than the threshold
are not processed, and the
smaller than the threshold
are set to zero, which will cause oscillations in the reconstruction of
. Although the soft threshold function is continuous and derivable at the threshold
, the wavelet coefficients processed by Equation (3) deviate from the actual wavelet coefficients. In order to compensate for the deficiencies of the above-mentioned threshold function, we proposed the improved threshold function as follows:
When
, Equation (4) can be written as:
When
, Equation (4) can be written as:
In Equations (2)–(4),
is the wavelet coefficient, and
λ is the threshold. Parameters
k and
m in Equation (4) include the adjustment parameters (
,
). The improved threshold function is not only continuous at the threshold but also high-order differentiable. The denoising effect of the wavelet packet threshold denoising algorithm is closely related to the choice of threshold and threshold function [
24]. In recent years, the following four threshold estimation criteria have been widely used: fixed threshold estimation (Sqtwolog), maximal minimum threshold estimation (Minimaxi), unbiased risk estimation (Rigsure), and heuristic threshold estimation (Heursure). In this paper, we select the Heursure estimation method to calculate the threshold
λ. The algorithm flow is as follows:
Square each element in the vector W, and then sort in ascending order to obtain a new vector (, , …, ). The length of the vector W is the integer N.
The threshold is the square root of the
i-th element of the vector
; then, the risk algorithm is as follows [
25]:
Calculate the
i value corresponding to the minimum risk (
i), as shown in Equation (8). The threshold λ calculation method is shown in Equation (9).
The procedure of the wavelet packet threshold denoising algorithm is as follows [
26]:
Select the appropriate wavelet function, and perform the wavelet packet decomposition calculation according to the best tree of the wavelet packet;
A threshold value is selected for the wavelet packet coefficients of each decomposition scale. Select the appropriate threshold function to process the wavelet packet coefficients ;
The processed wavelet packet coefficients are reconstructed by inverse wavelet packet transformation. The denoised signals are obtained.
The objective evaluation of the wavelet packet threshold denoising algorithm is described by the signal-to-noise ratio (SNR) and mean square error (MSE). The calculation methods of SNR and RMSE are shown in Equations (10) and (11).
In Equations (10) and (11), is the original signal, is the signal after denoising, and n is the signal length.
2.5. Feature Extraction
The sEMG signal is a 1D time-series signal of the neuromuscular system that is recorded on the skin surface. The analysis of sEMG concentrated on two main fields: the frequency domain and the time domain [
27]. In this study, we choose two time-domain features, Root Mean Square (RMS) and Integrated Electromyogram (IEMG), to describe the changes in the amplitude of the sEMG signals. The mathematical expressions of RMS and IEMG are expressed as:
The frequency-domain feature is extracted from the Fourier transform of the signal. We choose median frequency (MF) and mean power frequency (MPF) because of their abilities to reflect the fatigue-caused frequency changes of sEMG [
28,
29,
30]. The mathematical expressions of MF and MPF are expressed as:
In Equations (14) and (15),
f1 and
f2 determine the bandwidth, and
P(
f) is the power spectral density (PSD) of the sEMG signals.
P(
f) is expressed as:
In the above equations, L is the signal length, and x(f) is the sEMG signals in the frequency domain. The RMS, IEMG, MF, and MPF are extracted from the denoised sEMG signals by a moving window of 2 s. In order to construct a muscle fatigue feature data set, we divide the sEMG feature data into fatigue and non-fatigued sEMG data according to the time point of the AT. The V-Slope method was introduced to calculate the AT in the fatigue test.
2.6. Fatigue Recognition Model
The fatigue recognition model classifies the sEMG signal of muscle fatigue status and muscle non-fatigue status and was constructed based on LSTM networks. The proposed method is presented in
Figure 5. The LSTM network is a modified version of recurrent neural network (RNN). The control system of the LSTM unit consists of input, output and forget gates. The internal state
of the LSTM network records the historical information up to the current moment, and the three gates control the information transmission path [
31]. The forget gate
controls the amount of information that needs to be forgotten in the internal state
at the previous moment, and the activation state of the forget gate is computed as shown in Equation (17).
where
is the current input vector of the LSTM unit,
is the output of the previous LSTM unit, σ is the logistic sigmoid function,
is the weight vector, and
is the biased vector. The input gate
determines the amount of information stored in the current candidate state, and the calculation method is shown in Equations (18) and (19).
The data and information of the current LSTM unit are conveyed to the output gate, and the output calculations are shown in Equations (20) and (21).
The LSTM network used in this study consists of an LSTM layer, fully connected layer, rectified linear unit (ReLU) layer, dropout layer, and softmax layer. The fully connected layer connects all neurons of the previous and next layers. The dropout layer can make a certain neuron activation value stop working with a certain probability, which can make the model more generalized and prevent overfitting. Finally, the activation function of the softmax layer was used to classify the muscle state. The LSTM network used the initial hyperparameter configuration shown in
Table 1. We selected the Stochastic Momentum Gradient Descent (SGDM) algorithm to optimize the learnable parameters. During the LSTM network training process, these hyperparameters were reset many times until the optimal configuration was reached.
We constructed muscle fatigue recognition models based on SVM (support vector machine) and CNN (convolutional neural network) respectively and then compared the LSTM muscle fatigue recognition models. The kernel function has an important impact on the classification performance of the SVM. Numerous applications indicate that the Gaussian kernel function has good learning capability, so the SVM used in this study uses the Gaussian kernel function. The CNN network consists of 5 layers: an input layer, two convolutional layers, and two fully connected layers. The cross-entropy loss function is used to conduct back-propagation training based on the error. The initial learning rate is 0.1.
2.7. Evaluation of the Proposed Model
In order to ensure the generalization ability of the model, the muscle fatigue recognition model needs to be trained, verified, and tested on an independent data set [
32]. The verification method selected in this article is Holdout. This method consists of dividing the data into three independents subsets: training, validation, and test. The premise is that the ratio of fatigue data to non-fatigued data in each subset is approximately equal to the ratio of fatigue data to non-fatigued data in the total sEMG data set. The training data set was used for LSTM network training. The validation data set was used to evaluate performance during training through accurate measurements and error. The test data set was used for a final evaluation of the predictions performed by the model. Three datasets were independent, as
Figure 6 shows. The dataset division mode was divided into 3 types: (a) 70%, 10%, 20% (training, validation, testing); (b) 60%, 10%, 30% (training, validation, testing); (c) 50%, 10%, 40% (training, validation, testing).
We select four indicators of accuracy (
), sensitivity (
), specificity (
), and precision (
) to evaluate the performance of the model. The calculation methods of the above indicators are shown in Equations (22)–(25), respectively.
The parameters TP, TN, FF, and FN in Equations (22)–(25) are true fatigue (TF), true non-fatigue (TN), false fatigue (FF), and false non-fatigue (FN).
4. Discussion
In this section, the proposed methods and common methods will be discussed in terms of denoising performance and classification performance. The original sEMG signal was almost submerged by noise after adding Gaussian noise. The sEMG signal denoised by the hard threshold and soft threshold function still retains considerable noises. By comparing the images (
Figure 7 and
Figure 8) and denoising performance (
Table 2), it can be seen that the improved threshold function denoising method is more suitable for surface EMG signal denoising than the other two threshold functions. The signal denoised by the improved threshold function has better performance than the signal denoised by the hard threshold and soft threshold function in approximating the original signal.
The classification performancea of these methods are given in
Table 4 Wu et al. [
13] proposed a novel bacterial foraging algorithm (BFA)–Gaussian support vector classifier machine (GSVCM) model to improve the muscle fatigue classification accuracy. The RMS, IEMG, MPF, MF, and mean instantaneous frequency (MIF) features of the sEMG signals were extracted to evaluate the fatigue status during muscle contraction. With the training–testing rate of 80% 20%, the GSVCM model achieves the best accuracy of 93.94%. Khan et al. [
33] constructed three separate random forest models to classify muscle fatigue status based on eight sEMG features, with the best accuracy 87%. A CNN-SVM algorithm was proposed to identify muscle fatigue status [
33]. The RMS, IEMG, MPF, MF, and BSE features of the sEMG signals are extracted as the input dataset of the algorithm, and the best recognition accuracy is 86%.
This research proposes a new muscle fatigue recognition model based on LSTM, which is learned from scratch. In the preprocessing stage, the collected sEMG signal was processed by the wavelet packet threshold function denoising algorithm, and then the time-domain and frequency-domain features of the signal were extracted as the input of the LSTM algorithm. As shown in
Table 5, the wavelet packet threshold function denoising algorithm improves the performance of the classification algorithm, and the improved wavelet threshold function denoising algorithm has the best effect. With the same dataset as input, the LSTM algorithm performs better than CNN, SVM, and BFA–GSVCM in the accuracy of muscle fatigue recognition. The model proposed in this paper was trained, validated, and tested based on three different dataset ratios. The best accuracy of the model was obtained at type (a) 70%, 10%, 20% training, validation, testing. This indicates that the classification performance of the proposed model will improve with the increase in sEMG training data. However, the computation time for training process of the new muscle fatigue classification model is relatively long.
In addition, the model can be applied to the automatic detection of leg muscle fatigue during exercises other than cyclic resistance exercises. This model can be used as an auxiliary tool for coaches to monitor the muscle state of cyclists and long-distance runners during training. In the future, this model can be used to detect the state of muscles in other parts of the human body during exercise.
5. Conclusions
In this study, we analyzed the shortcomings of the traditional threshold function in the denoising of sEMG signal and proposed an improved threshold function. We compared the classification performance of LSTM with CNN, SVM, and other classification algorithms. In conclusion, the muscle fatigue recognition model constructed based on the improved wavelet packet threshold function denoising algorithm and LSTM network has excellent performance in the denoising of sEMG signal and the classification of muscle fatigue. The experimental results prove that the improved wavelet packet threshold function denoising algorithm is significantly better than the hard threshold and soft threshold functions in denoising the EMG signal. When the dataset was divided into training–verification–test 70% 10% 20%, the LSTM network achieved the best performance, including accuracy, sensitivity, specificity, and precision. Compared with other classification algorithms, the LSTM network achieved the best performance, including accuracy, sensitivity, and precision. However, we plan to improve the accuracy of the model in terms of feature extraction and algorithm optimization.