{ "nbformat_minor": 0, "nbformat": 4, "cells": [ { "execution_count": null, "cell_type": "code", "source": [ "%matplotlib inline" ], "outputs": [], "metadata": { "collapsed": false } }, { "source": [ "\n\n\nArtifact Correction with ICA\n============================\n\nICA finds directions in the feature space\ncorresponding to projections with high non-Gaussianity. We thus obtain\na decomposition into independent components, and the artifact's contribution\nis localized in only a small number of components.\nThese components have to be correctly identified and removed.\n\nIf EOG or ECG recordings are available, they can be used in ICA to\nautomatically select the corresponding artifact components from the\ndecomposition. To do so, you have to first build an Epoch object around\nblink or heartbeat event.\n\n" ], "cell_type": "markdown", "metadata": {} }, { "execution_count": null, "cell_type": "code", "source": [ "import numpy as np\n\nimport mne\nfrom mne.datasets import sample\n\nfrom mne.preprocessing import ICA\nfrom mne.preprocessing import create_eog_epochs, create_ecg_epochs\n\n# getting some data ready\ndata_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, 40, n_jobs=2) # 1Hz high pass is often helpful for fitting ICA\n\npicks_meg = mne.pick_types(raw.info, meg=True, eeg=False, eog=False,\n stim=False, exclude='bads')" ], "outputs": [], "metadata": { "collapsed": false } }, { "source": [ "Before applying artifact correction please learn about your actual artifacts\nby reading `tut_artifacts_detect`.\n\n" ], "cell_type": "markdown", "metadata": {} }, { "source": [ "Fit ICA\n-------\n\nICA parameters:\n\n" ], "cell_type": "markdown", "metadata": {} }, { "execution_count": null, "cell_type": "code", "source": [ "n_components = 25 # if float, select n_components by explained variance of PCA\nmethod = 'fastica' # for comparison with EEGLAB try \"extended-infomax\" here\ndecim = 3 # we need sufficient statistics, not all time points -> saves time\n\n# we will also set state of the random number generator - ICA is a\n# non-deterministic algorithm, but we want to have the same decomposition\n# and the same order of components each time this tutorial is run\nrandom_state = 23" ], "outputs": [], "metadata": { "collapsed": false } }, { "source": [ "Define the ICA object instance\n\n" ], "cell_type": "markdown", "metadata": {} }, { "execution_count": null, "cell_type": "code", "source": [ "ica = ICA(n_components=n_components, method=method, random_state=random_state)\nprint(ica)" ], "outputs": [], "metadata": { "collapsed": false } }, { "source": [ "we avoid fitting ICA on crazy environmental artifacts that would\ndominate the variance and decomposition\n\n" ], "cell_type": "markdown", "metadata": {} }, { "execution_count": null, "cell_type": "code", "source": [ "reject = dict(mag=5e-12, grad=4000e-13)\nica.fit(raw, picks=picks_meg, decim=decim, reject=reject)\nprint(ica)" ], "outputs": [], "metadata": { "collapsed": false } }, { "source": [ "Plot ICA components\n\n" ], "cell_type": "markdown", "metadata": {} }, { "execution_count": null, "cell_type": "code", "source": [ "ica.plot_components() # can you spot some potential bad guys?" ], "outputs": [], "metadata": { "collapsed": false } }, { "source": [ "Component properties\n--------------------\n\nLet's take a closer look at properties of first three independent components.\n\n" ], "cell_type": "markdown", "metadata": {} }, { "execution_count": null, "cell_type": "code", "source": [ "# first, component 0:\nica.plot_properties(raw, picks=0)" ], "outputs": [], "metadata": { "collapsed": false } }, { "source": [ "we can see that the data were filtered so the spectrum plot is not\nvery informative, let's change that:\n\n" ], "cell_type": "markdown", "metadata": {} }, { "execution_count": null, "cell_type": "code", "source": [ "ica.plot_properties(raw, picks=0, psd_args={'fmax': 35.})" ], "outputs": [], "metadata": { "collapsed": false } }, { "source": [ "we can also take a look at multiple different components at once:\n\n" ], "cell_type": "markdown", "metadata": {} }, { "execution_count": null, "cell_type": "code", "source": [ "ica.plot_properties(raw, picks=[1, 2], psd_args={'fmax': 35.})" ], "outputs": [], "metadata": { "collapsed": false } }, { "source": [ "Instead of opening individual figures with component properties, we can\nalso pass an instance of Raw or Epochs in ``inst`` arument to\n``ica.plot_components``. This would allow us to open component properties\ninteractively by clicking on individual component topomaps. In the notebook\nthis woks only when running matplotlib in interactive mode (``%matplotlib``).\n\n" ], "cell_type": "markdown", "metadata": {} }, { "execution_count": null, "cell_type": "code", "source": [ "# uncomment the code below to test the inteactive mode of plot_components:\n# ica.plot_components(picks=range(10), inst=raw)" ], "outputs": [], "metadata": { "collapsed": false } }, { "source": [ "Advanced artifact detection\n---------------------------\n\nLet's use a more efficient way to find artefacts\n\n" ], "cell_type": "markdown", "metadata": {} }, { "execution_count": null, "cell_type": "code", "source": [ "eog_average = create_eog_epochs(raw, reject=dict(mag=5e-12, grad=4000e-13),\n picks=picks_meg).average()\n\n# We simplify things by setting the maximum number of components to reject\nn_max_eog = 1 # here we bet on finding the vertical EOG components\neog_epochs = create_eog_epochs(raw, reject=reject) # get single EOG trials\neog_inds, scores = ica.find_bads_eog(eog_epochs) # find via correlation\n\nica.plot_scores(scores, exclude=eog_inds) # look at r scores of components\n# we can see that only one component is highly correlated and that this\n# component got detected by our correlation analysis (red).\n\nica.plot_sources(eog_average, exclude=eog_inds) # look at source time course" ], "outputs": [], "metadata": { "collapsed": false } }, { "source": [ "We can take a look at the properties of that component, now using the\ndata epoched with respect to EOG events.\nWe will also use a little bit of smoothing along the trials axis in the\nepochs image:\n\n" ], "cell_type": "markdown", "metadata": {} }, { "execution_count": null, "cell_type": "code", "source": [ "ica.plot_properties(eog_epochs, picks=eog_inds, psd_args={'fmax': 35.},\n image_args={'sigma': 1.})" ], "outputs": [], "metadata": { "collapsed": false } }, { "source": [ "That component is showing a prototypical average vertical EOG time course.\n\nPay attention to the labels, a customized read-out of the\n``mne.preprocessing.ICA.labels_``:\n\n" ], "cell_type": "markdown", "metadata": {} }, { "execution_count": null, "cell_type": "code", "source": [ "print(ica.labels_)" ], "outputs": [], "metadata": { "collapsed": false } }, { "source": [ "These labels were used by the plotters and are added automatically\nby artifact detection functions. You can also manually edit them to annotate\ncomponents.\n\nNow let's see how we would modify our signals if we removed this component\nfrom the data\n\n" ], "cell_type": "markdown", "metadata": {} }, { "execution_count": null, "cell_type": "code", "source": [ "ica.plot_overlay(eog_average, exclude=eog_inds, show=False)\n# red -> before, black -> after. Yes! We remove quite a lot!\n\n# to definitely register this component as a bad one to be removed\n# there is the ``ica.exclude`` attribute, a simple Python list\nica.exclude.extend(eog_inds)\n\n# from now on the ICA will reject this component even if no exclude\n# parameter is passed, and this information will be stored to disk\n# on saving\n\n# uncomment this for reading and writing\n# ica.save('my-ica.fif')\n# ica = read_ica('my-ica.fif')" ], "outputs": [], "metadata": { "collapsed": false } }, { "source": [ "Exercise: find and remove ECG artifacts using ICA!\n\n" ], "cell_type": "markdown", "metadata": {} }, { "execution_count": null, "cell_type": "code", "source": [ "ecg_epochs = create_ecg_epochs(raw, tmin=-.5, tmax=.5)\necg_inds, scores = ica.find_bads_ecg(ecg_epochs, method='ctps')\nica.plot_properties(ecg_epochs, picks=ecg_inds, psd_args={'fmax': 35.})" ], "outputs": [], "metadata": { "collapsed": false } }, { "source": [ "What if we don't have an EOG channel?\n-------------------------------------\n\nWe could either:\n\n1. make a bipolar reference from frontal EEG sensors and use as virtual EOG\n channel. This can be tricky though as you can only hope that the frontal\n EEG channels only reflect EOG and not brain dynamics in the prefrontal\n cortex.\n2. go for a semi-automated approach, using template matching.\n\nIn MNE-Python option 2 is easily achievable and it might give better results,\nso let's have a look at it.\n\n" ], "cell_type": "markdown", "metadata": {} }, { "execution_count": null, "cell_type": "code", "source": [ "from mne.preprocessing.ica import corrmap # noqa" ], "outputs": [], "metadata": { "collapsed": false } }, { "source": [ "The idea behind corrmap is that artefact patterns are similar across subjects\nand can thus be identified by correlating the different patterns resulting\nfrom each solution with a template. The procedure is therefore\nsemi-automatic. :func:`mne.preprocessing.corrmap` hence takes a list of\nICA solutions and a template, that can be an index or an array.\n\nAs we don't have different subjects or runs available today, here we will\nsimulate ICA solutions from different subjects by fitting ICA models to\ndifferent parts of the same recording. Then we will use one of the components\nfrom our original ICA as a template in order to detect sufficiently similar\ncomponents in the simulated ICAs.\n\nThe following block of code simulates having ICA solutions from different\nruns/subjects so it should not be used in real analysis - use independent\ndata sets instead.\n\n" ], "cell_type": "markdown", "metadata": {} }, { "execution_count": null, "cell_type": "code", "source": [ "# We'll start by simulating a group of subjects or runs from a subject\nstart, stop = [0, len(raw.times) - 1]\nintervals = np.linspace(start, stop, 4, dtype=int)\nicas_from_other_data = list()\nraw.pick_types(meg=True, eeg=False) # take only MEG channels\nfor ii, start in enumerate(intervals):\n if ii + 1 < len(intervals):\n stop = intervals[ii + 1]\n print('fitting ICA from {0} to {1} seconds'.format(start, stop))\n this_ica = ICA(n_components=n_components, method=method).fit(\n raw, start=start, stop=stop, reject=reject)\n icas_from_other_data.append(this_ica)" ], "outputs": [], "metadata": { "collapsed": false } }, { "source": [ "Remember, don't do this at home! Start by reading in a collection of ICA\nsolutions instead. Something like:\n\n``icas = [mne.preprocessing.read_ica(fname) for fname in ica_fnames]``\n\n" ], "cell_type": "markdown", "metadata": {} }, { "execution_count": null, "cell_type": "code", "source": [ "print(icas_from_other_data)" ], "outputs": [], "metadata": { "collapsed": false } }, { "source": [ "We use our original ICA as reference.\n\n" ], "cell_type": "markdown", "metadata": {} }, { "execution_count": null, "cell_type": "code", "source": [ "reference_ica = ica" ], "outputs": [], "metadata": { "collapsed": false } }, { "source": [ "Investigate our reference ICA:\n\n" ], "cell_type": "markdown", "metadata": {} }, { "execution_count": null, "cell_type": "code", "source": [ "reference_ica.plot_components()" ], "outputs": [], "metadata": { "collapsed": false } }, { "source": [ "Which one is the bad EOG component?\nHere we rely on our previous detection algorithm. You would need to decide\nyourself if no automatic detection was available.\n\n" ], "cell_type": "markdown", "metadata": {} }, { "execution_count": null, "cell_type": "code", "source": [ "reference_ica.plot_sources(eog_average, exclude=eog_inds)" ], "outputs": [], "metadata": { "collapsed": false } }, { "source": [ "Indeed it looks like an EOG, also in the average time course.\n\nWe construct a list where our reference run is the first element. Then we\ncan detect similar components from the other runs (the other ICA objects)\nusing :func:`mne.preprocessing.corrmap`. So our template must be a tuple like\n(reference_run_index, component_index):\n\n" ], "cell_type": "markdown", "metadata": {} }, { "execution_count": null, "cell_type": "code", "source": [ "icas = [reference_ica] + icas_from_other_data\ntemplate = (0, eog_inds[0])" ], "outputs": [], "metadata": { "collapsed": false } }, { "source": [ "Now we can run the CORRMAP algorithm.\n\n" ], "cell_type": "markdown", "metadata": {} }, { "execution_count": null, "cell_type": "code", "source": [ "fig_template, fig_detected = corrmap(icas, template=template, label=\"blinks\",\n show=True, threshold=.8, ch_type='mag')" ], "outputs": [], "metadata": { "collapsed": false } }, { "source": [ "Nice, we have found similar ICs from the other (simulated) runs!\nIn this way, you can detect a type of artifact semi-automatically for example\nfor all subjects in a study.\nThe detected template can also be retrieved as an array and stored; this\narray can be used as an alternative template to\n:func:`mne.preprocessing.corrmap`.\n\n" ], "cell_type": "markdown", "metadata": {} }, { "execution_count": null, "cell_type": "code", "source": [ "eog_component = reference_ica.get_components()[:, eog_inds[0]]\n\n# If you calculate a new ICA solution, you can provide this array instead of\n# specifying the template in reference to the list of ICA objects you want\n# to run CORRMAP on. (Of course, the retrieved component map arrays can\n# also be used for other purposes than artifact correction.)\n#\n# You can also use SSP to correct for artifacts. It is a bit simpler and\n# faster but also less precise than ICA and requires that you know the event\n# timing of your artifact.\n# See :ref:`tut_artifacts_correct_ssp`." ], "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" } } } }