1 Introduction↩︎

Brain-Computer Interfaces (BCIs) hold transformative potential across various domains, including the diagnosis of neurodegenerative diseases [1], advanced brain function mapping, robotic control, and the development of immersive gaming devices [2]. Among the various BCI paradigms, motor imagery (MI) BCIs decode intentionally modulated neural activity resulting from conscious cognitive effort, making them especially promising for developing assistive communication systems and neurorehabilitation protocols for patients with severe motor impairments [2]. Furthermore, BCI implementations leverage various techniques for recording brain activity, including electroencephalography (EEG), electrocorticography, functional magnetic resonance imaging, and functional near-infrared spectroscopy [3]. Among these modalities, BCIs based on EEG are particularly promising for real-world applications because of their high temporal resolution and superior portability.

In this study, we focus on brain activity occurring during motor imagery (MI) and tackle a classification task for imagined actions using EEG signals. Figure 1 presents a representative example of the task. In this case, the input comprises EEG signals recorded while the subject imagines an elbow extension. As its response, the model outputs predicted probabilities for each corresponding action. Despite intensive research, the accurate decoding of MI EEG signals is still non-trivial. For instance, a binary classification task using MI EEG signals from the OpenBMI [4] dataset, a representative model (e.g., [5]), demonstrates an error rate of approximately 22%.

Recent approaches have used Transformers to capture subject-agnostic dependencies in EEG signals (e.g., EEG Conformer [6] and Medformer [7]). However, to address the exponential growth in computational cost associated with longer time-series, these methods patchify and compress the input EEG signals along the temporal dimension, which may result in fine-grained temporal dependencies being lost. Moreover, while these approaches offer interpretability in spatio-temporal domains, they do not directly provide insights into the frequency domain, leaving the contributions of neurophysiologically significant frequency bands (e.g., mu band [8]) unclear.

Therefore, we propose Cortical-SSM, an extension of Deep SSM that explicitly captures integrated EEG dependencies across temporal, spatial, and frequency domains. This design enables the modeling of temporal dependencies in EEG signals without compressing them. Moreover, the proposed method provides direct visual explanations across the temporal, spatial, and frequency domains, in contrast to prior methods that operate on a subset of domains.

In Cortical-SSM, we extend a deep state space model (Deep SSM [9])—an architecture offering superior computational efficiency to Transformers for long sequence modeling—to effectively capture multi-scale temporal dependencies in EEG signals. Specifically, we introduce the Frequency-SSM module, which extracts spatio-temporal dependencies for each frequency component, and the Channel-SSM module, which captures temporal-frequency dependencies for each electrode. Furthermore, we introduce a Wavelet-Convolution module integrating both deterministic and adaptive frequency features for feature extraction from the frequency domain. The Frequency-SSM and Channel-SSM modules model temporal dependencies in EEG signals without requiring the patchification employed in prior studies, thereby facilitating the capture of fine-grained temporal variations. Moreover, by incorporating the Wavelet-Convolution, Cortical-SSM extracts frequency-analyzable features while simultaneously learning their representations. Frequency-SSM and Channel-SSM modules explicitly model frequency-wise and electrode-wise features derived from the Wavelet-Convolution module, enabling direct visual explanations across temporal, spatial, and frequency domains.

Our main contributions lie in the following aspects:

• We propose Cortical-SSM, an extension of Deep SSM that captures integrated EEG dependencies across temporal, spatial, and frequency domains.

• For frequency domain feature extraction, we introduce Wavelet-Convolution, which integrates deterministically obtained frequency components with adaptively derived frequency features. This approach enables the extraction of interpretable features while preserving their learnable representations.

2 Related Work↩︎

Deep learning-based methods for EEG decoding have been extensively studied, as reviewed by Abibullaev et al. [10] and Altaheri et al. [11]. Additionally, predictive methods for multivariate time series have been systematically summarized by Lara et al. [12] and Liang et al. [13]. Furthermore, deep state space models (Deep SSMs) have recently emerged as a promising architecture for sequence modeling, with advances documented by Patro et al. [14] and Wang et al. [15].

2.0.0.1 Deep learning-based EEG decoding.

Deep learning-based methods for decoding EEG signals have been extensively investigated [6], [16], and they demonstrate considerable potential for enhancing communication and rehabilitation in patients with physical paralysis [10], [11]. Early attempts rely on convolutional architectures [5], [17]. Nonetheless, the limited receptive field of convolutional layers has motivated the exploration of Transformer-based approaches [6], [7] to model temporal dependencies in EEG signals. Medformer [7], for instance, segments input signals into patches of varying temporal lengths and then feeds these patches into a Transformer. These Transformer-based approaches typically employ temporal patching strategies prior to Transformer processing to address computational complexity in long-sequence modeling. However, this preprocessing step may disrupt fine-grained temporal dependencies in EEG signals.

2.0.0.2 Multivariate time-series forecasting models.

Multivariate time-series forecasting has found broad applications across diverse domains, ranging from medical signal prediction (e.g., electromyography and electrocardiography analysis [18], [19]) as well as other applications [20], [21]. Notably, Transformer-based approaches are widely studied due to their demonstrated capability for long-range sequence modeling (e.g., [22], [23]). For example, Informer [22] proposed an efficient forecasting method based on a ProbSparse self-attention mechanism and a generative decoder. Other notable examples include PatchTST [24], which divides time-series signals into smaller patches for Transformer input, and Crossformer [23], which introduces a two-stage attention mechanism to capture inter-variable and temporal dependencies separately. In contrast, DLinear [25] achieves performance comparable to Transformer-based approaches by utilizing a simple yet effective MLP-based method that decomposes time-series into seasonal and trend components. Drawing on these findings, iTransformer [26] proposes a strategy in which the matrix operations in the Transformer are transposed, thereby capturing inter-variable dependencies through the Attention mechanism while modeling temporal dependencies with an MLP.

