Editing Bruun's FFT algorithm (section)

== The Bruun factorization ==

The basic Bruun algorithm for [[power of two|powers of two]] ''N''=''2''<sup>''n''</sup> factorizes ''z''<sup>''2''<sup>''n''</sup></sup>-''1'' recursively via the rules:

<math display="block">z^{2M}-1 = (z^M - 1) (z^M + 1) \,</math>
<math display="block">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 a real constant with |''a''| ≤ 2. If <math>a=2\cos(\phi)</math>, <math>\phi\in(0,\pi)</math>, then <math>\sqrt{2+a}=2\cos\tfrac\phi2</math> and <math>\sqrt{2-a}=2\cos(\tfrac\pi 2-\tfrac\phi 2)</math>.

At stage ''s'', ''s''=0,1,2,''n''-1, the intermediate state consists of ''2''<sup>''s''</sup> polynomials <math>p_{s,0},\dots,p_{s,2^s-1}</math> of degree ''2''<sup>''n''-''s''</sup> - ''1'' or less , where
<math display="block">\begin{align}
p_{s,0}(z)&= p(z) \mod \left(z^{2^{n-s}}-1\right)&\quad&\text{and}\\
p_{s,m}(z) &= p(z)\mod \left(z^{2^{n-s}}-2\cos\left(\tfrac{m}{2^s}\pi\right)z^{2^{n-1-s}}+1\right)&m&=1,2,\dots,2^s-1
\end{align}</math>

By the construction of the factorization of ''z''<sup>''2''<sup>''n''</sup></sup>-''1'', the polynomials ''p''<sub>''s'',''m''</sub>(''z'') each encode 2<sup>''n''-''s''</sup> values
<math display="block">X_k=p(e^{2\pi i\tfrac{k}{2^n}})</math>
of the Fourier transform, for ''m''=0, the covered indices are ''k''=''0'', 2<sup>''k''</sup>, 2∙2<sup>''s''</sup>, 3∙2<sup>''s''</sup>,..., (2<sup>''n''-''s''</sup>-1)∙2<sup>''s''</sup>, for ''m''>''0'' the covered indices are ''k''=''m'', 2<sup>''s''+1</sup>-''m'', 2<sup>''s''+1</sup>+''m'', 2∙2<sup>''s''+1</sup>-''m'', 2∙2<sup>''s''+1</sup>+''m'', ..., 2<sup>''n''</sup>-''m''.

During the transition to the next stage, the polynomial <math>p_{s,\ell}(z)</math> is reduced to the polynomials <math>p_{s+1,\ell}(z)</math> and <math>p_{s+1,2^s-\ell}(z)</math> via polynomial division. If one wants to keep the polynomials in increasing index order, this pattern requires an implementation with two arrays. An implementation in place produces a predictable, but highly unordered sequence of indices, for example for ''N''=16 the final order of the 8 linear remainders is (0, 4, 2, 6, 1, 7, 3, 5).

At the end of the recursion, for {{math|1=''s'' = ''n''-1}}, there remain {{math|2<sup>''n''-1</sup>}} linear polynomials encoding two Fourier coefficients {{math|''X''<sub>''0''</sub>}} and {{math|''X''<sub>2<sup>''n''-1</sup></sub>}} for the first and for the any other {{mvar|k}}th polynomial the coefficients {{math|''X''<sub>''k''</sub>}} and {{math|''X''<sub>2<sup>''n''</sup>-''k''</sub>}}.

At each recursive stage, all of the polynomials of the common degree {{math|4''M''-1}} are reduced to two parts of half the degree {{math|2''M''-1}}. The divisor of this polynomial remainder computation is a quadratic polynomial ''z''<sup>''m''</sup>, so that all reductions can be reduced to polynomial divisions of cubic by quadratic polynomials. There are {{math|1=''N''/2 = 2<sup>''n''−1</sup>}} of these small divisions at each stage, leading to an {{math|''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''<sub>''n''</sub> 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]] {{harv|Chen|Sorensen|1992}}.

=== 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 {{math|''φ''<sub>''N'',''α''</sub>(''z'')}} for positive integers {{mvar|N}} and for {{mvar|α}} in {{closed-open|0, 1}} by:
<math display="block">\phi_{N, \alpha}(z) = \begin{cases}
z^{2N} - 2 \cos (2 \pi \alpha) z^N + 1 & \text{if } 0 < \alpha < 1 \\ \\
z^{2N} - 1 & \text{if } \alpha = 0
\end{cases} </math>

Note that all of the polynomials that appear in the Bruun factorization above can be written in this form. The zeroes of these polynomials are <math>e^{2\pi i ( \pm\alpha + k ) / N}</math> for <math>k = 0, 1, \dots, N-1</math> in the <math>\alpha \neq 0</math> case, and <math>e^{2\pi i k / 2N}</math> for <math>k = 0, 1, \dots, 2N-1</math> in the <math>\alpha = 0</math> case. Hence these polynomials can be recursively factorized for a factor (radix) {{mvar|r}} via:

<math display="block">\phi_{rM, \alpha}(z) = \begin{cases}
\prod_{\ell=0}^{r-1} \phi_{M,(\alpha+\ell)/r} & \text{if } 0 < \alpha \leq 0.5 \\ \\
\prod_{\ell=0}^{r-1} \phi_{M,(1-\alpha+\ell)/r} & \text{if } 0.5 < \alpha < 1 \\ \\
\prod_{\ell=0}^{r-1} \phi_{M,\ell/(2r)} & \text{if } \alpha = 0
\end{cases} </math>