Epileptic signal classification method based on EWT and improved CSP

Abstract

The present invention proposes an epilepsy signal classification method based on EWT and improved CSP. The method comprises: preprocessing the signal, selecting the signal channel by correlation analysis, extracting the features of the three feature domains of time domain, frequency domain and space domain of the epileptic EEG signal in the interictal period and preictal period, and finally inputting the merged feature matrix of the three feature domains into the support vector machine classifier to achieve effective recognition of the interictal period and preictal period of the epileptic EEG signal. Among them, the time domain features include the root mean square value, the absolute mean value and the number of zero crossing points; the frequency domain features are extracted by using the empirical wavelet transform for multimodal decomposition and the Welch power spectrum for single mode selection, and then the Hilbert transform method is used for the single mode of each channel to extract the instantaneous amplitude and instantaneous frequency of the signal; the spatial domain features are extracted by using the improved CSP algorithm. In the present invention, the empirical wavelet transform overcomes the limitations of the EMD method, and in addition, the improved CSP algorithm has a high recognition rate and computational efficiency.

Description

Technical Field

[0001] The invention belongs to the field of pattern recognition, and relates to a method for recognizing an electroencephalogram (EEG) signal pattern, and in particular to a method for classifying epileptic EEG signals. Background Art

[0002] Epilepsy is one of the most common neurological diseases. It is a brain dysfunction disease caused by sudden abnormal discharge of brain neurons. Epilepsy attacks are generally very sudden, and the duration of each attack is different, so epilepsy attacks are highly uncertain. During an epileptic attack, the patient generally loses consciousness, convulses, foams at the mouth, and even threatens the patient's life. At present, according to a survey by the World Health Organization (WHO), there are about 50 million people suffering from epilepsy worldwide, of which 9 million are in China. Long-term and repeated epileptic attacks not only cause physical and mental torture to patients, but also bring huge economic losses. In addition, epileptic attacks can seriously affect patients' daily activities, such as driving outside, swimming, etc. At present, although most epilepsy patients can control and treat epilepsy through drugs and surgery, there are still some patients whose epilepsy cannot be cured. Therefore, if the characteristics of epilepsy signals before the attack can be discovered, epilepsy attacks can be predicted in advance, so that epilepsy patients can be reminded and protected in advance, which can not only avoid the harm caused by epileptic attacks to patients, but also help the prevention and treatment and research of epilepsy. At present, electroencephalogram (EEG) is one of the main tools for epilepsy prediction. In recent years, the use of EEG to process epileptic EEG signals has also been studied in depth. In 1997, Qu and Gotman proposed a nearest neighbor classification method based on EEG features in the time-frequency domain to detect epileptic seizures. In 2004, Gigola et al. used a method based on wavelet analysis and cumulative energy evolution to predict the occurrence of epileptic seizures in epileptic EEG data recorded by intracranial electrodes. Other researchers have proposed an epileptic seizure diagnosis system based on artificial neural networks (ANN). Weng and Khorasani proposed a method of using the average EEG amplitude and average EEG duration, variance coefficient, main frequency and average power spectrum as feature vectors to feed into an adaptive structured neural network. These methods mainly extract the features of the signal in the time-frequency domain, but the information of epileptic EEG signals is reflected in the three feature domains of time, frequency and space. Although the above methods have achieved good results in the recognition of epileptic EEG signals, they lack information about spatial features, and the lack of information will inevitably affect the results. Summary of the invention

[0003] In view of the shortcomings of the existing technology, in order to achieve the correct identification of EEG signals of epileptic patients, the present invention proposes an epilepsy signal classification method based on EWT and improved CSP. By extracting the fusion features of EEG signals in the time domain, frequency domain and spatial domain, various aspects of epileptic EEG signals are obtained to improve the recognition accuracy.

[0004] To achieve the above object of the invention, the method comprises the following steps:

[0005] Step 1: Data preprocessing;

