An extended binary subband canonical correlation analysis detection algorithm oriented to the radial contraction-expansion motion steady-state visual evoked paradigm

Abstract

The radial contraction-expansion motion paradigm is a novel steady-state visual evoked experimental paradigm, and the electroencephalography (EEG) evoked potential is different from the traditional luminance modulation paradigm. The signal energy is concentrated chiefly in the fundamental frequency, while the higher harmonic power is lower. Therefore, the conventional steady-state visual evoked potential recognition algorithms optimizing multiple harmonic response components, such as the extended canonical correlation analysis (eCCA) and task-related component analysis (TRCA) algorithm, have poor recognition performance under the radial contraction-expansion motion paradigm. This paper proposes an extended binary subband canonical correlation analysis (eBSCCA) algorithm for the radial contraction-expansion motion paradigm. For the radial contraction-expansion motion paradigm, binary subband filtering was used to optimize the weighting coefficients of different frequency response signals, thereby improving the recognition performance of EEG signals. The results of offline experiments involving 13 subjects showed that the eBSCCA algorithm exhibits a better performance than the eCCA and TRCA algorithms under the stimulation of the radial contraction-expansion motion paradigm. In the online experiment, the average recognition accuracy of 13 subjects was 88.68% ± 6.33%, and the average information transmission rate (ITR) was 158.77 ± 43.67 bits/min, which proved that the algorithm had good recognition effect signals evoked by the radial contraction-expansion motion paradigm.

Keywords

steady-state visual evoked potentials brain–computer interface radial contraction-expansion motion paradigm binary subband canonical correlation analysis extended binary subband canonical correlation analysis

1 Introduction

Brain–computer interface (BCI) system based on steady-state visual evoked potential (SSVEP) has attracted extensive attention in the field of BCI because of its high information transmission rate (ITR) and low training cost [1]. The widely used SSVEP paradigm evokes electroencephalography (EEG) signals in luminance modulation. However, this paradigm is easy to cause visual fatigue [2] and contains the risk of evoking photosensitive epilepsy. Additionally, the harmonic of the EEG signal evoked by the luminance modulation paradigm is prominent, which leads researchers to avoid the overlap between the fundamental frequency of a target and the harmonic of other targets when designing stimulation frequency, thus limiting the selection of stimulation frequency [3]. Han et al. studied a paradigm with radial contraction-expansion motion as the stimulus-evoked form and developed a high-interaction BCI based on flicker-free SSVEP [4]. This paradigm uses the feature that the human visual system is susceptible to motion perception to evoke EEG signals. The luminance of the stimuli in the paradigm remains constant, reducing the subjects’ discomfort caused by target flickering [5, 6]. In addition, the EEG evoked by this paradigm has low harmonic content and more frequencies in the range of sensitive frequencies [7] that can be used to stimulate the target than the luminance modulation paradigm. Therefore, the radial contraction-expansion motion paradigm is an excellent alternative to the traditional luminance modulation paradigm.

Majority of existing research on the radial contraction-expansion motion paradigm includes comparative studies of EEG characteristics evoked by different paradigms and practical research on the radial contraction-expansion motion paradigm. The identification algorithms used are mostly canonical correlation analysis (CCA) [8]. As a statistical method for measuring the potential correlation between two multidimensional variables, CCA is widely used in SSVEP-BCI [9]. However, in standard CCA-based methods, the canonical correlation value between SSVEP and the sin-cos reference signal tends to decrease as the flicker frequency increases, leading to a decrease in the accuracy of SSVEP detection at higher frequencies [10]. In addition, due to the different characteristics of evoked signals and the EEG signal intensity evoked by the radial contraction-expansion motion paradigm being weaker than that of the luminance modulation paradigm [11], the recognition effect of trained algorithms with excellent recognition performance under the luminance modulation paradigm is not satisfactory under the radial contraction-expansion motion paradigm, such as task-related component analysis (TRCA) [12], extended canonical correlation analysis (eCCA) [13].

The power spectral density of EEG signals is usually a decreasing function, and the power is low at higher frequencies [14]. Binary subband canonical correlation analysis (BsCCA) filters EEG data into two subbands. The first subband covers all target frequencies, and the second subband contains high-frequency components that need to be enhanced. It adjusts high- and low-frequency components through subband weights to improve the recognition accuracy of high-frequency stimuli [15]. Therefore, this study proposes an extended binary subband canonical correlation analysis (eBSCCA) detection algorithm oriented to the radial contraction-expansion motion paradigm in response to the above problems. The algorithm uses the BsCCA algorithm to improve the recognition accuracy of high-frequency stimuli and combines the dynamic window method to improve the ITR of the system. The experimental results show that, compared with the TRCA and eCCA algorithms, the eBSCCA algorithm proposed in this paper can effectively improve the recognition accuracy and ITR of EEG signals.

The structure of this paper is as follows: the second section introduces the characteristics of EEG signals evoked by radial contraction-expansion motion paradigm and the proposed algorithm; the third section introduces the experimental process, parameter selection, comparison algorithm, and algorithm performance evaluation criteria; the fourth section shows the parameter calculation results, offline experiment results, and online experiment results of the dataset; and in Section 5, the experimental results and paradigm performance are discussed; the sixth section is the conclusion of this paper.

2 Methods

2.1 Signal characteristics

In the radial contraction-expansion motion paradigm, the contraction-expansion motion of the target was realized by periodically reducing and increasing the radius of concentric circles with different radii. The cosine function adjusts the change process of the circle radius. The radial contraction-expansion process of a circle can be expressed by Eq. (1):

