Notas de Visión

Posgrado en Ingeniería Eléctrica. Maestría y Doctorado en Ingeniería

Eléctrica. Especialidades en Sistemas

Computacionales.

Dr. Félix Calderón Solorio

Septiembre 2009

Índice general

1. Introducción 1

1.1. Creación de Imágenes . . . . . . . . . . . . . . . . . . . . . . . 1 1.1.1. Manejo de colores RGB . . . . . . . . . . . . . . . . . 1 1.1.2. La escala de gris . . . . . . . . . . . . . . . . . . . . . 2 1.1.3. Interfase para desplegar imágenes en Java. . . . . . . . 3

1.2. Operadores Sobre imágenes . . . . . . . . . . . . . . . . . . . 3 1.2.1. Operadores puntuales sobre Imágenes . . . . . . . . . . 3 1.2.2. Operadores de Ventana sobre imágenes . . . . . . . . . 6

1.3. Transformaciones Geométricas . . . . . . . . . . . . . . . . . . 7 1.3.1. Transformación Afín . . . . . . . . . . . . . . . . . . . 7 1.3.2. Traslación . . . . . . . . . . . . . . . . . . . . . . . . . 7 1.3.3. Interpolación. . . . . . . . . . . . . . . . . . . . . . . . 9

2. Calibración de Camáras 19 2.1. El modelo de cámara pinhole . . . . . . . . . . . . . . . . . . 19 2.2. Cámaras CCD . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 2.3. Propiedades de la cámara Proyectiva . . . . . . . . . . . . . . 22 2.4. Calibración de Cámaras a partir de puntos correspondientes . 22

2.4.1. Centro de la Cámara . . . . . . . . . . . . . . . . . . . 23 2.4.2. Descomposición RQ . . . . . . . . . . . . . . . . . . . . 23 2.4.3. Ejemplo . . . . . . . . . . . . . . . . . . . . . . . . . . 24 2.4.4. Calculo del modelo utilizando puntos 3D - 2D . . . . . 25 2.4.5. Ejemplo . . . . . . . . . . . . . . . . . . . . . . . . . . 25

3. Procesamiento de señales 27 3.1. Definiciones. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 3.2. Señales pares e impares . . . . . . . . . . . . . . . . . . . . . . 27

3.2.1. Ejemplos. . . . . . . . . . . . . . . . . . . . . . . . . . 28 3.3. Señales básicas de tiempo continuo. . . . . . . . . . . . . . . . 31

3.3.1. Exponencial compleja y senoidal. . . . . . . . . . . . . 31 3.3.2. Función escalón unitario e impulso unitario. . . . . . . 32

4. Sistemas 33

4.1. Introducción . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 4.1.1. Sistemas inversos. . . . . . . . . . . . . . . . . . . . . . 34 4.1.2. Sistemas causales. . . . . . . . . . . . . . . . . . . . . . 34 4.1.3. Estabilidad. . . . . . . . . . . . . . . . . . . . . . . . . 35

4.2. Invariancia en el tiempo. . . . . . . . . . . . . . . . . . . . . . 35

4.2.1. Linealidad . . . . . . . . . . . . . . . . . . . . . . . . . 36 4.3. Ejemplos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36

4.3.1. Ejemplo . . . . . . . . . . . . . . . . . . . . . . . . . . 36 4.3.2. Ejemplo . . . . . . . . . . . . . . . . . . . . . . . . . . 36 4.3.3. Ejemplo . . . . . . . . . . . . . . . . . . . . . . . . . . 37

5. Convolución 39

5.1. Correlación . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 5.1.1. Propiedades . . . . . . . . . . . . . . . . . . . . . . . . 40 5.1.2. Ejemplo . . . . . . . . . . . . . . . . . . . . . . . . . . 41 5.1.3. Correlación en dos dimensiones . . . . . . . . . . . . . 41

5.2. Representación de señales en términos de impulsos. . . . . . . 42 5.3. Convolución . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43

5.3.1. Propiedades . . . . . . . . . . . . . . . . . . . . . . . . 43 5.3.2. Sucesión útil . . . . . . . . . . . . . . . . . . . . . . . . 45 5.3.3. Ejemplos . . . . . . . . . . . . . . . . . . . . . . . . . . 47

5.4. Algunos kerneles interesantes . . . . . . . . . . . . . . . . . . 49 5.4.1. Suavizadores . . . . . . . . . . . . . . . . . . . . . . . 49 5.4.2. Derivadas . . . . . . . . . . . . . . . . . . . . . . . . . 50

5.5. Respuesta de Sistemas lineales invariantes en el tiempo a exponenciales complejas. . . . . . . . . . . . . . . . . . . . . . . 54

6. Transformada de Fourier 56

6.1. Representación de señales periódicas. . . . . . . . . . . . . . . 56

6.1.1. Ejemplo . . . . . . . . . . . . . . . . . . . . . . . . . . 58 6.2. Propiedades de la Transformada de Discreta Fourier . . . . . . 61

6.2.1. Periodicidad . . . . . . . . . . . . . . . . . . . . . . . . 61 6.2.2. Linealidad . . . . . . . . . . . . . . . . . . . . . . . . . 61 6.2.3. Desplazamiento en tiempo . . . . . . . . . . . . . . . . 62 6.2.4. Desplazamiento en frecuencia . . . . . . . . . . . . . . 63 6.2.5. Conjugación . . . . . . . . . . . . . . . . . . . . . . . . 63 6.2.6. Inversión en Tiempo . . . . . . . . . . . . . . . . . . . 64 6.2.7. Escalamiento en tiempo . . . . . . . . . . . . . . . . . 64 6.2.8. Convolución . . . . . . . . . . . . . . . . . . . . . . . . 65 6.2.9. Multiplicación . . . . . . . . . . . . . . . . . . . . . . . 66 6.2.10. Diferenciación en Tiempo . . . . . . . . . . . . . . . . 66 6.2.11. Diferenciación en Frecuencia . . . . . . . . . . . . . . . 67 6.2.12. Propiedades de Simetría de la transformada de Fourier. 67 6.2.13. Relación de Parseval. . . . . . . . . . . . . . . . . . . . 69 6.2.14. Resumen de Propiedades . . . . . . . . . . . . . . . . . 69

6.3. Transformada de Fourier de algunas funciones interesantes . . 70 6.3.1. Exponencial Compleja . . . . . . . . . . . . . . . . . . 70 6.3.2. Función seno . . . . . . . . . . . . . . . . . . . . . . . 71 6.3.3. Función coseno . . . . . . . . . . . . . . . . . . . . . . 71 6.3.4. Función impulso unitario . . . . . . . . . . . . . . . . . 71 6.3.5. Constante . . . . . . . . . . . . . . . . . . . . . . . . . 72 6.3.6. Escalón Unitario . . . . . . . . . . . . . . . . . . . . . 72 6.3.7. Caja . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73 6.3.8. Ejemplos . . . . . . . . . . . . . . . . . . . . . . . . . . 74

6.4. Transformada de Fourier en dos dimensiones. . . . . . . . . . . 77 6.5. Convolución utilizando TF. . . . . . . . . . . . . . . . . . . . 81 6.6. Teorema del Muestreo . . . . . . . . . . . . . . . . . . . . . . 82

6.6.1. Transformada de Fourier de una señal periódica …. 82 6.6.2. Transformada de Furier de un Tren de Pulsos . . . . . 83 6.6.3. Señal muestreada . . . . . . . . . . . . . . . . . . . . . 84 6.6.4. Integración de la señal continua. . . . . . . . . . . . . . 86 6.6.5. Ejemplos . . . . . . . . . . . . . . . . . . . . . . . . . . 87

7. Filtros. 90 7.1. Filtros Pasa bajas. . . . . . . . . . . . . . . . . . . . . . . . . 90 7.2. Filtros Pasa Altas. . . . . . . . . . . . . . . . . . . . . . . . . 92 7.3. Filtro pasa bajas Butterworth. . . . . . . . . . . . . . . . . . . 93 7.4. Filtros de pasa banda. . . . . . . . . . . . . . . . . . . . . . . 94

7.4.1. Ejemplo . . . . . . . . . . . . . . . . . . . . . . . . . . 95 7.5. Filtro de Gabor . . . . . . . . . . . . . . . . . . . . . . . . . . 97

7.5.1. Filtro de Gabor en dos dimensiones . . . . . . . . . . . 97 7.6. Como entonar un filtro . . . . . . . . . . . . . . . . . . . . . . 98

7.6.1. Ejemplo . . . . . . . . . . . . . . . . . . . . . . . . . . 100 7.7. Filtro Elimina Banda . . . . . . . . . . . . . . . . . . . . . . . 100 7.8. Filtro de Mediana . . . . . . . . . . . . . . . . . . . . . . . . . 101 7.9. Filtro Binario . . . . . . . . . . . . . . . . . . . . . . . . . . . 101 7.10. Filtro de Membrana . . . . . . . . . . . . . . . . . . . . . . . 102

7.10.1. Probabilidad de un evento . . . . . . . . . . . . . . . . 102 7.10.2. Probabilidad condicional . . . . . . . . . . . . . . . . . 102 7.10.3. Independencia . . . . . . . . . . . . . . . . . . . . . . . 104 7.10.4. Regla de Bayes . . . . . . . . . . . . . . . . . . . . . . 104 7.10.5. Estimador de maxima verosimilitud. . . . . . . . . . . 104 7.10.6. Interpretación del filtro de Membrana en el dominio de

la Frecuencia. . . . . . . . . . . . . . . . . . . . . . . . 106

8. Aplicaciones 108

8.1. Detección de corrimientos utilizando técnicas de restauración . 108 8.1.1. Ejemplos . . . . . . . . . . . . . . . . . . . . . . . . . . 109

9. Tareas 113 9.1. Operadores de Ventana . . . . . . . . . . . . . . . . . . . . . . 113 9.2. Transformaciones Afín y Proyectivas . . . . . . . . . . . . . . 114 9.3. Cálculo de una Homografía . . . . . . . . . . . . . . . . . . . . 116 9.4. Modelo de la Cámara Pin-Hole . . . . . . . . . . . . . . . . . 119 9.5. Modelo Proyectivo de la cámara a partir de Imágenes de un

modelo 3D . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121 9.6. Convolución . . . . . . . . . . . . . . . . . . . . . . . . . . . . 127 9.7. Algunos Kerneles Interesantes . . . . . . . . . . . . . . . . . . 131 9.8. Cálculo de la Transformada Discreta de Fourier . . . . . . . . 136 9.9. Propiedades de la Transformada Discreta de Fourier . . . . . . 139 9.10. Transformada Discreta de Fourier de una señal muestreada . 148 9.11. Filtros . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 150

Capítulo 1

Introducción

1.1. Creación de Imágenes

El Puntillismo aparece por vez primera en 1884, encabezado por el pintor Neo - Impresionista Georges Seurat, y contando con entre sus más fieles seguidores tales como Henri - Edmond Cross y Vlaho Bokovac.Su procedimiento empleado por estos artistas, consistente en poner puntos de colores puros en vez de pinceladas sobre tela; logrará una técnica que fue el resultado de los estudios cromáticos llevados a cabo por Georges Seurat (1859 - 1891), pintor francés, quien en 1884 llegó a la división de tonos por la posición de toques de colores que, mirados a cierta distancia, crean en la retina las combinaciones deseadas [Wikipedia].

Figura 1.1: Imagen de puntos

1.1.1. Manejo de colores RGB

Definimos un pixel como la unidad básica de una imagen. Es un punto en cuyo valor se almacena la información conserniente a los colores. Los colores básicos que manejaremos son Rojo, Verde y Azul. Así pixel lo podemos ver como un arreglo P = [P_{r}, P_{g}, P_{b}]

En la siguiente figura podemos ver el como se almacena la información conserniente a los 32 bit o 4 Bytes

Figura 1.2: Formato de 32 bits para manejo de color

Una imagen discreta, la podemos definir como un arreglo bidimensional en el cual se almacena la información de cada punto de la imagen 256^{3} = 16, 777, 216diferentes colores por pixel. Dado que la información se maneja en Bytes para almacenar un byte necesitaremos de 4 bytes, tres correspondientes a los colores Rojo, verde y Azul y un cuarto que tiene información del nivel de transparencia. En la cámaras es común hablar de su resolución en Megabytes, ¿Que relación existe con esto?. Bueno una camara de 3 Mega pixeles, tiene aproximadamente 314 pixeles de lado. Esto es \sqrt{3 ∗ 1000 ∗ 1000} = 1732,05

En Forma matricial podemos hacer una imagen con colores verdes si definimos, un arreglo bidimensional donde ponemos un vector de colores {R,G,B}

( \begin{matrix}( \begin{matrix}0 \\ 255 \\ 0 \end{matrix} ) & ( \begin{matrix}0 \\ 0 \\ 0 \end{matrix} ) \\ ( \begin{matrix}0 \\ 0 \\ 0 \end{matrix} ) & ( \begin{matrix}0 \\ 255 \\ 0 \end{matrix} ) \end{matrix} )

Para este ejemplo si nuestro arreglo de imagen se llama A los colores almacenados en el renglón i columna j, la notación para referirnos a este punto es A_{i,j}

1.1.2. La escala de gris

En el caso de una imágen en tono de gris solamente se maneja una sola matriz de datos, así el color rojo será igual que el verde y el azul. La unión de los colores dara como resultado blanco y la ausencia dará negro. Nuestra imagen solamente tendrá colores que van desde el negro al blanco en 256 tonos diferentes.

Figura 1.3: Cubo de colores y escala de gris

)

Utilizando imágenes en tono de gris, solamente será necesario almacenar la información conserniente a un canal de color, o un promedio de los tres canales de color. Esto reduce en 1/3 la información que debe ser almacenada.

1.1.3. Interfase para desplegar imágenes en Java.

Ver Algoritmos y ejemplos en la página. ver inicial.java

1.2. Operadores Sobre imágenes

1.2.1. Operadores puntuales sobre Imágenes

Una imagen discreta, la podemos definir como un arreglo bidimensional en el cual se almacena la información de cada punto de la imagen. Como definimos en la sección anterior cada unidad es un pixel y en el caso de operadores puntuales consideraremos que cada uno de estos se aplica sobre cada pixel de la imagen.

Linea Recta

Dada la imagen que se muestra en la figura 1.4(a) y su histograma 1.4(b) en tono de gris, podemos calcular una transformación lineal que nos permita cambiar el contraste de la imagen.

(a) Imagen original

(b) Histograma

Figura 1.4: Imagen original de Lena y su histograma en tono de gris

Consideremos el caso de la ecuación de una línea recta y = mx + b donde m es la pendiente y b el cruze por cero.

g_{i,j} = (f_{i,j} − f_{min}) ∗ (g_{max} − g_{min})/(f_{max} − f_{min}) + g_{min};

donde el rango dinámico de la imagen dada es f_{min} = 0 y f_{max} = 255 y el rango de la nueva imagen es g_{min} y g_{max} los cuales fijaremos en el rango [10, 50]. La ecuación resultante es

g_{i,j} = 0,4f_{i,j} + 50
∀ < i, j >∈ Ω

Así para una imagen dadas estas condiciones podemos cambiar la apariencia haciendo una proyección del espacio original de la imagen a un nuevo espacio dado por una línea recta. El resultado de la aplicación a cada punto se presenta en la figura 1.2.1. En la figura 1.5(b) podemos ver como la aplicación de este operador reduce el rango dinámico de la imagen reduciendo el rango del histograma en la misma medida.

En el caso de una imagen con mayor información esta lucirá como:

Cambio de color a tono de Gris

Otro ejemplo de operadores puntuales sobre imagenes es la conversión de colores RGB a escala de gris. Así para convertir la imagen a colores a tono

(a) Imagen o (b) Histograma

Figura 1.5: Imagen de Lena modificada y su histograma en tono de gris

de gris, simplemente sacamos el promedio de cada uno de los pixeles en la 3 bandas aplicando la operación

A_{i,j} = (R_{i,j}+ G_{i,j} + B_{i,j})/ 3

Con esta operación tendremos una sola matriz con el promedio de los colores y cuando se crea la imagen a color las matrice R, G, B serán iguales a la matriz A.

Función Sigmoide

Otra manera de hacer esto es considerar una función diferente como la sigmoide, definida como

y = \frac{1}{1 + e^{−βx}}

donde β es un parámetro que controla la pendiente del escalon. Para hacer la aplicación de la función sigmoide, será necesario hacer un cambio al rango [−1, 1] de la imagen original ya que el dominio y rango de la función sigmoide es [−∞, ∞] y [0, 1] respectivamente.

Como luce una imagen en estas condiciones. En la figura 1.2.1, podemos ver la aplicación de la función sigmoide con β = 100 a la imagen original en la figura 1.4(a). En esta figura podemos ver la, la imagen procesada, el histogramas de la imagen y la gráfica de la sigmoide aplicada.

(a) Imagen o (b) Histograma (c) Sigmoide con β = 100

Figura 1.6: Imagen de Lena modificada por la función sigmoide

1.2.2. Operadores de Ventana sobre imágenes

El objetivo de un operador de ventana es reemplazar el valor de un pixel, por el contenido pesado en una ventana al rededor de esta. Así por ejemplo

I_{i,j} = ∑_{k=i−nk}^{i+nk} ∑_{l=j−nk}^{j+nk} w_{k,l}I_{k,l.}

donde w_{k,l} es un peso que puede tener varios valores con lo cual tenemos una posibilidad casi infinita de operadores. En la figura 1.7 tenemos un operador correspondiente a la media, donde cada uno de los pixeles es reemplazado por el promedio en una ventana de tamaña (2n_{k} + 1) × (2n_{k} + 1) y con w_{k,l} = \frac{1}{(2n_{k}+1)^{2}} .

Figura 1.7: Lena aplicando un operador de ventana que simula la media

1.3. Transformaciones Geométricas

En esta sección básicamente analizaremos dos transformaciones geométricas en imágenes. Estan son la transformación Afín y la transformación proyectiva.

1.3.1. Transformación Afín

Una transformación afin es una transformación geométrica que esta constitida por translación, escalamiento, rotación y cizallamiento.

1.3.2. Traslación

Una translación la podemos hacer simplemente asuminedo que las nuevas coordenadas x̂ = x + tx ŷ = y + ty les sumamos un valor tx o ty según corresponda. En coordenadas homogeneas queda como

( x̂ ; ŷ ; 1 ) = ( \begin{matrix}1 & 0 & tx \\ 0 & 1 & ty \\ 0 & 0 & 1 \end{matrix} ) ( \begin{matrix}x \\ y \\ 1 \end{matrix} )

Escalamiento

El escalamiento pude entenderse como hacer que una figura geométrica cambie su tamaño. Un escalamiento en x lo podemos representar por como x̂ = x ∗ s_{x} y en y como ŷ = y ∗ s_{y}. En coordenadas homogeneas se pude expresar como

S = ( \begin{matrix}s_{x} & 0 & 0 \\ 0 & s_{y} & 0 \\ 0 & 0 & 1 \end{matrix} )

Rotación

Consideremos el caso de un punto que rota respecto a un punto fijo. Las coordenadas x y y, en forma polar las podemos obtener como x = rcos(θ) y y = rseno(θ). Si consideramos que esta gira un ángulo ψ entonces podemos representar esta rotación en forma polar como

( x̂ ; ŷ ) = ( \begin{matrix}r cos(θ + ψ) \\ r seno(θ + ψ) \end{matrix} ) = ( \begin{matrix}r cos(θ) cos(ψ) − rseno(θ)seno(ψ) \\ r cos(θ)seno(ψ) + rseno(θ) cos(ψ) \end{matrix} )

( x̂ ; ŷ ) = ( \begin{matrix}x cos(ψ) − yseno(ψ) \\ xseno(ψ) + y cos(ψ) \end{matrix} )

Una transformación de rotación puede ser definida como

( x̂ ; ŷ ) = ( \begin{matrix}cos(ψ) & −seno(ψ) \\ seno(ψ) & cos(ψ) \end{matrix} ) ( \begin{matrix}x \\ y \end{matrix} )

y en coordenadas homogeneas puede representarse como

R(θ) = ( \begin{matrix}cos(θ) & −sin(θ) & 0 \\ sin(θ) & cos(θ) & 0 \\ 0 & 0 & 1 \end{matrix} )

Cizallamiento

El cizallamiento es una transformación dada por la matriz, donde c_{x} es el ángulo de cizallamiento respecto al eje x.

C_{x} = ( \begin{matrix}1 & tg(c_{x}^{∘}) & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{matrix} )

En caso de querer aplicar el cizallamiento en la dirección y, la matriz de transformación será.

C_{y} = ( \begin{matrix}1 & 0 & 0 \\ tg(c_{y}^{∘}) & 1 & 0 \\ 0 & 0 & 1 \end{matrix} )

Una vez que unimos todas las transformaciones tenemos una composición dada por una Rotacion R, un escalamiento S, un Cizallamiento C y una Translación T. Al conjunto de estas operaciones se le conoce como transformación afin.

T_{afin} = ( \begin{matrix}s_{x} cos(θ) & −s_{y} sin(θ) + s_{x} cos(θ) tan(c_{x}) & t_{x} \\ s_{x} sin(θ) & s_{y} cos(θ) + s_{x} sin(θ) tan(c_{x}) & t_{y} \\ 0 & 0 & 1 \end{matrix} )

1.3.3. Interpolación.

Todas las transformaciónes aplicadas anteriormente son para un punto en el espacio de coordenadas [x,y], pero que sucede cuando hablamos de imágenes. Una imagen esta compuesta de colores en una posición con coordenadas enteras [i,j]. Si aplicamos una transformación afín las coordenas enteras pasaran a ser coordenadas en números reales, así por ejemplo en el caso de una rotación tendremos

En estas condiciones en probable descidir hacer un truncamiento, sin embargo cuando se trata de un escalamiento con sx=2 y sy = 1, tendremos que ciertos pixeles no serán visitados.

Esta condición hará que no se tengán cubiertas todos los pixeles de una imagen ya que en el caso del escalamiento (s_{x} = 2) solamente las posiciones pares de la nueva imagen tendrán valor.¿Como corregir este efecto?.En lugar de barrer todos los puntos de la imagen origen, barreremos todos los puntos de la imagen destino y cualquier pixel que no este en posición entera lo calcularemos utilizando interpolación bilineal. Consideremos el caso de una matriz de rotación

T = ( \begin{matrix}0,951057 & −0,309017 & 0. \\ 0,309017 & 0,951057 & 0. \\ 0. & 0. & 1. \end{matrix} )

al aplicar la transformación para cada uno de los puntos de la imagen de acuerdo con el siguiente código tenemos, una solución dada por la figura 1.8

    for(i=0; i<nr; i++)
      for(j=0; j<nc; j++) {
        afin(i, j, x, T);
        if(x[0] >=0 && x[0] < nr && x[1] >=0 && x[1] < nc)
         I2[(int) x[1]][(int) x[0]] = I1[i][j];

}

Note los puntos negros generados por el redondeo. Esta condición es mas notoria cuando se escala una imagen con un factor superior a 1.

La función de interpolación bilineal en Java es

Figura 1.8: Lena aplicando una transformación

static public double bilineal(double f[][], double x0, double y0)

{

int ii, jj, nr, nc;
double a, b, y;
ii = (int) y0;
jj = (int) x0;
a = y0 - ii;
b = x0 - jj;
nr = f.length;
nc = f[0].length;
if (ii < 0 || jj < 0 || ii > nr - 2 || jj > nc - 2)
 return 0;
if ( (a == 0) && (b == 0))
 return f[ii][jj];
if (ii == nr - 1 && jj == nc - 1)
 return f[ii][jj];
if (b == 0 || jj == nc - 1) {
 y = (1 - a) * f[ii][jj] + a * f[ii + 1][jj];
 return y;

}

if (a == 0 || ii == nr - 1) {
 y = (1 - b) * f[ii][jj] + b * f[ii][jj + 1];
 return y;

}

y = (1 - a) * (1 - b) * f[ii][jj] + (1 - a) * b * f[ii][jj + 1]
+
  a * (1 - b) * f[ii + 1][jj] + a * b * f[ii + 1][jj + 1];
return y;

}

Una manera de resolver este problema es aplicar la transformación sobre los puntos de la imagen destino utilizando la inversa de la función y la función de interpolación de acuerdo con el código

 static public void afin(double f[][], double g[][], double T[])
c {
int i, j, nr, nc;
double in, jn;
nr = f.length;
nc = f[0].length;

for (i = 0; i < nr; i++)
 for (j = 0; j < nc; j++) { //System.out.println(T[0]);
  jn = T[0] * j  + T[1] * i  + T[2];
  in = T[3] * j  + T[4] * i  + T[5];
  g[i][j] = bilineal(f, jn, in);

} }

La diferencia la podemos ver en la imagen

(a) Aplicando T inversa con interpolación

(b) Aplicando T

Figura 1.9: Transformación proyectiva

ver (afin.java)

Transformación Log Polar

La geometría log polar en imagens fue motivada por la estrcutura de la retina de algunos sistema de vision y esta tiene algunas propiedades de compresion. Cuando comparamos el sistema usual de coordenadas cartesianas en imagenes, las imagenes log - polar permiten muestreo mas rápido en los sistema de visión artificial sin reducción en el tamaño de la vista y la parte centra de la retina (fovea).

La transformación log-polar es un mapeo de puntos en el plano cartesiano (x, y) a puntos en el plano log-polar (ξ, η):

ξ = log (\sqrt{x^{2} + y^{2}})

η = tg^{−1} (\frac{y}{x})

Aplicando esto tenemos que una imagen luce como la figura 1.10 ver (logpolar.java)

static public void LogPolar (double datos[][], double res[][]) {
  int nr, nc, i, j, ren, col;

Figura 1.10: Lena aplicando una transformación log polar

double alfa, r, x, y;
nr = datos.length;
nc = datos[0].length;
ren = res.length;
col = res[0].length;
System.out.println (ren + + col);
double rmax = Math.sqrt (ren*ren + col*col)/2.0;
for (i = 0; i < ren; i++)
 for (j = 0; j < col; j++) {
  alfa = (2.0*Math.PI*i)/(ren - 2);
  r = Math.exp (j*Math.log (rmax)/col);
  x = r*Math.cos (alfa) + nc/2.0;
  y = r*Math.sin (alfa) + nr/2.0;
  res[i][j] = bilineal (datos, x, y);

} } La tranformación log-polar inversa se calcula con el código

static public void InvLogPolar (double datos[][], double[][] res)

