Classified electroencephalogram artifact removal method and system based on multivariable empirical mode decomposition fusion

Abstract

The invention provides a classified electroencephalogram artifact removing method and system based on multivariable empirical mode decomposition fusion, and relates to the technical field of electroencephalogram signal denoising. Firstly, multi-channel original electroencephalogram signals are obtained, multivariable empirical mode decomposition is conducted on the multi-channel original electroencephalogram signals through an optimized multivariable empirical mode decomposition algorithm, and the multi-channel original electroencephalogram signals are obtained; according to the method, the intrinsic mode functions of the multi-channel original electroencephalogram signals are obtained through more accurate decomposition, the decomposition efficiency is improved, and the feature information of the multi-channel original electroencephalogram signals is reserved, and then the intrinsic mode functions of the multi-channel original electroencephalogram signals are classified based on a multi-dimensional interference classification threshold value. According to the method, electroencephalogram artifacts are removed from intrinsic mode functions of electroencephalogram signals under signals of different interference types by adopting different denoising methods, so that denoising is more targeted, the denoising effect is better, and finally, signal reconstruction fusion is performed on the intrinsic mode functions of the electroencephalogram signals from which the signals of different interference types are removed. The electroencephalogram signals without the electroencephalogram artifacts are obtained, and the accuracy and reliability of analyzing the electroencephalogram signals are improved.

Description

Technical Field

[0001] The present invention relates to the technical field of electroencephalogram (EEG) signal denoising, and more specifically, to a method and system for removing EEG artifacts by classification based on multivariate empirical mode decomposition and fusion. Background Art

[0002] With the advancement of science and the rapid development of biomedical technology, biological signals have been applied to research in various fields. Among them, EEG signals, as a biological signal, have been widely used in sleep classification, epilepsy, concentration analysis, and diagnosis of depression.

[0003] Electroencephalography (EEG) is a method of collecting brain electrical signals by placing electrodes on the head. It is widely used due to its high temporal resolution and non-invasive advantages. However, EEG signals themselves have the characteristics of low frequency and nonlinearity, and are easily interfered by electrooculogram and electromyography signals, which makes it difficult to obtain pure EEG, affecting clinical diagnosis and related experimental research. Therefore, it is of great significance to study an efficient denoising technology to remove the interference of electrooculogram and electromyography signals in EEG signals.

[0004] There are currently two methods for the study of interference removal of EEG signals: the first is the multivariate empirical mode decomposition method (MEEMD) based on the empirical mode decomposition method (EMD). This method decomposes the EEG signal to obtain a number of intrinsic mode functions (IMFs), which can effectively extract different frequency components. The decomposed intrinsic mode functions reflect the essential characteristics of the original EEG signal at different time scales and frequencies, and the decomposed intrinsic mode functions are arranged from high frequency to low frequency. In the decomposition process, the correlation between multi-channel EEG signals is fully utilized to enhance the accuracy of EEG signal decomposition. However, the existing MEEMD requires a reasonable selection of the number of decomposition layers and the penalty factor. The wrong selection of these two parameters will lead to unsatisfactory decomposition effects, which in turn leads to low efficiency and unsatisfactory denoising effects of EEG signals. The second is the independent component analysis (ICA) and classical component analysis (CCA) based on the blind source separation (BSS) method. These methods use the independent characteristics of source signals to separate noise signals from source signals, but have application limitations: ICA assumes that each component is independent of each other, but in actual applications, some source signals may be correlated, and this assumption is not always true. CCA focuses on the correlation between a single independent variable and a single dependent variable, without considering the correlation between the variables within the independent variable and dependent variable groups. This may cause the linear model to be unable to fully capture the true characteristics of the signal when processing complex EEG signals.

[0005] On December 12, 2023, China Invention Patent Publication No. (CN117216475A) announced a method and system for denoising EEG signals. In this scheme, the variational modal decomposition model is first optimized using the dung beetle optimization algorithm, and then the original EEG signal is decomposed using the optimized variational modal decomposition model. The modal component after decomposition is then calculated to determine whether the modal component is valid by calculating the correlation coefficient between the modal component after decomposition and the original EEG signal. Finally, a clean EEG signal is obtained by reconstructing the effective modal component and denoising with the wavelet threshold. However, this scheme uses the same denoising algorithm for different interference signals, and different interference signals have different characteristics. Using the same denoising method, the denoising is poorly targeted and the denoising effect is poor. Summary of the invention

[0006] In order to solve the problems of low denoising efficiency, poor denoising pertinence and poor denoising effect of current methods for EEG signal denoising, the present invention proposes a hierarchical EEG artifact removal method based on a multivariate empirical mode decomposition fusion algorithm, which combines the Grey Wolf optimization algorithm with the MEEMD to obtain the optimal number K of modal components and the optimal penalty factor α, sets a multidimensional classification threshold to accurately separate interference signals from EEG signals, uses different denoising algorithms for different types of interference signals, improves denoising pertinence and reduces computational complexity.

[0007] In order to achieve the above technical effects, the technical solution of the present invention is as follows:

[0008] In the first aspect, the present application proposes a classification and removal method of EEG artifacts by multivariate empirical mode decomposition fusion, comprising the following steps:

[0009] S1. Obtain multi-channel raw EEG signals;

[0010] S2. using the optimized multivariate empirical mode decomposition algorithm to perform multivariate empirical mode decomposition on the multi-channel original EEG signal to obtain the intrinsic mode function of the original EEG signal;

[0011] S3. Classifying the intrinsic mode function based on the multi-dimensional interference classification threshold to obtain the intrinsic mode function of the EEG signal under different interference type signals;

[0012] S4. using different denoising methods to remove EEG artifacts from the intrinsic mode functions of EEG signals under different interference types of signals, and obtaining the intrinsic mode functions of EEG signals from which different interference types of signals have been removed;

[0013] S5. Perform signal reconstruction and fusion on the intrinsic mode functions of the EEG signals from which signals of different interference types have been removed, to obtain the EEG signals from which EEG artifacts have been removed.

[0014] In the present technical scheme, firstly, a multi-channel original EEG signal is acquired, and then a multi-variable empirical mode decomposition is performed on the multi-channel original EEG signal using an optimized multivariate empirical mode decomposition algorithm, so as to more accurately decompose and obtain the intrinsic mode function of the multi-channel original EEG signal, thereby improving the decomposition efficiency and retaining the characteristic information of the multi-channel original EEG signal. Then, the intrinsic mode function of the multi-channel original EEG signal is classified based on the multi-dimensional interference classification threshold, and different denoising methods are used to remove EEG artifacts from the intrinsic mode function of the EEG signal under different interference type signals, so that the denoising is more targeted and the denoising effect is better. Finally, the intrinsic mode function of the EEG signal from which the signals of different interference types have been removed is reconstructed and fused to obtain the EEG signal from which the EEG artifacts have been removed, thereby improving the accuracy and reliability of the analysis of the EEG signal.

