2D Fourier Transforms

Preliminaries

import numpy as np
from numpy.fft import fft, ifft, fft2, ifft2, fftshift, ifftshift
import matplotlib.pyplot as plt
from scipy import ndimage, misc
from scipy.signal import gaussian

#from IPython.display import display, HTML
#display(HTML("<style>.container { width:98% !important; }</style>"))

np.set_printoptions(edgeitems=30, linewidth=250, formatter=dict(float=lambda x: "% .4f" % x))
np.set_printoptions(edgeitems=30, linewidth=250, formatter=dict(complex=lambda x: "% .4f" % x))

fftshift

# The DC (0) is at index 0. This is the order of the frequency indices
# for the Fourier coefficients that come out of the DFT.
k = list(range(16))
k

[0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15]

# The function fftshift moves (rotates) the Fourier coefs so that
# the DC is at (or beside) the centre, at index 8 in this case.
s = fftshift(k)
print(s)

[ 8  9 10 11 12 13 14 15  0  1  2  3  4  5  6  7]

f = np.arange(1,8, dtype=float)
print(f)

print('Raw DFT output')
F = np.fft.fft(f)
print(np.round(F, 2))

print('Shifted DFT output')
Fs = fftshift(F)
print(np.round(Fs, 2))

[1. 2. 3. 4. 5. 6. 7.]
Raw DFT output
[28. +0.j   -3.5+7.27j -3.5+2.79j -3.5+0.8j  -3.5-0.8j  -3.5-2.79j -3.5-7.27j]
Shifted DFT output
[-3.5-0.8j  -3.5-2.79j -3.5-7.27j 28. +0.j   -3.5+7.27j -3.5+2.79j -3.5+0.8j ]

Images

f = np.array(plt.imread('images/pd.jpg')[:,:,0], dtype=float)
plt.imshow(f,cmap='gray');

View the Fourier coefficients

F = fft2(f)

# Plot the modulus of the Fourier coefs
# DC is in the top-left corner
plt.figure(figsize=[6,6]);
plt.imshow(abs(F), cmap='gray');
plt.title(r"Modulus of the Fourier coefficients, $|F|$");

print(np.round(F[:6, :3], decimals=2))

[[ 2491737.       +0.j   -1504009.57 +25372.13j    39770.1    +859.05j]
 [ -972029.69 -40667.9j    471657.62 +44309.77j   160436.24  -9271.05j]
 [ -353612.05+102199.93j   285393.98-106532.76j   -68510.94 +41346.47j]
 [  -38838.71 +92651.41j    91690.36 -77397.41j  -115268.83 +32830.02j]
 [   43618.62 +42461.13j   -29564.08 -43405.79j   -46294.81 +50679.65j]
 [  100883.   -98530.08j   -77818.8  +86290.73j    31552.56 -50487.04j]]

plt.imshow(np.log(abs(F)), cmap='gray');
plt.title(r'$\log \left( |F| \right)$');

Shift it so the DC is in the centre

f = np.array(plt.imread('images/pd.jpg')[:,:,0], dtype=float)
F = fftshift(fft2(f))

plt.figure(figsize=[10,5]);
plt.subplot(1,2,1); plt.imshow(abs(F), cmap='gray'); plt.title(r'$|F|$')
plt.subplot(1,2,2); plt.imshow(np.log(abs(F)), cmap='gray'); plt.title(r'$\log (|F|)$');

What do the lines in the 2D DFT mean?

f = np.array(plt.imread('images/cardoor.jpg')[:,:,0], dtype=float)
plt.imshow(f, cmap='gray');

F = fftshift(fft2(f))
plt.imshow(np.log(abs(F)), cmap='gray');

2D Fourier Basis Vectors

# A 1D Fourier basis vector
G = np.zeros(128)
G[3] = 1.
g = ifft(G)
plt.plot(np.real(g), 'o');

# A 2D Fourier basis vector
F = np.zeros([256,256])
F1 = F.copy()
F1[3,0] = 1.
f1 = ifft2(F1)
plt.figure(figsize=[18,9])
plt.subplot(1,2,1)
plt.imshow(F1, cmap='gray')
plt.subplot(1,2,2)
plt.imshow(np.real(f1), cmap='gray');

# A 2D Fourier basis vector
F = np.zeros([256,256])
F1 = F.copy()
F1[4,12] = 1.0
f1 = ifft2(F1)
plt.figure(figsize=[18,9])
plt.subplot(1,2,1)
plt.imshow(F1, cmap='gray')
plt.subplot(1,2,2)
plt.imshow(np.real(f1), cmap='gray');

F1 = F.copy(); F1[3,5] = 1.0
F2 = F.copy(); F2[2,-3] = 2.0
f1 = ifft2(F1)
f2 = ifft2(F2)
plt.figure(figsize=[18,12]); MN = 256.**2; rg = 2.;
plt.subplot(2,3,1); plt.imshow(np.abs(F1), vmin=0, vmax=rg, cmap='gray'); plt.title('F1')
plt.subplot(2,3,4); plt.imshow(np.real(f1), vmin=-rg/MN, vmax=rg/MN, cmap='gray'); plt.title('f1')
plt.subplot(2,3,3); plt.imshow(np.abs(F2), vmin=0, vmax=rg, cmap='gray'); plt.title('F2')
plt.subplot(2,3,6); plt.imshow(np.real(f2), vmin=-rg/MN, vmax=rg/MN, cmap='gray'); plt.title('f2')
plt.subplot(2,3,2); plt.imshow(np.abs(F1+F2), vmin=0, vmax=rg, cmap='gray'); plt.title('F1+F2')
plt.subplot(2,3,5); plt.imshow(np.real(f1+f2), vmin=-rg/MN, vmax=rg/MN, cmap='gray'); plt.title('f1+f2');

What happens if we take the DFT twice in a row?

plt.imshow(f, cmap='gray'); plt.axis('off');

F = fft2(f)
plt.imshow(np.log(abs(F)), cmap='gray'); plt.axis('off');

G = fft2(fft2(fft2(F)))
plt.imshow(np.real(G), cmap='gray'); plt.axis('off');