Cerebral stroke patient upper limb rehabilitation evaluation method based on electroencephalogram and electromyographic signals

Abstract

The invention discloses a cerebral apoplexy patient upper limb rehabilitation evaluation method based on electroencephalogram and electromyographic signals. The method specifically comprises the following steps that signals are collected and preprocessed; electroencephalogram and electromyographic signal decomposition; performing multi-scale processing; calculating transmission spectrum entropy; and upper limb rehabilitation grade evaluation: mapping a fitting relationship between the MSTSE significance area value of each patient and a Brunnstrom evaluation scale based on a Spearman correlation coefficient, analyzing an FCMC characteristic index and a change trend which are significantly related to the grade of the Brunnstrom evaluation scale, and performing quantitative evaluation on the disease grade of the stroke patient. According to the technical scheme, a binary variational mode decomposition-multi-scale transfer spectrum entropy method is used for representing the FCMC coupling relation between the brain motor cortex and the corresponding muscular tissue, and in combination with a Brunnstrom rating scale, objective quantitative evaluation of the rehabilitation level of the upper limbs of the stroke patient is achieved.

Description

Technical Field

[0001] The present invention relates to the technical field of bioelectric signal processing, and in particular to an upper limb rehabilitation assessment method for stroke patients based on brain electromyographic signals. Background Art

[0002] In recent years, the risk of stroke patients ranks first in the world. Every year, there are more than 20 million patients with varying degrees of limb dysfunction due to cerebrovascular damage, which has a serious negative impact on people's normal life, work, family and society. Effective targeted rehabilitation training is conducive to the recovery of relevant motor functions of patients, but accurate rehabilitation assessment after stroke is the prerequisite for various rehabilitation training. At present, traditional clinical assessment methods are mostly completed by doctors based on clinical assessment scales and the degree of action execution of patients. This method has the following disadvantages: the assessment results are greatly affected by the doctor's subjective influence, and different doctors may give different judgment results for the same patient based on different assessment scales and their experience; some subtle manifestations of patients are not discovered by doctors in time, resulting in inaccurate assessment results.

[0003] With the development of medical technology, computer and electronic technology, clinical diagnosis and rehabilitation methods after stroke have gradually become more diverse. Electroencephalogram (EEG) and surface electromyography (sEMG) can quickly characterize the changes in the patient's damaged brain and muscle function. They are inexpensive and have been widely used in clinical practice. When the human body performs actions autonomously, the EEG generated by the cerebral cortex first reaches the motor unit through the central nervous system. The motor unit controls the muscle fiber to complete muscle contraction, generates action potentials on the muscle fiber, and realizes the execution of the action. This process reflects the brain's control of limb movement information and the muscle's response information to the brain's control of movement, which is called functional cortical muscle coupling (FCMC). A large number of studies have shown that synchronized EEG and sEMG information can characterize the functional coupling between the brain and its controlled muscles, reflecting the functional state of the motor nervous system. Stroke-induced damage to the cortex and neural pathways leads to local brain neurological dysfunction, which cannot send accurate and complete nerve impulses to the peripheral nervous system, thereby causing motor dysfunction.

[0004] There are many methods for studying FCMC, such as Granger Causality (GC), Transfer Entropy (TE), and generalized partial directed coherence (gPDC). GC can effectively calculate the causal relationship between signal sequences and reflect the directionality between cortical and muscle coupling, but GC can only measure the linear relationship between signals. Both EEG and sEMG are random signals, so the nonlinear relationship between signals cannot be well characterized. The TE algorithm quantitatively analyzes the nonlinear directed coupling relationship between signals and is widely used in FCMC research. However, TE can only capture the information flow relationship in one direction. For the bidirectional or multidirectional information flow in complex systems, TE lacks an effective explanation. The gPDC algorithm can distinguish the directional connection between signals and determine the direction of information transmission. It is also suitable for the causal relationship analysis of multiple time series in complex systems, but the algorithm is sensitive to noise and is greatly affected by parameters.

[0005] Since EMG signals are random signals and contain complex features of different scales, the use of decomposition algorithms is more common in FCMC, such as wavelet packets, empirical mode decomposition (EMD), variational mode decomposition (VMD), etc. The wavelet packet algorithm is often affected by the selection of basis functions. Improper selection will affect the decomposition effect and analysis results. The EMD algorithm has endpoint effects, is inherently uncertain, and lacks mathematical theoretical support. The VMD algorithm has high resolution and strong adaptability, but it is greatly affected by the initial parameters and can only process single-channel signals. Summary of the invention

[0006] The technical solution of the present invention to solve the above technical problem is to provide a method for evaluating upper limb rehabilitation of stroke patients based on brain electromyographic signals, comprising the following steps:

