Motor imagery electroencephalogram signal classification method combining function type data analysis and kernel extreme learning machine

Abstract

The invention discloses a motor imagery electroencephalogram signal classification method combining function type data analysis and a kernel extreme learning machine. The method is an MI-EEG signal four-classification problem identification method combining function type data analysis and a sparrow search algorithm optimization kernel extreme learning machine (SSAKELM). Firstly, collected MI-EEG signals are preprocessed, data samples are established according to the duration of single motor imagery, and then a multi-channel and multi-sample alpha-waveband power spectral density (PSD) and frequency relation curve is obtained through FIR filtering. According to the basic theory of functional data analysis (FDA), functional smoothing processing is carried out on all PSD-frequency curves one by one according to brain electrode channels by utilizing a B spline, then characteristic parameters of each sample are obtained by utilizing functional principal component analysis (FPCA) and principal differential analysis (PDA), and then the characteristic parameters of all the samples are merged for normalization processing. The method comprises the following steps: selecting a training sample proportion, randomly selecting a training sample, carrying out clustering training by using SSAHELM, and then carrying out category prediction on a test sample. The MI-EEG signal multivariate characteristic parameter extraction method provided by the invention has clear mathematical and physical significance, mainly relates to matrix operation, is small in calculation scale, and is easy to embed into a chip.

Description

Technical Field

[0001] The present invention relates to the technical field of motor imagery electroencephalogram (MI-EEG) classification problems in brain-computer interface technology, and specifically, to an accurate identification method for the four-classification problem of MI-EEG signals by combining functional principal component analysis (FPCA) and principal differential analysis (PDA) feature extraction with a sparrow search algorithm optimized kernel extreme learning machine (SSA_KELM). Background Art

[0002] A brain-computer interface (BCI) is a connection between the human brain and a computer (such as Figure 2 As shown in the figure), it realizes the information exchange between the brain and the device. This technology has profound theoretical research significance and broad application prospects. At present, my country has made plans and ideas for brain-computer interfaces at the strategic level.

[0003] Motor imagery electroencephalogram (MI-EEG) is an important research direction in the field of BCI research. It is a non-invasive brain-computer interface technology that records brain activity from the surface of the head through a scalp wearable device. Motor imagery brain-computer interface technology collects EEG signals, performs preprocessing, feature extraction, and pattern recognition, completes signal processing, deciphers the brain state or intention, encodes the processed information and converts it into specific instructions, so as to control the application of external devices, realize information exchange, movement assistance, and function recovery, etc., and may eventually transmit neural feedback signals to the human brain to complete the closed-loop processing of the communication system. Brain-computer interface technology based on motor imagery has broad application prospects in medical rehabilitation (such as intelligent prostheses and wheelchairs), auxiliary tool development, brain-controlled weapons (such as brain-controlled drones, brain-controlled robots, etc.), and other fields.

[0004] Motor imagery EEG signals have the following characteristics: (1) The amplitude of the EEG signal is weak, usually only a few microvolts to about 100 microvolts.

[0005] (2) The signal-to-noise ratio of EEG signals is low. EEG signals are prone to aliasing and interference between signals. At the same time, they are also interfered and attenuated when passing through the scalp and skull during acquisition. They are also easily interfered by factors such as the environment and body movement, such as electromagnetic interference from electrodes, cables, and power supplies, as well as physiological interference from muscle movement, blinking, and breathing. Using acquisition equipment with higher quality electrodes and cables, keeping the subjects in a relatively calm state and quiet environment, and using some advanced preprocessing and analysis methods can effectively improve the signal-to-noise ratio.

[0006] (3) EEG signals are highly non-stationary random signals.

[0007] (4) EEG signals have nonlinear characteristics, and their changes in time and space are nonlinear. The phase difference between different frequency components of electrical signals between different neurons changes irregularly over time, and nonlinear oscillations will also occur; the complex coupling relationship between neurons will show complex dynamic phenomena such as chaos, bifurcation and period doubling, which are all manifestations of the nonlinear characteristics of EEG signals.

[0008] (5) EEG signals are broadband multi-component signals with multiple frequency components in the frequency domain. The energy distribution of these frequency components is relatively wide and is not concentrated on a certain frequency point. There is a nonlinear coupling relationship between different frequency components, and there are artifacts and cross-interference in the signal.

[0009] (6) Multi-channel acquisition provides more spatial domain information for EEG signals. Multi-channel acquisition can remove common noise and cross-interference, provide higher signal resolution, and accurately locate the source of EEG activity. Conversely, more accurate signals can be acquired based on their spatial positioning.

[0010] The classification of motor imagery EEG signals is a complex and important research field, which involves multiple aspects such as the acquisition, processing, analysis and classification of EEG signals.

[0011] Brain waves are divided into five categories according to their frequency: alpha waves (8-14Hz), beta waves (14-30Hz), theta waves (4-8Hz), delta waves (below 4Hz) and gamma waves (above 30Hz).

[0012] The practice of EEG signal testing shows that the EEG signals of movement phenomena are mainly concentrated in the α band, so signal analysis in the α band is of great significance for MI-EEG signal classification.

[0013] Due to the low signal-to-noise ratio of the original MI-EEG signal and the high-dimensional characteristics of the data, it is almost impossible to directly use the original EEG data for classification. It must go through a series of preprocessing and feature parameter extraction processes to achieve a significant dimensionality reduction of the original data.

[0014] Preprocessing of motor imagery EEG signals is an important part of brain-computer interface technology, and it is crucial to improving the accuracy and efficiency of subsequent signal processing. The content of preprocessing mainly includes: (1) Filtering. Filtering is one of the most basic and important steps in the preprocessing process. The purpose of filtering is to remove noise and artifacts in EEG signals and improve the signal-to-noise ratio (SNR) of the signal. Filtering can be performed in multiple dimensions such as time domain, spatial domain and frequency domain. Common filtering methods include: time domain filtering, frequency domain filtering, and spatial domain filtering. (2) Artifact removal. EEG signals often contain artifacts that are not related to brain activity, such as electrooculography, electromyography, and electrocardiography. These artifacts will interfere with the analysis and classification of EEG signals, so they need to be removed during the preprocessing process. Artifact removal methods include independent component analysis (ICA), principal component analysis (PCA), etc., which can decompose the signal into multiple independent components and remove components related to artifacts. (3) Signal segmentation. Since EEG signals are continuous time series data, in order to facilitate subsequent processing and analysis, it is usually necessary to segment the signal into multiple shorter segments (also called epochs). Signal segmentation needs to determine the length and overlap of the segments according to the specific requirements of the experimental design. (4) Baseline correction. When collecting EEG signals, DC offset may be generated due to factors such as the contact resistance between the electrode and the scalp. This offset will affect the subsequent processing and analysis of the signal, so baseline correction is required. Baseline correction usually eliminates DC offset by setting a certain time point or time period of the signal as the reference point and adjusting the signal values ​​at other time points. (5) Other preprocessing steps. In addition to the above steps, it may also include signal quality assessment, detection and processing of bad channels (damaged electrode channels), etc.

