¿Implementando la potencia promedio de Matlab en Octave?

Question

Amigos,

Matlab 2007b (7.5.0) tiene una función avgpower.Ver aquí:

"El método AVGPower utiliza una aproximación de rectángulo a la integral para calcular la potencia promedio de la señal utilizando los datos de PSD almacenados en el objeto.

"El método AVGPower devuelve la potencia promedio de la señal que es el área debajo de la curva PSD".

Ejemplo de invocación:

    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]);

Estoy buscando replicar este tipo de funcionalidad en Octave.El función "pwelch" parece una posibilidad.Esto es:

    ...
    sig = amplitude * sin(2*pi*frequency*t);
    pwelch('R12+');
    [spectra, freq]=pwelch(signal, [], [], [], Fs, plot_type='dB');

Ahora yo pensar spectra tiene los valores y de la PSD y freq tiene la x valores.Por lo tanto, pude encontrar las muestras en frecuencia que caen entre "inferior" y "arriba" y ..Eh, ¿promediar los valores correspondientes en los espectros?Estoy bastante confuso en esto.

Por otra parte, los valores en "freq" no necesariamente corresponden a mi deseado arriba y abajo, y no estoy seguro de qué hacer al respecto.¿Qué pasa si el inferior o el superior caen justo en el medio del amplio rango de frecuencias?Por ejemplo, ¿tomo medio contenedor (es decir,interpolar linealmente)?

También podría ser posible obtener un único valor de algún tipo de FFT en lugar de usar pwelch.

¿Sugerencias?

Answer 1

Aparentemente estoy hablando conmigo mismo, pero aquí hay un código de octava propuesto para aquellos que deambulan por este camino.

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);