傅里叶空间滤波

Question

我有一个真实的向量时间序列 x，长度为 T，还有一个长度为 t 的滤波器 h，t 远小于 T。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. 选择一个小于 t/2 的填充大小 P，选择一个块大小 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，那就太好了。实际上，我的问题是三维的。此外，我不太关心因使用非因果滤波器而产生的伪影。

谢谢，
保罗。

Answer 1

你走在正确的道路上。这个技术叫做重叠保存处理。你提到的t长度够短吗？如果是的话，你可以选择一个块大小B，使得B > 2*min(length(x),length(h))，这样可以加快变换速度。在处理时，你只需要丢掉y_b的前半部分，而不是两头都丢。

为什么要丢掉前半部分呢？记住，频谱相乘在时间域中等同于循环卷积。用零填充的h进行卷积会在结果的前半部分产生奇怪的瞬态噪声，但到后半部分这些瞬态噪声就消失了，因为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的示例，如果你想了解更多细节。