Real-time (causal) phase estimation algo (#75)

README.md

@@ -153,3 +153,12 @@ and was adapted to python by [Giuseppe Ferraro](mailto:giuseppe.ferraro@isae-sup
     EEG Bad Channel Detection Using Local Outlier Factor (LOF). Sensors (Basel). 
     2022 Sep 27;22(19):7314. https://doi.org/10.3390/s22197314.
 ```
+
+### 6. Real-Time 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].
+
+```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
+
+```

doc/index.rst

@@ -38,6 +38,7 @@ Here is a list of the methods and techniques available in ``meegkit``:
    ~meegkit.dss
    ~meegkit.detrend
    ~meegkit.lof
+   ~meegkit.phase
    ~meegkit.ress
    ~meegkit.sns
    ~meegkit.star

examples/example_phase_estimation.ipynb

@@ -0,0 +1,122 @@
+{
+  "cells": [
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "\n# Causal phase estimation example\n\nThis example shows how to causally estimate the phase of a signal using two\noscillator models, as described in [1]_.\n\nUses `meegkit.phase.ResOscillator()` and `meegkit.phase.NonResOscillator()`.\n\n## References\n.. [1] Rosenblum, M., Pikovsky, A., K\u00fchn, A.A. et al. Real-time estimation\n       of phase and amplitude with application to neural data. Sci Rep 11, 18037\n       (2021). https://doi.org/10.1038/s41598-021-97560-5\n"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": null,
+      "metadata": {
+        "collapsed": false
+      },
+      "outputs": [],
+      "source": [
+        "import os\nimport sys\n\nimport matplotlib.pyplot as plt\nimport numpy as np\nfrom scipy.signal import hilbert\n\nfrom meegkit.phase import NonResOscillator, ResOscillator, locking_based_phase\n\nsys.path.append(os.path.join(\"..\", \"tests\"))\n\nfrom test_filters import generate_multi_comp_data, phase_difference  # noqa:E402\n\nrng = np.random.default_rng(5)"
+      ]
+    },
+    {
+      "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": [
+        "npt = 100000\nfs = 100\ns  = generate_multi_comp_data(npt, fs)  # Generate test data\ndt = 1 / fs\ntime = np.arange(npt) * dt"
+      ]
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "### Vizualize signal\nPlot the test signal's Fourier spectrum\n\n"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": null,
+      "metadata": {
+        "collapsed": false
+      },
+      "outputs": [],
+      "source": [
+        "f, ax = plt.subplots(2, 1)\nax[0].plot(time, s)\nax[0].set_xlabel(\"Time (s)\")\nax[0].set_title(\"Test signal\")\nax[1].psd(s, Fs=fs, NFFT=2048*4, noverlap=fs)\nax[1].set_title(\"Test signal's Fourier spectrum\")\nplt.tight_layout()"
+      ]
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "### Compute phase and amplitude\nWe compute the Hilbert phase and amplitude, as well as the phase and\namplitude obtained by the locking-based technique, non-resonant and\nresonant oscillator.\n\n"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": null,
+      "metadata": {
+        "collapsed": false
+      },
+      "outputs": [],
+      "source": [
+        "ht_ampl = np.abs(hilbert(s))  # Hilbert amplitude\nht_phase = np.angle(hilbert(s))  # Hilbert phase\n\nlb_phase = locking_based_phase(s, dt, npt)\nlb_phi_dif = phase_difference(ht_phase, lb_phase)\n\nosc = NonResOscillator(fs, 1.1)\nnr_phase, nr_ampl = osc.transform(s)\nnr_phase = nr_phase[:, 0]\nnr_phi_dif = phase_difference(ht_phase, nr_phase)\n\nosc = ResOscillator(fs, 1.1)\nr_phase, r_ampl = osc.transform(s)\nr_phase = r_phase[:, 0]\nr_phi_dif = phase_difference(ht_phase, r_phase)"
+      ]
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "## Results\nHere we reproduce figure 1 from the original paper [1]_.\n\n"
+      ]
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "The first row shows the test signal $s$ and its Hilbert amplitude $a_H$ ; one\ncan see that ah does not represent a good envelope for $s$. On the contrary,\nthe Hilbert-based phase estimation yields good results, and therefore we take\nit for the ground truth.\nRows 2-4 show the difference between the Hilbert phase and causally\nestimated phases ($\\phi_L$, $\\phi_N$, $\\phi_R$) are obtained by means of the\nlocking-based technique, non-resonant and resonant oscillator, respectively).\nThese panels demonstrate that the output of the developed causal algorithms\nis very close to the HT-phase. Notice that we show $\\phi_H - \\phi_N$\nmodulo $2\\pi$, since the phase difference is not bounded.\n\n"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": null,
+      "metadata": {
+        "collapsed": false
+      },
+      "outputs": [],
+      "source": [
+        "f, ax = plt.subplots(4, 2, sharex=True, sharey=True, figsize=(12, 8))\nax[0, 0].plot(time, s, time, ht_phase, lw=.75)\nax[0, 0].set_ylabel(r\"$s,\\phi_H$\")\nax[0, 0].set_title(\"Signal and its Hilbert phase\")\n\nax[1, 0].plot(time, lb_phi_dif, lw=.75)\nax[1, 0].axhline(0, color=\"k\", ls=\":\", zorder=-1)\nax[1, 0].set_ylabel(r\"$\\phi_H - \\phi_L$\")\nax[1, 0].set_ylim([-np.pi, np.pi])\nax[1, 0].set_title(\"Phase locking approach\")\n\nax[2, 0].plot(time, nr_phi_dif, lw=.75)\nax[2, 0].axhline(0, color=\"k\", ls=\":\", zorder=-1)\nax[2, 0].set_ylabel(r\"$\\phi_H - \\phi_N$\")\nax[2, 0].set_ylim([-np.pi, np.pi])\nax[2, 0].set_title(\"Nonresonant oscillator\")\n\nax[3, 0].plot(time, r_phi_dif, lw=.75)\nax[3, 0].axhline(0, color=\"k\", ls=\":\", zorder=-1)\nax[3, 0].set_ylim([-np.pi, np.pi])\nax[3, 0].set_ylabel(\"$\\phi_H - \\phi_R$\")\nax[3, 0].set_xlabel(\"Time\")\nax[3, 0].set_title(\"Resonant oscillator\")\n\nax[0, 1].plot(time, s, time, ht_ampl, lw=.75)\nax[0, 1].set_ylabel(r\"$s,a_H$\")\nax[0, 1].set_title(\"Signal and its Hilbert amplitude\")\n\nax[1, 1].axis(\"off\")\n\nax[2, 1].plot(time, s, time, nr_ampl, lw=.75)\nax[2, 1].set_ylabel(r\"$s,a_N$\")\nax[2, 1].set_title(\"Amplitudes\")\nax[2, 1].set_title(\"Nonresonant oscillator\")\n\nax[3, 1].plot(time, s, time, r_ampl, lw=.75)\nax[3, 1].set_xlabel(\"Time\")\nax[3, 1].set_ylabel(r\"$s,a_R$\")\nax[3, 1].set_title(\"Resonant oscillator\")\nplt.suptitle(\"Amplitude (right) and phase (left) estimation algorithms\")\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.11.6"
+    }
+  },
+  "nbformat": 4,
+  "nbformat_minor": 0
+}

examples/example_phase_estimation.py

@@ -0,0 +1,137 @@
+"""
+Causal phase estimation example
+===============================
+
+This example shows how to causally estimate the phase of a signal using two
+oscillator models, as described in [1]_.
+
+Uses `meegkit.phase.ResOscillator()` and `meegkit.phase.NonResOscillator()`.
+
+References
+----------
+.. [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
+
+"""
+import os
+import sys
+
+import matplotlib.pyplot as plt
+import numpy as np
+from scipy.signal import hilbert
+
+from meegkit.phase import NonResOscillator, ResOscillator, locking_based_phase
+
+sys.path.append(os.path.join("..", "tests"))
+
+from test_filters import generate_multi_comp_data, phase_difference  # noqa:E402
+
+rng = np.random.default_rng(5)
+
+###############################################################################
+# Build data
+# -----------------------------------------------------------------------------
+# First, we generate a multi-component signal with amplitude and phase
+# modulations, as described in the paper [1]_.
+
+###############################################################################
+npt = 100000
+fs = 100
+s  = generate_multi_comp_data(npt, fs)  # Generate test data
+dt = 1 / fs
+time = np.arange(npt) * dt
+
+###############################################################################
+# Vizualize signal
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+# Plot the test signal's Fourier spectrum
+f, ax = plt.subplots(2, 1)
+ax[0].plot(time, s)
+ax[0].set_xlabel("Time (s)")
+ax[0].set_title("Test signal")
+ax[1].psd(s, Fs=fs, NFFT=2048*4, noverlap=fs)
+ax[1].set_title("Test signal's Fourier spectrum")
+plt.tight_layout()
+
+###############################################################################
+# Compute phase and amplitude
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+# We compute the Hilbert phase and amplitude, as well as the phase and
+# amplitude obtained by the locking-based technique, non-resonant and
+# resonant oscillator.
+ht_ampl = np.abs(hilbert(s))  # Hilbert amplitude
+ht_phase = np.angle(hilbert(s))  # Hilbert phase
+
+lb_phase = locking_based_phase(s, dt, npt)
+lb_phi_dif = phase_difference(ht_phase, lb_phase)
+
+osc = NonResOscillator(fs, 1.1)
+nr_phase, nr_ampl = osc.transform(s)
+nr_phase = nr_phase[:, 0]
+nr_phi_dif = phase_difference(ht_phase, nr_phase)
+
+osc = ResOscillator(fs, 1.1)
+r_phase, r_ampl = osc.transform(s)
+r_phase = r_phase[:, 0]
+r_phi_dif = phase_difference(ht_phase, r_phase)
+
+
+###############################################################################
+# Results
+# -----------------------------------------------------------------------------
+# Here we reproduce figure 1 from the original paper [1]_.
+
+###############################################################################
+# The first row shows the test signal $s$ and its Hilbert amplitude $a_H$ ; one
+# can see that ah does not represent a good envelope for $s$. On the contrary,
+# the Hilbert-based phase estimation yields good results, and therefore we take
+# it for the ground truth.
+# Rows 2-4 show the difference between the Hilbert phase and causally
+# estimated phases ($\phi_L$, $\phi_N$, $\phi_R$) are obtained by means of the
+# locking-based technique, non-resonant and resonant oscillator, respectively).
+# These panels demonstrate that the output of the developed causal algorithms
+# is very close to the HT-phase. Notice that we show $\phi_H - \phi_N$
+# modulo $2\pi$, since the phase difference is not bounded.
+f, ax = plt.subplots(4, 2, sharex=True, sharey=True, figsize=(12, 8))
+ax[0, 0].plot(time, s, time, ht_phase, lw=.75)
+ax[0, 0].set_ylabel(r"$s,\phi_H$")
+ax[0, 0].set_title("Signal and its Hilbert phase")
+
+ax[1, 0].plot(time, lb_phi_dif, lw=.75)
+ax[1, 0].axhline(0, color="k", ls=":", zorder=-1)
+ax[1, 0].set_ylabel(r"$\phi_H - \phi_L$")
+ax[1, 0].set_ylim([-np.pi, np.pi])
+ax[1, 0].set_title("Phase locking approach")
+
+ax[2, 0].plot(time, nr_phi_dif, lw=.75)
+ax[2, 0].axhline(0, color="k", ls=":", zorder=-1)
+ax[2, 0].set_ylabel(r"$\phi_H - \phi_N$")
+ax[2, 0].set_ylim([-np.pi, np.pi])
+ax[2, 0].set_title("Nonresonant oscillator")
+
+ax[3, 0].plot(time, r_phi_dif, lw=.75)
+ax[3, 0].axhline(0, color="k", ls=":", zorder=-1)
+ax[3, 0].set_ylim([-np.pi, np.pi])
+ax[3, 0].set_ylabel("$\phi_H - \phi_R$")
+ax[3, 0].set_xlabel("Time")
+ax[3, 0].set_title("Resonant oscillator")
+
+ax[0, 1].plot(time, s, time, ht_ampl, lw=.75)
+ax[0, 1].set_ylabel(r"$s,a_H$")
+ax[0, 1].set_title("Signal and its Hilbert amplitude")
+
+ax[1, 1].axis("off")
+
+ax[2, 1].plot(time, s, time, nr_ampl, lw=.75)
+ax[2, 1].set_ylabel(r"$s,a_N$")
+ax[2, 1].set_title("Amplitudes")
+ax[2, 1].set_title("Nonresonant oscillator")
+
+ax[3, 1].plot(time, s, time, r_ampl, lw=.75)
+ax[3, 1].set_xlabel("Time")
+ax[3, 1].set_ylabel(r"$s,a_R$")
+ax[3, 1].set_title("Resonant oscillator")
+plt.suptitle("Amplitude (right) and phase (left) estimation algorithms")
+plt.tight_layout()
+plt.show()