Whole-note sound displays with Lyapunov exponent and energy; strike jitter - listening-for-chaos/double-pendulum-chaos

listening-for-chaos / double-pendulum-chaos

Whole-note sound displays with Lyapunov exponent and energy; strike jitter

- Faster stage solve: extrapolated initial guess, tolerance 1e-12 - Audio worker measures the note's Lyapunov exponent (shadow method) and energy; shown under the sound spectrogram on the same time axis - Sound spectrogram and charts show the whole note without scrolling - Energy charts start at 0; log charts use whole decades - Tiny random jitter (1e-9) on every strike

Jeremy Magland <jeremy.magland@gmail.com> committed 10 days ago commit 30a92c8ecfb3 parent 954de9b Browse files

10 changed files+372−48

README.mdmodified+4−4View file

@@ -2,7 +2,7 @@

22

33 Interactive double pendulum, step 1 of "Listening for chaos". Work in progress.

44

5-Current state: the pendulum hangs at rest until Strike, which starts it from hanging rest with the energy on the slider (default 50 J); Pause freezes the animation. The Ghost toggle adds a second pendulum (translucent orange) started from the same state with the outer arm's angle offset by 10^-12 to 10^-1 rad (default 10^-10), restarted on Strike or Reset, with a readout of the angular separation and the time since the ghost started: it grows exponentially when the motion is chaotic and stays tiny when it is regular. Charts "Separation" (ghost separation, log scale, last 30 s) and "Lyapunov" (largest Lyapunov exponent over a 10 s window and since the strike, last 30 s). Parameters (masses, lengths, g, decay time, presets) fold out underneath. Toggle chips choose the plots: a spectrogram of the selected signal (frequency axis auto-scaled, dashed lines at the normal modes) and live traces of energies, energy balance error, angles, momenta and bob positions. The sound card at the bottom has a Sound on/off toggle (with sound on, every Strike plays a note: a separate, sped-up copy of the motion from the same strike; a damped note ends once it has decayed by 60 dB), "Sound of" (the signal, also used by the animation spectrogram), volume, pitch (where the small-amplitude lower mode lands, 30 to 1000 Hz, default 100 Hz), and a spectrogram of the sound.

5+Current state: the pendulum hangs at rest until Strike, which starts it from hanging rest with the energy on the slider (default 50 J), plus a random relative jitter of up to 1e-9 in each angle and momentum (so chaotic strikes differ from one another while regular ones sound alike); Pause freezes the animation. The Ghost toggle adds a second pendulum (translucent orange) started from the same state with the outer arm's angle offset by 10^-12 to 10^-1 rad (default 10^-10), restarted on Strike or Reset, with a readout of the angular separation and the time since the ghost started: it grows exponentially when the motion is chaotic and stays tiny when it is regular. Charts "Separation" (ghost separation, log scale, last 30 s) and "Lyapunov" (largest Lyapunov exponent over a 10 s window and since the strike, last 30 s). Parameters (masses, lengths, g, decay time, presets) fold out underneath. Toggle chips choose the plots: a spectrogram of the selected signal (frequency axis auto-scaled, dashed lines at the normal modes) and live traces of energies, energy balance error, angles, momenta and bob positions. The sound card at the bottom has a Sound on/off toggle (with sound on, every Strike plays a note: a separate, sped-up copy of the motion from the same strike; a damped note ends once it has decayed by 60 dB), "Sound of" (the signal, also used by the animation spectrogram), volume, pitch (where the small-amplitude lower mode lands, 30 to 1000 Hz, default 100 Hz), and optional whole-note displays that do not scroll (strike at the left edge, end of the note at the right): a spectrogram of the sound, the largest Lyapunov exponent of the motion being heard (10 s window, physical time), and its energy (log scale). The time axis spans the predicted note length ((tau/2) ln 1e6 of physical time, divided by the speedup), is extended if the note runs longer, and is fitted to the note when it ends.

66

