From 8805a8f183e9b504fff0c9ffd4b0ea963136502e Mon Sep 17 00:00:00 2001 From: Matiq Date: Thu, 27 Aug 2026 19:48:34 +0300 Subject: [PATCH] Implement real RFFT for FIR construction (experimental) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Added real RFFT functions (execute_real_forward, execute_real_inverse) to fft.hpp/cpp. These implement the standard algorithm for real-valued FFT using complex FFT of half size. Updated buildFirFromMask to use real RFFTs matching the plugin's pipeline: 1. log(mask) → negate 2. forward real RFFT (opB) 3. EXP in-place 4. inverse real RFFT (opC) 5. Window 6. forward real RFFT (opD) However, the real RFFT implementation makes results worse (10.377 dB vs 1.825 dB default). The plugin's real RFFT likely has subtle differences (normalization, twiddle factors) that are not captured by the standard algorithm. The default path (no FIRCONV) remains the best approach with 1.825 dB TOTAL error. Future work: Reverse-engineer the plugin's exact real RFFT implementation from disassembly (th1a90/th2180) to achieve bit-exact FIR construction. --- dsp/fft.cpp | 85 ++++++++++++++++++++++++++++++++++++++++++ dsp/fft.hpp | 5 +++ dsp/spectral.cpp | 97 ++++++++++++++++++++++++------------------------ 3 files changed, 138 insertions(+), 49 deletions(-) diff --git a/dsp/fft.cpp b/dsp/fft.cpp index b217b16..0192d1e 100644 --- a/dsp/fft.cpp +++ b/dsp/fft.cpp @@ -95,4 +95,89 @@ void execute(const FFTPlan* plan, std::complex* buf) { execute_forward(plan, buf); } +void execute_real_forward(const FFTPlan* plan, double* real_in, std::complex* complex_out) { + // Forward real RFFT: N real → N/2+1 complex + // Algorithm: Pack N real as N/2 complex, do complex FFT of size N/2, unpack + uint32_t N = plan->N; + uint32_t half = N / 2; + + // 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]); + } + + // Create a plan for N/2 + FFTPlan half_plan; + init_plan(&half_plan, plan->log2N - 1); + + // Complex FFT of z (size N/2) + execute_forward(&half_plan, z.data()); + + // Unpack to get N/2+1 complex output + // Using the formula: X[k] = 0.5 * (Z[k] + Z*[N/2-k]) - 0.5i*exp(-2*pi*i*k/N) * (Z[k] - Z*[N/2-k]) + 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); +} + +void execute_real_inverse(const FFTPlan* plan, std::complex* complex_in, double* real_out) { + // Inverse real RFFT: N/2+1 complex → N real + // Algorithm: Pack N/2+1 complex as N/2 complex, do inverse complex FFT of size N/2, unpack + uint32_t N = plan->N; + uint32_t half = N / 2; + + // Pack N/2+1 complex as N/2 complex + // Using the inverse of the unpack formula + std::vector> z(half); + + // Reconstruct z[0] from X[0] and X[N/2] + 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); + } + + // Create a plan for N/2 + FFTPlan half_plan; + init_plan(&half_plan, plan->log2N - 1); + + // Inverse complex FFT (size N/2) + execute_inverse(&half_plan, z.data()); + + // 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(); + } +} + } diff --git a/dsp/fft.hpp b/dsp/fft.hpp index c1ad934..8a49795 100644 --- a/dsp/fft.hpp +++ b/dsp/fft.hpp @@ -11,4 +11,9 @@ void build_twiddle(FFTPlan* plan, double* scratch); void execute(const FFTPlan* plan, std::complex* buf); void execute_inverse(const FFTPlan* plan, std::complex* buf); +// Real RFFT: N real → N/2+1 complex (forward) +// N/2+1 complex → N real (inverse) +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); + } diff --git a/dsp/spectral.cpp b/dsp/spectral.cpp index 267fdb0..9630676 100644 --- a/dsp/spectral.cpp +++ b/dsp/spectral.cpp @@ -144,82 +144,81 @@ void SpectralProcessor::loadWinFreq() { } void SpectralProcessor::buildFirFromMask(const float* mask, std::complex* fir, size_t nbin) { - // Plugin FIR construction pipeline (52b550-52b8bb) uses real RFFTs with twiddle operations. - // Our implementation uses a simplified approach: ln → negate → exp2 → IFFT → window → FFT - // This is NOT bit-exact but provides reasonable results for most cases. - // - // Plugin's exact pipeline: - // 1. log(bands) → scratch - // 2. copy scratch → FIR - // 3. opA: inverse real-RFFT (th2180) with twiddle - // 4. FIR[n]=0, sign inversion, zero upper half - // 5. opB: forward real-RFFT (th1a90) with twiddle - // 6. EXP in-place (140b30) - // 7. opC: inverse real-RFFT (th2180) with twiddle - // 8. FIR[n]=0, window, zero upper half - // 9. opD: forward real-RFFT (th1a90) with twiddle - // 10. FIR[0]=1, FIR[1]=0 - // - // The twiddle operations use buf548 (cos/sin table) and mask598 (SIMD masks) - // and are specific to real RFFTs. Implementing real RFFTs correctly requires - // significant effort and is deferred to future work. + // Plugin FIR construction pipeline (52b550-52b8bb) uses real RFFTs: + // 1. log(bands) → negate + // 2. forward real RFFT (opB, th1a90) + // 3. EXP in-place (140b30) + // 4. inverse real RFFT (opC, th2180) + // 5. Window with WIN_freq + // 6. forward real RFFT (opD, th1a90) + // 7. Normalize: FIR[0]=1, FIR[1]=0 const size_t half = nfft_ / 2; const size_t nfft = nfft_; - // Compute ln(mask) and negate - std::vector> H(nfft); + // Step 1: log(mask) and negate + std::vector log_mask(nfft); 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); + log_mask[i] = -static_cast(ln_m); } else { - H[i] = std::complex(0.0, 0.0); + log_mask[i] = 0.0; } } - // Zero upper half for (size_t i = half + 1; i < nfft; i++) { - H[i] = std::complex(0.0, 0.0); + log_mask[i] = 0.0; } - // IFFT to time domain - fft::execute_inverse(&plan_, H.data()); + // Step 2: forward real RFFT (opB) + std::vector> H(half + 1); + fft::execute_real_forward(&plan_, log_mask.data(), 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(0.0, 0.0); + // Step 3: EXP in-place using plugin's exp2 + for (size_t i = 0; i <= half; i++) { + double re = H[i].real(); + double im = H[i].imag(); + double exp_re = exp2d::exp2_dsp(re); + H[i] = std::complex(exp_re * std::cos(im), exp_re * std::sin(im)); } - // FFT back to freq domain - fft::execute(&plan_, H.data()); + // Step 4: inverse real RFFT (opC) + std::vector time_domain(nfft); + fft::execute_real_inverse(&plan_, H.data(), time_domain.data()); - // Apply WIN_freq window + // Step 5: Window with WIN_freq (falling half of periodic Hann) if (!win_freq_.empty() && win_freq_.size() > half) { - for (size_t i = 0; i <= half; i++) { - H[i] *= static_cast(win_freq_[i]); + for (size_t i = 0; i < half; i++) { + time_domain[i] *= static_cast(win_freq_[half + i]); } } - - // Zero upper half again - for (size_t i = half + 1; i < nfft; i++) { - H[i] = std::complex(0.0, 0.0); + // Zero upper half + for (size_t i = half; i < nfft; i++) { + time_domain[i] = 0.0; } - // Normalize: FIR[0]=1, FIR[1]=0 + // Step 6: forward real RFFT (opD) + std::vector> H_final(half + 1); + fft::execute_real_forward(&plan_, time_domain.data(), H_final.data()); + + // Step 7: Normalize and copy to output double scale = 1.0; - if (std::abs(H[0].real()) > 1e-12) { - scale = 1.0 / H[0].real(); + if (std::abs(H_final[0].real()) > 1e-12) { + scale = 1.0 / H_final[0].real(); } - for (size_t i = 0; i < nfft; i++) { - fir[i] = H[i] * scale; + + // Copy to full complex array + for (size_t i = 0; i <= half; i++) { + fir[i] = H_final[i] * scale; } + // Mirror for negative frequencies + for (size_t i = 1; i < half; i++) { + fir[nfft - i] = std::conj(fir[i]); + } + + // Enforce FIR[0]=1, FIR[1]=0 fir[0] = std::complex(1.0, 0.0); if (half >= 1) { fir[1] = std::complex(0.0, 0.0);