[ENH] Add (poly)coherence functions (#76)

nbara · Jun 24, 2024 · cb2e817 · cb2e817
1 parent e220ac7
commit cb2e817
Show file tree

Hide file tree

Showing 29 changed files with 1,236 additions and 88 deletions.
diff --git a/README.md b/README.md
@@ -154,11 +154,19 @@ and was adapted to python by [Giuseppe Ferraro](mailto:giuseppe.ferraro@isae-sup
     2022 Sep 27;22(19):7314. https://doi.org/10.3390/s22197314.
 ```
 
-### 6. Real-Time Phase Estimation
+### 6. Phase Estimation
 
-This code is based on the Matlab implementation from [Michael Rosenblum](http://www.stat.physik.uni-potsdam.de/~mros), and its corresponding paper [1].
+The oscillator code is based on the Matlab implementation from [Michael
+Rosenblum](http://www.stat.physik.uni-potsdam.de/~mros), and its corresponding
+paper [1]. The Endpoint Corrected Hilbert Transform (ECHT) method was adapted
+from [2].
 
 ```sql
-[1] Rosenblum, M., Pikovsky, A., Kühn, A.A. et al. Real-time estimation of phase and amplitude with application to neural data. Sci Rep 11, 18037 (2021). https://doi.org/10.1038/s41598-021-97560-5
+[1] Rosenblum, M., Pikovsky, A., Kühn, A.A. et al. Real-time estimation of phase 
+    and amplitude with application to neural data. Sci Rep 11, 18037 (2021). 
+    https://doi.org/10.1038/s41598-021-97560-5
+[2] Schreglmann, S. R., Wang, D., Peach, R. L., Li, J., Zhang, X., Latorre, A., 
+    ... & Grossman, N. (2021). Non-invasive suppression of essential tremor via
+    phase-locked disruption of its temporal coherence. Nature communications, 12(1), 363.
 
 ```
diff --git a/doc/conf.py b/doc/conf.py
@@ -26,7 +26,7 @@
 # -- Project information -----------------------------------------------------
 
 project = "MEEGkit"
-copyright = "2023, Nicolas Barascud"
+copyright = "2024, Nicolas Barascud"
 author = "Nicolas Barascud"
 release = meegkit.__version__
 version = meegkit.__version__
@@ -63,7 +63,7 @@
     "show-inheritance": True,
     "exclude-members": "__weakref__"
 }
-numpydoc_show_class_members = True
+numpydoc_show_class_members = False
 
 # The suffix(es) of source filenames.
 # You can specify multiple suffix as a list of string:
@@ -129,3 +129,5 @@
     "ignore_pattern": "config.py",
     "run_stale_examples": False,
 }
+
+suppress_warnings = ["config.cache"]
diff --git a/doc/index.rst b/doc/index.rst
@@ -46,7 +46,6 @@ Here is a list of the methods and techniques available in ``meegkit``:
    ~meegkit.tspca
    ~meegkit.utils
 
-
 Examples gallery
 ----------------
 

diff --git a/doc/modules/meegkit.phase.rst b/doc/modules/meegkit.phase.rst
@@ -0,0 +1,23 @@
+meegkit.phase
+=============
+
+.. automodule:: meegkit.phase
+
+    .. rubric:: Classes
+
+    .. autosummary::
+
+       NonResOscillator
+       ResOscillator
+       Device
+       ECHT
+
+    .. rubric:: Functions
+
+    .. autosummary::
+
+       locking_based_phase
+       rk
+       init_coefs
+       one_step_oscillator
+       one_step_integrator
diff --git a/doc/modules/meegkit.utils.rst b/doc/modules/meegkit.utils.rst
@@ -6,11 +6,14 @@ meegkit.utils
    .. autosummary::
 
     auditory
+    buffer
+    coherence
     covariances
     denoise
     matrix
     sig
     stats
+    trca
 
 |
 
@@ -23,6 +26,28 @@ Auditory
    .. autosummary::
 
 
+|
+
+----
+
+Buffer
+------
+.. automodule:: meegkit.utils.buffer
+
+   .. autosummary::
+
+
+|
+
+----
+
+Coherence
+---------
+.. automodule:: meegkit.utils.coherence
+
+   .. autosummary::
+
+
 |
 
 ----

diff --git a/doc/sg_execution_times.rst b/doc/sg_execution_times.rst
@@ -0,0 +1,73 @@
+
+:orphan:
+
+.. _sphx_glr_sg_execution_times:
+
+
+Computation times
+=================
+**00:29.509** total execution time for 13 files **from all galleries**:
+
+.. container::
+
+  .. raw:: html
+
+    <style scoped>
+    <link href="https://cdnjs.cloudflare.com/ajax/libs/twitter-bootstrap/5.3.0/css/bootstrap.min.css" rel="stylesheet" />
+    <link href="https://cdn.datatables.net/1.13.6/css/dataTables.bootstrap5.min.css" rel="stylesheet" />
+    </style>
+    <script src="https://code.jquery.com/jquery-3.7.0.js"></script>
+    <script src="https://cdn.datatables.net/1.13.6/js/jquery.dataTables.min.js"></script>
+    <script src="https://cdn.datatables.net/1.13.6/js/dataTables.bootstrap5.min.js"></script>
+    <script type="text/javascript" class="init">
+    $(document).ready( function () {
+        $('table.sg-datatable').DataTable({order: [[1, 'desc']]});
+    } );
+    </script>
+
+  .. list-table::
+   :header-rows: 1
+   :class: table table-striped sg-datatable
+
+   * - Example
+     - Time
+     - Mem (MB)
+   * - :ref:`sphx_glr_auto_examples_example_trca.py` (``../examples/example_trca.py``)
+     - 00:11.533
+     - 0.0
+   * - :ref:`sphx_glr_auto_examples_example_dss_line.py` (``../examples/example_dss_line.py``)
+     - 00:07.115
+     - 0.0
+   * - :ref:`sphx_glr_auto_examples_example_phase_estimation.py` (``../examples/example_phase_estimation.py``)
+     - 00:07.064
+     - 0.0
+   * - :ref:`sphx_glr_auto_examples_example_mcca.py` (``../examples/example_mcca.py``)
+     - 00:01.182
+     - 0.0
+   * - :ref:`sphx_glr_auto_examples_example_mcca_2.py` (``../examples/example_mcca_2.py``)
+     - 00:00.535
+     - 0.0
+   * - :ref:`sphx_glr_auto_examples_example_asr.py` (``../examples/example_asr.py``)
+     - 00:00.512
+     - 0.0
+   * - :ref:`sphx_glr_auto_examples_example_ress.py` (``../examples/example_ress.py``)
+     - 00:00.459
+     - 0.0
+   * - :ref:`sphx_glr_auto_examples_example_detrend.py` (``../examples/example_detrend.py``)
+     - 00:00.259
+     - 0.0
+   * - :ref:`sphx_glr_auto_examples_example_star_dss.py` (``../examples/example_star_dss.py``)
+     - 00:00.246
+     - 0.0
+   * - :ref:`sphx_glr_auto_examples_example_star.py` (``../examples/example_star.py``)
+     - 00:00.197
+     - 0.0
+   * - :ref:`sphx_glr_auto_examples_example_dss.py` (``../examples/example_dss.py``)
+     - 00:00.146
+     - 0.0
+   * - :ref:`sphx_glr_auto_examples_example_echt.py` (``../examples/example_echt.py``)
+     - 00:00.132
+     - 0.0
+   * - :ref:`sphx_glr_auto_examples_example_dering.py` (``../examples/example_dering.py``)
+     - 00:00.129
+     - 0.0
diff --git a/examples/example_echt.ipynb b/examples/example_echt.ipynb
@@ -0,0 +1,97 @@
+{
+  "cells": [
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "\n# Endpoint-corrected Hilbert transform (ECHT) phase estimation\n\nThis example shows how to causally estimate the phase of a signal using\n\nUses `meegkit.phase.ECHT()`.\n\n## References\n.. [1] Schreglmann, S. R., Wang, D., Peach, R. L., Li, J., Zhang, X., Latorre,\n    A., ... & Grossman, N. (2021). Non-invasive suppression of essential tremor\n    via phase-locked disruption of its temporal coherence. Nature\n    communications, 12(1), 363.\n"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": null,
+      "metadata": {
+        "collapsed": false
+      },
+      "outputs": [],
+      "source": [
+        "import matplotlib.pyplot as plt\nimport numpy as np\nfrom scipy.signal import hilbert\n\nfrom meegkit.phase import ECHT\n\nrng = np.random.default_rng(38872)\n\nplt.rcParams[\"axes.grid\"] = True\nplt.rcParams[\"grid.linestyle\"] = \":\""
+      ]
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "## Build data\nFirst, we generate a multi-component signal with amplitude and phase\nmodulations, as described in the paper [1]_.\n\n"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": null,
+      "metadata": {
+        "collapsed": false
+      },
+      "outputs": [],
+      "source": [
+        "f0 = 2\n\nN = 500\nsfreq = 100\ntime = np.linspace(0, N / sfreq, N)\nX = np.cos(2 * np.pi * f0 * time - np.pi / 4)\nphase_true = np.angle(hilbert(X))\nX += rng.normal(0, 0.5, N)  # Add noise"
+      ]
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "### Compute phase and amplitude\nWe compute the Hilbert phase, as well as the phase obtained with the ECHT\nfilter.\n\n"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": null,
+      "metadata": {
+        "collapsed": false
+      },
+      "outputs": [],
+      "source": [
+        "phase_hilbert = np.angle(hilbert(X))  # Hilbert phase\n\n# Compute ECHT-filtered signal\nfilt_BW = f0 / 2\nl_freq = f0 - filt_BW / 2\nh_freq = f0 + filt_BW / 2\necht = ECHT(l_freq, h_freq, sfreq)\n\nXf = echt.fit_transform(X)\nphase_echt = np.angle(Xf)"
+      ]
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "### Visualize signal\nPlot the results\n\n"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": null,
+      "metadata": {
+        "collapsed": false
+      },
+      "outputs": [],
+      "source": [
+        "fig, ax = plt.subplots(3, 1, figsize=(8, 6))\nax[0].plot(time, X)\nax[0].set_xlabel(\"Time (s)\")\nax[0].set_title(\"Test signal\")\nax[0].set_ylabel(\"Amplitude\")\nax[1].psd(X, Fs=sfreq, NFFT=2048*4, noverlap=sfreq)\nax[1].set_ylabel(\"PSD (dB/Hz)\")\nax[1].set_title(\"Test signal's Fourier spectrum\")\nax[2].plot(time, phase_true, label=\"True phase\", ls=\":\")\nax[2].plot(time, phase_echt, label=\"ECHT phase\", lw=.5, alpha=.8)\nax[2].plot(time, phase_hilbert, label=\"Hilbert phase\", lw=.5, alpha=.8)\nax[2].set_title(\"Phase\")\nax[2].set_ylabel(\"Amplitude\")\nax[2].set_xlabel(\"Time (s)\")\nax[2].legend(loc=\"upper right\", fontsize=\"small\")\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.12.2"
+    }
+  },
+  "nbformat": 4,
+  "nbformat_minor": 0
+}
diff --git a/examples/example_echt.py b/examples/example_echt.py
@@ -0,0 +1,81 @@
+"""
+Endpoint-corrected Hilbert transform (ECHT) phase estimation
+============================================================
+
+This example shows how to causally estimate the phase of a signal using the
+Endpoint-corrected Hilbert transform (ECHT) [1]_.
+
+Uses `meegkit.phase.ECHT()`.
+
+References
+----------
+.. [1] Schreglmann, S. R., Wang, D., Peach, R. L., Li, J., Zhang, X., Latorre,
+    A., ... & Grossman, N. (2021). Non-invasive suppression of essential tremor
+    via phase-locked disruption of its temporal coherence. Nature
+    communications, 12(1), 363.
+
+"""
+import matplotlib.pyplot as plt
+import numpy as np
+from scipy.signal import hilbert
+
+from meegkit.phase import ECHT
+
+rng = np.random.default_rng(38872)
+
+plt.rcParams["axes.grid"] = True
+plt.rcParams["grid.linestyle"] = ":"
+
+###############################################################################
+# Build data
+# -----------------------------------------------------------------------------
+# First, we generate a multi-component signal with amplitude and phase
+# modulations, as described in the paper [1]_.
+f0 = 2
+
+N = 500
+sfreq = 100
+time = np.linspace(0, N / sfreq, N)
+X = np.cos(2 * np.pi * f0 * time - np.pi / 4)
+phase_true = np.angle(hilbert(X))
+X += rng.normal(0, 0.5, N)  # Add noise
+
+###############################################################################
+# Compute phase and amplitude
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+# We compute the Hilbert phase, as well as the phase obtained with the ECHT
+# filter.
+phase_hilbert = np.angle(hilbert(X))  # Hilbert phase
+
+# Compute ECHT-filtered signal
+filt_BW = f0 / 2
+l_freq = f0 - filt_BW / 2
+h_freq = f0 + filt_BW / 2
+echt = ECHT(l_freq, h_freq, sfreq)
+
+Xf = echt.fit_transform(X)
+phase_echt = np.angle(Xf)
+
+###############################################################################
+# Visualize signal
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+# Here we plot the original signal, its Fourier spectrum, and the phase obtained
+# with the Hilbert transform and the ECHT filter. The ECHT filter provides a
+# much smoother phase estimate than the Hilbert transform
+fig, ax = plt.subplots(3, 1, figsize=(8, 6))
+ax[0].plot(time, X)
+ax[0].set_xlabel("Time (s)")
+ax[0].set_title("Test signal")
+ax[0].set_ylabel("Amplitude")
+ax[1].psd(X, Fs=sfreq, NFFT=2048*4, noverlap=sfreq)
+ax[1].set_ylabel("PSD (dB/Hz)")
+ax[1].set_title("Test signal's Fourier spectrum")
+ax[2].plot(time, phase_true, label="True phase", ls=":")
+ax[2].plot(time, phase_echt, label="ECHT phase", lw=.5, alpha=.8)
+ax[2].plot(time, phase_hilbert, label="Hilbert phase", lw=.5, alpha=.8)
+ax[2].set_title("Phase")
+ax[2].set_ylabel("Amplitude")
+ax[2].set_xlabel("Time (s)")
+ax[2].legend(loc="upper right", fontsize="small")
+plt.tight_layout()
+plt.show()