API Reference ¶

fdscore.types.FDSTimePlan

Notes

The computation assumes a linear SDOF transfer from base acceleration to the selected metric and evaluates fatigue on the reconstructed response histories. For fixed slope_k, the absolute damage level scales globally as

\[\frac{p_{scale}^k}{N_{ref} S_{ref}^k}\]

Therefore p_scale, ref_stress, and ref_cycles change the magnitude of damage(f) but not its relative shape. When only the spectral shape and a compatible FDS-to-PSD inversion matter, a normalized S-N definition together with p_scale=1.0 is usually sufficient.

The choice of detrend can materially affect low-frequency damage, especially for displacement- and pseudo-velocity-based metrics, because offsets and slow drifts are amplified by the low-frequency dynamics of the oscillator bank.

If p_scale is omitted, p_scale=1.0 is assumed only for normalized S-N parameters with ref_stress=1 and ref_cycles=1. Physical workflows must pass p_scale explicitly.

References

ASTM E1049-85(2017). Standard Practices for Cycle Counting in Fatigue: Analysis.
Miner, M. A. (1945). “Cumulative Damage in Fatigue.” Journal of Applied: Mechanics, 12(3), A159-A164.
Crandall, S. H., & Mark, W. D. (1963). Random Vibrations in Mechanical: Systems. Academic Press.

fdscore.prepare_fds_time_plan(*, fs, n_samples, sdof, strict_nyquist=True)[source]¶

Precompute the FFT-domain transfer data for repeated FDS evaluations.

A time-domain FDS call repeatedly uses the same oscillator grid, damping, and FFT-bin transfer matrix when the sampling configuration is fixed. FDSTimePlan stores that reusable state so that repeated calls can skip transfer-matrix reconstruction.

Parameters:

fs (float) – Sampling rate in Hz.
n_samples (int) – Number of samples in the time histories that will reuse this plan.
sdof (SDOFParams) – Oscillator-grid definition and chosen response metric.
strict_nyquist (bool) – If True, frequencies above Nyquist raise ValidationError. If False, the frequency grid is truncated to the valid Nyquist range.

Returns:

Reusable transfer plan containing the validated oscillator grid, the implied damping ratio, and the complex FFT-domain transfer matrix.

Return type:

Notes

The stored matrix H is materialized as complex128 with shape (len(f0), n_fft_bins). Memory therefore scales approximately as

\[len(f_0) \times n_{fft\_bins} \times 16 \text{ bytes}\]

so plans trade memory for speed. For example, 400 oscillators and a 4 s signal sampled at 1 kHz require roughly 12 MB for the matrix alone.

fdscore.compute_fds_spectral_psd(f_psd_hz, psd_baseacc, *, duration_s, sn, sdof, p_scale=None)[source]¶

Compute a spectral Fatigue Damage Spectrum from an acceleration PSD.

This routine evaluates fatigue damage oscillator by oscillator using a PSD-domain response model and the library’s internal Dirlik spectral fatigue approximation.

Pipeline

The computation first builds the oscillator grid defined by sdof and the transfer matrix from base acceleration PSD to the selected response metric. For each oscillator, the response PSD is then computed as

\[P_{resp}(f; f_0) = p_{scale}^2 \left| H(f; f_0) \right|^2 P_{base}(f)\]

Dirlik’s method is applied to that response PSD to estimate life, and Miner damage is recovered through

\[D(f_0) = \frac{T}{life(f_0)}\]

Parameters:

f_psd_hz (numpy.ndarray) – One-sided input PSD frequency grid in Hz.
psd_baseacc (numpy.ndarray) – One-sided base-acceleration PSD defined on f_psd_hz.
duration_s (float) – Exposure duration associated with the PSD.
sn (fdscore.types.SNParams) – S-N curve definition used by the Dirlik life calculation.
sdof (fdscore.types.SDOFParams) – Oscillator-grid definition and response metric.
p_scale (float or None) – Additional scale factor applied to the response quantity before damage evaluation. This must match the fatigue convention used elsewhere in the workflow.

Returns:

Damage spectrum on the oscillator grid defined by sdof.

Return type:

Notes

Dirlik is a spectral fatigue approximation derived from response-PSD moments. It is not the same algorithm as rainflow counting on a realized time history, so absolute levels from compute_fds_spectral_psd(...) and compute_fds_time(...) should not be expected to match exactly.

Agreement with time-domain rainflow tends to improve when the response is well represented as a stationary Gaussian process and the record is long enough that the PSD is a stable descriptor. Differences grow for short records, strongly non-stationary environments, transient content, and non-Gaussian responses.

For fixed slope_k, p_scale, ref_stress, and ref_cycles act only as a global damage scaling factor. They change the magnitude of damage(f) but not its relative shape. Use SNParams.normalized(...) with p_scale=1.0 when a normalized workflow is sufficient.

Input PSD values are expected to be non-negative. Tiny negative values consistent with numerical noise are clamped to zero; materially negative values raise ValidationError.

References

Dirlik, T. (1985). Application of computers in fatigue analysis. Miner, M. A. (1945). “Cumulative Damage in Fatigue.” Journal of Applied Mechanics, 12(3), A159-A164.

fdscore.compute_fds_spectral_time(x, fs, *, sn, sdof, psd, duration_s=None, p_scale=None)[source]¶

Compute a spectral FDS from a time history through Welch plus Dirlik.

This convenience route first estimates a one-sided acceleration PSD from the input time history and then delegates to fdscore.compute_fds_spectral_psd().

Parameters:

x (numpy.ndarray) – One-dimensional base-acceleration time history.
fs (float) – Sampling rate in Hz.
sn (fdscore.types.SNParams) – S-N curve definition used for the fatigue calculation.
sdof (fdscore.types.SDOFParams) – Oscillator-grid definition and response metric.
psd (fdscore.types.PSDParams) – PSD-estimation configuration passed to compute_psd_welch(...).
duration_s (float or None) – Exposure duration associated with the damage estimate. If None, uses len(x) / fs.
p_scale (float or None) – Optional scale factor applied to the response quantity before damage evaluation.

Returns:

Spectral damage spectrum evaluated from the internally estimated PSD.

Return type:

Notes

This route combines two approximation layers:

Welch estimation of the PSD from a finite realization.
Dirlik spectral fatigue damage from the estimated PSD.

