High-frequency ringing: in-plane loss, band-limited listener sum, damping defaults - proof-of-concept/steelsim21

proof-of-concept / steelsim21

High-frequency ringing: in-plane loss, band-limited listener sum, damping defaults

Two sources of high-frequency content that outlasted the low modes: the stretching modes of the curved sheet (kHz, coupled to w by the bow) had only the sigma0 loss, and the listener sum sampled the instantaneous acceleration, which carries the per-step float32 rounding noise of the force evaluation. The ring buffer now stores the transverse velocity and the acceleration is the first difference of successive samples. Adds an in-plane loss parameter (default 20/s) and raises sigma1 to 5e-4 to match the estimated radiation damping.

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

7 changed files+32−18

THEORY.mdmodified+2−2View file

@@ -25,7 +25,7 @@ $$

2525

2626 with the bending stiffness $D = Eh^3/(12(1-\nu^2))$. The quadratic terms in the strains couple bending to stretching. They are what produce the gong-like behaviour: the pitch glide of a strongly struck plate, and the cascade of energy from the low modes that the push excites toward high frequencies. The in-plane inertia terms $\rho h u_{tt}$ and $\rho h v_{tt}$ are usually dropped in favour of an Airy stress function, on the grounds that in-plane waves are much faster than bending waves. We keep them, as requested, at the cost of a much smaller time step (Section 5).

2727

28-The loss terms follow the two-parameter model that is standard in plate sound synthesis. A mode of angular frequency $\omega$ decays at the rate $\sigma_0 + \sigma_1 k^2$ with $k^2 = \omega/\kappa$ and $\kappa = \sqrt{D/\rho h}$, so $\sigma_0$ sets the low-frequency decay and $\sigma_1$ makes high frequencies decay faster. The defaults ($\sigma_0 = 1.15\ \mathrm{s^{-1}}$, $\sigma_1 = 1.08\times10^{-4}\ \mathrm{m^2/s}$) give a 60 dB decay time of about 6 s at low frequencies and about 1 s near 10 kHz. These are not measured values; both parameters are exposed in the interface.

28+The loss terms follow the two-parameter model that is standard in plate sound synthesis. A bending mode of angular frequency $\omega$ decays at the rate $\sigma_0 + \sigma_1 k^2$ with $k^2 = \omega/\kappa$ and $\kappa = \sqrt{D/\rho h}$, so $\sigma_0$ sets the low-frequency decay and $\sigma_1$ makes high frequencies decay faster. The defaults ($\sigma_0 = 1.15\ \mathrm{s^{-1}}$, $\sigma_1 = 5\times10^{-4}\ \mathrm{m^2/s}$) give a 60 dB decay time of about 5 s at 100 Hz, 1.8 s at 1 kHz, and 0.5 s at 5 kHz. The value of $\sigma_1$ is chosen to match an estimate of the radiation damping of a thin steel sheet in air, which for this thickness corresponds to a loss factor of roughly $10^{-3}$ at 1 kHz; radiation, not the metal itself, is the main loss of a thin sheet at mid and high frequencies. The in-plane motion carries a separate frequency-independent loss $\sigma_{\mathrm{in}}$ (default $20\ \mathrm{s^{-1}}$). It matters because the stretching modes of the sheet lie at several kilohertz, couple linearly to transverse motion once the sheet is curved (Section 1, initial shape), and would otherwise ring for the whole low-frequency decay time; a material loss factor of about $10^{-3}$ at 5 kHz gives a rate of this order. None of these are measured values, and all three are exposed in the interface.

2929

3030 ### Initial shape (imperfection)

3131

@@ -79,7 +79,7 @@ $$

7979 p(\mathbf r, t) = \frac{\rho_{\mathrm{air}}}{2\pi} \sum_k \frac{\Delta A_k}{r_k}\, \ddot w_k(t - r_k/c),

8080 $$

8181

82-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.

82+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 velocity 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 acceleration is formed as the difference between successive delayed velocities, which is the average acceleration over one sample interval. This matters in single precision: the force evaluation carries a broadband rounding error of the order of $1\ \mathrm{m/s^2}$ per node, and sampling the instantaneous acceleration would put that noise straight into the audio, whereas the averaged form suppresses it by the square root of the number of steps per sample. 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.

8383

8484 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.

8585

src/main.tsmodified+4−3View file

@@ -58,6 +58,7 @@ const groups: { title: string; open: boolean; fields: FieldSpec[] }[] = [

5858 fields: [

5959 { key: "sigma0", label: "σ₀ (freq. independent)", unit: "1/s", scale: 1, step: 0.1, live: true },

6060 { key: "sigma1", label: "σ₁ (freq. dependent)", unit: "m²/s", scale: 1, step: 1e-5, live: true },

61+ { key: "sigmaIn", label: "In-plane loss", unit: "1/s", scale: 1, step: 5, live: true },

6162 ],

6263 },

