阅读说明：本技术 一种基于itd模态参数识别的脑电信号特征提取方法 (Electroencephalogram signal feature extraction method based on ITD modal parameter identification ) 是由 杨怀花 叶庆卫 于 2021-06-08 设计创作，主要内容包括：本发明公开了一种基于ITD模态参数识别的脑电信号特征提取方法,其根据脑电信号中的δ波节律信号、θ波节律信号、α波节律信号、β波节律信号各自的频率范围,设计各自对应的有限脉冲响应带通滤波器,再利用各自对应的有限脉冲响应带通滤波器对脑电信号进行滤波处理,提取得到δ波节律信号、θ波节律信号、α波节律信号、β波节律信号；对四个节律信号分别进行ITD模态参数识别,提取得到各自的四个ITD模态参数,进而得到四个节律信号各自的特征向量；将四个节律信号的特征向量构成脑电信号的特征向量；优点是其能够从运动想象的脑电信号中提取出有效的特征,从而能够极大地提高后续特征分类的正确率。(The invention discloses an electroencephalogram signal characteristic extraction method based on ITD modal parameter identification, which comprises the steps of designing respective corresponding finite impulse response band-pass filters according to respective frequency ranges of a delta nodal signal, a theta nodal signal, an alpha nodal signal and a beta nodal signal in an electroencephalogram signal, then carrying out filtering processing on the electroencephalogram signal by using the respective corresponding finite impulse response band-pass filters, and extracting to obtain the delta nodal signal, the theta nodal signal, the alpha nodal signal and the beta nodal signal; respectively carrying out ITD modal parameter identification on the four rhythm signals, extracting to obtain four respective ITD modal parameters, and further obtaining respective characteristic vectors of the four rhythm signals; forming feature vectors of the four rhythm signals into feature vectors of the electroencephalogram signals; the method has the advantages that effective features can be extracted from the electroencephalogram signals of the motor imagery, so that the accuracy of subsequent feature classification can be greatly improved.)

1. An EEG signal feature extraction method based on ITD modal parameter identification is characterized by comprising the following steps:

step 1: randomly acquiring an electroencephalogram signal in one channel at one time according to a given sampling frequency; then, the electroencephalogram signals are expressed in a vector form and are marked as Y, Y ═ Y₁…y_n…y_N](ii) a The given sampling frequency is F hertz, the dimension of Y is 1 XN, each element in Y is a sample value, N represents the total number of sample values in Y, N is not less than 5, N is F X t, t represents the acquisition time for acquiring one electroencephalogram signal of one channel, the unit of t is second, N is a positive integer, the initial value of N is 1, N is not less than 1 and not more than N, Y is a positive integer, N is a positive integer, the number of the sampling frequency is larger than 1X, the number of the sampling frequency is larger than N, each element in Y is one, the number of the sampling frequency is one, N represents the total number of the sampling values in Y, N is not less than 5, N is larger than 1, N is not less than F X, t represents the acquisition time for acquiring one electroencephalogram signal of one channel, the electroencephalogram signal, the unit of one channel, the unit of second, N is a unit, N is a positive integer, N is larger than 1, N is larger than N, N is larger than N, N is larger than N₁Represents the 1 st sample value in Y, Y_nRepresenting the n-th sample value in Y, Y_NDenotes the Nth sample value in Y, the symbol "[ alpha ]]"is a symbol representing a vector or matrix;

step 2: designing respective corresponding finite impulse response band-pass filters according to respective frequency ranges of a delta wave pitch signal, a theta wave pitch signal, an alpha wave pitch signal and a beta wave pitch signal in the electroencephalogram signals; then, respectively filtering Y by utilizing finite impulse response band-pass filters corresponding to delta wave pitch signals, theta wave pitch signals, alpha wave pitch signals and beta wave pitch signals in the electroencephalogram signals, correspondingly extracting the delta wave pitch signals, the theta wave pitch signals, the alpha wave pitch signals and the beta wave pitch signals in the Y, and correspondingly expressing the delta wave pitch signals, the theta wave pitch signals, the alpha wave pitch signals and the beta wave pitch signals in the Y in a vector form as Y_δ、Y_θ、Y_α、Y_β，Y_δ＝[y_δ,1…y_δ,n…y_δ,N]，Y_θ＝[y_θ,1…y_θ,n…y_θ,N]，Y_α＝[y_α,1…y_α,n…y_α,N]，Y_β＝[y_β,1…y_β,n…y_β,N](ii) a Wherein the frequency range of delta nodal signals in the electroencephalogram signals is 0.78-3.90 Hz, the frequency range of theta nodal signals in the electroencephalogram signals is 3.91-7.80 Hz, the frequency range of alpha nodal signals in the electroencephalogram signals is 7.81-13.28 Hz, the frequency range of beta nodal signals in the electroencephalogram signals is 13.29-30.47 Hz, and the frequency range of Y nodal signals in the electroencephalogram signals is 13.29-30.47 Hz_δ、Y_θ、Y_α、Y_βAre all of dimensions1×N，Y_δIs a sample value, y_δ,1Represents Y_δ1 st sample value of (1), y_δ,nRepresents Y_δN-th sample value of (1), y_δ,NRepresents Y_δSample value of Nth, Y_θIs a sample value, y_θ,1Represents Y_θ1 st sample value of (1), y_θ,nRepresents Y_θN-th sample value of (1), y_θ,NRepresents Y_θSample value of Nth, Y_αIs a sample value, y_α,1Represents Y_α1 st sample value of (1), y_α,nRepresents Y_αN-th sample value of (1), y_α,NRepresents Y_αSample value of Nth, Y_βIs a sample value, y_β,1Represents Y_β1 st sample value of (1), y_β,nRepresents Y_βN-th sample value of (1), y_β,NRepresents Y_βThe nth sample value;

