Implémentation de la puissance moyenne de Matlab dans Octave?

Question

Les gens,

Matlab 2007b (7.5.0) possède une fonction de puissance moyenne. Voir ici :

& "; La méthode avgpower utilise une approximation rectangle de l'intégrale à calculer la puissance moyenne du signal en utilisant les données PSD stockées dans le objet.

& "; La méthode avgpower renvoie la puissance moyenne du signal qui est l'aire sous la courbe PSD. "

Exemple d'appel:

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

Je cherche à reproduire ce type de fonctionnalité dans Octave. le fonction " pwelch " semble être une possibilité. À savoir:

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

Maintenant, je pense que le spectre a les valeurs y du PSD et que freq a le x valeurs. Donc, je pourrais trouver les échantillons dans la fréquence qui se situent entre & "; Inférieur &"; et " supérieur " et .. euh, faire la moyenne des valeurs correspondantes dans le spectre? Je suis assez flou à ce sujet.

De plus, les valeurs dans " freq " ne correspond pas nécessairement à mon désiré haut et bas, et je ne sais pas quoi faire à ce sujet. Que se passe-t-il si la partie inférieure ou supérieure tombe en plein milieu d'une large bande de fréquences? Par exemple, est-ce que je prends une demi-case (c'est-à-dire une interpolation linéaire)?

Il pourrait également être possible d’obtenir une valeur unique à partir d’une sorte de FFT au lieu d'utiliser pwelch.

Des suggestions?

Answer 1

Apparemment, je me parle à moi-même, mais voici un code Octave proposé pour ceux qui errent de cette façon.

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