The Fourier Transform

Jérémie Decock (www.jdhp.org)

Attention:

Veuillez noter que ce document est en cours de rédaction. Dans son état actuel, il n'est pas destiné à être lu par d'autres personnes que ses auteurs.

The following notations come from the following book: Toutes les mathématiques et les bases de l'informatique, Horst Stöcker (Dunod, 2002)

In [2]:
def plot(a, b):
    plt.axis('equal')
    
    plt.axis([-2.*np.pi, 2.*np.pi, -2.*np.pi, 2.*np.pi])
    plt.xticks(np.arange(-2.*np.pi, 2.*np.pi, 1))
    plt.yticks(np.arange(-2.*np.pi, 2.*np.pi, 1))
    
    plt.axhline(y=0, color='k')
    plt.axvline(x=0, color='k')
    
    t = np.arange(-math.pi, math.pi, 0.1)
    y = a * np.cos(t) + b * np.sin(t)
    plt.plot(t, y, "-r",
             t, a*np.cos(t), "--k",
             t, b*np.sin(t), "--k",
             a*np.cos(t), b*np.sin(t), "--b")
    
    plt.grid()
In [3]:
@interact(a=(-5., 5.), b=(-5., 5.))
def cos_sin(a, b):
    plot(a, b)

The Fast Fourier Transform (J.-L. Stark's notation)

The following notations come from the following book: Astronomical Image and Data Analysis, J.-L. Stark and F. Murtagh (Springer, 2002)

One dimension case

Fourier transform of a continuous function f(t) is defined by:

$$ \hat{f}(\nu) = \int_{-\infty}^{+\infty} f(t) e^{-i 2\pi\nu t} dt $$

the inverse Fourier transform is:

$$ f(t) = \int_{-\infty}^{+\infty} \hat{f}(\nu) e^{i 2\pi\nu t} du $$

The discrete Fourier transform is givent by:

$$ \hat{f}(u) = \frac{1}{N} \sum_{k=-\infty}^{+\infty} f(k) e^{-i 2\pi \frac{uk}{N}} $$

and the inverse discrete Fourier transform is:

$$ f(k) = \sum_{u=-\infty}^{+\infty} \hat{f}(u) e^{i 2\pi \frac{uk}{N}} $$

Two dimention case

In case of images (2 variables), the Fourier Transform is:

$$ \hat{f}(u, v) = \frac{1}{MN} \sum_{l=-\infty}^{+\infty} \sum_{k=-\infty}^{+\infty} f(k,l) e^{-i 2\pi \left( \frac{uk}{M} + \frac{vl}{N} \right) } $$

and the inverse discrete Fourier transform is:

$$ f(k,l) = \sum_{u=-\infty}^{+\infty} \sum_{v=-\infty}^{+\infty} \hat{f}(u,v) e^{i 2\pi \left( \frac{uk}{M} + \frac{vl}{N} \right) } $$

Complex notation

Since $\hat{f}(u, v)$ is generally complex, this can be written using the complex its real and imaginary parts:

$$ \hat{f}(u, v) = Re[\hat{f}(u, v)] + i Im[\hat{f}(u, v)] $$

with:

$$ Re[\hat{f}(u, v)] = \frac{1}{MN} \sum_{l=-\infty}^{+\infty} \sum_{k=-\infty}^{+\infty} f(k,l) \cos\left( 2\pi \left( \frac{uk}{M} + \frac{vl}{N} \right) \right) $$$$ Im[\hat{f}(u, v)] = \frac{-1}{MN} \sum_{l=-\infty}^{+\infty} \sum_{k=-\infty}^{+\infty} f(k,l) \sin\left( 2\pi \left( \frac{uk}{M} + \frac{vl}{N} \right) \right) $$

It can also be written using modulus and argument:

$$ \hat{f}(u, v) = |\hat{f}(u, v)| + e^{i ~ \text{arg} ~ \hat{f}(u, v)} $$

$|\hat{f}(u, v)|^2$ is called the power spectrum, and $\Theta(u,v) = \text{arg} ~ \hat{f}(u, v)$ the phase.

Fourier Transform with Python

In [4]:
import math

import numpy as np
import matplotlib.pyplot as plt
from matplotlib import cm

import PIL.Image as pil_img     # PIL.Image is a module not a class...
In [5]:
# %load %load https://raw.githubusercontent.com/jeremiedecock/snippets/master/python/numpy/fft_transform/fft2.py

# Fast Fourier Transform snippet
#
# Additional documentation:
# - Numpy implementation: http://docs.scipy.org/doc/numpy/reference/routines.fft.html
# - Scipy implementation: http://docs.scipy.org/doc/scipy/reference/fftpack.html

# SET OPTIONS ########################################################

shift = True                                      # Shift the zero to the center
threshold = 0.0001                                # The threshold value (between 0 and 1)
file_path = "./sample-images/julie_lebrun.jpeg"   # The file image to filter

# GET DATA ###########################################################

# Open the image and convert it to grayscale
signal = np.array(pil_img.open(file_path).convert('L'))

# Init plot
fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, figsize=(14, 8))

