A Least-Square Unified Framework for Spatial Filtering in SSVEP-Based BCIs

The steady-state visual evoked potential (SSVEP) has become one of the most prominent BCI paradigms with high information transfer rate, and has been widely applied in rehabilitation and assistive applications. This paper proposes a least-square (LS) unified framework to summarize the correlation analysis (CA)-based SSVEP spatial filtering methods from a machine learning perspective. Within this framework, the commonalities and differences between various spatial filtering methods appear apparent, the interpretation of computational factors becomes intuitive, and spatial filters can be determined by solving a generalized optimization problem with non-linear and regularization items. Moreover, the proposed LS framework provides the foundation of utilizing the knowledge behind these spatial filtering methods in further classification/regression model designs. Through a comparative analysis of existing representative spatial filtering methods, recommendations are made for the superior and robust design strategies. These recommended strategies are further integrated to fill the research gaps and demonstrate the ability of the proposed LS framework to promote algorithmic improvements, resulting in five new spatial filtering methods. This study could offer significant insights in understanding the relationships between various design strategies in the spatial filtering methods from the machine learning perspective, and would also contribute to the development of the SSVEP recognition methods with high performance.


I. INTRODUCTION
I N RECENT years, numerous brain-computer interface (BCI) paradigms have been developed to facilitate direct communications between human intentions and the external environment [1].Among these paradigms, the steady-state visual evoked potential (SSVEP)-based BCI has become a prominent BCI paradigm owing to its high performance, minimal training requirements, and small individual performance differences [2], [3], [4].The SSVEP-based BCI has been widely applied in rehabilitation and assistive applications [5], [6].The SSVEP-based BCI utilizes visual flickers with varying frequencies and phases to encode commands, while the command selections are identified by recognizing the SSVEP components from electroencephalogram (EEG) signals.To enhance the recognition performance, the spatial filtering technique is commonly employed [7], [8].It involves a weighted linear combination of EEG signals recorded from multiple channels, where the weights assigned to channels are collectively referred to as the spatial filter.
Currently, various correlation analysis (CA)-based methods have been proposed and have achieved state-of-the-art performance in SSVEP recognition [2], [8], [9], [10], [11], [12], [13], [14].Nevertheless, the relationships between these CA-based methods and the classification/regression models have been rarely explored.On the one hand, most existing reviews primarily focus on performance comparisons, with limited exploration of their computations, and thus cannot summarize such relationships [7], [15], [16].On the other hand, our previous study stands as the only systematic and methodological review, summarizing sixteen spatial filtering approaches within a generalized eigenvalue problem (GEP) framework and categorizing them into two types [8], while the GEP framework also fails to elucidate spatial filtering methods from the machine learning perspective.This limitation impedes the further utilization of the knowledge behind existing spatial filtering approaches in further SSVEP classification model designs, such as the projection matrix designs for extracting SSVEP features, the combination methods of inter-and intra-class features, and the transfer strategies of spatial filters.Additionally, the eigen-formulation of existing spatial filtering methods introduces critical analytical and computational drawbacks for the future developments, e.g., the difficulties of incorporating the regularization, integrating non-linear transforms, and extending to the knowledge-embedded artificial neural network (ANN) designs.Moreover, with the rapid development of SSVEP recognition algorithms, certain advanced spatial filtering methods, such as the task-discriminant component analysis (TDCA) that exhibits superior individual SSVEP recognition performance [13], cannot be integrated into the GEP framework.Therefore, an comprehensive review from the machine learning perspective is imperative to addresses these limitations.
This study applies the least square (LS) regression model, a commonly used machine learning model, to summarize various CA-based spatial filtering methods as a unified framework.In existing literature, component analysis methods and the common spatial pattern (CSP) method were reformulated as the LS problems, which effectively converts eigendecomposition problems into regression problems [17], [18].Although the frameworks of the component analysis and the regularization of the CSP method are not related to the SSVEP spatial filtering methods, these studies highlight that the LS format provides greater flexibility compared to the eigendecomposition format, making it more amenable to the incorporation of constraints, non-linearity, and ANNs, and thus is suitable for establishing a unified framework of the SSVEP spatial filtering method in this study.
The key contributions of this study are summarized as below: 1) This study comprehensively consolidates 11 existing CA-based SSVEP spatial filtering methods into a LS unified framework, leading to a better understanding of the spatial filtering methods from the machine learning perspective; 2) With the LS framework as a guide, this study summarizes the strategies of computing spatial filters into 5 distinct categories and performs a comparative analysis, facilitating a deeper exploration of relationships among existing spatial filtering approaches; 3) Beyond the spatial filtering computations, this study further summarizes the strategies of utilizing spatial filters into 4 distinct categories and compares these strategies, providing a comprehensive analysis of factors affecting recognition performance; 4) Based on these comparison results, the strategies with the superior and robust performance are recommended and integrated to fill the research gap and demonstrate the ability of the proposed LS framework to promote algorithmic improvements, leading to the development of 5 new spatial filtering methods.Overall, this study could provide significant contributes to the understanding of their computation strategies from the machine learning perspective, and holds the potential for future developments of SSVEP recognition methods.Codes are available at https://github.com/pikipity/SSVEP-Analysis-Toolbox

