Fourier (y III)

Esta entrada es una continuación de las entradas Fourier (I) y Fourier (II), donde hemos hablado sucesivamente de Fourier, las series y transformada que llevan su nombre y algunas aplicaciones.

La transformada discreta

Recordemos: la transformada de Fourier es un útil matemático que permite asociar a una función que depende de una variable física (del tiempo, de una o varias variables espaciales) otra función que depende de la frecuencia, es decir, que indica con qué frecuencia se presenta cada uno de los valores alcanzados.

Una versión finito-dimensional (y por tanto manejable desde un punto de vista computacional) de la transformada de Fourier es la transformada discreta.

Se trata de una aplicación de \({\bf R}^N\) en sí mismo que a cada \((x^0,x^1,\dots,x^{N-1})\) asigna los \(X^k\) siguientes:
\[
X^k = \sum_{n=0}^{N-1} x^n e^{-2 \pi i kn / N} , \quad k = 0,1, \dots N-1 .
\]

Se puede mirar como una versión aproximada (o discretizada) de la transformada habitual con un número de muestras igual a \(N\).

En efecto, consideremos un intervalo \([0,N]\) y una función \(f : [0,N] \to {\bf R}\) que es extendida por valores nulos a todo \({\bf R}\). Su transformada de Fourier \(\hat f\) está dada por
\[
\hat{f}(s) = \int_{-\infty}^{+\infty} f(t) e^{-2\pi i ts} \,dt \approx \sum_{n=0}^{N-1} f(n) e^{-2 \pi i s n / N}.
\]
Por tanto, los valores \(X^k\) asociados a los \(x^n = f(n)\) se pueden interpretar como aproximaciones de los \(\hat{f}(s_k)\) con \(s_k = k\).

En lo que sigue, identificaremos cada \(N\)-tupla \(\{ x^0, x^1, \dots, x^{N-1} \}\) de \({\bf R}^N\) con la función \(f_N : [0,N] \to {\bf R}\) constante a trozos que toma los valores \(f(n) = x^n\) en cada \(n = 0,1,\dots N-1\). Y denotaremos \({\mathcal F}_N\) la aplicación que a cada \(f_N\) de este tipo asigna la función análoga que, en cada \(k = 0,1,\dots N-1\), toma el valor \(X^k\).

He aquí algunas propiedades de la transformada discreta de Fourier cuyas demostraciones son sencillas:

• Linealidad: \({\mathcal F}_N (\lambda f_N + \mu g_N) = \lambda {\mathcal F}_N f_N + \mu {\mathcal F}_N g_N\) para todo \(\lambda\), \(f_N\), \(\mu\) y \(g_N\).

• Inversión: Si \(g_N = {\mathcal F}_N f_N\), entonces \(f_N(s) = N^{-1} {\mathcal F}_N g_N (-s)\).

• Identidad de Parseval: miradas \(f_N\), \(g_N\) y sus transformadas discretas como vectores de \({\bf R}^N\), se tiene que \({\mathcal F}_N f_N \cdot {\mathcal F}_N g_N = N f_N \cdot g_N\) para todas las \(f_N\) y \(g_N\).

• Convolución: dadas \(f_N\) y \(g_N\), llamaremos producto de convolución de \(f_N\) y \(g_N\) a la nueva función discreta \(h_N = f_N *_{\!N} g_N\), con\[h_N(n) = \sum_{m=0}^{N-1} f_N(m) g_N(n – m) \ \text{ para } \ n = 0,1,\dots N-1.\] Entonces se tiene siempre que \({\mathcal F}_N (f_N *_{\!N} g_N)(s) \equiv ({\mathcal F}_N f_N)(s) ({\mathcal F}_N g_N)(s)\). En otras palabras, la transformada discreta del producto de convolución de dos funciones discretas es el producto punto a punto (frecuencia a frecuencia) de sus respectivas transformadas discretas.

La transformada discreta de Fourier se puede utilizar con éxito en muchas aplicaciones:

1. En procesamiento de señales digitales, con \(f_N\) una señal que varía en el tiempo: la presión de una onda sonora, una señal radiofónica, los valores de la temperatura a lo largo de un intervalo de tiempo finito, etc.

2. Para describir la manera en que la luz viaja a través de múltiples sistemas ópticos. En particular, esto conduce a la reconstrucción tomográfica de objetos tridimensionales, con evidentes aplicaciones en medicina.

3. En operaciones complicadas, como puede ser el producto de dos polinomios de grado elevado. Para ello, se usan las propiedades de inversión y convolución enunciadas más arriba.

4. Un subproducto del uso que precede permite multiplicar con rapidez enteros grandes. En efecto, los enteros pueden ser tratados como valores tomados por adecuados polinomios, cuyos coeficientes están “bien preparados” para el cálculo. Combinado este método con el algoritmo descrito en la sección que sigue, aparecen los métodos más rápidos de los que se dispone para este fin.

5. Si aceptamos que los datos han sido “contaminados” por un ruido, esto es, son de la forma \(f_N *_{\!N} g_N\), donde \(g_N\) es una función regularizante convenientemente discretizada y convertida en una \(N\)-tupla \(\{ y_0, y_1, \dots, y_{N-1} \}\), la determinación de \(f_N\) puede realizarse fácilmente usando las propiedades de convolución e inversión: primero calculamos \({\mathcal F}_N (f_N *_{\!N} g_N) = {\mathcal F}_N f_N {\mathcal F}_N g_N\), después \({\mathcal F}_N f_N\) y finalmente \(f_N\).

Por tanto, queda claro que, dada \(f_N\), es decir, dados los \(x^n\), un método rápido de cálculo de \({\mathcal F}_N f_N\), es decir, de los \(X^k\), es muy deseable.

