Rapid gradient convolution kernel compensation surface electromyogram signal decomposition method based on stripping

Abstract

The invention provides a rapid gradient convolution kernel compensation surface electromyogram signal decomposition method based on stripping, and aims to improve the decomposition efficiency of high-density surface electromyogram signals. The method comprises the following steps of: after applying fast gradient convolution kernel compensation each time, acquiring an output motion unit issuing time sequence, and evaluating whether the sequence meets a preset standard or not; estimating a corresponding motion unit action potential through the identified motion unit distribution time sequence, and subtracting the motion unit action potential from the original electromyographic signal to generate a residual signal; then, the fast gradient convolution kernel compensation algorithm is applied to the residual signal again for further decomposition; and the process is continuously iterated until an effective motion unit issuing time sequence cannot be identified. According to the invention, the number of identified motion unit issuing time sequences can be increased.

Description

Technical Field

[0001] The present invention belongs to the field of biological signal processing, and relates to a method for decomposing surface electromyogram signals based on peeling-based fast gradient convolution kernel compensation. Background Art

[0002] Electromyogram (EMG) signals are electrical signals generated during muscle contraction or relaxation, and are widely used in fields such as physiology, medicine, sports science, and prosthetic control. By decomposing EMG signals, motor unit action potentials (MUAPs) and their firing time series can be extracted, providing important bases for disease diagnosis and evaluation of nerve dysfunction.

[0003] Surface electromyogram (sEMG) signals have received extensive attention due to their non-invasiveness. Usually, blind source separation (BSS) methods, such as convolution kernel compensation (CKC) or gradient convolution kernel compensation (gCKC), are used to decompose signals to obtain motor unit information. Fast gradient convolution kernel compensation (fgCKC) improves the optimization method of gCKC, and uses the information of past gradients and squared gradients to accelerate the convergence process of the current gradient, thereby significantly improving the decomposition efficiency. However, existing methods are difficult to extract a sufficient number of motor units in complex signal environments, especially for motor unit information with low signal amplitude or large noise influence, and the decomposition ability is limited. Summary of the Invention

[0004] In view of the deficiencies of the prior art, the present invention provides a method for decomposing surface electromyogram signals based on peeling-based fast gradient convolution kernel compensation. By introducing a peeling strategy, the identified motor unit signals are gradually subtracted to generate a simpler residual signal, thereby significantly increasing the number of decomposed motor units.

[0005] To achieve the above object, the technical means adopted by the present invention are as follows:

[0006] A method for decomposing surface electromyogram signals based on peeling-based fast gradient convolution kernel compensation, comprising the following steps:

[0007] Step 1: Mathematically model multi-channel surface electromyogram signals, and expand the modeled multi-channel surface electromyogram signals based on a preset expansion factor to obtain expanded surface electromyogram signals; obtain the firing time series and noise series corresponding to the expanded electromyogram signals;

[0008] Step 2: Based on the convolution mixture model of surface electromyogram signals, construct a cross-correlation matrix of surface electromyogram signals, and initialize the activity index through the cross-correlation matrix;

[0009] Step 3: Initialize the cross-correlation vector based on the activity index, set the gradient function and learning rate, introduce an exponentially weighted moving average model, and perform offline calculation and iterative update of the cross-correlation vector;

[0010] Step 4: Construct a sliding window for the original multi-channel surface electromyogram (EMG) signal, extract the surface EMG signal in the sliding window and expand it, calculate the expanded cross-correlation matrix, estimate the firing time series of each motor unit using the expanded cross-correlation matrix, and adjust the firing time series according to the updated cross-correlation vector. Repeat Step 4 until all signals are decomposed;

[0011] Step 5: Compare the firing time series with a preset threshold. If the coefficient of variation of the spike intervals or the average firing rate of the firing time series does not meet the requirements, discard it, stop peeling, jump to Step 8; otherwise, continue peeling and proceed to Step 6;

[0012] Step 6: Based on the firing time series, calculate the action potential of each motor unit, subtract the action potential from the original EMG signal to generate a residual signal, and re-enter the residual signal into Step 4 for the next round of decomposition;

