Motor imagery EEG signal processing method and system based on adaptive waveform features

Abstract

The present invention discloses a method and system for processing motor imagery EEG signals based on adaptive waveform features, relating to the fields of neuroscience and information technology. The method comprises the following steps: collecting EEG signals from subjects undergoing motor imagery training, preprocessing the collected raw multi-channel EEG signals, and extracting time series related to motor imagery training; determining the center frequency using an adaptive decomposition method to decompose neural oscillations in multiple frequency bands; calculating the nonlinearity, sharpness, and average power characteristics of the multi-frequency motor imagery signal components; and fusing the nonlinear waveform features and power features to identify effective motor imagery events. The present invention can identify whether a target user is performing effective motor imagery by acquiring the nonlinear waveform features of the EEG signal to be decoded, thereby improving the accuracy of subsequent EEG signal decoding.

Description

Technical Field

[0001] The present invention relates to the technical field of neuroscience and information technology, and in particular to a method and system for processing motor imagery electroencephalogram (EEG) signals based on adaptive waveform features. Background Art

[0002] Motor imagery is widely used in the pathological research and treatment of various movement disorders, including Parkinson's disease and stroke, as well as cognitive memory disorders such as attention deficit hyperactivity disorder and post-traumatic stress disorder. During motor imagery, the subject does not actually perform a movement, but instead repeatedly imagines themselves performing the movement in their mind, leveraging the brain's cognitive functions to improve actual motor performance. During this process, the subject's EEG activity is recorded in real time, and feedback information is conveyed to the subject through sound, images, and other forms, enabling them to self-perceive and regulate brain activity, thereby enhancing the effectiveness of motor imagery. This training method has been shown to help healthy subjects learn to self-regulate beta wave oscillations in the sensorimotor cortex. Therefore, it is worthwhile to further investigate the changes in neural oscillations during motor imagery-neurofeedback training.

[0003] Multiple studies have demonstrated that multi-band neural oscillations undergo corresponding changes during motor imagery. For example, amputees exhibited significantly reduced power in central mu waves (7-10 Hz) and low-frequency beta waves (13-20 Hz) during motor imagery. The power of alpha waves in the sensory cortex has been found to be positively correlated with learning efficiency. However, the extent to which motor imagery influences brain information transmission and activity patterns remains largely unexplored. Specifically, many studies use the average power of each frequency band as a key metric. This represents the average signal strength over time, but it neglects important temporal information and fails to accurately capture the instantaneous fluctuations of neural oscillations. In contrast, neural oscillations associated with motor imagery exhibit rapid changes, which are more pronounced in the signal's waveform fluctuations. Therefore, analyzing the nonlinear waveform characteristics of EEG signals provides an important opportunity to deepen our understanding of the neurophysiological mechanisms of motor imagery in the brain.

[0004] Neural oscillations are important vehicles for information exchange within the same brain region and between different brain regions. Numerous studies have shown that neural oscillations are not regular sinusoidal waveforms, and non-sinusoidal neural oscillation waveforms have been observed in a variety of organisms, including humans, mice, and rabbits. For example, theta oscillations in the hippocampus of mice and rabbits are sawtooth-shaped. Arched alpha oscillations have been observed in the sensory cortex of rats, and gamma cortical pyramidal interneurons can generate sawtooth or arched gamma oscillations. Furthermore, waveform characteristics can serve as indicators of behavioral and disease changes. Belluscio et al. found that hippocampal theta wave oscillations exhibited greater asymmetry during maze exploration and rapid eye movement (REM) sleep compared to resting mice. Ouedraogo et al. compared slow-wave oscillations in hippocampal dentate granule cells of epileptic and healthy rats. The study showed that while the frequency of slow-wave oscillations was essentially the same in the experimental and control groups, the waveforms in epileptic rats were significantly sharper than in the control group. Jackson et al. used multiple waveform features to analyze the beta frequency band of the motor cortex in Parkinson's patients and found that the ratio of the sharpness of the peak to the sharpness of the trough of the oscillation wave in patients with Parkinson's disease was increased compared to healthy controls. Furthermore, waveform asymmetry increased once medication was discontinued, and returned to symmetry with medication. These findings suggest that waveform features may serve as new indicators for monitoring the status of Parkinson's patients. However, Fourier transform-based EEG signal analysis methods inevitably smooth the signal, which can alter or even remove some local features in the oscillation waveform.

[0005] It can be seen that in the existing technology, the changes in the neural oscillation waveform cannot be retained when extracting EEG signal features, and the nonlinear characteristics of the neural oscillation waveform cannot be quantified, which leads to a decrease in the accuracy of waveform analysis and the inability to effectively identify motor imagery signals. Summary of the Invention

[0006] In view of this, the present invention provides a motor imagery EEG signal processing method and system based on adaptive waveform features, which can improve the accuracy and reliability of waveform analysis while retaining the changes in neural oscillation waveforms and successfully quantify the nonlinear characteristics of neural oscillation waveforms.

[0007] To achieve the above-mentioned object, the present invention provides a method for processing motor imagery EEG signals based on adaptive waveform features, the technical solution of which includes the following steps:

[0008] Step 1: Acquire original multi-channel EEG signals, which contain parts marked as being related to motor imagery.

