notebook.community

Edit and run

The Fourier Transform





In [ ]:



    


%matplotlib inline

import matplotlib
matplotlib.rcParams['figure.figsize'] = (6, 6)

import math
import cmath          # math functions for complex numbers
import numpy as np
import matplotlib.pyplot as plt

import ipywidgets
from ipywidgets import interact

TODO:

Remplacer http://www.jdhp.org/docs/fourier_transform_intro_fr/fourier_transform_intro_fr.html

  • ajouter notation simplifiée
    • fourier-transform-introduction
    • fourier-transform-slides
    • slides de Bologne
  • ajouter un plot interactif de f(t) (fonction 1D) pour les vecteurs (numpy) de coefs a et b donnés
  • ajouter un plot interactif pour calculer les vecteurs a et b d'une fonction 1D donnée
  • finir notation JL
    • sine transform
    • cosine transform
  • implementations python avec scipy

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

TODO...





In [ ]:



    


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 [ ]:



    


@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

TODO: ça doit être déplacé dans un notebook à part "python_fourier_transform" (?)





In [ ]:



    


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 [ ]:



    


# %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")



    











In [ ]:



    


# %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")



    











In [ ]:



    


# %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()



    











In [ ]:



    


# %load %load https://raw.githubusercontent.com/jeremiedecock/snippets/master/python/numpy/fft_transform/fft_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 ###########################################################

# Make the signal
t = np.arange(5. * -math.pi, 5. * math.pi, 0.01)
#t = np.arange(0., 10. * math.pi, 0.01)
#sig = np.sin(t) + 2. * np.sin(2. * t) + 4. * np.sin(3. * t)
sig = np.cos(t) + 2. * np.cos(2. * t) + 4. * np.cos(3. * t)

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

noised_sig = sig + 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.plot(t, sig)
ax1.set_title("Original signal")

ax2.plot(t, noise)
ax2.set_title("Noise")

ax3.plot(t, noised_sig)
ax3.set_title("Noised signal")


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

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

print(np.fft.rfft(sig))
print(type(np.fft.rfft(sig)[0]))

fourier_sig = np.fft.fft(sig)
fourier_noise = np.fft.fft(noise)
fourier_noised_sig = np.fft.fft(noised_sig)

if shift:
    fourier_sig = np.fft.fftshift(fourier_sig)
    fourier_noise = np.fft.fftshift(fourier_noise)
    fourier_noised_sig = np.fft.fftshift(fourier_noised_sig)

ax4.plot(t, abs(fourier_sig))
ax4.set_title("Fourier coefficients\nbefore filtering (signal)")

ax5.plot(t, abs(fourier_noise))
ax5.set_title("Fourier coefficients\nbefore filtering (noise)")

ax6.plot(t, abs(fourier_noised_sig))
ax6.set_title("Fourier coefficients\nbefore filtering (noised signal)")


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

max_value = np.max(abs(fourier_sig))

# 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_sig = fourier_sig * (abs(fourier_sig) > max_value*threshold)
filtered_fourier_noise = fourier_noise * (abs(fourier_noise) > max_value*threshold)
filtered_fourier_noised_sig = fourier_noised_sig * (abs(fourier_noised_sig) > max_value*threshold)

ax7.plot(t, abs(filtered_fourier_sig))
ax7.set_title("Fourier coefficients\nafter filtering (signal)")

ax8.plot(t, abs(filtered_fourier_noise))
ax8.set_title("Fourier coefficients\nafter filtering (noise)")

ax9.plot(t, abs(filtered_fourier_noised_sig))
ax9.set_title("Fourier coefficients\nafter filtering (noised signal)")


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

if shift:
    filtered_fourier_sig = np.fft.ifftshift(filtered_fourier_sig)
    filtered_fourier_noise = np.fft.ifftshift(filtered_fourier_noise)
    filtered_fourier_noised_sig = np.fft.ifftshift(filtered_fourier_noised_sig)

filtered_sig = np.fft.ifft(filtered_fourier_sig)
filtered_noise = np.fft.ifft(filtered_fourier_noise)
filtered_noised_sig = np.fft.ifft(filtered_fourier_noised_sig)

ax10.plot(t, abs(filtered_sig))
ax10.set_title("Filtered signal")

ax11.plot(t, abs(filtered_noise))
ax11.set_title("Filtered noise")

ax12.plot(t, abs(filtered_noised_sig))
ax12.set_title("Filtered noised signal")

# 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()



    











In [ ]:



    


# %load %load https://raw.githubusercontent.com/jeremiedecock/snippets/master/python/numpy/fft_transform/plot.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 ###########################################################

# Make the signal
t = np.arange(5. * -math.pi, 5. * math.pi, 0.01)
#t = np.arange(0., 10. * math.pi, 0.01)
#sig = np.sin(t) + 2. * np.sin(2. * t) + 4. * np.sin(3. * t)
sig = np.cos(t) + 2. * np.cos(2. * t) + 4. * np.cos(3. * t)

sig1 = np.cos(t)
sig2 = 2. * np.cos(2. * t)
sig3 = 4. * np.cos(3. * t)

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

noised_sig = sig + noise

## Init plots
#fig, (ax1, ax2, ax3, ax4) = plt.subplots(1, 4, figsize=(16, 5))
#
#N = 8
#
## Display images
#ax1.plot(t, sig)
#ax1.set_title(r"$f(t)$ = ")
#ax1.set_ylim([-N, N])
#
#ax2.plot(t, sig1)
#ax2.set_title(r"$a_1 \times \cos(t)$ +")
#ax2.set_ylim([-N, N])
#
#ax3.plot(t, sig2)
#ax3.set_title(r"$a_2 \times \cos(2t)$ +")
#ax3.set_ylim([-N, N])
#
#ax4.plot(t, sig3)
#ax4.set_title(r"$a_3 \times \cos(3t)$")
#ax4.set_ylim([-N, N])

# Init plots
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(16, 5))

# Display images
ax1.plot(t, sig)
ax1.set_title(r"Direct space")
#ax1.set_ylim([-N, N])

ax2.bar([1, 2, 3], [1, 2, 4], 0.01)
ax2.set_title(r"Fourier space")
ax2.set_xlim([0, 4])
ax2.set_ylim([0, 5])



    











In [ ]:



    


# %load %load https://raw.githubusercontent.com/jeremiedecock/snippets/master/python/numpy/fft_transform/rfft2.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.rfft2(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.irfft2(filtered_transformed_signal)

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

Content source: jdhp-docs/python-notebooks

Similar notebooks:

notebook.community | gallery | about