% 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[Transformadas]{Transformadas Coseno y Wavelet}

\section{Transformada Coseno a partir de la Transformada de Fourier}

Partiremos de una señal $x$ de tamaño $2N-1$ dada por \eqref{eq:sen}  la cual es real y par.  La transformada de Fourier de esta señal también será real y par de acuerdo con la propiedades de la Transformada de Fourier.  

\begin{equation}
x(n) = x(2N -1 -n) \label{eq:sen}
\end{equation}

\noindent 
de otra forma

\begin{eqnarray*}
x(0) = X(2N-1) \\
x(1) = X(2N-2) \\
x(2) = X(2N-3)\\
\vdots
\end{eqnarray*}

Un ejemplo de esta señal podría ser $x =[1,2,3,4,4,3,2]$. La transformada de Fourier de cualquier señal $x$ la con la forma dada por \eqref{eq:sen} la podemos calcular como:

\begin{equation*}
X(k) = \sum_{n=0}^{2N-2} x(n) e^{-j \frac{2 \pi n k}{2N-1}}
\end{equation*}

La sumatoria la podemos dividir en dos partes 

\begin{equation*}
X(k) = \sum_{n=0}^{N-1} x(n) e^{-j \frac{2 \pi n k}{2N-1}} + \sum_{n=N}^{2N-2} x(n) e^{-j \frac{2 \pi n k}{2N-1}}
\end{equation*}

Si en la segunda sumatoria hacemos el cambio de variable $n = 2N -1 -\hat n$ tenemos que $\hat n = 2 N -1 - n$ y podemos escribir

\begin{equation*}
X(k) = \sum_{n=0}^{N-1} x(n) e^{-j \frac{2 \pi n k}{2N-1}} + \sum_{\hat n=N-1}^{1} x(2N -1 -\hat n) e^{-j \frac{2 \pi (2N -1 -\hat n) k}{2N-1}}
\end{equation*}

Reorganizando términos tenemos

\begin{equation*}
X(k) = x(0) + \sum_{n=1}^{N-1} x(n) e^{-j \frac{2 \pi n k}{2N-1}} + \sum_{n=1}^{N-1} x(2N -1 - n) e^{-j \frac{2 \pi (- n) k}{2N-1}} \cancelto{1}{e^{-j \frac{2 \pi (2N -1) k}{2N-1}}}
\end{equation*}

Dado que $x(n) = x(2N -1 - n)$ podemos simplificar

\begin{equation*}
X(k) = x(0) + \sum_{n=1}^{N-1} x(n) \left ( e^{-j \frac{2 \pi n k}{2N-1}} +  e^{+j \frac{2 \pi n k}{2N-1}} \right)
\end{equation*}

\begin{equation}
X(k) = x(0) + 2\sum_{n=1}^{N-1} x(n) Cos \left(\frac{2 \pi n k}{2N-1}  \right) \label{eq:Tcos}
\end{equation}

La función que resulta ya es una sucesión de cosenos por lo cual podemos decir que es una Transformada Coseno implementada a partir de la transformada de Fourier.

El Código en Java para hacer la implementación de la Transformada Coseno es:

\begin{verbatim}
    static public void MiTCos(float x[], float X[]) {
        int N = x.length;
        int n, k;
        double w;
        
        w = (2.0*Math.PI/(2.0*N-1.0));
        
        for(k=0; k<N; k++){
            X[k] = x[0];
            for(n=1; n<N; n++)
                X[k] += x[n]*2.0*Math.cos(w*k*n);
        }
    }

\end{verbatim}

\subsection{Transformada Inversa}

Para llevar a cabo la transformada inversa consideraremos los mismos principios que la transformada Coseno. Por definición la transformada Inversa es:

\begin{equation*}
x(n) = \frac{1}{2N -1}\sum_{k=0}^{2N-2} X(k) e^{j \frac{2 \pi n k}{2N-1}}
\end{equation*}

La sumatoria la podemos dividir en dos partes 

\begin{equation*}
x(n) = \frac{1}{2N -1}\sum_{k=0}^{N-1} X(k) e^{j \frac{2 \pi n k}{2N-1}} + \frac{1}{2N -1}\sum_{k=N}^{2N-2} X(k) e^{j \frac{2 \pi n k}{2N-1}}
\end{equation*}

Si en la segunda sumatoria hacemos el cambio de variable $k = 2N -1 -\hat k$ tenemos que $\hat k = 2 N -1 - k$ y podemos escribir

\begin{equation*}
x(n) = \frac{1}{2N -1}\sum_{k=0}^{N-1} X(k) e^{j \frac{2 \pi n k}{2N-1}} + \frac{1}{2N -1}\sum_{\hat k=N-1}^{1} X(2N -1 -\hat k) e^{j \frac{2 \pi (2N -1 -\hat k) n}{2N-1}}
\end{equation*}

