72 lines
2.9 KiB
Python
72 lines
2.9 KiB
Python
#!/usr/bin/env python3
|
|
"""Phase B step 3: (rho, Delta) joint scan.
|
|
|
|
Hypothesis: content gap lives in IIR1 spike attenuation vs law level.
|
|
rho = IIR1 pole (DC-normalized: y = rho*acc + (1-rho)*x), canon rho=0.692.
|
|
Delta = additive shift of affine law c = max(1.8+D+0.11*dB, 0).
|
|
Anchored at canon old-law gains GO. Criterion: pred_err ~ 0 on ALL 4 anchors.
|
|
"""
|
|
import numpy as np
|
|
|
|
exec(open('/tmp/opencode/phaseA_grid_fast.py').read().split("# old-law reference gains")[0])
|
|
|
|
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]))
|
|
print('GO:', {k: round(v, 3) for k, v in GO.items()})
|
|
|
|
def gains_rho(dB, W, bm, rho, DLs):
|
|
"""dB [T,N]; law c=max(1.8+Dl+0.11*dB,0); IIR1 pole=rho (DC-norm).
|
|
returns [C,T] gain at bm for each Delta in DLs."""
|
|
C, T, N = len(DLs), dB.shape[0], dB.shape[1]
|
|
c = np.maximum(1.8 + DLs[:, None, None] + 0.11 * dB[None], 0.0)
|
|
acc = np.zeros((C, T)); y = np.empty_like(c)
|
|
for i in range(N):
|
|
acc = rho * acc + (1 - rho) * c[:, :, i]
|
|
y[:, :, i] = acc
|
|
acc = np.zeros((C, T))
|
|
for i in range(N):
|
|
acc = A2[i] * acc + B2[i] * y[:, :, i]
|
|
y[:, :, i] = 0.8 * np.exp2(-acc) * W[i]
|
|
for _ in range(2):
|
|
st = np.zeros((C, T))
|
|
for i in range(N):
|
|
st = y[:, :, i] * B3[i] + st * A3[i]
|
|
y[:, :, i] = st
|
|
st = y[:, :, -1].copy()
|
|
for i in range(N - 2, 0, -1):
|
|
st = y[:, :, i] * B3[i] + st * A3[i]
|
|
y[:, :, i] = st
|
|
return y[:, :, bm]
|
|
|
|
RHOS = np.linspace(0.30, 0.95, 131)
|
|
DLS = np.linspace(-1.5, 1.5, 121)
|
|
names = list(DATA)
|
|
best = []
|
|
for rho in RHOS:
|
|
ev = {}
|
|
for name, (dB, W, bm, _) in DATA.items():
|
|
g = gains_rho(dB, W, bm, rho, DLS)
|
|
agg = np.sqrt(np.mean(g ** 2, axis=1)) if name == 'res_500' else np.median(g, axis=1)
|
|
ev[name] = MEAS_OLD[name] + 20 * np.log10(agg / GO[name])
|
|
E = np.stack([ev[n] for n in names]) # [4, C]
|
|
rms = np.sqrt((E ** 2).mean(axis=0)) # per Delta
|
|
j = int(rms.argmin())
|
|
best.append((rms[j], rho, DLS[j], E[:, j]))
|
|
best.sort()
|
|
print('\ntop-10 (rms over 4 anchors):')
|
|
for rms, rho, dl, e in best[:10]:
|
|
print(f' rho={rho:.3f} D={dl:+.3f} rms={rms:.3f} ' +
|
|
' '.join(f'{n[:5]}:{v:+.2f}' for n, v in zip(names, e)))
|
|
print(f'\ncanon rho=0.692 D=0 reference:')
|
|
j0 = int(np.argmin(np.abs(DLS)))
|
|
for rho in [0.692]:
|
|
ev = {}
|
|
for name, (dB, W, bm, _) in DATA.items():
|
|
g = gains_rho(dB, W, bm, rho, DLS[j0:j0+1])
|
|
agg = np.sqrt(np.mean(g[0] ** 2)) if name == 'res_500' else np.median(g[0])
|
|
ev[name] = MEAS_OLD[name] + 20 * np.log10(agg / GO[name])
|
|
print(' ' + ' '.join(f'{n[:5]}:{v:+.2f}' for n, v in ev.items()))
|