# Common Spatio-Time-Frequency Patterns for Motor Imagery-Based Brain Machine Interfaces

Authors: Hiroshi Higashi; Toshihisa Tanaka
Journal: Computational Intelligence and Neuroscience
 (2013)

---

## Abstract

For efficient decoding of brain activities in analyzing brain function with an application to brain machine interfacing (BMI),
we address a problem of how to determine spatial weights (spatial patterns), bandpass filters (frequency patterns), and time windows
(time patterns) by utilizing electroencephalogram (EEG) recordings. To find these parameters, we develop a data-driven
criterion that is a natural extension of the so-called common spatial patterns (CSP) that are known to be effective features in BMI.
We show that the proposed criterion can be optimized by an alternating procedure to achieve fast convergence. Experiments
demonstrate that the proposed method can effectively extract discriminative features for a motor imagery-based BMI.


---

## Body

## 1. Introduction

Brain machine/computer interfacing (BMI/BCI) is a challenging technology of signal processing, machine learning, and neuroscience [1]. BMIs capture brain activities associated with mental tasks and external stimuli, realize nonmuscular communication, and control channel for conveying messages and commands to the external world [1–3]. Basically, noninvasively measured data such as electroencephalogram (EEG), magnetoencephalogram (MEG), and functional magnetic resonance imaging (fMRI) are widely used to observe brain activities. Among them, because of its simplicity and low cost, EEG is practical for use in engineering applications [4, 5].Efficient decoding around motor-cortex is a crucial technique for the realization of BMI associated with motor-imagery (MI-BMI) [6, 7] with the application to controlling external devices [7], prostheses [4], rehabilitation [8], and so forth. For instance, it is also known that the real and imaginary movements of hands and feet evoke the change of the so-called mu rhythm in different brain regions [2, 3]. Therefore, the accurate extraction of these changes from the measured EEG signals in the presence of measurement noise and spontaneous components which are related to other brain activities enables us to classify the EEG signal associated with the different motor (imagined) actions such as movement of the right hand, left hand, or feet.In classification of EEG signals in MI-BMI and analyzing of the brain activities during motor imagery, signal processing techniques such as bandpass filtering and spatial weighting are used [1]. For the processing, presuming the parameters such as coefficients of the filters and weights that extract the related components is a crucial issue. Moreover, the optimum parameters in classification are highly dependent on users and measurement environments [9].In order to determine the parameters, data-driven techniques that exploit observed data are widely used [1, 2]. The observed data essentially include class labels corresponding to the tasks. The techniques should find the parameters that extract discriminative features as much as possible. For example, the well-known common spatial pattern (CSP) method finds the spatial weights by using the observed signals [1, 9, 10] in such a way that the variances of the signals extracted by the linear combination of a multichannel signal and the spatial weights differ as much as possible between two classes. The standard CSP method has been extended to methods to estimate the other parameters, such as the frequency bands [11–16], and methods to select the CSP features extracted with various parameters [17, 18].Besides, one of the parameters to be decided by data-driven techniques is a time window because of the following reasons. A kind of BMIs is implemented based on cues which a user follows. In the BMI, the user begins to perform a task when the cue is given. Therefore, the time when the user begins to perform the task is known. However, the time when the brain activity associated with the task occurs is unknown. The time windows working to remove samples that do not contain the brain activity will not match the period when the cues are showed. For instance, the samples for a few hundreds of milliseconds after the cues are assumed not to be used to extract the features in previous works [13, 17, 19, 20], which heuristically determined the time window. Contrary to these works, this paper hypothesizes that an optimal observation period in classification depends on users. For example, reaction time defined as the elapsed time between the presentation of a sensory stimulus and the subsequent behavioral response is strongly associated with age [21]. The reaction time can be related to the time of response to the cues. Therefore, the time window should also be designed by using the observed signals or data-driven.In this paper, we propose a method for finding the time windows as well as the aforementioned parameters (CSP and temporal filters) by extending a framework proposed in [15]. The proposed method enables us to find these parameters that separately extract several components that are observed in different spatial patterns, frequency bands, and periods of time. We call these simultaneously designed parameters common spatio-time-frequency patterns (CSTFP). As compared to the methods, in which the parameters are selected out of the predefined candidates fixed in advance [18], the CSTFP method has a higher degree of freedom for the parameters reducing the computational costs. In CSTFP, the coefficients of the temporal filters and the spatial weights are searched in the set of real numbers. Moreover, although the time windows are selected out of a set of candidates, the computational cost of the CSTFP method does not grow rapidly even with a large number of candidates.The rest of this paper is organized as follows. Section2.1 reviews the CSP method. Next we, illustrate the CSTFP method from Section 2.2 to Section 2.4. We experimentally analyze the CSTFP method by using artificial data in Section 3. Section 4 presents experimental results of classification of EEG signals during motor imageries to show the performance of CSTFP. Finally, the conclusions of this paper are presented in Section 5.


## 2. Common Spatio-Time-Frequency Patterns (CSTFP)

We describe a novel method called the CSTFP method to simultaneously find the parameters for the spatial weights, the temporal filters, and the time windows. In this section, before introducing CSTFP, CSP is reviewed. Next, we define a signal extraction model using the spatial weighting, filtering, and time windowing. Then, we propose a criterion designing the parameters motivated by that for the CSP method. The proposed criterion evaluates not only the spatial weights but also the coefficients of the time windows and the temporal filters. Next, the optimization method for the proposed criterion with alternate updating is presented. Finally, we define feature vector with the feature extraction model for classification of unlabeled EEG signals.


### 2.1. Common Spatial Pattern (CSP): A Review

LetX∈ℝM×N be an observed signal, where M is the number of channels and N is the number of samples. In BMI application, we do not directly use X, but we use the filtered signal described as X^=ℋ(X) to find the CSP, where ℋ is a bandpass filter which passes the frequency components related to brain activity of motor imagery. Denote the components of X^ by X^=[x^1,…,x^N], where x^n∈ℝM and n is the time index. We assume the sets of the observed signals, 𝒞1 and 𝒞2, where 𝒞d contains the signals belonging to class d, d∈{1,2} is a class label, 𝒞1∩𝒞2=∅, and ∅ is a set having no elements. CSP is defined as the weight vector that maximizes the intra class variance in 𝒞c under the normalization of samples, where c is a class label. More specifically, for c being fixed, the weight vector is found by solving the following optimization problem [9, 10]:
(1)maxwdddsdEX∈𝒞c[1N∑n=1N‍|wT(x^n-μ)|2],subjectto∑d=1,2‍EX∈𝒞d[1N∑n=1N‍|wT(x^n-μ)|2]=1,
where EX∈𝒞d[·] denotes the expectation over 𝒞d, μ is the time average of X given by μ=(1/N)∑n=1N‍x^n, ·T is the transpose of a vector or a matrix, and |·| is the absolute value of a scalar. The solution of (1) is given by the generalized eigenvector corresponding to the largest generalized eigenvalue of the generalized eigenvalue problem described as
(2)Σcw=λ(Σ1+Σ2)w,
where Σd, for d=1,2, are defined as Σd=EX∈𝒞d[(1/N)∑n=1N‍(x^n-μ)(x^n-μ)T].


