Signal processing: AT2TS / VT2TS / DT2TS

# 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)
#>        t            s     ID   OCID
#>    <num>        <num> <char> <char>
#> 1: 0.000 0.000000e+00     AT     H1
#> 2: 0.008 1.003931e-05     AT     H1
#> 3: 0.016 1.802513e-01     AT     H1
#> 4: 0.024 2.961337e-01     AT     H1
#> 5: 0.032 3.912661e-01     AT     H1
#> 6: 0.040 4.817434e-01     AT     H1

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

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:

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.

tsw <- TSL2TSW(tsl)
roundtrip <- TSW2TSL(tsw)
head(roundtrip)
#>      OCID     ID     t            s
#>    <char> <char> <num>        <num>
#> 1:     H1     AT 0.000 0.000000e+00
#> 2:     H1     AT 0.008 1.003931e-05
#> 3:     H1     AT 0.016 1.802513e-01
#> 4:     H1     AT 0.024 2.961337e-01
#> 5:     H1     AT 0.032 3.912661e-01
#> 6:     H1     AT 0.040 4.817434e-01

Time regularisation

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

\[\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:

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:

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:

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)

References

mirror server hosted at Truenetwork, Russian Federation.