零填充奇数长度 FFT 时的实值振铃

信息处理 fft 插值 C
2021-12-29 07:14:44

所以我正在尝试编写一个频域插值器,它对信号的频率响应和逆变换进行零填充。我必须处理两种情况:

  1. 偶数响应 - 必须拆分Fs/2bin 因为它是模棱两可的。所以我复制了频谱的负部分,并n*(interp-1)-1在两者之间添加了零。
  2. 奇数响应 - 没有Fs/2bin 所以只需拆分正/负频率并n*(interp-1)在它们之间插入零。

可以在这里看到执行零填充的代码

// Copy negative frequency components to end of buffer and zero out middle
//  inp    - input buffer of complex floats
//    n    - transform size
//  interp - interpolation amount
void zero_pad_freq(cfloat_t *inp, size_t n, size_t interp) {
    if ((n % 2) == 0) {
        memmove(inp + n*interp - n/2, inp + n/2,     n/2*sizeof(cfloat_t));
        memset (inp + n/2 + 1, 0,       (n*(interp-1)-1)*sizeof(cfloat_t)); // Duplicate Fs/2 so we need one less zero

        inp[n/2]          /= 2.0;
        inp[n*interp-n/2] /= 2.0;
    } else {
        memmove(inp + n*interp - n/2, inp + (n+1)/2, n/2*sizeof(cfloat_t));
        memset (inp + (n+1)/2, 0,         (n*(interp-1))*sizeof(cfloat_t));
    }
}

第一种情况工作正常,我在啁啾信号上测试它,它插值很好,有一点数字噪声,但它通过 FFT 往返,所以你能做什么(首先50μs信号显示左右):

问题在于奇数长度变换,我只在真实样本上得到了非常令人发指的瞬态响应(50μs再次,真实的):

想象中的通道上有一个小涟漪,但没有那么糟糕:

这就像我搞砸了我的Fs/2bin 在奇怪的情况下,但没有Fs/2斌,所以我很不解。有人有什么想法吗?

2个回答

通过将高频区间归零,您可以有效地将信号频谱与矩形函数相乘。频率乘法是时间卷积,矩形的傅立叶对是 sinc。所以你真正做的是将时域信号与正弦波卷积,正弦波的主瓣宽度与矩形的长度成反比。这就是为什么众多滤波器设计技术(如Parks-McClellan)在所谓的“过渡区域”或“过渡”频带中设计,以使滤波器的频率响应不会发生瞬时变化。这些滤波器设计技术很重要,因为像您使用的“理想”滤波器在时域中具有如此不良的影响。

频域中的一个步骤将显示为时域中的波纹。如果您使用窗函数(例如汉明窗)平滑频率数据,它应该会显着减少波纹。