{
  "nbformat_minor": 0, 
  "nbformat": 4, 
  "cells": [
    {
      "execution_count": null, 
      "cell_type": "code", 
      "source": [
        "%matplotlib inline"
      ], 
      "outputs": [], 
      "metadata": {
        "collapsed": false
      }
    }, 
    {
      "source": [
        "\n# Brainstorm CTF phantom tutorial dataset\n\n\nHere we compute the evoked from raw for the Brainstorm CTF phantom\ntutorial dataset. For comparison, see [1]_ and:\n\n    http://neuroimage.usc.edu/brainstorm/Tutorials/PhantomCtf\n\nReferences\n----------\n.. [1] Tadel F, Baillet S, Mosher JC, Pantazis D, Leahy RM.\n       Brainstorm: A User-Friendly Application for MEG/EEG Analysis.\n       Computational Intelligence and Neuroscience, vol. 2011, Article ID\n       879716, 13 pages, 2011. doi:10.1155/2011/879716\n\n"
      ], 
      "cell_type": "markdown", 
      "metadata": {}
    }, 
    {
      "execution_count": null, 
      "cell_type": "code", 
      "source": [
        "# Authors: Eric Larson <larson.eric.d@gmail.com>\n#\n# License: BSD (3-clause)\n\nimport os.path as op\nimport numpy as np\nimport matplotlib.pyplot as plt\n\nimport mne\nfrom mne import fit_dipole\nfrom mne.datasets.brainstorm import bst_phantom_ctf\nfrom mne.io import read_raw_ctf\n\nprint(__doc__)"
      ], 
      "outputs": [], 
      "metadata": {
        "collapsed": false
      }
    }, 
    {
      "source": [
        "The data were collected with a CTF system at 2400 Hz.\n\n"
      ], 
      "cell_type": "markdown", 
      "metadata": {}
    }, 
    {
      "execution_count": null, 
      "cell_type": "code", 
      "source": [
        "data_path = bst_phantom_ctf.data_path()\n\n# Switch to these to use the higher-SNR data:\n# raw_path = op.join(data_path, 'phantom_200uA_20150709_01.ds')\n# dip_freq = 7.\nraw_path = op.join(data_path, 'phantom_20uA_20150603_03.ds')\ndip_freq = 23.\nerm_path = op.join(data_path, 'emptyroom_20150709_01.ds')\nraw = read_raw_ctf(raw_path, preload=True)"
      ], 
      "outputs": [], 
      "metadata": {
        "collapsed": false
      }
    }, 
    {
      "source": [
        "The sinusoidal signal is generated on channel HDAC006, so we can use\nthat to obtain precise timing.\n\n"
      ], 
      "cell_type": "markdown", 
      "metadata": {}
    }, 
    {
      "execution_count": null, 
      "cell_type": "code", 
      "source": [
        "sinusoid, times = raw[raw.ch_names.index('HDAC006-4408')]\nplt.figure()\nplt.plot(times[times < 1.], sinusoid.T[times < 1.])"
      ], 
      "outputs": [], 
      "metadata": {
        "collapsed": false
      }
    }, 
    {
      "source": [
        "Let's create some events using this signal by thresholding the sinusoid.\n\n"
      ], 
      "cell_type": "markdown", 
      "metadata": {}
    }, 
    {
      "execution_count": null, 
      "cell_type": "code", 
      "source": [
        "events = np.where(np.diff(sinusoid > 0.5) > 0)[1] + raw.first_samp\nevents = np.vstack((events, np.zeros_like(events), np.ones_like(events))).T"
      ], 
      "outputs": [], 
      "metadata": {
        "collapsed": false
      }
    }, 
    {
      "source": [
        "The CTF software compensation works reasonably well:\n\n"
      ], 
      "cell_type": "markdown", 
      "metadata": {}
    }, 
    {
      "execution_count": null, 
      "cell_type": "code", 
      "source": [
        "raw.plot()"
      ], 
      "outputs": [], 
      "metadata": {
        "collapsed": false
      }
    }, 
    {
      "source": [
        "But here we can get slightly better noise suppression, lower localization\nbias, and a better dipole goodness of fit with spatio-temporal (tSSS)\nMaxwell filtering:\n\n"
      ], 
      "cell_type": "markdown", 
      "metadata": {}
    }, 
    {
      "execution_count": null, 
      "cell_type": "code", 
      "source": [
        "raw.apply_gradient_compensation(0)  # must un-do software compensation first\nmf_kwargs = dict(origin=(0., 0., 0.), st_duration=10.)\nraw = mne.preprocessing.maxwell_filter(raw, **mf_kwargs)\nraw.plot()"
      ], 
      "outputs": [], 
      "metadata": {
        "collapsed": false
      }
    }, 
    {
      "source": [
        "Our choice of tmin and tmax should capture exactly one cycle, so\nwe can make the unusual choice of baselining using the entire epoch\nwhen creating our evoked data. We also then crop to a single time point\n(@t=0) because this is a peak in our signal.\n\n"
      ], 
      "cell_type": "markdown", 
      "metadata": {}
    }, 
    {
      "execution_count": null, 
      "cell_type": "code", 
      "source": [
        "tmin = -0.5 / dip_freq\ntmax = -tmin\nepochs = mne.Epochs(raw, events, event_id=1, tmin=tmin, tmax=tmax,\n                    baseline=(None, None))\nevoked = epochs.average()\nevoked.plot()\nevoked.crop(0., 0.)\ndel raw, epochs"
      ], 
      "outputs": [], 
      "metadata": {
        "collapsed": false
      }
    }, 
    {
      "source": [
        "To do a dipole fit, let's use the covariance provided by the empty room\nrecording.\n\n"
      ], 
      "cell_type": "markdown", 
      "metadata": {}
    }, 
    {
      "execution_count": null, 
      "cell_type": "code", 
      "source": [
        "raw_erm = read_raw_ctf(erm_path).apply_gradient_compensation(0)\nraw_erm = mne.preprocessing.maxwell_filter(raw_erm, coord_frame='meg',\n                                           **mf_kwargs)\ncov = mne.compute_raw_covariance(raw_erm)\ndel raw_erm\nsphere = mne.make_sphere_model(r0=(0., 0., 0.), head_radius=None)\ndip = fit_dipole(evoked, cov, sphere)[0]"
      ], 
      "outputs": [], 
      "metadata": {
        "collapsed": false
      }
    }, 
    {
      "source": [
        "Compare the actual position with the estimated one.\n\n"
      ], 
      "cell_type": "markdown", 
      "metadata": {}
    }, 
    {
      "execution_count": null, 
      "cell_type": "code", 
      "source": [
        "expected_pos = np.array([18., 0., 49.])\ndiff = np.sqrt(np.sum((dip.pos[0] * 1000 - expected_pos) ** 2))\nprint('Actual pos:     %s mm' % np.array_str(expected_pos, precision=1))\nprint('Estimated pos:  %s mm' % np.array_str(dip.pos[0] * 1000, precision=1))\nprint('Difference:     %0.1f mm' % diff)\nprint('Amplitude:      %0.1f nAm' % (1e9 * dip.amplitude[0]))\nprint('GOF:            %0.1f %%' % dip.gof[0])"
      ], 
      "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"
      }
    }
  }
}