and step 3: for Y_δIdentifying ITD modal parameters, and extracting to obtain Y_δFour ITD mode parameters of, will Y_δThe vector formed by the sequential arrangement of the four ITD modal parameters is used as Y_δFeature vector of (1), denoted fea_δ(ii) a In the same way for Y_θIdentifying ITD modal parameters, and extracting to obtain Y_θFour ITD mode parameters of, will Y_θThe vector formed by the sequential arrangement of the four ITD modal parameters is used as Y_θFeature vector of (1), denoted fea_θ(ii) a In the same way for Y_αIdentifying ITD modal parameters, and extracting to obtain Y_αFour ITD mode parameters of, will Y_αThe vector formed by the sequential arrangement of the four ITD modal parameters is used as Y_αFeature vector of (1), denoted fea_α(ii) a In the same way for Y_βIdentifying ITD modal parameters, and extracting to obtain Y_βFour ITD mode parameters of, will Y_βThe vector formed by the sequential arrangement of the four ITD modal parameters is used as Y_βFeature vector of (1), denoted fea_β(ii) a Wherein fea_δ、fea_θ、fea_α、fea_βThe dimensions of (A) are all 1 × 4;

and 4, step 4: will fea_δ、fea_θ、fea_α、fea_βThe vector formed by arranging the sequences is taken as a feature vector of Y and is denoted as fea, fea ═ fea_δ fea_θ fea_α fea_β](ii) a Wherein the dimension of fea is 1 × 16.

2. The electroencephalogram signal feature extraction method based on ITD modal parameter identification according to claim 1, characterized in that the specific process of the step 2 is as follows:

according to the frequency range of a delta-wave pitch signal in an electroencephalogram signal, designing a first finite impulse response band-pass filter, wherein the first finite impulse response band-pass filter is recorded as H1, the sampling frequency of H1 is F Hz, the order of H1 is 10, the upper edge frequency of a stop band of H1 is 0.7 Hz, the lower edge frequency of the stop band is 4 Hz, the upper edge frequency of a pass band is 0.8 Hz, and the lower edge frequency of the pass band is 3.5 Hz; then, H1 is used for carrying out filtering processing on Y, and delta wave pitch signals in Y are extracted;

similarly, according to the frequency range of a theta-wave pitch signal in the electroencephalogram signal, designing a second finite impulse response band-pass filter, which is recorded as H2, wherein the sampling frequency of H2 is F Hz, the order of H2 is 10, the upper edge frequency of a stop band of H2 is 3.9 Hz, the lower edge frequency of the stop band is 8 Hz, the upper edge frequency of a pass band is 4 Hz, and the lower edge frequency of the pass band is 7.5 Hz; then, filtering Y by using H2, and extracting to obtain a theta node law signal in Y;

according to the frequency range of alpha-nodal signals in the electroencephalogram signals, a third finite impulse response band-pass filter is designed and recorded as H3, the sampling frequency of H3 is F Hz, the order of H3 is 10, the upper edge frequency of a stop band of H3 is 7.8 Hz, the lower edge frequency of the stop band is 13 Hz, the upper edge frequency of a pass band is 7.9 Hz, and the lower edge frequency of the pass band is 12 Hz; then, filtering Y by using H3, and extracting to obtain an alpha wave pitch signal in Y;

designing a fourth finite impulse response band-pass filter according to the frequency range of beta-nodal signals in the electroencephalogram signals, wherein the fourth finite impulse response band-pass filter is recorded as H4, the sampling frequency of H4 is F Hz, the order of H4 is 10, the upper edge frequency of a stop band of H4 is 13 Hz, the lower edge frequency of the stop band is 31 Hz, the upper edge frequency of a pass band is 14 Hz, and the lower edge frequency of the pass band is 30 Hz; and then, filtering Y by using H4 to extract and obtain a beta wave law signal in Y.

3. The EEG signal feature extraction method based on ITD modal parameter identification according to claim 1 or 2, characterized in that in step 3, fea_δThe acquisition process comprises the following steps:

step 3_ 1: will Y_δ＝[y_δ,1…y_δ,n…y_δ,N]The discrete signal is regarded as a free response vibration signal for Y_δ＝[y_δ,1…y_δ,n…y_δ,N]Carrying out normal sampling, d-delay sampling and 2 d-delay sampling in sequence to obtain a normal sampling free response data vector, a d-delay sampling free response data vector and a 2 d-delay sampling free response data vector in sequence, and correspondingly marking as X_nor、X_{del_1}、X_{del_2}，X_nor＝[y_δ,1…y_δ,n…y_δ,m]，X_{del_1}＝[y_δ,d+1…y_δ,d+n…y_δ,d+m]，X_{del_2}＝[y_δ,2d+1…y_δ,2d+n…y_δ,2d+m](ii) a Then X is put in_nor＝[y_δ,1…y_δ,n…y_δ,m]Re-expressed as X_nor＝[x_nor,1…x_nor,n…x_nor,m]Is mixing X_{del_1}＝[y_δ,d+1…y_δ,d+n…y_δ,d+m]Re-expressed as X_{del_1}＝[x_{del_1,1}…x_{del_1,n}…x_{del_1,m}]Is mixing X_{del_2}＝[y_δ,2d+1…y_δ,2d+n…y_δ,2d+m]Re-expressed as X_{del_2}＝[x_{del_2,1}…x_{del_2,n}…x_{del_2,m}](ii) a Wherein d is a positive integer, d represents a sampling interval,m is a positive integer, m represents X_nor、X_{del_1}、X_{del_2}M is more than or equal to 1 and less than or equal to N-2d, X_nor、X_{del_1}、X_{del_2}All dimensions of (a) are 1 xm, y_δ,m、y_δ,d+1、y_δ,d+n、y_δ,d+m、y_δ,2d+1、y_δ,2d+n、y_δ,2d+mCorresponds to and represents Y_δM-th, d + 1-th, d + n-th, d + m-th, 2d + 1-th, 2d + n-th, 2d + m-th sample values, x_nor,1、x_nor,n、x_nor,mCorresponds to X_nor1 st element, nth element, mth element, x_{del_1,1}、x_{del_1,n}、x_{del_1,m}Corresponds to X_{del_1}1 st element, nth element, mth element, x_{del_2,1}、x_{del_2,n}、x_{del_2,m}Corresponds to X_{del_2}The 1 st element, the nth element, the mth element in (a);

step 3_ 2: according to X_norAnd X_{del_1}Constructing a first augmented matrix, denoted X_zg1，And according to X_{del_1}And X_{del_2}Constructing a second augmented matrix, denoted X_zg2，Then let A ═ X_zg2(X_zg1)^-1And using QR decomposition method to pair A ═ X_zg2(X_zg1)^-1Processing to obtain a characteristic value matrix of A, and recording the characteristic value matrix as delta; wherein, X_zg1And X_zg2All dimensions of (A) are 2X m, A represents a system feature matrix, (X)_zg1)^-1Is X_zg1The dimension of a is 2 × 2, and the dimension of Δ is 2 × 2;

