// SPDX-FileCopyrightText: © 2026 Jeffrey C. Ollie <jeff@ocjtech.us>
// SPDX-License-Identifier: MIT

//! An endless WOPR machine room hum.
//!
//! The model is additive and deliberately small: one sine oscillator per
//! measured partial (see `partials`), two noise generators for the floor
//! beneath them, and a slow random wander applied to the frequency and
//! amplitude of each oscillator so that the sound never settles into a loop.
//! There is no filter bank, no reverb and no tremolo, because the recording
//! does not call for any:
//!
//! * The throb is not a tremolo. The loudest partials sit about 3.7 Hz apart,
//!   so summing them beats at 3.7 Hz on its own. Modulating the sum as well
//!   would be doing the same job twice, audibly.
//! * The hum is not a filtered sawtooth, or anything else with a fundamental.
//!   Its partials are not harmonically related, so nothing subtractive
//!   reaches the same spectrum. `partials.zig` says more about why.
//!
//! What the wander buys is the absence of a period. Twenty-eight fixed sines
//! repeat exactly whenever their frequencies happen to share one, and the ear
//! finds that quickly and stops hearing a machine. Letting the whole bank
//! drift together by a fraction of a percent — the way an unregulated motor
//! does, and the way the recording does, whose strongest partial wanders
//! between 123.6 Hz and 125.2 Hz across twelve seconds — removes the period
//! without moving the pitch anywhere the ear can follow.
//!
//! ```zig
//! var hum: Hum = try .init(gpa, .{});
//! defer hum.deinit(gpa);
//!
//! var frames: [2048]f32 = undefined; // interleaved, 1024 stereo frames
//! while (true) {
//!     hum.render(&frames);
//!     try out.writeAll(std.mem.sliceAsBytes(&frames));
//! }
//! ```
//!
//! `render` is the only thing that runs per sample, it allocates nothing and
//! it never fails, so it is safe to call from an audio callback. `init` is
//! the only part that allocates, and it allocates once.

const Hum = @This();

const std = @import("std");
const Allocator = std.mem.Allocator;

const partials = @import("partials.zig");
const burst_table = @import("bursts.zig");

/// The block size at which the drift and shimmer are recomputed, in frames.
///
/// Wandering by a fraction of a percent every 1.5 ms is indistinguishable
/// from wandering every sample, and costs 1/64th as much; amplitudes are
/// ramped across the block rather than stepped, so nothing zippers.
const control_period = 64;

pub const Options = struct {
    /// Frames per second.
    sample_rate: u32 = 44_100,

    /// 1 for mono, 2 for stereo. Anything else is `error.UnsupportedChannelCount`.
    channels: u8 = 2,

    /// Chooses the random phases, pans and wander. The same seed and options
    /// give byte-identical output, which is what makes the tests possible;
    /// different seeds give a different instance of the same machine.
    seed: u64 = 0,

    /// The partials to sum. Defaults to the ones measured off the recording.
    partials: []const Partial = partials.wopr,

    /// Output level, in dBFS RMS. The default leaves about 8 dB of headroom
    /// above the peaks of a signal whose crest factor is around 10 dB; the
    /// reference recording sits at roughly -11 dBFS RMS at its loudest, which
    /// is louder than anything should be asked to render into a fixed-point
    /// format without a limiter in the way.
    level_dbfs: f64 = -18.0,

    /// The broadband floor under the hum, in dB relative to the tonal part,
    /// or `null` for no floor at all. `null` is the default.
    ///
    /// The recording has such a floor and this can reproduce it, but off is
    /// the better default and the reason is not acoustics. A room already
    /// has a floor -- fans, traffic, the building -- and a synthesiser that
    /// runs for hours as a background bed is heard *in* that room, not
    /// instead of it. Adding a second one only ever puts a steady, smooth
    /// hiss on top of the listener's own, which is the one thing here that
    /// people notice and dislike; the partials do not need it to sound like
    /// a machine.
    ///
    /// Turn it on with a level if the hum is going somewhere that has no
    /// floor of its own -- under dialogue in a dry mix, say. -28 dB is
    /// present without being audible as hiss; -16 dB matches the recording.
    hiss_level_db: ?f64 = null,

    /// The shape of that floor, when there is one.
    ///
    /// Fitted to the recording's own floor between 50 Hz and 800 Hz, which
    /// it follows to within 2 dB, and held well under it above 1 kHz.
    /// `BandShape` says why the fit is split that way.
    hiss: BandShape = .{ .high_hz = 160, .tilt_hz = 200, .cutoff_hz = 1000 },

    /// Rumble, in dB relative to the tonal part: the room and the air
    /// handling under everything else.
    ///
    /// Quiet, and deliberately so. The recording has almost nothing below
    /// 80 Hz -- its 20-40 Hz octave is 46 dB under the hum, which is to say
    /// absent -- so a synthesiser with a generous bottom end would be a
    /// different and less convincing sound however good it felt on a
    /// subwoofer.
    rumble_level_db: f64 = -42,

    /// The shape of the rumble: a band rather than everything below a
    /// cutoff, because the recording falls away below 40 Hz as steeply as it
    /// does above 80.
    rumble: BandShape = .{ .high_hz = 42, .tilt_hz = 78, .cutoff_hz = 300 },

    /// How far the whole bank drifts in frequency, as an RMS fraction. All
    /// the partials move together by this much, because they are one machine
    /// changing speed. 0.004 is ±0.4%, which matches the wander measured in
    /// the recording.
    drift: f64 = 0.004,

    /// How long the bank takes to wander, in seconds. Longer is slower and
    /// less noticeable.
    drift_seconds: f64 = 3.0,

    /// How far each partial drifts independently of the others, as an RMS
    /// fraction. Small, but it is what stops the beat pattern from repeating:
    /// the partials' spacing changes, so the 3.7 Hz throb breathes.
    voice_drift: f64 = 0.0008,

    /// How far each partial's level wanders, in dB RMS.
    shimmer_db: f64 = 2.0,

    /// How long a partial takes to wander in level, in seconds.
    shimmer_seconds: f64 = 4.0,

    /// The tone bursts heard over the hum -- the ping and the pong. Empty
    /// for a hum with nothing happening over it.
    bursts: []const Burst = burst_table.wopr,

    /// How often those arrive, and in what pattern.
    timing: Timing = .{},

    /// The room the bursts echo in, or `null` for no echo.
    echo: ?EchoShape = .{},

    /// Stereo spread, 0 to 1. Ignored when `channels` is 1.
    ///
    /// 0 puts the same signal in both channels. 1 is a pair of microphones
    /// `spacing_seconds` apart in a room; the default reproduces the
    /// reference recording's channel correlation of 0.97, which is a wide
    /// but overwhelmingly centred image.
    ///
    /// Correlation rather than the recording's side-to-mid ratio, which is
    /// 1.6 dB of plain level imbalance between its two channels on top of
    /// the image width. That is a property of how the recording was made,
    /// not of the room, and matching it would mean shipping a lopsided mix.
    width: f64 = 0.6,

    /// How far apart the two microphones are, expressed as the time it takes
    /// sound to travel between them. Scaled by `width`.
    ///
    /// This is the whole of the decorrelation that matters, and it is a
    /// delay rather than a random phase per partial for a reason worth
    /// recording: a random phase makes the channel correlation depend on
    /// where the *strongest* partial's draw happened to land, and one
    /// partial here is 8 dB above every other, so the image came out
    /// anywhere between 0.95 and 0.997 depending on the seed. A delay gives
    /// every partial a phase difference of `2*pi*f*tau`, which is what a
    /// spaced pair actually does, and makes the width a property of the
    /// options rather than of the seed.
    ///
    /// 450 microseconds is about 15 cm of air.
    spacing_seconds: f64 = 450e-6,
};

