CN104899436A - Electroencephalogram signal time-frequency analysis method based on multi-scale radial basis function and improved particle swarm optimization algorithm - Google Patents

Abstract

The present invention discloses an electroencephalogram signal time-frequency analysis method based on a multi-scale radial basis function and an improved particle swarm optimization algorithm. The method performs time-frequency characteristic extraction analysis on an electroencephalogram signal with a time-varying parameter modeling method, and introduces the multi-scale radial basis function for identifying time-varying parameters in an expansion manner. A time-varying parameter model is created first; then the time-varying parameters of the model are identified with a basis function expansion method, that is to say, the time-varying parameters are represented as linear weighted combination of a set of multi-scale radial basis functions, and the identification problem of the time-varying parameters is converted into identification of time-invariant parameters; the optimal scale of the radial basis function is determined by the particle swarm optimization algorithm; and finally, time-frequency distribution characteristics of the electroencephalogram signal are calculated according to a time-varying parameter estimation value and a power spectral density formula. Compared to an current time-frequency analysis method, the method is capable of obtaining relatively high time and frequency resolution at the same time, and accurately extracting the time-frequency distribution characteristics of the electroencephalogram signal, and is of great significance on application of the epileptic electroencephalogram signal and auxiliary diagnosis of epilepsy.

Description

Based on the EEG signals Time-Frequency Analysis Method of multiple dimensioned radial basis function and Modified particle swarm optimization algorithm

Technical field

The present invention relates to a kind of EEG signals Time-Frequency Analysis Method, particularly relate to a kind of Time-Frequency Analysis Method based on multiple dimensioned radial basis function and Modified particle swarm optimization algorithm, belong to Digital Signal Analysis and Processing technical field.

Background technology

Epilepsy is by brain nervous cell supersynchronous electric discharge repeatedly; the nervous system disease of the spontaneity caused, paroxysmal functional disorders of brain; there are some problems in the Main Means-excision focus of current epilepsy therapy and drug therapy; complication and some bad reactions may be brought to some patient; if can make a definite diagnosis in early days before epileptic attack; take safeguard measure in advance; just can reduce the risk of patient injury greatly, and facilitation can be played to the pathogenesis of understanding epilepsy and the new methods for the treatment of of research.Electroencephalogram (Electroencephalogram, EEG) a kind ofly assesses the state of brain and the important form of disease.By EEG(electrocardiogram) examination and immediate analysis EEG signals, from the brain electric information of epileptic, extract the characteristic parameter that can reflect cerebral function state, in advance clinical intervention is carried out to these people.Therefore, EEG(electrocardiogram) examination extract EEG signals feature the diagnosis of epilepsy and research are had great significance.

EEG signals is a kind of non-stationary signal.Because the frequency of non-stationary signal changes in time, time frequency analysis need be carried out to it, namely describe the analytical approach that signal spectrum content distributes in time, extract its time-frequency characteristics.EEG signals can be divided into 4 wave bands by frequency: δ ripple (0.5-4Hz), θ ripple (4-8Hz), α ripple (8-13Hz), β ripple (13-30Hz).

Current Time-Frequency Analysis Method mainly contains two classes: nonparametric technique and parametric technique.Nonparametric technique such as Short Time Fourier Transform represents based on the nonparametric distributed to signal time-frequency combination, and according to Heisenberg uncertainty principle, the major defect of the method is that temporal resolution and frequency resolution can not reach optimum simultaneously.Parametric technique carries out modeling to signal and parameter is estimated, can obtain higher time frequency resolution simultaneously.

The main task of time-varying parameter model method carries out identification to its time-varying parameter.The main basis function development method adopted at present.Time-varying parameter is expressed as the linear weighted combination (Chen Yu of one group of known basis function by the main thought of basis function development method, Chen Huaihai, Li Zancheng, Deng. based on the time-varying parameter identification [J] of monetary multiplier and wavelet transformation. external electronic measurement technique, 2011,30 (7): 20-23.), by time become problem and be converted into time invariant parameter identification problem about basis function, by pair time invariant parameter identification so that obtain time-varying parameter.Current alternative basis function has Fourier's base, radial basis function (Li Y, Wei H L, Billings S A, etal.Time-varying model identification for time – frequency feature extraction from EEG data [J] .Journal of Neuroscience Methods, 2011,196:151 – 158.) etc.Often kind of basis function has respective approximation properties, as Fourier basis functions and Legendre polynomial can time-varying parameters that effectively identification change is slow and level and smooth, and radial basis function can identification simultaneously level and smooth and change violent time-varying parameter.Thus the present invention adopts radial basis function to carry out time-varying parameter identification.