A. Notations
For a clear introduction of the proposed LS framework, Table I summarizes the key notations used in this paper.The variables with right subscripts n and m denote that they are corresponding to the n-th trial and the m-th stimulus.More details related to the SSVEP reference and template signals, the QR decomposition, and the singular value decomposition (SVD) can be found in Section S.I of the supplementary material.

B. GEP Framework of CA-Based Spatial Filtering Methods
Our previous study summarized 6 types of CA-based spatial filtering methods and then proposed 3 new CA-based spatial filtering methods [8].These 9 methods, i.e., standard canonical correlation analysis (sCCA), individual template based CCA (itCCA), transfer template CCA (ttCCA), extended CCA Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
(eCCA), multi-set CCA (MsetCCA), multi-stimulus CCA (ms-CCA), MsetCCA with reference signals (MsetCCA-R), task related component analysis (TRCA), and TRCA with reference signals (TRCA-R), were unified as a GEP that is where Z denotes the processed EEG signals, W denotes the eigenvectors and works as the spatial filter, denotes the corresponding eigenvalues, B is the covariance matrix of Z, and D denotes the temporal filter that is pre-defined and computed by the SSVEP reference or template signals.

C. Statistical Analysis Method
In the performance comparisons of following sections and the supplementary material, the paired t-test is utilized.The p-values are corrected by the Bonferroni correction [19].The significance of an observed difference in the N paired t-tests is denoted as "*" if the p-value is less than 0.05/N , "**" if the p-value is less than 0.01/N , "***" if the p-value is less than 0.001/N , and "n.s." if the p-value is larger than 0.05/N .

A. Reduced Rank Regression-Based Unified Framework
The reduced rank regression (RRR) is an extension of the LS regression, and was firstly introduced in early 1950s by Anderson [20].The RRR model has been utilized in various fields, such as the pattern recognition [21], [22], the time series analysis [23], [24], and the financial analysis [25].Compared to the conventional LS regression, the RRR aims to minimize the LS error subject to a rank constraint on the regression matrix, which can effectively reduces the number of free trainable parameters.With the rank constraint, the regression matrix in the LS regression can be parameterized as the outer product of two matrices of the rank constraint.This study proposes formulating various CA-based SSVEP spatial filtering methods as one kind of RRR problems: where W is the spatial filter.In principle, the matrix E presents the combined EEG features of inter-classes, i.e., where L E denotes the combination matrix of inter-class EEG features, P E denotes the orthogonal matrix applied for generating inter-class EEG features, and Z denotes the EEG data.
In addition, the matrix T is Similarly as E, K presents the combined EEG features of intraclasses, i.e., where L K denotes the combination matrix of intra-class EEG features, and P K denotes the orthogonal matrix applied for generating intra-class EEG features.Following [17] and [18], The RRR problem shown in (2) can be solved by the alternated least squares, which is illustrated in detail in Section S.III of the supplementary material.