[0009] Preprocess the original multi-channel EEG signals.

[0010] Data segments marked as related to motor imagery training are extracted from the preprocessed signals, and the data segments are spliced ​​by channels to obtain multidimensional EEG signals.

[0011] Step 2: For each channel signal in the multidimensional EEG signal, an adaptive decomposition method is used to decompose the IMF components and divide them into corresponding frequency bands according to the frequency characteristics of each IMF component.

[0012] Step 3: Calculate the nonlinearity, sharpness, and average power characteristics of the IMF components corresponding to each frequency band to obtain the average power, sharpness, and nonlinearity of each channel corresponding to each frequency band.

[0013] Step 4: Integrate the nonlinearity, sharpness, and average power features and use a machine learning classification model to identify effective motor imagery events.

[0014] Furthermore, step 1: pre-processing the original multi-channel EEG signal, the specific process is as follows:

[0015] Preprocess each channel data in the original multi-channel EEG signal one by one; the specific preprocessing process is as follows:

[0016] The average value of all electrodes used is used as the reference to eliminate errors introduced by variations in the original reference electrode.

[0017] Interpolation is performed on sample points whose absolute amplitude exceeds the set threshold.

[0018] The interpolated signal is band-pass filtered using a notch filter.

[0019] Furthermore, data segments marked as related to motor imagery training are extracted from the preprocessed signals, and the data segments are spliced ​​by channel to obtain multidimensional EEG signals, specifically:

[0020] The pre-processed signal is segmented, and the m-second data segments marked as related to motor imagery training are extracted and spliced ​​by channel to obtain the time series set x(t)={x i (t), i=1,2,…,I}, where I is the number of channels.

[0021] Furthermore, the adaptive decomposition method is an integrated empirical mode decomposition method, a multivariate empirical mode decomposition method, or a noise-assisted multivariate empirical mode decomposition method.

[0022] Furthermore, step 2 is specifically divided into the following steps:

[0023] S201: length is L, Represents K corresponding angles on the (I-1)-dimensional unit ball space The direction vector set of ; Represents the 1st to Ith direction vectors; Indicates the 1st to I-1th corresponding angles.

[0024] S202: Using the Hammersley sequence sampling method, a uniform sampling point set is obtained on the (I-1)-dimensional sphere to obtain an I-dimensional spatial direction vector.

[0025] S203: Initialization: Margin

[0026] S204: Extract the hth multivariate intrinsic mode function IMF h (t), the specific steps are:

[0027] a) Set y0(t) = r h-1 (t),p=1.

[0028] b) y p-1 (t) is the input signal, in each direction vector The mapping on Where · represents the vector inner product, * represents the vector scalar product, Represents a vector Modulus.

[0029] c) Confirm The instantaneous moment corresponding to the extreme point c is the number of extreme points.

[0030] d) Yes The extreme point series in Perform multivariate spline interpolation to obtain K corresponding multivariate envelopes

[0031] e) Calculate the mean of the extreme value envelope based on the center of gravity of the extreme value point

[0032] f) Calculate d(t) = y p-1 (t)-m(t), if the mean of d(t) is 0, the h-th multivariate intrinsic mode function IMF is obtained h (t)=d(t), otherwise, y p (t)=d(t), p increases by 1, and b) to f) are repeated.

[0033] S205: From signal r h-1 (t) minus IMF h (t) Get the new margin r h (t) = r h-1 (t)-IMF h(t); If the remainder is less than a given threshold, the algorithm ends and all multivariate intrinsic mode functions and residual components are obtained. Otherwise, h is incremented by 1 and S204) to S205) are repeated.

[0034] Finally, the original signal x(t) becomes the multivariate intrinsic mode function IMF h (t) and the remainder r(t) are of the form H is the number of eigenmode functions.

[0035] S206: Perform frequency division on the above IMF components, find the center frequency of each component, and classify it into six frequency bands: delta band, theta band, alpha band, low beta band, high beta band, and gamma band. The corresponding range of the delta band is [0.5 Hz, 4 Hz], the corresponding range of the theta band is (4 Hz, 8 Hz], the corresponding range of the alpha band is (8 Hz, 13 Hz], the corresponding range of the low beta band is (13 Hz, 22 Hz], the corresponding range of the high beta band is (22 Hz, 35 Hz], and the corresponding range of the gamma band is [60 Hz, 90 Hz].

[0036] Furthermore, in step 3, the nonlinearity, sharpness and average power characteristics of the IMF components corresponding to each frequency band are calculated, wherein the nonlinearity is calculated as follows:

[0037] A1) By performing Hilbert transform on the IMF signal, the phase function is obtained Then, the derivative of θ(t) is taken to calculate the instantaneous frequency IF of the IMF signal, that is, IF = dθ(t) / dt.

[0038] A2) Estimate the average frequency IF between adjacent zero points Z : Wherein T1 and T2 are adjacent zero-crossing time points.

[0039] A3) Considering the influence of local amplitude modulation on nonlinear estimation, calculate the average amplitude a between adjacent zero-crossing time points Z ,Right now

[0040] A4) Combining the standard deviation, which is commonly used to estimate the degree of data dispersion, we can obtain the degree of nonlinearity DoN as:

[0041]