### 2.2.  Signal Extraction Model

The target signal and signal extraction procedure are formulated in this section. The filtered signal of a target signal,X, denoted as x^=[x^1,…,x^K]T, is defined as
(3)x^n=bn∑m=1M‍wm∑p=1P‍hp[xn+P-p]m,
for n=1,…,K, K=N-P+1, where [·]i is the ith entry of a vector, P is the filter order of an FIR filter of which the coefficients are denoted by h1,…,hP, wm is a spatial weight for mth channel, and bn is a time window for nth sample that takes a binary value of either 0 or 1. The structure of temporally filtering, spatial weighting, and windowing for X is illustrated in Figure 1. In this model, wm is regarded as the spatial pattern, hp is regarded as the frequency pattern, and bn is regarded as the time pattern of the extracted signal.Figure 1
The filtering of an observed signal by a temporal filter,h, spatial weights, w, and a time window, b. Z-n is an operator of n samples delay.The sample variance ofx^n over time n=1,…,K is described as
(4)αX(w,h,b)=1∥b∥∑n=1K‍|x^n-1K∑i=1K‍x^i|2,
where w is defined as w=[w1,…,wM]T, h is defined as h=[h1,…,hP]T, b is defined as b=[b1,…,bN]T, and ∥·∥ is the Euclidean norm of a vector.The variance defined in (4) can be transformed to matrix-vector form as follows. We define An, n=1,…,K, whose elements are from X as
(5)[An]m,p=[X]m,n+P-p,
for m=1,…,M and p=1,…,P, where [·]i,j is the element at the ith row and the jth column of a matrix. Then, (4) can be modified to
(6)αX(w,h,b)=1∥b∥∑n=1K‍bn|wTA^nh|2,
where A^n is defined as A^n=An-∥b∥-1∑m=1K‍bmAm.


### 2.3.  Optimization for Sets of Parameters

We consider the problem of the design ofF sets of the FIR filter, the spatial weights, and the time window represented by {wi,hi,bi}i=1F. These sets of the parameters are designed in such a way that wi, hi, and bi maximize expectation of αX(wi,hi,bi) with respect to X∈𝒞c samples under the normalization of an expectation of αX(wi,hi,bi) over all of the observation. Additionally, we impose the orthonormality on hi, i=1,…,F to avoid the trivial solution. Moreover, the time windows are chosen from given candidates for efficient optimization. Therefore, we formulate the following maximization problem:
(7)max𝒫i,i=1,…,Fdd∑i=1F‍J^(𝒫i)+ϵK∥bi∥,subjecttohiThj∥hi∥∥hj∥=δij,i,j=1,…,F,ffffffffffffffhk∈𝒮k⊥,k=1,…,F,ffffffffffffffbl∈ℬ,l=1,…,F,
where 𝒫i represents set {wi,hi,bi} and J^(𝒫i) is the cost evaluating the ratio of the feature value in all samples defined as
(8)J^(𝒫i)=EX∈𝒞c[αX(wi,hi,bi)]∑d=1,2EX∈𝒞d[αX(wi,hi,bi)],where 𝒮i is any subspace in ℝP, ℬ is a candidate set for the time windows, defined as ℬ={bl}l=1L, c is a class label chosen from 1 and 2, ϵ is a regularization parameter, and δij is the Kronecker delta defined as 1 for i=j and 0, otherwise.Since it is difficult to simultaneously find all parameters, we consider sequential optimization to find the parameters with respect to each filter indexi. That is, we first find 𝒫1, and then find 𝒫2 under the constraint on h1. This sequential optimization is represented with respect to each i as
(9)max𝒫iddddJ^(𝒫i)+ϵK∥bi∥,subjecttohi∈𝒮^i⊥,b∈ℬ,
where 𝒮^i is a subspace defined as
(10)𝒮^i=Span{h1,…,hi-1}⊕𝒮i∈ℝP,
where Span(⋯) represents a subspace spanned by vectors and the operator denoted by ⊕ gives the direct sum of two subspaces. Methods for choosing 𝒮^i has been discussed in [15]. In (9), to optimize the parameters indexed with i, we adopt an alternating optimization procedure based on alternating least square (ALS). In the optimization, we separate the problem of (9) into three subproblems for wi, hi, and bi, respectively. Then, we update the parameters by alternating solving the subproblems. The three subproblems and these solutions are as follows.The first subproblem is to optimizewi. While fixing hi and bi, wi maximizing (9) is found as the generalized eigenvector corresponding to the largest generalized eigenvalue of the generalized eigenvalue problem [15] described as
(11)Rcwi=λ(R1+R2)wi,
where(12)Rd=EX∈𝒞d[1∥bi∥∑n=1K‍[b]nA^nhihiTA^nT],
for d=1,2, and λ is an eigenvalue.The second subprolem is to optimizehi. While fixing wi and bi, hi maximizing (9) is found as the generalized eigenvector corresponding to the largest generalized eigenvalue of the generalized eigenvalue problem [15] described as
(13)GQchi=ζ(Q1+Q2)hi,
where
(14)Qd=EX∈𝒞d[1∥bi∥∑n=1K‍[bi]nA^nTwiwiTA^n],
for d=1,2,
(15)G=IP-V(VT(Q1+Q2)-1V)-1VT(Q1+Q2)-1,V=[vi1,…,viDi]∈ℝP×Di,vi1,…,viDi are vectors spanning S^i described as
(16)Span{vi1,…,viDi}=𝒮^i,where IP is the P×P identify matrix and ζ is an eigenvalue.The third subproblem is to optimizebi while fixing wi and hi. The cost of (9) can be reduced to
(17)J3(bi∣wi,hi)=g1(bi)g1(bi)+g2(bi)+ϵK∥bi∥.
Then,
(18)gd(b)=EX∈𝒞d[1∥b∥∑n=1K‍(bnx~n-1∥b∥∑k=1K‍bkx~k)2]=1∥b∥EX∈𝒞d[bT(ξ+2μx~)+μ],
where x~n=wiTAnhi, x~=[x~1,…,x~K]T, ξ=[x~12,…,x~K2]T, and μ=∥b∥-1bTx~. Because bi are chosen out of ℬ, we calculate the values J3(b^i) for all candidates in ℬ and the optimal bi can be chosen as the candidates that maximize J3. This formulates that
(19)bi=arg maxb∈ℬJ3(b).The procedure to design the spatial weights, the temporal filters, and the time windows is summarized in Algorithm1 as a pseudocode.Algorithm 1:Design of the FIR filters, the spatial weights, and the time windows.
Input:  𝒞1 and 𝒞2: the sets of observed signals.
Parameters:  F:  the number of FIR filters, P: the filter order,
ℬ: the set of the candidates for the time window.
Output:  {wi}i=1F: the spatial weights, {hi}i=1F: the vectors of
FIR filter coefficients,{bi}i=1F: the time window.
for  i=1,…,F  do
Initializehi and bi.
Set the index of iteration ask=0.
repeat
k←k+1.
Updatewi by solving (11).
Updatehi by solving (13).
Updatebi by solving (19).
Calculate cost,Ck from the cost function, J^(𝒫i).
untilCk-Ck-1 is sufficiently small.
end for