B. Standard CCA
The sCCA is the most conventional CA-based spatial filtering method for the SSVEP recognition [10].For the m-th stimulus, the sCCA finds the corresponding spatial filter to maximize the correlation between the processed EEG signal X and the SSVEP reference signal Y m : According to [8], W m can be obtained by equivalently solving Utilizing the SVD, the solution of ( 7) can be found by solving the following eigenvalue problem: where T sCCA = P(Y m ) X S(X).According to [26], finding the eigenvectors in ( 8) is equivalent to minimizing Let M sCCA = XW m , the spatial filters of the sCCA can be obtained by minimizing C. Individual Template Based CCA The computation of spatial filters in the itCCA is similar as the sCCA [27].To avoid the interference from the spontaneous EEG activities, the itCCA applies the SSVEP template signals to replace the SSVEP reference signals used in the sCCA.Therefore, the spatial filters of the itCCA also can be calculated by minimizing where M itCCA = M sCCA , and T itCCA = P X m X S(X).

D. Transfer Template CCA
The core idea of the ttCCA is similar as the sCCA [28].The spatial filters of the ttCCA is designed to maximize the correlation between the SSVEP templates and the SSVEP reference signals.This idea was firstly proposed for the cross-subject SSVEP recognition by using signals collected from non-target subjects to compute the SSVEP template signals.The similar idea was also applied for the individual SSVEP recognition by using the signals from the target subject to compute the SSVEP template signals.Because the computation steps of these cross-subject and individual SSVEP recognition methods are the same, this study considers them together as the ttCCA.Following the derivation shown in Section III-B, the computation of the spatial filters in the itCCA can be obtained by minimizing where M ttCCA = X m W m , and

E. Extended CCA
The eCCA ensembles the sCCA, the itCCA and the ttCCA together, and computes 3 types of spatial filters [2].The computation processes of spatial filters are exactly the same as the sCCA, the itCCA, and the ttCCA, and thus will not be repeated.The spatial filter computation methods in the eCCA following the sCCA, the itCCA, and the ttCCA are termed as eCCA-1, eCCA-2, and eCCA-3, respectively.

F. Multi-Stimulus CCA
The ms-CCA extends the computation of the spatial filters in the ttCCA from the single stimulus case to the multi-stimulus case [12].It calculates the spatial filters by maximizing the correlations of the SSVEP templates and the SSVEP references for multiple stimuli simultaneously.In practice, the computation of the spatial filters in the ms-CCA can be expressed as arg max Following the derivation shown in Section III-B, the spatial filters of the ms-CCA can be computed by minimizing where

G. Task-Related Component Analysis Without and With Reference Signals
The TRCA does not use the averaged calibration signals as SSVEP templates and sine/cosine signals as reference signals, but directly maximizes the correlations between calibration signals in different trials [11].It can be achieved by maximizing the inter-trial covariance of calibration signals, i.e., arg max where inter m denotes the inter-trial covariance of EEG signals for m-th stimulus, i.e., inter and intra m denotes the intra-trial covariance of EEG signals for m-th stimulus, i.e., intra . Equivalently, the spatial filters of the TRCA can be obtained by solving the following eigenvalue problem: Following [26], finding the eigenvectors of ( 16) is equivalent to minimizing Then, the spatial filters of the TRCA can be obtained by minimizing where After the TRCA, Wong et al. introduced the orthogonal projection matrix to the TRCA, termed as TRCA-R [8].The TRCA-R maximizes the inter-trial covariance of EEG features to find the spatial filters, i.e., where F n,m , and Following the derivation from ( 16) to (18), the eigenvectors of ( 19) can be obtained by minimizing where M TRCA-R = M TRCA , and X n,m .

H. Multi-Stimulus Task-Related Component Analysis
The multi-stimulus TRCA (ms-TRCA) integrates the computation strategy used in the ms-CCA into the TRCA, which simultaneously considers calibration trials of multiple stimuli [12].The spatial filters of the ms-TRCA can be calculated by solving the following eigenvalue problem: Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
Following the derivation from ( 16) to (18), the eigenvectors of ( 21) can be calculated by minimizing where X n,m W, and .
The computations of the spatial filters in the ensemble TRCA (eTRCA), the eTRCA with reference signals (eTRCA-R), and the multi-stimulus eTRCA (ms-eTRCA) are the same as those in the TRCA, the TRCA-R, and the ms-TRCA, respectively.The major difference is how to use the spatial filters for the SSVEP recognition, which will be summarized in the Section IV-B.In this study, to simplify the expressions of these method names, (e)TRCA denotes TRCA and eTRCA, (e)TRCA-R denotes TRCA-R and eTRCA-R, and ms-(e)TRCA denotes ms-TRCA and ms-eTRCA.

