<jupyter_start><jupyter_text>Brainstorm Elekta phantom dataset tutorialHere we compute the evoked from raw for the Brainstorm Elekta phantomtutorial dataset. For comparison, see [1]_ and:    http://neuroimage.usc.edu/brainstorm/Tutorials/PhantomElektaReferences----------.. [1] <NAME>, <NAME>, <NAME>, <NAME>, <NAME>.       Brainstorm: A User-Friendly Application for MEG/EEG Analysis.       Computational Intelligence and Neuroscience, vol. 2011, Article ID       879716, 13 pages, 2011. doi:10.1155/2011/879716<jupyter_code># sphinx_gallery_thumbnail_number = 9

# Authors: <NAME> <<EMAIL>>
#
# License: BSD (3-clause)

import os.path as op
import numpy as np
import matplotlib.pyplot as plt

import mne
from mne import find_events, fit_dipole
from mne.datasets.brainstorm import bst_phantom_elekta
from mne.io import read_raw_fif

from mayavi import mlab
print(__doc__)<jupyter_output><empty_output><jupyter_text>The data were collected with an Elekta Neuromag VectorView system at 1000 Hzand low-pass filtered at 330 Hz. Here the medium-amplitude (200 nAm) dataare read to construct instances of :class:`mne.io.Raw`.<jupyter_code>data_path = bst_phantom_elekta.data_path(verbose=True)

raw_fname = op.join(data_path, 'kojak_all_200nAm_pp_no_chpi_no_ms_raw.fif')
raw = read_raw_fif(raw_fname)<jupyter_output><empty_output><jupyter_text>Data channel array consisted of 204 MEG planor gradiometers,102 axial magnetometers, and 3 stimulus channels. Let's get the eventsfor the phantom, where each dipole (1-32) gets its own event:<jupyter_code>events = find_events(raw, 'STI201')
raw.plot(events=events)
raw.info['bads'] = ['MEG2421']<jupyter_output><empty_output><jupyter_text>The data have strong line frequency (60 Hz and harmonics) and cHPI coilnoise (five peaks around 300 Hz). Here we plot only out to 60 secondsto save memory:<jupyter_code>raw.plot_psd(tmax=60., average=False)<jupyter_output><empty_output><jupyter_text>Let's use Maxwell filtering to clean the data a bit.Ideally we would have the fine calibration and cross-talk informationfor the site of interest, but we don't, so we just do:<jupyter_code>raw.fix_mag_coil_types()
raw = mne.preprocessing.maxwell_filter(raw, origin=(0., 0., 0.))<jupyter_output><empty_output><jupyter_text>We know our phantom produces sinusoidal bursts below 25 Hz, so let's filter.<jupyter_code>raw.filter(None, 40., fir_design='firwin')
raw.plot(events=events)<jupyter_output><empty_output><jupyter_text>Now we epoch our data, average it, and look at the first dipole response.The first peak appears around 3 ms. Because we low-passed at 40 Hz,we can also decimate our data to save memory.<jupyter_code>tmin, tmax = -0.1, 0.1
event_id = list(range(1, 33))
epochs = mne.Epochs(raw, events, event_id, tmin, tmax, baseline=(None, -0.01),
                    decim=3, preload=True)
epochs['1'].average().plot(time_unit='s')<jupyter_output><empty_output><jupyter_text>Let's use a sphere head geometry model and let's see the coordinatealignment and the sphere location. The phantom is properly modeled bya single-shell sphere with origin (0., 0., 0.).<jupyter_code>sphere = mne.make_sphere_model(r0=(0., 0., 0.), head_radius=0.08)

mne.viz.plot_alignment(raw.info, subject='sample', show_axes=True,
                       bem=sphere, dig=True, surfaces='inner_skull')<jupyter_output><empty_output><jupyter_text>Let's do some dipole fits. We first compute the noise covariance,then do the fits for each event_id taking the time instant that maximizesthe global field power.<jupyter_code># here we can get away with using method='oas' for speed (faster than "shrunk")
# but in general "shrunk" is usually better
cov = mne.compute_covariance(
    epochs, tmax=0, method='oas', rank=None)
