- home

Sept. 17, 2009 page compiled on November 8, 2015 at 9:21pm

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)|')

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)|')

Looking at this line Y = fft(y,NFFT)/L;

Why the code divides the FFT result by L?

And this line plot(f,2*abs(Y(1:NFFT/2+1)))

This scaling factor comes from Nyquist sampling theory.

Assuming the analog signal is with corresponding Fourier transform .

Sampling with some sampling frequency results in a sequence of numbers: the sampled sequence where is the sampling period and is the sequence number.

Let the sequence of complex numbers be the Fourier transform of the sampled sequence, i.e. the DFT of

If the sampling frequency is larger than twice the largest frequency in the signal then the magnitude of will be proportional to the magnitude of .

The proportionality factor turns out to be the sampling period .

This means on can write

Matlab returns back from the FFT() function when given a sequence of numbers.

Hence to scale and obtain the sampled version of then is multiplied by as per equation 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 .

That is why the code divided by .

It was actually multiplying by the sampling period which happened to be in this example.

I think it would have been clearer to write Y = fft(y,NFFT)*T

But in the example above it is same thing.

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 ﬀt() 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 ﬀt() result by the sampling period.

The proof of above equation 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 looking at equation 1-93 in the above.

Why does it multiply by ?

First, they should not have multiplied by 2 the entry (i.e. the DC entry at frequency 0). The correct way would have been to do , followed by to adjust the at DC back to what it was before plotting it.

But the reason they multiplied by is just normalization.

Matlab ﬀt() returns Y(k) in the range to (with a factor), and not from .

So, the code was only using from , i.e. ﬁrst half of the .

Hence it was throwing away the part from .

Since is symmetric, it was throwing away exactly half of the sequence .

To keep the things equal (or to keep the energy the same, if you well, and I am being loose in the terms used here) as before throwing half of away, we need to multiply by again to adjust it. (Similar to an area, which we cut in half, then to keep the area the same, we multiply the remaining area back by ).

The code above was plotting only half of . You do not need to multiply by if you are plotting the whole . Normally what is done is to use fftshift (which makes centered around 0 instead of how Matlab returns it)