[0042] Furthermore, in step 3, the nonlinearity, sharpness and average power characteristics of the IMF components corresponding to each frequency band are calculated, wherein the sharpness is calculated as follows:

[0043] B1) Find the zero point of the IMF component corresponding to each frequency band, and determine whether the signal is in the rising or falling phase based on the slope at the zero point.

[0044] B2) The time point corresponding to the maximum voltage between the zero point of the rising phase and the zero point of the subsequent falling phase is defined as the peak, and the time point corresponding to the minimum voltage between the zero point of the falling phase and the zero point of the subsequent rising phase is defined as the trough.

[0045] B3) If the oscillation reaches a peak or trough at time T, the voltage at the peak or trough is expressed as V T , the voltage t milliseconds before the peak or trough is expressed as V T-t The voltage t milliseconds after the peak or trough is expressed as V T+t , then the sharpness is expressed as:

[0046]

[0047] Furthermore, in step 3, the nonlinearity, sharpness, and average power characteristics of the IMF components corresponding to each frequency band are calculated, wherein the average power characteristics are calculated as follows:

[0048] C1) Use the Welch average power diagram method to obtain the power spectral density of the IMF component.

[0049] C2) For each IMF component, based on its power spectral density distribution in the frequency domain, integrate the power spectral density in the six frequency bands (delta band, theta band, alpha band, low beta band, high beta band, and gamma band) to obtain the total power value of the corresponding frequency band.

[0050] C3) Divide the total power value of each frequency band by the frequency band width to obtain the average power of the corresponding frequency band.

[0051] Furthermore, step 4: integrating the nonlinearity, sharpness and average power features, and using a machine learning classification model to identify effective motor imagery events, specifically:

[0052] For the EEG signal sequence x(t), the average power, sharpness and nonlinearity of I channels corresponding to six frequency bands were obtained.

[0053] The minimum redundancy-maximum relevance algorithm was used to sort the features extracted from all 6×3×I groups. The top K features in the weighted order were selected and input into various machine learning classification models to detect motor imagery events within a time span of m seconds.

[0054] Applicable machine learning classification models include linear regression model, boosted tree model, naive Bayes model or decision tree model.

[0055] Another embodiment of the present invention provides a motor imagery EEG signal processing system based on adaptive waveform features, including an EEG signal acquisition module, a frequency component extraction module, a feature calculation module, and a motor imagery time recognition module.

[0056] The EEG signal acquisition module is used to collect original multi-channel EEG signals, which include parts marked as related to motor imagery; preprocess the original multi-channel EEG signals; extract data segments marked as related to motor imagery training from the preprocessed signals, and splice the data segments by channel to obtain multidimensional EEG signals.

[0057] The frequency component extraction module uses an adaptive decomposition method to decompose each channel signal in the multidimensional EEG signal to obtain the IMF component, and divides each IMF component into the corresponding frequency band according to its frequency characteristics.

[0058] The feature calculation module calculates the nonlinearity, sharpness and average power characteristics of the IMF components corresponding to each frequency band, and obtains the average power, sharpness and nonlinearity of each channel corresponding to each frequency band.

[0059] The motor imagery time recognition module is used to integrate nonlinearity, sharpness and average power features and use a machine learning classification model to identify valid motor imagery events.

[0060] The specific processing process of the EEG signal acquisition module is as follows:

[0061] Preprocess each channel data in the original multi-channel EEG signal one by one; the specific preprocessing process is as follows:

[0062] The average value of all electrodes used is used as the reference to eliminate errors introduced by variations in the original reference electrode.

[0063] Interpolation is performed on sample points whose absolute amplitude exceeds the set threshold.

[0064] The interpolated signal is band-pass filtered using a notch filter.

[0065] The pre-processed signal is segmented, and the m-second data segments marked as related to motor imagery training are extracted and spliced ​​by channel to obtain the time series set x(t)={x i (t), i=1,2,…,I}, where I is the number of channels.

[0066] Frequency component extraction module, the specific processing process is:

[0067] S201: length is L, Represents K corresponding angles on the (I-1)-dimensional unit ball space The direction vector set of ; Represents the 1st to Ith direction vectors; Indicates the 1st to I-1th corresponding angles.

[0068] S202: Using the Hammersley sequence sampling method, a uniform sampling point set is obtained on the (I-1)-dimensional sphere to obtain an I-dimensional spatial direction vector.

[0069] S203: Initialization: Margin

[0070] S204: Extract the hth multivariate intrinsic mode function IMF h (t), the specific steps are:

[0071] a) Set y0(t) = r h-1 (t),p=1.

[0072] b) y p-1 (t) is the input signal, in each direction vector The mapping on Where · represents the vector inner product, * represents the vector scalar product, Represents a vector Modulus.

[0073] c) Confirm The instantaneous moment corresponding to the extreme point c is the number of extreme points.

[0074] d) Yes The extreme point series in Perform multivariate spline interpolation to obtain K corresponding multivariate envelopes

[0075] e) Calculate the mean of the extreme value envelope based on the center of gravity of the extreme value point

[0076] f) Calculate d(t) = y p-1 (t)-m(t), if the mean of d(t) is 0, the h-th multivariate intrinsic mode function IMF is obtained h (t)=d(t), otherwise, y p (t)=d(t), p increases by 1, and b) to f) are repeated.

