{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n\n# Cortical Signal Suppression (CSS) for removal of cortical signals\n\nThis script shows an example of how to use CSS\n:footcite:`Samuelsson2019` . CSS suppresses the cortical contribution\nto the signal subspace in EEG data using MEG data, facilitating\ndetection of subcortical signals. We will illustrate how it works by\nsimulating one cortical and one subcortical oscillation at different\nfrequencies; 40 Hz and 239 Hz for cortical and subcortical activity,\nrespectively, then process it with CSS and look at the power spectral\ndensity of the raw and processed data.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Author: John G Samuelsson <johnsam@mit.edu>\n# License: BSD-3-Clause\n# Copyright the MNE-Python contributors.\n\nimport matplotlib.pyplot as plt\nimport numpy as np\n\nimport mne\nfrom mne.datasets import sample\nfrom mne.simulation import simulate_evoked, simulate_sparse_stc"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Load sample subject data\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "data_path = sample.data_path()\nsubjects_dir = data_path / \"subjects\"\nmeg_path = data_path / \"MEG\" / \"sample\"\nfwd_fname = meg_path / \"sample_audvis-meg-eeg-oct-6-fwd.fif\"\nave_fname = meg_path / \"sample_audvis-no-filter-ave.fif\"\ncov_fname = meg_path / \"sample_audvis-cov.fif\"\ntrans_fname = meg_path / \"sample_audvis_raw-trans.fif\"\nbem_fname = subjects_dir / \"sample\" / \"bem\" / \"/sample-5120-bem-sol.fif\"\n\nraw = mne.io.read_raw_fif(meg_path / \"sample_audvis_raw.fif\")\nfwd = mne.read_forward_solution(fwd_fname)\nfwd = mne.convert_forward_solution(fwd, force_fixed=True, surf_ori=True)\nfwd = mne.pick_types_forward(fwd, meg=True, eeg=True, exclude=raw.info[\"bads\"])\ncov = mne.read_cov(cov_fname)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Find patches (labels) to activate\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "all_labels = mne.read_labels_from_annot(subject=\"sample\", subjects_dir=subjects_dir)\nlabels = []\nfor select_label in [\"parahippocampal-lh\", \"postcentral-rh\"]:\n    labels.append([lab for lab in all_labels if lab.name in select_label][0])\nhiplab, postcenlab = labels"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Simulate one cortical dipole (40 Hz) and one subcortical (239 Hz)\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "def cortical_waveform(times):\n    \"\"\"Create a cortical waveform.\"\"\"\n    return 10e-9 * np.cos(times * 2 * np.pi * 40)\n\n\ndef subcortical_waveform(times):\n    \"\"\"Create a subcortical waveform.\"\"\"\n    return 10e-9 * np.cos(times * 2 * np.pi * 239)\n\n\ntimes = np.linspace(0, 0.5, int(0.5 * raw.info[\"sfreq\"]))\nstc = simulate_sparse_stc(\n    fwd[\"src\"],\n    n_dipoles=2,\n    times=times,\n    location=\"center\",\n    subjects_dir=subjects_dir,\n    labels=[postcenlab, hiplab],\n    data_fun=cortical_waveform,\n)\nstc.data[np.where(np.isin(stc.vertices[0], hiplab.vertices))[0], :] = (\n    subcortical_waveform(times)\n)\nevoked = simulate_evoked(fwd, stc, raw.info, cov, nave=15)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Process with CSS and plot PSD of EEG data before and after processing\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "evoked_subcortical = mne.preprocessing.cortical_signal_suppression(evoked, n_proj=6)\nchs = mne.pick_types(evoked.info, meg=False, eeg=True)\n\npsd = np.mean(np.abs(np.fft.rfft(evoked.data)) ** 2, axis=0)\npsd_proc = np.mean(np.abs(np.fft.rfft(evoked_subcortical.data)) ** 2, axis=0)\nfreq = np.arange(\n    0, stop=int(evoked.info[\"sfreq\"] / 2), step=evoked.info[\"sfreq\"] / (2 * len(psd))\n)\n\nfig, ax = plt.subplots()\nax.plot(freq, psd, label=\"raw\")\nax.plot(freq, psd_proc, label=\"processed\")\nax.text(0.2, 0.7, \"cortical\", transform=ax.transAxes)\nax.text(0.8, 0.25, \"subcortical\", transform=ax.transAxes)\nax.set(ylabel=\"EEG Power spectral density\", xlabel=\"Frequency (Hz)\")\nax.legend()\n\n# References\n# ^^^^^^^^^^\n#\n# .. footbibliography::"
      ]
    }
  ],
  "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
}