---
title: "Signal processing: AT2TS / VT2TS / DT2TS"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Signal processing: AT2TS / VT2TS / DT2TS}
  %\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)
dt_acc <- data.table(t = seq(0, 5, by = 0.01), H1 = sin(2 * pi * 2 * seq(0, 5, by = 0.01)))
tsl <- AT2TS(dt_acc, units.source = "mm", Fmax = 10, audit = FALSE, output = "TSL")
head(tsl)
```

`gmsp` implements three workflows to build a consistent set of
**Acceleration (AT)**, **Velocity (VT)** and **Displacement (DT)** time
series from a single input quantity:

* `AT2TS()` starts from acceleration; produces VT and DT by integration.
* `VT2TS()` starts from velocity; produces AT by differentiation and DT
  by integration.
* `DT2TS()` starts from displacement; produces VT and AT by
  differentiation.

All three share the same internal pipeline:

1. Time regularisation if needed.
2. STFT parameter selection and resampling decision
   (internal strategy + optional `auditSTFT()`).
3. Optional anti-alias resampling.
4. Edge tapering to "turn off" low-amplitude pre/post segments.
5. Integration or differentiation (STFT-domain or finite differences).

References for STFT/windows and OLA: Allen & Rabiner (1977),
Harris (1978).

## Data formats

### Input (wide)

The workflows expect a `data.table` with:

* a time column (default `t`), and
* one or more signal columns (e.g. `H1`, `H2`, `UP`).

Use `time = "..."` when the input time column is not named `t`. The
pipeline canonicalizes time internally and `TSL` output is always `t`
and `s`; `COL.t` and `COL.s` are not part of this interface.

### Outputs

| `output` | Shape | Columns |
|---|---|---|
| `"ATo"`, `"VTo"`, `"DTo"` | wide | `ts` (time from 0), `Units`, channels |
| `"AT"`, `"VT"`, `"DT"` | wide | channels only |
| `"TSW"` | wide | `ts`, `AT.<OCID>`, `VT.<OCID>`, `DT.<OCID>` |
| `"TSL"` (default) | long | `t`, `s`, `ID` ∈ {AT,VT,DT}, `OCID` |

Canonical table projections are also available after construction. Use
`TSL2TSW()` to cast a long table to wide columns and `TSW2TSL()` to return to
the canonical long contract. `TSL2TSW()` writes canonical time as `t`;
`TSW2TSL()` accepts either `t` or constructor-style `ts`.

```{r tsl-tsw-converters, eval = TRUE}
tsw <- TSL2TSW(tsl)
roundtrip <- TSW2TSL(tsw)
head(roundtrip)
```

## Time regularisation

Given a time series with potentially irregular sampling, it builds a
uniform grid:

* Let $t_0 = \min(t)$ and $t_1 = \max(t)$.
* Let $\Delta t = \min(t_i - t_{i-1})$.
* "Rationalize" $\Delta t$ to force an integer sampling rate $F_s$:

$$\Delta t \leftarrow \mathrm{rationalize}(\Delta t),
  \qquad F_s = 1/\Delta t \in \mathbb{Z}.$$

Then it interpolates each channel using a monotone cubic Hermite
spline (`splinefun(..., method = "monoH.FC")`).

**Practical motivation:** many later steps (STFT bin quantization,
resampling to `TargetFs`) assume integer $F_s$.

Reference for monotone cubic interpolation: Fritsch & Carlson (1980).

Regularisation is an internal step of `AT2TS()`, `VT2TS()`, `DT2TS()`,
`TS2IMF()`, and `TSL2PS()`. It always works on canonical time `t`;
callers select a different input time column through the constructor
argument `time`.

## Internal STFT Strategy Selection

When `Fmax` is provided, the internal STFT strategy tries to choose `Fs_out` and `NW`
(from a candidate grid) to:

* ensure enough STFT windows (`MW >= MW.min`), and
* **minimise leakage** around `Fmax` by aligning `Fmax` to an STFT bin.

The leakage criterion used is

$$J = \left|\frac{F_{\max}}{\Delta f}
       - \mathrm{round}\!\left(\frac{F_{\max}}{\Delta f}\right)\right|
       + \text{penalty}(\Delta F_s),
  \qquad \Delta f = F_s / NW.$$

If the user passes `kNyq` explicitly, it forces
$F_{s,\mathrm{target}} \approx \mathrm{round}(kNyq\,F_{\max})$.

The function returns a list with `NW` (window length), `OVLP`
(overlap), `MW` (window count), `Fs` (target sampling rate after
resampling), the flags `Resample` / `Upsample` / `Downsample`, and a
`strategy` label.

## STFT audit: `auditSTFT()`

`auditSTFT()` is invoked from `AT2TS/VT2TS/DT2TS` when `audit = TRUE`.
It performs two checks.

**Spectrum and content above Fmax.** Estimates the fraction of
spectral amplitude above `Fmax` (via `buildFFT()`), and warns if the
spectral peak `Fpeak` is above `Fmax`.

**Strategy consistency.** Re-evaluates the anti-alias heuristics
(transition bins `bins.tr`, margin to Nyquist `nyq.margin`, window
count `MW`).

Returns `list(warnings = ..., info = ..., stft = ...)`.

## STFT / OLA filtering: `.ffilter()`

`.ffilter()` applies a frequency-domain filter defined by a `custom`
vector in the STFT domain:

1. Compute STFT with `seewave::stdft()` using a Hann window
   (`wn = "hanning"`).
2. Multiply each frequency-bin row by `custom`.
3. Reconstruct with `seewave::istft()`.

Notation: let $Z[k, m]$ be the STFT (bin $k$, frame $m$) and $H[k]$ be
the filter (`custom`). Then

$$\tilde{Z}[k,m] = H[k]\,Z[k,m]$$

and ISTFT/OLA yields $x[n]$.

## Integration: `.integrate()`

Integrates a signal $x[n]$ using a frequency-domain (STFT-frame-wise)
filter. The complex integration filter is

$$H_I(\omega_k) = \begin{cases}
  0, & k = 0 \\
  \dfrac{1}{i\,\omega_k}, & k > 0
\end{cases}$$

applied via `.ffilter(..., custom = H_I)`. The DC component is set to
zero (the integration constant is "fixed"), which helps reduce drift.

Implementation details:

* RMS normalisation
  $x \leftarrow x / \max(\mathrm{RMS}(x), \epsilon)$, then rescale.
* Padding to an STFT-compatible length via `.padZeros()`.

## Differentiation: `.derivate()`

Two methods are available via `method`.

**`method = "time"` (finite differences).** Four-point central
difference in the interior,

$$\dot{x}_i \approx
  \frac{-x_{i+2} + 8x_{i+1} - 8x_{i-1} + x_{i-2}}{12\,\Delta t}$$

with three-point formulas at the edges.

**`method = "freq"` (STFT-domain).** The derivative filter

$$H_D(\omega_k) = \begin{cases}
  0, & k = 0 \\
  i\,\omega_k, & k > 0
\end{cases}$$

optionally multiplied by an additional low-pass if `lowPass = TRUE` and
`Fmax` is provided.

## Anti-alias resampling: `.resample()`

Per channel:

1. Design a Butterworth-like low-pass with `Fpass`, `Fstop` (quantised
   to bins).
2. Filter via `.ffilter()` (anti-alias).
3. Resample to `TargetFs` using `signal::resample`.
4. Remove DC.
5. Rescale amplitudes to preserve peak (per-channel peak matching).

The Butterworth-like low-pass magnitude implemented is

$$\left|H_{LP}(f)\right|
   = \frac{1}
          {\sqrt{1
             + \left(\dfrac{1}{A_{stop}^2}-1\right)
               \left(\dfrac{f}{F_{stop}}\right)^{\!2N}}}$$

where the order $N$ is chosen to approximately satisfy
$|H(F_{pass})| = A_{pass}$.

## Edge tapering: `.taperA()`

`.taperA()` builds a window $w[n] \in [0,1]$ based on thresholds
normalised by the max amplitude:

* `Astop` and `Apass` are compared against $|x|/\max|x|$.
* It finds start/end indices where amplitude exceeds those thresholds.
* It creates smooth transitions using Butterworth-like filters in
  index-space (not frequency).

`AT2TS/VT2TS/DT2TS` apply this window at the beginning and the end of
each channel, and can optionally trim with `trimZeros`.

A companion `.taperI()` exists (energy-cumulative taper) but is not
used by the current workflows.

## Notes (parameter audit)

* `resample` is a parameter on `AT2TS/VT2TS/DT2TS`, but the effective
  decision is made by the internal STFT strategy (`STFTParams$Resample`).
* `time` selects the input time column for `AT2TS/VT2TS/DT2TS`; `TSL`
  output uses canonical `t` and `s`.

## References

* Allen, J. B., & Rabiner, L. R. (1977). A unified approach to
  short-time Fourier analysis and synthesis. *Proceedings of the
  IEEE*, 65(11), 1558–1564.
* Harris, F. J. (1978). On the use of windows for harmonic analysis
  with the discrete Fourier transform. *Proceedings of the IEEE*,
  66(1), 51–83.
* Fritsch, F. N., & Carlson, R. E. (1980). Monotone Piecewise Cubic
  Interpolation. *SIAM Journal on Numerical Analysis*, 17(2), 238–246.
* Oppenheim, A. V., & Schafer, R. W. (2010). *Discrete-Time Signal
  Processing* (3rd ed.). Pearson.