diff --git a/dsp/detect.cpp b/dsp/detect.cpp index a7f48fb..270ce3e 100644 --- a/dsp/detect.cpp +++ b/dsp/detect.cpp @@ -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* spectrum, float* mask) { - std::vector mag(nfft_); - for (size_t i = 0; i < nfft_; i++) { + size_t half = nfft_ / 2; + float bin_hz = sample_rate_ / static_cast(nfft_); + + // Compute magnitude per bin + std::vector mag(half + 1); + double total_energy = 0.0; + for (size_t i = 0; i <= half; i++) { mag[i] = static_cast(std::sqrt( spectrum[i].real() * spectrum[i].real() + spectrum[i].imag() * spectrum[i].imag())); + total_energy += static_cast(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(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]); diff --git a/dsp/fft.cpp b/dsp/fft.cpp index a26cf88..ece48d4 100644 --- a/dsp/fft.cpp +++ b/dsp/fft.cpp @@ -56,10 +56,6 @@ void execute_forward(const FFTPlan* plan, std::complex* buf) { void execute_inverse(const FFTPlan* plan, std::complex* 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++) { diff --git a/dsp/harness.cpp b/dsp/harness.cpp index cfbd900..f96b897 100644 --- a/dsp/harness.cpp +++ b/dsp/harness.cpp @@ -131,7 +131,7 @@ int main(int argc, char* argv[]) { std::vector 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 left_in(frames), right_in(frames); for (size_t i = 0; i < frames; i++) { diff --git a/dsp/soothe_constants.hpp b/dsp/soothe_constants.hpp new file mode 100644 index 0000000..77b0f72 --- /dev/null +++ b/dsp/soothe_constants.hpp @@ -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 diff --git a/dsp/spectral.cpp b/dsp/spectral.cpp index f1128d0..08255e2 100644 --- a/dsp/spectral.cpp +++ b/dsp/spectral.cpp @@ -41,11 +41,21 @@ void SpectralProcessor::stftFrame(const float* in, std::complex* out) { void SpectralProcessor::istftFrame(std::complex* in, float* out, float* overlap) { memcpy(tmp_buf_, in, nfft_ * sizeof(std::complex)); 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(wola_sum / hop_); + wola_computed = true; + } for (size_t i = 0; i < nfft_; i++) { overlap[i] += static_cast(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_]; diff --git a/dsp/twiddle_loader.hpp b/dsp/twiddle_loader.hpp new file mode 100644 index 0000000..6bbfc25 --- /dev/null +++ b/dsp/twiddle_loader.hpp @@ -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 +#include + +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(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(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