# Plot
ax1.imshow(signal, interpolation='nearest', cmap=cm.gray)
ax1.set_title("Original image")


# FOURIER TRANSFORM WITH NUMPY #######################################

# Do the fourier transform #############

transformed_signal = np.fft.fft2(signal)

if shift:
    transformed_signal = np.fft.fftshift(transformed_signal)

ax2.imshow(np.log10(abs(transformed_signal)),
           interpolation='nearest',
           cmap=cm.gray)
ax2.set_title("Fourier coefficients before filtering")


# Filter ###############################

max_value = np.max(abs(transformed_signal))
filtered_transformed_signal = transformed_signal * (abs(transformed_signal) > max_value*threshold)

ax3.imshow(np.log10(abs(filtered_transformed_signal)),
           interpolation='nearest',
           cmap=cm.gray)
ax3.set_title("Fourier coefficients after filtering")


# Do the reverse transform #############

if shift:
    filtered_transformed_signal = np.fft.ifftshift(filtered_transformed_signal)

filtered_signal = np.fft.ifft2(filtered_transformed_signal)

ax4.imshow(abs(filtered_signal), interpolation='nearest', cmap=cm.gray)
ax4.set_title("Filtered image")
/Users/jdecock/anaconda/lib/python3.5/site-packages/ipykernel_launcher.py:48: RuntimeWarning: divide by zero encountered in log10
Out[5]:
<matplotlib.text.Text at 0x11bcea470>
In [6]:
# %load %load https://raw.githubusercontent.com/jeremiedecock/snippets/master/python/numpy/fft_transform/fft2_basic_test.py

# Fast Fourier Transform snippets
# 
# Documentation:
# - Numpy implementation: http://docs.scipy.org/doc/numpy/reference/routines.fft.html
# - Scipy implementation: http://docs.scipy.org/doc/scipy/reference/fftpack.html

fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, figsize=(14, 8))

# MAKE DATA ##########################################################

pattern = np.zeros((4, 4))
pattern[1:3,1:3] = 1

signal = np.tile(pattern, (2, 2))

ax1.imshow(signal, interpolation='nearest', cmap=cm.gray)
ax1.set_title("Original image")


# FOURIER TRANSFORM WITH NUMPY #######################################
 

# Do the fourier transform #############

transformed_signal = np.fft.fft2(signal)

ax2.imshow(abs(transformed_signal), interpolation='nearest', cmap=cm.gray)
ax2.set_title("Fourier coefficients before filtering")

#shifted_transformed_signal = np.fft.fftshift(transformed_signal)
#shifted_filtered_signal = np.fft.ifftshift(transformed_signal)


# Filter ###############################

max_value = np.max(abs(transformed_signal))
filtered_transformed_signal = transformed_signal * (abs(transformed_signal) > max_value*0.5)

ax3.imshow(abs(filtered_transformed_signal), interpolation='nearest', cmap=cm.gray)
ax3.set_title("Fourier coefficients after filtering")


# Do the reverse transform #############

filtered_signal = np.fft.ifft2(filtered_transformed_signal)

ax4.imshow(abs(filtered_signal), interpolation='nearest', cmap=cm.gray)
ax4.set_title("Filtered image")
Out[6]:
<matplotlib.text.Text at 0x11beaadd8>
In [7]:
# %load %load https://raw.githubusercontent.com/jeremiedecock/snippets/master/python/numpy/fft_transform/fft2_with_noise.py

# Fast Fourier Transform snippet
# 
# Additional documentation:
# - Numpy implementation: http://docs.scipy.org/doc/numpy/reference/routines.fft.html
# - Scipy implementation: http://docs.scipy.org/doc/scipy/reference/fftpack.html

# SET OPTIONS ########################################################

shift = True                                      # Shift the zero to the center
threshold = 0.0001                                # The threshold value (between 0 and 1)
noise_coefficient = 0.5                           # The noise coefficient (between 0 and 1)
file_path = "./sample-images/julie_lebrun.jpeg"   # The file image to filter

# GET DATA ###########################################################

# Open the image and convert it to grayscale
img = np.array(pil_img.open(file_path).convert('L'))