[0006] The selected EEG data are divided into normal data and onset data based on the onset time; the data 30 to 120 minutes before the onset time is taken as normal data. The first half of the normal data is taken as the interictal data, and the second half is taken as the preictal data. After selecting the training set and the test set, the normal data in the set is segmented using a non-overlapping sliding window with a window size of a seconds, and the data of one window is taken as an experimental sample. The samples in the interictal period of the training set are marked as category 1, and the samples in the preictal period are marked as category 2.

[0007] Step 2: Select the EEG channel;

[0008] For the processed sample data, the EEG signals are calculated using the correlation calculation method, and the EEG channels in the sample are selected based on the calculation results. If the Pearson correlation coefficient between different channels is high, it means that there is information duplication between the two channels. The calculation formula of the Pearson correlation system is:

[0009]

[0010] cov(X,Y) is the covariance between two variables, σ X σ Y is the standard deviation between the two variables.

[0011] The correlation coefficient between two channels is obtained according to the calculation result of formula (1), and thus the correlation matrix is ​​obtained. A threshold b is set. If the correlation coefficient between two channels is greater than the threshold, it means they are correlated, otherwise it means they are uncorrelated. The EEG channels are selected according to the correlation coefficient. The selection process is:

[0012] (1) Select all channels that do not have related channels;

[0013] (2) For channels that are equivalent to each other, the channel with the smaller channel number is selected;

[0014] (3) Select the channel that contains the most channels from the remaining channels;

[0015] Step 3: Extract time domain features;

[0016] Extract the root mean square value, absolute mean value and number of zero crossings of the EEG signal of the channel selected in step 2 as the time domain features of the signal, and obtain the time domain training set matrix f 1train and the time domain test set matrix f 1test ; The calculation formula of the root mean square value is:

[0017]

[0018] The absolute mean value is calculated as:

[0019]

[0020] The formula for calculating the number of zero crossings is:

[0021]

[0022] l is the total number of sampling points, x j Represents the jth sampling point of the sample.

[0023] Step 4: Extract frequency domain features;

[0024] Extracting the frequency domain features of the EEG signal of the channel selected in step 2 specifically includes the following steps:

[0025] 4.1. For each selected channel, x modes are obtained using the EWT algorithm. The specific steps are as follows:

[0026] (1) Perform fast Fourier transform to obtain the Fourier spectrum in the range of 0 to π.

[0027] (2) Using empirical rules to determine a threshold T, for any local minimum point of the Fourier spectrum, as long as the scale space curve O generated by it i The length satisfies L i >T, then the local minimum point corresponds to the frequency ω i becomes the boundary frequency of its Fourier spectrum segmentation; after finding all the boundary frequencies of the Fourier spectrum, the Fourier spectrum is segmented into continuous N continuous segments.

[0028] (3) After determining the boundary frequency in the Fourier spectrum, the empirical wavelet function is constructed according to the construction ideas of Littlewood–Paley and Meyer wavelets. The empirical wavelet function is the bandpass filter on each continuous segment, and the expression is:

[0029]

[0030] ω n is the boundary frequency of the n-1th and nth filters, γ is the transition coefficient, which should satisfy β(x) is the indicative function, usually defined as:

[0031]

[0032] (4) Similar to the traditional wavelet transform, the detail coefficient W of the empirical wavelet function f (n,t) can be calculated as:

[0033]

[0034] Decompose to get the empirical mode f n (t) is:

[0035]

[0036] Where f is the decomposed function, t is time, τ is the integral variable, * is the convolution operation, ∧ is the complex conjugate operation, ∨ is the inverse Fourier transform, and - is the conjugate operation.

[0037] 4.2. Calculate the spectrum density of x modes by Welch power spectrum density analysis method, and obtain the mode with the largest spectrum density value as the single mode selected for each channel.

[0038] 4.3. For each channel single mode, the Hilbert transform method is used to extract the instantaneous amplitude and instantaneous frequency components of the signal as the frequency domain features of the signal. The frequency domain training set matrix f is obtained 2train , and the frequency domain test set matrix f 2test .

[0039] Step 5: Extract spatial features;