pub const Partial = partials.Partial;
pub const Burst = burst_table.Burst;
pub const Timing = burst_table.Timing;

/// The room the bursts echo in.
///
/// One delay line per channel with a little feedback, which is a poor
/// reverb and exactly the right amount of one: what the recording has is a
/// discrete repeat about 210 ms after each burst, not a tail. Its envelope
/// autocorrelation peaks at 200 to 215 ms.
///
/// The hum does not go through it, and not for want of realism. A
/// continuous sound convolved with its own echo is the same continuous
/// sound very slightly thicker; there is nothing to hear, and it would cost
/// a second delay line long enough to matter.
pub const EchoShape = struct {
    /// Delay for each channel, in seconds.
    ///
    /// Different per channel on purpose. A room is not symmetric, and two
    /// equal delays put the repeat in the middle of the head where the dry
    /// burst already is; a few tens of milliseconds between them is what
    /// makes it sound like a space rather than a copy.
    delay_s: [2]f64 = .{ 0.205, 0.232 },

    /// How much of each repeat feeds the next.
    ///
    /// Low, so that a burst gets *one* audible repeat rather than a train
    /// of them. `level_db` sets how loud the echo is and this sets how many
    /// there are, and the two are easy to confuse: the first version had
    /// 0.42 here, which put a third and fourth repeat 14 and 21 dB down --
    /// inaudible in the recording, where the noise floor covers them, and
    /// plainly audible here, where by default there is no floor. It made
    /// the whole thing sound like it was running fast when the burst rate
    /// was right.
    feedback: f64 = 0.15,

    /// The repeats lose their top end, as they would off real surfaces.
    /// A one-pole lowpass in the feedback path, at this frequency.
    damping_hz: f64 = 2600.0,

    /// Level of the first repeat, in dB relative to the dry burst.
    level_db: f64 = -7.0,
};

