Реализация средней мощности Matlab в Octave?

Question

Люди,

Matlab 2007b (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]);

Я хочу воспроизвести такого рода функциональность в Octave. Функция "pwelch" кажется вполне возможной.А именно:

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

Теперь я подумай spectra имеет значения y для PSD, а freq - значения x .Итак, я смог найти выборки по частоте, которые находятся между "нижним" и "верхним" и ..э-э, усреднить соответствующие значения в спектрах?Я довольно смутно представляю себе это.

Более того, значения в "freq" не обязательно соответствуют моим желаемым верхним и нижним значениям, и я не уверен, что с этим делать.Что делать, если нижняя или верхняя часть попадает прямо в середину ячейки wide 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);