傅立叶空间滤波

Question

我有一个长度为T的实矢量时间序列x和一个长度为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,选择块大小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连接
   • 选择下一个x_b与最后一个重叠P

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

我正在使用cuFFT来加快速度,因此如果大部分操作都是FFT,那将会很棒.实际问题是3D.此外,我不关心非因果过滤器的工件.

谢谢,保罗.

Answer 1

你走在正确的轨道上.该技术称为重叠保存处理.是t足够短,这是长度的FFT存放在内存？如果是这样,您可以选择块大小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例子.