Plot sensor denoising using oversampled temporal projection

This demonstrates denoising using the OTP algorithm [1] on data with with sensor artifacts (flux jumps) and random noise.

# Author: Eric Larson <>
# License: BSD (3-clause)

import os.path as op
import mne
import numpy as np

from mne import find_events, fit_dipole
from mne.datasets.brainstorm import bst_phantom_elekta
from import read_raw_fif


Plot the phantom data, lowpassed to get rid of high-frequency artifacts. We also crop to a single 10-second segment for speed. Notice that there are two large flux jumps on channel 1522 that could spread to other channels when performing subsequent spatial operations (e.g., Maxwell filtering, SSP, or ICA).

dipole_number = 1
data_path = bst_phantom_elekta.data_path()
raw = read_raw_fif(
    op.join(data_path, 'kojak_all_200nAm_pp_no_chpi_no_ms_raw.fif'))
raw.crop(40., 50.).load_data()
order = list(range(160, 170))
raw.copy().filter(0., 40.).plot(order=order, n_channels=10)


Opening raw data file /home/circleci/mne_data/MNE-brainstorm-data/bst_phantom_elekta/kojak_all_200nAm_pp_no_chpi_no_ms_raw.fif...
    Read a total of 13 projection items:
        planar-0.0-115.0-PCA-01 (1 x 306)  idle
        planar-0.0-115.0-PCA-02 (1 x 306)  idle
        planar-0.0-115.0-PCA-03 (1 x 306)  idle
        planar-0.0-115.0-PCA-04 (1 x 306)  idle
        planar-0.0-115.0-PCA-05 (1 x 306)  idle
        axial-0.0-115.0-PCA-01 (1 x 306)  idle
        axial-0.0-115.0-PCA-02 (1 x 306)  idle
        axial-0.0-115.0-PCA-03 (1 x 306)  idle
        axial-0.0-115.0-PCA-04 (1 x 306)  idle
        axial-0.0-115.0-PCA-05 (1 x 306)  idle
        axial-0.0-115.0-PCA-06 (1 x 306)  idle
        axial-0.0-115.0-PCA-07 (1 x 306)  idle
        axial-0.0-115.0-PCA-08 (1 x 306)  idle
    Range : 47000 ... 437999 =     47.000 ...   437.999 secs
Current compensation grade : 0
Reading 0 ... 10000  =      0.000 ...    10.000 secs...
Setting up low-pass filter at 40 Hz
h_trans_bandwidth chosen to be 10.0 Hz
Filter length of 331 samples (0.331 sec) selected

Now we can clean the data with OTP, lowpass, and plot. The flux jumps have been suppressed alongside the random sensor noise.

raw_clean = mne.preprocessing.oversampled_temporal_projection(raw)
raw_clean.filter(0., 40.)
raw_clean.plot(order=order, n_channels=10)


Processing MEG data using oversampled temporal projection
    Processing 1 data chunk of (at least) 10.0 sec with hann windowing
    The final 0.001 sec will be lumped into the final window
    Denoising     0.00 -    10.00 sec
Setting up low-pass filter at 40 Hz
h_trans_bandwidth chosen to be 10.0 Hz
Filter length of 331 samples (0.331 sec) selected

We can also look at the effect on single-trial phantom localization. See the Brainstorm Elekta phantom dataset tutorial for more information. Here we use a version that does single-trial localization across the 17 trials are in our 10-second window:

def compute_bias(raw):
    events = find_events(raw, 'STI201', verbose=False)
    events = events[1:]  # first one has an artifact
    tmin, tmax = -0.2, 0.1
    epochs = mne.Epochs(raw, events, dipole_number, tmin, tmax,
                        baseline=(None, -0.01), preload=True, verbose=False)
    sphere = mne.make_sphere_model(r0=(0., 0., 0.), head_radius=None,
    cov = mne.compute_covariance(epochs, tmax=0, method='shrunk',
    idx = epochs.time_as_index(0.036)[0]
    data = epochs.get_data()[:, :, idx].T
    evoked = mne.EvokedArray(data,, tmin=0.)
    dip = fit_dipole(evoked, cov, sphere, n_jobs=1, verbose=False)[0]
    actual_pos = mne.dipole.get_phantom_dipoles()[0][dipole_number - 1]
    misses = 1000 * np.linalg.norm(dip.pos - actual_pos, axis=-1)
    return misses

bias = compute_bias(raw)
print('Raw bias: %0.1fmm (worst: %0.1fmm)'
      % (np.mean(bias), np.max(bias)))
bias_clean = compute_bias(raw_clean)
print('OTP bias: %0.1fmm (worst: %0.1fmm)'
      % (np.mean(bias_clean), np.max(bias_clean),))


Raw bias: 2.5mm (worst: 5.2mm)
OTP bias: 1.3mm (worst: 1.5mm)


[1]Larson E, Taulu S (2017). Reducing Sensor Noise in MEG and EEG Recordings Using Oversampled Temporal Projection. IEEE Transactions on Biomedical Engineering.

Total running time of the script: ( 0 minutes 43.820 seconds)

Gallery generated by Sphinx-Gallery