[0015] At present, the traditional methods for extracting basic MI-EEG signal features mostly use time domain, frequency domain, time-frequency domain and spatial domain analysis methods. Currently, there are methods such as Common Spatial Pattern (CSP), Principal Component Analysis (PCA), Common Average Referencing (CAR), adaptive filtering, Independent Component Analysis (ICA) and digital filtering.

[0016] Feature extraction through frequency domain analysis often uses Fast Fourier Transform (FFT) and Auto Regressive (AR). Since the EEG signal is divided into five frequency bands, FFT can be used to convert the EEG signal from the time domain to the frequency domain to extract the signal's features. The AR model simulates the EEG signal by representing the signal on each channel as a linear combination of the signal at the previous point. In addition, the power spectral density (PSD) of the EEG signal can be calculated to characterize its features. Wavelet Transform (WT), Wavelet Packet Decomposition (Wavelet Packet Decomposition), Ensemble Empirical Mode Decomposition (EEMD) that can be applied to any time series, and Hilbert-huang Transform (HHT) using Hilbert spectrum analysis are also commonly used methods in time-frequency domain analysis.

[0017] Although the above-mentioned feature extraction methods based on time domain, frequency domain, time-frequency domain and space domain have clear mathematical and physical meanings, their classification accuracy is very low.

[0018] The above-mentioned traditional EEG signal feature extraction methods cannot extract more effective EEG features from high-dimensional EEG data. The ability of deep learning to automatically learn high-dimensional data and select features at the same time solves this problem. Through the powerful analysis ability of deep learning, effective representation of EEG signals can be achieved, providing new possibilities for accurately identifying different brain signals.

[0019] Deep learning models have the ability to extract features from high-dimensional data and learn from hierarchical representations, so MI-EEG feature extraction and classification methods based on deep learning have become a new research trend. Convolutional Neural Network (CNN), Deep Belief Network (DBN), Recurrent Neural Network (RNN), and EEGnet are the current mainstream algorithms for analyzing EEG signals in deep learning. The fusion of different deep learning models and their joint application in EEG classification and recognition has become a new research trend. Many researchers use a hybrid network that fuses CNN and RNN. This hybrid network is not a simple series relationship, but a hybrid structure constructed in parallel. The application of deep learning in the field of EEG signal classification and recognition will develop in the direction of multi-scale analysis of EEG signals in time, frequency and space.

[0020] Compared with traditional MI-EEG signal classification methods, various neural network feature extraction classification methods have achieved high accuracy. However, these methods have poor interpretability and lack clear mathematical-physical meaning. In addition, if these algorithms are written into chips, the price will be very expensive.

[0021] In summary, there are many MI-EEG signal classification algorithms, but a simple and easy-to-use algorithm with universal adaptability has not yet been established. In addition, more research on classification algorithms is needed in chip design for classification problems, which requires considering the timeliness and ease of implementation of the algorithm. Summary of the invention

[0022] The purpose of the present invention is to provide a MI-EEG signal classification problem identification method combining functional data analysis (FDA) feature extraction and sparrow search algorithm optimized kernel extreme learning machine (SSA_KELM) to achieve accurate prediction of motor imagery electroencephalogram (MI-EEG) categories.

[0023] The technical solution to achieve the purpose of the present invention is: a method for classifying motor imagery EEG signals combining functional data analysis and a kernel extreme learning machine, comprising the following steps:

[0024] 1. Use high-precision EEG signal acquisition equipment with more than 8 channels and select a sampling frequency of 500Hz to collect MI-EEG signals; pre-process the collected MI-EEG signals according to the EEG signal pre-processing technical requirements;

[0025] 2. Establish data samples according to the duration of a single motor imagery, and obtain the power spectral density PSD-frequency relationship curve of the alpha band of multiple channels and multiple samples through FIR filter;

[0026] 3. According to the basic theory of functional data analysis of FDA, all PSD-frequency relationship curves are functionally smoothed;

[0027] 4. Use functional principal component analysis (FPCA) to obtain the principal component score of each data sample, which is called the principal component score eigenvalue;

[0028] 5. Use functional principal differential analysis PDA to obtain the principal differential score of each sample, which is called the principal differential score eigenvalue;

[0029] 6. Combine the principal component score eigenvalue and principal differential score eigenvalue of each sample, and normalize the characteristic parameters of all samples;

[0030] 7. Randomly select training samples and test samples, and use the sparrow search algorithm to optimize the kernel extreme learning machine SSA_HELM for clustering training;

[0031] 8. Make category predictions for samples to be classified.

[0032] Furthermore, the preprocessing of the collected MI-EEG signal in step 1 specifically includes:

[0033] Signal filtering: bandpass filtering to retain the most valuable frequency range for motor imagery tasks, which is between 1 and 40 Hz; notch filtering, which uses notch filters to suppress noise at specific frequencies;

[0034] Artifact removal: by decomposing the independent components in the MI-EEG signal and removing the components related to eye movements and muscle activity;

[0035] Channel selection: Select the EEG channel that is most valuable for the motor imagery task, and select the MI-EEG signals collected by the movement-related electrodes C4, P4, Pz, Cz, P3, and C3;

[0036] Baseline correction: The MI-EE signal in the period before stimulus presentation was used as the baseline, and this baseline value was subtracted from the time series data in the period before and after stimulus presentation;

[0037] Signal amplification: amplify the MI-EEG signal.

[0038] Furthermore, the specific implementation of step 2 is as follows:

[0039] 2.1 Obtain the alpha wave time series in the range of 8-14 Hz of the preprocessed MI-EEG signal through the FIR filter;

[0040] 2.2 For the above time series, establish data samples according to the duration of a single motor imagery;

[0041] 2.3 For each sample time series, perform FFT operation to obtain the PSD-frequency relationship curve of each α-wave motor imagery EEG signal mentioned above.

[0042] Furthermore, the method for performing functional smoothing on all PSD-frequency relationship curves in step 3 is:

[0043] 3.1 According to the brain electrode channel, the discrete PSD-frequency data of each channel is fitted into a continuous function one by one; the sample data belonging to the C4, P4, Pz, Cz, P3, and C3 channels are fitted with functions respectively, that is, the four types of sample data of each channel are put together for function fitting;

[0044] MI-EEG signals are non-periodic, so we select B-spline basis functions suitable for non-periodic data, and take 6-8th order B-splines, so that the fitting curve has good differentiability and locality;

[0045] 3.2 After selecting the basis function, the coefficient vector c is estimated by minimizing the sum of squared errors using the weighted least squares criterion i ; The larger the number of solved time-based functions, the better the data fitting, but the worse the function smoothness, and the noise or abnormal data in the original data may be fitted into the function; in order to take into account both the data fitting and the smoothness of the function, the a-order roughness penalty function is introduced.

