CN101869477A - Self-adaptive EEG signal ocular artifact automatic removal method - Google Patents

Abstract

The invention provides a self-adaptive EEG signal ocular artifact automatic removal method, which belongs to the technical field of biological information and is mainly used in a pretreatment process for acquiring an EEG signal. The method comprises: performing real-time empirical mode decomposition (EMD) of collected EEG data having ocular artifacts; performing Hilbert transform of all obtained mode components to obtain a instantaneous frequency; according to the time-frequency property of the ocular artifacts in the EEG signal and the statistical property of the empirical mode components, performing the threshold filtering of all obtained mode components; and performing data reconstruction by using all mode components obtained after filtration. The method solves the manual screening problem of the empirical mode components having the ocular artifacts, thereby automatically removing the ocular artifacts from the EEG signal.

Description

The automatic removal method of the electric artefact of eye in a kind of self adaptation EEG signals

Technical field:

The present invention relates to the biology information technology field, particularly EEG signals (Electroencephalography, EEG) collection and preconditioning technique.Specifically relate to based on Hilbert-Huang transform (Hilbert-Huang Transform, HHT) the electric artefact of eye (Electrooculography, automatic removal technology EOG) in the EEG signals.

Background technology:

The analysis of EEG signals and research are being brought into play important effect in fields such as neuroscience, psychology, biomedicines.Yet the EEG signal is a kind of very faint electricity physiological signal, is vulnerable to the influence in various types of interference source in gatherer process.The brain electrical interference has caused very big difficulty for the analysis EEG signals, has especially hindered the computer automatic analysis and the diagnosis of EEG signals.Therefore, how to detect and eliminate the various interference components among the EEG effectively, extracting true and reliable EEG information, the analysis and the research of EEG signals is played crucial effects.

Artefact in the EEG signals mainly contains the electric artefact of eye, tongue electricity artefact, myoelectricity artefact, pulse artefact and perspiration artefact etc.With respect to other interference signal sources, the electric artefact of eye is a kind of topmost interference component that is present in the EEG signal.It results from human body self, when nictation or ocular movement, can cause bigger potential change.When gathering EEG, the electric artefact of eye sends from its source, permeates into whole scalp, its amplitude can reach 100mV, and the EEG signals amplitude is very faint, and general scalp EEG signals has only about 50 μ V, the EEG signal that therefore will cause collecting produces obviously distortion, forms the electric artefact of eye.And the frequency band of the electric artefact of eye has also covered the frequency band of EEG signals, is difficult to by the method for filtering it be removed.In experiment, the experimenter blinks, eye movement occurs at random, and avoids nictation and eye movement for a long time, is difficult to accomplish.Therefore, the research of eye electrical interference removing method is the important content in the EEG signals pretreatment always.

Along with the continuous development of signal processing technology, the method that artefact is removed has also had considerable improvement.Early stage method mainly is to control by experiment to reduce electric artefact of eye and artefact subtraction (Prtifact Subtraction, PS) remove artefact, principal component analysis (Principle Component Analysis had appearred again afterwards, PCA), independent component analysis (Independent Component Analysis, ICA) and wavelet transformation (WaveletTransform, new method such as WT).

So-called experiment control method is meant to require the subjects to close one's eyes in experiment or avoid nictation and eye movement as far as possible, reaches the purpose of the electric artefact of restriction eye.This method is difficult to make it to be used widely to subjects's harsh requirement.

The artefact subtraction is a kind of method of using early, and its feature is visual and understandable, explicit physical meaning.The basic assumption of this method is that the observation EEG signal that measures is the linear combination of real source EEG signal and artefact, and true EEG source is uncorrelated with artefact, and artefact can record by measurement means.

y _t(i)＝y(i)-kx(i)....................................(1)

(1) in the formula, y _t(i) EEG behind the artefact is removed in expression; The EEG that y (i) expression measures; The electric artefact of x (i) expression eye; K represents proportionality coefficient, and the EEG after promptly correcting removes a certain proportion of artefact from the EEG signal of measuring.The regressive method of the permanent really usefulness of k value is estimated.But all depending on, regressive method sets up a correct recurrence lead (EOG), and the activated diffusion of electro-ocular signal and EEG signals all has the amphicheirality, be the composition that has also comprised EEG among the actual EOG of recording, so homing method is removed some EEG signal when removing artefact inevitably.

