Understanding Fourier transform example in Matlab I'm studying about Fourier series and transform and I get confused with the following Matlab example of Fourier transformation: 
Fs = 1000;                    % Sampling frequency
T = 1/Fs;                     % Sample time
L = 1000;                     % Length of signal
t = (0:L-1)*T;                % Time vector
% Sum of a 50 Hz sinusoid and a 120 Hz sinusoid
x = 0.7*sin(2*pi*50*t) + sin(2*pi*120*t); 
y = x + 2*randn(size(t));     % Sinusoids plus noise
plot(Fs*t(1:50),y(1:50))
title('Signal Corrupted with Zero-Mean Random Noise')
xlabel('time (milliseconds)') 

This gives:

I have no problem with this part, but the Fourier transform is where I get lost:
NFFT = 2^nextpow2(L); % Next power of 2 from length of y
Y = fft(y,NFFT)/L;
f = Fs/2*linspace(0,1,NFFT/2+1);

% Plot single-sided amplitude spectrum.
plot(f,2*abs(Y(1:NFFT/2+1))) 
title('Single-Sided Amplitude Spectrum of y(t)')
xlabel('Frequency (Hz)')
ylabel('|Y(f)|')

, which gives:

This is the code which confuses me: 
NFFT = 2^nextpow2(L); % Next power of 2 from length of y
Y = fft(y,NFFT)/L;
f = Fs/2*linspace(0,1,NFFT/2+1);

% Plot single-sided amplitude spectrum.
plot(f,2*abs(Y(1:NFFT/2+1)))

I don't understand all the confusing multiplications and divisions done here....can someone explain why we are for example here: plot(f,2*abs(Y(1:NFFT/2+1))) multiplying by two etc. etc. Why not just do:  plot(f, abs(Y(1:NFFT/2+1))). Why are we dividing by L in this line: Y = fft(y,NFFT)/L;. What are we doing in this line: f = Fs/2*linspace(0,1,NFFT/2+1);??
The definitions for the Matlab functions are the following: 
p = nextpow2(A) returns the smallest power of two that is greater than or equal to the absolute value of A. (That is, p that satisfies 2^p >= abs(A)). 

abs(X) returns an array Y such that each element of Y is the absolute value of the corresponding element of X.

Y = fft(X,n) returns the n-point DFT. fft(X) is equivalent to fft(X, n) where n is the size of X in the first nonsingleton dimension. If the length of X is less than n, X is padded with trailing zeros to length n. If the length of X is greater than n, the sequence X is truncated. When X is a matrix, the length of the columns are adjusted in the same manner.

y = linspace(a,b,n) generates a row vector y of n points linearly spaced between and including a and b. For n < 2, linspace returns b.

In other words I would just want someone to better explain the four given code lines to me. What are we doing an why :) 
Thank you for any help!
P.S.
here is the reference:
http://www.mathworks.se/help/matlab/ref/fft.html
 A: Read up on the differences between the DFT (i.e., discrete fourier transform), the FFT (a fast version of DFT), the *Fourier Series and the Fourier Transform (i.e. the time-continuous case).
Basically, fft computes the DFT. The DFT is simply an invertible linear map from $\mathbb{C}^n$ to itself, i.e. you may imagine it as a change of base. For real-valued inputs, the right half of the result of the DFT is always the conjugated mirror-image of the left part, i.e. for a $N$-point DFT you have that $F(n) = \overline{F(N-n)}$. This is why the code ignores the right half of Y.
The DFT also loses all information about the time scale, since the input is just a vector of real or complex-valued samples. To nevertheless display actual frequencies in Hz, the code needs to recover those - that's what the linspacecode is about. It produces the labels for the $x$-axis, nothing more.
A basic invariant of the DFT (and also other kinds of fourier transforms) is parseval's theorem (http://en.wikipedia.org/wiki/Parseval%27s_theorem). Scaling the amplitudes by a factor of two ensure that this theorem holds for the plotted value - the scaling factor compensates for ignoring the right half of the result. 
A: Based on Parseval's theorem
$${\displaystyle \sum _{n=0}^{N-1}|x[n]|^{2}={\frac {1}{N}}\sum _{k=0}^{N-1}|X[k]|^{2}}$$
This expresses that the energies in the time- and frequency-domain are the same. It means the magnitude $|X[k]|$ of each frequency bin $k$ is contributed by $N$ samples. In order to find out the average contribution by each sample, the magnitude is normalized as $|X[k]|/N$, which leads to
$${\displaystyle \frac{1}{N}\sum _{n=0}^{N-1}|x[n]|^{2}=\sum _{k=0}^{N-1}\left|\frac{X[k]}{N}\right|^{2}},$$
where the LHS is the power of the signal.
The multiplication by two is due to the conversion from the two-sided to single-side spectrum. However, it was wrong in the old version as it should exclude Y(1) which corresponds to the zero frequency and is fixed in the newer version as Y(2:end-1) = 2*Y(2:end-1);
The answer by OmidS was misleading as the MATLAB example is actually NOT wrong, at least for the normalization part Y = fft(y,NFFT)/L; In the newer version of Matlab (I am using 2020a), the example has been modified so that L is not equal to Fs and Y=fft(y)/L.
However, he explained the conversion from the CTFT and DTFT.
A: because sometimes for fast fourier multiplication,scientist use length which is power of two,or some length $k=2^l$,this is idea of  fast fourier transform,
http://www.mathworks.com/matlabcentral/newsreader/view_thread/243154
http://www.mathworks.com/company/newsletters/articles/faster-finite-fourier-transforms-matlab.html
A: I think this question deserves a more satisfactory answer. The truth is that the MATLAB example is actually wrong in dividing the fft by the signal length in the time domain (which is L):
Y = fft(y,NFFT)/L; % The MATLAB example which is actually wrong

