Optimal Fractal Feature and Neural Network: EEG Based BCI Applications

1. Introduction

This chapter provides the theoretical principle of fractal dimension analysis applied for a brain-computer interface (BCI). Fractal geometry is a mathematical tool for dealing with complex systems. A method of estimating its dimension has been widely used to describe objects in space, since it has been found useful for the analysis of biological data. Moreover, fractal dimension (FD) is one of most popular fractal features. The term waveform applies to the shape of a wave, usually drawn as instantaneous values of a periodic quantity versus time. Any waveform is an infinite series of points. Aside of classical methods such as moment statistics and regression analysis, properties such as the Kolmogorov-Sinai entropy [1], the apparent entropy [2] and the FD [3] have been proposed to deal with the problem of pattern analysis of waveforms. The FD may convey information on spatial extent and self similarity and self affinity [4]. Unfortunately, although precise methods to determine the FD have already been proposed [5,6], their usefulness is severely limited since they are computer intensive and their evaluation is time consuming. Recently, the FD is relatively intensive to data scaling and shows a strong correlation with the human movement of EEG data [7-9]. The time series with fractal nature are to be describable by the functions of fractional Brownian motion (fBm), for which the FD can easily be set. Waveform FD values indicate the complexity of a pattern or the quantity of information embodied in a waveform pattern in terms of morphology, spectra, and variance. This chapter, we investigate most widely popular six algorithms to be feature patterns for a BCI in which algorithms are in time and frequency domain approaches. However, many methods in evaluating FDs have been developed. As the existing algorithms, they present different characteristics therefore the optimum condition respects to the requirements of BCI application should be seriously considered.

2. Experimental paradigm

By NLAB, under supervision of Prof.Masahiro Nakagawa, Department of Electrical Engieering, Nagaoka University of Technology, 1603-1 Kamitomioka, Nagaoka-shi 940-2188 Japan. E-mail: [e-mail] Http://pelican.nagaokaut.ac.jp/

In this study, we tested 7-male and 3-female. Their ages were 21-32 years (mean age, 25.3 years; SD ±3.6), with 180 trials per subject (30 trails per task). At the beginning of each EEG recording session, the subjects were instructed to relax with their eyes open for 30 s. The EEG data from this period was used as the baseline for the tasks. The duration of one trial was 6 s throughout the experiment. Each trail in the experiment was divided into two periods, a relaxing period of 0-3 s and an imaging period of 3-6 s. When the subjects had their eyes open on a blank screen, they would be in a relaxed state of mind (relaxing period). In the imaging period, the subjects imagined with eyes open, as represented by arrows. This study focuses on the motor movement functions in the brain. Therefore, imagination of motors movements are used as the tasks throughout this experiment. To assist the subject imagine the task without difficulty, we provide the acting stimuli in accordance to the given tasks. Explicit indicators are also shown on a monitor for informing subject during the experiment, all selected tasks are listed in Table 1.

Table 1.

Selected tasks in this study

3. Fractal analysis and methods

In this chapter, most widely popular six algorithms in time and frequency domain analysis have been addressed; the box-counting method (BCM), Higuchi‘s method (HM), variance fractal method (VFD), detrended fluctuation analysis (DFA), power spectral density analysis (PSDA), and critical exponent method (CEM). To measure the fractality in short time intervals of time-sequential data from one end of the waveform to the other sequentially, the dynamical changes in the FDs with respect the time series based on the time-dependent fractal dimensions (TDFD) were observed. For the classification process, two new feature parameters on the basis of fractal dimension; Kullback–Leibler (K-L) divergence and the different expected values were presented. In experimental results, DFA was selected to evaluate fraction dimensions on the basis of TDFDs [11, 36]. The reason is that, we found that the DFA provides fast computation time and also presents reasonable values in terms of accuracy and variability when comparison with each other.

3.1. Time series-based analysis

3.1.1. Box-Counting Method (BCM)

One of the most common methods for calculating the FD of a self-similar fractal is the box-counting method (BCM). The definition of BCM is a bounded set in the Euclidian $n$ -space that it composes self-similar property [25], $N_r\propto {\left(1/r\right)}^D\mathrm{ }$ by covering a structure with boxes of radii, $r$, the FD can be determined by

$$D_{BCM}=\underset{r\to 0}{\lim }\frac{{\log }_2\left(N_r\right)}{{\log }_2\left(1/r\right)}$$E1

We then repeat this process with several different radii. To implement this algorithm, number of contained box, $n_r$, is computed from the difference between the maximum and minimum amplitudes of the data divide by the changed radius, as

