The Nyquist frequency \(f_N = f_s/2 = 500\) Hz ensures adequate representation of all relevant neurophysiological frequencies. The anti-aliasing filter was implemented as a zero-phase Butterworth filter of order 8:

where \(\omega _c = 2\pi \cdot 450\) rad/s is the cutoff frequency and \(n = 8\) is the filter order.

Motor imagery paradigm design

The experimental paradigm consisted of four distinct motor imagery classes: left hand (LH), right hand (RH), both feet (BF), and tongue (T) movements. The trial structure followed a carefully designed temporal sequence to optimize neural response capture:

where: – \(t_{\text {fix}} \in [0, 2]\) s: fixation cross presentation – \(t_{\text {cue}} \in [2, 3]\) s: visual cue presentation – \(t_{\text {MI}} \in [3, 7]\) s: motor imagery execution – \(t_{\text {rest}} \in [7, 9]\) s: inter-trial interval

The probability of each class appearing was balanced according to:

Signal preprocessing pipeline

Artifact detection and removal

The preprocessing pipeline begins with automated artifact detection using a multivariate approach⁴⁸. For each channel c and time window w, we compute the standardized amplitude:

where \(\mu _{c,w}\) and \(\sigma _{c,w}\) are the local mean and standard deviation. Channels are marked as artifactual if:

where \(\tau _z = 5\), \(\tau _k = 10\), and \(\tau _p\) are empirically determined thresholds, and Kurt denotes kurtosis.

Attention network training

Independent component analysis for artifact removal

We employ extended Infomax ICA to decompose the multichannel EEG into statistically independent components:

where \({X} \in 97.2477{R}^{C \times T}\) is the observed EEG data, \({A} \in 97.2477{R}^{C \times C}\) is the mixing matrix, and \({S} \in 97.2477{R}^{C \times T}\) contains the independent components. The unmixing matrix \({W} = {A}^{-1}\) is obtained by maximizing the entropy:

where \({y} = {W}{x}\) and \(p_i\) is the probability density function of the i-th component.

The weight update rule follows the natural gradient approach:

where \(\eta\) is the learning rate and \({\phi }({y})\) is the score function vector. Artifactual components are identified using multiple criteria:

where MI denotes mutual information with EOG channels, and \(P_{s_i}(f)\) is the power spectral density of component \(s_i\).

Spectral filtering and band-power extraction

For motor imagery analysis, we focus on specific frequency bands known to exhibit task-related modulations. The filter bank approach employs multiple bandpass filters⁴⁹:

where the coefficients are designed for the following frequency bands: – \(\theta\): [4, 8] Hz – \(\alpha\): [8, 13] Hz – \(\beta _{\text {low}}\): [13, 20] Hz – \(\beta _{\text {high}}\): [20, 30] Hz – \(\gamma _{\text {low}}\): [30, 50] Hz

The instantaneous band power is computed using the Hilbert transform:

where \(\mathscr {H}\) denotes the Hilbert transform operator.

Spatial filtering using optimized common spatial patterns

The traditional CSP algorithm is enhanced through regularization and multi-class extensions. For the two-class case, the optimization problem becomes:

This is solved via the generalized eigenvalue problem:

For improved robustness, we employ Tikhonov regularization:

where \(\gamma \in [0,1]\) is the regularization parameter determined through cross-validation.

The multi-class extension uses the one-versus-rest (OVR) approach:

Feature extraction and dimensionality reduction

After spatial filtering, we extract logarithmic band power features:

where \(T_w\) is the time window length and \({w}_i\) is the i-th spatial filter. Principal Component Analysis (PCA) is applied for dimensionality reduction:

where \({U}_d\) contains the top d eigenvectors of the feature covariance matrix:

The number of components d is selected to retain 95

Deep learning architectures

Convolutional neural network for spatial feature learning

The convolutional architecture is specifically designed for EEG spatial pattern extraction. The first layer performs spatial convolution across electrode locations:

where C is the number of channels, \(T_k\) is the temporal kernel size, and \(\sigma\) is the activation function.

The spatial attention mechanism is incorporated as:

The attended feature representation becomes:

Long short-term memory networks for temporal dynamics

The LSTM architecture captures the temporal evolution of motor imagery patterns. The complete LSTM dynamics are governed by:

where \(\odot\) denotes element-wise multiplication, and \({V}_f, {V}_i, {V}_o\) are the peephole connection weights.

Schematic diagram of CNN-LSTM model structure.

The gradient flow through LSTM is carefully managed to prevent vanishing gradients:

This product of forget gates allows long-term dependencies to be preserved when \({f}_t \approx 1\).

Bidirectional processing for complete temporal context

The bidirectional LSTM processes the sequence in both forward and backward directions:

The combined representation is:

This bidirectional processing captures both causal and anti-causal dependencies in the EEG signal.

