傅里叶空间滤波

Posted 2023-03-12

【中文标题】傅里叶空间滤波

【英文标题】：Fourier space filtering

【发布时间】：2011-04-16 02:39:19

【问题描述】：

我有一个长度为 T 的实向量时间序列 x 和一个长度为 t

我想用 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 连接 选择下一个 x_b 与最后一个重叠的 P

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

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

谢谢， 保罗。

【参考方案1】：

你在正确的轨道上。该技术称为overlap-save processing。 t 是否足够短以至于该长度的 FFT 适合内存？如果是这样，你可以选择你的块大小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 示例。

【讨论】：

感谢重叠保存参考。我在 Press et al., Numerical Recipes 中读到了关于时域滤波的内容，但我不确定如何将其映射到频域。我不确定是否要放弃：1）为什么 y_b 的后半部分而不是结尾，2）在您的其他 SO 帖子中，您放弃了前半部分。

我的过滤器 h 是从原始数据的平均值得出的，h(f) ~ 1/f 和相位设置为 0。我正在用这个过滤器过滤合成信号以提供更多的频谱就像我的原始数据一样。我不确定如何在时域中设计此过滤器。您指出的一件事是 ifft(h) 最好在一端为零以避免振铃伪影。我会检查这个，因为它很可能不是。是否可以将时域中的汉明窗应用于频域中的某个窗口方法（您在其他 SO 帖子中的第一个示例）？

是的 - 很抱歉搞砸了第一半/第二半的问题。我更新了该更正，以及一些关于生成行为良好的 h 的想法。

您是否有机会尝试使用1/f 频率分布创建随机噪声？有一些技术可以直接生成此类信号 - 例如，请参阅 local.wasp.uwa.edu.au/~pbourke/fractals/noise。

非常感谢代码和插图。这很有帮助。我将不得不阅读使用你的提示。我借了 Proakis 的书。关于 1/f 分布的生成，我正在做与您链接中的示例“景观频率合成”类似的事情作为测试用例，除了我正在制作电影。确定性生成器是最有趣的。