% Activate the following line by filling in the right side. If for example the name of the root file is Main.tex, write
% "...root = Main.tex" if the chapter file is in the same directory, and "...root = ../Main.tex" if the chapter is in a subdirectory.
 
%!TEX root =  ../vision.tex

\chapter{Filtros.}

\section{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

\begin{eqnarray*}
z_{1}\left( n\right) =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]
\end{eqnarray*}

si volvemos a convolucionar con el mismos kernel tenemos

\[
z_{2}\left( n\right) =[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

\begin{center}
1,1

1,2,1

1,3,3,1

1,4,6,4,1
\end{center}

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

\begin{eqnarray*}
H(k)=\sum_{k=0}^{N}h(n)e^{-j\left( \frac{2\pi }{N}\right) nk} \\
=\left[ 1\times e^{-j\left( \frac{2\pi }{N}\right) 0k}+1\times
e^{-j\left( \frac{2\pi }{N}\right) 1k}\right] \\
=\left[ 1+e^{-j\left( \frac{2\pi }{N}\right) k}\right]
\end{eqnarray*}

la magnitud de este filtro es

\begin{eqnarray*}
\left| H(k)\right| =\sqrt{\left( 1+\cos \frac{2\pi k}{N}\right)
^{2}+\left( sen\frac{2\pi k}{N}\right) ^{2}} \\
=\sqrt{1+2\cos \frac{2\pi k}{N}+\cos ^{2}\frac{2\pi k}{N}+sen^{2}%
\frac{2\pi k}{N}} \\
=\sqrt{2+2\cos \frac{2\pi k}{N}} \\
=2 \cos \frac{\pi k}{N}
\end{eqnarray*}

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úa las altas frecuencias. En la figura podemos ver el comportamiento del filtro pasa bajas en el dominio de la frecuencia (ver figura \ref{fig:pasabajas}).

\begin{figure}[h!]
{\normalsize \centering
\centerline{\psfig{figure = Filtros/figuras/pasabajas.eps,width=8.0cm}} }
\caption{Respuesta a la frecuencia del filtro de pasabajas}
\label{fig:pasabajas}
\end{figure}

\section{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

\begin{eqnarray*}
z_{1}\left( n\right) =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]
\end{eqnarray*}

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

\begin{center}
+1,-1

+1,-2,+1

+1,-3,+3,-1

+1,-4,+6,-4,-1
\end{center}

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

\begin{eqnarray*}
H(k)=\sum_{k=0}^{N}h(n)e^{-j\left( \frac{2\pi }{N}\right) nk} \\
=\left[ 1\times e^{-j\left( \frac{2\pi }{N}\right) 0k}-1\times
e^{-j\left( \frac{2\pi }{N}\right) 1k}\right] \\
=\left[ 1-e^{-j\left( \frac{2\pi }{N}\right) k}\right]
\end{eqnarray*}

la magnitud de este filtro es

\begin{eqnarray*}
\left| H(k)\right| =\sqrt{\left( 1-\cos \frac{2\pi k}{N}\right)
^{2}+\left( sen\frac{2\pi k}{N}\right) ^{2}} \\
=\sqrt{1-2\cos \frac{2\pi k}{N}+\cos ^{2}\frac{2\pi k}{N}+sen^{2}\frac{2\pi k}{N}} \\
=\sqrt{2-2\cos \frac{2\pi k}{N}} \\
=2 sen\frac{\pi k}{N}
\end{eqnarray*}

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 \ref{fig:pasaaltas}).

\begin{figure}[h!]
{\normalsize \centering
\centerline{\psfig{figure = Filtros/figuras/pasaaltas.eps,width=8.0cm}} }
\caption{Respuesta a la frecuencia del filtro de pasa altas}
\label{fig:pasaaltas}
\end{figure}

\section{Filtro pasa bajas Butterworth.}

La magnitud al cuadrado de la respuesta de un filtro Butterworth esta dada por la expresión

\[
\left| H\left( k\right) \right| ^{2}=\frac{1}{1+\left( \frac{2\pi k}{N}%
\right) ^{2n}}
\]

Note que $\left| H\left( 0\right) \right| ^{2}=1$ y que $\left| H\left(\frac{N}{2\pi }\right) \right| ^{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 \ref{fig:buther}.

\begin{figure}[h!]
{\normalsize \centering
\centerline{\psfig{figure = Filtros/figuras/buther.eps,width=8.0cm}} }
\caption{Respuesta a la frecuencia del filtro de Butherword}
\label{fig:buther}
\end{figure}

Podemos notar que para valores bajos de $n$ la ventana es suave y para valores grandes tiende a una ventana cuadrada.

\section{Filtros de pasa banda.}

La ecuación de una función Gaussiana $g$ la podemos representar por:

\[
g\left( \mu ,\sigma \right) =\frac{1}{\sqrt{2\pi }\sigma }e^{-\frac{\left(
x-\mu \right) ^{2}}{2\sigma ^{2}}}
\]

y podemos demostrar que la transformada de Fourier es otra Gaussiana $G$ dada como

\[
g\left( \mu ,\sigma _{t}\right) \stackrel{\mathcal{F}}{\leftrightarrow }
G\left( \mu ,\sigma _{\omega }\right)
\]

y que $\sigma _{t}\sigma _{\omega }\geq 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 $\mu =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 $\mu $

\[
G\left( k\right) =\frac{1}{\sqrt{2\pi }\sigma }e^{-\frac{k^{2}}{2\sigma ^{2}}%
}
\]

si queremos hacer el filtrado en una frecuencia $w_{k}=\frac{2\pi }{N}k$ simplemente situamos una Gaussiana en $G\left( k-\mu \right) $ y otra en $G\left( k+\mu \right) $ la transformada de Fourier de estas será

\begin{eqnarray*}
G\left( k\right) \stackrel{\mathcal{F}}{\leftrightarrow }g(n) \\
\frac{1}{2}G\left( k-\mu \right) \stackrel{\mathcal{F}}{\leftrightarrow }%
\frac{1}{2}g(n)e^{j\left( \frac{2\pi }{N} \mu n\right) } \\
\frac{1}{2}G\left( k+\mu \right) \stackrel{\mathcal{F}}{\leftrightarrow }%
\frac{1}{2}g(n)e^{-j\left( \frac{2\pi }{N} \mu n \right)}
\end{eqnarray*}

el filtro en la frecuencia es la suma de $G\left( k-\mu \right) $ y $G\left(k+\mu \right) $ por lo que el filtro de convolución es la suma de estos dos

\begin{eqnarray*}
h_{1}(n)=\frac{1}{2}g(n)\left[ \cos \left( \frac{2\pi }{N} \mu
n \right) +jsen\left( \frac{2\pi }{N}  \mu n \right)\right] \\
h_{2}(n)=\frac{1}{2}g(n)\left[ \cos \left( \frac{2\pi }{N}\mu
n \right) -jsen\left( \frac{2\pi }{N} \mu n \right) \right] \\
h(n) = h_{1}(n)+h_{2}(n)=g(n)\cos \left( \frac{2\pi }{N} \mu n\right)
\end{eqnarray*}

\subsection{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 \pi}{N}$ recordando que la frecuencia de la señal muestreada es $\omega_s = \omega t$. Con esto tenemos que

\begin{eqnarray*}
k &=& \frac{\omega T N}{2 \pi} \\
k &=& \frac{10 \times 0.05 \times 64}{2 \pi} \\
k &=& 5.0929
\end{eqnarray*}

Haciendo esto tenemos que nuestra Gaussiana la colocaremos en $\mu = 5.0929$ y dado que no tenemos una frecuencia entera utilizaremos una varianza $\sigma = 0.5$. El la figura \ref{fig:filtro_orig}, se muestra la señal original con ruido y en la figuras \ref{fig:filtro_real} y \ref{fig:filtro_imag} la parte real e imaginaria de su transformada Discreta de Fourier. El la figura \ref{fig:filtro_gauss} se muestra las Gaussiana entonadas para filtrar la señal y en la figura \ref{fig:filtro_sal} el resultado de haber aplicado el filtro.

\begin{figure}[h!]
\begin{center}
{\subfigure[Señal con ruido \label{fig:filtro_orig}]{\psfig{figure = Filtros/figuras/filtro_orig.eps,width=5.0cm}}
\subfigure[Parte real DFT \label{fig:filtro_real}]{\psfig{figure = Filtros/figuras/filtro_real.eps,width=5.0cm}} 
\subfigure[Parte Imaginaria DFT \label{fig:filtro_imag}]{\psfig{figure = Filtros/figuras/filtro_imag.eps,width=5.0cm}}
\subfigure[Filtro entonado \label{fig:filtro_gauss}]{\psfig{figure = Filtros/figuras/filtro_gauss.eps,width=5.0cm}}
\subfigure[Señal filtrada \label{fig:filtro_sal}]{\psfig{figure = Filtros/figuras/filtro_sal.eps,width=5.0cm}}\label{fig:filtro} }
\end{center}
\caption{Ejemplo de filtrado de una señal con ruido}
\end{figure}


\section{Filtro de Gabor}

Pero que ocurre si en lugar de colocar un filtro en $\frac{1}{2}G\left(k-\mu \right) $ y otro en $\frac{1}{2}G(k+\mu )$, pusiéramos uno solo en $G\left( k-\mu \right) $

\[
h(n)=g(n)\left[ \cos \left( \frac{2\pi }{N} \mu n \right) +jsen\left( \frac{%
2\pi }{N}  \mu n \right) \right]
\]

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

\begin{eqnarray*}
h_{R}(n)=g(n)\cos \left( \frac{2\pi }{N} \mu n \right) \\
h_{I}(n)=g(n)sen\left( \frac{2\pi }{N}  \mu n \right)
\end{eqnarray*}

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.

\subsection{Filtro de Gabor en dos dimensiones}

Para el caso de dos dimensiones el Kernel de Gabor que utilizaremos es

\[
h(n,m)=g(n,m)e^{j \left(\frac{2\pi }{N}nk_{0}+\frac{2\pi }{M}ml_{0} \right)}
\]
La parte real del filtro es
\[
h_{R}(n,m)=g(n,m)\cos \left(\frac{2\pi }{N}nk_{0}+\frac{2\pi }{M}ml_{0}\right)
\]
y la parte imaginaria es
\[
h_{I}(n,m)=g(n,m)\sin \left(\frac{2\pi }{N}nk_{0}+\frac{2\pi }{M}ml_{0}\right)
\]

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 y aplicando identidades trigonométricas podemos escribir:

\begin{eqnarray*}
h_{R}(n,m) &=&g(n)*g(m)\left[ \cos \left(\frac{2\pi }{N}nk_{0} \right)\cos \left(\frac{2\pi }{M}ml_{0} \right)-\sin \left(\frac{2\pi }{N}nk_{0} \right)\sin \left(\frac{2\pi }{M}ml_{0}\right)\right]
\\
h_{R}(n,m) &=&\left[ g(n)\cos \left(\frac{2\pi }{N}nk_{0} \right)\right] *\left[
g(m)\cos \left(\frac{2\pi }{M}ml_{0} \right)\right] -\left[ g(n)\sin \left(\frac{2\pi }{N}nk_{0} \right)\right] *\left[ g(m)\sin \left(\frac{2\pi }{M}ml_{0} \right)\right]
\end{eqnarray*}

Si definimos

\begin{eqnarray*}
g_{1}(n)=g(n)\cos \left(\frac{2\pi }{N}nk_{0} \right) \\
g_{2}(n)=g(n)\sin \left(\frac{2\pi }{N}nk_{0} \right) \\
g_{3}(m)=g(m)\cos \left(\frac{2\pi }{M}ml_{0} \right) \\
g_{4}(m)=g(m)\sin \left(\frac{2\pi }{M}ml_{0} \right)
\end{eqnarray*}

podemos escribir la parte real como

\[
h_{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.

\[
h_{I}(n,m)=g_{1}(n)*g_{4}(m)+g_{2}(n)*g_{3}(m)
\]

En la figura \ref{fig:gabor}, podemos ver unos ejemplos del filtro de Gabor para diferentes valores de frecuencia espacial.

\begin{figure}[h!]
\begin{center}
{\normalsize 
\subfigure[Parte real, entonado en la frecuencia espacial (0,0)\label{fig:gabor_r00}]{\psfig{figure = Filtros/figuras/gabor_r00.eps,width=5.0cm}} \hfill %
\subfigure[Parte imaginaria, entonado en la frecuencia espacial (0,0)\label{fig:gabor_i00}]{\psfig{figure = Filtros/figuras/gabor_i00.eps,width=5.0cm}} 
\subfigure[Parte real, entonado en la frecuencia espacial (2,2)\label{fig:gabor_r22}]{\psfig{figure = Filtros/figuras/gabor_r22.eps,width=5.0cm}} \hfill 
\subfigure[Parte imaginaria, entonado en la frecuencia espacial (2,2)\label{fig:gabor_i22}]{\psfig{figure = Filtros/figuras/gabor_i22.eps,width=5.0cm}} \label{fig:gabor} }
\end{center}
\caption{Filtros de Gabor}
\end{figure}

\begin{algorithm}[H]
\begin{algorithmic}[1]
	\STATE{Dadas $k_0$,  $l_0$, $\sigma_n$, $\sigma_m$ y $I$}\\
	\FOR {$n = -3\sigma_n$  \TO  $3\sigma_n$}
		\STATE{$g[n] = \frac{1}{\sqrt{2 \pi} \sigma} exp^{-n^2/(2 \sigma_n^2)}$}
		\STATE {$g_1[n] = g[n] \times cos(2 \pi n k_0/N)$}\\
		\STATE {$g_2[n] = g[n] \times sin(2 \pi n k_0/N)$}\\
	\ENDFOR
	\FOR {$m = -3\sigma_m$  \TO  $3\sigma_m$}
		\STATE{$g[m] = \frac{1}{\sqrt{2 \pi} \sigma} exp^{-m^2/(2 \sigma_m^2)}$}
		\STATE {$g_3[m] = g[m] \times cos(2 \pi m l_0/M)$}\\
		\STATE {$g_4[m] = g[m] \times sin(2 \pi m l_0/M)$}\\
	\ENDFOR
	\STATE{Calcular $f_r^{(1)}[n,m] = I[n,m]*g_1[n]*g_3[m]$}
	\STATE{Calcular $f_r^{(2)}[n,m] = I[n,m]*g_2[n]*g_4[m]$}
	\STATE{Calcular $f_i^{(1)}[n,m] = I[n,m]*g_1[n]*g_4[m]$}
	\STATE{Calcular $f_i^{(2)}[n,m] = I[n,m]*g_2[n]*g_3[m]$}	
	\FOR {$n = 0$  \TO  $N-1$}
		\FOR {$m=0$ \TO $M-1$}
			\STATE {Calcular $f_r[n,m] = f_r^{(1)}[n,m] - f_r^{(2)}[n,m]$ }\\
			\STATE {Calcular $f_i[n,m] = f_i^{(1)}[n,m] + f_i^{(2)}[n,m]$ }\\
		\ENDFOR
	\ENDFOR
	\RETURN {$f_r$ y $f_i$}
\end{algorithmic}
\caption{Filtro de Gabor}\label{alg:Gabor}
\end{algorithm}

\section{Filtro de Wiener}

Dadas dos imágenes $f^{\prime }$ y $f$  y un kernel $h[n]$ podemos decir:

\[
f^{\prime }[n]=f*h[n]
\]

lo cual indica que si realizamos convolución con un kernel $h[n]$  de $f[n]$, tendremos una imagen $f\prime[n]$. Aplicando transformada de Fourier podemos intentar calcular el valor del kernel $h(n)$ haciendo:

\begin{eqnarray*}
f^{\prime }[n]=f*h[n] \\
F^{\prime }[k]=F[k]H[k]
\end{eqnarray*}

de esta expresión podemos despejar el valor de $H[k] = F^\prime[k]/F[k]$, 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

\begin{eqnarray*}
\left| f^{\prime }-f*g\right| ^{2} \\
\text{sugeto a }\left| l*g\right| ^{2}
\end{eqnarray*}

donde $f^{\prime }[n]$ es la imagen filtrada, $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 el kernel que intentamos encontrar. Para nuestro calculo procedemos

\begin{eqnarray*}
U(g)=\sum_{n=0}^{N-1}\left[ f^\prime [n] - f*g[n] \right]^{2}+\beta
\sum_{n=0}^{N-1}\left| l*g[n]\right| ^{2} \\
=\sum_{n=0}^{N-1}\left[ f^{\prime }[n]-\sum_{i=0}^{N-1}f[n-i]g[i]\right]
^{2}+\beta \sum_{n=0}^{N-1}\left[ \sum_{i=0}^{N-1}l[n-i]g[i]\right] ^{2}
\end{eqnarray*}

Para encontrar el valor de g que minimiza esta expresión hacemos

\begin{eqnarray*}
\frac{\partial U(g)}{\partial g[i]}=-2\sum_{n=0}^{N-1}\left[ f^{\prime
}[n]-\sum_{i=0}^{N-1}f[n-i]g[i]\right] f[n-i] \\
+2\beta \sum_{n=0}^{N-1}\left[ \sum_{i=0}^{N-1}l[n-i]g[i]\right] l[n-i] \\
=0
\end{eqnarray*}

reorganizado términos tenemos

\begin{eqnarray*}
\sum_{n=0}^{N-1}\left[ \sum_{i=0}^{N-1}f[n-i]g[i]\right] f[n-i]+\beta
\sum_{n=0}^{N-1}\left[ \sum_{i=0}^{N-1}l[n-i]g[i]\right]
l[n-i]=\sum_{n=0}^{N-1}f^{\prime }[n]f[n-i] \\
\left[ f*f^{r}+\beta l*l^{r}\right] *g[n]=f^{\prime }*f^{r}[n]
\end{eqnarray*}

Al aplicar transformada de Fourier tenemos:

\begin{eqnarray*}
G[k]=\frac{\left[ F^{\prime }[k]\right] ^{*}F[k]}{\left[ F[k]\right]
^{*}F[k]+\beta \left[ L[k]\right] ^{*}L[k]} \\
G[k]=\frac{\left[ F^{\prime }[k]\right] ^{*}F[k]}{\left| F[k])\right|
^{2}+\beta \left| L[k]\right| ^{2}}
\end{eqnarray*}

al sacar la antitrasformada de Fourier calculamos el kernel $g[n]$  que aproxima al kernel original $h[n]$.

En dos dimensiones el filtro queda

\begin{equation}
G[k,l]=\frac{\left[ F^{\prime }[k,l]\right] ^{*}F[k,l]}{\left| F[k,l]\right|
^{2}+\beta \left| L[k,l]\right| ^{2}}
\end{equation}

donde $ F^{\prime }[k,l] ^{*} $ es el conjugado de la imagen.

En general podemos considerar que $L[k,l] = cte$ y simplificar la expresión 

\begin{equation}
G[k,l]=\frac{\left[ F^{\prime }[k,l]\right] ^{*}F[k,l]}{\left| F[k,l]\right|
^{2}+\beta }
\end{equation}

\begin{algorithm}[H]
\begin{algorithmic}[1]
	\STATE{Dadas $f$,  $f'$ y $\beta$}\\
	\STATE{Calcular $F = TDF[f]$ y $F' = TDF[f']$} \\
	\FOR {$k = 0$  \TO  $N-1$}
		\FOR {$l=0$ \TO $M-1$}
			\STATE {Calcular $G[k,l] = \frac{F'[k,l]^*F[k,l]}{|F[k,l]|^2 + \beta}$ }\\
		\ENDFOR
	\ENDFOR
	\STATE{Calcular $g = TDF^{-1}[G]$}\\
	\RETURN g
\end{algorithmic}
\caption{Filtro de Wiener}\label{alg:Wiener}
\end{algorithm}

\section{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.

Consideremos una función muestreada con tamaño $N$, en el dominio de Fourier en la posición $k=N$ tenemos que $N$ representa la frecuencia angular $2\pi $, por lo tanto, podemos representar todas las frecuencia como múltiplos de $\frac{2\pi }{N}$\ como $\omega =\frac{2\pi }{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=\left( 1/f_{m}\right) seg$, sustituyendo esto tenemos:

\begin{eqnarray*}
f_{t}=\frac{k}{N}\text{ }[ciclos/\left( 1/f_{m}\right) seg] \\
f_{t}=\frac{f_{m}k}{N}\text{ }[ciclos/seg] \\
f_{t}=f_{m}f_{k}\text{ }[ciclos/seg]
\end{eqnarray*}
finalmente $k=\frac{f_{t}}{f_{m}}N.$

\subsection{Ejemplo}

Consideremos una señal $x(t)=sen(200\pi 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ño de la muestra $N=128$.

La frecuencia de la señal es $f_{t}=200\pi /(2\pi )=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\pi t)=sen(200\pi \frac{n}{800})=sen(\frac{2\pi }{128}16n)
\]

\section{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

\[
\widehat{X}_{k_{0}}(k)=X(k)\frac{1}{\sqrt{2\pi }\sigma _{k}}e^{\frac{-\left(
k-k_{0}\right) ^{2}}{2\sigma _{k}^{2}}}
\]

Para eliminar la señal hacemos

\[
\widehat{X}_{-k_{0}}(k)=X(k)\frac{1}{\sqrt{2\pi }\sigma _{k}}\left[ 1-e^{%
\frac{-\left( k-k_{0}\right) ^{2}}{2\sigma _{k}^{2}}}\right]
\]

y aplicamos la transformada de Fourier

\begin{eqnarray*}
\widehat{x}_{-k_{0}}(n) &=&x(n)*cte\left[ \delta \left( n\right) -g\left(
n\right) e^{\frac{2\pi k_{0}}{N}n}\right] \\
&=&cte\left[ x(n)-x(n)*g\left( n\right) e^{\frac{2\pi k_{0}}{N}n}\right] \\
&=&x(n)-\widehat{x}_{k_{0}}(n)
\end{eqnarray*}

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.

\section{Filtro de Mediana}

Un filtro de mediana consiste en buscar en una vecindad de tamaño $N\times 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\times M$ pixeles de menor a mayor y la mediana será el valor a la mitad del arreglo \cite{Parker97}. 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ían se remplazar la media o moda por los valores min o max de la ventana de convolución.

\section{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\left( n,m\right) =\bigvee_{k=-\infty }^{\infty }\bigvee_{l=-\infty
}^{\infty }x\left( k,l\right) h\left( n-k,m-l\right) =x\vee h
\]

donde $\bigvee $ representa la operación or y esta operación es conocida como la or convolución o dilatación. La implementación en OpenCV se muestra a continuación.

\begin{verbatim}
void Dilacion(CvMat *ent, CvMat *sal, CvMat *h) {
    int i, j, k, l;
    int nr, nc, w, s;
    
    nr = ent->rows;
    nc = sal->cols;
    w = h->rows;
    
    for(i=0; i<nr-w; i++)
        for(j=0; j<nc-w; j++){
            s = 0;
            for(k=0; k<w; k++)
                for(l=0; l<w; l++) {
                    s += punto(ent, i+k, j+l)*punto(h, k, l);
                }
            punto(sal, i, j) = s >=1 ? 1 : 0;
        }
}
\end{verbatim}

Por otro lado podemos calcular la and convolución como

\[
y\left( n,m\right) =\bigwedge_{k=-\infty }^{\infty }\bigwedge_{l=-\infty
}^{\infty }x\left( k,l\right) h\left( n-k,m-l\right) =x\wedge h
\]

La or convolución también es conocida como operación de dilación y la and convolución como contracción o erosión. 

\begin{verbatim}
void Erosion(CvMat *ent, CvMat *sal, CvMat *h) {
    int i, j, k, l;
    int nr, nc, w, s;
    
    nr = ent->rows;
    nc = sal->cols;
    w = h->rows;
    
    for(i=0; i<nr-w; i++)
        for(j=0; j<nc-w; j++){
            s = 1;
            for(k=0; k<w; k++)
                for(l=0; l<w; l++) {
                        s *= punto(ent, i+k, j+l)*punto(h, k, l);
                    }
                punto(sal, i, j) = s;
            }
}
\end{verbatim}

Si hacemos la aplicación de ambos $\left(x\vee h\right) \wedge 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 $\left( x\vee h\right) \wedge h$ y la apertura a la operación $\left( x\wedge h\right) \vee h$

Para obtener los bordes de la imagen podemos hacer la diferencia de la imagen original con la imagen dilatada haciendo $x-\left( x\vee h\right) .$

\section{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.

\subsection{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}
\]

\subsection{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=\left\{ 1,2,3,4,5,6,7\right\} $, donde cada uno de los eventos son igualmente probables y estos eventos se agrupan en dos conjuntos $A=\left\{ 1,2,3,4\right\} $ y $B=\left\{ 3,4,5,6,7\right\}$, tal que $S=A\cup 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 \ref{fig:conjunto}).

\begin{figure}[h!]
{\normalsize \centering
\centerline{\psfig{figure = Filtros/figuras/conjunto.eps,width=8.0cm}} }
\caption{Ejemplo de dos conjuntos}
\label{fig:conjunto}
\end{figure}

Note que existe una intersección entre los conjuntos $A$ y $B$, en notación de conjuntos este conjunto es $A\cap 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\left\{ A\cap B\right\} }{n\left\{ B\right\} }
\]

si dividimos esta ecuación entre el total de elementos del espacio muestral tenemos

\begin{eqnarray*}
P(A|B)=\frac{n\left\{ A\cap B\right\} /N}{n\left\{ B\right\} /N} \\
=\frac{P(A\cap B)}{P(B)}
\end{eqnarray*}

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\cap B)}{P(B)}
\]

\subsection{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\cap B)=P(A)P(B).$ 

\subsection{Regla de Bayes}

Si los eventos $B_{1},B_{2},...,B_{k}$ constituyen una partición del espacio muestral S, donde $P\left( B_{i}\right) \neq 0$, entonces, para cualquier evento A de S

\[
P\left( B_{r}|A\right) =\frac{P\left( B_{r}\right) P\left( A|B_{r}\right) }{%
\sum_{i=1}^{k}P\left( B_{i}\right) P\left( A|B_{i}\right) }
\]
ver \cite{Walpole87}

\subsection{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 $\sigma ^{2}=1/2$, la cual se expresa como:

\[
p(g_{i}|f_{i})=\frac{1}{C_{1}}e^{-\left( f_{i}-g_{i}\right) ^{2}}
\]

considerando que los datos son generados de manera independiente tenemos que la distribución de probabilidades conjunta es

\[
p(g|f)=\prod_{i=0}^{N-1}\frac{1}{C_{1}}e^{-\left( f_{i}-g_{i}\right) ^{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)=\prod_{i=0}^{N-1}\frac{1}{C_{2}}e^{-\lambda \nabla
^{2}f_{i}}=\prod_{i=0}^{N-1}\frac{1}{C_{2}}e^{-\lambda \left(
f_{i}-f_{i-1}\right) ^{2}}
\]

aplicando la regla de Bayes, encontramos la expresión para la probabilidad a posteriori dada por:

\begin{eqnarray}
p(f|g)=\frac{p(g|f)*p(f)}{p(g)}  \label{poste} \\
p(f|g)=\prod_{i=0}^{N-1}\frac{1}{C}e^{-\left[ \left( f_{i}-g_{i}\right)
^{2}+\lambda \left( f_{i}-f_{i-1}\right) ^{2}\right] }  \nonumber \\
p(f|g)=\prod_{i=0}^{N-1}\frac{1}{C}e^{-U(f)}  \nonumber
\end{eqnarray}

La ecuación \ref{poste} 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)=\sum_{i=0}^{N-1}\left[ \left( f[i]-g[i]\right) ^{2}+\lambda \left(
f[i]-f[i-1]\right) ^{2}\right]
\]

la cual representa la energía potencial almacenada en un sistema de resortes acoplados, tal como se muestra en la figura \ref{fig:resortes}.

\begin{figure}[h!]
{\normalsize \centering
\centerline{\psfig{figure = Filtros/figuras/resortes.eps,width=8.0cm}} }
\caption{Sistema de resortes equivalente al filtro de membrana}
\label{fig:resortes}
\end{figure}

Para calcular el valor de $f$ que minimiza la función $U(f)$, hacemos $\frac{ \partial U(f)}{\partial f[i]}=0$

\begin{equation}
\frac{\partial U(f)}{\partial f[i]}=2\left( f[i]-g[i]\right) +2\lambda
\left( f[i]-f[i-1]\right) -2\lambda \left( f[i]-f[i]\right) =0
\label{mem}
\end{equation}

lo cual nos da el siguiente sistema lineal de ecuaciones :

\[
\left[
\begin{array}{llllllll}
1+\lambda & -\lambda &  &  &  &  &  &  \\
-\lambda & 1+2\lambda & -\lambda &  &  &  &  &  \\
&  & \ddots & \ddots & \ddots &  &  &  \\
&  & -\lambda & 1+2\lambda & -\lambda &  &  &  \\
&  &  & -\lambda & 1+2\lambda & -\lambda &  &  \\
&  &  &  & \ddots & \ddots & \ddots &  \\
&  &  &  &  & -\lambda & 1+2\lambda & -\lambda \\
&  &  &  &  &  & -\lambda & 1+\lambda
\end{array}
\right] \left[
\begin{array}{l}
f[0] \\
f[1] \\
\vdots \\
f[i] \\
f[i+1] \\
\vdots \\
f[N-2] \\
f[N-1]
\end{array}
\right] =\left[
\begin{array}{l}
g[0] \\
g[1] \\
\vdots \\
g[i] \\
g[i+1] \\
\vdots \\
g[N-2] \\
g[N-1]
\end{array}
\right]
\]

\subsection{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 \ref{mem} y agrupando los términos obtenemos:

\[
F[k]+\lambda \left( F[k]-F[k]e^{-\left( \frac{2\pi }{N}\right) k}\right)
-\lambda \left( F[k]e^{\left( \frac{2\pi }{N}\right) k}-F[k]\right) =G[k]
\]

agrupando términos obtenemos

\begin{eqnarray*}
\left[ 1+2\lambda -\lambda e^{-\left( \frac{2\pi }{N}\right) k}-\lambda
e^{\left( \frac{2\pi }{N}\right) k}\right] F[k]=G[k] \\
\left[ 1+2\lambda -2\lambda \cos \left( \frac{2\pi }{N}k\right) \right]
F[k]=G[k]
\end{eqnarray*}

finalmente tenemos que

\begin{eqnarray*}
F[k]=\frac{1}{\left[ 1+2\lambda -2\lambda \cos \left( \frac{2\pi }{N} k\right) \right] }G[k]=H[k]G[k] \\
H[k]=\frac{1}{\left[ 1+2\lambda -2\lambda \cos \left( \frac{2\pi }{N}k\right) \right] }
\end{eqnarray*}

donde la función $H[k]$ es un filtro pasa bajas equivalente al filtro de membrana. Note que cuando el valor de $\lambda $ aumenta tenemos que el ancho de banda se disminuye, en las figuras \ref{fig:lambda1} y \ref{fig:lambda1000} podemos ver la respuesta a la frecuencia de este filtro para valores de $\lambda $ 1 y 1000 respectivamente.

\begin{figure}[h!]
{\normalsize \centering
\centerline{\psfig{figure = Filtros/figuras/lambda1.eps,width=8.0cm}} }
\caption{Respuesta a la frecuencia para el filtro de membrana con $\lambda
=1 $}
\label{fig:lambda1}
\end{figure}

\begin{figure}[h!]
{\normalsize \centering
\centerline{\psfig{figure = Filtros/figuras/lambda1000.eps,width=8.0cm}} }
\caption{Respuesta a la frecuencia para el filtro de membrana con $\lambda=1000$}
\label{fig:lambda1000}
\end{figure}

En dos dimensiones tenemos

\[
U(f)=\sum_{n=0}^{N-1}\sum_{m=0}^{M-1}\left[ \left( f[n,m]-g[n,m]\right) ^{2}+\lambda \left(
f[n,m]-f[n-1,m]\right) ^{2} + \lambda \left(f[n,m]-f[n,m-1]\right) ^{2}\right]
\]

La derivada es

\begin{eqnarray*}
\frac{\partial U(f)}{\partial f[n,m]}=2\left( f[n,m]-g[n,m]\right) + \\
2\lambda \left( f[n,m]-f[n-1, m]\right) -2\lambda \left( f[n+1,m]-f[n,m]\right) \\
2\lambda \left( f[n,m]-f[n, m-1]\right) -2\lambda \left( f[n,m]-f[n,m+1]\right) =0
\label{mem}
\end{eqnarray*}

La Transformada Discreta de Fourier queda

\begin{eqnarray*}
\left[ 1+4\lambda -\lambda e^{\left( \frac{2\pi }{N}\right) k} - \lambda e^{-\left( \frac{2\pi }{M}\right) k}-\lambda
e^{\left( \frac{2\pi }{M}\right) l} -\lambda e^{\left( -\frac{2\pi }{M}\right) l}\right] F[k,l]=G[k,l] \\
\left[ 1+4\lambda -2\lambda \cos \left( \frac{2\pi }{N}k\right)   -2\lambda \cos \left( \frac{2\pi }{M}l\right) \right]
F[k,l]=G[k,l]
\end{eqnarray*}

El Filtro en el dominio de Fourier queda

\begin{equation*}
H[k,l] = \frac{1}{\left[ 1+4\lambda -2\lambda \cos \left( \frac{2\pi }{N}k\right)   -2\lambda \cos \left( \frac{2\pi }{M}l\right) \right]} 
\end{equation*}

\begin{algorithm}[H]
\begin{algorithmic}[1]
	\STATE{Dadas $g$ y $\lambda$}\\
	\STATE{Calcular  $G = TDF[g]$ } \\
	\FOR {$k = 0$  \TO  $N-1$}
		\FOR {$l=0$ \TO $M-1$}
			\STATE {Calcular $F[k,l] = \frac{G[k,l]}{1+4\lambda - 2\lambda cos( 2\pi k/N) -2\lambda cos(2 \pi l/M)}$ }\\
		\ENDFOR
	\ENDFOR
	\STATE{Calcular $f = TDF^{-1}[F]$}\\
	\RETURN f
\end{algorithmic}
\caption{Filtro de Membrana}\label{alg:membrana}
\end{algorithm}