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(ks)
[ 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)
F = np.fft.fft(f)
print(np.round(F, 2))
Fs = fftshift(F)
print(np.round(Fs, 2))
[ 1.0000 2.0000 3.0000 4.0000 5.0000 6.0000 7.0000] [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] [-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 ]
f = np.array(plt.imread('pd.jpg')[:,:,0], dtype=float)
plt.imshow(f,cmap='gray');
F = fft2(f)
# Plot the modulus of the Fourier coefs
# DC is in the top-left corner
plt.figure(figsize=[8,8]);
plt.imshow(abs(F), cmap='gray');
plt.title(r"Modulus of the Fourier coefficients, $|F|$");
print(np.round(F[:3, :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]]
plt.imshow(np.log(abs(F)), cmap='gray');
plt.title(r'$\log \left( |F| \right)$');
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|)$');
# A 1D Fourier basis vector
G = np.zeros(128)
G[3] = 1.
g = ifft(G)
plt.plot(np.real(g), 'o');
# A 2D Fourier vasis vector
F = np.zeros([256,256])
F1 = F.copy()
F1[5,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[2,6] = 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');
plt.imshow(f, cmap='gray'); plt.axis('off');
F = fft2(f)
plt.imshow(np.log(abs(F)), cmap='gray'); plt.axis('off');
G = fft2(F)
plt.imshow(np.real(G), cmap='gray'); plt.axis('off');