The EEG signals collected by the device are generally two-dimensional data after pre-processing. N_c is denoted as the number of channels and N_t as temporal points. Input data X contains N samples x₁,x₂,x₃ …… x_N ∈ R^D, where D = N_tN_c represents the number of features in each sample. Then, X = [x₁,x₂,x₃, … x_N]^T ∈ R^{N × D} and y = [y₁,y₂,y₃, …, y_N]^T ∈ R^N represent the corresponding labels, where y_i ∈ {1,−1} is the class label. Its mathematical model can be expressed linearly as follows:

where w = [w₁,w₂,w₃ …… w_D]^T is a learnable weight vector, ε is noise, and X can be replaced by Φ(X) expressed in the form of a kernel function. Assume ε ∼ 𝒩 (0, σ²I_N), then y ∼ 𝒩 (Xw, σ²I_N) and its probabilistic framework is.

The RSBSBL adds the symmetric positive definite matrix in the variance term of the distribution that w obeys. The input data of one channel are regarded as a block. So, for the mathematical model (1), assume that w_b(∀b) is mutually independent and Gaussian distributed.

where w_b containing several w_i is bth block of w, γ_b is a non-negative scalar that controls the variance of w_b, B_b is a positive definite matrix reflecting the intra-block correlation, and N_b is the number of blocks. Since the features of a channel are considered to be a block, N_b = N_c.

In our case of EEG signal, b is the index of channels. In a channel of EEG signal with corresponding weight w_b, it is assumed that all its feature weights share the same γ_b to control the variance of their distribution, and B_b controls the intra-block correlation.

In this case, we express the prior of w as p(w|γ, B) ∼ 𝒩 (0, Σ₀), where Σ₀ is

the posterior probability has been calculated by the Bayesian rule,

and the corresponding variance and mean of the posterior probability density p(w|y, σ², γ, B) ∼ 𝒩 (μ_w, Σ_w) can be described as

When N ≥ D, the Eqs (6) and (7) are suitable because the maximum size of the inverse matrix is D in this case. Now, we give the iterative ways when N < D. According to the matrix inversion formula and the matrix identity.

we replace the Eqs (6) and (7) with the following equations:

To find the iterative equation of the parameters Θ = {γ, B, σ²}, the expectation–maximization (EM) is used to maximize log p (y|Θ). The Q function is.

The first term of the Q function is related to σ² and the second term is related to γ and B. Then, we can get the parameters iteratively by maximizing the Q function.

where *old represents the hyperparameter in the previous iteration, and the superscript b of μwb and Σwb indicates the bth block in μ_w and Σ_w with the size of d_b × 1 and d_b × d_b (d_b is the number of elements in w_b).

The potential similarity exists in the adjacent electrode signals considering the volume conduction effects in the brain (Hassan and Wendling, 2018). We assign the same B_re for channels with close locations for regional smoothing, and the region G_re contains g_re channels. As shown in Figure 1, all the channels are divided into 13 regions by position, and each region contains at least three channels. B_re is the average of blocks in region re (re ∈ [1, 13]).

We use a first-order Auto-Regressive (AR) process to model the intra-block correlation. Many applications have used the AR process to express it (Zhang and Rao, 2011; Zhang et al., 2013; Yin et al., 2020). Thus, to find a symmetric positive definite matrix to approximate B, it can be constrained to the following form of the Toeplitz matrix.

Empirically calculate r=m1m0, where m₀ is the average of the main diagonal of B_re and m₁ is the average of the main sub-diagonal.

Regarding the feature extracted from the same channel as a block, we perform RSBSBL to get the weight vector of features and design a channel selection based on the weight vector as Algorithm 1.

As shown in Algorithm 1, the parameters are initialized, and the shear threshold τ is set. Then, from Line 3 to Line 12, the algorithm iteratively solves BSBL and prunes the γ. Line 4 to Line 8 decide the calculation of Σ_w, μ_w, so that the large time cost caused by finding the inverse matrix of a large-size matrix can be alleviated. The parameters are updated on Line 9 and Line 10. Figure 2 illustrates the relationship between the parameters in a single iteration, where the parameters calculated simultaneously have the same color. The solid line indicates the passing relationship between the parameters of this iteration, and the dashed line indicates the passing relationship between the parameters of this iteration and the next iteration. After the parameters are calculated, in order to achieve the sparse block effect, make γ_b to 0 when γ_b is less than the threshold τ. Then, it comes into the next iteration until the convergence criterion is satisfied. Line 13 automatically selects the channels with γ_b greater than the shear threshold τ. Finally, the algorithm returns the selected channel and the corresponding weight vector.

