Getting started with gmsp

This vignette walks through the signal-processing core of gmsp on a synthetic acceleration record. The four entry points exercised are:

Synthetic input

We build a 5-second, two-channel acceleration record sampled at 200 Hz, with two narrowband harmonic components per channel and a mild amplitude envelope. Amplitudes are in mm/s² (the package’s canonical base unit).

NP <- 1000L
dt <- 1 / 200
t  <- seq.int(0, by = dt, length.out = NP)

env <- exp(-((t - 2.5) / 1.5)^2)

DT <- data.table(
  t  = t,
  H1 = env * (3.0 * sin(2 * pi *  5 * t) + 0.6 * sin(2 * pi * 20 * t)),
  H2 = env * (2.0 * sin(2 * pi *  7 * t) + 0.4 * sin(2 * pi * 17 * t))
)

str(DT)
#> Classes 'data.table' and 'data.frame':   1000 obs. of  3 variables:
#>  $ t : num  0 0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04 0.045 ...
#>  $ H1: num  0 0.0517 0.0952 0.1242 0.1375 ...
#>  $ H2: num  0 0.0402 0.0764 0.1045 0.1221 ...
#>  - attr(*, ".internal.selfref")=<pointer: 0x103679470>

Run AT2TS

AT2TS() takes wide input (one time column plus one or more signal columns). With output = "TSL" (the default), it returns a long table keyed by (t, s, ID, OCID) where ID ∈ {AT, VT, DT} is the kinematic quantity and OCID is the channel identifier inherited from the input columns.

TSL <- AT2TS(
  DT,
  units.source = "mm",
  units.target = "mm",
  Fmax        = 25,
  output      = "TSL",
  audit       = FALSE,
  verbose     = FALSE
)

TSL[, .N, by = .(ID, OCID)]
#>        ID   OCID     N
#>    <char> <char> <int>
#> 1:     AT     H1   998
#> 2:     AT     H2   998
#> 3:     VT     H1   998
#> 4:     VT     H2   998
#> 5:     DT     H1   998
#> 6:     DT     H2   998
head(TSL[ID == "VT" & OCID == "H1"], 6)
#>        t     s     ID   OCID
#>    <num> <num> <char> <char>
#> 1: 0.000     0     VT     H1
#> 2: 0.005     0     VT     H1
#> 3: 0.010     0     VT     H1
#> 4: 0.015     0     VT     H1
#> 5: 0.020     0     VT     H1
#> 6: 0.025     0     VT     H1

A quick visual sanity check on the H1 channel:

op <- par(mfrow = c(3, 1), mar = c(4, 4, 2, 1))
plot(TSL[ID == "AT" & OCID == "H1", .(t, s)], type = "l",
     main = "AT (mm/s^2)", xlab = "", ylab = "")
plot(TSL[ID == "VT" & OCID == "H1", .(t, s)], type = "l",
     main = "VT (mm/s)",   xlab = "", ylab = "")
plot(TSL[ID == "DT" & OCID == "H1", .(t, s)], type = "l",
     main = "DT (mm)",     xlab = "t [s]", ylab = "")

AT, VT and DT for channel H1

par(op)

Intensity measures

getIntensity() consumes a long TSL with columns RSN, OCID, ID, t, s. It computes 20 + scalar measures per (RSN, OCID, ID) group: PGA / PGV / PGD, RMS, zero crossings, Arias intensity and its positive / negative variants, three significant-duration measures (D5–95, D5–75, D20–80), CAV and CAV5, mean period Tm, and the derived indices EPI and PDI.

TSL[, RSN := "demo"]
IM <- getIntensity(TSL, units.source = "mm", units.target = "mm")

dcast(IM[IM %in% c("PGA", "PGV", "PGD", "AI", "D0595", "CAV")],
      OCID + IM ~ ., value.var = "value")
