# Evaluation of Methods for Estimating Fractal Dimension in Motor Imagery-Based Brain Computer Interface

Authors: Chu Kiong Loo; Andrews Samraj; Gin Chong Lee
Journal: Discrete Dynamics in Nature and Society
 (2011)

---

## Abstract

A brain computer interface BCI enables direct communication between a brain and a computer translating brain activity into computer commands using preprocessing, feature extraction, and classification operations. Feature extraction is crucial, as it has a substantial effect on the classification accuracy and speed. While fractal dimension has been successfully used in various domains to characterize data exhibiting fractal properties, its usage in motor imagery-based BCI has been more recent. In this study, commonly used fractal dimension estimation methods to characterize time series Katz's method, Higuchi's method, rescaled range method, and Renyi's entropy were evaluated for feature extraction in motor imagery-based BCI by conducting offline analyses of a two class motor imagery dataset. Different classifiers fuzzy k-nearest neighbours FKNN, support vector machine, and linear discriminant analysis were tested in combination with these methods to determine the methodology with the best performance. This methodology was then modified by implementing the time-dependent fractal dimension TDFD, differential fractal dimension, and differential signals methods to determine if the results could be further improved. Katz's method with FKNN resulted in the highest classification accuracy of 85%, and further improvements by 3% were achieved by implementing the TDFD method.


---

## Body

## 1. Introduction

A brain computer interface (BCI) enables direct communication between a brain and a computer translating brain activity into computer commands, thus providing nonmuscular interaction with the environment. Sensorimotor rhythms (SMRs) are rhythmic brain waves found in the frequency range of 8 to 12 Hz over the left and right sensorimotor cortices. Movement, movement preparation, and motor imagery desynchronize SMRs, whereas during relaxation or postmovement, they are synchronized [1]. Since motor imagery does not require any muscular activity, motor imagery-regulated SMRs are commonly utilized in BCI [2, 3]. This is particularly beneficial for people with neurological disorders, since their voluntary muscular activities might be impaired. Another advantage of the utilization of motor imagery-regulated SMRs in BCI is the short training period required [4]. Motor imagery tasks are identified by detecting the synchronization and desynchronization of SMRs. The most common motor imagery tasks are imagery hand [5], foot [4], and tongue [3] movements. Once acquired, SMRs are analyzed using preprocessing, feature extraction, and classification operations.Feature extraction is the process of accurately simplifying the representation of data by reducing its dimensionality while extracting its relevant characteristics for the desired task. It has a substantial effect on the classification accuracy and speed, since classification carried out without a successful feature extraction process on a high dimensional and redundant data would be computationally complex and would overfit the training data. Fractal dimension is a statistical measure indicating the complexity of an object or a quantity that is self-similar over some region of space or time interval. It has been successfully used in various domains to characterize such objects and quantities [6, 7], but its usage in motor imagery-based BCI has been more recent [8, 9]. There are several fractal dimension estimation methods, some of which are not applicable to all types of data exhibiting fractal properties. In order to achieve a higher classification accuracy and speed, the fractal dimension estimation method that is most suitable to the data at hand should be chosen.In this study, Katz’s method [10], Higuchi’s method [11], the rescaled range (R/S) method [12], and Renyi’s entropy [13] were evaluated for feature extraction in motor imagery-based BCI by conducting offline analyses of a two-class motor imagery dataset. Fuzzy k nearest neighbors (FKNN), support vector machine (SVM), and linear discriminant analysis (LDA) were tested in combination with these methods to determine the methodology with the best performance. This methodology was then modified by implementing time-dependent fractal dimension (TDFD) [14], differential fractal dimension, (DFD) and differential signals (DS) [15].


## 2. Materials and Methods

### 2.1. Dataset

