FFT Y TWIDDLES FACTORS

La primera aparición de un método eficiente para el cálculo de la transformada de Fourier (TDF), como muchos otros algoritmos, se debe a Gauss (Heideman et al., 1985), Gauss estaba interesado en ciertos cálculos astronómicos –un área recurrente para la aplicación de la FFT– para obtener la interpolación de las órbitas de algunos asteroides a partir de un conjunto de observaciones equiespaciadas.

Gauss notó que una serie de Fourier de N = N1N2 observaciones se podría descomponer como N2 TDF de largo N1 que se combinaban como N1 TDF de largo N2. El algoritmo de Gauss, desarrollado alrededor de 1805, pasó desapercibido porque fue publicado en latín, después de su muerte, en la recopilación de todas sus obras.

John W. Tukey concibió la idea básica del algoritmo para el cálculo rápido de la transformada de Fourier durante una reunión del Comité Científico Asesor del presidente Kennedy. Entre otros temas, el Comité estudiaba el desarrollo de técnicas para detección de pruebas nucleares en la Unión Soviética, ya que la ratificación de un tratado para la eliminación de las mismas dependía de la obtención de un método que permitiese detectarlas sin necesidad de visitar las plantas nucleares rusas. Una idea era analizar los registros sismológicos obtenidos mediante sismógrafos ubicados en el lecho marino. La cantidad de sensores y la longitud de los datos demandaban algoritmos que permitieran calcular las TDF rápidamente.

Richard Garwin, de IBM, era otro de los asistentes a la reunión; cuando Tukey le mostró su idea inmediatamente vislumbró el amplio rango de aplicaciones. Decidido a implementarla, se dirigió al centro de cómputos de IBM Research en Yorktown Heights, Nueva York, donde James W. Cooley era un miembro relativamente reciente del plantel. Según algunas versiones (Brigham, 1974), Cooley fue el encargado de programar el algoritmo porque, según sus propias palabras, “no tenía nada importante que hacer”. Según otras versiones (Burk et al., 2004), Garwin tuvo que insistir para que Cooley dejase de lado otros proyectos en los que estaba trabajando, y se dedicase a programar ese algoritmo.

Debido al secreto impuesto por el Comité, Garwin hizo creer a Cooley que el programa se utilizaría para calcular la periodicidad de las orientaciones del spin de los átomos en un cristal tridimensional de He3 (Cooley, 1992). En cualquier caso, Cooley lo implementó muy rápidamente, con la esperanza de sacarse ese proyecto de encima y poder seguir con su trabajo. Sin embargo, los pedidos de copias del programa y de instrucciones para ejecutarlo se empezaron a acumular.

En algún momento se vislumbró la posibilidad de patentar el algoritmo, pero finalmente se decidió que fuese de dominio público. En 1965 apareció en Mathematics of Computation el ahora famoso “An algorithm for the machine calculation of complex Fourier series” por Cooley y Tukey; fue publicado rápidamente (en menos de ocho meses), y el artículo junto con el proselitismo de Garwin hicieron que el algoritmo se difundiera muy rápidamente. Sin la tenacidad de Garwin, es posible que la FFT siguiese siendo desconocida aún hoy en día. Aunque según reza la introducción al algoritmo de la FFT en Numerical Recipes –un libro clásico sobre algoritmos computacionales– “si se acelera cualquier algoritmo no trivial por un factor de un millón todo el mundo buscará la forma de encontrarle aplicaciones útiles”.

Este trabajo fue un hito en la literatura sobre transformadas de Fourier, y se le ha atribuido (Oppenheim y Schafer, 1975; Rabiner y Rader, 1972) el tremendo incremento en las aplicaciones y el desarrollo del procesamiento digital de señales que comenzó alrededor de 1970.

Luego que Cooley y Tukey publicaran sus descubrimientos, comenzaron a aparecer reportes de otros investigadores que habían utilizado técnicas similares (Cooleyet al., 1967).

A su vez, este artículo citaba un par de trabajos de Runge (1903 y 1905); estos artículos, junto con las notas de curso de Runge y Konig (1964) describen esencialmente las ventajas computacionales del algoritmo de la FFT tal como se lo conoce hoy en día. El algoritmo de 1903 se concentraba en el caso en que N era potencia de 2, y el artículo de 1905 extendía la idea para potencias de 3. Fue utilizado ampliamente alrededor de 1940, de manera que fue bastante bien conocido. De todos modos, desapareció después de la Segunda Guerra Mundial.

