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

Содержание
Введение
В предыдущих разделах мы рассмотрели дискретное преобразование Фурье (ДПФ) и его свойства, а также алгоритмы быстрого преобразования Фурье (БПФ). Мы выяснили, что алгоритмы БПФ позволяют существенно экономить вычислительные ресурсы при расчете ДПФ. Однако иногда требуется произвести расчет не полного ДПФ, а лишь фиксированного количества спектральных отсчетов. Например, такая задача встает при декодировании DTMF сигналов. В этом случае применение алгоритмов БПФ может оказаться неэффективным, поскольку для получения одного спектрального отсчета потребуется произвести расчет полного ДПФ.
В данном разделе мы рассмотрим алгоритм Гёрцеля [1], используемого для расчета фиксированного набора спектральных отсчетов дискретного преобразования Фурье, который нашел широкое применение в задачах декодирования DTMF сигналов. Как будет показано ниже, алгоритм Гёрцеля может быть реализован в виде БИХ-фильтра второго порядка.
Рекуррентное соотношение для расчета фиксированного спектрального отсчета сигнала
Выражение для $N-$точечного дискретного преобразования Фурье сигнала $s(n)$, $n = 0 \ldots N-1$ имеет вид:
$$ S(k) = \sum_{n = 0}^{N-1} s(n) \cdot W_N^{n \cdot k}; $$ $$ W_N^{n \cdot k} = \exp \left( - j \cdot \frac{2 \pi}{N} \cdot n \cdot k \right); \textrm{ } k = 0 \ldots N-1. $$
(1)
Поворотные коэффициенты $W_N^{n \cdot k}$ обладают следующим свойством:
$$ W_N^{- k \cdot N} = \exp \left( j \cdot \frac{2 \pi}{N} \cdot N \cdot n \right) = 1; \textrm{ } k = 0 \ldots N-1. $$
(2)
Таким образом, умножение выражения (1) на $W_N^{- k \cdot N}$ не приведет к изменению результата:
$$ S(k) = W_N^{-k \cdot N} \sum_{n = 0}^{N-1} s(n) \cdot W_N^{n \cdot k} = \sum_{n = 0}^{N-1} s(n) \cdot W_N^{k \cdot (n-N)}. $$
(3)
Распишем сумму (3):
$$ S(k) =s(0) \cdot W_N^{-k \cdot N} + s(1) \cdot W_N^{-k \cdot (N-1)} + s(2) \cdot W_N^{-k \cdot (N-2)} + \ldots $$ $$ \ldots + s(N-2) \cdot W_N^{-k \cdot 2} + s(N-1) \cdot W_N^{-k}. $$
(4)
Обозначим $S(k) = y_{N-1}(k)$ для фиксированного номера $k$ спектрального отсчета ДПФ. Также вынесем в (4) $W_N^{-k}$ за скобки и получим:
$$ y_{N-1}(k) = W_N^{-k} \cdot \Big( s(N-1) + \underbrace{W_N^{-k} \cdot s(N-2) + \ldots + s(0) \cdot W_N^{-k \cdot (N-1)} }_{=y_{N-2}(k)}\Big) = $$ $$ = W_N^{-k} \cdot \big( s(N-1) + y_{N-2}(k)\big) . $$
(5)
Выразим $y_{N-2}(k)$ через $y_{N-3}(k)$. Для этого снова вынесем $W_N^{-k}$ за скобки и получим:
$$ y_{N-2}(k) = W_N^{-k} \cdot s(N-2) + W_N^{-2\cdot k} \cdot s(N-3) + \ldots + s(0) \cdot W_N^{-k \cdot (N-1)} = $$ $$ = W_N^{-k} \cdot \Big( s(N-2) + \underbrace{W_N^{-k} \cdot \ s(N-3) + y_{N-2}(k) \ldots + s(0) \cdot W_N^{-k \cdot (N-2)}}_{=y_{N-3}(k)} \Big) . $$
(6)
Таким образом, мы получили рекуррентное соотношение для вычисления $y_{r}(k)$ на любом шаге $r = 0 \ldots N-1$ через $y_{r-1}(k)$, полученных на предыдущем шаге для фиксированного номера спектрального отсчета $k$:
Рисунок 1. Структурная схема БИХ-фильтра,
реализующего расчет спектрального
отсчета с номером $k$
$$ y_r(k) = W_N^{-k} \cdot \Big( s(r) + y_{r-1}(k) \Big), $$
(7)
где $s(r)$ – отсчет входного сигнала с номером $r$.
Данное рекуррентное соотношение на шаге с номером $r = N-1$ приводит нас к $S(k) = y_{N-1}(k)$, т.е. на $N-1$ шаге получим спектральный отсчет с номером $k$.
Проанализировав (7), обратим внимание, что рекуррентное соотношение можно трактовать как разностное уравнение БИХ-фильтра первого порядка с комплексным коэффициентом $W_N^{-k}$, структурная схема которого приведена на рисунке 1.
Обозначим через $X(z)$ и $Y(z)$ $z-$образы $x(r)$ и $y_r(k)$ соответственно. Тогда $z^{-1} \cdot Y(z) -$это $z-$образ $y_{r-1}(k),$ и разностное уравнение (7) в операторном виде равно
$$ Y(z) = W_N^{-k} \cdot \left( X(z) + z^{-1} \cdot Y(z) \right). $$
(8)
Передаточная характеристика полученного БИХ-фильтра имеет вид:
$$ H(z) = \frac{Y(z)}{X(z)} = \frac{W_N^{-k}}{1-W_N^{-k} \cdot z^{-1}}. $$
(9)
Таким образом, мы зафиксировали номер спектрального отсчета $k$ и преобразовали выражение ДПФ для одного фиксированного спектрального отсчета к рекуррентному соотношению, которое позволяет нам на шаге $N-1$ получить значение искомого $k$-го спектрального отсчета $S(k)$. При этом все промежуточные значения рекуррентного соотношения нас не интересуют, а интересует только $y_{N-1}(k) = S(k).$
Затем мы трактовали рекуррентное соотношение как разностное уравнение БИХ-фильтра с передаточной характеристикой (9). В результате мы получили фильтр первого порядка с комплексным коэффициентом $W_N^{-k},$ применение которого на $N-1$-ой итерации на выходе дает $S(k)$. При фиксированном $k$ значение комплексного коэффициента $W_N^{-k}$ всегда постоянно.
Алгоритм Гёрцеля
Для вычисления одного спектрального отсчета $S(k)$ при использовании БИХ-фильтра требуется $N$ комплексных умножений и сложений, что не дает никаких преимуществ по сравнению с прямым вычислением ДПФ, согласно выражению (1).
Однако сокращение вычислений можно получить, если домножить числитель и знаменатель передаточной характеристики БИХ-фильтра (9) на $1-W_N^k \cdot z^{-1}$:
$$ H(z) = \frac{W_N^{-k}}{1-W_N^{-k} \cdot z^{-1}} = \frac{W_N^{-k} \cdot \left(1-W_N^k \cdot z^{-1} \right)}{ \left(1-W_N^{-k} \cdot z^{-1} \right) \cdot \left(1-W_N^k \cdot z^{-1} \right)} = \frac{W_N^{-k} - z^{-1} }{1- z^{-1} \cdot \left( W_N^{-k} + W_N^{k} \right) + z^{-2}}. $$
(10)
Рассмотрим более подробно сумму:
$$ W_N^{-k} + W_N^{k} = \exp \left( j \cdot \frac{2\pi}{N} \cdot k \right) + \exp \left( -j \cdot \frac{2\pi}{N} \cdot k \right) = 2 \cdot \cos \left( \frac{2\pi}{N} \cdot k \right). $$
(11)
Тогда (10) с учетом(11) можно представить следующим образом:
$$ H(z) = \frac{W_N^{-k} - z^{-1} }{1- 2 \cdot \cos \left( \frac{2\pi}{N} \cdot k \right) \cdot z^{-1} + z^{-2}}. $$
(12)
Передаточная характеристика (12) соответствует БИХ-фильтру второго порядка, стуктурная схема которого показана на рисунке 2, где $\alpha = 2 \cdot \cos \left( 2\pi \cdot \frac{k}{N}\right)$.
Рисунок 2. Структурная схема БИХ-фильтра
реализующего расчет спектрального
отсчета с номером $k$
Обратим внимание, что $\alpha = 2 \cdot \cos \left( 2\pi \cdot \frac{k}{N}\right)$ – вещественный коэффициент, соответственно умножения в рекурсивной ветви фильтра стали вещественными, причем умножения на $-1$ можно не учитывать. У нас остался один комплексный коэффициент $W_N^{-k}$ в числителе передаточной характеристики. Но мы выяснили, что промежуточные значения $y_r(k)$ нам неважны, соответственно умножать на $W_N^{-k}$ можно только на $N-1$ итерации, когда рассчитывается непосредственно $S(k)$.
Таким образом, применение фильтра второго порядка позволяет рассчитать $S(k)$ для фиксированного $k$ с использованием $N$ умножений, но не комплексных, а вещественных. Также требуется одно комплексное умножение на $W_N^{-k}$ на последней итерации.
Внимание! В литературе [2,3] приведена структурная схема, отличная от рисунка 2. Проверка показала, что приведенный в литературе [2,3] фильтр второго порядка выдает ошибочный результат на выходе. Фильтр, приведенный на рисунке 2, выдает точное значения отсчета ДПФ. Вы можете самостоятельно удостовериться в правильности результата фильтра, представленного на рисунке 2, выполнив проверочный скрипт goertzel_test.m в системе Matlab или в GNU Octave, или goertzel_test.py на языке Python.
								
