A Study on the Relationship between RPE and sEMG in Dynamic Contraction Based on the GPR Method

Abstract

The rating of perceived exertion (RPE) and surface electromyography (sEMG) describe exercise intensity subjectively and objectively, while there has been a lack of research on the relationship between them during dynamic contractions to predict exercise intensity, comprehensively. The purpose of this study was to establish a model of the relationship between sEMG and RPE during dynamic exercises. Therefore, 20 healthy male subjects were organized to perform an incremental load test on a cycle ergometer, and the subjects’ RPEs (Borg Scale 6–20) were collected every minute. Additionally, the sEMGs of the subjects’ eight lower limb muscles were collected. The sEMG features based on time domain, frequency domain and time–frequency domain methods were extracted, and the relationship model was established using Gaussian process regression (GPR). The results show that the sEMG and RPE of the selected lower limb muscles are significantly correlated (p < 0.05) but that they have different monotonic correlation degrees. The model that was established with all three domain features displayed optimal performance and when the RPE was 13, the prediction error was the smallest. The study is significant for lower limb muscle training strategy and quantification of training intensity from both subjective and objective aspects, and lays a foundation for sEMG further applications in rehabilitation medicine and sports training.

1. Introduction

Exercise intensity is a crucial factor in sports training. Exercising at intensities that are too low will not lead to the desired effects, while over-intensity may result in local muscular injuries [1,2]. Hence, exercise intensity should be monitored effectively through the use of scientific training methods and guidance contribute in order to optimize intensity and to avoid injuries caused by exercising at intensity levels that are too high as much as possible [3].

At present, scholars generally evaluate exercise intensity using measurable mechanical indicators outside of the body (including speed, power and force) as well as with interior indicators, which can be divided into objective indicators (such as surface electromyography (sEMG), heart rate, maximum oxygen uptake and blood lactate [4]) and subjective ones (such as the rating of perceived exertion (RPE) [3,5,6,7]). Subjectively, the RPE represents exercise intensity, or, more specifically, comprehensive feedback regarding physiological–psychological intensity [8].

Some scholars are committed to establishing relationships between RPE and objective indicators [3,9,10,11] during various exercises. Compared with the objective indicators that emphasize the evaluation of participants’ general exercise intensity [12], sEMG indicates the exercise intensity of the local muscles [13]. sEMG reflects the weighted sum of the electrical activity of the active motor units (Mus) in a specific type of exercise. Mu is the minimum motor functional unit of the body, so the sEMG contains information such as muscle fiber type, size and neural strategy [14]. Originally, RPE was used to evaluate general intensity, and later, other studies applied it to the evaluation of local muscles. Mosher et al. measured the overall and the local RPE in order to explore the effects of nitrate intake on resistance training [15]. Mathew et al. measured the working muscles and the overall RPE of the subjects in their to compare the cardiorespiratory and perceptual responses to exercises using self-regulated and imposed power outputs distributed between the arms and legs [16]. Moreover, through an intermittent grip strength experiment, Guo et al. discovered that RPE is not only able to indicate central fatigue, but that it can also indicate peripheral local muscle fatigue [17]. It can be seen that both RPE and sEMG are able to evaluate the local muscles, and therefore, it is of great significance to establish the relationship between RPE and sEMG, which would contribute to correct and integral instructions for the exercise intensity required for local muscles as well as in personalized exercise programs.

Although sEMG contains abundant body information, there is usually no explicit relationship between the muscle activity state, neural strategy and sEMG, so it is necessary to use appropriate data analysis methods to extract features that can better reflect the above information from sEMG. At present, sEMG feature extraction methods can be divided into three categories: time domain methods, frequency domain methods and time–frequency domain methods [18]. The features of time domain methods mainly include the root mean square value (RMS), the integrated electromyogram value (IEMG) and the zero-crossing rate (ZC). The time domain features have strong real-time performance, and amplitude-based time domain features are usually used to quantify muscle activation and the strength level [14]. The features of frequency domain methods mainly include the median (MF) and mean frequency of the power spectrum (MNF). Previous studies have shown that muscle fiber fatigue will be accompanied by a decrease in MF and MNF [3,19]. The features of the time–frequency domain methods include the wavelet transform (WT) method and short-time Fourier transform (STFT) method. Karthick et al. found that several time–frequency methods, such as wavelet analysis, can handle nonstationary and multi-component sEMG, and that time–frequency domain features are able to clearly identify muscle fatigue status during dynamic contractions [20]. However, using time–frequency domain methods for feature extraction has the disadvantage of high algorithm complexity.

There are some drawbacks indicated in existing studies on the relationship between RPE and sEMG in dynamic contractions. When establishing the relationship between sEMG and RPE while running, Chai Guozhong et al. only extracted two indicators: the RMS and MNF [3]. During the running process, the measured muscles were either stretched or contracted, and the joint angles and muscle power changed rapidly during recruitment de-recruitment of Mus, which are unable to stabilize the sEMG [21]. Therefore, traditional time domain and the frequency domain methods are not applicable for sEMG feature extraction, and a more complex method that is suitable for non-stable signals is needed, such as a time–frequency domain method.

