sigal average power estimation

Question

Sylvain Rousseau 2013-12-9

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/109222-sigal-average-power-estimation

编辑： Sylvain Rousseau 2013-12-10

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

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Suneesh 2013-12-9

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/109222-sigal-average-power-estimation#answer_117888

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

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

Sylvain Rousseau 2013-12-10

编辑：Sylvain Rousseau 2013-12-10

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

请先登录，再进行评论。