Note renderer: step at the stability step and read the signal between steps by cubic interpolation (no longer at least 4 steps per audio sample); plate notes 1.4-2.3x faster - listening-for-chaos/chaos-apps

listening-for-chaos / chaos-apps

Note renderer: step at the stability step and read the signal between steps by cubic interpolation (no longer at least 4 steps per audio sample); plate notes 1.4-2.3x faster

Jeremy Magland <jmagland@flatironinstitute.org> committed 9 days ago commit 0394e117a2dd parent ad9866b Browse files

2 changed files+50−16

src/common/note-render.tsmodified+49−15View file

@@ -1,11 +1,21 @@

11 // Renders a note in a worker: a struck system run from a given state, sped

22 // up, sampled as audio. Each audio sample advances the physics by

3-// speedup / sampleRate seconds, in equal substeps no longer than the

3+// speedup / sampleRate seconds. The signal is read OVERSAMPLE times per

4+// audio sample, then low-passed and decimated (./antialias.ts), so motion

5+// above half the sample rate does not alias. The integrator step is the

46 // system's step for the current energy (re-chosen every DIAG_EVERY samples,

5-// so the steps lengthen as a damped note decays). The signal is read

6-// OVERSAMPLE times per audio sample (the substeps per sample are a multiple

7-// of it), then low-passed and decimated (./antialias.ts), so motion above

8-// half the sample rate does not alias. On request a Lyapunov

7+// so the steps lengthen as a damped note decays): when that is shorter than

8+// the interval between reading points, an equal division of it, so the

9+// steps land on the reading points; when longer, the step itself, and the

10+// signal is read between steps by cubic interpolation through the four

11+// nearest steps (Lagrange). A step resolves the fastest motion (omega h <=

12+// 0.2 for the stepFor of every app), so the interpolation error is about

13+// 0.02 (omega h)^4 of a mode's amplitude, -88 dB at the top mode and far

14+// less for lower ones; and the long steps occur only when the highest mode

15+// is below about 6 kHz in sound (0.2 x OVERSAMPLE x rate / 2 pi), far under

16+// the filter's passband edge. Until 2026-10-09 the steps were never longer

17+// than the reading interval, which cost up to several times more steps for

18+// notes whose highest mode is low. On request a Lyapunov

919 // shadow runs along (LYAP_WINDOW of physical time; it doubles the cost). The

1020 // energy and the windowed exponent (NaN without the shadow) are recorded

1121 // once per DIAG_EVERY samples. A damped note ends once its energy has fallen

@@ -26,7 +36,7 @@ export const LYAP_WINDOW = 10; // s of physical time