Differences relative to compute_fds_time(...) or to compute_fds_spectral_psd(...) with a reference PSD are therefore expected for finite-length signals.

References

Dirlik, T. (1985). Application of computers in fatigue analysis.

fdscore.scale_fds(fds, factor)[source]¶

Scale all damage ordinates of a fatigue damage spectrum.

Parameters:

fds (fdscore.types.FDSResult) – Input fatigue damage spectrum to be scaled.
factor (float) – Positive finite multiplicative factor applied uniformly to all damage values.

Returns:

New fatigue damage spectrum with the same frequency grid as the input and damage values multiplied by factor.

Return type:

Notes

This operation is purely algebraic. It is appropriate when a known post hoc scaling of damage is required, for example to account for an external calibration factor or to compare normalized spectra on a common basis. The function does not revisit the underlying stress cycles, SDOF response model, or S-N assumptions used to generate the original result.

The returned object stores provenance metadata indicating that the spectrum was produced by scale_fds(), together with the applied factor and the upstream provenance chain of the input spectrum.

fdscore.sum_fds(fds_list, weights=None)[source]¶

Form a weighted sum of mutually compatible fatigue damage spectra.

Parameters:

fds_list (collections.abc.Sequence of fdscore.types.FDSResult) – Sequence of fatigue damage spectra defined on compatible frequency grids and generated under the same compatibility contract.
weights (collections.abc.Sequence of float or None, optional) – Non-negative weights applied to each spectrum before summation. When omitted, unit weights are used for all inputs.

Returns:

Weighted sum of the input fatigue damage spectra on the reference frequency grid.

Return type:

fdscore.types.PSDResult

Notes

Compatibility is enforced internally before summation so that the spectra can be meaningfully combined without mixing damage metrics, oscillator assumptions, or incompatible S-N definitions. This function therefore acts as a safe superposition utility for spectra that already satisfy a shared engineering interpretation.

The summation is performed directly on the damage ordinates. No normalization is applied to the weights, so the caller controls whether the result represents a simple sum, a convex combination, or another weighted aggregate.

Provenance metadata is preserved and extended to record the number of inputs, the applied weights, and the provenance chain of each contributing spectrum.

fdscore.invert_fds_closed_form(fds, *, test_duration_s, strict_metric=True)[source]¶

Invert an FDS to an equivalent acceleration PSD using the closed-form Henderson-Piersol method.

Overview

The inversion proceeds in two steps. First, Miner damage is converted to a Damage Potential spectrum using the proportionality factor K derived from the S-N parameters and p_scale:

\[DP(f) = \frac{Damage(f)}{K}\]

Second, the closed-form inverse of H&P (1995), Eq. 12 is applied:

\[G_{aa}(f) = f \zeta \left[\frac{DP(f)}{f T}\right]^{2 / b}\]

The round-trip reconstruction error stored in meta["reconstruction"]["med_abs_log10"] is the median absolute value of

\[\log_{10}\left(\frac{D_{recon}}{D_{target}}\right)\]

Metric restriction

This method is valid only for metric="pv" (pseudo-velocity).

The physical justification follows from two independent lines of argument that converge on pseudo-velocity as the stress-relevant quantity:

Empirical / experimental: Gaberson & Chalmers (1969) established that modal velocity is the best single predictor of shock and vibration severity across a wide range of structures. Their work is explicitly cited by Henderson & Piersol (1995) as the basis for using velocity rather than acceleration or displacement to estimate stress.
Theoretical: Crandall (1962) derived analytically that, for a structure vibrating at resonance, peak stress is proportional to peak velocity, not peak acceleration or displacement. This result holds for geometrically simple structures (beams, plates under bending) and is the theoretical foundation cited by H&P (1995) in their Eq. 9.

For a lightly damped SDOF oscillator (\(\zeta < 0.1\)), pseudo-velocity

\[PV = 2 \pi f_0 x_{rel}\]

and relative velocity \(v_{rel}\) are approximately equal at resonance, because the phase between displacement and velocity at the resonant peak makes

\[|v_{rel}| \approx \omega_0 |x_{rel}|\]

This equivalence justifies using metric="pv" as a numerically convenient proxy for relative velocity in the closed-form derivation.

Using metric="acc", "vel", or "disp" would require a different transfer function at resonance and a different proportionality to stress, yielding a structurally different inversion equation that is not implemented here. Passing strict_metric=False suppresses the guard but does not make the derivation valid for those metrics.

Global damage scaling cancellation¶

When fds was computed with compatible settings, the absolute damage scaling carried by p_scale, ref_stress, and ref_cycles cancels exactly in the inversion: those parameters affect the magnitude of damage(f) and of K in equal proportion, leaving the inverted PSD unchanged. This invariance is verified by test_closed_form_psd_is_invariant_to_global_damage_scaling.

Requirements

fds.meta["compat"] must exist and contain metric, q, p_scale, and sn.
metric must be "pv" unless strict_metric=False.
test_duration_s must be the intended test duration, not the original signal duration.

type fds:: FDSResult
param fds:: Target FDS result. It must carry meta["compat"] produced by any compute_fds_* function in this library.
type fds:: fdscore.types.FDSResult
type test_duration_s:: float
param test_duration_s:: Duration of the equivalent test \(T\) in seconds. This is the denominator in the DP-to-PSD step and directly controls the amplitude of the inverted PSD.
type test_duration_s:: float
type strict_metric:: bool
param strict_metric:: If True (default), raise ValidationError when fds.meta["compat"]["metric"] != "pv".
type strict_metric:: bool
returns:: Equivalent acceleration PSD on the same frequency grid as fds.f. meta["reconstruction"] contains the round-trip log10 error used for quality assessment.
rtype:: fdscore.types.PSDResult

References

Henderson, G. R. & Piersol, A. G. (1995). “Fatigue Damage Related: Descriptor for Random Vibration Test Environments.” Sound and Vibration, October 1995, 20-24. Equations 11-12.
Gaberson, H. A. & Chalmers, R. H. (1969). “Modal Velocity as a: Criterion of Shock Severity.” Shock and Vibration Bulletin, No. 40, Pt. 2, 31-49.
Crandall, S. H. (1962). “Relationship between Stress and Velocity in: Resonant Vibration.” Journal of the Acoustical Society of America, 34(12), 1960-1961.
Crandall, S. H. & Mark, W. D. (1963). Random Vibrations in Mechanical: Systems. Academic Press, New York. pp. 113-117.

