License

Nonlinear block for WinEEG (further – Software) is distributed as a plug-in for WinEEG program. This means:

1. All copyrights to Software are exclusively owned by the author – Alexey A. Mekler.

2. Software may be used only in conjunction with the only copy of WinEEG program.

3. No person or company may distribute Software without written permission from the copyright holder.

4. SOFTWARE IS DISTRIBUTED "AS IS". NO WARRANTY OF ANY KIND IS EXPRESSED OR IMPLIED. YOU USE AT YOUR OWN RISK. THE AUTHOR WILL NOT BE LIABLE FOR DATA LOSS, DAMAGES AND LOSS OF PROFITS OR ANY OTHER KIND OF LOSS WHILE USING OR MISUSING THIS SOFTWARE.

5. Using Software signifies acceptance of these terms and conditions of the license.


Additional copyrighting.

Nonlinear block for WinEEG includes some copyrighted programs as a part of distribution:

1. Gnuplot. Copyright © 1986 - 1993, 1998, 2004 Thomas Williams, Colin Kelley. WWW: http://www.gnuplot.info/

2. Cygwin (cygwin1.dll). Copyright © 2000, 2001,2002, 2003, 2004, 2005 Red Hat, Inc. WWW: http://www.cygwin.com/ .

Software works in conjunction with some programs from TISEAN project. These programs are subject of separate licensing and must be downloaded from TISEAN home page: http://www.mpipks-dresden.mpg.de/~tisean/TISEAN_2.1/index.html . TISEAN license terms could be found in the downloaded package or at TISEAN homepage as well. TISEAN at this time is freeware, but must be downloaded personally by the end user. We do not distribute TISEAN.


Citations

Algorithms, used in software are described in:

1. A. Mekler. Calculation of EEG correlation dimension: Large massifs of experimental data.// Computer Methods and Programs in Biomedicine, 92, (2008), 154–160

And also in:

2. Mekler A.A. Applications of the methods of nonlinear dynamical systems analysis for EEG signal processing.// In: Actual problems of modern mathematics: sci. notes., Vol. 13 (issue 2), ed. prof. Kalashnikov E.V., edition of LGU named by A.S. Pushkin, St.-Petersburg, 2004 (153 p.), pp. 112 – 140. (In Russian – Меклер А.А Применение аппарата нелинейного анализа динамических систем для обработки сигналов ЭЭГ. //В: сб. Актуальные проблемы современной математики: учёные записки. Т. 13 (вып. 2). п/ред. проф. Калашникова Е.В., изд. ЛГУ им. А.С. Пушкина, С.-Пб., 2004 г., 153 стр., стр. 112 – 140.)

3. Mekler A.A. About choosing of lag parameters for calculating of correlation dimension of reconstructed attractor of EEG. // In: Theoretical and methodical problems of secondary and high school education. Mathematics and informatics. Vol. 15, ed. prof. Hamov G.G., S.-Petersburg – Murmansk, pp. 42-45. (In Russian: Меклер А.А. К вопросу о выборе лаговых параметров при вычислении корреляционной размерности восстановленного аттрактора ЭЭГ.//In: Теоретические и методические проблемы обучения в школе и вузе. Математика и информатика. Межвузовский сборник научных трудов. Вып. 15, п/ред. проф. Хамова Г.Г., С.-Петербург – Мурманск, 2005, стр. 42-47.)

4. Mekler A.A. the use of the apparatus of nonliner analysis of dynamic systems for data processing of electroencephalogram. // Bulletin of New Medical Technologies. 2007, Vol. XIV, #1, pp. 73-77. (In Russian – А.А. Меклер. Применение аппарата нелинейного анализа динамических систем для обработки сигналов ЭЭГ. // Вестник новых медицинских технологий. 2007, вып. XIV, №1, стр. 73 – 77)

5. This documentation.

In [1] general view is given. In [3] it is shown, why do are default lag parameters (lag and Theiler window) chosen as they were set in software.