EEG signals is a kind of faint bio signal, and because its noise is large, the features such as the low and non-stationary of energy is strong, are thus difficult to extract the potential temporal characteristics with biological nature.In this context, study a kind of Time-Frequency Analysis Method based on time-varying parameter modelling, and adopt radial basis function identification time-varying parameter can obtain higher time frequency resolution simultaneously, effectively can extract brain electricity time-frequency characteristics, the accurate time-frequency characteristics for epileptic EEG Signal extracts and assists the diagnosis of epileptic condition significant.

Summary of the invention

The invention provides a kind of brain based on time-varying parameter model electricity Time-Frequency Analysis Method, adopt multiple dimensioned radial basis function (Multi-scale radial basis function, MRBF) time-varying parameter is launched, then pair time-varying uncertainty problem is converted into time-invariant model, and time, invariant parameter carries out identification.Wherein, radial basis function has multiple dimensioned and multi-resolution characteristics, and to change, fast and slowly time-varying parameter can effective recognition and tracking, has been widely used in the time-varying parameter identification with multiple behavioral characteristics.The optimal scale of radial basis function is determined by Modified particle swarm optimization algorithm (Particle swarm optimization, PSO).PSO algorithm is a kind of Swarm Intelligence Algorithm, optimum solution (Li Xiuying can be sought in numerous particle according to enlightening information-fitness value, Han Zhigang. a kind of Nonlinear System Identification based on particle group optimizing [J]. control and decision-making, 2011,26 (11): 1627-1631.).The present invention, by carrying out time-frequency characteristics extraction and analysis to epileptic EEG Signal and normal brain activity electric signal, is accurately extracted the time-frequency characteristics of epileptic EEG Signal and normal brain activity electric signal, can provide quantitative technical support to the auxiliary diagnosis of epilepsy cerebral disease.

For achieving the above object, the invention provides:

Based on a non-stationary signal Time-Frequency Analysis Method for multiple dimensioned radial basis function and Modified particle swarm optimization algorithm, comprise the steps:

1. time-varying parameter modelling;

2. time-varying parameter launches, and time-varying model is converted into time-invarying parameter model;

3. select center and the yardstick of radial basis function;

4. time, invariant parameter is estimated;

5. time-varying uncertainty;

6. time frequency analysis, asks signal time-frequency distributions feature by time-varying uncertainty value.

Wherein, in described step 3, by the optimal scale of Modified particle swarm optimization algorithms selection radial basis function.

In described step 3, the inertia weight of Modified particle swarm optimization algorithm medium velocity more new formula is the changing value relevant with fitness value.

In described step 3, weigh the fitness value of the index of correlation as particle of models fitting effect.

The advantage of the non-stationary signal Time-Frequency Analysis Method based on multiple dimensioned radial basis function and particle swarm optimization algorithm provided by the present invention comprises:

1. time frequency resolution is high, can accurate recognition very noisy, non-stationary, nonlinear signal, accurately extracts its time-frequency characteristics;

2. can select radial basis function yardstick flexibly according to signal, adaptability is good;

3. accurately extract epileptic EEG Signal and normal brain activity electric signal time-frequency characteristics, the auxiliary diagnosis for epilepsy cerebral disease provides quantitative technological guidance.

Accompanying drawing explanation

Fig. 1 is the schematic flow sheet of the Time-Frequency Analysis Method according to the embodiment of the present invention.

Fig. 2 is the schematic flow sheet that particle group optimizing method selects multiple dimensioned radial basis function optimal scale.

Fig. 3 (a) shows true epileptic EEG Signal.Fig. 3 (b)-3 (c) shows according to the Time-Frequency Analysis Method of the embodiment of the present invention and classical Time-Frequency Analysis Method true epileptic EEG Signal time frequency analysis Comparative result; Wherein Fig. 3 (b) is the time frequency analysis result of fourier methods in short-term, the time frequency analysis result that Fig. 3 (c) is the inventive method.