### 2.4.  Feature Vector Definition

With designedwi, hi, and bi, the feature vector of an EEG signal, X, for classification is defined as
(20)y=[αX(w^1(1),h1,b1),…,αX(w^1(2r),h1,b1),…,dαX(w^F(1),hF,bF),…,αX(w^F(2r),hF,bF)]T,
where the set of {w^i(1),…,w^i(2r)} is the CSPs corresponding to hi and bi. w^i(m), i=1,…,F, m=1,…,2r, are decided as follows. By solving (11) with hi and bi, we obtain MF spatial patterns as w~i(m) for i=1,…,F and m=1,…,M, where w~i(m) is the unit-length eigenvector corresponding to the mth largest eigenvalue of (11). Then, w^i(m) are defined as w^i(j)=w~i(j) and w^i(r+j)=w~i(M-j+1) for i=1,…,F, j=1,…,r.


## 2.1. Common Spatial Pattern (CSP): A Review

LetX∈ℝM×N be an observed signal, where M is the number of channels and N is the number of samples. In BMI application, we do not directly use X, but we use the filtered signal described as X^=ℋ(X) to find the CSP, where ℋ is a bandpass filter which passes the frequency components related to brain activity of motor imagery. Denote the components of X^ by X^=[x^1,…,x^N], where x^n∈ℝM and n is the time index. We assume the sets of the observed signals, 𝒞1 and 𝒞2, where 𝒞d contains the signals belonging to class d, d∈{1,2} is a class label, 𝒞1∩𝒞2=∅, and ∅ is a set having no elements. CSP is defined as the weight vector that maximizes the intra class variance in 𝒞c under the normalization of samples, where c is a class label. More specifically, for c being fixed, the weight vector is found by solving the following optimization problem [9, 10]:
(1)maxwdddsdEX∈𝒞c[1N∑n=1N‍|wT(x^n-μ)|2],subjectto∑d=1,2‍EX∈𝒞d[1N∑n=1N‍|wT(x^n-μ)|2]=1,
where EX∈𝒞d[·] denotes the expectation over 𝒞d, μ is the time average of X given by μ=(1/N)∑n=1N‍x^n, ·T is the transpose of a vector or a matrix, and |·| is the absolute value of a scalar. The solution of (1) is given by the generalized eigenvector corresponding to the largest generalized eigenvalue of the generalized eigenvalue problem described as
(2)Σcw=λ(Σ1+Σ2)w,
where Σd, for d=1,2, are defined as Σd=EX∈𝒞d[(1/N)∑n=1N‍(x^n-μ)(x^n-μ)T].


## 2.2.  Signal Extraction Model

The target signal and signal extraction procedure are formulated in this section. The filtered signal of a target signal,X, denoted as x^=[x^1,…,x^K]T, is defined as
(3)x^n=bn∑m=1M‍wm∑p=1P‍hp[xn+P-p]m,
for n=1,…,K, K=N-P+1, where [·]i is the ith entry of a vector, P is the filter order of an FIR filter of which the coefficients are denoted by h1,…,hP, wm is a spatial weight for mth channel, and bn is a time window for nth sample that takes a binary value of either 0 or 1. The structure of temporally filtering, spatial weighting, and windowing for X is illustrated in Figure 1. In this model, wm is regarded as the spatial pattern, hp is regarded as the frequency pattern, and bn is regarded as the time pattern of the extracted signal.Figure 1
The filtering of an observed signal by a temporal filter,h, spatial weights, w, and a time window, b. Z-n is an operator of n samples delay.The sample variance ofx^n over time n=1,…,K is described as
(4)αX(w,h,b)=1∥b∥∑n=1K‍|x^n-1K∑i=1K‍x^i|2,
where w is defined as w=[w1,…,wM]T, h is defined as h=[h1,…,hP]T, b is defined as b=[b1,…,bN]T, and ∥·∥ is the Euclidean norm of a vector.The variance defined in (4) can be transformed to matrix-vector form as follows. We define An, n=1,…,K, whose elements are from X as
(5)[An]m,p=[X]m,n+P-p,
for m=1,…,M and p=1,…,P, where [·]i,j is the element at the ith row and the jth column of a matrix. Then, (4) can be modified to
(6)αX(w,h,b)=1∥b∥∑n=1K‍bn|wTA^nh|2,
where A^n is defined as A^n=An-∥b∥-1∑m=1K‍bmAm.


## 2.3.  Optimization for Sets of Parameters