Calculation procedures, performed in Software, are based on algorithms, proposed in TISEAN package. These algorithms are described in: R. Hegger, H. Kantz, and T. Schreiber, Practical implementation of nonlinear time series methods: The TISEAN package, CHAOS 9, 413 (1999).


Hardware requirements:

The same as for WinEEG, but fast computer is preferable. For fast calculating processor frequency must be no less then 1 GHz.

Operation system requirements:

Program was developed using operation system Windows 2000 and tested with operation systems Windows 98, 2000 and XP.

Input signal requirements:

Duration:

There are contradictory requirements for duration of signals to be processed:

– On one hand they must be long enough, because most of theorems in this area are formulated and proved for endless signals. Routines performed in this software can’t work at all if signal is too short. There are many criteria of sufficiency duration of signal, but they have only recommendation character.

– On the other hand, analyzing of long time series takes too much time. Also, EEG, registered during long time, may differ from one fragment to another.

According to these reasons, there is lower limit of EEG duration:

– For plotting of attractor – 10 samples (1/25 second at sampling rate 250 Hz);

– For calculating [image] (any step) – 1000 samples (4 seconds at sampling rate 250 Hz).

If fragment, which is selected for EEG analysis, is shorter than required, respectable options are disabled.

Also, analyzing of signals, which duration is more than 50000 samples (200 seconds) may take a lot of time. Usually they take for analysis an EEG samples which duration is between 4000 and 10000 samples (16 ÷ 40 seconds at sampling rate 250 Hz).

Band range:

There are many publications, in which they use 30-Hz low-pass filter. However, there may be a reasons to use as wide band range as possible. Ideally – without 50 Hz notch filter. In the last case it is necessary to record EEG without any hum, that is very hard to do. So in most cases notch filter must be enabled.

Artifacts:

Before EEG processing all artifacts in it must be deleted or marked (with the corresponding tool in WinEEG).


General issues

Here will be described the main ideas of calculation of correlation dimension of reconstructed attractor.

I. General notes.

Let there is a dynamical system, described by several variables [image] and its behavior is described by some system of differential equations.

We can build a phase space for this system; it will be a space, defined by basis vectors – variables, which describe this system. So, at each moment state of this system will be determined by the set of values of respective variables and this set may be considered as coordinates of a point in the phase space. We say that this point represents state of this system. In next moment state of system changes and new point in phase space appears to represent state of system. We say that system moves in phase space. Set of points, through which system passes is called trajectory of system in phase space. It has coordinates [image] .

If the system has chaotic behavior, its trajectory is concentrated in limited region of phase space and never crosses itself [1]. This trajectory is called attractor of the system. This attractor has no integer dimensionality, or, as they say, fractal dimensionality (Latin – fracture).

There is a Takens theorem [2] and its main proposition is following:

We can take single realization of chaotic process [image] and we can set so called lag space. Lag space is [image] dimensional space with the following basis: , where [image] is time delay in time series [image] and [image] some constant. So, we can build a trajectory in lag space – [image] . This trajectory can be considered as reconstructed attractor of the system.

In the case of digitized signals we deal with discrete time series [image] Every point [image] of time series has its projection in the lag space [image] .

If [image] is large enough, then reconstructed attractor of the system will have the same metrical properties (dimension as well) as original one. As far as we are trying to embed reconstructed attractor in lag space, we call it embedding space and [image] embedding dimension [image] .

More detailed information can be obtained in [3, 4] and huge amount of other literature.

There are different methods of calculating fractal dimension of reconstructed attractor. Dimension, calculated by one of these methods is called “Correlation dimension” – [image] .


Reconstructed attractor plot.

Reconstructed attractor can be plotted in time-delay space. Of course it can be plotted only in 2- or 3-dimension space. The projection of multidimensional attractor into 2- or 3-dimentional space will not give its perfect reproduction. Yet it may be useful for conclusions about optimal choice of time delay [image] .


Main theoretical issues

1) Main theoretical issues.

