Question about fft2(2D Fourier Transform) of a Gaussian function
39 views (last 30 days)
Show older comments
Zhu Li
on 19 Jul 2019
Commented: David Goodmanson
on 26 Jul 2019
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.
0 Comments
Accepted Answer
David Goodmanson
on 21 Jul 2019
Edited: David Goodmanson
on 21 Jul 2019
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
2 Comments
David Goodmanson
on 26 Jul 2019
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.
More Answers (0)
See Also
Categories
Find more on Spectral Measurements in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!