"""Post-process a simulated sheet sound so that its frequency-dependent decay follows a recording. Usage: python3 match.py REFERENCE.wav SIMULATED.wav OUT_PREFIX Both signals are aligned at their onsets and analysed with the same STFT. Band envelopes (dB versus time) are formed on a logarithmic frequency grid. Two outputs are written: OUT_PREFIX_decay.wav per-band decay rates of the simulation are replaced by those of the recording; the initial spectral balance of the simulation is kept. OUT_PREFIX_envelope.wav the whole band envelope (initial balance and decay) of the recording is imposed on the simulation. The recording's noise floor is estimated from the segment before the onset; where a band of the recording has decayed into that floor, the fitted decay rate is used instead of the measured envelope so that noise is not copied. """ import sys, wave, numpy as np from scipy.signal import stft, istft FS = 48000 NPERSEG, HOP = 4096, 512 # Log-spaced band edges, about a third of an octave, 25 Hz .. 16 kHz. EDGES = 25 * 2 ** (np.arange(0, 28.1) / 3) CENTERS = np.sqrt(EDGES[:-1] * EDGES[1:]) def load(path): if path.endswith(".f32"): return np.fromfile(path, dtype=" 1: x = x.reshape(-1, w.getnchannels()).mean(1) if fs != FS: from scipy.signal import resample_poly from math import gcd g = gcd(FS, fs) x = resample_poly(x, FS // g, fs // g) return x def save(path, x): x = np.clip(x, -1, 1) w = wave.open(path, "wb") w.setnchannels(1); w.setsampwidth(2); w.setframerate(FS) w.writeframes((x * 32767).astype(" 0.3 * e.max())) return max(0, i - int(0.005 * FS)) def band_envelopes(P, f): """dB power per band and frame, smoothed over about 100 ms.""" E = np.zeros((len(CENTERS), P.shape[1])) for b in range(len(CENTERS)): m = (f >= EDGES[b]) & (f < EDGES[b + 1]) E[b] = 10 * np.log10(P[m].sum(0) + 1e-30) k = max(1, int(0.1 * FS / HOP)) ker = np.ones(k) / k for b in range(len(CENTERS)): E[b] = np.convolve(E[b], ker, mode="same") return E def fit_slopes(E, t, floor, t0=0.15, tmax=3.0, sim_slopes=None): """Decay slope (dB/s) per band and how it was obtained. "fit": fitted from t0 until the band first comes within 6 dB of its noise floor (or tmax), when that leaves at least 0.5 s. "bound": the band reaches the floor sooner but has dropped by at least 10 dB from t0: the slope over that drop is taken. The true decay is at least this fast. "unknown": the band sits at the floor from the start (the low bands of the recording are within a few decibels of the room noise): the simulation's own slope is kept, i.e. nothing is imposed. """ slopes = np.zeros(len(CENTERS)) status = ["unknown"] * len(CENTERS) i0 = np.searchsorted(t, t0) for b in range(len(CENTERS)): i1 = min(np.searchsorted(t, tmax), len(t) - 1) low = np.where(E[b, i0:] < floor[b] + 6)[0] if len(low): i1 = min(i1, i0 + low[0]) if t[i1] - t[i0] >= 0.5: slopes[b] = np.polyfit(t[i0:i1], E[b, i0:i1], 1)[0] status[b] = "fit" elif i1 > i0 + 2 and E[b, i0] - E[b, i1] >= 10: slopes[b] = (E[b, i1] - E[b, i0]) / (t[i1] - t[i0]) status[b] = "bound" # Bands without information borrow the slope of the nearest determined band; if # none exists, they keep the simulation's own slope. known = [b for b in range(len(CENTERS)) if status[b] != "unknown"] for b in range(len(CENTERS)): if status[b] == "unknown": 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) return slopes, status def main(ref_path, sim_path, out, debug=False): ref = load(ref_path) sim = load(sim_path) o_ref, o_sim = onset(ref), onset(sim) noise = ref[: max(int(0.2 * FS), o_ref - int(0.05 * FS))] ref = ref[o_ref:] sim = sim[o_sim:] ref_full = ref n = min(len(ref), len(sim)) ref, sim = ref[:n], sim[:n] f, t, Zr = stft(ref, FS, nperseg=NPERSEG, noverlap=NPERSEG - HOP) _, _, Zs = stft(sim, FS, nperseg=NPERSEG, noverlap=NPERSEG - HOP) Er = band_envelopes(np.abs(Zr) ** 2, f) Es = band_envelopes(np.abs(Zs) ** 2, f) # Noise floor of the recording per band (median over the pre-onset segment). _, _, Zn = stft(noise, FS, nperseg=NPERSEG, noverlap=NPERSEG - HOP) En = band_envelopes(np.abs(Zn) ** 2, f) # The pre-onset segment gives the floor; the last half second of the recording, # where the high bands have long since died away, gives a second estimate. _, _, Zt = stft(ref_full[-int(0.5 * FS) :], FS, nperseg=NPERSEG, noverlap=NPERSEG - HOP) Et = band_envelopes(np.abs(Zt) ** 2, f) floor = np.maximum(np.median(En, axis=1) + 3, np.median(Et, axis=1) - 3) prefloor = np.median(En, axis=1) ks, _ = fit_slopes(Es, t, np.full(len(CENTERS), -400.0)) kr, status = fit_slopes(Er, t, floor, sim_slopes=ks) print("band centre T60 ref T60 sim ref fit ref level@0.1s pre-onset noise") for b in range(len(CENTERS)): t60 = lambda k: (-60 / k if k < 0 else float("inf")) 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}") # Gain in dB per band and frame. G_decay = (kr - ks)[:, None] * t[None, :] i_ref = np.searchsorted(t, 0.1) # Envelope mode: follow the recording while it is above its floor, then continue # along its fitted slope from the last trusted frame. G_env = Er - Es for b in range(len(CENTERS)): if status[b] == "unknown": # Only the balance at 0.1 s is known; apply the borrowed decay after that. G_env[b, i_ref:] = G_env[b, i_ref] + (kr[b] - ks[b]) * (t[i_ref:] - t[i_ref]) continue # Follow the recording until the band first comes within 10 dB of its floor # (later frames may contain room noises), then continue along its fitted decay. below = np.where(Er[b, i_ref:] < floor[b] + 10)[0] if len(below) == 0: continue last = i_ref + below[0] - 1 G_env[b, last + 1 :] = G_env[b, last] + (kr[b] - ks[b]) * (t[last + 1 :] - t[last]) # The attack (first 100 ms) is scaled like the early ring rather than frame by frame: # the simulated strike has a much sharper click than the recording, and the STFT # frames cannot separate the two anyway. G_env[:, :i_ref] = G_env[:, [i_ref]] G_decay = np.clip(G_decay, -60, 60) G_env = np.clip(G_env, -60, 60) def apply(G): # Interpolate the band gains onto the STFT bins (linear in log frequency). lf = np.log(np.maximum(f, 1.0)) lc = np.log(CENTERS) Gf = np.zeros(Zs.shape) for j in range(Zs.shape[1]): Gf[:, j] = np.interp(lf, lc, G[:, j]) Y = Zs * 10 ** (Gf / 20) _, y = istft(Y, FS, nperseg=NPERSEG, noverlap=NPERSEG - HOP) return y[:n] y_decay = apply(G_decay) y_env = apply(G_env) # Output level: the body of each sound (RMS over 0.1 to 1 s) is set to that of the # recording, and the attack is soft-clipped like the recording's own clipped attack. i_a, i_b = int(0.1 * FS), int(1.0 * FS) body = lambda x: np.sqrt(np.mean(x[i_a:i_b] ** 2)) target = body(ref) def finish(x): x = x * (target / body(x)) return np.tanh(x) for name, y in (("decay", y_decay), ("envelope", y_env), ("sim_aligned", sim)): save(f"{out}_{name}.wav", finish(y)) if debug: np.save(f"{out}_{name}_scaled.npy", y * (target / body(y))) save(f"{out}_reference_aligned.wav", ref) print(f"body RMS of the recording (0.1-1 s): {20 * np.log10(target):.1f} dBFS; outputs level-matched there, attack soft-clipped") np.savez(f"{out}_analysis.npz", f=f, t=t, Er=Er, Es=Es, floor=floor, kr=kr, ks=ks, centers=CENTERS) if __name__ == "__main__": main(*sys.argv[1:4], debug="--debug" in sys.argv)