{
  "nbformat_minor": 0, 
  "nbformat": 4, 
  "cells": [
    {
      "execution_count": null, 
      "cell_type": "code", 
      "source": [
        "%matplotlib inline"
      ], 
      "outputs": [], 
      "metadata": {
        "collapsed": false
      }
    }, 
    {
      "source": [
        "\n# Brainstorm Elekta phantom tutorial dataset\n\n\nHere we compute the evoked from raw for the Brainstorm Elekta phantom\ntutorial dataset. For comparison, see [1]_ and:\n\n    http://neuroimage.usc.edu/brainstorm/Tutorials/PhantomElekta\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\n\nimport mne\nfrom mne import find_events, fit_dipole\nfrom mne.datasets.brainstorm import bst_phantom_elekta\nfrom mne.io import read_raw_fif\n\nprint(__doc__)"
      ], 
      "outputs": [], 
      "metadata": {
        "collapsed": false
      }
    }, 
    {
      "source": [
        "The data were collected with an Elekta Neuromag VectorView system at 1000 Hz\nand low-pass filtered at 330 Hz. Here the medium-amplitude (200 nAm) data\nare read to construct instances of :class:`mne.io.Raw`.\n\n"
      ], 
      "cell_type": "markdown", 
      "metadata": {}
    }, 
    {
      "execution_count": null, 
      "cell_type": "code", 
      "source": [
        "data_path = bst_phantom_elekta.data_path()\n\nraw_fname = op.join(data_path, 'kojak_all_200nAm_pp_no_chpi_no_ms_raw.fif')\nraw = read_raw_fif(raw_fname)"
      ], 
      "outputs": [], 
      "metadata": {
        "collapsed": false
      }
    }, 
    {
      "source": [
        "Data channel array consisted of 204 MEG planor gradiometers,\n102 axial magnetometers, and 3 stimulus channels. Let's get the events\nfor the phantom, where each dipole (1-32) gets its own event:\n\n"
      ], 
      "cell_type": "markdown", 
      "metadata": {}
    }, 
    {
      "execution_count": null, 
      "cell_type": "code", 
      "source": [
        "events = find_events(raw, 'STI201')\nraw.plot(events=events)\nraw.info['bads'] = ['MEG2421']"
      ], 
      "outputs": [], 
      "metadata": {
        "collapsed": false
      }
    }, 
    {
      "source": [
        "The data have strong line frequency (60 Hz and harmonics) and cHPI coil\nnoise (five peaks around 300 Hz). Here we plot only out to 60 seconds\nto save memory:\n\n"
      ], 
      "cell_type": "markdown", 
      "metadata": {}
    }, 
    {
      "execution_count": null, 
      "cell_type": "code", 
      "source": [
        "raw.plot_psd(tmax=60.)"
      ], 
      "outputs": [], 
      "metadata": {
        "collapsed": false
      }
    }, 
    {
      "source": [
        "Let's use Maxwell filtering to clean the data a bit.\nIdeally we would have the fine calibration and cross-talk information\nfor the site of interest, but we don't, so we just do:\n\n"
      ], 
      "cell_type": "markdown", 
      "metadata": {}
    }, 
    {
      "execution_count": null, 
      "cell_type": "code", 
      "source": [
        "raw.fix_mag_coil_types()\nraw = mne.preprocessing.maxwell_filter(raw, origin=(0., 0., 0.))"
      ], 
      "outputs": [], 
      "metadata": {
        "collapsed": false
      }
    }, 
    {
      "source": [
        "We know our phantom produces sinusoidal bursts below 25 Hz, so let's filter.\n\n"
      ], 
      "cell_type": "markdown", 
      "metadata": {}
    }, 
    {
      "execution_count": null, 
      "cell_type": "code", 
      "source": [
        "raw.filter(None, 40., h_trans_bandwidth='auto', filter_length='auto',\n           phase='zero', fir_window='hamming')\nraw.plot(events=events)"
      ], 
      "outputs": [], 
      "metadata": {
        "collapsed": false
      }
    }, 
    {
      "source": [
        "Now we epoch our data, average it, and look at the first dipole response.\nThe first peak appears around 3 ms. Because we low-passed at 40 Hz,\nwe can also decimate our data to save memory.\n\n"
      ], 
      "cell_type": "markdown", 
      "metadata": {}
    }, 
    {
      "execution_count": null, 
      "cell_type": "code", 
      "source": [
        "tmin, tmax = -0.1, 0.1\nevent_id = list(range(1, 33))\nepochs = mne.Epochs(raw, events, event_id, tmin, tmax, baseline=(None, -0.01),\n                    decim=5, preload=True)\nepochs['1'].average().plot()"
      ], 
      "outputs": [], 
      "metadata": {
        "collapsed": false
      }
    }, 
    {
      "source": [
        "Let's do some dipole fits. The phantom is properly modeled by a single-shell\nsphere with origin (0., 0., 0.). We compute covariance, then do the fits.\n\n"
      ], 
      "cell_type": "markdown", 
      "metadata": {}
    }, 
    {
      "execution_count": null, 
      "cell_type": "code", 
      "source": [
        "t_peak = 60e-3  # ~60 MS at largest peak\nsphere = mne.make_sphere_model(r0=(0., 0., 0.), head_radius=None)\ncov = mne.compute_covariance(epochs, tmax=0)\ndata = []\nfor ii in range(1, 33):\n    evoked = epochs[str(ii)].average().crop(t_peak, t_peak)\n    data.append(evoked.data[:, 0])\nevoked = mne.EvokedArray(np.array(data).T, evoked.info, tmin=0.)\ndel epochs, raw\ndip = fit_dipole(evoked, cov, sphere, n_jobs=1)[0]"
      ], 
      "outputs": [], 
      "metadata": {
        "collapsed": false
      }
    }, 
    {
      "source": [
        "Now we can compare to the actual locations, taking the difference in mm:\n\n"
      ], 
      "cell_type": "markdown", 
      "metadata": {}
    }, 
    {
      "execution_count": null, 
      "cell_type": "code", 
      "source": [
        "actual_pos = mne.dipole.get_phantom_dipoles()[0]\ndiffs = 1000 * np.sqrt(np.sum((dip.pos - actual_pos) ** 2, axis=-1))\nprint('Differences (mm):\\n%s' % diffs[:, np.newaxis])\nprint('\u03bc = %s' % (np.mean(diffs),))"
      ], 
      "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"
      }
    }
  }
}