#> Key: <OCID, IM>
#>       OCID     IM            .
#>     <char> <char>        <num>
#>  1:     H1     AI 0.0014181966
#>  2:     H1    CAV 4.9525029281
#>  3:     H1  D0595 2.4750000000
#>  4:     H1    PGA 3.4213696194
#>  5:     H1    PGD 0.0032915973
#>  6:     H1    PGV 0.1035611225
#>  7:     H2     AI 0.0006313606
#>  8:     H2    CAV 3.3362851682
#>  9:     H2  D0595 2.4600000000
#> 10:     H2    PGA 2.3665485261
#> 11:     H2    PGD 0.0011334207
#> 12:     H2    PGV 0.0501314469

Response spectrum

TSL2PS() integrates the SDOF equation of motion for a canonical TSL table and reports PSA / PSV / SD. It derives grouping keys from TSL metadata instead of requiring a public BY argument.

Tn  <- 10 ^ seq(log10(0.05), log10(3), length.out = 60)

PS <- TSL2PS(TSL[OCID == "H1"], xi = 0.05, Tn = Tn, output = "PSL")

head(PS)
#>       RSN   OCID         Tn     ID        S
#>    <char> <char>      <num> <char>    <num>
#> 1:   demo     H1 0.00000000    PSA 3.421370
#> 2:   demo     H1 0.05000000    PSA 8.831002
#> 3:   demo     H1 0.05359301    PSA 6.370405
#> 4:   demo     H1 0.05744422    PSA 4.882943
#> 5:   demo     H1 0.06157217    PSA 4.297133
#> 6:   demo     H1 0.06599676    PSA 4.016459

plot(PS[ID == "PSA", .(Tn, S)], log = "x",
     type = "l", lwd = 2,
     xlab = "Tn [s]", ylab = "PSA [mm/s^2]",
     main = "5%-damped pseudo-acceleration spectrum (H1)")
#> Warning in xy.coords(x, y, xlabel, ylabel, log): 1 x value <= 0 omitted from
#> logarithmic plot

Pseudo-acceleration response spectrum

IMF decomposition

TS2IMF() decomposes one signal into intrinsic mode functions (EMD / EEMD / VMD) and a residue, optionally returning a reconstruction filtered by frequency content. Default engine is VMD.

AT_H1 <- TSL[ID == "AT" & OCID == "H1", .(t, s)]
IMFs  <- TS2IMF(AT_H1, method = "vmd", K = 6, output = "TSW")

names(IMFs)
#> [1] "t"       "signal"  "IMF1"    "IMF2"    "IMF3"    "IMF4"    "IMF5"   
#> [8] "IMF6"    "residue"
head(IMFs[, 1:5])
#>        t signal         IMF1       IMF2      IMF3
#>    <num>  <num>        <num>      <num>     <num>
#> 1: 0.000      0 -0.003560409 -0.3520249 0.2301532
#> 2: 0.004      0 -0.003556392 -0.3458191 0.2259957
#> 3: 0.008      0 -0.003548982 -0.3336690 0.2181172
#> 4: 0.012      0 -0.003538805 -0.3157821 0.2066475
#> 5: 0.016      0 -0.003526167 -0.2924627 0.1917765
#> 6: 0.020      0 -0.003511544 -0.2641076 0.1737500

op <- par(mfrow = c(2, 1), mar = c(4, 4, 2, 1))
plot(IMFs$t, IMFs$IMF1, type = "l",
     main = "IMF1", xlab = "", ylab = "")
plot(IMFs$t, IMFs$IMF2, type = "l",
     main = "IMF2", xlab = "t [s]", ylab = "")

First two IMFs of the AT H1 signal

par(op)

Where to go next

Four vignettes ship with the package:

vignette("signal-processing",  package = "gmsp")  # AT2TS / VT2TS / DT2TS math
vignette("imfs",               package = "gmsp")  # TS2IMF decomposition (EMD / EEMD / VMD)
vignette("spectra",            package = "gmsp")  # TSL2PS elastic SDOF spectra
vignette("intensity-measures", package = "gmsp")  # getIntensity output details
vignette("database",           package = "gmsp")  # optional file-based indexing layer

In-session help: ?AT2TS, ?TS2IMF, ?TSL2PS, ?getIntensity, or help(package = "gmsp").

Rendered documentation: https://averriK.github.io/gmsp/.