{
  "nbformat_minor": 0, 
  "nbformat": 4, 
  "cells": [
    {
      "execution_count": null, 
      "cell_type": "code", 
      "source": [
        "%matplotlib inline"
      ], 
      "outputs": [], 
      "metadata": {
        "collapsed": false
      }
    }, 
    {
      "source": [
        "\n# Compute power and phase lock in label of the source space\n\n\nCompute time-frequency maps of power and phase lock in the source space.\nThe inverse method is linear based on dSPM inverse operator.\n\nThe example also shows the difference in the time-frequency maps\nwhen they are computed with and without subtracting the evoked response\nfrom each epoch. The former results in induced activity only while the\nlatter also includes evoked (stimulus-locked) activity.\n\n"
      ], 
      "cell_type": "markdown", 
      "metadata": {}
    }, 
    {
      "execution_count": null, 
      "cell_type": "code", 
      "source": [
        "# Authors: Alexandre Gramfort <alexandre.gramfort@telecom-paristech.fr>\n#\n# License: BSD (3-clause)\n\nimport numpy as np\nimport matplotlib.pyplot as plt\n\nimport mne\nfrom mne import io\nfrom mne.datasets import sample\nfrom mne.minimum_norm import read_inverse_operator, source_induced_power\n\nprint(__doc__)"
      ], 
      "outputs": [], 
      "metadata": {
        "collapsed": false
      }
    }, 
    {
      "source": [
        "Set parameters\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_raw.fif'\nfname_inv = data_path + '/MEG/sample/sample_audvis-meg-oct-6-meg-inv.fif'\nlabel_name = 'Aud-rh'\nfname_label = data_path + '/MEG/sample/labels/%s.label' % label_name\n\ntmin, tmax, event_id = -0.2, 0.5, 2\n\n# Setup for reading the raw data\nraw = io.read_raw_fif(raw_fname)\nevents = mne.find_events(raw, stim_channel='STI 014')\ninverse_operator = read_inverse_operator(fname_inv)\n\ninclude = []\nraw.info['bads'] += ['MEG 2443', 'EEG 053']  # bads + 2 more\n\n# Picks MEG channels\npicks = mne.pick_types(raw.info, meg=True, eeg=False, eog=True,\n                       stim=False, include=include, exclude='bads')\nreject = dict(grad=4000e-13, mag=4e-12, eog=150e-6)\n\n# Load epochs\nepochs = mne.Epochs(raw, events, event_id, tmin, tmax, picks=picks,\n                    baseline=(None, 0), reject=reject,\n                    preload=True)\n\n# Compute a source estimate per frequency band including and excluding the\n# evoked response\nfrequencies = np.arange(7, 30, 2)  # define frequencies of interest\nlabel = mne.read_label(fname_label)\nn_cycles = frequencies / 3.  # different number of cycle per frequency\n\n# subtract the evoked response in order to exclude evoked activity\nepochs_induced = epochs.copy().subtract_evoked()\n\nplt.close('all')\n\nfor ii, (this_epochs, title) in enumerate(zip([epochs, epochs_induced],\n                                              ['evoked + induced',\n                                               'induced only'])):\n    # compute the source space power and the inter-trial coherence\n    power, itc = source_induced_power(\n        this_epochs, inverse_operator, frequencies, label, baseline=(-0.1, 0),\n        baseline_mode='percent', n_cycles=n_cycles, n_jobs=1)\n\n    power = np.mean(power, axis=0)  # average over sources\n    itc = np.mean(itc, axis=0)  # average over sources\n    times = epochs.times\n\n    ##########################################################################\n    # View time-frequency plots\n    plt.subplots_adjust(0.1, 0.08, 0.96, 0.94, 0.2, 0.43)\n    plt.subplot(2, 2, 2 * ii + 1)\n    plt.imshow(20 * power,\n               extent=[times[0], times[-1], frequencies[0], frequencies[-1]],\n               aspect='auto', origin='lower', vmin=0., vmax=30., cmap='RdBu_r')\n    plt.xlabel('Time (s)')\n    plt.ylabel('Frequency (Hz)')\n    plt.title('Power (%s)' % title)\n    plt.colorbar()\n\n    plt.subplot(2, 2, 2 * ii + 2)\n    plt.imshow(itc,\n               extent=[times[0], times[-1], frequencies[0], frequencies[-1]],\n               aspect='auto', origin='lower', vmin=0, vmax=0.7,\n               cmap='RdBu_r')\n    plt.xlabel('Time (s)')\n    plt.ylabel('Frequency (Hz)')\n    plt.title('ITC (%s)' % title)\n    plt.colorbar()\n\nplt.show()"
      ], 
      "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"
      }
    }
  }
}