{

int nr, nc, i, j, ii, jj, ren, col;
double alfa, r;
double fac = 65.0;
nr = datos.length;
nc = datos[0].length;
ren = res.length;

 col = res[0].length;
 double phimax = Math.log (Math.sqrt (nr*nr + nc*nc)/2.0);
 for (i = 0; i < ren; i++)
  for (j = 0; j < col; j++) {
   ii = ren/2 - i;
   jj = col/2 - j;
   alfa = (Math.atan2 (ii, jj)/Math.PI + 1.0)*(double) (nr - 2)/2.0;
   r = Math.log (Math.sqrt (ii*ii + jj*jj))/phimax*nc;
   res[i][j] = bilineal (datos, r, alfa);

} }

Transformación Proyectiva

Una transformación proyectiva es la generalización de las transformaciones R^{2} → R^{2}, en la que las líneas paralelas no son transformadas necesariamene en tales. Se puede expresar de la forma

( x̂ ; ŷ ; ẑ ) = ( \begin{matrix}A & t \\ v^{T} & 1 \end{matrix} ) ( \begin{matrix}x \\ y \\ 1 \end{matrix} )

donde v^{T} es la diferencia entre la transformación proyectiva y la afín. Note para el caso de la transformación proyectiva la coordenada ẑ no necesariamente es unitaria lo cual implica que la imagen transfromada se fue a un plano con coordenadas diferentes a esta. Una manera de regresar al sistema de coordenadas es dividir entre ẑ. Así la matriz de transformación queda

( x̃ ; ỹ ; 1 ) = \frac{1}{ẑ} ( \begin{matrix}A & t \\ v^{T} & 1 \end{matrix} ) ( \begin{matrix}x \\ y \\ 1 \end{matrix} )

Esta ecuación es conocida tambien como homógrafía. La aplicación a una imagen de una transformación

T = ( \begin{matrix}1 & 0 & 0 \\ 0 & 1 & 0 \\ 1e − 3 & 0 & 1 \end{matrix} )

Aplicando esta transformación tenemos una proyeccion de la imagen de Lena que se muestra en la figura 1.11

Figura 1.11: Transformación Proyectiva en la imagen de Lena

El código Java para hacer esta transformación es:

static public void proyecccion (double f[][], double g[][], double
T[])

{

int i, j, nr, nc;
double in, jn, kn;
nr = f.length;
nc = f[0].length;
for (i = 0; i < nr; i++)
for (j = 0; j < nc; j++) {
// System.out.println (T[0]);
jn = T[0]*j + T[1]*i + T[2];
in = T[3]*j + T[4]*i + T[5];
kn = T[6]*j + T[7]*i + 1.0;
jn /= kn;
in /= kn;
g[i][j] = bilineal (f, jn, in);

} }

Parámetros a partir de puntos correspondientes

En esta sección vamos a calcular los parametros de una homografía, la cual representa un modelo proyectivo, dada por la ecuación

( \begin{matrix}x_{2} \\ y_{2} \\ z_{2} \end{matrix} ) = ( \begin{matrix}h_{1} & h_{2} & h_{3} \\ h_{4} & h_{5} & h_{6} \\ h_{7} & h_{8} & h_{9} \end{matrix} ) ( \begin{matrix}x_{1} \\ y_{1} \\ z_{1} \end{matrix} ) = ( \begin{matrix}A & t \\ v^{T} & k \end{matrix} ) ( \begin{matrix}x_{1} \\ y_{1} \\ z_{1} \end{matrix} )

donde el vector P_{1} = [x_{1}, y_{1}, z_{1}]^{T} es un punto en el espacio tridimensional, el vector P_{2} = [x_{2}, y_{2}, z_{2}]^{T} es la proyección en un plazo con z_{2} = 1y ambos puntos son correspondientes P_{1} ⟷ P_{2}

El problema consiste en determinar el vector H = [h_{1}, h_{2}, h_{3}, h_{4}, h_{5}, h_{6}, h_{7}, h_{8}]^{T} a partir de un conjunto de puntos correspondientes. Una manera de resolver este problema es utilizando mínimos cuadrados, sin embargo tiene problema de estabilidad numérica por lo que es preferible utilizar el método siguiente.

Dado los puntos y que estos son correspondientes, entonces podemos decir que P_{2} × HP_{1} = 0, dado que ambos vectores on paralelos. Esto no lleva a la siguiente representación

( \begin{matrix}0^{T} & −z_{2,i}P_{1,i}^{T} & y_{2,i}P_{1,i}^{T} \\ z_{2,i}P_{1,i}^{T} & 0^{T} & −x_{2,i}P_{1,i}^{T} \end{matrix} ) ( \begin{matrix}H_{1} \\ H_{2} \\ H_{3} \end{matrix} ) = 0

donde H_{1} = [h_{1}, h_{2}, h_{3}]^{T}, H_{2} = [h_{4}, h_{5}, h_{6}]^{T} y H_{3} = [h_{7}, h_{8}, h_{9}]^{T} .

La solución del sistema de ecuaciones puede llevarse a cabo utilizando la descompisición en valores singulares SVD. La descomposición en valores singulares calcula tres matrices USV = SVD(A) donde U y V son matrices y S una matriz diagonal. La solución del sistema es el último vector de la matriz V (ver Harley-Zisserman pp 90).

Ejemplo

Dados los puntos P1 = { [10.0 10.0 1.0], [10.0 50.0 1.0], [50.0 50.0 1.0 ], [50.0 10.0 1.0 ], [25.0 25.0 1.0]} y P2 = {[8.9108, 10.89108, 1.0], [4.95049, 46.53465, 1.0], [39.04761, 48.57142, 1.0], [42.85714, 14.28571, 1.0 ], [20.48780, 25.36585, 1.0 ], determinar la Homografía correspondiente. La siguiente imagen muestra los puntos P1 a la izquierda y los puntos P2 a la derecha.

Figura 1.12: Puntos correspondientes

De acuerdo con la ecuación

( \begin{matrix}0^{T} & −z_{2,i}P_{1,i}^{T} & y_{2,i}P_{1,i}^{T} \\ z_{2,i}P_{1,i}^{T} & 0^{T} & −x_{2,i}P_{1,i}^{T} \end{matrix} ) ( \begin{matrix}H_{1} \\ H_{2} \\ H_{3} \end{matrix} ) = BH = 0

Tenemos

B = ( \begin{matrix}0. & 0. & 0. & −10. & −10. & −1. & 108,911 & 108,911 & 10,891 \\ 10. & 10. & 1. & 0. & 0. & 0. & −89,109 & −89,109 & −8,911 \\ 0. & 0. & 0. & −10. & −50. & −1. & 465,347 & 2326,73 & 46,535 \\ 10. & 50. & 1. & 0. & 0. & 0. & −49,505 & −247,525 & −4,95 \\ 0. & 0. & 0. & −50. & −50. & −1. & 2428,57 & 2428,57 & 48,571 \\ 50. & 50. & 1. & 0. & 0. & 0. & −1952,38 & −1952,38 & −39,048 \\ 0. & 0. & 0. & −50. & −10. & −1. & 714,286 & 142,857 & 14,286 \\ 50. & 10. & 1. & 0. & 0. & 0. & −2142,86 & −428,571 & −42,857 \\ 0. & 0. & 0. & −25. & −25. & −1. & 634,146 & 634,146 & 25,366 \\ 25. & 25. & 1. & 0. & 0. & 0. & −512,195 & −512,195 & −20,488 \end{matrix} )

Haciendo la SVD(B) tenemos que la matriz v es

v = ( \begin{matrix}−0,008799 & −0,01804 & 0,4915 & 0,0774 & −0,4276 & −0,5010 & −0,07493 & 0,3712 & −0,4178 \\ −0,006549 & −0,001930 & 0,6014 & −0,7039 & 0,06790 & 0,3676 & −0,009472 & −0,02489 & 0,04643 \\ −0,0001982 & −0,0003293 & 0,01553 & −0,01201 & −0,01818 & 0,007560 & 0,7942 & −0,3909 & −0,4641 \\ −0,008647 & −0,002579 & −0,5090 & −0,7008 & −0,03908 & −0,4922 & 0,005096 & 0,05909 & −0,04643 \\ −0,01059 & 0,0176 & −0,3700 & −0,07726 & −0,4068 & 0,6088 & −0,06827 & 0,3757 & −0,41782 \\ −0,0002489 & 0,0002550 & −0,01247 & −0,01494 & −0,02256 & −0,005127 & −0,5945 & −0,6558 & −0,4642 \\ 0,7018 & 0,7119 & 0,01084 & −0,00820 & −0,01566 & −0,01418 & 0,00008070 & 0,0008968 & −0,0004646 \\ 0,7119 & −0,7018 & −0,01073 & −0,007749 & −0,01538 & 0,01458 & −0,0002675 & 0,0004953 & 2,6997e − 7 \\ 0,01740 & −0,000201 & −0,001353 & 0,02651 & 0,8025 & −0,01339 & −0,07222 & 0,3660 & −0,4642 \end{matrix} )

Tomando la última columna de la matriz v y normalizando respecto a v_{9,9} tenemos que la Homografía es

H = ( \begin{matrix}0,900036 \\ −0,100016 \\ 0,999866 \\ 0,100021 \\ 0,899986 \\ 0,99995 \\ 0,00100087 \\ −5,8153e − 7 \\ 1. \end{matrix} )

Calibración de Camáras

En esta sección daremos el modelo de una cámara y la manera de calcular el módelo de la misma. El modelo esta diseñando principalmente para sensores como CCD pero también es aplicable a otras cámaras.

2.1. El modelo de cámara pinhole

Consideremos el punto central de proyección del espacio en un plano. Centraremos la proyección del espacio Euclidiano de coordenadas y consideraremos que el plano es Z = f, el cual es llamado el plano focal de la cámara. Bajo el modelo de camara pinhole, un punto de coordenadas en el espacio X = [X, Y, Z]^{T} es mapeado en el punto del plano de la imagen donde las lineas se unen en el centro de proyeción. Esto puede mostrarse por triángulos similares en las figuras 2.1(a) y 2.1(b).

En la figura, por triangulos similares podemos ver que el punto [X, Y, Z]^{T} es mapeado a un punto en el plano de la cámara como [fX/Z, fY/Z, f]^{T}. Esto lo podems escribir de forma matricial como

( \begin{matrix}fX \\ fY \\ Z \end{matrix} ) = ( \begin{matrix}f & & & 0 \\ & f & & 0 \\ & & 1 & 0 \end{matrix} )( \begin{matrix}X \\ Y \\ Z \\ 1 \end{matrix} )

esta matriz puede ser escrita como diag(f, f, 1)[I|0] así la matriz de proyección para una cámara con modelo pinhole es P = diag(f, f, 1)[I|0]. La

(a) Tres dimensiones  (b) dos dimensiones

Figura 2.1: Modelo de la cámara de pinhole

expresión anterior asume que el origen de coordenadas en el plano de la imagen es el punto principal. En la práctica no es así por lo cual debemos aplicar un mapeo

(X, Y, Z)^{T} | → (fX/Z + p_{x}, fX/Z + p_{y})^{T}

esta ecuación puede expresarse como

( \begin{matrix}fX + Zp_{x} \\ fY + Zp_{y} \\ Z \end{matrix} ) = ( \begin{matrix}f & & p_{x} & 0 \\ & f & p_{y} & 0 \\ & & 1 & 0 \end{matrix} )( \begin{matrix}X \\ Y \\ Z \\ 1 \end{matrix} )

donde

K = ( \begin{matrix}f & & p_{x} \\ & f & p_{y} \\ & & 1 \end{matrix} )

y x = K[I|0]X_{cam} donde la matriz K es llamada la matriz de calibración de la cámara.

En general puntos en el espacio son expresados en términos de coordenadas euclidianas diferentes, conocidas como marco de coordenadas del mundo. Estos dos marco son relacionados por una rotación y una translación. Si X̃son las coordenadas homogeneas del mundo estas pueden ser representadas en el marco de la camara como como X̃_{cam} = R (X̃ − C̃), donde C̃ representa las coordenadas de la cámara centradas en las coordenadas del mundo. El modelo puede ser escrito en coordenadas homogenas como

x_{cam} = ( \begin{matrix}R & −RC̃ \\ 0 & 1 \end{matrix} ) ( \begin{matrix}X \\ Y \\ Z \\ 1 \end{matrix} )

Poniendo todo junto tenemos x = KR [I |−C̃ ] X donde la matriz de la

cámara es P = K[R|t] donde t = −RC̃

2.2. Cámaras CCD

El modelo de cámara pinhole asume que las coordenadas de la imagen tienen igual escala en ambos ejes. En el caso de la cámaras CCD existe una posibilidad adicional de no tener pixeles cuadrados. Por lo tanto manejaremos en lugar de usar un factor de escala f en la diagonal consideraremos que existe un número de pixeles en la direción x llamado m_{x} y uno en la direción y llamado m_{y} por lo tanto nuestra matriz de parámetros queda como

K = ( \begin{matrix}fm_{x} & & x_{0} \\ & fm_{y} & y_{0} \\ & & 1 \end{matrix} ) = ( \begin{matrix}α_{x} & & x_{0} \\ & α_{y} & y_{0} \\ & & 1 \end{matrix} )

Podemos consideran un modelo mas general si agregamos un parametro s al que denominaremos sesgo. En la mayoría de la cámaras normales este valor es cero, sin embargo lo consideraremos en el sentido de tener un modelo más general.

K = ( \begin{matrix}α_{x} & s & x_{0} \\ & α_{y} & y_{0} \\ & & 1 \end{matrix} )

El modelo de la cámara proyectiva es una matriz no cuadrada con tres renglones y cuatro columnas, por lo tanto el modelo tendrá 11 grados de libertad. Este modelo lo podemos escribir como P = M [I |M^{−1}p_{4} ] = KR [I |−C̃ ]donde p_{4} es la última columna de la matriz P. El modelo completo es conocido como el módelo de la cámara proyetiva.

2.3. Propiedades de la cámara Proyectiva

Algunas propiedades de la cámara proyectiva son

Centro de la camara. El centro de la cámara es el vector C el cual puede ser calculado para una cámara finita como el vector en el espacio nulo es decir PC=0. Esto se resuelve utilizando descomposición en Valores singulares y

Puntos columna. Para i =1,2,3 los vectores columna p_{i}son los puntos de fuga en la imagen correpondientes a los ejes X, Y y Z. Columna p_{4} es la imagen de las coordenadas del origen.

Plano principal. El plano principal de la cámara is p^{3}, el último renglón de P.

Ejes del plano. Los planos p^{1}y p^{2} representan los planos a través del centro de la cámara.

Punto principal. Punto de imagen x_{0} = Mm^{3}es el punto principal de la cámar, donde m^{3T} es el tercer renglón de la matriz M.

Rayo principal. El rayo principal de la camara es el rayo que pasa a través de la cámara con centro C con dirección m^{3T}. El eje principal es el vector v = det(M)m^{3}es la dirección hacia el frente de la cámara.

2.4. Calibración de Cámaras a partir de puntos correspondientes

La calibración de cámaras la definiremos como el proceso de calcular los parámetros internos a partir de la observación de una imagen y sus puntos correspondiente en el mundo real. Comenzaremos por dada la matriz de proyección, calcular la matriz de parametros internos K, el centro de la cámara C y la matriz de rotación de la cámara P = KR[I|C].

2.4.1. Centro de la Cámara

Podemos ver de lo expuesto en la sección anterior que el centro de la cámara dado la matriz de proyección se puede calcular como P C = 0 cuya solución es el vector con el menor eigenvalor una vez calculado la descomposición en valores singulares (espacio nulo).

2.4.2. Descomposición RQ

Definamos la matriz M = KR la cual es una matrix de 3 renglones por 3 columnas. Para ello tendremos que hacer una factorización matricial. Recordemos que la matriz R puede calcularse como la multiplicación de tres matrices de rotación en cada eje como R = R_{z}^{T} R_{y}^{T} R_{x}^{T} dado que todas estas matrices son ortogonales, entonces podemos escibir

M = KR = KR_{z}^{T} R_{y}^{T} R_{x}^{T} MR_{x}R_{y}R_{z} = K

Las matrices de rotación estan definidas como

R_{x} = ( \begin{matrix}1 & 0 & 0 \\ 0 & c_{x} & −s_{x} \\ 0 & s_{x} & c_{x} \end{matrix} ) , R_{y} = ( \begin{matrix}c_{y} & 0 & s_{y} \\ 0 & 1 & 0 \\ −s_{y} & 0 & c_{y} \end{matrix} ) y R_{x} = ( \begin{matrix}c_{y} & −s_{z} & 0 \\ s_{z} & 1 & 0 \\ 0 & 0 & c_{y} \end{matrix} )

donde c_{i} = cos (θ_{i}) y s_{i} = seno (θ_{i})

Dada la matriz M Los pasos son

1. Multiplicar M′ = MR_{x} y calcular c_{x} y s_{x}para que el elemento M′_{3,2} = 0 y c_{x}^{2} + s_{x}^{2} = 1. La solución es

c_{x} = \frac{M_{3,3}}{\sqrt{M_{3,2}^{2}+M_{3,3}^{2}}} ys_{x} = −\frac{M_{3,2}}{\sqrt{M_{3,2}^{2}+M_{3,3}^{2}}}

2. Multiplicar M” = M′R_{y} y calcular c_{y} y s_{y}para que el elemento M”_{3,1} = 0 y c_{y}^{2} + s_{y}^{2} = 1. La solución es

c_{y} = \frac{M_{3,3}}{\sqrt{M_{3,1}^{2}+M_{3,3}^{2}}} ys_{y} = −\frac{M_{3,1}}{\sqrt{M_{3,1}^{2}+M_{3,3}^{2}}}

3. Finalmente multiplicar K = M”R_{z} y calcular c_{z} y s_{z}para que el elemento M”_{2,1} = 0 y c_{z}^{2} + s_{z}^{2} = 1. La solución es

c_{z} = −\frac{M_{2,2}}{\sqrt{M_{2,1}^{2}+M_{2,2}^{2}}} ys_{z} = \frac{M_{2,1}}{\sqrt{M_{2,1}^{2}+M_{2,2}^{2}}}

ver calibra camaras.java

2.4.3. Ejemplo

Dada la matriz de proyección de la cámara

P = ( \begin{matrix}353,553 & 339,645 & 277,744 & −1,44946 × 10^{6} \\ −103,528 & 23,3212 & 459,607 & −632525. \\ 0,707107 & −0,353553 & 0,612372 & −918,559 \end{matrix} )

determinar el centro de la cámara y la matriz de parámetros internos. La matriz M de acuerdo con la definición es:

M = ( \begin{matrix}353,553 & −103,528 & 0,707107 \\ 339,645 & 23,3212 & −0,353553 \\ 277,744 & 459,607 & 0,612372 \end{matrix} )

Para encontrar el centro de la cámara hacemos la descomposición en Valores singulares de la matriz [u, s, v] = svd(P) donde

u = ( \begin{matrix}−0,916531 & −0,399962 & −0,000424686 \\ −0,399962 & 0,916532 & −0,000479022 \\ −0,000580829 & 0,00026918 & 1. \end{matrix} )

s = ( \begin{matrix}1,58146 × 10^{6} & 0. & 0. & 0. \\ 0. & 406,369 & 0. & 0. \\ 0. & 0. & 0,837954 & 0. \end{matrix} )

v = ( \begin{matrix}−0,000178718 & −0,581477 & 0,723847 & 0,371391 \\ −0,000202738 & −0,281691 & −0,607393 & 0,742782 \\ −0,000277203 & 0,76324 & 0,327293 & 0,557086 \\ 1. & 0,0000505429 & 0,0000969494 & 0,00037139 \end{matrix} )

El vector v_{4} es el que tiene el menor eigenvalue así que dividiendo este entre v_{4,4} tenemos que el centro de la camára es C = {1000., 2000., 1500., 1.}

De la descomposición QR tenemos que la matriz M puede representarse como:

K×R = ( \begin{matrix}468,1648 & 91,2251 & 300,0001 \\ 0,0 & 427,201 & 199,9999 \\ 0,0 & 0,0 & 1,0 \end{matrix} )×( \begin{matrix}0,4138 & 0,9091 & 0,0471 \\ −0,5733 & 0,2201 & 0,7892 \\ 0,7071 & −0,3535 & 0,6124 \end{matrix} )

2.4.4. Calculo del modelo utilizando puntos 3D - 2D

Así dado un conjunto de puntos en 3D Q = {Q_{1}, Q_{2},..., Q_{N}} con Q_{i} = (X_{i}, Y_{i}, Z_{i})^{T} y otro conjunto de puntos correspondientes en el plano de la imagen q = {q_{1}, q_{2},..., q_{N}} con q_{i} = (x_{i}, y_{i}, w_{i})^{T} podemos determinar el modelo de la cámara con al menos 6 puntos correspondientes X_{i} ⟷ x_{i}resolviendo

( \begin{matrix}x_{i} \\ y_{i} \\ w_{i} \end{matrix} ) = ( \begin{matrix}P_{1,1} & P_{1,2} & P_{1,3} \\ P_{2,1} & P_{2,2} & P_{2,3} \\ P_{3,1} & P_{3,2} & P_{3,3} \end{matrix} ) ( \begin{matrix}X_{i} \\ Y_{i} \\ Z_{i} \end{matrix} ) = ( \begin{matrix}Q_{i}^{T} P_{1} \\ Q_{i}^{T} P_{2} \\ Q_{i}^{T} P_{3} \end{matrix} )

Para resolver calculamos un vector paralelo, el cual tiene un producto vectorial igual a cero q_{i} × P Q_{i} = 0. El sistema resultante lo podemos calcular encontrando el vector en el espacio nulo con descomposición en valores singulares. El sistema resultante a resolver es

q_{i} × P Q_{i} = ( \begin{matrix}î & ĵ & k̂ \\ x_{i} & y_{i} & w_{i} \\ Q_{i}^{T} P_{1} & Q_{i}^{T} P_{2} & Q_{i}^{T} P_{3} \end{matrix} ) = ( \begin{matrix}y_{i}Q_{i}^{T} P_{3} − w_{i}Q_{i}^{T} P_{2} \\ x_{i}Q_{i}^{T} P_{3} − w_{i}Q_{i}^{T} P_{1} \\ x_{i}Q_{i}^{T} P_{3} − w_{i}Q_{i}^{T} P_{1} \end{matrix} ) = 0

de lo cual podemos desprender

( \begin{matrix}0 & −w_{i}Q_{i}^{T} & y_{i}Q_{i}^{T} \\ −w_{i}Q_{i}^{T} & 0 & x_{i}Q_{i}^{T} \end{matrix} ) ( \begin{matrix}P_{1} \\ P_{2} \\ P_{3} \end{matrix} ) = AP = 0

2.4.5. Ejemplo

Dados los los puntos correpondientes Q = { [0.0, 0.0, 0.0], [ 100.0, 0.0, 0.0], [100.0, 100.0, 0.0], [0.0, 100.0 , 0.0], [0.0, 0.0, 100.0], [100.0, 0.0 100.0], [100.0, 100.0, 100.0], [0.0, 100.0, 100.0 ]} y q = {[1577.9715837523775, 688.6057400776651], [1667.8746657863205, 758.2462570249891], [1562.6523714350803, 725.2525691697816], [1483.8812040033367, 660.6388854847863], [1658.2870049496003, 684.1821822330893], [1762.4087684498734, 758.8465253032915], [1645.2811209801278, 723.3689600937458], [1554.5611061379307, 654.4720145727945]}, calcular

a) La matriz de proyección de la cámara,

b) La posición del centro de la cámara,

c) La matriz de Rotación y

d) la matriz K de parámetros internos.

A_{1.,6} ( \begin{matrix}0. & 0. & 0. & 0. & 0. & 0. \\ 0. & 0. & 0. & 1. & 0. & 0. \\ 0. & 0. & 0. & 0. & −99,9999 & 0. \\ 100. & 0. & 0. & 1. & 0. & 0. \\ 0. & 0. & 0. & 0. & −99,9999 & −99,9999 \\ 100. & 100. & 0. & 1. & 0. & 0. \\ 0. & 0. & 0. & 0. & 0. & −99,9999. \\ 0. & 100. & 0. & 1. & 0. & 0. \\ 0. & 0. & 0. & 0. & 0. & 0. \\ 0. & 0. & 100. & 1. & 0. & 0. \\ 0. & 0. & 0. & 0. & −99,9999 & 0. \\ 0. & 0. & 0. & 0. & 0. & 0. \\ 100. & 100. & 100. & 1. & 0. & 0. \\ 0. & 0. & 0. & 0. & 0. & 0. \end{matrix} )

A_{7.,12} ( \begin{matrix}0. & −0,9999 & 0. & 0. & 0. & 688,606 \\ 0. & 0. & 0. & 0. & 0. & −1577,97 \\ 0. & −0,9999 & 75824,6 & 0. & 0. & 758,246 \\ 0. & 0. & −166787. & 0. & 0. & −1667,87 \\ 0. & −0,9999 & 72525,3 & 72525,3 & 0. & 725,253 \\ 0. & 0. & −156265. & −156265. & 0. & −1562,65 \\ 0. & −0,9999 & 0. & 66063,9 & 0. & 660,639 \\ 0. & 0. & 0. & −148388. & 0. & −1483,88 \\ −99,9999 & −0,9999 & 0. & 0. & 68418,2 & 684,182 \\ 0. & 0. & 0. & 0. & −165829. & −1658,29 \\ −99,9999 & −0,9999 & 75884,7 & 0. & 75884,7 & 758,846 \\ 0. & 0. & −1,27487 × 10^{10} & 0. & −1,27487 × 10^{10} & −1,27487 × 10^{6} \\ 0. & 0. & −164528. & −164528. & −164528. & −1645,28 \\ 0. & 0. & −1,01742 × 10^{10} & −1,01742 × 10^{10} & −1,01742 × 10^{6} \end{matrix} )

a) Comenzamos por hacer la SVD. El último vector de la matriz v después de la descomposición es V_{12} = [2,2356E−4, 2,1476E−4, 1,7562E−4, −0,9165, −6,5463E−5, 1,4746E− 5, 2,9062E−4, −0,3999, 4,4712E−7, −2,2356E−7, 3,8721E−7, −5,8082E−4] La matriz de proyección

P = [ \begin{matrix}2,2356E − 4 & 2,1476E − 4 & 1,7562E − 4 & −0,9165 \\ −6,5463E − 5 & 1,4746E − 5 & 2,9062E − 4 & −0,3999 \\ 4,4712E − 7 & −2,2356E − 7 & 3,8721E − 7 & −5,8082E − 4 \end{matrix} ]

b) El centro de la cámara es Para calcular el centro de la cámara calculamos svd de P. La matriz v es

