{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import math as m\n",
    "import sys\n",
    "import time\n",
    "import io\n",
    "from urllib.request import urlopen\n",
    "import numpy as np\n",
    "import scipy as sp\n",
    "import scipy.io.wavfile as wavfile\n",
    "import scipy.signal as sig\n",
    "import soundfile as sf\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"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Exemplo da seção 4.2.1\n",
    "### Resposta do filtro $$w(n) = \\frac{1}{2}x(n)+\\frac{1}{2}x(n-1)$$ às entradas da forma $$x(k) = \\sin(2\\pi qk/N), \\forall k=0,...,\\frac{N}{2}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib notebook\n",
    "N=100\n",
    "k=np.arange(0,N)\n",
    "W=[1]+[0]*(N//2)\n",
    "fig = plt.figure()\n",
    "#plt.ylim([-2, 2])\n",
    "ax1 = fig.add_subplot(121)\n",
    "ax1.set_title(\"entrada\")\n",
    "ax2 = fig.add_subplot(122)\n",
    "ax2.set_title(\"saída\")\n",
    "w = 0*k\n",
    "line1, = ax1.plot(k, w, 'r-')\n",
    "ax1.set_ylim([-1.2,1.2])\n",
    "ax2.set_ylim([-1.2,1.2])\n",
    "line2, = ax2.plot(k, w, 'r-')\n",
    "for q in range(1,N//2+1):\n",
    "  # x(k) é a q-ésima função \"básica\"\n",
    "  x = np.sin(2*m.pi*q*k/N)\n",
    "  w = 0.5*np.roll(x,1)+0.5*x\n",
    "  # Mostra a saída do filtro para cada função básica.\n",
    "  # Observe como a amplitude diminui e a fase inicial\n",
    "  # também varia.\n",
    "  line1.set_ydata(x)\n",
    "  line2.set_ydata(w)\n",
    "  fig.canvas.draw()\n",
    "  #fig.canvas.flush_events()\n",
    "  plt.show()\n",
    "  time.sleep(0.5)\n",
    "  # W contém a resposta em magnitude (ou seja, o fator de\n",
    "  # escala) do filtro para as funções \"básicas\".\n",
    "  W[q] = np.linalg.norm(w)/np.linalg.norm(x)\n",
    "# corrige a última magnitude (é 0 independentemente da fase da entrada)\n",
    "W[N//2] = 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "# Mostra o gráfico da resposta em magnitude \n",
    "plt.figure()\n",
    "plt.plot(W)\n",
    "plt.title(\"relação entre entrada e saída\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Exemplo 4.1 e figuras 4.1 e 4.2\n",
    "# Retoma a função do exemplo na seção 2.3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "T=1\n",
    "N=128\n",
    "t=np.linspace(0,T,N)\n",
    "x = 2*np.cos(2*m.pi*5*t) + 0.8*np.sin(2*m.pi*12*t) + 0.3*np.cos(2*m.pi*47*t)\n",
    "# Calcula a saída do filtro da média.\n",
    "w = 0.5*(x+np.roll(x,1))\n",
    "# Figura 4.1\n",
    "plt.figure()\n",
    "plt.plot(t,w,'-',t,x,'-')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calcula as FFTs do sinal original e do sinal filtrado\n",
    "X = np.fft.fft(x)\n",
    "W = np.fft.fft(w)\n",
    "c = X/N\n",
    "d = W/N\n",
    "E = N*abs(c)**2\n",
    "F = N*abs(d)**2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Figura 4.2\n",
    "plt.figure()\n",
    "plt.plot(E[:N//2+1],'-',F[:N//2+1],'-')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Como é difícil observar as diferenças dos espectros,\n",
    "# as 3 figuras a seguir fazem um zoom dos 3 picos dos\n",
    "# espectros de x (sinal original) e w (sinal filtrado)\n",
    "plt.figure()\n",
    "plt.plot(np.arange(3,8),E[3:8],np.arange(3,8),F[3:8])\n",
    "plt.figure()\n",
    "plt.plot(np.arange(10,15),E[10:15],np.arange(10,15),F[10:15])\n",
    "plt.figure()\n",
    "plt.plot(np.arange(44,53),E[44:53],np.arange(44,53),F[44:53])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Exemplo 4.3 e Figura 4.3:\n",
    "### O filtro da média $w(n) = \\frac{1}{2}x(n)+\\frac{1}{2}x(n-1)$ pode ser\n",
    "### escrito como $w = x*l$, onde $*$ é a convolução. Os\n",
    "### coeficientes do filtro são $l(0)=\\frac{1}{2}$ e $l(1)=\\frac{1}{2}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 128\n",
    "l = np.zeros(N)\n",
    "l[0] = l[1] = 0.5\n",
    "# Mostra os espectros de magnitude e fase de L (Figura 4.3)  \n",
    "L = np.fft.fft(l)\n",
    "plt.figure()\n",
    "plt.plot(np.arange(-N//2+1,N//2+1),np.fft.fftshift(np.abs(L)))\n",
    "plt.figure()\n",
    "plt.plot(np.arange(-N//2+1,N//2+1),np.fft.fftshift(np.angle(L)))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Exemplo 4.4: filtro passa-alta projetado a partir de uma resposta em frequência ideal."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 128\n",
    "# H deverá cortar todas as frequências até k=29 (inclusive)\n",
    "H = np.concatenate((np.zeros(30),np.ones(69),np.zeros(29)))\n",
    "# Calcula a resposta impulsiva do filtro\n",
    "h = np.real(np.fft.ifft(H))\n",
    "# Figura 4.4: resposta impulsiva\n",
    "plt.figure()\n",
    "plt.plot(h)\n",
    "plt.show()\n",
    "# Resposta em frequência do filtro\n",
    "plt.figure()\n",
    "plt.plot(np.arange(-N//2+1,N//2+1),np.fft.fftshift(H))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Versões computacionalmente mais \"baratas\" (com menos taps) do mesmo filtro passa-alta.\n",
    "### Os limiares abaixo são usados para \"cortar\" coeficientes do filtro e com isso diminuir o número de termos na equação de convolução"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "EPSVALS =  [ 0.01, 0.05, 0.1 ]\n",
    "for k in range(len(EPSVALS)):\n",
    "  epsilon = EPSVALS[k]\n",
    "  # seleciona os coeficientes acima do limiar\n",
    "  mask = (abs(h)>=epsilon)\n",
    "  # calcula o número de coeficientes do filtro\n",
    "  ntaps = sum(mask)\n",
    "  # heps é o filtro simplificado com limiar epsilon\n",
    "  heps = h*mask\n",
    "  # Mostra a resposta impulsiva com os coeficientes eliminados\n",
    "  plt.figure()\n",
    "  plt.plot(heps)\n",
    "  plt.title(\"Epsilon = \"+str(epsilon)+\", Filtro com \"+str(ntaps)+\" taps\")\n",
    "  # Compara a resposta do filtro simplificado com a resposta ideal\n",
    "  # do filtro original (pontilhada)\n",
    "  plt.figure()\n",
    "  plt.plot(np.arange(-N//2+1,N//2+1),np.fft.fftshift(np.real(np.fft.fft(heps))),'-',np.arange(-N//2+1,N//2+1),np.fft.fftshift(H),\".\")\n",
    "  plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Colocando em perspectiva a \"perfeição\" da resposta em frequência do filtro original (com banda de corte nula)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Este exemplo (que não está no livro) se contrapõe à idéia de que o filtro desenhado a partir da resposta em frequência H tenha exatamente o comportamento de filtro ideal que o gráfico de H poderia sugerir.\n",
    "\n",
    "#### Para este exemplo, construiremos \"na mão\" a resposta em magnitude do filtro H para funções \"básicas\" do tipo $\\cos(2\\pi f\\frac{n}{N})$, inicialmente para frequências do tipo $f=0,1,...,N-1$, e posteriormente para outras frequências intermediárias."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "FREQS = np.arange(-N//2+1,N//2+1,1)\n",
    "R = [ 0 ] * len(FREQS)\n",
    "t = np.arange(N)\n",
    "for k in range(N):\n",
    "  # cria função básica de frequência k\n",
    "  x = np.cos(2*np.pi*FREQS[k]*t/N)\n",
    "  # filtra por h\n",
    "  y = np.real(np.fft.ifft(np.fft.fft(x)*H))\n",
    "  # constrói gráfico de magnitude da resposta\n",
    "  # (fator de escala da saída do filtro)\n",
    "  R[k] = np.linalg.norm(y)/np.linalg.norm(x)\n",
    "# Mostra o gráfico da resposta de magnitude.\n",
    "plt.figure()\n",
    "plt.plot(FREQS,R)\n",
    "plt.ylim([-0.5, 1.5])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Repete a mesma construção com frequências não-inteiras\n",
    "### Observe que a única diferença do código abaixo é que agora as frequências são f=0,0.1,0.2,...,N-1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "FREQS = np.arange(-N//2+1,N//2+1,0.1)\n",
    "R = [ 0 ] * len(FREQS)\n",
    "t = np.arange(N)\n",
    "for k in range(len(FREQS)):\n",
    "  # cria função básica de frequência k\n",
    "  x = np.cos(2*np.pi*FREQS[k]*t/N)\n",
    "  # filtra por h\n",
    "  y = np.real(np.fft.ifft(np.fft.fft(x)*H))\n",
    "  # constrói gráfico de magnitude da resposta\n",
    "  # (fator de escala da saída do filtro)\n",
    "  R[k] = np.linalg.norm(y)/np.linalg.norm(x)\n",
    "# Mostra o gráfico da resposta de magnitude.\n",
    "plt.figure()\n",
    "plt.plot(FREQS,R)\n",
    "plt.ylim([-0.5, 1.5])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Figura 4.6: Filtro da média 2D\n",
    "### Magnitude da resposta em frequência da \"máscara\" de suavização de 9 pontos"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "M = N = 100\n",
    "h = np.zeros((M,N))\n",
    "# Cria um bloco 3x3 centralizado em (0,0) \n",
    "for j in range(-1,2):\n",
    "  for k in range(-1,2):\n",
    "    h[j%M,k%N] = 1/9\n",
    "# Duas visualizações alternativas para a DFT.\n",
    "H = np.fft.fft2(h)\n",
    "fig = plt.figure()\n",
    "ax = Axes3D(fig)\n",
    "x = np.arange(-M//2+1,M//2+1)\n",
    "y = np.arange(N//2,-N//2,-1)\n",
    "x, y = np.meshgrid(x,y)\n",
    "f = np.abs(np.fft.fftshift(H))\n",
    "ax.plot_surface(x,y,f,cmap=matplotlib.cm.binary_r)\n",
    "plt.figure()\n",
    "plt.imshow(abs(np.fft.fftshift(H)),cmap=matplotlib.cm.binary_r)\n",
    "plt.colorbar()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Figuras 4.7 e 4.8: Uma imagem \"original\", uma versão ruidosa, e duas versões filtradas"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# função auxiliar\n",
    "def rgb2gray(rgb):\n",
    "    fil = [0.299, 0.587, 0.144]\n",
    "    return np.dot(rgb, fil)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "M=rgb2gray(imread(\"double_ferris.jpg\"))\n",
    "S=M.shape\n",
    "# Imagem original\n",
    "plt.figure()\n",
    "plt.imshow(M, cmap='gray')\n",
    "# Versão com ruído adicionado:\n",
    "Mnoisy = M + 20*(np.random.rand(*S)-0.5)\n",
    "plt.figure()\n",
    "plt.imshow(Mnoisy, cmap='gray')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Filtragem usando média de 9 pontos."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# versão 1: calculando as convoluções no domínio do tempo\n",
    "# Define a resposta impulsiva do filtro\n",
    "d = np.zeros((3,3))\n",
    "for j in range(3):\n",
    "  for k in range(3):\n",
    "    d[j,k] = 1/9\n",
    "# O código abaixo faz a filtragem no domínio do tempo.\n",
    "antes = time.clock()\n",
    "Mdenoised = sig.convolve2d(Mnoisy,d,mode='same', boundary='wrap')\n",
    "depois = time.clock()\n",
    "print(\"Tempo da filtragem espacial =\",depois-antes,\"segundos\")\n",
    "plt.figure()\n",
    "plt.imshow(Mdenoised, cmap='gray')\n",
    "plt.show()\n",
    "# Passa a imagem ruidosa 5 vezes no mesmo filtro.\n",
    "Mdenoised10 = Mnoisy.copy()\n",
    "antes = time.clock()\n",
    "for i in range(10):\n",
    "  Mdenoised10 = sig.convolve2d(Mdenoised10,d,mode='same', boundary='wrap')\n",
    "depois = time.clock()\n",
    "print(\"Tempo de 10 filtragens espaciais =\",depois-antes,\"segundos\")\n",
    "plt.figure()\n",
    "plt.imshow(Mdenoised10, cmap='gray')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# versão 2: calculando as convoluções no domínio da frequência\n",
    "# Define a resposta impulsiva do filtro\n",
    "d = np.zeros(S)\n",
    "for j in range(-1,2):\n",
    "  for k in range(-1,2):\n",
    "    d[j%S[0],k%S[1]] = 1/9\n",
    "# O código abaixo faz a filtragem no domínio da frequência.\n",
    "antes = time.clock()\n",
    "Mdenoised = np.real(np.fft.ifft2(np.fft.fft2(Mnoisy)*np.fft.fft2(d)))\n",
    "depois = time.clock()\n",
    "print(\"Tempo da filtragem espectral =\",depois-antes,\"segundos\")\n",
    "plt.figure()\n",
    "plt.imshow(Mdenoised, cmap='gray')\n",
    "plt.show()\n",
    "# Passa a imagem ruidosa 10 vezes no mesmo filtro.\n",
    "antes = time.clock()\n",
    "Mdenoised10 = np.real(np.fft.ifft2(np.fft.fft2(Mnoisy)*np.fft.fft2(d)**10))\n",
    "depois = time.clock()\n",
    "print(\"Tempo de 10 filtragens espectrais =\",depois-antes,\"segundos\")\n",
    "plt.figure()\n",
    "plt.imshow(Mdenoised10, cmap='gray')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Figuras 4.9 e 4.10: detecção de bordas na horizontal e vertical"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### imagem original"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "plt.imshow(M,cmap='gray')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### filtragem para detecção de bordas na vertical"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "v = np.zeros((1,2))\n",
    "v[0,0], v[0,1] = 1, -1\n",
    "Mver = np.abs(sig.convolve2d(M,v,mode='same', boundary='wrap'))\n",
    "maximo = np.amax(Mver)\n",
    "# mostra a imagem filtrada\n",
    "plt.figure()\n",
    "plt.subplot(1,2,1)\n",
    "plt.imshow(1-Mver/maximo, cmap='gray_r')\n",
    "plt.subplot(1,2,2)\n",
    "plt.imshow(1-Mver/maximo, cmap='gray')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### filtragem para detecção de bordas na horizontal"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "h = np.zeros((2,1))\n",
    "h[0,0], h[1,0] = 1, -1\n",
    "Mhor = np.abs(sig.convolve2d(M,h,mode='same', boundary='wrap'))\n",
    "maximo = np.amax(Mhor)\n",
    "# mostra a imagem filtrada\n",
    "plt.figure()\n",
    "plt.subplot(1,2,1)\n",
    "plt.imshow(1-Mhor/maximo,cmap='gray_r')\n",
    "plt.subplot(1,2,2)\n",
    "plt.imshow(1-Mhor/maximo,cmap='gray')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### combina os resultados das filtragens horizontal e vertical"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Mbordas = np.sqrt(Mhor**2+Mver**2)\n",
    "maximo = np.amax(Mbordas)\n",
    "# bordas combinadas\n",
    "plt.figure()\n",
    "plt.subplot(1,2,1)\n",
    "plt.imshow(1-Mbordas/maximo, cmap='gray_r')\n",
    "plt.subplot(1,2,2)\n",
    "plt.imshow(1-Mbordas/maximo, cmap='gray')\n",
    "plt.show()"
   ]
  }
 ],
 "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.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}