Zero-padding Hilbert

[tristan-deep]

@vincentvdschaft already made a great blog on IQ demodulation. I here want to add a small part showcasing the need for zero-padding the signal before taking the Hilbert transform. We do this in our implementation of envelope detecting ultrasound data. It mitigates the cyclic convolution effect that arises from the FFT's implicit assumption that the signal is periodic.

import matplotlib.pyplot as plt
import numpy as np
from keras import ops

from zea.ops import hilbert
from zea.visualize import set_mpl_style

set_mpl_style()


def envelope_detect(data, axis=-1, pad=True):
    """Adaptation of zea.ops.envelope_detect.

    Envelope detection using Hilbert transform with optional zero padding.

    Args:
        data (np.ndarray): Real-valued signal.
        axis (int): Axis along which to apply the Hilbert transform.
        pad (bool): Whether to apply zero-padding to reduce edge effects.

    Returns:
        np.ndarray: Envelope of the signal (same shape as input).
    """
    n_ax = ops.shape(data)[axis]
    n_ax_float = ops.cast(n_ax, "float32")
    if pad:
        # Calculate next power of 2: M = 2^ceil(log2(n_ax))
        log2_n_ax = ops.log(n_ax_float) / ops.log(2.0)
        M = ops.cast(2 ** ops.ceil(log2_n_ax), "int32")
    else:
        M = n_ax

    data = hilbert(data, N=M, axis=axis)
    indices = ops.arange(n_ax)
    data = ops.take(data, indices, axis=axis)

    # Compute magnitude
    real = ops.real(data)
    imag = ops.imag(data)
    data = ops.sqrt(real**2 + imag**2)
    data = ops.cast(data, "float32")
    return data


fs = 1000
t = np.linspace(0, 1, fs, endpoint=False)
x = np.sin(2 * np.pi * 50 * t)
x[int(0.3 * fs) :] = 0.0  # abrupt cutoff

env_no_pad = envelope_detect(x, pad=False)
env_padded = envelope_detect(x, pad=True)


fig, axs = plt.subplots(3, 1, figsize=(7, 6), sharex=True)

axs[0].plot(t, x, color="gray")
axs[0].set_title("Truncated sine wave")

axs[1].plot(t, env_no_pad, color="orange")
axs[1].set_title("Envelope (no padding)")

axs[2].plot(t, env_padded, color="cyan")
axs[2].set_title("Envelope (with padding)")

axs[-1].set_xlabel("time [s]")
plt.savefig("test_hilbert.png", bbox_inches="tight", dpi=150)