Note

This page was generated from a Jupyter Notebook. Download the notebook (.ipynb)

# Skipped in CI: Colab/bootstrap dependency install cell.

HHT: Analysis

This tutorial explains the procedure for Hilbert-Huang Transform (HHT) analysis using gwexpy. HHT is well-suited for analyzing nonlinear and non-stationary signals, and is a powerful method for extracting local frequency variations (instantaneous frequency) that cannot be captured by traditional Fourier or wavelet transforms. The sample code in this tutorial uses TimeSeries.hht().

HHT consists of two main steps:

1. Empirical Mode Decomposition (EMD): Decomposes the signal into Intrinsic Mode Functions (IMFs).

2. Hilbert Spectral Analysis: Applies the Hilbert transform to each IMF to calculate instantaneous amplitude and instantaneous frequency.

Note: To use this feature, the PyEMD (EMD-signal) package must be installed.

pip install EMD-signal

import astropy.units as u
import matplotlib.pyplot as plt
import numpy as np

from gwexpy.timeseries import TimeSeries

# Check PyEMD installation
try:
    import PyEMD

    print("PyEMD is installed and ready.")
except ImportError:
    raise ImportError(
        "This tutorial requires 'PyEMD' (EMD-signal). Please run: pip install EMD-signal"
    )

PyEMD is installed and ready.

1. Theoretical Background and Advantages of HHT

What is HHT (Hilbert-Huang Transform)?

HHT is a powerful method for analyzing nonlinear and non-stationary signal data. It consists of two main processes:

1. EMD (Empirical Mode Decomposition): Decomposes the original signal into a sum of Intrinsic Mode Functions (IMFs) with different frequency scales.

2. HSA (Hilbert Spectral Analysis): Applies the Hilbert transform to each IMF to calculate Instantaneous Frequency and instantaneous amplitude.

Differences from Traditional STFT

Short-Time Fourier Transform (STFT) and wavelet transforms have a trade-off between time resolution and frequency resolution (limited by the uncertainty principle). In contrast, HHT adaptively generates basis functions based on local time scales, so it is not subject to this limitation and can achieve high time and frequency resolution simultaneously.

EEMD (Ensemble EMD)

In standard EMD, “mode mixing” can occur where different frequency components are mixed within a single IMF due to intermittent signals or noise. To solve this, EEMD (Ensemble EMD) has been proposed.

• Generate multiple ensemble members by adding white noise to the original signal.

• Perform EMD on each signal.

• Take the ensemble average of the resulting IMFs to cancel out the noise and stably extract physically meaningful modes.

Physical intent: in CCSN-like signals we care about whether an extracted ridge tracks a real mode family, such as SASI or a proto-neutron-star oscillation, rather than a basis artefact created by a fixed window.

2. Simulation: Reproducing SASI Modes

In gravitational wave analysis of core-collapse supernovae (CCSN), the instability called SASI (Standing Accretion Shock Instability) is an important observation target. Here, we create a Sine-Gaussian signal simulating gravitational waves from SASI and analyze it (Reference: Takeda et al. 2021).

The Sine-Gaussian signal is expressed by the following equation:

\[h(t) = A \exp\left( - \left( \frac{2\pi (t - t_0)}{t_{\text{width}}} \right)^2 \right) \sin(2\pi f_0 t)\]

import warnings
warnings.filterwarnings("ignore", category=UserWarning)
warnings.filterwarnings("ignore", category=DeprecationWarning)

import numpy as np
import matplotlib.pyplot as plt
from gwexpy.timeseries import TimeSeries


def sine_gaussian(t, t0, t_width, f0, amp=1.0):
    exponent = -(((2 * np.pi * (t - t0)) / t_width) ** 2)
    return amp * np.exp(exponent) * np.sin(2 * np.pi * f0 * t)


# Parameter settings (mimicking SASI mode)
duration = 0.5  # seconds
fs = 4096  # sampling rate
t = np.linspace(0, duration, int(duration * fs))

# Generate signal
t0 = 0.25  # peak time
t_width = 0.4  # width
f0 = 130.0  # frequency (Hz)
clean_signal = sine_gaussian(t, t0, t_width, f0, amp=1e-21)

# Add noise
np.random.seed(42)
noise_amp = 0.1e-21
noise = np.random.normal(0, noise_amp, len(t))
data = clean_signal + noise

# Create TimeSeries object
ts_data = TimeSeries(
    data, t0=0, sample_rate=fs, name="Simulated SASI Signal", unit="strain"
)