Hierarchical CNN-LSTM architecture

The hybrid CNN-LSTM model leverages the complementary strengths of convolutional and recurrent architectures in Fig. 2. The hierarchical processing pipeline can be formalized as:

where \(\mathscr {S}\) represents spatial filtering, \(\mathscr {C}_{\text {CNN}}\) denotes convolutional feature extraction, \(\mathscr {P}\) is pooling, \(\mathscr {R}\) is reshaping, and \(\mathscr {L}_{\text {LSTM}}\) represents temporal modeling.

The convolutional block consists of multiple layers with progressively increasing receptive fields:

where \(*\) denotes the convolution operation, \(K_l\) is the number of feature maps in layer l, and \(\phi\) is the activation function.

We employ depthwise separable convolutions to reduce computational complexity:

The depthwise convolution applies a single filter per input channel:

Followed by pointwise convolution to combine channels:

Multi-head attention mechanism

The attention mechanism enables the model to focus on task-relevant spatiotemporal patterns. We implement multi-head attention with H parallel attention heads:

where each head is computed as:

The scaled dot-product attention is:

where \(d_k\) is the dimension of the key vectors, and the scaling factor \(1/\sqrt{d_k}\) prevents gradient vanishing in the softmax function. For temporal attention over LSTM hidden states, we compute:

The attention weights are normalized:

The context vector aggregates the attended features:

Self-attention for global dependencies

To capture long-range dependencies in EEG signals, we incorporate self-attention layers:

The position-aware self-attention includes positional encodings:

The input with positional encoding becomes:

Real-time motor imagery classification

Layer normalization placement: We employ a pre-normalization architecture following the recent advances in transformer design. The complete attention block implements the following sequence:

The LayerNorm operation is applied with learnable parameters \(\gamma\) and \(\beta\):

where \(\mu\) and \(\sigma ^2\) are computed across the feature dimension, and \(\epsilon = 1e^{-6}\) for numerical stability.

Attention dropout configuration: We implement dropout at three specific locations within the attention mechanism:

1. 1.

   Attention weight dropout (\(p_{attn} = 0.3\)): Applied to the attention probability matrix after softmax normalization:

   $$\begin{aligned} \text {Attention}(Q,K,V) = \text {Dropout}\left( \text {softmax}\left( \frac{QK^T}{\sqrt{d_k}}\right) , 0.3\right) V \end{aligned}$$

   (69)

2. 2.

   Output projection dropout (\(p_{proj} = 0.2\)): Applied to the concatenated multi-head output before the final linear projection:

   $$\begin{aligned} \text {MultiHead}(Q,K,V) = \text {Dropout}(\text {Concat}(\text {head}_1, …, \text {head}_H), 0.2) W^O \end{aligned}$$

   (70)

3. 3.

   Residual connection dropout (\(p_{res} = 0.1\)): Applied to the attention output before the residual addition:

   $$\begin{aligned} H^{(l)} = H^{(l-1)} + \text {Dropout}(\text {Attn}^{(l)}, 0.1) \end{aligned}$$

   (71)

L2 regularization on attention weights: We apply L2 regularization specifically to the attention projection matrices with coefficient \(\lambda _{attn} = 0.001\):

The total loss function becomes:

Attention head initialization: Each attention head is initialized using Xavier uniform initialization:

where \(d_{model} = 256\) and \(d_k = 32\).

Gradient scaling for attention: To stabilize training, we apply gradient scaling to attention parameters:

where \(\alpha _{scale} = 0.5\) prevents attention gradients from dominating the optimization.

Training methodology and optimization

Loss function design

For multi-class motor imagery classification, we employ a weighted cross-entropy loss to handle class imbalance:

where \(w_k = \frac{N}{K \cdot N_k}\) is the class weight, \(N_k\) is the number of samples in class k.

To improve feature discrimination, we add a center loss term:

where \({f}_i\) is the feature representation and \({c}_{y_i}\) is the center of class \(y_i\).

The total loss combines multiple objectives:

where \(\mathscr {L}_{\text {reg}}\) is the regularization term:

combining L2 (Frobenius norm) and L1 regularization for weight sparsity.

Schematic diagram of RF model structure.

Optimization algorithm

We employ the Adam optimizer with modifications for improved convergence:

where \({g}_t = \nabla _{\varvec{\theta }} \mathscr {L}(\varvec{\theta }_{t-1})\) is the gradient at time t.

The learning rate follows a cosine annealing schedule with warm restarts:

where \(T_{\text {cur}}\) is the number of epochs since the last restart.

Gradient clipping and normalization

To prevent gradient explosion in deep networks, we apply gradient clipping:

Layer normalization is applied to stabilize training:

where \(\mu\) and \(\sigma ^2\) are the mean and variance computed across the feature dimension.

Data augmentation strategies

