Algoritmo FFT de Bruun

O algoritmo de Bruun é um algoritmo de transformada rápida de Fourier (FFT) baseado numa abordagem invulgar de fatorização recursiva de polinómios, proposto para potências de dois por G. Bruun em 1978 e generalizado para tamanhos compostos pares arbitrários por H. Murakami em 1996. Como as suas operações envolvem apenas coeficientes reais até à última etapa de cálculo, foi inicialmente proposto como uma forma de calcular eficientemente a transformada discreta de Fourier (DFT) de dados reais. O algoritmo de Bruun não tem tido um uso generalizado, contudo, uma vez que as abordagens baseadas no algoritmo FFT de Cooley-Tukey comum foram adaptadas com sucesso para dados reais com pelo menos a mesma eficiência. Além disso, existem evidências de que o algoritmo de Bruun possa ser intrinsecamente menos preciso do que o de Cooley-Tukey face à precisão numérica finita (Storn 1993).

No entanto, o algoritmo de Bruun ilustra uma estrutura algorítmica alternativa que pode expressar tanto a si próprio como ao algoritmo de Cooley-Tukey e, por conseguinte, fornece uma perspetiva interessante sobre as FFTs que permite misturas dos dois algoritmos e outras generalizações.

Recorde-se que a DFT é definida pela fórmula: $${\displaystyle X_{k}=\sum _{n=0}^{N-1}x_{n}e^{-{\frac {2\pi i}{N}}nk}\qquad k=0,\dots ,N-1.}$$

Por conveniência, denotemos as ${\displaystyle N}$ raízes da unidade por ${\displaystyle \omega _{N}^{n}}$ (${\displaystyle n=0,\dots ,N-1}$): $${\displaystyle \omega _{N}^{n}=e^{-{\frac {2\pi i}{N}}n}}$$ e defina-se o polinómio ${\displaystyle x(z)}$ cujos coeficientes são ${\displaystyle x_{n}}$: $${\displaystyle x(z)=\sum _{n=0}^{N-1}x_{n}z^{n}.}$$

A DFT pode então ser compreendida como uma redução deste polinómio; isto é, ${\displaystyle X_{k}}$ é dado por: $${\displaystyle X_{k}=x(\omega _{N}^{k})=x(z)\mod (z-\omega _{N}^{k})}$$ onde mod denota a operação de resto polinomial (teorema do resto polinomial). A chave para algoritmos rápidos como o de Bruun ou Cooley-Tukey advém do facto de que é possível executar este conjunto de ${\displaystyle N}$ operações de resto em fases recursivas.

Para calcular a DFT, precisamos de avaliar o resto de ${\displaystyle x(z)}$ módulo ${\displaystyle N}$ polinómios de grau 1, tal como descrito acima. Avaliar estes restos um a um é equivalente a avaliar a fórmula habitual da DFT diretamente, e requer ${\displaystyle O(N^{2})}$ operações. Contudo, é possível combinar estes restos recursivamente para reduzir o custo, utilizando o seguinte truque: se quisermos avaliar ${\displaystyle x(z)}$ módulo dois polinómios ${\displaystyle U(z)}$ e ${\displaystyle V(z)}$, podemos primeiro obter o resto módulo o seu produto ${\displaystyle U(z)V(z)}$, o que reduz o grau do polinómio ${\displaystyle x(z)}$ e torna as operações de módulo subsequentes menos dispendiosas computacionalmente.

O produto de todos os monómios ${\displaystyle (z-\omega _{N}^{k})}$ para ${\displaystyle k=0\dots N-1}$ é simplesmente ${\displaystyle z^{N}-1}$ (cujas raízes são claramente as ${\displaystyle N}$ raízes da unidade). Pretende-se então encontrar uma fatorização recursiva de ${\displaystyle z^{N}-1}$ em polinómios de poucos termos e de grau cada vez menor. Para calcular a DFT, toma-se ${\displaystyle x(z)}$ módulo cada nível desta fatorização, recursivamente, até se chegar aos monómios e ao resultado final. Se cada nível da fatorização dividir cada polinómio num número ${\displaystyle O(1)}$ (limitado por uma constante) de polinómios menores, cada um com um número ${\displaystyle O(1)}$ de coeficientes não nulos, então as operações de módulo para esse nível demoram um tempo ${\displaystyle O(N)}$; uma vez que haverá um número logarítmico de níveis, a complexidade global é ${\displaystyle O(N\log N)}$.