[0077] S205: From signal r h-1 (t) minus IMF h (t) Get the new margin r h (t) = r h-1 (t)-IMF h(t); If the remainder is less than a given threshold, the algorithm ends and all multivariate intrinsic mode functions and residual components are obtained. Otherwise, h is incremented by 1 and S204) to S205) are repeated.

[0078] Finally, the original signal x(t) becomes the multivariate intrinsic mode function IMF h (t) and the remainder r(t) are of the form H is the number of eigenmode functions.

[0079] S206: Perform frequency division on the above IMF components, find the center frequency of each component, and classify it into six frequency bands: delta band, theta band, alpha band, low beta band, high beta band, and gamma band. The corresponding range of the delta band is [0.5 Hz, 4 Hz], the corresponding range of the theta band is (4 Hz, 8 Hz], the corresponding range of the alpha band is (8 Hz, 13 Hz], the corresponding range of the low beta band is (13 Hz, 22 Hz], the corresponding range of the high beta band is (22 Hz, 35 Hz], and the corresponding range of the gamma band is [60 Hz, 90 Hz].

[0080] Beneficial effects:

[0081] The present invention provides a method and system for processing motor imagery EEG signals based on adaptive waveform features. By incorporating an adaptive decomposition method into the nonlinear waveform and average power of motor imagery EEG signals, the system uses these features as key features for evaluating effective motor imagery. This method aims to improve the accuracy and reliability of waveform analysis while preserving the variability of neural oscillation waveforms, successfully quantifying the nonlinear characteristics of neural oscillation waveforms. This helps improve the identification of effective motor imagery signals. BRIEF DESCRIPTION OF THE DRAWINGS

[0082] Figure 1 This is a block diagram of the adaptive motor imagery EEG signal processing method and system principle based on nonlinear waveform characteristics provided by the present invention. DETAILED DESCRIPTION

[0083] The present invention is described in detail below with reference to the accompanying drawings and embodiments.

[0084] In order to quantify the waveform characteristics of brain waves during motor imagery, this embodiment proposes a method and system for adaptive motor imagery EEG signal processing based on nonlinear waveform characteristics. Figure 1 This is a technical flow chart of the present invention, wherein steps 1 to 3 are the EEG signal acquisition module, step 4 is the frequency component extraction module, steps 5 to 7 are the nonlinear waveform and power feature calculation module, and steps 8 and 9 are the motor imagery event recognition module. The detailed process is as follows:

[0085] Step 1: Acquire original multi-channel EEG signals, which contain parts marked as being related to motor imagery.

[0086] Preprocess the original multi-channel EEG signals.

[0087] Data segments marked as related to motor imagery training are extracted from the preprocessed signals, and the data segments are spliced ​​by channels to obtain multidimensional EEG signals.

[0088] In an embodiment of the present invention, the subject performs motor imagination by watching video prompts and records corresponding labels (such as imagining right hand movement or imagining left hand movement). At the same time, the EEG device collects the subject's multi-channel cortical EEG or deep EEG signals, with the sampling frequency set to fs, and stores them in a computer.

[0089] Preprocess the data for each channel of the above EEG signal one by one. The average value used by all electrodes is used as a reference to eliminate the error caused by changes in the original reference electrode. Through visual inspection, delete signal segments with particularly large and overly dense noise interference, and interpolate sample points with absolute amplitudes exceeding a threshold (e.g., 100 μV). Use a 50 / 60 Hz notch filter to eliminate power frequency interference, and perform bandpass filtering between 0.5 and 120 Hz to remove EEG signal drift caused by factors such as scalp sweat, while retaining information in the frequency range of interest.

[0090] The pre-processed signal is segmented, and the m-second data segments marked as motor imagery training are extracted and spliced ​​according to the channels to obtain the multi-dimensional EEG signal x(t) = {x i (t), i=1,2,…,I}, where I is the number of channels.

[0091] Step 2: For each channel signal in the multidimensional EEG signal, an adaptive decomposition method is used to decompose the IMF components and divide them into corresponding frequency bands according to the frequency characteristics of each IMF component.

[0092] The present invention uses adaptive decomposition methods, such as integrated empirical mode decomposition, multivariate empirical mode decomposition, and noise-assisted multivariate empirical mode decomposition, for each channel signal in the multidimensional EEG signal x(t), and divides it into corresponding frequency bands according to the frequency characteristics of each component. Taking the multivariate empirical mode decomposition method as an example, the specific steps are as follows:

[0093] S201: I-dimensional vector sequence Represents the EEG signal of I channels, with a length of L, Represents K corresponding angles on the (I-1)-dimensional unit ball space Direction Vector Set.

[0094] S202: Using the Hammersley sequence sampling method, a uniform sampling point set is obtained on the (I-1)-dimensional sphere to obtain an I-dimensional spatial direction vector.

[0095] S203: Initialization: Margin

[0096] S204: Extract the hth multivariate intrinsic mode function IMF h (t):

[0097] a) Set y0(t) = r h-1 (t),p=1.

[0098] b) y p-1 (t) is the input signal, in each direction vector Mapping on Where · represents the vector inner product, * represents the vector scalar product, Represents a vector Modulus.

