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

Алгоритм Гёрцеля (англ. Goertzel algorithm) — это специальная реализация дискретного преобразования Фурье (ДПФ) в форме рекурсивного фильтра. Данный алгоритм был предложен Джеральдом Гёрцелем в 1958 году. В отличие от быстрого преобразования Фурье, вычисляющего все частотные компоненты ДПФ, алгоритм Гёрцеля позволяет эффективно вычислить значение одного частотного компонента.

Алгоритм Гёрцеля является популярным алгоритмом для решения задачи детектирования и декодирования тональных сигналов в телефонии.

Варианты названия алгоритма

В русскоязычной литературе нет устоявшегося варианта транскрипции фамилии автора алгоритма. Распространены варианты «Алгоритм Герцеля», «Алгоритм Гертцеля», «Алгоритм Горцеля» и другие.

Описание алгоритма

Пусть [math]\displaystyle{ x_n,\ n = 0, \dots, N - 1 }[/math] — измеренные значения сигнала, которые являются входными данными для дискретного преобразования Фурье, а [math]\displaystyle{ X_k,\ k = 0, \dots, N - 1 }[/math] — частотные компоненты дискретного преобразования Фурье, по определению равные [math]\displaystyle{ X_k = \sum_{n=0}^{N-1} x_n e^{-\frac{2 \pi i}N k n} }[/math]. Для расчёта [math]\displaystyle{ X_k }[/math] с помощью алгоритма Гёрцеля:

1. Последовательно вычисляются члены последовательности [math]\displaystyle{ s_n }[/math] для [math]\displaystyle{ n = 0, \dots, N - 1 }[/math] по рекуррентной формуле [math]\displaystyle{ s_n = 2 \cos\left(\frac{2 \pi k}N\right) s_{n-1} - s_{n-2} + x_n }[/math], где [math]\displaystyle{ s_{-1} = s_{-2} = 0 }[/math].
2. Искомое значение частотного компонента получается как [math]\displaystyle{ X_k = e^{\frac{2 \pi i}N k} s_{N-1} - s_{N-2} }[/math].

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

[math]\displaystyle{ |X_k|^2 = s_{N-1}^2 - 2 \cos\left(\frac{2 \pi k}N\right) s_{N-1} s_{N-2} + s_{N-2}^2. }[/math]

Вывод

Для удобства записи введём обозначение [math]\displaystyle{ W = e^{-\frac{2 \pi i}N} }[/math]. Определение k-го частотного компонента дискретного преобразования Фурье тогда примет вид [math]\displaystyle{ X_k = \sum_{n=0}^{N-1} x_n W^{kn} }[/math]. Также заметим, что [math]\displaystyle{ W^{-N} = e^{2 \pi i} = 1 }[/math], и для любого действительного [math]\displaystyle{ k }[/math] выполняется равенство

[math]\displaystyle{ W^k + W^{-k} = e^{-\frac{2 \pi i}N k} + e^{\frac{2 \pi i}N k} = 2 \cos\left(\frac{2 \pi k}N\right). }[/math]

Преобразуем исходное определение [math]\displaystyle{ X_k }[/math]:

[math]\displaystyle{ \begin{align} X_k & = \sum_{n=0}^{N-1} {x_n W^{kn}} = \\ & = W^{-kN} \sum_{n=0}^{N-1} {x_n W^{kn}} = \\ & = \sum_{n=0}^{N-1} {x_n [W^{-k}]^{N-n}} = \\ & = x_0 [W^{-k}]^N + x_1 [W^{-k}]^{N-1} + \dots + x_{N-2} [W^{-k}]^2 + x_{N-1} W^{-k} = \\ & = \bigg( \Big( \big( \dots (x_0 W^{-k} + x_1) W^{-k} + \dots \big) W^{-k} + x_{N-2} \Big) W^{-k} + x_{N-1} \bigg) W^{-k}. \end{align} }[/math]

Процесс последовательного вычисления стоящих во вложенных скобках выражений мы теперь можем записать с помощью рекуррентной формулы:

[math]\displaystyle{ y_n = W^{-k} y_{n-1} + x_n,\ n = 0, 1, 2, \dots, N - 1. }[/math]

