python - 傅里叶空间滤波

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

由于条件不成立，我尝试了以下 hack:

1. 选择填充大小 P < t/2，选择 block 大小 B，使得 B + 2P 是一个很好的 FFT 大小
2. 通过样条插值将 h 缩放为大小 B + 2P > t (h_scaled)
3. y = [];环形:
   • 从x中取出长度为B+2P的 block (称为x_b)
   • 执行 y_b = ifft(fft( x_b ) * h_scaled)
   • 从 y_b 的任一侧删除填充 P 并与 y 连接
   • 通过P选择下一个与上一个重叠的x_b

这是一个好的策略吗？怎样选择padding P比较好？这样做的正确方法是什么？我不太了解信号处理。这是学习的好机会。

我正在使用 cuFFT 来加快速度，所以如果大部分操作都是 FFT，那就太好了。实际问题是 3D。此外，我不担心来自非因果过滤器的伪影。

谢谢， 保罗。

Answer 1

您走在正确的轨道上。该技术称为 overlap-save processing . t 是否足够短，以至于该长度的 FFT 适合内存？如果是这样，您可以选择 block 大小 B 使 B > 2*min(length(x),length(h)) 并进行快速转换。然后，当您处理时，您将丢弃 y_b 的前半部分，而不是从两端丢弃。

要了解为什么要删除前半部分，请记住频谱乘法与时域中的循环卷积相同。与零填充的 h 卷积会在结果的前半部分产生奇怪的瞬变，但到后半部分所有瞬变都消失了，因为 x 中的循环环绕点与 h 的零部分对齐。在 "Theory and Application of Digital Signal Processing" by Lawrence Rabiner and 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 值来延长直线段的长度。

我写了关于不连续性问题的文章，还在 another SO question 中写了一个 Matlab/Octave 示例如果您想查看更多详细信息。