P4: port framed_render.py model to C++ (FramedDetector: twin-res -> level -> Pchip LUT -> warp -> res^rp gain), replace empirical Detector; tt_base 49.93%->61.45% (single-band); multi-band refs still 100% (model is single-band)

2026-08-20 10:31:43 +03:00
parent 45dfe2b8c2
commit 0cdc57972c
6 changed files with 190 additions and 9 deletions
@@ -0,0 +1,136 @@
+#include "framed_model.hpp"
+#include "twin.hpp"
+#include "freqpath.hpp"
+#include <cmath>
+#include <cstring>
+#include <algorithm>
+
+namespace {
+
+constexpr float G_FIT = 0.9963f;
+constexpr float W_FIT = 0.3335f;
+constexpr float A_FIT = 0.9807f;
+constexpr float GAIN = 4.132f;              // sqrtf(param_5), sens_lin = GAIN^2
+
+// framed_render.py LUT nodes (Pchip): C = G_FIT*LUT(xv) + W_FIT*warp^A_FIT
+constexpr double LX[12] = {-0.75, -0.5012, -0.5, -0.2012, 0.0988, 0.2488,
+                           0.3988, 0.5488, 0.574, 0.61, 0.75, 1.0};
+constexpr double LY[12] = {0.4402, 0.366, 0.4552, 0.459, 0.541, 0.576,
+                           0.608, 0.636, 0.5645, 0.6471, 0.6562, 0.6670};
+constexpr int LN = 12;
+
+// PchipInterpolator (monotone cubic Hermite, scipy semantics) eval at one x.
+double pchip_eval(double x) {
+    if (x <= LX[0]) return LY[0];
+    if (x >= LX[LN - 1]) return LY[LN - 1];
+    // slopes
+    double m[LN];
+    for (int i = 0; i < LN - 1; i++) m[i] = (LY[i + 1] - LY[i]) / (LX[i + 1] - LX[i]);
+    double d[LN];
+    d[0] = m[0];
+    for (int i = 1; i < LN - 1; i++) {
+        if (m[i - 1] * m[i] <= 0.0) d[i] = 0.0;
+        else {
+            double w1 = 2 * LX[i + 1] - LX[i] - LX[i - 1];
+            double w2 = LX[i + 1] - LX[i - 1];
+            d[i] = (w1 + w2) / ((w1 / m[i - 1]) + (w2 / m[i]));
+        }
+    }
+    d[LN - 1] = m[LN - 2];
+    int i = 0;
+    while (i < LN - 2 && x > LX[i + 1]) i++;
+    double h = LX[i + 1] - LX[i];
+    double t = (x - LX[i]) / h;
+    double y0 = LY[i], y1 = LY[i + 1], d0 = d[i], d1 = d[i + 1];
+    double h00 = 2 * t * t * t - 3 * t * t + 1;
+    double h10 = t * t * t - 2 * t * t + t;
+    double h01 = -2 * t * t * t + 3 * t * t;
+    double h11 = t * t * t - t * t;
+    double v = h00 * y0 + h10 * h * d0 + h01 * y1 + h11 * h * d1;
+    return std::max(LY[LN - 2] < LY[LN - 1] ? 0.0 : -1e9,
+                    std::min(v, *std::max_element(LY, LY + LN) * 1.0));
+}
+
+double lut_eval(double x) {
+    double v = pchip_eval(x);
+    double lo = *std::min_element(LY, LY + LN);
+    double hi = *std::max_element(LY, LY + LN);
+    return std::max(lo, std::min(hi, v));
+}
+
+}  // namespace
+
+FramedDetector::FramedDetector(size_t nfft, float sample_rate)
+    : nfft_(nfft), sample_rate_(sample_rate), fc_(0), q_(1.0), wsum_(0), inited_(false) {
+    res_.resize(nfft / 2 + 1, 1.0f);
+    warp_.resize(nfft, 0.0f);
+    am_.resize(nfft / 2 + 1, 0.0f);
+}
+
+FramedDetector::~FramedDetector() {}
+
+void FramedDetector::setParams(float fc, float q) {
+    fc_ = fc;
+    q_ = std::max(1.0f, q);
+    size_t half = nfft_ / 2;
+
+    // twin resonance per bin (res = |2B(z)/A(z)|)
+    float sens_lin = GAIN * GAIN;  // param_5
+    detkernel::twin_coeff c = detkernel::build_twin_coeff(
+        static_cast<double>(sample_rate_), static_cast<double>(fc),
+        static_cast<double>(q_), sens_lin);
+    std::vector<detkernel::cplxf> z(half + 1);
+    std::vector<detkernel::cplxf> out(half + 1);
+    for (size_t k = 0; k <= half; k++) {
+        double theta = 2.0 * M_PI * static_cast<double>(k) / static_cast<double>(nfft_);
+        z[k].re = static_cast<float>(std::cos(theta));
+        z[k].im = static_cast<float>(std::sin(theta));
+    }
+    detkernel::twin_apply(c, z.data(), half + 1, out.data());
+    for (size_t k = 0; k <= half; k++) {
+        res_[k] = std::sqrt(out[k].re * out[k].re + out[k].im * out[k].im);
+        res_[k] = std::max(res_[k], 1e-12f);
+    }
+
+    // warp (freqpath 0x5406a8)
+    detkernel::build_warp(static_cast<float>(sample_rate_),
+                          static_cast<int>(nfft_), warp_.data());
+
+    inited_ = true;
+}
+
+void FramedDetector::processFrame(const std::complex<double>* spectrum, float* mask) {
+    size_t half = nfft_ / 2;
+    // wsum = sum of sqrt-hann window (compute once via nfft)
+    if (wsum_ == 0.0) {
+        double s = 0.0;
+        for (size_t i = 0; i < nfft_; i++) {
+            s += std::sqrt(0.5 * (1.0 - std::cos(2.0 * M_PI * i / (nfft_ - 1))));
+        }
+        wsum_ = s;
+    }
+
+    double tatt = 0.011, trel = 0.08;
+    double att = std::exp(-1.0 * (nfft_ / 4) / (tatt * sample_rate_));
+    double rel = std::exp(-1.0 * (nfft_ / 4) / (trel * sample_rate_));
+
+    double rp = 0.0275 * std::pow(static_cast<double>(q_), 0.2159);
+
+    for (size_t k = 0; k <= half; k++) {
+        double a_cur = 2.0 * std::abs(spectrum[k]) / wsum_;
+        double am = am_[k];
+        if (a_cur > am) am = att * am + (1.0 - att) * a_cur;
+        else am = rel * am + (1.0 - rel) * a_cur;
+        am_[k] = static_cast<float>(am);
+
+        double xv = std::log10(std::max(am / static_cast<double>(res_[k]), 1e-9));
+        double C = G_FIT * lut_eval(xv) +
+                   W_FIT * std::pow(static_cast<double>(warp_[k]), A_FIT);
+        double gain = std::max(1.0 - C, 1e-9) *
+                      std::pow(static_cast<double>(res_[k]), rp);
+        mask[k] = static_cast<float>(gain);
+    }
+    for (size_t k = half + 1; k < nfft_; k++) {
+        mask[k] = mask[nfft_ - k];
+    }
+}