5200b59Write raw float pressure from the render command; add post-processing scripts for matching a recordingJeremy Magland 1"""Post-process a simulated sheet sound so that its frequency-dependent decay follows a recording.
98b9eb0Matching scripts: bands-per-octave option, linear mode frequencies, partial listingJeremy Magland 3Usage: python3 match.py REFERENCE.wav SIMULATED.wav OUT_PREFIX [--bands N] [--debug]
5--bands N sets the number of bands per octave (default 3, i.e. third-octave bands;
612 gives bands narrow enough to weight individual partials).
5200b59Write raw float pressure from the render command; add post-processing scripts for matching a recordingJeremy Magland 7
8Both signals are aligned at their onsets and analysed with the same STFT. Band
9envelopes (dB versus time) are formed on a logarithmic frequency grid. Two
10outputs are written:
12 OUT_PREFIX_decay.wav per-band decay rates of the simulation are replaced by
13 those of the recording; the initial spectral balance
14 of the simulation is kept.
15 OUT_PREFIX_envelope.wav the whole band envelope (initial balance and decay) of
16 the recording is imposed on the simulation.
18The recording's noise floor is estimated from the segment before the onset;
19where a band of the recording has decayed into that floor, the fitted decay
20rate is used instead of the measured envelope so that noise is not copied.
21"""
22import sys, wave, numpy as np
23from scipy.signal import stft, istft
25FS = 48000
26NPERSEG, HOP = 4096, 512
98b9eb0Matching scripts: bands-per-octave option, linear mode frequencies, partial listingJeremy Magland 27# Log-spaced band edges, 25 Hz .. 16 kHz, BANDS per octave (set by set_bands).
28EDGES = CENTERS = None
31def 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:])
37set_bands(3)
5200b59Write raw float pressure from the render command; add post-processing scripts for matching a recordingJeremy Magland 38
40def load(path):
41 if path.endswith(".f32"):
42 return np.fromfile(path, dtype="<f4").astype(float)
43 w = wave.open(path)
44 fs = w.getframerate()
45 x = np.frombuffer(w.readframes(w.getnframes()), dtype="<i2").astype(float) / 32768
46 if w.getnchannels() > 1:
47 x = x.reshape(-1, w.getnchannels()).mean(1)
48 if fs != FS:
49 from scipy.signal import resample_poly
50 from math import gcd
51 g = gcd(FS, fs)
52 x = resample_poly(x, FS // g, fs // g)
53 return x
56def save(path, x):
57 x = np.clip(x, -1, 1)
58 w = wave.open(path, "wb")
59 w.setnchannels(1); w.setsampwidth(2); w.setframerate(FS)
60 w.writeframes((x * 32767).astype("<i2").tobytes())
61 w.close()
64def onset(x):
65 e = np.abs(x)
66 i = int(np.argmax(e > 0.3 * e.max()))
67 return max(0, i - int(0.005 * FS))
70def band_envelopes(P, f):
71 """dB power per band and frame, smoothed over about 100 ms."""
72 E = np.zeros((len(CENTERS), P.shape[1]))
73 for b in range(len(CENTERS)):
74 m = (f >= EDGES[b]) & (f < EDGES[b + 1])
75 E[b] = 10 * np.log10(P[m].sum(0) + 1e-30)
76 k = max(1, int(0.1 * FS / HOP))
77 ker = np.ones(k) / k
78 for b in range(len(CENTERS)):
79 E[b] = np.convolve(E[b], ker, mode="same")
80 return E
83def fit_slopes(E, t, floor, t0=0.15, tmax=3.0, sim_slopes=None):
84 """Decay slope (dB/s) per band and how it was obtained.
86 "fit": fitted from t0 until the band first comes within 6 dB of its noise
87 floor (or tmax), when that leaves at least 0.5 s.
88 "bound": the band reaches the floor sooner but has dropped by at least 10 dB
89 from t0: the slope over that drop is taken. The true decay is at
90 least this fast.
91 "unknown": the band sits at the floor from the start (the low bands of the
92 recording are within a few decibels of the room noise): the
93 simulation's own slope is kept, i.e. nothing is imposed.
94 """
95 slopes = np.zeros(len(CENTERS))
96 status = ["unknown"] * len(CENTERS)
97 i0 = np.searchsorted(t, t0)
98 for b in range(len(CENTERS)):
99 i1 = min(np.searchsorted(t, tmax), len(t) - 1)
100 low = np.where(E[b, i0:] < floor[b] + 6)[0]
101 if len(low):
102 i1 = min(i1, i0 + low[0])
103 if t[i1] - t[i0] >= 0.5:
104 slopes[b] = np.polyfit(t[i0:i1], E[b, i0:i1], 1)[0]
105 status[b] = "fit"
106 elif i1 > i0 + 2 and E[b, i0] - E[b, i1] >= 10:
107 slopes[b] = (E[b, i1] - E[b, i0]) / (t[i1] - t[i0])
108 status[b] = "bound"
109 # Bands without information borrow the slope of the nearest determined band; if
110 # none exists, they keep the simulation's own slope.
111 known = [b for b in range(len(CENTERS)) if status[b] != "unknown"]
112 for b in range(len(CENTERS)):
113 if status[b] == "unknown":
114 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)
115 return slopes, status
118def main(ref_path, sim_path, out, debug=False):
119 ref = load(ref_path)
120 sim = load(sim_path)
121 o_ref, o_sim = onset(ref), onset(sim)
122 noise = ref[: max(int(0.2 * FS), o_ref - int(0.05 * FS))]
123 ref = ref[o_ref:]
124 sim = sim[o_sim:]
125 ref_full = ref
126 n = min(len(ref), len(sim))
127 ref, sim = ref[:n], sim[:n]
129 f, t, Zr = stft(ref, FS, nperseg=NPERSEG, noverlap=NPERSEG - HOP)
130 _, _, Zs = stft(sim, FS, nperseg=NPERSEG, noverlap=NPERSEG - HOP)
131 Er = band_envelopes(np.abs(Zr) ** 2, f)
132 Es = band_envelopes(np.abs(Zs) ** 2, f)
133 # Noise floor of the recording per band (median over the pre-onset segment).
134 _, _, Zn = stft(noise, FS, nperseg=NPERSEG, noverlap=NPERSEG - HOP)
135 En = band_envelopes(np.abs(Zn) ** 2, f)
136 # The pre-onset segment gives the floor; the last half second of the recording,
137 # where the high bands have long since died away, gives a second estimate.
138 _, _, Zt = stft(ref_full[-int(0.5 * FS) :], FS, nperseg=NPERSEG, noverlap=NPERSEG - HOP)
139 Et = band_envelopes(np.abs(Zt) ** 2, f)
140 floor = np.maximum(np.median(En, axis=1) + 3, np.median(Et, axis=1) - 3)
142 prefloor = np.median(En, axis=1)
143 ks, _ = fit_slopes(Es, t, np.full(len(CENTERS), -400.0))
144 kr, status = fit_slopes(Er, t, floor, sim_slopes=ks)
146 print("band centre T60 ref T60 sim ref fit ref level@0.1s pre-onset noise")
147 for b in range(len(CENTERS)):
148 t60 = lambda k: (-60 / k if k < 0 else float("inf"))
149 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}")
151 # Gain in dB per band and frame.
152 G_decay = (kr - ks)[:, None] * t[None, :]
153 i_ref = np.searchsorted(t, 0.1)
154 # Envelope mode: follow the recording while it is above its floor, then continue
155 # along its fitted slope from the last trusted frame.
156 G_env = Er - Es
157 for b in range(len(CENTERS)):
158 if status[b] == "unknown":
159 # Only the balance at 0.1 s is known; apply the borrowed decay after that.
160 G_env[b, i_ref:] = G_env[b, i_ref] + (kr[b] - ks[b]) * (t[i_ref:] - t[i_ref])
161 continue
162 # Follow the recording until the band first comes within 10 dB of its floor
163 # (later frames may contain room noises), then continue along its fitted decay.
164 below = np.where(Er[b, i_ref:] < floor[b] + 10)[0]
165 if len(below) == 0:
166 continue
167 last = i_ref + below[0] - 1
168 G_env[b, last + 1 :] = G_env[b, last] + (kr[b] - ks[b]) * (t[last + 1 :] - t[last])
169 # The attack (first 100 ms) is scaled like the early ring rather than frame by frame:
170 # the simulated strike has a much sharper click than the recording, and the STFT
171 # frames cannot separate the two anyway.
172 G_env[:, :i_ref] = G_env[:, [i_ref]]
173 G_decay = np.clip(G_decay, -60, 60)
174 G_env = np.clip(G_env, -60, 60)
176 def apply(G):
177 # Interpolate the band gains onto the STFT bins (linear in log frequency).
178 lf = np.log(np.maximum(f, 1.0))
179 lc = np.log(CENTERS)
180 Gf = np.zeros(Zs.shape)
181 for j in range(Zs.shape[1]):
182 Gf[:, j] = np.interp(lf, lc, G[:, j])
183 Y = Zs * 10 ** (Gf / 20)
184 _, y = istft(Y, FS, nperseg=NPERSEG, noverlap=NPERSEG - HOP)
185 return y[:n]
187 y_decay = apply(G_decay)
188 y_env = apply(G_env)
190 # Output level: the body of each sound (RMS over 0.1 to 1 s) is set to that of the
191 # recording, and the attack is soft-clipped like the recording's own clipped attack.
192 i_a, i_b = int(0.1 * FS), int(1.0 * FS)
193 body = lambda x: np.sqrt(np.mean(x[i_a:i_b] ** 2))
194 target = body(ref)
195 def finish(x):
196 x = x * (target / body(x))
197 return np.tanh(x)
198 for name, y in (("decay", y_decay), ("envelope", y_env), ("sim_aligned", sim)):
199 save(f"{out}_{name}.wav", finish(y))
200 if debug:
201 np.save(f"{out}_{name}_scaled.npy", y * (target / body(y)))
202 save(f"{out}_reference_aligned.wav", ref)
203 print(f"body RMS of the recording (0.1-1 s): {20 * np.log10(target):.1f} dBFS; outputs level-matched there, attack soft-clipped")
204 np.savez(f"{out}_analysis.npz", f=f, t=t, Er=Er, Es=Es, floor=floor, kr=kr, ks=ks, centers=CENTERS)
207if __name__ == "__main__":
98b9eb0Matching scripts: bands-per-octave option, linear mode frequencies, partial listingJeremy Magland 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)