{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n\n# Frequency and time-frequency sensor analysis\n\nThe objective is to show you how to explore the spectral content\nof your data (frequency and time-frequency). Here we'll work on Epochs.\n\nWe will use this dataset: `somato-dataset`. It contains so-called event\nrelated synchronizations (ERS) / desynchronizations (ERD) in the beta band.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Authors: Alexandre Gramfort <alexandre.gramfort@inria.fr>\n#          Stefan Appelhoff <stefan.appelhoff@mailbox.org>\n#          Richard H\u00f6chenberger <richard.hoechenberger@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 matplotlib.pyplot as plt\nimport numpy as np\n\nimport mne\nfrom mne.datasets import somato"
      ]
    },
    {
      "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# Setup for reading the raw data\nraw = mne.io.read_raw_fif(raw_fname)\n# crop and resample just to reduce computation time\nraw.crop(120, 360).load_data().resample(200)\nevents = mne.find_events(raw, stim_channel=\"STI 014\")\n\n# picks MEG gradiometers\npicks = mne.pick_types(raw.info, meg=\"grad\", eeg=False, eog=True, stim=False)\n\n# Construct Epochs\nevent_id, tmin, tmax = 1, -1.0, 3.0\nbaseline = (None, 0)\nepochs = mne.Epochs(\n    raw,\n    events,\n    event_id,\n    tmin,\n    tmax,\n    picks=picks,\n    baseline=baseline,\n    reject=dict(grad=4000e-13, eog=350e-6),\n    preload=True,\n)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Frequency analysis\n\nWe start by exploring the frequency content of our epochs.\n\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Let's first check out all channel types by averaging across epochs.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "epochs.compute_psd(fmin=2.0, fmax=40.0).plot(\n    average=True, amplitude=False, picks=\"data\", exclude=\"bads\"\n)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Now, let's take a look at the spatial distributions of the PSD, averaged\nacross epochs and frequency bands.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "epochs.compute_psd().plot_topomap(ch_type=\"grad\", normalize=False, contours=0)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Alternatively, you can also create PSDs from `~mne.Epochs` methods directly.\n\n<div class=\"alert alert-info\"><h4>Note</h4><p>In contrast to the methods for visualization, the ``compute_psd`` methods\n   do **not** scale the data from SI units to more \"convenient\" values. So\n   when e.g. calculating the PSD of gradiometers via\n   :meth:`~mne.Epochs.compute_psd`, you will get the power as\n   ``(T/m)\u00b2/Hz`` (instead of ``(fT/cm)\u00b2/Hz`` via\n   :meth:`~mne.Epochs.plot_psd`).</p></div>\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "_, ax = plt.subplots()\nspectrum = epochs.compute_psd(fmin=2.0, fmax=40.0, tmax=3.0, n_jobs=None)\n# average across epochs first\nmean_spectrum = spectrum.average()\npsds, freqs = mean_spectrum.get_data(return_freqs=True)\n# then convert to dB and take mean & standard deviation across channels\npsds = 10 * np.log10(psds)\npsds_mean = psds.mean(axis=0)\npsds_std = psds.std(axis=0)\n\nax.plot(freqs, psds_mean, color=\"k\")\nax.fill_between(\n    freqs,\n    psds_mean - psds_std,\n    psds_mean + psds_std,\n    color=\"k\",\n    alpha=0.5,\n    edgecolor=\"none\",\n)\nax.set(\n    title=\"Multitaper PSD (gradiometers)\",\n    xlabel=\"Frequency (Hz)\",\n    ylabel=\"Power Spectral Density (dB)\",\n)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Notably, :meth:`mne.Epochs.compute_psd` supports the keyword argument\n``average``, which specifies how to estimate the PSD based on the individual\nwindowed segments. The default is ``average='mean'``, which simply calculates\nthe arithmetic mean across segments. Specifying ``average='median'``, in\ncontrast, returns the PSD based on the median of the segments (corrected for\nbias relative to the mean), which is a more robust measure.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Estimate PSDs based on \"mean\" and \"median\" averaging for comparison.\nkwargs = dict(fmin=2, fmax=40, n_jobs=None)\npsds_welch_mean, freqs_mean = epochs.compute_psd(\n    \"welch\", average=\"mean\", **kwargs\n).get_data(return_freqs=True)\npsds_welch_median, freqs_median = epochs.compute_psd(\n    \"welch\", average=\"median\", **kwargs\n).get_data(return_freqs=True)\n\n# Convert power to dB scale.\npsds_welch_mean = 10 * np.log10(psds_welch_mean)\npsds_welch_median = 10 * np.log10(psds_welch_median)\n\n# We will only plot the PSD for a single sensor in the first epoch.\nch_name = \"MEG 0122\"\nch_idx = epochs.info[\"ch_names\"].index(ch_name)\nepo_idx = 0\n\n_, ax = plt.subplots()\nax.plot(\n    freqs_mean,\n    psds_welch_mean[epo_idx, ch_idx, :],\n    color=\"k\",\n    ls=\"-\",\n    label=\"mean of segments\",\n)\nax.plot(\n    freqs_median,\n    psds_welch_median[epo_idx, ch_idx, :],\n    color=\"k\",\n    ls=\"--\",\n    label=\"median of segments\",\n)\n\nax.set(\n    title=f\"Welch PSD ({ch_name}, Epoch {epo_idx})\",\n    xlabel=\"Frequency (Hz)\",\n    ylabel=\"Power Spectral Density (dB)\",\n)\nax.legend(loc=\"upper right\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Lastly, we can also retrieve the unaggregated segments by passing\n``average=None`` to :meth:`mne.Epochs.compute_psd`. The dimensions of\nthe returned array are ``(n_epochs, n_sensors, n_freqs, n_segments)``.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "welch_unagg = epochs.compute_psd(\"welch\", average=None, **kwargs)\nprint(welch_unagg.shape)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n## Time-frequency analysis: power and inter-trial coherence\n\nWe now compute time-frequency representations (TFRs) from our Epochs.\nWe'll look at power and inter-trial coherence (ITC).\n\nTo this we'll use the function :func:`mne.time_frequency.tfr_morlet`\nbut you can also use :func:`mne.time_frequency.tfr_multitaper`\nor :func:`mne.time_frequency.tfr_stockwell`.\n\n<div class=\"alert alert-info\"><h4>Note</h4><p>The ``decim`` parameter reduces the sampling rate of the time-frequency\n      decomposition by the defined factor. This is usually done to reduce\n      memory usage. For more information refer to the documentation of\n      :func:`mne.time_frequency.tfr_morlet`.</p></div>\n\ndefine frequencies of interest (log-spaced)\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "freqs = np.logspace(*np.log10([6, 35]), num=8)\nn_cycles = freqs / 2.0  # different number of cycle per frequency\npower, itc = epochs.compute_tfr(\n    method=\"morlet\",\n    freqs=freqs,\n    n_cycles=n_cycles,\n    average=True,\n    return_itc=True,\n    decim=3,\n)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Inspect power\n\n<div class=\"alert alert-info\"><h4>Note</h4><p>The generated figures are interactive. In the topo you can click\n    on an image to visualize the data for one sensor.\n    You can also select a portion in the time-frequency plane to\n    obtain a topomap for a certain time-frequency region.</p></div>\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "power.plot_topo(baseline=(-0.5, 0), mode=\"logratio\", title=\"Average power\")\npower.plot(picks=[82], baseline=(-0.5, 0), mode=\"logratio\", title=power.ch_names[82])\n\nfig, axes = plt.subplots(1, 2, figsize=(7, 4), layout=\"constrained\")\ntopomap_kw = dict(\n    ch_type=\"grad\", tmin=0.5, tmax=1.5, baseline=(-0.5, 0), mode=\"logratio\", show=False\n)\nplot_dict = dict(Alpha=dict(fmin=8, fmax=12), Beta=dict(fmin=13, fmax=25))\nfor ax, (title, fmin_fmax) in zip(axes, plot_dict.items()):\n    power.plot_topomap(**fmin_fmax, axes=ax, **topomap_kw)\n    ax.set_title(title)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Joint Plot\nYou can also create a joint plot showing both the aggregated TFR\nacross channels and topomaps at specific times and frequencies to obtain\na quick overview regarding oscillatory effects across time and space.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "power.plot_joint(\n    baseline=(-0.5, 0), mode=\"mean\", tmin=-0.5, tmax=2, timefreqs=[(0.5, 10), (1.3, 8)]\n)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Inspect ITC\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "itc.plot_topo(title=\"Inter-Trial coherence\", vmin=0.0, vmax=1.0, cmap=\"Reds\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "<div class=\"alert alert-info\"><h4>Note</h4><p>Baseline correction can be applied to power or done in plots.\n    To illustrate the baseline correction in plots, the next line is\n    commented::\n\n    # power.apply_baseline(baseline=(-0.5, 0), mode='logratio')</p></div>\n\n## Exercise\n\n   - Visualize the inter-trial coherence values as topomaps as done with\n     power.\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
}