El impacto de una revelación inesperada

¿A poco no cambian las cosas cuando en el momento menos esperado descubres de una u otra manera algo que por algún tiempo había vagado sin rumbo aparente por los pasillos del inconsciente? Siempre sucede de la misma manera. Al estar buscando  respuestas a los grandes enigmas universales, uno se encuentra en el mundo de las tinieblas ya que la claridad nunca llega en el momento indicado y más aún, siempre surgen dudas adicionales en el camino. Pero, también todo es diferente cuando en un solo instante las cosas parecen aclararse, no me refiero a resolver los enigmas, sino que a raíz de dicha epifanía uno cree comprender el porqué de las cosas, tampoco significa que tienes que estar de acuerdo con ellas…, pero así sucede, desgraciadamente .

Si tenemos una imagen en formato PGM, que según sus autores “está designado para ser extremadamente fácil de aprender”, ya que cada uno de los valores de los datos representa un pixel en escala de grises de 8 bits, no representa en mayor problema hacer un programa para obtener la imagen del archivo. Afortunadamente no es este el caso, en MATLAB leemos un archivo pgm de la siguiente manera:

imagen = imread('peppers.pgm','pgm');

Ahora, imagen es una variable que contiene una matriz de MxN pixeles, donde N y M dependen del tamaño de la imagen. A esto lo vemos como nuestra señal, y le aplicamos la Transformada de Fourier Discreta, pero al ser una imagen y como lo mencionaba antes, es una matriz de 2 dimensiones, lo que quiere decir que tenemos que realizar la transformada a través de los renglones y después de las columnas. El código que se presenta a continuación fue sacado de aquí, aunque lo modifique ligeramente para solamente pasarle los datos y calcular automáticamente el tiempo y las frecuencias a partir de esto, lo que realiza esta función es calcular la DFT a los datos de entrada:

1. function X=dft(x)
2. [sx1 sx2] =size(x); %Obtener el tamaño de la matriz
3. t = [0:sx2-1];     %Calcular la base de tiempo
4. f = [0:(1/sx2):(1-(1/sx2))];    %Calcular las frecuencias
5. t = t(:); % Formatear 't' en vector columna
6. x = x(:); % Formatear 'x' en vector columna
7. f = f(:); % Formatear 'f' en vector columna
8. W = exp(-2*pi*j * f*t');
9. X = W * double(x);

Código 1. Función que realiza la DFT en MATLAB

Las últimas dos líneas de código representan en esencia la transformada de fourier. Tengo que hacer referencia al post que escribí iniciando este tema, ya que no recuerdo de que fuente tomé la ecuación en el post anterior, pero noto que es diferente en cuanto X(k)  y x(n) está representadas como X(k+1) y x(n+1), creo que la diferencia radica en donde se tomen los coeficientes.

Ok, otra vez, la ecuación siguiente nos muestra la transformada de fourier discreta (Tomada de: Practical Digital Signal Processing for Engineers and Technicians, pp 63)

Ecuación 1

Ahora, viendo el programa relacionamos la línea 8 con lo siguiente:

   Ecuación 2

Podemos observar de la ecuación anterior que k es el equivalente a la líneas 3 y  n/N es la línea 4. Esto es, k es el ‘tiempo’ o mas bien dicho cada muestra en la imagen  y las frecuencias en las que contienen son n/N. (Si, esto me sigue causando problemas emocionales todavía). En la última línea del código anterior (la nueve) tenemos las muestras x(n) con los componentes de frecuencia (lo representado en la ecuación 2).

1. imagen = imread('peppers.pgm','pgm');
2. imfinfo('peppers.pgm')
3. figure(1); imshow(imagen);
4. [s1 s2]= size(imagen);
5. dftimagen = zeros(s1);
6. for i=1:s2
7.     dftimagen(i,:) = dft(imagen(i,:));
8. end
9. dftimagen= dftimagen.';
10. for i=1:s2
11.     dftimagen(i,:) = dft(dftimagen(i,:));
12. end
13. dftimagen = dftimagen.';
14. figure(2);imshow(uint8(ifft2(dftimagen)));

Código 2. Programa que realiza la DFT a una imagen de 8 bits escala de grises en MATLAB

De las líneas 1 a la 3 se lee la imagen y muestra en pantalla, en las líneas 4 y 5 se crea una nueva variable para almacenar el resultado de la transformada (este código tiene el inconveniente que da por hecho que es una matriz cuadrada).  Las líneas 6 a 8 realizan la DFT a cada renglón de la imagen, posteriormente se realiza una traspuesta (en la línea 9) a la imagen para realizar la DFT a lo largo de las columnas en las lineas 10 a 12 y se regresa con la traspuesta de la línea 13.

Por último, se toma la variable que tiene almacenada la imagen transformada y usando la transformada inversa que viene incluida en MATLAB se recupera la imagen original y se muestra en pantalla (línea 14) para verificar que se halla realizado correctamente la DFT.

Solo para terminar esta imagen en la que hice pruebas es una imagen de 512 x 512 pixeles y es bastante lento, no tengo los datos de cuanto se tarda, pero no se compara con la función que trae integrada MATLAB para realizar la FFT (no de en vano se llama Fast Fourier Transform, claro)

Isn’t ironic? Don’t you think?