Artifact Correction#

MNE-RT provides five real-time artifact correction methods that can be selected via RTStream (artifact_correction= parameter) or the CLI (--artifact-correction). All methods expose a common fit / transform interface and operate sample-by-sample (or chunk-by-chunk) during the closed-loop session.


Adaptive LMS filter#

Class: mne_rt.tools.AdaptiveLMSFilter

The Least Mean Squares (LMS) algorithm [1] is a stochastic gradient-descent adaptive filter designed to cancel a reference signal (e.g., an EOG or ECG channel) from the M/EEG data channels.

Model. For channel \(i\) at time \(t\), the cleaned signal is:

\[y_i(t) = x_i(t) - \mathbf{w}_i(t)^\top \mathbf{r}(t)\]

where \(x_i(t)\) is the raw sample, \(\mathbf{r}(t) \in \mathbb{R}^q\) is the reference signal vector (possibly delayed copies of the artefact channel), and \(\mathbf{w}_i(t)\) is the adaptive weight vector.

Weight update (Widrow–Hoff rule):

\[\mathbf{w}_i(t+1) = \mathbf{w}_i(t) + \mu\, y_i(t)\, \mathbf{r}(t)\]

where \(\mu > 0\) is the step size (learning rate). The algorithm converges when the cross-correlation between \(\mathbf{r}(t)\) and the residual \(y_i(t)\) reaches zero — i.e., when the filter has learned to predict the artifact contribution in channel \(i\) from the reference.

Stability condition. The step size must satisfy \(0 < \mu < 2/(\lambda_{\max}\, q)\), where \(\lambda_{\max}\) is the largest eigenvalue of the reference autocorrelation matrix. MNE-RT uses a normalised variant (NLMS) that adapts \(\mu\) automatically:

\[\mu_{\mathrm{norm}}(t) = \frac{\mu_0}{\epsilon + \|\mathbf{r}(t)\|^2}\]

When to use. Best suited for ocular (EOG) or cardiac (ECG) artifacts when a clean reference channel is available. Requires no calibration baseline — the filter adapts online from sample 1.


ORICA — Online Recursive ICA#

Class: mne_rt.tools.ORICA

ORICA [2] is an online (sample-recursive) variant of Independent Component Analysis (ICA) [3] that continuously updates the unmixing matrix \(\mathbf{W}\) as new data arrive.

ICA model. EEG is assumed to be a linear mixture of independent sources:

\[\mathbf{x}(t) = \mathbf{A}\, \mathbf{s}(t), \quad \mathbf{s}(t) = \mathbf{W}\, \mathbf{x}(t)\]

where \(\mathbf{A} \in \mathbb{R}^{p \times p}\) is the mixing matrix, \(\mathbf{W} = \mathbf{A}^{-1}\) is the unmixing matrix, and \(\mathbf{s}(t)\) are the independent components (ICs).

Recursive update. ORICA updates \(\mathbf{W}\) with each new sample using an exponentially weighted natural-gradient rule:

\[\mathbf{W}(t) \leftarrow \mathbf{W}(t) + \lambda(t)\, \bigl[\mathbf{I} - \mathbf{f}\!\left(\mathbf{s}(t)\right)\mathbf{s}(t)^\top\bigr]\, \mathbf{W}(t)\]

where \(\mathbf{f}(\cdot)\) is a nonlinear score function (e.g., \(f(s) = \tanh(s)\) for super-Gaussian sources) and \(\lambda(t)\) is a decreasing forgetting factor that controls how fast old data are discarded.

Artifact removal. ICs whose kurtosis, focal power, or spatial map match known artifact signatures (blinks, muscle, ECG) are identified and zeroed before back-projection:

\[\hat{\mathbf{x}}(t) = \mathbf{A}_{\text{clean}}\, \mathbf{W}_{\text{clean}}\, \mathbf{x}(t)\]

where the subscript “clean” denotes the subsets of columns / rows with artifact ICs removed.

When to use. Suitable for EEG with diverse artifact types. No reference channel is required. Because \(\mathbf{W}\) adapts continuously, ORICA tracks slow non-stationarities (electrode drift, changing impedances).


GEDAI — GED-based Artifact Isolation#

Class: mne_rt.tools.GEDAIDenoiser

Important