P = [ \begin{matrix}−1,0E − 4 & 0,5815 & −0,7237 & 0,3714 \\ −1,0E − 4 & 0,2817 & 0,6074 & 0,7428 \\ −2,0E − 4 & −0,7631 & −0,3272 & 0,5571 \\ 1,0 & 0,0 & 0,0 & 4,0E − 4 \end{matrix} ]

El vector con el eigenvalor menor es v_{4} = [0,3714, 0,7128, 0,5571, 4,0E−4] por lo tanto el centro de la cámara es C = [1000,0007, 2000,002, 1500,0003, 1,0]

De la descomposición QR tenemos que la matriz M puede representarse como:

K×R = ( \begin{matrix}468,1648 & 91,2251 & 300,0001 \\ 0,0 & 427,201 & 199,9999 \\ 0,0 & 0,0 & 1,0 \end{matrix} )×( \begin{matrix}0,4138 & 0,9091 & 0,0471 \\ −0,5733 & 0,2201 & 0,7892 \\ 0,7071 & −0,3535 & 0,6124 \end{matrix} )

Procesamiento de señales

3.1. Definiciones.

Representaremos a una señal discreta por x(t) la cual pude ser función de una o varias variables. Esta señal discreta resulta de hacer la multiplicación de una señal continua por un tren de pulsos representados por

x (t) = f (t) δ(t − kT)

donde f (t) es una función continua y δ(t − kT) es un tren de pulsos unitarios con periodo T. La función δ(t) es la delta de dirac la cual toma el valor unitario cuando su argumento es cero (ver figura 3.1).

3.2. Señales pares e impares

Podemos decir que una señal es par si es idéntica a su reflexién alrededor del origen, esto es

x(−t) = x(t)

y que es impar si

x(−t) = −x(t)

Una característica importante de cualquier señal, es que esta pude ser representada por la suma de una señal par y una señal impar

x (t) = P(x(t)) + I(x(t))

Figura 3.1: Multiplicación de una función por un tren de pulsos donde la parte par la calculamos

P(x(t)) = \frac{1}{2} [x (t) + x(−t)]

y la impar por

I(x(t)) = \frac{1}{2} [x (t) − x(−t)]

ver [Oppenheim, 1997]

3.2.1. Ejemplos.

(a) Demostrar que si x (n) es discreta impar entonces

∑_{n=−∞}^{∞} x (n) = 0

Demostración: si, x (n) es impar, entonces

x (n) = −x (−n)

es decir

x (n) + x (−n) = 0

considerando

∑_{n=−∞}^{∞} x (n) = ∑_{n=−∞}^{−1} x (n) + x (0) + ∑_{n=1}^{∞} x (n)

reordenando

∑_{−∞}^{∞} x (n) = x (0) + ∑_{n=−∞}^{−1} x (n) + ∑_{n=1}^{∞} x (n)

= x (0) + ∑_{n=1}^{∞} x (−n) + ∑_{n=1}^{∞} x (n)

= x (0) − ∑_{n=1}^{∞} x (n) + x (−n)

con x (0) = 0 y sustituyendo la definición de señal par en la ecuación anterior tenemos :

∑_{−∞}^{∞} x (n) = 0 + ∑_{n=1}^{∞} 0

∑_{−∞}^{∞} x (n) = 0

(b) Comprobar, si x_{1} (n) es impar y x_{2} (n) es par, entonces

x_{1} (n) · x_{2} (n) ⟹ impar

Comprobación: Hacemos que

z (n) = x_{1} (n) · x_{2} (n)

y

z (−n) = x_{1} (−n) · x_{2} (−n)

Además, sabemos que

x_{1} (n) = −x_{1} (−n)

y

x_{2} (n) = x_{2} (−n)

sustituyendo las ecuaciones tenemos

z (n) = −x_{1} (−n) · x_{2} (−n)

resulta

z (n) = −z (−n)

Lo cual implica que z (n) es impar, i.e. x_{1} (n) · x_{2} (n) es impar.

(c) Considere que x (n) es una señal con parte par P (x (n)) y parte impar I (x (n)) , demostrar que

∑_{n=−∞}^{∞} x^{2} (n) = ∑_{n=−∞}^{∞} P^{2} (x (n)) + ∑_{n=−∞}^{∞} I^{2} (x (n))

Demostración: Hacemos

P (x (n)) = x_{p} (n)
I (x (n)) = x_{i} (n)

y

∑_{n=−∞}^{∞} x^{2} (n) = ∑_{n=−∞}^{∞} [x_{p} (n) + x_{i} (n)]^{2}

= ∑_{n=−∞}^{∞} (x_{p}^{2} (n) + 2x_{p} (n) x_{i} (n) + x_{i}^{2} (n))

sustituyendo el término x_{p} (n) x_{i} (n) por z (n) y reordenando

∑_{n=−∞}^{∞} x^{2} (n) = ∑_{n=−∞}^{∞} (x_{p}^{2} (n) + x_{i}^{2} (n)) + 2 ∑_{n=−∞}^{∞} z (n)

Como z (n) es el producto de una señal par y una impar, usamos el resultado del inciso b) y, concluimos que z (n) es impar. Ahora, usando el resultado del inciso a) sabemos que

2 ∑_{n=−∞}^{∞} z (n) = 0

Si sustituimos este resultado en la ecuación tenemos que

∑_{n=−∞}^{∞} x^{2} (n) = ∑_{n=−∞}^{∞} (x_{p}^{2} (n) + x_{i}^{2} (n))

∑_{n=−∞}^{∞} x^{2} (n) = ∑_{n=−∞}^{∞} x_{p}^{2} (n) + ∑_{n=−∞}^{∞} x_{i}^{2} (n)

Una señal es periódica si tiene la propiedad de tener un valor positivo T para el cual

x (t) = x (t + T) ∀t

en este caso diremos que la señal es periódica con periodo T

3.3. Señales básicas de tiempo continuo.

3.3.1. Exponencial compleja y senoidal.

La señal exponencial de tiempo continuo es de la forma

x (t) = Ce^{at}

donde C y a son, en general números complejos. Dependiendo de los valores de estos parámetros, la exponencial puede adoptar varias características diferentes. Una clase de exponenciales complejas que nos interesa es

x (t) = e^{jw_{o}t}

donde w_{0} = 2πf_{0} es definida como la velocidad angular en radianes por segundo y f_{0} es la frecuencia en ciclos por segundo o Hertz. Utilizando la relación de Euler esta señal puede expresarse como

x (t) = cos (w_{0}t) + jsen (w_{o}t)

Una propiedad importante de esta señal, es su periodicidad. Una señal es periódica si cada intervalo de tiempo T, tenemos un repetición de la señal. Para verificar esta propiedad hacemos

e^{jw_{o}t} = e^{jw_{o}(t+T)}
e^{jw_{o}t} = e^{jw_{o}t}e^{jw_{o}T}

para que esta ecuación se cumpla debemos tener que

e^{w_{o}T} = 1

Existen dos posibilidades para que esta condición se cumpla: la primera cuando w_{o} = 0, la cual es periódica para cualquier valor de T, pero si w_{o} ≠ 0, entonces tenemos que existe un valor T_{o} al cual llamamos periodo fundamental y esta dado por la ecuación

T_{o} = \frac{2π}{|w_{0}|}

3.3.2. Función escalón unitario e impulso unitario.

Otra señal de interés es la función escalón unitario la cual esta dada por