[0007] Signal acquisition and preprocessing;

[0008] Brain and myoelectric signal decomposition:

[0009] Construct a complex signal, combine the EEG signal x(t) and the EMG signal y(t) into a complex form, defined as: z(t) = x(t) + jy(t), where j is an imaginary unit;

[0010] Direction selection, select a set of directions to project z(t), select multiple different directions θ d , where θ dis a uniformly distributed angle on the unit circle. When D directions are selected, each direction θ d It can be expressed as: d=0,1,2,…,D-1;

[0011] Projection calculation, project the complex signal z(t) to a certain direction θ, and calculate the projection p of z(t) in this direction θ (t), p θ (t) = Re (z (t) e -jθ ), where Re represents the real part and θ is the angle of the projection direction;

[0012] Parameter initialization, initialization of the parameters required for VMD, including the number of modes K, the center frequency w of each mode k , Lagrange multiplier λ, modal function u k , penalty parameter α and convergence tolerance tol;

[0013] Define the objective function, requiring each mode u k (t) have finite bandwidth and center frequency, and the constraint is that the sum of all modes is equal to the input signal p θ (t), the sum of the estimated bandwidths of each mode is minimized, and the constrained variational model is as follows:

[0014]

[0015] st∑ k u k =p θ (t);

[0016] In the formula, {w k} represents the frequency center of each intrinsic mode function (IMFs), {u k} are the K finite bandwidth IMFs components obtained after decomposition, To find the partial derivative, * is the convolution operator, δ(t) is the Dirac function, Indicates that the spectrum is shifted to the baseband;

[0017] Variational solution, introduce Lagrange multiplication operator λ and penalty coefficient α, convert the constrained variational problem into an unconstrained problem, and then use alternating direction multiplier (ADMM) to continuously update the center frequency w of each mode k , Lagrange multiplier λ and modal function u k , to obtain the optimal solution of the constraint model;

[0018] Determine the eigenmode function, and through the continuous iteration of the sub-step variational solution, stop the iteration when the constraint conditions are met to obtain each mode u k (t), and then the real data P is obtained by performing inverse Fourier transform on it.θ (t), each P θ (t) contains the real part X(t) and the imaginary part Y(t), then the real part X(t) is the IMFs set of the EEG signal, and the imaginary part Y(t) is the IMFs set of the EMG signal.

[0019] Multi-scale processing, using the moving average method to analyze IMFs at multiple scales, extract the local change characteristics of brain electromyographic signals, and choose to use the sliding window to average the IMFs sequence:

[0020]

[0021] Among them, L(P θ (t),τ) represents the moving average under scale τ, P θ (t) represents each brain and myoelectric IMFs sequence, T is the length of the IMFs sequence, and the brain and myoelectric IMFs after multi-scale processing are defined as X S With Y S ;

[0022] Transfer spectral entropy calculation: After completing multi-scale processing, the EMG IMFs sequence is reconstructed in phase space, and then converted into a two-dimensional matrix through two-dimensional Fourier transform to calculate the transfer spectral entropy.

[0023] The upper limb rehabilitation grade was assessed by mapping the fitting relationship between the MSTSE significance area value of each patient and the Brunnstrom rating scale based on the Spearman correlation coefficient, analyzing the FCMC characteristic indicators and change trends that were significantly correlated with the Brunnstrom rating scale grade, and quantitatively assessing the disease grade of stroke patients.

[0024] Furthermore, the signal acquisition and preprocessing steps include:

[0025] Combined with the Brunnstrom rating scale, according to the clinical manifestations of different stroke patients, an arm flexion experimental paradigm was designed under the guidance of clinicians. According to the international 10-20 system, EEG of stroke patients at leads Fp1, Fp2, F3, F4, C3, C4, Cz, P3, P4, A1, and A2 were collected respectively. Surface electromyography sensors were used to collect sEMG at the anterior bundle of the deltoid muscle, middle bundle of the deltoid muscle, biceps brachii, triceps brachii, and extensor digitorum of the forearm. The EEG was preprocessed, and the specific operations included: re-referencing, baseline correction, filtering, removal of oculoculography, downsampling, and EEG data segmentation; the sEMG was preprocessed, including: baseline correction, filtering, EMG data segmentation, etc.

[0026] Furthermore, the step of calculating the transfer spectral entropy includes:

[0027] Perform phase reconstruction on the multi-scale processed sequence and convert X S With Y S Expanded to a two-dimensional matrix P S With Q S ;

[0028] The time domain matrix P S With Q S Through the two-dimensional Fourier transform to the frequency domain, we get the matrix W S and V S ;

