|
|
- Алгоритм Гертцеля
Пусть наша задача состоит в вычислении параметров отдельных спектральных составляющих входного сигнала, причем выборка сигнала представляет собой массив вещественных чисел. Для решения этой задачи существует так называемый алгоритм Гертцеля. Описать этот алгоритм можно так:
(9) Sk[n] = x[n] + Mk*Sk[n-1] – Sk[n-2]
(10) Sk[-1]=Sk[-2]=0
(11) Mk=2*cos(2*p*k/N)
(12) Yk[n]=Sk[n]- WkN *Sk[n-1]
Как вы помните,
(13) WkN= exp(-j*2*p*k/N).
Вычисление спектра с помощью алгоритма Гертцеля начинается с задания начальных условий (10) и заполнения таблицы коэффициентов Mk (11). Здесь индекс “k” соответствует номеру вычисляемой гармоники. Далее, используя отсчет входного сигнала x[n] при n=0, вычисляют Sk[0] (9). Теперь по формуле (12) можно вычислить Yk[n]. Все значения Yk[n] при этом считать не нужно, так как амплитуде интересующей нас гармоники соответствует лишь значение Yk[N-1] (почему происходит именно так, в двух словах объяснить не берусь). Таким образом, для вычисления амплитуды некоторой гармоники нужно сначала найти N значений коэффициента Sk[n], последним из которых будет Sk[N-1]. Для этого нам потребуется N умножений на Mk=2*cos(2*p*k/N). Последний штрих – с помощью (12) нужно вычислить Yk[n], а это еще одно умножение. Итак, для расчета амплитуды одной гармоники по алгоритму Гертцеля нужно примерно N+1 умножений, в то же время ДПФ для решения этой же задачи потребуется около N*2 умножений. Поэтому алгоритм Гертцеля в данном случае более эффективен. Сразу заметим, что применять рассматриваемый алгоритм для расчета параметров всех N/2 гармоник бессмысленно, ведь БПФ вычисляется гораздо быстрее. Алгоритм Гертцеля хорош там, где нужно “выловить” из входного сигнала лишь несколько спектральных составляющих. Такая задача возникает, например, при декодировании сигналов тонального набора телефонного номера (DTMF). Правда следует заметить, что данный алгоритм требует более-менее мощного вычислительного устройства. Приведу пример: Пусть у нас есть 84 отсчета входного сигнала, и пусть мы хотим оценить амплитуды всех 42-х гармоник. Процессор 80386 DX-40 справился с этой задачей примерно за 50 миллисекунд. Поэтому для однокристальных микроконтроллеров этот алгоритм, скорее всего, не очень подходит.
В заключение приведу простейшую программу, написанную на Turbo Pascal 7.0. Эта программа рассчитывает амплитуды гармоник некоторого сигнала по его выборке и выводит их на экран. Сама выборка перед расчетом помещается в файл Viborka.bin:
{--- Вычисление спектра с использованием алгоритма Гертцеля. Этот алгоритм полезен, если нужно вычислить спектр только в нескольких точках. Так, при вычислении с помощью ДПФ спектра в одной точке требуется примерно N*2 умножений. При вычислении же с помощью алгоритма Гертцеля нужно всего лишь N+1 умножений. Для вычисления спектра в каждой точке алгоритм использовать бессмысленно, гораздо лучше использовать БПФ. ---}
program Goertzel_Algorithm;
uses crt;
{--- "N_Big" - количество входных отсчетов ---}
const n_big = 16;
pii = 3.1415926535897932385;
{--- "How_Many_Garmonics" - количество вычисляемых гармоник спектра ---}
how_many_garmonics = n_big div 2;
var viborka : file of byte;
entries : array [0..n_big - 1] of real;
s_array : array [-2..n_big - 1] of real;
m_coeff : array [0..n_big - 1] of real;
amplitudes : array [0..how_many_garmonics] of real;
{--- "K" - это номер вычисляемой гармоники; ---}
{--- "N" - это номер отсчета входного сигнала ---}
i,j,k : integer;
transformer : byte;
{--- Самая главная процедура в этой программе. Процедура вычисляет спектр входного сигнала, используя для этого алгоритм Гертцеля. В результате получаем массив из "N_Big" элементов, первая половина которого соответствует амплитудам гармоник входного сигнала ---}
procedure goertzel;
begin
for k:=0 to how_many_garmonics do begin
for i:=0 to n_big - 1 do
s_array[i]:=entries[i] + (m_coeff[k] * s_array[i-1]) -
s_array[i-2];
amplitudes[k]:=sqrt(s_array[n_big-1]*s_array[n_big-1]+
s_array[n_big - 2]*s_array[n_big - 2] -
m_coeff[k]*s_array[n_big - 1]*
s_array[n_big - 2]);
end;
end;
{--- Процедура формирует массив элементов M(K) для выполнения алгоритма Гертцеля ---}
procedure m_initial;
begin
k:=0;
for k:=0 to n_big - 1 do
m_coeff[k]:=2*cos(2*pii*k/n_big);
end;
begin
clrscr;
{--- Заполним таблицу множителей M(K) ---}
m_initial;
{--- Заполним первые два элемента S[-2] и S[-1], иначе алгоритм не сможет стартовать ---}
s_array[-2]:=0;
s_array[-1]:=0;
{***Сопоставим файлу переменную VIBORKA***}
assign (viborka,'viborka.bin');
{***Проверим наличие ошибки***}
if ioresult<>0 then writeln ('Ошибка при открытии файла !');
{***Откроем файл. Если файл не существует, мы получим ошибку***}
reset (viborka);
{--- Получим выборку входного сигнала ---}
for i:=0 to n_big - 1 do begin
read (viborka,transformer);
entries[i]:=transformer;
end;
{--- Теперь можно заняться вычислением спектра: ---}
goertzel;
{--- Выведем результаты расчетов на экран ---}
writeln ('Амплитуды гармоник равны:');
for k:=0 to how_many_garmonics do
writeln (amplitudes[k]:5:3);
writeln ('Нажмите <ПРОБЕЛ> ...');
while readkey<>' ' do;
end.
E-mail: 
info@telesys.ru