u (t) = \{ \begin{matrix}0, t < 0 \\ 1, t > 0 \end{matrix}

Otra señal bastante ótil es la función impulso unitario de tiempo continuo la cual estó relacionada con el escalón unitario por la ecuación

u (t) = ∫_{−∞}^{t} δ (τ) dτ

Sistemas

4.1. Introducción

Un sistema se puede ver como cualquier proceso que produce una transformación de señales. Un sistema tiene una señal de entrada y una señal de salida la cual estó relacionada con la entrada a través de la transformación del sistema.

y(t) = T [x(t)]

donde x(t) es la señal de entrada, y (t) es la señal de salida y T [] es la transformación del sistema.

Entre los sistemas podemos tener interconexiones serie y paralelo dadas como los mostrados en la figura 4.1.

Sistemas con y sin memoria.

Si la salida de un sistema para cada valor de la variable independiente depende solo de la entrada en ese mismo instante de tiempo se dice que el sistema no tiene memoria.

y(t) = Rx(t)

Un ejemplo de un sistema con memoria es

y(n) = ∑_{k=−∞}^{n} x(k)

Figura 4.1: Sistemas serie y paralelo

4.1.1. Sistemas inversos.

Decimos que un sistema es invertible si dada una transformación T podemos encontrar la transformación T^{−1} tal que:

y(t) = T [x(t)]
x(t) = T^{−1}[y(t)]

Un ejemplo de un sistema que no es invertible es y (t) = 0.

4.1.2. Sistemas causales.

Un sistema es causal si su salida en cualquier instante de tiempo depende sólo de los valores en el tiempo presente y en el pasado. Estos sistemas también son llamados no anticipativo.

y(t) = x(t − 1)
y(t) = ∫_{−∞}^{t} x(τ)dτ

4.1.3. Estabilidad.

Un sistema es llamado estable si se produce una salida acotada para una entrada acotada. Una señal es acotada si

|x(k)| ≤ M < ∞

y la señal de salida y(k) = T [x(k)] , es acotada si

|y(k)| ≤ N < ∞

Ejemplo. Consideremos la secuencia de los nómeros de Fibonacci donde x(k) = {1, 2, 3, 5, 8, 13, ...} ó x(k) = x(k −1)+x(k −2). En este caso tenemos un sistema no acotado ya que la sucesión crece indefinidamente para valores de k.

Un ejemplo de señal estable lo tenemos en la sucesión

x (n) = ∑_{k=1}^{n} \frac{1}{n^{2}}

para esta sucesión podemos ver que

x (n) − x (n − 1) = \frac{1}{n^{2}}

si tomamos el límite podemos verificar que la sucesión converge en un valor estable

lím_{n→∞} [x (n) − x (n − 1)] = lím_{n→∞} \frac{1}{n^{2}} = 0

4.2. Invariancia en el tiempo.

Para que un sistema sea invariante en el tiempo se debe cumplir que para un desplazamiento en la señal de entrada se produzca el mismo desplazamiento en la señal de salida.

T [D^{n} [x (t)]] = D^{n} [T [x (t)]]

Ejemplo. Considere la señal y (t) = sen [x (t)]

D^{k} [x (t)] = x(t − k)
= sen [x(t − k)]
D^{k}y (t) = sen [x(t − k)]

4.2.1. Linealidad

La característica principal de los sistemas lineales es

T [αx(t) + βy(t)] = T [αx(t)] + T [βy(t)]

esta propiedad es conocida como el principio de superposición.

Ejemplo. Considere el sistema y(t) = mx(t) + b. Que valores debe tener el sistema para ser lineal.

Consideremos dos señales x_{1}(t) y x_{2}(t)

y_{1}(k) = mx_{1}(t) + b
y_{2}(k) = mx_{2}(t) + b
y_{1}(k) + y_{2}(k) = m (x_{1}(t) + x_{2}(t)) + 2b

Si aplicamos la transformación a la suma de x_{1}(t) y x_{2}(t) tendremos

y_{3}(k) = m (x_{1}(t) + x_{2}(t)) + b

note que y_{3}(k) es diferente de y_{1}(k) + y_{2}(k), la ónica posibilidad es que la constante b sea igual a cero.

4.3. Ejemplos

4.3.1. Ejemplo

Demuestre que y(n) = x(n) − x(n − 1) es invariante en el tiempo. Solución. Sea, v(n) = D^{n_{0}} {T [x(n)]} y v̂(n) = T {D^{n_{0}}[x(n)]}. Un sistema es invariante si v(n) = v̂(n)

Para nuestro ejemplo v(n) = D^{n_{0}}[x(n)−x(n−1)] = x(n−n_{0})−x(n−n_{0}−1) por otro lado v̂(n) = T [x(n−n_{0})] = x(n−n_{0})−x(n−n_{0}−1), como v(n) = v̂(n) el sistema es invariante en el tiempo.

Nota: En este ejemplo, hacemos una aproximación de la derivada utilizando diferencias finitas y(n) ≃ \frac{dx(n)}{dn} = x(n) − x(n − 1) .

4.3.2. Ejemplo

Demuestre que y(n) = ∑_{k=n−n_{0}}^{n} x(k) es invariante en el tiempo.

Solución. En óste caso, v(t) = D^{n_{0}} [∑_{k=k_{0}}^{k_{f}} x(k)] = ∑_{k=k_{0}}^{k_{f}} x(k − n_{0}) y

v̂(k) = T [D^{n_{0}} [x(k)]] = ∑_{k=k_{0}}^{k_{f}} x(k − n_{0}). Como v(k) = v̂(k) el sistema es invariante en el tiempo.

4.3.3. Ejemplo

Dada la sucesión

y(n) = \frac{1}{4}y(n − 1) + x(n)  (4.1)

1. Probar que es un sistema Lineal invariante en el Tiempo (LIT) y

2. Determinar la salida y(n) si x(n) = δ(n − 1)

Comenzaremos por dar una solución de la recurrencia, haciendo sustituciones sucesivas

y(n) = x(n) + \frac{1}{4}y(n − 1)

= x(n) + \frac{x(n − 1)}{4} + \frac{1}{4}y(n − 2)

= x(n) + \frac{x(n − 1)}{4} + \frac{x(n − 2)}{4^{2}} + \frac{1}{4^{2}} y(n − 3)

= x(n) + \frac{x(n − 1)}{4} + \frac{x(n − 2)}{4^{2}} + \frac{x(n − 3)}{4^{3}} + ⋯ + \frac{x(1)}{4^{n−1}} + \frac{x(0)}{4^{n}}

donde de manera general la solución de la recursión es

y(n) = ∑_{k=0}^{n} \frac{x(n − k)}{4^{k}}  (4.2)

Linealidad Por otra parte tenemos que verificar que el sistema dado por 4.1 es Lineal,

primero calculamos utilizando 4.2

T [αx_{1}(n)] + T [βx_{2}(n)] = α ∑_{k=0}^{n} \frac{x_{1}(n − k)}{4^{k}} + β ∑_{k=0}^{n} \frac{x_{2}(n − k)}{4^{k}}

y segundo utilizando 4.2 calculamos

T [αx_{1}(n) + βx_{2}(n)] = ∑_{k=0}^{n} \frac{αx_{1}(n − k) + βx_{2}(n − k)}{4^{k}}

= α ∑_{k=0}^{n} \frac{x_{1}(n − k)}{4^{k}} + β ∑_{k=0}^{n} \frac{x_{2}(n − k)}{4^{k}}

lo cual demuestra que el sistema es lineal.

Invariancia en el tiempo

En este caso tenemos que mostrar D^{n_{0}} [T [x(n)]] = T [D^{n_{0}} [x(n)]] , aplicando la solución de la recursión dada por 4.2, tenemos

D^{n_{0}} [T [x(n)]] = D^{n_{0}} [∑_{k=0}^{n} \frac{x(n − k)}{4^{k}}] = ∑_{k=0}^{n} \frac{x(n − n_{0} − k)}{4^{k}}

Si aplicamos la traslación a la señal x tenemos

D^{n_{0}} [x(n)] = x(n − n_{0}), x(n − 1 − n_{0}), x(n − 2 − n_{0}), ⋯, x(−n_{0}),

y posteriormente calculamos

T [D^{n_{0}} [x(n)]] = ∑_{k=0}^{n} \frac{x(n − n_{0} − k)}{4^{k}}

De lo anterior podemos concluir que el sistema es invariante en el tiempo. Respuesta al impulso Para determinar la salida y(n) = \frac{1}{4}y(n − 1) + x(n) y suponiendo que

y(−1) = 0, se obtiene la siguiente tabla al evaluar recursivamente:

\begin{matrix}n & x(n) & y(n) \\ 0 & 1 & 1 \\ 1 & 0 & \frac{1}{4} \\ 2 & 0 & \frac{1}{4^{2}} \\ 3 & 0 & \frac{1}{4^{3}} \end{matrix}

Se observa que para cualquier n la salida seró y(n) = (\frac{1}{4})^{n}, note que esta solución puede verificarse utilizando la solución de la recursión.

Convolución

5.1. Correlación

Consideremos dos señales x y y, donde x = [y_{i}, y_{i+1}, ⋯ y_{i+N}] es un subconjunto de los valores de la señal y = [y_{0}, y_{1}, y_{2}, ⋯, y_{M}]. óEn que punto la señoal x empatará mejor con la señal y. Comenzaremos por utilizar como medida de similitud el producto escalar de vectores, así el ángulo entre estas dos señales lo podemos calcular como

cos(θ) = \frac{x^{T} ŷ}{|x||ŷ|}

donde ŷ = [y_{j}, y_{j+1} ⋯, y_{j+N}] y θ es el ángulo entre los dos vectores. Note que cuando los dos vectores son iguales tendremos que el ángulo es cero, el valor máximo será dado por |x|^{2} y que el producto x^{T} ŷ tenderá a este máximo.

El producto de x^{T} ŷ lo podemos calcular como

r_{yx} = ∑_{k=0}^{N} x(k)ŷ(k)

dado que ŷ es un subconjunto de y, si aplicamos un desplazamiento podemos calcular el producto x^{T} ŷ_{n} en las n posiciones posibles de la serie con la siguiente expresión

r_{xy}(n) = ∑_{k=0}^{N} x(k)y(k + n)

La cross–correlación de las señales x (n) y y (n) es la secuencia r_{xy} (n) definida por

r_{xy} (n) = ∑_{k=−∞}^{∞} x (k) y (k + n)

r_{xy} (t) = ∫_{−∞}^{∞} x (τ) y (τ + t) dτ

La función de correlación suministra una medida de la similitud o inter– dependencia entre las funciones x y y.

5.1.1. Propiedades

Conmutación

Demostrar que la función de correlación y auto–correlación no son conmutativas, es decir

r_{xy} (t) = r_{yx} (−t)
r_{xx} (t) = r_{xx} (−t)

para ello hacemos

r_{xy} = ∫_{−∞}^{∞} x (τ) y (τ + t) dτ

y sustituimos τ̂ = τ + t

r_{yx}(−t) = ∫_{−∞}^{∞} x (τ̂ − t) y \binom{}{τ} dτ̂

en el caso de la auto–correlación tenemos

r_{xx} = ∫_{−∞}^{∞} x (τ) x (τ + t) dτ

y sustituimos τ̂ = τ + t

r_{xx}(−t) = ∫_{−∞}^{∞} x (τ̂ − t) x \binom{}{τ} dτ̂

Lo cual significa que la correlación no es conmutativa.

Invariancia a traslación

Es la correlación invariantes a traslación.

r_{xy} = ∫_{−∞}^{∞} x (τ + a) y (τ + t) dτ

si sustituimos

τ̂ = τ + a
τ = τ̂ − a

∫_{−∞}^{∞} x \binom{}{τ} y (τ̂ + (t − a)) dτ = r_{xy} (t − a)

dado lo anterior la correlación no es invariante a traslación.

Linealidad

Pero, Será lineal la correlación?

r_{xy} (t) = ∫_{−∞}^{∞} [x_{1} (τ) + x_{2} (τ)] y (τ + t) dτ

= ∫_{−∞}^{∞} x_{1} (τ) y (τ + t) dτ + ∫_{−∞}^{∞} x_{2} (τ) y (τ + t) dτ
= r_{x_{1}y} + r_{x_{2}}y

si es una transformación lineal la correlación.

5.1.2. Ejemplo

Suponga una señal x(t) = [0, 1, 2, 3, 4, 0, 1, 2, 3, 4, …] (Fig. 6.5(a)) y una señal y(t) = [0, 1, 2, 3, 4] (Fig. 5.1(b)). La correlación entre ambas señales, se muestra en la Fig. 6.5(b). Note que los máximos se obtienen cuando la señal y empata perfectamente con x.

5.1.3. Correlación en dos dimensiones

Si la señal esta en dos dimensiones, como es el caso de imágenes la correlación se representa como

r_{xy}(n, m) = ∑_{k=−n_{k}}^{n_{k}} ∑_{l=−n_{l}}^{n_{l}} x(k, l) ∗ y(k + n, l + m)

(a) x(t)  (b) y(t)

(c) z(t) correlación entre x(t) y y(t)

Figura 5.1: Correlación

5.2. Representación de señales en términos de impulsos.

La función impulso unitario, puede utilizarse para construir una clase amplia de señales. Para ilustrar como funciona consideremos que tenemos una señal discreta dada por x(n). Si queremos ver un parte de la señal en cierto instante de tiempo multiplicamos la señal de entrada por un impulso en el instante que deseamos analizar, así por ejemplo:

x (−1) δ(n + 1) = \{ \begin{matrix}x (−1) si n=-1 \\ 0 si no \end{matrix}

x (0) δ(n + 0) = \{ \begin{matrix}x (0) si n=0 \\ 0 si no \end{matrix}

x (1) δ(n − 1) = \{ \begin{matrix}x (1) si n=1 \\ 0 si no \end{matrix}

por lo tanto la suma suma de estos términos me da la señal x. Matemáticamente podemos expresar esta suma como

y (n) = ∑_{k=−∞}^{∞} x (k) δ (n − k)

esta ecuación es llamada la propiedad de escudriñamiento del impulso unitario.

5.3. Convolución

Podemos extender el concepto sustituyendo la función impulso por cualquier otra función, así obtenemos la expresión de la convolución.

y (n) = ∑_{k=−∞}^{∞} x (k) h (n − k)

note que también es un sistema LIT. La convolución la podemos representar por y = x ∗ h

5.3.1. Propiedades

Conmutación

La convolución es conmutativa es decir

f_{1} ∗ f_{2} = f_{2} ∗ f_{1}

demostración, para probar la propiedad conmutativa hacemos

f_{1} ∗ f_{2} = ∑_{k=−∞}^{∞} f_{1} (k) f_{2} (n − k)

haciendo el cambio de variable k̂ = n − k tenemos

∑_{k=−∞}^{∞} f_{1} (n − k̂) f_{2} \binom{}{k} = f_{2} ∗ f_{1}

Asociativa

La propiedad asociativa consiste en

(f_{1} ∗ f_{2}) ∗ f_{3} = f_{1} ∗ (f_{2} ∗ f_{3})

para demostrar la propiedad asociativa hacemos f_{1} ∗ f_{2} = g y f_{2} ∗ f_{3} = h

g ∗ f_{3} = f_{1} ∗ h

puesto que

g (n) = ∑_{k=−∞}^{∞} f_{1} (k) f_{2} (n − k)

se tiene que

g (n) ∗ f_{3} (n) = ∑_{j=−∞}^{∞} g (j) f_{3} (n − j)

= ∑_{j=−∞}^{∞} [ ∑_{k=−∞}^{∞} f_{1} (k) f_{2} (j − k)] f_{3} (n − j)

sustituyendo l = j − k y cambiando el orden de las sumatorias

= ∑_{k=−∞}^{∞} f_{1} (k) ∑_{l=−∞}^{∞} f_{2} (l) f_{3} (n − (l + k))

= ∑_{k=−∞}^{∞} f_{1} (k) ∑_{l=−∞}^{∞} f_{2} (l) f_{3} ((n − k) − l)

= ∑_{k=−∞}^{∞} f_{1} (k) h (n − k)

= f_{1} ∗ h Invariancia a traslación

Para probar que la convolución es un sistema invariante a traslación hacemos la convolución de la señal x(k) trasladada un valor k_{0}

y(n) = ∑_{k=−∞}^{∞} x (k − k_{0}) h (n − k)

haciendo k̂ = k − k_{0} tenemos:

∑_{k̂=−∞}^{∞} x \binom{}{k} h (n − (k̂ + k_{0}))

∑_{k̂=−∞}^{∞} x \binom{}{k} h ((n − k_{0}) − k̂)

= y(n − k_{0})

Linealidad

Para probar que se trata de un sistema lineal, calculamos la convolución de la suma de dos señales x_{1} y x_{2}

y(n) = ∑_{k=−∞}^{∞} [x_{1} (k) + x_{2} (k)] h (n − k)

= ∑_{k=−∞}^{∞} x_{1} (k) h (n − k) + ∑_{k=−∞}^{∞} x_{2} (k) h (n − k)

= y_{1}(n) + y_{2}(n)

5.3.2. Sucesión útil

Una sucesión que nos será especialmente útil para realizar los cálculos de convolución y correlación es la sucesión geométrica

s_{N} = ∑_{n=0}^{N} a^{n} = 1 + a + a^{2} + ... + a^{N}

para esta serie podemos ver

s_{N} = 1 + a(1 + a + a^{2} + ... + a^{N−1})
= 1 + as_{N−1}

que pasa si multiplicamos s_{N} por (1 − a)

(1 − a) s_{N} = (1 + a + a^{2} + ... + a^{N}) −

= (a + a^{2} + a^{3}... + a^{N+1})
= 1 − a^{N+1}

de lo cual concluimos que

s_{N} = \frac{1 − a^{N+1}}{1 − a}

Esta serie será convergente en el caso de que |a| < 1 y divergente en el caso de que |a| > 1. Lo cual lo podemos verificar haciendo

lím_{N→∞} a^{N} = 0

y el valor de convergencia lo podemos calcular con

s_{N} = 1 + as_{N−1}

en el lómite s_{N} = s_{N−1} = r

r = 1 + ar r(1 − a) = 1

r = \frac{1}{1 − a}

En el caso de que |r| = 1, la serie converge al valor de

s_{N} = ∑_{n=0}^{N} 1^{n} = 1 + 1 + 1^{2} + ... + 1^{N}

= N + 1

y en el caso de |r| > 1 la serie diverge

lím_{N→∞} a^{N} ≠ 0

5.3.3. Ejemplos

Ejemplo 1

Considere una señal x(n) = α^{n} y un kernel dado por h (n) = β^{n}. Calcular la convolución de estas dos señales.

y (n) = ∑_{k=0}^{N−1} α^{k}β^{(n−k)}

y (n) = ∑_{k=0}^{N−1} α^{k}β^{n}β^{−k}

y (n) = β^{n} ∑_{k=0}^{N−1} (\frac{α}{β})^{k}

y (n) = β^{n} \frac{1 − (\frac{α}{β})^{N}}{1 − \frac{α}{β}}

si quisieramos calcular la correlación hacemos

y (n) = ∑_{k=0}^{N−1} α^{k}β^{(k−n)}

y (n) = ∑_{k=0}^{N−1} α^{k}β^{−n}β^{k}

y (n) = β^{−n} ∑_{k=0}^{N−1} (αβ)^{k}

y (n) = β^{−n} \frac{1 − (αβ)^{N}}{1 − (αβ)}

Ejemplo 2

Mostrar que la convolución de una señal impulso unitario con un kernel cualquiera es el mismo kernel

δ(n) ∗ h(n) = h(n)

prueba Por definición tenemos

y(n) = ∑_{k=−∞}^{∞} x(k)h(n − k)

para nuestro caso sustituimos x(n) por la función impulso δ(n) dando lugar a

y(n) = ∑_{k=−∞}^{∞} δ(k)h(n − k)

recordemos que la función impulso, será igual a 1 solo cuando su argumento es cero, por lo cual

y(n) = ⋯ + 0 × h(n − 1) + 1 × h(n) + 0 × h(n + 1) + ⋯

y(n) = h(n)

Por esta propiedad al kernel se le conoce también como respuesta al impulso.

Ejemplo 3

Mostrar que dado y(n) = x(n) ∗ h(n) podemos hacer y′(n) = x(n) ∗ h′(n) Por definición tenemos

y(n) = ∑_{k=−∞}^{∞} x(k) ∗ h(n − k)

y que la derivada de una función discretizada puede ser representada por y′(n) = y(n) − y(n − 1) tenemos

y(n) − y(n − 1) = ∑_{k=−∞}^{∞} x(k) ∗ h(n − k) − ∑_{k=−∞}^{∞} x(k) ∗ h((n − 1) − k)

y(n) − y(n − 1) = ∑_{k=−∞}^{∞} x(k) ∗ [h(n − k) − h((n − 1) − k)]

y′(n) = ∑_{k=−∞}^{∞} x(k) ∗ h′(n − k)

5.4. Algunos kerneles interesantes

5.4.1. Suavizadores

Caja

Para implementar este suavizador, en una dimension, utilizamos una señal dada como h(n) = [0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0]

En dos dimensiones, podemos implementar este kernel de manera separada si hacemos la siguiente operación

( \begin{matrix}0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 \\ 0 & 1 & 1 & 1 & 0 \\ 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 \end{matrix} ) ∗ ( \begin{matrix}0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 \end{matrix} ) = ( \begin{matrix}0 & 0 & 0 & 0 & 0 \\ 0 & 1 & 1 & 1 & 0 \\ 0 & 1 & 1 & 1 & 0 \\ 0 & 1 & 1 & 1 & 0 \\ 0 & 0 & 0 & 0 & 0 \end{matrix} )

En la figura 5.2(a) y 5.2(d), se presentan gráficamente, los kerneles para suavizadores de caja en una y dos dimensiones

Binomial

Este kernel esta basado en los coeficientes binomiales. Una manera fácil de calcularlo es utilizar un triángulo de Pascal, de la siguiente forma.

1,1 1,2,1 1,3,3,1 1,4,6,4,1 1,5,10,10,5,1 1,6,15,20,15,6,1

Dado b_{1} = [1, 1] podemos calcular

b_{2} = b_{1} ∗ b_{1}
b_{3} = b_{2} ∗ b_{1}
b_{4} = b_{3} ∗ b_{1}
⋮
b_{k+1} = b_{k} ∗ b_{1}

Para su implementación este kernel es normalizado para que su suma de 1. En dos dimensiones de manera separable se puede implementar como

B_{3} = b_{3} ∗ b_{3}

Note que este kernel es recursivo.

\frac{1}{4} ( \begin{matrix}0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 \\ 0 & 1 & 2 & 1 & 0 \\ 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 \end{matrix} ) ∗ \frac{1}{4} ( \begin{matrix}0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 2 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 \end{matrix} ) = \frac{1}{16} ( \begin{matrix}0 & 0 & 0 & 0 & 0 \\ 0 & 1 & 2 & 1 & 0 \\ 0 & 2 & 4 & 2 & 0 \\ 0 & 1 & 2 & 1 & 0 \\ 0 & 0 & 0 & 0 & 0 \end{matrix} )

Gaussiano

Este kernel es creado utilizando una campana de Gauss. Par calcularla utilizamos la expresión

g(x) = \frac{1}{\sqrt{2π}σ} exp^{\frac{(x−μ)^{2}}{2σ^{2}}}

donde: μ es la media de la distribución y σ es la varianza. En nuestro caso la media la consideraremos cero y la varianza la utilizaremos como una estimación del tamaóo del kernel.

En dos dimensiones calcularemos nuestro kernel de manera desacoplada haciendo g(x, y) = g(x) ∗ g(y). En la figura 5.2(c) se muestra el kernel Gaussiano en una dimensión y 5.2(f) en dos dimensiones.

5.4.2. Derivadas

Derivada en x y y

La derivada de una función esta definida como

\frac{df(x)}{dx} = lím_{h→0} \frac{f(x) − f(x − h)}{h}

considerando un incremento unitario, que para el caso discreto es el más pequeño y que la señal fue capturada con una frecuencia de muestreo adecuada tenemos que la derivada la podemos aproximar por

\frac{df(x)}{dx} ≈ f(x) − f(x − h)

(a) Caja en 1d  (b) Binomial 1d  (c) Gaussiano 1d

(d) Caja 2d  (e) Binomial 2d  (f) Gaussiano 2d

Figura 5.2: Kerneles para suavizado de señales en una y dos dimensiones

De la expresión anterior podemos ver que el kernel de derivadas seró d(n) = [1 − 1], asó la derivada la calcularemos como y(n) ∗ d(n).

Para calcular las derivadas de orden N, solamente debemos recordar que se debe aplicar sucesivamente la derivada de orden uno, con ello damos lugar a la siguiente familia de kerneles.

+1,-1 +1,-2,+1 +1,-3,+3,-1 +1,-4,+6,-4,-1

Derivadas de Gaussianas

Un problema que presenta el cólculo de las derivadas utilizando convolución es que son amplificadores de ruido. Esto quiere decir que si la señal presenta ruido, la señal resultante de la convolución, el ruido será más notorio. Una manera de eliminar el ruido de una señal es aplicar un suavizador, el mós apropiado es el suavizado Gaussiano [Jahne, 2002]. Así pues, para eliminar el ruido, es deseable primero convolucionamos con un kernel gaussiano y luego aplicamos un kernel de derivadas de la siguiente forma

y(n) = (x(n) ∗ h(n)) ∗ d(n) = x(n) ∗ (h(n) ∗ d(n)) = x(n) ∗ \frac{dg(x)}{dx}

Resulta que tanto h(n) y d(n) son aproximaciones de una Gaussiana y una derivada respectivamente, y que ambos dan lugar a un nuevo kernel, pero es más inteligente calcular este kernel como la derivada de una función gaussiana. El kernel de derivada Gaussiano es

\frac{dg(x)}{dx} = \frac{−x}{σ^{2}} [\frac{1}{\sqrt{2π}σ} exp^{\frac{−x^{2}}{2σ^{2}}} ]

El kernel de segunda derivada es

\frac{d^{2}g(x)}{dx} = [\frac{x^{2}}{σ^{4}} − \frac{1}{σ^{2}} ] [\frac{1}{\sqrt{2π}σ} exp^{\frac{−x^{2}}{2σ^{2}}} ]

(a) Primer derivada 1d  (b) Segunda derivada 1d

Figura 5.3: Kerneles de derivadas de gaussianas con σ = 10

En dos dimensiones, la implementación de estos kerneles lo haremos de manera separable de la siguiente manera. Así la derivada en la dirección x será

g(x, y) = g(x) ∗ g(y)

\frac{∂g(x, y)}{∂x} = \frac{∂g(x)}{∂x} ∗ g(y)

\frac{∂g(x, y)}{∂y} = g(x) ∗ \frac{∂g(y)}{∂y}

Las segundas derivadas se calcularán como:

\frac{∂^{2}g(x, y)}{∂x^{2}} = \frac{∂^{2}g(x)}{∂x^{2}} ∗ g(y)

\frac{∂^{2}g(x, y)}{∂x∂y} = \frac{∂g(x)}{∂x} ∗ \frac{∂g(y)}{∂y}

\frac{∂^{2}g(x, y)}{∂y^{2}} = g(x) ∗ \frac{∂^{2}g(y)}{∂y^{2}}

En la figuras 5.4(a), 5.4(b) y 5.4(c), se presentan los kerneles de derivadas de Gaussianas en dos dimensiones.

(a) Kernel \frac{∂^{2}g(x,y)}{∂x^{2}} en 2d (b) Kernel \frac{∂^{2}g(x,y)}{∂x∂y} en 2d (c) Kernel \frac{∂^{2}g(x,y)}{∂y^{2}} en 2d

Figura 5.4: Kerneles de derivadas de gaussianas con σ = 10 en dos dimensiones

Laplaciano

El Laplaciano resulta de sumar las segundas derivadas en x y y de una señal bidimensional y esta dado por la siguiente expresión

Δ = \frac{∂^{2}g(x, y)}{∂x^{2}} + \frac{∂^{2}g(x, y)}{∂y^{2}}

Este kernel lo podemos construir sumando los kerneles de segundas derivadas de la siguiente forma

( \begin{matrix}0 & 0 & 0 \\ 1 & −2 & 1 \\ 0 & 0 & 0 \end{matrix} ) + ( \begin{matrix}0 & 1 & 0 \\ 0 & −2 & 0 \\ 0 & 1 & 0 \end{matrix} ) = ( \begin{matrix}0 & 1 & 0 \\ 1 & −4 & 1 \\ 0 & 1 & 0 \end{matrix} )

De manera robusta lo podemos calcular utilizando derivadas de Gaussiana, lo cual da lugar al kernel denominado ”Mexican Hat”. Este kernel se utiliza para determinar bordes en una imagen y podemos ver este kernel en la figura 5.5

Figura 5.5: Laplaciano. Mexican Hat

5.5. Respuesta de Sistemas lineales invariantes en el tiempo a exponenciales complejas.

La importancia de las exponenciales complejas en el estudios de sistemas LTI proviene del hecho, de que la respuesta de un sistema LTI a una entrada exponencial compleja es la misma exponencial compleja modificada solo en amplitud.

T [z^{n}] = H(z)z^{n}

donde el factor complejo de la amplitud H(z) seró en general una función de la variable compleja s. Esto lo podemos mostrar haciendo

y(n) = ∑_{k=−∞}^{∞} h (k) x (n − k)

y(n) = ∑_{k=−∞}^{∞} h (k) z^{(n−k)}

y(n) = ∑_{k=−∞}^{∞} h (k) z^{n}z^{−k}

y(n) = z^{n} ∑_{k=−∞}^{∞} h (k) z^{−k}

y(n) = H(z)z^{n}

donde

H(z) = ∑_{k=−∞}^{∞} h (k) z^{−k}

hemos demostrado que cualquier exponencial compleja es una función característica de un sistema LIT.

Si la entrada de un sistema LIT de tiempo discreto se presenta como una combinación de exponenciales complejas, esto es, si

x (n) = ∑_{k} a_{k}z_{k}^{n}

entonces la salida es

y (n) = ∑_{k} a_{k}H (z_{k}) z_{k}^{n}

Transformada de Fourier

6.1. Representación de señales periódicas.

Recordemos que la exponencial compleja e^{j(2π/N)n} es periódica con periodo N y existe una familia de exponenciales complejas también con periodo N dadas por

ϕ_{k} (n) = e^{jk(2π/N)n}

todas estas señales tienen frecuencias que son múltiplos de las misma frecuencia fundamental 2π/N. Así cuando k = 0 llamaremos a este el componente de CD de la señal, cuando k = 1 como armónico de frecuencia fundamenta, k = 2 segundo armónico y así sucesivamente.

Con esto queremos llegar a hacer la representación de una señal periódica como la combinación de exponenciales complejas, tal que

x (n) = ∑_{k=0}^{N−1} a_{k}ϕ_{k} (n) = ∑_{k=0}^{N−1} a_{k}e^{jk(2π/N)n}

esta ecuación es conocida como la serie de Fourier de tiempo discreto y los coeficientes a_{k} como los coeficientes de esta.

Para determinar los coeficientes de la serie de Fourier procedemos

x (0) = a_{0}ϕ_{0} (0) + a_{1}ϕ_{1} (0) + a_{2}ϕ_{2} (0) + ... + a_{N−1}ϕ_{N−1} (0)
x (1) = a_{0}ϕ_{0} (1) + a_{1}ϕ_{1} (1) + a_{2}ϕ_{2} (1) + ... + a_{N−1}ϕ_{N−1} (1)
x (2) = a_{0}ϕ_{0} (2) + a_{1}ϕ_{1} (2) + a_{2}ϕ_{2} (2) + ... + a_{N−1}ϕ_{N−1} (2)

x (N − 1) = a_{0}ϕ_{0} (N − 1) + a_{1}ϕ_{1} (N − 1) + ... + a_{N−1}ϕ_{N−1} (N − 1)

En forma matricial

[ \begin{matrix}ϕ_{0} (0) & ϕ_{1} (0) & … & ϕ_{N−1} (0) \\ ϕ_{0} (1) & ϕ_{1} (1) & … & ϕ_{N−1} (1) \\ ⋮ & ⋮ & ... & ⋮ \\ ϕ_{0} (N − 1) & ϕ_{1} (N − 1) & … & ϕ_{N−1} (N − 1) \end{matrix} ] [ a_{0} ; a_{1} ; ⋮ ; a_{N−1} ] = [ x_{0} ; x_{1} ; ⋮ ; x_{N−1} ]

Tenemos un sistema Φa = x para el cual la solución puede ser calculada si encontramos la inversa de la matriz Φ. Pero una alternativa es multiplicar por e^{−jr(2π/N)n} ambos lados de la representación en serie de Fourier y sumar para todos los N tórminos de la serie

∑_{n=0}^{N−1} x (n) e^{−jr(2π/N)n} = ∑_{n=0}^{N−1} ∑_{k=0}^{N−1} a_{k}e^{jk(2π/N)n}e^{−jr(2π/N)n}

∑_{n=0}^{N−1} x (n) e^{−jr(2π/N)n} = ∑_{k=0}^{N−1} a_{k} ∑_{n=0}^{N−1} e^{j(k−r)(2π/N)n}

la suma interior para n en el lado derecho de la ecuación es cero a menos que k − r sea cero o un móltiplo entero de N . Por tanto si escogemos valores de r sobre el mismo rango de k, tendremos que esta suma es igual a N si k = r y cero de lo contrario, por lo que tenemos

Na_{r} = ∑_{n=0}^{N−1} x (n) e^{−jr(2π/N)n}

a_{r} = \frac{1}{N} ∑_{n=0}^{N−1} x (n) e^{−jr(2π/N)n}

= \frac{1}{N} ∑_{n=0}^{N−1} x (n) ϕ_{r}^{∗} (n)

En lo general representaremos la transformada de Fourier como

X (k) = \frac{1}{N} ∑_{n=0}^{N−1} x (n) ϕ_{r}^{∗} (n)

Figura 6.1: Señal a descomponer en Fourier

6.1.1. Ejemplo

Encontrar los coeficientes de Fourier de la señal mostrada en la Fig. 6.1, cuyos elementos son x(n) = [1, 1, 1, 0, 0, 0, 0, 1, 1]^{T}.

Para encontrar los coeficientes de Fourier de la señal x, es necesario calcular el siguiente sistema de ecuaciones:

Φa_{k} = x_{n}

donde Φ es la matriz de exponenciales complejas exp^{j \frac{2π}{N}} kn (N es la cantidad de muestras de la señal, en este caso es 9) y es igual,

Φ = ( \begin{matrix}ϕ_{0}(0) & ϕ_{1}(0) & ϕ_{2}(0) & ϕ_{3}(0) & ϕ_{4}(0) & ϕ_{5}(0) & ϕ_{6}(0) & ϕ_{7}(0) & ϕ_{8}(0) \\ ϕ_{0}(1) & ϕ_{1}(1) & ϕ_{2}(1) & ϕ_{3}(1) & ϕ_{4}(1) & ϕ_{5}(1) & ϕ_{6}(1) & ϕ_{7}(1) & ϕ_{8}(1) \\ ϕ_{0}(2) & ϕ_{1}(2) & ϕ_{2}(2) & ϕ_{3}(2) & ϕ_{4}(2) & ϕ_{5}(2) & ϕ_{6}(2) & ϕ_{7}(2) & ϕ_{8}(2) \\ ϕ_{0}(3) & ϕ_{1}(3) & ϕ_{2}(3) & ϕ_{3}(3) & ϕ_{4}(3) & ϕ_{5}(3) & ϕ_{6}(3) & ϕ_{7}(3) & ϕ_{8}(3) \\ ϕ_{0}(4) & ϕ_{1}(4) & ϕ_{2}(4) & ϕ_{3}(4) & ϕ_{4}(4) & ϕ_{5}(4) & ϕ_{6}(4) & ϕ_{7}(4) & ϕ_{8}(4) \\ ϕ_{0}(5) & ϕ_{1}(5) & ϕ_{2}(5) & ϕ_{3}(5) & ϕ_{4}(5) & ϕ_{5}(5) & ϕ_{6}(5) & ϕ_{7}(5) & ϕ_{8}(5) \\ ϕ_{0}(6) & ϕ_{1}(6) & ϕ_{2}(6) & ϕ_{3}(6) & ϕ_{4}(6) & ϕ_{5}(6) & ϕ_{6}(6) & ϕ_{7}(6) & ϕ_{8}(6) \\ ϕ_{0}(7) & ϕ_{1}(7) & ϕ_{2}(7) & ϕ_{3}(7) & ϕ_{4}(7) & ϕ_{5}(7) & ϕ_{6}(7) & ϕ_{7}(7) & ϕ_{8}(7) \\ ϕ_{0}(8) & ϕ_{1}(8) & ϕ_{2}(8) & ϕ_{3}(8) & ϕ_{4}(8) & ϕ_{5}(8) & ϕ_{6}(8) & ϕ_{7}(8) & ϕ_{8}(8) \end{matrix} )

además,

a_{k} = ( \begin{matrix}a_{0} \\ a_{1} \\ a_{2} \\ a_{3} \\ a_{4} \\ a_{5} \\ a_{6} \\ a_{7} \\ a_{8} \end{matrix} )

y finalmente,

x_{n} = ( \begin{matrix}x(0) \\ x(1) \\ x(2) \\ x(3) \\ x(4) \\ x(5) \\ x(6) \\ x(7) \\ x(1) \end{matrix} )

Al resolver el sistema de ecuaciones se obtuvieron los coeficientes:

a_{k} = ( \begin{matrix}0,5556 \\ 0,3199 \\ −0,0591 \\ −0,1111 \\ 0,0725 \\ 0,0725 \\ −0,1111 \\ −0,0591 \\ 0,3199 \end{matrix} )

Note que a_{k} no contiene elementos imaginario. Esto se debe a que la señal es par y por lo tanto, se puede representar con suma de cosenos.

A continuación se va a reconstruir la señal original usando los coeficientes de Fourier obtenidos calculando,

y_{n} = Φa_{k}

Para poder observar la importancia de cada uno de los coeficientes encontrados, primero se va a reconstruir la señal usando un solo coeficiente de manera que a_{k} = [a_{0}, 0, 0, 0, 0, 0, 0, 0, 0]^{T}, después se van a usar dos coeficientes, de modo que a_{k} = [a_{0}, a1, 0, 0, 0, 0, 0, 0, 0]^{T}, y así sucesivamente hasta utilizar los 9 coeficientes. Los resultados de este ejercicio se muestran en las Figs. 6.2(a)-6.2(i).

(a)  (b)  (c)

(d)  (e)  (f)

(g)  (h)  (i)

Figura 6.2: Reconstrucción de la señal usando los coeficientes de Fourier

6.2. Propiedades de la Transformada de Discreta Fourier

La transformada de Fourier discreta, presenta algunas propiedades un cuanto diferentes que la transformada continua. A continuación se mencionan algunas y se da la demostración de tales

6.2.1. Periodicidad

Sea una señal x(n), mostrar que la transformada de Fourier, X(k) tiene periodo N.

Si la función es periódica entonces

X(k) = X(k + N)

La transformada de Fourier para ambas señales es

X(k) = \frac{1}{N} ∑_{n=0}^{N−1} x(n)e^{−j(\frac{2π}{N})nk}

X(k + N) = \frac{1}{N} ∑_{n=0}^{N−1} x(n)e^{−j(\frac{2π}{N})n(k+N)}

X(k + N) = \frac{1}{N} ∑_{n=0}^{N−1} x(n)e^{−j(\frac{2π}{N})nk}e^{−j(\frac{2π}{N})nN}

X(k + N) = \frac{1}{N} ∑_{n=0}^{N−1} x(n)e^{−j(\frac{2π}{N})nk}

6.2.2. Linealidad

Si

x_{1}(n) ↔^{F} X_{1}(k)

x_{2}(n) ↔^{F} X_{2}(k)

entonces

ax_{1}(n) + bx_{2}(n) ↔^{F} aX_{1}(k) + bX_{2}(k)

prueba

X(k) = \frac{1}{N} ∑_{n=0}^{N−1} x(n)e^{−j(\frac{2π}{N})nk}

para la suma de ax_{1}(n) + bx_{2}(n) queda

= \frac{1}{N} ∑_{n=0}^{N−1} [ax_{1}(n) + bx_{2}(n)] e^{−j(\frac{2π}{N})nk}

= \frac{a}{N} ∑_{n=0}^{N−1} x_{1}(n)e^{−j(\frac{2π}{N})nk} + \frac{b}{N} ∑_{n=0}^{N−1} x_{2}(n)e^{−j(\frac{2π}{N})nk}

= aX_{1}(k) + bX_{2}(k)

6.2.3. Desplazamiento en tiempo

Si

x(n) ↔^{F} X(k)

x(n − n_{0}) ↔^{F} e^{−j(\frac{2π}{N})n_{0}k}X(k)

demostración

X(k) = \frac{1}{N} ∑_{n=0}^{N} x(n − n_{0})e^{j(\frac{2π}{N})nk}

haciendo en cambio de variable n̂ = n − n_{0} tenemos

X(k) = \frac{1}{N} ∑_{n̂=−n_{0}}^{N−1−n_{0}} x(n̂)e^{−j(\frac{2π}{N})(n̂+n_{0})k}
X(k) = \frac{1}{N} ∑_{n̂=0}^{N−1} x(n̂)e^{−j(\frac{2π}{N})n̂k}e^{−j(\frac{2π}{N})n_{0}k}
X(k) = e^{−j(\frac{2π}{N})n_{0}k} \frac{1}{N} ∑_{n̂=0}^{N−1} x(n̂)e^{−j(\frac{2π}{N})n̂k}
X(k) = e^{−j(\frac{2π}{N})n_{0}k}X(k)

6.2.4. Desplazamiento en frecuencia

Ahora si hacemos un desplazamiento en frecuencia tenemos que

e^{j(\frac{2π}{N})nk_{0}}x(n) ↔^{F} X(k − k_{0})

demostración

X_{2}k = ∑_{k=0}^{N−1} X(k − k_{0})e^{j(\frac{2π}{N})nk}

haciendo el cambio de variable k̂ = k − k_{0} tenemos

X_{2}k = ∑_{k̂=−k_{0}}^{N−1−k_{0}} X(k̂)e^{j(\frac{2π}{N})n(k̂+k_{0})}
X_{2}k = = ∑_{k̂=0}^{N−1} X(k̂)e^{j(\frac{2π}{N})nk̂}e^{j(\frac{2π}{N})nk_{0}}
X_{2}k = = e^{j(\frac{2π}{N})nk_{0}} ∑_{k̂=0}^{N−1} X(k̂)e^{j(\frac{2π}{N})nk̂}
X_{2}k = = e^{j(\frac{2π}{N})nk_{0}}x(n)

6.2.5. Conjugación

La definición del conjugado de un número complejo es x = R + jI

x^{∗} = R − jI

Mostrar que:

x^{∗}(n) ↔^{F} X^{∗}(−k)

comenzaremos por

X_{2}(k) = \frac{1}{N} ∑_{n=0}^{N−1} x^{∗}(n)e^{−j(\frac{2π}{N})nk}
X_{2}(k) = \frac{1}{N} [∑_{n=0}^{N−1} x(n)e^{−j(\frac{2π}{N})n(−k)}]^{∗}
X_{2}(k) = X^{∗}(−k)

6.2.6. Inversión en Tiempo

Mostrar que

x(−n) ↔^{F} X(−k)

Comenzaremos por calcular la transformada de Fourier de x(−n)

X_{2}(k) = \frac{1}{N} ∑_{n=0}^{N−1} x(−n)e^{−j(\frac{2π}{N})nk}

hacemos el cambio de variable m = −n

X_{2}(k) = \frac{1}{N} ∑_{m=0}^{−N+1} x(m)e^{−j(\frac{2π}{N})(−m)k}
X_{2}(k) = \frac{1}{N} ∑_{m=0}^{N−1} x(m)e^{−j(\frac{2π}{N})m(−k)}
X_{2}(k) = X(−k)

6.2.7. Escalamiento en tiempo

Sea x(an) una señal multiplicada por un entero, calcular la transformada de Fourier para este caso

X_{2}(k) = \frac{1}{N} ∑_{n=0}^{N−1} x(an)e^{−j(\frac{2π}{N})nk}

haciendo m = an tenemos

X_{2}(k) = \frac{1}{N} ∑_{n=0}^{N−1} x(m)e^{−j(\frac{2π}{N})(\frac{m}{a})k}
X_{2}(k) = \frac{1}{N} ∑_{m=0}^{a(N−1)} x(m)e^{−j(\frac{2π}{N})m(\frac{k}{a})}
= X(\frac{k}{a})

Note que el rango de la TDF, también se escala.

6.2.8. Convolución

El teorema de la convolución afirma que si F [x (n)] = X(k) y F [y (n)] = Y (k) entonces

x (n) ∗ y (n) = NX(k)Y (k)

demostración:

Consideremos que la convolución g(n) = x (n) ∗ y (n) esté dada por

g(n) = ∑_{m=0}^{N−1} x (m) y (n − m)

y la transformada de fourier de g(n) es

G(k) = \frac{1}{N} ∑_{n=0}^{N−1} [g(n)] e^{−j(\frac{2π}{N})nk}

G(k) = \frac{1}{N} ∑_{n=0}^{N−1} [∑_{m=0}^{N−1} x (m) y (n − m)] e^{−j(\frac{2π}{N})nk}

cambiando el orden de la sumatoria

G(k) = \frac{1}{N} ∑_{m=0}^{N−1} x (m) [∑_{n=0}^{N−1} y (n − m) e^{−j(\frac{2π}{N})nk}]

haciendo n̂ = n − m

G(k) = \frac{1}{N} ∑_{m=0}^{N−1} x (m) [∑_{n=0}^{N−1} y \binom{}{n} e^{−j(\frac{2π}{N})(n̂+m)k}]

G(k) = \frac{1}{N} ∑_{m=0}^{N−1} x (m) e^{−j(\frac{2π}{N})mk} [∑_{n̂=0}^{N−1} y \binom{}{n} e^{−j(\frac{2π}{N})n̂k}]

G(k) = N [\frac{1}{N} ∑_{m=0}^{N−1} x (m) e^{−j(\frac{2π}{N})mk}] [\frac{1}{N} ∑_{n̂=0}^{N−1} y \binom{}{n} e^{−j(\frac{2π}{N})n̂k}]
G(k) = NX(k)Y (k)

6.2.9. Multiplicación

De manera inversa podemos ver que si F [x (n)] = X(k) y F [y (n)] = Y (k) entonces

x (n) y (n) = X(k) ∗ Y (k)

demostración:

Consideremos que la convolución G(k) = X(k) ∗ Y (k) está dada por

G(k) = ∑_{l=0}^{N−1} X (l) Y (k − l)

y su transformada de Fourier la podemos calcular como:

g(n) = ∑_{k=0}^{N−1} [G(k)] e^{j(\frac{2π}{N})nk}
g(n) = ∑_{k=0}^{N−1} [∑_{l=0}^{N−1} X (l) Y (k − l)] e^{j(\frac{2π}{N})nk}

cambiando el orden de la sumatoria

g(n) = ∑_{l=0}^{N−1} X (l) [∑_{k=0}^{N−1} Y (k − l) e^{j(\frac{2π}{N})nk}]

haciendo k̂ = k − l

g(n) = ∑_{l=0}^{N−1} X (l) [ ∑_{k̂=−l}^{N−1−l} Y \binom{}{k} e^{j(\frac{2π}{N})n(k̂+l)}]
g(n) = [∑_{l=0}^{N−1} X (l) e^{j(\frac{2π}{N})nl}] [∑_{k̂=0}^{N−1} Y \binom{}{k} e^{j(\frac{2π}{N})nk̂}]
g(n) = x(n)y(n)

6.2.10. Diferenciación en Tiempo

Demostrar que

x(n) − x(n − 1) ↔^{F} (1 − e^{−j(\frac{2π}{N})k}) X(k)

tenemos que

x(n) ↔^{F} X(k)

y

x(n − 1) ↔^{F} e^{−j(\frac{2π}{N})k}X(k)X(k)

por superposición demostramos la primera.

6.2.11. Diferenciación en Frecuencia

Demostrar que

X(k) − X(k − 1) ↔^{F} (1 − e^{j(\frac{2π}{N})k}) x(n)

tenemos que

X(k) ↔^{F} x(n)

y

X(n − 1) ↔^{F} e^{j(\frac{2π}{N})n}x(n)

por superposición demostramos la primera.

6.2.12. Propiedades de Simetría de la transformada de

Fourier.

Para cualquier secuencia x(n) la parte real la podemos calcular haciendo

R [x (n)] = \frac{1}{2} [x(n) + x^{∗}(n)]

y la parte imaginaria como

I [x (n)] = \frac{1}{2} [x(n) − x^{∗}(n)]

Si x(n) es una secuencia real entonces:

\frac{1}{2} [x(n) − x^{∗}(n)] = 0

lo que es equivalente a

x(n) = x^{∗}(n)

Al aplicar la transformada de Fourier tendremos que

X(k) = X^{∗}(−k)

A partir de esto se observa que la parte real de X (k) es una función par de k y la parte imaginaria de X (k) es una función impar de k. De manera similar, la magnitud de X (k) es una función par y el ángulo de de fase es una función impar. Además

E {x (n)} ↔^{F} R {X (n)}

y

O {x (n)} ↔^{F} I {X (n)}

demostración:

En el caso de una se estéal real y par se cumple que x(n) = x(−n) y que x(n) = x^{∗}(n) por lo tanto en el dominio de Fourier

X(k) = X(−k) = R [X(−k)] + jI [X(−k)]
X(k) = X^{∗}(−k) = R [X(−k)] − jI [X(−k)]

la única posibilidad de que esto sea cierto es que la parte imaginaria sea cero. Así la transformada de Fourier se transforma en

= \frac{1}{N} ∑_{n=0}^{N} x(n) cos (\frac{2π}{N} nk)

a esta ecuación es común que se le conozca como la transformada coseno.

En el caso de una señal real e impar, se cumple que x(n) = −x(−n) x(n) = x^{∗}(n), entonces en el dominio de Fourier

X(k) = −X(−k) = −R [X(−k)] − jI [X(−k)]
X(k) = X^{∗}(−k) = R [X(−k)] − jI [X(−k)]

la única posibilidad de que esto ocurra es que la parte real sea cero. La transformada de Fourier en este caso queda como

= j \frac{1}{N} ∑_{n=0}^{N} x(n)sen (\frac{2π}{N} nk)

6.2.13. Relación de Parseval.

Si x (n) y X (k) son la se estéal y su transformada respectivamente, entonces tenemos que:

∑_{n=0}^{N−1} |x(n)|^{2} ↔^{F} N ∑_{k=0}^{N−1} |X(k)|^{2}

podemos representar como

∑_{n=0}^{N−1} x(n)x^{∗}(n) = N ∑_{k=0}^{N−1} X(k)X(k)^{∗}

= N ∑_{k=0}^{N−1} X(k) [\frac{1}{N} ∑_{n=0}^{N−1} x(n)e^{−j(\frac{2π}{n})nk}]^{∗}

= N ∑_{k=0}^{N−1} X(k)\frac{1}{N} ∑_{n=0}^{N−1} x^{∗}(n)e^{j(\frac{2π}{n})nk}

cambiando el orden en que se hacen las sumatorias

= ∑_{n=0}^{N−1} x^{∗}(n) [∑_{k=0}^{N−1} X(k)e^{j(\frac{2π}{n})nk}]

= ∑_{n=0}^{N−1} x^{∗}(n)x(n)

= ∑_{n=0}^{N−1} |x(n)|^{2}

6.2.14. Resumen de Propiedades

Dado lo anterior podemos resumir en la siguiente tabla las propiedades de la transformada discreta de Fourier.

Propiedad  Señal en  Señal en
Propiedad  el tiempo  la frecuencia
Periodicidad  x(n)  X(k) periodo N
Linealidad  ax_{1}(n) + bx_{2}(n)  aX_{1}(k) + bX_{2}(k)
Desp. en tiempo  x(n − n_{0})  e^{−j(\frac{2π}{N})n_{0}k}X(k)
Desp. en frecuencia  e^{j(\frac{2π}{N})nk_{0}}x(n)  X(k − k_{0})
Conjugación  x^{∗}(n)  X^{∗}(−k)
Inversión en tiempo  x(−n)  X(−k)
Expansión en Tiempo  x(an)  X(k/a)
Convolución  x(n) ∗ y(n)  NX(k)Y (k)
Multiplicación  x(n)y(n)  X(k) ∗ Y (k)
Dif. en tiempo  x(n) − x(n − 1)  (1 − e^{−j(\frac{2π}{N})k}) X(k)
Dif. en Frecuencia (1 − e^{j(\frac{2π}{N})k}) x(n)  X(k) − X(k − 1)
Simetría conjugada x(n) real X(k) = X^{∗}(−k)

X_{r}(k) = X_{r}(−k) par X_{i}(k) = −X_{i}(−k) impar

Simetría x(n) real y par X(k) real y par
Simetría  x(n) real e impar X(k) imaginaria e impar
Parseval  ∑_{n=0}^{N−1} |x(n)|^{2}  N ∑_{k=0}^{N−1} |X(k)|^{2}

6.3. Transformada de Fourier de algunas funciones interesantes

6.3.1. Exponencial Compleja

Dado x(n) = e^{j(\frac{2π}{N})k_{0}n} la transformada de fourier esta dada como

X(k) = \frac{1}{N} ∑_{n=0}^{N−1} x(n)e^{−j(\frac{2π}{N})kn}

= \frac{1}{N} ∑_{n=0}^{N−1} e^{j(\frac{2π}{N})k_{0}n}e^{−j(\frac{2π}{N})kn}

= \frac{1}{N} ∑_{n=0}^{N−1} e^{j(\frac{2π}{N})(k_{0}−k)n}

esta suma tendrá solución N en el caso de que k = k_{0} y cero en los demás casos. La transformada de Fourier es

e^{j(\frac{2π}{N})k_{0}n} ↔^{F} δ (k − k_{0})

6.3.2. Función seno

Dada la función seno, podemos hacer la representación de esta función como

sen (\frac{2π}{N} k_{0}n) = \frac{1}{2j} [e^{j(\frac{2π}{N})k_{0}n} − e^{−j(\frac{2π}{N})k_{0}n}]

dado que tenemos la suma de dos exponenciales complejas, la transformada de Fourier es

sen (\frac{2π}{N} k_{0}n) ↔^{F} −\frac{j}{2}δ (k − k_{0}) + \frac{j}{2}δ (k + k_{0})

6.3.3. Función coseno

La representación de la función coseno la podemos dar como

cos (\frac{2π}{N} k_{0}n) = \frac{1}{2} [e^{j(\frac{2π}{N})k_{0}n} + e^{−j(\frac{2π}{N})k_{0}n}]

de manera similar que en la función seno tenemos

cos (\frac{2π}{N} k_{0}n) ↔^{F} \frac{1}{2}δ (k − k_{0}) + \frac{1}{2}δ (k + k_{0})

6.3.4. Función impulso unitario

Recordemos que la función impulso unitario es

δ (n) = \{ \begin{matrix}1 si n = 0 \\ 0 en caso contrario \end{matrix}

la transformada de Fourier para esta la calculamos

X(k) = \frac{1}{N} ∑_{n=0}^{N−1} δ (n − n_{0}) e^{−j(\frac{2π}{N})kn}

= \frac{1}{N} e^{−j(\frac{2π}{N})kn_{0}}

podemos comprobar que |e^{−j(\frac{2π}{N})kn_{0}} | = 1 para cualquier valor que tome k por lo que finalmente la magnitud de la transformada de Fourier es

δ (n − n_{0}) ↔^{F} \frac{1}{N}

6.3.5. Constante

Consideremos una función x (n) = a la transformada de Fourier es

X(k) = \frac{1}{N} ∑_{n=0}^{N−1} ae^{−j(\frac{2π}{N})kn}

= \frac{a}{N} ∑_{n=0}^{N−1} e^{−j(\frac{2π}{N})kn}

la sumatoria tendró la solución N para k = 0 y 0 para k ≠ 0 por lo cual

a ↔^{F} aδ (k)

6.3.6. Escalón Unitario

La función escalón unitario esta definida como u (n_{0}) = \{ \begin{matrix}1 si n >= n_{0} \\ 0 en caso contrario \end{matrix}

la transformada de Fourier estará dada como

U(k) = \frac{1}{N} ∑_{n=0}^{N−1} u(n_{0})e^{−j(\frac{2π}{N})kn}

U(k) = \frac{1}{N} ∑_{n=n_{0}}^{N−1} e^{−j(\frac{2π}{N})kn}

U(k) = \frac{1}{N} ∑_{n=0}^{N−1} e^{−j(\frac{2π}{N})kn} − \frac{1}{N} ∑_{n=0}^{n_{0}−1} e^{−j(\frac{2π}{N})kn}

Lo cual da como resultado U(k) = \{ \begin{matrix}\frac{N−n_{0}}{N} si k = 0 \\ Û(k) en caso contrario \end{matrix} con

Û(k) = −\frac{1 − e^{−j(\frac{2π}{N})kn_{0}}}{N(1 − e^{−j(\frac{2π}{N})k})}
Û(k) = −\frac{e^{−j(\frac{π}{N})kn_{0}} (e^{j(\frac{π}{N})kn_{0}} − e^{−j(\frac{π}{N})kn_{0}})}{Ne^{−j(\frac{π}{N})k} (e^{j(\frac{π}{N})k} − e^{−j(\frac{π}{N})k})}
Û(k) = −\frac{1}{N} \frac{sin(\frac{πn_{0}k}{N})}{sin(\frac{πk}{N})}e^{−j(\frac{π}{N})k(n_{0}−1)}

6.3.7. Caja

Consideremos una función x (n) dada por la siguiente expresión

c(n) = \{ \begin{matrix}1 si -d < n < d \\ 0 si no \end{matrix}

Aplicando la formula podemos ver

C(k) = \frac{1}{N} ∑_{n=−N/2}^{N/2} c(n)e^{−j(\frac{2π}{N})kn}
C(k) = \frac{1}{N} ∑_{n=−d}^{d} e^{−j(\frac{2π}{N})kn}

Podemos organizar los términos de la siguiente forma

C(k) = \frac{1}{N} ∑_{n=−d}^{−1} e^{−j(\frac{2π}{N})kn} + \frac{1}{N} ∑_{n=0}^{d} e^{−j(\frac{2π}{N})kn}
C(k) = \frac{1}{N} ∑_{n=1}^{d} e^{−j(\frac{2π}{N})k(−n)} + \frac{1}{N} ∑_{n=0}^{d} e^{−j(\frac{2π}{N})kn}
C(k) = \frac{1}{N} ∑_{n=0}^{d} e^{j(\frac{2π}{N})kn} − \frac{1}{N} + \frac{1}{N} ∑_{n=0}^{d} e^{−j(\frac{2π}{N})kn}

Lo cual da como resultado

C(k) = \{ \begin{matrix}si k = 0 \\ ...... & \frac{2d+1}{N} \\ si no \\ ...... & \frac{1}{N} & (\frac{1−e^{\frac{2π}{N} k(d+1)}}{1−e^{\frac{2π}{N} k}} + \frac{1−e^{− \frac{2π}{N} k(d+1)}}{1−e^{− \frac{2π}{N} k}} − 1) \end{matrix}

En la figura 6.3 se presenta la parte real de esta función ya que la imaginaria es cero.

La transformada de Fourier continua de esta función, esta dada por la siguiente expresión

C(k) = \frac{1}{N} ∫_{−∞}^{∞} c(t)e^{−j( \frac{2π}{N})kt}dt

al sustituir nuestra función obtenemos

C(k) = \frac{1}{N} ∫_{−d}^{d} e^{−j( \frac{2π}{N})kt}dt
C(k) = \frac{1}{N} \frac{1}{[−j(\frac{2π}{N})k]} e^{−j( \frac{2π}{N})kt}|_{−d}^{d}
C(k) = \frac{1}{N} \frac{1}{[−j(\frac{2π}{N})k]} (−2j) sen [\frac{2π}{N} kd]
C(k) = \frac{sen [\frac{2π}{N}kd]}{πk}

la cual luce como 6.4

6.3.8. Ejemplos

En un programa en Java que utilice las clases funciones y gráfica, obtener y explicar brevemente la Transformada de Fourier Discreta (TFD) de las siguientes señales:

1. x_{1}(n) = cos(\frac{2π}{N}n) + sin(\frac{10π}{N}n)

2. x_{2}(n) = δ(n) + δ(n − 3)

3. x_{3}(n) = 1 + δ(n + 5) + cos(\frac{6π}{N}n) + exp^{−j( \frac{2π}{N})n} Figura 6.3: Función Sinc con d = 10 y N = 128.

Figura 6.4: Función Sinc con k = 10 y N = 128.

Ejercicio 1

La figura 6.3.8 muestra las cuatro gráficas que se obtienen al aplicar la TFD en la señal 1. Se puede observar en la figura 6.5(a) la señal que contiene parte real pero no contiene parte imaginaria, además la parte real se forma por un seno y un coseno. El resultado de aplicar la TFD a esta señal equivale a aplicar la TFD a las dos señales que la componen por separado y luego sumarlas. En las figuras 6.5(b) y 6.5(c) se puede observar el resultado de la aplicación de la TFD a la señal 1, cuya solución analítica es:

X_{1}(k) = \frac{1}{2}δ(k − 1) + \frac{1}{2}δ(k + 1) − \frac{j}{2}δ(k − 5) + \frac{j}{2}δ(k + 5)

(a) x(n) 1 (b) Parte real de la FFT(x(n))

(c) Parte imaginaria de la FFT(x(n))

Figura 6.5: Señal 1 y su transformada de fourier

Con las figuras podemos constatar que el coseno efectivamente contribuye con dos impulsos positivos en k = 1 y k = −1. Asimismo, el seno contribuye con dos pulsos con diferente signo en la parte imaginaria de la señal, el positivo en k = 5 y el negativo en k = −5.

Ejercicio 2

La solución analítica de la señal 2 es:

X_{2}(k) = \frac{1}{N} + \frac{1}{N} cos(\frac{6π}{N} k) − \frac{j}{N} sin(\frac{6π}{N} k)

Las figura 6.3.8 muestra las gróficas de la señal original su TFD. Observando las gráficas 6.6(b) y 6.6(c) se puede concluir que existe un coseno con frecuencia 3 en la parte real y un seno con igual frecuencia en la parte imaginaria. La contribución del impulso en el origen casi no se percibe en las gróficas, ya que es es una constante que se suma al resto de la señal.

Ejercicio 3

Las figuras 6.7(a) y 6.7(b) muestra la parte real e imaginaria respectivamente de la señal 3. Al aplicar la TFD de manera analótica a la señal se obtiene:

(a) Señal real (b) Parte real de la FFT(x(n))

(c) Parte imaginaria de la FFT(x(n))

Figura 6.6: Señal 2 y su transformada de fourier

X_{3}(k) = δ(k) + \frac{1}{N} cos(\frac{6π}{N} k) − \frac{j}{N} sin(\frac{6}{N}) + δ(k + 1)

La contribución de la constante se puede observar en la grófica de la figura6.7(c) como un impulso en 1 de la parte real. La contribución del impulso es coseno con frecuencia 5 en la parte real y un seno negativo con frecuencia 5 en la parte imaginaria (figura 6.7(d)). El coseno contribuye con un par de impulsos en la parte real, cuando k = 3 y k = −3, y finalmente el exponencial complejo contribuye con un pulso en k = −1.

6.4. Transformada de Fourier en dos dimensiones.

En dos dimensiones la transformada de Fourier se expresa como

X(k, l) = \frac{1}{NM} ∑_{n=0}^{N−1} ∑_{m=0}^{M−1} x(n, m)e^{−j(\frac{2π}{NM})(nk+ml)}

reorganizando términos tenemos:

X(k, l) = \frac{1}{NM} ∑_{m=0}^{M−1} [∑_{n=0}^{N−1} x(n, m)e^{−j(\frac{2π}{NM})nk}] e^{−j(\frac{2π}{NM})ml}

de esta expresión podemos ver que la transformada de Fourier en dos dimensiones, es equivalente a realizar la transformada de Fourier en una dirección del arreglo bidimensional y despuós realizarlo en la otra dirección.

(a) Señal real en el tiempo (b) Señal imaginaria en el tiempo (c) Señal real en la frecuencia

(d) Señal imaginaria en la imaginaria

Figura 6.7: Señal 3 y su transformada de fourier

La trasformada Discreta de Fourier TDF es O(N^{2}) y existe la Transformada Rópida de Fourier FFT la cual es O(N log N), para una señal en una sola dimensión. Para el caso de señales en dos dimensiones, la TDF es O((NM)^{2}) y la FFT es O(NM log NM), razón por la cual se utiliza la FFT. Sin embargo, la FFT requiere que el tamaño de la señales sea potencia de 2, es decir N = 2^{r} y M = 2^{s}. Todas la propiedades de la transformada de Fourier son validas para el caso multidimensional y para hacer la demostración solo hay que tener en cuenta que se trata de arreglos bidimensionales.

Transformada de Fourier de una señal bidimensional

Dada la señal x (n.m) = cos [\frac{2∗π}{N} (k_{0}n + l_{0}m)] determinar su transformada de Fourier.

Esta señal puede ser representada como la suma de dos exponenciales complejas

x(n, m) = cos [\frac{2 ∗ π}{N} (k_{0}n + l_{0}m)]

= \frac{1}{2}e^{j(\frac{2π}{N})(k_{0}n+l_{0}m)} + \frac{1}{2}e^{−j(\frac{2π}{N})(k_{0}n+l_{0}m)}

por lo tanto

X(k, l) = \frac{1}{2}δ (k − k_{0}, l − l_{0}) + \frac{1}{2}δ (k + k_{0}, l + l_{0})

En la figura 6.8 podemos ver la señal bidimensional correspondiente x(n, m) = cos [\frac{2∗π}{N} (10n + 20m))] y su transformada de Fourier. Note que en la imagen de la transformada de Fourier aparecen dos picos en las coordenadas [10, 20] y [118, 108], lo cual, corresponde con la deducción anterior.

Figura 6.8: Señal coseno en dos dimensiones

Ejemplo

El primer renglón de una imagen es creado utilizando una distribución normal N(0, 1) (ver 6.10) y los renglones subsiguientes se calculan de acuerdo con la siguiente sucesión f(n, m) = f (n − 1, m − 1) tal que se produce un patron como el mostrado en la figura 6.9 (a la izquierda). Explique a que se debe que la magnitud de la transformada de Fourier esta dominado por una linea recta a 45 grados. La transformada de Fourier del ruido gaussiano la podemos ver a la derecha de la figura 6.9 .

Figura 6.9: A la izquierda Patron regular creado con una señal uni– dimensional trasladada (ruido Gaussiano) y a la derecha Magnitud de su transformada de Fourier

Figura 6.10: Valores aleatorios calculado con una distribución Normal con media cero y varianza unitaria

Demostración:

El primer renglón de imagen fue generado utilizando una distribución normal, y los renglones subsecuentes utilizando la recursión

f(n, m) = f (n − 1, m − 1)

la transformada de Fourier para f(n, m) es F (k, l) si aplicamos la propiedad de translación de la TF, podemos ver:

f(n, m) ↔^{F} F (k, l)

f(n − n_{0}, m − m_{0}) ↔^{F} e^{−j(\frac{2π}{NM})(n_{0}k+m_{0}l)}F (k, l)

lo cual nos da como resultado la transformada de fourier F (x), esta multiplicada por una exponencial compleja e^{−j(\frac{2π}{NM})(n_{0}k+m_{0}l)} cuyo móximo se localiza sobre la recta n_{0}k + m_{0}l = 0, dado que n_{0} = m_{0} = 1, tendremos que la pendiente de esta linea recta es 45 grados.

6.5. Convolución utilizando TF.

En una imagen tenemos que la referencia se encuentra en la esquina superior izquierda. Cuando estamos utilizando la transformada de Fourier debemos recordar que estó, considera que las señales tienen periodos N en la dirección de x y M en la dirección de y. Así, una imagen estará representada virtualmente, por un conjunto de copias infinitas es un espacio infinito.

Para realizar una translación de la imagen, una manera es convolucionar la imagen con la función delta de Dirac δ (n − n_{0}, m − m_{0}) donde n_{0} y m_{0} es el desplazamiento, hay que recordar que la convolución es Lineal e invariante a la translación, por lo que, la convolución de dos señales no se vera afectada cuando hacemos la translación del kernel.

Cuando aplicamos el teorema de la convolución

x(n, m) ∗ y(n, m) ↔^{F} X(k, l)Y (K, l)

debemos tomar en cuenta que si el kernel se desplaza, la transformada de Fourier del kernel estará multiplicada por una exponencial compleja

x(n − n_{0}, y − m_{0}) ↔^{F} X(k, l)e^{−j( \frac{2π}{N} n_{0}k+ \frac{2π}{M} m_{0}l)}

x(n, m) ∗ δ (n − n_{0}, m − m_{0}) ↔^{F} X(k, l)e^{−j( \frac{2π}{N} n_{0}k+ \frac{2π}{M} m_{0}l)}

Cuando aplicamos el teorema de la convolución tendremos

x(n, m) ∗ [y(n, m) ∗ δ (n − n_{0}, m − m_{0})] ↔^{F} X(k, l) [Y (k, l)e^{−j( \frac{2π}{N} n_{0}k+ \frac{2π}{M} m_{0}l)}]

[x(n, m) ∗ δ (n − n_{0}, m − m_{0})] ∗ y(n, m) ↔^{F} [X(k, l)e^{−j( \frac{2π}{N} n_{0}k+ \frac{2π}{M} m_{0}l)}] Y (k, l)

x(n − n_{0}, m − m_{0}) ∗ y(n, n) ↔^{F} [X(k, l)e^{−j( \frac{2π}{N} n_{0}k+ \frac{2π}{M} m_{0}l)}] Y (k, l)

lo cual significa que x(n − n_{0}, m − m_{0}), estará desplazada a las coordenadas [n_{0}, m_{0}]

6.6. Teorema del Muestreo

Consideramos el caso de una señal continua en el tiempo dada por x(t) = cos(15πt). Es claro que la señal tiene una frecuencia angular w = 15π, una frecuencia f = 7,5 y un periodo T = \frac{1}{7,5}. Esto significa que si graficamos la señal tendremos siete ciclos y medio por segundo y la grófica para esta lucirá como la figura 6.11(a). Si esta misma función la submuestreamos con incrementos de \frac{π}{9} lucirá como la fig 6.11(b) y esta misma a intervalos de 0.01 lucirá como 6.11(c). Cual es el valor de incremento al que debo discretizar mi señal sin perder información?.

(a) Señal original (b) Señal muestreada a \frac{π}{9} (c) Señal muestreada a 0.01

Figura 6.11: Señal cos(15πt) con diferentes muestreos

6.6.1. Transformada de Fourier de una señal periódica

Supongamos que una señal con periodo T_{0} (con ω_{0} = \frac{2π}{T_{0}}), pudiera representarse como la serie

x(t) = ∑_{k=−∞}^{k=∞} a_{k}e^{jω_{0}kt}  (6.1)

multiplicando ambos miembros de la ecuación e integrando en el intervalo

∫− \frac{T_{0}}{2}^{\frac{T_{0}}{2}} x(t)e^{−jnω_{0}t}dt = ∫− \frac{T_{0}}{2}^{\frac{T_{0}}{2}} ∑_{k=−∞}^{k=∞} a_{k}e^{jω_{0}kt}e^{−jnω_{0}t}dt

reagrupando los términos tenemos

∫− \frac{T_{0}}{2}^{\frac{T_{0}}{2}} x(t)e^{−jnω_{0}t}dt = ∑_{k=−∞}^{k=∞} a_{k} [∫− \frac{T_{0}}{2}^{\frac{T_{0}}{2}} e^{jω_{0}(k−n)t}] dt

La evaluación de la integral en los corchetes es, considerando que k ≠ n:

∫− \frac{T_{0}}{2}^{\frac{T_{0}}{2}} e^{jω_{0}(k−n)t}dt = \frac{1}{j(k − n)ω_{0}} e^{jω_{0}(k−n)t}|− \frac{T_{0}}{2}^{\frac{T_{0}}{2}}
= \frac{1}{j(k − n)ω_{0}} [e^{j(k−n)π} − e^{−j(k−n)π}]
= \frac{1}{j(k − n)ω_{0}} [−1 + 1]
= 0

en el caso de que k = n, la integral se transforma en

∫− \frac{T_{0}}{2}^{\frac{T_{0}}{2}} dt = T_{0}

Por lo tanto la transformada de Fourier de una señal periódica ecuación (6.1) es

a_{n} = ∫− \frac{T_{0}}{2}^{\frac{T_{0}}{2}} x(t)e^{−jnω_{0}t}dt  (6.2)

6.6.2. Transformada de Furier de un Tren de Pulsos

Un tren de pulsos lo podemos definir por

p(t) = ⋯ + δ (t + 2T) + δ (t + T) + δ (t) +
δ (t − T) + δ (t + 2T) + ⋯

la cual podemos representar en forma compacta como

p(t) = ∑_{k=−∞}^{∞} δ (t − kT)

donde T = \frac{2π}{ω}. Note que esta ecuación es una función periódica y su transformada puede ser calculada utilizando el resultado de la ecuación (6.2)

En este caso hablamos de una transformada de Fourier continua por lo cual haremos la transformada de Fourier utilizando la definición continua

P (ω) = ∫_{−∞}^{+∞} p(t)e^{−jωt}dt

sustituyendo tenemos

P (ω) = ∑_{k=−∞}^{∞} ∫_{−∞}^{+∞} δ (t − kT) e^{−jωt}dt  (6.3)

Aplicando el resultado de la ecuación (6.2) tenemos

P (ω) = ∑_{k=−∞}^{∞} \frac{1}{T} ∫− \frac{T_{0}}{2}^{\frac{T_{0}}{2}} δ (t − kT) e^{−jωt}dt

de donde se desprende que si k = 0

\frac{1}{T} ∫− \frac{T_{0}}{2}^{\frac{T_{0}}{2}} δ (t − kT) e^{−jωt}dt = \frac{1}{T} = \frac{1}{T} δ(ω)

entonces da como resultado un tren de pulsos equiespaciados en la frecuencia con ω_{s} = \frac{2π}{T}.

P (ω) = ∑_{k=−∞}^{∞} δ (ω − kω_{s})

6.6.3. Señal muestreada

Una señal muestreada se calcula como

x_{p}(t) = x(t)p(t)

donde p(t) es un tren de pulsos. Aplicando el teorema de la convolución tenemos que la transformada de Fourier de esta señal es

X_{p}(ω) = X(ω) ∗ P (ω)

X_{p}(ω) = X(ω) ∗ ∑_{k=−∞}^{∞} δ (ω − ω_{s})

= ∑_{k=−∞}^{∞} X(ω − kω_{s})

lo cual significa que tendremos copias de nuestra señal X(ω) repetidas cada \frac{2π}{T}. Pero que ocurre si el ancho de banda de la señal X(ω) sobre pasa el valor \frac{2π}{T}?.

Llamemos ω_{s} la frecuencia de muestreo la cual esta dada como

ω_{s} = \frac{2π}{T}

y ω_{m} el ancho de banda de la señal x(t), podemos ver que si

ω_{s} > 2ω_{m}

podemos reconstruir exactamente la señal continua x (t) a partir de la señal muestreada x_{p} (t). La frecuencia de muestreo ω_{s} tambión se conoce como la frecuencia de Nyquist.

Podemos notar que

x_{p}(n) = x(nT) = \frac{1}{2π} ∫_{−∞}^{∞} X (ω) e^{jωnT} dω

dado que la señal en el dominio de la frecuencia es periódica podemos considerar que

x_{p}(n) = \frac{1}{2π} ∑_{r=−∞}^{∞} ∫_{(2r−1)π/T}^{(2r+1)π/T} X (ω) e^{jωnT} dω

y cada tórmino de la suma puede ser reducido a una integral sobre el rango −\frac{π}{T} a \frac{π}{T} haciendo el cambio de variable ω = Ω + \frac{2πr}{T}

x_{p}(n) = \frac{1}{2π} ∑_{r=−∞}^{∞} ∫_{−π/T}^{π/T} X (Ω + \frac{2πr}{T}) e^{jΩnT} e^{j2πrn}dΩ

Si intercambiamos el orden de integración y considerando que e^{j2πrn} = 1 para cualquier valor de n y r, tenemos que:

x_{p}(n) = \frac{1}{2π} ∫_{−π/T}^{π/T} [ ∑_{r=−∞}^{∞} X (Ω + \frac{2πr}{T}) e^{jΩnT} ] dΩ

haciendo un cambio de variable Ω = \frac{ω}{T}

x_{p}(n) = \frac{1}{2π} ∫_{−π}^{π} [\frac{1}{T} ∑_{r=−∞}^{∞} X (\frac{ω}{T} + \frac{2πr}{T}) e^{jωn}] dω

de esta ecuación tenemos que

X_{p} (ω) = \frac{1}{T} ∑_{r=−∞}^{∞} X (\frac{ω}{T} + \frac{2πr}{T})

o si hacemos ω̂ = \frac{ω}{T} tenemos

X_{p} (ω̂T) = \frac{1}{T} ∑_{r=−∞}^{∞} X (ω̂ + \frac{2πr}{T})

6.6.4. Integración de la señal continua.

Ahora que hacemos para recuperar la señal continua a partir de la señal muestreada. Puesto que la señal X_{p}(ω) solo esta definida en el rango [−\frac{π}{T} , \frac{π}{T} ], calcularemos la antitrasformada en estos límites, además X_{p} (ωT) = \frac{1}{T} X (ω)

x(t) = \frac{1}{2π} ∫_{− \frac{π}{T}}^{+ \frac{π}{T}} T X_{p} (ωT) e^{jωt}dω

= \frac{T}{2π} ∫_{− \frac{π}{T}}^{+ \frac{π}{T}} [ ∑_{n=−∞}^{∞} x (nT) e^{−jωnT} ] e^{jωt}dω

= ∑_{n=−∞}^{∞} x (nT) [\frac{T}{2π} ∫_{− \frac{π}{T}}^{+ \frac{π}{T}} e^{jω(t−nT)}dω]

= ∑_{n=−∞}^{∞} x (nT) \frac{sen (\frac{π}{T} ) (t − nT)}{(\frac{π}{T} ) (t − nT)}

El código para realizar esta integración es:

  public static double Interpola(double x[], double t, double T)
  {
    int n, N = x.length, i;
    double val =0;
    for(i=-3*N; i<3*N; i++)
    {
      n = (i%N);
      if(n<0) n = N + n;

      val += x[n] * Sinc( (t - (double) i * T) * Math.PI / T);
    }
    return val;
  }

6.6.5. Ejemplos

Ejemplo 1

Aquella frecuencia que de acuerdo con el teorema de muestreo, debe ser excedida por la frecuencia de muestreo se llama razón de Nyquist. Determine la razón de Nyquist para las siguientes señales

a) x(t) = 1 + cos(2000πt) + sen(4000πt).

La transformada de Fourier para esta señal es

x(k) = δ(0)+\frac{1}{2} [δ(w − 2000π) + δ(w + 2000π)]+\frac{j}{2} [δ(w − 4000π) − δ(w + 4000π)]

lo cual indica que el espectro de frecuencia tendrá valores en el intervalo [−4000π, 4000π], o bien que ω_{b} = 4000π. De acuerdo con esto la razón de Nyquist es ω_{s} = 8000π.

b) x(t) = \frac{sen(200πt)}{πt}. Podemos ver que la antitransformada de Fourier de una caja es

x(t) = A ∫_{−d}^{d} e^{jωt}dω

x(t) = \frac{A}{jt} e^{jωt} |_{−d}^{d}

x(t) = \frac{1}{jt} (2j) sen [dt]

x(t) = \frac{2Asen [dt]}{t}

de donde podemos ver que la transformada de Fourier de esta señal es una caja de ancho d = 200π y altura A = \frac{1}{2π}, por lo tanto su ancho de banda seró de ω_{b} = 200π y su razón de Nyquist seró ω_{s} = 400π.

Dada la señal continua x(t) = sen^{2}(10πt), calcular:

a) El ancho de banda correspondiente a esta señal

b) La transformada discreta de Fourier si la señal es muestreada en con periodo T_{s} = \frac{5}{2N} segundos, con un nómero de muestras N = 128.

c) Cual es la transformada discreta de fourier si la señal es muestreada con un periodo del doble del anterior.

Para esta señal aplicando identidades trigonomótricas tenemos que:

x(t) = sen^{2}(10πt) = \frac{1}{2} (1 − cos(20πt))

Respuestas

a) La frecuencia del primer termino es ω_{0} = 0 y del segundo termino es ω_{1} = 20π, con lo cual tenemos que el ancho de banda es ω_{b} = 20π

b)En este caso la frecuencia de muestreo es ω_{s} = \frac{2π}{\frac{5}{2N}} = \frac{4Nπ}{5} considerando 128 muestras la frecuencia de muestreo es ω_{s} = 321,6991, note que

ω_{s} >= 2 ∗ ω_{b}3
321,6921 >= 40π

Nuestra señal cumple con el teorema del muestreo, ahora es necesario calcula k de la frecuencia \frac{2π}{N} para ello

x_{p}[n] = x(nT)

x_{p}[n] = \frac{1}{2} (1 − cos(20πnT))

x_{p}[n] = \frac{1}{2} (1 − cos(\frac{2π}{N} kn))

Por comparación podemos ver que 20πT = \frac{2π}{N}k con lo cual tenemos

\frac{2π}{N} k = 20πT
k = \frac{20πT N}{2π}
k = \frac{20π \frac{5}{2N} N}{2π}
k = 25

Con esto tenemos que la transformada Discreta de Fourier es

X[k] = \frac{1}{2}δ(k) − \frac{1}{4} (δ(k − 25) + δ(k + 25)))

