1"""Post-process a simulated sheet sound so that its frequency-dependent decay follows a recording.
3Usage: python3 match.py REFERENCE.wav SIMULATED.wav OUT_PREFIX
5Both signals are aligned at their onsets and analysed with the same STFT. Band
6envelopes (dB versus time) are formed on a logarithmic frequency grid. Two
7outputs are written:
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.
15The recording's noise floor is estimated from the segment before the onset;
16where a band of the recording has decayed into that floor, the fitted decay
17rate is used instead of the measured envelope so that noise is not copied.
18"""
19import sys, wave, numpy as np
20from scipy.signal import stft, istft
22FS = 48000
23NPERSEG, HOP = 4096, 512
24# Log-spaced band edges, about a third of an octave, 25 Hz .. 16 kHz.
25EDGES = 25 * 2 ** (np.arange(0, 28.1) / 3)
26CENTERS = np.sqrt(EDGES[:-1] * EDGES[1:])
29def 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
45def 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()
53def 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))
59def 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
72def 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.
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
107def 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]
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)
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)
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}")
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)
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]
176 y_decay = apply(G_decay)
177 y_env = apply(G_env)
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)
196if __name__ == "__main__":
197 main(*sys.argv[1:4], debug="--debug" in sys.argv)