Understanding Fourier transform example in Matlab

Question

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

Answer 1

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

Answer 2

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.

Answer 3

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.