I. Task-Discriminant Component Analysis
The TDCA applies three strategies of computing the spatial filters to improve the recognition performance [13].Firstly, EEG signals are augmented across channels by including delayed signals.Secondly, the EEG signals are extended across samples by combining the original EEG signals and the corresponding EEG features.For the m-th stimulus, the extended EEG signal in the n-th trial can be expressed as Thirdly, the TDCA applies the two-dimensional linear discriminant analysis to find the common spatial filter of multiple stimuli, i.e., where Compared with the eigenformulation of the ms-TRCA, the TDCA further removes the averaged trends of calibration signals.Therefore, the TDCA can be regarded as an extension of the ms-TRCA.Following the derivation from ( 21) to (22), the eigenvectors of ( 24) can be computed by minimizing where W, and .
It is worth noting that, although (24) is also an eigendecomposition problem, it cannot be covered by the GEP framework proposed in [8] due to the projection matrix in the right side and the differences of the data combination way between the left and right sides.

J. Multi-Set CCA Without and With Reference Signals
The optimization target of computing the spatial filters in the MsetCCA is to maximize the inter-trial covariance of calibration signals, which is the same as the TRCA [29].However, the MsetCCA assumes that the spatial filters are different across trials, which is different from the spatial filtering methods introduced in above sections.The spatial filters in the MsetCCA can be obtained by solving the following eigenvalue problem: Compared with ( 16), the only difference is the combination way of calibration signals in different trials.Therefore, following the derivation from ( 16) to (18), the eigenvectors of ( 26) can be computed by minimizing where X n,m W m , and Similarly as the TRCA-R, Wong et al. also incorporated the orthogonal matrix into the MsetCCA to maximize the inter-trial covariance of EEG features with different spatial filters across trials, termed as MsetCCA-R [8].Following the derivation from (19) to (20), the spatial filters of the MsetCCA-R can be obtained by where M MsetCCA-R = M MsetCCA , and IV. CATEGORIZATIONS, COMPARISONS, AND RECOMMENDATIONS OF SPATIAL FILTERING COMPUTATION AND USAGE STRATEGIES With the help of the proposed LS unified framework, the commonalities and differences among the spatial filtering computations in these CA-based methods are summarized by grouping them into 5 categories.Additionally, besides the computation of spatial filters, the strategies related to their utilization are further summarized into 4 types.Finally, the prominent strategies are recommended through simulations.

A. Categorizations of Spatial Filtering Computation Strategies
In the LS framework, the matrices L E and P E are related to the inter-class features.The matrices L K and P K are related Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.

TABLE II PARAMETERS IN PROPOSED LS UNIFIED FRAMEWORK OF CA-BASED SPATIAL FILTERING METHODS
to the intra-class features.L E and L K present how to combine EEG features of inter-and intra-classes, respectively.P E and P K are the orthogonal projection matrices that extract EEG features of inter-and intra-classes, respectively.The combined intra-class features K can be regarded as the normalization item of the combined inter-class features E in the output of the LS regression T. With the help of the proposed LS framework shown in (2), the computation strategies of spatial filters in the CA-based methods is revealed, especially for the data combination strategies and the feature extraction strategies between inter-and intra-classes.The parameters in the propsoed LS framework of the CA-based methods are summarized in Table II.
Except L E that has same functions for all methods, other parameters in the LS framework have different forms and functions.As shown in Fig. 1, the spatial filtering computation strategies in the existing CA-based methods can be grouped by 5 categories based on the differences of parameters in the LS framework: 1) How to extract inter-class features, which is related to P E ; 2) How to concatenate data in Z; 3) Which data is applied for the spatial filtering computation, which is related to the data types included in Z; 4) How to combine intra-class features, which is related to L K ; 5) How to extract intra-class features, which is related to P K ;

