Calculating energy

Question

Raf K 2012-3-20

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/32894-calculating-energy

回答： Mouna Bouhachem 2021-9-15

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 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

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

Answer 1

Dr. Seis 2012-3-21

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/32894-calculating-energy#answer_41619

在 MATLAB Online 中打开

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 个评论
显示 4更早的评论隐藏 4更早的评论

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 ".*"?

请先登录，再进行评论。

Answer 2

Wayne King 2012-3-20

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/32894-calculating-energy#answer_41466

在 MATLAB Online 中打开

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.

请先登录，再进行评论。

Answer 3

Raf K 2012-3-20

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/32894-calculating-energy#answer_41470

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