Update documentation

.buildinfo

@@ -0,0 +1,4 @@
+# Sphinx build info version 1
+# This file hashes the configuration used when building these files. When it is not found, a full rebuild will be done.
+config: b086e516bdfb25948dbfc05f6c2139c4
+tags: 645f666f9bcd5a90fca523b33c5a78b7

_downloads/127999bfe25ff667ca6d2d8ac8525a33/example_dss_line.ipynb

@@ -0,0 +1,140 @@
+{
+  "cells": [
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "\n# Remove line noise with ZapLine\n\nFind a spatial filter to get rid of line noise [1]_.\n\nUses meegkit.dss_line().\n\n## References\n.. [1] de Cheveign\u00e9, A. (2019). ZapLine: A simple and effective method to\n    remove power line artifacts [Preprint]. https://doi.org/10.1101/782029\n"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": null,
+      "metadata": {
+        "collapsed": false
+      },
+      "outputs": [],
+      "source": [
+        "# Authors: Maciej Szul <maciej.szul@isc.cnrs.fr>\n#          Nicolas Barascud <nicolas.barascud@gmail.com>\nimport os\n\nimport matplotlib.pyplot as plt\nimport numpy as np\nfrom scipy import signal\n\nfrom meegkit import dss\nfrom meegkit.utils import create_line_data, unfold"
+      ]
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "# Line noise removal\n\n"
+      ]
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "## Remove line noise with dss_line()\nWe first generate some noisy data to work with\n\n"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": null,
+      "metadata": {
+        "collapsed": false
+      },
+      "outputs": [],
+      "source": [
+        "sfreq = 250\nfline = 50\nnsamples = 10000\nnchans = 10\ndata = create_line_data(n_samples=3 * nsamples, n_chans=nchans,\n                        n_trials=1, fline=fline / sfreq, SNR=2)[0]\ndata = data[..., 0]  # only take first trial\n\n# Apply dss_line (ZapLine)\nout, _ = dss.dss_line(data, fline, sfreq, nkeep=1)"
+      ]
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "Plot before/after\n\n"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": null,
+      "metadata": {
+        "collapsed": false
+      },
+      "outputs": [],
+      "source": [
+        "f, ax = plt.subplots(1, 2, sharey=True)\nf, Pxx = signal.welch(data, sfreq, nperseg=500, axis=0, return_onesided=True)\nax[0].semilogy(f, Pxx)\nf, Pxx = signal.welch(out, sfreq, nperseg=500, axis=0, return_onesided=True)\nax[1].semilogy(f, Pxx)\nax[0].set_xlabel(\"frequency [Hz]\")\nax[1].set_xlabel(\"frequency [Hz]\")\nax[0].set_ylabel(\"PSD [V**2/Hz]\")\nax[0].set_title(\"before\")\nax[1].set_title(\"after\")\nplt.show()"
+      ]
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "## Remove line noise with dss_line_iter()\nWe first load some noisy data to work with\n\n"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": null,
+      "metadata": {
+        "collapsed": false
+      },
+      "outputs": [],
+      "source": [
+        "data = np.load(os.path.join(\"..\", \"tests\", \"data\", \"dss_line_data.npy\"))\nfline = 50\nsfreq = 200\nprint(data.shape)  # n_samples, n_chans, n_trials\n\n# Apply dss_line(), removing only one component\nout1, _ = dss.dss_line(data, fline, sfreq, nremove=1, nfft=400)"
+      ]
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "Now try dss_line_iter(). This applies dss_line() repeatedly until the\nartifact is gone\n\n"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": null,
+      "metadata": {
+        "collapsed": false
+      },
+      "outputs": [],
+      "source": [
+        "out2, iterations = dss.dss_line_iter(data, fline, sfreq, nfft=400)\nprint(f\"Removed {iterations} components\")"
+      ]
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "Plot results with dss_line() vs. dss_line_iter()\n\n"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": null,
+      "metadata": {
+        "collapsed": false
+      },
+      "outputs": [],
+      "source": [
+        "f, ax = plt.subplots(1, 2, sharey=True)\nf, Pxx = signal.welch(unfold(out1), sfreq, nperseg=200, axis=0,\n                      return_onesided=True)\nax[0].semilogy(f, Pxx, lw=.5)\nf, Pxx = signal.welch(unfold(out2), sfreq, nperseg=200, axis=0,\n                      return_onesided=True)\nax[1].semilogy(f, Pxx, lw=.5)\nax[0].set_xlabel(\"frequency [Hz]\")\nax[1].set_xlabel(\"frequency [Hz]\")\nax[0].set_ylabel(\"PSD [V**2/Hz]\")\nax[0].set_title(\"dss_line\")\nax[1].set_title(\"dss_line_iter\")\nplt.tight_layout()\nplt.show()"
+      ]
+    }
+  ],
+  "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.10.13"
+    }
+  },
+  "nbformat": 4,
+  "nbformat_minor": 0
+}

