ESP32-S3 FFT计算输入信号最大频率结果异常求助

问题背景 #

在ESP32-S3项目中，基于ESP-IDF的FFT DSP示例实现输入信号最大频率计算。生成包含3Hz、5Hz的正弦合成信号（5Hz幅度更大），采样率1000Hz，采样时长3秒（3000个采样点），FFT大小设为4096，但最终计算得到的最大频率约为612Hz，远高于预期的5Hz。

核心代码 #

信号生成代码 #

double composite_signal(double t, SignalComponent *component) {
    double signal = 0;
    
    // Summing sine components
    for (int i = 0; i < component->num_sin_components; i++) {
        signal += component->sin_amplitudes[i] * sin(2 * M_PI * component->sin_frequencies[i] * t);
    }
    
    return signal;
}

SignalComponent component = {
    2, // Number of sine components
    {3,5}, // Sine frequencies in HZ
    {2,4} // Sine amplitudes
};

// 在main函数中
for (int i = 0; i < NUM_SAMPLES; i++) {
    double t = (double)i / SAMPLING_RATE;
    wave[i] = composite_signal(t, &component);
}

参数定义 #

#define SAMPLING_RATE 1000 // Hz
#define DURATION 3 // seconds
#define NUM_SAMPLES (SAMPLING_RATE * DURATION)
#define TH 3 //threshold for magnitude, max frequency of the fft

FFT处理代码 #

dsps_wind_hann_f32(wind, NUM_SAMPLES);
// Convert two input vectors to one complex vector
for (int i = 0 ; i < NUM_SAMPLES ; i++) {
    y_cf[i * 2] = wave[i] * wind[i];
    y_cf[i * 2 + 1] = 0; // I don't have a second vector as in the example
}
// FFT
unsigned int start_b = dsp_get_cpu_cycle_count();
dsps_fft2r_fc32(y_cf, NUM_SAMPLES);
unsigned int end_b = dsp_get_cpu_cycle_count();
// Bit reverse
dsps_bit_rev_fc32(y_cf, NUM_SAMPLES);
// Convert one complex vector to two complex vectors
dsps_cplx2reC_fc32(y_cf, NUM_SAMPLES);

//RETRIEVE THE MAX FREQ OF THE SIGNAL
int maxI = 0;
double maxM = 0;
for (int i = 0; i < NUM_SAMPLES; i++)
{
    if(y_cf[i]>TH && y_cf[i]>maxM){
        maxI = i;
        maxM = y_cf[i];
    }
}
ESP_LOGI(TAG, "THE MAX Index is %i", maxI);
float maxF = maxI;
maxF *= SAMPLING_RATE / (float) NUM_SAMPLES; //max frequency retrival

实际输出 #

I (22594) main: THE MAX Index is 1838
I (22944) main: The maximum frequency of the signal is 612.666687, and the relative magnitude value is 564.219421.

问题根源分析 #

1. FFT输入长度不符合要求
   ESP-IDF的dsps_fft2r_fc32函数要求输入的FFT长度必须是2的幂次（如1024、2048、4096）。你的NUM_SAMPLES=3000不是2的幂次，强行用这个长度执行FFT会导致计算逻辑混乱，输出错误的频率点。

2. 幅度计算逻辑错误
   dsps_cplx2reC_fc32执行后，y_cf数组中存储的是每个频率点的实部和虚部交替数据，并非直接是幅度值。你直接用y_cf[i]比较大小，实际上是拿实部或虚部作为幅度，完全错误，导致选中了错误的索引。

3. 频率索引范围超出有效区间
   根据奈奎斯特采样定理，采样率1000Hz时，有效频率范围是0_{500Hz，对应FFT索引范围是0到`NUM_FFT/2`。你遍历整个`NUM_SAMPLES`，会进入镜像频率区间（500Hz}1000Hz），这部分是无效的镜像数据，所以得到了612Hz这种错误结果。

修正方案 #

步骤1：调整FFT长度为2的幂次 #

将FFT长度设为2的幂次，对3000个采样点补零到4096：

#define NUM_SAMPLES 3000
#define NUM_FFT 4096  // 2的幂次，符合FFT函数要求

步骤2：正确计算幅度值 #

遍历FFT结果时，计算每个频率点的幅度平方（避免开方，提升效率）来比较大小。

步骤3：限制索引到有效频率区间 #

只遍历0到NUM_FFT/2的索引，避免进入镜像频率区域。

修正后的FFT处理代码 #

// 初始化汉宁窗，长度为NUM_SAMPLES
dsps_wind_hann_f32(wind, NUM_SAMPLES);

// 初始化复数数组，长度为2*NUM_FFT（每个复数占两个float）
float y_cf[2 * NUM_FFT] = {0};

// 填充采样数据并加窗，剩余部分补零
for (int i = 0 ; i < NUM_SAMPLES ; i++) {
    y_cf[i * 2] = wave[i] * wind[i];
    y_cf[i * 2 + 1] = 0;
}
// 超出采样点的部分已经初始化为0，无需额外处理

// 执行FFT（使用NUM_FFT作为长度，而非NUM_SAMPLES）
unsigned int start_b = dsp_get_cpu_cycle_count();
dsps_fft2r_fc32(y_cf, NUM_FFT);
unsigned int end_b = dsp_get_cpu_cycle_count();

// 位反转
dsps_bit_rev_fc32(y_cf, NUM_FFT);

// 转换为实部虚部分开的格式
dsps_cplx2reC_fc32(y_cf, NUM_FFT);

// 查找最大幅度的频率点
int maxI = 0;
double maxM = 0;
// 只遍历有效频率区间：0到NUM_FFT/2
for (int i = 0; i < NUM_FFT / 2; i++)
{
    // 获取实部和虚部
    float real = y_cf[i];
    float imag = y_cf[i + NUM_FFT];
    // 计算幅度平方（比开方更快，不影响大小比较）
    double mag_sq = real * real + imag * imag;
    if(mag_sq > TH * TH && mag_sq > maxM){
        maxI = i;
        maxM = mag_sq;
    }
}
ESP_LOGI(TAG, "THE MAX Index is %i", maxI);
// 计算对应频率：索引 * (采样率 / FFT长度)
float maxF = (float)maxI * SAMPLING_RATE / NUM_FFT;
ESP_LOGI(TAG, "The maximum frequency of the signal is %f, and the relative magnitude value is %f", maxF, sqrt(maxM));

额外优化建议 #

• 如果不需要3秒的采样时长，可以直接将NUM_SAMPLES设为4096，这样不需要补零，计算效率更高。
• 汉宁窗的长度要和采样点长度一致，补零部分不需要加窗，避免引入额外误差。

内容的提问来源于stack exchange，提问作者Mimmo