[0015] Preferably, the process of obtaining multi-channel original EEG signals in step S1 is:

[0016] S11. Place multiple electrodes at different locations on the subject's scalp;

[0017] S12. Use electrodes to record data on the electrical activities of multiple areas or multiple neurons in the subject's brain, and organize them into multi-channel original EEG signals of the subject's brain.

[0018] Preferably, the multivariable empirical mode decomposition algorithm is optimized. Specifically, the multivariable empirical mode decomposition algorithm is optimized using the gray wolf optimization algorithm. The process is:

[0019] The ranges of the number of modal components K and the penalty factor α in the multivariate empirical mode decomposition algorithm are defined as:

[0020] K min ≤K≤K max

[0021] α min ≤α≤α max

[0022] Among them, K min Indicates the minimum value of the number of modal components, 1, K max represents the maximum number of modal components, α min Represents the minimum value of the penalty factor, α max Indicates the maximum value of the penalty factor;

[0023] The number of modal components K and the penalty factor α are combined to generate P candidate solutions, expressed as:

[0024] Y i =[K i ,α i ],i=1,2,…,P

[0025] Among them, Y i Represents the combination of the number of modal components K and the penalty factor α;

[0026] The objective function is constructed to measure the kurtosis of each modal component, expressed as:

[0027]

[0028] Where f(α, K) represents the objective function used to calculate the kurtosis of each modal component, x j is the modal component IMF i The jth sample point of is the modal component IMF i The mean of the modal component is iThe standard deviation of , N is the number of samples of the modal component, and K is the number of modal components;

[0029] like If the minimum value is obtained, that is, the optimal kurtosis is obtained, then the candidate solution is the optimal solution. According to the expression of the optimal solution, the optimal number of modal components K and the optimal penalty factor α are determined. Otherwise, the candidate solution is updated and iterated. The update and iteration expression is:

[0030] Y i new =Y α -A·|C·Y α -Y i |

[0031] Among them, Y i new represents the updated candidate solution, Y α represents the reference solution under the combination of the number of current modal components K and the penalty factor α, A represents the coefficient of the control update amplitude, and C represents a coefficient of the distance from the current solution, which affects the direction and moving speed of the candidate solution. The calculation expressions of A and C are:

[0032] A=2a·r1-a,C=2·r2

[0033] Among them, r1 and r2 are random numbers in the range of [0,1], and a decreases linearly from 2 to 0 as the current number of iterations t increases. The expression of a is:

[0034]

[0035] Among them, a0 is the initial value of a, T max is the maximum number of iterations. When the iterative update reaches the maximum number of iterations, the candidate solution after iterative update is obtained and taken as the optimal solution. The optimal number K of modal components and the optimal penalty factor α are determined according to the expression of the optimal solution.

[0036] According to the above technical means, the multivariate empirical mode decomposition algorithm is optimized, the decomposition efficiency of multi-channel original EEG signals is improved, and the optimal number K of modal components and the optimal penalty factor α are obtained.

[0037] Preferably, the specific process of performing multivariate integrated empirical mode decomposition on multi-channel original EEG signals using the optimized multivariate empirical mode decomposition algorithm is as follows:

[0038] S21. Decompose each channel signal in the multi-channel EEG signal into the corresponding intrinsic mode function IMF ij , the expression is:

[0039]

[0040] Among them, x i (t) is the i-th channel signal to be decomposed; N i is the number of IMFs of the ith channel, that is, the optimal number of modal components K obtained after optimization by the Grey Wolf Optimization Algorithm, R i (t) is the residual signal associated with the current intrinsic mode function, α is the optimal penalty factor α obtained after optimization by the gray wolf optimization algorithm;

[0041] S22. The intrinsic mode functions of all channels are statistically combined to obtain the expression of the set of each intrinsic mode function:

[0042]

[0043] Where M represents the number of channels in the multi-channel EEG signal.

[0044] According to the above technical means, the efficiency of decomposing the intrinsic mode function of the multi-channel original EEG signal is improved, and the accuracy and robustness of the decomposition of the multi-channel original EEG signal are improved.

[0045] Preferably, before classifying the intrinsic mode function based on the multidimensional interference classification threshold, a multidimensional interference classification threshold is set, and the specific process is:

[0046] Perform fast Fourier transform on each intrinsic mode function to obtain the spectrum representation within the intrinsic mode function, which is expressed as:

[0047] S(f)=FFT(IMF i )

[0048] Among them, S(f) is the spectrum representation after fast Fourier transform;

[0049] Calculate the frequency at which the spectrum representation reaches its maximum value. The expression is:

[0050] f max =argmax|S(f)|

[0051] The spectrum of the eigenmode function is calculated to represent the energy E in the frequency range of the electrooculogram signal. eye , the expression is:

[0052]

[0053] The spectrum of the eigenmode function is calculated to represent the energy E in the frequency range of the electromyographic signal. muscle , the expression is:

[0054]

[0055] The spectrum of the eigenmode function is calculated to represent the energy E in the frequency range of the electrooculogram signal. eye The spectrum of the eigenmode function represents the energy E in the frequency range of the electromyographic signal muscle The ratio G is expressed as:

[0056]

[0057] According to the obtained frequency f max and energy ratio G, setting the classification threshold parameters, the classification threshold parameters including setting the energy ratio threshold as G th and the frequency threshold is f th If f max <f th , and G>G th , then the frequency f max The corresponding intrinsic mode function contains the electrooculographic interference signal that needs to be removed; if f max >f th , and G <G th , then the frequency f max The corresponding intrinsic mode function contains the electromyographic interference signal that needs to be removed.

[0058] According to the above technical means, the accuracy of classifying the intrinsic mode function is improved.

[0059] Preferably, after calculating the frequency corresponding to the maximum value of the spectrum representation, the calculated frequency is preprocessed to remove signals with a frequency of more than 50 Hz.

[0060] Preferably, the electrooculographic interference signal is removed using an independent component analysis method, and the removal process is as follows:

[0061] The intrinsic mode function IMF containing the electrooculographic interference signal eye It is expressed as:

[0062] IMF eye =[imf1(t),imf2(t),…,imf n (t)] T

[0063] Where n is the number of intrinsic mode functions classified as containing electrooculographic interference signals that need to be removed;

[0064] IMF is analyzed using independent component analysis algorithm. eye The signal matrix is ​​decomposed into independent components S and mixing matrix A:

[0065] IMF eye =A·S

[0066] Among them, the independent component S contains the electrooculographic interference signal and the EEG signal;