% N-точечное ДПФ
N = 100;
% Номер спектрального отсчета (бина)
k = 32;

%% Исходный сигнал и расчет его ДПФ    
s = sin(2*pi*(0:N-1)*k/N+pi/6);
X = fft(s);
S_dft = X(k+1); % k+1 потому что Matlab и GNU Octave индексируют массивы начиная с единицы 


%%  Алгоритм Герцеля с DSPLIB.ORG. Формула (12) и рисунок 2.
%
%                W - z^-1
% H(z) = ---------------------------
%          1 - alpha * z^-1 + z^-2
%
%
%	W = W_N^(-k) = exp(2i*pi*k/N).
%	alpha = 2*cos(2*pi*k/N).
%

w = exp(2i*pi*k/N);

%alpha
alpha = 2*cos(2*pi*k/N);

% Выход БИХ-фильтра с DSPLIB.ORG 
b_dsplib = [w, -1, 0];
a_dsplib = [1, -alpha, 1];
X = filter(b_dsplib, a_dsplib, s);
S_dsplib = X(end);


%% Алгоритм Герцеля из книги  
% А. Оппенгейм, Р. Шафер Цифровая обработка сигналов. Москва, Техносфера, 2006.
% страница 634 формула (9.9) рисунок 9.2
%
%              1 - W * z^-1
% H(z) = ---------------------------
%          1 - alpha * z^-1 + z^-2
%
%
%	W = W_N^(k) = exp(-2i*pi*k/N).
%	alpha = 2*cos(2*pi*k/N).
%

