{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n\n# Plot sensor denoising using oversampled temporal projection\n\nThis demonstrates denoising using the OTP algorithm :footcite:`LarsonTaulu2018`\non data with with sensor artifacts (flux jumps) and random noise.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Author: Eric Larson <larson.eric.d@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 numpy as np\n\nimport mne\nfrom mne import find_events, fit_dipole\nfrom mne.datasets.brainstorm import bst_phantom_elekta\nfrom mne.io import read_raw_fif\n\nprint(__doc__)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Plot the phantom data, lowpassed to get rid of high-frequency artifacts.\nWe also crop to a single 10-second segment for speed.\nNotice that there are two large flux jumps on channel 1522 that could\nspread to other channels when performing subsequent spatial operations\n(e.g., Maxwell filtering, SSP, or ICA).\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "dipole_number = 1\ndata_path = bst_phantom_elekta.data_path()\nraw = read_raw_fif(data_path / \"kojak_all_200nAm_pp_no_chpi_no_ms_raw.fif\")\nraw.crop(40.0, 50.0).load_data()\norder = list(range(160, 170))\nraw.copy().filter(0.0, 40.0).plot(order=order, n_channels=10)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Now we can clean the data with OTP, lowpass, and plot. The flux jumps have\nbeen suppressed alongside the random sensor noise.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "raw_clean = mne.preprocessing.oversampled_temporal_projection(raw)\nraw_clean.filter(0.0, 40.0)\nraw_clean.plot(order=order, n_channels=10)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We can also look at the effect on single-trial phantom localization.\nSee the `tut-brainstorm-elekta-phantom`\nfor more information. Here we use a version that does single-trial\nlocalization across the 17 trials are in our 10-second window:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "def compute_bias(raw):\n    events = find_events(raw, \"STI201\", verbose=False)\n    events = events[1:]  # first one has an artifact\n    tmin, tmax = -0.2, 0.1\n    epochs = mne.Epochs(\n        raw,\n        events,\n        dipole_number,\n        tmin,\n        tmax,\n        baseline=(None, -0.01),\n        preload=True,\n        verbose=False,\n    )\n    sphere = mne.make_sphere_model(r0=(0.0, 0.0, 0.0), head_radius=None, verbose=False)\n    cov = mne.compute_covariance(epochs, tmax=0, method=\"oas\", rank=None, verbose=False)\n    idx = epochs.time_as_index(0.036)[0]\n    data = epochs.get_data(copy=False)[:, :, idx].T\n    evoked = mne.EvokedArray(data, epochs.info, tmin=0.0)\n    dip = fit_dipole(evoked, cov, sphere, n_jobs=None, verbose=False)[0]\n    actual_pos = mne.dipole.get_phantom_dipoles()[0][dipole_number - 1]\n    misses = 1000 * np.linalg.norm(dip.pos - actual_pos, axis=-1)\n    return misses\n\n\nbias = compute_bias(raw)\nprint(f\"Raw bias: {np.mean(bias):0.1f}mm (worst: {np.max(bias):0.1f}mm)\")\nbias_clean = compute_bias(raw_clean)\nprint(f\"OTP bias: {np.mean(bias_clean):0.1f}mm (worst: {np.max(bias_clean):0.1f}m)\")"
      ]
    },
    {
      "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
}