THE FAST FOURIER TRANSFORM

When we write the expression  

$$B(Z) \eq b_0 + b_1Z + b_2Z^2 + \cdots$$ (23)

we have both a time function and its Fourier transform. If we plot the coefficients $(b_0, b_1, \ldots )$, we plot the time function. If we evaluate and plot (23) at numerous real $\omega$,then we have plotted the transform. (Note that for real $\omega$, Z is of unit magnitude; i.e., on the unit circle.) Since $\omega$ is a continuous variable and everything in a computer is finite, how do we select a finite number of values $\omega_k$for plotting? The usual choice is to take evenly spaced frequencies. The lowest frequency can be zero. [Note $Z(\omega = 0) = e^{io} = 1$.] A frequency as high as $\omega = 2\pi$[Note $Z(\omega = 2\pi) = e^{i2\pi} = 1$ also] need not be considered, since (23) gives the same value for it as for zero frequency. Choosing uniformly spaced frequencies between these limits we have  

$$\omega_k \eq { (0, 1, 2, \ldots , M- 1) \, 2\pi \over M}$$ (24)

where M is some integer. Now let us abbreviate $B(Z(\omega_k))$ as B_k.

For the special case of an N-point time function where N = 4, (23) may be expressed by the matrix multiplication  

$$\left[ \begin{array} {c} B_0 \\ B_1 \\ B_2 \\ B_3 \end{a... ...gin{array} {c} b_0 \\ b_1 \\ b_2 \\ b_3 \end{array} \right]$$ (25)

where  

$$W \eq e^{2\pi i/N}$$ (26)

It is not essential to choose N = M as we have done in (25), but it is a convenience. There is no loss of generality because one may always append zeros to a time function before inserting it into (25). A convenience of the choice N = M is that the matrix in (25) will then be square and there will be an exact inverse. In fact, the inverse to (25) may be easily shown to be  

$$\left[ \begin{array} {c} b_0 \\ b_1 \\ b_2 \\ b_3 \end{a... ...gin{array} {c} B_0 \\ B_1 \\ B_2 \\ B_3 \end{array} \right]$$ (27)

Since 1/W is the complex conjugate of W, the matrices of (25) and (27) are just complex conjugates of one another. In fact, one observes no fundamental mathematical difference between time functions and frequency functions. This ``duality'' would be even more complete if we had used a scale factor of N<sup>-1/2</sup> in each of (25) and (27) rather than 1 in (25) and N<sup>-1</sup> in (27). Note also that time functions and frequency functions could be interchanged in the mnemonic table describing symmetries. In fact, our earlier observation that the product of two frequency functions amounts to a convolution of two time functions corresponds to the convolution of the corresponding two frequency functions. We will not ``provide'' this duality as it is standard fare in both mathematics and systems theory books. However we will occasionally call upon the reader to realize that in any theorem the meanings of ``time'' and ``frequency'' may be interchanged.

In making a plot of the transform B_k for $(k = 0, 1, \ldots , M - 1)$,the frequency axis ranges as $0 \leq \omega_k < 2\pi$.It is often more natural to display the interval $-\pi \leq \omega < \pi$.Since the transform is periodic with period $2\pi$,values of B_k on the interval $\pi \leq \omega <2\pi$may simply be moved to the interval $-\pi \leq \omega < 0$ for display.

Thus, for N = 8 one might plot successively

$$B_4 \ \ B_5 \ \ B_6 \ \ B_7 \ \ B_0 \ \ B_1 \ \ B_2 \ \ B_3$$ (28)

corresponding to values of $\omega$ equal to

$$-\pi, -{3\pi \over 4}, -{\pi \over 2}, -{\pi \over 4}, 0, {\pi \over 4}, {\pi \over 2}, {3\pi \over 4}$$ (29)

One advantage of this display interval is that for continuous time series which are sampled sufficiently densely in time the transform values B_k get small on both ends. If the time series is real, the real part of B_k has even symmetry about B₀; the imaginary part has odd symmetry about B₀. Then, one need not bother to display half the values. Choice of an odd value of N would enable us to put $\omega = 0$ exactly in the middle of the interval, but the reader will soon see why we stick to an even number of data points.

The matrix times vector operation in (25) requires N² multiplications and additions. The rest of this section describes a trick method, called the fast Fourier transform, of accomplishing the matrix multiplication in N log₂ N multiplications and additions. Since, for example, log₂ 1024 is 10, this is a tremendous saving in effort.

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})$,one finds 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})$$ (30)

The process of doubling is used many times during the process of computing 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 of one point each. 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}$$ (31) (32)

The transforms of two N-point series are by definition

$$\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}$$ (33) (34)

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

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

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

First

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

We 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}$$ (36) (37)

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, \... ...(V^N)}^l \\ &= & \sum^{2N- 1}_{l = 0} z_l V^{lm} (-1)^l \nonumber\end{eqnarray}$$ (38) (39) (40)

 

1-8
Figure 8 A program to do fast Fourier transform. Modified from Brenner. Calling this program twice returns the original data. SIGNI should be +1. on one call and -1. on the other. LX must be a power of 2.