Question about fft2(2D Fourier Transform) of a Gaussian function

Hello,

I have a question about fft2 of a Gaussian fuction. The following is my code:

zr = pi*12.5^2/0.85;
w = 12.5*sqrt(1+(10/zr)^2);
R = 10+zr^2/10;

N0 = 622;
n0 = linspace(-250,250,N0);
[x0,y0] = meshgrid(n0,n0);
r0 = sqrt(x0.^2+y0.^2);
Gau0 = exp(-1*r0.^2/w^2).*exp(1i*2*pi/0.85*10).*exp(-1i*((r0.^2)/2*R));
% figure,imagesc(n0,n0,abs(Gau0)),axis square,colormap hot

Fm = fftshift(fft2(Gau0));

figure,imagesc(n0,n0,abs(Fm).^2),axis square,colormap hot
figure,imagesc(n0,n0,angle(Fm)),axis square,colormap hot

Q: The problem is that the FT(Fourier transform) of a Gaussian function should be another Gaussian function. But if I change the sample number (N0 in the code), the amplitude of the FT is strange, a kind of periodic function sometimes. That is to say, the FT of a Gaussian function is related to the sampling number? It is very hard for me to understand such things. Maybe it is because it is undersampling in my code?

Thank you very much for your help (I hope my question is clear), and any help will be much appreciated.

Answer 1

[David Goodmanson]

Hi Zhu Li,

A couple of things are going on here. First, there is the checkerboard pattern in the angle plot. I added an extra plot (fig 3) to show the checkerboard pattern there as well. The same is also true for imag(Fm). The reason is that the gaussian Gau0 does not have its center at the origin of coordinates. The center is shifted, which leads to a multiplicative phase factor in the fft. The shift is as large as can be done for an fft, and the corresponding phase oscillation is as fast as possible (nyquist frequency). So the fft has a mulitplicative factor of [1 -1 1 -1 1 -1 ...] in each spacial frequency coordinate, leading to the checkerboard. (In plot 2 the abs function rather misleadingly got rid of the checkerboard, as happens sometimes with abs).

If Gau0 is shifted to the origin before doing the fft, the checkerboard goes away as in plots 6 and 7.

Your gaussian is of the form exp(-r0^2*(-a+i*b)). The transform looks like

exp(-kr^2*(-a-i*b)/(4*(a^2+b^2))      [1]

where kr is the radial spacial frequency coordinate. Consequently you have the radial oscillations of plot 6 and 7.

A more traditional, real gaussian has b=0 and the fft is shown in plots 8 and 9. The result should be purely real, and indeed the rings are not present in those plots. However, when the size of the fft is down around numerical machine precision, the angle calculation can turn out to be anything. These angles are meaningless. Hence the solid red disc in plot 10 which is angle = 0, surrounded by numerical confetti when the gaussian [1] above is too small.

zr = pi*12.5^2/0.85;
w = 12.5*sqrt(1+(10/zr)^2);
R = 10+zr^2/10;
N0 = 622;
n0 = linspace(-250,250,N0);
[x0,y0] = meshgrid(n0,n0);
r0 = sqrt(x0.^2+y0.^2);
Gau0 = exp(-1*r0.^2/w^2).*exp(1i*2*pi/0.85*10).*exp(-1i*((r0.^2)/2*R));
figure(1),imagesc(n0,n0,abs(Gau0)),axis square,colormap hot

Fm = fftshift(fft2(Gau0));
figure(2),imagesc(n0,n0,abs(Fm)),axis square,colormap hot
figure(3),imagesc(n0,n0,real(Fm)),axis square,colormap hot
figure(4),imagesc(n0,n0,angle(Fm)),axis square,colormap hot

% shift gaussian center to origin
Fmiff = fftshift(fft2(ifftshift(Gau0)));
figure(5),imagesc(n0,n0,abs(Fmiff)),axis square,colormap hot
figure(6),imagesc(n0,n0,real(Fmiff)),axis square,colormap hot
figure(7),imagesc(n0,n0,angle(Fmiff)),axis square,colormap hot

% real gaussian
Gaureal = exp(-1*r0.^2/w^2).*exp(1i*2*pi/0.85*10);
Fmiffreal = fftshift(fft2(ifftshift(Gaureal)));

figure(8),imagesc(n0,n0,abs(Fmiffreal)),axis square,colormap hot
figure(9),imagesc(n0,n0,real(Fmiffreal)),axis square,colormap hot
figure(10),imagesc(n0,n0,angle(Fmiffreal)),axis square,colormap hot

[Zhu Li]

Thank you very much, David.

As you said, it's probably because the phase changes (oscillate) too fast to sample enough for the fft function. But I try changing the sample number (N0) into an odd one. The abs of this function is still the checkboard pattern. So I guess the right "pattern" from your method is just a coincidence? Maybe the reason is the second phase term in the function, which triggers the phase of the function change too fast.

[David Goodmanson]

Hi Zhu Li,

In the case of odd N, the shift in the spacial domain is still quite large, leading to an oscillation in the frequency domain that is still quite high. Not quite as fast as alternating [1 -1 1 -1 ...], but almost, so you still see a checkerboard pattern. The fix is the same as for even N, inserting an ifftshift command.