The off-diagonal matrix B makes the weights w in the same block relevant in distribution. It means that the correlation of the features from the same channel can be reflected during the process. Moreover, the components of the temporal correlation of different channels in close locations are the same because the B_re of channels in the same region are shared. The sparsity of weights will form the units of channels. The features from one channel share the same weight distribution whose variance is controlled by γ. For practicality, up to five channels are removed in a single iteration when making a channel selection.

Three datasets were used in this study to validate the proposed method. DS1 is BCI Competition II dataset IIb (one participant) (Blankertz et al., 2004) and DS2 is BCI Competition III dataset II (two participants) (Blankertz et al., 2006). DS3 is the EEG signal collected in our lab (12 participants). The stimulus numbers for each participant of the above three datasets are shown in Table 1.

DS1 and DS2 provided by the BCI Competition are public datasets and follow the same experimental paradigm of Farwell and Donchin, as shown in Figure 3. In a six-by-six character matrix containing 26 characters and 10 numbers, participants were asked to focus on a specified character in each trial (a trial is a set of stimuli that can support the output of a recognized character). They could do this by mentally counting the target stimuli’ number of flashes (intensifications). The paradigm continuously intensified and randomly scanned all rows and columns of the matrix at a rate of 5.7 Hz. Each row and column in the matrix was randomly intensified for 100 ms and was left blank for 75 ms. DS1 contained 42 training characters and 31 testing characters. The training set of DS2 contained 85 characters, and the testing set contained 100 characters. A trial for each character had 15 epochs to apply reliable spelling, and each epoch was comprised of 12 intensifications. Both datasets were collected using a 64-channel cap, filtered by 0.1–60 Hz, and digitized at a sampling rate of 240 Hz. DS1 and DS2 can be downloaded from the websites: http://www.bbci.de/competition/ii/ and http://www.bbci.de/competition/iii/.

DS3 was collected in our lab. Its paradigm was similar to the BCI Competition. It contained 26 characters and 10 numbers. DS3 consisted of 12 participants who were graduate students between the ages of 20 and 26 years, with normal or corrected-to-normal vision. The experiments used a 64-channel wireless EEG acquisition system (Neuracle, NeuSen W series, 59 EEG, 4EOG, 1ECG) to acquire data at the sampling rate of 1,000 Hz. In the paradigm, each row and column in the six-by-six matrix was randomly intensified for 80 ms and kept extinguished for 80 ms. A trial for each target character included four epochs, and each epoch had 12 intensifications. Participants were required to spell 36 characters. We randomly selected 18 characters as the training dataset and the rest as the test dataset.

Considering that some channels contain less task-relevant information but more noise, it is vital to use a reasonable method to select the most effective channels. This study compares the proposed RSBSBL with two empirical channel sets (Set 1 and Set 2) (Krusienski et al., 2008), LASSO, GLASSO, and SBL in the case of using the same pre-processing process and classifier. Set 1 includes Fz, Cz, Pz, Oz, PO7, and PO8. Set 2 includes Fz, FCz, Cz, C3, C4, CPz, Pz, P3, P4, P7, P8, POz, PO3, PO4, PO7, PO8, Oz, O1, and O2.

Figure 3 shows the diagram of the data processing framework, which includes three main parts: (1) pre-processing, (2) channel selection, and (3) classification. DS1 and DS2 shared the same pre-processing: bandpass filtering of data from 0.5 to 20 Hz and downsampling by a factor of 5. Then, the sampling rate of the data was 48 Hz. We intercepted 0–667 ms after each stimulus as the primary analysis objective was to obtain 32 sampling points for each stimulus. For the DS3, the 59-channel dataset that went through 0.5–20 Hz bandpass filtering was down-sampled to 50 Hz and the data segment from 0 to 600 ms was taken after stimulation to obtain 30 sampling points for each stimulus. Thus, denoting the number of channels as N_c and number of signal sampling points as N_t, a 1 × D feature matrix was obtained for each stimulus, where D = N_tN_c. A feature matrix was labeled “1” only if the corresponding stimulus belongs to the row or column of the target characters. Otherwise, the label was assigned to “0.”

