2009-03-19 11 views
1

フォーク、MatlabのavgpowerをOctaveで実装していますか?

Matlab 2007b(7.5.0)にはavgpower関数があります。 hereを参照してください

「avgpower方法は オブジェクトに格納されているPSDデータを用いて、信号の平均電力を計算 と一体に矩形近似を使用

」avgpower方法は、信号の平均電力を返しますこれ PSD曲線下面積である「

例の呼び出し:。

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

私が探していますこの種の機能をOctaveで再現しています。 関数 "pwelch"は可能性のようです。ウィットに:

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

は、今私はスペクトルはPSDのyの値を持っており、FREQは、x 値を持っていると思う。だから、私は "下" と "上"との間にある周波数のサンプルを見つけることができます。そして、平均、スペクトルの対応する値ですか? 私はこれをかなりあいまいにしています。

さらに、「freq」の値は、私の 希望の上限と下限に対応しているとは限りません。そのために何をすべきかわかりません。下段または上段が広い周波数のビンの真ん中に落ちたらどうなりますか?たとえば、半分のビンを取るか(すなわち、直線的に補間する)か?

pwelchを使用する代わりに、ある種のFFT から単一の値を取得することもできます。

提案?

答えて

2

明らかに私は自分自身に話していますが、ここでこのようにさまよっている人のためにオクターブコードを提案しています。

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

私はPythonで作ったが、それでも良い例がある – PatriceG

関連する問題