Something went wrong. Try again.
A synthesiser for the low machine-room hum of the WOPR computer in WarGames, built from partials measured off a recording of it.
Something went wrong. Try again.
1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162636465666768697071727374757677787980818283848586878889909192939495969798991001011021031041051061071081091101111121131141151161171181191201211221231241251261271281291301311321331341351361371381391401411421431441451461471481491501511521531541551561571581591601611621631641651661671681691701711721731741751761771781791801811821831841851861871881891901911921931941951961971981992002012022032042052062072082092102112122132142152162172182192202212222232242252262272282292302312322332342352362372382392402412422432442452462472482492502512522532542552562572582592602612622632642652662672682692702712722732742752762772782792802812822832842852862872882892902912922932942952962972982993003013023033043053063073083093103113123133143153163173183193203213223233243253263273283293303313323333343353363373383393403413423433443453463473483493503513523533543553563573583593603613623633643653663673683693703713723733743753763773783793803813823833843853863873883893903913923933943953963973983994004014024034044054064074084094104114124134144154164174184194204214224234244254264274284294304314324334344354364374384394404414424434444454464474484494504514524534544554564574584594604614624634644654664674684694704714724734744754764774784794804814824834844854864874884894904914924934944954964974984995005015025035045055065075085095105115125135145155165175185195205215225235245255265275285295305315325335345355365375385395405415425435445455465475485495505515525535545555565575585595605615625635645655665675685695705715725735745755765775785795805815825835845855865875885895905915925935945955965975985996006016026036046056066076086096106116126136146156166176186196206216226236246256266276286296306316326336346356366376386396406416426436446456466476486496506516526536546556566576586596606616626636646656666676686696706716726736746756766776786796806816826836846856866876886896906916926936946956966976986997007017027037047057067077087097107117127137147157167177187197207217227237247257267277287297307317327337347357367377387397407417427437447457467477487497507517527537547557567577587597607617627637647657667677687697707717727737747757767777787797807817827837847857867877887897907917927937947957967977987998008018028038048058068078088098108118128138148158168178188198208218228238248258268278288298308318328338348358368378388398408418428438448458468478488498508518528538548558568578588598608618628638648658668678688698708718728738748758768778788798808818828838848858868878888898908918928938948958968978988999009019029039049059069079089099109119129139149159169179189199209219229239249259269279289299309319329339349359369379389399409419429439449459469479489499509519529539549559569579589599609619629639649659669679689699709719729739749759769779789799809819829839849859869879889899909919929939949959969979989991000100110021003100410051006100710081009101010111012// 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);}