We consider the problem of the design ofF sets of the FIR filter, the spatial weights, and the time window represented by {wi,hi,bi}i=1F. These sets of the parameters are designed in such a way that wi, hi, and bi maximize expectation of αX(wi,hi,bi) with respect to X∈𝒞c samples under the normalization of an expectation of αX(wi,hi,bi) over all of the observation. Additionally, we impose the orthonormality on hi, i=1,…,F to avoid the trivial solution. Moreover, the time windows are chosen from given candidates for efficient optimization. Therefore, we formulate the following maximization problem:
(7)max𝒫i,i=1,…,Fdd∑i=1F‍J^(𝒫i)+ϵK∥bi∥,subjecttohiThj∥hi∥∥hj∥=δij,i,j=1,…,F,ffffffffffffffhk∈𝒮k⊥,k=1,…,F,ffffffffffffffbl∈ℬ,l=1,…,F,
where 𝒫i represents set {wi,hi,bi} and J^(𝒫i) is the cost evaluating the ratio of the feature value in all samples defined as
(8)J^(𝒫i)=EX∈𝒞c[αX(wi,hi,bi)]∑d=1,2EX∈𝒞d[αX(wi,hi,bi)],where 𝒮i is any subspace in ℝP, ℬ is a candidate set for the time windows, defined as ℬ={bl}l=1L, c is a class label chosen from 1 and 2, ϵ is a regularization parameter, and δij is the Kronecker delta defined as 1 for i=j and 0, otherwise.Since it is difficult to simultaneously find all parameters, we consider sequential optimization to find the parameters with respect to each filter indexi. That is, we first find 𝒫1, and then find 𝒫2 under the constraint on h1. This sequential optimization is represented with respect to each i as
(9)max𝒫iddddJ^(𝒫i)+ϵK∥bi∥,subjecttohi∈𝒮^i⊥,b∈ℬ,
where 𝒮^i is a subspace defined as
(10)𝒮^i=Span{h1,…,hi-1}⊕𝒮i∈ℝP,
where Span(⋯) represents a subspace spanned by vectors and the operator denoted by ⊕ gives the direct sum of two subspaces. Methods for choosing 𝒮^i has been discussed in [15]. In (9), to optimize the parameters indexed with i, we adopt an alternating optimization procedure based on alternating least square (ALS). In the optimization, we separate the problem of (9) into three subproblems for wi, hi, and bi, respectively. Then, we update the parameters by alternating solving the subproblems. The three subproblems and these solutions are as follows.The first subproblem is to optimizewi. While fixing hi and bi, wi maximizing (9) is found as the generalized eigenvector corresponding to the largest generalized eigenvalue of the generalized eigenvalue problem [15] described as
(11)Rcwi=λ(R1+R2)wi,
where(12)Rd=EX∈𝒞d[1∥bi∥∑n=1K‍[b]nA^nhihiTA^nT],
for d=1,2, and λ is an eigenvalue.The second subprolem is to optimizehi. While fixing wi and bi, hi maximizing (9) is found as the generalized eigenvector corresponding to the largest generalized eigenvalue of the generalized eigenvalue problem [15] described as
(13)GQchi=ζ(Q1+Q2)hi,
where
(14)Qd=EX∈𝒞d[1∥bi∥∑n=1K‍[bi]nA^nTwiwiTA^n],
for d=1,2,
(15)G=IP-V(VT(Q1+Q2)-1V)-1VT(Q1+Q2)-1,V=[vi1,…,viDi]∈ℝP×Di,vi1,…,viDi are vectors spanning S^i described as
(16)Span{vi1,…,viDi}=𝒮^i,where IP is the P×P identify matrix and ζ is an eigenvalue.The third subproblem is to optimizebi while fixing wi and hi. The cost of (9) can be reduced to
(17)J3(bi∣wi,hi)=g1(bi)g1(bi)+g2(bi)+ϵK∥bi∥.
Then,
(18)gd(b)=EX∈𝒞d[1∥b∥∑n=1K‍(bnx~n-1∥b∥∑k=1K‍bkx~k)2]=1∥b∥EX∈𝒞d[bT(ξ+2μx~)+μ],
where x~n=wiTAnhi, x~=[x~1,…,x~K]T, ξ=[x~12,…,x~K2]T, and μ=∥b∥-1bTx~. Because bi are chosen out of ℬ, we calculate the values J3(b^i) for all candidates in ℬ and the optimal bi can be chosen as the candidates that maximize J3. This formulates that
(19)bi=arg maxb∈ℬJ3(b).The procedure to design the spatial weights, the temporal filters, and the time windows is summarized in Algorithm1 as a pseudocode.Algorithm 1:Design of the FIR filters, the spatial weights, and the time windows.
Input:  𝒞1 and 𝒞2: the sets of observed signals.
Parameters:  F:  the number of FIR filters, P: the filter order,
ℬ: the set of the candidates for the time window.
Output:  {wi}i=1F: the spatial weights, {hi}i=1F: the vectors of
FIR filter coefficients,{bi}i=1F: the time window.
for  i=1,…,F  do
Initializehi and bi.
Set the index of iteration ask=0.
repeat
k←k+1.
Updatewi by solving (11).
Updatehi by solving (13).
Updatebi by solving (19).
Calculate cost,Ck from the cost function, J^(𝒫i).
untilCk-Ck-1 is sufficiently small.
end for


## 2.4.  Feature Vector Definition

With designedwi, hi, and bi, the feature vector of an EEG signal, X, for classification is defined as
(20)y=[αX(w^1(1),h1,b1),…,αX(w^1(2r),h1,b1),…,dαX(w^F(1),hF,bF),…,αX(w^F(2r),hF,bF)]T,
where the set of {w^i(1),…,w^i(2r)} is the CSPs corresponding to hi and bi. w^i(m), i=1,…,F, m=1,…,2r, are decided as follows. By solving (11) with hi and bi, we obtain MF spatial patterns as w~i(m) for i=1,…,F and m=1,…,M, where w~i(m) is the unit-length eigenvector corresponding to the mth largest eigenvalue of (11). Then, w^i(m) are defined as w^i(j)=w~i(j) and w^i(r+j)=w~i(M-j+1) for i=1,…,F, j=1,…,r.


## 3. Experimental Analysis of Artificial Signal

