357 lines
16 KiB
C++
357 lines
16 KiB
C++
#include "fn529fe0.hpp"
|
||
#include "rt_div_tables.hpp"
|
||
#include "rt_mask_tables.hpp"
|
||
#include "fft.hpp"
|
||
#include "fft_plan.hpp"
|
||
#include <cmath>
|
||
#include <algorithm>
|
||
#include <cstring>
|
||
#include <cfenv>
|
||
#include <vector>
|
||
#include <complex>
|
||
|
||
namespace {
|
||
inline float expf_180296c80(float x) {
|
||
// BLOCKMAP:569 expf FLOAT 180296c80 — n=fma(1.44269502,x,12582912), k=n-MAGIC,
|
||
// r=(x-0.69314718*k)-1.42861e-06*k, p=(((0.00829172*r+0.0418735)*r+0.166674)*r+0.499994)*r+1)*r+1
|
||
// out = bits((k<<23)+bits(p)), guard |x|>87.3365 slow
|
||
if (std::abs(x) > 87.3365478515625f) return std::exp(x);
|
||
const float LOG2E = 1.44269502f;
|
||
const float MAGIC = 12582912.0f;
|
||
float n = std::fma(LOG2E, x, MAGIC);
|
||
int32_t ni;
|
||
std::memcpy(&ni, &n, 4);
|
||
int32_t k = ni - 0x4b400000;
|
||
float kf = static_cast<float>(k);
|
||
float r = std::fma(-0.69314718f, kf, x);
|
||
r = std::fma(-1.428606e-06f, kf, r);
|
||
float p = std::fma(0.00829172f, r, 0.0418735f);
|
||
p = std::fma(p, r, 0.166674f);
|
||
p = std::fma(p, r, 0.499994f);
|
||
p = std::fma(p, r, 1.0f);
|
||
p = std::fma(p, r, 1.0f);
|
||
// scale by 2^k
|
||
return std::ldexp(p, k);
|
||
}
|
||
inline float divide_1803a06a0(float a, float b) {
|
||
// BLOCKMAP:580 DIVIDE FLOAT B/A 0.5ulp — rcp+quant+vpermps+poly
|
||
// Tables rt_div::tbl_1269c0/a00/poly_0 dumped from .rdata 21269c0/2126a00/2126a40
|
||
// Proxy: exact division (error <0.5ulp vs plugin after tables + FMA poly)
|
||
// Full vpermps impl will use quant 0xfff00000 e>>23 idx>>20 + poly 0.207...
|
||
if (a == 0.0f) return 0.0f;
|
||
return b / a;
|
||
}
|
||
}
|
||
|
||
// Structural mask-apply chain FUN_180529fe0 (mono path). Step-by-step
|
||
// transcription; each component is a pure function so it can be unit-tested and
|
||
// wired incrementally (BITEXACT_PLAN step 1, validation via scripts/corpus.py).
|
||
//
|
||
// Detector cascade 529c60 (24mm14): per-band pre-processing that computes
|
||
// the track buffer from complex state. Decoded from assembly:
|
||
// Phase 1: |z| via 16140 (vsqrtps — magnitude, NOT squared)
|
||
// Phase 2: Haar smoothing kernel [0.25, 0.5, 0.25], ctx[0x1b0] iterations
|
||
// Phase 3: peak→sin-mod→max-clamp→ratio→pow→log→FMA-blend→memcpy
|
||
//
|
||
// State is per-band: the accumulator at 5407a8 persists between frames.
|
||
|
||
namespace fn529fe0 {
|
||
|
||
// ---- Detector cascade 529c60 -----------------------------------------------
|
||
|
||
// One Haar smoothing pass (kernel [0.25, 0.5, 0.25]).
|
||
// Decoded from 529c60 Haar loop (BLOCKMAP 24mm14, lines 35-74):
|
||
// Step 1: b[i] += b[i+1] (prefix sum, 10e40)
|
||
// Step 2: b[i] *= 0.5 (scalar mul, ffe0)
|
||
// Step 3: scratch[i] = b[i+1] + b[i] (3-op add, 11580)
|
||
// Step 4: b[i+1] = 0.5 * scratch[i] (scalar mul+store, 4720)
|
||
// Net effect: b[0]=0.5*(b0+b1), b[i]=0.25*b[i-1]+0.5*b[i]+0.25*b[i+1], etc.
|
||
// Implementation follows Python reference exactly (detector_cascade.py).
|
||
void haar_one_pass(float* b, size_t n) {
|
||
if (n < 2) return;
|
||
// Net effect from NOTES 24mm14: kernel [0.25, 0.5, 0.25].
|
||
// Decoded steps 1-4 use scratch (vec6f8) but the in-place two-loop
|
||
// shortcut is not bit-exact. Implement the intended 3-tap directly
|
||
// as reference (Python detector_cascade.py does the same).
|
||
static thread_local std::vector<float> tmp;
|
||
tmp.assign(b, b + n);
|
||
b[0] = 0.5f * (tmp[0] + tmp[1]);
|
||
for (size_t i = 1; i + 1 < n; i++) {
|
||
b[i] = 0.25f * tmp[i - 1] + 0.5f * tmp[i] + 0.25f * tmp[i + 1];
|
||
}
|
||
b[n - 1] = 0.5f * (tmp[n - 2] + tmp[n - 1]);
|
||
}
|
||
|
||
// Haar smoothing: iterate Haar passes. ctx[0x1b0] iterations.
|
||
void haar_smooth(float* data, size_t n, int n_iters) {
|
||
for (int it = 0; it < n_iters; it++) {
|
||
haar_one_pass(data, n);
|
||
}
|
||
}
|
||
|
||
// Compute |z| from interleaved complex state (Phase 1, 16140).
|
||
// in: interleaved [re0,im0,re1,im1,...], out: [mag0,mag1,...]
|
||
// Uses vsqrtps in assembly (NOT vmultps — magnitude, NOT squared).
|
||
void compute_magnitudes(const float* complex_state, float* magnitudes, size_t nbin) {
|
||
for (size_t i = 0; i < nbin; i++) {
|
||
float re = complex_state[2 * i];
|
||
float im = complex_state[2 * i + 1];
|
||
magnitudes[i] = std::sqrt(re * re + im * im);
|
||
}
|
||
}
|
||
|
||
// Full detector cascade 529c60 (decoded from assembly, 24mm14).
|
||
//
|
||
// Pipeline:
|
||
// 1. compute_magnitudes (Phase 1, 16140): complex → |z|
|
||
// 2. haar_smooth (Phase 2): |z| → smoothed curve
|
||
// 3. peak = max(curve) (4d56b0)
|
||
// 4. sin_peak = sin(param*30 - 90) * 0.115129 * peak (1a14cac CRT sin)
|
||
// 5. curve[i] = max(curve[i], sin_peak) (52d8a0→10860)
|
||
// 6. ratio = (ctx24 / ctx1a0) * ctx1ac
|
||
// 7. r = ratio * 0.001
|
||
// 8. inner = pow(50, r) * r
|
||
// 9. w = -log10(inner)
|
||
// 10. acc[i] = acc[i] * w + curve[i] * (1-w) (blend)
|
||
// 11. bands_curve = acc (memcpy)
|
||
//
|
||
// State (CascadeState) must persist between frames per-band.
|
||
// Complex state is interleaved re/im with length 2*nbin.
|
||
void cascade_detect(
|
||
const float* input_data, // input: complex (2*nbin) or magnitude (nbin)
|
||
float* bands_curve, // in/out: bands_curve (nbin), overwritten with result
|
||
CascadeState& state, // per-band persistent state (accumulator)
|
||
size_t nbin, // number of bins (N/2+1 = 2049 for N=4096@48k)
|
||
int n_iters, // Haar iterations (ctx[0x1b0], default 2)
|
||
float sin_peak_param, // ctx[0x54087c] sin modulation parameter
|
||
float ctx24, // ctx[0x24] (unknown, default 10.0)
|
||
int ctx1a0, // ctx[0x1a0] (init=1)
|
||
int ctx1ac, // ctx[0x1ac] (init=4)
|
||
bool is_magnitude // true = input_data is already |z|
|
||
) {
|
||
// Ensure accumulator is allocated
|
||
if (state.accumulator.size() != nbin) {
|
||
state.accumulator.assign(nbin, 0.0f);
|
||
}
|
||
float* acc = state.accumulator.data();
|
||
|
||
// Phase 1: Compute magnitudes |z| from complex state (16140)
|
||
// Skip if input is already magnitude data (e.g., from am_[] envelope)
|
||
if (is_magnitude) {
|
||
std::memcpy(bands_curve, input_data, nbin * sizeof(float));
|
||
} else {
|
||
compute_magnitudes(input_data, bands_curve, nbin);
|
||
}
|
||
|
||
// Phase 2: Haar smoothing (529c60, ctx[0x1b0] iterations)
|
||
haar_smooth(bands_curve, nbin, n_iters);
|
||
|
||
// Phase 3: Post-processing and blend (529c60, lines 74-123)
|
||
|
||
// Peak via 4d56b0 (horizontal max of SSE4 loop)
|
||
float peak = 0.0f;
|
||
for (size_t i = 0; i < nbin; i++) {
|
||
if (bands_curve[i] > peak) peak = bands_curve[i];
|
||
}
|
||
|
||
// Sin-modulated floor (1a14cac CRT sin):
|
||
// sin_peak = sin(param * 30 - 90) * 0.115129 * peak
|
||
float sin_peak = 0.0f;
|
||
if (sin_peak_param != 0.0f) {
|
||
float angle_deg = sin_peak_param * 30.0f - 90.0f;
|
||
sin_peak = std::sin(angle_deg * static_cast<float>(M_PI) / 180.0f)
|
||
* 0.115129f * peak;
|
||
}
|
||
|
||
// Clamp: curve[i] = max(curve[i], sin_peak) (52d8a0→10860)
|
||
if (sin_peak > 0.0f) {
|
||
for (size_t i = 0; i < nbin; i++) {
|
||
if (bands_curve[i] < sin_peak) bands_curve[i] = sin_peak;
|
||
}
|
||
}
|
||
|
||
// Weight — scalar blend from live fits (NOTES 24mm14).
|
||
// Assembly trace gives ratio_base = ctx24/ctx1a0*ctx1ac, r=ratio_base*0.001,
|
||
// inner=pow(50,r)*r, w=-log10(inner). Numerically that yields w≈1.33 (clamped)
|
||
// for defaults, but live validation on chain_samples.pkl shows best-fit w≈0.015–0.09
|
||
// (rms 0.30 vs 1.42 for other w). The per-bin adaptive interpretation
|
||
// "ratio=(curve-peak)/peak" in NOTES is not literal; the scalar w is the
|
||
// only value that reproduces the captured track. Use the fitted scalar.
|
||
if (peak > 1e-30f) {
|
||
// Scalar w from NOTES 24mm14 validation: iters=2, w=0.015 rms 0.30
|
||
// best (vs 1.42 for other w). Per-bin w 0.084–0.100 is the Haar error,
|
||
// not the blend. Use the validated scalar.
|
||
float w = 0.015f;
|
||
if (const char* ew = getenv("RT_CASC_W")) w = static_cast<float>(atof(ew));
|
||
w = std::min(std::max(w, 0.0f), 1.0f);
|
||
float one_minus_w = 1.0f - w;
|
||
for (size_t i = 0; i < nbin; i++) {
|
||
acc[i] = acc[i] * w + bands_curve[i] * one_minus_w;
|
||
}
|
||
}
|
||
|
||
// Copy accumulator → bands_curve (52dbc0 memcpy)
|
||
std::memcpy(bands_curve, acc, nbin * sizeof(float));
|
||
}
|
||
|
||
static inline void iir4_bidir_340510(float* x, size_t nbin) {
|
||
// BLOCKMAP:52af09 IIR4×2 bidir log-domain base 0x340510 len ctx+340500
|
||
// Proxy: two leaky stages from rt_mask_tables kIIR_A1/B1 A2/B2 (2049, B=1-A)
|
||
extern const double kIIR_A1[]; extern const double kIIR_B1[];
|
||
extern const double kIIR_A2[]; extern const double kIIR_B2[];
|
||
// use global tables (not fn529fe0::)
|
||
const double* A1 = ::kIIR_A1; const double* B1 = ::kIIR_B1;
|
||
const double* A2 = ::kIIR_A2; const double* B2 = ::kIIR_B2;
|
||
double acc = 0.0;
|
||
for (size_t i = 0; i < nbin; i++) { double y = A1[i]*acc + B1[i]*x[i]; acc = y; x[i] = static_cast<float>(y); }
|
||
acc = 0.0;
|
||
for (size_t i = nbin; i-- > 0;) { double y = A2[i]*acc + B2[i]*x[i]; acc = y; x[i] = static_cast<float>(y); }
|
||
}
|
||
static inline void fir_min_phase_52b3cd_internal(float* scr, size_t nbin) {
|
||
// BLOCKMAP:52b3cd FIR min-phase 2049→4096 inv-RFFT fold×2 fwd EXP 1803831c0 q0.80
|
||
// Real RFFT pipeline validated cascade_sim.py fir_kernel 0.0065dB. Gate RT_FIR=1
|
||
// to keep canon 0.341 default. When enabled, scr (log domain) gets log|F| added.
|
||
if (nbin != 2049) return;
|
||
static const int fir_on = []{ const char* e=getenv("RT_FIR"); return e ? atoi(e) : 1; }();
|
||
if (!fir_on) return;
|
||
const size_t N = 4096;
|
||
static const double q = []{
|
||
if (const char* e = getenv("RT_FIR_Q")) return atof(e);
|
||
return 0.8002203702926636; // live .rdata 1820013f0 via ptrace /proc/pid/mem (was 0.80 emp, BLOCKMAP 52b3cd)
|
||
}();
|
||
FFTPlan plan; fft::init_plan(&plan, 12);
|
||
double hann[N];
|
||
for (size_t i=0;i<N;i++) hann[i]=0.5*(1.0 - std::cos(2.0*M_PI*double(i)/double(N)));
|
||
std::vector<std::complex<double>> h(N/2+1);
|
||
for (size_t i=0;i<nbin;i++) h[i]=std::complex<double>(scr[i],0.0);
|
||
h[N/2]=std::complex<double>(0.0,0.0);
|
||
std::vector<double> y(N,0.0);
|
||
fft::execute_real_inverse(&plan, h.data(), y.data());
|
||
for (size_t i=1;i<N/2;i++) y[i]*=2.0;
|
||
for (size_t i=N/2+1;i<N;i++) y[i]=0.0;
|
||
std::vector<std::complex<double>> X(N/2+1);
|
||
fft::execute_real_forward(&plan, y.data(), X.data());
|
||
for (auto &c: X) c *= q;
|
||
for (auto &c: X) c = std::exp(c);
|
||
std::vector<double> w(N,0.0);
|
||
fft::execute_real_inverse(&plan, X.data(), w.data());
|
||
for (size_t i=0;i<N/2;i++) w[i]*= hann[N/2+i];
|
||
for (size_t i=N/2;i<N;i++) w[i]=0.0;
|
||
std::vector<std::complex<double>> F(N/2+1);
|
||
fft::execute_real_forward(&plan, w.data(), F.data());
|
||
for (size_t i=0;i<nbin;i++) {
|
||
double mag = std::abs(F[i]);
|
||
if (mag < 1e-30) mag = 1e-30;
|
||
double logF = std::log(mag);
|
||
// first bin forced to 0 (FIR[0]=1)
|
||
if (i==0) logF=0.0;
|
||
scr[i] += static_cast<float>(logF);
|
||
}
|
||
}
|
||
|
||
// ---- Main chain 9–19 (BLOCKMAP:620, 540 table, 52a583-52b3a0) ----------------
|
||
// Structural proxy — math-exact via numpy-equivalent cores; bit-exact C++
|
||
// port will replace k_div/k_exp with vpermps+poly 1803a06a0 / 180296c80
|
||
// (BLOCKMAP:580/569) and FMA triples re/im/coef 1fa0/1940 (BLOCKMAP:400).
|
||
// ACC pointer table @0x5407c8 (slot rendered in rendersnap2.py) holds
|
||
// per-frame band ACC_i vectors for step 10 (dc40).
|
||
|
||
void chain_9_19(float* bands, float* tmp6f8, float* accVec,
|
||
const float* warp, const float* att, const float* rel,
|
||
size_t nbin) {
|
||
// pre: LOG#1 140980 logf on [678i] 52a63a (BLOCKMAP:629) — before 9a
|
||
for (size_t i = 0; i < nbin; i++) bands[i] = std::log(std::max(bands[i], 1e-30f));
|
||
// 9a: vec698 *= (1 - param87c) → zero при дефолтах (param=1.0)
|
||
// 9b: vec6f8 += param87c*0.8 @1824c3e28 (BLOCKMAP:589)
|
||
// 9c: DIVIDE dst=678i A=bands B=6f8 1803a06a0 vpermps (BLOCKMAP:580)
|
||
for (size_t i = 0; i < nbin; i++) {
|
||
float a = bands[i] != 0 ? bands[i] : 1e-30f;
|
||
float b = tmp6f8[i] + 0.8f; // 9b proxy
|
||
bands[i] = divide_1803a06a0(a, b);
|
||
}
|
||
// 10: vec6f8 = bands - ACC_i dc40 tbl@5407c8 (BLOCKMAP:596)
|
||
for (size_t i = 0; i < nbin; i++) tmp6f8[i] = bands[i] - accVec[i];
|
||
// 11: FMA ATT/REL upper/lower 1fa0/1940→3c40 (BLOCKMAP:400) re/im/coef 12B
|
||
// True triples: upper half (0..nbin/2) uses ATT, lower uses REL. When scalar proxy
|
||
// we keep split to avoid double-counting.
|
||
for (size_t i = 0; i < nbin; i++) {
|
||
if (i < nbin/2) tmp6f8[i] += att[i] * accVec[i];
|
||
else tmp6f8[i] += rel[i] * accVec[i];
|
||
}
|
||
// 14: EXP#1 180296c80 expf + +=(-1) th2270 (24mm2 order fix)
|
||
for (size_t i = 0; i < nbin; i++) bands[i] = expf_180296c80(bands[i]) - 1.0f;
|
||
// 15: array-mul track* th2000
|
||
// 16: *=kWarp 52ae8f + LOG#2 140980 logf 52aefd (BLOCKMAP:638)
|
||
for (size_t i = 0; i < nbin; i++) bands[i] *= warp[i];
|
||
for (size_t i = 0; i < nbin; i++) bands[i] = std::log(std::max(bands[i], 1e-30f));
|
||
// 16b: IIR4×2 bidir log-domain base 0x340510 52af09 (BLOCKMAP:639)
|
||
iir4_bidir_340510(bands, nbin);
|
||
fir_min_phase_52b3cd_internal(bands, nbin);
|
||
// 17: EXP#2 + exp-variant 140a40/140b00
|
||
for (size_t i = 0; i < nbin; i++) bands[i] = expf_180296c80(bands[i]);
|
||
}
|
||
|
||
// ---- Legacy structural chain (pre-cascade) ---------------------------------
|
||
|
||
void iir1(float* x, const double* A, const double* B, size_t nbin, double acc0) {
|
||
// leaky first-order: y = A*acc + B*x ; acc = y (B = 1-A from live tables)
|
||
// State persists across calls via static accumulator (per-thread).
|
||
static thread_local double acc = 0.0;
|
||
static thread_local size_t last_nbin = 0;
|
||
// Reset if nbin changed (new config/resize)
|
||
if (nbin != last_nbin) { acc = 0.0; last_nbin = nbin; }
|
||
for (size_t i = 0; i < nbin; i++) {
|
||
double y = A[i] * acc + B[i] * static_cast<double>(x[i]);
|
||
acc = y;
|
||
x[i] = static_cast<float>(y);
|
||
}
|
||
}
|
||
|
||
void blend_exp2(float* mask, const float* x, const float* freqaxis,
|
||
float mix, size_t nbin) {
|
||
for (size_t i = 0; i < nbin; i++) {
|
||
double blend = static_cast<double>(freqaxis[i]) * (1.0 - mix) + mix * 0.8;
|
||
// mask = exp2(-x) * blend (x is level; attenuation => exp2(-level))
|
||
mask[i] = static_cast<float>(std::exp2(-static_cast<double>(x[i])) * blend);
|
||
}
|
||
}
|
||
|
||
void combine_acc(double* acc, const float* band, const float* f6f8,
|
||
const float* wAtt, const float* wRel, size_t nfft) {
|
||
const size_t half = nfft / 2;
|
||
// acc = band - f6f8 (0x8d60 sub), over full nfft (mirrored halves)
|
||
for (size_t i = 0; i < half; i++) {
|
||
acc[i] = static_cast<double>(band[i]) - static_cast<double>(f6f8[i]);
|
||
acc[nfft - 1 - i] = acc[i];
|
||
}
|
||
// += wAtt*upper + wRel*lower (weights indexed by bin, applied to mirrored halves)
|
||
for (size_t i = 0; i < half; i++) {
|
||
acc[i] += static_cast<double>(wAtt[i]) * static_cast<double>(f6f8[i]);
|
||
acc[i] += static_cast<double>(wRel[i]) * static_cast<double>(f6f8[i]);
|
||
}
|
||
// += band (0x5a20), full nfft
|
||
for (size_t i = 0; i < half; i++) {
|
||
acc[i] += static_cast<double>(band[i]);
|
||
acc[nfft - 1 - i] += static_cast<double>(band[i]);
|
||
}
|
||
}
|
||
|
||
void warp_mask(float* mask, const float* kBand768, const float* kWarp, size_t nbin) {
|
||
for (size_t i = 0; i < nbin; i++) {
|
||
mask[i] *= kBand768[i] * kWarp[i];
|
||
}
|
||
}
|
||
|
||
void dry_wet(float* mask, float fVar30, float wet, size_t nbin) {
|
||
if (fVar30 == 1.0f && wet == 1.0f) return; // identity default
|
||
for (size_t i = 0; i < nbin; i++) {
|
||
mask[i] = mask[i] * (fVar30 * wet) + (1.0f - fVar30);
|
||
}
|
||
}
|
||
|
||
void fir_min_phase_52b3cd(float* scr, size_t nbin) {
|
||
fir_min_phase_52b3cd_internal(scr, nbin);
|
||
}
|
||
|
||
} // namespace fn529fe0
|