## Friday, August 5, 2011

### Fourier Transforms and image filtering

In the previous post we have seen how to compute the fft of a one dimensional signal. In this post we will see how to compute and plot the Fourier spectrum of an image (we can think at images as a bidimensional signals) and how to use the spectrum in order to make an high pass filter in the frequency domain. The following function implements an high pass filter in three steps:
• computes the frequency spectrum
• cuts the low frequencies from the spectrum
• and finally, the filtered time domain signal is calculated from the modified frequency spectrum
```from pylab import *
from numpy import *

def freqHighPass(img, r):
"""
Applies an High Pass filter to img in the frequency domain.
In:  img, image to filter
r, percentage of spectrum to keep (in [0 1])
Out: ifiltered, filtered version of the image
P, filtered spectrum
"""
I = fftshift(fft2(img)) # entering to frequency domain
# fftshift moves zero-frequency component
# to the center of the array
P = zeros(I.shape,dtype=complex)

r = int((r*min(img.shape))/2);
c1 = I.shape[0]/2 # spectrum center
c2 = I.shape[1]/2

for i in range(c1-r,c1+r):  # frequency cutting
for j in range(c2-r,c2+r): # around the center
P[i,j] = I[i,j]

ifiltered = real(ifft2(ifftshift(P)))
return ifiltered,P # back to the spatial domain
```
This is the first application of the function freqHighPass():
```i = imread('diane.jpg') # load an image
i = mean(i,2) # to get a 2-D array

ifiltered,P = freqHighPass(i,0.4) # 40% of spectrum preserved

subplot(2,2,1)
title('Original - Fourier Specturm')
imshow(log(abs(fftshift(fft2(i))))**2)
subplot(2,2,2)
title('Filtered - Fourier Specturm')
imshow(log(abs(P))**2)
subplot(2,2,3)
title('Original image')
imshow(flipud(i))
gray()
subplot(2,2,4)
title('Filtered image')
imshow(flipud(ifiltered))
show()
```
The result of the code above is the following plot. We can see the spectrum of the original image, the filtered spectrum (the matrix P returned by freqHighPass), the original image and the filtered image.

And this is another application that preserve more spectrum than the first:
```i = imread('castle.png')

ifiltered,P = freqHighPass(i,0.7) # 70% of spectrum preserved
figure()
subplot(1,2,1)
title('Original image')
imshow(i)
gray()
subplot(1,2,2)
title('Filtered image')
imshow(ifiltered)
show()
```
In this case the original image is affected by noise. As we can see the filter is able to reduce noise.

1. Can you compare this to the DCT and how to implement that?

2. Yes, maybe in the next post. DCT rule but the my favorite are always Wavelets :)

3. That's cool. Each method is ideal to some situation. My favorite is PCA!

4. Thx for great post. But how are you avoiding ringing in the output image? Wouldn't the sharp cut-off in I cause ringing? (I am a new student to signal processing.)

5. Hi testcluster, I don't care about the ringing effect in the post. If you need to avoid the ringing effect you have to choose the cut-off frequency equal to (or very close to) Nyquist frequency.
You can find out more about the Nyquist frequency in this post:

http://glowingpython.blogspot.com/2011/09/sampling-theorem-explained-with-numpy.html

6. Isn't that a low pass filter?

1. It is a high pass filter because it passes high frequencies signals but blocks frequencies lower than the cutoff frequency. Just look at the filtered Fourier specturm, it as only the highest frequencies of the original one.