Time-frequency analysis of surface EMG signals for human movement identification

The assessment of muscle-recruitment timing from surface EMG signal (sEMG) is relevant in different fields, including clinical gait analysis and robotic systems to interpret user’s motion intention. However, available methods typically provide only information in time domain without evaluating muscle-activation frequency content. This study aims to propose a novel adaptative algorithm for detecting muscle activation in time-frequency domain based on continuous wavelet transform (CWT) analysis. Precisely, the novel contribution of the proposed algorithm consists of evaluating the frequency range of every muscle activations detected in time domain. Performances are evaluated on a test bench of 720 simulated sEMG signals and 105 real sEMG signals, stratified for signal-to-noise ratio (SNR), and then validated against different reference algorithms. Outcomes indicate that the proposed approach can provide an accurate prediction of muscle onset and offset timing in both simulated (mean absolute error, MAE ≈ 10 ms) and real datasets (MAE < 30 ms), minimally affected by the SNR variability and compatible with the timing of sEMG-driven assistive devices. Concomitantly, the maximum frequency of the activations is computed, ranging from around 100 Hz up to almost 500 Hz. This suggest a large within-muscle between-muscle variability of the frequency range. In conclusion, the current study introduces a novel reliable wavelet-based algorithm to detect both time and frequency content of muscle activation, suitable in different conditions of signal quality.


