Muscle fatigue threshold detection method based on surface electromyogram signals

Abstract

The invention relates to the technical field of biomedical signal processing, and discloses a muscle fatigue threshold detection method based on surface electromyogram signals, which comprises the following steps: acquiring synchronous surface electromyogram and angular velocity signals, and segmenting to generate a motion period sequence; establishing a baseline kinematics envelope template based on the initial angular velocity data; performing morphological verification on the current angular velocity segment, and calculating actual electromechanical delay after verification is passed; comparing the reference delay to obtain a time drift distance, and defining an electromyographic signal interception time boundary in combination with an initial matching window; intercepting a signal segment in the calibration boundary, calculating a frequency band energy distribution index and generating a time sequence; the sequence is reconstructed into a two-dimensional phase space state vector, the maximum Lyapunov index is calculated, when a continuous counting mechanism is met, it is judged that dynamic topology bifurcation occurs, and a control intervention instruction is output. According to the invention, electromechanical phase dislocation caused by fatigue can be dynamically compensated, and accurate detection and safe intervention of the muscle fatigue state are realized.

Description

Technical Field

[0001] This invention relates to the field of biomedical signal processing technology, specifically to a method for detecting muscle fatigue threshold based on surface electromyography signals. Background Technology

[0002] Surface electromyography (EMG) signals, as an objective indicator that can non-invasively reflect the electrophysiological activity of the neuromuscular system, are widely used in rehabilitation medicine, sports biomechanical analysis, and exoskeleton robot control. During periodic continuous movement in the human body, it is usually necessary to extract EMG signal segments from each movement cycle and assess the fatigue state of the target muscle group by extracting time-frequency domain features or nonlinear dynamic indices.

[0003] Existing methods for detecting muscle fatigue typically rely on sampling signals within a fixed time window set at the initial, fatigue-free state. However, during continuous dynamic work, the physiological functions and kinematic characteristics of the human body undergo dynamic changes as muscle fatigue accumulates. On one hand, fatigue leads to a decrease in muscle contraction rate, and the intervention of neural compensatory mechanisms deforms the movement trajectory, resulting in changes in the actual duration and rhythm of a single movement. On the other hand, the accumulation of metabolic products (such as lactic acid) within muscles reduces the efficiency of calcium ion channel release and reabsorption within muscle fibers, leading to a prolonged lag time between the transmission of neural excitation commands and the generation of actual mechanical force, i.e., increased electromechanical delay.

[0004] Because existing methods lack dynamic compensation mechanisms for the evolution of the aforementioned time-domain parameters, continuing to use fixed time boundaries to extract electromyographic signals will cause the sampling interval to deviate from the actual electrophysiological force exertion phase. This phase misalignment means that the extracted signal segments cannot correspond to the actual muscle contraction phase, thus introducing a large amount of irrelevant data from non-target phases. If nonlinear phase space reconstruction and Lyapunov exponent calculation are directly performed on signal sequences with phase misalignment, it will cause discontinuity in the phase space trajectory, making the final output dynamic characteristics unable to reflect the physiological decline pattern of muscles, leading to misjudgment or underjudgment of the critical threshold of muscle fatigue, and making it difficult to provide timely and safe control intervention. Summary of the Invention

[0005] To address the shortcomings of existing technologies, this invention provides a method for detecting muscle fatigue threshold based on surface electromyography (EMG) signals. This method solves the problem that in continuous dynamic movement, the increased muscle fatigue leads to changes in movement rhythm and increased electromechanical delay, resulting in force phase misalignment when using traditional fixed time windows to capture EMG signals, thus causing inaccurate assessment of muscle fatigue status.

[0006] To achieve the above objectives, the present invention provides the following technical solution:

[0007] This invention provides a method for detecting muscle fatigue threshold based on surface electromyography signals, comprising the following steps:

[0008] Acquire surface electromyography (EMG) signal sequences and angular velocity signal sequences, and segment them to generate motion cycle sequences;

[0009] A baseline kinematic envelope template is established based on the angular velocity data of the initial number of motion cycles in the motion cycle sequence, defining the initial kinematic phase matching window and the reference delay time;

[0010] Extract the surface electromyography (EMG) signal fragments and angular velocity signal fragments of the current movement cycle, perform morphological calibration of the angular velocity signal fragments with the baseline kinematic envelope template, and calculate the actual electromechanical delay after the morphological calibration is passed.

[0011] The electromechanical delay time drift is obtained by subtracting the actual electromechanical delay from the reference delay time, and the electromyographic signal interception time boundary is defined by combining the initial kinematic phase matching window.

[0012] Within the time boundary of electromyography signal interception, surface electromyography signal segments are intercepted to calculate the frequency band energy distribution index and generate a time series of the frequency band energy distribution index.

[0013] The frequency band energy distribution index time series is reconstructed into a two-dimensional phase space state vector, the maximum Lyapunov exponent is calculated, and a control intervention command is output when the maximum Lyapunov exponent satisfies the continuous counting mechanism.

[0014] Specifically, after filtering and preprocessing the synchronously acquired surface electromyography (SEM) signal sequence and angular velocity (Angular velocity) signal sequence, the data points of the filtered Angular velocity signal sequence are traversed according to the time sequence, and the moments when the signal values ​​cross zero level are marked as zero-crossing points. Two adjacent in-phase zero-crossing points on the time axis are used as the start and end boundaries of the motion interval. The timestamps of the start and end zero-crossing points of the motion interval are read, and the filtered SEM and Angular velocity signal sequences are extracted within the start and end boundaries to obtain the corresponding SEM and Angular velocity signal segments. The difference between the end and start zero-crossing point timestamps is calculated to obtain the corresponding actual duration. The SEM and Angular velocity signal segments containing timestamp information are bound and stored with their corresponding period indices and actual durations, thus segmenting the continuous signal into discrete motion period sequences.

[0015] Based on this, angular velocity signal segments and surface electromyography (EMG) signal segments within a continuous initial number of motion cycles are extracted as sample data. The angular velocity signal sequence is time-axis normalized, and the average amplitude of the angular velocity waveform is calculated at each relative time point to generate a baseline kinematic envelope template. The peak angular velocity of the baseline kinematic envelope template is extracted, and continuous data segments with amplitudes between the peak angular velocity multiplied by a set lower amplitude limit coefficient and the peak angular velocity multiplied by an upper amplitude limit coefficient are identified. The interval covered by these continuous data segments on the relative cycle progress axis is established as the initial kinematic phase matching window, and its relative start time percentage and relative end time percentage are recorded. The starting point of EMG activity and the extreme points of angular velocity signals in the sample data are extracted respectively, their absolute time differences are calculated, and the arithmetic mean is taken as the reference delay time.

[0016] To ensure the validity of the analyzed data, the angular velocity waveform of the current cycle is resampled using an interpolation algorithm to make its data length the same as the baseline kinematic envelope template, and the Pearson correlation coefficient between the two is calculated. If the correlation coefficient is less than the isomorphism threshold, the frequency band energy distribution index of the previous motion cycle is directly output; if the correlation coefficient is greater than or equal to the isomorphism threshold, the action is deemed valid. After verification, the Teager-Kaiser energy operator is used to perform a nonlinear transformation on the surface electromyography signal to obtain the instantaneous energy sequence. The mean and standard deviation of the background noise within the resting data window are calculated to establish an adaptive trigger threshold. The timestamps corresponding to the data points in the instantaneous energy sequence that first continuously exceed the threshold for a set duration are used as the neural excitation pacing points. Simultaneously, the derivative of the original angular velocity signal is used to obtain the timestamps corresponding to the extreme points of angular acceleration as the extreme points of mechanical force. The absolute time difference between the two is calculated to obtain the actual electromechanical delay of the current motion cycle.

[0017] Subsequently, the difference between the current actual electromyographic delay and the baseline delay time is calculated to obtain the electromyographic delay time drift. Using the start zero-crossing timestamp and the actual duration combined with the percentage of the relative start and end times, the start and end baseline time points are calculated respectively. The electromyographic delay time drift is subtracted from each of the two baseline time points to define the time boundary for electromyographic signal extraction. The essence of this process is to reverse-shift the time extraction window of the electromyographic signal by measuring the actual neural-mechanical force exertion delay change in a single movement, establishing a dynamic time compensation mechanism to ensure that the extracted frequency band features are strictly aligned with the actual force exertion phase of the human body.

[0018] After truncating the aligned signal segment, zeros are padded and a Discrete Fourier Transform (DFT) is performed to obtain the power spectral density sequence. Using the initial power spectral median frequency as the critical segmentation frequency, the frequency interval is divided into low-frequency and high-frequency bands, and discrete summation is performed on each. The ratio of the energy integral values ​​of the low-frequency and high-frequency bands is calculated to obtain the frequency band energy distribution index sequence. The optimal delay time is calculated using the mutual information method, and the current data point and the lagging data point are extracted from the sequence to form a two-dimensional phase space state vector.