[0099] c) Confirm The instantaneous moment corresponding to the extreme point c is the number of extreme points.

[0100] d) Yes The extreme point series in Perform multivariate spline interpolation to obtain K corresponding multivariate envelopes

[0101] e) Calculate the mean of the extreme value envelope based on the center of gravity of the extreme value point

[0102] f) Calculate d(t) = y p-1 (t)-m(t), if the mean of d(t) is 0, the h-th multivariate intrinsic mode function IMF is obtained h (t)=d(t), otherwise, y p (t)=d(t), p=p+1, repeat steps b) to f).

[0103] S205: From signal r h-1 (t) minus IMF h (t) Get the new margin r h (t) = r h-1 (t)-IMF h (t); If the residual is less than a given threshold (such as 0.075), the algorithm ends and all multivariate intrinsic mode functions and residual components are obtained. Otherwise, let h = h + 1 and repeat steps 4) to 5). Finally, the original signal x(t) becomes a multivariate intrinsic mode function IMF. h (t) and the remainder r(t) are of the form H is the number of eigenmode functions.

[0104] S206: Frequency division is performed on the above-mentioned IMF components to find the center frequency of each component and classify it into six frequency bands: delta band, theta band, alpha band, low beta band, high beta band, and gamma band. The corresponding range of the delta band is [0.5 Hz, 4 Hz], the corresponding range of the theta band is (4 Hz, 8 Hz], the corresponding range of the alpha band is (8 Hz, 13 Hz], the corresponding range of the low beta band is (13 Hz, 22 Hz], the corresponding range of the high beta band is (22 Hz, 35 Hz], and the corresponding range of the gamma band is [60 Hz, 90 Hz]. For example, the power spectrum S(f) of each IMF component is calculated, the energy of each frequency component is multiplied by its frequency, the products are added, and finally divided by the total energy to obtain the weighted average center frequency.

[0105] Step 3: Calculate the nonlinearity, sharpness, and average power characteristics of the IMF components corresponding to each frequency band to obtain the average power, sharpness, and nonlinearity of each channel corresponding to each frequency band.

[0106] The specific steps for calculating the nonlinearity of the frequency IMF corresponding to each channel signal are as follows:

[0107] A1) By performing Hilbert transform on the IMF signal, the phase function is obtained Then, the derivative of θ(t) is taken to calculate the instantaneous frequency IF of the IMF signal, that is, IF = dθ(t) / dt.

[0108] A2) Estimate the average frequency between adjacent zero points (IF Z ): Use the zero-crossing method to estimate the average frequency between adjacent zero points Wherein T1 and T2 are adjacent zero-crossing time points.

[0109] A3) Consider the impact of local amplitude modulation on nonlinear estimation: Similar to the zero-point method, calculate the average amplitude a between adjacent zero-crossing time points Z ,Right now

[0110] A4) Combining the standard deviation, which is commonly used to estimate the degree of data dispersion, we obtain the formula for calculating the degree of nonlinearity (DoN):

[0111]

[0112] The specific steps for calculating the sharpness of the IMF corresponding to the frequency of each channel signal are as follows:

[0113] B1) Find the zero point of the IMF component corresponding to each frequency band, and determine whether the signal is in the rising or falling phase based on the slope at the zero point.

[0114] B2) The time point corresponding to the maximum voltage between the zero point of the rising phase and the zero point of the subsequent falling phase is defined as the peak. Similarly, the time point corresponding to the minimum voltage between the zero point of the falling phase and the zero point of the subsequent rising phase is defined as the trough.

[0115] B3) If the oscillation reaches a peak or trough at time T, the voltage at the peak or trough is expressed as V T The voltage t milliseconds (t defaults to 5) before the peak or trough is expressed as V T-t The voltage t milliseconds after the peak or trough is expressed as V T+t Since the extreme sharpness is independent of polarity, and polarity is affected by the polarization direction related to the electrode placement, the average sharpness of the trough and peak is calculated. Then the sharpness can be expressed as

[0116]

[0117] The specific steps for calculating the average power characteristics of the frequency IMF corresponding to each channel signal are as follows:

[0118] C1) Using the Welch average power diagram method (existing method), the power spectral density of the IMF component is obtained.

[0119] C2) For each IMF component, based on its power spectral density distribution in the frequency domain, integrate the power spectral density in the six frequency bands (delta band, theta band, alpha band, low beta band, high beta band, and gamma band) to obtain the total power value of the corresponding frequency band.

[0120] C3) Divide the total power value of each frequency band by the frequency band width to obtain the average power of the corresponding frequency band.

[0121] Step eight, after the calculations in steps five to seven, the average power, sharpness, and nonlinearity of the six frequency bands corresponding to the I channel are obtained for the EEG signal sequence x(t). Next, the minimum redundancy-maximum correlation algorithm is used to sort the features extracted from all 6*3*I (frequency band*average power / sharpness / nonlinearity*channel) groups, and the top K features (default is 20) in the weighted order are selected and input into a variety of machine learning classification models to detect motor imagery events within a time length of m seconds. Applicable machine learning classification models include linear regression models, boosting tree models, naive Bayes models, decision trees, and other models. These models will be used to identify motor imagery events.