[0013] Step 7: Repeat Steps 5 - 6, continue to extract the sliding window signal, estimate the firing time series, and peel off the motor unit components from the residual signal until no effective motor unit firing time series can be identified;

[0014] Step 8: Complete the decomposition of all signals and finally output the effective motor unit firing time series.

[0015] Further, Step 1 specifically includes:

[0016] Step 1.1: Represent the EMG signal obtained from M channels generated by the activities of N motor units as:

[0017]

[0018] where, x i (n) represents the nth sample of the ith channel, h ij is the action potential of the jth motor unit with length L in the ith channel, and s j (n) represents the firing time series of the jth motor unit;

[0019] Add noise to the electromechanical signal and represent the EMG signal as:

[0020] y i (n) = x i (n) + ω i (n)

[0021] where, y i (n) represents the EMG signal containing noise, and ω i (n) represents the noise in the signal;

[0022] Step 1.2: Expand the EMG signal based on a preset expansion factor of K-1, and the obtained EMG signal is expressed as:

[0023]

[0024] Then the corresponding firing time series and noise are respectively:

[0025]

[0026] Step 1.3: The surface EMG signal convolution mixture model after expansion is expressed as where H is the mixing matrix:

[0027]

[0028] where, h ij represents the action potential of the jth motor unit with length L in the ith channel.

[0029] Furthermore, Step 2 specifically includes:

[0030] Construct a cross-correlation matrix based on the convolution mixture model of the surface EMG signal:

[0031]

[0032] where, E() represents calculating the mathematical expectation;

[0033] Initialize the activity index:

[0034]

[0035] Furthermore, Step 3 specifically includes:

[0036] Initialize the cross-correlation vector based on the activity index:

[0037] n 0 = maxarg n (γ(n))

[0038] γ(n 0 ) = 0

[0039]

[0040] where, represents the initial cross-correlation vector, n 0 represents the index when γ(n) reaches the maximum value, represents the index n 0 at which the signal after expansion

[0041] Set the gradient function and learning rate:

[0042] Set as the initial gradient function, and at the same time select the learning rate η = 0.01,

[0043] Introduce the exponentially weighted moving average model:

[0044] y t =(1 - λ)x t +(1 - λ)·λx t-1 +(1 - λ)·λ 2 ×x t-2 +L+(1 - λ)·λ t-1 x 1 +λ t y 0

[0045] where λ is a weight variable with a value ranging from 0 to 1, t is the current time, and y t is the gradient weighted average at the current time t, and x t is the gradient at the current time t;

[0046] Introduce the exponentially weighted moving average model and recalculate the gradient function as:

[0047]

[0048] where, is the firing time series of the jth motion unit at the mth sample, λ and β are weight variables with values ranging from 0 to 1, k is the number of iterations. When k = 1, initialize the gradient

[0049] The update rule for the cross - correlation vector in the offline decomposition process is:

[0050]

[0051] where k is the number of iterations, is the cross - correlation vector at the (k + 1)th time, is the cross - correlation vector at the kth time, η represents the learning rate at the kth time, ε = 10 -8 , represents the mth sample of the surface electromyogram signal.

[0052] Furthermore, in step 4, the cross - correlation vector and the cross - correlation matrix are updated according to the following formula:

[0053]

[0054] where Ψ j ={n 1 ,n 2,...} is the firing time series decomposed from the cross-correlation vector calculated according to the offline process in the sliding window data, l is the learning rate, is the cross-correlation matrix, is the cross-correlation vector, represents the surface electromyogram signal after delay spread within the sliding window.

[0055] Furthermore, in step 5, the coefficient of variation of the interspike interval is calculated according to the following formula:

[0056]

[0057] where SD(ISI) is the standard deviation of the interspike interval, μ(ISI) is the mean of the interspike interval, and ISI represents the time interval between two consecutive discharge events;

[0058] The average discharge rate is calculated according to the following formula:

[0059]

[0060] where N is the number of discharges and T is the observation time.

[0061] Furthermore, step 6 specifically includes:

[0062] Step 6.1: Calculate the action potential of each motor unit using the least mean square error method:

