sigal average power estimation

Hi !

I've created a snippet to test my insight of power the density spectra computed from the fft matlab function. Unfortunately, I have a factor 4 between average power computed from the signal in the time domain and the average power estimated from the FFT. May I ask you some hints ? Regards. Here my snippet :

N=1024;
Fs=1e3;
Ts=1/Fs;
timeSlot = N*Ts;
tVect = (0:N-1)*Ts;
F_axis = (0:N-1)*Fs/N;
F_axis_Single = F_axis(1:length(F_axis)/2);

A2=2;A=sqrt(A2);
% number of harmonics
nH = 10;
% build a set of nH random frequencies ranged from 0 to Fs/2
freqIDs = sort(round((N/2)*rand(1,nH )));
setOfFreqs =  F_axis_Single(freqIDs);
% buid a set of harmonics of frequencies defined above
noiseHarmos = A*sin(2*pi*tVect'*setOfFreqs);
noiseSig = sum(noiseHarmos,2);

% figure;plot(tVect,sigs);

fftSig_Db = fft(noiseSig); 
fftSig_Sg = fftSig_Db(1:N/2);   % FFT single sided coeff 
fftSig_Sg(2:end) = 2*fftSig_Sg(2:end);

% norm of FFT
absFftSig_Sg = abs(fftSig_Sg); % /N  to get FFT in V
absFftSig_Sg_V = abs(fftSig_Sg)/N; % /N  to get FFT in V
% DSP
dspSigSg = absFftSig_Sg.^2/N;

% Average powa from time domain
powFromSig = mean(noiseSig.^2)
%  Average powa from frequency domain
powFromDsp = sum(dspSigSg)/N

Answer 1

[Suneesh]

I guess what you are trying to verify is the validity of Parseval's theorem. Since FFT calculates the Discrete Fourier Transform (but not an (infinitely long) discrete-time Fourier transform ) a scaling is required. Please refer to the following links to clear this up:

http://www.mathworks.com/help/signal/ug/psd-estimate-using-fft.html

http://www.mathworks.us/matlabcentral/answers/98132-how-does-the-output-of-the-psd-function-compare-with-the-output-of-the-fft-function-in-matlab-7-11

[Sylvain Rousseau]

Dear Suneesh,

you're right, I'd like to verify Parseval Theorem. I went through the links you suggested, but I still do not understand properly the way the PSD is computed from the matlab FFT.

Parseval tells that if you perform the integral of the squared modules of the Fourier coefficients, then divide this area with the time windows T during which you oberved the signal, then you retrieve the signal average power.

If I first consider the discrete signal x made of T/Ts = N samples, then the average power is simply : avPowFromSig = sum(x^2)*Ts/T = sum(x^2)/N

In the frequency space, I will use xdft = fft(x) to retrieve the Fourier coefs. To get psd, first, I take squared modules of the dft coefs : sqModOfDft = abs(xdft ).^2 But I have to remove the scaling factor embeded into the xdft vector to get the actual square of the dft coefs magnitude : sqModOfDft = sqModOfDft/N^2 Now, the integral of sqModOfDft in the frequency domain is : powFromPsd = Fs/N*sum(sqModOfDft ) To get the average power I just have to divide by the window time duration : avPowFromPsd = powFromPsd/T Unfortunately, this reasoning leads to a bad result.

May I as you some help ?

Regards