diff --git a/dsp/CMakeLists.txt b/dsp/CMakeLists.txt index a80bdb2..de21559 100644 --- a/dsp/CMakeLists.txt +++ b/dsp/CMakeLists.txt @@ -21,6 +21,7 @@ add_library(soothe2_dsp SHARED fftconv.cpp vlog.cpp leveltrack.cpp + framed_model.cpp ) add_executable(soothe2_harness harness.cpp) diff --git a/dsp/framed_model.cpp b/dsp/framed_model.cpp new file mode 100644 index 0000000..a3b1dc5 --- /dev/null +++ b/dsp/framed_model.cpp @@ -0,0 +1,136 @@ +#include "framed_model.hpp" +#include "twin.hpp" +#include "freqpath.hpp" +#include +#include +#include + +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(sample_rate_), static_cast(fc), + static_cast(q_), sens_lin); + std::vector z(half + 1); + std::vector out(half + 1); + for (size_t k = 0; k <= half; k++) { + double theta = 2.0 * M_PI * static_cast(k) / static_cast(nfft_); + z[k].re = static_cast(std::cos(theta)); + z[k].im = static_cast(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(sample_rate_), + static_cast(nfft_), warp_.data()); + + inited_ = true; +} + +void FramedDetector::processFrame(const std::complex* 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(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(am); + + double xv = std::log10(std::max(am / static_cast(res_[k]), 1e-9)); + double C = G_FIT * lut_eval(xv) + + W_FIT * std::pow(static_cast(warp_[k]), A_FIT); + double gain = std::max(1.0 - C, 1e-9) * + std::pow(static_cast(res_[k]), rp); + mask[k] = static_cast(gain); + } + for (size_t k = half + 1; k < nfft_; k++) { + mask[k] = mask[nfft_ - k]; + } +} diff --git a/dsp/framed_model.hpp b/dsp/framed_model.hpp new file mode 100644 index 0000000..b23c4bc --- /dev/null +++ b/dsp/framed_model.hpp @@ -0,0 +1,36 @@ +#pragma once +#include +#include +#include + +// FramedDetector — C++ port of framed_render.py (Phase 5 pilot, mean=0.175 dB). +// Per-frame per-bin mask: +// res[f] = |2B(z)/A(z)| (twin resonance, detkernel::twin) +// am = 2|X_k|/wsum (per-bin input amplitude, smoothed attack/release) +// xv = log10(am / res) +// C = G_FIT*LUT(xv) + W_FIT*warp(f)^A_FIT +// gain = max(1-C, eps) * res^(rp0*Q^drp) +class FramedDetector { +public: + FramedDetector(size_t nfft, float sample_rate); + ~FramedDetector(); + + // fc = band center freq, q = resonance Q (>=1). Precomputes res[] + warp[]. + void setParams(float fc, float q); + + // spectrum = forward FFT of one windowed frame (nfft/2+1 bins used). + // mask output = per-bin complex gain (nfft entries, mirror applied). + void processFrame(const std::complex* spectrum, float* mask); + +private: + size_t nfft_; + float sample_rate_; + float fc_; + float q_; + double wsum_; + + std::vector res_; // per-bin |2B/A| + std::vector warp_; // per-bin warp tilt + std::vector am_; // smoothed per-bin amplitude state + bool inited_; +}; diff --git a/dsp/harness.cpp b/dsp/harness.cpp index a4ae922..db5f058 100644 --- a/dsp/harness.cpp +++ b/dsp/harness.cpp @@ -156,10 +156,18 @@ int main(int argc, char* argv[]) { } SpectralProcessor sp(2048, 512); - sp.setDetectorParams( - static_cast(params.sharpness), - static_cast(params.selectivity), - static_cast(params.depth)); + float fc = 500.0f; + float q = 1.0f; + if (!params.bands.empty()) { + for (const auto& b : params.bands) { + if (b.on > 0.5f && b.freq > 1.0f) { + fc = static_cast(b.freq); + q = static_cast(b.q > 0.0f ? b.q : 1.0f); + break; + } + } + } + sp.setDetectorParams(fc, q); std::vector left(frames, 0.0f), right(frames, 0.0f); encode_ms(left_in.data(), right_in.data(), frames); diff --git a/dsp/spectral.cpp b/dsp/spectral.cpp index 08255e2..7632b86 100644 --- a/dsp/spectral.cpp +++ b/dsp/spectral.cpp @@ -21,8 +21,8 @@ SpectralProcessor::~SpectralProcessor() { delete[] tmp_buf_; } -void SpectralProcessor::setDetectorParams(float sharpness, float selectivity, float depth) { - detector_.setParams(sharpness, selectivity, depth); +void SpectralProcessor::setDetectorParams(float fc, float q) { + detector_.setParams(fc, q); } void SpectralProcessor::computeWindow() { diff --git a/dsp/spectral.hpp b/dsp/spectral.hpp index 1dbb941..18d0aa7 100644 --- a/dsp/spectral.hpp +++ b/dsp/spectral.hpp @@ -4,7 +4,7 @@ #include #include #include "fft.hpp" -#include "detect.hpp" +#include "framed_model.hpp" constexpr size_t DEFAULT_NFFT = 2048; constexpr size_t DEFAULT_HOP = 512; @@ -14,7 +14,7 @@ public: SpectralProcessor(size_t nfft = DEFAULT_NFFT, size_t hop = DEFAULT_HOP); ~SpectralProcessor(); - void setDetectorParams(float sharpness, float selectivity, float depth); + void setDetectorParams(float fc, float q); void processBlock(float* in, float* out, size_t num_samples, size_t num_channels = 1); private: @@ -26,7 +26,7 @@ private: std::complex* tmp_buf_; std::vector overlap_; std::vector mask_; - Detector detector_; + FramedDetector detector_; size_t frame_count_; size_t output_pos_;