Compare commits

...
3 Commits
Author SHA1 Message Date
Matiq f689023089 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.
2026-08-27 12:18:57 +03:00
Matiq 575d26a771 STFT partitioned conv: FIR construction pipeline (52b550-52b8bb)
Implement minimum-phase FIR design from BLOCKMAP:
- buildFirFromMask: mask → 1/mask (reciprocal via log→negate→exp) →
  IFFT → causal window → FFT → normalize → complex multiply
- RT_FIRCONV=2 activates the new path
- RT_FIRCONV=1 preserved as simple mask × audio (legacy)

Results (tone1kq single band):
  default (pointwise):  500Hz=-25.35 dB, 1kHz=-50.60 dB
  FIRCONV=1 (mask mul): 500Hz=-1.31 dB, 1kHz=-25.92 dB
  FIRCONV=2 (min-phase): same as FIRCONV=1

The twiddle stages (ops B/C/D with cos/sin tables) are the missing
piece for bit-exact FIR construction. They perform FMA operations
with twiddle factors that modify the mask shape.

Note: dual_b1q_0.5.wav reference is empty (0 bytes) — corpus can't run.
Needs regeneration.
2026-08-27 06:13:18 +03:00
Matiq 4d78785f0b cascade integration: sin-peak floor, complex twin resp storage, per-band cascade
- Add sin-peak floor mechanism (529c60): RT_CASC_SINPEAK param
  Formula: sin_peak = sin(param*30-90) * 0.115129 * peak_level
  Floor active for param in [3,9], max at param=6 (ln10/20=0.115129)
  Prevents over-reduction by clamping level curve from below

- Store complex twin filter responses in FramedDetector::setParams()
  for cascade 529c60 per-band processing

- Add cascade state persistence (fn529fe0::CascadeState per band)

- ctx[0x24] = 48000 (sample rate, from commit 0e90918)
  With init values ctx[0x1a0]=1, ctx[0x1ac]=4, cascade w=0 (passthrough)

