FFT Low Pass Filtering methods

Biftheunderstudy

Senior member
Aug 15, 2006
375
1
81
First off, let me announce that I know very little about DSP, and I'm coming at this from a mathematical perspective so my language may be....different....

I have a set of data stored in a 3D array of size 64x128x128 which I would like to decompose into "large scale" and "small scale" components. Practically, this means I would like to do a low and high pass filter respectively. The cut-off frequency I'm interested in, corresponds to about the first 10 Fourier components (whatever that works out to).

My first attempt involved a simple brick wall filter, I fft'ed the data, then set everything but the first 10 (and last 10) components along each axis to zero and called this a low-pass. To get the high-pass, I took the original dataset and subtracted the low-pass result from it.

After some research, I came to the realization that this was introducing sinc ringing in my time domain signal and I should do some sort of windowing.

So, instead of the brick wall filter, I made a Tukey window (like a tunable Hanning window) and applied it in the same way as the brick wall filter. This seems to give better results on a test array made of noise (it looks like a lowpass filter).

My question is, can I do better and are there some "gotcha" type scenarios that I need to watch out for? I'm programming this all in python (numpy, scipy, matplotlib etc).

I'm investigating the use of time domain filtering, but I can't wrap my head around the DSP language and the nuances involved. One of the things which may complicate everything is that the data is periodic (produced from a simulation with periodic boundary conditions).

Notably, using numpy.signal and firwin along with lfilter is there a way to take advantage of the periodicity to prime the filter, and how would I do it?

P.S.
I haven't included any code, since it's a mess and I kind of want a fresh approach.
 

esun

Platinum Member
Nov 12, 2001
2,214
0
0
I think you have the idea for frequency-domain filtering. Window, take the FFT, zero components, and inverse FFT. Choosing the window function is pretty subjective based on what you want. You can play around with different ones until you get a result you like.

One thing you could do since your signal is periodic is create a periodic repetition of it before processing it. This will make the resolution of the FFT better. Of course, this assumes that the signal is actually supposed to be periodic and that its periodicity is not just a simulation artifact (otherwise, you don't want to use the periodic extension since that will distort the result).

As for time-domain filtering, you'll have to play around with some filter design tools probably. I recall scipy has some filter functions that you can basically specify passband and stopband parameters to and it will construct a filter for you. If you are having trouble determining the passband, you can convert the 10th FFT bin into a relative frequency and use that value.

As for taking advantage of the periodicity for time-domain filtering, you can again repeat your signal so that the filter processes a periodic representation before getting to your actual signal. This may be more desirable than a zero boundary condition. You'd only need to buffer up to the delay of the filter (so if you have a 100-tap filter, then you'd need 100 points leading up to your starting point). You of course have to remember not to use the outputs prior to the start of your "real" data (i.e., up to the number of taps in your filter).
 
Last edited:

slugg

Diamond Member
Feb 17, 2002
4,723
80
91
I came in here thinking this was an image processing problem, AKA a 2D spacial domain. Nope. Now I'm reminded that I gotta hit the books and learn something about DSP.

Not to derail, but how does the Sci-Py stack fare against Matlab?
 

esun

Platinum Member
Nov 12, 2001
2,214
0
0
I much prefer SciPy (and related tools) over Matlab. I gave up Matlab at my last job (where Matlab was pretty much entrenched) and went to Python cold turkey and had little trouble doing everything I needed to do in Python.

The only way I'd use Matlab in the future is if I absolutely needed Simulink or if I needed a large subset of a specific toolbox that Python doesn't currently have an equivalent of and didn't have time to roll my own (since, honestly, if you spend a couple weeks making your own it'll save you in the long run vs a $4000 license + recurring costs + lock-in to Matlab).

One thing I must recommend is installing Numpy binaries compiled against the MKL. You can get them for Windows here:

http://www.lfd.uci.edu/~gohlke/pythonlibs/#numpy

There are some instructions here for compiling with the Intel compiler and MKL for Linux:

http://www.scipy.org/Installing_SciPy/Linux
 

slugg

Diamond Member
Feb 17, 2002
4,723
80
91
Glad to hear that the changeover is smooth. I'll give it a shot. I've learned a lot from your posts in the past, esun, so I hope you're right again! :)
 

Biftheunderstudy

Senior member
Aug 15, 2006
375
1
81
After a bit of fiddling with scipy.signal, I found the filtfilt and butter functions and did a lowpass filter and compared to my method. The time-domain filtering appears to be pretty bad at the edges (probably because I'm a noob though).

In the end, I convinced myself that the fft method was the better choice for my case.

Now I'm curious about your comment of adding copies of the data on the edges...I've heard of this, but in my case it would make the fft rediculous. (I lied a little, the data cube is 128x512x512 which already pushes my 16GB of RAM to its limit). The data should be perfectly periodic, so the fft is well suited to this sort of thing.

So this got me thinking, what if I copied the data and then used an overlap-add method? Although this is probably more trouble than its worth...


I also switched from Matlab cold turkey. Python is free, and fixes some annoying bugs with Matlab. The mutable types and making basically everything a pointer is very frustrating though.

I also switched from writing it in Fortran90 -- in a funny way, aside from the row-major vs column major and indexing from 0, the two are very similar.