Bruun's FFT algorithm: Difference between revisions

Revision as of 19:20, 12 March 2006

Revision as of 21:06, 9 August 2006

Line 1:

'''Bruun's algorithm''' is a [[fast Fourier transform]] (FFT) algorithm based on an unusual recursive [[polynomial]]-factorization approach, proposed for powers of two by G. Bruun in 1978 and generalized to arbitrary even composite sizes by H. Murakami in 1996. Because its operations involve only real coefficients until the last computation stage, it was initially proposed as a way to efficiently compute the [[discrete Fourier transform]] (DFT) of real data. Bruun's algorithm has not seen widespread use, however, as approaches based on the ordinary [[Cooley-Tukey FFT algorithm|Cooley-Tukey FFT algorithm]] have been successfully adapted to real data with at least as much efficiency. Furthermore, there is evidence that Bruun's algorithm may be intrinsically less accurate than Cooley-Tukey in the face of finite numerical precision (Storn, 1993).

Nevertheless, Bruun's algorithm illustrates an alternative algorithmic framework that can express both itself and the Cooley-Tukey algorithm, and thus provides an interesting perspective on FFTs that permits mixtures of the two algorithms and other generalizations.

Line 11:

k = 0,\dots,N-1. </math>

For convenience, let us denote the ''N'' [[root of unity|roots of unity]] by &omega;ω_''N''^''n'' (''n''=0..''N''-1):

:<math>\omega_N^n = e^{-\frac{2\pi i}{N} n }</math>

Line 29:

In order to compute the DFT, we need to evaluate the remainder of ''x''(''z'') modulo ''N'' [[monomial|monomials]] as described above. Evaluating these remainders one by one is equivalent to the evaluating the usual DFT formula directly, and requires O(''N''²) operations. However, one can ''combine'' these remainders recursively to reduce the cost, using the following trick: if we want to evaluate ''x''(''z'') modulo two polynomials ''U''(''z'') and ''V''(''z''), we can first take the remainder modulo their product ''U''(''z'') ''V''(''z''), which reduces the [[degree (mathematics)|degree]] of the polynomial ''x''(''z'') and makes subsequent modulo operations less computationally expensive.

The product of all of the monomials (''z'' - &omega;ω_''N''^''k'') for ''k''=0..''N''-1 is simply ''z''^''N''-1 (whose roots are clearly the ''N'' roots of unity). One then wishes to find a recursive factorization of ''z''^''N''-1 into polynomials of few terms and smaller and smaller degree. To compute the DFT, one takes ''x''(''z'') modulo each level of this factorization in turn, recursively, until one arrives at the monomials and the final result. If each level of the factorization splits every polynomial into an O(1) (constant-bounded) number of smaller polynomials, each with an O(1) number of nonzero coefficients, then the modulo operations for that level take O(''N'') time; since there will be a logarithmic number of levels, the overall complexity is O (''N'' log ''N'').

More explicitly, suppose for example that ''z''^''N''-1 = ''F''₁(''z'') ''F''₂(''z'') ''F''₃(''z''), and that ''F''_''k''(''z'') = ''F''_''k'',1(''z'') ''F''_''k'',2(''z''), and so on. The corresponding FFT algorithm would consist of first computing ''x''_''k''(''z'') = ''x''(''z'') mod

Line 39:

====Cooley-Tukey as polynomial factorization====

The standard decimation-in-frequency (DIF) radix-''r'' Cooley-Tukey algorithm corresponds closely to a recursive factorization. For example, radix-2 DIF Cooley-Tukey factors ''z''^''N''-1 into ''F''₁ = (''z''^''N''/2-1) and ''F''₂ = (''z''^''N''/2+1). These modulo operations reduce the degree of ''x''(''z'') by 2, which corresponds to dividing the problem size by 2. Instead of recursively factorizing ''F''₂ directly, though, Cooley-Tukey instead first computes ''x''₂(''z'' &omega;ω_''N''), shifting all the roots (by a ''twiddle factor'') so that it can apply the recursive factorization of ''F''₁ to both subproblems. That is, Cooley-Tukey ensures that all subproblems are also DFTs, whereas this is not generally true for an arbitrary recursive factorization (such as Bruun's, below).

== The Bruun factorization ==

Line 48:

:<math>z^{4M} + az^{2M} + 1 = (z^{2M} + \sqrt{2-a}z^M+1) (z^{2M} - \sqrt{2-a}z^M + 1)</math>

where ''a'' is ana real constant with |''a''| &le;≤ 2. At the end of the recursion, for ''M''=1, you are left with degree-2 polynomials that can then be evaluated modulo two roots (''z'' - &omega;ω_''N''^''k'') for each polynomial. Thus, at each recursive stage, all of the polynomials are factorized into two parts of half the degree, each of which has at most three nonzero terms, leading to an O (''N'' log ''N'') algorithm for the FFT.

Moreover, since all of these polynomials have purely real coefficients (until the very last stage), they automatically exploit the special case where the inputs ''x''_''n'' are purely real to save roughly a factor of two in computation and storage. One can also take straightforward advantage of the case of real-symmetric data for computing the [[discrete cosine transform]] (Chen and Sorensen, 1992).

Line 54:

==== Generalization to arbitrary radices ====

The Bruun factorization, and thus the Bruun FFT algorithm, was generalized to handle arbitrary ''even'' composite lengths, i.e. dividing the polynomial degree by an arbitrary ''radix'' (factor), as follows. First, we define a set of polynomials &phi;φ_{''N'',&alpha;α}(''z'') for positive integers ''N'' and for &alpha;α in <nowiki>[0,1)</nowiki> by:

:<math>\phi_{N, \alpha}(z) =