|
| 1 | +""" |
| 2 | +Endpoint-corrected Hilbert transform (ECHT) phase estimation |
| 3 | +============================================================ |
| 4 | +
|
| 5 | +This example shows how to causally estimate the phase of a signal using the |
| 6 | +Endpoint-corrected Hilbert transform (ECHT) [1]_. |
| 7 | +
|
| 8 | +Uses `meegkit.phase.ECHT()`. |
| 9 | +
|
| 10 | +References |
| 11 | +---------- |
| 12 | +.. [1] Schreglmann, S. R., Wang, D., Peach, R. L., Li, J., Zhang, X., Latorre, |
| 13 | + A., ... & Grossman, N. (2021). Non-invasive suppression of essential tremor |
| 14 | + via phase-locked disruption of its temporal coherence. Nature |
| 15 | + communications, 12(1), 363. |
| 16 | +
|
| 17 | +""" |
| 18 | +import matplotlib.pyplot as plt |
| 19 | +import numpy as np |
| 20 | +from scipy.signal import hilbert |
| 21 | + |
| 22 | +from meegkit.phase import ECHT |
| 23 | + |
| 24 | +rng = np.random.default_rng(38872) |
| 25 | + |
| 26 | +plt.rcParams["axes.grid"] = True |
| 27 | +plt.rcParams["grid.linestyle"] = ":" |
| 28 | + |
| 29 | +############################################################################### |
| 30 | +# Build data |
| 31 | +# ----------------------------------------------------------------------------- |
| 32 | +# First, we generate a multi-component signal with amplitude and phase |
| 33 | +# modulations, as described in the paper [1]_. |
| 34 | +f0 = 2 |
| 35 | + |
| 36 | +N = 500 |
| 37 | +sfreq = 100 |
| 38 | +time = np.linspace(0, N / sfreq, N) |
| 39 | +X = np.cos(2 * np.pi * f0 * time - np.pi / 4) |
| 40 | +phase_true = np.angle(hilbert(X)) |
| 41 | +X += rng.normal(0, 0.5, N) # Add noise |
| 42 | + |
| 43 | +############################################################################### |
| 44 | +# Compute phase and amplitude |
| 45 | +# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
| 46 | +# We compute the Hilbert phase, as well as the phase obtained with the ECHT |
| 47 | +# filter. |
| 48 | +phase_hilbert = np.angle(hilbert(X)) # Hilbert phase |
| 49 | + |
| 50 | +# Compute ECHT-filtered signal |
| 51 | +filt_BW = f0 / 2 |
| 52 | +l_freq = f0 - filt_BW / 2 |
| 53 | +h_freq = f0 + filt_BW / 2 |
| 54 | +echt = ECHT(l_freq, h_freq, sfreq) |
| 55 | + |
| 56 | +Xf = echt.fit_transform(X) |
| 57 | +phase_echt = np.angle(Xf) |
| 58 | + |
| 59 | +############################################################################### |
| 60 | +# Visualize signal |
| 61 | +# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
| 62 | +# Here we plot the original signal, its Fourier spectrum, and the phase obtained |
| 63 | +# with the Hilbert transform and the ECHT filter. The ECHT filter provides a |
| 64 | +# much smoother phase estimate than the Hilbert transform |
| 65 | +fig, ax = plt.subplots(3, 1, figsize=(8, 6)) |
| 66 | +ax[0].plot(time, X) |
| 67 | +ax[0].set_xlabel("Time (s)") |
| 68 | +ax[0].set_title("Test signal") |
| 69 | +ax[0].set_ylabel("Amplitude") |
| 70 | +ax[1].psd(X, Fs=sfreq, NFFT=2048*4, noverlap=sfreq) |
| 71 | +ax[1].set_ylabel("PSD (dB/Hz)") |
| 72 | +ax[1].set_title("Test signal's Fourier spectrum") |
| 73 | +ax[2].plot(time, phase_true, label="True phase", ls=":") |
| 74 | +ax[2].plot(time, phase_echt, label="ECHT phase", lw=.5, alpha=.8) |
| 75 | +ax[2].plot(time, phase_hilbert, label="Hilbert phase", lw=.5, alpha=.8) |
| 76 | +ax[2].set_title("Phase") |
| 77 | +ax[2].set_ylabel("Amplitude") |
| 78 | +ax[2].set_xlabel("Time (s)") |
| 79 | +ax[2].legend(loc="upper right", fontsize="small") |
| 80 | +plt.tight_layout() |
| 81 | +plt.show() |
0 commit comments