An attempt on explaining 2 scaling factors seen in the Matlab help on FFT
By Nasser Abbasi
Updated 9/17/09
Matlab help on FFT (doc FFT) shows this example
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)')
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)|')
We need to talk above the lines with bold face in the code above. Why in the line
Y = fft(y,NFFT)/L;
The code divides the FFT result by L?
And in the second bold line
plot(f,2*abs(Y(1:NFFT/2+1)))
Why does it multiply by 2?
Y = fft(y,NFFT)/L;
This scaling factor L comes from Nyquist sampling theory.
Assume the analog signal is x(t) with corresponding Fourier transform X(w). Let
us now sample x(t) with some sampling frequency. This will result in a sequence
of numbers: the sampled sequence x(n*T) where T is the sampling period and
n=0,1,…
Now let the sequence of complex numbers Y(k) represent the Fourier transform of
the sampled sequence, i.e. the DFT of x(n*T).
If the sampling frequency is larger than twice the largest frequency in the
signal x(t) then the magnitude of X(w) will be proportional to the magnitude of
Y(k).
The proportionality factor turns out to be the sampling period T.
i.e. we can write
X(w) = T * Y(k) ------(1)
Matlab returns back Y(k) from the FFT() function when given a sequence of
numbers.
Hence to scale Y(k) and obtain the sampled version of X(w), we need to multiply
Y(k) by T as per equation (1) above.
The example from Matlab help above was using one second for the duration of the
data and it sampled the data at a sampling frequency such that T = 1/L .
That is why the code divided by L.
It was actually multiplying by the sampling period T
which happened to be 1/L in this one example.
I think it would have been clearer to write
Y = fft(y,NFFT)*T
But in the example above it is same thing.
So, as a rule, if you know the signal is being sampled at a frequency larger
than twice the largest frequency embedded in the signal, then multiply the DFT
you obtain from Matlab fft() function by the sampling period.
Since in practice most of time we make sure we sample at
large enough frequency, then normal practice is to multiply the fft() result by
the sampling period.
The proof of equation (1) requires a bit of algebra, I scanned 2 pages from a
book I have on this (DSP by Liu and Peled) pages 38 and 39. Here it is
And look at question 1-93 in the above.
plot(f,2*abs(Y(1:NFFT/2+1)))
Why does it multiply by 2?
First, they should not have multiplied by 2 the Y(1) entry (i.e. the DC entry, ie. Frequency 0). The correct way would have been to do 2*Y, followed by Y(1)=Y(1)/2 to adjust the Y at DC back to what it was before plotting it.
But the reason they multiply Y by 2 is just normalization.
Matlab fft() returns Y(k) in the range 0 to fs (with a 2*pi factor there), and not from –fs/2..fs/2.
So, the code was only using Y(k) from 0..fs/2, i.e. first half of the Y(k). Hence it was throwing away the part from fs/2..fs. Since X(k) is symmetric, it was throwing away exactly half of the sequence Y(k) content.
To keep the ‘things’ equal (to keep the ‘energy’ if you well, and I am being loose in the terms I am using here) the same as before throwing half of Y away, we need to multiply by 2 again to ‘adjust’ it. (Think of an area, which we cut in half, then to keep the area the same, we multiply the remaining area back by 2.
The code above was plotting only half of Y(k). You do not need to multiply by 2 if you are plotting the whole Y(k). Normally what is done is to use fftshift (which make Y(k) centered around 0 instead of how Matlab returns it..