The motor imagery dataset from the BCI Competition II (Data set III) provided by the Department of Medical Informatics, Institute for Biomedical Engineering, University of Technology Graz was analyzed. The data was acquired over seven runs from a healthy 25-year-old female subject during imagery left and right hand movements. The signals were recorded with a sampling rate of 128 Hz from three electrodes placed at the standard positions of the 10–20 international system (C3, Cz, and C4) and filtered between 0.5 and 30 Hz. Each run consisted of 40 trials and each trial was nine seconds long. During the first two seconds of each trial, neither a stimulus was presented nor did the subject perform any motor imagery task. After this period, an acoustic and a visual stimulus indicating the beginning of the motor imagery task were presented. Then, for six seconds, a cue (a left or right arrow) indicating the required motor imagery task was presented (in a random order for each trial), and the subject performed this task. During this period, a feedback bar was displayed. Both the training and testing sets consisted of 140 samples.


### 2.2. Preprocessing

The samples from each electrode were zero phase filtered using a 6th-order bandpass digital Butterworth filter with cutoff frequencies of 0.5 and 30 Hz in both the forward and reverse directions. The last six seconds of each trial were extracted to discard the period without any motor imagery. Two different electrode configurations (C3 and C4, and C3, Cz, and C4) were tested.


### 2.3. Feature Extraction

In Katz’s method, Higuchi’s method, and theR/S method, the fractal dimension of the samples from selected electrodes were concatenated into feature vectors. In the TDFD, DFD, and DS methods, the fractal dimensions were estimated using the fractal dimension estimation method of the methodology with the best performance.


#### 2.3.1. Katz’s Method

Katz’s method [10] calculates the fractal dimension of a sample as follows: the sum and average of the Euclidean distances between the successive points of the sample (L and a, resp.) are calculated as well as the maximum distance between the first point and any other point of the sample (d). The fractal dimension of the sample (D) then becomes(2.1)D=log⁡(L/a)log⁡(d/a)=log⁡(n)log⁡(n)+log⁡(d/L),
where n is L divided by a.


#### 2.3.2. Higuchi’s Method

Higuchi’s method [11] calculates the fractal dimension of a sample as follows: first, subsample sets (Xk) are constructed from the sample (X) as(2.2)Xkm={X(m+ik)}i=0[(N-m)/k],
where k∈ [1, kmax], m∈ [1, k] and N is the sample size. Then, the length of each Xk(Lm) is calculated as (2.3)Lm(k)=(∑i=1[(N-m)/k]|X(m+ik)-X(m+(i-1)k)|(N-1)/([(N-m)/k]k))k.
Finally, the fractal dimension of the sample (D) is solved from(2.4)〈L(k)〉∝k-D,
where 〈L〉 is the average of Lm. Three kmax values from the range of 8 to 18 [16] (8, 13, and 18) were tested.


#### 2.3.3. R/S Method

TheR/S method [12] calculates the fractal dimension of a sample by iteratively dividing it into nonoverlapping subsamples with decreasing subsample size and performing the following operations at each iteration: for each subsample, a new subsample (X) is constructed from its zero mean (ξ) such that the nth point of X is the cumulative sum of the first n points of ξ. Then, the difference between the maximum and the minimum values, and the standard deviation of X (R and S, resp.) are calculated in order to obtain their ratio (R/S). Finally, R/S of each X is averaged ((R/S)avg). After obtaining (R/S)avg at each iteration, the Hurst exponent (H) becomes the slope of the log-log plot of (R/S)avg versus subsample size. The fractal dimension then becomes 2-H.


#### 2.3.4. Renyi’s Entropy

Renyi’s entropy [13] is generalization of Shannon’s entropy. Renyi’s entropy is defined as(2.5)Rq=11-qlog2(∑i=1npiq),
where q>0,   q≠1. However, R1 exists and the value is subset of Shannon’s entropy. Therefore, Shannon’s entropy is limit case to Rq when q=1(2.6)lim⁡q→1Rq=S.
This is the fractal dimension on the basic of Renyi’s entropy(2.7)Dq=lim⁡δx→011-qlog2log2∑i=1npiqlog2δx,
can be replaced with faster algorithm(2.8)Dq=lim⁡δx→0log2∑i=1nNiq-qlog2Nlog2δx.


