model: диагностический мост FUN_180529fe0 (rmse 0.236 dB): C(f)=1.221*LUT(log10(L0/res)) + 0.358*warp(f)^3.143 — dual2000-глубина через АДДИТИВНЫЙ аккумулятор 0x5407c8 (не warp*LUT, тот >10 dB), res(500)=0.117 Q-независим (C500-const), warp^3.14~pi => кратный каскад

This commit is contained in:
2026-08-18 18:02:31 +03:00
parent 111d47530e
commit b1066f31e0
+128
View File
@@ -0,0 +1,128 @@
#!/usr/bin/env python3
"""model_fir.py — РЕЗУЛЬТАТ диагностического моста (2026-08-18).
Вопрос: даёт ли РЕАЛЬНАЯ цепочка FUN_180529fe0 (masкa -> warp -> FIR/DC) глубину
dual_b1q с реDD-константами из декомпа? Ответ: ДА, через АДДИТИВНЫЙ аккумулятор.
МОДЕЛЬ (rmse=0.236 dB на 36 точках):
C(f) = g·LUT(log10(L0/res_band(f))) + w·warp(f)^a
red = -20·log10(1-C)
g=1.221 w=0.358 a=3.143 (фит, bounds [0.9..1.5],[0.2..1.2],[2..6])
LUT = PCHIP B.11 (заморожен), warp = 0.87·Kx/(K+x), K=exp(2.0723), x=f/2000
КЛЮЧЕВЫЕ ВЫВОДЫ:
1. dual500-константа объяснена БЕЗ tilt: res_band(500; fc=500)=0.117 не зависит от Q
(твин в центре полосы Q-независим); C500 = g·LUT(0.574)·warp^a(500≈0) = 0.692 -> 10.23 dB.
2. Глубина dual2000 НЕ мультипликативный warp·LUT (ratio warp 2000/500=3.66 -> rmse>10 dB);
её даёт добавка w·warp(f)^3.14: на 2000 Гц вклад 0.358·0.773^3.14=0.149, на 1000<=0.02
(иначе t1k/t1kq рушатся) => точка фита фиксирована единственной. Экзотика a≈π: кандидат —
кратный каскад warp (freq-axis 0x540698 / ∏0x540688 / двойной FFT-проход).
3. Остаток dual2000 при Q=0.1 (14.4 vs 15.2) = тонкая форма LUT-колена 0.574 (0.7 dB) —
не закрывается per-bin формами; требует точного ур-я импульса (FFT-уровень 0x535a70).
4. Порядок измерений в фите: DUAL.ravel() (500,2000 попарно) | T1KQ | T1K.
5. res(2000) растёт 0.216(Q=0.1)..1.988(Q=10) и входит ЧЕРЕЗ LUT(xv), где xv=log10(L0/res):
xv(2000)=0.308(Q=0.1)..-0.656(Q=10) -> C2000-Q-зависимость совпадает.
"""
import numpy as np
from scipy.optimize import least_squares
from scipy.interpolate import PchipInterpolator
FS = 44100.0
GAIN = 4.132
QS = [0.1, 0.2, 0.3, 0.5, 0.7, 1.0, 1.5, 2.0, 3.0, 5.0, 10.0]
DUAL = np.array([(10.220, 15.224), (10.219, 13.963), (10.219, 12.947),
(10.219, 11.725), (10.218, 11.119), (10.216, 10.689),
(10.210, 10.412), (10.203, 10.305), (10.182, 10.225),
(10.113, 10.183), (9.822, 10.165)])
FCS = [800.0, 900.0, 950.0, 1000.0, 1050.0, 1100.0, 1200.0]
T1KQ = np.array([7.868, 8.536, 8.726, 8.788, 8.729, 8.575, 8.115])
T1K = np.array([14.548, 15.332, 15.553, 15.626, 15.557, 15.378, 14.840])
L0_DUAL = 10 ** (-7.142 / 20)
L0_T1KQ = 10 ** (-18.063 / 20)
L0_T1K = 1.0
LUT_X = np.array([-0.750, -0.500, -0.250, 0.000, 0.250, 0.500, 0.574, 0.610, 0.750, 1.000])
LUT_Y = np.array([0.4402, 0.4552, 0.4813, 0.5072, 0.5329, 0.5332, 0.5645, 0.6471, 0.6562, 0.6670])
_LUT = PchipInterpolator(LUT_X, LUT_Y)
def lut(x):
return np.float64(np.clip(_LUT(float(x)), LUT_Y[0], LUT_Y[-1]))
def res_band(f, fc, Q):
w0 = fc * 2 * np.pi / FS
c, s = np.cos(w0), np.sin(w0)
p = (s * 0.5) / Q
a, a2 = p * GAIN, p / GAIN
A = [a + 1, -2 * c, 1 - a]
B = [a2 + 1, -2 * c, 1 - a2]
w = 2 * np.pi * f / FS
z = np.exp(-1j * w)
return np.abs(2.0 * (B[0] + B[1] * z + B[2] * z * z) /
(A[0] + A[1] * z + A[2] * z * z))
def warp(f):
x = f / 2000.0
return 0.87 * 7.942 * x / (7.942 + x)
def make_combos():
combos = []
for q in QS:
combos += [(500.0, 500.0, q, L0_DUAL), (2000.0, 500.0, q, L0_DUAL)]
for fc in FCS:
combos.append((1000.0, fc, 0.9999978, L0_T1KQ))
for fc in FCS:
combos.append((1000.0, fc, 0.9999978, L0_T1K))
return combos
COMBOS = make_combos()
def pred(P):
g, w, a = P
out = []
for ft, fc, Q, L0 in COMBOS:
xv = np.log10(L0 / res_band(ft, fc, Q))
m = g * lut(xv) + w * warp(ft) ** a
out.append(-20 * np.log10(1 - np.clip(m, 0, 0.999)))
return np.array(out)
def run():
meas = np.concatenate([DUAL.ravel(), T1KQ, T1K])
rng = np.random.default_rng(2)
best = None
for trial in range(40):
try:
p0 = rng.uniform([0.9, 0.2, 2.0], [1.5, 1.2, 6.0])
r = least_squares(lambda p: pred(p) - meas, p0,
bounds=([0.9, 0.2, 2.0], [1.5, 1.2, 6.0]),
max_nfev=400, ftol=1e-11, x_scale='jac')
pr = pred(r.x)
rm = np.sqrt(np.mean((pr - meas) ** 2))
if best is None or rm < best[0]:
best = (rm, pr, r.x)
except Exception:
pass
rm, pr, P = best
seg = lambda i0, i1: np.sqrt(np.mean((pr[i0:i1] - meas[i0:i1]) ** 2))
print(f'rmse={rm:.4f} dual={seg(0,22):.3f} t1kq={seg(22,29):.3f} t1k={seg(29,36):.3f} '
f'P(g,w,a)={np.round(P,4)} [warp^a: 2000={warp(2000)**P[2]:.3f}, 1000={warp(1000)**P[2]:.3f}]')
print('d500=' + ' '.join(f'{pr[2*i]:5.2f}/{meas[2*i]:5.2f}' for i in range(11)))
print('d2k =' + ' '.join(f'{pr[2*i+1]:5.2f}/{meas[2*i+1]:5.2f}' for i in range(11)))
print('t1kq=' + ' '.join(f'{pr[22+j]:5.2f}/{meas[22+j]:5.2f}' for j in range(7)))
print('t1k =' + ' '.join(f'{pr[29+j]:5.2f}/{meas[29+j]:5.2f}' for j in range(7)))
# числовая справка res/warp
print('\nres(500;fc500)=%.3f (Q-независим, объясняет C500-константу)' % res_band(500, 500, 1.0))
print('res(2000): ' + ' '.join('%.3f' % res_band(2000, 500, q) for q in QS))
print('xv(2000)=log10(L0D/res): ' +
' '.join('%.3f' % np.log10(L0_DUAL / res_band(2000, 500, q)) for q in QS))
if __name__ == '__main__':
run()