_downloads/1d7c08ebfaa26dd44535867493fad39a/example_star_dss.py

@@ -0,0 +1,111 @@
+"""
+Example demonstrating STAR + DSS
+================================
+
+This example shows how one can effectively combine STAR and DSS to recover
+signal components which would not have been discoverable with either these
+two techniques alone, due to the presence of strong artifacts.
+
+This example replicates figure 1 in [1]_.
+
+References
+----------
+.. [1] de Cheveigné A (2016) Sparse Time Artifact Removal, Journal of
+   Neuroscience Methods, 262, 14-20, doi:10.1016/j.jneumeth.2016.01.005
+
+"""
+import matplotlib.pyplot as plt
+import numpy as np
+from scipy.optimize import leastsq
+
+from meegkit import dss, star
+from meegkit.utils import demean, normcol, tscov
+
+# import config  # noqa
+
+rng = np.random.default_rng(9)
+
+###############################################################################
+# Create simulated data
+# -----------------------------------------------------------------------------
+# Simulated data consist of N channels, 1 sinusoidal target, N-3 noise sources,
+# with temporally local artifacts on each channel.
+
+# Source
+n_chans, n_samples = 10, 1000
+f = 2
+target = np.sin(np.arange(n_samples) / n_samples * 2 * np.pi * f)
+target = target[:, np.newaxis]
+noise = rng.standard_normal((n_samples, n_chans - 3))
+
+# Create artifact signal
+SNR = np.sqrt(1)
+x0 = normcol(np.dot(noise, rng.standard_normal((noise.shape[1], n_chans)))) + \
+    SNR * target * rng.standard_normal((1, n_chans))
+x0 = demean(x0)
+artifact = np.zeros(x0.shape)
+for k in np.arange(n_chans):
+    artifact[k * 100 + np.arange(20), k] = 1
+x = x0 + 10 * artifact
+
+
+def _sine_fit(x):
+    """Fit a sinusoidal trend."""
+    guess_mean = np.mean(x)
+    guess_std = np.std(x)
+    guess_phase = 0
+    t = np.linspace(0, 4 * np.pi, x.shape[0])
+
+    # Optimization function, in this case, we want to minimize the difference
+    # between the actual data and our "guessed" parameters
+    def func(y):
+        return np.mean(x - (y[0] * np.sin(t + y[1]) + y[2])[:, None], 1)
+
+    est_std, est_phase, est_mean = leastsq(
+        func, [guess_std, guess_phase, guess_mean])[0]
+    data_fit = est_std * np.sin(t + est_phase) + est_mean
+    return np.tile(data_fit, (x.shape[1], 1)).T
+
+
+###############################################################################
+# 1) Apply STAR
+# -----------------------------------------------------------------------------
+y, w, _ = star.star(x, 2)
+
+###############################################################################
+# 2) Apply DSS on raw data
+# -----------------------------------------------------------------------------
+c0, _ = tscov(x)
+c1, _ = tscov(x - _sine_fit(x))
+[todss, _, pwr0, pwr1] = dss.dss0(c0, c1)
+z1 = normcol(np.dot(x, todss))
+
+###############################################################################
+# 3) Apply DSS on STAR-ed data
+# -----------------------------------------------------------------------------
+# Here the bias function is the original signal minus the sinusoidal trend.
+c0, _ = tscov(y)
+c1, _ = tscov(y - _sine_fit(y))
+[todss, _, pwr0, pwr1] = dss.dss0(c0, c1)
+z2 = normcol(np.dot(y, todss))
+
+###############################################################################
+# Plots
+# -----------------------------------------------------------------------------
+f, (ax0, ax1, ax2, ax3) = plt.subplots(4, 1, figsize=(7, 9))
+ax0.plot(target, lw=.5)
+ax0.set_title("Target")
+
+ax1.plot(x, lw=.5)
+ax1.set_title(f"Signal + Artifacts (SNR = {SNR})")
+
+ax2.plot(z1[:, 0], lw=.5, label="Best DSS component")
+ax2.set_title("DSS")
+ax2.legend(loc="lower right")
+
+ax3.plot(z2[:, 0], lw=.5, label="Best DSS component")
+ax3.set_title("STAR + DSS")
+ax3.legend(loc="lower right")
+
+f.set_tight_layout(True)
+plt.show()

