How to Actually Use FFT to Implement Low Pass Filter

Question

have searched high and low for a practical guide to implement a simple 50Hz low pass filter on a signal.I understand the theory but its the practical but that's missing. If I was using MATLAB it would be trivial because there are so many examples but I am coding this by code.

This is where I am at:

1: I have a series of 8192 bytes sampled at 1000Hz in an array (0-256)

2: I have generated the frequency domain representation of a simple step filter in a symmetric array (+/-F, Ire, Iim) :

0,1,0
0.1220703125,1,0
0.244140625,1,0
0.3662109375,1,0
0.48828125,1,0
0.6103515625,1,0
.....
-0.6103515625,1,0
-0.48828125,1,0
-0.3662109375,1,0
-0.244140625,1,0
-0.1220703125,1,0

3: I have used Inverse FFT of the filter to obtain the time representations of the filter as (t, Ire, Iim) [as expected the im parts are ~0 and the filter gives the expected decaying ringing]:

0.000,819,0
0.001,805.600848087199,4.11892742135933E-12
0.002,766.191083986647,7.80009390410896E-12
0.003,703.07869137741,1.07704956064936E-11
0.004,619.930866652693,1.26768595620774E-11
0.005,521.519198753447,1.32877042702262E-11
0.006,413.390428478234,1.26576527037514E-11
0.007,301.48763797423,1.08400510789863E-11
8,0.008,191.750273447113,7.79765141345479E-12
....
8.186,413.390428478235,7.79385628257856E-09
8.187,521.519198753448,9.84392145575441E-09
8.188,619.930866652693,1.17116023545805E-08
8.189,703.07869137741,1.32913150485692E-08
8.19,766.191083986645,1.44921696865197E-08
8.191,805.600848087196,1.52439926792702E-08

OK so what do I do with this data???

From what I have read in the literature the thoery is that in the freq domain you multiply the spectra but because this is the DFT the usual course is to convolve in the time domain but the question remains:

1: Do I simply multiply the two time series together (doesn't feel right)?

2: If not how is the weighting applied?

3: How do I normalise any output back to 0-256 values?

I realise this is a pretty basic question but there doesn't seem to be a really basic answer or example out there :-(

Thanks

iced

Answer 1

this is not really a question but is an answer to a question I asked here and the answer is a little long winded so I thought I would post the answer here as a guide to others. I have done this because I found it demonically difficult to find a real "how to" for do basic filtering. I found a few guides to this helpful here and here So here goes, I hope it helps some others:

I had a data set of biometric data that was full of hum and noise, most of which was >30Hz and I had tried various "averaging" techniques to get rid of this with varying success so I turned to FFT to do the job. While I had studied transforms in general at uni and understood the basic theory and maths I needed to code this for a small app for a tri club (no its not commercial!).

So here is the psuedo code with comments:

SF: Sample Frequency
Samples: is array with samples 0..255
N: is number of samples (8192)
NP: is root of 2^NP = N (13)
CF: cutoff freq of filter (30)
W: number of coefficients in filter (40)
dare: real data input
daim: imaginary data input

for i = 0 To N – 1
{
// time axis for graphing results
time(i) = i / SF

// freq axis for graphing, note symetric FFT so need -ve freq
if i <= N / 2 then // Positive frequencies
freq(i) = i * SF / N
else // Negative frequencies
freq(i) = i * SF / N - SF

// copy input data to local array 
// and normalise to 127 = 0, centering the input/transform

dare(i) = Samples(i) – 127
daim(i) = 0
}

// fft the input data
FFT(NP, N, dare, daim)

// Create the impulse response of the filter in the frequency domain
// this could be any filter but I use a simple box here: 1 in the passband and 0 elsewhere. Could use any filter here:

// Note that I used symmetric complex FFT: it is necessary
// fill in the negative frequencies as well.

for i = 0 To N - 1
{
re(i) = 0
if i <= (CF * N / SF) or i >= (N - CF * N / SF) 
re(i) = 1
im(i) = 0
}

// Inverse filter into the time domain 
// Note: imaginary parts should be ~ 0 after iFFT
IFFT(NP, N, re, im)

// Hann window to truncate the filter series at say +/-20
// can use other window formulas (e.g. bartlett) if desired
for i = 0 To N – 1
{
window(i) = 0
if i <= Won2 then
window(i) = (0.5 + 0.5 * Cos(i * 2pi / (W + 1)))
else if i >= N - Won2 then
window(i) = (0.5 + 0.5 * Cos((N - i) * 2pi / (W + 1)))
}    

// Apply the Hann window to the filter
for i = 0 To N – 1
{
//careful here if complex multiply required
re(i) = re(i) * window(i)
im(i) = 0
}

// can IFFT and check for freq response and back again

// Shift the filter to make it i.e. +ve times only
//shift lower coefs up
for i = Won2 To 0 Step -1
{
re(i + Won2) = re(i) // 
}
//wrap -ve coeffs around
for i = 0 To Won2 – 1
{
re(i) = re(N - Won2 + i) //wrap N -> 0, N - 1 -> 1 etc
re(N - Won2 + i) = 0
}

// Transform filter back into the frequency domain 
FFT(NP, N, re, im)

//multiply input data by filter COMPLEX MULT REQUIRED!!
// may be able to shorten this step if desired
for i = 0 To N – 1
{
tmpre = dare(i) * re(i) - daim(i) * im(i)
tmpim = dare(i) * im(i) + daim(i) * re(i)
dare(i) = tmpre
daim(i) = tmpim
}

//inverse back to time domain
IFFT(NP, N, dare, daim)

// copy back to input normalised and re centered
// may need to normalise/scale at this stage
for i = 0 To N – 1
{
Samples(i) = dare(i) + 127
}

Answer 2

1. You convolve the two time series together. Think of the filter as a kind of weighted moving average. Imagine how you would implement a moving average, as it is also a convolution operation. E.g. a five sample moving average is convolution with a filter with response [0.2 0.2 0.2 0.2 0.2]

Or you can multiply the Fourier transform of the signal and the Fourier transform of the filter together, then perform the inverse Fourier transform. It is (almost) the same as convolution, with the exception that the Fourier transform assumes the signal is periodic, therefore filtering at the start of the signal used values from the end as well, and vice versa. Convolution can do this as well, but usually it is implemented by padding the signal with zeros at the start and end.

3. $$f'(x) = \frac{f(x) - \min(f)}{\max(f) - \min(f)} \times 255$$

(bytes are 0 - 255)