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);