Mode Decomposition with EMD, VMD, and SSA: Python Implementation and Comparison for Non-Stationary Signals

Why Mode Decomposition? Limits of Fourier and Wavelet

The Fourier transform assumes a superposition of stationary sinusoids and struggles with non-stationary signals whose frequency changes momentarily, or non-linear signals whose amplitude and frequency vary simultaneously. Short-Time Fourier Transform (STFT) and Wavelet analysis improve time-frequency localization, but they fix the basis (window or mother wavelet) up front, limiting accuracy when the basis does not match the data.

Mode decomposition extracts Intrinsic Mode Functions (IMFs) in a data-driven way from the signal itself. The three established methods are:

• EMD (Empirical Mode Decomposition) — Huang (1998); removes local means iteratively via the sifting algorithm.
• VMD (Variational Mode Decomposition) — Dragomiretskiy (2014); solves a constrained variational problem with ADMM.
• SSA (Singular Spectrum Analysis) — separates periodic and trend components via SVD of the trajectory matrix.

All three handle non-stationary, non-linear signals well, and combined with the Hilbert transform (https://yuhi-sa.github.io/en/posts/20260318_hilbert_transform/1/) form the Hilbert-Huang Transform (HHT) for instantaneous amplitude and frequency.

1. EMD: Empirical Mode Decomposition

1.1 IMF Conditions

An IMF must satisfy:

1. The number of extrema and the number of zero crossings differ by at most one.
2. The mean of the upper and lower envelopes (cubic-spline interpolation of maxima and minima) is zero at every point.

Locally this represents an “amplitude-modulated single-frequency component”.

1.2 The Sifting Algorithm

From \(x(t)\) , iterate

where \(m_k(t)\) is the mean of the upper and lower envelopes of \(h_k\) . Stop when a Cauchy SD criterion (typically 0.2–0.3) is satisfied, set \(\mathrm{IMF}_1 = h_K\) , and recurse on the residual \(r_1 = x - \mathrm{IMF}_1\) .

1.3 EMD in Python with PyEMD

import numpy as np
import matplotlib.pyplot as plt
from PyEMD import EMD

fs = 1000
t = np.arange(0, 2, 1 / fs)
x = (
    np.sin(2 * np.pi * (5 + 10 * t) * t)
    + (1 + 0.5 * np.sin(2 * np.pi * 1.0 * t)) * np.sin(2 * np.pi * 60 * t)
    + 0.1 * np.random.randn(len(t))
)

emd = EMD()
imfs = emd(x)
print(f"Number of IMFs: {len(imfs)}")

EMD suffers from mode mixing (similar frequencies appearing in different IMFs, or different frequencies sharing an IMF). EEMD (Ensemble EMD) mitigates this by averaging many trials with added white noise.

2. VMD: Variational Mode Decomposition

2.1 Constrained Variational Problem

VMD decomposes \(x(t)\) into \(K\) AM-FM modes \(u_k(t)\) designed so each mode has tight bandwidth around a center frequency \(\omega_k\) :

subject to \(\sum_k u_k(t) = x(t)\) .

2.2 ADMM Solution

Form the augmented Lagrangian and iterate via the Alternating Direction Method of Multipliers:

1. Update \(u_k\) as a Wiener-like filter in the frequency domain.
2. Update \(\omega_k\) as the centroid of \(|u_k|^2\) .
3. Dual ascent on the multiplier \(\lambda\) .

You must pre-specify \(K\) and the bandwidth penalty \(\alpha\) ; EMD is adaptive.

2.3 VMD in Python with vmdpy

from vmdpy import VMD

K = 4
alpha = 2000
tau = 0.0
DC = False
init = 1
tol = 1e-7

u, u_hat, omega = VMD(x, alpha, tau, K, DC, init, tol)
print(f"VMD center frequencies (Hz): {omega[-1] * fs}")

VMD suppresses EMD’s mode mixing and gives stable frequency separation. Choosing \(K\) via the elbow of energy contribution is common.

3. SSA: Singular Spectrum Analysis

3.1 Trajectory Matrix and SVD

Build the trajectory (Hankel) matrix with window length \(L\) :

SVD \(X = U \Sigma V^\top\) yields principal component modes from the largest singular values \(\sigma_i\) .

3.2 Grouping and Diagonal Averaging

Group SVD outputs into trend (large \(\sigma\) ), periodic (adjacent pairs), and noise (small \(\sigma\) ). Apply diagonal averaging (Hankelization) to recover each group as a time series.

3.3 SSA in Python (numpy only)

def ssa_decompose(x, L=100, n_components=4):
    N = len(x)
    K = N - L + 1
    X = np.array([x[i:i + L] for i in range(K)]).T
    U, s, Vt = np.linalg.svd(X, full_matrices=False)
    components = []
    for i in range(n_components):
        Xi = s[i] * np.outer(U[:, i], Vt[i, :])
        recon = np.zeros(N)
        counts = np.zeros(N)
        for ii in range(L):
            for jj in range(K):
                recon[ii + jj] += Xi[ii, jj]
                counts[ii + jj] += 1
        components.append(recon / counts)
    return np.array(components), s

modes, sigmas = ssa_decompose(x, L=200, n_components=4)

SSA has only two parameters (\(L\) and the number of components), is numerically stable thanks to linear algebra, and is often more robust than EMD/VMD on short or noisy series.

4. Comparison

5. Hilbert-Huang Transform: Instantaneous Frequency

Apply the Hilbert transform to each IMF \(c_k(t)\) to form the analytic signal \(z_k = c_k + j \mathcal{H}\{c_k\} = a_k(t) e^{j \phi_k(t)}\) . Instantaneous frequency is \(f_k(t) = (2\pi)^{-1} d\phi_k / dt\) . See https://yuhi-sa.github.io/en/posts/20260318_hilbert_transform/1/ for details.

from scipy.signal import hilbert

def hht(imfs, fs):
    analytic = hilbert(imfs, axis=-1)
    amp = np.abs(analytic)
    phase = np.unwrap(np.angle(analytic), axis=-1)
    freq = np.diff(phase, axis=-1) / (2 * np.pi) * fs
    return amp, freq

amp, freq = hht(imfs, fs)

HHT is used in bearing fault tracking, acoustic feature extraction, ECG R-wave analysis, and more.

6. Side-by-Side Comparison

from PyEMD import EMD as EMD_lib
from vmdpy import VMD

emd = EMD_lib()
imfs_emd = emd(x)
u_vmd, _, _ = VMD(x, alpha=2000, tau=0.0, K=4, DC=False, init=1, tol=1e-7)
ssa_modes, _ = ssa_decompose(x, L=200, n_components=4)

def dominant_freq(signal, fs):
    spec = np.abs(np.fft.rfft(signal))
    freqs = np.fft.rfftfreq(len(signal), 1 / fs)
    return freqs[np.argmax(spec)]

print("EMD dominant Hz:", [f"{dominant_freq(c, fs):.1f}" for c in imfs_emd[:4]])
print("VMD dominant Hz:", [f"{dominant_freq(c, fs):.1f}" for c in u_vmd])
print("SSA dominant Hz:", [f"{dominant_freq(c, fs):.1f}" for c in ssa_modes])

Practical guideline:

• Known mode count / stability → VMD
• Unknown / adaptive exploration → EMD (or EEMD)
• Short series / trend extraction → SSA

7. Applications

• Bearing diagnostics: extract modes near BPFO/BPFI fault frequencies, detect anomalies via instantaneous amplitude (VMD + HHT).
• ECG R-wave extraction: separate respiratory drift and pulse with SSA.
• Climate trend separation: yearly, seasonal, short-term variations into three SSA components.
• Acoustic scene analysis: VMD band-separates ambient and voice; feed into ML (see https://yuhi-sa.github.io/en/posts/20260525_ml_timeseries_guide/1/).

8. Learning Checklist

1. Two IMF conditions
2. Sifting stopping criteria (SD / S-number)
3. Difference between EEMD and CEEMDAN
4. Center-frequency update in VMD
5. Role of ADMM dual variable
6. Choosing \(L\) in SSA from data period
7. Implementing diagonal averaging
8. Relation between Hilbert transform and instantaneous frequency
9. Typical mode-mixing scenarios (close frequencies, intermittent signals)
10. Minimal code with PyEMD / vmdpy / numpy.linalg.svd

Related Articles

• https://yuhi-sa.github.io/en/posts/20260524_time_frequency_guide/1/ — Time-Frequency Analysis Selection Hub
• https://yuhi-sa.github.io/en/posts/20260318_hilbert_transform/1/ — Hilbert Transform and Analytic Signal
• https://yuhi-sa.github.io/en/posts/20260226_wavelet/1/ — Wavelet Transform
• https://yuhi-sa.github.io/en/posts/20260522_wavelet_packet/1/ — Wavelet Packet Decomposition
• https://yuhi-sa.github.io/en/posts/20260429_stft/1/ — Short-Time Fourier Transform
• https://yuhi-sa.github.io/en/posts/20260225_fft/1/ — Fast Fourier Transform
• https://yuhi-sa.github.io/en/posts/20260228_fft_window_psd/1/ — Window Functions and PSD
• https://yuhi-sa.github.io/en/posts/20260228_timeseries_anomaly/1/ — Time-Series Anomaly Detection
• https://yuhi-sa.github.io/en/posts/20260525_ml_timeseries_guide/1/ — ML Time-Series Hub