w = exp(-2i*pi*k/N);

%alpha
alpha = 2*cos(2*pi*k/N);

b_opp = [1, -w,     0];
a_opp = [1, -alpha, 1];

X = filter(b_opp, a_opp, s);
S_opp = X(end);


%%  PRINT 
fprintf('DFT ....................................%9.4f %+9.4f i\n', real(S_dft), imag(S_dft));
fprintf('dspib.org Algorithm.....................%9.4f %+9.4f i\n', real(S_dsplib), imag(S_dsplib));
fprintf('Oppenheim Algorithm.....................%9.4f %+9.4f i\n', real(S_opp), imag(S_opp));
		
		
								
# -*- coding: utf-8 -*-

import numpy as np
import scipy.signal as signal
						
# N-точечное ДПФ
N = 100
# Номер спектрального отсчета (бина)
k = 32

# Исходный сигнал и расчет его ДПФ 
t = np.linspace(0, N, N, endpoint = False, dtype = 'float64') 
s = np.sin(2*np.pi*t*k/N+np.pi/6)
X = np.fft.fft(s)
S = X[k]

"""
Алгоритм Герцеля с DSPLIB.ORG. Формула (12) и рисунок 2.

               W - z^-1
H(z) = ---------------------------
         1 - alpha * z^-1 + z^-2


W = W_N^(-k) = exp(2i*pi*k/N).
alpha = 2*cos(2*pi*k/N).
"""

w = np.exp(2j * np.pi * k/N);

#alpha
alpha = 2*np.cos(2*np.pi*k/N);

# Выход БИХ-фильтра с DSPLIB.ORG 
b0 = np.array([w, -1, 0],     dtype = 'complex128')
a0 = np.array([1, -alpha, 1], dtype = 'complex128')
X = signal.lfilter(b0, a0, s)
Y = X[N-1]