Reorganizando términos tenemos

\begin{equation*}
x(n) =\frac{1}{2N -1} \left ( X(k) + \sum_{k=1}^{N-1} X(k) e^{j \frac{2 \pi n k}{2N-1}} + \sum_{k=1}^{N-1} X(2N -1 - k) e^{j \frac{2 \pi n (-k)}{2N-1}} \cancelto{1}{e^{-j \frac{2 \pi (2N -1) n}{2N-1}}} \right)
\end{equation*}

Dado que $X(k) = X(2N -1 - k)$ podemos simplificar

\begin{equation*}
x(n) = \frac{1}{2N -1} \left( X(0) + \sum_{k=1}^{N-1} X(k) \left ( e^{j \frac{2 \pi n k}{2N-1}} +  e^{-j \frac{2 \pi n k}{2N-1}} \right) \right )
\end{equation*}

Finalmente

\begin{equation}
x(n) = \frac{1}{2N -1} \left(X(0) + 2\sum_{k=1}^{N-1} X(k) Cos \left(\frac{2 \pi n k}{2N-1}  \right) \label{eq:Tcosinv} \right)
\end{equation}

\subsection{Ejemplo}

Consideremos una señal $x_1 = [1,2,3,4,4,3,2]$  y otra señal $x_2 = [1,2,3,4]$ calcular:

\begin{enumerate}
\item Calcular la transformada de Fourier $x_1$,
\item Calcular la transformada Coseno dada por \eqref{eq:Tcos} para $x_2$ ,
\item Calcular la transformada Coseno Inversa dada por \eqref{eq:Tcosinv} para $X_2(k)$ y 
\item Concluir
\end{enumerate}

La Transformada de Fourier de $x_1$ es  

\begin{equation*}
X_1(k) = [19,  -5.04892 , -0.307979, -0.643104, -0.643104, -0.307979 , -5.04892]
\end{equation*}

La Transformada Coseno de $x_2$  es

\begin{equation*}
X_2(k) = [19.0, -5.04892, -0.307979, -0.643104]
\end{equation*}


La transformada Coseno inversa para $X_2(k)$ es 
\begin{equation*}
\hat x_2(n) = [1,2,3,4]
\end{equation*}

Note que los términos de $X_2(k)$ corresponden exactamente a los primero 4 términos de $X1(k)$.

\section{Definición de la DCT}

Para derivar la DCT de una señal real de con $N$ muestras $x'= [x'(0), x'(1), \cdots, x'(N-1)]$, primeramente construiremos una señal $x'(n)$ con $2N$ puntos 

\begin{equation}
x(n) = 
\left \{
\begin{array} {ll}
x'(n)     & (  0 \leq n \leq N-1) \\
x'(-n-1) & (-N\leq n \leq -1) 
\end{array}
\right .
\label{eq:senalCos}
\end{equation}


\begin{figure}[h!]
\begin{center}
   \subfigure[$I_1$ ]{\includegraphics[scale=0.7]{Transformadas/figuras/senal1.eps} \label{fig:Cosgraf0}}
   \subfigure[$I_1$ ]{\includegraphics[scale=0.7]{Transformadas/figuras/senal2.eps} \label{fig:Cosgraf1}}
\label{fig:Cosgraf} 
\end{center}
\caption{Reconstrucción de la señal usando los coeficientes de Fourier}
\end{figure}

En la Figura \ref{fig:Cosgraf0} se muestra la imagen correspondiente a la señal generada con \eqref{eq:senalCos}.
Pare esta señal $x(n)$ asumiremos que tiene periodo $2N$ y es par con respecto al punto $n=-1/2$.  Si aplicamos una translación de $1/2$  hacia la derecha de la señal $x(n)$ haciendo el cambio de variable $n' = n + 1/2$, entonces tenemos que $x(n) = x(n'-1/2)$ y la señal resultante es par con respecto a $n'=0$ como se muestra en la Figura \ref{fig:Cosgraf1}.

