Compute ICA on MEG data and remove artifactsΒΆ

ICA is fit to MEG raw data. The sources matching the ECG and EOG are automatically found and displayed. Subsequently, artifact detection and rejection quality are assessed.

# Authors: Denis Engemann <denis.engemann@gmail.com>
#          Alexandre Gramfort <alexandre.gramfort@telecom-paristech.fr>
#
# License: BSD (3-clause)

import numpy as np

import mne
from mne.preprocessing import ICA
from mne.preprocessing import create_ecg_epochs, create_eog_epochs
from mne.datasets import sample

Setup paths and prepare raw data

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, add_eeg_ref=False)
raw.filter(1, 45, n_jobs=1, l_trans_bandwidth=0.5, h_trans_bandwidth=0.5,
           filter_length='10s', phase='zero-double')

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...
Band-pass filtering from 1 - 45 Hz
fir_window in 0.13 is "hann" but will change to "hamming" in 0.14
  1. Fit ICA model using the FastICA algorithm
# Other available choices are `infomax` or `extended-infomax`
# We pass a float value between 0 and 1 to select n_components based on the
# percentage of variance explained by the PCA components.

ica = ICA(n_components=0.95, method='fastica')

picks = mne.pick_types(raw.info, meg=True, eeg=False, eog=False,
                       stim=False, exclude='bads')

ica.fit(raw, picks=picks, decim=3, reject=dict(mag=4e-12, grad=4000e-13))

# maximum number of components to reject
n_max_ecg, n_max_eog = 3, 1  # here we don't expect horizontal EOG components

Out:

Fitting ICA to data using 305 channels.
Please be patient, this may take some time
Inferring max_pca_components from picks.
    Rejecting  epoch based on MAG : [u'MEG 1711']
Artifact detected in [4242, 4343]
    Rejecting  epoch based on MAG : [u'MEG 1711']
Artifact detected in [5858, 5959]
    Rejecting  epoch based on MAG : [u'MEG 1411']
Artifact detected in [13433, 13534]
Selection by explained variance: 123 components
  1. identify bad components by analyzing latent sources.
title = 'Sources related to %s artifacts (red)'

# generate ECG epochs use detection via phase statistics

ecg_epochs = create_ecg_epochs(raw, tmin=-.5, tmax=.5, picks=picks)

ecg_inds, scores = ica.find_bads_ecg(ecg_epochs, method='ctps')
ica.plot_scores(scores, exclude=ecg_inds, title=title % 'ecg', labels='ecg')

show_picks = np.abs(scores).argsort()[::-1][:5]

ica.plot_sources(raw, show_picks, exclude=ecg_inds, title=title % 'ecg')
ica.plot_components(ecg_inds, title=title % 'ecg', colorbar=True)

ecg_inds = ecg_inds[:n_max_ecg]
ica.exclude += ecg_inds

# detect EOG by correlation

eog_inds, scores = ica.find_bads_eog(raw)
ica.plot_scores(scores, exclude=eog_inds, title=title % 'eog', labels='eog')

show_picks = np.abs(scores).argsort()[::-1][:5]

ica.plot_sources(raw, show_picks, exclude=eog_inds, title=title % 'eog')
ica.plot_components(eog_inds, title=title % 'eog', colorbar=True)

eog_inds = eog_inds[:n_max_eog]
ica.exclude += eog_inds
  • ../_images/sphx_glr_plot_ica_from_raw_001.png
  • ../_images/sphx_glr_plot_ica_from_raw_002.png
  • ../_images/sphx_glr_plot_ica_from_raw_003.png
  • ../_images/sphx_glr_plot_ica_from_raw_004.png
  • ../_images/sphx_glr_plot_ica_from_raw_005.png
  • ../_images/sphx_glr_plot_ica_from_raw_006.png

Out:

Reconstructing ECG signal from Magnetometers
Band-pass filtering from 8 - 16 Hz
Number of ECG events detected : 284 (average pulse 61 / min.)
Creating RawArray with float64 data, n_channels=1, n_times=41700
    Range : 0 ... 41699 =      0.000 ...   277.709 secs
Ready.
284 matching events found
No baseline correction applied
Created an SSP operator (subspace dimension = 3)
Loading data for 284 events and 151 original time points ...
0 bad epochs dropped
Reconstructing ECG signal from Magnetometers
... filtering ICA sources
Band-pass filtering from 1 - 10 Hz
Multiple deprecated filter parameters were used:
phase in 0.13 is "zero-double" but will change to "zero" in 0.14
fir_window in 0.13 is "hann" but will change to "hamming" in 0.14
lower transition bandwidth in 0.13 is 0.5 Hz but will change to "auto" in 0.14
upper transition bandwidth in 0.13 is 0.5 Hz but will change to "auto" in 0.14
The default filter length in 0.13 is "10s" but will change to "auto" in 0.14
... filtering target
Band-pass filtering from 1 - 10 Hz
Multiple deprecated filter parameters were used:
phase in 0.13 is "zero-double" but will change to "zero" in 0.14
fir_window in 0.13 is "hann" but will change to "hamming" in 0.14
lower transition bandwidth in 0.13 is 0.5 Hz but will change to "auto" in 0.14
upper transition bandwidth in 0.13 is 0.5 Hz but will change to "auto" in 0.14
The default filter length in 0.13 is "10s" but will change to "auto" in 0.14
  1. Assess component selection and unmixing quality
# estimate average artifact
ecg_evoked = ecg_epochs.average()
ica.plot_sources(ecg_evoked, exclude=ecg_inds)  # plot ECG sources + selection
ica.plot_overlay(ecg_evoked, exclude=ecg_inds)  # plot ECG cleaning

eog_evoked = create_eog_epochs(raw, tmin=-.5, tmax=.5, picks=picks).average()
ica.plot_sources(eog_evoked, exclude=eog_inds)  # plot EOG sources + selection
ica.plot_overlay(eog_evoked, exclude=eog_inds)  # plot EOG cleaning

# check the amplitudes do not change
ica.plot_overlay(raw)  # EOG artifacts remain
  • ../_images/sphx_glr_plot_ica_from_raw_007.png
  • ../_images/sphx_glr_plot_ica_from_raw_008.png
  • ../_images/sphx_glr_plot_ica_from_raw_009.png
  • ../_images/sphx_glr_plot_ica_from_raw_010.png
  • ../_images/sphx_glr_plot_ica_from_raw_011.png

Out:

Transforming to ICA space (123 components)
Zeroing out 3 ICA components
EOG channel index for this subject is: [375]
Filtering the data to remove DC offset to help distinguish blinks from saccades
Band-pass filtering from 2 - 45 Hz
Band-pass filtering from 1 - 10 Hz
Now detecting blinks and generating corresponding events
Number of EOG events detected : 46
46 matching events found
No baseline correction applied
Created an SSP operator (subspace dimension = 3)
Loading data for 46 events and 151 original time points ...
0 bad epochs dropped
Transforming to ICA space (123 components)
Zeroing out 3 ICA components
Transforming to ICA space (123 components)
Zeroing out 3 ICA components
# To save an ICA solution you can say:
# ica.save('my_ica.fif')

# You can later load the solution by saying:
# from mne.preprocessing import read_ica
# read_ica('my_ica.fif')

# Apply the solution to Raw, Epochs or Evoked like this:
# ica.apply(epochs)

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

Generated by Sphinx-Gallery