$$n_r\left(i\right)=⌊\frac{\text{max}\left(x_r\right)-\text{min}\left(x_r\right)}{r\left(i\right)}⌋,\mathrm{ } \text{for}\mathrm{ }\left\{r\in 2^k|k=\mathrm{1,2},\dots ,{\text{log}}_2\left(L\right)-1\right\}$$E2

$$N_r=\sum _i{n}_r\left(i\right)$$E3

where $N_r\mathrm{ }$ is a total number of contained boxes, $x_r$ represents the EEG time series with length $L$, $r\left(i\right)$ is a radius by changing a step of $k$ within the $i$-th subdivision window, and integer-part function denoted by $⌊\bullet ⌋$. To obtain the FD, least-square linear fitted line corresponds to the slope of the plot ${\text{log}}_2\left(N_r\right)\mathrm{ }$ versus ${\text{log}}_2\left(1/r\right)$ is applied.

3.1.2. Variance Fractal Dimension (VFD)

This method is determined by the Hurst exponent, $H$, whose calculation was divided from the properties of the fBm data. The basic idea of calculation is based on the power law relationship between the variance of the amplitude increments of the input time series, which was produced by a dynamic process over time. The main advantage of VFD was its support of the real-time computation [17]. We also selected the VFD for estimating the FD of EEG data. The amplitude increments of a datum over a time interval $∆t$ adhere to the following power law relationship $\text{Var[}x\text{(}{t}_2\text{-}{t}_1\text{)]}\text{∝}{\left|t_2\text{-}{t}_1\right|}^{2H}$, the Hurst exponent can be calculated by using a log-log plot then given by

$$H=\underset{∆t\to 0 }{\lim }\left(\frac{1}{2}\frac{{\log }_2\left(\mathrm{V}\mathrm{a}\mathrm{r}\left[{\left(∆x\right)}_{∆t}\right]\right)}{{\log }_2\left(∆t\right)} \right)$$E4

The variance in each window per stage k can be calculated as follows

$$\mathrm{V}\mathrm{a}\mathrm{r}\left[∆x_{∆t}\right]=\frac{1}{(N_k-1)}\left[\sum _{j=1}^{N_k}{\left(∆x\right)}_{jk}^2-\frac{1}{N_k}{\left(\sum _{j=1}^{N_k}{\left(∆x\right)}_{jk}\right)}^2\right]$$E5

$${\left(∆x\right)}_{jk}=x\left({jn}_k\right)-x\left(\left(j-1\right)n_k\right), \mathrm{f}\mathrm{o}\mathrm{r} j=\mathrm{1,2},\dots ,N_k$$E6

The least-square linear fitted line corresponds to the slope of the plot ${\text{log}}_2\left(n_k\right)$ and ${\text{log}}_2\left(\mathrm{V}\mathrm{a}\mathrm{r}{\left[∆x\right]}_k\right)$, the Hurst exponent is computed as $H=0.5s$ where s is the obtained slope then the FD can be estimated as

$$D_{VFD}=2-H$$E7

The process of calculating the FD essentially involves segmenting the entire input time series data into numerous subsequence (or window). The values $k$ represents the integer range chosen such that each window of size $N_T$ contains a number $n_k=2^k$ of smaller windows of size $N_k=⌊N_T/n_k⌋$.

3.1.3. Higuchi’s Method (HM)

FD is another measure of data complexity, generally evaluated in phase space by means of correlation dimension. Higuchi proposed an algorithm for the estimation of FD directly in time domain without reconstructing the strange attractor [18]. The HM also gives reasonable estimate of the FD in the case of short time segment. The HM algorithm based on the given finite time series $y=\left\{y\left(1\right), y\left(2\right), \dots ,y\left(N\right)\right\}$, a new time series, $y_m^k$, are constructed by the following equation

$$y_m^k=\left\{y\left(m\right),\mathrm{ }y\left(m+k\right),y\left(m+2k\right),\dots ,y\left(m+⌊\frac{N-m}{k}⌋\right)\bullet k\right\}\mathrm{ }\mathrm{f}\mathrm{o}\mathrm{r}\mathrm{ }m=\mathrm{1,2},\dots ,k$$E8

where both m and k are integers and they indicate the initial time and the time interval, respectively. The length, $L_m\left(k\right)$, of each curves is computed as

$$L_m\left(k\right)=\frac{1}{k}\left\{\left(\sum _{i=1}^{⌊\frac{N-m}{k}⌋}\left|y\left(m+ik\right)-y(m+(i-1\left)k\right)\right|\right)\bullet \frac{N-1}{⌊\frac{N-m}{k}⌋\bullet k}\right\}$$E9