License notice. mne-rt itself is BSD-3-Clause, but AdaptiveMultibandGEDAIDenoiser (multiband=True / "causal" below) depends on the gedai package, which is distributed under the PolyForm Noncommercial License 1.0.0. That license restricts use to noncommercial purposes, and the underlying GEDAI algorithm is patent-pending with no commercial patent grant included. gedai is an opt-in extra — pip install "mne-rt[gedai]" — and is not installed by default or by mne-rt[full]; GEDAIDenoiser does not use it. Commercial use of the adaptive multiband GEDAI feature is not covered by mne-rt’s permissive license. See the “License” section of the project README for details.

GEDAI uses Generalized Eigendecomposition (GED) [4] to find spatial filters that maximally separate brain signal from broadband activity, optionally anchored to a leadfield-derived forward model [5].

GED formulation. MNE-RT uses a band-vs-broadband formulation. Let \(\mathbf{R}_\text{band}\) be the covariance of the band-limited signal (e.g. 8–30 Hz for alpha/beta) and \(\mathbf{R}_\text{broad}\) be the broadband covariance. The GED solves:

\[\mathbf{R}_\text{band}\,\mathbf{W} = \mathbf{R}_\text{broad}\,\mathbf{W}\,\boldsymbol{\Lambda}\]

The columns of \(\mathbf{W}\) are spatial filters and the eigenvalues \(\lambda_k = \mathbf{w}_k^\top \mathbf{R}_\text{band}\, \mathbf{w}_k \;/\; \mathbf{w}_k^\top \mathbf{R}_\text{broad}\, \mathbf{w}_k\) express the band-to-broadband power ratio along each filter direction.

Eigenvalue interpretation. Filters are sorted by \(\lambda_k\) in descending order:

  • Large \(\lambda_k\) — high band/broadband ratio → brain-like (targeted oscillation dominant).

  • Small \(\lambda_k\) — low band/broadband ratio → non-brain (broadband, artifact-like).

Artifact isolation. find_noise_components() returns the \(n_\text{noise}\) components with the smallest eigenvalues (least band-specific), which are then zeroed before back-projection:

\[\hat{\mathbf{X}} = \mathbf{W}_{\text{brain}}\, \mathbf{W}_{\text{brain}}^\dagger\, \mathbf{X}\]

where \(\mathbf{W}_{\text{brain}}\) is the subset of \(\mathbf{W}\) with the largest eigenvalues (brain-like directions retained).

Leadfield incorporation. When a forward model is available (fit_from_leadfield()), GEDAI replaces \(\mathbf{R}_\text{band}\) with the leadfield outer product \(\mathbf{L}\mathbf{L}^\top\), directly constraining filters to be consistent with neural generators.

Delayed Haar multiband mode. mne_rt.tools.AdaptiveMultibandGEDAIDenoiser (fit_gedai(multiband="causal") or multiband=True) uses FIR kernels obtained from pyGEDAI’s actual Haar MODWT reconstruction, rather than Butterworth crossovers. Every retained band is shifted to the same causal delay, preserving their relative phase, and only retained, cleaned bands are summed. Bands whose upper edge is at or below wavelet_low_cutoff are excluded entirely: no low-frequency residual is added back. This dyadic wavelet cutoff is not a brick-wall frequency filter.

The default is six cycles per band. Automatic level selection includes baseline length, and epoch sizes reuse pyGEDAI’s adaptive parameter helper. Each active band’s covariance uses a finite trailing epoch rather than an EMA with an unbounded tail. A chunk updates the covariance and its spatial projection; single-sample chunks do not erase the covariance. refit=True can re-solve the spatial filters independently of update_covariance. At a 0.5 Hz cutoff and a 20-second baseline, the slowest retained epoch is 15.36 seconds at 200 Hz or 12 seconds at 256 Hz.

The default broadband pass includes the input 0.1 Hz Haar prefilter, the low-cut prefilter, and the broadband model’s own 0.1 Hz prefilter, matching pyGEDAI’s unfiltered-input path. Broadband SENSAI calibration uses overlapping one-second epochs; per-wavelet calibration uses non-overlapping epochs. broadband_pass=False disables broadband cleaning and its additional prefilters, but retains the initial 0.1 Hz prefilter and band exclusion. In RTStream.fit_gedai, use causal_broadband_pass=False to select this wavelet-only configuration and recalibrate its per-band SENSAI sensitivities.