[0029] By calculating the M-dimensional vector w S (f) and v S (f) The transfer spectral entropy reveals the information transfer between different frequency components.

[0030] Furthermore, the variational solution step includes:

[0031] Introduce the Lagrange multiplication operator λ and penalty coefficient α:

[0032]

[0033] In the formula, α is the penalty function to ensure the reconstruction accuracy of the signal, and λ is the lagrange multiplier to ensure the strictness of the constraint conditions;

[0034] Based on Plancherel Fourier transform, the above formula is converted to the frequency domain for calculation, and the center frequency w of each mode is continuously updated. k , Lagrange multiplier λ and modal function u k , to obtain the optimal solution of the constraint model:

[0035]

[0036] Where ∧ represents Fourier transform, n is the number of iterations, and η is the fidelity coefficient.

[0037] In combination with the above technical solutions and the technical problems solved, the advantages and positive effects of the technical solutions to be protected by the present invention are as follows:

[0038] First, the binary variational mode decomposition-multiscale transfer spectral entropy method is used to characterize the FCMC coupling relationship between the motor cortex of the brain and the corresponding muscle tissue, and combined with the Brunnstrom rating scale to achieve an objective quantitative assessment of the upper limb rehabilitation level of stroke patients.

[0039] (1) The technical solution of the present invention solves a technical problem that people have been eager to solve but have never been able to solve successfully:

[0040] Based on bivariate variational mode decomposition (mBVMD), this paper synchronously decomposes brain electromyographic signals, obtains brain electromyographic signals of two same frequency bands, realizes the same-frequency coupling analysis of brain electromyographic signals, and explores the neuromuscular coupling characteristics. Compared with the VMD algorithm, BVMD expands the one-dimensional signal to two dimensions, ensures the matching of intrinsic mode function components in quantity and scale, and solves the problem of modal calibration between signals. At the same time, based on the phase space reconstruction technology, the frequency domain information of brain electromyographic signals is extracted, and the transfer spectral entropy of patients with different levels of stroke is calculated, thereby realizing the quantitative and objective assessment of the disease level of stroke patients.

[0041] (2) The technical solution of the present invention overcomes technical prejudice:

[0042] In the clinical assessment process of stroke disease grade, clinicians complete the assessment of stroke patients' disease grade based on the grade rating scale and clinical subjective experience, which is somewhat subjective. The present invention collects synchronous brain and myoelectric signals of stroke patients when they perform rehabilitation movements, mines the coupling degree between signals, maps the relationship between different assessment grades and signal coupling, and realizes the assessment of stroke patients' disease grade. To a certain extent, it overcomes the current technical bias of subjective clinical assessment of stroke patients' disease grade. BRIEF DESCRIPTION OF THE DRAWINGS

[0043] In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the drawings required for use in the embodiments or the description of the prior art will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on the structures shown in these drawings without paying creative work.

[0044] Figure 1 This is a flow chart of the upper limb rehabilitation assessment method for stroke patients based on brain myoelectric signals according to the present invention;

[0045] Figure 2 Schematic diagram of the experimental paradigm;

[0046] Figure 3 Execute a flow chart for the experimental paradigm;

[0047] Figure 4 This is a schematic diagram of the EEG leads used in the present invention;

[0048] Figure 5 A schematic diagram of the myoelectric tissue used in the present invention;

[0049] Figure 6 Decomposition of IMFs waveform for complex EMG signals;

[0050] Figure 7 It is the MSTSE heat map of brain electromyographic signal;

[0051] Figure 8 This is a mapping diagram between BVMD-MSTSE and Brunnstrom rating scale in stroke patients. DETAILED DESCRIPTION

[0052] The present invention proposes an upper limb rehabilitation assessment method for stroke patients based on brain electromyography signals, aiming to design a rehabilitation assessment method for objectively and quantitatively assessing the upper limb rehabilitation level of stroke patients.

[0053] The following is a description of the upper limb rehabilitation assessment method for stroke patients based on brain myoelectric signals proposed by the present invention in a specific embodiment:

[0054] In the technical solution of this embodiment, a method for evaluating upper limb rehabilitation of stroke patients based on brain electromyographic signals comprises the following steps:

[0055] Step 1: Signal acquisition and preprocessing;

[0056] Step 2: Decomposition of brain and myoelectric signals:

[0057] Construct a complex signal, combine the EEG signal x(t) and the EMG signal y(t) into a complex form, defined as: z(t) = x(t) + jy(t), where j is an imaginary unit;

[0058] Direction selection, select a set of directions to project z(t), select multiple different directions θ d , where θ d is a uniformly distributed angle on the unit circle. When D directions are selected, each direction θ d It can be expressed as: d=0,1,2,…,D-1;

