傅立叶空间滤波

Question

我有一个长度为T的实向量时间序列x和一个长度为t << T的滤波器h。h是傅立叶空间中的一个实且对称的滤波器。大约是1 / f。

我想用h过滤x以得到y。

假设t == T，长度为T的FFT可以放入存储器中（都不成立）。要在python中获取经过过滤的x，我可以这样做：

import numpy as np
from scipy.signal import fft, ifft

y = np.real( np.ifft( np.fft(x) * h ) ) )

由于条件不成立，我尝试了以下技巧：

1. 选择一个填充大小P <t / 2，选择一个块大小B，使B + 2P是一个好的FFT大小
2. 通过样条插值将h缩放为B + 2P> t（h_scaled）
3. y = []; 循环：
   • 从x取长度为B + 2P的块（称为x_b）
   • 执行y_b = ifft（fft（x_b）* h_scaled）
   • 从y_b的任一侧放下填充P，并与y连接
   • 选择下一个与P重叠的下一个x_b

这是个好策略吗？如何选择一个好的填充P？正确的方法是什么？我不太了解信号处理。这是学习的好机会。

我正在使用cuFFT来加快速度，因此，如果大部分操作是FFT，那将是很好的。实际的问题是3D。另外，我也不关心因果滤波器的伪影。

谢谢，保罗。

Answer 1

您走在正确的轨道上。该技术称为重叠保存处理。是否t足够短，以至于该长度的FFT可以存储在内存中？如果是的话，你可以选择你的块大小B，从而`B

2*min(length(x),length(h))为快速变换，使。然后，当您处理时，您将放下的前一半y_b`，而不是从两端都放下。

要了解为什么要丢弃前半部分，请记住，频谱乘法与时域中的圆形卷积相同。与填充零的卷积h会在结果的上半部分产生怪异的毛刺瞬变，但是到下半部分，所有瞬变都消失了，因为in中的圆形环绕点x与的零部分对齐h。劳伦斯·拉宾纳（Lawrence
Rabiner）和伯纳德·戈德（Bernard
Gold）

在“数字信号处理的理论和应用”中对图片进行了很好的解释。

重要的是，您的时域过滤器至少一端会逐渐减小到0，以免出现振铃失真。您提到h在频域中这是真实的，这意味着它具有全0相位。通常，这样的信号将仅以循环的方式连续，并且当用作滤波器时会在整个频带上产生失真。创建合理滤波器的一种简单方法是在相位为0的频域中对其进行设计，逆变换和旋转。例如：

def OneOverF(N):
    import numpy as np
    N2 = N/2; #N has to be even!
    x = np.hstack((1, np.arange(1, N2+1), np.arange(N2-1, 0, -1)))
    hf = 1/(2*np.pi*x/N2)
    ht = np.real(np.fft.ifft(hf)) # discard tiny imag part from numerical error
    htrot = np.roll(ht, N2)
    htwin = htrot * np.hamming(N)
    return ht, htrot, htwin

（我对Python还是很陌生，请告诉我是否有更好的编码方法）。

如果你比较的频率响应ht，htrot和htwin，您将看到以下（X轴是归一化频率高达pi）：

ht顶部有很多涟漪。这是由于边缘不连续。 htrot，在中间，效果更好，但仍有涟漪。
htwin平滑流畅，但以稍高的频率展平为代价。请注意，可以通过为N使用更大的值来延长直线部分的长度。

我写了有关不连续性问题的文章，如果想了解更多细节，还写了另一个SO问题中的Matlab / Octave示例。