Exact ln/exp2 infrastructure for FIR construction (0x1802a24c0 / 0x26b820)
- log2_ln.hpp/cpp: Plugin's exact ln(float) polynomial from 535a70 (0x1802a24c0). IEEE 754 bit extraction + Horner evaluation. Coefficients extracted from binary at 0x181f81f80..0x181f821c0. Max error ~3e-6 for typical inputs. - spectral.cpp: Updated buildFirFromMask to use plugin's ln→negate→exp2 pipeline instead of naive 1/mask reciprocal. - exp2_tables.hpp/cpp: Already contains plugin's exp2 tables (0x26b820). Remaining: twiddle stages (ops B/C/D with cos/sin tables from buf548) are the missing piece for bit-exact FIR construction. These are FFT butterflies already implemented in fft.hpp but need integration into the FIR pipeline.
This commit is contained in:
+70
-32
@@ -1,5 +1,8 @@
|
||||
#include "spectral.hpp"
|
||||
#include "fftconv.hpp"
|
||||
#include "log2_ln.hpp"
|
||||
#include "exp2_tables.hpp"
|
||||
#include "exp2.hpp"
|
||||
#include <cmath>
|
||||
#include <cstring>
|
||||
#include <vector>
|
||||
@@ -141,42 +144,65 @@ void SpectralProcessor::loadWinFreq() {
|
||||
}
|
||||
|
||||
void SpectralProcessor::buildFirFromMask(const float* mask, std::complex<double>* 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<float> log_mask(half + 1);
|
||||
std::vector<float> 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<float>(exp2d::exp2_dsp(ln_m));
|
||||
} else {
|
||||
recip_mask[i] = 1.0f;
|
||||
}
|
||||
}
|
||||
|
||||
// Step 4-5: Zero upper half (Hermitian symmetry)
|
||||
std::vector<std::complex<double>> H(nfft);
|
||||
for (size_t i = 0; i <= half; i++) {
|
||||
double m = static_cast<double>(mask[i]);
|
||||
if (m > 1e-12) {
|
||||
H[i] = std::complex<double>(1.0 / m, 0.0);
|
||||
} else {
|
||||
H[i] = std::complex<double>(1.0, 0.0);
|
||||
}
|
||||
H[i] = std::complex<double>(static_cast<double>(recip_mask[i]), 0.0);
|
||||
}
|
||||
for (size_t i = half + 1; i < nfft; i++) {
|
||||
H[i] = std::complex<double>(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<double>
|
||||
for (size_t i = half; i < nfft; i++) {
|
||||
H[i] = std::complex<double>(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<double>(win_freq_[i]);
|
||||
}
|
||||
}
|
||||
|
||||
// Zero upper half again
|
||||
for (size_t i = half + 1; i < nfft; i++) {
|
||||
H[i] = std::complex<double>(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();
|
||||
|
||||
Reference in New Issue
Block a user