[0019] Finally, the state vectors in the two-dimensional phase space are traversed to find the reference state vector and its nearest neighbor state vector. Both vectors are simultaneously moved backward, and the Euclidean distance between the state vectors at different evolution step sizes is calculated. The average of the logarithmic distances is extracted to generate a local divergence evolution curve. A least-squares straight line is fitted to the linearly rising region of the initial stage of the curve, and the slope of the fitted line is used as the maximum Lyapunov exponent. The average exponent of the initial period is multiplied by a proportionality margin to generate a bifurcation threshold. An anomaly counter is used to continuously check for exceeding the current maximum Lyapunov exponent. When the anomaly count reaches the number of tolerance cycles, a dynamic topological bifurcation is determined in the target muscle tissue, and a control intervention command is generated. The command is encapsulated into a data frame and sent to the rehabilitation exoskeleton device for resistance adjustment, to the front-end interface for state updates, or to the electrical stimulation device for shutdown protection.

[0020] This invention provides a method for detecting muscle fatigue threshold based on surface electromyography signals. It has the following beneficial effects:

[0021] 1. This invention obtains the electromyographic delay time drift by calculating the difference between the actual electromyographic delay and the reference delay time in the current movement cycle, and uses this to dynamically compensate for the interception time boundary of the electromyographic signal. This method solves the problem that the traditional fixed time window cannot adapt to the decrease in neuromuscular conduction rate during fatigue, enabling the intercepted electromyographic signal segments to correspond to the actual contraction stage of the muscle, thus improving the reliability of effective data extraction.

[0022] 2. This invention performs morphological calibration of the angular velocity signal segment of the current motion cycle with a baseline kinematic envelope template and calculates the Pearson correlation coefficient. When morphological distortion is detected, a state-preserving strategy is implemented. This mechanism can identify abnormal motion cycles caused by compensatory force exertion or movement deformation during muscle fatigue detection and uses the frequency band energy distribution index of the previous cycle as a replacement output, avoiding interference from invalid data in the dynamic state analysis and enhancing the method's anti-interference capability in real-world motion scenarios.

[0023] 3. This invention reconstructs the time series of the frequency band energy distribution index into a two-dimensional phase space state vector and calculates the maximum Lyapunov exponent. Simultaneously, a continuous counting mechanism is introduced for fatigue threshold determination. This method utilizes nonlinear dynamic indicators to assess the state of the muscle system maintaining stable contraction and eliminates occasional physiological noise fluctuations through continuous counting conditions, reducing the system's false judgment rate and improving the accuracy of fatigue detection and control intervention command output. Attached Figure Description

[0024] Figure 1 This is a schematic diagram of the method flow of the present invention;

[0025] Figure 2 This is a schematic diagram of the system architecture of the present invention;

[0026] Figure 3 The diagram shows a comparison and verification of the evolution trends of the traditional median frequency index and the maximum Lyapunov exponent of the present invention. (a) is the evolution trend of the median frequency of the traditional power spectrum; (b) is the evolution trend of the maximum Lyapunov exponent. Detailed Implementation

[0027] The technical solutions in the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0028] Please see the appendix Figure 1 This invention provides a method for detecting muscle fatigue threshold based on surface electromyography signals, comprising the following steps:

[0029] S1, acquire surface electromyography signal sequence and angular velocity signal sequence, segment to generate motion cycle sequence, and extract the current motion cycle and corresponding angular velocity signal segment;

[0030] S2, Based on the angular velocity data of the previous few motion cycles, establish a baseline kinematic envelope template and define the initial kinematic phase matching window and reference delay time;

[0031] S3, perform morphological calibration on the angular velocity signal segment and the baseline kinematic envelope template, and calculate the actual electromechanical delay of the current motion cycle after the calibration is passed;

[0032] S4, the electromechanical delay time drift is obtained by subtracting the actual electromechanical delay from the reference delay time, and the electromyographic signal interception time boundary is defined by combining the initial kinematic phase matching window;

[0033] S5, extract surface electromyography signal segments within the time boundary of electromyography signal extraction, calculate the frequency band energy distribution index of the current movement cycle and generate a frequency band energy distribution index time series;

[0034] S6 reconstructs the frequency band energy distribution index time series into a two-dimensional phase space state vector, calculates the maximum Lyapunov exponent of the current motion cycle, and determines the occurrence of dynamic topological bifurcation when the continuous counting mechanism is satisfied, and outputs control intervention commands.

[0035] Please see the appendix Figure 2 This invention provides a muscle fatigue threshold detection system based on surface electromyography signals, comprising: a data acquisition hardware environment and a core data processing logic module.

[0036] The data acquisition hardware environment includes a surface electromyography (SEMG) sensor, an inertial measurement unit (IMU), and a microprocessor. The SEMG sensor and IMU are worn on the surface of the target muscle and the corresponding movable joint, respectively. Both the SEMG sensor and IMU are connected to the microprocessor via a communication bus. The microprocessor has a unified hardware clock source internally. Using this hardware clock source, the microprocessor assigns timestamps to the received SEMG and angular velocity signals, achieving synchronous alignment of the acquired data on the time axis.

[0037] The core data processing logic module is deployed inside the microprocessor. This core data processing logic module includes a data preprocessing module, a morphological verification module, a time axis compensation module, a feature extraction module, and a dynamics determination module.

[0038] The data preprocessing module receives synchronization signals acquired by the surface electromyography (EMG) sensor and the inertial measurement unit (IMU). Based on the envelope characteristics of the angular velocity signals, the data preprocessing module segments the continuous data into a sequence of motion cycles, using the zero-crossing point as the boundary. In the initial stage of system operation, the data preprocessing module establishes a baseline kinematic envelope template based on the angular velocity data of the previous few motion cycles and defines the initial kinematic phase matching window.

[0039] The morphological calibration module is connected to the data preprocessing module. Starting from the period after the initial stage, the morphological calibration module calculates the rhythm scaling factor between the actual duration and the baseline duration. The morphological calibration module resamples and aligns the angular velocity waveform of the current period with the baseline kinematic envelope template and calculates the Pearson correlation coefficient. When the correlation coefficient is below a preset threshold, the morphological calibration module executes a state preservation strategy, stopping data updates for the current period. When the correlation coefficient reaches or exceeds the preset threshold, the morphological calibration module transmits the current period data to the time axis compensation module.

[0040] The time axis compensation module receives the validated motion cycle data. It uses the Teager-Kaiser energy operator to extract the neural excitation pacemaker of the surface electromyography (EMG) signal and differentiates the angular velocity signal to extract the extreme points of angular acceleration as the extreme points of mechanical force. The time axis compensation module calculates the time difference between these two points as the actual electromechanical delay. Based on the difference between this actual electromechanical delay and the reference delay, the time axis compensation module calculates the time drift compensation amount and, in conjunction with the rhythm scaling factor, calculates the sampling window time boundary for the next motion cycle.

[0041] The feature extraction module is connected to the time axis compensation module. The feature extraction module extracts a segment of surface electromyography (EMG) signal within the calibrated time boundary and performs a Fast Fourier Transform (FFT) on this signal segment. The feature extraction module calculates the ratio of the energy integral of the preset high-frequency sensitive band to the energy integral of the low-frequency sensitive band, generating a frequency band energy distribution index. The feature extraction module outputs a feature time series along with the motion cycle.

[0042] The kinetics determination module receives the feature time series output by the feature extraction module. It maps the frequency band energy distribution index and its first-order difference to a two-dimensional phase space state vector. The kinetics determination module calculates the maximum Lyapunov exponent of the trajectory evolution of adjacent state phase points in the two-dimensional phase space. When the kinetics determination module detects that the value of the maximum Lyapunov exponent is equal to or greater than zero, and this state is maintained for a preset number of consecutive cycles, the kinetics determination module outputs a trigger signal representing the arrival of the muscle fatigue threshold.

[0043] The data preprocessing module receives the surface electromyography (EMG) signal sequence with timestamps assigned by the microprocessor. and angular velocity signal sequence The continuous data stream is then denoised and segmented.

[0044] The data preprocessing module processes the synchronously acquired surface electromyography signal sequences. and angular velocity signal sequence Filtering preprocessing is performed on the surface electromyography signal sequence. The data preprocessing module uses bandpass filters and notch filters for frequency domain truncation. The low-frequency cutoff and high-frequency cutoff frequencies of the bandpass filters are configured to 20Hz and 500Hz, respectively, to extract the effective frequency band of muscle electrophysiological activity and filter out low-frequency motion artifacts and high-frequency thermal noise. The center frequency of the notch filter is configured to 50Hz or 60Hz to eliminate power grid interference in the environment. This is applied to angular velocity signal sequences. The data preprocessing module uses a low-pass filter for smoothing. The cutoff frequency of this low-pass filter is configured with a preset value between 5Hz and 10Hz to attenuate high-frequency mechanical vibration noise during limb movement, while retaining the low-frequency fundamental component characterizing the macroscopic movement trajectory of the limb. Generally speaking, the effective energy of muscle electrophysiological signals is mainly concentrated in the mid-to-high frequency range, while the frequency of macroscopic limb mechanical movement in the human body is relatively low. Through the above-mentioned frequency range division based on the characteristics of different physical quantities, the data preprocessing module can suppress interference from the external environment and the device itself while preserving the true physiological and motor characteristics. Regarding the specific digital implementation methods of the bandpass filter, notch filter, and low-pass filter, such as using a Butterworth filter or a finite-length unit impulse response filter, those skilled in the art can perform conventional configurations based on the computing power overhead of the controller. The specific discretization difference equation solution is a well-known technique in the field and will not be elaborated upon here.

