56 lines
2.3 KiB
Python
56 lines
2.3 KiB
Python
#!/usr/bin/env python3
|
|
"""Phase B step 1c: neighborhood MEAN/RMS pooling scan (max already refuted).
|
|
|
|
Rationale: slide_max raises lvl at tone bins via sidelobes -> over-reduction
|
|
(tones broke -2.7..-4.9). Mean/RMS pooling does the opposite for an isolated
|
|
narrow peak among quiet neighbours -> less reduction on tones, ~neutral on
|
|
wide noise. Pooling in LINEAR lvl domain (am ~ lvl near flat res).
|
|
"""
|
|
import numpy as np
|
|
|
|
exec(open('/tmp/opencode/phaseA_grid_fast.py').read().split("# old-law reference gains")[0])
|
|
|
|
def slide(x, w, op):
|
|
if w <= 1: return x
|
|
h = w // 2
|
|
xp = np.pad(x, ((0, 0), (h, h)), mode='edge')
|
|
win = np.lib.stride_tricks.sliding_window_view(xp, w, axis=1)
|
|
if op == 'mean': return win.mean(axis=-1)
|
|
return np.sqrt((win ** 2).mean(axis=-1))
|
|
|
|
RAW = {}
|
|
for name, traj in [('res_500','traj_res500.bin'), ('al_12','traj_al12.bin'),
|
|
('al_24','traj_al24.bin'), ('t1k_1000','traj_t1k.bin')]:
|
|
RAW[name] = load_traj('/tmp/opencode/'+traj)[WIN[name][0]:WIN[name][1]]
|
|
|
|
GO = {}
|
|
for name,(dB,W,bm,_) in DATA.items():
|
|
c_old = np.clip((dB+13.78)/82.07, 0, 1)**0.344*4.2
|
|
g = gains_batch(dB, W, bm, np.array([0.]), np.array([0.]), np.array([99.]), C_pre=c_old)
|
|
GO[name] = float(g[0].mean()) if name=='res_500' else float(np.median(g[0]))
|
|
|
|
def eval_variant(op, w):
|
|
saved = {n: DATA[n] for n in DATA}
|
|
for name,(dB,W,bm,_) in DATA.items():
|
|
T = slide(RAW[name], w, op)
|
|
dBp = np.log10(np.maximum(T, 1e-12))*20.0
|
|
if name != 'res_500':
|
|
lv = dBp[:, bm]; keep = dBp[lv >= lv.max()-6]
|
|
else:
|
|
keep = dBp
|
|
DATA[name] = (np.ascontiguousarray(keep), W, bm, None)
|
|
out = {}
|
|
for name,(dB,W,bm,_) in DATA.items():
|
|
g = gains_batch(dB, W, bm, np.array([1.8]), np.array([0.11]), np.array([99.]))
|
|
agg = np.sqrt(np.mean(g**2, axis=1)) if name=='res_500' else np.median(g, axis=1)
|
|
out[name] = MEAS_OLD[name] + 20*np.log10(float(agg[0])/GO[name])
|
|
DATA.update(saved)
|
|
return out
|
|
|
|
print(f'{"op":>4} {"w":>3} ' + ' '.join(f'{n:>9}' for n in DATA) + ' rms')
|
|
for op in ['mean', 'rms']:
|
|
for w in [3, 5, 9, 17, 33]:
|
|
e = eval_variant(op, w)
|
|
tot = np.sqrt(sum(v*v for v in e.values())/len(e))
|
|
print(f'{op:>4} {w:>3} ' + ' '.join(f'{v:+9.2f}' for v in e.values()) + f' {tot:.2f}')
|