Algoritmo FFT de Bruun

Algoritmo de transformada rápida de Fourier

El algoritmo de Bruun es un algoritmo de transformada rápida de Fourier (FFT) basado en un enfoque inusual de factorización polinómica recursiva , 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 implican 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, existe 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 expresarse a sí mismo y al algoritmo de Cooley-Tukey y, por lo tanto, proporciona una perspectiva interesante sobre las FFT que permite combinaciones de los dos algoritmos y otras generalizaciones.

Un enfoque polinomial para la DFT

Recuerde que la DFT está definida por 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.}$$

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

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

Factorizaciones recursivas y FFT

Para calcular la DFT, necesitamos evaluar el resto de los polinomios de grado 1 de módulo N como se describió anteriormente. Evaluar estos restos uno por uno equivale a evaluar directamente la fórmula DFT habitual y requiere operaciones O ( N ² ). Sin embargo, se pueden combinar estos restos de forma recursiva para reducir el costo, usando el siguiente truco: si queremos evaluar polinomios en módulo dos y , primero podemos tomar el resto en módulo su producto , lo que reduce el grado del polinomio y realiza operaciones de módulo posteriores. menos costoso computacionalmente. ${\displaystyle x(z)}$ ${\displaystyle x(z)}$ ${\displaystyle U(z)}$ ${\displaystyle V(z)}$ ${\displaystyle U(z)}$ ${\displaystyle V(z)}$ ${\displaystyle x(z)}$

El producto de todos los monomios para k =0.. N -1 es simple (cuyas raíces son claramente las N raíces de la unidad). Entonces se desea encontrar una factorización recursiva de en polinomios de pocos términos y de grado cada vez menor. Para calcular la DFT, se toma módulo de cada nivel de esta factorización por turno, de forma recursiva, hasta llegar a los monomios y al resultado final. Si cada nivel de factorización divide cada polinomio en un número O(1) (con límites constantes) 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 O( N ) tiempo ; dado que habrá un número logarítmico de niveles, la complejidad general es O ( N log N ). ${\displaystyle (z-\omega _{N}^{k})}$ ${\displaystyle z^{N}-1}$ ${\displaystyle z^{N}-1}$ ${\displaystyle x(z)}$

Más explícitamente, supongamos, por ejemplo , que y eso , etc. 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. ${\displaystyle z^{N}-1=F_{1}(z)F_{2}(z)F_{3}(z)}$ ${\displaystyle F_{k}(z)=F_{k,1}(z)F_{k,2}(z)}$

Además, siempre que los factores polinomiales 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 polinomial

El algoritmo estándar de diezmado en frecuencia (DIF) radix- r Cooley-Tukey se corresponde estrechamente con una factorización recursiva. Por ejemplo, radix-2 DIF Cooley–Tukey factoriza en y . Estas operaciones de módulo reducen el grado de por 2, lo que corresponde a dividir el tamaño del problema por 2. Sin embargo, en lugar de factorizar directamente de forma recursiva, Cooley-Tukey primero calcula x ₂ ( z ω _N ), desplazando todas las raíces (por un giro) . factor ) para que pueda aplicar la factorización recursiva de a ambos subproblemas. Es decir, Cooley-Tukey asegura que todos los subproblemas también sean DFT, mientras que esto generalmente no es cierto para una factorización recursiva arbitraria (como la de Bruun, más adelante). ${\displaystyle z^{N}-1}$ ${\displaystyle F_{1}=(z^{N/2}-1)}$ ${\displaystyle F_{2}=(z^{N/2}+1)}$ ${\displaystyle x(z)}$ ${\displaystyle F_{2}}$ ${\displaystyle F_{1}}$

La factorización de Bruun