[0067] Calculate the information entropy H(S) of each independent component S i ), the expression is:

[0068]

[0069] Among them, P(S i =imf) is the ith independent component S containing the intrinsic mode function of the electrooculographic interference signal to be removed i The probability distribution of

[0070] Sort the information entropy of all independent components to obtain the sorted list H of information entropy sorted The expression is:

[0071] H sorted =sort(H(S1),H(S2),…,H(S m ))

[0072] Set the information entropy threshold H threshold , if the value of information entropy is greater than the information entropy threshold, the corresponding independent component is removed; if the information entropy is less than or equal to the information entropy threshold, the corresponding independent component is retained:

[0073]

[0074] S selected1 ={S i |H(S i ) <H threshold}

[0075] Among them, m is the total number of independent components, S selected1 To remove the independent components after electrooculographic interference signals;

[0076] After removing the electrooculographic interference signal, the independent component S selected1 Reconstructed with the mixing matrix into the intrinsic mode function with the electrooculographic interference signal removed The expression is:

[0077]

[0078] Among them, A -1 is the inverse matrix of the mixing matrix.

[0079] According to the above technical means, the electrooculographic interference signals in the multi-channel original EEG signals are removed in a targeted manner, thereby improving the denoising efficiency and accuracy.

[0080] Preferably, a typical component analysis method is used to remove the electromyographic interference signal, and the specific removal process is as follows:

[0081] The intrinsic mode function IMF containing the electromyographic interference signal miscle It is expressed as:

[0082] IMF muscle =[IMF1(t),IMF2(t),…,IMF p (t)] T

[0083] Where p is the total number of eigenmode functions containing myoelectric interference signals;

[0084] Compute the correlation coefficient for each canonical correlation component:

[0085]

[0086] Among them, Cov(IMF j ,M q ) represents the covariance between the jth intrinsic mode function containing myoelectric noise and the reference myoelectric signal M(t), Var(IMF j ) represents the variance of the jth intrinsic mode function containing myoelectric noise, Var(M q ) represents the variance of the reference electromyographic signal M(t), r q represents the correlation coefficient of the qth typical correlation component;

[0087] Calculate the mean μ of the canonical correlation coefficient R and standard deviation σ R , the expression is:

[0088]

[0089] Based on the mean and standard deviation of the typical correlation coefficient, the threshold T for removing the electromyographic interference signal is set. If the typical correlation coefficient of the intrinsic mode function containing the electromyographic interference signal is greater than the threshold, it is the intrinsic mode function containing the electromyographic interference signal that needs to be removed:

[0090] T=μ R +α·σ R

[0091] Among them, α is the adjustment factor;

[0092] Remove the intrinsic mode function containing the electromyographic interference signal and reconstruct the signal into the intrinsic mode function with the electromyographic interference signal removed

[0093]

[0094] in, is a set of intrinsic mode functions containing the electromyographic interference signal that needs to be removed.

[0095] According to the above technical means, the myoelectric interference signals in the multi-channel original EEG signals are removed in a targeted manner, thereby improving the denoising efficiency and accuracy.

[0096] Preferably, the intrinsic mode function of the electrooculographic interference signal and the electromyographic interference signal is subjected to signal reconstruction fusion to obtain the electroencephalogram signal X from which the electroencephalogram artifacts have been removed. clean , whose expression is:

[0097]

[0098] Where M is the number of channels.

[0099] In the second aspect, the present application also proposes a classification and removal of EEG artifacts system based on multivariate empirical mode decomposition fusion, the system comprising:

[0100] A signal acquisition module is used to acquire multi-channel raw EEG signals;

[0101] A multivariate empirical mode decomposition module is used to perform multivariate empirical mode decomposition on multi-channel original EEG signals using an optimized multivariate empirical mode decomposition algorithm to obtain the intrinsic mode function of the original EEG signals;

[0102] Interference classification module, which classifies the intrinsic mode function based on the multi-dimensional interference classification threshold to obtain the intrinsic mode function of the EEG signal under different interference type signals;

[0103] A denoising module is used to remove EEG artifacts by using different denoising methods on the Eigenmode functions of EEG signals under different interference types, and obtain the Eigenmode functions of EEG signals from which the different interference types of signals have been removed respectively;

[0104] The reconstruction fusion module is used to reconstruct and fuse the intrinsic mode functions of the EEG signals from which signals of different interference types have been removed, so as to obtain the EEG signals from which EEG artifacts have been removed.

[0105] Compared with the prior art, the technical solution of the present invention has the following beneficial effects:

[0106] The present invention proposes a method and system for removing EEG artifacts by classification based on multivariate empirical mode decomposition and fusion. First, a multi-channel original EEG signal is obtained, and a multivariate empirical mode decomposition is performed on the multi-channel original EEG signal using an optimized multivariate empirical mode decomposition algorithm to more accurately decompose and obtain the intrinsic mode function of the multi-channel original EEG signal, thereby improving the decomposition efficiency and retaining the characteristic information of the multi-channel original EEG signal. Then, the intrinsic mode function of the multi-channel original EEG signal is classified based on a multidimensional interference classification threshold, and different denoising methods are used to remove EEG artifacts from the intrinsic mode function of the EEG signal under different interference type signals, so that denoising is more targeted and the denoising effect is better. Finally, the intrinsic mode function of the EEG signal from which the different interference type signals have been removed is reconstructed and fused to obtain the EEG signal from which the EEG artifacts have been removed, thereby improving the accuracy and reliability of analyzing the EEG signal. BRIEF DESCRIPTION OF THE DRAWINGS

[0107] Figure 1 A schematic diagram showing a flow chart of a method for removing EEG artifacts by classification based on multivariate empirical mode decomposition fusion proposed in Embodiment 1 of the present invention;

[0108] Figure 2 A schematic diagram showing a flow chart of a multivariable empirical mode decomposition optimization algorithm proposed in Embodiment 2 of the present invention;

[0109] Figure 3 A schematic diagram showing a process of setting a multi-dimensional interference classification threshold value proposed in Embodiment 2 of the present invention;

[0110] Figure 4 A structural diagram showing the system for classifying and removing EEG artifacts based on multivariate empirical mode decomposition fusion proposed in Example 3 of the present invention. DETAILED DESCRIPTION

[0111] The drawings are for illustrative purposes only and should not be construed as limiting the present patent;

[0112] In order to better illustrate the present embodiment, some parts of the drawings may be omitted, enlarged or reduced, and do not represent the actual size;

[0113] It is understandable to those skilled in the art that descriptions of certain well-known contents in the drawings may be omitted.

[0114] The technical solution of the present invention is further described below in conjunction with the accompanying drawings and embodiments.

