diff options
Diffstat (limited to 'src/audio')
| -rw-r--r-- | src/audio/fft.cc | 25 | ||||
| -rw-r--r-- | src/audio/fft.h | 3 |
2 files changed, 16 insertions, 12 deletions
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 <cstddef> +// 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 |