[0045] The data preprocessing module is based on the filtered angular velocity signal sequence. Zero-crossing point detection and motion cycle segmentation are performed. In periodic variable load or variable posture motion tasks, the reciprocating motion of the target joint will result in a sequence of angular velocity signals. The values ​​of angular velocity alternate between positive and negative. When the human body performs periodic movements (such as gait or continuous squats), the joints move in two opposite directions. The moment when the angular velocity is zero usually corresponds to the switching point of the joint's movement direction or the extreme position of the movement. The data preprocessing module iterates through the filtered angular velocity data points according to the time series, extracting the moments when the signal value crosses the zero level and marking them as zero-crossing points. The data preprocessing module uses two adjacent in-phase zero-crossing points on the time axis—that is, two consecutive zero-crossing points that change from negative to positive or from positive to negative—as the start and end boundaries of the movement interval. The data preprocessing module defines all data segments within these start and end boundaries as a complete movement cycle. Using these physical extreme positions as dividing boundaries allows mapping the actual mechanical movement cycle of the human body.

[0046] The data preprocessing module assigns a cycle index to each segmented motion cycle. ,in The data preprocessing module reads the first... The starting zero-crossing timestamp of each cycle Intersection timestamp with end zero And calculate the number according to the following formula. The actual duration of each cycle :

[0047] ;

[0048] in, This is the actual duration. To end the zero-crossing timestamp, The starting zero-crossing timestamp is used. The data preprocessing module combines the surface electromyography (EMG) signal fragments and angular velocity signal fragments containing the timestamp information with their corresponding period indices. and actual duration The data is then bound and stored. Continuously acquired signals are segmented into discrete motion cycle sequences, which form the corresponding discrete data segments for subsequent calculations of period-based feature variables.

[0049] After segmenting the motion cycle sequence, the data preprocessing module extracts baseline features for the initial stage of system operation as reference data in a non-fatigue state.

[0050] The data preprocessing module calculates the average baseline cycle time during the initial phase of system operation. In the initial stages of a cyclical movement, the muscle exertion patterns and duration exhibit high repetitiveness. The data preprocessing module extracts continuous data from the preceding movements. Using the motion cycles as sample data for establishing the baseline model, the above calculations were performed. Actual duration of each cycle The arithmetic mean of the initial number of motion cycles. The value range is typically configured to be 3 to 5 consecutive periods, which can represent the motion characteristics of the system in its initial state. The calculation process uses the following formula:

[0051] ;

[0052] in, This represents the average baseline period duration. This represents the initial number of motion cycles; The summation symbol; Index for motion cycles; The starting value of the motion cycle index is 1; Indicates the index of the motion cycle From 1 to The data are summed up. This refers to the actual duration.

[0053] The data preprocessing module generates a baseline kinematic envelope template. Due to the actual duration of continuous human movement Fluctuations will occur between each cycle, directly affecting the angular velocity signal sequence. Overlaying data can cause misalignment on the timeline. The data preprocessing module... angular velocity signal sequence within one period Time axis normalization is performed. The data preprocessing module maps the absolute time intervals of each cycle to a relative cycle progress axis of 0 to 100%, and performs time axis normalization at each relative time node. The amplitudes of each waveform are averaged. After superposition and averaging, a baseline kinematic envelope template representing the force exertion law of the target action is generated. For data resampling during waveform time axis normalization, those skilled in the art can use linear interpolation or cubic spline interpolation. The specific interpolation mapping calculations are well-known techniques in the field and will not be elaborated here.

[0054] The data preprocessing module is based on the baseline kinematic envelope template. Define the initial kinematic phase-matching window. From the perspective of human muscle dynamics, the isotonic contraction phase of muscles during dynamic continuous movement is the observation interval for extracting electrophysiological features for fatigue assessment. This contraction phase corresponds, macroscopically, to a specific waveform interval before and after the limb's velocity accelerates to its peak. The data preprocessing module traverses the baseline kinematic envelope template. Data points, extract peak angular velocity from waveform The data preprocessing module is internally configured with an amplitude lower limit proportional coefficient. and amplitude upper limit ratio coefficient The range of values ​​for these two coefficients satisfies In practical applications, the amplitude lower limit proportionality coefficient The upper limit of amplitude can be configured to be 0.3 to 0.5. Configurable to 0.8 to 0.9, used to extract the main force waveband of the motion. The data preprocessing module uses the baseline kinematic envelope template. Within the peak range, find the amplitude located at and The data preprocessing module defines the interval covered by this data segment on the relative periodic progress axis as the initial kinematic phase matching window. The data preprocessing module records the percentage of the window's relative start time on the normalized time axis. and relative end time percentage Using a relative time percentage as a boundary variable allows for synchronous scaling of the intercepted interval when the overall motion speeds up or slows down, maintaining the consistency of the sampled physical phase.

[0055] The data preprocessing module calculates the system's baseline delay time. Establishing the baseline kinematic envelope template During the sample period, the data preprocessing module synchronously reads the previous data. Surface electromyography signal sequence within one cycle The data preprocessing module extracts these respectively. The starting points of surface electromyography (EMG) activity within each cycle and the corresponding extreme points of mechanical force are used to calculate the absolute time difference between the two points. The method for extracting and calculating this absolute time difference is the same as the calculation method for a single cycle in the time axis compensation module described later. The data preprocessing module will measure the... The arithmetic mean of the absolute time differences is calculated, and this arithmetic mean is set as the base delay time. Reference delay time The baseline electromechanical coupling time, during which electrophysiological commands are transmitted to muscles and trigger limb mechanical movements, is recorded when the system is in a non-fatigue state. This parameter is cached within the microprocessor as a benchmark for subsequent detection of the increase in electromechanical delay and calculation of drift compensation.

[0056] After baseline feature extraction is completed, the system enters the continuous monitoring phase. This begins from the initial phase, i.e., the [number]th cycle... At the start of each motion cycle, the morphological verification module performs verification operations on each subsequent received motion cycle to prevent non-standard movements from introducing abnormal data. Among these... This represents the initial number of motion cycles.

[0057] During continuous work or rehabilitation training, the duration of a single movement dynamically changes as muscle fatigue deepens or external loads alter. From the perspective of muscle physiology and kinesiology, as motor units fatigue, the rate of muscle contraction decreases, leading to a prolonged absolute time for the limb to complete the prescribed movement; or, due to compensatory force mechanisms, the movement trajectory deforms, causing the movement to end prematurely. If the system continues to use a fixed absolute time length from the baseline state to extract surface electromyographic signals, the actual extracted electrophysiological signal segments will be misaligned with the phase of the muscle's actual physical force exertion. To address this issue, the morphological verification module extracts a rhythm scaling factor to quantify the degree of stretching or compression on the time axis.

[0058] The morphological school verification module obtains the actual duration. The current motion cycle index is... And satisfy The actual duration The acquisition method is the same as in the initial stage, that is, the data preprocessing module calculates the difference between the end zero crossover time stamp and the start zero crossover time stamp of the current cycle, and transmits the difference data to the morphological verification module.

[0059] The morphological school module calculates the rhythm scaling factor of the current motion cycle. The morphological school verification module will determine the actual duration. Divide by the average baseline cycle time of the system's internal cache The calculation process uses the following formula:

[0060] ;

[0061] in, This is the rhythm scaling factor for the current motion cycle; This refers to the actual duration. This represents the average baseline period duration.

[0062] Rhythm scaling factor As a dimensionless parameter, it characterizes the degree of variation of the current limb movement relative to a baseline movement over time. When calculated... When this occurs, it indicates that the current movement cycle is longer than the baseline state, and the body's current movement rhythm is slower; when the calculation shows... When the time is short, it indicates that the duration of the current movement cycle is shorter than the baseline state, and the current movement rhythm of the human body is faster. In practical applications, normal rhythm changes in the human body are usually within a certain physical range, and the rhythm scaling factor... The effective value range is generally configured between 0.5 and 2.0. If the value exceeds this range, it usually corresponds to extreme situations such as sensor detachment or a fall. This parameter enables the system to establish a time mapping relationship between the current action and the reference action, providing a data basis for the subsequent time axis compensation module to adjust the absolute time boundary of the sampling matching window on a relative time scale. Through the above extraction and calculation process, the specific lower-level operation of dynamic tracking of human movement rhythm is realized.

[0063] After calculating the rhythm scaling factor for the current motion cycle, the morphological verification module performs a morphological similarity assessment on the angular velocity data of the current motion cycle to filter out distorted data introduced by non-standard movements.

[0064] Because the actual duration of each motion varies, the number of data points on the time axis of the current cycle's angular velocity waveform differs from that of the baseline kinematic envelope template, making direct point-by-point comparison impossible. The morphological verification module resamples the current cycle's angular velocity waveform using an interpolation algorithm to ensure its data length matches that of the baseline kinematic envelope template. Since the baseline kinematic envelope template maps the absolute time interval to a relative cycle progress axis from 0% to 100% during construction, the total number of data points after resampling... The number of data points is typically fixed at 100 or 1000 equal division points to ensure that each data point corresponds to the same action phase. For the interpolation algorithm in the resampling process, those skilled in the art can use linear interpolation or cubic spline interpolation. The specific numerical calculations are well-known techniques in the field and will not be elaborated here.

