Write raw float pressure from the render command; add post-processing scripts for matching a recording
6 changed files+316−1
.gitignoremodified+9−0View file
@@ -2,3 +2,12 @@ node_modules/
22 dist/
33 dist-cli/
44 docs/gallery/raw/
5+
6+# Post-processing experiments: keep the scripts, not the audio and figures
7+matching/*.wav
8+matching/*.f32
9+matching/*.npy
10+matching/*.npz
11+matching/*.png
12+matching/*.json
13+matching/*.log
matching/README.mdadded+17−0View file
@@ -0,0 +1,17 @@
1+# Matching a simulated sound to a recording by post-processing
2+
3+This directory holds scripts for imposing the frequency-dependent decay of a recorded sheet on a simulated sound, as a quick way to hear how far the loss model is from the real sheet before changing the model itself. The generated audio and figures are not committed.
4+
5+## Files
6+
7+- `analyze.py WAV ...`: per-band levels, fitted decay times and the strongest partials of one or more clips.
8+- `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.
9+- `plot.py OUT`: spectrograms of the four clips and band decay curves of the recording against the simulation.
10+
11+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.
12+
13+## Method and its limits
14+
15+For each band the recording's decay rate is fitted from 0.15 s until the band comes within 6 dB of its noise floor (estimated from the pre-onset segment and from the end of the recording). A band that reaches the floor within half a second gives only a lower bound on its decay rate, which is used as is; a band that sits at the floor from the start (the lowest bands of the recording, which are within a few decibels of the room noise) carries no information, and the decay rate of the nearest determined band is borrowed. The gain applied to the simulation is the difference of the two decay rates times the time since the onset (decay mode), or the difference of the two envelopes while the recording is above its floor and the extrapolated slope afterwards (envelope mode). The gains are interpolated across STFT bins linearly in log frequency and the sound is resynthesized by inverse STFT.
16+
17+The procedure changes envelopes, not partial frequencies: the modes of the simulated sheet stay where they are. It also cannot separate the attack from the early ring, so the first 100 ms are scaled like the early ring.
matching/analyze.pyadded+32−0View file
@@ -0,0 +1,32 @@
1+import sys, wave, numpy as np
2+from scipy.signal import stft
3+def load(path):
4+ w=wave.open(path); fs=w.getframerate(); n=w.getnframes()
5+ x=np.frombuffer(w.readframes(n),dtype='<i2').astype(float)/32768
6+ if w.getnchannels()>1: x=x.reshape(-1,w.getnchannels()).mean(1)
7+ return fs,x
8+def onset(x,fs):
9+ e=np.abs(x); i=np.argmax(e>0.3*e.max()); return max(0,i-int(0.005*fs))
10+bands=[(30,60),(60,120),(120,250),(250,500),(500,1000),(1000,2000),(2000,4000),(4000,8000),(8000,16000)]
11+def analyze(path):
12+ fs,x=load(path); o=onset(x,fs); x=x[o:]
13+ print(f'== {path}: fs {fs}, {len(x)/fs:.2f} s after onset, peak {np.abs(x).max():.3f} FS')
14+ f,t,Z=stft(x,fs,nperseg=4096,noverlap=4096-512)
15+ P=np.abs(Z)**2
16+ print(' band lvl@0.1s T60 fit(0.2-2s) T60 fit(0.2-end) level@2s')
17+ for lo,hi in bands:
18+ m=(f>=lo)&(f<hi); e=10*np.log10(P[m].sum(0)+1e-20)
19+ i0=np.searchsorted(t,0.1); e0=e[i0]
20+ def fit(a,b):
21+ ia,ib=np.searchsorted(t,a),np.searchsorted(t,b); ib=min(ib,len(t)-1)
22+ if ib-ia<4: return float('nan')
23+ k=np.polyfit(t[ia:ib],e[ia:ib],1)[0]; return -60/k if k<0 else float('inf')
24+ i2=min(np.searchsorted(t,2.0),len(t)-1)
25+ print(f' {lo:5d}-{hi:5d} {e0:7.1f} {fit(0.2,2.0):8.2f} {fit(0.2,t[-1]):8.2f} {e[i2]-e0:7.1f}')
26+ # strongest partials in 0.3-1.5 s
27+ s=x[int(0.3*fs):int(1.5*fs)]; S=np.abs(np.fft.rfft(s*np.hanning(len(s))))**2; fr=np.fft.rfftfreq(len(s),1/fs)
28+ idx=[i for i in range(2,len(S)-2) if S[i]>S[i-1] and S[i]>S[i+1] and fr[i]>25]
29+ idx=sorted(idx,key=lambda i:-S[i])[:14]
30+ print(' strongest partials (Hz, dB rel max):', ', '.join(f'{fr[i]:.0f}({10*np.log10(S[i]/S[idx[0]]):.0f})' for i in sorted(idx)))
31+ return fs,x
32+for p in sys.argv[1:]: analyze(p)
matching/match.pyadded+197−0View file
@@ -0,0 +1,197 @@
1+"""Post-process a simulated sheet sound so that its frequency-dependent decay follows a recording.
2+
3+Usage: python3 match.py REFERENCE.wav SIMULATED.wav OUT_PREFIX
4+
5+Both signals are aligned at their onsets and analysed with the same STFT. Band
6+envelopes (dB versus time) are formed on a logarithmic frequency grid. Two
7+outputs are written:
8+
9+ OUT_PREFIX_decay.wav per-band decay rates of the simulation are replaced by
10+ those of the recording; the initial spectral balance
11+ of the simulation is kept.
12+ OUT_PREFIX_envelope.wav the whole band envelope (initial balance and decay) of
13+ the recording is imposed on the simulation.
14+
15+The recording's noise floor is estimated from the segment before the onset;
16+where a band of the recording has decayed into that floor, the fitted decay
17+rate is used instead of the measured envelope so that noise is not copied.
18+"""
19+import sys, wave, numpy as np
20+from scipy.signal import stft, istft
21+
22+FS = 48000
23+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+
28+
29+def load(path):
30+ if path.endswith(".f32"):
31+ return np.fromfile(path, dtype="<f4").astype(float)
32+ w = wave.open(path)
33+ fs = w.getframerate()
34+ x = np.frombuffer(w.readframes(w.getnframes()), dtype="<i2").astype(float) / 32768
35+ if w.getnchannels() > 1:
36+ x = x.reshape(-1, w.getnchannels()).mean(1)
37+ if fs != FS:
38+ from scipy.signal import resample_poly
39+ from math import gcd
40+ g = gcd(FS, fs)
41+ x = resample_poly(x, FS // g, fs // g)
42+ return x
43+
44+
45+def save(path, x):
46+ x = np.clip(x, -1, 1)
47+ w = wave.open(path, "wb")
48+ w.setnchannels(1); w.setsampwidth(2); w.setframerate(FS)
49+ w.writeframes((x * 32767).astype("<i2").tobytes())
50+ w.close()
51+
52+
53+def onset(x):
54+ e = np.abs(x)
55+ i = int(np.argmax(e > 0.3 * e.max()))
56+ return max(0, i - int(0.005 * FS))
57+
58+
59+def band_envelopes(P, f):
60+ """dB power per band and frame, smoothed over about 100 ms."""
61+ E = np.zeros((len(CENTERS), P.shape[1]))
62+ for b in range(len(CENTERS)):
63+ m = (f >= EDGES[b]) & (f < EDGES[b + 1])
64+ E[b] = 10 * np.log10(P[m].sum(0) + 1e-30)
65+ k = max(1, int(0.1 * FS / HOP))
66+ ker = np.ones(k) / k
67+ for b in range(len(CENTERS)):
68+ E[b] = np.convolve(E[b], ker, mode="same")
69+ return E
70+
71+
72+def fit_slopes(E, t, floor, t0=0.15, tmax=3.0, sim_slopes=None):
73+ """Decay slope (dB/s) per band and how it was obtained.
74+
75+ "fit": fitted from t0 until the band first comes within 6 dB of its noise
76+ floor (or tmax), when that leaves at least 0.5 s.
77+ "bound": the band reaches the floor sooner but has dropped by at least 10 dB
78+ from t0: the slope over that drop is taken. The true decay is at
79+ least this fast.
80+ "unknown": the band sits at the floor from the start (the low bands of the
81+ recording are within a few decibels of the room noise): the
82+ simulation's own slope is kept, i.e. nothing is imposed.
83+ """
84+ slopes = np.zeros(len(CENTERS))
85+ status = ["unknown"] * len(CENTERS)
86+ i0 = np.searchsorted(t, t0)
87+ for b in range(len(CENTERS)):
88+ i1 = min(np.searchsorted(t, tmax), len(t) - 1)
89+ low = np.where(E[b, i0:] < floor[b] + 6)[0]
90+ if len(low):
91+ i1 = min(i1, i0 + low[0])
92+ if t[i1] - t[i0] >= 0.5:
93+ slopes[b] = np.polyfit(t[i0:i1], E[b, i0:i1], 1)[0]
94+ status[b] = "fit"
95+ elif i1 > i0 + 2 and E[b, i0] - E[b, i1] >= 10:
96+ slopes[b] = (E[b, i1] - E[b, i0]) / (t[i1] - t[i0])
97+ status[b] = "bound"
98+ # Bands without information borrow the slope of the nearest determined band; if
99+ # none exists, they keep the simulation's own slope.
100+ known = [b for b in range(len(CENTERS)) if status[b] != "unknown"]
101+ for b in range(len(CENTERS)):
102+ if status[b] == "unknown":
103+ slopes[b] = slopes[known[np.argmin(np.abs(np.array(known) - b))]] if known else (sim_slopes[b] if sim_slopes is not None else 0.0)
104+ return slopes, status
105+
106+
107+def main(ref_path, sim_path, out, debug=False):
108+ ref = load(ref_path)
109+ sim = load(sim_path)
110+ o_ref, o_sim = onset(ref), onset(sim)
111+ noise = ref[: max(int(0.2 * FS), o_ref - int(0.05 * FS))]
112+ ref = ref[o_ref:]
113+ sim = sim[o_sim:]
114+ ref_full = ref
115+ n = min(len(ref), len(sim))
116+ ref, sim = ref[:n], sim[:n]
117+
118+ f, t, Zr = stft(ref, FS, nperseg=NPERSEG, noverlap=NPERSEG - HOP)
119+ _, _, Zs = stft(sim, FS, nperseg=NPERSEG, noverlap=NPERSEG - HOP)
120+ Er = band_envelopes(np.abs(Zr) ** 2, f)
121+ Es = band_envelopes(np.abs(Zs) ** 2, f)
122+ # Noise floor of the recording per band (median over the pre-onset segment).
123+ _, _, Zn = stft(noise, FS, nperseg=NPERSEG, noverlap=NPERSEG - HOP)
124+ En = band_envelopes(np.abs(Zn) ** 2, f)
125+ # The pre-onset segment gives the floor; the last half second of the recording,
126+ # where the high bands have long since died away, gives a second estimate.
127+ _, _, Zt = stft(ref_full[-int(0.5 * FS) :], FS, nperseg=NPERSEG, noverlap=NPERSEG - HOP)
128+ Et = band_envelopes(np.abs(Zt) ** 2, f)
129+ floor = np.maximum(np.median(En, axis=1) + 3, np.median(Et, axis=1) - 3)
130+
131+ prefloor = np.median(En, axis=1)
132+ ks, _ = fit_slopes(Es, t, np.full(len(CENTERS), -400.0))
133+ kr, status = fit_slopes(Er, t, floor, sim_slopes=ks)
134+
135+ print("band centre T60 ref T60 sim ref fit ref level@0.1s pre-onset noise")
136+ for b in range(len(CENTERS)):
137+ t60 = lambda k: (-60 / k if k < 0 else float("inf"))
138+ print(f" {CENTERS[b]:7.0f} Hz {t60(kr[b]):7.2f} {t60(ks[b]):7.2f} {status[b]:8s} {Er[b, np.searchsorted(t, 0.1)]:7.1f} {prefloor[b]:7.1f}")
139+
140+ # Gain in dB per band and frame.
141+ G_decay = (kr - ks)[:, None] * t[None, :]
142+ i_ref = np.searchsorted(t, 0.1)
143+ # Envelope mode: follow the recording while it is above its floor, then continue
144+ # along its fitted slope from the last trusted frame.
145+ G_env = Er - Es
146+ for b in range(len(CENTERS)):
147+ if status[b] == "unknown":
148+ # Only the balance at 0.1 s is known; apply the borrowed decay after that.
149+ G_env[b, i_ref:] = G_env[b, i_ref] + (kr[b] - ks[b]) * (t[i_ref:] - t[i_ref])
150+ continue
151+ # Follow the recording until the band first comes within 10 dB of its floor
152+ # (later frames may contain room noises), then continue along its fitted decay.
153+ below = np.where(Er[b, i_ref:] < floor[b] + 10)[0]
154+ if len(below) == 0:
155+ continue
156+ last = i_ref + below[0] - 1
157+ G_env[b, last + 1 :] = G_env[b, last] + (kr[b] - ks[b]) * (t[last + 1 :] - t[last])
158+ # The attack (first 100 ms) is scaled like the early ring rather than frame by frame:
159+ # the simulated strike has a much sharper click than the recording, and the STFT
160+ # frames cannot separate the two anyway.
161+ G_env[:, :i_ref] = G_env[:, [i_ref]]
162+ G_decay = np.clip(G_decay, -60, 60)
163+ G_env = np.clip(G_env, -60, 60)
164+
165+ def apply(G):
166+ # Interpolate the band gains onto the STFT bins (linear in log frequency).
167+ lf = np.log(np.maximum(f, 1.0))
168+ lc = np.log(CENTERS)
169+ Gf = np.zeros(Zs.shape)
170+ for j in range(Zs.shape[1]):
171+ Gf[:, j] = np.interp(lf, lc, G[:, j])
172+ Y = Zs * 10 ** (Gf / 20)
173+ _, y = istft(Y, FS, nperseg=NPERSEG, noverlap=NPERSEG - HOP)
174+ return y[:n]
175+
176+ y_decay = apply(G_decay)
177+ y_env = apply(G_env)
178+
179+ # Output level: the body of each sound (RMS over 0.1 to 1 s) is set to that of the
180+ # recording, and the attack is soft-clipped like the recording's own clipped attack.
181+ i_a, i_b = int(0.1 * FS), int(1.0 * FS)
182+ body = lambda x: np.sqrt(np.mean(x[i_a:i_b] ** 2))
183+ target = body(ref)
184+ def finish(x):
185+ x = x * (target / body(x))
186+ return np.tanh(x)
187+ for name, y in (("decay", y_decay), ("envelope", y_env), ("sim_aligned", sim)):
188+ save(f"{out}_{name}.wav", finish(y))
189+ if debug:
190+ np.save(f"{out}_{name}_scaled.npy", y * (target / body(y)))
191+ save(f"{out}_reference_aligned.wav", ref)
192+ print(f"body RMS of the recording (0.1-1 s): {20 * np.log10(target):.1f} dBFS; outputs level-matched there, attack soft-clipped")
193+ np.savez(f"{out}_analysis.npz", f=f, t=t, Er=Er, Es=Es, floor=floor, kr=kr, ks=ks, centers=CENTERS)
194+
195+
196+if __name__ == "__main__":
197+ main(*sys.argv[1:4], debug="--debug" in sys.argv)
matching/plot.pyadded+58−0View file
@@ -0,0 +1,58 @@
1+"""Spectrogram comparison of the recording, the simulation and the matched outputs.
2+
3+Usage: python3 plot.py OUT_PREFIX (reads the wavs written by match.py)
4+"""
5+import sys, wave, numpy as np
6+import matplotlib
7+matplotlib.use("Agg")
8+import matplotlib.pyplot as plt
9+from scipy.signal import stft
10+
11+FS = 48000
12+
13+
14+def load(path):
15+ w = wave.open(path)
16+ return np.frombuffer(w.readframes(w.getnframes()), dtype="<i2").astype(float) / 32768
17+
18+
19+out = sys.argv[1]
20+files = [("recording", f"{out}_reference_aligned.wav"), ("simulation", f"{out}_sim_aligned.wav"),
21+ ("matched: decay rates", f"{out}_decay.wav"), ("matched: full envelope", f"{out}_envelope.wav")]
22+fig, axes = plt.subplots(2, 2, figsize=(13, 8), sharex=True, sharey=True)
23+for ax, (title, path) in zip(axes.flat, files):
24+ x = load(path)
25+ f, t, Z = stft(x, FS, nperseg=4096, noverlap=4096 - 512)
26+ S = 20 * np.log10(np.abs(Z) + 1e-9)
27+ S -= S.max()
28+ ax.pcolormesh(t, f, S, vmin=-90, vmax=0, cmap="magma", shading="auto")
29+ ax.set_yscale("log")
30+ ax.set_ylim(25, 12000)
31+ ax.set_title(title)
32+ ax.set_xlabel("time [s]")
33+ ax.set_ylabel("frequency [Hz]")
34+fig.suptitle("Reference hit versus simulated strike (dB relative to each clip's peak, 90 dB range)")
35+fig.tight_layout()
36+fig.savefig(f"{out}_spectrograms.png", dpi=110)
37+print(f"wrote {out}_spectrograms.png")
38+
39+# Band decay curves for a few bands.
40+a = np.load(f"{out}_analysis.npz")
41+t, Er, Es, centers, floor = a["t"], a["Er"], a["Es"], a["centers"], a["floor"]
42+fig, axes = plt.subplots(2, 4, figsize=(14, 6), sharex=True)
43+for ax, fc in zip(axes.flat, [35, 89, 178, 356, 713, 1425, 2851, 5702]):
44+ b = int(np.argmin(np.abs(centers - fc)))
45+ ax.plot(t, Er[b] - Er[b, np.searchsorted(t, 0.1)], label="recording")
46+ ax.plot(t, Es[b] - Es[b, np.searchsorted(t, 0.1)], label="simulation")
47+ ax.axhline(floor[b] - Er[b, np.searchsorted(t, 0.1)], color="gray", ls=":", lw=1, label="recording noise floor")
48+ ax.set_title(f"{centers[b]:.0f} Hz band")
49+ ax.set_ylim(-80, 10)
50+ ax.grid(alpha=0.3)
51+axes[0, 0].legend(fontsize=8)
52+for ax in axes[1]:
53+ ax.set_xlabel("time after onset [s]")
54+for ax in axes[:, 0]:
55+ ax.set_ylabel("level rel. 0.1 s [dB]")
56+fig.tight_layout()
57+fig.savefig(f"{out}_band_decay.png", dpi=110)
58+print(f"wrote {out}_band_decay.png")
src/cli/main.tsmodified+3−1View file
@@ -1,6 +1,6 @@
11 // Command-line renderer.
22 //
3-// steelsim render [--out DIR] [--name ID] [--gain DB] key=value ...
3+// steelsim render [--out DIR] [--name ID] [--gain DB] key=value ... (also writes ID.f32, the raw pressure)
44 // steelsim gallery MANIFEST.json [--out DIR] [--only SWEEP_ID] [--force] [--html-only] [key=value ...]
55 //
66 // Common options: --engine chrome|dawn (default chrome), --chrome PATH (Chrome
@@ -100,6 +100,8 @@ async function render(): Promise<void> {
100100 const r = await runSimulation(overrides as Partial<PhysParams>, console.log);
101101 console.log(" simulation finished");
102102 const rec = writeClip(dir, id, id, overrides, r, gainDb);
103+ // Raw listener pressure [Pa] as little-endian float32, for post-processing without the 16-bit floor.
104+ writeFileSync(join(dir, `${id}.f32`), Buffer.from(r.audio.buffer, r.audio.byteOffset, r.audio.byteLength));
103105 console.log(" files written");
104106 console.log(` ${r.grid.nx}x${r.grid.ny} nodes, ${r.steps} steps in ${r.seconds.toFixed(1)} s (${(r.steps / r.seconds / 1e3).toFixed(0)} k steps/s)`);
105107 console.log(` peak ${rec.metrics.peakPa.toFixed(3)} Pa (${rec.metrics.peakDbSpl.toFixed(1)} dB SPL), centroid ${rec.metrics.centroidHz.toFixed(0)} Hz, wav peak ${rec.peakDbfs.toFixed(1)} dBFS`);