R = r + A cos (2 π f_{c} t + φ_{0})

where r represents the initial radius of the circle; R represents the actual radius of the circle in the process of radial contraction expansion; A represents the amplitude of the contraction-expansion movement of the circle; f _c represents the movement frequency; and φ ₀ represents the initial phase.

As shown in Fig. 1, when the phase φ = 2πf _c t + φ ₀ gradually increases from 0 to π, the circle’s radius changes from large to small, and the target contracts; when the phase gradually increases from π to 2π, the circle’s radius increases from small, and the target expands.

The frequency at which the movement direction changes is called the flip frequency denoted as f, f = 2 f _c [16]. The frequency of the EEG signal evoked by the radial contraction-expansion motion paradigm corresponds to the flip frequency f. Only the signal’s fundamental frequency was significant in the spectrum, and the second and higher harmonics were insignificant [17]. Fig. 2 was derived from the offline data analysis of each subject in this study. As shown in Fig. 2(A), the EEG signal energy evoked by the radial contraction-expansion motion paradigm is concentrated in the fundamental frequency component. As shown in Fig. 2(B), the average signal-to-noise ratio (SNR) of the EEG signal evoked by high-frequency is lower than that of the low-frequency part and showed a downward trend. The SNR of each frequency in Fig. 2(B) is calculated using Eq. (2) [18]:

Fig. 1

Schematic diagram of the contraction-expansion process of the stimulus target.

SNR = \frac{tr ({\hat{X}}_{k} \cdot Φ^{Η} \cdot Φ \cdot {\hat{X}}_{k}^{H})}{tr [{\hat{X}}_{k} \cdot (I_{N} - Φ^{Η} \cdot Φ) \cdot {\hat{X}}_{k}^{H}]}

where the constraint $Φ \cdot Φ^{Η} = I_{2 N_{h}} \cdot Φ \in ℝ^{2 N_{h} \times N_{s}}$ is defined as the 5-fold complex sinusoidal template of the stimulation frequency, N _h represents the number of harmonics (N _h = 5), and N _s represents the number of sampling points. ${\hat{X}}_{k} \in ℝ^{N_{c} \times N_{s}}$ represents the EEG template corresponding to the kth target frequency, obtained by averaging the offline data by block; and N _c represents the number of channels. The template data ${\hat{X}}_{k}$ contains 9 EEG channels of the occipital area (PO3、PO4、PO5、PO6、PO7、PO8、O_z、O1、O2). In Fig. 2(B), the ordinate value is the average SNR of the 13 subjects.

Fig. 2

Signal characterization analysis of offline data. (A) The average unilateral amplitude spectrogram (11, 17, and 19 Hz), the abscissa represents the response frequency, and the ordinate represents the unilateral amplitude. (B) The average SNR of each stimulation frequency. The abscissa represents the stimulation frequency, and the ordinate represents the average SNR of the EEG response at this frequency.

2.2 Proposed algorithm

This paper proposes an extended BsCCA detection algorithm for the radial contraction-expansion motion paradigm. The CCA algorithm is improved based on the characteristics of the EEG evoked by this paradigm, which is used in training scenarios. The algorithm uses the BsCCA algorithm to improve the recognition accuracy of high-frequency stimuli [19] and uses the dynamic window method to improve the ITR further.

2.2.1 BsCCA

During the experiment, the SNR of the EEG signal in the high-frequency part was lower than that in the low-frequency part, which decreased the recognition accuracy of high-frequency components of some subjects. The BsCCA algorithm proposed by Islam et al. can improve the recognition accuracy of high-frequency stimuli [19]. As shown in Fig. 3(A), the BsCCA algorithm proposed in Ref. [19] was applied to the scene without training. However, this study further applied it to the scene with training. Figs. 3(B) and (C) show the flow of the BsCCA algorithm with training.

Fig. 3

Schematic diagram of BsCCA algorithm. (A) Schematic diagram of the flowchart of the untrained BsCCA algorithm. (B) Schematic diagram of the training process of the trained BsCCA algorithm. (C) Schematic diagram of the detection process of the trained BsCCA algorithm.

In Fig. 3, $Y_{f} \in ℝ^{2 N_{h} \times N_{s}}$ represents the sine-cosine reference signal, $X \in ℝ^{N_{c} \times N_{s}}$ represents the test data, α ₁ and α ₂ represent the weights of the all-pass subband and the high-frequency optimized subband, respectively. The sin-cos reference signal Y _f is set to [9]:

ϒ_{f} = [\begin{matrix} \sin (2 π f \frac{0}{f_{s}}) & \dots & \sin (2 π f \frac{N_{s} - 1}{f_{s}}) \\ \cos (2 π f \frac{0}{f_{s}}) & \dots & \cos (2 π f \frac{N_{s} - 1}{f_{s}}) \\ ⋮ & ⋱ & ⋮ \\ \sin (2 π N_{h} f \frac{0}{f_{s}}) & \dots & \sin (2 π N_{h} f \frac{N_{s} - 1}{f_{s}}) \\ \cos (2 π N_{h} f \frac{0}{f_{s}}) & \dots & \cos (2 π N_{h} f \frac{N_{s} - 1}{f_{s}}) \end{matrix}]

where f represents the stimulation frequency; and f _s represents the 1000 Hz sampling rate.