Latency is explicit. Symmetric Haar reconstruction needs future context relative to its output timestamp. The streaming implementation obtains that context by delaying output, never by reading future acquisition chunks. The fixed delay is available as delay_samples and latency. Input sample index t produces output corresponding to t - delay_samples. Startup output includes delayed calibration history. With the default broadband pass, 0.5 Hz cutoff, and a 20-second baseline, the delay is 12.78 seconds at 200 Hz and 17.984375 seconds at 256 Hz. This excludes computation time and is not near-zero-latency processing.

RTStream.record_main feeds causal GEDAI only previously unseen timestamped samples in chunks of at most 50 ms. Overlapping analysis windows are assembled from cleaned output, not passed repeatedly through the stateful filters. Startup calibration-history output is discarded; feedback begins once a full cleaned analysis window is available. Window and marker timestamps refer to the delayed source samples. Missing input samples stop acquisition with an explicit error rather than silently joining discontinuous history.

For example:

nf.fit_gedai(
    multiband="causal",
    cycles_per_wavelet=6,
    causal_wavelet_low_cutoff=0.5,
)
print(nf.gedai.latency)

For the low-delay, wavelet-only configuration:

nf.fit_gedai(
    multiband="causal",
    shrinkage=0.05,
    cycles_per_wavelet=6,
    causal_wavelet_low_cutoff=1.0,
    causal_prefilter_mode="none",
    causal_broadband_pass=False,
)

Optional shorter Haar cascades. The default causal_prefilter_mode="reference" retains the original prefilter cascade. "lowcut" removes the two 0.1 Hz stages and keeps one Haar highpass at the requested cutoff. "none" skips standalone prefilters entirely; ignored wavelet bands are still excluded from output. In the latter case, broadband GEDAI sees low-frequency activity before wavelet exclusion. Calibration follows the selected cascade, so changing modes requires refitting. These modes are not claims of offline pyGEDAI equivalence.

At a 1 Hz cutoff, the "lowcut" / "none" fixed delays are respectively 1.270 / 0.635 seconds at 200 Hz, 0.9921875 / 0.49609375 seconds at 256 Hz, and 1.020 / 0.510 seconds at 500 Hz (for a sufficiently long baseline). Changing cycles from six to four or two shortens covariance epochs, not these kernel delays when the retained bands are unchanged. Upsampling adds no signal information and requires its own resampling delay and computation to be accounted for.

Causal SENSAI calibration. The default n_artifact="auto" uses pyGEDAI’s SENSAI optimizer, not a fixed component count or an eigenvalue heuristic. During fit, each Haar baseline band is divided into epochs using cycles_per_wavelet at the band’s lower frequency (half the upper edge for the lowest band), with at least as many samples as channels when the baseline permits. Only complete epochs are used, or the available baseline if shorter than one epoch. Physical data and reference covariance scales are retained without unit-power normalization.

For every calibration epoch, the generalized eigenproblem is \(C_b w = \lambda R w\), with the same regularized reference used online. SENSAI chooses the threshold maximizing signal subspace similarity minus three times noise subspace similarity, using pyGEDAI’s EEG scoring convention and up to three reference principal components. It searches sensitivity values using pyGEDAI’s frequency-dependent bounds and default tolerance. Scoring and optimization are delegated to the installed pyGEDAI version; this module does not modify its SENSAI objective or optimizer. Record that dependency version for reproducibility and refit after upgrading. Historical benchmark results obtained with a different SENSAI implementation are not accuracy guarantees for this pipeline. The fitted sensitivity \(s_b\) stays fixed until refitting. For each local covariance spectrum, \(p_{98,b}(t)\) is recalculated as the 98th percentile of \(\log|\lambda| + 100\) over nonzero absolute eigenvalues. The current cutoff is

\[\theta_b(t) = \exp\left(\frac{105-s_b}{100}p_{98,b}(t) - 100\right).\]

Streaming rejects absolute eigenvalues at or above \(\theta_b(t)\). The eigenvalue cutoff is not frozen to its calibration value. The calibration cutoffs are available through eigenvalue_thresholds. The last used adaptive cutoffs are available through current_eigenvalue_thresholds. There is no clean ground-truth input to calibration. A zero-variance baseline band cannot be calibrated; a warning is emitted and automatic rejection is disabled for that band until refitting.

