Task

Calculate the   FFT   (Fast Fourier Transform)   of an input sequence.

The most general case allows for complex numbers at the input
and results in a sequence of equal length, again of complex numbers.
If you need to restrict yourself to real numbers, the output should
be the magnitude   (i.e.:   sqrt(re2 + im2))   of the complex result.

The classic version is the recursive Cooley–Tukey FFT. Wikipedia has pseudo-code for that.
Further optimizations are possible but not required.




 
PROGRAM RDFT;
 
(*)
 
        Free Pascal Compiler version 3.2.0 [2020/06/14] for x86_64
        The free and readable alternative at C/C++ speeds
        compiles natively to almost any platform, including raspberry PI *
        Can run independently from DELPHI / Lazarus
 
        For debian Linux: apt -y install fpc
        It contains a text IDE called fp
 
 
(*)
 
USES
 
    crt,
    math,
    sysutils,
    ucomplex;
 
 
 
TYPE
 
    table = array  of complex;
 
 
 
PROCEDURE Split ( T: table ; EVENS: table; ODDS:table ) ;
 
    VAR
 
        k:  integer ; 
 
    BEGIN
 
        FOR k := 0 to Length ( T ) - 1 DO
 
            IF Odd ( k ) THEN
 
                ODDS [ k DIV 2 ]    := T [ k ]
 
            ELSE
 
                EVENS [ k DIV 2 ]   := T [ k ]
 
    END;
 
 
 
PROCEDURE WriteCTable ( L: table ) ;
 
    VAR
 
        x   :integer ;
 
    BEGIN
 
        FOR x := 0  to length ( L ) - 1 DO
 
            BEGIN
 
                Write   ( Format ('%3.3g ' , [ L [ x ].re ] ) ) ;
 
                IF ( L [ x ].im >= 0.0 ) THEN Write ( '+' ) ;
 
                WriteLn ( Format ('%3.5gi' , [ L [ x ].im ] ) ) ;
 
            END ;
 
    END;
 
 
 
FUNCTION FFT ( L : table ): table ;
 
    VAR
 
        k       :   integer ;
        N       :   integer ;
        halfN   :   integer ;
        E       :   table   ;
        Even    :   table   ;
        O       :   table   ;
        Odds    :   table   ;
        T       :   complex ;
 
    BEGIN
 
        N   :=  length ( L )        ;
 
        IF N < 2 THEN
 
            EXIT ( L )  ;
 
        halfN := ( N DIV 2 )        ;
 
        SetLength ( E,    halfN )   ;
 
        SetLength ( O,    halfN )   ;                
 
        Split     ( L, E, O )       ;
 
        SetLength ( L, 0 )   	    ;
 
        SetLength ( Even, halfN )   ;
 
        Even :=     FFT ( E )       ;
 
        SetLength ( E   , 0 )       ;
 
        SetLength ( Odds, halfN )   ;
 
        Odds :=     FFT ( O )       ;
 
        SetLength ( O   , 0 )       ;
 
        SetLength ( L,    N )       ;
 
        FOR k := 0 to halfN - 1 DO
 
            BEGIN
 
                T               := Cexp ( -2 * i * pi * k / N ) * Odds [ k ];
 
                L [ k ]         := Even [ k ] + T                      	    ;
 
                L [ k + halfN ] := Even [ k ] - T                      	    ;
 
            END ;
 
        SetLength ( Even, 0 )   ;
 
        SetLength ( Odds, 0 )   ;
 
        FFT :=  L               ;
 
    END ;
 
 
 
VAR
 
    Ar  :   array of complex ;
 
    x   :   integer ;
 
BEGIN
 
 
 
    SetLength ( Ar, 8 ) ;
 
    FOR x := 0 TO 3 DO
 
	BEGIN
 
	    Ar [ x ]        :=  1.0 ;
 
	    Ar [ x + 4 ]    :=  0.0 ;
 
	END;
 
    WriteCTable ( FFT ( Ar ) )  ;
 
    SetLength ( Ar, 0 ) ;
 
 
 
END.
(*)
    Output:
 
        4 +  0i
        1 -2.4142i
        0 +  0i
        1 -0.41421i
        0 +  0i
        1 +0.41421i
        0 +  0i
        1 +2.4142i
 
 
 
