Variational Mode Decomposition

Decomposes a sampled signal into a given number of modes, each concentrated around a center frequency, by alternating updates of the mode spectra, the centers and a multiplier.

Usage

decomposeVMD(signal, samplingInterval, initialCenters, alpha, tau,
             fixedZeroCenter, tolerance, maxIterations)

Arguments

signal: Non-empty numeric vector of finite values without dim, in the units of the quantity.
samplingInterval: Sampling interval $\Delta t$ in seconds, a finite number greater than zero.
initialCenters: Numeric vector of initial center frequencies in Hz, each in $[0, 1/(2\Delta t)]$; its length is the number of modes $K$. No default.
alpha: Penalty $\alpha > 0$ of the bandwidth, in $s^2$ because it multiplies squared angular frequencies in rad/s. No default. sktime and VMDecomp take $\alpha$ on cycles per sample without the factor 2: their value $\alpha_\nu$ corresponds to $\alpha = \alpha_\nu \Delta t^2 / (8 \pi^2)$.
tau: Dual ascent step $\tau \ge 0$, dimensionless. 0 solves the penalized problem with a zero multiplier; a positive value the constrained problem. No default; the package tests use 0 and 0.5.
fixedZeroCenter: TRUE keeps the first center at 0 and requires initialCenters[1] == 0; FALSE updates every center. No default.
tolerance: Threshold $\epsilon > 0$ of the relative change of the mode spectra between sweeps, dimensionless. No default; the 1e-7 of sktime belongs to an absolute criterion, not to this relative change.
maxIterations: Whole number $\ge 1$ of sweeps. No default; sktime and VMDecomp allocate 500 states.

Details

The record of $N$ samples is extended to $T = 2N$ samples by reflecting $\lceil N/2 \rceil$ samples at the start and the rest at the end, and transformed with the unnormalized discrete Fourier transform; the grid holds $|\omega|$ in rad/s. Sweeps start from zero mode spectra and multiplier. In each sweep mode $k$, updated after modes $1, \ldots, k - 1$, is $$\hat u_k = \frac{\hat x - \sum_{i \ne k} \hat u_i + \hat\lambda / 2}{1 + 2 \alpha (\omega - \omega_k)^2},$$ its center becomes the power-weighted mean of $|\omega|$, and the sweep closes with $\hat\lambda \leftarrow \hat\lambda + \tau (\hat x - \sum_k \hat u_k)$. The iteration stops with "tolerance" when $\sum_k \|\hat u_k^{(n+1)} - \hat u_k^{(n)}\|^2 / \|\hat u_k^{(n)}\|^2 < \epsilon$, which cannot happen before the second sweep, or with "iterationLimit" after maxIterations sweeps; the last sweep is returned. The modes are the inverse transforms restricted to the record. With tau = 0 the penalized formulation permits a nonzero fit residual. VMD defines its modes through a bandwidth objective, without imposing the geometrical IMF conditions of EMD.

Value

A list with

modes: Double matrix, samples by modes, on the input grid.
fitResidual: Numeric vector $x - \sum_k u_k$.
centers: Numeric vector of the $K$ centers of the last sweep, in Hz.
iterations: Integer vector repeating, for every mode, the number of sweeps.
converged: Logical vector repeating, for every mode, the result of the joint stopping test.
stopReason: "tolerance" or "iterationLimit".

Conditions

A mode without spectral energy leaves its center undefined and stops the function with an error of class xemdZeroEnergy.

References

Dragomiretskiy, K. and Zosso, D. (2014). Variational Mode Decomposition. IEEE Transactions on Signal Processing, 62(3), 531–544. doi:10.1109/TSP.2013.2288675

Examples

SamplingInterval <- 1 / 1000
Time <- seq_len(1000) * SamplingInterval
Signal <- cos(2 * pi * 2 * Time) + 0.25 * cos(2 * pi * 24 * Time) +
    0.0625 * cos(2 * pi * 288 * Time)
OUT <- decomposeVMD(Signal, samplingInterval = SamplingInterval,
                    initialCenters = (0:2) / (6 * SamplingInterval),
                    alpha = 2000 * SamplingInterval^2 / (8 * pi^2),
                    tau = 0.5, fixedZeroCenter = FALSE, tolerance = 1e-6,
                    maxIterations = 98)
OUT$centers
#> [1]   2.000004  23.999532 287.962079
OUT$stopReason
#> [1] "tolerance"
max(abs(OUT$fitResidual))
#> [1] 0.02030723