I. INTRODUCTION
Muscle activity is typically quantified during human movement using surface electromyography (sEMG). sEMG is a non-invasive approach that permits the attainment of many different tasks, such as estimating muscle function, evaluating muscular fatigue, identifying movement patterns, explaining neuromotor-system strategies in complex tasks, controlling sEMG-driven assistive devices, and so on [1][2][3]. Expressly, quantifying muscular activity from sEMG signals is acknowledged to be particularly relevant in clinical gait analysis [4] and robotic systems to interpret user's motion intention [3,5,6].
Most of the literature in the field is mainly focused on the analysis of the muscular activity in time domain to extract the timing of activation onset and offset. This approach is extensively adopted in clinics to quantitatively characterize different neuromotor disorders [7,8], patient follow-up, and improve rehabilitation strategies [9]. Furthermore, myoelectric non-pattern-recognition-based prosthetic control is acknowledged to include an activation-timing algorithm as one of the fundamental steps of the control procedure [10]. To this aim, the assessment of the onset event is more crucial than identifying the offset event. On the other hand, the approach of myoelectric pattern-recognitioncontrolled devices typically mimics the functionality of an actual limb, employing EMG signals from residual muscles leftover after amputation or congenital disability [11]. This approach presupposes a reliable acquisition and preprocessing of sEMG signals, a focused feature-extraction procedure, an accurate signal-classification process to predict a subset of designated motor tasks, and eventually a multifunction prosthesis control. Thus, the capability of these systems to decode the motor intention is determined by the quality of system performance at each of the phases mentioned above [12]. Implementing an onset-detection algorithm could help since it has been shown that the performance of the classification phase could be improved by windowing sEMG signals using a reliable assessment of activation events [13].
In the same way, it has been suggested that the onsetdetection-based interpretation of sEMG signal could be adopted together with or as a possible alternative to a more complex EEG-based approach to improve movement intention detection in stroke patients without residual movements [14]. More generally, onset-offset detection algorithms are used to develop increasingly sophisticated and reliable EMG-based human-machine interfaces [15]. Further applications of the assessment of activation onset in rehabilitation robotics are also reported [6,16].
The onset (and sometimes the offset) timing is one of the most frequently analyzed parameters in the characterization of sEMG signals. Different approaches are introduced to assess it. Visual identification is typically adopted to detect the onset and offset timings of muscle activation [17,18]. This approach consists of inspecting EMG signals by trained human-movement specialists who visually identify the onset event through a preset criterion. Unfortunately, this empirical procedure is very time-consuming, mainly when identifying activation events in a large dataset or long-lasting sEMG signals [19]. Moreover, it could be susceptible to inter-expert variability, specifically for the low value of signal-to-noise ratio (SNR) [20].
To lower procedure duration, limit inter-expert variability, and improve objectivity of the analysis, computer-based techniques are introduced. A classical category of automatic methods is the one representing the threshold-based algorithms [21,22], including the double-threshold (DT) statistical algorithm [23,24]. DT is even now considered among the most reliable detection techniques; DT is indeed broadly used both in research and in clinics, and it is also implemented in commercial systems [25,26]. A more recent approach adopts an adaptive threshold computed by the constant false alarm rate (CFAR) algorithm, where the EMG signal is refined by morphological hole filling used to close up and fill out missing information [27].
Despite the great effort in developing of techniques for EMG-signal processing in time domain, only a small number of studies are focused on sEMG analysis in frequency domain. Frequency analysis is mainly performed to compute the power spectrum of the sEMG signal by applying the Fourier transform, quantify the frequency content of the whole signal and/or provide information on the muscle fatigue process [28,29].
It is acknowledged that muscle fibers can be broken down into two main types: slow-twitch fibers, which respond slowly to stimulus but are resistant to fatigue, and fast-twitch fibers, which respond but also wear out more quickly. The different fibers are characterized by different sEMG-signal frequency content [30][31][32]. Thus, for a deeper comprehension of the electrophysiological processes behind neuromuscular activations, sEMG signals should be analyzed in time as well as frequency domain. This is successfully achieved, for example, by Lauer et al. [33]. They introduced a time-frequency index that is significantly correlated to gait kinetics, kinematics, and clinical measures of motor impairment in children with cerebral palsy and is sensitive to walking ability according to Gross Motor Functional Classification Scale. Another example of the significant contribution of time-frequency analysis is presented by Trigili et al. [34], who adopted a cognitive human-robot interface approach to detect the users' motion intention fast through a machine-learning interpretation of a set of subject-independent sEMG time-domain features.
Moreover, Sacco et al. [35] investigated the spectral properties and muscle energy patterns in diabetic subjects with neuropathy during gait using wavelet transform, showing that lower energy content could be found at a lower frequency (slow-twitch fibers) due to impairments caused by neuropathy. The time-frequency approach has also been adopted in complex classification tasks. Wang et al. [36] presented a wavelet-based correlation dimension approach to classify SEMG signals in a prosthetic control system. Recently, wavelet transform has become an efficient tool to assess pertinent information from the EMG signal [37]. As shown by both above-mentioned and other studies, indeed, wavelet transform analysis has been proven to be a reliable tool for the simultaneous characterization of EMG data in time and frequency domain, providing a suitable approach for the identification of motor-unit/fiber recruitment patterns [38], the quantification of motor strategy patterns [35], the project of rehabilitation devices [36], and the assessments of clinical features [33,39,40].
Thus, the present study aims to propose a novel adaptative algorithm for detecting muscle activation in time-frequency domain, based on continuous wavelet transform (CWT) analysis. Precisely, the approach's novel contribution consists of evaluating the frequency range of each of the muscle activations detected in time domain by the proposed algorithm. Acknowledged techniques, such as Fourier transform, are widely adopted to assess the power spectrum of the whole sEMG signal. However, this approach typically provides a single frequency content for all the possible different activations and the silent (only noise) portion of the signal [28,29]. The time-frequency approach of the proposed algorithm goes further, allowing to identify muscle activity intervals in time domain and concomitantly assessing the frequency range for each one of these activations, providing specific information in frequency domain about every single activation of the analyzed task. To our knowledge, the present analysis is the first to provide a reliable tool able to do this.
Moreover, it is acknowledged that detecting onset and offset events is hard to achieve in sEMG signals characterized by high background-noise levels, reflecting in low SNR values [19,41]. This is true for all the methods mentioned above. To address or at least improve this drawback, the Teager-Kaiser energy operator (TKEO) and its extended version (ETKEO) are introduced [42,43]. The aim is to filter the sEMG signal and, thus, increase SNR value to limit the errors in detecting activation events. Despite this, worsening detection performances could also be observed when TKEO is implemented to improve signal quality [41,44]. A further issue to consider is that SNR could assume different values in an other portion of the sEMG signal, especially in prolonged tasks like walking, running, cycling, and in specific conditions as muscle fatigue assessment. Variability of SNR within the same recording is likely due to the alteration of electrode-skin contact features, the change of noise power, and the variations in the ground reference level [24]. This can influence the assessment of onset-offset events just in the signal range where SNR degrades. For all these reasons, although many EMG-based methodologies are developed to identify muscle recruitment, a gold standard is not available yet.
Thus, the further goal of the present study is to attempt to limit the deterioration of algorithm performance associated with very noisy signals. To this aim, the CWT approach is chosen since it involves decomposing sEMG signal into many multiresolution components and applying a series of low-and high-pass filters, making this approach very suitable for handling signals affected by low SNR values [45]. Briefly, the algorithm consists of windowing the sEMG signal in short segments where SNR value is supposed to remain constant, filtering sEMG signal by CWT decomposition, estimating the local time-frequency energy density of the signal through the scalogram function [39], and then assessing onset and offset events as the begin and the end of the time interval where the scalogram is exceeding the 1% of the peak value of energy density. The windowing procedure is performed in order to handle the expected within-signal variability of SNR values.
Thus, in the present paper, we intend to first show that our algorithm is reliable in time domain since it can provide a quantitative time characterization of muscle activity (onset and offset events), at least comparable with the reference algorithm reported in the literature [23,39,46,47]. Then, we show that it can further provide new information, that is, the frequency content of every single activation detected in time domain.

