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

Содержание

DSPL-2.0 — свободная библиотека алгоритмов цифровой обработки сигналов

Распространяется под лицензией LGPL v3

Страница библиотеки на GitHub

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

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

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

Рекуррентное соотношение для расчета фиксированного спектрального отсчета сигнала

Выражение для N-точечного дискретного преобразования Фурье сигнала s(n), n = 0 \ldots N-1 имеет вид:

equation 1
(1)
Поворотные коэффициенты W_N^{n k} обладают следующим свойством:

equation 2
(2)
Таким образом, умножение выражения (1) на W_N^{- k N} не приведет к изменению результата:

equation 3
(3)
Распишем сумму (3):
equation 4
(4)
Обозначим S(k) = y_{N-1}(k) для фиксированного номера k спектрального отсчета ДПФ. Также вынесем в (4) W_N^{-k} за скобки и получим:
equation 5
(5)
Выразим y_{N-2}(k) через y_{N-3}(k). Для этого снова вынесем W_N^{-k} за скобки и получим:
equation 6
(6)
Таким образом, мы получили рекуррентное соотношение для вычисления y_{r}(k) на любом шаге r = 0 \ldots N-1 через y_{r-1}(k), полученных на предыдущем шаге для фиксированного номера спектрального отсчета k:

equation 7
(7)
где s(r) — отсчет входного сигнала с номером r.

Данное рекуррентное соотношение на шаге с номером r = N-1 приводит нас к S(k) = y_{N-1}(k), т.е. на N-1 шаге получим спектральный отсчет с номером k.

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

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

Обозначим через X(z) и Y(z) z-образы x(r) и y_r(k) соответственно.

Тогда z^{-1} Y(z) это z-образ y_{r-1}(k), и разностное уравнение (7) в операторном виде равно

equation 8
(8)

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

equation 9
(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 z^{-1}:

equation 10
(10)
Рассмотрим более подробно сумму:
equation 11
(11)
Тогда (10) с учетом (11) можно представить следующим образом:
equation 12
(12)
Передаточная характеристика (12) соответствует БИХ-фильтру второго порядка, стуктурная схема которого показана на рисунке 2, где \alpha = 2 \cos  \left( 2\pi  \frac{k}{N}\right).

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

Обратим внимание, что \alpha = 2\cos \left( 2\pi  \frac{k}{N}\right) — вещественный коэффициент, соответственно умножения в рекурсивной ветви фильтра стали вещественными, причем умножения на -1 можно не учитывать. У нас остался один комплексный коэффициент W_N^{-k} в числителе передаточной характеристики[1]. Но мы выяснили, что промежуточные значения y_r(k) нам неважны, соответственно умножать на W_N^{-k} можно только на N-1 итерации, когда рассчитывается непосредственно S(k).

Таким образом, применение фильтра второго порядка позволяет рассчитать S(k) для фиксированного k с использованием N умножений, но не комплексных, а вещественных. Также требуется одно комплексное умножение на W_N^{-k} на последней итерации.

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

equation 13
(13)
где v(r)-промежуточные значения, которые рассчитываются итерационно:
equation 14
(14)

Пример использования алгоритма Гёрцеля

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

Исходный сигнал
Рисунок 3. Исходный сигнал

Рассчитаем при помощи алгоритма Герцеля спектральный отсчет S(1) с номером k=1.

Предварительно рассчитаем коэффициенты \alpha и W_N^{-k}:

equation 15
(15)
В качестве начальных условий расчета зададим v(-2) = v(-1) = 0.

Произведем итерационный расчет v(r), r = 0 \ldots N-1 согласно (14):

equation 16
(16)
equation 17
(17)
equation 18
(18)
equation 19
(19)
equation 20
(20)
equation 21
(21)
equation 22
(22)
equation 23
(23)
Тогда, согласно выражению (13) можно получить значения для спектрального отсчета S(1):
equation 24
(24)

Использование алгоритма Гёрцеля для декодирования DTMF сигнала

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

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

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

Кодирование символов DTMF
Рисунок 4. Кодирование символов DTMF

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

При частоте дискретизации F_{\textrm{s}} = 8 кГц, для 205-точечного ДПФ (N = 205), номера спектральных отсчетов k приведены в таблице 1.

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

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

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

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

Осциллограммы сигналов DTMF
Рисунок 5. Осциллограммы сигналов DTMF

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

Амплитуды спектральных отсчетов DTMF, рассчитанные при помощи алгоритма Гёрцеля
Рисунок 6. Амплитуды спектральных отсчетов DTMF, рассчитанные при помощи алгоритма Гёрцеля

Выводы

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

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

Программная реализация в библиотеке DSPL

Данные для построения рисунков данного раздела были просчитаны при использовании библиотеки 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 выдает верный результат.

Список литературы
[1] Goertzel G. An Algorithm for the Evaluation of Finite Trigonometric Series. The American Mathematical Monthly, Vol. 65, No. 1, Jan., 1958, pp. 34-35.

[2] Оппенгейм А., Шаффер Р. Цифровая обработка сигналов. Москва, Техносфера, 2012. 1048 с. ISBN 978-5-94836-329-5

[3] Блейхут Р. Быстрые алгоритмы цифровой обработки сигналов. Москва, Мир, 1989, 448 c.

Последнее изменение страницы: 14.11.2020 (11:20:43)
Страница создана Latex to HTML translator ver. 5.20.11.14