#### 2.3.5. TDFD Method

In TDFD method, a window (with sizes) is slid over a sample by a time step, and the fractal dimension of the part of the sample inside the window is estimated. The fractal dimensions were concatenated into feature vectors. Different window sizes were tested using a time step of one second.


#### 2.3.6. DFD and DS Methods

The DFD method is a variation of the DS method. In the DFD method, first, the fractal dimensions of the samples from selected electrodes are estimated, and then, the pairwise differences of the fractal dimensions are calculated. However, in the DS method [15], first, the pairwise differences of the samples from selected electrodes are calculated, and then, the fractal dimensions of the pairwise differences are estimated. In both methods, the resultant values were concatenated into feature vectors. Only the three electrode configuration was tested, since the two electrode configuration results in one-dimensional feature vectors.


### 2.4. Classification

After constructing the feature vectors, the test samples were classified as imagery left or right hand movements using different classifiers. FKNN, SVM, and LDA were tested.FKNN is a variation of KNN. The main difference between the two is that KNN assigns a class label to a sample that is most frequent among the k nearest neighbors of that sample, whereas FKNN assigns a membership value for each class in this neighborhood and classifies the sample as the class with the highest membership value. The membership value for a class was calculated by dividing the sum of the distances between the samples belonging to this class and the test sample by the sum of the distances between all the samples in the neighborhood and the testing sample. Number of nearest neighbors between one and the square root of the sample length were tested.SVM separates the samples using a hyperplane that maximizes the margin between those belonging to different classes. SVM with a linear kernel was used.LDA finds a linear combination of features that best separates the samples belonging to different classes and can be used as a classifier. To assign a class label to a sample, the probabilities of the sample belonging to each class were estimated using LDA. The label of the class with the highest probability was then assigned to the sample.


## 2.1. Dataset

The motor imagery dataset from the BCI Competition II (Data set III) provided by the Department of Medical Informatics, Institute for Biomedical Engineering, University of Technology Graz was analyzed. The data was acquired over seven runs from a healthy 25-year-old female subject during imagery left and right hand movements. The signals were recorded with a sampling rate of 128 Hz from three electrodes placed at the standard positions of the 10–20 international system (C3, Cz, and C4) and filtered between 0.5 and 30 Hz. Each run consisted of 40 trials and each trial was nine seconds long. During the first two seconds of each trial, neither a stimulus was presented nor did the subject perform any motor imagery task. After this period, an acoustic and a visual stimulus indicating the beginning of the motor imagery task were presented. Then, for six seconds, a cue (a left or right arrow) indicating the required motor imagery task was presented (in a random order for each trial), and the subject performed this task. During this period, a feedback bar was displayed. Both the training and testing sets consisted of 140 samples.


## 2.2. Preprocessing

The samples from each electrode were zero phase filtered using a 6th-order bandpass digital Butterworth filter with cutoff frequencies of 0.5 and 30 Hz in both the forward and reverse directions. The last six seconds of each trial were extracted to discard the period without any motor imagery. Two different electrode configurations (C3 and C4, and C3, Cz, and C4) were tested.


## 2.3. Feature Extraction

In Katz’s method, Higuchi’s method, and theR/S method, the fractal dimension of the samples from selected electrodes were concatenated into feature vectors. In the TDFD, DFD, and DS methods, the fractal dimensions were estimated using the fractal dimension estimation method of the methodology with the best performance.


### 2.3.1. Katz’s Method

Katz’s method [10] calculates the fractal dimension of a sample as follows: the sum and average of the Euclidean distances between the successive points of the sample (L and a, resp.) are calculated as well as the maximum distance between the first point and any other point of the sample (d). The fractal dimension of the sample (D) then becomes(2.1)D=log⁡(L/a)log⁡(d/a)=log⁡(n)log⁡(n)+log⁡(d/L),
where n is L divided by a.


### 2.3.2. Higuchi’s Method

