Реализация способности Маглаба в Октаве?

Question

Люди,

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

Теперь я думаю спектров имеют значение Y СДП и частота имеет Й- значения. Итак, я мог найти образцы по частоте, которые находятся между «нижними» и «верхними» и ... эр, средними значениями в спектрах? Я довольно нечеткий на этом.

Кроме того, значения в «freq» не обязательно соответствуют моим желаемым верхним и нижним, и я не уверен, что с этим делать. Что делать, если нижний или верхний падает прямо в середине широкого диапазона частот? Например, беру ли я половину бина (т. Е. Линейно интерполировать)?

Возможно также получить одно значение из своего рода 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);