---
title: "IMF decomposition: TS2IMF (EMD / EEMD / VMD)"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{IMF decomposition: TS2IMF (EMD / EEMD / VMD)}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(collapse = TRUE, comment = "#>", eval = FALSE)
```

```{r minimal-example, eval = TRUE}
# 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)
```

`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:

1. the number of zero crossings and extrema differ at most by 1, and
2. 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:

```r
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:

1. apply a Hann window $w[n]$,
2. compute FFT with zero-padding to a power of 2,
3. compute a power spectrum $P(f)$,
4. 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:

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

```r
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"]
```

## References

* Huang, N. E., Shen, Z., Long, S. R., Wu, M. C., Shih, H. H., Zheng,
  Q., Yen, N.-C., Tung, C. C., & 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 A*, 454(1971), 903–995.
* Wu, Z., & Huang, N. E. (2009). Ensemble Empirical Mode
  Decomposition: A Noise-Assisted Data Analysis Method. *Advances in
  Adaptive Data Analysis*, 1(1), 1–41.
* Dragomiretskiy, K., & Zosso, D. (2014). Variational Mode
  Decomposition. *IEEE Transactions on Signal Processing*, 62(3),
  531–544.