Comparing FFT of Function to Analytical FT Solution in Matlab

Question

I am trying to compare the FFT of exp(-t^2) to the function's analytical fourier transform, exp(-(w^2)/4)/sqrt(2), over the frequency range -3 to 3.

I have written the following matlab code and have iterated on it MANY times now with no success.

fs = 100; %sampling frequency
dt = 1/fs;
t = 0:dt:10-dt; %time vector
L = length(t); %number of sample points
%N = 2^nextpow2(L);  %necessary?

y = exp(-(t.^2));
Y=dt*ifftshift(abs(fft(y)));

freq = (-L/2:L/2-1)*fs/L; %freq vector

F = (exp(-(freq.^2)/4))/sqrt(2); %analytical solution

%Y_valid_pts = Y(W>=-3 & W<=3); %compare for freq = -3 to 3
%npts = length(Y_valid_pts);
% w = linspace(-3,3,npts);
% Fe = (exp(-(w.^2)/4))/sqrt(2);

error = norm(Y - F)  %L2 Norm for error

hold on;
plot(freq,Y,'r');
plot(freq,F,'b');
xlabel('Frequency, w');
legend('numerical','analytic');
hold off;

You can see that right now, I am simply trying to get the two plots to look similar. Eventually, I would like to find a way to do two things: 1) find the minimum sampling rate, 2) find the minimum number of samples, to reach an error (defined as the L2 norm of the difference between the two solutions) of 10^-4.

I feel that this is pretty simple, but I can't seem to even get the two graphs visually agree. If someone could let me know where I'm going wrong and how I can tackle the two points above (minimum sampling frequency and minimum number of samples) I would be very appreciative.

Thanks

Why do you multiply by dt in the line "Y=dt*ifftshift(abs(fft(y)));" Is it because you are trying to compare the result to the continuous Fourier case, which has a dt in the integral? — teeeeee
In the past, I have always divided the output of the fft by the number of samples - what is the difference, and when would you use each method? — teeeeee

SleuthEye SleuthEye · Accepted Answer · 2014-10-11T00:50:19

A first thing to note is that the Fourier transform pair for the function exp(-t^2) over the +/- infinity range, as can be derived from tables of Fourier transforms is actually:

Gaussian FT: exp(-t^2) <==> exp(-(w^2)/4)*sqrt(pi)

Finally, as you are generating the function exp(-t^2), you are limiting the range of t to positive values (instead of taking the whole +/- infinity range). For the relationship to hold, you would thus have to generate exp(-t^2) with something such as:

t = 0:dt:10-dt;     %time vector
t = t - 0.5*max(t); %center around t=0
y = exp(-(t.^2));

Then, the variable w represents angular frequency in radians which is related to the normalized frequency freq through:

w = 2*pi*freq;

Thus,

F = (exp(-((2*pi*freq).^2)/4))*sqrt(pi); %analytical solution

Comparing FFT of Function to Analytical FT Solution in Matlab

1 Answers