Higuchi’s method [11] calculates the fractal dimension of a sample as follows: first, subsample sets (Xk) are constructed from the sample (X) as(2.2)Xkm={X(m+ik)}i=0[(N-m)/k],
where k∈ [1, kmax], m∈ [1, k] and N is the sample size. Then, the length of each Xk(Lm) is calculated as (2.3)Lm(k)=(∑i=1[(N-m)/k]|X(m+ik)-X(m+(i-1)k)|(N-1)/([(N-m)/k]k))k.
Finally, the fractal dimension of the sample (D) is solved from(2.4)〈L(k)〉∝k-D,
where 〈L〉 is the average of Lm. Three kmax values from the range of 8 to 18 [16] (8, 13, and 18) were tested.


### 2.3.3. R/S Method

TheR/S method [12] calculates the fractal dimension of a sample by iteratively dividing it into nonoverlapping subsamples with decreasing subsample size and performing the following operations at each iteration: for each subsample, a new subsample (X) is constructed from its zero mean (ξ) such that the nth point of X is the cumulative sum of the first n points of ξ. Then, the difference between the maximum and the minimum values, and the standard deviation of X (R and S, resp.) are calculated in order to obtain their ratio (R/S). Finally, R/S of each X is averaged ((R/S)avg). After obtaining (R/S)avg at each iteration, the Hurst exponent (H) becomes the slope of the log-log plot of (R/S)avg versus subsample size. The fractal dimension then becomes 2-H.


### 2.3.4. Renyi’s Entropy

Renyi’s entropy [13] is generalization of Shannon’s entropy. Renyi’s entropy is defined as(2.5)Rq=11-qlog2(∑i=1npiq),
where q>0,   q≠1. However, R1 exists and the value is subset of Shannon’s entropy. Therefore, Shannon’s entropy is limit case to Rq when q=1(2.6)lim⁡q→1Rq=S.
This is the fractal dimension on the basic of Renyi’s entropy(2.7)Dq=lim⁡δx→011-qlog2log2∑i=1npiqlog2δx,
can be replaced with faster algorithm(2.8)Dq=lim⁡δx→0log2∑i=1nNiq-qlog2Nlog2δx.


### 2.3.5. TDFD Method

In TDFD method, a window (with sizes) is slid over a sample by a time step, and the fractal dimension of the part of the sample inside the window is estimated. The fractal dimensions were concatenated into feature vectors. Different window sizes were tested using a time step of one second.


### 2.3.6. DFD and DS Methods

The DFD method is a variation of the DS method. In the DFD method, first, the fractal dimensions of the samples from selected electrodes are estimated, and then, the pairwise differences of the fractal dimensions are calculated. However, in the DS method [15], first, the pairwise differences of the samples from selected electrodes are calculated, and then, the fractal dimensions of the pairwise differences are estimated. In both methods, the resultant values were concatenated into feature vectors. Only the three electrode configuration was tested, since the two electrode configuration results in one-dimensional feature vectors.


## 2.3.1. Katz’s Method

Katz’s method [10] calculates the fractal dimension of a sample as follows: the sum and average of the Euclidean distances between the successive points of the sample (L and a, resp.) are calculated as well as the maximum distance between the first point and any other point of the sample (d). The fractal dimension of the sample (D) then becomes(2.1)D=log⁡(L/a)log⁡(d/a)=log⁡(n)log⁡(n)+log⁡(d/L),
where n is L divided by a.


## 2.3.2. Higuchi’s Method

Higuchi’s method [11] calculates the fractal dimension of a sample as follows: first, subsample sets (Xk) are constructed from the sample (X) as(2.2)Xkm={X(m+ik)}i=0[(N-m)/k],
where k∈ [1, kmax], m∈ [1, k] and N is the sample size. Then, the length of each Xk(Lm) is calculated as (2.3)Lm(k)=(∑i=1[(N-m)/k]|X(m+ik)-X(m+(i-1)k)|(N-1)/([(N-m)/k]k))k.
Finally, the fractal dimension of the sample (D) is solved from(2.4)〈L(k)〉∝k-D,
where 〈L〉 is the average of Lm. Three kmax values from the range of 8 to 18 [16] (8, 13, and 18) were tested.


