Anti-aliased pickup: solve at R fs and low-pass filter the pickup down to fs (both solvers); CACHE_VERSION 3; chaos-check analysis scripts
10 changed files+478−10
README.mdmodified+3−2View file
@@ -17,7 +17,7 @@ npm run build # static site in dist/
1717 - Ends: clamped (fixed point), sliding (massless ring on a rod: x fixed, free across, holds tension), or free (no support; only consistent with `T = 0`). Any bead can be pinned.
1818 - Excitation: pluck (triangle between the nearest supports, released at rest), strike (Gaussian initial velocity), hammer (raised-cosine force pulse). Direction is an angle in the y-z plane.
1919 - Damping from two decay times: `sigma(omega) = sigma0 + sigma1 omega^2`, fitted to T60 at two frequencies for transverse modes. `sigma0` is mass-proportional; `sigma1` is a dashpot on each spring's relative velocity.
20-- Audio: pickup acceleration (forces / mass, at the full-step velocity), summed over all beads (the default; on a uniform string clamped at both ends this cancels the even modes) or from one bead, along y, z or x. The `.wav` is normalized to -1 dBFS.
20+- Audio: pickup acceleration (forces / mass, at the full-step velocity), summed over all beads (the default; on a uniform string clamped at both ends this cancels the even modes) or from one bead, along y, z or x. The `.wav` is normalized to -1 dBFS. The pickup is anti-aliased (see Solvers).
2121 - Solver: velocity Verlet (explicit, symplectic without damping). `dt = cfl * 2 / omega_max`, rounded down to a whole number of steps per audio sample; `omega_max` is the highest frequency of the system linearized at rest (power iteration, axial `k` and transverse `T/l` per spring, bending, plus a mallet's contact stiffness). See Solvers.
2222 - 3D model view (three.js): Setup shows the beads at rest with supports (wall = clamped, translucent plane + ring = sliding), pinned beads (cubes), excitation (pluck shape, strike velocity arrows, hammer force arrows) and pickup (bead + axis arrow). Replay animates the saved frames ("Replay frames" per second, default 2000) with adjustable exaggeration, on a timeline shared with the 2D side/top views and space-time image.
2323 - Checks shown per run: energy balance (energy + dissipated - work in - initial) relative to peak energy, largest strain. Linear transverse and longitudinal modes (Sturm bisection on the tridiagonal stiffness) are marked on the spectrum.
@@ -49,6 +49,7 @@ String tests (`test/string.test.ts`): modes vs the discrete-string formula; smal
4949
5050 Both solvers integrate the same scheme (`src/sim/solver.ts` on the CPU in double precision, `src/sim/gpu.ts` with WebGPU in single precision); the WebGPU toggle next to Run picks one, and a failed GPU run falls back to the CPU with the reason in the run summary. The toggle is on by default for the drum and cymbal and off for the string, whose few dozen beads run faster on the CPU.
5151
52+- Anti-aliased pickup (`src/sim/antialias.ts`, applied in the worker to both solvers): the solvers record the pickup once per sample, without a filter, so motion above fs/2 folded back into the audible band (a 10 mm pluck of a 5-bead string put 28% of its pickup power above 24 kHz). A run is now solved at R fs, with R = ceil((2 f_max + fs/2) / fs) capped at the steps per sample (f_max = omega_max / 2 pi, the factor 2 leaving room for nonlinear combination tones), and the pickup is low-pass filtered (Kaiser-windowed sinc, flat to 0.45 fs, at least 86 dB down from 0.5 fs) and kept every R-th sample. The contact force is subsampled. Presets get R = 1 (low f_max) to 6; the steps per sample round up to a multiple of R (11 to 12 for the steel string). The cost is the per-sample work done R times as often: about 25% for the strings, and 68% on the CPU and 36% with WebGPU for the default cymbal (R = 2), where recomputing the velocity-dependent forces at every internal sample dominates. Recording the extra pickup samples inside the solvers without that recomputation would recover most of it. CPU and WebGPU audio still agree to 0.11% RMS over the first 50 ms of the default cymbal. Tests: `test/antialias.test.ts`.
5253 - Time step: `omega_max` comes from power iteration on the linearized stiffness rather than a Gershgorin bound, which overestimated it by 1.4 to 1.5x on the cymbal (so steps were that much shorter than needed).
5354 - CPU: at each audio sample the solver recomputes only the velocity-dependent forces (damping, friction) at the full-step velocity; spring and hinge elasticity depend on positions, which have not changed. Together with the time step, this made the cymbal presets 1.4 to 1.6x faster.
5455 - WebGPU: each phase of a time step (kick and drift, hinge angles and contact, hinge moments, bead forces, kick with rigid-motion removal) is a small kernel; a command buffer holds the dispatches for many samples, and step and sample counters live on the GPU. Forces are gathered per bead from fixed-width adjacency lists, so there are no atomics. To keep small vibrations from vanishing in f32, the state is the displacement from rest, and spring extensions and hinge angles are computed from precomputed rest differences. Sums over beads (rigid motion, energy, audio) go through per-workgroup partials. On the first 50 ms of every preset, the audio agrees with the CPU to within 0.5% RMS (0.1% for most), and the energy-balance check and largest strain match.
@@ -61,7 +62,7 @@ Both solvers integrate the same scheme (`src/sim/solver.ts` on the CPU in double
6162 Runs are cached on [artifact-cache](https://artifact-cache.tempory.net) (app `toy-instruments`, `src/sim/cache.ts`). Reads are public, so anyone opening the site gets a result someone has already computed; writes need the app's passcode, entered once per browser under "Shared cache" in the panel.
6263
6364 - Key: `v<CACHE_VERSION>/<instrument>/<cpu|gpu>/<SHA-256>`, hashed over everything the solver receives (bead system, excitation, options), not the UI parameters. A model change that alters the built system therefore misses old entries on its own; bump `CACHE_VERSION` when the solver's output changes for the same inputs. CPU and WebGPU results are cached separately; a WebGPU run that fell back to the CPU is stored as a CPU result.
64-- Value: the whole `SimResult` in a small binary format, so a cached result shows everything a fresh run does, split into two entries: `<key>/light` (audio, energy, contact force, scalars; a few hundred kB) and `<key>/frames` (the replay frames; the default cymbal's are about 35 MB). The frames are stored first, so a light entry always has its frames. Results whose frames are over the 100 MB entry limit are not stored. Version 1 kept both in one entry.
65+- Value: the whole `SimResult` in a small binary format, so a cached result shows everything a fresh run does, split into two entries: `<key>/light` (audio, energy, contact force, scalars; a few hundred kB) and `<key>/frames` (the replay frames; the default cymbal's are about 35 MB). The frames are stored first, so a light entry always has its frames. Results whose frames are over the 100 MB entry limit are not stored. Version 1 kept both in one entry; version 3 has the anti-aliased pickup.
6566 - As the settings (or the WebGPU toggle) change, the page checks whether the cache has a result for them (a HEAD request, debounced) and, if so, shows its date and size with a Load cached button next to Run. Run always simulates; when the passcode is set, every completed run is stored, replacing any existing entry.
6667 - Because the key hashes floating-point inputs, a browser whose `Math` functions round differently in the last bit (possible in Safari) builds a different key and misses entries computed elsewhere. It never gets a wrong result.
6768
analysis/chaos-check.tsadded+209−0View file
@@ -0,0 +1,209 @@
1+// Is the large-pluck motion of the 5-bead string chaotic, or a solver artifact?
2+//
3+// Runs the post's large pluck (n=5, tension 5.001 N, bead 2, 10 mm, 2 s) with the app's own
4+// model and CPU solver (double precision), and measures
5+// 1. step size: runs at cfl 0.1, 0.05, ..., halving; each against the next finer one.
6+// For chaos with a convergent solver, the runs agree up to a time that grows by about
7+// ln(4)/lambda per halving (velocity Verlet is second order, so the error drops 4x).
8+// 2. perturbations: pluck amplitude times (1 + eps), at two step sizes. The separation
9+// should grow exponentially at the same rate lambda (the largest Lyapunov exponent)
10+// at both step sizes, and its divergence time should grow like ln(1/eps)/lambda.
11+// 3. energy balance error of every run (velocity Verlet nearly conserves energy).
12+// 4. statistics: band levels of the sound (octave bands) at every step size.
13+// 5. the gentle pluck (0.1 mm) against the linear modes.
14+// 6. aliasing: the solver point-samples the pickup at fs (48 kHz) without a low-pass
15+// filter, so sound above fs/2 folds back into the spectrogram. Runs at 16 fs and
16+// reports the fraction of the pickup's power above 24 kHz.
17+//
18+// Positions are saved at every audio sample (frameRate = fs), since the runs separate
19+// within milliseconds; the replay rate does not affect the dynamics.
20+//
21+// npx tsx analysis/chaos-check.ts [out.json] (about a minute per fine run)
22+
23+import { writeFileSync } from 'node:fs';
24+import { CONFIGS } from '../src/instruments/configs';
25+import { parseSettings, simOptions } from '../src/ui/settings';
26+import { simulate } from '../src/sim/solver';
27+import { spectrum } from '../src/ui/fft';
28+import type { SimResult } from '../src/sim/types';
29+
30+const BASE = 'n=5&tension=5.001&excBead=2';
31+const AMP = 0.01;
32+
33+interface Run {
34+ label: string;
35+ cfl: number;
36+ amp: number;
37+ r: SimResult;
38+}
39+
40+function run(label: string, amp: number, cfl: number, over: { fs?: number; frameRate?: number } = { frameRate: 48000 }): Run {
41+ const cfg = CONFIGS['bead-string'];
42+ const s = parseSettings(cfg, `${BASE}&pluckAmp=${amp}&cfl=${cfl}`);
43+ // set directly: a link only accepts the app's sample rates
44+ Object.assign(s.params, over);
45+ const built = cfg.setup(s.params).built!;
46+ const r = simulate(built.sys, built.exc, simOptions(s.params, built));
47+ const emax = Math.max(...r.energy);
48+ let bal = 0;
49+ for (const b of r.balance) bal = Math.max(bal, Math.abs(b));
50+ console.error(`${label}: dt ${r.dt.toExponential(3)} s (${r.stepsPerSample}/sample), ${r.wallSeconds.toFixed(1)} s wall, energy error ${((bal / emax) * 100).toExponential(2)}% of peak`);
51+ return { label, cfl, amp, r };
52+}
53+
54+// RMS over beads of the displacement difference at each frame (m)
55+function separation(a: SimResult, b: SimResult): Float64Array {
56+ const nf = a.frameTimes.length;
57+ const w = a.frames.length / nf;
58+ const d = new Float64Array(nf);
59+ for (let f = 0; f < nf; f++) {
60+ let s = 0;
61+ for (let i = 0; i < w; i++) {
62+ const x = a.frames[f * w + i] - b.frames[f * w + i];
63+ s += x * x;
64+ }
65+ d[f] = Math.sqrt(s / (w / 3));
66+ }
67+ return d;
68+}
69+
70+// first time the separation exceeds thr (s), or null
71+function crossing(t: Float64Array, d: Float64Array, thr: number): number | null {
72+ for (let i = 0; i < d.length; i++) if (d[i] > thr) return t[i];
73+ return null;
74+}
75+
76+// growth rate of the running maximum of ln d between lo and hi (1/s): least squares
77+function growthRate(t: Float64Array, d: Float64Array, lo: number, hi: number): { lambda: number; t0: number; t1: number } | null {
78+ const xs: number[] = [];
79+ const ys: number[] = [];
80+ let m = 0;
81+ for (let i = 0; i < d.length; i++) {
82+ m = Math.max(m, d[i]);
83+ if (m >= lo && m <= hi) {
84+ xs.push(t[i]);
85+ ys.push(Math.log(m));
86+ }
87+ }
88+ if (xs.length < 10) return null;
89+ const n = xs.length;
90+ const mx = xs.reduce((s, x) => s + x, 0) / n;
91+ const my = ys.reduce((s, y) => s + y, 0) / n;
92+ let sxy = 0;
93+ let sxx = 0;
94+ for (let i = 0; i < n; i++) {
95+ sxy += (xs[i] - mx) * (ys[i] - my);
96+ sxx += (xs[i] - mx) ** 2;
97+ }
98+ return { lambda: sxy / sxx, t0: xs[0], t1: xs[n - 1] };
99+}
100+
101+// sound level in octave bands (dB re the loudest band), from 1 s on (the settled motion)
102+function bands(r: SimResult): { f: number; db: number }[] {
103+ const start = Math.round(1 * r.fs);
104+ const s = spectrum(r.audio.subarray(start), r.fs);
105+ const out: { f: number; db: number }[] = [];
106+ for (let fc = 125; fc <= 16000; fc *= 2) {
107+ let p = 0;
108+ for (let k = 0; k < s.freq.length; k++) if (s.freq[k] >= fc / Math.SQRT2 && s.freq[k] < fc * Math.SQRT2) p += 10 ** (s.db[k] / 10);
109+ out.push({ f: fc, db: 10 * Math.log10(p) });
110+ }
111+ const top = Math.max(...out.map((b) => b.db));
112+ return out.map((b) => ({ f: b.f, db: b.db - top }));
113+}
114+
115+const downsample = (x: Float64Array, k: number) => Array.from(x.filter((_, i) => i % k === 0));
116+
117+const out: Record<string, unknown> = { settings: `${BASE}&pluckAmp=${AMP}` };
118+
119+// 1. step size
120+const CFLS = [0.1, 0.05, 0.025, 0.0125, 0.00625, 0.003125];
121+const byCfl = CFLS.map((c) => run(`cfl ${c}`, AMP, c));
122+const t = byCfl[0].r.frameTimes;
123+const DS = 24; // keep every 24th frame (0.5 ms) in the output
124+out.time = downsample(t, DS);
125+out.stepSize = byCfl.slice(0, -1).map((a, i) => {
126+ const b = byCfl[i + 1];
127+ const d = separation(a.r, b.r);
128+ return {
129+ cfl: a.cfl,
130+ against: b.cfl,
131+ dt: a.r.dt,
132+ divergeAt1um: crossing(t, d, 1e-6),
133+ divergeAt1mm: crossing(t, d, 1e-3),
134+ d: downsample(d, DS),
135+ };
136+});
137+
138+// 2. perturbations at two step sizes
139+const EPS = [1e-5, 1e-7, 1e-9, 1e-11];
140+const pert: Record<string, unknown>[] = [];
141+for (const [ci, cfl] of [[0, 0.1], [2, 0.025]] as const) {
142+ const base = byCfl[ci];
143+ for (const eps of EPS) {
144+ const p = run(`cfl ${cfl}, eps ${eps}`, AMP * (1 + eps), cfl);
145+ const d = separation(base.r, p.r);
146+ pert.push({
147+ cfl,
148+ eps,
149+ divergeAt1mm: crossing(t, d, 1e-3),
150+ growth: growthRate(t, d, 1e-8, 1e-4),
151+ bands: bands(p.r),
152+ d: downsample(d, DS),
153+ });
154+ }
155+}
156+out.perturbation = pert;
157+
158+// 3. energy balance and 4. sound statistics per step size
159+out.runs = byCfl.map(({ cfl, r }) => {
160+ const emax = Math.max(...r.energy);
161+ let bal = 0;
162+ for (const b of r.balance) bal = Math.max(bal, Math.abs(b));
163+ return { cfl, dt: r.dt, stepsPerSample: r.stepsPerSample, wallSeconds: r.wallSeconds, energyError: bal / emax, maxStretch: r.maxStretch, bands: bands(r) };
164+});
165+
166+// 5. gentle pluck vs the linear modes: f_k = (1/pi) sqrt(T / (m a)) sin(k pi / 12)
167+const gentle = run('gentle 0.1 mm', 0.0001, 0.1);
168+{
169+ const s = spectrum(gentle.r.audio, gentle.r.fs);
170+ const T = 5.001;
171+ const m = 5.08e-5;
172+ const a = 0.0127;
173+ const linear = [1, 2, 3, 4, 5].map((k) => ({ k, hz: (1 / Math.PI) * Math.sqrt(T / (m * a)) * Math.sin((k * Math.PI) / 12) }));
174+ const top = Math.max(...s.db);
175+ out.gentle = linear.map(({ k, hz }) => {
176+ // strongest bin within 3% of the predicted frequency
177+ let best = -1;
178+ for (let i = 0; i < s.freq.length; i++) if (Math.abs(s.freq[i] - hz) < 0.03 * hz && (best < 0 || s.db[i] > s.db[best])) best = i;
179+ return { k, predictedHz: hz, peakHz: s.freq[best], levelDb: s.db[best] - top };
180+ });
181+}
182+
183+// 6. aliasing
184+out.aliasing = [0.001, AMP].map((amp) => {
185+ const { r } = run(`${amp * 1000} mm at 16 fs`, amp, 0.1, { fs: 768000, frameRate: 2000 });
186+ const s = spectrum(r.audio, r.fs);
187+ let total = 0;
188+ let above = 0;
189+ let above20 = 0;
190+ for (let k = 0; k < s.freq.length; k++) {
191+ const p = 10 ** (s.db[k] / 10);
192+ total += p;
193+ if (s.freq[k] > 24000) above += p;
194+ if (s.freq[k] > 20000) above20 += p;
195+ }
196+ // octave bands up to 384 kHz, dB re the loudest
197+ const b: { f: number; db: number }[] = [];
198+ for (let fc = 125; fc <= 256000; fc *= 2) {
199+ let p = 0;
200+ for (let k = 0; k < s.freq.length; k++) if (s.freq[k] >= fc / Math.SQRT2 && s.freq[k] < fc * Math.SQRT2) p += 10 ** (s.db[k] / 10);
201+ b.push({ f: fc, db: 10 * Math.log10(p) });
202+ }
203+ const top = Math.max(...b.map((x) => x.db));
204+ return { amp, fractionAbove24k: above / total, fractionAbove20k: above20 / total, bands: b.map((x) => ({ f: x.f, db: x.db - top })) };
205+});
206+
207+const file = process.argv[2] ?? 'chaos-check.json';
208+writeFileSync(file, JSON.stringify(out));
209+console.error(`wrote ${file}`);
analysis/plot_chaos.pyadded+98−0View file
@@ -0,0 +1,98 @@
1+"""Plots for analysis/chaos-check.ts: python3 analysis/plot_chaos.py chaos-check.json out_dir"""
2+
3+import json
4+import sys
5+from pathlib import Path
6+
7+import matplotlib
8+
9+matplotlib.use("Agg")
10+import matplotlib.pyplot as plt
11+import numpy as np
12+
13+d = json.load(open(sys.argv[1]))
14+out = Path(sys.argv[2] if len(sys.argv) > 2 else ".")
15+out.mkdir(parents=True, exist_ok=True)
16+t = np.array(d["time"]) * 1000 # ms
17+TMAX = 30
18+
19+
20+def sep_axes(ax, title):
21+ ax.set_yscale("log")
22+ ax.set_xlim(0, TMAX)
23+ ax.set_ylim(1e-12, 3e-2)
24+ ax.set_xlabel("time (ms)")
25+ ax.set_ylabel("RMS separation of the beads (m)")
26+ ax.set_title(title, fontsize=10)
27+ ax.grid(alpha=0.3)
28+
29+
30+# 1. step size
31+fig, ax = plt.subplots(figsize=(6.5, 4))
32+for s in d["stepSize"]:
33+ y = np.maximum(np.array(s["d"]), 1e-13)
34+ ax.plot(t, y, lw=1, label=f"cfl {s['cfl']} vs {s['against']}")
35+sep_axes(ax, "Same pluck (10 mm), time step halved: runs agree for longer")
36+ax.legend(fontsize=8)
37+fig.tight_layout()
38+fig.savefig(out / "step-size.png", dpi=150)
39+
40+# 2. perturbations
41+fig, axs = plt.subplots(1, 2, figsize=(11, 4), sharey=True)
42+for ax, cfl in zip(axs, sorted({p["cfl"] for p in d["perturbation"]}, reverse=True)):
43+ for p in [p for p in d["perturbation"] if p["cfl"] == cfl]:
44+ y = np.maximum(np.array(p["d"]), 1e-13)
45+ g = p["growth"]
46+ lab = f"eps {p['eps']:g}" + (f", growth {g['lambda']:.0f}/s" if g else "")
47+ ax.plot(t, y, lw=1, label=lab)
48+ sep_axes(ax, f"Pluck amplitude times (1 + eps), cfl {cfl}")
49+ ax.legend(fontsize=8)
50+fig.tight_layout()
51+fig.savefig(out / "perturbation.png", dpi=150)
52+
53+# 3. energy error vs time step
54+fig, ax = plt.subplots(figsize=(5, 4))
55+dt = np.array([r["dt"] for r in d["runs"]])
56+err = np.array([r["energyError"] for r in d["runs"]])
57+ax.loglog(dt, err, "o-", label="largest energy balance error / peak energy")
58+ax.loglog(dt, err[0] * (dt / dt[0]) ** 2, "k--", lw=0.8, label="slope 2")
59+ax.set_xlabel("time step (s)")
60+ax.set_ylabel("relative energy error")
61+ax.legend(fontsize=8)
62+ax.grid(alpha=0.3, which="both")
63+fig.tight_layout()
64+fig.savefig(out / "energy.png", dpi=150)
65+
66+# 4. sound statistics: octave bands, step sizes vs perturbations
67+fig, ax = plt.subplots(figsize=(6.5, 4))
68+for r in d["runs"]:
69+ ax.plot([b["f"] for b in r["bands"]], [b["db"] for b in r["bands"]], "o-", color="C0", alpha=0.6, lw=1)
70+for p in d["perturbation"]:
71+ ax.plot([b["f"] for b in p["bands"]], [b["db"] for b in p["bands"]], "s--", color="C1", alpha=0.6, lw=1)
72+ax.plot([], [], "o-", color="C0", label="6 time steps")
73+ax.plot([], [], "s--", color="C1", label="8 tiny perturbations")
74+ax.set_xscale("log")
75+ax.set_xlabel("octave band centre (Hz)")
76+ax.set_ylabel("level, 1 s to 2 s (dB re loudest band)")
77+ax.set_title("Sound statistics: time step vs perturbation", fontsize=10)
78+ax.legend(fontsize=8)
79+ax.grid(alpha=0.3)
80+fig.tight_layout()
81+fig.savefig(out / "bands.png", dpi=150)
82+
83+# 5. aliasing
84+fig, ax = plt.subplots(figsize=(6.5, 4))
85+for a in d["aliasing"]:
86+ ax.plot([b["f"] for b in a["bands"]], [b["db"] for b in a["bands"]], "o-",
87+ label=f"{a['amp'] * 1000:g} mm: {a['fractionAbove24k'] * 100:.1f}% of power above 24 kHz")
88+ax.axvline(24000, color="k", lw=0.8, ls="--")
89+ax.text(25000, -5, "24 kHz", fontsize=8)
90+ax.set_xscale("log")
91+ax.set_xlabel("octave band centre (Hz)")
92+ax.set_ylabel("pickup level (dB re loudest band)")
93+ax.set_title("Pickup at 768 kHz: content above the 48 kHz Nyquist limit", fontsize=10)
94+ax.legend(fontsize=8)
95+ax.grid(alpha=0.3)
96+fig.tight_layout()
97+fig.savefig(out / "aliasing.png", dpi=150)
98+print("wrote", out)
src/sim/antialias.tsadded+86−0View file
@@ -0,0 +1,86 @@
1+import { timeStep } from './solver';
2+import type { BeadSystem, Excitation, SimOptions, SimResult } from './types';
3+
4+// Anti-aliased pickup. The solvers record the pickup once per output sample (point
5+// sampling), so motion above fs/2 would fold back into the audible band: a 10 mm pluck of
6+// the 5-bead string puts about 28% of its pickup power above 24 kHz. Instead, a run is
7+// solved with an internal sample rate R fs, and the pickup is low-pass filtered (a Kaiser-
8+// windowed sinc, passband to 0.45 fs, stopband from 0.5 fs, about 86 dB down) and kept
9+// every R-th sample. R is chosen so that the internal Nyquist frequency is above twice the
10+// system's highest linear frequency plus fs/2 (room for nonlinear combination tones),
11+// capped at the solver's steps per sample, so the time step does not shrink.
12+
13+export function oversampling(sys: BeadSystem, exc: Excitation, opt: SimOptions): number {
14+ const { omegaMax, stepsPerSample } = timeStep(sys, exc.mallet, opt.fs, opt.cfl);
15+ const fmax = omegaMax / (2 * Math.PI);
16+ return Math.max(1, Math.min(stepsPerSample, Math.ceil((2 * fmax + opt.fs / 2) / opt.fs)));
17+}
18+
19+// the options the solver runs with: R times the sample rate, delays in internal samples
20+export function oversampledOptions(opt: SimOptions, R: number): SimOptions {
21+ return { ...opt, fs: opt.fs * R, pickupDelay: opt.pickupDelay?.map((d) => d * R) };
22+}
23+
24+const BETA = 8.6; // Kaiser window: about 86 dB stopband
25+
26+// filter taps at the internal rate, symmetric about the middle, summing to 1
27+const kernels = new Map<number, Float64Array>();
28+export function lowpass(R: number): Float64Array {
29+ let h = kernels.get(R);
30+ if (h) return h;
31+ const df = 0.05 / R; // transition width, cycles per internal sample (0.45 to 0.5 of fs)
32+ const fc = 0.475 / R; // cutoff, cycles per internal sample
33+ let N = Math.ceil((86 - 7.95) / (14.36 * df)) + 1;
34+ if (N % 2 === 0) N++;
35+ const M = (N - 1) / 2;
36+ h = new Float64Array(N);
37+ const i0b = besselI0(BETA);
38+ let sum = 0;
39+ for (let i = 0; i < N; i++) {
40+ const x = i - M;
41+ const sinc = x === 0 ? 2 * fc : Math.sin(2 * Math.PI * fc * x) / (Math.PI * x);
42+ const r = x / M;
43+ h[i] = sinc * (besselI0(BETA * Math.sqrt(Math.max(0, 1 - r * r))) / i0b);
44+ sum += h[i];
45+ }
46+ for (let i = 0; i < N; i++) h[i] /= sum;
47+ kernels.set(R, h);
48+ return h;
49+}
50+
51+function besselI0(x: number): number {
52+ let s = 1;
53+ let term = 1;
54+ for (let k = 1; k < 50; k++) {
55+ term *= (x / (2 * k)) ** 2;
56+ s += term;
57+ if (term < 1e-12 * s) break;
58+ }
59+ return s;
60+}
61+
62+// Filter and decimate an internal-rate result to fs. complete: the whole run (the signal
63+// is zero after its end); otherwise a partial result, of which only the output samples
64+// whose filter window is covered are kept.
65+export function decimateResult(r: SimResult, R: number, complete: boolean): SimResult {
66+ if (R === 1) return { ...r, oversampling: 1 };
67+ const h = lowpass(R);
68+ const M = (h.length - 1) / 2;
69+ const x = r.audio;
70+ const nOut = complete ? Math.ceil(x.length / R) : Math.max(0, Math.floor((x.length - 1 - M) / R) + 1);
71+ const audio = new Float32Array(nOut);
72+ for (let i = 0; i < nOut; i++) {
73+ const c = i * R; // centre, in internal samples; the signal is zero before t = 0
74+ const a = Math.max(0, c - M);
75+ const b = Math.min(x.length - 1, c + M);
76+ let s = 0;
77+ for (let j = a; j <= b; j++) s += h[j - c + M] * x[j];
78+ audio[i] = s;
79+ }
80+ let contactForce: Float32Array | undefined;
81+ if (r.contactForce) {
82+ contactForce = new Float32Array(nOut);
83+ for (let i = 0; i < nOut; i++) contactForce[i] = r.contactForce[Math.min(i * R, r.contactForce.length - 1)];
84+ }
85+ return { ...r, audio, contactForce, fs: r.fs / R, stepsPerSample: r.stepsPerSample * R, oversampling: R };
86+}
src/sim/cache.tsmodified+2−2View file
@@ -7,9 +7,9 @@ import type { BeadSystem, Excitation, SimOptions, SimResult } from './types';
77 // scalars: everything but the replay frames, small) and <key>/frames (the replay frames,
88 // usually most of the size). The hash covers everything the solver sees, so model
99 // changes miss old entries by themselves; bump CACHE_VERSION when the solver's output
10-// changes for the same inputs.
10+// changes for the same inputs. Version 3: the pickup is anti-aliased (antialias.ts).
1111
12-export const CACHE_VERSION = 2;
12+export const CACHE_VERSION = 3;
1313 const BASE = 'https://artifact-cache.tempory.net/api/apps/toy-instruments';
1414 const PASS_KEY = 'toy-instruments.cache-passcode';
1515 export const MAX_ENTRY_BYTES = 100 * 1024 * 1024;
src/sim/types.tsmodified+3−0View file
@@ -98,4 +98,7 @@ export interface SimResult {
9898 malletPos?: Float32Array; // mallet position q (3 per frame) at frame times (m)
9999 backend?: string; // which solver ran, e.g. 'CPU' or 'WebGPU (intel gen-12lp)'
100100 fallback?: string; // why WebGPU was asked for but not used
101+ // the pickup was computed at this many times fs and low-pass filtered (antialias.ts);
102+ // stepsPerSample counts steps per output sample
103+ oversampling?: number;
101104 }
src/sim/worker.tsmodified+11−2View file
@@ -1,5 +1,6 @@
11 import { simulate } from './solver';
22 import { simulateGPU } from './gpu';
3+import { decimateResult, oversampledOptions, oversampling } from './antialias';
34 import type { BeadSystem, Excitation, SimOptions, SimResult } from './types';
45
56 export interface WorkerRequest {
@@ -13,9 +14,16 @@ const buffers = (r: SimResult) =>
1314 [r.audio, r.frames, r.frameTimes, r.energy, r.balance, r.contactForce, r.malletPos].filter((a) => a).map((a) => a!.buffer) as ArrayBuffer[];
1415
1516 self.onmessage = async (ev: MessageEvent<WorkerRequest>) => {
16- const { sys, exc, opt, gpu } = ev.data;
17+ const { sys, exc, gpu } = ev.data;
18+ const out = ev.data.opt;
19+ // solve at R times the sample rate; the pickup is filtered back down to fs (antialias.ts)
20+ const R = oversampling(sys, exc, out);
21+ const opt = oversampledOptions(out, R);
1722 const progress = (f: number) => self.postMessage({ type: 'progress', fraction: f });
18- const partial = (result: SimResult) => self.postMessage({ type: 'partial', result }, { transfer: buffers(result) });
23+ const partial = (r: SimResult) => {
24+ const result = decimateResult(r, R, false);
25+ self.postMessage({ type: 'partial', result }, { transfer: buffers(result) });
26+ };
1927 try {
2028 let result: SimResult | undefined;
2129 let fallback: string | undefined;
@@ -28,6 +36,7 @@ self.onmessage = async (ev: MessageEvent<WorkerRequest>) => {
2836 }
2937 }
3038 result ??= { ...simulate(sys, exc, opt, progress, (p) => partial({ ...p, backend: 'CPU', fallback })), backend: 'CPU', fallback };
39+ result = decimateResult(result, R, true);
3140 self.postMessage({ type: 'done', result }, { transfer: buffers(result) });
3241 } catch (e) {
3342 self.postMessage({ type: 'error', message: String(e) });
src/ui/instrumentPage.tsmodified+2−1View file
@@ -339,7 +339,8 @@ export function mountInstrumentPage<P extends SimParams>(root: HTMLElement, cfg:
339339 [`Peak ${quantity.label}`, `${fmtNum(peak)} ${quantity.unit}`],
340340 ]);
341341 setAudio(encodeWav(r.audio, r.fs), update);
342- q('audio-note').textContent = `16-bit mono, ${r.fs} Hz, normalized so the peak is −1 dBFS${partial ? ' (of the sound so far)' : ''}.`;
342+ const aa = (r.oversampling ?? 1) > 1 ? ` Anti-aliased: the pickup was computed at ${r.oversampling}× the sample rate and low-pass filtered below ${fmtNum(r.fs / 2000)} kHz.` : '';
343+ q('audio-note').textContent = `16-bit mono, ${r.fs} Hz, normalized so the peak is −1 dBFS${partial ? ' (of the sound so far)' : ''}.${aa}`;
343344 drawWave();
344345 spec.setSignal(r.audio, r.fs);
345346 drawSpectrum();
src/ui/settings.tsmodified+8−3View file
@@ -1,4 +1,5 @@
11 import { timeStep } from '../sim/solver';
2+import { oversampling } from '../sim/antialias';
23 import type { SimOptions } from '../sim/types';
34 import type { Built, InstrumentConfig, SimParams } from './instrumentPage';
45
@@ -48,11 +49,15 @@ export function simOptions(p: SimParams, built: Built): SimOptions {
4849 return { duration: p.duration, fs: p.fs, cfl: p.cfl, frameRate: p.frameRate, pickup: built.pickup, pickupDelay: built.pickupDelay };
4950 }
5051
51-// Size of a run: time step and element updates (a bending hinge costs about four
52+// Size of a run: time step (the solver runs at R times fs, see antialias.ts),
53+// steps per output sample and element updates (a bending hinge costs about four
5254 // springs). A desktop CPU does roughly CPU_RATE updates per second.
5355 export const CPU_RATE = 2.2e7;
5456 export function runCost(p: SimParams, built: Built) {
55- const { omegaMax, stepsPerSample } = timeStep(built.sys, built.exc.mallet, p.fs, p.cfl);
57+ const opt = simOptions(p, built);
58+ const R = oversampling(built.sys, built.exc, opt);
59+ const { omegaMax, stepsPerSample: inner } = timeStep(built.sys, built.exc.mallet, p.fs * R, p.cfl);
60+ const stepsPerSample = inner * R;
5661 const elements = built.sys.si.length + 4 * (built.sys.bending?.w0.length ?? 0);
57- return { omegaMax, stepsPerSample, work: stepsPerSample * p.fs * p.duration * elements };
62+ return { omegaMax, stepsPerSample, oversampling: R, work: stepsPerSample * p.fs * p.duration * elements };
5863 }
test/antialias.test.tsadded+56−0View file
@@ -0,0 +1,56 @@
1+import { describe, expect, it } from 'vitest';
2+import { decimateResult, lowpass } from '../src/sim/antialias';
3+import type { SimResult } from '../src/sim/types';
4+
5+const fs = 48000;
6+
7+// a result holding only a pickup signal at R fs
8+function result(x: Float32Array, R: number): SimResult {
9+ return {
10+ audio: x,
11+ fs: fs * R,
12+ frames: new Float32Array(0),
13+ frameTimes: new Float64Array(0),
14+ energy: new Float64Array(0),
15+ balance: new Float64Array(0),
16+ dt: 1 / (fs * R),
17+ stepsPerSample: 1,
18+ omegaMax: 0,
19+ wallSeconds: 0,
20+ maxStretch: 0,
21+ };
22+}
23+
24+const sine = (f: number, R: number, n: number) => Float32Array.from({ length: n * R }, (_, i) => Math.sin((2 * Math.PI * f * i) / (fs * R)));
25+
26+// RMS of the middle half (away from the edges)
27+function rms(x: Float32Array): number {
28+ const a = Math.floor(x.length / 4);
29+ const b = Math.floor((3 * x.length) / 4);
30+ let s = 0;
31+ for (let i = a; i < b; i++) s += x[i] * x[i];
32+ return Math.sqrt(s / (b - a));
33+}
34+
35+describe('anti-aliased pickup', () => {
36+ for (const R of [2, 5]) {
37+ it(`R = ${R}: passes the audible band, removes what would fold back`, () => {
38+ const h = lowpass(R);
39+ expect(h.reduce((s, v) => s + v, 0)).toBeCloseTo(1, 12);
40+ const n = 8000;
41+ // 1 kHz and 20 kHz pass (to 0.1 dB)
42+ for (const f of [1000, 20000]) expect(Math.abs(20 * Math.log10(rms(decimateResult(result(sine(f, R, n), R), R, true).audio) / Math.SQRT1_2))).toBeLessThan(0.1);
43+ // 26 and 36 kHz would fold back to 22 and 12 kHz: at least 80 dB down
44+ for (const f of [26000, 36000]) expect(20 * Math.log10(rms(decimateResult(result(sine(f, R, n), R), R, true).audio) / Math.SQRT1_2)).toBeLessThan(-80);
45+ });
46+ }
47+
48+ it('a partial result is a prefix of the complete one', () => {
49+ const R = 4;
50+ const x = Float32Array.from({ length: 6000 * R }, (_, i) => Math.sin(i * 0.01) * Math.exp(-i / 9000) + 0.3 * Math.sin(i * 1.7));
51+ const full = decimateResult(result(x, R), R, true).audio;
52+ const part = decimateResult(result(x.slice(0, 3000 * R), R), R, false).audio;
53+ expect(part.length).toBeGreaterThan(2800);
54+ expect(Array.from(part)).toEqual(Array.from(full.subarray(0, part.length)));
55+ });
56+});