[0065] The morphological verification module calculates the Pearson correlation coefficient between the resampled angular velocity waveform of the current period and the baseline kinematic envelope template. The Pearson correlation coefficient measures the degree of linear correlation between two time series in terms of morphological trends. Mathematically, because the mean of each data point is subtracted from the calculation and the denominator undergoes variance normalization, this calculation method can eliminate amplitude differences caused by the absolute magnitude of limb force exertion, retaining only the similarity characteristics of waveform fluctuations, thus truly reflecting the consistency of the movement's force trajectory. The calculation process uses the following formula:

[0066] ;

[0067] in, The Pearson correlation coefficient for the current motion cycle; The summation symbol; Index the data points; This indicates that the starting value of the data point index is 1; The total number of data points for the baseline kinematic envelope template; Indicates the index of data points From 1 to The data are summed up. The angular velocity waveform of the current period after resampling The value of each data point; This is the average value of the angular velocity waveform for the current period after resampling; The first of the baseline kinematic envelope template The value of each data point; The mean of the data for the baseline kinematic envelope template; Index for motion cycles.

[0068] The morphological school verification module has an internal isomorphism determination threshold. In practical engineering applications, even standard actions will exhibit minute physiological jitter with each execution. Isomorphism determination threshold. Typically configured between 0.75 and 0.90, this is used to intercept severely deformed movements while tolerating normal physiological fluctuations. The morphological schooling module uses the Pearson correlation coefficient of the current movement cycle. Isomorphism threshold Perform numerical comparisons and execute the corresponding branch control logic based on the comparison results.

[0069] When the Pearson correlation coefficient of the current exercise cycle Less than the isomorphism threshold At this point, the morphological verification module determines that morphological distortion has occurred in the current motion cycle. Morphological distortion is usually caused by interference from non-target movements, external collisions to the limbs, or momentary sensor slippage. In this state, the electromyography (EMG) data of the current cycle cannot reflect the true working state of the target muscle. The morphological verification module directly intercepts the data of the current cycle and executes a state preservation strategy. Specifically, the morphological verification module stops further processing of the current cycle data and extracts the final output feature of the previous cycle stored in the system memory, namely the frequency band energy distribution index of the previous cycle, and directly outputs it as the frequency band energy distribution index of the current cycle. Nonlinear phase space reconstruction depends on the continuous evolution of time series features in the topological structure. If abnormal feature values ​​generated by distorted movements are mixed into the time series, it will cause jump-like discontinuities in the phase space trajectory, causing deviations in the calculation of the maximum Lyapunov exponent. The state preservation strategy can mask this type of distortion noise and maintain the continuity of the feature data stream topology.

[0070] When the Pearson correlation coefficient of the current exercise cycle Equal to or greater than the isomorphism threshold At this point, the morphological verification module determines that the current action is valid. The morphological verification module then transmits the current cycle's data stream and the corresponding rhythm scaling factor to the time axis compensation module, initiating subsequent electromechanical delay calculations and time window calibration processes. Through the aforementioned waveform alignment, correlation calculations, and threshold-based branch determination, the system effectively screens and isolates abnormal actions.

[0071] After the data stream of the current motion cycle passes the isomorphic evaluation of the morphological school module, the time axis compensation module receives the surface electromyography signal segment and angular velocity signal segment of that cycle and extracts the electromechanical delay physical features.

[0072] Electromechanical delay is the time difference between a muscle receiving neurophysiological stimulation and generating a physically observable mechanical force. The time axis compensation module locates the true starting point of muscle electrophysiological activity. This module uses the Teager-Kaiser energy operator to perform a nonlinear transformation on the surface electromyography (EMG) signal segment of the current movement cycle. From a signal processing perspective, the Teager-Kaiser energy operator can simultaneously reflect the instantaneous amplitude and frequency of a time series, exhibiting high sensitivity to high-frequency energy jumps generated during muscle contraction, thereby amplifying the initial jump characteristics of the EMG signal envelope. The calculation process uses the following core formula:

[0073] ;

[0074] in, For the Teager-Kaiser energy operator; This is a segment of surface electromyography (EMG) signals during the current movement cycle; The first segment of the surface electromyography signal in the current movement cycle Values ​​of discrete sampling points; The first segment of the surface electromyography signal in the current movement cycle The square of the values ​​of each discrete sampling point; For discrete sampling point index; The first segment of the surface electromyography signal in the current movement cycle Values ​​of discrete sampling points; The first segment of the surface electromyography signal in the current movement cycle The values ​​of discrete sampling points.

[0075] After the aforementioned nonlinear transformation, the time-axis compensation module acquires the instantaneous energy sequence of the electromyographic signal. The time-axis compensation module extracts the resting data window at the beginning of the movement cycle. In practical applications, this resting data window is typically taken as the first 100 to 200 milliseconds before the start of the cycle or at the very beginning of the cycle, when the target muscle is usually in a relaxed state. The time-axis compensation module calculates the mean and standard deviation of the background noise within this resting data window, and adds a specific multiple of the standard deviation to the mean as an adaptive trigger threshold, typically configured as 3 to 5 times. In the instantaneous energy sequence, the time-axis compensation module records the timestamps of the data points that first continuously exceed the adaptive trigger threshold for a preset duration as the neural excitation pacemaker timestamps. The preset duration is typically configured to be between 15 and 25 milliseconds to eliminate interference from single high-frequency noise spikes, ensuring that the extracted commands are genuine and continuous muscle contractions.

[0076] The time-axis compensation module locates the extreme points of mechanical force application in the kinematic dimension. Since limb angular velocity reflects the macroscopic motion process, and the mechanical force generated by muscle contraction directly corresponds to the angular acceleration of limb movement in physics, the time-axis compensation module performs first-order differential calculation on the current cycle angular velocity signal before resampling to obtain the angular acceleration time series of the current motion cycle. For the discrete differential algorithm of the time series, those skilled in the art can use forward differencing or central differencing methods; the specific numerical derivations are well-known techniques in the field and will not be elaborated here.

[0077] The time-axis compensation module iterates through the generated angular acceleration time series to find the point of maximum amplitude within the force-generating phase of the cycle. In periodic reciprocating motion, the system extracts the point of maximum angular acceleration amplitude in the force-generating direction based on the specific motion direction of the current cycle. This extreme point corresponds to the moment when the muscle group generates maximum mechanical explosive force. The time-axis compensation module extracts the timestamp corresponding to this extreme point on the global time axis and records it as the mechanical force-generating extreme point timestamp. .

[0078] The time axis compensation module calculates the timestamps of the extreme points of mechanical force application. With neural excitation pacemaker timestamp The absolute time difference is used to obtain the actual electromechanical delay of the current motion cycle. The calculation process uses the following formula:

[0079] ;

[0080] in, The actual electromechanical delay of the current motion cycle; This is a timestamp for the extreme point of mechanical force exertion. For the time stamp of the neural excitation pacemaker; Index for motion cycles.

[0081] Actual electromechanical delay This records the actual physical time consumed by the muscle to complete the electromechanical coupling transition in the current state. As fatigue gradually accumulates, the efficiency of calcium ion release and reabsorption within muscle fibers decreases, leading to a numerically prolonged trend in the actual electromechanical delay. This parameter directly reflects the current microscopic physiological state of the muscle and, as a control variable for calculating the time axis drift, is passed to the next computational node within the system. Through the aforementioned feature point extraction and time difference calculation, the function of extracting electromechanical delay features is supported by complete lower-level steps.

[0082] After obtaining the actual electromechanical delay of the current motion cycle, the time axis compensation module calculates the time drift and performs feedforward compensation of the time window to correct the phase misalignment caused by dynamic fatigue.

[0083] The time axis compensation module retrieves the reference latency time from the microprocessor's internal cache. The time axis compensation module calculates the actual electromechanical delay of the current motion cycle. Compared with the reference delay time The difference is used to obtain the electromechanical delay time drift for the current cycle. The calculation process uses the following formula:

[0084] ;

[0085] in, This is the difference operator, representing the difference in the change of a physical quantity; The symbol for the basic variable of electromechanical delay time; This represents the electromechanical delay time shift during the current motion cycle. The actual electromechanical delay of the current motion cycle; The base delay time; Index for motion cycles.

[0086] Electromechanical delay time drift Physically, this represents the increased lag between neural control commands and actual mechanical contraction due to accumulated muscle fatigue. As fatigue deepens, the electrical activity of the nervous system advances on the time axis in order to output the same mechanical force as the baseline state. From a physiological perspective, since the generation of electrophysiological signals necessarily precedes mechanical movement, when the electrophysiological delay increases, the sampling window must open earlier on the time axis to capture the true electrophysiological precursor phase that triggers a specific mechanical force exertion. Therefore, this increased lag time needs to be subtracted from the original baseline time point, shifting the sampling interval to the left on the absolute time axis, i.e., towards historical time. If no compensation is performed during signal interception, the electromyographic signal intercepted according to the fixed mechanical movement phase will lag behind the true electrophysiological precursor phase. The time axis compensation module uses this drift amount as a feedforward control parameter to correct the absolute time boundary of the sampling window by shifting it forward.