Due to the complex relationship between RPE and sEMG, it is difficult to establish an accurate prediction model, so a data-driven approach is needed to establish the model. The Gaussian process regression (GPR) method is able to obtain a prediction model by learning historical samples and by inputting test samples into the model to obtain the prediction results, which have probabilistic significance [22]. Compared with support vector machines (SVM) and other methods, GPR is more suitable for nonlinear systems with higher dimensions, and the hyperparameters of the kernel function are easier to obtain.

The purpose of the study was to establish the relationship between RPE and sEMG based on cycle ergometer exercise. It was assumed that there would be some consistency between RPE and sEMG during dynamic exercises and that there would be some quantitative means that could be used to describe the relationship. This type of exercise is dynamic, so three methods, time domain, frequency domain and time–frequency domain methods were used. GPR was adopted in machine learning to establish the relationship model because it can overcome the noise that is present in neuromuscular measurements [23] and copes well with sEMG’s nonlinearity and high-complexity [22]. This study is important for the accurate judgement of exercise intensity and is potentially applicable for instruction during cycling-based exercise.

2. Materials and Methods

A. Subjects

A total of 20 healthy male skiers (height: 176.0 ± 4.8 cm; weight: 73.4 ± 9.1 kg; age: 19.9 ± 2.0 years) are selected to take part in the incremental load experiment on power bicycles. According to a survey that the participants completed regarding their exercise habits, the participants indicated that they generally exercise three times a week (non-cycling-related exercises). All the skiers did not have any neuromuscular diseases, plastic surgery history or any physical conditions that might affect the movement of their lower limbs. They were required not to train hard, consume alcohol or caffeine within 48 h before the experiment, and had at least 7 h of sleep before the experiment to ensure they were in good physical condition and their muscles were not fatigued and in good physical condition. The Ethics Committee of Physical Education College of Jilin University approved this experiment, which was conducted in compliance with the Declaration of Helsinki, and each subject submitted written informed consent before the experiment.

B. Procedures

The dominant leg of each subject was determined by inquiring about which leg the participant preferred to use to kick when playing soccer [24]. We made preparations to collect sEMG signals from their dominant legs. First, in order to reduce skin impedance and to ensure the accuracy of the measurement signal, we shaved the body hair from the subjects’ legs and used medical alcohol cotton balls to clean the newly exposed skin [25]. Due to the strong influence that innervation zones have on sEMG and to improve the standards and repeatability of the signals, we referred to the optimal electrode locations [26] chosen by Marco et al. Trigno Avanti Sensor electrodes (DELSYSY, Massachusetts, United States) were attached to the following sites: 20 mm on the anatomical landmark frames (ALF) vastus medialis (VM), 30% of the ALF rectus femoris (RF); 20 mm on the AFL vastus lateralis (VL); 80% of the ALF semitendinosus (S); 20% of the ALF biceps femoris (BF); 80% of the ALF tibialis anterior (TA); 80% of the ALF gastrocnemius lateralis (GL); 90% of the ALF gastrocnemius medialis (GM) [27].

Table 1. Muscles and the electrode locations. VM represents vastus medialis. RF represents rectus femoris. VL represents vastus lateralis. S represents semitendinosus. BF represents biceps femoris. TA represents tibialis anterior. GL represents gastrocnemius lateralis. GM represents gastrocnemius medialis.

Name of the Target Muscle	VM, RF, VL	S, BF	TA	GL, GM
Electrode location

shows the specific positions of each of the electrodes. In order to reduce the influence of motion artifact noise and to not affect the subjects’ normal movements, medical adhesive tapes were used to fix the electrodes to the skin. At the same time, the subjects’ heart rates were monitored with a heart rate belt (Polar H9, Polar Electro Oy, Kempele, Finland), which was worn on the chest across the nipples, with the sensor placed between the two nipples and closer to the left nipple.

The subjects were then asked to sit on the cycle ergometer (Monark 839 e, Vansbro, Sweden), and the seat and grip were adjusted to the proper height in order to carry out step incremental load test. The experiment consisted of a 5 min warm-up at a zero-load pace of 60 rpm, during which the participants became familiar with the equipment and listened to the instructions for Borg RPE Scale (6–20) test read by the professional conducting the test. The subjects described their self-perceived strenuousness and tiredness using the numbers from 6 to 15, with large numbers representing higher levels of fatigue [28].