A. ALGORITHM FOR ONSET AND OFFSET DETETCTION
Simulated and real sEMG signals are preliminary bandpass filtered (2nd order Butterworth filter, cut-off frequency 20-450 Hz). Later, all signals are processed to achieve the wavelet scalogram, according to the following procedure. Continuous Wavelet Transform (CWT) permits identifying time variations of non-stationary-signal frequency content, maintaining the resolution of signal processing in both frequency and time domain [48]. CWT of sEMG(t) is defined as the inner product shown in the successive Equation 1: where ψa,b(t) is the mother wavelet defined in Equation 2: where a is the scale parameter, b is the translation parameter (time shifting), and ψa,b(t) is achieved as the mother wavelet function ψ(t) computed at scale a in time b. ψa,b(t) must satisfy mathematical criteria like finite energy and no zero-frequency component to be admissible [49]. The scale is inversely proportional to spectral components. Low scales (i.e., high frequencies) provide more local information. High scales (i.e., low frequencies) supply more global information about the sEMG signal.
In the present study, wavelet transform is implemented by adopting Daubechis of order 4 with 8 levels of decomposition (db4) as the mother wavelet. This choice is made since Daubechies is acknowledged as a suitable mother wavelet for detecting signal changes and because its shape is similar to the shape of motor unit action potentials [50]. CWT decomposes a signal into several multiresolution components (coefficients) and performs a series of low-and high-pass filter operations followed by down-sampling. In this way, the sEMG signal is decomposed into its frequency content form and then reconstructed. CWT was used for removing noise from the sEMG signal and providing energy localization in time-frequency. CWT denoising is based on the decomposition of the signal, modification of detail coefficients concerning a defined soft threshold and then signal reconstruction concerning the new coefficients. This is achieved in the present study, applying the denoising algorithm included in Wavelet Matlab Toolbox. The soft threshold (Donoho threshold) is adopted for the denoising procedure. sEMG signal was reconstructed by revised CWT coefficients [51].
Energy localization in time-frequency of sEMG signal is identified using CWT scalogram function [45,51]. Scalogram function, PsEMG, is defined as the square of the absolute value of CWT coefficients, WsEMG, as reported in Equation 3: For each sEMG signal (simulated or real), the scalogram of denoised sEMG is computed in windows of a suitable portion of the sEMG signal (number of samples), depending on the signal characteristics (signal duration, muscle-activity duration). Details are reported in the following sections.
Timing of onset and offset events is assessed as the interval where the scalogram exceeds the 1% of the peak value of energy density in each specific portion of the sEMG signal. It is acknowledged that muscle recruitments lasting less than 30 ms have no effect in controlling joint motion [23]. Thus, activation intervals detected by the algorithm but lasting less than 60 samples are discarded. In the same way, activation intervals distant less than 60 samples are put together in a single activation interval. Once an activation interval is detected in time domain, the correspondent activation interval in frequency domain is computed as the frequency range associated with the scalogram function to that specific activation time interval. Maximum and minimum of the frequency content is, thus, quantified for each one of the activations identified in time domain.

B. ALGORITHM VALIDATION
Simulated and real data are considered to evaluate the appropriateness of the proposed approach to assess onsetoffset events of muscular recruitment in time-frequency domain. To this purpose, the following sub-section describes the generation of a test bench of simulated sEMG signals. The next sub-section introduces the experimental protocol adopted to design two different datasets of real sEMG signals. The current approach is validated by direct comparison with reference approaches on both real and simulated signals.

