Thread Subject:

Why does FFT(x) must be divided by number of samples, i.e. length(x)?

The following code is extracted from Matlab help file and has been modified as follows:

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
x = 10*sin(2*pi*50*t)
y = x
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)|')


May I know why do we need to divide fft(y,NFFT) by L? If we don't do that, then the peak at 50Hz is not equal to 10, but equal to 10000.

Thank you very much

onemilimeter wrote:
...

> May I know why do we need to divide fft(y,NFFT) by L? If we don't
> dot hat, then the peak at 50Hz is not equal to 10, but equal to 10000....

Simple answer is "because"... :)

Actually, it's a result of the integration process and a choice TMW made
of what to return as the value of the FFT function as implemented.
Other implementations do choose to normalize; TMW didn't...

--

On Dec 8, 6:05 pm, "onemilimeter " <on...@example.com> wrote:
> The following code is extracted from Matlab help file and has been modified as follows:
>
> 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
> x = 10*sin(2*pi*50*t)
> y = x
> 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)|')
>
> May I know why do we need to divide fft(y,NFFT) by L? If we don't do that, then the peak at 50Hz is not equal to 10, but equal to 10000.
>
> Thank you very much

If y(t) is periodic with period T and fundamental frequency
f0 = 1/T, the Fourier series coefficients (m=1,2,...) are given
by

Ym = (1/T)*INTEGRAL(0,T){ dt exp(-j*2*pi*f(m)*t) * y(t) ) }

where f(m) = (m-1)*f0 = (m-1)/T.

Since the integrand is periodic with period T, the approximation
to the integral represented by the sum of N-1 interior rectangles
of width dt and 2 end rectangles of width dt/2 centered at times
t(n) = (n-1)*dt , n = 1:N+1 can be reduced to the sum of N
rectangles of width dt centered at times t(n) = (n-1)*dt, n=1:N.

Therefore

Ym = (1/T)*SUM(n=1:N){ dt exp(-j*2*pi*f(m)*t(n)) * y(n) ) }

The factor dt/T = 1/N is missing in the MATLAB implementation.

Hope that helps.

Greg

On 9 Des, 00:05, "onemilimeter " <on...@example.com> wrote:

> May I know why do we need to divide fft(y,NFFT) by L?

First of all, it's the convention just about everybody use.
Second, the reason why everybody use this particular
convention is that it saves a few computations. If you
implement a filter, say, y[n] = x[n] (*) h[n] in frequency
domain as

y = ifft(fft(x).*fft(h));

there are obviously two forward DFTs and one inverse DFT.
With this convention one scales only one transform, the IFFT,
instead of scaling three places. It might not look like much
of a saving these days, but once upon a time it was significant
enough that the idea made it to the textbooks, where it first
of all became the standard method, and second has stuck ever
since.

Rune

Rune Allnor <allnor@tele.ntnu.no> wrote in message <cf337b2e-c731-4ae1-bd09-37fff3b63908@u7g2000yqm.googlegroups.com>...
> On 9 Des, 00:05, "onemilimeter " <on...@example.com> wrote:
>
> > May I know why do we need to divide fft(y,NFFT) by L?
>
> First of all, it's the convention just about everybody use.
> Second, the reason why everybody use this particular
> convention is that it saves a few computations. If you
> implement a filter, say, y[n] = x[n] (*) h[n] in frequency
> domain as
>
> y = ifft(fft(x).*fft(h));
>
> there are obviously two forward DFTs and one inverse DFT.
> With this convention one scales only one transform, the IFFT,
> instead of scaling three places. It might not look like much
> of a saving these days, but once upon a time it was significant
> enough that the idea made it to the textbooks, where it first
> of all became the standard method, and second has stuck ever
> since.
>
> Rune


Hello all, interesting thread,

I'd like to know why in

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

you multiply the absolute value by two?

thanks

Tom

>
> Hello all, interesting thread,
>
> I'd like to know why in
>
> plot(f,2*abs(Y(1:NFFT/2+1)))
>
> you multiply the absolute value by two?
>
> thanks
>
> Tom

I would like to know why there is the multiply by two also ...

Thx,
Forrest

