FFT

[=?UNKNOWN?Q?Kov=E1cs_Andr=E1s?=]

> Kerem , akinek van, kuldjon nekem egy FFT algoritmust (vagy linket),
> legjobb Pascal-ban lenne de minden mas megoldas is erdekel.

Ez Delphiben van, felig meddig en irtam, ha ilyesmire gondoltal,
tudok hozza nemi magyarazatot fuzni, ha talaltal jobbat akkor
nem nezek utana.



unit ufft;

interface
uses math;
type

float=single;
tdata=array[0..1000000] of float;

pdata=^tdata;

tkompl=object
 x,y:float;
 function absval:float;
 function fi:float;
end;


tcdata=array[0..1000000] of tkompl;
pcdata=^tcdata;
function fft(adat:pdata;n:integer):pcdata;
implementation
function wn(n,k:integer):tkompl;
 begin
  with result do begin
   x:=cos(2*pi*k/n);
   y:=-sin(2*pi*k/n);
  end;
 end;

function fft(adat:pdata;n:integer):pcdata;
 var eredm,eprs,eptl:pcdata;
   aprs,aptl:pdata;
   i:integer;
 begin
   eredm:=nil;
   getmem(eredm,sizeof(tkompl)*n);


   if (n>2) then begin
       eprs:=nil;
       getmem(eprs,sizeof(tkompl)*n);
       eptl:=nil;
       getmem(eptl,sizeof(tkompl)*n);
       getmem(aprs,sizeof(float)*n);
       getmem(aptl,sizeof(float)*n);
       for i:=0 to n div 2-1 do begin
          aprs[i]:=adat[2*i];
   aptl[i]:=adat[2*i+1];
       end;
       freemem(adat);
       eprs:=fft(aprs,n div 2);
       eptl:=fft(aptl,n div 2);
       for i:=0 to n div 2-1 do begin
   eredm[i].x:=eprs[i].x+wn(n,i).x*eptl[i].x-wn(n,i).y*eptl[i].y;
   eredm[i].y:=eprs[i].y+wn(n,i).x*eptl[i].y+wn(n,i).y*eptl[i].x;
   eredm[i+n div 2].x:=eprs[i].x-wn(n,i).x*eptl[i].x+wn(n,i).y*eptl[i].y;
   eredm[i+n div 2].y:=eprs[i].y-wn(n,i).x*eptl[i].y-wn(n,i).y*eptl[i].x;
       end;
       freemem(eprs);
       freemem(eptl);
   end else begin
     eredm[0].x:=adat[0]+adat[1];
 eredm[1].x:=adat[0]-adat[1];
 eredm[0].y:=0;
 eredm[1].y:=0;
 freemem(adat);
   end;
   result:=eredm;
 end;
{ tkompl }

function tkompl.absval: float;
begin
 result:=sqrt(sqr(x)+sqr(y));
end;

function tkompl.fi: float;
begin
 result:=ArcTan2(y,x);
end;

end.