Matiq/soothe2-re
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

B.11: LUT-кривая уровень->маска извлечена (rmse 0.072 dB, 36 точек, 3 уровня): коллапс dual+t1kq+t1k на ОДНУ кривую x=log10(L0/res); степенной закон C~L0^p НЕ работает на 0dB (p_eff 0.19 vs 0.085) - LUT насыщающая gamma из FUN_180563440 (param_1+0x188); форма резонанса стабильна; model_lut.py канонический + roadmap B.11

Browse Source
This commit is contained in:
Matiq
2026-08-18 13:02:56 +03:00
parent 9abecb60d5
commit 32a63955f8
13 changed files with 1021 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files
extract_lut.py
+73
View File
@@ -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}),')
fit_level.py
+86
View File
@@ -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()
fit_lut.py
+87
View File
@@ -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()
fit_lut2.py
+82
View File
@@ -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()
fit_lut3.py
+84
View File
@@ -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()
fit_lut4.py
+80
View File
@@ -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()
fit_lut5.py
+84
View File
@@ -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()
model_lut.py
+108
View File
@@ -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()
roadmap.md
+32
View File
@@ -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-фит).
verify_lut.py
+72
View File
@@ -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()
verify_lut2.py
+76
View File
@@ -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()
verify_lut3.py
+66
View File
@@ -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()
verify_lut4.py
+91
View File
@@ -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()