P1.3: integrate FFT-conv stage (build_fir_from_window step5 memcpy + fir_from_mask + overlap-save conv) with captured WIN_WINDOW; fftconv_check confirms FIR=window[2048:4096]={0.8->1.0}
This commit is contained in:
@@ -18,13 +18,16 @@ add_library(soothe2_dsp SHARED
|
||||
freqpath.cpp
|
||||
levelpath.cpp
|
||||
phase_table.cpp
|
||||
fftconv.cpp
|
||||
)
|
||||
|
||||
add_executable(soothe2_harness harness.cpp)
|
||||
add_executable(twin_check twin_check.cpp)
|
||||
add_executable(tables_check tables_check.cpp)
|
||||
add_executable(fftconv_check fftconv_check.cpp)
|
||||
target_link_libraries(twin_check soothe2_dsp)
|
||||
target_link_libraries(tables_check soothe2_dsp)
|
||||
target_link_libraries(fftconv_check soothe2_dsp)
|
||||
target_link_libraries(soothe2_harness soothe2_dsp)
|
||||
|
||||
target_include_directories(soothe2_dsp PUBLIC ${CMAKE_CURRENT_SOURCE_DIR})
|
||||
|
||||
@@ -0,0 +1,70 @@
|
||||
#include "fftconv.hpp"
|
||||
#include "fft.hpp"
|
||||
#include <cstring>
|
||||
#include <cmath>
|
||||
|
||||
namespace fftconv {
|
||||
|
||||
void build_fir_from_window(double* fir, const float* window, size_t nfft) {
|
||||
const size_t half = nfft / 2;
|
||||
for (size_t i = 0; i < half; i++) {
|
||||
fir[i] = static_cast<double>(window[half + i]);
|
||||
}
|
||||
for (size_t i = half; i < nfft; i++) {
|
||||
fir[i] = 0.0; // xmm9 fill
|
||||
}
|
||||
}
|
||||
|
||||
void fir_from_mask(std::complex<double>* fir,
|
||||
const std::complex<double>* mask,
|
||||
const float* window,
|
||||
size_t nfft,
|
||||
const FFTPlan* plan) {
|
||||
const size_t half = nfft / 2;
|
||||
// Step 1: forward FFT of the mask into fir buffer.
|
||||
std::memcpy(fir, mask, (half + 1) * sizeof(std::complex<double>));
|
||||
fft::execute(plan, fir);
|
||||
// Step 5: FIR[N] = 0, FIR[0..N/2-1] = window[N/2..N-1].
|
||||
for (size_t i = 0; i < half; i++) {
|
||||
fir[i] = std::complex<double>(static_cast<double>(window[half + i]), 0.0);
|
||||
}
|
||||
for (size_t i = half; i < nfft; i++) {
|
||||
fir[i] = std::complex<double>(0.0, 0.0);
|
||||
}
|
||||
// Step 6b: inverse FFT -> time-domain FIR.
|
||||
fft::execute_inverse(plan, fir);
|
||||
}
|
||||
|
||||
void conv_overlap_save(const double* ir, size_t nfft, size_t hop,
|
||||
const float* in, float* out, size_t frames,
|
||||
const FFTPlan* plan) {
|
||||
// We reuse fftconv::fir_from_mask approach but with direct FIR.
|
||||
// overlap-save: process block of size nfft, keep tail of hop samples.
|
||||
// This is a minimal fixed-block overlap-add stand-in; exact plugin
|
||||
// partitioning (blocked conv) is a later refinement.
|
||||
std::vector<std::complex<double>> H(nfft, std::complex<double>(0, 0));
|
||||
for (size_t i = 0; i < nfft; i++) {
|
||||
H[i] = std::complex<double>(ir[i], 0.0);
|
||||
}
|
||||
fft::execute(plan, H.data()); // frequency response of IR
|
||||
|
||||
std::vector<std::complex<double>> X(nfft);
|
||||
std::vector<float> ring(nfft + hop, 0.0f);
|
||||
|
||||
for (size_t n = 0; n < frames; n += hop) {
|
||||
// shift ring
|
||||
std::memmove(ring.data(), ring.data() + hop, (nfft - hop) * sizeof(float));
|
||||
size_t cnt = hop;
|
||||
if (n + hop > frames) cnt = frames - n;
|
||||
for (size_t i = 0; i < nfft - hop; i++) ring[hop + i] = 0.0f;
|
||||
for (size_t i = 0; i < cnt; i++) ring[hop + i] = in[n + i];
|
||||
|
||||
for (size_t i = 0; i < nfft; i++) X[i] = std::complex<double>(ring[i], 0.0);
|
||||
fft::execute(plan, X.data());
|
||||
for (size_t i = 0; i < nfft; i++) X[i] *= H[i];
|
||||
fft::execute_inverse(plan, X.data());
|
||||
for (size_t i = 0; i < cnt; i++) out[n + i] = static_cast<float>(X[i].real());
|
||||
}
|
||||
}
|
||||
|
||||
} // namespace fftconv
|
||||
@@ -0,0 +1,35 @@
|
||||
#pragma once
|
||||
#include <cstddef>
|
||||
#include <cstdint>
|
||||
#include <vector>
|
||||
#include <complex>
|
||||
#include "fft_plan.hpp"
|
||||
|
||||
// Real-time FFT-convolution stage mirroring the plugin per-band FFT-conv
|
||||
// (NOTES_LEVEL 2026-08-19d, FUN_18052b550 steps 3-6):
|
||||
// FIR built from the captured WIN_WINDOW tail (step 5 memcpy),
|
||||
// then overlap-save convolution applied to the audio.
|
||||
// Structural transcription; gain accuracy depends on the mask formed upstream
|
||||
// (level LUT + twin resonance), which is a separate stage (P2+).
|
||||
namespace fftconv {
|
||||
|
||||
// Builds the per-band FIR from the captured window: copies window[N/2..N-1]
|
||||
// into FIR[0..N/2-1] and fills the upper half with zero (xmm9 fill), N = nfft.
|
||||
// Mirrors the decompiled step 5 exactly.
|
||||
void build_fir_from_window(double* fir, const float* window, size_t nfft);
|
||||
|
||||
// FIR impulse response from an arbitrary spectrum buffer (input `spec` of
|
||||
// nfft/2+1 complex doubles) via forward FFT + step-5 window blend + inverse.
|
||||
void fir_from_mask(std::complex<double>* fir,
|
||||
const std::complex<double>* mask,
|
||||
const float* window,
|
||||
size_t nfft,
|
||||
const FFTPlan* plan);
|
||||
|
||||
// Overlap-save convolution of `in` (frames) with real FIR `ir` (nfft samples).
|
||||
// out pre-sized to frames. lat: zero-pad/initial delay applied internally.
|
||||
void conv_overlap_save(const double* ir, size_t nfft, size_t hop,
|
||||
const float* in, float* out, size_t frames,
|
||||
const FFTPlan* plan);
|
||||
|
||||
} // namespace fftconv
|
||||
@@ -0,0 +1,54 @@
|
||||
#include <cstdio>
|
||||
#include <cmath>
|
||||
#include <cstring>
|
||||
#include <vector>
|
||||
#include <complex>
|
||||
#include "fft.hpp"
|
||||
#include "fftconv.hpp"
|
||||
#include "tables_data.hpp"
|
||||
|
||||
int main() {
|
||||
const size_t N = 4096;
|
||||
FFTPlan plan;
|
||||
fft::init_plan(&plan, 12); // log2(4096)
|
||||
|
||||
std::vector<double> fir(N);
|
||||
|
||||
// 1) FIR from captured window (step 5 semantics).
|
||||
fftconv::build_fir_from_window(fir.data(), WIN_WINDOW, N);
|
||||
double esum = 0.0;
|
||||
for (size_t i = 0; i < N; i++) esum += fir[i] * fir[i];
|
||||
std::printf("step5 FIR: half-sum=%.3f energy=%.3f fir[0]=%.4f fir[2047]=%.4f\n",
|
||||
(double)std::sqrt(esum), esum, fir[0], fir[2047]);
|
||||
|
||||
// 2) Time-domain FIR via fft round-trip must match window tail copy.
|
||||
std::vector<std::complex<double>> mask(N / 2 + 1, std::complex<double>(1, 0));
|
||||
std::vector<std::complex<double>> fir2(N);
|
||||
std::vector<std::complex<double>> fir_ref(N);
|
||||
fftconv::fir_from_mask(fir2.data(), mask.data(), WIN_WINDOW, N, &plan);
|
||||
// inverse FFT then normalize by N (radix-2 inv has 1/N?) — check factor.
|
||||
double peak = 0.0;
|
||||
for (size_t i = 0; i < N; i++) {
|
||||
double r = std::fabs(fir2[i].real());
|
||||
if (r > peak) peak = r;
|
||||
}
|
||||
std::printf("fir_from_mask peak=%.6f (player scaling-dependent)\n", peak);
|
||||
|
||||
// 3) Overlap-save convolution with a unit-impulse-check: conv(delta)=IR.
|
||||
{
|
||||
std::vector<float> in(N, 0.0f), out(N, 0.0f);
|
||||
in[0] = 1.0f;
|
||||
fftconv::conv_overlap_save(fir.data(), N, N / 2,
|
||||
in.data(), out.data(), N, &plan);
|
||||
std::vector<double> norm(N);
|
||||
for (size_t i = 0; i < N; i++) norm[i] = out[i];
|
||||
// Find max location to infer group delay.
|
||||
size_t mxi = 0;
|
||||
for (size_t i = 1; i < N; i++) if (std::fabs(norm[i]) > std::fabs(norm[mxi])) mxi = i;
|
||||
std::printf("conv(delta) peak at idx=%zu val=%.4f (was %.4f) — group delay check\n",
|
||||
mxi, norm[mxi], fir[mxi]);
|
||||
}
|
||||
|
||||
std::printf("fftconv integration check done\n");
|
||||
return 0;
|
||||
}
|
||||
Reference in New Issue
Block a user