2.0.0.3 Deep state space models.

Although the Transformer architecture has been widely adopted across various domains [27], [28], the quadratic computational complexity \(O(N^2)\) of its attention mechanism with respect to sequence length \(N\) introduces fundamental scalability bottlenecks. This limitation has motivated extensive research into efficient alternative architectures [29], [30]. Among these alternatives, Deep SSMs [9], [31] have emerged as a promising framework for efficiently capturing long-range dependencies. In line with this trend, EEG classification models [32], [33] have also leveraged Deep SSMs to capture temporal dependencies in signals. However, current models predominantly adopt Mamba [31] as their foundational Deep SSM, leaving the optimal choice for the EEG classification task ambiguous. In contrast, our proposed method extends S5 [9], a Deep SSM explicitly designed to capture inter-variable dependencies within the state space, to effectively address the intrinsic multivariate nature of EEG signals.

3 Preliminaries↩︎

3.0.0.1 Deep state space models

Recent advancements in Deep SSMs [9], [34] have demonstrated their remarkable advantages over predominant architectures (including Transformer [35]) across various sequence modeling tasks. Inspired by classical SSM [36], Deep SSMs establish a principled framework in which input signals \(\mathbf{x}(t)\in\mathbb{R}^P\) are mapped to output signals \(\mathbf{y}(t)\in\mathbb{R}^P\) via latent states \(\mathbf{h}(t)\in\mathbb{R}^Q\), as follows:

\[\begin{align} \label{eq:mimo95ssm} \frac{d\mathbf{h}(t)}{dt} = \mathbf{A}\mathbf{h}(t) + \mathbf{B}\mathbf{x}(t) ,\qquad \mathbf{y}(t) = \mathbf{C}\mathbf{h}(t) + \mathbf{D}\mathbf{x}(t), \end{align}\tag{1}\]

where \(P\) and \(Q\) denote the number of variables for the input/output signals and the latent states, respectively. Moreover, \(\mathbf{A}\in\mathbb{R}^{Q\times Q}\) represents the state matrix, while \(\mathbf{B}\in\mathbb{R}^{Q\times P}\), \(\mathbf{C}\in\mathbb{R}^{P\times Q}\), and \(\mathbf{D}\in\mathbb{R}^{P\times P}\) denote the projection matrices. Notably, the Deep SSM variant S5 [9] has demonstrated remarkable effectiveness in capturing sequential relationships for continuous signals. Here, we detail the processing steps employed in S5.

In S5, the HiPPO-N matrix [37] is adopted as \(\mathbf{A}\) to effectively capture long-range dependencies in sequential signals. Since the HiPPO-N matrix is real symmetric, it can be diagonalized as \(\mathbf{A}=\mathbf{V\Lambda V}^{-1}\), thereby transforming Equation (1 ) into the following form: \[\begin{align} \label{eq:mimo95ssm2} \frac{d\tilde{\mathbf{h}}(t)}{dt} = \mathbf{\Lambda}\mathbf{\tilde{h}}(t) + \mathbf{\tilde{B}}\mathbf{x}(t),\qquad\mathbf{y}(t)=\mathbf{\tilde{C}}\mathbf{\tilde{h}}(t) + \mathbf{D}\mathbf{x}(t), \end{align}\tag{2}\]

where \(\tilde{\mathbf{h}}(t)=\mathbf{V}^{-1}\mathbf{h}(t)\), \(\tilde{\mathbf{B}}=\mathbf{V}^{-1}\mathbf{B}\), and \(\tilde{\mathbf{C}}=\mathbf{CV}\). Moreover, by introducing a timescale parameter \(\mathbf{\Delta}\in\mathbb{R}_+\), Equation (2 ) is discretized using the zero-order hold (ZOH) method [38] as follows: \[\begin{align} \label{eq:mimo95ssm3} \mathbf{\tilde{h}}_t = \mathbf{\bar{\Lambda}}\mathbf{\tilde{h}}_{t-1} + \mathbf{\bar{B}}\mathbf{x}_t,\qquad\mathbf{y}_t=\mathbf{\bar{C}}\mathbf{\tilde{h}}_t + \mathbf{\bar{D}}\mathbf{x}_t, \end{align}\tag{3}\] where \(\bar{\mathbf{\Lambda}}=\mathrm{exp}(\mathbf{\Lambda}\mathbf{\Delta})\), \(\bar{\mathbf{B}}=\mathbf{\Lambda}^{-1}\left(\bar{\mathbf{\Lambda}}-\mathbf{I}\right)\tilde{\mathbf{B}}\), \(\bar{\mathbf{C}}=\tilde{\mathbf{C}}\), \(\bar{\mathbf{D}}=\mathbf{D}\). In practice, \(\mathbf{\Delta}\in\mathbb{R}^Q\) is used for the timescale parameter, and \(\mathbf{D}\) is restricted as a diagonal matrix. Under these conditions, the learnable parameters consist of \(\mathrm{diag}(\mathbf{\Lambda})\), \(\tilde{\mathbf{B}}\), \(\tilde{\mathbf{C}}\), \(\mathrm{diag}(\mathbf{D})\), and \(\mathbf{\Delta}\). Furthermore, S5 achieves efficient modeling of Equation (3 ) through the introduction of parallel scanning.

4 Method↩︎

4.1 Cortical-SSM↩︎

