An Algorithm for the Machine Calculation of Complex Fourier Series

複素フーリエ級数の機械計算のためのアルゴリズム

James W. Cooley and John W. Tukey による

$2^m$要因実験の交互作用を計算する効率的な方法はYatesによって導入され、広く彼の名で知られている。$3^m$への一般化はBox et al. [1]によって与えられた。Good [2]はこれらの方法を一般化し、その応用の一分野としてフーリエ級数の計算に用いることのできる優れたアルゴリズムを与えた。Goodの方法はその最も一般的な形では、$N$-ベクトルに、$m$個の疎行列に因数分解できる$N \times N$行列を乗算しなければならないある種の問題に適用できる。ただし、$m$は$\log N$に比例する。この結果、$N^2$ではなく$N \log N$に比例する数の演算を必要とする手続きが得られる。ここではこれらの方法を複素フーリエ級数の計算に適用する。これはデータ点の数が非常に合成数である場合、またはそのように選択できる場合に有用である。ここではアルゴリズムをかなり異なる形で導出し、提示する。$N$の選択についても考察する。また、$N = 2^m$とした二進コンピューターを使用することで特別な利点が得られること、および与えられたフーリエ係数に用いる$N$個のデータ記憶場所からなる配列内だけで計算全体を実行できることを示す。

複素フーリエ級数を計算する問題を考える。

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.