Calculating of correlation dimension of reconstructed attractor is based on Grassberger-Procaccia algorithm [5]. This algorithm is used in nonlinear module for WinEEG with some modifications.

The central part of this algorithm is calculation of correlation integral:

,

where [image] – radius of some small sphere, [image] – number of points in reconstructed attractor

and [image] – Heaviside function: .

and [image] – [image] th and [image] th elements of time series and [image] – distance between projections of these elements in [image] dimensioned lag space.

Actually, correlation integral shows share of pairs of points of reconstructed attractor, that are located within distances less than [image] .

It is proved [5], that [image] , where [image] – dimension of reconstructed attractor.

So we can assume, that

.

Dimension of reconstructed attractor, calculated in this way is called correlation dimension of reconstructed attractor [image] . Dimension of lag space, in which attractor was reconstructed is called embedding dimension [image] .

So, we can set the following basic algorithm of calculations:

1) Setting [image] as 1 and calculating [image] ;

2) plotting [image] in log-log scale;

3) finding the linear fragment of this plot (so called scaling region);

4) Calculating of [image] as the slope of fragment that was found.

5) Increasing [image] and repeating all procedure. If [image] exceeds [image] , then after increasing of [image] [image] will increase too, other way – not.

On fig. 1 one can see log-log plots of [image] , calculated for [image] changing from 1 to 10[1].

Fig. 1. , plotted in log-log scale. changes from 1 (upper line) to 10 (lower line).

The classic way to find linear fragment of [image] , plotted in log-log scale, is to go so called Rapp graph. This plot shows the local slopes of [image] – fig. 2.

To estimate [image] at some value of [image] one has to find plateau on this plot and calculate the average value of [image] at this plateau.

Fig. 2. Rapp graphs – local slopes of curves, shown at fig. 1.

Thus one can calculate [image] for any [image] . To estimate proper [image] it is necessary to plot [image] . Value of [image] at which plot will go saturated can be considered enough for embedding and [image] at this value – real [image] (fig. 3).

Fig. 3. Saturating of curve at [image] shows that it is optimal value of embedding dimension and [image] at this value equals approx. 1.251.

[1] All sample calculations and graphs are given for Henon map – one of the classic chaotic processes with low-dimensional attractor.


Lag

Estimating of lag

If time series is endless, value of lag ( [image] ) has no matter [2]. If time series is finite, value of [image] depends on [image] . Exact algorithm of determining of [image] doesn’t exist, but there are some recommendations. According them delay coordinates must be independent. So, the simplest way to estimate lag is to set it as the time of autocorrelation (x-coordinate of the first zero or minimum of autocorrelation function). Instead of autocorrelation function one can use a function of mutual information. It takes into account not only linear, but also nonlinear correlations [6].


Theiler window

If time series is finite, then correlation integral must be calculated only using independent points. In this case, especially if original signal is represented with high sampling frequency, calculation of [image] may give spurious results, because points, placed in time series close to each other are temporally correlated [7]. To avoid it, it was proposed to exclude from calculation of correlation integral near-by-placed in time series points as temporarily correlated. I.e. in correlation integral must be no less than some constant. Minimal distance in time series between points to be included in calculation of correlation integral is called Theiler window.

There are different methods of estimating of Theiler window. It may be set as the first minimum or zero of autocorrelation function or time delayed mutual information function. It was proposed, that there is more correct method of estimating width of Theiler window with the help of so called space-time separation plot [8].

This plot is a curve of a constant probability density of two points of reconstructed attractor to be in [image] neighborhood, if distance between them in time series is [image] . Usually a set of curves are being plotted – for different values of probability density. This plot shows, how large the temporal distance between two points should be to consider them as independent [9]. See. fig. 4.

In this case we shall consider width of Theiler window as time of the first maximum, that is common for all curves (on fig. 4: ).

Fig. 4. Space-time separation plot. Probability densities are increasing from bottom to top.


Embedding dimension.

In the case of classic chaotic process value of embedding dimension can be set as value of saturation of the curve [image] . As far as EEG is noisy signal, this curve doesn’t go perfectly saturated. In this case one can evaluate proper [image] by the means of searching False Nearest Neighbors [12].

