Compute ICA components on epochs

ICA is fit to MEG raw data. We assume that the non-stationary EOG artifacts have already been removed. The sources matching the ECG are automatically found and displayed.

Note

This example does quite a bit of processing, so even on a fast machine it can take about a minute to complete.

# Authors: Denis Engemann <denis.engemann@gmail.com>
#
# License: BSD (3-clause)

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

print(__doc__)

Read and preprocess the data. Preprocessing consists of:

  • MEG channel selection

  • 1-30 Hz band-pass filter

  • epoching -0.2 to 0.5 seconds with respect to events

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, eeg=False, exclude='bads', stim=True)
raw.filter(1, 30, fir_design='firwin')

# longer + more epochs for more artifact exposure
events = mne.find_events(raw, stim_channel='STI 014')
epochs = mne.Epochs(raw, events, event_id=None, tmin=-0.2, tmax=0.5)

Out:

Opening raw data file /home/circleci/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...
Filtering raw data in 1 contiguous segment
Setting up band-pass filter from 1 - 30 Hz

FIR filter parameters
---------------------
Designing a one-pass, zero-phase, non-causal bandpass filter:
- Windowed time-domain design (firwin) method
- Hamming window with 0.0194 passband ripple and 53 dB stopband attenuation
- Lower passband edge: 1.00
- Lower transition bandwidth: 1.00 Hz (-6 dB cutoff frequency: 0.50 Hz)
- Upper passband edge: 30.00 Hz
- Upper transition bandwidth: 7.50 Hz (-6 dB cutoff frequency: 33.75 Hz)
- Filter length: 497 samples (3.310 sec)

319 events found
Event IDs: [ 1  2  3  4  5 32]
319 matching events found
Applying baseline correction (mode: mean)
Not setting metadata
Created an SSP operator (subspace dimension = 3)
4 projection items activated

Fit ICA model using the FastICA algorithm, detect and plot components explaining ECG artifacts.

ica = ICA(n_components=0.95, method='fastica').fit(epochs)

ecg_epochs = create_ecg_epochs(raw, tmin=-.5, tmax=.5)
ecg_inds, scores = ica.find_bads_ecg(ecg_epochs)

ica.plot_components(ecg_inds)
../../_images/sphx_glr_plot_run_ica_001.png

Out:

Fitting ICA to data using 305 channels (please be patient, this may take a while)
Inferring max_pca_components from picks
Loading data for 319 events and 106 original time points ...
0 bad epochs dropped
Selection by explained variance: 126 components
/home/circleci/.local/lib/python3.7/site-packages/sklearn/decomposition/fastica_.py:119: ConvergenceWarning: FastICA did not converge. Consider increasing tolerance or the maximum number of iterations.
  ConvergenceWarning)
Loading data for 319 events and 106 original time points ...
Fitting ICA took 30.0s.
Reconstructing ECG signal from Magnetometers
Setting up band-pass filter from 8 - 16 Hz

FIR filter parameters
---------------------
Designing a two-pass forward and reverse, zero-phase, non-causal bandpass filter:
- Windowed frequency-domain design (firwin2) method
- Hann window
- Lower passband edge: 8.00
- Lower transition bandwidth: 0.50 Hz (-12 dB cutoff frequency: 7.75 Hz)
- Upper passband edge: 16.00 Hz
- Upper transition bandwidth: 0.50 Hz (-12 dB cutoff frequency: 16.25 Hz)
- Filter length: 2048 samples (13.639 sec)

Number of ECG events detected : 284 (average pulse 61 / min.)
284 matching events found
No baseline correction applied
Not setting metadata
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

Plot properties of ECG components:

ica.plot_properties(epochs, picks=ecg_inds)
  • ../../_images/sphx_glr_plot_run_ica_002.png
  • ../../_images/sphx_glr_plot_run_ica_003.png

Out:

Loading data for 319 events and 106 original time points ...
Loading data for 319 events and 106 original time points ...
    Using multitaper spectrum estimation with 7 DPSS windows
319 matching events found
No baseline correction applied
Not setting metadata
0 projection items activated
0 bad epochs dropped
319 matching events found
No baseline correction applied
Not setting metadata
0 projection items activated
0 bad epochs dropped

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

Estimated memory usage: 358 MB

Gallery generated by Sphinx-Gallery