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

Алгоритм Герцеля — это метод цифровой обработки сигналов (DSP) для эффективной оценки отдельных членов дискретного преобразования Фурье (DFT). Это полезно в некоторых практических приложениях, таких как распознавание тонов двухтональной многочастотной сигнализации (DTMF), производимых кнопками клавиатуры традиционного аналогового телефона . Алгоритм был впервые описан Джеральдом Герцелем в 1958 году.

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

Алгоритм Герцеля также можно использовать «обратно» в качестве функции синтеза синусоиды, которая требует только 1 умножения и 1 вычитания на сгенерированную выборку.

Алгоритм

Основной расчет в алгоритме Герцеля имеет форму цифрового фильтра , и по этой причине алгоритм часто называют фильтром Герцеля . Фильтр работает с входной последовательностью $x[n]$ в каскаде из двух каскадов с параметром $\omega _{0}$, что дает анализируемую частоту, нормализованную в радианах на выборку.

На первом этапе рассчитывается промежуточная последовательность, $s[n]$:

$s[n]=x[n]+2\cos(\omega _{0})s[n-1]-s[n-2].$ (1)

На втором этапе применяется следующий фильтр к $s[n]$, создавая выходную последовательность $y[n]$:

$y[n]=s[n]-e^{-j\omega _{0}}s[n-1].$ (2)

второго порядка Можно наблюдать, что первая ступень фильтра представляет собой БИХ-фильтр со структурой прямой формы . Эта конкретная структура обладает тем свойством, что ее внутренние переменные состояния равны прошлым выходным значениям этого этапа. Входные значения $x[n]$ для $n<0$ предполагаются равными 0. Чтобы установить начальное состояние фильтра, чтобы можно было начать оценку с выборки $x[0]$состояниям фильтра присваиваются начальные значения $s[-2]=s[-1]=0$. Чтобы избежать опасности наложения псевдонимов , частота $\omega _{0}$ часто ограничивается диапазоном от 0 до π (см. теорему выборки Найквиста – Шеннона ); использование значения вне этого диапазона не бессмысленно, но эквивалентно использованию совмещенной частоты внутри этого диапазона, поскольку экспоненциальная функция является периодической с периодом 2π в $\omega _{0}$.

Фильтр второй ступени можно рассматривать как КИХ-фильтр , поскольку в его расчетах не используются какие-либо из его прошлых выходных данных.

Методы Z-преобразования могут применяться для изучения свойств каскада фильтров. Z-преобразование первой ступени фильтра, заданное в уравнении (1), равно

${\begin{aligned}{\frac {S(z)}{X(z)}}&={\frac {1}{1-2\cos(\omega _{0})z^{-1}+z^{-2}}}\\&={\frac {1}{(1-e^{+j\omega _{0}}z^{-1})(1-e^{-j\omega _{0}}z^{-1})}}.\end{aligned}}$ (3)

Z-преобразование второй ступени фильтра, заданное в уравнении (2), равно

${\frac {Y(z)}{S(z)}}=1-e^{-j\omega _{0}}z^{-1}.$ (4)

Тогда объединенная передаточная функция каскада двух ступеней фильтра равна

${\begin{aligned}{\frac {S(z)}{X(z)}}{\frac {Y(z)}{S(z)}}={\frac {Y(z)}{X(z)}}&={\frac {(1-e^{-j\omega _{0}}z^{-1})}{(1-e^{+j\omega _{0}}z^{-1})(1-e^{-j\omega _{0}}z^{-1})}}\\&={\frac {1}{1-e^{+j\omega _{0}}z^{-1}}}.\end{aligned}}$ (5)

Это можно преобразовать обратно в эквивалентную последовательность во временной области, а термины развернуть обратно до первого входного термина по индексу. $n=0$:

${\begin{aligned}y[n]&=x[n]+e^{+j\omega _{0}}y[n-1]\\&=\sum _{k=-\infty }^{n}x[k]e^{+j\omega _{0}(n-k)}\\&=e^{j\omega _{0}n}\sum _{k=0}^{n}x[k]e^{-j\omega _{0}k}\qquad {\text{since }}\forall k<0,x[k]=0.\end{aligned}}$ (6)

Численная стабильность

Можно заметить, что полюса фильтра Z-преобразования расположены в точках $e^{+j\omega _{0}}$ и $e^{-j\omega _{0}}$, на окружности единичного радиуса с центром в начале плоскости комплексного Z-преобразования. Это свойство указывает на то, что процесс фильтрации минимально стабилен и уязвим к накоплению числовых ошибок при вычислении с использованием арифметических вычислений низкой точности и длинных входных последовательностей. Численно устойчивая версия была предложена Кристианом Рейншем .

Вычисления ДПФ

В важном случае вычисления члена ДПФ применяются следующие специальные ограничения.

• Фильтрация завершается по индексу $n=N$, где $N$ — количество членов во входной последовательности ДПФ.
• Частоты, выбранные для анализа Герцеля, ограничены специальной формой

$\omega _{0}=2\pi {\frac {k}{N}}.$ (7)

• Индексный номер $k$ указывающий, что «частотный элемент» ДПФ выбран из набора индексных номеров

$k\in \{0,1,2,...,N-1\}.$ (8)

Сделав эти замены в уравнении (6) и заметив, что член $e^{+j2\pi k}=1$, уравнение (6) тогда примет следующий вид:

$y[N]=\sum _{n=0}^{N}x[n]e^{-j2\pi {\frac {nk}{N}}}.$ (9)

Мы можем заметить, что правая часть уравнения (9) чрезвычайно похожа на определяющую формулу для члена ДПФ. $X[k]$, член ДПФ для индексного номера $k$, но не совсем то же самое. Суммирование, показанное в уравнении (9), требует $N+1$ входные параметры, но только $N$ входные условия доступны при оценке ДПФ. Простой, но неэлегантный способ — расширить входную последовательность $x[n]$ с еще одним искусственным значением $x[N]=0$. Из уравнения (9) мы видим, что математическое влияние на конечный результат такое же, как и удаление члена $x[N]$ от суммирования, таким образом получая предполагаемое значение ДПФ.

Однако существует более элегантный подход, позволяющий избежать дополнительного прохода фильтра. Из уравнения (1) мы можем отметить, что когда расширенный входной член $x[N]=0$ используется на последнем этапе,

$s[N]=2\cos(\omega _{0})s[N-1]-s[N-2].$ (10)

Таким образом, алгоритм можно завершить следующим образом:

• завершить БИХ-фильтр после обработки входного термина $x[N-1]$,
• примените уравнение (10) для построения $s[N]$ из предыдущих результатов $s[N-1]$ и $s[N-2]$,
• применить уравнение (2) с рассчитанным $s[N]$ ценность и с $s[N-1]$ производится окончательным прямым расчетом фильтра.

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

${\begin{aligned}y[N]&=s[N]-e^{-j2\pi {\frac {k}{N}}}s[N-1]\\&=(2\cos(\omega _{0})s[N-1]-s[N-2])-e^{-j2\pi {\frac {k}{N}}}s[N-1]\\&=e^{j2\pi {\frac {k}{N}}}s[N-1]-s[N-2].\end{aligned}}$ (11)

Обратите внимание, что остановка обновлений фильтра в установленный срок $N-1$ и немедленное применение уравнения (2), а не уравнения (11), пропускает окончательные обновления состояния фильтра, что дает результат с неправильной фазой.

Конкретная структура фильтрации, выбранная для алгоритма Герцеля, является ключом к его эффективным расчетам ДПФ. Мы можем заметить, что только одно выходное значение $y[N]$ используется для расчета ДПФ, поэтому расчеты для всех остальных выходных членов опускаются. Поскольку КИХ-фильтр не рассчитывается, расчеты этапа БИХ $s[0],s[1]$и т. д. можно отбросить сразу после обновления внутреннего состояния первого этапа.