Simulated sEMG signals
Performances of the present approach in assessing onset and offset events are firstly evaluated on a test bench of simulated sEMG signals. Simulated sEMG signal is modelled as the overlapping of two stationary processes: the original signal generated by muscular activity and the background noise, fundamentally provoked by the crosstalk of neighbouring muscles and circuit noise [52]. Specifically, the background noise is simulated as a gaussian process with zero mean and variance 2 . The portion of the sEMG signal where the muscle is active is modelled as the background uncorrelated background noise plus the bandlimited stochastic process with zero-mean Gaussian distribution of amplitude and fixed power level, according to the acknowledged method introduced in [23]. This distribution is accomplished by band-pass filtering (cut-off frequencies = 80-120 Hz) a Gaussian series of uncorrelated samples [23]. This Gaussian distribution is truncated to simulate the sEMG activity due to muscle activation (Fig. 1).
Each simulated sEMG signal is generated with a sampling frequency fs = 2000 Hz, a time window of 1 s as suggested in [23], and a randomly variable value of the median of the gaussian µ, with 0 ≤ µ < 1. The physiological variability of sEMG signals is reproduced regulating the standard deviation, σ, and the time support, 2×α×σ, of the Gaussian distribution. Suitable values of σ are set in order to achieve the wanted value of signal-to-noise ratio (SNR), where: Randomly varying the value of parameter µ, 8 simulated sEMG signals are generated for all the different combinations of the values adopted for σ (50, 100, and 150 ), (1, 1.5, and 2), and SNR, (from 3 dB to 30 dB, with step = 3), considering the suggestions reported in [23]. A total of 720 simulated signals are obtained. CWT scalogram is computed separately for each one of those 2000-sample simulated signals.
The start and the end of truncated gaussian function used to model the simulated signal are adopted as the ground-truth events for the onset and offset, respectively.

