{
  "nbformat_minor": 0, 
  "nbformat": 4, 
  "cells": [
    {
      "execution_count": null, 
      "cell_type": "code", 
      "source": [
        "%matplotlib inline"
      ], 
      "outputs": [], 
      "metadata": {
        "collapsed": false
      }
    }, 
    {
      "source": [
        "\n# Temporal whitening with AR model\n\n\nHere we fit an AR model to the data and use it\nto temporally whiten the signals.\n\n\n"
      ], 
      "cell_type": "markdown", 
      "metadata": {}
    }, 
    {
      "execution_count": null, 
      "cell_type": "code", 
      "source": [
        "# Authors: Alexandre Gramfort <alexandre.gramfort@telecom-paristech.fr>\n#\n# License: BSD (3-clause)\n\nimport numpy as np\nfrom scipy import signal\nimport matplotlib.pyplot as plt\n\nimport mne\nfrom mne.time_frequency import fit_iir_model_raw\nfrom mne.datasets import sample\n\nprint(__doc__)\n\ndata_path = sample.data_path()\n\nraw_fname = data_path + '/MEG/sample/sample_audvis_raw.fif'\nproj_fname = data_path + '/MEG/sample/sample_audvis_ecg-proj.fif'\n\nraw = mne.io.read_raw_fif(raw_fname)\nproj = mne.read_proj(proj_fname)\nraw.info['projs'] += proj\nraw.info['bads'] = ['MEG 2443', 'EEG 053']  # mark bad channels\n\n# Set up pick list: Gradiometers - bad channels\npicks = mne.pick_types(raw.info, meg='grad', exclude='bads')\n\norder = 5  # define model order\npicks = picks[:1]\n\n# Estimate AR models on raw data\nb, a = fit_iir_model_raw(raw, order=order, picks=picks, tmin=60, tmax=180)\nd, times = raw[0, 10000:20000]  # look at one channel from now on\nd = d.ravel()  # make flat vector\ninnovation = signal.convolve(d, a, 'valid')\nd_ = signal.lfilter(b, a, innovation)  # regenerate the signal\nd_ = np.r_[d_[0] * np.ones(order), d_]  # dummy samples to keep signal length"
      ], 
      "outputs": [], 
      "metadata": {
        "collapsed": false
      }
    }, 
    {
      "source": [
        "Plot the different time series and PSDs\n\n"
      ], 
      "cell_type": "markdown", 
      "metadata": {}
    }, 
    {
      "execution_count": null, 
      "cell_type": "code", 
      "source": [
        "plt.close('all')\nplt.figure()\nplt.plot(d[:100], label='signal')\nplt.plot(d_[:100], label='regenerated signal')\nplt.legend()\n\nplt.figure()\nplt.psd(d, Fs=raw.info['sfreq'], NFFT=2048)\nplt.psd(innovation, Fs=raw.info['sfreq'], NFFT=2048)\nplt.psd(d_, Fs=raw.info['sfreq'], NFFT=2048, linestyle='--')\nplt.legend(('Signal', 'Innovation', 'Regenerated signal'))\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"
      }
    }
  }
}