unit Fourier;
{=================================================================}
{  Версия файла 3.0.000                                           }
{  от 4 декабря 2002 г.                                           }
{==========================} interface {==========================}

uses GlType;

  //*************************************************************************//
  //Процедура быстрого преобразования Фурье (алгоритм Кули-Тьюки)            //
  //Data - комплексный массив исходных данных                                //
  //ISign - направление преобразования                                       //
  // sgn_Fw  =  1 - прямое преобразование                                    //
  // sgn_Rew = -1 - обратное преобразование                                  //
  //*************************************************************************//
  // Примечание: результаты обратного преобразования Фурье дожны быть        //
  // поделены на количество данных NN для того чтобьы восстановить           //
  // исходный комплексный массив данных !!!!                                 //
  //*************************************************************************//

procedure fft(Data:PComplexArr;NN:integer;ISign:integer);

 //Процедура прямого преобразования Фурье                                  //
 // двух действительных последовательностей                                //
 //************************************************************************//
 //Data - массив входных сигналов , Re - первый сигнал ,Im - второй сигнал //
 //X - преобразование Фурье первого сигнала                                //
 //Y - преобразование Фурье второго сигнала                                //
 //************************************************************************//
procedure DoubleFFT(Data,X,Y:PComplexArr;N:integer);

const
 sgn_Fw=1;
 sgn_Rew=-1;

implementation

procedure fft;
 var
  temp,wp,w:TComplex;
  i,j,m,mmax,istep,ii,jj:integer;
  theta,wtemp:double;
begin
   j := 0;
   FOR i := 0 to nn-1 DO BEGIN
      IF j > i THEN BEGIN
       temp:=Data^[j];
       Data^[j]:=Data^[i];
       Data^[i]:=temp
       END;
      m := nn div 2;
      WHILE (m >= 1) AND (j+1 > m)  DO BEGIN
       j := j-m;
       m:=m div 2;
       END;
      j := j+m
   END;
 mmax := 1;
   WHILE (nn > mmax) DO BEGIN
      istep:=2*mmax;
      theta:=Pi/(isign*mmax);
      wp.Re:=-2.0*sqr(sin(0.5*theta));
      wp.Im:=sin(theta);
      w.Re:=1;
      w.Im:=0;
      FOR ii := 0 to mmax-1 DO BEGIN
         FOR jj := ((2*(nn-ii)-1) DIV (istep*2)) DOWNTO 0 DO BEGIN  (***)
            i := ii+jj*istep;
            j := i+mmax;
            temp.Re := w.Re*data^[j].Re-w.Im*data^[j].Im;
            temp.Im := w.Re*data^[j].Im+w.Im*data^[j].Re;
            data^[j].Re := data^[i].Re-temp.Re;
            data^[j].Im := data^[i].Im-temp.Im;
            data^[i].Re := data^[i].Re+temp.Re;
            data^[i].Im := data^[i].Im+temp.Im
         END;
         wtemp:=w.Re;
         w.Re:=w.Re*wp.Re-w.Im*wp.Im+w.Re;
         w.Im:=w.Im*wp.Re+wtemp*wp.Im+w.Im
      END;
      mmax := istep
   END;
end;

procedure DoubleFFT;
 var i:integer;
  Z,Z1:TComplex;
begin
 fft(Data,N,sgn_Fw);
 X^[0].Re:=Data^[0].Re;X^[0].Im:=0;
 Y^[0].Re:=Data^[0].Im;Y^[0].Im:=0;
 for i:=1 to N-1 do begin
  Z:=Data^[i];
  Z1:=Data^[N-i];
  X^[i].Re:=0.5*(Z.Re+Z1.Re);
  X^[i].Im:=0.5*(Z.Im-Z1.Im);
  Y^[i].Re:=0.5*(Z.Im+Z1.Im);
  Y^[i].im:=0.5*(Z1.Re-Z.Re)
 end
end;

end.