R – Implementing Matlab’s avgpower in Octave

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?