step 3_ 3: calculating Y from Delta_δHas a natural frequency and an attenuation coefficient, corresponding to f_δAnd λ_δ，Wherein the content of the first and second substances,to representThe imaginary part of (a) is,to representReal part of, Δ_1,1Line 1, column 1 elements representing Δ, t' at Y_δ＝[y_δ,1…y_δ,n…y_δ,N]A time interval of d sample values is sampled,

step 3_ 4: applying free response theory to X_norThe description is as follows: x_nor＝ψ_norE_nor，Wherein psi_norRepresents X_norMode vector of ψ_norDimension of (1X 2, E)_norFor the introduction of intermediate variables, E_norDimension of (d) is 2 xm, e represents a natural base, and μ ═ λ_δ+j2πf_δ，μ^*Is a conjugate of mu, mu^*＝-λ_δ-j2πf_δJ represents an imaginary unit of number,then solving X by using a least square method_nor＝ψ_norE_norTo obtain psi_nor，ψ_nor＝X_norE_nor ^T(E_norE_nor ^T)^-1Wherein, the upper mark is 'T' tableRepresenting transposes of vectors or matrices, (E)_norE_nor ^T)^-1Is E_norE_nor ^TThe inverse of (1); recalculate psi_norIs the amplitude and phase of (1), the correspondence is expressed as xi_norAndξ_nor＝abs(ψ_nor(1))+abs(ψ_nor(2))，wherein abs () is the absolute value function, ψ_nor(1) To indicate psi_norElement # 1 of_nor(2) To indicate psi_norThe number 2 element of (a) is,to representThe real part of (2); finally xi is_norAs Y_δAmplitude of (1), re-noted as xi_δWill beAs Y_δIs newly recorded as

Step 3_ 5: f. of_δ、λ_δ、ξ_δ、Is Y_δFour ITD modal parameters of, will f_δ、λ_δ、ξ_δ、Sequentially arranging the constructed vectors as Y_δFeature vector fea of_δ，

In the step 3, fea_θ、fea_α、fea_βIs acquired in the same manner as the procedure of step 3_1 to step 3_5, wherein f is_θ、λ_θ、ξ_θ、Corresponds to Y_θNatural frequency, attenuation coefficient, amplitude, phase, f_α、λ_α、ξ_α、Corresponds to Y_αNatural frequency, attenuation coefficient, amplitude, phase, f_β、λ_β、ξ_β、Corresponds to Y_βNatural frequency, attenuation coefficient, amplitude, phase.

4. The EEG signal feature extraction method based on ITD modal parameter identification as claimed in claim 3, characterized in that EEG signals are randomly acquired for many times in the same channel according to a given sampling frequency, feature vectors of the EEG signals acquired for each time in the same channel are acquired in the same way according to the processes of step 1 to step 4, the feature vectors of the EEG signals randomly acquired for many times in the same channel form a feature matrix, and each row of vectors in the feature matrix represents the feature vector of the EEG signals randomly acquired for one time; and then classifying the characteristic matrix to obtain a motor imagery type corresponding to the electroencephalogram signals randomly acquired each time.

5. The EEG signal feature extraction method based on ITD modal parameter identification as claimed in claim 4, characterized in that for a plurality of different channels, EEG signals are collected at a plurality of channels at the same time according to a given sampling frequency, and feature vectors of the EEG signals collected at each channel are obtained in the same manner according to the processes of step 1 to step 4, and the feature vectors of the EEG signals collected at a plurality of channels at the same time are arranged in sequence to form a row vector; after multiple times of acquisition, obtaining a plurality of row vectors, arranging the row vectors in sequence to form a characteristic matrix, wherein each row vector in the characteristic matrix represents a row vector formed by arranging the characteristic vectors of the electroencephalogram signals acquired at multiple channels at one time in sequence; and then classifying the characteristic matrix to obtain a motor imagery type corresponding to the electroencephalogram signal acquired each time.

Technical Field

The invention relates to a signal feature extraction technology, in particular to an electroencephalogram signal feature extraction method based on ITD (Ibrahim Time domain) modal parameter identification.

Background

The Brain-Computer Interface (BCI) is a communication control system that can realize direct communication between the human Brain and the external environment, and the communication control system can make the human Brain independent of the normal output pathway of the Brain composed of peripheral nerves and muscles of the Brain. In recent years, due to the widespread development of many technologies such as computer science, biomedicine, signal processing, and artificial intelligence, attention has been focused on brain-computer interface technology as a cross-discipline thereof.

The brain-computer interface technology comprises five main parts: the method comprises the steps of electroencephalogram signal acquisition, preprocessing, feature extraction, feature classification and identification and control equipment. The key point of the current technical research is how to effectively extract the features of the electroencephalogram signals and improve the accuracy of classification of the electroencephalogram signals. Hitherto, people extract features of electroencephalogram signals by using methods of a plurality of signal processing technologies, such as a CSP (Common Spatial Pattern) algorithm, a sparse self-coding method, and the like, so that quantitative analysis of electroencephalogram signals is greatly advanced, but whether dynamic behaviors exist in the brain is always a hot topic. From a kinetic point of view, the neuronal network system of the human brain is an extremely intricate, ultra-high-dimensional, nonlinear system formed by many neuronal connections. In view of the multiple coupling property of the life phenomenon, together with various nonlinear characteristics such as the activity and the large damping characteristic of the biological material, the living system has strong nonlinearity, so that the deep research can be carried out from the nonlinear characteristics, and the deep research has profound significance on the signal processing technology and the life phenomenon research.

Disclosure of Invention

The invention aims to solve the technical problem of providing an electroencephalogram characteristic extraction method based on the electroencephalogram signal nonlinear vibration theory, which can extract effective characteristics from electroencephalogram signals of motor imagery, thereby greatly improving the accuracy of subsequent characteristic classification.

The technical scheme adopted by the invention for solving the technical problems is as follows: an EEG signal feature extraction method based on ITD modal parameter identification is characterized by comprising the following steps:

