Self-equilibrated release force and loudness normalization - proof-of-concept/steelsim21

proof-of-concept / steelsim21

Self-equilibrated release force and loudness normalization

The recorded finger reaction had a rigid-body component, so applying it as an external force during the release accelerated the whole plate and the Rayleigh sum turned that into a 20 ms pressure pulse (a click). The rigid components are now projected out at release. Playback and WAV export gain a loudness mode (RMS to -16 dBFS with a soft clipper), the new default.

Jeremy Magland <jmagland@flatironinstitute.org> committed 1 month ago commit b9a5a7eeaa3d parent 6f03b39 Browse files

7 changed files+102−16

THEORY.mdmodified+2−2View file

@@ -37,7 +37,7 @@ The same projection defines what "holding" the sheet means during the push. The

3737

3838 ## 3. Excitation

3939

40-The simulation has three phases. During the push (default 0.1 s) the transverse displacement of the nodes within the finger radius follows a smooth ramp to the push amount, and heavy artificial damping ($\sigma_0 = 300\ \mathrm{s^{-1}}$) removes the waves that the ramp launches. During the hold the constraint stays in place until the non-rigid kinetic energy drops below a threshold or a maximum settle time elapses. At the release the clock is set to zero, the damping returns to its physical value, and the recording of audio and movie frames begins. The finger does not let go instantaneously: the reaction force that held each node of the patch is recorded at the end of the hold and then applied as an external force that falls smoothly to zero over the release time (20 ms by default). This matters for the sound. An instantaneous release is a step in the finger force, which excites all modes at once, and because the radiated pressure is proportional to the plate acceleration, the result is a nearly flat pressure spectrum, a cymbal-like click. A release over 20 ms confines the excitation to roughly the first hundred hertz, as a fingertip does. The release time can be set to zero to recover the instantaneous case.

40+The simulation has three phases. During the push (default 0.1 s) the transverse displacement of the nodes within the finger radius follows a smooth ramp to the push amount, and heavy artificial damping ($\sigma_0 = 300\ \mathrm{s^{-1}}$) removes the waves that the ramp launches. During the hold the constraint stays in place until the non-rigid kinetic energy drops below a threshold or a maximum settle time elapses. At the release the clock is set to zero, the damping returns to its physical value, and the recording of audio and movie frames begins. The finger does not let go instantaneously: the reaction force that held each node of the patch is recorded at the end of the hold and then applied as an external force that falls smoothly to zero over the release time (20 ms by default). This matters for the sound. An instantaneous release is a step in the finger force, which excites all modes at once, and because the radiated pressure is proportional to the plate acceleration, the result is a nearly flat pressure spectrum, a cymbal-like click. A release over 20 ms confines the excitation to roughly the first hundred hertz, as a fingertip does. The rigid-body components of the recorded reaction are removed before it is applied, so that the finger force and the distributed holding reaction relax together and the plate is not accelerated as a whole. The release time can be set to zero to recover the instantaneous case.

4141

4242 ## 4. Spatial discretization

4343

@@ -69,7 +69,7 @@ $$

6969

7070 where the sum runs over the nodes, $r_k$ is the distance from node $k$ to the listener, and $c = 343$ m/s. The transverse acceleration of every node is recorded in a ring buffer at the audio rate, and each output sample is a delayed, weighted sum over the plate, evaluated on the GPU with linear interpolation between samples. The common delay is discarded; only the spread of delays across the plate (about 5 samples for a listener 1 m in front) is kept. The plate is unbaffled in reality, so the sound from the two faces partly cancels at low frequencies. The "dipole" option is a heuristic for this: each element is weighted by the cosine of its angle to the listener, and the signal is high-pass filtered with a first-order filter at $c/(2\sqrt{L_x L_y})$, about 460 Hz for the default sheet. An accurate treatment would require solving the exterior acoustic problem, which is outside the scope of this simulator.

7171

72-The pressure is reported in pascals, and the peak sound pressure level is displayed. Playback and the downloaded WAV file are normalized to their peak.

