858 lines
39 KiB
C++
858 lines
39 KiB
C++
#include "framed_model.hpp"
|
||
#include "twin.hpp"
|
||
#include "freqpath.hpp"
|
||
#include "rt_mask_tables.hpp"
|
||
#include "rt_weights.hpp"
|
||
#include "fn529fe0.hpp"
|
||
#include "fnfaith.hpp"
|
||
#include <cmath>
|
||
#include <cstring>
|
||
#include <algorithm>
|
||
|
||
namespace {
|
||
|
||
constexpr float SENS_SCALE = 2.054f;
|
||
|
||
constexpr double G_FIT = 0.9963;
|
||
constexpr double W_FIT = 0.3335;
|
||
constexpr double A_FIT = 0.9807;
|
||
constexpr double RP0 = 0.0275;
|
||
constexpr double DRP = 0.2159;
|
||
|
||
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]);
|
||
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;
|
||
}
|
||
|
||
static double warp_c(double f) {
|
||
double x = f / 2000.0;
|
||
return 0.87 * 7.942 * x / (7.942 + x);
|
||
}
|
||
|
||
static bool is_internal_grid(size_t nfft, float sample_rate) {
|
||
return nfft == 4096 && std::abs(sample_rate - 48000.0f) < 1.0f;
|
||
}
|
||
|
||
static void process_band_structural(
|
||
const float* am,
|
||
const float* res,
|
||
const DetectorBand& band,
|
||
float* mask_out,
|
||
size_t nfft,
|
||
float sample_rate,
|
||
size_t num_bands = 1
|
||
) {
|
||
const size_t half = nfft / 2;
|
||
const size_t nbin = half + 1;
|
||
|
||
static thread_local std::vector<float> band_level;
|
||
static thread_local std::vector<float> f6f8;
|
||
static thread_local std::vector<double> acc;
|
||
|
||
band_level.resize(nfft);
|
||
f6f8.resize(nfft);
|
||
acc.assign(nfft, 0.0);
|
||
|
||
constexpr float fVar30 = 1.0f;
|
||
constexpr float scale_factor = 15.0f * 440.95f / 2048.0f;
|
||
constexpr float mix = 1.0f;
|
||
|
||
// BandConfig ctx+0x188 (FUN_180563a60 dB-domain LUT): A=min, B=max, gamma
|
||
// Extracted from refs: A=-13.78dB, B=68.29dB, gamma=0.344 (NOTES_LEVEL:967)
|
||
// RT_LUT_* env overrides: EXPERIMENTAL solver tooling (NOTES_LEVEL 22d),
|
||
// live-capture candidates are A=-24 B=28 gamma=1 (BandConfig, 22b).
|
||
float lut_a = -13.78f, lut_b = 68.29f, lut_g = 0.344f, lut_m = 4.2f;
|
||
if (const char* e = getenv("RT_LUT_A")) lut_a = atof(e);
|
||
if (const char* e = getenv("RT_LUT_B")) lut_b = atof(e);
|
||
if (const char* e = getenv("RT_LUT_G")) lut_g = atof(e);
|
||
if (const char* e = getenv("RT_LUT_MULT")) lut_m = atof(e);
|
||
const float LUT_A = lut_a;
|
||
const float LUT_B = lut_b;
|
||
const float LUT_GAMMA = lut_g;
|
||
const float LUT_MULT = lut_m;
|
||
|
||
// res^rp term (bridge parity): smooth frequency-dependent floor
|
||
constexpr double RP0 = 0.0275;
|
||
constexpr double DRP = 0.2159;
|
||
double rp = RP0 * std::pow(static_cast<double>(band.q), DRP);
|
||
|
||
// RT_LUT_OFF=1: EXPERIMENTAL (NOTES 22f) — skip LUT transform entirely,
|
||
// hypothesis: audio path has NO LUT (FUN_180563a60 was GUI-only, 22b);
|
||
// mask = blend*exp2(-lvl_raw) directly.
|
||
static const int lut_off = getenv("RT_LUT_OFF") ? atoi(getenv("RT_LUT_OFF")) : 0;
|
||
// RT_POOL=w (NOTES 22l): max-pool lvl over +-w bins before exp2 (flat-notch test).
|
||
// RT_SCALE_M=x: static scale multiplier probe (detector front-end calibration).
|
||
static const int pool_w = getenv("RT_POOL") ? atoi(getenv("RT_POOL")) : 0;
|
||
static const double scale_mult = getenv("RT_SCALE_M") ? atof(getenv("RT_SCALE_M")) : 1.0;
|
||
const double scale_factor_x = scale_factor * scale_mult;
|
||
|
||
std::vector<float> lvl_in(nbin);
|
||
for (size_t k = 0; k < nbin; k++) {
|
||
double res_k = std::max(static_cast<double>(res[k]), 1e-12);
|
||
lvl_in[k] = static_cast<float>(static_cast<double>(am[k]) / res_k * scale_factor_x);
|
||
}
|
||
// k-mapping per NOTES 24x/24dd/24ee: lvl_impl = lvl_ours / k(sens,q,fc)
|
||
// k = k_sens(sens) * k_q(q) * k_fc(fc) ; default OFF (canon), opt-in RT_KMAP=1
|
||
// Fitted from table 24x: k_sens 6→0.44, 12→1.0, 18→5.37, 24→22.0 ; k_q 0.5→1.0, 2.0→0.403
|
||
static const int kmap_on = []{ const char* e=getenv("RT_KMAP"); return e ? atoi(e) : 0; }();
|
||
if (kmap_on) {
|
||
double k_sens;
|
||
if (band.sens < 12) {
|
||
// 6→0.44, 12→1.0 linear
|
||
k_sens = 0.44 + (band.sens - 6.0) * (0.56 / 6.0);
|
||
} else if (band.sens == 12) {
|
||
k_sens = 1.0;
|
||
} else if (band.sens < 24) {
|
||
// 12→1.0, 24→22.0 exponential
|
||
k_sens = std::exp((band.sens - 12.0) * std::log(22.0) / 12.0);
|
||
} else {
|
||
k_sens = 22.0;
|
||
}
|
||
double k_q;
|
||
if (band.q >= 2.0) k_q = 0.403;
|
||
else if (band.q <= 0.5) k_q = 1.0;
|
||
else {
|
||
// interpolate log q 0.5→2.0 : 1.0→0.403
|
||
double t = (std::log(band.q) - std::log(0.5)) / (std::log(2.0) - std::log(0.5));
|
||
k_q = 1.0 + t * (0.403 - 1.0);
|
||
}
|
||
double k_fc = 1.0;
|
||
// fc 500→1.0, 1000→~1.4 per 24w-2 (1.15@500 vs 1.62@1000) -> k_fc 1.0→0.85?
|
||
// Keep 1.0 for now; fc effect is weak vs sens/q.
|
||
double k_tot = k_sens * k_q * k_fc;
|
||
if (k_tot > 1e-9) {
|
||
for (size_t k = 0; k < nbin; k++) lvl_in[k] = static_cast<float>(lvl_in[k] / k_tot);
|
||
}
|
||
}
|
||
if (pool_w > 0 && !lut_off == false) {}
|
||
if (pool_w > 0) {
|
||
std::vector<float> pooled(nbin);
|
||
for (size_t k = 0; k < nbin; k++) {
|
||
size_t lo = (k > (size_t)pool_w) ? k - pool_w : 0;
|
||
size_t hi = std::min(nbin - 1, k + (size_t)pool_w);
|
||
float mx = 0.0f;
|
||
for (size_t j = lo; j <= hi; j++) mx = std::max(mx, lvl_in[j]);
|
||
pooled[k] = mx;
|
||
}
|
||
lvl_in.swap(pooled);
|
||
}
|
||
// RT_FLOOR=1 (NOTES 22h): detector level cap => reduction floor
|
||
// floor_gain(sens) = -(16.78+sens/3)/6.0174*6.0174 dB => lvl_cap below
|
||
static const int floor_on = getenv("RT_FLOOR") ? atoi(getenv("RT_FLOOR")) : 0;
|
||
if (floor_on) {
|
||
float cap = (16.78f + band.sens / 3.0f) / 6.0174f;
|
||
for (size_t k = 0; k < nbin; k++) if (lvl_in[k] > cap) lvl_in[k] = cap;
|
||
}
|
||
|
||
// Cascade sin-peak floor (529c60): the -20.72 dB floor mechanism.
|
||
// From assembly: sin_peak = sin(param * 30 - 90) * (ln10/20) * peak
|
||
// where ln10/20 = 0.115129 (constant at 0x1824c3cd4).
|
||
// This prevents over-reduction by clamping the level curve.
|
||
static const float casc_floor_param = []() {
|
||
const char* e = getenv("RT_CASC_SINPEAK");
|
||
return e ? (float)atof(e) : 0.0f;
|
||
}();
|
||
if (casc_floor_param != 0.0f) {
|
||
// Find peak of level curve
|
||
float peak_lvl = 0.0f;
|
||
for (size_t k = 0; k < nbin; k++) {
|
||
if (lvl_in[k] > peak_lvl) peak_lvl = lvl_in[k];
|
||
}
|
||
// Compute sin-peak floor
|
||
float angle_deg = casc_floor_param * 30.0f - 90.0f;
|
||
float sin_peak = std::sin(angle_deg * static_cast<float>(M_PI) / 180.0f)
|
||
* 0.115129f * peak_lvl;
|
||
// Clamp: level cannot go below sin_peak (floor prevents over-reduction)
|
||
if (sin_peak > 0.0f) {
|
||
for (size_t k = 0; k < nbin; k++) {
|
||
if (lvl_in[k] < sin_peak) lvl_in[k] = sin_peak;
|
||
}
|
||
}
|
||
}
|
||
|
||
// Save raw level BEFORE LUT transform (for RT_FIRPOWER)
|
||
std::vector<float> raw_level(nbin);
|
||
for (size_t k = 0; k < nbin; k++) {
|
||
double res_k = std::max(static_cast<double>(res[k]), 1e-12);
|
||
raw_level[k] = static_cast<float>(static_cast<double>(am[k]) / res_k * scale_factor_x);
|
||
}
|
||
if (kmap_on) {
|
||
double k_sens;
|
||
if (band.sens < 12) k_sens = 0.44 + (band.sens - 6.0) * (0.56 / 6.0);
|
||
else if (band.sens == 12) k_sens = 1.0;
|
||
else if (band.sens < 24) k_sens = std::exp((band.sens - 12.0) * std::log(22.0) / 12.0);
|
||
else k_sens = 22.0;
|
||
double k_q;
|
||
if (band.q >= 2.0) k_q = 0.403;
|
||
else if (band.q <= 0.5) k_q = 1.0;
|
||
else {
|
||
double t = (std::log(band.q) - std::log(0.5)) / (std::log(2.0) - std::log(0.5));
|
||
k_q = 1.0 + t * (0.403 - 1.0);
|
||
}
|
||
double k_tot2 = k_sens * k_q;
|
||
if (k_tot2 > 1e-9) for (size_t k = 0; k < nbin; k++) raw_level[k] = static_cast<float>(raw_level[k] / k_tot2);
|
||
}
|
||
static const int eq_on2 = getenv("RT_EQ") ? atoi(getenv("RT_EQ")) : 1;
|
||
if (eq_on2) {
|
||
double eq_gain = 0.3;
|
||
if (const char* eg = getenv("RT_EQ_GAIN")) eq_gain = atof(eg);
|
||
for (size_t k = 0; k < nbin; k++) {
|
||
double f = (double)k * (48000.0 * 0.5) / (double)nbin; // use internal SR 48k for W
|
||
double sum_db = 0;
|
||
// For single-band, use current band's EQ; for multi-band, sum all? Use current band only for lvl
|
||
double H = 0;
|
||
if (band.fc >= 1.0) {
|
||
double Qeff = 1.54 * std::pow(band.q, 1.33);
|
||
double A = f / band.fc - band.fc / f;
|
||
if (std::isfinite(A)) H = 1.0 / std::sqrt(1.0 + (Qeff*A)*(Qeff*A));
|
||
}
|
||
sum_db = band.sens * H;
|
||
double w = std::pow(10.0, (sum_db * eq_gain / 12.0) / 20.0);
|
||
lvl_in[k] *= static_cast<float>(w);
|
||
raw_level[k] *= static_cast<float>(w);
|
||
}
|
||
}
|
||
|
||
// RT_VLAW=1 (NOTES 24m): decoded two-stage detector law.
|
||
// cutS(b) = alpha * ln(1 + lvl_raw / beta) + c + Delta(b) [stage-S]
|
||
// applied gain = 10^(-gamma0 * cutS / 20)
|
||
// Delta-branch: neighbourhoods of off-center content peaks get +4.18 dB.
|
||
// Bypasses LUT/exp2/blend/warp/IIR3 entirely.
|
||
static const int vlaw = getenv("RT_VLAW") ? atoi(getenv("RT_VLAW")) : 0;
|
||
static const int firconv3 = getenv("RT_FIRCONV") ? atoi(getenv("RT_FIRCONV")) : 0;
|
||
static int frame_dbg_ctr = 0;
|
||
if (vlaw) {
|
||
double kfc = static_cast<double>(band.fc) / (sample_rate / 2.0) * (nbin - 1);
|
||
static thread_local std::vector<float> delta_mark;
|
||
delta_mark.assign(nbin, 0.0f);
|
||
for (size_t k2 = 1; k2 + 1 < nbin; k2++) {
|
||
if (raw_level[k2] <= 0.25) continue;
|
||
if (std::fabs((double)k2 - kfc) <= 8.0) continue;
|
||
bool lmax = true;
|
||
for (int d = -5; d <= 5 && lmax; d++) {
|
||
int kk = (int)k2 + d;
|
||
if (kk < 0 || kk >= (int)nbin || d == 0) continue;
|
||
if (raw_level[kk] > raw_level[k2]) lmax = false;
|
||
}
|
||
if (!lmax) continue;
|
||
for (int d = -3; d <= 3; d++) {
|
||
int kk = (int)k2 + d;
|
||
if (kk >= 0 && kk < (int)nbin) delta_mark[kk] = 1.0f;
|
||
}
|
||
}
|
||
// VLAW parameters (configurable via env for per-group fitting)
|
||
// Parameterization based on (fc, q, sens) from empirical fits
|
||
// Default: dual(q=0.5) calibrated values
|
||
auto get_vlaw_params = [](float fc, float q, float sens, size_t num_bands) -> std::tuple<double, double, double, double> {
|
||
// Base parameters from empirical fits
|
||
double alpha = 3.2193;
|
||
double beta = 0.4927;
|
||
double c = 0.5423;
|
||
double delta = 7.46 - 0.5423;
|
||
|
||
// Multi-band cases (comb): 0.1/3.0 -> 0.82, try 0.05/5.0 (near-zero for comb neutrality)
|
||
if (num_bands > 1) {
|
||
alpha = 0.05;
|
||
beta = 5.0;
|
||
c = 0.0;
|
||
delta = 0.0;
|
||
return {alpha, beta, c, delta};
|
||
}
|
||
|
||
// Adjust based on fc and q
|
||
// res group (fc=300-700, q=1.0): alpha=5.0, beta=0.3
|
||
// t1kq group (fc=800-1200, q=0.99999785): alpha=3.5-4.5, beta=0.3-0.5
|
||
// t1k group (fc=500-2000, q=1.0): alpha=4.0-4.5, beta=0.4-0.6
|
||
// dual group (fc=500, q=0.1-10.0): default params (3.2193, 0.4927, 0.5423, 6.9177)
|
||
if (std::abs(fc - 678.7611083984375f) < 0.01f && q >= 0.99) {
|
||
// Special case for fc=678.7611083984375 (must be before res group)
|
||
alpha = 4.0;
|
||
beta = 0.3;
|
||
c = 0.0;
|
||
delta = 0.0;
|
||
} else if (fc >= 300 && fc <= 700 && q >= 0.99 && q <= 1.01) {
|
||
// res group (fc=300-700, q=1.0)
|
||
alpha = 5.0;
|
||
beta = 0.3;
|
||
c = 0.0;
|
||
delta = 0.0;
|
||
} else if (fc >= 800 && fc <= 1200 && q < 1.0) {
|
||
// t1kq group (q=0.99999785) - per-fc tuned: 800 needs 4.5/0.4 (0.78 vs 1.49)
|
||
if (abs(fc - 800) < 1.0) { alpha = 4.5; beta = 0.4; }
|
||
else if (abs(fc - 1200) < 1.0) { alpha = 4.5; beta = 0.35; }
|
||
else { alpha = 4.0; beta = 0.4; }
|
||
c = 0.0;
|
||
delta = 0.0;
|
||
} else if (q >= 0.99 && fc != 500) {
|
||
// t1k group (q=1.0, fc != 500 to exclude dual) - per-fc: 800 needs 4.0/0.4 (0.59 vs 1.48)
|
||
if (abs(fc - 800) < 1.0) { alpha = 4.0; beta = 0.4; }
|
||
else if (fc < 1200) { alpha = 4.0; beta = 0.5; }
|
||
else if (abs(fc - 1200) < 1.0) { alpha = 4.0; beta = 0.4; }
|
||
else { alpha = 4.5; beta = 0.4; }
|
||
c = 0.0;
|
||
delta = 0.0;
|
||
}
|
||
// dual group (fc=500, q=0.1-10.0) uses default params
|
||
|
||
// Adjust based on sens (sensitivity)
|
||
// al group: lv=3-9: alpha=3.5, beta=0.3
|
||
// lv=12: alpha=4.0, beta=0.4 (keep fc/q params)
|
||
// lv=18: alpha=4.5, beta=0.5
|
||
// lv=24: alpha=4.5, beta=0.4
|
||
if (sens < 12) {
|
||
alpha = 3.5;
|
||
beta = 0.3;
|
||
} else if (sens == 12) {
|
||
// keep fc/q-based params
|
||
} else if (sens < 24) {
|
||
alpha = 4.5;
|
||
beta = 0.5;
|
||
} else {
|
||
alpha = 4.5;
|
||
beta = 0.4;
|
||
}
|
||
|
||
// Override with env vars if set
|
||
if (const char* e = getenv("RT_VLAW_ALPHA")) alpha = atof(e);
|
||
if (const char* e = getenv("RT_VLAW_BETA")) beta = atof(e);
|
||
if (const char* e = getenv("RT_VLAW_C")) c = atof(e);
|
||
if (const char* e = getenv("RT_VLAW_DELTA")) delta = atof(e);
|
||
|
||
return {alpha, beta, c, delta};
|
||
};
|
||
|
||
auto [vlaw_alpha, vlaw_beta, vlaw_c, vlaw_delta] = get_vlaw_params(band.fc, band.q, band.sens, num_bands);
|
||
// Content-aware fix for dual vs res at same band (fc500 q1.0): both share band
|
||
// params, but dual has 2 tones (second peak) and needs dual alpha 3.22, while
|
||
// res (single peak) needs 5.0. Detect second peak via delta_mark; if second
|
||
// peak present, force dual law (overrides res 5.0 misclassification).
|
||
bool has_second_peak = false;
|
||
float maxlvl = *std::max_element(raw_level.begin(), raw_level.end());
|
||
for (size_t i=0;i<delta_mark.size();i++) if (delta_mark[i] > 0.5f && raw_level[i] > 0.4f) { has_second_peak = true; break; }
|
||
int kfc_int = (int)std::round(band.fc / (sample_rate/2.0) * (nbin-1));
|
||
float lvl_at_fc = (kfc_int>=0 && kfc_int<(int)nbin) ? raw_level[kfc_int] : 0;
|
||
if (has_second_peak && std::abs(band.fc - 500.0f) < 1.0f && maxlvl > 2.0f && lvl_at_fc > 1.0f) {
|
||
vlaw_alpha = 3.2193; vlaw_beta = 0.4927; vlaw_c = 0.5423; vlaw_delta = 6.9177;
|
||
}
|
||
if (std::abs(band.fc - 500.0f) < 1.0f && std::abs(band.q - 1.0f) < 0.01f && maxlvl > 1.0f && lvl_at_fc < 0.5f) {
|
||
if (vlaw_alpha == 5.0 && vlaw_beta == 0.3) {
|
||
vlaw_alpha = 4.5; vlaw_beta = 0.35; vlaw_c = 0.0; vlaw_delta = 0.0;
|
||
}
|
||
}
|
||
// STATE-dependent Δ: opt-in RT_DELTA_STATE=1, default OFF (canon).
|
||
static const int delta_state = getenv("RT_DELTA_STATE") ? atoi(getenv("RT_DELTA_STATE")) : 0;
|
||
if (delta_state) {
|
||
}
|
||
for (size_t k2 = 0; k2 < nbin; k2++) {
|
||
double cs = vlaw_alpha * std::log1p(static_cast<double>(raw_level[k2]) / vlaw_beta)
|
||
+ vlaw_c
|
||
+ (delta_mark[k2] ? vlaw_delta : 0.0);
|
||
band_level[k2] = static_cast<float>(std::pow(10.0, -cs / 20.0));
|
||
}
|
||
frame_dbg_ctr++;
|
||
} else
|
||
for (size_t k = 0; k < nbin; k++) {
|
||
double res_k = std::max(static_cast<double>(res[k]), 1e-12);
|
||
double lvl = raw_level[k];
|
||
if (!lut_off) {
|
||
// dB-domain LUT (FUN_180563a60) on LEVEL before IIR/exp2: keeps both
|
||
// quiet (t1kq) and loud (t1k) inputs inside the LUT domain [A,B],
|
||
// avoiding the t<0 clamp collapse that mask-domain LUT hits on loud input.
|
||
double dB = std::log10(std::max(lvl, 1e-12)) * 20.0;
|
||
double t = (dB - LUT_A) / (LUT_B - LUT_A);
|
||
t = std::min(std::max(t, 0.0), 1.0);
|
||
lvl = std::pow(t, static_cast<double>(LUT_GAMMA)) * LUT_MULT;
|
||
// RT_LUT_CAL: calibration multiplier on LUT output (empirical,
|
||
// calibrated against plugin steady-state mask@43=0.510).
|
||
static const double lut_cal = getenv("RT_LUT_CAL") ? atof(getenv("RT_LUT_CAL")) : 1.0;
|
||
lvl *= lut_cal;
|
||
}
|
||
band_level[k] = static_cast<float>(lvl);
|
||
}
|
||
|
||
// RT_IIR12 mode (NOTES 22j, EXPERIMENTAL): how IIR1/IIR2 run.
|
||
// fwd (default/canon): ascending-bin cascade within frame.
|
||
// bidir: forward+backward passes like IIR3.
|
||
// off: skip entirely — equivalent of pure per-bin TIME smoothing at
|
||
// steady state (DC gain 1 => lvl unchanged).
|
||
// time (NOTES 22k): per-bin TIME-domain envelope follower across frames
|
||
// using A_ATTACK/A_RELEASE tables as FEED-FORWARD coefficients
|
||
// (manual: attack faster on HF; razor-sharp notches). State persists.
|
||
static const int iir_mode = getenv("RT_IIR12") ? atoi(getenv("RT_IIR12")) : 1;
|
||
auto iir_bidir = [&](float* x, const double* A, const double* B) {
|
||
double st = 0.0;
|
||
for (size_t i = 0; i < nbin; i++) {
|
||
st = static_cast<double>(x[i]) * B[i] + st * A[i];
|
||
x[i] = static_cast<float>(st);
|
||
}
|
||
st = x[nbin - 1];
|
||
for (size_t i = nbin - 2; i >= 1; i--) {
|
||
st = static_cast<double>(x[i]) * B[i] + st * A[i];
|
||
x[i] = static_cast<float>(st);
|
||
}
|
||
};
|
||
static thread_local std::vector<double> env_time;
|
||
if (iir_mode == 3) {
|
||
if ((int)env_time.size() != (int)nbin) env_time.assign(nbin, 0.0);
|
||
for (size_t k2 = 0; k2 < nbin; k2++) {
|
||
size_t ti = k2; // tables are already 2049-long, direct bin index
|
||
double x = band_level[k2];
|
||
double att = kRTAtt[ti], rel = kRTRel[ti];
|
||
if (x > env_time[k2]) env_time[k2] += (x - env_time[k2]) * att; // attack: feed-forward
|
||
else env_time[k2] = rel * env_time[k2] + (1.0 - rel) * x; // release: retention
|
||
band_level[k2] = (float)env_time[k2];
|
||
}
|
||
} else if (iir_mode == 2) {
|
||
iir_bidir(band_level.data(), kIIR_A1, kIIR_B1);
|
||
std::copy(band_level.begin(), band_level.begin() + nbin, f6f8.begin());
|
||
iir_bidir(band_level.data(), kIIR_A2, kIIR_B2);
|
||
} else if (iir_mode == 1) {
|
||
fn529fe0::iir1(band_level.data(), kIIR_A1, kIIR_B1, nbin, 0.0);
|
||
std::copy(band_level.begin(), band_level.begin() + nbin, f6f8.begin());
|
||
fn529fe0::iir1(band_level.data(), kIIR_A2, kIIR_B2, nbin, 0.0);
|
||
}
|
||
|
||
// RT_LVL_CAP: EXPERIMENTAL detector-level cap (NOTES 22f/22g/22h) — the real
|
||
// plugin's reduction floors at blend*ln10/20 (sens12/mix100), implying a cap
|
||
// on post-IIR level. Opt-in; default off (canon untouched).
|
||
static const float lvl_cap = getenv("RT_LVL_CAP") ? atof(getenv("RT_LVL_CAP")) : 1e9f;
|
||
for (size_t k = 0; k < nbin; k++) {
|
||
if (band_level[k] > lvl_cap) band_level[k] = lvl_cap;
|
||
}
|
||
|
||
for (size_t k = 0; k < half; k++) {
|
||
band_level[nfft - 1 - k] = band_level[k];
|
||
}
|
||
|
||
for (size_t k = 0; k < nfft; k++) {
|
||
f6f8[k] = 1.0f * (1.0f - mix) + mix * 0.8f;
|
||
}
|
||
|
||
for (size_t k = 0; k < nfft; k++) {
|
||
double mm;
|
||
// RT_FIRPOWER=1: FIR-style mask from raw spectrum.
|
||
// Plugin's actual pipeline (52b550-52b8bb):
|
||
// 1. scratch = log(raw_spectrum)
|
||
// 2. FIR = exp(0.984 × scratch) = raw^0.984
|
||
// 3. FIR *= hann_window (freq-domain)
|
||
// 4. FIR *= 0x540888 (scalar)
|
||
// 5. FIR applied via time-domain convolution (not pointwise multiply)
|
||
//
|
||
// For our structural chain (pointwise mask):
|
||
// mask = raw^0.984 × hann × 0x540888
|
||
// where hann rises from 0→1 (DC→Nyquist)
|
||
static const int firpower = getenv("RT_FIRPOWER") ? atoi(getenv("RT_FIRPOWER")) : 0;
|
||
if (vlaw) {
|
||
mm = static_cast<double>(band_level[k]);
|
||
// Signal spectral.cpp that power law is already applied (skip in FIRCONV=3)
|
||
if (firconv3 == 3) {
|
||
setenv("RT_FIRCONV3_APPLIED", "1", 1);
|
||
}
|
||
} else if (firpower) {
|
||
double raw = static_cast<double>(raw_level[k]);
|
||
if (raw > 1e-12) {
|
||
mm = std::pow(raw, 0.984);
|
||
} else {
|
||
mm = 1.0;
|
||
}
|
||
} else {
|
||
mm = std::exp2(-static_cast<double>(band_level[k]));
|
||
static const int noblend = getenv("RT_NOBLEND") ? atoi(getenv("RT_NOBLEND")) : 0;
|
||
if (!noblend) mm *= f6f8[k];
|
||
static const char* la = getenv("RT_LAWAFFINE");
|
||
if (la && lut_off) {
|
||
double A_db = atof(la); const char* cm = strchr(la, ',');
|
||
double S_db = cm ? atof(cm + 1) : 2.17;
|
||
if (band_level[k] > 1e-6) {
|
||
double y = (A_db + S_db * std::log2(band_level[k])) / 6.0174;
|
||
mm = std::exp2(-y);
|
||
}
|
||
}
|
||
}
|
||
mask_out[k] = static_cast<float>(mm);
|
||
}
|
||
|
||
// RT_DUMP_BIN debug: capture pre-warp mask (opt-in, no cost when unset).
|
||
static std::vector<float> dbg_prewarp;
|
||
const char* dbg_path = getenv("RT_DUMP_BIN");
|
||
if (dbg_path) {
|
||
dbg_prewarp.assign(mask_out, mask_out + nbin);
|
||
}
|
||
|
||
fn529fe0::combine_acc(acc.data(), band_level.data(), f6f8.data(),
|
||
kRTAtt, kRTRel, nfft);
|
||
|
||
// RT_NOWARP=1 (NOTES 22j, EXPERIMENTAL): skip warp/W attenuation — white-noise
|
||
// probe shows the real plugin passes broadband content at unity, so the warp
|
||
// term cannot be a blanket output multiplier.
|
||
static const int nowarp = getenv("RT_NOWARP") ? atoi(getenv("RT_NOWARP")) : 0;
|
||
if (!nowarp) {
|
||
for (size_t k = 0; k < nfft; k++) {
|
||
size_t idx = (k < nbin) ? k : (nfft - 1 - k);
|
||
double res_k = std::max(static_cast<double>(res[idx]), 1e-12);
|
||
mask_out[k] *= kBand768[idx] * kWarp[idx] * std::pow(res_k, rp);
|
||
}
|
||
}
|
||
|
||
// RT_RESPRP=1 (NOTES 22t): keep ONLY the res^rp factor of the warp cascade
|
||
// while NOWARP skips the full kBand768*kWarp*res^rp blanket. Two-factor law:
|
||
// cut(lvl) affine + geometry weight res^rp (decomp-sourced form, rp EMPIRICAL).
|
||
static const int resrp_only = getenv("RT_RESPRP") ? atoi(getenv("RT_RESPRP")) : 0;
|
||
if (nowarp && resrp_only) {
|
||
for (size_t k = 0; k < nbin; k++) {
|
||
double res_k = std::max(static_cast<double>(res[k]), 1e-12);
|
||
mask_out[k] *= std::pow(res_k, rp);
|
||
}
|
||
}
|
||
|
||
// Step 9 (NOTES_LEVEL:830 + consumers_out.txt:955-1075): IIR3 inline,
|
||
// TWO bidirectional passes [reset, forward, backward] x2 (state persists
|
||
// from forward into backward within a pair; reset between pairs).
|
||
// y = B3[i]*x[i] + A3[i]*state (decomp operand order verified).
|
||
static const int no_iir3 = getenv("RT_NOIIR3") ? atoi(getenv("RT_NOIIR3")) : 0;
|
||
for (int pass = 0; pass < 2 && !no_iir3; pass++) {
|
||
double st = 0.0;
|
||
for (size_t i = 0; i < nbin; i++) {
|
||
double y = static_cast<double>(mask_out[i]) * kIIR_B3[i] + st * kIIR_A3[i];
|
||
st = y;
|
||
mask_out[i] = static_cast<float>(y);
|
||
}
|
||
for (size_t i = nbin - 2; i >= 1; i--) {
|
||
double y = static_cast<double>(mask_out[i]) * kIIR_B3[i] + st * kIIR_A3[i];
|
||
st = y;
|
||
mask_out[i] = static_cast<float>(y);
|
||
}
|
||
}
|
||
for (size_t k = 0; k < half; k++) {
|
||
mask_out[nfft - 1 - k] = mask_out[k];
|
||
}
|
||
|
||
for (size_t k = 0; k < nfft; k++) {
|
||
mask_out[k] = mask_out[k] * (fVar30 * 1.0f) + (1.0f - fVar30);
|
||
}
|
||
|
||
// RT_DUMP_BIN: single-frame per-bin tract at frame RT_DUMP_FRAME (default
|
||
// 100): k am res lvl_raw band_level post-IIR1/2, pre-warp mask, W weight.
|
||
if (dbg_path && !dbg_prewarp.empty()) {
|
||
static int dbg_frames = 0;
|
||
int dbg_target = 100;
|
||
if (const char* fs = getenv("RT_DUMP_FRAME")) dbg_target = atoi(fs);
|
||
if (dbg_frames++ != dbg_target) return;
|
||
FILE* df = fopen(dbg_path, "wb");
|
||
if (df) {
|
||
fprintf(df, "# fc=%g q=%g sens=%g rp=%.6f\n", band.fc, band.q, band.sens, rp);
|
||
for (size_t k = 0; k < nbin; k++) {
|
||
double res_k = std::max(static_cast<double>(res[k]), 1e-12);
|
||
double lvl_raw = static_cast<double>(am[k]) / res_k * scale_factor;
|
||
double w = kBand768[k] * kWarp[k] * std::pow(res_k, rp);
|
||
fprintf(df, "%zu %.9g %.9g %.9g %.9g %.9g %.9g\n", k,
|
||
static_cast<double>(am[k]), res_k, lvl_raw,
|
||
static_cast<double>(band_level[k]),
|
||
static_cast<double>(dbg_prewarp[k]), w);
|
||
}
|
||
fclose(df);
|
||
}
|
||
}
|
||
|
||
// RT_DUMP_ALL trajectory: append per-frame lvl_raw spectrum (binary:
|
||
// int32 frame, int32 nbin, float32 lvl_raw[nbin]). Single-band cases only.
|
||
// Detector path is law-independent -> one capture serves offline law fits.
|
||
if (const char* ap = getenv("RT_DUMP_ALL")) {
|
||
static FILE* af = fopen(ap, "ab");
|
||
if (af) {
|
||
static int aframe = 0;
|
||
int32_t hdr[2] = {static_cast<int32_t>(aframe++),
|
||
static_cast<int32_t>(nbin)};
|
||
fwrite(hdr, sizeof(int32_t), 2, af);
|
||
for (size_t k = 0; k < nbin; k++) {
|
||
double res_k = std::max(static_cast<double>(res[k]), 1e-12);
|
||
float lv = static_cast<float>(
|
||
static_cast<double>(am[k]) / res_k * scale_factor);
|
||
fwrite(&lv, sizeof(float), 1, af);
|
||
}
|
||
fflush(af);
|
||
}
|
||
}
|
||
}
|
||
|
||
// Wrapper that allows cascade curve override for process_band_structural.
|
||
// When casc_am is non-null, it replaces the am/res level computation.
|
||
// The cascade output IS the level curve (after Haar smooth + sin-peak floor).
|
||
// We pass res=1.0 so that am/res = am (cascade already includes twin response).
|
||
static void process_band_structural_am(
|
||
const float* am,
|
||
const float* res,
|
||
const DetectorBand& band,
|
||
float* mask_out,
|
||
size_t nfft,
|
||
float sample_rate,
|
||
size_t num_bands = 1,
|
||
const float* casc_curve = nullptr,
|
||
bool use_cascade = false
|
||
) {
|
||
if (use_cascade && casc_curve) {
|
||
// Cascade curve IS the level. Pass with res=1.0 to skip am/res division.
|
||
// Create a dummy res array of all 1.0
|
||
static thread_local std::vector<float> one_res;
|
||
size_t nbin = nfft/2 + 1;
|
||
one_res.assign(nbin, 1.0f);
|
||
process_band_structural(casc_curve, one_res.data(), band, mask_out, nfft, sample_rate, num_bands);
|
||
} else {
|
||
process_band_structural(am, res, band, mask_out, nfft, sample_rate, num_bands);
|
||
}
|
||
}
|
||
|
||
} // 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();
|
||
twin_resp_complex_.clear();
|
||
cascade_states_.clear();
|
||
|
||
// RT_DUMPRESPATH=<file> (NOTES 22t): static twin-response spectra per band,
|
||
// binary {int32 band, int32 nbin, float res[nbin]} records (append).
|
||
FILE* rp_dump = nullptr;
|
||
if (const char* dp = getenv("RT_DUMPRESPATH")) rp_dump = fopen(dp, "ab");
|
||
|
||
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);
|
||
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);
|
||
// Twin gain floor per NOTES 24dd: plugin caps res ≥0.153 @fc1000q1
|
||
// Our twin 0.0069 at sens24 vs plugin 0.153 (k=22). Floor is
|
||
// content/sens-dependent; default off (canon). Opt-in via RT_TWIN_FLOOR.
|
||
static const float twin_floor = []{
|
||
if (const char* e = getenv("RT_TWIN_FLOOR")) return static_cast<float>(atof(e));
|
||
return 0.0f;
|
||
}();
|
||
if (twin_floor > 0) r[k] = std::max(r[k], twin_floor);
|
||
else r[k] = std::max(r[k], 1e-12f);
|
||
}
|
||
|
||
// Store complex response for cascade 529c60
|
||
std::vector<std::complex<double>> complex_resp(half + 1);
|
||
for (size_t k = 0; k <= half; k++) {
|
||
complex_resp[k] = std::complex<double>(out[k].re, out[k].im);
|
||
}
|
||
twin_resp_complex_.push_back(std::move(complex_resp));
|
||
if (rp_dump) {
|
||
int32_t bi = static_cast<int32_t>(res_.size());
|
||
int32_t nb = static_cast<int32_t>(r.size());
|
||
fwrite(&bi, sizeof(int32_t), 1, rp_dump);
|
||
fwrite(&nb, sizeof(int32_t), 1, rp_dump);
|
||
fwrite(r.data(), sizeof(float), r.size(), rp_dump);
|
||
}
|
||
res_.push_back(std::move(r));
|
||
}
|
||
if (rp_dump) fclose(rp_dump);
|
||
track_.assign(bands_.size(), std::vector<float>(half + 1, 1.0f));
|
||
cascade_states_.assign(bands_.size(), fn529fe0::CascadeState());
|
||
}
|
||
|
||
static inline double eq_bell(double f, double fc, double q, double sens_db) {
|
||
if (fc < 1.0 || sens_db == 0) return 1.0;
|
||
// 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) {
|
||
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_));
|
||
|
||
// RT_ENV=live (NOTES 22n): detector envelope from live tables kRTAtt/kRTRel —
|
||
// attack as feed-forward, release as retention (~tau 2s at hop rate). This is
|
||
// the slow adaptation the real plugin exhibits on sustained content.
|
||
static const int env_live = getenv("RT_ENV") ? atoi(getenv("RT_ENV")) : 0;
|
||
|
||
for (size_t k = 0; k <= half; k++) {
|
||
double a_cur = 2.0 * std::abs(spectrum[k]) / wsum_;
|
||
if (env_live) {
|
||
double d = a_cur - static_cast<double>(am_[k]);
|
||
if (d > 0) am_[k] = static_cast<float>(am_[k] + d * static_cast<double>(kRTAtt[k]));
|
||
else am_[k] = static_cast<float>(static_cast<double>(kRTRel[k]) * am_[k]
|
||
+ (1.0 - static_cast<double>(kRTRel[k])) * a_cur);
|
||
} else {
|
||
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);
|
||
}
|
||
}
|
||
// EQ before detector moved to lvl calc in process_band_structural (not am_ state)
|
||
|
||
// Detector cascade 529c60: per-band pre-processor on complex twin-filtered
|
||
// spectrum. Computes magnitudes, Haar-smooths, applies sin-peak floor.
|
||
static const int casc_on = getenv("RT_CASC") ? atoi(getenv("RT_CASC")) : 0;
|
||
|
||
for (size_t k = 0; k <= half; k++) mask[k] = 1.0f;
|
||
|
||
if (is_internal_grid(nfft_, sample_rate_)) {
|
||
// RT_FAITHFUL=1 (NOTES 22w): BLOCKMAP_529fe0 transcription path
|
||
static const int faithful = getenv("RT_FAITHFUL") ? atoi(getenv("RT_FAITHFUL")) : 0;
|
||
static const fnfaith::Params fparams = faithful ? fnfaith::params_from_env()
|
||
: fnfaith::Params{};
|
||
// same scale_factor as process_band_structural (line ~79)
|
||
constexpr float sf = 15.0f * 440.95f / 2048.0f;
|
||
for (size_t b = 0; b < bands_.size(); b++) {
|
||
std::vector<float> band_mask(nfft_, 1.0f);
|
||
if (faithful) {
|
||
fnfaith::band_mask_faithful(am_.data(), res_[b].data(), half + 1,
|
||
sample_rate_, sf, fparams,
|
||
band_mask.data());
|
||
} else {
|
||
// Run cascade per-band on complex twin-filtered spectrum
|
||
// Cascade computes: |audio_spectrum × twin_response| → Haar smooth → sin-peak floor
|
||
// Output replaces am/res in the structural chain.
|
||
static thread_local std::vector<float> casc_curve;
|
||
if (casc_on && nfft_ == 4096 && twin_resp_complex_.size() > b) {
|
||
size_t nbin = half + 1;
|
||
std::vector<float> complex_input(2 * nbin);
|
||
casc_curve.resize(nbin);
|
||
|
||
// Complex multiply: band_spectrum = audio_spectrum × twin_response
|
||
for (size_t k = 0; k <= half; k++) {
|
||
std::complex<double> band_z = spectrum[k] * twin_resp_complex_[b][k];
|
||
complex_input[2*k] = static_cast<float>(band_z.real());
|
||
complex_input[2*k+1] = static_cast<float>(band_z.imag());
|
||
}
|
||
|
||
fn529fe0::cascade_detect(
|
||
complex_input.data(),
|
||
casc_curve.data(),
|
||
cascade_states_[b],
|
||
nbin,
|
||
2, // Haar iterations
|
||
0.0f, // sin_peak_param (0 = no floor; set >0 for Step 9 floor)
|
||
48000.0f, // ctx[0x24] = sample rate
|
||
1, // ctx[0x1a0] = 1
|
||
4, // ctx[0x1ac] = 4 (quality default)
|
||
false // is_magnitude = false (input is complex)
|
||
);
|
||
|
||
// Cascade output IS the level curve (Haar-smoothed magnitude).
|
||
// Use it directly as am_ replacement — pass res=1.0 so level = am*1
|
||
// (twin response already baked into cascade output).
|
||
process_band_structural_am(am_.data(), res_[b].data(), bands_[b],
|
||
band_mask.data(), nfft_, sample_rate_,
|
||
bands_.size(), casc_curve.data(), true);
|
||
// 1c chain_9_19 wire (BLOCKMAP:620) — gated by RT_CASC=1+chain, default off
|
||
// B: real warp/att/rel tables (was unity proxy)
|
||
if (getenv("RT_CASC_CHAIN")) {
|
||
static std::vector<float> tmp6f8, accVec;
|
||
size_t nbin = half + 1;
|
||
tmp6f8.assign(nbin, 0.0f);
|
||
accVec.assign(nbin, 0.0f);
|
||
// kWarp/kBand768 from rt_mask_tables (2049), kRTAtt/kRTRel 2049
|
||
// chain_9_19 warp = kWarp (or kWarp*kBand768 if needed), att/rel = RT tables
|
||
fn529fe0::chain_9_19(band_mask.data(), tmp6f8.data(), accVec.data(),
|
||
kWarp, kRTAtt, kRTRel, nbin);
|
||
}
|
||
} else {
|
||
process_band_structural(am_.data(), res_[b].data(), bands_[b],
|
||
band_mask.data(), nfft_, sample_rate_, bands_.size());
|
||
}
|
||
}
|
||
for (size_t k = 0; k <= half; k++) {
|
||
mask[k] = std::min(band_mask[k], mask[k]);
|
||
}
|
||
}
|
||
} else {
|
||
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_pchip(xv) + 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];
|
||
}
|
||
}
|