Использую алгоритм Герцеля для выделения одной составляющей из сигнала, амплитуда вычисляется верно.
Понадобилось получить также фазу. Вычисляю ее так 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 ?
Полная версия этой страницы: Как получить правильную фазу из алгоритма Герцеля
Форум разработчиков электроники ELECTRONIX.ru > Цифровая обработка сигналов - ЦОС (DSP) > Алгоритмы ЦОС (DSP)
Герцель - Герцелю рознь.
Вот этот правильно считает http://www.dsplib.ru/content/goertzel/goertzel.html
Вот этот правильно считает http://www.dsplib.ru/content/goertzel/goertzel.html
У алгоритма Герцеля, даже с плавающей запятой, при больших N (N>200) накапливается большая погрешность.
Так что -18 может быть и нормально.
Так что -18 может быть и нормально.
Цитата
Герцель - Герцелю рознь.
Вот этот правильно считает http://www.dsplib.ru/content/goertzel/goertzel.html
Вот этот правильно считает http://www.dsplib.ru/content/goertzel/goertzel.html
Самое забавное что я использовал алгоритм именно с этого источника !
Цитата
У алгоритма Герцеля, даже с плавающей запятой, при больших N (N>200) накапливается большая погрешность.
Так что -18 может быть и нормально.
Так что -18 может быть и нормально.
У меня 500 выборок.
Если бы накапливалась погрешность тогда бы и после добавления начальной фазы в сигнал была бы неправильная дельта по фазе, но ведь у меня в сумме получается -18+70=88 градусов, вопрос в том откуда берется изначальная фаза -18 градусов, причем она зависит еще от частоты которую хотим выделить из сигнала, правда пока не понял какая зависимость
Дело не в накопленной погрешности, дело скорее всего в растекании спектра. Ну не попадают 50Гц входного сигнала ровно на бины ДФТ/БПФ/Герцеля для частоты дискретизации 500 Гц и N=512.
если резулатат вашего алгоритма Герцеля отличается от 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.
PS по поводу накопления ошибок поставьте в исходных данных
и все равно ошибка будет на уровне 1E-6 при абсолютном значении спектра 1E6 (т.е. на 12 порядков ниже)
Код
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);
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;
N = 500000;
t = 0:N-1;
Fs = 500;
k = 50000;
и все равно ошибка будет на уровне 1E-6 при абсолютном значении спектра 1E6 (т.е. на 12 порядков ниже)
Цитата
Получаем фазу -90 градусов, поскольку фаза 0 будет если взять s = сos(2*pi*50*t/Fs).
Попробовал получилось так, но теперь непонятно почему если cos то фаза 0, а если sin то 90 Ведь все должно быть наоборот ???
Цитата(kumle @ Nov 14 2012, 13:20) 
Попробовал получилось так, но теперь непонятно почему если cos то фаза 0, а если sin то 90 Ведь все должно быть наоборот ???
Комплексная экспонента есть exp(-jwt) = cos(wt)-j*sin(wt). Фаза есть atan(imag(S)/real(S)) . Соответсвенно для косинуса получим фазу 0 или pi поскольку спектр будет чисто вещественным, а для синуса спектр будет чисто комплексным и фаза будет +-pi/2
Да все правильно, это я туплю
Еще вопросик, я сделал по алгоритму Герцеля и сделал ДПФ для одного и того же к.
Тестовый сигнал взял тот же самый в обоих случаях.
Результат получился такой:
Для Герцеля: Re=0; Im=-699 Фаза=-90 градусов
Для ДПФ: Re=0; Im=699 Фаза=90 градусов
Непонятно почему получились разные знаки у мнимой части и соотв. у фазы, ведь я подавал на вход один и тот же сигнал ?
Сам разобрался, забыл домножить мнимую часть на -1 !!!
Тестовый сигнал взял тот же самый в обоих случаях.
Результат получился такой:
Для Герцеля: Re=0; Im=-699 Фаза=-90 градусов
Для ДПФ: Re=0; Im=699 Фаза=90 градусов
Непонятно почему получились разные знаки у мнимой части и соотв. у фазы, ведь я подавал на вход один и тот же сигнал ?
Цитата
Еще вопросик, я сделал по алгоритму Герцеля и сделал ДПФ для одного и того же к.
Тестовый сигнал взял тот же самый в обоих случаях.
Результат получился такой:
Для Герцеля: Re=0; Im=-699 Фаза=-90 градусов
Для ДПФ: Re=0; Im=699 Фаза=90 градусов
Непонятно почему получились разные знаки у мнимой части и соотв. у фазы, ведь я подавал на вход один и тот же сигнал ?
Тестовый сигнал взял тот же самый в обоих случаях.
Результат получился такой:
Для Герцеля: Re=0; Im=-699 Фаза=-90 градусов
Для ДПФ: Re=0; Im=699 Фаза=90 градусов
Непонятно почему получились разные знаки у мнимой части и соотв. у фазы, ведь я подавал на вход один и тот же сигнал ?
Сам разобрался, забыл домножить мнимую часть на -1 !!!
Для просмотра полной версии этой страницы, пожалуйста, пройдите по ссылке.