Fig. 4 (a) shows true normal brain electric signal.Fig. 4 (b)-4 (c) shows according to the Time-Frequency Analysis Method of the embodiment of the present invention and classical Time-Frequency Analysis Method the time frequency analysis Comparative result of true normal brain electric signal; Wherein Fig. 4 (b) is the time frequency analysis result of fourier methods in short-term, the time frequency analysis result that Fig. 4 (c) is the inventive method.

Embodiment

Below in conjunction with the drawings and specific embodiments, the present invention is described in further detail.

The time frequency analysis of true EEG signals verifies time frequency analysis effect of the present invention as embodiments of the invention, and the present invention and classical Time-Frequency Analysis Method is contrasted.The time frequency analysis of true EEG signals comprises and comparing the time frequency analysis of epileptic EEG Signal and normal brain electric signal.

Experimental data is from the eeg data storehouse of Bonn, Germany epilepsy research department clinical acquisitions, and sample frequency is 173.6Hz, and the duration is 23.6 seconds.As shown in Fig. 3 (a) He Fig. 4 (a), be respectively epileptic EEG Signal and normal brain electric signal.

Fig. 1 illustrates the process flow diagram according to embodiment of the present invention method, comprising:

First time-varying parameter modelling is carried out to signal, carried out determining rank (step 1) to model by final predicated error (Final prediction error, FPE) model selection criteria; Then time-varying parameter MRBF is launched, invariant parameter form (step 2) when being expressed as; Then select MRBF center and by the optimal scale (step 3) of PSO algorithms selection MRBF; Time again pair, invariant parameter is estimated, is matrix form by model representation, pair time invariant parameter carry out identification (step 4), obtain time-varying uncertainty value; By time invariant parameter estimated value and MRBF estimate time-varying parameter (step 5); Finally carry out time frequency analysis, calculate time-frequency spectrum by time-varying parameter, obtain time frequency distribution map (step 6).

Lower mask body introduction is according to the concrete grammar of the Time-Frequency Analysis Method based on multiple dimensioned radial basis function and Modified particle swarm optimization algorithm provided by the present invention:

1. time-varying parameter modelling

Time-varying parameter model parameter is time dependent, and a p rank time-varying parameter model is as follows:

$y\left(t\right)=\underset{i=1}{\overset{p}{\mathrm{\Σ}}}{a}_i\left(t\right)y(t-i)+e\left(t\right)---\left(1\right)$

In formula: p is model order, a _i(n) (i=1,2 ..., p) be time-varying parameter, e (n) is average is 0, and variance is σ ²white Gaussian noise.

Rank determined by model: the present invention adopts final predicated error (FPE) order selection criteria preference pattern exponent number, and its expression formula is as follows:

In formula, N is data length, and p is model order, for predicated error.

2. time-varying parameter launches, and time-varying parameter is expressed as the weighted linear combination of one group of multiple dimensioned radial basis function, is shown below:

$a_i\left(t\right)=\underset{m=1}{\overset{M}{\mathrm{\Σ}}}{c}_{i,m}{f}_m\left(t\right)---\left(3\right)$

In formula, c _i,mfor the time constant weight coefficient of basis function, f _mt () is Gaussian radial basis function.M is the dimension of basis function.Above formula is substituted into formula (1) to obtain:

$y\left(t\right)=\underset{i=1}{\overset{p}{\mathrm{\Σ}}}\underset{m=1}{\overset{M}{\mathrm{\Σ}}}{c}_{i,m}{f}_m\left(t\right)y(t-i)+e\left(t\right)---\left(4\right)$

Formula (4) is time-invarying parameter model.

3. the determination of radial basis function center and yardstick

The present invention adopts Gaussian radial basis function to expand time-varying parameter as basis function.Due to Gaussian radial basis function have radially symmetrical, approximation capability is strong, effective, the advantages such as pace of learning is fast, and less by the impact of signal characteristic, so the present invention selects Gaussian radial basis function to launch time-varying parameter as basis function, its expression formula is as follows:

$\mathrm{\φ}\left(\right||x-c|\left|\right)=\exp [-\frac{{\left|\right|x-c\left|\right|}^2}{2{\mathrm{\σ}}^2}]---\left(5\right)$

