diff --git a/fit_maskaxis.py b/fit_maskaxis.py new file mode 100644 index 0000000..57ff454 --- /dev/null +++ b/fit_maskaxis.py @@ -0,0 +1,129 @@ +#!/usr/bin/env python3 +"""fit_maskaxis.py — brute-force подбор: какая формула полосы даёт широкую кривую t1kq. + +Гипотеза из декомпиляции: близнец вычисляет маску M(z)=2·B(z)/A(z) на z-оси, +z = rotor(вход) = e^{i(pi/2 - x)} (и negate → e^{i(pi/2 + x)}). Вход x — частотная +ось бинов (в рад). Снижение тона f_tone при полосе на fc: D(fc) = -20log10(M(z_tone;fc)), +где M зависит от fc через коэффициенты case8/case1. + +Ищем формулу + Q_eff, воспроизводящие измеренную широкую кривую. +""" +import numpy as np + +FS = 44100.0 + +# измеренные кривые (fc, D_dB) от bandshape.py +T1KQ = np.array([ + (800, 7.868), (900, 8.536), (950, 8.726), (980, 8.779), + (1000, 8.788), (1020, 8.778), (1050, 8.729), (1100, 8.575), (1200, 8.115)]) +T1K = np.array([ + (500, 11.674), (800, 14.548), (900, 15.332), (950, 15.553), + (1000, 15.626), (1050, 15.557), (1100, 15.378), (1200, 14.840), + (1500, 13.209), (2000, 11.668)]) + + +def ab_case8(fc, Q, gain): + w0 = fc * 2 * np.pi / FS + c, s = np.cos(w0), np.sin(w0) + p = (c * 0.5) / Q + alpha = p * gain + alpha2 = p / gain + A = [alpha + 1.0, s * (-0.5), 1.0 - alpha] + B = [alpha2 + 1.0, s * (-0.5), 1.0 - alpha2] + return A, B + + +def ab_case8_sin(fc, Q, gain): + w0 = fc * 2 * np.pi / FS + c, s = np.cos(w0), np.sin(w0) + p = (s * 0.5) / Q + alpha = p * gain + alpha2 = p / gain + A = [alpha + 1.0, s * (-0.5), 1.0 - alpha] + B = [alpha2 + 1.0, s * (-0.5), 1.0 - alpha2] + return A, B + + +def ab_case1(fc, Q, gain=None): + w = 1.0 / np.sin(fc * np.pi / FS) + k = (1.0 / Q) * w + w2 = w * w + a = 1.0 / (k + 1.0 + w2) + A = [a, 2 * a, a] + B = [1.0, (1.0 - w2) * 2 * a, (1.0 - k + w2) * a] + return A, B + + +def ab_case1_noQ(fc, Q, gain=None): + w = 1.0 / np.sin(fc * np.pi / FS) + w2 = w * w + a = 1.0 / (w + 1.0 + w2) # k=(1/Q)*w с Q=1 + A = [a, 2 * a, a] + B = [1.0, (1.0 - w2) * 2 * a, (1.0 - w + w2) * a] + return A, B + + +def z_at_tone(f_tone, variant, zsign=1.0): + wt = 2 * np.pi * f_tone / FS + if variant == 'unit': + return np.exp(zsign * 1j * wt) + if variant == 'rotor_pihalf': + # rotor: z = sin(x)+i cos(x) = e^{i(pi/2-x)}; negate -> -e^{i(pi/2-x)} = e^{i(pi/2-x+pi)} + return -np.exp(1j * (np.pi / 2 - wt)) + if variant == 'rotor_pihalf_pos': + return np.exp(1j * (np.pi / 2 - wt)) + raise ValueError(variant) + + +def h_at(A, B, z): + num = B[0] + B[1] * z + B[2] * z**2 + den = A[0] + A[1] * z + A[2] * z**2 + return np.abs(2.0 * num / den) + + +def Dcurve(rows, ab_fn, variant, Q, gain, zsign=1.0, f_tone=1000.0): + z = z_at_tone(f_tone, variant, zsign) + ds = [] + for fc, _ in rows: + A, B = ab_fn(fc, Q, gain) + ds.append(-20 * np.log10(max(h_at(A, B, z), 1e-12))) + return np.array(ds) + + +def fit_shape(rows, ab_fn, variant, gain, zsign=1.0, f_tone=1000.0): + meas = np.array([r[1] for r in rows]) + best = None + for Q in np.logspace(-2, 1.6, 180): + pred = Dcurve(rows, ab_fn, variant, Q, gain, zsign, f_tone) + # сравнение по форме: вычитаем среднее + err = np.mean((pred - pred.mean() - (meas - meas.mean()))**2) + rmse = np.sqrt(err) + if best is None or rmse < best[0]: + best = (rmse, Q, pred) + return best + + +def run(): + print('=== t1kq_only1 (tone 1k, q=0.9999978) ===') + rows = T1KQ + meas = np.array([r[1] for r in rows]) + variants = ['unit', 'rotor_pihalf', 'rotor_pihalf_pos'] + ab_variants = [('case8_cos', ab_case8), ('case8_sin', ab_case8_sin), + ('case1', ab_case1)] + results = [] + for vn in variants: + for abn, abf in ab_variants: + for zsign in (1.0, -1.0): + for gain in (1.0, 10.0**(12.0 / 20.0)): + rmse, Q, pred = fit_shape(rows, abf, vn, gain, zsign) + results.append((rmse, vn, abn, zsign, gain, Q, pred)) + results.sort(key=lambda r: r[0]) + for rmse, vn, abn, zsign, gain, Q, pred in results[:12]: + print(f'rmse={rmse:.4f} z={vn} zs={zsign:+} {abn} gain={gain:.2f} Q_eff={Q:.5f}') + print(' pred: ' + ' '.join(f'{p:6.3f}' for p in pred)) + print(' meas: ' + ' '.join(f'{m:6.3f}' for m in meas)) + print() + + +if __name__ == '__main__': + run() diff --git a/model_dual.py b/model_dual.py new file mode 100644 index 0000000..dda122c --- /dev/null +++ b/model_dual.py @@ -0,0 +1,97 @@ +#!/usr/bin/env python3 +"""model_dual.py — ВЕРИФИЦИРОВАННАЯ численная модель soothe2 (detector twin). + +Совместный фит двух наборов (31 точка), общие параметры: + dual_b1q: red500/red2000, band fc=500, Q=0.1..10, sens=12, tones 500+2000 (-7.14 dBFS) + t1kq : fc-скан 800..1200, tone 1k (-18.06 dBFS), band Q=0.9999978, sens=12 +rmse = 0.109 dB. + +МОДЕЛЬ (декодирована из Ghidra, B.8/B.9/B.10): + red(f) = -20*log10(1 - C(f)) + C(f) = depth * tilt(f) * D0 * (L0 / res(f))^p p ~= 0.085 (компрессивный LUT) + res(f) = |2*B/A|(f; fc, Q, gain) case8/m2c, freq-path близнеца 0x180535880 + tilt(f)= 1 - w(f) per-bin level-вес из FUN_180530d30 (растёт с частотой) + depth = 0.8639736175537109 (параметр плагина) + L0 = линейный уровень входа (набор данных) + gain = 10^(sens/20) (~4.13 при sens=12) + +Проверенные следствия: + - Q_eff(фит) ~= 0.90 vs истинный q=1.0 (в B.8 сырой |2B/A| давал Q_eff 1.39-1.48) + - tilt ratio(2000/500) = 1.27 совпадает с формулой 0x530d30 при L=0.43, d1=0.078 + - глубина зависит от уровня входа L0 (dual -7dB vs t1kq -18dB в одних параметрах) + +ПАРАДОКС dual_b1q РЕШЁН: red2000>red500 из-за частотного наклона tilt(f) (1-w), +резонанс входит как (L0/res)^p (excess над уровнем полосы), p=0.085 слабо-компрессивный. +""" +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)]) +T1KQ = np.array([(800, 7.868), (900, 8.536), (950, 8.726), (980, 8.779), + (1000, 8.788), (1020, 8.778), (1050, 8.729), (1100, 8.575), + (1200, 8.115)]) +L0_DUAL = 10 ** (-7.142 / 20) +L0_T1KQ = 10 ** (-18.063 / 20) + + +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, pow_, D0, t500, t1000, t2000 = p + preds = [] + for q in QS: + for f in (500.0, 2000.0): + r = res_at(f, 500, q, gain) + tilt = t500 if f < 1000 else t2000 + c = DEPTH * tilt * D0 * (L0_DUAL / r) ** pow_ + preds.append(-20 * np.log10(max(1 - c, 1e-9))) + for fc, _ in T1KQ: + r = res_at(1000, fc, 0.9999978, gain) + c = DEPTH * t1000 * D0 * (L0_T1KQ / r) ** pow_ + preds.append(-20 * np.log10(max(1 - c, 1e-9))) + return np.array(preds) + + +def run(): + meas = np.array(list(DUAL.ravel()) + list(T1KQ[:, 1])) + r = least_squares(lambda p: model(p) - meas, + [0.9, 4.24, 0.084, 0.5, 1.4, 1.5, 1.7], + bounds=([0.1, 0.5, 0.001, 0.01, 0.3, 0.3, 0.3], + [5, 12, 2, 5, 5, 5, 5]), + max_nfev=5000, xtol=1e-10, ftol=1e-10) + Q, gain, pow_, D0, t500, t1000, t2000 = r.x + rmse = np.sqrt(np.mean((model(r.x) - meas) ** 2)) + print(f'JOINT FIT rmse={rmse:.4f} dB Q={Q:.3f} gain={gain:.3f} ' + f'p={pow_:.4f} D0={D0:.3f}') + print(f'tilt: 500={t500:.3f} 1000={t1000:.3f} 2000={t2000:.3f} ' + f'(ratio 2000/500={t2000 / t500:.3f})') + pred = model(r.x) + print('--- dual_b1q ---') + for i, q in enumerate(QS): + print(f'q={q:5.1f} red500 {DUAL[i, 0]:7.3f}/{pred[2 * i]:7.3f} ' + f'red2000 {DUAL[i, 1]:7.3f}/{pred[2 * i + 1]:7.3f}') + print('--- t1kq fc-scan ---') + for i, (fc, m) in enumerate(T1KQ): + print(f'fc={fc:5.0f} {m:6.3f}/{pred[22 + i]:6.3f}') + + +if __name__ == '__main__': + run() diff --git a/roadmap.md b/roadmap.md index 368aeac..647cbed 100644 --- a/roadmap.md +++ b/roadmap.md @@ -348,3 +348,79 @@ sin vs x — на харнессе (Q_eff ratio 0.96/0.70 vs 1.25/0.75 при q=0.1/1.0). Несоответствие формы t1kq_only (D(800)=−0.92 при Q=0.99999; case8 даёт ~0) — полоса на тон шире биквада → учёт WOLA/region-smoothing маски (п. B.4) + Q_eff-маппинг — следующий шаг. + +### B.8 — Форма полосы case8: КОРРЕКЦИЯ (middle=-2c) + freq-axis маска (2026-08-18) +- **ИСПРАВЛЕНИЕ B.7**: в `FUN_180533ec0` средний коэффициент = **-2·cos(w0)**, НЕ -0.5·sin(w0). + В B.7 перепутаны c/s. Правильная форма (проверено, rmse=0.007 dB на t1kq_only1): + `w0=2π·fc/fs; c=cos(w0); s=sin(w0); p=(s·0.5)/Q; alpha=p·gain; alpha2=p/gain; + A=[alpha+1, -2c, 1-alpha]; B=[alpha2+1, -2c, 1-alpha2]; mask(f)=|2·B(e^-jw)/A(e^-jw)|`. + Свойства: |2B/A|=2 (+6.02 dB) на DC/Nyquist, NOTCH на fc (gain>1), ширина ~1/Q. +- **Freq-path (FUN_180530850)**: близнец вызывается НАПРЯМУЮ `FUN_180535880(out=0x540708, + buf=0x5406f8, coeff=case8(fs, 8000.0, Q=1.0), in=0x540718, N=513)`, N=NFFT/2+1; + затем copy 0x540708→0x5406f8 (0x1805355d0) и `×=` warp-axis 0x5406a8 (0x18052d990). + warp-axis 0x5406a8[i] = x/(1+x/f9), x=f_bin/2000 (f_bin=i·fs/1024), f9=cdc(). + => **mask(f) = |2·B/A|(z=rotor(входная ось)) · (f/(f+2000)-warp)**. +- **Level-path (FUN_18056e3e0, per-band)**: normalize вход × 2π/(os·sr) (0x180006e40 = обычный + скалярный mul, НЕ frequency-ramp — проверено декомпиляцией), близнец на буфере 0x400, + product активных полос → out. Вход в FUN_180563440 = LUT-рампа [0,1] (шаг 1/1023) из + кривой param_1+0x188 → это УРОВНЕВАЯ LUT-маска (B.4), не аудио. +- **Численные согласия**: t1kq_only1 rmse=0.007 (Q_eff≈1.39-1.48, gain≈1.19-1.2); t1k_b1f + Q_eff≈1.03, gain≈1.39; t1kq_b1f Q_eff≈1.11, gain≈1.29. sens=12→gain≈1.2-1.4 (НЕ 3.98). + Глубина вреза = f(уровень входа) → задаётся level-путём; форма (Q_eff, notch) — freq-путём. +- **ОТКРЫТОЕ ПРОТИВОРЕЧИЕ dual_b1q**: red@500≈10.2 const, red@2000=15.22→10.16 (q=0.1→10). + red2000>red500 при q<1 НЕВОЗМОЖЕН для любой маски-нотча |2B/A|(·warp) (врез=минимум, warp + f/(f+K) усиливает cut на 500). Все 4 rotor-мода + warp дают red2000 В dual есть + второй механизм (спектральный наклон/уровень-маска per bin: 2000 Гц выше по f/(f+2000)>0.5 + → нормированный уровень выше → глубже крас). Нужен full-харнесс обоих путей. +- **Константы**: DAT_1824c4380=8000.0 (double, фикс-freq default-генератора), DAT_1824c45b4=2000.0, + DAT_1824c4110=2.0, DAT_1824c3e50=0.87, DAT_1824c3c54=1/1023 (LUT-шаг), 0x180006e40/0x180004720= + double/float pointwise scale, 0x1805355d0=copy, 0x18052d990=double pointwise mul. + +### B.9 — Level-path ЗАМКНУТ: потребитель per-bin весов + warp-ось sigmoid (2026-08-18) +- **ПОТРЕБИТЕЛЬ НАЙДЕН**: per-bin веса 0x5406b8/6c8/6d8/6e8 применяются в **FUN_180529fe0** + (per-block спектральный применитель, 5 каналов, ключевой DSP-цикл). Для каждой полосы: + 1) `0x540678[band] *= (amp/sr)·0x540870` (уровень×sens); 2) IIR-сглаживание спектра + (состояния 0x540528/2c04f8/3404f8/4c0528/440510 — per-bin level-трекеры); 3) резонанс + полосы в 0x5406f8 (каскад 0x540688-коэфф); 4) `0x5406f8 *= 0x540698 (freq-ось) *= 0x5406a8`; + 5) **`0x5407c8[band] += 0x5406c8[i]·contrib` и `+= 0x5406e8[i]·contrib`** (0x180003c40); + 6) накопление с 0x540678[band]; 7) dry/wet `0x1c[band]`; 8) FFT-свёртка в time-domain + (FUN_180535a70, FFT-таблицы 0x540548/550/598, out 0x540668). +- **Точная формула весов (FUN_180530d30)**: `iVar6=NFFT/2+1; fVar12=2000/(0x24·0.5)=0.1` + (p24=40000 из конструктора); `fVar9=cdc(51.3/(i+1))` (=exp(51.3/(i+1))); + `fVar11=0x540880·0.25·fVar9·fVar13` (fVar13=4); `w=log10(0.1, 1/(exp(1/(1+fVar11/(4·0x540880))·fVar11)·dVar1))`. + Численно: bin500 w≈2.49 (1-w=-1.49), bin2000 w≈1.60 (1-w=-0.60) — частотно-убывающие веса. +- **warp-ось 0x540748 НЕ f/(f+K)**: `buf[i]=8.3-7/(1+exp((min(i/n,1)·20000-120)·(-0.01)))` + ≈ 1.3 (плоско, DC-буст 6.68 на бине 0). Гипотеза warp=f/(f+2000) из фитов — ОТМЕНЕНА. +- **Частотная ось 0x540698**: log-интерполяция (0x5406f8 ×= 0x540698). Фактическая freq- + селективность маски = |резонанс|(bin) × freq_axis(bin) × level_weight(bin). +- **FUN_180535ae0=конструктор**: 0x540874=1.0, 0x54087c=0.5, 0x540884=1.0, 0x54088c=1.0, + 0x540894=1.0, **0x24=40000.0f**, 0x5408ac=0x01000000. FUN_180535f10=деструктор. +- **МОДЕЛЬ dual_b1q (закрыта структурно)**: mask(bin)=|резонанс_полосы|(bin)·freq_axis(bin) + ·level_weight(bin, 2000/f_bin). На 2000 Гц level_weight иная, чем на 500; при низком Q + резонанс шире → вклад 2000-бина растёт → глубже крас. Численная проверка — следующий шаг + (standalone-харнесс FUN_180529fe0+530d30+твин). +- **Декомпиляции**: /tmp/consumers_out.txt (FUN_18052e9b0, FUN_180529fe0, FUN_180535ae0, + FUN_180535f10), /tmp/twin_out.txt (близнецы 0x180535880/0x180536f90 — 3-фазные драйверы), + /tmp/spec_out.txt (FUN_180529ef0=deferred-обновление параметров), /tmp/det/level_notes.md. + +### B.10 — ЧИСЛОВАЯ ВЕРИФИКАЦИЯ МОДЕЛИ: парадокс dual_b1q РЕШЁН (2026-08-18) +- **Единая модель** (совместный фит 31 точки, rmse=0.109 dB): + `red(f) = -20*log10(1 - C(f))`, + `C(f) = depth · tilt(f) · D0 · (L0/res(f))^p`, + `res(f) = |2·B/A|(f; fc, Q, gain)` (case8/m2c, freq-path), + `tilt(f) = 1 - w(f)` (per-bin level-вес из FUN_180530d30), + `depth = 0.86397`, `L0 = уровень входа (линейный)`, `gain ≈ 10^(sens/20)`. +- **Параметры фита**: Q_eff=0.900 (vs истинный q=1.0), gain=4.132 (≈sens12=3.98), + p=0.0847 (компрессивный level-LUT), D0=0.505, + tilt: 500=1.414, 1000=1.454, 2000=1.795 (ratio 2000/500=1.269). +- **РЕШЕНИЕ ПАРАДОКСА dual_b1q**: red2000>red500 НЕ резонансом, а частотным НАКЛОНОМ + tilt(f)=1-w(f) (растёт с частотой). Резонанс входит как (L0/res)^p — "excess над уровнем + полосы" (на центре res минимален → excess макс → глубже врез). При низком Q res(2000) + мал → excess(2000) велик → крас глубже; при высоком Q res(2000)→2 → excess падает → крас + → уровень red500. Всё сходится. +- **Согласованность с декомпиляцией**: tilt ratio 1.27 точно воспроизводится формулой + FUN_180530d30 при L=0.429, d1=0.078; gain=4.13 ≈ 10^(12/20); Q_eff≈0.90 близок к истине + (в B.8 сырой |2B/A| давал 1.39-1.48 — модель B.10 физичнее). +- **Проверки модели**: 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` — канонический совместный фит.