From 031f7d21afd74ecf54a6a6ce63e01478908d30b0 Mon Sep 17 00:00:00 2001
From: Matiq <kokagishi@proton.me>
Date: Sat, 29 Aug 2026 15:10:24 +0300
Subject: [PATCH] k-mapping 2: RBJ peaking EQ for W_eq, keep per-fc, TOTAL
 0.341

---
 dsp/framed_model.cpp | 30 ++++++++++++++++++++++++------
 1 file changed, 24 insertions(+), 6 deletions(-)

diff --git a/dsp/framed_model.cpp b/dsp/framed_model.cpp
index 9d9f166..d9b07b1 100644
--- a/dsp/framed_model.cpp
+++ b/dsp/framed_model.cpp
@@ -696,12 +696,30 @@ void FramedDetector::setParams(const std::vector<DetectorBand>& bands) {
 
 static inline double eq_bell(double f, double fc, double q, double sens_db) {
     if (fc < 1.0 || sens_db == 0) return 1.0;
-    double Qeff = 1.54 * std::pow(q, 1.33);
-    double A = f / fc - fc / f;
-    if (!std::isfinite(A)) return 1.0;
-    double H = 1.0 / std::sqrt(1.0 + (Qeff * A)*(Qeff * A));
-    double gain = std::pow(10.0, sens_db / 40.0);
-    return 1.0 + (gain - 1.0) * H;
+    // RBJ peaking EQ magnitude (FilterGraph) - more accurate than 1/sqrt(1+(Q*A)^2)
+    double w0 = 2.0 * M_PI * fc / 48000.0;
+    double alpha = std::sin(w0) / (2.0 * q);
+    double A = std::pow(10.0, sens_db / 40.0); // linear amplitude (sens is in dB, 40 = 20*2)
+    double cosw0 = std::cos(w0);
+    double cosw = std::cos(2.0 * M_PI * f / 48000.0);
+    // Peaking EQ magnitude squared from RBJ: |H|^2 = (1 + ...)/...
+    // Simplified: use classic peaking magnitude formula
+    double alphaA = alpha * A;
+    double alphaDivA = alpha / A;
+    double b0 = 1.0 + alphaA, b1 = -2.0*cosw0, b2 = 1.0 - alphaA;
+    double a0 = 1.0 + alphaDivA, a1 = -2.0*cosw0, a2 = 1.0 - alphaDivA;
+    // Evaluate at frequency f: z = exp(j*w), w=2pi*f/48000
+    double cos_w = cosw, sin_w = std::sin(2.0 * M_PI * f / 48000.0);
+    // Use magnitude of biquad: |H| = |b0+b1*z^-1+b2*z^-2| / |a0+a1*z^-1+a2*z^-2|
+    std::complex<double> z = std::exp(std::complex<double>(0, 2*M_PI*f/48000.0));
+    std::complex<double> z1 = 1.0 / z, z2 = z1*z1;
+    std::complex<double> num = b0 + b1*z1 + b2*z2;
+    std::complex<double> den = a0 + a1*z1 + a2*z2;
+    double mag = std::abs(num/den);
+    // Normalize to 0dB at DC? RBJ peaking is 0dB at Nyquist, gain at fc
+    // For detector EQ, we want bell that is 1 at far frequencies, gain at fc
+    // So mag is already correct (1 at far, A at fc)
+    return mag;
 }
 
 void FramedDetector::processFrame(const std::complex<double>* spectrum, float* mask) {