Fig. 3(A) shows the flow of the untrained BsCCA algorithm: after the test data X passes through two band-pass filters, an all-pass subband X ⁽¹⁾ and a high-frequency optimized subband X ⁽²⁾ are obtained. The two subbands perform CCA with Y _f respectively to obtain the maximum correlation coefficient $ρ_{k_{1}}^{(1)}$ and $ρ_{k_{1}}^{(2)}$ , then $ρ_{k_{1}}^{(1)}$ and $ρ_{k_{1}}^{(2)}$ are weighted by subband weights α ₁ and α ₂ to obtain the weighted sum $ρ_{k_{1}}$ .

Fig. 3(B) shows the training process of the trained BsCCA algorithm: after the EEG template ${\hat{X}}_{k}$ passes through two band-pass filters, an all-pass subband ${\hat{X}}_{k}^{(1)}$ and a high-frequency optimized subband ${\hat{X}}_{k}^{(2)}$ are obtained. ${\hat{X}}_{k}^{(1)}$ and Y _f perform CCA to get the spatial filter $w_{y_{f} {\hat{X}}_{k}}^{(1)} \in ℝ^{1 \times 2 N_{h}}$ and $w_{{\hat{X}}_{k} y_{f}}^{(1)} \in ℝ^{1 \times N_{c}}$ , ${\hat{X}}_{k}^{(2)}$ and Y _f perform CCA to get the spatial filter $w_{y_{f} {\hat{X}}_{k}}^{(2)} \in ℝ^{1 \times 2 N_{h}}$ and $w_{{\hat{X}}_{k} y_{f}}^{(2)} \in ℝ^{1 \times N_{c}}$ . According to the CCA algorithm, the formula for solving the spatial filter is shown in Eq. (4) [9]:

[w_{ϒ_{f} {\hat{X}}_{k}}^{(n)}, w_{{\hat{X}}_{k} ϒ_{f}}^{(n)}] = \underset{w_{ϒ_{f} {\hat{X}}_{k}}^{(n)}, w_{{\hat{X}}_{k} ϒ_{f}}^{(n)}}{argmax} (\frac{E [{(w_{ϒ_{f} {\hat{X}}_{k}}^{(n)} \cdot ϒ_{f})}^{T} \cdot (w_{{\hat{X}}_{k} ϒ_{f}}^{(n)} \cdot X^{(n)})]}{\sqrt{E [{(w_{ϒ_{f} {\hat{X}}_{k}}^{(n)} \cdot ϒ_{f})}^{T} \cdot (w_{ϒ_{f} {\hat{X}}_{k}}^{(n)} \cdot ϒ_{f})] E [{(w_{{\hat{X}}_{k} ϒ_{f}}^{(n)} \cdot X^{(n)})}^{T} \cdot (w_{{\hat{X}}_{k} ϒ_{f}}^{(n)} \cdot X^{(n)})]}})

where n represents the subband index.

Fig. 3(C) shows the detection process of the trained BsCCA algorithm: the test data X passes through two band-pass filters to obtain an all-pass subband X ⁽¹⁾ and a high-frequency optimized subband X ⁽²⁾. For X ⁽¹⁾, the Pearson correlation coefficient $ρ_{k_{3}}^{(1)}$ of vectors $w_{y_{f} {\hat{X}}_{k}}^{(1)} \cdot Y_{f}$ and $w_{{\hat{X}}_{k} y_{f}}^{(1)} \cdot X^{(1)}$ is obtained from Eq. (5) [20]. For X ⁽²⁾, the Pearson correlation coefficient $ρ_{k_{3}}^{(2)}$ of vectors $w_{y_{f} {\hat{X}}_{k}}^{(2)} \cdot Y_{f}$ and $w_{{\hat{X}}_{k} y_{f}}^{(2)} \cdot X^{(2)}$ is obtained from (5). $ρ_{k_{3}}^{(1)}$ and $ρ_{k_{3}}^{(2)}$ are weighted by subband weights α ₁ and α ₂ to obtain the weighted sum $ρ_{k_{3}}$ .

τ_{(w_{ϒ_{f} {\hat{X}}_{k}}^{(n)} \cdot ϒ_{f}, w_{{\hat{X}}_{k} ϒ_{f}}^{(n)} \cdot X^{(n)})} = \frac{cov (w_{ϒ_{f} {\hat{X}}_{k}}^{(n)} \cdot ϒ_{f}, w_{{\hat{X}}_{k} ϒ_{f}}^{(n)} \cdot X^{(n)})}{σ_{w_{ϒ_{f} {\hat{X}}_{k}}^{(n)} \cdot ϒ_{f}} σ_{w_{{\hat{X}}_{k} ϒ_{f}}^{(n)} \cdot X^{(n)}}}

where cov(·) represents the covariance between two variables; and σ represents the variable’s standard deviation.

2.2.2 eBSCCA

Fig. 4 shows the flowchart of the eBSCCA algorithm. The decision feature consists of three parts: $ρ_{k_{1}}$ , $ρ_{k_{2}}$ and $ρ_{k_{3}}$ . The calculation method is as follows:

$ρ_{k_{1}}$ : The calculation method is shown in Fig. 3(A).

$ρ_{k_{2}}$ : The spatial filter calculation method is shown in Eq. (4), and the Pearson correlation coefficient calculation method is shown in Eq. (5).

$ρ_{k_{3}}$ : The spatial filter calculation method is shown in Fig. 3(B), and the calculation method of the decision result is shown in Fig. 3(C).

2.2.3 Dynamic window