[0059] Projection calculation, project the complex signal z(t) to a certain direction θ, and calculate the projection p of z(t) in this direction θ (t), p θ (t) = Re (z (t) e -jθ ), where Re represents the real part and θ is the angle of the projection direction;

[0060] Parameter initialization, initialization of the parameters required for VMD, including the number of modes K, the center frequency w of each mode k , Lagrange multiplier λ, modal function u k , penalty parameter α and convergence tolerance tol;

[0061] Define the objective function, requiring each mode u k(t) has a finite bandwidth and center frequency, with the constraint that the sum of the modes is equal to the input signal p θ (t), the sum of the estimated bandwidths of each mode is minimized, and the constrained variational model is as follows:

[0062]

[0063] st∑ k u k =p θ (t);

[0064] In the formula, {w k} represents the frequency center of each IMFs component, {u k} are the K finite bandwidth IMFs components obtained after decomposition, To find the partial derivative, * is the convolution operator, δ(t) is the Dirac function, Indicates that the spectrum is shifted to the baseband;

[0065] Variational solution, introduce Lagrange multiplication operator λ and penalty coefficient α, convert the constrained variational problem into an unconstrained problem, and then use alternating direction multiplier (ADMM) to continuously update the center frequency w of each mode k , Lagrange multiplier λ and modal function u k , to obtain the optimal solution of the constraint model;

[0066] Determine the eigenmode function, and through the continuous iteration of the sub-step variational solution, stop the iteration when the constraint conditions are met to obtain each mode u k (t), and then the real data P is obtained by performing inverse Fourier transform on it. θ (t), each P θ (t) contains the real part X(t) and the imaginary part Y(t), then the real part X(t) is the IMFs set of the EEG signal, and the imaginary part Y(t) is the IMFs set of the EMG signal.

[0067] Step 3: Multi-scale processing, use the moving average method to analyze the IMF multi-scale, extract the local change characteristics of the brain and myoelectric signals, and choose to use the sliding window to average the IMFs sequence:

[0068]

[0069] Among them, L(P θ (t),τ) represents the moving average under scale τ, P θ (t) represents each brain and myoelectric IMFs sequence, T is the length of the IMFs sequence, and the brain and myoelectric IMFs after multi-scale processing are defined as X S With Y S ;

[0070] Step 4: Transfer spectral entropy calculation: perform phase space reconstruction on the EMG IMFs sequence after multi-scale processing, and then convert the reconstructed sequence into a two-dimensional matrix through two-dimensional Fourier transform to calculate the transfer spectral entropy;

[0071] Step 5: Assessment of upper limb rehabilitation level. Based on the Spearman correlation coefficient, the fitting relationship between the MSTSE significance area value of each patient and the Brunnstrom rating scale was mapped. The FCMC characteristic indicators and change trends that were significantly correlated with the Brunnstrom rating scale level were analyzed to quantitatively assess the disease level of stroke patients.

[0072] Furthermore, the signal acquisition and preprocessing steps include:

[0073] Combined with the Brunnstrom rating scale, according to the clinical manifestations of different stroke patients, an arm flexion experimental paradigm was designed under the guidance of clinicians. According to the international 10-20 system, EEG of stroke patients at leads Fp1, Fp2, F3, F4, C3, C4, Cz, P3, P4, A1, and A2 were collected respectively. Surface electromyography sensors were used to collect sEMG at the anterior bundle of the deltoid muscle, middle bundle of the deltoid muscle, biceps brachii, triceps brachii, and extensor digitorum of the forearm. The EEG was preprocessed, and the specific operations included: re-referencing, baseline correction, filtering, removal of oculoculography, downsampling, and EEG data segmentation; the sEMG was preprocessed, including: baseline correction, filtering, EMG data segmentation, etc.

[0074] Furthermore, the step of calculating the transfer spectral entropy includes:

[0075] Perform phase reconstruction on the multi-scale processed sequence and convert X S With Y S Expanded to a two-dimensional matrix P S With Q S ;

[0076] The time domain matrix P S With Q S Through the two-dimensional Fourier transform to the frequency domain, we get the matrix W S and V S ;

[0077] By calculating the M-dimensional vector w S (f) and v S (f) The transfer spectral entropy reveals the information transfer between different frequency components.

[0078] Furthermore, the variational solution step includes:

[0079] Introduce the Lagrange multiplication operator λ and penalty coefficient α:

[0080]

[0081] In the formula, α is the penalty function to ensure the reconstruction accuracy of the signal, and λ is the lagrange multiplier to ensure the strictness of the constraint conditions;

