{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n\n# Compute source level time-frequency timecourses using a DICS beamformer\n\nIn this example, a Dynamic Imaging of Coherent Sources (DICS)\n:footcite:`GrossEtAl2001` beamformer is used to transform sensor-level\ntime-frequency objects to the source level. We will look at the event-related\nsynchronization (ERS) of beta band activity in the `somato dataset\n<somato-dataset>`.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Authors: Marijn van Vliet <w.m.vanvliet@gmail.com>\n#          Alex Rockhill <aprockhill@mailbox.org>\n#\n# License: BSD-3-Clause\n# Copyright the MNE-Python contributors.\n\nimport numpy as np\n\nimport mne\nfrom mne.beamformer import apply_dics_tfr_epochs, make_dics\nfrom mne.datasets import somato\nfrom mne.time_frequency import csd_tfr\n\nprint(__doc__)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Organize the data that we will use for this example.\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\"\nfname_fwd = (\n    data_path / \"derivatives\" / f\"sub-{subject}\" / f\"sub-{subject}_task-{task}-fwd.fif\"\n)\nsubjects_dir = data_path / \"derivatives\" / \"freesurfer\" / \"subjects\""
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "First, we load the data and compute for each epoch the time-frequency\ndecomposition in sensor space.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Load raw data and make epochs.\nraw = mne.io.read_raw_fif(raw_fname)\nevents = mne.find_events(raw)\nepochs = mne.Epochs(\n    raw,\n    events[:22],  # just for execution speed of the tutorial\n    event_id=1,\n    tmin=-1,\n    tmax=2.5,\n    reject=dict(\n        grad=5000e-13,  # unit: T / m (gradiometers)\n        mag=5e-12,  # unit: T (magnetometers)\n        eog=250e-6,  # unit: V (EOG channels)\n    ),\n    preload=True,\n)\n\n# We are mostly interested in the beta band since it has been shown to be\n# active for somatosensory stimulation\nfreqs = np.linspace(13, 31, 5)\n\n# Use Morlet wavelets to compute sensor-level time-frequency (TFR)\n# decomposition for each epoch. We must pass ``output='complex'`` if we wish to\n# use this TFR later with a DICS beamformer. We also pass ``average=False`` to\n# compute the TFR for each individual epoch.\nepochs_tfr = epochs.compute_tfr(\n    \"morlet\", freqs, n_cycles=5, return_itc=False, output=\"complex\", average=False\n)\n\n# crop either side to use a buffer to remove edge artifact\nepochs_tfr.crop(tmin=-0.5, tmax=2)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Now, we build a DICS beamformer and project the sensor-level TFR to the\nsource level.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Compute the Cross-Spectral Density (CSD) matrix for the sensor-level TFRs.\n# We are interested in increases in power relative to the baseline period, so\n# we will make a separate CSD for just that period as well.\ncsd = csd_tfr(epochs_tfr, tmin=-0.5, tmax=2)\nbaseline_csd = csd_tfr(epochs_tfr, tmin=-0.5, tmax=-0.1)\n\n# use the CSDs and the forward model to build the DICS beamformer\nfwd = mne.read_forward_solution(fname_fwd)\n\n# compute scalar DICS beamfomer\nfilters = make_dics(\n    epochs.info,\n    fwd,\n    csd,\n    noise_csd=baseline_csd,\n    pick_ori=\"max-power\",\n    reduce_rank=True,\n    real_filter=True,\n)\n\n# project the TFR for each epoch to source space\nepochs_stcs = apply_dics_tfr_epochs(epochs_tfr, filters, return_generator=True)\n\n# average across frequencies and epochs\ndata = np.zeros((fwd[\"nsource\"], epochs_tfr.times.size))\nfor epoch_stcs in epochs_stcs:\n    for stc in epoch_stcs:\n        data += (stc.data * np.conj(stc.data)).real\n\nstc.data = data / len(epochs) / len(freqs)\n\n# apply a baseline correction\nstc.apply_baseline((-0.5, -0.1))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Let's visualize the source time course estimate. We can see the\nexpected activation of the two gyri bordering the central sulcus, the\nprimary somatosensory and motor cortices (S1 and M1).\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "fmax = 4500\nbrain = stc.plot(\n    subjects_dir=subjects_dir,\n    hemi=\"both\",\n    views=\"dorsal\",\n    initial_time=1.2,\n    brain_kwargs=dict(show=False),\n    add_data_kwargs=dict(\n        fmin=fmax / 10,\n        fmid=fmax / 2,\n        fmax=fmax,\n        scale_factor=0.0001,\n        colorbar_kwargs=dict(label_font_size=10),\n    ),\n)"
      ]
    }
  ],
  "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
}