## 2.3.3. R/S Method

TheR/S method [12] calculates the fractal dimension of a sample by iteratively dividing it into nonoverlapping subsamples with decreasing subsample size and performing the following operations at each iteration: for each subsample, a new subsample (X) is constructed from its zero mean (ξ) such that the nth point of X is the cumulative sum of the first n points of ξ. Then, the difference between the maximum and the minimum values, and the standard deviation of X (R and S, resp.) are calculated in order to obtain their ratio (R/S). Finally, R/S of each X is averaged ((R/S)avg). After obtaining (R/S)avg at each iteration, the Hurst exponent (H) becomes the slope of the log-log plot of (R/S)avg versus subsample size. The fractal dimension then becomes 2-H.


## 2.3.4. Renyi’s Entropy

Renyi’s entropy [13] is generalization of Shannon’s entropy. Renyi’s entropy is defined as(2.5)Rq=11-qlog2(∑i=1npiq),
where q>0,   q≠1. However, R1 exists and the value is subset of Shannon’s entropy. Therefore, Shannon’s entropy is limit case to Rq when q=1(2.6)lim⁡q→1Rq=S.
This is the fractal dimension on the basic of Renyi’s entropy(2.7)Dq=lim⁡δx→011-qlog2log2∑i=1npiqlog2δx,
can be replaced with faster algorithm(2.8)Dq=lim⁡δx→0log2∑i=1nNiq-qlog2Nlog2δx.


## 2.3.5. TDFD Method

In TDFD method, a window (with sizes) is slid over a sample by a time step, and the fractal dimension of the part of the sample inside the window is estimated. The fractal dimensions were concatenated into feature vectors. Different window sizes were tested using a time step of one second.


## 2.3.6. DFD and DS Methods

The DFD method is a variation of the DS method. In the DFD method, first, the fractal dimensions of the samples from selected electrodes are estimated, and then, the pairwise differences of the fractal dimensions are calculated. However, in the DS method [15], first, the pairwise differences of the samples from selected electrodes are calculated, and then, the fractal dimensions of the pairwise differences are estimated. In both methods, the resultant values were concatenated into feature vectors. Only the three electrode configuration was tested, since the two electrode configuration results in one-dimensional feature vectors.


## 2.4. Classification

After constructing the feature vectors, the test samples were classified as imagery left or right hand movements using different classifiers. FKNN, SVM, and LDA were tested.FKNN is a variation of KNN. The main difference between the two is that KNN assigns a class label to a sample that is most frequent among the k nearest neighbors of that sample, whereas FKNN assigns a membership value for each class in this neighborhood and classifies the sample as the class with the highest membership value. The membership value for a class was calculated by dividing the sum of the distances between the samples belonging to this class and the test sample by the sum of the distances between all the samples in the neighborhood and the testing sample. Number of nearest neighbors between one and the square root of the sample length were tested.SVM separates the samples using a hyperplane that maximizes the margin between those belonging to different classes. SVM with a linear kernel was used.LDA finds a linear combination of features that best separates the samples belonging to different classes and can be used as a classifier. To assign a class label to a sample, the probabilities of the sample belonging to each class were estimated using LDA. The label of the class with the highest probability was then assigned to the sample.


## 3. Results