6364 {

@@ -149,7 +150,7 @@ function buildForm(): void {

149150 note.textContent = "Real sheet metal is never flat. A bow of a few mm makes the sheet a shallow shell with an asymmetric stiffness; push with a negative amount to push against the bow.";

150151 det.appendChild(note);

151152 }

152- if (f.key === "sigma1") {

153+ if (f.key === "sigmaIn") {

153154 const note = document.createElement("p");

154155 note.className = "sub";

155156 note.id = "dampNote";

@@ -172,9 +173,9 @@ function onFieldChange(f: FieldSpec, inp: HTMLInputElement): void {

172173 const v = parseFloat(inp.value);

173174 if (!isFinite(v)) return;

174175 (params as unknown as Record<string, number>)[f.key] = v * f.scale;

175- if (f.key === "sigma0" || f.key === "sigma1") {

176+ if (f.key === "sigma0" || f.key === "sigma1" || f.key === "sigmaIn") {

176177 updateDampingNote();

177- sim?.setDamping(params.sigma0, params.sigma1);

178+ sim?.setDamping(params.sigma0, params.sigma1, params.sigmaIn);

178179 }

179180 if (f.key.startsWith("listener")) applyListener();

180181 }

src/params.tsmodified+3−1View file

@@ -9,6 +9,7 @@ export interface PhysParams {

99 rho: number; // density [kg/m^3]

1010 sigma0: number; // frequency-independent loss [1/s]

1111 sigma1: number; // frequency-dependent loss [m^2/s]

12+ sigmaIn: number; // additional loss of the in-plane motion [1/s]

1213 holeInset: number; // distance of the thread hole from the two edges of the top-left corner [m]

1314 pushAmount: number; // corner displacement [m]

1415 fingerRadius: number; // radius of the pushed patch [m]

@@ -39,7 +40,8 @@ export const defaultParams: PhysParams = {

3940 nu: 0.29,

4041 rho: 7900,

4142 sigma0: 1.15,

42- sigma1: 1.08e-4,

43+ sigma1: 5e-4,

44+ sigmaIn: 20,

4345 holeInset: 0.012,

4446 pushAmount: 0.01,

4547 fingerRadius: 0.01,

src/sim/shaders/audio.wgslmodified+10−3View file

@@ -1,6 +1,10 @@

11 // Retarded-time sum over the plate for one output sample (Rayleigh integral).

22 // lw[k] = (weight, delay in samples) for node k; delay is relative to the

33 // closest node so that the ring buffer only has to cover the spread.

4+// The ring buffer holds the transverse velocity at the sample instants; the

5+// acceleration is the first difference of the delayed velocity, which is the

6+// average acceleration over one sample interval (band-limited, and immune to

7+// the per-step rounding noise of the force evaluation).

48

59 @group(0) @binding(0) var<uniform> P: Params;

610 @group(0) @binding(1) var<uniform> SP: StepParams;

@@ -22,9 +26,12 @@ fn main(@builtin(local_invocation_id) lid: vec3u) {

2226 let f = d - f32(k0);

2327 let s0 = ((slot - k0) % H + H) % H;

2428 let s1 = ((slot - k0 - 1) % H + H) % H;

25- let a0 = hist[u32(s0) * P.n + k];

26- let a1 = hist[u32(s1) * P.n + k];

27- acc += lw[k].x * mix(a0, a1, f);

29+ let s2 = ((slot - k0 - 2) % H + H) % H;

30+ let v0 = hist[u32(s0) * P.n + k];

31+ let v1 = hist[u32(s1) * P.n + k];

32+ let v2 = hist[u32(s2) * P.n + k];

33+ let a = (mix(v0, v1, f) - mix(v1, v2, f)) * P.fs;

34+ acc += lw[k].x * a;

2835 }

2936 sh[t] = acc;

3037 workgroupBarrier();

src/sim/shaders/common.wgslmodified+2−2View file

@@ -20,7 +20,7 @@ struct Params {

2020 Ix: f32, // sum m y^2

2121 Iy: f32, // sum m x^2

2222 Iz: f32, // sum m (x^2 + y^2)

23- pad0: f32,

23+ fs: f32, // audio sample rate [Hz]

2424 };

2525

2626 // Per-step parameters, addressed with a dynamic uniform offset.

@@ -34,7 +34,7 @@ struct StepParams {

3434 sigma0: f32,

3535 sigma1: f32,

3636 releaseGain: f32, // fraction of the recorded reaction still applied (finger letting go)

37- pad0: f32,

37+ sigmaIn: f32, // additional frequency-independent loss of the in-plane motion

3838 pad1: f32,

3939 pad2: f32,

4040 };

src/sim/shaders/step.wgslmodified+5−3View file

@@ -12,7 +12,7 @@

1212 @group(0) @binding(3) var<storage, read> velIn: array<f32>;

1313 @group(0) @binding(4) var<storage, read_write> posOut: array<f32>;

1414 @group(0) @binding(5) var<storage, read_write> velOut: array<f32>;

15-@group(0) @binding(6) var<storage, read_write> hist: array<f32>; // [slot * n + node] acceleration of w

15+@group(0) @binding(6) var<storage, read_write> hist: array<f32>; // [slot * n + node] velocity of w at the audio sample instants

1616 @group(0) @binding(7) var<storage, read_write> react: array<f32>; // constraint reaction (acceleration of w) per node

1717 @group(0) @binding(8) var<storage, read> w0: array<f32>; // stress-free initial shape

1818

@@ -160,13 +160,15 @@ fn main(@builtin(global_invocation_id) gid: vec3u, @builtin(local_invocation_id)

160160 if (j < i32(P.ny) - 1) { lap += (DW(i, j + 1) - vel.z) * ihy2; }

161161 acc -= 2.0 * SP.sigma0 * vel;

162162 acc.z += 2.0 * SP.sigma1 * lap;

163+ // In-plane (stretching) waves lie at several kHz; they get their own loss.

164+ acc.x -= 2.0 * SP.sigmaIn * vel.x;

165+ acc.y -= 2.0 * SP.sigmaIn * vel.y;

163166

164167 // Finger letting go: the reaction recorded at the end of the hold, ramped to zero.

165168 if ((SP.flags & 4u) != 0u) { acc.z += react[k] * SP.releaseGain; }

166169

167- if (SP.audioSlot >= 0) { hist[u32(SP.audioSlot) * n + k] = acc.z; }

168-

169170 var vNew = vel + P.dt * acc;

171+ if (SP.audioSlot >= 0) { hist[u32(SP.audioSlot) * n + k] = vNew.z; }

170172 let pOld = vec3f(posIn[k], posIn[n + k], posIn[2u * n + k]);

171173 var pNew = pOld + P.dt * vNew;

172174 // Compensated update of w: the increment dt * v can be far below the

src/sim/simulation.tsmodified+6−4View file

@@ -162,7 +162,7 @@ export class Simulation {

162162 pf[4] = g.hx; pf[5] = g.hy; pf[6] = this.dt; pf[7] = m.nu;

163163 pf[8] = m.D; pf[9] = m.C; pf[10] = m.G; pf[11] = m.rhoh;

164164 pf[12] = params.Lx / 2; pf[13] = -params.Ly / 2; pf[14] = params.fingerRadius ** 2; pf[15] = mass;

165- pf[16] = Ix; pf[17] = Iy; pf[18] = Ix + Iy; pf[19] = 0;

165+ pf[16] = Ix; pf[17] = Iy; pf[18] = Ix + Iy; pf[19] = params.fs;

166166 this.paramsBuf = device.createBuffer({ size: pf.byteLength, usage: GPUBufferUsage.UNIFORM | GPUBufferUsage.COPY_DST });

167167 device.queue.writeBuffer(this.paramsBuf, 0, pf);

168168

@@ -307,8 +307,8 @@ export class Simulation {

307307 this.delayClamped = false;

308308 for (let k = 0; k < g.n; k++) {

309309 let d = data[2 * k + 1] - dmin;

310- if (d > HIST_LEN - 2) {

311- d = HIST_LEN - 2;

310+ if (d > HIST_LEN - 3) {

311+ d = HIST_LEN - 3;

312312 this.delayClamped = true;

313313 }

314314 data[2 * k + 1] = d;

@@ -320,9 +320,10 @@ export class Simulation {

320320 this.params.radiation = radiation;

321321 }

322322

323- setDamping(sigma0: number, sigma1: number): void {

323+ setDamping(sigma0: number, sigma1: number, sigmaIn: number): void {

324324 this.params.sigma0 = sigma0;

325325 this.params.sigma1 = sigma1;

326+ this.params.sigmaIn = sigmaIn;

326327 }

327328

328329 setDuration(duration: number): void {

@@ -452,6 +453,7 @@ export class Simulation {

452453 f32[o + 6] = sigma0;

453454 f32[o + 7] = p.sigma1;

454455 f32[o + 8] = releaseGain;

456+ f32[o + 9] = p.sigmaIn;

455457 }

456458 this.device.queue.writeBuffer(this.stepBuf, 0, this.stepData, 0, steps * this.stepStride);

457459