2d Fourier Transforms: FFT vs Fourier Optics

Question

I am trying to use programming to increase my understanding of Fourier optics. I know that physically and mathematically the Fourier transform of a Fourier transform is inverted -> F{F{f(x)} = f(-x). I am having two problems 1) The second transform doesn't return anything like the original function except in the simple gaussian case (which makes it even more confusing), and 2) there seems to be some scaling factor that requires me to "zoom in" and distort the transformed image to a point that it is much less helpful (as illustrated below). **Editted with suggestions from @Cris Luengo

#%% Playing with 2d Fouier Transform

import numpy as np
from scipy import fftpack
import matplotlib.pyplot as plt
import LightPipes

wavelength = 792*nm
size = 15*mm
N = 600
w0=3*mm

# Fields
sq = np.zeros([100,100])
sq[25:75, 25:75] = 1
F=Begin(size,wavelength,N)
I0 = Intensity(0,GaussBeam(F, w0, LG=True, n=0, m=0))
I1 = Intensity(0,GaussBeam(F, w0, LG=False, n=0, m=1))+Intensity(0,GaussBeam(F, w0, LG=False, n=1, m=0))

# Plot transforms
f = sq
F = np.fft.fftshift(fftpack.fft2(f))
F_F = fftpack.fft2((F))

plt.subplot(331), plt.imshow(f)
plt.title(r'f'), plt.xticks([]), plt.yticks([])
plt.subplot(332), plt.imshow(np.abs(F))
plt.title(r'F\{f\}'), plt.xticks([]), plt.yticks([])
plt.subplot(333), plt.imshow(np.abs(F_F))
plt.title('F\{F\{f\}\}'), plt.xticks([]), plt.yticks([])

# plt.subplot(331), plt.imshow(f)
# plt.title(r'f'), plt.xticks([]), plt.yticks([])
# plt.subplot(332), plt.imshow(np.abs(F))
# plt.title(r'F\{f\}'), plt.xticks([]), plt.yticks([])
# plt.subplot(333), plt.imshow(np.abs(F_F))
# plt.title('F\{F\{f\}\}'), plt.xticks([]), plt.yticks([])

f = I0
F = np.fft.fftshift(fftpack.fft2(f))
F_F = fftpack.fft2((F))

plt.subplot(334), plt.imshow(f)
plt.title(r'f'), plt.xticks([]), plt.yticks([])
plt.subplot(335), plt.imshow(np.abs(F))
plt.title(r'F\{f\}'), plt.xticks([]), plt.yticks([])
plt.subplot(336), plt.imshow(np.abs(F_F))
plt.title('F\{F\{f\}\}'), plt.xticks([]), plt.yticks([])

f = I1
F = fftpack.fft2(f)
F_F = fftpack.fft2(F)

plt.subplot(337), plt.imshow(f)
plt.title(r'f'), plt.xticks([]), plt.yticks([])
plt.subplot(338), plt.imshow(np.abs(F))
plt.title(r'F\{f\}'), plt.xticks([]), plt.yticks([])
plt.subplot(339), plt.imshow(np.abs(F_F))
plt.title('F\{F\{f\}\}'), plt.xticks([]), plt.yticks([])

plt.tight_layout()
plt.show()

Answer 1

After chatting with Cris, it seems that there is no scaling factor, this type of DFT just works this way it seems. So the solution I have found is to increase the pixels to the point that I can zoom in and have a clear enough image. It's not a great solution but paired with LightPipes it is now possible to get an idea of what the transform of light modes will look like, as well as illustrate that at the image plane of a lens system they will appear as they did in the front focal field.

#%% Playing with 2d Fourier Transform

import numpy as np
import matplotlib.pyplot as plt
import LightPipes
from scipy.fftpack import fft2 as fft
from numpy.fft import fftshift, ifftshift

wavelength = 792*nm
size = 100*mm
N = 1000
w0=3*mm

# Fields
sq = np.zeros([100,100])
sq[25:75, 25:75] = 1
F=Begin(size,wavelength,N)
I0 = Intensity(0,GaussBeam(F, w0, LG=True, n=0, m=0))
I1 = Intensity(0,GaussBeam(F, w0, LG=False, n=0, m=1))+Intensity(0,GaussBeam(F, w0, LG=False, n=1, m=0))

# Plot transforms
f = sq
F = fftshift(fft(ifftshift(f)))
F_F = fftshift(fft(ifftshift(F)))

plt.subplot(331), plt.imshow(f,cmap = cmap)
plt.title(r'f'), plt.xticks([]), plt.yticks([])
plt.subplot(332), plt.imshow(np.abs(F),cmap = cmap)
plt.title(r'F\{f\}'), plt.xticks([]), plt.yticks([])
plt.subplot(333), plt.imshow(np.abs(F_F),cmap = cmap)
plt.title('F\{F\{f\}\}'), plt.xticks([]), plt.yticks([])

f = I0
F = fftshift(fft(ifftshift(f)))
F_F = fftshift(fft(ifftshift(F)))

plt.subplot(334), plt.imshow(f[450:550,450:550],cmap = cmap)
plt.title(r'f'), plt.xticks([]), plt.yticks([])
plt.subplot(335), plt.imshow(np.abs(F)[450:550,450:550],cmap = cmap)
plt.title(r'F\{f\}'), plt.xticks([]), plt.yticks([])
plt.subplot(336), plt.imshow(np.abs(F_F)[450:550,450:550],cmap = cmap)
plt.title('F\{F\{f\}\}'), plt.xticks([]), plt.yticks([])

f = I1
F = fftshift(fft(ifftshift(f)))
F_F = fftshift(fft(ifftshift(F)))

plt.subplot(337), plt.imshow(f[450:550,450:550],cmap = cmap)
plt.title(r'f'), plt.xticks([]), plt.yticks([])
plt.subplot(338), plt.imshow(np.abs(F)[450:550,450:550],cmap = cmap)
plt.title(r'F\{f\}'), plt.xticks([]), plt.yticks([])
plt.subplot(339), plt.imshow(np.abs(F_F)[450:550,450:550],cmap = cmap)
plt.title('F\{F\{f\}\}'), plt.xticks([]), plt.yticks([])

plt.tight_layout()
plt.show()