Traditional SSVEP recognition algorithms only analyze fixed-length data, such as CCA and filter bank canonical correlation analysis (FBCCA) [21]. These algorithms usually estimate the optimal data window length in offline datasets based on the performance of the BCI system. The fixed window length is then applied to online experiments [18]. However, due to the complexity and non-stationarity of spontaneous EEG signals, the optimal window length is different among different subjects [22]. In addition, the psychological status and environmental changes of subjects are also the influencing factors of the optimal window length. The dynamic window method can effectively reduce the above adverse effects and dynamically adjust the data length required by the algorithm while maintaining the high precision of the decision to achieve a higher ITR, thereby improving the performance of the BCI [23].

Fig. 4

Flowchart of eBSCCA algorithm.

In this study, the dynamic window method sets a minimum calculation window length T _low, a maximum calculation window length T _high, and a dynamic window threshold d. Each calculation of the eBSCCA algorithm will obtain a 35-dimensional correlation coefficient vector ρ = {ρ ₁, ρ ₂ ,…, ρ ₃₅ }, and the dynamic window threshold is the kurtosis value of the correlation coefficient vector ρ. The threshold reflects the steepness of the distribution of vector values. The greater the kurtosis value, the more significant the correlation coefficient values. The kurtosis calculation formula of vector ρ is shown in Eq. (6) [24, 25]:

d = \frac{Q \sum_{i = 1}^{Q} {(ρ_{i} - \bar{ρ})}^{4}}{{[\sum_{i = 1}^{Q} {(ρ_{i} - \bar{ρ})}^{2}]}^{2}}

where ρ_i represents the ith correlation coefficient; $\bar{ρ}$ represents the average value of the vector ρ ; and Q represents the number of stimulus targets. Fig. 5 shows the calculation flow of the dynamic window method.

Fig. 5

Flow chart of dynamic window method.

Algorithm 1 shows the pseudo-code of the eBSCCA algorithm.

Algorithm 1

eBSCCA algorithm.

3 Experimental design

3.1 SSVEP dataset

In this study, the subband parameters in the BsCCA algorithm are calculated using the SSVEP dataset of the finals of the 2020 World Robotics Competition – BCI Brain Control Robot Competition Technical Competition (Dataset I). The dataset contains the EEG data of 10 subjects, and each subject participated in 6 block experiments. The first three blocks were offline experiments, and the last three blocks were online experiments. The stimulation paradigm contains 35 targets, the stimulation frequency range was 3–20 Hz, the interval was 0.5 Hz, and all the initial phases were π / 2. The experimental data was in the block, and each block was the EEG data collected continuously. A single trial lasted 5 s, including 3 s stimulation time and 2 s rest time. In the stimulation process, 35 targets were presented simultaneously, and the motion amplitude of each target changed cosine according to its predetermined frequency. Subjects looked at the prompted target during stimulation to produce the steady-state visual evoked response in their EEG signal. The trigger was recorded at the beginning of the stimulation phase in each trial. The experimental data were 65 channels of EEG data with a sampling frequency of 1000 Hz without other filtering processing.

3.2 Paradigm design

As shown in Fig. 6(A), the paradigm consists of 35 radial contraction-expansion targets with a gray background outside the stimulus targets. The target stimulation frequency was from 3 to 20 Hz, the interval was 0.5 Hz, and the initial phase was π / 2. The static image of the stimulus pattern was composed of 8 concentric rings. Each concentric ring was divided into 24 equal parts, composed of white and black lattice alternately, and the area of the white area was equal to that of the black area. The radius of the outermost circle of the stimulus target was 100 pixels. The radii of the 8 concentric rings ranged from 12.5 pixels to 100 pixels and were evenly distributed according to the difference value of 12.5 pixels. A gray area with a radius of 6.25 pixels was set in the center of the circle, and a black dot with a radius of 1.25 pixels was set in the center of the gray area to focus on the subjects’ attention. Fig. 6(B) shows the distribution of the target stimulation frequency.

Fig. 6

Paradigm interface design and single target flip frequency. (A) Paradigm interface. (B) Flip frequency of each target.

In the contraction-expansion motion paradigm, the display’s high refresh rate is an essential factor in reducing flicker perception [4]. This study sets the screen refresh rate to 240 Hz, consistent with the 2020 World Robotics Competition—BCI Brain Control Robot Competition SSVEP Training Competition.

3.3 Data collection

Neusen W series wireless EEG acquisition system produced by NEURACLE company was used in this study. In the experiment, 64 channels EEG cap arranged according to the international 10–20 system was used. The display model used was Dell AW2518H (1920 × 1080 resolution). The Tsinghua Institutional Review Board approved the experiment, and all subjects signed an informed consent form before the experiment and received a certain reward after the experiment.

Thirteen subjects participated in the experiment, including five women and eight men. All subjects had normal or corrected to normal vision. The subjects sat on a chair and relaxed during the experiment, keeping the distance between the eyes and the screen about 60 cm. After the experiment began, the subjects were required to minimize blinking and physical activity. The experiment was divided into two parts: offline and online. The experiment took the block as the unit for data collection. Four block data were collected for offline and online experiments, and the sampling rate was 1000 Hz. Each block contained 35 trials in the offline experiment, and each target was prompted once in random order. Each trial was stimulated for 4 s, and there was 1 s to find the target between the two trials and 2 min of rest between the two blocks. The offline data of each subject was used as the training data of the subject’s online experiment. Each block of the online experiment also contained 35 trials, and each trial was stimulated for 4 s. There was a 1 s search time between the two trials and a 2 min rest time between the two blocks.

3.4 Parameter selection