72+The pressure is reported in pascals, and the peak sound pressure level is displayed. Playback and the downloaded WAV file are normalized either to their peak or, by default, to the loudness of the sustained sound (the RMS after the first 50 ms is scaled to -16 dBFS and a soft clipper bounds the onset), since a signal dominated by the lowest modes has a large peak-to-RMS ratio.

7373

7474 ## 7. Limitations

7575

src/audio/level.tsadded+25−0View file

@@ -0,0 +1,25 @@

1+// Playback level: peak normalization, or loudness normalization with a soft clipper.

2+

3+export type LevelMode = "peak" | "loudness";

4+

5+/** Returns a copy of `samples` scaled for playback or export, in [-1, 1]. */

6+export function normalize(samples: Float32Array, fs: number, mode: LevelMode, peak = 0.95): Float32Array {

7+ const out = new Float32Array(samples.length);

8+ if (samples.length === 0) return out;

9+ if (mode === "peak") {

10+ let max = 0;

11+ for (let k = 0; k < samples.length; k++) max = Math.max(max, Math.abs(samples[k]));

12+ const g = max > 0 ? peak / max : 0;

13+ for (let k = 0; k < samples.length; k++) out[k] = samples[k] * g;

14+ return out;

15+ }

16+ // Loudness: scale the RMS of the sustained part (after the first 50 ms) to -16 dBFS,

17+ // then soft-clip so that transients cannot exceed the peak.

18+ const start = Math.min(samples.length - 1, Math.round(0.05 * fs));

19+ let sum = 0;

20+ for (let k = start; k < samples.length; k++) sum += samples[k] * samples[k];

21+ const rms = Math.sqrt(sum / Math.max(1, samples.length - start));

22+ const g = rms > 0 ? 10 ** (-16 / 20) / rms : 0;

23+ for (let k = 0; k < samples.length; k++) out[k] = peak * Math.tanh((samples[k] * g) / peak);

24+ return out;

25+}

src/audio/player.tsmodified+4−6View file

@@ -1,5 +1,7 @@

11 // Web Audio playback of a Float32Array signal with position reporting.

22

3+import { normalize, type LevelMode } from "./level.ts";

4+