[0087] The time axis compensation module combines the relative time percentage defined by the data preprocessing module with the actual duration of the current cycle to calculate the adaptively compensated electromyographic signal truncation time boundary on the absolute time axis. The time axis compensation module obtains the start-zero crossover timestamp of the current movement cycle and calculates the base time point under dynamic rhythm scaling based on the relative start time percentage and relative end time percentage. The time axis compensation module subtracts the electromyographic delay time drift from the base time point to obtain the final truncation boundary. In specific calculations, the relative start time percentage... and relative end time percentage All values ​​are represented as decimals between 0 and 1 to ensure consistency of time units. The calculation process uses the following formula:

[0088] ;

[0089] ;

[0090] in, The start time of the compensated electromyographic signal extraction; This is the timestamp of the starting zero-crossing point of the current motion cycle; Percentage of the relative start time; This refers to the actual duration. This represents the electromechanical delay time shift during the current motion cycle. The end time of the compensated electromyographic signal interception; Percentage of relative end time; Index for motion cycles.

[0091] The time axis compensation module extracts the start time based on the compensated electromyographic signal. and the end time of the compensated electromyographic signal interception The system extracts target signal segments from a continuous surface electromyography (EMG) signal data stream. Using the aforementioned compensation mapping mechanism, the system, on a macroscopic scale, maps the actual duration... Changes in motor rhythm were tracked, and electrophysiological changes were measured using electromechanical delay time drift. The electromechanical phase shift induced by fatigue was corrected. The function of adaptively adjusting the sampling window and compensating for phase misalignment was supported by specific lower-level features through this numerical calculation process. The extracted surface electromyography signal fragments were transmitted to the phase space reconstruction module, ensuring the data phase consistency for subsequent nonlinear fatigue feature extraction.

[0092] After the time axis compensation module completes the extraction of the surface electromyography signal segment, the feature extraction module receives the time segment and converts it from the time domain to the frequency domain to extract the frequency band energy distribution index.

[0093] From the perspective of muscle electrophysiology, as muscle fatigue intensifies, the accumulation of metabolic byproducts such as lactic acid leads to a decrease in the intramuscular pH, which in turn slows down the conduction velocity of action potentials on muscle fibers. At the macroscopic signal level, this microscopic physiological change manifests as a general shift in the power spectral density of surface electromyography (EMG) signals from high-frequency to low-frequency regions. To quantify this energy transfer phenomenon, the power spectral density of the compensated signal was calculated.

[0094] In continuous dynamic motion, the actual duration of each motion cycle is constantly changing, and the time axis compensation module performs feedforward translation correction on the absolute time boundary of the sampling window. This results in the number of discrete data points contained in each captured surface electromyography signal segment being unequal. If the data length is not standardized and frequency domain transformation is directly performed on signal segments of different lengths, it will lead to differences in the frequency resolution of each cycle, causing a physical misalignment when performing energy integration within a fixed frequency band.

[0095] To address the aforementioned resolution inconsistency issue, the feature extraction module performs zero-padding on the surface electromyography (EMG) signal segments of the current cycle. The feature extraction module sets the number of analysis points to a power of 2 that is greater than the maximum segment length found in the system. In practical configurations, for electromyography (EMG) front-end acquisition devices with a sampling frequency of 1000 Hz, the number of analysis points... The value is typically set to 1024 or 2048. The feature extraction module continuously adds zero-value data points to the end of the current signal segment until the total number of data points in the sequence equals the set number of analysis points. This zero-padding operation unifies the spectral resolution of all motion cycles to a fixed scale without altering the original signal's frequency components.

[0096] After data alignment, the feature extraction module performs a Discrete Fourier Transform on the zero-padded electromyography (EMG) signal sequence. The calculation process uses the following core formula:

[0097] ;

[0098] in, The first complex number sequence in the frequency domain of the current motion period Values ​​for each frequency point; For discrete sampling point index; The initial value of the discrete sampling point index is 0; To analyze the number of points; Indicates the index of discrete sampling points From 0 to The data are summed up. The first zero-filled surface electromyography signal sequence of the current movement cycle The value of each data point; is the base of the natural logarithm; The imaginary unit; Pi is a constant. This is the index of the discrete frequency sampling point, with a value ranging from 0 to... ; This serves as the index for the motion cycle. For the underlying butterfly operation implementation of fast calculation of the Discrete Fourier Transform, those skilled in the art can use fast algorithms based on radix-2 or radix-4. The specific numerical derivation and array shifting are well-known techniques in the field and will not be elaborated here.

[0099] After obtaining the frequency domain complex sequence, in order to establish the discrete frequency sampling point index The feature extraction module calculates the frequency resolution and derives the physical frequency corresponding to each index by establishing a correspondence between the index and the actual physical frequency. The core formula for calculating this physical frequency involves indexing the discrete frequency sampling points. Multiply by the sampling frequency of the hardware acquisition front end Divide by the number of analysis points This mapping relationship enables the system to address specific real physical frequency bands within a discrete array sequence.

[0100] The feature extraction module calculates the power spectral density of the current motion cycle. Power spectral density reflects the signal energy per unit frequency band; this calculation eliminates interference from differences in the absolute value of total energy caused by variations in the duration of different actions. The calculation process uses the following formula:

[0101] ;

[0102] in, The power spectral density sequence of the current motion period is the first... Values ​​for each frequency point; To analyze the number of points; The sampling frequency of the hardware acquisition front end; The first complex number sequence in the frequency domain of the current motion period Values ​​for each frequency point; This is the index of the discrete frequency sampling point, with a value ranging from 0 to... ; Index for motion cycles.

[0103] After generating the complete power spectral density sequence, the feature extraction module truncates it according to the effective physiological frequency band of the surface electromyography (EMG) signal. The effective energy of the human surface EMG signal is mainly concentrated in the frequency range of 20 Hz to 500 Hz. Frequency bands below 20 Hz are usually mixed with motion artifact noise caused by electrode cable swaying or skin physical deformation; frequency bands above 500 Hz are mainly composed of environmental high-frequency electromagnetic interference and equipment background thermal noise. The feature extraction module retains the numerical points in the power spectral density sequence corresponding to the 20 Hz to 500 Hz frequency range and removes data outside this physical range. Specifically, based on the aforementioned physical frequency mapping relationship, the feature extraction module calculates the starting frequency index corresponding to 20 Hz and the ending frequency index corresponding to 500 Hz, and extracts the array segment in the power spectral density sequence located between these starting and ending frequency indices. This truncation operation filters out noise interference from non-physiological frequency bands, providing a high signal-to-noise ratio frequency domain data foundation for subsequent calculation of the band energy distribution index, thus providing clear lower-level technical support for frequency domain transformation and feature extraction operations.

[0104] After obtaining the truncated effective physiological frequency band power spectral density sequence, the feature extraction module performs frequency band division and energy integration on the sequence to quantify the spectral shift phenomenon during muscle contraction.

[0105] The feature extraction module sets a critical segmentation frequency, dividing the effective physiological frequency band from 20 Hz to 500 Hz into low-frequency and high-frequency bands. Due to the differences in physiological characteristics among different muscle groups, the critical segmentation frequency is usually configured to be between 80 Hz and 100 Hz. The specific value of the critical segmentation frequency is usually determined based on the median frequency of the initial power spectrum of the target muscle group in a non-fatigue state, to ensure that the energy distribution in the high and low frequency bands is at a relatively balanced level in the early stages of exercise. Based on the aforementioned established physical frequency mapping relationship, the feature extraction module calculates the discrete frequency sampling point indices corresponding to 20 Hz, the critical segmentation frequency, and 500 Hz. The specific calculation method is to multiply the target physical frequency by the number of analysis points, divide by the sampling frequency, and round the calculation result to obtain the corresponding integer index.

[0106] The feature extraction module performs discrete summation on the power spectral density sequence in the low-frequency band to obtain the low-frequency energy integral value of the current movement cycle; it also performs discrete summation in the high-frequency band to obtain the high-frequency energy integral value of the current movement cycle. This discrete integration operation is mathematically equivalent to calculating the area under the power spectral density curve in a specific frequency band, and physically represents the total surface electromyographic energy contained within that specific frequency band. For numerical integration algorithms of discrete sequences, those skilled in the art can use the rectangular method or the trapezoidal rule; the specific numerical derivation and program loop implementation are well-known techniques in the field and will not be elaborated upon here.

[0107] The feature extraction module calculates the ratio of the low-frequency band energy integral value to the high-frequency band energy integral value of the current motion cycle to obtain the frequency band energy distribution index of the current motion cycle. Since the numerator and denominator in the ratio calculation both belong to the accumulation of discrete sequences, their frequency resolution (i.e., the physical bandwidth corresponding to a single frequency point) cancels each other out during the division process. Therefore, in the specific algorithm implementation, a simple summation of the power spectral density amplitudes can be directly used to equivalently perform the actual energy integration calculation. The calculation process uses the following core formula:

[0108] ;

[0109] in, This represents the frequency band energy distribution index for the current motion cycle; This is the index of the discrete frequency sampling point, with a value ranging from 0 to... ; This is the starting frequency index corresponding to 20 Hz; This is the index of the segmentation frequency corresponding to the critical segmentation frequency; This is the cutoff frequency index for the low-frequency band. Indicates the index of discrete frequency sampling points from arrive The data in the low-frequency band are summed up. The power spectral density sequence of the current motion period is the first... Values ​​for each frequency point; This is the end frequency index corresponding to 500 Hz; Indicates the index of discrete frequency sampling points from arrive The data in the high-frequency band are summed up. Index for motion cycles.

