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

[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 ?

[ivan219]

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

[анатолий]

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

[kumle]

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

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

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

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

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

[Самурай]

Дело не в накопленной погрешности, дело скорее всего в растекании спектра. Ну не попадают 50Гц входного сигнала ровно на бины ДФТ/БПФ/Герцеля для частоты дискретизации 500 Гц и N=512.

[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

[kumle]

Получаем фазу -90 градусов, поскольку фаза 0 будет если взять s = сos(2*pi*50*t/Fs).

Попробовал получилось так, но теперь непонятно почему если cos то фаза 0, а если sin то 90 Ведь все должно быть наоборот ???

[bahurin]

Попробовал получилось так, но теперь непонятно почему если cos то фаза 0, а если sin то 90 Ведь все должно быть наоборот ???

Комплексная экспонента есть exp(-jwt) = cos(wt)-j*sin(wt). Фаза есть atan(imag(S)/real(S)) . Соответсвенно для косинуса получим фазу 0 или pi поскольку спектр будет чисто вещественным, а для синуса спектр будет чисто комплексным и фаза будет +-pi/2

[kumle]

Да все правильно, это я туплю

[kumle]

Еще вопросик, я сделал по алгоритму Герцеля и сделал ДПФ для одного и того же к.
Тестовый сигнал взял тот же самый в обоих случаях.
Результат получился такой:

Для Герцеля: Re=0; Im=-699 Фаза=-90 градусов

Для ДПФ: Re=0; Im=699 Фаза=90 градусов

Непонятно почему получились разные знаки у мнимой части и соотв. у фазы, ведь я подавал на вход один и тот же сигнал ?

Еще вопросик, я сделал по алгоритму Герцеля и сделал ДПФ для одного и того же к.
Тестовый сигнал взял тот же самый в обоих случаях.
Результат получился такой:

Для Герцеля: Re=0; Im=-699 Фаза=-90 градусов

Для ДПФ: Re=0; Im=699 Фаза=90 градусов

Непонятно почему получились разные знаки у мнимой части и соотв. у фазы, ведь я подавал на вход один и тот же сигнал ?

Сам разобрался, забыл домножить мнимую часть на -1 !!!