The typical classification methods of P300 include traditional machine learning methods and neural network-based methods. Traditional machine learning can achieve outstanding performance with less complexity. This study regarded BLDA as a unified classifier for different channel selection algorithms.

The optimal combination of parameters was determined by a 10-fold cross-validation. There were two modes of the selected channel number in the experiment for the channel selections: automatic and fixed. When the channel number was determined automatically, we used a threshold to determine the channel number. For LASSO and SBL, the absolute values of the feature weights in one channel were summed up to represent the importance of the channel. The threshold equaled the mean minus 0.5 times the standard deviation of the channel importance values, and the channels with importance values higher than the threshold were selected. As for GLASSO and RSBSBL, automatic channel selection had been enabled in the methods. When the number of selected channels is fixed (M channels were selected), we used the same way to evaluate each channel. For all the four methods, the absolute values of the feature weights w of each channel were summed, and the top M channels were selected in descending order.

We used character recognition accuracy to evaluate the performance of a classification. The character recognition accuracy is defined as follows:

where C_{test_total} represents the total number of characters in the test dataset, and C_{test_correct} is the sum of all the correctly predicted characters. Besides, to evaluate the significance of performance difference, we introduced a non-parametric statistical hypothesis test, the Wilcoxon signed-rank test. The Wilcoxon signed-rank test can be used as an alternative to the paired t-test for matched pairs when the population cannot be assumed to be normally distributed. The significance of the pairs can be confirmed when the corresponding p-value is less than 0.05.

Channel selection is supposed to reserve channels with more helpful information and exclude the channels with more noise. According to the data processing, we chose a unified classifier to verify the performance of different methods for a fair comparison. In Table 2, we compared the character recognition accuracy of each method on the three datasets, and the number of selected channels was automatically determined as described in section “Parameter Setting.” Set 1 and Set 2 are empirical subsets of channels (Set 1 contains 6 channels and Set 2 contains 19 channels). The best results were marked in bold, and the number of channels selected for each participant is presented in the corresponding parentheses.

For DS1, RSBSBL selected the minimum number of channels when the classification accuracy of all the methods was 100%. For DS2, RSBSBL had the highest average accuracy, 97.50%, which was 1.00% higher than the second-ranked GLASSO. Although SBL selected fewer channels than others, the average recognition accuracy was 95.00%.

For DS3, RSBSBL as a channel selection method could bring higher accuracy with BLDA in 11 participants among 12 and got 88.43% average accuracy by eliminating insufficient data than using all channels. It outperformed the second-ranked GSBL on an average by 7.88% and selected the fewest channels as 17. We evaluated the significance of the classification performance of DS3 via the Wilcoxon signed-rank test and found that the proposed method performed significantly better than others (RSBSBL vs. LASSO: p = 0.005 < 0.05; RSBSBL vs. GLASSO: p = 0.013 < 0.05; RSBSBL vs. SBL: p = 0.003 < 0.05).

To further compare the effectiveness of the four methods, we compared the recognition results of the algorithms when M channels were selected (M = [4, 8, 12, 16]). Top M channels were selected by ranking the corresponding channels according to the sum of the absolute values of the feature weights. The classifiers were retrained with the data with the selected channel. It was supposed that the number of channels M′ automatically selected by the method was less than the value of M. In that case, the latest deleted M-M′ channels are added according to the order in which they were deleted during the iteration of the method.

Figure 4 shows the accuracy of each method on DS1, DS2, and DS3, with the horizontal coordinates of the bars indicating the selection of the top M channels. For DS1, the accuracy of all the methods was the same except that the accuracy of SBL was 96.77% when eight channels were selected, and it was lower than others. For DS2, SBL and RSBSBL obtained better performance with 80% average recognition accuracy when four channels were selected. When 8, 12, and 16 channels were selected, GLASSO obtained an average recognition accuracy of 78.5, 84.5, 91, and 92%, respectively, and RSBSBL obtained a better performance of 80, 85.5, 91.5, and 93.5%, respectively. For DS3, GLASSO obtained average recognition accuracy of 73.61, 75.93, 75.46, and 79.63% when 4, 8, 12, and 16 channels were selected, respectively. Moreover, RSBSBL obtained the best performance of 74.07, 82.87, 80.09, and 80.56%, respectively. The average recognition accuracies of LASSO, GLASSO, and RSBSBL on DS3 with M = 16 were 77.31, 79.63, and 80.56%, respectively. The results of experiments with the fixed number of selected channels revealed that the feature weights generated by RSBSBL could provide more reasonable guidelines for the channel selection.