[0040] An improved CSP algorithm, namely the weighted regularized common spatial pattern algorithm, is used to extract the spatial features of the EEG information of the channel selected in step 2. Specifically, the following steps are included:

[0041] 5.1. Obtaining Regularized Co-spatial Pattern Algorithm

[0042] In order to overcome the sensitivity of the CSP algorithm to noise and overfitting, the Tikhonov regularization term is added to the objective function to obtain a regularized co-space mode algorithm to penalize solutions that do not satisfy the prior information. The objective function is expressed as:

[0043]

[0044] Where C1 and C2 are the covariance matrices of category 1 and category 2 respectively, and P(w) represents the penalty function, which is used to measure the degree to which the spatial filter w satisfies the given prior conditions. The more the filter w satisfies the conditions, the smaller the corresponding penalty function value P(w). Therefore, in order to maximize It is necessary to minimize P(w) to ensure that the spatial filter satisfies the prior information. α is the defined regularization parameter. The larger its value, the more it can meet the conditions. Where P(w) is a quadratic penalty function, expressed as:

[0045]

[0046] Where I represents the identity matrix. Adding the quadratic penalty term, formula (9) can be rewritten as:

[0047]

[0048] The corresponding Lagrangian function expression is:

[0049]

[0050] Using the Lagrangian operator λ to derive w and make its derivative equal to 0, the eigenvalue expression is shown in formula (13):

[0051] (C2+αI) -1 C1w=λw (13)

[0052] By maximizing the objective function in formula (14), we obtain a filter that maximizes C2 while minimizing C1:

[0053]

[0054] Get the eigenvector corresponding to the largest eigenvalue M2 = (C2 + αI) -1 C2, as the final filter w.

[0055] 5.2. Obtaining the Weighted Tikhonov Regularized Co-spatial Pattern Algorithm

[0056] Since the Tikhonov regularized common spatial pattern algorithm penalizes all high-weight channels equally, there is no specific distinction between the contribution of channels. However, for epileptic EEG signals, the features extracted from different channels contribute differently to the spatial filter. Therefore, different penalty items need to be set for different channels. That is, if a given channel contributes greatly to the spatial filter, then this channel should be given a relatively small penalty, that is, to ensure that the CSP algorithm assigns it a large weight. On the contrary, a larger penalty is given to this channel.

[0057] To this end, a weighted Tikhonov regularized CSP algorithm is proposed. The specific expression of its penalty term is:

[0058] P d (w) = w T D w w (15)

[0059] Where D w is a diagonal matrix, which can be expressed as:

[0060] D w =diag(w G ) (16)

[0061] In the formula, w G (i) is the penalty level assigned to the i-th channel. G The calculation formula for (i) is:

[0062]

[0063] In the formula, It represents the fth spatial filter obtained by using the CSP algorithm for the i-th sample, Ω is the total number of samples, 2*N f is the total number of filters. The penalty level of a channel is set to the inverse of the average weight of the channel in the CSP filter obtained from other samples. That is, the contribution of a given channel will be reflected by its weight in the filter obtained by the CSP algorithm from other samples. By taking different values ​​of the number of spatial filters m, different classification accuracies are obtained. The value m that maximizes the classification accuracy is used as the final number of filters to obtain the spatial domain training set matrix f 3train , and the test set matrix f 3test .

[0064] Step 6: Merge the time domain features, frequency domain features and spatial domain feature matrices obtained in steps 3, 4 and 5 to obtain the training feature matrix f train =[f 1train f 2train f 3train ] and the test feature matrix f test =[f 1test f 2test f 3test ].

[0065] Step 7: Input the feature matrix obtained in step 6 into the support vector machine-based classifier according to the training set and the test set to obtain the recognition result.

[0066] The classification process is as follows: after obtaining the training set and test set data, the support vector machine is used to establish a training model. The LIBSVM toolbox is used to select the radial basis kernel function to train the training set matrix and classification labels. The cross-validation method is used to globally optimize the kernel function parameters and error penalty factors, and the optimal parameters are input into the support vector machine for actual classification of the test set.