[0046] Furthermore, in step 4, the functional principal component analysis FPCA is used to obtain the principal component score of each sample, which is specifically implemented as follows:

[0047] 4.1 Covariance function estimation; Estimating the covariance function of functional data, which describes the correlation between data changes at different time points or variable values;

[0048] 4.2 Solving eigenvalues ​​and eigenfunctions: By solving the eigenvalue problem of the covariance function, a series of eigenvalues ​​and corresponding eigenfunctions are obtained; these eigenvalues ​​and eigenfunctions represent the main change patterns in the data and the importance of these patterns in the data respectively;

[0049] 4.3 Sorting and selection: sort the eigenvalues ​​from large to small, and select the eigenfunctions corresponding to the first M largest eigenvalues, where M is the number of dimensions after dimensionality reduction;

[0050] 4.4 Projection and dimensionality reduction; Projection is to project the original functional data onto a low-dimensional subspace defined by the selected M characteristic functions, which is achieved by calculating the inner product of the original data and each characteristic function. The projection is the eigenvalue used for future classification, called the principal component score eigenvalue; the result after projection is the dimensionality reduction data, which represents the main change pattern in the original data in a lower dimension. Furthermore, in step 5, the functional principal differential analysis PDA is used to obtain the principal differential score of each sample, which is specifically implemented as follows:

[0051] (1) Construction of linear differential operators and differential equations

[0052] The K-order homogeneous differential equation is expressed as:

[0053]

[0054] In the formula, is the coefficient parameter, which can be a variable function or a constant;

[0055] Define a linear differential operator:

[0056]

[0057] Where I is the identity operator; the homogeneous differential equation associated with it is:

[0058]

[0059] The primary task of PDA is to estimate the linear differential operator L, hoping that the linear differential operator L can eliminate the given N functions x(t) = [x 1 (t),x 2 (t),…,x N (t)] T , that is, we hope that the fitting residual Lxi(t) provided by the differential equation (17) is as close to 0 as possible; the residual norm sum of squares is defined as the fitting criterion using the least squares method:

[0060]

[0061] If the linear differential equation is non-homogeneous, that is, there is a forcing function δ(t) representing the external input or disturbance of the system, then the differential equation is expressed as:

[0062]

[0063] At this time, Lx(t) = δ(t), and the least squares fitting criterion of the linear differential operator L becomes:

[0064]

[0065] (2) Estimation of weight function

[0066] Based on the least squares fitting criterion of formula (6), the coefficient parameters of the linear differential operator L are estimated using the point-by-point minimization method. The fitting criterion of the point-by-point estimation is:

[0067]

[0068] In the formula, When t is fixed, the fitting criterion is just a simple least squares fitting criterion, which can be expressed in matrix form:

[0069]

[0070] In the formula, τ(t)=[D K x 1 (t),D K x 2 (t),…,D K x N (t)] T ; Z(t) is the point-by-point estimation design matrix, whose i-th row element is z i (t) = [-x i (t),-Dx i (t),…,-D K-1 x i (t),δ i (t)];

[0071] Keeping t constant, the coefficient vector can be obtained by minimizing point by point The least squares solution of is:

[0072]

[0073] Pointwise estimates of coefficient vectors The existence of depends on Z(t) T The determinant of Z(t) is non-zero for any t;

[0074] The estimate of the linear differential operator L is obtained as follows:

[0075]

[0076] RSQ is a commonly used tool in PDA to measure the goodness of fit of differential equations:

[0077]

[0078] The value range of RSQ is [0,1]. The closer its value is to 1, the better the linear differential operator L fits the function.

[0079] (3) Obtaining the principal differential score eigenvalue

[0080] Theorem: Let Λ be the Hilbert space, Y be the linear closed subspace of Λ, Y ⊥ is the orthogonal complement of Y. Any element belonging to Λ can be uniquely decomposed into the component belonging to Y and the component belonging to Y ⊥ The sum of the weights;

[0081] According to the theorem, the function The function space can be orthogonally decomposed into two parts, denoted as and in express The null space of for orthogonal complement of i (t) can be orthogonally decomposed into in It is called the structural component and satisfies It combines the common functional differential characteristics of functional data; It is called the residual component and satisfies It reflects the linear differential operator The differential characteristics of functions cannot be described;

[0082] When PDA is used as a feature extraction technique, the coefficient parameter of the differential operator L is considered to be a constant, and it is assumed that there is no forcing function; for a K-order non-degenerate linear differential equation, any solution can be expressed as a linear combination of K linearly independent complex exponential functions; construct the characteristic equation:

[0083]

[0084] The characteristic equation has K linearly independent complex roots The corresponding complex exponential function is:

[0085]

[0086] The space spanned by K linearly independent complex exponential functions as a linear differential operator The null space of:

[0087]

[0088] function A set of basis functions are formed in the solution space of the differential equation, which are called principal differential functions or differential characteristic functions. It is expressed as:

[0089]

[0090] In the formula, is the function x i The kth principal differential score (PDS) of (t), which is the function x extracted by PDA i (t) is the characteristic parameter of

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

[0092] (1) The PSD-frequency curve of the MI-EEG signal is functionalized to process the discrete signal into a continuous signal. Functional principal component analysis (FPCA) and principal differential analysis (PDA) are used to extract signal features. Its characteristic is that after the data is functionalized, it will be regarded as a whole, rather than just a sequential arrangement of individual observations, and its internal structure can be mathematically mined. In addition, the functionalized data can be subjected to multi-order differential analysis to mine the dynamic information of the signal frequency change rate, which is not considered by almost all current classification algorithms.

[0093] (2) Mathematically rigorous, with clear physical meaning, subjective cognition has little influence on feature extraction. Both FPCA and PDA have strict mathematical definitions, and the extracted features are the projections of the sample functions on the basis functions, which have clear mathematical and physical meanings.

[0094] (3) The classification process is simple and clear, the solution is fast, and it is easy to implement in the chip. BRIEF DESCRIPTION OF THE DRAWINGS

[0095] Figure 1 Flowchart of the motor imagery EEG signal classification method combining functional data analysis feature extraction and kernel extreme learning machine.

[0096] Figure 2 Schematic diagram of the brain-computer interface (BCI) system.

[0097] Figure 3 For high-quality EEG acquisition process.

[0098] Figure 4 This is an example of the original EEG signal.

[0099] Figure 5 This is an example of 8-channel EEG data after preprocessing.

[0100] Figure 6 Examples of FFT and power spectral density analysis.

[0101] Figure 7 It is a 6-channel right foot movement imagination time series after FIR8-15Hz filtering.

[0102] Figure 8 This is an example of the EEG time series data of the right foot C4 channel after FIR8-15Hz filtering.