c) Si duplicamos el periodo de Muestreo T_{s} = \frac{5}{N} aplicando los pasos anteriores tenemos

\frac{2π}{N} k = 20πT
k = \frac{20πT N}{2π}
k = \frac{20π \frac{5}{N} N}{2π}
k = 50

y la TDF es

X[k] = \frac{1}{2}δ(k) − \frac{1}{4} (δ(k − 50) + δ(k + 50)))

Note como la frecuencia de muestreo afecta la transformada discreta de Fourier.

Ejemplo 3

Para el ejemplo anterior considerando un periodo de muestreo T_{s} = 0,0390625, calcular utilizando en código dado y los valores muestreados x_{p}[n], los valores de la función en t = T, t = 1,5T y t = 2,0.

Solución

Para llenar el arreglo dada la función x(t) = sen^{2}(10π) se hizo

x_{p}[n] = sen^{2}(10πnT)

Al aplicar

t calculados valores reales
0.0390625 0.9415440651830205 0.9415440651830208
0.05859375 0.964359871520714 0.9637760657954398
0.078125 0.6343932841636447 0.6343932841636455

Filtros.

7.1. Filtros Pasa bajas.

La idea de un filtro es permitir el paso de solamente un cierto ancho de banda de un espectro de frecuencias. En el caso de un filtro pasa bajas la frecuencias que este dejará pasar son aquellas que se encuentran cerca de la frecuencia cero o componente de C.D.

