#!/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()