We counted the selected channels at the same location and used it to describe the number of times a channel has been selected in the dataset. If 6 of the 12 participants’ selected channels contain Pz, then the contribution value of the channel corresponding to the Pz electrode is 6. Figure 5 indicates the scalp distributions of the contribution value of channels on DS1, DS2, and DS3. The color changes from red to blue, indicating that the channel was selected less often. As shown in Figure 5, when the number of selected channels was small (M = 4, 8), RSBSBL selected the occipital and parietal electrodes more often. It shows that, in addition to the P300 potential, the early visual components also contribute to a classification in the paradigm (Blankertz et al., 2011).

In RSBSBL, γ_b smaller than the threshold τ was set to zero, indicating that τ determines the pruning strength. We analyzed the change in the number of channels selected and the recognition results when τ is assigned different values in the range 10⁻⁸ to 10⁻¹. The recognition accuracy of each participant varying with τ was normalized to highlight the location of the optimal threshold. Figure 6 illustrates the effect of the threshold on the proposed method. The x-axis indicates the number of selected channels, the y-axis indicates the value of τ, and the z-axis indicates the participant ID. The color changes from red to blue, indicating that the point corresponds to a higher to lower normalized accuracy.

As shown in Figure 6, the number of channels selected by each participant increased as the threshold value decreased. When the threshold was less than or equal to 10⁻⁶, the number of selected channels was the original number in the dataset, and the algorithm loses the ability to select the channels automatically. Therefore, 10-fold cross-validation can be used to select the optimal parameter values in the range of 10⁻⁶ to 10⁻¹. From the curves corresponding to P3.2, P3.3, P3.7, and P3.12, using selected channels can obtain better recognition accuracy than using all the channels, which proves that channel selection can remove weak task-relevant and noisy channels to improve the classification accuracy.

Fabiani et al. (1987) confirmed that the visual P300 paradigm should at least include Fz, Cz, Pz electrodes signed as the 10–20 international electrode system. Krusienski et al. (2008) and McCann et al. (2015) made sure that Fz, Cz, Pz, Oz, PO7, and PO8 corresponded to the parietal and occipital regions of the brain that take a significant part in the recognition of P300 signals. In Table 2, Set 1 and Set 2 represent two empirical channel subsets. Set 1 includes Fz, Cz, Pz, Oz, PO7, and PO8. Set 2 includes Fz, FCz, Cz, C3, C4, CPz, Pz, P3, P4, P7, P8, POz, PO3, PO4, PO7, PO8, Oz, O1, and O2. It can be seen that for many participants (P2.1, P2.2, P3.1, P3.2, P3.5, P3.6), the character recognition accuracy was lower when the empirical channel subsets were used. The empirical selection may not include some channels that contribute to the classification. The channels assumed to reflect visual components and also some frontal channels contribute to the classification for some participants. It also indicates the lower robustness of the empirical channel subset. In Figure 5, the scalp mapped according to channel selection of RSBSBL could be observed with high values in Pz, P3, P4, O1, O2, Oz, PO7, PO8, and POz regions. These electrodes are very similar to the abovementioned electrodes, which are closely related to the visually induced ERPs. The P1, N1, and N2 components are mainly concentrated in the parietal and occipital regions. And the central distribution of P2 and P3 is elongated along the midline electrodes (Blankertz et al., 2011). It can be assumed that a multitude of ERP components is affected by attention to the target and utilized by classifiers rather than just the P300 (Treder and Blankertz, 2010). In addition, it can be found from Figures 5, 6 that many participants in DS3 had poorer classification using full-channel data compared to DS1 and DS2, and their topographic maps select more frontal channels when M = 8, 12,16. This phenomenon may be due to the effect of eye artifacts and noise during the experiment.

