lunes, 12 de enero de 2015

Fourier y sus senos

Hoy es #lunesTetas, así que supongo que ya sabréis lo que toca hablar, ¿no? Efectivamente, de esas curvas tan sensuales, siempre bamboleándose, pa'rriba y pa'bajo, pa'rriba y pa'bajo. Lo habéis adivinado: hoy toca senos
Pero no ese tipo de senos. Más bien estos otros
Sí, ya lo sé, me he ido por la tangente (qué chispa que tengo, eh?). Pero ya que has hecho el esfuerzo de llegar hasta aquí, sigue leyendo. Te advierto que es condición sine qua non para un rato de diversión con la trigonometría (tetas y chistes malos, qué mejor combinación).

A diferencia de lo que ocurre con la señorita tan maja de ahí arriba, los senos de las matemáticas (y los cosenos de los matemáticos, no vayan a decir que discriminamos a nadie) tienen la curiosa propiedad de que se repiten exactamente de forma periódica e indefinida.
La Maja Periódica es seguramente la obra más desconocida de la serie
y es que \(\sin (\theta + 2 \pi n) = \sin \theta\) y \(\cos(\theta + 2\pi n) = \cos \theta\) para cualquier número entero \(n\).


Desarrollador en serie

La importancia que tiene esto es que muchos fenómenos exhiben algún tipo de periodicidad (desde las horas de luz al día a la estructura de un cristal), y nos gustaría modelarlo con unas funciones tan sencillas como los senos y los cosenos.
Así que, si tenemos una función que se repita cada \(2\pi\), será o \(\sin(\theta)\) o \(\cos(\theta)\). ¡Pues no! Porque también tienen la misma periodicidad \(\sin(2\theta)\), \(\sin(3\theta)\), y en general \(\sin(k\theta)\) (y también \(\cos(k\theta)\) con \(k\) entero. Así pues, la forma más general para una función \(f(\theta)\) tal que \(f(\theta+2\pi n) = f(\theta)\) está dado por una serie de Fourier:\[f(\theta) = \sum_{k=0}^\infty[a_k \sin(k\theta) + b_k \cos(k\theta)]\]
Como lo de llevar a cuestas los senos y los cosenos por separado es bastante incómodo, vamos a unificarlos aunque sea algo más complejo (estoy sembrao), usando:\[\sin\theta = \frac{1}{2i}\left(e^{i\theta} - e^{-i\theta}\right) \qquad\qquad \cos\theta = \frac{1}{2}\left(e^{i\theta} + e^{-i\theta}\right)\]
Con lo cual el desarrllo de Fourier es: \[f(\theta) = \sum_{k=-\infty}^\infty c_k e^{ik\theta}\]
Como buen desarrollo en serie, esperaremos que para las funciones razonables, los coeficientes \(c_k\) vayan tendiendo a cero al aumentar \(|k|\), de modo que podamos truncar en algún momento la serie.

¿Y cómo calculamos los \(c_k\)? Pues vamos a valernos de la siguiente propiedad: la ortogonalidad de las exponenciales complejas. ¿Orto-qué? Muy fácil, tomando dos enteros distintos, \(m\) y \(n\), hagamos la siguiente integral:
\[\int_0^{2\pi} e^{im\theta}e^{-in\theta}d\theta = \int_0^{2\pi} e^{i(m-n)\theta}d\theta=\left.\frac{e^{i(m-n)\theta}}{i(m-n)}\right|_0^{2\pi} = \frac{1-1}{i(m-n)}=0\]
mientras que si tenemos \(m=n\):
\[\int_0^{2\pi}e^{im\theta}e^{-im\theta}d\theta = \int_0^{2\pi}d\theta = 2\pi\]
Aplicándola al desarrollo de Fourier tenemos:
\[\int_0^{2\pi}e^{-im\theta} f(\theta)d\theta = \int_0^{2\pi}e^{-im\theta}\sum_{k=-\infty}^\infty c_k e^{ik\theta} d\theta = \sum_{k=-\infty}^\infty c_k \int_0^{2\pi} e^{-im\theta}e^{ik\theta}d\theta = 2\pi c_m\]
con lo que concluimos que
\[c_m = \frac{1}{2\pi}\int_0^{2\pi} e^{-im\theta} f(\theta)d\theta\]

Comprimidos

El desarrollo de Fourier tiene innumerables aplicaciones. Podría hablaros, por ejemplo, de su relación con el principio de incertidumbre, o de su papel para determinar los gaps de energía en sólidos cristalinos. Pero vamos a ver un ejemplo más cotidiano, que seguro que todos habéis utilizado en alguna ocasión. Solamente tenéis que mirar al principio de la entrada. 
- ¿Una mujer desnuda?¿Puedo usar el desarrollo de Fourier para desnudar a una mujer? 
- No creo que sea una técnica para ligar muy efectiva. 
- Ceci n'est pas une femme
- Vale, un cuadro.  
- Mais non! Ceci n'est pas un tableau
- En fin, tú ganas, una foto de un cuadro de una mujer desnuda. 
Y más concretamente, una foto en formato JPEG.


Al efectuar el procesamiento de señales digitales, (como es el caso de las imágenes) en vez de una función \(f(\theta)\), se trabaja con puntos muestreados uniformemente. Con un número suficiente de muestras (determinado por el teorema de Nyquist) es posible reconstruir su serie de Fourier. Para ello se utiliza la transformada discreta de Fourier (DFT). Este algoritmo es tan fundamental que resulta importante hacerlo de la manera más eficiente posible: la transformada rápida de Fourier (FFT), para \(N\) muestras, necesita \(N \log N\) operaciones frente a las \(N^2\) de los algoritmos originales. Y si en vez de ordenadores convencionales usamos ordenadores cuánticos, el poder actuar simultáneamente sobre los estados cuánticos que están superpuestos en un qubit acelera aún más el proceso.

Una imagen de tipo mapa de bits es un conjunto de valores muestreados (la unidad de muestreo es el píxel) de tres magnitudes (rojo/verde/azul, tono/brillo/saturación).
El canal de brillo de la imagen de la maja

Para guardarla como jpg, la imagen es dividida en cuadrados de 8 píxeles de lado. Con esos 64 puntos, para cada una de las magnitudes del color se realiza su FFT (en realidad una transformada discreta de coseno, pero es básicamente lo mismo) y se obtienen 64 coeficientes del desarrollo de Fourier. 
La maja tras pasarse una FFT

Hasta aquí parece que no hemos ganado nada, de unos 64 números hemos conseguido otros 64. El truco está en que, como hemos visto, el valor de los coeficientes \(c_k\) va decreciendo a medida que aumenta \(|k|\): por ello, los últimos términos del desarrollo en cada bloque se pueden despreciar. 
 
Nos hemos cargado el 25% de los datos con menor \(|c_k|\)

Al no considerarlos, la imagen ocupa menos espacio, pero también se pierde calidad. 
 
Resultado tras deshacer la FFT

Al guardar una imagen como jpg, los programas te suelen pedir la calidad de la imagen, que indica cómo de pequeños deben ser los valores de \(c_k\) para que los cuente como ceros.
La maja con una calidad del 15% (izquierda) y del 90% (derecha)
No sé si el contenido del post es adecuado para el #lunesTetas, pero seguro que sí lo es para el #lunes\(\theta\)s (leído #lunesThetas)

Para saber más 

E. Roberts: Lossy algorithms - JPEG. University of Stanford
K. Cabeen y P. Gent: Image compression and the discrete cosine transform (pdf). College of the Redwoods
Las imágenes están realizadas con GIMP y su plugin para hacer FFT.