Real sEMG signals
Two different datasets of real sEMG signals are used: the first one is adopted to validate the current approach; the second one is utilized to stress the potentiality of the present approach in time-frequency domain. The first one consists of a dataset of sEMG signals measured in 18 adult subjects carrying out knee extension and elbow flexion. It is made available by Tenan et al. at the following link (https://github.com/TenanATC/EMG) [18]. Correspondent ground truth is also provided. Two motor tasks are considered to evaluate the dependability of the approach under various conditions and in signals from different muscles. Participants carried out each exercise three times at a self-selected exercise pace. A mass (2.3 kg) is applied to subject's right ankle to improve the efficacy of the knee extension. To the same aim, elbow-flexion task is also carried out with the same mass applied to the subject's right wrist. During both tasks, the volunteer is seated in a stationary chair. sEMG probes are attached on biceps brachii (BB) for monitoring elbow flexion and vastus lateralis (VL) for knee extension. sEMG signals are analogue-to-digital converted at 2048 Hz. In total, 105 sEMG signals are acquired and then band-pass filtered with cut-off frequencies of 10-1000 Hz. sEMG signals are visually inspected by three experts to write down all the activation-onset in a randomized and double-blind fashion (i.e., specific task and study identifier are not known). Each expert checks every trial twice within a week. Average over the six onset values is adopted as ground truth for the experiments in [18], and is also used in the present study with the same aim. Additional details are reported [18]. In this dataset, for each signal CWT scalogram is computed in a window of 2000 samples around the ground-truth event.
Gait data characterize the second dataset of real data adopted in the present study. Specifically, sEMG and footfloor contact signals during 30-control-subject walking are retrospectively taken from the database assembled at Movement Analysis Lab, Università Politecnica delle Marche, Ancona, Italy and used for preceding studies of the present group of researchers [53,54]. The dataset is available by contacting the corresponding author of present and previous works. Obese and overweight volunteers (body mass index, BMI > 25) and subjects affected by any pathological condition, joint pain, or undergone orthopaedic surgery are excluded. Gait data are captured using multichannel recording system Step32 (Medical Technology, Italy; sampling rate: 2 kHz; resolution: 12 bit). sEMG signals are acquired in each leg by single differential probes placed over gastrocnemius lateralis (GL) and vastus lateralis (VL), respecting SENIAM guidelines for EMGsensor positioning [55]. Foot-floor-contact signals are measured by three footswitches placed under the heel, the first and the fifth metatarsal head of each foot and processed to segment the different gait cycles, according to [56]. Each subject walked barefoot at a self-selected pace for about 5 minutes, following an eight-shaped path, which involves natural deceleration, acceleration, and reversing. Additional details about data acquisition are available in [53]. The present research has been undertaken following the ethical principles of the Helsinki Declaration and was approved by the local ethical committee. In this dataset, CWT scalogram is computed separately in every single stride.

Algorithm performances
Performances of the proposed approach are evaluated in terms of precision, recall, F1-score, bias, and mean average error (MAE), assessed in all the true positives. Bias is assessed as the mean time distance between the predicted event and the correspondent ground-truth event. MAE is computed as the absolute value of the mean time distance between the predicted event and the correspondent groundtruth event. A predicted event occurring at time tp is acknowledged as true positive if a ground-truth event of the same kind occurs at time tg, such that | tg -tp | < T. T is a time tolerance set to 100 ms according to [40]. Differently, the predicted event is considered as false positive. After that, a further processing procedure is implemented to identify and then remove activation intervals that are too short to be physiologically plausible. According to [23], activation intervals ≤ 30 ms are discarded. For the same reason, consecutive activation intervals distanced less than 30 ms are united.

C. STATISTICS
Statistical difference of data distributions is evaluated. Firstly, Shapiro-Wilk test is adopted to appraise the normality of each data distribution. Next, a two-tailed, nonpaired Student's t-test is applied to verify the significance of

A. SIMULATED sEMG SIGNALS
Mean performances of onset and offset detection timing in simulated signals are shown in Table I; MAE and bias (on only  TPs), precision, recall, and F1-score are computed. Mean precision, recall, and F1-score are > 99% for the onset and > 98% for the offset (Table I). This means that the algorithm fails in identifying TPs in no more than 1% of the onset events. Nevertheless, MAE is computed in all the simulated signals to investigate if reporting MAE for TPs only could give false indication on comparing the two methods. Differences in MAE values are minimal (< 4%) compared with what was reported in Table I (TP only) and do not change the results of comparison between methods. Variability of MAE in function of α, σ, and SNR is quantified in Table II and III for onset and offset, respectively. A colour-level-coded portrayal is utilized to provide a useful visual picture of MAE variability. The worst values of MAE are located in the table region characterized by the highest value of α (2) and σ (100 and 150 ms) and the lowest value of SNR (3-12 dB), concomitantly. Nevertheless, MAE value never exceeds 50 ms. Results stratified for different SNR are highlighted in Fig. 2, in terms of mean MAE score (orange bars in panel A for onset and B for offset) and F1-score (orange bars in panel C for onset and D for offset). Specifically, mean F1-score never falls below 96% for onset and 94% for offset values, for the proposed approach. In the same way, mean MAE never exceeds 18 ms for onset and 22 ms for offset values. In the same panels, performances are directly compared with those obtained by DT algorithm in the same dataset. No significant differences are identified between F1-score values achieved by the two approaches in the whole SNR range (p > 0.05, Fig 2, panels C  and D). Similarly, no change is measured in MAE values (p > 0.05, Fig 2, panels A and B), except for SNR = 15 dB for the onset (p = 0.046) and for SNR = 12 dB for the offset (p = 0.031). Offset -MAE (ms) σ = 50 ms σ = 100 ms σ = 150 ms

B. REAL sEMG SIGNALS
The collection of sEMG signals assembled by Tenan et al. [18] is adopted to validate the current approach on real data. An example of raw and denoised sEMG signal in a representative subject is depicted in Fig 3. To test the effect of different SNR values, Tenan's dataset is firstly split, creating two subsets characterized by SNR ≤ 8 dB (first subset, 52 signals) and 8 < SNR ≤ 16 dB (second subset, 35 signals), respectively. Outcomes are presented in Table IV (first two rows) in terms of mean, standard deviation (SD), median, 25-percentile, and 75-percentile of the absolute error of onset prediction. Accuracy in onset prediction in both subsets is 100% (precision = recall = F1-score = 100%). No significant differences in the absolute error of onset prediction are detected between the two subsets (p > 0.05). Then, the first subset (SNR ≤ 8 dB) is further divided into four ranges of signals characterized by increasing SNR values with a step of 2 dB, as highlighted in the first column of Table IV. Outcomes are reported in the lower four rows of Table IV. Results in this four-range subsets are directly compared with results achieved by the following four different onset-detection algorithms in the same 52-signal dataset and reported in [44]: a method based on the doublethreshold statistical algorithm, DT [23]; the wavelet approach introduced in [39], WLT; the procedure founded on the CUSUM logic [46], CUSUM; and the recent approach grounded on the profile-likelihood maximization, employing  discrete Fibonacci search, PROLIFIC [47]. Two filtering procedures are also considered (TKEO and ETKEO). To facilitate the interpretation, results are shown all together in Table V. The present approach achieves the lowest absoluteerror values for all the metrics and SNR ranges (Table V), except for SD in signals with SNR ≤ 2.

C. TIME-FREQUENCY ANALYSIS
To demonstrate the capability of the proposed method to provide a time-frequency characterization of muscle activation, the present study adopts the dataset mentioned above composed of basographic and sEMG signals captured during walking of 30 control subjects. An example of CWT scalogram computed in a representative stride for GL and VL is reported in Fig. 4. In a total of 200 strides taken randomly from this dataset, three main activations for GL and four main activations for VL are identified in the gait cycle (GC). The mean value (± SD) over the 200-stride dataset of these activations is summarized in Table VI, for both GL and VL. A graphical representation is depicted in Fig. 5 to help visualize each main activation's position and duration within GC. For each activation identified in the time domain, the frequency range (Hz) is computed as the difference between the maximum and minimum frequency values and depicted in Fig. 6. In the upper panel of Fig. 6, frequency-range bars are split into three groups since three main activations are detected for GL. The first group (GL1) is characterized by the frequency-range bars of all the activations detected in early stance (mean value onset-offset = 4-14 %GC, as reported in Table VI). The second group (GL2) is characterized by the frequency-range bars of all the activations detected during push-off (mean value onset-offset = 21-45 %GC). The third group (GL3) is characterized by the frequency-range bars of all the activations detected during swing (mean value onsetoffset = 83-89 %GC). In the same way, in the lower panel, frequency-range bars are split into four groups since four main activations are detected for VL. The first group (VL1) is characterized by all the frequency-range bars of all the activations detected in early stance (mean value onset-offset = 0-15 %GC, as reported in Table VI). The second group (VL2) is characterized by the frequency-range bars of all the activations detected during push-off (mean value onset-offset = 33-42 %GC) The third (VL3) and the fourth (VL4) groups are characterized by the frequency-range bars of all the activations detected during early and late swing, respectively (mean value onset-offset = 72-78 %GC and 84-100 %GC). A large variability of the frequency range is detected for both muscles, ranging from 85 Hz to 489 Hz.

IV. DISCUSSION
The current study proposes a novel adaptative algorithm for assessing onset and offset timing of muscle activation in the time-frequency domain, based on CWT analysis. One of the added values of the current approach is to provide a concomitant assessment of the onset/offset timing and the frequency content of the muscle activation. Acknowledged techniques, such as Fourier transform, are widely adopted to assess power spectrum of the whole sEMG signal, including all the possible different activations and the silent (only noise) portion of the signal [28,29]. The time-frequency approach of the proposed algorithm allows identifying muscle activity intervals in time domain and concomitantly assessing the frequency range for each of these activations, providing specific information in frequency domain about every single activation of the analyzed task.
Three datasets are involved in testing the time-frequency  performances of the proposed algorithm: 1) a test bench of 720 simulated sEMG signals; 2) 105 real sEMG signals measured in vastus lateralis during knee extension and in biceps brachii during elbow flexion; and 3) real sEMG signals from gastrocnemius lateralis, and vastus lateralis collected in thirty healthy subjects during ground walking.