In this study, we address a classification task for EEG signals recorded while subjects engage in motor imagery (MI) tasks. In this task, it is desirable to accurately classify the corresponding EEG signals recorded while subjects engage in MI involving \(N\) types of movement. The input is EEG signals recorded while the subjects performed MI, and the output is a predicted probability corresponding to the MI.

Unlike previous studies [32], [33], our method jointly models spatio-temporal features within individual frequency components and temporal-frequency features within individual electrodes. This parallel design yields a comprehensive multi-domain representation for EEG classification. Because the architecture is not specific to MI, it can also be applied to decoding other endogenous and exogenous stimuli from EEG signals.

Figure 2 shows the overall architecture of our proposed Cortical-SSM. The proposed method consists of three main modules, Wavelet-Convolution, Frequency-SSM, and Channel-SSM. The input \(\mathbf{X}\in\mathbb{R}^{M\times T}\) to our model is the recorded EEG signal, where \(M\) and \(T\) denote the number of electrodes and the sequence length, respectively.

4.2 Wavelet-Convolution↩︎

In the Wavelet-Convolution module, we integrate deterministically derived frequency components with learned frequency features. Previous models for EEG classification widely employed one-dimensional convolutional layers to learn frequency features through their kernels [5], [6]. Despite the promising results achieved by these models, their black-box nature poses challenges for direct interpretability and impedes frequency domain analysis of the extracted features. In contrast, the conventional short-time Fourier transform (STFT) and wavelet transform extract frequency components in a deterministic manner, thereby yielding inherently interpretable features. However, these fixed representations constrain the model’s performance (see Subsection 6.3 for details).

To address this trade-off between learnability and interpretability, we propose a module integrating deterministically derived frequency components from an Explanation Branch (E-Branch) with trainable frequency features from an Adaptation Branch (A-Branch). This architecture maintains both the learnability of these representations and the extraction of non-black-box features. Specifically, we employ the continuous wavelet transform (CWT) for the E-Branch (first term of Equation (4 )), and a one-dimensional convolutional layer for the A-Branch (second term of Equation (4 )). Unlike existing methods that integrate the wavelet transform within convolutional layers [39], [40], our module processes CWT and 1D convolutional features in parallel, thereby enabling effective feature fusion.

Let \(\mathbf{x}_m\in\mathbb{R}^{T}\) denote the EEG signal recorded from the \(m\)-th electrode (\(m=1\ldots M\)). The module’s output \(\mathbf{\tilde{x}}_m\in\mathbb{R}^{F\times T}\) for each \(\mathbf{x}_m\) is formulated as shown below, yielding \(\mathbf{\tilde{X}} = \left\{\mathbf{\tilde{x}}_m\;\middle|\;m=1\ldots M\right\}\in\mathbb{R}^{M\times F\times T}\): \[\begin{align} \tilde{\mathbf{x}}_m = \frac{1}{2}\mathrm{LayerNorm}\left(\mathrm{CWT}\left(\mathbf{x}_m\right)\right) + \frac{1}{2}\mathrm{LayerNorm}\left(\mathrm{Conv1D}\left(\mathbf{x}_m\right)\right), \label{eq:wavelet-conv} \end{align}\tag{4}\] where \(\mathrm{Conv1D}(\cdot)\), \(\mathrm{CWT}(\cdot)\), and \(\mathrm{LayerNorm}(\cdot)\) represent one-dimensional convolutional layers, CWT, and layer normalization [41], respectively. Additionally, \(F\) represents the dimension in the frequency domain. Details of the processing are described below.

First, in \(\mathrm{Conv1D}(\cdot)\), the kernel length \(K\) is set to \(K=f_\mathrm{sample}/2\) following EEGNet [5], and frequency features are extracted. Here, \(f_\mathrm{sample}\) represents the sampling frequency.

In \(\mathrm{CWT}(\cdot)\), a filter bank \(\mathbf{\Psi}\in\mathbb{R}^{F\times T}\) is constructed to extract features for individual frequency components. For the extraction of frequency components, the Morlet wavelet [42] is employed as the standard mother wavelet. Let \(\boldsymbol{\psi}_f=\left[\psi_{f,-\frac{T}{2}},\psi_{f,-\frac{T}{2}+1},\;\ldots,\;\psi_{f,\frac{T}{2}}\right]\in\mathbb{R}^T\) denote the Morlet wavelet used to extract a frequency component \(f\). \(\psi_{f,t}\) is defined as follows: \[\begin{align} \psi_{f,t}=\sqrt{\frac{1}{s}}\pi^{-\frac{1}{4}}e^{i\omega_{0}\frac{t}{s}}e^{-\frac{1}{2}{\left(\frac{t}{s}\right)}^2},\;s = \frac{\omega_0 f_\mathrm{sample}}{2\pi f}, \end{align}\] where \(\omega_0\) and \(s\) denote the center frequency and scaling factor, respectively. Based on the above, \(\mathbf{\Psi}\) is expressed by the following equation: \[\begin{align} \mathbf{\Psi}=\left\{\psi_f\;\middle|\;f=f_{\mathrm{min}}+\alpha (f_{\mathrm{max}} - f_{\mathrm{min}})/F,\;\alpha=1\ldots F\right\}, \end{align}\] where \(f_\mathrm{min}\) and \(f_\mathrm{max}\) denote the minimum and maximum \(f\) targeted by the CWT, respectively.