The main idea of this method is following.

Let attractor is reconstructed in m–dimensional space. Let’s find pairs of points, located close to each other. Then we can reconstruct this attractor in m+1– dimensional space. If dimension of attractor is less then m, then closest neighborhoods of points in m–dimensional space will project to the closest neighborhoods of the same points in m+1– dimensional space. But if dimension of attractor exceeds m, then distances between some points, that where close to each other in m–dimensional space, will rapidly increase in m+1– dimensional space. These points are called false nearest neighbors (FNN). When [image] increases relative amount of FNN decreases. Value of [image] when fraction of FNN is zero is sufficient for embedding (fig. 6).

Fig. 6. Fractions of false nearest neighbors for different signals.

This method sometimes gives satisfactory results, but not always. For example, if the signal is too noisy, fraction of FNN may not reduce to zero at all. Also, if time series is too short, fraction of FNN may be calculated not for all values of [image] . Thus, this method is quite useful, but not perfectly reliable. On fig.6 FNN fracture measured for EEG signal increases at [image] higher then 5. It means that original time series is too short (zero at [image] – spurious, meaningless value because effect of lack of time series length takes place already at [image] ).

Due to all these considerations there is a reason to make calculations for a some fixed value of [image] . As far as [image] of EEG signal is usually varies between 5 and 8, one can set [image] as, say, 10 by default and calculate [image] . This may have reason only if researcher goes to find changes in value of correlation dimension, calculated for EEG, recorded in different conditions. Ofcourse, it does not give us absolute value of correlation dimension. But often this is enough.


correlation dimension

Estimating of correlation dimension.

The simplest way to estimate [image] is to find plateau at Rapp graph (fig. 2).One can see that these graphs are not smooth enough to estimate dimension precisely. If we deal with real EEG signal – not classic chaotic process – then we can hardly find plateau and estimate [image] from these graphs. There is alternative estimator of [image] – Takens-Theiler estimator [10, 11]. It reads:

Graphs, plotted by such means are presented on fig. 5.

Fig. 5. Dimension estimation by means of Takens-Theiler estimator.


Fast calculation of correlation integral.

Calculating of distances between all points of reconstructed attractor may take a lot of time (this time is proportional to the second power of number of points in reconstructed attractor). It is possible to reduce amount of calculations by the following means.

We can take for each [image] only limited number of pairs of points (say, 1000) and make all calculations only for them. If we choose these points randomly, then the proportion of calculated distances (greater then [image] and less then [image] ) would be the same as if we’ll take all points. If [image] is too small for reliable statistics, whole point will be chosen.


Resume.

There is no theoretically proven methods of estimating parameters for calculating of [image] ( [image] , Theiler window, [image] ). One can use them in different combinations. To our opinion the best combination of these parameters is following:

– unit;

Theiler window – derived from space-time separating plot;

– fixed value (we set it as 10, but it may be something between 10 and 15).

Researchers have to mind that setting of embedding parameters also depends on concrete tasks to be solved.


Installation notes.

  1. Run setup file.
  2. Download from TISEAN homepage - http://www.mpipks-dresden.mpg.de/~tisean/TISEAN_2.1/index.html - TISEAN distribution, unpack it and extract from subdirectory .\exe_files\ following program files:
  • av-d2.exe
  • c2t.exe
  • cygwin.dll
  • d2.exe
  • delay.exe
  • false_nearest.exe
  • mutual.exe
  • rescale.exe
  • stp.exe
  3. Place these files into WinEEG directory.

After that you can run WinEEG and use nonlinear module, running it from menu “Analysis”.


This block performs nonlinear analysis of EEG

General notes.

Block of nonlinear analysis of EEG performs following types of EEG signal processing:

– Plotting of attractor, reconstructed in 2- or 3-D space;

– Calculating of correlation dimension of reconstructed attractor ( [image] ).

