Implementazione in Matlab del avgpower in Ottava?

Question

La gente,

Matlab 2007b (7.5.0) ha un avgpower funzione.Vedere qui:

"Il avgpower metodo utilizza un rettangolo approssimazione dell'integrale di calcolare il segnale di media potenza con il PSD dati memorizzati nel oggetto.

"Il avgpower metodo restituisce la potenza media del segnale l'area sotto la curva di densità spettrale."

Esempio di invocazione:

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

Sto cercando di replicare questo tipo di funzionalità in Ottava.Il funzione "pwelch" sembra una possibilità.Vale a dire:

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

Ora Ho pensare spectra ha i valori di y del PSD e freq è la x valori.Così, ho potuto trovare i campioni in freq che rientrano tra "basso" e "superiore" e ..er, in media i valori corrispondenti in spectra?Io sono abbastanza fuzzy su questo.

Inoltre, i valori in "freq" non necessariamente corrisponde al mio desiderato superiore e inferiore, e io non sono sicuro di cosa fare al riguardo.Se è per quello, inferiore o superiore, cade proprio nel mezzo di una vasta freq bin?Per esempio, posso prendere la metà di un bin (es.interpolare linearmente)?

Potrebbe anche essere possibile ottenere un singolo valore di una sorta di FFT invece di utilizzare pwelch.

Suggerimenti?

Answer 1

A quanto pare, sto parlando a me stesso, ma qui è stato proposto di Ottava codice per coloro che vagano in questo modo.

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