Signal processing: AT2TS / VT2TS / DT2TS

In gmsp: Ground Motion Signal Processing

knitr::opts_chunk$set(collapse = TRUE, comment = "#>", eval = FALSE)

# 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.

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.

gmsp documentation built on July 18, 2026, 5:07 p.m.