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 de factorización polinómica recursiva inusual , propuesto para potencias de dos por G. Bruun en 1978 y generalizado a tamaños compuestos arbitrarios 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 se ha utilizado ampliamente, ya que los enfoques basados ​​en el algoritmo FFT de Cooley-Tukey ordinario 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 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 Cooley-Tukey y, por lo tanto, proporciona una perspectiva interesante sobre las FFT que permite mezclas de los dos algoritmos y otras generalizaciones.

Un enfoque polinomial para la DFT

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): y definamos 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 de resto del polinomio . La clave para algoritmos rápidos como el de Bruun o el de Cooley–Tukey proviene del hecho de que se puede realizar este conjunto de N operaciones de resto 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 polinomios módulo N de grado 1 como se describió anteriormente. Evaluar estos restos uno por uno es equivalente a evaluar 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 módulo dos polinomios y , primero podemos tomar el resto módulo su producto , lo que reduce el grado del polinomio y hace que las operaciones módulo posteriores sean menos costosas 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 simplemente (cuyas raíces son claramente las N raíces de la unidad). Luego 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 cada nivel de esta factorización por turno, 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) (constantemente 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 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 que , 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 de resto 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 del resto chino .

Cooley–Tukey como factorización polinómica

El algoritmo estándar de diezmado en frecuencia (DIF) de Cooley-Tukey de base r corresponde estrechamente a una factorización recursiva. Por ejemplo, el algoritmo de Cooley-Tukey de base 2 DIF factoriza en y . Estas operaciones módulo reducen el grado de en 2, lo que corresponde a dividir el tamaño del problema por 2. Sin embargo, en lugar de factorizar recursivamente de manera directa, Cooley-Tukey primero calcula x ₂ ( z ω _N ), desplazando todas las raíces (por un factor de giro ) de modo 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 no suele ser cierto para una factorización recursiva arbitraria (como la de Bruun, a continuación). ${\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 a través de 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 | a | ≤ 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 los polinomios y a través de la división polinómica. Si se desea mantener los polinomios en orden de índice creciente, este patrón requiere una implementación con dos matrices. Una implementación en su lugar produce una secuencia de índices predecible, pero altamente desordenada; por ejemplo, para N = 16, el orden final de los 8 residuos 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 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 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 cuadráticos. Hay N /2 = 2 ^{n −1} de estas pequeñas divisiones en cada etapa, lo que conduce 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 en el que las entradas x _n son puramente reales para ahorrar aproximadamente un factor de dos en el cálculo y el almacenamiento. También se puede aprovechar directamente el caso de datos reales simétricos para calcular la transformada de coseno discreta (Chen y Sorensen, 1992).

Generalización a bases arbitrarias

La factorización de Bruun, y por lo tanto el algoritmo FFT de Bruun, se generalizó para manejar longitudes compuestas pares arbitrarias, es decir, dividiendo el grado del polinomio por un radio arbitrario (factor), de la siguiente manera. Primero, definimos un conjunto de polinomios φ _{N , α} ( z ) para números enteros positivos N y para α en [0, 1) por: $${\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 de forma recursiva para un factor (raíz) 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). " Filtros DFT y FFT de transformada z ". IEEE Trans. on Acoustics, Speech and Signal Processing (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". IEEE Trans. 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 ]". IEEE Trans. Signal Process . 41 (7): 2371–2375.
• Murakami, Hideo (1994). "Algoritmos de decimación en tiempo y decimación en frecuencia con valores reales". IEEE Trans. Circuits Syst. II: Analog and Digital Sig. Proc . 41 (12): 808–816.
• Murakami, Hideo (1996). "Algoritmos de transformada de Fourier discreta rápida y convolución cíclica de valores reales 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. 7th Intl. Workshop, SAMOS 2007 (Samos, Grecia, 16-19 de julio de 2007).