In data processing, the algorithm starts to intercept data with a delay of 140 ms from the start of each trial to avoid using the low SNR data at the start of the trial [26]. In the preprocessing stage, the algorithm used a comb filter to suppress the interference of 50 Hz power frequency components [27]. In brain topography, responses evoked by the radial contraction-expansion motion paradigm were concentrated in the occipital region [4]. Therefore, the algorithm used nine channels in the occipital region for calculation (PO3、PO4、PO5、PO6、PO7、PO8、O_z、O1、O2). The all-pass subband covers all target frequencies in the range of [3, 20] Hz, and the weight was set to 1; the high-frequency optimization subband covers the high-frequency components that need to be enhanced, in the range of [13, 20] Hz, and the weight was set to 0.25.

The offline algorithm comparison used the leave-one-out cross-validation method. Each group retained one block as the test sample, and the remaining blocks were used as the training sample. The maximum calculation time of the dynamic window of the eBSCCA algorithm was set to 2 s, the dynamic window threshold value ranged from 1 to 10, the value interval was 1, and there were 10 data points in total.

The online experiment used the four blocks of the offline experiment as training data, and the leave-one-out cross-validation method was used in the training process. One block was reserved for each group as a verification sample, and the rest of the blocks were used as training samples. During the training process, the maximum calculation time of the dynamic window of the eBSCCA algorithm was set to 2 s, the dynamic window threshold value ranged from 1 to 10, and the value interval was 1, with a total of 10 data points. Finally, the dynamic window threshold with the highest average ITR was used as the threshold used by the online algorithm.

3.5 Contrast algorithm

This study used the eCCA and TRCA algorithms as comparison algorithms. The eCCA algorithm is composed of four correlation coefficients [13]. The TRCA algorithm used filter banks to decompose the SSVEP into subband components to efficiently extract the independent information embedded in the harmonic components [12]. In the literature [12], the subband setting method of the TRCA algorithm is consistent with the FBCCA algorithm, and the upper and lower cutoff frequencies of the mth subband are m× 8Hz and 90 Hz, respectively. Among them, the upper cutoff frequency was set to 90 Hz because the EEG signals evoked by the SSVEP paradigm showed high SNR in the frequency band of the upper limit frequency of 90 Hz. The FBCCA algorithm was proposed for the luminance modulation paradigm with the stimulus frequency of 8 to 15.8 Hz, and the lower cutoff frequency corresponds to each harmonic frequency band [21]. According to the stimulation frequency range of the radial contraction-expansion motion paradigm from 3 to 20 Hz, the subband frequency bands of the TRCA algorithm were set as [3, 90] Hz, [20, 90] Hz, [40, 90] Hz, [60, 90] Hz, and [80,90] Hz, the weights were set using Eq. (7):

δ (m) = m^{- 1.25} + 0.25

where m is the subband index, m∊[1,5].

In BCI research, average ITR and accuracy are common performance evaluation indicators of algorithms. The block accuracy rate was obtained by dividing the number of correct predictions for each block by the total number of pre-tests, and the accuracy rate for each subject was the average of the block accuracy rates. The formula for calculating the ITR is [28]:

ITR = \frac{1}{T} [\log_{2} Q + (1 - P) \log_{2} \frac{1 - P}{Q - 1} + P \log_{2} P]

where T represents the decision time, and in this study, the decision time T of each subject was added with 0.5 s of sight shift and focused time on the basis of the time window used by the algorithm [18]; Q represents the number of stimulus targets; and P represents the decision accuracy. The T value in each block is the mean value of the decision time of each trial in the block, and the decision time T of each subject is the average value of the decision time of each block.

4 Results

4.1 Subband and weight parameter calculations

Subbands and weights were obtained according to Dataset I. The all-pass subband was set to [3, 20] Hz, covering all target frequencies and the weight α ₁ was 1. The high-frequency optimization subband was set to [x, 20] Hz, focusing on enhancing high-frequency components, and the weight was α ₂. The lower cutoff frequency x and weighting coefficient α ₂ of the high-frequency optimization subband were calculated using the grid search method. The value range of x was 9–16, and the interval was 1. The algorithm used nine channels in the occipital region for calculation (PO3、PO4、PO5、PO6、PO7、PO8、O_z、O1、O2). The maximum calculation time of the dynamic window was set to 2 s, and the threshold was 20.

Fig. 7 shows the grid search results. The results show that the optimal ITR (77.84 bits/min) can be obtained when the frequency range of the high-frequency optimized subband and weight is set to [13, 20] Hz and 0.25, respectively.

4.2 Offline experiment results

The offline analysis was derived from the offline experimental dataset in this study. Fig. 2(B) shows the SNR of SSVEP responses evoked by different stimulus targets in the offline dataset. From Fig. 2(B), it can be observed that the SNR of the SSVEP component increases with the frequency increase at the stimulation frequency of 3–13 Hz. The SNR of SSVEP components decreased with frequency increase under stimulation at 13–20 Hz. The high-frequency optimized subband calculated according to Dataset I in this study was [13, 20] Hz, which can completely cover the frequency range where the SSVEP SNR decreases.

Fig. 7

Grid search result. The abscissa is the weight of the high-frequency optimized subband (subband2), the ordinate is the lower cutoff frequency of the high-frequency optimized subband (subband2), and the color blocks represent the average ITR value of 10 subjects.

In Fig. 8, each brain topography map corresponds to a stimulation frequency, and the energy of each channel was calculated from the offline data corresponding to the stimulation frequency. This study normalized each subject using its maximum channel energy as a standard to eliminate the effect of energy differences between subjects. It averaged the normalized values of 13 subjects. Fig. 8 shows that, except for the response of the lowest stimulation frequency of 3 Hz, which was relatively scattered, the rest of the target responses were concentrated in the occipital area.

Fig. 8

