Comparing Steady-State Visually Evoked Potentials Frequency Estimation Methods in Brain-Computer Interface With the Minimum Number of EEG Channels

Introduction: Brain-Computer Interface (BCI) systems provide a communication pathway between users and systems. BCI systems based on Steady-State Visually Evoked Potentials (SSVEP) are widely used in recent decades. Different feature extraction methods have been introduced in the literature to estimate SSVEP responses to BCI applications. Methods: In this study, the new algorithms, including Canonical Correlation Analysis (CCA), Least Absolute Shrinkage and Selection Operator (LASSO), L1-regularized Multi-way CCA (L1-MCCA), Multi-set CCA (MsetCCA), Common Feature Analysis (CFA), and Multiple Logistic Regression (MLR) are compared using proper statistical methods to determine which one has better performance with the least number of EEG electrodes. Results: It was found that MLR, MsetCCA, and CFA algorithms provided the highest performances and significantly outperformed CCA, LASSO, and L1-MCCA algorithms when using 8 EEG channels. However, when using only 1 or 2 EEG channels d, CFA method provided the highest F-scores. This algorithm not only outperformed MLR and MsetCCA when applied on different electrode montages but also provided the fastest computation time on the test set. Conclusion: Although MLR method has already demonstrated to have higher performance in comparison with other frequency recognition algorithms, this study showed that in a practical SSVEP-based BCI system with 1 or 2 EEG channels and short-time windows, CFA method outperforms other algorithms. Therefore, it is proposed that CFA algorithm is a promising choice for the expansion of practical SSVEP-based BCI systems.


