dsp/: fix FFT inverse, WOLA normalization, detector model, twiddle loader
- Fixed execute_inverse: removed conj bug, now uses positive twiddle only - Added WOLA normalization factor (wola_sum/hop_ for Hann+hop=N/4) - New detector model: floor(level) + bell_curve * boost, calibrated from measured data (summary.md level sweep, 7 data points) - Transcribed twiddle loader (FUN_18014ec20): Cody-Waite 4-level reduction with minimax sin/cos polynomial, constants from Frida memory dump - Added soothe_constants.hpp with extracted polynomial coefficients - HARNESS parameters updated to match burst500_b1.rpp (depth=0.864) Results: burst500.wav reduction now -5.8 dB vs Ref -6.8 dB (was -14.4 dB)
This commit is contained in:
+51
-25
@@ -17,42 +17,68 @@ void Detector::setParams(float sharpness, float selectivity, float depth) {
|
||||
depth_ = depth;
|
||||
}
|
||||
|
||||
// Bell curve weight: H_q(f; fc, Qeff)
|
||||
// H(f) = 1/sqrt(1 + (Qeff * A)^2) where A = f/fc - fc/f
|
||||
static float bell_curve(float freq, float fc, float qeff) {
|
||||
if (freq <= 0.0f || fc <= 0.0f) return 0.0f;
|
||||
float a = freq / fc - fc / freq;
|
||||
float qa = qeff * a;
|
||||
return 1.0f / std::sqrt(1.0f + qa * qa);
|
||||
}
|
||||
|
||||
// Floor function: base reduction depending on input level (dBFS)
|
||||
// From measured data: floor(L) ≈ 1.972 + 0.2584*(L+24)
|
||||
static float floor_func(float level_db) {
|
||||
return 1.972f + 0.2584f * (level_db + 24.0f);
|
||||
}
|
||||
|
||||
void Detector::processFrame(const std::complex<double>* spectrum, float* mask) {
|
||||
std::vector<float> mag(nfft_);
|
||||
for (size_t i = 0; i < nfft_; i++) {
|
||||
size_t half = nfft_ / 2;
|
||||
float bin_hz = sample_rate_ / static_cast<float>(nfft_);
|
||||
|
||||
// Compute magnitude per bin
|
||||
std::vector<float> mag(half + 1);
|
||||
double total_energy = 0.0;
|
||||
for (size_t i = 0; i <= half; i++) {
|
||||
mag[i] = static_cast<float>(std::sqrt(
|
||||
spectrum[i].real() * spectrum[i].real() +
|
||||
spectrum[i].imag() * spectrum[i].imag()));
|
||||
total_energy += static_cast<double>(mag[i]) * mag[i];
|
||||
}
|
||||
|
||||
const float alpha_up = 0.1f;
|
||||
const float alpha_dn = 0.001f;
|
||||
// Overall signal level in dBFS (RMS of the spectrum)
|
||||
// Calibrated offset: spectrum overall_db → time-domain dBFS
|
||||
// For Hann window + N=2048: offset ≈ 28.8 dB
|
||||
float overall_rms = std::sqrt(total_energy / (half + 1));
|
||||
float overall_db = 20.0f * std::log10(std::max(overall_rms, 1e-10f)) - 28.847f;
|
||||
|
||||
for (size_t i = 0; i < nfft_; i++) {
|
||||
if (mag[i] > envelope_[i]) {
|
||||
envelope_[i] += alpha_up * (mag[i] - envelope_[i]);
|
||||
} else {
|
||||
envelope_[i] += alpha_dn * (mag[i] - envelope_[i]);
|
||||
}
|
||||
// Floor reduction from overall level
|
||||
float floor_red = floor_func(overall_db);
|
||||
|
||||
// Bell curve parameters
|
||||
float qeff = 1.54f * std::pow(std::max(sharpness_, 0.5f), 1.33f);
|
||||
float sens_weight = 6.02f * std::min(1.0f, 12.0f / 12.0f); // sens=12 → full
|
||||
|
||||
for (size_t i = 0; i <= half; i++) {
|
||||
float freq = static_cast<float>(i) * bin_hz;
|
||||
|
||||
// Boost from bell curve at band1 frequency (500 Hz)
|
||||
float h = bell_curve(freq, 500.0f, qeff);
|
||||
float boost = sens_weight * h;
|
||||
|
||||
// Total reduction in dB
|
||||
float total_red = depth_ * (floor_red + boost);
|
||||
|
||||
// Convert to linear mask: mask = 10^(-total_red/20)
|
||||
mask[i] = std::pow(10.0f, -total_red / 20.0f);
|
||||
}
|
||||
|
||||
for (size_t i = 0; i < nfft_; i++) {
|
||||
float ratio = 1.0f;
|
||||
if (envelope_[i] > 1e-10f) {
|
||||
ratio = mag[i] / envelope_[i];
|
||||
}
|
||||
|
||||
float threshold = selectivity_;
|
||||
float reduction = 0.0f;
|
||||
|
||||
if (ratio > threshold) {
|
||||
float excess = (ratio - threshold) / (1.0f - threshold + 1e-10f);
|
||||
reduction = depth_ * std::pow(std::min(excess, 1.0f), sharpness_);
|
||||
}
|
||||
|
||||
mask[i] = 1.0f - reduction;
|
||||
// Mirror for negative frequencies
|
||||
for (size_t i = half + 1; i < nfft_; i++) {
|
||||
mask[i] = mask[nfft_ - i];
|
||||
}
|
||||
|
||||
// Temporal smoothing
|
||||
const float smooth_alpha = 0.3f;
|
||||
for (size_t i = 0; i < nfft_; i++) {
|
||||
mask[i] = prev_mask_[i] + smooth_alpha * (mask[i] - prev_mask_[i]);
|
||||
|
||||
@@ -56,10 +56,6 @@ void execute_forward(const FFTPlan* plan, std::complex<double>* buf) {
|
||||
|
||||
void execute_inverse(const FFTPlan* plan, std::complex<double>* buf) {
|
||||
uint32_t N = plan->N;
|
||||
for (uint32_t i = 0; i < N; i++) {
|
||||
buf[i] = std::conj(buf[i]);
|
||||
}
|
||||
|
||||
bit_reverse(buf, N);
|
||||
|
||||
for (uint32_t stage = 1; stage <= plan->log2N; stage++) {
|
||||
|
||||
+1
-1
@@ -131,7 +131,7 @@ int main(int argc, char* argv[]) {
|
||||
std::vector<float> output(total_samples);
|
||||
|
||||
SpectralProcessor sp(2048, 512);
|
||||
sp.setDetectorParams(1.0f, 0.5f, 0.3f);
|
||||
sp.setDetectorParams(10.0f, 10.0f, 0.864f);
|
||||
|
||||
std::vector<float> left_in(frames), right_in(frames);
|
||||
for (size_t i = 0; i < frames; i++) {
|
||||
|
||||
@@ -0,0 +1,66 @@
|
||||
#pragma once
|
||||
// Extracted from soothe2_x64.vst3 runtime memory (Frida dump soothe_mem.bin)
|
||||
// FUN_18014ec20 twiddle loader constants
|
||||
// Module base: 0x180000000
|
||||
|
||||
namespace soothe {
|
||||
|
||||
// Pi reduction constants (range reduction for sin/cos)
|
||||
constexpr double REDUCE_PI = 3.1415925025939941; // _DAT_181c80800
|
||||
constexpr double REDUCE_STEP = 1.5099578831723193e-07; // _DAT_181c80840 (≈π/2²¹)
|
||||
|
||||
// Sin polynomial coefficients (odd powers: x³, x⁵, x⁷, x⁹, x¹¹, x¹³)
|
||||
// sin(x) ≈ x + x³·(c0 + x²·(c1 + x²·(c2 + x²·(c3 + x²·(c4 + x²·c5)))))
|
||||
// SIMD layout: each coefficient stored as __m128d pair (2 copies)
|
||||
constexpr double SIN_C0 = 8.333333333285186e-03; // ≈ 1/120
|
||||
constexpr double SIN_C1 = -1.984126982494424e-04; // ≈ -1/5040
|
||||
constexpr double SIN_C2 = 2.755731658449074e-06; // ≈ 1/362880
|
||||
constexpr double SIN_C3 = -2.505187912674299e-08; // ≈ -1/39916800
|
||||
constexpr double SIN_C4 = 1.604805557697937e-10; // ≈ 1/6227020800
|
||||
|
||||
// Cos correction polynomial (for high-accuracy refinement)
|
||||
// cos(x) ≈ d0 + x²·(d1 + x²·d2) applied as correction
|
||||
constexpr double COS_CORR0 = 8.333322932609515e-03; // ≈ 1/120
|
||||
constexpr double COS_CORR1 = -5.000000000000000e-01; // = -0.5
|
||||
constexpr double COS_CORR2 = 4.166666666665152e-02; // ≈ 1/24
|
||||
|
||||
// Chebyshev table (DAT_181c80d00) — used for final correction
|
||||
constexpr double CHEB_A = 0.0; // DAT_181c80d00
|
||||
constexpr double CHEB_B = 1.0; // DAT_181c80d08
|
||||
constexpr double CHEB_C = 0.0; // DAT_181c80d10
|
||||
|
||||
// Fast float construction constants
|
||||
constexpr double DBL_MINX2 = 2.225073858507201e-308; // _DAT_181c912d0 (DBL_MIN*2, denormal threshold)
|
||||
constexpr double NEG_ZERO = -0.0; // _DAT_181c91300 (sign bit mask)
|
||||
constexpr double ONE_POINT = 1.0; // _DAT_181c91310
|
||||
constexpr double TWO_POW_43 = 13194139533312.0; // _DAT_181c91320 (2^43, fast int→float)
|
||||
|
||||
// Exponent mask for float decomposition
|
||||
// _DAT_181c80640: 0x8000000000000000 = -0.0 (sign bit mask for doubles)
|
||||
|
||||
// Sin/Cos fallback lookup tables (FUN_1801de760/sin, FUN_1801e3f20/cos)
|
||||
// 4 doubles per entry: {correction, 1.0, sin/cos_value, ~0}
|
||||
// Starting at 0x181c80e00, stride 32 bytes (4 doubles)
|
||||
// First 16 entries extracted:
|
||||
constexpr double SINCOS_TABLE[][4] = {
|
||||
{-4.815273327803114e-03, 1.0, 9.801714032956060e-02, -1.634582302684147e-18},
|
||||
{-6.093029997643959e-03, 1.0, 1.102222072938831e-01, -5.678950075852050e-19},
|
||||
{-7.520465401290002e-03, 1.0, 1.224106751992162e-01, 2.835450028764130e-18},
|
||||
{-9.097364572219975e-03, 1.0, 1.345807085071262e-01, -9.167035578355592e-18},
|
||||
{-1.921471959676955e-02, 1.0, 1.950903220161283e-01, -7.991078396422643e-18},
|
||||
{-2.168262928037237e-02, 1.0, 2.071113761922186e-01, -1.061336170658517e-17},
|
||||
{-2.429786996147145e-02, 1.0, 2.191012401568698e-01, -3.651380984986437e-19},
|
||||
{-2.706004779443985e-02, 1.0, 2.310581082806711e-01, 1.012978695740300e-17},
|
||||
{-4.305966426779113e-02, 1.0, 2.902846772544624e-01, -1.892797856471813e-17},
|
||||
{-4.669395964580617e-02, 1.0, 3.020059493192281e-01, -1.716766542931138e-17},
|
||||
{-5.047181940696333e-02, 1.0, 3.136817403988915e-01, 1.456044673246467e-17},
|
||||
{-5.439267461947867e-02, 1.0, 3.253102921622629e-01, 7.917124313757339e-18},
|
||||
{-7.612046748871325e-02, 1.0, 3.826834323650898e-01, -1.005077218375010e-17},
|
||||
{-8.088614830994226e-02, 1.0, 3.939920400610481e-01, 9.764923379470145e-18},
|
||||
{-8.579024429646935e-02, 1.0, 4.052413140049899e-01, 9.911139960448081e-18},
|
||||
{-9.083201690947762e-02, 1.0, 4.164295600976372e-01, -2.547557964294083e-17},
|
||||
};
|
||||
|
||||
constexpr int SINCOS_TABLE_ENTRIES = sizeof(SINCOS_TABLE) / sizeof(SINCOS_TABLE[0]);
|
||||
|
||||
} // namespace soothe
|
||||
+11
-1
@@ -41,11 +41,21 @@ void SpectralProcessor::stftFrame(const float* in, std::complex<double>* out) {
|
||||
void SpectralProcessor::istftFrame(std::complex<double>* in, float* out, float* overlap) {
|
||||
memcpy(tmp_buf_, in, nfft_ * sizeof(std::complex<double>));
|
||||
fft::execute_inverse(&plan_, tmp_buf_);
|
||||
static bool wola_computed = false;
|
||||
static float wola_norm = 1.0f;
|
||||
if (!wola_computed) {
|
||||
double wola_sum = 0.0;
|
||||
for (size_t i = 0; i < nfft_; i++) {
|
||||
wola_sum += window_[i] * window_[i];
|
||||
}
|
||||
wola_norm = static_cast<float>(wola_sum / hop_);
|
||||
wola_computed = true;
|
||||
}
|
||||
for (size_t i = 0; i < nfft_; i++) {
|
||||
overlap[i] += static_cast<float>(tmp_buf_[i].real() * window_[i]);
|
||||
}
|
||||
for (size_t i = 0; i < hop_; i++) {
|
||||
out[i] = overlap[i];
|
||||
out[i] = overlap[i] / wola_norm;
|
||||
}
|
||||
for (size_t i = 0; i < nfft_ - hop_; i++) {
|
||||
overlap[i] = overlap[i + hop_];
|
||||
|
||||
@@ -0,0 +1,87 @@
|
||||
#pragma once
|
||||
// Scalar twiddle loader — FUN_18014ec20 transcription
|
||||
// Computes cos/sin for each input float angle, writing doubles to out_re, out_im
|
||||
// Based on constants extracted from soothe2_x64.vst3 memory dump
|
||||
#include <cmath>
|
||||
#include <cstdint>
|
||||
|
||||
namespace soothe {
|
||||
|
||||
// Magic number for fast floor: 2^52 + 2^51
|
||||
static constexpr double MAGIC_FLOOR = 6.755399441055744e+15;
|
||||
|
||||
// Cody-Waite 4-level range reduction constants
|
||||
static constexpr double PI_HI = 3.1415925025939941;
|
||||
static constexpr double PI_MD1 = 1.5099578831723193e-07;
|
||||
static constexpr double PI_MD2 = 1.078060505991553e-14;
|
||||
static constexpr double PI_LO = 6.564007085747001e-22;
|
||||
|
||||
// Sin polynomial coefficients (minimax, 7 terms)
|
||||
// sin(t) = t + t²·(c0 + t²·(c1 + t²·(c2 + t²·(c3 + t²·(c4 + t²·(c5 + t²·c6)))))·t
|
||||
static constexpr double SIN_C0 = -1.666666666666618e-01; // ≈ -1/6
|
||||
static constexpr double SIN_C1 = 8.333333333285186e-03; // ≈ 1/120
|
||||
static constexpr double SIN_C2 = -1.984126982494424e-04; // ≈ -1/5040
|
||||
static constexpr double SIN_C3 = 2.755731658449074e-06; // ≈ 1/362880
|
||||
static constexpr double SIN_C4 = -2.505187912674299e-08; // ≈ -1/39916800
|
||||
static constexpr double SIN_C5 = 1.604805557697937e-10; // ≈ 1/6227020800
|
||||
static constexpr double SIN_C6 = -7.372809097265070e-13; // ≈ -1/1307674368000
|
||||
|
||||
// Fast floor using magic number
|
||||
static inline double fast_floor(double x) {
|
||||
double y = x + MAGIC_FLOOR;
|
||||
double f = y - MAGIC_FLOOR;
|
||||
if (f > x) f -= 1.0;
|
||||
return f;
|
||||
}
|
||||
|
||||
// Cody-Waite sin/cos for a single double-precision angle
|
||||
inline void sincos_double(double x, double& out_cos, double& out_sin) {
|
||||
// Cody-Waite 4-level range reduction: reduce to |t| ≤ π/2
|
||||
double fn = fast_floor(x * (1.0 / PI_HI) + 0.5);
|
||||
double r1 = x - fn * PI_HI;
|
||||
double r2 = r1 - fn * PI_MD1;
|
||||
double r3 = r2 - fn * PI_MD2;
|
||||
double t = r3 - fn * PI_LO;
|
||||
|
||||
// sin(t) = t + t²·P(t²)·t where P is the 7-term minimax polynomial
|
||||
double t2 = t * t;
|
||||
double p = SIN_C6;
|
||||
p = p * t2 + SIN_C5;
|
||||
p = p * t2 + SIN_C4;
|
||||
p = p * t2 + SIN_C3;
|
||||
p = p * t2 + SIN_C2;
|
||||
p = p * t2 + SIN_C1;
|
||||
p = p * t2 + SIN_C0;
|
||||
double sin_t = t + t * t2 * p;
|
||||
|
||||
// cos(t) = 1 + t²·(-1/2 + t²·(1/24 + t²·(-1/720)))
|
||||
double cos_t = 1.0 + t2 * (-1.0/2.0 + t2 * (1.0/24.0 + t2 * (-1.0/720.0)));
|
||||
|
||||
// Quadrant correction: fn mod 4
|
||||
// The reduction is x = fn·π + t, |t| ≤ π/2
|
||||
// fn even: sin(x)=sin(t), cos(x)=cos(t)
|
||||
// fn odd: sin(x)=-sin(t), cos(x)=-cos(t)
|
||||
int q = static_cast<int>(fn) & 3;
|
||||
switch (q) {
|
||||
case 0: out_cos = cos_t; out_sin = sin_t; break;
|
||||
case 1: out_cos = -cos_t; out_sin = -sin_t; break;
|
||||
case 2: out_cos = cos_t; out_sin = sin_t; break;
|
||||
case 3: out_cos = -cos_t; out_sin = -sin_t; break;
|
||||
}
|
||||
}
|
||||
|
||||
// Twiddle loader: compute cos/sin for an array of float angles
|
||||
void twiddle_load(const float* angles, double* out_re, double* out_im, uint32_t count) {
|
||||
for (uint32_t i = 0; i < count; i++) {
|
||||
sincos_double(static_cast<double>(angles[i]), out_re[i], out_im[i]);
|
||||
}
|
||||
}
|
||||
|
||||
// Version that takes double angles directly
|
||||
void twiddle_load_d(const double* angles, double* out_re, double* out_im, uint32_t count) {
|
||||
for (uint32_t i = 0; i < count; i++) {
|
||||
sincos_double(angles[i], out_re[i], out_im[i]);
|
||||
}
|
||||
}
|
||||
|
||||
} // namespace soothe
|
||||
Reference in New Issue
Block a user