API

Normalized two-sided FFT spectrum of an AccelerationTrace.

spectrum = AccelerationSpectrum(trace)

computes the FFT of trace.acceleration, normalizes so that ∫|F(f)|²df = rms², reorders into a two-sided layout, and holds the result. Requires trace.tlist to be uniform. See write_acceleration_spectrum and check_acceleration_spectrum.

Fields

• freq_axis: frequency axis of two-sided spectrum (Hz), length n_samples + 1; uniform with spacing df = freqaxis[2] − freqaxis[1]
• spectrum: complex spectral amplitude (m/s²/√Hz), length n_samples + 1
• rms: RMS computed from the spectrum as √(∑|F(f)|²·df) (m/s²); should match the source trace RMS — checked by check_acceleration_spectrum

Time-domain acceleration waveform generated by superimposing random sinusoids.

trace = AccelerationTrace(
    psd;
    tlist, f_band_min, f_band_max, n_oscillators, rng, target_rms,
    show_progress)

superimposes n_oscillators random sinusoids with physically derived amplitudes Aₖ = √(2·S(fₖ)·Δf/M) and holds the resulting time-domain waveform and oscillator parameters. When target_rms is given, all arrays are rescaled by eta = targetrms/actualrms to enforce the RMS exactly. Use AccelerationSpectrum to compute the FFT spectrum; write with write_acceleration_trace.

Keyword Arguments

• tlist: time grid (s); typically a uniform range or Vector{Float64}, but any AbstractVector is accepted — oscillators are evaluated at each tlist[it], so non-uniform grids work. For AccelerationSpectrum the grid must be uniform (asserted there).
• f_band_min=psd.f_min, f_band_max=psd.f_max: frequency band for oscillator sampling (Hz). Must be in the frequency range supported by psd.
• n_oscillators: number of random sinusoids to superimpose
• rng: random number generator; any AbstractRNG; defaults to Random.default_rng() (non-fixed seed). Pass a Ran2State for legacy reproducible output.
• target_rms: desired RMS acceleration (m/s²); if nothing (default), no rescaling is applied and the raw Monte Carlo RMS is kept.
• show_progress: if true, display a progress bar while drawing oscillators; default true

Fields

• tlist: time grid (s), length n_samples; typically uniform but may be arbitrary — AccelerationSpectrum asserts uniformity
• acceleration: acceleration waveform (m/s²), length n_samples
• velocity: velocity waveform (m/s), length n_samples; satisfies $v̇ = a$ entry-by-entry
• displacement: displacement waveform (m), length n_samples; satisfies $ẍ = a$ entry-by-entry
• oscillators: 3 × n_oscillators matrix of (amplitude Aₖ (m/s²), frequency fₖ (Hz), phase φₖ (rad))
• f_band_min, f_band_max: noise band limits (Hz)
• rms: actual RMS of acceleration (m/s²)
• eta: rescaling factor applied to all arrays; 1.0 if target_rms=nothing, otherwise targetrms / rawrms

Acceleration PSD S(f) loaded from a two-column ASCII file.

psd = PSDData(file; acceleration_unit=1.0, interpolation_method=:linear)

reads the power spectral density (PSD) from file. The file must contain two whitespace-separated columns; lines starting with # are ignored as comments:

• Column 1: frequency (Hz), must be sorted in ascending order
• Column 2: one-sided acceleration PSD in a²/Hz

Calling PSDData() without positional arguments instead prompts for the file path and acceleration_unit interactively from stdin.

The acceleration_unit is the conversion factor from the file's acceleration unit $a$ to m/s²: acceleration_unit=9.81 if the PSD is in g²/Hz; acceleration_unit=1 if already in (m/s²)²/Hz.

By the Wiener–Khinchin theorem, √∫ S(f) df is the RMS √⟨a²⟩ of a stationary random process with one-sided PSD S(f), so psd.rms is the physically expected RMS acceleration.

The resulting psd object is a callable psd(f) → S(f) in (m/s²)²/Hz, regardless of the file's unit, returning interpolated values for the PSD.

A PSDData object can be plotted directly with any Plots-compatible backend, showing S(f) in (m/s²)²/Hz vs frequency in Hz.

Keyword Arguments

• acceleration_unit: conversion factor from file units to m/s²; default 1.0
• interpolation_method: :linear (default, reproduces trapezoid-rule results exactly) or :pchip (monotone-preserving cubic, smoother S(f) between nodes)

Fields

• f_min: lowest frequency in the PSD file (Hz)
• f_max: highest frequency in the PSD file (Hz)
• n_samples: number of data points in the PSD file
• rms: $√∫S(f)df$ in (m/s²), the expected root-mean-square of the acceleration trace described by the PSD.

State for the ran2! long-period (> 2×10¹⁸) random-number generator.

rng = Ran2State(-76989321)

initialises the generator from Numerical Recipes in Fortran, 2nd ed., §7.1 (© 1986–92 Numerical Recipes Software), translated faithfully with Int32 arithmetic. Pass a negative seed; advance with ran2! or rand(rng).

Implements Random.AbstractRNG so it can be passed anywhere a Julia RNG is expected.

Keyword Arguments

• verbose: if true, print "Initialize RAN2" on the first draw; default false

Check that the spectrum RMS is consistent with the source acceleration trace.

rms_err = check_acceleration_spectrum(trace, spectrum; verbose=true)

