50 lines
1.1 KiB
Plaintext
50 lines
1.1 KiB
Plaintext
test: PROCEDURE OPTIONS (MAIN, REORDER); /* Derived from Fortran Rosetta Code */
|
|
|
|
/* In-place Cooley-Tukey FFT */
|
|
FFT: PROCEDURE (x) RECURSIVE;
|
|
DECLARE x(*) COMPLEX FLOAT (18);
|
|
DECLARE t COMPLEX FLOAT (18);
|
|
DECLARE ( N, Half_N ) FIXED BINARY (31);
|
|
DECLARE ( i, j ) FIXED BINARY (31);
|
|
DECLARE (even(*), odd(*)) CONTROLLED COMPLEX FLOAT (18);
|
|
DECLARE pi FLOAT (18) STATIC INITIAL ( 3.14159265358979323E0);
|
|
|
|
N = HBOUND(x);
|
|
|
|
if N <= 1 THEN return;
|
|
|
|
allocate odd((N+1)/2), even(N/2);
|
|
|
|
/* divide */
|
|
do j = 1 to N by 2; odd((j+1)/2) = x(j); end;
|
|
do j = 2 to N by 2; even(j/2) = x(j); end;
|
|
|
|
/* conquer */
|
|
call fft(odd);
|
|
call fft(even);
|
|
|
|
/* combine */
|
|
half_N = N/2;
|
|
do i=1 TO half_N;
|
|
t = exp(COMPLEX(0, -2*pi*(i-1)/N))*even(i);
|
|
x(i) = odd(i) + t;
|
|
x(i+half_N) = odd(i) - t;
|
|
end;
|
|
|
|
FREE odd, even;
|
|
|
|
END fft;
|
|
|
|
|
|
DECLARE data(8) COMPLEX FLOAT (18) STATIC INITIAL (
|
|
1, 1, 1, 1, 0, 0, 0, 0);
|
|
DECLARE ( i ) FIXED BINARY (31);
|
|
|
|
call fft(data);
|
|
|
|
do i=1 TO 8;
|
|
PUT SKIP LIST ( fixed(data(i), 25, 12) );
|
|
end;
|
|
|
|
END test;
|