A. SIMULATED sEMG SIGNALS
The algorithm is tested in simulated data according to what was performed in previous studies, which have introduced novel algorithms for onset-offset detection [23,24,39,40]. The main reason is the need for a ground truth that provides an error-free assessment of onset-offset events for validating the algorithm. It is known that in real sEMG signals, the ground truth for onset and offset events is provided by trained human-movement specialists who inspect the signals and visually identify the onset-offset events. This procedure could be very susceptible to interexpert variability [20], specifically for the low signal-tonoise ratio value, as in the present study. Otherwise, the muscle activation is generated in the simulated sEMG signal truncating the Gaussian distribution ( Figure 1). Thus the onset and offset events are available with no error.
A total of 720 simulated sEMG signals are generated to quantify the ability of the algorithm to detect the onset and offset timing and assess its sensitivity to SNR variability. The proposed algorithm can achieve mean performance parameter (precision, recall, and F1-score) values very close to 100%, for both onset and offset detection, as shown in Table I. On average, correspondent values of bias and MAE are 6.0 ± 16.4 ms and 9.9 ± 14.4 ms for the onset and -1.1 ± 19.6 ms and 12.9 ± 14.8 ms for the offset, respectively. Vannozzi et al. reported higher mean bias (≈ 7 ms for both onset and offset) for their wavelet-based algorithm on a dataset of simulated sEMG signals with SNRs ranging from 8 to 20 dB [40]. Even higher bias values are achieved by Merlo et al. in applying a different wavelet-based method to simulated sEMG signals with 2 ≤ SNR ≤ 8 [39].
The analysis of MAE variability in function of SNR, α, and σ can be appreciated in Table II for onset and Table III for offset prediction, thanks to a colour-level-coded representation. As described in table legend, the darkest green zones highlight the lowest MAE values. Symmetrically, the darkest red zones show the highest MAE values. 10 ms is the threshold value adopted to discriminate red from green zones, since it is compliant with the mean MAE values displayed in Table I. Concomitance of highest values of α and σ seems to affect onset/offset detection more than SNR variability. This outcome suggests that detection capability gets worse with augmenting the activation timeduration, being the time support (i.e., the duration of a single activation) defined as 2×α×σ. Despite this, MAE never overcomes 50 ms, a value that seems to be compliant with the maximum acceptable response time of 300 ms, essential for natural control of an exoskeleton [57]. Furthermore, the mean worsening of onset/offset detection identified for the signals with the lowest SNR value (SNR = 3) compared to the signals with the highest value (SNR = 30) is very small, in terms of both F1-score (≈ 3% decrease for onset and ≈ 6% decrease for the offset) and MAE (10-ms increase for the onset and 11-ms increase for the offset). This suggests good robustness of the algorithm to SNR variability. A direct comparison with the results achieved by DT algorithm on the same simulated signals is adopted to evaluate the reliability of the proposed approach. The two methods show no significant difference (p > 0.05) between correspondent F1-score and MAE values in the whole SNR range (with a single exception, Fig. 2). Fig. 2 seems to indicate that DT reported a lower MAE value for onset detection for many SNR values. However, the statistical analysis shows that the possible differences between MAE values provided by DT and CWT are not significant (p > 0.05), except for one single SNR value. Otherwise, average bias values are lower for CWT prediction, even if not statistically (p > 0.05, Table I). For offset detection, the proposed CWT algorithm seems to work better than DT for the lowest SNR values, while the trend seems to invert for higher SNR values. Also, the possible differences are not significant (p > 0.05), except for one single SNR value. Thus, the concomitance of these results allows asserting that the capability of the two algorithms in detecting onset/offset events in simulated sEMG signals is comparable. Moreover, it is worth reminding that the analysis of the algorithm performances in time domain is achieved not to demonstrate that the proposed algorithm is able to outperform the other approaches in assessing onset and offset events in time domain, but just to show that the CWT algorithm is reliable in time domain since it can provide an assessment of onset and offset events at least comparable with reference algorithm reported in the literature [23,27,39,46,47]. The actual aim of the present study is, indeed, to propose an algorithm able to provide further and new information, that is, the frequency content of every single activation detected in time domain, as highlighted in the following section C.
All these findings contribute to conclude that proposed approach is able to provide a reliable prediction of muscle onset and offset events in a large dataset of simulated sEMG signals, being minimally affected by large SNR variability (3-30 dB).