We give an analysis of CSTFP by a toy experiment with an artificial signal in this section. We assume a 2-class BMI where the observed EEG signals are modeled by a mixture of narrow-band signals (see Figure2). In this model, a trial signal belonging to class d is given by
(21)x[n]=∑i=1Ns‍si[n]ai(d)+η,n=0,…,N-1,
where d is a class label taking either 1 or 2, x[n]∈ℝM is a vector representing a signal at a discrete time point n, N is the number of time samples for a trial, M is the number of channels, si[n]∈ℝ is an ith source signal of feature components, Ns is the number of the source signals, ai(d)∈ℝM is a vector defined as ai(d)=[ai1(d),…,aiM(d)]T, aim(d)∈ℝ is the amplitude of si[n] in the mth channel for class d, and η∈ℝM is a stochastic noise. The source signals, si[n], are generated as
(22)si[n]=ti[n]ℜ[∑k=0N-1‍Si[k]ejθej(2πk/N)n],
where Si[k]∈ℝ represents a discrete spectrum, ti[n] represents the time window that decides the period when the ith source signal is generated, θ∈ℝ is a stochastic phase of the source signals, and the operator denoted by ℜ takes a real part of a complex number.Frequency patterns of the spectra of the source signals ((a)–(d)) and the amplitude characteristics of the designed frequency patterns ((e)–(h)). The horizontal axis represents frequency and its unit is Hz.
(a)
1st source,s1
(b)
2nd source,s2
(c)
3rd source,s3
(d)
4th source,s4
(e)
1st filter,h1
(f)
2nd filter,h2
(g)
3rd filter,h3
(h)
4th filter,h4In particular, 200 artificial signals that we used in the experiment were generated with the conditions shown in Table1 where 𝒩(m,σ2) is a Gaussian distribution with a mean, m, and a variance, σ2, and 𝒰(a,b) is a uniform distribution whose minimum and maximum values are denoted by a and b, respectively.Table 1
The conditions for generating the artificial signals.
Parameter
Value or distributions
Number of channel,M
10
Number of samples,N
100
Number of trials for each class
100
Sampling frequency
100 Hz
Number of sources,Ns
4
Spectra of sources,|Si[k]|
Figures2(a)–2(d)
Stochastic phase,θ
𝒰
(
0,2
π
)
Time windows,ti[n]
Figures3(a)–3(d)
Amplitudes,ai(d)
Figures4(a)–4(d)
Stochastic noise,[η]m
𝒩
(
0,0.1
)We applied the CSTFP method to the artificial signals as follows. The class label represented byc in (8) was set to 1. The number of the sets of the parameters, F, was set to 4. The order of the temporal filters was set to 41, yielding that P=41; the length of a filtered signal, K, is 60. For the given candidates for the time windows, we define the following set. First, we define ten K-dimensional vectors as
(23)dj=[0,…,0︸D(j-1),1,…,1︸D,0,…,0]T,j=1,…,10,
where D=6. Then, we use all combinations of {dj}j=110 represented by
(24)b^l=p1d1+⋯+p10d10
and pj∈{0,1}, j=1,…,10, as the given candidate set. Therefore, the number of the candidates, L, was 1023. Moreover, for 𝒮^i that determines the search space for hi, we used 𝒮^i=Span{h1,…,hi-1,hi-1′,…,hi-1′} [15], where each element of hj′ is defined as [hj′]k=[hj]((k+1)modP), k=1,…,P, and the operator denoted by amodb takes a residue of dividing a by b. In addition, the regularization parameter, ϵ, was set to 0.05.The optimization resulted in Figures2(e)–2(h) for the amplitude characteristics of the FIR filters, Figures 3(e)–3(h) for the time windows, and Figures 4(e)–4(h) for the normalized spatial weights. The centers of passbands of the filters shown in Figures 2(e)–2(h) coincide with the centers of the source signals shown in Figures 2(a)–2(d). Moreover, in the spectrum of the source and the amplitude characteristic of the designed FIR filter that have similar center frequencies, the spatial amplitude corresponding to the source and the spatial weight vector corresponding to the FIR filter are similar to each other. For instance, the correlation coefficient between a1(1) shown in Figure 4(a) with circles and w2(1) shown in Figure 4(f) with circles is 0.971.Time patterns as the time windows of the source signals ((a)–(d)) and the designed time windows ((e)–(h)). The horizontal axis represents time and its unit is the second.
(a)
1st source,s1
(b)
2nd source,s2
(c)
3rd source,s3
(d)
4th source,s4
(e)
1st window,b1
(f)
2nd window,b2
(g)
3rd window,b3
(h)
4th window,b4Spatial patterns as the spatial amplitudes of the source signals ((a)–(d)) and the designed spatial weights ((e)–(h)). The amplitudes plotted by circles areai(1) and the amplitudes plotted by cross are ai(2) for i=1,…,4. The weights plotted by circles are wi(1) and the weights plotted by cross are wi(10) for i=1,…,4. The horizontal axis represents the channel number.
(a)
1st source,a1(1) and a1(2)
(b)
2nd source,a2(1) and a2(2)
(c)
3rd source,a3(1) and a3(2)
(d)
4th source,a4(1) and a4(2)
(e)
1st weights,w^1(1) and w^1(10)
(f)
2nd weights,w^2(1) and w^2(10)
(g)
3rd weights,w^3(1) and w^3(10)
(h)
4th weights,w^4(1) and w^4(10)The results also suggest that the time windows designed by the CSTFP method can remove samples observed in the periods of time that do not contain the source signals. For instance, hence, the source signal,s1[n], is not observed in the first 25 samples according to Figure 3(a), the time window for extracting s1[n] is expected to remove the first 25 samples. Because |S1[k]| and the amplitude characteristics shown in Figure 2(f) have similar center frequency, we can decide that the time window for extracting s1[n] is the time window shown in Figure 3(f). Although the designed time windows do not coincide with ti[n] because b is applied to an FIR-filtered signal that is shorter than that of the original one, the time window removes samples in the first 10 samples, as we expected. Moreover, we can decide that the time window shown in Figure 3(c) is for extracting s3[n] due to the same reason. As the observed signals do not have s3[n] in the last 25 samples (see Figure 3(c)), the time window shown in Figure 3(c) removes the samples in the last 20 samples.


## 4. Experiment of EEG Signal Classification

A comprehensive comparative study was performed to illustrate the ability of the CSTFP method to produce more accurate classification of EEG signals during motor imagery over several conventional methods (CSP [10], common sparse spectral spatial patterns (CSSSP) [12], filter bank CSP (FBCSP) [17], and discriminative filter bank CSP (DFBCSP) [15]).


### 4.1. Data Description

We used dataset IVa from BCI competition III [22], which was provided by Fraunhofer FIRST (Intelligent Data Analysis Group) and Campus Benjamin Franklin of the Charité - University Medicine Berlin (Department of Neurology, Neurophysics Group) [23] and dataset 1 from BCI competition IV, which was provided by Berlin Institute of Technology (Machine Learning Laboratory), Fraunhofer FIRST (Intelligent Data Analysis Group), and Campus Benjamin Franklin of the Charité - University Medicine Berlin (Department of Neurology, Neurophysics Group) [24]. The condition for each dataset is shown in Table 2. They have two classes of motor imagery. The signals in the provided datasets were recorded with the sampling rate of 1000 Hz.Table 2
Description of the datasets.
Dataset IVa
Dataset 1
Classes
Right hand and right foot
2 tasks from foot, left hand, and right hand
Subject labels
a
a
,
a
l
,
a
v
,
a
w
,
a
y
a
,
b
,
f
,
g
Number of channels
118
59
Signal length
3.5 secs
4 secs
Sampling rate
100 Hz
100 Hz
Number of the trials per class
140
100We furthermore applied to this dataset a Butterworth low-pass filter whose cutoff frequency is 50 Hz and the filter order is 4, and downsampled to 100 Hz.


### 4.2.  Result

