Как получить правильную фазу из алгоритма Герцеля Опции

[kumle]

Использую алгоритм Герцеля для выделения одной составляющей из сигнала, амплитуда вычисляется верно.
Понадобилось получить также фазу. Вычисляю ее так atan(Im/Re)
Для проверки прогоняю алгоритм через массив выборок синуса
for (i=0; i<N; i++)
{
S[i]=150*sin(2*3.1415926*50*i*0.002); //частота дискретизации 500 Гц, N=512 выборок
}
у которого начальная фаза 0, но когда вычисляю арктангенс, то получается -18 градусов ?? (хотя должен быть ноль по идее)

Причем если добавить фазу в массив:
for (i=0; i<N; i++)
{
S[i]=150*sin(2*3.1415926*50*i*0.002+(88*pi/180)); //добавил 88 градусов
}
то в результате вычисления фаза получается 70 градусов, что как раз составляет -18+88.
Отсюда вопрос что это за -18 ?

[kumle]

Герцель - Герцелю рознь.
Вот этот правильно считает http://www.dsplib.ru/content/goertzel/goertzel.html

Самое забавное что я использовал алгоритм именно с этого источника !

У алгоритма Герцеля, даже с плавающей запятой, при больших N (N>200) накапливается большая погрешность.
Так что -18 может быть и нормально.

У меня 500 выборок.

Если бы накапливалась погрешность тогда бы и после добавления начальной фазы в сигнал была бы неправильная дельта по фазе, но ведь у меня в сумме получается -18+70=88 градусов, вопрос в том откуда берется изначальная фаза -18 градусов, причем она зависит еще от частоты которую хотим выделить из сигнала, правда пока не понял какая зависимость

[bahurin]

если резулатат вашего алгоритма Герцеля отличается от fft то значит где-то ошибка. Результат должен совпадать. Теперь по существу. Вот матлаб код для проверки. Берем 500 отсчетов синусоиды s = sin(2*pi*50*t/Fs) на частоте 50 Гц и рассчитываем 50-ый отсчет FFT используя алгоритм герцеля. Получаем фазу -90 градусов, поскольку фаза 0 будет если взять s = сos(2*pi*50*t/Fs). Если посмотрите спектр FFT, то нет эффекта растекания поскольку я взял 500 отсчетов. Если вы берете 512 отсчетов то эффект растекания будет и фаза изменится. При этом если посмотрите из примера разница между отсчетом FFT и герцелем ничтожно мала (1E-12). Так и должно быть поскольку алгоритм Герцеля и рассчитывает один отсчет FFT.

Код

clear all;
close all;
clc;

%% input data
N = 500;
t = 0:N-1;
Fs = 500;
k = 50;

s = sin(2*pi*50*t/Fs);

%% Goertzel algorithm
W = exp(2i*pi*k/N);
v = zeros(1,N);
a = 2*cos(2*pi*k/N);

v(1) = s(1);
v(2) = s(2) + a*v(1);

for n = 3:N
    v(n) = s(n) + a*v(n-1) - v(n-2);
end
G = W*v(n) - v(n-1);

%% fft
S = fft(s);
figure; plot(abs(S))

%%
fprintf('разница между fft и Герцелем: %e\n', G - S(k+1));
fprintf('фаза: %.4f\n', angle(G)*180/pi);

PS по поводу накопления ошибок поставьте в исходных данных

Код

%% input data
N = 500000;
t = 0:N-1;
Fs = 500;
k = 50000;

и все равно ошибка будет на уровне 1E-6 при абсолютном значении спектра 1E6 (т.е. на 12 порядков ниже)

Сообщение отредактировал bahurin - Nov 5 2012, 13:47