[0122] To execute the above method, another embodiment of the present invention further provides a motor imagery EEG signal processing system based on adaptive waveform features, which is characterized in that it includes an EEG signal acquisition module, a frequency component extraction module, a feature calculation module and a motor imagery time recognition module.

[0123] The EEG signal acquisition module is used to collect original multi-channel EEG signals, which include parts marked as related to motor imagery; preprocess the original multi-channel EEG signals; extract data segments marked as related to motor imagery training from the preprocessed signals, and splice the data segments by channel to obtain multidimensional EEG signals.

[0124] The frequency component extraction module uses an adaptive decomposition method to decompose each channel signal in the multidimensional EEG signal to obtain the IMF component, and divides each IMF component into the corresponding frequency band according to its frequency characteristics.

[0125] The feature calculation module calculates the nonlinearity, sharpness and average power characteristics of the IMF components corresponding to each frequency band, and obtains the average power, sharpness and nonlinearity of each channel corresponding to each frequency band.

[0126] The motor imagery time recognition module is used to integrate nonlinearity, sharpness and average power features and use a machine learning classification model to identify valid motor imagery events.

[0127] The specific processing process of the EEG signal acquisition module is as follows:

[0128] Preprocess each channel data in the original multi-channel EEG signal one by one; the specific preprocessing process is as follows:

[0129] The average value of all electrodes used is used as the reference to eliminate errors introduced by variations in the original reference electrode.

[0130] Interpolation is performed on sample points whose absolute amplitude exceeds the set threshold.

[0131] The interpolated signal is band-pass filtered using a notch filter.

[0132] The pre-processed signal is segmented, and the m-second data segments marked as related to motor imagery training are extracted and spliced ​​by channel to obtain the time series set x(t)={x i (t), i=1,2,…,I}, where I is the number of channels.

[0133] Frequency component extraction module, the specific processing process is:

[0134] S201: length is L, Represents K corresponding angles on the (I-1)-dimensional unit ball space The direction vector set of ; Represents the 1st to Ith direction vectors; Indicates the 1st to I-1th corresponding angles.

[0135] S202: Using the Hammersley sequence sampling method, a uniform sampling point set is obtained on the (I-1)-dimensional sphere to obtain an I-dimensional spatial direction vector.

[0136] S203: Initialization: Margin

[0137] S204: Extract the hth multivariate intrinsic mode function IMF h (t), the specific steps are:

[0138] a) Set y0(t) = r h-1 (t),p=1.

[0139] b) y p-1 (t) is the input signal, in each direction vector The mapping on Where · represents the vector inner product, * represents the vector scalar product, Represents a vector Modulus.

[0140] c) Confirm The instantaneous moment corresponding to the extreme point c is the number of extreme points.

[0141] d) Yes The extreme point series in Perform multivariate spline interpolation to obtain K corresponding multivariate envelopes

[0142] e) Calculate the mean of the extreme value envelope based on the center of gravity of the extreme value point

[0143] f) Calculate d(t) = y p-1 (t)-m(t), if the mean of d(t) is 0, the h-th multivariate intrinsic mode function IMF is obtained h (t)=d(t), otherwise, y p (t)=d(t), p increases by 1, and b) to f) are repeated.

[0144] S205: From signal r h-1 (t) minus IMF h (t) Get the new margin r h (t) = r h-1 (t)-IMF h(t); If the remainder is less than a given threshold, the algorithm ends and all multivariate intrinsic mode functions and residual components are obtained. Otherwise, h is incremented by 1 and S204) to S205) are repeated.

[0145] Finally, the original signal x(t) becomes the multivariate intrinsic mode function IMF h (t) and the remainder r(t) are of the form H is the number of eigenmode functions.

[0146] S206: Perform frequency division on the above IMF components, find the center frequency of each component, and classify it into six frequency bands: delta band, theta band, alpha band, low beta band, high beta band, and gamma band. The corresponding range of the delta band is [0.5 Hz, 4 Hz], the corresponding range of the theta band is (4 Hz, 8 Hz], the corresponding range of the alpha band is (8 Hz, 13 Hz], the corresponding range of the low beta band is (13 Hz, 22 Hz], the corresponding range of the high beta band is (22 Hz, 35 Hz], and the corresponding range of the gamma band is [60 Hz, 90 Hz].

[0147] The results of this system show that by obtaining the EEG signals of the target user within the corresponding time window and identifying that the EEG signals to be decoded are all valid motor imagery, the accuracy of subsequent EEG signal decoding can be improved.

[0148] In summary, the above are only preferred embodiments of the present invention and are not intended to limit the scope of protection of the present invention. Any modifications, equivalent replacements, improvements, etc. made within the spirit and principles of the present invention should be included in the scope of protection of the present invention.

Claims