Finally, we employ \(\mathrm{LayerNorm}(\cdot)\) to normalize the features obtained from \(\mathrm{Conv1D}(\cdot)\) and \(\mathrm{CWT}(\cdot)\). The learnable affine parameters in \(\mathrm{LayerNorm}(\cdot)\) adaptively modulate the feature distributions of both branches, effectively acting as an implicit gating mechanism. This enables channel-wise and frequency-wise adaptive weighting prior to fusion, rather than enforcing a fixed combination. Furthermore, \(\mathrm{LayerNorm}(\cdot)\) is specifically applied along the temporal dimension, as explained below. Normalization techniques widely applied in deep learning include batch normalization [43], group normalization [44], and layer normalization across dimensions at a given time step. These methods normalize features across variables. However, in multivariate time-series, when an arbitrary event occurs across variables and its effects appear at different time steps in the sequence, normalization across variables is known to introduce mutual noise [45], [46]. Similarly, in EEG classification tasks, when internal stimuli derived from motor activity are recorded as signal sources via electrodes, these stimuli influence different temporal points within each electrode’s signal. Therefore, we adopt layer normalization along the temporal dimension to address these issues.

4.3 Frequency-SSM↩︎

The Frequency-SSM module independently captures spatio-temporal feature interactions within each individual frequency component. Previous EEG studies [47], [48] have reported that MI tasks elicit frequency-specific power variations localized to functionally relevant cortical regions. Motivated by these neurophysiological findings, Frequency-SSM explicitly models the spatio-temporal dependencies for each frequency component in an independent manner. This design effectively tracks the power variations of those frequency bands associated with MI. The proposed module comprises \(L\) hierarchically organized blocks, where each block includes layer normalization, a feed-forward network, and a Deep SSM. In the following, we provide a detailed description of the processing steps applied to the input \(\mathbf{U}^{(l)}\) within the \(l\)-th block (\(l=0,\ldots ,L\)), where \(\mathbf{U}^{(0)}=\tilde{\mathbf{X}}\).

We first normalize \(\mathbf{U}^{(l)}\) along the temporal dimension for each frequency component \(f\in\{1,\ldots,F\}\): \[\begin{align} \tilde{\mathbf{u}}_f^{(l)} = \mathrm{LayerNorm}(\mathbf{u}_f^{(l)}). \end{align}\] Next, we capture the temporal dependencies in \(\tilde{\mathbf{u}}_f^{(l)}\). In EEG classification tasks, Transformer-based approaches are widely employed to model temporal dependencies [7], [16]. However, for long-sequence modeling, Deep SSMs have demonstrated superior performance to Transformer [31], [37]. Therefore, our module incorporates a Deep SSM to capture temporal dependencies in EEG signals. Deep SSMs can be broadly categorized into time-invariant [9], [49] and time-varying [31], [34] systems. Although time-varying Deep SSMs are prevalent in EEG classification [32], [33], prior work [31] suggests that their inherent selection mechanisms can be detrimental for continuous signals, and they have exhibited inferior performance in certain speech synthesis tasks. Therefore, we opt for time-invariant Deep SSMs. Furthermore, Deep SSMs can be classified by their input-output configuration into Single-Input Single-Output (SISO) [34], [37] or Multi-Input Multi-Output (MIMO) [9], [50]. Given the multi-electrode nature of EEG signals, a MIMO configuration that preserves inter-variable dependencies within the state space is deemed appropriate. For these reasons, we extend S5 [9] as the Deep SSM with a time-invariant and MIMO configuration.

By defining the operation of S5 as \(\mathrm{SSM}(\cdot)\), the feature \(\mathbf{u}_f^{(l+1)}\) that captures temporal dependencies in \(\tilde{\mathbf{u}}_f^{(l)}\) is obtained as the following equation: \[\begin{align} \mathbf{u}_f^{(l+1)} = \mathrm{FFN}\left(\mathrm{SSM}\left(\tilde{\mathbf{u}}_f^{(l)}\right)\right) + \tilde{\mathbf{u}}_f^{(l)}, \end{align}\] where \(\mathrm{FFN}(\cdot)\) denotes the feed-forward network. Subsequently, the output \(\mathbf{{U}}^{(l+1)} = \left[ \mathbf{{u}}^{(l+1)}_1,\;\mathbf{{u}}^{(l+1)}_2,\;\ldots\;\mathbf{{u}}^{(l+1)}_F\right]\in\mathbb{R}^{M\times F\times T}\) of the module at the \(l\)-th block is obtained, capturing temporal dependencies independently for each frequency component.

4.4 Channel-SSM↩︎

In the Channel-SSM, temporal-frequency features are extracted independently for each electrode. By explicitly modeling electrode-specific dependencies of temporal-frequency features, we capture localized variations in signal intensity associated with MI. This module comprises \(L\) stacked blocks, each incorporating layer normalization, a feed-forward network, and Deep SSM layers hierarchically. The input \(\mathbf{V}^{(l)}\) within the \(l\)-th block (\(l=0,\ldots ,L\)) is modeled as follows, where \(\mathbf{V}^{(0)}=\tilde{\mathbf{X}}\): \[\begin{align} \mathbf{v}_m^{(l+1)} = \mathrm{FFN}\left(\mathrm{SSM}\left(\mathrm{LayerNorm}\left({\mathbf{v}}_m^{(l)}\right)\right)\right) + \mathrm{LayerNorm}\left({\mathbf{v}}_m^{(l)}\right), \end{align}\] yielding the block’s final output \(\mathbf{{V}}^{(l+1)} = \left[ \mathbf{{v}}^{(l+1)}_1,\;\mathbf{{v}}^{(l+1)}_2,\;\ldots\;\mathbf{{v}}^{(l+1)}_M\right]\in\mathbb{R}^{M\times F\times T}\).

