Implementing Matlab's avgpower in Octave?

Question

Folks,

Matlab 2007b (7.5.0) has an avgpower function. See here:

"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. 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? I'm pretty fuzzy on this.

Moreover, the values in "freq" don't necessarily correspond to my desired upper and lower, and I'm not sure what to do about that. What if the lower or upper falls right in the middle of wide freq bin? For example, do I take half a bin (i.e. linearly interpolate)?

It might also be possible to get a single value from some sort of FFT instead of using pwelch.

Suggestions?

Answer 1

Apparently I'm talking to myself, but here is some proposed Octave code for those who wander this way.

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

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

idx1=max(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);