Fourier Transform in Python 2D

Question

I want to perform numerically Fourier transform of Gaussian function using fft2. Under this transformation the function is preserved up to a constant.

I create 2 grids: one for real space, the second for frequency (momentum, k, etc.). (Frequencies are shifted to zero). I evaluate functions and eventually plot the results.

Here is my code

import numpy as np
import matplotlib.pyplot as plt
from scipy.fftpack import fft2, ifft2
from mpl_toolkits.mplot3d import Axes3D

"""CREATING REAL AND MOMENTUM SPACES GRIDS"""
N_x, N_y = 2 ** 11, 2 ** 11
range_x, range_y = np.arange(N_x), np.arange(N_y)
dx, dy = 0.005, 0.005
# real space grid vectors
xv, yv = dx * (range_x - 0.5 * N_x), dy * (range_y - 0.5 * N_y)
dk_x, dk_y = np.pi / np.max(xv), np.pi / np.max(yv)
# momentum space grid vectors, shifted to center for zero frequency
k_xv, k_yv = dk_x * np.append(range_x[:N_x//2], -range_x[N_x//2:0:-1]), \
            dk_y * np.append(range_y[:N_y//2], -range_y[N_y//2:0:-1])

# create real and momentum spaces grids
x, y = np.meshgrid(xv, yv, sparse=False, indexing='ij')
kx, ky = np.meshgrid(k_xv, k_yv, sparse=False, indexing='ij')

"""FUNCTION"""
f = np.exp(-0.5 * (x ** 2 + y ** 2))
F = fft2(f)
f2 = ifft2(F)
"""PLOTTING"""
fig = plt.figure()
ax = Axes3D(fig)
surf = ax.plot_surface(x, y, np.abs(f), cmap='viridis')
# for other plots I changed to
# surf = ax.plot_surface(kx, ky, np.abs(F), cmap='viridis')
# surf = ax.plot_surface(x, y, np.abs(f2), cmap='viridis')
plt.show()

So, the plots for gaussian, fourier(gaussian), inverse_fourier(fourier(gaussian)) are the following:Initial, Fourier, Inverse Fourier

Using plt.imshow(), I additionally plot fourier of gaussian:

   plt.imshow(F)
   plt.colorbar()
   plt.show()

The result is as follows: imshow

That doesn't make sense. I expect see the same gaussian function as the initial up to some constant order of unity.

I would be very glad if someone could clarify this for me.

Solution

I think you are a bit puzzled by the shape of your output F. Especially, you might wonder why you see such a sharp peak and not a wide-spread gaussian.

I changed your code a little bit:

 import numpy as np
 import matplotlib.pyplot as plt
 from scipy.fftpack import fft2, ifft2
 from mpl_toolkits.mplot3d import Axes3D

 """CREATING REAL AND MOMENTUM SPACES GRIDS"""
 N_x, N_y = 2 ** 10, 2 ** 10
 range_x, range_y = np.arange(N_x), np.arange(N_y)
 dx, dy = 0.005, 0.005
 # real space grid vectors
 xv, yv = dx * (range_x - 0.5 * N_x), dy * (range_y - 0.5 * N_y)
 dk_x, dk_y = np.pi / np.max(xv), np.pi / np.max(yv)
 # momentum space grid vectors, shifted to center for zero frequency
 k_xv, k_yv = dk_x * np.append(range_x[:N_x//2], -range_x[N_x//2:0:-1]), \
             dk_y * np.append(range_y[:N_y//2], -range_y[N_y//2:0:-1])

 # create real and momentum spaces grids
 x, y = np.meshgrid(xv, yv, sparse=False, indexing='ij')
 kx, ky = np.meshgrid(k_xv, k_yv, sparse=False, indexing='ij')

 """FUNCTION"""
 sigma=0.05
 f = 1/(2*np.pi*sigma**2) * np.exp(-0.5 * (x ** 2 + y ** 2)/sigma**2)
 F = fft2(f)
 """PLOTTING"""
 fig = plt.figure()
 ax = Axes3D(fig)
 surf = ax.plot_surface(x, y, np.abs(f), cmap='viridis')
 # for other plots I changed to
 fig2 = plt.figure()
 ax2 =Axes3D(fig2)
 surf = ax2.plot_surface(kx, ky, np.abs(F)*dx*dy, cmap='viridis')
 plt.show()

Notice that I introduced a sigma parameter to control the width of the gaussian. I now invite you to play with the following parameters: N_x and N_y, d_x and d_y and sigma.

You should then see the inverse behaviour of gaussian in real-space and in fourier space: The larger the gaussian in real-space, the narrower in fourier-space and vice-versa.

So with the currently set parameters in my code, you get the following plots:

Real space:

Fourier Space:

Comments

[Merlin1896]

Duplicate of stackoverflow.com/questions/19975030/… ?

[Max Borovkov]

@Merlin1896 seems like multiplication (discussed in this link) changes nothing. So, it's not a point, I think

[Merlin1896]

Well, multiplying F by dy*dx/2/np.pi certainly does change something: Now all values of F lie between zero and one.

[pluto]

How can I use the 'surf' plot to the 2D DFT in the link stackoverflow.com/questions/70768384/…