(Principle Component Analysis is a kind of multiple tracks signal processing method PCA) to PCA, is a kind of common method that PARAMETERS IN THE LINEAR MODEL is estimated.Basic thought is to utilize orthogonality principle that original relevant independent variable is transformed to the separate variable of another group.This method can effectively be found out in the data the element and the structure of " mainly ", removes noise and redundancy, with original complex data dimensionality reduction, discloses and is hidden in complex data simple structure behind.PCA is separate composition respectively leading signal decomposition on EEG respectively leads the basis that distributes, and removes unwanted composition, and reconstruct EEG again is to reach the purpose of removing artefact.Studies show that PCA significantly is better than based on regressive artefact subtraction on effect, yet PCA can not separate from EEG fully and the artefact of the current potential of its waveform similarity, effect is not as independent component analysis.Because PCA comes decomposed signal by the quadrature principle, the independence that does not have consideration respectively to decompose component.

(Independent Component Analysis is to separate (Blind Source Separation, BSS) the next multiple tracks signal processing method of technical development by blind source signal ICA) to independent component analysis.Its ultimate principle is independently to pass the multiple tracks observation signal in principle the optimization algorithm according to statistics to be decomposed into some independent elements, thereby realizes the enhancing and the analysis of signal.The thinking of ICA comes from central limit theorem: a class mean and variance are the stochastic variable of the same order of magnitude, and coefficient result must be near Gauss distribution.Therefore the non-Gauss of the separating resulting of one group of mixed signal that the mutual statistical independent information source is produced through the linearity combination measures, and when its non-Gauss reaches maximum, can think that mixed signal has realized optimal separation.The basic thought of ICA can be described below: establish time series X (t)=[x ₁(t), x ₂(t) ... x _n(t)] be the observation signal of n dimension, time series S (t)=[s ₁(t), s ₂(t) ... s _m(t)] be m the mutual statistical independent source signal of observation signal X (t), and observation signal X (t) is that source signal S (t) mixes and generation, i.e. X (t)=AS (t) through a unknown matrix A linearity.Under mixed matrix A and source signal S (t) condition of unknown, only utilize independently hypothesis of observation signal X (t) and source signal statistics, seek a linear transformation separation matrix W, wish that output signal U (t)=WX (t)=WAS (t) approaches real source signal S (t) as far as possible.

Think that in theory the interfering signal that electrocardio artefact, the electric artefact of eye, myoelectricity artefact and other interference sources in the EEG signals produced all is to be produced by separate information source.And EEG signals, electromyographic signal belong to inferior gaussian signal, and electrocardiosignal, electro-ocular signal belong to this signal of superelevation.Do not belong to gaussian signal, thereby the ICA method is suitable for the EEG Signal Processing.According to the criterion difference of independent degree between each component of tolerance, the ICA method has various ways, minimum as mutual information, information is very big, negentropy is maximum, based on nuclear, maximum likelihood estimation etc.In actual applications, the ICA method does not need the artefact reference electrode, can be used for the separation of various artefacts, and it is higher to remove the artefact precision, has shown its superiority.But with the ICA algorithm observation signal that contains Gaussian noise is directly carried out the independent variable analysis, produce bigger error sometimes.In addition, the order of each isolated component is a random alignment in the ICA decomposition result, need carry out manual screening, how to use ICA to carry out automatic artefact and remove, and also be an open question.

Wavelet transformation (Wavelet Transform) is a breakthrough development of Fourier transform, enjoys the good reputation of " school microscop ".Compare with Fourier transform, wavelet transformation has good time-frequency characteristic, has higher frequency resolution and lower temporal resolution in low frequency part, has higher temporal resolution and lower frequency resolution at HFS.This specific character just makes wavelet transformation have adaptivity to signal, and the processing of non-stationary signal is had good effect with analyzing.

