傅立叶空间滤波

问题描述

我有一个长度为T的实矢量时间序列x和一个长度为t≤的滤波器h. 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. 选择填充尺寸P
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，那就太好了.实际的问题是3D.另外，我也不关心因果滤波器的伪影.

谢谢，保罗.

解决方案

您处在正确的轨道上.该技术称为重叠保存处理. 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使用更大的值来延长直线部分的长度.

我写了有关不连续性的问题，并且还在，如果您想了解更多详细信息.

I have a real vector time series x of length T and a filter h of length t << T. h is a filter in fourier space, real and symmetric. It is approximately 1/f.

I would like to filter x with h to get y.

Suppose t == T and FFT's of length T could fit into memory (neither of which are true). To get my filtered x in python, I would do:

import numpy as np
from scipy.signal import fft, ifft

y = np.real( np.ifft( np.fft(x) * h ) ) )

Since the conditions don't hold, I tried the following hack:

1. Select a padding size P < t/2, select a block size B such that B + 2P is a good FFT size
2. Scale h via spline interpolation to be of size B + 2P > t (h_scaled)
3. y = []; Loop:
   • Take block of length B + 2P from x (called x_b)
   • Perform y_b = ifft(fft( x_b ) * h_scaled)
   • Drop padding P from either side of y_b and concatenate with y
   • Select next x_b overlapping with last by P

Is this a good strategy? How do I select the padding P in a good way? What is the proper way to do this? I don't know much signal processing. This is a good chance to learn.

I am using cuFFT to speed things up so it would be great if the bulk of the operations are FFTs. The actual problem is 3D. Also, I am not concerned about artifacts from an acausal filter.

Thanks,Paul.

解决方案

You're on the right track. The technique is called overlap-save processing. Is t short enough that FFTs of that length fit in memory? If so, you can pick your block size B such that B > 2*min(length(x),length(h)) and makes for a fast transform. Then when you process, you drop the first half of y_b, rather than dropping from both ends.

To see why you drop the first half, remember that the spectral multiplication is the same as circular convolution in the time domain. Convolving with the zero-padded h creates weird glitchy transients in the first half of the result, but by the second half all the transients are gone because the circular wrap point in x is aligned with the zero part of h. There's a good explanation of this, with pictures, in "Theory and Application of Digital Signal Processing" by Lawrence Rabiner and Bernard Gold.

It's important that your time domain filter taper to 0 at least on one end so that you don't get ringing artifacts. You mention that h is real in the frequency domain, which implies that it has all 0 phase. Usually, such a signal will be continuous only in a circular fashion, and when used as a filter will create distortion all through the frequency band. One easy way to create a reasonable filter is to design it in the frequency domain with 0 phase, inverse transform, and rotate. For instance:

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

(I'm pretty new to Python, please let me know if there's a better way to code this).

If you compare the frequency responses of ht, htrot, and htwin, you see the following (x-axis is normalized frequency up to pi):

ht, at the top, has lots of ripple. This is due to the discontinuity at the edge. htrot, in the middle, is better, but still has ripple. htwin is nice and smooth, at the expense of flattening out at a slightly higher frequency. Note that you can extend the length of the straight-line section by using a bigger value for N.

I wrote about the discontinuity issue, and also wrote a Matlab/Octave example in another SO question if you want to see more detail.

这篇关于傅立叶空间滤波的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持！