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.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
Current compensation grade : 0
    Range : 6450 ... 48149 =     42.956 ...   320.665 secs
Ready.
Reading 0 ... 41699  =      0.000 ...   277.709 secs...

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
Number of ECG events detected : 285 (average pulse 61 / min.)
Computing projector
Filtering a subset of channels. The highpass and lowpass values in the measurement info will not be updated.
Band-pass filtering from 1 - 35 Hz
285 matching events found
No baseline correction applied
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
Now detecting blinks and generating corresponding events
Number of EOG events detected : 46
Computing projector
Filtering a subset of channels. The highpass and lowpass values in the measurement info will not be updated.
Band-pass filtering from 1 - 35 Hz
46 matching events found
No baseline correction applied
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
Applying baseline correction (mode: mean)
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
Applying baseline correction (mode: mean)
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
Applying baseline correction (mode: mean)
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 16.008 seconds)

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