In formula, c is the center of radial basis function, σ ²for the yardstick of basis function, which determine the distance of basis function around central point.X is input amendment.When radial basis function is as basis function, x (x=1,2 ..., N, N are sampling length) and be discrete sampling time series.

Suitable radial basis function center and yardstick is selected to be the keys of accurate recognition time-varying parameter.Be distributed in whole time-varying parameter to make radial basis function, to guarantee that all local of radial basis function to time-varying parameter are accurately estimated, the center of radial basis function is evenly distributed in time-varying parameter by the present invention, if the center of a kth radial basis function and radial basis function is expressed as:

${\mathrm{\φ}}_k\left(\right||x-c_k|\left|\right)=\exp [-\frac{{\left|\right|x-c_k\left|\right|}^2}{2{\mathrm{\σ}}_k^2}],k=\mathrm{1,2},...,M---\left(6\right)$

$c_k=\frac{k\×N}{M}---\left(7\right)$

In formula, M is the dimension of radial basis function, and x is sample input, c _kfor the center of a kth radial basis function, for the yardstick of a kth radial basis function, which determine the distance of this basis function around central point, N is sampled data length.

In the present invention, the optimal scale of radial basis function is calculated by Modified particle swarm optimization algorithm.

Fig. 2 is the optimal scale process flow diagram that Modified particle swarm optimization algorithm calculates radial basis function, and concrete steps are as follows:

(a) initialization: initialization particle, particle number, fitness value, local optimum particle, global optimum's particle, particle maximal value and minimum value, speed maximal value, iterations etc.

Candidate's yardstick of radial basis function is:

${\mathrm{\σ}}_k^2=\frac{N^2\×2^{-S_k}}{M}---\left(8\right)$

In formula, N is sampled data length, and M is the dimension of radial basis function, s _kfor needing the arbitrary integer carrying out regulating, and maximal value is within 10, and particle is s _kthe vector formed is u _i=[s ₁, s ₂..., s _m].

Getting random integers is particle initialize, and initial fitness value is 0, and population and iterations can adjust according to characteristic change parameter.

B () calculates the fitness value of each particle

The fitness value of particle is the enlightening information helping to find optimal particle, and in the present invention, adopt the index of correlation of models fitting effect as fitness function, index of correlation formula is as follows:

$R^2=1-\frac{\underset{t=1}{\overset{N}{\mathrm{\Σ}}}{(Y\left(t\right)-\hat{Y}\left(t\right))}^2}{\underset{t=1}{\overset{N}{\mathrm{\Σ}}}{(Y\left(t\right)-\stackrel{\‾}{Y})}^2}---\left(9\right)$

In formula, R ²for the index of correlation, N is sampled data length, for the model prediction of data Y exports.The index of correlation larger explanation models fitting effect is better.

Particle is substituted into yardstick formula (8), obtain the yardstick of radial basis function, then according to radial basis function center, obtain multiple dimensioned radial basis function.Again by multiple dimensioned radial basis function method of deploying to the identification of Model Parameters of time-varying system, obtain the model predication value of measurement data Y finally calculate the index of correlation according to formula (9), obtain the fitness value of particle.

C () upgrades local optimum particle.Local optimum particle is the particle that in previous cycle, adaptive value is maximum, finds the particle that fitness value is maximum, and this particle is as local optimum particle.

D () upgrades global optimum's particle.Global optimum's particle is the particle that in all circulations, fitness value is maximum.If local optimum particle fitness value is larger than global optimum particle, then assignment is to global optimum's particle.Otherwise global optimum's particle is constant.

E () upgrades the speed of each particle, by particle rapidity more new formula upgrade the speed of each particle.

Particle rapidity more new formula is as follows:

$v_i^{(h+1)}=\mathrm{\ω}\×v_i^{\left(h\right)}+\mathrm{\θ}\×c_1\×(C_{\mathrm{best}}^{\left(h\right)}-u_i^{\left(h\right)})+\mathrm{\θ}\×c_2\×(g_{\mathrm{best}}^{\left(h\right)}-u_i^{\left(h\right)})---\left(10\right)$