step 1: randomly acquiring an electroencephalogram signal in one channel at one time according to a given sampling frequency; then, the electroencephalogram signals are expressed in a vector form and are marked as Y, Y ═ Y₁…y_n…y_N](ii) a Wherein the given sampling frequency is F Hz, and the dimension of YIs 1 multiplied by N, each element in Y is a sample value, N represents the total number of sample values in Y, N is more than or equal to 5, N is F multiplied by t, t represents the acquisition time of one electroencephalogram signal of one channel, the unit of t is second, N is a positive integer, the initial value of N is 1, N is more than or equal to 1 and less than or equal to N, Y₁Represents the 1 st sample value in Y, Y_nRepresenting the n-th sample value in Y, Y_NDenotes the Nth sample value in Y, the symbol "[ alpha ]]"is a symbol representing a vector or matrix;

step 2: designing respective corresponding finite impulse response band-pass filters according to respective frequency ranges of a delta wave pitch signal, a theta wave pitch signal, an alpha wave pitch signal and a beta wave pitch signal in the electroencephalogram signals; then, respectively filtering Y by utilizing finite impulse response band-pass filters corresponding to delta wave pitch signals, theta wave pitch signals, alpha wave pitch signals and beta wave pitch signals in the electroencephalogram signals, correspondingly extracting the delta wave pitch signals, the theta wave pitch signals, the alpha wave pitch signals and the beta wave pitch signals in the Y, and correspondingly expressing the delta wave pitch signals, the theta wave pitch signals, the alpha wave pitch signals and the beta wave pitch signals in the Y in a vector form as Y_δ、Y_θ、Y_α、Y_β，Y_δ＝[y_δ,1…y_δ,n…y_δ,N]，Y_θ＝[y_θ,1…y_θ,n…y_θ,N]，Y_α＝[y_α,1…y_α,n…y_α,N]，Y_β＝[y_β,1…y_β,n…y_β,N](ii) a Wherein the frequency range of delta nodal signals in the electroencephalogram signals is 0.78-3.90 Hz, the frequency range of theta nodal signals in the electroencephalogram signals is 3.91-7.80 Hz, the frequency range of alpha nodal signals in the electroencephalogram signals is 7.81-13.28 Hz, the frequency range of beta nodal signals in the electroencephalogram signals is 13.29-30.47 Hz, and the frequency range of Y nodal signals in the electroencephalogram signals is 13.29-30.47 Hz_δ、Y_θ、Y_α、Y_βAll dimensions of (A) are 1 XN, Y_δIs a sample value, y_δ,1Represents Y_δ1 st sample value of (1), y_δ,nRepresents Y_δN-th sample value of (1), y_δ,NRepresents Y_δSample value of Nth, Y_θIs a sample value, y_θ,1Represents Y_θ1 st sample value of (1), y_θ,nRepresents Y_θN-th sample value of (1), y_θ,NRepresents Y_θSample value of Nth, Y_αIs a sample value, y_α,1Represents Y_α1 st sample value of (1), y_α,nRepresents Y_αN-th sample value of (1), y_α,NRepresents Y_αSample value of Nth, Y_βIs a sample value, y_β,1Represents Y_β1 st sample value of (1), y_β,nRepresents Y_βN-th sample value of (1), y_β,NRepresents Y_βThe nth sample value;

and step 3: for Y_δIdentifying ITD modal parameters, and extracting to obtain Y_δFour ITD mode parameters of, will Y_δThe vector formed by the sequential arrangement of the four ITD modal parameters is used as Y_δFeature vector of (1), denoted fea_δ(ii) a In the same way for Y_θIdentifying ITD modal parameters, and extracting to obtain Y_θFour ITD mode parameters of, will Y_θThe vector formed by the sequential arrangement of the four ITD modal parameters is used as Y_θFeature vector of (1), denoted fea_θ(ii) a In the same way for Y_αIdentifying ITD modal parameters, and extracting to obtain Y_αFour ITD mode parameters of, will Y_αThe vector formed by the sequential arrangement of the four ITD modal parameters is used as Y_αFeature vector of (1), denoted fea_α(ii) a In the same way for Y_βIdentifying ITD modal parameters, and extracting to obtain Y_βFour ITD mode parameters of, will Y_βThe vector formed by the sequential arrangement of the four ITD modal parameters is used as Y_βFeature vector of (1), denoted fea_β(ii) a Wherein fea_δ、fea_θ、fea_α、fea_βThe dimensions of (A) are all 1 × 4;

and 4, step 4: will fea_δ、fea_θ、fea_α、fea_βThe vector formed by arranging the sequences is taken as a feature vector of Y and is denoted as fea, fea ═ fea_δfea_θfea_αfea_β](ii) a Wherein the dimension of fea is 1 × 16。

The specific process of the step 2 is as follows:

according to the frequency range of a delta-wave pitch signal in an electroencephalogram signal, designing a first finite impulse response band-pass filter, wherein the first finite impulse response band-pass filter is recorded as H1, the sampling frequency of H1 is F Hz, the order of H1 is 10, the upper edge frequency of a stop band of H1 is 0.7 Hz, the lower edge frequency of the stop band is 4 Hz, the upper edge frequency of a pass band is 0.8 Hz, and the lower edge frequency of the pass band is 3.5 Hz; then, H1 is used for carrying out filtering processing on Y, and delta wave pitch signals in Y are extracted;

similarly, according to the frequency range of a theta-wave pitch signal in the electroencephalogram signal, designing a second finite impulse response band-pass filter, which is recorded as H2, wherein the sampling frequency of H2 is F Hz, the order of H2 is 10, the upper edge frequency of a stop band of H2 is 3.9 Hz, the lower edge frequency of the stop band is 8 Hz, the upper edge frequency of a pass band is 4 Hz, and the lower edge frequency of the pass band is 7.5 Hz; then, filtering Y by using H2, and extracting to obtain a theta node law signal in Y;

according to the frequency range of alpha-nodal signals in the electroencephalogram signals, a third finite impulse response band-pass filter is designed and recorded as H3, the sampling frequency of H3 is F Hz, the order of H3 is 10, the upper edge frequency of a stop band of H3 is 7.8 Hz, the lower edge frequency of the stop band is 13 Hz, the upper edge frequency of a pass band is 7.9 Hz, and the lower edge frequency of the pass band is 12 Hz; then, filtering Y by using H3, and extracting to obtain an alpha wave pitch signal in Y;

