visualdynamics.core.shock_synthesis

shock_synthesis

A shock a shaker can play, synthesized to meet a shock response spectrum: a sum of wavelets, their amplitudes iterated until every point of the spectrum is inside the requirement's band, then read against the shaker's shock ratings.

An SRS requirement leaves the waveform free. Many time histories have one spectrum, and which one is played decides the event's length, its peak acceleration and the travel it asks of the shaker. So a synthesis is two questions at once: does the waveform's spectrum meet the requirement, and can the shaker play the waveform. Both are answered here, with the package's own readings: the spectrum by core.srs and judged by compliance.srs_tolerance, the travel, velocity and force by core.stroke.read_stroke against a core.shaker.ShakerRating.

The wavelet. Each component is a tone of frequency f windowed by one half cycle of a slower sine, N half cycles of the tone long, with N odd and at least three:

w(t) = A sin(πt/T) sin(2πft),   0 ≤ t ≤ T = N / (2f),

and zero outside. It starts and ends at zero acceleration, and both its integrals come back to zero at its end: written as a difference of two cosines of (N - 1)π t/T and (N + 1)π t/T, the velocity is a difference of their sines, zero at T for any whole N, and the displacement a difference of one minus their cosines, zero at T when N ± 1 are even, which is N odd. So every wavelet leaves the table where it found it, at rest and at center, and so does any sum of them: no compensation pulse is needed, and none is added. On the samples the same holds exactly when the wavelet spans a whole number of samples, so each wavelet's length is rounded to whole samples and its tone set from that length, N times the sample rate over twice the samples, a fraction of a percent off the point it is for at the top of the band. Its samples then sum to zero (the velocity's trapezoid rule ends at zero) and are symmetric about the middle (so the displacement's ends at zero too).

The iteration. There is one wavelet per point of the requirement, at that point's frequency, all of one polarity and centered on one instant. Each pass computes the spectrum of the sum at the requirement's points, of the requirement's type and at its Q, and scales each wavelet's amplitude by the requirement over the spectrum at its point. The closest pass whose every point lies inside the requirement's inner band (its warning limits, read as compliance.srs_tolerance reads them, or compliance.SRS_INNER_DB without them) is kept, and the passes go on while they bring the spectrum closer, until it is within ACCURACY_DB everywhere or a pass gains less than STALL_DB, so the record has room inside the band for what the shaker adds to it. Wavelets a twelfth of an octave apart overlap in frequency: in opposite polarity they cancel at the oscillators between them, and the iteration chases the cancellation rather than closing on the requirement, so one polarity is kept. A requirement no pass meets in MAX_PASSES is refused, naming its worst point.

The shaker. The record is read raw, from rest at its first sample, by read_stroke, against the ratings given. A wavelet of more half cycles builds its oscillator's response over more cycles, so it meets the same spectrum with a smaller amplitude and asks for less velocity and travel, and a longer event. The synthesis starts at the fewest half cycles asked for and, while the shaker cannot hold the result, tries two more at a time up to MOST_HALF_CYCLES, keeping the first that fits: the shortest event the shaker can play. When none fits, the requirement is refused, naming the limit that runs out at the closest attempt and the band of the requirement whose wavelets alone take that limit: the low end for the travel and the velocity, which are spent at low frequency, the high end for the force. The spectrum is never split into a low and a high pulse to fit, since the article's modes are then loaded separately, and a test that splits it is a different test.

References

The wavelet, its integrals and the iteration are worked out above; this is where the rules the synthesis keeps are set out.

1. MIL-STD-810H (2019), Environmental Engineering Considerations and Laboratory Tests, Method 516.8, Shock. Paragraph 2.3.3.1: no synthesized transient is unique to a spectrum, a pulse that exceeds the exciter calls for a compromise, and splitting a spectrum into a low- and a high-frequency pulse is not acceptable. Paragraph 4.2.2.2: the spectrum is judged point by point against a band on a share of the points. Annex A, paragraph 1.1: a record whose spectrum is computed is sampled at ten times its highest frequency or more.

Classes:

Name	Description
SynthesizedChannel	One channel's synthesis: its wavelets, how it was reached, and
ShockSynthesis	A synthesized shock and what it was synthesized to.

Functions:

Name	Description
wavelet	One wavelet of unit amplitude on its own samples, both ends
synthesize	A shock a shaker can play, whose spectrum meets specification.
suggested_rate	A sample rate for a synthesis: SUGGESTED_RATIO times the

Classes

SynthesizedChannel dataclass

SynthesizedChannel(dof: str, points: ndarray, target: ndarray, frequency: ndarray, samples: ndarray, start: ndarray, amplitude: ndarray, half_cycles: int, passes: int, tried: tuple[tuple[int, str, float], ...], spectrum: ndarray, reading: StrokeReading)

One channel's synthesis: its wavelets, how it was reached, and what it asks of the shaker.

The wavelets are one per point of the requirement, in its order: each one's tone as played (frequency, Hz), its length and its start in samples, and its amplitude in m/s^2.

Attributes:

Name	Type	Description
worst_db	float	The spectrum's largest distance from the requirement, in dB

Attributes

worst_db property

worst_db: float

The spectrum's largest distance from the requirement, in dB either way.

ShockSynthesis dataclass

ShockSynthesis(history: TimeHistory, specification: ShockSpecification, rating: ShakerRating, channels: list[SynthesizedChannel] = list())

A synthesized shock and what it was synthesized to.

history is the record, one acceleration channel per DOF of the requirement synthesized, in SI, from rest to rest with a quiet tail of a period of the lowest point so the spectrum's low oscillators reach their peaks inside it. channels holds each one's wavelets and its reading against the shaker.

Methods:

Name	Description
tolerance	The spectrum judged against the requirement, every channel,
describe	The synthesis in lines, a channel at a time: its wavelets and

Attributes:

Name	Type	Description
srs	Srs	The record's spectrum at the requirement's points, of its

Attributes

srs property

srs: Srs

The record's spectrum at the requirement's points, of its type and at its Q, one record per channel; NaN at a point outside the band a channel's requirement is written on.

Methods:

tolerance

tolerance() -> list[tuple[str, str | None, dict[str, Any]]]

The spectrum judged against the requirement, every channel, as compliance.srs_tolerance_all judges a measured one: the requirement's own bands and the default share rule.

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

def tolerance(self) -> list[tuple[str, str | None, dict[str, Any]]]:
    """The spectrum judged against the requirement, every channel,
    as `compliance.srs_tolerance_all` judges a measured one: the
    requirement's own bands and the default share rule."""
    from .compliance import srs_tolerance_all

    return srs_tolerance_all(self.specification, self.srs)

describe

describe(unit_system: Any = None) -> list[str]

The synthesis in lines, a channel at a time: its wavelets and passes, its event length, the band it met, then its stroke reading's judgments (StrokeReading.describe).

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

def describe(self, unit_system: Any = None) -> list[str]:
    """The synthesis in lines, a channel at a time: its wavelets and
    passes, its event length, the band it met, then its stroke
    reading's judgments (`StrokeReading.describe`)."""
    from ..units import DEFAULT_SYSTEM

    us = DEFAULT_SYSTEM if unit_system is None else unit_system
    judged = {dof: result for dof, _event, result in self.tolerance()}
    lines = []
    for c in self.channels:
        rate = self.history.sample_rate
        length = float(c.samples.max()) / rate
        lines.append(f'{c.dof}: {c.points.size} wavelets of '
                     f'{c.half_cycles} half cycles, {c.passes} '
                     f'pass{"es" * (c.passes != 1)}, event '
                     f'{1e3 * length:.4g} ms')
        result = judged.get(c.dof)
        if result is not None:
            lines.append(f'{result["inside_count"]}/{result["points"]} '
                         'points inside the band, within '
                         f'{c.worst_db:.2f} dB of the requirement')
        lines.extend(line for line in c.reading.describe(us)
                     if ' against ' in line)
    return lines

Functions:

wavelet

wavelet(samples: int, half_cycles: int) -> ndarray

One wavelet of unit amplitude on its own samples, both ends included: sin(πn/M) sin(πNn/M) for n from 0 to M, M = samples, N = half_cycles.

Its tone is N half cycles over M samples, N times the sample rate over 2M in Hz. Its samples sum to zero and are symmetric about the middle when N is odd, so integrated twice by the trapezoid rule its velocity and its displacement both end at zero.

Raises:

Type	Description
ValueError	If half_cycles is not an odd whole number of at least three, or the wavelet is shorter than two samples a half cycle.

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

def wavelet(samples: int, half_cycles: int) -> np.ndarray:
    """One wavelet of unit amplitude on its own samples, both ends
    included: sin(πn/M) sin(πNn/M) for n from 0 to M, M = `samples`,
    N = `half_cycles`.

    Its tone is N half cycles over M samples, N times the sample rate
    over 2M in Hz. Its samples sum to zero and are symmetric about the
    middle when N is odd, so integrated twice by the trapezoid rule its
    velocity and its displacement both end at zero.

    Raises
    ------
    ValueError
        If `half_cycles` is not an odd whole number of at least three,
        or the wavelet is shorter than two samples a half cycle.
    """
    n = _half_cycles(half_cycles)
    m = int(samples)
    if m < 2 * n:
        raise ValueError(f'a wavelet of {n} half cycles over {m} samples '
                         'is under two samples a half cycle')
    k = np.arange(m + 1)
    return np.sin(np.pi * k / m) * np.sin(np.pi * n * k / m)

synthesize

synthesize(specification: ShockSpecification, *, sample_rate: float, rating: ShakerRating | None = None, half_cycles: int = HALF_CYCLES, dofs: Sequence[str] | None = None, progress: Callable[[int, int], None] | None = None) -> ShockSynthesis

A shock a shaker can play, whose spectrum meets specification.

made = synthesize(spec, sample_rate=20480.0,
                  rating=ShakerRating(stroke=0.0508, velocity=1.78,
                                      force=35600.0,
                                      moving_mass=45.0))
made.history, made.tolerance(), made.channels[0].reading

One wavelet per point of the requirement, their amplitudes iterated until every point of the record's spectrum lies inside the requirement's inner band; then the record is read against rating raw (core.stroke.read_stroke). While the shaker cannot hold it the wavelets are lengthened two half cycles at a time, up to MOST_HALF_CYCLES, and the first that fits is kept. The module docstring has the method.

Parameters:

Name	Type	Description	Default
specification	ShockSpecification	The requirement, an acceleration SRS in declared units, with its Q and type; its warning limits are the band every point is put inside, compliance.SRS_INNER_DB without them.	required
sample_rate	float	Samples per second, above RATE_RATIO times the highest point.	required
rating	ShakerRating	The shaker's shock ratings and the masses it carries; a rating left out is not judged, and with none the first synthesis that meets the requirement is kept.	None
half_cycles	int	The fewest half cycles each wavelet has: odd, from three to MOST_HALF_CYCLES.	`HALF_CYCLES`
dofs	sequence of str	Which of the requirement's DOFs, in this order; all by default.	None
progress	callable	Called with (done, total) a pass at a time, the total raised as each attempt begins; the app's progress bar uses it.	None

Returns:

Type	Description
ShockSynthesis	The record and each channel's wavelets, spectrum and reading.

Raises:

Type	Description
TypeError	If specification is not a ShockSpecification.
ValueError	If the requirement is not an acceleration SRS in declared units, names no DOF asked for, is not a level above zero at a point inside its band, or runs to RATE_RATIO times the sample rate's fifth; if the half cycles are not odd and from three to MOST_HALF_CYCLES; if no pass meets the requirement (naming the worst point); or if no synthesis fits the shaker (naming the limit that runs out and the band that asks for it).

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

def synthesize(specification: ShockSpecification, *, sample_rate: float,
               rating: ShakerRating | None = None,
               half_cycles: int = HALF_CYCLES,
               dofs: Sequence[str] | None = None,
               progress: Callable[[int, int], None] | None = None
               ) -> ShockSynthesis:
    """A shock a shaker can play, whose spectrum meets `specification`.

        made = synthesize(spec, sample_rate=20480.0,
                          rating=ShakerRating(stroke=0.0508, velocity=1.78,
                                              force=35600.0,
                                              moving_mass=45.0))
        made.history, made.tolerance(), made.channels[0].reading

    One wavelet per point of the requirement, their amplitudes iterated
    until every point of the record's spectrum lies inside the
    requirement's inner band; then the record is read against `rating`
    raw (`core.stroke.read_stroke`). While the shaker cannot hold it the
    wavelets are lengthened two half cycles at a time, up to
    `MOST_HALF_CYCLES`, and the first that fits is kept. The module
    docstring has the method.

    Parameters
    ----------
    specification : ShockSpecification
        The requirement, an acceleration SRS in declared units, with its
        Q and type; its warning limits are the band every point is put
        inside, `compliance.SRS_INNER_DB` without them.
    sample_rate : float
        Samples per second, above `RATE_RATIO` times the highest point.
    rating : ShakerRating, optional
        The shaker's shock ratings and the masses it carries; a rating
        left out is not judged, and with none the first synthesis that
        meets the requirement is kept.
    half_cycles : int, default `HALF_CYCLES`
        The fewest half cycles each wavelet has: odd, from three to
        `MOST_HALF_CYCLES`.
    dofs : sequence of str, optional
        Which of the requirement's DOFs, in this order; all by default.
    progress : callable, optional
        Called with ``(done, total)`` a pass at a time, the total raised
        as each attempt begins; the app's progress bar uses it.

    Returns
    -------
    ShockSynthesis
        The record and each channel's wavelets, spectrum and reading.

    Raises
    ------
    TypeError
        If `specification` is not a `ShockSpecification`.
    ValueError
        If the requirement is not an acceleration SRS in declared units,
        names no DOF asked for, is not a level above zero at a point
        inside its band, or runs to `RATE_RATIO` times the sample rate's
        fifth; if the half cycles are not odd and from three to
        `MOST_HALF_CYCLES`; if no pass meets the requirement (naming the
        worst point); or if no synthesis fits the shaker (naming the
        limit that runs out and the band that asks for it).
    """
    from ..units import UNKNOWN
    from .data import ShockSpecification
    from .progress import Ticker
    from .signals import _history

    if not isinstance(specification, ShockSpecification):
        raise TypeError('a shock is synthesized to a shock specification, '
                         f'not a {type(specification).__name__}')
    rate = float(sample_rate)
    if not (np.isfinite(rate) and rate > 0.0):
        raise ValueError('the sample rate has to be above zero')
    first = _half_cycles(half_cycles)
    if first > MOST_HALF_CYCLES:
        raise ValueError(f'{first} half cycles is past the most tried, '
                         f'{MOST_HALF_CYCLES}')
    rating = ShakerRating() if rating is None else rating
    own = [str(d) for d in specification.response_dof]
    wanted = (list(dict.fromkeys(own)) if dofs is None
              else [str(d) for d in dofs])
    if not wanted:
        raise ValueError('name at least one DOF to synthesize')
    twice = sorted({d for d in wanted if wanted.count(d) > 1})
    if twice:
        raise ValueError(f'{", ".join(twice)} is named more than once; a '
                         'record holds one channel a DOF')
    missing = [d for d in wanted if d not in own]
    if missing:
        raise ValueError(f'the requirement names no {", ".join(missing)}; '
                         f'it names {", ".join(own)}')
    ticker = Ticker(progress)
    rows, channels, records = [], [], []
    for dof in wanted:
        # the acceleration record at the DOF: one at a DOF can sit beside
        # another quantity's, and two of acceleration are two answers
        at = [i for i, d in enumerate(own) if d == dof]
        held = [i for i in at
                if specification.ordinate_dim[i] == 'acceleration']
        if len(held) > 1:
            raise ValueError(f'the requirement holds {len(held)} '
                             f'acceleration records at {dof}; keep one')
        if not held:
            dimension = specification.ordinate_dim[at[0]]
            if dimension == UNKNOWN:
                raise ValueError(f'the requirement at {dof} does not say '
                                 'what it is a spectrum of; declare its '
                                 'units first')
            raise ValueError(f'a shaker plays acceleration, and the '
                             f'requirement at {dof} is a spectrum of '
                             f'{dimension}')
        record = held[0]
        records.append(record)
        row, a = _channel(specification, record, rate, rating, first, ticker)
        rows.append(row)
        channels.append(SynthesizedChannel(
            dof=dof, points=a['points'], target=a['target'],
            frequency=a['tone'], samples=a['samples'], start=a['start'],
            amplitude=a['amplitude'], half_cycles=a['half_cycles'],
            passes=a['passes'], tried=a['tried'], spectrum=a['spectrum'],
            reading=a['reading']))
    # one clock for every channel: a shorter one keeps its quiet tail
    # longer, still at rest
    longest = max(row.size for row in rows)
    block = np.zeros((len(rows), longest))
    for i, row in enumerate(rows):
        block[i, :row.size] = row
    what = ', '.join(sorted({f'{c.half_cycles} half cycles'
                             for c in channels}))
    history = _history(block, rate, wanted, ['acceleration'] * len(rows),
                       f'SRS synthesis, wavelets of {what}')
    # shown in the requirement's units where it names them, held in SI
    # as every record is
    history.ordinate_unit = [specification.ordinate_unit[r] or shown
                             for r, shown in zip(records,
                                                 history.ordinate_unit,
                                                 strict=True)]
    return ShockSynthesis(history=history, specification=specification,
                          rating=rating, channels=channels)

suggested_rate

suggested_rate(specification: ShockSpecification) -> float

A sample rate for a synthesis: SUGGESTED_RATIO times the requirement's highest written point, rounded up to a multiple of 1024 Hz.

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

def suggested_rate(specification: ShockSpecification) -> float:
    """A sample rate for a synthesis: `SUGGESTED_RATIO` times the
    requirement's highest written point, rounded up to a multiple of
    1024 Hz."""
    f = np.asarray(specification.abscissa, dtype=float)
    level = np.abs(np.asarray(specification.ordinate, dtype=float))
    written = f[np.any(np.isfinite(level) & (level > 0.0), axis=0)]
    top = float(written.max()) if written.size else 2000.0
    return float(1024 * int(np.ceil(SUGGESTED_RATIO * top / 1024.0)))