El filtro pasa bajas más simple es h(n) = [\frac{1}{2}, \frac{1}{2}] consideremos una señal dada por x(n) = [1, 1, 1, 2, 2, 2, 2, 1, 1, 1] al realizar la convolución de esta señal con el kernel tenemos

z_{1} (n) = x(n) ∗ h(n)

= [1, 1, 1, 2, 2, 2, 2, 1, 1, 1] ∗ [\frac{1}{2}, \frac{1}{2}]

= [1, 1, 1, 1,5, 2, 2, 1,5, 1, 1, 1]

si volvemos a convolucionar con el mismos kernel tenemos

z_{2} (n) = [1, 1, 1, 1,25, 1,75, 1,75, 1,25, 1, 1, 1]

el aplicar sucesivamente el filtrado pasa–bajas dará lugar a una señal plana. En lugar de convolucionar varias veces la señal x(n) con el kernel h(n) resulta mós practico hacer la convolución del kernel consigo mismo dando lugar a

1,1 1,2,1 1,3,3,1 1,4,6,4,1

Pero veamos la transformada de Fourier del kernel propuesto para ser filtro pasa–bajas.

H(k) = \frac{1}{N} ∑_{k=0}^{N} h(n)e^{−j(\frac{2π}{N})nk}

= \frac{1}{N} [1 × e^{−j(\frac{2π}{N})0k} + 1 × e^{−j(\frac{2π}{N})1k}]

= \frac{1}{N} [1 + e^{−j(\frac{2π}{N})k}]

la magnitud de este filtro es

|H(k)| = \frac{1}{N} \sqrt{(1 + cos \frac{2πk}{N} )^{2} + (sen\frac{2πk}{N} )^{2}}

= \frac{1}{N} \sqrt{1 + 2 cos \frac{2πk}{N} + cos^{2} \frac{2πk}{N} + sen^{2} \frac{2πk}{N}}

= \frac{1}{N} \sqrt{2 + 2 cos \frac{2πk}{N}}

= \frac{1}{N} cos \frac{πk}{N}

note que esta función tiene un móximo en cero y un mínimo en N, razón por la cual, deja pasar las bajas frecuencias y atenúna las altas frecuencias. En la figura podemos ver el comportamiento del filtro pasa bajas en el dominio de la frecuencia (ver figura 7.1).

Figura 7.1: Respuesta a la frecuencia del filtro de pasabajas

7.2. Filtros Pasa Altas.

El filtro pasa bajas más simple es h(n) = [−\frac{1}{2}, \frac{1}{2}] consideremos una señal dada por x(n) = [1, 1, 1, 2, 2, 2, 2, 1, 1, 1] al realizar la convolución de esta señal con el kernel tenemos

z_{1} (n) = x(n) ∗ h(n)

= [1, 1, 1, 2, 2, 2, 2, 1, 1, 1] ∗ [−\frac{1}{2}, \frac{1}{2}]

= [1, 0, 0, −0,5, 0, 0, 0,5, 0, 0, 0]

El aplicar sucesivamente el filtrado pasabajas dará lugar a una señal plana. En lugar de convolucionar varias veces la señal x(n) con el kernel h(n) resulta mós práctico hacer la convolución del kernel consigo mismo dando lugar a

+1,-1 +1,-2,+1 +1,-3,+3,-1 +1,-4,+6,-4,-1

Pero veamos la transformada de Fourier del kernel propuesto para ser filtro pasa altas.

H(k) = \frac{1}{N} ∑_{k=0}^{N} h(n)e^{−j(\frac{2π}{N})nk}

= \frac{1}{N} [1 × e^{−j(\frac{2π}{N})0k} − 1 × e^{−j(\frac{2π}{N})1k}]

= \frac{1}{N} [1 − e^{−j(\frac{2π}{N})k}]

la magnitud de este filtro es

|H(k)| = \frac{1}{N} \sqrt{(1 − cos \frac{2πk}{N} )^{2} + (sen\frac{2πk}{N} )^{2}}

= \frac{1}{N} \sqrt{1 − 2 cos \frac{2πk}{N} + cos^{2} \frac{2πk}{N} + sen^{2} \frac{2πk}{N}}

= \frac{1}{N} \sqrt{2 − 2 cos \frac{2πk}{N}}

= \frac{1}{N} sen\frac{πk}{N}

note que esta función tiene un móximo en \frac{N}{2} y un mínimo en 0, razón por la cual, deja pasar las altas frecuencias y atenña las bajas frecuencias. En la figura podemos ver el comportamiento del filtro pasa bajas en el dominio de la frecuencia (ver figura 7.2).

Figura 7.2: Respuesta a la frecuencia del filtro de pasa altas

7.3. Filtro pasa bajas Butterworth.

La magnitud al cuadrado de la respuesta de un filtro Butterworth esta dada por la expresión

|H (k)|^{2} = \frac{1}{1 + (\frac{2πk}{N})^{2n}}

Note que |H (0)|^{2} = 1 y que |H (\frac{N}{2π})|^{2} = 0,5, lo cual le da su característica de ser un filtro pasa bajas, pero adicionalmente este filtro tiene otro parámetro de control que es el exponente n al cual esta elevado la frecuencia. Si cambiamos este valor tendremos que la ventana del filtro se modifica tal como se observa en la figura 7.3.

Podemos notar que para valores bajos de n la ventana es suave y para valores grandes tiende a una ventana cuadrada.

Figura 7.3: Respuesta a la frecuencia del filtro de Butherword

7.4. Filtros de pasa banda.

La ecuación de una función gaussiana g la podemos representar por:

g (μ, σ) = \frac{1}{\sqrt{2π}σ} e^{− \frac{(x−μ)^{2}}{2σ^{2}}}

y podemos demostrar que la transformada de fourier es otra gaussiana G dada como

g (μ, σ_{t}) ↔^{F} G (μ, σ_{ω})

y que σ_{t}σ_{ω} ≥ K. Esta expresión es la formulación del principio de incertidumbre y el caso de una gaussiana se tiene la igualdad. El principio de incertidumbre establece que no podemos tener definición simultáneamente en el dominio de la frecuencia y del tiempo.

Esta función puede ser utilizada como filtro pasa bajas cuando μ = 0, pero con valores diferentes cambiaremos la banda de la señal que queremos filtrar. Así por ejemplo para filtrar una señal consideraremos que la gaussiana tiene una media cero y que aplicaremos una translación en el tiempo y/o en la frecuencia de la señal dada por μ

G (k) = \frac{1}{\sqrt{2π}σ} e^{− \frac{k^{2}}{2σ^{2}}}

si queremos hacer el filtrado en una frecuencia w_{k} = \frac{2π}{N}k simplemente situamos una gaussiana en G (k − μ) y otra en G (k + μ) la transformada de fourier de estas será

G (k) ↔^{F} g(n)

\frac{1}{2}G (k − μ) ↔^{F} \frac{1}{2}g(n)e^{j(\frac{2π}{N})μn}

\frac{1}{2}G (k + μ) ↔^{F} \frac{1}{2}g(n)e^{−j(\frac{2π}{N})μn}

el filtro en la frecuencia es la suma de G (k − μ) y G (k + μ) por lo que el filtro de convolución es la suma de estos dos

f_{1}(n) = \frac{1}{2}g(n) [cos (\frac{2π}{N} ) μn + jsen (\frac{2π}{N} ) μn]

f_{2}(n) = \frac{1}{2}g(n) [cos (\frac{2π}{N} ) μn − jsen (\frac{2π}{N} ) μn]

f_{1}(n) + f_{2}(n) = g(n) cos (\frac{2π}{N} ) μn

7.4.1. Ejemplo

Una señal x(t) = seno(10t) + N(0, 0,3), esta corrompida con ruido de Naturaleza gaussiana con media cero y varianza 0.3. Entonar calcular la media y la varianza de un filtro gaussiano pasa banda dada una frecuencia de muestreo T_{s} = 0,05 seg y un número total de muestras de N = 64

Podemos calcular el móltiplo de la frecuencia \frac{2π}{N} recordando que la frecuencia de la señal muestreada es ω_{s} = ωt. Con esto tenemos que

k = \frac{ωT N}{2π}
k = \frac{10 × 0,05 × 64}{2π}
k = 5,0929

Haciendo esto tenemos que nuestra Gaussiana la colocaremos en μ = 5,0929 y dado que no tenemos una frecuencia entera utilizaremos una varianza σ = 0,5. El la figura 7.4(a), se muestra la señal original con ruido y en la figuras 7.4(b) y 7.4(c) la parte real e imaginaria de su transformada Discreta de Fourier. El la figura 7.4(d) se muestra las Gaussiana entonadas para filtrar la señal y en la figura 7.4(e) el resultado de haber aplicado el filtro.

(a) Señal con ruido  (b) Parte real DFT

(c) Parte Imaginaria DFT  (d) Filtro entonado

(e) Señal filtrada

Figura 7.4: Ejemplo de filtrado de una señal con ruido

7.5. Filtro de Gabor

Pero que ocurre si en lugar de colocar un filtro en \frac{1}{2}G (k − μ) y otro en \frac{1}{2}G(k + μ), pusióramos uno solo en G (k − μ)

f(n) = g(n) [cos (\frac{2π}{N} ) μn + jsen (\frac{2π}{N} ) μn]

note que la parte real es equivalente a la encontrada anteriormente, pero además tenemos una parte compleja. El hecho de tener solamente un filtro G en el dominio de la frecuencia, nos da mós información acerca de la banda que andamos buscando

f_{R}(n) = g(n) cos (\frac{2π}{N} ) μn

f_{I}(n) = g(n)sen (\frac{2π}{N} ) μn

a este tipo de filtro se le conoce como filtro de Gabor. La magnitud del filtro nos daró los puntos de la señal donde la frecuencia para la cual entonamos el filtro produce esa señal, siendo mas alta en los puntos de dicha frecuencia.

7.5.1. Filtro de Gabor en dos dimensiones

Para el caso de dos dimensiones el Kernel de Gabor que utilizaremos es

f(n, m) = g(n, m)e^{j( \frac{2π}{N} nk_{0}+ \frac{2π}{M} ml_{0})}

La parte real del filtro es

f_{R}(n, m) = g(n, m) cos(\frac{2π}{N} nk_{0} + \frac{2π}{M} ml_{0})

y la parte imaginaria es

f_{I}(n, m) = g(n, m) sin(\frac{2π}{N} nk_{0} + \frac{2π}{M} ml_{0})

Dado que la convolución requiere de (NM)^{2} operaciones, resulta mas rápido implementar el kernel de Gabor de manera separable, así el número de operaciones seró N^{2}M + NM^{2}. De manera separable la parte real queda

f_{R}(n, m) = g(n, m) cos(\frac{2π}{N} nk_{0} + \frac{2π}{M} ml_{0})
f_{R}(n, m) = g(n) ∗ g(m) [cos(\frac{2π}{N} nk_{0}) cos(\frac{2π}{M} ml_{0})) − sin(\frac{2π}{N} nk_{0}) sin(\frac{2π}{M} ml_{0}))]
f_{R}(n, m) = [g(n) cos(\frac{2π}{N} nk_{0})] ∗ [g(m) cos(\frac{2π}{M} ml_{0}))] − [g(n) sin(\frac{2π}{N} nk_{0})] ∗ [g(m) sin(\frac{2π}{M} ml_{0}))]

Si definimos

g_{1}(n) = g(n) cos(\frac{2π}{N} nk_{0})

g_{2}(n) = g(n) sin(\frac{2π}{N} nk_{0})

g_{3}(m) = g(m) cos(\frac{2π}{M} ml_{0})

g_{4}(m) = g(m) sin(\frac{2π}{M} ml_{0})

podemos escribir la parte real como

f_{R}(n, m) = g_{1}(n) ∗ g_{3}(m) − g_{2}(n) ∗ g_{4}(m)

De forma similar podemos calcular la parte imaginaria y llegar al siguiente kernel separable.

f_{I}(n, m) = g_{1}(n) ∗ g_{4}(m) + g_{2}(n) ∗ g_{3}(m)

En la figura 7.5.1, podemos ver unos ejemplos del filtro de Gabor para diferentes valores de frecuencia espacial.

7.6. Como entonar un filtro

Resulta más natural al momento de entonar un filtro, hablar de frecuencias cuyas unidades son ciclos por segundo o Hertz. Dado el proceso de discretización, perdemos la noción del tiempo y nuestra frecuencia en el dominio discreto ahora esta dada en ciclos/muestra. La pregunta es como convertir Hertz a ciclos/muestra.

(a) Parte real, entonado en la frecuencia espacial (0,0)

(b) Parte imaginaria, entonado en la frecuencia espacial (0,0)

(c) Parte real, entonado en la frecuencia espacial (2,2)

(d) Parte imaginaria, entonado en la frecuencia espacial (2,2)

Figura 7.5: Filtros de Gabor

Consideremos una función muestreada con tama´ no N, en el dominio de Fourier en la posición k = N tenemos que N representa la frecuencia angular 2π, por lo tanto, podemos representar todas las frecuencia como móltiplos de \frac{2π}{N} como ω = \frac{2π}{N}k [rad/muestra] y la frecuencia es f_{k} = \frac{k}{N} [ciclos/muestra]. Dada una frecuencia de muestreo f_{m} en [muestras/seg] podemos calcular el tiempo necesario para una muestra, así : 1muestra = (1/f_{m}) seg, sustituyendo esto tenemos:

f_{t} = \frac{k}{N} [ciclos/ (1/f_{m}) seg]

f_{t} = \frac{f_{m}k}{N} [ciclos/seg]

f_{t} = f_{m}f_{k} [ciclos/seg]

finalmente k = \frac{f_{t}}{f_{m}} N.

7.6.1. Ejemplo

Consideremos una señal x(t) = sen(200πt)+ruido, la cual es muestreada con una f_{s} = 800 Hz, calcular el valor al cual debe ser entonado el filtro pasa bandas para un tama´ no de la muestra N = 128.

La frecuencia de la señal es f_{t} = 200π/(2π) = 100 hz, por lo tanto el valor al que debe entonarse es k = \frac{100}{800}128 = 16. De otra forma

x(t) = sen(200πt) = sen(200π \frac{n}{800}) = sen(\frac{2π}{128}16n)

7.7. Filtro Elimina Banda

Dada una señal x(n) el filtro pasa banda nos permitirá para la señal dada, solo mostrar la señal para la cual fue entonado el filtro. Pero como hacer si queremos mostrar la señal eliminando una banda dada.

En el domino de la frecuencia multiplicamos la transformada de Fourier de la señal por una gaussiana posicionada en k_{0}, esto es

X̂_{k_{0}}(k) = X(k)\frac{1}{\sqrt{2π}σ_{k}} e^{\frac{−(k−k_{0})^{2}}{2σ^{2}_{k}}}

Para eliminar la señal hacemos

X̂_{−k_{0}}(k) = X(k)\frac{1}{\sqrt{2π}σ_{k}} [1 − e^{\frac{−(k−k_{0})^{2}}{2σ^{2}_{k}}} ]

y aplicamos la transformada de Fourier

x̂_{−k_{0}}(n) = x(n) ∗ cte [δ (n) − g (n) e^{\frac{2πk_{0}}{N} n}]
= cte [x(n) − x(n) ∗ g (n) e^{\frac{2πk_{0}}{N} n}]
= x(n) − x̂_{k_{0}}(n)

De la expresión anterior podemos deducir que el filtro elimina banda consiste en restar a la señal original, la señal después de aplicar el filtro pasa banda.

7.8. Filtro de Mediana

