MATLAB Answers

Calculating energy

497 views (last 30 days)
Raf K
Raf K on 20 Mar 2012
Commented: Benjamin Shafer on 4 Aug 2017
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.

  0 Comments

Sign in to comment.

Accepted Answer

Dr. Seis
Dr. Seis on 21 Mar 2012
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 Comments

Show 3 older comments
Ray Lee
Ray Lee on 5 Dec 2014
why would you normalize fft by fs?
Dr. Seis
Dr. Seis on 29 Dec 2014
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
Benjamin Shafer on 4 Aug 2017
Why do you square y with its complex conjugate and not just itself? Shouldn't we multiplying element wise with ".*"?

Sign in to comment.

More Answers (3)

Wayne King
Wayne King on 20 Mar 2012
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 Comments

Dr. Seis
Dr. Seis on 20 Mar 2012
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
Wayne King on 20 Mar 2012
Exactly, I was assuming the OP was really interested in average power.

Sign in to comment.


Raf K
Raf K on 20 Mar 2012
Hi thank you for your input. Actually I'm interested of the Energy (sum) not in the average power. Here is my example:
x contains data of 114 second of a signal sampled at 500 Hz.
my code:
t=0:1:114000; F=1:1:90; [y,f,t,p] = spectrogram(x,250,125,F,500,'yaxis');
Now I have a p matrix ( rows frequencies (1:90) and columns are time windows as generated by the spectrogram command above) and the p values reflect the power of certain frequency in certain window.
My question is simply: to calculate the ENERGY of the frequency band 10-90 in the first window ... should i simply sum p for frequencies 10-90 or should i multiply each frequency by its corresponding power then sum it up
Thanks a lot ...

  3 Comments

Wayne King
Wayne King on 20 Mar 2012
Then you can certainly just sum the powers at at the respective frequencies.
Raf K
Raf K on 20 Mar 2012
Thank you
Dr. Seis
Dr. Seis on 21 Mar 2012
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).

Sign in to comment.


Raf K
Raf K on 24 Mar 2012
Thanks a lot

  0 Comments

Sign in to comment.

Sign in to answer this question.