L. H. Thomas, un estadístico del IBM Watson Laboratory también había desarrollado una técnica eficiente para el cálculo de la TDF (Thomas, 1963), que había encontrado en un libro de Stumpf (1939). Thomas generalizó los resultados de Stumpf y derivó una técnica muy similar a la de la FFT.

También Good, en 1960, extendió una idea de Yates (1937) para calcular la interacción de 2a experimentos, y obtuvo un procedimiento para calcular las transformadas de Fourier de N puntos que es esencialmente equivalente al de Thomas.

Respecto a los desarrollos posteriores, Yavne (1968) presentó un algoritmo –en un artículo poco difundido– que requiere el menor número conocido de multiplicaciones y de sumas para FFTs de longitud 2a , tanto para datos reales como complejos. Esta marca se mantuvo hasta hace muy poco tiempo, como se comenta más adelante, al menos para algoritmos prácticos. El mismo número de operaciones fue obtenido posteriormente por otros algoritmos más sencillos; debido al estilo críptico de Yavne pocos investigadores pudieron aplicar sus ideas en la época de su publicación.

En 1976 Rader y Brenner presentaron un algoritmo en el cual todos los factores eran reales o imaginarios puros (es decir, no necesitaba de los factores de rotación o twiddle factors), lo que permite disminuir sensiblemente el número de operaciones, a costa de aumentar el número de sumas y de una mayor sensitividad al ruido de redondeo.

Tomó casi 15 años la aparición de nuevos algoritmos para FFTs de largo 2a que requirieran tan pocas operaciones como el de Yavne. En 1984 cuatro trabajos se presentaron casi simultáneamente (Duhamel y Hollmann, 1984; Martens, 1984; Stasinski, 1984; Vetterli y Nussbaumer, 1984) que presentaban las bases de los métodos de “raíces partidas”, los que se basan en utilizar un algoritmo de raíz 2 para la parte par de la señal y uno de base 4 para la parte impar.

El trabajo de Winograd en teoría de complejidad (1977), particularmente sobre el número de multiplicaciones requeridas para calcular el producto o la convolución de polinomios, permitió aplicar los resultados previos de Good (1960) y de Rader (1968), quien mostró cómo calcular una FFT de largo N primo en una convolución circular de largo N  1. Si bien estos resultados fueron considerados como meras curiosidades al momento de su publicación, la combinación de los mismos, primero por Winograd (1976) y luego por Kolba y Parks (1977) despertó un gran interés en este tipo del algoritmos, que se conocen colectivamente como Winograd Fourier Transforms Algorithm (WFTA), y que requieren el menor número de multiplicaciones para TDFs de longitud moderada. La desventaja es que este procedimiento, lo mismo que uno muy similar conocido como algoritmo de factores primos, se basan en conceptos matemáticos muy avanzados, ocasionando cierta demora y esfuerzo para trasladar estos resultados teóricos a algoritmos computacionales eficientes.

La búsqueda de mejores prestaciones sigue en la actualidad. El último algoritmo conocido fue presentado en 2006 (Johnson y Frigo, 2006): es una variación del algoritmo de raíces partidas y según sus autores utiliza un 6 % menos de operaciones que el algoritmo de Yavne.

Resulta llamativo que los algoritmos anteriores al de Cooley y Tukey hayan desaparecido o permaneciesen ocultos, mientras que éste tuvo enorme repercusión. Una explicación posible es que el interés en los aspectos teóricos del procesamiento digital de señales fue creciendo junto con las mejoras tecnológicas en la industria de los semiconductores, ya que la disponibilidad de mayor potencia de cómputo a costo razonable permite desarrollar muchas aplicaciones novedosas.

Esquema de  sistema que incluyen un módulo con FFT:

Ell avance de la tecnología ha influido en el diseño del algoritmo: mientras que en los años 60 y principios del 70 el interés se centraba en disminuir el número de multiplicaciones, el progreso reveló que también son de importancia fundamental el número de sumas y de accesos a memoria (desde el punto de vista del software) y los costos de comunicación (desde el punto de vista del hardware).