For the experiments, as a sample for each trial, we used a signal observed in the period fromT1 to T2 [sec] after the cue that directs the subject to perform the task. In the experiments, T1 was tuned by a method we mention later. T2 was set to 3.5 and 4 seconds for BCI competition III dataset IVa and IV dataset 1, respectively.In order to compare the classification abilities for the methods, we obtained the classification accuracy rates by5×5 cross-validation (CV). In each classification in the CV, we separated learning samples for selecting the parameters of the feature extraction and a linear discriminant analysis (LDA) classifier and test samples for obtaining classification accuracy rates.For the methods to be compared (CSP, CSP-Exh, CSSSP, FBCSP, DFBCSP, and CSTFP), the parameters for the feature extraction were obtained as follows.(i)
CSP: the parameters determined in this method are spatial weights. Before obtaining the spatial weights by the CSP method, we applied the Butterworth bandpass filter with the passband of 7–30 Hz. In the CSP method, we minimized the variance cost of the right hand class in (1). The eigenvectors corresponding to the r largest and r smallest eigenvalues of the eigenvalue problem (2) were given as the spatial weights.
(ii)
CSP-Exh: the parameters determined in this method are spatial weights and a passband of the Butterworth filter. The passband of the Butterworth filters was tuned asfl-fu Hz by an exhaustive search by the CSP method and the learning samples. After the filtering with the passband, the spatial weights were given by the same manner as that of the CSP method.
(iii)
CSSSP: the parameters determined in this method are spatial weights and a bandpass filter. The bandpass filter between 7 and 30 Hz was applied as preprocessing [12]. CSSSP was applied with regularization parameters, C, and the parameter for the number of the spatial weights, r. The order of the filter was fixed to 16 [12].
(iv)
FBCSP: the parameters determined in this method arer bandpass filters out of a filter bank and associated spatial weights. FBCSP was applied with the mutual information based best individual feature and a naïve Bayesian Parzen window (NBPW) classifier [17]. The filter bank comprising 9 bandpass filters covering 4–40 Hz was used. All filters were Chebyshev type II filters with a bandwidth of 4 Hz each. In FBCSP, the number of the spatial weights, NM, in each band was set to 8. These parameters were decided by referring to [17].
(v)
DFBCSP: the parameters determined in this method areF FIR filters and r spatial weights associated to each FIR filter. DFBCSP was applied with the FIR filter order of 41 as done in [15]. In optimization, we stopped iteration when error of the cost function between successive iterations becomes under 10-5.
(vi)
CSTFP: the parameters determined in this method areF FIR filters, the corresponding F time windows, and r spatial weights associated with each FIR filter. We fixed T1=0 to observe the behavior of the resulting time window. CSTFP is applied with the FIR filter order, 41 [15], and the following candidate set for the time windows. The candidate set for the time windows, ℬ, consists of the vectors defined as
(25)b^l=[0,…,0︸D,1,…,1︸O,0,…,0]T,
where we choose D and O out of a set {0,5,10,…,K} such that O>50 and D+O≤K. The regularization parameter, ϵ, was set to 0.1. In alternating optimization, we initialize hi as a random vector which is orthonormalized from v1,…,vDi and bi as a vector all the elements of which are one. We stopped iteration when error of the cost function between successive iterations becomes under 10-5.For the obscure parameters such asr in the above list, we furthermore tuned them by 5 × 5 CV using the learning samples as done in [25]. We conducted the nested CV [25] with all of the combinations of these parameters and obtained the classification accuracy rates. We adopted the combination that performed the highest rates as the parameters. The parameters tuned by the nested CV in the learning data and the candidates for them are summarized in Table 3.Table 3
The parameters decided by the nested CV in the learning data in the classification experiments.
Method
Parameters and candidates
CSP
T
1
∈
{
0,0.25,0.5,0.75
}
r
∈
{
1,2
,
…
,
10
}
CSP-Exh
T
1
∈
{
0,0.25,0.5,0.75
}
r
∈
{
1,2
,
…
,
10
}
f
l
∈
{
1,2
,
…
,
48
}
f
u
∈
{
f
l
+
1
,
f
l
+
2
,
…
,
49
}
CSSSP
T
1
∈
{
0,0.25,0.5,0.75
}
r
∈
{
1,2
,
…
,
10
}
C
∈
{
0,0.01,0.1,0.2,0.5,1
,
2,5
}
FBCSP
T
1
∈
{
0,0.25,0.5,0.75
}
r
∈
{
1,2
,
…
,
10
}
DFBCSP
T
1
∈
{
0,0.25,0.5,0.75
}
r
∈
{
1,2
,
…
,
10
}
F
∈
{
1,2
,
…
,
5
}
CSTFP
r
∈
{
1,2
,
…
,
10
}
F
∈
{
1,2
,
…
,
5
}After we obtained the feature vectors extracted by the filters, spatial weights, and the time windows that are designed by each listed method, we calculated the logarithm of the feature vectors. Then, the LDA for the learning samples was used to obtain a projector onto the 1-dimensional space. The threshold for classification was determined as the middle point of two class averages over the projected learning samples. The feature vectors from the test samples were classified by the projection and the threshold and we obtained the classification accuracy rate in each CV. In Tables4 and 5, CSTFP results in the highest accuracy rate in the average over whole subjects.Table 4
Classification accuracy [%] given by5×5 CV in dataset IVa from BCI competition III. The figure with ± represents the standard deviation (S.D.) over CV.
Method
Subject
Average
aa
al
av
aw
ay
CSP
79.9
±
4.6
98.4
±
1.6
71.9
±
3.5
96.8
±
2.6
92.5
±
3.1
87.9
CSP-Exh
90.8
±
3.6
99.1
±
1.4
74.4
±
5.5
99.1
±
1.4
93.0
±
4.2
91.3
CSSSP
91.5
±
3.7
99.1
±
1.6
71.1
±
6.9
98.8
±
1.4
93.1
±
3.5
90.7
FBCSP
91.3
±
3.4
99.3
±
1.5
70.7
±
7.4
98.2
±
1.6
88.4
±
4.8
89.6
DFBCSP
91.9
±
2.6
98.9
±
1.5
74.8
±
4.0
98.9
±
1.6
96.0
±
3.1
92.1
CSTFP
92.6
±
2.1
98.9
±
1.5
75.4
±
5.3
99.0
±
1.3
96.0
±
2.4
92.5Table 5
Classification accuracy [%] given by5×5 CV in dataset 1 from BCI competition IV. The figure with ± represents the standard deviation (S.D.) over CV.
Method
Subject
Average
a
b
f
g
CSP
89.6
±
5.1
68.8
±
6.7
78.1
±
6.5
93.9
±
3.4
82.6
CSP-Exh
92.3
±
2.5
85.3
±
6.9
78.1
±
6.5
93.9
±
3.6
90.2
CSSSP
89.8
±
5.5
70.4
±
9.5
85.3
±
5.3
93.4
±
4.1
84.7
FBCSP
87.3
±
8.6
63.9
±
14.1
78.6
±
6.9
94.1
±
4.1
81.0
DFBCSP
89.7
±
5.0
85.6
±
6.9
92.5
±
4.4
93.7
±
4.3
90.4
CSTFP
91.4
±
4.0
90.6
±
7.2
93.3
±
3.6
93.5
±
4.1
92.2Moreover, we conducted the classification experiments, in which the parameters shown in Table3 were fixed to show the effects of the changes of the parameters on the classification accuracy rates. Figure 5 shows the accuracy rates of CSP and the rates of each method averaged over all subjects with each T1. For each subject, the parameters shown in Table 3 except for T1 were fixed to the combinations of the parameters that performed the highest accuracy rates with each T1. In Figure 5(a), T1 that perform the highest accuracy rates are different among the subjects. Moreover, Figure 5(b) shows that the accuracy rates highly depend on T1 in the conventional methods. Figure 6 shows the variation of the accuracy rates by the various regularization parameters, ϵ, in CSTFP. For each subject, the parameters, r and F, were fixed to the combinations of the parameters that performed the highest accuracy rates with each ϵ. CSTFP performs higher accuracy as ϵ goes higher than 0.1.The variation of the accuracy rates by variousT1.
(a)
CSP
(b)
Average over all subjectsThe variation of the accuracy rates by various regularization parameters,ϵ.
(a)
BCI competition III dataset IVa
(b)
BCI competition IV dataset 1We show examples of the spatial patterns, the frequency patterns (the amplitude characteristics of the FIR filters), and the time patterns designed by CSTFP in Figures7 and 8. As we can observe in Figures 6, 7(a), and 8(a), the short time windows caused by small ϵ result in poor classification accuracy rates. All of the time windows shown in Figures 7 and 8 remove the samples observed within a few hundreds of milliseconds after the cue. The results suggest that the brain activities related to the task cannot be observed just after the cue. In Figures 7(c)–7(e) and 8(c)–8(e), the time windows do not significantly change in around 0.1–0.4 of ϵ. This result of the relations between ϵ and the time windows corresponds with the result shown in Figure 6 in which the classification accuracy rates do not strongly depend on the regularization parameter, ϵ, by more than 0.12. Moreover, the FIR filters do not highly differ even with various ϵ. However, the spatial weights with the small windows (Figures 7(a) and 8(a)) differ from the other weights with the longer windows. Furthermore, there are slight differences in the time patterns and the frequency patterns between subject aa and subject av.Examples of spatio-time-frequency patterns as the spatial weights,w^1(1) (top), the amplitude response of the filter, h1,…,h4 (center), and the time windows, b1,…,b4 (bottom), designed by the CSTFP method with various ϵ in the result of subject aa. The vertical axis represents normalized amplitude.
(a)
ϵ
=
0.04
(b)
ϵ
=
0.08
(c)
ϵ
=
0.12
(d)
ϵ
=
0.16
(e)
ϵ
=
0.4Examples of spatio-time-frequency patterns as the spatial weights,w^1(1) (top), the amplitude response of the filter, h1,…,h4 (center), and the time windows, b1,…,b4 (bottom), designed by the CSTFP method with various ϵ in the result of subject av. The vertical axis represents normalized amplitude.
(a)
ϵ
=
0.04
(b)
ϵ
=
0.08
(c)
ϵ
=
0.12
(d)
ϵ
=
0.16
(e)
ϵ
=
0.4


