Diskrete Fourier-Transformation

Da in PASCAL die komplexe Arithmetik nicht unterstützt wird, muß die e-Funktion mit der Eulerschen Formel in Real- und Imaginärteil zerlegt werden:

program DiskreteFourierTransformation;

        {Berechnung einer diskreten Fouriertransformation am Beispiel
        y(x) = cos(3t) + sin(10t)}

    uses
        stdfile;

    const
        MaxArray = 1024; {<= 32767}

    var
        k, m, n:      integer;
        omega, angle: double;
        funkt:        array [0..MaxArray] of double; {Vorgegebene Funktion}
        realt:        array [0..MaxArray] of double; {Realteil}
        imagt:        array [0..MaxArray] of double; {Imaginaerteil}

    procedure init;
        var
            fName: string;
            f:     text;
        begin {init}
            writeln('Darstellung der Berechnung einer diskreten Fouriertransformation');
            writeln;
            write('In welcher Datei liegt die Funktion?');
            writeln('In dft.dat liegt als Beispiel: y(x) = cos(3t) + sin(10t) in 32 Intervallen von delta t = 2Pi/32');
            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 StdRead(f, funkt[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, ' Werte eingelesen.');
        end; {init}

    procedure result; 
	var
	    i: Integer;
        begin {result}
	    for i:=0 to n-1 do begin
		writeln(i:2,'   ', realt[i]:7:3,' +i* ', imagt[i]:6:3)
	    end {for}
	end; {result}

    begin {DFT}
        init;
        omega := 2.0*Pi/n;
        for k:=0 to n-1 do begin
            realt[k] := 0.0;
            imagt[k] := 0.0;
            for m:=0 to n-1 do begin
                angle := k*omega*m;
                realt[k] := realt[k] + funkt[m]*cos(angle);
                imagt[k] := imagt[k] - funkt[m]*sin(angle);
            end {for}
        end; {for}
        result
    end. {DFT}