Wherein, C _bestfor local optimum particle, g _bestfor global optimum's particle, be i-th particle of the h time circulation, be the speed of i-th particle of the h time circulation, be the local optimum particle of the h time circulation, it is current global optimum particle after the h time circulation.θ is the random number between 0 ~ 1, c ₁and c ₂effect be prevent particle to be absorbed in local optimum:

$c_1=2.5-\frac{2\×h}{1\×H}---\left(11\right)$

$c_2=0.5-\frac{2\×h}{1\×H}---\left(12\right)$

In formula, H is maximum cycle, h=1,2 ..., H.

Particle rapidity is more in new formula, and ω is inertia weight.In conventional particle colony optimization algorithm, ω is generally constant, is unfavorable for finding optimal particle flexibly.In the present invention, adopt a kind of Modified particle swarm optimization algorithm, make inertia weight value be associated with fitness value, to find optimal particle more flexibly, its expression formula is as follows:

$\mathrm{\ω}=1-\mathrm{\θ}\×\frac{F_i-F_{\mathrm{best}}}{F_{\mathrm{best}}}---\left(13\right)$

In formula, θ is the random number between 0 ~ 1, F _ithe fitness value of i-th particle, F _bestit is the fitness value of global optimum's particle.

When speed close to 0 time, be particle rapidity assignment again by following formula:

$v_i^{[h+1]}{|}_j=\±\mathrm{\θ}\×\mathrm{\γ}\×V_{\max }---\left(14\right)$

In formula, be in a jth element, the γ random number that to be constant value 0.1, θ be between 0 ~ 1.

F () upgrades the position of each particle.Particle position more new formula is as follows:

$u_i^{(h+1)}=u_i^{\left(h\right)}+v_i^{(h+1)}---\left(15\right)$

G () returns step (b), repeated execution of steps (b), (c), (d), (e) and (f), until obtain optimal particle.Finally obtain global optimum's particle and be optimal particle, obtain the optimal scale of radial basis function further according to yardstick formula.

4. time, invariant parameter is estimated.By least square method to the time invariant parameter c in time-invarying parameter model _i,mestimate, first time-invariant model is written as matrix form.

Parameter in formula (4) is defined as follows:

f(t)＝[f ₁(t),f ₂(t),…,f _M(t)],

H _i(t)＝y(t-i)f(t),

H(t)＝[H ₁(t),H ₂(t),…,H _p(t)],

C _i＝[c _i,1,c _i,2,…,c _i,M],

C＝[C ₁,C ₂,…,C _p],

Then write formula (4) as matrix form:

y(t)＝H(t)C ^T+e(t) (16)

In formula, T is transpose of a matrix.C is tried to achieve by least-squares algorithm, namely

In formula, H=[H (1), H (2) ..., H (N) ,] Y=[y (1), y (2) ..., y (N)].

5. time-varying uncertainty.When obtaining after invariant parameter estimated value, by time invariant parameter and multiple dimensioned radial basis function substitute into formula (3) estimated value of time-varying parameter can be obtained.

6. time frequency analysis.According to the time-varying parameter estimated, applied power spectrum formula carries out time frequency analysis to signal, and its computing formula is as follows:

$S(f,t)=\frac{{\widehat{\mathrm{\σ}}}_e^2}{{|1-{\mathrm{\Σ}}_{i=1}^p{\hat{a}}_i{e}^{-j2\mathrm{\πif}/f_s}|}^2},---\left(18\right)$

In formula, f _sfor sample frequency, for the estimated value of time-varying parameter, for the variance of observational error.

Obtain the time frequency distribution map of signal, namely obtain the time-frequency distributions feature of signal.

The inventive method and basis function development method and classical Time-Frequency Analysis Method and Short Time Fourier Transform method compare, and analyze, disclose the instantaneous time varying characteristic that EEG signals is potential to true epileptic EEG Signal and normal brain activity electric signal time-frequency result.Time frequency analysis result is respectively as shown in Fig. 3 (b)-3 (c) He Fig. 4 (b)-4 (c).

From epileptic EEG Signal and normal brain activity electric signal time frequency analysis result, two kinds of methods obtain more consistent frequecy characteristic distribution, but obviously the time frequency resolution of Short Time Fourier Transform method is low, result is as shown in Fig. 3 (b) He Fig. 4 (b), and the distribution of time-frequency figure medium frequency is unintelligible.And the time frequency resolution of the inventive method is high, result is as shown in Fig. 3 (c) He Fig. 4 (c).

