From d88e8ab7bfa4cfebf1eb5dcc48f9877647fa900a Mon Sep 17 00:00:00 2001 From: Matiq Date: Sat, 29 Aug 2026 21:33:54 +0300 Subject: [PATCH] FIR via FFT RFFT fast + FMA half-split exposed, keep TOTAL 0.341 (q0.80 gated) --- dsp/fn529fe0.cpp | 62 ++++++++++-------------------------------------- dsp/fn529fe0.hpp | 7 ++++-- 2 files changed, 18 insertions(+), 51 deletions(-) diff --git a/dsp/fn529fe0.cpp b/dsp/fn529fe0.cpp index 38ce3e7..c13329d 100644 --- a/dsp/fn529fe0.cpp +++ b/dsp/fn529fe0.cpp @@ -1,6 +1,8 @@ #include "fn529fe0.hpp" #include "rt_div_tables.hpp" #include "rt_mask_tables.hpp" +#include "fft.hpp" +#include "fft_plan.hpp" #include #include #include @@ -205,7 +207,7 @@ 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(y); } } -static inline void fir_min_phase_52b3cd(float* scr, size_t nbin) { +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. @@ -214,68 +216,26 @@ static inline void fir_min_phase_52b3cd(float* scr, size_t nbin) { 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 + FFTPlan plan; fft::init_plan(&plan, 12); double hann[N]; for (size_t i=0;i> h(N/2+1); for (size_t i=0;i(scr[i],0.0); h[N/2]=std::complex(0.0,0.0); - // y = irfft(h) std::vector 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 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 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(std::cos(-angle), std::sin(-angle)); - } - y[n]= sum.real() / double(N); - } + fft::execute_real_inverse(&plan, h.data(), y.data()); for (size_t i=1;i> X(N/2+1); - for (size_t k=0;k<=N/2;k++) { - std::complex sum(0,0); - for (size_t n=0;n(std::cos(angle), std::sin(angle)); - } - X[k]=sum; - } + fft::execute_real_forward(&plan, y.data(), X.data()); for (auto &c: X) c *= q; for (auto &c: X) c = std::exp(c); - // w = irfft(X) std::vector w(N,0.0); - for (size_t n=0;n 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 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(std::cos(-angle), std::sin(-angle)); - } - w[n]= sum.real() / double(N); - } + fft::execute_real_inverse(&plan, X.data(), w.data()); for (size_t i=0;i> F(N/2+1); - for (size_t k=0;k<=N/2;k++) { - std::complex sum(0,0); - for (size_t n=0;n(std::cos(angle), std::sin(angle)); - } - F[k]=sum; - } + fft::execute_real_forward(&plan, w.data(), F.data()); for (size_t i=0;i