带有或不带有窗口的KISS FFT输出

Jin*_*Jin 4 signal-processing windowing dft kissfft

我目前正在尝试使用Kiss FFT将FFT实施到AVR32微控制器中,以进行信号处理。而且我的输出有一个奇怪的问题。基本上,我会将ADC样本(使用函数发生器进行测试)传递到fft(真实输入,256 n大小)中,并且检索到的输出对我来说很有意义。但是,如果我将汉明窗应用于ADC样本,然后将其传递给FFT,则峰值幅度的频率仓是错误的(并且与之前没有开窗的结果不同)。ADC样本具有DC偏移,因此我消除了偏移,但仍不适用于窗口样本。

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

    488260 //直流仓
    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);        
Run Code Online (Sandbox Code Playgroud)

Sle*_*Eye 5

频域输出说明

在频域中,矩形和汉明窗如下所示:

频域中的矩形和汉明窗

如您所知,时域乘以一个窗口对应于频域中的卷积,这实际上将信号的能量散布在多个频点上,这通常被称为频谱泄漏。对于您选择的特定窗口(上面在频域中进行了描述),您可能会注意到汉明窗在主瓣外散布的能量要少得多,但该主瓣的宽度比矩形窗宽一些。

结果,使用汉明窗时,DC能量的尖峰最终会扩散到bin 0上方,并扩散到bin 1中。并不是说在bin 1中有一个很强的峰值。实际上,如果对提供的数据进行绘图,您应该看到在索引6处看到的峰值实际上仍然与和处于相同的频率。不使用汉明窗:

数据

如果要消除周围信号仓中的DC尖峰和泄漏,请消除数据中的偏差(本质上应用陷波滤波器),或者在寻找“第一个强峰值”。

过滤问题

最后,需要注意的是还有与IIR滤波器来实现,由此阵列的方式的问题x,并y会界被索引出来的时候ctr==0ctr==1(除非你已经做了一些特殊规定,并xy实际上指针从一开始的偏移分配的缓冲区数)。无论有无窗口,这都可能影响结果。如果仅过滤单个数据块,则通常的假设是较早的样本为零。在这种情况下,您可以使用以下方法避免超出范围的索引编制:

// 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]);
}
Run Code Online (Sandbox Code Playgroud)

另一方面,如果您要过滤多个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;
Run Code Online (Sandbox Code Playgroud)

然后,您可以在这些缓冲区中使用偏移量。索引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];
Run Code Online (Sandbox Code Playgroud)

最后,在每个块的末尾,您必须在索引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];
Run Code Online (Sandbox Code Playgroud)