  ## 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

```