77 ```bash

88 npm install

@@ -17,8 +17,8 @@ npm run dev

1717 - Damping: linear drag on each bob, force -gamma m_i v_i (proportional to its mass; equal for equal bobs), so the generalized forces are -gamma p and every mode's amplitude decays with time constant tau = 2/gamma. Default tau = 540 s: about 10% amplitude per minute in the animation, about 2 s in the sound at x262. The state carries the dissipated energy W, so E + W is checked for conservation. With damping the energy slider sets the strike energy; with tau = 0 it sets the energy live.

1818 - Live energy change: both momenta are scaled by one factor at fixed angles (a push along the current motion). A target below the current potential energy can't be reached at once; until it can, kinetic energy is removed at a rate of 6 per second.

1919 - Lyapunov exponent: two-trajectory (Benettin) method in `src/lyapunov.ts`. A hidden shadow pendulum is kept 1e-8 away (angles in rad, momenta in units of sqrt((m1 + m2) l1^2 E_top)), moved back after every step, and the log growth averaged. The shadow keeps its own energy (forcing it onto the pendulum's energy by rescaling momenta inflates the separation near turning points). Checked (no damping, 300 s): 2 J and 10 J fall toward 0 (0.01/s), 30 J ~1.0/s, 50 J ~1.4/s; unchanged for shadow distances 1e-6 and 1e-10. With damping (tau 540 s) a 50 J strike's 10 s value falls from ~1.4 to ~0.1/s by 300 s.

20-- Integrator: 2-stage Gauss-Legendre (order 4, symplectic), fixed step h chosen per parameter set by a short calibration run so the relative energy error stays below 1e-6 up to the top of the energy slider (twice the both-inverted energy); 3.13 ms for the defaults, fixed-point stage solve. Angles are wrapped to (-pi, pi] after each step. The physics costs well under 1 ms of CPU per second; drawing dominates, so traces draw at most a min/max pair per pixel column and the spectrogram scrolls in new columns instead of redrawing.

20+- Integrator: 2-stage Gauss-Legendre (order 4, symplectic), fixed step h chosen per parameter set by a short calibration run so the relative energy error stays below 1e-6 up to the top of the energy slider (twice the both-inverted energy); 3.13 ms for the defaults, fixed-point stage solve (started from the previous step's stage derivatives extrapolated in time; tolerance 1e-12, below which the energy error does not change). Angles are wrapped to (-pi, pi] after each step. The physics costs well under 1 ms of CPU per second; drawing dominates, so traces draw at most a min/max pair per pixel column and the spectrogram scrolls in new columns instead of redrawing.

2121

22-Audio: a Web Worker (`src/audio-worker.ts`) integrates the audio copy with substeps of at most 3.2 ms of simulated time per sample, high-passes at 20 Hz, and streams 1024-sample chunks through a MessageChannel to an AudioWorklet (`src/audio-worklet.ts`) that plays them from a ring buffer, keeping about 0.25 s queued. The page shows how many times faster than real time the worker computes, and any dropouts. The worker also sends each chunk to the page, which feeds it to the audio spectrogram as the audio clock reaches it (2048-sample window, new column every 512 samples).

22+Audio: a Web Worker (`src/audio-worker.ts`) integrates the audio copy with substeps of at most 3.2 ms of simulated time per sample, high-passes at 20 Hz, and streams 1024-sample chunks through a MessageChannel to an AudioWorklet (`src/audio-worklet.ts`) that plays them from a ring buffer, keeping about 0.25 s queued. The page shows how many times faster than real time the worker computes, and any dropouts. The worker also sends each chunk to the page, which feeds it to the audio spectrogram as the audio clock reaches it (2048-sample window, new column every 512 samples). The worker runs the same Lyapunov shadow on the note and sends one exponent and one energy value per 256 samples with each chunk; the page plots them in step with the spectrogram. The whole-note spectrogram keeps every column (pairs are averaged past 4096 columns) and draws only new pixel columns unless the axis or color scale changes.

2323

24-Code: `src/pendulum.ts` (physics), `src/sim.ts` (running simulation with energy target), `src/audio.ts` (audio setup), `src/signals.ts` (audible quantities), `src/lyapunov.ts` (Lyapunov exponent), `src/traces.ts` (strip charts), `src/spectrogram.ts` (spectrogram), `src/main.ts` (page).

24+Code: `src/pendulum.ts` (physics), `src/sim.ts` (running simulation with energy target), `src/audio.ts` (audio setup), `src/signals.ts` (audible quantities), `src/lyapunov.ts` (Lyapunov exponent), `src/audio-shared.ts` (constants shared with the worker), `src/traces.ts` (strip charts), `src/spectrogram.ts` (spectrogram), `src/main.ts` (page).

index.htmlmodified+6−2View file

@@ -61,6 +61,7 @@

6161 button.sound.on { background: var(--accent-soft); border-color: #b9cde0; color: var(--accent); }

6262

6363 #controls { margin-top: 14px; }

64+ #energy { min-width: 70px; }

6465 #energy-value { min-width: 5.5em; text-align: right; }

6566 #ghost-row { margin-top: 10px; min-height: 28px; }

6667 #ghost-controls { flex: 1; }

@@ -92,7 +93,7 @@

9293 figure { margin: 14px 0 0; }

9394 figure canvas { width: 100%; height: auto; border-radius: 4px; }

9495 figcaption { margin-top: 4px; }

95- #traces canvas { width: 100%; height: auto; margin-top: 12px; }

96+ #traces canvas, #audio-traces canvas { width: 100%; height: auto; margin-top: 12px; }

9697 #plots-empty { margin: 14px 0 2px; }

9798

9899 #audio { margin-top: 16px; }

@@ -174,7 +175,9 @@

174175 <button id="sound" type="button" class="sound" title="With sound on, every Strike plays a note">Sound off</button>

175176 <label><span class="label">Sound of</span> <select id="signal"></select></label>

176177 <label><span class="label">Volume</span> <input id="volume" type="range" min="0" max="1" step="0.01" value="0.5" /></label>

177- <label class="chip"><input id="audio-spec-toggle" type="checkbox" /><span>Spectrogram</span></label>

178+ <span id="audio-chips" class="chips">

179+ <label class="chip"><input id="audio-spec-toggle" type="checkbox" /><span>Spectrogram</span></label>

180+ </span>

178181 </div>

179182 <div class="row">

180183 <label class="grow"><span class="label">Pitch</span> <input id="pitch" type="range" min="0" max="1000" step="1" aria-label="Pitch of the lower mode" /></label>

@@ -185,6 +188,7 @@

185188 <div id="audio-spectrogram"></div>

186189 <figcaption id="audio-spec-label" class="muted"></figcaption>

187190 </figure>

191+ <div id="audio-traces"></div>

188192 </section>

189193 </main>

190194 <script type="module" src="/src/main.ts"></script>

src/audio-shared.tsadded+7−0View file

@@ -0,0 +1,7 @@

1+// Constants shared by the audio worker and the page.

2+

3+// Samples per diagnostic value (Lyapunov exponent, energy) sent with each

4+// chunk; divides the chunk size.

5+export const LAMBDA_EVERY = 256;

6+

7+export const END_RATIO = 1e-6; // energy ratio at which a damped note ends (-60 dB amplitude)

src/audio-worker.tsmodified+21−5View file

@@ -5,9 +5,13 @@

55 // seconds queued ahead of playback. Also sends a copy of each chunk to the

66 // page (for the audio spectrogram) and reports how fast it computes relative

77 // to real time. With damping, a note ends (and computing stops) once its

8-// energy has fallen by END_RATIO (-60 dB in amplitude).

8+// energy has fallen by END_RATIO (-60 dB in amplitude). The note's largest

9+// Lyapunov exponent (10 s window, physical time) is measured alongside; it and

10+// the energy are sent with the samples, one value per LAMBDA_EVERY samples.

911

1012 import { energy, normalModes, stepSize, topEnergy, type Params } from "./pendulum";

13+import { END_RATIO, LAMBDA_EVERY } from "./audio-shared";

14+import { Lyapunov } from "./lyapunov";

1115 import { signalByKey, type Signal } from "./signals";

1216 import { Sim } from "./sim";

1317

@@ -15,7 +19,6 @@ const CHUNK = 1024; // samples per message

1519 const BUFFER = 0.25; // seconds of audio kept queued

1620 const HIGHPASS_HZ = 20; // DC blocker cutoff

1721 const REPORT_MS = 500;

18-const END_RATIO = 1e-6; // energy ratio at which a damped note ends

1922

2023 export type ToWorker =

2124 | { type: "init"; P: Params; alpha: number; sampleRate: number; port: MessagePort }

@@ -29,7 +32,7 @@ export type ToWorker =

2932 export type FromWorker =

3033 | { type: "stats"; speed: number; underruns: number }

3134 | { type: "ended"; session: number }

32- | { type: "samples"; session: number; data: Float32Array };

35+ | { type: "samples"; session: number; data: Float32Array; lambda: Float32Array; energy: Float32Array };

3336

3437 let sim: Sim;

3538 let port: MessagePort;

@@ -43,6 +46,7 @@ let session = 0;

4346 let running = false;

4447 let produced = 0, consumed = 0, underruns = 0;

4548 let eStart = 0; // energy at the start of the note

49+const lyap = new Lyapunov(10);

4650 let computeMs = 0, computedSamples = 0, lastReport = 0;

4751 // One-pole DC blocker.

4852 let hpR = 0.997, hpX = 0, hpY = 0;

@@ -83,6 +87,7 @@ onmessage = (e: MessageEvent<ToWorker>) => {

8387 session++;

8488 sim.setState(msg.y, msg.E0, msg.target);

8589 eStart = energy(sim.P, sim.y);

90+ lyap.start(sim);

8691 produced = consumed = underruns = 0;

8792 computeMs = computedSamples = 0;

8893 lastReport = performance.now();

@@ -115,19 +120,30 @@ function pump(): void {

115120 const t0 = performance.now();

116121 while (produced - consumed < BUFFER * fs) {

117122 const chunk = new Float32Array(CHUNK);

123+ const lambda = new Float32Array(CHUNK / LAMBDA_EVERY);

124+ const noteEnergy = new Float32Array(CHUNK / LAMBDA_EVERY);

118125 // Simulated time per sample; the speedup follows the current parameters.

119126 const dt = pitchHz / fLower / fs;

120127 const n = Math.ceil(dt / hMax);

121128 const h = dt / n;

122129 for (let i = 0; i < CHUNK; i++) {

123- for (let k = 0; k < n; k++) sim.step(h);

130+ for (let k = 0; k < n; k++) {

131+ sim.step(h);

132+ lyap.step(sim, h);

133+ }

134+ if ((i + 1) % LAMBDA_EVERY === 0) {

135+ lambda[(i + 1) / LAMBDA_EVERY - 1] = lyap.windowed();

136+ noteEnergy[(i + 1) / LAMBDA_EVERY - 1] = energy(sim.P, sim.y);

137+ }

124138 const x = value();

125139 hpY = x - hpX + hpR * hpY;

126140 hpX = x;

127141 chunk[i] = gain * hpY;

128142 }

129143 const copy = chunk.slice();

130- postMessage({ type: "samples", session, data: copy } satisfies FromWorker, { transfer: [copy.buffer] });

144+ postMessage({ type: "samples", session, data: copy, lambda, energy: noteEnergy } satisfies FromWorker, {

145+ transfer: [copy.buffer, lambda.buffer, noteEnergy.buffer],

146+ });

131147 port.postMessage({ session, data: chunk }, [chunk.buffer]);

132148 produced += CHUNK;

133149 computedSamples += CHUNK;

src/audio.tsmodified+4−2View file

@@ -15,7 +15,9 @@ type Stats = Extract<FromWorker, { type: "stats" }>;

1515

1616 export interface AudioCallbacks {

1717 onStats: (s: Stats) => void;

18- onSamples: (data: Float32Array) => void; // samples of the current note

18+ // Samples of the current note, with its Lyapunov exponent (1/s, physical)

19+ // and energy (J), one value each per LAMBDA_EVERY samples.

20+ onSamples: (data: Float32Array, lambda: Float32Array, energy: Float32Array) => void;

1921 onEnded: () => void; // the current note decayed and stopped

2022 }

2123

@@ -66,7 +68,7 @@ export class AudioOut {

6668 const msg = e.data;

6769 if (msg.type === "stats") this.cb.onStats(msg);

6870 else if (msg.session !== this.session || !this.playing) return;

69- else if (msg.type === "samples") this.cb.onSamples(msg.data);

71+ else if (msg.type === "samples") this.cb.onSamples(msg.data, msg.lambda, msg.energy);

7072 else if (msg.type === "ended") {

7173 this.playing = false;

7274 this.cb.onEnded();

src/main.tsmodified+69−17View file

@@ -1,14 +1,16 @@

11 import { AudioOut } from "./audio";

2+import { END_RATIO, LAMBDA_EVERY } from "./audio-shared";

23 import { Lyapunov } from "./lyapunov";

34 import {

45 energy, kinetic, normalModes, potential, stepSize, topEnergy, type Params,

56 } from "./pendulum";

67 import { SIGNALS, signalByKey } from "./signals";

78 import { Sim } from "./sim";

8-import { Spectrogram } from "./spectrogram";

9+import { NoteSpectrogram, Spectrogram } from "./spectrogram";

910 import { TracePanel } from "./traces";

1011

1112 const ALPHA = 0.5; // fraction of energy in the lower mode

13+const STRIKE_JITTER = 1e-9; // random relative perturbation of each strike (see Sim.strike)

1214 // Parameters; changed with the parameter form (Apply restarts everything).

1315 // tau: amplitude decay time from drag on the bobs; 540 s barely shows in the

1416 // animation (~10% amplitude per minute) but makes each struck note ring down

@@ -101,6 +103,7 @@ const traces = new TracePanel(

101103 {

102104 label: "Energy",

103105 unit: "J",

106+ zero: true,

104107 series: [

105108 { label: "E", value: () => energy(P, sim.y) },

106109 { label: "T", value: () => kinetic(P, sim.y) },

@@ -179,41 +182,81 @@ const showPitch = () => {

179182 `${lo.toFixed(0)} Hz · upper ${hi.toFixed(0)} Hz · ×${S.toFixed(0)}` +

180183 (damped() ? ` · decay ${(P.tau / S).toFixed(1)} s` : "");

181184 audioSpecLabel.textContent =

182- `Spectrogram of the sound (${signal.label}, sped up ×${(lo / modeHz[0]).toFixed(0)}). Dashed: normal modes.`;

185+ `Spectrogram of the whole note, strike (left) to end (right): ${signal.label}, sped up ×${(lo / modeHz[0]).toFixed(0)}. Dashed: normal modes.`;

183186 };

184187 showPitch();

185188

186-// Spectrogram of the audio, fed with the worker's samples as the audio clock

187-// reaches them (2048-sample window, new column every 512 samples, 960

188-// columns: ~10 s at 48 kHz). Rebuilt if the audio device's sample rate

189-// differs.

189+// Displays of the current note, from strike (left) to end (right), fed with

190+// the worker's samples as the audio clock reaches them. The time axis spans

191+// the predicted note length (energy decays at about 2/tau, so a drop by

192+// END_RATIO takes ~(tau/2) ln(1/END_RATIO) of physical time), is extended

193+// if the note runs longer, and fitted to the note when it ends.

194+// Spectrogram: 2048-sample window, new column every 512 samples. Rebuilt if

195+// the audio device's sample rate differs.

190196 const makeAudioSpec = (rate: number) => {

191- const s = new Spectrogram({ dt: 1 / rate, n: 2048, nfft: 4096, hop: 512, refFreqs: audioModeHz(), width: 960 });

197+ const s = new NoteSpectrogram({ dt: 1 / rate, n: 2048, nfft: 4096, hop: 512, refFreqs: audioModeHz(), width: 960 });

192198 document.getElementById("audio-spectrogram")!.replaceChildren(s.canvas);

193199 return s;

194200 };

201+const predictedNoteLength = () => {

202+ const S = pitchHz() / modeHz[0];

203+ return damped() ? ((P.tau / 2) * Math.log(1 / END_RATIO)) / S : 10;

204+};

205+let noteSpan = 10; // s of sound across the displays

206+let received = 0; // samples received for the current note

207+const setNoteSpan = (span: number) => {

208+ noteSpan = span;

209+ audioSpec.setSpan(span);

210+ audioTraces.setFixed({ t0: 0, span });

211+};

195212 let audioSpecRate = 48000;

196213 let audioSpec = makeAudioSpec(audioSpecRate);

197214 const audioSpecToggle = document.getElementById("audio-spec-toggle") as HTMLInputElement;

198215 const audioSpecFig = document.getElementById("audio-spec-fig")!;

199216 audioSpecToggle.addEventListener("change", () => (audioSpecFig.hidden = !audioSpecToggle.checked));

200-const queue: Float32Array[] = []; // sample chunks not yet shown

217+// Lyapunov exponent and energy of the note, under the spectrogram on the same

218+// time axis.

219+let noteLambda = NaN, noteEnergy = NaN;

220+const audioTraces = new TracePanel(

221+ document.getElementById("audio-chips")!,

222+ document.getElementById("audio-traces")!,

223+ [

224+ { label: "Lyapunov", unit: "1/s", series: [{ label: "λ 10 s", value: () => noteLambda }] },

225+ { label: "Energy", unit: "J", zero: true, series: [{ label: "E", value: () => noteEnergy }] },

226+ ],

227+ 10,

228+ 960,

229+);

230+audioTraces.setFixed({ t0: 0, span: noteSpan });

231+

232+// Chunks of the current note not yet shown, fed to the displays as the audio

233+// clock reaches them.

234+const queue: { data: Float32Array; lambda: Float32Array; energy: Float32Array }[] = [];

201235 let queueOffset = 0; // samples of queue[0] already shown

202236 let fed = 0; // samples shown this session

203237

204238 function feedAudioSpec() {

205239 const played = audio.playedSamples();

206240 while (fed < played && queue.length > 0) {

207- const chunk = queue[0];

208- const n = Math.min(chunk.length - queueOffset, Math.ceil(played - fed));

209- for (let i = 0; i < n; i++) audioSpec.push(chunk[queueOffset + i]);

241+ const { data, lambda, energy: en } = queue[0];

242+ const n = Math.min(data.length - queueOffset, Math.ceil(played - fed));

243+ for (let i = 0; i < n; i++) {

244+ const j = queueOffset + i;

245+ audioSpec.push(data[j]);

246+ if ((j + 1) % LAMBDA_EVERY === 0) {

247+ noteLambda = lambda[(j + 1) / LAMBDA_EVERY - 1];

248+ noteEnergy = en[(j + 1) / LAMBDA_EVERY - 1];

249+ audioTraces.sample((fed + i + 1) / audio.sampleRate);

250+ }

251+ }

210252 queueOffset += n;

211253 fed += n;

212- if (queueOffset === chunk.length) {

254+ if (queueOffset === data.length) {

213255 queue.shift();

214256 queueOffset = 0;

215257 }

216258 }

259+ if (audio.playing && fed / audio.sampleRate > noteSpan) setNoteSpan(1.25 * noteSpan);

217260 }

218261

219262 const audio = new AudioOut(P, ALPHA, {

@@ -224,16 +267,24 @@ const audio = new AudioOut(P, ALPHA, {

224267 (s.speed < 1.5 ? " (too slow for live audio)" : "") +

225268 (s.underruns > 0 ? ` · ${s.underruns} dropouts` : "");

226269 },

227- onSamples: (data) => queue.push(data),

228- onEnded: () => (audioStatus.textContent = "Note finished (decayed by 60 dB)."),

270+ onSamples: (data, lambda, energy) => {

271+ queue.push({ data, lambda, energy });

272+ received += data.length;

273+ },

274+ onEnded: () => {

275+ audioStatus.textContent = "Note finished (decayed by 60 dB).";

276+ setNoteSpan(received / audio.sampleRate);

277+ },

229278 });

230279

231280 // Start a note from the animation's current state.

232281 function startNote() {

233282 queue.length = 0;

234- queueOffset = fed = 0;

283+ queueOffset = fed = received = 0;

284+ audioTraces.clear();

235285 audio.play(sim.y, sim.E0, sim.target, pitchHz(), signal.key);

236- audioSpec.clear();

286+ audioSpec.start(predictedNoteLength());

287+ setNoteSpan(predictedNoteLength());

237288 audioStatus.textContent = "";

238289 }

239290

@@ -320,7 +371,7 @@ function showEnergy() {

320371

321372 document.getElementById("controls")!.addEventListener("submit", (e) => {

322373 e.preventDefault();

323- sim.strike(sliderEnergy());

374+ sim.strike(sliderEnergy(), STRIKE_JITTER);

324375 resetGhost();

325376 lyap.start(sim);

326377 setRunning(true);

@@ -497,6 +548,7 @@ function frame(now: number) {

497548 if (!specFig.hidden) spectrogram.draw();

498549 feedAudioSpec();

499550 if (!audioSpecFig.hidden) audioSpec.draw();

551+ audioTraces.draw();

500552 requestAnimationFrame(frame);

501553 }

502554 requestAnimationFrame(frame);

src/pendulum.tsmodified+28−5View file

@@ -137,9 +137,15 @@ export function strike(P: Params, E: number, alpha: number): Float64Array {

137137

138138 // 2-stage Gauss-Legendre (order 4; symplectic and time-symmetric when there

139139 // is no damping). The implicit stage equations are solved by fixed-point

140-// iteration.

140+// iteration, started from the previous step's stage derivatives extrapolated

141+// linearly in time (or from f(y) if y was changed from outside since the

142+// last step, or h changed).

141143 const S3 = Math.sqrt(3) / 6;

142144 const A11 = 0.25, A12 = 0.25 - S3, A21 = 0.25 + S3, A22 = 0.25;

145+// Stage times c1 = 1/2 - S3, c2 = 1/2 + S3. Extrapolating the previous step's

146+// stage derivatives (at c1, c2) linearly to this step's stage times

147+// (1 + c1, 1 + c2): K_i = K1 + (K2 - K1) * (1 + c_i - c1) / (c2 - c1).

148+const E1 = 1 / (2 * S3), E2 = (1 + 2 * S3) / (2 * S3);

143149

144150 export class GaussLegendre2 {

145151 private k1 = new Float64Array(STATE_SIZE);

@@ -150,6 +156,8 @@ export class GaussLegendre2 {

150156 private y2 = new Float64Array(STATE_SIZE);

151157 iterations = 0; // fixed-point iterations used by the last step

152158 private P: Params;

159+ private lastY = new Float64Array(STATE_SIZE); // state after the last step

160+ private lastH = NaN;

153161 private gamma: number; // drag rate

154162

155163 constructor(P: Params) {

@@ -158,9 +166,20 @@ export class GaussLegendre2 {

158166 }

159167

160168 step(y: Float64Array, h: number): void {

161- const { k1, k2, n1, n2, y1, y2, P, gamma } = this;

162- rhs(P, gamma, y, k1);

163- k2.set(k1);

169+ const { k1, k2, n1, n2, y1, y2, P, gamma, lastY } = this;

170+ let continued = h === this.lastH;

171+ for (let i = 0; continued && i < STATE_SIZE; i++) continued = y[i] === lastY[i];

172+ if (continued) {

173+ // k1, k2 still hold the previous step's stage derivatives.

174+ for (let i = 0; i < STATE_SIZE; i++) {

175+ const a = k1[i], d = k2[i] - k1[i];

176+ k1[i] = a + E1 * d;

177+ k2[i] = a + E2 * d;

178+ }

179+ } else {

180+ rhs(P, gamma, y, k1);

181+ k2.set(k1);

182+ }

164183 let it = 0;

165184 for (; it < 50; it++) {

166185 for (let i = 0; i < STATE_SIZE; i++) {

@@ -176,7 +195,9 @@ export class GaussLegendre2 {

176195 }

177196 k1.set(n1);

178197 k2.set(n2);

179- if (diff <= 1e-14 * scale) break;

198+ // 1e-12: the energy error (from the step size) is unchanged for

199+ // tolerances from 1e-14 to 1e-11 over 600 s at 50 and 100 J.

200+ if (diff <= 1e-12 * scale) break;

180201 }

181202 this.iterations = it + 1;

182203 for (let i = 0; i < STATE_SIZE; i++) y[i] += 0.5 * h * (k1[i] + k2[i]);

@@ -184,6 +205,8 @@ export class GaussLegendre2 {

184205 // angle without bound and lose precision in sin/cos.

185206 y[0] = wrap(y[0]);

186207 y[1] = wrap(y[1]);

208+ lastY.set(y);

209+ this.lastH = h;

187210 }

188211 }

189212

src/sim.tsmodified+14−1View file

@@ -35,8 +35,21 @@ export class Sim {

3535 this.target = E;

3636 }

3737

38- strike(E: number): void {

38+ // Strike from hanging rest with energy E. With jitter > 0 each angle and

39+ // momentum gets a random relative change of up to `jitter` (angles in rad,

40+ // momenta relative to the strike momentum), so chaotic strikes differ from

41+ // one another while regular ones stay alike; the first step then restores

42+ // the energy to exactly E.

43+ strike(E: number, jitter = 0): void {

3944 this.y = strike(this.P, E, this.alpha);

45+ if (jitter > 0) {

46+ const pScale = Math.hypot(this.y[2], this.y[3]);

47+ const r = () => jitter * (2 * Math.random() - 1);

48+ this.y[0] += r();

49+ this.y[1] += r();

50+ this.y[2] += r() * pScale;

51+ this.y[3] += r() * pScale;

52+ }

4053 this.E0 = energy(this.P, this.y);

4154 this.target = E;

4255 this.t = 0;

src/spectrogram.tsmodified+192−0View file

@@ -265,3 +265,195 @@ class FFT {

265265 }

266266 }

267267 }

268+

269+// Spectrogram of one whole note, from its start (left edge) to its end (right

270+// edge), without scrolling. The time axis spans `span` seconds, set at the

271+// start (predicted length) and changed with setSpan (extended if the note

272+// runs longer, fitted when it ends). All columns are kept; past MAX_COLS,

273+// adjacent columns are averaged in pairs (halving the time resolution) so

274+// memory stays bounded. Colors are relative to the loudest column so far, so

275+// the note's decay shows. Only new pixel columns are drawn, except when the

276+// axis, the frequency range, or the color scale changes.

277+const MAX_COLS = 4096;

278+

279+export class NoteSpectrogram {

280+ canvas: HTMLCanvasElement;

281+ private width: number;

282+ private dt: number;

283+ private n: number;

284+ private hop: number;

285+ private bins: number;

286+ private refFreqs: number[];

287+ private buf: Float64Array;

288+ private nSamples = 0;

289+ private hann: Float64Array;

290+ private fft: FFT;

291+ private re: Float64Array;

292+ private im: Float64Array;

293+ private cols: Float32Array[] = []; // stored columns (power)

294+ private merge = 1; // raw columns per stored column

295+ private pending: Float32Array; // raw columns being averaged into the next stored one

296+ private pendingCount = 0;

297+ private sum: Float64Array; // sum of all raw columns (time-averaged spectrum)

298+ private peak = 0;

299+ private fmax = 1;

300+ private span = 10; // s shown across the width

301+ private drawnX = 0; // pixel columns already drawn

302+ private dirty = true; // full redraw needed

303+ private image: ImageData;

304+

305+ constructor(opts: SpectrogramOptions) {

306+ this.width = opts.width ?? 960;

307+ this.dt = opts.dt;

308+ this.n = opts.n;

309+ this.hop = opts.hop;

310+ this.bins = opts.nfft / 2 + 1;

311+ this.refFreqs = opts.refFreqs;

312+ this.buf = new Float64Array(this.n);

313+ this.hann = new Float64Array(this.n);

314+ for (let i = 0; i < this.n; i++) this.hann[i] = 0.5 - 0.5 * Math.cos((2 * Math.PI * i) / (this.n - 1));

315+ this.fft = new FFT(opts.nfft);

316+ this.re = new Float64Array(opts.nfft);

317+ this.im = new Float64Array(opts.nfft);

318+ this.pending = new Float32Array(this.bins);

319+ this.sum = new Float64Array(this.bins);

320+ this.canvas = document.createElement("canvas");

321+ this.canvas.width = this.width;

322+ this.canvas.height = HEIGHT;

323+ this.image = this.canvas.getContext("2d")!.createImageData(this.width, HEIGHT);

324+ }

325+

326+ // Start a new note shown over `span` seconds.

327+ start(span: number): void {

328+ this.nSamples = 0;

329+ this.cols = [];

330+ this.merge = 1;

331+ this.pending.fill(0);

332+ this.pendingCount = 0;

333+ this.sum.fill(0);

334+ this.peak = 0;

335+ this.span = span;

336+ this.dirty = true;

337+ }

338+

339+ setSpan(span: number): void {

340+ this.span = span;

341+ this.dirty = true;

342+ }

343+

344+ setRefFreqs(f: number[]): void {

345+ this.refFreqs = f;

346+ this.dirty = true;

347+ }

348+

349+ push(x: number): void {

350+ this.buf[this.nSamples % this.n] = x;

351+ this.nSamples++;

352+ if (this.nSamples >= this.n && this.nSamples % this.hop === 0) this.addColumn();

353+ }

354+

355+ // Seconds of note covered by the stored columns.

356+ private covered(): number {

357+ return (this.cols.length * this.merge * this.hop + this.n) * this.dt;

358+ }

359+

360+ private addColumn(): void {

361+ const { buf, re, im, hann, n, bins } = this;

362+ let mean = 0;

363+ for (let i = 0; i < n; i++) mean += buf[i];

364+ mean /= n;

365+ const start = this.nSamples % n;

366+ re.fill(0);

367+ im.fill(0);

368+ for (let i = 0; i < n; i++) re[i] = (buf[(start + i) % n] - mean) * hann[i];

369+ this.fft.transform(re, im);

370+ for (let k = 0; k < bins; k++) {

371+ const p = re[k] * re[k] + im[k] * im[k];

372+ this.pending[k] += p / this.merge;

373+ this.sum[k] += p;

374+ }

375+ if (++this.pendingCount < this.merge) return;

376+ const col = this.pending;

377+ this.pending = new Float32Array(bins);

378+ this.pendingCount = 0;

379+ this.cols.push(col);

380+ let colPeak = 0;

381+ for (let k = 1; k < bins; k++) if (col[k] > colPeak) colPeak = col[k];

382+ if (colPeak > 2 * this.peak) this.dirty = true; // 3 dB louder: rescale colors

383+ this.peak = Math.max(this.peak, colPeak);

384+ if (this.cols.length === MAX_COLS) {

385+ const merged: Float32Array[] = [];

386+ for (let c = 0; c < MAX_COLS; c += 2) {

387+ const a = this.cols[c], b = this.cols[c + 1], m = new Float32Array(bins);

388+ for (let k = 0; k < bins; k++) m[k] = 0.5 * (a[k] + b[k]);

389+ merged.push(m);

390+ }

391+ this.cols = merged;

392+ this.merge *= 2;

393+ this.dirty = true;

394+ }

395+ }

396+

397+ draw(): void {

398+ const { bins, cols } = this;

399+ const df = 1 / (2 * (bins - 1) * this.dt);

400+ let sumPeak = 0;

401+ for (let k = 1; k < bins; k++) if (this.sum[k] > sumPeak) sumPeak = this.sum[k];

402+ let top = bins - 1;

403+ while (top > 1 && this.sum[top] < sumPeak * 10 ** (-CONTENT_DB / 10)) top--;

404+ if (cols.length > 0) {

405+ const fmax = Math.min(niceCeil(1.2 * top * df), (bins - 1) * df);

406+ if (fmax !== this.fmax) this.dirty = true;

407+ this.fmax = fmax;

408+ }

409+ // Pixel columns whose time is covered by stored columns.

410+ const ready = Math.min(this.width, Math.floor((this.covered() / this.span) * this.width));

411+ if (this.dirty) this.drawnX = 0;

412+ if (!this.dirty && ready <= this.drawnX) return;

413+

414+ const data = this.image.data;

415+ const peakDb = 10 * Math.log10(this.peak + 1e-30);

416+ const colDt = this.merge * this.hop * this.dt; // s per stored column

417+ const t0 = this.n * this.dt; // the first column ends here

418+ const from = this.dirty ? 0 : this.drawnX;

419+ const to = this.dirty ? this.width : ready;

420+ for (let x = from; x < to; x++) {

421+ const t = ((x + 0.5) / this.width) * this.span;

422+ const c = Math.round((t - t0) / colDt);

423+ const col = x < ready && c >= 0 && c < cols.length ? cols[c] : null;

424+ for (let y = 0; y < HEIGHT; y++) {

425+ const p = 4 * (y * this.width + x);

426+ if (!col) {

427+ data[p] = data[p + 1] = data[p + 2] = 246; // no data (yet)

428+ } else {

429+ const f = (this.fmax * (HEIGHT - 1 - y)) / (HEIGHT - 1);

430+ const k = Math.min(bins - 1, Math.round(f / df));

431+ const db = 10 * Math.log10(col[k] + 1e-30);

432+ colormap((db - (peakDb - RANGE_DB)) / RANGE_DB, data, p);

433+ }

434+ data[p + 3] = 255;

435+ }

436+ }

437+ this.drawnX = ready;

438+ this.dirty = false;

439+

440+ const ctx = this.canvas.getContext("2d")!;

441+ ctx.putImageData(this.image, 0, 0);

442+ const py = (f: number) => (HEIGHT - 1) * (1 - f / this.fmax);

443+ ctx.strokeStyle = "rgba(255,255,255,0.7)";

444+ ctx.setLineDash([4, 4]);

445+ for (const f of this.refFreqs) {

446+ if (f >= this.fmax) continue;

447+ ctx.beginPath();

448+ ctx.moveTo(0, py(f));

449+ ctx.lineTo(this.width, py(f));

450+ ctx.stroke();

451+ }

452+ ctx.setLineDash([]);

453+ ctx.font = "12px ui-monospace, monospace";

454+ label(ctx, `${this.fmax} Hz`, 2, 2);

455+ label(ctx, "0", 2, HEIGHT - 18);

456+ const end = `${this.span.toFixed(1)} s`;

457+ label(ctx, end, this.width - ctx.measureText(end).width - 8, HEIGHT - 18);

458+ }

459+}

src/traces.tsmodified+27−12View file

@@ -17,11 +17,12 @@ export interface Group {

1717 series: Series[];

1818 shown?: boolean; // initially shown

1919 window?: number; // visible time span for this chart (default: the panel's)

20- log?: boolean; // log10 y axis (non-positive values leave a gap)

20+ log?: boolean; // log10 y axis (non-positive values leave a gap); range in whole decades

21+ zero?: boolean; // non-negative quantity: the y range starts at 0

2122 }

2223

2324 const CAPACITY = 16384; // samples kept per series (ring buffer)

24-const WIDTH = 480, HEIGHT = 92; // CSS pixels; drawn at devicePixelRatio

25+const HEIGHT = 92; // CSS pixels; drawn at devicePixelRatio

2526 const TOP = 18; // legend row height; the plot is below it

2627 const COLORS = ["#2f5d8a", "#c8742c", "#4a8c5c", "#a8476e"];

2728 const FONT = "11px ui-monospace, SFMono-Regular, Menlo, monospace";

@@ -36,18 +37,21 @@ export class TracePanel {

3637 private head = 0; // next write index

3738 private canvases: HTMLCanvasElement[];

3839 private dpr = window.devicePixelRatio || 1;

40+ private width: number; // chart width, CSS pixels

41+ private fixed: { t0: number; span: number } | null = null; // fixed x range instead of scrolling

3942 // Per draw: buffer index and pixel column of each visible sample, oldest first.

4043 private idx = new Int32Array(CAPACITY);

4144 private col = new Int32Array(CAPACITY);

4245

4346 // Toggle chips go into `toggles`, charts into `charts`.

44- constructor(toggles: HTMLElement, charts: HTMLElement, groups: Group[], window: number) {

47+ constructor(toggles: HTMLElement, charts: HTMLElement, groups: Group[], window: number, width = 480) {

48+ this.width = width;

4549 this.groups = groups;

4650 this.window = window;

4751 this.v = groups.map((g) => g.series.map(() => new Float64Array(CAPACITY)));

4852 this.canvases = groups.map((g) => {

4953 const canvas = document.createElement("canvas");

50- canvas.width = WIDTH * this.dpr;

54+ canvas.width = this.width * this.dpr;

5155 canvas.height = HEIGHT * this.dpr;

5256 canvas.hidden = !g.shown;

5357 const box = document.createElement("input");

@@ -69,6 +73,15 @@ export class TracePanel {

6973 return this.canvases.some((c) => !c.hidden);

7074 }

7175

76+ setWindow(window: number): void {

77+ this.window = window;

78+ }

79+

80+ // Show times t0 .. t0 + span across the width (no scrolling); null to scroll.

81+ setFixed(range: { t0: number; span: number } | null): void {

82+ this.fixed = range;

83+ }

84+

7285 clear(): void {

7386 this.count = 0;

7487 this.head = 0;

@@ -87,14 +100,14 @@ export class TracePanel {

87100 // Indices and pixel columns of the samples in the last `win` seconds.

88101 private layout(win: number): number {

89102 const newest = (this.head - 1 + CAPACITY) % CAPACITY;

90- const tStart = this.t[newest] - win;

103+ const tStart = this.fixed ? this.fixed.t0 : this.t[newest] - win;

91104 let n = 0;

92105 while (n < this.count && this.t[(newest - n + CAPACITY) % CAPACITY] >= tStart) n++;

93106 const first = (newest - n + 1 + CAPACITY) % CAPACITY;

94107 for (let m = 0; m < n; m++) {

95108 const i = first + m < CAPACITY ? first + m : first + m - CAPACITY;

96109 this.idx[m] = i;

97- this.col[m] = Math.floor(((this.t[i] - tStart) / win) * WIDTH);

110+ this.col[m] = Math.floor(((this.t[i] - tStart) / win) * this.width);

98111 }

99112 return n;

100113 }

@@ -108,7 +121,7 @@ export class TracePanel {

108121 this.groups.forEach((g, j) => {

109122 const canvas = this.canvases[j];

110123 if (canvas.hidden) return;

111- const win = g.window ?? this.window;

124+ const win = this.fixed ? this.fixed.span : (g.window ?? this.window);

112125 if (win !== laidOut) n = this.layout((laidOut = win));

113126 const tv = g.log ? (x: number) => (x > 0 ? Math.log10(x) : NaN) : (x: number) => x;

114127 const ctx = canvas.getContext("2d")!;

@@ -123,12 +136,14 @@ export class TracePanel {

123136 }

124137 const empty = !(hi >= lo);

125138 if (empty) { lo = 0; hi = 1; }

139+ if (g.zero) lo = 0;

140+ if (g.log) { lo = Math.floor(lo); hi = Math.ceil(hi); }

126141 if (hi - lo <= 1e-300) { lo -= 1; hi += 1; } // constant trace

127- const pad = 0.1 * (hi - lo);

128- const yMin = lo - pad, yMax = hi + pad;

142+ const pad = g.log ? 0 : 0.1 * (hi - lo);

143+ const yMin = g.zero || g.log ? lo : lo - pad, yMax = hi + pad;

129144 const py = (x: number) => HEIGHT - ((x - yMin) / (yMax - yMin)) * (HEIGHT - TOP);

130145

131- ctx.clearRect(0, 0, WIDTH, HEIGHT);

146+ ctx.clearRect(0, 0, this.width, HEIGHT);

132147 ctx.lineWidth = 1.25;

133148 ctx.lineJoin = "round";

134149 g.series.forEach((s, k) => {

@@ -195,10 +210,10 @@ export class TracePanel {

195210 ctx.fillText(g.unit, x, 3);

196211 }

197212 if (!empty) {

198- const show = g.log ? (y: number) => fmt(10 ** y) : fmt;

213+ const show = g.log ? (y: number) => `1e${y}` : (y: number) => (y === 0 ? "0" : fmt(y));

199214 ctx.fillStyle = MUTED;

200215 ctx.textAlign = "right";

201- ctx.fillText(`${show(lo)} … ${show(hi)}${g.log ? " (log)" : ""}`, WIDTH - 4, 3);

216+ ctx.fillText(`${show(lo)} … ${show(hi)}${g.log ? " (log)" : ""}`, this.width - 4, 3);

202217 ctx.textAlign = "left";

203218 }

204219 });