"""
Алгоритм Герцеля из книги  
А. Оппенгейм, Р. Шафер Цифровая обработка сигналов. Москва, Техносфера, 2006.
страница 634 формула (9.9) рисунок 9.2

             1 - W * z^-1
H(z) = ---------------------------
         1 - alpha * z^-1 + z^-2

W = W_N^(k) = exp(-2i*pi*k/N).
alpha = 2*cos(2*pi*k/N).
"""

w = np.exp(-2j * np.pi * k/N);

#alpha
alpha = 2*np.cos(2*np.pi*k/N);

b1 = np.array([1, -w,     0],     dtype = 'complex128')
a1 = np.array([1, -alpha, 1],     dtype = 'complex128')

X = signal.lfilter(b1, a1, s)
Z = X[N-1]


# PRINT 
print('FFT......................', S, 'модуль', np.abs(S))
print('DSPLIB.ORG...............', Y, 'модуль', np.abs(Y))
print('А. Оппенгейм.............', Z, 'модуль', np.abs(Z))
		
		
Также ошибочная структура приведена в справке пакета MATLAB (хотя встроенная функция пакета выдает верный результат), а также в других источниках.
Таким образом, спектральный отсчет $S(k)$ равен:
$$ S(k) = y_{N-1}(k) = W_N^{-k} \cdot v(N-1) - v(N -2), $$
(13)
где $v(r)-$промежуточные значения, которые рассчитываются итерационно:
$$ v(r) = s(r) + 2 \cdot \cos \left(2\pi\cdot\frac{k}{N} \right) \cdot v(r-1) - v(r-2). $$
(14)
Пример использования алгоритма Гёрцеля
Рассмотрим пример использования алгоритма Гёрцеля. Пусть $N=8$, исходный сигнал $s(n)$ показан на рисунке 3.
Рисунок 3. Исходный сигнал
Рассчитаем при помощи алгоритма Герцеля спектральный отсчет $S(1)$ с номером $k=1$.
Предварительно рассчитаем коэффициенты $\alpha$ и $W_N^{-k}$:
$$ \alpha = 2 \cdot \cos \left(2\pi\cdot\frac{1}{8} \right) = 1.4142; $$ $$ W_N^{-k} = W_8^{-1} = \exp\left( j \cdot \frac{2\pi}{8} \right) = 0.7071 + j \cdot 0.7071. $$
(15)
В качестве начальных условий расчета зададим $v(-2) = v(-1) = 0$.
Произведем итерационный расчет $v(r)$, $r = 0 \ldots N-1$ согласно (14):
\begin{align} v(0) =& s(0) + \alpha \cdot v(-1) -v(-2) = 3; \\ v(1) =& s(1) + \alpha \cdot v(0) - v(-1)= 6.2426;\\ v(2) =& s(2) + \alpha \cdot v(1) -v(0) = 6.8284; \\ v(3) =& s(3) + \alpha \cdot v(2) - v(1)= 2.4142;\\ v(4) =& s(4) + \alpha \cdot v(3) -v(2) = -2.4142; \\ v(5) =& s(5) + \alpha \cdot v(4) - v(3)= -7.8284;\\ v(6) =& s(6) + \alpha \cdot v(5) -v(4) = -11.6569; \\ v(7) =& s(7) + \alpha \cdot v(6) - v(5)= -10.6569. \end{align}
(16)
Тогда, согласно выражению (13) можно получить значения для спектрального отсчета $S(1)$:
$$ S(1) = W_8^{-1} \cdot v(7) - v(6) = 4.1213 - j \cdot 7.5355. $$
(17)
Использование алгоритма Гёрцеля для декодирования DTMF сигнала
Рисунок 4. Кодирование
символов DTMF
DTMF (англ. Dual-Tone Multi-Frequency) сигнал представляет собой двухтональный сигнал, который используется для набора номера цифровых АТС.
Каждый символ DTMF кодируется двумя тонами, как это показано на рисунке 4. Один тон берется из нижней группы частот, а второй из верхней группы.
Таким образом, для декодирования DTMF нам необходимо рассчитать восемь спектральных отсчетов ДПФ с индексами, соответствующими частотам восьми тональных сигналов. При частоте дискретизации $F_{\textrm{s}} = 8$ кГц, для $205-$точечного ДПФ ($N = 205$), номера спектральных отсчетов $k$ приведены в таблице 1.

