# 13.4. The Fourier Transform¶

The 2D discrete Fourier Transform decomposes a 2D array into a sum of complex exponentials:

$f(x,y) = \frac{1}{mn} \sum_{u=0}^{m-1} \sum_{v=0}^{n-1} F(u,v) \exp \left[2\pi i \left(\frac{xu}{m} + \frac{yv}{n}\right) \right]$

where the Fourier coefficients are defined by

$F(u,v) = \sum_{x=0}^{m-1} \sum_{y=0}^{n-1} f(x,y) \exp \left[-2\pi i \left(\frac{xu}{m} + \frac{yv}{n}\right) \right]$

They are useful in image processing for performing fast filtering, compression, detection of periodic features, etc.

First let us plot the Fourier coefficients of some simple images:

# Code from previous section
using PyPlot
using Statistics
B = mean(A, dims=3)[:,:,1]
function imshow_scale(A)
# Like imshow(A) but scales the values to [0,1] and supports grayscale

A .-= minimum(A)            # Scale and shift to [0,1]
A ./= maximum(A)
if ndims(A) < 3
A = reshape(A, size(A,1), size(A,2), 1)
end
if size(A,3) == 1
A = repeat(A, 1, 1, 3)  # Set R=G=B for grayscale
end
imshow(A)
end;

# Uncomment below if the package is not already installed

using FFTW

function imagefft_demo(A)
AF = fftshift(fft(A))
subplot(1,2,1); imshow_scale(A);
subplot(1,2,2); imshow_scale(log.(1 .+ abs.(AF)));
return
end

imagefft_demo (generic function with 1 method)

G = [zeros(256,128) ones(256,128)]
imagefft_demo(G) G = zeros(256, 256)
G[78:178, 78:178] .= 1.0
imagefft_demo(G) G = Float32[ (i+j<329) && (i+j>182) && (i-j>-67) && (i-j<73) for i = 1:256, j = 1:256 ]
imagefft_demo(G) G = Float32[ sin(2pi*(10i + 20j)/256) for i = 1:256, j = 1:256 ]
imagefft_demo(G) Many of these patterns can be understood from the underlying Fourier expansion. However, for a general image the pattern of the Fourier spectrum is less clear. We would expect the coefficients from a relatively smooth image to decay away from the center though:

imagefft_demo(B) ## 13.4.1. Removing periodic noise¶

One application of image processing using the Fourier transform is to remove periodic noise. Below we demonstrate this using a made-up example with a given frequency and direction of the noise, but it can be made more general.

Bpernoise = copy(B)
Bpernoise = B + 0.5*Float32[sin(2π*10j / size(B,2)) for i = 1:size(B,1), j = 1:size(B,2) ]
imagefft_demo(Bpernoise) Now compute the Fourier transform, and set the (known) noise frequencies to zero:

# Filter
BF = fftshift(fft(Bpernoise))
mid = size(B) .÷ 2 .+ 1
BF[mid, mid + 10] = 0
BF[mid, mid - 10] = 0
Bfiltered = real.(ifft(ifftshift(BF)))

imagefft_demo(Bfiltered) 