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.

Answer 1

[Dr. Seis]

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.

[Raf K]

Thanks a lot - what u did makes sense however I got this error:

" Error using *

Inner matrix dimensions must agree./"

I believe the problem is in this line

y_squared = (y/Fs)*conj(y/Fs)

but isn't y_squared in your equations proportional to p values i got from the spectrogram command.

[Dr. Seis]

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]

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]

why would you normalize fft by fs?

[Dr. Seis]

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]

Why do you square y with its complex conjugate and not just itself? Shouldn't we multiplying element wise with ".*"?

Answer 2

[Wayne King]

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!

[Dr. Seis]

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]

Exactly, I was assuming the OP was really interested in average power.

Answer 3

[Raf K]

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 ...

[Wayne King]

Then you can certainly just sum the powers at at the respective frequencies.

[Raf K]

Thank you

[Dr. Seis]

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).

Answer 4

[Raf K]

Thanks a lot

Answer 5

[Mouna Bouhachem]

Thanks Very Much!