[0067] Preferably, the time of normal data in step 1 is 60 minutes, and the non-overlapping sliding window size a used is 10 seconds.

[0068] Preferably, the threshold b set in step 2 is 0.7.

[0069] Compared with the prior art, the present invention has the following beneficial effects:

[0070] (1) Select the EEG channels in the sample and only select the channels without repeated information, which reduces the amount of data and improves the running speed.

[0071] (2) The EWT method is used to overcome the limitations of the EMD method in the process of signal decomposition and analysis. By using a filter bank composed of adaptive wavelet filters, the EWT method is able to decompose any non-stationary signal into different modes, each of which is adjusted around a clear frequency to fully utilize the characteristics of the intrinsic mode function.

[0072] (3) A weighted regularized common spatial pattern algorithm is proposed to extract the spatial features of epileptic EEG signals. Compared with the traditional CSP algorithm, the weighted regularized common spatial pattern algorithm has a higher recognition rate and shorter calculation time.

[0073] (4) The characteristics of the EEG signal in the time domain, frequency domain and spatial domain are extracted simultaneously to obtain more comprehensive information.

[0074] (5) The time domain feature extraction method has low computational complexity, is easy to obtain, and can obtain better classification features. BRIEF DESCRIPTION OF THE DRAWINGS

[0075] Figure 1 It is a flow chart of the implementation of the present invention;

[0076] Figure 2 The time division of the interictal period and preictal period of the epileptic EEG signal in the embodiment;

[0077] Figure 3 A list of related channels of each channel obtained by sorting out in the embodiment; DETAILED DESCRIPTION

[0078] The embodiments of the present invention are described in detail below in conjunction with the accompanying drawings: Figure 1 As shown, the present invention comprises the following steps:

[0079] Step 1: Data preprocessing.

[0080] The epilepsy EEG data used in this example are selected from the CHB-MIT public database, and the first patient's data is used as the experimental data. The data channel is 23, and the data is divided into normal data and onset data. The 60 minutes before the onset of the onset are selected as normal data, and the normal data is further divided into interictal period and preictal period. Figure 2 As shown in the figure, the 30 minutes before the onset of the disease are selected as the preictal period, and the 30 to 60 minutes before the onset of the disease are selected as the interictal period. The normal data of the patient before the first three onsets are taken as the training set, and the normal data before the last two onsets are taken as the test set. The normal data before the onset are segmented using a non-overlapping sliding window with a window size of 10 seconds, and a window of 10 seconds of data is used as a test sample. The training set after segmentation has a total of 1080 samples, the first 540 of which are interictal samples, marked as category 1, and the last 540 are preictal samples, marked as category 2; the test set has a total of 720 samples, the first 360 of which are interictal samples, and the last 360 are preictal samples. The sampling frequency of the data is 256Hz, and the number of sample points for each sample is 2560.

[0081] Step 2: Select the EEG channel.

[0082] A sample is randomly selected from step 1, and then the EEG data in the sample is calculated using the correlation calculation method, and the EEG channels are selected based on the calculation results. If the Pearson correlation coefficient between different channels is high, it means that there is information duplication between the two channels.

[0083] The Pearson correlation coefficient is:

[0084]

[0085] Where: cov(X,Y) is the covariance between two variables, σ X σ Y is the standard deviation between the two variables.

[0086] According to formula (18), the correlation coefficients between two channels are calculated to obtain the correlation matrix, in which the correlation coefficients on the left diagonal of the matrix are all 1, indicating that the correlation between a channel and itself is completely positively correlated. The threshold is set to 0.7, and a correlation coefficient greater than 0.7 indicates that they are correlated; a correlation coefficient less than 0.7 indicates that they are not correlated. After sorting, the correlation channels of each channel can be obtained, such as Figure 3 The selection process of EEG channels is as follows:

[0087] (1) Select all channels that do not have related channels;

[0088] The selected channels are channels 4, 6, 7, 8, 10, 11, 17, 18, and 20;

[0089] (2) For channels that are equivalent to each other, the channel with the smaller channel number is selected;

