{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# MAC0317/MAC5920\n",
    "## Introdução ao Processamento de Sinais Digitais\n",
    "### Seção 2.4: A DFT unidimensional"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "import math as m\n",
    "import io\n",
    "from urllib.request import urlopen\n",
    "import numpy as np\n",
    "import scipy.io.wavfile as wavfile\n",
    "import soundfile as sf\n",
    "import matplotlib.pyplot as plt\n",
    "import IPython.display as ipd"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Seção 2.4.1 Definição da DFT\n",
    "\n",
    "No exemplo (1.23) vimos que um vetor $x\\in\\mathbb{C}^N$ pode ser escrito na base das exponenciais complexas $\\{E_k|k=0,\\ldots,N-1\\}$ como $x=\\displaystyle\\sum_{k=0}^{N-1}c_kE_k$ onde $$c_k=\\frac{(x,E_k)}{(E_k,E_k)}=\\frac{1}{N}(x,E_k).$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Definição 2.1\n",
    "\n",
    "A DFT de um vetor $x\\in\\mathbb{C}^N$ qualquer é o vetor  $X=DFT(x)\\in\\mathbb{C}^N$ cujas componentes são $$X_k = (x,E_k) = \\sum_{n=0}^{N-1}x_ne^{-i2\\pi k\\frac{n}{N}},\\ k=0,\\ldots,N-1.$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Observe que $X_k = Nc_k$ na definição da DFT: isso significa que, normalmente, ao computar a DFT por uma função de biblioteca, por exemplo $\\texttt{numpy.fft.fft(x)}$, os coeficientes $X_k$ são $N$ vezes maiores do que os pesos $c_k$ associados às componentes $E_k$ na combinação linear $x=\\displaystyle\\sum_{k=0}^{N-1}c_kE_k$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Observação 2.1\n",
    "\n",
    "O coeficiente $c_0 = \\frac{X_0}{N}$ mede a contribuição da forma de onda constante $E_0=(1,1,\\ldots,1)$. Em outras palavras, $$c_0 = \\frac{X_0}{N} = \\frac{1}{N}\\sum_{n=0}^{N-1}x_n,$$ ou seja, $c_0$ é o valor médio do sinal $x$, também chamado de componente dc (*direct current*)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Observação 2.2\n",
    "Na expressão da definição 2.4.1, o vetor $X$ pode ser indexado em qualquer componente $k$, mesmo fora do intervalo $0\\le k<N$, e sua expressão é periódica em $k$ com período $N$:\n",
    "$$X_{N+k} = \\sum_{n=0}^{N-1}x_ne^{-i2\\pi(N+k)\\frac{n}{N}} = \\sum_{n=0}^{N-1}x_n\\overbrace{e^{-i2\\pi n}}^{=1}e^{-i2\\pi k\\frac{n}{N}} = X_k,\\ \\forall k\\in\\mathbb{Z}.$$\n",
    "Frequentemente visualizamos as componentes $X_k$ no intervalo $-\\frac{N}{2}<k\\le +\\frac{N}{2}$, que correspondem às frequências que respeitam a condição de Shannon-Nyquist (frequências que não produzem rebatimento)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Definição 2.2\n",
    "\n",
    "A IDFT de um vetor $X\\in\\mathbb{C}^N$ qualquer é o vetor  $x=IDFT(X)\\in\\mathbb{C}^N$ dado pela expressão $$x=\\frac{1}{N}\\sum_{k=0}^{N-1}X_kE_k,$$ de tal forma que suas componentes podem ser escritas como $$x_n = \\frac{1}{N}\\sum_{k=0}^{N-1}X_ke^{i2\\pi k\\frac{n}{N}} = \\frac{1}{N}(X,\\overline{E_n}),\\ n=0,\\ldots,N-1.$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "**Observação:** Existem textos e softwares que podem definir o par DFT/IDFT de maneira diferente. Poderíamos ter definido o par DFT/IDFT pelas equações $c_k=\\displaystyle\\frac{1}{N}\\sum_{n=0}^{N-1}x_ne^{-i2\\pi k\\frac{n}{N}}$ e $x=\\displaystyle\\sum_{k=0}^{N-1}c_kE_k$, onde a constante $\\frac{1}{N}$ apareceria na DFT ao invés de na IDFT, como é mais comum."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Também seria possível definir a DFT/IDFT com a mesma constante nas duas fórmulas:\n",
    "$$\\tilde{X}_k = \\frac{1}{\\sqrt{N}}(x,E_k)$$\n",
    "e\n",
    "$$x = \\frac{1}{\\sqrt{N}}\\sum_{k=0}^{N-1}\\tilde{X}_kE_k.$$\n",
    "Essa versão corresponde a usar a base ortonormal das exponenciais complexas $\\left\\{\\frac{E_k}{\\|E_k\\|}\\mid k=0,\\ldots,N-1\\right\\}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Em qualquer caso, as definições da DFT e da IDFT devem ser \"compatíveis\", ou seja, precisam satisfazer $$IDFT(DFT(x))=x,\\ \\forall x.$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Exemplo de sinal com rebatimento (aliasing)\n",
    "\n",
    "Vamos considerar outra vez o sinal da seção 2.3\n",
    "$$x(t)=2\\cos(2\\pi 5t)+0.8\\sin(2\\pi 12 t)+0.3\\cos(2\\pi 47t)$$\n",
    "amostrado no intervalo $t\\in[0,1)$ usando duas taxas de amostragem diferentes: $N=128$ e $N=64$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "N = 128; t = np.arange(0, 1, 1/N); N2 = N/2; t_sub = np.arange(0, 1, 1/N2)\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",
    "x_sub = 2*np.cos(2* m.pi *5*t_sub) + 0.8*np.sin(2*m.pi*12*t_sub) + 0.3*np.cos(2*m.pi*47*t_sub)\n",
    "f, ax = plt.subplots(1,2,figsize=(15,5))\n",
    "ax[0].set_title(\"Sinal original\");ax[0].plot(t,x)\n",
    "ax[1].plot(t_sub,x_sub); ax[1].set_title(\"Sinal amostrado a 64 Hz\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Espectro de energia das duas versões do sinal\n",
    "\n",
    "Lembrando da expressão $\\|x\\|^2 = \\displaystyle\\sum_{k=0}^{N-1}N|c_k|^2$, e do fato de que as duas componentes $E_k$ e $E_{-k}$ juntas representam uma componente real (seno ou cosseno), vamos usar como medida de energia os valores\n",
    "$$N(|c_{-k}|^2+|c_k|^2),$$\n",
    "ou equivalentemente, lembrando que $c_k=\\frac{X_k}{N}$,\n",
    "$$\\frac{|X_{-k}|^2+|X_k|^2}{N}.$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "X= np.fft.fft(x); c = X/ (N); E = N * abs(c)**2\n",
    "X_sub= np.fft.fft(x_sub); c_sub = X_sub / (N/2); E_sub = (N/2) * abs(c_sub)**2\n",
    "f, ax = plt.subplots(1,2,figsize=(15,5))\n",
    "ax[0].plot(E[0:int(N/2)]); ax[0].set_title(\"Energia do sinal original\");\n",
    "ax[1].plot(E_sub[0:int(N/4)]); ax[1].set_title(\"Energia do sinal reamostrado (Figura 2.7)\");\n",
    "plt.show(); f = 47; 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": [
    "Nesse exemplo o rebatimento da componente 47 Hz em uma componente de 17 Hz ocorre pois, usando uma taxa de amostragem de 64 Hz,\n",
    "$$\\begin{array}{ll}\n",
    "\\cos\\left(2\\pi 47\\frac{n}{64}\\right) \n",
    "& = \\cos\\left(2\\pi 47\\frac{n}{64} - 2\\pi n\\frac{64}{64}\\right)\\\\\n",
    "& = \\cos\\left(2\\pi(47-64)\\frac{n}{64}\\right)\\\\\n",
    "& = \\cos\\left(2\\pi(-17)\\frac{n}{64}\\right)\\\\\n",
    "& = \\cos\\left(2\\pi 17\\frac{n}{64}\\right),\\ \\forall n\\in\\mathbb{Z},\n",
    "\\end{array}$$\n",
    "ou seja, as componentes de 47 Hz e de 17 Hz são indistinguíveis nessa taxa de amostragem."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Exemplo com pulsos quadrados (Figuras 2.8 - 2.10)\n",
    "\n",
    "Considere um pulso quadrado de largura $R$ e duração $N$, ou seja, um sinal $x\\in\\mathbb{R}^N$ tal que\n",
    "$$x_n = \\left\\{\\begin{array}{ll}\n",
    "1&n<R\\\\\n",
    "0&n\\ge R,\n",
    "\\end{array}\\right.$$\n",
    "cuja DFT é (exercício 2.11 resolvido):\n",
    "$$X_k = \\sum_{n=0}^{N-1}x_ne^{-i2\\pi k\\frac{n}{N}} = \\sum_{n=0}^{R-1}(\\overbrace{e^{-i2\\pi k\\frac{1}{N}}}^z)^n = \\frac{1-z^R}{1-z} = \\frac{1-e^{-i2\\pi k\\frac{R}{N}}}{1-e^{-i2\\pi k\\frac{1}{N}}}.$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Cada componente $X_k = \\frac{1-e^{-i2\\pi k\\frac{R}{N}}}{1-e^{-i2\\pi k\\frac{1}{N}}}$ possui energia\n",
    "$$\\begin{array}{ll}\n",
    "|X_k|^2\n",
    "& = \\frac{|1-e^{-i2\\pi k\\frac{R}{N}}|}{|1-e^{-i2\\pi k\\frac{1}{N}}|}\n",
    "= \\frac{\\left(1-\\cos\\left(2\\pi k\\frac{R}{N}\\right)\\right)^2+\\left(\\sin\\left(2\\pi k\\frac{R}{N}\\right)\\right)^2}{\\left(1-\\cos\\left(2\\pi k\\frac{1}{N}\\right)\\right)^2+\\left(\\sin\\left(2\\pi k\\frac{1}{N}\\right)\\right)^2}\\\\\n",
    "& = \\frac{1-2\\cos\\left(2\\pi k\\frac{R}{N}\\right)+\\left(\\cos\\left(2\\pi k\\frac{R}{N}\\right)\\right)^2+\\left(\\sin\\left(2\\pi k\\frac{R}{N}\\right)\\right)^2}{1-2\\cos\\left(2\\pi k\\frac{1}{N}\\right)+\\left(\\cos\\left(2\\pi k\\frac{1}{N}\\right)\\right)^2+\\left(\\sin\\left(2\\pi k\\frac{1}{N}\\right)\\right)^2}\\\\\n",
    "& = \\frac{2-2\\cos\\left(2\\pi k\\frac{R}{N}\\right)}{2-2\\cos\\left(2\\pi k\\frac{1}{N}\\right)}\\\\\n",
    "& = \\frac{1-\\cos\\left(2\\pi k\\frac{R}{N}\\right)}{1-\\cos\\left(2\\pi k\\frac{1}{N}\\right)}.\n",
    "\\end{array}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "def pulso(R,title):\n",
    "    N = 256; t = np.arange(N);\n",
    "    f, ax = plt.subplots(len(R),2,figsize=(15,10)); f.subplots_adjust(hspace=0.35)\n",
    "    for i in range(len(R)):\n",
    "        r = R[i]; x = np.zeros(N); x[0:r:1] = 1; X = np.fft.fft(x); c = X / N; E = N * abs(c)**2\n",
    "        ax[i,0].set_title(f\"Pulso Quadrado com R = {r}\"); ax[i,0].plot(x)\n",
    "        ax[i,1].set_title(f\"FFT do pulso quadrado com R = {r}\"); ax[i,1].plot(np.arange(N) - N/2, np.roll(E, int(N/2)))\n",
    "        f.suptitle(title)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "pulso([50,20],\"Figura 2.8\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "pulso([10,5],\"Figura 2.9\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "pulso([2,1],\"Figura 2.10\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Exemplo de sinal aleatório e sua DFT (Figura 2.11)\n",
    "\n",
    "Considere $x\\in\\mathbb{R}^N$ obtido por $N$ sorteios aleatórios de uma mesma distribuição normal, ou seja,\n",
    "$$x_n = \\mathcal{N}(0,\\sigma^2),$$\n",
    "onde $\\mu=0$ é a média e $\\sigma^2$ é a variância da distribuição normal."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "É possível provar (exercício 2.12) que as componentes de Fourier $X_k$, sendo elas próprias combinações lineares das amostras aleatórias $x_n$, também possuem distribuição normal, sendo \n",
    "$$X_0 = \\mathcal{N}(0,N\\sigma^2)$$\n",
    "e\n",
    "$$X_k = \\mathcal{N}\\left(0,\\frac{N}{2}\\sigma^2\\right).$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "N = 256; t = np.arange(N); x = np.random.rand(N) * 2 - 1\n",
    "X = np.fft.fft(x); c = X / N; E = N * abs(c)**2\n",
    "f, ax = plt.subplots(1,2,figsize=(15,5))\n",
    "ax[0].set_title(\"Sinal aleatório\"); ax[0].plot(t, x)\n",
    "ax[1].set_title(\"Energia da DFT\"); ax[1].plot(np.arange(N) - N/2, np.roll(E, int(N/2))); plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Sugestões para plotar DFTs\n",
    "\n",
    "- É usual plotarmos $|X_k|$ (magnitude) ou $|X_k|^2$ (energia), ou ainda alguma versão com mudança de escala ($|c_k|=|X_k|/N$ ou $N|c_k|^2 = |X_k|^2/N$). Apenas em aplicações específicas estaremos interessados no espectro de fase ($\\measuredangle X_k$)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "- Quando os valores de $X_k$ variam em uma faixa muito grande de valores, ou quando são muito baixos, é comum plotarmos espectros em escala logarítmica ($\\log|X_k|$, ou $\\log(1+|X_k|)$ para evitarmos $\\log 0 = -\\infty$)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "- Escolha a faixa de frequências mais conveniente para o que se quer ilustrar. Frequentemente estaremos interessados em combinar as componentes $-k$ e $k$ para observar energias associadas a componentes reais (senos/cossenos), através de expressões como $(|X_{-k}|^2+|X_k|^2)/N$, para $k=1,\\ldots,\\frac{N}{2}-1$ (observe que as componentes $k=0$ e $k=\\frac{N}{2}$ devem ser consideradas separadamente)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Observação 2.3: Índices da DFT e Frequências em Hz\n",
    "\n",
    "Os índices $k$ das componentes $X_k$ da DFT de um sinal $x\\in\\mathbb{C}^N$ representam frequências medidas em **número de ciclos completos no intervalo $[0,N)$**. Para traduzirmos esses índices em frequências em Hz, precisamos considerar **quantos segundos** dura o sinal $x$, através da expressão $T = \\frac{N}{R}$. Uma componente que percorre $k$ ciclos em $N$ amostras percorrerá $k$ ciclos em $T$ segundos, o que corresponde a uma frequência de $\\frac{k}{T} = \\frac{kR}{N}$ Hz."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "Observe que os índices $k=0,1,\\ldots,\\frac{N}{2}$ corresponderão dessa forma às frequências\n",
    "$$\\frac{0R}{N}=0,\\frac{1R}{N}=\\frac{R}{N},\\ldots,\\frac{NR}{2N}=\\frac{R}{2}\\quad\\mbox{Hz}.$$\n",
    "Os índices $k=\\frac{N}{2}+1,\\ldots,N-1$, que equivalem (por rebatimento) aos índices $k=-\\frac{N}{2}+1,\\ldots,-1$, correspondem às frequências\n",
    "$$\\frac{(-N/2+1)R}{N}=-\\frac{R}{2}+\\frac{R}{N},\\ldots,\\frac{-1R}{N}=-\\frac{R}{N}\\quad\\mbox{Hz}.$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Exemplo de áudio: Sinal de sino e sua DFT (Figura 2.12)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "url=\"https://www.ime.usp.br/~mqz/bells004.wav\" # abre arquivo wav de url:\n",
    "x, rate = sf.read(io.BytesIO(urlopen(url).read()))\n",
    "x_mono = x[:,0] + x[:,1] # converte sinal estereo para mono\n",
    "x_mono = x_mono-np.mean(x_mono,0) # corrige média do sinal\n",
    "plt.figure(figsize=(15,5)); plt.plot(x_mono); 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_mono); N = len(x); E = 2*abs(X[0:int(N/2):1])**2/N\n",
    "f, ax = plt.subplots(1,2,figsize=(15,5))\n",
    "ax[0].set_title(\"Espectro de Energia\"); ax[0].plot(np.arange(N/2) *rate/N, E)\n",
    "ax[1].set_title(\"Espectro de Energia (detalhe)\"); ax[1].plot(np.arange(6000) *rate/N, E[0:6000:1]); 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": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "Y = X.copy(); 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)); EY = 2*abs(Y[0:int(N/2):1])**2/N\n",
    "f, ax = plt.subplots(1,2,figsize=(15,5))\n",
    "ax[0].set_title(\"Sinal ressintetizado\"); ax[0].plot(rex)\n",
    "ax[1].set_title(\"Espectro ressintetizado\"); ax[1].plot(np.arange(N/2) *rate/N, EY); plt.show()\n",
    "ipd.Audio(rex, rate=rate)"
   ]
  }
 ],
 "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.7.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}