/// How one band of the noise floor is shaped.
///
/// Three stages, because the floor has to do two different things. Through
/// the mid-band it has to *tilt*, following the recording's own slope, and a
/// pair of one-poles does that. Above the band it has to actually stop, and
/// a cascade of one-poles does not: a one-pole digital lowpass flattens out
/// towards Nyquist at `(1-a)/(1+a)`, so a corner low enough to shape this
/// floor leaves a broadband tail only about 30 dB down per pole. Two
/// Butterworth biquads have a double zero at Nyquist and keep falling,
/// which is what lets the band actually end.
///
/// None of this is a detail, and two rounds of listening went into the
/// numbers. A version shaped with two one-poles alone matched the
/// recording's octave band energies to within 2 dB everywhere below 1 kHz
/// and hissed audibly at 4 to 8 kHz -- where the ear is at its most
/// sensitive, and where a band total, dominated by the octave below it,
/// showed nothing wrong. Adding the biquads at 2.6 kHz fixed that and left
/// a quieter hiss still, this time because a *smooth, steady* floor at
/// 1 to 3 kHz is heard as hiss whatever its level: the recording's content
/// up there has a spectral flatness of 0.11 to 0.33 at its peaks and swings
/// 27 dB from frame to frame, because it is discrete and intermittent
/// rather than a floor at all. Matching its energy there with smooth noise
/// matches the measurement and not the sound.
///
/// So the corners below are fitted to the recording between 200 Hz and
/// 800 Hz, where it really does have a floor and where the floor's audible
/// job is -- filling in between the partials so the hum sounds like a room
/// -- and constrained to stay well under the recording above 1 kHz rather
/// than to match it.
pub const BandShape = struct {
    /// Two one-pole highpasses: the floor stops below the hum rather than
    /// continuing down to DC, and stops as steeply as the recording does.
    high_hz: f64,
    /// Two one-pole lowpasses, which set the slope through the mid-band.
    tilt_hz: f64,
    /// Two Butterworth biquads, which end the band.
    cutoff_hz: f64,
};

pub const InitError = error{
    UnsupportedChannelCount,
    NoPartials,
    /// A partial at or above half the sample rate has nowhere to be.
    PartialAboveNyquist,
    /// A burst with more than `burst_table.max_partials` of them, or none.
    BadBurst,
    /// An echo delay that is not a positive number of seconds, or is longer
    /// than `max_echo_seconds`.
    BadEchoDelay,
} || Allocator.Error;

/// The longest echo that will be allocated, which bounds what a bad option
/// can ask the allocator for.
pub const max_echo_seconds = 4.0;

/// How many bursts may sound at once. The recording's flurries put them
/// 80 ms apart and each lasts about 100 ms, so two or three overlap; this
/// is enough that a burst is never cut short to make room for the next.
const max_burst_voices = 6;

voices: []Voice,
/// The shared floor, then one private to each channel. Three rather than
/// two because a filter has state: stepping one generator twice per frame to
/// get two channels runs its poles at twice the sample rate they were
/// designed for, which moves every cutoff up an octave and makes the hiss
/// audibly brighter than it was tuned to be.
noise: [3]NoiseFloor,
drift: Wander,
rng: std.Random.DefaultPrng,
sample_rate: f64,
channels: u8,
/// Applied to the whole mix, folding in the RMS normalisation and `level_dbfs`.
gain: f64,
hiss_gain: f64,
rumble_gain: f64,
/// How much of each channel's noise is its own rather than shared.
noise_spread: f64,
/// Frames left before the control values are recomputed.
countdown: u32 = 0,

// -- the bursts over the hum -------------------------------------------------

/// What may be fired. Empty means nothing is.
bursts: []const Burst,
timing: Timing,
/// `timing`'s intervals, already multiplied by `Timing.scale` and the
/// sample rate, so firing one costs no division.
gap_frames: [2]f64,
rest_frames: [2]f64,
/// The bursts draw from a generator of their own, never from `rng`.
///
/// Not tidiness: `render` runs the hum's control rate and the burst
/// scheduler over the same buffer, and a caller is free to ask for the
/// frames in any sized pieces it likes. Sharing one generator would
/// interleave the two streams of draws differently for a 64 frame call than
/// for a 4096 frame one, and the same seed would stop meaning the same
/// audio.
burst_rng: std.Random.DefaultPrng,
sounding: [max_burst_voices]BurstVoice = @splat(.{}),
/// Frames until the next burst starts.
until_burst: u64 = 0,
/// Bursts left in the flurry being played, and which kind it is made of.
flurry_left: u32 = 0,
flurry_kind: usize = 0,
/// Peak amplitude for each burst kind, already including the output level.
burst_amp: []f64,
echo: ?Echo,

