Skip to content

Commit e220ac7

Browse files
nbaranbara
authored
Real-time (causal) phase estimation algo (nbara#75)
1 parent 9204cfd commit e220ac7

File tree

8 files changed

+1317
-0
lines changed

8 files changed

+1317
-0
lines changed

‎README.md‎

Lines changed: 9 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -153,3 +153,12 @@ and was adapted to python by [Giuseppe Ferraro](mailto:giuseppe.ferraro@isae-sup
153153
EEG Bad Channel Detection Using Local Outlier Factor (LOF). Sensors (Basel).
154154
2022 Sep 27;22(19):7314. https://doi.org/10.3390/s22197314.
155155
```
156+
157+
### 6. Real-Time Phase Estimation
158+
159+
This code is based on the Matlab implementation from [Michael Rosenblum](http://www.stat.physik.uni-potsdam.de/~mros), and its corresponding paper [1].
160+
161+
```sql
162+
[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
163+
164+
```

‎doc/index.rst‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -38,6 +38,7 @@ Here is a list of the methods and techniques available in ``meegkit``:
3838
~meegkit.dss
3939
~meegkit.detrend
4040
~meegkit.lof
41+
~meegkit.phase
4142
~meegkit.ress
4243
~meegkit.sns
4344
~meegkit.star

‎examples/example_phase_estimation.ipynb‎

Lines changed: 122 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,122 @@
1+
{
2+
"cells": [
3+
{
4+
"cell_type": "markdown",
5+
"metadata": {},
6+
"source": [
7+
"\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"
8+
]
9+
},
10+
{
11+
"cell_type": "code",
12+
"execution_count": null,
13+
"metadata": {
14+
"collapsed": false
15+
},
16+
"outputs": [],
17+
"source": [
18+
"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)"
19+
]
20+
},
21+
{
22+
"cell_type": "markdown",
23+
"metadata": {},
24+
"source": [
25+
"## Build data\nFirst, we generate a multi-component signal with amplitude and phase\nmodulations, as described in the paper [1]_.\n\n"
26+
]
27+
},
28+
{
29+
"cell_type": "code",
30+
"execution_count": null,
31+
"metadata": {
32+
"collapsed": false
33+
},
34+
"outputs": [],
35+
"source": [
36+
"npt = 100000\nfs = 100\ns = generate_multi_comp_data(npt, fs) # Generate test data\ndt = 1 / fs\ntime = np.arange(npt) * dt"
37+
]
38+
},
39+
{
40+
"cell_type": "markdown",
41+
"metadata": {},
42+
"source": [
43+
"### Vizualize signal\nPlot the test signal's Fourier spectrum\n\n"
44+
]
45+
},
46+
{
47+
"cell_type": "code",
48+
"execution_count": null,
49+
"metadata": {
50+
"collapsed": false
51+
},
52+
"outputs": [],
53+
"source": [
54+
"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()"
55+
]
56+
},
57+
{
58+
"cell_type": "markdown",
59+
"metadata": {},
60+
"source": [
61+
"### 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"
62+
]
63+
},
64+
{
65+
"cell_type": "code",
66+
"execution_count": null,
67+
"metadata": {
68+
"collapsed": false
69+
},
70+
"outputs": [],
71+
"source": [
72+
"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)"
73+
]
74+
},
75+
{
76+
"cell_type": "markdown",
77+
"metadata": {},
78+
"source": [
79+
"## Results\nHere we reproduce figure 1 from the original paper [1]_.\n\n"
80+
]
81+
},
82+
{
83+
"cell_type": "markdown",
84+
"metadata": {},
85+
"source": [
86+
"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"
87+
]
88+
},
89+
{
90+
"cell_type": "code",
91+
"execution_count": null,
92+
"metadata": {
93+
"collapsed": false
94+
},
95+
"outputs": [],
96+
"source": [
97+
"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()"
98+
]
99+
}
100+
],
101+
"metadata": {
102+
"kernelspec": {
103+
"display_name": "Python 3",
104+
"language": "python",
105+
"name": "python3"
106+
},
107+
"language_info": {
108+
"codemirror_mode": {
109+
"name": "ipython",
110+
"version": 3
111+
},
112+
"file_extension": ".py",
113+
"mimetype": "text/x-python",
114+
"name": "python",
115+
"nbconvert_exporter": "python",
116+
"pygments_lexer": "ipython3",
117+
"version": "3.11.6"
118+
}
119+
},
120+
"nbformat": 4,
121+
"nbformat_minor": 0
122+
}

‎examples/example_phase_estimation.py‎

Lines changed: 137 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,137 @@
1+
"""
2+
Causal phase estimation example
3+
===============================
4+
5+
This example shows how to causally estimate the phase of a signal using two
6+
oscillator models, as described in [1]_.
7+
8+
Uses `meegkit.phase.ResOscillator()` and `meegkit.phase.NonResOscillator()`.
9+
10+
References
11+
----------
12+
.. [1] Rosenblum, M., Pikovsky, A., Kühn, A.A. et al. Real-time estimation
13+
of phase and amplitude with application to neural data. Sci Rep 11, 18037
14+
(2021). https://doi.org/10.1038/s41598-021-97560-5
15+
16+
"""
17+
import os
18+
import sys
19+
20+
import matplotlib.pyplot as plt
21+
import numpy as np
22+
from scipy.signal import hilbert
23+
24+
from meegkit.phase import NonResOscillator, ResOscillator, locking_based_phase
25+
26+
sys.path.append(os.path.join("..", "tests"))
27+
28+
from test_filters import generate_multi_comp_data, phase_difference # noqa:E402
29+
30+
rng = np.random.default_rng(5)
31+
32+
###############################################################################
33+
# Build data
34+
# -----------------------------------------------------------------------------
35+
# First, we generate a multi-component signal with amplitude and phase
36+
# modulations, as described in the paper [1]_.
37+
38+
###############################################################################
39+
npt = 100000
40+
fs = 100
41+
s = generate_multi_comp_data(npt, fs) # Generate test data
42+
dt = 1 / fs
43+
time = np.arange(npt) * dt
44+
45+
###############################################################################
46+
# Vizualize signal
47+
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
48+
# Plot the test signal's Fourier spectrum
49+
f, ax = plt.subplots(2, 1)
50+
ax[0].plot(time, s)
51+
ax[0].set_xlabel("Time (s)")
52+
ax[0].set_title("Test signal")
53+
ax[1].psd(s, Fs=fs, NFFT=2048*4, noverlap=fs)
54+
ax[1].set_title("Test signal's Fourier spectrum")
55+
plt.tight_layout()
56+
57+
###############################################################################
58+
# Compute phase and amplitude
59+
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
60+
# We compute the Hilbert phase and amplitude, as well as the phase and
61+
# amplitude obtained by the locking-based technique, non-resonant and
62+
# resonant oscillator.
63+
ht_ampl = np.abs(hilbert(s)) # Hilbert amplitude
64+
ht_phase = np.angle(hilbert(s)) # Hilbert phase
65+
66+
lb_phase = locking_based_phase(s, dt, npt)
67+
lb_phi_dif = phase_difference(ht_phase, lb_phase)
68+
69+
osc = NonResOscillator(fs, 1.1)
70+
nr_phase, nr_ampl = osc.transform(s)
71+
nr_phase = nr_phase[:, 0]
72+
nr_phi_dif = phase_difference(ht_phase, nr_phase)
73+
74+
osc = ResOscillator(fs, 1.1)
75+
r_phase, r_ampl = osc.transform(s)
76+
r_phase = r_phase[:, 0]
77+
r_phi_dif = phase_difference(ht_phase, r_phase)
78+
79+
80+
###############################################################################
81+
# Results
82+
# -----------------------------------------------------------------------------
83+
# Here we reproduce figure 1 from the original paper [1]_.
84+
85+
###############################################################################
86+
# The first row shows the test signal $s$ and its Hilbert amplitude $a_H$ ; one
87+
# can see that ah does not represent a good envelope for $s$. On the contrary,
88+
# the Hilbert-based phase estimation yields good results, and therefore we take
89+
# it for the ground truth.
90+
# Rows 2-4 show the difference between the Hilbert phase and causally
91+
# estimated phases ($\phi_L$, $\phi_N$, $\phi_R$) are obtained by means of the
92+
# locking-based technique, non-resonant and resonant oscillator, respectively).
93+
# These panels demonstrate that the output of the developed causal algorithms
94+
# is very close to the HT-phase. Notice that we show $\phi_H - \phi_N$
95+
# modulo $2\pi$, since the phase difference is not bounded.
96+
f, ax = plt.subplots(4, 2, sharex=True, sharey=True, figsize=(12, 8))
97+
ax[0, 0].plot(time, s, time, ht_phase, lw=.75)
98+
ax[0, 0].set_ylabel(r"$s,\phi_H$")
99+
ax[0, 0].set_title("Signal and its Hilbert phase")
100+
101+
ax[1, 0].plot(time, lb_phi_dif, lw=.75)
102+
ax[1, 0].axhline(0, color="k", ls=":", zorder=-1)
103+
ax[1, 0].set_ylabel(r"$\phi_H - \phi_L$")
104+
ax[1, 0].set_ylim([-np.pi, np.pi])
105+
ax[1, 0].set_title("Phase locking approach")
106+
107+
ax[2, 0].plot(time, nr_phi_dif, lw=.75)
108+
ax[2, 0].axhline(0, color="k", ls=":", zorder=-1)
109+
ax[2, 0].set_ylabel(r"$\phi_H - \phi_N$")
110+
ax[2, 0].set_ylim([-np.pi, np.pi])
111+
ax[2, 0].set_title("Nonresonant oscillator")
112+
113+
ax[3, 0].plot(time, r_phi_dif, lw=.75)
114+
ax[3, 0].axhline(0, color="k", ls=":", zorder=-1)
115+
ax[3, 0].set_ylim([-np.pi, np.pi])
116+
ax[3, 0].set_ylabel("$\phi_H - \phi_R$")
117+
ax[3, 0].set_xlabel("Time")
118+
ax[3, 0].set_title("Resonant oscillator")
119+
120+
ax[0, 1].plot(time, s, time, ht_ampl, lw=.75)
121+
ax[0, 1].set_ylabel(r"$s,a_H$")
122+
ax[0, 1].set_title("Signal and its Hilbert amplitude")
123+
124+
ax[1, 1].axis("off")
125+
126+
ax[2, 1].plot(time, s, time, nr_ampl, lw=.75)
127+
ax[2, 1].set_ylabel(r"$s,a_N$")
128+
ax[2, 1].set_title("Amplitudes")
129+
ax[2, 1].set_title("Nonresonant oscillator")
130+
131+
ax[3, 1].plot(time, s, time, r_ampl, lw=.75)
132+
ax[3, 1].set_xlabel("Time")
133+
ax[3, 1].set_ylabel(r"$s,a_R$")
134+
ax[3, 1].set_title("Resonant oscillator")
135+
plt.suptitle("Amplitude (right) and phase (left) estimation algorithms")
136+
plt.tight_layout()
137+
plt.show()

0 commit comments

Comments
 (0)