The length of curve for time interval $k$, $L_m\left(k\right)$ is computed as the average of the $m$ curves. A relationship of this algorithm is $L_m\left(k\right)\propto k^{-D_{HM}}$ therefore we apply the least-squares fitting line of ${\text{log}}_2\left(L\left(k\right)\right)$ versus ${\text{log}}_2\left(k\right)$, the negative slope of the obtained line is calculated giving the estimate of the FD, $D_{HM}$.

3.1.4. Detrended Fluctuation Analysis (DFA)

The idea of DFA was invented originally to investigate the long-range dependence in coding and noncoding DNA nucleotide sequence [19]. Due to simplicity in implementation, the DFA is now becoming the most important method in the field of fractal analysis [20]. Therefore the DFA method was also applied to FD estimation in this study.

$$X\left(k\right)=\sum _{i=1}^k\left[x\left(i\right)-\stackrel{-}{x}\right]$$E10

This integrated series is divided into non-overlapping intervals of length n. In each interval, a least squares line is fit to the data (representing the trend in the interval). The series $X\left(k\right)$ is then locally detrended by subtracting the theoretical values $X_n\left(k\right)$ given by the regression. For a given interval length $n$, the characteristic size of fluctuation for this integrated and detrended series is calculated by

$$F\left(n\right)=\sqrt{\frac{1}{N}\sum _{i=1}^N{\left[X\left(k\right)-X_n\left(k\right)\right]}^2}$$E11

This computation is repeated over all possible interval lengths (in practice, the minimum length is around 10, and the maximum is a half-length of input data, giving two adjacent intervals). In this experiment, we use the power of 2 based length for input EEG data therefore in this case we can range $n=2^k$ for $k=\mathrm{4,5},\dots ,{\log }_2\left(L\right)-1$. A relationship is expected, as $F\left(n\right)\propto n^{\alpha }$ where $\alpha$is expressed as the slope of a double logarithmic plot ${\log }_2\left(F\left(n\right)\right)$ versus ${\log }_2\left(n\right)$. As PSD, Then $\alpha$ can be converted into the Hurst exponent $H=\alpha -1$ and the estimated FD according as$D_{DFA}=3-\alpha$.

3.2. Frequency-based analysis

3.2.1. Power Spectral Density Analysis (PSDA)

This method is widely used for assessing the fractal properties of time series, the method based on frequency domain analysis, and also works on the basis of the 1/f-like power scaling which can be obtained by the fast Fourier transform (FFT) algorithm. It is called a power-law relationship of Mandelbrot and van Ness [25] can be expressed as follows

$$f\left(s\right)\propto \frac{1}{s^{-\beta }}$$E12

where $s$ is the frequency, $f\left(s\right)$ is the correspondence power at a given frequency, and $\beta$ is the scaling exponent which is estimated by calculating the negative slope of the linear fit relating $\log \left(f\left(s\right)\right)$ to $\log \left(s\right).$ The PSDA was applied to raw EEG data to classify the EEG characteristics. This study used the enhanced version of PSDA that initially proposed by Fougère [26] and modified by Eke [27]. To improve the stability of spectral estimates, raw EEG data of each individual trial were first detrending by subtracted the mean of the data series from each value, and then each value is multiplied by a parabolic window is given as below

$$\psi \left(i\right)=1-{\left(\frac{2i}{M+1}-1\right)}^2,\mathrm{ } \mathrm{f}\mathrm{o}\mathrm{r}\mathrm{ }i=\mathrm{1,2},\dots ,M$$E13

Finally the scaling exponent was estimated by the least-square fitting of log-log domain for the frequencies lower than 1/8 of the maximum sampling rate, i.e., less than 64 Hz in this study. The Hurst exponent can be determined as $H=\left(\beta -1\right)/2$ where $1<\beta <3$ and the estimated FD can be computed by the following equation

$$D_{PSDA}=2+\frac{\left(1-\beta \right)}{2}$$E14

3.2.2. Critical Exponent Method (CEM)

This method was initially proposed by Nakagawa [21]. The CEM has been applied in the physiological data analysis and featuring [22-24]. The PSD of observed fBm data, $P_H\left(\upsilon \right)$, in the frequency domain is determined as

$$P_H\left(\upsilon \right)\propto {\upsilon }^{2H+1}={\upsilon }^{-\beta }$$E15

The moment, $I_{\alpha }$, and the moment exponent, $\alpha$, of the PSD is determined as

$$I_{\alpha }={\int }_1^{\mathrm{\Omega }}{P}_H\left(\upsilon \right){\upsilon }^{\alpha }d\upsilon ,\mathrm{ } (-\mathrm{\infty }<\alpha <+\infty )$$E16