Normalized EEG energy distribution topographic map.

As shown in Fig. 9, under the radial contraction-expansion motion paradigm, the eCCA algorithm can obtain the best recognition effect by using the reference signal with a harmonic number of 1. This optimization parameter will be used for algorithm comparison.

Fig. 9

The influence of harmonic number of reference signal on the performance of eCCA algorithm. (A) Average recognition accuracy. The abscissa is the time window length, and the ordinate is the average recognition accuracy. (B) Average ITR. The abscissa is the time window length, and the ordinate is the average ITR.

As shown in Fig. 10, under the radial contraction-expansion motion paradigm, the TRCA algorithm can get the best recognition effect by setting the number of subbands to 2. This optimization parameter will be used for algorithm comparison.

Fig. 11 shows the comparison results of average accuracy and ITR under different algorithms. As shown in Fig. 11(A), the recognition accuracy of each algorithm increases with the increase in decision time. Fig. 11(B) shows that the TRCA algorithm reached the highest average ITR of 75.67 ± 49.85 bits/min in the 1.5 s time window. The eCCA algorithm reached the highest average ITR of 107.66 ± 61.79 bits/min in the 1.4 s time window. Furthermore, the eBSCCA algorithm reached the highest average ITR of 151.34 ± 42.45 bits/min in the average decision time of 1.59 s.

Each algorithm’s average ITR and accuracy under the highest average ITR were taken for a paired t-test with Bonferroni correction. Table 1 shows the results. The test results show that the p-value of average ITR and accuracy between the eBSCCA and TRCA algorithms was less than 0.01. The p-value of average ITR and accuracy between the eBSCCA algorithm and the eCCA algorithm was less than 0.01. The experimental results show that the proposed eBSCCA algorithm significantly improves the average ITR and accuracy (p < 0.01).

Fig. 10

The effect of the number of subbands on the performance of the TRCA algorithm. (A) Average recognition accuracy. The abscissa is the time window length, and the ordinate is the average recognition accuracy. (B) Average ITR. The abscissa is the time window length, and the ordinate is the average ITR.

Fig. 11

Comparison of average accuracy and ITR under different algorithms. (A) Average accuracy comparison. The errorbar represents the standard deviation of the average accuracy of each subject. (B) Average ITR comparison. The errorbar represents the standard deviation of the average ITR of each subject. The TRCA and teCCA algorithms contain 16 data points corresponding to the decision window length in the range of 1–2.5 s, and the interval is 0.1. The eBSCCA algorithm contains 10 data points, which correspond to dynamic window thresholds in the range of 1–10, with an interval of 0.1. Since the decision window length when the dynamic window threshold is similar to that when the dynamic window threshold is 2, the first and second data points of the eBSCCA algorithm in the figure almost coincide. represents TRCA algorithm, represents eCCA algorithm, represents eBSCCA algorithm. The hollow point corresponds to the highest point of ITR of the algorithm.

Table 1

Paired t-test (with Bonferroni correction) for the mean ITR and mean accuracy of the optimal performance of each algorithm.

	(TRCA, eBSCCA)	(eCCA, eBSCCA)	(eCCA, TRCA)
p-value (ITR)	1.92×10⁻⁸	0.00071	0.0004
p-value (accuracy)	0.00021	0.0091	0.00068

The offline experiment included 1820 trials. As shown in Fig. 12(A), in the eBSCCA algorithm, the calculation workload of the computer increases with the increase in the calculation window length of the trial, and the average calculation workload under the maximum calculation window length was less than 40 ms. In Fig. 12(B), the calculation window length distribution of the trials changes with the subjects’ signal quality. If the signal quality of the subject is good, the length of the calculation window is small. If the signal quality of the subject is poor, the calculation window is large. As shown in Fig. 12(C), the eBSCCA algorithm uses the dynamic window method and the computer calculation workload changes with the signal quality. The better the quality of the tested signal, the less the computer calculation workload. In Figs. 12(D) and (E), the TRCA and ECCA algorithms used fixed windows for the calculation, so the computer calculation workload of each subject was less different.

4.3 Online experiment results

Table 2 shows the specific values of the online experiment results of each subject.

As shown in Table. 2, the average discrimination accuracy of the subjects reached 88.68%, and the discrimination accuracy of all subjects was greater than 80%. The average ITR of the subjects reached 158.77 bits/min. Except for subject S3, the ITRs of the other subjects were all greater than 100 bits/min. Among them, the ITR of 4 subjects reached 200 bits/min, and the average ITR of subject S8 was the highest, reaching 237.72 bits/min. Except that the average decision time of subject S3 was longer than 2 s, which was 2.22 s. The decision times of the other subjects were all within the range of 1–2 s.

Fig. 12

Different algorithms to calculate time-consuming statistics. (A) Computer calculation workload of different calculation window lengths in the eBSCCA algorithm. The abscissa represents the calculation window length of the trial (including 0.5 s of sight transfer time), and the ordinate represents the average calculation workload of the computer corresponding to each calculation window length range. (B) Statistics on the length of calculated data used in trials in the eBSCCA algorithm. The abscissa represents the calculation window length of the trial (including 0.5 s of sight transfer time), and the ordinate represents the number of trials corresponding to each calculation window length range. (C) The time-consuming statistics of the eBSCCA algorithm with optimal parameters under the dynamic stop window length condition. The abscissa represents the computer calculation workload of trials, and the ordinate represents the number of trials corresponding to the computer calculation workload range. (D) Time-consuming statistics of eCCA algorithm with optimal parameters. (E) Time-consuming statistics of TRCA algorithm with optimal parameters.