All procedures begin, of course, from selecting interval of input data. One can select for analysis some selection of EEG, some fragment, or whole EEG recording. If one use selection of EEG, it must be done before going nonlinear analysis (see main textbook on WinEEG).

After time interval of EEG is chosen, one can select a channel for analysis. This selection has no matter, if one is going to calculate [image] for all channels at once (see further description).

There are two text boxes, showing information that may be useful.

“Number of samples” shows, how many points of digitized EEG are selected.

“Sampling rate” – shows sampling rate that was set for digitizing of EEG.

EEG duration is related to number of samples in following way:

.

So, if it is shown that 1000 samples are selected and sampling rate is 250 Hz, it means, that four seconds of EEG will be processed.

It is important, because there are requirements for minimal duration of signal to be processed (see technical requirements).

NB! Calculation of [image] requires some additional parameters. It is strongly recommended to read literature on this subject. Yet, if user is in doubt, how to set options, [image] can be calculated using default settings. Our experience shows, that these settings are optimal for EEG processing.


This block performs nonlinear analysis of EEG

Step one. For both procedures it is necessary to set time delay [image] (lag). It may be unit or any other. It can be set manually or automatically – as the first minimum of time delayed mutual information function (see technical notes). If lag is set automatically, it is possible to preview it by pressing button “Estimate”. If it is not previewed, it would be still calculated automatically. Default lag – unit, i.e. (sampling rate)-1.

Time delay is shown in milliseconds, but one ought to mind, that its size can be only divisible to unit lag. If any other size of lag is set, it goes automatically enlarged to the closest size that is divisible to unit lag. For example, if sampling rate is 250 Hz, then unit lag is 1/250=4 milliseconds and possible lag values are 4, 8, 12, 16, and so on milliseconds. If it is set lag of 3 milliseconds, it will be automatically corrected to 4 milliseconds; 5 milliseconds will be corrected to 8 milliseconds and so on.

Time delay cannot be set as zero. In this case it is being automatically corrected to unit size.

After time delay is set, one can either plot reconstructed attractor in 2- or 3-D time delayed space or continue procedure of calculating of correlation dimension. Obviously, if researcher would like to get attractor plot and respective value of [image] , in both cases the same time delay must be used.


This block performs nonlinear analysis of EEG

Attractor plot.

Attractor, plotted in 2-D space, can be displayed at the plotting tool or stored in file. 3-D plot can be only displayed in plotting tool. After 3-D plot is displayed, it can be rotated by click-and-drugging.

Displayed attractor may be copied to clipboard by clicking left upper corner of display window and selecting respective menu item.

When the attractor is displayed in plotting tool, there are some digits are shown at the lower left corner. If it is 2-D plot, these digits mean x- and y-coordinates of cursor; if it is 3-D plot, they show value of angle of rotation around vertical and horizontal axes (degrees).


This block performs nonlinear analysis of EEG

Correlation dimension.

To perform calculating of [image] some more steps must be done.

Step two is setting of size of Theiler window. It can be set manually or automatically – derived from space-time separation plot (see technical notes). If Theiler window is set automatically, it is possible to preview it by pressing button “Estimate”. If it is not previewed, it would be still calculated automatically.

Like delay value (lag), size of Theiler window is expressed in milliseconds. So, it can be set according the same rules. Unlike delay value, size of Theiler window can be set as zero.

Step three is setting of dimension of lag space – embedding dimension ( [image] ). Embedding dimension can be set in diapason 8÷15 – optimal for EEG processing. Also, there is a tool for estimating the value of [image] – calculating of false nearest neighbors (FNN) in attractor, reconstructed in spaces of different dimensionality.

NB! Calculating of [image] may take a lot of time – depending on length of the selected for analysis time series. This time increases proportionally second power of the length of selected time series.

At this step one can set the regime of fast calculation of [image] – by checking box “Fast” (see technical notes). If this regime is set, program works much faster and results in most cases are precise enough. Yet, if researcher is going to find small fluctuations in [image] values (of different EEG’s for example), there is a reason to calculate [image] in common way.

