{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n\n# Transform EEG data using current source density (CSD)\n\nThis script shows an example of how to use CSD\n:footcite:`PerrinEtAl1987,PerrinEtAl1989,Cohen2014,KayserTenke2015`.\nCSD takes the spatial Laplacian of the sensor signal (derivative in both\nx and y). It does what a planar gradiometer does in MEG. Computing these\nspatial derivatives reduces point spread. CSD transformed data have a sharper\nor more distinct topography, reducing the negative impact of volume conduction.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Authors: Alex Rockhill <aprockhill@mailbox.org>\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 sample\n\nprint(__doc__)\n\ndata_path = sample.data_path()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Load sample subject data\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "meg_path = data_path / \"MEG\" / \"sample\"\nraw = mne.io.read_raw_fif(meg_path / \"sample_audvis_raw.fif\")\nraw = raw.pick(picks=[\"eeg\", \"eog\", \"ecg\", \"stim\"], exclude=\"bads\").load_data()\nevents = mne.find_events(raw)\nraw.set_eeg_reference(projection=True).apply_proj()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Plot the raw data and CSD-transformed raw data:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "raw_csd = mne.preprocessing.compute_current_source_density(raw)\nraw.plot()\nraw_csd.plot()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Also look at the power spectral densities:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "raw.compute_psd().plot(picks=\"data\", exclude=\"bads\", amplitude=False)\nraw_csd.compute_psd().plot(picks=\"data\", exclude=\"bads\", amplitude=False)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "CSD can also be computed on Evoked (averaged) data.\nHere we epoch and average the data so we can demonstrate that.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "event_id = {\n    \"auditory/left\": 1,\n    \"auditory/right\": 2,\n    \"visual/left\": 3,\n    \"visual/right\": 4,\n    \"smiley\": 5,\n    \"button\": 32,\n}\nepochs = mne.Epochs(raw, events, event_id=event_id, tmin=-0.2, tmax=0.5, preload=True)\nevoked = epochs[\"auditory\"].average()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "First let's look at how CSD affects scalp topography:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "times = np.array([-0.1, 0.0, 0.05, 0.1, 0.15])\nevoked_csd = mne.preprocessing.compute_current_source_density(evoked)\nevoked.plot_joint(title=\"Average Reference\", show=False)\nevoked_csd.plot_joint(title=\"Current Source Density\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "CSD has parameters ``stiffness`` and ``lambda2`` affecting smoothing and\nspline flexibility, respectively. Let's see how they affect the solution:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "fig, ax = plt.subplots(4, 4, layout=\"constrained\")\nfig.set_size_inches(10, 10)\nfor i, lambda2 in enumerate([0, 1e-7, 1e-5, 1e-3]):\n    for j, m in enumerate([5, 4, 3, 2]):\n        this_evoked_csd = mne.preprocessing.compute_current_source_density(\n            evoked, stiffness=m, lambda2=lambda2\n        )\n        this_evoked_csd.plot_topomap(\n            0.1, axes=ax[i, j], contours=4, time_unit=\"s\", colorbar=False, show=False\n        )\n        ax[i, j].set_title(f\"stiffness={m}\\n\u03bb\u00b2={lambda2}\")"
      ]
    },
    {
      "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
}