带窗函数或不带窗函数的KISS FFT输出

Question

带窗函数或不带窗函数的KISS FFT输出

5

我目前正在尝试将fft应用于avr32微控制器进行信号处理，使用的是kiss fft。但在输出方面遇到了奇怪的问题。基本上，我将ADC样本（用函数发生器测试）传递给fft（实际输入，256n大小），并且检索到的输出对我来说是有意义的。然而，如果我将汉明窗应用于ADC样本，然后将其传递给FFT，峰值幅度的频率处于错误位置（与未应用窗口的先前结果不同）。 ADC样本具有DC偏移，因此我消除了该偏移量，但仍无法使用带窗口的样本。以下是通过rs485输出的前几个输出值。第一列是未使用窗口的fft输出，而第二列是使用窗口的输出。从第1列开始，峰值位于第6行（6 x fs（10.5kHz）/0.5N），为我提供了正确的输入频率结果，而第2列在第2行具有峰值幅度（除DC bin外），这对我来说没有意义。任何建议都将有所帮助。提前致谢。

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);

- Jin

1个回答

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- SleuthEye · Accepted Answer

频域输出澄清

在频域中，矩形窗口和Hamming窗口的样子如下：

频域中的矩形和Hamming窗口

您可能已经知道，在时间域中通过窗口进行乘法运算对应于在频率域中进行卷积运算，这实际上会将信号的能量分散到多个频率 bin 中，通常被称为spectral leakage。对于您选择的特定窗口（如上所示），您可能会注意到 Hamming 窗口在主瓣外部分散的能量要少得多，但其主瓣略微比矩形窗口宽一些。

因此，使用 Hamming 窗口时直流能量的峰值最终会扩散到第 0 个 bin 之外，并进入第 1 个 bin。实际上，并没有出现第 1 个 bin 中很强的峰值。事实上，如果绘制您提供的数据，您应该看到在使用和不使用 Hamming 窗口时，您观察到的在索引 6 处的尖峰实际上仍然存在于相同的频率处：

如果要消除直流尖峰和周围 bin 中的泄漏，可以删除数据中的偏置（实质上应用陷波滤波器），或者在查找“第一个强峰”时还需要忽略更多低频 bin。

过滤问题

最后，请注意 IIR 滤波器实现方式存在问题，其中当 ctr==0 和 ctr==1 时会超出数组 x 和 y 的索引范围（除非您已经作出了某些特殊规定，并且 x 和 y 实际上是指针，其与分配缓冲区的开始处有一个偏移量）。这可能会影响使用和不使用窗口的结果。如果您只过滤单个数据块，则通常假定之前的样本为零。在这种情况下，您可以通过以下方式避免越界索引：

// 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处的“过去样本”以包含当前块的最后一个样本，以便准备处理下一个输入块：

// 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];