{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n\n# Compute Power Spectral Density of inverse solution from single epochs\n\nCompute PSD of dSPM inverse solution on single trial epochs restricted\nto a brain label. The PSD is computed using a multi-taper method with\nDiscrete Prolate Spheroidal Sequence (DPSS) windows.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Author: Martin Luessi <mluessi@nmr.mgh.harvard.edu>\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\n\nimport mne\nfrom mne.datasets import sample\nfrom mne.minimum_norm import compute_source_psd_epochs, read_inverse_operator\n\nprint(__doc__)\n\ndata_path = sample.data_path()\nmeg_path = data_path / \"MEG\" / \"sample\"\nfname_inv = meg_path / \"sample_audvis-meg-oct-6-meg-inv.fif\"\nfname_raw = meg_path / \"sample_audvis_raw.fif\"\nfname_event = meg_path / \"sample_audvis_raw-eve.fif\"\nlabel_name = \"Aud-lh\"\nfname_label = meg_path / \"labels\" / f\"{label_name}.label\"\nsubjects_dir = data_path / \"subjects\"\n\nevent_id, tmin, tmax = 1, -0.2, 0.5\nsnr = 1.0  # use smaller SNR for raw data\nlambda2 = 1.0 / snr**2\nmethod = \"dSPM\"  # use dSPM method (could also be MNE or sLORETA)\n\n# Load data\ninverse_operator = read_inverse_operator(fname_inv)\nlabel = mne.read_label(fname_label)\nraw = mne.io.read_raw_fif(fname_raw)\nevents = mne.read_events(fname_event)\n\n# Set up pick list\ninclude = []\nraw.info[\"bads\"] += [\"EEG 053\"]  # bads + 1 more\n\n# pick MEG channels\npicks = mne.pick_types(\n    raw.info, meg=True, eeg=False, stim=False, eog=True, include=include, exclude=\"bads\"\n)\n# Read epochs\nepochs = mne.Epochs(\n    raw,\n    events,\n    event_id,\n    tmin,\n    tmax,\n    picks=picks,\n    baseline=(None, 0),\n    reject=dict(mag=4e-12, grad=4000e-13, eog=150e-6),\n)\n\n# define frequencies of interest\nfmin, fmax = 0.0, 70.0\nbandwidth = 4.0  # bandwidth of the windows in Hz"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Compute source space PSD in label\n\n..note:: By using \"return_generator=True\" stcs will be a generator object\n         instead of a list. This allows us so to iterate without having to\n         keep everything in memory.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "n_epochs_use = 10\nstcs = compute_source_psd_epochs(\n    epochs[:n_epochs_use],\n    inverse_operator,\n    lambda2=lambda2,\n    method=method,\n    fmin=fmin,\n    fmax=fmax,\n    bandwidth=bandwidth,\n    label=label,\n    return_generator=True,\n    verbose=True,\n)\n\n# compute average PSD over the first 10 epochs\npsd_avg = 0.0\nfor i, stc in enumerate(stcs):\n    psd_avg += stc.data\npsd_avg /= n_epochs_use\nfreqs = stc.times  # the frequencies are stored here\nstc.data = psd_avg  # overwrite the last epoch's data with the average"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Visualize the 10 Hz PSD:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "brain = stc.plot(\n    initial_time=10.0,\n    hemi=\"lh\",\n    views=\"lat\",  # 10 HZ\n    clim=dict(kind=\"value\", lims=(20, 40, 60)),\n    smoothing_steps=3,\n    subjects_dir=subjects_dir,\n)\nbrain.add_label(label, borders=True, color=\"k\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Visualize the entire spectrum:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "fig, ax = plt.subplots()\nax.plot(freqs, psd_avg.mean(axis=0))\nax.set_xlabel(\"Freq (Hz)\")\nax.set_xlim(stc.times[[0, -1]])\nax.set_ylabel(\"Power Spectral Density\")"
      ]
    }
  ],
  "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
}