При этом для вычисления [math]\displaystyle{ y_0 }[/math] примем [math]\displaystyle{ y_{-1} = 0 }[/math], а искомое значение k-го частотного компонента выражается как [math]\displaystyle{ X_k = W^{-k} y_{N-1} }[/math]. Полученное выражение для [math]\displaystyle{ y_n }[/math] представляет собой БИХ-фильтр первого порядка:

[math]\displaystyle{ y_n - W^{-k} y_{n-1} = x_n. }[/math]

(1)

Применив z-преобразование к обеим частям выражения (1), получим

[math]\displaystyle{ \big(1 - W^{-k} z^{-1}\big) Y(z) = X(z), }[/math]

где [math]\displaystyle{ X(z) }[/math] и [math]\displaystyle{ Y(z) }[/math] — образы дискретных временных сигналов [math]\displaystyle{ x_n }[/math] и [math]\displaystyle{ y_n }[/math] соответственно. Из этого выражения выразим передаточную функцию фильтра (1) и преобразуем её, домножив числитель и знаменатель на [math]\displaystyle{ 1 - W^k z^{-1} }[/math]:

[math]\displaystyle{ \begin{align} H(z) & = \frac{Y(z)}{X(z)} =\\ & = \frac 1 {1 - W^{-k} z^{-1}} = \\ & = \frac{1 - W^k z^{-1}}{(1 - W^k z^{-1})(1 - W^{-k} z^{-1})} = \\ & = \frac{1 - W^k z^{-1}}{1 - (W^k +W^{-k})z^{-1} + z^{-2}} = \\ & = \frac{1 - W^k z^{-1}}{1 - 2 \cos \left( \frac{2 \pi k}N \right) z^{-1} + z^{-2}}. \end{align} }[/math]

Передаточную функцию [math]\displaystyle{ H(z) }[/math] можно представить в виде произведения функций [math]\displaystyle{ P(z) = \frac 1 {1 - 2 \cos\left(\frac{2 \pi k}N\right) z^{-1} + z^{-2}} }[/math] и [math]\displaystyle{ Q(z) = 1 - W^k z^{-1} }[/math], что во временно́й области соответствует последовательному соединению двух фильтров с передаточными функциями [math]\displaystyle{ P(z) }[/math] и [math]\displaystyle{ Q(z) }[/math]. Таким образом, исходный фильтр (1) может быть представлен в виде комбинации двух фильтров:

[math]\displaystyle{ s_n - 2 \cos \left(\frac{2 \pi k}N\right) s_{n-1} + s_{n-2} = x_n, }[/math]

(2)

[math]\displaystyle{ y_n = s_n - W^k s_{n-1}. }[/math]

(3)

Поскольку для вычисления [math]\displaystyle{ X_k }[/math] не требуется всех выходных значений фильтра (3), а только значение [math]\displaystyle{ y_{N-1} }[/math], то вычисление [math]\displaystyle{ y_0, y_1, \dots, y_{N-2} }[/math] можно опустить благодаря тому, что фильтр (3) является КИХ-фильтром, выход которого [math]\displaystyle{ y_n }[/math] зависит только от входов [math]\displaystyle{ s_n }[/math] и [math]\displaystyle{ s_{n-1} }[/math]. Значение [math]\displaystyle{ X_k }[/math] можем тогда выразить следующим образом:

[math]\displaystyle{ X_k = W^{-k} y_{N-1} = W^{-k} \big(s_{N-1} - W^k s_{N-2}\big) = W^{-k} s_{N-1} - s_{N-2}. }[/math]

Мощность сигнала тогда может быть вычислена как

[math]\displaystyle{ \begin{align} |X_k|^2 & = X_k \cdot \overline X_k = \\ & = (W^{-k} s_{N-1} - s_{N-2})(\overline{[W^{-k}]} s_{N-1} - s_{N-2}) = \\ & = (W^{-k} s_{N-1} - s_{N-2})(W^k s_{N-1} - s_{N-2}) = \\ & = s_{N-1}^2 - (W^k + W^{-k}) s_{N-1} s_{N-2} + s_{N-2}^2 = \\ & = s_{N-1}^2 - 2 \cos \left( \frac{2 \pi k}N \right) s_{N-1} s_{N-2} + s_{N-2}^2. \end{align} }[/math]

Устойчивость алгоритма

