(*  FFT COMPUTES THE (FAST) FOURIER TRANSFORM OF THE VECTOR X
    (A COMPLEX ARRAY OF DIMENSION N). SOURCE: Ferziger; Numerical
        methods for engineering applications.
 
    X = DATA TO BE TRANSFORMED; ON RETURN IT CONTAINS THE TRANSFORM.
    N = SIZE OF VECTOR. MUST BE A POWER OF 2 (<4294967296).
    K = 1 FOR FORWARD TRANSFORM.
    K = -1 FOR INVERSE TRANSFORM. *)

unit fourier;

interface 

uses ucomplex;

procedure fft(var x:array of complex;n,k:longint);

implementation

function log2n(n:longint):longint;
var m,y:longint;
begin
    m:=1;
    y:=0;
    while(m<n) do begin
        m+=m; y+=1
    end;
    if(m<>n) then begin
        writeln(stderr,'N= ',n,' is not a power of 2.');
        halt(1);
    end;
    log2n:=y
end;

procedure fft(var x:array of complex;n,k:longint);
var u:array [0..31] of complex;
    w,v,t:complex;
    gain:double;
    l2n:longint;
    ii,jj,stage,sby2,s:longint;
    p,q:longint;
    jndex,index,itemp:longint;
begin
    l2n:=log2n(n);
    gain:=1.0/n;
    u[0]:=cexp(-k*2.0*pi*i*gain);
    for ii:=1 to l2n-1 do u[ii]:=u[ii-1]*u[ii-1];
    sby2:=n;
    for stage:=0 to l2n-1 do begin
        w:=1.0;
        v:=u[stage];
        s:=sby2;
        sby2:=sby2 div 2;
        for jj:=0 to sby2-1 do begin
            ii:=0;
            while(ii<n) do begin
                p:=ii+jj;
                q:=p+sby2;
                t:=x[p]+x[q];
                x[q]:=(x[p]-x[q])*w;
                x[p]:=t;
                ii+=s
            end;
            w*=v;
        end
    end;
    for ii:=0 to n-1 do begin
        index:=ii;
        jndex:=0;
        for jj:=0 to l2n-1 do begin
            jndex+=jndex;
            itemp:=index div 2;
            if(itemp+itemp<>index) then jndex+=1;
            index:=itemp
        end;
        jj:=jndex;
        if(jj>=ii) then begin
            t:=x[jj];
            x[jj]:=x[ii];
            x[ii]:=t
        end
    end;
    if(k<0) then begin
        for ii:=0 to n-1 do x[ii]*=gain;
    end
end;

end.