傅里叶域滤波。

Question

傅里叶域滤波。

5

我有一个长度为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：

选择填充大小P
通过样条插值将h缩放到大小B + 2P>t（h_scaled）
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的。此外，我不关心来自非因果滤波器的伪影。

谢谢，保罗。

- Paul

1个回答

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- mtrw · Accepted Answer

您正在正确的道路上。这种技术称为重叠保存处理。如果t足够短，以至于该长度的FFT适合内存，那么您可以选择块大小B，使得B > 2 * min(length(x),length(h))，从而实现快速变换。然后在处理时，您会丢弃y_b的前半部分，而不是从两端丢弃。

要了解为什么要丢弃前半部分，请记住，频谱乘法与时间域中的循环卷积相同。用零填充的h进行卷积会在结果的前半部分创建奇怪的故障传输，但到第二半部分时，所有传输都消失了，因为x中的循环包装点与h的零部分对齐。在"数字信号处理的理论和应用" by Lawrence Rabiner and Bernard Gold中有一个好的解释，附有图片。

重要的是，你的时域滤波器至少在一端渐变为0，以避免出现振铃伪影。你提到h在频域中是实数，这意味着它具有全零相位。通常，这样的信号只会以循环方式连续存在，并且在用作滤波器时会在整个频带中产生失真。创建一个合理的滤波器的一种简单方法是在频域中设计一个具有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）： 1/f filters with length 64

ht在顶部有很多纹波。这是由于边缘处的不连续性引起的。htrot在中间更好，但仍然有纹波。htwin非常平滑，但代价是在略高的频率处变平。请注意，您可以通过使用更大的N值来延长直线段的长度。

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