1

I have Gaussian beam in 2D:

2D Gaussian

After doing fft2 and angle I get strange results:

magnitude of FFT

phase of FFT

def finite2D(x,y, N, M, a, hx):
    f = np.array([[0.0]*N]*N)
    for i in range(len(x)):
        for k in range(len(y)):
            f[i][k] = np.exp(-(x[i]*x[i] + y[k]*y[k]))

    D1 = fftpack.fft2(f)
    D2 = fftpack.fftshift(D1)

    b = N*N/(4*a*M)
    x = np.linspace(-b, b, N)
    y = np.linspace(-b, b, N)
    xx, yy = np.meshgrid(x, y)
    plt.imshow(np.abs(D2))
    plt.show()

    plt.imshow(np.angle(D2))
    plt.show(True)
    return D2, phas

a = 5
N = 128
M = 256
b = N*N/(4*a*M)
hx = 2*a/N
x = np.linspace(-a, a, N)
y = np.linspace(-a, a, N)
finite2D(x,y, N, M, a, hx)

It should be phase 0 or close to 0. Why is this not the case, and how do I fix this?

///Updated:

def finite2D(x,y, N, M, a, hx):
    f = np.array([[0.0]*N]*N)
    for i in range(len(x)):
        for k in range(len(y)):
            f[i][k] = np.exp(-(x[i]*x[i] + y[k]*y[k]))

    f = fftpack.ifftshift(f)
    D1 = fftpack.fft2(f)
    D2 = fftpack.fftshift(D1)

    b = N*N/(4*a*M)
    x = np.linspace(-b, b, N)
    y = np.linspace(-b, b, N)
    xx, yy = np.meshgrid(x, y)
    plt.imshow(np.abs(D2))
    plt.show()

    plt.imshow(np.angle(D2))
    plt.show(True)
    return D2

a = 5
N = 128
M = 256
b = N*N/(4*a*M)
hx = 2*a/N
x = np.linspace(-a, a, N, endpoint=False)
y = np.linspace(-a, a, N, endpoint=False)
finite2D(x,y, N, M, a, hx)

Phase: Phase

SoH
  • 49
  • 7
  • Possible duplicate of [Analytical Fourier transform vs FFT of functions in Matlab](https://stackoverflow.com/questions/49317834/analytical-fourier-transform-vs-fft-of-functions-in-matlab) – Cris Luengo Apr 07 '19 at 13:08
  • @CrisLuengo I do not quite understand how this can solve my problem. Could you help me fix my code? – SoH Apr 07 '19 at 13:57

1 Answers1

1

The FFT asumes that the origin is in the top-left corner of the image. Thus, you are computing the FFT of a Gaussian shifted by half the image size. This shift leads to a high-frequency phase shift in the frequency domain.

To solve the problem, you need to shift the origin of your Gaussian signal to the top-left corner of the image. ifftshift does this:

f = fftpack.ifftshift(f)
D1 = fftpack.fft2(f)
D2 = fftpack.fftshift(D1)

Note that where the magnitude is very small, the phase is defined by rounding errors, don’t expect zero phase there.

The updated result looks good, but there still is a very small gradient in the central region. This is caused by the half-pixel shift of the Gaussian. This shift is given by the definition of the x and y coordinates:

N = 128
x = np.linspace(-a, a, N)
y = np.linspace(-a, a, N)

For an even-sized N, do

x = np.linspace(-a, a, N, endpoint=False)
y = np.linspace(-a, a, N, endpoint=False)

such that there is a sample where x==0.

Cris Luengo
  • 55,762
  • 10
  • 62
  • 120
  • @SoH: Please check updated answer. I didn’t pay attention to details after finding the big obvious issue. :) – Cris Luengo Apr 07 '19 at 14:35
  • @SoH: Sorry for not being explicit, both the `x` and the `y` need the same treatment. Notice how your change of `x` removed one of the gradients. Changing `y` in the same way will remove the other gradient. – Cris Luengo Apr 07 '19 at 15:33