visualdynamics.core.srs

srs

Shock response spectra.

An SRS is not a spectrum of the shock. It is the answer to a question asked of it: hang a single-degree-of-freedom oscillator off the base, shake the base with the measured transient, and record the largest response the oscillator ever reaches. Sweep the oscillator's natural frequency across the band of interest and the curve of those peaks is the shock response spectrum.

That is why two shocks with the same SRS can look nothing alike, and why an SRS cannot be inverted back to a time history: every peak is one number taken from one whole run of one filter, and the phase that produced it is gone.

The filter is Smallwood's ramp-invariant recursion (Smallwood 1981), which is the method the shock community settled on. Multiplying by the continuous transfer function in the frequency domain is the obvious alternative and it is wrong near Nyquist: the sampled input is a ramp between samples, not a train of impulses, and the ramp-invariant filter is the one that gets that right. The recursion below is exact for that assumption rather than an approximation of it.

Damping is quoted either as a ratio or as Q = 1/(2*zeta). Q = 10 — 5% of critical — is the near-universal default, and the value a specification is assumed to be written at unless it says otherwise.

Functions:

Name	Description
damping_for	The damping ratio an amplification factor means.
q_for	The amplification factor a damping ratio means.
octave_frequencies	Natural frequencies from low to high, per_octave to an octave.
ramp_invariant	Smallwood's filter coefficients, one row per natural frequency.
peaks	(highest, lowest) the oscillators reach, one per frequency.
maximax	The peak absolute response, whichever way it went.

Functions:

damping_for

damping_for(q: float = DEFAULT_Q) -> float

The damping ratio an amplification factor means.

Source code in src/visualdynamics/core/srs.py

def damping_for(q: float = DEFAULT_Q) -> float:
    """The damping ratio an amplification factor means."""
    return 1.0 / (2.0 * float(q))

q_for

q_for(damping: float) -> float

The amplification factor a damping ratio means.

Source code in src/visualdynamics/core/srs.py

def q_for(damping: float) -> float:
    """The amplification factor a damping ratio means."""
    return 1.0 / (2.0 * float(damping))

octave_frequencies

octave_frequencies(
    low: float, high: float, per_octave: int = PER_OCTAVE
) -> ndarray

Natural frequencies from low to high, per_octave to an octave.

Geometric, because an SRS is read on a log axis and a linear grid would crowd the top of the band and starve the bottom.

Source code in src/visualdynamics/core/srs.py

def octave_frequencies(low: float, high: float,
                       per_octave: int = PER_OCTAVE) -> np.ndarray:
    """Natural frequencies from `low` to `high`, `per_octave` to an octave.

    Geometric, because an SRS is read on a log axis and a linear grid
    would crowd the top of the band and starve the bottom.
    """
    low, high = float(low), float(high)
    if not (low > 0.0 and high > low):
        raise ValueError(f'octave band {low} to {high} is not a band')
    octaves = np.log2(high / low)
    return low * 2.0 ** (np.arange(round(octaves * per_octave) + 1)
                         / float(per_octave))

ramp_invariant

ramp_invariant(
    frequencies: ArrayLike,
    sample_rate: float,
    damping: float,
) -> tuple[ndarray, ndarray]

Smallwood's filter coefficients, one row per natural frequency.

Returns (b, a), each (F, 3), for the absolute acceleration of the oscillator given the acceleration of its base.

The gain at DC is exactly one — b.sum() / a.sum() == 1 — which is the statement that an infinitely stiff oscillator goes wherever its base goes. It is the cheapest check that the coefficients are the right ones, and test_srs makes it.

Source code in src/visualdynamics/core/srs.py