To improve model generalization, we implement several EEG-specific data augmentation techniques:

1. 1.

   **Temporal Shifting**: Random temporal shifts within physiologically plausible ranges in Fig. 3:

   $$\begin{aligned} \tilde{{x}}
   (88)

2. 2.

   **Amplitude Scaling**: Channel-wise amplitude modulation:

   $$\begin{aligned} \tilde{x}_c
   (89)

3. 3.

   **Gaussian Noise Addition**:

   $$\begin{aligned} \tilde{{x}}
   (90)

   where \(\sigma _n = 0.01 \cdot \text {std}({x})\).

4. 4.

   **Mixup Regularization**: Linear interpolation between training samples:

   $$\begin{aligned} \tilde{{x}} = \lambda {x}_i + (1-\lambda ) {x}_j, \quad \tilde{y} = \lambda y_i + (1-\lambda ) y_j \end{aligned}$$

   (91)

   where \(\lambda \sim \text {Beta}(\alpha , \alpha )\) with \(\alpha = 0.2\).

Evaluation metrics and statistical analysis

Performance metrics

Beyond accuracy, we employ multiple metrics to comprehensively evaluate model performance:

1. 1.

   **Cohen’s Kappa**: Accounts for chance agreement:

   $$\begin{aligned} \kappa = \frac{p_o – p_e}{1 – p_e} \end{aligned}$$

   (92)

   where \(p_o\) is the observed agreement and \(p_e\) is the expected agreement by chance.

2. 2.

   **Matthews Correlation Coefficient**: Balanced measure for multi-class:

   $$\begin{aligned} \text {MCC} = \frac{\sum _{k,l,m} C_{kk}C_{lm} – C_{kl}C_{mk}}{\sqrt{\sum _k(\sum _l C_{kl})(\sum _{l’} C_{l’k})}\sqrt{\sum _k(\sum _l C_{lk})(\sum _{l’} C_{kl’})}} \end{aligned}$$

   (93)

   where C is the confusion matrix.

3. 3.

   **Area Under the ROC Curve (AUC)**: For multi-class, we use one-vs-rest:

   $$\begin{aligned} \text {AUC}_{\text {macro}} = \frac{1}{K}\sum _{k=1}^K \text {AUC}_k \end{aligned}$$

   (94)

Statistical significance testing

To assess the statistical significance of performance differences, we employ:

1. 1.

   **Wilcoxon Signed-Rank Test**: For paired comparisons across subjects:

   $$\begin{aligned} W = \sum _{i=1}^n \text {sgn}(x_{2,i} – x_{1,i}) \cdot R_i \end{aligned}$$

   (95)

   where \(R_i\) is the rank of \(|x_{2,i} – x_{1,i}|\).

2. 2.

   **Friedman Test**: For comparing multiple models:

   $$\begin{aligned} \chi _F^2 = \frac{12n}{k(k+1)}\sum _{j=1}^k R_j^2 – 3n(k+1) \end{aligned}$$

   (96)

   where n is the number of subjects, k is the number of models, and \(R_j\) is the sum of ranks for model j.

3. 3.

   **Post-hoc Nemenyi Test**: For pairwise comparisons after Friedman test:

   $$\begin{aligned} CD = q_{\alpha } \sqrt{\frac{k(k+1)}{6n}} \end{aligned}$$

   (97)

   where \(q_{\alpha }\) is the critical value from the Studentized range distribution.

Cross-validation strategy

We implement a nested cross-validation scheme to ensure unbiased performance estimation:

where \(\varvec{\theta }_i^*\) is obtained from inner cross-validation:

Nested cross-validation protocol

To ensure unbiased performance estimation and robust hyperparameter optimization, we implemented a rigorous nested cross-validation strategy with explicit leave-two-subjects-out (L2SO) outer loop validation. The complete protocol consists of the following components:

Outer loop (model evaluation): We employed a leave-two-subjects-out cross-validation scheme, systematically holding out 2 participants for testing while using the remaining 13 participants for training and validation. This resulted in \(\left( {\begin{array}{c}15\\ 2\end{array}}\right) = 105\) outer folds, ensuring that every possible pair of subjects served as the test set exactly once. The outer loop provides an unbiased estimate of model generalization to unseen subjects, which is critical for BCI applications where cross-subject variability is substantial.

Inner loop (hyperparameter optimization): For each outer fold, we performed hyperparameter optimization using the 13 training subjects through stratified 5-fold cross-validation. The inner training set was split using an 80%-20% ratio, where 80% (approximately 10.4 subjects worth of data) was used for model training and 20% (approximately 2.6 subjects worth of data) was reserved for validation during hyperparameter search.

Parameter grid specification: The hyperparameter search space encompassed the following parameters:

• Learning rates: \(\alpha \in \{0.0001, 0.0005, 0.001, 0.005, 0.01\}\)

