HoloHilbert.jl — Holo-Hilbert Spectral Analysis

Why you need this: Fourier analysis cannot detect multiplicative interactions in signals. When a slow modulation multiplies a carrier (not adds), classical spectrograms show only noise. HHSA reveals the hidden modulation.

The Problem: Fourier Fails on Multiplicative Signals

Consider a hybrid signal that combines a multiplicative component (noise whose amplitude is modulated) with an additive one (a carrier with AM), both driven by the same modulation frequency $\Omega = 0.3$ Hz:

$$x(t) = \underbrace{[1 + \sin(2\pi\Omega t)],\varepsilon(t)}_{\text{multiplicative}} + \underbrace{[1 + 0.7\sin(2\pi\Omega t)]\sin(2\pi\omega_c t)}_{\text{additive, }\omega_c = 5\text{ Hz}}$$

Fourier cannot report $\Omega$ as a spectral line. For the additive part the identity

$$[1 + m\sin\Omega t]\sin\omega t = \sin\omega t + \tfrac{m}{2}\left[\cos(\omega-\Omega)t - \cos(\omega+\Omega)t\right]$$

puts the modulation into sidebands at $\omega \pm \Omega$, never at $\Omega$ itself; for the multiplicative part it smears across the whole noise floor. HHSA instead applies a second Hilbert transform to the EMD amplitude envelope (the "holo" step), which measures $\Omega$ directly:

Figure: Panel (b): the Fourier spectrum shows the 5 Hz carrier and its 4.7 / 5.3 Hz sidebands, but the noise floor at the marked $\Omega = 0.3$ Hz is flat — the modulation frequency itself never appears. Panel (c): the holo-spectrum resolves the energy in (carrier $\omega$ × modulation $\Omega$) space. Note the triangular support: the physical constraint $\omega > \Omega$ leaves the upper-left half-plane empty.

What is HHSA?

Holo-Hilbert Spectral Analysis extends the classical Hilbert-Huang Transform with a second layer:

1. Stage 1 (Standard HHT):

   • Decompose signal into Intrinsic Mode Functions (IMFs) via Empirical Mode Decomposition.
   • Extract instantaneous amplitude $A_j(t)$ and carrier frequency $\omega_j(t)$ via Hilbert transform.

2. Stage 2 (The "Holo" Step):

   • Decompose each amplitude envelope $A_j(t)$ again via EMD.
   • Extract modulation frequency $\Omega_{j,k}(t)$ from the second-layer IMFs.

3. Result:

   • A 2-D holo-spectrum showing energy density in (carrier frequency ω × modulation frequency Ω) space, not just (time × frequency).

Mathematical Summary

Each Stage-1 IMF $c_j(t)$ is converted to an analytic signal via the Hilbert transform:

$$z_j(t) = c_j(t) + i \mathcal{H}[c_j(t)] = A_j(t) e^{i\phi_j(t)}$$

The amplitude envelope $A_j(t)$ then undergoes EMD to reveal modulation:

$$A_j(t) = \sum_k IA_{j,k}(t) + r_j(t)$$

where $IA_{j,k}$ are modulation IMFs. Applying the Hilbert transform again yields $\Omega_{j,k}(t)$.

The holo-spectrum is the 2-D energy map:

$$H(\omega, \Omega) = \sum_{j,k} \int A_{j,k}^2(t) , \delta(\omega - \omega_j(t)) , \delta(\Omega - \Omega_{j,k}(t)) , dt$$

For a detailed treatment, see Huang et al. (2016) "On Holo-Hilbert spectral analysis: a full informational spectral representation for nonlinear and non-stationary data," Phil. Trans. R. Soc. A 374, 20150196. (doi:10.1098/rsta.2015.0196)

Features

• Pure Julia EMD — natural cubic spline envelopes, plateau-aware extrema detection (clipped or quantised data would otherwise lose over half its extrema), mirror boundary extension, and the Rilling et al. (2003) stopping criterion by default (stop_criterion = :rilling). Unlike Huang's SD criterion, it tests the IMF condition before subtracting the mean envelope, so a mono-component AM signal is not over-sifted. This matters: over-sifting erodes the amplitude modulation that the second HHSA layer exists to measure. The legacy SD criterion is available as stop_criterion = :sd.
• Minimal dependencies — only DSP.jl (Hilbert transform) and Plots.jl (visualization). No legacy packages, no SPM12.
• Two-layer HHSA pipeline — automatically extracts both carrier and modulation frequencies.
• Masked sifting (mask = :auto) — the mode-mixing cure of Deering & Kaiser (2005), applied in both layers as the reference HHSA toolbox does. Without it a carrier that shares an IMF with another component has a wandering ω(t), and its modulation energy smears along the ω axis until the feature disappears. On the nested cross-scale signal of Quinn et al. (2021) it is the difference between recovering the second modulation and missing it entirely: without a mask the five strongest density peaks are all the same 5 Hz feature smeared along ω, while with one the two ground-truth blobs come out as peaks #1 and #2, at (ω, Ω) = (4.72, 0.51) and (36.2, 4.57) Hz against a truth of (5, 0.5) and (37, 5).
• Physically constrained — the paper's $\omega > \Omega$ condition is enforced by default, and non-oscillatory envelope components (the DC level of $A_j(t)$) are excluded from the spectrum instead of swamping it at $\Omega \approx 0$.
• Reproducible — 25 passing tests, each annotated with the paper section it checks. On the clean AM benchmark $(1+0.8\sin 2\pi 6t)\sin 2\pi 60t$ the holo-spectrum peak lands at $(\omega, \Omega) = (59.4, 5.9)$ Hz against a ground truth of $(60, 6)$.