## 4.1. Data Description

We used dataset IVa from BCI competition III [22], which was provided by Fraunhofer FIRST (Intelligent Data Analysis Group) and Campus Benjamin Franklin of the Charité - University Medicine Berlin (Department of Neurology, Neurophysics Group) [23] and dataset 1 from BCI competition IV, which was provided by Berlin Institute of Technology (Machine Learning Laboratory), Fraunhofer FIRST (Intelligent Data Analysis Group), and Campus Benjamin Franklin of the Charité - University Medicine Berlin (Department of Neurology, Neurophysics Group) [24]. The condition for each dataset is shown in Table 2. They have two classes of motor imagery. The signals in the provided datasets were recorded with the sampling rate of 1000 Hz.Table 2
Description of the datasets.
Dataset IVa
Dataset 1
Classes
Right hand and right foot
2 tasks from foot, left hand, and right hand
Subject labels
a
a
,
a
l
,
a
v
,
a
w
,
a
y
a
,
b
,
f
,
g
Number of channels
118
59
Signal length
3.5 secs
4 secs
Sampling rate
100 Hz
100 Hz
Number of the trials per class
140
100We furthermore applied to this dataset a Butterworth low-pass filter whose cutoff frequency is 50 Hz and the filter order is 4, and downsampled to 100 Hz.


## 4.2.  Result

