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:

(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
(2.19) |

and

(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:

^{8}F_{(1:4)} |
= | ||

^{8}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_fftThere 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 NONEstatement, 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 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 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, ")", /))') yThe 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.