Finally, the outputs \(\mathbf{U}^{(L)}\) and \(\mathbf{V}^{(L)}\) from Frequency-SSM and Channel-SSM, respectively, are integrated as follows, yielding the predicted probability \(p(\hat{\mathbf{y}})\in\mathbb{R}^{N}\) of action corresponding to \(\mathbf{X}\): \[\begin{align} p(\mathbf{\hat{y}}) = \mathrm{FFN}\left(\left[\mathrm{AvgPool}\left(\mathbf{U}^{(L)}\right);\mathrm{AvgPool}\left(\mathbf{V}^{(L)}\right)\right]\right), \end{align}\] where \(N\) represents the number of action types. Moreover, \(\mathrm{AvgPool}(\cdot)\) denotes the average pooling layer, which aggregates the input features along the temporal dimension following the approach proposed in S4 [37]. Furthermore, we use the cross-entropy loss as the loss function.

4.5 Visual Explanations↩︎

In our proposed method, we generate visual explanations in both the spatio-temporal and temporal-frequency domains through the following procedures. Our approach extends Grad-CAM [51] to generate visual explanations tailored for time-series signals. Grad-CAM is formulated as follows, generating a feature map \(\mathbf{Z}^n\in\mathbb{R}^{I\times J}\) that represents the explanation for the \(n\)-th output class from the gradient of the loss with respect to \(\hat{y}_n\): \[\begin{align} \alpha_r= \frac{1}{IJ}\sum_{i\in I}\sum_{j\in J}\frac{\partial \hat{y}^{(n)}}{\partial o_{r,i,j}}\;,\qquad \mathbf{Z}^{(n)}=\mathrm{ReLU}\left(\sum_{r\in R}\alpha_r \mathbf{o}_r\right), \end{align}\] where \(\hat{y}_n\) denotes the \(n\)-th element of \(\hat{\mathbf{y}}\). \(I\), \(J\), and \(R\) denote the number of vertical and horizontal pixels, and the number of dimensions of the feature map, respectively. Furthermore, \(\mathbf{o}_r\in\mathbb{R}^{I\times J}\) and \(o_{r,i,j}\in\mathbb{R}\) denote the feature map of the \(r\)-th dimension and the feature indexed by position \(\left(i,j\right)\) within \(\mathbf{o}_r\), respectively. In this context, \(\mathbf{o}_r\) is required to retain the spatial relationships of the input.

Next, we describe our method for generating visual explanations in the spatio-temporal domain. While Grad-CAM operates on 2D feature maps \(\mathbf{o}_r\) and performs weighting across \(\left\{\mathbf{o}_r\right\}_{r=1}^R\), our method deals with time-series signals and therefore performs weighting across 1D time-series vectors. Here, we generate visual explanations using \(\mathbf{V}^{(L)}\), the feature closest to the output layer, which preserves the spatio-temporal relationships of the input signal. As detailed in Subsection 4.4, the features within \(\mathbf{V}^{(L)}\) are processed independently for each electrode. Therefore, we generate a separate visual explanation for each electrode by weighting across \(\left\{\mathbf{v}_{f,m}^{(L)}\right\}_{f=1}^F\). This yields the feature map \(\mathbf{Z}_\mathrm{ch}^{n}\in\mathbb{R}^{M\times T}\), representing the spatio-temporal visual explanation for the \(n\)-th output class, obtained as follows: \[\begin{align} \tag{5} \alpha^n_{f,m}&=&\frac{1}{T}\sum_{t\in T}\frac{\partial \hat{y}_n}{\partial {v}^{(L)}_{f,m,t}},\\ \tag{6} \mathbf{Z}^n_{\mathrm{ch}}&=&\left\{\mathrm{ReLU}\left(\sum_{f\in F}\alpha^n_{f,m}\mathbf{v}^{(L)}_{f,m}\right)\;\middle|\;m=1\ldots M\right\}. \end{align}\]

Finally, we describe our method for generating visual explanations in the temporal-frequency domain. Here, we generate visual explanations using \(\mathbf{U}^{(L)}\), the feature closest to the output layer, which preserves the temporal-frequency relationships of the input signal. As detailed in Subsection 4.3, in contrast to \(\mathbf{V}^{(L)}\), the features within \(\mathbf{U}^{(L)}\) are processed independently for each frequency component. Therefore, we generate a separate visual explanation for each frequency component by weighting across \(\left\{\mathbf{u}_{f,m}^{(L)}\right\}_{m=1}^M\). This yields the feature map \(\mathbf{Z}_\mathrm{freq}^{n}\in\mathbb{R}^{F\times T}\), representing the temporal-frequency visual explanation for the \(n\)-th output class, as in Equations (5 ) and (6 ).

5 Experiments↩︎

5.1 Datasets and Data Pre-Processing↩︎

In our experiments, we employed two publicly available MI EEG datasets, OpenBMI [4] and Stieger2021 [61]. These datasets were chosen because previous studies have demonstrated substantial domain shifts across subjects in MI EEG recordings [62], [63], making them suitable for evaluating the generalizability and robustness of decoding models. Both datasets include multiple sessions collected from more than fifty participants, enabling systematic subject-independent evaluation. Accordingly, we adopted a cross-subject \(k\)-fold cross-validation scheme (\(k=8\)).

For preprocessing, we followed the minimal pipeline described in [64]. EEG signals were downsampled to 250 Hz and used directly as model inputs without explicit noise or artifact removal. Training, validation, and test sets were constructed following the procedures described below and were used for model training, hyperparameter tuning, and performance evaluation, respectively. During training, we monitored classification accuracy on the validation set at each epoch, and model parameters achieving the highest validation accuracy were retained for final evaluation on the held-out test set. Further details regarding each dataset are provided below.

5.1.0.1 OpenBMI.

