2026-04-29 06:57:46 +03:00

83 lines
2.2 KiB
Plaintext

//
// Tiny in-place radix-2 Cooley-Tukey FFT. Module-level scratch buffers so we
// don't allocate per frame.
//
// FFT_SIZE is a compile-time constant (power of two). 1024 at 44.1 kHz gives
// ~23 ms windows and ~43 Hz bin resolution — a good balance for the
// spectrum visualizer.
//
FFT_SIZE :: 1024;
HALF_FFT :: FFT_SIZE / 2;
fft_window: [FFT_SIZE] float;
fft_re: [FFT_SIZE] float;
fft_im: [FFT_SIZE] float;
fft_initted: bool;
fft_init :: () {
if fft_initted return;
// Hann window: smooth attenuation at edges, no spectral leakage from
// arbitrary-phase sample boundaries.
for 0..FFT_SIZE-1 {
t := cast(float) it / cast(float) (FFT_SIZE - 1);
fft_window[it] = 0.5 - 0.5 * cos(2 * PI * t);
}
fft_initted = true;
}
//
// In-place complex FFT on (fft_re, fft_im).
//
fft :: () {
n := FFT_SIZE;
// Bit-reverse permutation.
j := 0;
for i: 0..n-2 {
if i < j {
t_re := fft_re[i]; t_im := fft_im[i];
fft_re[i] = fft_re[j]; fft_im[i] = fft_im[j];
fft_re[j] = t_re; fft_im[j] = t_im;
}
bit := n >> 1;
while j & bit {
j ^= bit;
bit >>= 1;
}
j ^= bit;
}
// Cooley-Tukey butterflies.
seg_len := 2;
while seg_len <= n {
ang := -2.0 * PI / cast(float) seg_len;
w_re := cos(ang);
w_im := sin(ang);
i := 0;
while i < n {
tw_re := 1.0;
tw_im := 0.0;
half := seg_len / 2;
for k: 0..half-1 {
a := i + k;
b := i + k + half;
v_re := fft_re[b] * tw_re - fft_im[b] * tw_im;
v_im := fft_re[b] * tw_im + fft_im[b] * tw_re;
fft_re[b] = fft_re[a] - v_re;
fft_im[b] = fft_im[a] - v_im;
fft_re[a] = fft_re[a] + v_re;
fft_im[a] = fft_im[a] + v_im;
next_re := tw_re * w_re - tw_im * w_im;
next_im := tw_re * w_im + tw_im * w_re;
tw_re = next_re;
tw_im = next_im;
}
i += seg_len;
}
seg_len *= 2;
}
}