[0103] Fig. 9 This is an example of PSD-frequency samples for the C4 channel.

[0104] Fig.10 Example of PSD-frequency sample function for C4 channel.

[0105] Fig.11 This is an example of a PSD-frequency sample of the Cz channel.

[0106] Fig.12 Example of Cz channel PSD-frequency sample functionalization.

[0107] Fig.13 is the 10th-order principal component function (PCF) of the C3 channel samples.

[0108] Fig.14 is the 10th-order principal component function (PCF) of the Cz channel samples.

[0109] Fig.15 It is the differential characteristic function of the P4 channel sample function.

[0110] Fig.16 The differential characteristic function of the Pz channel sample function.

[0111] Fig.17 Examples of unnormalized principal component score eigenvalues.

[0112] Fig.18 is an example of an unnormalized principal differential score eigenvalue.

[0113] Fig.19 This is an example of normalized training sample feature values ​​(when the training samples account for 70%, there are 336 training samples in total).

[0114] Fig. 20 This is an example of normalized eigenvalues ​​of test samples (when the training samples account for 70%, there are 144 test samples in total).

[0115] Fig.21 It is the KELM network structure.

[0116] Fig. 22 This is the SSA-KELM algorithm process.

[0117] Fig.23 This is an example of a confusion matrix.

[0118] Fig.24 This is a comparison chart between the true category and the predicted category (when the training accounts for 70%). DETAILED DESCRIPTION

[0119] In view of the fact that the current deep learning-based motor image electroencephalogram (MI-EEG) feature extraction methods lack interpretability and are difficult to implement in chips, the present invention proposes a multivariate feature parameter extraction method with clear mathematical and physical meanings, which is easy to implement in chips.

[0120] The present invention will be further described below in conjunction with the accompanying drawings and embodiments.

[0121] S1: High-quality MI-EEG signal acquisition.

[0122] For the four-category problem, you can imagine the movements of the left hand, right hand, left foot, and right foot respectively. In order to facilitate subsequent machine learning training, try to collect more samples.

[0123] When choosing test equipment, try to use professional equipment software, at least 8 acquisition channels: Use professional EEG acquisition equipment and analysis software to ensure the accuracy and consistency of data acquisition and processing. These devices usually have the necessary filtering, amplification and noise reduction functions to help improve signal quality.

[0124] When subjects are measuring motor imagery EEG signals, in order to ensure the accuracy and reliability of the data, they need to pay attention to the following matters:

[0125] (1) Stay relaxed. Before the experiment, the subjects should try to keep their bodies relaxed and avoid strenuous exercise or emotional fluctuations to reduce the interference of physiological factors on EEG signals. At the same time, the experimental environment should be quiet and comfortable, avoiding external noise and light influences.

[0126] (2) Reduce eye movements. During the experiment, the subjects should try to reduce blinking and eye movements as much as possible, because these movements will produce obvious artifacts and affect the quality of EEG signals. Eye fatigue can be relieved by taking a proper rest between experiments.

[0127] (3) Maintain focused imagination. The subjects need to enter the motor imagination state immediately after the prompt appears, and try to stay focused without being disturbed by the outside world. The imagined movements should be clear and coherent, such as imagining clenching a fist or moving toes, which helps to enhance the motor imagination component in the EEG signal.

[0128] (4) Maintain a stable posture. During the experiment, the subject should keep the body and head posture stable to avoid introducing interference signals due to changes in body position. A headrest or other support can be used to help fix the head position.

[0129] (5) Conduct multiple repeated experiments. Since EEG signals have large randomness and individual differences, it is recommended to conduct multiple repeated experiments to obtain more reliable and stable data results. An appropriate rest time can be set between each experiment to avoid fatigue.

[0130] (6) Pay attention to the electrode impedance. Before and after the experiment, make sure that all electrodes are in good contact with the scalp, and the electrode impedance should be controlled within a reasonable range (usually less than 20 kΩ). Conductive paste or saline can be used to improve the contact quality between the electrode and the skin.

[0131] (7) Record basic information. The experimental record should include the basic information of the subjects, such as age, gender, health status, etc. This information will help with subsequent data analysis and interpretation.

[0132] Figure 3 It shows the acquisition process of an advanced EEG device. Figure 4This is an example of an original EEG signal. In the example of this patent, the EEG device measures and obtains 6 channels of 4×120 MI-EEG sample data of four categories (left hand, right hand, left foot, right foot motor imagery), and each motor imagery lasts for 5 seconds.

[0133] S2: Preprocess the collected MI-EEG signals.

[0134] The preprocessing of motor imagery EEG signals is a key link in brain-computer interface technology, which is crucial for improving signal quality, reducing noise interference, and subsequent feature extraction and classification recognition. Motor imagery EEG signal preprocessing includes signal filtering, artifact removal, channel selection, etc.

[0135] (1) Signal filtering. The frequency range most valuable for motor imagery tasks is retained through a bandpass filter, usually between 0.5 and 50 Hz. Notch filtering may also be included, that is, noise of a specific frequency (such as power frequency noise) is suppressed using a notch filter. This helps to remove low-frequency and high-frequency noise and improve the signal-to-noise ratio of the signal.

[0136] (2) Artifact removal: Independent component analysis is the most commonly used method, which reduces data noise by decomposing the independent components in the EEG signal and removing components related to eye movements and muscle activities.

[0137] (3) Channel selection. Select the EEG channels that are most valuable for motor imagery tasks to reduce the dimensionality and complexity of data processing. Related studies often select electrode positions related to movement, such as C4, P4, Pz, Cz, P3, and C3.

[0138] (4) Baseline correction: The data in the period before stimulus presentation is used as the baseline, and this baseline value is subtracted from the subsequent data processing to eliminate differences between subjects.

[0139] (5) Signal amplification: Since EEG signals are weak bioelectric signals (generally in the range of 5-100 μV), they usually need to be amplified before feature extraction.

[0140] (6) Detrending: Remove the long-term trend term in the signal to eliminate low-frequency interference such as baseline drift.

[0141] (7) Spatial filtering: For example, applying a common spatial pattern (CSP) filter to enhance task-related EEG components and suppress irrelevant signals.

[0142] Figure 5 shows an example of 8-channel EEG data after preprocessing.

[0143] S3: Perform FFT analysis and power spectral density analysis on the preprocessed MI-EEG signal for channel selection.

[0144] Figure 6 This is an example of FFT and power spectral density analysis of MI-EEG signals, where obvious α-band features can be seen.

[0145] S4: Obtain the alpha wave of the preprocessed MI-EEG signal through the FIR filter.

[0146] Alpha waves play an important role in motor imagery EEG signals. They reflect the state of the brain, promote motor imagery, and affect the accuracy and reliability of subsequent processing. Their role is mainly reflected in the following aspects:

