An Algorithm for the Machine Calculation of Complex ...

1 An Algorithm for the Machine Calculation ofComplex Fourier SeriesBy James W. Cooley and John W. TukeyAn efficient method for the Calculation of the interactions of a 2m factorial ex-periment was introduced by Yates and is widely known by his name. The generaliza-tion to 3m was given by Box et al. [1]. Good [2] generalized these methods and gaveelegant algorithms for which one class of applications is the Calculation of Fourierseries. In their full generality, Good's methods are applicable to certain problems inwhich one must multiply an JV-vector by an JV X N matrix which can be factoredinto m sparse matrices, where m is proportional to log JV.

2 This results in a procedurerequiring a number of operations proportional to JV log JV rather than JV2. Thesemethods are applied here to the Calculation of Complex Fourier series. They areuseful in situations where the number of data points is, or can be chosen to be, ahighly composite number. The Algorithm is here derived and presented in a ratherdifferent form. Attention is given to the choice of JV. It is also shown how specialadvantage can be obtained in the use of a binary computer with JV = 2m and howthe entire Calculation can be performed within the array of JV data storage locationsused for the given Fourier the problem of calculating the Complex Fourier series(1) XiJ)=ZAik)-Wjk, 0,1.

3 ,iV-l,*-0where the given Fourier coefficients Aik) are Complex and W is the principalJVth root of unity,(2) W = straightforward Calculation using ( 1 ) would require JV2 operations where "opera-tion" means, as it will throughout this note, a Complex multiplication followed by acomplex Algorithm described here iterates on the array of given Complex Fourieramplitudes and yields the result in less than 2JV Iog2 JV operations without requiringmore data storage than is required for the given array A. To derive the Algorithm ,suppose JV is a composite, , JV = rvr2.

4 Then let the indices in (1) be expressed/ = in + o , o - 0, 1, , rx 1, i = 0, 1, , r2 - 1,( * Jk = hr2 + fcc, h = 0, 1, - , r2 1, h = 0, 1, , n , one can write(4) Xijijo) = Aih , ko)-W*">Wik\*0 *1 Received August 17, 1964. Research in part at Princeton University under the sponsorshipof the Army Research Office (Durham). The authors wish to thank Richard Garwin for hisessential role in communication and encouragement. 297 License or copyright restrictions may apply to redistribution; see W.

5 COOLEY AND JOHN W. TUKEYS incethe inner sum, over fci, depends only on 0 and k0 and can be defined as a new array,(6) Aiijo, h) = Aiki, ko)-Wlok r>.kiThe result can then be written(7) Xiji, jo) = Aiijo ,ko)- W(hri+Mk\There are JV elements in the array Ai, each requiring rx operations, giving a totalof JVri operations to obtain Ai. Similarly, it takes Nr2 operations to calculate Xfrom Ai. Therefore, this two-step Algorithm , given by (6) and (7), requires a totalof(8) T = Nin + r2) is easy to see how successive applications of the above procedure, starting withits application to (6), give an m-step Algorithm requiring(9) T = Niri + r2+ + rm)operations, where(10) N = rvr2--- r = Sjtj with Sj ,tj > 1, then s, + t < r3- unless y = t = 2, when s + t,- = r.

6 In general, then, using as many factors as possible provides a minimum to (9), butfactors of 2 can be combined in pairs without loss . If we are able to choose JV to behighly composite, we may make very real gains. If all r are equal to r, then, from( 10 ) we have(11) m = logrJVand the total number of operations is(12) Tir) = JV = rmsnl" , then we find thatT(13)so thatm-r + n-s + p-t + ,log2 JV = m-log2 r + n-log2 s + p-log21 +N log2 JVis a weighted mean of the quantitiesr s tlog2 r ' log2 s ' log21 'License or copyright restrictions may apply to redistribution.

7 See Calculation OF Complex FOURIER SERIES 299whose values run as followsr log r2 use of r = 3 is formally most efficient, but the gain is only about 6% overthe use of 2 or 4, which have other advantages. If necessary, the use of r, up to 10can increase the number of computations by no more than 50%. Accordingly, wecan find "highly composite" values of JV within a few percent of any given possible, the use of JV = rm with r = 2 or 4 offers important advantagesfor computers with binary arithmetic, both in addressing and in Algorithm with r = 2 is derived by expressing the indices in the form, N - -r2~* + + i-2+ o,(14)k = fcm_1-2m~1 + + fci-2 + fco,where / and kv are equal to 0 or 1 and are the contents of the respective bit positionsin the binary representation of / and k.

8 All arrays will now be written as functionsof the bits of their indices. With this convention (1 ) is written(15) Xijm-i, ,/ ) - Aikm-i, , ko)-wik-1-2m'l+-+ik",kQ kx km_\where the sums are over kv = 0, 1. Since(I ) T^J*m-l-2m_1 _ Jp-J0*m-l-2m-1the innermost sum of (15), over km^i, depends only on 0, fcm_2, , k0 and canbe written(17) Aiijo, km_2, , ko) = YiAikm-i, , ko)-Whkm-'-2m'\*m-lProceeding to the next innermost sum, over fcm_2, and so on, and usingng\ -pp'-*m-|-2' -< _ p^OJ_i-2'-l+.

9 +jo)*m-i-2' -'one obtains successive arrays,Aiijo, , jl-l , fcm-!-l , , ko)(19) = ( 0, * ,3,-*,h -i, ..,ko)-wl ->-2'-1+-+M-k-'-2m-1km-lfor I = 1, 2, , or copyright restrictions may apply to redistribution; see JAMES W. COOLEY AND JOHN W. TUKEYW riting out the sum this appears as l( 0 , , I-l , fcm-i-l , * , h)= Ai-iijo , , ji-2, 0, -i-i, ,ko)+ (-l)i,-'i"-,4w(i,, ,jt-t, 1, fcm_t_i, , *o) ,-+ -+*>- -,> jw . o, to the indexing convention, this is stored in a location whose index is(21) /0-2"-1 + + j,-i-2m~l + km+i 2m-'-1 + + can be seen in (20) that only the two storage locations with indices having 0 and1 in the 2m~l bit position are involved in the computation.

10 Parallel computation ispermitted since the operation described by (20 ) can be carried out with all values of o, , ji-2, and ko, , fcm-j-i simultaneously. In some applications* it is con-venient to use (20) to express A in terms of A _2, giving what is equivalent to analgorithm with r = last array calculated gives the desired Fourier sums,(22) -X"( m-i, , o) = Amijo, , /m_i)in such an order that the index of an X must have its binary bits put in reverseorder to yield its index in the array Am .In some applications, where Fourier sums are to be evaluated twice, the aboveprocedure could be programmed so that no bit-inversion is necessary.