35 export class Player {

46 private ctx: AudioContext | null = null;

57 private source: AudioBufferSourceNode | null = null;

@@ -18,17 +20,13 @@ export class Player {

1820 return Math.min(this.length, this.offset + this.ctx.currentTime - this.startedAt);

1921 }

2022

21- play(samples: Float32Array, fs: number, offset = 0, peak = 0.9): void {

23+ play(samples: Float32Array, fs: number, offset = 0, mode: LevelMode = "peak"): void {

2224 this.stop();

2325 if (samples.length === 0) return;

2426 if (!this.ctx) this.ctx = new AudioContext({ sampleRate: fs });

2527 const ctx = this.ctx;

2628 const buffer = ctx.createBuffer(1, samples.length, fs);

27- const ch = buffer.getChannelData(0);

28- let max = 0;

29- for (let k = 0; k < samples.length; k++) max = Math.max(max, Math.abs(samples[k]));

30- const gain = max > 0 ? peak / max : 0;

31- for (let k = 0; k < samples.length; k++) ch[k] = samples[k] * gain;

29+ buffer.getChannelData(0).set(normalize(samples, fs, mode, 0.9));

3230 const src = ctx.createBufferSource();

3331 src.buffer = buffer;

3432 src.connect(ctx.destination);

src/audio/wav.tsmodified+5−5View file

@@ -1,9 +1,9 @@

11 // 16-bit PCM mono WAV encoder.

22

3-export function encodeWav(samples: Float32Array, fs: number, peak = 0.99): Blob {

4- let max = 0;

5- for (let k = 0; k < samples.length; k++) max = Math.max(max, Math.abs(samples[k]));

6- const gain = max > 0 ? peak / max : 0;

3+import { normalize, type LevelMode } from "./level.ts";

4+

5+export function encodeWav(samples: Float32Array, fs: number, mode: LevelMode = "peak"): Blob {

6+ const scaled = normalize(samples, fs, mode, 0.99);

77 const buf = new ArrayBuffer(44 + samples.length * 2);

88 const v = new DataView(buf);

99 const str = (o: number, s: string) => { for (let k = 0; k < s.length; k++) v.setUint8(o + k, s.charCodeAt(k)); };

@@ -21,7 +21,7 @@ export function encodeWav(samples: Float32Array, fs: number, peak = 0.99): Blob

2121 str(36, "data");

2222 v.setUint32(40, samples.length * 2, true);

2323 for (let k = 0; k < samples.length; k++) {

24- const s = Math.max(-1, Math.min(1, samples[k] * gain));

24+ const s = Math.max(-1, Math.min(1, scaled[k]));

2525 v.setInt16(44 + 2 * k, Math.round(s * 32767), true);

2626 }

2727 return new Blob([buf], { type: "audio/wav" });

src/main.tsmodified+10−3View file

@@ -6,6 +6,7 @@ import { Renderer, type ViewMode } from "./render/renderer.ts";

66 import { Spectrogram } from "./audio/spectrogram.ts";

77 import { Player } from "./audio/player.ts";

88 import { encodeWav } from "./audio/wav.ts";

9+import { type LevelMode } from "./audio/level.ts";

910

1011 const $ = <T extends HTMLElement>(id: string) => document.getElementById(id) as T;

1112

@@ -251,6 +252,7 @@ const timeLabel = document.createElement("span");

251252 const btnLive = mkButton("Live");

252253 const btnPlay = mkButton("Play");

253254 const btnWav = mkButton("Download .wav");

255+const levelSel = document.createElement("select");

254256 let exagInput: HTMLInputElement;

255257 let rangeInput: HTMLInputElement;

256258

@@ -288,21 +290,26 @@ function buildMovieToolbar(): void {

288290 if (!sim || sim.samples === 0) return;

289291 const offset = live ? 0 : scrubFrame / sim.fps;

290292 live = false;

291- player.play(sim.processedAudio(), sim.fs, offset);

293+ player.play(sim.processedAudio(), sim.fs, offset, levelSel.value as LevelMode);

292294 btnPlay.textContent = "Stop";

293295 updateLiveButton();

294296 });

295297 player.onEnded = () => { btnPlay.textContent = "Play"; };

296298 btnWav.addEventListener("click", () => {

297299 if (!sim || sim.samples === 0) return;

298- const blob = encodeWav(sim.processedAudio(), sim.fs);

300+ const blob = encodeWav(sim.processedAudio(), sim.fs, levelSel.value as LevelMode);

299301 const a = document.createElement("a");

300302 a.href = URL.createObjectURL(blob);

301303 a.download = `sheet_${sim.samples}samples_${sim.fs}Hz.wav`;

302304 a.click();

303305 setTimeout(() => URL.revokeObjectURL(a.href), 5000);

304306 });

305- tb.append(btnPlay, btnLive, scrub, timeLabel, btnWav);

307+ for (const [v, t] of [["loudness", "Loudness (RMS, soft clip)"], ["peak", "Peak"]]) {

308+ const o = document.createElement("option");

309+ o.value = v; o.textContent = t;

310+ levelSel.appendChild(o);

311+ }

312+ tb.append(btnPlay, btnLive, scrub, timeLabel, labeled("Level", levelSel), btnWav);

306313 updateLiveButton();

307314 }

308315

src/sim/shaders/rigid.wgslmodified+28−0View file

@@ -66,3 +66,31 @@ fn apply(@builtin(global_invocation_id) gid: vec3u) {

6666 vel[n + k] -= stats[7] + stats[9] * x;

6767 vel[2u * n + k] -= stats[8] + stats[10] * y - stats[11] * x;

6868 }

69+

70+// Removes the three rigid components (translation and two rotations) from a

71+// per-node acceleration field so that it is self-equilibrated. One workgroup.

72+@group(0) @binding(4) var<storage, read_write> field: array<f32>;

73+var<workgroup> sh3: array<vec3f, 256>;

74+

75+@compute @workgroup_size(256)

76+fn projectField(@builtin(local_invocation_id) lid: vec3u) {

77+ let t = lid.x;

78+ var s = vec3f(0.0);

79+ for (var k = t; k < P.n; k += 256u) {

80+ let i = i32(k % P.nx); let j = i32(k / P.nx);

81+ let m = P.rhoh * P.hx * P.hy * trapWeight(i, j);

82+ let f = field[k];

83+ s += vec3f(m * f, m * nodeY(j) * f, -m * nodeX(i) * f);

84+ }

85+ sh3[t] = s;

86+ workgroupBarrier();

87+ for (var stride = 128u; stride > 0u; stride >>= 1u) {

88+ if (t < stride) { sh3[t] += sh3[t + stride]; }

89+ workgroupBarrier();

90+ }

91+ let c = vec3f(sh3[0].x / P.mass, sh3[0].y / P.Ix, sh3[0].z / P.Iy);

92+ for (var k = t; k < P.n; k += 256u) {

93+ let i = i32(k % P.nx); let j = i32(k / P.nx);

94+ field[k] -= c.x + c.y * nodeY(j) - c.z * nodeX(i);

95+ }

96+}

src/sim/simulation.tsmodified+28−0View file

@@ -96,6 +96,8 @@ export class Simulation {

9696 private stepPipe: GPUComputePipeline;

9797 private sumsPipe: GPUComputePipeline;

9898 private applyPipe: GPUComputePipeline;

99+ private projectPipe: GPUComputePipeline;

100+ private projectBG: GPUBindGroup;

99101 private audioPipe: GPUComputePipeline;

100102 private extractPipe: GPUComputePipeline;

101103 private stepBG: GPUBindGroup[];

@@ -202,6 +204,16 @@ export class Simulation {

202204 const rbPipeLayout = device.createPipelineLayout({ bindGroupLayouts: [rbLayout] });

203205 this.sumsPipe = device.createComputePipeline({ layout: rbPipeLayout, compute: { module: rigidModule, entryPoint: "sums" } });

204206 this.applyPipe = device.createComputePipeline({ layout: rbPipeLayout, compute: { module: rigidModule, entryPoint: "apply" } });

207+ const projectLayout = device.createBindGroupLayout({

208+ entries: [

209+ { binding: 0, visibility: GPUShaderStage.COMPUTE, buffer: { type: "uniform" } },

210+ { binding: 4, visibility: GPUShaderStage.COMPUTE, buffer: { type: "storage" } },

211+ ],

212+ });

213+ this.projectPipe = device.createComputePipeline({

214+ layout: device.createPipelineLayout({ bindGroupLayouts: [projectLayout] }),

215+ compute: { module: rigidModule, entryPoint: "projectField" },

216+ });

205217 this.audioPipe = mk(audioWgsl, "main", device.createPipelineLayout({ bindGroupLayouts: [audioLayout] }));

206218 this.extractPipe = mk(extractWgsl, "main");

207219

@@ -232,6 +244,13 @@ export class Simulation {

232244 ],

233245 }),

234246 );

247+ this.projectBG = device.createBindGroup({

248+ layout: projectLayout,

249+ entries: [

250+ { binding: 0, resource: { buffer: this.paramsBuf } },

251+ { binding: 4, resource: { buffer: this.react } },

252+ ],

253+ });

235254 this.audioBG = device.createBindGroup({

236255 layout: audioLayout,

237256 entries: [

@@ -515,6 +534,15 @@ export class Simulation {

515534 this.phase = "recording";

516535 this.stepCount = 0;

517536 this.totalSteps = 0;

537+ // The finger force and the distributed holding reaction relax together:

538+ // remove the rigid-body components of the recorded reaction.

539+ const enc = this.device.createCommandEncoder();

540+ const pass = enc.beginComputePass();

541+ pass.setPipeline(this.projectPipe);

542+ pass.setBindGroup(0, this.projectBG);

543+ pass.dispatchWorkgroups(1);

544+ pass.end();

545+ this.device.queue.submit([enc.finish()]);

518546 // Frame 0 is the held, settled shape.

519547 if (this.lastPrepFrame) this.frames.push(this.lastPrepFrame);

520548 }