def ramp_invariant(frequencies: ArrayLike, sample_rate: float,
                   damping: float) -> tuple[np.ndarray, np.ndarray]:
    """Smallwood's filter coefficients, one row per natural frequency.

    Returns (b, a), each (F, 3), for the absolute acceleration of the
    oscillator given the acceleration of its base.

    The gain at DC is exactly one — b.sum() / a.sum() == 1 — which is
    the statement that an infinitely stiff oscillator goes wherever its
    base goes. It is the cheapest check that the coefficients are the
    right ones, and `test_srs` makes it.
    """
    frequencies = np.asarray(frequencies, dtype=float)
    dt = 1.0 / float(sample_rate)
    zeta = float(damping)
    if not 0.0 < zeta < 1.0:
        raise ValueError(f'damping {zeta} is not under critical')

    wn = 2.0 * np.pi * frequencies
    wd = wn * np.sqrt(1.0 - zeta ** 2)
    decay = np.exp(-zeta * wn * dt)
    turned = wd * dt
    cosine = decay * np.cos(turned)
    # sin(wd dt)/(wd dt) -> 1 as the frequency goes to nothing, which is
    # where the ramp-invariant form is otherwise 0/0
    shaped = decay * np.sinc(turned / np.pi)

    b = np.stack([1.0 - shaped,
                  2.0 * (shaped - cosine),
                  decay ** 2 - shaped], axis=-1)
    a = np.stack([np.ones_like(cosine),
                  -2.0 * cosine,
                  decay ** 2], axis=-1)
    return b, a

peaks

peaks(
    signal: ArrayLike,
    frequencies: ArrayLike,
    sample_rate: float,
    damping: float | None = None,
    q: float = DEFAULT_Q,
) -> tuple[ndarray, ndarray]

(highest, lowest) the oscillators reach, one per frequency.

The filter runs over the record once for every natural frequency at once — the state is a vector across frequencies and the loop is over time — so nothing of length (frequencies x samples) is ever held. Only the running extremes are kept, which is all an SRS is.

Source code in src/visualdynamics/core/srs.py

def peaks(signal: ArrayLike, frequencies: ArrayLike, sample_rate: float,
          damping: float | None = None,
          q: float = DEFAULT_Q) -> tuple[np.ndarray, np.ndarray]:
    """(highest, lowest) the oscillators reach, one per frequency.

    The filter runs over the record once for every natural frequency at
    once — the state is a vector across frequencies and the loop is over
    time — so nothing of length (frequencies x samples) is ever held.
    Only the running extremes are kept, which is all an SRS is.
    """
    x = np.asarray(signal, dtype=float).ravel()
    zeta = damping_for(q) if damping is None else float(damping)
    b, a = ramp_invariant(frequencies, sample_rate, zeta)
    b0, b1, b2 = b[:, 0], b[:, 1], b[:, 2]
    a1, a2 = a[:, 1], a[:, 2]

    y1 = np.zeros(b0.shape)
    y2 = np.zeros(b0.shape)
    highest = np.zeros(b0.shape)
    lowest = np.zeros(b0.shape)
    x1 = x2 = 0.0
    for sample in x:
        y = b0 * sample + b1 * x1 + b2 * x2 - a1 * y1 - a2 * y2
        np.maximum(highest, y, out=highest)
        np.minimum(lowest, y, out=lowest)
        y2, y1 = y1, y
        x2, x1 = x1, sample
    return highest, lowest

maximax

maximax(
    signal: ArrayLike,
    frequencies: ArrayLike,
    sample_rate: float,
    damping: float | None = None,
    q: float = DEFAULT_Q,
) -> ndarray

The peak absolute response, whichever way it went.

The reading a shock specification is written against: a shock that only ever pushes one way is no gentler for it.

Source code in src/visualdynamics/core/srs.py

def maximax(signal: ArrayLike, frequencies: ArrayLike, sample_rate: float,
            damping: float | None = None,
            q: float = DEFAULT_Q) -> np.ndarray:
    """The peak absolute response, whichever way it went.

    The reading a shock specification is written against: a shock that
    only ever pushes one way is no gentler for it.
    """
    highest, lowest = peaks(signal, frequencies, sample_rate, damping, q)
    return np.maximum(highest, -lowest)