Differences between FFT and analytical Fourier Transform

In MATLAB, if I perform an FFT on a simple Gaussian spot with a phase grating:

S = exp(-2*((X.^2+Y.^2).^0.5/sigma_x).^2);
Ein = S.*exp(1i*X);
Eout = fftshift(fft2(Ein));
phase = mod(pi+angle(Eout),2*pi/1);
Eout = abs(Eout).*exp(1i*phase);


I find a strange grid like phase in the Fourier plane. The output amplitude is as I'd expect from Fourier optics, but the phase seems unphysical.

I tried a 1D analogue of this case in Mathematica with the analytical Fourier transform and found a flat phase in the Fourier plane:

FTgrat[p_] =FourierTransform[Exp[-((2 ((x)^2) )/\[Sigma])]*(Exp[- I x]), x, p];
a = Abs[FTgrat[p]];
b = Arg[FTgrat[p]];
Grid[{{Plot[a, {p, -6, 6}, PlotRange -> All],
Plot[b, {p, 0, 2}, PlotRange -> All]}}]


Is the grid-like phase a mathematical artifact of an FFT, and can it be avoided?

• Please provide the code (inlcuding plots) you use to facilitate answering. – Yves Klett Jun 9 '15 at 11:38
• ^i.e. the Mathematica code – dr.blochwave Jun 9 '15 at 11:38
• I can't see the titles for the figures I'm uploading after I post them I'm afraid, hope it's not too unclear. (First figures are the Matlbab output / Second Mathematica) – user12800 Jun 9 '15 at 12:12
• As a reminder, Mathematica uses a different normalization for the FFT from MATLAB. Set FourierParameters appropriately. – J. M.'s discontentment Jun 9 '15 at 12:17
• Without more details it's hard to say, but I note that your 1D case has a Gaussian centered at the origin, while your 2D case is centered far from the origin. – John Doty Jun 9 '15 at 14:08

How to make the numerical Fourier transform match up exactly with the analytical one? I'll just work a 1-dimensional example here, but it should apply in 2D as well. We take a Gaussian pulse whose phase is zero at time zero, and it should have a completely real and positive Fourier transformation

pulse[t_] := Exp[-t^2/ (2 σ^2)] Exp[-I ω0 t];
pulseω = FourierTransform[pulse[t], t, ω]
ω0 = 2;
σ = 10;

Grid[
{{Plot[Re@pulse[t], {t, -50, 50}, PlotLabel -> "Re E(t)"]}, {Plot[
Abs@pulseω, {ω, 1.5, 2.5},
PlotLabel -> "Abs E(ω)"],
Plot[Arg@pulseω, {ω, 1.5, 2.5},
PlotLabel -> "Arg E(ω)"]}}]


Now let's do it numerically and make sure we get the same thing.

So we sample in the time domain at a given rate, this corresponds to an inverse sampling rate in the frequency domain, as described in detail here, in less detail here, and probably more definitively on another site. But most of the signal-processing literature I find on the web is opaque to me. I just remember these functions and use them every time,

dt = 0.5;
num = 2^12;
df = 2 \[Pi]/(num dt);
fftshift[flist_] := RotateRight[flist, (Length@flist)/2 - 1];
Print["Frequency Resolution = " <> ToString[df]];
Print["Frequency Range = +/-" <> ToString[num/2 df]];
timevalues := Table[t, {t, -dt num/2 + dt, num/2 dt, dt}];
pulse[t_] := Exp[-t^2/ 100] Exp[-I 2 t];
timelist = pulse /@ timevalues;

(*Frequency Resolution = 0.00306796*)

(*Frequency Range = +/-6.28319 *)


Now we plot the Fourier transform, keeping in mind that we need to shift the zero frequency to the center of the data using fftshift

ftlist = fftshift[Fourier[timelist]];

Grid[{{ListLinePlot[Transpose[{timevalues, Re@timelist}],
PlotRange -> {{-50, 50}, All}, PlotLabel -> "Re E(t)"]},
{ListLinePlot[Abs@ftlist, PlotLabel -> "Abs E(\[Omega])",
DataRange -> df {-num/2 + 1, num/2},
PlotRange -> {{1.5, 2.5}, All}, GridLines -> Automatic]
, ListLinePlot[Arg@ftlist, PlotLabel -> "Arg E(\[Omega])",
PlotRange -> {{1.5, 2.5}, All},
DataRange -> df {-num/2, num/2}]}}]


So clearly we've introduced this phase that is linear in ω, which is what you see in your 2D data. The answer is that you also need to shift your time values so that the zero time is at the beginning. I never see this written about almost anywhere, but without doing it you don't get the right result.

To repeat that, we have to apply a RotateLeft operation on the timevalues, take the FFT, and then apply a RotateRight on the resulting frequency values.

timevalues :=
RotateLeft[Table[t, {t, -dt num/2 + dt, num/2 dt, dt}], num/2 - 1];
timelist = pulse /@ timevalues;

ftlist = fftshift[Fourier[timelist]];

Grid[{{ListLinePlot[Sort@Transpose[{timevalues, Re@timelist}],
PlotRange -> {{-50, 50}, All}, PlotLabel -> "Re E(t)"]},
{ListLinePlot[Abs@ftlist, PlotLabel -> "Abs E(\[Omega])",
DataRange -> df {-num/2 + 1, num/2},
PlotRange -> {{1.5, 2.5}, All}, GridLines -> Automatic]
, ListLinePlot[Arg@ftlist, PlotLabel -> "Arg E(\[Omega])",
PlotRange -> {{1.5, 2.5}, All},
DataRange -> df {-num/2, num/2}]}}]


And now we have the right answer!