Filtering in the domain of frequency

To make a Fourier-discrete transform to an MxN matrix, MATLAB uses the algorithm Fast Fourier Transform:

F = fft2(f)

where f is the image, and F is the transformed image. The image has to be in double format. You can show the Fourier spectrum with:


S = abs(F);
imshow(S, []);

To show the Fourier transform in more detail, it is better to apply it a logarithmic transform:

figure, imshow(log(1+abs(F)),[])

Fourier transform of an image after a logartimic transform

We need to make a translation to place the centre of the transform in the position (M/2,N/2).

Fc = fftshift(F)

Fourier transform of an image with its center shifted

Just be aware that this:

fftshift(fft2(f))

is not the same as this:

fft2(fftshift(f))

To get the inverse transform:

f = ifft2(F)

In theory, this should give us a real image, but due to issues rounding doubles, it has a small imaginary part. So it's better to use:

f = real(ifft2(F))

If we had already translated the transform, we can invert the process with:

F = ifftshift(Fc)

To make a filtering in the domain of frequency, we use the convolution theorem:

f(x,y)*g(x,y) ⇔ F(u,v)G(u,v)
f(x,y)g(x,y) ⇔ F(u,v)*G(u,v)

This equation shows that we can apply filters in the domain of frequency. But we have again a small problem with borders when we have to apply filters. The discrete-Fourier transform returns a periodic function where each period is the transformed image. This can make strange phenomenons appear when we have to filter the image.

When applying the Fourier transform, we turn the image into a periodic one, and then, do the Fourier transform over that periodic function.

We would have to extend the image borders to avoid problems when doing the transformations. How much? If we suppose that the image f(x,y) has a size of AxB, and the mask to apply h(x,y) has a size of CxD, then:


P ≥ A + C - 1
Q ≥ B + D - 1

To add those borders to our image before we apply a filter, we use the following command, which applies a Fourier transform to our image:

F = fft2(f, P, Q)

Normally, we will work with filters of the same size as the image (MxN), so we use as sizes:


P ≥ 2M - 1
Q ≥ 2N - 1

Steps to follow in order to filter an image:

  1. Obtain the right size to extend an image:
    
    [M N] = size(f);
    P = 2 * M;
    Q = 2 * N;
  2. Obtain the Fourier transform, with extended borders:
    F = fft2(f, P, Q)
  3. Generate a filter function H, with a size of PxQ. The filter should be in the format shown in points 4, 5 and 6. If it is centred, use fftshift(H) to de-centre it.
  4. Multiply the transform by the filter:
    G = H.*F
  5. Obtain the real part of the inverse transform:
    g = real(ifft2(G))
  6. Reduce the image size to its original size:
    g = g(1:m,1:n)

We can create frequency filters from spatial filters:

H = freqz2(h, R, C)

where h is a spatial filter, R is the size in rows we want it to have, C is the size in columns, and H is the filter that works in the domain of frequency. For example, a Sobel filter:

h = fspecial('sobel')';
H = freqz2(h, 500, 500);

Sobel transform

ifftshift(H)

Sobel transform shifted