[0110] Band energy distribution index This study quantifies the trend of electromyographic signal energy shifting from high to low frequencies at both physical and physiological levels. With increasing repetitions and accumulated muscle micro-fatigue, the energy of nerve discharges in the high-frequency band gradually decreases, while the energy in the low-frequency band gradually increases. This leads to an increase in the numerator and a decrease in the denominator, resulting in a monotonically increasing value. This index eliminates amplitude interference caused by individual absolute strength and the intensity of a single exertion, objectively reflecting the degree of decline in muscle function.

[0111] The feature extraction module stores and arranges the frequency band energy distribution index calculated for each consecutive motion cycle in chronological order of the motion cycle occurrence, generating a frequency band energy distribution index time series. This series, as a discrete one-dimensional nonlinear time series, constitutes the basic input data stream for subsequent nonlinear dynamics analysis and phase space reconstruction. Through the above frequency band division, energy integration, and ratio calculation, the function of extracting the frequency band energy distribution index and generating a time series has obtained complete lower-level feature support.

[0112] After receiving the frequency band energy distribution index time series generated by the feature extraction module, the phase space reconstruction module performs high-dimensional spatial projection on the one-dimensional time series to reconstruct the phase space trajectory that can reflect the nonlinear dynamic evolution characteristics inside the muscle system.

[0113] From a kinetic perspective of human physiology, muscle motion control during fatigue accumulation is a nonlinear dynamic system. A simple one-dimensional frequency band energy distribution exponential time series can only reflect the overall transfer trend of muscle energy along the frequency axis, failing to reveal the interactions between the various state variables within the system. The phase space reconstruction module, based on Takens' delayed embedding theorem, maps the one-dimensional time series to a multi-dimensional phase space. Considering the balance between the degrees of freedom characteristics of the muscle fatigue evolution process and the computational resource requirements, the phase space reconstruction module configures the system's embedding dimension to two dimensions.

[0114] The phase space reconstruction module calculates the optimal delay time required to reconstruct the phase space. The setting of the delay time directly determines the quality of the phase space reconstruction. If the delay time is too small, the reconstructed adjacent coordinate components will be too close in value, causing the phase trajectory to be squeezed near the diagonal of the coordinate system, making it impossible to unfold the dynamic characteristics of the system; if the delay time is too large, the adjacent coordinate components will lose their dynamic correlation, causing the reconstructed phase space to exhibit random divergence characteristics. For the solution process of the optimal delay time, those skilled in the art can use the mutual information method or the autocorrelation function method. The specific probability distribution calculation and local minimum threshold determination are well-known techniques in the field and will not be elaborated here. In the actual evaluation of continuous dynamic repetitive motion, this optimal delay time usually corresponds to a data span of 1 to 3 motion cycles.

[0115] The phase space reconstruction module constructs a two-dimensional state vector for the frequency band energy distribution exponential time series based on the acquired optimal delay time. Starting from the current data point in the sequence, the module extracts later data points on the time axis according to the set optimal delay time span, combining both to form the coordinate components of a two-dimensional vector. Since the frequency band energy distribution exponential time series is a discrete sequence based on motion periods, the optimal delay time... In numerical calculations, it is defined as a positive integer representing the number of delayed motion cycles. The calculation process uses the following core formula:

[0116] ;

[0117] in, For the first The two-dimensional phase space state vector corresponding to each motion cycle; This represents the frequency band energy distribution index for the current motion cycle; The frequency band energy distribution index corresponding to the motion period after the optimal delay time is used here as the second coordinate component of the two-dimensional phase space state vector. The optimal delay time for phase space reconstruction represents the number of delayed motion cycles; Index for motion cycles; This is the index for the delayed motion cycle.

[0118] The phase space reconstruction module transforms the original one-dimensional scalar data sequence into a series of two-dimensional phase space state vectors arranged continuously on the time axis by traversing the entire frequency band energy distribution exponential time series. To prevent array out-of-bounds errors, the motion cycle index is used during actual traversal and calculation. The maximum effective value is the total length of the original time series minus the optimal delay time. These state vectors are projected and connected sequentially on a two-dimensional coordinate plane, forming a phase space trajectory describing the evolution of muscle fatigue. As the human musculoskeletal system progresses from a stable, non-fatigued contraction state to a fatigued state, the topological structure of this phase space trajectory expands outward from a clustered state. Regarding the functional overview of phase space reconstruction, the specific data mapping process involving the introduction of time delay calculations and the combination of two-dimensional vectors, as described above, provides underlying technical support with clear mathematical logic and geometric meaning. This reconstruction operation ensures that subsequent calculations of nonlinear dynamic exponents can obtain input conditions with complete topological properties.

[0119] After obtaining the reconstructed two-dimensional phase space state vector sequence, the nonlinear dynamics exponent calculation module calculates the local divergence rate of adjacent trajectories in the phase space and extracts the maximum Lyapunov exponent to evaluate the dynamic bifurcation characteristics of the system.

[0120] The nonlinear dynamics exponent calculation module traverses the reconstructed two-dimensional phase space state vectors. For each selected reference state vector, it searches for the state vector with the closest spatial distance in the entire phase space and takes it as the nearest neighbor state vector. To avoid mistaking vectors with similar temporal distances on the same time trajectory as spatially adjacent points, the nonlinear dynamics exponent calculation module sets a time separation window parameter. This time separation window parameter is typically configured to be greater than the average number of motion periods of the frequency band energy distribution exponent time series, to ensure that the two selected vectors have topological independence in the dynamic space, rather than simple temporal autocorrelation. The nonlinear dynamics exponent calculation module stipulates that the index difference between the reference state vector and the nearest neighbor state vector on the time axis must be greater than this window parameter. The nonlinear dynamics exponent calculation module calculates the initial Euclidean distance between the reference state vector and its nearest neighbor state vector.

[0121] After determining the initial adjacent vector pairs, the nonlinear dynamics exponent calculation module tracks the evolution of these two state vectors in phase space over time. An evolution step size is set, which is a discrete integer index in numerical calculations, representing the number of motion cycles to be shifted backward. Its value starts from 0 and gradually increases to a preset maximum evolution step size. The nonlinear dynamics exponent calculation module simultaneously shifts the reference state vector and the nearest neighbor state vector backward in the time series by the number of cycles corresponding to this step size, and calculates the Euclidean distance between the two new state vectors after the evolutionary shift, i.e., the evolutionary distance. For the nearest neighbor search algorithm in phase space and the underlying numerical calculation of the Euclidean distance, those skilled in the art can use the kd-tree algorithm or spatial grid partitioning method. The specific data structures and distance comparisons are well-known techniques in the field and will not be elaborated upon here.

[0122] To quantify the exponential divergence characteristics of the phase trajectory over time, the nonlinear dynamics exponential calculation module calculates the average of the natural logarithms of the evolution distances of all reference vector pairs, yielding the overall local divergence evolution curve of the system. The calculation process uses the following core formula:

[0123] ;

[0124] in, For the evolution step size is The mean logarithmic divergence value at time; The total number of reference state vectors; Index for motion cycles; The starting value of the motion cycle index is 1; Indicates the index of the motion cycle From 1 to The data are summed up. The operator for natural logarithms; For the first Each reference state vector and its nearest neighbor state vector pass through the phase space. Euclidean distance after one evolutionary step; This is the evolution step size index, and its value ranges from 0 to the preset maximum number of evolution steps.

[0125] Obtain different evolution step sizes Corresponding mean logarithmic divergence value Subsequently, an evolution curve is generated with the evolution step size as the x-axis and the mean logarithmic divergence value as the y-axis. The nonlinear dynamics exponent calculation module extracts the linearly rising region of this evolution curve in the initial evolution stage. From the perspective of dynamic evolution, after the initial linear rise, the divergence curve, constrained by the finite boundary of the phase space, will eventually tend towards a gentle saturation state. In specific calculations, the nonlinear dynamics exponent calculation module calculates the difference slope of the evolution curve, and determines the inflection point where the slope value decreases and the change amplitude is close to zero as the cutoff point of the linear region, thus intercepting the corresponding linearly rising interval. After extracting this interval, the nonlinear dynamics exponent calculation module uses the least squares method to perform a straight line fitting on the data points of this linear region. The slope of the calculated fitted line is the maximum Lyapunov exponent of the system.

[0126] The maximum Lyapunov exponent, in a physical sense, characterizes the exponential divergence rate of adjacent trajectories in phase space of a system. From the perspective of nonlinear dynamics in muscle physiology, when muscles are in a stable, fatigue-free working state, the neuromuscular system maintains a relatively fixed contraction rhythm through the alternating recruitment of muscle fibers. The phase space trajectories remain convergent or diverge very slowly, and the maximum Lyapunov exponent remains at a relatively small value. As the execution of movements leads to the accumulation of metabolic products, the physiological function of the muscle's micro-force-generating units declines, and the nervous system needs to constantly change its motor unit recruitment strategies to compensate for the loss of force. When this compensatory mechanism reaches its limit, the dynamic equilibrium within the system is disrupted, and even small perturbations can cause the phase trajectories to diverge rapidly, altering the system's dynamic characteristics. At this point, the value of the maximum Lyapunov exponent exhibits a nonlinear surge, indicating that the system has undergone dynamic bifurcation and entered an unstable state.