[0063]

[0064] where, is the estimated motor unit action potential matrix, S is the extended peak column matrix, X is the multi-channel sEMG signal matrix, S T is the transpose of matrix S, S T S is the covariance matrix of the peak column, S T X is the correlation matrix between the signal and the peak column;

[0065] Step 6.2: Calculate the residual signal using the reconstructed signal:

[0066] First, calculate the reconstructed signal of the identified motor unit:

[0067]

[0068] where, S is the extended peak column matrix, is the estimated motor unit action potential matrix,

[0069] Subtract the reconstructed signal from the original signal to obtain the residual signal:

[0070]

[0071] Among them, X residual is the residual signal matrix for the next iteration.

[0072] Compared with the prior art, the present invention has the following advantages:

[0073] A fast gradient convolution kernel compensation surface electromyogram signal decomposition method combining a peeling strategy provided by the present invention generates a simpler residual signal by introducing a peeling strategy and gradually subtracting the decomposed motor unit signals, thereby significantly improving the subsequent decomposition ability. The number of decomposed motor units is effectively increased. On the basis of maintaining the decomposition efficiency and accuracy of the fgCKC method, the present invention further expands the decomposition ability and is more suitable for applications in complex signal environments. BRIEF DESCRIPTION OF THE DRAWINGS

[0074] In order to more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the following will briefly introduce the drawings required for the description of the embodiments or the prior art. Obviously, the drawings in the following description are some embodiments of the present invention. For those of ordinary skill in the art, other drawings can be obtained based on these drawings without creative efforts.

[0075] Figure 1 It is a flowchart of a peeling-based fast gradient convolution kernel compensation surface electromyogram signal decomposition method in an embodiment of the present invention.

[0076] Figure 2 It is a decomposition result diagram of signals with different signal-to-noise ratios and excitation levels in an embodiment of the present invention.

[0077] Figure 3 It is a decomposition result diagram of signals with different expansion factors and excitation levels in an embodiment of the present invention. DETAILED DESCRIPTION OF THE EMBODIMENTS

[0078] In order to enable those skilled in the art to better understand the solution of the present invention, the following will clearly and completely describe the technical solutions in the embodiments of the present invention with reference to the drawings in the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, rather than all the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those of ordinary skill in the art without creative efforts shall fall within the protection scope of the present invention.

[0079] As Figure 1 shown, the present invention provides a peeling-based fast gradient convolution kernel compensation surface electromyogram signal decomposition method, including the following steps:

[0080] Step 1. Mathematically model the multi-channel surface electromyogram (sEMG) signals and take the expansion factor as K-1 to expand the EMG signals. Generally, K is taken as 10.

[0081]

[0082] Among them, x i (n) represents the nth sample of the ith channel, and h ij is the action potential of the jth motor unit with length L in the ith channel, and s j (n) represents the firing time series of the jth motor unit;

[0083] Considering the noise in the EMG signals, the above EMG signals can be expressed as:

[0084] y i (n) = x i (n) + ω i (n)

[0085] Among them, y i (n) represents the EMG signal containing noise, and ω i (n) represents the noise in the signal.

[0086] Expanded EMG signals:

[0087]

[0088] Then the corresponding firing time series and noise are respectively:

[0089]

[0090] The convolution mixture model of the expanded surface EMG signals can be expressed as where H is the mixing matrix:

[0091]

[0092] Among them, h ij still represents the action potential of the jth motor unit with length L in the ith channel.

[0093] Step 2. Calculate the cross-correlation matrix and initialize the activity index where E() represents calculating the mathematical expectation.

[0094] Step 3. Initialize the cross-correlation vector and set the gradient function and learning rate. Introduce the exponentially weighted moving average model, perform offline calculation and iteratively update the cross-correlation vector.

[0095] The cross-correlation vector can be initialized as:

[0096] n0 = maxarg n (γ(n))

[0097] γ(n 0 ) = 0

[0098]

[0099] wherein, represents the initial cross - correlation vector, and n 0 represents the index when γ(n) reaches the maximum value, represents the index n 0 when the extended signal. The activity index with subscript n 0 is set to 0 to prevent this value from being reused.

