Why is FFT result divided by NFFT instead of the root of NFFT?

Hi everyone,

if

yfft = fft(y)

According to parvela's theorem, the equation below must be achieved.

sum(y.^2) == sum(abs(yfft).^2/length(yfft))

On both sides of the equation, actually, the power values are summed, right?.

So why is yfft divided by length of yfft to find the amplitude spectrum in the MATLAB FFT example? Shouldn't it be divided by the root of yfft length to ensure the equation?

Do I think wrong? If you could enlighten me on this matter, I would appreciate it.

Answer 1

[Paul]

Hi Miktat,

The energy spectrum and the amplitude spectrum are two different things. Using the conventions of fft, the output yfft is divided by numel(y) so as to match the amplitude of the complex sinusoids that make up the underlying signal.

For example, define a finite duration signal as the sum of two complex sinusoids

n = 0:99;
y = 1.5*exp(1j*2*pi*.3*n) + 3*exp(1j*2*pi*.7*n);

Compute the FFT

Y = fft(y);

Plot the result in the frequency domain, dividing by numel(Y)

N = numel(y);
figure
stem((0:N-1)/N,abs(Y)/N)
ylabel('abs(Y)')

We see that dividing the FFT, as defined in fft, by N yields the amplitude spectrum in the frequency domain with peaks at the same amplitudes as the constituent complex sinusoids.

You are, of course, correct that we can compute the energy in the signal in either domain

format long e
[sum(y.*conj(y)) sum(Y.*conj(Y))/N]

ans = 1×2

1.0e+00 * 1.125000000000000e+03 1.125000000000000e+03

We can plot the energy spectrum if that's what's desired

figure
stem((0:N-1)/N,Y.*conj(Y)/N)
ylabel('Y*conj(Y)/N')

with the understanding that the spectrum is to be summed to get the energy in the signal.

If we zero pad the FFT, then, as @dpb pointed out, we still divide by the length of the time domain signal to get the proper amplitude (not the length of the FFT)

Y = fft(y,1024);
Nfft = numel(Y);
figure
stem((0:Nfft-1)/Nfft,abs(Y)/N,'.')
ylabel('abs(Y)')

The peaks at 0.3 and 0.7 are still very close to the amplitudes of the complex sinusoids; the rest of the non-zero elements are filling in the true amplitude spectrum of the windowed signal.

However, the energy computation uses the length of the fft, i.e., Nfft

[sum(y.*conj(y)) sum(Y.*conj(Y))/Nfft]

ans = 1×2

1.0e+00 * 1.125000000000000e+03 1.124999999999999e+03

So you'd have to scale by Nfft if wanting to plot the energy spectrum after zero padding.

[Miktat Vural]

Thanks for your explanation. I guess I confused energy spectrum with amplitude spectrum.

Answer 2

[dpb]

Well, let's try it and see...

Fs = 1000;            % Sampling frequency
T = 1/Fs;             % Sampling period
L = 1500;             % Length of signal
t = (0:L-1)*T;
S = 0.7*sin(2*pi*50*t) + sin(2*pi*120*t);
X = S + 2*randn(size(t));
Y = fft(X);
sum(X.^2)==sum(abs(Y.^2))/L

ans = logical
   0

sum(X.^2)-sum(abs(Y.^2))/L

ans = 9.0949e-13

shows it is equal to almost the double-precision.

NOTA BENE: Remember that if pad the signal so that NFFT is >L, the divisor is still L, NOT NFFT since there's no power in the zero-padded elements.

[Miktat Vural]

Thanks for your explanation.