fdscore.invert_fds_iterative_spectral(target, *, f_psd_hz, psd_seed, duration_s, sn, sdof, p_scale, params=IterativeInversionParams(iters=30, gamma=0.8, gain_min=0.2, gain_max=5.0, alpha_sharpness=1.0, floor=1e-30, smooth_enabled=True, smooth_window_bins=11, smooth_every_n_iters=1, prior_blend=0.0, prior_power=1.0, edge_anchor_hz=0.0, edge_anchor_blend=0.0, tail_cap_start_hz=0.0, tail_cap_ratio=0.0, low_cap_ratio=0.0, post_smooth_window_bins=0, post_smooth_blend=1.0, post_refine_iters=0, post_refine_gamma=0.5, post_refine_min=0.7, post_refine_max=2.2))[source]¶

Iteratively synthesize a PSD that matches a target FDS with a spectral predictor.

This routine solves the inverse problem “find an acceleration PSD whose predicted FDS matches target” by repeatedly calling the spectral predictor fdscore.compute_fds_spectral_psd().

Algorithm

Let \(F_{target}\) be the target damage spectrum and let \(F(P)\) be the spectral predictor evaluated at candidate PSD \(P\). The method builds a PSD-to-oscillator influence matrix from the SDOF transfer model and converts it to a normalized redistribution matrix \(\alpha\).

Oscillator-wise multiplicative gains are computed as

\[s_i = \left(\frac{F_{target, i}}{F_i(P)}\right)^{2 / k}\]

and clipped to the interval defined by gain_min and gain_max. Those gains are projected back to PSD bins through

\[u = \exp\left(\alpha^T \log(s)\right)\]

The candidate PSD is then updated multiplicatively as

\[P \leftarrow P \, u^{\gamma}\]

followed by optional smoothing, prior blending, and edge capping. The updated PSD is re-evaluated and scored by the median absolute log10-domain mismatch.

The function always outputs a one-sided acceleration PSD on the user supplied f_psd_hz grid.

Parameters:

target (fdscore.types.FDSResult) – Target FDS result, expressed as damage versus oscillator frequency.
f_psd_hz (numpy.ndarray) – Frequency grid in Hz for the synthesized acceleration PSD.
psd_seed (numpy.ndarray) – Strictly positive seed PSD defined on f_psd_hz.
duration_s (float) – Exposure duration used by the spectral predictor.
sn (fdscore.types.SNParams) – S-N curve definition used by the predictor.
sdof (fdscore.types.SDOFParams) – Oscillator-grid definition and response metric used by the predictor.
p_scale (float) – Response scale factor used by the predictor. It must be compatible with the way target was computed.
params (fdscore.types.IterativeInversionParams) – Iteration and regularization parameters.

Returns:

Synthesized acceleration PSD on f_psd_hz. The result metadata includes convergence diagnostics, reconstructed FDS traces, and explicit per-engine parameter usage.

Return type:

Notes

The update strategy implemented here is a library-specific multiplicative inversion heuristic. It is not claimed as a closed-form method from the fatigue literature.

The convergence metric stored in meta["diagnostics"]["best_err"] is the median absolute log10-domain mismatch over bins where both target and predicted damage are positive:

\[median\left(\left|\log_{10} F_{pred} - \log_{10} F_{target}\right|\right)\]

The main loop performs two predictor evaluations per iteration: one on the current PSD to derive oscillator-wise gains, and one after the multiplicative update to score the candidate that may become the new best solution.

If post-smoothing and post-refinement are enabled, an additional spectral inversion pass is launched from the smoothed candidate using a reduced parameter set. The metadata distinguishes whether the best result came from the main loop, the post-smooth stage, or the post-refine stage.

References

Dirlik, T. (1985). Application of computers in fatigue analysis. Miner, M. A. (1945). “Cumulative Damage in Fatigue.” Journal of Applied Mechanics, 12(3), A159-A164.

fdscore.invert_fds_iterative_time(target, *, f_psd_hz, psd_seed, fs, duration_s, sn, sdof, p_scale, params=IterativeInversionParams(iters=30, gamma=0.8, gain_min=0.2, gain_max=5.0, alpha_sharpness=1.0, floor=1e-30, smooth_enabled=True, smooth_window_bins=11, smooth_every_n_iters=1, prior_blend=0.0, prior_power=1.0, edge_anchor_hz=0.0, edge_anchor_blend=0.0, tail_cap_start_hz=0.0, tail_cap_ratio=0.0, low_cap_ratio=0.0, post_smooth_window_bins=0, post_smooth_blend=1.0, post_refine_iters=0, post_refine_gamma=0.5, post_refine_min=0.7, post_refine_max=2.2), n_realizations=1, seed=0, nfft=None, target_duration_s=None)[source]¶

Iteratively synthesize a PSD that matches a target FDS with a time predictor.

This routine solves the inverse problem “find an acceleration PSD whose time-domain FDS matches target” by repeatedly synthesizing time histories from the candidate PSD and re-evaluating fdscore.compute_fds_time().

Algorithm

Let \(F_{target}\) be the target damage spectrum and let \(F(P)\) be the average time-domain predictor response for candidate PSD \(P\). The method builds a PSD-to-oscillator influence matrix from the SDOF transfer model and converts it to a normalized redistribution matrix \(\alpha\).

For each predictor call, the routine synthesizes n_realizations random-phase time histories from \(P\), evaluates their FDS with fdscore.compute_fds_time(), and averages the resulting damage spectra.

Oscillator-wise multiplicative gains are computed as

\[s_i = \left(\frac{F_{target, i}}{F_i(P)}\right)^{2 / k}\]

and clipped to the interval defined by gain_min and gain_max. Those gains are projected back to PSD bins through

\[u = \exp\left(\alpha^T \log(s)\right)\]

The candidate PSD is then updated multiplicatively as

\[P \leftarrow P \, u^{\gamma}\]

followed by optional smoothing and prior blending. The updated PSD is re-evaluated and scored by the median absolute log10-domain mismatch.

Parameters:

target (fdscore.types.FDSResult) – Target FDS result. It must carry compatibility metadata.
f_psd_hz (numpy.ndarray) – Frequency grid in Hz for the synthesized acceleration PSD.
psd_seed (numpy.ndarray) – Strictly positive seed PSD defined on f_psd_hz.
fs (float) – Sampling rate in Hz used during time synthesis and FDS evaluation.
duration_s (float) – Synthetic duration in seconds used for each predictor realization.
sn (fdscore.types.SNParams) – S-N curve definition used by the time-domain predictor.
sdof (fdscore.types.SDOFParams) – Oscillator-grid definition and response metric used by the predictor.
p_scale (float) – Response scale factor used by the predictor. It must be compatible with the way target was computed.
params (fdscore.types.IterativeInversionParams) – Iteration and regularization parameters.
n_realizations (int) – Number of random-phase synthesized histories averaged per predictor call.
seed (int or None) – Seed for reproducible synthesis. Realization r uses seed + r.
nfft (int or None) – FFT length used during synthesis. If None, the synthesis routine chooses the next power of two.
target_duration_s (float or None) –
Optional duration to which the predictor damage is rescaled through

\[D_{scaled} = D_{synth} \frac{T_{target}}{T_{synth}}\]

If None, no duration scaling is applied.

Returns:

Synthesized acceleration PSD on f_psd_hz. The result metadata includes convergence diagnostics, predictor configuration, and explicit per-engine parameter usage.

Return type:

fdscore.types.PSDResult

Notes

The inversion heuristic implemented here is library-specific. It combines a stochastic time-history predictor with a multiplicative PSD update and does not correspond to a published closed-form inversion formula.

The convergence metric stored in meta["diagnostics"]["best_err"] is the median absolute log10-domain mismatch over bins where both target and predicted damage are positive.

The main loop performs two predictor evaluations per iteration: one on the current PSD to derive oscillator-wise gains, and one after the update to score the candidate that may become the best solution.

Because the predictor is stochastic, convergence depends on n_realizations, seed, duration_s, and nfft. The fixed internal predictor policy is recorded in meta["diagnostics"]["predictor_config"] and currently uses synthesize_time_from_psd(remove_mean=True) together with compute_fds_time(detrend="none", batch_size=64).

References

ASTM E1049-85(2017). Standard Practices for Cycle Counting in Fatigue Analysis. Miner, M. A. (1945). “Cumulative Damage in Fatigue.” Journal of Applied Mechanics, 12(3), A159-A164.

PSD Analysis and Signal Synthesis¶

These functions estimate acceleration PSDs, derive scalar PSD metrics, and synthesize stationary Gaussian realizations.

fdscore.compute_psd_welch(x, fs, psd)[source]¶

Compute a one-sided PSD using Welch’s method.

Parameters:

x (numpy.ndarray) – One-dimensional input time history.
fs (float) – Sampling rate in Hz.
psd (fdscore.types.PSDParams) – PSD-estimation settings. The current implementation supports only method="welch".

Return type:

tuple[ndarray, ndarray]

Returns:

f_hz (numpy.ndarray) – One-dimensional frequency grid in Hz.
psd_values (numpy.ndarray) – One-dimensional PSD values in units squared per hertz.

Notes

This helper is purely numerical and performs no file I/O.

The output is the one-sided PSD returned by Welch’s method under the requested detrending, window, overlap, and segment-length settings. When psd.fmin or psd.fmax is provided, the spectrum is cropped after estimation.

Tiny negative PSD values caused by floating-point noise are clipped to zero. Materially negative values are rejected because they would invalidate downstream log-domain operations and stochastic synthesis.

References

Welch, P. D. (1967). The use of the fast Fourier transform for the estimation of power spectra: A method based on time averaging over short, modified periodograms. IEEE Transactions on Audio and Electroacoustics, 15(2), 70-73.

fdscore.compute_psd_metrics(psd, *, f_hz=None, duration_s=None, acc_unit=None, acc_to_m_s2=None, bands_hz=((5.0, 20.0), (20.0, 80.0), (80.0, 200.0), (200.0, 400.0)))[source]¶

Compute scalar summary metrics from an acceleration PSD.

The routine derives RMS, Gaussian peak estimates, velocity and displacement metrics, and band-limited RMS values from a one-sided acceleration PSD.

Parameters:

psd (fdscore.types.PSDResult or numpy.ndarray) – Either a PSDResult or a one-dimensional PSD array.
f_hz (numpy.ndarray or None) – Frequency vector in Hz when psd is passed as a raw array. Must be omitted when psd is a PSDResult.
duration_s (float or None) – Exposure duration used for Gaussian peak estimation. If None, peak-related outputs are returned as nan.
acc_unit (str or None) – Acceleration unit of the PSD values. Supported values are "g" and "m/s2".
acc_to_m_s2 (float or None) – Explicit conversion factor from the PSD acceleration unit to \(m/s^2\). If provided, overrides acc_unit.
bands_hz (collections.abc.Sequence) – Sequence of (f_lo, f_hi) frequency-band pairs used to compute band-limited RMS acceleration in g.

Returns:

Structured set of scalar metrics including acceleration RMS and peak, zero-upcrossing rate, effective cycles, velocity metrics, displacement metrics, and per-band RMS values.

Return type:

fdscore.types.PSDMetricsResult

Notes

RMS values are obtained by integrating the one-sided PSD. Velocity and displacement PSDs are derived from acceleration through division by \(\omega^2\) and \(\omega^4\), respectively, for nonzero angular frequency.

When duration_s is provided, the peak factor is estimated with a Gaussian extreme-value approximation of Davenport form:

\[\nu_0 = \frac{1}{2 \pi} \sqrt{\frac{m_2}{m_0}}\]

\[n_{eff} = \nu_0 T\]

\[u = \sqrt{2 \ln(n_{eff})}\]

\[peak\_factor \approx u + \frac{\gamma_E}{u}\]

where \(m_0\) and \(m_2\) are spectral moments and \(\gamma_E\) is the Euler-Mascheroni constant. The implementation floors \(n_{eff}\) at \(e\) to avoid a singular or negative argument in the logarithm for very short exposures or very low zero-upcrossing rates.

These peak estimates assume a stationary Gaussian response process. They should be interpreted cautiously for strongly non-stationary environments, deterministic transients, or markedly non-Gaussian signals.

