Algoritmo FFT de Bruun

El algoritmo de Bruun es un algoritmo de transformada rápida de Fourier (FFT) basado en un enfoque de factorización polinomial recursiva inusual , propuesto para potencias de dos por G. Bruun en 1978 y generalizado a tamaños compuestos pares arbitrarios por H. Murakami en 1996. Debido a que sus operaciones involucran solo coeficientes reales hasta la última etapa de cálculo, se propuso inicialmente como una forma de calcular eficientemente la transformada discreta de Fourier (DFT) de datos reales. Sin embargo, el algoritmo de Bruun no ha tenido un uso generalizado, ya que los enfoques basados ​​en el algoritmo FFT ordinario de Cooley-Tukey se han adaptado con éxito a datos reales con al menos la misma eficiencia. Además, hay evidencia de que el algoritmo de Bruun puede ser intrínsecamente menos preciso que Cooley-Tukey frente a una precisión numérica finita ( Storn 1993 ) .

Sin embargo, el algoritmo de Bruun ilustra un marco algorítmico alternativo que puede expresar tanto el algoritmo de Cooley-Tukey como a sí mismo, y por lo tanto proporciona una perspectiva interesante sobre las FFT que permite combinaciones de ambos algoritmos y otras generalizaciones.

Recordemos que la DFT se define mediante la 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.}$$

Para mayor comodidad, denotemos las N raíces de la unidad por ω _N ⁿ ( n  =  0,  ..., N − 1):    $${\displaystyle \omega _ {N}^{n}=e^{-{\frac {2\pi i}{N}}n}}$$ y definimos el polinomio x ( z ) cuyos coeficientes son x _n : $${\displaystyle x(z)=\sum _ {n=0}^{N-1}x_ {n}z^{n}.}$$

La DFT puede entenderse entonces como una reducción de este polinomio; es decir, X _k viene dada por: $${\displaystyle X_{k}=x(\omega _{N}^{k})=x(z)\mod (z-\omega _{N}^{k})}$$ donde mod denota la operación de resto polinomial . La clave de algoritmos rápidos como el de Bruun o el de Cooley-Tukey reside en el hecho de que este conjunto de N operaciones de resto se puede realizar en etapas recursivas.

Para calcular la DFT, necesitamos evaluar el resto de${\displaystyle x(z)}$módulo N polinomios de grado 1 como se describió anteriormente. Evaluar estos restos uno por uno es equivalente a evaluar directamente la fórmula DFT usual y requiere O( N² ) operaciones. ^Sin embargo, se pueden combinar estos restos recursivamente para reducir el costo, utilizando el siguiente truco: si queremos evaluar${\displaystyle x(z)}$módulo dos polinomios${\displaystyle U(z)}$y${\displaystyle V(z)}$, primero podemos tomar el resto módulo su producto${\displaystyle U(z)}$${\displaystyle V(z)}$lo cual reduce el grado del polinomio${\displaystyle x(z)}$y hace que las operaciones de módulo posteriores sean menos costosas desde el punto de vista computacional.

El producto de todos los monomios${\displaystyle (z-\omega _ {N}^{k})}$para k =0.. N -1 es simplemente${\displaystyle z^{N}-1}$(cuyas raíces son claramente las N raíces de la unidad). Entonces se desea encontrar una factorización recursiva de${\displaystyle z^{N}-1}$en polinomios de pocos términos y de grado cada vez menor. Para calcular la DFT, se toma${\displaystyle x(z)}$módulo cada nivel de esta factorización a su vez, recursivamente, hasta llegar a los monomios y al resultado final. Si cada nivel de la factorización divide cada polinomio en un número O(1) (constante-acotado) de polinomios más pequeños, cada uno con un número O(1) de coeficientes distintos de cero, entonces las operaciones de módulo para ese nivel toman un tiempo O( N ); dado que habrá un número logarítmico de niveles, la complejidad general es O( N log N ).