1. A method for processing motor imagery EEG signals based on adaptive waveform features, characterized in that: The steps include: Step 1: collecting original multi-channel EEG signals, wherein the original multi-channel EEG signals include a portion marked as being related to motor imagery; Preprocessing the original multi-channel EEG signal; Extracting data segments marked as relevant to motor imagery training from the preprocessed signals, and splicing the data segments by channel to obtain multidimensional EEG signals; Step 2: For each channel signal in the multidimensional EEG signal, an adaptive decomposition method is used to decompose the IMF components, and the IMF components are divided into corresponding frequency bands according to their frequency characteristics. Step 2 is specifically divided into the following steps: S201: Length L , Represented on the (I-1)-dimensional unit ball space K Corresponding angles The direction vector set of ; Represents the 1st to Ith direction vector; Indicates the 1st to I-1th corresponding angles; S202: Using the Hammersley sequence sampling method, a uniform sampling point set is obtained on the (I-1)-dimensional sphere, and we get I dimensional space direction vector; S203: Initialization: Margin ; S204: Extract h Multivariate Intrinsic Mode Function , the specific steps are: a) Settings ; b) is the input signal, in each direction vector The mapping on , where · represents the vector inner product, * represents the vector scalar product, Represents a vector Take the mold; c) Determine The instantaneous moment corresponding to the extreme point , c is the number of extreme points; d) Yes The extreme point series in Perform multivariate spline interpolation to obtain K The corresponding multivariate envelope ; e) Calculate the mean of the extreme value envelope based on the center of gravity of the extreme value point ; f) Calculation ,like The mean is 0, and the h Multivariate Intrinsic Mode Function ,otherwise, , p increases by 1, repeat b)~f); S205: From signal Subtract Get a new margin ; If the residual is less than the given threshold, the algorithm ends and all multivariate intrinsic mode functions and residual components are obtained. Otherwise, let h Increment by 1, repeat S204)~S205); Finally, the original signal x(t) Become a multivariate intrinsic mode function IMF h ( t ) and margin r ( t ) and the form , H is the number of eigenmode functions; S206: performing frequency division on the above IMF components, finding the center frequency of each component, and classifying it into six frequency bands: delta band, theta band, alpha band, low beta band, high beta band, and gamma band. The corresponding range of the delta band is [0.5 Hz, 4 Hz], the corresponding range of the theta band is (4 Hz, 8 Hz], the corresponding range of the alpha band is (8 Hz, 13 Hz], the corresponding range of the low beta band is (13 Hz, 22 Hz], the corresponding range of the high beta band is (22 Hz, 35 Hz], and the corresponding range of the gamma band is [60 Hz, 90 Hz]; Step 3: Calculate the nonlinearity, sharpness, and average power characteristics of the IMF components corresponding to each frequency band to obtain the average power, sharpness, and nonlinearity of each channel corresponding to each frequency band; Step 4: Integrate the nonlinearity, sharpness, and average power features and use a machine learning classification model to identify effective motor imagery events. Step 4 specifically includes: For EEG signal sequences x ( t ), the average power, sharpness and nonlinearity of the six frequency bands corresponding to one channel were obtained; Use the minimum redundancy-maximum correlation algorithm to sort all 6×3×I groups of extracted features, select the top K features in the weight sorting, input multiple machine learning classification models, and perform m Detection of motor imagery events within a time span of seconds; Applicable machine learning classification models include linear regression model, boosted tree model, naive Bayes model or decision tree model.

2. The method for processing motor imagery EEG signals based on adaptive waveform features according to claim 1, wherein: The step 1: pre-processing the original multi-channel EEG signal, the specific process is: Preprocess each channel data in the original multi-channel EEG signal one by one; the specific preprocessing process is: The average value used by all electrodes is used as the reference to eliminate the error caused by the variation of the original reference electrode; Perform interpolation processing on sample points whose absolute amplitude exceeds the set threshold; The interpolated signal is band-pass filtered using a notch filter.

3. The method for processing motor imagery EEG signals based on adaptive waveform features according to claim 1, wherein: The data segments marked as related to motor imagery training are extracted from the pre-processed signals, and the data segments are spliced ​​according to channels to obtain multi-dimensional EEG signals, specifically: The pre-processed signal is segmented to extract the signals marked as relevant to motor imagery training. m The data segments of 1 second are collected and spliced ​​by channel to obtain the time series set x(t)={xi(t), i = 1, 2, …,I} marked as related to motor imagery training. I is the number of channels.

4. The method for processing motor imagery EEG signals based on adaptive waveform features according to claim 3, wherein: The adaptive decomposition method is an integrated empirical mode decomposition method and a multivariate empirical mode decomposition method.