The classification accuracies (Table1) and the computation times (Table 2) were evaluated for each fractal dimension calculation method and classifier combination. Katz’s method was the fastest method, and combining it with FKNN, the highest classification accuracy of 85% (the three electrode configuration and k=9) as well as the second highest classification accuracy of 83% (the two electrode configuration and k=9) were achieved. R/S method with any classifier performed the slowest with the classification accuracies and the computation times ranging from 69% to 71% and 7.32 to 11.07 s, respectively. On the other hand, Renyi’s entropy with any classifier performed the worst with the classification accuracies and the computation times ranging from 55% to 66% and 1.84 to 4.87 s, respectively. The performances of the rest of the combinations were similar (Tables 1 and 2). The classification accuracies (except for the R/S method and Renyi’s entropy) and computation times (except for the R/S method) increased with the number of the number of selected electrodes.Table 1
Maximum classification accuracies (percentage of correctly classified testing samples) obtained by the combination of fractal dimension estimation methods and classifiers with the two and three electrode configurations (and the parameters (k and kmax) used to obtain these values).
Classification accuracy (%)Katz’s methodHiguchi’s methodR/S methodRenyi’s entropyC3, C4C3, Cz, C4C3, C4C3, Cz, C4C3, C4C3, Cz, C4C3, C4C3, Cz, C4FKNN83 (k=9)85 (k=9)77(k=7, kmax=18)79 (k=5, kmax=18)71 (k=9)69 (k=9)66 (k=9, q=4)65 (k=9, q=5)SVM777978 (kmax=13)81 (kmax=13)717059 (q=2)55 (q=4)LDA788178 (kmax=13)79 (kmax=13)717059 (q=3)57 (q=18)Table 2
Computation times (time it took for feature extraction and classification) corresponding to the maximum classification accuracies obtained by the combination of fractal dimension estimation methods and classifiers with the two and three electrode configurations (and the parameters (k and kmax) used to obtain these values).
Computation time (s)Katz’s methodHiguchi’s methodR/S methodRenyi’s entropyC3, C4C3, Cz, C4C3, C4C3, Cz, C4C3, C4C3, Cz, C4C3, C4C3, Cz, C4FKNN0.17 (k=9)0.23 (k=9)1.03 (k=7, kmax=18)1.5 (k=5, kmax=18)7.37 (k=9)11.07 (k=9)2.35 (k=9, q=4)2.60 (k=9, q=5)SVM0.340.341.07 (kmax=13)1.4 (kmax=13)7.3610.991.90 (q=2)4.87 (q=4)LDA0.120.210.83 (kmax=13)1.26 (kmax=13)7.3210.991.84 (q=3)3.01 (q=18)Table3 shows the computation times and classification accuracies obtained by modifying the best performing methodology. Although all the modifications increased the computation time, further improvements in the classification accuracy (by 3%) were achieved only by implementing TDFD method (the two channel configuration, k=5 and s=10). However, implementing the DFD and DS methods resulted in lower classification accuracies.Table 3
Computation times and classification accuracies obtained by modifying the highest performing methodology (Katz’s Method with FKNN) (and the parameters (k and s) used to obtain these values).
TDFD methodDFD methodDS methodC3, C4C3, Cz, C4C3, Cz, C4C3, Cz, C4Classification accuracy (%)88 (k=5, s=10)85 (k=7, s=64)84 (k=11)71 (k=11)Computation time (s)3.47 (k=5, s=10)0.94 (k=7, s=64)0.41 (k=11)0.26 (k=11)Mental activity may modulate FD of EEG signal which implies that it is timed-ependent in nature. By implementing TDFD method in Katz’s Method with FKNN, we may measure the fractality in short time intervals of time-sequential data from one end of the waveform to the other sequentially, and we may observe the dynamical changes in the FDs with respect the time series. These FDs, namely, are referred to the time-dependent fractal dimensions (TDFD) [17, 18].Katz’s algorithm is the most consistent method due to its exponential transformation of FD values and relative insensitivity to noise. Hiaguchi’s method, however, yields a more accurate estimation of signal FD, when tested on synthetic data, but it is more sensitive to noise. In the experiment, EEG datasets used are real data sets which contain noise, hence Katz’s method exhibits better result [19].


## 4. Conclusions

Since all fractal dimension estimation methods are not applicable to all types of data exhibiting fractal properties, commonly used fractal dimension estimation methods to characterize time series with different classifiers were evaluated to find the most suitable method for motor imagery data. Katz’s method with FKNN was determined to be the best methodology, and the results were further improved by implementing the TDFD method. The results warrant further research to use this methodology in online analysis of motor imagery data and analysis of other signals.



---

*Source: 724697-2011-11-30.xml*