Using FFT for 3D array representation of 2D field

Question

I need to obtain the fourier transform of a complex field. I'm using python.

My input is a 2D snapshot of the electric field in the xy-plane.

I currently have a 3D array F[x][y][z] where F[x][y][0] contains the real component and F[x][y]1 contains the complex component of the field.

My current code is very simple and does this:

result=np.fft.fftn(F)
result=np.fft.fftshift(result)

I have the following questions:

1) Does this correctly compute the fourier transform of the field, or should the field be entered as a 2D matrix with each element containing both the real and imaginary component instead?

2) I entered the complex component values of the field using the real multiple only (i.e if the complex value is 6i I entered 6), is this correct or should this be entered as a complex value instead (i.e. entered as '6j')?

3) As this is technically a 2D input field, should I use np.fft.fft2 instead? Doing this means the output is not centered in the middle.

4) The output does not look like what I'd expect the fourier transform of F to look like, and I'm unsure what I'm doing wrong.

Full example code:

import numpy as np
import matplotlib.pyplot as plt
from matplotlib import cm

x, y = np.meshgrid(np.linspace(-1,1,100), np.linspace(-1,1,100))
d = np.sqrt(x*x+y*y)
sigma, mu = .35, 0.0
g1 = np.exp(-( (d-mu)**2 / ( 2.0 * sigma**2 ) ) )

F=np.empty(shape=(300,300,2),dtype=complex)
for x in range(0,300):
        for y in range(0,300):
            if y<50 or x<100 or y>249 or x>199:
                    F[x][y][0]=g1[0][0]
                    F[x][y][1]=0j
            elif y<150:
                    F[x][y][0]=g1[x-100][y-50]
                    F[x][y][1]=0j
            else:
                    F[x][y][0]=g1[x-100][y-150]
                    F[x][y][1]=0j

F_2D=np.empty(shape=(300,300))
for x in range(0,300):
    for y in range(0,300):
            F_2D[x][y]=np.absolute(F[x][y][0])+np.absolute(F[x][y][1])

plt.imshow(F_2D)
plt.show()

result=np.fft.fftn(F)
result=np.fft.fftshift(result)

result_2D=np.empty(shape=(300,300))
for x in range(0,300):
    for y in range(0,300):
            result_2D[x][y]=np.absolute(result[x][y][0])+np.absolute(result[x][y][1])

plt.imshow(result_2D)
plt.show()

plotting F gives this:

With np.fft.fftn, the image shown at the end is:

And with np.fft.fft2:

Neither of these look like what I would expect the fourier transform of F to look like.

Answer 1

I add here another answer, suitable to the added code.

The answer is still np.fft.fft2(). Here's an example. I modified the code slightly. To verify that we need fft2 I discarded one of the blobs, and then we know that a single Gaussian blob should transform into a Gaussian blob (with a certain phase, that's not shown when plotting absolute value). I also decreased the standard deviation so that the frequency response will widen a little.

Code:

import numpy as np
import matplotlib.pyplot as plt

x, y = np.meshgrid(np.linspace(-1,1,100), np.linspace(-1,1,100))
d = np.sqrt(x**2+y**2)
sigma, mu = .1, 0.0
g1 = np.exp(-( (d-mu)**2 / ( 2.0 * sigma**2 ) ) )
N = 300
positions = [ [150,100] ]#, [150,200] ]
sz2 = [int(x/2) for x in g1.shape]
F_2D = np.zeros([N,N])
for x0,y0 in positions:
    F_2D[ x0-sz2[0]: x0+sz2[0], y0-sz2[1]:y0+sz2[1] ] = g1 + 1j*0.

result = np.fft.fftshift(np.fft.fft2(F_2D))

plt.subplot(211); plt.imshow(F_2D)
plt.subplot(212); plt.imshow(np.absolute(result))
plt.title('$\sigma$=.1')
plt.show()

Result:

To get back to the original problem, we need only change

positions = [ [150,100] , [150,200] ] and sigma=.35 instead of sigma=.1.

Answer 2

You should use complex numpy variables (by using 1j) and use fft2. For example:

N = 16
x0 = np.random.randn(N,N,2)
x = x0[:,:,0] + 1j*x0[:,:,1]
X = np.fft.fft2(x)

Using fftn on x0 will do a 3D FFT, and using fft will do vector-wise 1D FFT.