designing a fourth finite impulse response band-pass filter according to the frequency range of beta-nodal signals in the electroencephalogram signals, wherein the fourth finite impulse response band-pass filter is recorded as H4, the sampling frequency of H4 is F Hz, the order of H4 is 10, the upper edge frequency of a stop band of H4 is 13 Hz, the lower edge frequency of the stop band is 31 Hz, the upper edge frequency of a pass band is 14 Hz, and the lower edge frequency of the pass band is 30 Hz; and then, filtering Y by using H4 to extract and obtain a beta wave law signal in Y.

In the step 3, fea_δThe acquisition process comprises the following steps:

step 3_ 1: will Y_δ＝[y_δ,1…y_δ,n…y_δ,N]Discrete signals are treated as free-response vibrationsSignal, to Y_δ＝[y_δ,1…y_δ,n…y_δ,N]Carrying out normal sampling, d-delay sampling and 2 d-delay sampling in sequence to obtain a normal sampling free response data vector, a d-delay sampling free response data vector and a 2 d-delay sampling free response data vector in sequence, and correspondingly marking as X_nor、X_{del_1}、X_{del_2}，X_nor＝[y_δ,1…y_δ,n…y_δ,m]，X_{del_1}＝[y_δ,d+1…y_δ,d+n…y_δ,d+m]，X_{del_2}＝[y_δ,2d+1…y_δ,2d+n…y_δ,2d+m](ii) a Then X is put in_nor＝[y_δ,1…y_δ,n…y_δ,m]Re-expressed as X_nor＝[x_nor,1…x_nor,n…x_nor,m]Is mixing X_{del_1}＝[y_δ,d+1…y_δ,d+n…y_δ,d+m]Re-expressed as X_{del_1}＝[x_{del_1,1}…x_{del_1,n}…x_{del_1,m}]Is mixing X_{del_2}＝[y_δ,2d+1…y_δ,2d+n…y_δ,2d+m]Re-expressed as X_{del_2}＝[x_{del_2,1}…x_{del_2,n}…x_{del_2,m}](ii) a Wherein d is a positive integer, d represents a sampling interval,m is a positive integer, m represents X_nor、X_{del_1}、X_{del_2}M is more than or equal to 1 and less than or equal to N-2d, X_nor、X_{del_1}、X_{del_2}All dimensions of (a) are 1 xm, y_δ,m、y_δ,d+1、y_δ,d+n、y_δ,d+m、y_δ,2d+1、y_δ,2d+n、y_δ,2d+mCorresponds to and represents Y_δM-th, d + 1-th, d + n-th, d + m-th, 2d + 1-th, 2d + n-th, 2d + m-th sample values, x_nor,1、x_nor,n、x_nor,mCorresponds to X_nor1 st element, nth element, mth element, x_{del_1,1}、x_{del_1,n}、x_{del_1,m}Correspond toRepresents X_{del_1}1 st element, nth element, mth element, x_{del_2,1}、x_{del_2,n}、x_{del_2,m}Corresponds to X_{del_2}The 1 st element, the nth element, the mth element in (a);

step 3_ 2: according to X_norAnd X_{del_1}Constructing a first augmented matrix, denoted X_zg1，And according to X_{del_1}And X_{del_2}Constructing a second augmented matrix, denoted X_zg2，Then let A ═ X_zg2(X_zg1)^-1And using QR decomposition method to pair A ═ X_zg2(X_zg1)^-1Processing to obtain a characteristic value matrix of A, and recording the characteristic value matrix as delta; wherein, X_zg1And X_zg2All dimensions of (A) are 2X m, A represents a system feature matrix, (X)_zg1)^-1Is X_zg1The dimension of a is 2 × 2, and the dimension of Δ is 2 × 2;

step 3_ 3: calculating Y from Delta_δHas a natural frequency and an attenuation coefficient, corresponding to f_δAnd λ_δ，Wherein the content of the first and second substances,to representThe imaginary part of (a) is,to representReal part of, Δ_1,1Line 1, column 1 elements representing Δ, t' atY_δ＝[y_δ,1…y_δ,n…y_δ,N]A time interval of d sample values is sampled,

step 3_ 4: applying free response theory to X_norThe description is as follows: x_nor＝ψ_norE_nor，Wherein psi_norRepresents X_norMode vector of ψ_norDimension of (1X 2, E)_norFor the introduction of intermediate variables, E_norDimension of (d) is 2 xm, e represents a natural base, and μ ═ λ_δ+j2πf_δ，μ^*Is a conjugate of mu, mu^*＝-λ_δ-j2πf_δJ represents an imaginary unit of number,then solving X by using a least square method_nor＝ψ_norE_norTo obtain psi_nor，ψ_nor＝X_norE_nor ^T(E_norE_nor ^T)^-1Wherein the superscript "T" denotes the transpose of the vector or matrix, (E)_norE_nor ^T)^-1Is E_norE_nor ^TThe inverse of (1); recalculate psi_norIs the amplitude and phase of (1), the correspondence is expressed as xi_norAndξ_nor＝abs(ψ_nor(1))+abs(ψ_nor(2))，wherein abs () is the absolute value function, ψ_nor(1) To indicate psi_norElement # 1 of_nor(2) To indicate psi_norThe number 2 element of (a) is,to representThe real part of (2); finally xi is_norAs Y_δAmplitude of (1), re-noted as xi_δWill beAs Y_δIs newly recorded as

Step 3_ 5: f. of_δ、λ_δ、ξ_δ、Is Y_δFour ITD modal parameters of, will f_δ、λ_δ、ξ_δ、Sequentially arranging the constructed vectors as Y_δFeature vector fea of_δ，

In the step 3, fea_θ、fea_α、fea_βIs acquired in the same manner as the procedure of step 3_1 to step 3_5, wherein f is_θ、λ_θ、ξ_θ、Corresponds to Y_θNatural frequency, attenuation coefficient, amplitude, phase, f_α、λ_α、ξ_α、Corresponds to Y_αNatural frequency, attenuation coefficient, amplitude, phase, f_β、λ_β、ξ_β、Corresponds to Y_βNatural frequency, attenuation coefficient, amplitude, phase.

