Phases and envelopes

In this tutorial, we will introduce Pynapple’s functionality for computing signal phases and envelopes. Most of this functionality is part of the pynapple.process.signal module and is built on the Hilbert transform.

Show code cell source

Hide code cell source

# we'll import the packages we're going to use
import numpy as np
import pynapple as nap
import matplotlib.pyplot as plt
import seaborn as sns

# some configuration, you can ignore this
custom_params = {"axes.spines.right": False, "axes.spines.top": False}
sns.set_theme(
    style="ticks", palette="colorblind", font_scale=1.5, rc=custom_params
)

Extracting the analytic signal

Let us start by simulating a noisy oscillatory signal containing two frequencies:

sampling_rate_hz = 1000
times = np.arange(0, 5, 1 / sampling_rate_hz)

# Low frequency (8 Hz)
low_freq_hz = 8
signal = np.cos(2 * np.pi * low_freq_hz * times)

# High-frequency (40 Hz)
high_freq_hz = 40
segments = [(1.5, 2.0), (3.0, 3.5)]
for start, end in segments:
    mask = (times >= start) & (times <= end)
    signal[mask] = np.cos(2 * np.pi * high_freq_hz * times[mask])

# Add noise
signal = signal + 0.3 * np.random.normal(size=len(times))

# Convert to Tsd
signal = nap.Tsd(t=times, d=signal)

Let’s visualize that:

Show code cell source

Hide code cell source

segment = nap.IntervalSet(1, 4)
plt.figure(figsize=(10,3))
plt.plot(signal.restrict(segment))
plt.xlabel("time (s)")
plt.ylabel("amplitude (a.u.)")
plt.tight_layout();

Now, imagine that we are interested in the low frequency part of the signal. We can start by applying a bandpass filter using apply_bandpass_filter to keep only the relevant frequencies (5-10Hz):

filtered_signal = nap.apply_bandpass_filter(
    signal, (5, 10), fs=sampling_rate_hz, mode="butter"
)
filtered_signal

Time (s)
----------  -----------
0.0         0.431097
0.001       0.424725
0.002       0.417263
0.003       0.408704
0.004       0.399043
0.005       0.388277
0.006       0.376407
...
4.993       7.69622e-05
4.994       0.000141857
4.995       0.000188552
4.996       0.000219716
4.997       0.000237831
4.998       0.000245183
4.999       0.000243858
dtype: float64, shape: (5000,)

Let’s visualize that together with the original signal:

Show code cell source

Hide code cell source

plt.figure(figsize=(10,3))
plt.plot(signal.restrict(segment), label="signal")
plt.plot(filtered_signal.restrict(segment), label="filtered signal")
plt.xlabel("time (s)")
plt.ylabel("amplitude (a.u.)")
plt.legend(loc="upper left", bbox_to_anchor=(1, 1))
plt.tight_layout();

We can now use apply_hilbert_transform to extract the analytic signal:

analytic_signal = nap.apply_hilbert_transform(filtered_signal)
analytic_signal

Time (s)
----------  ----------------------------------------------
0.0         (0.4310970099596455-0.389180608806092j)
0.001       (0.42472509556786964-0.0945978539377379j)
0.002       (0.41726340335953166-0.07713514412161379j)
0.003       (0.4087042274931062+0.032269935027983934j)
0.004       (0.39904266461746524+0.04937495440398126j)
0.005       (0.38827668999479675+0.12155693572414608j)
0.006       (0.37640722631227685+0.13825270206051443j)
...
4.993       (7.696222489530556e-05-0.06078391984351132j)
4.994       (0.00014185676659299133-0.050423101980494246j)
4.995       (0.00018855154093598686-0.09463313216781304j)
4.996       (0.00021971563980446263-0.08290336509174884j)
4.997       (0.0002378306303482077-0.16221293763392533j)
4.998       (0.00024518288525214306-0.14836692668958085j)
4.999       (0.00024385847552243833-0.40822925879990607j)
dtype: complex128, shape: (5000,)

If we visualize the analytic signal with the input signal, you will notice that the analytic signal appears identical to the original signal:

Show code cell source

Hide code cell source

plt.figure(figsize=(10,3))
plt.plot(filtered_signal.restrict(segment), label="filtered signal")
plt.plot(analytic_signal.restrict(segment), label="analytic signal")
plt.xlabel("time (s)")
plt.ylabel("amplitude (a.u.)")
plt.legend(loc="upper left", bbox_to_anchor=(1, 1))
plt.title("nap.apply_hilbert_transform(filtered_signal)")
plt.tight_layout();

/home/runner/.local/lib/python3.12/site-packages/matplotlib/cbook.py:1719: ComplexWarning: Casting complex values to real discards the imaginary part
  return math.isfinite(val)
/home/runner/.local/lib/python3.12/site-packages/matplotlib/cbook.py:1355: ComplexWarning: Casting complex values to real discards the imaginary part
  return np.asarray(x, float)

