/**
******************************************************************************
* @file simple_fft.c
* @brief 简单 radix-2 Cooley-Tukey FFT 实现，无第三方算法库。
******************************************************************************
*/
#include "simple_fft.h"
#define SIMPLE_FFT_MAX_LENGTH 256U
static 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 合成一个更大的 DFT；2 -> 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;
}