OpenBMI is a public EEG dataset containing recordings from 54 healthy participants engaging in a two-class MI task involving left- and right-hand grasping [4]. The experimental protocol follows [47], where visual cues indicate the movement to be imagined. Each subject performed two sessions consisting of 400 MI trials per session, resulting in 21,600 samples recorded at 1000 Hz using 62 Ag–AgCl electrodes positioned according to the international 10–20 system [65]. We extract the 4-second MI period for each trial. For cross-validation (\(k=8\)), subjects are randomly partitioned into training (44 subjects), validation (5 subjects), and test (5 subjects), corresponding to 17,600, 2,000, and 2,000 samples, respectively.

5.1.0.2 Stieger2021.

Stieger2021 contains EEG recordings collected from 62 participants performing MI-based cursor control tasks [61]. In this work, we focus on the left/right (LR) task, which involves imagined left- and right-hand grasping. Cursor movement was controlled using alpha power extracted from electrodes C3 and C4, and visual cues indicating the intended movement direction were presented on a screen. Across 7–11 sessions, each participant contributed between 1,050 and 1,650 MI trials, recorded using 64 electrodes (10–10 system) at 1000 Hz. Following prior work, we focus on 202,950 samples from 41 participants who completed all sessions, and use 4-second windows centered around MI onset. Subject-independent \(k\)-fold cross-validation was conducted with splits of 33, 4, and 4 subjects for training, validation, and testing, yielding 54,450, 6,600, and 6,600 samples, respectively.

5.2 Baseline Methods↩︎

This task involves (i) classification of EEG signals, and (ii) handling of multivariate time-series signals recorded from multiple electrodes. Accordingly, we selected baseline methods from the following perspectives:

• We chose baselines that have been successfully applied to EEG classification, similar to this task, including Shallow ConvNet [59], Deep ConvNet [59], EEGNet [5], TSception [60], EEG Conformer [6], and Medformer [7].

• Considering the multivariate nature of EEG signals, we selected baselines that have demonstrated effectiveness in multivariate time-series forecasting, such as Informer [22], Autoformer [52], FEDformer [53], Crossformer [23], DLinear [25], TimesNet [54], PatchTST [24], TimeMixer [55], iTransformer [26], S-Mamba [58], and TimeMachine [57].

5.3 Implementation Details↩︎

Table ¿tbl:tab:hyperparams? summarizes the experimental settings of the proposed method. The hyperparameters were fixed a priori such that the number of trainable parameters was comparable to that of the baselines.

Table ¿tbl:tab:model95param95comp? further reports the number of trainable parameters and the inference time for the proposed method and all baselines. The proposed model contained approximately 0.83 million parameters, with an inference time of 2.20 ms per sample, demonstrating computational cost comparable to the baselines.

All experiments were conducted on an NVIDIA GeForce RTX 4090 GPU (24 GB VRAM) and an Intel Core i9-13900KF CPU with 64 GB RAM. The training time per task for the proposed method was approximately 157 min on OpenBMI and 50 min on Stieger2021.

6 Results↩︎

6.1 Quantitative Results↩︎

We conducted experiments to compare the performance of our model with baselines on the OpenBMI [4] and Stieger2021 [61] datasets. Quantitative results are summarized in Table ¿tbl:tab:quantitative95main?. Values reported in the table represent the mean and standard deviation obtained across \(k\)-fold cross-validation (\(k=8\)). For the evaluation metrics, we employed accuracy, Macro-F1, AUROC (macro-averaged), AUPRC (macro-averaged), and Cohen’s Kappa. We used these metrics because they are standard for the classification of EEG signals.

Table ¿tbl:tab:quantitative95main? presents the performance of our proposed method, which achieved the highest accuracy, macro-F1, AUROC, AUPRC, and Cohen’s Kappa on OpenBMI, with values of 81.62%, 81.56%, 90.14%, 90.17%, and 0.63 points, respectively. The proposed method outperformed the second-best model by 3.72, 3.78, 4.33, 4.81, and 0.07 points for these respective metrics. Moreover, the proposed method achieved state-of-the-art results on Stieger2021 with accuracy, macro-F1, AUROC, AUPRC, and Cohen’s Kappa of 87.12%, 87.11%, 94.66%, 94.82%, and 0.74 points, respectively. Compared with the second-best model, our method exhibited improvements of 4.49, 4.51, 1.86, 2.73, and 0.09 points for accuracy, macro-F1, AUROC, AUPRC, and Cohen’s Kappa, respectively.

Furthermore, Figure 4 reports subject-wise decoding performance on both benchmarks. Specifically, for each dataset, we compare the proposed method with the top two baseline models identified in Table 1 in terms of classification accuracy. Each point corresponds to an individual subject, while violin and box plots summarize the distribution across subjects, enabling a comprehensive assessment of inter-subject variability beyond aggregated statistics such as mean and standard deviation.

As shown in Figure 4, the proposed method consistently outperforms the baselines for the majority of subjects across both datasets. To quantitatively assess these differences, we apply the Wilcoxon signed-rank test, a nonparametric test for paired samples, and observe statistically significant improvements across all evaluation metrics and both benchmarks (\(p < 0.05\)). The use of a nonparametric test is justified by a normality analysis: Shapiro–Wilk tests conducted on fold-wise accuracy values reject the normality assumption (\(p < 0.05\)), likely due to subject-specific variability and domain shifts inherent in EEG signals [62], [63]. These deviations from normality violate the assumptions of parametric tests such as the \(t\)-test, thereby necessitating a nonparametric alternative. Overall, these results indicate that the observed performance gains are not driven by a small subset of subjects but are consistently achieved across the population.

