An Algorithm for the Machine Calculation of Complex Fourier Series

By James W. Cooley and John W. Tukey

An efficient method for the calculation of the interactions of a $2^m$ factorial experiment was introduced by Yates and is widely known by his name. The generalization to $3^m$ was given by Box et al. [1]. Good [2] generalized these methods and gave elegant algorithms for which one class of applications is the calculation of Fourier series. In their full generality, Good’s methods are applicable to certain problems in which one must multiply an $N$-vector by an $N \times N$ matrix which can be factored into $m$ sparse matrices, where $m$ is proportional to $\log N$. This results in a procedure requiring a number of operations proportional to $N \log N$ rather than $N^2$. These methods are applied here to the calculation of complex Fourier series. They are useful in situations where the number of data points is, or can be chosen to be, a highly composite number. The algorithm is here derived and presented in a rather different form. Attention is given to the choice of $N$. It is also shown how special advantage can be obtained in the use of a binary computer with $N = 2^m$ and how the entire calculation can be performed within the array of $N$ data storage locations used for the given Fourier coefficients.

Consider the problem of calculating the complex Fourier series

where the given Fourier coefficients $A(k)$ are complex and $W$ is the principal $N$th root of unity,

A straightforward calculation using (1) would require $N^2$ operations where “operation” means, as it will throughout this note, a complex multiplication followed by a complex addition.

The algorithm described here iterates on the array of given complex Fourier amplitudes and yields the result in less than $2N \log_2 N$ operations without requiring more data storage than is required for the given array $A$. To derive the algorithm, suppose $N$ is a composite, i.e., $N = r_1 \cdot r_2$. Then let the indices in (1) be expressed

Then, one can write

Received August 17, 1964. Research in part at Princeton University under the sponsorship of the Army Research Office (Durham). The authors wish to thank Richard Garwin for his essential role in communication and encouragement.

Since

the inner sum, over $k_1$, depends only on $j_0$ and $k_0$ and can be defined as a new array,

The result can then be written

There are $N$ elements in the array $A_1$, each requiring $r_1$ operations, giving a total of $Nr_1$ operations to obtain $A_1$. Similarly, it takes $Nr_2$ operations to calculate $X$ from $A_1$. Therefore, this two-step algorithm, given by (6) and (7), requires a total of

operations.

It is easy to see how successive applications of the above procedure, starting with its application to (6), give an $m$-step algorithm requiring

operations, where

If $r_j = s_j t_j$ with $s_j, t_j > 1$, then $s_j + t_j < r_j$ unless $s_j = t_j = 2$, when $s_j + t_j = r_j$.

In general, then, using as many factors as possible provides a minimum to (9), but factors of 2 can be combined in pairs without loss. If we are able to choose $N$ to be highly composite, we may make very real gains. If all $r_j$ are equal to $r$, then, from (10) we have

and the total number of operations is

If $N = r^m s^n t^p \cdots$, then we find that

so that

is a weighted mean of the quantities

whose values run as follows

The use of $r_j = 3$ is formally most efficient, but the gain is only about 6% over the use of 2 or 4, which have other advantages. If necessary, the use of $r_j$ up to 10 can increase the number of computations by no more than 50%. Accordingly, we can find "highly composite" values of $N$ within a few percent of any given large number.

Whenever possible, the use of $N = r^m$ with $r = 2$ or 4 offers important advantages for computers with binary arithmetic, both in addressing and in multiplication economy.

The algorithm with $r = 2$ is derived by expressing the indices in the form

where $j_v$ and $k_v$ are equal to 0 or 1 and are the contents of the respective bit positions in the binary representation of $j$ and $k$. All arrays will now be written as functions of the bits of their indices. With this convention (1) is written

where the sums are over $k_v = 0, 1$. Since

the innermost sum of (15), over $k_{m-1}$, depends only on $j_0, k_{m-2}, \cdots, k_0$ and can be written

Proceeding to the next innermost sum, over $k_{m-2}$, and so on, and using

one obtains successive arrays,

for $l = 1, 2, \cdots, m$.

Writing out the sum this appears as

According to the indexing convention, this is stored in a location whose index is

It can be seen in (20) that only the two storage locations with indices having 0 and 1 in the $2^{m-l}$ bit position are involved in the computation. Parallel computation is permitted since the operation described by (20) can be carried out with all values of $j_0, \cdots, j_{l-2}$, and $k_0, \cdots, k_{m-l-1}$ simultaneously. In some applications* it is convenient to use (20) to express $A_l$ in terms of $A_{l-2}$, giving what is equivalent to an algorithm with $r = 4$.

The last array calculated gives the desired Fourier sums,

in such an order that the index of an $X$ must have its binary bits put in reverse order to yield its index in the array $A_m$.

In some applications, where Fourier sums are to be evaluated twice, the above procedure could be programmed so that no bit-inversion is necessary. For example, consider the solution of the difference equation,

The present method could be first applied to calculate the Fourier amplitudes of $F(j)$ from the formula

The Fourier amplitudes of the solution are, then,

The $B(k)$ and $A(k)$ arrays are in bit-inverted order, but with an obvious modification of (20), $A(k)$ can be used to yield the solution with correct indexing.

A computer program for the IBM 7094 has been written which calculates three-dimensional Fourier sums by the above method. The computing time taken for computing three-dimensional $2^a \times 2^b \times 2^c$ arrays of data points was as follows:

• A multiple-processing circuit using this algorithm was designed by R. E. Miller and S. Winograd of the IBM Watson Research Center. In this case $r = 4$ was found to be most practical.

IBM Watson Research Center Yorktown Heights, New York

Bell Telephone Laboratories, Murray Hill, New Jersey

Princeton University Princeton, New Jersey

1. G. E. P. Box, L. R. Connor, W. R. Cousins, O. L. Davies (Ed.), F. R. Hirnsworth & G. P. Silitto, The Design and Analysis of Industrial Experiments, Oliver & Boyd, Edinburgh, 1954.
2. I. J. Good, "The interaction algorithm and practical Fourier series," J. Roy. Statist. Soc. Ser. B., v. 20, 1958, p. 361-372; Addendum, v. 22, 1960, p. 372-375. MR 21 1674; MR 23 A4231.

一种用于机器计算复傅里叶级数的算法

作者：James W. Cooley和John W. Tukey

Yates提出了一种用于计算$2^m$阶乘实验交互作用的高效方法，并以其名字广为人知。Box等人[1]给出了向$3^m$的推广。Good[2]进一步推广了这些方法，并给出了优雅的算法，其中一类应用是计算傅里叶级数。在最一般的形式下，Good的方法适用于这样一类问题：必须将一个$N$-向量与一个$N \times N$矩阵相乘，而该矩阵可以分解为$m$个稀疏矩阵，其中$m$与$\log N$成正比。由此得到的过程所需的运算次数与$N \log N$成正比，而不是与$N^2$成正比。本文将这些方法应用于复傅里叶级数的计算。当数据点数是或可以选择为一个高度合成数时，这些方法很有用。本文以一种颇为不同的形式推导并给出了该算法。文中还讨论了$N$的选择。此外，本文说明了在$N = 2^m$的二进制计算机上如何获得特殊的优势，以及如何在用于给定傅里叶系数的$N$个数据存储位置组成的数组内完成整个计算。

考虑计算复傅里叶级数的问题