151 lines
6.2 KiB
C++
151 lines
6.2 KiB
C++
#include "framed_model.hpp"
|
||
#include "twin.hpp"
|
||
#include "freqpath.hpp"
|
||
#include "rt_mask_tables.hpp"
|
||
#include "rt_weights.hpp"
|
||
#include <cmath>
|
||
#include <cstring>
|
||
#include <algorithm>
|
||
|
||
namespace {
|
||
|
||
// sens XML -> internal sens_stored = sens * 2.054 (NOTES_TWIN:74: XML 12 -> 24.65 dB).
|
||
constexpr float SENS_SCALE = 2.054f;
|
||
|
||
// Live-captured BandConfig parameters from DSP snapshot (2026-08-20).
|
||
// +0x180 (level LUT curve, FUN_180563a60): A = -24.0, B = +28.0, gamma = 1.0, flag = 0.
|
||
// +0x188 (freq-range shaper, FUN_180563440): A = 16.0, B = 20000.0, gamma = 1.0, flag = 0.
|
||
// These values are identical for both render_long.rpp and t1kq_only1_1000 configs.
|
||
// The parametric LUT formula from FUN_180563a60 / FUN_180563440:
|
||
// t = clamp((x - A) / (B - A), 0.0, 1.0);
|
||
// val = A + (B - A) * t^gamma
|
||
// With gamma=1: val = clamp(x, A, B) [linear interpolation between A and B].
|
||
// The x input is the mask-dependent dB-scaled value (mask * 8.6859 from 0x24c43e0).
|
||
|
||
// ctx+0x188 A/B/gamma) evaluated at the measured (xv, C) nodes (al_* dataset +
|
||
// B.11 anchors). Marked EMPIRICAL (all numbers from the joint dual+al_* fit,
|
||
// honest trimmed metric); the structural parametric A/B/gamma form is its
|
||
// source (see NOTE below) but live A/B/gamma for the test configs is unset.
|
||
constexpr double G_FIT = 0.9963;
|
||
constexpr double W_FIT = 0.3335;
|
||
constexpr double A_FIT = 0.9807;
|
||
constexpr double RP0 = 0.0275; // res^rp(Q) gain term, rp = RP0·Q^drp
|
||
constexpr double DRP = 0.2159;
|
||
|
||
// LUT knots (xv = log10(level), level = am/res):
|
||
static constexpr double kLX[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 };
|
||
static constexpr double kLY[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 };
|
||
|
||
static double lut_pchip(double x) {
|
||
int n = 12;
|
||
x = std::min(std::max(x, kLX[0]), kLX[n - 1]);
|
||
// Monotone cubic Hermite (Fritsch–Carlson), matching scipy PchipInterpolator.
|
||
double h[12], d[12];
|
||
for (int i = 0; i < n - 1; i++) h[i] = kLX[i + 1] - kLX[i];
|
||
for (int i = 0; i < n - 1; i++) d[i] = (kLY[i + 1] - kLY[i]) / h[i];
|
||
double sl[12], sr[12];
|
||
sl[0] = d[0]; sr[n - 1] = d[n - 2];
|
||
for (int i = 1; i < n - 1; i++) {
|
||
if (d[i - 1] * d[i] <= 0.0) { sl[i] = sr[i - 1] = 0.0; continue; }
|
||
double w1 = 2 * h[i] + h[i - 1], w2 = h[i] + 2 * h[i - 1];
|
||
sl[i] = (w1 + w2) / (w1 / d[i - 1] + w2 / d[i]);
|
||
sr[i - 1] = sl[i];
|
||
}
|
||
int i = std::upper_bound(kLX, kLX + n, x) - kLX - 1;
|
||
i = std::max(0, std::min(i, n - 2));
|
||
double hh = h[i], t = (x - kLX[i]) / hh;
|
||
double t2 = t * t, t3 = t2 * t;
|
||
double h00 = 2 * t3 - 3 * t2 + 1, h10 = t3 - 2 * t2 + t;
|
||
double h01 = -2 * t3 + 3 * t2, h11 = t3 - t2;
|
||
double y = h00 * kLY[i] + h10 * hh * sr[i] + h01 * kLY[i + 1] + h11 * hh * sl[i + 1];
|
||
return y;
|
||
}
|
||
|
||
// freq-path warp 0x5406a8 (NOTES_LEVEL:181; build_warp): 0.87·K·x/(K+x), K=exp(2.0723).
|
||
static double warp_c(double f) {
|
||
double x = f / 2000.0;
|
||
return 0.87 * 7.942 * x / (7.942 + x);
|
||
}
|
||
|
||
} // namespace
|
||
|
||
FramedDetector::FramedDetector(size_t nfft, float sample_rate)
|
||
: nfft_(nfft), sample_rate_(sample_rate), wsum_(0) {
|
||
am_.resize(nfft / 2 + 1, 0.0f);
|
||
}
|
||
|
||
FramedDetector::~FramedDetector() {}
|
||
|
||
void FramedDetector::setParams(const std::vector<DetectorBand>& bands) {
|
||
bands_ = bands;
|
||
size_t half = nfft_ / 2;
|
||
res_.clear();
|
||
track_.clear();
|
||
|
||
for (const auto& b : bands_) {
|
||
std::vector<float> r(half + 1, 1.0f);
|
||
float sens_lin = std::pow(10.0f, b.sens * SENS_SCALE / 20.0f); // param_5
|
||
detkernel::twin_coeff c = detkernel::build_twin_coeff(
|
||
static_cast<double>(sample_rate_), static_cast<double>(b.fc),
|
||
static_cast<double>(b.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++) {
|
||
r[k] = std::sqrt(out[k].re * out[k].re + out[k].im * out[k].im);
|
||
r[k] = std::max(r[k], 1e-12f);
|
||
}
|
||
res_.push_back(std::move(r));
|
||
}
|
||
track_.assign(bands_.size(), std::vector<float>(half + 1, 1.0f));
|
||
}
|
||
|
||
void FramedDetector::processFrame(const std::complex<double>* spectrum, float* mask) {
|
||
size_t half = nfft_ / 2;
|
||
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_));
|
||
|
||
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);
|
||
}
|
||
|
||
for (size_t k = 0; k <= half; k++) mask[k] = 1.0f;
|
||
|
||
for (size_t b = 0; b < bands_.size(); b++) {
|
||
double rp = RP0 * std::pow(static_cast<double>(bands_[b].q), DRP);
|
||
double fk = 0.0;
|
||
double fstep = (sample_rate_ * 0.5) / static_cast<double>(half);
|
||
for (size_t k = 0; k <= half; k++) {
|
||
double res_k = std::max(static_cast<double>(res_[b][k]), 1e-12);
|
||
double lvl = static_cast<double>(am_[k]) / res_k;
|
||
double xv = std::log10(std::max(lvl, 1e-9));
|
||
double C = G_FIT * lut_parametric(xv, CAP_A_LEVEL, CAP_B_LEVEL, CAP_GAMMA) + W_FIT * std::pow(warp_c(fk), A_FIT);
|
||
double g = std::max(1.0 - C, 1e-9) * std::pow(res_k, rp);
|
||
mask[k] = std::min(static_cast<float>(g), mask[k]);
|
||
fk += fstep;
|
||
}
|
||
}
|
||
for (size_t k = half + 1; k < nfft_; k++) {
|
||
mask[k] = mask[nfft_ - k];
|
||
}
|
||
} |