Table 2 and Figure 3 show the superiority of RSBSBL in channel selection. When the number of channels was determined automatically, the proposed method achieved the highest average recognition accuracy of 100, 97.5, and 88.43% for DS1, DS2, and DS3, with the lowest average number of channels on DS1 and DS3. The RSBSBL achieved better performance than the compared methods when selecting the channels with the fixed number, and the average accuracies of 90.21, 80, and 74.07% were obtained with the top four selected channels on the three datasets.

To verify the performance of RSBSBL, we compared the proposed method with the state-of-the-art developments in recent years on DS2, as shown in Table 3. Most of them are based on evolutionary computational algorithms (Kee et al., 2015; Khairullah et al., 2020; Tang et al., 2020; Martinez-Cagigal et al., 2022). The channel selection methods and classifiers used in each study are shown in the table.

The shear threshold τ significantly impacted the final results, so cross-validation was required to determine the optimal parameters. According to the analysis of parameter sensitivity, as shown in Figure 6, the recommended threshold selection range was [10⁻⁶, 10⁻¹]. Besides, Figure 6 reflects the variation of character recognition accuracy with the shear threshold for each participant. Compared with others, P3.2, P3.3, P3.9, P3.10, and P3.12 cannot achieve the best recognition accuracies with the full channels, which implies that the EEG signals of these participants have more channels with noise, and these channels are not conducive to signal classification. As shown in Table 2, when determining the number of channels automatically, RSBSBL can achieve the best recognition accuracies of them with the corresponding number of selected channels of 14, 13, 15, 15, and 13, respectively. It confirms that RSBSBL can remove unfavorable channels and improve the recognition accuracies.

To verify the effectiveness of regional smoothing, we conducted further controlled experiments on the three datasets, and the results are shown in Table 4. Case 1 represents that B is a unit matrix, implying that no temporal correlation is considered. Case 2 has the same B for all blocks, indicating that all channels share the same B. Case 3 has a different B matrix for each block, showing that regional smoothing is no longer done. The comparison between Case 3 and Case 1 in Table 4 illustrates the improvement of the model due to temporal correlation. The comparison between our algorithm and Cases 3 and 1 indicates the improvement brought by region smoothing. The“*” in Table 4 represents a significant difference in our method after Wilcoxon signed-rank test (RSBSBL vs. Case 1: p = 0.015 < 0.05; RSBSBL vs. Case 2: p = 0.031 < 0.05; RSBSBL vs. Case 3: p = 0.124).

As described in sections “Data Descriptions” and “The Framework of Data Processing,” for DS1 and DS2, N_t = 32 and N_c = 64 after pre-processing, then we can get a 1 × D (D = N_tN_c = 2048) vector for each stimulus. As shown in Table 1, in the training datasets of DS1 and DS2, the total number of stimuli was 7,560 and 15,300, which is larger than the number of features D. For DS3, N_t = 30 and N_c = 59 after pre-processing, then the feature is a 1 × D (D = N_tN_c = 1770) vector. In Table 1, in the training datasets of DS3, the total number of stimuli was 864, which is smaller than the number of its features.

In a preliminary study, we found that inappropriate iterations can make the algorithm to have a large time cost [e.g., using equations (10) and (11) on DS1 and DS2]. Therefore, a strategy of automatic selection of the iteration method is used to avoid this problem. In Figure 7, we analyze the variation of the matrix inversion run-time when the size of the matrix increases (the matrix is a square matrix). In the left part, the horizontal axis represents the size of the square matrix. The vertical axis is the value after taking the logarithm of the time, and the actual time (s) is also indicated in the figure. It can be noticed that the time spent on matrix inversion is more than 1 s when the matrix size is larger than 3,000 × 3,000. Therefore, we consider that the method may not be suitable for data with numbers of features and samples larger than 3,000. Of course, this problem can be solved by reducing the number of features and optimizing the iteration steps. The right bar in Figure 7 indicates the average time cost of the proposed method on the three data sets, which is acceptable.

The sparse Bayesian algorithm can make the sparsity of the algorithm change by changing the prior distribution of w (Tipping, 2001). Zhang et al. (2015) used the Laplace distribution instead of the traditional Gaussian distribution for the classification of P300 signals using SBL. Therefore, RSBSBL can change the prior of the weights to make the sparsity stronger in the future, such as the Gamma distribution. The proposed method used the EM algorithm for iteration, and there is still room for improvement in the computational speed. In the future, we will also explore the suitability of the proposed method for other ERPs.