$k$1820222431343842
$f(k)$, (Гц)6977708529411209133614771633
Таблица 1. Индексы спектральных отсчетов $k$, соответствующих тонам DTMF,
$F_{\textrm{s}} = 8$ кГц и $N = 205$
Таким образом, для декодирования DTMF нам необходимо расчитать всего 8 спектральных отсчетов вместо расчета $205-$точечного ДПФ.
Моделирование DTMF сигналов и расчет спектральных отсчетов производились с использованием библиотеки DSPL. Исходный код программы приведен в разделе примеры библиотеки DSPL.
Временные осциллограммы DTMF сигналов, полученные в результате работы программы моделирования, показаны на рисунке 5.
Рисунок 5. Осциллограммы сигналов DTMF
Амплитуды спектральных отсчетов, рассчитанные при помощи функции dspl_goertzel, для индексов, соответствующих частотам DTMF (смотри таблицу 1), показаны на рисунке 6.
Рисунок 6. Амплитуды спектральных отсчетов DTMF, рассчитанные при помощи алгоритма Гёрцеля
Пользователи Matlab и GNU Octave могут выполнить скрипт goertzel_dtmf.m для построения рисунков 5 и 6.
	
					
% Пример использования алгоритма Гёрцеля для анализа DTMF сигналов
clear all;
close all;
clc;

addpath('functions');

Fs = 8000;
N  = 205;
frq = [697, 770, 852, 941, 1209, 1336, 1477, 1633];
ind = [18,  20,  22,  24,  31,   34,   38,   42]+1;

t = (0:N-1) / Fs;

for(k = 1:4)  
  for(m = 1:4)
    s = sin(2*pi*frq(k)*t) + sin(2*pi*frq(m+4)*t);
    S = goertzel(s, ind);    
    
    figure(1); subplot(4,4,(k-1)*4+m); 
    plot(t,s); axis([0, 0.025, -2, 2]);
    
    figure(2); subplot(4,4,(k-1)*4+m); 
    stem(frq,abs(S), '.'); axis([500, 1800, 0, 130]);
  end 
end

	
function S = goertzel(s, ind)
% S = goertzel(s, ind)
% Function returns DFT samples indexes ind, calculated by Goertzel algorithm
% for input signal s.
%
% Input parameters
%  s   - input signal vector [N x 1];
%  ind - DFT samples indexes vector (indexation starts from 0)
%
% Ouptut parameters
%  S   - DFT samples corresponds to indexes ind
%
% Author: Sergey Bakhurin (dsplib.org)
%

  N = length(s);
  for k = 1:length(ind);
    
    w = exp(2i*pi*(ind(k)-1)/N);
    
    %alpha
    alpha = 2*cos(2*pi*(ind(k)-1)/N);
    
    % second-order IIR 
    b = [w, -1, 0];
    a = [1, -alpha, 1];
    X = filter(b, a, s);
    S(k) = X(end);
  end

Выводы
В данном разделе мы рассмотрели алгоритм Гёрцеля, используемый для расчета отдельных спектральных отсчетов при помощи БИХ-фильтра. Данный алгоритм позволяет достичь высокой производительности при расчете небольшого количества спектральных отсчетов. При этом алгоритм Гёрцеля уступает в производительности алгоритмам быстрого преобразования Фурье для расчета полного ДПФ.
Алгоритм Гёрцеля находит применение при декодировании DTMF сигналов тонового набора цифровой АТС, ввиду необходимости анализа амплитуд восьми спектральных отсчетов. В данном разделе приведен результат моделирования DTMF сигналов с использованием функции dspl_goertzel, входящей в состав библиотеки DSPL, а также m-скриптов Matlab и GNU Octave.

Ваши комментарии, вопросы и пожелания вы можете оставить на форуме или в гостевой книге.
Список литературы
[1] Goertzel Gerald An Algorithm for the Evaluation of Finite Trigonometric Series. The American Mathematical Monthly, Vol. 65, No. 1, Jan., 1958, pp. 34-35

[2] Оппенгейм А., Шафер Р. Цифровая обработка сигналов. Москва, Техносфера, 2006.

[3] Голд Б., Рэйдер Ч. Цифровая обрботка сигналов. С приложением работы Д. Кайзера "Цифровые фильтры" Москва, Советское радио, 1973.


Oбнаружили ошибку в тексте? Выделите ее мышкой и нажмите