References

Davenport, A. G. (1964). “Note on the Distribution of the Largest Value of: a Random Function with Application to Gust Loading.” Proceedings of the Institution of Civil Engineers, 28(2), 187-196.

fdscore.synthesize_time_from_psd(*, f_psd_hz, psd, fs, duration_s, seed=None, nfft=None, remove_mean=True)[source]¶

Synthesize a stationary Gaussian time history from a one-sided PSD.

Parameters:

f_psd_hz (numpy.ndarray) – One-dimensional one-sided PSD frequency grid in Hz.
psd (numpy.ndarray) – One-dimensional one-sided PSD values in units squared per hertz.
fs (float) – Sampling rate in Hz.
duration_s (float) – Desired duration in seconds. The output length is N = round(duration_s * fs).
seed (int or None) – Random seed used for reproducible phase generation.
nfft (int or None) – FFT length used for synthesis. If None, the next power of two greater than or equal to N is used. If provided, it must satisfy nfft >= N.
remove_mean (bool) – If True, subtract the sample mean after synthesis.

Returns:

Synthesized time history of length N.

Return type:

numpy.ndarray

Notes

This routine is typically used by iterative inversion predictors that map a PSD to a synthetic time history.

The generated realization is stationary and Gaussian by construction and preserves the target PSD only in the statistical, ensemble sense. It does not reproduce transient structure, non-stationarity, or non-Gaussian tails.

The PSD is first interpolated onto the FFT frequency bins associated with nfft. Random phases are then assigned to the interior one-sided bins, while the singleton DC and Nyquist bins are handled separately to preserve the Hermitian structure required by the real inverse FFT.

For finite durations, Welch estimates from the synthesized signal will not match the input PSD exactly point by point.

Deterministic Harmonic Workflows¶

These APIs cover single-tone, dwell-profile, and dwell-discretized sweep analyses.

fdscore.compute_ers_sine(*, freq_hz, amp, sdof, input_motion='acc', peak_mode='abs')[source]¶

Compute deterministic ERS for a single harmonic base excitation.

The returned spectrum is the response-amplitude spectrum of the selected SDOF metric. For a pure sinusoid, the absolute peak response is equal to the response amplitude, so the ERS is obtained directly from the harmonic transfer magnitude.

Parameters:

freq_hz (float) – Excitation frequency of the sinusoid in Hz.
amp (float) – Input amplitude expressed in the units implied by input_motion.
sdof (SDOFParams) – Oscillator-grid definition and response metric.
input_motion (Literal['acc', 'vel', 'disp']) –
Type of the specified harmonic amplitude. The routine converts the input to equivalent base-acceleration amplitude through

\[\begin{split}a_{base,amp} = \begin{cases} amp, & \text{if input\_motion = "acc"} \\ \omega \, amp, & \text{if input\_motion = "vel"} \\ \omega^2 \, amp, & \text{if input\_motion = "disp"} \end{cases}\end{split}\]

where \(\omega = 2 \pi f\).
peak_mode (Literal['abs']) – Peak convention for the ERS. Only "abs" is currently supported.

Returns:

Response-amplitude spectrum on the oscillator grid.

Return type:

Notes

The result is deterministic and contains no cycle counting or statistical assumptions. It is obtained from the magnitude of the linear SDOF transfer function at the single forcing frequency.

fdscore.compute_fds_sine(*, freq_hz, amp, duration_s, sn, sdof, input_motion='acc', p_scale=None)[source]¶

Compute deterministic FDS for a single harmonic base excitation.

This routine evaluates the damage caused by a sinusoidal base input of frequency freq_hz acting for duration_s seconds. Each oscillator is driven at the same forcing frequency and its steady-state response amplitude is converted directly to Miner damage.

Parameters:

freq_hz (float) – Excitation frequency of the sinusoid in Hz.
amp (float) – Input amplitude expressed in the units implied by input_motion.
duration_s (float) – Exposure duration in seconds.
sn (fdscore.types.SNParams) – S-N curve definition used for Miner damage accumulation.
sdof (fdscore.types.SDOFParams) – Oscillator-grid definition and response metric.
input_motion (str) – Optional input-motion type. Accepted values are "acc", "vel", and "disp". Internally the routine converts velocity and displacement amplitudes to equivalent base-acceleration amplitude using \(\omega amp\) and \(\omega^2 amp\), respectively.
p_scale (float or None) – Optional response scale factor applied before damage evaluation.

Returns:

Deterministic damage spectrum on the oscillator grid.

Return type:

Notes

For a single harmonic input, the damage model reduces to

\[D(f_0) = N_{cycles} \frac{\left(p_{scale} S(f_0)\right)^k}{C}\]

where \(N_{cycles} = f \, T\), \(S(f_0)\) is the response load amplitude derived from the chosen metric, and \(C\) is the S-N intercept. When sn.amplitude_from_range is False, the full range 2 S(f_0) is used instead of alternating amplitude.

References

Miner, M. A. (1945). “Cumulative Damage in Fatigue.” Journal of Applied: Mechanics, 12(3), A159-A164.

fdscore.compute_ers_dwell_profile(segments, *, sdof, peak_mode='abs')[source]¶

Compute mission-level ERS as the envelope of harmonic dwell segments.

Parameters:

segments (Sequence[SineDwellSegment]) – Sequence of deterministic harmonic dwell segments.
sdof (SDOFParams) – Oscillator-grid definition and response metric.
peak_mode (Literal['abs']) – Peak convention for the ERS. Only "abs" is currently supported.

Returns:

Pointwise envelope of the ERS results produced by the individual segments.

Return type:

Notes

Each segment is evaluated independently with compute_ers_sine(), and the final result is the pointwise maximum across segments on a compatible oscillator grid.

fdscore.compute_fds_dwell_profile(segments, *, sn, sdof, p_scale=None)[source]¶

Compute mission-level FDS as the sum of harmonic dwell-segment damage.

Parameters:

segments (collections.abc.Sequence) – Sequence of fdscore.types.SineDwellSegment dwell segments.
sn (fdscore.types.SNParams) – S-N curve definition used for Miner damage accumulation.
sdof (fdscore.types.SDOFParams) – Oscillator-grid definition and response metric.
p_scale (float or None) – Optional response scale factor applied before damage evaluation.