"Forrest " <forrestlundstrom@conmed.com> wrote in message
news:je7794$st3$1@newscl01ah.mathworks.com...
>>
>> Hello all, interesting thread,
>>
>> I'd like to know why in plot(f,2*abs(Y(1:NFFT/2+1))) you multiply the
>> absolute value by two?
>>
>> thanks Tom
>
> I would like to know why there is the multiply by two also ...

For context, this is from the second example in the documentation for FFT:

http://www.mathworks.com/help/techdoc/ref/fft.html

--
Steve Lord
slord@mathworks.com
To contact Technical Support use the Contact Us link on
http://www.mathworks.com

"Forrest" wrote in message <je7794$st3$1@newscl01ah.mathworks.com>...
>
>
> I would like to know why there is the multiply by two also ...
=================

Because

A*exp(-j*(2*pi*f+phi)) +A*exp(j*(2*pi*f+phi)) = 2*A cos(2*pi*f+phi)

In other words, if you want to convert FFT coefficients, which are coefficients of exponentials, to cosine coefficients, they are related by a factor of 2.

On Jan 6, 8:41 am, "Matt J " <mattjacREM...@THISieee.spam> wrote:
> "Forrest" wrote in message <je7794$st...@newscl01ah.mathworks.com>...
>
> > I would like to know why there is the multiply by two also ...
>
>================
>
> Because
>
> A*exp(-j*(2*pi*f+phi)) +A*exp(j*(2*pi*f+phi)) = 2*A cos(2*pi*f+phi)
>
> In other words, if you want to convert FFT coefficients, which are coefficients of exponentials, to cosine coefficients, they are related by a factor of 2.

In addition, note that when f=0 and f=Fs/2, the amplitude is 1, not 2.

Hope this helps.

Greg

1. could anybody please explain to me about the length of signal L (it may sound like a stupid question, but I'm brandnew with FFT) ? or show me where I can get the information about it?

2. In the following code (I removed the noise), if I change L (even L=5000), I couldn't get the real amplitudes.

3. In this code, if I change NFFT=L, I will get the exact values of amplitude, no matter what L is.

Thanks very much in advance.

Fs = 1000; % Sampling frequency
T = 1/Fs; % Sample time
L = 5000; % Length of signal
t = (0:L-1)*T; % Time vector
x = 0.7*sin(2*pi*50*t) + sin(2*pi*120*t);
NFFT = 2^nextpow2(L); % Next power of 2 from length of y
X = fft(x,NFFT)/L;
f = Fs/2*linspace(0,1,NFFT/2+1);
plot(f,2*abs(X(1:NFFT/2+1)))

"Mai Le" <lemai289@yahoo.com.vn> wrote in message <k031km$pka$1@newscl01ah.mathworks.com>...
> 1. could anybody please explain to me about the length of signal L (it may sound like a stupid question, but I'm brandnew with FFT) ? or show me where I can get the information about it?
>
> 2. In the following code (I removed the noise), if I change L (even L=5000), I couldn't get the real amplitudes.
>
> 3. In this code, if I change NFFT=L, I will get the exact values of amplitude, no matter what L is.
===============

It's because your frequency sampling axis doesn't contain 50 or 120. The frequency sample spacing in your code is dF=1/NFFT/T
which doesn't divide evenly into 50 or 120 unless L=NFFT.

The bottom line is that to sample particular frequencies, like 50 or 120, you have to start the design by choosing your frequency sampling interval dF first, so that it divides evenly into those desired frequencies. Then, you have to design your T and NFFT based on dF to satisfy

T*NFFT=1/dF


I don't know why you think NFFT needs to be a power of 2, but if you insist that it must be, one way to do this is as follows:

dF=0.2;%Pre-select frequency sample spacing
NFFT = 2^nextpow2(5000); % Next power of 2 from length of y
T=1/NFFT/dF;

t = (0:NFFT-1)*T; % Time vector
x = 0.7*sin(2*pi*50*t) + sin(2*pi*120*t);

 X = fft(x)/NFFT;
 f = Fs/2*linspace(0,1,NFFT/2+1);
 plot(f,2*abs(X(1:NFFT/2+1)))