[0147] (1) Reflecting the state of the brain. Alpha waves usually appear in a state of relaxation and preparation, such as closing your eyes to rest, meditating, or just waking up from sleep. In the motor imagery task, when the subject enters a relaxed state and prepares for motor imagery, the alpha wave will be significantly enhanced. The appearance of this waveform indicates that the brain is in a relatively quiet but alert state, which is conducive to motor imagery. Alpha waves are one of the basic rhythms in adult EEG, and their frequency range is between 8-14Hz. It reflects the rhythmic activity of the cerebral cortex and is one of the important indicators of brain activity.

[0148] (2) Promote motor imagery. The appearance of alpha waves helps reduce the impact of external interference on the brain, allowing the subjects to focus more on the task of motor imagery. This improvement in concentration helps motor imagery become clearer and more accurate.

[0149] (3) Affects subsequent feature extraction and classification. In brain-computer interface (BCI) technology, feature extraction and classification of α waves is one of the key steps in subsequent processing. By analyzing and processing α waves, feature information related to motor imagery can be extracted, providing an important basis for subsequent decision-making and control.

[0150] (4) Improve recognition accuracy. The stability and regularity of α waves help improve the recognition accuracy of motor imagery EEG signals.

[0151] Figure 7 It shows the 6-channel right foot movement imagination time series after FIR filter 8-15Hz filtering. Figure 8 It shows an example of EEG time series data of the right foot C4 channel after FIR filter 8-15Hz filtering.

[0152] S5: Obtain the power spectral density PSD-frequency curve of the above α-wave motor imagery EEG signal.

[0153] Fig. 9 shows an example of C4 channel PSD-frequency samples, Fig.11 An example of Cz channel PSD-frequency samples is shown.

[0154] S6: Using the basic principle of functional data analysis (FDA), the PSD-frequency sample curves of all channels are functionalized, and the functionalization process is performed according to the channels.

[0155] When using FDA to analyze PSD-frequency sample curve data, the first step is to fit the discrete observation data into a continuous function. Let the i-th (i=1,2,…,N) observation sample be {(y ij ,t j )|j=1,2,…,n}, then:

[0156] y ij =x i (t j )+e ij ,j=1,2,…,n(1)

[0157] In the formula, x i (t j ) is t j Function x of frequency i (t) (PSD-frequency function), in the present invention, t j It represents the frequency value, e ij is the corresponding fitting error, n is the number of points selected in the frequency coordinate, in the present invention n=1251, one point is sampled every 0.2 Hz (the sampling rate of the EEG device is 500 Hz), and the frequency range is 0-250 Hz.

[0158] Use the basis function expansion method to convert x i (t) represents the linear combination of basis functions

[0159]

[0160] In the formula, c i =[c i1 ,c i2 ,…,c iU ] T is the coefficient vector of the basis function, φ=[φ 1 (t),φ 2 (t),…,φ U (t)] T is a functional vector composed of U basis functions. MI-EEG signals are non-periodic, so the B-spline basis function suitable for non-periodic data is selected, which has good differentiability and locality.

[0161] After selecting the basis function, the coefficient vector c is estimated by minimizing the sum of squared errors using the weighted least squares criterion. iThe larger the number of solved time-based functions, the better the data fit, but the worse the function smoothness, and the noise or abnormal data in the original data may be fitted into the function. In order to take into account both the data fit and the function smoothness, the a-order roughness penalty function is introduced

[0162]

[0163] Where: D a represents the a-order differential operator, J = ∫D a φD a φ T dt.

[0164] After the introduction of the rough penalty function, the estimated coefficient vector c i The weighted least squares criterion is

[0165]

[0166] Where: y i is the vector form of the observed sample value, Φ is the element φ u (t ij ), W is the weighting matrix, and λ is the smoothing parameter.

[0167] Coefficient vector c i The estimated value is

[0168]

[0169] This results in a functional representation of the sample data, and then the function analysis method can be used for the next step of data mining.

[0170] In the example of the present invention, there are 6 channels with 4×120 four-category (left hand, right hand, left foot, right foot motor imagination) MI-EEG functionalized sample data, and the B-spline order of each channel is 6. 20% to 80% of them are randomly selected as training samples, and the remaining samples are classified and predicted.

[0171] Fig.10 It shows an example of C4 channel PSD-frequency sample functionalization, Fig.12 An example of Cz channel PSD-frequency sample functionalization is shown.

[0172] S7: Perform functional principal component analysis (FPCA) on the functionalized PSD-frequency sample curve to obtain the corresponding principal component score eigenvalue.

[0173] When performing FPCA, for the convenience of representation, the N sample functions are still recorded as functions {x i(t), i = 1, 2, ..., N}. The first principal component function (PCF) of the sample function refers to the characteristic function ξ with the L2 norm of 1 that maximizes the variance of the projection of the sample function on it. 1 (t). The sample function is 1 The projection onto (t) is called the first functional principal component score (FPCS).

[0174]

[0175] So solving ξ 1 The problem (t) is expressed as

[0176]

[0177] Solve for the mth characteristic function ξ m (t) is to satisfy ξ m (t) maximizes the mean square of the m-th order FPCS when it is orthogonal to the previous m-1 characteristic functions, that is,

[0178]

[0179] The characteristic function ξ(t) satisfies the characteristic equation (7)

[0180] ∫v(s,t)ξ(t)dt=ρξ(s) (9)

[0181] In the formula: s, t are independent variables, is the covariance function and ρ is the eigenvalue.

[0182] N sample functions x = [x 1 (t),x 2 (t),…,x N (t)] T The basis function expansion of is represented by x=Cφ, where C is the N×U-order basis function coefficient matrix. The basis function expansion of the characteristic function ξ(t) is ξ(t)=b T φ,b=[b 1 ,b 2 ,…,b U ] T is the coefficient vector. Define the elements as the U-order symmetric matrix whose inner products are basis functions, i.e. Γ=∫φφ T dt; use A to represent the coefficient vector c i The covariance matrix of -1 C T C, then formula (9) can be expressed as

[0183] AΓb=ρb (10)

[0184] Solving the coefficient vector b gives the characteristic function ξ(t). In order to prevent the principal component function from being too rough, the same roughness penalty function as in the data functionization process is used to smooth the principal component function, that is, to maximize the sample variance with roughness penalty:

[0185]

[0186] The basis function expansion form of the above formula is

[0187]

[0188] The generalized eigenvalue problem is

[0189] ΓAΓb=ρ(Γ+λJ)b (13)

[0190] Perform Cholesky decomposition Γ+λJ=SS T , transformed into a symmetric matrix eigenvalue analysis problem:

[0191] (S -1 ΓAΓ(S -1 ) T )(S T b)=ρ(S T b) (14)

