{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "\n\n# Brainstorm CTF phantom dataset tutorial\n\nHere we compute the evoked from raw for the Brainstorm CTF phantom\ntutorial dataset. For comparison, see :footcite:`TadelEtAl2011` and:\n\n https://neuroimage.usc.edu/brainstorm/Tutorials/PhantomCtf\n\n## References\n.. footbibliography::\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Authors: Eric Larson \n#\n# License: BSD-3-Clause\n# Copyright the MNE-Python contributors." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import warnings\n\nimport matplotlib.pyplot as plt\nimport numpy as np\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__)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The data were collected with a CTF system at 2400 Hz.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "data_path = bst_phantom_ctf.data_path(verbose=True)\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 = data_path / \"phantom_20uA_20150603_03.ds\"\ndip_freq = 23.0\nerm_path = data_path / \"emptyroom_20150709_01.ds\"\nraw = read_raw_ctf(raw_path, preload=True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The sinusoidal signal is generated on channel HDAC006, so we can use\nthat to obtain precise timing.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "sinusoid, times = raw[raw.ch_names.index(\"HDAC006-4408\")]\nplt.figure()\nplt.plot(times[times < 1.0], sinusoid.T[times < 1.0])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's create some events using this signal by thresholding the sinusoid.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "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" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The CTF software compensation works reasonably well:\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "raw.plot()" ] }, { "cell_type": "markdown", "metadata": {}, "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": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "raw.apply_gradient_compensation(0) # must un-do software compensation first\nmf_kwargs = dict(origin=(0.0, 0.0, 0.0), st_duration=10.0)\nraw = mne.preprocessing.maxwell_filter(raw, **mf_kwargs)\nraw.plot()" ] }, { "cell_type": "markdown", "metadata": {}, "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": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "tmin = -0.5 / dip_freq\ntmax = -tmin\nepochs = mne.Epochs(\n raw, events, event_id=1, tmin=tmin, tmax=tmax, baseline=(None, None)\n)\nevoked = epochs.average()\nevoked.plot(time_unit=\"s\")\nevoked.crop(0.0, 0.0)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\nLet's use a `sphere head geometry model `\nand let's see the coordinate alignment and the sphere location.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "sphere = mne.make_sphere_model(r0=(0.0, 0.0, 0.0), head_radius=0.08)\n\nmne.viz.plot_alignment(\n raw.info, subject=\"sample\", meg=\"helmet\", bem=sphere, dig=True, surfaces=[\"brain\"]\n)\ndel raw, epochs" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To do a dipole fit, let's use the covariance provided by the empty room\nrecording.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "raw_erm = read_raw_ctf(erm_path).apply_gradient_compensation(0)\nraw_erm = mne.preprocessing.maxwell_filter(raw_erm, coord_frame=\"meg\", **mf_kwargs)\ncov = mne.compute_raw_covariance(raw_erm)\ndel raw_erm\n\nwith warnings.catch_warnings(record=True):\n # ignore warning about data rank exceeding that of info (75 > 71)\n warnings.simplefilter(\"ignore\")\n dip, residual = fit_dipole(evoked, cov, sphere, verbose=True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Compare the actual position with the estimated one.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "expected_pos = np.array([18.0, 0.0, 49.0])\ndiff = np.sqrt(np.sum((dip.pos[0] * 1000 - expected_pos) ** 2))\nprint(f\"Actual pos: {np.array_str(expected_pos, precision=1)} mm\")\nprint(f\"Estimated pos: {np.array_str(dip.pos[0] * 1000, precision=1)} mm\")\nprint(f\"Difference: {diff:0.1f} mm\")\nprint(f\"Amplitude: {1e9 * dip.amplitude[0]:0.1f} nAm\")\nprint(f\"GOF: {dip.gof[0]:0.1f} %\")" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.12.2" } }, "nbformat": 4, "nbformat_minor": 0 }