B. Categorizations of Spatial Filtering Utilization Strategies
Except the computations of spatial filters shown in Section III, how to use the spatial filters also has large effects on the SSVEP recognition.All CA-based methods follow the same computation paradigm.Firstly, they measure the correlation coefficients of the processed EEG signals and the calibration signals or the SSVEP reference signals for all stimuli.The correlation coefficients of the m-th stimulus can be computed by where Z denotes the reorganized processed EEG signal, W m is the spatial filter of the m-th stimulus in the recognition process, Y m is the calibration signals or the SSVEP reference signals of the m-th stimulus in the recognition process, V m contains the harmonic weights of the m-th stimulus in the recognition process, and D m is applied for the normalization, i.e., Then, the CA-based methods combine all elements in r m for each stimulus.Finally, the stimulus with the highest combined Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.correlation coefficient is regarded as the target.The major differences of the recognition between different CA-based methods are the strategies of constructing W m , Z, V m , Y m and the combined correlation coefficients.These strategies can be summarized as four types.
In the first type of strategies, most CA-based methods, including the sCCA, itCCA, ttCCA, TRCA, ms-CCA, ms-TRCA and TDCA, recognize the processed EEG signal by directly applying the obtained spatial filters and the SSVEP template/reference signals.Because these methods only contain one element in r m , they directly use r m to identify the recognition target.
As mentioned in Section III-E, the eCCA integrates spatial filters of the sCCA, itCCA, and the ttCCA together to improve the recognition performance, which is regarded as the second type of strategies.The eCCA utilizes 3 , and where W m,1 , W m,2 and W m,3 denote the spatial filters W m obtained from the eCCA-1, the eCCA-2, and the eCCA-3, respectively, and U m denotes the harmonic weights obtained from the eCCA-1.Four correlation coefficients obtained from different spatial filters are combined by sign {r m } • r T m × r m where denotes the summation of all elements, "sign" denotes the element-wise signum function.This ensemble strategy can also be used as a general method of combining the recognition results of multiple CA-based methods, and has shown prominent improvements of the recognition performance in many approaches, such as the combination of the ms-CCA and the ms-eTRCA [12] and the stimulus-stimulus transfer method [30].
The eTRCA, the eTRCA-R and the ms-eTRCA introduces the third type of strategies.They combine spatial filters of multiple stimuli together for the recognition, i.e., In addition, r m also only contains one element and can be directly compared for the recognition.It should be noticed that the ms-eTRCA computes spatial filers by only considering the calibration signals of the stimuli with close frequencies, resulting in different spatial filers for different stimuli [12].
The MsetCCA and the MsetCCA-R assume that the spatial filters of different trials are different, and combine these spatial filters of various trials together for the recognition, which is regarded as the fourth type of strategies.The spatial filters W m obtained from the computation shown in Section III-J contain the spatial filters of all trials, i.e., W m = Then, the MsetCCA and the MsetCCA-R re-calculate the spatial filters and recognize the target following the itCCA with the utilization of the new template signal Y m .

C. Strategy Comparisons
To compare various strategies of computing and utilizing spatial filters in existing CA-based spatial filtering methods, Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.

TABLE III METHODS WITH HIGHEST AVERAGED ACCURACY AT 0.25S SIGNAL LENGTH FOR FIVE DATASETS TABLE IV METHODS WITH HIGHEST AVERAGED ITR IN FIVE PUBLIC SSVEP DATASETS
simulations are carried out on 5 public with various subject numbers, stimulus numbers, stimulus designs, and devices, i.e., Benchmark Dataset BETA Dataset [32], eldBETA Dataset [33], UCSD Dataset [15], and Wearable SSVEP BCI Dataset [34].More details related to these datasets can be found in Section S.II of the supplementary material.
The performance comparisons are based on the individual classification with the leave-one-block-out validation.More spcifically, for each stimulus, because 1 block only contains 1 trial, 1 trial is used for testing, and remaining trials are adopted to computing spatial filters and SSVEP templates.This entire evaluation is repeated to test all trials.The recognition performance is evaluated by the recognition accuracy and the information transfer rate (ITR).The ITR is computed by , where P is the recognition accuracy, and T is the summation of the break time designed for shifting visual attention, the visual latency time, the signal length, and the computation time of the recognition process [2].The recognition performance is evaluated on 4 signal lengths: 0.25s, 0.50s, 0.75s, and 1.00s.The detailed results are illustrated from Section S.IV to Section S.VIII in the supplementary material.The methods with highest averaged accuracy at short signal length (0.25s) and with highest ITR are summarized in Table III and Table IV.