But the wavelet transformation denoising requires the frequency band of signal and noise can not aliasing.But the frequency band of EEG and myoelectricity artefact, pulse and the electric artefact of eye is aliasing mutually, so researcher begins the wavelet transformation of classics is combined with existing denoising method.Classical ICA model is not considered the existence of noise, therefore with the ICA algorithm observation signal that contains Gaussian noise is directly carried out the independent variable analysis, produces bigger error sometimes.Some researcheres are attempted combining with wavelet analysis and ICA, with the signal to noise ratio of the soft threshold method raising of small echo brain electricity, utilize ICA to isolate source signal again, remove a noise and the electric artefact of eye in the brain electricity effectively.But the order of each isolated component is a random alignment in the ICA decomposition result, need carry out manual screening to the composition after decomposing, and how to use ICA to carry out automatic artefact and removes, and still fails to solve well.In addition, when carrying out wavelet analysis, the selection of wavelet basis is an important problems, but does not also have good method at present, mainly is to judge the quality of wavelet basis by the result with wavelet analysis with theoretical error, and selectes wavelet basis thus.

In order to study transient state and non-stationary phenomenon, people such as U.S. scientist Huang of Chinese origin had proposed a kind of brand-new non-linear and non-stationary signal processing method in 1998---and Hilbert-Huang transform (Hilbert-Huang Transform, HHT).

(Hilbert-Huang Transform HHT), is a kind of brand-new analytical method to Hilbert-Huang transform, and it decomposes (EMD) by empirical modal and Hilbert conversion two parts are formed, and its core is EMD.The HHT method is a kind ofly to have more adaptive time-frequency localization analytical method than Fourier transformation and wavelet transformation etc., is more suitable for being used to study non-stationary signal.

In order to study transient state and non-stationary phenomenon, frequency must be the function of time.Thus, draw the instantaneous frequency definition of (Instantaneous Frequency is called for short IF).But the notion of instantaneous frequency is disputable always, and most of viewpoints are thought that it does not exist or only thinks and existed under given conditions.And the non-stationary signal that definition of overall importance changes constantly for frequency in order to obtain significant instantaneous frequency, must be revised as partial restrictive condition without any meaning based on restrictive condition of overall importance.People such as Norden E Huang have proposed physically to define the essential condition of a significant instantaneous frequency IF: function is symmetrical in local zero-mean; And the number of zero crossing and extreme point equates.On these condition bases, people such as Norden E Huang have proposed the notion of intrinsic mode functions (Intrinsic Mode Function IMF).An intrinsic mode functions must satisfy following two conditions: (1) at whole data length, the number of extreme value and zero crossing must equate or differ from one at the most; (2) in any data point, on average being necessary for of the envelope of local maximum and the envelope of local minimum is zero, and promptly signal is about the local symmetry of time shaft.First qualifications is very tangible, is similar to the distribution of traditional stationary Gaussian process.Second condition is the place of innovation, and it becomes traditional qualification of overall importance the qualification of locality.This qualification can be removed the fluctuation of the instantaneous frequency that causes owing to waveform is asymmetric.Satisfy the mode component of above two conditions, have only an extreme point between its continuous two zero crossings, promptly include only the vibration of a basic model, do not have complicated stack ripple to exist.The basic model component is not restricted to narrow band signal thus defined, can be the non-stationary signal with certain broadband.

There is instantaneous frequency in IMF, and it can be tried to achieve by the Hilbert conversion.But general signal usually is a sophisticated signal, does not satisfy the IMF condition, also can't obtain instantaneous frequency.So people's creativeness such as Norden E Huang propose to suppose: any sophisticated signal all is made up of many IMF components; Each IMF mutual superposition just forms composite signal.People such as Norden E Huang have proposed Empirical mode decomposition (EMD) based on this hypothesis again

$X\left(t\right)=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{c}_i\left(t\right)+r\left(t\right)...\left(2\right)$

Wherein, n is the number of IMF component.This decomposition algorithm mainly acts on: the one, remove the stack ripple, and the 2nd, make waveform symmetry more.

Utilize the EMD method sophisticated signal can be resolved into the combination of a plurality of IMF, these IMF are carried out the Hilbert conversion, can obtain the instantaneous spectrum of each IMF component.The Hilbert conversion is a kind of linear transformation, and it emphasizes local property, and the instantaneous frequency that is obtained by it is best definition.To random time sequence X (t), the Hilbert conversion Y (t) that can obtain it is:

$Y\left(t\right)=\frac{1}{\mathrm{\π}}\underset{-\∞}{\overset{+\∞}{\∫}}\frac{X\left(\mathrm{\τ}\right)}{t-\mathrm{\τ}}\mathrm{d\τ}...\left(3\right)$

