Empirical Mode Decomposition

Decomposes a signal by successive extraction: each component is obtained by sifting the current residue and is subtracted from that residue. Its numerical acceptance depends on the selected stopping rule.

Usage

decomposeEMD(signal, sifting)

Arguments

signal: Numeric vector of finite real values without dim, in the units of the quantity. Positions are sample indices, so no time vector or sampling interval is needed. Attributes such as names, tsp or class do not reach the outputs.
sifting: Named list of sifting controls: stopRule and the fields of that rule, described in section ‘Sifting controls’. Fields with a default may be omitted; every other field is required, and a field that does not belong to the rule is an error.

Details

With $r_0 = x$, extraction $k$ computes the mode $c_k = E_1(r_{k-1})$ by sifting and updates the residue $$r_k = r_{k-1} - c_k,$$ so that $x = \sum_{k=1}^{K} c_k + r_K$ up to rounding. Before extraction $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, and with "insufficientExtrema" when the boundary rule cannot build the envelopes of $r_{k-1}$.

Sifting starts from $h_0 = r_{k-1}$ and repeats $h_{j+1} = h_j - m_j$, where $m_j$ is the mean of the upper and lower envelopes of $h_j$. The envelopes interpolate the maxima and the minima, extended at both borders by sifting$boundary: a cubic spline with not-a-knot end conditions through four or more nodes, the parabola through three and the line through two. A run of equal interior samples is one extremum at the exact centre of the run. Each state $h_j$ is examined once, in this order: the stopping rule accepts it (converged TRUE); $j$ equals fixedSift under "fixed" (converged NA); $j$ equals maxSift under an adaptive rule, or the envelopes of $h_j$ cannot be built (converged FALSE); otherwise one more subtraction follows. Non-convergence is reported in converged, without a warning.

Sifting controls

The same list configures the extractor of decomposeEEMD, decomposeCEEMD, decomposeCEEMDAN, decomposeICEEMDAN, resiftModes, decomposeMaskedEMD, decomposeIteratedMaskedEMD, estimateMaskFrequency and, with the restrictions stated in their help pages, decomposeBEMD, decomposeMEMD and decomposeNAMEMD. Counts are accepted as doubles with a whole value.

stopRule: Categorical, without default: "SD", "Snumber", "doubleThreshold" or "fixed". The rule "consecutive" belongs to the vector sifting of decomposeMEMD and is rejected here.
tolerance: Under "SD": dimensionless threshold, greater than zero, without default. From $j = 1$ the state is accepted when $\sum_n (h_{j-1}[n] - h_j[n])^2 / h_{j-1}[n]^2$ is finite and below it; the sum runs over all samples, without the factor $1/N$, and a null denominator is not regularized. The 1998 formulation proposes values from 0.2 to 0.3.
stableSift: Under "Snumber": the count $S \ge 1$, without default. From $j = 1$ the state is accepted when the counts of extrema $E$ and zero crossings $Z$ have stayed equal over $S$ consecutive states and $|E - Z| \le 1$. The 2003 formulation reports typical values from three to five.
thresholds: Under "doubleThreshold": two dimensionless numbers $0 < \theta_1 < \theta_2$; default c(0.05, 0.5).
fraction: Under "doubleThreshold": the admitted fraction $\alpha$ in $[0, 1)$; default 0.05. With $\sigma = |m / a|$, $m$ the half sum and $a$ the half difference of the envelopes, the state, the signal included, is accepted when $|E - Z| \le 1$, the fraction of samples with $\sigma \ge \theta_1$ is at most $\alpha$ and no sample has $\sigma \ge \theta_2$. Where $m = a = 0$ the ratio is undefined and counts as no exceedance; where only $a = 0$ it is infinite.
fixedSift: Under "fixed": the count $J \ge 1$ of subtractions of every extraction, without default.
maxSift: Under "SD", "Snumber" and "doubleThreshold": the count $M \ge 1$ of subtractions after which an extraction ends without acceptance; under "Snumber" $M \ge S$. Default 2000 under "doubleThreshold", none under the other rules. The rule "fixed" has no maxSift.
maxModes: Count $\ge 1$ or Inf, the limit of extractions; default Inf, which extracts until the residue ends the decomposition.
boundary: Categorical: "reflection" or "extrapolation"; default "reflection" under "doubleThreshold" and "fixed", none under "SD" and "Snumber". "reflection" mirrors two extrema of each family per side about the first extremum or the border sample and needs three extrema; "extrapolation" anchors the upper envelope at each border at the larger of the border sample and the line through the two nearest maxima, and the lower envelope at the smaller of the border sample and the line through the two nearest minima.

Value

A list with

modes: Double matrix, samples by modes, in the units of signal; it has zero columns when no mode is extracted. No names.
residue: Numeric vector $r_K$ of the input length; equal to signal when no mode is extracted.
iterations: Integer vector with the number of subtractions of each mode.
converged: Logical vector by mode: TRUE when the stopping rule accepted the mode, FALSE when maxSift was reached or the envelopes could not be built, NA under "fixed".
stopReason: "terminalResidue", "modeLimit" or "insufficientExtrema".

References

Huang, N. E., Shen, Z., Long, S. R., Wu, M. C., Shih, H. H., Zheng, Q., Yen, N.-C., Tung, C. C. and Liu, H. H. (1998). The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis. Proceedings of the Royal Society of London. Series A, 454, 903–995. doi:10.1098/rspa.1998.0193

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 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 = Inf)
OUT <- decomposeEMD(Signal, sifting = Sifting)
dim(OUT$modes)
#> [1] 3000    3
OUT$iterations
#> [1] 2 4 2
OUT$stopReason
#> [1] "insufficientExtrema"
max(abs(rowSums(OUT$modes) + OUT$residue - Signal))
#> [1] 1.332268e-15

# The rule alone: the other fields take their defaults.
identical(decomposeEMD(Signal, sifting = list(stopRule = "doubleThreshold")),
          OUT)
#> [1] TRUE

# Ten fixed subtractions per mode, at most three modes.
Fixed <- decomposeEMD(Signal, sifting = list(stopRule = "fixed",
                                             fixedSift = 10, maxModes = 3))
Fixed$converged
#> [1] NA NA NA
Fixed$stopReason
#> [1] "modeLimit"