IMF decomposition: TS2IMF (EMD / EEMD / VMD)

# Minimal executable example
library(gmsp)
library(data.table)
t_vec <- seq(0, 5, by = 0.01)
dt_sig <- data.table(t = t_vec, s = sin(2 * pi * 2 * t_vec) + 0.3 * sin(2 * pi * 8 * t_vec))
imf_res <- TS2IMF(dt_sig, method = "vmd", K = 4, output = "IMF")
print(imf_res)
#>       IMF    Fo[Hz]     BW[Hz] E[k]/Eo [%]
#>    <char>     <num>      <num>       <num>
#> 1:   IMF1 0.8121344 0.48023152    1.479176
#> 2:   IMF2 1.9518203 0.11210220   43.727802
#> 3:   IMF3 2.0687453 0.08820797   27.534890
#> 4:   IMF4 7.9992331 0.10765251    7.920870

TS2IMF() decomposes a time series into Intrinsic Mode Functions (IMFs) and a residue. A common goal is to separate oscillatory components by scale / band and then reconstruct a filtered signal by removing or keeping selected IMFs.

Canonical references: EMD (Huang et al., 1998), EEMD (Wu & Huang, 2009), VMD (Dragomiretskiy & Zosso, 2014).

What is an IMF?

In the EMD framework, a component $c(t)$ is an IMF if it satisfies the classic definition:

the number of zero crossings and extrema differ at most by 1, and
the local mean of the upper and lower envelopes is approximately zero.

This enables a meaningful instantaneous frequency (via the Hilbert transform) for non-stationary signals.

Input and grouping

.x: one signal as a data.table with canonical columns t and s.

Note: internally TS2IMF() regularizes the time grid and rebuilds a time vector ts from 0. The workflow requires t / s as the working column names.

TS2IMF() is a worker and does not detect groups. For a canonical TSL, group outside the helper:

TSL[ID == "AT",
    TS2IMF(.SD, method = "vmd", K = 6, output = "TSL"),
    by = .(RecordID, OCID),
    .SDcols = c("t", "s")]

Available engines

Parameter method selects one of three engines.

VMD (`method = "vmd"`, default)

Calls VMDecomp::vmd(...). Returns a matrix $U \in \mathbb{R}^{N \times K}$ ($K$ modes). The residue is defined as

\[r[n] = s[n] - \sum_{k=1}^{K} U_k[n].\]

EMD (`method = "emd"`)

Calls EMD::emd(...) with boundary and stop.rule controlling sifting. The output AUX$imf and AUX$residue are taken as-is.

EEMD (`method = "eemd"`)

Calls hht::EEMD(...) followed by hht::EEMDCompile(...). Runs an ensemble with noise (noise.amp, noise.type, trials). Uses a tempdir() as a working folder and then compiles the averaged IMFs.

Output

If output = NULL, TS2IMF() returns a list with three elements.

Element	Shape	Description
`TSW`	wide	`t`, `signal`, `IMF1..IMFk`, `residue`
`TSL`	long	melt of `TSW` with columns `t`, `IMF`, `s`
`IMF`	metrics	per-IMF center frequency / bandwidth / relative energy

If output is "TSW", "TSL" or "IMF", only that component is returned.

Per-IMF metrics

For each IMF the code estimates the center frequency Fo[Hz], the bandwidth BW[Hz], and the relative energy E[k]/Eo [%]. The procedure is:

apply a Hann window $w[n]$,
compute FFT with zero-padding to a power of 2,
compute a power spectrum $P(f)$,
apply a −10 dB threshold relative to the maximum and compute Fo / BW as moments of $P(f)$ over “active” bins.

Reconstruction filters

Inside TS2IMF() a local helper .applyIMFFilter() adds a reconstructed signal signal.R when any of the following arguments is set.

imf.remove — if character, removes those columns (IMF1, IMF2, …); if numeric, positive values remove IMFs counted from the first mode, negative values remove IMFs counted from the last mode, and both selections are unioned. For example, c(1L, -1L, -2L) removes IMF1 plus the last two IMFs. Zero, non-finite, and out-of-range numeric values are ignored.
remove.Fo = c(f1, f2) — removes IMFs whose frequency interval $[F_o - BW,\,F_o + BW]$ lies entirely within $[f_1, f_2]$.
keep.Fo = c(f1, f2) — keeps only IMFs whose frequency interval $[F_o - BW,\,F_o + BW]$ lies entirely within $[f_1, f_2]$.
keep.Residue — if TRUE, signal.R adds the residue on top of the kept IMFs.

Composing with AT2TS / VT2TS / DT2TS

In practice the typical pattern is:

Start from a TSL (e.g. output of AT2TS(..., output = "TSL")).
Filter a specific ID + OCID (e.g. ID == "DT" & OCID == "H1").
Run TS2IMF() on that slice, or inside a data.table grouped call.
Use RES$TSL[IMF == "signal.R"] as the reconstructed (filtered) signal.

AT_H1 <- TSL[ID == "AT" & OCID == "H1", .(t, s)]
RES   <- TS2IMF(AT_H1, method = "vmd", K = 6,
                keep.Fo = c(2, 15), keep.Residue = TRUE)

# Reconstructed signal (only IMFs whose [Fo-BW, Fo+BW] fits in [2, 15] Hz):
RES$TSL[IMF == "signal.R"]