Obtain analytic signal Z (t) by X (t) and Y (t):

Z(t)＝X(t)+iY(t)＝a(t)e ^jθ(t)............................(4)

$a\left(t\right)=\sqrt{X^2\left(t\right)+Y^2\left(t\right)}...\left(5\right)$

$\mathrm{\θ}\left(t\right)={\mathrm{tg}}^{-1}\frac{Y\left(t\right)}{X\left(t\right)}...\left(6\right)$

(5) and (6) formula, clearly express instantaneous amplitude and instantaneous phase, reflected the temporal properties of data well.On this basis of (6) formula, defined instantaneous frequency and be:

$f\left(t\right)=\frac{1}{2\mathrm{\π}}\·\frac{\mathrm{d\θ}\left(t\right)}{\mathrm{dt}}...\left(7\right)$

(7) formula shows, instantaneous frequency is the function of time, and it has disclosed the tolerance of a certain moment signal energy in the frequency intensity.Making that the HHT conversion is more enough successfully is applied to processing non-linear, non-stationary signal.

Summary of the invention

The present invention is based on the electric artefact method of removaling of eye in the EEG signals of HHT, and its feature mainly comprises: according to the time-frequency characteristic of in the EEG signals electric artefact, and in conjunction with the statistical property of empirical modal component, the whole modal components that obtain are carried out threshold filter; Utilize filtered whole modal components to carry out data reconstruction.The invention solves manual screening problem, thereby reach the purpose of from EEG signals, removing the electric artefact of eye automatically the empirical modal component that comprises the electric artefact of eye.

The automatic removal method of the electric artefact of eye is characterized in that in a kind of self adaptation EEG signals provided by the invention, may further comprise the steps:

1) empirical modal of EEG signals decomposes

EEG signals X (t) is a non-stationary signal, carries out empirical modal and decomposes EMD, obtains each IMF component c _i(t) and remainder r (t);

Wherein, i is for decomposing the number of the IMF component that obtains, and the value of i is determined by X (t) fully;

2) EEG signals EMD is decomposed all IMF component c of gained _i(t) carry out Hilbert transform, try to achieve instantaneous frequency f _i(t),, utilize following formula to carry out bandpass filtering, extract the signal g in the EEG signals frequency band based on instantaneous frequency _i(t);

$g_i\left(t\right)=\left\{\begin{array}{ccc}0& \mathrm{if}& f_i\left(t\right)>f_h\\ 0& \mathrm{if}& f_i\left(t\right)<f_l\\ c_i\left(t\right)& \mathrm{if}& \mathrm{other}\end{array}\right.$ (I IV VI VII VIII); Here f _lBe 1Hz, f _hBe 30Hz;

3) utilize formula threshold function table formula (II) and threshold filter function formula (IV), carry out threshold filter, remove electro-ocular signal, obtain signal w _i(t);

τ _i＝mean|M _i-m _i|+std|M _i-m _i|.............................(II)

Wherein, M _iBe i IMF component g _iThe value of each extreme point (t); Mean represents to average, and std represents to ask standard deviation; m _iBe i IMF component g _i(t) time average;

Suppose time series g _i(t) length is N, then

$m_i=\frac{1}{N}\underset{t=1}{\overset{N}{\mathrm{\&Sigma;}}}{g}_i\left(t\right)...$ (Ⅲ)

Threshold function table τ _i, characterized g _i(t) signal is with average m _iThe amplitude upper limit for datum line; On this basis, at the big characteristics of eye electricity artefact amplitude, designed threshold filter function formula (IV);

$w_i\left(t\right)=\left\{\begin{array}{ccc}{m}_i& \mathrm{if}& |g_i\left(t\right)-m_i|>{\mathrm{\&tau;}}_i\\ g_i\left(t\right)& \mathrm{if}& |g_i\left(t\right)-m_i|\&le;{\mathrm{\&tau;}}_i\end{array}\right....$ （Ⅳ）

4) utilize filtered whole modal components w according to public formula V _i(t) carry out data reconstruction, do not contained the more purified EEG signals Y (t) of eye electricity;

$Y\left(t\right)=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{w}_i\left(t\right)...$ （Ⅴ）。

