Calculating energy
显示 更早的评论
Hi,
I would appreciate your help and suggestions. I know this is a silly question for most of the readers, but this has been bugging me for sometime. I generated a PSD of a signal using spectrogram command (matlab signal processing toolbox). I used a 250 ms window. to calculate the energy for a frequency band in a certain time window, can i simply add up the power values of the individual frequencies or should i multiply power values * frequency then add up the results .. or ... for example: WINDOW1: (freq 1 Hz power 20) (freq 2Hz power 20) (freq 3 Hz power 15) is the energy of 1-3 Hz band in WINDOW1 is 20 + 20 + 15 or is it 1 X 20 + 2 X 20 + 3 X 15 ... Thanks a lot.
采纳的回答
Dr. Seis
2012-3-21
In order to obtain energy, do not even worry about the power. Just do this:
Using "*y*" and "*f*" from the output of spectrogram -
N = 250; % window length
Fs = 500; % sample frequency
df = Fs/N; % frequency increment
y_squared = (y/Fs)*conj(y/Fs); % Fs is used to normalize the FFT amplitudes
energy_10Hz_to_90Hz = 2*sum(y_squared( f>=10 & f<=90,: ))*df;
If at any point the range of frequencies you want to calculate the energy over includes either 0Hz or the Nyquist frequency, then you will not be multiplying those amplitudes by 2. Otherwise the above should work. You need to multiply the sum of squared amplitudes by the frequency increment (df) as I do above because you are essentially performing discrete integration across the interval between 10Hz and 90Hz.
6 个评论
Raf K
2012-3-22
Dr. Seis
2012-3-22
Sorry about that...
y_squared = (y/Fs).*conj(y/Fs);
Yes, they are most certainly proportional. However, the power values in "p" are the result of multiplying the squared-amplitude values (similar to mine) with another constant "k" (see documentation for spectrogram). When I tried to use the power amplitudes to finagle an estimate for energy (in my tests) my amplitudes were off by some constant factor, which is why I decided not to even bother with the power values anymore and just stick to what I know.
annnne
2013-3-15
energy_10Hz_to_90Hz = 2*sum(y_squared( f>=10 & f<=90,: ))*df; must be change?? or fixed? if i just use delta,theta alpha and beta?? what the value for my energy?
Ray Lee
2014-12-5
why would you normalize fft by fs?
Dr. Seis
2014-12-29
The Fourier transform of a continuous time signal is an integral where the variable of integration is dt. Since we do not record data continuously, the continuous Fourier transform turns into a discrete Fourier transform (DFT) and our variable of integration ( dt ) ends up as a scalar multiplier (delta_t = 1/Fs) whose product with the DFT (which in Matlab is just a summation of the product of the time domain amplitudes with a time-&frequency-dependent exponential term - i.e., no normalization) approximates the continuous Fourier transform over a finite interval. The more samples you have over a time interval the more closely you will approximate a continuous result.
If you want the amplitude values in the frequency domain to have a physical meaning, then you multiply the DFT by delta_t (i.e., normalize by Fs). If all you want is a display to look at, then normalizing by Fs doesn't matter because you are mostly interested in relative amplitude changes versus frequency.
Benjamin Shafer
2017-8-4
Why do you square y with its complex conjugate and not just itself? Shouldn't we multiplying element wise with ".*"?
更多回答（4 个）
Wayne King
2012-3-20
Hi Raf, you should do this:
should i multiply power values * frequency then add up the result
BUT you want to multiply by the width in frequency, or your delta f. What is the different in Hz between adjacent DFT bins? You can use diff() on your frequency vector to get that info. Multiplying your power estimates by that delta f and summing the result approximates the integral under the PSD.
I don't know how long this vector is in samples (250 msec), but you can check yourself by just using the avgpower method on a single 250 msec sample
x = randn(1e3,1);
psdestx = psd(spectrum.periodogram,x,'Fs',1,'NFFT',length(x));
% power in frequency interval 1/4 to 1/2
pwr = avgpower(psdestx,[1/4 1/2]);
It's probably more meaningful to take the frequency interval width into account as you suggest in your second option.
For example:
t = 0:0.001:1-0.001;
x = 2*cos(2*pi*100*t);
psdestx = psd(spectrum.periodogram,x,'Fs',1e3,'NFFT',length(x));
pwr = avgpower(psdestx);
Answer is 4/2 as you would expect, now just look in the range around 100 Hz.
pwr = avgpower(psdestx,[95 105])
Still 2, because all the power is there!
2 个评论
Dr. Seis
2012-3-20
But that is average power... the energy in a cosine wave is dependent on how long the cosine wave signal is. The energy (in the time domain) is equal to the integral of the amplitudes squared... therefore the energy in your signal will increase more and more as you make your "t" longer and longer. That means that the frequency domain amplitude at 100Hz (in your example) gets bigger and bigger... remember Parseval's theorem!
Wayne King
2012-3-20
Exactly, I was assuming the OP was really interested in average power.
Raf K
2012-3-20
0 个投票
3 个评论
Wayne King
2012-3-20
Then you can certainly just sum the powers at at the respective frequencies.
Raf K
2012-3-20
Dr. Seis
2012-3-21
Energy is NOT the sum of power. Power has units of energy per second (in the time domain) and energy per Hertz (in the frequency domain).
类别
在 帮助中心 和 File Exchange 中查找有关 Parametric Spectral Estimation 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