2636 export interface NoteSystem {

2737 energy(): number;

2838 stepFor(E: number): number; // integrator step (s) at energy E

29- step(h: number, n: number): void; // advance by h during sample n

39+ step(h: number, n: number): void; // advance by h; n is the sample at the end of the step

3040 signal(): number; // value heard now

3141 shadow(): Shadow; // a fresh Lyapunov shadow of the current state

3242 diag?(k: number): void; // called at diagnostic sample k (sample k DIAG_EVERY)

@@ -47,17 +57,37 @@ export interface Note {

4757 }

4858

4959 export function renderNote(sys: NoteSystem, o: NoteOptions): Note {

50- const dtSample = o.speedup / o.sampleRate;

51- let h = dtSample;

60+ const dtSample = o.speedup / o.sampleRate, dtRead = dtSample / OVERSAMPLE;

61+ let h = dtRead;

5262 const E0 = sys.energy();

5363 const nMax = Math.round((o.damped ? o.maxLength : UNDAMPED_LENGTH) * o.sampleRate);

5464 const over = new Float32Array(OVERSAMPLE * nMax); // signal at OVERSAMPLE points per sample

55- let per = 1; // substeps per oversampled point

5665 const en = new Float32Array(Math.ceil(nMax / DIAG_EVERY));

5766 const lam = new Float32Array(en.length);

5867 const lyap = new Lyapunov(LYAP_WINDOW);

59- let t = 0; // physical time (s)

68+ let t = 0; // physical time (s) of the system's state

6069 if (o.lyapunov) lyap.start(sys.shadow(), t);

70+ // The last four steps' times and signal values (a ring), for reading

71+ // between them.

72+ const ht = new Float64Array(4), hy = new Float64Array(4);

73+ let hn = 0; // values pushed so far

74+ const push = (tt: number, y: number) => {

75+ ht[hn & 3] = tt;

76+ hy[hn & 3] = y;

77+ hn++;

78+ };

79+ push(0, sys.signal());

80+ // Signal at time tau, from the last (up to) four steps around it.

81+ const read = (tau: number): number => {

82+ const k0 = Math.max(0, hn - 4);

83+ let y = 0;

84+ for (let i = k0; i < hn; i++) {

85+ let w = 1;

86+ for (let j = k0; j < hn; j++) if (j !== i) w *= (tau - ht[j & 3]) / (ht[i & 3] - ht[j & 3]);

87+ y += w * hy[i & 3];

88+ }

89+ return y;

90+ };

6191 let n = 0;

6292 for (; n < nMax; n++) {

6393 if (n % DIAG_EVERY === 0) {

@@ -67,16 +97,20 @@ export function renderNote(sys: NoteSystem, o: NoteOptions): Note {

6797 lam[k] = o.lyapunov ? lyap.windowed() : NaN;

6898 sys.diag?.(k);

6999 if (o.damped && E <= END_RATIO * E0) break;

70- per = Math.ceil(dtSample / (OVERSAMPLE * sys.stepFor(E)));

71- h = dtSample / (OVERSAMPLE * per);

100+ const step = sys.stepFor(E);

101+ h = step >= dtRead ? step : dtRead / Math.ceil(dtRead / step);

72102 }

73103 for (let q = 0; q < OVERSAMPLE; q++) {

74- over[OVERSAMPLE * n + q] = sys.signal();

75- for (let i = 0; i < per; i++) {

76- sys.step(h, n);

104+ const tau = (OVERSAMPLE * n + q) * dtRead;

105+ // Step until the second-to-last step is at or past tau, so tau lies

106+ // between the middle two of the last four.

107+ while (hn < 2 || ht[(hn - 2) & 3] < tau) {

108+ sys.step(h, Math.floor((t + h) / dtSample));

77109 t += h;

78110 if (o.lyapunov) lyap.step(h, t);

111+ push(t, sys.signal());

79112 }

113+ over[OVERSAMPLE * n + q] = read(tau);

80114 }

81115 }

82116 const nDiag = Math.ceil(n / DIAG_EVERY);

src/plate/README.mdmodified+1−1View file

@@ -24,6 +24,6 @@ Clamped edges (step 5, 2026-10-09). Edges select: simply supported, clamped, fre

2424

2525 A membrane tension T for held edges (drum-like: springs between neighbours and to the supports, energy T/2 uᵀLu, exact sine modes with ω² = (T l + κ l²)/m when simply supported) was added in step 5 and removed the same day: the suggested drum setting (κ 2, T 200, h 1 mm) did not sound like a drum. Drums to be revisited later.

2626

27-Known cost issue (2026-10-09): the 4x oversampled note forces at least 4 steps per audio sample, and the step stability needs is longer than that whenever the speed-up is small, that is when the lowest mode is close in frequency to the top of the band: clamped (lowest mode higher; a tensioned membrane, tried briefly, was far worse, about 8x more steps than needed). Clamped default note: 12 s to render in the browser. Fix planned: read the signal between steps by interpolation (shared note renderer).

27+Cost issue, fixed (2026-10-09, `../common/note-render.ts`): the 4x oversampled note used to force at least 4 steps per audio sample, more than stability needs whenever the highest mode is low in sound (below about 6 kHz at 48 kHz: clamped above all). The renderer now steps at the stability step and reads the signal between steps by cubic interpolation through the four nearest steps. Checks: a linear undamped note (simply supported, 1 J, 10 s) has nothing above −140 dB away from the modes either way (old −144 dB, new −141 dB), with 1.9x fewer steps; damped nonlinear notes (2.51 J, h 5 cm) keep level and spectral centroid over time within the spread a 10% change of step gives anyway (chaotic notes: a different step gives a different, equivalent realization; for clamped, up to 2.7 dB and 100 Hz from a 10% step change, more than old against new). Browser, same conditions, default note: simply supported 7.2 → 5.1 s, clamped 11.9 → 5.2 s, free 9.3 → 6.1 s. Beam and FPUT chain unchanged in practice (their top modes sit above 6 kHz at the default pitch, so their steps were already shorter than the reading interval).

2828

2929 Lyapunov exponent (step 6, 2026-10-09): as the beam's (two-trajectory method, `../common/lyapunov.ts`), with the plate's shadow in `lyapunov.ts`: a second plate D0 = 1e-8 away, distances with displacements in units of A/√(N²) (A the lowest mode's amplitude at which stretching equals bending, `nonlinearAmplitude`; 1 m for a linear plate) and momenta in m ω₁ times that; the initial offset is a fixed pseudo-random pattern (free: its rigid part removed). Layout now as the beam's: Parameters folded under the plate, a Real time card on the right with Lyapunov (λ 10 s and since the strike, first 5 s left out of the y range) and Energy error charts. Sound panel: Lyapunov chart and λ vs E scatter, both off at first (the shadow doubles the time to compute a note: default note 4.3 s without, about 12 s with, in the browser). Exponents (2026-10-09, N = 8, kick at (1.3, 1.1), undamped, 300 s, λ since the strike, 1/s): linear plate 0.007 at every energy (finite-time; the last 100 s, −0.001: zero); simply supported, h 5 cm: 0.008, 0.065, 0.26, 0.53, 1.01 at 0.001, 0.1, 1, 10, 100 J; free, h 5 cm: 0.006, 0.015, 0.16, 1.54, 8.4. The free plate's exponent grows much faster with energy.