pub fn init(gpa: Allocator, options: Options) InitError!Hum {
    if (options.channels != 1 and options.channels != 2) return error.UnsupportedChannelCount;
    if (options.partials.len == 0) return error.NoPartials;

    const rate: f64 = @floatFromInt(options.sample_rate);
    const nyquist = rate / 2;
    // Checked against the drifted extreme rather than the nominal frequency,
    // so that a partial cannot alias only sometimes, minutes into a run.
    const headroom = 1 + 4 * (options.drift + options.voice_drift);
    for (options.partials) |p| {
        if (p.freq <= 0 or p.freq * headroom >= nyquist) return error.PartialAboveNyquist;
    }

    for (options.bursts) |b| {
        if (b.partials.len == 0 or b.partials.len > burst_table.max_partials) return error.BadBurst;
        if (!(b.weight > 0)) return error.BadBurst;
        for (b.partials) |p| if (p.freq <= 0 or p.freq >= nyquist) return error.BadBurst;
    }
    if (options.echo) |e| for (e.delay_s) |d| {
        if (!(d > 0) or d > max_echo_seconds) return error.BadEchoDelay;
    };

    var prng: std.Random.DefaultPrng = .init(options.seed);
    const rng = prng.random();

    const control_rate = rate / control_period;
    const width = std.math.clamp(options.width, 0, 1);

    const voices = try gpa.alloc(Voice, options.partials.len);
    errdefer gpa.free(voices);

    const burst_amp = try gpa.alloc(f64, options.bursts.len);
    errdefer gpa.free(burst_amp);

    var echo: ?Echo = if (options.echo) |shape|
        try Echo.init(gpa, shape, rate, options.channels)
    else
        null;
    errdefer if (echo) |*e| e.deinit(gpa);

    var power: f64 = 0;
    for (voices, options.partials) |*v, p| {
        power += p.amp * p.amp;

        // A pan either side of centre, for texture, and the phase
        // difference the microphone spacing gives this partial, which is
        // what actually sets the width.
        const placed = pan(width * 0.24, rng);
        const skew = p.freq * width * options.spacing_seconds;

        v.* = .{
            .freq = p.freq,
            .amp = p.amp,
            .phase = .{ rng.float(f64), 0 },
            .gain = placed,
            .freq_wander = .init(options.drift_seconds, control_rate, options.voice_drift),
            .amp_wander = .init(options.shimmer_seconds, control_rate, options.shimmer_db),
            .inc = 0,
            .level = p.amp,
            .level_step = 0,
        };
        v.phase[1] = wrap(v.phase[0] + skew);
    }

    // The bank is normalised to unit RMS before the output level is applied,
    // so that changing the partial table -- or handing in a different one --
    // changes the timbre and not the loudness. Sines contribute half their
    // squared amplitude to the mean square.
    // The pans are random, so nothing makes the two channels come out at the
    // same level on their own -- and with one partial 8 dB above the rest,
    // where that one happens to sit decides the balance of the whole mix.
    // Measured over a dozen seeds it strays past a decibel, which is an
    // audibly lopsided image and one that changes every run. Equalising the
    // two sides here keeps the width the pans bought and throws the
    // imbalance away.
    var side_power: [2]f64 = .{ 0, 0 };
    for (voices) |v| for (0..2) |ch| {
        side_power[ch] += v.amp * v.amp * v.gain[ch] * v.gain[ch];
    };
    const balanced = (side_power[0] + side_power[1]) / 2;
    for (voices) |*v| for (0..2) |ch| {
        if (side_power[ch] > 0) v.gain[ch] *= @sqrt(balanced / side_power[ch]);
    };

    const tonal_rms = @sqrt(power / 2);
    const level = std.math.pow(f64, 10, options.level_dbfs / 20);

    // The shimmer is symmetric in decibels, which is not symmetric in power:
    // a partial spends as long 3.5 dB up as 3.5 dB down, and 3.5 dB up is the
    // larger change. Left alone that makes the output louder than
    // `level_dbfs` asked for, by e^(2 c^2 sigma^2) in power for c = ln(10)/20,
    // so it is taken back out here rather than left as a surprise that grows
    // with the shimmer setting.
    const c = @log(@as(f64, 10)) / 20;
    const shimmer_bias = @exp(-c * c * options.shimmer_db * options.shimmer_db);

    // The bursts are specified relative to the hum, and `gain` has already
    // been folded into the hum, so their amplitudes carry the output level
    // directly and are added after it.
    for (burst_amp, options.bursts) |*a, b| {
        a.* = level * std.math.pow(f64, 10, b.level_db / 20);
    }

    const interval_scale = options.timing.scale() * rate;

    return .{
        .voices = voices,
        .bursts = options.bursts,
        .timing = options.timing,
        .gap_frames = .{
            options.timing.gap_s[0] * interval_scale,
            options.timing.gap_s[1] * interval_scale,
        },
        .rest_frames = .{
            options.timing.rest_s[0] * interval_scale,
            options.timing.rest_s[1] * interval_scale,
        },
        // Somewhere inside the first rest, so that two instances started
        // together do not ping in unison.
        .until_burst = @intFromFloat(rng.float(f64) * options.timing.rest_s[1] * interval_scale),
        .burst_rng = .init(rng.int(u64)),
        .burst_amp = burst_amp,
        .echo = echo,
        .noise = .{
            .init(rng.int(u64), options, rate),
            .init(rng.int(u64), options, rate),
            .init(rng.int(u64), options, rate),
        },
        .drift = .init(options.drift_seconds, control_rate, options.drift),
        .rng = prng,
        .sample_rate = rate,
        .channels = options.channels,
        .gain = level * shimmer_bias / tonal_rms,
        .hiss_gain = if (options.hiss_level_db) |db| tonal_rms * std.math.pow(f64, 10, db / 20) else 0,
        .rumble_gain = tonal_rms * std.math.pow(f64, 10, options.rumble_level_db / 20),
        .noise_spread = if (options.channels == 1) 0 else width,
    };
}

pub fn deinit(hum: *Hum, gpa: Allocator) void {
    gpa.free(hum.voices);
    gpa.free(hum.burst_amp);
    if (hum.echo) |*e| e.deinit(gpa);
    hum.* = undefined;
}