# Add noise
# TODO: lets the user choose the noise method : uniform, normal, poisson, ...
#noise = np.random.rand(*img.shape) * 255. * noise_coefficient
noise = np.random.poisson(255. * noise_coefficient, size=img.shape)

noised_img = img + noise

# Init plots
fig, ((ax1, ax4, ax7, ax10), (ax2, ax5, ax8, ax11), (ax3, ax6, ax9, ax12)) = plt.subplots(3, 4, figsize=(16, 10))

# Display images
ax1.imshow(img, interpolation='nearest', cmap=cm.gray)
ax1.set_title("Original image")

ax2.imshow(noise, interpolation='nearest', cmap=cm.gray)
ax2.set_title("Noise")

ax3.imshow(noised_img, interpolation='nearest', cmap=cm.gray)
ax3.set_title("Noised image")


# FOURIER TRANSFORM WITH NUMPY #######################################

# Do the fourier transform #############

fourier_img = np.fft.fft2(img)
fourier_noise = np.fft.fft2(noise)
fourier_noised_img = np.fft.fft2(noised_img)

if shift:
    fourier_img = np.fft.fftshift(fourier_img)
    fourier_noise = np.fft.fftshift(fourier_noise)
    fourier_noised_img = np.fft.fftshift(fourier_noised_img)

ax4.imshow(np.log10(abs(fourier_img)),
           interpolation='nearest',
           cmap=cm.gray)
ax4.set_title("Fourier coefficients\nbefore filtering (image)")

ax5.imshow(np.log10(abs(fourier_noise)),
           interpolation='nearest',
           cmap=cm.gray)
ax5.set_title("Fourier coefficients\nbefore filtering (noise)")

ax6.imshow(np.log10(abs(fourier_noised_img)),
           interpolation='nearest',
           cmap=cm.gray)
ax6.set_title("Fourier coefficients\nbefore filtering (noised image)")


# Filter ###############################

max_value = np.max(abs(fourier_img))

# TODO: lets the user choose between a threshold based on max, mean, median, ... or a geometric mask (square or circle to the center)
filtered_fourier_img = fourier_img * (abs(fourier_img) > max_value*threshold)
filtered_fourier_noise = fourier_noise * (abs(fourier_noise) > max_value*threshold)
filtered_fourier_noised_img = fourier_noised_img * (abs(fourier_noised_img) > max_value*threshold)

ax7.imshow(np.log10(abs(filtered_fourier_img)),
           interpolation='nearest',
           cmap=cm.gray)
ax7.set_title("Fourier coefficients\nafter filtering (image)")

ax8.imshow(np.log10(abs(filtered_fourier_noise)),
           interpolation='nearest',
           cmap=cm.gray)
ax8.set_title("Fourier coefficients\nafter filtering (noise)")

ax9.imshow(np.log10(abs(filtered_fourier_noised_img)),
           interpolation='nearest',
           cmap=cm.gray)
ax9.set_title("Fourier coefficients\nafter filtering (noised image)")


# Do the reverse transform #############

if shift:
    filtered_fourier_img = np.fft.ifftshift(filtered_fourier_img)
    filtered_fourier_noise = np.fft.ifftshift(filtered_fourier_noise)
    filtered_fourier_noised_img = np.fft.ifftshift(filtered_fourier_noised_img)

filtered_img = np.fft.ifft2(filtered_fourier_img)
filtered_noise = np.fft.ifft2(filtered_fourier_noise)
filtered_noised_img = np.fft.ifft2(filtered_fourier_noised_img)

ax10.imshow(abs(filtered_img), interpolation='nearest', cmap=cm.gray)
ax10.set_title("Filtered image")

ax11.imshow(abs(filtered_noise), interpolation='nearest', cmap=cm.gray)
ax11.set_title("Filtered noise")

ax12.imshow(abs(filtered_noised_img), interpolation='nearest', cmap=cm.gray)
ax12.set_title("Filtered noised image")


# SAVE FILES ######################

ax1.set_axis_off()
ax2.set_axis_off()
ax3.set_axis_off()
ax4.set_axis_off()
ax5.set_axis_off()
ax6.set_axis_off()

ax7.set_axis_off()
ax8.set_axis_off()
ax9.set_axis_off()

ax10.set_axis_off()
ax11.set_axis_off()
ax12.set_axis_off()
/Users/jdecock/anaconda/lib/python3.5/site-packages/ipykernel_launcher.py:80: RuntimeWarning: divide by zero encountered in log10
/Users/jdecock/anaconda/lib/python3.5/site-packages/ipykernel_launcher.py:85: RuntimeWarning: divide by zero encountered in log10
/Users/jdecock/anaconda/lib/python3.5/site-packages/ipykernel_launcher.py:90: RuntimeWarning: divide by zero encountered in log10