{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "import math as m\n",
    "import sys\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 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": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Seção 2.3: Um exemplo motivacional"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Figura 2.2 - Sinal analógico no domínio do tempo"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "T = 1\n",
    "N = 128\n",
    "t = np.arange(0, T, 1/N)\n",
    "doispi = 2 * m.pi\n",
    "x = 2*np.cos(doispi*5*t) + 0.8*np.sin(doispi*12*t) + 0.3*np.cos(doispi*47*t)\n",
    "print(\"Primeiras 10 amostras do sinal: \",np.round(x[0:10],3))\n",
    "plt.plot(t,x); plt.xlabel(\"t\"); plt.ylabel(\"x(t)\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "X = np.fft.fft(x)\n",
    "c = X / N\n",
    "print(np.round(c,2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Mostra coeficientes não-nulos"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "eps = 1e-8\n",
    "for k, ck in np.ndenumerate(c):\n",
    "    if abs(ck) > eps: print (\"c({0}) = {1:.2f}\".format(k,ck)) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Mostra energia em cada componente"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "print(\"E(x)\\t=\\t{0:.2f}\".format(np.linalg.norm(x)**2))\n",
    "print(\"--------------------------\")\n",
    "E = N * abs(c)**2\n",
    "for (k,), Ek in np.ndenumerate(E):\n",
    "    if abs(Ek) > eps: print(\"E({0})\\t=\\t{1:.2f}\".format(k,Ek))\n",
    "print(\"--------------------------\")\n",
    "print(\"E(c)*N\\t=\\t{0:.2f}\".format(sum(E)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Figuras 2.3 e 2.4 - Espectro do sinal de exemplo"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "### Espectro centralizado na frequência 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "plt.plot(np.arange(N) - N/2, np.roll(E, int(N/2)))\n",
    "plt.xlabel(\"k\"); plt.ylabel(\"E(k)\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Espectro no intervalo [0, N-1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "plt.plot(E); plt.xlabel(\"k\"); plt.ylabel(\"E(k)\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Figura 2.5 - Espectro com pares de frequências $\\pm\\ f_k$ agrupadas"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "E[1:int(N/2)-1:1] *= 2\n",
    "plt.plot(E[0:int(N/2)]); plt.xlabel(\"k\"); plt.ylabel(\"E(k)+E(-k)\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Figura 2.6 - Sinal ruidoso e sinal filtrado"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "# Elimina \"altas frequências\" (acima de 40 Hz)\n",
    "Y = X\n",
    "Y[40:N-40:1] *= 0\n",
    "y = np.real(np.fft.ifft(Y))\n",
    "plt.subplot(1, 2, 1);plt.title(\"Sinal original\");plt.plot(t,x)\n",
    "plt.subplot(1, 2, 2);plt.title(\"Sinal filtrado\");plt.plot(t,y)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Seção 2.4: A DFT unidimensional"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Figura 2.7 - Espectro do sinal com rebatimento (aliasing)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "N2 = N/2\n",
    "tt = np.arange(0, 1, 1/N2)\n",
    "xx = 2*np.cos(doispi*5*tt) + 0.8*np.sin(doispi*12*tt) + 0.3*np.cos(doispi*47*tt)\n",
    "plt.plot(tt,xx); plt.title(\"Sinal amostrado a 64 Hz\"); plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "XX = np.fft.fft(xx)\n",
    "cc = XX / (N2)\n",
    "EE = (N/2) * abs(cc)**2\n",
    "plt.plot(EE[0:int(N/4)]); plt.title(\"Energia do sinal reamostrado\"); plt.show()\n",
    "f = 47\n",
    "print(\"Componente de {0} Hz rebatida = {1} Hz\".format(f,f%N2 if f%N2<=N2/2 else N2-f%N2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Figuras 2.8 - 2.10 - Pulsos quadrados"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "N = 256\n",
    "t = np.arange(N)\n",
    "R = [50, 10, 2 ,1]\n",
    "for r in R:\n",
    "    x = np.zeros(N)\n",
    "    x[0:r:1] = 1\n",
    "    plt.subplot(1, 2, 1)\n",
    "    plt.plot(x)\n",
    "    X = np.fft.fft(x)\n",
    "    c = X / N\n",
    "    E = N * abs(c)**2\n",
    "    plt.subplot(1, 2, 2)\n",
    "    plt.plot(np.arange(N) - N/2, np.roll(E, int(N/2)))\n",
    "    plt.show()\n",
    "fig = plt.figure()\n",
    "plt.plot()\n",
    "plt.axis('off')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Figura 2.11 - Sinal aleatório e sua DFT"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "N = 256; t = np.arange(N)\n",
    "x = np.random.rand(N) * 2 - 1\n",
    "plt.subplot(1, 2, 1); plt.title(\"Sinal temporal\"); plt.plot(t, x)\n",
    "X = np.fft.fft(x)\n",
    "c = X / N; E = N * abs(c)**2\n",
    "plt.subplot(1, 2, 2); plt.title(\"Energia da DFT\"); plt.plot(np.arange(N) - N/2, np.roll(E, int(N/2)))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Figura 2.12: Sinal de sino e sua DFT"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "# para abrir arquivo wav de url:\n",
    "url=\"http://static1.grsites.com/archive/sounds/bells/bells004.wav\"\n",
    "x, rate = sf.read(io.BytesIO(urlopen(url).read()))\n",
    "# esse arquivo tem uma codificação estranha, com uma média muito maior do que 0: corrige\n",
    "x = x-np.mean(x,0)\n",
    "# para abrir arquivo wav local:\n",
    "# rate, x = wavfile.read('bell.wav')\n",
    "plt.plot(x)\n",
    "plt.show()\n",
    "ipd.Audio(x.T, rate=rate)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "X = np.fft.fft(x[:,0])\n",
    "N = len(x)\n",
    "E = 2*abs(X[0:int(N/2):1])**2/N\n",
    "plt.subplot(1, 2, 1); plt.plot(np.arange(N/2) *rate/N, E); plt.title(\"Espectro de Energia\")\n",
    "plt.subplot(1, 2, 2); plt.plot(np.arange(6000) *rate/N, E[0:6000:1]); plt.title(\"Espectro de Energia (detalhe)\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Exemplo de remoção de componentes mais \"fracas\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "Y = X.copy()\n",
    "C = 0.1 * abs(Y).max()\n",
    "for val in np.nditer(Y, op_flags=['readwrite']):\n",
    "    if abs(val) < C:\n",
    "        val[...] = 0\n",
    "rex = np.real(np.fft.ifft(Y))\n",
    "EY = 2*abs(Y[0:int(N/2):1])**2/N\n",
    "plt.subplot(1, 2, 1); plt.plot(rex); plt.title(\"Sinal ressintetizado\")\n",
    "plt.subplot(1, 2, 2); plt.plot(np.arange(N/2) *rate/N, EY); plt.title(\"Espectro ressintetizado\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "ipd.Audio(rex, rate=rate)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Seção 2.5: Propriedades da DFT"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Seção 2.6: Transformada Rápida de Fourier"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Seção 2.7: A DFT bidimensional"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "# Figura 2.13 - Imagem e várias apresentações da DFT"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "url = \"http://sutherncharm.files.wordpress.com/2009/09/double-ferris.jpg\"\n",
    "M = imread(urlopen(url).read())\n",
    "plt.subplot(1, 2, 1); plt.title(\"Imagem original\"); plt.imshow(M)\n",
    "N = np.fft.fft2(M,axes=(0,1)); absN = np.abs(N)\n",
    "plt.subplot(1, 2, 2); plt.title(\"Espectro\"); plt.imshow(255 * absN / absN.max())\n",
    "plt.show()\n",
    "halfheight, halfwidth = int(absN.shape[0]/2), int(absN.shape[1]/2)\n",
    "shiftN = np.roll(absN, (halfheight,halfwidth), axis=(0,1))\n",
    "plt.subplot(1, 2, 1); plt.title(\"Espectro centralizado em (0,0)\"); plt.imshow(255 * shiftN / shiftN.max() )\n",
    "Nlog = np.log(1 + shiftN)\n",
    "maxi, mini = Nlog.max(), Nlog.min(); Nshow = (Nlog - mini) / (maxi - mini)\n",
    "plt.subplot(1, 2, 2); plt.title(\"Log-Espectro normalizado\"); plt.imshow(Nshow)\n",
    "plt.show()\n",
    "C = 0.001 * absN.max()\n",
    "for val in np.nditer(N, op_flags=['readwrite']):\n",
    "    if abs(val) < C:\n",
    "        val[...] = 0\n",
    "reM = np.round(np.real(np.fft.ifft2(N, axes=(0,1))))\n",
    "reMin, reMax = reM.min(), reM.max()\n",
    "reMshow = (reM - reM.min())/ (reM.max() - reM.min())\n",
    "plt.title(\"Ressíntese após compressão\"); plt.imshow(reMshow)\n",
    "plt.show()\n",
    "fig = plt.figure()\n",
    "plt.plot()\n",
    "plt.axis('off')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "# Figura 2.14: Imagem artificial e visualizações da DFT"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "O = np.zeros((300,300,3)); O[100:150,100:110,0:3] += 1\n",
    "plt.subplot(1, 2, 1); plt.title(\"Imagem original\"); plt.imshow(O)\n",
    "\n",
    "N = np.fft.fft2(O,axes=(0,1)); absN = abs(N)\n",
    "plt.subplot(1, 2, 2); plt.title(\"Espectro\"); plt.imshow(absN/ absN.max())\n",
    "plt.show()\n",
    "\n",
    "halfheigth, halfwidth = int(absN.shape[0]/2), int(absN.shape[1]/2)\n",
    "shiftN = np.roll(absN, (halfheigth, halfwidth), axis=(0,1))\n",
    "plt.subplot(1, 2, 1); plt.title(\"Espectro centralizado em (0,0)\"); plt.imshow(shiftN / shiftN.max())\n",
    "\n",
    "Nlog = np.log(1 + shiftN)\n",
    "maxi, mini = Nlog.max(), Nlog.min()\n",
    "Nshow = (Nlog - mini) / (maxi - mini)\n",
    "plt.subplot(1, 2, 2); plt.title(\"Log-Espectro normalizado\"); plt.imshow(Nshow)\n",
    "plt.show()\n",
    "\n",
    "fig = plt.figure()\n",
    "plt.plot()\n",
    "plt.axis('off')\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Slideshow",
  "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
}