Winograd's Short DFT Algorithms

In 1976, S. Winograd presented a new DFT algorithm which had significantly fewer multiplications than the Cooley-Tukey FFT which had been published eleven years earlier. This new Winograd Fourier Transform Algorithm (WFTA) is based on the type- one index map from Multidimensional Index Mapping with each of the relatively prime length short DFT's calculated by very efficient special algorithms. It is these short algorithms that this section will develop. They use the index permutation of Rader described in the another module to convert the prime length short DFT's into cyclic convolutions. Winograd developed a method for calculating digital convolution with the minimum number of multiplications. These optimal algorithms are based on the polynomial residue reduction techniques of Polynomial Description of Signals: Equation 1 to break the convolution into multiple small ones , , , , , . The operation of discrete convolution defined by y ( n ) = ∑ k h ( n - k ) x ( k ) is called a bilinear operation because, for a fixed h(n), y(n) is a linear function of x(n) and for a fixed x(n) it is a linear function of h(n). The operation of cyclic convolution is the same but with all indices evaluated modulo N. Recall from Polynomial Description of Signals: Equation 3 that length-N cyclic convolution of x(n) and h(n) can be represented by polynomial multiplication Y ( s ) = X ( s ) H ( s ) mod ( s N - 1 ) This bilinear operation of and can also be expressed in terms of linear matrix operators and a simpler bilinear operator denoted by o which may be only a simple element-by-element multiplication of the two vectors , , . This matrix formulation is Y = C [ A X o B H ] where X, H and Y are length-N vectors with elements of x(n), h(n) and y(n) respectively. The matrices A and B have dimension M x N , and C is N x M with M≥N. The elements of A, B, and C are constrained to be simple; typically small integers or rational numbers. It will be these matrix operators that do the equivalent of the residue reduction on the polynomials in . In order to derive a useful algorithm of the form to calculate , consider the polynomial formulation again. To use the residue reduction scheme, the modulus is factored into relatively prime factors. Fortunately the factoring of this particular polynomial, sN-1, has been extensively studied and it has considerable structure. When factored over the rationals, which means that the only coefficients allowed are rational numbers, the factors are called cyclotomic polynomials , , . The most interesting property for our purposes is that most of the coefficients of cyclotomic polynomials are zero and the others are plus or minus unity for degrees up to over one hundred. This means the residue reduction will generally require no multiplications. The operations of reducing X(s) and H(s) in are carried out by the matrices A and B in . The convolution of the residue polynomials is carried out by the o operator and the recombination by the CRT is done by the C matrix. More details are in , , , , but the important fact is the A and B matrices usually contain only zero and plus or minus unity entries and the C matrix only contains rational numbers. The only general multiplications are those represented by o. Indeed, in the theoretical results from computational complexity theory, these real or complex multiplications are usually the only ones counted. In practical algorithms, the rational multiplications represented by C could be a limiting factor. The h(n) terms are fixed for a digital filter, or they represent the W terms from Multidimensional Index Mapping: Equation 1 if the convolution is being used to calculate a DFT. Because of this, d=BH in can be precalculated and only the A and C operators represent the mathematics done at execution of the algorithm. In order to exploit this feature, it was shown , that the properties of allow the exchange of the more complicated operator C with the simpler operator B. Specifically this is given by Y = C [ A X o B H ] Y ' = B T [ A X o C T H ' ] where H' has the same elements as H, but in a permuted order, and likewise Y' and Y. This very important property allows precomputing the more complicated CTH' in rather than BH as in . Because BH or CTH' can be precomputed, the bilinear form of and can be written as a linear form. If an M x M diagonal matrix D is formed from d=CTH, or in the case of , d=BH, assuming a commutative property for o, becomes Y ' = B T D A X and becomes Y = C D A X In most cases there is no reason not to use the same reduction operations on X and H, therefore, B can be the same as A and then becomes Y ' = A T D A X In order to illustrate how the residue reduction is carried out and how the A matrix is obtained, the length-5 DFT algorithm started in The DFT as Convolution or Filtering: Matrix 1 will be continued. The DFT is first converted to a length-4 cyclic convolution by the index permutation from The DFT as Convolution or Filtering: Equation 3 to give the cyclic convolution in The DFT as Convolution or Filtering. To avoid confusion from the permuted order of the data x(n) in The DFT as Convolution or Filtering, the cyclic convolution will first be developed without the permutation, using the polynomial U(s) U ( s ) = x ( 1 ) + x ( 3 ) s + x ( 4 ) s 2 + x ( 2 ) s 3 U ( s ) = u ( 0 ) + u ( 1 ) s + u ( 2 ) s 2 + u ( 3 ) s 3 and then the results will be converted back to the permuted x(n). The length-4 cyclic convolution in terms of polynomials is Y ( s ) = U ( s ) H ( s ) mod ( s 4 - 1 ) and the modulus factors into three cyclotomic polynomials s 4 - 1 = ( s 2 - 1 ) ( s 2 + 1 ) = ( s - 1 ) ( s + 1 ) ( s 2 + 1 ) = P 1 P 2 P 3 Both U(s) and H(s) are reduced modulo these three polynomials. The reduction modulo P1 and P2 is done in two stages. First it is done modulo (s2-1), then that residue is further reduced modulo (s-1) and (s+1). U ( s ) = u 0 + u 1 s + u 2 s 2 + u 3 s 3 U ' ( s ) = ( ( U ( s ) ) ) ( s 2 - 1 ) = ( u 0 + u 2 ) + ( u 1 + u 3 ) s U 1 ( s ) = ( ( U ' ( s ) ) ) P 1 = ( u 0 + u 1 + u 2 + u 3 ) U 2 ( s ) = ( ( U ' ( s ) ) ) P 2 = ( u 0 - u 1 + u 2 - u 3 ) U 3 ( s ) = ( ( U ( s ) ) ) P 3 = ( u 0 - u 2 ) + ( u 1 - u 3 ) s The reduction in of the data polynomial can be denoted by a matrix operation on a vector which has the data as entries. 1 0 1 0 0 1 0 1 u 0 u 1 u 2 u 3 = u 0 + u 2 u 1 + u 3 and the reduction in is 1 0 - 1 0 0 1 0 - 1 u 0 u 1 u 2 u 3 = u 0 - u 2 u 1 - u 3 Combining and gives one operator 1 0 1 0 0 1 0 1 1 0 - 1 0 0 1 0 - 1 u 0 + u 2 u 1 + u 3 u 0 - u 2 u 1 - u 3 = u 0 + u 2 u 1 + u 3 u 0 - u 2 u 1 - u 3 = w 0 w 1 v 0 v 1 Further reduction of v0+v1s is not possible because P3=s2+1 cannot be factored over the rationals. However s2-1 can be factored into P1P2=(s-1)(s+1) and, therefore, w0+w1s can be further reduced as was done in and by 1 1 w 0 w 1 = w 0 + w 1 = u 0 + u 2 + u 1 + u 3 1 - 1 w 0 w 1 = w 0 - w 1 = u 0 + u 2 - u 1 - u 3 Combining , and gives 1 1 0 0 1 - 1 0 0 0 0 1 0 0 0 0 1 1 0 1 0 0 1 0 1 1 0 - 1 0 0 1 0 - 1 u 0 u 1 u 2 u 3 = r 0 r 1 v 0 v 1 The same reduction is done to H(s) and then the convolution of is done by multiplying each residue polynomial of X(s) and H(s) modulo each corresponding cyclotomic factor of P(s) and finally a recombination using the polynomial Chinese Remainder Theorem (CRT) as in Polynomial Description of Signals: Equation 9 and Polynomial Description of Signals: Equation 13. Y ( s ) = K 1 ( s ) U 1 ( s ) H 1 ( s ) + K 2 ( s ) U 2 ( s ) H 2 ( s ) + K 3 ( s ) U 3 ( s ) H 3 ( s ) mod (s4-1) where U1(s)=r1 and U2(s)=r2 are constants and U3(s)=v0+v1s is a first degree polynomial. U1 times H1 and U2 times H2 are easy, but multiplying U3 time H3 modulo (s2+1) is more difficult. The multiplication of U3(s) times H3(s) can be done by the Toom-Cook algorithm , , which can be viewed as Lagrange interpolation or polynomial multiplication modulo a special polynomial with three arbitrary coefficients. To simplify the arithmetic, the constants are chosen to be plus and minus one and zero. The details of this can be found in , , . For this example it can be verified that ( ( v 0 + v 1 s ) ( h 0 + h 1 s ) ) ) s 2 + 1 = ( v 0 h 0 - v 1 h 1 ) + ( v 0 h 1 + v 1 h 0 ) s which by the Toom-Cook algorithm or inspection is 1 - 1 0 - 1 - 1 1 1 0 0 1 1 1 v 0 v 1 o 1 0 0 1 1 1 h 0 h 1 = y 0 y 1 where o signifies point-by-point multiplication. The total A matrix in is a combination of and giving A X = A 1 A 2 A 3 X = 1 0 0 0 0 1 0 0 0 0 1 0 0 0 0 1 0 0 1 1 1 1 0 0 1 - 1 0 0 0 0 1 0 0 0 0 1 1 0 1 0 0 1 0 1 1 0 - 1 0 0 1 0 - 1 u 0 u 1 u 2 u 3 = r 0 r 1 v 0 v 1 where the matrix A3 gives the residue reduction s2-1 and s2+1, the upper left-hand part of A2 gives the reduction modulo s-1 and s+1, and the lower right-hand part of A1 carries out the Toom-Cook algorithm modulo s2+1 with the multiplication in . Notice that by calculating in the three stages, seven additions are required. Also notice that A1 is not square. It is this “expansion" that causes more than N multiplications to be required in o in or D in . This staged reduction will derive the A operator for The method described above is very straight-forward for the shorter DFT lengths. For N=3, both of the residue polynomials are constants and the multiplication given by o in is trivial. For N=5, which is the example used here, there is one first degree polynomial multiplication required but the Toom-Cook algorithm uses simple constants and, therefore, works well as indicated in . For N=7, there are two first degree residue polynomials which can each be multiplied by the same techniques used in the N=5 example. Unfortunately, for any longer lengths, the residue polynomials have an order of three or greater which causes the Toom-Cook algorithm to require constants of plus and minus two and worse. For that reason, the Toom-Cook method is not used, and other techniques such as index mapping are used that require more than the minimum number of multiplications, but do not require an excessive number of additions. The resulting algorithms still have the structure of . Blahut and Nussbaumer have a good collection of algorithms for polynomial multiplication that can be used with the techniques discussed here to construct a wide variety of DFT algorithms. The constants in the diagonal matrix D can be found from the CRT matrix C in using d=CTH' for the diagonal terms in D. As mentioned above, for the smaller prime lengths of 3, 5, and 7 this works well but for longer lengths the CRT becomes very complicated. An alternate method for finding D uses the fact that since the linear form or calculates the DFT, it is possible to calculate a known DFT of a given x(n) from the definition of the DFT in Multidimensional Index Mapping: Equation 1 and, given the A matrix in , solve for D by solving a set of simultaneous equations. The details of this procedure are described in . A modification of this approach also works for a length which is an odd prime raised to some power: N=PM. This is a bit more complicated , but has been done for lengths of 9 and 25. For longer lengths, the conventional Cooley-Tukey type- two index map algorithm seems to be more efficient. For powers of two, there is no primitive root, and therefore, no simple conversion of the DFT into convolution. It is possible to use two generators , , to make the conversion and there exists a set of length 4, 8, and 16 DFT algorithms of the form in in . In an operation count of several short DFT algorithms is presented. These are practical algorithms that can be used alone or in conjunction with the index mapping to give longer DFT's as shown in The Prime Factor and Winograd Fourier Transform Algorithms. Most are optimized in having either the theoretical minimum number of multiplications or the minimum number of multiplications without requiring a very large number of additions. Some allow other reasonable trade-offs between numbers of multiplications and additions. There are two lists of the number of multiplications. The first is the number of actual floating point multiplications that must be done for that length DFT. Some of these (one or two in most cases) will be by rational constants and the others will be by irrational constants. The second list is the total number of multiplications given in the diagonal matrix D in . At least one of these will be unity ( the one associated with X(0)) and in some cases several will be unity ( for N=2M ). The second list is important in programming the WFTA in The Prime Factor and Winograd Fourier Transform Algorithm: The Winograd Fourier Transform Algorithm. Length Mult Mult Adds N Non-one Total 2 0 4 4 3 4 6 12 4 0 8 16 5 10 12 34 7 16 18 72 8 4 16 52 9 20 22 84 11 40 42 168 13 40 42 188 16 20 36 148 17 70 72 314 19 76 78 372 25 132 134 420 32 68 - 388

Because of the structure of the short DFTs, the number of real multiplications required for the DFT of real data is exactly half that required for complex data. The number of real additions required is slightly less than half that required for complex data because (N-1) of the additions needed when N is prime add a real to an imaginary, and that is not actually performed. When N=2m, there are (N-2) of these pseudo additions. The special case for real data is discussed in , , . The structure of these algorithms are in the form of X'=CDAX or BTDAX or ATDAX from and . The A and B matrices are generally M by N with M≥N and have elements that are integers, generally 0 or ±1. A pictorial description is given in .

The flow graph in should be compared with the matrix description of and , and with the programs in , , , and the appendices. The shape in illustrates the expansion of the data by A. That is to say, AX has more entries than X because M>N. The A operator consists of additions, the D operator gives the M multiplications (some by one) and AT contracts the data back to N values with additions only. M is one half the second list of multiplies in . An important characteristic of the D operator in the calculation of the DFT is its entries are either purely real or imaginary. The reduction of the W vector by (s(N-1)/2-1) and (s(N-1)/2+1) separates the real and the imaginary constants. This is discussed in , . The number of multiplications for complex data is only twice those necessary for real data, not four times. Although this discussion has been on the calculation of the DFT, very similar results are true for the calculation of convolution and correlation, and these will be further developed in Algorithms for Data with Restrictions. The ATDA structure and the picture in are the same for convolution. Algorithms and operation counts can be found in , , .

The Bilinear Structure The bilinear form introduced in and the related linear form in are very powerful descriptions of both the DFT and convolution. Bilinear: Y = C [ A X o B H ] Linear: Y = C D A X Since is a bilinear operation defined in terms of a second bilinear operator o , this formulation can be nested. For example if o is itself defined in terms of a second bilinear operator @, by X o H = C ' [ A ' X @ B ' H ] then becomes Y = C C ' [ A ' A X @ B ' B H ] For convolution, if A represents the polynomial residue reduction modulo the cyclotomic polynomials, then A is square (e.g. and o represents multiplication of the residue polynomials modulo the cyclotomic polynomials. If A represents the reduction modulo the cyclotomic polynomials plus the Toom-Cook reduction as was the case in the example of , then A is NxM and o is term-by- term simple scalar multiplication. In this case AX can be thought of as a transform of X and C is the inverse transform. This is called a rectangular transform because A is rectangular. The transform requires only additions and convolution is done with M multiplications. The other extreme is when A represents reduction over the N complex roots of sN-1. In this case A is the DFT itself, as in the example of (43), and o is point by point complex multiplication and C is the inverse DFT. A trivial case is where A, B and C are identity operators and o is the cyclic convolution. This very general and flexible bilinear formulation coupled with the idea of nesting in gives a description of most forms of convolution.

Winograd's Complexity Theorems Because Winograd's work , , , , , has been the foundation of the modern results in efficient convolution and DFT algorithms, it is worthwhile to look at his theoretical conclusions on optimal algorithms. Most of his results are stated in terms of polynomial multiplication as Polynomial Description of Signals: Equation 3 or . The measure of computational complexity is usually the number of multiplications, and only certain multiplications are counted. This must be understood in order not to misinterpret the results. This section will simply give a statement of the pertinent results and will not attempt to derive or prove anything. A short interpretation of each theorem will be given to relate the result to the algorithms developed in this chapter. The indicated references should be consulted for background and detail. Theorem 1 Given two polynomials, x(s) and h(s), of degree N and M respectively, each with indeterminate coefficients that are elements of a field H, N+M+1 multiplications are necessary to compute the coefficients of the product polynomial x(s)h(s). Multiplication by elements of the field G (the field of constants), which is contained in H, are not counted and G contains at least N+M distinct elements. The upper bound in this theorem can be realized by choosing an arbitrary modulus polynomial P(s) of degree N+M+1 composed of N+M+1 distinct linear polynomial factors with coefficients in G which, since its degree is greater than the product x(s)h(s), has no effect on the product, and by reducing x(s) and h(s) to N+M+1 residues modulo the N+M+1 factors of P(s). These residues are multiplied by each other, requiring N+M+1 multiplications, and the results recombined using the Chinese remainder theorem (CRT). The operations required in the reduction and recombination are not counted, while the residue multiplications are. Since the modulus P(s) is arbitrary, its factors are chosen to be simple so as to make the reduction and CRT simple. Factors of zero, plus and minus unity, and infinity are the simplest. Plus and minus two and other factors complicate the actual calculations considerably, but the theorem does not take that into account. This algorithm is a form of the Toom-Cook algorithm and of Lagrange interpolation , , , . For our applications, H is the field of reals and G the field of rationals. Theorem 2 If an algorithm exists which computes x(s)h(s) in N+M+1 multiplications, all but one of its multiplication steps must necessarily be of the form m k = ( g k ' + x ( g k ) ) ( g k " + h ( g k ) ) for k = 0 , 1 , . . . , N + M where gk are distinct elements of G; and gk' and gk" are arbitrary elements of G This theorem states that the structure of an optimal algorithm is essentially unique although the factors of P(s) may be chosen arbitrarily. Theorem 3 Let P(s) be a polynomial of degree N and be of the form P(s)=Q(s)k, where Q(s) is an irreducible polynomial with coefficients in G and k is a positive integer. Let x(s) and h(s) be two polynomials of degree at least N-1 with coefficients from H, then 2N-1 multiplications are required to compute the product x(s)h(s) modulo P(s). This theorem is similar to Theorem 1 with the operations of the reduction of the product modulo P(s) not being counted. Theorem 4 Any algorithm that computes the product x(s)h(s) modulo P(s) according to the conditions stated in Theorem 3 and requires 2N-1 multiplications will necessarily be of one of three structures, each of which has the form of Theorem 2 internally. As in Theorem 2, this theorem states that only a limited number of possible structures exist for optimal algorithms. Theorem 5 If the modulus polynomial P(s) has degree N and is not irreducible, it can be written in a unique factored form P(s)=P1m1(s)P2m2(s)...Pkmk(s) where each of the Pi(s) are irreducible over the allowed coefficient field G. 2N-k multiplications are necessary to compute the product x(s)h(s) modulo P(s) where x(s) and h(s) have coefficients in H and are of degree at least N-1. All algorithms that calculate this product in 2N-k multiplications must be of a form where each of the k residue polynomials of x(s) and h(s) are separately multiplied modulo the factors of P(s) via the CRT. Corollary: If the modulus polynomial is P(s)=sN-1, then 2N-t(N) multiplications are necessary to compute x(s)h(s) modulo P(s), where t(N) is the number of positive divisors of N. Theorem 5 is very general since it allows a general modulus polynomial. The proof of the upper bound involves reducing x(s) and h(s) modulo the k factors of P(s). Each of the k irreducible residue polynomials is then multiplied using the method of Theorem 4 requiring 2Ni-1 multiplies and the products are combined using the CRT. The total number of multiplies from the k parts is 2N-k. The theorem also states the structure of these optimal algorithms is essentially unique. The special case of P(s)=sN-1 is interesting since it corresponds to cyclic convolution and, as stated in the corollary, k is easily determined. The factors of sN-1 are called cyclotomic polynomials and have interesting properties , , . Theorem 6 , Consider calculating the DFT of a prime length real-valued number sequence. If G is chosen as the field of rational numbers, the number of real multiplications necessary to calculate a length-P DFT is u(DFT(N))=2P-3-t(P-1) where t(P-1) is the number of divisors of P-1. This theorem not only gives a lower limit on any practical prime length DFT algorithm, it also gives practical algorithms for N=3,5, and 7. Consider the operation counts in to understand this theorem. In addition to the real multiplications counted by complexity theory, each optimal prime-length algorithm will have one multiplication by a rational constant. That constant corresponds to the residue modulo (s-1) which always exists for the modulus P(s)=sN-1-1. In a practical algorithm, this multiplication must be carried out, and that accounts for the difference in the prediction of Theorem 6 and count in . In addition, there is another operation that for certain applications must be counted as a multiplication. That is the calculation of the zero frequency term X(0) in the first row of the example in The DFT as Convolution or Filtering: Matrix 1. For applications to the WFTA discussed in The Prime Factor and Winograd Fourier Transform Algorithms: The Winograd Fourier Transform Algorithm, that operation must be counted as a multiply. For lengths longer than 7, optimal algorithms require too many additions, so compromise structures are used. Theorem 7 , If G is chosen as the field of rational numbers, the number of real multiplications necessary to calculate a length-N DFT where N is a prime number raised to an integer power: N=Pm , is given by u ( D F T ( N ) ) = 2 N - ( ( m 2 + m ) / 2 ) t ( P - 1 ) - m - 1 where t(P-1) is the number of divisors of (P-1). This result seems to be practically achievable only for N=9, or perhaps 25. In the case of N=9, there are two rational multiplies that must be carried out and are counted in but are not predicted by Theorem 7. Experience indicates that even for N=25, an algorithm based on a Cooley-Tukey FFT using a type 2 index map gives an over-all more balanced result. Theorem 8 If G is chosen as the field of rational numbers, the number of real multiplications necessary to calculate a length-N DFT where N=2m is given by u ( D F T ( N ) ) = 2 N - m 2 - m - 2 This result is not practically useful because the number of additions necessary to realize this minimum of multiplications becomes very large for lengths greater than 16. Nevertheless, it proves the minimum number of multiplications required of an optimal algorithm is a linear function of N rather than of NlogN which is that required of practical algorithms. The best practical power-of-two algorithm seems to the Split-Radix FFT discussed in The Cooley-Tukey Fast Fourier Transform Algorithm: The Split-Radix FFT Algorithm. All of these theorems use ideas based on residue reduction, multiplication of the residues, and then combination by the CRT. It is remarkable that this approach finds the minimum number of required multiplications by a constructive proof which generates an algorithm that achieves this minimum; and the structure of the optimal algorithm is, within certain variations, unique. For shorter lengths, the optimal algorithms give practical programs. For longer lengths the uncounted operations involved with the multiplication of the higher degree residue polynomials become very large and impractical. In those cases, efficient suboptimal algorithms can be generated by using the same residue reduction as for the optimal case, but by using methods other than the Toom-Cook algorithm of Theorem 1 to multiply the residue polynomials. Practical long DFT algorithms are produced by combining short prime length optimal DFT's with the Type 1 index map from Multidimensional Index Mapping to give the Prime Factor Algorithm (PFA) and the Winograd Fourier Transform Algorithm (WFTA) discussed in The Prime Factor and Winograd Fourier Transform Algorithms. It is interesting to note that the index mapping technique is useful inside the short DFT algorithms to replace the Toom-Cook algorithm and outside to combine the short DFT's to calculate long DFT's.

The Automatic Generation of Winograd's Short DFTs by Ivan Selesnick, Polytechnic Institute of New York University

Introduction Efficient prime length DFTs are important for two reasons. A particular application may require a prime length DFT and secondly, the maximum length and the variety of lengths of a PFA or WFTA algorithm depend upon the availability of prime length modules. This , , , discusses automation of the process Winograd used for constructing prime length FFTs , for N<7 and that Johnson and Burrus extended to N<19. It also describes a program that will design any prime length FFT in principle, and will also automatically generate the algorithm as a C program and draw the corresponding flow graph. Winograd's approach uses Rader's method to convert a prime length DFT into a P-1 length cyclic convolution, polynomial residue reduction to decompose the problem into smaller convolutions , , and the Toom-Cook algorithm , . The Chinese Remainder Theorem (CRT) for polynomials is then used to recombine the shorter convolutions. Unfortunately, the design procedure derived directly from Winograd's theory becomes cumbersome for longer length DFTs, and this has often prevented the design of DFT programs for lengths greater than 19. Here we use three methods to facilitate the construction of prime length FFT modules. First, the matrix exchange property , , is used so that the transpose of the reduction operator can be used rather than the more complicated CRT reconstruction operator. This is then combined with the numerical method for obtaining the multiplication coefficients rather than the direct use of the CRT. We also deviate from the Toom-Cook algorithm, because it requires too many additions for the lengths in which we are interested. Instead we use an iterated polynomial multiplication algorithm . We have incorporated these three ideas into a single structural procedure that automates the design of prime length FFTs.

Matrix Description It is important that each step in the Winograd FFT can be described using matrices. By expressing cyclic convolution as a bilinear form, a compact form of prime length DFTs can be obtained. If y is the cyclic convolution of h and x, then y can be expressed as y = C [ A x . * B h ] where, using the Matlab convention, .* represents point by point multiplication. When A,B, and C are allowed to be complex, A and B are seen to be the DFT operator and C, the inverse DFT. When only real numbers are allowed, A, B, and C will be rectangular. This form of convolution is presented with many examples in . Using the matrix exchange property explained in and this form can be written as y = R B T [ C T R h . * A x ] where R is the permutation matrix that reverses order. When h is fixed, as it is when considering prime length DFTs, the term CTRh can be precomputed and a diagonal matrix D formed by D=diag{CTRh}. This is advantageous because in general, C is more complicated than B, so the ability to “hide" C saves computation. Now y=RBTDAx or y=RATDAx since A and B can be the same; they implement a polynomial reduction. The form y=RTDAxT can also be used for the prime length DFTs, it is only necessary to permute the entries of x and to ensure that the DC term is computed correctly. The computation of the DC term is simple, for the residue of a polynomial modulo a-1 is always the sum of the coefficients. After adding the x0 term of the original input sequence, to the s-l residue, the DC term is obtained. Now DFT{x}=RATDAx. In Johnson observes that by permuting the elements on the diagonal of D, the output can be permuted, so that the R matrix can be hidden in D, and DFT{x}=ATDAx. From the knowledge of this form, once A is found, D can be found numerically .

Programming the Design Procedure Because each of the above steps can be described by matrices, the development of a prime length FFTs is made convenient with the use of a matrix oriented programming language such as Matlab. After specifying the appropriate matrices that describe the desired FFT algorithm, generating code involves compiling the matrices into the desired code for execution. Each matrix is a section of one stage of the flow graph that corresponds to the DFT program. The four stages are: Permutation Stage: Permutes input and output sequence. Reduction Stage: Reduces the cyclic convolution to smaller polynomial products. Polynomial Product Stage: Performs the polynomial multiplications. Multiplication Stage: Implements the point-by-point multiplication in the bilinear form. Each of the stages can be clearly seen in the flow graphs for the DFTs. shows the flow graph for a length 17 DFT algorithm that was automatically drawn by the program.

The programs that accomplish this process are written in Matlab and C. Those that compute the appropriate matrices are written in Matlab. These matrices are then stored as two ASCII flies, with the dimensions in one and the matrix elements in the second. A C program then reads the flies and compiles them to produce the final FFT program in C

The Reduction Stage The reduction of an Nth degree polynomial, X(s), modulo the cyclotomic polynomial factors of (sN-1) requires only additions for many N, however, the actual number of additions depends upon the way in which the reduction proceeds. The reduction is most efficiently performed in steps. For example, if N=4 and ((X(s))s-1,((X(s))s+1and ((X(s))s2+1 where the double parenthesis denote polynomial reduction modulo (s-1), s+1, and s2+1), then in the first step ((X(s)))s2-1, and ((Xs)))s2+1 should be computed. In the second step, ((Xs)))s-1 and ((Xs)))s+1 can be found by reducing ((X(s)))s2-1 This process is described by the diagram in .

When N is even, the appropriate first factorization is (SN/2-1)(sN/2+1), however, the next appropriate factorization is frequently less obvious. The following procedure has been found to generate a factorization in steps that coincides with the factorization that minimizes the cumulative number of additions incurred by the steps. The prime factors of N are the basis of this procedure and their importance is clear from the useful well-known equation sN-1=∏n|NCn(s) where Cn(s) is the nth cyclotomic polynomial. We first introduce the following two functions defined on the positive integers, ψ ( N ) = the smallest prime factor of N for N > 1 and ψ(1)=1. Suppose P(s) is equal to either (sN-1) or an intermediate noncyclotomic polynomial appearing in the factorization process, for example, (a2-1), above. Write P(s) in terms of its cyclotomic factors, P ( s ) = C k 1 ( s ) C k 2 ( s ) ⋯ C k L define the two sets, G and G, by G = { k 1 , ⋯ , k L } and G = { k / g c d ( G ) : k ∈ G } and define the two integers, t and T, by t = min { ψ ( k ) : k ∈ G , k > 1 } and T = max n u ( k , t ) : k ∈ G } Then form two new sets, A = { k ∈ G : T ∣ k } and B = { k ∈ G : T | k } The factorization of P(s), P ( s ) = ( ∏ k ∈ A C k ( s ) ) ( ∏ k ∈ B C k ( s ) ) has been found useful in the procedure for factoring (sN-1). This is best illustrated with an example. Example: N=36 Step 1. Let P(s)=s36-1. Since P=C1C2C3C4C6C9C12C18C36 G = G = { 1 , 2 , 3 , 4 , 6 , 9 , 12 , 18 , 36 } t = min { 2 , 3 } = 2 A = { k ∈ G : 4 | k } = { 1 , 2 , 3 , 6 , 9 , 18 } B = { k ∈ G : 4 | k } = { 4 , 12 , 36 } Hence the factorization of s36-1 into two intermediate polynomials is as expected, ∏ k ∈ A C k ( s ) = s 18 - 1 , ∏ k ∈ B C k ( s ) = s 18 + 1 If a 36th degree polynomial, X(s), is represented by a vector of coefficients, X=(x35,⋯,x0)', then ((X(s))s18-1 (represented by X') and ((X(s))s18+1 (represented by X") is given by t e s t which entails 36 additions. Step 2. This procedure is repeated with P(s)=s18-1 and P(s)=s18+1. We will just show it for the later. Let P(s)=s18+1. Since P=C4C12C36 G = { 4 , 12 , 36 } , G ' = { l , 3 , 9 } t = min 3 = 3 T = max ν ( k , 3 ) : k ∈ G = max l , 3 , 9 = 9 A = k ∈ G : 9 | k } = { 4 , 12 } B = k ∈ G : 9 | k } = { 36 } This yields the two intermediate polynomials, s 6 + 1 , and s 12 - s 6 + 1 In the notation used above, X ' X ' ' = I 6 - I 6 I 6 I 6 I 6 - I 6 I 6 X entailing 24 additions. Continuing this process results in a factorization in steps In order to see the number of additions this scheme uses for numbers of the form N=P-1 (which is relevant to prime length FFT algorithms) figure 4 shows the number of additions the reduction process uses when the polynomial X(s) is real. Figure 4: Number of Additions for Reduction Stage

The Polynomial Product Stage The iterated convolution algorithm can be used to construct an N point linear convolution algorithm from shorter linear convolution algorithms . Suppose the linear convolution y, of the n point vectors x and h (h known) is described by y = E n T D E n x where En is an “expansion" matrix the elements of which are ±l's and 0's and D is an appropriate diagonal matrix. Because the only multiplications in this expression are by the elements of D, the number of multiplications required, M(n), is equal to the number of rows of En. The number of additions is denoted by A(n). Given a matrix En1 and a matrix En2, the iterated algorithm gives a method for combining En1 and En2 to construct a valid expansion matrix, En, for N≤n1n2. Specifically, E n 1 , n 2 = ( I M ( n 2 ) ⊗ E n 1 ) ( E n 2 × I n 1 ) The product n1n2 may be greater than N, for zeros can be (conceptually) appended to x. The operation count associated with En1,n2 is A ( n 1 , n 2 ) = n ! A ( n 2 ) + A ( n 1 ) M n 2 ) M ( n 1 , n 2 ) = M ( n 1 ) M ( n 2 ) Although they are both valid expansion matrices, En1,n2≠En2,n1 and An1,n2≠An2,n1 Because Mn1,n2≠Mn2,n1 it is desirable to chose an ordering of factors to minimize the additions incurred by the expansion matrix. The following , follows from above.

Multiple Factors Note that a valid expansion matrix, EN, can be constructed from En1,n2 and En3, for N≤n1n2n3. In general, any number of factors can be used to create larger expansion matrices. The operation count associated with En1,n2,n3 is A ( n 1 , n 2 , n 3 ) = n 1 n 2 A ( n 3 ) + n 1 A ( n 2 ) M ( n 3 ) + A ( n 1 ) M ( n 2 ) M ( n 3 ) M ( n 1 , n 2 , n 3 ) = M ( n 1 ) M ( n 2 ) M ( n 3 ) These equations generalize in the predicted way when more factors are considered. Because the ordering of the factors is relevant in the equation for A(.) but not for M(.), it is again desirable to order the factors to minimize the number of additions. By exploiting the following property of the expressions for A(.) and M(.), the optimal ordering can be found . reservation of Optimal Ordering. Suppose A(n1,n2,n3)≤min{A(nk1,nk2,nk3):k1,k2,k3∈{1,2,3} and distinct}, then A ( n 1 , n 2 ) ≤ A ( n 2 , n 1 ) A ( n 2 , n 3 ) ≤ A ( n 3 , n 2 ) A ( n 1 , n 3 ) ≤ A ( n 3 , n 1 ) The generalization of this property to more than two factors reveals that an optimal ordering of {n1,⋯,nL-i} is preserved in an optimal ordering of {n1,⋯nL}. Therefore, if (n1,⋯nL) is an optimal ordering of {n1,⋯nL}, then (nk,nk+1) is an optimal ordering of {nk,nk+1} and consequently A ( n k ) M ( n k ) - n k ≤ A ( n k + 1 ) M ( n k + 1 ) - n k + 1 for all k=1,2,⋯,L-1. This immediately suggests that an optimal ordering of {n1,⋯nL} is one for which A ( n 1 ) M ( n 1 ) - n 1 ⋯ A ( n L ) M ( n L ) - n L is nondecreasing. Hence, ordering the factors, {n1,⋯nL}, to minimize the number of additions incurred by En1,⋯,nL simply involves computing the appropriate ratios.

Discussion and Conclusion We have designed prime length FFTs up to length 53 that are as good as the previous designs that only went up to 19. Table 1 gives the operation counts for the new and previously designed modules, assuming complex inputs. It is interesting to note that the operation counts depend on the factorability of P-1. The primes 11, 23, and 47 are all of the form 1+2P1 making the design of efficient FFTs for these lengths more difficult. Further deviations from the original Winograd approach than we have made could prove useful for longer lengths. We investigated, for example, the use of twiddle factors at appropriate points in the decomposition stage; these can sometimes be used to divide the cyclic convolution into smaller convolutions. Their use means, however, that the 'center* multiplications would no longer be by purely real or imaginary numbers. N Mult Adds 7 16 72 11 40 168 13 40 188 17 82 274 19 88 360 23 174 672 29 190 766 31 160 984 37 220 920 41 282 1140 43 304 1416 47 640 2088 53 556 2038

The approach in writing a program that writes another program is a valuable one for several reasons. Programming the design process for the design of prime length FFTs has the advantages of being practical, error-free, and flexible. The flexibility is important because it allows for modification and experimentation with different algorithmic ideas. Above all, it has allowed longer DFTs to be reliably designed.More details on the generation of programs for prime length FFTs can be found in the 1993 Technical Report.