Introduction

Introduction

The development of algorithms for the fast computation of the Discrete FourierTransform in the last 30 years originated with the radix 2 Cooley-Tukey FFTand the theory and variety of FFTs has grown significantly since then. Most of the work has focused onFFTs whose sizes are composite, for the algorithms depend on the ability to factorthe length of the data sequence, so that the transform can be found by takingthe transform of smaller lengths. For this reason, algorithms for prime length transforms arebuilding blocks for many composite length FFTs -the maximum length and the variety of lengths of a PFA or WFTA algorithm depend upon theavailability of prime length FFT modules. As such, prime length Fast Fourier Transformsare a special, important and difficult case.

Fast algorithms designed for specific short prime lengths have been developed and have been written as straight line code [link] , [link] . These dedicated programs rely upon an observation made in Rader's paper [link] in which he shows that a prime $p$ length DFT can be found by performing a $p-1$ length circular convolution.Since the publication of that paper, Winograd had developed a theory of multiplicative complexity for transforms and designed algorithms for convolution that attainthe minimum number of multiplications [link] . Although Winograd's algorithms are very efficient for small prime lengths,for longer lengths they require a large number of additions and the algorithms become very cumbersome to design.This has prevented the design of useful prime length FFT programs for lengths greater than 31.Although we have previously reported the design of programs for prime lengths greater than 31 [link] those programs required more additions than necessary and were long.Like the previously existing ones, they simply consisted of a long list of instructionsand did not take advantage of the attainable common structures.

In this paper we describe a set of programs for circular convolution and prime length FFTsthat are are short, possess great structure, share many computational procedures,and cover a large variety of lengths. Because the underlying convolution is decomposed into aset of disjoint operations they can be performed in paralleland this parallelism is made clear in the programs. Moreover, each of these independent operations is made up of a sequenceof sub-operations of the form $I\otimes A\otimes I$ where $\otimes$ denotes the Kronecker product. These operations can be implemented asvector/parallel operations [link] . Previous programs for prime length FFTs do nothave these features: they consist of straight line code and are not amenable tovector/parallel implementations.

We have also developed a program that automatically generates these programs for circular convolution and prime length DFTs.This code generating program requires information only about a set of modules for computing cyclotomic convolutions.We compute these non-circular convolutions by computing a linear convolution and reducing the result.Furthermore, because these linear convolution algorithms can be built from smaller ones, the only modulesneeded are ones for the linear convolution of prime length sequences. It turns out that with linear convolution algorithms foronly the lengths 2 and 3, we can generate a wide variety of prime length FFT algorithms.In addition, the code we generate is made up of calls to a relatively small set of functions.Accordingly, the subroutines can be designed and optimized to specifically suit a given architecture.

The programs we describe use Rader's conversion of a prime point DFT into a circular convolution, but this convolutionwe compute using the split nesting algorithm [link] . As Stasinski notes [link] , this yields algorithms possessing greater structure and simpler programs anddoesn't generally require more computation.

On the row-column method

In computing the DFT of an $n=n_1{n}_2$ point sequence where $n_1$ and $n_2$ are relatively prime, a row-column method can be employed.That is, if an $n_1\times n_2$ array is appropriately formed from the $n$ point sequence, then its DFT can be computed by computing the DFTof the rows and by then computing the DFT of the columns. The separability of the DFT makes this possible.It should be mentioned, however, that in at least two papers [link] , [link] it is mistakenly assumed that the row-column method can also be applied to convolution.Unfortunately, the convolution of two sequences can not be found by forming two arrays, byconvolving their rows, and by then convolving their columns. This misunderstanding about the separability of convolutionalso appears in [link] where the author incorrectly writes a diagonal matrix of a bilinear form as a Kronecker product.If it were a Kronecker product, then there would indeed exist a row-column method for convolution.

Earlier reports on this work were published in the conference proceedings [link] , [link] , [link] and a fairly complete report was published in the IEEE Transaction on Signal Processing [link] . Some parts of this approach appear in the Connexions book, Fast Fourier Transforms . This work is built on and an extension of that in [link] which is also in the Connexions Technical Report .