We will consider the frequency bands as finite value and substitute Eq. (15) to Eq. (16) thus the equation was given as

$$I_{\alpha }\propto {\int }_1^{\mathrm{\Omega }}{\nu }^{\alpha -\beta }d\nu =\mathrm{ }{\int }_1^{\mathrm{\Omega }}{\nu }^{A-1}d\nu =\frac{1}{A}\left({\mathrm{\Omega }}^A-1\right)=\frac{2}{A}\text{exp}\left(\frac{uA}{2}\right)\text{sinh}\left(\frac{uA}{2}\right)$$E17

where $A=\alpha -\beta +1$, $u=\log \left(\mathrm{\Omega }\right)$, and $\mathrm{\Omega }$ is the upper frequency variable which was normalized to the lower bound of the integration region as 1. Hence taking the logarithm of moment and the third order derivative can be written as

$$\frac{d^3\text{log}{I}_{\alpha }}{{d\alpha }^3}=\frac{2}{8}{u}^3{\text{cosech}}^3\left(\frac{uA}{2}\right)\text{cosh}\left(\frac{uA}{2}\right)-\frac{1}{A^3}$$E18

We then determine the critical value $\alpha ={\alpha }_c$ which satisfies for Eq. (18) is equal to zero. Finally, the FD can be estimated by

$$D_{CEM}=2-\frac{{\alpha }_c}{2}$$E19

The PSDA and CEM used the FFT algorithm for transforming the EEG time series data into frequency components.

4. Performance assessment

Since high computation time and low variability are the requirement of the BCI system, the optimal algorithm for evaluating fractal dimension should be examined carefully. In this experiment, we selected two signals that fractal dimension value can be easily set for assessing. The usefulness of output results can be helped us know which algorithm is suited for applying it as the feature extractor.

4.1. Assessing with artificially generated signals

1. Fractional Brownian motion (fBm) signal

   fBm data were simulated by NLAB, Department of Electrical Engieering, Nagaoka University of Technology, 1603-1 Kamitomioka, Nagaoka-shi 940-2188 Japan.

   . General model of the fBm signal can be written as the following fractional differential equation [32]

$$\frac{{\mathrm{d}}^{H+1/2}{B}_H\left(t\right)}{{\mathrm{d}t}^{H+1/2}}\propto w\left(t\right)$$E20

The mathematical model for solving the fBm signal in Eq.(20) was proposed by Mandelbrot and Van Ness [25], which is defined by

$$B_H\left(t\right)=\frac{1}{\text{Γ}(H+1/2)}\left({\int }_{-\mathrm{\infty }}^0\left[{\left(t-s\right)}^{H-1/2}-{\left(-s\right)}^{H-1/2}\right]dB\left(s\right)+{\int }_0^t\left[{\left(t-s\right)}^{H-1/2}\right]dB\left(s\right)\right)$$E21

where $\text{Γ}\text{(∙)}$ is the Gamma function $\text{Γ}\left(a\right)={\int }_0^{+\infty }{x}^{a-1}{e}^{-x}dx$ and $H$ is the Hurst exponent where $0<H<1$. $\text{B}$ is a stochastic process. To simulate the fBm signal, we can simply implement the above definition to discrete time series at n-th index by Riemann’s type sums as follows:

$$B_H\left(n\right)=\frac{1}{\text{Γ}(h+1/2)}\left(\sum _{k=-b}^0\left[{\left(n-k\right)}^{H-\frac{1}{2}}-{\left(-k\right)}^{H-\frac{1}{2}}\right]G_1\left(k\right)+\sum _{k=0}^n\left[{\left(n-k\right)}^{H-\frac{1}{2}}\right]G_2\left(k\right)\right)$$E22

In case $H=1/2$, Eq.(4.22) can be reduced, we then obtain

$$B_{1/2\mathrm{ }}\left(n\right)=\sum _{k=0}^{n}G\left(k\right)$$E23

For $n=\mathrm{1,2},\dots ,N$ and $b=\sqrt{N^3}$. $G$ is a normal distribution space.

1. Synthetic (Syn) signal. This signal was produced using the deterministic Weierstrauss cosine function [28]. To simulate the Syn signal, we can simply implement the above definition to discrete time series at n-th index by [33]

$$W_h\left(n\right)=\mathrm{ }\sum _{i=0}^M{\gamma }^{-iH}\mathit{cos}\left(2\pi {\gamma }^{i}n\right)$$E24

