真正的 FFT 输出

Question

我已经使用 Kiss fft 库在 at32ucb 系列 ucontroller 中实现了 fft，目前正在努力处理 fft 的输出。 我的目的是分析来自压电扬声器的声音。 目前，发声器的频率是 420Hz，我成功地从 fft 输出（用示波器交叉检查）得到。但是，如果我将函数发生器波形放入系统，则输出频率仅为预期的一半。 我怀疑是我弄错了频率仓计算公式；当前使用，fft_peak_magnitude_index*采样频率/fft_size。 我的输入是真实的，并且是真实的 fft。 （输出样本 = N/2） 并且还在 fft 之前进行 iir 过滤和窗口化。 任何建议都会有很大帮助！

            // IIR filter calculation, n = 256 fft points       
        for (ctr=0; ctr<n; ctr++)
        {       
            // filter calculation
            y[ctr] = num_coef[0]*x[ctr];
            y[ctr] += (num_coef[1]*x[ctr-1]) - (den_coef[1]*y[ctr-1]);
            y[ctr] += (num_coef[2]*x[ctr-2]) - (den_coef[2]*y[ctr-2]);
            y1[ctr] = y[ctr] - 510; //eliminate dc offset

            // hamming window
            hamming[ctr] = (0.54-((0.46) * cos(2*M_PI*ctr/n)));
            window[ctr] = hamming[ctr]*y1[ctr];

            fft_input[ctr].r = window[ctr];
            fft_input[ctr].i = 0;
            fft_output[ctr].r = 0;
            fft_output[ctr].i = 0;
        }


        kiss_fftr_cfg fftConfig = kiss_fftr_alloc(n,0,NULL,NULL);
        kiss_fftr(fftConfig, (kiss_fft_scalar * )fft_input, fft_output);

        peak = 0;
        freq_bin = 0;       
        for (ctr=0; ctr<n1; ctr++)
            {   
                fft_mag[ctr] = 10*(sqrt((fft_output[ctr].r * fft_output[ctr].r) + (fft_output[ctr].i * fft_output[ctr].i)))/(0.5*n);

                if(fft_mag[ctr] > peak)
                {
                    peak = fft_mag[ctr];
                    freq_bin = ctr;
                }

            frequency = (freq_bin*(10989/n)); // 10989 is the sampling freq
                //************************************
                //Usart write
                char filtResult[10];
                //sprintf(filtResult, "%04d %04d %04d\n", (int)peak, (int)freq_bin, (int)frequency);
                sprintf(filtResult, "%04d  %04d  %04d\n", (int)x[ctr], (int)fft_mag[ctr], (int)frequency);
                char c;
                char *ptr = &filtResult[0];
                do
                {
                    c = *ptr;
                    ptr++;
                    usart_bw_write_char(&AVR32_USART2, (int)c);
                    // sendByte(c);

                } while (c != '\n');
            }

Answer 1

主要问题可能是您声明fft_input 的方式。 基于your previous question，您将fft_input 分配为kiss_fft_cpx 的数组。另一方面，函数kiss_fftr 需要一个标量数组。通过将输入数组转换为 kiss_fft_scalar ：

kiss_fftr(fftConfig, (kiss_fft_scalar * )fft_input, fft_output);

KissFFT 本质上会看到一个实值数据数组，其中每隔一个样本包含零（您填写为虚部的内容）。这实际上是原始信号的上采样版本（尽管没有插值），即采样率有效两倍的信号（在您的 freq_bin 到 frequency 转换中未考虑）。要解决这个问题，我建议您将数据打包到 kiss_fft_scalar 数组中：

kiss_fft_scalar fft_input[n];
...
for (ctr=0; ctr<n; ctr++)
{       
    ...
    fft_input[ctr] = window[ctr];
    ...
}
kiss_fftr_cfg fftConfig = kiss_fftr_alloc(n,0,NULL,NULL);
kiss_fftr(fftConfig, fft_input, fft_output);

另请注意，在寻找峰值幅度时，您可能只对最终的最大峰值感兴趣，而不是运行最大值。因此，您可以将循环限制为仅计算峰值（如果需要，在以下 sprintf 语句中使用 freq_bin 而不是 ctr 作为数组索引）：

for (ctr=0; ctr<n1; ctr++)
{   
  fft_mag[ctr] = 10*(sqrt((fft_output[ctr].r * fft_output[ctr].r) + (fft_output[ctr].i * fft_output[ctr].i)))/(0.5*n);

  if(fft_mag[ctr] > peak)
  {
    peak = fft_mag[ctr];
    freq_bin = ctr;
  }
} // close the loop here before computing "frequency"

最后，在计算与幅度最大的 bin 相关的频率时，您需要确保使用浮点运算完成计算。如果我怀疑n 是一个整数，那么您的公式将使用整数运算执行10989/n 因子，从而导致截断。这可以通过以下方式简单地解决：

frequency = (freq_bin*(10989.0/n)); // 10989 is the sampling freq