[0090] The selected channels are channels 3, 5, 12, 14, and 15;

[0091] (3) Select the channel that contains the most channels from the remaining channels;

[0092] The remaining channels include 3 related channels, namely, channels 2, 13, and 22, and channels 13 and 22 are both related to channel 2, so channel 2 is selected; the remaining channels include 2 related channels, namely, channels 1 and 21, and channels 1 is related to the selected channel 5, and channels 21 is related to the selected channel 2, so neither of them is selected;

[0093] According to the calculation results of the Pearson correlation coefficient, 15 channels were obtained after sifting out the 23 channels, namely channels 2, 3, 4, 5, 6, 7, 8, 10, 11, 12, 14, 15, 17, 18, and 20.

[0094] Step 3: Extract time domain features;

[0095] The time domain features of the 15-channel epilepsy EEG data obtained in step 2 are extracted. The three time domain features are: root mean square value x RMS , absolute mean value x MAV and the number of zero crossings x ZC .

[0096] RMS value x RMS It reflects the power of epileptic EEG signal and can measure the energy of the signal:

[0097]

[0098] Absolute mean x MAV It reflects the average amplitude of epileptic EEG signal deviating from zero moment:

[0099]

[0100] Zero crossing number x ZC Reflects the degree of fluctuation of epileptic EEG signals:

[0101]

[0102] l is the total number of sampling points, x j Represents the jth sampling point of the sample.

[0103] Get the time domain training set matrix f with a size of 1080×(15×3) 1train , and the time domain test set matrix f of size 720×(15×3) 1test .

[0104] Step 4: Extract frequency domain features;

[0105] The frequency domain feature extraction is performed on the 15-channel epileptic EEG data obtained in step 2. The following steps are included:

[0106] 4.1. Use the EWT algorithm to obtain x modes for each of the 15 channels. The specific steps are as follows:

[0107] (1) Perform fast Fourier transform to obtain the Fourier spectrum in the range of 0 to π.

[0108] (2) Using empirical rules to determine a threshold T, for any local minimum point of the Fourier spectrum, as long as the scale space curve O generated by it i The length satisfies L i >T, then the local minimum point corresponds to the frequency ω i becomes the boundary frequency of its Fourier spectrum segmentation; after finding all the boundary frequencies of the Fourier spectrum, the Fourier spectrum is segmented into continuous N continuous segments.

[0109] (3) After determining the boundary frequency in the Fourier spectrum, the empirical wavelet function is constructed according to the construction ideas of Littlewood–Paley and Meyer wavelets. The empirical wavelet function is the bandpass filter on each continuous segment, and the expression is:

[0110]

[0111] In the formula, ω n is the boundary frequency of the n-1th and nth filters, γ is the transition coefficient, which should satisfy β(x) is the indicative function, usually defined as:

[0112]

[0113] (4) Similar to the traditional wavelet transform, the detail coefficient W of the empirical wavelet function f (n,t) can be calculated as:

[0114]

[0115] Decompose to get the empirical mode f n (t) is:

[0116]

[0117] Where f is the decomposed function, t is time, τ is the integral variable, * is the convolution operation, ∧ is the complex conjugate operation, ∨ is the inverse Fourier transform, and - is the conjugate operation.

[0118] 4.2. Calculate the spectrum density of x modes by Welch power spectrum density analysis method, and obtain the mode with the largest spectrum density value as the single mode selected for each channel.

[0119] 4.3. For each single mode of the channel, the Hilbert transform method is used to extract the instantaneous amplitude and instantaneous frequency components of the signal as the frequency domain features of the signal.

[0120] The frequency domain training set matrix f with a size of 1080×(15×2) is obtained 2train , and the frequency domain test set matrix f of size 720×(15×2) 2test .

[0121] Step 5: Extract spatial features;

[0122] An improved CSP algorithm, namely the weighted regularized common spatial pattern algorithm, is used to extract the spatial features of the 15-channel epileptic EEG data obtained in step 2.

[0123] 5.1. Obtaining Regularized Co-spatial Pattern Algorithm

