{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n\n# Compute power and phase lock in label of the source space\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"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Authors: Alexandre Gramfort <alexandre.gramfort@inria.fr>\n#\n# License: BSD-3-Clause\n# Copyright the MNE-Python contributors."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import matplotlib.pyplot as plt\nimport numpy as np\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__)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Set parameters\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "data_path = sample.data_path()\nmeg_path = data_path / \"MEG\" / \"sample\"\nraw_fname = meg_path / \"sample_audvis_raw.fif\"\nfname_inv = meg_path / \"sample_audvis-meg-oct-6-meg-inv.fif\"\nlabel_names = [\"Aud-lh\", \"Aud-rh\"]\nfname_labels = [meg_path / \"labels\" / f\"{ln}.label\" for ln in label_names]\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(\n    raw.info, meg=True, eeg=False, eog=True, stim=False, include=include, exclude=\"bads\"\n)\nreject = dict(grad=4000e-13, mag=4e-12, eog=150e-6)\n\n# Load epochs\nepochs = mne.Epochs(\n    raw,\n    events,\n    event_id,\n    tmin,\n    tmax,\n    picks=picks,\n    baseline=(None, 0),\n    reject=reject,\n    preload=True,\n)\n\n# Compute a source estimate per frequency band including and excluding the\n# evoked response\nfreqs = np.arange(7, 30, 2)  # define frequencies of interest\nlabels = [mne.read_label(fl) for fl in fname_labels]\nlabel = labels[0]\nn_cycles = freqs / 3.0  # 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\nfig, axes = plt.subplots(2, 2, layout=\"constrained\")\nfor ii, (this_epochs, title) in enumerate(\n    zip([epochs, epochs_induced], [\"evoked + induced\", \"induced only\"])\n):\n    # compute the source space power and the inter-trial coherence\n    power, itc = source_induced_power(\n        this_epochs,\n        inverse_operator,\n        freqs,\n        label,\n        baseline=(-0.1, 0),\n        baseline_mode=\"percent\",\n        n_cycles=n_cycles,\n        n_jobs=None,\n    )\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    ax = axes[ii, 0]\n    ax.imshow(\n        20 * power,\n        extent=[times[0], times[-1], freqs[0], freqs[-1]],\n        aspect=\"auto\",\n        origin=\"lower\",\n        vmin=0.0,\n        vmax=30.0,\n        cmap=\"RdBu_r\",\n    )\n    ax.set(xlabel=\"Time (s)\", ylabel=\"Frequency (Hz)\", title=f\"Power ({title})\")\n\n    ax = axes[ii, 1]\n    ax.imshow(\n        itc,\n        extent=[times[0], times[-1], freqs[0], freqs[-1]],\n        aspect=\"auto\",\n        origin=\"lower\",\n        vmin=0,\n        vmax=0.7,\n        cmap=\"RdBu_r\",\n    )\n    ax.set(xlabel=\"Time (s)\", ylabel=\"Frequency (Hz)\", title=f\"ITC ({title})\")\n    fig.colorbar(ax.images[0], ax=axes[ii])"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "In the example above, we averaged power across vertices after calculating\npower because we provided a single label for power calculation and therefore\npower of all sources within the single label were returned separately. When\nwe provide a list of labels, power is averaged across sources within each\nlabel automatically. With a list of labels, averaging is performed before\nrescaling, so choose a baseline method appropriately.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Get power from multiple labels\nmulti_label_power = source_induced_power(\n    epochs,\n    inverse_operator,\n    freqs,\n    labels,\n    baseline=(-0.1, 0),\n    baseline_mode=\"mean\",\n    n_cycles=n_cycles,\n    n_jobs=None,\n    return_plv=False,\n)\n\n# visually compare evoked power in left and right auditory regions\nfig, axes = plt.subplots(ncols=2, layout=\"constrained\")\nfor l_idx, l_power in enumerate(multi_label_power):\n    ax = axes[l_idx]\n    ax.imshow(\n        l_power,\n        extent=[epochs.times[0], epochs.times[-1], freqs[0], freqs[-1]],\n        aspect=\"auto\",\n        origin=\"lower\",\n        vmin=multi_label_power.min(),\n        vmax=multi_label_power.max(),\n        cmap=\"RdBu_r\",\n    )\n    title = f\"{labels[l_idx].hemi.upper()} Evoked Power\"\n    ax.set(xlabel=\"Time (s)\", ylabel=\"Frequency (Hz)\", title=title)\n    fig.colorbar(ax.images[0], ax=ax)"
      ]
    }
  ],
  "metadata": {
    "kernelspec": {
      "display_name": "Python 3",
      "language": "python",
      "name": "python3"
    },
    "language_info": {
      "codemirror_mode": {
        "name": "ipython",
        "version": 3
      },
      "file_extension": ".py",
      "mimetype": "text/x-python",
      "name": "python",
      "nbconvert_exporter": "python",
      "pygments_lexer": "ipython3",
      "version": "3.12.2"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}