returns |spectrum.rms − trace.rms| / trace.rms, the relative difference between the RMS computed independently from the FFT (√(∑|F(f)|²·df)) and the time-domain RMS stored in trace. By the normalization convention this should be zero up to floating-point precision (typically < 1e-14); larger values indicate a bug in the spectrum normalization.

Keyword Arguments

• verbose: if true, print both RMS values and the relative error; default false

Check that an AccelerationTrace RMS is consistent with Monte-Carlo statistics.

n_sigma = check_acceleration_trace(psd, trace; verbose=true)

returns the deviation of the oscillator-amplitude-based RMS from the band-limited PSD RMS, expressed in units of the Monte-Carlo standard deviation σ. The band limits and oscillator count are read directly from trace. Values ≤ 4 are expected ~99.994% of the time for correct generation; larger values suggest a bug. Whether target_rms was given during generation does not matter: the raw RMS is recovered from the stored oscillator amplitudes.

Keyword Arguments

• verbose: if true, print expected/raw RMS and σ summary; default false

Monte-Carlo convergence

With M = n_oscillators sinusoids with amplitudes Aₖ = √(2 S(fₖ) Δf/M) and frequencies fₖ ~ Uniform[fbandmin, fbandmax] (Δf = fbandmax − fbandmin), the raw squared RMS is

raw_rms² = Σₖ Aₖ²/2 = (Δf/M) Σₖ S(fₖ)

— Δf times the sample mean of S(fₖ). By the law of large numbers this converges to

E[raw_rms²] = Δf · E[S] = ∫_band S(f) df = expected_rms²

where expected_rms = rms(psd; f_min=f_band_min, f_max=f_band_max).

Note: raw_rms here is computed from the oscillator amplitudes, not from the time-domain waveform. The time-domain RMS equals √(Σₖ Aₖ²/2) only when every oscillator completes many cycles ((1/N) Σₙ cos²(ωₖtₙ + φₖ) → 1/2). For oscillators with fₖ < df = 1/(N·dt), this time average depends heavily on the random phase and can deviate far from 1/2. Using the amplitude-based formula keeps raw_rms a pure MC quantity so that σ_rms below gives the correct z-score regardless of window length.

The M PSD samples S(fₖ) = Aₖ²·M/(2·Δf) are recovered from the stored oscillator amplitudes, and their sample variance estimates Var[S(f)] directly. The variance of the sample mean then gives

Var[raw_rms²] ≈ (Δf²/M) · var(S(fₖ))

and first-order error propagation (δ√x ≈ δx / 2√x) yields

σ_rms ≈ Δf · √(var(S(fₖ)) / M) / (2 · expected_rms)

Run the interactive noise generation program, reading all inputs from stdin.

main()

queries for sampling parameters, a PSD file, oscillator count, random seed, target RMS, and output file paths, then writes a time trace and spectrum file.

Draw the next uniform random Float64 in (0, 1) from a Ran2State.

x = ran2!(rng)

advances rng in-place and returns the next value.

Read one stdin line and parse n values of the given type.

(f_min, f_max) = read_stdin(Float64, 2)
x              = read_stdin(Int)

strips any label up to and including the first colon, so "tstep (s): 1e-4" and "1e-4" parse identically. Float64 inputs accept Fortran syntax ("1d-4" → 1e-4). Returns a Tuple for n > 1, a scalar for n = 1.

Resample an AccelerationTrace at an arbitrary time grid.

resampled = resample_acceleration_trace(trace; tlist=0:1e-3:10.0)

re-evaluates the oscillator superposition from trace.oscillators at each point in tlist and returns a new AccelerationTrace with the same oscillator parameters (eta = 1.0, rms is the actual RMS of the resampled waveform). Write the result with write_acceleration_trace.

Keyword Arguments

• tlist: new time grid (s); may be uniform or non-uniform
• show_progress: if true, display a progress bar while evaluating time points; default true

RMS acceleration from a PSDData, AccelerationTrace, or AccelerationSpectrum.

band_rms  = rms(psd; f_min=20.0, f_max=1000.0)
trace_rms = rms(trace)
spec_rms  = rms(spectrum; f_min=20.0, f_max=1000.0)

For PSDData, integrates S(f) over [fmin, fmax] using the interpolant; without windowing (f_min=0, f_max=Inf), returns psd.rms. For AccelerationTrace, returns trace.rms directly (no keyword arguments). For AccelerationSpectrum, sums |F(f)|²·df over all f with |f| ∈ [fmin, fmax]; without windowing, recovers spectrum.rms.

Comparing rms(psd; f_min=f_band_min, f_max=f_band_max) to rms(trace) when target_rms=nothing isolates Monte-Carlo fluctuations from the systematic band-limitation bias.

Keyword Arguments

• f_min: lower frequency bound (Hz); default 0.0
• f_max: upper frequency bound (Hz); default Inf

Write an AccelerationSpectrum to a spectrum ASCII file.

write_acceleration_spectrum(spectrum, "spectrum.dat")

writes four columns: frequency (Hz), Re (m/s²/√Hz), Im (m/s²/√Hz), two-sided PSD ((m/s²)²/Hz).

Write an AccelerationTrace to a time-trace ASCII file.

write_acceleration_trace(trace, "time_trace.dat")

writes four columns: time (s), acceleration (m/s²), velocity (m/s), displacement (m). The time column is taken directly from trace.tlist.