{
  "nbformat_minor": 0, 
  "nbformat": 4, 
  "cells": [
    {
      "execution_count": null, 
      "cell_type": "code", 
      "source": [
        "%matplotlib inline"
      ], 
      "outputs": [], 
      "metadata": {
        "collapsed": false
      }
    }, 
    {
      "source": [
        "\n# Compute Power Spectral Density of inverse solution from single epochs\n\n\nCompute PSD of dSPM inverse solution on single trial epochs restricted\nto a brain label. The PSD is computed using a multi-taper method with\nDiscrete Prolate Spheroidal Sequence (DPSS) windows.\n\n\n"
      ], 
      "cell_type": "markdown", 
      "metadata": {}
    }, 
    {
      "execution_count": null, 
      "cell_type": "code", 
      "source": [
        "# Author: Martin Luessi <mluessi@nmr.mgh.harvard.edu>\n#\n# License: BSD (3-clause)\n\nimport numpy as np\nimport matplotlib.pyplot as plt\n\nimport mne\nfrom mne.datasets import sample\nfrom mne.minimum_norm import read_inverse_operator, compute_source_psd_epochs\n\nprint(__doc__)\n\ndata_path = sample.data_path()\nfname_inv = data_path + '/MEG/sample/sample_audvis-meg-oct-6-meg-inv.fif'\nfname_raw = data_path + '/MEG/sample/sample_audvis_raw.fif'\nfname_event = data_path + '/MEG/sample/sample_audvis_raw-eve.fif'\nlabel_name = 'Aud-lh'\nfname_label = data_path + '/MEG/sample/labels/%s.label' % label_name\n\nevent_id, tmin, tmax = 1, -0.2, 0.5\nsnr = 1.0  # use smaller SNR for raw data\nlambda2 = 1.0 / snr ** 2\nmethod = \"dSPM\"  # use dSPM method (could also be MNE or sLORETA)\n\n# Load data\ninverse_operator = read_inverse_operator(fname_inv)\nlabel = mne.read_label(fname_label)\nraw = mne.io.read_raw_fif(fname_raw)\nevents = mne.read_events(fname_event)\n\n# Set up pick list\ninclude = []\nraw.info['bads'] += ['EEG 053']  # bads + 1 more\n\n# pick MEG channels\npicks = mne.pick_types(raw.info, meg=True, eeg=False, stim=False, eog=True,\n                       include=include, exclude='bads')\n# Read epochs\nepochs = mne.Epochs(raw, events, event_id, tmin, tmax, picks=picks,\n                    baseline=(None, 0), reject=dict(mag=4e-12, grad=4000e-13,\n                                                    eog=150e-6))\n\n# define frequencies of interest\nfmin, fmax = 0., 70.\nbandwidth = 4.  # bandwidth of the windows in Hz\n\n# compute source space psd in label\n\n# Note: By using \"return_generator=True\" stcs will be a generator object\n# instead of a list. This allows us so to iterate without having to\n# keep everything in memory.\n\nstcs = compute_source_psd_epochs(epochs, inverse_operator, lambda2=lambda2,\n                                 method=method, fmin=fmin, fmax=fmax,\n                                 bandwidth=bandwidth, label=label,\n                                 return_generator=True)\n\n# compute average PSD over the first 10 epochs\nn_epochs = 10\nfor i, stc in enumerate(stcs):\n    if i >= n_epochs:\n        break\n\n    if i == 0:\n        psd_avg = np.mean(stc.data, axis=0)\n    else:\n        psd_avg += np.mean(stc.data, axis=0)\n\npsd_avg /= n_epochs\nfreqs = stc.times  # the frequencies are stored here\n\nplt.figure()\nplt.plot(freqs, psd_avg)\nplt.xlabel('Freq (Hz)')\nplt.ylabel('Power Spectral Density')\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"
      }
    }
  }
}