Más explícitamente, supongamos, por ejemplo, que${\displaystyle z^{N}-1=F_{1}(z)F_{2}(z)F_{3}(z)}$y que${\displaystyle F_{k}(z)=F_{k,1}(z)F_{k,2}(z)}$y así sucesivamente. El algoritmo FFT correspondiente consistiría en calcular primero x _k ( z ) = x ( z ) mod F _k ( z ), luego calcular x _{k , j} ( z ) = x _k ( z ) mod F _{k , j} ( z ), y así sucesivamente, creando recursivamente más y más polinomios restantes de grado cada vez menor hasta llegar a los resultados finales de grado 0.

Además, siempre que los factores polinómicos en cada etapa sean relativamente primos (lo que para los polinomios significa que no tienen raíces comunes), se puede construir un algoritmo dual invirtiendo el proceso con el teorema chino del resto .

El algoritmo básico de Bruun para potencias de dos N = 2 ⁿ factoriza z ^{2 n} - 1 recursivamente mediante las siguientes reglas:

$${\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)}$$

donde a es una constante real con | a | ≤ 2. Si${\displaystyle a=2\cos(\phi )}$,${\displaystyle \phi \in (0,\pi )}$, entonces${\displaystyle {\sqrt {2+a}}=2\cos {\tfrac {\phi }{2}}}$y${\displaystyle {\sqrt {2-a}}=2\cos({\tfrac {\pi }{2}}-{\tfrac {\phi }{2}})}$.

En la etapa s , s = 0, 1, 2, n - 1, el estado intermedio consta de 2 ^s polinomios${\displaystyle p_{s,0},\dots ,p_{s,2^{s}-1}}$de grado 2 ^{n - s} - 1 o menor, donde $${\displaystyle {\begin{aligned}p_{s,0}(z)&=p(z)\mod \left(z^{2^{ns}}-1\right)&\quad &{\text{y}}\\p_{s,m}(z)&=p(z)\mod \left(z^{2^{ns}}-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}}}$$

Mediante la construcción de la factorización de z ^{2 n} - 1 , los polinomios p _{s , m} ( z ) codifican cada uno 2 ^{n - s} valores. $${\displaystyle X_{k}=p(e^{2\pi i{\tfrac {k}{2^{n}}}})}$$ de la transformada de Fourier, para m =0, los índices cubiertos son k = 0 , 2 ^k , 2∙2 ^s , 3∙2 ^s ,..., (2 ^{n - s} -1)∙2 ^s , para m > 0 los índices cubiertos son k = m , 2 ^{s +1} - m , 2 ^{s +1} + m , 2∙2 ^{s +1} - m , 2∙2 ^{s +1} + m , ..., 2 ⁿ - m .

Durante la transición a la siguiente etapa, el polinomio${\displaystyle p_{s,\ell }(z)}$se reduce a los polinomios${\displaystyle p_{s+1,\ell }(z)}$y${\displaystyle p_{s+1,2^{s}-\ell }(z)}$mediante división polinómica. Si se desea mantener los polinomios en orden ascendente de índice, este patrón requiere una implementación con dos arreglos. Una implementación in situ produce una secuencia de índices predecible, pero muy desordenada; por ejemplo, para N = 16, el orden final de los 8 restos lineales es (0, 4, 2, 6, 1, 7, 3, 5).

Al final de la recursión, para s = n -1 , quedan 2 ^{n -1} polinomios lineales que codifican dos coeficientes de Fourier X ₀ y X _{2 ^{n -1}} para el primero y para cualquier otro k -ésimo polinomio los coeficientes X _k y X _{2 ⁿ - k} .

En cada etapa recursiva, todos los polinomios de grado común 4 M -1 se reducen a dos partes de la mitad del grado 2 M -1 . El divisor de este cálculo del resto polinómico es un polinomio cuadrático z ^m , de modo que todas las reducciones se pueden reducir a divisiones polinómicas de polinomios cúbicos por polinomios cuadráticos. Hay N /2 = 2 ^{n −1} de estas pequeñas divisiones en cada etapa, lo que da como resultado un algoritmo O ( N log N ) para la FFT.

Además, dado que todos estos polinomios tienen coeficientes puramente reales (hasta la última etapa), aprovechan automáticamente el caso especial en el que las entradas x _n son puramente reales para ahorrar aproximadamente un factor de dos en cálculo y almacenamiento. También se puede aprovechar directamente el caso de datos reales simétricos para calcular la transformada discreta del coseno ( Chen y Sorensen, 1992 ) .