avgpower de Implementação Matlab em Octave?

Question

Gente,

Matlab 2007b (7.5.0) tem uma função avgpower. Consulte aqui :

"O método avgpower usa uma aproximação retângulo para a integral para calcular a potência média do sinal usando os dados PSD armazenados na objeto.

"O método avgpower retorna a potência média do sinal que é a área sob a curva de PSD. "

Exemplo invocação:

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

Eu estou olhando para replicar esse tipo de funcionalidade em Octave. o função "pwelch" parece ser uma possibilidade. A saber:

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

Agora eu pensar espectros tem os valores de y do PSD e freq tem a x valores. Então, eu poderia encontrar as amostras em freq que se situam entre "baixar" e "superior" e .. er, a média dos valores correspondentes no espectro? Estou muito confuso sobre isso.

Além disso, os valores "freq" não corresponde necessariamente ao meu desejado superior e inferior, e eu não sei o que fazer sobre isso. E se o inferior ou superior cai bem no meio de bin freq ampla? Por exemplo, posso tomar metade de um bin (ou seja linearmente interpolate)?

Também pode ser possível obter um único valor de algum tipo de FFT em vez de usar pwelch.

Sugestões?

Answer 1

Aparentemente eu estou falando para mim mesmo, mas aqui é algum código Octave proposto para aqueles que vagueiam desta forma.

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