Typical parameter of this regime – N = 1000. This parameter shows number of pairs of points, used for calculation of [image] . See technical notes for details.

Step four – calculating of [image] itself. There are two possible ways of output – value of [image] , calculated for established [image] , or plot [image] . In the last case largest value of [image] is established at the step three.

Estimating of [image] is modified. Here Takens-Theiler estimator is used (see technical notes).

Step five – calculating of [image] for all channels at once. At the output user will get a brain map of [image] and table with results. This procedure may take up to an hour or even more – depending on regime, time series length and computer productivity.

One can see that calculating of [image] needs to choose many optional parameters. A default set of parameters may be used by pressing button “Default”. The default parameters are following:

– lag ( [image] ) – unit ( [image] ; f – sampling rate of analogue-digital conversion of EEG signal; typically – 250 Hz);

– Theiler window – auto (see technical notes);

– Embedding dimension [image] - 10.

Calculating with these parameters will give satisfactory results in most cases.

In some cases calculations may take a lot of time and there may be need to interrupt them. It can be done by pressing button “Terminate”.


This block performs nonlinear analysis of EEG

Known bugs

  1. If window, which displays reconstructed attractor is, being moved (by drug and drop), cursor coordinates in it may disappear, or change their look. It is a bug of graphic tool “Gnuplot v.4.0” that is used here. To override it, just minimize this window and restore it back. Also one can slightly resize it.


Literature

1. Nicolis G., Prigogine I. Exploring complexity. An introduction. NY.,1989.

2. Takens F. Detecting strange attractors in turbulence. In: D.A. Rand and L.S. Young, editors, Dynamical Systems and Turbulence, Warwick 1980, Lecture Notes in Mathematics, 898, 366 – 381. Springer, Berlin, 1981.

3. Schuster H.G., Deterministic chaos. An introduction. Phisic-Verlag, Weinheim, 1984.

4. Schroeder M. Fractals, Chaos, Power Laws. Minutes from an Infinite Paradise. 1991.

5. Grassberger P., Procaccia I., Characterization of strange attractors. In: Physical Review Letters, 1983, v.50, pp. 346-349.

6. Fraiser A.M., Swinney H.L. Independent Coordinates for Strange Attractors From Mutual Information. Physical Rev. A, 33, #2, 1134 – 1140, 1986.

7. Theiler J. Spurious Dimension from Correlation Algorithms Applied to Limited Time Series Data. In: Physical Rev. A, 34, #3, 2427, 1986.

8. Provenzale A., Smith L. A., Vio R., and Murante G. Distinguishing Between Low-Dimensional Dynamics and Randomness in Measured Time Series. In: Physica D, 58, 31, 1992.

9. Hegger R., Kantz H., and Schreiber T., Practical Implementation of Nonlinear Time Series Methods. The TISEAN package. In: CHAOS 9, 413, 1999.

10. F. Takens, in: B. L. J. Braaksma, H. W. Broer, and F. Takens, eds., ``Dynamical Systems and Bifurcations'', Lecture Notes in Math. Vol. 1125, Springer, Heidelberg (1985).

11. J. Theiler, Lacunarity in a best estimator of fractal dimension, Phys. Lett. A 135, 195 (1988).

12. Kennel M. B., R. Brown, and H. D. I. Abarbanel. Determining embedding dimension for phase-space reconstruction using a geometrical construction, Phys. Rev. A, 45, 3403 (1992).

Additional literature:

Thomas Elbert; William J. Ray; Zbigniew J. Kowalik; James E. Skinner; Karl Eugen Graf; Niels Birbaumer. Chaos and physiology: deterministic chaos in excitable cell assemblies. Physiological Reviews, Jan 1994 v74 n1 p1(47).

Haken H. Principles of brain functioning. A synergetic approach to brain activity, behaviour and cognition. 1995. Springer.

Chaos in brain function. Ed.E.Basar-Berlin, Heidelberg, London. Springer-Verlag. 42–47. 1990.