• Dropout rates: \(p_1 \in \{0.15, 0.2, 0.25, 0.3\}\), \(p_2 \in \{0.4, 0.5, 0.6\}\)

• LSTM units: \(N_{lstm} \in \{64, 96, 128, 192, 256\}\)

• Attention heads: \(H \in \{4, 6, 8, 12, 16\}\)

• Regularization weights: \(\lambda _1 \in \{0.001, 0.003, 0.005\}\), \(\lambda _2 \in \{0.0001, 0.0005, 0.001\}\)

This resulted in \(5 \times 4 \times 3 \times 5 \times 5 \times 3 \times 3 = 13,500\) hyperparameter combinations evaluated through Bayesian optimization using Gaussian Process surrogate models.

Early stopping criteria: To prevent overfitting during hyperparameter search, we implemented early stopping with the following criteria: (1) patience of 20 epochs with no improvement in validation loss, (2) minimum delta of 0.001 for considering an improvement, and (3) restoration of best weights upon early termination. The maximum training duration was capped at 200 epochs per configuration.

Performance aggregation: The final performance estimate was computed as:

where \(E_i(\theta _i^*)\) represents the test performance of the optimal hyperparameters \(\theta _i^*\) found for the i-th outer fold. The reported standard deviations reflect the variability across the 105 outer folds, providing a robust estimate of model stability across different subject combinations.

Statistical validation: To ensure the reliability of our nested CV results, we computed 95% confidence intervals using the t-distribution with 104 degrees of freedom (105 folds – 1). Additionally, we performed the Shapiro-Wilk test to verify the normality of performance distributions across folds, confirming that parametric statistical tests were appropriate for our analysis.

Advanced signal processing techniques

Riemannian geometry-based features

EEG covariance matrices lie on the manifold of symmetric positive definite (SPD) matrices. We compute the Riemannian mean of trial covariance matrices:

where \(d_R\) is the Riemannian distance:

The tangent space mapping at the Riemannian mean provides vectorized features:

where \(\text {vec}_u\) is the upper triangular vectorization and \(\log _{\bar{{C}}}\) is the Riemannian logarithm.

Wavelet packet decomposition

For multi-resolution analysis, we employ wavelet packet decomposition:

where \(p \in \{0,1\}\) indicates low-pass (\(h_0\)) or high-pass (\(h_1\)) filtering, and j is the decomposition level.

The best basis selection uses the Shannon entropy criterion:

The optimal wavelet packet tree \(\mathscr {T}^*\) minimizes:

Phase-locking value and connectivity analysis

To assess functional connectivity during motor imagery, we compute the phase-locking value (PLV):

where \(\phi _i^n(f)\) is the instantaneous phase of channel i at frequency f for trial n.

The weighted phase lag index (wPLI) provides a more robust connectivity measure:

where \({S}_{ij}\) is the cross-spectral density between channels i and j.

Deep learning implementation details

Network architecture specifications

The complete CNN-LSTM-Attention architecture consists of the following layers: Spatial convolution block

where \(F_1 = 8\) initial filters, \(K_t = 64\) temporal kernel size, \(D = 2\) depth multiplier, and \(P_t = 8\) pooling size.

Separable convolution block

where \(F_2 = 16\), \(K_t’ = 16\), \(p_1 = 0.25\), and \(P_t’ = 8\).

LSTM block

where \(N_{\text {lstm}} = 128\) LSTM units.

Attention block

*Classification head

where \(N_{\text {dense}} = 64\) and \(p_2 = 0.5\).

Initialization strategies

Proper weight initialization is crucial for deep network training. For convolutional layers, we use He initialization:

For LSTM weights, we employ orthogonal initialization:

Forget gate biases are initialized to 1.0 to encourage information flow:

Ensemble methods and model averaging

To improve robustness, we employ ensemble techniques:

Snapshot ensembling: Save models at different training epochs with cosine annealing:

Weighted averaging: Based on validation performance:

where \(\beta\) is a temperature parameter.

Spatiotemporal attention fusion mechanism

Our spatiotemporal attention employs a parallel dual-stream architecture that simultaneously computes spatial and temporal attention components, followed by adaptive fusion. Given bidirectional LSTM outputs \(H \in 97.2477{R}^{T_{seq} \times d_{hidden}}\), we reshape to \(H_{spatial} \in 97.2477{R}^{T_{seq} \times C \times d_{electrode}}\) to preserve spatial electrode structure, where \(d_{electrode} = d_{hidden}/C\).

The spatial attention operates across electrode channels at each temporal position:

Parallel temporal attention operates across time steps:

The adaptive fusion mechanism dynamically weights spatial and temporal information:

where \(c_{spatial,avg} = \frac{1}{T_{seq}} \sum _{t=1}^{T_{seq}} c_{spatial}^{
(129)