[0082] Based on Plancherel Fourier transform, the above formula is converted to the frequency domain for calculation, and the center frequency w of each mode is continuously updated. k , Lagrange multiplier λ and modal function u k , to obtain the optimal solution of the constraint model:

[0083]

[0084] Where ∧ represents Fourier transform, n is the number of iterations, and η is the fidelity coefficient.

[0085] Embodiment 1:

[0086] A method for evaluating upper limb rehabilitation of stroke patients based on brain electromyographic signals. The specific process is as follows: Figure 1 As shown, the following steps are included:

[0087] Step 1: Signal acquisition and preprocessing. Under the guidance of clinicians, arm flexion was set as the experimental paradigm in combination with the Brunnstrom rating scale, and the synchronous EEG signals and surface electromyographic signals of several stroke patients were collected when they performed arm flexion, ensuring that at least 5 subjects of each level of Brunnstrom were included in the stroke subjects.

[0088] The experimental paradigm includes two actions: arm flexion and rest, such as Figure 2 As shown in the figure. When the subject is performing arm flexion, a prompter is placed in front of the subject to give the subject a picture and sound prompt to prompt the subject to perform the corresponding action. During the signal collection process, the subject needs to sit upright on the chair, with the arm hanging naturally, keep the body stable, and reduce unnecessary movements, including but not limited to: turning the head, swallowing, blinking, etc.; Figure 3The figure shows the specific action execution process. The preparation stage is from 0 to 3 seconds. The subject is given a 1-second voice prompt and a "+" is displayed on the prompter to focus the subject's attention. During this period, the subject's arm naturally hangs down. In the action execution stage, the subject needs to slowly bend his arm forward until his arm is parallel to the ground and remains stable. The action execution stage lasts for 8 seconds. In the rest stage, the subject needs to slowly lower his arm until his arm naturally hangs down and remains vertical to the ground. The rest stage lasts for 12 seconds. Each subject needs to complete 3 sessions of actions. Each session includes 10 arm flexion actions and 10 rest actions. After completing 1 session of actions, the subject will get a 10-minute rest time to prevent mental and muscle fatigue.

[0089] In the EEG acquisition process, the international 10-20 system is strictly followed. The subjects need to clean their scalp in advance and wear the EEG cap correctly. The reference point is the central lead Cz. To avoid distortion of the collected EEG signal, the impedance between the electrode and the scalp must be less than 5kΩ before the EEG signal is collected. Figure 4 As shown, EEG is collected from leads Fp1, Fp2, F3, F4, C3, C4, Cz, P3, P4, A1, and A2 respectively.

[0090] During the sEMG acquisition process, alcohol was used to clean the skin of the subject's upper arm, and scrub was used to polish the subject's local skin when necessary to reduce skin impedance. Figure 5 As shown, sEMGs were collected from the anterior deltoid, middle deltoid, biceps, triceps and extensor digitorum of the forearm.

[0091] The EEG preprocessing operations include but are not limited to re-referencing, baseline correction, filtering, electrooculogram removal, downsampling, and EEG data segmentation. In the re-referencing operation, A1 and A2 are selected as new reference points; in the baseline correction operation, the baseline drift of EEG is reduced by the de-averaging method; in the filtering operation, 0.1-200Hz bandpass filtering and 50Hz notch processing are performed on the EEG; in the electrooculogram removal operation, the ICA algorithm is used to remove the subject's electrooculogram artifacts; in the EEG signal downsampling operation, the sampling rate of the EEG signal needs to be the same as the sampling rate of the electromyography signal; in the EEG data segmentation operation, the time range of the action execution stage is determined based on the marks in the EEG acquisition system, and the corresponding data segment is extracted.

[0092] The preprocessing operation of the sEMG includes but is not limited to baseline correction, filtering, electromyographic data segmentation, etc. In the baseline correction operation, the baseline drift of the sEMG is reduced by taking the average method; in the filtering operation, the sEMG is subjected to 0.1-200Hz bandpass filtering and 50Hz notch processing; in the electromyographic data segmentation operation, the time range of the action execution stage is determined based on the marks in the electromyographic acquisition system, and the corresponding data segment is extracted.

[0093] Step 2: Decomposition of brain and myoelectric signals. This embodiment uses binary variational mode decomposition (BVMD) to achieve synchronous decomposition of brain and myoelectric signals, which specifically includes the following contents:

[0094] 2.1 Constructing complex signals: The preprocessed EEG signal x(t) and surface electromyography signal y(t) are combined into a complex signal z(t)=x(t)+jy(t), where j is an imaginary unit;