mne.viz.plot_evoked_white(epochs['1'].average(), cov)

data = []
t_peak = 0.036  # true for Elekta phantom
for ii in event_id:
    evoked = epochs[str(ii)].average().crop(t_peak, t_peak)
    data.append(evoked.data[:, 0])
evoked = mne.EvokedArray(np.array(data).T, evoked.info, tmin=0.)
del epochs, raw
dip, residual = fit_dipole(evoked, cov, sphere, n_jobs=1)<jupyter_output><empty_output><jupyter_text>Do a quick visualization of how much variance we explained, putting thedata and residuals on the same scale (here the "time points" are the32 dipole peak values that we fit):<jupyter_code>fig, axes = plt.subplots(2, 1)
evoked.plot(axes=axes)
for ax in axes:
    ax.texts = []
    for line in ax.lines:
        line.set_color('#98df81')
residual.plot(axes=axes)<jupyter_output><empty_output><jupyter_text>Now we can compare to the actual locations, taking the difference in mm:<jupyter_code>actual_pos, actual_ori = mne.dipole.get_phantom_dipoles()
actual_amp = 100.  # nAm

fig, (ax1, ax2, ax3) = plt.subplots(nrows=3, ncols=1, figsize=(6, 7))

diffs = 1000 * np.sqrt(np.sum((dip.pos - actual_pos) ** 2, axis=-1))
print('mean(position error) = %0.1f mm' % (np.mean(diffs),))
ax1.bar(event_id, diffs)
ax1.set_xlabel('Dipole index')
ax1.set_ylabel('Loc. error (mm)')

angles = np.rad2deg(np.arccos(np.abs(np.sum(dip.ori * actual_ori, axis=1))))
print(u'mean(angle error) = %0.1f°' % (np.mean(angles),))
ax2.bar(event_id, angles)
ax2.set_xlabel('Dipole index')
ax2.set_ylabel(u'Angle error (°)')

amps = actual_amp - dip.amplitude / 1e-9
print('mean(abs amplitude error) = %0.1f nAm' % (np.mean(np.abs(amps)),))
ax3.bar(event_id, amps)
ax3.set_xlabel('Dipole index')
ax3.set_ylabel('Amplitude error (nAm)')

fig.tight_layout()
plt.show()<jupyter_output><empty_output><jupyter_text>Let's plot the positions and the orientations of the actual and the estimateddipoles<jupyter_code>def plot_pos_ori(pos, ori, color=(0., 0., 0.), opacity=1.):
    x, y, z = pos.T
    u, v, w = ori.T
    mlab.points3d(x, y, z, scale_factor=0.005, opacity=opacity, color=color)
    q = mlab.quiver3d(x, y, z, u, v, w,
                      scale_factor=0.03, opacity=opacity,
                      color=color, mode='arrow')
    q.glyph.glyph_source.glyph_source.shaft_radius = 0.02
    q.glyph.glyph_source.glyph_source.tip_length = 0.1
    q.glyph.glyph_source.glyph_source.tip_radius = 0.05


mne.viz.plot_alignment(evoked.info, bem=sphere, surfaces='inner_skull',
                       coord_frame='head', meg='helmet', show_axes=True)

# Plot the position and the orientation of the actual dipole
plot_pos_ori(actual_pos, actual_ori, color=(0., 0., 0.), opacity=0.5)
# Plot the position and the orientation of the estimated dipole
plot_pos_ori(dip.pos, dip.ori, color=(0.2, 1., 0.5))
mlab.view(70, 80, distance=0.5)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.12431360585723002, 'fraction_numerical': 0.02989627821842587, 'mean_word_length': 5.197542533081285, 'pattern_counts': {'":': 0, '<': 57, '<?xml version=': 0, '>': 57, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '12', 'dir': 'jupyter-structured-clean-dedup', 'id': '196198', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 2552, 'n_tokens_neox': 2316, 'n_words': 614, 'path': '0.17/_downloads/2aba6a5c9f79fe16cdce1a232bc5e327/plot_brainstorm_phantom_elekta.ipynb', 'repo_name': 'drammock/mne-tools.github.io', 'size': '6590'}