[ZhekSooN 0]

Всем добрый день.

Решил скопировать алгоритм БПФ из книги Юкио Сато "Обработка сигналов" как самый простой и прозрачный из всех мною виданных.

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

Прикрепляю страницу из книги и свой исходник

FFT.zip

[maior 0]

В свое время сделал то же самое с алгоритмом из "Рабинер-Гоулд" - все заработало сразу, в разных видах, очень быстро и было очень компактно по коду.

[ZhekSooN 0]

maior, не поможите ли с кодом? Может ваш сохранился...

[ivan219 0]

http://alglib.sources.ru/fasttransforms/

[ZhekSooN 0]

http://alglib.sources.ru/fasttransforms/

Спасибо, алглибовские сырцы уже смотрел, но в первую очередь инетерсуюсь прозрачными "монолитными" алгоритмами без лишнего мусора в виде кучи функций и заточки под определенный язык. :smile3046:

[ivan219 0]

Это?

 

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;

 

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

[bahurin 0]

Спасибо, алглибовские сырцы уже смотрел, но в первую очередь инетерсуюсь прозрачными "монолитными" алгоритмами без лишнего мусора в виде кучи функций и заточки под определенный язык. :smile3046:

 

если вам нужно использовать бпф в своем приложении возьмите бибилиотеку там есть БПФ. Если хотите потренироваться, то флаг вам в руки. Когда напишете алгоритм опять таки возьмите указанную библиотеку и сравните быстродействие если ваш алгоритм будет быстрее считайте, что вы добились успехов в написании БПФ.

 

PS в библиотеке алгоритм кули тьюки по основанию 2. Все у кого есть свои собственные функции БПФ предлагаю сравнить быстродействие с указанной библиотекой

Изменено 29 августа, 2010 пользователем bahurin

[ZhekSooN 0]

если вам нужно использовать бпф в своем приложении возьмите бибилиотеку там есть БПФ. Если хотите потренироваться, то флаг вам в руки. Когда напишете алгоритм опять таки возьмите указанную библиотеку и сравните быстродействие если ваш алгоритм будет быстрее считайте, что вы добились успехов в написании БПФ.

 

PS в библиотеке алгоритм кули тьюки по основанию 2. Все у кого есть свои собственные функции БПФ предлагаю сравнить быстродействие с указанной библиотекой

Ну, в первую очередь меня интересует алгоритм, а не результат. Во-вторых, он будет переписываться на ассембелере для AVR и экстримально оптимизироваться, так что мне ну никак не подходят реализации с любой оптимизацией (как, например, от ivan219) под конкретные условия (в алгоритме ivan219 задумывалось, что взятие синуса будет долгим, вот и оптимизировалось). В алгоритме от Сато нету ни капли лишней воды, но он не работает :( Вот поэтому и пытаюсь найти ошибку, а не беру готовую либу (кстати, такие есть и для AVR, но хочется сделать лучше :) ).

[ZhekSooN 0]

Урррааа!Я нашел проблему - она крылась в коде двоично-инверсной перестановки. Поменял этот участок на заведомо рабочий - все заработало как надо! Благодарю всех за участие в обсуждении.

[L++ 0]

Урррааа!Я нашел проблему - она крылась в коде двоично-инверсной перестановки. Поменял этот участок на заведомо рабочий - все заработало как надо! Благодарю всех за участие в обсуждении.

 

в строке 300 должно быть : J=J-N2 и т.д.