5. The method for processing motor imagery EEG signals based on adaptive waveform features according to claim 1, wherein: Step 3: Calculate the nonlinearity, sharpness and average power characteristics of the IMF components corresponding to each frequency band, wherein the nonlinearity is calculated in the following manner: A1) By performing Hilbert transform on the IMF signal, we can obtain the phase function , and then θ ( t ) to calculate the instantaneous frequency of the IMF signal IF ,Right now IF = d θ ( t ) / d t ; A2) Estimate the average frequency between adjacent zero points IFZ : ;in T 1 and T 2 is the adjacent zero-crossing time point; A3) Considering the impact of local amplitude modulation on nonlinear estimation, calculate the average amplitude between adjacent zero-crossing time points ,Right now ; A4) Combined with the standard deviation commonly used to estimate the degree of data dispersion, the degree of nonlinearity DoN is obtained as: 。 6. The method for processing motor imagery EEG signals based on adaptive waveform features according to claim 1, wherein: Step 3: Calculate the nonlinearity, sharpness, and average power characteristics of the IMF components corresponding to each frequency band, wherein the method for calculating the sharpness is specifically as follows: B1) Find the zero point of the IMF component corresponding to each frequency band and determine whether the signal is in the rising or falling phase based on the slope at the zero point; B2) The time point corresponding to the maximum voltage between the zero point of the rising phase and the zero point of the subsequent falling phase is defined as the peak, and the time point corresponding to the minimum voltage between the zero point of the falling phase and the zero point of the subsequent rising phase is defined as the trough; B3) If the oscillation reaches a peak or trough at time 𝑇, the voltage at the peak or trough is expressed as VT , before the peak or trough t The voltage in milliseconds is expressed as VT−t , the voltage t milliseconds after the peak or trough is expressed as VT+t , then the sharpness is expressed as: （2）。 7. The method for processing motor imagery EEG signals based on adaptive waveform features according to claim 5, wherein: Step 3: Calculate the nonlinearity, sharpness, and average power characteristics of the IMF components corresponding to each frequency band, wherein the average power characteristics are calculated as follows: C1) Use the Welch average power diagram method to obtain the power spectral density of the IMF component; C2) For each IMF component, based on the distribution of its power spectral density in the frequency domain, integrate the power spectral density in the six frequency bands (delta band, theta band, alpha band, low beta band, high beta band, and gamma band) to obtain the total power value of the corresponding frequency band; C3) Divide the total power value of each frequency band by the frequency band width to obtain the average power of the corresponding frequency band.

8. A motor imagery EEG signal processing system based on adaptive waveform features, characterized in that: It includes an EEG signal acquisition module, a frequency component extraction module, a feature calculation module, and a motor imagery time recognition module; The EEG signal acquisition module is used to collect original multi-channel EEG signals, wherein the original multi-channel EEG signals include parts marked as being related to motor imagery; and pre-process the original multi-channel EEG signals; Extracting data segments marked as relevant to motor imagery training from the preprocessed signals, and splicing the data segments by channel to obtain multidimensional EEG signals; The frequency component extraction module decomposes each channel signal in the multidimensional EEG signal using an adaptive decomposition method to obtain IMF components, and divides each IMF component into corresponding frequency bands according to its frequency characteristics; The feature calculation module calculates the nonlinearity, sharpness and average power characteristics of the IMF components corresponding to each frequency band to obtain the average power, sharpness and nonlinearity of each channel corresponding to each frequency band; The motor imagery time recognition module is used to integrate nonlinearity, sharpness and average power features and use a machine learning classification model to identify effective motor imagery events; The EEG signal acquisition module specifically processes as follows: Preprocess each channel data in the original multi-channel EEG signal one by one; the specific preprocessing process is: The average value used by all electrodes is used as the reference to eliminate the error caused by the variation of the original reference electrode; Perform interpolation processing on sample points whose absolute amplitude exceeds the set threshold; The interpolated signal is band-pass filtered using a notch filter; The pre-processed signal is segmented to extract the signals marked as relevant to motor imagery training. m Seconds data segments are spliced ​​by channel to obtain a set of time series marked as related to motor imagery training x(t)= { xi ( t ), i = 1, 2, …, I }, I is the number of channels; The specific processing process of the frequency component extraction module is as follows: S201: Length L , Indicates that ( I -1) dimensional unit sphere space K Corresponding angles The direction vector set of ; Represents the 1st to Ith direction vector; Indicates the 1st to I-1th corresponding angles; S202: Using Hammersley sequence sampling method, in ( I -1) Get a uniform sampling point set on the spherical surface, and get I dimensional space direction vector; S203: Initialization: Margin ; S204: Extract h Multivariate Intrinsic Mode Function , the specific steps are: a) Settings ; b) is the input signal, in each direction vector The mapping on , where · represents the vector inner product, * represents the vector scalar product, Represents a vector Take the mold; c) Determine The instantaneous moment corresponding to the extreme point , c is the number of extreme points; d) Yes The extreme point series in Perform multivariate spline interpolation to obtain K The corresponding multivariate envelope ; e) Calculate the mean of the extreme value envelope based on the center of gravity of the extreme value point ; f) Calculation ,like The mean is 0, and the h Multivariate Intrinsic Mode Function ,otherwise, , p increases by 1, repeat b)~f); S205: From signal Subtract Get a new margin ; If the residual is less than the given threshold, the algorithm ends and all multivariate intrinsic mode functions and residual components are obtained. Otherwise, let h Increment by 1, repeat S204)~S205); Finally, the original signal x(t) Become a multivariate intrinsic mode function IMF h ( t ) and margin r ( t ) and the form , H is the number of eigenmode functions; S206: Frequency division is performed on the above IMF components to find the center frequency of each component and classify it into six frequency bands: delta band, theta band, alpha band, low beta band, high beta band, and gamma band. The corresponding range of the delta band is [0.5 Hz, 4 Hz], the corresponding range of the theta band is (4 Hz, 8 Hz], the corresponding range of the alpha band is (8 Hz, 13 Hz], the corresponding range of the low beta band is (13 Hz, 22 Hz], the corresponding range of the high beta band is (22 Hz, 35 Hz], and the corresponding range of the gamma band is [60 Hz, 90 Hz].

Images