Fast Fourier transform: Difference between revisions

=={{header|Pascal}}==

(*)

</lang pascal>

(*)

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 ;

R : table ;

T : table ;

BEGIN

N := length ( L ) ;

IF N < 2 THEN

EXIT ( L ) ;

halfN := ( N DIV 2 ) ;

SetLength ( R, N ) ;

SetLength ( T, halfN ) ;

SetLength ( E, halfN ) ;

SetLength ( O, halfN ) ;

SetLength ( Even, halfN ) ;

SetLength ( Odds, halfN ) ;

Split ( L, E, O ) ;

Even := FFT ( E ) ;

Odds := FFT ( O ) ;

FOR k := 0 to halfN - 1 DO

BEGIN

T [ k ] := Cexp ( -2 * i * pi * k / N ) * Odds [ k ];

R [ k ] := Even [ k ] + T [ k ] ;

R [ k + halfN ] := Even [ k ] - T [ k ] ;

END ;

SetLength ( T , 0 ) ;

SetLength ( E , 0 ) ;

SetLength ( O , 0 ) ;

SetLength ( Even, 0 ) ;

SetLength ( Odds, 0 ) ;

FFT := R ;

SetLength ( R, 0 ) ;

END ;

VAR

Ar : array [ 0..7 ] of complex ;

x : integer ;

BEGIN

FOR x := 0 TO 3 DO

Ar [ x ] := 1.0 ;

Ar [ x + 4 ] := 0.0 ;

WriteCTable ( FFT ( Ar ) ) ;

END.

(*)

Output:

4 + 0i

1 -2.4142i

0 + 0i

1 -0.41421i

0 + 0i

1 +0.41421i

0 + 0i

1 +2.4142i

JPD 2021/12/24

(*)