A Recursive Implementation of the Fast Fourier Transform

Before we can embark on implementing the algorithm discussed above in Octave one word of caution. Throughout this course I have been referring to Octave as ``interactive Fortran'', but, although in many respects Octave's notation is quite similar to that of Fortran, sometimes it is annoyingly different. The situation is a little like with differences between Bourne Shell and C-shell - those little differences can be quite enough to derail one's shell scripts and one's composure too.

And the difference that I need to point to you right now relates to the way ranges and strides are described in both environments. And so, in Fortran (1:10:2) means ``take integers from 1 through 10 with a stride of 2'', i.e., 1, 3, 5, 7, 9. In Octave the same is accomplished by (1:2:10):

octave:5> (1:2:10)
ans =

   1   3   5   7   9

octave:6>

Now let us have a look at the top FFT formula:

$${}^8\boldsymbol{F}(\boldsymbol{f}_{[1:8]}) = \left( \begin{... ...4\boldsymbol{F}(\boldsymbol{f}_{[2:8:2]}) \end{array} \right)$$ (2.18)

How would we implement this in Octave? First let us assume that our input data lives on an 8-entries long vector x, and that we have a function called my_fft that can evaluate a discrete Fourier Transform from any regularly spaced subset of x, for example x(1:2:8) - remember that here we already use Octave notation, which means that the second, not the third, integer in the sequence (1:2:8) defines the stride. So

$${}^4\boldsymbol{F}(\boldsymbol{f}_{[1:8:2]}) = \hbox{\tt my\_fft}\left(x(1:2:8)\right)$$ (2.19)

and

$${}^4\boldsymbol{F}(\boldsymbol{f}_{[2:8:2]}) = \hbox{\tt my\_fft}\left(x(2:2:8)\right)$$ (2.20)

What the function is going to deliver in this case is an array of 4 complex numbers that represents a 4-point discrete Fourier transform. To get the 8-point transform:

$$\hbox{\tt my\_fft}\left(x(1:2:8)\right) + [1, \omega_8, \omega_8^2, \omega_8^3] \cdot \hbox{\tt my\_fft}\left(x(2:2:8)\right)$$ ⁸F_(1:4) =

$$\hbox{\tt my\_fft}\left(x(1:2:8)\right) - [1, \omega_8, \omega_8^2, \omega_8^3] \cdot \hbox{\tt my\_fft}\left(x(2:2:8)\right)$$ ⁸F_(5:8) =

And the beauty of recursive algorithms that this is really the only thing that we have to code. Plus the condition that stops the recursion, of course.

So, here is the Octave code:

octave:1> function y = my_fft(x)
> n = length(x);
> if n == 1
>    y = x;
> else
>    m = n/2;
>    y_top = my_fft(x(1:2:n));
>    y_bottom = my_fft(x(2:2:n));
>    d = exp(-2 * pi * i / n) .^ (0:m-1);
>    z = d .* y_bottom;
>    y = [ y_top + z , y_top - z ];
> end
> endfunction
octave:2> x = 2 * pi / 8 * (0:7);
octave:3> x = sin(x);
octave:4> 2 * pi / 8 * my_fft(x)
ans =

 Columns 1 through 3:

   0.00000 + 0.00000i  -0.00000 - 3.14159i   0.00000 - 0.00000i

 Columns 4 through 6:

   0.00000 + 0.00000i   0.00000 + 0.00000i   0.00000 + 0.00000i

 Columns 7 and 8:

   0.00000 + 0.00000i  -0.00000 + 3.14159i

octave:35>

It is easy enough to translate this code to Fortran-90, remembering, of course, that (1:2:n) needs to be replaced with (1:n:2). Here is the translation not only of function my_fft itself, which in the Fortran program shown below is implemented as a subroutine fft, but also of computations that lead to setting up vector x:

PROGRAM try_fft

  IMPLICIT NONE

  DOUBLE PRECISION, PARAMETER :: pi = 3.141592653589793d0
  INTEGER :: i
  COMPLEX(kind=8), DIMENSION(8) :: x, y

  x = (/ (i, i=0,7) /)
  x = 2 * pi / 8 * x
  x = SIN(x)
  CALL fft(x, y); y = 2 * pi / 8 * y
  WRITE(*,'(8 ("(", en12.3, ",", en12.3, ")", /))') y
  STOP 0

CONTAINS

  RECURSIVE SUBROUTINE fft(x, y)

    IMPLICIT NONE

    COMPLEX(kind=8), DIMENSION(:) :: x
    COMPLEX(kind=8), DIMENSION(SIZE(x)) :: y

    COMPLEX(kind=8), DIMENSION(:), ALLOCATABLE :: y_top, y_bottom, z, d
    DOUBLE PRECISION, PARAMETER :: pi = 3.141592653589793d0

    INTEGER :: n, m, i

    n = SIZE(x)
    IF (n .EQ. 1) THEN
       y = x
    ELSE
       m = n/2
       ALLOCATE(y_top(m), y_bottom(m), d(m), z(m))
       CALL fft(x(1:n:2), y_top)
       CALL fft(x(2:n:2), y_bottom)
       d = (/ (EXP(-2. * pi * (0, 1.) / n) ** i, i = 0, m-1) /)
       z = d * y_bottom
       y(1:m) = y_top + z; y(m+1:n) = y_top - z
    END IF
  END SUBROUTINE fft

END PROGRAM try_fft

There isn't really much more writing here than in the Octave program, although we have to be more formal, and, in particular, all variables must be declared, especially following the

IMPLICIT NONE

statement, and we have to remember to allocate space for y_top, y_bottom, d and z.

Here is how this program compiles and runs:

gustav@blanc:../src 18:01:22 !569 $ f90 -o fft fft.f90
gustav@blanc:../src 18:04:08 !570 $ ./fft
(   8.987E-18,   0.000E+00)
(-174.393E-18,  -3.142E+00)
(  96.184E-18, -87.197E-18)
( 418.010E-18,   0.000E+00)
( 183.380E-18,   0.000E+00)
( 174.393E-18,   0.000E+00)
(  96.184E-18,  87.197E-18)
(-802.744E-18,   3.142E+00)

 STOP 0
gustav@blanc:../src 18:04:11 !571 $

The algorithm presented in this section is perhaps the simplest of all FFT algorithms. It suffers from two inefficiencies. The first one is the statement:

d = (/ (EXP(-2. * pi * (0, 1.) / n) ** i, i = 0, m-1) /)

which results in calculating more terms of $\omega_n^k$ than is really necessary: some will be unnecessarily recalculated on climbing up from the recursion. This can be addressed relatively trivially: for example we could precompute only as many $\omega_n^k$ as are needed and then we would simply look them up.

The other inefficiency is recursion itself, which although it can be made nearly as efficient as iterations in some cases (e.g., for the so called ``tail recursion''), it will usually take more resources and run slower than a simple DO loop.

Regarding Fortran itself there are two new elements in the code, which you haven't encountered yet. The first one is the so called engineering notation, which is used when formatting the complex numbers:

  WRITE(*,'(8 ("(", en12.3, ",", en12.3, ")", /))') y

The format is flagged by en. The other two formats that can be used for floating point numbers are e and f.

The other novelty is the keyword RECURSIVE that appears in front of SUBROUTINE fft. A subroutine or a function must be declared RECURSIVE if it is such. Fortran-77 and earlier Fortrans did not support recursion, so this is still a new and seldom used feature to most Fortran programmers.