{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n\n# Compute source power spectral density (PSD) of VectorView and OPM data\n\nHere we compute the resting state from raw for data recorded using\na Neuromag VectorView system and a custom OPM system.\nThe pipeline is meant to mostly follow the Brainstorm :footcite:`TadelEtAl2011`\n[OMEGA resting tutorial pipeline](https://neuroimage.usc.edu/brainstorm/Tutorials/RestingOmega)_.\nThe steps we use are:\n\n1. Filtering: downsample heavily.\n2. Artifact detection: use SSP for EOG and ECG.\n3. Source localization: dSPM, depth weighting, cortically constrained.\n4. Frequency: power spectral density (Welch), 4 s window, 50% overlap.\n5. Standardize: normalize by relative power for each source.\n\n## Preprocessing\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Authors: Denis Engemann <denis.engemann@gmail.com>\n#          Luke Bloy <luke.bloy@gmail.com>\n#          Eric Larson <larson.eric.d@gmail.com>\n#\n# License: BSD-3-Clause\n# Copyright the MNE-Python contributors."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import mne\nfrom mne.filter import next_fast_len\n\nprint(__doc__)\n\ndata_path = mne.datasets.opm.data_path()\nsubject = \"OPM_sample\"\n\nsubjects_dir = data_path / \"subjects\"\nbem_dir = subjects_dir / subject / \"bem\"\nbem_fname = bem_dir / f\"{subject}-5120-5120-5120-bem-sol.fif\"\nsrc_fname = bem_dir / f\"{subject}-oct6-src.fif\"\nvv_fname = data_path / \"MEG\" / \"SQUID\" / \"SQUID_resting_state.fif\"\nvv_erm_fname = data_path / \"MEG\" / \"SQUID\" / \"SQUID_empty_room.fif\"\nvv_trans_fname = data_path / \"MEG\" / \"SQUID\" / \"SQUID-trans.fif\"\nopm_fname = data_path / \"MEG\" / \"OPM\" / \"OPM_resting_state_raw.fif\"\nopm_erm_fname = data_path / \"MEG\" / \"OPM\" / \"OPM_empty_room_raw.fif\"\nopm_trans = mne.transforms.Transform(\"head\", \"mri\")  # use identity transform\nopm_coil_def_fname = data_path / \"MEG\" / \"OPM\" / \"coil_def.dat\""
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Load data, resample. We will store the raw objects in dicts with entries\n\"vv\" and \"opm\" to simplify housekeeping and simplify looping later.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "raws = dict()\nraw_erms = dict()\nnew_sfreq = 60.0  # Nyquist frequency (30 Hz) < line noise freq (50 Hz)\nraws[\"vv\"] = mne.io.read_raw_fif(vv_fname, verbose=\"error\")  # ignore naming\nraws[\"vv\"].load_data().resample(new_sfreq, method=\"polyphase\")\nraws[\"vv\"].info[\"bads\"] = [\"MEG2233\", \"MEG1842\"]\nraw_erms[\"vv\"] = mne.io.read_raw_fif(vv_erm_fname, verbose=\"error\")\nraw_erms[\"vv\"].load_data().resample(new_sfreq, method=\"polyphase\")\nraw_erms[\"vv\"].info[\"bads\"] = [\"MEG2233\", \"MEG1842\"]\n\nraws[\"opm\"] = mne.io.read_raw_fif(opm_fname)\nraws[\"opm\"].load_data().resample(new_sfreq, method=\"polyphase\")\nraw_erms[\"opm\"] = mne.io.read_raw_fif(opm_erm_fname)\nraw_erms[\"opm\"].load_data().resample(new_sfreq, method=\"polyphase\")\n# Make sure our assumptions later hold\nassert raws[\"opm\"].info[\"sfreq\"] == raws[\"vv\"].info[\"sfreq\"]"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Explore data\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "titles = dict(vv=\"VectorView\", opm=\"OPM\")\nkinds = (\"vv\", \"opm\")\nn_fft = next_fast_len(int(round(4 * new_sfreq)))\nprint(f\"Using n_fft={n_fft} ({n_fft / raws['vv'].info['sfreq']:0.1f} s)\")\nfor kind in kinds:\n    fig = (\n        raws[kind]\n        .compute_psd(n_fft=n_fft, proj=True)\n        .plot(picks=\"data\", exclude=\"bads\", amplitude=True)\n    )\n    fig.suptitle(titles[kind])"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Alignment and forward\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Here we use a reduced size source space (oct5) just for speed\nsrc = mne.setup_source_space(subject, \"oct5\", add_dist=False, subjects_dir=subjects_dir)\n# This line removes source-to-source distances that we will not need.\n# We only do it here to save a bit of memory, in general this is not required.\ndel src[0][\"dist\"], src[1][\"dist\"]\nbem = mne.read_bem_solution(bem_fname)\n# For speed, let's just use a 1-layer BEM\nbem = mne.make_bem_solution(bem[\"surfs\"][-1:])\nfwd = dict()\n\n# check alignment and generate forward for VectorView\nkwargs = dict(azimuth=0, elevation=90, distance=0.6, focalpoint=(0.0, 0.0, 0.0))\nfig = mne.viz.plot_alignment(\n    raws[\"vv\"].info,\n    trans=vv_trans_fname,\n    subject=subject,\n    subjects_dir=subjects_dir,\n    dig=True,\n    coord_frame=\"mri\",\n    surfaces=(\"head\", \"white\"),\n)\nmne.viz.set_3d_view(figure=fig, **kwargs)\nfwd[\"vv\"] = mne.make_forward_solution(\n    raws[\"vv\"].info, vv_trans_fname, src, bem, eeg=False, verbose=True\n)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "And for OPM:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "with mne.use_coil_def(opm_coil_def_fname):\n    fig = mne.viz.plot_alignment(\n        raws[\"opm\"].info,\n        trans=opm_trans,\n        subject=subject,\n        subjects_dir=subjects_dir,\n        dig=False,\n        coord_frame=\"mri\",\n        surfaces=(\"head\", \"white\"),\n    )\n    mne.viz.set_3d_view(figure=fig, **kwargs)\n    fwd[\"opm\"] = mne.make_forward_solution(\n        raws[\"opm\"].info, opm_trans, src, bem, eeg=False, verbose=True\n    )\n\ndel src, bem"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Compute and apply inverse to PSD estimated using multitaper + Welch.\nGroup into frequency bands, then normalize each source point and sensor\nindependently. This makes the value of each sensor point and source location\nin each frequency band the percentage of the PSD accounted for by that band.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "freq_bands = dict(alpha=(8, 12), beta=(15, 29))\ntopos = dict(vv=dict(), opm=dict())\nstcs = dict(vv=dict(), opm=dict())\n\nsnr = 3.0\nlambda2 = 1.0 / snr**2\nfor kind in kinds:\n    noise_cov = mne.compute_raw_covariance(raw_erms[kind])\n    inverse_operator = mne.minimum_norm.make_inverse_operator(\n        raws[kind].info, forward=fwd[kind], noise_cov=noise_cov, verbose=True\n    )\n    stc_psd, sensor_psd = mne.minimum_norm.compute_source_psd(\n        raws[kind],\n        inverse_operator,\n        lambda2=lambda2,\n        n_fft=n_fft,\n        dB=False,\n        return_sensor=True,\n        verbose=True,\n    )\n    topo_norm = sensor_psd.data.sum(axis=1, keepdims=True)\n    stc_norm = stc_psd.sum()  # same operation on MNE object, sum across freqs\n    # Normalize each source point by the total power across freqs\n    for band, limits in freq_bands.items():\n        data = sensor_psd.copy().crop(*limits).data.sum(axis=1, keepdims=True)\n        topos[kind][band] = mne.EvokedArray(100 * data / topo_norm, sensor_psd.info)\n        stcs[kind][band] = 100 * stc_psd.copy().crop(*limits).sum() / stc_norm.data\n    del inverse_operator\ndel fwd, raws, raw_erms"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Now we can make some plots of each frequency band. Note that the OPM head\ncoverage is only over right motor cortex, so only localization\nof beta is likely to be worthwhile.\n\n## Alpha\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "def plot_band(kind, band):\n    \"\"\"Plot activity within a frequency band on the subject's brain.\"\"\"\n    lf, hf = freq_bands[band]\n    title = f\"{titles[kind]} {band}\\n({lf:d}-{hf:d} Hz)\"\n    topos[kind][band].plot_topomap(\n        times=0.0,\n        scalings=1.0,\n        cbar_fmt=\"%0.1f\",\n        vlim=(0, None),\n        cmap=\"inferno\",\n        time_format=title,\n    )\n    brain = stcs[kind][band].plot(\n        subject=subject,\n        subjects_dir=subjects_dir,\n        views=\"cau\",\n        hemi=\"both\",\n        time_label=title,\n        title=title,\n        colormap=\"inferno\",\n        time_viewer=False,\n        show_traces=False,\n        clim=dict(kind=\"percent\", lims=(70, 85, 99)),\n        smoothing_steps=10,\n    )\n    brain.show_view(azimuth=0, elevation=0, roll=0)\n    return fig, brain\n\n\nfig_alpha, brain_alpha = plot_band(\"vv\", \"alpha\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Beta\nHere we also show OPM data, which shows a profile similar to the VectorView\ndata beneath the sensors. VectorView first:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "fig_beta, brain_beta = plot_band(\"vv\", \"beta\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Then OPM:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "fig_beta_opm, brain_beta_opm = plot_band(\"opm\", \"beta\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## References\n.. footbibliography::\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
}