[0115] The positional relationships described in the drawings are only for illustrative purposes and should not be construed as limiting the present patent;

[0116] Example 1

[0117] This embodiment proposes a classification method for removing EEG artifacts based on multivariate empirical mode decomposition fusion. The flowchart of this method is shown in Figure 1 , including the following steps:

[0118] S1. Obtain multi-channel raw EEG signals;

[0119] S2. using the optimized multivariate empirical mode decomposition algorithm to perform multivariate empirical mode decomposition on the multi-channel original EEG signal to obtain the intrinsic mode function of the original EEG signal;

[0120] S3. Classifying the intrinsic mode function based on the multi-dimensional interference classification threshold to obtain the intrinsic mode function of the EEG signal under different interference type signals;

[0121] S4. Different denoising methods are used to remove EEG artifacts from the intrinsic mode functions of EEG signals under different interference types, and the intrinsic mode functions of EEG signals with different interference types removed are obtained respectively.

[0122] In this embodiment, firstly, a multi-channel original EEG signal is obtained, and then the multi-channel original EEG signal is subjected to multi-variable empirical mode decomposition using an optimized multi-variable empirical mode decomposition algorithm, so that the intrinsic mode function of the multi-channel original EEG signal can be more accurately decomposed, thereby improving the efficiency of decomposition and retaining the characteristic information of the multi-channel original EEG signal. Then, the intrinsic mode function of the multi-channel original EEG signal is classified based on the multi-dimensional interference classification threshold, and different denoising methods are used to remove EEG artifacts from the intrinsic mode function of the EEG signal under different interference type signals, so that denoising is more targeted and the denoising effect is better. Finally, the intrinsic mode function of the EEG signal from which the signals of different interference types have been removed is subjected to signal reconstruction and fusion to obtain an EEG signal from which EEG artifacts have been removed. The entire technical process scheme adopted in this embodiment improves the accuracy and reliability of analyzing EEG signals.

[0123] Example 2

[0124] In this embodiment, the process of obtaining multi-channel original EEG signals is as follows:

[0125] S11. Place multiple electrodes at different locations on the subject's scalp;

[0126] S12. Use electrodes to record data on the electrical activities of multiple areas or multiple neurons in the subject's brain, and organize them into multi-channel original EEG signals of the subject's brain.

[0127] Specifically, the process of step S12 is:

[0128] After recording is completed, the recorded data is exported from the recording device and saved in a computer;

[0129] Use professional EEG signal processing software to pre-process and analyze the data, including filtering, artifact removal, segmentation and other operations;

[0130] The preprocessed data is organized into a multi-channel raw EEG signal format.

[0131] In this embodiment, the gray wolf optimization algorithm is used to optimize the multivariable empirical mode decomposition algorithm. For the optimization flow chart, see Figure 2 , the process is:

[0132] The ranges of the number of modal components K and the penalty factor α in the multivariate empirical mode decomposition algorithm are defined as:

[0133] K min ≤K≤K max

[0134] α min ≤α≤α max

[0135] Among them, K min Represents the minimum number of modal components, K max represents the maximum number of modal components, α min Represents the minimum value of the penalty factor, α max Indicates the maximum value of the penalty factor;

[0136] Specifically, K min The value of is generally 2; K max The value of is generally 10; min The value of is usually 1000; max The value is generally 3000.

[0137] The number of modal components K and the penalty factor α are combined to generate P candidate solutions, expressed as:

[0138] Y i =[K i ,α i ],i=1,2,…,P

[0139] Among them, Y i Represents the combination of the number of modal components K and the penalty factor α;

[0140] The objective function is constructed to measure the kurtosis of each modal component, expressed as:

[0141]

[0142] Where f(α, K) represents the objective function used to calculate the kurtosis of each modal component, x j is the modal component IMF iThe jth sample point of is the modal component IMF i The mean of the modal component is i The standard deviation of , N is the number of samples of the modal component, and K is the number of modal components;

[0143] like The minimum value 1 is obtained, that is, the optimal kurtosis is obtained, then the candidate solution is the optimal solution, and the optimal number K of modal components and the optimal penalty factor α are determined according to the expression of the optimal solution. Otherwise, the candidate solution is updated and iterated, and the update and iteration expression is:

[0144] Y i new =Y α -A·|C·Y α -Y i |

[0145] Among them, Y i new represents the updated candidate solution, Y α represents the reference solution under the combination of the number of current modal components K and the penalty factor α, A represents the coefficient of the control update amplitude, and C represents a coefficient of the distance from the current solution, which affects the direction and moving speed of the candidate solution. The calculation expressions of A and C are:

[0146] A=2a·r1-a,C=2·r2

[0147] Among them, r1 and r2 are random numbers in the range of [0,1], and a decreases linearly from 2 to 0 as the current number of iterations t increases. The expression of a is:

[0148]

[0149] Among them, a0 is the initial value of a, T max is the maximum number of iterations. When the iterative update reaches the maximum number of iterations, the candidate solution after iterative update is obtained and taken as the optimal solution. The optimal number K of modal components and the optimal penalty factor α are determined according to the expression of the optimal solution.

[0150] In this embodiment, the specific process of performing multivariate integrated empirical mode decomposition on multi-channel original EEG signals using the optimized multivariate empirical mode decomposition algorithm is as follows:

[0151] S21. Decompose each channel signal in the multi-channel EEG signal into the corresponding intrinsic mode function IMF ij , the expression is:

[0152]

[0153] Among them, xi (t) is the i-th channel signal that needs to be decomposed; N i is the number of IMFs of the ith channel, that is, the optimal number of modal components K obtained after optimization by the Grey Wolf Optimization Algorithm, R i (t) is the residual signal associated with the current intrinsic mode function, α is the optimal penalty factor α obtained after optimization by the gray wolf optimization algorithm;

[0154] S22. The intrinsic mode functions of all channels are statistically combined to obtain the expression of the set of each intrinsic mode function:

[0155]

[0156] Where M represents the number of channels in the multi-channel EEG signal.

[0157] In this embodiment, before classifying the intrinsic mode function based on the multidimensional interference classification threshold, a multidimensional interference classification threshold is set. The flowchart of this step is shown in FIG. Figure 3 , the specific process is:

[0158] Perform fast Fourier transform on each intrinsic mode function to obtain the spectrum representation within the intrinsic mode function, which is expressed as:

[0159] S(f)=FFT(IMF i )

[0160] Among them, S(f) is the spectrum representation after fast Fourier transform;

[0161] Calculate the frequency at which the spectrum representation reaches its maximum value. The expression is:

[0162] f max =argmax|S(f)|

[0163] The spectrum of the eigenmode function is calculated to represent the energy E in the frequency range of the electrooculogram signal. eye , the expression is:

[0164]