[0192] Solving the above equation, we can get the principal component function (PCF), and retain the first M, denoted as ξ(t) = [ξ 1 (t),ξ 2 (t),…,ξ M (t)], ξ(t) can be regarded as a set of M-dimensional orthogonal basis functions generated inside the functional data and capable of representing the original data. In the example of the present invention, M=10.

[0193] The sample function x i (t) is projected onto ξ(t) to obtain the characteristic parameters We call it the principal component score eigenvalue.

[0194] For the example of the present invention, each channel collects 4 types of MI-EEG data, each type has 120 sample data, and each sample PSD-frequency function has 10 principal component score eigenvalues, so the principal component scores of all sample curve data of 6 channels are 480×60 matrices. Here, the original data has been greatly reduced in dimension.

[0195] Fig.13 represents the 10th-order principal component function (PCF) of the C3 channel sample, Fig.14The 10th-order principal component function (PCF) of the Cz channel samples is shown. Fig.17 shows examples of unnormalized principal component score eigenvalues.

[0196] S8: Perform functional principal differential analysis (PDA) on the functionalized PSD-frequency sample curve to obtain the corresponding principal differential score eigenvalue (PDS).

[0197] Principal Differential Analysis (PDA) is a concise functional analysis method that models the observed data through a linear differential equation. Using differential equations as modeling objects involves related concepts in functional linear models and principal component analysis. For a single function, the differential equation can capture the shape characteristics of the function and its derivatives; for multiple functions, since the solution set of the differential equation is a complete function space with the characteristics of principal component analysis, a function subspace can be found from it to capture the main change characteristics in the sample data. Based on this feature, the present invention expands PDA into a feature extraction technology.

[0198] S8.1: Linear differential operators and differential equations

[0199] In general, the K-order homogeneous differential equation is expressed as:

[0200]

[0201] In the formula, is the coefficient parameter, which can be a variable function or a constant.

[0202] Define a linear differential operator:

[0203]

[0204] Where I is the identity operator. The homogeneous differential equation associated with it is:

[0205]

[0206] The primary task of PDA is to estimate the linear differential operator L, hoping that the linear differential operator L can eliminate the given N functions x(t) = [x 1 (t),x 2 (t),…,x N (t)] T , that is, we hope that the fitting residual Lx provided by the differential equation (17) i (t) is as close to 0 as possible. The residual norm sum of squares is defined as the fitting criterion using the least squares method:

[0207]

[0208] If the linear differential equation is non-homogeneous, that is, there is a forcing function δ(t) representing the external input or disturbance of the system, then the differential equation is expressed as:

[0209]

[0210] At this time, Lx(t) = δ(t), and the least squares fitting criterion of the linear differential operator L becomes:

[0211]

[0212] S8.2: Estimation of weight function

[0213] Based on the least squares fitting criterion of formula (20), the coefficient parameters of the linear differential operator L are estimated using the point-by-point minimization method. The fitting criterion of the point-by-point estimation is:

[0214]

[0215] In the formula, When t is fixed, the fitting criterion is just a simple least squares fitting criterion. It can be expressed in matrix form:

[0216]

[0217] In the formula, τ(t)=[D K x 1 (t),D K x 2 (t),…,D K x N (t)] T ; Z(t) is the point-by-point estimation design matrix, whose i-th row element is z i (t) = [-x i (t),-Dx i (t),…,-D K-1 x i (t),δ i (t)].

[0218] Keeping t constant, the coefficient vector can be obtained by minimizing point by point The least squares solution of is:

[0219]

[0220] It should be noted that the coefficient vector is estimated pointwise The existence of depends on Z(t) T The determinant of Z(t) is non-zero for any t.

[0221] The estimate of the linear differential operator L is obtained as follows:

[0222]

[0223] RSQ is a commonly used tool in PDA to measure the goodness of fit of differential equations:

[0224]

[0225] The value range of RSQ is [0,1]. The closer its value is to 1, the better the linear differential operator L fits the function.

[0226] S8.3: Principal differential score eigenvalue calculation

[0227] Theorem: Let Λ be the Hilbert space, Y be the linear closed subspace of Λ, Y ⊥ is the orthogonal complement of Y. Any element belonging to Λ can be uniquely decomposed into the component belonging to Y and the component belonging to Y ⊥ The sum of the weights.

[0228] According to the theorem, the function The function space can be orthogonally decomposed into two parts, denoted as and in express The null space of for This shows that the function x i (t) can be orthogonally decomposed into in It is called the structural component and satisfies It combines the common functional differential characteristics of functional data; It is called the residual component and satisfies It reflects the linear differential operator The differential characteristics of functions cannot be described.

[0229] When PDA is used as a feature extraction technique, the coefficient parameter of the differential operator L is considered to be a constant, and it is assumed that there is no forcing function. For a K-order non-degenerate linear differential equation, any solution can be expressed as a linear combination of K linearly independent complex exponential functions. Construct the characteristic equation:

[0230]

[0231] The characteristic equation has K linearly independent complex roots The corresponding complex exponential function is:

[0232]

[0233] The space spanned by K linearly independent complex exponential functions as a linear differential operator The null space of:

[0234]

[0235] function A set of basis functions are formed in the solution space of the differential equation, which are called principal differential functions or differential characteristic functions. It can be expressed as:

[0236]

[0237] In the formula, is the function x i The kth principal differential score (PDS) of (t), which is the function x extracted by PDA i (t) is the characteristic parameter of

[0238] In Hilbert space, inner product calculation also satisfies the distributive law and the associative law, so:

[0239]

[0240] The principal differential score eigenvalue PDS is solved by the following formula:

[0241] f i PDA =C -1 d i (31) In the formula, f i PDA is the principal differential score vector, C is the matrix whose elements are the inner products of complex exponential functions, d i The element is a function x i The vector of the inner product of (t) and the complex exponential function is expressed as:

[0242]

[0243] Comparing PDA and FPCA, we can find that they are somewhat similar in form, both involving minimizing the sum of squared errors, but they are fundamentally different. FPCA is centered on the mean function and generates a projection space from the first C orthogonal characteristic functions of the characteristic equation of the covariance function, achieving data dimensionality reduction from a geometric perspective, and the lost information is the missing component in the selection of comprehensive indicators. PDA, on the other hand, fits the function by estimating the K-order linear differential operator L, and achieves data dimensionality reduction by constructing the null space of the linear differential operator L. The lost information is the fitting residual of the linear differential operator.

[0244] As an example, Fig.15It represents the differential characteristic function of the P4 channel sample function. Fig.16 The differential characteristic function of the Pz channel sample function is shown. In the example of the present invention, a third-order or fourth-order linear differential operator is calculated for the sample function.

[0245] For the example of the present invention, each channel (a total of 6 channels) collects 4 types of MI-EEG data, each type has 120 sample data, and each sample PSD-frequency function has 3 or 4 main differential score values, so the main differential score of all sample curve data of 6 channels is a 480×23 matrix. Here, the original data is greatly reduced in dimension.