Returns:

Damage spectrum equal to the sum of the per-segment deterministic FDS results.

Return type:

Notes

The function assumes linear cumulative damage and therefore sums the independent damage spectra computed for each dwell segment.

References

Miner, M. A. (1945). “Cumulative Damage in Fatigue.” Journal of Applied: Mechanics, 12(3), A159-A164.

fdscore.compute_ers_sine_sweep(*, f_start_hz, f_stop_hz, amp, duration_s, sdof, input_motion='acc', peak_mode='abs', spacing='log', n_steps=200)[source]¶

Approximate deterministic ERS for a sine sweep via dwell discretization.

Parameters:

f_start_hz (float) – Sweep start frequency in Hz.
f_stop_hz (float) – Sweep stop frequency in Hz.
amp (float) – Input amplitude expressed in the units implied by input_motion.
duration_s (float) – Total sweep duration in seconds.
sdof (SDOFParams) – Oscillator-grid definition and response metric.
input_motion (Literal['acc', 'vel', 'disp']) – Type of the specified harmonic amplitude.
peak_mode (Literal['abs']) – Peak convention for the ERS. Only "abs" is currently supported.
spacing (Literal['linear', 'log']) – Frequency spacing of the dwell discretization, either "linear" or "log".
n_steps (int) – Number of dwell segments used to approximate the sweep.

Returns:

Deterministic ERS approximation of the sweep.

Return type:

Notes

The sweep is approximated as a sequence of constant-frequency dwell segments of equal duration. Convergence improves as n_steps increases, especially near sharp resonances and for logarithmic sweeps spanning large frequency ranges.

fdscore.compute_fds_sine_sweep(*, f_start_hz, f_stop_hz, amp, duration_s, sn, sdof, input_motion='acc', p_scale=None, spacing='log', n_steps=200)[source]¶

Approximate deterministic FDS for a sine sweep via dwell discretization.

Parameters:

f_start_hz (float) – Sweep start frequency in Hz.
f_stop_hz (float) – Sweep stop frequency in Hz.
amp (float) – Input amplitude expressed in the units implied by input_motion.
duration_s (float) – Total sweep duration in seconds.
sn (fdscore.types.SNParams) – S-N curve definition used for Miner damage accumulation.
sdof (fdscore.types.SDOFParams) – Oscillator-grid definition and response metric.
input_motion (str) – Optional input-motion type. Accepted values are "acc", "vel", and "disp".
p_scale (float or None) – Optional response scale factor applied before damage evaluation.
spacing (str) – Optional spacing mode for the dwell discretization. Accepted values are "linear" and "log".
n_steps (int) – Number of dwell segments used to approximate the sweep.

Returns:

Deterministic FDS approximation of the sweep.

Return type:

Notes

The result is obtained by discretizing the sweep into harmonic dwell segments and summing their individual damage spectra. Convergence with respect to n_steps should be checked whenever the sweep crosses narrow resonances or when the requested damage estimate is used quantitatively.

References

Miner, M. A. (1945). “Cumulative Damage in Fatigue.” Journal of Applied: Mechanics, 12(3), A159-A164.

ERS Workflows¶

These functions compute and combine generic extreme-response spectra.

fdscore.compute_ers_time(x, fs, sdof, *, detrend='linear', strict_nyquist=True, batch_size=64, peak_mode='abs', plan=None, engine='incremental', zoh_r_max=0.2)[source]¶

Compute time-domain ERS by evaluating an SDOF bank on an input history.

Parameters:

x (numpy.ndarray) – One-dimensional input time history.
fs (float) – Sampling rate in Hz.
sdof (fdscore.types.SDOFParams) – Oscillator-grid definition and response metric.
detrend ({"linear", "mean", "none"}, optional) – Preprocessing mode applied before the FFT.
strict_nyquist (bool, optional) – Whether oscillator frequencies at or above Nyquist should raise an error instead of being clipped.
batch_size (int, optional) – Number of oscillator responses reconstructed simultaneously in each FFT batch. Only used when engine="fft".
peak_mode ({"abs"}, optional) – Peak convention for the ERS. The current implementation supports only absolute peaks.
plan (fdscore.types.FDSTimePlan or None, optional) – Optional reusable transfer plan matching the current sampling and oscillator configuration. Only used when engine="fft".
engine ({"fft", "incremental"}, optional) –
Internal computation engine.

"fft"
Original FFT-domain engine. Applies the continuous SDOF transfer function to rfft(x) and reconstructs each oscillator response with irfft in batches.

"incremental" (default)
Sample-by-sample SDOF integration using exact ZOH state- transition matrices. Peaks are tracked online, so the full (n_osc, n_samples) response matrix is never materialized. Oscillators near Nyquist are integrated with adaptive upsampling controlled by zoh_r_max.
zoh_r_max (float, optional) – Maximum tolerated f0 / Nyquist_effective ratio for the incremental ZOH engine. Smaller values reduce the high-frequency ZOH attenuation error at the cost of larger upsample factors. Ignored when engine="fft".

Returns:

Extreme-response spectrum evaluated on the validated oscillator grid.

Return type:

Notes

The ERS is tied to the selected sdof.metric and therefore may represent absolute acceleration, relative displacement, relative velocity, or pseudo-velocity depending on the chosen SDOF setup.

A compatible FDSTimePlan can be reused by the FFT engine because it stores only the transfer data for a fixed (fs, n_samples, metric, q, f-grid) configuration. Reuse avoids rebuilding the transfer matrix on repeated calls with the same sampling contract.

The computational pipeline is:

Validate the SDOF definition and requested frequency grid.
Enforce or apply Nyquist handling.
Preprocess the input signal.
Evaluate the oscillator bank with the selected engine.
Extract the peak absolute response for each oscillator.

The incremental engine is designed to match the FFT baseline while avoiding the cost of inverse FFT reconstruction for every oscillator. Small differences near the top of the frequency grid remain possible because the two engines use different discretization schemes (continuous FRF on discrete FFT bins versus causal ZOH state integration).

fdscore.compute_ers_spectral_psd(f_psd_hz, psd_baseacc, *, duration_s, sdof, nyquist_hz=None, edge_correction=True)[source]¶

