{
 "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"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Exemplo 5.1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 1000\n",
    "t = np.linspace(0, 1-1/N, N)\n",
    "f = 0.5*np.sin(2*m.pi*96*t)+0.5*np.sin(2*m.pi*235*t)\n",
    "g = np.concatenate( (np.sin(2*m.pi*96*t[0:N//2]), np.sin(2*m.pi*235*t[N//2:N])) , axis=0) \n",
    "\n",
    "plt.rcParams['figure.figsize'] = [20, 12]\n",
    "fig, axes = plt.subplots(2,2)\n",
    "axes[0,0].set_title('Espectro da função f')\n",
    "axes[0,0].plot( range(N//2) , abs(np.fft.fft(f)[0:N//2]) )\n",
    "axes[0,0].axis([0, 500, 0, 300])\n",
    "axes[0,1].set_title('Função f')\n",
    "axes[0,1].plot(f)\n",
    "\n",
    "axes[1,0].set_title('Espectro da função g')\n",
    "axes[1,0].plot( range(N//2) , abs(np.fft.fft(g)[0:N//2]) )\n",
    "axes[1,0].axis([0, 500, 0, 300])\n",
    "axes[1,1].set_title('Função g')\n",
    "axes[1,1].plot(g)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Exemplo 5.2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# mesmo N, t, f\n",
    "w = np.zeros(N)\n",
    "w[99:149] = np.ones(50)\n",
    "ftil = w * f\n",
    "\n",
    "plt.rcParams['figure.figsize'] = [15, 12]\n",
    "fig, axes = plt.subplots(2,1)\n",
    "axes[0].plot( range(N//2), abs(np.fft.fft(ftil)[0:N//2]) )\n",
    "axes[0].axis([0, 500, 0, 14])\n",
    "axes[0].set_title(\"Espectro do sinal multiplicado pela janela\")\n",
    "\n",
    "axes[1].plot( range(25), abs(np.fft.fft(ftil[99:149])[0:25]))\n",
    "axes[1].axis([0, 25, 0, 14])\n",
    "axes[1].set_title(\"Espectro do sinal recortado\");\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Exemplo 5.3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.rcParams['figure.figsize'] = [15, 6]\n",
    "fig, axes = plt.subplots()\n",
    "axes.plot( range(-N//2, N//2), np.roll(abs(np.fft.fft(w)), N//2 ))\n",
    "axes.set_title(\"Espectro da janela (retangular) do exemplo 5.2\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Exemplo 5.4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.rcParams['figure.figsize'] = [15, 18]\n",
    "M = np.array([3, 10, 100, 250])\n",
    "fig, axes = plt.subplots(len(M), 1)\n",
    "\n",
    "for j in range(len(M)):\n",
    "    w = np.concatenate( (np.ones(M[j]), np.zeros(N-M[j])) )\n",
    "    axes[j].plot(range(-N//2, N//2), np.roll(abs(np.fft.fft(w)), N//2))\n",
    "\n",
    "    axes[j].set_xlim([-N/2+1, N/2])\n",
    "    axes[j].set_title(\"Espectro da janela retangular com M={}\".format(M[j]))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Exemplo adicional: Figura 5.2 com diferentes tamanhos de janela"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.rcParams['figure.figsize'] = [15, 18]\n",
    "M = np.array([3, 10, 100, 250])\n",
    "fig, axes = plt.subplots(len(M), 2)\n",
    "\n",
    "for j in range(len(M)):\n",
    "    w = np.concatenate( (np.ones(M[j]), np.zeros(N-M[j])) )\n",
    "    ftil = w * f\n",
    "    fcut = ftil[0:M[j]]\n",
    "    axes[j,0].plot(range(N//2+1), abs(np.fft.fft(ftil)[0:N//2+1]))\n",
    "    axes[j,0].set_xlim([0, N//2])\n",
    "    axes[j,0].set_title(\"Espectro de w*f com M={}\".format(M[j]))\n",
    "    axes[j,1].plot(range(M[j]//2+1), abs(np.fft.fft(fcut)[0:M[j]//2+1]))\n",
    "    axes[j,1].set_xlim([0, M[j]//2])\n",
    "    axes[j,1].set_title(\"Espectro do sinal recortado com M={}\".format(M[j]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Exemplo 5.5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def show_spectra(N, t, g, M, O):\n",
    "    for j in range(M.size):\n",
    "        Snrows = M[j] // 2\n",
    "        Sncols = m.floor( (N-M[j]) / O[j] ) + 1\n",
    "        S = np.zeros( (Snrows, Sncols) )\n",
    "        for k in range(Sncols):\n",
    "            inicio = int(k * O[j])\n",
    "            fim = int(k * O[j] + M[j])\n",
    "            magspec = abs(np.fft.fft(g[inicio:fim]))\n",
    "            S[0:Snrows, k] = magspec[0:Snrows]\n",
    "        S = np.flipud(S)\n",
    "        S = S / np.amax(S)\n",
    "        plt.imshow(S, cmap='gray', aspect='auto', interpolation='none', extent=[0, 1, 0, 500])\n",
    "        plt.title(\"Espectrograma com M={}, m={}\".format(M[j],O[j]))\n",
    "        plt.ylabel(\"Frequência [Hertz]\")\n",
    "        plt.xlabel(\"Tempo [s]\")\n",
    "        plt.show()\n",
    "\n",
    "\n",
    "plt.rcParams['figure.figsize'] = [6, 6]\n",
    "N = 1000\n",
    "t = np.linspace(0, 1-1/N, N)\n",
    "f = 0.5*np.sin(2*m.pi*96*t)+0.5*np.sin(2*m.pi*235*t)\n",
    "g = np.concatenate( (np.sin(2*m.pi*96*t[0:N//2]), np.sin(2*m.pi*235*t[N//2:N])) , axis=0) \n",
    "\n",
    "M = np.array([20, 50, 100, 200])\n",
    "# saltos de 60% do valor da janela\n",
    "O = 0.6 * M\n",
    "# a figura 4 é um exemplo de que mesmo diminuindo o salto\n",
    "# não se ganha muita informação temporal sobre o instante\n",
    "# da transição\n",
    "O[3] = 10;\n",
    "show_spectra(N, t, g, M, O)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Exemplo 5.6"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.rcParams['figure.figsize'] = [12, 9]\n",
    "N = 1000\n",
    "t = np.linspace(0, 1-1/N, N)\n",
    "omega = 150 + 50 * np.cos(2*m.pi * t)\n",
    "f = np.sin(2*m.pi* 111 * t)+0.5*np.sin(2*m.pi* 123 * t)+0.5*np.sin(2*m.pi* omega * t)\n",
    "\n",
    "plt.plot(range(N//2),np.log(1+abs(np.fft.fft(f)[1:N//2+1])))\n",
    "plt.xlim([0, N/2])\n",
    "plt.title(\"Espectro do sinal com 3 senoides, 1 com frequencia variavel\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.rcParams['figure.figsize'] = [6, 6]\n",
    "M = np.array([200, 100, 50, 20])\n",
    "O = np.array([20, 60, 50, 12])\n",
    "\n",
    "show_spectra(N, t, f, M, O)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Figura 5.10 e outros exemplos de janelas"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 1000\n",
    "wr = np.zeros(N)\n",
    "wr[99:149] = np.ones(50)\n",
    "    \n",
    "plt.rcParams['figure.figsize'] = [15, 6]\n",
    "fig, axes = plt.subplots(1,2)\n",
    "axes[0].set_title('Janela retangular')\n",
    "axes[0].plot(range(-2, 52), wr[97:151])\n",
    "axes[1].set_title('Espectro da janela retangular')\n",
    "freqs = np.linspace(-N/2,N/2,N)\n",
    "Fmag = abs(np.fft.fft(wr))\n",
    "axes[1].plot(freqs, np.fft.fftshift(Fmag));"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_window_func(f, N, windowname):\n",
    "    plt.rcParams['figure.figsize'] = [15, 6]\n",
    "    fig, axes = plt.subplots(1,2)\n",
    "    axes[0].set_title(\"Janela {}\".format(windowname))\n",
    "    axes[0].plot(f)\n",
    "    axes[1].set_title(\"Espectro da janela {}\".format(windowname))\n",
    "    freqs = np.linspace(-N/2,N/2, f.size)\n",
    "    Fmag = abs(np.fft.fft(f))\n",
    "    axes[1].plot(freqs, np.fft.fftshift(Fmag))\n",
    "\n",
    "win = np.bartlett(50)\n",
    "plot_window_func(win, N, 'Bartlett')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "win = np.hamming(50)\n",
    "plot_window_func(win, N, 'Hamming')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sigma = 10\n",
    "win = sp.signal.gaussian(50, sigma)\n",
    "plot_window_func(win, N, 'Gaussiana')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.rcParams['figure.figsize'] = [15, 10]\n",
    "wsize = 50\n",
    "ns = range(wsize)\n",
    "win_bartlett = np.bartlett(wsize)\n",
    "win_hamming = np.hamming(wsize)\n",
    "sigma = 10\n",
    "win_gaussian = sp.signal.gaussian(wsize, sigma)\n",
    "\n",
    "\n",
    "fig, axes = plt.subplots(1,1)\n",
    "axes.plot(ns, win_bartlett, 'g-', label='Bartlett')\n",
    "axes.plot(ns, win_hamming, 'r-', label='Hamming')\n",
    "axes.plot(ns, win_gaussian, 'b-', label='Gaussiana')\n",
    "axes.legend(loc='best')\n",
    "plt.show()\n",
    "\n",
    "freqs = np.linspace(-N/2,N/2, wsize)\n",
    "spec_bartlett = abs(np.fft.fft(win_bartlett))\n",
    "spec_hamming = abs(np.fft.fft(win_hamming))\n",
    "spec_gaussian = abs(np.fft.fft(win_gaussian))\n",
    "\n",
    "fig, axes = plt.subplots(1,1)\n",
    "axes.plot(freqs, np.fft.fftshift(spec_bartlett), 'g-', label='Bartlett')\n",
    "axes.plot(freqs, np.fft.fftshift(spec_hamming), 'r-', label='Hamming')\n",
    "axes.plot(freqs, np.fft.fftshift(spec_gaussian), 'b-', label='Gaussiana')\n",
    "axes.legend(loc='best')\n",
    "plt.show()\n"
   ]
  }
 ],
 "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
}