Installation

Requirements

• Julia 1.10 or later (tested with 1.12)

Quick Setup

Once registered in Julia General Registry:

using Pkg
Pkg.add("HoloHilbert")

Or from the repo:

cd /path/to/HoloHilbert
julia --project=.

Then from the Julia REPL:

using Pkg
Pkg.instantiate()   # Install DSP.jl and Plots.jl

Usage

Interactive Notebook Demo

For an interactive walkthrough of all three benchmark cases (additive AM, multiplicative noise, two-carrier AM):

cd notebooks
jupyter notebook hhsa_paper_demo.ipynb

Basic Example

using HoloHilbert

# Your signal: 1D Float64 vector, sampling rate fs in Hz
x = your_signal_vector  
fs = 1000.0              # Hz

# Run HHSA (Holo-Hilbert Spectral Analysis)
result = hhsa(x, fs; max_imfs=8, max_sift=200)

# Visualize
plot_imfs(result)                    # Stage-1 IMF decomposition
plot_hilbert_spectrum(result)        # Time × carrier frequency
plot_holo_spectrum(result)           # Carrier frequency × modulation frequency

Run the Demo

julia --project=. examples/demo.jl

This generates synthetic AM signals and produces four PNGs in output/:

• 01_signal.png — raw input
• 02_stage1_imfs.png — EMD decomposition
• 03_hilbert_spectrum.png — classical time-frequency map
• 04_holo_spectrum.png — the holo-spectrum (ω × Ω)

The Multiplicative / Additive Test

julia --project=. examples/paper_test_multiplicative.jl

This builds the hybrid signal shown at the top of this README (multiplicative noise + additive AM carrier, both at $\Omega = 0.3$ Hz) and produces the three-panel figure. It follows the argument of Section 2 of Huang et al. (2016) rather than reproducing their Figure 2b numerically — the paper's own example uses $f_s = 1$ Hz with a 0.01 Hz sine, which is exercised in the test suite instead (test/runtests.jl, "Paper core demonstration").

Main Functions

emd(signal; max_imfs=10, max_sift=200, stop_criterion=:rilling, θ1=0.05, θ2=0.5, α=0.05, sd_threshold=0.1)

Empirical Mode Decomposition. Returns a vector of IMFs (last element is the residue).

• max_imfs: Maximum number of IMFs to extract.
• max_sift: Maximum iterations per sifting loop.
• stop_criterion: :rilling (default) or :sd.
• θ1, θ2, α: Rilling thresholds — sifting stops when $\sigma(t) = |m(t)/a(t)| < \theta_2$ everywhere and $< \theta_1$ on a fraction $1-\alpha$ of samples.
• sd_threshold: Huang SD threshold; used only when stop_criterion = :sd.

instantaneous_params(signal, fs) → (amplitude, frequency)

Instantaneous amplitude and frequency from the Hilbert analytic signal. The frequency is signed: a negative value flags a sample where the input is not a well-formed IMF. Taking abs of it would fold that onto a high positive frequency and fabricate energy, so consumers discard those samples instead.

hhsa(signal, fs; envelope=:hilbert, mask=nothing, emd_kw...)

Two-layer Holo-Hilbert Spectral Analysis. Remaining keywords pass through to emd.

