Algoritmo FFT de Bruun

El algoritmo de Bruun es un algoritmo de transformada rápida de Fourier (FFT) basado en un enfoque inusual de factorización polinomial recursiva, propuesto para potencias de dos por G. Bruun en 1978 y generalizado a arbitrariamente incluso tamaños compuestos por H. Murakami en 1996. Debido a que sus operaciones involucran solo coeficientes reales hasta la última etapa de cálculo, inicialmente se propuso como una forma de calcular de manera eficiente 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, existe evidencia de que el algoritmo de Bruun puede ser intrínsecamente menos preciso que el de Cooley-Tukey frente a la precisión numérica finita (Storn, 1993).

Sin embargo, el algoritmo de Bruun ilustra un marco algorítmico alternativo que puede expresarse tanto a sí mismo como al algoritmo de Cooley-Tukey y, por lo tanto, brinda una perspectiva interesante sobre las FFT que permite mezclas de los dos algoritmos y otras generalizaciones.

Un enfoque polinomial de la DFT

Recuerde que la DFT se define mediante la fórmula:

Xk=.. n=0N− − 1xne− − 2π π iNnkk=0,...... ,N− − 1.{displaystyle X_{k}=sum ¿Qué?

Por comodidad, denotemos las raíces N de la unidad con ω_Nⁿ (n = 0,..., N − 1):

⋅ ⋅ Nn=e− − 2π π iNn{displaystyle omega ¿Qué? {2pi} {fn}n}

y define el polinomio x(z) cuyos coeficientes son x_n:

x()z)=.. n=0N− − 1xnzn.{displaystyle x(z)=sum ¿Qué?

La DFT puede entenderse entonces como una reducción de este polinomio; es decir, X_k viene dado por:

Xk=x()⋅ ⋅ Nk)=x()z)mod()z− − ⋅ ⋅ Nk){displaystyle X_{k}=x(omega) ¿Qué?

donde mod denota la operación de resto polinomial. La clave para algoritmos rápidos como Bruun's o Cooley–Tukey proviene del hecho de que uno puede realizar este conjunto de operaciones de resto N en etapas recursivas.

Factorizaciones recursivas y FFT

Para calcular el DFT, necesitamos evaluar el resto del x()z){displaystyle x(z)} modulo N polinomios grado-1 como se describe anteriormente. Evaluar estos restos uno por uno es equivalente a la evaluación de la fórmula DFT habitual directamente, y requiere O(N²) operaciones. Sin embargo, uno puede combinar estos restos recursivamente para reducir el costo, utilizando el siguiente truco: si queremos evaluar x()z){displaystyle x(z)} modulo dos polinomios U()z){displaystyle U(z)} y V()z){displaystyle V(z)}, primero podemos tomar el resto modulo su producto U()z){displaystyle U(z)} V()z){displaystyle V(z)}, que reduce el grado del polinomio x()z){displaystyle x(z)} y hace que las operaciones de modulo posteriores sean menos costosas.

El producto de todos los monomiales ()z− − ⋅ ⋅ Nk){displaystyle (z-omega) ¿Qué? para k=0..N-1 es simplemente zN− − 1{displaystyle z^{N}-1} (cuyas raíces son claramente N raíces de la unidad). Uno desea entonces encontrar una factorización recurrente de zN− − 1{displaystyle z^{N}-1} en polinomios de pocos términos y menor grado. Para calcular el DFT, se necesita x()z){displaystyle x(z)} modulo cada nivel de esta factorización a su vez, recursivamente, hasta que uno llega a los monomiales y el resultado final. Si cada nivel de la factorización divide cada polinomio en un número O(1) (continuado) de polinomios más pequeños, cada uno con un número O(1) de coeficientes no cero, entonces las operaciones de modulo para ese nivel toman O(N) tiempo; ya 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 zN− − 1=F1()z)F2()z)F3()z){displaystyle z^{N}-1=F_{1}(z)F_{2}(z)F_{3}(z)}, y eso Fk()z)=Fk,1()z)Fk,2()z){displaystyle F_{k}(z)=F_{k,1}(z)F_{k,2}(z)}, y así sucesivamente. El algoritmo FFT correspondiente consistía en el primer cálculo x_k()z) x()zMod F_k()zEntonces computación x_k,j()z) x_k()zMod F_k,j()z), y así sucesivamente, la creación de polinomios más y más restantes de menor y menor grado hasta que uno llega a los resultados finales de grado-0.

Además, siempre que los factores del polinomio 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.

Cooley–Tukey como factorización de polinomios

El radio estándar decimación en frecuencia (DIF)r El algoritmo Cooley-Tukey corresponde estrechamente a una factorización recursiva. Por ejemplo, radix-2 DIF Cooley–Tukey factors zN− − 1{displaystyle z^{N}-1} en F1=()zN/2− − 1){displaystyle F_{1}=(z^{N/2}-1)} y F2=()zN/2+1){displaystyle F_{2}=(z^{N/2}+1)}. Estas operaciones de modulo reducen el grado de x()z){displaystyle x(z)} por 2, que corresponde a dividir el tamaño del problema por 2. En lugar de la factorización recursiva F2{displaystyle F_{2} directamente, sin embargo, Cooley–Tukey primero computes x₂()z ⋅_N), cambiar todas las raíces (por twiddle factor) para que pueda aplicar la factorización recursiva de F1{displaystyle F_{1} a ambos subproblemas. Es decir, Cooley-Tukey asegura que todos los subproblemas son también DFT, mientras que esto no es generalmente cierto para una factorización recursiva arbitraria (como Bruun's, abajo).

La factorización de Bruun

Did you mean:

The basic Bruun algorithm for powers of two N=2^{n factorise z2n-1 recursively via the rules:}

z2M− − 1=()zM− − 1)()zM+1){displaystyle z^{2M}-1=(z^{M}-1)(z^{M}+1),}

