{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n\n# Compute a cross-spectral density (CSD) matrix\n\nA cross-spectral density (CSD) matrix is similar to a covariance matrix, but in\nthe time-frequency domain. It is the first step towards computing\nsensor-to-sensor coherence or a DICS beamformer.\n\nThis script demonstrates the three methods that MNE-Python provides to compute\nthe CSD:\n\n1. Using short-term Fourier transform: :func:`mne.time_frequency.csd_fourier`\n2. Using a multitaper approach: :func:`mne.time_frequency.csd_multitaper`\n3. Using Morlet wavelets: :func:`mne.time_frequency.csd_morlet`\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Author: Marijn van Vliet <w.m.vanvliet@gmail.com>\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.datasets import sample\nfrom mne.time_frequency import csd_fourier, csd_morlet, csd_multitaper\n\nprint(__doc__)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "In the following example, the computation of the CSD matrices can be\nperformed using multiple cores. Set ``n_jobs`` to a value >1 to select the\nnumber of cores to use.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "n_jobs = 1"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Loading the sample dataset.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "data_path = sample.data_path()\nmeg_path = data_path / \"MEG\" / \"sample\"\nfname_raw = meg_path / \"sample_audvis_raw.fif\"\nfname_event = meg_path / \"sample_audvis_raw-eve.fif\"\nraw = mne.io.read_raw_fif(fname_raw)\nevents = mne.read_events(fname_event)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "By default, CSD matrices are computed using all MEG/EEG channels. When\ninterpreting a CSD matrix with mixed sensor types, be aware that the\nmeasurement units, and thus the scalings, differ across sensors. In this\nexample, for speed and clarity, we select a single channel type:\ngradiometers.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "picks = mne.pick_types(raw.info, meg=\"grad\")\n\n# Make some epochs, based on events with trigger code 1\nepochs = mne.Epochs(\n    raw,\n    events,\n    event_id=1,\n    tmin=-0.2,\n    tmax=1,\n    picks=picks,\n    baseline=(None, 0),\n    reject=dict(grad=4000e-13),\n    preload=True,\n)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Computing CSD matrices using short-term Fourier transform and (adaptive)\nmultitapers is straightforward:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "csd_fft = csd_fourier(epochs, fmin=15, fmax=20, n_jobs=n_jobs)\ncsd_mt = csd_multitaper(epochs, fmin=15, fmax=20, adaptive=True, n_jobs=n_jobs)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "When computing the CSD with Morlet wavelets, you specify the exact\nfrequencies at which to compute it. For each frequency, a corresponding\nwavelet will be constructed and convolved with the signal, resulting in a\ntime-frequency decomposition.\n\nThe CSD is constructed by computing the correlation between the\ntime-frequency representations between all sensor-to-sensor pairs. The\ntime-frequency decomposition originally has the same sampling rate as the\nsignal, in our case ~600Hz. This means the decomposition is over-specified in\ntime and we may not need to use all samples during our CSD computation, just\nenough to get a reliable correlation statistic. By specifying ``decim=10``,\nwe use every 10th sample, which will greatly speed up the computation and\nwill have a minimal effect on the CSD.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "frequencies = [16, 17, 18, 19, 20]\ncsd_wav = csd_morlet(epochs, frequencies, decim=10, n_jobs=n_jobs)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "The resulting :class:`mne.time_frequency.CrossSpectralDensity` objects have a\nplotting function we can use to compare the results of the different methods.\nWe're plotting the mean CSD across frequencies.\n:func:`mne.time_frequency.CrossSpectralDensity.plot()` returns a list of\ncreated figures; in this case, each returned list has only one figure\nso we use a Python trick of including a comma after our variable name\nto assign the figure (not the list) to our ``fig`` variable:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "plot_dict = {\n    \"Short-time Fourier transform\": csd_fft,\n    \"Adaptive multitapers\": csd_mt,\n    \"Morlet wavelet transform\": csd_wav,\n}\nfor title, csd in plot_dict.items():\n    (fig,) = csd.mean().plot()\n    fig.suptitle(title)"
      ]
    }
  ],
  "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
}