[0246] Fig.18 An example of unnormalized principal differential score eigenvalues ​​is shown.

[0247] S9: Merge the principal component score eigenvalues ​​and the principal differential score eigenvalues ​​to form an eigenvalue sequence for each sample, and normalize all eigenvalues ​​to facilitate subsequent machine learning.

[0248] In the example introduced in the present invention, each channel collects 4 types of MI-EEG data, each type has 120 sample data, and the original total data set size is (4×120×2500)×6 dimensions. After the above feature extraction, each sample PSD-frequency function curve has only 10 principal component score eigenvalues ​​and 3 or 4 principal differential score eigenvalues, so the principal component score (eigenvalue) of all sample curve data of 6 channels is 480×83 dimensions. In this way, the original data is greatly reduced in dimension.

[0249] Fig.19 It shows an example of normalized training sample feature values ​​(when the training sample accounts for 70%, there are 336 training samples in total). Fig. 20 It shows an example of normalized feature values ​​of test samples (when the training samples account for 70%, there are 144 test samples in total).

[0250] S10: Clustering training and classification prediction of sparrow search algorithm optimized kernel extreme learning machine (SSA_HELE).

[0251] Extreme Learning Machine (ELM) is a special single hidden layer feedforward neural network. Different from the traditional single hidden layer feedforward neural network, ELM randomly generates the connection weights between the input layer and the hidden layer and the bias of the hidden layer. It only needs to set the number of hidden layer nodes to obtain the only optimal solution, which improves the learning speed and generalization performance while ensuring the learning accuracy. In the nearly 20 years since the birth of ELM theory, it has been widely studied and applied. Applying the kernel function idea to ELM, the kernel extreme learning machine (KELM) algorithm is formed, which performs better than ELM in many problems. Its network structure is as follows: Fig.21 shown.

[0252] The sparrow search algorithm takes the foraging process of sparrow groups as a reference and idealizes it into a finder-joiner model. The main rules are as follows:

[0253] (1) The discoverer has a higher energy reserve and is responsible for finding food. It provides the foraging area and direction for the population. The joiners follow the discoverer to get food. The energy reserve of the sparrow depends on the fitness value of the individual sparrow.

[0254] (2) Once a sparrow discovers a predator, the individual will send out an alarm signal. When the alarm value is greater than the safety value, the discoverer will take the joiner to other safe areas to forage for food.

[0255] (3) The identities of discoverers and joiners change dynamically, but the ratio between them remains constant.

[0256] (4) The lower the energy reserves of the joiner, the worse its foraging position in the population, and the more likely it is to fly to other places to forage for more energy.

[0257] (5) During the foraging process, the joiners can always search for the discoverer that provides the best food, and then forage around the discoverer. At the same time, some joiners will monitor the discoverer to obtain a better foraging area.

[0258] (6) When aware of danger, sparrows at the edge of the group will move quickly to a safe area, while sparrows in the middle of the group will move around randomly.

[0259] According to the above rules, suppose that in the M-dimensional search space, the population consists of W sparrows, of which the proportions of discoverers and joiners are c and c respectively. P and c S =1-c P , the proportion of vigilant people is c V , the lower and upper boundaries of the population are p L and p U, the position of the wth sparrow is p w =[p w1 ,p w2 ,...,p wM ] T , the fitness value is f(p i ).

[0260] During each iteration, the finder's position is updated as follows:

[0261]

[0262] Where, τ is the current iteration step; τ max is the maximum number of iterations; is the current position of the wth sparrow; R 0 ∈(0,1] is a random number; R 1 is a random value that follows a normal distribution; L is an M-dimensional column vector whose elements are all 1; R 2 ∈[0,1] is the warning value; ST∈[0.5,1] is the safety threshold.

[0263] The joiner position is updated as follows:

[0264]

[0265] In the formula, The best position currently occupied by the discoverer; is the current global worst position; A + =(A T A) - 1 A T , A is an M-dimensional column vector whose elements are randomly assigned 1 or -1.

[0266] The positions of the Sentinels have been updated as follows:

[0267]

[0268] In the formula, is the current global optimal position; α 0 is the step size control parameter, which is a random number that obeys the standard normal distribution; R 3 ∈[-1,1] is a random number; f w is the current individual fitness value of the sparrow, f best and f worst is the best and worst fitness value of the current group; ε 0 is a very small constant.

[0269] In the present invention, the sparrow search algorithm (SSA) is used to optimize the regularization coefficient and kernel function parameters of the kernel extreme learning machine (KELM). The optimization process is as follows: Fig. 22 shown.

[0270] In the example of the present invention, 20-80% of the samples are randomly selected as classification training, and the remaining samples are used as classification tests. Table 1 shows the recognition accuracy of the test set under different proportions of training samples.

[0271] Table 1 Test set recognition accuracy under different training sample ratios

[0272]

[0273] It can be seen that the MI-EEG signal classification method provided by the present invention is concise and efficient. When the training sample is only 20%, the minimum classification accuracy can reach 71.61%. When the training sample is 80%, the minimum classification accuracy can reach 85.42%.

[0274] Fig.23 Given an example confusion matrix, Fig.24 This is an example of a comparison chart between the true category and the predicted category (when the training accounts for 70%).

Claims

1. A motor imagery EEG signal classification method combining functional data analysis and kernel extreme learning machine, characterized in that: The following steps are involved:

1. Use high-precision EEG signal acquisition equipment with more than 8 channels to collect MI-EEG signals; According to the technical requirements of EEG signal preprocessing, the collected MI-EEG signals are preprocessed; 2. Establish data samples according to the duration of a single motor imagery, and obtain the power spectral density PSD-frequency relationship curve of the alpha band of multiple channels and multiple samples through FIR filter; 3. According to the basic theory of functional data analysis of FDA, all PSD-frequency relationship curves are functionally smoothed; 4. Use functional principal component analysis (FPCA) to obtain the principal component score of each data sample, which is called the principal component score eigenvalue; 5. Use functional principal differential analysis PDA to obtain the principal differential score of each sample, which is called the principal differential score eigenvalue; 6. Combine the principal component score eigenvalue and principal differential score eigenvalue of each sample, and normalize the characteristic parameters of all samples; 7. Randomly select training samples and test samples, and use the sparrow search algorithm to optimize the kernel extreme learning machine SSA_HELM for clustering training; 8. Make category predictions for samples to be classified.

2. The method for classifying motor imagery EEG signals combining functional data analysis and kernel extreme learning machine according to claim 1, characterized in that: The preprocessing of the collected MI-EEG signal in step 1 specifically includes: Signal filtering: bandpass filtering to retain the most valuable frequency range for motor imagery tasks, which is between 1 and 40 Hz; notch filtering, which uses notch filters to suppress noise at specific frequencies; Artifact removal: by decomposing the independent components in the MI-EEG signal and removing the components related to eye movements and muscle activity; Channel selection: Select the EEG channel that is most valuable for the motor imagery task, and select the MI-EEG signals collected by the movement-related electrodes C4, P4, Pz, Cz, P3, and C3; Baseline correction: The MI-EE signal in the period before stimulus presentation was used as the baseline, and this baseline value was subtracted from the data at the subsequent time. Signal amplification: amplify the MI-EEG signal.