[0127] Regarding the functional generalization of calculating local divergence and extracting the maximum Lyapunov exponent, specific lower-level technical features are obtained through the detailed mathematical calculation process described above, which involves finding the spatial nearest neighbor, tracking the evolutionary distance, and fitting the slope of the logarithmic divergence curve. The extraction of this exponent transforms the microscopic fatigue evolution process within the muscle system into a quantifiable numerical indicator, providing a computational basis for determining the critical bifurcation point of muscle fatigue in the system.

[0128] After obtaining the maximum Lyapunov index of the current movement cycle, the system controller compares it with a threshold to determine the dynamic state of the target muscle group and outputs control intervention commands.

[0129] The system controller retrieves the initial Lyapunov index baseline value stored internally. This baseline value is obtained during the initial movement phase when the test subject is not fatigued. To ensure the validity of the baseline value, the system controller first evaluates the actual duration variance of the first few movement cycles before starting the calculation. After confirming that the continuous movement rhythm is stable and the test subject has entered a standard exertion state, the system records the maximum Lyapunov index of the first few movement cycles (e.g., the first 5 cycles) and calculates their arithmetic mean, using this average as the baseline value. The system controller calculates the bifurcation determination threshold based on the baseline value. The calculation process uses the following core formula:

[0130] ;

[0131] in, The threshold for determining bifurcation; This serves as the initial Lyapunov index baseline value; The scaling factor margin is set to quantify the critical point at which the system abruptly changes from stable convergence to divergence. The value is typically set between 1.5 and 2.0. This parameter is based on the fact that when human muscles enter a state of deep fatigue, the neural compensation mechanism reaches its limit, and the local divergence rate in phase space usually rises to more than 1.5 times that of the initial steady state.

[0132] The system controller compares the maximum Lyapunov exponent of the current motion cycle with the bifurcation determination threshold. The judgment logic uses the following core inequality:

[0133] ;

[0134] in, The maximum Lyapunov index during the current motion cycle; The threshold for determining bifurcation; Index for motion cycles.

[0135] To prevent misjudgments caused by sudden changes in a single electrophysiological signal or computational noise, the system controller employs a continuous counting mechanism for status confirmation. The system controller determines the maximum Lyapunov exponent of the current motion cycle. Is it greater than or equal to the bifurcation determination threshold? If the value is greater than or equal to the threshold, the system controller increments the internal anomaly counter by one; if it is less than the threshold, the system controller resets the anomaly counter to zero. The system controller then determines whether the current value of the anomaly counter has reached the preset number of fault tolerance cycles. To balance real-time monitoring with anti-interference capability, this number of fault tolerance cycles is typically set to 3 to 5 consecutive motion cycles.

[0136] When the fault counter reaches the number of fault tolerance cycles, the system controller determines that a dynamic topological bifurcation has occurred in the target muscle system, and the muscle has entered an unstable fatigue state. Based on this, the system controller generates a control intervention command and sends it to the peripheral hardware device via the communication interface. In specific engineering implementations, the system controller encapsulates the control intervention command into a standard data frame format and sends it to the actuator via the Controller Area Network (CAN bus) or Serial Peripheral Interface (SPI). Specifically, the system controller sends a resistance adjustment command to the rehabilitation exoskeleton device via the industrial communication bus, instructing the underlying motor driver to reduce damping torque or increase auxiliary thrust to protect the joints from overload damage caused by muscle exhaustion; the system controller sends a status update data packet to the front-end graphical user interface, changing the fatigue status indicator on the screen from green to red and driving a buzzer to sound, reminding the operator to terminate the current training; the system controller sends a shutdown command to the functional electrical stimulation device, cutting off the pulse output of the external stimulation electrodes to avoid inducing muscle rigidity and spasms.

[0137] Please see the appendix Figure 3 To verify the technical effect of the present invention, the system underwent continuous dynamic motion fatigue assessment test in a controlled physical experimental environment.

[0138] In the experimental setup, subjects wore a lower limb rehabilitation exoskeleton device with surface electromyography (EMG) sensors attached to the surface of the vastus lateralis muscle. They continuously performed constant-load knee extension and flexion movements at a set rhythm until subjective exhaustion prevented them from maintaining the intended movement. The sampling frequency of the EMG sensors was configured at 1000 Hz. The system simultaneously recorded and calculated the traditional median frequency (MDF) index and the maximum Lyapunov exponent derived from the reconstructed phase space based on the band energy distribution index of this invention.

[0139] From the appendix Figure 3Subplot (a), the evolution trend of the median frequency (MDF) in the traditional power spectrum, shows that the MDF calculated value for each motion cycle, represented by the solid line, exhibits an overall decreasing trend throughout the entire motion process. The dashed line represents the decreasing trend of the calculated linear fit. This linear index is highly susceptible to interference from changes in muscle length, relative electrode displacement, and fluctuations in force exertion during continuous dynamic contraction, leading to severe oscillations in the median frequency calculated from adjacent motion cycles. Because this curve exhibits a progressive linear decay with high variance noise, it is difficult for the system to determine a clear and objective fatigue critical point on the time axis. Forcibly setting a fixed value as the threshold can easily lead to premature false triggering or delayed judgment by the system.

[0140] From the appendix Figure 3 Subgraph (b), which is the evolution trend diagram of the maximum Lyapunov exponent of this invention, shows the technical advantages of this invention. In the initial and middle stages of the motion (i.e., the motion cycle index represented by the horizontal axis in the diagram)... Within the range of 0 to 60, before the subject's muscles enter an unstable state, the maximum Lyapunov exponent curve, represented by the solid line in this invention, fluctuates within a low and stable numerical range. The system extracts the arithmetic mean of the first 5 effective calculation cycles of this stable phase as the initial baseline value, and multiplies it by a proportionality margin of 1.5 to generate the bifurcation determination threshold. (i.e., the horizontal dotted line in the diagram).

[0141] As the number of exercise cycles increases, the muscle's micro-neural compensatory mechanisms reach their limit. Around the 65th exercise cycle (as shown in the figure), the maximum Lyapunov exponent curve breaks through its previous plateau, exhibiting a non-linear, abrupt increase. When the curve exceeds the bifurcation threshold for three consecutive exercise cycles, the system's internal anomaly counter reaches the tolerance cycle limit. At the control command trigger point marked by the black dot and vertical auxiliary dashed line (as shown in the figure), the system determines that the target muscle group has undergone a dynamic topological bifurcation. At this point, the system controller immediately issues a control intervention command. Upon receiving the command, the exoskeleton device reduces the motor damping torque, effectively preventing compensatory joint sprains caused by complete muscle exhaustion in the subject.

[0142] The comparison results demonstrate that the nonlinear dynamic index calculation method based on phase space reconstruction provided by this invention transforms the fuzzy progressive fatigue process in continuous dynamic motion into a system state bifurcation determination with clear mathematical boundaries.

Claims

1. A method for detecting muscle fatigue threshold based on surface electromyography signals, characterized in that, Includes the following steps: Acquire surface electromyography (EMG) signal sequences and angular velocity signal sequences, and segment them to generate motion cycle sequences; Based on the angular velocity data of the initial number of motion cycles in the motion cycle sequence, a baseline kinematic envelope template is established to define the initial kinematic phase matching window and the reference delay time. Extract the surface electromyography signal segment and angular velocity signal segment of the current movement cycle, perform morphological calibration between the angular velocity signal segment and the baseline kinematic envelope template, and calculate the actual electromechanical delay after the morphological calibration is passed. The electromechanical delay time drift is obtained by subtracting the actual electromechanical delay from the reference delay time, and the electromyographic signal interception time boundary is defined by combining the initial kinematic phase matching window. Within the time boundary of the electromyography signal interception, the surface electromyography signal segment is intercepted to calculate the frequency band energy distribution index and generate a frequency band energy distribution index time series. The frequency band energy distribution index time series is reconstructed into a two-dimensional phase space state vector, the maximum Lyapunov exponent is calculated, and a control intervention command is output when the maximum Lyapunov exponent satisfies the continuous counting mechanism.

2. The method for detecting muscle fatigue threshold based on surface electromyography signals according to claim 1, characterized in that, The process of acquiring surface electromyography (EMG) signal sequences and angular velocity (Angular velocity) signal sequences and segmenting them to generate motion cycle sequences specifically includes: The synchronously acquired surface electromyography signal sequence and angular velocity signal sequence were filtered and preprocessed. By iterating through the data points of the filtered angular velocity signal sequence in time series, the moment when the signal value crosses the zero level is extracted and marked as the zero crossing point; Two adjacent in-phase zero-crossing points on the time axis are used as the start and end boundaries of the motion interval. The start and end zero-crossing point timestamps of the motion interval are read, and the filtered surface electromyography (EMG) signal sequence and angular velocity signal sequence are extracted within the start and end boundaries to obtain the corresponding EMG signal segments and angular velocity signal segments. Calculate the difference between the end zero-crossing point timestamp and the start zero-crossing point timestamp to obtain the corresponding actual duration; The surface electromyography signal segment containing timestamp information and the angular velocity signal segment are bound and stored with the corresponding period index and the actual duration, thereby cutting the filtered signal sequence into discrete motion cycle sequences.