/// Fill `out` with the next samples, interleaved by channel.
///
/// `out.len` must be a multiple of the channel count. Nothing here allocates,
/// blocks or can fail, so this is safe to call from an audio callback.
pub fn render(hum: *Hum, out: []f32) void {
    std.debug.assert(out.len % hum.channels == 0);
    const rng = hum.rng.random();

    var i: usize = 0;
    while (i < out.len) {
        if (hum.countdown == 0) {
            hum.retune(rng);
            hum.countdown = control_period;
        }
        const frames_left = (out.len - i) / hum.channels;
        const run = @min(hum.countdown, frames_left);
        hum.countdown -= @intCast(run);

        if (hum.channels == 1) hum.renderMono(out[i..][0..run]) else hum.renderStereo(out[i..][0 .. run * 2]);
        i += run * hum.channels;
    }

    // Added over the finished hum rather than mixed into it, because their
    // levels are relative to it and because they go through the echo and it
    // does not.
    hum.renderBursts(out);
}

/// Recompute the per-partial frequency and level. Called once per
/// `control_period` frames.
fn retune(hum: *Hum, rng: std.Random) void {
    const speed = 1 + hum.drift.next(rng);
    const inv_rate = 1 / hum.sample_rate;
    const inv_period: f64 = 1.0 / @as(f64, control_period);
    for (hum.voices) |*v| {
        v.inc = v.freq * speed * (1 + v.freq_wander.next(rng)) * inv_rate;
        const target = v.amp * std.math.pow(f64, 10, v.amp_wander.next(rng) / 20);
        v.level_step = (target - v.level) * inv_period;
    }
}

fn renderMono(hum: *Hum, out: []f32) void {
    @memset(out, 0);
    for (hum.voices) |*v| {
        var level = v.level;
        var phase = v.phase[0];
        for (out) |*s| {
            s.* += @floatCast(level * sine(phase));
            phase = wrap(phase + v.inc);
            level += v.level_step;
        }
        v.phase[0] = phase;
        v.level = level;
    }
    const floor = &hum.noise[0];
    for (out) |*s| {
        const n = floor.next(hum.hiss_gain, hum.rumble_gain);
        s.* = @floatCast((s.* + n) * hum.gain);
    }
}

fn renderStereo(hum: *Hum, out: []f32) void {
    @memset(out, 0);
    for (hum.voices) |*v| {
        for (0..2) |ch| {
            var level = v.level;
            var phase = v.phase[ch];
            const g = v.gain[ch];
            var k: usize = ch;
            while (k < out.len) : (k += 2) {
                out[k] += @floatCast(g * level * sine(phase));
                phase = wrap(phase + v.inc);
                level += v.level_step;
            }
            v.phase[ch] = phase;
            // Both channels walked the same ramp; keep the one they agree on.
            if (ch == 1) v.level = level;
        }
    }

    // A shared noise source plus a private one per channel. `noise_spread`
    // decides how much of each channel is its own, which is what keeps the
    // floor from either collapsing to the centre or throwing the channel
    // correlation away entirely.
    const common = @sqrt(1 - hum.noise_spread * hum.noise_spread);
    var k: usize = 0;
    while (k < out.len) : (k += 2) {
        const shared = hum.noise[0].next(hum.hiss_gain, hum.rumble_gain);
        const l = hum.noise[1].next(hum.hiss_gain, hum.rumble_gain);
        const r = hum.noise[2].next(hum.hiss_gain, hum.rumble_gain);
        out[k] = @floatCast((out[k] + common * shared + hum.noise_spread * l) * hum.gain);
        out[k + 1] = @floatCast((out[k + 1] + common * shared + hum.noise_spread * r) * hum.gain);
    }
}

/// Render the bursts, and their echo, over what is already in `out`.
fn renderBursts(hum: *Hum, out: []f32) void {
    if (hum.bursts.len == 0 and hum.echo == null) return;
    const rng = hum.burst_rng.random();
    const channels = hum.channels;
    const frames = out.len / channels;

    for (0..frames) |i| {
        if (hum.until_burst == 0) hum.fire(rng);
        hum.until_burst -= 1;

        var dry: [2]f64 = .{ 0, 0 };
        for (&hum.sounding) |*v| v.add(&dry, channels);

        if (hum.echo) |*e| e.process(&dry, channels);

        for (0..channels) |ch| {
            const k = i * channels + ch;
            out[k] += @floatCast(dry[ch]);
        }
    }
}

/// Start a burst, and decide when the next one is.
fn fire(hum: *Hum, rng: std.Random) void {
    if (hum.bursts.len == 0) {
        // Nothing to fire, but the countdown still has to go somewhere
        // other than zero or this is called every frame.
        hum.until_burst = @intFromFloat(@max(1, hum.rest_frames[1]));
        return;
    }

    if (hum.flurry_left == 0) {
        hum.flurry_kind = hum.pickBurst(rng);
        hum.flurry_left = geometric(rng, hum.timing.flurry);
    }

    // Nearly always the flurry's own kind; occasionally the other, which is
    // what the recording does at 7.34 s where a ping and a pong land within
    // three milliseconds of each other.
    const kind = if (hum.bursts.len > 1 and rng.float(f64) < hum.timing.stray)
        hum.pickBurst(rng)
    else
        hum.flurry_kind;
    hum.start(kind, rng);

    hum.flurry_left -= 1;
    const range = if (hum.flurry_left == 0) hum.rest_frames else hum.gap_frames;
    const span = range[0] + rng.float(f64) * (range[1] - range[0]);
    hum.until_burst = @intFromFloat(@max(1, span));
}

