Bruun's FFT algorithm

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 size 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 of real data. Bruun's algorithm has not seen widespread use, however, as approaches based on the ordinary Cooley-Tukey FFT algorithm (q.v.) 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 alternative perspective on FFTs that permits mixtures of the two algorithms and other generalizations.

A Polynomial Approach to the DFT

Recall that the DFT is defined by the formula:

${\displaystyle f_{j}=\sum _{k=0}^{n-1}x_{k}e^{-{\frac {2\pi i}{n}}jk}\qquad j=0,\dots ,n-1.}$ 

For convenience, let us denote the n roots of unity by ω_n^k:

${\displaystyle \omega _{n}^{k}=e^{-{\frac {2\pi i}{n}}k}}$ 

and define the the polynomial X(z) whose coefficients are x_k:

${\displaystyle X(z)=\sum _{k=0}^{n-1}x_{k}z^{k}}$ 

The DFT can then be understood as a reduction of this polynomial: f_j is given by:

${\displaystyle f_{j}=X(\omega _{n}^{j})=X(z)\mod (z-\omega _{n}^{j})}$ 

where mod (modulo) denotes the polynomial remainder operation. The key to fast algorithms like Bruun's or Cooley-Tukey comes from the fact that one can perform this set of n remainder operations in recursive stages.

Recursive Factorizations and FFTs

In order to compute the DFT, we need to evaluate the remainder of X(z) modulo n 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: X(z) mod U(z) = [X(z) mod U(z)*V(z)] mod U(z). That is, 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 of the polynomial X(z) and makes subsequent modulo operations less computationally expensive.

The product of all of the monomials (z - ω_n^j) for j=0..n-1 is simply zⁿ-1 (whose roots are clearly the n roots of unity). One then wishes to find a recursive factorization of zⁿ-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 each polynomials into a bounded number of smaller polynomials, each of bounded degree, 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).

Moreover, as long as the factorized polynomials are relatively prime (which for polynomials means that they have no common roots), one can construct a dual algorithm by reversing the process with the Chinese Remainder Theorem.

Cooley-Tukey as polynomial factorization

The standard decimation-in-frequency radix-r Cooley-Tukey algorithm corresponds to the recursive factorization:

${\displaystyle z^{rm}-\omega _{p}^{q}=\prod _{\ell =0}^{r-1}\left(z^{m}-\omega _{pr}^{q\ell }\right)}$ 

where w is a root of unity.

Thus, radix-2 Cooley-Tukey factors zⁿ-1 into (z^n/2-1) (z^n/2+1), and these into (z^n/4-1) (z^n/4+1) and (z^n/4-i) (z^n/4+i), respectively, and so on. In this way, each stage of modulo operations reduces the degree of X(z) by two, which corresponds to dividing the problem size by two.

The Bruun Factorization

The basic Bruun factorization for powers of two factorizes z^2m-1 into (z^m-1) (z^m+1), and factorizes:

${\displaystyle z^{4m}+az^{2m}+1=(z^{2m}+{\sqrt {2-a}}z^{m}+1)(z^{2}m-sqrt{2-a}z^{m}+1)}$ 

where a is a real constant with |a| ≤ 2; this factorization can be done recursively. 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 - ω_n^j) 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_k 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).

Generalization to arbitrary radix

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 φ_{n, α}(z) for positive integers n and for α in [0,1) by:

${\displaystyle \phi _{n,\alpha }(z)=\left\{{\begin{matrix}z^{2n}-2\cos(2\pi \alpha )z^{n}+1&{\mbox{if }}0<\alpha <1\\\\z^{2n}-1&{\mbox{if }}\alpha =0\end{matrix}}\right.}$ 

Note that all of the polynomials that appear in the Bruun factorization above can be written in this form. These polynomials can be recursively factorized for a factor (radix) r as follows:

${\displaystyle \phi _{rm,\alpha }(z)=\left\{{\begin{matrix}\prod _{\ell =0}^{r-1}\phi _{m,(\alpha +\ell )/r}&{\mbox{if }}0<\alpha \leq 0.5\\\\\prod _{\ell =0}^{r-1}\phi _{m,(1-\alpha +\ell )/r}&{\mbox{if }}0.5<\alpha <1\\\\\prod _{\ell =0}^{r-1}\phi _{m,\ell /(2r)}&{\mbox{if }}\alpha =0\end{matrix}}\right.}$ 

---

References:

• Georg Bruun, "z-Transform DFT filters and FFTs," IEEE Trans. on Acoustics, Speech and Signal Processing (ASSP) 26 (1), 56-63 (1978).
• H. J. Nussbaumer, Fast Fourier Transform and Convolution Algorithms (Springer-Verlag: Berlin, 1990).
• Yuhang Wu, "New FFT structures based on the Bruun algorithm," IEEE Trans. ASSP 38 (1), 188-191 (1990)
• Jianping Chen and Henrik Sorensen, "An efficient FFT algorithm for real-symmetric data," Proc. ICASSP 5, 17-20 (1992).
• Rainer Storn, "Some results in fixed point error analysis of the Bruun-FTT [sic] algorithm," IEEE Trans. ASSP 41 (7), 2371-2375 (1993).
• Hideo Murakami, "Real-valued decimation-in-time and decimation-in-frequency algorithms," IEEE Trans. Circuits Syst. II: Analog and Digital Sig. Proc. 41 (12), 808-816 (1994).
• Hideo Murakami, "Real-valued fast discrete Fourier transform and cyclic convolution algorithms of highly composite even length," Proc. ICASSP 3, 1311-1314 (1996).