Процесс рекуррентного вычисления членов последовательности [math]\displaystyle{ s_n }[/math] является цифровым БИХ-фильтром второго порядка. Как и любой БИХ-фильтр, он чувствителен к ошибкам, которые возникают в результате квантования и использования арифметических операций со словами конечной длины. Более того, поскольку оба полюса фильтра ([math]\displaystyle{ z=e^{ -\frac{2 \pi i}N} }[/math] и [math]\displaystyle{ z=e^{ \frac{2 \pi i}N} }[/math]) лежат на единичной окружности, ошибки округления могут привести и к неустойчивости фильтра. В связи с этим, алгоритм Гёрцеля следует применять с осторожностью при больших длинах окон (большие значения [math]\displaystyle{ N }[/math]), особенно при использовании арифметики с низкой разрядностью.

Вычислительная эффективность алгоритма

Для вычисления одного частотного компонента ДПФ комплексной последовательности отсчётов длины [math]\displaystyle{ N }[/math] с помощью алгоритма Гёрцеля требуется [math]\displaystyle{ 2N+4 }[/math] умножений и [math]\displaystyle{ 4N+4 }[/math] сложений / вычитаний (для действительной последовательности — [math]\displaystyle{ N+2 }[/math] умножения и [math]\displaystyle{ 2N+1 }[/math] сложений), не считая затрат на вычисление постоянных коэффициентов [math]\displaystyle{ 2 \cos \tfrac{2 \pi k}N }[/math] и [math]\displaystyle{ e^{ \frac{2 \pi i}N k} }[/math]. При этом метод не требует хранения каких-либо таблиц коэффициентов, а основной объём арифметических вычислений метода может производиться по мере поступления входных отсчётов [math]\displaystyle{ x_n }[/math].

Вычисление одного частотного компонента непосредственно по формуле-определению дискретного преобразования Фурье требует большего числа арифметических действий (в комплексном случае [math]\displaystyle{ 4N }[/math] умножений и [math]\displaystyle{ 4N }[/math] сложений, а в действительном [math]\displaystyle{ 2N }[/math] умножений и [math]\displaystyle{ 4N }[/math] сложений), чем алгоритм Гёрцеля, хотя асимптотически требует также [math]\displaystyle{ O(N) }[/math] действий. Кроме того, для эффективного проведения вычислений по прямой формуле требуется хранить таблицу коэффициентов.

Другой альтернативой алгоритму Гёрцеля является быстрое преобразование Фурье. БПФ-алгоритм Кули — Тьюки по основанию 2 требует [math]\displaystyle{ 2 \log_2 N }[/math] умножений и [math]\displaystyle{ 3 \log_2 N }[/math] сложений на вычисление одного частотного компонента (а в случае действительной входной последовательности количество арифметических действий уменьшается вдвое относительно приведённых значений), однако все частотные компоненты должны вычисляться одновременно. Когда число частотных компонентов, подлежащих вычислению, невелико, то применение алгоритма Гёрцеля эффективнее, чем применение БПФ. Если необходимо вычислить M частотных компонентов, то применение алгоритма Гёрцеля выгоднее при условии

[math]\displaystyle{ M \lt \frac{5}{6} \log_2 N. }[/math]

Кроме того, алгоритмы БПФ должны применяться к полным блокам данных длины N и могут требовать хранения таблиц коэффициентов для эффективной реализации. Семейство алгоритмов быстрого преобразования Фурье включает алгоритмы с различными свойствами и вычислительной эффективностью, в том числе и так называемые усечённые алгоритмы БПФ, позволяющие вычислять подмножество набора частотных компонентов. В каждом случае решение вопроса о предпочтительности использования алгоритма Гёрцеля или БПФ зависит от выбора конкретного варианта БПФ, длины блока входных данных и количества частотных компонентов, которые необходимо вычислить.

Литература

• Эммануил С. Айфичер, Барри У. Джервис. Цифровая обработка сигналов: практический подход = Digital Signal Processing: A Practical Approach. — 2-е издание. — М.: Издательский дом «Вильямс», 2004. — 992 с. — ISBN 5-8459-0710-1.
• Блейхут Р. Быстрые алгоритмы цифровой обработки сигналов = Fast Algorithms for Digital Signal Processing. — М.: Мир, 1989. — 448 с. — ISBN 5-09-001009-2.

Ссылки

• Алгоритм Герцеля (Goertzel algorithm)
• Динамический пересчет спектральных отсчетов на каждом такте дискретизации. Модифицированный алгоритм Герцеля