/// Choose a burst kind, weighted by `Burst.weight`.
fn pickBurst(hum: *Hum, rng: std.Random) usize {
    var total: f64 = 0;
    for (hum.bursts) |b| total += b.weight;
    var pick = rng.float(f64) * total;
    for (hum.bursts, 0..) |b, i| {
        pick -= b.weight;
        if (pick <= 0) return i;
    }
    return hum.bursts.len - 1;
}

/// Hand a burst to a voice.
fn start(hum: *Hum, kind: usize, rng: std.Random) void {
    const burst = &hum.bursts[kind];

    // The oldest sounding voice if none is free. Stealing is better than
    // dropping: a flurry that runs past the polyphony should crowd, which
    // is what a room does, rather than go quiet.
    var slot: *BurstVoice = &hum.sounding[0];
    var oldest: u32 = 0;
    for (&hum.sounding) |*v| {
        if (v.burst == null) {
            slot = v;
            break;
        }
        if (v.frame >= oldest) {
            oldest = v.frame;
            slot = v;
        }
    }

    // The recording's fundamentals vary by about half a percent from one
    // occurrence to the next, so these do too: identical pitch every time
    // is the sound of a sample being retriggered.
    const detune = 1 + 0.004 * rng.floatNorm(f64);
    const inv_rate = 1 / hum.sample_rate;

    slot.* = .{
        .burst = burst,
        .amp = hum.burst_amp[kind],
        .attack = @max(1, @as(u32, @intFromFloat(burst.attack_s * hum.sample_rate))),
        .hold = @intFromFloat(@max(0, burst.hold_s *
            (1 + burst.hold_spread * (rng.float(f64) * 2 - 1)) * hum.sample_rate)),
        .release = @max(1, @as(u32, @intFromFloat(burst.release_s * hum.sample_rate))),
        .gain = pan(hum.noise_spread * 0.9, rng),
    };
    for (burst.partials, 0..) |p, k| {
        slot.inc[k] = p.freq * detune * inv_rate;
        // A phase of its own, so two overlapping bursts of the same kind do
        // not reinforce into one twice as loud.
        slot.phase[k] = rng.float(f64);
    }
}

/// Equal-power gains for a source placed `spread` either side of centre.
fn pan(spread: f64, rng: std.Random) [2]f64 {
    const angle = std.math.pi / 4.0 + spread * 0.5 * rng.floatNorm(f64);
    return .{ @cos(angle) * std.math.sqrt2, @sin(angle) * std.math.sqrt2 };
}

/// A geometric draw with the given mean, at least 1 and at most 32.
fn geometric(rng: std.Random, mean: f64) u32 {
    if (!(mean > 1)) return 1;
    const p = 1 / mean;
    const u = @max(rng.float(f64), 1e-12);
    const n = 1 + @floor(@log(u) / @log(1 - p));
    return @intFromFloat(std.math.clamp(n, 1, 32));
}

/// One sounding burst.
const BurstVoice = struct {
    /// `null` when the voice is free.
    burst: ?*const Burst = null,
    phase: [burst_table.max_partials]f64 = @splat(0),
    inc: [burst_table.max_partials]f64 = @splat(0),
    amp: f64 = 0,
    gain: [2]f64 = .{ 1, 1 },
    /// Frames since it started.
    frame: u32 = 0,
    attack: u32 = 0,
    hold: u32 = 0,
    release: u32 = 0,

    fn add(v: *BurstVoice, dst: *[2]f64, channels: u8) void {
        const burst = v.burst orelse return;

        var sample: f64 = 0;
        for (burst.partials, 0..) |p, k| {
            sample += p.amp * sine(v.phase[k]);
            v.phase[k] = wrap(v.phase[k] + v.inc[k]);
        }
        sample *= v.amp * v.envelope();

        if (channels == 1) {
            dst[0] += sample;
        } else {
            dst[0] += sample * v.gain[0];
            dst[1] += sample * v.gain[1];
        }

        v.frame += 1;
        if (v.frame >= v.attack + v.hold + v.release) v.burst = null;
    }

    /// A raised cosine up, a flat hold, a raised cosine down.
    ///
    /// Flat-topped rather than struck-and-decaying because that is what the
    /// recording shows: the ping reaches full level within 5 ms, holds
    /// within 3 dB for 70 ms, and is gone 10 ms after that. Raised cosine at
    /// both ends rather than linear because a corner in the envelope is a
    /// click, and at this level a click is the only thing anybody would
    /// hear.
    fn envelope(v: BurstVoice) f64 {
        if (v.frame < v.attack) {
            const u = @as(f64, @floatFromInt(v.frame)) / @as(f64, @floatFromInt(v.attack));
            return 0.5 - 0.5 * @cos(std.math.pi * u);
        }
        const held = v.frame - v.attack;
        if (held < v.hold) return 1;
        const out = held - v.hold;
        if (out >= v.release) return 0;
        const u = @as(f64, @floatFromInt(out)) / @as(f64, @floatFromInt(v.release));
        return 0.5 + 0.5 * @cos(std.math.pi * u);
    }
};