[0100] Set the initial gradient function and the learning rate η = 0.01;

[0101] Introduce the exponentially weighted moving average model, calculate offline and iteratively update the cross - correlation vector

[0102] y t = (1 - λ)x t +(1 - λ)·λx t-1 +(1 - λ)·λ 2 ×x t-2 +L+(1 - λ)·λ t-1 x 1 +λ t y 0

[0103] where λ is a weight variable with a value ranging from 0 to 1. At the current time t, the value of y t is affected by both x t and y t-1 . This formula is called the exponentially weighted moving average model of the stochastic process y;

[0104] The gradient function recalculated by introducing the exponentially weighted moving average model is:

[0105]

[0106] wherein, is the spike train of the j - th motor unit at the m - th sample, λ and β are weight variables with values ranging from 0 to 1, k is the number of iterations. When k = 1, initialize the gradient

[0107] The update rule of the cross - correlation vector in the offline decomposition process is:

[0108]

[0109] where k is the number of iterations, is the cross-correlation vector at the (k + 1)-th time, is the cross-correlation vector at the k-th time, η represents the learning rate at the k-th time, and ε = 10 -8 , represents the m-th sample of the extended EMG signal.

[0110] Step 4. Construct a sliding window for the original multi-channel surface EMG signal, extract the EMG signal in the sliding window and expand it into Update the cross-correlation matrix using the formula and estimate the firing time series of the j-th motor unit Combine the cross-correlation vector updated in Step 3 with the cross-correlation vector updated in this step and the cross-correlation matrix

[0111] The cross-correlation vector and the cross-correlation matrix are updated according to the following formula:

[0112]

[0113] where Ψ j ={n 1 ,n 2 ,...} is the firing time series decomposed from the cross-correlation vector calculated in the offline process in the sliding window data, l is the learning rate, taking 0.1, is the cross-correlation matrix, is the cross-correlation vector, represents the surface EMG signal after delay expansion within the sliding window.

[0114] Step 5. Verify the validity of the firing time series:

[0115] Calculate the coefficient of variation of the interspike interval:

[0116]

[0117] where SD(ISI) is the standard deviation of the interspike interval, μ(ISI) is the mean of the interspike interval, and ISI represents the time interval between two consecutive discharge events (spikes or pulses).

[0118] Calculate the average discharge rate:

[0119]

[0120] where N is the number of discharges and T is the observation time.

[0121] It is required that CoV ISI≤0.3 and ADR ≤ 35 Hz, keep this time series, otherwise discard it.

[0122] Step 6. Introduce a peeling strategy to generate a residual signal:

[0123] Calculate the motor unit action potential (MUAP):

[0124]

[0125] where is the estimated motor unit action potential matrix, S is the extended peak column matrix, X is the multi-channel sEMG signal matrix, S T is the transpose of matrix S, S T S is the covariance matrix of the peak columns, S T X is the correlation matrix between the signal and the peak columns.

[0126] Reconstruct the identified signal:

[0127]

[0128] where S is the extended peak column matrix, is the estimated motor unit action potential matrix.

[0129] Subtract the reconstructed signal from the original signal to generate a residual signal:

[0130]

[0131] where X residual is the residual signal matrix for the next iteration.

[0132] Step 7. Iteratively decompose the residual signal

[0133] Use the residual signal as the input and re-enter the fgCKC algorithm. Repeat Step 5 and Step 6 to gradually peel off the motor unit components. This process is iterated until no valid firing time series can be identified.

[0134] Step 8. Output the valid motor unit information

[0135] After all signal decompositions are completed, output the final motor unit spike train (MUST), significantly increasing the number of decomposed MUSTs.

[0136] The present invention decomposes from three dimensions: signal-to-noise ratio (SNR), excitation level (MVC), and expansion factor (R), and compares and strips the number results of MUST decomposed by fgckc and fgckc. Two sets of comparative experiments are mainly carried out: for simulated signals with an expansion factor of 10, signal-to-noise ratios (SNR) of 5 dB, 10 dB, 20 dB, 30 dB, and excitation levels (MVC) of 10%, 30%, 50%, 25 experiments are carried out for each group, with a total of 4×3×25 experiments; for simulated signals with a signal-to-noise ratio (SNR) of 20 dB, expansion factors (R) of 4, 6, 10, 16, and excitation levels (MVC) of 10%, 30%, 50%, 25 experiments are carried out for each group, with a total of 4×3×25 experiments.