El algoritmo básico de Bruun para potencias de dos N = 2 ⁿ factoriza z ^{2 n} - 1 recursivamente mediante las 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 | un | ≤ 2. Si , , entonces y . ${\displaystyle a=2\cos(\phi )}$ ${\displaystyle \phi \in (0,\pi )}$ ${\displaystyle {\sqrt {2+a}}=2\cos {\tfrac {\phi }{2}}}$ ${\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 de grado 2 ^{n - s} - 1 o menos, donde ${\displaystyle p_{s,0},\dots ,p_{s,2^{s}-1}}$ $${\displaystyle {\begin{aligned}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{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 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 . $${\displaystyle X_{k}=p(e^{2\pi i{\tfrac {k}{2^{n}}}})}$$

Durante la transición a la siguiente etapa, el polinomio se reduce a polinomios y mediante división polinomial. Si uno quiere mantener los polinomios en orden de índice creciente, este patrón requiere una implementación con dos matrices. Una implementación implementada 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). ${\displaystyle p_{s,\ell }(z)}$ ${\displaystyle p_{s+1,\ell }(z)}$ ${\displaystyle p_{s+1,2^{s}-\ell }(z)}$

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 del 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 del polinomio es un polinomio cuadrático z ^m , de modo que todas las reducciones pueden reducirse a divisiones polinómicas de polinomios cúbicos por 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), automáticamente explotan 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 simétricos reales para calcular la transformada discreta del coseno (Chen y Sorensen 1992).

Generalización a radices arbitrarias.

La factorización de Bruun, y por tanto el algoritmo FFT de Bruun, se generalizó para manejar longitudes compuestas pares arbitrarias, es decir, dividir el grado del polinomio por una base (factor) arbitraria , como sigue. Primero, definimos un conjunto de polinomios φ _{N , α} ( z ) para enteros positivos N y para α en [0, 1) mediante: $${\displaystyle \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}}}$$

Tenga en cuenta que todos los polinomios que aparecen en la factorización de Bruun anterior se pueden escribir de esta forma. Los ceros de estos polinomios son para en el caso y para en el caso. Por lo tanto, estos polinomios se pueden factorizar recursivamente para un factor (base) r mediante: ${\displaystyle e^{2\pi i(\pm \alpha +k)/N}}$ ${\displaystyle k=0,1,\dots ,N-1}$ ${\displaystyle \alpha \neq 0}$ ${\displaystyle e^{2\pi ik/2N}}$ ${\displaystyle k=0,1,\dots ,2N-1}$ ${\displaystyle \alpha =0}$

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

Referencias

• Bruun, Georg (1978). " z -Transformar filtros DFT y FFT". Traducción IEEE . de Acústica, Habla y Procesamiento de Señales (ASSP) . 26 (1): 56–63.
• Nussbaumer, HJ (1990). Algoritmos de convolución y transformada rápida de Fourier . Berlín: Springer-Verlag.
• Wu, Yuhang (1990). "Nuevas estructuras FFT basadas en el algoritmo de Bruun". Traducción IEEE. ASSP . 38 (1): 188-191.
• Chen, Jianping; Sorensen, Henrik (1992). "Un algoritmo FFT eficiente para datos simétricos reales". Proc. ICASSP . 5 : 17–20.
• Storn, Rainer (1993). "Algunos resultados en el análisis de errores de punto fijo del algoritmo Bruun-FTT [ sic ]". Traducción IEEE. Proceso de señal . 41 (7): 2371–2375.
• Murakami, Hideo (1994). "Algoritmos de diezmación en tiempo y frecuencia de valor real". Traducción IEEE. Sistema de circuitos. II: Señal Analógica y Digital. Proc . 41 (12): 808–816.
• Murakami, Hideo (1996). "Algoritmos de convolución cíclica y transformada de Fourier discreta rápida de valor real de longitud uniforme altamente compuesta". Proc. ICASSP . 3 : 1311-1314.
• Shashank Mittal, Md. Zafar Ali Khan, MB Srinivas, "Un estudio comparativo de diferentes arquitecturas FFT para radio definida por software", Lecture Notes in Computer Science 4599 ( Sistemas informáticos integrados: arquitecturas, modelado y simulación ), 375-384 (2007 ). Proc. Séptimo Internacional Taller, SAMOS 2007 (Samos, Grecia, 16 al 19 de julio de 2007).