Successive Variational Mode Decomposition

Extracts variational modes one at a time from the whole signal, keeping the centers of the modes already extracted fixed, and stops after a known number of modes or when the residual power meets a known noise power.

Usage

decomposeSVMD(signal, samplingInterval, initialCenter = 0, alpha, maxAlpha,
              tau, tolerance, maxIterations = 299, stopping)

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.

initialCenter

Initial center of every mode in Hz, in $[0, 1/(2\Delta t)]$; default 0.

alpha

Initial penalty $\alpha > 0$ in $s^2$. No default.

maxAlpha

Final penalty, at least alpha, in $s^2$; equal to alpha gives a fixed penalty. No default. PySDKit uses 10 and 20000 for its denominator $1 + 2 \alpha_\nu (\nu - \nu_L)^2$ on cycles per sample $\nu$; a value $\alpha_\nu$ there corresponds to $\alpha_\nu \Delta t^2 / (4 \pi^2)$.

tau

Dual ascent step $\tau \ge 0$. No default; PySDKit uses 0.

tolerance

Threshold $\epsilon_1 > 0$ of the relative change of the mode spectrum between inner iterations. No default; PySDKit uses 1e-6.

maxIterations

Whole number $\ge 1$ of inner iterations for each value of the penalty; default 299. The defaults of initialCenter and maxIterations are those of the authors' svmd.m 1.1.1.

stopping

Named list with the field stopRule and exactly the fields of that rule:

"modeCount": modeCount, the whole number of modes to extract, without default.
"noisePower": noisePower, the noise power $\sigma^2 > 0$ in squared units of the signal; noiseTolerance, the dimensionless $\epsilon_2 > 0$; maxModes, a whole number $\ge 1$ or Inf, default Inf.

Details

The signal is extended and transformed as in decomposeVMD. Mode $L$ is extracted from the whole signal with the centers $\omega_1, \ldots, \omega_{L-1}$ of the previous modes fixed. Its penalty starts at alpha and doubles while below maxAlpha; when doubling would exceed maxAlpha, one last computation runs at maxAlpha. Each penalty restarts the inner iteration from the spectrum and center reached, with a zero multiplier and the penalty $H_{L-1}(\omega) = \sum_{i<L} 1 / (\alpha^2 (\omega - \omega_i)^4)$. With $A = \alpha^2 (\omega - \omega_L)^4$, the inner update is $$\hat u_L \leftarrow \frac{\hat x + A \hat u_L + \hat\lambda / 2}{(1 + A)\left(1 + 2 \alpha (\omega - \omega_L)^2 + H_{L-1}\right)},$$ followed by the new center and the dual ascent, and it stops when the relative change of $\hat u_L$ is below tolerance or after maxIterations. Extraction stops with "modeCount" after modeCount modes; with "noisePower" when $|\sigma^2 - P_L| / \sigma^2 < \epsilon_2$, $P_L$ the mean square of the residual on the record; or with "modeLimit" after maxModes modes without meeting it. With maxModes = Inf the extraction under "noisePower" continues until the condition holds. The modes keep the extraction order. Their bandwidth formulation does not impose the geometrical IMF conditions of EMD.

Value

A list with

modes: Double matrix, samples by modes, in extraction order.
fitResidual: Numeric vector $x - \sum_L u_L$.
centers: Numeric vector of the centers in Hz.
iterations: Integer vector by mode, the inner iterations of all its penalties.
converged: Logical vector by mode, TRUE when every inner loop of the mode met the change test.
stopReason: "modeCount", "noisePower" or "modeLimit".

Conditions

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

References

Nazari, M. and Sakhaei, S. M. (2020). Successive variational mode decomposition. Signal Processing, 107610. doi:10.1016/j.sigpro.2020.107610

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 <- decomposeSVMD(Signal, samplingInterval = SamplingInterval,
                     initialCenter = 0,
                     alpha = 10 * SamplingInterval^2 / (4 * pi^2),
                     maxAlpha = 20000 * SamplingInterval^2 / (4 * pi^2),
                     tau = 0, tolerance = 1e-6, maxIterations = 299,
                     stopping = list(stopRule = "modeCount", modeCount = 5))
OUT$centers
#> [1]   1.999983 287.997040 278.041353  23.999702 297.789812
OUT$stopReason
#> [1] "modeCount"