[0095] 2.2 Direction selection: Select a set of directions to project z(y). Generally, you can select multiple different directions θ d , where θ d is a uniformly distributed angle on the unit circle. When D directions are selected, each direction θ d It can be expressed as: d=0,1,2,…,D-1, the value range of D is set between 10-20;

[0096] 2.3 Projection calculation: Project the complex signal z(t) to a certain direction θ and calculate the projection signal p of z(t) in this direction θ (t), p θ (t) = Re (z (t) e -jθ ), where Re represents the real part and θ is the angle of the projection direction;

[0097] 2.4 Parameter initialization: Initialize the parameters required for VMD, including the number of modes K and the center frequency w of each mode k , Lagrange multiplier λ, modal function u k , penalty parameter α and convergence tolerance tol;

[0098] 2.5 Define the objective function: require each mode u k (t) has a finite bandwidth and center frequency, with the constraint that the sum of the modes is equal to the input signal p θ (t), the sum of the estimated bandwidths of each mode is minimized, and the constrained variational model is as follows:

[0099]

[0100] st∑k u k =p θ (t);

[0101] In the formula, {w k} represents the frequency center of each IMFs component, {u k} are the K finite bandwidth IMFs components obtained after decomposition, To find the partial derivative, * is the convolution operator, δ(t) is the Dirac function, Indicates that the spectrum is shifted to the baseband;

[0102] 2.6 Variational solution. To solve the above constrained optimal problem, the Lagrange multiplication operator λ and penalty coefficient α are introduced to convert the constrained variational problem into an unconstrained problem:

[0103]

[0104] In the formula, α is the penalty function to ensure the reconstruction accuracy of the signal, and λ is the lagrange multiplier to ensure the strictness of the constraint conditions;

[0105] The variational problem is further solved by the alternating direction multiplier (ADMM). Based on the Plancherel Fourier transform, the above equation is converted to the frequency domain for calculation, and the center frequency w of each mode is continuously updated. k , Lagrange multiplier λ and modal function u k , to obtain the optimal solution of the constraint model:

[0106]

[0107] Where ∧ represents Fourier transform, n is the number of iterations, and η is the fidelity coefficient.

[0108] 2.7 Determine the eigenmode function. Set the constraints:

[0109]

[0110] Through the continuous iteration of sub-step 2.6, when the constraint conditions are met, the iteration is stopped to obtain each mode u k (t), and then the real data P is obtained by performing inverse Fourier transform on it. θ (t), each P θ (t) contains the real part X(t) and the imaginary part Y(t), then the real part X(t) is the IMFs set of the EEG signal, and the imaginary part Y(t) is the IMFs set of the EMG signal.

[0111] Step 3: Multi-scale processing. Use the moving average method to perform multi-scale analysis on IMF and extract the local change characteristics of brain electromyographic signals:

[0112]

[0113] Among them, L(P θ (t),τ) represents the moving average under scale τ, P θ (t) represents each brain and myoelectric IMFs sequence, T is the length of the IMFs sequence, and the brain and myoelectric IMFs after multi-scale processing are defined as X S With Y S .

[0114]

[0115] Step 4: Calculation of transfer spectral entropy: After completing multi-scale processing, the EMG IMFs sequence is reconstructed in phase space, and then converted into a two-dimensional matrix through two-dimensional Fourier transform to calculate the transfer spectral entropy.

[0116] 4.1 In order to better capture the complex dynamic characteristics and frequency domain characteristics of the signal, the multi-scale processed sequence is firstly reconstructed. S With Y S Expanded to a two-dimensional matrix P S With Q S :

[0117]

[0118] In the formula, σ = 1, 2, …, N-(M-1)Δt, where M is the embedding dimension, which is usually set to 15, and Δt is the embedding delay, which is generally set to 2.

[0119] 4.2 The time domain matrix P S With Q S Through the two-dimensional Fourier transform to the frequency domain, we get the matrix W S and V S :

[0120]

[0121] Wherein, the frequency resolution is Δf=fs / N, fs is the sampling frequency of the brain electromyographic signal, and N=T-τ+1.

[0122] 4.3 By calculating the M-dimensional vector w S (f) and v S The transfer spectrum entropy of (f) reveals the information transfer between different frequency components:

[0123]

[0124] Step 5: Assessment of upper limb rehabilitation level.

[0125] In the upper limb rehabilitation grade assessment, the consistency relationship between the transfer spectral entropy and the clinical assessment scale was explored based on the Spearman method. The calculation process of the Spearman correlation coefficient is as follows:

[0126]

