From 32a63955f8b96fed306c9c65822f88d41b385971 Mon Sep 17 00:00:00 2001 From: Matiq Date: Tue, 18 Aug 2026 13:02:56 +0300 Subject: [PATCH] =?UTF-8?q?B.11:=20LUT-=D0=BA=D1=80=D0=B8=D0=B2=D0=B0?= =?UTF-8?q?=D1=8F=20=D1=83=D1=80=D0=BE=D0=B2=D0=B5=D0=BD=D1=8C->=D0=BC?= =?UTF-8?q?=D0=B0=D1=81=D0=BA=D0=B0=20=D0=B8=D0=B7=D0=B2=D0=BB=D0=B5=D1=87?= =?UTF-8?q?=D0=B5=D0=BD=D0=B0=20(rmse=200.072=20dB,=2036=20=D1=82=D0=BE?= =?UTF-8?q?=D1=87=D0=B5=D0=BA,=203=20=D1=83=D1=80=D0=BE=D0=B2=D0=BD=D1=8F)?= =?UTF-8?q?:=20=D0=BA=D0=BE=D0=BB=D0=BB=D0=B0=D0=BF=D1=81=20dual+t1kq+t1k?= =?UTF-8?q?=20=D0=BD=D0=B0=20=D0=9E=D0=94=D0=9D=D0=A3=20=D0=BA=D1=80=D0=B8?= =?UTF-8?q?=D0=B2=D1=83=D1=8E=20x=3Dlog10(L0/res);=20=D1=81=D1=82=D0=B5?= =?UTF-8?q?=D0=BF=D0=B5=D0=BD=D0=BD=D0=BE=D0=B9=20=D0=B7=D0=B0=D0=BA=D0=BE?= =?UTF-8?q?=D0=BD=20C~L0^p=20=D0=9D=D0=95=20=D1=80=D0=B0=D0=B1=D0=BE=D1=82?= =?UTF-8?q?=D0=B0=D0=B5=D1=82=20=D0=BD=D0=B0=200dB=20(p=5Feff=200.19=20vs?= =?UTF-8?q?=200.085)=20-=20LUT=20=D0=BD=D0=B0=D1=81=D1=8B=D1=89=D0=B0?= =?UTF-8?q?=D1=8E=D1=89=D0=B0=D1=8F=20gamma=20=D0=B8=D0=B7=20FUN=5F1805634?= =?UTF-8?q?40=20(param=5F1+0x188);=20=D1=84=D0=BE=D1=80=D0=BC=D0=B0=20?= =?UTF-8?q?=D1=80=D0=B5=D0=B7=D0=BE=D0=BD=D0=B0=D0=BD=D1=81=D0=B0=20=D1=81?= =?UTF-8?q?=D1=82=D0=B0=D0=B1=D0=B8=D0=BB=D1=8C=D0=BD=D0=B0;=20model=5Flut?= =?UTF-8?q?.py=20=D0=BA=D0=B0=D0=BD=D0=BE=D0=BD=D0=B8=D1=87=D0=B5=D1=81?= =?UTF-8?q?=D0=BA=D0=B8=D0=B9=20+=20roadmap=20B.11?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- extract_lut.py | 73 +++++++++++++++++++++++++++++++++ fit_level.py | 86 +++++++++++++++++++++++++++++++++++++++ fit_lut.py | 87 +++++++++++++++++++++++++++++++++++++++ fit_lut2.py | 82 +++++++++++++++++++++++++++++++++++++ fit_lut3.py | 84 ++++++++++++++++++++++++++++++++++++++ fit_lut4.py | 80 ++++++++++++++++++++++++++++++++++++ fit_lut5.py | 84 ++++++++++++++++++++++++++++++++++++++ model_lut.py | 108 +++++++++++++++++++++++++++++++++++++++++++++++++ roadmap.md | 32 +++++++++++++++ verify_lut.py | 72 +++++++++++++++++++++++++++++++++ verify_lut2.py | 76 ++++++++++++++++++++++++++++++++++ verify_lut3.py | 66 ++++++++++++++++++++++++++++++ verify_lut4.py | 91 +++++++++++++++++++++++++++++++++++++++++ 13 files changed, 1021 insertions(+) create mode 100644 extract_lut.py create mode 100644 fit_level.py create mode 100644 fit_lut.py create mode 100644 fit_lut2.py create mode 100644 fit_lut3.py create mode 100644 fit_lut4.py create mode 100644 fit_lut5.py create mode 100644 model_lut.py create mode 100644 verify_lut.py create mode 100644 verify_lut2.py create mode 100644 verify_lut3.py create mode 100644 verify_lut4.py diff --git a/extract_lut.py b/extract_lut.py new file mode 100644 index 0000000..d311cc8 --- /dev/null +++ b/extract_lut.py @@ -0,0 +1,73 @@ +#!/usr/bin/env python3 +"""extract_lut.py — извлечение LUT-кривой level->mask из коллапса. + +Форма B.10 ЗАФИКСИРОВАНА (Q=0.9, gain=4.13, tilt 1.414/1.454/1.795). +По каждой из 36 точек вычисляем x=log10(L0/res) и y_lut=(1-10^(-red/20))/(depth*tilt). +Все точки должны лечь на ОДНУ монотонную кривую y_lut(x). Печатаем кривую. +""" +import numpy as np + +FS = 44100.0 +DEPTH = 0.8639736175537109 +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]) +L_DUAL = 10 ** (-7.142 / 20) +L_T1KQ = 10 ** (-18.063 / 20) +L_T1K = 1.0 +TILT = {500: 1.414, 1000: 1.454, 2000: 1.795} +Q, G = 0.900, 4.132 + + +def res_at(ft, fc, Q, g): + w0 = fc * 2 * np.pi / FS + c, s = np.cos(w0), np.sin(w0) + p = (s * 0.5) / Q + a, a2 = p * g, p / g + A = [a + 1, -2 * c, 1 - a] + B = [a2 + 1, -2 * c, 1 - a2] + w = 2 * np.pi * ft / 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 add(pts, L, res, red, tilt, tag): + y = (1 - 10 ** (-red / 20)) / (DEPTH * tilt) + pts.append((np.log10(L / res), y, tag)) + + +pts = [] +for q in QS: + for f in (500.0, 2000.0): + add(pts, L_DUAL, res_at(f, 500, q, G), DUAL[QS.index(q)][0 if f < 1000 else 1], + TILT[f], 'dual') +for i, fc in enumerate(FCS): + r = res_at(1000, fc, 0.9999978, G) + add(pts, L_T1KQ, r, T1KQ[i], TILT[1000], 't1kq') + add(pts, L_T1K, r, T1K[i], TILT[1000], 't1k') +pts.sort() +# печать всех точек: x, y_lut +print('x=log10(L0/res) y_lut tag') +for x, y, tag in pts: + print(f'{x:+.3f} {y:.4f} {tag}') +# биннинг для кривой +import collections +bins = collections.defaultdict(list) +for x, y, tag in pts: + bins[round(x * 4) / 4].append(y) +print('\nкривая (бин 0.25 по x):') +xs, ys = [], [] +for bx in sorted(bins): + v = np.mean(bins[bx]) + xs.append(bx); ys.append(v) + print(f'x={bx:+.2f} y={v:.4f} (n={len(bins[bx])}, spread={np.std(bins[bx]):.4f})') +# подгонка свободной кривой к (xs, ys): монотонная интерполяция +print('\nНЕПАРАМЕТРИЧЕСКАЯ КРИВАЯ (для model):') +for bx, v in zip(xs, ys): + print(f' ({bx:+.3f}, {v:.4f}),') diff --git a/fit_level.py b/fit_level.py new file mode 100644 index 0000000..2848007 --- /dev/null +++ b/fit_level.py @@ -0,0 +1,86 @@ +#!/usr/bin/env python3 +"""fit_level.py — расширенная B.10: три уровня входа (0, -7, -18 dBFS). + +Проверяем, что резонансная форма |2B/A| одна для всех уровней, а глубина +масштабируется законом C(dB). Вопрос: степенной ли закон C ~ L0^p или +экспоненциальный по dB C ~ 10^(dB/k) с k=dB-независимым? +""" +import numpy as np +from scipy.optimize import least_squares + +FS = 44100.0 +DEPTH = 0.8639736175537109 + +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] +# измеренные (компонент-корреляция 1к, надёжно): +T1KQ = np.array([7.868, 8.536, 8.726, 8.788, 8.729, 8.575, 8.115]) # -18.06 dB +T1K = np.array([14.548, 15.332, 15.553, 15.626, 15.557, 15.378, 14.840]) # 0 dB +L_DUAL = 10 ** (-7.142 / 20) +L_T1KQ = 10 ** (-18.063 / 20) +L_T1K = 1.0 + + +def res_at(f_tone, fc, Q, gain): + 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_tone / 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 model(p): + Q, gain, p0, p1, D0, t1000 = p # p(dB) = p0 + p1*dB + out = [] + for q in QS: + for f in (500.0, 2000.0): + r = res_at(f, 500, q, gain) + tilt = t1000 if f == 2000 else 1.45 + c = DEPTH * tilt * D0 * (L_DUAL / r) ** (p0 + p1 * (-7.142)) + out.append(-20 * np.log10(max(1 - c, 1e-9))) + for i, fc in enumerate(FCS): + r = res_at(1000, fc, 0.9999978, gain) + c = DEPTH * t1000 * D0 * (L_T1KQ / r) ** (p0 + p1 * (-18.063)) + out.append(-20 * np.log10(max(1 - c, 1e-9))) + c = DEPTH * t1000 * D0 * (L_T1K / r) ** p0 + out.append(-20 * np.log10(max(1 - c, 1e-9))) + return np.array(out) + + +def run(): + meas = np.array(list(DUAL.ravel()) + list(T1KQ) + list(T1K)) + x0 = [0.9, 4.24, 0.085, 0.003, 0.5, 1.454] + r = least_squares(lambda p: model(p) - meas, x0, + bounds=([0.1, 0.5, 0.0, -0.01, 0.01, 0.3], + [5, 12, 0.3, 0.02, 5, 5]), + max_nfev=10000, xtol=1e-12, ftol=1e-12) + Q, gain, p0, p1, D0, t1000 = r.x + rmse = np.sqrt(np.mean((model(r.x) - meas) ** 2)) + print(f'FIT rmse={rmse:.4f} dB Q={Q:.3f} gain={gain:.3f} ' + f'p(dB)={p0:.4f}{p1:+.4f}*dB D0={D0:.3f} t1000={t1000:.3f}') + print('p при 0/-7/-18 dB:', p0, p0+p1*-7.142, p0+p1*-18.063) + pred = model(r.x) + n = 0 + print('--- dual_b1q ---') + for i, q in enumerate(QS): + print(f'q={q:5.1f} {DUAL[i,0]:7.3f}/{pred[2*i]:7.3f} ' + f'{DUAL[i,1]:7.3f}/{pred[2*i+1]:7.3f}') + print('--- t1kq -18dB ---') + for i, fc in enumerate(FCS): + print(f'fc={fc:5.0f} {T1KQ[i]:6.3f}/{pred[22+i]:6.3f}') + print('--- t1k 0dB ---') + for i, fc in enumerate(FCS): + print(f'fc={fc:5.0f} {T1K[i]:6.3f}/{pred[29+i]:6.3f}') + + +if __name__ == '__main__': + run() diff --git a/fit_lut.py b/fit_lut.py new file mode 100644 index 0000000..83df9ff --- /dev/null +++ b/fit_lut.py @@ -0,0 +1,87 @@ +#!/usr/bin/env python3 +"""fit_lut.py — полная LUT «уровень->маска» (B.11). + +Все 3 набора (0/-7/-18 dBFS) ложатся на ОДНУ кривую C_norm = f(L0/res), +f = min+(max-min)*(k*z)^(1/c) — gamma-LUT из декомпиляции (FUN_180563440, +кривая param_1+0x188, линейный флаг: val=min+(max-min)*x^(1/c)). +""" +import numpy as np +from scipy.optimize import least_squares + +FS = 44100.0 +DEPTH = 0.8639736175537109 +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]) +L_DUAL = 10 ** (-7.142 / 20) +L_T1KQ = 10 ** (-18.063 / 20) +L_T1K = 1.0 +TILT = {500: 1.414, 1000: 1.454, 2000: 1.795} + + +def res_at(ft, fc, Q, g): + w0 = fc * 2 * np.pi / FS + c, s = np.cos(w0), np.sin(w0) + p = (s * 0.5) / Q + a, a2 = p * g, p / g + A = [a + 1, -2 * c, 1 - a] + B = [a2 + 1, -2 * c, 1 - a2] + w = 2 * np.pi * ft / 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 lut(z, mn, mx, c, k): + v = mn + (mx - mn) * (k * z) ** (1.0 / c) + return np.minimum(v, mx) + + +def model(p): + Q, g, mn, mx, c, k = p + out = [] + for q in QS: + for f in (500.0, 2000.0): + r = res_at(f, 500, q, g) + C = DEPTH * TILT[f] * lut(L_DUAL / r, mn, mx, c, k) + out.append(-20 * np.log10(max(1 - C, 1e-9))) + for i, fc in enumerate(FCS): + r = res_at(1000, fc, 0.9999978, g) + C = DEPTH * TILT[1000] * lut(L_T1KQ / r, mn, mx, c, k) + out.append(-20 * np.log10(max(1 - C, 1e-9))) + C = DEPTH * TILT[1000] * lut(L_T1K / r, mn, mx, c, k) + out.append(-20 * np.log10(max(1 - C, 1e-9))) + return np.array(out) + + +def run(): + meas = np.array(list(DUAL.ravel()) + list(T1KQ) + list(T1K)) + x0 = [0.9, 4.13, 0.0, 1.0, 6.0, 1.0] + r = least_squares(lambda p: model(p) - meas, x0, + bounds=([0.1, 0.5, 0.0, 0.5, 1.0, 1e-4], + [5, 12, 2.0, 5.0, 30.0, 50.0]), + max_nfev=20000, xtol=1e-12, ftol=1e-12) + Q, g, mn, mx, c, k = r.x + rmse = np.sqrt(np.mean((model(r.x) - meas) ** 2)) + print(f'LUT-FIT rmse={rmse:.4f} dB Q={Q:.3f} gain={g:.3f}') + print(f'LUT: min={mn:.3f} max={mx:.3f} gamma={c:.3f} k={k:.4f}') + pred = model(r.x) + print('--- dual_b1q ---') + for i, q in enumerate(QS): + print(f'q={q:5.1f} {DUAL[i,0]:7.3f}/{pred[2*i]:7.3f} ' + f'{DUAL[i,1]:7.3f}/{pred[2*i+1]:7.3f}') + print('--- t1kq -18dB ---') + for i, fc in enumerate(FCS): + print(f'fc={fc:5.0f} {T1KQ[i]:6.3f}/{pred[22+i]:6.3f}') + print('--- t1k 0dB ---') + for i, fc in enumerate(FCS): + print(f'fc={fc:5.0f} {T1K[i]:6.3f}/{pred[29+i]:6.3f}') + + +if __name__ == '__main__': + run() diff --git a/fit_lut2.py b/fit_lut2.py new file mode 100644 index 0000000..bc023a1 --- /dev/null +++ b/fit_lut2.py @@ -0,0 +1,82 @@ +#!/usr/bin/env python3 +"""fit_lut2.py — логистическая LUT (B.11): C_norm = mn + (mx-mn)/(1+exp(-(x-x0)/w)), +x = log10(L0/res). Коллапс всех 36 точек на одну кривую.""" +import numpy as np +from scipy.optimize import least_squares + +FS = 44100.0 +DEPTH = 0.8639736175537109 +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]) +L_DUAL = 10 ** (-7.142 / 20) +L_T1KQ = 10 ** (-18.063 / 20) +L_T1K = 1.0 +TILT = {500: 1.414, 1000: 1.454, 2000: 1.795} + + +def res_at(ft, fc, Q, g): + w0 = fc * 2 * np.pi / FS + c, s = np.cos(w0), np.sin(w0) + p = (s * 0.5) / Q + a, a2 = p * g, p / g + A = [a + 1, -2 * c, 1 - a] + B = [a2 + 1, -2 * c, 1 - a2] + w = 2 * np.pi * ft / 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 lut(x, mn, mx, x0, w): + return mn + (mx - mn) / (1.0 + np.exp(-(x - x0) / w)) + + +def model(p): + Q, g, mn, mx, x0, w = p + out = [] + for q in QS: + for f in (500.0, 2000.0): + r = res_at(f, 500, q, g) + C = DEPTH * TILT[f] * lut(np.log10(L_DUAL / r), mn, mx, x0, w) + out.append(-20 * np.log10(max(1 - C, 1e-9))) + for i, fc in enumerate(FCS): + r = res_at(1000, fc, 0.9999978, g) + C = DEPTH * TILT[1000] * lut(np.log10(L_T1KQ / r), mn, mx, x0, w) + out.append(-20 * np.log10(max(1 - C, 1e-9))) + C = DEPTH * TILT[1000] * lut(np.log10(L_T1K / r), mn, mx, x0, w) + out.append(-20 * np.log10(max(1 - C, 1e-9))) + return np.array(out) + + +def run(): + meas = np.array(list(DUAL.ravel()) + list(T1KQ) + list(T1K)) + x0 = [0.9, 4.13, 0.0, 1.0, 0.5, 0.4] + r = least_squares(lambda p: model(p) - meas, x0, + bounds=([0.1, 0.5, 0.0, 0.4, -1.0, 0.01], + [5, 12, 0.5, 5.0, 3.0, 5.0]), + max_nfev=30000, xtol=1e-12, ftol=1e-12) + Q, g, mn, mx, x0, w = r.x + rmse = np.sqrt(np.mean((model(r.x) - meas) ** 2)) + print(f'LOGISTIC LUT-FIT rmse={rmse:.4f} dB Q={Q:.3f} gain={g:.3f}') + print(f'LUT: mn={mn:.3f} mx={mx:.3f} x0={x0:.3f} w={w:.3f}') + pred = model(r.x) + print('--- dual_b1q ---') + for i, q in enumerate(QS): + print(f'q={q:5.1f} {DUAL[i,0]:7.3f}/{pred[2*i]:7.3f} ' + f'{DUAL[i,1]:7.3f}/{pred[2*i+1]:7.3f}') + print('--- t1kq -18dB ---') + for i, fc in enumerate(FCS): + print(f'fc={fc:5.0f} {T1KQ[i]:6.3f}/{pred[22+i]:6.3f}') + print('--- t1k 0dB ---') + for i, fc in enumerate(FCS): + print(f'fc={fc:5.0f} {T1K[i]:6.3f}/{pred[29+i]:6.3f}') + + +if __name__ == '__main__': + run() diff --git a/fit_lut3.py b/fit_lut3.py new file mode 100644 index 0000000..d7a4d6a --- /dev/null +++ b/fit_lut3.py @@ -0,0 +1,84 @@ +#!/usr/bin/env python3 +"""fit_lut3.py — LUT-кривая как свободная функция (B.11). + +Форма резонанса (Q,gain,tilt) ЗАФИКСИРОВАНА по B.10. Фитится только +LUT g(x), x=log10(L0/res): крас=-20log10(1 - DEPTH*tilt*g(x)). +""" +import numpy as np +from scipy.optimize import least_squares + +FS = 44100.0 +DEPTH = 0.8639736175537109 +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]) +L_DUAL = 10 ** (-7.142 / 20) +L_T1KQ = 10 ** (-18.063 / 20) +L_T1K = 1.0 +TILT = {500: 1.414, 1000: 1.454, 2000: 1.795} +QG, GG = 0.900, 4.132 + + +def res_at(ft, fc, Q, g): + w0 = fc * 2 * np.pi / FS + c, s = np.cos(w0), np.sin(w0) + p = (s * 0.5) / Q + a, a2 = p * g, p / g + A = [a + 1, -2 * c, 1 - a] + B = [a2 + 1, -2 * c, 1 - a2] + w = 2 * np.pi * ft / 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 lut(x, mn, mx, x0, w): + return mn + (mx - mn) / (1.0 + np.exp(-(x - x0) / w)) + + +def model(p): + mn, mx, x0, w = p + out = [] + for q in QS: + for f in (500.0, 2000.0): + r = res_at(f, 500, q, GG) + C = DEPTH * TILT[f] * lut(np.log10(L_DUAL / r), mn, mx, x0, w) + out.append(-20 * np.log10(max(1 - C, 1e-9))) + for i, fc in enumerate(FCS): + r = res_at(1000, fc, 0.9999978, GG) + C = DEPTH * TILT[1000] * lut(np.log10(L_T1KQ / r), mn, mx, x0, w) + out.append(-20 * np.log10(max(1 - C, 1e-9))) + C = DEPTH * TILT[1000] * lut(np.log10(L_T1K / r), mn, mx, x0, w) + out.append(-20 * np.log10(max(1 - C, 1e-9))) + return np.array(out) + + +def run(): + meas = np.array(list(DUAL.ravel()) + list(T1KQ) + list(T1K)) + r = least_squares(lambda p: model(p) - meas, [0.44, 0.67, -0.1, 0.3], + bounds=([0.3, 0.5, -0.8, 0.02], [0.6, 1.0, 0.5, 2.0]), + max_nfev=30000, xtol=1e-12, ftol=1e-12) + mn, mx, x0, w = r.x + rmse = np.sqrt(np.mean((model(r.x) - meas) ** 2)) + print(f'LUT-FIT (shape fixed) rmse={rmse:.4f} dB') + print(f'LUT logistic: mn={mn:.3f} mx={mx:.3f} x0={x0:.3f} w={w:.3f}') + pred = model(r.x) + print('--- dual_b1q ---') + for i, q in enumerate(QS): + print(f'q={q:5.1f} {DUAL[i,0]:7.3f}/{pred[2*i]:7.3f} ' + f'{DUAL[i,1]:7.3f}/{pred[2*i+1]:7.3f}') + print('--- t1kq -18dB ---') + for i, fc in enumerate(FCS): + print(f'fc={fc:5.0f} {T1KQ[i]:6.3f}/{pred[22+i]:6.3f}') + print('--- t1k 0dB ---') + for i, fc in enumerate(FCS): + print(f'fc={fc:5.0f} {T1K[i]:6.3f}/{pred[29+i]:6.3f}') + + +if __name__ == '__main__': + run() diff --git a/fit_lut4.py b/fit_lut4.py new file mode 100644 index 0000000..1197a99 --- /dev/null +++ b/fit_lut4.py @@ -0,0 +1,80 @@ +#!/usr/bin/env python3 +"""fit_lut4.py — LUT-кривая логистическая (B.11), форма B.10 зафиксирована.""" +import numpy as np +from scipy.optimize import least_squares + +FS = 44100.0 +DEPTH = 0.8639736175537109 +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]) +L_DUAL = 10 ** (-7.142 / 20) +L_T1KQ = 10 ** (-18.063 / 20) +L_T1K = 1.0 +TILT = {500: 1.414, 1000: 1.454, 2000: 1.795} +GG = 4.132 + + +def res_at(ft, fc, Q, g): + w0 = fc * 2 * np.pi / FS + c, s = np.cos(w0), np.sin(w0) + p = (s * 0.5) / Q + a, a2 = p * g, p / g + A = [a + 1, -2 * c, 1 - a] + B = [a2 + 1, -2 * c, 1 - a2] + w = 2 * np.pi * ft / 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 lut(x, mn, mx, x0, w): + return mn + (mx - mn) / (1.0 + np.exp(-(x - x0) / w)) + + +def model(p): + mn, mx, x0, w = p + out = [] + for q in QS: + for f in (500.0, 2000.0): + r = res_at(f, 500, q, GG) + C = DEPTH * TILT[f] * lut(np.log10(L_DUAL / r), mn, mx, x0, w) + out.append(-20 * np.log10(max(1 - C, 1e-9))) + for i, fc in enumerate(FCS): + r = res_at(1000, fc, 0.9999978, GG) + C = DEPTH * TILT[1000] * lut(np.log10(L_T1KQ / r), mn, mx, x0, w) + out.append(-20 * np.log10(max(1 - C, 1e-9))) + C = DEPTH * TILT[1000] * lut(np.log10(L_T1K / r), mn, mx, x0, w) + out.append(-20 * np.log10(max(1 - C, 1e-9))) + return np.array(out) + + +def run(): + meas = np.array(list(DUAL.ravel()) + list(T1KQ) + list(T1K)) + r = least_squares(lambda p: model(p) - meas, [0.44, 0.67, -0.1, 0.3], + bounds=([0.3, 0.5, -0.8, 0.02], [0.6, 1.0, 0.5, 2.0]), + max_nfev=30000, xtol=1e-12, ftol=1e-12) + mn, mx, x0, w = r.x + rmse = np.sqrt(np.mean((model(r.x) - meas) ** 2)) + print(f'LUT-FIT (shape fixed) rmse={rmse:.4f} dB') + print(f'LUT logistic: mn={mn:.3f} mx={mx:.3f} x0={x0:.3f} w={w:.3f}') + pred = model(r.x) + print('--- dual_b1q ---') + for i, q in enumerate(QS): + print(f'q={q:5.1f} {DUAL[i,0]:7.3f}/{pred[2*i]:7.3f} ' + f'{DUAL[i,1]:7.3f}/{pred[2*i+1]:7.3f}') + print('--- t1kq -18dB (idx 22+2i) ---') + for i, fc in enumerate(FCS): + print(f'fc={fc:5.0f} {T1KQ[i]:6.3f}/{pred[22+2*i]:6.3f}') + print('--- t1k 0dB (idx 23+2i) ---') + for i, fc in enumerate(FCS): + print(f'fc={fc:5.0f} {T1K[i]:6.3f}/{pred[23+2*i]:6.3f}') + + +if __name__ == '__main__': + run() diff --git a/fit_lut5.py b/fit_lut5.py new file mode 100644 index 0000000..790cbdd --- /dev/null +++ b/fit_lut5.py @@ -0,0 +1,84 @@ +#!/usr/bin/env python3 +"""fit_lut5.py — gamma-LUT из декомпиляции (FUN_180563440): +val = min + (max-min)*(norm_level)^(1/gamma), norm_level = (L0/res)/k.""" +import numpy as np +from scipy.optimize import least_squares + +FS = 44100.0 +DEPTH = 0.8639736175537109 +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]) +L_DUAL = 10 ** (-7.142 / 20) +L_T1KQ = 10 ** (-18.063 / 20) +L_T1K = 1.0 + + +def res_at(ft, fc, Q, g): + w0 = fc * 2 * np.pi / FS + c, s = np.cos(w0), np.sin(w0) + p = (s * 0.5) / Q + a, a2 = p * g, p / g + A = [a + 1, -2 * c, 1 - a] + B = [a2 + 1, -2 * c, 1 - a2] + w = 2 * np.pi * ft / 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 lut(z, mn, mx, gamma, k): + v = mn + (mx - mn) * (np.minimum(z * k, 1.0)) ** (1.0 / gamma) + return np.minimum(v, mx) + + +def model(p): + Q, g, t500, t1000, t2000, mn, mx, gamma, k = p + out = [] + for q in QS: + for f in (500.0, 2000.0): + r = res_at(f, 500, q, g) + tilt = t500 if f < 1000 else t2000 + C = DEPTH * tilt * lut(L_DUAL / r, mn, mx, gamma, k) + out.append(-20 * np.log10(max(1 - C, 1e-9))) + for i, fc in enumerate(FCS): + r = res_at(1000, fc, 0.9999978, g) + C = DEPTH * t1000 * lut(L_T1KQ / r, mn, mx, gamma, k) + out.append(-20 * np.log10(max(1 - C, 1e-9))) + C = DEPTH * t1000 * lut(L_T1K / r, mn, mx, gamma, k) + out.append(-20 * np.log10(max(1 - C, 1e-9))) + return np.array(out) + + +def run(): + meas = np.array(list(DUAL.ravel()) + list(T1KQ) + list(T1K)) + x0 = [0.9, 4.13, 1.4, 1.45, 1.8, 0.44, 0.67, 2.0, 1.0] + r = least_squares(lambda p: model(p) - meas, x0, + bounds=([0.1, 0.5, 0.3, 0.3, 0.3, 0.3, 0.5, 0.5, 1e-3], + [5, 12, 5, 5, 5, 5, 5, 30, 100]), + max_nfev=50000, xtol=1e-13, ftol=1e-13) + Q, g, t500, t1000, t2000, mn, mx, gamma, k = r.x + rmse = np.sqrt(np.mean((model(r.x) - meas) ** 2)) + print(f'GAMMA-LUT FIT rmse={rmse:.4f} dB') + print(f'Q={Q:.3f} gain={g:.3f} tilt:500={t500:.3f} 1000={t1000:.3f} 2000={t2000:.3f}') + print(f'LUT: mn={mn:.4f} mx={mx:.4f} gamma={gamma:.3f} k={k:.4f}') + pred = model(r.x) + print('--- dual_b1q ---') + for i, q in enumerate(QS): + print(f'q={q:5.1f} {DUAL[i,0]:7.3f}/{pred[2*i]:7.3f} ' + f'{DUAL[i,1]:7.3f}/{pred[2*i+1]:7.3f}') + print('--- t1kq -18dB ---') + for i, fc in enumerate(FCS): + print(f'fc={fc:5.0f} {T1KQ[i]:6.3f}/{pred[22+2*i]:6.3f}') + print('--- t1k 0dB ---') + for i, fc in enumerate(FCS): + print(f'fc={fc:5.0f} {T1K[i]:6.3f}/{pred[23+2*i]:6.3f}') + + +if __name__ == '__main__': + run() diff --git a/model_lut.py b/model_lut.py new file mode 100644 index 0000000..0d1c95d --- /dev/null +++ b/model_lut.py @@ -0,0 +1,108 @@ +#!/usr/bin/env python3 +"""model_lut.py — КАНОНИЧЕСКАЯ модель B.11 (2026-08-18). + +Единая непараметрическая LUT-кривая "уровень -> маска" для ВСЕХ уровней входа. + +МОДЕЛЬ (verify, rmse=0.072 dB на 36 точках — dual_b1q 22 + t1kq fc-скан 7 ++ t1k fc-скан 7 при 3 уровнях входа: -7.14 / -18.06 / 0 dBFS): + + red(f) = -20*log10(1 - C(f)) + C(f) = depth * tilt(f) * LUT(log10(L0 / res(f))) + res(f) = |2·B/A|(f; fc, Q, gain) case8/m2c (freq-path, близнец 0x180535880) + LUT = PCHIP-узлы (x=log10(L0/res), y=нормированная маска), таблица ниже + depth = 0.8639736175537109 + L0 = линейный уровень входа (dual 10^(-7.142/20), t1kq 10^(-18.063/20), t1k 1.0) + tilt(f)= 1 - w(f) per-bin level-вес FUN_180530d30 (500/1000/2000: 1.414/1.454/1.795) + +СВОЙСТВА (открытие B.11): + - Степенной закон C = D0·(L0/res)^p (B.10, p=0.0847) НЕ описывает высокий уровень + (0 dBFS): наклон d(ln C)/d(dB) падает с уровнем -> LUT компрессивная (насыщается). + - Кривая монотонна, с "коленом" при x≈0.58 (скачок 0.56 -> 0.65) и асимптотами + y->0.44 (низкий уровень) / y->0.67 (высокий уровень). + - Форма резонанса (Q=0.900, gain=4.132, tilt) ОДИНАКОВА для всех уровней; + различие между dual/t1kq/t1k = только положение x = L0/res на LUT-кривой. + - Подтверждает структурную модель level-path: per-bin уровень (0x540678 IIR-трекеры + 0x540528..) -> LUT-кривая param_1+0x188 (FUN_180563440: mn+(mx-mn)·x^(1/gamma)) + -> маска, применяемая с per-bin весами (0x530d30) в FUN_180529fe0. +""" +import numpy as np +from scipy.interpolate import PchipInterpolator + +FS = 44100.0 +DEPTH = 0.8639736175537109 + +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 + +TILT = {500: 1.414, 1000: 1.454, 2000: 1.795} +Q_FIT, GAIN_FIT = 0.900, 4.132 + +# узлы LUT (x=log10(L0/res), y=норм. маска) — фит по 36 точкам, rmse=0.072 +LUT_KNOTS_X = np.array([-0.750, -0.500, -0.250, 0.000, 0.250, 0.500, 0.574, 0.610, 0.750, 1.000]) +LUT_KNOTS_Y = np.array([0.4402, 0.4552, 0.4813, 0.5072, 0.5329, 0.5332, 0.5645, 0.6471, 0.6562, 0.6670]) + + +def lut(x): + p = PchipInterpolator(LUT_KNOTS_X, LUT_KNOTS_Y) + v = p(np.asarray(x)) + return np.clip(v, LUT_KNOTS_Y[0], LUT_KNOTS_Y[-1]) + + +def res_at(ft, fc, Q, gain): + 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 * ft / 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 red(f_tone, fc, Q, gain, L0, tilt): + r = res_at(f_tone, fc, Q, gain) + C = DEPTH * tilt * lut(np.log10(L0 / r)) + return -20 * np.log10(max(1 - C, 1e-9)) + + +def run(): + preds, meas = [], [] + print('--- dual_b1q (tones 500+2000, band fc=500, -7.142 dBFS) ---') + for i, q in enumerate(QS): + for f, m in ((500.0, DUAL[i, 0]), (2000.0, DUAL[i, 1])): + p = red(f, 500.0, q, GAIN_FIT, L0_DUAL, TILT[f]) + preds.append(p); meas.append(m) + print(f'q={q:5.1f} 500 {DUAL[i,0]:7.3f}/{preds[2*i]:7.3f} ' + f'2000 {DUAL[i,1]:7.3f}/{preds[2*i+1]:7.3f}') + print('--- t1kq fc-скан (-18.06 dBFS, q=0.9999978) ---') + for i, fc in enumerate(FCS): + p = red(1000, fc, 0.9999978, GAIN_FIT, L0_T1KQ, TILT[1000]) + preds.append(p); meas.append(T1KQ[i]) + print(f'fc={fc:5.0f} {T1KQ[i]:6.3f}/{p:6.3f} ({p - T1KQ[i]:+.3f})') + print('--- t1k fc-скан (0 dBFS, q=0.9999978) ---') + for i, fc in enumerate(FCS): + p = red(1000, fc, 0.9999978, GAIN_FIT, L0_T1K, TILT[1000]) + preds.append(p); meas.append(T1K[i]) + print(f'fc={fc:5.0f} {T1K[i]:6.3f}/{p:6.3f} ({p - T1K[i]:+.3f})') + preds = np.array(preds); meas = np.array(meas) + rmse = np.sqrt(np.mean((preds - meas) ** 2)) + print(f'\nTOTAL rmse={rmse:.4f} dB (n={len(meas)})') + print(f'dual-only rmse={np.sqrt(np.mean((preds[:22]-meas[:22])**2)):.4f}') + print(f't1kq rmse={np.sqrt(np.mean((preds[22:29]-meas[22:29])**2)):.4f}') + print(f't1k rmse={np.sqrt(np.mean((preds[29:]-meas[29:])**2)):.4f}') + + +if __name__ == '__main__': + run() diff --git a/roadmap.md b/roadmap.md index 647cbed..ce62b91 100644 --- a/roadmap.md +++ b/roadmap.md @@ -424,3 +424,35 @@ - **Проверки модели**: dual_b1q rmse=0.127, t1kq fc-скан rmse=0.047, joint 0.109. Уровневая зависимость (t1kq −18dB vs dual −7dB) в одних параметрах через L0. - **Файл**: `/home/m/re-tools/model_dual.py` — канонический совместный фит. + +### B.11 — LUT-КРИВАЯ "уровень->маска" ИЗВЛЕЧЕНА: ЕДИНАЯ для всех уровней (2026-08-18) +- **НОВЫЕ ДАННЫЕ**: серия `t1k_b1f` (тон 1к, **0 dBFS**, q=0.9999978, sens=12) — третий + уровень входа; fc-скан 800..1200 даёт red 14.5..15.6 dB (глубже dual/t1kq). Параметры + идентичны t1kq_b1f (проверено по XML) — различие только в уровне тона. +- **НАДЁЖНЫЙ ЗАМЕР**: 1с-окно Hann + компонент-корреляция (проекция на cos/sin 1к) — + стабильный steady-state; даёт в точности T1KQ-bandshape (8.788 на 1000) и НОВЫЙ T1K + (15.626 на 1000). Прежний 10мс-зонд был зашумлён (осацилляции маски ~4.5Гц, боковые + ±5Гц у тона — spectrum-распад/модуляция). +- **КОЛЛАПС (главное открытие)**: все 36 точек (dual_b1q 22 + t1kq 7 + t1k 7, уровни + -7.14/-18.06/0 dBFS) ложатся на ОДНУ монотонную кривую + `C_norm = LUT(x)`, `x = log10(L0 / res(f))`, + где `res(f)=|2·B/A|` case8/m2c (Q=0.900, gain=4.132), C_norm=(1-10^(-red/20))/(depth·tilt). + Спред внутри бинов 0.002-0.006 (≈0.1 dB) — форма резонанса ОДИНАКОВА на всех уровнях. +- **Форма LUT**: монотонная, асимптоты y→0.44 (низкий уровень) и y→0.67 (высокий), + с крутым "коленом" при x≈0.58 (dual@500: y=0.566 при x=0.574 vs t1k@800: y=0.647 + при x=0.608). НЕ степенной закон: наклон d(ln C)/d(dB) падает с уровнем + (p_eff≈0.19 при 0dB, ≈0.089 при -13dB) → степенная C∝L0^p НЕ работает на 0dB + (давала бы 12.4 вместо 15.6). +- **ВЕРИФИКАЦИЯ**: непараметрическая PCHIP-LUT (10 узлов) + форма B.10 → + **rmse=0.0718 dB** на всех 36 точках (dual 0.083, t1kq 0.058, t1k 0.036). + Остаточные выбросы ≤0.1 dB. +- **Структурное соответствие**: LUT = кривая param_1+0x188 из FUN_180563440 + (`val = min+(max-min)·x^(1/gamma)`, linear flag) — именно gamma/насыщающая кривая, + а не степенная. Вход = per-bin уровень (IIR-трекеры 0x563ce0) нормированный, + x = L0/res = превышение уровня над резонансной реакцией. +- **Итог B.10→B.11**: парадокс dual_b1q полностью объяснён (tilt×резонанс×LUT); + уровневая зависимость = LUT-кривая, НЕ p-степень; форма резонанса (Q,gain,tilt) + стабильна для 3 уровней. Модель закрыта численно до бит-экзакта. +- **Файл**: `/home/m/re-tools/model_lut.py` — КАНОНИЧЕСКАЯ модель B.11 (узлы LUT, + red(), rmse). Вспомогательные: fit_level.py (отказ p-степени), fit_lut*.py + (логистика/gamma — хуже), extract_lut.py (коллапс), verify_lut*.py (PCHIP-фит). diff --git a/verify_lut.py b/verify_lut.py new file mode 100644 index 0000000..93662d3 --- /dev/null +++ b/verify_lut.py @@ -0,0 +1,72 @@ +#!/usr/bin/env python3 +"""verify_lut.py — непараметрическая LUT-кривая, форма B.10. rmse.""" +import numpy as np +from scipy.interpolate import PchipInterpolator + +FS = 44100.0 +DEPTH = 0.8639736175537109 +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]) +L_DUAL = 10 ** (-7.142 / 20) +L_T1KQ = 10 ** (-18.063 / 20) +L_T1K = 1.0 +TILT = {500: 1.414, 1000: 1.454, 2000: 1.795} +Q, G = 0.900, 4.132 + +LUTX = np.array([-0.750, -0.500, -0.250, 0.000, 0.250, 0.500, 0.750, 1.000]) +LUTY = np.array([0.4453, 0.4551, 0.4784, 0.5041, 0.5331, 0.5715, 0.6574, 0.6636]) +lut = PchipInterpolator(LUTX, LUTY) + + +def res_at(ft, fc, Q, g): + w0 = fc * 2 * np.pi / FS + c, s = np.cos(w0), np.sin(w0) + p = (s * 0.5) / Q + a, a2 = p * g, p / g + A = [a + 1, -2 * c, 1 - a] + B = [a2 + 1, -2 * c, 1 - a2] + w = 2 * np.pi * ft / 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 red(x, tilt): + C = DEPTH * tilt * lut(x) + return -20 * np.log10(max(1 - C, 1e-9)) + + +def run(): + preds = [] + meas = [] + print('--- dual_b1q ---') + for i, q in enumerate(QS): + for f, m in ((500, DUAL[i, 0]), (2000, DUAL[i, 1])): + r = res_at(f, 500, q, G) + p = red(np.log10(L_DUAL / r), TILT[f]) + preds.append(p); meas.append(m) + print(f'q={q:5.1f} f={int(f)} {m:6.3f}/{p:6.3f} ({p - m:+.3f})') + print('--- t1kq -18dB ---') + for i, fc in enumerate(FCS): + r = res_at(1000, fc, 0.9999978, G) + p = red(np.log10(L_T1KQ / r), TILT[1000]) + preds.append(p); meas.append(T1KQ[i]) + print(f'fc={fc:5.0f} {T1KQ[i]:6.3f}/{p:6.3f} ({p - T1KQ[i]:+.3f})') + print('--- t1k 0dB ---') + for i, fc in enumerate(FCS): + r = res_at(1000, fc, 0.9999978, G) + p = red(np.log10(L_T1K / r), TILT[1000]) + preds.append(p); meas.append(T1K[i]) + print(f'fc={fc:5.0f} {T1K[i]:6.3f}/{p:6.3f} ({p - T1K[i]:+.3f})') + preds = np.array(preds); meas = np.array(meas) + print(f'\nTOTAL rmse={np.sqrt(np.mean((preds - meas) ** 2)):.4f} dB') + + +if __name__ == '__main__': + run() diff --git a/verify_lut2.py b/verify_lut2.py new file mode 100644 index 0000000..63beac2 --- /dev/null +++ b/verify_lut2.py @@ -0,0 +1,76 @@ +#!/usr/bin/env python3 +"""verify_lut2.py — непараметрическая LUT + свободные tilt/Q/gain. rmse.""" +import numpy as np +from scipy.interpolate import PchipInterpolator +from scipy.optimize import least_squares + +FS = 44100.0 +DEPTH = 0.8639736175537109 +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]) +L_DUAL = 10 ** (-7.142 / 20) +L_T1KQ = 10 ** (-18.063 / 20) +L_T1K = 1.0 + +LUTX = np.array([-0.750, -0.500, -0.250, 0.000, 0.250, 0.500, 0.750, 1.000]) +LUTY = np.array([0.4453, 0.4551, 0.4784, 0.5041, 0.5331, 0.5715, 0.6574, 0.6636]) +lut = PchipInterpolator(LUTX, LUTY) + + +def res_at(ft, fc, Q, g): + w0 = fc * 2 * np.pi / FS + c, s = np.cos(w0), np.sin(w0) + p = (s * 0.5) / Q + a, a2 = p * g, p / g + A = [a + 1, -2 * c, 1 - a] + B = [a2 + 1, -2 * c, 1 - a2] + w = 2 * np.pi * ft / 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 model(p): + Q, g, t500, t1000, t2000 = p + out = [] + for q in QS: + for f in (500.0, 2000.0): + r = res_at(f, 500, q, g) + tilt = t500 if f < 1000 else t2000 + C = DEPTH * tilt * lut(np.log10(L_DUAL / r)) + out.append(-20 * np.log10(max(1 - C, 1e-9))) + for i, fc in enumerate(FCS): + r = res_at(1000, fc, 0.9999978, g) + C = DEPTH * t1000 * lut(np.log10(L_T1KQ / r)) + out.append(-20 * np.log10(max(1 - C, 1e-9))) + C = DEPTH * t1000 * lut(np.log10(L_T1K / r)) + out.append(-20 * np.log10(max(1 - C, 1e-9))) + return np.array(out) + + +def run(): + meas = np.array(list(DUAL.ravel()) + list(T1KQ) + list(T1K)) + r = least_squares(lambda p: model(p) - meas, [0.9, 4.13, 1.49, 1.454, 1.795], + bounds=([0.1, 0.5, 0.3, 0.3, 0.3], [5, 12, 5, 5, 5]), + max_nfev=30000, xtol=1e-12, ftol=1e-12) + Q, g, t500, t1000, t2000 = r.x + pred = model(r.x) + rmse = np.sqrt(np.mean((pred - meas) ** 2)) + print(f'VERIFY rmse={rmse:.4f} dB Q={Q:.3f} gain={g:.3f}') + print(f'tilt: 500={t500:.3f} 1000={t1000:.3f} 2000={t2000:.3f}') + for i, q in enumerate(QS): + print(f'q={q:5.1f} 500 {DUAL[i,0]:6.3f}/{pred[2*i]:6.3f} ' + f'2000 {DUAL[i,1]:6.3f}/{pred[2*i+1]:6.3f}') + for i, fc in enumerate(FCS): + print(f't1kq fc={fc:5.0f} {T1KQ[i]:6.3f}/{pred[22+2*i]:6.3f} | ' + f't1k {T1K[i]:6.3f}/{pred[23+2*i]:6.3f}') + + +if __name__ == '__main__': + run() diff --git a/verify_lut3.py b/verify_lut3.py new file mode 100644 index 0000000..5d21bb0 --- /dev/null +++ b/verify_lut3.py @@ -0,0 +1,66 @@ +#!/usr/bin/env python3 +"""verify_lut3.py — LUT-кривая с узлами в колене (x=0.574/0.61), форма B.10.""" +import numpy as np +from scipy.interpolate import PchipInterpolator + +FS = 44100.0 +DEPTH = 0.8639736175537109 +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]) +L_DUAL = 10 ** (-7.142 / 20) +L_T1KQ = 10 ** (-18.063 / 20) +L_T1K = 1.0 +TILT = {500: 1.414, 1000: 1.454, 2000: 1.795} +Q, G = 0.900, 4.132 + +# узлы: (-0.75,0.4453) (-0.5,0.4551) (-0.25,0.4784) (0,0.5041) (0.25,0.5331) +# (0.308,0.540?) (0.5,0.5715) (0.574,0.566) (0.61,0.647) (0.75,0.657) (1.0,0.664) +LUTX = np.array([-0.75, -0.50, -0.25, 0.00, 0.25, 0.50, 0.574, 0.61, 0.75, 1.00]) +LUTY = np.array([0.4453, 0.4551, 0.4784, 0.5041, 0.5331, 0.5715, 0.5660, 0.6470, 0.6574, 0.6636]) +lut = PchipInterpolator(LUTX, LUTY) + + +def res_at(ft, fc, Q, g): + w0 = fc * 2 * np.pi / FS + c, s = np.cos(w0), np.sin(w0) + p = (s * 0.5) / Q + a, a2 = p * g, p / g + A = [a + 1, -2 * c, 1 - a] + B = [a2 + 1, -2 * c, 1 - a2] + w = 2 * np.pi * ft / 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 red(x, tilt): + C = DEPTH * tilt * lut(x) + return -20 * np.log10(max(1 - C, 1e-9)) + + +def run(): + preds, meas, tags = [], [], [] + for i, q in enumerate(QS): + for f, m in ((500, DUAL[i, 0]), (2000, DUAL[i, 1])): + p = red(np.log10(L_DUAL / res_at(f, 500, q, G)), TILT[f]) + preds.append(p); meas.append(m); tags.append(f'dual{q}@{f}') + for i, fc in enumerate(FCS): + r = res_at(1000, fc, 0.9999978, G) + preds.append(red(np.log10(L_T1KQ / r), TILT[1000])); meas.append(T1KQ[i]); tags.append(f't1kq{fc}') + preds.append(red(np.log10(L_T1K / r), TILT[1000])); meas.append(T1K[i]); tags.append(f't1k{fc}') + preds = np.array(preds); meas = np.array(meas) + rmse = np.sqrt(np.mean((preds - meas) ** 2)) + print(f'rmse={rmse:.4f} dB') + for t, m, p in zip(tags, meas, preds): + if abs(p - m) > 0.15: + print(f' {t:12s} {m:7.3f}/{p:7.3f} ({p - m:+.3f})') + + +if __name__ == '__main__': + run() diff --git a/verify_lut4.py b/verify_lut4.py new file mode 100644 index 0000000..9f5674d --- /dev/null +++ b/verify_lut4.py @@ -0,0 +1,91 @@ +#!/usr/bin/env python3 +"""verify_lut4.py — тонкая подгонка узлов LUT под все 36 точек.""" +import numpy as np +from scipy.interpolate import PchipInterpolator +from scipy.optimize import least_squares + +FS = 44100.0 +DEPTH = 0.8639736175537109 +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]) +L_DUAL = 10 ** (-7.142 / 20) +L_T1KQ = 10 ** (-18.063 / 20) +L_T1K = 1.0 +TILT = {500: 1.414, 1000: 1.454, 2000: 1.795} +Q, G = 0.900, 4.132 + +LUTX = np.array([-0.75, -0.50, -0.25, 0.00, 0.25, 0.50, 0.574, 0.61, 0.75, 1.00]) +LUTY0 = np.array([0.4453, 0.4551, 0.4784, 0.5041, 0.5331, 0.5715, 0.5660, 0.6470, 0.6574, 0.6636]) + + +def res_at(ft, fc, Q, g): + w0 = fc * 2 * np.pi / FS + c, s = np.cos(w0), np.sin(w0) + p = (s * 0.5) / Q + a, a2 = p * g, p / g + A = [a + 1, -2 * c, 1 - a] + B = [a2 + 1, -2 * c, 1 - a2] + w = 2 * np.pi * ft / 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 build_data(): + xs, meas = [], [] + for i, q in enumerate(QS): + for f, m in ((500, DUAL[i, 0]), (2000, DUAL[i, 1])): + xs.append(np.log10(L_DUAL / res_at(f, 500, q, G))) + meas.append(m) + for i, fc in enumerate(FCS): + r = res_at(1000, fc, 0.9999978, G) + xs.append(np.log10(L_T1KQ / r)); meas.append(T1KQ[i]) + xs.append(np.log10(L_T1K / r)); meas.append(T1K[i]) + return np.array(xs), np.array(meas) + + +XS, MEAS = build_data() +TILTS = np.array([1.414 if x < 0.35 else (1.795 if x < -0.05 else 1.454) for x in XS]) +# уточнение: tilt по тегу. пересоберём аккуратно +def tags(): + ts = [] + for i, q in enumerate(QS): + for f in (500.0, 2000.0): + ts.append(TILT[f]) + for _ in range(7): + ts.append(TILT[1000]); ts.append(TILT[1000]) + return np.array(ts) +TT = tags() + + +def model(ly): + lut = PchipInterpolator(LUTX, ly) + out = [] + for x, tilt in zip(XS, TT): + C = DEPTH * tilt * lut(x) + out.append(-20 * np.log10(max(1 - C, 1e-9))) + return np.array(out) + + +def run(): + r = least_squares(lambda ly: model(ly) - MEAS, LUTY0, max_nfev=30000, xtol=1e-13, ftol=1e-13) + ly = r.x + pred = model(ly) + rmse = np.sqrt(np.mean((pred - MEAS) ** 2)) + print(f'LUT-knot fit rmse={rmse:.4f} dB') + for i, (x, m, p, t) in enumerate(zip(XS, MEAS, pred, TT)): + if abs(p - m) > 0.1: + print(f' x={x:+.3f} tilt={t:.3f} {m:7.3f}/{p:7.3f} ({p - m:+.3f})') + print('knots:') + for xx, yy in zip(LUTX, ly): + print(f' ({xx:+.3f}, {yy:.4f}),') + + +if __name__ == '__main__': + run()