Table 2

Online experiment results of each subject (including 0.5 s eye shift time).

Subject	Age (years), gender	Average ITR (bits/min)	Average accuracy (%)	Average time (s)
S1	24, male	142.97	88.57	1.70
S2	26, male	202.80	92.86	1.30
S3	26, male	92.18	80.00	2.22
S4	23, female	158.19	92.14	1.65
S5	25, female	161.17	91.43	1.59
S6	26, female	178.46	81.43	1.18
S7	19, male	201.01	95.71	1.39
S8	23, male	237.72	95.71	1.18
S9	25, male	123.28	82.86	1.76
S10	22, male	201.37	97.86	1.46
S11	22, male	111.20	81.43	1.89
S12	24, female	142.18	90.71	1.78
S13	21,female	111.45	82.14	1.91
Average	23.46 ± 2.15	158.77 ± 43.67	88.68 ± 6.33	1.6 ± 0.31

5 Discussion

In the current work we propose an extended BsCCA detection algorithm for the radial contraction-expansion motion paradigm. This algorithm improves the CCA algorithm according to the characteristics of the high fundamental frequency component, less doubling frequency component, and the solid low-frequency response of EEG signals evoked by the radial contraction-expansion motion paradigm. The BsCCA algorithm was used to improve the recognition accuracy of high-frequency stimuli. The dynamic window method was used to minimize the required data length to achieve higher ITR and improve the performance of BCI. The offline experimental results showed that the proposed algorithm significantly outperforms the eCCA and TRCA algorithms in the signals evoked by the radial contraction-expansion motion paradigm (p < 0.01). In the online experiment, the average recognition accuracy of 13 subjects reached 88.68% ± 6.33%, and the average ITR reached 158.77 ± 43.67 bits/min, which proved that the algorithm has a good recognition effect in the signal evoked by the radial contraction-expansion motion paradigm. Moreover, in the 2020 World Robotics Competition – BCI Brain Control Robot Competition, the algorithm won first place in the final of the SSVEP training competition, which further proved its usability.

As shown in Fig. 13, the recognition performance of the three algorithms was improved with the increase in the amount of training data. Among them, the TRCA algorithm was most affected by the amount of training data. In the case of less training data, the advantages of the proposed ebscca algorithm were more obvious. Research [12] has shown that both the eCCA and TRCA algorithms are greatly affected by training data. When the number of training blocks was 3, the prediction accuracy of the algorithm was less than 60%. In the offline algorithm comparison in this paper, the number of training blocks was 3. Therefore, the lack of training data may be one of the reasons for the poor performance of the eCCA and TRCA algorithms.

For subjects S3, S9, S11, and S13, the recognition accuracy of 30 target frequencies in the range of 5.5–20 Hz has been improved to varying degrees without 3–5 Hz target frequency. In the offline data set, the time window length of the CCA algorithm was set to 3 s. After removing the target frequency of 3–5 Hz, the average recognition accuracy of subject S3, S9, S11, and S13 increased from 61.67%–68.33%, 90.83%–99.17%, 60%–70.83%, and 55.83%–84.17%, respectively. The recognition accuracy was dramatically improved, especially in subject S13, by nearly 30%. Therefore, to enhance the universality of this paradigm, the stimulation frequency range of the paradigm should be improved. Although, the flip frequency of 3 Hz is too low, and some subjects have a poor recognition effect.

6 Conclusion

The radial contraction-expansion motion paradigm is different from the traditional luminance modulation paradigm. This paradigm has a weak sense of flicker. Additionally, the subjects are not prone to visual fatigue during use, reducing the risk of photosensitive epilepsy. Therefore, this paradigm is an excellent choice to replace the luminance modulation paradigm. However, the current research on this paradigm mainly focuses on the comparative study of the EEG performance evoked by the paradigm and the practical study of the paradigm. There is a lack of research on applicable algorithms. Thus, the extended BsCCA algorithm proposed in this paper for the radial contraction-expansion motion paradigm can effectively improve EEG signal recognition performance and provide a good algorithm choice for the practical application of this paradigm in the future.

Fig. 13

The influence of the number of training data blocks on the performance of each algorithm. (A) The abscissa represents the number of training data blocks, and the ordinate represents the average recognition accuracy. (B) The abscissa represents the number of training data blocks, and the ordinate represents the average ITR.

The radial contraction-expansion motion paradigm has been developed for spellers at the application level, and it works well. In addition, this paradigm has also been used to develop attention-driven games, and researchers believe that this paradigm can be used to train participants’ attention, which can be used in special education schools or centers dealing with attention disorders in the future [29].

Footnotes

Consent

The experiment has been approved by the Tsinghua Institutional Review Board, and all subjects signed an informed consent form before the experiment.

Conflict of interests

All contributing authors report no conflict of interests in this work.

Funding

This work is granted by National Natural Science Foundation of China (Grant Nos. 62006024, 62071057), the Fundamental Research Funds for the Central Universities (BUPT Project No. 2019XD17), Aeronautical Science Foundation of China (NO. 2019ZG073001).

Acknowledgements

All the authors thank the subjects who worked for many hours on this project, and also thank the director of the 2020 World Robotics Competition –BCI Brain Control Robot Competition Technical Competition for providing the dataset.

Authors’ contribution

Yuxue Zhao completed the experimental design, performed the data analyses and wrote the manuscript. Hongxin Zhang performed analysis on the algorithm and directed experiments. Chen Yang contributed to the conception of the study and designed experimental systems. Yuanzhen Wang, Chenxu Li and Ruilin Xu provided technical support for the construction of the experimental platform.