Introduction
rain-Computer Interface (BCI) can provide a new way of communication and control between the human brain and external devices. Although this communication method could be used by everyone, its applications have been mostly used for the patients with disabilities so far (Resalat & Saba, 2016) to improve their quality of life (Tello, Müller, Bastos-Filho, & Ferreira, 2014). Recently, many new BCI applications have been validated for disabled people, including neuroprostheses, spelling for communication, control of artificial hands and legs, and playing with games (Zhang, Xu, Cheng, & Yao, 2014a).
There are different techniques to detect brain activities for the realization of a BCI system, including Electroencephalography (EEG), functional near-infrared spectroscopy, Magnetoencephalography (MEG) and functional Magnetic Resonance Imaging (fMRI). Recording EEG signals from the scalp is a common approach with the neural data signals for BCI devices since it can be easily recorded in most of the circumstances by rather simple equipment. Also, it has a high temporal resolution (Zhang, Zhou, Jin, Wang, & Cichocki, 2014b).
A variety of EEG-based BCI paradigms have been proposed based on different brain responses to transfer the user intent to the computer, such as Sensorimotor Rhythms (SMRs), Slow Cortical Potentials (SCPs), P300 Event-Related Potentials (ERPs), and Steady-State Visually Evoked Potentials (SSVEP) (Zhang et al., 2014b). Recently research on SSVEP-based BCI systems has been extended because of their high Information Transfer Rate (ITR) and easy training phase (Ge et al., 2017).
SSVEP is a neurophysiological reaction excited in the occipital and occipitoparietal region of the brain by a flickering visual stimulus at a specified frequency. Responses of SSVEP contain the basic frequency of the visual stimulus and some of its harmonics are sometimes accompanied with a sub-harmonic. Therefore, SSVEP-based BCI systems understand the user commands when the user looks at one stimulus by identifying the corresponding frequency components in the EEG (Wang et al., 2016). Different approaches have been introduced in the literature to estimate SSVEP responses for BCI applications (Liu, Chen, Ai, & Xie, 2014). A conventional method is the Power Spectral Density Analysis (PSDA) (Wei, Xiao, & Lu, 2011). However, when using this algorithm, the length of processing time window has to be typically more than 3 s to have acceptable frequency resolution, resulting in low ITR (Ming, Xiaorong, Shangkai, & Dingfeng, 2002).
Canonical Correlation Analysis (CCA) (Lin, Zhang, Wu, & Gao, 2007) is a method for exploring the relationships between two multivariate sets of vectors. One of these variables is the EEG signal from various channels and the other is the artificial sine-cosine signal. Researchers have demonstrated that CCA usually outperforms PSDA (Wang et al., 2014). Moreover, high accuracy BCI systems could be developed using CCA in which their required data frame is as short as 2 s (Zhang, Jin, Qing, Wang, & Wang, 2012). Zhang et al. (2011) (Cited byWang et al., 2016 proposed a multi-way development of standard CCA (MCCA) by inspecting the correlation among various variables, including space and trial modes of multidimensional EEG data and sine-cosine signals. MCCA and its L1-regularized development (L1-MCCA) (Zhang et al., 2013) have been introduced to supply refined SS-VEP recognition fulfillment in comparison with standard CCA (Wang et al., 2016).
Least Absolute Shrinkage and Selection Operator (LASSO) algorithm (Zhang et al., 2012) is another method proposed in the literature for computing the contribution of different stimulus frequencies and their harmonics in the recorded EEG signal. The frequency with maximum contribution degree was identified as the goal frequency. LASSO algorithm has been shown to improve frequency identification performance with a shorter time window compared with CCA algorithm (Liu et al., 2014). Zhang et al. (2014b) introduced a Multi-set CCA method (MsetCCA) for frequency identification in SSVEP. In this method, training information is used for reference signals (Zhang et al., 2014b). Accuracy of MsetCCA approach is more than CCA and MCCA methods when using time windows shorter than 2 s (Zhang et al., 2014b). Zhang et al. (2015) considered that a set of EEG signals would share certain common components in response to a specified stimulus on different subjects. These EEG signals contain common components which may have characteristics of SSVEP responses. Therefore, these components could be further efficient reference data for SSVEP identification in using correlation methods. So they proposed a Common Feature Analyzes (CFA) method (Zhang, Zhou, Jin, Wang, & Cichocki, 2015). It outperformed SSVEP identification accuracy compared with those of the CCA and the MCCA methods in a 0.5-s time window (Zhang et al., 2015).
Another approach that has been recently proposed for distinguishing features of SSVEP is Multivariate Linear Regression (MLR) (Wang et al., 2016). This algorithm outperformed CCA and MCCA algorithms and especially provided higher classification accuracies when using a short-time window of 1 s (Wang et al., 2016). The abovementioned methods and their properties are listed in Table 1.
With the advancement in BCI research, these systems are coming out of the lab for practical applications. In addition to high accuracy and information transfer rate, user comfort is critical. It must be met at least in a level that a potential user accepts and continues to work with the system (Qing, Zheng, Yue, Yuankui, & Ge, 2015).
One of the significant factors to ensure user comfort is that the BCI system works with only a few EEG electrodes. A large number of electrodes increases the preparation time beyond the acceptable level for practical use and increases the size and final cost of the system. On the other hand, achieving high and or enough accuracy and information transfer rate with a minimum number of electrodes is challenging due to the reduction of their available information. As shown in Table 1, most of the studies that proposed and compared algorithms for SS-VEP detection used a large number of EEG channels.
The goal of this study was to ascertain which of the newly-proposed algorithms (reviewed above) is more appropriate for the development of a practical SSVEPbased BCI system with high accuracy and information transfer rate when using a minimum number of electrodes. Here CCA (Lin et al., 2007), LASSO (Zhang et al., 2012), L1MCCA (Zhang et al., 2013), MsetCCA (Zhang et al., 2014b), CFA (Zhang et al., 2015) and MLR (Wang et al., 2016) methods are compared using proper statistical methods to determine their performance in real applications.
The rest of the paper is organized as follows. In the next section, information about the experimental protocol and the BCI methods used in this study as well as the statistical methods are presented. Section 3 supplies the results of the performance assessment. Ultimately, the discussion is provided in Section 4.

Methods
Numerous methods have been proposed for SSVEP frequency recognition in the literature. Of those, LASSO (Zhang et al., 2012), L1-MCCA (Zhang et al., 2013), MsetCCA (Zhang et al., 2014b), CFA (Zhang et al., 2015), and MLR (Wang et al., 2016) have demonstrated better performances and therefore we selected them for this study. Also, CCA method is a benchmark method and thus we selected it for this study. These algorithms are described as follows: SSVEP detection based on CCA; SSVEP detection based on LASSO; SSVEP detection based on L1MCCA; SSVEP detection based on MsetCCA; SSVEP detection based on CFA; SSVEP detection based on MLR.
CCA is a multivariate statistical method to explore the underlying correlation between two sets of data. The first set (X) is the EEG recorded from several channels while the second one (Y) is the sine-cosine reference signals reconstructed as below: (1)

Y=
; t= , , ..., , where N refers to the number of harmonics, f i denotes the i th stimulus frequency (fundamental frequency), n represents the number of sampling points, and F s is the sampling rate.
CCA tries to seek a pair of linear transforms, W x and W y that maximizes the correlation between x=X T W x and y=Y T W y . The following optimization problem is solved for each frequency (Zhang et al., 2011). (2) where p i is the association between the recorded signals and the synthetic waveform of the i th frequency. The frequency with the maximum correlation coefficient is selected as the target frequency (Lin et al., 2007). May, June 2019, Volume 10, Number 3 In LASSO algorithm, each EEG trial assumes that SS-VEPs are standard linear regression models (Eq. 1) for the response y ∈ R n .
To create the model of SSVEP identification, a symmetric square-wave signal X, corresponding to the stimulus frequencies, is considered as the reference signal shown in Eq. (3).

X=
;t= , ,..., , where N refers to the number of harmonics, f i denotes the i th stimulus frequency (fundamental frequency), n represents the number of sampling points, and F S is the sampling rate.
The LASSO estimator β ̂ among the EEG signals y and the artificial reference set X were calculated by Eq. (2).
Then the contribution degree of each stimulus frequency is calculated for all the recorded channels of EEG as: , where M equals the number of channels, K denotes the number of harmonics, and CD i refers to the contribution degree of the i th square-wave in the signal. The maximum contribution implies the target frequency which the subject gazing at (Zhang et al., 2012). L1-regularization (L1-MCCA) was proposed in the literature to give a function which can automatically select features for optimizing reference data in SSVEP-BCI detection (Zhang et al., 2017).
To construct the SSVEP recognition model, consider a three-way tensor X ∈ R I×n×K (channel×time×experiment) formed by EEG signals recorded from some channels out of numerous experiments with a particular stimulus frequency and a signal collection Y ∈ R 2N×n shown in Eq. (3).
The optimization problem in L1-MCCA is formulated as: , where w 1 ∈ R I , w 3 ∈ R K , v ∈ R 2N are projection vectors and λ 1 , λ 2 . λ 3 are adjustment parameters. The LASSO estimation is equivalent to (5) (Tibshirani, 2011;Zhang et al., 2012) when any two of w 1 , w 3 , and v are constant. This problem would be resolved by an alternating LASSO method (Zhang et al., 2013).
The Multiset CCA (MsetCCA) method was recently proposed for reference data optimization from a common component in numerous calibration experiments (Nakanishi, Wang, Wang, & Jung, 2015). MsetCCA was extended for correlation maximization among canonical variables from numerous collections of random variates with distinguishing multiple linear transforms (Zhang et al., 2014b). Assume numerous groups of random variables X i ∈ R I ×J i (i=1,2,…,N). In order to maximize the throughout correlation across canonical variables, the MAXVAR objective function is characterized as: , where C ij =X i X j T is the between-set covariance matrix, and ρ is the correlation coefficient. The objective function in (6), may be transformed into the eigenvalue problem with the approach of Lagrange multipliers as: Assume, X 1,m , X 2,m , ..., X N,m ∈ R C×P (C channels×P points) demonstrate EEG signals groups, including experiments at the m-th stimulus frequency f m . The MsetCCA is performed to distinguish numerous spatial filters w 1,m , w 2,m , ..., w N,m , to maximize throughout correlation between the canonical variables z 1,m , z 2,m , ...,z Nm ̃ ̃ ̃ with the linked spatial filtering z i,m =w T i,m X i,m ̃ (i=1,2,…,N). These canonical variables show the common components between numerous training signals considered to have better accuracy of the real SSVEP specifications in comparison with sine-cosine reference signals. Canonical variables were combined to construct the reference signal optimization at frequency f m . These optimized reference signals are defined as: For each stimulus frequency f m , the corresponding reference signal Y m is consumed for computing the maximum correlation coefficients with the EEG signals (Zhang et al., 2014b).
Another approach is to explicitly model the common and distinct components. This approach has a similar concept with the multiset approach (Wang et al., 2016) to exploit the common features which multiple EEG signals share at the same stimulus frequency.
In this method, a set of matrices X={X k ∈X I ×T :k=1,2,…,K} k , share at least one common dimension T. They can be, for example, a set of multichannel EEG signals (channel×time×point) recorded for the same visual stimulus but from different subjects. These data matrices can be factorized in a linked way shown in Eq. (9): (11) is the number of latent components with R k <I k and C is the number of common components. Finally, Ā k and Ă k are the partitions of the mixing coefficients A k corresponding to B and B k ̌. For M stimulus frequencies, the common features B k (k=1,2,..., K) at each frequency would be exploited for a new test signal x∈X R ̌T ; the target frequency is identified as (Zhang et al., 2015): (12) LR is a technique for modeling the relationship between a scalar dependent vector and one or more independent vectors. Consider EEG training data X=[x 1 ,x 2 ,…,x N ]∈R D×N , where D demonstrates the dimension of feature (D=C channels×P temporal points) and N is the number of signal points. Corresponding to the training points, a matrix of the label is created, using one of M class coding. The goal of MLR is to find distinguished subspaces by minimizing the objective function as follows: , where W=[w 1 ,w 2 ,…,w C ]∈R S×C demonstrates the projection matrix and b is the model separation.
The model of MLR is then considered as:

Data Description
EEG data from an online available dataset (Wang et al., 2016) were used in our study. The dataset contained data from ten subjects (all males, aged from 21 to 27 years) and recorded from 8 channels (P7, P3, Pz, P4, P8, O1, Oz and O2). The sampling rate was fixed at 250 Hz. During the experiment, the participants were seated in a comfortable armchair 60 cm away from the center of the monitor. The experiment was performed in a shielded room. The EEG data were then bandpass filtered from 4 to 45 Hz. Four frequencies (6, 8, 9, and 10 Hz) were adopted in a recording session. The subjects were asked to gaze at each stimulus frequency for 4 s. Overall, there are 80 trials in this dataset. This dataset has been widely used in the literature as the benchmark (Wang et al., 2016;Zhang et al., 2012;Zhang et al., 2014b;Zhang et al., 2013;Zhang et al., 2015).

Experimental evaluation
In this study, the new algorithms, including MLR and CFA are compared with MsetCCA, L1MCCA, LASSO, and CCA algorithms to assess their efficacy for target frequency detection in a practical BCI system. The number of harmonics needed to define reference signals was set to 2 for CCA, LASSO, and L1-MCCA algorithms.
In LASSO and L1-MCCA methods, the lambda parameters were set to 0.5 and 0.02, respectively (Zhang et al., 2012;Zhang et al., 2013). The leave-one-run-out cross-validation was used to assess the average detection accuracy in the entire analyzed methods. The signals from 19 runs were used as training signals whereas the signals from the left-out runs were used for validation (Zhang et al., 2014b).

Validation
The performance of the methods was assessed using the F-score (F1 score). It was shown in the literature that this criterion is more suitable than the overall accuracy of multi-class problems (Sokolova & Lapalme, 2009) because the latter index overestimates the performance of the analyzed methods. Moreover, by using F-score, it is possible to analyze the performance of the system for each class better than the traditional accuracy measurements. For the calculation of F-score, first, the multiclass confusion matrix of four frequencies is created. Then, precision, sensitivity, and F-score are calculated for each frequency class as below:

Pr= TP TP+FP
, where Pr is the precision TP and FP and are the number of true and false positive predictions for each class. (16)

Se= TP TP+FN
, where Se is the sensitivity and FN is the number of false negative predictions for each class. (17)

F-score= 2×Pr×Se
Pr+Se , where F-score is the harmonic mean of the precision and sensitivity (=recall) of the corresponding class. This parameter is not dependent on TN (True Negative), and thus the performance of the analyzed system is not overestimated (Marateb, Mansourian, Adibi, & Farina, 2014).

Statistical methods
It is important to use proper statistical tests for rigorous comparison between different approaches, otherwise random or insignificant differences are considered as significant (Bossuyt et al., 2015;Karimimehr et al., 2017). The McNemar's test (Webb & Copsey, 2011) was used to identify SSVEP frequency estimation algorithms outperforming the other methods. Pairwise comparisons were performed using the McNemar's test. Then, the methods with higher and significant performance were shown as significantly outperforming methods.
For further analysis, GEE (Generalized Estimating Equation) (Hardin & Hilbe, 2007) was applied to determine any significant differences among the selected algorithms when different electrode montages (monopolar or bipolar) and the number of channels (1 or 2 channels) were used with repeated measurements in the time windows of 1 second. The level of statistical significance was considered 0.1 to be on the safe side and not to miss significant differences in this small sample size data (Lavrakas, 2008).

Results
We compared the performance of different SSVEP frequency estimation algorithms. Figure 1 demonstrates the resulting F-scores for time windows of 0.5 to 4 seconds when all 8 channels of recorded EEG were used. MLR method outperformed the other 6 methods in terms of mean F-score performance for the time windows of 0.5 to 1.5 s. MsetCCA and CFA methods were ranked second, while the performance of the CFA method was significantly higher than that of MsetCCA method for the time window of 0.5 s (P<0.05). Overall, LASSO had the lowest performance compared with the others.
Overall, pairwise McNemar's test of these results demonstrated that MLR, MsetCCA, and CFA algorithms significantly outperformed the other analyzed methods in the detection of gazed frequencies using 8 EEG channels (P<0.05). In fact, in the entire analysis windows, MLR, MsetCCA, CFA, L1-MCCA, CCA, and LASSO significantly outperformed 5, 3, 3, 2, 1, and 0 times in pairwise comparisons. Thus, only those three methods were studied in the rest of the paper. In a practical SSVEP-based BCI system, when only one channel of EEG is used, the performance of the algorithms might change compared with when using more channels. Therefore, the F-score measures of MLR, MsetCCA, and CFA algorithms were compared when only one channel of EEG (Oz-Pz channel) is used (Diez, Mut, Laciar, & Avila, 2010) (Figure 2).
For one-channel EEG, McNemar's test was implemented to determine significant differences between CFA, MLR, and MsetCCA methods in different time windows (Table 2). In time windows smaller than 2 s, the CFA method had better F-score in target frequency recognition than MLR and MsetCCA methods. To do further extensive analysis, the F-score values for MLR, CFA, and MsetCCA methods were calculated for a time window of 1 s, when different electrode montages (monopolar or bipolar) and the number of channels (1 or 2 channels) were used ( Table 2).
The results are based on leave-one-out cross-validation test on all 20 trials of the data. CFA and MLR methods outperformed MsetCCA in the entire scenarios. Moreover, on average, CFA outperformed MLR method. Generalized Estimating Equation (GEE) test was applied to determine any significant differences between these algorithms when different scenarios (Table 3) were used. The result showed that CFA method significantly outperformed the other two algorithms (P<0.07). Also, MLR significantly outperformed MsetCCA method (P=0.001).
The other important factor in the development of an online BCI system is the computational cost of the algorithms. The average running time for 19 runs of algorithms was obtained on a laptop with Windows 8.1 Operating System and 2.6 GHz Intel Core i5 CPU with 2 GB RAM and all of the algorithms were implemented in MATLAB 2014a (Table 4). It should also be stated that the MLR and CFA methods were implemented in the vectorized form.  In the training phase, MLR is the fastest algorithm and L1-MCCA is the most time-consuming algorithm. However, in applying the tuned algorithm to frequency detection in the test set, CFA provided the fastest computations, while LASSO took the longest computational time.

Discussion
In a practical SSVEP-based BCI system, it is important to use accurate frequency detection algorithms when a short-time window with only a few EEG channels are used for analysis. In this study, the state-of-the-art algorithms for SSVEP detection, including CCA, LASSO, L1-MCCA, MsetCCA, CFA, and MLR were compared on a benchmark database. The results demonstrated higher F-scores for MsetCCA, CFA and MLR methods in comparison with CCA, LASOO, and L1-MCCA methods for entire time windows when 8 channels of EEG are used. These results are in accordance with the results reported by Wang et al. (2016), demonstrating a higher performance for MLR method in comparison with CFA, MCCA, and CCA methods. Zhang et al. (2014b) suggested that MsetCCA provided a higher accuracy in comparison with MCCA and CCA methods.
Here it is demonstrated that MsetCCA also outperforms L1-MCCA. In another study, Zhang et al. (2015) showed that the CFA method provided higher accuracies than MCCA and CCA methods. Here it is demonstrated that CFA method also provides higher F-scores in comparison with L1-MCCA.
Further analysis demonstrated that while MLR method outperforms the other two algorithms when several EEG channels are used for analysis, it is the CFA method that provides the highest F-scores when only 1 or 2 EEG channels are used. We also compared the performance of CCA with a set of monopolar signals such as Oz and Pz or bipolar (Oz-Pz) signal with that of the other methods in Figure 2 (namely as CFA, MLR, MsetCCA). The results did not change and CFA and MLR methods still outperformed the other two algorithms for time windows shorter than 1 s. Table 3 can also be examined by another aspect. The 8th row in Table 3 represents the F-scores for Oz-Pz EEG channel, while in rows 1 to 7, the F-scores are provided when an additional bipolar channel was used. In most cases, the MLR method could effectively use additional information provided by the extra channels to improve the F-score of the system. In the CFA method, although the F-scores are higher, using an extra EEG channel does not have a consistently positive effect on the F-scores and the same can be observed for MsetCCA method. This capability of the MLR technique can also be observed when comparing the F-scores for Oz montage with Oz, O1 and Oz, O 2 montages. Therefore, the MLR method could use information in multichannel EEG towards providing higher F-scores.

The results presented in
In practical applications of an SSVEP-based BCI system without the help of experts, the proper positioning of the electrodes is also a challenging issue. Therefore, another question is how much the algorithms should be The best result for the testing phase is displayed in bold.
robust to correct the position of electrodes. The results provided in Table 3 demonstrated that the performance of frequency detection algorithms is quite sensitive to the placement of Oz electrode. In fact, changing the position of this electrode to O1 or O2 considerably degrades the F-score of the system. This result can be observed when Oz or Oz-Pz montages were used. However, the results show that the CFA method was less sensitive to the placement of Oz electrode compared with other methods. The comparison of F-score values between monopolar and bipolar montages in Table 3 indicates that Oz-Pz montage has better F-score value than other monopolar montages. This result is in accordance with the result reported by Diez et al. (2010).
The CFA algorithm also has the lowest computational time which is important for online implementation of a BCI system. Therefore, it is proposed that the CFA algorithm may be a proper choice in the development of practical SSVEP-based BCI systems. In the literature, different methods were only compared with CCA and MCCA methods (Zhang et al., 2012;Zhang et al., 2014b;Zhang et al., 2013;Zhang et al., 2015). Our study suggests that newly-developed algorithms could also be compared with CFA due to its higher performance for short-time windows when few EEG channels are used.
The limitations of the study are the small dataset, a small number of participants, and the limited number of frequencies. Therefore it is suggested to test these methods on a larger dataset. Moreover, more detailed analysis of sensitivity on the electrode placement can be performed when a large number of electrodes are used.

Compliance with ethical guidelines
There was no ethical considerations to be considered in this research.

Funding
This research did not receive any specific grant from funding agencies in the public, commercial, or not-forprofit sectors.