带或不带窗口的 KISS FFT 输出

Question

我目前正在尝试使用 kiss fft 将 fft 实现到 avr32 微控制器中以进行信号处理。我的输出有一个奇怪的问题。基本上，我将 ADC 样本（使用函数发生器进行测试）传递到 fft（实际输入，256 n 大小）并且检索到的输出对我来说很有意义。但是，如果我将 Hamming window 应用于 ADC 样本，然后将它们传递给 FFT，则峰值幅度的频率仓是错误的（并且与之前没有 windowing 的结果不同）。 ADC 样本具有 DC 偏移，因此我消除了偏移，但它仍然不适用于 windowed 样本。

下面是通过rs485输出的前几个值。第一列是没有 window 的 fft 输出，而第二列是带有 window 的输出。从第 1 列开始，峰值位于第 6 行 (6 x fs (10.5kHz) / 0.5N) 给了我正确的输入频率结果，其中第 2 列在第 2 行（直流仓除外）有一个峰值幅度，这对我来说没有意义. 任何建议都会有所帮助。提前致谢。

    488         260   //dc bin
    5            97
    5            41
    5            29  
    4            26
    10           35
    133          76
    33           28
    21           6
    17           3

kiss_fft_scalar zero;
memset(&zero,0,sizeof(zero));
kiss_fft_cpx fft_input[n];
kiss_fft_cpx fft_output[n];

for(ctr=0; ctr<n; ctr++)
{
    fft_input[ctr].r = zero;
    fft_input[ctr].i = zero;
    fft_output[ctr].r =zero;
    fft_output[ctr].i = zero;
}

// IIR filter calculation

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] - 500;
    // hamming window
    hamming[ctr] = (0.54-((0.46) * cos(2*PI*ctr/256)));
    window[ctr] = hamming[ctr]*y[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);


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

    //Usart write
    char filtResult[10];
    sprintf(filtResult, "%04d %04d\n", (int)x[ctr], (int)fft_mag[ctr]);
    //sprintf(filtResult, "%04d %04d\n", (int)x[ctr], (int)window[ctr]);
    char c;
    char *ptr = &filtResult[0];
    do
    {   
        c = *ptr;
        ptr++;
        usart_bw_write_char(&AVR32_USART2, (int)c);
        // sendByte(c);

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

kiss_fft_cleanup();
free(fftConfig);

Answer 1

频域输出说明

在频域中，矩形和汉明windows看起来像：

如您所知，在时域中乘以 window 对应于频域中的卷积，这实质上是将信号的能量分布在通常称为 spectral leakage。对于您选择的特定 windows（如上频域所示），您可能会注意到汉明 window 在主瓣外传播的能量要少得多，但主瓣稍微宽一点比矩形 window.

因此，当使用 Hamming window 时，DC 能量的尖峰最终会传播到 bin 0 并进入 bin 1。与其说你在 bin 1 中有一个强峰。事实上，如果你绘制你提供的数据，你应该看到你在索引 6 看到的尖峰实际上仍然以相同的频率存在并且不使用汉明 window:

如果您想消除周围 bin 中的 DC 尖峰和泄漏，要么消除数据中的偏差（基本上应用陷波滤波器），要么您将不得不忽略更多的低频 bin在寻找你的 "first strong spike".

时

过滤问题

最后，请注意 IIR 滤波器的实现方式也存在一个问题，即数组 x 和 y 将在 ctr==0 和 ctr==1 （除非你做了一些特殊的规定，并且 x & y 实际上是指针，从分配的缓冲区开始处有一个偏移量）。这可能会影响使用和不使用 window 的结果。如果您只过滤单个数据块，一个常见的假设是较早的样本为零。在这种情况下，您可以使用以下方法避免越界索引：

// filter calculation
y[ctr] = num_coef[0]*x[ctr];

if (ctr>=1)
{
  y[ctr] += (num_coef[1]*x[ctr-1]) - (den_coef[1]*y[ctr-1]);
}
if (ctr>=2)
{
  y[ctr] += (num_coef[2]*x[ctr-2]) - (den_coef[2]*y[ctr-2]);
}

另一方面，如果您想过滤 n 个样本的多个块，则必须记住前一个块的最后几个样本。这可以通过分配比块大小稍大的缓冲区来完成：

x = malloc((n+2)*sizeof(kiss_fft_scalar));
y = malloc((n+2)*sizeof(kiss_fft_scalar));
// initialize "past samples" for the first block, assuming data was zero
x[0] = x[1] = 0;
y[0] = y[1] = 0;

然后你可以在这些缓冲区中使用偏移量。索引 0 和 1 代表过去的样本，而索引 2 的缓冲区的其余部分填充了当前输入数据块。这导致以下稍微修改过的过滤代码：

// filter calculation
y[ctr+2] = num_coef[0]*x[ctr+2];

y[ctr+2] += (num_coef[1]*x[ctr+1]) - (den_coef[1]*y[ctr+1]);
y[ctr+2] += (num_coef[2]*x[ctr]) - (den_coef[2]*y[ctr]);

// hamming window
hamming[ctr] = (0.54-((0.46) * cos(2*PI*ctr/256)));
window[ctr] = hamming[ctr]*y[ctr+2];

最后，在每个块的末尾，您必须更新索引 0 和 1 处的 "past samples"，当前块的最后样本准备好处理下一个输入块：

// remember last 2 samples of block
x[0] = x[n-2];
x[1] = x[n-1];
y[0] = y[n-2];
y[1] = y[n-1];

带或不带窗口的 KISS FFT 输出

KISS FFT output with or without windowing

signal-processing

windowing

dft

kissfft

频域输出说明

过滤问题