Randomly acquiring electroencephalogram signals in the same channel for multiple times according to a given sampling frequency, acquiring the characteristic vectors of the electroencephalogram signals acquired in the same channel every time in the same mode according to the processes of the step 1 to the step 4, and forming a characteristic matrix by the characteristic vectors of the electroencephalogram signals acquired in the same channel for multiple times, wherein each row of vectors in the characteristic matrix represents the characteristic vector of the electroencephalogram signals acquired at one time; and then classifying the characteristic matrix to obtain a motor imagery type corresponding to the electroencephalogram signals randomly acquired each time.

Acquiring electroencephalogram signals at a plurality of channels according to a given sampling frequency aiming at a plurality of different channels, acquiring the characteristic vector of the electroencephalogram signals acquired at each channel in the same mode according to the processes from step 1 to step 4, and arranging the characteristic vectors of the electroencephalogram signals acquired at the plurality of channels at the same time in sequence to form a row vector; after multiple times of acquisition, obtaining a plurality of row vectors, arranging the row vectors in sequence to form a characteristic matrix, wherein each row vector in the characteristic matrix represents a row vector formed by arranging the characteristic vectors of the electroencephalogram signals acquired at multiple channels at one time in sequence; and then classifying the characteristic matrix to obtain a motor imagery type corresponding to the electroencephalogram signal acquired each time.

Compared with the prior art, the invention has the advantages that:

1) the method provided by the invention starts from a nonlinear vibration theory, analyzes the electroencephalogram signal dynamic model, researches the electroencephalogram signal from the aspect of vibration and engineering, has an important scientific theoretical value, and has an important significance for electroencephalogram signal characteristic extraction.

2) The characteristics of the electroencephalogram signals extracted by the method are used for classification, and simulation comparison is carried out on the existing motor imagery electroencephalogram signal database, so that the vibration mode characteristics acquired by the method have higher classification accuracy and better stability.

Drawings

Fig. 1 is a block diagram of the overall implementation of the method of the present invention.

Detailed Description

The invention is described in further detail below with reference to the accompanying examples.

The invention provides an electroencephalogram characteristic extraction method based on ITD modal parameter identification, the overall implementation block diagram of which is shown in figure 1, and the method comprises the following steps:

step 1: randomly acquiring an electroencephalogram signal in one channel at one time according to a given sampling frequency; then, the electroencephalogram signals are expressed in a vector form and are marked as Y, Y ═ Y₁…y_n…y_N](ii) a The given sampling frequency is F hertz, the given sampling frequency is a sampling frequency parameter of an acquisition system, the dimension of Y is 1 xN, each element in Y is a sample value, N represents the total number of the sample values in Y, N is not less than 5, N is F x t, t represents the acquisition time for acquiring one electroencephalogram signal of one channel, the unit of t is second, N is a positive integer, the initial value of N is 1, N is not less than 1 and not more than N, Y is not less than 1 and not more than N, and₁represents the 1 st sample value in Y, Y_nRepresenting the n-th sample value in Y, Y_NDenotes the Nth sample value in Y, the symbol "[ alpha ]]"is a symbol representing a vector or a matrix.

Step 2: designing respective corresponding Finite Impulse Response (FIR) band-pass filters according to respective frequency ranges (namely frequency bands) of a delta wave pitch signal, a theta wave pitch signal, an alpha wave pitch signal and a beta wave pitch signal in the electroencephalogram signal; then, respectively filtering Y by utilizing respective corresponding finite impulse response band-pass filters of a delta wave pitch signal, a theta wave pitch signal, an alpha wave pitch signal and a beta wave pitch signal in the electroencephalogram signals, and correspondingly extracting the delta wave pitch signal, the theta wave pitch signal, the alpha wave pitch signal and the beta wave pitch signal in the YSignal, represented correspondingly in vector form as Y_δ、Y_θ、Y_α、Y_β，Y_δ＝[y_δ,1…y_δ,n…y_δ,N]，Y_θ＝[y_θ,1…y_θ,n…y_θ,N]，Y_α＝[y_α,1…y_α,n…y_α,N]，Y_β＝[y_β,1…y_β,n…y_β,N](ii) a Wherein the frequency range of delta nodal signals in the electroencephalogram signals is 0.78-3.90 Hz, the frequency range of theta nodal signals in the electroencephalogram signals is 3.91-7.80 Hz, the frequency range of alpha nodal signals in the electroencephalogram signals is 7.81-13.28 Hz, the frequency range of beta nodal signals in the electroencephalogram signals is 13.29-30.47 Hz, and the frequency range of Y nodal signals in the electroencephalogram signals is 13.29-30.47 Hz_δ、Y_θ、Y_α、Y_βAll dimensions of (A) are 1 XN, Y_δIs a sample value, y_δ,1Represents Y_δ1 st sample value of (1), y_δ,nRepresents Y_δN-th sample value of (1), y_δ,NRepresents Y_δSample value of Nth, Y_θIs a sample value, y_θ,1Represents Y_θ1 st sample value of (1), y_θ,nRepresents Y_θN-th sample value of (1), y_θ,NRepresents Y_θSample value of Nth, Y_αIs a sample value, y_α,1Represents Y_α1 st sample value of (1), y_α,nRepresents Y_αN-th sample value of (1), y_α,NRepresents Y_αSample value of Nth, Y_βIs a sample value, y_β,1Represents Y_β1 st sample value of (1), y_β,nRepresents Y_βN-th sample value of (1), y_β,NRepresents Y_βThe nth sample value in (1).

In this embodiment, the specific process of step 2 is:

according to the frequency range of a delta-wave pitch signal in an electroencephalogram signal, designing a first finite impulse response band-pass filter, wherein the first finite impulse response band-pass filter is recorded as H1, the sampling frequency of H1 is F Hz, the order of H1 is 10, the upper edge frequency of a stop band of H1 is 0.7 Hz, the lower edge frequency of the stop band is 4 Hz, the upper edge frequency of a pass band is 0.8 Hz, and the lower edge frequency of the pass band is 3.5 Hz; and then, filtering Y by using H1, and extracting to obtain a delta wave law signal in Y.

