HOW FAST FOURIER TRANSFORM WORKS

A basic building block in the fast Fourier transform is called ``doubling.'' Given a series $(x_0, x_1, \ldots , x_{N - 1})$and its sampled Fourier transform $(X_0, X_1, \ldots , X_{N - 1})$,and another series $(y_0, y_1, \ldots , y_{N - 1})$and its sampled Fourier transform $(Y_0, Y_1, \ldots , Y_{N - 1})$,there is a trick to find easily the transform of the interlaced double-length series

$$z_t \eq (x_0, y_0, x_1, y_1, \ldots , x_{N - 1}, y_{N-1})$$ (20)

The process of doubling is used many times during the computing of a fast Fourier transform. As the word ``doubling" might suggest, it will be convenient to suppose that N is an integer formed by raising 2 to some integer power. Suppose N = 8 = 2³. We begin by dividing our eight-point series into eight separate series, each of length one. The Fourier transform of each of the one-point series is just the point. Next, we use doubling four times to get the transforms of the four different two-point series (x₀, x₄), (x₁, x₅), (x₂, x₆), and (x₃, x₇). We use doubling twice more to get the transforms of the two different four-point series (x₀, x₂, x₄, x₆) and (x₁, x₃, x₅, x₇). Finally, we use doubling once more to get the transform of the original eight-point series $(x_0, x_1, x_2, \ldots , x_7)$.It remains to look into the details of the doubling process. Let

$$\begin{eqnarray} V & = & e^{i2\pi/2N} = W^{1/2} \ V^N & = & e^{i\pi} = -1\end{eqnarray}$$ (21) (22)

By definition, the transforms of two N-point series are

$$\begin{eqnarray} X_k &= & \sum^{N- 1}_{j = 0} x_j V^{2jk} \quad (k = 0, 1, \ldot... ...& \sum^{N- 1}_{j = 0} y_j V^{2jk} \quad (k = 0, 1, \ldots , N -1) \end{eqnarray}$$ (23) (24)

Likewise, the transform of the interlaced series $z_j = (x_0, y_0, x_1, y_1, \ldots , x_{N - 1}, y_{N- 1})$is

$$Z_k \eq \sum^{2N- 1}_{l = 0} z_l V^{lk} \quad (k = 0, 1, \ldots , 2N -1)$$ (25)

To make Z_k from X_k and Y_k, we require two separate formulas, one for k = 0, 1, $\ldots$, N -1, and the other for k = N, N + 1, $\ldots$, 2N - 1. Start from the sum

$$Z_k \eq \sum^{2N- 1}_{l = 0} z_l V^{lk} \quad (k = 0, 1, \ldots , N -1)$$ (26)

and then split the sum into two parts, noting that x_j multiplies even powers of V, and y_j multiplies odd powers:

$$\begin{eqnarray} Z_k &= & \sum^{N- 1}_{j = 0} x_j V^{2jk} + V^k \sum^{N - 1}_{j = 0} y_j V^{2jk} \ &= & X_k + V^k Y_k\end{eqnarray}$$ (27) (28)

We obtain the last half of the Z_k by

$$\begin{eqnarray} Z_k &= & \sum^{2N- 1}_{l = 0} z_l V^{lk} \quad (k = N, N + 1, \... ...= & X_{k-N} - V^{k-N} Y_{k-N} \quad (k = N, N + 1, \ldots , 2N -1)\end{eqnarray}$$ (29) (30) (31) (32) (33) (34) (35)

The subroutine ftu() does not follow this analysis in detail.

If you would like some industrial grade FFT programs, search the web for "prime factor FFT".