La transformada discreta de Fourier para la señal $x(n'-1/2)$ la podemos calcular con

\begin{equation*}
X(k) = \sum_{n'=-N+1/2}^{N-1/2} x\left(n' -\frac{1}{2}\right) e^{-j 2 \pi n' k/2N}
\end{equation*}

Podemos hacer la extensión de la exponencial en su forma rectangular y dado que $x$ es real expresar

\begin{equation*}
X(k) =\sum_{n'=-N+1/2}^{N-1/2} x\left(n' -\frac{1}{2}\right) cos( 2\pi n' k/2N) - j \sum_{n'=-N+1/2}^{N-1/2} x\left(n' -\frac{1}{2}\right) sin(2 \pi n' k/2N)
\end{equation*}

Dado que la señal $x$ es par y que el $cos(.)$ y el $sin(.)$ son par e impar respectivamente con respecto a $n'=0$ entonces la segunda sumatoria será cero.

\begin{equation*}
X(k) = \sum_{n'=-N+1/2}^{N-1/2} x\left(n' -\frac{1}{2}\right) cos( 2\pi n' k/2N) 
\end{equation*}

Separando la sumatoria en dos términos

\begin{equation*}
X(k) =   \sum_{n'=-N+1/2}^{-1/2} x\left(n' -\frac{1}{2}\right) cos( 2\pi n' k/2N)  + \sum_{n'=1/2}^{N-1/2} x\left(n' -\frac{1}{2}\right) cos( 2\pi n' k/2N) 
\end{equation*}

\begin{equation*}
X(k) =  \sum_{n'=1/2}^{N-1/2} x\left(-n' +\frac{1}{2}\right) cos( -2\pi n' k/2N)  + \sum_{n'=1/2}^{N-1/2} x\left(n' -\frac{1}{2}\right) cos(2 \pi n' k/2N) 
\end{equation*}

Dado que la señal es par entonces $x(n' -1/2) = x(-n' + 1/2)$ podemos simplificar la expresión y obtener

\begin{equation*}
X(k) =  2 \sum_{n'=1/2}^{N-1/2} x\left(n' -\frac{1}{2}\right) cos( 2\pi n' k/2N) 
\end{equation*}

Finalmente si hacemos el cambio de variable $n= n'-1/2$ obtenemos 

\begin{equation}
X(k) = 2 \sum_{n=0}^{N-1} x(n) cos \left( \frac{\pi (n+1/2) k}{N} \right) 
\end{equation}

\subsection{Transformada Coseno Inversa}

Para derivar la DCT inversa de una señal real $X'$  con $N$ muestras $X'= [X'(0), X'(1), \cdots, X'(N-1)]$, primeramente construiremos una  nueva señal $X$ con $2N$ puntos dada por \eqref{eq:senalCosInv}

\begin{equation}
X(k) = 
\left \{
\begin{array} {ll}
X'(k)     & (  0 \leq k \leq N-1) \\
X'(-k-1) & (-N\leq k \leq -1) 
\end{array}
\right .
\label{eq:senalCosInv}
\end{equation}

Pare esta señal $X(k)$ asumiremos que tiene periodo $2N$ y es par con respecto al punto $k=-1/2$.  Si aplicamos una translación de $1/2$  hacia la derecha de la señal $X(k)$ haciendo el cambio de variable $k' = k + 1/2$, entonces tenemos que $X(k) = X(k'-1/2)$ y la señal resultante es par con respecto a $k'=0$.

La  transformada discreta de Fourier inversa para la señal $X(k'-1/2)$ la podemos calcular con

\begin{equation*}
x(n) = \frac{1}{2N} \sum_{k'=-N+1/2}^{N-1/2} X\left(k' -\frac{1}{2}\right) e^{j 2 \pi n k'/2N}
\end{equation*}

Podemos hacer la extensión de la exponencial en su forma rectangular y dado que $X$ es real podemos expresar

\begin{eqnarray*}
x(n) &=&\frac{1}{2N}  \sum_{k'=-N+1/2}^{N-1/2} X\left(k' -\frac{1}{2}\right) cos( 2\pi n k'/2N) \\
&+& \frac{j}{2N}  \sum_{k'=-N+1/2}^{N-1/2} X\left(k' -\frac{1}{2}\right) sin(2 \pi n k'/2N) 
\end{eqnarray*}

Dado que la señal $X$ es par y que el $cos(.)$ y el $sin(.)$ son par e impar respectivamente con respecto a $n'=0$ entonces la segunda sumatoria será cero.

\begin{equation*}
x(n) = \frac{1}{2N} \sum_{k'=-N+1/2}^{N-1/2} X\left(k' -\frac{1}{2}\right) cos( 2\pi n k'/2N) 
\end{equation*}

Separando la sumatoria en dos términos

\begin{equation*}
x(n) =   \frac{1}{2N} \left( \sum_{k'=-N+1/2}^{-1/2} X\left(k' -\frac{1}{2}\right) cos( 2\pi n k'/2N)  + \sum_{k'=1/2}^{N-1/2} X\left(k' -\frac{1}{2}\right) cos( 2\pi n k'/2N) \right)
\end{equation*}

\begin{equation*}
x(n) =   \frac{1}{2N} \left( \sum_{k'=1/2}^{N-1/2} X\left(-k' +\frac{1}{2}\right) cos( -2\pi n k'/2N)  + \sum_{k'=1/2}^{N-1/2} X\left(k' -\frac{1}{2}\right) cos(2 \pi n k'/2N)  \right)
\end{equation*}

Dado que la señal es par entonces, es decir $X(k' -1/2) = X(-k' + 1/2)$, podemos simplificar la expresión y obtener

\begin{equation*}
x(n) =  \frac{2}{2N}  \sum_{k'=1/2}^{N-1/2} X\left(k' -\frac{1}{2}\right) cos( 2\pi n k'/2N) 
\end{equation*}

Finalmente si hacemos el cambio de variable $k= k'-1/2$ obtenemos 

\begin{equation}
x(n) = \frac{1}{N} \sum_{k=0}^{N-1} X(k) cos \left( \frac{\pi n(k+1/2) }{N} \right) 
\end{equation}


\subsection{Transformadas Coseno Utilizadas}

Se pude tener un sin número de definiciones de la transformada coseno, sin embargo, centraremos en cuatro principalmente:

\begin{equation*}
\begin{array}{|c|l|}
	\hline
	\text{tipo} & \text{Definición}\\
	\hline
	&\\
	1 & X_1(k) = \sqrt{\frac{2}{N-1}} \left( \frac{x(0)}{2} + (-1)^s \frac{x(n-1)}{2} \sum_{n=1}^{N-2} x(n) cos \left( \frac{\pi n k}{N-1} \right) \right )\\
	& \\
	\hline
	&\\
	2 & X_2(k) = \frac{1}{\sqrt{N}} \sum_{n=0}^{N-1} x(n) cos \left( \frac{\pi (n+1/2) k}{N} \right) \\
	&\\
	\hline
	&\\
	3 & X_3(k) = \frac{1}{\sqrt{N}} \left(  x(0) + 2\sum_{n=1}^{N-1} x(n) cos \left( \frac{\pi n (k+1/2)}{N} \right)  \right)\\
	&\\
	\hline
	&\\
	4 & X_4(k) = \sqrt{\frac{2}{N}} \sum_{n=0}^{N-1} x(n) cos \left( \frac{\pi (n+1/2) k}{N} \right) \\
	&\\
	\hline
\end{array}
\end{equation*}

Las transformadas discreta de coseno inversa para los tipos 1,2,3 y 4 es la tipo 1, 3, 2 y 4 respectivamente.

\subsection{Ejemplo de aplicación}

Dada una señal $x[n] = sin(0.5 n) + N(0, 0.3)$ queremos obtener una señal filtrada eliminando el $75\%$ de los coeficientes correspondientes a las altas frecuencias. Considere que se tienen $N=64$ muestras.

En la Figura \ref{fig:eje1d} se muestran los resultados en cuatro sub figuras. En la Figura \ref{fig:orig} se muestra la señal original la cual tiene una frecuencia $w = 0.5$. La Figura \ref{fig:origT} muestra la transformada de Coseno de la señal $x$, note como cerca de la frecuencia fundamental se muestra el máximo de la señal. La Figura \ref{fig:final} muestra la Transformada Coseno una vez que se eliminaron el $75\%$ de las altas frecuencias y finalmente en la Figura \ref{fig:finalT} sobre impuestas la señal original y la señal muestreada.

La frecuencia fundamental la podemos calcular haciendo

\begin{equation*}
0.5 n = \frac{\pi (k + 1/2) n}{64}
\end{equation*}

despejando tenemos

\begin{equation*}
k = \frac{64 \times 0.5}{\pi} - \frac{1}{2} = 9.6858
\end{equation*}

\begin{figure}[h!]
\begin{center}
   \subfigure[$x$ ]{\includegraphics[scale=0.35]{Transformadas/figuras/Eje_1d_X0.jpg} \label{fig:orig}}
   \subfigure[$X$ ]{\includegraphics[scale=0.35]{Transformadas/figuras/Eje_1d_X.jpg} \label{fig:origT}}
   \subfigure[$x_f$ ]{\includegraphics[scale=0.35]{Transformadas/figuras/Eje_1d_Xf0.jpg} \label{fig:final}}
   \subfigure[$X_f$ ]{\includegraphics[scale=0.35]{Transformadas/figuras/Eje_1d_Xf.jpg} \label{fig:finalT}}
\end{center}
\caption{Ejemplo de filtrado y compactación de una señal utilizando Transformada Coseno}
\label{fig:eje1d}
\end{figure} 

\subsection{Transformada Discreta Coseno en dos Dimensiones}

En forma similar al cálculo de la TDC en una dimensión donde se definió una señal par \eqref{eq:senalCos},  en dos dimensiones partimos de la transformada discreta de Fourier para la señal bidimensional par desplazada $x(n'-1/2, m'-1/2)$.

\begin{equation*}
X(k, l) = \sum_{n'=-N+1/2}^{N-1/2} \sum_{m'=-M+1/2}^{M-1/2}x\left(n' -\frac{1}{2}, m' -\frac{1}{2}\right) e^{-j 2 \pi n' k/2N - j 2 \pi m' k/2M}
\end{equation*}

Aplicando la propiedad asociativa de la suma y la propiedad de los exponentes podemos escribir como \eqref{eq:DCT_2D}

\begin{equation}
X(k, l) = \sum_{n'=-N+1/2}^{N-1/2} \left (\sum_{m'=-M+1/2}^{M-1/2}x\left(n' -\frac{1}{2}, m' -\frac{1}{2}\right) e^{ - j 2 \pi m' k/2M} \right) e^{-j 2 \pi n' k/2N }
\label{eq:DCT_2D}
\end{equation}

Si definimos  $X_c$ como \eqref{eq:DCT_2D_paso} por podemos escribir \eqref{eq:DCT_2D}. Esto significa que podemos calcular la transformada coseno primero por columnas $X_c$  \eqref{eq:DCT_2D} y posteriormente por renglones \eqref{eq:DCT_2D_final}.

\begin{equation}
X_c(k, l) = \sum_{m'=-M+1/2}^{M-1/2}x\left(n' -\frac{1}{2}, m' -\frac{1}{2}\right) e^{ - j 2 \pi m' k/2M} 
\label{eq:DCT_2D_paso}
\end{equation}

\begin{equation}
X(k, l) = \sum_{n'=-N+1/2}^{N-1/2} X_c\left(n' -\frac{1}{2}, m' -\frac{1}{2}\right) e^{-j 2 \pi n' k/2N }
\label{eq:DCT_2D_final}
\end{equation}

\subsection{Compresión de Imágenes utilizando TDC}

Una aplicación de la TDC es la compresión de  imágenes. Para mostrar comó se lleva  cabo la compresión vamos a considerar una imagen a la cual le aplicaremos la TDC y nos quedaremos con el 100, 80, 25 y 10\% de las bajas frecuencias, a estas imágenes recortadas se les aplica la TDC inversa y el resultado se muestra en la Figura \ref{fig:eje2d}.  En el primer renglón Figuras \ref{fig:lena100}, \ref{fig:lena80}, \ref{fig:lena25} y \ref{fig:lena10}, se muestra las imágenes,  resultado de dejar pasar 100, 80, 25 y 10\% de las bajas frecuencias. En el segundo renglón Figuras \ref{fig:mask100}, \ref{fig:mask80}, \ref{fig:mask25} y \ref{fig:mask10} se muestra en blanco los valores de la TDC que se conservan y en negro los que se eliminan. Note que se puede tener un recorte significativo de altas frecuencias y la señal que se obtiene mantiene mucha información de la imagen original, condición que permite hacer la compresión de la imagen. 

\begin{figure}[h!]
\begin{center}
   \subfigure[Lena 100 ]{\includegraphics[scale=0.35]{Transformadas/figuras/salida_100.jpg} \label{fig:lena100}}
   \subfigure[Lena 80 ]{\includegraphics[scale=0.35]{Transformadas/figuras/salida_80.jpg} \label{fig:lena80}}
   \subfigure[Lena 25 ]{\includegraphics[scale=0.35]{Transformadas/figuras/salida_25.jpg} \label{fig:lena25}}
   \subfigure[Lena 10 ]{\includegraphics[scale=0.35]{Transformadas/figuras/salida_10.jpg} \label{fig:lena10}}
   \subfigure[Mascara 100 ]{\includegraphics[scale=0.35]{Transformadas/figuras/mascara_100.jpg} \label{fig:mask100}}
   \subfigure[Mascara 80]{\includegraphics[scale=0.35]{Transformadas/figuras/mascara_80.jpg} \label{fig:mask80}}
   \subfigure[Mascara 25 ]{\includegraphics[scale=0.35]{Transformadas/figuras/mascara_25.jpg} \label{fig:mask25}}
   \subfigure[Mascara 10 ]{\includegraphics[scale=0.35]{Transformadas/figuras/mascara_10.jpg} \label{fig:mask10}}
\end{center}
\caption{Ejemplo de filtrado y compactación de una imagen utilizando Transformada Coseno}
\label{fig:eje2d}
\end{figure} 

\section{Transformada Wavelet}

La Transformada Discreta Wavelet TDW de una señal $x$ se calcula pasándolo a través de una serie de filtros. Primero, la señal $x$ se pasan a través de un filtro de paso bajas con respuesta al impulso $g$  para posteriormente muestrearla a la mitad. Esta operación se puede representar como la convolución dada por \eqref{w_baja}.

\begin{equation}
	y_{low} [n] = \sum_{k = \infty}^\infty g[k] x[2n-k] 
	\label{w_baja}
\end{equation}

Posteriormente la señal $x$ es convolucionada con un filtro pasa altas y mestreada a la mitad de acuerdo con \eqref{w_alta}. El proceso completo se muestra en la Figura \ref{fig:ej_wavelet}

\begin{equation}
	y_{high} [n] = \sum_{k= \infty}^\infty h[k] x[2n-k]
	\label{w_alta}
\end{equation}

\begin{figure}
  \centering
    \includegraphics{Transformadas/figuras/wavelet.jpg}
  \caption{Ejemplo de aplicación de la transformada Wavelet}
  \label{fig:ej_wavelet}
\end{figure}


El código Java para realizar la TDW se muestra a continuación:

\begin{verbatim}
       static void pwt(double a[], double cc[], double cr[], int n, int isign) {
        double ai,ai1,wksp[];
        int i, ii, j, k, nh, ncof;
        
        if (n < 4) return;
        
        ncof = cc.length; // Numero de coeficientes de los kerneles
        //System.out.println("Numero de coeficientes " + ncof);
 
        nh=n >> 1;        // Calcula la mitad del arreglo
        n = nh << 1;      // forza a que n sea par
        wksp = new double[n];
        
        for (j=0;j<n;j++) wksp[j]=0.0;
        
        if (isign >= 0) {
            for (ii=0,i=0;i<n;i+=2,ii++) {
                for (k=0;k<ncof;k++) {
                    j = (i+k)< n ? i+k : i+k-n;
                    wksp[ii] += cc[k]*a[j];
                    wksp[ii+nh] += cr[k]*a[j];
                }
            }
        }
        else {
            for (ii=0,i=0;i<n;i+=2,ii++) {
                ai=a[ii];
                ai1=a[ii+nh];
                for (k=0;k<ncof;k++) {
                    j= (i+k) < n ? i+k : i+k-n;
                    wksp[j] += cc[k]*ai;
                    wksp[j] += cr[k]*ai1;
                }
            }
        }
        
        for (j=0;j<n;j++) a[j]=wksp[j];
    }
\end{verbatim}

\subsection{Ejemplo}

Dada la señal $x =  [ 56, 40, 8, 24, 48, 48, 40, 10 ]$ y los filtros  $g=[1/2, 1/2]$ y $h=[1/2, -1/2]$ determinar la transformada discreta wavelet de $x$.

Los resultados de la TDW se muestran en la tabla \ref{t:DWT}. En el primer renglón se muestra la señal original y en los siguiente renglones la TDW para diferentes niveles $k = 1, 2, 3, \cdots$.  Note que el primer elemento del último renglón corresponde al promedio de la señal.

\begin{table}[]
\caption{Transformada Wavelet para una señal discreta}
\begin{center}
\begin{tabular}{|c|cccccccc|}
        \hline
        $x$ & 56 & 40 & 8 & 24 & 48 & 48 & 40 & 16 \\
        \hline
        \hline
        $X^{(k)}$ & \multicolumn{4}{c||}{$d_{low}$} & \multicolumn{4}{c|}{\color{red}{$d_{high}$}}\\
        \hline
	$X^{(1)}$ & 48 & 16 & 48 & 28 & \textcolor{red}{8} & \textcolor{red}{-8} & \textcolor{red}{0} & \textcolor{red}{12}\\
	\hline
	$X^{(2)}$ & 32 & 38 & \textcolor{red}{16} & \textcolor{red}{10} & \textcolor{red}{8} & \textcolor{red}{-8} & \textcolor{red}{0} & \textcolor{red}{12}\\
	\hline
	$X^{(3)}$ & 35 & \textcolor{red}{-3} & \textcolor{red}{16} & \textcolor{red}{10} & \textcolor{red}{8} & \textcolor{red}{-8} & \textcolor{red}{0} & \textcolor{red}{12}\\
	\hline
\end{tabular}
\end{center}
\label{t:DWT}
\end{table}%

\subsection{Antitransformada Wavelet}

Consideremos el caso en que se utilizan dos kerneles $g = [c, c]$ y $h = [c, -c]$. Con esta condición podemos ver que:

\begin{eqnarray*}
   d_{low}[i] = x[i] \times c + x[i+1] \times c \\
   d_{high}[i] = x[i] \times c - x[i+1] \times c
\end{eqnarray*}

Si sumamos y restamos $d_{low}[i]$ y $d_{high}[i]$ tenemos

\begin{eqnarray*}
   d_{low}[i] + d_{high}[i] = 2 x[i] \times c  \\
   d_{low}[i] - d_{high}[i] = 2 x[i+1] \times c 
\end{eqnarray*}

Finalmente obtenemos

\begin{eqnarray*}
  x[i]  = \frac{d_{low}[i] + d_{high}[i]}{ 2 c}  \\
  x[i+1] =  \frac{d_{low}[i] - d_{high}[i]}{ 2  c} 
\end{eqnarray*}

Esto nos permite ver qué la transformada wavelet es invertible y se puede mostrar que la invertibilidad no es exclusiva de kerneles con dos valores solamente. En la tabla \ref{t:IDWT} se muestra el procedimiento inverso para calcular la inversa de la TDW.

\begin{table}[]
\caption{Transformada Discreta Wavelet Inversa para una señal }
\begin{center}
\begin{tabular}{|c|cccccccc|}
        \hline
        $X$ & 35.0 & \textcolor{red}{-3.0} & \textcolor{red}{16.0} & \textcolor{red}{10.0} & \textcolor{red}{8.0} & \textcolor{red}{-8.0} & \textcolor{red}{0.0} & \textcolor{red}{12.0} \\
        \hline
        \hline
        $X^{(k)}$ & \multicolumn{4}{c||}{$d_{low}$} & \multicolumn{4}{c|}{\color{red}{$d_{high}$}}\\
        \hline
        \hline
 	$X^{(3)}$ & 35.0 & \textcolor{red}{-3.0} &  &  &  &  &  &  \\
 	$X^{(2)}$ & 32.0 & 38.0 & \textcolor{red}{16.0} & \textcolor{red}{10.0} &   &   &   &   \\
 	$X^{(1)}$ & 48.0 & 16.0 & 48.0 & 28.0 & \textcolor{red}{8.0} & \textcolor{red}{-8.0} & \textcolor{red}{0.0} & \textcolor{red}{12.0} \\
	$x$ & 56.0 & 40.0 & 8.0 & 24.0 & 48.0 & 48.0 & 40.0 & 16.0 \\ 
	\hline
\end{tabular}
\end{center}
\label{t:IDWT}
\end{table}%

\subsection{Transformada Wavelet para compresión de señales}

Una aplicación de la TDW es la compresión de señales. Para ello se toma como criterio poner en cero los coeficientes con los valores menos significativos. Así por ejemplo dada la señal $x =  [56.0, 40.0, 8.0, 24.0, 48.0, 48.0, 40.0, 16.0  ]$  y su transformada Wavelet es $X = [35,-3,16,10,8,-8,0,12]$ podemos construir dos señales. La primera una señal donde se hacen cero todos los valores con valor absoluto menor o igual a 4 $X_1 = [35, 0,16,10,8,-8,0,12]$ y una segunda señal donde se hacen cero los valores cuyo valor absoluto es menor igual a 9 $X_1 = [35, 0,16,10,0,0,0,12]$. Las antitransformadas de estás señales son $x_1 = [59.0, 43.0,$  $11.0,  27.0,  45.0,  45.0,  37.0,  13.0]$ y $x_2 = [51.0,  51.0,  19.0,$ $  19.0,  45.0,  45.0,  37.0,  13.0] $ respectivamente. En la Figura \ref{fig:wave4} se muestra la señal $x$ en negro y la señal $x_1$ en rojo y de manera similar para la Figura \ref{fig:wave9}.  

\begin{figure}[h!]
\begin{center}
   \subfigure[Señal original y filtrada con $|X(k)| \geq 4$]{\includegraphics[scale=0.35]{Transformadas/figuras/wavelet_4.eps} \label{fig:wave4}}
   \subfigure[Señal original y filtrada con $|X(k)| \geq 9$ ]{\includegraphics[scale=0.35]{Transformadas/figuras/wavelet_9.eps} \label{fig:wave9}}
\end{center}
\caption{Ejemplo de filtrado y compactación de una señal utilizando Transformada Wavelet}
\label{fig:wave1d}
\end{figure} 

\subsection{Transformada Wavelet en dos dimensiones}

La transformada Discreta Wavelet pude realizarse en mas de una dimensión al igual que la convolución. Sin embargo si consideramos kerneles separables en cada dimensión la TDW puede realizarse de manera optima. Así dado un filtro pasa bajas $g$ y un pasa altas $h$ la transformada Wavelet puede llevarse a cabo mediante \eqref{eq:w2_baja}
y \eqref{eq:w2_alta}

\begin{equation}
	y_{low} [n,m] = \sum_{k=\infty}^\infty  \left[ \sum_{l=\infty}^\infty x[2n-k,2m-l] g[l] \right ] g[k] 
	\label{eq:w2_baja}
\end{equation}

\begin{equation}
	y_{high} [n,m] = \sum_{k=\infty}^\infty  \left[ \sum_{l=\infty}^\infty x[2n-k, 2m-l] h[l] \right ] h[k] 
	\label{eq:w2_alta}
\end{equation}

El código en Java correspondiente se muestra a continuación. La función $pwd\_T$ recibe una imagen $I$, el nivel al cual se va aplicar la TDW $p$ y una variable $ising = 1 $ para transformada directa y $ising = -1$ para transformada inversa. La función $pwd2d$ realiza la transformada Wavelet utilizando los kerneles $g = [1/\sqrt{2}, 1/\sqrt{2}]$ y $h = [1/\sqrt{2}, -1/\sqrt{2}]$

\begin{verbatim}
   static public void pwd_T(double I[][], int p, int ising) {
        int i, j, nr = I.length, nc = I[0].length, fact;
        
        Imagen a;
        a = new Imagen(I, "datos");
        System.out.println("Datos " + p + " " + ising);
        
        if(ising > 0 ) fact = 1;
        else fact = 1 << (p-1);
        
        for(i=0; i < p; i++) {
            //System.out.println(fact + " " + p);
            pwd2d(I, nr/fact, nc/fact, ising);
            if(ising > 0) fact *= 2;
            else fact /=2;
        }
    }
    
    static public void pwd2d(double I[][], int nr, int nc, int ising){
        int i, j, N;
        
        N = nr > nc ? nr : nc;
        // Haar Wavelet Transform
        
        double c = 1.0/Math.sqrt(2);
        double cc[] = {c, c};
        double cr[] = new double[2];
        cr[0] = -cc[1]; cr[1] = cc[0];
        
        
        //Daubechies D4 Wavelet Transform

        double aux[] = new double[N];
        
        for(i=0; i<nr; i++) {
            for(j=0; j<nc; j++)
                aux [j] = I[i][j];
            pwt(aux, cc, cr, nc, ising);    
            for(j=0; j<nc; j++)
                I[i][j] = aux [j];
        }
        
        for(j=0; j<nc; j++) {
            for(i=0; i<nr; i++)
                aux [i] = I[i][j];
            pwt(aux, cc, cr, nr, ising);    

            for(i=0; i<nr; i++)
                I[i][j] = aux [i] ;
        }
    }
\end{verbatim}

\begin{verbatim}
    static void pwt(double a[], double cc[], double cr[], int n, int isign) {
        double ai,ai1,wksp[];
        int i, ii, j, k, nh, ncof;
        
        //if (n < 4) return;
        
        ncof = cc.length; // Numero de coeficientes de los kerneles
        //System.out.println("Numero de coeficientes " + ncof);
 
        nh=n >> 1;        // Calcula la mitad del arreglo
        n = nh << 1;      // forza a que n sea par
        wksp = new double[n];
        
        for (j=0;j<n;j++) wksp[j]=0.0;
        
        if (isign >= 0) {
            for (ii=0,i=0;i<n;i+=2,ii++) {
                for (k=0;k<ncof;k++) {
                    j = (i+k)< n ? i+k : i+k-n;
                    wksp[ii] += cc[k]*a[j];
                    wksp[ii+nh] += cr[k]*a[j];
                }
            }
        }
        else {
            for (ii=0,i=0;i<n;i+=2,ii++) {
                ai=a[ii];
                ai1=a[ii+nh];
                for (k=0;k<ncof;k++) {
                    j= (i+k) < n ? i+k : i+k-n;
                    wksp[j] += cc[k]*ai;
                    wksp[j] += cr[k]*ai1;
                }
            }
        }
        
        for (j=0;j<n;j++) a[j]=wksp[j];
    }
\end{verbatim}

\subsection{Ejemplo de compresión de Imágenes}

Una forma de llevar a cabo la compresión de imágenes utilizando la Transformada Wavelet es hacer cero los valores de frecuencia alta. El la Figura \ref{fig:wave2d} se muestra un ejemplo de la imagen de Lena. La primer columna corresponde a la Transformada Wavelest de la imagen, la segunda columna a la mascara utilizada para eliminar frecuencias y la tercer columna la anti transformada Wavelet. Por renglones tenemos diferentes niveles de la transformada Wavelet.  

\begin{figure}[h!]
\begin{center}
   \subfigure[Transformada Wavelet ]{\includegraphics[scale=0.50]{Transformadas/figuras/lenaT.jpg} \label{fig:lenaT}}
\subfigure[Mascara]{\includegraphics[scale=0.50]{Transformadas/figuras/lenaM.jpg} \label{fig:lenaM}}
   \subfigure[Señal comprimida]{\includegraphics[scale=0.50]{Transformadas/figuras/lenaS.jpg} \label{fig:lenaS}}
    \subfigure[Transformada Wavelet ]{\includegraphics[scale=0.50]{Transformadas/figuras/lenaT2.jpg} \label{fig:lenaT2}}
\subfigure[Mascara]{\includegraphics[scale=0.50]{Transformadas/figuras/lenaM2.jpg} \label{fig:lenaM2}}
   \subfigure[Señal comprimida]{\includegraphics[scale=0.50]{Transformadas/figuras/lenaS2.jpg} \label{fig:lenaS2}}
   \subfigure[Transformada Wavelet ]{\includegraphics[scale=0.50]{Transformadas/figuras/lenaT3.jpg} \label{fig:lenaT3}}
\subfigure[Mascara]{\includegraphics[scale=0.50]{Transformadas/figuras/lenaM3.jpg} \label{fig:lenaM3}}
   \subfigure[Señal comprimida]{\includegraphics[scale=0.50]{Transformadas/figuras/lenaS3.jpg} \label{fig:lenaS3}}   
\end{center}
\caption{Ejemplo de filtrado y compactación de una imagen utilizando Transformada Wavelet}
\label{fig:wave2d}
\end{figure} 