• mask: nothing (default), :auto, or a Stage-1 mask frequency in cycles per sample. Enables masked sifting in both layers (Stage 2 always auto-estimates its own, since an envelope's frequency content is unrelated to the Stage-1 mask). Note the trade-off: on a signal that plain sifting already resolves cleanly, masking can fragment a mode across several low-energy IMFs. Use it when mode mixing is the problem, not by reflex.
• envelope: how the Stage-1 amplitude $A_j(t)$ is obtained — :hilbert (default, $|z_j(t)|$) or :nht (iterative normalisation, Huang et al. 2009). A 12-point sweep over $f_c \in {40,60,80}$ Hz × $f_m \in {2,4,6,10}$ Hz found :hilbert better on both accuracy (mean $|\Omega$ error$|$ 0.146 vs 0.178 Hz) and spurious energy (0.17% vs 0.37%), so :nht is offered for robustness checks rather than as an improvement.

Returns a HoloResult struct containing:

• imfs, carrier_amp, carrier_freq — Stage 1 results
• mod_imfs, mod_amp, mod_freq — Stage 2 results
• fs — sampling rate

plot_imfs(result)

Stacked plot of Stage-1 IMFs and residue.

plot_hilbert_spectrum(result; n_freq=200, f_max=fs/2)

Time × carrier frequency heatmap (classical Hilbert spectrum).

plot_holo_spectrum(result; n_carrier=200, n_mod=100, enforce_constraint=true, dc_tol=0.5, ...)

Carrier frequency × modulation frequency heatmap (the holo-spectrum).

holo_spectrum(result; n_carrier=128, n_mod=64, f_carrier_max=fs/2, f_mod_max=fs/4, enforce_constraint=true, dc_tol=0.5, density=false, spacing=:linear, f_carrier_min=1.0, f_mod_min=0.01)

Compute the binned 2-D energy matrix directly (without plotting).

• enforce_constraint: discard $(\omega, \Omega)$ pairs with $\omega \le \Omega$ (paper Section 3). Set false only for signals with no genuine carrier, where $\omega(t)$ is not a meaningful frequency.

• dc_tol: a Stage-2 mode enters the spectrum only if $|\text{mean}| < \texttt{dc_tol} \cdot \text{rms}$, i.e. it satisfies the zero-local-mean IMF condition. This excludes the DC level of the envelope, which otherwise carries more energy than the modulation itself and pins the peak at $\Omega \approx 0$.

• density: divide by bin area and sample interval so that $\iint H , d\omega , d\Omega$ equals the total energy at any grid resolution. Note this makes the integral grid-invariant, not the peak: a near-monochromatic feature concentrates into fewer bins as the grid refines, so its peak density grows — correct behaviour for a density.

• spacing: :linear (default) or :log. Log axes match the convention of the Python emd package, but they are not an improvement here: a 128-bin log axis over 1–100 Hz is wider at 37 Hz (1.34 Hz) and much narrower at 5 Hz (0.18 Hz) than the linear grid, so it biases the density toward slow carriers. On the test signal above the (37, 5) feature falls from the 2nd-ranked bin to the 72nd. Use it to compare against published figures drawn on log axes, and pair it with density=true, without which bins of wildly different width are not comparable at all.

Energy outside the grid is discarded, not clamped onto the edge bins.

All entry points (emd, hhsa, instantaneous_params) accept any real vector — ranges, views, Float32, Int — and convert once at the boundary.

Repository Structure

Source code:

• src/HHSA.jl — module entry point
• src/emd.jl — EMD sifting loop, spline envelopes, extrema detection
• src/holo_spectrum.jl — Hilbert transform, two-layer HHSA pipeline, visualization

Examples & demos:

• examples/demo.jl — synthetic AM signal (automatic figure generation)
• examples/paper_test_multiplicative.jl — reproduce README Figure (hybrid signal)
• examples/parameter_search.jl — sweep over carrier/modulation frequencies

Notebooks:

• notebooks/hhsa_paper_demo.ipynb — interactive paper benchmarks (Jupyter)

Tests & metadata:

• test/runtests.jl — 25 unit tests with section references
• Project.toml, Manifest.toml — dependency management
• LICENSE — MIT
• References/ — papers and documents

Notes & Caveats

• Performance: This implementation prioritizes clarity and transparency over optimization. For very large signals (>100k samples), consider parallelizing the inner EMD loop or using compiled spline libraries.
• Parameter tuning: EMD is sensitive to sd_threshold and max_sift. Lower sd_threshold extracts finer modes but increases computation. Adjust for your data.
• End effects: The code reflects the outermost extremum past each boundary to reduce end artifacts. This is simpler than the multi-extremum extension of Rilling et al. (2003); for periodic signals, consider wrapping or extending the signal before analysis.
• Over-sifting: Damage to amplitude modulation is monotonic in the number of sifts. On a mono-component AM signal, a single superfluous sift injects ~16 spurious extrema into the Hilbert envelope and splits the modulation across three Stage-2 modes. This is why :rilling is the default; if you switch to :sd, expect the modulation peak to be biased low (≈4.7 Hz instead of 6 Hz on the benchmark above).
• Bin resolution: The holo-spectrum is computed on a fixed 2-D grid. Adjust n_carrier, n_mod, f_carrier_max, f_mod_max to match your frequency scales of interest.

Reference

Original paper & foundational works:

• Huang, N. E., Shen, Z., Long, S. R., Wu, M. C., Shih, H. H., Zheng, Q., Yen, N.-C., Tung, C. C., & Liu, H. H. (1998). The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis. Proc. R. Soc. Lond. A 454(1971), 903–935.

• Rilling, G., Flandrin, P., & Gonçalves, P. (2003). On empirical mode decomposition and its algorithms. IEEE-EURASIP Workshop on Nonlinear Signal and Image Processing (NSIP-03).

• Huang, N. E., Hu, K., Yang, A. C. C., Chang, H.-C., Jia, D., Liang, W.-K., Yeh, J.-R., Kao, C.-L., & Liu, C. (2016). On Holo-Hilbert spectral analysis: a full informational spectral representation for nonlinear and non-stationary data. Phil. Trans. R. Soc. A 374, 20150196. (link)

Citation

If you use this package in your research, please cite the original authors of the algorithm (Huang et al.).

Feedback

For questions, bug reports, or feature requests, please open an issue on GitHub or contact the maintainer at akejja@estudiantes.fisica.unam.mx.