Similarly, according to the frequency range of a theta-wave pitch signal in the electroencephalogram signal, designing a second finite impulse response band-pass filter, which is recorded as H2, wherein the sampling frequency of H2 is F Hz, the order of H2 is 10, the upper edge frequency of a stop band of H2 is 3.9 Hz, the lower edge frequency of the stop band is 8 Hz, the upper edge frequency of a pass band is 4 Hz, and the lower edge frequency of the pass band is 7.5 Hz; and then, filtering Y by using H2 to extract and obtain a theta node law signal in Y.

According to the frequency range of alpha-nodal signals in the electroencephalogram signals, a third finite impulse response band-pass filter is designed and recorded as H3, the sampling frequency of H3 is F Hz, the order of H3 is 10, the upper edge frequency of a stop band of H3 is 7.8 Hz, the lower edge frequency of the stop band is 13 Hz, the upper edge frequency of a pass band is 7.9 Hz, and the lower edge frequency of the pass band is 12 Hz; and then, filtering Y by using H3, and extracting to obtain an alpha wave law signal in Y.

Designing a fourth finite impulse response band-pass filter according to the frequency range of beta-nodal signals in the electroencephalogram signals, wherein the fourth finite impulse response band-pass filter is recorded as H4, the sampling frequency of H4 is F Hz, the order of H4 is 10, the upper edge frequency of a stop band of H4 is 13 Hz, the lower edge frequency of the stop band is 31 Hz, the upper edge frequency of a pass band is 14 Hz, and the lower edge frequency of the pass band is 30 Hz; and then, filtering Y by using H4 to extract and obtain a beta wave law signal in Y.

And step 3: for Y_δIdentifying ITD modal parameters, and extracting to obtain Y_δFour ITD mode parameters of, will Y_δThe vector formed by the sequential arrangement of the four ITD modal parameters is used as Y_δFeature vector of (1), denoted fea_δ(ii) a In the same way for Y_θIdentifying ITD modal parameters, and extracting to obtain Y_θFour ITD mode parameters of, will Y_θThe vector formed by the sequential arrangement of the four ITD modal parameters is used as Y_θFeature vector of (1), denoted fea_θ(ii) a In the same way for Y_αIdentifying ITD modal parameters, and extracting to obtain Y_αFour ITD mode parameters of, will Y_αThe vector formed by the sequential arrangement of the four ITD modal parameters is used as Y_αFeature vector of (1), denoted fea_α(ii) a In the same way for Y_βIdentifying ITD modal parameters, and extracting to obtain Y_βFour ITD mode parameters of, will Y_βThe vector formed by the sequential arrangement of the four ITD modal parameters is used as Y_βFeature vector of (1), denoted fea_β(ii) a Wherein fea_δ、fea_θ、fea_α、fea_βAre all 1 x 4.

In this embodiment, fea in step 3_δThe acquisition process comprises the following steps:

step 3_ 1: analyzing from the induction process of the motor imagery electroencephalogram signals, when motor imagery induction signal stimulation occurs, the brain starts to respond, at the moment, the total number of synaptic transmission signals of neurons of the brain is increased sharply, a corresponding kinetic equation shows strong nonlinearity, the electroencephalogram signals collected in the period of time are transient nonlinear signals, and the transient nonlinear signals have chaotic characteristics and are very difficult to classify and process; and then, after adapting to the evoked stimulus signals, the total number of synaptic transmission signals of the neurons of the brain is reduced in the inhibition stage, the number of the synaptic transmission signals is continuously reduced in the corresponding kinetic equation, the electroencephalogram signals acquired in the period are free response vibration signals, and the kinetic equation is degenerated into a weak nonlinear vibration equation. Thus, Y will be_δ＝[y_δ,1…y_δ,n…y_δ,N]The discrete signal is regarded as a free response vibration signal for Y_δ＝[y_δ,1…y_δ,n…y_δ,N]Carrying out normal sampling, d-delay sampling and 2 d-delay sampling in sequence to obtain a normal sampling free response data vector, a d-delay sampling free response data vector and a 2 d-delay sampling free response data vector in sequence, and correspondingly marking as X_nor、X_{del_1}、X_{del_2}，X_nor＝[y_δ,1…y_δ,n…y_δ,m]，X_{del_1}＝[y_δ,d+1…y_δ,d+n…y_δ,d+m]，X_{del_2}＝[y_δ,2d+1…y_δ,2d+n…y_δ,2d+m](ii) a Then X is put in_nor＝[y_δ,1…y_δ,n…y_δ,m]Re-expressed as X_nor＝[x_nor,1…x_nor,n…x_nor,m]Is mixing X_{del_1}＝[y_δ,d+1…y_δ,d+n…y_δ,d+m]Re-expressed as X_{del_1}＝[x_{del_1,1}…x_{del_1,n}…x_{del_1,m}]Is mixing X_{del_2}＝[y_δ,2d+1…y_δ,2d+n…y_δ,2d+m]Re-expressed as X_{del_2}＝[x_{del_2,1}…x_{del_2,n}…x_{del_2,m}](ii) a Wherein d is a positive integer, d represents a sampling interval,m is a positive integer, m represents X_nor、X_{del_1}、X_{del_2}M is more than or equal to 1 and less than or equal to N-2d, X_nor、X_{del_1}、X_{del_2}All dimensions of (a) are 1 xm, y_δ,m、y_δ,d+1、y_δ,d+n、y_δ,d+m、y_δ,2d+1、y_δ,2d+n、y_δ,2d+mCorresponds to and represents Y_δM-th, d + 1-th, d + n-th, d + m-th, 2d + 1-th, 2d + n-th, 2d + m-th sample values, x_nor,1、x_nor,n、x_nor,mCorresponds to X_nor1 st element, nth element, mth element, x_{del_1,1}、x_{del_1,n}、x_{del_1,m}Corresponds to X_{del_1}1 st element, nth element, mth element, x_{del_2,1}、x_{del_2,n}、x_{del_2,m}Corresponds to X_{del_2}The 1 st element, the nth element, the mth element in (1).

Step 3_ 2: according to X_norAnd X_{del_1}Constructing a first augmented matrix, denoted X_zg1，And according to X_{del_1}And X_{del_2}Constructing a second augmented matrix, denoted X_zg2，Then let A ═ X_zg2(X_zg1)^-1And using QR (orthogonal trigonometric) decomposition method to determine A ═ X_zg2(X_zg1)^-1Processing to obtain a characteristic value matrix of A, and recording the characteristic value matrix as delta; wherein, X_zg1And X_zg2All dimensions of (A) are 2X m, A represents a system feature matrix, (X)_zg1)^-1Is X_zg1The dimension of a is 2 × 2 and the dimension of Δ is 2 × 2.