Compute a spectral/random ERS from a one-sided acceleration PSD.

This workflow interprets the ERS as the expected extreme response of a stationary Gaussian process over duration_s. It is therefore a different quantity from fdscore.compute_ers_time(), which extracts the maximum observed in a realized time history.

For generic metrics, the response PSD is built as

\[P_{resp}(f; f_n) = \left| H(f; f_n) \right|^2 P_{base}(f)\]

and the response moments

\[m_0 = \int P_{resp}(f) \, df\]

\[m_2 = \int (2\pi f)^2 P_{resp}(f) \, df\]

are used to estimate the peak rate and the expected Gaussian maximum.

For sdof.metric="acc", the implementation uses a Lalanne-style relative-displacement/random-peak formulation that is better aligned with classical random-vibration ERS practice for acceleration spectra. Other metrics use the exact response-PSD moment route.

Parameters:

f_psd_hz (numpy.ndarray) – One-sided PSD frequency grid in Hz.
psd_baseacc (numpy.ndarray) – One-sided base-acceleration PSD defined on f_psd_hz.
duration_s (float) – Exposure duration associated with the expected extreme response.
sdof (fdscore.types.SDOFParams) – Oscillator-grid definition and response metric.
nyquist_hz (float or None, optional) – Original Nyquist limit of the underlying time-history sampling. This is used only when edge_correction=True and the PSD has been cropped below that limit.
edge_correction (bool, optional) – Whether to apply an automatic high-frequency edge correction for cropped PSDs. The correction is a per-oscillator raised-cosine taper over a bandwidth proportional to zeta * f_n and is a no-op when the PSD already reaches nyquist_hz.

Returns:

Spectral/random ERS evaluated on the oscillator grid defined by sdof.

Return type:

Notes

This route assumes that the PSD is an adequate descriptor of the environment and that Gaussian extreme-value statistics are an acceptable approximation for the response process.

The high-frequency edge correction exists because PSD exports are often cropped below the original Nyquist frequency for plotting or exchange. Near the top of the oscillator grid, that cropping can artificially suppress the response moments unless the missing tail is reconstructed.

References

Davenport, A. G. (1964). Note on the Distribution of the Largest Value of a Random Function with Application to Gust Loading. Lalanne, C. Mechanical Vibration and Shock Analysis, Volume 3 (Random Vibration).

fdscore.compute_ers_spectral_time(x, fs, *, sdof, psd, duration_s=None, edge_correction=True)[source]¶

Compute a spectral/random ERS from a time history through Welch.

This convenience route first estimates a one-sided acceleration PSD with fdscore.compute_psd_welch() and then delegates to fdscore.compute_ers_spectral_psd().

Parameters:

x (numpy.ndarray) – One-dimensional base-acceleration time history.
fs (float) – Sampling rate in Hz.
sdof (fdscore.types.SDOFParams) – Oscillator-grid definition and response metric.
psd (fdscore.types.PSDParams) – PSD-estimation settings passed to compute_psd_welch(...).
duration_s (float or None, optional) – Exposure duration associated with the expected extreme response. If None, uses len(x) / fs.
edge_correction (bool, optional) – Whether to enable the automatic high-frequency edge correction when the internally estimated PSD is cropped below Nyquist.

Returns:

Spectral/random ERS evaluated from the internally estimated PSD.

Return type:

Notes

This function is a convenience wrapper. It is not equivalent to fdscore.compute_ers_time(), which operates on the realized time history directly and extracts the maximum observed response.

fdscore.envelope_ers(results)[source]¶

Compute a pointwise envelope across compatible ERS results.

Parameters:: results (collections.abc.Sequence of fdscore.types.ERSResult) – Response spectra to combine. All inputs must share the same frequency grid and ERS compatibility signature.
Returns:: Envelope spectrum whose response at each oscillator frequency is the maximum of the corresponding input values.
Return type:: fdscore.types.ERSResult

Notes

Compatibility is enforced internally, so the function will reject inputs that mix different response metrics, peak conventions, oscillator assumptions, or incompatible frequency grids.

The returned result preserves the reference grid and stores provenance metadata recording the number and origin of contributing spectra.

Shock Workflows¶

These functions cover transient shock analysis, event-based workflows, rolling spectra, sided envelopes, and half-sine reduction.

fdscore.compute_srs_time(x, fs, sdof, *, detrend='mean', strict_nyquist=True, peak_mode='abs')[source]¶

Compute a shock response spectrum (SRS) from a base-acceleration time history.

Parameters:

x (numpy.ndarray) – One-dimensional base-acceleration time history.
fs (float) – Sampling rate in Hz.
sdof (fdscore.types.SDOFParams) – Oscillator-grid definition. For this wrapper, sdof.metric must be "acc".
detrend ({"linear", "mean", "median", "none"}, optional) – Detrending mode applied before the shock-spectrum calculation.
strict_nyquist (bool, optional) – Whether oscillator frequencies at or above Nyquist should raise an error instead of being clipped.
peak_mode ({"abs", "pos", "neg", "both"}, optional) – Requested peak convention. "both" returns separate negative and positive spectra.

Returns:

One-sided results are returned as ERSResult. When peak_mode="both", the function returns a ShockSpectrumPair containing the negative and positive sides.

Return type:

Notes

This wrapper uses the dedicated recursive shock engine implemented in fdscore._shock_iir.

The returned spectrum is tagged with ers_kind="shock_response_spectrum" so that downstream envelope and compatibility operations can distinguish it from generic ERS or PVSS results.

fdscore.compute_pvss_time(x, fs, sdof, *, detrend='mean', strict_nyquist=True, peak_mode='abs')[source]¶

Compute a pseudo-velocity shock spectrum (PVSS) from a base-acceleration time history.

Parameters:

x (numpy.ndarray) – One-dimensional base-acceleration time history.
fs (float) – Sampling rate in Hz.
sdof (fdscore.types.SDOFParams) – Oscillator-grid definition. For this wrapper, sdof.metric must be "pv".
detrend ({"linear", "mean", "median", "none"}, optional) – Detrending mode applied before the shock-spectrum calculation.
strict_nyquist (bool, optional) – Whether oscillator frequencies at or above Nyquist should raise an error instead of being clipped.
peak_mode ({"abs", "pos", "neg", "both"}, optional) – Requested peak convention. "both" returns separate negative and positive spectra.

