MATLAB: FFT of Gaussian Pulse in Time Domain

energy spectrum densityfftMATLAB

Referring to the Gaussian Pulse in Time Domain (2nd example of https://uk.mathworks.com/help/matlab/ref/fft.html ) I want to highlight 2 things:
1) In this example you write:
Fs = 100;           % Sampling frequency

t = -0.5:1/Fs:0.5;  % Time vector 

L = length(t);      % Signal length

X = 1/(4*sqrt(2*pi*0.01))*(exp(-t.^2/(2*0.01)));
n = 2^nextpow2(L);
Y = fft(X,n);
f = Fs*(0:(n/2))/n;
P = abs(Y/n);
plot(f,P(1:n/2+1))
but this is wrong. It should instead be:
Fs = 100;           % Sampling frequency
t = -0.5:1/Fs:0.5;  % Time vector 
L = length(t);      % Signal length
X = 1/(4*sqrt(2*pi*0.01))*(exp(-t.^2/(2*0.01)));
n = 2^nextpow2(L);
Y = fft(X,n)
f = Fs*(0:(n/2))/n;
P = abs(Y/L);
plot(f,P(1:n/2+1))
2) I do not get why in this example, differently from the other two examples, you do not multiply by 2 the vector P. In other words, I would expect something like this:
P = abs(Y/L);
P = P(1:n/2+1);
P(2:end-1) = 2*P(2:end-1);
I made a comparison with the analytical solution and indeed it turns out that it is correct not to multiply by 2, but why is that?
Thanks

Best Answer

Hello Francesco,
It appears that you are looking for the kind fft normalization that approximates a Fourier transform. If that is not what you are looking for, then this answer is not germane, but if it is, then
Y(f) = Integral {-inf,inf}  X(t) exp(-i 2 pi f t) dt
which you can approximate with
dt = 1/Fs;
Y = fft(X)*dt;
In not every case is dt is defined by dt = 1/Fs but it is here.
For reasons known only to themselves, Mathworks has a gaussian function whose integral as a continuous function = 1/4. For f=0 the expression for Y above approximates this integral so as you know, the answer should be Y(0) = 1/4.
They chose a different normalization fft(X)/n (less appropriate to this situation) so you don't get that answer. Your expression fft(X)/L does give the answer, but that's a coincidence having to do with the specific limits of the t array. To get the answer of 1/4, the factor should be dt. The code below illustrates that, where the limits of the t array have been changed to +-.7 to remove the coincidence.
Rather than plot out half of the frequency domain, the code plots the entire domain and uses fftshift to put zero frequency at the center of the freq array.
As to the question of doubled amplitudes, you can do that any time including this one, but it only makes sense for a set of real signals, each of the form A_m cos(2 pi f_m t + phi) [maybe times a decay factor exp(-g_m t) ]. In that case the fft(X)/n normalization gives peaks at frequency +-f_m, of height A_m / 2, each with a complex phase factor. People take abs of that, double the size of the peak (the f=0 value should not be doubled but that mistake gets made often) and plot out positive frequencies only. That gives you the amplitudes A_m of the cos functions so you can see a spectrum. But at that point you have thrown away phase information so you can’t do an ifft to get back to the original signal.
In the case of a continuous function like a Gaussian there is no real incentive or purpose to do all that.
Fs = 100;           % Sampling frequency
t = -0.7:1/Fs:0.7;  % Time vector 
L = length(t);      % Signal length
X = 1/(4*sqrt(2*pi*0.01))*(exp(-t.^2/(2*0.01)));
%Plot the pulse in the time domain.
figure(1)
plot(t,X)
title('Gaussian Pulse in Time Domain')
xlabel('Time (t)')
ylabel('X(t)')
n = 2^nextpow2(L);
%Convert the Gaussian pulse to the frequency domain.
Y = fftshift(fft(X,n));
%Define the frequency domain and plot the unique frequencies.
f = Fs*(-n/2:n/2-1)/n;
P1 = abs(Y)/n;         % blue
P2 = abs(Y)/L;         % red
P3 = abs(Y)*(1/Fs);    % yellow
figure(2)
plot(f,P1,f,P2,f,P3) 
title('Gaussian Pulse in Frequency Domain')
xlabel('Frequency (f)')
ylabel('|P(f)|')
Related Question