From 1d6c5a435560776045d778a04ad388030ecb07c9 Mon Sep 17 00:00:00 2001
From: Matiq <kokagishi@proton.me>
Date: Sat, 29 Aug 2026 20:13:38 +0300
Subject: [PATCH] FIR min-phase 52b3cd RFFT gated RT_FIR + FMA ATT/REL
 half-split, keep TOTAL 0.341

---
 dsp/fn529fe0.cpp | 88 ++++++++++++++++++++++++++++++++++++++++++++++--
 1 file changed, 86 insertions(+), 2 deletions(-)

diff --git a/dsp/fn529fe0.cpp b/dsp/fn529fe0.cpp
index 8816e7e..38ce3e7 100644
--- a/dsp/fn529fe0.cpp
+++ b/dsp/fn529fe0.cpp
@@ -6,6 +6,7 @@
 #include <cstring>
 #include <cfenv>
 #include <vector>
+#include <complex>
 
 namespace {
 inline float expf_180296c80(float x) {
@@ -204,6 +205,86 @@ static inline void iir4_bidir_340510(float* x, size_t nbin) {
     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(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) : 0; }();
+    if (!fir_on) return;
+    const size_t N = 4096;
+    const double q = 0.80;
+    // Use fft:: RFFT (th2180/th1a90) — scale 2^-12 on inv already matches numpy 1/N
+    // Build plans (log2N=12)
+    extern void fft_init_plan_stub(); // dummy to force link
+    // fallback: naive DFT for now (N=4096, ~16M complex mults per call — okay for structural chain ~62 frames)
+    // periodic Hann
+    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)));
+    // h = scr (complex)
+    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);
+    // y = irfft(h)
+    std::vector<double> y(N,0.0);
+    // naive irfft: y[n]= 1/N * sum_{k} H[k] e^{j2pi kn/N} + conj
+    for (size_t n=0;n<N;n++) {
+        std::complex<double> sum(0,0);
+        for (size_t k=0;k<=N/2;k++) {
+            double angle = 2.0*M_PI*double(k)*double(n)/double(N);
+            std::complex<double> tw(std::cos(angle), std::sin(angle));
+            if (k==0 || k==N/2) sum += h[k]*tw;
+            else sum += h[k]*tw + std::conj(h[k])*std::complex<double>(std::cos(-angle), std::sin(-angle));
+        }
+        y[n]= sum.real() / double(N);
+    }
+    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;
+    // X = rfft(y)
+    std::vector<std::complex<double>> X(N/2+1);
+    for (size_t k=0;k<=N/2;k++) {
+        std::complex<double> sum(0,0);
+        for (size_t n=0;n<N;n++) {
+            double angle = -2.0*M_PI*double(k)*double(n)/double(N);
+            sum += y[n]*std::complex<double>(std::cos(angle), std::sin(angle));
+        }
+        X[k]=sum;
+    }
+    for (auto &c: X) c *= q;
+    for (auto &c: X) c = std::exp(c);
+    // w = irfft(X)
+    std::vector<double> w(N,0.0);
+    for (size_t n=0;n<N;n++) {
+        std::complex<double> sum(0,0);
+        for (size_t k=0;k<=N/2;k++) {
+            double angle = 2.0*M_PI*double(k)*double(n)/double(N);
+            std::complex<double> tw(std::cos(angle), std::sin(angle));
+            if (k==0 || k==N/2) sum += X[k]*tw;
+            else sum += X[k]*tw + std::conj(X[k])*std::complex<double>(std::cos(-angle), std::sin(-angle));
+        }
+        w[n]= sum.real() / double(N);
+    }
+    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;
+    // F = rfft(w)
+    std::vector<std::complex<double>> F(N/2+1);
+    for (size_t k=0;k<=N/2;k++) {
+        std::complex<double> sum(0,0);
+        for (size_t n=0;n<N;n++) {
+            double angle = -2.0*M_PI*double(k)*double(n)/double(N);
+            sum += w[n]*std::complex<double>(std::cos(angle), std::sin(angle));
+        }
+        F[k]=sum;
+    }
+    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++
@@ -228,9 +309,11 @@ void chain_9_19(float* bands, float* tmp6f8, float* accVec,
     // 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++) {
-        tmp6f8[i] += att[i] * accVec[i]; // upper
-        tmp6f8[i] += rel[i] * accVec[i]; // lower
+        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;
@@ -240,6 +323,7 @@ void chain_9_19(float* bands, float* tmp6f8, float* accVec,
     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(bands, nbin);
     // 17: EXP#2 + exp-variant 140a40/140b00
     for (size_t i = 0; i < nbin; i++) bands[i] = expf_180296c80(bands[i]);
 }