D. Strategy Recommendations
As shown in Fig. 1, the entire evolution of CA-based spatial filtering methods encompasses the transition from utilizing single-trial data with a single stimulus to employing multi-trial data with multiple stimuli.Since the spatial distribution variances of different trials in single stimulus are small, the early attempts, such as the (e)TRCA and the eCCA, integrate calibration signals of multiple trials in one stimulus to find the spatial filters.Due to the costly EEG collection procedures, the availability of calibration data pre-recorded for a single stimulus is severely restricted, which has become a well-known practical challenge in SSVEP-based BCIs, particularly when dealing with a large number of stimuli [30], [35], [36], [37], [38].One natural idea of addressing this issue is transferring calibration data from a non-target domain to the target domain.As the concept of leveraging the commonality of spatial filters across different stimuli is introduced in [12], an increasing number of approaches utilize calibration signals from multiple stimuli to find the common spatial filters, and achieve the data transfer cross classes.According to the comparison results shown in Table III and Table IV, this strategy demonstrates notable performance in most cases.However, as pointed in [12], only signals of stimuli with close frequencies have the similar spatial distributions.To achieve the stimulus-stimulus transfer methods with broader applicability, the signal aliments are required [30], [37].In addition to the stimulus-stimulus transfer, several cross-subject, cross-device and cross-dataset transfer approaches have been proposed recently and have achieved prominent performance [28], [35], [36], [39].Although achieving the general transfer learning across multiple domains remains an open issue for further investigation, incorporating calibration data from multiple domains is a valuable strategy of obtaining the reliable spatial filters.
Most CA-based spatial filtering methods recommend utilizing the SSVEP features extracted by the SSVEP reference/template signals.In light of the frequency-and phase-locking characteristics of the SSVEP signals, the utilization of SSVEP reference/template signals is prevalent in these methods.Nevertheless, (e)TRCA and the ms-(e)TRCA deviate from this approach and directly employing the sampling points in EEG signals.They outperform most methods that rely on SSVEP features as shown from Section S.IV to Section S.VIII in the supplementary material.This superior performance can likely be attributed to the fact that certain characteristics of SSVEP signals elude the extraction through the use of SSVEP reference/template signals, such as the delayed responses of the second-order linear system reported in [40].By considering both the sampling points as the temporal features and the features extracted by SSVEP reference signals as the spectral features, the TDCA integrates features from different domains together.As shown in Table III and Table IV, this approach yields significant and reliable performance.While further studies are needed to determine effective SSVEP feature extraction methods, the projection matrix design strategy employed in TDCA proves benefits in enhancing recognition performance.
For the strategy of utilizing spatial filters, directly applying them can yield good ITR performance due to the prominent performance of the TDCA.However, combining multi-stimulus spatial filters can offer enhanced stability, especially for shorter signal lengths, as demonstrated in Table III.By considering spatial filters as different recognition models, combining multi-stimulus spatial filters can be treated as an ensemble technique that is fundamental and commonly employed in the machine learning [41], [42].Although the optimal methods of integrating multiple spatial filters require further investigation, aggregating spatial filters from different trials and stimuli is recommended for achieving robust performance.

