{
 "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": {},
   "outputs": [],
   "source": [
    "N = 2**8; t = np.arange(N)\n",
    "a = -1; b = 1\n",
    "x = np.random.rand(N) * (b - a) + a  #amostra valores de uma distribuição uniforme em [a, b)\n",
    "sigma = (b - a) / m.sqrt(12)   #desvio padrão da variável aleatória com distribuição uniforme em [a,b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "#referência: Time Series Analysis and Its Applications: With R Examples (Robert H. Shumway, David S. Stoffer)\n",
    "plt.subplot(4, 2, 1); plt.title(\"Sinal temporal\"); plt.plot(t, x)\n",
    "shifted_range = np.arange(N) - N/2\n",
    "\n",
    "#uma sequência de variáveis aleatórias descorrelacionadas com variância sigma^2 terá densidade espectral constante igual a sigma^2 (ex 4.4 do Sumway & Stoffer)\n",
    "plt.subplot(4, 2, 2); plt.title(\"Densidade espectral\"); plt.plot(shifted_range, (sigma ** 2) * np.ones(N) )\n",
    "\n",
    "X = np.fft.fft(x)\n",
    "c = X / m.sqrt(N) # coeficentes da DFT (eq. 4.18 do Shumway & Stoffer)\n",
    "E = abs(c)**2  #peridograma (energia) do sinal  (eq. 4.20 do Shumway & Stoffer): seria o mesmo obtido fazendo c=X/N; E = N*abs(c)**2\n",
    "shifted_E = np.roll(E, int(N/2)) #faz um \"shift à direita\" de  int(N/2) posições em E\n",
    "\n",
    "#Ponto positivo:quando N tende a infinito, a esperança (valor médio) do periodograma tende à densidade espectral (eq. 4.20 do Shumway & Stoffer)\n",
    "#-> Logo, o periodograma é um estimador assintoticamente não viesado da densidade espectral\n",
    "#Problema: a variância do periodograma não diminui conforme N aumenta (ver prop. 4.4 e exemplo 4.10 do Shumway & Stoffer)\n",
    "#-> Logo, o periodograma não é um estimador assintoticamente \"estável\" densidade espectral\n",
    "plt.subplot(4, 2, 3); plt.title(\"Periodograma (energia)\"); plt.plot(shifted_range, shifted_E)\n",
    "\n",
    "#Alternativa: tentar diminuir a variância através da suavização do periodograma (ver seção 4.5 do Shumway & Stoffer)\n",
    "#-> Uma opção é usar médias móveis, com uma janela deslizante de largura L \n",
    "def smooth(a, L=3, kernel='daniell', repeat=False):\n",
    "    if L % 2 != 1 or L <= 1 or not isinstance(L, int): #L deve ser ímpar, inteiro e >= 3\n",
    "        raise ValueError(\"Parameter 'L' (%d) must be an odd integer greater or equal to 3\" % (L))  \n",
    "    if L > len(a):\n",
    "        raise ValueError(\"Parameter 'L' (%d) is too large to apply kernel again\" % (L))              \n",
    "        \n",
    "    weights = False #pesos da media\n",
    "    output = False  #saida (sinal suavizado)\n",
    "    mode = 'valid' #faz a media ponderada apenas quando a janela estiverem totalmente contida no sinal \n",
    "    if kernel == 'daniell':\n",
    "        weights = np.ones(L)/L #vetor de pesos de tamanho L, com cada elemento igual a 1/L \n",
    "    elif kernel == 'modified.daniell':\n",
    "        w1 = 0.5 / (L-1) #peso das pontas\n",
    "        w2 = 2 * w1 #peso no interior da janela\n",
    "        weights = np.concatenate([[w1],np.ones(L-2)*w2,[w1]])\n",
    "    else:\n",
    "        raise ValueError(\"Parameter 'kernel' (%s) must be '' or ''\" % (kernel))  \n",
    "    \n",
    "    output = np.convolve(a, weights, 'valid')\n",
    "    \n",
    "    if repeat:\n",
    "        if L > len(output):\n",
    "            raise ValueError(\"Parameter 'L' (%d) is too large to apply the kernel again\" % (L))              \n",
    "        output = np.convolve(output, weights, 'valid')\n",
    "        \n",
    "    return output\n",
    "\n",
    "L = 3\n",
    "smoothed_E = smooth(shifted_E,L,kernel='daniell', repeat=False)\n",
    "plt.subplot(4, 2, 4); plt.title(\"Per. suavizado (Daniell)\"); plt.plot(shifted_range[int(np.floor(L/2)):-int(np.floor(L/2))], smoothed_E)\n",
    "\n",
    "smoothed_E = smooth(shifted_E,L,kernel='modified.daniell', repeat=False)\n",
    "plt.subplot(4, 2, 5); plt.title(\"Per. suavizado (Daniell modif.)\"); plt.plot(shifted_range[int(np.floor(L/2)):-int(np.floor(L/2))], smoothed_E)\n",
    "\n",
    "smoothed_E = smooth(shifted_E,L,kernel='daniell', repeat=True)\n",
    "plt.subplot(4, 2, 6); plt.title(\"Per. suavizado (Dan. repetido)\"); plt.plot(shifted_range[int(2*np.floor(L/2)):-int(2*np.floor(L/2))], smoothed_E)\n",
    "\n",
    "smoothed_E = smooth(shifted_E,L,kernel='modified.daniell', repeat=True)\n",
    "plt.subplot(4, 2, 7); plt.title(\"Per. suavizado (Dan. modif. repetido)\"); plt.plot(shifted_range[int(2*np.floor(L/2)):-int(2*np.floor(L/2))], smoothed_E)\n",
    "\n",
    "plt.tight_layout()\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
}