From 994c8e29bac4a5969cffb9eb913d2e74692bc71c Mon Sep 17 00:00:00 2001 From: skal Date: Tue, 24 Mar 2026 07:52:25 +0100 Subject: fix(fft): replace iterative twiddle with direct cosf/sinf, add tests A-E fft_radix2 now computes wr=cosf(angle*k)/wi=sinf(angle*k) directly per k, eliminating float drift over long iteration runs. Iterative approach documented in comment for reference. Tests A-E added (bit-reverse, small-N DFT, twiddle drift, DCT small/large N). arrays_match tolerance reverted to 5e-3. TODO.md updated. handoff(Gemini): fft twiddle fix complete, 38/38 tests passing. --- src/audio/fft.cc | 25 +++++++++++++------------ src/audio/fft.h | 3 +++ 2 files changed, 16 insertions(+), 12 deletions(-) (limited to 'src/audio') diff --git a/src/audio/fft.cc b/src/audio/fft.cc index 6b8ba8e..6512fdc 100644 --- a/src/audio/fft.cc +++ b/src/audio/fft.cc @@ -9,7 +9,7 @@ // Bit-reversal permutation (in-place) // Reorders array elements by reversing their binary indices -static void bit_reverse_permute(float* real, float* imag, size_t N) { +void bit_reverse_permute(float* real, float* imag, size_t N) { size_t temp_bits = N; size_t num_bits = 0; while (temp_bits > 1) { @@ -48,13 +48,19 @@ static void fft_radix2(float* real, float* imag, size_t N, int direction) { const size_t half_stage = stage_size / 2; const float angle = direction * 2.0f * PI / stage_size; - // Precompute twiddle factors for this stage - float wr = 1.0f; - float wi = 0.0f; - const float wr_delta = cosf(angle); - const float wi_delta = sinf(angle); - for (size_t k = 0; k < half_stage; k++) { + // Direct twiddle factor: numerically stable, no drift. + // Faster iterative alternative (~4 mul+add vs 2 transcendentals) but + // accumulates float drift over many iterations (e.g. 256 steps for + // N=512 final stage). Shown here for reference: + // float wr = 1.0f, wi = 0.0f; + // const float wr_delta = cosf(angle), wi_delta = sinf(angle); + // // per k: const float wr_old = wr; + // // wr = wr_old * wr_delta - wi * wi_delta; + // // wi = wr_old * wi_delta + wi * wr_delta; + const float wr = cosf(angle * (float)k); + const float wi = sinf(angle * (float)k); + // Apply butterfly to all groups at this stage for (size_t group_start = k; group_start < N; group_start += stage_size) { const size_t i = group_start; @@ -70,11 +76,6 @@ static void fft_radix2(float* real, float* imag, size_t N, int direction) { real[i] = real[i] + temp_real; imag[i] = imag[i] + temp_imag; } - - // Update twiddle factor for next k (rotation) - const float wr_old = wr; - wr = wr_old * wr_delta - wi * wi_delta; - wi = wr_old * wi_delta + wi * wr_delta; } } } diff --git a/src/audio/fft.h b/src/audio/fft.h index df37ad5..6a54742 100644 --- a/src/audio/fft.h +++ b/src/audio/fft.h @@ -8,6 +8,9 @@ #include +// Bit-reversal permutation (in-place). Exposed for testing. +void bit_reverse_permute(float* real, float* imag, size_t N); + // Forward FFT: Time domain → Frequency domain // Input: real[] (length N), imag[] (length N, can be zeros) // Output: real[] and imag[] contain complex frequency bins -- cgit v1.2.3