[0124] In order to overcome the sensitivity of the CSP algorithm to noise and overfitting, the Tikhonov regularization term is added to the objective function to obtain a regularized co-space mode algorithm to penalize solutions that do not satisfy the prior information. The objective function is expressed as:

[0125]

[0126] Where C1 and C2 are the covariance matrices of category 1 and category 2 respectively, and P(w) represents the penalty function, which is used to measure the degree to which the spatial filter w satisfies the given prior conditions. The more the filter w satisfies the conditions, the smaller the corresponding penalty function value P(w). Therefore, in order to maximize It is necessary to minimize P(w) to ensure that the spatial filter satisfies the prior information. α is a defined regularization parameter. The larger its value, the more it can meet the conditions. Using the regularization parameter, the optimization process can be directed to a good spatial filter when the training data is limited or noisy.

[0127] Where P(w) is a quadratic penalty function, expressed as:

[0128]

[0129] In the formula, I represents the unit matrix.

[0130] By adding the quadratic penalty term, formula (26) can be rewritten as:

[0131]

[0132] The corresponding Lagrangian function expression is:

[0133]

[0134] Using the Lagrangian operator λ to derive w and make its derivative equal to 0, the eigenvalue expression is shown in formula (30):

[0135] (C2+αI) -1 C1w=λw (30)

[0136] By maximizing the objective function in formula (14), we obtain a filter that maximizes C2 while minimizing C1:

[0137]

[0138] Get the eigenvector corresponding to the largest eigenvalue M2 = (C2 + αI) -1 C2, as the final filter w.

[0139] 5.2. Obtaining the Weighted Tikhonov Regularized Co-spatial Pattern Algorithm

[0140] Since the Tikhonov regularized common spatial pattern algorithm penalizes all high-weight channels equally, there is no specific distinction between the contribution of channels. However, for epileptic EEG data, the features extracted from different channels contribute differently to the spatial filter. Therefore, different penalty items need to be set for different channels. That is, if a given channel contributes greatly to the spatial filter, then this channel should be given a relatively small penalty, that is, to ensure that the CSP algorithm assigns it a large weight. On the contrary, a larger penalty is given to this channel.

[0141] To this end, a weighted Tikhonov regularized CSP algorithm is proposed. The specific expression of its penalty term is:

[0142] P d (w) = w T D w w (32)

[0143] Where D w is a diagonal matrix, which can be expressed as:

[0144] D w =diag(w G ) (33)

[0145] In the formula, w G (i) is the penalty level assigned to the i-th channel. G The calculation formula for (i) is:

[0146]

[0147] In the formula, It represents the fth spatial filter obtained by using the CSP algorithm for the i-th sample, Ω is the total number of samples, 2*N f is the total number of filters. That is, the penalty level of a channel is set to the inverse of the average weight of the channel in the CSP filter obtained from other samples. Therefore, the contribution of a given channel will be reflected by its weight in the filter obtained by the CSP algorithm from other samples. In the spatial feature extraction based on improved CSP, the number of spatial filters m directly affects the experimental results. After a large number of repeated experiments, different classification accuracies are obtained by taking different m values. The m value that maximizes the classification accuracy is used as the final number of filters, denoted as M. Therefore, the spatial training set matrix f with a size of 1080×M is obtained. 3train , and the spatial domain test set matrix f of size 720×M 3test .

[0148] Step 6: Merge the time domain features, frequency domain features and spatial domain feature matrices obtained in steps 3, 4 and 5 to obtain the training feature matrix f train =[f 1train f 2train f 3train ] and the test feature matrix f test =[f 1test f 2test f 3test ].

[0149] Step 7: Input the feature matrix obtained in step 6 into a classifier based on a support vector machine (SVM) according to the training set and the test set to obtain the recognition result.

[0150] Support vector machine classification process: After obtaining the training set and test set data, use SVM to establish a training model. SVM is trained on the training set matrix and classification labels, using the LIBSVM toolbox and the radial basis kernel function. The cross-validation method is used to globally optimize the kernel function parameters and error penalty factors, and the optimal parameters are input into the support vector machine for actual classification of the test set.