[0137] The decomposition results of signals with different signal-to-noise ratios and excitation levels are as Figure 2 shown, and the decomposition results of signals with different expansion factors and excitation levels are as Figure 3 shown. The number results of stripping MUST decomposed by fgckc and fgckc show that a fast gradient convolution kernel compensation surface electromyogram signal decomposition method combining a stripping strategy proposed by the present invention can decompose more motor unit discharge time series, with the average decomposition number increased by 20%, and basically maintaining the same accuracy as fgckc.

[0138] The present invention proposes a surface electromyogram signal stripping and decomposition method based on fast gradient convolution kernel compensation. After each application of fast gradient convolution kernel compensation, the output motor unit discharge time series is obtained, and it is evaluated whether the series meets a predetermined standard; the corresponding motor unit action potential is estimated through the identified motor unit discharge time series, and it is subtracted from the original electromyogram signal to generate a residual signal; subsequently, the fast gradient convolution kernel compensation algorithm is applied to the residual signal again for further decomposition; this process is iterated continuously until no effective motor unit discharge time series can be identified. The present invention can increase the number of identified motor unit discharge time series.

[0139] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention, rather than to limit them; although the present invention has been described in detail with reference to the foregoing embodiments, those of ordinary skill in the art should understand that they can still modify the technical solutions described in the foregoing embodiments, or perform equivalent replacements for some or all of the technical features; and these modifications or replacements do not make the essence of the corresponding technical solutions deviate from the scope of the technical solutions of the embodiments of the present invention.

Claims

1. A fast gradient convolution kernel compensated surface electromyography signal decomposition method based on peeling, characterized in that: The following steps are involved: Step 1, mathematically modeling the multi-channel surface electromyographic signal, and expanding the modeled multi-channel surface electromyographic signal based on a preset expansion factor to obtain an expanded surface electromyographic signal; Obtaining the release time series and noise series corresponding to the extended electromyographic signal; Step 2: constructing a cross-correlation matrix of the surface electromyography signal based on the convolution mixture model of the surface electromyography signal, and initializing the activity index through the cross-correlation matrix; Step 3: Initialize the cross-correlation vector based on the activity index and set the gradient function and learning rate, introduce the exponentially weighted moving average model, perform offline calculation and iteratively update the cross-correlation vector; Step 4, constructing a sliding window for the original multi-channel surface electromyographic signal, extracting and expanding the surface electromyographic signal in the sliding window, calculating the expanded cross-correlation matrix, using the expanded cross-correlation matrix to estimate the firing time series of each motor unit, and adjusting the firing time series according to the updated cross-correlation vector, repeating step 4 until the decomposition of all signals is completed; Step 5, comparing the discharge time series with a preset threshold value, if the peak interval variation coefficient or the average discharge rate of the discharge time series does not meet the requirements, then discard it, stop stripping, and jump to step 8, otherwise, continue stripping and proceed to step 6; Step 6: Based on the firing time series, the action potential of each motor unit is calculated, and the action potential is subtracted from the original electromyographic signal to generate a residual signal, and the residual signal is re-input into step 4 for the next round of decomposition; Step 7, repeating steps 5-6, continuing to extract the sliding window signal, estimating the release time series, and stripping the components of the motion unit from the residual signal until no valid motion unit release time series can be identified; Step 8: Complete the decomposition of all signals and finally output the effective motor unit release time series.

