Analytical Fourier transform vs FFT of functions in Matlab

The plots at the bottom of the question are not mirrored. If you plot those using lines instead of dots you’ll see the numeric results have very high frequencies. The absolute component matches, but the phase doesn’t. When this happens, it’s almost certainly a case of a shift in the time domain.

And indeed, you define the time domain function with the origin in the middle. The FFT expects the origin to be at the first (leftmost) sample. This is what ifftshift is for:

Y = dt*fftshift(fft(ifftshift(y)));

ifftshift moves the origin to the first sample, in preparation for the fft call, and fftshift moves the origin of the result to the middle, for display.


Edit

Your t does not have a 0:

>> t(L/2+(-1:2))
ans =
  -1.5000e-05  -5.0000e-06   5.0000e-06   1.5000e-05

The sample at t(floor(L/2)+1) needs to be 0. That is the sample that ifftshift moves to the leftmost sample. (I use floor there in case L is odd in size, not the case here.)

To generate a correct t do as follows:

fs = 1e5; % sampling frequency
L = 30 * fs;
t = -floor(L/2):floor((L-1)/2);
t = t / fs;

I first generate an integer t axis of the right length, with 0 at the correct location (t(floor(L/2)+1)==0). Then I convert that to seconds by dividing by the sampling frequency.

With this t, the Y as I suggest above, and the rest of your code as-is, I see this for the Gaussian example:

>> max(abs(F-Y))
ans =    4.5254e-16

For the other function I see larger differences, in the order of 6e-6. This is due to the inability to sample the Heaviside function. You need t=0 in your sampled function, but H doesn’t have a value at 0. I noticed that the real component has an offset of similar magnitude, which is caused by the sample at t=0.

Typically, the sampled Heaviside function is set to 0.5 for t=0. If I do that, the offset is removed completely, and max difference for the real component is reduced by 3 orders of magnitude (largest errors happen for values very close to 0, where I see a zig-zag pattern). For the imaginary component, the max error is reduced to 3e-6, still quite large, and is maximal at high frequencies. I attribute these errors to the difference between the ideal and sampled Heaviside functions.

You should probably limit yourself to band-limited functions (or nearly-band-limited ones such as the Gaussian). You might want to try to replace the Heaviside function with an error function (integral of Gaussian) with a small sigma (sigma = 0.8 * fs is the smallest sigma I would consider for proper sampling). Its Fourier transform is known.

Leave a Comment