An Algorithm for the Machine Calculation of Complex Fourier Series

Thuật toán cho việc tính toán chuỗi Fourier phức bằng máy

Bởi James W. Cooley và John W. Tukey

Một phương pháp hiệu quả để tính toán các tương tác của một thí nghiệm giai thừa $2^m$ đã được Yates đưa ra và được biết rộng rãi theo tên ông. Sự tổng quát hóa cho $3^m$ được Box và cộng sự đưa ra [1]. Good [2] đã tổng quát hóa các phương pháp này và đưa ra các thuật toán thanh lịch, trong đó một lớp ứng dụng là tính toán chuỗi Fourier. Ở dạng tổng quát đầy đủ, các phương pháp của Good áp dụng cho một số bài toán trong đó cần nhân một vectơ $N$ với một ma trận $N \times N$ có thể được phân tích thành $m$ ma trận thưa, trong đó $m$ tỉ lệ với $\log N$. Điều này dẫn đến một thủ tục yêu cầu số phép toán tỉ lệ với $N \log N$ thay vì $N^2$. Các phương pháp này được áp dụng ở đây cho việc tính toán chuỗi Fourier phức. Chúng hữu ích trong các tình huống khi số điểm dữ liệu là, hoặc có thể được chọn là, một số có tính hợp thành cao. Thuật toán được suy ra và trình bày ở đây dưới một dạng khá khác biệt. Sự chú ý được dành cho việc lựa chọn $N$. Cũng chỉ ra rằng có thể thu được lợi thế đặc biệt khi sử dụng máy tính nhị phân với $N = 2^m$ và toàn bộ phép tính có thể được thực hiện bên trong mảng gồm $N$ vị trí lưu trữ dữ liệu được sử dụng cho các hệ số Fourier đã cho.

Xét bài toán tính toán chuỗi Fourier phức

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.