Fourier space filtering

Question

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.

Answer 1

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.