The time-frequency distributions feature of epileptic EEG Signal and normal brain activity electric signal clearly can be obtained by the time-frequency figure of the inventive method.Especially, in the time frequency distribution map of epileptic EEG Signal, clearly can observe two kinds of frequency range ripples: θ ripple and β ripple.As shown in Fig. 3 (c), θ ripple (4-8Hz) occurs respectively between 0 ~ 6 second, between 9 ~ 11 seconds and 14 ~ 22 seconds; β ripple (13-30Hz) appears at the time period between 1 ~ 4 second and 12 ~ 22 seconds.Similarly, in the time frequency distribution map of normal brain activity electric signal, the brain wave of two kinds of different frequency ranges can be observed: δ ripple and α ripple.As shown in Fig. 4 (c), the time period of δ ripple (0-4Hz) between 4 ~ 5 seconds, 10 ~ 12 seconds, 14 ~ 16 seconds and 18 ~ 22 seconds occurs, α ripple (8-12Hz) mainly appears at the time period between 17 ~ 20 seconds and 14 ~ 15 seconds.

From the time frequency analysis result of true EEG signals, compared with classical Short Time Fourier Transform Time-frequency method, the inventive method accurately can obtain the time-frequency distributions feature of true EEG signals.From epileptic EEG Signal and normal brain activity electric signal time frequency analysis Comparative result, normal brain activity electric signal on the whole time period, the electrical activity of brain that occurrence frequency is lower, and the brain wave frequency motion frequency of epileptic EEG Signal obviously increases.And epileptic's event of showing effect all is detected, the outbreak event be detected, although the outbreak duration very short (2s), the inventive method still can detect, and the brain electrical feature extracted has obvious change in the epileptic attack phase.The inventive method accurately detects in the EEG signals of epileptic attack to there is the brain wave with normal brain activity electric signal different frequency, and providing the instantaneous frequency distributed intelligence of different frequency, these information can provide quantitative analysis tool for doctor to the auxiliary diagnosis that epilepsy cerebral disease is shown effect.

EEG signals Time-Frequency Analysis Method based on multiple dimensioned radial basis function and Modified particle swarm optimization algorithm provided by the present invention is mainly accurately extract the instantaneous time-frequency distributions feature of EEG signals, and the diagnosis of auxiliary epilepsy cerebral disease proposes.But obviously, the Time-Frequency Analysis Method described in this instructions is also applicable to the time frequency analysis of other non-stationary signals.

Above the EEG signals Time-Frequency Analysis Method based on multiple dimensioned radial basis function and Modified particle swarm optimization algorithm provided by the present invention has been described in detail, but obvious scope of the present invention is not limited thereto.When not departing from the protection domain that appended claims limits, to the various changes of above-described embodiment all within the scope of the present invention.

Claims (4)

1., based on an EEG signals Time-Frequency Analysis Method for multiple dimensioned radial basis function and Modified particle swarm optimization algorithm, it is characterized in that comprising:

Step 1. time-varying parameter modelling;

Step 2. time-varying parameter launches, and time-varying model is converted into time-invarying parameter model;

Step 3. selects center and the yardstick of radial basis function;

During step 4., invariant parameter is estimated;

Step 5. time-varying uncertainty;

Step 6. time frequency analysis, asks signal time-frequency distributions feature by time-varying uncertainty value.

2., as claimed in claim 1 based on the EEG signals Time-Frequency Analysis Method of multiple dimensioned radial basis function and Modified particle swarm optimization algorithm, it is characterized in that:

Described step 3 comprises, by the optimal scale of Modified particle swarm optimization algorithms selection radial basis function.

3., as claimed in claim 2 based on the EEG signals Time-Frequency Analysis Method of multiple dimensioned radial basis function and Modified particle swarm optimization algorithm, it is characterized in that:

Described step 3 comprises, and the inertia weight of Modified particle swarm optimization algorithm medium velocity more new formula is the changing value relevant with fitness value.

4., as claimed in claim 2 based on the EEG signals Time-Frequency Analysis Method of multiple dimensioned radial basis function and Modified particle swarm optimization algorithm, it is characterized in that:

Described step 3 comprises, and weighs the fitness value of the index of correlation as particle of models fitting effect.