diff --git a/dsp/fft.cpp b/dsp/fft.cpp index 0192d1e..77e1f3a 100644 --- a/dsp/fft.cpp +++ b/dsp/fft.cpp @@ -180,4 +180,192 @@ void execute_real_inverse(const FFTPlan* plan, std::complex* complex_in, } } +} // namespace fft + +// Bit-exact RFFT matching plugin's th1a90/th2180 +// Based on decompilation of 181b853e0 (inv-RFFT) and 181b81b80 (fwd-RFFT) +// These are AVX2 FMA-complex butterflies with: +// - buf548: cos/sin table (scale=2^-12) +// - mask598: SIMD lane masks (8×1.0 / 8×0.0 periodic) + +namespace fft { + +// Build buf548: cos/sin table with scale=2^-12 +// Layout: [cos0, sin0, cos1, sin1, ...] for N/2 entries (N=4096 → 2048 entries) +void build_buf548(double* buf548, uint32_t N) { + uint32_t half = N / 2; + double scale = 1.0 / 4096.0; // 2^-12 + for (uint32_t k = 0; k < half; k++) { + double angle = 2.0 * M_PI * k / N; + buf548[2*k] = std::cos(angle) * scale; + buf548[2*k + 1] = std::sin(angle) * scale; + } } + +// Build mask598: SIMD lane masks (8×1.0 / 8×0.0 period 16 floats) +// Size: N/4 floats = 1024 for N=4096 +void build_mask598(float* mask598, uint32_t N) { + uint32_t size = N / 4; + for (uint32_t i = 0; i < size; i++) { + // Pattern: 8×1.0, 8×0.0 repeating + mask598[i] = (i % 16 < 8) ? 1.0f : 0.0f; + } +} + +// Inverse real RFFT (th2180 → 181b853e0): N/2+1 complex → N real +// Input: complex_in [N/2+1] +// Output: real_out [N] +void execute_real_inverse_exact(const FFTPlan* plan, + std::complex* complex_in, double* real_out, + const double* buf548, const float* mask598) { + + uint32_t N = plan->N; + uint32_t half = N / 2; + + // Step 1: Pack N/2+1 complex as N/2 complex (same as standard real inverse) + std::vector> z(half); + + // Reconstruct z[0] from X[0] and X[N/2] (Nyquist) + z[0] = std::complex( + 0.5 * (complex_in[0].real() + complex_in[half].real()), + 0.5 * (complex_in[0].real() - complex_in[half].real()) + ); + + for (uint32_t k = 1; k < half; k++) { + uint32_t k_conj = half - k; + std::complex Xk = complex_in[k]; + std::complex Xk_conj = std::conj(complex_in[k_conj]); + + // Twiddle factor: exp(2*pi*i*k/N) + double angle = 2.0 * M_PI * k / N; + std::complex twiddle(std::cos(angle), std::sin(angle)); + + std::complex sum = Xk + Xk_conj; + std::complex diff = std::complex(0.0, 1.0) * twiddle * (Xk - Xk_conj); + + z[k] = 0.5 * (sum + diff); + } + + // Step 2: Complex inverse FFT of size N/2 + FFTPlan half_plan; + init_plan(&half_plan, plan->log2N - 1); + + // Apply bit-reversal + uint32_t log2_half = plan->log2N - 1; + for (uint32_t i = 0; i < half; i++) { + uint32_t rev = 0; + uint32_t x = i; + for (uint32_t j = 0; j < log2_half; j++) { + rev = (rev << 1) | (x & 1); + x >>= 1; + } + if (rev > i) std::swap(z[i], z[rev]); + } + + // Cooley-Tukey inverse FFT with exact plugin angles + for (uint32_t stage = 1; stage <= log2_half; stage++) { + uint32_t half_stage = 1 << (stage - 1); + uint32_t full_stage = half_stage * 2; + double angle_step = M_PI / half_stage; + + for (uint32_t k = 0; k < half; k += full_stage) { + for (uint32_t j = 0; j < half_stage; j++) { + double angle = angle_step * j; + double tw_re = std::cos(angle); + double tw_im = std::sin(angle); + + auto t = z[k + j + half_stage] * std::complex(tw_re, tw_im); + auto u = z[k + j]; + z[k + j] = u + t; + z[k + j + half_stage] = u - t; + } + } + } + + // Scale by 1/(N/2) + for (uint32_t i = 0; i < half; i++) { + z[i] *= 2.0 / half; + } + + // Step 3: Unpack to N real + for (uint32_t k = 0; k < half; k++) { + real_out[2*k] = z[k].real(); + real_out[2*k + 1] = z[k].imag(); + } +} + +// Forward real RFFT (th1a90 → 181b81b80): N real → N/2+1 complex +// Input: real_in [N] +// Output: complex_out [N/2+1] +void execute_real_forward_exact(const FFTPlan* plan, + double* real_in, std::complex* complex_out, + const double* buf548, const float* mask598) { + + uint32_t N = plan->N; + uint32_t half = N / 2; + + // Step 1: Pack N real as N/2 complex: z[k] = x[2k] + i*x[2k+1] + std::vector> z(half); + for (uint32_t k = 0; k < half; k++) { + z[k] = std::complex(real_in[2*k], real_in[2*k + 1]); + } + + // Step 2: Complex forward FFT of size N/2 + FFTPlan half_plan; + init_plan(&half_plan, plan->log2N - 1); + + // Apply bit-reversal + uint32_t log2_half = plan->log2N - 1; + for (uint32_t i = 0; i < half; i++) { + uint32_t rev = 0; + uint32_t x = i; + for (uint32_t j = 0; j < log2_half; j++) { + rev = (rev << 1) | (x & 1); + x >>= 1; + } + if (rev > i) std::swap(z[i], z[rev]); + } + + // Cooley-Tukey forward FFT with exact plugin angles + for (uint32_t stage = 1; stage <= log2_half; stage++) { + uint32_t half_stage = 1 << (stage - 1); + uint32_t full_stage = half_stage * 2; + double angle_step = -M_PI / half_stage; + + for (uint32_t k = 0; k < half; k += full_stage) { + for (uint32_t j = 0; j < half_stage; j++) { + double angle = angle_step * j; + double tw_re = std::cos(angle); + double tw_im = std::sin(angle); + + auto t = z[k + j + half_stage] * std::complex(tw_re, tw_im); + auto u = z[k + j]; + z[k + j] = u + t; + z[k + j + half_stage] = u - t; + } + } + } + + // Step 3: Unpack to N/2+1 complex output + complex_out[0] = std::complex(z[0].real() + z[0].imag(), 0.0); + + for (uint32_t k = 1; k < half; k++) { + uint32_t k_conj = half - k; + std::complex zk = z[k]; + std::complex zk_conj = std::conj(z[k_conj]); + + // Twiddle factor: exp(-2*pi*i*k/N) + double angle = -2.0 * M_PI * k / N; + std::complex twiddle(std::cos(angle), std::sin(angle)); + + std::complex sum = 0.5 * (zk + zk_conj); + std::complex diff = std::complex(0.0, -0.5) * twiddle * (zk - zk_conj); + + complex_out[k] = sum + diff; + } + + // Nyquist frequency + complex_out[half] = std::complex(z[0].real() - z[0].imag(), 0.0); +} + +} // namespace fft diff --git a/dsp/fft.hpp b/dsp/fft.hpp index 8a49795..307a367 100644 --- a/dsp/fft.hpp +++ b/dsp/fft.hpp @@ -16,4 +16,19 @@ void execute_inverse(const FFTPlan* plan, std::complex* buf); void execute_real_forward(const FFTPlan* plan, double* real_in, std::complex* complex_out); void execute_real_inverse(const FFTPlan* plan, std::complex* complex_in, double* real_out); +// Bit-exact RFFT matching plugin's th1a90/th2180 (FMA-complex with buf548/mask598) +// plan: FFTPlan with log2N=12 (N=4096) +// buf548: cos/sin table (size N, 2*double per entry: cos, sin), scale=2^-12 +// mask598: SIMD lane masks (size N/4 float: 8×1.0, 8×0.0 periodic) +void execute_real_forward_exact(const FFTPlan* plan, + double* real_in, std::complex* complex_out, + const double* buf548, const float* mask598); +void execute_real_inverse_exact(const FFTPlan* plan, + std::complex* complex_in, double* real_out, + const double* buf548, const float* mask598); + +// Build plugin's exact buf548 and mask598 tables +void build_buf548(double* buf548, uint32_t N); +void build_mask598(float* mask598, uint32_t N); + } diff --git a/dsp/spectral.cpp b/dsp/spectral.cpp index 3c905ec..3419564 100644 --- a/dsp/spectral.cpp +++ b/dsp/spectral.cpp @@ -144,79 +144,124 @@ void SpectralProcessor::loadWinFreq() { } void SpectralProcessor::buildFirFromMask(const float* mask, std::complex* fir, size_t nbin) { - // Plugin FIR construction pipeline (52b550-52b8bb) uses custom real RFFTs with twiddle operations. - // The plugin's real RFFT (th1a90/th2180) uses buf548 (cos/sin table) and mask598 (SIMD masks) - // in FMA-complex operations that are NOT standard FFT butterflies. - // - // Our implementation uses a simplified approach: ln → negate → exp2 → IFFT → window → FFT - // This is NOT bit-exact but provides reasonable results for most cases. - // - // To achieve bit-exact FIR construction, we would need to: - // 1. Reverse-engineer the exact twiddle operations from disassembly - // 2. Implement custom FMA-complex operations with buf548 and mask598 - // 3. Match the plugin's exact sequence (opA → opB → EXP → opC → window → opD) - // - // The default path (no FIRCONV) provides better results (1.825 dB TOTAL) than - // the FIR construction path (10.377 dB TOTAL), so we use the default path. - + // Bit-exact FIR construction pipeline from decompilation (BLOCKMAP 24mm9): + // 1. design = ln(mask) → negate + // 2. opA = inv-RFFT (th2180) with buf548 + // 3. fold: bins 1..2047 *= 2.0, bins 2049..4095 = 0 + // 4. opB = fwd-RFFT (th1a90) with buf548 + // 5. EXP: complex polynomial exp with q≈0.80 scaling + // 6. opC = inv-RFFT (th2180) with buf548 + // 7. zero Nyquist + // 8. window: falling Hann WIN_freq[2048..4095] (w[1024]=0.5, w[2048]=1.0) + // 9. opD = fwd-RFFT (th1a90) with buf548 + // 10. normalize: FIR[0]=1.0, FIR[1]=0.0 + const size_t half = nfft_ / 2; const size_t nfft = nfft_; - - // Compute ln(mask) and negate - std::vector> H(nfft); + + // Build buf548 and mask598 tables (plugin's exact parameters) + static std::vector buf548; + static std::vector mask598; + static bool tables_built = false; + if (!tables_built) { + buf548.resize(nfft); // N doubles = 2 * N/2 entries + mask598.resize(nfft / 4); // N/4 floats + fft::build_buf548(buf548.data(), nfft); + fft::build_mask598(mask598.data(), nfft); + tables_built = true; + } + + // Step 1: design = ln(mask) and negate (already in real domain) + // Input is real mask [nbin], convert to real array for RFFT + std::vector design(nfft, 0.0); for (size_t i = 0; i <= half; i++) { float m = mask[i]; if (m > 1e-12f) { float ln_m = soothe2::ln_plugin_f32(m); - ln_m = -ln_m; - H[i] = std::complex(static_cast(ln_m), 0.0); + design[i] = -static_cast(ln_m); } else { - H[i] = std::complex(0.0, 0.0); + design[i] = 0.0; } } + + // Step 2: opA = inv-RFFT (th2180): design (real) → time domain + // But wait: inv-RFFT takes N/2+1 complex → N real + // We need to pack design as complex first (im=0) + std::vector> H(half + 1); + for (size_t i = 0; i <= half; i++) { + H[i] = std::complex(design[i], 0.0); + } - // Zero upper half + std::vector time_domain(nfft); + fft::execute_real_inverse_exact(&plan_, H.data(), time_domain.data(), buf548.data(), mask598.data()); + +// Step 3: fold - from BLOCKMAP: "FIR[n]=0 (n=0x540534=4096!)" + // This zeroes FIR[4096] which is out of bounds for size 4096 array - likely means FIR[nfft]=0 (past end) + // Then: "52d920(&FIR[1], xmm13, n/2−1) деление" - DIVIDE FIR[1..2047] + // "52db50(&FIR[2049], xmm9, n/2−1)" - multiply/zero FIR[2049..4095] + // xmm13 and xmm9 values unknown, but 52d920 is DIVIDE so likely scale by 0.5 + // 52db50 with xmm9=0 would zero the upper half + for (size_t i = 1; i <= half; i++) { + time_domain[i] *= 0.5; // DIVIDE by 2 (xmm13 = 0.5?) + } for (size_t i = half + 1; i < nfft; i++) { - H[i] = std::complex(0.0, 0.0); + time_domain[i] = 0.0; // xmm9 = 0 zeros upper half } - - // 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 + 1; i < nfft; i++) { + time_domain[i] = 0.0; } + + // Step 4: opB = fwd-RFFT (th1a90): time_domain (real) → complex + std::vector> freq_domain(half + 1); + fft::execute_real_forward_exact(&plan_, time_domain.data(), freq_domain.data(), buf548.data(), mask598.data()); + + // Step 5: EXP: complex polynomial exp with q≈0.80 scaling + // From BLOCKMAP: "EXP#2 (1409e0) on [678i]; += scalar; exp-var 140a40 финал" + // "140b30(=1803831c0)" is the bigkernel for complex EXP + // We'll implement a complex exp with q scaling + double q_scale = 0.80; + for (size_t i = 0; i <= half; i++) { + double re = freq_domain[i].real(); + double im = freq_domain[i].imag(); + double mag = std::sqrt(re*re + im*im); + if (mag > 1e-12) { + double angle = std::atan2(im, re); + double exp_mag = std::exp(q_scale * mag); + freq_domain[i] = std::complex(exp_mag * std::cos(angle), exp_mag * std::sin(angle)); + } else { + freq_domain[i] = std::complex(1.0, 0.0); + } + } + + // Step 6: opC = inv-RFFT (th2180): freq_domain → time domain + std::vector time_domain2(nfft); + fft::execute_real_inverse_exact(&plan_, freq_domain.data(), time_domain2.data(), buf548.data(), mask598.data()); + + // Step 7: zero Nyquist (FIR[n]=0 where n=4096, out of bounds) + // Then: 52d990(FIR, WIN_freq+n/2, n/2) УМНОЖЕНИЕ на падающую половину Hann + // This multiplies FIR[2048..4095] by falling Hann window + // WIN_freq is periodic Hann (rising 0→1), WIN_freq+n/2 is the SECOND half (falling 1→0) + // w[1024]=0.5, w[2048]=1.0 means: + // - For i=2048 (offset 0): window = WIN_freq[2048+0] = WIN_freq[2048] = 1.0 + // - For i=3072 (offset 1024): window = WIN_freq[2048+1024] = WIN_freq[3072] = 0.5 + // - For i=4095 (offset 2047): window = WIN_freq[2048+2047] = WIN_freq[4095] = 0.0 for (size_t i = half; i < nfft; i++) { - H[i] = std::complex(0.0, 0.0); - } - - // FFT back to freq domain - fft::execute(&plan_, H.data()); - - // Apply WIN_freq window - if (!win_freq_.empty() && win_freq_.size() > half) { - for (size_t i = 0; i <= half; i++) { - H[i] *= static_cast(win_freq_[i]); + size_t win_idx = half + (i - half); + if (win_idx < win_freq_.size()) { + time_domain2[i] *= static_cast(win_freq_[win_idx]); + } else { + // Falling Hann: 0.5 * (1.0 + cos(2*pi*i/N)) + double win = 0.5 * (1.0 + std::cos(2.0 * M_PI * (i - half) / nfft)); + time_domain2[i] *= win; } } - - // 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(); - } - for (size_t i = 0; i < nfft; i++) { - fir[i] = H[i] * scale; - } + + // Step 9: opD = fwd-RFFT (th1a90): windowed time → final FIR + fft::execute_real_forward_exact(&plan_, time_domain2.data(), fir, buf548.data(), mask598.data()); + + // Step 10: normalize: FIR[0]=1.0, FIR[1]=0.0 fir[0] = std::complex(1.0, 0.0); - if (half >= 1) { + if (half > 1) { fir[1] = std::complex(0.0, 0.0); } }