Алгоритм Гёрцеля

[Бахурин Сергей]

Пожалуйста, оставляйте ваши вопросы, пожелания и комментарии к статье
"Алгоритм Гёрцеля"


Для быстрого комментирования нажмите кнопку [Ответить] и введите свое имя и сообщение. Также вы можете пройти регистрацию на форуме, после которой вам будет доступна форма быстрого ответа.

[Гость]

Приветствую

По поводу ошибки в структурной схеме, распространенной в литературе. Ошибки на самом деле нет, есть нюанс.
Небольшая модификация скрипта из статьи:

Код: Выделить всё

--- goertzel_test.m	Thu Oct 21 19:56:44 2021
+++ goertzel_test_new.m	Thu Oct 21 19:56:48 2021
@@ -56 +56 @@
-X = filter(b_opp, a_opp, s);
+X = filter(b_opp, a_opp, [s 0]);

...и все заколосилось:

Код: Выделить всё

DFT ....................................  25.0000  -43.3013 i
dspib.org Algorithm.....................  25.0000  -43.3013 i
Oppenheim Algorithm.....................  25.0000  -43.3013 i

Если посмотреть пример в упомянутой справке Matlab, который расположен прямо под структурной схемой, то можно обнаружить ту же фишку.
В английской википедии есть разъяснение ситуации.

[Grot]

Еще один нюанс.
При выводе алгоритма Гёрцеля опираются на свойство (2) в статье или (9.4) в Оппенгейме, поэтому в литературе декларируется использование алгоритма для целочисленных бинов.
По факту, его можно использовать и для нецелых бинов, но при этом вносится погрешность в значение фазы, что для некоторых приложений не существенно.
Более того, если обратить внимание на вышеупомянутую структурную схему справки Matlab, можно обнаружить «лишний» блок на выходе. Этот блок как раз и занимается коррекцией погрешности фазы при нецелых бинах — успешно занимается. Можно сравнить с результатом отсчета DFT. А в конце справки приводится ссылка на статью, в которой этот результат был получен.

[Grot]

Кстати, Лайонс тоже этими вопросами интересовался:
A Simpler Goertzel Algorithm
Goertzel Algorithm for a Non-integer Frequency Index

[Бахурин Сергей]

Во-первых огромнейшее спасибо за очень ценные замечания. Алгортим для нецелых k довольно интересный. Надо будет дописать замечания по этому поводу. И очень редко встречаются такие исчерпывающие комментарии!

Но тем не менее, есть у меня одно такое большое НО!
Смотрите действительно если добавить в вектор входного сигнала 0, то, как вы верно сказали все заколосится:

X = filter(b_opp, a_opp, [s 0]);

А теперь представим: идет поток данных на вход рекурсвного фильтра (скажем фильтр реализован в FPGA и вся схема работает синхронно как конвейер, т.е. каждый такт на входе появляется новый отсчет, и на выходе фильтра отсчет данных).

Каждые N отсчетов на выходе получается одно значение спектра с точностью до фазы.

Теперь для того чтобы строку X = filter(b_opp, a_opp, [s 0]) реализовать, нужно каким то образом пропустить N+1 первый отсчет и впихнуть вместо него 0. Т.е. мы должны потерять N+1 отсчет данных из потока, потому что для оценки согласно приведенной строки, размерность вектора s стала уже не N а N+1. Все это видится как-то кривенько. Почему я не могу подать на вход фильтра поток и каждые N тактов получать выход DFT?

При этом в структуре приведенны в моей статье ничего пропускать (или впихивать дополнительный ноль) не придется, каждые N тактов на выходе будет корректный отсчет спетра.

При этом надо заметить, что умножение выхода фильтра из моей статьи на сайте на коррекцию exp(-j*2*pi*k) также возвращает валидное значение для дробного k, проверил через интерполяционную формулу Дирихле, а также с функцией
goertzel_non_integer_k(x,k) из статьи Лайноса. Тогда при целых k этот корректор вырождается в 1 и перестает влиять.

[Grot]

Все это видится как-то кривенько.

Согласен, непрактичный вариант. Но это не ошибка.
Интересно бы узнать, чего добивались отцы-теоретики таким вариантом, но сильно копать нет мотива.

[Grot]

Есть еще одна тема для алгоритма Гёрцеля.
Вот он эквивалентен сумме отсчетов произведения сигнала на экспоненту, и начальная фаза у экспоненты — нулевая. Если бы удалось обобщить алгоритм Гёрцеля на произвольную начальную фазу, то это открыло бы возможность когерентного сложения субинтервалов накопления.

[Бахурин Сергей]

Честно говоря, мотивы мне тоже не ясны.

Есть предположение, что однажды получили рекуррентное соотношение в том виде, как они везде приводят (вывод кстати приведён в русской Википедии, а в английской нет, хотя раньше он вроде бы был, но с ошибкой). После этого это раскопировалось по вмем учебникам, по пути забыв про выходной фазовый множитель (его нету у оппенгейма и других) авторов.

Потом при реализации нашелся костыль в виде добавления 0 к массиву и всё так и живёт.

Но это лишь предположение, не могу подтвердить, а наоборот надеюсь кто-то обоснует, что есть смысл в этой схеме.

[Бахурин Сергей]

Есть еще одна тема для алгоритма Гёрцеля.
Вот он эквивалентен сумме отсчетов произведения сигнала на экспоненту, и начальная фаза у экспоненты — нулевая. Если бы удалось обобщить алгоритм Гёрцеля на произвольную начальную фазу, то это открыло бы возможность когерентного сложения субинтервалов накопления.

Не понял произвольная начальная фаза чего? Комплексных экспонент в ДПФ? Или самого сигнала?

[Grot]

Не понял произвольная начальная фаза чего? Комплексных экспонент в ДПФ?

Да.
Если мы разбиваем интервал на два субинтервала, то во втором — фаза экспоненты уже не нулевая в общем случае.