Claims

1. An epilepsy signal classification method based on EWT and improved CSP, characterized by: The method comprises the following steps: Step 1: Data preprocessing; The selected EEG data are divided into normal data and onset data based on the onset time; the data 30 to 120 minutes before the onset time are taken as normal data; the first half of the normal data are taken as interictal data, and the second half as preictal data; after selecting the training set and the test set, the normal data in the set are segmented using a non-overlapping sliding window with a window size of a seconds, and the data in one window is taken as an experimental sample; the samples in the interictal period in the training set are marked as category 1, and the samples in the preictal period are marked as category 2; Step 2: Select the EEG channel; For the processed sample data, the EEG signals are calculated using the correlation calculation method, and the EEG channels in the sample are selected according to the calculation results; the calculation formula of the Pearson correlation system is: cov(X,Y) is the covariance between two variables, σ X σ Y is the standard deviation between the two variables; The correlation coefficient between two channels is obtained according to the calculation result of formula (1), and thus the correlation matrix is ​​obtained; a threshold b is set, and the correlation coefficient between two channels is greater than the threshold, indicating that they are correlated, otherwise it indicates that they are uncorrelated; the EEG channel is selected according to the correlation coefficient, and the selection process is as follows: (1) Select all channels that do not have related channels; (2) For channels that are equivalent to each other, the channel with the smaller channel number is selected; (3) Select the channel that contains the most channels from the remaining channels; Step 3: Extract time domain features; Extract the root mean square value, absolute mean value and number of zero crossings of the EEG signal of the channel selected in step 2 as the time domain features of the signal, and obtain the time domain training set matrix f 1train and the time domain test set matrix f 1test ; Step 4: Extract frequency domain features; Extracting the frequency domain features of the EEG signal of the channel selected in step 2 specifically includes the following steps: 4.

1. For each selected channel, x modes are obtained using the EWT algorithm. The specific steps are as follows: (1) Perform fast Fourier transform to obtain the Fourier spectrum in the range of 0 to π; (2) Using empirical rules to determine a threshold T, for any local minimum point of the Fourier spectrum, as long as the scale space curve O generated by it i The length satisfies L i >T, then the local minimum point corresponds to the frequency ω i becomes the boundary frequency of its Fourier spectrum segmentation; after finding all the boundary frequencies of the Fourier spectrum, the Fourier spectrum is segmented into N consecutive segments; (3) After determining the boundary frequency in the Fourier spectrum, the empirical wavelet function is constructed according to the construction ideas of Littlewood-Paley and Meyer wavelets. The empirical wavelet function is the bandpass filter on each continuous segment, and the expression is: In the formula, ω n is the boundary frequency of the n-1th and nth filters, γ is the transition band coefficient, which should satisfy β(x) is the indicative function, usually defined as: (4) Similar to the traditional wavelet transform, the detail coefficient W of the empirical wavelet function f (n,t) can be calculated as: Decompose to get the empirical mode f n (t) is: In the formula, f is the decomposed function, t is time, τ is the integral variable, * is the convolution operation, ∧ is the complex conjugate operation, ∨ is the inverse Fourier transform, is the conjugate operation; 4.

2. Calculate the spectrum density of x modes by Welch power spectrum density analysis method, and obtain the mode with the largest spectrum density value as the single mode selected for each channel; 4.

3. For each channel, the Hilbert transform method is used to extract the instantaneous amplitude and instantaneous frequency components of the signal as the frequency domain features of the signal; the frequency domain training set matrix f is obtained 2train , and the frequency domain test set matrix f 2test ; Step 5: Extract spatial features; An improved CSP algorithm, namely the weighted regularized common spatial pattern algorithm, is used to extract the spatial domain features of the EEG information of the channel selected in step 2; specifically, the following steps are included: 5.

