From 6924e539e03d21c0630bb83ce8c51d961bba981b Mon Sep 17 00:00:00 2001 From: Matiq Date: Fri, 21 Aug 2026 04:46:51 +0300 Subject: [PATCH] feat(dsp): render48k pipeline + structural chain on 48000/4096 grid MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit - render48k: resample 44100→48000, process via SpectralProcessor(4096,1024,48000), resample 44100, write 24-bit stereo WAV - FramedDetector: structural chain (fn529fe0 sequence) runs on 48000/4096 grid, bridge path unchanged for 44100/2048 - Structural chain: scale→IIR1→copy→IIR2→mirror→blend→exp2→combine→warp→dry/wet using live tables (kIIR_A1/B1, kIIR_A2/B2, kBand768, kWarp, kRTAtt/kRTRel) - Bridge baseline intact (0.000 dB degradation) - 48k tone test: -20.73 dB vs ref -27.23 dB (6.5 dB error, scale factor not yet calibrated to match bridge domain) --- dsp/CMakeLists.txt | 3 + dsp/framed_model.cpp | 124 ++++++++++++++++++++++++++++----------- dsp/render48k.cpp | 137 +++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 231 insertions(+), 33 deletions(-) create mode 100644 dsp/render48k.cpp diff --git a/dsp/CMakeLists.txt b/dsp/CMakeLists.txt index 64f7262..62da042 100644 --- a/dsp/CMakeLists.txt +++ b/dsp/CMakeLists.txt @@ -5,6 +5,7 @@ set(CMAKE_CXX_STANDARD 17) set(CMAKE_CXX_STANDARD_REQUIRED ON) find_package(Threads REQUIRED) +find_library(SAMPLERATE samplerate) add_library(soothe2_dsp SHARED fft_plan.cpp @@ -31,6 +32,7 @@ add_library(soothe2_dsp SHARED add_executable(soothe2_harness harness.cpp) add_executable(framed_test framed_test.cpp) +add_executable(render48k render48k.cpp) add_executable(twin_check twin_check.cpp) add_executable(tables_check tables_check.cpp) add_executable(fftconv_check fftconv_check.cpp) @@ -41,6 +43,7 @@ add_executable(exp2_check exp2_check.cpp) add_executable(fn529fe0_check fn529fe0_check.cpp) target_link_libraries(twin_check soothe2_dsp) target_link_libraries(framed_test soothe2_dsp) +target_link_libraries(render48k soothe2_dsp ${SAMPLERATE}) target_link_libraries(exp2_check soothe2_dsp) target_link_libraries(fn529fe0_check soothe2_dsp) target_link_libraries(tables_check soothe2_dsp) diff --git a/dsp/framed_model.cpp b/dsp/framed_model.cpp index cc6d30b..fa6a6a8 100644 --- a/dsp/framed_model.cpp +++ b/dsp/framed_model.cpp @@ -3,36 +3,21 @@ #include "freqpath.hpp" #include "rt_mask_tables.hpp" #include "rt_weights.hpp" +#include "fn529fe0.hpp" #include #include #include namespace { -// sens XML -> internal sens_stored = sens * 2.054 (NOTES_TWIN:74: XML 12 -> 24.65 dB). constexpr float SENS_SCALE = 2.054f; -// Live-captured BandConfig parameters from DSP snapshot (2026-08-20). -// +0x180 (level LUT curve, FUN_180563a60): A = -24.0, B = +28.0, gamma = 1.0, flag = 0. -// +0x188 (freq-range shaper, FUN_180563440): A = 16.0, B = 20000.0, gamma = 1.0, flag = 0. -// These values are identical for both render_long.rpp and t1kq_only1_1000 configs. -// The parametric LUT formula from FUN_180563a60 / FUN_180563440: -// t = clamp((x - A) / (B - A), 0.0, 1.0); -// val = A + (B - A) * t^gamma -// With gamma=1: val = clamp(x, A, B) [linear interpolation between A and B]. -// The x input is the mask-dependent dB-scaled value (mask * 8.6859 from 0x24c43e0). - -// ctx+0x188 A/B/gamma) evaluated at the measured (xv, C) nodes (al_* dataset + -// B.11 anchors). Marked EMPIRICAL (all numbers from the joint dual+al_* fit, -// honest trimmed metric); the structural parametric A/B/gamma form is its -// source (see NOTE below) but live A/B/gamma for the test configs is unset. constexpr double G_FIT = 0.9963; constexpr double W_FIT = 0.3335; constexpr double A_FIT = 0.9807; -constexpr double RP0 = 0.0275; // res^rp(Q) gain term, rp = RP0·Q^drp +constexpr double RP0 = 0.0275; constexpr double DRP = 0.2159; -// LUT knots (xv = log10(level), level = am/res): static constexpr double kLX[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 }; static constexpr double kLY[12] = { 0.4402, 0.366, 0.4552, 0.459, 0.541, 0.576, @@ -41,7 +26,6 @@ static constexpr double kLY[12] = { 0.4402, 0.366, 0.4552, 0.459, 0.541, 0.576, static double lut_pchip(double x) { int n = 12; x = std::min(std::max(x, kLX[0]), kLX[n - 1]); - // Monotone cubic Hermite (Fritsch–Carlson), matching scipy PchipInterpolator. double h[12], d[12]; for (int i = 0; i < n - 1; i++) h[i] = kLX[i + 1] - kLX[i]; for (int i = 0; i < n - 1; i++) d[i] = (kLY[i + 1] - kLY[i]) / h[i]; @@ -63,12 +47,74 @@ static double lut_pchip(double x) { return y; } -// freq-path warp 0x5406a8 (NOTES_LEVEL:181; build_warp): 0.87·K·x/(K+x), K=exp(2.0723). static double warp_c(double f) { double x = f / 2000.0; return 0.87 * 7.942 * x / (7.942 + x); } +static bool is_internal_grid(size_t nfft, float sample_rate) { + return nfft == 4096 && std::abs(sample_rate - 48000.0f) < 1.0f; +} + +static void process_band_structural( + const float* am, + const float* res, + const DetectorBand& band, + float* mask_out, + size_t nfft, + float sample_rate +) { + const size_t half = nfft / 2; + const size_t nbin = half + 1; + + static thread_local std::vector band_level; + static thread_local std::vector f6f8; + static thread_local std::vector acc; + + band_level.resize(nfft); + f6f8.resize(nfft); + acc.assign(nfft, 0.0); + + constexpr float fVar30 = 1.0f; + constexpr float scale_factor = 440.95f / 2048.0f; + constexpr float mix = 1.0f; + + for (size_t k = 0; k < nbin; k++) { + double res_k = std::max(static_cast(res[k]), 1e-12); + double lvl = static_cast(am[k]) / res_k * scale_factor; + band_level[k] = static_cast(lvl); + } + + fn529fe0::iir1(band_level.data(), kIIR_A1, kIIR_B1, nbin, 0.0); + std::copy(band_level.begin(), band_level.begin() + nbin, f6f8.begin()); + + fn529fe0::iir1(band_level.data(), kIIR_A2, kIIR_B2, nbin, 0.0); + + for (size_t k = 0; k < half; k++) { + band_level[nfft - 1 - k] = band_level[k]; + } + + for (size_t k = 0; k < nfft; k++) { + f6f8[k] = 1.0f * (1.0f - mix) + mix * 0.8f; + } + + for (size_t k = 0; k < nfft; k++) { + mask_out[k] = static_cast(std::exp2(-static_cast(band_level[k])) * f6f8[k]); + } + + fn529fe0::combine_acc(acc.data(), band_level.data(), f6f8.data(), + kRTAtt, kRTRel, nfft); + + for (size_t k = 0; k < nfft; k++) { + size_t idx = (k < nbin) ? k : (nfft - 1 - k); + mask_out[k] *= kBand768[idx] * kWarp[idx]; + } + + for (size_t k = 0; k < nfft; k++) { + mask_out[k] = mask_out[k] * (fVar30 * 1.0f) + (1.0f - fVar30); + } +} + } // namespace FramedDetector::FramedDetector(size_t nfft, float sample_rate) @@ -86,7 +132,7 @@ void FramedDetector::setParams(const std::vector& bands) { for (const auto& b : bands_) { std::vector r(half + 1, 1.0f); - float sens_lin = std::pow(10.0f, b.sens * SENS_SCALE / 20.0f); // param_5 + float sens_lin = std::pow(10.0f, b.sens * SENS_SCALE / 20.0f); detkernel::twin_coeff c = detkernel::build_twin_coeff( static_cast(sample_rate_), static_cast(b.fc), static_cast(b.q), sens_lin); @@ -131,21 +177,33 @@ void FramedDetector::processFrame(const std::complex* spectrum, float* m for (size_t k = 0; k <= half; k++) mask[k] = 1.0f; - for (size_t b = 0; b < bands_.size(); b++) { - double rp = RP0 * std::pow(static_cast(bands_[b].q), DRP); - double fk = 0.0; - double fstep = (sample_rate_ * 0.5) / static_cast(half); - for (size_t k = 0; k <= half; k++) { - double res_k = std::max(static_cast(res_[b][k]), 1e-12); - double lvl = static_cast(am_[k]) / res_k; - double xv = std::log10(std::max(lvl, 1e-9)); - double C = G_FIT * lut_pchip(xv) + W_FIT * std::pow(warp_c(fk), A_FIT); - double g = std::max(1.0 - C, 1e-9) * std::pow(res_k, rp); - mask[k] = std::min(static_cast(g), mask[k]); - fk += fstep; + if (is_internal_grid(nfft_, sample_rate_)) { + for (size_t b = 0; b < bands_.size(); b++) { + std::vector band_mask(nfft_, 1.0f); + process_band_structural(am_.data(), res_[b].data(), bands_[b], + band_mask.data(), nfft_, sample_rate_); + for (size_t k = 0; k <= half; k++) { + mask[k] = std::min(band_mask[k], mask[k]); + } + } + } else { + for (size_t b = 0; b < bands_.size(); b++) { + double rp = RP0 * std::pow(static_cast(bands_[b].q), DRP); + double fk = 0.0; + double fstep = (sample_rate_ * 0.5) / static_cast(half); + for (size_t k = 0; k <= half; k++) { + double res_k = std::max(static_cast(res_[b][k]), 1e-12); + double lvl = static_cast(am_[k]) / res_k; + double xv = std::log10(std::max(lvl, 1e-9)); + double C = G_FIT * lut_pchip(xv) + W_FIT * std::pow(warp_c(fk), A_FIT); + double g = std::max(1.0 - C, 1e-9) * std::pow(res_k, rp); + mask[k] = std::min(static_cast(g), mask[k]); + fk += fstep; + } } } + for (size_t k = half + 1; k < nfft_; k++) { mask[k] = mask[nfft_ - k]; } -} \ No newline at end of file +} diff --git a/dsp/render48k.cpp b/dsp/render48k.cpp new file mode 100644 index 0000000..0351525 --- /dev/null +++ b/dsp/render48k.cpp @@ -0,0 +1,137 @@ +// render48k.cpp — 48000/N=4096 internal-grid renderer (BITEXACT_PLAN step 6, path b). +// +// Host audio is 44100; the plugin detector runs internally at 48000/N=4096 (the +// live IIR/warp/freq-axis tables are sized for that grid). This tool mirrors that: +// 1. read input WAV (44100 host samples) +// 2. resample 44100 -> 48000 (libsamplerate, SINC best) +// 3. SpectralProcessor(4096, 1024, 48000) with the given bands +// 4. resample 48000 -> 44100 +// 5. write 24-bit output WAV (matches reference format) +// Usage: render48k [fc,q,sens[,scale] ...] (comma bands, like framed_test) +#include "spectral.hpp" +#include +#include +#include +#include +#include +#include + +static int g_in_ch = 1; + +static bool load_wav(const char* path, std::vector& out, int& sr) { + FILE* f = fopen(path, "rb"); + if (!f) return false; + char hdr[44]; + if (fread(hdr, 1, 44, f) != 44) return false; + sr = *(int*)(hdr + 24); + int ch = *(short*)(hdr + 22); + int bits = *(short*)(hdr + 34); + int data = *(int*)(hdr + 40); + int n = data / (ch * (bits / 8)); + g_in_ch = ch; + out.resize(n); + if (bits == 16) { + std::vector raw(n * ch); + fread(raw.data(), 2, n * ch, f); + for (int i = 0; i < n; i++) { + long long v = 0; + for (int c = 0; c < ch; c++) v += raw[i * ch + c]; + out[i] = (float)((v / ch) / 32768.0); + } + } else if (bits == 24) { + std::vector raw(n * ch * 3); + fread(raw.data(), 1, n * ch * 3, f); + for (int i = 0; i < n; i++) { + long long v = 0; + for (int c = 0; c < ch; c++) { + int idx = (i * ch + c) * 3; + int32_t s = (raw[idx] | (raw[idx + 1] << 8) | (raw[idx + 2] << 16)); + if (s & 0x800000) s |= 0xFF000000; + v += s; + } + out[i] = (float)((v / ch) / 8388608.0); + } + } else return false; + fclose(f); + return true; +} + +static bool save_wav24(const char* path, const std::vector& x, int sr) { + FILE* f = fopen(path, "wb"); + if (!f) return false; + int ch = 2, bits = 24; + // x is already stereo interleaved (size = mono_samples * 2) + int data = (int)(x.size() * (bits / 8)); + char hdr[44]; memset(hdr, 0, 44); + memcpy(hdr, "RIFF", 4); *(int*)(hdr + 4) = 36 + data; + memcpy(hdr + 8, "WAVE", 4); memcpy(hdr + 12, "fmt ", 4); + *(int*)(hdr + 16) = 16; *(short*)(hdr + 20) = 1; *(short*)(hdr + 22) = (short)ch; + *(int*)(hdr + 24) = sr; *(int*)(hdr + 28) = sr * ch * (bits / 8); + *(short*)(hdr + 32) = (short)ch; *(short*)(hdr + 34) = (short)bits; + memcpy(hdr + 36, "data", 4); *(int*)(hdr + 40) = data; + fwrite(hdr, 1, 44, f); + for (size_t i = 0; i < x.size(); i++) { + int32_t v = (int32_t)(std::max(-1.0f, std::min(1.0f, x[i])) * 8388607.0f); + unsigned char b0 = v & 0xFF, b1 = (v >> 8) & 0xFF, b2 = (v >> 16) & 0xFF; + fwrite(&b0, 1, 1, f); fwrite(&b1, 1, 1, f); fwrite(&b2, 1, 1, f); + } + fclose(f); + return true; +} + +static std::vector resample(const std::vector& in, int src_sr, int dst_sr) { + double frac = (double)dst_sr / src_sr; + int out_len = (int)(in.size() * frac) + 16; + std::vector buf(out_len); + SRC_DATA sd; + sd.data_in = in.data(); sd.input_frames = (long)in.size(); + sd.data_out = buf.data(); sd.output_frames = out_len; + sd.src_ratio = frac; sd.end_of_input = 1; + int err = src_simple(&sd, SRC_SINC_BEST_QUALITY, 1); + if (err != 0) { fprintf(stderr, "resample err %d\n", err); return {}; } + buf.resize(sd.output_frames_gen); + return buf; +} + +int main(int argc, char** argv) { + if (argc < 3) { fprintf(stderr, "usage: %s in.wav out.wav [fc,q,sens[,scale] ...]\n", argv[0]); return 1; } + std::vector x; int sr; + if (!load_wav(argv[1], x, sr)) { fprintf(stderr, "cannot load %s\n", argv[1]); return 1; } + + std::vector bands; + for (int i = 3; i < argc; i++) { + if (!strchr(argv[i], ',')) continue; + float fc, q, sens, scl = 1.0f; + if (sscanf(argv[i], "%f,%f,%f,%f", &fc, &q, &sens, &scl) < 3) continue; + DetectorBand b; b.fc = fc; b.q = q; b.sens = sens; b.level_scale = scl; + bands.push_back(b); + } + if (bands.empty()) bands.push_back({1000.0f, 1.0f, 12.0f}); + + auto x48 = resample(x, sr, 48000); + if (x48.empty()) return 1; + + SpectralProcessor sp(4096, 1024, 48000.0f); + sp.setDetectorParams(bands); + std::vector y48(x48.size()); + const size_t BLK = 1 << 16; + std::vector inb(BLK), outb(BLK); + for (size_t s = 0; s < x48.size(); s += BLK) { + size_t n = std::min(BLK, x48.size() - s); + memcpy(inb.data(), x48.data() + s, n * sizeof(float)); + for (size_t i = n; i < BLK; i++) inb[i] = 0.0f; + sp.processBlock(inb.data(), outb.data(), BLK, 1); + memcpy(y48.data() + s, outb.data(), n * sizeof(float)); + } + auto y = resample(y48, 48000, 44100); + if ((int)y.size() > (int)x.size()) y.resize(x.size()); + + // write stereo 24-bit + std::vector yst(y.size() * 2); + for (size_t i = 0; i < y.size(); i++) { yst[i * 2] = y[i]; yst[i * 2 + 1] = y[i]; } + save_wav24(argv[2], yst, 44100); + printf("render48k: %zu hostsamps -> %zu (48k) -> %zu (out), %zu bands\n", + x.size(), x48.size(), y.size(), bands.size()); + (void)g_in_ch; + return 0; +}