For the experiments, as a sample for each trial, we used a signal observed in the period fromT1 to T2 [sec] after the cue that directs the subject to perform the task. In the experiments, T1 was tuned by a method we mention later. T2 was set to 3.5 and 4 seconds for BCI competition III dataset IVa and IV dataset 1, respectively.In order to compare the classification abilities for the methods, we obtained the classification accuracy rates by5×5 cross-validation (CV). In each classification in the CV, we separated learning samples for selecting the parameters of the feature extraction and a linear discriminant analysis (LDA) classifier and test samples for obtaining classification accuracy rates.For the methods to be compared (CSP, CSP-Exh, CSSSP, FBCSP, DFBCSP, and CSTFP), the parameters for the feature extraction were obtained as follows.(i)
CSP: the parameters determined in this method are spatial weights. Before obtaining the spatial weights by the CSP method, we applied the Butterworth bandpass filter with the passband of 7–30 Hz. In the CSP method, we minimized the variance cost of the right hand class in (1). The eigenvectors corresponding to the r largest and r smallest eigenvalues of the eigenvalue problem (2) were given as the spatial weights.
(ii)
CSP-Exh: the parameters determined in this method are spatial weights and a passband of the Butterworth filter. The passband of the Butterworth filters was tuned asfl-fu Hz by an exhaustive search by the CSP method and the learning samples. After the filtering with the passband, the spatial weights were given by the same manner as that of the CSP method.
(iii)
CSSSP: the parameters determined in this method are spatial weights and a bandpass filter. The bandpass filter between 7 and 30 Hz was applied as preprocessing [12]. CSSSP was applied with regularization parameters, C, and the parameter for the number of the spatial weights, r. The order of the filter was fixed to 16 [12].
(iv)
FBCSP: the parameters determined in this method arer bandpass filters out of a filter bank and associated spatial weights. FBCSP was applied with the mutual information based best individual feature and a naïve Bayesian Parzen window (NBPW) classifier [17]. The filter bank comprising 9 bandpass filters covering 4–40 Hz was used. All filters were Chebyshev type II filters with a bandwidth of 4 Hz each. In FBCSP, the number of the spatial weights, NM, in each band was set to 8. These parameters were decided by referring to [17].
(v)
DFBCSP: the parameters determined in this method areF FIR filters and r spatial weights associated to each FIR filter. DFBCSP was applied with the FIR filter order of 41 as done in [15]. In optimization, we stopped iteration when error of the cost function between successive iterations becomes under 10-5.
(vi)
CSTFP: the parameters determined in this method areF FIR filters, the corresponding F time windows, and r spatial weights associated with each FIR filter. We fixed T1=0 to observe the behavior of the resulting time window. CSTFP is applied with the FIR filter order, 41 [15], and the following candidate set for the time windows. The candidate set for the time windows, ℬ, consists of the vectors defined as
(25)b^l=[0,…,0︸D,1,…,1︸O,0,…,0]T,
where we choose D and O out of a set {0,5,10,…,K} such that O>50 and D+O≤K. The regularization parameter, ϵ, was set to 0.1. In alternating optimization, we initialize hi as a random vector which is orthonormalized from v1,…,vDi and bi as a vector all the elements of which are one. We stopped iteration when error of the cost function between successive iterations becomes under 10-5.For the obscure parameters such asr in the above list, we furthermore tuned them by 5 × 5 CV using the learning samples as done in [25]. We conducted the nested CV [25] with all of the combinations of these parameters and obtained the classification accuracy rates. We adopted the combination that performed the highest rates as the parameters. The parameters tuned by the nested CV in the learning data and the candidates for them are summarized in Table 3.Table 3
The parameters decided by the nested CV in the learning data in the classification experiments.
Method
Parameters and candidates
CSP
T
1
∈
{
0,0.25,0.5,0.75
}
r
∈
{
1,2
,
…
,
10
}
CSP-Exh
T
1
∈
{
0,0.25,0.5,0.75
}
r
∈
{
1,2
,
…
,
10
}
f
l
∈
{
1,2
,
…
,
48
}
f
u
∈
{
f
l
+
1
,
f
l
+
2
,
…
,
49
}
CSSSP
T
1
∈
{
0,0.25,0.5,0.75
}
r
∈
{
1,2
,
…
,
10
}
C
∈
{
0,0.01,0.1,0.2,0.5,1
,
2,5
}
FBCSP
T
1
∈
{
0,0.25,0.5,0.75
}
r
∈
{
1,2
,
…
,
10
}
DFBCSP
T
1
∈
{
0,0.25,0.5,0.75
}
r
∈
{
1,2
,
…
,
10
}
F
∈
{
1,2
,
…
,
5
}
CSTFP
r
∈
{
1,2
,
…
,
10
}
F
∈
{
1,2
,
…
,
5
}After we obtained the feature vectors extracted by the filters, spatial weights, and the time windows that are designed by each listed method, we calculated the logarithm of the feature vectors. Then, the LDA for the learning samples was used to obtain a projector onto the 1-dimensional space. The threshold for classification was determined as the middle point of two class averages over the projected learning samples. The feature vectors from the test samples were classified by the projection and the threshold and we obtained the classification accuracy rate in each CV. In Tables4 and 5, CSTFP results in the highest accuracy rate in the average over whole subjects.Table 4
Classification accuracy [%] given by5×5 CV in dataset IVa from BCI competition III. The figure with ± represents the standard deviation (S.D.) over CV.
Method
Subject
Average
aa
al
av
aw
ay
CSP
79.9
±
4.6
98.4
±
1.6
71.9
±
3.5
96.8
±
2.6
92.5
±
3.1
87.9
CSP-Exh
90.8
±
3.6
99.1
±
1.4
74.4
±
5.5
99.1
±
1.4
93.0
±
4.2
91.3
CSSSP
91.5
±
3.7
99.1
±
1.6
71.1
±
6.9
98.8
±
1.4
93.1
±
3.5
90.7
FBCSP
91.3
±
3.4
99.3
±
1.5
70.7
±
7.4
98.2
±
1.6
88.4
±
4.8
89.6
DFBCSP
91.9
±
2.6
98.9
±
1.5
74.8
±
4.0
98.9
±
1.6
96.0
±
3.1
92.1
CSTFP
92.6
±
2.1
98.9
±
1.5
75.4
±
5.3
99.0
±
1.3
96.0
±
2.4
92.5Table 5
Classification accuracy [%] given by5×5 CV in dataset 1 from BCI competition IV. The figure with ± represents the standard deviation (S.D.) over CV.
Method
Subject
Average
a
b
f
g
CSP
89.6
±
5.1
68.8
±
6.7
78.1
±
6.5
93.9
±
3.4
82.6
CSP-Exh
92.3
±
2.5
85.3
±
6.9
78.1
±
6.5
93.9
±
3.6
90.2
CSSSP
89.8
±
5.5
70.4
±
9.5
85.3
±
5.3
93.4
±
4.1
84.7
FBCSP
87.3
±
8.6
63.9
±
14.1
78.6
±
6.9
94.1
±
4.1
81.0
DFBCSP
89.7
±
5.0
85.6
±
6.9
92.5
±
4.4
93.7
±
4.3
90.4
CSTFP
91.4
±
4.0
90.6
±
7.2
93.3
±
3.6
93.5
±
4.1
92.2Moreover, we conducted the classification experiments, in which the parameters shown in Table3 were fixed to show the effects of the changes of the parameters on the classification accuracy rates. Figure 5 shows the accuracy rates of CSP and the rates of each method averaged over all subjects with each T1. For each subject, the parameters shown in Table 3 except for T1 were fixed to the combinations of the parameters that performed the highest accuracy rates with each T1. In Figure 5(a), T1 that perform the highest accuracy rates are different among the subjects. Moreover, Figure 5(b) shows that the accuracy rates highly depend on T1 in the conventional methods. Figure 6 shows the variation of the accuracy rates by the various regularization parameters, ϵ, in CSTFP. For each subject, the parameters, r and F, were fixed to the combinations of the parameters that performed the highest accuracy rates with each ϵ. CSTFP performs higher accuracy as ϵ goes higher than 0.1.The variation of the accuracy rates by variousT1.
(a)
CSP
(b)
Average over all subjectsThe variation of the accuracy rates by various regularization parameters,ϵ.
(a)
BCI competition III dataset IVa
(b)
BCI competition IV dataset 1We show examples of the spatial patterns, the frequency patterns (the amplitude characteristics of the FIR filters), and the time patterns designed by CSTFP in Figures7 and 8. As we can observe in Figures 6, 7(a), and 8(a), the short time windows caused by small ϵ result in poor classification accuracy rates. All of the time windows shown in Figures 7 and 8 remove the samples observed within a few hundreds of milliseconds after the cue. The results suggest that the brain activities related to the task cannot be observed just after the cue. In Figures 7(c)–7(e) and 8(c)–8(e), the time windows do not significantly change in around 0.1–0.4 of ϵ. This result of the relations between ϵ and the time windows corresponds with the result shown in Figure 6 in which the classification accuracy rates do not strongly depend on the regularization parameter, ϵ, by more than 0.12. Moreover, the FIR filters do not highly differ even with various ϵ. However, the spatial weights with the small windows (Figures 7(a) and 8(a)) differ from the other weights with the longer windows. Furthermore, there are slight differences in the time patterns and the frequency patterns between subject aa and subject av.Examples of spatio-time-frequency patterns as the spatial weights,w^1(1) (top), the amplitude response of the filter, h1,…,h4 (center), and the time windows, b1,…,b4 (bottom), designed by the CSTFP method with various ϵ in the result of subject aa. The vertical axis represents normalized amplitude.
(a)
ϵ
=
0.04
(b)
ϵ
=
0.08
(c)
ϵ
=
0.12
(d)
ϵ
=
0.16
(e)
ϵ
=
0.4Examples of spatio-time-frequency patterns as the spatial weights,w^1(1) (top), the amplitude response of the filter, h1,…,h4 (center), and the time windows, b1,…,b4 (bottom), designed by the CSTFP method with various ϵ in the result of subject av. The vertical axis represents normalized amplitude.
(a)
ϵ
=
0.04
(b)
ϵ
=
0.08
(c)
ϵ
=
0.12
(d)
ϵ
=
0.16
(e)
ϵ
=
0.4


## 5. Conclusion

We have proposed a novel method called CSTFP for the classification of EEG signals during motor imagery. Our objective is to design the time windows that are adopted in signal processing for a cue-based BMI. Incorporating the idea into DFBCSP has allowed us to simultaneously design the parameters for the time windows, the spatial weights, and the FIR filters. These parameters are optimized in a single criterion based on the CSP method. We have shown the optimization procedure for the problem of the CSTFP method that is conducted by sequentially and alternatively solving the subproblems into which the original problem is divided. Through the experiments of the artificial signals and actual EEG signals, we have shown the performance of CSTFP. In the experiment, we have demonstrated that CSTFP achieves high classification accuracy rate. Our experimental results also suggest that the CSTFP method can find the frequency bands and the time periods in which brain activities associated with a mental task can be observed.Methods for finding time intervals in which the feature components associated to some brain activities are observed are needed for accurate classification in an asynchronous BCI/BMI [26]. We would like to develop such methods by applying the proposed CSTFP algorithm to our future works.



---

*Source: 537218-2013-11-03.xml*