Selection overrides. A supplied artifact template takes precedence over automatic selection in the direct denoiser API. RTStream uses SENSAI for multiband GEDAI and does not implicitly pass an existing blink template. An explicit integer n_artifact retains fixed-count selection: the largest eigenvalues are rejected in leadfield mode (excess power relative to the model), and the smallest are rejected in band-filter mode (least band-specific). n_artifact=0 without a template disables component rejection, but not the prefilters, output delay, or sub-cutoff band exclusion. Fixed positive counts can remove clean activity. Without a leadfield, automatic SENSAI uses the baseline broadband covariance as its signal reference, so its validity depends on the baseline being representative; it is not a theoretical brain reference.

The uncleaned Haar reconstruction matches pyGEDAI after alignment, away from recording boundaries. Spatial cleaning still differs: offline pyGEDAI blends two half-epoch-shifted, cosine-weighted reconstructions, whereas this streaming implementation uses one projection from a trailing covariance epoch per acquired chunk. This remaining difference is intentional; an offline future-containing epoch is not available at the current output time. SENSAI does not guarantee preservation of every clean source, and processing deadlines and denoising quality should be checked for the intended montage and hardware. A target Pearson correlation is not guaranteed.

When to use. Most powerful when a resting-state baseline is available and the artifact has a stable spatial signature (cardiac, powerline). Requires a calibration recording (fit_from_raw() uses the first N seconds).


ASR — Artifact Subspace Reconstruction#

Class: mne_rt.tools.ASRDenoiser

ASR [6] [7] learns the covariance statistics of a clean baseline segment, then projects out components that deviate beyond a threshold in subsequent data windows.

Calibration. The baseline recording is divided into overlapping windows of length \(T_{\mathrm{win}}\). The sample covariance per window is:

\[\mathbf{C}_j = \frac{1}{n_j} \mathbf{X}_j \mathbf{X}_j^\top\]

Windows with the highest total power (top max_dropout_fraction fraction by \(\operatorname{tr}(\mathbf{C}_j)\)) are discarded as likely artifactual. The mean of the remaining clean covariances is eigendecomposed:

\[\mathbf{C}_0 = \mathbf{U}\,\mathrm{diag}(d_1,\dots,d_p)\,\mathbf{U}^\top\]

Per-component RMS thresholds are set as:

\[T_k = \sigma \cdot \sqrt{d_k}, \qquad k = 1,\dots,p\]

where \(\sigma\) = cutoff (default 5).

Cleaning (online). For each incoming window \(\mathbf{X} \in \mathbb{R}^{p \times n}\):

  1. Project to calibration-component space: \(\mathbf{Z} = \mathbf{U}^\top(\mathbf{X} - \bar{\mathbf{x}})\).

  2. Compute per-component RMS: \(r_k = \sqrt{\operatorname{mean}(Z_k^2)}\).

  3. Zero artifact-dominated components: \(Z_k \leftarrow 0\) if \(r_k > T_k\).

  4. Back-project and restore the mean: \(\hat{\mathbf{X}} = \mathbf{U}\mathbf{Z} + \bar{\mathbf{x}}\).

Because steps 1–4 reduce to two matrix multiplies, the per-chunk runtime is \(\mathcal{O}(p^2 n)\).

When to use. General-purpose EEG artifact suppressor. Robust to transient high-amplitude artifacts (muscle bursts, electrode pops). Requires a brief clean baseline (fit_asr() records 60 s by default).


RTMaxwellFilter — Real-Time SSS / tSSS#

Class: mne_rt.tools.RTMaxwellFilter

The Signal Space Separation (SSS) method [8] and its temporal extension (tSSS) [9] exploit the physics of quasi-static magnetic fields to separate brain-origin signals from external interference in MEG recordings.

Physical basis. Inside the sensor array the magnetic field satisfies the Laplace equation \(\nabla^2 \phi = 0\). Solutions are decomposed into a basis of spherical harmonic multipoles:

\[\phi(\mathbf{r}) = \sum_{\ell=1}^{L_{\mathrm{in}}} \sum_{m=-\ell}^{\ell} \alpha_{\ell m}\, \phi_{\ell m}^{\mathrm{in}}(\mathbf{r}) \;+\; \sum_{\ell=1}^{L_{\mathrm{ex}}} \sum_{m=-\ell}^{\ell} \beta_{\ell m}\, \phi_{\ell m}^{\mathrm{ex}}(\mathbf{r})\]

where \(\phi^{\mathrm{in}}_{\ell m}\) are the internal basis functions (sources inside the array — brain signals) and \(\phi^{\mathrm{ex}}_{\ell m}\) are the external basis functions (sources outside — environmental noise).

SSS projector. Let \(\mathbf{S}_{\mathrm{in}}\) and \(\mathbf{S}_{\mathrm{ex}}\) be the basis matrices evaluated at all sensor positions. The SSS projector onto the internal subspace is:

\[\mathbf{P}_{\mathrm{SSS}} = \mathbf{S}_{\mathrm{in}} \mathbf{S}_{\mathrm{in}}^\dagger\]

where \((\cdot)^\dagger\) denotes the Moore–Penrose pseudoinverse. Applied to raw MEG data \(\mathbf{X}\):

\[\hat{\mathbf{X}}_{\mathrm{SSS}} = \mathbf{P}_{\mathrm{SSS}}\, \mathbf{X}\]

This single matrix multiply (cached at calibration time) suppresses all signal components that cannot originate inside the sensor array.

tSSS — temporal extension. Internal components that correlate with external ones over a sliding time window indicate that interference has “leaked” into the internal space (e.g., due to a nearby ferromagnetic object). tSSS [9] removes these correlated components by projecting the internal coefficients onto the subspace orthogonal to the external ones:

\[\hat{\boldsymbol{\alpha}} = \bigl(\mathbf{I} - \Pi_{\mathrm{corr}}\bigr)\, \boldsymbol{\alpha}\]

where \(\Pi_{\mathrm{corr}}\) is the projector onto the directions of \(\boldsymbol{\alpha}\) that correlate above threshold with the external coefficient matrix over the buffer.

Real-time implementation in MNE-RT. RTMaxwellFilter caches \(\mathbf{P}_{\mathrm{SSS}}\) once during fit() using mne.preprocessing.compute_maxwell_basis(), then applies it as a matrix multiply per incoming chunk. This gives:

  • Zero additional latency beyond a matrix multiply per chunk.

  • Numerical equivalence with offline MNE maxwell_filter() (verified to machine-epsilon level).

tSSS runs periodically on a rolling buffer with maxwell_filter() (st_only=True, skipping the already-applied spatial step).

Empty-room noise covariance. When an empty-room recording is provided, RTMaxwellFilter uses a system-identification approach: MNE’s full Maxwell filter (with noise-informed regularisation) is applied to a Gaussian test signal, and the effective operator is recovered by least squares:

\[\mathbf{P}_{\mathrm{SSS}}^{\mathrm{ER}} = \mathbf{Y}\, \mathbf{X}^\dagger\]

where \(\mathbf{X}\) is the test input (good-channel MEG) and \(\mathbf{Y}\) is the filtered output. This transparently incorporates fine calibration, cross-talk compensation, and noise-informed regularisation into the single cached matrix — no re-implementation of SSS is required.

Baseline requirement. Unlike ASR or GEDAI, RTMaxwellFilter requires no baseline recording — \(\mathbf{P}_{\mathrm{SSS}}\) depends only on sensor geometry, which is fixed at scanner installation.

When to use. MEG data only. Mandatory first denoising step in any MEG pipeline. Use SSS mode for environmental noise suppression; add tSSS when nearby ferromagnetic objects or implants are present.


Bad Channel Detection#

Class: mne_rt.tools.BadChannelDetector

Artifact correction assumes all channels are functional. A disconnected or noisy electrode contaminating adjacent channels can degrade every downstream method. BadChannelDetector evaluates each incoming data window against up to four independent criteria and uses a rolling-window majority vote to flag persistently bad channels.

Criteria

"flat"

Channels whose RMS amplitude falls below flat_threshold. Catches disconnected leads and dead electrodes.

"variance"

Channels whose RMS is a robust outlier across all channels, measured by a MAD-based z-score: \(z = (\text{RMS} - \text{median})\;/\;(\text{MAD} \times 1.4826)\). Flags both excessively noisy and suspiciously quiet channels.