Un filtro de mediana consiste en buscar en una vecindad de tama´ no N × M la mediana de un conjunto de pixeles y reemplazar el valor del pixel centrada por la mediana. Para calcular la mediana será necesario, al igual que la convolución, hacer un barrido sobre toda la imagen y ordenar los N × M pixeles de menor a mayor y la mediana seró el valor a la mitad del arreglo [Parker, 1997]. El filtro de mediana tiene la caracterástica de preservar los bordes.

En base a esto podemos proponer el filtro de moda, en el cual a diferencia de la mediana de la muestra se calcula la moda. Otras alternativas podría se remplazar la media o moda por los valores min o max de la ventana de convolución.

7.9. Filtro Binario

En el caso de tener una imagen binaria aplicar un filtro de convolución dar´¡ como resultado una imagen con números reales. En este caso la convolución la reemplazaremos por

y (n, m) = ⋁_{k=−∞}^{∞} ⋁_{l=−∞}^{∞} x (k, l) h (n − k, m − l) = x ∨ h

donde ⋁ representa la operación or y esta operación es conocida como la or convolución.

Por otro lado podemos calcular la and convolución como

y (n, m) = ⋀_{k=−∞}^{∞} ⋀_{l=−∞}^{∞} x (k, l) h (n − k, m − l) = x ∧ h

La or convolución también es conocida como operación de dilación y la and convolución como contracción. Si hacemos la aplicación de ambos (x ∨ h) ∧ h, tenemos primero una dilación y posteriormente una contracción, lo cual permite eliminar ruido de la imagen binaria. Se definen cerradura a la operación (x ∨ h) ∧ h y la apertura a la operación (x ∧ h) ∨ h

Para obtener los bordes de la imagen podemos hacer la diferencia de la imagen original con la imagen dilatada haciendo x − (x ∨ h) .

7.10. Filtro de Membrana

En esta sección describiremos la aplicación de la regularización al cálculo de un filtro de membrana. La regularización es una técnica de adaptable de filtrado de señales que permite estimar componentes de baja frecuencia. Abordaremos el problema desde un punto de vista estadístico, para ello, calcularemos primero el estimador de móxima verosimilitud, agregaremos información a priori y finalmente con regla de Bayes calcularemos el filtro.

7.10.1. Probabilidad de un evento

Si un experimento puede dar como resultado cualquiera N resultado diferentes igualmente probables, y si exactamente n de estos resultado corresponden al evento A, entonces la probabilidad del evento A es

P (A) = \frac{n}{N}

7.10.2. Probabilidad condicional

La probabilidad de que ocurra un evento B cuando se sabe que ha ocurrido algún otro evento A se denomina probabilidad condicional y se denota como P (B|A).

Consideremos el espacio muestral S = {1, 2, 3, 4, 5, 6, 7}, donde cada uno de los eventos son igualmente probables y estos eventos se agrupan en dos conjuntos A = {1, 2, 3, 4} y B = {3, 4, 5, 6, 7}, tal que S = A ∪ B. De acuerdo a la definición de probabilidad de un evento, la probabilidad de A es P (A) = 4/7 y la probabilidad del conjunto B es P (B) = 5/7. En la figura podemos ver, de manera mas clara, como se agrupan estos conjuntos (ver fig 7.6).

Figura 7.6: Ejemplo de dos conjuntos

Note que existe una intersección entre los conjuntos A y B, en notación de conjuntos este conjunto es A ∩ B y la probabilidad de estos elementos es 2/7.

Ahora queremos determinar la probabilidad de que ocurra un elemento del conjunto A dado que solo tengo solamente elementos del conjunto B, la cual denotaremos por P (A|B). Aplicando nuestra definición podemos ver que

P (A|B) = \frac{n \{A ∩ B\}}{n \{B\}}

si dividimos esta ecuación entre el total de elementos del espacio muestral tenemos

P (A|B) = \frac{n \{A ∩ B\} /N}{n \{B\} /N}

= \frac{P (A ∩ B)}{P (B)}

Para nuestros datos tenemos

P (A|B) = \frac{2}{5}

Por otro lado la probabilidad de P (B|A) = 1/2. Note que estos valores siguen cumpliendo nuestra definición de probabilidad muestral.

Dado lo anterior, podemos definir la probabilidad condicional de B dado A, que se denota P (B|A) com

P (A|B) = \frac{P (A ∩ B)}{P (B)}

7.10.3. Independencia

Dos eventos A y B son independientes si y solo si P (B|A) = P (B) y P (A|B) = P (A), tambión se cumple que P (A ∩ B) = P (A)P (B).

7.10.4. Regla de Bayes

Si los eventos B_{1}, B_{2}, ..., B_{k} constituyen una partición del espacio muestral S, donde P (B_{i}) ≠ 0, entonces, para cualquier evento A de S

P (B_{r}|A) = \frac{P (B_{r}) P (A|B_{r})}{∑_{i=1}^{k} P (B_{i}) P (A|B_{i})}

ver [Walpole and Myers, 1987]

7.10.5. Estimador de maxima verosimilitud.

Para encontrar el estimador de máxima verosimilitud f, dado un conjunto de observaciones g, suponemos que estos tienen una distribución de probabilidad normal con varianza σ^{2} = 1/2, la cual se expresa como:

p(g_{i}|f_{i}) = \frac{1}{C_{1}} e^{−(f_{i}−g_{i})^{2}}

considerando que los datos son generados de manera independiente tenemos que la distribución de probabilidades conjunta es

p(g|f) = ∏_{i=0}^{N−1} \frac{1}{C_{1}} e^{−(f_{i}−g_{i})^{2}}

Note que la probabilidad conjunta es máxima cuando g_{i} = f_{i}. Queremos calcular una señal suave, por lo que, nuestra información a priori será

p(f) = ∏_{i=0}^{N−1} \frac{1}{C_{2}} e^{−λ∇^{2}f_{i}} = ∏_{i=0}^{N−1} \frac{1}{C_{2}} e^{−λ(f_{i}−f_{i−1})^{2}}

aplicando la regla de Bayes, encontramos la expresión para la probabilidad a posteriori dada por:

p(f|g) = \frac{p(g|f) ∗ p(f)}{p(g)}  (7.1)

p(f|g) = ∏_{i=0}^{N−1} \frac{1}{C} e^{−[(f_{i}−g_{i})^{2}+λ(f_{i}−f_{i−1})^{2}]}

p(f|g) = ∏_{i=0}^{N−1} \frac{1}{C} e^{−U(f)}

La ecuación 7.1 es la distribución de probabilidad a posteriori, el mínimo de esta función, lo encontramos cuando la función de energía U(f) es mínima. Nuestro problema lo traducimos en calcular el mínimo de la siguiente función:

U(f) = ∑_{i=0}^{N−1} [(f_{i} − g_{i})^{2} + λ (f_{i} − f_{i−1})^{2}]

la cual representa la energía potencial almacenada en un sistema de resortes acoplados, tal como se muestra en la figura 7.7.

Para calcular el valor de f que minimiza la función U(f), hacemos \frac{∂U(f)}{∂f_{i}} = 0

\frac{∂U(f)}{∂f_{i}} = 2 (f_{i} − g_{i}) + 2λ (f_{i} − f_{i−1}) − 2λ (f_{i+1} − f_{i}) = 0  (7.2)

lo cual nos da el siguiente sistema lineal de ecuaciones :

[ \begin{matrix}1 + λ & −λ \\ −λ & 1 + 2λ & −λ \\ & & ... & ... & ... \\ & & −λ & 1 + 2λ & −λ \\ & & & −λ & 1 + 2λ & −λ \\ & & & & ... & ... & ... \\ & & & & & −λ & 1 + 2λ & −λ \\ & & & & & & −λ & 1 + λ \end{matrix} ] [ f_{1} ; f_{2} ; ⋮ ; f_{i} ; f_{i+1} ; ⋮ ; f_{N−1} ; f_{N} ] = [ g_{1} ; g_{2} ; ⋮ ; g_{i} ; g_{i+1} ; ⋮ ; g_{N−1} ; g_{N} ]

Figura 7.7: Sistema de resortes equivalente al filtro de membrana

7.10.6. Interpretación del filtro de Membrana en el dominio de la Frecuencia.

Tomando la transformada discreta de Fourier a ambos miembros de la ecuación 7.2 y agrupando los términos obtenemos:

F_{k} + λ (F_{k} − F_{k}e^{−(\frac{2π}{N})k}) − λ (F_{k}e^{(^{\frac{2π}{N}})k} − F_{k}) = G_{k}

agrupando tórminos obtenemos

[1 + 2λ − λe^{−(\frac{2π}{N})k} − λe^{(^{\frac{2π}{N}})k}] F_{k} = G_{k}

[1 + 2λ − 2λ cos (\frac{2π}{N} k)] F_{k} = G_{k}

finalmente tenemos que

F_{k} = \frac{1}{[1 + 2λ − 2λ cos (\frac{2π}{N}k)]}G_{k} = H_{k}G_{k}

H_{k} = \frac{1}{[1 + 2λ − 2λ cos (\frac{2π}{N}k)]}

donde la función H_{k} es un filtro pasa bajas equivalente al filtro de membrana. Note que cuando el valor de λ aumenta tenemos que el ancho de banda se disminuye, en las figuras 7.8 y 7.9 podemos ver la respuesta a la frecuencia de este filtro para valores de λ 1 y 1000 respectivamente.

Figura 7.8: Respuesta a la frecuencia para el filtro de membrana con λ = 1 Figura 7.9: Respuesta a la frecuencia para el filtro de membrana con λ = 1000

Aplicaciones

8.1. Detección de corrimientos utilizando técnicas de restauración

Dadas dos imágenes f′ y f, la delta de dirac δ y el operador de corrimiento S_{d} podemos decir:

f′ = S_{d}f = S_{d}f ∗ δ = f ∗ S_{d}δ = f ∗ δ_{d}

lo cual indica que si realizamos convolución con un impulso desplazado, tendremos una imagen desplazada a un nuevo punto de coordenadas dado por d. Aplicando transformada de Fourier podemos intentar calcular el valor de desplazamiento haciendo

f′(n) = f ∗ δ_{d}(n)
F′(k) = F (k)e^{j(\frac{2π}{N})kd}

de esta expresión podemos despejar el valor de d pero si un valor de F (k) = 0, no tenemos manera de estimar el valor del mismo.

En su lugar utilizaremos técnicas de Regularización, para ello, minimizamos la función

|f′ − f ∗ g|^{2}
sugeto a |l ∗ g|^{2}

donde f′(n) es la imagen desplazada, f(n) es la imagen origen, l(n) es un operador lineal que escogeremos de la mejor manera para facilitar la solución y g(n) es la función impulso que intentamos encontrar. Para nuestro calculo procedemos

U(g) = ∑_{n=0}^{N−1} [f′(n) − f ∗ g(n)]^{2} + β ∑_{n=0}^{N−1} |l ∗ g(n)|^{2}

= ∑_{n=0}^{N−1} [f′(n) − ∑_{i=0}^{N−1} f(n − i)g(i)]^{2} + β ∑_{n=0}^{N−1} [∑_{i=0}^{N−1} l(n − i)g(i)]^{2}

Para encontrar el valor de g que minimiza esta expresión hacemos

\frac{∂U(g)}{∂g(i)} = −2 ∑_{n=0}^{N−1} [f′(n) − ∑_{i=0}^{N−1} f(n − i) ∗ g(i)] f(n − i)

+2β ∑_{n=0}^{N−1} [∑_{i=0}^{N−1} l(n − i) ∗ g(i)] l(n − i)

= 0

reorganizado términos tenemos

∑_{n=0}^{N−1} [∑_{i=0}^{N−1} f(n − i) ∗ g(i)] f(n − i) + β ∑_{n=0}^{N−1} [∑_{i=0}^{N−1} l(n − i) ∗ g(i)] l(n − i) = ∑_{n=0}^{N−1} f′(n)f(n − i)

[f ∗ f^{r} + βl ∗ l^{r}] ∗ g(n)+ = f′ ∗ f^{r}(n)

Al aplicar transformada de Fourier tenemos:

G(k) = \frac{[F′(k)]^{∗} F (k)}{[F (k)]^{∗} F (k) + β [L(k)]^{∗} L(k)}

G(k) = \frac{[F′(k)]^{∗} F (k)}{|F (k)|^{2} + β |L(k)|^{2}}

al sacar la antitrasformada de Fourier calculamos la función impulso g(n).

8.1.1. Ejemplos

Se utilizaron como entrada la imagen f mostrada en 8.1(a), y la imagen f′ de la figura 8.1(b). Al parámetro β se le dio un valor de 0.1.

Después de calcular g(n) al computar la TFD^{−1} de G(k), se localizó una magnitud móxima igual a 1926.3965, en las coordenadas (2,179), lo que significa que el segmento de la imagen de f′ las coordenadas cuya columna es mayor que 179 no es visible desde f, y que los renglones de f′ menores que 2 no son visibles desde f. Una vez ensamblada la imagen panorámica con la información anterior, el resultado se muestra en la Fig. 8.2.

(a) Imagen f  (b) Imagen f′

Figura 8.1: Imógenes de entrada para el sistema.

Figura 8.2: Imagen panorómica que se formó con la imagen f y f′.

Pudiera parecer que el ensamble no fue correcto, puesto que las líneas rectas no parecen ajustar perfectamente. Sin embargo, este problema se debe a tres factores principales. El primero, es que la cámara que se utilizó cuenta con una distorsión considerable, que se puede observar en la curvatura de las lóneas que se suponen rectas. En segundo lugar, al tomar lafotografía pudo haber existido una transformación más alla de una traslación horizontal. Finalmente, la iluminación de la escena captada por la cámara puede ser diferente para cada fotografía lo cual perjudicaría el comportamiento del algoritmo A pesar de estos tres factores, se encuentra una aproximación bastante aceptable que puede servir incluso como proveedor de valores iniciales de corrección proyectiva. En el caso ideal, en donde de estos factores no se presentan, el resultado es también ideal. Para ejemplificar esta situación, podemos observar las imágenes iniciales en la Fig. 8.1.1, y el resultado se observa en la Fig. 8.4. En este caso, el empalme es perfecto puesto que no existe otra variante más que traslación, ni cambio de iluminación entre las imágenes, ni distorsión de la cámara. De hecho, este ejemplo se realizó fragmentando una imagen en tres más pequeñas con diferente traslación.

(a) Imagen 1  (b) Imagen 2  (c) Imagen 3

Figura 8.3: Imógenes de entrada para la creación de una panorómica final.

Figura 8.4: Imagen panorómica que se formó con la imagen 1, 2 y 3. 112

Tareas

9.1. Operadores de Ventana

Dada la imagen origen f (ver figura 9.1), establecer una ventana W y obtener la imagen destino g mediante las siguientes operaciones:

1. Calcular el valor máximo de la ventana en (i, j) y colocarlo en g_{ij}.

2. Sobre la ventana en (i,j) calcular la mediana y colocar el valor de la mediana en g_{ij}.

3. Para la ventana en (i,j) calcular g_{ij} = ∑_{k=i−n_{k}}^{i+n_{k}} ∑_{l=j−n_{k}}^{j+n_{k}} w_{kl}f_{kl} considerando:

w_{kl} = \begin{cases}−1, & si l < j \\ 0, & si l = j \\ 1, & si l > j \end{cases}

Desarrollo A continuación se muestra el resultado depués de aplicar las operaciones. Las figuras 9.2(a), 9.2(b) y 9.2(c) muestran las soluciones a cada uno de los incisos.

Figura 9.1: Imagen original.

(a) Inciso 1.

(b) Inciso 2.  (c) Inciso 3.

Figura 9.2: Operadores de ventana.

9.2. Transformaciones Afín y Proyectivas

Dada la transformación:

T = [ \begin{matrix}1.3 & −0.1 & 10 \\ −0.5 & 1 & 5 \\ 0 & 1 ∗ 10^{−4} & 1 \end{matrix} ]

calcular:

a) Una imagen destino barriendo todos los puntos de la imagen original.

b) Una imagen destino barriendo todos los puntos de la imagen destino y aproximando coordenadas al entero más próximo.

c) Igual que b) pero usando interpolación bilineal.

Figura 9.3: Imagen original.

Desarrollo En las figuras 9.4(a), 9.4(b) y 9.4(a) podemos observar los resultados de aplicar las transformaciones indicadas en cada uno de los incisos.

(a) Inciso a.

(b) Inciso b.  (c) Inciso c.

Figura 9.4: Aplicacion de la Transformación afín aplicando varios procedimientos

9.3. Cálculo de una Homografía

Dado el conjunto de puntos:

P_{1} = [ \begin{matrix}10, 10 \\ 10, 50 \\ 50, 50 \\ 50, 10 \end{matrix} ]

P_{2} = [ \begin{matrix}8.9108, 10.89108 \\ 4.95049, 46.53465 \\ 39.04761, 48.57142 \\ 42.85714, 14.28571 \end{matrix} ]

que se muestran en figura 9.5, realizar el cálculo de los parámetros de una homografía que representa un modelo proyectivo utilzando mínimos cuadrados y comparar con le m’etodo de descomposici’on en valores singulares (SVD). Desarrollo Del cálculo de las tranformaciones respectivas en x y y proyectadas en el plano z = 1, se obtienen las siguientes ecuaciones:

h_{1}x + h_{2}y + h_{3} − h_{7}x̂y − h_{8}x̂y = x̂  (9.1)

Figura 9.5: Conjuntos P_{1} y P_{2}.
h_{4}x + h_{5}y + h_{6} − h_{7}ŷx − h_{8}ŷy = ŷ  (9.2)

donde x̂, ŷ pertenecen al conjunto de puntos P_{2} y x, y pertenecen al conjunto de puntos P_{1}. Las cuales se pueden expresar en forma matricial de la siguiente manera:

[ \begin{matrix}x & y & 1 & 0 & 0 & 0 & −x̂x & −x̂y \\ 0 & 0 & 0 & x & y & 1 & −ŷx & −ŷy \end{matrix} ] [ \begin{matrix}h_{1} \\ h_{2} \\ h_{3} \\ h_{4} \\ h_{5} \\ h_{6} \\ h_{7} \\ h_{8} \end{matrix} ] = [ x̂ ; ŷ ]

como tenemos más incógnitas que ecuaciones podemos multiplicar ambos lados de la ecuacion por la transpuesta de la matriz de términos que multilpican a las h_{i}:

[ \begin{matrix}x & 0 \\ y & 0 \\ 1 & 0 \\ 0 & x \\ 0 & y \\ 0 & 1 \\ −x̂x & −ŷx \\ −x̂y & −ŷy \end{matrix} ] [ \begin{matrix}x & y & 1 & 0 & 0 & 0 & −x̂x & −x̂y \\ 0 & 0 & 0 & x & y & 1 & −ŷx & −ŷy \end{matrix} ] [ \begin{matrix}h_{1} \\ h_{2} \\ h_{3} \\ h_{4} \\ h_{5} \\ h_{6} \\ h_{7} \\ h_{8} \end{matrix} ] = [ \begin{matrix}x & 0 \\ y & 0 \\ 1 & 0 \\ 0 & x \\ 0 & y \\ 0 & 1 \\ −x̂x & −ŷx \\ −x̂y & −ŷy \end{matrix} ] [ x̂ ; ŷ ]

sin embargo, al despejar el vector de incógnitas se observa que la matriz resultante de multiplicar la matriz de coeficientes por su transpuesta se trata de una matriz singular sin importar cual sea el valor de x̂, ŷ, x o y, lo cual, no no permite calcular su inversa y por lo tanto no es posible resolver el sistema.

Para poder obtener los parámetros se observó que un par de puntos correspondientes generan dos ecuaciones del sistema a resolver, por lo que 4 pares de puntos correspondientes nos generarian igual número de ecuacions que incógnitas.

[ \begin{matrix}x_{1} & y_{1} & 1 & 0 & 0 & 0 & −x̂_{1}x_{1} & −x̂_{1}y_{1} \\ 0 & 0 & 0 & x_{1} & y_{1} & 1 & −ŷ_{1}x_{1} & −ŷ_{1}y_{1} \\ x_{2} & y_{2} & 1 & 0 & 0 & 0 & −x̂_{2}x_{2} & −x̂_{2}y_{2} \\ 0 & 0 & 0 & x_{2} & y_{2} & 1 & −ŷ_{2}x_{2} & −ŷ_{2}y_{2} \\ x_{3} & y_{3} & 1 & 0 & 0 & 0 & −x̂_{3}x_{3} & −x̂_{3}y_{3} \\ 0 & 0 & 0 & x_{3} & y_{3} & 1 & −ŷ_{3}x_{3} & −ŷ_{3}y_{3} \\ x_{4} & y_{4} & 1 & 0 & 0 & 0 & −x̂_{4}x_{4} & −x̂_{4}y_{4} \\ 0 & 0 & 0 & x_{4} & y_{4} & 1 & −ŷ_{4}x_{4} & −ŷ_{4}y_{4} \end{matrix} ] [ \begin{matrix}h_{1} \\ h_{2} \\ h_{3} \\ h_{4} \\ h_{5} \\ h_{6} \\ h_{7} \\ h_{8} \end{matrix} ] = [ x̂_{1} ; ŷ_{1} ; x̂_{2} ; ŷ_{2} ; x̂_{3} ; ŷ_{3} ; x̂_{4} ; ŷ_{4} ]

El vector solución se muestra a continuación:

[ \begin{matrix}h_{1} \\ h_{2} \\ h_{3} \\ h_{4} \\ h_{5} \\ h_{6} \\ h_{7} \\ h_{8} \end{matrix} ] = [ \begin{matrix}0.9000 \\ −0.1000 \\ 0.9999 \\ 0.1000 \\ 0.9000 \\ 1.0000 \\ 0.0010 \\ 0.0000 \end{matrix} ]

el cual es semejante al obtenido por el método de SVD:

[ \begin{matrix}h_{1} \\ h_{2} \\ h_{3} \\ h_{4} \\ h_{5} \\ h_{6} \\ h_{7} \\ h_{8} \end{matrix} ] = [ \begin{matrix}0.8999999999999954 \\ −0.09999999999999934 \\ 0.9999999999999948 \\ 0.09999999999999791 \\ 0.8999999999999946 \\ 1.0000000000001192 \\ 0.0010000000000000475 \\ 1.1034825056055406E − 16 \end{matrix} ]

9.4. Modelo de la Cámara Pin-Hole

Dada la geometría y el modelo de la c’amara Pin-Hole, figura 9.6, elaborar un demo en el que se muestre la ubicación y orientación del centro de la cámara en el mundo real y, de acuerdo a dichos parámetros, que muestre la respectiva vista de la cámara.

Figura 9.6: Geometría y diagrama de la cámara Pin-Hole.

Desarrollo Con el objeto de mostrar la correcta obtención del modelo de la camara, el demo sugiere un objeto en el mundo real con forma de caja, ver figura 9.7. Las lńeas de color rojo son el contorno de la cara frontal, mientras que las líneas punteadas son el contorno de la cara posterior. el origen de la flecha verde en la figura anterior muestra el centro de la cámara y la direcci’on de la misma muestra el rayo principal. El punto negro muestra nuestro punto de referencia en el mundo real.

En la siguiente figura se muestra la vista correspondiente a la ubicación y orientaci’on del centro de la cámara mostrada en la figura 9.7 la cual se encuentra ubicada al frente del objeto centrada respecto a la cara frontal y

Figura 9.7: Ejemplo de posicionamiento y orientaci’on de la c’amara.

sin rotación en algún eje.

Figura 9.8: Vista de la cámara.

A continiación se muestra una imagen del demo en el cual se pueden variar diversos parámetros referentes a la ubicación e internos. También se muestran las variaciones de las matrices P y M de acuerdo a la manipulación de los parámetros.

Figura 9.9: Demo Pin-Hole.

9.5. Modelo Proyectivo de la cámara a partir

de Imágenes de un modelo 3D

Dado un conjunto de puntos R_{i} ⟷ S_{i} del mundo real 3D a la imágen 2D R_{i} = [X_{i}, Y_{i}, Z_{i}, 1]^{T}

S_{i} = [x_{i}, y_{i}, w_{i}, 1]^{T}
S_{i} = P R_{i}

[ \begin{matrix}x_{i} \\ y_{i} \\ w_{i} \end{matrix} ] = [ \begin{matrix}P_{11} & P_{12} & P_{13} & P_{14} \\ P_{21} & P_{22} & P_{23} & P_{24} \\ P_{31} & P_{32} & P_{33} & P_{34} \\ P_{41} & P_{42} & P_{43} & P_{33} \end{matrix} ] [ \begin{matrix}X_{i} \\ Y_{i} \\ Z_{i} \\ 1 \end{matrix} ]

dado que si dos vectores son paralelos entre s’i, el producto vectorial es = 0:

S_{i} × P R_{i} = 0

Si suponemos,

[ \begin{matrix}P_{11} & P_{12} & P_{13} & P_{14} \\ P_{21} & P_{22} & P_{23} & P_{24} \\ P_{31} & P_{32} & P_{33} & P_{34} \end{matrix} ] = [ \begin{matrix}P 1 \\ P 2 \\ P 3 \end{matrix} ]

P R_{i} = [ \begin{matrix}P 1R_{i} \\ P 2R_{i} \\ P 3R_{i} \end{matrix} ]

S_{i} × P R_{i} = 0 ⟹ [ \begin{matrix}i & j & k \\ x_{i} & y_{i} & w_{i} \\ P 1R_{i} & P 2R_{i} & P 3R_{i} \end{matrix} ]

para cada punto R_{i} del mundo real mapeado en la imagen 2D se obtiene un par de ecuaciones:

y_{i}P 3R_{i} − wiP 2R_{i} = 0
−x_{i}P 3R_{i} + w_{i}P 1R_{r} = 0

por lo tanto necesitamos al menos seis puntos para tener un número igual de ecuaciones e incógnitas.

Pasando nuestro sistema a la forma:

[ \begin{matrix}. & . & . \\ . & . & . \\ . & . & . \\ 0^{T} & −wiR_{i} & y_{i}R_{i} \\ w_{i}R_{i} & 0^{T} & −x_{i}R_{i} \\ . & . & . \\ . & . & . \\ . & . & . \end{matrix} ] [ \begin{matrix}P 1 \\ P 2 \\ P 3 \end{matrix} ] = 0

tenemos un sistema lineal que resolvemos utilzando el m’etodo de descomposición en valores singulares (SVD):

[ \begin{matrix}P 1 \\ P 2 \\ P 3 \end{matrix} ] = [ \begin{matrix}P_{11} & P_{12} & P_{13} & P_{14} \\ P_{21} & P_{22} & P_{23} & P_{24} \\ P_{31} & P_{32} & P_{33} & P_{34} \end{matrix} ]

