Sande-Tukey-FFT-Algorithmus

program SchnelleFourierTransformation;

        {Berechnung einer diskreten Fouriertransformation mittels des
        Sande-Tukey-FFT-Algorithmus}

    uses
        stdfile;

    const
        MaxArray = 1024; {<= 32767}


    var
        i, j, k, l, m, n, n1, n2:          integer;
        angle, argument, xdum, ydum, c, s: double;
        x:                                 array [0..MaxArray] of double;
        y:                                 array [0..MaxArray] of double;

    procedure init;
        var
            fName: string;
            f:     text;
        begin {init}
            writeln('Darstellung der Berechnung einer diskreten Fouriertransformation');
            writeln('mittels des Sande-Tukey-FFT-Algorithmus.)');
            writeln;
            writeln('In welcher Datei liegt die Funktion? (In fft.dat liegt als Beispiel y(x)=cos(3t)-sin(10t.)');
            readln(fName);
            if StdOpen(f, fName) then begin
                n := -1;
                repeat
                    n := n + 1;
                    if n > MaxArray then begin
                        writeln('Die Datei "', fName, '" hat zu viele Punkte (mehr als ', MaxArray, ').');
                        halt
                    end; {if}
                until not StdRead2(f, x[n], y[n]);
                n := n - 1
            end else begin {if}
                writeln('Die Datei "', fName, '" kann nicht geoeffnet werden.');
                halt
            end; {if}
            close(f);
            writeln('Insgesamt wurden ', n + 1, ' Werte eingelesen.');
            n := n + 1;
            m := trunc(ln(n)/ln(2.0));
        end; {init}

    procedure result; 
	var 
	    i: integer;
	begin {result}
	    for i:=1 to n-1 do begin
		writeln(y[i]/n:7:3,' ',x[i]/n:7:3)
	    end {for}
	end; {result}

    begin {FFT}
        init;
        n2 := n;
        for k:=1 to m do begin
            n1 := n2;
            n2 := n2 div 2;
            angle := 0.0;
            argument := 2.0*Pi/n1;
            for j:=0 to n2-1 do begin
                c := cos(angle);
                s := sin(angle);
                i := j;
                while i < n-1 do begin
                    l    := i + n2;
                    xdum := x[i] - x[l];
                    x[i] := x[i] + x[l];
                    ydum := y[i] - y[l];
                    y[i] := y[i] + y[l];
                    x[l] := xdum*c - ydum*s;
                    y[l] := ydum*c + xdum*s;
                    i := i+n1
                end; {while}
                angle := (j+1)*argument
            end {for}
        end; {for}
        j := 0;
        for i:=0 to n-2 do begin
            if i < j then begin
                xdum := x[j];
                x[j] := x[i];
                x[i] := xdum;
                ydum := y[j];
                y[j] := y[i];
                y[i] := ydum
            end; {if}
            k := n div 2;
            while k < j+1 do begin
               j := j - k;
               k := k div 2
            end; {while}
            j := j + k
        end; {for}
        result
    end. {FFT}