3. The method for detecting muscle fatigue threshold based on surface electromyography signals according to claim 2, characterized in that, The process of establishing a baseline kinematic envelope template based on the angular velocity data of an initial number of motion cycles in the motion cycle sequence, defining the initial kinematic phase matching window and the reference delay time, specifically includes: Extract angular velocity signal segments and surface electromyography signal segments within the initial number of consecutive motion cycles in the motion cycle sequence as sample data; The angular velocity signal sequence within the initial number of motion cycles is normalized on the time axis. The average waveform amplitude of the angular velocity signal segment of the initial number of motion cycles is calculated at each relative time node to generate the baseline kinematic envelope template characterizing the force exertion law of the target action. Extract the peak angular velocity of the baseline kinematic envelope template, find continuous data segments whose amplitude is between the value obtained by multiplying the peak angular velocity by a set lower amplitude ratio coefficient and the value obtained by multiplying the peak angular velocity by a set upper amplitude ratio coefficient, establish the interval covered by the continuous data segments on the relative period progress axis as the initial kinematic phase matching window, and record the relative start time percentage and relative end time percentage of the initial kinematic phase matching window on the normalized time axis; The starting point of the surface electromyography (SEMG) activity and the mechanical force extreme point of the corresponding angular velocity signal segment of the initial number of motion cycles are extracted respectively. The absolute time difference between the starting point of the SEMG activity and the mechanical force extreme point is calculated. The arithmetic mean of the measured absolute time difference is obtained and the arithmetic mean is used as the reference delay time.

4. The method for detecting muscle fatigue threshold based on surface electromyography signals according to claim 3, characterized in that, The step of performing morphological calibration of the angular velocity signal segment with the baseline kinematic envelope template specifically includes: The angular velocity waveform of the current period is resampled using an interpolation algorithm so that the data length of the resampled angular velocity waveform of the current period is the same as the data length of the baseline kinematic envelope template. Calculate the Pearson correlation coefficient between the angular velocity waveform of the current period after resampling and the baseline kinematic envelope template; When the Pearson correlation coefficient is less than the set isomorphism determination threshold, the frequency band energy distribution index of the previous motion cycle is read and the frequency band energy distribution index of the previous motion cycle is directly output as the frequency band energy distribution index of the current motion cycle. When the Pearson correlation coefficient is equal to or greater than the isomorphism determination threshold, the current action is determined to be valid and the verification is passed.

5. The method for detecting muscle fatigue threshold based on surface electromyography signals according to claim 4, characterized in that, After the morphological verification is passed, the actual electromechanical delay is calculated, specifically including: The Teager-Kaiser energy operator is used to perform nonlinear transformation on the surface electromyography signal segment of the current exercise cycle to obtain the instantaneous energy sequence. The resting data window of the initial stage is extracted, and the mean and standard deviation of the background noise in the resting data window are calculated. The mean plus the standard deviation by a set multiple is used as the adaptive trigger threshold. In the instantaneous energy sequence, find the data point that first continuously exceeds the adaptive trigger threshold and lasts for a set duration, and record the timestamp corresponding to the data point as the neural excitation pacing point timestamp; The first-order differential of the current period angular velocity signal before resampling is used to obtain the angular acceleration time series. The timestamp corresponding to the maximum angular acceleration amplitude point in the angular acceleration time series is extracted and recorded as the timestamp of the mechanical force extreme point. Calculate the absolute time difference between the timestamp of the mechanical force excitation extreme point and the timestamp of the neural excitation pacing point to obtain the actual electromechanical delay of the current motion cycle.

6. The method for detecting muscle fatigue threshold based on surface electromyography signals according to claim 5, characterized in that, The step of obtaining the electromechanical delay time drift by subtracting the actual electromechanical delay from the reference delay time, and defining the electromyographic signal interception time boundary in conjunction with the initial kinematic phase matching window, specifically includes: Calculate the difference between the actual electromechanical delay of the current motion cycle and the reference delay time to obtain the electromechanical delay time drift. Obtain the starting zero-crossing point timestamp and actual duration of the current exercise cycle, calculate the product of the actual duration and the percentage of the relative starting time, and add it to the starting zero-crossing point timestamp to obtain the starting base time point; Calculate the product of the actual duration and the percentage of the relative end time, and add it to the start zero-crossing timestamp to obtain the end base time point; The electromyographic delay time drift is subtracted from the starting base time point and the ending base time point respectively to obtain the compensated electromyographic signal capture start time and compensated electromyographic signal capture end time, thereby defining the electromyographic signal capture time boundary.

7. The method for detecting muscle fatigue threshold based on surface electromyography signals according to claim 6, characterized in that, The step of extracting surface electromyographic signal segments within the time boundary of the electromyographic signal extraction, calculating the frequency band energy distribution index, and generating a time series of the frequency band energy distribution index specifically includes: Zero-filling is performed on the truncated surface electromyography signal fragments until the total number of data points equals the number of analysis points, where the number of analysis points is a power of 2 greater than the maximum fragment length of the signal sequence. The zero-padding electromyography signal sequence is subjected to discrete Fourier transform to obtain a frequency domain complex sequence for calculating the power spectral density. The numerical points corresponding to the frequency range of 20 Hz to 500 Hz are retained to generate a power spectral density sequence. Extract the surface electromyography signal segments of the initial number of exercise cycles, calculate the median frequency of the initial power spectrum as the critical segmentation frequency, and divide the frequency range of 20 Hz to 500 Hz into a low-frequency band and a high-frequency band according to the critical segmentation frequency. Discretely accumulate and sum the power spectral density sequence in the low-frequency band and the high-frequency band respectively to obtain the low-frequency band energy integral value and the high-frequency band energy integral value of the current motion cycle; The frequency band energy distribution index is obtained by calculating the ratio of the low-frequency band energy integral value to the high-frequency band energy integral value, and the frequency band energy distribution index of each consecutive motion cycle is arranged in chronological order to generate the frequency band energy distribution index time series.

8. The method for detecting muscle fatigue threshold based on surface electromyography signals according to claim 7, characterized in that, The step of reconstructing the frequency band energy distribution index time series into a two-dimensional phase space state vector specifically includes: The optimal delay time required for phase space reconstruction of the frequency band energy distribution index time series is calculated using the mutual information method. Select the current data point from the frequency band energy distribution index time series, and extract the lag data point on the time axis that differs from the current data point by the optimal delay time span; The coordinate components of a two-dimensional vector are formed by combining the current data point and the lagged data point. By traversing the time series, the original one-dimensional scalar data sequence is transformed into a two-dimensional phase space state vector that is continuously arranged on the time axis.

9. The method for detecting muscle fatigue threshold based on surface electromyography signals according to claim 8, characterized in that, The calculation of the maximum Lyapunov exponent specifically includes: Traverse the reconstructed two-dimensional phase space state vectors. For the selected reference state vector, find the nearest neighbor state vector in the phase space that is spatially closest to the reference state vector, and limit the index difference between the reference state vector and the nearest neighbor state vector on the time axis to be greater than the set time separation window parameter. Obtain the evolution step size, which is gradually increased from zero to the set maximum evolution step size. Simultaneously shift the reference state vector and the nearest neighbor state vector backward in the time series by the number of cycles corresponding to the evolution step size. Calculate the Euclidean distance between the two new state vectors after the evolution shift and obtain the evolution distance. Calculate the average natural logarithm of the evolution distance for all reference state vector pairs, obtain the average logarithmic divergence value corresponding to different evolution step sizes, and generate a local divergence evolution curve with the evolution step size on the horizontal axis and the average logarithmic divergence value on the vertical axis. The linearly rising region of the local divergence evolution curve in the initial evolution stage is extracted. The least squares method is used to fit the data points in the linearly rising region to a straight line, and the slope of the fitted line is calculated. The slope is used as the maximum Lyapunov exponent.

10. The method for detecting muscle fatigue threshold based on surface electromyography signals according to claim 9, characterized in that, The step of outputting a control intervention command when the maximum Lyapunov exponent satisfies the continuous counting mechanism specifically includes: Calculate the maximum Lyapunov exponent for each of the initial number of motion cycles, and use the arithmetic mean as the initial Lyapunov exponent benchmark value. Multiply the initial Lyapunov exponent benchmark value by a set proportional coefficient margin to generate a bifurcation determination threshold. Set an anomaly counter and compare the maximum Lyapunov exponent of the current motion cycle with the bifurcation determination threshold. If the maximum Lyapunov exponent is greater than or equal to the bifurcation determination threshold, the abnormality counter value is incremented by one; If the value is less than the bifurcation determination threshold, the anomaly counter will be cleared to zero. Determine whether the current value of the anomaly counter has reached the set number of fault tolerance cycles. When the number of fault tolerance cycles is reached, it is determined that the continuous counting mechanism is satisfied, and the target muscle tissue is determined to have undergone dynamic topological bifurcation, thereby generating the control intervention command. The control intervention command is encapsulated into a standard data frame format and sent to a receiving terminal device through a communication interface. The receiving terminal device can be a rehabilitation exoskeleton device, a front-end graphical user interface device, or a functional electrical stimulation device. The control intervention command sent to the rehabilitation exoskeleton device is a resistance adjustment command, the control intervention command sent to the front-end graphical user interface device is a status update data packet, and the control intervention command sent to the functional electrical stimulation device is a shutdown command.

Images