Dada la matriz P calculamos los parámetros de la c’amara,

P = KR[T | −C]

donde: C = posición del punto principal de la cámara en el mundo real, esto lo podemos obtener al encontrar el vector en la dirección del espacio nulo P C = 0.

R = rotación de la cámara en las tres dimensiones.

K = matriz de parámetros internos de la cámara:

K = [ \begin{matrix}α_{x} & s & P_{x} \\ & α_{y} & P_{y} \\ & & 1 \end{matrix} ]

M = KR ; R = R_{z}^{T} R_{y}^{T} R_{x}^{T}

Para calcular K dado M y dado que una matriz de rotación es ortogonal:

M = KR = R_{z}^{T} R_{y}^{T} R_{x}^{T}
MR_{x} = R_{z}^{T} R_{y}^{T} R_{x}^{T} R_{x}
MR_{x} = R_{z}^{T} R_{y}^{T} I
MR_{x} = M′
M′R_{y} = M′′
M′′R_{z} = K

Desarrollo

Para la obtención de los parámetros de la cámara se sacaron dos fotos tomadas desde diferentes ángulos, de acuerdo a la densidad de pixeles/cm^{2} de las imágenes, se sabe que cada una mide aproximadamente 106×141 cm.

Para el cálculo de los parámetros con la primer figura se tomaron los puntos que aparecen en verde y como centro de referencia el punto en color azul cuyas coordenadas est’an dadas en cent’imetros:

Figura 9.10: Puntos tomados de la Imagen (a).

R_{i} = [ \begin{matrix}14.4 & 0 & 0 & 1 \\ 0 & 14.4 & 0 & 1 \\ 14.4 & 14.4 & 0 & 1 \\ 14.4 & 14.4 & 14.4 & 1 \\ 0 & 0 & 0 & 1 \\ 1.8 & 0 & 12.6 & 1 \\ 14.4 & 0 & 14.4 & 1 \\ 7.2 & 7.2 & 0 & 1 \\ 14.4 & 7.2 & 7.2 & 1 \\ 7.2 & 0 & 7.2 & 1 \end{matrix} ] ; S_{i} = [ \begin{matrix}94.3 & 58.84 & 1 \\ 38.52 & 1.11 & 1 \\ 100.44 & 9.1 & 1 \\ 137.27 & 32.49 & 1 \\ 44.77 & 63.08 & 1 \\ 53.23 & 99.69 & 1 \\ 115.57 & 98.53 & 1 \\ 69.99 & 38.1 & 1 \\ 109.22 & 42.83 & 1 \\ 75.14 & 78.74 & 1 \end{matrix} ]

P = [ \begin{matrix}−0.0487473 & 0.0133769 & 0.0144258 & −0.578005 \\ 0.000857938 & 0.0487553 & −0.00771372 & −0.812739 \\ −0,0000531202 & 0.000196116 & 0.000303633 & −0.0128544 \end{matrix} ]

K = [ \begin{matrix}127.1457 & −4.0322 & 68.476 \\ 0.0 & 127.0717 & 49.1362 \\ 0.0 & 0.0 & 1.0 \end{matrix} ]

C = [ \begin{matrix}2.4616 \\ 21.2699 \\ 29.6608 \end{matrix} ]

Para el cálculo de los par’ametros con la primer figura se tomaron los puntos que aparecen en verde y como centro de referencia el punto en color azul cuyas coordenadas est’an dadas en cent’imetros:

Figura 9.11: Puntos tomados de la Imagen (a).

R_{i} = [ \begin{matrix}14.4 & 0 & 0 & 1 \\ 1.8 & 12.6 & 0 & 1 \\ 14.4 & 14.4 & 0 & 1 \\ 14.4 & 14.4 & 14.4 & 1 \\ 0 & 0 & 0 & 1 \\ 0 & 0 & 14.4 & 1 \\ 14.4 & 0 & 14.4 & 1 \\ 7.2 & 7.2 & 0 & 1 \\ 14.4 & 7.2 & 7.2 & 1 \\ 7.2 & 0 & 7.2 & 1 \end{matrix} ] ; S_{i} = [ \begin{matrix}65.19 & 39.3 & 1 \\ 45.79 & 5.15 & 1 \\ 99.69 & 23.35 & 1 \\ 99.55 & 87.49 & 1 \\ 19.83 & 25.33 & 1 \\ 10.16 & 81.14 & 1 \\ 61.52 & 85.58 & 1 \\ 56.37 & 23.99 & 1 \\ 78.25 & 57.93 & 1 \\ 41.49 & 57.29 & 1 \end{matrix} ]

P = [ \begin{matrix}−0.126195 & −0.0154735 & 0.0224553 & −0.598582 \\ −0.047378 & 0.0594941 & −0.100067 & −0.779932 \\ −0.000451098 & 0.000736939 & 0.000222169 & −0v0306284 \end{matrix} ]

K = [ \begin{matrix}−130.0555 & 0.8049 & 63.4632 \\ 0.0 & 130.1217 & 54.0046 \\ 0.0 & 0.0 & 1.0 \end{matrix} ]

C = [ \begin{matrix}−6.1702 \\ 33.2872 \\ 14.9179 \end{matrix} ]

Conculsiones

Durante los primeros experimentos, con pocos puntos (7 pts.), se observó que la exactitud en la obtenci’on de las medidas tanto en el mundo real tanto como en la imagen en 2D es indispensable para la correcta obtención de los parámetros internos y posici’on de la cámara. El centro de la imagen se aproxima al valor ideal y la posici’on de la c’amara es congruente con la imagen. En experimentos posteriores, tomando m’as puntos, el error no mejoró considerablemente, se observ’o que no importa cuantos puntos se tomen como referencia para al cálculo de los par’ametros sino que m’as bien es importante la exactitud en las mediciones, puesto que al tomar m’as puntos siempre se va tener un margen de error en cada uno de ellos al igual que si tomamos sólo unos pocos. Al final se tomaron 10 puntos para este ejercicio.

9.6. Convolución

Dada la expresi’on de la convoluci’on:

y(n) = ∑_{k=−∞}^{∞} ∑_{l=−∞}^{∞} x(k, l)h(n − k, m − l)

demuestre a trav’es de la convolución de dos o m’as im’agenes que esta es:

a) Conmutativa.

b) Asociativa.

c) Invariante a la traslaci’on.

d) Lineal.

Desarrollo

Dadas las imagenes: cada una de 201×201 pixeles se obtuvieron los re-

(a) Im  (b) K_{1}^{(c) K_{2}  (d) K_{3}}

sultados que se muestran en las siguientes secciones.

Comnutatividad

f_{1} ⊗ f_{2} = f_{2} ⊗ f_{1}

Los resultados son:

(e) Im ⊗K_{1}  (f) K_{1}⊗ Im

Figura 9.12: Conmutatividad.

Asociatividad

(f_{1} ⊗ f_{2}) ⊗ f_{3} = f_{1} ⊗ (f_{2} ⊗ f_{3})

Los resultados son:

(a) (Im⊗K_{1}) ⊗ K_{3} (b) Im⊗(K_{1} ⊗ K_{3})

Figura 9.13: Asiciatividad.

Invarriabilidad a traslaci’on

Hata ahora hemos convolucionado im’agenes considerando que estas se inician en el origen en ambas dimensiones. Para probar que la convoluci’on es invariate a la traslaci’on trasladaremos a la imagen Im 50 pixeles en ambas dimensiones obteniendo: al convlolucionar la imagen trasladada con K_{3}

Figura 9.14: Imagen trasladada.

obtenemos los siguientes resultados:

(a) Translate(Im)

⊗K_{1}

(b) Im ⊗K_{1}

Figura 9.15: Invarianza a traslaci’on.

Linealidad

y_{1}(n) = ∑_{k=−∞}^{∞} ∑_{l=−∞}^{∞} x_{1}(k, l)h(n − k, m − l)

y_{2}(n) = ∑_{k=−∞}^{∞} ∑_{l=−∞}^{∞} x_{2}(k, l)h(n − k, m − l)

∑_{k=−∞}^{∞} ∑_{l=−∞}^{∞}[x_{1}(k, l) + x_{2}(k, l)]h(n − k, m − l) = y_{1}(n) + y_{2}(n)

(a) Im+K_{3}  (b) Im⊗K_{1}  (c) K_{3} ⊗ K_{1}

(d) (Im+K_{3}) ⊗ K_{1} {(e) (Im⊗K_{1}) + (K_{3} ⊗ \atop K_{1})}

9.7. Algunos Kerneles Interesantes

a) Generar los kerneles Caja Binomial Gaussiano Primeras derivadas en X y Y Segundas derivadas

b) Aplicar cada kernel a una imagen

Desarrollo

La imagen a la cual le aplicaremos los kerneles generados se muestra a continuaci’on:

Figura 9.16: Imagen original .

Para todos los kerneles generados se tomaron los mismos par’ametros. ’Estos son de 31 × 31 pixeles y de un n’umero de tama no de caja igual a 15, coeficientes binomiale = 10 y sigma igual a 15 respectivamente.

Generaci’on de kerneles

Caja

(a) Ker- nel de caja. (b) Imagen resultante.

Binomial

(c) Kernel bino- mial.  (d) Imagen resultante.

Gaussiano

(e) Kernel gaus- siano  (f) Imagen resultante.

Primeras derivadas

(g) ; Kernel ; D_{x}  (h) Imagen resultante.

{(i) Ker- \atop nel D_{y}} (j) Imagen resultante.

Segundas derivadas

(k) ; Kernel ; D_{xx}  (l) Imagen resultante.

(m) ; Kernel ; D_{xy}  (n) Imagen resultante.

(ñ) ; Kernel ; D_{yy}  (o) Imagen resultante.

9.8. Cálculo de la Transformada Discreta de Fourier

a) Dada la señal

X[n] = sen (\frac{6π}{N} ) n

con N=10 calcular la TDF resolviendo el sistema Φa = X.

b) Dada la señal

X[n] = \begin{cases}0 & si n < n_{0} \\ 1 & si n_{0} ≤ x < n_{1} \\ 0 & si n ≥ n_{0} \end{cases}

c) Dado a de 2 calcular la anti-transformada discreta de fourier activando los coeficientes [0], [0, 1], [0, 1, 2], ..., [0 − (N − 1)].

Desarrollo Existe una familia de exponenciales complejas con periodo N dadas por:

Φ_{k}(n) = e^{jk(\frac{2π}{N})n} = cos (\frac{k2πn}{N}) + jsen^{(\underline{k2πn})}

Con esto podemos llegar a hacer la representaci’ on de una señal peri’odica como la combinaci’on de señales complejas tal que:

X[N] = ∑_{k=0}^{N−1} a_{k}Φ_{k}(n)

Para determinar los coeficientes de Fourier tenemos un sistema lineal Φa = X en forma matricial:

Resolviendo el sistema ΦA = X, encontrando la inversa de la matriz Φ calculamos a:

A = [ \begin{matrix}0,0000 + 0,0000i \\ 0,0000 + 0,0000i \\ 0,0000 + 0,0000i \\ 0,0000 + 0,5000i \\ 0,0000 − 0,0000i \\ 0,0000 − 0,0000i \\ −0,0000 − 0,0000i \\ −0,0000 − 0,5000i \\ −0,0000 + 0,0000i \\ 0,0000 + 0,0000i \end{matrix} ]

2.

Tomando N = 10, n_{0} = 4 y n_{1} = 6 obtenemos la señal:

Resolviendo el sistema ΦA = X:

A = [ \begin{matrix}0,2000 + 0,0000i \\ −0,1118 + 0,1539i \\ −0,0500 − 0,1539i \\ 0,1118 + 0,0363i \\ −0,0500 + 0,0363i \\ −0,0000 − 0,0000i \\ −0,0500 − 0,0363i \\ 0,1118 − 0,0363i \\ −0,0500 + 0,1539i \\ −0,1118 − 0,1539i \end{matrix} ]

3.

Dado que de 2 se obtuvo:

A = [ \begin{matrix}0,2000 + 0,0000i \\ −0,1118 + 0,1539i \\ −0,0500 − 0,1539i \\ 0,1118 + 0,0363i \\ −0,0500 + 0,0363i \\ −0,0000 − 0,0000i \\ −0,0500 − 0,0363i \\ 0,1118 − 0,0363i \\ −0,0500 + 0,1539i \\ −0,1118 − 0,1539i \end{matrix} ]

Procedemos a calcular la anti-transformada discreta de Fourier como: X = ΦA, graficando la suma de las componentes en orden ascendente se muestra a continuación la evolución de la señal resulante pero primeramente la imagen original.

Figura 9.17: Suma de las componentes activando los coeficientes del 0 al 4.

9.9. Propiedades de la Transformada Discreta de Fourier

a) Demostrar las propiedades de la transformada de Fourer mediante el uso de im ágenes.

b) Demostrar la transformada de Fourier en dos dimensiones para las señales

Impulso unitario Constante

Propiedades

Linealidad

Esta propiedad nos dice que es equivalente multiplicar dos señales por determinado factor cada una y después sumarlas a multiplicar las transfor-Figura 9.18: Suma de las componentes activando los coeficientes del 5 al 9 a parte de los anteriormente activados.

madas de Fourier por los mismos factores respectivamente y llevar a cabo la suma en el dominio de la frecuencia, es decir:

ax(n) + by(n) ⟶ aX(k) + bY (k)

En las im ágenes se muestra primero las imágenes originales y la compa-

Figura 9.19: Im ágenes originales.

ración de las operaciones en ambos dominios. Las imagenes x y y fueron multiplicadas respectivamente por los factores a = 2,0 y b = 5,0.

(a) Suma en el tiempo.

(b) Transformada inversa de la suma en frecuencia.

Figura 9.20: Comparación de resultados en tiempo y frecuencia.

Desplazamiento en tiempo

Esta propiedad nos nos dice que desplazar una imagen en el dominio del tiempo es equivalente a multiplicar su transformada de Fourier por la exponencial compleja que correspondiente a n_{0} y m_{0}, es decir:

x(n − n_{0}, m − m_{0}) ⟶ e^{−j(\frac{2π}{N})(n_{0}k+m_{0}l)}X(k)

En las im ágenes anteriores se muestra la equivalencia de las operaciones en los dos dominios con un desplazamiento en el tiempo de n_{0} = 50 y m_{0} = 80.

Desplazamiento en frecuencia

Desplazar una señal, en este caso una imagen, en el dominio de la frecuencia equivale a multiplicar la imagen original en el dominio del tiempo por la exponencial compleja correspondiente a n_{0} y m_{0}, es decir:

e^{j(\frac{2π}{N})(n_{0}k+m_{0}l)}x(n) ⟶ X(k − k_{0}, l − l_{0})

En las imágenes se muestra la equivalencia de estas operaciones en ambos dominios con un desplazamiento en frecuencia de n_{0} = 2 y m_{0} = 5.

Simetría de señal real y par

Esta propiedad nos dice que la transformada de Fourier de una señal real y par es una señal con parte real únicamente y que también es par, es decir:

x(n, m) real y par ⟶ X(k, l) real y par

(a) Imagen original.

(b) Imagen desplazada en tiempo.

(c) Transformada inversa de operaci ón en frec .

Figura 9.21: Comparación de resultados en tiempo y frecuencia.

Con el objeto de visualizar mejor lo que se explica acerca de esta propiedad, en las imágenes se muestra ambas partes (real e imaginaria) al aplicar transformada inversa y con lo cual se observa que sólo hay la parte real.

(b) Operacion en dominio del tiempo.

(c) Transformada inversa de despl en frec.

Figura 9.22: Comparacion de resultados en tiempo y frecuencia.

Simetría de señal real e impar

Esta propiedad nos dice que la transformada de Fourier de una señal real e impar es una señal con parte imaginaria únicamente y que es impar, es decir:

x(n) real y par ⟶ X(k) real y par

Para visualizar mejor lo que se explica acerca de esta propiedad, en las imágenes se muestra ambas partes (real e imaginaria) al aplicar transformada inversa y con lo cual se observa que sólo hay parte imaginaria.

(a) Imagen real-par.

(b) Parte real de la transformada inversa.

(c) Parte imag. de la transformada inversa.

Figura 9.23: Operadores de ventana.

Parseval

Esta propiedad nos habla de que la transformada de Fourier es una operación conservativa ya que al aplicarla, la magnitud del resultado es igual a la magnitud de la señal original.

∑_{n=0}^{N−1} ∑_{m=0}^{M−1} | x(n, m) |⟶ NM ∑_{k=0}^{N−1} ∑_{l=0}^{M−1} | X(k, l) |

Al obtener la magnitud de ambas imágenes se obtuvieron los siguientes resultados:

Magnitud(Imagen) = 1.875586515 × 10^{9}
Magnitud(T F (Imagen)) = 1.8755865149999826 × 10^{9}

con lo cual se comprueba la relaci ón de Parseval.

(b) Parte real de la transformada inversa.

(c) Parte imag de la transformada inversa.

Figura 9.24: Operadores de ventana.

(a) Imagen original.

(b) Magnitud de la transformada de Fourier.

Figura 9.25: Imagen en tiempo y frecuencia.

Demostraciones

La transformada de Fourer en dos dimensiones se expresa:

X(k, l) = \frac{1}{NM} ∑_{n=0}^{N−1} ∑_{m=0}^{M−1} x(n, m)e^{−j(\frac{2π}{NM})(nk+ml)}

Impulso un itario

δ(n − n_{o}, m − m_{o}) = e^{−j(\frac{2π}{N} kn_{0}+ \frac{2π}{M} lm_{0})}

δ(n, m) = \begin{cases}1 & si n = 0 y m = 0 \\ 0 & en caso contrario \end{cases}

La transformada de Fourier, de acuerdo a su expresión, la calculamos:

X(k, l) = \frac{1}{NM} ∑_{n=0}^{N−1} ∑_{m=0}^{M−1} δ(n − n_{0}, m − m_{0})e^{−j(\frac{2π}{NM})(nk+ml)}

Tenemos entonces que para cualquier valor de m ≠ m_{0} y n ≠ n_{0} la funci ón impulso unitario δ(n, m) = 0. Por lo tanto la sumatoria se reduce a:

X(k, l) = \frac{1}{NM} e^{−j(\frac{2π}{NM})(n_{0}k+m_{0}l)} = \frac{1}{NM} e^{−j(\frac{2π}{NM} n_{0}k+ \frac{2π}{NM} m_{0}l)}

Constante

Consideramos:

x(n, m) = C

La transformada de Fourier, de acuerdo a su expresi ón, es:

X(k, l) = \frac{1}{NM} ∑_{n=0}^{N−1} ∑_{m=0}^{M−1} Ce^{−j(\frac{2π}{NM})(nk+ml)}

X(k, l) = \frac{C}{NM} ∑_{n=0}^{N−1} ∑_{m=0}^{M−1} e^{−j(\frac{2π}{NM})(nk+ml)} = \frac{C}{NM} ∑_{n=0}^{N−1} [∑_{m=0}^{M−1} e^{−j(\frac{2π}{NM})(ml)}] e^{−j(\frac{2π}{NM})(nk)}

La sumatoria interna equivale a M para l = 0 y 0 para l ≠ 0 y la externa eqivale a MN para k = 0 y 0 para k ≠ 0 por lo cual:

F F T (C) = Cδ(k, l)

Conculsiones Con el desarrollo de esta tarea se logró una mejor comprensión de la transformada de Fourier la cual tiene aplicaciones importantes en el tratamiento deseñales y diseño de filtros digitales gracias a que, al aplicarla, nos permite representar una señal en el dominio de la frecuencia y, a su vez, al aplicar la transformada inversa, regresar a la señal original en el dominio del tiempo. Principalmente se observó que gracias a sus propiedades es posible realizar ciertas operaciones en el dominio del tiempo y obtener los mismos resultados al hacerlas en el dominio de la frecuencia, algo que podemos aprovechar dependiendo de las características del problema que queramos resolver ya que podemos deducir en que dominio conviene más realizar la operación.

9.10. Transformada Discreta de Fourier de una señal muestreada

Dada la señal x(t) = cos(300t) + sen(450t):

Figura 9.26: x(t) = cos(300t) + sen(450t)

a) Discretizar la señal con un periodo de muestreo w_{0} = 200\frac{rad}{seg}.

b) Graficar la señal discretizada.

c) Obtener la TDF y la TDF inversa de la señal discretizada y graficar los resultados.

Desarrollo

Dado que: w_{0} = 200\frac{rad}{seg}, f_{0} = \frac{200}{2π}Hz ∴ T_{0} = \frac{2π}{200}segs. La señal discretizada será:

x[n] = x(nT_{0}) = cos [300 (\frac{2π}{200}) n] + sen [450 (\frac{2π}{200}) n]

x[n] = cos (3πn) + sen (4.5πn)

como sabemos, la frecuencia más grande que podemos ver en nuestro espectro relativamente a nuestro número de muestras es π, la cual equivale a \frac{N}{2} , ∴

x[n] = cos (3\frac{N}{2} n) + sen (\frac{9N}{4} n)

cuya transformada de Fourier es:

X(k) = \frac{1}{2} [δ (k − \frac{3N}{2} ) + δ (k + \frac{3N}{2} )]− \frac{j}{2} [δ (k − \frac{9N}{4} ) − δ (k + \frac{9N}{4} )]

Lo cual por considerar nuestra señal periódica queda:

X(k) = \frac{1}{2} [δ (k − \frac{N}{2} ) + δ (k + \frac{N}{2} )] − \frac{j}{2} [δ (k − \frac{N}{4} ) − δ (k + \frac{N}{4} )]

A continuación se muestra un acercamiento de la imagen original y la imagen muestreada, se ve claramente una gran diferencia.

Figura 9.27: Comparación de señal original y señal muestreada.

Tomando en consideración un número de muestras N = 128, a continuación se muestra la transformada de Fourier graficada de [0 − 2π], se observa efectivamente las frecuencias calculadas con anterioridad y la transformada inversa la cual no se trata de la señal original (debido a la frecuencia de muestreo).

Figura 9.28: Transformada de Fourier y señal reconstruida con TDFI.

9.11. Filtros

a) Dada la señal:

I(n, m) = cos (\frac{2π}{N} 20n) + N(0, σ)

N(0, σ) = ruido gaussiano, con media μ = 0 y varianza σ = 0.1, multiplicar la señal en el dominio de la frecuencia por una caja centrada en la frecuencia correspondiente para eliminar el ruido.

b) Para la misma señal I(n, m) calcular un kernel de convolución separable, multiplicarlo por una exponencial compleja e^{j(θk_{0}n+θl_{0}m)} para hacer un filtrado del ruido en la señal en el dominio del tiempo.

Desarrollo

Un filtro permite el paso solamente de ciertas frecuencias. Un filtro pasa bajas se caracteriza por permitir pasar las frecuencias cercanas a la componente de directa y atenuar las frecuencias más altas. Por otra parte un filtro

Figura 9.29: I(n, m).

pasa altas, es un tipo de filtro en cuya respuesta en frecuencia se atenúan las componentes de baja frecuencia pero no las de alta frecuencia. Un filtro paso banda es un tipo de filtro electrónico que deja pasar un determinado rango de frecuencias de una señal y atenúa el paso del resto.

1.

La transformada de Fourier de la señal I(n, m) se muestra en la figura 9.30, para hacer el filtrado en el dominio de la frecuencia es necesario hacer una multiplicación punto a punto de la trasformada de Fourier por nuestra caja (mostrada en la figura 9.31), esto es:

[Re(k, l) + Im(k, l)] · Caja(k, l)

la cual, debemos ubicar o centrar en la posición de la frecuencia principal que queremos dejar pasar, en este caso se requiere que solo pasen las frecuencias correspondientes a (k_{0}, l_{0}) = (20, 0) por lo que se centrará la caja en estas coordenadas. Se debe considerar también que el ancho de banda de nuestro filtro, es decir, las frecuencias que éste permitira pasar, depende del tamaño de nuestra caja. El resultado del filtrado se muestra en la figura 9.32.

Figura 9.30: Transformada de Fourier de I(m,n).

Figura 9.31: Caja centrada en (k_{0}, l_{0}) = (20, 0).

Figura 9.32: Resultado del filtrado de I(n, m).

2.

Para implementar un filtro que brinde un resultado semejante al ejercicio anterior pero mediante convolución en el dominio del tiempo y a un bajo costo computacional relativo, es conveniente tener kerneles de convolución separables y aplicar el desplazamiento en frecuencia necesario, multiplicando por la exponencial compleja correspondiente, en este caso (k_{0}, l_{0}) = (20, 0), obteniendo:

Figura 9.33: Kerneles de caja separables desplazados en fercuencia.

desarrollando la operación, esto es:

Caja(n, m) → [caja(n)caja(m)]e^{j(\frac{2π}{N} k_{0}n+ \frac{2pi}{M} l_{0}m)}

Considerando θ = \frac{2π}{N}, tenemos entonces:

Caja = g(n)g(m)e^{jθk_{0}n}e^{jθl_{0}m}
= [g(n)e^{jθk_{0}n}]g(m)e^{jθl_{0}m}
= g(n)[cos(θk_{0}n) + jsen(θk_{0}n)]g(m)[cos(θl_{0}m) + jsen(θl_{0}m)]

y separando la parte real y la imaginaria:

Real(n, m) = g(n)cos(θk_{0}n)g(m)cos(θl_{0}m) − g(n)sen(θk_{0}n)g(m)sen(θl_{0}m)

Imag(n, m) = g(n)sen(θk_{0}n)g(m)cos(θl_{0}m) + g(n)cos(θk_{0}n)g(m)sen(θl_{0}m)

cuyos resultados se pueden apreciar en la figura 9.34, la cual muestra el senal original filtrada.

Figura 9.34: Kerneles de caja separables desplazados en fercuencia.

Conculsiones

Para la realización de esta tarea se aprovecharon las propiedades de la convlución para la implementación de dos filtros, uno en el dominio del tiempo y el segundo en el dominio de la frecuencia, obteniéndose resultados similares en ambos casos al aplicarlos a una señal ruidosa para limpiarla. Lo cual comprobó la posiblidad de implementar filtros tanto en el espacio geométrico como en el espacio de la frecuencia

Bibliografía

[Jahne, 2002] Jahne, B. (2002). Digital Image Processing. Springer, New York, N.Y.

[Oppenheim, 1997] Oppenheim, A. V. (1997). Señales y Sistemas. Pearson Education, México.

[Parker, 1997] Parker, J. R. (1997). Algorithms for Image Processing and Computer Vision. Wiley, USA.

[Walpole and Myers, 1987] Walpole, R. E. and Myers, R. H. (1987). Probabilidad y Estadistica para Ingenieros. Interamericana, México D.F.