where $\gamma >1$, and we set $M=26$ and $\gamma =5$. The fractal dimension of this signal is given by $D=2-H$.

4.2. Conditions for assessing

A summary of the selected parameter values in this experiment is shown in Table 2. The selection of signal lengths that are powers of two was motivated by the requirements of frequency-based algorithms (PSDA and CEM). Output performances of each algorithm were shown in Fig. 1.

4.3. Results of performance assessing

In order to test the effect of signal lengths on $H$ estimation and each algorithm was applied on the entire series. The estimated fractal dimension values are then averaged at the change of $H$. Before utilizing those of algorithms, we should regard as the following things:

1. Computation time. Short computation time is the requirement of the BCI applications. The fBm and Syn signals are performed to complete the entire series on a laptop 1.6 GHz Pentium with memory of 512 MB. Table 3 shows the results of the average computation time.

2. Accuracy. To compare with the theoretical FD value (true FD value), we used a mean-squared error (MSE) value, which can be defined by

$${\mathcal{E}}_{\mathrm{D}}=\frac{1}{L}{\sum }_{i=1}^L{\left(D_i-{\widehat{D}}_i\right)}^2$$E25

where $D$ is a theoretical FD value (true value), $\mathrm{ }\widehat{D}$ is an estimated FD value at ith step of $H$, and L is a signal length. Results are shown in Tables 4 and 5.

1. Variability (or robustness). The standard deviation (SD) value is computed to indicate that that variability of algorithms. Results are shown in Tables 6 and 7.

Algorithm	Parameter	Value
BCM	Step size of radius	$k=\mathrm{1,2},\dots ,{\log }_2\left(p\right)-1$
VFD	Step size of window	$k=\mathrm{2,3},\dots ,{\log }_2\left(p\right)-1$
HM	Maximum interval length	${ k}_{max}=p-3$
DFA	Step size of interval length	$k=\mathrm{4,5},\dots ,{\log }_2\left(p\right)-1$
PSDA	FFT point; Maximum frequency	$N_{FFT}=p$; $f_{max}=64\mathrm{H}\mathrm{z}$
CEM	FFT point; Step size of moment exponent	$N_{FFT}=p$; ${\alpha }_{∆}=0.001$

Table 2.

Selected parameters in the experiment.

Figure 1.

The results for fractal analysis, performed on a series of test EEG data. (a) BCM algorithm: a log-log plot of a total number of contained boxes versus radius. (b) VFD algorithm: a log-log plot of variance versus window length. (c) HM algorithm: a log-log plot of a curve’s length versus interval length. (d) DFA algorithm: a log-log plot of a total length $\mathrm{F}\left(\mathrm{n}\right)$ versus interval length. (e) PSDA algorithm: a log-log plot of power spectra versus frequency. (f) CEM algorithm: a log-log plot of moment function versus moment exponent value [16].

Signal length	${\mathit{t}}_{\mathit{a}\mathit{v}\mathit{r}}$ [ms]
BCM	VFD	HM	DFA	PSDA	CEM
L = 2^7	2.60	0.77	14.90	10.05	2.07	140.40
L = 2^8	3.50	1.10	148.60	20.20	2.90	227.50
L = 2^9	6.80	2.80	385.20	40.20	4.90	1,767.80
L = 2^10	12.60	4.30	697.90	82.50	10.80	3,351.20
L = 2^11	26.00	7.10	1,400.70	168.80	20.30	6,538.50
L = 2^12	54.00	12.40	2,490.50	329.70	46.40	12,964.80
L = 2^13	114.90	61.10	4,160.70	895.10	137.80	25,318.10
L = 2^14	265.50	164.70	14,630.80	2,971.50	478.70	50,157.60
L = 2^15	676.10	607.90	47,082.50	11,210.20	1,882.30	101,809.70
Average	129.11	95.79	7,890.20	1,747.53	287.35	22,475.06

Table 3.

Average computation time in millisecond (best value is marked in bold).

Signal length	${\mathcal{E}}_{\mathbf{D}}$(fBm signal)
BCM	VFD	HM	DFA	PSDA	CEM
L = 2^7	0.091	0.018	0.003	0.004	0.027	0.039
L = 2^8	0.069	0.013	0.031	0.002	0.031	0.047
L = 2^9	0.054	0.008	0.033	0.001	0.096	0.040
L = 2^10	0.044	0.007	0.001	0.001	0.099	0.045
L = 2^11	0.038	0.005	0.019	0.002	0.102	0.042
L = 2^12	0.032	0.003	0.015	0.001	0.102	0.041
L = 2^13	0.023	0.001	0.001	0.004	0.101	0.046
L = 2^14	0.022	0.004	0.001	0.001	0.105	0.046
L = 2^15	0.019	0.007	0.002	0.001	0.103	0.053
Average	0.044	0.007	0.012	0.002	0.085	0.044

