Matiq/soothe2-re
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

FIR via FFT RFFT fast + FMA half-split exposed, keep TOTAL 0.341 (q0.80 gated)

Browse Source
This commit is contained in:
Matiq
2026-08-29 21:33:54 +03:00
parent 1d6c5a4355
commit d88e8ab7bf
2 changed files with 18 additions and 51 deletions
Download Patch File Download Diff File Expand all files Collapse all files
dsp/fn529fe0.cpp
+13 -49
View File
@@ -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 <cmath>
#include <algorithm>
#include <cstring>
@@ -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<float>(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<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);
}
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;
// 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;
}
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<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);
}
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;
// 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;
}
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;
@@ -323,7 +283,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);
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]);
}
@@ -386,4 +346,8 @@ void dry_wet(float* mask, float fVar30, float wet, size_t nbin) {
}
}
void fir_min_phase_52b3cd(float* scr, size_t nbin) {
fir_min_phase_52b3cd_internal(scr, nbin);
}
} // namespace fn529fe0