The right scaling needed to adhere to Parseval's theorem would be dividing the Fourier transform by the sampling frequency:
Y = fft(y,NFFT)/Fs; % The Correct Scaling

Incidentally, these two are the same in this MATLAB example! (L = Fs = 1000)

It's important to see why the contradiction to Parseval's theorem is caused in the first place. The MATLAB fft function is not to be blamed here! Actually, it adheres to the Parseval's theorem in the simple case of calculating the FFT with the same number of points as the signal length:
x = rand(10000,1); % For any vector (of any length)
A = sum( abs(x).^2 ) - sum( abs(fft(x)).^2 )/length(x); 

For any vector (of any length), A in the above code would be extremely close to 0 (any difference is due to numerical computation errors) as stated by the Parseval's theorem for discrete Fourier transform (DFT):
$$
\sum_{n=0}^{N-1} | x[n] |^2  =   \frac{1}{N} \sum_{k=0}^{N-1} | X[k] |^2
$$
So who is the culprit?! The discrepancy happens only when fft is calculated with a length other than the original input vector's length, which is specified via the second argument of the fft function. But still, the fft function is not doing anything wrong. It simply pads the input vector with enough zeros (or truncates it if the vector length is more than the demanded fft length) and then performs the fft. The Parseval's theorem holds for the padded (truncated) signal but we shouldn't expect it to hold for the original signal (because of dependency on N in Parseval's theorem). 
So what can we do now in this case where the signal and its DFT are of different lengths? We can emulate the continuous case of Parseval's theorem: 
$$
\int_{-\infty}^\infty | x(t) |^2 \, dt   =   \int_{-\infty}^\infty | X(f) |^2 \, df
$$
We can approximate the integral by a sum. The differential elements should be:
dt = 1/Fs;
df = Fs/NFFT; % NFFT is the length of FFT; 
% df is the above because the last point in fft corresponds 
% to Fs [Hz] in the continuous Fourier transform 
% so each fft step is equivalent to Fs/NFFT [Hz]

B = sum( abs(y).^2 )*dt - sum( abs(fft(y, NFFT)/ScalingFactor).^2 )*df

Assuming that we want Parseval's theorem to hold, we need to set the ScalingFactor so that B would be zero. Setting it as Fs would exactly do that:
ScalingFactor = Fs; % The output of fft has to be divided by this

The output of the fft function has to be divided by the sampling frequency of the input vector (Fs) (or equivalently multiplied by dt = 1/Fs) for the Parseval's theorem to hold!

Credits go to Elige Grant for writing about this way of scaling here
and to Rick Rosson for this
A: I also think the example is incorrect, but for a different reason. I'm OK with scaling by L for the following reason. The average signal "power" is given by
$$
P_\mathrm{av}=\frac{1}{T}\sum_{n=0}^{L-1}|x(n)|^2\tau_s \\
=\frac{1}{f_sT}\sum_{n=0}^{L-1}|x(n)|^2 \\
=\frac{1}{L}\sum_{n=0}^{L-1}|x(n)|^2
$$
Combining this with the DFT version of Parseval's theorem yields
$$
P_\mathrm{av}=\frac{1}{L^2}\sum_{k=0}^{L-1}|X(k)|^2 \\
=\sum_{k=0}^{L-1}|X(k)/L|^2
$$
where $X$ is the DFT of $x$. Thus scaling by L seems convenient. I have a problem, however, with the example of how to determine the single-sided version. Since $X(L-k)=X^*(k)$ we have
$$
\sum_{k=0}^{L-1}|X(k)/L|^2=\sum_{k=0}^{L/2}|X(k)/L|^2+\sum_{k=L/2+1}^{L-1}|X(k)/L|^2 \\
=\sum_{k=0}^{L/2}|X(k)/L|^2+\sum_{k=L/2+1}^{L-1}|X^*(L-k)/L|^2 \\
=\sum_{k=0}^{L/2}|X(k)/L|^2+\sum_{k=1}^{L/2-1}|X(k)/L|^2 \\
=|X(0)|^2+2\sum_{k=1}^{L/2-1}|X(k)/L|^2+|X(L/2)|^2
$$
This suggests that in the one-sided version the "interior" components should by multiplied by $\sqrt{2}$ rather than $2$.
My DSP bible does not address this issue and I cannot quote any other authority, so I have to submit this as only my opinion.
