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

Введение Рекуррентное соотношение для расчета фиксированного спектрального отсчета сигнала Алгоритм Гёрцеля Пример использования алгоритма Гёрцеля Использование алгоритма Гёрцеля для декодирования DTMF сигнала Выводы Программная реализация в библиотеке DSPL Приложения Примечания Список литературы

DSPL-2.0 — свободная библиотека алгоритмов цифровой обработки сигналов Распространяется под лицензией LGPL v3 Страница проекта на SourceForge

Обнаружили ошибку? Выделите ее мышью и нажмите

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

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

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

Данное рекуррентное соотношение на шаге с номером приводит нас к , т.е. на шаге получим спектральный отсчет с номером .

Проанализировав (7), обратим внимание, что рекуррентное соотношение можно трактовать как разностное уравнение БИХ-фильтра первого порядка с комплексным коэффициентом , структурная схема которого приведена на рисунке 1.

Рисунок 1. Структурная схема БИХ-фильтра, реализующего расчет спектрального отсчета с номером

Обозначим через и -образы и соответственно.

Тогда это -образ и разностное уравнение (7) в операторном виде равно

Передаточная характеристика полученного БИХ-фильтра имеет вид:

Таким образом, мы зафиксировали номер спектрального отсчета и преобразовали выражение ДПФ для одного фиксированного спектрального отсчета к рекуррентному соотношению, которое позволяет нам на шаге получить значение искомого -го спектрального отсчета . При этом все промежуточные значения рекуррентного соотношения нас не интересуют, а интересует только

Затем мы трактовали рекуррентное соотношение как разностное уравнение БИХ-фильтра с передаточной характеристикой (9).

В результате мы получили фильтр первого порядка с комплексным коэффициентом применение которого на -ой итерации на выходе дает .

При фиксированном значение комплексного коэффициента всегда постоянно.

Для вычисления одного спектрального отсчета при использовании БИХ-фильтра требуется комплексных умножений и сложений, что не дает никаких преимуществ по сравнению с прямым вычислением ДПФ, согласно выражению (1).

Однако сокращение вычислений можно получить, если домножить числитель и знаменатель передаточной характеристики БИХ-фильтра (9) на :

Рисунок 2. Структурная схема БИХ-фильтра реализующего расчет спектрального отсчета с номером

Обратим внимание, что  — вещественный коэффициент, соответственно умножения в рекурсивной ветви фильтра стали вещественными, причем умножения на можно не учитывать. У нас остался один комплексный коэффициент в числителе передаточной характеристики[1]. Но мы выяснили, что промежуточные значения нам неважны, соответственно умножать на можно только на итерации, когда рассчитывается непосредственно .

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

Таким образом, спектральный отсчет равен:

Рассмотрим пример использования алгоритма Гёрцеля. Пусть , исходный сигнал показан на рисунке 3.

Рассчитаем при помощи алгоритма Герцеля спектральный отсчет с номером .

Предварительно рассчитаем коэффициенты и :

Произведем итерационный расчет , согласно (14):

DTMF (англ. Dual-Tone Multi-Frequency) сигнал представляет собой двухтональный сигнал, который используется для набора номера цифровых АТС.

Каждый символ DTMF кодируется двумя тонами, как это показано на рисунке 4.

Один тон берется из нижней группы частот, а второй из верхней группы.

Таким образом, для декодирования DTMF нам необходимо рассчитать восемь спектральных отсчетов ДПФ с индексами, соответствующими частотам восьми тональных сигналов.

При частоте дискретизации кГц, для -точечного ДПФ (), номера спектральных отсчетов приведены в таблице 1.

Таблица 1. Индексы спектральных отсчетов DTMF для -точечного ДПФ

Таким образом, для декодирования DTMF нам необходимо расчитать всего 8 спектральных отсчетов вместо расчета -точечного ДПФ.

Моделирование DTMF сигналов и расчет спектральных отсчетов производились с использованием библиотеки DSPL, в которой реализованы функции goertzel и goertzel_cmplx для анализа вещественных и комплексных данных. Исходный код программы симуляции сигналов DTMF и расчета спектральных компонент согласно таблицы 1 приведен в файле goertzel_dtmf.c .