This happens because the analytic signal is complex-valued. The real part is exactly the input signal. The imaginary part is the Hilbert transform (a 90° phase-shifted version). When you pass a complex signal to matplotlib, it automatically plots only the real part (see the warnings above).

To actually see what’s going on, you can plot the real and imaginary parts separately:

Show code cell source

Hide code cell source

plt.figure(figsize=(10,3))
plt.plot(np.real(analytic_signal).restrict(segment), label="real part")
plt.plot(
    np.imag(analytic_signal).restrict(segment),
    label="imaginary part",
)
plt.xlabel("time (s)")
plt.ylabel("amplitude (a.u.)")
plt.legend(loc="upper left", bbox_to_anchor=(1, 1));

Computing the signal envelope

From the analytic signal, it is often the case that we will compute other things. For one, we can extract the envelope of a signal, by taking the absolute value. To make things easy, Pynapple provides compute_hilbert_envelope to compute the envelope in one go:

envelope = nap.compute_hilbert_envelope(filtered_signal)
envelope

Time (s)
----------  ---------
0.0         0.580781
0.001       0.435132
0.002       0.424333
0.003       0.409976
0.004       0.402086
0.005       0.40686
0.006       0.400994
...
4.993       0.060784
4.994       0.0504233
4.995       0.0946333
4.996       0.0829037
4.997       0.162213
4.998       0.148367
4.999       0.408229
dtype: float64, shape: (5000,)

Visualizing the envelope over the signal:

Show code cell source

Hide code cell source

plt.figure(figsize=(10,3))
plt.plot(filtered_signal.restrict(segment), label="filtered signal")
plt.plot(envelope.restrict(segment), label="envelope", color="red")
plt.xlabel("time (s)")
plt.ylabel("amplitude (a.u.)")
plt.legend(loc="upper left", bbox_to_anchor=(1, 1))
plt.title("nap.compute_hilbert_envelope(filtered_signal)")
plt.tight_layout();

Computing the signal phase

We can also estimate the signal’s phase, by taking angle and wrapping. To make things easy, Pynapple provides compute_hilbert_phase to compute the phase in one go:

phase = nap.compute_hilbert_phase(filtered_signal)
phase

Time (s)
----------  ---------
0.0         5.54884
0.001       6.06404
0.002       6.10039
0.003       0.0787932
0.004       0.123108
0.005       0.303402
0.006       0.351999
...
4.993       4.71366
4.994       4.7152
4.995       4.71438
4.996       4.71504
4.997       4.71386
4.998       4.71404
4.999       4.71299
dtype: float64, shape: (5000,)

Visualizing the phase with the signal:

Show code cell source

Hide code cell source

fig, ax_sig = plt.subplots(figsize=(10, 3))
signal_line = ax_sig.plot(filtered_signal.restrict(segment))
ax_sig.set_ylabel("amplitude (a.u.)")
ax_phase = ax_sig.twinx()
ax_phase.spines["right"].set_visible(True)
phase_line = ax_phase.plot(phase.restrict(segment), color="red", linewidth=0.5)
ax_phase.set_ylabel("phase (rad)")
ax_sig.set_xlabel("time (s)")
ax_sig.legend(
    signal_line + phase_line,
    ["filtered signal", "phase"],
    loc="upper left",
    bbox_to_anchor=(1.15, 1),
)
plt.title("nap.compute_hilbert_phase(filtered_signal)")
plt.tight_layout();

Detecting oscillatory events

Having looked at the low frequency part of our signal, we might also be interested in the high frequency part. To start with, we might simply be interested in finding the epochs where the signal is oscillating at high frequencies. Pynapple provides the detect_oscillatory_events exactly for such a goal.

To get it to work nicely, you will have to the tune the following detection parameters:

• frequency band: the band of frequencies you are interested in

• threshold band: minimum and maximum thresholds to apply to the z-scored envelope of the squared signal

• duration band: minimum and maximum duration of the events

• minimum interval between events

events = nap.detect_oscillatory_events(
    data=signal,
    epochs=signal.time_support,
    frequency_band=(35, 45),
    threshold_band=(1,5),
    duration_band=(0.4,0.6),
    min_interval=0.05,
)
events

  index    start    end    power    amplitude    peak_time
      0    1.503  1.985    -0.28         1.11          1.6
      1    3.011  3.493    -0.39         1.05          3.4
shape: (2, 2), time unit: sec.

We can then visualize the found events on top of the original signal as validation:

Show code cell source

Hide code cell source

plt.figure(figsize=(10, 3))
plt.plot(signal.restrict(segment), label="signal")
first = True
for s, e in events.intersect(segment).values:
    if first:
        plt.axvspan(s, e, color="orange", alpha=0.3, label="event")
        first = False
    else:
        plt.axvspan(s, e, color="orange", alpha=0.3)
plt.xlabel("time (s)")
plt.ylabel("amplitude (a.u.)")
plt.legend(loc="upper left", bbox_to_anchor=(1, 1))
plt.title("nap.detect_oscillatory_events(signal, ...)")
plt.tight_layout();