Add an optional k^4 loss term (sigma2), off by default
A recording of the real sheet showed its high bands decaying faster than
the linear-in-frequency law allows. The new term applies the graph
Laplacian of the transverse velocity twice, is dissipative, and can be
changed while a simulation runs like the other loss parameters.
8 changed files+62−14
THEORY.mdmodified+2−2View file
@@ -20,12 +20,12 @@ $$
2020 $$
2121
2222 $$
23-\rho h\, w_{tt} = -D\nabla^4 w + \partial_x(N_{xx} w_x + N_{xy} w_y) + \partial_y(N_{xy} w_x + N_{yy} w_y) - 2\sigma_0\rho h\, w_t + 2\sigma_1 \rho h\, \nabla^2 w_t ,
23+\rho h\, w_{tt} = -D\nabla^4 w + \partial_x(N_{xx} w_x + N_{xy} w_y) + \partial_y(N_{xy} w_x + N_{yy} w_y) - 2\sigma_0\rho h\, w_t + 2\sigma_1 \rho h\, \nabla^2 w_t - 2\sigma_2 \rho h\, \nabla^4 w_t ,
2424 $$
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 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.
28+The loss terms extend the two-parameter model that is standard in plate sound synthesis by an optional third term. A bending mode of angular frequency $\omega$ decays at the rate $\sigma_0 + \sigma_1 k^2 + \sigma_2 k^4$ with $k^2 = \omega/\kappa$ and $\kappa = \sqrt{D/\rho h}$, so $\sigma_0$ sets the low-frequency decay, $\sigma_1$ makes high frequencies decay faster in proportion to frequency, and $\sigma_2$ (zero by default) lets the loss curve upward above a few kilohertz, as a recording of a real sheet suggested it does; it is implemented as the biharmonic of the velocity, the graph Laplacian applied twice, and is dissipative for the same reason as the $\sigma_1$ term. 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
matching/README.mdmodified+1−0View file
@@ -9,6 +9,7 @@ This directory holds scripts for imposing the frequency-dependent decay of a rec
99 - `plot.py OUT`: spectrograms of the four clips and band decay curves of the recording against the simulation.
1010 - `eq.py REFERENCE.wav SIMULATED.f32 OUT`: splits the envelope-mode gain table into a static per-band equalizer (the balance at 0.1 s) and a time-varying remainder, reports both, and writes the simulation through the static equalizer (`OUT_eq.wav`, 28 band gains) and through a cubic-in-log-frequency approximation of it (`OUT_eq4.wav`, 4 coefficients).
1111 - `partials.py CLIP ...`: the dominant partials of several clips side by side.
12+- `hfcheck.py OUT ...`: from a `match.py` analysis, the high-band levels of recording and simulation at a few times after the onset, and the high-band excess of the simulation at the attack.
1213 - `modes.ts OUT.bin [key=value ...]` and `modes.py OUT.bin OUT.json [fmax]`: linearized mode frequencies of the model sheet (all three displacement fields, including the stiffening by the initial bow), for comparing the model's partials with a recording without running the simulator. Any parameter can be overridden, for example `bowAmplitude=0` or `E=178e9`.
1314
1415 `match.py` takes `--bands N` to set the number of bands per octave (default 3). Twelfth-octave bands are narrow enough to weight individual partials, which makes the result closer to the recording but says less about the loss model.
matching/hfcheck.pyadded+26−0View file
@@ -0,0 +1,26 @@
1+"""High-frequency comparison from a match.py analysis file.
2+
3+Usage: python3 hfcheck.py OUT_PREFIX ...
4+For each prefix, prints the level of the recording and the simulation in the bands
5+from 2.3 to 5.7 kHz at a few times after the onset, relative to 0.1 s, and the static
6+balance (recording minus simulation at 0.1 s) of the 3.6 to 9 kHz bands relative to
7+the 356 to 2263 Hz mids.
8+"""
9+import sys, numpy as np
10+
11+for out in sys.argv[1:]:
12+ a = np.load(f"{out}_analysis.npz")
13+ t, Er, Es, C, floor = a["t"], a["Er"], a["Es"], a["centers"], a["floor"]
14+ i0 = np.searchsorted(t, 0.1)
15+ times = [0.3, 0.6, 1.0, 1.5]
16+ idx = [np.searchsorted(t, x) for x in times]
17+ G0 = (Er - Es)[:, i0]
18+ mid = G0[(C >= 356) & (C <= 2263)].mean()
19+ hf = G0[(C >= 3592) & (C <= 9051)].mean()
20+ print(f"== {out}: high band excess of the simulation at the attack {mid - hf:+.1f} dB (0 = balanced)")
21+ print(" band recording drop (dB) at 0.3/0.6/1.0/1.5 s simulation drop recording floor")
22+ for c in (2263, 2851, 3592, 4525, 5702):
23+ b = int(np.argmin(np.abs(C - c)))
24+ r = " ".join(f"{Er[b, i] - Er[b, i0]:6.0f}" for i in idx)
25+ s = " ".join(f"{Es[b, i] - Es[b, i0]:6.0f}" for i in idx)
26+ print(f" {C[b]:5.0f} Hz {r} {s} {floor[b] - Er[b, i0]:6.0f}")
src/main.tsmodified+3−2View file
@@ -71,6 +71,7 @@ const groups: { title: string; open: boolean; fields: FieldSpec[] }[] = [
7171 fields: [
7272 { key: "sigma0", label: "σ₀ (freq. independent)", unit: "1/s", scale: 1, step: 0.1, live: true },
7373 { key: "sigma1", label: "σ₁ (freq. dependent)", unit: "m²/s", scale: 1, step: 1e-5, live: true },
74+ { key: "sigma2", label: "σ₂ (k⁴, 0 = off)", unit: "×10⁻⁹ m⁴/s", scale: 1e-9, step: 1, live: true },
7475 { key: "sigmaIn", label: "In-plane loss", unit: "1/s", scale: 1, step: 5, live: true },
7576 ],
7677 },
@@ -233,9 +234,9 @@ function onFieldChange(f: FieldSpec, inp: HTMLInputElement): void {
233234 if (!isFinite(v)) return;
234235 (params as unknown as Record<string, number>)[f.key] = v * f.scale;
235236 if (f.key === "Lx" || f.key === "Ly" || f.key === "spacing") updateGridNote();
236- if (f.key === "sigma0" || f.key === "sigma1" || f.key === "sigmaIn") {
237+ if (f.key === "sigma0" || f.key === "sigma1" || f.key === "sigma2" || f.key === "sigmaIn") {
237238 updateDampingNote();
238- sim?.setDamping(params.sigma0, params.sigma1, params.sigmaIn);
239+ sim?.setDamping(params.sigma0, params.sigma1, params.sigmaIn, params.sigma2);
239240 }
240241 if (f.key.startsWith("listener")) applyListener();
241242 }
src/params.tsmodified+4−2View file
@@ -10,6 +10,7 @@ export interface PhysParams {
1010 rho: number; // density [kg/m^3]
1111 sigma0: number; // frequency-independent loss [1/s]
1212 sigma1: number; // frequency-dependent loss [m^2/s]
13+ sigma2: number; // loss rising with the fourth power of the wavenumber [m^4/s]; 0 disables it
1314 sigmaIn: number; // additional loss of the in-plane motion [1/s]
1415 holeInset: number; // distance of the thread hole from the two edges of the top-left corner [m]
1516 pushAmount: number; // corner displacement [m]
@@ -56,6 +57,7 @@ export const defaultParams: PhysParams = {
5657 rho: 7900,
5758 sigma0: 1.15,
5859 sigma1: 5e-4,
60+ sigma2: 0,
5961 sigmaIn: 20,
6062 holeInset: 0.012,
6163 pushAmount: 0.01,
@@ -257,10 +259,10 @@ export function nodeWeights(g: Grid): Float32Array<ArrayBuffer> {
257259 return w;
258260 }
259261
260-/** Decay time T60 [s] of a mode with angular frequency omega under the two-parameter loss model. */
262+/** Decay time T60 [s] of a bending mode with angular frequency omega under the loss model sigma0 + sigma1 k^2 + sigma2 k^4. */
261263 export function t60(m: Material, p: PhysParams, omega: number): number {
262264 const kappa = Math.sqrt(m.D / m.rhoh);
263265 const k2 = omega / kappa;
264- const sigma = p.sigma0 + p.sigma1 * k2;
266+ const sigma = p.sigma0 + p.sigma1 * k2 + (p.sigma2 ?? 0) * k2 * k2;
265267 return (3 * Math.LN10) / sigma; // amplitude factor 1000 = 60 dB
266268 }
src/sim/shaders/common.wgslmodified+1−1View file
@@ -44,7 +44,7 @@ struct StepParams {
4444 releaseGain: f32, // fraction of the recorded reaction still applied (finger letting go)
4545 sigmaIn: f32, // additional frequency-independent loss of the in-plane motion
4646 parity: u32, // which of the two stick-state slots holds the current state
47- pad2: f32,
47+ sigma2: f32, // loss rising with k^4 (biharmonic of the velocity); 0 disables it
4848 };
4949
5050 fn nodeX(i: i32) -> f32 { return (f32(i) - 0.5 * f32(P.nx - 1u)) * P.hx; }
src/sim/shaders/step.wgslmodified+22−6View file
@@ -53,6 +53,17 @@ fn V(i: i32, j: i32) -> f32 { return posIn[P.n + id(i, j)]; }
5353 fn DW(i: i32, j: i32) -> f32 { return velIn[2u * P.n + id(i, j)]; }
5454 fn W0(i: i32, j: i32) -> f32 { return w0[id(i, j)]; }
5555
56+// Graph Laplacian of the transverse velocity at node (i, j), over neighbours in the sheet.
57+fn lapDW(i: i32, j: i32) -> f32 {
58+ let c = DW(i, j);
59+ var lap = 0.0;
60+ if (INS(i - 1, j)) { lap += (DW(i - 1, j) - c) * ihx2; }
61+ if (INS(i + 1, j)) { lap += (DW(i + 1, j) - c) * ihx2; }
62+ if (INS(i, j - 1)) { lap += (DW(i, j - 1) - c) * ihy2; }
63+ if (INS(i, j + 1)) { lap += (DW(i, j + 1) - c) * ihy2; }
64+ return lap;
65+}
66+
5667 // Weighted bending "moments" at nodes; zero where the second difference is undefined.
5768 fn Ma(i: i32, j: i32) -> f32 {
5869 if (!(INS(i - 1, j) && INS(i, j) && INS(i + 1, j))) { return 0.0; }
@@ -180,15 +191,20 @@ fn main(@builtin(global_invocation_id) gid: vec3u, @builtin(local_invocation_id)
180191 let mass = P.rhoh * P.hx * P.hy * weight;
181192 var acc = F * (1.0 / mass);
182193
183- // Loss: frequency-independent on all components, frequency-dependent (graph Laplacian) on w.
194+ // Loss: frequency-independent on all components; on w also the graph Laplacian of
195+ // the velocity (sigma1, ~k^2) and, when enabled, its biharmonic (sigma2, ~k^4).
184196 let vel = vec3f(velIn[k], velIn[n + k], velIn[2u * n + k]);
185- var lap = 0.0;
186- if (INS(i - 1, j)) { lap += (DW(i - 1, j) - vel.z) * ihx2; }
187- if (INS(i + 1, j)) { lap += (DW(i + 1, j) - vel.z) * ihx2; }
188- if (INS(i, j - 1)) { lap += (DW(i, j - 1) - vel.z) * ihy2; }
189- if (INS(i, j + 1)) { lap += (DW(i, j + 1) - vel.z) * ihy2; }
197+ let lap = lapDW(i, j);
190198 acc -= 2.0 * SP.sigma0 * vel;
191199 acc.z += 2.0 * SP.sigma1 * lap;
200+ if (SP.sigma2 > 0.0) {
201+ var bih = 0.0;
202+ if (INS(i - 1, j)) { bih += (lapDW(i - 1, j) - lap) * ihx2; }
203+ if (INS(i + 1, j)) { bih += (lapDW(i + 1, j) - lap) * ihx2; }
204+ if (INS(i, j - 1)) { bih += (lapDW(i, j - 1) - lap) * ihy2; }
205+ if (INS(i, j + 1)) { bih += (lapDW(i, j + 1) - lap) * ihy2; }
206+ acc.z -= 2.0 * SP.sigma2 * bih;
207+ }
192208 // In-plane (stretching) waves lie at several kHz; they get their own loss.
193209 acc.x -= 2.0 * SP.sigmaIn * vel.x;
194210 acc.y -= 2.0 * SP.sigmaIn * vel.y;
src/sim/simulation.tsmodified+3−1View file
@@ -365,10 +365,11 @@ export class Simulation {
365365 this.params.radiation = radiation;
366366 }
367367
368- setDamping(sigma0: number, sigma1: number, sigmaIn: number): void {
368+ setDamping(sigma0: number, sigma1: number, sigmaIn: number, sigma2 = 0): void {
369369 this.params.sigma0 = sigma0;
370370 this.params.sigma1 = sigma1;
371371 this.params.sigmaIn = sigmaIn;
372+ this.params.sigma2 = sigma2;
372373 }
373374
374375 setDuration(duration: number): void {
@@ -516,6 +517,7 @@ export class Simulation {
516517 f32[o + 8] = releaseGain;
517518 f32[o + 9] = p.sigmaIn;
518519 u32[o + 10] = parity;
520+ f32[o + 11] = p.sigma2 ?? 0;
519521 }
520522 this.device.queue.writeBuffer(this.stepBuf, 0, this.stepData, 0, steps * this.stepStride);
521523