V. NEW SPATIAL FILTERING METHODS
To demonstrate that the proposed LS framework could facilitate the development of new methods and address research gaps, 5 new spatial filtering methods are designed by integrating the strategies recommended in Section IV-D.The parameters of the proposed LS framework employed in these newly developed methods have been summarized in Table II.
A. New TRCA-Related Methods 2 extensions of the ms-TRCA can be obtained by introducing the projection matrices of the SSVEP reference signals into P E and P K .One method, termed as ms-TRCA-R-2, only includes the projection matrix in P E , which also can be regarded as the multi-stimulus extension of the TRCA-R.The other method, termed as ms-TRCA-R-1, includes the projection matrix in both P E and P K .Specifically, P E and P K of these two new methods can be expressed as Because the TDCA can provide the best performance in most datasets as listed in Table III and Table IV, L E , L K and Z of the ms-TRCA-R-1 and the ms-TRCA-R-2 are designed following the TDCA, i.e., After computing the spatial filters, the ms-TRCA-R-1 and the ms-TRCA-R-2 can follow the strategy shown in (31) to utilize the spatial filtering for the recognition.Moreover, Table III and Table IV show that combining spatial filters of multiple stimuli during the recognition can provide prominent performance in the TRCA-related methods.Therefore, the proposed TRCA-related methods could also follow this strategy of utilizing spatial filters shown in (32) for the recognition.Following the abbreviation rule of the TRCA-related methods mentioned in Section III-H, the methods of combining spatial filters are denoted as the corresponding ensemble methods, i.e., ms-eTRCA-R-1 and ms-eTRCA-R-2, respectively.In this study, to simplify the expressions of method names, ms-(e)TRCA-R-1 denotes ms-TRCA-R-1 and ms-eTRCA-R-1, and ms-(e)TRCA-R-2 denotes ms-TRCA-R-2 and ms-eTRCA-R-2.

B. New MsetCCA-Related Methods
Since the strategy of utilizing multi-stimulus calibration signals can improve the recognition performance as recommended in Section IV-D, III new methods, termed as Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
ms-MsetCCA-R-2, ms-MsetCCA-R-1 and ms-MsetCCA-R-3, can be obtained by extending the MsetCCA-R to its multi-stimulus cases, and applying the augmented calibration signals, i.e., The projection matrices and the combination matrices of the ms-MsetCCA-R- , and In addition, the ms-MsetCCA-R-3 further extends the ms-MsetCCA-R-1 following the TDCA, i.e., , and Due to the poor recognition performance of the MsetCCA and the MsetCCA-R, the proposed ms-MsetCCA-R-2, ms-MsetCCA-R-1 and ms-MsetCCA-R-3 does not use the strategy of utilizing spatial filters in the MsetCCA and the MsetCCA-R.Following the strategy of utilizing spatial filters in the ms-eTRCA, the proposed three MsetCCA-related methods firstly computes common spatial filters by utilizing calibration signals of stimuli that have close stimulus frequencies, and then combine all spatial filters of all trials and all stimuli for the recognition.More specifically, the correlation coefficients are calculated by (29) with C. Performance Evaluation of New TRCA-and MsetCCA-Related Methods The recognition performance of the proposed TRCA-and MsetCCA-related methods is verified by the Benchmark Dataset.All settings of simulations and evaluations are the same as those shown in Section IV-C.Fig. 2 and Fig. 3 illustrate the performance comparisons of the TRCA-and MsetCCA-related methods, respectively.More performance  comparisons of the existing and new multi-stimulus methods are illustrated in Section S.IX of the supplementary material.
In Fig. 2, we can observe that the performance of both the proposed ms-eTRCA-R-1 and ms-eTRCA-R-2 methods is better than or equal to that of existing TRCA-related methods, especially in the case of short signal length.In the new TRCA-related methods, the proposed ms-eTRCA-R-2 can provide significantly better performance than existing TRCA-related methods on all signal lengths.Fig. 3 shows that all new MsetCCA-related methods, especially the proposed ms-MsetCCA-R-3, are significantly better than existing MsetCCA methods.The prominent performance of the new two TRCA-related methods and the new three MsetCCA-related methods indicates that the proposed LS framework is very helpful for finding the prominent strategies of computing and using spatial filters, and achieving the further improvements of existing methods.