After the warm-up, the initial load was 50 W, and the load increased by 50 W every 3 min. In the last 10 s of every minute, the RPE values of the lower limb muscles were reported. The experiment stopped when the subjects could not reach the prescribed rhythm for 5 s, or, in order to protect their health, when their heart rate exceeded 90% of the subject’s expected maximum heart rate (220 minus age) or RPE ≥ 18. Using the Trigno Avanti Platform (DELSYSY, Natick, MA, USA), the sEMGs of the target muscles were measured at a frequency of 1926 Hz.

C. Signal Processing and Data Analysis

First, the sEMGs that were obtained were processed visually to ensure that there was no manual intervention or power wire interference, as shown in Figure 1. The sEMGs were processed using MATLAB 2020b (Math Works, Inc., Natick, MA, USA). Then, in order to reduce the impact of noise, Butterworth bandpass filter (20–500 Hz) was used to filter the signal. The subjects were required to ride at 60 r/min with a contraction cycle of 1 s. Providing instructions in the moment ensured that muscle contraction occurred within the first half second. The signals were intercepted into smaller segments of 3 s. It is assumed that the contractions in the segment were almost the same, and they were processed evenly as one contraction in order to reduce the data size. During this time period, the signal could be approximately regarded as stationary signal. In each segment, time-domain features, frequency domain features based on Fourier transform and time–frequency domain features were extracted from the signal. These features are: IEMG, RMS and ZC; MF and MNF; the RMS (WT_RMS), MF (WT_MF) and MNF (WT_MNF) of the effective frequency band wavelet coefficients and MF (WT_MF(HF)) and MNF (WT_MNF(HF)) of the high-frequency portion of the wavelet coefficients; the RMS (STFT_RMS), MF (STFT_MF) and MNF (STFT_MNF) of the effective frequency band based on STFT; the and MF (STFT_MF(HF)) and MNF (STFT_MNF(HF)) of the high-frequency portion based on STFT.

The IEMG, RMS, and ZC can be expressed as follows:

$$\mathrm{IEMG}=\frac{1}{\mathrm{N}}{\displaystyle \sum }_{\mathrm{i}=1}^{\mathrm{N}}\left|{\mathrm{x}}_{\mathrm{i}}\right|,$$

(1)

$$\mathrm{RMS}=\sqrt{\frac{1}{\mathrm{N}}{\displaystyle \sum }_{\mathrm{i}=1}^{\mathrm{N}}{\mathrm{x}}_{\mathrm{i}}^2},$$

(2)

$$\mathrm{ZC}={\displaystyle \sum }_{\mathrm{i}=1}^{\mathrm{N}-1}\left[\mathrm{sgn}\left({\mathrm{x}}_{\mathrm{i}}{\times \mathrm{x}}_{\mathrm{i}+1}\right)\cap \left|{\mathrm{x}}_{\mathrm{i}}{-\mathrm{x}}_{\mathrm{i}+1}\right|\ge 0\right],$$

(3)

