Artifact Correction with SSP

import numpy as np

import mne
from mne.datasets import sample
from mne.preprocessing import compute_proj_ecg, compute_proj_eog

# getting some data ready
data_path = sample.data_path()
raw_fname = data_path + '/MEG/sample/sample_audvis_filt-0-40_raw.fif'

raw = mne.io.read_raw_fif(raw_fname, preload=True)
raw.set_eeg_reference()
raw.pick_types(meg=True, ecg=True, eog=True, stim=True)

Out:

Opening raw data file /home/ubuntu/mne_data/MNE-sample-data/MEG/sample/sample_audvis_filt-0-40_raw.fif...
    Read a total of 4 projection items:
        PCA-v1 (1 x 102)  idle
        PCA-v2 (1 x 102)  idle
        PCA-v3 (1 x 102)  idle
        Average EEG reference (1 x 60)  idle
    Range : 6450 ... 48149 =     42.956 ...   320.665 secs
Ready.
Current compensation grade : 0
Reading 0 ... 41699  =      0.000 ...   277.709 secs...
An average reference projection was already added. The data has been left untouched.

Compute SSP projections

projs, events = compute_proj_ecg(raw, n_grad=1, n_mag=1, average=True)
print(projs)

ecg_projs = projs[-2:]
mne.viz.plot_projs_topomap(ecg_projs)

# Now for EOG

projs, events = compute_proj_eog(raw, n_grad=1, n_mag=1, average=True)
print(projs)

eog_projs = projs[-2:]
mne.viz.plot_projs_topomap(eog_projs)
  • ../_images/sphx_glr_plot_artifacts_correction_ssp_001.png
  • ../_images/sphx_glr_plot_artifacts_correction_ssp_002.png

Out:

Including 4 SSP projectors from raw file
Running ECG SSP computation
Reconstructing ECG signal from Magnetometers
Setting up band-pass filter from 5 - 35 Hz
Number of ECG events detected : 285 (average pulse 61 / min.)
Computing projector
Setting up band-pass filter from 1 - 35 Hz
285 matching events found
Created an SSP operator (subspace dimension = 3)
4 projection items activated
Loading data for 285 events and 91 original time points ...
    Rejecting  epoch based on MAG : [u'MEG 1421']
    Rejecting  epoch based on MAG : [u'MEG 1411', u'MEG 1421']
    Rejecting  epoch based on MAG : [u'MEG 1711']
    Rejecting  epoch based on EOG : [u'EOG 061']
    Rejecting  epoch based on EOG : [u'EOG 061']
    Rejecting  epoch based on MAG : [u'MEG 1711']
    Rejecting  epoch based on MAG : [u'MEG 1411', u'MEG 1421']
    Rejecting  epoch based on MAG : [u'MEG 1411']
    Rejecting  epoch based on MAG : [u'MEG 1421']
    Rejecting  epoch based on EOG : [u'EOG 061']
    Rejecting  epoch based on EOG : [u'EOG 061']
    Rejecting  epoch based on EOG : [u'EOG 061']
    Rejecting  epoch based on MAG : [u'MEG 1411']
    Rejecting  epoch based on MAG : [u'MEG 1441']
    Rejecting  epoch based on EOG : [u'EOG 061']
15 bad epochs dropped
No EEG channels found. Forcing n_eeg to 0
Adding projection: planar--0.200-0.400-PCA-01
Adding projection: axial--0.200-0.400-PCA-01
Done.
[<Projection  |  PCA-v1, active : False, n_channels : 102>, <Projection  |  PCA-v2, active : False, n_channels : 102>, <Projection  |  PCA-v3, active : False, n_channels : 102>, <Projection  |  Average EEG reference, active : False, n_channels : 60>, <Projection  |  ECG-planar--0.200-0.400-PCA-01, active : False, n_channels : 203>, <Projection  |  ECG-axial--0.200-0.400-PCA-01, active : False, n_channels : 102>]
Including 4 SSP projectors from raw file
Running EOG SSP computation
EOG channel index for this subject is: [314]
Filtering the data to remove DC offset to help distinguish blinks from saccades
Setting up band-pass filter from 2 - 45 Hz
Setting up band-pass filter from 1 - 10 Hz
Now detecting blinks and generating corresponding events
Number of EOG events detected : 46
Computing projector
Setting up band-pass filter from 1 - 35 Hz
46 matching events found
Created an SSP operator (subspace dimension = 3)
4 projection items activated
Loading data for 46 events and 61 original time points ...
    Rejecting  epoch based on MAG : [u'MEG 1421']
    Rejecting  epoch based on MAG : [u'MEG 1411', u'MEG 1421']
    Rejecting  epoch based on MAG : [u'MEG 1411', u'MEG 1421']
    Rejecting  epoch based on MAG : [u'MEG 1411']
    Rejecting  epoch based on MAG : [u'MEG 1411', u'MEG 1421']
5 bad epochs dropped
No EEG channels found. Forcing n_eeg to 0
Adding projection: planar--0.200-0.200-PCA-01
Adding projection: axial--0.200-0.200-PCA-01
Done.
[<Projection  |  PCA-v1, active : False, n_channels : 102>, <Projection  |  PCA-v2, active : False, n_channels : 102>, <Projection  |  PCA-v3, active : False, n_channels : 102>, <Projection  |  Average EEG reference, active : False, n_channels : 60>, <Projection  |  EOG-planar--0.200-0.200-PCA-01, active : False, n_channels : 203>, <Projection  |  EOG-axial--0.200-0.200-PCA-01, active : False, n_channels : 102>]

Apply SSP projections

MNE is handling projections at the level of the info, so to register them populate the list that you find in the ‘proj’ field