Table 4.

MSE value performed with the fBm signal (best value is marked in bold).

Signal length	$\mathit{S}\mathit{D}$(fBm signal)
BCM	VFD	HM	DFA	PSDA	CEM
L = 2^7	0.083	0.009	0.054	0.003	0.059	0.038
L = 2^8	0.067	0.023	0.004	0.002	0.048	0.048
L = 2^9	0.044	0.002	0.003	0.009	0.173	0.026
L = 2^10	0.030	0.003	0.002	< 0.001	0.107	0.033
L = 2^11	0.020	0.057	0.003	< 0.001	0.098	0.041
L = 2^12	0.018	0.083	0.002	0.005	0.089	0.036
L = 2^13	0.026	0.105	0.004	0.022	0.027	0.044
L = 2^14	0.037	0.182	0.004	0.042	0.051	0.035
L = 2^15	0.049	0.183	0.003	0.065	0.029	0.029
Average	0.034	0.072	0.009	0.016	0.076	0.037

Table 5.

MSE value performed with the Syn signal (best value is marked in bold).

Signal length	$\mathit{S}\mathit{D}$(fBm signal)
BCM	VFD	HM	DFA	PSDA	CEM
L = 2^7	0.056	0.201	0.109	0.112	0.248	0.069
L = 2^8	0.047	0.150	0.089	0.095	0.154	0.056
L = 2^9	0.061	0.116	0.040	0.069	0.090	0.047
L = 2^10	0.034	0.088	0.032	0.049	0.056	0.042
L = 2^11	0.032	0.078	0.019	0.039	0.036	0.037
L = 2^12	0.023	0.079	0.016	0.036	0.028	0.031
L = 2^13	0.023	0.062	0.014	0.028	0.019	0.027
L = 2^14	0.018	0.054	0.013	0.024	0.014	0.026
L = 2^15	0.015	0.021	0.011	0.022	0.013	0.014
Average	0.034	0.094	0.038	0.053	0.073	0.039

Table 6.

Deviation performed with the fBm signal (best value is marked in bold).

Signal length	$\mathit{S}\mathit{D}$(Syn signal)
BCM	VFD	HM	DFA	PSDA	CEM
L = 2^7	0.039	0.151	0.056	0.085	0.278	0.036
L = 2^8	0.031	0.116	0.026	0.044	0.236	0.025
L = 2^9	0.023	0.074	0.022	0.043	0.141	0.046
L = 2^10	0.017	0.068	0.021	0.016	0.148	0.038
L = 2^11	0.018	0.070	0.019	0.011	0.165	0.030
L = 2^12	0.019	0.091	0.018	0.010	0.175	0.022
L = 2^13	0.022	0.063	0.017	0.009	0.127	0.020
L = 2^14	0.019	0.049	0.015	0.014	0.121	0.017
L = 2^15	0.008	0.011	0.009	0.012	0.119	0.012
Average	0.022	0.077	0.022	0.027	0.168	0.027

Table 7.

Deviation performed with the Syn signal (best value is marked in bold).

4.4. Time-Dependent Fractal Dimensions (TDFD)

One of common measure of irregularity in fluctuations as time-sequential data may be the fractal dimension. Fractal dimensions estimate the degree of freedom of fluctuations in a signal. To measure the fractality in short time intervals of time-sequential data from one end of the waveform to the other sequentially, we may observe the dynamical changes in the fractal dimensions with respect the time series. These fractal dimensions, namely, the time-dependent fractal dimensions (TDFD) [22].

In this study, EEG signal was sampled at 512 Hz (or $L=\mathrm{3,072}$ points) whose duration is 6 seconds. We set a window function is a rectangular type, window size = 512 points (1 s), moving window with intervals = 10 points (19.5 ms) because the temporal resolution of EEG is millisecond or even better [29]. The number of obtained points can compute from $N_{FD}=⌊(L-L_w)/∆t ⌋+1$ where $L$ is a signal length, $L_w$ is a window length ($L_w\le L$), and $∆t$ is an interval. Thus, we then obtained 257 points of fractal dimension. The time series-based algorithms of fractal dimension present more effective than frequency-based algorithms [16].

Figure 2.

Computation process of TDFDs.

4.5. Evaluation of features by Kullback-Leibler divergence

In probability theory and information theory, the Kullback–Leibler divergence [34] (or called K-L divergence) is a measure of the difference between two probability distributions: from an actual probability distribution $P$ to an arbitrary probability distribution $Q$.

