Matching scripts: bands-per-octave option, linear mode frequencies, partial listing
6 changed files+127−5
.gitignoremodified+2−0View file
@@ -11,3 +11,5 @@ matching/*.npz
1111 matching/*.png
1212 matching/*.json
1313 matching/*.log
14+matching/*.bin
15+matching/*.txt
matching/README.mdmodified+3−0View file
@@ -7,6 +7,9 @@ This directory holds scripts for imposing the frequency-dependent decay of a rec
77 - `analyze.py WAV ...`: per-band levels, fitted decay times and the strongest partials of one or more clips.
88 - `match.py REFERENCE.wav SIMULATED.{wav,f32} OUT`: aligns both sounds at their onsets, measures band envelopes on a third-octave grid with a 4096-sample STFT, and writes `OUT_decay.wav` (the simulation with the recording's per-band decay rates, its own initial balance kept), `OUT_envelope.wav` (the recording's whole band envelope, balance and decay, imposed on the simulation), and the aligned inputs `OUT_reference_aligned.wav` and `OUT_sim_aligned.wav`. All outputs are level-matched on the body of the sound (RMS over 0.1 to 1 s) and the attack is soft-clipped.
99 - `plot.py OUT`: spectrograms of the four clips and band decay curves of the recording against the simulation.
10+- `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`.
11+
12+`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.
1013
1114 The simulated input should be the raw float pressure that `steelsim render` writes next to the wav (`NAME.f32`); a 16-bit wav loses the tail of the simulated sound after a few seconds.
1215
matching/match.pymodified+20−5View file
@@ -1,6 +1,9 @@
11 """Post-process a simulated sheet sound so that its frequency-dependent decay follows a recording.
22
3-Usage: python3 match.py REFERENCE.wav SIMULATED.wav OUT_PREFIX
3+Usage: python3 match.py REFERENCE.wav SIMULATED.wav OUT_PREFIX [--bands N] [--debug]
4+
5+--bands N sets the number of bands per octave (default 3, i.e. third-octave bands;
6+12 gives bands narrow enough to weight individual partials).
47
58 Both signals are aligned at their onsets and analysed with the same STFT. Band
69 envelopes (dB versus time) are formed on a logarithmic frequency grid. Two
@@ -21,9 +24,17 @@ from scipy.signal import stft, istft
2124
2225 FS = 48000
2326 NPERSEG, HOP = 4096, 512
24-# Log-spaced band edges, about a third of an octave, 25 Hz .. 16 kHz.
25-EDGES = 25 * 2 ** (np.arange(0, 28.1) / 3)
26-CENTERS = np.sqrt(EDGES[:-1] * EDGES[1:])
27+# Log-spaced band edges, 25 Hz .. 16 kHz, BANDS per octave (set by set_bands).
28+EDGES = CENTERS = None
29+
30+
31+def set_bands(per_octave):
32+ global EDGES, CENTERS
33+ EDGES = 25 * 2 ** (np.arange(0, 9.34 * per_octave + 0.5) / per_octave)
34+ CENTERS = np.sqrt(EDGES[:-1] * EDGES[1:])
35+
36+
37+set_bands(3)
2738
2839
2940 def load(path):
@@ -194,4 +205,8 @@ def main(ref_path, sim_path, out, debug=False):
194205
195206
196207 if __name__ == "__main__":
197- main(*sys.argv[1:4], debug="--debug" in sys.argv)
208+ args = [a for a in sys.argv[1:] if not a.startswith("--")]
209+ if "--bands" in sys.argv:
210+ set_bands(int(sys.argv[sys.argv.index("--bands") + 1]))
211+ args = [a for a in args if a != sys.argv[sys.argv.index("--bands") + 1]]
212+ main(*args[:3], debug="--debug" in sys.argv)
matching/modes.pyadded+17−0View file
@@ -0,0 +1,17 @@
1+"""Mode frequencies from the stiffness written by modes.ts.
2+
3+Usage: python3 modes.py K.bin META.json [fmax]
4+Prints the frequencies (Hz) of the modes below fmax, excluding the six rigid modes.
5+"""
6+import sys, json, numpy as np
7+
8+K = np.fromfile(sys.argv[1], dtype="<f8")
9+meta = json.load(open(sys.argv[2]))
10+n = meta["n"]
11+K = K.reshape(n, n)
12+K = 0.5 * (K + K.T)
13+ev = np.linalg.eigvalsh(K)
14+f = np.sqrt(np.clip(ev, 0, None)) / (2 * np.pi)
15+fmax = float(sys.argv[3]) if len(sys.argv) > 3 else 1500
16+modes = f[(f > 5) & (f < fmax)]
17+print(json.dumps({"params": meta["params"], "modes": [round(float(v), 1) for v in modes]}))
matching/modes.tsadded+47−0View file
@@ -0,0 +1,47 @@
1+// Linearized mode frequencies of the simulated sheet (with its initial bow), for
2+// comparing the model's partials with a recording without running the simulator.
3+//
4+// Usage: node modes.ts OUT.bin [key=value ...]
5+// Writes the mass-symmetrized stiffness of all three displacement fields over the
6+// nodes of the sheet as float64, and prints JSON with the degrees of freedom and the
7+// grid; eigenvalues are then omega^2 (see modes.py).
8+import { writeFileSync } from "node:fs";
9+import { accelerations, type Fields } from "../src/sim/cpuRef.ts";
10+import { defaultParams, initialShape, makeGrid, makeMaterial, nodeWeight, type PhysParams } from "../src/params.ts";
11+
12+const out = process.argv[2];
13+const p: PhysParams = { ...defaultParams };
14+for (const a of process.argv.slice(3)) {
15+ const [k, v] = a.split("=");
16+ const d = (p as unknown as Record<string, unknown>)[k];
17+ (p as unknown as Record<string, unknown>)[k] = typeof d === "number" ? parseFloat(v) : v;
18+}
19+const g = makeGrid(p);
20+const m = makeMaterial(p);
21+const w0 = Float64Array.from(initialShape(p, g));
22+const nodes: number[] = [];
23+for (let k = 0; k < g.n; k++) if (nodeWeight(g, k % g.nx, Math.floor(k / g.nx)) > 0) nodes.push(k);
24+const N = nodes.length;
25+const n = 3 * N;
26+const mass = (d: number) => m.rhoh * g.hx * g.hy * nodeWeight(g, nodes[d % N] % g.nx, Math.floor(nodes[d % N] / g.nx));
27+const x: Fields = { u: new Float64Array(g.n), v: new Float64Array(g.n), w: new Float64Array(g.n) };
28+const y: Fields = { u: new Float64Array(g.n), v: new Float64Array(g.n), w: new Float64Array(g.n) };
29+const y0: Fields = { u: new Float64Array(g.n), v: new Float64Array(g.n), w: new Float64Array(g.n) };
30+accelerations(g, m, x, y0, w0); // rest state (zero for a stress-free bow)
31+const eps = 1e-9;
32+const K = new Float64Array(n * n);
33+const arr = (f: Fields, c: number) => (c === 0 ? f.u : c === 1 ? f.v : f.w);
34+for (let col = 0; col < n; col++) {
35+ const c = Math.floor(col / N), k = nodes[col % N];
36+ arr(x, c)[k] = eps;
37+ accelerations(g, m, x, y, w0);
38+ arr(x, c)[k] = 0;
39+ const sc = Math.sqrt(1 / mass(col)) / eps;
40+ for (let row = 0; row < n; row++) {
41+ const rc = Math.floor(row / N), rk = nodes[row % N];
42+ const acc = arr(y, rc)[rk] - arr(y0, rc)[rk];
43+ K[row * n + col] = -acc * Math.sqrt(mass(row)) * sc;
44+ }
45+}
46+writeFileSync(out, Buffer.from(K.buffer));
47+console.log(JSON.stringify({ n, N, nx: g.nx, ny: g.ny, params: { E: p.E, nu: p.nu, rho: p.rho, h: p.h, bowShape: p.bowShape, bowAmplitude: p.bowAmplitude, spacing: p.spacing } }));
matching/partials.pyadded+38−0View file
@@ -0,0 +1,38 @@
1+"""Dominant partials (0.3 to 1.5 s after the onset) of several clips, side by side.
2+
3+Usage: python3 partials.py CLIP.wav|CLIP.f32 ...
4+"""
5+import sys, wave, numpy as np
6+from scipy.signal import find_peaks
7+
8+FS = 48000
9+
10+
11+def load(p):
12+ if p.endswith(".f32"):
13+ return np.fromfile(p, dtype="<f4").astype(float)
14+ w = wave.open(p)
15+ return np.frombuffer(w.readframes(w.getnframes()), dtype="<i2").astype(float) / 32768
16+
17+
18+def onset(x):
19+ e = np.abs(x)
20+ return max(0, int(np.argmax(e > 0.3 * e.max())) - int(0.005 * FS))
21+
22+
23+def partials(x, a=0.3, b=1.5, fmin=100, fmax=1300, n=12):
24+ s = x[int(a * FS):int(b * FS)]
25+ S = 20 * np.log10(np.abs(np.fft.rfft(s * np.hanning(len(s)))) + 1e-12)
26+ f = np.fft.rfftfreq(len(s), 1 / FS)
27+ m = (f > fmin) & (f < fmax)
28+ S = S - S[m].max()
29+ idx, _ = find_peaks(S, prominence=8, distance=int(3 / (FS / len(s))))
30+ idx = [i for i in idx if m[i]]
31+ idx = sorted(idx, key=lambda i: -S[i])[:n]
32+ return " ".join(f"{f[i]:.0f}({S[i]:.0f})" for i in sorted(idx))
33+
34+
35+for p in sys.argv[1:]:
36+ x = load(p)
37+ x = x[onset(x):]
38+ print(f"{p:36s} {partials(x)}")