{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "\n\n# Explore event-related dynamics for specific frequency bands\n\nThe objective is to show you how to explore spectrally localized\neffects. For this purpose we adapt the method described in\n:footcite:`HariSalmelin1997` and use it on the somato dataset.\nThe idea is to track the band-limited temporal evolution\nof spatial patterns by using the :term:`global field power` (GFP).\n\nWe first bandpass filter the signals and then apply a Hilbert transform. To\nreveal oscillatory activity the evoked response is then subtracted from every\nsingle trial. Finally, we rectify the signals prior to averaging across trials\nby taking the magnitude of the Hilbert.\nThen the :term:`GFP` is computed as described in\n:footcite:`EngemannGramfort2015`, using the sum of the\nsquares but without normalization by the rank.\nBaselining is subsequently applied to make the :term:`GFP` comparable\nbetween frequencies.\nThe procedure is then repeated for each frequency band of interest and\nall :term:`GFPs` are visualized. To estimate uncertainty, non-parametric\nconfidence intervals are computed as described in :footcite:`EfronHastie2016`\nacross channels.\n\nThe advantage of this method over summarizing the Space \u00d7 Time \u00d7 Frequency\noutput of a Morlet Wavelet in frequency bands is relative speed and, more\nimportantly, the clear-cut comparability of the spectral decomposition (the\nsame type of filter is used across all bands).\n\nWe will use this dataset: `somato-dataset`\n\n## References\n.. footbibliography::\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Authors: Denis A. Engemann \n# Stefan Appelhoff \n#\n# License: BSD-3-Clause\n# Copyright the MNE-Python contributors." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import matplotlib.pyplot as plt\nimport numpy as np\n\nimport mne\nfrom mne.baseline import rescale\nfrom mne.datasets import somato\nfrom mne.stats import bootstrap_confidence_interval" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Set parameters\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "data_path = somato.data_path()\nsubject = \"01\"\ntask = \"somato\"\nraw_fname = data_path / f\"sub-{subject}\" / \"meg\" / f\"sub-{subject}_task-{task}_meg.fif\"\n\n# let's explore some frequency bands\niter_freqs = [(\"Theta\", 4, 7), (\"Alpha\", 8, 12), (\"Beta\", 13, 25), (\"Gamma\", 30, 45)]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We create average power time courses for each frequency band\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# set epoching parameters\nevent_id, tmin, tmax = 1, -1.0, 3.0\nbaseline = None\n\n# get the header to extract events\nraw = mne.io.read_raw_fif(raw_fname)\nevents = mne.find_events(raw, stim_channel=\"STI 014\")\n\nfrequency_map = list()\n\nfor band, fmin, fmax in iter_freqs:\n # (re)load the data to save memory\n raw = mne.io.read_raw_fif(raw_fname)\n raw.pick(picks=[\"grad\", \"eog\"]) # we just look at gradiometers\n raw.load_data()\n\n # bandpass filter\n raw.filter(\n fmin,\n fmax,\n n_jobs=None, # use more jobs to speed up.\n l_trans_bandwidth=1, # make sure filter params are the same\n h_trans_bandwidth=1,\n ) # in each band and skip \"auto\" option.\n\n # epoch\n epochs = mne.Epochs(\n raw,\n events,\n event_id,\n tmin,\n tmax,\n baseline=baseline,\n reject=dict(grad=4000e-13, eog=350e-6),\n preload=True,\n )\n # remove evoked response\n epochs.subtract_evoked()\n\n # get analytic signal (envelope)\n epochs.apply_hilbert(envelope=True)\n frequency_map.append(((band, fmin, fmax), epochs.average()))\n del epochs\ndel raw" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we can compute the Global Field Power\nWe can track the emergence of spatial patterns compared to baseline\nfor each frequency band, with a bootstrapped confidence interval.\n\nWe see dominant responses in the Alpha and Beta bands.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Helper function for plotting spread\ndef stat_fun(x):\n \"\"\"Return sum of squares.\"\"\"\n return np.sum(x**2, axis=0)\n\n\n# Plot\nfig, axes = plt.subplots(4, 1, figsize=(10, 7), sharex=True, sharey=True)\ncolors = plt.colormaps[\"winter_r\"](np.linspace(0, 1, 4))\nfor ((freq_name, fmin, fmax), average), color, ax in zip(\n frequency_map, colors, axes.ravel()[::-1]\n):\n times = average.times * 1e3\n gfp = np.sum(average.data**2, axis=0)\n gfp = mne.baseline.rescale(gfp, times, baseline=(None, 0))\n ax.plot(times, gfp, label=freq_name, color=color, linewidth=2.5)\n ax.axhline(0, linestyle=\"--\", color=\"grey\", linewidth=2)\n ci_low, ci_up = bootstrap_confidence_interval(\n average.data, random_state=0, stat_fun=stat_fun\n )\n ci_low = rescale(ci_low, average.times, baseline=(None, 0))\n ci_up = rescale(ci_up, average.times, baseline=(None, 0))\n ax.fill_between(times, gfp + ci_up, gfp - ci_low, color=color, alpha=0.3)\n ax.grid(True)\n ax.set_ylabel(\"GFP\")\n ax.annotate(\n f\"{freq_name} ({fmin:d}-{fmax:d}Hz)\",\n xy=(0.95, 0.8),\n horizontalalignment=\"right\",\n xycoords=\"axes fraction\",\n )\n ax.set_xlim(-1000, 3000)\n\naxes.ravel()[-1].set_xlabel(\"Time [ms]\")" ] } ], "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 }