Fourier space filtering

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.

5
задан Paul 23 September 2010 в 06:28
поделиться