z4M+az2M+1=()z2M+2− − azM+1)()z2M− − 2− − azM+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 a=2#⁡ ⁡ ()φ φ ){displaystyle a=2cos(phi)}, φ φ ▪ ▪ ()0,π π ){displaystyle phi in (0,pi)}, entonces 2+a=2#⁡ ⁡ φ φ 2{fnMicroc {fnK}=2cos {fnMicroc} } {2}} y 2− − a=2#⁡ ⁡ ()π π 2− − φ φ 2){fnMicroc {fnMicroc}=2cos {fnMicroc {fnMicroc} } {2}-{tfrac {phi }{2}}}.

En el escenario s, s=0,1,2,n-1, el estado intermedio consta de 2^s polinomios ps,0,...... ,ps,2s− − 1{fnMicrosoft Sans Serif} grado 2^n-s - 1 o menos donde

ps,0()z)=p()z)mod()z2n− − s− − 1)yps,m()z)=p()z)mod()z2n− − s− − 2#⁡ ⁡ ()m2sπ π )z2n− − 1− − s+1)m=1,2,...... ,2s− − 1{fnMicrosoft Sans Serif} {fnMicrosoft Sans Serif} {fnMicrosoft Sans Serif}

Por la construcción de la factorización de z²ⁿ-1, los polinomios p_s,m(z) cada uno codifica 2^n-s valores

Xk=p()e2π π ik2n){displaystyle ¿Qué? }

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 ps,l l ()z){displaystyle p_{s,ell }(z)} se reduce a los polinomios ps+1,l l ()z){displaystyle p_{s+1,ell }(z)} y ps+1,2s− − l l ()z){displaystyle p_{s+1,2^{s}-ell }(z)} vía división polinomio. Si uno quiere mantener los polinomios en orden índice creciente, este patrón requiere una implementación con dos arrays. Una implementación en su lugar produce una secuencia predecible, pero altamente no ordenada de índices, por ejemplo, N=16 la orden final del 8 restos lineales es (0, 4, 2, 6, 1, 7, 3, 5).

Al final de la recursividad, 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 4M-1 se reducen a dos partes de la mitad del grado 2M -1. El divisor de este cálculo del resto polinomial 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 lleva a un algoritmo O (N log N) para la FFT.

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

Generalización a radices arbitrarias

(feminine)

La factorización de Bruun, y por lo tanto el algoritmo FFT de Bruun, se generalizó para manejar longitudes compuestas pares arbitrarias, es decir, dividir el grado del polinomio por una base arbitraria (factor), como sigue. Primero, definimos un conjunto de polinomios φ_N,α(z) para enteros positivos N y para α en [0,1) por:

<math alttext="{displaystyle phi _{N,alpha }(z)=left{{begin{matrix}z^{2N}-2cos(2pi alpha)z^{N}+1&{mbox{if }}0<alpha φ φ N,α α ()z)={}z2N− − 2#⁡ ⁡ ()2π π α α )zN+1si0.α α .1z2N− − 1siα α =0{displaystyle phi _{N,alpha }(z)=left{begin{matrix}z^{2N}-2cos(2pialpha)z^{N}+1 {mbox{if }0 seleccionóalpha <1\\z^{2N}-1} {mbox{if}}f}f}f}f}f}f}f}fn0}f}fn0}f}<img alt="phi _{{N,alpha }}(z)=left{{begin{matrix}z^{{2N}}-2cos(2pi alpha)z^{N}+1&{mbox{if }}0<alpha

Tenga en cuenta que todos los polinomios que aparecen en la factorización Bruun arriba se pueden escribir en esta forma. Los ceros de estos polinomios son e2π π i()± ± α α +k)/N{displaystyle e^{2pi i(pm alpha +k)/N} para k=0,1,...... ,N− − 1{displaystyle k=0,1,dotsN-1} en el α α ل ل 0{displaystyle alpha neq 0} caso, y e2π π ik/2N{displaystyle e^{2pi ik/2N} para k=0,1,...... ,2N− − 1{displaystyle k=0,1,dots2N-1} en el α α =0{displaystyle alpha =0} caso. Por lo tanto estos polinomios pueden ser factorizados recursivamente para un factor (radix) r via:

<math alttext="{displaystyle phi _{rM,alpha }(z)=left{{begin{array}{ll}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 φ φ rM,α α ()z)={}∏ ∏ l l =0r− − 1φ φ M,()α α +l l )/rsi0.α α ≤ ≤ 0.5∏ ∏ l l =0r− − 1φ φ M,()1− − α α +l l )/rsi0.5.α α .1∏ ∏ l l =0r− − 1φ φ M,l l /()2r)siα α =0{displaystyle phi _{rM,alpha }(z)=left{begin{array}{ll}prod ################################################################################################################################################################################################################################################################ ################################################################################################################################################################################################################################################################ ################################################################################################################################################################################################################################################################<img alt="phi _{{rM,alpha }}(z)=left{{begin{array}{ll}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