simple_fft.c

/********************************************************************************* file simple_fft.c* brief 简单 radix-2 Cooley-Tukey FFT 实现无第三方算法库。*******************************************************************************/#include simple_fft.h#define SIMPLE_FFT_MAX_LENGTH 256Ustatic uint8_t SimpleFFT_IsPowerOfTwo(uint16_t value);static uint16_t SimpleFFT_ReverseBits(uint16_t value, uint8_t bit_count);static uint8_t SimpleFFT_Log2(uint16_t value);uint8_t SimpleFFT_ForwardReal(const float *input, SimpleFFT_Complex_t *output, uint16_t length){uint16_t i;uint16_t block_size;uint8_t bits;if ((input 0) || (output 0) || (length SIMPLE_FFT_MAX_LENGTH) ||(SimpleFFT_IsPowerOfTwo(length) 0U)){return 0U;}/* FFT 首先做位倒序重排输入第 i 个样本写入输出的 bit_reverse(i) 位置。 */bits SimpleFFT_Log2(length);for (i 0U; i length; i){uint16_t reversed SimpleFFT_ReverseBits(i, bits);output[reversed].real input[i];output[reversed].imag 0.0f; /* 实数输入的虚部固定为 0。 */}/* 每一级把两个小 DFT 合成一个更大的 DFT2 - 4 - ... - length。 */for (block_size 2U; block_size length; block_size 1U){float step_real;float step_imag;uint16_t half (uint16_t)(block_size 1U);uint16_t start;uint16_t j;/** 本级旋转因子 exp(-j*2*pi/block_size)。只支持最大 256 点* 所有需要的 cos/sin 常量均列出避免调用 math 库。*/switch (block_size){case 2U: step_real -1.0f; step_imag 0.0f; break;case 4U: step_real 0.0f; step_imag -1.0f; break;case 8U: step_real 0.70710678f; step_imag -0.70710678f; break;case 16U: step_real 0.92387953f; step_imag -0.38268343f; break;case 32U: step_real 0.98078528f; step_imag -0.19509032f; break;case 64U: step_real 0.99518473f; step_imag -0.09801714f; break;case 128U: step_real 0.99879546f; step_imag -0.04906767f; break;case 256U: step_real 0.99969882f; step_imag -0.02454123f; break;default: return 0U;}for (start 0U; start length; start (uint16_t)(start block_size)){float w_real 1.0f;float w_imag 0.0f;for (j 0U; j half; j){uint16_t top (uint16_t)(start j);uint16_t bottom (uint16_t)(top half);float t_real w_real * output[bottom].real - w_imag * output[bottom].imag;float t_imag w_real * output[bottom].imag w_imag * output[bottom].real;float top_real output[top].real;float top_imag output[top].imag;float next_w_real;output[top].real top_real t_real;output[top].imag top_imag t_imag;output[bottom].real top_real - t_real;output[bottom].imag top_imag - t_imag;/* w w * step得到本级下一个蝶形运算的旋转因子。 */next_w_real w_real * step_real - w_imag * step_imag;w_imag w_real * step_imag w_imag * step_real;w_real next_w_real;}}}return 1U;}float SimpleFFT_GetBinEnergy(const SimpleFFT_Complex_t *spectrum, uint16_t bin){if (spectrum 0){return 0.0f;}return spectrum[bin].real * spectrum[bin].real spectrum[bin].imag * spectrum[bin].imag;}static uint8_t SimpleFFT_IsPowerOfTwo(uint16_t value){return ((value 2U) ((value (uint16_t)(value - 1U)) 0U)) ? 1U : 0U;}static uint8_t SimpleFFT_Log2(uint16_t value){uint8_t bits 0U;while (value 1U){value 1U;bits;}return bits;}static uint16_t SimpleFFT_ReverseBits(uint16_t value, uint8_t bit_count){uint16_t result 0U;while (bit_count 0U){result (uint16_t)((result 1U) | (value 1U));value 1U;bit_count--;}return result;}