Временные осциллограммы DTMF сигналов, полученные в результате работы программы моделирования goertzel_dtmf.c , показаны на рисунке 5.

Амплитуды спектральных отсчетов, рассчитанные при помощи функции алгоримта Герцеля, для индексов, соответствующих частотам DTMF (смотри таблицу 1), показаны на рисунке 6.

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

Алгоритм Гёрцеля находит применение при декодировании DTMF сигналов тонового набора цифровой АТС, ввиду необходимости анализа амплитуд восьми спектральных отсчетов. В данном разделе приведен результат моделирования DTMF сигналов с использованием функции написанных на языке Python, а также m-скриптов Matlab и GNU Octave.

Данные для построения рисунков данного раздела были просчитаны при использовании библиотеки DSPL-2.0

Ниже приведён исходный код программы расчета данных для построения рисунков 5 и 6 :

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include "dspl.h"



#define FS   8000.0
#define N    205




void scale_data(double* x, double* y, int n, 
                double xmin, double xmax,
                double ymin, double ymax,
                double w, double h)
{
  double kx = w / (xmax - xmin);
  double ky = h / (ymax - ymin);
  
  int k;
  
  for(k = 0; k < n; k++)
  {
    x[k] = (x[k] - xmin) * kx;
    y[k] = (y[k] - ymin) * ky;
  }
}


int main(int argc, char* argv[])
{
  double t[N];         /* время                      */
  double s[N];         /* исходный сигнал            */
  void* hdspl; /* DSPL handle        */
  void* hplot; /* GNUPLOT handle     */
  
  complex_t cS[N];
  double    rS[N];
  
  double frq[8]  = {697.0, 770.0, 852.0, 941.0, 1209.0, 1336.0, 1477.0, 1633.0}; 
  int    ind[8]=  {    18,    20,    22,    24,     31,     34,     38,     42};
  char   smb[16] = "123a456b789cs0rd";
  double frq8[8];
  int k, m, n;
  char fname[64];  /* имя файла для сохранения рассчитанных данных */


  /* Загружаем libdspl.dll */
  hdspl = dspl_load();
  if(!hdspl)
  {
    printf("Cannot to load libdspl!\n");
    return 0;
  }
  
  
  for(k = 0; k < 4; k++)
  {
    for(m = 0; m < 4; m++)
    {
      linspace(0, (double)N/FS, N, DSPL_PERIODIC, t);
      for(n = 0; n < N; n++)
      {
        s[n] = sin(M_2PI * frq[k] * t[n]) + sin(M_2PI * frq[m+4] * t[n]);
      }
      goertzel(s, N, ind, 8, cS);
      for(n = 0; n < 8; n++)
        rS[n] = ABS(cS[n]);
      scale_data(t, s, N, 0.0, t[N-1], -2.0, 2.0, 4.5, 2.0);
      memcpy(frq8, frq, 8*sizeof(double));
      scale_data(frq8, rS, 8, 500.0, 1800.0, 0.0, 120.0, 4.5, 2.0);
      sprintf(fname, "dat/dtmf_sym_%c_time.txt", smb[k*4+m]);
      writetxt(t,s,N,fname);
      sprintf(fname, "dat/dtmf_sym_%c_freq.txt", smb[k*4+m]);
      writetxt(frq8,rS,8,fname);
      
    }  
  }
  
  


  /* Выгружаем libdspl.dll из памяти */
  dspl_free(hdspl);

  return 0;
}

Проверочный скрипт goertzel_test.m для проверки структурной схемы в системе Matlab или в GNU Octave:

		
						
% 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));
		

Проверочный скрипт goertzel_test.py для проверки структурной схемы в Python3:

# -*- 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))

[1] Внимание! В литературе приведена структурная схема, отличная от рисунка 2. Проверка показала, что приведенный в литературе фильтр второго порядка выдает ошибочный результат на выходе. Фильтр, приведенный на рисунке 2, выдает точное значения отсчета ДПФ. Вы можете самостоятельно удостовериться в правильности результата фильтра, представленного на рисунке  2, выполнив проверочный скрипт goertzel_test.m в системе Matlab или в GNU Octave, или goertzel_test.py на языке Python. Также ошибочная структура приведена в справке пакета MATLAB, а также в других источниках, хотя встроенная функция пакета MATLAB выдает верный результат.