# Plot
fig, ax = plt.subplots(figsize=(10, 4))
ax.plot(t, ts_data.value, label="Noisy Signal", color="lightgray")
ax.plot(t, clean_signal, label="Target Signal (SASI)", color="tab:blue")
ax.set_xlabel("Time [s]")
ax.set_ylabel("Strain")
ax.set_xlim(0, duration)
ax.legend()
plt.title("Simulated SASI Gravitational Wave Signal")
plt.show()

3. Comparative Analysis: STFT vs HHT

We apply STFT (spectrogram) and HHT (Hilbert spectrum) to the same signal and compare their resolution differences.

Note: In the following results, STFT may appear more robust to noise and the signal may look cleaner. This is because standard EMD is susceptible to noise and causes “mode mixing.” This problem is dramatically improved by using EEMD, which we will introduce in the next section, but let’s first check the basic behavior of HHT.

Common mistake: reading every sharp HHT ridge as a physical mode. If the ridge fragments when you perturb the padding, decomposition method, or IMF count, it is more likely a decomposition artefact than a robust signal feature.

# --- 1. STFT (Spectrogram) using gwexpy API ---
# Convert sample counts to time-based parameters for gwexpy API:
# nperseg=128 samples → fftlength = 128/fs seconds
# noverlap=120 samples → overlap = 120/fs seconds
# spectrogram2 computes an STFT-like spectrogram (no stride needed)
spec = ts_data.spectrogram2(fftlength=128/fs, overlap=120/fs)
f_stft = spec.frequencies.value
t_stft = spec.times.value
Sxx = spec.value

# --- 2. HHT (EMD + Hilbert) ---
# Important: EMD has numerical tolerance, so it may not decompose well for very small amplitudes (e.g., 1e-21 for GW).
# Therefore, we normalize the data first, then scale back the results.
data_std = np.std(ts_data.value)
ts_norm = TimeSeries(
    ts_data.value / data_std,
    t0=ts_data.t0,
    sample_rate=ts_data.sample_rate,
    name="Normalized SASI Signal",
    unit=ts_data.unit,
)

hht_result = ts_norm.hht(
    emd_method="emd",
    emd_kwargs={
        "sift_max_iter": 200,
        "stopping_criterion": 0.2,
    },
    hilbert_kwargs={
        "pad": 200,
        "if_smooth": 11,
    },
    output="dict",
)

imfs = hht_result["imfs"]
if_dict = hht_result["if"]
ia_dict = hht_result["ia"]
n_imfs = len(imfs)

print(f"Extracted {n_imfs} IMFs from normalized data.")

# --- Plot ---
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(16, 6))

# STFT
mesh = ax1.pcolormesh(t_stft, f_stft, Sxx.T, shading="gouraud", cmap="viridis")
ax1.set_ylabel("Frequency [Hz]")
ax1.set_xlabel("Time [sec]")
ax1.set_title("Short-Time Fourier Transform (STFT)")
ax1.set_ylim(0, 400)
ax1.set_xlim(0, duration)
plt.colorbar(mappable=mesh, ax=ax1, label="Power")

# HHT Plot
# Filter out low-amplitude points (noise) for visibility
amp_list = [np.max(ia_dict[key].value) for key in imfs.keys()]
max_amp = (np.max(amp_list) * data_std) if amp_list else 1.0
threshold = max_amp * 0.1  # Display only points above 10% of max amplitude

sc = None

for key in imfs.keys():
    t_vals = if_dict[key].times.value
    inst_freq = if_dict[key].value
    inst_amp = ia_dict[key].value * data_std
    mask = inst_amp > threshold
    if np.any(mask):
        sc = ax2.scatter(
            t_vals[mask],
            inst_freq[mask],
            c=inst_amp[mask],
            cmap="viridis",
            s=10,
            alpha=0.8,
            vmin=0,
            vmax=max_amp,
        )

ax2.set_ylabel("Instantaneous Frequency [Hz]")
ax2.set_xlabel("Time [sec]")
ax2.set_title("Hilbert-Huang Transform (HHT)")
ax2.set_ylim(0, 400)
ax2.set_xlim(0, duration)
ax2.grid(alpha=0.3)

if sc is not None:
    cbar = plt.colorbar(mappable=sc, ax=ax2)
    cbar.set_label("Instantaneous Amplitude")

plt.tight_layout()
plt.show()