6.2 Interpretability↩︎

We further examine the proposed method by generating visual explanations, following the procedure described in Section 4.5. Figure 3 presents sample-agnostic visual explanations on the OpenBMI dataset [4]. Column (a) corresponds to left-hand grasp imagery, while Column (b) corresponds to right-hand grasp imagery. Rows (i)–(iv) show representative subjects (S2, S9, S15, and S25, respectively). For each subject and class, the upper panels depict spatio-temporal topographic visualizations at 0.5 s, 1.5 s, 2.5 s, and 3.5 s relative to MI onset, obtained by projecting \(\mathbf{Z}_\mathrm{ch}^n\) (introduced in Section 4.5) onto scalp maps aligned with the 10–20 system [65]. The lower panels illustrate the corresponding temporal–frequency visualizations, representing \(\mathbf{Z}_\mathrm{freq}^n\), where 0 seconds denotes the onset of motor imagery. Across subjects and both imagery classes, the temporal–frequency maps consistently highlight activity in the mu band (approximately 10 Hz), which is well known to be associated with motor imagery in EEG [8]. Moreover, the topographic maps show strong attention around electrodes C3 and C4, which are located over the motor cortex and are neurophysiologically linked to hand motor control [8]. These observations indicate that the proposed method reliably attends to neurophysiologically meaningful frequency bands and spatial regions, in a sample-agnostic manner.

Furthermore, to quantitatively assess the explanations, we evaluate their causal faithfulness using insertion and deletion curves [66], as shown in Fig 5. Insertion and deletion curves are standard metrics for evaluating the faithfulness of visual explanations. They quantify how the model’s prediction changes as regions identified as important by a saliency map are progressively inserted or removed. Formally, given an input \(\boldsymbol{X}\), we define a sequence of perturbed inputs \(\{\boldsymbol{X}_k\}\) by inserting (or deleting) the top \(k~\%\) of features ranked by the saliency map. For each \(\boldsymbol{X}_k\), we compute the predicted probability of the ground-truth class, \(p(\hat{y}_k=y)\). The insertion (or deletion) curve is then defined as the trajectory of \(p(\hat{y}_k=y)\) as a function of \(k\). A faithful explanation is expected to yield a steep increase in the insertion curve and a steep decrease in the deletion curve for small \(k\), indicating that the identified features are indeed critical for the model’s prediction.

In Fig. 5, panels (a) and (b) show the insertion (blue) and deletion (red) curves, respectively, while panels (i) and (ii) correspond to results on 100 samples from OpenBMI and Stieger2021. As a control, we also reported curves obtained by inserting or deleting randomly selected features (gray). The shaded regions indicate variability across samples, corresponding to ±3 standard deviations. Performance was quantified using the area under the curve (AUC). Since all benchmark tasks considered in this study are binary classification tasks, the prediction probability was evaluated within the range [0.5,1.0], where 0.5 corresponds to random performance. The results showed that the insertion curves consistently achieved higher AUC than the random baseline, while the deletion curves yielded lower AUC, demonstrating a sharper rise (or drop) in prediction confidence for small \(k\). Furthermore, the differences between the proposed method and the random baseline were statistically significant under the Wilcoxon signed-rank test (\(p<0.001\)). These findings provide quantitative evidence that the features highlighted by the proposed explanations are not merely correlational but are causally relevant to the model’s predictions.

6.3 Ablation Study↩︎

To investigate the effectiveness of each module, we conducted ablation studies on the following three conditions.

6.3.0.1 Wavelet-Convolution ablation.

Table 1 presents the performance of different frequency domain feature extraction methods in the E-Branch and A-Branch of the Wavelet-Convolution module. We compared five model variants (1-i)-(1-v) with combinations of short-time Fourier transform (STFT) and continuous wavelet transform (CWT) in the E-Branch, along with 1D convolutional layers in the A-Branch. The results reveal that models employing CWT in the E-Branch exhibited markedly different performance on OpenBMI [4], with the classification accuracy of Model (1-ii) being 2.76 points lower than that of Model (1-v). Similarly, when models using 1D convolutions in the A-Branch were compared, Model (1-iii) underperformed Model (1-v) by 1.37 points in the corresponding metric.

To explain this observation, we consider the intrinsic properties of MI EEG signals, which exhibit both well-defined frequency characteristics and pronounced non-stationarity. Regarding frequency characteristics, MI EEG signals contain discriminative patterns in specific frequency bands, particularly the \(\mu\)-band (8–12 Hz), as widely reported in prior work [68] and confirmed by our visual explanations (Section 6.2). Accordingly, employing the Continuous Wavelet Transform (CWT) in the E-Branch enables deterministic extraction of such physiologically meaningful components. Moreover, the superiority of the CWT-based configuration (Model (1-v)) over the STFT-based variant (Model (1-iv)) can be attributed to the transient nature of \(\mu\)-band modulations, which are not well captured by the fixed temporal resolution of STFT. Furthermore, in terms of non-stationarity, EEG signals are affected by temporal electrode shifts, device characteristics, and subject-specific variability [69]. To account for these factors, the A-Branch employs one-dimensional convolutional layers that function as adaptive filters, enabling the model to learn task-relevant frequency patterns while compensating for such variations. Taken together, these results suggest that the integration of deterministic frequency decomposition (via CWT) and adaptive feature learning (via one-dimensional convolution) provides a complementary inductive bias that is well aligned with the characteristics of MI EEG signals.

6.3.0.2 Frequency-SSM and Channel-SSM ablation.