/// A delay line per channel with damped feedback: the room the bursts are
/// heard in.
const Echo = struct {
    line: [2][]f64,
    pos: [2]usize,
    damping: [2]OnePole,
    feedback: f64,
    level: f64,
    channels: u8,

    fn init(gpa: Allocator, shape: EchoShape, rate: f64, channels: u8) Allocator.Error!Echo {
        var line: [2][]f64 = .{ &.{}, &.{} };
        errdefer for (line) |l| gpa.free(l);
        for (&line, shape.delay_s) |*l, seconds| {
            const frames: usize = @intFromFloat(@max(1, seconds * rate));
            l.* = try gpa.alloc(f64, frames);
            @memset(l.*, 0);
        }
        return .{
            .line = line,
            .pos = .{ 0, 0 },
            .damping = .{ .init(shape.damping_hz, rate), .init(shape.damping_hz, rate) },
            .feedback = std.math.clamp(shape.feedback, 0, 0.95),
            .level = std.math.pow(f64, 10, shape.level_db / 20),
            .channels = channels,
        };
    }

    fn deinit(e: *Echo, gpa: Allocator) void {
        for (e.line) |l| gpa.free(l);
        e.* = undefined;
    }

    /// Add each channel's repeats to it, and feed the dry signal in.
    fn process(e: *Echo, dry: *[2]f64, channels: u8) void {
        for (0..channels) |ch| {
            const tap = e.line[ch][e.pos[ch]];
            e.line[ch][e.pos[ch]] = dry[ch] + e.damping[ch].step(tap) * e.feedback;
            e.pos[ch] = (e.pos[ch] + 1) % e.line[ch].len;
            dry[ch] += tap * e.level;
        }
    }
};

/// One partial.
const Voice = struct {
    freq: f64,
    amp: f64,
    /// Turns, one per channel. Stereo gives the two channels different
    /// phases; mono uses only the first.
    phase: [2]f64,
    gain: [2]f64,
    freq_wander: Wander,
    amp_wander: Wander,
    /// Phase increment per frame, in turns. Set by `retune`.
    inc: f64,
    /// Current amplitude, ramped towards the wander's target across a block.
    level: f64,
    level_step: f64,
};

/// A bounded random walk: white noise through a one-pole lowpass, scaled so
/// that it settles at a known RMS.
///
/// A plain random walk is unbounded and would eventually take a partial
/// anywhere; a sine LFO is periodic, which is the thing being avoided. This
/// is neither: it wanders, and it stays.
const Wander = struct {
    state: f64 = 0,
    pole: f64,
    gain: f64,

    fn init(tau_seconds: f64, rate: f64, rms: f64) Wander {
        const a = @exp(-1 / (tau_seconds * rate));
        // y = a*y + g*x with x of unit variance settles at var(y) =
        // g^2/(1-a^2), so this is the g that makes var(y) come out at rms^2.
        return .{ .pole = a, .gain = rms * @sqrt(1 - a * a) };
    }

    fn next(w: *Wander, rng: std.Random) f64 {
        w.state = w.state * w.pole + rng.floatNorm(f64) * w.gain;
        return w.state;
    }
};

/// The broadband floor: hiss and rumble, each a band of shaped white noise.
///
/// Two lowpass poles on the hiss rather than one, because the recording's
/// floor falls at about 12 dB an octave above a few hundred hertz and one
/// pole falls at 6.
const NoiseFloor = struct {
    rng: std.Random.DefaultPrng,
    hiss: Band,
    rumble: Band,

    fn init(seed: u64, options: Options, rate: f64) NoiseFloor {
        return .{
            .rng = .init(seed),
            .hiss = .init(options.hiss, rate),
            .rumble = .init(options.rumble, rate),
        };
    }

    fn next(f: *NoiseFloor, hiss_gain: f64, rumble_gain: f64) f64 {
        // Its own generator rather than the caller's: the noise floor is then
        // independent of how many partials there are, so adding one does not
        // change the hiss.
        const rng = f.rng.random();
        // Skipped rather than multiplied by zero when the hiss is off, which
        // is the default: six filter stages a sample, per channel, for a
        // signal nobody asked for.
        const hiss = if (hiss_gain != 0) f.hiss.next(rng) * hiss_gain else 0;
        return hiss + f.rumble.next(rng) * rumble_gain;
    }
};