2. The method for decomposing surface electromyographic signals based on stripping fast gradient convolution kernel compensation according to claim 1, characterized in that: Step 1 specifically includes: Step 1.1: The electromyographic signals generated by the activities of N motor units acquired by M channels are expressed as: Among them, x i (n) represents the nth sample of the ith channel, h ij is the action potential of the jth motor unit in the ith channel with length L, s j (n) represents the firing time series of the jth motor unit; Adding noise to the electromechanical signal, the electromyographic signal is expressed as: y i (n)=x i (n)+ω i (n) Among them, y i (n) represents the myoelectric signal containing noise, ω i (n) represents the noise in the signal; Step 1.2: Expand the electromyographic signal based on the preset expansion factor K-1, and the obtained electromyographic signal is expressed as: The corresponding emission time series and noise are: Step 1.3: The expanded surface electromyography signal convolution mixture model is expressed as Where H is the mixing matrix: Among them, h ij represents the action potential of the jth motor unit with length L in the ith channel.

3. The method for decomposing surface electromyographic signals based on stripping fast gradient convolution kernel compensation according to claim 2, characterized in that: Step 2 specifically includes: Based on the convolution mixture model of surface electromyography signals, the cross-correlation matrix is ​​constructed: Among them, E() represents the calculation of mathematical expectation; Initialize the activity index:

4. The method for decomposing surface electromyographic signals based on stripping fast gradient convolution kernel compensation according to claim 3, characterized in that: Step 3 specifically includes: Initialize the cross-correlation vector based on the activity index: n0=maxarg n (γ(n)) γ(n0)=0 in, represents the initial cross-correlation vector, n0 represents the index when γ(n) reaches its maximum value, represents the expanded signal at index n0, Set the gradient function and learning rate: set up is the initial gradient function, and the learning rate is selected as η = 0.

01. Introducing the exponentially weighted moving average model: y t =(1-λ)x t +(1-λ)·λx t-1 +(1-λ)·λ 2 ×x t-2 +L+(1-λ)·λ t-1 x1+λ t y0 Among them, λ is a weight variable with a value between 0 and 1, t is the current time, and y t is the weighted average of the gradient at the current time t, x t is the gradient at the current time t; The gradient function recalculated by introducing the exponentially weighted moving average model is: in, is the firing time series of the jth motor unit at the mth sample, λ and β are weight variables ranging from 0 to 1, k is the number of iterations, when k = 1, the gradient is initialized The update rule of the cross-correlation vector in the offline decomposition process is: Where k is the number of iterations, is the k+1th cross-correlation vector, is the k-th cross-correlation vector, η represents the k-th learning rate, ε=10 -8 , Represents the mth sample of the extended EMG signal.

5. The method for decomposing surface electromyographic signals based on stripping fast gradient convolution kernel compensation according to claim 4, characterized in that: In step 4, the cross-correlation vector and cross-correlation matrix are updated according to the following formula: Among them, j ={n1,n2,...} is the release time series decomposed from the cross-correlation vector calculated in the offline process in the sliding window data, l is the learning rate, is the cross-correlation matrix, is the cross-correlation vector, Represents the surface electromyographic signal after delay expansion within the sliding window.

6. The method for decomposing surface electromyographic signals based on stripping fast gradient convolution kernel compensation according to claim 5, characterized in that: In step 5, the peak interval variation coefficient is calculated according to the following formula: Where SD(ISI) is the standard deviation of the spike interval, μ(ISI) is the mean of the spike interval, and ISI is expressed as the time interval between two consecutive discharge events; The average discharge rate is calculated according to the following formula: Where N is the number of discharges and T is the observation time.

7. The method for decomposing surface electromyographic signals based on stripping fast gradient convolution kernel compensation according to claim 6, characterized in that: Step 6 specifically includes: Step 6.

1. Calculate the action potential of each motor unit using the minimum mean square error method: in, is the estimated motor unit action potential matrix, S is the expanded peak column matrix, X is the multi-channel sEMG signal matrix, S T is the transpose of matrix S, S T S is the covariance matrix of the peak column, S T X is the correlation matrix between the signal and the peak column; Step 6.2: Calculate the residual signal using the reconstructed signal: First, the reconstructed signal of the identified motion unit is calculated: Among them, S is the expanded peak column matrix, is the estimated motor unit action potential matrix, Subtract the reconstructed signal from the original signal to obtain the residual signal: Among them, X residual is the residual signal matrix, which is used for the next iteration.

Images