{
  "nbformat_minor": 0, 
  "nbformat": 4, 
  "cells": [
    {
      "execution_count": null, 
      "cell_type": "code", 
      "source": [
        "%matplotlib inline"
      ], 
      "outputs": [], 
      "metadata": {
        "collapsed": false
      }
    }, 
    {
      "source": [
        "\n\nCompute ICA on MEG data and remove artifacts\n============================================\n\nICA is fit to MEG raw data.\nThe sources matching the ECG and EOG are automatically found and displayed.\nSubsequently, artifact detection and rejection quality are assessed.\n\n"
      ], 
      "cell_type": "markdown", 
      "metadata": {}
    }, 
    {
      "execution_count": null, 
      "cell_type": "code", 
      "source": [
        "# Authors: Denis Engemann <denis.engemann@gmail.com>\n#          Alexandre Gramfort <alexandre.gramfort@telecom-paristech.fr>\n#\n# License: BSD (3-clause)\n\nimport numpy as np\n\nimport mne\nfrom mne.preprocessing import ICA\nfrom mne.preprocessing import create_ecg_epochs, create_eog_epochs\nfrom mne.datasets import sample"
      ], 
      "outputs": [], 
      "metadata": {
        "collapsed": false
      }
    }, 
    {
      "source": [
        "Setup paths and prepare raw data.\n\n"
      ], 
      "cell_type": "markdown", 
      "metadata": {}
    }, 
    {
      "execution_count": null, 
      "cell_type": "code", 
      "source": [
        "data_path = sample.data_path()\nraw_fname = data_path + '/MEG/sample/sample_audvis_filt-0-40_raw.fif'\n\nraw = mne.io.read_raw_fif(raw_fname, preload=True)\nraw.filter(1, 45, n_jobs=1, l_trans_bandwidth=0.5, h_trans_bandwidth=0.5,\n           filter_length='10s', phase='zero-double')\nraw.annotations = mne.Annotations([1], [10], 'BAD')\nraw.plot(block=True)\n\n# For the sake of example we annotate first 10 seconds of the recording as\n# 'BAD'. This part of data is excluded from the ICA decomposition by default.\n# To turn this behavior off, pass ``reject_by_annotation=False`` to\n# :meth:`mne.preprocessing.ICA.fit`.\nraw.annotations = mne.Annotations([0], [10], 'BAD')"
      ], 
      "outputs": [], 
      "metadata": {
        "collapsed": false
      }
    }, 
    {
      "source": [
        "1) Fit ICA model using the FastICA algorithm.\n\n"
      ], 
      "cell_type": "markdown", 
      "metadata": {}
    }, 
    {
      "execution_count": null, 
      "cell_type": "code", 
      "source": [
        "# Other available choices are `infomax` or `extended-infomax`\n# We pass a float value between 0 and 1 to select n_components based on the\n# percentage of variance explained by the PCA components.\n\nica = ICA(n_components=0.95, method='fastica')\n\npicks = mne.pick_types(raw.info, meg=True, eeg=False, eog=False,\n                       stim=False, exclude='bads')\n\nica.fit(raw, picks=picks, decim=3, reject=dict(mag=4e-12, grad=4000e-13))\n\n# maximum number of components to reject\nn_max_ecg, n_max_eog = 3, 1  # here we don't expect horizontal EOG components"
      ], 
      "outputs": [], 
      "metadata": {
        "collapsed": false
      }
    }, 
    {
      "source": [
        "2) identify bad components by analyzing latent sources.\n\n"
      ], 
      "cell_type": "markdown", 
      "metadata": {}
    }, 
    {
      "execution_count": null, 
      "cell_type": "code", 
      "source": [
        "title = 'Sources related to %s artifacts (red)'\n\n# generate ECG epochs use detection via phase statistics\n\necg_epochs = create_ecg_epochs(raw, tmin=-.5, tmax=.5, picks=picks)\n\necg_inds, scores = ica.find_bads_ecg(ecg_epochs, method='ctps')\nica.plot_scores(scores, exclude=ecg_inds, title=title % 'ecg', labels='ecg')\n\nshow_picks = np.abs(scores).argsort()[::-1][:5]\n\nica.plot_sources(raw, show_picks, exclude=ecg_inds, title=title % 'ecg')\nica.plot_components(ecg_inds, title=title % 'ecg', colorbar=True)\n\necg_inds = ecg_inds[:n_max_ecg]\nica.exclude += ecg_inds\n\n# detect EOG by correlation\n\neog_inds, scores = ica.find_bads_eog(raw)\nica.plot_scores(scores, exclude=eog_inds, title=title % 'eog', labels='eog')\n\nshow_picks = np.abs(scores).argsort()[::-1][:5]\n\nica.plot_sources(raw, show_picks, exclude=eog_inds, title=title % 'eog')\nica.plot_components(eog_inds, title=title % 'eog', colorbar=True)\n\neog_inds = eog_inds[:n_max_eog]\nica.exclude += eog_inds"
      ], 
      "outputs": [], 
      "metadata": {
        "collapsed": false
      }
    }, 
    {
      "source": [
        "3) Assess component selection and unmixing quality.\n\n"
      ], 
      "cell_type": "markdown", 
      "metadata": {}
    }, 
    {
      "execution_count": null, 
      "cell_type": "code", 
      "source": [
        "# estimate average artifact\necg_evoked = ecg_epochs.average()\nica.plot_sources(ecg_evoked, exclude=ecg_inds)  # plot ECG sources + selection\nica.plot_overlay(ecg_evoked, exclude=ecg_inds)  # plot ECG cleaning\n\neog_evoked = create_eog_epochs(raw, tmin=-.5, tmax=.5, picks=picks).average()\nica.plot_sources(eog_evoked, exclude=eog_inds)  # plot EOG sources + selection\nica.plot_overlay(eog_evoked, exclude=eog_inds)  # plot EOG cleaning\n\n# check the amplitudes do not change\nica.plot_overlay(raw)  # EOG artifacts remain"
      ], 
      "outputs": [], 
      "metadata": {
        "collapsed": false
      }
    }, 
    {
      "execution_count": null, 
      "cell_type": "code", 
      "source": [
        "# To save an ICA solution you can say:\n# ica.save('my_ica.fif')\n\n# You can later load the solution by saying:\n# from mne.preprocessing import read_ica\n# read_ica('my_ica.fif')\n\n# Apply the solution to Raw, Epochs or Evoked like this:\n# ica.apply(epochs)"
      ], 
      "outputs": [], 
      "metadata": {
        "collapsed": false
      }
    }
  ], 
  "metadata": {
    "kernelspec": {
      "display_name": "Python 2", 
      "name": "python2", 
      "language": "python"
    }, 
    "language_info": {
      "mimetype": "text/x-python", 
      "nbconvert_exporter": "python", 
      "name": "python", 
      "file_extension": ".py", 
      "version": "2.7.13", 
      "pygments_lexer": "ipython2", 
      "codemirror_mode": {
        "version": 2, 
        "name": "ipython"
      }
    }
  }
}