2-sensors_time_frequency.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Frequency and time-frequency sensors analysis\n",
    "\n",
    "\n",
    "The objective is to show you how to explore the spectral content\n",
    "of your data (frequency and time-frequency). Here we'll work on Epochs.\n",
    "\n",
    "    Authors: Alexandre Gramfort <alexandre.gramfort@inria.fr>\n",
    "             Stefan Appelhoff <stefan.appelhoff@mailbox.org>\n",
    "             Richard Höchenberger <richard.hoechenberger@gmail.com>\n",
    "             Denis A. Engemann <denis.engemann@gmail.com>\n",
    "    License: BSD (3-clause)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib qt\n",
    "import os\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "import mne\n",
    "from mne.time_frequency import tfr_morlet, psd_multitaper, psd_welch"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Set parameters\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "mne.set_log_level('WARNING')\n",
    "\n",
    "# Change the following path to where the folder ds000117-practical is on your disk\n",
    "data_path = os.path.expanduser(\"~/work/data/ds000117-practical/\")\n",
    "\n",
    "raw_fname = os.path.join(data_path,\n",
    "    'derivatives/meg_derivatives/sub-01/ses-meg/meg/sub-01_ses-meg_task-facerecognition_run-01_proc-sss_meg.fif')\n",
    "\n",
    "epochs_fname = raw_fname.replace('_meg.fif', '-epo.fif')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<EpochsFIF  |   141 events (all good), -0.5 - 2 sec, baseline [-0.2, 0], ~313.0 MB, data loaded,\n",
       " 'face/famous/first': 24\n",
       " 'face/famous/immediate': 10\n",
       " 'face/famous/long': 14\n",
       " 'face/unfamiliar/first': 24\n",
       " 'face/unfamiliar/immediate': 10\n",
       " 'face/unfamiliar/long': 10\n",
       " 'scrambled/first': 25\n",
       " 'scrambled/immediate': 13\n",
       " 'scrambled/long': 11>"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "epochs = mne.read_epochs(epochs_fname)\n",
    "epochs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# epochs.resample(200., npad='auto')  # resample to reduce computation time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<EpochsFIF  |   141 events (all good), -0.5 - 2 sec, baseline [-0.2, 0], ~313.0 MB, data loaded,\n",
       " 'face/famous/first': 24\n",
       " 'face/famous/immediate': 10\n",
       " 'face/famous/long': 14\n",
       " 'face/unfamiliar/first': 24\n",
       " 'face/unfamiliar/immediate': 10\n",
       " 'face/unfamiliar/long': 10\n",
       " 'scrambled/first': 25\n",
       " 'scrambled/immediate': 13\n",
       " 'scrambled/long': 11>"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "epochs.apply_proj()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Frequency analysis\n",
    "------------------\n",
    "\n",
    "We start by exploring the frequence content of our epochs.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's first check out all channel types by averaging across epochs.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "epochs.plot_psd(fmin=2., fmax=40., average=False, bandwidth=2);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-info\">\n",
    "    <b>REMARK</b>:\n",
    "     <ul>\n",
    "    <li> Select a frequency range in the plot to inspect topographies</li>\n",
    "\n",
    "<li> The \"bandwidth\" parameter controls the spectral resolution of the multitaper. You can increase the resolution by chosing a narrower bandwidth at the cost of longer computation time.</li>\n",
    "    </ul>\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's take a look at the spatial distributions of the PSD.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "epochs.plot_psd_topomap(ch_type='eeg', normalize=False, cmap='viridis');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "epochs.plot_psd_topomap(ch_type='mag', normalize=False, cmap='viridis');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "epochs.plot_psd_topomap(ch_type='grad', normalize=False, cmap='viridis');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-info\">\n",
    "    <b>REMARK</b>:\n",
    "     <ul>\n",
    "    <li>Sometimes it can be interesting  to consider the relative power, defined as the power in a given band divided by the total power. To explore this option, have a look at the \"normalize\" keyword. </li>\n",
    "    </ul>\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Alternatively, you can also create PSDs from Epochs objects with functions\n",
    "that start with ``psd_`` such as `mne.time_frequency.psd_multitaper` and\n",
    "`mne.time_frequency.psd_welch`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "psd_multitaper?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(141, 376, 95)"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "psds, freqs = psd_multitaper(epochs, fmin=2, fmax=40, n_jobs=1, bandwidth=2)\n",
    "psds.shape   # n_epochs x n_channels x n_frequencies"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "psds = np.mean(10. * np.log10(psds), axis=0)  # use dB and average over epochs\n",
    "picks_grad = mne.pick_types(epochs.info, meg=\"mag\", eeg=False)\n",
    "\n",
    "f, ax = plt.subplots()\n",
    "ax.plot(freqs, psds[picks_grad].T, color='k', alpha=0.3)\n",
    "ax.set(title='Multitaper PSD (gradiometers)', xlabel='Frequency (Hz)',\n",
    "       ylabel='Power Spectral Density (dB)')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can clearly see 2 peaks."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notably, `mne.time_frequency.psd_welch` supports the keyword argument\n",
    "``average``, which specifies how to estimate the PSD based on the individual\n",
    "windowed segments. The default is ``average='mean'``, which simply calculates\n",
    "the arithmetic mean across segments. Specifying ``average='median'``, in\n",
    "contrast, returns the PSD based on the median of the segments (corrected for\n",
    "bias relative to the mean), which is a more robust measure."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Lastly, we can also retrieve the unaggregated segments by passing\n",
    "``average=None`` to `mne.time_frequency.psd_welch`. The dimensions of\n",
    "the returned array are ``(n_epochs, n_sensors, n_freqs, n_segments)``.\n",
    "This can be interesting when one is interested in computing statistics across segments or use custom functions for aggregation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Time-frequency analysis: power and inter-trial coherence\n",
    "\n",
    "We now compute time-frequency representations (TFRs) from our Epochs.\n",
    "We'll look at power and inter-trial coherence (ITC).\n",
    "\n",
    "To this we'll use the function `mne.time_frequency.tfr_morlet`\n",
    "but you can also use `mne.time_frequency.tfr_multitaper`\n",
    "or `mne.time_frequency.tfr_stockwell`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "# define frequencies of interest (log-spaced)\n",
    "freqs = np.logspace(*np.log10([2, 30]), num=20)\n",
    "n_cycles = freqs / 2.  # different number of cycle per frequency\n",
    "power, itc = tfr_morlet(epochs, freqs=freqs, n_cycles=n_cycles, use_fft=True,\n",
    "                        return_itc=True, decim=3, n_jobs=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<AverageTFR  |  time : [-0.100000, 1.600000], freq : [2.000000, 30.000000], nave : 141, channels : 376, ~17.4 MB>"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "power.crop(-0.1, 1.6)  # crop to remove edge artifacts"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<AverageTFR  |  time : [-0.100000, 1.600000], freq : [2.000000, 30.000000], nave : 141, channels : 376, ~17.4 MB>"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "itc.crop(-0.1, 1.6)  # crop to remove edge artifacts"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Inspect power\n",
    "-------------\n",
    "\n",
    "<div class=\"alert alert-info\"><h4>Note</h4><p>The generated figures are interactive. In the topo you can click\n",
    "    on an image to visualize the data for one sensor.\n",
    "    You can also select a portion in the time-frequency plane to\n",
    "    obtain a topomap for a certain time-frequency region.</p></div>\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/dengeman/anaconda3/lib/python3.7/site-packages/matplotlib/cbook/__init__.py:216: RuntimeWarning: `mask` is None, not masking the plot ...\n",
      "  func(*args, **kwargs)\n",
      "/Users/dengeman/anaconda3/lib/python3.7/site-packages/matplotlib/cbook/__init__.py:216: RuntimeWarning: `mask` is None, not adding contour to the plot ...\n",
      "  func(*args, **kwargs)\n",
      "/Users/dengeman/anaconda3/lib/python3.7/site-packages/matplotlib/cbook/__init__.py:216: RuntimeWarning: `mask` is None, not masking the plot ...\n",
      "  func(*args, **kwargs)\n",
      "/Users/dengeman/anaconda3/lib/python3.7/site-packages/matplotlib/cbook/__init__.py:216: RuntimeWarning: `mask` is None, not adding contour to the plot ...\n",
      "  func(*args, **kwargs)\n"
     ]
    }
   ],
   "source": [
    "baseline_mode = 'logratio'\n",
    "baseline = (None, 0)\n",
    "power.plot_topo(baseline=baseline, mode=baseline_mode, title='Average power');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "power.plot([82], baseline=baseline, mode=baseline_mode, title=power.ch_names[82]);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axis = plt.subplots(1, 3, figsize=(7, 4))\n",
    "power.plot_topomap(ch_type='grad', tmin=0.5, tmax=1.5, fmin=4, fmax=7,\n",
    "                   baseline=baseline, mode=baseline_mode, axes=axis[0],\n",
    "                   title='Theta', show=False, contours=1)\n",
    "power.plot_topomap(ch_type='grad', tmin=0.5, tmax=1.5, fmin=8, fmax=12,\n",
    "                   baseline=baseline, mode=baseline_mode, axes=axis[1],\n",
    "                   title='Alpha', show=False, contours=1)\n",
    "power.plot_topomap(ch_type='grad', tmin=0.5, tmax=1.5, fmin=15, fmax=30,\n",
    "                   baseline=baseline, mode=baseline_mode, axes=axis[2],\n",
    "                   title='Beta', show=False, contours=1)\n",
    "mne.viz.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Joint Plot\n",
    "----------\n",
    "You can also create a joint plot showing both the aggregated TFR\n",
    "across channels and topomaps at specific times and frequencies to obtain\n",
    "a quick overview regarding oscillatory effects across time and space.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "power.plot_topo?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<Figure size 1280x960 with 4 Axes>,\n",
       " <Figure size 1280x960 with 4 Axes>,\n",
       " <Figure size 1280x960 with 4 Axes>]"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "power.plot_joint(baseline=baseline, mode='mean', tmin=None, tmax=None,\n",
    "                 timefreqs=[(0.25, 2.), (1., 11.)])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Inspect ITC\n",
    "-----------\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/dengeman/anaconda3/lib/python3.7/site-packages/matplotlib/cbook/__init__.py:216: RuntimeWarning: `mask` is None, not masking the plot ...\n",
      "  func(*args, **kwargs)\n",
      "/Users/dengeman/anaconda3/lib/python3.7/site-packages/matplotlib/cbook/__init__.py:216: RuntimeWarning: `mask` is None, not adding contour to the plot ...\n",
      "  func(*args, **kwargs)\n"
     ]
    }
   ],
   "source": [
    "itc.plot_topo(title='Inter-Trial coherence', vmin=0., vmax=1., cmap='Reds');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-info\"><h4>Note</h4><p>Baseline correction can be applied to power or done in plots.\n",
    "    To illustrate the baseline correction in plots, the next line is\n",
    "    commented power.apply_baseline(baseline=(-0.5, 0), mode='logratio')</p></div>\n",
    "\n",
    "<div class=\"alert alert-success\">\n",
    "    <b>EXERCISE</b>:\n",
    "     <ul>\n",
    "      <li>Visualize the inter-trial coherence values as topomaps as done with\n",
    "     power</li>\n",
    "    </ul>\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}