{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import math as m\n",
    "import sys\n",
    "import numpy as np\n",
    "import scipy as sp\n",
    "import scipy.io.wavfile as wavfile\n",
    "import scipy.fftpack as spfft\n",
    "import scipy.signal as signal\n",
    "from imageio import imread\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.animation as anim\n",
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "import IPython\n",
    "import IPython.display as ipd\n",
    "import ipywidgets as ipw\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%run helpers.py\n",
    "%run dwts.py"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Exemplo 6.1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 128\n",
    "k = np.linspace(0, N, N)\n",
    "x = 0.5 * np.sin(2*m.pi*3*k/N) + 0.5 * np.sin(2*m.pi*49*k/N)\n",
    "\n",
    "plt.rcParams['figure.figsize'] = [15, 6]\n",
    "plt.plot(k,x)\n",
    "plt.ylim([-1, 1])\n",
    "plt.title(\"Sinal original\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# resposta impulsiva do primeiro filtro de Haar (passa-baixas)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "l = [ 0.5, 0.5 ];\n",
    "# recorta a convolução linear para manter o tamanho do vetor\n",
    "xl = np.convolve(x,l)[0:N]\n",
    "plt.plot(k,xl)\n",
    "plt.ylim([-1, 1])\n",
    "plt.title(\"Sinal filtrado (frequencias baixas)\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# resposta impulsiva do segundo filtro de Haar (passa-altas)\n",
    "h = [ 0.5, -0.5 ]\n",
    "xh = np.convolve(x,h)[0:N]\n",
    "plt.plot(k,xh)\n",
    "plt.ylim([-1, 1])\n",
    "plt.title(\"Sinal filtrado (frequencias altas)\")\n",
    "# Observe que x = xl+xh:\n",
    "print(\"Observe que ||xl+xh-x|| = {}\".format(np.linalg.norm(xh+xl-x)));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Construção dos vetores de aproximação de detalhe"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(xl[0:N:2])\n",
    "plt.ylim([-1, 1])\n",
    "plt.title(\"Coeficientes de aproximacao\")\n",
    "plt.show()\n",
    "\n",
    "plt.plot(xh[1:N:2])\n",
    "plt.ylim([-1, 1])\n",
    "plt.title(\"Coeficientes de detalhe\");\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Exemplo 6.2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def exemplo_62(perfect=False):\n",
    "    \n",
    "    N = 1024\n",
    "    k = np.linspace(0, N, N+1) if perfect else np.linspace(0, N-1, N)\n",
    "    \n",
    "    some1s =  np.ones(N//4)\n",
    "    # sinal com vários trechos independentes\n",
    "    x = np.concatenate( (np.sin(2*m.pi*12* k[0:N//4]/N), 0.5*some1s, 0.1*some1s, -0.2*some1s))\n",
    "    # sentinela para permitir a reconstrução perfeita\n",
    "    if perfect:\n",
    "        x = np.append(x, x[N-1]) \n",
    "    return N, k, x \n",
    "\n",
    "N, k, x = exemplo_62(perfect=True)\n",
    "plt.plot(k,x)\n",
    "plt.ylim([-1, 1])\n",
    "plt.title(\"Sinal original\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# constroi a componente de frequencias baixas\n",
    "l = [ 0.5, 0.5 ]\n",
    "xl = np.convolve(x, l)\n",
    "\n",
    "plt.plot(xl[0:N:2])\n",
    "plt.ylim([-1, 1])\n",
    "plt.title(\"Coeficientes de aproximacao\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# constroi a componente de frequencias altas\n",
    "h = [ 0.5, -0.5 ]\n",
    "xh = np.convolve(x,h)\n",
    "\n",
    "plt.plot(xh[0:N:2])\n",
    "plt.ylim([-1, 1])\n",
    "plt.title(\"Coeficientes de detalhe\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Continuação do exemplo 6.2: compressão"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Reconstrução usando transformada de Haar, zerando o vetor xh:\n",
    "Uxl = np.zeros(N+1)\n",
    "Uxl[0:N+1:2] = xl[0:N+1:2]\n",
    "vl = np.convolve(Uxl,[ 1, 1 ])[1:N+1]\n",
    "\n",
    "plt.plot(vl)\n",
    "plt.axis([0, N, -1, 1])\n",
    "plt.title(\"Compactacao usando a transformada de Haar\")\n",
    "plt.show()\n",
    "print(\"Erro da compactação por Haar (50%) = {0:.2f}\".format( 100*distortion(x[0:N], vl) ))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Reconstrução usando DFT, mantendo 50% dos coeficientes (mais altos)\n",
    "X = np.fft.fft(x[0:N])\n",
    "\n",
    "Xcompactado = np.copy(X)\n",
    "limiar = np.median(abs(X))\n",
    "Xcompactado[abs(Xcompactado) < limiar] = 0\n",
    "\n",
    "xnovo = np.real(np.fft.ifft(Xcompactado))\n",
    "\n",
    "plt.plot(xnovo)\n",
    "plt.axis([0, N, -1, 1])\n",
    "plt.title(\"Compactacao usando DFT\")\n",
    "plt.show()\n",
    "# Calcula erro relativo da reconstrução\n",
    "print(\"Erro da compactação por Fourier (50%) = {0:.4f}\".format(100*distortion(x[0:N], xnovo)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Outra tentativa"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ao invés de anular o vetor xh,\n",
    "# vamos manter 50% dos coeficientes de xh (os mais altos).\n",
    "# A representação terá no total 75% do tamanho de x\n",
    "# (50% correspondente a xl e 25% correspondente a xh).\n",
    "\n",
    "xhnovo = xh.copy()\n",
    "limiar = np.median(abs(xh))\n",
    "xhnovo[abs(xhnovo) < limiar] = 0\n",
    "Uxh = np.zeros(N+1)\n",
    "Uxh[0:N+1:2] = xhnovo[0:N+1:2]\n",
    "vh = np.convolve(Uxh,[ -1, 1 ])[1:N+1]\n",
    "\n",
    "plt.plot(vl+vh);\n",
    "plt.axis([0, N, -1, 1])\n",
    "plt.title(\"Compactacao usando a transformada de Haar\")\n",
    "plt.show()\n",
    "# Calcula erro relativo da reconstrução\n",
    "print(\"Erro da compactação por Haar (75%) = {0:f}\".format(100*distortion(x[0:N], vl+vh)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Reconstrução usando DFT"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# mantendo 75%\n",
    "# dos coeficientes (mais altos)\n",
    "# A expressao statistics(abs(X))(2) corresponde ao\n",
    "# primeiro quartil da distribuição de abs(X).\n",
    "limiar = np.percentile(abs(X).T, 25)\n",
    "Xcompactado = X.copy()\n",
    "Xcompactado[abs(Xcompactado) < limiar] = 0\n",
    "xnovo = np.real(np.fft.ifft(Xcompactado))\n",
    "\n",
    "plt.plot(xnovo)\n",
    "plt.axis([0, N, -1, 1])\n",
    "plt.title(\"Compactacao usando DFT\")\n",
    "plt.show()\n",
    "# Calcula erro relativo da reconstrução\n",
    "print(\"Erro da compactação por Fourier (75%) = {0:.4f}\".format(100*distortion(x[0:N], xnovo)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "N, k, x = exemplo_62()\n",
    "plt.rcParams['figure.figsize'] = [15, 5]\n",
    "plt.plot(k,x)\n",
    "plt.title(\"Sinal original\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "plt.rcParams['figure.figsize'] = [15, 16]\n",
    "nsteps = 4\n",
    "fig, axes = plt.subplots(nsteps,3)\n",
    "for j in range(1, nsteps+1):\n",
    "    y = dwt(x, L=j)\n",
    "    axes[j-1,0].plot( k[0:N//2**j], y[0:N//2**j] )\n",
    "    axes[j-1,0].set_title(\"Coefs de aproximacao de {}a etapa(s)\".format(j))\n",
    "    \n",
    "    axes[j-1,1].plot( k[N//2**j:N//2**(j-1)], y[N//2**j:N//2**(j-1)] )\n",
    "    axes[j-1,1].set_title(\"Coefs de detalhes de {}a etapa(s)\".format(j))\n",
    "    \n",
    "    # zera todos os coeficientes de detalhes de todas as etapas\n",
    "    y[N//2**j:N] = 0.0;\n",
    "    axes[j-1,2].plot( k, idwt(y,j) );\n",
    "    axes[j-1,2].set_title(\"Reconstrucao da {}a etapa(s)\".format(j))\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Exemplo 6.9"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "N, k, x = exemplo_62()\n",
    "plt.rcParams['figure.figsize'] = [15, 5]\n",
    "plt.plot(k,x)\n",
    "plt.title(\"Sinal original\")\n",
    "plt.show()\n",
    "nsteps = 4\n",
    "for j in range(1, nsteps+1):\n",
    "    plt.plot(k, dwt(x,j))\n",
    "    plt.title(\"DWT (Haar) de {} etapa(s)\".format(j))\n",
    "    plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Figuras 6.14 e 6.15"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N, k, x = exemplo_62()\n",
    "plt.rcParams['figure.figsize'] = [15, 5]\n",
    "plt.plot(k,x)\n",
    "plt.title(\"Sinal original\")\n",
    "plt.show()\n",
    "nsteps = 4\n",
    "for j in range(1, nsteps+1):\n",
    "    M = N // (2**j)\n",
    "    coeffilter = np.concatenate( (np.ones(M), np.zeros(N-M)))\n",
    "    # Daubechies \n",
    "    #y = idwt( dwt(x, j, \"Daubechies 4-tap\") * coeffilter, j, \"Daubechies 4-tap\")\n",
    "    # Haar filter\n",
    "    #y = idwt( (dwt(x, j) * coeffilter), j )\n",
    "    # Le Gall 5/3\n",
    "    y = idwt( dwt(x, j, \"Le Gall 5/3\") * coeffilter, j, \"Le Gall 5/3\" )\n",
    "    plt.plot(k,y);\n",
    "    # aqui outra medida de erro?\n",
    "    plt.title(\"Recontrucao a partir da etapa {0:d}: erro = {1:2.2f}%\".format(j, 100*distortion(-x,-y)))\n",
    "    plt.show()\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Exemplo 6.10, figuras 6.16 e 6.17"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "M = rgb2gray(imread(\"double_ferris.jpg\"))[0:896,0:672]\n",
    "\n",
    "plt.rcParams['figure.figsize'] = [9, 9]\n",
    "plt.imshow(M, cmap='gray')\n",
    "plt.title(\"Imagem original\")\n",
    "plt.show()\n",
    "\n",
    "for j in range(1,4):\n",
    "    N = dwt2 (M.copy(), j, \"Le Gall 5/3\")\n",
    "    Nlog = np.log( 1+abs(N) )\n",
    "    Nlog /= np.amax(Nlog)\n",
    "    plt.rcParams['figure.figsize'] = [9, 9]\n",
    "    plt.imshow(Nlog, cmap='gray')\n",
    "    plt.title(\"DWT (Le Gall 5/3) de {} etapa(s)\".format(j))\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Exemplo 6.10, figuras 6.16 e 6.17"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "M = rgb2gray(imread(\"double_ferris.jpg\"))[0:896,0:672]\n",
    "my, mx = M.shape\n",
    "\n",
    "plt.rcParams['figure.figsize'] = [12, 12]\n",
    "plt.imshow(M, cmap='gray')\n",
    "plt.title(\"Imagem original\")\n",
    "plt.show()\n",
    "\n",
    "for j in range(1,4):\n",
    "    N = dwt2(M.copy(), j, \"Le Gall 5/3\")\n",
    "    Mask = np.zeros(N.shape)\n",
    "    Mask[0:my//2, 0:mx//2] += 1\n",
    "    O = idwt2 (N * Mask, j, \"Le Gall 5/3\")\n",
    "    plt.imshow(np.uint8(O), cmap='gray')\n",
    "    plt.title(\"Reconstrucao a partir da etapa {0:d}: erro = {1:2.2f}%\".format(j,100*distortion(-np.float64(M), -O)))\n",
    "    plt.show()\n",
    "    my //= 2\n",
    "    mx //= 2"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}