Extracted 8 IMFs from normalized data.

4. Application: Mode Extraction and Parameter Estimation with EEMD

Applying EEMD

In the simple EMD results above, the frequency may oscillate due to the influence of noise, or the original 130Hz signal may be split into multiple IMFs (mode mixing). Using EEMD (Ensemble EMD) can improve this.

In Sasaoka et al. (2024), it was shown that in gravitational wave analysis from CCSN, using EEMD (or CEEMD) can accurately extract the oscillation modes of the proto-neutron star (PNS) and estimate PNS physical quantities (mass-to-radius ratio \(M/R^2\)) from the frequency evolution.

The following code is an example of running EEMD using gwexpy’s TimeSeries.hht().

Failure mode to watch for: increasing eemd_trials can stabilise IMFs, but it does not rescue a poor band selection or a signal dominated by broadband transients. If the target physics is a narrow drifting mode, crop to a plausible time/frequency region first instead of hoping EEMD will separate everything automatically.

# Run EEMD (gwexpy API)

print("Running EEMD on normalized data via TimeSeries.hht()...")
eemd_result = ts_norm.hht(
    emd_method="eemd",
    emd_kwargs={
        "eemd_trials": 10,
        "eemd_noise_std": 0.2,
        "random_state": 42,
        "sift_max_iter": 200,
        "stopping_criterion": 0.2,
    },
    hilbert_kwargs={
        "pad": 200,
        "if_smooth": 11,
    },
    output="dict",
)

eimfs = eemd_result["imfs"]
n_eimfs = len(eimfs)
print(f"EEMD completed. Extracted {n_eimfs} IMFs.")

# Visualize all IMFs
fig, axes = plt.subplots(n_eimfs, 1, figsize=(10, 2 * n_eimfs), sharex=True)
if n_eimfs == 1:
    axes = [axes]

for i, key in enumerate(eimfs.keys()):
    imf_ts = eimfs[key]
    axes[i].plot(imf_ts.times.value, imf_ts.value * data_std, "k")
    axes[i].set_ylabel(f"IMF {i + 1}")
    axes[i].grid(alpha=0.3)
    if i == n_eimfs - 1:
        axes[i].set_xlabel("Time [s]")
        axes[i].set_xlim(0, duration)

plt.suptitle("EEMD Decomposed Modes")
plt.tight_layout()
plt.show()

Running EEMD on normalized data via TimeSeries.hht()...
EEMD completed. Extracted 9 IMFs.

5. Hilbert Spectrum (HHTSpectrogram)

By specifying output='spectrogram', you can obtain a time-frequency map. By default, HHTSpectrogram.plot() uses a logarithmic scale.

spec = ts_norm.hht(
    output="spectrogram",
    emd_method="eemd",
    emd_kwargs={
        "eemd_trials": 10,
        "eemd_noise_std": 0.2,
        "random_state": 42,
        "sift_max_iter": 200,
        "stopping_criterion": 0.2,
    },
    hilbert_kwargs={
        "pad": 200,
        "if_smooth": 11,
    },
    fmin=0,
    fmax=400,
    n_bins=120,
    weight="ia2",
)

plot = spec.plot()
plt.ylim(0, 400)
plt.show()

6. Low-level APIs (Optional)

For more fine-grained control, you can directly use the following methods:

• emd() for decomposition.

• hilbert_analysis() for IA/IF calculation.

• instantaneous_frequency() / envelope() for individual IMFs.

However, using hht() is recommended for general analysis.

Hints

• Edge effects are common in Hilbert analysis. Consider using hilbert_kwargs={'pad': N} to mitigate this, and ignore the very ends of the data.

• If you need smoother IMFs, increase eemd_trials (this will increase computation time).

• For strict reproducibility, use emd_method='emd'.

• Do not interpret negative or rapidly sign-flipping instantaneous frequency samples near amplitude nulls; mask low-amplitude regions before turning IF into physics.

Summary

Using TimeSeries.hht() allows you to practically apply the theoretical power of HHT:

• Tracking instantaneous frequencies of non-stationary signals like SASI (Section 2).

• Confirming superior resolution compared to STFT (Section 3).

• Resolving mode mixing and estimating physical quantities using EEMD (Section 4).

• Obtaining time-frequency maps with output='spectrogram' (Section 5).

In practice, the reliable workflow is: isolate the physics you care about, verify that the IMF/ridge is stable under modest analysis changes, and only then interpret the instantaneous frequency as a mode tracker.