EEG signals X (t) decomposes all modal components c that obtain by EMD _i(t), because c _i(t) pass through bandpass filtering and threshold filter and handled the i that an obtains w _i(t) signal, filtering the eye electric artefact composition, so i w _i(t) signal can participate in data reconstruction, thereby has obtained purified relatively EEG signals Y (t).Solved comprising the empirical modal component c of the electric artefact of eye _i(t) manual screening problem, thus reach the purpose of from EEG signals, removing the electric artefact of eye automatically.

Description of drawings

The electric artefact of eye is removed structured flowchart in Fig. 1 EEG signals of the present invention;

The program flow diagram of the existing EMD decomposition method of Fig. 2;

The placement sketch map of the existing international standard ten-twenty electrode system of Fig. 3;

Fig. 4 adopts existing EMD method, to the result of EEG signals decomposition.

Fig. 5 adopts method of the present invention to experimentize, to the design sketch of the electric artefact removal of eye in the EEG signals.

The specific embodiment

Based on the electric artefact removal method of eye in the EEG signals of HHT, the complete procedure of its signal processing comprises the individual part of following four (1.2.3.4).Wherein, part 1 is existing method, and the feature of the present patent application comprises the individual part of three (2.3.4):

The signal processing complete procedure as shown in Figure 1.

1. the empirical modal of EEG signals decomposes (EMD).

EEG signals X (t) is a non-stationary signal, in order to study its transient state characteristic, must carry out empirical modal and decompose (EMD), obtains each IMF component c _i(t) and remainder r (t).

Wherein, i is for decomposing the number of the IMF component that obtains, and i is complete in data-driven, and the value of i is determined by X (t) fully, different X (t) signal, and the number i that decomposes the IMF component that obtains is different, this point has also embodied the adaptivity that EMD decomposes.This is different from wavelet transformation, and which floor carries out decompose, be to set by the people in advance.

The EMD step of decomposition:

The program circuit that EMD decomposes as shown in Figure 2.

The first step: determine all Local Extremum of time series X (t), with the cubic spline function curve all maximum points are carried out interpolation then, simulate the coenvelope line of primary signal X (t).In like manner, obtain the lower envelope line of X (t).All data points of signal are between two envelopes, and the meansigma methods of upper and lower envelope is designated as m (t).

Second step: the signal X (t) makes h with the difference note of m (t) ₁(t), ideally, h ₁(t) should be a basic model component.In general, it does not satisfy two required conditions of IMF (1. at whole data length, the number of extreme value and zero crossing must equate or differ from one at the most; 2. in any data point, on average being necessary for of the envelope of local maximum and the envelope of local minimum is zero, and promptly signal is about the local symmetry of time shaft), this seasonal X (t)=h ₁(t), carry out the first step, until h ₁(t) be a basic model component, note c ₁(t)=h ₁(t).

Two qualificationss of basic model component are (in any data point, on average being necessary for of the envelope of local maximum and the envelope of local minimum is zero, be that signal is about the local symmetry of time shaft) be a kind of theoretic requirement, in the screening process of reality, be difficult to guarantee that the local mean value of signal definitely is zero, therefore must provide the criterion that a screening process stops, whether it can utilize the value of two variance SD between the successive result less than preset threshold δ, is used as criterion:

$\mathrm{SD}=\underset{t=0}{\overset{T}{\mathrm{\&Sigma;}}}\frac{{(h_{1(k-1)}\left(t\right)-h_{1k}\left(t\right))}^2}{{h^2}_{1(k-1)}\left(t\right)}...\left(8\right)$

When SD＜δ, think two qualificationss that satisfied the basic model component, otherwise think dissatisfied.Experience shows: δ gets and was advisable in 0.2～0.3 o'clock, both can guarantee the linear and stable of IMF, can make IMF have corresponding physical significance again.

The 3rd step: decomposite a basic model component h from original series X (t) ₁(t) afterwards, deduct h with X (t) ₁(t), obtain the surplus value sequence X ₁(t)=X (t)-h ₁(t).

The 4th step: make X (t)=X ₁(t), be used as new " original " sequence, repeat above-mentioned one two three steps, extract the second, the three successively, until n basic model component c _n(t).At this moment, X _N+1(t) become a monotonic sequence, wherein no longer comprise the information of any pattern, the remainder r of primary signal that Here it is (t)=X _N+1(t).