[0165] The spectrum of the eigenmode function is calculated to represent the energy E in the frequency range of the electromyographic signal. muscle , the expression is:

[0166]

[0167] The spectrum of the eigenmode function is calculated to represent the energy E in the frequency range of the electrooculogram signal. eye The spectrum of the eigenmode function represents the energy E in the frequency range of the electromyographic signal muscle The ratio G is expressed as:

[0168]

[0169] According to the obtained frequency f max and energy ratio G, setting the classification threshold parameters, the classification threshold parameters including setting the energy ratio threshold as G th =1 and the frequency threshold is f th =10Hz; if f max <10Hz, and G>1, then the frequency f max The corresponding intrinsic mode function contains the electrooculographic interference signal that needs to be removed; if f max >10Hz, and G<1, then the frequency f max The corresponding intrinsic mode function contains the electromyographic interference signal that needs to be removed.

[0170] In this embodiment, after calculating the frequency corresponding to the maximum value of the spectrum representation, the calculated frequency is preprocessed to remove signals with a frequency of more than 50 Hz.

[0171] In this embodiment, different denoising methods are used to remove EEG artifacts from the intrinsic mode functions of EEG signals under different interference types. Specifically, the independent component analysis method is used to remove the electrooculographic interference signal. The removal process is as follows:

[0172] The intrinsic mode function IMF containing the electrooculographic interference signal eye It is expressed as:

[0173] IMF eye =[imf1(t),imf2(t),…,imf n (t)] T

[0174] Where n is the number of intrinsic mode functions classified as containing electrooculographic interference signals that need to be removed;

[0175] IMF eye The signal matrix is ​​decomposed into independent components S and mixing matrix A:

[0176] IMF eye =A·S

[0177] Among them, the independent component S contains the electrooculographic interference signal and the EEG signal;

[0178] Calculate the information entropy H(S) of each independent component S i ), the expression is:

[0179]

[0180] Among them, P(S i=imf) is the ith independent component S containing the intrinsic mode function of the electrooculographic interference signal to be removed i The probability distribution of

[0181] Sort the information entropy of all independent components to obtain the sorted list H of information entropy sorted The expression is:

[0182] H sorted =sort(H(S1),H(S2),…,H(S m ))

[0183] Set the information entropy threshold H threshold , if the value of information entropy is greater than the information entropy threshold, the corresponding independent component is removed; if the information entropy is less than or equal to the information entropy threshold, the corresponding independent component is retained:

[0184]

[0185] S selected1 ={S i |H(S i ) <H threshold}

[0186] Among them, m is the total number of independent components, S selected1 To remove the independent components after electrooculographic interference signals;

[0187] After removing the electrooculographic interference signal, the independent component S selected1 Reconstructed with the mixing matrix into the intrinsic mode function with the electrooculographic interference signal removed The expression is:

[0188]

[0189] Among them, A -1 is the inverse matrix of the mixing matrix.

[0190] In this embodiment, different denoising methods are used to remove EEG artifacts from the intrinsic mode functions of EEG signals under different interference types. Specifically, the typical component analysis method is used to remove the EMG interference signal. The removal process is as follows:

[0191] The intrinsic mode function IMF containing the electromyographic interference signal muscle It is expressed as:

[0192] IMF muscle =[IMF1(t),IMF2(t),…,IMF p (t)] T

[0193] Where p is the total number of eigenmode functions containing myoelectric interference signals;

[0194] Compute the correlation coefficient for each canonical correlation component:

[0195]

[0196] Among them, Cov(IMF j ,M q ) represents the covariance between the jth intrinsic mode function containing myoelectric noise and the reference myoelectric signal M(t), Var(IMF j ) represents the variance of the jth intrinsic mode function containing myoelectric noise, Var(M q ) represents the variance of the reference electromyographic signal M(t), r q represents the correlation coefficient of the qth typical correlation component;

[0197] Calculate the mean μ of the canonical correlation coefficient R and standard deviation σ R , the expression is:

[0198]

[0199] Based on the mean and standard deviation of the typical correlation coefficient, the threshold T for removing the electromyographic interference signal is set to 0.5. If the typical correlation coefficient of the intrinsic mode function containing the electromyographic interference signal is greater than the threshold, it is the intrinsic mode function containing the electromyographic interference signal that needs to be removed:

[0200] T=μ R +α·σ R

[0201] Among them, α is the adjustment factor;

[0202] Remove the intrinsic mode function containing the electromyographic interference signal and reconstruct the signal into the intrinsic mode function with the electromyographic interference signal removed

[0203]

[0204] in, is a set of intrinsic mode functions containing the electromyographic interference signal that needs to be removed.

[0205] In this embodiment, the intrinsic mode functions of the EEG signals from which different types of interference signals have been removed are reconstructed and fused to obtain the EEG signal X from which EEG artifacts have been removed. clean , whose expression is:

[0206]

[0207] Where M is the number of channels.

[0208] Example 3

[0209] This embodiment proposes a classification and removal of EEG artifacts system based on multivariate empirical mode decomposition fusion. Figure 4 , the system comprising:

[0210] A signal acquisition module is used to acquire multi-channel raw EEG signals;

[0211] A multivariate empirical mode decomposition module is used to perform multivariate empirical mode decomposition on multi-channel original EEG signals using an optimized multivariate empirical mode decomposition algorithm to obtain the intrinsic mode function of the original EEG signals;

[0212] Interference classification module, which classifies the intrinsic mode function based on the multi-dimensional interference classification threshold to obtain the intrinsic mode function of the EEG signal under different interference type signals;

[0213] A denoising module is used to remove EEG artifacts by using different denoising methods on the Eigenmode functions of EEG signals under different interference types, and obtain the Eigenmode functions of EEG signals from which the different interference types of signals have been removed respectively;

[0214] The reconstruction fusion module is used to reconstruct and fuse the intrinsic mode functions of the EEG signals from which signals of different interference types have been removed, so as to obtain the EEG signals from which EEG artifacts have been removed.

[0215] In this embodiment, in the signal acquisition module, the process of acquiring multi-channel original EEG signals is as follows:

[0216] S11. Place multiple electrodes at different locations on the subject's scalp;

[0217] S12. Use electrodes to record data on the electrical activities of multiple areas or multiple neurons in the subject's brain, and organize them into multi-channel original EEG signals of the subject's brain.

[0218] Specifically,

[0219] After recording is completed, the recorded data is exported from the recording device and saved in a computer;

[0220] Use professional EEG signal processing software to pre-process and analyze the data, including filtering, artifact removal, segmentation and other operations;

[0221] The preprocessed data is organized into a multi-channel raw EEG signal format.

