Re: avgpower in Octave?

[Dan F]

Thanks for your comments everyone. For others who come this way, I
attach below some code that seems to agree with Matlab, more or less,
under equivalent circumstances (welch, same window size, etc.).

Comments welcome.

Dan

function[avgp] = oavgpower(signal, sampling_freq, lowfreq, highfreq, window)

[spectra, freq]=pwelch(signal, window, [], [], sampling_freq);

idx1=max(find(freq <= lowfreq));
idx2=min(find(freq >= highfreq));

% Index and actual frequency of lower and upper bins
%idx1
%freq(idx1)
%idx2
%freq(idx2)

% 0: don't include the last bin
width = [diff(freq); 0];

pvec = width(idx1:idx2).*spectra(idx1:idx2);
avgp = sum(pvec);


On Thu, Mar 19, 2009 at 12:13 PM, Francesco Potorti`
<address@hidden> wrote:
>>Matlab 2007b (7.5.0) has an avgpower function. From
>>http://www.mathworks.com/products/signal/demos.html?file=/products/demos/shipping/signal/spectralanalysisobjsdemo.html:
>>
>>"The avgpower method uses a rectangle approximation to the integral to
>>calculate the signal's average power using the PSD data stored in the
>>object.
>>
>>"The avgpower method returns the average power of the signal which is
>>the area under the PSD curve."
>>
>>Example invocation:
>>
>>    numSamples = 10000
>>    frequency = 20
>>    amplitude = 10
>>    Fs = 1000
>>    t = [0:1/numSamples:1];
>>    sig = amplitude * sin(2*pi*frequency*t);
>>    h = spectrum.periodogram('rectangular');
>>    hopts = psdopts(h, signal);
>>    set(hopts,'Fs',Fs);
>>    p = psd(h,signal,hopts);
>>    lower = 12
>>    upper = 30
>>    beta_power = p.avgpower([lower upper]);
>>
>>I am looking to replicate this sort of functionality in Octave. The
>>function "pwelch" seems like a possibility.
>
> You use pwelch to compute the PSD.  It is compatible with Matlab's.
>
>>                                            To wit:
>>
>>    ...
>>    sig = amplitude * sin(2*pi*frequency*t);
>>    pwelch('R12+');
>>    [spectra, freq]=pwelch(signal, [], [], [], Fs, plot_type='dB');
>>
>>Now I *think* spectra has the y values of the PSD and freq has the x
>>values. So, I could find the samples in freq that fall between "lower"
>>and "upper" and .. er, average the corresponding values in spectra?
>
> You should compute the integral, so I think that you can average them
> and multiply by upper-lower, if upper and lower are frequencies as they
> appear to be.
>
>>Moreover, the values in "spectra" don't necessarily correspond to my
>>desired upper and lower, and I'm not sure what to do about that.
>
> I do not follow you here.
>
>>It might also be possible to get a single value from some sort of FFT
>>instead of using pwelch.
>
> You want to compute the power is a given range of frequencies, if I
> understand correctly.  So you must first estimate the PSD, which is why
> you need pwelch.
>
> --
> Francesco Potortì (ricercatore)        Voice: +39 050 315 3058 (op.2111)
> ISTI - Area della ricerca CNR          Fax:   +39 050 315 2040
> via G. Moruzzi 1, I-56124 Pisa         Email: address@hidden
> (entrance 20, 1st floor, room C71)     Web:   http://fly.isti.cnr.it/