[0127] Where ρ is the Spearman correlation coefficient, h i is the rank difference of each group of data. In the present invention, the result of each group of data is recorded as The Brunnstrom rank correlation, num is the total number of observations, and the mean value ρ is obtained after calculating the MSTSE significant area value between each muscle tissue and the brain mean , to characterize the fitting relationship between the MSTSE significance area value and the Brunnstrom rating scale of each patient, analyze the FCMC characteristic indicators and change trends that are significantly correlated with the Brunnstrom rating scale level, and complete the quantitative assessment of the disease level of stroke patients.

[0128] Step 6: Example Testing

[0129] The present invention collected brain electromyographic signals of 6 stroke patients when they performed the experimental paradigm. All patients had right upper limb hemiplegia and no cognitive impairment, no history of any mental illness, all patients were informed of the experimental protocol before the experiment, and had never had similar experimental experience. Among all stroke patients, subject 1 corresponds to stage 6 patients in the Brunnstrom upper limb assessment scale, subjects 2 and 5 correspond to stage 5 patients in the Brunnstrom upper limb assessment scale, subjects 3 and 6 correspond to stage 4 patients in the Brunnstrom upper limb assessment scale, and subject 4 corresponds to stage 3 patients in the Brunnstrom upper limb assessment scale.

[0130] The pre-processed action segment EEG and myoelectric signals are combined into complex signals, and the complex signals are decomposed based on BVMD. Since the most significant brain area that controls human limb movement is located in the contralateral motor area of ​​the brain, the present invention selects the EEG at the C3 lead of the motor area and each EMG signal for coupling characteristic analysis. Figure 6 The figure shows the IMFs waveform diagram after decomposing a complex signal, and 5 groups of IMFs signal waveform diagrams are obtained.

[0131] After completing the multi-scale processing of the IMFs signal using the moving average method, the coupling strength between the multi-scale signals is calculated based on the transfer spectral entropy, such as Figure 7The figure shows the distribution diagram of the average value of the spectral entropy transferred between the brain and myoelectric signals of the six subjects, where the ordinate represents each EEG IMFs, and the abscissa represents the IMFs of the anterior deltoid muscle, middle deltoid muscle, biceps brachii, triceps brachii and extensor digitorum of the forearm respectively. As can be seen from the figure, the MSTSE coupling intensity of Subject 1 is the strongest, and the MSTSE coupling intensity of Subject 4 is the weakest. Since the Brunnstrom stage of Subject 1 is stage 6, the muscle contraction state of Subject 1 when performing the paradigm action is similar to that of normal people. The coupling intensity between the anterior bundle of the deltoid muscle and the middle bundle of the deltoid muscle and the brain is the strongest, the coupling intensity between the biceps and the triceps and the brain is second, and the coupling intensity between the forearm finger extensor and the brain is the lowest; the Brunnstrom stage of Subject 4 is stage 3. When performing the paradigm action, the phenomenon of shoulder adduction and elbow flexion occurs, which is the biggest difference from the muscle contraction state of normal people. The coupling intensity between the biceps and the triceps and the brain is the strongest, the coupling intensity between the anterior bundle of the deltoid muscle and the middle bundle of the deltoid muscle and the brain is second, and the coupling intensity between the forearm finger extensor and the brain is the lowest.

[0132] By calculating the significant area value of MSTSE in patients with various grades of stroke, the correlation between each MSTSE and the Brunnstrom rating scale was mapped based on Spearman correlation analysis, such as Figure 8 The figure shows the mapping relationship diagram. It can be seen from the figure that there is a significant correlation between the MSTSE significant area value between the brain and myoelectric signals and the Brunnstrom rating scale, and there is a positive correlation, indicating that the FCMC characteristics calculated based on BVMD-MSTSE can be used to evaluate the rehabilitation status of stroke patients.

[0133] The above is only a preferred specific embodiment of the present invention, but the protection scope of the present invention is not limited thereto. Any changes or substitutions that can be easily thought of by a person skilled in the art within the technical scope disclosed by the present invention should be included in the protection scope of the present invention. Therefore, the protection scope of the present invention should be based on the protection scope of the claims.

Claims

