parallel-fft.a68
1 COMMENT
2
3 @section Synopsis
4
5 Parallel Fast Fourier Transform in recursive form.
6
7 Unnormalised Fast Fourier Transform in recursive form.
8 F (k) = F even (k) + exp (2 pi i k / n) * F odd (k).
9 Parameter dir = +- 1 determines direction of the transform.
10 Assume that the lower bound of f is zero, and that its length is a power of 2.
11
12 COMMENT
13
14 PROC fft = (REF () LONG COMPLEX f, INT dir) VOID:
15
16 IF INT length = UPB f + 1;
17 length > 1
18 THEN INT middle = length % 2;
19 # Calculate transforms at sublevels #
20 (0 .. middle - 1) LONG COMPLEX f even, f odd;
21 FOR i FROM 0 TO middle - 1
22 DO f even (i) := f (2 * i);
23 f odd (i) := f (2 * i + 1)
24 OD;
25 PAR (fft (f even, dir), fft (f odd, dir));
26 # Calculate transform at this level #
27 FOR k FROM 0 TO middle - 1
28 DO LONG REAL phi = dir * 2 * long pi * k / length;
29 LONG COMPLEX factor = long cos (phi) I long sin (phi) * f odd (k);
30 f (k) := f even (k) + factor;
31 f (k + middle) := f even (k) - factor
32 OD
33 FI;
34
35 SKIP
