252 lines
10 KiB
Python
252 lines
10 KiB
Python
#!/usr/bin/env python3
|
||
"""cascade_sim.py — структурный симулятор тракта маски soothe2 (float-путь).
|
||
|
||
Цель (24mm3): воспроизвести scr(f)=дизайн-сигнал детектора; применённая
|
||
маска = exp(γ·scr), γ=1.760561 (точно), trk@688=exp(scr@628) бит-в-бит.
|
||
|
||
Структура канонической цепи (BLOCKMAP 24hh/24ii + 24mm2):
|
||
шаг 9a: vec698 *= (1 − [54087c]) ; zero при дефолтах
|
||
шаг 9b: vec6f8 += [54087c]·0.8 ; xmm10=0.8 @1824c3e28
|
||
шаг 9c: bands_curve_i /= ... divide-ядро ; dst=678i, A/B уточняются
|
||
шаг 10: vec6f8 = bands_curve_i − ACC_i ; dc40, ACC @таблицы 0x5407c8
|
||
шаг 11: fma att/rel (тройки re/im/coef) ; коэф @6c8/6e8
|
||
шаг 12: COPY ; memcpy
|
||
шаг 13: зеркало 9
|
||
шаг 14: expf(bands_curve); bands_curve += (−1) ; ПОРЯДОК исправлен 24mm2
|
||
шаг 15: bands_curve *= track_i ; th2000 array-mul
|
||
шаг 16: bands_curve *= kWarp@[5406a8]
|
||
шаг 17: expf ещё раз ; call-site 52b32c
|
||
пост-17: exp-вариант(140a40) + pow?(140b00)
|
||
FIR-секция: кривая-float(140b30→1803831c0) + sincos-twiddle(140aa0)
|
||
|
||
ЯДРА (структурная фаза — математически точные numpy-эквиваленты;
|
||
канонический C++ порт = инструкци-точная транскрипция, см. BLOCKMAP 24mm2):
|
||
"""
|
||
import numpy as np
|
||
import glob
|
||
import os
|
||
|
||
GAMMA = 1.760561 # 24mm3: показатель степени, rms фита 0 на чистых кадрах
|
||
N = 2049 # число бинов полной сетки
|
||
|
||
|
||
# ---------------------------------------------------------------- ядра ----
|
||
def k_exp(x):
|
||
"""expf-ядро 180296c80. Структурная фаза: np.exp.
|
||
Каноническая формула (для C++ порта, FMA-точно):
|
||
n = fma(log2e_hi=1.4427f, x, 12582912.0f); k = n - MAGIC
|
||
r = (x - 0.693146f*k) - 1.42861e-06f*k
|
||
p = (((0.00829172f*r+0.0418735f)*r+0.166674f)*r+0.499994f)*r+1)*r+1
|
||
out = bits((k<<23) + bits(p)); guard |x|>87.3365 -> slow path
|
||
"""
|
||
return np.exp(x)
|
||
|
||
|
||
def k_exp_exact(x):
|
||
"""Bit-exact expf 180296c80 (BLOCKMAP:569) — float32 FMA poly, guard 87.3365.
|
||
Python структурный прокси: np.exp (float64) — точный C++ порт в dsp/exp2.cpp
|
||
через fmaf+бит-манипуляции (1824c...), ошибка <0.5 ulp vs плагин."""
|
||
return np.exp(np.asarray(x, dtype=np.float64)).astype(np.float32)
|
||
|
||
|
||
def k_div(a, b):
|
||
"""divide-ядро 1803a06a0: dst = B/A (~0.5 ulp, rcp+таблицы+полином).
|
||
Структурная фаза: точное деление."""
|
||
return b / a
|
||
|
||
|
||
def k_div_exact(a, b):
|
||
"""Bit-exact DIVIDE 1803a06a0 (BLOCKMAP:580) — rcp+quant+vpermps+poly.
|
||
Структурный прокси: точное деление; C++ порт копирует 0.5ulp полином
|
||
0.207515 -0.241687 0.288535 -0.360671 ... 0.240264 через vpermps tbl@1821269c0."""
|
||
return np.asarray(b, dtype=np.float64) / np.asarray(a, dtype=np.float64)
|
||
|
||
|
||
# ------------------------------------------------------------ данные -----
|
||
def load_tract(path):
|
||
"""tract_*.txt: k am res lvl_raw band_level prewarp w"""
|
||
t = np.loadtxt(path)
|
||
return {'am': t[:, 1], 'res': t[:, 2], 'lvl': t[:, 3]}
|
||
|
||
|
||
def load_frame(npz):
|
||
"""Слоты кадра rendersnap2 → dict[int, np.ndarray]."""
|
||
d = np.load(npz)
|
||
out = {}
|
||
for k in d.keys():
|
||
if k.startswith('0x'):
|
||
out[int(k[2:], 16)] = d[k]
|
||
return out, d['t_snap']
|
||
|
||
|
||
def pick_clean_frame(ds_dir, min_bins=8):
|
||
"""Отбор чистых стационарных кадров по фазам (24mm3):
|
||
возвращает лучший на фазе γ* (~1.7606, маска применена)
|
||
и лучший на фазе γ=1 (степень ещё не применена)."""
|
||
classes = {'gamma': None, 'identity': None}
|
||
for f in sorted(glob.glob(os.path.join(ds_dir, 'ph*.npz'))):
|
||
try:
|
||
S, ts = load_frame(f)
|
||
except Exception:
|
||
continue
|
||
if not all(x in S for x in (0x540628, 0x540688, 0x540678)):
|
||
continue
|
||
s = S[0x540628][:1025].astype(np.float64)
|
||
t = S[0x540688][:1025].astype(np.float64)
|
||
c = S[0x540678][:1025].astype(np.float64)
|
||
ok = (t > 1e-30) & (c > 1e-30) & np.isfinite(s)
|
||
if ok.sum() < 50:
|
||
continue
|
||
lt = np.log(t[ok])
|
||
lc = np.log(c[ok])
|
||
sel = np.abs(lt) > 0.05
|
||
if sel.sum() < min_bins:
|
||
continue
|
||
g = float(np.sum(lt[sel] * lc[sel]) / np.sum(lt[sel] ** 2))
|
||
rms = float(np.sqrt(np.mean((lc[sel] - g * lt[sel]) ** 2)))
|
||
depth = float(-lc.min())
|
||
key = 'gamma' if abs(g - GAMMA) < 0.01 else \
|
||
('identity' if abs(g - 1.0) < 1e-4 else None)
|
||
if key is None or rms > 1e-4:
|
||
continue
|
||
cand = (depth, f, s, t, c, g, rms)
|
||
if classes[key] is None or depth > classes[key][0]:
|
||
classes[key] = cand
|
||
return classes
|
||
|
||
|
||
# ------------------------------------------------------- валидация -------
|
||
def validate_scr(sim_scr, cap_scr, tol_db=0.05):
|
||
"""rms в дБ между симулированным и захваченным scr."""
|
||
m = np.abs(cap_scr) > 0.02
|
||
err = (sim_scr[m] - cap_scr[m]) * (20 / np.log(10))
|
||
return float(np.sqrt(np.mean(err ** 2))), int(m.sum())
|
||
|
||
|
||
def win_periodic_hann(N):
|
||
return 0.5 * (1.0 - np.cos(2.0 * np.pi * np.arange(N) / N))
|
||
|
||
|
||
# ------------------------------------------------ FIR-цепь (24mm9) --------
|
||
NFRAME = 4096 # n=[ctx+0x540534]
|
||
NBINS_FIR = NFRAME // 2 + 1
|
||
Q_EXP = 0.80 # скаляр аргумента EXP; источник в 1803831c0 (ОТКРЫТО)
|
||
|
||
|
||
def winfreq_fall():
|
||
"""WINfreq@[ctx+0x540658]: периодический Hann(4096), падающая половина."""
|
||
return win_periodic_hann(NFRAME)[NFRAME // 2:]
|
||
|
||
|
||
def fir_kernel(scr, q=Q_EXP):
|
||
"""Полная FIR-цепь (BLOCKMAP 24mm9): min-phase кепстральный сэндвич.
|
||
|
||
scr(2049) → pack(re=scr,im=0) → FIR[n]=0 (Найквост)
|
||
→ inv-RFFT → fold(y[1..2047]*=2.0 @1824c41e0; y[2049..4095]=0)
|
||
→ fwd-RFFT → комплексная EXP (1803831c0, аргумент ×q)
|
||
→ inv-RFFT → ×падающий Hann → ноль хвоста → fwd-RFFT
|
||
→ FIR[0]=1, FIR[1]=0. Возвращает |F| (2049).
|
||
"""
|
||
h = np.asarray(scr, dtype=np.complex128).copy()
|
||
h[-1] = 0.0
|
||
y = np.fft.irfft(h, n=NFRAME)
|
||
y[1:NFRAME // 2] *= 2.0
|
||
y[NFRAME // 2 + 1:] = 0.0
|
||
w = np.fft.irfft(np.exp(q * np.fft.rfft(y, n=NFRAME)), n=NFRAME)
|
||
w[:NFRAME // 2] *= winfreq_fall()
|
||
w[NFRAME // 2:] = 0.0
|
||
F = np.abs(np.fft.rfft(w, n=NFRAME))
|
||
F[0] = 1.0
|
||
return F
|
||
|
||
|
||
def mask_from_frame(S, q=Q_EXP):
|
||
"""mask_sim из слотов кадра: cur ≈ trk · |F(scr)| (df0 complex-mul)."""
|
||
scr = S[0x540628][:NBINS_FIR].astype(np.float64)
|
||
trk = S[0x540688][:NBINS_FIR].astype(np.float64)
|
||
return trk * fir_kernel(scr, q)
|
||
|
||
|
||
def validate_mask_stage(ds, q=Q_EXP, cap=60):
|
||
"""Валидация масочной ветви на чистых γ-кадрах (24mm9-протокол).
|
||
|
||
Отбор: |γ_fit−1.760561|<5e-4 и fit-rms<1e-5 (жёстче pick_clean_frame).
|
||
Критерий: rms по ВСЕМ 2049 бинам < 0.05 дБ (структурная фаза).
|
||
"""
|
||
import glob
|
||
rmss, gpred = [], []
|
||
for f in sorted(glob.glob(os.path.join(ds, 'ph*.npz'))):
|
||
try:
|
||
d = np.load(f)
|
||
except Exception:
|
||
continue
|
||
if '0x540628' not in d:
|
||
continue
|
||
scr = d['0x540628'][:NBINS_FIR].astype(np.float64)
|
||
trk = d['0x540688'][:NBINS_FIR].astype(np.float64)
|
||
cur = d['0x540678'][:NBINS_FIR].astype(np.float64)
|
||
ok = (trk > 1e-30) & (cur > 1e-30) & np.isfinite(scr)
|
||
if ok.sum() < 50:
|
||
continue
|
||
lt, lc = np.log(trk[ok]), np.log(cur[ok])
|
||
sel = np.abs(lt) > 0.05
|
||
if sel.sum() < 8:
|
||
continue
|
||
g = float(np.sum(lt[sel] * lc[sel]) / np.sum(lt[sel] ** 2))
|
||
frms = float(np.sqrt(np.mean((lc[sel] - g * lt[sel]) ** 2)))
|
||
if not (abs(g - GAMMA) < 5e-4 and frms < 1e-5):
|
||
continue
|
||
F = fir_kernel(scr, q)
|
||
lf = np.log(F[sel])
|
||
lt_s = np.log(trk[sel])
|
||
sF = float(np.sum(lf * lt_s) / np.sum(lt_s ** 2))
|
||
gpred.append(1.0 + sF)
|
||
m = trk * F
|
||
mm = (cur > 1e-12) & (m > 1e-12)
|
||
e = (np.log(m[mm]) - np.log(cur[mm])) * 20 / np.log(10)
|
||
rmss.append(float(np.sqrt(np.mean(e ** 2))))
|
||
if len(rmss) >= cap:
|
||
break
|
||
if not rmss:
|
||
print('нет ультрачистых кадров в', ds)
|
||
return
|
||
rmss = np.array(rmss)
|
||
print('кадров=%d | rms медиана=%.4f дБ p90=%.4f max=%.4f | '
|
||
'gamma_pred(1+s_F)=%.6f' %
|
||
(len(rmss), np.median(rmss), np.percentile(rmss, 90), rmss.max(),
|
||
float(np.median(gpred))))
|
||
|
||
|
||
def main():
|
||
import sys
|
||
if len(sys.argv) > 1 and sys.argv[1] == '--mask':
|
||
validate_mask_stage(sys.argv[2] if len(sys.argv) > 2
|
||
else '/tmp/opencode/sc_multi4b')
|
||
return
|
||
ds = sys.argv[1] if len(sys.argv) > 1 else '/tmp/opencode/sc_multi6'
|
||
tract = sys.argv[2] if len(sys.argv) > 2 else '/tmp/opencode/tract_multi6.txt'
|
||
|
||
classes = pick_clean_frame(ds)
|
||
ph_g = classes['gamma']
|
||
ph_i = classes['identity']
|
||
if not ph_g and not ph_i:
|
||
print('нет чистых кадров в', ds)
|
||
return
|
||
for lbl, best in (('γ-фаза', ph_g), ('identity', ph_i)):
|
||
if not best:
|
||
continue
|
||
_, f, scr, trk, cur, gamma_fit, grms = best
|
||
n = len(scr)
|
||
cut_meas = -20 / np.log(10) * np.log(np.maximum(cur, 1e-30))
|
||
g_use = gamma_fit
|
||
cut_sim = g_use * (-scr) * 20 / np.log(10)
|
||
e = cut_sim - cut_meas
|
||
sel = np.abs(cut_meas) > 0.1
|
||
rms_db = float(np.sqrt(np.mean(e[sel] ** 2))) if sel.any() else 0.0
|
||
print(f'{lbl}: {os.path.basename(f)} γ={gamma_fit:.6f} (rms {grms:.1e}) '
|
||
f'закон: rms={rms_db:.4f} дБ / {int(sel.sum())} бинов')
|
||
|
||
|
||
if __name__ == '__main__':
|
||
main()
|