The fast Fourier transform (FFT) is an algorithm which can take the discrete Fourier transform
of a array of size *n* = 2^{N} in **Θ**(*n* ln(*n*)) time.
This algorithm is generally performed in place and this implementation continues in that
tradition. Two implementations are provided:

- The first implementation,
`FFT.basic.h`, emphasizes the algorithm; consequently,
to simplify the presentation, it allocates additional memory when dividing the vector
into even and odd entries and explicit recursion is used, and
- The second implementation,
`FFT.array.h`, demonstrates how the FFT algorithm
can be performed with only **O**(1) additional memory and without recursive function
calls. Reference, Maple 8, Waterloo Maple Inc., Waterloo, Ontario.

The website Fastest Fourier Transform in the West (see
wikipedia for a summary) provides the fastest portable
free software implementation of the FFT.

Some comments on the code:

The value π is found using `double const PI = 4.0*std::atan( 1.0 );` The inverse tangent function
of 1 returns π/4 and and multiplication by a power of 2 does not affect the relative error of a floating
point number: in this case, multiplication by four simply adds two to the exponent.

The testing function applies the fast Fourier transform to the vector
(2, 2, 0, 0, 0, 0, 0, 0)^{T}. It prints the
vector, the fast Fourier transform of the vector, and the inverse fast Fourier transform
of that result. The output
is

(1,0) (1,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0)
(2,0) (1.70711,-0.707107) (1,-1) (0.292893,-0.707107) (0,0) (0.292893,0.707107) (1,1) (1.70711,0.707107)
(1,0) (1,0) (0,0) (0,0) (0,0) (0,0) (0,0) (0,0)

You will note that 1/√2 ≈ 0.707107, and 1 − 1/√2 ≈ 0.292893 .

The function printing the complex array will print 0 in place of either component
being less than 10^{-15}. This is to make aid in the visualization of the
result.

# The Discrete Fourier Transform

The discrete Fourier transform maps a vector from the space domain to the complex frequency domain. To describe
this, recall that a vector of *n* dimensions represents coordinates in space: the basis vectors
are the *n* unit vectors (1, 0, 0, ···), (0, 1, 0, ···), *etc*. These are shown in Figure 1.

Figure 1. The unit basis vectors for a vector of dimension 8.

While these unit vectors are exceelent for representing the specific location of
a vector, they do not suggest any periodicity. The discrete Fourier transforms is a linear transformation
from a basis of unit vectors to a basis of vectors which are capture the periodic behaviour of the
unit vector: if the copies of the vector were laid end-to-end, what periodicity would there be
in the patterns and what are the periods thereof.

Figure 2. The basis vectors of the frequency space of dimension 8.

`template <typename T> class complex<T>`

This code uses the STL complex class which has a constructor which takes the real and imaginary parts
as arguments. The following functions are friends of the class and for the argument *z* = *a* + *jb*:

Function | Returns |

`double abs(z)` | |*z*| = √(*a*^{2} + *b*^{2}) |

`double arg(z)` | atan2(ℑ(*z*), ℜ(*z*)) = atan2(*b*, *a*) |

`complex conj(z)` | *z*^{*} = *a* − *jb* |

`double imag(z)` | ℑ(*z*) = *b* |

`double norm(z)` | |*z*|^{2} = *a*^{2} + *b*^{2} |

`complex polar(r, θ)` | *r**e*^{jθ} |

`double real(z)` | ℜ(*z*) = *a* |

`complex exp(z)` | *e*^{z} |

`complex log(z)` | ln(*z*) |

`complex pow(w, z)` | *w*^{z} |

`complex sqrt(z)` | √*z* |

`complex sin(z)` | sin(*z*) |

`complex cos(z)` | cos(*z*) |

`complex sinh(z)` | sinh(*z*) |

`complex cosh(z)` | cosh(*z*) |

A complex number *a* + *jb* is printed as the ordered pair `(a,b)`.

# Determining Powers of Two

The assertion `assert( n > 0 && (n & (~n + 1)) == n );` checks to
ensure that the argument *n* is both positive and a power of two. The operation
`n & (~n + 1)` selects the least significant one in the integer. If the least
significant one is also the most significant one, then the bitwise and will return
that value; otherwise, the result will not equal the original value *n*. This is shown
in Figure 3 with both a power of two (64) and a number which is not a power of two (84).

Figure 3. Selecting the least
significant one (1) bit in a binary number.