- All tests pass: fn529fe0_check, render48k build OK
2026-08-27 03:12:56 +03:00
7 changed files with 317 additions and 8 deletions
+1
View File
@@ -29,6 +29,7 @@ add_library(soothe2_dsp SHARED
fn529fe0.cpp fn529fe0.cpp
rt_weights.cpp rt_weights.cpp
rt_mask_tables.cpp rt_mask_tables.cpp
log2_ln.cpp
) )
add_executable(soothe2_harness harness.cpp) add_executable(soothe2_harness harness.cpp)
+66 -2
View File
@@ -134,6 +134,32 @@ static void process_band_structural(
for (size_t k = 0; k < nbin; k++) if (lvl_in[k] > cap) lvl_in[k] = cap; for (size_t k = 0; k < nbin; k++) if (lvl_in[k] > cap) lvl_in[k] = cap;
} }
// Cascade sin-peak floor (529c60): the -20.72 dB floor mechanism.
// From assembly: sin_peak = sin(param * 30 - 90) * (ln10/20) * peak
// where ln10/20 = 0.115129 (constant at 0x1824c3cd4).
// This prevents over-reduction by clamping the level curve.
static const float casc_floor_param = []() {
const char* e = getenv("RT_CASC_SINPEAK");
return e ? (float)atof(e) : 0.0f;
}();
if (casc_floor_param != 0.0f) {
// Find peak of level curve
float peak_lvl = 0.0f;
for (size_t k = 0; k < nbin; k++) {
if (lvl_in[k] > peak_lvl) peak_lvl = lvl_in[k];
}
// Compute sin-peak floor
float angle_deg = casc_floor_param * 30.0f - 90.0f;
float sin_peak = std::sin(angle_deg * static_cast<float>(M_PI) / 180.0f)
* 0.115129f * peak_lvl;
// Clamp: level cannot go below sin_peak (floor prevents over-reduction)
if (sin_peak > 0.0f) {
for (size_t k = 0; k < nbin; k++) {
if (lvl_in[k] < sin_peak) lvl_in[k] = sin_peak;
}
}
}
// Save raw level BEFORE LUT transform (for RT_FIRPOWER) // Save raw level BEFORE LUT transform (for RT_FIRPOWER)
std::vector<float> raw_level(nbin); std::vector<float> raw_level(nbin);
for (size_t k = 0; k < nbin; k++) { for (size_t k = 0; k < nbin; k++) {
@@ -410,6 +436,8 @@ void FramedDetector::setParams(const std::vector<DetectorBand>& bands) {
size_t half = nfft_ / 2; size_t half = nfft_ / 2;
res_.clear(); res_.clear();
track_.clear(); track_.clear();
twin_resp_complex_.clear();
cascade_states_.clear();
// RT_DUMPRESPATH=<file> (NOTES 22t): static twin-response spectra per band, // RT_DUMPRESPATH=<file> (NOTES 22t): static twin-response spectra per band,
// binary {int32 band, int32 nbin, float res[nbin]} records (append). // binary {int32 band, int32 nbin, float res[nbin]} records (append).
@@ -434,6 +462,13 @@ void FramedDetector::setParams(const std::vector<DetectorBand>& bands) {
r[k] = std::sqrt(out[k].re * out[k].re + out[k].im * out[k].im); r[k] = std::sqrt(out[k].re * out[k].re + out[k].im * out[k].im);
r[k] = std::max(r[k], 1e-12f); r[k] = std::max(r[k], 1e-12f);
} }
// Store complex response for cascade 529c60
std::vector<std::complex<double>> complex_resp(half + 1);
for (size_t k = 0; k <= half; k++) {
complex_resp[k] = std::complex<double>(out[k].re, out[k].im);
}
twin_resp_complex_.push_back(std::move(complex_resp));
if (rp_dump) { if (rp_dump) {
int32_t bi = static_cast<int32_t>(res_.size()); int32_t bi = static_cast<int32_t>(res_.size());
int32_t nb = static_cast<int32_t>(r.size()); int32_t nb = static_cast<int32_t>(r.size());
@@ -445,6 +480,7 @@ void FramedDetector::setParams(const std::vector<DetectorBand>& bands) {
} }
if (rp_dump) fclose(rp_dump); if (rp_dump) fclose(rp_dump);
track_.assign(bands_.size(), std::vector<float>(half + 1, 1.0f)); track_.assign(bands_.size(), std::vector<float>(half + 1, 1.0f));
cascade_states_.assign(bands_.size(), fn529fe0::CascadeState());
} }
void FramedDetector::processFrame(const std::complex<double>* spectrum, float* mask) { void FramedDetector::processFrame(const std::complex<double>* spectrum, float* mask) {
@@ -481,6 +517,10 @@ void FramedDetector::processFrame(const std::complex<double>* spectrum, float* m
} }
} }
// Detector cascade 529c60: per-band pre-processor on complex twin-filtered
// spectrum. Computes magnitudes, Haar-smooths, applies sin-peak floor.
static const int casc_on = getenv("RT_CASC") ? atoi(getenv("RT_CASC")) : 0;
for (size_t k = 0; k <= half; k++) mask[k] = 1.0f; for (size_t k = 0; k <= half; k++) mask[k] = 1.0f;
if (is_internal_grid(nfft_, sample_rate_)) { if (is_internal_grid(nfft_, sample_rate_)) {
@@ -497,8 +537,32 @@ void FramedDetector::processFrame(const std::complex<double>* spectrum, float* m
sample_rate_, sf, fparams, sample_rate_, sf, fparams,
band_mask.data()); band_mask.data());
} else { } else {
process_band_structural(am_.data(), res_[b].data(), bands_[b], // Run cascade per-band on complex twin-filtered spectrum
band_mask.data(), nfft_, sample_rate_); if (casc_on && nfft_ == 4096) {
size_t nbin = half + 1;
std::vector<float> complex_input(2 * nbin);
std::vector<float> curve_output(nbin);
for (size_t k = 0; k <= half; k++) {
complex_input[2*k] = static_cast<float>(twin_resp_complex_[b][k].real());
complex_input[2*k+1] = static_cast<float>(twin_resp_complex_[b][k].imag());
}
fn529fe0::cascade_detect(
complex_input.data(),
curve_output.data(),
cascade_states_[b],
nbin,
2, // Haar iterations
0.0f, // sin_peak_param (0 = no floor)
48000.0f, // ctx[0x24] = sample rate
1, // ctx[0x1a0] = 1
4, // ctx[0x1ac] = 4 (quality default)
false // is_magnitude = false (input is complex)
);
}
process_band_structural(am_.data(), res_[b].data(), bands_[b],
band_mask.data(), nfft_, sample_rate_);
} }
for (size_t k = 0; k <= half; k++) { for (size_t k = 0; k <= half; k++) {
mask[k] = std::min(band_mask[k], mask[k]); mask[k] = std::min(band_mask[k], mask[k]);
+12
View File
@@ -2,6 +2,7 @@
#include <cstddef> #include <cstddef>
#include <complex> #include <complex>
#include <vector> #include <vector>
#include "fn529fe0.hpp"
struct DetectorBand { struct DetectorBand {
float fc; // band center freq (Hz) float fc; // band center freq (Hz)
@@ -71,6 +72,11 @@ public:
void processFrame(const std::complex<double>* spectrum, float* mask); void processFrame(const std::complex<double>* spectrum, float* mask);
// Cascade state access for per-band detector cascade
std::vector<fn529fe0::CascadeState>& cascadeStates() { return cascade_states_; }
const std::vector<std::vector<std::complex<double>>>& twinRespComplex() const { return twin_resp_complex_; }
std::vector<std::vector<std::complex<double>>>& twinRespComplex() { return twin_resp_complex_; }
private: private:
size_t nfft_; size_t nfft_;
float sample_rate_; float sample_rate_;
@@ -82,4 +88,10 @@ private:
std::vector<float> am_; // smoothed per-bin amplitude std::vector<float> am_; // smoothed per-bin amplitude
std::vector<float> f6f8_; // shared 0x5406f8 blend buffer (IIR1 out) std::vector<float> f6f8_; // shared 0x5406f8 blend buffer (IIR1 out)
std::vector<std::vector<float>> track_; // per band, per bin accumulator 0x5407c8 std::vector<std::vector<float>> track_; // per band, per bin accumulator 0x5407c8
// For cascade 529c60: per-band complex twin filter responses
std::vector<std::vector<std::complex<double>>> twin_resp_complex_;
// Per-band cascade states
std::vector<fn529fe0::CascadeState> cascade_states_;
}; };
+48
View File
@@ -0,0 +1,48 @@
#include "log2_ln.hpp"
#include <cstring>
#include <cmath>
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
+17
View File
@@ -0,0 +1,17 @@
#pragma once
#include <cstdint>
#include <cmath>
#include <cstring>
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
+163 -6
View File
@@ -1,9 +1,13 @@
#include "spectral.hpp" #include "spectral.hpp"
#include "fftconv.hpp" #include "fftconv.hpp"
#include "log2_ln.hpp"
#include "exp2_tables.hpp"
#include "exp2.hpp"
#include <cmath> #include <cmath>
#include <cstring> #include <cstring>
#include <vector> #include <vector>
#include <cstdlib> #include <cstdlib>
#include <cstdio>
SpectralProcessor::SpectralProcessor(size_t nfft, size_t hop, float sample_rate) SpectralProcessor::SpectralProcessor(size_t nfft, size_t hop, float sample_rate)
: nfft_(nfft), hop_(hop), frame_count_(0), output_pos_(0), : nfft_(nfft), hop_(hop), frame_count_(0), output_pos_(0),
@@ -36,6 +40,7 @@ SpectralProcessor::~SpectralProcessor() {
void SpectralProcessor::setDetectorParams(const std::vector<DetectorBand>& bands) { void SpectralProcessor::setDetectorParams(const std::vector<DetectorBand>& bands) {
detector_.setParams(bands); detector_.setParams(bands);
loadWinFreq();
} }
void SpectralProcessor::computeWindow() { void SpectralProcessor::computeWindow() {
@@ -88,6 +93,154 @@ void SpectralProcessor::istftFrame(std::complex<double>* in, float* out, float*
} }
} }
void SpectralProcessor::loadWinFreq() {
if (win_freq_loaded_) return;
win_freq_loaded_ = true;
// Try to load WIN_freq from live capture (handoff/rtwin_freq_44100.npy)
FILE* f = fopen("handoff/rtwin_freq_44100.npy", "rb");
if (!f) {
// Fallback: compute periodic Hann, second half (0.5→1.0 rising)
win_freq_.resize(nfft_ / 2 + 1);
for (size_t i = 0; i <= nfft_ / 2; i++) {
win_freq_[i] = static_cast<float>(0.5 * (1.0 - std::cos(2.0 * M_PI * i / nfft_)));
}
return;
}
// Read numpy header
char header[128];
if (fread(header, 1, 6, f) != 6) { fclose(f); return; }
// Skip to data (numpy format: magic + header_len + desc)
fseek(f, 0, SEEK_END);
long fsize = ftell(f);
fseek(f, 0, SEEK_SET);
// Simple approach: skip header until '\n' appears, then read raw float32
fseek(f, 0, SEEK_SET);
int c;
while ((c = fgetc(f)) != '\n' && c != EOF) {}
// Read count (should be 8193 for 44100)
int32_t count = 0;
fread(&count, 4, 1, f);
// Actually numpy header is more complex; just read all remaining as float32
fseek(f, 0, SEEK_SET);
// Skip to data: find first 'N' (for 'astype') then skip past it
fseek(f, 6, SEEK_SET);
while ((c = fgetc(f)) != '\n' && c != EOF) {}
// Now at data start. Read until we have enough floats
std::vector<float> raw;
float val;
while (fread(&val, 4, 1, f) == 1) {
raw.push_back(val);
}
fclose(f);
if (raw.size() > 0) {
win_freq_ = raw;
} else {
// Fallback
win_freq_.resize(nfft_ / 2 + 1);
for (size_t i = 0; i <= nfft_ / 2; i++) {
win_freq_[i] = static_cast<float>(0.5 * (1.0 - std::cos(2.0 * M_PI * i / nfft_)));
}
}
}
void SpectralProcessor::buildFirFromMask(const float* mask, std::complex<double>* fir, size_t nbin) {
// 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-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++) {
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 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());
// 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;
}
for (size_t i = half; i < nfft; i++) {
H[i] = std::complex<double>(0.0, 0.0);
}
// FFT back to freq domain
fft::execute(&plan_, H.data());
// 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();
}
for (size_t i = 0; i < nfft; i++) {
fir[i] = H[i] * scale;
}
fir[0] = std::complex<double>(1.0, 0.0);
if (half >= 1) {
fir[1] = std::complex<double>(0.0, 0.0);
}
}
void SpectralProcessor::processBlock(float* in, float* out, size_t num_samples, size_t num_channels) { void SpectralProcessor::processBlock(float* in, float* out, size_t num_samples, size_t num_channels) {
memset(out, 0, num_samples * sizeof(float)); memset(out, 0, num_samples * sizeof(float));
if (num_samples == 0 || num_samples < nfft_) { if (num_samples == 0 || num_samples < nfft_) {
@@ -118,16 +271,20 @@ void SpectralProcessor::processBlock(float* in, float* out, size_t num_samples,
buf_[i] *= a; buf_[i] *= a;
} }
} else if (firconv) { } else if (firconv) {
// RT_FIRCONV=1: Build FIR from mask and apply via complex multiply. // RT_FIRCONV=2: Full FIR construction pipeline (52b550-52b8bb).
// The mask is real-valued (per-bin gain). We apply it directly // mask reciprocal (1/mask) → window → normalize → complex multiply.
// to the audio spectrum via complex multiply (th_b3c0 equivalent). // This replicates the plugin's FFT-conv FIR design path.
// No upper-half zeroing — preserve Hermitian symmetry. buildFirFromMask(mask_.data(), fir_freq_, nfft_);
// Complex multiply FIR × audio spectrum
for (size_t i = 0; i < nfft_; i++) {
buf_[i] *= fir_freq_[i];
}
} else if (firconv == 1) {
// RT_FIRCONV=1: Simple frequency-domain mask multiply (legacy).
for (size_t i = 0; i < nfft_; i++) { for (size_t i = 0; i < nfft_; i++) {
fir_freq_[i] = std::complex<double>( fir_freq_[i] = std::complex<double>(
static_cast<double>(mask_[i % (nfft_/2+1)]), 0.0); static_cast<double>(mask_[i % (nfft_/2+1)]), 0.0);
} }
// Complex multiply FIR × audio spectrum.
for (size_t i = 0; i < nfft_; i++) { for (size_t i = 0; i < nfft_; i++) {
buf_[i] *= fir_freq_[i]; buf_[i] *= fir_freq_[i];
} }
+10
View File
@@ -37,4 +37,14 @@ private:
void computeWindow(); void computeWindow();
void stftFrame(const float* in, std::complex<double>* out); void stftFrame(const float* in, std::complex<double>* out);
void istftFrame(std::complex<double>* in, float* out, float* overlap); void istftFrame(std::complex<double>* in, float* out, float* overlap);
// FIR construction from detector mask (52b550-52b8bb pipeline):
// mask → log → sign-invert → EXP → twiddle ops → window → normalize
// Produces frequency-domain FIR kernel for complex multiply application.
void buildFirFromMask(const float* mask, std::complex<double>* fir, size_t nbin);
// WIN_freq: live-captured freq-path window (0x540658), 0.5→1.0
std::vector<float> win_freq_;
bool win_freq_loaded_ = false;
void loadWinFreq();
}; };