Кажется, это создает парадокс: для завершения алгоритма этап КИХ-фильтра должен быть оценен один раз с использованием двух последних выходных данных этапа БИХ-фильтра, в то время как для эффективности вычислений итерация БИХ-фильтра отбрасывает его выходные значения. Здесь применяются свойства структуры фильтра прямой формы. Две внутренние переменные состояния БИХ-фильтра предоставляют два последних значения выходного сигнала БИХ-фильтра, которые являются условиями, необходимыми для оценки ступени КИХ-фильтра.

Приложения

Условия спектра мощности

Анализ уравнения (6): окончательный проход БИХ-фильтра для расчета члена $y[N]$ используя дополнительное входное значение $x[N]=0$ применяет комплексный множитель величины 1 к предыдущему члену $y[N-1]$. Следовательно, $y[N]$ и $y[N-1]$ представляют собой эквивалентную мощность сигнала. В равной степени справедливо применить уравнение (11) и вычислить мощность сигнала из члена $y[N]$ или применить уравнение (2) и вычислить мощность сигнала из члена $y[N-1]$. Оба случая приводят к следующему выражению для мощности сигнала, представленному термином ДПФ $X[k]$:

${\begin{aligned}X[k]\,X'[k]&=y[N]\,y'[N]=y[N-1]\,y'[N-1]\\&=s^{2}[N-1]+s^{2}[N-2]-2\cos \left(2\pi {\frac {k}{N}}\right)\,s[N-1]\,s[N-2].\end{aligned}}$ (12)

В приведенном ниже псевдокоде входные данные с действительным знаком хранятся в массиве x и переменные sprev и sprev2 временно сохранить историю вывода БИХ-фильтра. Nterms количество выборок в массиве, а Kterm соответствует интересующей частоте, умноженной на период выборки.

Nterms defined here
Kterm selected here
ω = 2 × π × Kterm / Nterms;
coeff := 2 × cos(ω)

sprev := 0
sprev2 := 0
for each index n in range 0 to Nterms-1 do
    s := x[n] + coeff × sprev - sprev2
    sprev2 := sprev
    sprev := s
end

power := sprev2 + sprev22 - (coeff × sprev × sprev2)

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

Один термин ДПФ с вещественной арифметикой

Случай вещественных входных данных возникает часто, особенно во встроенных системах, где входные потоки являются результатом прямых измерений физических процессов. Когда входные данные имеют вещественное значение, внутренние переменные состояния фильтра sprev и sprev2 Также можно наблюдать, что они имеют действительные значения, следовательно, на первом этапе БИХ не требуется никаких сложных арифметических действий. Оптимизация для арифметики с действительными значениями обычно так же проста, как применение соответствующих типов данных с действительными значениями для переменных.

После вычислений с использованием входного термина $x[N-1]$и итерации фильтра завершаются, для оценки члена ДПФ необходимо применить уравнение (11). В окончательном вычислении используется комплексная арифметика, но ее можно преобразовать в вещественную арифметику, разделив действительные и мнимые члены:

${\begin{aligned}c_{r}&=\cos(2\pi {\tfrac {k}{N}}),\\c_{i}&=\sin(2\pi {\tfrac {k}{N}}),\\y[N]&=c_{r}s[N-1]-s[N-2]+jc_{i}s[N-1].\end{aligned}}$ (13)

По сравнению с приложением анализа спектра мощности, единственная разница заключается в расчете, используемом для завершения:

(Same IIR filter calculations as in the signal power implementation)
XKreal = sprev * cr - sprev2;
XKimag = sprev * ci;

Обнаружение фазы

Это приложение требует такой же оценки термина ДПФ. $X[k]$, как обсуждалось в предыдущем разделе, с использованием входного потока с действительным или комплексным знаком. Тогда фазу сигнала можно оценить как

$\phi =\tan ^{-1}{\frac {\Im (X[k])}{\Re (X[k])}},$ (14)

принимая соответствующие меры предосторожности для сингулярностей, квадрантов и т. д. при вычислении функции обратного тангенса.

Сложные сигналы в реальной арифметике

Поскольку сложные сигналы линейно разлагаются на действительную и мнимую части, алгоритм Герцеля можно вычислять в вещественной арифметике отдельно по последовательности действительных частей, что дает $y_{\text{r}}[n]$, и по последовательности мнимых частей, что дает $y_{\text{i}}[n]$. После этого два частичных результата с комплексными значениями можно повторно объединить:

$y[n]=y_{\text{r}}[n]+jy_{\text{i}}[n].$ (15)

Вычислительная сложность

• Согласно теории сложности вычислений , вычисление набора $M$ Условия ДПФ с использованием $M$ применение алгоритма Герцеля к набору данных с $N$ значения с «стоимостью за операцию» $K$ имеет сложность $O(KNM)$.

Чтобы вычислить один ДПФ бин $X(f)$ для сложной входной последовательности длины $N$, алгоритм Герцеля требует $2N$ умножения и $4\ N$ сложения/вычитания в цикле, а также 4 умножения и 4 конечных сложения/вычитания, всего $2N+4$ умножения и $4N+4$ добавления/вычитания. Это повторяется для каждого из $M$ частоты.

• Напротив, использование БПФ для набора данных с $N$ значения имеют сложность $O(KN\log _{2}(N))$.

Это сложнее применить напрямую, поскольку это зависит от используемого алгоритма БПФ, но типичным примером является БПФ по основанию 2, для которого требуется $2\log _{2}(N)$ умножения и $3\log _{2}(N)$ сложения/вычитания на бин ДПФ для каждого из $N$ контейнеры.

В выражениях порядка сложности, когда количество вычисляемых термов $M$ меньше, чем $\log N$, преимущество алгоритма Герцеля очевидно. Но поскольку код БПФ сравнительно сложен, коэффициент «затраты на единицу работы» $K$ часто больше для БПФ, и практическое преимущество в пользу алгоритма Герцеля даже для $M$ в несколько раз больше, чем $\log _{2}(N)$.

В качестве практического правила для определения того, является ли БПФ по основанию 2 или алгоритм Герцеля более эффективным, отрегулируйте количество членов. $N$ в наборе данных до ближайшей точной степени 2, называя это $N_{2}$, и алгоритм Герцеля, вероятно, будет быстрее, если

$M\leq {\frac {5N_{2}}{6N}}\log _{2}(N_{2})$

Реализации БПФ и платформы обработки оказывают значительное влияние на относительную производительность. Некоторые реализации БПФ выполнять внутренние вычисления комплексных чисел для генерации коэффициентов «на лету», значительно увеличивая их «стоимость K на единицу работы». Алгоритмы БПФ и ДПФ могут использовать таблицы предварительно вычисленных значений коэффициентов для повышения числовой эффективности, но это требует большего доступа к значениям коэффициентов, буферизованным во внешней памяти, что может привести к увеличению конкуренции за кэш, что сводит на нет некоторое числовое преимущество.

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

См. также

• Алгоритм Блюстейна БПФ (chirp-Z)
• Частотная манипуляция (FSK)
• Фазовая манипуляция (PSK)

Дальнейшее чтение

• Проакис, Дж.Г.; Манолакис, Д.Г. (1996), Цифровая обработка сигналов: принципы, алгоритмы и приложения , Аппер-Сэддл-Ривер, Нью-Джерси: Прентис-Холл, стр. 480–481, Bibcode : 1996dspp.book.....P

Внешние ссылки

• Алгоритм Герцеля в Wayback Machine (архивировано 28 июня 2018 г.)
• Алгоритм DSP для частотного анализа
• Алгоритм Герцеля Кевина Бэнкса