B. REAL sEMG SIGNALS
The reliability of onset detection in real sEMG signals is tested on Tenan's database [18]. Detection of offset timing is not computed, given that the ground-truth offset is not available. The choice of this dataset has been driven by the following considerations: it is a large and complete sEMG dataset available for free; it includes a reliable ground truth for onset detection, which is already used and tested in previous studies [19,44,47]; the peculiarities of knee extension and elbow flexion allow a simple and accurate visual detection of the muscle onset and, consequently, reliable identification of the ground truth; results of the application to this dataset of four onset-detection techniques and two filters are available for a direct comparison [44]. A recent and exciting study proposes an onset and offset detection algorithm for muscle activation validated during hand-close movement of 10 subjects [27]. It is worth noticing that the outcomes of the present study are at least comparable, in terms of MAE ± SD, with those reported in [27], although the algorithm is tested on a different dataset and in a different experimental condition. These results further support the reliability of the present algorithm in detecting muscle onset in time domain.
On the selected dataset composed by the 52 signals with SNR ≤ 8 dB, the proposed algorithm globally outperforms each one of the approaches employed in [44], achieving the lowest onset-detection absolute error, in terms of all error parameters (mean, SD, median, 25th percentile, and 75th percentile; bold red characters in Table V). This encouraging outcome is associated to a perfect detection accuracy of the proposed algorithm (precision = recall = F1-score = 100%). With the aim of a more detailed comparison, the SNR range is split in further four subsets. The proposed approach provides the lowest absolute error values also in these SNRdependent subsets, as highlighted with bold black characters in Table V (with a single exception for SD when SNR ≤ 2 dB). These results are also confirmed when the four onsetdetection reference algorithms are combined with the advanced filtering technique ETKEO (Table V).
The effect of different SNR values is evaluated by comparing the detection performances obtained by the proposed algorithm in signals characterized by low SNR values (SNR ≤ 8 dB) vs. signals with high SNR value (8 dB < SNR ≤ 16 dB). The absence of significant differences (p > 0.05) between the two subsets strongly supports the algorithm's robustness to SNR between-signal variability suggested by the analysis of simulated signals. Moreover, to handle the expected within-signal variability of SNR values, the proposed approach adopts a windowing procedure, splitting the sEMG signal in 2000-sample segments where SNR value is supposed to remain constant. Thus, the main procedure of the algorithm (denoising, estimation of scalogram function, removal of the portion of the signal < 1% of the peak value of energy density) is applied in an adaptative way to each one of the 2000-sample windows, for both real and simulated signals. Since each signal is split into 2000-sample segments, the denoising, the scalogram function, and the peak value of energy density are computed for each segment (i.e., every 2000 samples). This means that the threshold adopted for computing onset and offset timing (1% of the peak value of energy density) is not a fixed value, but it changes from signal to signal and from segment to segment within the same signal, adapting to the peculiar characteristics of the sEMG signal in that specific segment of the signal. sEMG signals are split into 2000-sample segments for three main reasons: 1) for considering sEMG-signal segments of the same (comparable) duration in the three datasets adopted in the present study (2000 samples for simulated signal, 2000-2100 samples for a signal during walking, and thus 2000 samples also for signals coming from Tenan's dataset); 2) for implementing an adaptative approach, as stressed above; 3) because SNR value is supposed to remain constant in such a small segment.