References

Wang

Gao

Hong

, et al. Brain-computer interfaces based on visual evoked potentials. IEEE Eng Med Biol Mag 2008, 27(5): 64–71.

Yao

Tang

, et al. Amplitude modulation of steady-state visual evoked potentials by event-related potentials in a working memory task. J Biol Phys 2010, 36(3): 261–271.

Zhang

Xie

, et al. Brain response to luminance-based and motion-based stimulation using inter-modulation frequencies. PLoS One 2017, 12(11):e0188073.

Han

Xie

, et al. Highly interactive brain-computer interface based on flicker-free steady-state motion visual evoked potential. Sci Rep 2018, 8(1): 5835.

Yan

Chen

, et al. Steady-state motion visual evoked potential (SSMVEP) enhancement method based on time-frequency image fusion. Comput Intell Neurosci 2019, 2019: 9439407.

Chai

Zhang

Guan

, et al. A radial zoom motion-based paradigm for steady state motion visual evoked potentials. Front Hum Neurosci 2019, 13: 127.

Galloway

. Human brain electrophysiology: evoked potentials and evoked magnetic fields in science and medicine. Br J Ophthalmol 1990, 74(4): 255

Zhang

Hou

, et al. Enhancing detection of SSMVEP induced by action observation stimuli based on task-related component analysis. Sensors (Basel) 2021, 21(16): 5269.

Bin

Gao

Yan

, et al. An online multi-channel SSVEP-based brain-computer interface using a canonical correlation analysis method. J Neural Eng 2009, 6(4): 046002.

10.

Lin

Zao

, et al. SNR analysis of high-frequency steady-state visual evoked potentials from the foveal and extrafoveal regions of human retina. Annu Int Conf IEEE Eng Med Biol Soc 2012, 2012: 1810–1814.

11.

Dang

, et al. MHLCNN: Multi-harmonic linkage CNN model for SSVEP and SSMVEP signal classification. IEEE Trans Circuits Syst II-Express Briefs 2022, 69(1): 244–248.

12.

Nakanishi

Wang

Chen

, et al. Enhancing detection of SSVEPs for a high-speed brain speller using task-related component analysis. IEEE Trans Biomed Eng 2018, 65(1): 104–112.

13.

Nakanishi

Wang

, et al. A high-speed brain speller using steady-state visual evoked potentials. Int J Neural Syst 2014, 24(6): 1450019.

14.

Parameshwaran

Thiagarajan

. Characterizing peaks in the EEG power spectrum. Biomed Phys Eng Express 2019, 5(4): 045023.

15.

Islam

Tanaka

Nakanishi

, et al. Frequency recognition of steady-state visually evoked potentials using binary subband canonical correlation analysis with reduced dimension of reference signals. In 2016 IEEE International Conference on Acoustics, Speech and Signal Processing. Shanghai, China, 2016, pp 769–773.

16.

Yan

, et al. Steady-state motion visual evoked potential (SSMVEP) based on equal luminance colored enhancement. PLoS One 2017, 12(1): e0169642.

17.

Yan

Xie

, et al. Four novel motion paradigms based on steady-state motion visual evoked potential. IEEE Trans Biomed Eng 2018, 65(8): 1696–1704.

18.

Yang

Han

Wang

, et al. A dynamic window recognition algorithm for SSVEP-based brain-computer interfaces using a spatio-temporal equalizer. Int J Neural Syst 2018, 28(10): 1850028.

19.

Rabiul Islam

Khademul Islam Molla

Nakanishi

, et al. Unsupervised frequency-recognition method of SSVEPs using a filter bank implementation of binary subband CCA. J Neural Eng 2017, 14(2): 026007.

20.

Correlation analysis of visual verbs’ subcategorization based on Pearson’s correlation coefficient. In 2010 International Conference on Machine Learning and Cybernetics. Qingdao, China, 2010, pp 2042–2046.

21.

Chen

Wang

Gao

, et al. Filter bank canonical correlation analysis for implementing a high-speed SSVEP-based brain-computer interface. J Neural Eng 2015, 12(4): 046008.

22.

Chen

Yang

Chen

, et al. A novel training-free recognition method for SSVEP-based BCIs using dynamic window strategy. J Neural Eng 2021, 18(3), DOI 10.1088/1741-2552/ab914e.

23.

Chen

Zhao

, et al. Adaptive asynchronous control system of robotic arm based on augmented reality-assisted brain-computer interface. J Neural Eng 2021, 18(6), DOI 10.1088/1741-2552/ ac3044.

24.

Joanes

Gill

. Comparing measures of sample skewness and kurtosis. J Royal Stat Soc Ser D Stat 1998, 47(1): 183–189.

25.

Balanda

Macgillivray

. Kurtosis: A critical review. Am Stat 1988. 42(2): 111–119.

26.

Chen

Wang

Nakanishi

, et al. High-speed spelling with a noninvasive brain-computer interface. Proc Natl Acad Sci USA 2015, 112(44): E6058–E6067.

27.

Wolpaw

Birbaumer

McFarland

, et al. Brain-computer interfaces for communication and control. Clin Neurophysiol 2002, 113(6): 767–791.

28.

Obermaier

Neuper

Guger

, et al. Information transfer rate in a five-classes brain-computer interface. IEEE Trans Neural Syst Rehabil Eng 2001, 9(3): 283–288.

29.

Perez-Valero

Lopez-Gordo

Vaquero-Blasco

. An attention-driven videogame based on steady-state motion visual evoked potentials. Expert Syst 2021, 38(4): e12682.