Table 2 shows the performance of different architectures in the Frequency-SSM and Channel-SSM. We compared models using the following architectures for capturing temporal dependencies: (2-i) Attention [35], (2-ii) S4-LegS [37], (2-iii) Mega [67], (2-iv) Mamba-2 [34], and (2-v) S5 [9]. Table 2 indicates that the classification accuracy of models (2-i), (2-ii), (2-iii), and (2-iv) on OpenBMI underperformed Model (2-v) by 2.88, 1.93, 1.76, and 1.33 points, respectively. These results suggest that S5, a Deep SSM employing a time-invariant system and MIMO configuration, is effective for capturing temporal dependencies in EEG signals.

6.3.0.3 Module-wise ablation.

Table 3 presents the performance of the three main modules: Wavelet-Convolution, Frequency-SSM, and Channel-SSM. We compared four model configurations: (3-i) exclusion of Wavelet-Convolution, (3-ii) exclusion of Frequency-SSM, (3-iii) exclusion of Channel-SSM, and (3-iv) the complete model incorporating all three modules. Regarding classification accuracy on OpenBMI, Models (3-i), (3-ii), and (3-iii) underperformed Model (3-iv) by 5.74, 2.48, and 0.98 points, respectively. These findings indicate that each module contributes to improving overall model performance, with the Wavelet-Convolution module exerting the most significant impact.

6.4 Sensitivity Analysis↩︎

6.4.0.1 Sequence length analysis

To investigate the effect of sequence length on decoding performance, we performed a sensitivity analysis on the OpenBMI dataset by resampling EEG signals originally recorded at 1000 Hz to three sampling rates (250 Hz, 500 Hz, and 1000 Hz), as illustrated in Figure 6. The 250 Hz condition corresponds to the setting used in the original experimental setup, and all conditions employ 4-second EEG segments as input. We compared our method against two representative EEG decoding baselines: EEG Conformer [6] and Medformer [7]. While EEG Conformer and Medformer revealed accuracy drops of up to 7.34 and 3.73 points, respectively, our method shows only a marginal decrease of 1.06 points. These results suggest that the proposed method scales with sequence length.

6.4.0.2 Signal-to-noise ratio ablation

To evaluate the impact of noise and signal degradation, we conducted an additional experiment in which the signal-to-noise ratio (SNR) of EEG signals from the OpenBMI dataset was progressively degraded by -0 dB, -5 dB, -10 dB, and -15 dB as shown in Figure 7. Here, we compared the performance of our proposed method with the second-best baseline, EEG Conformer [6]. The signal bandwidth was set to 1-100 Hz. As shown in the figure, our method consistently outperformed EEG Conformer across all SNR levels in terms of classification accuracy and macro-F1, demonstrating its effectiveness even in the presence of substantial noise and signal degradation common in clinical settings.

7 Discussion↩︎

While the results presented in Table ¿tbl:tab:quantitative95main? indicate the potential advantages of the proposed method, it is important to note that the current formulation assumes a synchronous setting. In this study, synchronous data with explicitly controlled motor imagery timing were adopted to enable principled interpretation of the visual explanations. Because the onset of motor imagery is known in this setting, the temporal correspondence between model attributions and the underlying neural processes can be clearly interpreted. In practical applications, however, MI EEG decoding is typically required to operate in an asynchronous setting, where the onset of motor imagery is unknown. Extending the proposed method to support interpretable visual explanations under such asynchronous conditions therefore remains an important direction for future work.

In addition to this experimental constraint, as with the existing baselines, Cortical-SSM does not explicitly address subject-level domain shifts [62], [63]. As shown in Section 6.1, the Shapiro–Wilk test rejects normality across folds, suggesting substantial variance. Incorporating domain adaptation techniques, such as learning invariant representations offers a promising path forward.

From the architectural perspective, Cortical-SSM primarily fuses multi-domain features at a later stage of the model, following largely independent processing in earlier stages. While this design enhances interpretability, it may in some cases encourage greater reliance on a subset of domains (e.g., temporal) when they are sufficient for the task, potentially limiting the utilization of complementary information from other domains. Future work could explore more tightly coupled or progressive integration strategies to further promote balanced feature utilization.

Moreover, although the visual explanations of the proposed method highlight brain regions that are broadly consistent with established neurophysiological findings, caution is required when interpreting these patterns in causal terms. The visualization reflects gradients of the prediction error and does not directly measure neural activation itself. Accordingly, the generated explanations should be understood as emphasizing discriminative patterns that contribute to accurate MI classification, rather than as providing direct evidence of underlying neurophysiological mechanisms. Given this limitation, future work may also explore integrating the proposed method with Electrophysiological Source Imaging (ESI) approaches [70], [71], enabling a more direct linkage to cortical sources and providing insights into the causal structure of the decoding process.

Despite these considerations, the present findings suggest that jointly modeling the spatial, temporal, and spectral structure of EEG is a useful strategy for MI EEG decoding. Moreover, providing visual explanations appears particularly important for EEG-based brain functional imaging, where the inherently low signal-to-noise ratio makes physiologically plausible interpretation essential for reliable decoding.

8 Conclusion↩︎

In this study, we focused on a classification task based on EEG signals recorded during MI tasks. We proposed Cortical-SSM, an extension of Deep SSMs designed to capture integrated dependencies across temporal, spatial, and frequency domains. For frequency feature extraction, we proposed the Wavelet-Convolution, which extracts non-black-box frequency-analyzable features while maintaining a learnable representation. In comprehensive evaluations on two MI EEG benchmarks, our method consistently outperformed the comparison baseline methods. Furthermore, we demonstrated that for both EEG signals, neurophysiologically significant regions were attended to in the visual explanations generated by our proposed method.

References↩︎