3. The method for classifying motor imagery EEG signals combining functional data analysis and kernel extreme learning machine according to claim 1, characterized in that: In step 1, a sampling frequency of 500 Hz is selected to collect MI-EEG signals.

4. The method for classifying motor imagery EEG signals combining functional data analysis and kernel extreme learning machine according to claim 1, characterized in that: The specific implementation of step 2 is: 2.1 Obtain the alpha wave time series in the range of 8-14 Hz of the preprocessed MI-EEG signal through the FIR filter; 2.2 For the above time series, establish data samples according to the duration of a single motor imagery; 2.3 For each sample time series, perform FFT operation to obtain the PSD-frequency relationship curve of each α-wave motor imagery EEG signal mentioned above.

5. The method for classifying motor imagery EEG signals combining functional data analysis and kernel extreme learning machine according to claim 1, characterized in that: The method for functional smoothing of all PSD-frequency relationship curves described in step 3 is: 3.1 According to the brain electrode channel, the discrete PSD-frequency data of each channel is fitted into a continuous function one by one; the sample data belonging to the C4, P4, Pz, Cz, P3, and C3 channels are fitted with functions respectively, that is, the four types of sample data of each channel are put together for function fitting; MI-EEG signals are non-periodic, so we select B-spline basis functions suitable for non-periodic data, and take 6-8th order B-splines, so that the fitting curve has good differentiability and locality; 3.2 After selecting the basis function, the coefficient vector c is estimated by minimizing the sum of squared errors using the weighted least squares criterion i ; The larger the number of solved time-based functions, the better the data fitting, but the worse the function smoothness, and the noise or abnormal data in the original data may be fitted into the function; in order to take into account both the data fitting and the smoothness of the function, the a-order roughness penalty function is introduced.

6. The method for classifying motor imagery EEG signals combining functional data analysis and kernel extreme learning machine according to claim 1, characterized in that: In step 4, functional principal component analysis FPCA is used to obtain the principal component score of each sample, which is specifically implemented as follows: 4.1 Covariance function estimation; Estimating the covariance function of functional data, which describes the correlation between data changes at different time points or variable values; 4.2 Solving eigenvalues ​​and eigenfunctions: By solving the eigenvalue problem of the covariance function, a series of eigenvalues ​​and corresponding eigenfunctions are obtained; these eigenvalues ​​and eigenfunctions represent the main change patterns in the data and the importance of these patterns in the data respectively; 4.3 Sorting and selection: sort the eigenvalues ​​from large to small, and select the eigenfunctions corresponding to the first M largest eigenvalues, where M is the number of dimensions after dimensionality reduction; 4.4 Projection and dimensionality reduction; Projection is to project the original functional data onto a low-dimensional subspace defined by the selected M characteristic functions. This is achieved by calculating the inner product of the original data and each characteristic function. The projection is the eigenvalue used for future classification, called the principal component score eigenvalue; the result after projection is the dimensionality reduced data, which represents the main change pattern in the original data in a lower dimension.

7. The method for classifying motor imagery EEG signals combining functional data analysis and kernel extreme learning machine according to claim 1, characterized in that: In step 5, the functional principal differential analysis PDA is used to obtain the principal differential score of each sample, which is specifically implemented as follows: (1) Construction of linear differential operators and differential equations The K-order homogeneous differential equation is expressed as: In the formula, is the coefficient parameter, which can be a variable function or a constant; Define a linear differential operator: Where I is the identity operator; the homogeneous differential equation associated with it is: The primary task of PDA is to estimate the linear differential operator L, hoping that the linear differential operator L can eliminate the given N functions x(t) = [x1(t), x2(t), …, x N (t)] T , that is, we hope that the fitting residual Lx provided by the differential equation (17) i (t) is as close to 0 as possible; the residual norm sum of squares is defined as the fitting criterion using the least squares method: If the linear differential equation is non-homogeneous, that is, there is a forcing function δ(t) representing the external input or disturbance of the system, then the differential equation is expressed as: At this time, Lx(t) = δ(t), and the least squares fitting criterion of the linear differential operator L becomes: (2) Estimation of weight function Based on the least squares fitting criterion of formula (6), the coefficient parameters of the linear differential operator L are estimated using the point-by-point minimization method. The fitting criterion of the point-by-point estimation is: In the formula, When t is fixed, the fitting criterion is just a simple least squares fitting criterion, which can be expressed in matrix form: In the formula, τ(t)=[D K x1(t),D K x2(t),…,D K x N (t)] T ; Z(t) is the point-by-point estimation design matrix, whose i-th row element is z i (t) = [-x i (t),-Dx i (t),…,-D K-1 x i (t),δ i (t)]; Keeping t constant, the coefficient vector can be obtained by minimizing point by point The least squares solution of is: Pointwise estimates of coefficient vectors The existence of depends on Z(t) T The determinant of Z(t) is non-zero for any t; The estimate of the linear differential operator L is obtained as follows: RSQ is a commonly used tool in PDA to measure the goodness of fit of differential equations: The value range of RSQ is [0,1]. The closer its value is to 1, the better the linear differential operator L fits the function. (3) Obtaining the principal differential score eigenvalue Theorem: Let Λ be the Hilbert space, Y be the linear closed subspace of Λ, Y ⊥ is the orthogonal complement of Y. Any element belonging to Λ can be uniquely decomposed into a component belonging to Y and a component belonging to Y. ⊥ The sum of the weights; According to the theorem, the function The function space can be orthogonally decomposed into two parts, denoted as and in express The null space of for orthogonal complement of i (t) can be orthogonally decomposed into in It is called the structural component and satisfies It combines the common functional differential characteristics of functional data; It is called the residual component and satisfies Linear differential operator The differential characteristics of functions cannot be described; When PDA is used as a feature extraction technique, the coefficient parameter of the differential operator L is considered to be a constant, and it is assumed that there is no forcing function; for a K-order non-degenerate linear differential equation, any solution can be expressed as a linear combination of K linearly independent complex exponential functions; construct the characteristic equation: The characteristic equation has K linearly independent complex roots The corresponding complex exponential function is: The space spanned by K linearly independent complex exponential functions as a linear differential operator The null space of: function A set of basis functions are formed in the solution space of the differential equation, which are called principal differential functions or differential characteristic functions. It is expressed as: In the formula, is the function x i The k-th principal differential score PDS of (t), that is, the function x extracted by PDA i (t) is the characteristic parameter of

Images