2. all IMF component c that EEG signals EMD decomposed gained _i(t) utilize formula (3) to carry out Hilbert transform, utilize formula (7) to try to achieve instantaneous frequency f _i(t),, utilize formula (9) to carry out bandpass filtering, extract the signal g in the EEG signals frequency band based on instantaneous frequency _i(t).

$g_i\left(t\right)=\left\{\begin{array}{ccc}0& \mathrm{if}& f_i\left(t\right)>f_h\\ 0& \mathrm{if}& f_i\left(t\right)<f_l\\ g_i\left(t\right)& \mathrm{if}& \mathrm{other}\end{array}\right....\left(9\right)$

According to medical research, the EEG signals frequency band range is 1-30Hz, wherein can be divided into 6 frequency ranges, δ again ₂Wave band (1-2Hz), δ ₁Wave band (2-4Hz), θ wave band (4-8Hz), α wave band (8-13Hz), β ₁Wave band (3-20Hz), β ₂Wave band (20-30Hz).Here f _lBe set at 1Hz, f _hBe set at 30Hz.

The frequency domain filter that the bandpass filtering of this paper is different from the past, frequency definition in the past is of overall importance, is the frequency on the whole time shaft, has lost temporal information fully.The frequency of this paper is an instantaneous frequency, is the function of time, not only can reflect frequency information, and can reflect temporal information simultaneously.Therefore for the unstable signal filtering problem, adopting the wave filter based on instantaneous frequency, is only.

3. according to the statistical property and the big characteristics of the electric artifact amplitude of eye of IMF component, design threshold function table formula (10) and threshold filter function formula (12), utilized formula (10) (12) to carry out threshold filter, removed electro-ocular signal, obtained signal w _i(t).

τ _i＝mean|M _i-m _i|+std|M _i-m _i|........................(10)

Wherein, M _iBe i IMF component g _iThe value of each extreme point (t); Mean represents to average, and std represents to ask standard deviation.m _iBe i IMF component g _i(t) time average.Suppose time series g _i(t) length is N, then

$m_i=\frac{1}{N}\underset{t=1}{\overset{N}{\mathrm{\&Sigma;}}}{g}_i\left(t\right)...\left(11\right)$

In theory, IMF component g _i(t) time average m _iBe 0, but because when the EMD of reality decomposed, the time SD＜δ (δ gets 0.2～0.3) of employing was a judgment criterion, so can not guarantee IMF component g _i(t) time average m _iDefinitely be 0, therefore, formula (11) is of practical significance.

Threshold function table τ _i, characterized g _i(t) signal is with average m _iThe amplitude upper limit for datum line.On this basis, at the big characteristics of eye electricity artefact amplitude, designed threshold filter function formula (12).

$w_i\left(t\right)=\left\{\begin{array}{ccc}{m}_i& \mathrm{if}& |g_i\left(t\right)-m_i|>{\mathrm{\&tau;}}_i\\ g_i\left(t\right)& \mathrm{if}& |g_i\left(t\right)-m_i|\&le;{\mathrm{\&tau;}}_i\end{array}\right....\left(12\right)$

Threshold function table formula (10) (11) and threshold filter function formula (12) are to be determined by the statistical property of data, have realized the adaptive thresholding value filtering that complete data drives in some sense.

4. utilize filtered whole modal components w according to formula (13) _i(t) carry out data reconstruction, obtain

The more purified EEG signals Y (t) that does not contain the eye electricity.

$Y\left(t\right)=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{w}_i\left(t\right)...\left(13\right)$

The process (program circuit as shown in Figure 2) of being decomposed by EMD and two conditions of IMF are as can be known, remainder r (t) is not about the local symmetric signal of time shaft, can think the baseline drift of EEG signals, therefore remainder r (t) does not participate in reorganization, has solved the problem of the baseline drift of EEG signals well.

In order to provide concrete embodiment, and further specify effect of the present invention, the thinking eeg data that will provide in conjunction with Colorado state university EEG research center provides concrete implementation step and experimental result, verifies that this method goes the effect of artefact.

The data brief introduction:

Eeg data derives from Colorado state university EEG research center.The EEG signal collects as follows, and electrode leads the 10-20 system rest in C3, C4, P3, P4, these 6 positions of O1, O2 according to international standard, in addition, also synchronous acquisition one lead the EOG signal, reference electrode is placed on A1 and A2, as shown in Figure 3.The signals sampling frequency is 250Hz, and the analog filtering scope is 0.1～100Hz.Each EEG writing time is 10s, and data length N is 2500.It is EOG signals that O2 leads signal and synchronous acquisition that data are adopted in this experiment, as shown in Figure 5.

Implementation step:

1. according to EMD decomposition process (as Fig. 2) O2 is led EEG signals X (t) and carry out the EMD decomposition, obtain the IMF component c of each layer _i(t) and remainder r (t).Wherein δ gets 0.2, and decomposition result has obtained 8 layers of IMF component (c as shown in Figure 4 ₁(t), c ₂(t) ... c ₈And remainder r (t) (t)).

2. to 8 layers of IMF component c ₁(t), c ₂(t) ... c ₈(t) utilize formula (3) to carry out Hilbert transform, utilize formula (7) to try to achieve instantaneous frequency f ₁(t), f ₂(t) ... f ₈(t),, utilize formula (9) to carry out bandpass filtering, extract the signal g in the EEG signals frequency band based on instantaneous frequency ₁(t), g ₂(t) ... g ₈(t).

3. to signal g ₁(t), g ₂(t) ... g ₈(t) utilize formula (10) (12) to carry out threshold filter, remove electro-ocular signal, obtain signal w ₁(t), w ₂(t) ... w ₈(t).

Utilize filtered 8 modal components w according to formula (13) ₁(t), w ₂(t) ... w ₈(t) carry out data reconstruction, do not contained the more purified EEG signals Y (t) of eye electricity.Experimental result as shown in Figure 5.

Claims (1)

1. the automatic removal method of the electric artefact of eye in the self adaptation EEG signals is characterized in that, may further comprise the steps:

1) empirical modal of EEG signals decomposes

EEG signals X (t) is a non-stationary signal, carries out empirical modal and decomposes EMD, obtains each IMF component c _i(t) and remainder r (t);

Wherein, i is for decomposing the number of the IMF component that obtains, and the value of i is determined by X (t) fully;

2) EEG signals EMD is decomposed all IMF component c of gained _i(t) carry out Hilbert transform, try to achieve instantaneous frequency f _i(t),, utilize following public formula I to carry out bandpass filtering, extract the signal g in the EEG signals frequency band based on instantaneous frequency _i(t);