1. Obtaining Regularized Co-spatial Pattern Algorithm In order to overcome the sensitivity of the CSP algorithm to noise and overfitting, the Tikhonov regularization term is added to the objective function to obtain a regularized co-space mode algorithm to penalize solutions that do not satisfy the prior information. The objective function is expressed as: Where C1 and C2 are the covariance matrices of category 1 and category 2 respectively, and P(w) represents the penalty function, which is used to measure the degree to which the spatial filter w satisfies the given prior conditions; the more the filter w satisfies the conditions, the smaller the corresponding penalty function value P(w); therefore, in order to maximize It is necessary to minimize P(w) to ensure that the spatial filter satisfies the prior information; α is the defined regularization parameter, the larger its value, the more it can meet the conditions; Where P(w) is a quadratic penalty function, expressed as: In the formula, I represents the unit matrix; By adding the quadratic penalty term, formula (6) can be rewritten as: The corresponding Lagrangian function expression is obtained as follows: Using the Lagrangian operator λ to derive w and make its derivative equal to 0, the eigenvalue expression is shown in formula (10): (C2+αI) -1 C1w=λw (10) By maximizing the objective function in formula (11), we obtain a filter that maximizes C2 while minimizing C1: Get the eigenvector corresponding to the largest eigenvalue M2 = (C2 + αI) -1 C2, as the final filter w; 5.

2. Obtaining the Weighted Tikhonov Regularized Co-spatial Pattern Algorithm Since the Tikhonov regularized common spatial pattern algorithm imposes equal penalties on all high-weight channels, there is no specific distinction between the contribution of channels; however, for epileptic EEG signals, the features extracted from different channels contribute differently to the spatial filter; therefore, different penalty items need to be set for different channels; To this end, a weighted Tikhonov regularized CSP algorithm is proposed; the specific expression of its penalty term is: P d (w)=w T D w w (12) Where D w is a diagonal matrix, which can be expressed as: D w =diag(w G ) (13) In the formula, w G (i) is the penalty level assigned to the i-th channel; w G The calculation formula for (i) is: In the formula, It represents the fth spatial filter obtained by using the CSP algorithm for the i-th sample, Ω is the total number of samples, 2*N f is the total number of filters; the penalty level of a channel is set to the inverse of the average weight of the channel in the CSP filter obtained from other samples; that is, the contribution of a given channel will be reflected by its weight in the filter obtained by the CSP algorithm from other samples; by taking different values ​​of the number of spatial filters m, different classification accuracies are obtained, and the value m that maximizes the classification accuracy is used as the final number of filters to obtain the spatial domain training set matrix f 3train , and the test set matrix f 3test ; Step 6: Merge the time domain features, frequency domain features and spatial domain feature matrices obtained in steps 3, 4 and 5 to obtain the training feature matrix f train =[f 1train f 2train f 3train ] and the test feature matrix f test =[f 1test f 2test f 3test ]; Step 7: Input the feature matrix obtained in step 6 into the support vector machine-based classifier according to the training set and the test set to obtain the recognition result.

2. The epilepsy signal classification method based on EWT and improved CSP as claimed in claim 1, characterized in that: In step 1, the normal data time is 60 minutes, and the non-overlapping sliding window size a used is 10 seconds.

3. The epilepsy signal classification method based on EWT and improved CSP as claimed in claim 1, characterized in that: In the step 2, the threshold b is set to 0.

7.

4. The epilepsy signal classification method based on EWT and improved CSP as claimed in claim 1, characterized in that: In step 3, the calculation formula of the root mean square value is: The absolute mean value is calculated as: The formula for calculating the number of zero crossings is: In the formula, l is the total number of sampling points, x j Represents the jth sampling point of the sample.

5. The epilepsy signal classification method based on EWT and improved CSP as claimed in claim 1, characterized in that: In step seven, the support vector machine classification process is as follows: after obtaining the training set and test set data, a training model is established using the support vector machine; the LIBSVM toolbox is used to select the radial basis kernel function to train the training set matrix and classification labels; the cross-validation method is used to globally optimize the kernel function parameters and the error penalty factor, and the optimal parameters are input into the support vector machine for actual classification of the test set.

Images