Algol 68 Genie

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

Deze website is gearchiveerd door de KB, nationale bibliotheek.

This website is archived by the National Library of the Netherlands.