La transformada rápida de Fourier

Desgraciadamente, para grandes valores de \(N\), el cálculo directo de los \(X^k\) exige un número excesivo de operaciones, del orden de \(N \times N\).

En 1965, J.M. Cooley y J. Tuckey [1] introdujeron un algoritmo de cálculo que, explotando la forma particular que tienen las sumas anteriores, permite reducir considerablemente esta cantidad. Más precisamente, su aplicación necesita del orden de \(N \log N\) operaciones. Este algoritmo recibe desde entonces en la literatura la denominación “transformada rápida de Fourier” (fast Fourier transform).

Para comprender bien cómo funciona la transformada rápida, veamos qué se puede hacer en el caso particular en que \(N = 4\), esto es,
\[
X^k = \sum_{n=0}^{N-1} x^n w^n[k] , \ \text{ con } \ w^n[k] := e^{- \pi i kn / 2} .
\]

En principio, el número de operaciones necesarias es del orden de \(N \times N\). Re-escribiendo estas igualdades en forma matricial, tenemos:
\[ \displaystyle
\left(\begin{array}{c} X^0 \\ X^1 \\ X^2 \\ X^3 \end{array}\right)
= \left(\begin{array}{cccc}
w^0[0] & w^1[0] & w^2[0] & w^3[0] \\
w^0[1] & w^1[1] & w^2[1] & w^3[1] \\
w^0[2] & w^1[2] & w^2[2] & w^3[2] \\
w^0[3] & w^1[3] & w^2[3] & w^3[3]
\end{array}\right)
\left(\begin{array}{c} x^0 \\ x^1 \\ x^2 \\ x^3 \end{array}\right) .
\]

Pero la matriz precedente coincide con
\[
\left(\begin{array}{cccc}
1 & 1 & 1 & 1 \\
1 & w^1[1] & w^2[1] & w^3[1] \\
1 & w^2[1] & 1 & w^2[1] \\
1 & w^3[1] & w^2[1] & w^1[1]
\end{array}\right) .
\]

Así, después de intercambiar la segunda y tercera fila, el sistema se puede escribir en la forma
\[ \displaystyle
\left(\begin{array}{c} X^0 \\ X^2 \\ X^1 \\ X^3 \end{array}\right)
= \left(\begin{array}{cccc}
1 & 1 & 1 & 1 \\
1 & w^2[1] & 1 & w^2[1] \\
1 & w^1[1] & w^2[1] & w^3[1] \\
1 & w^3[1] & w^2[1] & w^1[1]
\end{array}\right)
\left(\begin{array}{c} x^0 \\ x^1 \\ x^2 \\ x^3 \end{array}\right) .
\]
y la estructura tan particular de la matriz de coeficientes \(A\) permite la factorización que sigue:
\[
A
= \left(\begin{array}{cccc}
1 & 1 & 0 & 0 \\
1 & w^2[1] & 0 & 0 \\
0 & 0 & 1 & w^1[1] \\
0 & 0 & 1 & w^3[1]
\end{array}\right)
\left(\begin{array}{cccc}
1 & 0 & 1 & 0 \\
0 & 1 & 0 & 1 \\
1 & 0 & w^2[1] & 0 \\
0 & 1 & 0 & w^2[1]
\end{array}\right).
\]

Usando esta factorización, vemos que el número de operaciones que precisa el cálculo de las \(X^k\) es \(4\) productos y \(8\) sumas, es decir, \(12\) operaciones, en vez de \(N \times N = 16\).

Más generalmente, si (por ejemplo) \(N = 2^p\), un razonamiento similar pero más complejo permite (tras reordenar las filas del sistema) factorizar la matriz de coeficientes en \(p\) matrices huecas, de manera que se reduce considerablemente el número de operaciones. Un análisis detallado muestra que este es
\[
\frac{N^2}{2^p} + pN = N(1 + \log_2 N)
\]
y, por tanto, como hemos anunciado, es del orden de \(N \log N\).

Por ejemplo, si \(N=1\,024\), la determinación directa de la transformada discreta requeriría \(N^2=1\,048\,576\) operaciones; mientras que la transformada rápida sólo precisa \(N \log_2 N = 10\,240\) (aproximadamente la centésima parte).

Cuando \(N\) no coincide exactamente con una potencia de \(2\), modificando adecuadamente el método, se llega a una conclusión análoga; véase [2].

Naturalmente, es posible hablar de la transformada discreta de Fourier en dimensión \(d \geq 2\) y de la correspondiente transformada rápida. En este caso, aceptando de nuevo que el número de muestras es \(N\), se puede probar que el número de operaciones necesarias es del orden de \(Nd \log (Nd)\).

Se ha dicho con frecuencia que la transformada rápida ha tenido tanto impacto que puede ser considerada una de las herramientas numéricas más útiles jamás inventadas. De acuerdo con lo que se ha mencionado más arriba, es muy útil en análisis espectral; convolución rápida; codificación de fuentes, habla e imagen; búsqueda de raíces de ecuaciones polinómicas, cálculo del producto de grandes enteros, etc.

Para saber más: algunas referencias

[1] J.M. Cooley y J. Tuckey, An algorithm for the machine calculation of complex Fourier series. Mathematics of Computation. 19 (90): 297?301 (1965).

[2] B. Osgood, Lectures on the Fourier transform and its applications. Pure and Applied Undergraduate Texts, 33. American Mathematical Society, Providence, RI, 2019.

[3] V. Serov, Fourier Series, Fourier Transform and their Applications to Mathematical Physics, Applied Mathematical Sciences, Vol. 197, Springer International Publishing, 2017.

[4] https://vibrationresearch.com/blog/fast-fourier-transform-fft-analysis/