{
  "nbformat_minor": 0, 
  "nbformat": 4, 
  "cells": [
    {
      "execution_count": null, 
      "cell_type": "code", 
      "source": [
        "%matplotlib inline"
      ], 
      "outputs": [], 
      "metadata": {
        "collapsed": false
      }
    }, 
    {
      "source": [
        "\n========================================================\nTime-frequency on simulated data (Multitaper vs. Morlet)\n========================================================\n\nThis examples demonstrates on simulated data the different time-frequency\nestimation methods. It shows the time-frequency resolution trade-off\nand the problem of estimation variance.\n\n"
      ], 
      "cell_type": "markdown", 
      "metadata": {}
    }, 
    {
      "execution_count": null, 
      "cell_type": "code", 
      "source": [
        "# Authors: Hari Bharadwaj <hari@nmr.mgh.harvard.edu>\n#          Denis Engemann <denis.engemann@gmail.com>\n#\n# License: BSD (3-clause)\n\nimport numpy as np\n\nfrom mne import create_info, EpochsArray\nfrom mne.time_frequency import tfr_multitaper, tfr_stockwell, tfr_morlet\n\nprint(__doc__)"
      ], 
      "outputs": [], 
      "metadata": {
        "collapsed": false
      }
    }, 
    {
      "source": [
        "Simulate data\n\n"
      ], 
      "cell_type": "markdown", 
      "metadata": {}
    }, 
    {
      "execution_count": null, 
      "cell_type": "code", 
      "source": [
        "sfreq = 1000.0\nch_names = ['SIM0001', 'SIM0002']\nch_types = ['grad', 'grad']\ninfo = create_info(ch_names=ch_names, sfreq=sfreq, ch_types=ch_types)\n\nn_times = int(sfreq)  # 1 second long epochs\nn_epochs = 40\nseed = 42\nrng = np.random.RandomState(seed)\nnoise = rng.randn(n_epochs, len(ch_names), n_times)\n\n# Add a 50 Hz sinusoidal burst to the noise and ramp it.\nt = np.arange(n_times, dtype=np.float) / sfreq\nsignal = np.sin(np.pi * 2. * 50. * t)  # 50 Hz sinusoid signal\nsignal[np.logical_or(t < 0.45, t > 0.55)] = 0.  # Hard windowing\non_time = np.logical_and(t >= 0.45, t <= 0.55)\nsignal[on_time] *= np.hanning(on_time.sum())  # Ramping\ndata = noise + signal\n\nreject = dict(grad=4000)\nevents = np.empty((n_epochs, 3), dtype=int)\nfirst_event_sample = 100\nevent_id = dict(sin50hz=1)\nfor k in range(n_epochs):\n    events[k, :] = first_event_sample + k * n_times, 0, event_id['sin50hz']\n\nepochs = EpochsArray(data=data, info=info, events=events, event_id=event_id,\n                     reject=reject)"
      ], 
      "outputs": [], 
      "metadata": {
        "collapsed": false
      }
    }, 
    {
      "source": [
        "Consider different parameter possibilities for multitaper convolution\n\n"
      ], 
      "cell_type": "markdown", 
      "metadata": {}
    }, 
    {
      "execution_count": null, 
      "cell_type": "code", 
      "source": [
        "freqs = np.arange(5., 100., 3.)\n\n# You can trade time resolution or frequency resolution or both\n# in order to get a reduction in variance\n\n# (1) Least smoothing (most variance/background fluctuations).\nn_cycles = freqs / 2.\ntime_bandwidth = 2.0  # Least possible frequency-smoothing (1 taper)\npower = tfr_multitaper(epochs, freqs=freqs, n_cycles=n_cycles,\n                       time_bandwidth=time_bandwidth, return_itc=False)\n# Plot results. Baseline correct based on first 100 ms.\npower.plot([0], baseline=(0., 0.1), mode='mean', vmin=-1., vmax=3.,\n           title='Sim: Least smoothing, most variance')\n\n\n# (2) Less frequency smoothing, more time smoothing.\nn_cycles = freqs  # Increase time-window length to 1 second.\ntime_bandwidth = 4.0  # Same frequency-smoothing as (1) 3 tapers.\npower = tfr_multitaper(epochs, freqs=freqs, n_cycles=n_cycles,\n                       time_bandwidth=time_bandwidth, return_itc=False)\n# Plot results. Baseline correct based on first 100 ms.\npower.plot([0], baseline=(0., 0.1), mode='mean', vmin=-1., vmax=3.,\n           title='Sim: Less frequency smoothing, more time smoothing')\n\n\n# (3) Less time smoothing, more frequency smoothing.\nn_cycles = freqs / 2.\ntime_bandwidth = 8.0  # Same time-smoothing as (1), 7 tapers.\npower = tfr_multitaper(epochs, freqs=freqs, n_cycles=n_cycles,\n                       time_bandwidth=time_bandwidth, return_itc=False)\n# Plot results. Baseline correct based on first 100 ms.\npower.plot([0], baseline=(0., 0.1), mode='mean', vmin=-1., vmax=3.,\n           title='Sim: Less time smoothing, more frequency smoothing')\n\n# #############################################################################\n# Stockwell (S) transform\n\n# S uses a Gaussian window to balance temporal and spectral resolution\n# Importantly, frequency bands are phase-normalized, hence strictly comparable\n# with regard to timing, and, the input signal can be recoverd from the\n# transform in a lossless way if we disregard numerical errors.\n\nfmin, fmax = freqs[[0, -1]]\nfor width in (0.7, 3.0):\n    power = tfr_stockwell(epochs, fmin=fmin, fmax=fmax, width=width)\n    power.plot([0], baseline=(0., 0.1), mode='mean',\n               title='Sim: Using S transform, width '\n                     '= {:0.1f}'.format(width), show=True)\n\n# #############################################################################\n# Finally, compare to morlet wavelet\n\nn_cycles = freqs / 2.\npower = tfr_morlet(epochs, freqs=freqs, n_cycles=n_cycles, return_itc=False)\npower.plot([0], baseline=(0., 0.1), mode='mean', vmin=-1., vmax=3.,\n           title='Sim: Using Morlet wavelet')"
      ], 
      "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"
      }
    }
  }
}