[0222] In this embodiment, in the multivariate empirical mode decomposition module, the gray wolf optimization algorithm is used to optimize the multivariate empirical mode decomposition algorithm. For the optimization flow chart, see Figure 2 , the process is:

[0223] The ranges of the number of modal components K and the penalty factor α in the multivariate empirical mode decomposition algorithm are defined as:

[0224] K min ≤K≤K max

[0225] α min ≤α≤α max

[0226] Among them, K min Represents the minimum number of modal components, K max Represents the maximum number of modal components, α min Represents the minimum value of the penalty factor, α max Indicates the maximum value of the penalty factor;

[0227] Specifically, K min The value of is generally 2; K max The value of is usually 10; min The value of is usually 1000; max The value is generally 3000.

[0228] The number of modal components K and the penalty factor α are combined to generate P candidate solutions, expressed as:

[0229] Y i =[K i ,α i ],i=1,2,…,P

[0230] Among them, Y i Represents the combination of the number of modal components K and the penalty factor α;

[0231] The objective function is constructed to measure the kurtosis of each modal component, expressed as:

[0232]

[0233] Where f(α, K) represents the objective function used to calculate the kurtosis of each modal component, x j is the modal component IMF i The jth sample point of is the modal component IMF i The mean of the modal component is i The standard deviation of , N is the number of samples of the modal component, and K is the number of modal components;

[0234] like The minimum value 1 is obtained, that is, the optimal kurtosis is obtained, then the candidate solution is the optimal solution, and the optimal number K of modal components and the optimal penalty factor α are determined according to the expression of the optimal solution. Otherwise, the candidate solution is updated and iterated, and the update and iteration expression is:

[0235] Y i new =Y α -A·|C·Y α -Y i |

[0236] Among them, Y i new represents the updated candidate solution, Y α represents the reference solution under the combination of the number of current modal components K and the penalty factor α, A represents the coefficient of the control update amplitude, and C represents a coefficient of the distance from the current solution, which affects the direction and moving speed of the candidate solution. The calculation expressions of A and C are:

[0237] A=2a·r1-a,C=2·r2

[0238] Among them, r1 and r2 are random numbers in the range of [0,1], and a decreases linearly from 2 to 0 as the current number of iterations t increases. The expression of a is:

[0239]

[0240] Among them, a0 is the initial value of a, T max is the maximum number of iterations. When the iterative update reaches the maximum number of iterations, the candidate solution after iterative update is obtained and taken as the optimal solution. The optimal number K of modal components and the optimal penalty factor α are determined according to the expression of the optimal solution.

[0241] In this embodiment, in the multivariate empirical mode decomposition module, the specific process of using the optimized multivariate empirical mode decomposition algorithm to perform multivariate integrated empirical mode decomposition on the multi-channel original EEG signal is as follows:

[0242] Decompose each channel signal in the multi-channel EEG signal into the corresponding intrinsic mode function IMF ij , the expression is:

[0243]

[0244] Among them, x i (t) is the i-th channel signal that needs to be decomposed; N i is the number of IMFs of the ith channel, that is, the optimal number of modal components K obtained after optimization by the Grey Wolf Optimization Algorithm, R i(t) is the residual signal associated with the current intrinsic mode function, α is the optimal penalty factor α obtained after optimization by the gray wolf optimization algorithm;

[0245] The intrinsic mode functions of all channels are statistically combined to obtain the set expression of each intrinsic mode function:

[0246]

[0247] Where M represents the number of channels in the multi-channel EEG signal.

[0248] In this embodiment, in the interference classification module, before classifying the intrinsic mode function based on the multidimensional interference classification threshold, a multidimensional interference classification threshold is set. The flowchart of this step is shown in FIG. Figure 3 , the specific process is:

[0249] Perform fast Fourier transform on each intrinsic mode function to obtain the spectrum representation within the intrinsic mode function, which is expressed as:

[0250] S(f)=FFT(IMF i )

[0251] Among them, S(f) is the spectrum representation after fast Fourier transform;

[0252] Calculate the frequency at which the spectrum representation reaches its maximum value. The expression is:

[0253] f max =argmax|S(f)|

[0254] The spectrum of the eigenmode function is calculated to represent the energy E in the frequency range of the electrooculogram signal. eye , the expression is:

[0255]

[0256] The spectrum of the eigenmode function is calculated to represent the energy E in the frequency range of the electromyographic signal. muscle , the expression is:

[0257]

[0258] The spectrum of the eigenmode function is calculated to represent the energy E in the frequency range of the electrooculogram signal. eye The spectrum of the eigenmode function represents the energy E in the frequency range of the electromyographic signal muscle The ratio G is expressed as:

[0259]

[0260] According to the obtained frequency f maxand energy ratio G, set the classification threshold parameters, the classification threshold parameters include setting the energy ratio threshold to G th =1 and the frequency threshold is f th =10Hz; if f max <10Hz, and G>1, then the frequency f max The corresponding intrinsic mode function contains the electrooculographic interference signal that needs to be removed; if f max >10Hz, and G<1, then the frequency f max The corresponding intrinsic mode function contains the electromyographic interference signal that needs to be removed.

[0261] In this embodiment, in the denoising module, after calculating the frequency corresponding to the maximum value of the spectrum representation, the calculated frequency is preprocessed to remove signals with a frequency of more than 50 Hz.

[0262] In this embodiment, in the denoising module, different denoising methods are used to remove EEG artifacts for the intrinsic mode functions of EEG signals under different interference types of signals. Specifically, the independent component analysis method is used to remove the electrooculographic interference signal. The removal process is as follows:

[0263] The intrinsic mode function IMF containing the electrooculographic interference signal eye It is expressed as:

[0264] IMF eye =[imf1(t),imf2(t),…,imf N (t)] T

[0265] Where n is the number of intrinsic mode functions classified as containing electrooculographic interference signals that need to be removed;

[0266] IMF eye The signal matrix is ​​decomposed into independent components S and mixing matrix A:

[0267] IMF eye =A·S

[0268] Among them, the independent component S contains the electrooculographic interference signal and the EEG signal;

[0269] Calculate the information entropy H(S) of each independent component S i ), the expression is:

[0270]

[0271] Among them, P(S i =imf) is the ith independent component S containing the intrinsic mode function of the electrooculographic interference signal to be removed i The probability distribution of

[0272] Sort the information entropy of all independent components to obtain the sorted list H of information entropy sorted The expression is:

[0273] H sorted =sort(H(S1),H(S2),…,H(S m ))

[0274] Set the information entropy threshold H threshold , if the value of information entropy is greater than the information entropy threshold, the corresponding independent component is removed; if the information entropy is less than or equal to the information entropy threshold, the corresponding independent component is retained:

[0275]

[0276] S selected1 ={S i |H(S i ) <H threshold}