$$D_{KL}\left(P\parallel Q\right)={\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}p\left(x\right)\log \frac{p\left(x\right)}{q\left(x\right)}dx$$E26

For probability distributions $P$ and $Q$ of a discrete random variable is defined to be

$$D_{KL}\left(P\parallel Q\right)=\sum _{i}P\left(i\right)\log \frac{P\left(i\right)}{Q\left(i\right)}$$E27

Important properties of the K-L divergence are

1. The K–L divergence is always positive $D_{KL}\left(P\parallel Q\right)\ge 0$,

2. $D_{KL}\left(P\parallel Q\right)=0\mathrm{ }$ if and only if $P=\mathrm{Q}$,

3. The K-L divergence is not symmetric $D_{KL}\left(P\parallel Q\right)\ne D_{KL}\left(Q\parallel P\right)$.

Figure 3.

Feature extraction concepts on the basis of the K-L divergence.

The different expected values between imaging and relaxing periods are also proposed to use as the featuring parameter together with the K-L divergence. In general, if $X$ is random variable defined on a probability space, then the different expected value of two random variables can be defined by

$${∆E\left\{X\right\}}_{PQ}=\sum _i{p}_i{x}_i-\sum _j{q}_i{x}_j$$E28

In this study, we applied the K-L divergence and different expected values for finding the features between relaxing and imaging periods. According to Eq.(28), $P\left(i\right)$ is regarded as relaxing period, and $Q\left(i\right)$ is regarded as imaging period. The probability distributions can be approximated by normalizing histogram of the data. The output patterns of K-L divergences versus the different expected value are also shown in Fig. 4. These obtained results can be classified by the neural network. More explanations of this classification process are in chapter 5.

Figure 4.

Left-hand movement imagination (LH-MI).

5. Neural network classifier

This study is concerned with the most common class of neural networks (NNs) called backpropagation algorithm. Backpropagation is an abbreviation for the backwards propagation of error. The learning process is processed based on the Levenberg–Marquardt (LM) backpropagation method [30] the Hessian matrix can be approximated as the following equation

$$\text{H}\approx 2{\mathbf{J}}^T\mathbf{J}$$E29

and the gradient can be determined as

$$\text{∇}\text{J}=2J^T\mathbf{e}$$E30

where $J$ is the Jacobian matrix that contains first derivatives of the network error with respect to weights and biases, and $\mathbf{ }\mathbf{e}$ is a vector of network errors.

$$J\left(w\right)=\left(\begin{array}{ccc}\frac{\partial e_1\left(w\right)}{\text{∂}{w}_1}& \cdots & \frac{\partial e_1\left(w\right)}{\text{∂}{w}_n}\\ ⋮& \ddots & ⋮\\ \frac{\partial e_N\left(w\right)}{\text{∂}{w}_1}& \cdots & \frac{\partial e_N\left(w\right)}{\text{∂}{w}_n}\end{array}\right)$$E31

The LM algorithm can use for updating the vector of weights, $\mathbf{w}$, by

$${\mathbf{w}}_{n+1}={\mathbf{w}}_n-{\left({\text{J}}^T\text{J}+\mu \text{I}\right)}^{-1}{\text{J}}^T\text{e}$$E32

$$∆\text{w}=-{\left({\text{J}}^T\text{J}+\mu \text{I}\right)}^{-1}{\text{J}}^T\text{e}$$E33

where $\text{I}$ is an identity matrix, $\mu$ is a learning rate, and $n$ is an iteration step. This algorithm is particularly well suited for training NN with MSE. Moreover, the LM algorithm has been shown to be much faster than backpropagation in a variety of applications [31].

Figure 5.

Classification diagram of the proposed method.

Task	Desired vector
LH-MI	${\left[\begin{array}{cc}1& 0\end{array}\mathrm{ } \begin{array}{cc}0& 0\end{array}\right]}^T$
RH-MI	${\left[\begin{array}{cc}0& 1\end{array}\mathrm{ } \begin{array}{cc}0& 0\end{array}\right]}^T$
FT-MI	${\left[\begin{array}{cc}0& 0\end{array} \begin{array}{cc}1& 0\end{array}\right]}^T$
TG-MI	${\left[\begin{array}{cc}0& 0\end{array} \begin{array}{cc}0& 1\end{array}\right]}^T$

Table 8.

Desired vectors from the proposed network

5.1. Cross-validation procedure

