ru.algorithms

 
 - RU.ALGORITHMS ----------------------------------------------------------------
 From : Evgenij Masherov                     2:5020/175.2   17 May 2001  15:02:57
 To : Vladimir Siltchenko
 Subject : БПФ
 -------------------------------------------------------------------------------- 
 

 Thu May 17 2001 08:57, Vladimir Siltchenko wrote to All:
  VS> Сорри за назойливость, но не мог бы ктонибудь поделиться толковым
  VS> описанием(+теория)  алгоритма сабжа? Ссылки только приветствуются.
 
  VS>                 C уважением, Vladimir Siltchenko.
 
 //БПФ для 512 точек.
 
 procedure fft512(var Ar:array of double; var Ai:array of double);
 // Действительные и мнимые части коэффициентов хранятся в отдельных массивах.
 var
 Le,Le1,jt,it,ip,lt,kt:Integer;
 Tr,Ti,Ur,Ui,Wr,Wi:Double;
 
 begin
 // Основная идея БПФ состоит в том, что если прямой расчет 
 // требует квадратичного по отношению к размерности исходных массивов
 // времени, то для половины длины вектора нужно вчетверо меньше операций.
 // Последовательно деля массив пополам, причем полученные на каждом шаге
 // вектора коэффициентов объединяются в единый домножением на (комплексный)
 // множитель поворота и сложением/вычитанием, обнаруживаем, что делений таких
 // у нас логарифм от длины вектора, а на каждом шаге число операций
 // пропорционально оной длине.
 // (Можно делить не на два, а на три, четыре и т.п. - в зависимости от
 // этого можно придти к иным алгоритмам, скажем, алгоритм Винограда
 // делит на части, являющиеся простыми числами.
 // В данном алгоритме мы делим отрезок на четные и нечетные отсчеты -
 // при этом целесообразно точки предварительно переставить, иначе
 // перепутаются частоты
 
 // А вот и перестановка
      for it:=0 to 511 do
      begin
      kt:=256;
      jt:=0;
      Le:=1;
      Le1:=it;
      for lt:=1 to 9 do
      begin
       if kt<=Le1 then begin
       jt:=jt+Le;
       Le1:=Le1-kt;
       end;
       Le:=Le*2;
       kt:=kt div 2;
      end;
      if it<jt then begin
         Tr:=Ar[jt];Ti:=Ai[jt];
         Ar[jt]:=Ar[it];Ai[jt]:=Ai[it];
         Ar[it]:=Tr;Ai[it]:=Ti;
         end;
      end;
 // Порядок точек - инверсный битовый (т.е. если записать номер числа 
 // в исходном векторе в обратном порядке битов - получим номер в выходном
 
 // Основной цикл
      for lt:=1 to 9 do
 // Hомер очередного деления
      begin
       Le:=1;
       for jt:=1 to lt do Le:=Le*2;
       Le1:=Le div 2;
       Ur:=1.0;
       Ui:=0.0;
 // Поворачивающий множитель (начальное значение)
       Wr:=cos(Pi/Le1);
       Wi:=sin(Pi/Le1);
       for jt:=0 to Le1-1 do
       begin
        it:=jt;
        while it<512-Le1 do
        begin
         ip:=it+Le1;
 // "Бабочка" - умножаем на поворачиватель и складываем/вычитаем
         Tr:=Ar[ip]*Ur-Ai[ip]*Ui;
         Ti:=Ar[ip]*Ui+Ai[ip]*Ur;
         Ar[ip]:=Ar[it]-Tr;
         Ai[ip]:=Ai[it]-Ti;
         Ar[it]:=Ar[it]+Tr;
         Ai[it]:=Ai[it]+Ti;
         it:=it+Le;
        end;
 
 // Hовое значение поворачивателя
       Tr:=Ur*Wr-Ui*Wi;
       Ti:=Ur*Wi+Ui*Wr;
       Ur:=Tr;
       Ui:=Ti;
       end;
      end;
 end;
 
 // Результат записан на месте исходного комплексного вектора
 // Для обратного преобразования - воспользуйтесь тем, что
 // обратное БПФ есть комплексно сопряженное БПФ от комплексно сопряженного
 // исходного вектора
 // Примечание: нормировка не проводилась, так что обратное БПФ от БПФ будет
 // в N раз больше исходного вектора. Hе забывайте делить!
 
 Евгений Машеров АКА СанитарЖеня
 
 --- ifmail v.2.15
  * Origin: FidoNet Online - http://www.fido-online.com (2:5020/175.2)
 
 

Вернуться к списку тем, сортированных по: возрастание даты  уменьшение даты  тема  автор