Step 3_ 3: calculating Y from Delta_δHas a natural frequency and an attenuation coefficient, corresponding to f_δAnd λ_δ，Wherein the content of the first and second substances,to representThe imaginary part of (a) is,to representReal part of, Δ_1,1Line 1, column 1 elements representing Δ, t' at Y_δ＝[y_δ,1…y_δ,n…y_δ,N]A time interval of d sample values is sampled,

step 3_ 4: applying free response theory to X_norThe description is as follows: x_nor＝ψ_norE_nor，Wherein psi_norRepresents X_norMode vector of ψ_norDimension of (1X 2, E)_norFor the introduction of intermediate variables, E_norDimension of (d) is 2 × m, e represents a natural base, e is 2.17 …, μ is λ_δ+j2πf_δ，μ^*Is a conjugate of mu, mu^*＝-λ_δ-j2πf_δJ represents an imaginary unit of number, then solving X by using a least square method_nor＝ψ_norE_norTo obtain psi_nor，ψ_nor＝X_norE_nor ^T(E_norE_nor ^T)^-1Wherein the superscript "T" denotes the transpose of the vector or matrix, (E)_norE_nor ^T)^-1Is E_norE_nor ^TThe inverse of (1); recalculate psi_norIs the amplitude and phase of (1), the correspondence is expressed as xi_norAndξ_nor＝abs(ψ_nor(1))+abs(ψ_nor(2))，wherein abs () is the absolute value function, ψ_nor(1) To indicate psi_norElement # 1 of_nor(2) To indicate psi_norThe number 2 element of (a) is,to representThe real part of (2); finally xi is_norAs Y_δAmplitude of (1), re-noted as xi_δWill beAs Y_δIs newly recorded as

Step 3_ 5: f. of_δ、λ_δ、ξ_δ、Is Y_δFour ITD modal parameters of, will f_δ、λ_δ、ξ_δ、Sequentially arranging the constructed vectors as Y_δFeature vector fea of_δ，

In this embodiment, fea in step 3_θ、fea_α、fea_βIs acquired in the same manner as the procedure of step 3_1 to step 3_5, wherein f is_θ、λ_θ、ξ_θ、Corresponds to Y_θNatural frequency, attenuation coefficient, amplitude, phase, f_α、λ_α、ξ_α、Corresponds to Y_αNatural frequency, attenuation coefficient, amplitude, phase, f_β、λ_β、ξ_β、Corresponds to Y_βNatural frequency, attenuation coefficient, amplitude, phase.

And 4, step 4: will fea_δ、fea_θ、fea_α、fea_βThe vector formed by arranging the sequences is taken as a feature vector of Y and is denoted as fea, fea ═ fea_δ fea_θ fea_α fea_β](ii) a Wherein the dimension of fea is 1 × 16.

On the basis of the above steps 1 to 4, subsequent classification may be performed, such as: randomly acquiring electroencephalogram signals in the same channel for multiple times according to a given sampling frequency, acquiring the characteristic vectors of the electroencephalogram signals acquired in the same channel every time in the same mode according to the processes of the step 1 to the step 4, and forming a characteristic matrix by the characteristic vectors of the electroencephalogram signals acquired in the same channel for multiple times, wherein each row of vectors in the characteristic matrix represents the characteristic vector of the electroencephalogram signals acquired at one time; and then classifying the characteristic matrix to obtain a motor imagery type corresponding to the electroencephalogram signals randomly acquired each time.

For another example: acquiring electroencephalogram signals at a plurality of channels according to a given sampling frequency aiming at a plurality of different channels, acquiring the characteristic vector of the electroencephalogram signals acquired at each channel in the same mode according to the processes from step 1 to step 4, and arranging the characteristic vectors of the electroencephalogram signals acquired at the plurality of channels at the same time in sequence to form a row vector; after multiple times of acquisition, obtaining a plurality of row vectors, arranging the row vectors in sequence to form a characteristic matrix, wherein each row vector in the characteristic matrix represents a row vector formed by arranging the characteristic vectors of the electroencephalogram signals acquired at multiple channels at one time in sequence; and then classifying the characteristic matrix to obtain a motor imagery type corresponding to the electroencephalogram signal acquired each time.

In this embodiment, when a specific arrangement mode is not given, the sequence arrangement referred to in multiple places may be arranged according to the same self-defined sequence when the same kind of information is acquired, and may be arranged according to a sequence in general.

Let Y be [ 9.729.386.017.467.235.077.925.735.185.71 ]]F is 200, then Y_δ＝[0.7029 1.4001 1.8677 2.4437 3.0069 3.4138 4.0195 4.4626 4.8542 5.2688]，Y_θ＝[0.5803 1.2786 1.8847 2.6101 3.3521 3.9654 4.6916 5.2505 5.6846 6.0344]，Y_α＝[-0.0364 0.3748 1.1926 2.1650 3.2642 4.3742 5.2636 5.9892 6.3732 6.3498]，Y_β＝[-1.2230 -2.5862 -2.9156 -2.1651 -0.4555 1.6578 2.6811 2.5778 1.5536 0.1579](ii) a And when d is equal to 2, f_δ＝2.3755，λ_δ＝0.0828，ξ_δ＝7.8459，f_θ＝5.4966，λ_θ＝-14.5290，ξ_θ＝3.2526，f_α＝8.7462，λ_α＝-32.9419，ξ_α＝-1.9559，f_β＝21.3646，λ_β＝2.9961，ξ_β＝-3.1203，

In order to verify the feasibility and effectiveness of the method of the invention, experiments were carried out on the method of the invention.

The method is used for extracting the characteristics of the CLA left and right hand motor imagery data set in the current maximum Slow Cortical Potential (SCP) data set of the international standard electroencephalogram signal database, the extracted characteristic matrix is classified through the adaboost algorithm, the classification accuracy is as high as 92.7%, and compared with most characteristic extraction algorithms such as a CSP (Common Spatial Pattern) algorithm, a sparse self-coding method and the like, the accuracy is obviously improved, and the feasibility and the effectiveness of the method are fully demonstrated.