1. A method for evaluating upper limb rehabilitation of stroke patients based on brain electromyographic signals, characterized in that: The following steps are involved: Signal acquisition and preprocessing; Brain and myoelectric signal decomposition: Construct a complex signal, combine the EEG signal x(t) and the EMG signal y(t) into a complex form, defined as: z(t) = x(t) + jy(t), where j is an imaginary unit; Direction selection, select a set of directions to project the complex signal z(t), select multiple different directions θ d , where θ d is a uniformly distributed angle on the unit circle, select D directions, each direction θ d It can be expressed as: d=0,1,2,…,D-1; Projection calculation, project the complex signal z(t) to a certain direction θ, and calculate the projection p of z(t) in this direction θ (t), p θ (t) = Re (z (t) e -jθ ), where Re represents the real part and θ is the angle of the projection direction; Parameter initialization, initialization of the parameters required for VMD, including the number of modes K, the center frequency w of each mode k , Lagrange multiplier λ, modal function u k , penalty parameter α and convergence tolerance tol; Define the objective function, requiring each mode u k (t) has a finite bandwidth and center frequency, with the constraint that the sum of the modes is equal to the input signal p θ (t), the sum of the estimated bandwidths of each mode is minimized, and the constrained variational model is as follows: In the formula, {w k } represents the frequency center of each intrinsic mode function (IMFs), {u k } are the K finite bandwidth IMFs components obtained after decomposition, To find the partial derivative, * is the convolution operator, δ(t) is the Dirac function, Indicates that the spectrum is shifted to the baseband; Variational solution, introduce Lagrange multiplication operator λ and penalty coefficient α, alternating direction multipliers, and continuously update the center frequency w of each mode k , Lagrange multiplier λ and modal function u k , used to obtain the optimal solution of the constraint model; Determine the intrinsic mode function, set the constraints through continuous iteration of step variational solution: When the constraints are met, stop the iteration and get each mode u k (t); perform inverse Fourier transform to obtain the real part data P θ (t), each P θ (t) contains the real part X(t) and the imaginary part Y(t), where the real part X(t) is the IMFs set of the EEG signal and the imaginary part Y(t) is the IMFs set of the EMG signal; Multi-scale processing: Use the moving average method to analyze IMFs at multiple scales, extract the local change characteristics of brain and myoelectric signals, and choose to use a sliding window to average the IMFs sequence: Among them, L(P θ (t),τ) represents the moving average under scale τ, P θ (t) represents each brain and myoelectric IMFs sequence, T is the length of the IMFs sequence, and the brain and myoelectric IMFs after multi-scale processing are defined as X S With Y S ; Transfer spectral entropy calculation: After completing multi-scale processing, the EMG IMFs sequence is reconstructed in phase space, and then converted into a two-dimensional matrix through two-dimensional Fourier transform to calculate the transfer spectral entropy; Assessment of upper limb rehabilitation level: Based on the Spearman correlation coefficient, the fitting relationship between the MSTSE significance area value of each patient and the Brunnstrom rating scale was mapped, and the FCMC characteristic indicators and change trends that were significantly correlated with the Brunnstrom rating scale level were analyzed to quantitatively assess the disease level of stroke patients.

2. The method for evaluating upper limb rehabilitation of stroke patients based on brain myoelectric signals according to claim 1, characterized in that: The signal acquisition and preprocessing steps include: Combined with the Brunnstrom rating scale, according to the clinical manifestations of different stroke patients, an arm flexion experimental paradigm was designed under the guidance of clinicians. According to the international 10-20 system, EEG of stroke patients was collected at leads Fp1, Fp2, F3, F4, C3, C4, Cz, P3, P4, A1, and A2, and surface electromyography sensors were used to collect sEMG at the anterior bundle of the deltoid muscle, middle bundle of the deltoid muscle, biceps brachii, triceps brachii, and extensor digitorum of the forearm. The EEG was preprocessed, including re-referencing, baseline correction, filtering, removal of oculoculography, downsampling, and EEG data segmentation; the sEMG was preprocessed, including baseline correction, filtering, and EMG data segmentation.

3. The method for evaluating upper limb rehabilitation of stroke patients based on brain myoelectric signals according to claim 1, characterized in that: The step of calculating the transfer spectrum entropy comprises: Perform phase reconstruction on the multi-scale processed sequence and convert X S With Y S Expanded to a two-dimensional matrix P S With Q S ; The time domain matrix P S With Q S Through the two-dimensional Fourier transform to the frequency domain, we get the matrix W S and V S ; By calculating the M-dimensional vector w S (f) and v S (f) The transfer spectral entropy is used to reveal the information transfer between different frequency components.

4. The method for evaluating upper limb rehabilitation of stroke patients based on brain myoelectric signals according to claim 1, characterized in that: The variational solution step includes: Introduce the Lagrange multiplication operator λ and penalty coefficient α: In the formula, α is the penalty function to ensure the reconstruction accuracy of the signal, and λ is the lagrange multiplier; Based on Plancherel Fourier transform, the above formula is converted to the frequency domain for calculation, and the center frequency w of each mode is continuously updated. k , Lagrange multiplier λ and modal function u k , and obtain the optimal solution of the constraint model: Where ∧ represents Fourier transform, n is the number of iterations, and η is the fidelity coefficient.

Images