497 views (last 30 days)

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

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
on 4 Aug 2017

Sign in to comment.

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!

Dr. Seis
on 20 Mar 2012

Sign in to comment.

Raf K
on 20 Mar 2012

Wayne King
on 20 Mar 2012

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

Dr. Seis
on 21 Mar 2012

Sign in to comment.

Sign in to answer this question.

Opportunities for recent engineering grads.

Apply Today
## 0 Comments

Sign in to comment.