raw.info['projs'] += eog_projs + ecg_projs

Yes this was it. Now MNE will apply the projs on demand at any later stage, so watch out for proj parmeters in functions or to it explicitly with the .apply_proj method

Demonstrate SSP cleaning on some evoked data

events = mne.find_events(raw, stim_channel='STI 014')
reject = dict(grad=4000e-13, mag=4e-12, eog=150e-6)
# this can be highly data dependent
event_id = {'auditory/left': 1}

epochs_no_proj = mne.Epochs(raw, events, event_id, tmin=-0.2, tmax=0.5,
                            proj=False, baseline=(None, 0), reject=reject)
epochs_no_proj.average().plot(spatial_colors=True)


epochs_proj = mne.Epochs(raw, events, event_id, tmin=-0.2, tmax=0.5, proj=True,
                         baseline=(None, 0), reject=reject)
epochs_proj.average().plot(spatial_colors=True)
  • ../_images/sphx_glr_plot_artifacts_correction_ssp_003.png
  • ../_images/sphx_glr_plot_artifacts_correction_ssp_004.png

Out:

319 events found
Events id: [ 1  2  3  4  5 32]
72 matching events found
Created an SSP operator (subspace dimension = 7)
    Rejecting  epoch based on EOG : [u'EOG 061']
    Rejecting  epoch based on EOG : [u'EOG 061']
    Rejecting  epoch based on EOG : [u'EOG 061']
    Rejecting  epoch based on EOG : [u'EOG 061']
    Rejecting  epoch based on EOG : [u'EOG 061']
    Rejecting  epoch based on MAG : [u'MEG 1711']
    Rejecting  epoch based on EOG : [u'EOG 061']
    Rejecting  epoch based on EOG : [u'EOG 061']
    Rejecting  epoch based on EOG : [u'EOG 061']
    Rejecting  epoch based on EOG : [u'EOG 061']
    Rejecting  epoch based on EOG : [u'EOG 061']
    Rejecting  epoch based on EOG : [u'EOG 061']
    Rejecting  epoch based on EOG : [u'EOG 061']
    Rejecting  epoch based on EOG : [u'EOG 061']
    Rejecting  epoch based on EOG : [u'EOG 061']
    Rejecting  epoch based on EOG : [u'EOG 061']
    Rejecting  epoch based on EOG : [u'EOG 061']
72 matching events found
Created an SSP operator (subspace dimension = 7)
8 projection items activated
    Rejecting  epoch based on EOG : [u'EOG 061']
    Rejecting  epoch based on EOG : [u'EOG 061']
    Rejecting  epoch based on EOG : [u'EOG 061']
    Rejecting  epoch based on EOG : [u'EOG 061']
    Rejecting  epoch based on EOG : [u'EOG 061']
    Rejecting  epoch based on MAG : [u'MEG 1711']
    Rejecting  epoch based on EOG : [u'EOG 061']
    Rejecting  epoch based on EOG : [u'EOG 061']
    Rejecting  epoch based on EOG : [u'EOG 061']
    Rejecting  epoch based on EOG : [u'EOG 061']
    Rejecting  epoch based on EOG : [u'EOG 061']
    Rejecting  epoch based on EOG : [u'EOG 061']
    Rejecting  epoch based on EOG : [u'EOG 061']
    Rejecting  epoch based on EOG : [u'EOG 061']
    Rejecting  epoch based on EOG : [u'EOG 061']
    Rejecting  epoch based on EOG : [u'EOG 061']
    Rejecting  epoch based on EOG : [u'EOG 061']

Looks cool right? It is however often not clear how many components you should take and unfortunately this can have bad consequences as can be seen interactively using the delayed SSP mode:

evoked = mne.Epochs(raw, events, event_id, tmin=-0.2, tmax=0.5,
                    proj='delayed', baseline=(None, 0),
                    reject=reject).average()

# set time instants in seconds (from 50 to 150ms in a step of 10ms)
times = np.arange(0.05, 0.15, 0.01)

evoked.plot_topomap(times, proj='interactive')
  • ../_images/sphx_glr_plot_artifacts_correction_ssp_005.png
  • ../_images/sphx_glr_plot_artifacts_correction_ssp_006.png

Out:

72 matching events found
Entering delayed SSP mode.
Created an SSP operator (subspace dimension = 7)
    Rejecting  epoch based on EOG : [u'EOG 061']
    Rejecting  epoch based on EOG : [u'EOG 061']
    Rejecting  epoch based on EOG : [u'EOG 061']
    Rejecting  epoch based on EOG : [u'EOG 061']
    Rejecting  epoch based on EOG : [u'EOG 061']
    Rejecting  epoch based on MAG : [u'MEG 1711']
    Rejecting  epoch based on EOG : [u'EOG 061']
    Rejecting  epoch based on EOG : [u'EOG 061']
    Rejecting  epoch based on EOG : [u'EOG 061']
    Rejecting  epoch based on EOG : [u'EOG 061']
    Rejecting  epoch based on EOG : [u'EOG 061']
    Rejecting  epoch based on EOG : [u'EOG 061']
    Rejecting  epoch based on EOG : [u'EOG 061']
    Rejecting  epoch based on EOG : [u'EOG 061']
    Rejecting  epoch based on EOG : [u'EOG 061']
    Rejecting  epoch based on EOG : [u'EOG 061']
    Rejecting  epoch based on EOG : [u'EOG 061']

now you should see checkboxes. Remove a few SSP and see how the auditory pattern suddenly drops off

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

Download Python source code: plot_artifacts_correction_ssp.py
Download Jupyter notebook: plot_artifacts_correction_ssp.ipynb

Generated by Sphinx-Gallery