Then, the featured parameters typically perform in terms of fractal properties. The NN was selected as classifier and it includes adaptive self-learning that was learnt to produce the correct classifications based on a set of training examples via cross validation process. A tenfold cross validation method was used to increase the reliability of the results [13]. For each experiment, we divided all the data into three sets, namely, the training set (50% of the entire data), cross validation set (20% of the entire data), and testing set (30% of the entire data). According to the tasks, a set of desired values for network training are shown in Table 8. The MSE rate is succeeded versus the iteration step is shown in Fig. 6. In Table 9, we compared the proposed method with the conventional feature, namely, band power estimation (BPE) [14]. The training process makes the network to define its decision boundaries. Figures 7 to 8 show a tight decision boundary after training the network.

6. Discussions

This chapter discusses the major contributions of this study, 12 electrode channels were used to measure EEG signal over the sensorimotor area of the cortex. The accuracy rates of motor imagery tasks are satisfactory. The error in this classification is due to the ambiguous data between the relaxing and imaging periods, particularly when the subjects have imagined the movements. However, in this experiment, throughout the period of recoding raw EEG signals without training and also without rejecting a bad trial. On contrary, they can be different for different individuals. This makes the EEG signal depend on the measured subject. Even if different persons perform the same tasks and the measurement was done under identical conditions the resulting signals could differ [35].

Figure 6.

Log-log plot of learning curve.

Method	Accuracy [%]
LH-MI	RH-MI	FT-MI	TG-MI
BPE	77.6	79.5	76.4	77.3
Proposed	82.5	84.3	83.6	82.6

Table 9.

A comparison result of average accuracy rates.

Figure 7.

Decision boundary of the trained network.

Figure 8.

Corresponding 3-D surface of decision boundaries for each task.

Figure 9.

Updated weight values of the NN after learning; (Top panel) Proposed method, (Bottom panel) BPE Method.

The proposed method provides 2880 features for each subject, whereas BPE method provides 21600 features. The current features obtained from the proposed method present an easy practicable and reliable method for training the classification algorithms. Nevertheless, the number of features from BPE method can be reduced by changing the number of electrodes. We used the NN to classify the features processed by BPE method. After learning the NN, learning curves and decision boundaries of BPE method is shown in Fig. 9. As the results, output decision functions of classifying features based on BPE method does not perform well and cannot made clearly separation boundaries.

The six algorithms present specific properties in terms of capability and variability. Some algorithms appeared inapplicable for the EEG time series data in this study, such as PSDA and CEM algorithms, for example. The main reason was that two algorithms based on frequency analysis since estimating the Hurst exponent can be utilized by means of the FFT method where the FFT can be substantially applicable in the long-time series data, and it easily declines when the EEG data has been analyzed in a short-time series. However, we can enhance the variability and avoid such drawbacks by increasing the sampling rate at the same period of experiment. In terms of variability, the HM algorithm presents the lowest average values of SD. The DFA provides the fast computation time and also presents reasonable values in terms of accuracy and variability when comparison with each other. Although, VDF gives extremely fast computation time but the VFD itself has the drawbacks in terms of accuracy and robustness. The VFD algorithm appropriates for real-time application, since its consumed computation time was less than 10 ms. Moreover, as proposed algorithms we suggest that the most appropriate algorithm depends on an application purpose; for instance, the HM algorithm was suited for evaluating the FD of EEG data with high precision, whereas HM algorithm extremely consumed the processing time in comparison with the other algorithms. To show the performance of all algorithms, we then compared all estimated FDs with the theoretical value especially computed by short interval, i.e., 2⁷, 2⁸ and 2⁹ points; these results were shown in Figs. 10 to 15.

Figure 10.

Fractal dimension evaluated by BCM. (Right) fBm signal; (left) Syn signal.

Figure 11.

Fractal dimension evaluated by VFD. (Right) fBm signal; (left) Syn signal.

Figure 12.

Fractal dimension evaluated by HM. (Right) fBm signal; (left) Syn signal.

Figure 13.

Fractal dimension evaluated by DFA. (Right) fBm signal; (left) Syn signal.

Figure 14.

Fractal dimension evaluated by PSDA. (Right) fBm signal; (left) Syn signal.

Figure 15.

Fractal dimension evaluated by CEM. (Right) fBm signal; (left) Syn signal.

7. Conclusion

By selecting the best performance of fractal algorithm, we can obviously improve the rate of classification and can develop the novel methods in terms of fractal properties. Fractal features and the experimental framework can be applied not only binary-command BCIs, but also multi-command BCIs. The measurement of FD from EEG data can be capable not only in particular works, but also in publicity data. The FD characterizes the self-affine property of EEG data and has a direct relation to the different tasks.