Юкио Сато "Обработка сигналов", Подозреваю ошибку в алгоритме БПФ Опции

[ZhekSooN]

Всем добрый день.
Решил скопировать алгоритм БПФ из книги Юкио Сато "Обработка сигналов" как самый простой и прозрачный из всех мною виданных.
В книге он написан на языке Basic, я же переписал его на Pascal`е. Переписал точно, ошибок вроде не нашел, но.. алгоритм не работает. В тестовой программе на выходе получается бред вместо исходной волны. Может ли кто-нибудь хорошо понимающий в алгоритмах БПФ сказать, где же ошибка в книге и/или у меня?
Прикрепляю страницу из книги и свой исходник

Эскизы прикрепленных изображений

 

Прикрепленные файлы

 FFT.zip ( 929 байт ) Кол-во скачиваний: 44

 

[ivan219]

Это?

Код

Type
TReal1DArray = array of Double;

procedure FastFourierTransform(var a : TReal1DArray; nn : Integer;
InverseFFT : Boolean); // БПФ комплексной функции

Быстрое преобразование Фурье

Алгоритм проводит быстрое преобразование Фурье комплексной
функции, заданной nn отсчетами на действительной оси.

В зависимости от переданных параметров, может выполняться
как прямое, так и обратное преобразование.

Входные параметры:
nn - Число значений функции. Должно быть степенью
   двойки. Алгоритм не проверяет правильность
   переданного значения.
a - array [0 .. 2*nn-1] of Real
   Значения функции. I-ому значению соответствуют
   элементы a[2*I]  (вещественная  часть)
   и a[2*I+1] (мнимая часть).
InverseFFT - направление преобразования. True, если обратное, False, если прямое.

procedure FastFourierTransform(var a : TReal1DArray; nn : Integer;
InverseFFT : Boolean);
var
  Jj, n ,Mmax ,m , j, istep, i, isign: Integer;
  wtemp, wr, wpr, wpi, wi, theta, tempr, tempi: Double;
begin
if InverseFFT then isign := -1
else isign := 1;

n := 2 * nn;
j := 1;
I := 1;
mmax := 2;

while I <= n do   // Реверс
  begin
   if j > i then
    begin
     tempr := a[j - 1];
     tempi := a[j];
     a[j - 1] := a[i - 1];
     a[j] := a[i];
     a[i - 1] := tempr;
     a[i] := tempi;
    end;
   m := nn;
   while (m >= 2) and (j > m) do
    begin
     j := j - m;
     m := m div 2;
    end;
   j := j + m;
   I := I + 2;
  end;              // Реверс

while n > mmax do   // FFT цикл mmax = 2, 4, 8, 16..2 * NN
  begin
   istep := 2 * mmax;                 // 4, 8, 16..2 * NN
   theta := 2 * Pi * isign  / mmax;   // Прямое 2 * Pi * (-1) / (2, 4, 8, 16..2 * NN) Обратное 2 * Pi * 1 / (2, 4, 8, 16..2 * NN)
   wpr := -2 * sqr(sin(0.5 * theta)); // Cos(X) - 1
   wpi := sin(theta);                 // Sin(X)
   wr := 1;
   wi := 0;
   M := 1;
   while M <= mmax do
    begin
     for Jj := 0 to (N - M) div istep do
      begin
       i := m + Jj * istep;
       j := i + mmax;
       tempr := wr * a[j - 1] - wi * a[j];
       tempi := wr * a[j] + wi * a[j - 1];
       a[j - 1] := a[i - 1] - tempr;
       a[j] := a[i] - tempi;
       a[i - 1] := a[i - 1] + tempr;
       a[i] := a[i] + tempi;
      end;
     wtemp := wr;
     wr := wr * wpr - wi * wpi + wr;
     wi := wi * wpr + wtemp * wpi + wi;
     M := M + 2;
    end;
   mmax := istep;
  end;                   //FFT

if InverseFFT then      // Обратное FFT
  for I := 0 to N - 1 do
   a[I] := a[I] / nn;
end;

P.S. единственно что замени строки

Код

wpr := -2 * sqr(sin(0.5 * theta)); // Cos(X) - 1
wr := wr * wpr - wi * wpi + wr;
wi := wi * wpr + wtemp * wpi + wi;

на

Код

wpr := Cos(theta); // Cos(X)
wr := wr * wpr - wi * wpi;
wi := wi * wpr + wtemp * wpi;

Тогда получится прозрачней некуда.