Masked Empirical Mode Decomposition

Decomposes a signal by successive extraction in which each mode is the average first mode of the residue plus sinusoidal masks of several phases, under one of three mask rules: the opposite sine pair of 2005, dyadic masks (EMD-DMS) and uniformly distributed phases (UPEMD).

Usage

decomposeMaskedEMD(signal, samplingInterval, masking, sifting)

Arguments

signal

Numeric vector of finite real values without dim, in the units of the quantity.

samplingInterval

Sampling interval in seconds, a finite number greater than zero, without default. Mask frequencies are in Hz.

masking

Named list with the field rule and exactly the fields of that rule, without defaults:

"pair": frequency, in Hz, in $(0, 1/(2\Delta t)]$, for instance the value of estimateMaskFrequency; amplitude, greater than zero, in the units of the signal.
"dyadic": firstFrequency, in Hz, in $(0, 1/(2\Delta t)]$; amplitudeRatio, the dimensionless $V_m > 0$. The example of the EMD-DMS article uses $V_m = 10$.
"uniformPhase": amplitudeRatio, the dimensionless $\varepsilon_0 > 0$; phaseCount, the whole number $P \ge 2$ of phases. The package tests use $\varepsilon_0 = 1$ and $P = 4$.

sifting

Named list of sifting controls, as in decomposeEMD; the rule "consecutive" is rejected.

Details

From $r_0 = x$, before stage $k$ the decomposition ends, in this order, with "terminalResidue" when $r_{k-1}$ has at most one extremum, with "modeLimit" when $k - 1$ equals sifting$maxModes or the schedule of masks is exhausted, and with "insufficientExtrema" when the envelopes of $r_{k-1}$ cannot be built. Otherwise the masks of stage $k$ are $a_k \cos(2 \pi f_k n \Delta t + \theta_p)$, $n = 0, \ldots, N - 1$, with the first sample as phase origin, and the mode is the average of the first sifted component of $r_{k-1}$ plus each mask, $$c_k = \frac{1}{P} \sum_{p=1}^{P} E_1\left(r_{k-1} + w_p\right), \qquad r_k = r_{k-1} - c_k.$$ A stage is retained unless one of its phases has no first mode or loses its envelopes during sifting; then the decomposition stops before adding the stage mean or subtracting it from the residue. A phase that reaches the iteration limit is kept and flags the stage FALSE. The phases and divisor are unchanged. The rules fix the phases, frequencies and amplitudes:

"pair": Phases $-\pi/2$ and $\pi/2$, so the masks are $\pm a \sin(2 \pi f n \Delta t)$ and the mode is $(E_1(x + w) + E_1(x - w)) / 2$. The schedule has one stage: the decomposition ends with "modeLimit" after one extraction unless the residue is terminal.
"dyadic": Phases 0 and $\pi$; $f_k = 2^{1-k} f_1$; amplitude $V_m$ times the mean absolute value of the maxima of the residue.
"uniformPhase": Phases $2 \pi p / P$, $p = 0, \ldots, P - 1$; $f_k = 1 / (2^k \Delta t)$; amplitude $\varepsilon_0$ times the standard deviation of the residue with divisor $N$.

The masks are deterministic: the function draws no random numbers.

Value

A list with

modes: Double matrix, samples by modes.
residue: Numeric vector, the final residue.
iterations: Integer vector by mode, the sum of the subtractions of the extractions of the stage.
converged: Logical vector by mode: TRUE when every extraction of the stage accepted; NA when all completed fixed sifting. Interrupted stages are reported in termination.
stopReason: "terminalResidue", "modeLimit", "insufficientExtrema" or "incompleteExtraction", the last when at least one phase fails.
maskFrequencies: Numeric vector, the mask frequency of each mode in Hz.
termination: Data frame with mode, the next stage index after the retained modes; phase, its one-based index in the phase sequence; iterations, subtractions of that failed phase; and stopReason, "terminalResidue", "modeLimit", "insufficientExtrema" or "iterationLimit". One row is returned for each phase of the interrupted stage that did not accept, including the phases that reached the iteration limit. If the decomposition stops before a stage starts, one row carries its stopReason, phase = NA and zero iterations.

References

Deering, R. and Kaiser, J. F. (2005). The use of a masking signal to improve empirical mode decomposition. In IEEE International Conference on Acoustics, Speech, and Signal Processing (ICASSP), 4, 485–488.

The selection of the reflected extrema and the spline interpolation of the envelope nodes are adapted from emd.m (G. Rilling, March 2007), distributed with CEEMDAN revision af9a0f83, copyright (c) 2020 Marcelo A. Colominas, MIT License; see the NOTICE file of the installed package.

Examples

Time <- seq(0, 9, length.out = 3000)
Signal <- sin(pi * Time) + sin(2 * pi * Time) + sin(6 * pi * Time) +
    0.5 * Time
Sifting <- list(stopRule = "doubleThreshold", thresholds = c(0.05, 0.5),
                fraction = 0.05, maxSift = 2000, boundary = "reflection",
                maxModes = 3)
# Uniform phases, with one sample per second.
OUT <- decomposeMaskedEMD(Signal, samplingInterval = 1,
                          masking = list(rule = "uniformPhase",
                                         amplitudeRatio = 1, phaseCount = 4),
                          sifting = Sifting)
OUT$maskFrequencies
#> [1] 0.500 0.250 0.125
OUT$stopReason
#> [1] "modeLimit"

# The pair of 2005 at the estimated mask frequency.
Frequency <- estimateMaskFrequency(Signal, samplingInterval = 1,
                                   sifting = Sifting)
Pair <- decomposeMaskedEMD(Signal, samplingInterval = 1,
                           masking = list(rule = "pair", frequency = Frequency,
                                          amplitude = sqrt(mean((Signal -
                                              mean(Signal))^2))),
                           sifting = Sifting)
dim(Pair$modes)
#> [1] 3000    1
Pair$stopReason
#> [1] "modeLimit"