VI. CONCLUSION AND FUTURE PERSPECTIVES
This study proposes a LS unified framework for SSVEP spatial filtering methods and conducts a systematic and methodological review of CA-based methods within this framework.The proposed LS framework reveals the interconnections among these spatial filtering methods and the regression model.Subsequently, this study also provides a comprehensive summary of the strategies employed in spatial filter computations with the help of the proposed LS framework.Based on their commonalities and the differences, the strategies of computing spatial filters are grouped into 5 distinct categories.Moreover, beyond the LS framework, the strategies of utilizing spatial filters are further condensed into 4 types.By comparing the performance of these strategies through simulations and integrating prominent strategies to design new methods, this study offers valuable insights into the strengths and weaknesses of these strategies.Overall, this study provides a comprehensive analysis of the CA-based SSVEP spatial filtering methods within the proposed LS framework, contributing a deeper understanding of their underlying relationships and future advancements of SSVEP recognition methods.We outline 3 potential avenues below for future investigations.
Section V demonstrates the efficacy of the proposed LS framework in reassembling recommended existing strategies to fill the research gap and promote algorithmic improvements.Besides using the existing strategies, the proposed LS framework also offers insights into the development of novel computation strategies.For instance, one potential modification involves the combination matrices, L E and L K .Their current designs treat all trials in all stimuli equally, despite the temporal variability characteristics of EEG signals, such as different spatial distributions and the signal qualities across trials.Investigating the inclusion of trial weights in the combination matrices requires further studies.Another potential modification relates to the projection matrices, P E and P K .Fig. 1 illustrates that the strategy of utilizing P X m in P E is only applied in the itCCA and the eCCA, indicating unexplored potential.Specifically, there remains a research gap in terms of how to construct and employ SSVEP template signals for feature extraction in SSVEP recognition.
The proposed LS framework describes the relationships between the spatial filtering computations and the regression problem.Owing to the flexibility of the proposed LS framework in comparison to the existing GEP framework, the non-linearity and the regularization on the spatial filtering computations can be easily further investigated in the future studies.According to the generalized RRR model, the proposed LS framework shown in (2) can be modified to extend the CA-based spatial filtering methods into their non-linear form with regularization, i.e., arg min W,M (E) WM T − (T) 2 F + ρR (W) , where (•) and (•) denote the non-linear kernel methods of EEG features, ρ is the regularization parameter, and R (•) presents the regularization function.These studies would contribute to the explorations of the non-linear characteristics of EEG signals and address the channel selection issue in EEG signal analysis.
Recently, there has been a growing interest in employing the ANN for the SSVEP recognition [3], [38].However, the performance of existing ANN models in SSVEP recognition remains limited.This limitation can be attributed to a lack of sufficient considerations of the knowledge and specific characteristics associated with SSVEP signals within these ANN models.Although incorporating knowledge and signal characteristics from conventional methods into ANN models has yielded promising results in other BCI paradigms [43], [44], such approaches are rarely utilized in SSVEP recognition.By establishing a connection between the knowledge-based spatial filtering methods and the data-driven regression methods, the proposed LS framework demonstrates great potential for incorporating strategies in spatial filtering methods into ANN models, and further envisions a promising future for developing the knowledge-embedded ANN-based models in the SSVEP recognition.

Fig. 1 .
Fig. 1.Categorizations of the spatial filtering computations in the CA-based methods.Regions with green, red, purple and blue colors denotes the categories of computation strategies based on how to concatenate data in Z, which data is applied for the spatial filtering computation, how to combine intra-class features, and how to extract intra-class features.The arrows show the developments from simple to complex computations in the category of the computation strategy based on how to extract inter-class features, where the solid arrows show the developments of each classes, and the dotted arrows show the developments cross classes.The methods with the right superscipt a are proposed in this study.
m where W n,m denotes the spatial filter of the n-th trial for the m-th stimulus.The MsetCCA and the MsetCCA-R firstly apply the spatial filters to the calibration signals of each trial, and combine the filtered calibration signals together as the new template signal, which can be expressed as 2 and the ms-MsetCCA-R-1 are designed following ms-(e)TRCA-R-2 and the ms-(e)TRCA-R-1, i.e., (e)TRCA-R-2) E

Fig. 2 .
Fig. 2. Recognition performance comparisons between the existing prominent and new TRCA-related methods in the Benchmark Dataset.(a) Recognition accuracy performance.(b) ITR performance.

Fig. 3 .
Fig. 3. Recognition performance comparisons between the existing and new MsetCCA-related methods in the Benchmark Dataset.(a) Recognition accuracy performance.(b) ITR performance.

TABLE I TABLE OF KEY
NOTATIONS