Returns:

One-sided results are returned as ERSResult. When peak_mode="both", the function returns a ShockSpectrumPair containing the negative and positive sides.

Return type:

Notes

This wrapper uses the dedicated recursive shock engine implemented in fdscore._shock_iir.

The returned spectrum is tagged with ers_kind="pseudo_velocity_shock_spectrum" so that downstream tooling can distinguish PVSS results from generic ERS and classical acceleration-based SRS.

fdscore.detect_shock_events(x, fs, *, detrend='median', polarity='abs', threshold_reference='rms', threshold_multiplier=5.0, threshold_value=None, min_separation_s=0.05, window_s=None)[source]¶

Detect discrete shock events in a 1D acceleration time history.

Parameters:

x (numpy.ndarray) – One-dimensional acceleration time history.
fs (float) – Sampling rate in Hz.
detrend ({"linear", "mean", "median", "none"}, optional) – Preprocessing mode applied before peak detection.
polarity ({"abs", "pos", "neg"}, optional) – Whether detection is performed on absolute, positive-only, or negative-only peaks.
threshold_reference ({"rms", "std", "peak"}, optional) – Reference statistic used when threshold_value is not provided.
threshold_multiplier (float, optional) – Positive multiplier applied to the selected threshold reference.
threshold_value (float or None, optional) – Explicit absolute threshold. When provided, it overrides the reference-based threshold calculation.
min_separation_s (float, optional) – Minimum temporal separation enforced between detected peaks.
window_s (float or None, optional) – Total extraction-window duration centered on each detected peak. When omitted, the detector derives the window from min_separation_s.

Returns:

Detected events returned as a ShockEventSet, together with sampling metadata and detector settings.

Return type:

Notes

This is an axis-first detector intended for event-oriented shock workflows. It is most useful when a long recording contains multiple isolated transients and each transient should be analyzed with a dedicated local spectrum.

The stored event windows are index based. They can therefore be used directly to slice the original signal without recomputing detection state.

fdscore.compute_rolling_srs_time(x, fs, sdof, events, *, detrend='median', strict_nyquist=True, peak_mode='abs')[source]¶

Compute SRS on each detected shock-event window.

Parameters:

x (numpy.ndarray) – One-dimensional source signal containing the detected events.
fs (float) – Sampling rate in Hz.
sdof (fdscore.types.SDOFParams) – Oscillator-grid definition. For this wrapper, sdof.metric must be "acc".
events (object) – Detected event windows that define the rolling segments. The input must be a ShockEventSet instance.
detrend ({"linear", "mean", "median", "none"}, optional) – Preprocessing mode applied independently to each event window.
strict_nyquist (bool, optional) – Whether oscillator frequencies at or above Nyquist should raise an error instead of being clipped.
peak_mode ({"abs", "pos", "neg"}, optional) – Peak convention used for each event-local spectrum.

Returns:

Rolling result returned as a RollingERSResult. Each row corresponds to one detected event and each column corresponds to one oscillator frequency.

Return type:

Notes

Rolling shock spectra are currently event-window based rather than fixed-hop or sliding-window based. The time coordinate stored in the result corresponds to the detected peak time of each event.

If no events are present, the function still returns a valid empty matrix with the validated oscillator frequency grid.

fdscore.compute_rolling_pvss_time(x, fs, sdof, events, *, detrend='median', strict_nyquist=True, peak_mode='abs')[source]¶

Compute PVSS on each detected shock-event window.

Parameters:

x (numpy.ndarray) – One-dimensional source signal containing the detected events.
fs (float) – Sampling rate in Hz.
sdof (fdscore.types.SDOFParams) – Oscillator-grid definition. For this wrapper, sdof.metric must be "pv".
events (object) – Detected event windows that define the rolling segments. The input must be a ShockEventSet instance.
detrend ({"linear", "mean", "median", "none"}, optional) – Preprocessing mode applied independently to each event window.
strict_nyquist (bool, optional) – Whether oscillator frequencies at or above Nyquist should raise an error instead of being clipped.
peak_mode ({"abs", "pos", "neg"}, optional) – Peak convention used for each event-local spectrum.

Returns:

Rolling result returned as a RollingERSResult. Each row corresponds to one detected event and each column corresponds to one oscillator frequency.

Return type:

Notes

Rolling shock spectra are currently event-window based rather than fixed-hop or sliding-window based. The time coordinate stored in the result corresponds to the detected peak time of each event.

If no events are present, the function still returns a valid empty matrix with the validated oscillator frequency grid.

fdscore.envelope_srs(results)[source]¶

Compute an envelope across compatible shock-response spectra.

Parameters:: results (collections.abc.Sequence) – Sequence containing either one-sided ERSResult values or two-sided ShockSpectrumPair values. Mixing the two forms in a single call is not allowed.
Returns:: Enveloped one-sided or two-sided shock-response spectrum, matching the representation of the input sequence.
Return type:: object

Notes

This function accepts only spectra whose ERS compatibility metadata identifies them as "shock_response_spectrum". For pair-valued inputs, the negative and positive sides are enveloped independently.

fdscore.envelope_pvss(results)[source]¶

Compute an envelope across compatible pseudo-velocity shock spectra.

Parameters:: results (collections.abc.Sequence) – Sequence containing either one-sided ERSResult values or two-sided ShockSpectrumPair values. Mixing the two forms in a single call is not allowed.
Returns:: Enveloped one-sided or two-sided PVSS result, matching the representation of the input sequence.
Return type:: object

Notes

This function accepts only spectra whose ERS compatibility metadata identifies them as "pseudo_velocity_shock_spectrum". For pair-valued inputs, the negative and positive sides are enveloped independently.

fdscore.fit_half_sine_to_pvss(pvss, *, polarity='pos')[source]¶

Fit a half-sine acceleration pulse whose PVSS envelopes the input PVSS.

Parameters:

pvss (fdscore.types.ERSResult) – One-sided pseudo-velocity shock spectrum compatible with the library’s PVSS contract.
polarity ({"pos", "neg"}, optional) – Polarity assigned to the fitted half-sine pulse.

Returns:

Fitted pulse returned as a HalfSinePulse, expressed in acceleration amplitude and pulse duration.

Return type: