{
  "nbformat_minor": 0, 
  "nbformat": 4, 
  "cells": [
    {
      "execution_count": null, 
      "cell_type": "code", 
      "source": [
        "%matplotlib inline"
      ], 
      "outputs": [], 
      "metadata": {
        "collapsed": false
      }
    }, 
    {
      "source": [
        "\n\n# Frequency and time-frequency sensors analysis\n\n\nThe objective is to show you how to explore the spectral content\nof your data (frequency and time-frequency). Here we'll work on Epochs.\n\nWe will use the somatosensory dataset that contains so\ncalled event related synchronizations (ERS) / desynchronizations (ERD) in\nthe beta band.\n\n"
      ], 
      "cell_type": "markdown", 
      "metadata": {}
    }, 
    {
      "execution_count": null, 
      "cell_type": "code", 
      "source": [
        "import numpy as np\nimport matplotlib.pyplot as plt\n\nimport mne\nfrom mne.time_frequency import tfr_morlet, psd_multitaper\nfrom mne.datasets import somato"
      ], 
      "outputs": [], 
      "metadata": {
        "collapsed": false
      }
    }, 
    {
      "source": [
        "Set parameters\n\n"
      ], 
      "cell_type": "markdown", 
      "metadata": {}
    }, 
    {
      "execution_count": null, 
      "cell_type": "code", 
      "source": [
        "data_path = somato.data_path()\nraw_fname = data_path + '/MEG/somato/sef_raw_sss.fif'\n\n# Setup for reading the raw data\nraw = mne.io.read_raw_fif(raw_fname)\nevents = mne.find_events(raw, stim_channel='STI 014')\n\n# picks MEG gradiometers\npicks = mne.pick_types(raw.info, meg='grad', eeg=False, eog=True, stim=False)\n\n# Construct Epochs\nevent_id, tmin, tmax = 1, -1., 3.\nbaseline = (None, 0)\nepochs = mne.Epochs(raw, events, event_id, tmin, tmax, picks=picks,\n                    baseline=baseline, reject=dict(grad=4000e-13, eog=350e-6),\n                    preload=True)\n\nepochs.resample(150., npad='auto')  # resample to reduce computation time"
      ], 
      "outputs": [], 
      "metadata": {
        "collapsed": false
      }
    }, 
    {
      "source": [
        "Frequency analysis\n------------------\n\nWe 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": "markdown", 
      "metadata": {}
    }, 
    {
      "execution_count": null, 
      "cell_type": "code", 
      "source": [
        "epochs.plot_psd(fmin=2., fmax=40.)"
      ], 
      "outputs": [], 
      "metadata": {
        "collapsed": false
      }
    }, 
    {
      "source": [
        "Now let's take a look at the spatial distributions of the PSD.\n\n"
      ], 
      "cell_type": "markdown", 
      "metadata": {}
    }, 
    {
      "execution_count": null, 
      "cell_type": "code", 
      "source": [
        "epochs.plot_psd_topomap(ch_type='grad', normalize=True)"
      ], 
      "outputs": [], 
      "metadata": {
        "collapsed": false
      }
    }, 
    {
      "source": [
        "Alternatively, you can also create PSDs from Epochs objects with functions\nthat start with ``psd_`` such as\n:func:`mne.time_frequency.psd_multitaper` and\n:func:`mne.time_frequency.psd_welch`.\n\n"
      ], 
      "cell_type": "markdown", 
      "metadata": {}
    }, 
    {
      "execution_count": null, 
      "cell_type": "code", 
      "source": [
        "f, ax = plt.subplots()\npsds, freqs = psd_multitaper(epochs, fmin=2, fmax=40, n_jobs=1)\npsds = 10 * np.log10(psds)\npsds_mean = psds.mean(0).mean(0)\npsds_std = psds.mean(0).std(0)\n\nax.plot(freqs, psds_mean, color='k')\nax.fill_between(freqs, psds_mean - psds_std, psds_mean + psds_std,\n                color='k', alpha=.5)\nax.set(title='Multitaper PSD (gradiometers)', xlabel='Frequency',\n       ylabel='Power Spectral Density (dB)')\nplt.show()"
      ], 
      "outputs": [], 
      "metadata": {
        "collapsed": false
      }
    }, 
    {
      "source": [
        "Time-frequency analysis: power and intertrial coherence\n-------------------------------------------------------\n\nWe now compute time-frequency representations (TFRs) from our Epochs.\nWe'll look at power and intertrial coherence (ITC).\n\nTo this we'll use the function :func:`mne.time_frequency.tfr_morlet`\nbut you can also use :func:`mne.time_frequency.tfr_multitaper`\nor :func:`mne.time_frequency.tfr_stockwell`.\n\n"
      ], 
      "cell_type": "markdown", 
      "metadata": {}
    }, 
    {
      "execution_count": null, 
      "cell_type": "code", 
      "source": [
        "# define frequencies of interest (log-spaced)\nfreqs = np.logspace(*np.log10([6, 35]), num=8)\nn_cycles = freqs / 2.  # different number of cycle per frequency\npower, itc = tfr_morlet(epochs, freqs=freqs, n_cycles=n_cycles, use_fft=True,\n                        return_itc=True, decim=3, n_jobs=1)"
      ], 
      "outputs": [], 
      "metadata": {
        "collapsed": false
      }
    }, 
    {
      "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 censor.\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": "markdown", 
      "metadata": {}
    }, 
    {
      "execution_count": null, 
      "cell_type": "code", 
      "source": [
        "power.plot_topo(baseline=(-0.5, 0), mode='logratio', title='Average power')\npower.plot([82], baseline=(-0.5, 0), mode='logratio')\n\nfig, axis = plt.subplots(1, 2, figsize=(7, 4))\npower.plot_topomap(ch_type='grad', tmin=0.5, tmax=1.5, fmin=8, fmax=12,\n                   baseline=(-0.5, 0), mode='logratio', axes=axis[0],\n                   title='Alpha', vmax=0.45, show=False)\npower.plot_topomap(ch_type='grad', tmin=0.5, tmax=1.5, fmin=13, fmax=25,\n                   baseline=(-0.5, 0), mode='logratio', axes=axis[1],\n                   title='Beta', vmax=0.45, show=False)\nmne.viz.tight_layout()\nplt.show()"
      ], 
      "outputs": [], 
      "metadata": {
        "collapsed": false
      }
    }, 
    {
      "source": [
        "Inspect ITC\n-----------\n\n"
      ], 
      "cell_type": "markdown", 
      "metadata": {}
    }, 
    {
      "execution_count": null, 
      "cell_type": "code", 
      "source": [
        "itc.plot_topo(title='Inter-Trial coherence', vmin=0., vmax=1., cmap='Reds')"
      ], 
      "outputs": [], 
      "metadata": {
        "collapsed": false
      }
    }, 
    {
      "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"
      ], 
      "cell_type": "markdown", 
      "metadata": {}
    }, 
    {
      "source": [
        "Exercise\n--------\n\n   - Visualize the intertrial coherence values as topomaps as done with\n     power.\n\n"
      ], 
      "cell_type": "markdown", 
      "metadata": {}
    }
  ], 
  "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"
      }
    }
  }
}