Die Umsetzung Matlab avgpower in Octave?

Question

Leute,

Matlab 2007b (7.5.0) eine avgpower Funktion. Siehe hier :

"The avgpower Verfahren verwendet ein Rechteck Annäherung an dem integralen Bestandteil Berechnen der Durchschnittsleistung des Signals in der PSD-Daten unter Verwendung der gespeicherten Objekt.

"Die avgpower Methode gibt die mittlere Leistung des Signals, das ist die Fläche unter der PSD-Kurve. "

Beispiel Aufruf:

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

Ich suche diese Art von Funktionalität in Octave zu replizieren. Das Funktion „pwelch“ scheint wie eine Möglichkeit. Nämlich:

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

Jetzt I denken Spektren, die die y-Werte der PSD hat und Freq hat die x Werte. So konnte ich die Proben in Freq finden, die zwischen „unten“ fallen und „oben“ und .. äh, den Mittelwert der entsprechenden Werte in den Spektren? Ich bin mir ziemlich unscharf zu diesem Thema.

Darüber hinaus sind die Werte in „Freq“ entsprechen nicht notwendigerweise meiner gewünscht obere und untere, und ich bin nicht sicher, was das zu tun. Was, wenn die untere oder obere fällt genau in der Mitte des breiten Freq ist? Zum Beispiel nehme ich eine halbes ist (das heißt linear interpoliert)?

Es könnte auch möglich sein, einen einzelnen Wert aus einer Art von FFT zu erhalten anstelle der Verwendung pwelch.

Verbesserungsvorschläge?

Answer 1

Offenbar sprechen ich mir selbst, aber hier ist einige vorgeschlagene Octave Code für diejenigen, die auf diese Weise wandern.

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