{
  "nbformat_minor": 0, 
  "nbformat": 4, 
  "cells": [
    {
      "execution_count": null, 
      "cell_type": "code", 
      "source": [
        "%matplotlib inline"
      ], 
      "outputs": [], 
      "metadata": {
        "collapsed": false
      }
    }, 
    {
      "source": [
        "\n# Compute the power spectral density of raw data\n\n\nThis script shows how to compute the power spectral density (PSD)\nof measurements on a raw dataset. It also show the effect of applying SSP\nto the data to reduce ECG and EOG artifacts.\n\n"
      ], 
      "cell_type": "markdown", 
      "metadata": {}
    }, 
    {
      "execution_count": null, 
      "cell_type": "code", 
      "source": [
        "# Authors: Alexandre Gramfort <alexandre.gramfort@telecom-paristech.fr>\n#          Martin Luessi <mluessi@nmr.mgh.harvard.edu>\n#          Eric Larson <larson.eric.d@gmail.com>\n# License: BSD (3-clause)\n\nimport numpy as np\nimport matplotlib.pyplot as plt\n\nimport mne\nfrom mne import io, read_proj, read_selection\nfrom mne.datasets import sample\nfrom mne.time_frequency import psd_multitaper\n\nprint(__doc__)"
      ], 
      "outputs": [], 
      "metadata": {
        "collapsed": false
      }
    }, 
    {
      "source": [
        "Set parameters\n\n"
      ], 
      "cell_type": "markdown", 
      "metadata": {}
    }, 
    {
      "execution_count": null, 
      "cell_type": "code", 
      "source": [
        "data_path = sample.data_path()\nraw_fname = data_path + '/MEG/sample/sample_audvis_raw.fif'\nproj_fname = data_path + '/MEG/sample/sample_audvis_eog-proj.fif'\n\ntmin, tmax = 0, 60  # use the first 60s of data\n\n# Setup for reading the raw data (to save memory, crop before loading)\nraw = io.read_raw_fif(raw_fname).crop(tmin, tmax).load_data()\nraw.info['bads'] += ['MEG 2443', 'EEG 053']  # bads + 2 more\n\n# Add SSP projection vectors to reduce EOG and ECG artifacts\nprojs = read_proj(proj_fname)\nraw.add_proj(projs, remove_existing=True)\n\n\nfmin, fmax = 2, 300  # look at frequencies between 2 and 300Hz\nn_fft = 2048  # the FFT size (n_fft). Ideally a power of 2\n\n# Let's first check out all channel types\nraw.plot_psd(area_mode='range', tmax=10.0, show=False)\n\n# Now let's focus on a smaller subset:\n# Pick MEG magnetometers in the Left-temporal region\nselection = read_selection('Left-temporal')\npicks = mne.pick_types(raw.info, meg='mag', eeg=False, eog=False,\n                       stim=False, exclude='bads', selection=selection)\n\n# Let's just look at the first few channels for demonstration purposes\npicks = picks[:4]\n\nplt.figure()\nax = plt.axes()\nraw.plot_psd(tmin=tmin, tmax=tmax, fmin=fmin, fmax=fmax, n_fft=n_fft,\n             n_jobs=1, proj=False, ax=ax, color=(0, 0, 1),  picks=picks,\n             show=False)\n\n# And now do the same with SSP applied\nraw.plot_psd(tmin=tmin, tmax=tmax, fmin=fmin, fmax=fmax, n_fft=n_fft,\n             n_jobs=1, proj=True, ax=ax, color=(0, 1, 0), picks=picks,\n             show=False)\n\n# And now do the same with SSP + notch filtering\n# Pick all channels for notch since the SSP projection mixes channels together\nraw.notch_filter(np.arange(60, 241, 60), n_jobs=1)\nraw.plot_psd(tmin=tmin, tmax=tmax, fmin=fmin, fmax=fmax, n_fft=n_fft,\n             n_jobs=1, proj=True, ax=ax, color=(1, 0, 0), picks=picks,\n             show=False)\n\nax.set_title('Four left-temporal magnetometers')\nplt.legend(['Without SSP', 'With SSP', 'SSP + Notch'])\n\n# Alternatively, you may also create PSDs from Raw objects with ``psd_*``\nf, ax = plt.subplots()\npsds, freqs = psd_multitaper(raw, low_bias=True, tmin=tmin, tmax=tmax,\n                             fmin=fmin, fmax=fmax, proj=True, picks=picks,\n                             n_jobs=1)\npsds = 10 * np.log10(psds)\npsds_mean = psds.mean(0)\npsds_std = psds.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', xlabel='Frequency',\n       ylabel='Power Spectral Density (dB)')\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"
      }
    }
  }
}