_downloads/22a83c59bf23910872ecdf61ccc22837/example_dss.ipynb

@@ -0,0 +1,97 @@
+{
+  "cells": [
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "\n# DSS example\n\nFind the linear combinations of multichannel data that maximize repeatability\nover trials.\n\nUses meegkit.dss0().\n"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": null,
+      "metadata": {
+        "collapsed": false
+      },
+      "outputs": [],
+      "source": [
+        "import matplotlib.pyplot as plt\nimport numpy as np\n\nfrom meegkit import dss\nfrom meegkit.utils import fold, rms, tscov, unfold\n\nrng = np.random.default_rng(5)"
+      ]
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "## Create simulated data\n\n"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": null,
+      "metadata": {
+        "collapsed": false
+      },
+      "outputs": [],
+      "source": [
+        "# Data are time * channel * trials.\nn_samples = 100 * 3\nn_chans = 30\nn_trials = 100\nnoise_dim = 20  # dimensionality of noise\n\n# Source signal\nsource = np.hstack((\n    np.zeros((n_samples // 3,)),\n    np.sin(2 * np.pi * np.arange(n_samples // 3) / (n_samples / 3)).T,\n    np.zeros((n_samples // 3,))))[np.newaxis].T\ns = source * rng.standard_normal((1, n_chans))  # 300 * 30\ns = s[:, :, np.newaxis]\ns = np.tile(s, (1, 1, 100))\n\n# Noise\nnoise = np.dot(\n    unfold(rng.standard_normal((n_samples, noise_dim, n_trials))),\n    rng.standard_normal((noise_dim, n_chans)))\nnoise = fold(noise, n_samples)\n\n# Mix signal and noise\nSNR = 0.1\ndata = noise / rms(noise.flatten()) + SNR * s / rms(s.flatten())"
+      ]
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "## Apply DSS to clean them\n\n"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": null,
+      "metadata": {
+        "collapsed": false
+      },
+      "outputs": [],
+      "source": [
+        "# Compute original and biased covariance matrices\nc0, _ = tscov(data)\n\n# In this case the biased covariance is simply the covariance of the mean over\n# trials\nc1, _ = tscov(np.mean(data, 2))\n\n# Apply DSS\n[todss, _, pwr0, pwr1] = dss.dss0(c0, c1)\nz = fold(np.dot(unfold(data), todss), epoch_size=n_samples)\n\n# Find best components\nbest_comp = np.mean(z[:, 0, :], -1)"
+      ]
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "## Plot results\n\n"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": null,
+      "metadata": {
+        "collapsed": false
+      },
+      "outputs": [],
+      "source": [
+        "f, (ax1, ax2, ax3) = plt.subplots(3, 1)\nax1.plot(source, label=\"source\")\nax2.plot(np.mean(data, 2), label=\"data\")\nax3.plot(best_comp, label=\"recovered\")\nplt.legend()\nplt.show()"
+      ]
+    }
+  ],
+  "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.10.13"
+    }
+  },
+  "nbformat": 4,
+  "nbformat_minor": 0
+}