1"""How much of the envelope match is a static equalizer?
3Usage: python3 eq.py REFERENCE.wav SIMULATED.f32 OUT_PREFIX
5Recomputes the envelope-mode gain table of match.py, splits it into a static part
6(the per-band gain at 0.1 s) and a time-varying remainder, reports the size of the
7remainder, and writes
8 OUT_eq.wav the simulation through the static per-band equalizer (one gain per band)
9 OUT_eq4.wav the same with the equalizer smoothed to a cubic in log frequency (4 numbers)
10"""
11import sys, numpy as np
12from scipy.signal import stft, istft
13import match as M
15ref_path, sim_path, out = sys.argv[1:4]
16FS, NPERSEG, HOP = M.FS, M.NPERSEG, M.HOP
17ref = M.load(ref_path)
18sim = M.load(sim_path)
19o_ref, o_sim = M.onset(ref), M.onset(sim)
20noise = ref[: max(int(0.2 * FS), o_ref - int(0.05 * FS))]
21ref, sim = ref[o_ref:], sim[o_sim:]
22ref_full = ref
23n = min(len(ref), len(sim))
24ref, sim = ref[:n], sim[:n]
25f, t, Zr = stft(ref, FS, nperseg=NPERSEG, noverlap=NPERSEG - HOP)
26_, _, Zs = stft(sim, FS, nperseg=NPERSEG, noverlap=NPERSEG - HOP)
27Er, Es = M.band_envelopes(np.abs(Zr) ** 2, f), M.band_envelopes(np.abs(Zs) ** 2, f)
28_, _, Zn = stft(noise, FS, nperseg=NPERSEG, noverlap=NPERSEG - HOP)
29En = M.band_envelopes(np.abs(Zn) ** 2, f)
30_, _, Zt = stft(ref_full[-int(0.5 * FS):], FS, nperseg=NPERSEG, noverlap=NPERSEG - HOP)
31Et = M.band_envelopes(np.abs(Zt) ** 2, f)
32floor = np.maximum(np.median(En, axis=1) + 3, np.median(Et, axis=1) - 3)
33i_ref = np.searchsorted(t, 0.1)
34C = M.CENTERS
36G0 = (Er - Es)[:, i_ref] # static part: balance at 0.1 s
37print("static equalizer (dB gain applied to the simulation per third-octave band):")
38for c, g in zip(C, G0):
39 print(f" {c:7.0f} Hz {g:+6.1f}")
40# Time-varying remainder over the frames where the recording is above its floor.
41print("\nremainder after the static part (recording above floor + 10 dB):")
42print(" band trusted until remainder RMS remainder at end of trusted range")
43for b in range(len(C)):
44 below = np.where(Er[b, i_ref:] < floor[b] + 10)[0]
45 last = i_ref + (below[0] if len(below) else len(t) - i_ref) - 1
46 if last <= i_ref + 5:
47 continue
48 r = (Er - Es)[b, i_ref:last] - G0[b]
49 print(f" {C[b]:6.0f} Hz {t[last]:5.2f} s {np.sqrt(np.mean(r ** 2)):5.1f} dB {r[-1]:+6.1f} dB")
51# Static equalizer applied to the simulation, full and smoothed.
52lf = np.log(np.maximum(f, 1.0))
53lc = np.log(C)
54usable = (C > 100) & (C < 6000) # where the recording informs the balance
55coef = np.polyfit(lc[usable], G0[usable], 3)
56G4 = np.polyval(coef, np.clip(lc, lc[usable].min(), lc[usable].max()))
57print("\ncubic in log frequency (4 coefficients) versus the 28 band gains, RMS difference over 100 Hz to 6 kHz:",
58 f"{np.sqrt(np.mean((G4 - G0)[usable] ** 2)):.1f} dB")
59i_a, i_b = int(0.1 * FS), int(1.0 * FS)
60body = lambda x: np.sqrt(np.mean(x[i_a:i_b] ** 2))
61target = body(ref)
62for name, G in (("eq", G0), ("eq4", G4)):
63 Gf = np.interp(lf, lc, G)
64 _, y = istft(Zs * 10 ** (Gf / 20)[:, None], FS, nperseg=NPERSEG, noverlap=NPERSEG - HOP)
65 y = y[:n]
66 M.save(f"{out}_{name}.wav", np.tanh(y * (target / body(y))))
67print(f"wrote {out}_eq.wav and {out}_eq4.wav")