$$\mathrm{sgn}\left(\mathrm{x}\right)=\{\begin{array}{c}1, \mathrm{if} \mathrm{x}\ge 0\\ 0, \mathrm{otherwise}\end{array},$$

(4)

where ${\mathrm{x}}_{\mathrm{i}}$ represents the sEMG value at sampling point i and N represents the total number of sEMG signal sampling points.

MF and MNF can be expressed as follows:

$${{\displaystyle \int }}_0^{\mathrm{MF}}\mathrm{p}\left(\mathrm{f}\right)\mathrm{df}={{\displaystyle \int }}_{\mathrm{MF}}^{\infty }\mathrm{p}\left(\mathrm{f}\right)\mathrm{df}=\frac{1}{2}{{\displaystyle \int }}_0^{\infty }\mathrm{p}\left(\mathrm{f}\right)\mathrm{df},$$

(5)

$$\mathrm{MNF}={{\displaystyle \int }}_0^{\infty }\frac{\mathrm{fp}\left(\mathrm{f}\right)\mathrm{df}}{{{\displaystyle \int }}_0^{\infty }\mathrm{p}\left(\mathrm{f}\right)\mathrm{df}},$$

(6)

where f represents the frequency and p(f) represents the value of the signal power spectrum at frequency f.

To obtain the features of the wavelet coefficient sEMG, we first took the sym6 function as the parent wavelet to carry out discrete wavelet transformation for sEMG [29]. Since the available sEMG energy is distributed in the range of 20–500 Hz and mainly in the range of 50–150 Hz and because the time scale calculation needs to be simplified to avoid the omission of useful information, a wavelet coefficient that was calculated from the effective frequency band was selected to calculate WT_RMS, WT_MF and WT_MNF. Additionally, WT_MF (HF) and WT_MNF (HF) were also calculated as sEMG features in the high frequency section (45–250 Hz) [30].

Similarly, in order to obtain the sEMG features based on STFT, we performed STFT on sEMG and calculated each feature in the range of 20–250 Hz and 45–250 Hz.

In order to smooth the sEMG features, a sliding window filtering algorithm (10 s data window and 5 s sliding window) was used to weaken the influence of data jitter [31].

To minimize the impact of individual differences, all of the data were linearly normalized.

D. Statistical Analysis

SPSS Statistics 21.0.0 (IBM, Armonk, NY, USA) was used to conduct the statistical analysis. After the evaluation of data distribution using the Kolmogorov–Smirnov test, the Spearman coefficient was used to analyze the monotone correlation between the RPE and sEMG features of all of the muscles. Statistical significance determined at p < 0.05.

E. Machine Learning

Considering that sEMG is a highly complex and nonlinear signal with a small amount of data [22], the GPR method was used to establish the regression model between sEMG features and RPE. GPR is a non-parametric probabilistic model that is the kernel function based [32] and that does not seek to build a certain type of optimal fitting model (such as linear or exponential form) between data input and output, and instead calculates the posterior test distribution of new input data to predict output data. The advantage of this method is that it fits well nonlinear functions, meaning that it is beneficial for small datasets [33].

Gaussian process (GP) can be regarded as a set of random variables. Each random variable is defined by $\mathrm{f}\left(\mathrm{x}\right)(\mathrm{x}\in \mathrm{X})$ its mean function $\mathsf{\mu }\left(\mathrm{x}\right)$ and covariance function k(x,x’), which can be expressed as:

$$\left[\begin{array}{c}\mathrm{f}\left({\mathrm{x}}_1\right)\\ ⋮\\ \mathrm{f}\left({\mathrm{x}}_{\mathrm{n}}\right)\end{array}\right]~\mathrm{N}\left(\begin{array}{cc}\left[\begin{array}{c}\mathsf{\mu }\left({\mathrm{x}}_1\right)\\ ⋮\\ \mathsf{\mu }\left({\mathrm{x}}_{\mathrm{n}}\right)\end{array}\right]& \left[\begin{array}{ccc}\mathrm{k}\left({\mathrm{x}}_1{,\mathrm{x}}_1\right)& \cdots & \mathrm{k}\left({\mathrm{x}}_1{,\mathrm{x}}_{\mathrm{n}}\right)\\ ⋮& \ddots & ⋮\\ \mathrm{k}\left({\mathrm{x}}_{\mathrm{n}}{,\mathrm{x}}_1\right)& \cdots & \mathrm{k}\left({\mathrm{x}}_{\mathrm{n}}{,\mathrm{x}}_{\mathrm{n}}\right)\end{array}\right]\end{array}\right).$$

(7)

It can also be written as:

$$\mathrm{f}\left(\mathrm{x}\right)~\mathrm{GP}\left(\mathsf{\mu }\left(\mathrm{x}\right),\mathrm{k}\left(\mathrm{x},\mathrm{x}’\right)\right),$$

(8)

where each dimension of the Gaussian function corresponds to an element x of the index set X and the respective components of the random vector represent the values of $\mathrm{f}\left(\mathrm{x}\right)$.

MATLAB was used for all of the GP modeling. We took the average sEMG features from each minute as the independent variables ${\mathrm{x}}_{\mathrm{i}}$, measured the RPE value during this period as the dependent variable ${\mathrm{y}}_{\mathrm{i}}$ and established the data set ${\left\{\left({\mathrm{x}}_{\mathrm{i}},{\mathrm{y}}_{\mathrm{i}}\right)\right\}}_{\mathrm{i}=1}^n$, where ${\mathrm{x}}_{\mathrm{i}}\in {\mathrm{R}}^{\mathrm{d}},{\mathrm{y}}_{\mathrm{i}}\in \mathrm{R}$. The response of $\mathrm{y}$ can be expressed as [34] follows:

$$\mathrm{y}=\mathrm{f}\left(\mathrm{x}\right)+\mathsf{\epsilon },$$

(9)

where f(x) is a GP with a mean of 0 and a covariance of $\mathrm{k}(\mathrm{x},\mathrm{x}’)$ and ε is Gaussian noise.

If the output prediction is made for the new input variable Z (z = [z₁,…,z_m]^T ), based on the properties of Gaussian distribution, the joint distribution of training points and prediction points is as follows:

$$\left[\begin{array}{c}\mathrm{y}\\ \mathrm{f}\left(\mathrm{Z}\right)\end{array}\right]~\mathrm{N}\left(\begin{array}{cc}\left[\begin{array}{c}\mathsf{\mu }\left(\mathrm{X}\right)\\ \mathsf{\mu }\left(\mathrm{Z}\right)\end{array}\right]& \left[\begin{array}{cc}\mathrm{K}\left(\mathrm{X},\mathrm{X}\right){+\sigma }_{\mathrm{n}}^2\mathrm{I}& \mathrm{K}{\left(\mathrm{Z},\mathrm{X}\right)}^{\mathrm{T}}\\ \mathrm{K}\left(\mathrm{Z},\mathrm{X}\right)& \mathrm{K}\left(\mathrm{Z},\mathrm{Z}\right)\end{array}\right]\end{array}\right),$$

(10)

where $\mathrm{y}= {\left[{\mathrm{y}}_1,\cdots {,\mathrm{y}}_{\mathrm{n}}\right]}^{\mathrm{T}}$, $\mathsf{\mu }\left(\mathrm{X}\right)$ = ${\left[{\mathsf{\mu }(\mathrm{x}}_1),\cdots {,\mathsf{\mu }(\mathrm{x}}_{\mathrm{n}})\right]}^{\mathrm{T}}$ and $\mathsf{\mu }\left(\mathrm{Z}\right)$ = ${\left[{\mathsf{\mu }(\mathrm{z}}_1),\cdots {,\mathsf{\mu }(\mathrm{z}}_{\mathrm{m}})\right]}^{\mathrm{T}}$$,\mathrm{K}\left(\mathrm{X},\mathrm{X}\right)$ is the $\mathrm{n}\times \mathrm{n}$ matrix, in which the $(\mathrm{i},\mathrm{j})$ element ${[\mathrm{K}\left(\mathrm{X},\mathrm{X}\right)]}_{\mathrm{ij}}{ =\mathrm{k}(\mathrm{x}}_{\mathrm{i}}{,\mathrm{x}}_{\mathrm{j}})$, $\mathrm{K}\left(\mathrm{Z},\mathrm{X}\right)$ is the matrix $\mathrm{m}\times \mathrm{n}$ and ${[\mathrm{K}\left(\mathrm{Z},\mathrm{X}\right)]}_{\mathrm{ij}}{ =\mathrm{k}(\mathrm{z}}_{\mathrm{i}}{,\mathrm{x}}_{\mathrm{j}})$. $\mathrm{K}\left(\mathrm{Z},\mathrm{Z}\right)$ is the matrix $\mathrm{m}\times \mathrm{m}$, and ${[\mathrm{K}\left(\mathrm{Z},\mathrm{Z}\right)]}_{\mathrm{ij}}{ =\mathrm{k}(\mathrm{z}}_{\mathrm{i}}{,\mathrm{z}}_{\mathrm{j}})$.

Using the conditional distribution property of Gaussian distribution, the following can be deduced:

$$\mathrm{y}(\mathrm{Z})|\mathrm{y},\mathrm{X},\mathrm{Z}~\mathrm{N}(\hat{\mathrm{u}}\hat{\sum }+{\sigma }^2\mathrm{I}),$$

(11)

where the mean function and covariance function are defined as:

$$\hat{\mathrm{u}}=\mathrm{K}{(\mathrm{X},\mathrm{Z})}^{\mathrm{T}}{\left(\mathrm{K}\right(\mathrm{X},\mathrm{X})+{\sigma }^2\mathrm{I})}^{-1}\mathrm{y},$$

(12)

$$\hat{\mathrm{\sum }}=\mathrm{K}(\mathrm{Z},\mathrm{Z})-\mathrm{K}{(\mathrm{X},\mathrm{Z})}^{\mathrm{T}}{\left(\mathrm{K}\right(\mathrm{X},\mathrm{X})+{\sigma }^2\mathrm{I})}^{-1}\mathrm{K}(\mathrm{X},\mathrm{Z})$$

(13)

In addition, the definition of y(Z) is similar to y.

From this, any new input can be derived as the average of the posterior prediction distribution.

The selection of the covariance function (kernel function) is generally based on the assumption of the main function to be modeled. In this study, we chose the exponential kernel function and the rational quadratic kernel function. The exponential kernel function is closely related to the Gaussian kernel function, with the exception of the square being left out of the norm. It is also a radial basis function kernel. It can be described as follows:

$$\mathrm{k}\left(\mathrm{x},\mathrm{x}’\right){=\mathrm{e}}^{-\frac{||\mathrm{x}-\mathrm{x}’||}{{2\sigma }^2}},$$

(14)

where ${\sigma }^2$ represents the variance.

The rational quadratic kernel function is less computationally intensive than the Gaussian kernel function and can be used as an alternative when using the Gaussian becomes too expensive. It can be described as:

$$\mathrm{k}\left(\mathrm{x},\mathrm{x}’\right)=1-\frac{{||\mathrm{x}-\mathrm{x}’||}^2}{{||\mathrm{x}-\mathrm{x}’||}^2+\mathrm{c}},$$

(15)

where c represents an optional constant.

Then, 20% of the data set were selected randomly as the test set, k-fold cross-validation (k = 5) and leave one out validation to evaluate the generalization ability of the model, root mean square error (RMSE) and R-squared (R²) as the evaluation index of the model, and input the remaining 80% into the model for training.

3. Results

There was a significant correlation (p < 0.05) between 15 sEMG from all of the eight lower limb muscles and RPE, but each demonstrated different monotone correlations. Not considering the positive and the negative directions, we determined that the RPE and sEMG features are strongly correlated when the absolute value of Spearman’s rank correlation coefficient (represented by |R| later) is more than 0.7; when 0.4 < |R| < 0.7, a moderate correlation is observed; when 0.1 < |R| < 0.4, a weak correlation is observed; when |R| < 0.1, a negligible correlation is observed [35]. As shown in Figure 2, the 15 different sEMG features of GM showed a moderate monotonic correlation with RPE. The monotone correlation between the time domain features and RPE was stronger than the frequency domain features. The degree of the monotone correlation between the RMS and RPE of the eight muscles was similar to that of the RMS and RPE of the time–frequency features. The monotone correlation between the MF and MNF of eight muscles, and the RPE was slightly weaker than that of the MF and MNF in the effective frequency band and in the high-frequency time–frequency features. In general, the monotone correlation between the features obtained by the two time–frequency domain methods and RPE is very similar.

According to Spearman’s analysis, the sEMG features of some of the selected lower limb muscles had strong correlations with RPE (

Table 2. sEMG features and the corresponding muscles show a strong correlation with RPE. 95% CI represents the range of the corresponding correlation coefficient when p < 0.05; ^# Subjects p < 0.05 refers to the probability of the number of samples that meet the condition (p < 0.05).

Feature	Muscle	|R|	95%CI	# Subjects p < 0.05
IEMG	RF	0.85	0.41~0.98	20/20
IEMG	VM	0.92	0.77~0.96	20/20
IEMG	VL	0.90	0.62~0.97	20/20
IEMG	BF	0.90	0.74~0.97	20/20
IEMG	S	0.86	0.56~0.96	20/20
IEMG	GL	0.73	0.12~0.95	19/20
IEMG	TA	0.76	0.29~0.94	20/20
RMS	RF	0.82	0.32~0.97	20/20
RMS	VM	0.90	0.75~0.96	20/20
RMS	VL	0.88	0.54~0.96	20/20
RMS	BF	0.87	0.67~0.97	20/20
RMS	S	0.83	0.47~0.96	20/20
RMS	TA	0.72	0.29~0.91	20/20
ZC	VM	0.75	0.38~0.94	20/20
ZC	VL	0.72	0.39~0.95	20/20
ZC	BF	0.81	0.42~0.94	19/20
ZC	S	0.80	0.39~0.94	20/20
WT_MNF(HF)	BF	0.71	0.21~0.92	19/20
WT_RMS	RF	0.84	0.22~0.97	20/20
WT_RMS	VM	0.91	0.8~0.96	20/20
WT_RMS	VL	0.89	0.65~0.96	20/20
WT_RMS	BF	0.89	0.8~0.97	20/20
WT_RMS	S	0.85	0.61~0.95	20/20
WT_RMS	TA	0.71	0.3~0.91	20/20
STFT_MF(HF)	BF	0.73	0.42~0.91	18/20
STFT_MNF(HF)	BF	0.72	0.14~0.92	19/20
STFT_RMS	RF	0.84	0.22~0.97	20/20
STFT_RMS	VM	0.90	0.79~0.96	20/20
STFT_RMS	VL	0.88	0.62~0.96	20/20
STFT_RMS	BF	0.88	0.71~0.97	20/20
STFT_RMS	S	0.84	0.61~0.95	20/20
STFT_RMS	TA	0.71	0.31~0.93	20/20

). The IEMG of seven muscles had strong correlations with the RPE. RMS, WT_RMS and STFT_RMS of six muscles had strong correlations with the RPE and ZC of four muscles had strong correlations with the RPE. All of the RMS, ZC, WT_RMS and STFT_RMS of the VM, VL, BF and S had strong correlations with the RPE. At the same time, for the same muscle, the |R| of the WT_RMS and STFT_RMS was a little more than that observed for the RMS. Among them, the IEMG of the VM had the strongest monotone correlation with the RPE.

We assessed the stability of all the sEMG features selected, and by plotting all the sEMG features of each muscle across all subjects, we found that the variation trend of the features of each subject was basically consistent and within a certain range (Figure 3). We believed that all the features we selected had certain stability.

In order to demonstrate the advantages of the GPR method during sEMG processing, we also adopted the SVM method to build a model of the relationship between sEMG and RPE. For the SVM model, the linear kernel function and cubic kernel function were selected for modeling, and two validation methods were used to verify the model (

Table 3. Shows regression performance of support vector machines (SVM) and Gaussian process regression (GPR) models with the different kernel functions. ^# represents that the goodness of fit value of the model comes from leave one out validation, and * represents that the goodness of fit value of the model comes from 5-fold cross validation. WT indicates that the model is built based on all of the time domain features, frequency domain features and wavelet features. The STFT representative model was built based on all of the time domain features, frequency domain features, and STFT features. R² represents R-squared and RMSE represents root mean squared error.

	RMSE #	R2 #	RMSE *	R2 *
Linear Kernel SVM (WT)	1.7205	0.63	1.6650	0.63
Cubic Kernel SVM (WT)	1.5653	0.70	1.5693	0.67
Rational Quadratic Kernel GPR (WT)	1.5411	0.71	1.4324	0.73
Exponential Kernel GPR (WT)	1.4635	0.74	1.4285	0.73
Linear Kernel SVM (STFT)	1.7661	0.61	1.6879	0.62
Cubic Kernel SVM (STFT)	1.6239	0.67	1.5288	0.69
Rational Quadratic Kernel GPR (STFT)	1.4644	0.74	1.4103	0.74
Exponential Kernel GPR (STFT)	1.4830	0.73	1.4126	0.75

). The results show that there is a large gap between the goodness of fit values obtained by leave one out validation and 5-fold cross validation, and the model fitting effect based on 5-fold cross validation is better. Both model validation methods show that the performance of the GPR model is better than that of the SVM model. The GPR model constructed by STFT in the time–frequency domain is superior to that constructed by WT. In general, the GPR model constructed based on the rational quadratic kernel function is optimal when the STFT features are selected from the time–frequency domain features.

Apart from basic evaluation of regression performance using metrics including R² and RMSE, further evaluation of classification performance is added for each model. Raw data with Borg Scale < 12 is classified as negative, the rest is positive, and all regression results are rounded up to the nearest integer, as the pre-process of evaluation. The results show that there is little difference between classification performance of different models by two validation methods (

Table 4. Shows classification performance of SVM and GPR models with the different kernel functions.

	Accuracy	Kappa Score	Sensitivity	Specificity
(a) The Performance of the Model Comes from Leave One Out Validation
Linear Kernel SVM (WT)	0.846154	0.828979	0.942149	0.839324
Cubic Kernel SVM (WT)	0.833333	0.815837	0.939394	0.825219
Rational Quadratic Kernel GPR (WT)	0.846154	0.828791	0.939394	0.839295
Exponential Kernel GPR (WT)	0.846154	0.828697	0.922865	0.841335
Linear Kernel SVM (STFT)	0.846154	0.829104	0.922865	0.841335
Cubic Kernel SVM (STFT)	0.858974	0.844142	0.922865	0.855411
Rational Quadratic Kernel GPR (STFT)	0.884615	0.871499	0.936639	0.881497
Exponential Kernel GPR (STFT)	0.884615	0.871734	0.936639	0.881497
(b) The Performance of the Model Comes from 5-Fold Cross Validation
Linear Kernel SVM (WT)	0.846154	0.828916	0.942149	0.839324
Cubic Kernel SVM (WT)	0.846154	0.829756	0.955923	0.837256
Rational Quadratic Kernel GPR (WT)	0.846154	0.828791	0.939394	0.839295
Exponential Kernel GPR (WT)	0.846154	0.828697	0.922865	0.841335
Linear Kernel SVM (STFT)	0.846154	0.829135	0.922865	0.841335
Cubic Kernel SVM (STFT)	0.846154	0.829880	0.906336	0.843374
Rational Quadratic Kernel GPR (STFT)	0.884615	0.871499	0.936639	0.881497
Exponential Kernel GPR (STFT)	0.884615	0.871734	0.936639	0.881497

). Both model validation methods show that the performance of the GPR model is better than that of the SVM model. The GPR model constructed by STFT in the time–frequency domain is superior to that constructed by WT in the time–frequency domain. In general, the GPR model constructed based on the exponential kernel function is optimal when the STFT features are selected from the time–frequency domain features.

The fitting degree was the best when time domain features (IEMG, RMS and ZC), frequency domain features (MF and MNF) and high-frequency portion features based on STFT (STFT_MF(HF) and STFT_MNF(HF)) were used for modeling and when the RMSE was 1.4029 and R2 was 0.74 (Figure 4). The performance was better with the model that used all of the time domain features than the one using all of the time–frequency domain features (features based on WT and STFT) and was relatively poor with the one using only frequency domain features. The performance of the model using all of the time domain features was better than that of the one using only single time domain feature, and this rule also applied to the frequency domain features and the time–frequency domain features. The model based on STFT features is superior to that based on WT features. Whether to extract the high frequency proportion of time–frequency domain features for modeling had a certain influence on the model performance.

The prediction error was different in actual RPE measurement obtained when using the GPR models based on all of the time domain features, frequency domain features and STFT features (Figure 5), which functioned the best when the RPE measured as 13, and when further away from 13, a less desirable effect can be observed. At the same time, when the measured RPE was less than 13, the predicted value was relatively low, and when the measured RPE was more than 13, the predicted value was relatively high. Different kernel functions have little effect on model performance.

4. Discussion

The purpose of this study was to establish a quantitative relationship between RPE and sEMG through machine learning. We assumed that the relationship could be described by a certain quantitative method during dynamic exercise. Our research supports this assumption. The results show that the relationship can be described by the GPR model and that it achieves a desirable goodness of fit.

In this study, statistical analysis methods were used to analyze the correlation between sEMG features and RPE, and there were significant correlations between the RPE and the respective sEMG features of the selected lower limb muscles (p < 0.05). Therefore, we believe that during the step incremental load test, it is meaningful to establish relationships based on the VM, RF, VL, S, BF, TA, GL and GM through the use of the three sEMG pre-processing methods: time domain, frequency domain and the time–frequency domain methods. When cycling, coordination is achieved by muscle recruitment and via coordination in the lower limb tendon structures [36]. During a pedal cycle, different muscles are believed to make different functional contributions [37]. The results of this study show that despite the different contributions of the eight selected muscles during the exercise, all of their neuromuscular progressive changes are correlated with the RPE. Guo et al. discovered that the RPE of the subjects in their study was significantly correlated with the RMS of sEMG in the intermittent grip strength test, which is consistent with our results [17]. However, Chai et al. discovered that in two running exercises, the RMS of the BF and sEMG did not have any significant correlation with one another, nor did MNF and RPE [3]. We believe that the different exercises caused these differences.

It has been suggested that WT and STFT should be suitable for analyzing the nonstationary myoelectric signals that are recorded during dynamic contractions [38,39]. Our study confirms that because of the Spearman’s rank correlation coefficient between WT_RMS and STFT_RMS and RPE is higher than RMS (Table 2). Similarly, during the modeling process (Figure 4), using the time–frequency domain features were better than using the corresponding frequency domain indicators.

The Spearman correlation is a non-parametric measure of the monotonicity of the relationship between the two data sets. Since different the sEMG features of different muscles and RPE were not all strongly correlated, the monotonic correlation model cannot be used to determine the relationship. Therefore, machine learning model was used to set up a model for the RPE based on multiple sEMG features. During dynamic contractions, sEMG can be described as a complex, highly nonlinear and nonstationary signal. According to the principles of the two algorithms, GPR method is more suitable for processing nonlinear signals than SVM method. Using two different validation methods, we confirmed that GPR models were always performing better than SVM models in our study (Table 3 and Table 4).

The model was only able to show optimal performance when all of the three domain features were used for modeling (Figure 4). That is, all of the three domain features that were calculated included typical neuromuscular information during the step incremental load test and could be identified by the model established in this study. Furthermore, each feature carried distinctive information regarding its correlation with RPE. Plus, the study also shows that STFT features provide more efficient information than WT features in using sEMG features to predict RPE. We believe that WT has the characteristic of adaptive frequency resolution, and compared with STFT, it may be less efficient to extract effective information from sEMG in the adaptive process.

When the model shows a score of 13, which is “somewhat hard”, then the RPE is the most accurate, and when it is further away from 13, the error is larger (Figure 5). Lee et al. believed that the evaluation of the exercise intensity with RPE (6–20) was moderately reliable [40]. Yang et al. believed that there might be some errors during the RPE (6–20) quantification of athletes with a low resting heart rate and a high reserve heart rate [8]. Based on the experiments conducted above, we hold that model errors probably result from the improper scale division selected for RPE. In other words, when the subjects hardly feel fatigue (RPE ≤ 11) or exhausted (RPE > 16), it is difficult to determine RPE the value or to reflect their actual condition. Therefore, extremely high or low RPE value cannot reflect the real condition of the subjects’ fatigue. So, in the future, scholars can redefine the Borg scale for a specific sport to improve the accuracy of the RPE value.

There are some limitations associated with the present study. Only a small number of healthy young males were selected to participate in the study, so considering the gender and age differences in neuro-muscle and muscle metabolism [41], the model depicting the relationship between RPE and sEMG that was established in this study may not apply well to other gender or age groups. The relationship that was established in this study was determined based on the results of a cycling experiment, so people should be cautious when applying this model to other exercises.

5. Conclusions

In summary, this study established a relationship model between the sEMG and RPE of the lower limb muscles during dynamic contractions with the use of machine learning to predict exercise intensity. Compared with the present related studies, the model extracted the sEMG features using time–frequency domain methods, and it turns effective to provide muscle information. What is more, we adopted GPR algorithm to build the predictive model innovatively and the established model performs well, since GPR is suitable for handling nonstationary signals such as sEMG. The model can quantitatively analyze degree 6 to 18 in RPE, especially can accurately correlate degree 12 to 15 to their corresponding sEMG. The future research direction lies in how to design a more reasonable subjective scale to describe human exercise intensity and explore a better model and better sEMG features of quantifying exercise intensity. This study contributes to instructions provided during cycling exercises and will help coaches or rehabilitation physicians to control lower limb conditions accurately, enabling them adjust subsequent dynamic exercise strategies.