$g_i\left(t\right)=\left\{\begin{array}{ccc}0& \mathrm{if}& f_i\left(t\right)>f_h\\ 0& \mathrm{if}& f_i\left(t\right)<f_l\\ c_i\left(t\right)& \mathrm{if}& \mathrm{ohter}\end{array}\right.---\left(I\right);$ Here f _lBe 1Hz, f _hBe 30Hz;

3) utilize public formula II of formula threshold function table and threshold filter function formula (IV), carry out threshold filter, remove electro-ocular signal, obtain signal w _i(t);

τ _i＝mean|M _i-m _i|+std|M _i-m _i|.............................(Ⅱ)

Wherein, M _iBe i IMF component g _iThe value of each extreme point (t); Mean represents to average, and std represents to ask standard deviation; m _iBe i IMF component g _i(t) time average;

Suppose time series g _i(t) length is N, then

$m_i=\frac{1}{N}\underset{t=1}{\overset{N}{\mathrm{\&Sigma;}}}{g}_i\left(t\right)...\left(\mathrm{III}\right)$

Threshold function table τ _i, characterized g _i(t) signal is with average m _iThe amplitude upper limit for datum line; On this basis, at the big characteristics of eye electricity artefact amplitude, designed threshold filter function formula (IV);

$w_i\left(t\right)=\left\{\begin{array}{ccc}{m}_i& \mathrm{if}& |g_i\left(t\right)-m_i|>{\mathrm{\&tau;}}_i\\ g_i\left(t\right)& \mathrm{if}& |g_i\left(t\right)-m_i|\&le;{\mathrm{\&tau;}}_i\end{array}\right....\left(\mathrm{IV}\right)$

Utilize filtered whole modal components w according to formula (V) _i(t) carry out data reconstruction, do not contained the more purified EEG signals Y (t) of eye electricity;

$Y\left(t\right)=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{w}_i\left(t\right)...\left(V\right).$

Images