/// White noise shaped by a `BandShape` and scaled back to unit RMS.
///
/// The scaling is what lets the corners be moved without also having to
/// re-tune the level: a `_level_db` option means the same loudness whatever
/// band it is spread across.
const Band = struct {
    high: [2]OnePole,
    tilt: [2]OnePole,
    top: [2]Biquad,
    scale: f64,

    fn init(shape: BandShape, rate: f64) Band {
        const high: OnePole = .init(shape.high_hz, rate);
        const tilt: OnePole = .init(shape.tilt_hz, rate);
        const top: Biquad = .lowpass(shape.cutoff_hz, rate);
        return .{
            .high = @splat(high),
            .tilt = @splat(tilt),
            .top = @splat(top),
            .scale = 1 / responseRms(high, tilt, top),
        };
    }

    fn next(b: *Band, rng: std.Random) f64 {
        var v = rng.floatNorm(f64) * b.scale;
        // A one-pole highpass is what the matching lowpass does not pass.
        for (&b.high) |*p| v -= p.step(v);
        for (&b.tilt) |*p| v = p.step(v);
        for (&b.top) |*q| v = q.step(v);
        return v;
    }

    /// The RMS this chain produces from unit-variance white noise.
    ///
    /// Integrated over the spectrum rather than composed from each stage's
    /// own figure, because no stage after the first is being fed white noise
    /// and so none of them attenuates it by the factor it would.
    fn responseRms(high: OnePole, tilt: OnePole, top: Biquad) f64 {
        const steps = 4096;
        var sum: f64 = 0;
        for (0..steps) |i| {
            const w = std.math.pi * (@as(f64, @floatFromInt(i)) + 0.5) / steps;
            const highed = high.highpassPower(w);
            const tilted = tilt.lowpassPower(w);
            const topped = top.power(w);
            sum += highed * highed * tilted * tilted * topped * topped;
        }
        return @sqrt(sum / steps);
    }
};

const OnePole = struct {
    state: f64 = 0,
    pole: f64,

    fn init(cutoff_hz: f64, rate: f64) OnePole {
        return .{ .pole = @exp(-2 * std.math.pi * cutoff_hz / rate) };
    }

    fn step(p: *OnePole, in: f64) f64 {
        p.state = in * (1 - p.pole) + p.state * p.pole;
        return p.state;
    }

    /// |H(w)|^2 for `step` itself: y = (1-a)x + a*y[-1].
    fn lowpassPower(p: OnePole, w: f64) f64 {
        const a = p.pole;
        return (1 - a) * (1 - a) / (1 - 2 * a * @cos(w) + a * a);
    }

    /// |H(w)|^2 for `in - step(in)`, which is the matching highpass:
    /// H = a(1 - e^-jw) / (1 - a e^-jw).
    fn highpassPower(p: OnePole, w: f64) f64 {
        const a = p.pole;
        return a * a * (2 - 2 * @cos(w)) / (1 - 2 * a * @cos(w) + a * a);
    }
};

/// A second-order section, in the direct form and with the Audio EQ Cookbook
/// coefficients.
///
/// Only the lowpass is here, because only the lowpass is needed. What it has
/// that a pair of one-poles does not is a double zero at Nyquist, so its
/// stopband keeps going down instead of settling on a floor.
const Biquad = struct {
    b0: f64,
    b1: f64,
    b2: f64,
    a1: f64,
    a2: f64,
    x1: f64 = 0,
    x2: f64 = 0,
    y1: f64 = 0,
    y2: f64 = 0,

    fn lowpass(cutoff_hz: f64, rate: f64) Biquad {
        // Clamped below Nyquist: the bilinear transform puts the corner at
        // infinity there, and a `--rate` low enough to reach it should give
        // a dull noise floor rather than a division by zero.
        const cutoff = @min(cutoff_hz, rate * 0.45);
        const w0 = 2 * std.math.pi * cutoff / rate;
        // Butterworth: the flattest passband a single section has.
        const q = std.math.sqrt1_2;
        const alpha = @sin(w0) / (2 * q);
        const cos_w0 = @cos(w0);
        const a0 = 1 + alpha;
        return .{
            .b0 = (1 - cos_w0) / 2 / a0,
            .b1 = (1 - cos_w0) / a0,
            .b2 = (1 - cos_w0) / 2 / a0,
            .a1 = -2 * cos_w0 / a0,
            .a2 = (1 - alpha) / a0,
        };
    }

    fn step(q: *Biquad, x: f64) f64 {
        const y = q.b0 * x + q.b1 * q.x1 + q.b2 * q.x2 - q.a1 * q.y1 - q.a2 * q.y2;
        q.x2 = q.x1;
        q.x1 = x;
        q.y2 = q.y1;
        q.y1 = y;
        return y;
    }

    fn power(q: Biquad, w: f64) f64 {
        const c1 = @cos(w);
        const s1 = @sin(w);
        const c2 = @cos(2 * w);
        const s2 = @sin(2 * w);
        const num_re = q.b0 + q.b1 * c1 + q.b2 * c2;
        const num_im = -(q.b1 * s1 + q.b2 * s2);
        const den_re = 1 + q.a1 * c1 + q.a2 * c2;
        const den_im = -(q.a1 * s1 + q.a2 * s2);
        return (num_re * num_re + num_im * num_im) / (den_re * den_re + den_im * den_im);
    }
};

inline fn sine(turns: f64) f64 {
    return @sin(turns * 2 * std.math.pi);
}

inline fn wrap(turns: f64) f64 {
    return turns - @floor(turns);
}