[0277] Among them, m is the total number of independent components, S selected1 To remove the independent components after electrooculographic interference signals;

[0278] After removing the electrooculographic interference signal, the independent component S selected1 Reconstructed with the mixing matrix into the intrinsic mode function with the electrooculographic interference signal removed The expression is:

[0279]

[0280] Among them, A -1 is the inverse matrix of the mixing matrix.

[0281] In this embodiment, in the denoising module, different denoising methods are used to remove EEG artifacts for the intrinsic mode functions of EEG signals under different interference types of signals. Specifically, the typical component analysis method is used to remove the EMG interference signal. The removal process is as follows:

[0282] The intrinsic mode function IMF containing the electromyographic interference signal muscle It is expressed as:

[0283] IMF muscle =[IMF1(t),IMF2(t),…,IMF p (t)] T

[0284] Where p is the total number of eigenmode functions containing myoelectric interference signals;

[0285] Compute the correlation coefficient for each canonical correlation component:

[0286]

[0287] Among them, Cov(IMF j ,M q ) represents the covariance between the jth intrinsic mode function containing myoelectric noise and the reference myoelectric signal M(t), Var(IMF j ) represents the variance of the jth intrinsic mode function containing myoelectric noise, Var(M q ) represents the variance of the reference electromyographic signal M(t), r q represents the correlation coefficient of the qth typical correlation component;

[0288] Calculate the mean μ of the canonical correlation coefficient R and standard deviation σ R , the expression is:

[0289]

[0290] Based on the mean and standard deviation of the typical correlation coefficient, the threshold T for removing the electromyographic interference signal is set to 0.5. If the typical correlation coefficient of the intrinsic mode function containing the electromyographic interference signal is greater than the threshold, it is the intrinsic mode function containing the electromyographic interference signal that needs to be removed:

[0291] T=μ R +α·σ R

[0292] Among them, α is the adjustment factor;

[0293] Remove the intrinsic mode function containing the electromyographic interference signal and reconstruct the signal into the intrinsic mode function with the electromyographic interference signal removed

[0294]

[0295] in, is a set of intrinsic mode functions containing the electromyographic interference signal that needs to be removed.

[0296] In this embodiment, in the reconstruction fusion module, the intrinsic mode functions of the EEG signals from which the signals of different interference types have been removed are subjected to signal reconstruction fusion to obtain the EEG signal X from which the EEG artifacts have been removed. clean , whose expression is:

[0297]

[0298] Where M is the number of channels.

[0299] Obviously, the above embodiments of the present invention are only examples for clearly illustrating the present invention, and are not intended to limit the embodiments of the present invention. For those skilled in the art, other different forms of changes or modifications can be made based on the above description. It is not necessary and impossible to list all the embodiments here. Any modifications, equivalent substitutions and improvements made within the spirit and principles of the present invention should be included in the protection scope of the claims of the present invention.

Claims

1. A classification method for removing EEG artifacts based on multivariate empirical mode decomposition fusion, characterized in that: The following steps are involved: S1. Obtain multi-channel raw EEG signals; S2. using the optimized multivariate empirical mode decomposition algorithm to perform multivariate empirical mode decomposition on the multi-channel original EEG signal to obtain the intrinsic mode function of the original EEG signal; S3. Classifying the intrinsic mode function based on the multi-dimensional interference classification threshold to obtain the intrinsic mode function of the EEG signal under different interference type signals; S4. using different denoising methods to remove EEG artifacts from the intrinsic mode functions of EEG signals under different interference types of signals, and obtaining the intrinsic mode functions of EEG signals from which different interference types of signals have been removed; S5. Reconstruct and fuse the intrinsic mode functions of the EEG signals from which signals of different interference types have been removed to obtain the EEG signals from which EEG artifacts have been removed.

2. The method for removing EEG artifacts by classification based on multivariate empirical mode decomposition fusion according to claim 1 is characterized in that: The process of obtaining multi-channel original EEG signals in step S1 is as follows: S11. Place multiple electrodes at different locations on the subject's scalp; S12. Use electrodes to record data on the electrical activities of multiple areas or multiple neurons in the subject's brain, and organize them into multi-channel original EEG signals of the subject's brain.

3. The method for removing EEG artifacts by classification based on multivariate empirical mode decomposition fusion according to claim 1, characterized in that: The Grey Wolf Optimization Algorithm is used to optimize the multivariable empirical mode decomposition algorithm. The process is as follows: The ranges of the number of modal components K and the penalty factor α in the multivariate empirical mode decomposition algorithm are defined as: K min ≤K≤K max α min ≤α≤α max Among them, K min Represents the minimum number of modal components, K max represents the maximum number of modal components, α min Represents the minimum value of the penalty factor, α max Indicates the maximum value of the penalty factor; The number of modal components K and the penalty factor α are combined to generate P candidate solutions, expressed as: Y i =[K i ,α i ],i=1,2,...,P Among them, Y i Represents the combination of the number of modal components K and the penalty factor α; The objective function is constructed to measure the kurtosis of each modal component, expressed as: Among them, f(α, K) represents the objective function to measure the kurtosis of each modal component, x j is the modal component IMF i The jth sample point of is the modal component IMF i The mean of the modal component is i The standard deviation of , N is the number of samples of the modal component, and K is the number of modal components; like If the minimum value is obtained, that is, the optimal kurtosis is obtained, then the candidate solution is the optimal solution. According to the expression of the optimal solution, the optimal number of modal components K and the optimal penalty factor α are determined. Otherwise, the candidate solution is updated and iterated. The update and iteration expression is: in, represents the updated candidate solution, Y α represents the reference solution under the combination of the number of current modal components K and the penalty factor α, A represents the coefficient of the control update amplitude, and C represents a coefficient of the distance from the current solution, which affects the direction and moving speed of the candidate solution. The calculation expressions of A and C are: A=2a·r1-a,C=2·r2 Among them, r1 and r2 are random numbers in the range of [0, 1], and a decreases linearly from 2 to 0 as the current number of iterations t increases. The expression of a is: Among them, a0 is the initial value of a, T max is the maximum number of iterations. When the iterative update reaches the maximum number of iterations, the candidate solution after iterative update is obtained and taken as the optimal solution. The optimal number K of modal components and the optimal penalty factor α are determined according to the expression of the optimal solution.

4. The method for removing EEG artifacts by classification based on multivariate empirical mode decomposition fusion according to claim 3 is characterized in that: The specific process of using the optimized multivariate empirical mode decomposition algorithm to perform multivariate integrated empirical mode decomposition on multi-channel original EEG signals is as follows: S21. Decompose each channel signal in the multi-channel EEG signal into the corresponding intrinsic mode function IMF ij , the expression is: Among them, x i (t) represents the i-th channel signal to be decomposed; N i is the number of IMFs of the ith channel, that is, the optimal number of modal components K obtained after optimization by the Grey Wolf Optimization Algorithm, R i (t) is the residual signal associated with the current intrinsic mode function, α is the optimal penalty factor α obtained after optimization by the gray wolf optimization algorithm; S22. The intrinsic mode functions of all channels are statistically combined to obtain the expression of the set of each intrinsic mode function: Where M represents the number of channels in the multi-channel EEG signal.

5. The method for removing EEG artifacts by classification based on multivariate empirical mode decomposition fusion according to claim 4 is characterized in that: Before classifying the intrinsic mode function based on the multidimensional interference classification threshold, the multidimensional interference classification threshold is set. The specific process is as follows: Perform a fast Fourier transform on each intrinsic mode function to obtain the spectrum representation within the intrinsic mode function, which is expressed as: S(f)=FFT(IMF i ) Among them, S(f) is the spectrum representation after fast Fourier transform; Calculate the frequency at which the spectrum representation reaches its maximum value. The expression is: f max =arg max|S(f)| The spectrum of the eigenmode function is calculated to represent the energy E in the frequency range of the electrooculogram signal. eye , the expression is: The spectrum of the eigenmode function is calculated to represent the energy E in the frequency range of the electromyographic signal. muscle , the expression is: The spectrum of the eigenmode function is calculated to represent the energy E in the frequency range of the electrooculogram signal. eye The spectrum of the eigenmode function represents the energy E in the frequency range of the electromyographic signal muscle The ratio G is expressed as: According to the obtained frequency f max and energy ratio G, setting the classification threshold parameters, the classification threshold parameters including setting the energy ratio threshold as G th and the frequency threshold is f th If f max <f th , and G>G th , then the frequency f max The corresponding intrinsic mode function contains the electrooculographic interference signal that needs to be removed; if f max >f th , and G <G th , then the frequency f max The corresponding intrinsic mode function contains the electromyographic interference signal that needs to be removed.

6. The method for removing EEG artifacts by classification based on multivariate empirical mode decomposition fusion according to claim 5, characterized in that: After calculating the frequency corresponding to the maximum value of the spectrum representation, the calculated frequency is preprocessed to remove signals with a frequency of more than 50 Hz.

7. The method for removing EEG artifacts by classification based on multivariate empirical mode decomposition fusion according to claim 4, characterized in that: The independent component analysis method is used to remove the electrooculographic interference signal. The removal process is as follows: The intrinsic mode function IMF containing the electrooculographic interference signal eye It is expressed as: IMF eye =[imf1(t),imf2(t),...,imf n (t)] T Where n is the number of intrinsic mode functions classified as containing electrooculographic interference signals that need to be removed; IMF is analyzed using independent component analysis algorithm. eye The signal matrix is ​​decomposed into independent components S and mixing matrix A: IMF eye =A·S Among them, the independent component S contains the electrooculographic interference signal and the EEG signal; Calculate the information entropy H(S) of each independent component S i ), the expression is: Among them, P(S i =imf) is the ith independent component S containing the intrinsic mode function of the electrooculographic interference signal to be removed i The probability distribution of Sort the information entropy of all independent components to obtain the sorted list H of information entropy sorted The expression is: H sorted =sort(H(S1),H(S2),...,H(S m )) Let the information entropy threshold be H threshold , if the value of information entropy is greater than the information entropy threshold, the corresponding independent component is removed; if the information entropy is less than or equal to the information entropy threshold, the corresponding independent component is retained: S selected1 ={S i |H(S i )<H threshold } Among them, m is the total number of independent components, S selected1 To remove the independent components after electrooculographic interference signals; After removing the electrooculographic interference signal, the independent component S selected1 Reconstructed with the mixing matrix into the intrinsic mode function with the electrooculographic interference signal removed The expression is: Among them, A -1 is the inverse matrix of the mixing matrix.

8. The method for removing EEG artifacts by classification based on multivariate empirical mode decomposition fusion according to claim 7, characterized in that: The typical component analysis method is used to remove the electromyographic interference signal. The specific removal process is as follows: The intrinsic mode function IMF containing the electromyographic interference signal muscle It is expressed as: IMF muscle =[IMF1(t),IMF2(t),...,IMF p (t)] T Where p is the total number of eigenmode functions containing myoelectric interference signals; Compute the correlation coefficient for each canonical correlation component: Among them, Cov(IMF j , M q ) represents the covariance between the jth intrinsic mode function containing myoelectric noise and the reference myoelectric signal M(t), Var(IMF j ) represents the variance of the jth intrinsic mode function containing myoelectric noise, Var(M q ) represents the variance of the reference electromyographic signal M(t), r q represents the correlation coefficient of the qth typical correlation component; Calculate the mean μ of the canonical correlation coefficient R and standard deviation σ R , the expressions are: Set the threshold T for removing the electromyographic interference signal. If the typical correlation coefficient of the intrinsic mode function containing the electromyographic interference signal is greater than T, the intrinsic mode function containing the electromyographic interference signal to be removed is: T=μ R +a·s R Among them, α is the adjustment factor; Remove the intrinsic mode function containing the electromyographic interference signal and reconstruct the signal into the intrinsic mode function with the electromyographic interference signal removed in, is a set of intrinsic mode functions containing the electromyographic interference signal that needs to be removed.

9. The method for removing EEG artifacts by classification based on multivariate empirical mode decomposition fusion according to claim 8, characterized in that: The intrinsic mode functions of the removed electrooculographic interference signals and electromyographic interference signals are reconstructed and fused to obtain the EEG signal X with the EEG artifacts removed. clean , whose expression is: Where M is the number of channels.

10. A classification and removal of EEG artifacts system based on multivariate empirical mode decomposition fusion, characterized in that: The system is used to implement the method according to any one of claims 1 to 9, comprising: A signal acquisition module is used to acquire multi-channel raw EEG signals; A multivariate empirical mode decomposition module is used to perform multivariate empirical mode decomposition on multi-channel original EEG signals using an optimized multivariate empirical mode decomposition algorithm to obtain the intrinsic mode function of the original EEG signals; Interference classification module, which classifies the intrinsic mode function based on the multi-dimensional interference classification threshold to obtain the intrinsic mode function of the EEG signal under different interference type signals; A denoising module is used to remove EEG artifacts by using different denoising methods on the Eigenmode functions of EEG signals under different interference types, and obtain the Eigenmode functions of EEG signals from which the different interference types of signals have been removed respectively; The reconstruction fusion module is used to reconstruct and fuse the intrinsic mode functions of the EEG signals from which signals of different interference types have been removed, so as to obtain the EEG signals from which EEG artifacts have been removed.

Images