The hardware and bandwidth for this mirror is donated by dogado GmbH, the Webhosting and Full Service-Cloud Provider. Check out our Wordpress Tutorial.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]dogado.de.

Getting started with gmsp

library(gmsp)
library(data.table)

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

These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.
Health stats visible at Monitor.