Mais explicitamente, suponha-se, por exemplo, que ${\displaystyle z^{N}-1=F_{1}(z)F_{2}(z)F_{3}(z)}$, e que ${\displaystyle F_{k}(z)=F_{k,1}(z)F_{k,2}(z)}$, e assim sucessivamente. O algoritmo FFT correspondente consistiria em primeiro calcular ${\displaystyle x_{k}(z)=x(z)\mod F_{k}(z)}$, depois calcular ${\displaystyle x_{k,j}(z)=x_{k}(z)\mod F_{k,j}(z)}$, e assim por diante, criando recursivamente cada vez mais polinómios de resto de grau cada vez menor até se chegar aos resultados finais de grau 0.

Além disso, desde que os fatores polinomiais em cada etapa sejam primos entre si (o que para polinómios significa que não têm raízes comuns), é possível construir um algoritmo dual revertendo o processo com o Teorema chinês dos restos.

O algoritmo básico de Bruun para potências de dois ${\displaystyle N=2^{n}}$ fatora ${\displaystyle z^{2^{n}}-1}$ recursivamente através das regras:

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

onde ${\displaystyle a}$ é uma constante real com ${\displaystyle |a|\leq 2}$. Se ${\displaystyle a=2\cos(\phi )}$, ${\displaystyle \phi \in (0,\pi )}$, então ${\displaystyle {\sqrt {2+a}}=2\cos {\tfrac {\phi }{2}}}$ e ${\displaystyle {\sqrt {2-a}}=2\cos({\tfrac {\pi }{2}}-{\tfrac {\phi }{2}})}$.

Na etapa ${\displaystyle s}$, ${\displaystyle s=0,1,2,n-1}$, o estado intermédio consiste em ${\displaystyle 2^{s}}$ polinómios ${\displaystyle p_{s,0},\dots ,p_{s,2^{s}-1}}$ de grau ${\displaystyle 2^{n-s}-1}$ ou inferior, em que: $${\displaystyle {\begin{aligned}p_{s,0}(z)&=p(z)\mod \left(z^{2^{n-s}}-1\right)&\quad &{\text{e}}\\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{aligned}}}$$

Pela construção da fatorização de ${\displaystyle z^{2^{n}}-1}$, os polinómios ${\displaystyle p_{s,m}(z)}$ codificam cada um ${\displaystyle 2^{n-s}}$ valores $${\displaystyle X_{k}=p(e^{2\pi i{\tfrac {k}{2^{n}}}})}$$ da transformada de Fourier; para ${\displaystyle m=0}$, os índices cobertos são ${\displaystyle k=0,2^{k},2\cdot 2^{s},3\cdot 2^{s},\dots ,(2^{n-s}-1)\cdot 2^{s}}$; para ${\displaystyle m>0}$ os índices cobertos são ${\displaystyle k=m,2^{s+1}-m,2^{s+1}+m,2\cdot 2^{s+1}-m,2\cdot 2^{s+1}+m,\dots ,2^{n}-m}$.

Durante a transição para a etapa seguinte, o polinómio ${\displaystyle p_{s,\ell }(z)}$ é reduzido aos polinómios ${\displaystyle p_{s+1,\ell }(z)}$ e ${\displaystyle p_{s+1,2^{s}-\ell }(z)}$ através de divisão polinomial. Se se pretender manter os polinómios por ordem crescente de índice, este padrão requer uma implementação com dois arrays (matrizes). Uma implementação in-place produz uma sequência de índices previsível, mas altamente desordenada; por exemplo, para ${\displaystyle N=16}$, a ordem final dos 8 restos lineares é (0, 4, 2, 6, 1, 7, 3, 5).

No final da recursividade, para ${\displaystyle s=n-1}$, restam ${\displaystyle 2^{n-1}}$ polinómios lineares que codificam dois coeficientes de Fourier ${\displaystyle X_{0}}$ e ${\displaystyle X_{2^{n-1}}}$ para o primeiro e, para qualquer outro ${\displaystyle k}$-ésimo polinómio, os coeficientes ${\displaystyle X_{k}}$ e ${\displaystyle X_{2^{n}-k}}$.

Em cada etapa recursiva, todos os polinómios do grau comum ${\displaystyle 4M-1}$ são reduzidos a duas partes de metade do grau ${\displaystyle 2M-1}$. O divisor deste cálculo de resto polinomial é um polinómio quadrático ${\displaystyle z^{m}}$, pelo que todas as reduções podem ser convertidas em divisões polinomiais de polinómios cúbicos por quadráticos. Existem ${\displaystyle N/2=2^{n-1}}$ destas pequenas divisões em cada etapa, originando um algoritmo ${\displaystyle O(N\log N)}$ para a FFT.

Além disso, como todos estes polinómios têm coeficientes puramente reais (até à última etapa), exploram automaticamente o caso especial em que as entradas ${\displaystyle x_{n}}$ são puramente reais para poupar sensivelmente um fator de dois no cálculo e no armazenamento. Também é possível tirar partido direto do caso de dados simétricos reais para calcular a transformada discreta de cosseno (Chen & Sorensen 1992).