C. FREQUENCY ANALYSIS
Walking is one of the most suitable motor tasks to get insights into human motion. Thus, the capability of the proposed algorithm to characterize muscular activity in timefrequency domain is assessed in a dataset of sEMG data collected from GL and VL during 30-healthy-adult walking.
As reported in Fig. 5 (upper panel), three main activations are identified for GL. Mean activation intervals depicted in Fig. 5 match what reported in reference [58] and more recent [59] literature. In particular, the activation starting after the end of the heel rocker and finishing at the heel off is the main GL activity (GL2), interpreted as the active participation of GL in restraining the forward progression of the tibia over the talus during the second rocker, hence controlling dorsiflexion [58]. The other two activations (GL1 and GL3) occur less frequently for control purposes, as discussed in [59]. Four main activations are detected for VL during GC (Fig. 5, lower panel). The two activations at the beginning (VL1) and the end (VL4) of GC are acknowledged as the typical vastii recruitment to generate tension in the terminal swing in preparation for weight-bearing at initial contact and to control knee flexion during weight acceptance [58]. The other two activations highlighted in Fig. 5 are consistent with what reported in [60]: the first one (VL2) is likely adopted to stabilize the patella before the pre-swing phase, and the second one (very sporadic, VL3) is connected to the activity of vastii at the end of GC.
As shown in previous chapters, the performance of the proposed CWT method in time domain is tested against the double threshold (DT) algorithm [23] in the same dataset of simulated sEMG signals. Then, CWT method is also tested against four reference algorithms, the same DT used for simulated data, WLT, CUSUM, and PROLIFIC [23,39,46,47] in the same dataset of real sEMG signals available in the literature [18]. Besides the promising outcomes achieved in time domain, the present algorithm's actual novel contribution consists of providing a characterization of the frequency content of each single muscle activation detected in time domain. To our knowledge, this is achieved here for the first time (Fig. 6). Thus, a direct comparison with other approaches is not possible. However, some exciting considerations could be made. As expected, every activation presents a frequency content consistent with the typical frequency spectrum of the sEMG signal (10-500 Hz) [61], supporting the robustness of the proposed approach.
Furthermore, the large variability of the frequency range is detected within the muscle and between muscles. In particular, the frequency range associated with the GL3 group seems to be reduced compared to the other activation groups (Fig. 6). In the same way, the frequency range associated with the VL2 group decreases compared to VL1 and VL4. A further difference between muscles is evident. The reduction of the frequency range is detected mainly in groups of activations located in the late swing for GL (GL3); differently, the decrease of VL frequency range is in midstance (VL2). A possible explanation of these phenomena could lie in the different functional tasks of each activation related to muscle recruitment during GC. Although quantitative analysis of the frequency content associated with different activations is beyond the goal of the present study, these findings raise a new issue that deserves to be investigated. Novel and specific tools, like the one proposed here, could be beneficial to this aim.

V. CONCLUSION
The high accuracy of outcomes and the successful validation performed with a direct comparison against acknowledged detection algorithms endorse the reliability of continuous wavelet transform and, expressly, of the proposed algorithm in detecting muscular recruitment in the time-frequency domain. The involvement of both simulated data and real sEMG signals (from two different motor tasks: elbow flexion and knee extension) supports the general validity of the present findings.
The novel contribution of the current approach is the quantification not only of the onset-offset timing but mainly the frequency content of each one of the activations detected in the time domain. Moreover, the low sensitivity to between-subject and within-signal SNR variability makes this approach suitable for all experimental conditions. Limited values of detection error make the current approach compatible with the timing of sEMG-driven assistive devices. However, further developments should focus on evaluating the suitability of the present approach for on-field applications.