meegkit.phase#

Real-time phase and amplitude estimation using resonant oscillators.

Notes

While the original paper [1] only describes the algorithm for a single channel, this implementation can handle multiple channels. However, the algorithm is not optimized for this purpose. Therefore, this code is likely to be slow for large input arrays (n_channels >> 10), since an individual oscillator is instantiated for each channel.

[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

Classes

NonResOscillator([fs, nu, alpha_a, alpha_p, ...])

Real-time measurement of phase and amplitude using non-resonant oscillator.

ResOscillator([fs, nu, update_factor, ...])

Real-time measurement of phase and amplitude using a resonant oscillator.

Device(om0, dt, damping, type[, mu])

Measurement device.

ECHT(l_freq, h_freq, sfreq[, n_fft, filt_order])

Endpoint Corrected Hilbert Transform (ECHT).

Functions

locking_based_phase(s, dt, npt)

Compute the locking-based phase.

rk(phi, dt, n_steps, omega, epsilon, s_prev, ...)

Runge-Kutta method for phase calculation.

init_coefs(om0, dt, alpha)

Compute coefficients for solving oscillator's equations.

one_step_oscillator(x, xd, alpha, eta, ...)

Perform one step of the oscillator equations.

one_step_integrator(z, edelmu, mu, dt, spp, ...)

Perform one step of the integrator in the oscillator equations.
class meegkit.phase.Device(om0, dt, damping, type, mu=500)#

Bases: object

Measurement device.

This class implements a measurement device, as describes by Rosenblum et al.

Parameters:

  • om0 (float) – Oscillator frequency (estimation).

  • dt (float) – Sampling interval.

  • damping (float) – Damping parameter.

  • type (str in {"integrator", "oscillator"}) – Type of the device.

  • mu (float) – Integrator parameter (only for type=”integrator”). Should be >> target frequency of interest.

__init__(om0, dt, damping, type, mu=500)#
init_coefs(om0, dt, damping)#

Initialize coefficients for solving oscillator’s equations.

step(sprev, s, snew)#

Perform one step of the oscillator equations.

Parameters:

  • x (float) – State variable x.

  • y (float) – State variable y.

  • sprev (float) – Previous measurement.

  • s (float) – Current measurement.

  • snew (float) – New measurement.

Returns:

  • x (float) – Updated state variable x.

  • y (float) – Updated state variable y.

class meegkit.phase.ECHT(l_freq, h_freq, sfreq, n_fft=None, filt_order=2)#

Bases: object

Endpoint Corrected Hilbert Transform (ECHT).

See [1] for details.

Parameters:

  • X (ndarray, shape=(n_samples, n_channels)) – Time domain signal.

  • l_freq (float | None) – Low-cutoff frequency of a bandpass causal filter. If None, the data is only low-passed.

  • h_freq (float | None) – High-cutoff frequency of a bandpass causal filter. If None, the data is only high-passed.

  • sfreq (float) – Sampling rate of time domain signal.

  • n_fft (int, optional) – Length of analytic signal. If None, it defaults to the length of X.

  • filt_order (int, optional) – Order of the filter. Default is 2.

Notes

One common implementation of the Hilbert Transform uses a DFT (aka FFT) as part of its computation. Inherent to the DFT is the assumption that a finite sample of a signal is replicated infinitely in time, effectively abutting the end of a sample with its replicated start. If the start and end of the sample are not continuous with each other, distortions are introduced by the DFT. Echt effectively smooths out this ‘discontinuity’ by selectively deforming the start of the sample. It is hence most suited for real-time applications in which the point/s of interest is/are the most recent one/s (i.e. last) in the sample window.

We found that a filter bandwidth (BW=h_freq-l_freq) of up to half the signal’s central frequency works well.

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.

Examples

>>> f0 = 2
>>> filt_BW = f0 / 2
>>> N = 1000
>>> sfreq = N / (2 * np.pi)
>>> t = np.arange(-2 * np.pi, 0, 1 / sfreq)
>>> X = np.cos(2 * np.pi * f0 * t - np.pi / 4)
>>> l_freq = f0 - filt_BW / 2
>>> h_freq = f0 + filt_BW / 2
>>> Xf = echt(X, l_freq, h_freq, sfreq)
__init__(l_freq, h_freq, sfreq, n_fft=None, filt_order=2)#
fit(X, y=None)#

Fit the ECHT transform to the input signal.

Parameters:

X (ndarray, shape=(n_samples, n_channels)) – The input signal to be transformed.

fit_transform(X, y=None)#

Fit the ECHT transform to the input signal and transform it.

transform(X)#

Apply the ECHT transform to the input signal.

Parameters:

X (ndarray, shape=(n_samples, n_channels)) – The input signal to be transformed.

Returns:

Xf – The transformed signal (complex-valued).

Return type:

ndarray, shape=(n_samples, n_channels)

class meegkit.phase.NonResOscillator(fs=250, nu=1.1, alpha_a=6.0, alpha_p=0.2, update_factor=5)#

Bases: object

Real-time measurement of phase and amplitude using non-resonant oscillator.

This estimator relies on the resonance effect. The measuring “device” consists of two linear damped oscillators. The oscillators’ frequency is much larger than the frequency of the signal, i.e., the system is far from resonance. We choose the damping parameters to ensure that (i) the phase of the first linear oscillator equals that of the input and that (ii) amplitude of the second one and the input relate by a known constant multiplicator. The technique yields both phase and amplitude of the input signal.

This estimator includes an automated frequency-tuning algorithm to adjust to the a priori unknown signal frequency.

For a detailed description, refer to [1].

Parameters:

  • fs (float) – Sampling frequency in Hz.

  • nu (float) – Rough estimate of the main frequency of interest.

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

__init__(fs=250, nu=1.1, alpha_a=6.0, alpha_p=0.2, update_factor=5)#
transform(X)#

Transform the input signal into phase and amplitude estimates.

Parameters:

X (ndarray, shape=(n_samples, n_channels)) – The input signal to be transformed.

Returns:

  • phase (ndarray, shape=(n_samples, n_channels)) – Current phase estimate.

  • ampl (ndarray, shape=(n_samples, n_channels)) – Current amplitude estimate.

class meegkit.phase.ResOscillator(fs=1000, nu=4.5, update_factor=5, freq_adaptation=True, assume_centered=False)#

Bases: object

Real-time measurement of phase and amplitude using a resonant oscillator.

In this version, the corresponding “device” consists of a linear oscillator in resonance with the measured signal and of an integrating unit. It also provides both phase and amplitude of the input signal. The method exploits the known relation between the resonant oscillator’s phase and amplitude and those of the input. Additionally, the resonant oscillator acts as a bandpass filter for experimental data.

This filter also includes a frequency-adaptation algorithm, and removes low-frequency trend (baseline fluctuations). For a detailed description, refer to [1].

Parameters:

  • npt (int) – Number of points to be measured and processed

  • fs (float) – Sampling frequency in Hz.

  • nu (float) – Rough estimate of the tremor frequency.

  • freq_adaptation (bool) – Whether to use the frequency adaptation algorithm (default=True).

Returns:

  • Dphase (numpy.ndarray) – Array containing the phase values.

  • Dampl (numpy.ndarray) – Array containing the amplitude values.

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

__init__(fs=1000, nu=4.5, update_factor=5, freq_adaptation=True, assume_centered=False)#
transform(X)#

Transform the input signal into phase and amplitude estimates.

Parameters:

X (ndarray, shape (n_samples, n_channels)) – The input signal to be transformed.

Returns:

  • phase (float) – Current phase estimate.

  • ampl (float) – Current amplitude estimate.

meegkit.phase.init_coefs(om0, dt, alpha)#

Compute coefficients for solving oscillator’s equations.

Parameters:

  • om0 (float) – Oscillator frequency (estimation).

  • dt (float) – Sampling interval.

  • alpha (float) – Half of the damping.

Returns:

  • C1 (float) – Coefficient C1.

  • C2 (float) – Coefficient C2.

  • C3 (float) – Coefficient C3.

  • eetadel (complex) – Exponential term for amplitude device.

  • ealdel (float) – Exponential term for amplitude device.

  • eta (float) – Square root of the difference of oscillator frequency squared and damping squared.

meegkit.phase.locking_based_phase(s, dt, npt)#

Compute the locking-based phase.

This technique exploits the ideas from the synchronization theory. It is well-known that a force s(t) acting on a limit-cycle oscillator can entrain it. It means that the oscillator’s frequency becomes equal to that of the force, and their phases differ by a constant. Thus, the phase of the locked limit-cycle oscillator will correspond to the phase of the signal. For our purposes, it is helpful to use the simplest oscillator model, the so-called phase oscillator. To ensure the phase-locking to the force, we have to adjust the oscillator’s frequency to the signal’s frequency. We assume that we do not know the latter a priori, but can only roughly estimate it. We propose a simple approach that starts with this estimate and automatically tunes the “device’s” frequency to ensure the locking and thus provides the instantaneous phase.

Note that the amplitude is not determined with this approach. Technically, the algorithm reduces to solving differential equation incorporating measured data given at discrete time points; for details, see [1].

Parameters:

  • s (array-like) – Input signal.

  • dt (float) – Time step.

  • npt (int) – Number of points.

Returns:

lb_phase – Locking-based phase.

Return type:

ndarray

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

meegkit.phase.one_step_integrator(z, edelmu, mu, dt, spp, sp, s)#

Perform one step of the integrator in the oscillator equations.

Parameters:

  • z (float) – State variable z.

  • edelmu (float) – Exponential term for integrator.

  • mu (float) – Integrator parameter.

  • dt (float) – Sampling interval.

  • spp (float) – Second previous measurement.

  • sp (float) – Previous measurement.

  • s (float) – Current measurement.

Returns:

z – Updated state variable z.

Return type:

float

meegkit.phase.one_step_oscillator(x, xd, alpha, eta, eeta_del, eal_del, C1, C2, C3, spp, sp, s)#

Perform one step of the oscillator equations.

Parameters:

  • x (float) – State variable x.

  • xd (float) – Derivative of state variable x.

  • alpha (float) – Damping parameter.

  • eta (float) – Square root of the difference of oscillator frequency squared and damping squared.

  • eeta_del (complex) – Exponential term for time step.

  • eal_del (float) – Exponential term for time step.

  • C1 (float) – Coefficient C1.

  • C2 (float) – Coefficient C2.

  • C3 (float) – Coefficient C3.

  • spp (float) – Second orevious measurement.

  • sp (float) – Previous measurement.

  • s (float) – Current measurement.

Returns:

  • x (float) – Updated state variable x.

  • xd (float) – Updated derivative of state variable x.

meegkit.phase.rk(phi, dt, n_steps, omega, epsilon, s_prev, s, s_new)#

Runge-Kutta method for phase calculation.

Parameters:

  • phi (float) – Initial phase value.

  • dt (float) – Time step size.

  • n_steps (int) – Number of steps to iterate.

  • omega (float) – Angular frequency.

  • epsilon (float) – Amplitude of the oscillation.

  • s_prev (float) – Previous phase value.

  • s (float) – Current phase value.

  • s_new (float) – New phase value.

Returns:

phi – The calculated phase value.

Return type:

float