"correlation"

Channels whose mean Pearson correlation with their \(K\) nearest spatial neighbours falls below corr_threshold. Requires channel positions to be set in info (montage applied).

"hf_noise"

Channels with an abnormally high ratio of HF power (> hf_cutoff Hz) to broadband power, using the same MAD z-score. Catches EMG, loose cable noise, and intermittent contacts.

Rolling-window vote. A channel is declared bad only when it is flagged by at least one criterion in ≥ min_bad_frac of the last history_windows windows. With the default min_bad_frac=0.5 this is a majority vote over the past 30 windows, suppressing false positives from transient artifacts.

Usage example:

detector = BadChannelDetector(raw.info, method="all")
while streaming:
    window = stream.get_data(1.0)        # 1-second chunk
    bad_channels = detector.update(window)
    raw.info["bads"] = bad_channels

Riemannian Potato — Online Artifact Detection#

Class: mne_rt.tools.RiemannianPotatoDetector

The Riemannian Potato [10][11] is a covariance-based method for detecting — rather than correcting — artifactual data windows in real time. It operates on the manifold of symmetric positive-definite (SPD) matrices by computing the Riemannian distance from each incoming covariance matrix to the geometric mean of a clean calibration set. When this distance exceeds a z-score threshold, the window is flagged as an artifact and can be discarded or buffered for later ICA clean-up.

Why Riemannian geometry? Standard Euclidean distances on covariance matrices are sensitive to the arbitrary scaling of individual channels. The Riemannian (affine-invariant) metric is scale-invariant: multiplying all channels by any positive constant leaves the distance unchanged. This makes the potato robust to amplitude drift and electrode-impedance changes across sessions.

Algorithm. Let \(\Sigma_t\) be the sample covariance matrix of window \(t\) and \(\bar{\Sigma}\) be the geometric mean of calibration covariances. The z-score is:

\[z_t = \frac{\delta(\Sigma_t, \bar{\Sigma}) - \mu_\delta}{\sigma_\delta}\]

where \(\delta\) is the Riemannian distance and \(\mu_\delta\), \(\sigma_\delta\) are the mean and standard deviation of distances computed over the calibration set. Window \(t\) is declared an artifact when \(z_t > \theta\) (default \(\theta = 3\)).

Note

The Riemannian Potato is a detection method, not a correction method. It identifies and rejects artifactual windows; it does not reconstruct clean signal from contaminated data. Pair it with ORICA or ASRDenoiser when online correction is needed.

Usage example:

from mne_rt.tools import RiemannianPotatoDetector

# Calibrate on 60 clean 1-second windows
detector = RiemannianPotatoDetector(threshold=3.0)
detector.fit(clean_windows)       # shape (n_windows, n_ch, n_samples)

while streaming:
    window = stream.get_data(1.0)
    is_clean, z_score = detector.detect(window)
    if is_clean:
        compute_nf_feature(window)
    else:
        logger.debug("Artifact detected (z=%.2f) — window skipped", z_score)

When to use. Any EEG or MEG setup where transient artifacts (blinks, muscle bursts, electrode pops) need to be excluded from the feature computation. Particularly effective when reliable reference channels for LMS are unavailable. Requires a short clean calibration phase (≥ 30 windows recommended) before streaming begins.


Method comparison#

Method

Modality

Baseline needed

Reference ch.

Adapts online

Primary target

LMS

EEG / MEG

No

Yes (EOG/ECG)

Yes

Ocular, cardiac

ORICA

EEG / MEG

No

No

Yes

Mixed artifacts

GEDAI

EEG / MEG

Yes (N s)

No

No

Cardiac, powerline

ASR

EEG / MEG

Yes (N s)

No

No

Transient high-amplitude

RTMaxwellFilter

MEG only

No

No

No (periodic tSSS)

Environmental noise

BadChannelDetector

EEG / MEG

No

No

Yes (rolling)

Disconnected / noisy channels

RiemannianPotatoDetector

EEG / MEG

Yes (N windows)

No

No (detection only)

Transient artifacts (blinks, muscle, pops)


See also the artifact comparison example for a quantitative benchmark of LMS, ASR, and GEDAI on simulated EEG, and the RTMaxwell example for a demonstration of real-time SSS on the MNE sample MEG dataset.


References#