实施Matlab的在八度avgpower？

Question

民间，

Matlab的2007年b（7.5.0）具有avgpower功能。请参见这里：

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

我期待复制这种在八度的功能。该 功能“pwelch”似乎是一个可能性。即：

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

现在我认为光谱具有PSD的y值和频率具有在x 值。所以，我能找到的频率落在之间的“低”的样本 和“上”和..呃，平均在光谱相应的值？ 我对这个很模糊。

此外，在“频率”的值不一定对应于我的 所需的上，下，和我不知道该怎么做了一番。如果上或下跌倒就在广泛频率斌中间？例如，是否取一半的仓（即线性插入）？

这也可能是可能从某种FFT的得到一个单一的价值 代替使用pwelch。

建议？

Answer 1

显然，我对自己说，但这里是为那些谁徘徊这样一些建议倍频代码。

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