scipy-fft.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Scipy FFTPack"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "#### One dimensional DFT"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The FFT y[k] of length NN of the length-NN sequence x[n] is defined as\n",
    "\n",
    "\\begin{align}\n",
    "y[k] = \\sum_{n=0}^{N-1} e^{-2 \\pi j \\frac{k n}{N} } x[n] \\, ,\n",
    "\\end{align}\n",
    "<br>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "and the inverse transform is defined as follows\n",
    "\\begin{align}\n",
    "x[n] = \\frac{1}{N} \\sum_{k=0}^{N-1} e^{2 \\pi j \\frac{k n}{N} } y[k] \\, .\n",
    "\\end{align}\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from scipy.fftpack import fft, ifft\n",
    "import numpy as np\n",
    "\n",
    "x = np.array([1.0, 2.0, 1.0, -1.0, 1.5])\n",
    "y = fft(x)\n",
    "y\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "yinv = ifft(y)\n",
    "yinv"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For N even, the elements y[1]...y[N/2−1] contain the positive-frequency terms, and the elements y[N/2]...y[N−1] contain the negative-frequency terms, in order of decreasingly negative frequency. For N odd, the elements y[1]...y[(N−1)/2] contain the positive- frequency terms, and the elements y[(N+1)/2]...y[N−1] contain the negative- frequency terms, in order of decreasingly negative frequency."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In case the sequence x is real-valued, the values of y[n] for positive frequencies is the conjugate of the values y[n] for negative frequencies (because the spectrum is symmetric). Typically, only the FFT corresponding to positive frequencies is plotted."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from scipy.fftpack import fft\n",
    "N = 600\n",
    "T = 1.0/800.0\n",
    "x = np.linspace(0.0, N*T,N)\n",
    "y = np.sin(50.0 * 2.0 * np.pi *x) + 0.5 * np.sin(80.0 * 2.0 * np.pi * x)\n",
    "yf = fft(y)\n",
    "xf = np.linspace(0.0, 1.0/(2.0*T), N//2)\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "plt.plot(xf, 2.0/N*np.abs(yf[0:N//2]))\n",
    "plt.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The FFT input signal is inherently truncated. This truncation can be modelled as multiplication of an infinite signal with a rectangular window function. In the spectral domain this multiplication becomes convolution of the signal spectrum with the window function spectrum, being of form \\sin(x)/x. This convolution is the cause of an effect called spectral <b>leakage</b> (see [WPW]). Windowing the signal with a dedicated window function helps mitigate spectral leakage. The example below uses a Blackman window from scipy.signal and shows the effect of windowing (the zero component of the FFT has been truncated for illustrative purposes)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAD8CAYAAAB+UHOxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xd4VGXa+PHvTQgdgrRQEqmCIigQRETBYMUCWBfQZXdf\nFF71dVfd4qq41vW3rm6xoa6CvSBYAVmRdRPpCigdEURKAOklgVAS7t8fzwwZwiSZybST5P5c17lm\nzskpdw5k7nnKeR5RVYwxxlQ91RIdgDHGmMSwBGCMMVWUJQBjjKmiLAEYY0wVZQnAGGOqKEsAxhhT\nRVkCMMaYKsoSgDHGVFGWAIwxpoqqnugAStOkSRNt06ZN2Mft37+funXrRj+gCFlc4fNqbBZXeCyu\n8EQa18KFC3eoatMyd1RVzy4ZGRlaHllZWeU6LtYsrvB5NTaLKzwWV3gijQtYoCF8xloVkDHGVFFx\nqwISkbrA88BhIFtV347XtY0xxpwoohKAiLwiIttEZFmx7QNEZJWIrBGRe3ybrwHeV9WRwKBIrmuM\nMSZykZYAXgOeA97wbxCRJGAMcDGQA8wXkUlAGrDUt1thhNc1xsTAkSNHyMnJ4eDBg1E9b0pKCitX\nrozqOaOhosdVq1Yt0tLSSE5OLtd1IkoAqjpDRNoU29wLWKOqawFEZDwwGJcM0oBFWPdTYzwpJyeH\n+vXr06ZNG0QkaufNzc2lfv36UTtftFTkuFSVnTt3kpOTQ9u2bct1HdEIJ4TxJYApqtrFt34dMEBV\nb/atDwfOBv6IKy0cBGaV1AYgIqOAUQCpqakZ48ePDzumvLw86tWrF/ZxsWZxhc+rsVXWuFJSUmjf\nvn1UP/wBCgsLSUpKiuo5o6Gix6Wq/PDDD+zdu/e47f3791+oqj3LOj4WjcDB/ueoqu4H/qesg1X1\nJeAlgJ49e2pmZmbYAWRnZ1Oe42LN4gqfV2OrrHGtXLmSBg0aRC8gn4r8TTsRwomrVq1adO/evVzX\niUVVTA6QHrCeBmwO5wQiMlBEXiqe1Sq6hQsbsmRJoqMwxhgnFiWA+cApItIW2AQMBW4I5wSqOhmY\n3LNnz5ExiC8h9u2D3/++G3XqwP79iY7GGO9KSkqia9eux9Y//vhj1q1bx+DBg4/VdTdp0oT+/fsz\nceJEAJYuXXrsmBEjRvCb3/wm/oFXQBElABF5F8gEmohIDvCgqo4TkduBaUAS8IqqLo840gruDV8/\nqcOHExuHMV5Xu3ZtFi1adNy2devW0bdvX6ZMmXLc9tGjRwNQr169E44xZYu0F9CwErZPBaaW97wi\nMhAY2KFDh/KewnNee829tm6d0DCMMeYYTw4GVxmrgNasca/bt0f/3Bs2wJ//DGPGQDm7Axtzgjvv\nhGh9qS4srE1SEnTrBk89Vfq++fn5dOvWDYC2bdvy0UcfATBz5sxj26+//vpj3/5N+XkyAVQ2Bw7A\n3r1Qq1Yh+/YlcegQ1KwZvfNPnw4vvwx33QWnnRa98xqTCMGqgICgVUAmMpYA4mDLFvfasd0elq9I\nYceO6rRqFb3zHzjgXvPyondOY8r6ph6O3Nx8T3a3rOo8+URuZesGunkzDGQSi1c04R/8lm3bont+\nf68i611kjAmHJxOAqk5W1VEpKSmJDiUqtmyBp7gTgPP5MurtAP4PfisBGGPCYVVAcbBls9ICVw9U\nhwOWAIwpRV6Q/8iZmZmlPuEc7BhTNk+WACqbHev3Uxs3umILtkQ9AVgbgDGmPCwBxMGBda7SP69N\nG+qxn705uVE9v7UBGGPKwxJAHBzZ5BLA/nbtACjYENbQSGWyKiBjTHlYAoiDwp9cnU+e78nmgo1b\nonr+tE1fsYzTObp9Z1TPa4yp3DyZACpTN9DDh6Fgs68E0L69e10T3QRwxk+fczoraLnqv1E9rzGm\ncvNkAqhM3UCXLIFGhb42AF8CqL5jC3v2RO8arXLd1HGt18+I3kmNMZWeJxNAZTJ/PjRjG0fr1uNw\no0YU1qhFC7ZEdV6A1gdcAuiw2RKAqfiSkpLo1q3bsWXdunVkZ2eTkpJybNtFF13EY489dmw98Jhn\nnnkmbrHefPPNrFixIqxjvDSTnD0HEGPz58NlNbcjqc1ABG3eghYbtrB4MfTrF4ULHD1K28OrOEhN\nTt67FHbtgkaNonBiYxKjIg0HPXbs2LhfM5riVgIQkXYiMk5E3o/XNb3g66+hfb1tSLNmACQ1bURq\n8m4WL47SBTZsoLbm818uoBoK69dH6cTGVA5du3Zlz549qCqNGzfmDd/kHMOHDycrK6vE4yZMmMBv\nf/tbAJ5++mna+Xrx/fDDD5x33nmAe0BtwYIFgEtCo0eP5swzz6R3795s3boVgB9//JFzzjmHrl27\ncv/99x87v6ryhz/8gS5dutC1a1fee+89AG677TamTnWj6V999dWMGDECgHHjxh13fDSEVAIQkVeA\nK4Ft/snffdsHAE/jJn4Zq6qPl3QOVV0L3FSVEsCGDbB8ObRqsQ2angyA1K5Ns/r5zJoVpYusdNU/\nczmHy/k3UW1cMFVbFMeDrl1YSKjjQUd7OOhzzz2X2bNn07p1a9q1a8fMmTP5xS9+wbx583jiiSdK\nPK5fv348+eSTx67duHFjNm3axKxZs+jbt+8J++/fv5/evXvz2GOPcffdd/Pyyy9z//33c8cdd3Dr\nrbfyi1/8gjFjxhzb/8MPP2TRokUsXryYHTt2cNZZZ9GvXz/69evHnDlzGDJkCJs2bWKLbzTJWbNm\nMXTo0JB+51CFWgJ4DRgQuEFEkoAxwGVAZ2CYiHQWka4iMqXY0iyqUVcQH3/sXhsRUC3jSwCrVsG6\ndZFfo/DHDQAswv1hUAl6TpmqzV8FtGjRomMf/uCGg/ZvD2cugL59+zJjxgxmzJjBrbfeytKlS9m0\naRONGjUqtT6+efPm5OXlkZuby8aNG7nhhhuYMWMGM2fODJoAatSowZVXXglARkYG63x/4LNnz2bY\nMDd31vDhw4/tP2vWLIYNG0ZSUhKpqamcf/75zJ8/n759+zJnzhxWrFhB586dSU1NZcuWLcydO5c+\nffqE/HuHIqQSgKrOEJE2xTb3Atb4vtkjIuOBwar6F1xpoVxEZBQwCiA1NZXs7Oywz5GXl1eu46Lt\n1VfPpE2bGuieXDbt3k1eXh479u+nprr++s88s4pBgyLrEpq6cDGnARtwJYyV8+axtWHDsM7hlfsV\njFdjq6xxpaSkkJvre1L90UejExRQWFhIUlKSW8kt+0n43GL7HDhwgIKCghO2l3aMX48ePXjmmWdI\nT0/ngQceYOLEibz11lucffbZFBYWlnrOnj178sILL9CuXTt69OjBm2++yezZs3nooYfIzc2lsLCQ\n/fv3k5ubS3Jy8rExiQ4fPkx+fj65ubmoKnl5eVSvXv3YtXJzczl06BAHDx48tu3IkSPk5+fToEED\n9uzZwyeffEKvXr3YvXs3b7zxBnXq1An6ex48eLD8/+aqGtICtAGWBaxfh6v28a8PB54r5fjGwIvA\nD8C9oVwzIyNDyyMrK6tcx0XT6tWqIqoPPaSqDRqo3nWXi2vIED3asaO2bq06cGDk19l31wOqoGk1\nt6mC6lNPhX0OL9yvkng1tsoa14oVK6ITSDH79u0Led+6deuesC0rK0uvuOKKsI4JdMopp6j/8+Tx\nxx/XtLQ0/fjjj8uM69VXX9X09HR9+eWXtaCgQE899VTt3r37sZ+ff/75On/+/BNimDhxov7yl79U\nVdWBAwfqm2++qaqqzz///LH9PvjgA73kkku0oKBAt23bpieffLJu2bJFVVVvuOEGTU9P19WrV+u8\nefM0LS1N77zzzqAxBvs3AxZoCJ+xkTQCS7B8Ukqi2amqt6hqe3WlhEptzBhX5TlqFHDwYNEUYLVr\nI/n5XHstfPYZ7Izw4d2CvXnkUZfkZicBoHusCsiY4s4++2w6duwIuCqhTZs2HWvI3bx5M5dffnnQ\n4/r27cvGjRvp168fSUlJpKenHzsuVE8//TRjxoyha9eubNq06dj2q6++mjPOOIMzzzyTCy64gCee\neILmzZsD0KdPHwoKCujQoQM9evRg165dQaudIhZKltDgJYBzgGkB6/cS4jf7EK41EHipQ4cOwdNy\nGRL97Wz7dtX69VWHDVPVwkL3zfyhh1xct96q2qSJLl7sNj/7bGTX2nHNSN1Mcz3rLNVc6uqR3/w2\n7HMk+n6VxquxVda4vFACiKfKEFeiSgDzgVNEpK2I1ACGApMiON8xWsGfBP7zn90AbfffDxw65DbW\nquVea9eG/HzOOAO6d4eXXgItsdxUNt2XSy71adYM9pJCwQ7rBWSMCU1ICUBE3gXmAp1EJEdEblLV\nAuB2YBqwEpigqstjF2rFsGSJq/4ZMQI6d8ZV/8AJCQBV7rgDli6FadPKfz3NyyOPejRrBntoSOEu\nqwIyxoQmpASgqsNUtYWqJqtqmqqO822fqqod1dXrPxatoCrqYHAHD8KvfgWNG8PjjwdshOMTwNGj\ncOQIw4ZBq1bw2GPlLwVIrisBpKf7SgCWAEyENJIiqYmrSP+tPDkWUEWsAlKFX/8avv3WVes0buz7\nQbAEAJCfT40aMHo0zJoFkyeX77rV9rsEcOaZrgRgVUAmErVq1WLnzp2WBCoAVWXnzp3U8n+2lION\nBRQFqnDffTB2rHsdNCjgh6UkAFJSuPlmePpp+O1v4cILoW7d8K4tB1wV0DkZMIcU2Lsm4t/HVF1p\naWnk5OSwPcrzlh48eDCiD6pYqehx1apVi7S0tHJfxxJAhFThwQddlc///q9rAD5OaQkASE6GF1+E\n/v3hnnvg2WfDu77k5XK4Zn1OPhlm1G1Icp6VAEz5JScn07Zt26ifNzs7m+7du0f9vJGq6nF5sgqo\norQBHDgAw4e7ByZHjIDnnwcp/nSEvxdQwHMAwLEEAJCZCXfcAc89B198EV4MyQdzqdWkPiJQvXEK\ntQ57+54ZY7zDkwmgIrQBLF4MffrAO++4b/0vvwzVgt3NMkoAfv/v/0HHji6h5OSEFoMWHqX20f3U\nb1kfgJrNG1JDDxdd0xhjSuHJBOBl+/a5qpqMDNiyBaZMcQ25QT/8IeQEUKcOTJzoJna/8sqQhkth\nw3cHqIZyUrob0Kp+K5cw96yzaiBjTNksAYRo71548klo1w7++lf45S/dSMwlPEFeJMQEAHDGGS4J\nLFsG118fdJfjrFrgskSz9q4E0KS9SwDffWXVQMaYslkCKMOiRa5xt2VLuPtu6NnTzfI1blyIE2+F\nkQAALr0U/vUv+PxzuOKK0ksCs6e5kQdbdnIJ4NQebrTAr7LLyBzGGIP1AjpBYSF88w189JFbvvvO\nfWbfcAPcequr+glLmAkA4KabXJvxr34FF1zgSgVt2hy/T34+ZE922aFWE1cFVLtBMgBfzzkSZpDG\nmKqoSicAVVePv3y5+9CfOdM9lLV3rxvJMzPTPdw1bBicdFI5L1KOBADw859DSop77d7dlQquv76o\nl9HYsUCer3hQ35UASHYJYO33BWze7EotxhhTEk8mABEZCAzs0KFDuY7Pza3ODz+4Adn8y549sHlz\n0bJunatr37276LhTT4Wf/cxN1n7ZZQFP80bCnwBK6QZakoEDXRXU0KEwZIh7YGzIEJegHnkE7ume\nB99SlACqu3/OZI4wbhz86U9RiN8YU2l5MgGo6mRgcs+ePUeW5/hnn+3A9OnBf1ajhvtmnJ7uPuy7\ndClamjSJIOiSlLME4Ne2rSuVjB0L//yne14A4JJL4N4huXATJ5QAzu1VwJgxrs3Cn3eMMaY4TyaA\nSF1++RZuvLE5deu6oRXq1XPVKS1bum/1JzysFUshPAhWluRk1/5wyy2uykoEWrQAxrpGYPzzmvpK\nAD+7+giP3+vm3v7jH6PwOxhjKqW4JgARuQq4AmgAjFPVz2NxnW7d9pKZGYszl8PBg67Y4X9QIDnZ\nvQ8jAfiJFKvXP3zYvfqTi68E0L3LEa66Ch5+GK65Bk45JYL4jTGVVsjdQEXkFRHZJiLLim0fICKr\nRGSNiNxT2jlU9WNVHQncAgwpX8gVzMGDRdU/4D7F/XMCRKqgwL36vvkfey0o4Nln3cNlV13l2gyM\nMaa4cJ4DeA0YELhBRJKAMcBlQGdgmIh0FpGuIjKl2NIs4ND7fcdVfsUTAEQ/ASQluVdfCYAjR0hL\ngwkT4PvvYcCA4xu7jTEGwkgAqjoD2FVscy9gjaquVdXDwHhgsKouVdUriy3bxPkr8G9V/SZ6v4aH\nxSMBBCkBQNEzBAsXugfYliyJ/JLGmMoj0jaAVsDGgPUc4OxS9v81cBGQIiIdVPXF4juIyChgFEBq\nairZ2dlhB5WXl1eu42LhtA0bqK/K19nZx+LqBeStX8+KCGM8efVq2gFfzp6NJidTa/NmegMrlyxh\na4sWADRsCE891YAHHzydXr2qM3LkWq66atOxQgN4634V59XYLK7wWFzhiVtcocwc71+ANsCygPXr\ngbEB68OBZ8M5Z2lLRkbGCbPdhyIrK6tcx8XEVVepdu2qqgFxnXmm6qBBkZ/7kUdUQbWw0K2vX+/W\nx407YdefflIdMMD9OCND9csvi37mqftVjFdjs7jCY3GFJ9K4gAUawmdspGMB5QDpAetpwOYIz1lh\n5gMISTyqgAJ7GAEcOXEoiNRUmDoV3nvPPQh3/vlw8cUwb17kYRhjKqZIE8B84BQRaSsiNYChwKRI\ng9IKMB9AyGKdAKoH1OIVawMoTsQ9/LZmDfztb+4p43POgbvuOpOPPnLjIBljqo5wuoG+C8wFOolI\njojcpKoFwO3ANGAlMEFVl8cm1AoqlgmgsPD4BFBKCSBQnTrwu9/Bjz/CE0/Ali21ueYa6NDBTUwT\n6oQ0xpiKLZxeQMNUtYWqJqtqmqqO822fqqodVbW9qj4WjaAqVRXQoUOxLQEEtuaWUQIorl49+MMf\n4O23v+KDD9yIo6NHw8knu6Em3nnHTXtpjKmcPDkfgFUBhah4FVCIJYDikpKUa66BrCz44Qd44AFY\nvRpuvBGaNnWjkI4fH9osZcaYisOTCaBSiWcVUJglgGDatYOHHnKJICvLzUkwa5YbErtpU/dk8Wuv\nwU8/RRK4McYLLAHEWn5+/EoASUmupTfMEkAw1aq5+RDGjHFtAjNnusHoFi6E//kfNxhdRgbcfz/M\nmWMNyMZURJ5MAJWqDWD//qLROv1i1QYALiFEUAIIJikJzjvPjS66YQN8+61rLK5TBx5/HM4915UO\nfvYzeOEFN4uaeyzEGONlnkwAlaYNQBXy8oIngMOHI//aXLwEAK4dIAolgJKIQLducO+9rlSwfbt7\ntmDQIFcSuO02OO00aNXKTaP50kuuPcESgjHeUynnA/CMAwfcJ1+wBACufaBu3fKfv3gbALj1GCaA\n4k46yX3z/9nP3K/qbzvwL+++6/Zr2dI9c9Cnj3vt0cMmqzEm0SwBxFJesQlb/AInhYkkAQSrAkpO\njnoVUKhE3LMEHTrAyJEuIXz/PWRnw4wZMHcufPCB27dGDZcEzjmnaGnVKs6T9RhTxVkCiKVQEkAk\nglUBxbkEUBoR6NTJLf/7v27bTz+5ROBfnn/eTXUJbriKHj1c43JGhntvVUfGxI4lgFhKRAJIYAkg\nFM2bw9VXuwVcU8iiRfDVV66H0cKFMG0aHD3qfp6S0ofevYsSwxlnuK6qxQs+xpjwWQKIpVgngGBt\nADFuBI62GjWgVy+3+B044OYuWLgQpkzZyZYtLXjyyaK8Vrs2nH46dO1atHTp4koQVoVkTOg8mQBE\nZCAwsEOHDokOJTLxKAHEoRtovNWpA717u+X001eRmdmCgwdh6dLjl08/hVdfLTquSZOihNC5M5x6\nqluaNbPEYEwwnkwAqjoZmNyzZ8+RiY4lIomqAqpAJYBQ1aoFZ53llkDbthUlhGXL3OvYscePYdSw\noUsEnToVJYVOnaB9e1cCMaaq8mQCqDQSUQVUCUoA4WjWDC680C1+R4/Cxo2wapV7KO2779z76dPh\n9deL9ktKckmgY0f36l86dHAD41lyMJWdJYBYSkQVUCUtAYSjWjVo3dotl1xy/M/27XPJIDA5rF7t\nnlnYv//4c6SnH58U/O/btYvv72NMrMQtAYjIacAdQBPgC1V9IV7XTph4JIDi4wxVsRJAuBo0CF6V\npApbt7oH2fzLmjXu9aOPYMeO4/evX/9c2rUrSjTFl6ZNrd3BeF9ICUBEXgGuBLapapeA7QOAp4Ek\n3NzAj5d0DlVdCdwiItWAN4CqkQBEij7w/awNwHNEXBfV5s3d2EbF7d0La9e6pLB2Lcydu42Cglas\nXetKD8WHyq5d282r0Lp10WtamnvYzb80aGBJwiRWqCWA14DncB/cAIhIEjAGuBg3N/B8EZmESwZ/\nKXb8CFXdJiKDgFuBNyOMu2LwjwNU/K/c2gAqnJQU6N7dLQDZ2avJzGwFuNLDnj2wfn3w5dtv3ZhJ\nxdWte3xCCLY0b37iP7Ex0RLSfy1VnSEibYpt7gWsUdW1ACIyHhisqn/BlRaCnWcSMElEPgXeKW/Q\nFUawgeAg9m0A/qonExcibkykk05yA+UFk58PmzaVvMyYAZs3n5i7q1Vz1UnNm7vnHPyvwd43bhz7\n39VULpF8t2gFbAxYzwHOLmlnEckErgFqAlNL2W8UMAogNTWV7OzssAPLy8sr13HRdtratdRPSuJr\nXyzH4jp6lExg3cqVrIsgzp5795K/ezfLA87Rdd8+kvfs4ZswzuuV+xWMV2OLJK6WLd1SvB3i6FHY\nsyeZHTtqBiw12LXLLRs21GDxYvf+yJETB/KtVk1p2LA3jRrlctJJR2jU6DANGx4mJeXICUuDBkeo\nX7+AanEaD7gy/jvGUrziiiQBBKu9LHHkFlXNBrLLOqmqvgS8BNCzZ0/NzMwMO7Ds7GzKc1zU/f3v\n0KzZsViOi6tmTdqkptImkjhr1qReixbH/66pqXDoUFi/v2fuVxBejS2Rcam6NomffnIN11u3+t8L\nixbtRqQFW7fCihXuZ4cPBz9PtWrQqJF7gM6/NG584rq/dNOwoXutXTv8tgv7dwxPvOKKJAHkAOkB\n62nA5sjCcSrVk8DBqoAgOpPCFBYGfxLYGoErNRH3Yex/wC1QdrZ7ctpP1XVv3bHDLTt3Fr0vvu2H\nH9yYTDt2lP5fKDn5+IRQ2mvDhq79ZP36OuTkQP367k/CxnLyhkgSwHzgFBFpC2wChgI3RCOoSvUk\ncLNmwX8WjQRQAQeDM/El4j5w69VzD7eFwj+PkT9B7N7tGrn37Cl6H/i6e7frGeVfD/7fr9dxa3Xq\nuF5Q9esHX4L9rE4d92dTp07wpWZN61UVrlC7gb4LZAJNRCQHeFBVx4nI7cA0XM+fV1R1ecwirYhy\nc0t+aihWCcBKACZCIkUfum3bhnesqhuGIzA57NsHX321gvT0zuTmuvXc3BOXnJzj18P98/D3uA6W\nHPzba9Z0S40a7nXbtvZ89tnx20p6739NTnYlmOrV3Wt53vun7050wgq1F9CwErZPpZQG3fKyKqAQ\nlTQaqJUATIKIuO6t/i6ufnXrbiMzs3NY5yooKEoGeXnuz+XAgROXkrYH/iw317WHHDrklsOH3Wt+\nfgsKCkpuJ4kXfzIoSgj9yMqCvn1je11P9jCuFFVAy5a5/n3t2wf/ee3abiSzo0cpd1eMkkYDtRKA\nqQSqVy9qgI6V7OxZZGZmour+bIonCP9r4PsjR9x3r8JC9ydY3veqJy7gXtet20haWuvY/eI+nkwA\nlcIjj7gy9C23BP/5+ee7XkK33OJmTi8PawMwJipEXBVPjRruzzbRsrN/pG3b2CeAOPUCrmIOHHAD\nyIwc6frZBfPkkzB0qJsk1z/9VbisDcAYEwFPJgARGSgiL+3duzfRoZTP/PnuwzlwjOLiRODii2HX\nLjdzenkE6wZqJQBjTIg8mQBUdbKqjkpJSUl0KOUzZ4577d279P369Dl+/3BZCcAYEwFPJoAKb/Zs\nOO20kqt//Dp2dC1c0UwAVgIwxoTIEkC0HT0Kc+cWfbsvTbVqcM45bv9wqZbcDdTfxcAYY0phCSDa\nvv/e1euHkgDA7bdihXtiJhz+huNg3UDBSgHGmDJZAog2f3VOsFlFgvEninnzwruO/wM+WAkArB3A\nGFMmSwDRNmeOq/vv2DG0/c86y1UFhdsOUFIC8K9bAjDGlMGTCaBCdwOdM8d9qw91kI969eDMM13D\ncTgKC91rsG6gYFVAxpgyeTIBVNhuoFu3wsqVoVf/+J1/vksc4YwNZCUAY0yEPJkAKqz//Me9Xnxx\neMdddJEbZCScaqCy2gCsBGCMKUNcE4CI1BWRhSISdM7gCm/6dDeFkn/m8FD16+c+yL/4IvRjrARg\njIlQSAlARF4RkW0isqzY9gEiskpE1ojIPSGc6o/AhPIE6nmqLgFceGH4o3vWrw9nn11UggiFtQEY\nYyIU6ifVa8CAwA0ikgSMAS4DOgPDRKSziHQVkSnFlmYichGwAtgaxfi9Y8UK2Lw5/Oofv4suggUL\nQn8ewEoAxpgIhZQAVHUGsKvY5l7AGlVdq6qHgfHAYFVdqqpXFlu2Af2B3rhpI0eKSOVqf5g+3b1G\nkgBUITs7tP2tDcAYE6FI5gNoBWwMWM8Bzi5pZ1UdDSAivwJ2qGrQMZBFZBQwCiA1NZXsUD8QA+Tl\n5ZXruEh0HT+e2mlpfP3jj/Djj2HHJQUFnFu7Nltff53VIcyAUWf9enoBK1atYlvAORutWcMZwObR\no/n+zjtDmn07EfcrVF6NzeIKj8UVnrjFpaohLUAbYFnA+vXA2ID14cCzoZ4vlCUjI0PLIysrq1zH\nlVt+vmrduqq33VbqbmXGdcUVqu3aqR49WvY1ly51kwhNnHj89iNHVEeMcD/76KOyzxNKXAnk1dgs\nrvBYXOGJNC5ggYbwGRtJNUwOkB6wngZsjuB8x1S4B8E+/xz274eBAyM7z6BBsHYtLFlS+n7Ll8Od\nd7r3wdoAXnjBPWA2bVpk8RhjKrVIEsB84BQRaSsiNYChwKRoBKUV7UGwiRPdsM6lTQATiquucj2I\nPvig9P0hTaBRAAAa0UlEQVSef76oy2jxBABuXrsLL4TPPrNRQY0xJQq1G+i7wFygk4jkiMhNqloA\n3A5MA1YCE1R1eexC9ahDh2DSJPfh7W+ALa9mzdwzAWUlgBkzit6XVMc/YACsW1f+2caMMZVeqL2A\nhqlqC1VNVtU0VR3n2z5VVTuqantVfSxaQVWoKqDp02HfPrj++uic79prXZfSlSuD/3zXLlgW8DhG\nsBIAFJVGPNjAZYzxBk92xaxQVUATJ0LDhpFX//hdfbV7LakUMGvW8eslJYAOHSA1FWbOjE5cxphK\nx5MJoMI4dAg++QQGD3b17tHQqpUbTXTixOA/L/6BXlIVkAj07XtiwjDGGB9LAJH46CPYuxduuCG6\n5x02zPUEWrz4xJ/NmAHt2xetl1QCAJcA1q+HjRtL3scYU2V5MgFUmDaAl16Ctm3dU7zRdMMNrkTx\n6qvHb8/Lg4ULj+9uWlYCAKsGMsYE5ckEUCHaAFavhqwsuPnm8Ad/K0ujRq5X0ZtvwoEDRdvnzXOD\nwF16adG20hLAGWdAgwaWAIwxQXkyAVQIY8e6+vf/+Z/YnP///s/1+Hn77aJt06e7awZOOF/aUA/+\nfS0BGGOCsARQHocPu+qZgQOhRYvYXKNvXzevwD/+4b71Hz0K770Hl1zivtX7lVYC8J9n+XLYuTM2\ncRpjKixLAOXx9tuwfTvcemvsriEC990H330H48bB3LmuQbd4g3NZCaBfP/ca+PCYMcZgCSB8R4/C\nE09At27lH/o5VNdeC+edB3/4A9x0kxvfZ/Dg4/cpa7TPs8+GlBSYOjV2cRpjKiRLAOGaNMl9K7/7\nbvctPZZE4K234JRTXFfOTz5xs4cFKqsEkJzsqo2mTrVxgYwxx7EEEA5V+OtfXdfPaA39UJbWrV3v\nn5wcuOCCou3+b/5lJQCAyy93s5UtWhSbGI0xFZInE4BnnwP4z3/ch/HvfhfaB2+0VK/uRhsN5C8J\nhDDhC5dd5l4//TS6cRljKjRPJgBPPgdw9Cjcc4/7Rn7zzYmOxrUHgOuRVJbUVDjrLEsAxpjjeDIB\neNLEifDNN/Doo1CzZqKjgSlT3DMIaWmh7X/FFfDVV7BjR2zjMsZUGHFLACKSKSIzReRFEcmM13Wj\n4sgRGD0aunaN/rg/5XXmmfDKK6FVAYFLAKoweXJs4zLGVBihTgjziohsE5FlxbYPEJFVIrJGRO4p\n4zQK5AG1cNNJVhz//Cf88AM8/njoH7hek5Hhhoh+661ER2KM8YhQSwCvAQMCN4hIEjAGuAzoDAwT\nkc4i0lVEphRbmgEzVfUy4I/Aw9H7FWJs/Xp4+GHX//7yyxMdTfmJwM9/7sYvyqlY+dcYExuhzgg2\nA9hVbHMvYI2qrlXVw8B4YLCqLlXVK4st21T1qO+43YAHKtFDdMcd7vWZZxIbRzTceKOrBnrnnURH\nYozxANEQHw4SkTbAFFXt4lu/Dhigqjf71ocDZ6vq7SUcfw1wKdAQeEFVs0vYbxQwCiA1NTVj/Pjx\nYfw6Tl5eHvX8vWQi0OTLL+ny0EP8MGoUG4cNi/h80YorEt1vv52kAwdYMG7csQfZvBBXSbwam8UV\nHosrPJHG1b9//4Wq2rPMHVU1pAVoAywLWL8eGBuwPhx4NtTzhbJkZGRoeWRlZZXruONs3qzaqJHq\nWWepHj4c+fk0SnFFaswYVVD99ttjmzwRVwm8GpvFFR6LKzyRxgUs0BA+YyPpBZQDpAespwGbIzjf\nMQl/EEwVRoyA/Hw3Jn9ycmLiiIUhQ1w31hdfTHQkxpgEiyQBzAdOEZG2IlIDGApMikZQmugHwf7+\nd/jsM3jySejUKTExxErjxq4x+I03bIhoY6q4ULuBvgvMBTqJSI6I3KSqBcDtwDRgJTBBVZfHLtQ4\n+eIL+OMf4brr4LbbEh1NbNxxhyvdvPxyoiMxxiRQSAPaqGrQFlBVnQpEfZxhERkIDOzQoUO0T126\ndetcFcmpp7qHrGI92meidO0KF14Izz3nxjUyxlRJnhwKIiFVQNu3u7l2Cwvho49OHHa5srnzTti0\nCSZMSHQkxpgE8WQCiLu8PDdUwoYNboydjh0THVHsXX45dOkCjzyCFBYmOhpjTAJYAtizx33zX7jQ\nzbl77rmJjig+qlVzA9t9/z2pn3+e6GiMMQngyQQQt26g27e7uvD5811VyKBBsb2e1wweDGedRZvX\nX4dDhxIdjTEmzjyZAOLSBvDNN9CzJ6xY4aZavPba2F3Lq0Tgz3+m1tat9lyAMVWQJxNATB096j7s\nzj3XPfA1a1bRjFlV0cUXs6tnT3jwQdi6NdHRGGPiqGolgOXL3by6t94K550HCxa4YZKrMhFW/+Y3\ncOCAm+jeGFNlVP4EcPQozJzpJnHv0sVNjD52LHz+OTRrlujoPCE/PR3+8Af3dPCMGYkOxxgTJ3Gc\n2Tx+UpYsgdmzYdky+PJL2LLF9esfPRruussNh2COd999bpjoESNckvTgCInGmOiqlAmg+bRpMHUq\ntGkDffq4YR2uvNI+1EpTty68/jpkZsLvf2+NwsZUAZ5MAJEOBbH25ptpMX585X+aN9r69XNDQ/zt\nbzBwoHs4zhhTaXmyDSDSbqBHTjrJPvzL69FH4Ywz4Be/cGMjGWMqLU8mAJNAtWrBBx+4MZGuuw4O\nHkx0RMaYGLEEYE7UoYPrEbRwIYwa5Z6XMMZUOnFLACJSTUQeE5FnReSX8bquKadBg+Dhh92MaI88\nkuhojDExEOqEMK+IyDYRWVZs+wARWSUia0TknjJOMxhoBRzBTSdpvO5Pf4Jf/QoeesjNj2CMqVRC\n7QX0GvAc8IZ/g4gkAWOAi3Ef6PNFZBKQBPyl2PEjgE7AXFX9l4i8D3wRWegm5kTgX/+CzZvh5pvd\nXMI33pjoqIwxUSIaYv2uiLQBpqhqF9/6OcBDqnqpb/1eAFUt/uHvP/7nwGFVnSAi76nqkBL2GwWM\nAkhNTc0YP358WL8QQF5eHvU82Oe/osZV7eBBut53Hw0XL2bF6NFsv+ACz8SWKBZXeCyu8EQaV//+\n/Reqas8yd1TVkBagDbAsYP06YGzA+nDguVKOrwOMA54F/i+Ua2ZkZGh5ZGVlleu4WKvQceXlqfbt\nq5qUpDphQsxj8qvQ9ywBLK7wVNa4gAUawmdsJA+CBZswt8TihKoeAG6K4HomkerWhU8/dTOJDRkC\nO3a4QfWMMRVWJL2AcoD0gPU0YHNk4ThxmxDGhKd+fZg2zQ2rcdttcP/91kXUmAoskgQwHzhFRNqK\nSA1gKDApGkFpIiaFN6GpUwc+/BBGjoTHHnO9hOxhMWMqpFC7gb4LzAU6iUiOiNykqgXA7cA0YCUw\nQVWXxy5U4xnVq7veQY884h4Y69cPcqxnrzEVTUhtAKo6rITtU4GpUY2IyAeDM3Eg4p4T6NoVhg93\nE+u8/z707ZvoyIwxIfLkUBBWBVSBXHUVfPUVNGgA/fu7weQKCxMdlTEmBJ5MAKaC6dzZTa85ZAg8\n8IBLBBs2JDoqY0wZLAGY6EhJgbfecm0C334LZ54JEyYkOipjTCk8mQCsG2gFJeLaAxYtgo4dXYng\n2mvdlJzGGM/xZAKwNoAKrn17mDUL/vIX9/BY585uMDl7ZsAYT/FkAjCVQHIy3HMPLFniZhi76Sa4\n6CJYuTLRkRljfCwBmNjq2BGysuCFF+Cbb1wyuOsu2LMn0ZEZU+VZAjCxV60a3HILfP89jBgBTz/t\nEsPYsdZl1JgEsgRg4qdpU/cE8YIFLgGMHAm9esF//pPoyIypkiwBmPjr0QNmzoS333ajil58sWsf\nmD8/0ZEZU6V4MgFYN9AqQARuuMFVCz31FCxe7EoD110H332X6OiMqRI8mQCsG2gVUrMm3HEHrF3r\n5h7+/HM4/XSXHJYtK/NwY0z5eTIBmCqofn148EH44Qf43e9g8mTo2pXTH3gAFi5MdHTGVEpxSwAi\n0ldEXhSRsSIyJ17XNRVM06bwxBOwbh088AANv/0WevZ0M5HNnp3o6IypVEKdD+AVEdkmIsuKbR8g\nIqtEZI2I3FPaOVR1pqreAkwBXi9/yKZKaNwYHn6YeePHuyeKFyyA886DPn3csNMFBYmO0JgKL9QS\nwGvAgMANIpIEjAEuAzoDw0Sks4h0FZEpxZZmAYfeALwbhdhNFVBYt657onjdOnjmGdi2Da6/Hjp0\ngH/8A6yjgDHlFlICUNUZwK5im3sBa1R1raoeBsYDg1V1qapeWWzZBiAiJwN7VXVfNH8JUwXUqQO/\n/jWsWgUffQStW7u2gvR0uPNOWLMm0REaU+GIhjhAl4i0Aaaoahff+nXAAFW92bc+HDhbVW8v5RwP\nA9NUtcQ2ABEZBYwCSE1NzRg/fnxov0mAvLw86tWrF/ZxsWZxha+02OqtWkX6++/TNCuLaoWF7MrI\nYPOgQezs0wetHtJkdzGJK5EsrvBU1rj69++/UFV7lrmjqoa0AG2AZQHr1wNjA9aHA8+Ger5QloyM\nDC2PrKysch0XaxZX+EKKbdMm1UcfVU1PVwXVFi1U//Qn1Q0bEhtXAlhc4amscQELNITP2Eh6AeUA\n6QHracDmCM53jD0IZsLSsiXcf797lmDSJOjeHf78Z2jTBgYPhilT4MiRREdpjOdEkgDmA6eISFsR\nqQEMBSZFIyi1B8FMeVSvDgMHujkI1q51jcdffeW2pae7NoOlSxMdpTGeEWo30HeBuUAnEckRkZtU\ntQC4HZgGrAQmqOry2IVqTBjatIHHHoONG+GTT1z30WefdcNRZ2S4HkU7diQ6SmMSKtReQMNUtYWq\nJqtqmqqO822fqqodVbW9qj4WraCsCshETXIyDBoEH34Imze7D35ww0+0bAlXX+16FR08mNg4jUkA\nTw4FYVVAJiaaNHFdSRcudDOV/eY3MHcuXHMNpKbCr34Fn31m7QWmyvBkAjAm5rp2hb/9DXJyYNo0\nN3n9xx/DZZdBixZuApvsbJuwxlRqlgBM1Va9OlxyiZu0futWlwQuvhjefBP694eTT3bVRV9+acnA\nVDqeTADWBmASomZN12303XfdkBPjx8NZZ7lZzDIzXclg5EiYOhU5fDjR0RoTMU8mAGsDMAlXty4M\nGeJKBNu3w3vvwYUXutcrruDcq6+GYcNgwgTIzU10tMaUiycTgDGeUr8+/OxnrmSwfTt8+inb+vd3\ncxkPGeKGsB440JUUNm5MdLTGhCy2A6YYU9nUrAmXX873derQsm9fN0fBhx+6ksKUKW6fM86AK65w\nS+/ekJSU2JiNKYGVAIwpr6Qk6NfPzWn844+wfDn89a/QsKGb1Oa886BZM7jxRnjnHdhVfEBdYxLL\nSgDGRIMIdO7slrvvht273fzGn34K//63SwDVqrknki+91PU06tnTSgcmoawEYEwsnHSSax944w34\n6Sf3wNl998GBA/CnP7mqoaZN4brrXNvBjz8mOmJTBXmyBCAiA4GBHTp0SHQoxkQuKcl94PfuDY8+\n6hqSv/gCpk93ywcfuP3at3clg4svhgsucFVJxsSQJ0sA1g3UVGpNm8LQoTBuHKxfDytXujGKOneG\nt95yTyU3bgxnn+2qk6ZOhX02iZ6JPk+WAIypMkTg1FPd8utfu3GIvvrKtR9kZ7sG5iefdO0HGRnu\ngbTzz3cNzPYFyUTIEoAxXpKc7D7czzvPrR84APPmuWSQnQ1PP12UEHr0cAkhM9MSgimXuCUA34Tw\nz+Aml/9eVR+P17WNqbDq1HHtARdc4Nbz849PCM884wa1E4HTT6dj27buYbRzz4W2bd12Y0oQUgIQ\nkVeAK4Ft6psU3rd9APA0kISbH7i0D/WuwPuq+paIvBdBzMZUXbVru0Hq+vd36/n5rofRrFkwezbN\nvvgCJk92P2ve3HU7Pfdct3TvDjVqJC524zmhlgBeA54D3vBvEJEkYAxwMW5+4PkiMgmXDP5S7PgR\nwDzgfREZAbwZWdjGGMAlhIASwqwvviCzaVP3hPKcOUVPKgPUquUGt+vTxzUw9+oFrVolMHiTaOIm\nkA9hR5E2wBR/CUBEzgEeUtVLfev3Aqhq8Q9///G/B75W1Rki8r6qXlfCfqOAUQCpqakZ48ePD+sX\nAsjLy6NevXphHxdrFlf4vBpbRYqrxs6dNFi2jBTfUm/1aqr5hrY+1Lgxuaeeyr5TTyX31FPJ7dSJ\ngvr14xKXF1TWuPr3779QVXuWuaOqhrQAbYBlAevX4ap9/OvDgedKOb4L8D7wIvC3UK6ZkZGh5ZGV\nlVWu42LN4gqfV2Or0HHl56vOm6f6zDOqP/+5aseOqlC0dOzotj/9tOqcOar798cnrgSorHEBCzSE\nz9hIGoGDtS6VWJxQ1WW+pGGMSaRatVwV0NlnF23bswcWLICvv3bLF1+4ZxLA9Tjq1Mm1IfTo4V67\ndYNGjRITv4maSBJADpAesJ4GbI4sHMeeBDYmzho2hIsucovfpk0uGXz7rVu+/NKNaeTXurVLBv6l\nWzdIS7OeRxVIJAlgPnCKiLQFNgFDgRuiEZSqTgYm9+zZc2Q0zmeMKYdWreDqq93it307LFpUlBS+\n/RY++cRVHoEbA6lLlxMX40mhdgN9F8gEmohIDvCgqo4TkduBabieP6+o6vKYRWqMSbymTYvGK/LL\ny4MlS1wyWLbMLe+8AwFTup7TuLGrPgpMCqed5ibbMQkTUgJQ1WElbJ8KTI1qRFgVkDEVSr16rmtp\nnz5F21Rh8+ZjCWH39Ok037EDXnzRPbvg17Kla18ovrRubUNlx4Enh4KwKiBjKjgRV4XUqhVceinf\nZWTQPDMTCgvd0NdLl8J338GqVW557z03h4JfzZrQoYNLBh07uqV9e/d0c8uWlhyixJMJwBhTSSUl\nuQ/24qV7Vdixoygh+Jfly2HSJCgoKNq3Rg1XQmjXzi1t2x7/PiXFGqJDZAnAGJN4Iq59oWnTooHw\n/I4cccNm//gjrF1b9Lp2reulFFhyAFcllZYG6enHvwa+tyQBeDQBWBuAMeaY5OTgpQa/PXuOTw45\nOW7ZuNGVILZsKeql5FevHqSm0r1WLVe91KwZpKae+JqaWqmThScTgLUBGGNC1rBh0bMIwRw54pKA\nPyn4E8TWrRz97jv4/ns3mN6OHScmCnAPwjVs6Lq4Bi6B2xo0gLp1S1/q1HHJrJp35uHyZAIwxpio\nSU6Gk092SzGLs7PJzMx0KwUFsHMnbN0K27a5161bYdcuV80UuKxf70oeu3e7BBOOatVcO0ZycvCl\nWjXOys+H99+HnmUP5xMJSwDGGANQvXpRtU+oVN2kPfv2wf79pS/5+S5ZHD7sXktaVNm/dSt169aN\n3e/qYwnAGGPKS6SoiieKVmRn0+y006J6zmC8UxlljDEmriwBGGNMFeXJBCAiA0Xkpb0BY4kYY4yJ\nLk8mAFWdrKqjUlJSEh2KMcZUWp5MAMYYY2LPEoAxxlRRlgCMMaaKsgRgjDFVlGiwsS88QkS2A+vL\ncWgTYEeUw4kGiyt8Xo3N4gqPxRWeSONqrapNy9rJ0wmgvERkgarGdhCNcrC4wufV2Cyu8Fhc4YlX\nXFYFZIwxVZQlAGOMqaIqawJ4KdEBlMDiCp9XY7O4wmNxhScucVXKNgBjjDFlq6wlAGOMMWWodAlA\nRAaIyCoRWSMi9yQ4lnUislREFonIAt+2RiIyXURW+15PikMcr4jINhFZFrAtaBziPOO7f0tEpEec\n43pIRDb57tkiEbk84Gf3+uJaJSKXxjCudBHJEpGVIrJcRO7wbU/oPSslLi/cs1oi8rWILPbF9rBv\ne1sR+cp3z94TkRq+7TV962t8P28T57heE5EfA+5ZN9/2uP3/910vSUS+FZEpvvX43i9VrTQLkAT8\nALQDagCLgc4JjGcd0KTYtieAe3zv7wH+Goc4+gE9gGVlxQFcDvwbEKA38FWc43oI+H2QfTv7/j1r\nAm19/85JMYqrBdDD974+8L3v+gm9Z6XE5YV7JkA93/tk4CvfvZgADPVtfxG41ff+NuBF3/uhwHtx\njus14Log+8ft/7/ver8F3gGm+Nbjer8qWwmgF7BGVdeq6mFgPDA4wTEVNxh43ff+deCqWF9QVWcA\nu0KMYzDwhjrzgIYi0iKOcZVkMDBeVQ+p6o/AGty/dyzi2qKq3/je5wIrgVYk+J6VEldJ4nnPVFXz\nfKvJvkWBC4D3fduL3zP/vXwfuFBEJI5xlSRu//9FJA24AhjrWxfifL8qWwJoBWwMWM+h9D+QWFPg\ncxFZKCKjfNtSVXULuD9ooFmCYispDi/cw9t9xe9XAqrIEhKXr6jdHffN0TP3rFhc4IF75qvOWARs\nA6bjShx7VLUgyPWPxeb7+V6gcTziUlX/PXvMd8/+KSI1i8cVJOZoewq4GzjqW29MnO9XZUsAwTJi\nIrs5nauqPYDLgP8TkX4JjCVUib6HLwDtgW7AFuDvvu1xj0tE6gEfAHeq6r7Sdg2yLWaxBYnLE/dM\nVQtVtRuQhitpBJvU1n/9uMVWPC4R6QLcC5wKnAU0Av4Yz7hE5Epgm6ouDNxcyrVjEldlSwA5QHrA\nehqwOUGxoKqbfa/bgI9wfxRb/UVK3+u2BIVXUhwJvYequtX3B3sUeJmiKou4xiUiybgP2bdV9UPf\n5oTfs2BxeeWe+anqHiAbV4feUESqB7n+sdh8P08h9OrASOMa4KtOU1U9BLxK/O/ZucAgEVmHq6q+\nAFciiOv9qmwJYD5wiq8lvQausWRSIgIRkboiUt//HrgEWOaL55e+3X4JfJKI+EqJYxLwC19viN7A\nXn+1RzwUq2+9GnfP/HEN9fWGaAucAnwdoxgEGAesVNV/BPwoofespLg8cs+aikhD3/vawEW4Noos\n4DrfbsXvmf9eXgf8V30tnHGI67uARC64evbAexbzf0tVvVdV01S1De5z6r+qeiPxvl/Ras32yoJr\nxf8eV/84OoFxtMP1wFgMLPfHgqu3+wJY7XttFIdY3sVVDRzBfZO4qaQ4cEXNMb77txToGee43vRd\nd4nvP32LgP1H++JaBVwWw7jOwxWvlwCLfMvlib5npcTlhXt2BvCtL4ZlwAMBfwdf4xqgJwI1fdtr\n+dbX+H7eLs5x/dd3z5YBb1HUUyhu//8DYsykqBdQXO+XPQlsjDFVVGWrAjLGGBMiSwDGGFNFWQIw\nxpgqyhKAMcZUUZYAjDGmirIEYIwxVZQlAGOMqaIsARhjTBX1/wGunXEwN1S3MwAAAABJRU5ErkJg\ngg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fc7fe55b450>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from scipy.fftpack import fft\n",
    "N = 600 # number of sample points\n",
    "T = 1.0/800.0 # sample spacing\n",
    "x = np.linspace(0.0, N*T, N)\n",
    "y = np.sin(50.0 * 2.0 * np.pi * x) + 0.5 * np.sin(80.0 * 2.0 * np.pi * x)\n",
    "yf = fft(y)\n",
    "\n",
    "from scipy.signal import blackman\n",
    "w = blackman(N)\n",
    "ywf = fft(y*w)\n",
    "xf = np.linspace(0.0, 1.0/(2.0*T), N/2)\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "plt.semilogy(xf[1:N//2], 2.0/N * np.abs(yf[1:N//2]), '-b')\n",
    "plt.semilogy(xf[1:N//2], 2.0/N * np.abs(ywf[1:N//2]), '-r')\n",
    "plt.legend(['FFT', 'FFT w. window'])\n",
    "plt.grid()\n",
    "plt.show()\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In case the sequence x is complex-valued, the spectrum is no longer symmetric. To simplify working with the FFT functions, scipy provicdes the flollowing two helper functions.\n",
    "<br>The function <b>fftfreq</b> returns the FFT sample frequency points\n",
    "<br>The function <b>fftshift</b> returns the swapping the lower and upper halves of a vector."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The example below plots the FFT of two complex exponentials; note the asymmetric spectrum."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHORJREFUeJzt3X+UXOV93/H3d3/qxwrQD1hkSSARK8SyjQ0SwhziWqpJ\nECQHNT0mhbaO7dhRTm2c9NA4kUoPSclpfWzcYzctKebEbu0e2zImjqsS2TLF2jpOCkbilyVAZpFB\nEpIQSCC0kvbXzLd/zJ3du6PZndndufPce/fzOmfO3nvnmbnf2R195tFz7zzX3B0REcmXltAFiIhI\n4yncRURySOEuIpJDCncRkRxSuIuI5JDCXUQkhxTuIiI5pHAXEckhhbuISA61hdrxokWLfPny5VN6\n7OnTp5k7d25jC2qAtNYF6a1NdU2O6pqcPNa1e/fu1939wpoN3T3IbfXq1T5VO3funPJjk5TWutzT\nW5vqmhzVNTl5rAvY5XVkrIZlRERySOEuIpJDCncRkRyqGe5m9lUzO2Zme8a538zsL8ys18yeMbOr\nGl+miIhMRj099/8BbJjg/huBldFtE/Dfpl+WiIhMR81wd/cfAycmaLIR+Hp0IPdR4AIzW9yoAkVE\nZPIaMea+BDgYWz8UbRMRkUDM67jMnpktBx5y93dVue9vgc+6+0+i9UeAP3b33VXabqI0dEN3d/fq\nrVu3Tqnovr4+urq6pvTYJKW1LkhvbXmv6+dvFJjTZiyd15hzF/L++2q0PNa1fv363e6+pmbDek6G\nB5YDe8a578vAbbH1fcDiWs+pLzE1V1pry3tdG770Y//kN3Y35Lnc8//7arQ81kUTv8S0Dfid6KyZ\n9wEn3f1IA55XJPOGC0WGC8XQZcgMVHNuGTP7FrAOWGRmh4A/BdoB3P0+YDtwE9ALnAE+llSxIlnj\nQB0jnyINVzPc3f22Gvc78KmGVSSSI+6Osl1C0DdURRKknruEonAXSZJTPtFApKkU7iIJKmpYRgJR\nuIskqDQso3iX5lO4iyTIHfXcJQiFu0iCHNcBVQlC4S6SIPXcJRSFu0iCXGfLSCAKdxGRHFK4iySo\nNIlT6CpkJlK4iyTIKR1UFWk2hbtIgkpj7qGrkJlI4S6SIJ0KKaEo3EUSVDoVUukuzadwF0mQZoWU\nUBTuIgnSl5gkFIW7SKKU7hKGwl0kQe6laX9Fmk3hLpIgRx13CUPhLpKg0jdUFe/SfAp3kQSp5y6h\nKNxFEqRvqEooCneRBOkaqhKKwl0kSfoWkwSicBdJkMbcJRSFu0iCNJ+7hKJwF0mQ5nOXUBTuIgnS\n2TISisJdJEGaz11CUbiLJEizQkooCneRBJXOhFS8S/PVFe5mtsHM9plZr5ltrnL/JWa208yeNLNn\nzOymxpcqkkHKdQmkZribWStwL3AjsAq4zcxWVTT7d8AD7n4lcCvwl40uVCSLNOYuodTTc18L9Lr7\nfncfBLYCGyvaOHBetHw+cLhxJYpkl+Zzl1Da6mizBDgYWz8EXFPR5s+AH5rZp4G5wPUNqU4k4/QN\nVQnFah3sMbNbgBvc/RPR+oeBte7+6VibO6Ln+k9mdi3wFeBd7l6seK5NwCaA7u7u1Vu3bp1S0X19\nfXR1dU3psUlKa12Q3tryXtfHfnCai+can33/nAZUlf/fV6Plsa7169fvdvc1NRuWLyYw3g24FtgR\nW98CbKlosxdYFlvfD1w00fOuXr3ap2rnzp1TfmyS0lqXe3pry3tdyzc/5Ou/0Jjncs//76vR8lgX\nsMtr5La71zXm/jiw0sxWmFkHpQOm2yraHAA+CGBm7wBmAa/V8dwiueYal5FAaoa7uw8DtwM7gOco\nnRWz18zuNrObo2b/Bvg9M3sa+Bbw0egTRmTG0z8ECaGeA6q4+3Zge8W2u2LLzwLXNbY0kWwr92/U\nz5EQ9A1VkYSUM13RLiEo3EUSUg51ddwlBIW7SELKX17SfO4SgsJdJCEjwzLKdglA4S6SkHKPXeEu\nISjcRRKiUJeQFO4iCdOpkBKCwl0kIToVUkJSuIskRGPuEpLCXSQh5VDXfO4SgsJdJCFe8VOkmRTu\nIgkZnVsmcCEyIyncRRLiVZZEmkXhLpIQfUNVQlK4iyRFp0JKQAp3kYSMngqpeJfmU7iLJERfYpKQ\nFO4iCdF87hKSwl0kIbrMnoSkcBdJiL7EJCEp3EUS4kp3CUjhLpKQkbNlAtchM5PCXSQho19iUrxL\n8yncRRKiUyElJIW7SEI0n7uEpHAXSYjmc5eQFO4iCdHJMhKSwl0kIa5BdwlI4S6SkNFsV7pL8ync\nRRKmIXcJQeEukhCNykhIdYW7mW0ws31m1mtmm8dp89tm9qyZ7TWzbza2TJHs0XzuElJbrQZm1grc\nC/wacAh43My2ufuzsTYrgS3Ade7+hpldlFTBIlmhnruEVE/PfS3Q6+773X0Q2ApsrGjze8C97v4G\ngLsfa2yZItmj+dwlpHrCfQlwMLZ+KNoW98vAL5vZ35vZo2a2oVEFimSVhmMkJKv1BjSzW4Ab3P0T\n0fqHgbXu/ulYm4eAIeC3gaXA3wHvcvc3K55rE7AJoLu7e/XWrVunVHRfXx9dXV1TemyS0loXpLe2\nPNd1pK/Ilp+cBeC/3zAHM0tFXUlQXZMznbrWr1+/293X1Gzo7hPegGuBHbH1LcCWijb3AR+NrT8C\nXD3R865evdqnaufOnVN+bJLSWpd7emvLc10vvHrKL/2Th/zSP3nIC4Xi9IvyfP++kpDHuoBdXiO3\n3b2uYZnHgZVmtsLMOoBbgW0Vbb4HrAcws0WUhmn21/HcIjnmVZZEmqNmuLv7MHA7sAN4DnjA3fea\n2d1mdnPUbAdw3MyeBXYCn3H340kVLZIF8RFP1/i7NFnNUyEB3H07sL1i212xZQfuiG4iwtjeuqJd\nmk3fUBVJyNiee7g6ZGZSuIskJD5hmOZ0l2ZTuIskRHkuISncRRIS760r6KXZFO4iCRkz5q5DqtJk\nCneRJlDPXZpN4S6SkLE9d5HmUriLJCQ+FKMvMUmzKdxFEqKeu4SkcBdJyJhvqCrdpckU7iIJcXXd\nJSCFu0hCxs4to3SX5lK4iyREc8tISAp3kcRoPncJR+EukhDN5y4hKdwllV549RR3PPAUw4Vi6FKm\nTPO5S0gKd0mlR/cf57tPvMKJ04OhS5kyjblLSAp3SaWij/2ZRfGhGJ0tI82mcJdUKk+Xm+WLXOhL\nTBKSwl1SqVDMQbhrWEYCUrhLKpXDsJjd46ljJw7TsIw0mcJdUikPwzKo5y4BKdwllUYPqGY3FXUq\npISkcJdUGu25By5kGsZeQzXDL0QySeEuqVSMUj3LoagDqhKSwl1SqdxjL2Q4FbNbueSBwl1SaWRY\nJstny4wZlglYiMxICndJJc/B2TKaz11CUrhLKpWHYzKc7ToVUoJSuEsq5eNUSM3nLuEo3CWVyqGe\n6QOqms9dAqor3M1sg5ntM7NeM9s8QbsPmZmb2ZrGlSgzUTkLsxyKuj62hFQz3M2sFbgXuBFYBdxm\nZquqtJsH/AHwWKOLlJlndOKwwIVMg2aFlJDq6bmvBXrdfb+7DwJbgY1V2v058Hmgv4H1yQw1eipk\ndlNx7P86svs6JJva6mizBDgYWz8EXBNvYGZXAsvc/SEz+6PxnsjMNgGbALq7u+np6Zl0wQB9fX1T\nfmyS0loXpLe28eo6eHAAgCeefIqzB1qbXFVjfl8/e3V4ZPmxxx7n0LzpH+LK2t8xtJlcVz3hblW2\njXRDzKwF+CLw0VpP5O73A/cDrFmzxtetW1dXkZV6enqY6mOTlNa6IL21jVfXj07ugQMvc8V73sN1\nb1+Umromo3/PUXhyNwBrrr6ayy+el4q6kqC6JqcZddXTlTgELIutLwUOx9bnAe8CeszsJeB9wDYd\nVJXpyMPFOtB87hJQPeH+OLDSzFaYWQdwK7CtfKe7n3T3Re6+3N2XA48CN7v7rkQqlhkhH9dQrb4s\n0gw1w93dh4HbgR3Ac8AD7r7XzO42s5uTLlBmptxNP5DdlyEZVc+YO+6+Hdhese2ucdqum35ZMtMV\nR6YfyG4qjj3PPbuvQ7JJ31CVVCoUx/7MojHTDyjbpckU7pJKuRiWyW7pkgMKd0mlXAzLxJez+zIk\noxTukkr5OFtGp0JKOAp3SaVizoZlMvwyJKMU7pJKo+EeuJBp0HzuEpLCXVKpfO3UbE8cFl/O7uuQ\nbFK4SyrlblgmXBkyQyncJZVycUA1vpzh1yHZpHCXVMpHz13zuUs4CndJpbyd557l/4FINincJZXK\nYZjl6QfQqZASkMJdUikX0w+MmVsmu69DsknhLqlUvlhHlkNRZ8tISAp3SaV8fIkptpzh1yHZpHCX\nVBo9FTK7qaj53CUkhbukUnk4ppDhrvuYQM/uy5CMUrhLKo2OuQcuZBo05i4hKdwllXIxLBNfzu7L\nkIxSuEsqeQ4OqMYTXWPu0mwKd0kl9dxFpkfhLqk0cipkhrvuGnOXkBTukkrlA6oZzvaxl9lT112a\nTOEuqeQ5GJYpqucuASncJZXyNiuk0l2aTeEuqZSL6Qd0towEpHCXVBqZ8jfDPfe4YpanLpZMUrhL\nKuXjSkyx5XBlyAylcJdUGh1zD1zINGg+dwlJ4S6pVB7G0HnuIlNTV7ib2QYz22dmvWa2ucr9d5jZ\ns2b2jJk9YmaXNr5UmUlycUA1vpzh1yHZVDPczawVuBe4EVgF3GZmqyqaPQmscfcrgAeBzze6UJlZ\n8jbmrr67NFs9Pfe1QK+773f3QWArsDHewN13uvuZaPVRYGljy5SZJh9zy8TH3AMWIjNSPeG+BDgY\nWz8UbRvPx4HvT6cokVxcIFtj7hJQWx1trMq2qu9VM/uXwBrgA+PcvwnYBNDd3U1PT099VVbo6+ub\n8mOTlNa6IL21jVdX/8AgAIdeOUxPz/EmV9WY39f+/YMjy3v27mXO8X3TrCp7f8fQZnJd9YT7IWBZ\nbH0pcLiykZldD9wJfMDdB6o9kbvfD9wPsGbNGl+3bt1k6wWgp6eHqT42SWmtC9Jb23h1tf7fH8LQ\nEIsXL2bduitSU9dk7Cm+AC/8HIBVq1ax7oq3paKuJKiuyWlGXfUMyzwOrDSzFWbWAdwKbIs3MLMr\ngS8DN7v7scaXKTNNeTgm09dQ9erLIs1QM9zdfRi4HdgBPAc84O57zexuM7s5anYP0AV8x8yeMrNt\n4zydSF1GZ4UMW8d0+DjLIs1Qz7AM7r4d2F6x7a7Y8vUNrktmuLydCqlvqEqz6Ruqkkrl4ZgsZ6Jm\ngpSQFO6SSnm4WIfG3CUkhbukUi4OqI5Zzu7rkGxSuEsq5WFWyHjxms9dmk3hLqnj7jmZfqD6skgz\nKNwldeJ5nuVwj9eus2Wk2RTukjrxUMzwkLvmlpGgFO6SOvFAz/TFOsZdEUmewl1SZ2zPPbupOLbn\nnt3XIdmkcJfUyc2wjOZzl4AU7pI6xZwcUEVj7hKQwl1SZ+xZJgELmSZdQ1VCUrhL6sQPomb6G6rx\nDyn13aXJFO6SOnkZltHcMhKSwl1SJ5fDMsGqkJlK4S6pk8dTITP9KSWZpHCX1MnL9ANjToUMWIfM\nTAp3SZ34QdRChlNRY+4SksJdUiePE27l5XVIdijcJXVyMyyTk2/aSjYp3CV1xhxQzfBFLhwwG10W\naSaFu6ROfMw92z13aInSXcMy0mwKd0mdYk4ORDpOV2cbrS3Gm2eGQpcjM0xb6AJEKsV7uYUMp/tb\nZ4c5b3YbczpaOXzybOhyZIZRuEvqnBoYBmBOR2umh2VePnGGSxfM5czgMEdP9ocuR2YYDctI6hw8\ncQaASxfOzfTEYQeOn2bZgjksvmA2RxTu0mQKd0mdcrivXT6fw2+eZaiQvVNm3uof4o0zQ1y6cA6L\nz5vFkZNndVBVmkrhLqlz4MQZLprXyZWXzGeo4Ox/7XTokibtwPHSB9QlUc+9f6iog6rSVAp3SZ0D\nJ85wyYI5vGPxeQA8f/StwBVN3ouv9QGlcF9ywSwAXjqevQ8pyS6Fu6SKu/PS66Vwv+zCubS3GnsP\nZy/ctz11mAvndXL5xfNYu2IhbS3Gjr2vhi5LZhCFu6TKrpff4Ohb/Vy9YgHtrS1cs2Ih333iFc4O\nFkKXVre9h0+yc98xPrR6Ke2tLSyY28H7Vy7ie0++wpnB4dDlyQxRV7ib2QYz22dmvWa2ucr9nWb2\n7ej+x8xseaMLlfwbHC7yhR37OH92Oxvf+zYA/vD6lbzeN8CW7z5D/1C6A97d+Yfe1/nE13Zx4bxO\nPv6rK0bu+/0P/BKvnurnMw+m/3VIPtQ8z93MWoF7gV8DDgGPm9k2d3821uzjwBvu/nYzuxX4HPDP\nkihY8mdwuMiPnj/Gl3/8Ik8eeJMv3PIe5nSU3ppXL1/AZ264nHt27OMfXjzOze95G2tXLOBXLj6P\nJfNn09piweruGxjm4IkzHDhxhp8dOsmPnj/Gs0feYskFs/nKR65mUVfnSNv3XbaQzRt+hc9+/3l2\nvXSC37l2Oesvv4jLL54X9DVIftXzJaa1QK+77wcws63ARiAe7huBP4uWHwT+q5mZ69yvXCkWnaI7\nBXeKRWLLTtFLc8J4tG1gqEj/cIGzgwXODhXoHypwdrDI2aECJ04P8NqpAY6dGuCJF89y9OEfMFRw\nLprXyX++9b1sfO+SMfv91Pq3c+WyC/jq3/+Cr/2/l/irn/wCgBaD+XM6WNjVwXmz2pnV3hrdWpjV\n3kpnWwstZrQYmBktZpiVHldaLq0Xis5wwSkUiwwVnULBOXR4gO8dfZKhojNcKHJmsMCp/mFO9Q/R\nNzBMX/8wp2NDRa0txruXnM9//K1380+vWsKs9tZzfn+//4Ff4t1Lz+e/PNLLPTv2cc+OfXS0tnDp\nwjksnT+b+XM7WDCng/lzO+jqbGN2eyud7S3MHnldrex/s8CiV07S2mK0tRgt0c/SegstLdDW0jJ6\nf/QazcCIlin9Pko/S8uSP1Yrf83sQ8AGd/9EtP5h4Bp3vz3WZk/U5lC0/mLU5vXxnnfNmjW+a9eu\nSRe89acH+OIP9jB79uyRbdVeQbWXVXkF+qpt6vw4qvZ76x8YoLOzs6LdxDWMW0fVfVatpK7nGhoc\npL2jI9auSh0wJqiLHoV5tK2ROttauHBeJ/NbB7nuncu5ZsUC3r9yEW2tE48Unh4YZt+rp9h39BRH\n3jzL66cHOd43wKn+YfqHCvQPFaOfBQaGixTdR16XO6PrXnpN7k5bS0spIFtLgdjW0sLQ4ABdc0v/\nM2hvaWFOZytdnW2cN6udrs425s1qY2FXJ5csmMMlC+aw4sK5dHXW/4Xvw2+e5dH9x9n36in2v3aa\noyf7OXF6kDfODHIm0PGFcvCPfCBQ2jDyITDBh8Pw8BAd7e1VPygqt5zbpMpjrFaLam3snPsHKv5N\nnltL7Q+2c/ZT5SHV9n1um1EblhTY/M+vr7nv6vXYbndfU6tdPe/Gaq++8p96PW0ws03AJoDu7m56\nenrq2P1YR48Ns2xukba2gZoFVNt4zh+hzv3W+mMBDHUWaW8/94BZPR2jqm/eOhvW89ihYae9vThh\nG4h6tZR7ui2lXm95e9QDbIFYb7jcPtYGaG81OlqgoxU6Wo2OVuhsMdpbYV6HMau19Pi+vmG6Zh+F\no0f5ydFqL7i6xcDiDqADmF+tRWt0m5q+viJdXeXfUDG6DQEVc8ScgOMn4Hjv5PexALh2Nlx7SXlL\nCzCLwYLTPwyDRWewAIMFZ6gIAwXoO3OWzlmzKBRLE6wVfPSDqrRMxXLUnfDSK4DRDoBXLkcL5e3l\nf8Cj7Xxsm5F2ztCQ09Y+8iyjJl6tPhVyjceM+7j4/VGD4VlF2tqqH8SuvwNVq5baBVduaisOTyn/\nJsXdJ7wB1wI7YutbgC0VbXYA10bLbcDrRP8rGO+2evVqn6qdO3dO+bFJSmtd7umtTXVNjuqanDzW\nBezyGrnt7nWdLfM4sNLMVphZB3ArsK2izTbgI9Hyh4AfRUWIiEgANYdl3H3YzG6n1DtvBb7q7nvN\n7G5KnyDbgK8A/9PMeoETlD4AREQkkLqOALn7dmB7xba7Ysv9wC2NLU1ERKZK31AVEckhhbuISA4p\n3EVEckjhLiKSQwp3EZEcqjn9QGI7NnsNeHmKD19E6YtSaZPWuiC9tamuyVFdk5PHui519wtrNQoW\n7tNhZru8jrkVmi2tdUF6a1Ndk6O6Jmcm16VhGRGRHFK4i4jkUFbD/f7QBYwjrXVBemtTXZOjuiZn\nxtaVyTF3ERGZWFZ77iIiMoHMhLuZ/ZGZuZktitbNzP4iuij3M2Z2VaztR8zshej2kfGfdVr1/Hm0\n36fM7Idm9raU1HWPmT0f7ftvzOyC2H1borr2mdkNse0TXgC9QXXdYmZ7zaxoZmsq7gtWV5U6m77P\niv1/1cyORVc3K29bYGYPR++bh81sfrR93Pdag2taZmY7zey56G/4h2moK9rXLDP7qZk9HdX276Pt\nK8zssai2b0fTlWNmndF6b3T/8gRrazWzJ83soSA11TPpe+gbsIzSlMMvA4uibTcB36d04Z/3AY9F\n2xcA+6Of86Pl+QnUdF5s+Q+A+1JS168DbdHy54DPRcurgKeBTmAF8CKjlyt6EbiM0nWNngZWJVDX\nO4DLgR5gTWx70Loqamz6PqvU8I+Aq4A9sW2fBzZHy5tjf9Oq77UEaloMXBUtzwN+Hv3dgtYV7cuA\nrmi5HXgs2ucDwK3R9vuAfxUtfzL2b/VW4NsJ1nYH8E3goWi9qTVlpef+ReCPGXu1qo3A173kUeAC\nM1sM3AA87O4n3P0N4GFgQ6MLcve3YqtzY7WFruuH7l6+rtijwNJYXVvdfcDdfwH0Urr4+cgF0N19\nEChfAL3RdT3n7vuq3BW0rgoh9jmGu/+Y0jUR4jYCX4uWvwb8k9j2au+1Rtd0xN2fiJZPAc8BS0LX\nFdXj7t4XrbZHNwf+MfDgOLWVa34Q+KBZ468QbmZLgd8A/ipat2bXlPpwN7ObgVfc/emKu5YAB2Pr\nh6Jt421Porb/YGYHgX8BlOe3D15XzO9S6kGlra64NNUV+ncxnm53PwKloAUuirY3vd5oyOBKSj3k\nVNQVDX88BRyj1Gl6EXgz1smJ73+ktuj+k8DCBMr6EqUOafnStQubXVP9l2tPkJn9H+DiKnfdCfxb\nSkMN5zysyjafYHtD63L3/+XudwJ3mtkW4HbgT9NQV9TmTmAY+Eb5YePsv9oHfGJ1VXtY0nVNQsP+\nRk3S1HrNrAv4a+Bfu/tbE3Qum1qXuxeA90bHl/6G0hDgePtPvDYz+03gmLvvNrN1dew3kZpSEe7u\nfn217Wb2bkrjsE9Hb6SlwBNmtpbSJ9+yWPOlwOFo+7qK7T2NrKuKbwJ/Syncg9cVHaz9TeCDHg3k\nTVAXE2xvaF3jSLyuBtUS0qtmttjdj0TDG8ei7U2r18zaKQX7N9z9u2mpK87d3zSzHkpj7heYWVvU\nE47vv1zbITNrA87n3GGw6boOuNnMbgJmAedR6sk3t6akDiYkcQNeYvSA6m8w9qDNT6PtC4BfUDpo\nOT9aXpBALStjy58GHkxJXRuAZ4ELK7a/k7EHLvdTOoDYFi2vYPQg4jsT/Bv2MPaAairqimpp+j7H\nqWM5Yw+o3sPYA5efn+i9lkA9Bnwd+FLF9qB1Rfu6ELggWp4N/B2ljs13GHvw8pPR8qcYe/DygYT/\nlusYPaDa1Jqa+qZtwC/qJUbD3YB7KY2v/awiMH6X0oG5XuBjCdXy18Ae4BngfwNLUlJXL6Xxu6ei\n232x++6M6toH3BjbfhOlMyBepDSEkkRdv0WphzIAvArsSENdVeps+j4r9v8t4AgwFP2+Pk5p/PUR\n4IXo54Ja77UG1/SrlIYJnom9r24KXVe0ryuAJ6Pa9gB3RdsvA34a/Xv4DtAZbZ8VrfdG91+W8N9z\nHaPh3tSa9A1VEZEcSv3ZMiIiMnkKdxGRHFK4i4jkkMJdRCSHFO4iIjmkcBcRySGFu4hIDincRURy\n6P8DvUqWfxsNT98AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fc7f35cb410>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from scipy.fftpack import fft, fftfreq, fftshift\n",
    "N = 400\n",
    "T = 1.0/800.0\n",
    "x = np.linspace(0.0, N*T, N)\n",
    "y = np.exp(50.0 * 1.j * 2.0 * np.pi * x) + 0.5 * np.exp(-80.0 * 1.j * 2.0 * np.pi * x)\n",
    "yf = fft(y)\n",
    "xf = fftfreq(N, T)\n",
    "xf = fftshift(xf)\n",
    "yplot = fftshift(yf)\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "plt.plot(xf, 1.0/N * np.abs(yplot))\n",
    "plt.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The function <b>rfft</b> calculates the FFT of a real sequence and outputs the FFT coefficients y[n] with separate real and imaginary parts. In case of N being even: [y[0], Re(y[1]), Im(y[1]),..., Re(y[N/2])]; in case N being odd [y[0], Re(y[1]),Im(y[1]),..., Re(y[N/2]),\n",
    "Im(y[N/2])].\n",
    "\n",
    "The corresponding function <b>irfft</b> calculates the IFFT of the FFT coefficients with this special ordering."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 5.50+0.j        ,  2.25-0.4330127j , -2.75-1.29903811j,\n",
       "        1.50+0.j        , -2.75+1.29903811j,  2.25+0.4330127j ])"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from scipy.fftpack import fft, rfft, irfft\n",
    "x = np.array([1.0, 2.0, 1.0, -1.0, 1.5, 1.0])\n",
    "fft(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "yr = rfft(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 1. ,  2. ,  1. , -1. ,  1.5])"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "irfft(yr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 4.50000000+0.j        ,  2.08155948-1.65109876j,\n",
       "       -1.83155948+1.60822041j, -1.83155948-1.60822041j,\n",
       "        2.08155948+1.65109876j])"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x = np.array([1.0, 2.0, 1.0, -1.0, 1.5])\n",
    "fft(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 4.5       ,  2.08155948, -1.65109876, -1.83155948,  1.60822041])"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rfft(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Two and n-dimensional DFT"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The functions <b>fft2, ifft2</b> provide 2-dimensional FFT and IFFT, respectively. Similar, <b>fftn, ifftn</b> provide n-dimensional FFT and IFFT respectively.\n",
    "<br>\n",
    "The example below demonstrates a 20dimensional IFFT and plots the resulting time-domain signals."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD8CAYAAABn919SAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAD0tJREFUeJzt3U+IXGW+xvHnuVE3JnAT84f8M+0iNLoQcmlmvMyAEQnk\nziazGTCLoReBMOCAgovJzN24dOXugjQYzEIcBhSShSAhJA7CIOl7FcdMiMkdjDZpTbcXNO0qgd9d\n9Ek4KdPVp+qcqnP6198PhDqn6q3z/qre8HDy5tR5HRECAKx9/9J2AQCAZhDoAJAEgQ4ASRDoAJAE\ngQ4ASRDoAJAEgQ4ASdQKdNuHbV+xfc32iaaKAgAMzsP+sMj2BklfSDokaU7SRUlHI+IfzZUHAKjq\noRrv/ZmkaxHxT0my/WdJRyStGOhbtz4WE48/PnSH1z/57L79fQeeHvpYa1nd7+HLr77S4uJ3bqqe\nuuOKZjCuef33J58uRsS21drVCfTdkr4u7c9J+nm/N0w8/rhmP7owdIe/e3TPfftv1DjWWlb3e5j6\n5cHmilH9cUUzGNe8/Oi/Xq/Srs4c+oPOBH4yf2P7uO1Z27MLi9/V6A5dwrjmxLiubXXO0Ock7S3t\n75F0o7dRRMxImpGkqX87UOtOYG/8OFfn7Wl04XtoclzRHYzr2lbnDP2ipP22n7D9iKQXJJ1ppiwA\nwKCGPkOPiDu2fy/pA0kbJJ2MiEuNVQYAGEidKRdFxPuS3m+oFgBADfxSFACSINABIAkCHQCSINAB\nIAkCHQCSINABIAkCHQCSINABIAkCHQCSqPVLUUBavj/73Vv6duHGYcijfKto/m6tjjN0AEiCQAeA\nJAh0AEiCOXTUtu/A0+t2OUCMFvPmg+EMHQCSINABIAkCHQCSINABIAkCHQCSINABIAkCHQCSINAB\nIAkCHQCSINABIIlVA932Sds3bX9eem6L7bO2rxaPm0dbJgBgNVXO0N+SdLjnuROSzkXEfknnin0A\nQItWDfSI+Kuk/+t5+oikU8X2KUm/brguAMCAhp1D3xER85JUPG5vriQAwDBG/p+ito/bnrU9u7D4\n3ai7w5gwrjkxrmvbsIH+re2dklQ83lypYUTMRMRURExt2/rYkN2haxjXnBjXtW3YQD8jabrYnpZ0\nuplyAADDqnLZ4juS/iZp0vac7WOSXpN0yPZVSYeKfQBAi1Zdgi4ijq7w0vMN1wIAqIFfigJAEgQ6\nACRBoANAEgQ6ACRBoANAEgQ6ACRBoANAEgQ6ACRBoANAEgQ6ACRBoANAEgQ6ACRBoANAEgQ6ACRB\noANAEgQ6ACRBoANAEgQ6ACRBoANAEgQ6ACRBoANAEgQ6ACRBoANAEgQ6ACSxaqDb3mv7vO3Lti/Z\nfql4fovts7avFo+bR18uAGAlVc7Q70h6JSKelPSMpBdtPyXphKRzEbFf0rliHwDQklUDPSLmI+J/\niu1bki5L2i3piKRTRbNTkn49qiIBAKsbaA7d9oSkA5I+lrQjIual5dCXtL3p4gAA1VUOdNsbJb0r\n6eWI+GGA9x23PWt7dmHxu2FqRAcxrjkxrmtbpUC3/bCWw/ztiHivePpb2zuL13dKuvmg90bETERM\nRcTUtq2PNVEzOoBxzYlxXduqXOViSW9KuhwRr5deOiNputielnS6+fIAAFU9VKHNLyT9VtLfbX9a\nPPcnSa9J+ovtY5K+kvSb0ZQIAKhi1UCPiI8keYWXn2+2HADAsPilKAAkQaADQBIEOgAkQaADQBIE\nOgAkQaADQBIEOgAkQaADQBIEOgAkUeWn/0Bf1z/5TL97dI8k6Y0f51quBmjW3b/bUvf/fnOGDgBJ\nEOgAkASBDgBJMIeO2vYdeFpvfHSh7TKAkej6vHkZZ+gAkASBDgBJEOgAkASBDgBJEOgAkASBDgBJ\nEOgAkASBDgBJEOgAkASBDgBJEOgAkASBDgBJOCLG15m9IOlHSYtj6/TBtq7zGvZFxLamDlaM63W1\n/7223X/bNWQdV3Wghrb7rzS2Yw10SbI9GxFTY+2UGsai7c/Udv9dqaFpXfhMbdfQdv9VMeUCAEkQ\n6ACQRBuBPtNCn72oYTTa/kxt9y91o4amdeEztV1D2/1XMvY5dADAaDDlAgBJEOgAkASBDgBJEOgA\nkASBDgBJEOgAkASBDgBJEOgAkASBDgBJEOgAkASBDgBJEOgAkEStQLd92PYV29dsn2iqKADA4Ia+\n26LtDZK+kHRI0pyki5KORsQ/misPAFDVQzXe+zNJ1yLin5Jk+8+SjkhaMdC3bt0aExMTkqRbt27d\n99qNGzfubS8tLd3b3rhx473tXbt23feeTZs23dsuH698rKaPVz5Wv+OVj9X08Yb57srH++abb/T9\n999bDSmPK9rz5ZdfanFxcSTjeuXKlRXbTU5O3tvuWrtym7baNfFZl5aWFqusKVon0HdL+rq0Pyfp\n5/3eMDExodnZWUnShx9+eN9rr7766r3tCxcu3Nuempp6YBtJevbZZ+9tl4/X267J45WP1e945WM1\nfbxhvrtyu+PHj6tJ5XFFe3rHu67yuD733HMrtjt//vy97a61K7dpq10Tn/XChQvXV3yxpM4c+oPO\nBH4yf2P7uO1Z27MLCws1ukOXMK45Ma5rW50z9DlJe0v7eyTd6G0UETMqlm+amppieaQkyuM6OTkZ\nvf9qwPj1TsUNY6Vx7f2XZHm/PPbl53vf00a73r+XbbQb5rvr3T948KCqqHOGflHSfttP2H5E0guS\nztQ4HgCghqHP0CPiju3fS/pA0gZJJyPiUmOVAQAGUmfKRRHxvqT3G6oFAFBDrUAHpOXLJnuvGML4\n9V6u28Tx7o5r71Vb5f3y1Rnlud7e95TblV/r167u8XqvHFnpeL3tmjzeMN/dg95XBT/9B4AkCHQA\nSIIpF9S2tLT0k0uusPaVx3U9XzJYt10Tl2VWxRk6ACRBoANAEky5oLaNGzc2fh8RDK7p++mUx3WY\nK1GGvcKkyXb9rhwZV7smruKpijN0AEiCQAeAJJhyQW27du3ih0Ud0PRtkcvjup6vMOnCDcKq4gwd\nAJIg0AEgCQIdAJJgDh21bdq0aagbCaFZvevYNnG8u+O6ni8ZrNuuicsyq86nc4YOAEkQ6ACQBFMu\nqO3WrVtD3UgIzWpiTdHe47GmaDduEDaONUUBAB1CoANAEky5oDaWoOsGlqDrX1u/47EEHQCgUwh0\nAEiCQAeAJJhDR22sKZoTa4o2065Ta4raPmn7pu3PS89tsX3W9tXicfPAPQMAGlVlyuUtSYd7njsh\n6VxE7Jd0rtgHALRo1SmXiPir7Ymep49IOlhsn5J0QdIfGqwLawhrinYDa4r2r62tdmthTdEdETEv\nScXj9iGPAwBoyMivcrF93Pas7dmFhYVRd4cxKY/r7du32y4HDWFc17Zhr3L51vbOiJi3vVPSzZUa\nRsSMpBlJmpqaiiH7Q8eUx3VycjL4pWj7mlhTdKVxXc9XmHThBmFVDXuGfkbSdLE9Len0kMcBADSk\nymWL70j6m6RJ23O2j0l6TdIh21clHSr2AQAtqnKVy9EVXnq+4VoAADXwS1HUxpqi3cCaov1ra6sd\na4oCAAZGoANAEky5oDbWFO0G1hTtX1tb7VhTFAAwMAIdAJJgygW1saZoN7CmaP/a+h2PNUUBAJ1C\noANAEky5oDaWoMuJJeiaadepJegAAGsDgQ4ASRDoAJAEc+iojTVFu4E1RfvX1la7tbCmKACgYwh0\nAEiCKRfUtmvXLn4p2gFNrClaVh7X9XzJYBduEFYVZ+gAkASBDgBJMOWC2liCrhtYgq5/bW21Ywk6\nAMDACHQASIJAB4AkmENHbawp2g2sKdq/trbadWpNUdt7bZ+3fdn2JdsvFc9vsX3W9tXicXOlHgEA\nI1FlyuWOpFci4klJz0h60fZTkk5IOhcR+yWdK/YBAC1ZdcolIuYlzRfbt2xflrRb0hFJB4tmpyRd\nkPSHkVSJTmNN0W5gTdH+tfU73rpcU9T2hKQDkj6WtKMI+7uhv33g3gEAjakc6LY3SnpX0ssR8cMA\n7ztue9b27MLCwjA1ooPK43r79u22y0FDGNe1rdJVLrYf1nKYvx0R7xVPf2t7Z0TM294p6eaD3hsR\nM5JmJGlqaioaqBkdUB5X28GaojmsNK7r+QqTLtwgrKoqV7lY0puSLkfE66WXzkiaLranJZ0euHcA\nQGOqnKH/QtJvJf3d9qfFc3+S9Jqkv9g+JukrSb8ZTYkAgCqqXOXykSSv8PLzzZYDABgWvxRFbawp\n2g2sKdq/trbasaYoAGBgBDoAJMGUC2pjTdFuYE3R/rW11Y41RQEAAyPQASAJplxQG2uKdgNrivav\nra12rCkKABgYgQ4ASTDlgtpYgq4bWIKuf21ttevUEnQAgLWBQAeAJAh0AEiCOXTUxpqi3cCaov1r\n63e8dbmmKACguwh0AEiCKRfUtrS0NNSNhNBt5XFdz5cMduEGYVVxhg4ASRDoAJAEUy6ojSXouoEl\n6PrX1lY7lqADAAyMQAeAJAh0AEiCOXTUxpqi3cCaov1ra6sda4oCAAZGoANAEo6I8XVmL0j6UdLi\n2Dp9sK3rvIZ9EbGtqYMV43pd7X+vbfffdg1Zx1UdqKHt/iuN7VgDXZJsz0ZEqxctU8NotP2Z2u6/\nKzU0rQufqe0a2u6/KqZcACAJAh0Akmgj0Gda6LMXNYxG25+p7f6lbtTQtC58prZraLv/SsY+hw4A\nGA2mXAAgibEGuu3Dtq/Yvmb7xJj6PGn7pu3PS89tsX3W9tXicfOIa9hr+7zty7Yv2X6pjTpGhXFl\nXBvsk3GtYWyBbnuDpP+S9B+SnpJ01PZTY+j6LUmHe547IelcROyXdK7YH6U7kl6JiCclPSPpxeKz\nj7uOxjGujGvD3hLjOryIGMsfSf8u6YPS/h8l/XFMfU9I+ry0f0XSzmJ7p6Qr4/oeij5PSzrUdh2M\nK+PKuOYa13FOueyW9HVpf654rg07ImJekorH7ePq2PaEpAOSPm6zjgYxrmJcR4xxrWicge4HPLeu\nLrGxvVHSu5Jejogf2q6nIYwr45rSWhzXcQb6nKS9pf09km6Msf+yb23vlKTi8eaoO7T9sJb/crwd\nEe+1VccIMK6M66gxrhWNM9AvStpv+wnbj0h6QdKZMfZfdkbSdLE9reU5spGxbUlvSrocEa+3VceI\nMK6M66gxrlWN+T8XfiXpC0n/K+k/x9TnO5LmJd3W8lnHMUmPafl/qa8Wj1tGXMMvtfzP1c8kfVr8\n+dW462BcGVfGNfe48ktRAEiCX4oCQBIEOgAkQaADQBIEOgAkQaADQBIEOgAkQaADQBIEOgAk8f8+\nQvxA4rcYVQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fc7f3141090>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from scipy.fftpack import ifftn\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.cm as cm\n",
    "N= 30\n",
    "f, ((ax1, ax2, ax3), (ax4, ax5, ax6)) = plt.subplots(2, 3, sharex='col', sharey='row')\n",
    "xf = np.zeros((N,N))\n",
    "xf[0, 5] = 1\n",
    "xf[0, N-5] = 1\n",
    "Z = ifftn(xf)\n",
    "ax1.imshow(xf, cmap=cm.Reds)\n",
    "ax4.imshow(np.real(Z), cmap=cm.gray)\n",
    "xf = np.zeros((N, N))\n",
    "xf[5, 0] = 1\n",
    "xf[N-5, 0] = 1\n",
    "Z = ifftn(xf)\n",
    "ax2.imshow(xf, cmap=cm.Reds)\n",
    "ax5.imshow(np.real(Z), cmap=cm.gray)\n",
    "xf = np.zeros((N, N))\n",
    "xf[5, 10] = 1\n",
    "xf[N-5, N-10] = 1\n",
    "Z = ifftn(xf)\n",
    "ax3.imshow(xf, cmap=cm.Reds)\n",
    "ax6.imshow(np.real(Z), cmap=cm.gray)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Discrete Cosine Transforms"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Scipy provides a DCT with the function dct and a corresponding IDCT with the function idct. There are 8 types of the DCT [WPC], [Mak]; however, only the first 3 types are implemented in scipy. “The” DCT generally refers to DCT type 2, and “the” Inverse DCT generally refers to DCT type 3. In addition, the DCT coefficients can be normalized differently (for most types, scipy provides None and ortho). Two parameters of the dct/idct function calls allow setting the DCT type and coefficient normalization.\n",
    "\n",
    "For a single dimension array x, dct(x, norm=’ortho’) is equal to MATLAB dct(x).\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Type I DCT"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Scipy uses the following definition of the unnormalized DCT-I (norm='None'):\n",
    "\\begin{align}\n",
    "y[k] = x_0 + (-1)^k x_{N-1} + 2\\sum_{n=1}^{N-2} x[n]\n",
    "\\cos\\left({\\pi nk\\over N-1}\\right),\n",
    "\\qquad 0 \\le k < N.\n",
    "\\end{align}\n",
    "Only None is supported as normalization mode for DCT-I. Note also that the DCT-I is only supported for input size > 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Type II DCT"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Scipy uses the following definition of the unnormalized DCT-II (norm='None'):\n",
    "\\begin{align}\n",
    "y[k] = 2 \\sum_{n=0}^{N-1} x[n] \\cos \\left({\\pi(2n+1)k \\over 2N} \\right)\n",
    "\\qquad 0 \\le k < N.\n",
    "\\end{align}\n",
    "<br>\n",
    "In case of the normalized DCT (norm='ortho'), the DCT coefficients y[k] are multiplied by a scaling factor f:\n",
    "\\begin{split}f = \\begin{cases} \\sqrt{1/(4N)}, & \\text{if $k = 0$} \\\\    \\sqrt{1/(2N)},\n",
    "& \\text{otherwise} \\end{cases} \\, .\\end{split}\n",
    "\n",
    "In this case, the DCT “base functions” \n",
    "\\begin{align} \n",
    "\\phi_k[n] = 2 f \\cos \\left({\\pi(2n+1)k \\over 2N} \\right) \n",
    "\\end{align} become orthonormal:\n",
    "\\begin{align}\n",
    "\\sum_{n=0}^{N-1} \\phi_k[n] \\phi_l[n] = \\delta_{lk}\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Type III DCT"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Scipy uses the following definition of the unnormalized DCT-III (norm='None'):\n",
    "\\begin{align}\n",
    "y[k] = x_0 + 2 \\sum_{n=1}^{N-1} x[n] \\cos\\left({\\pi n(2k+1) \\over 2N}\\right)\n",
    "\\qquad 0 \\le k < N,\n",
    "\\end{align}\n",
    "or, for norm='ortho':\n",
    "\\begin{align}\n",
    "y[k] = {x_0\\over\\sqrt{N}} + {2\\over\\sqrt{N}} \\sum_{n=1}^{N-1} x[n]\n",
    "\\cos\\left({\\pi n(2k+1) \\over 2N}\\right) \\qquad 0 \\le k < N.\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The (unnormalized) DCT-III is the inverse of the (unnormalized) DCT-II, up to a factor 2N. The orthonormalized DCT-III is exactly the inverse of the orthonormalized DCT- II. The function idct performs the mappings between the DCT and IDCT types.\n",
    "\n",
    "The example below shows the relation between DCT and IDCT for different types and normalizations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 1. ,  2. ,  1. , -1. ,  1.5])"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from scipy.fftpack import dct, idct\n",
    "x = np.array([1.0, 2.0, 1.0, -1.0, 1.5])\n",
    "dct(dct(x, type=2, norm='ortho'), type=3, norm='ortho')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 10.,  20.,  10., -10.,  15.])"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "idct(dct(x, type=2), type=2) # scaling factor 2 * N = 10"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 1. ,  2. ,  1. , -1. ,  1.5])"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "idct(dct(x, type=3, norm='ortho'), type=3, norm='ortho') # no scaling"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([  8.,  16.,   8.,  -8.,  12.])"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "idct(dct(x, type=1), type = 1) # scaling factor 2*(N-1) = 8"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The DCT exhibits the “energy compaction property”, meaning that for many signals only the first few DCT coefficients have significant magnitude. Zeroing out the other coefficients leads to a small reconstruction error, a fact which is exploited in lossy signal compression (e.g. JPEG compression).\n",
    "\n",
    "The example below shows a signal x and two reconstructions (x20x20 and x15x15)from the signal’s DCT coefficients. The signal x20x20 is reconstructed from the first 20 DCT coefficients, x15x15 is reconstructed from the first 15 DCT coefficients. It can be seen that the relative error of using 20 coefficients is still very small (~0.1%), but provides a five-fold compression rate."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.0010901402257005154"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from scipy.fftpack import dct, idct\n",
    "import matplotlib.pyplot as plt\n",
    "N = 100\n",
    "t = np.linspace(0, 20, N)\n",
    "x = np.exp(-t/3) * np.cos(2*t)\n",
    "y = dct(x, norm='ortho')\n",
    "window = np.zeros(N)\n",
    "window[:20] = 1\n",
    "yr = idct(y*window, norm='ortho')\n",
    "sum(abs(x-yr)**2)/sum(abs(x)**2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX8AAAD8CAYAAACfF6SlAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnXl8lNX1/983CwlhJ4HIIpMgqLgrorUqJorirrVKxYi7\nqVqr1lYEoyFg0yot1l9t1aLWLVFEv1apRVHMRNsKihYXVmULRhBkJ4QASc7vjzuTTIaZZPYJmfN+\nveY193nmznPPPDPzee5z7r3nGBFBURRFSSyS4m2AoiiKEntU/BVFURIQFX9FUZQERMVfURQlAVHx\nVxRFSUBU/BVFURIQFX9FUZQERMVfURQlAVHxVxRFSUBS4m2AP7KysiQnJyfk9+/atYsuXbpEzqAI\noXYFh9oVHGpXcHREuz777LNNItKnzYoi0i4fw4cPl3BwOp1hvT9aqF3BoXYFh9oVHB3RLuBTCUBj\n1e2jKIqSgKj4K4qiJCAq/oqiKAlIux3wVRRFiTT79u2jurqauro6AHr06MHSpUvjbNX+BGJXeno6\nAwcOJDU1NaQ2VPwVRUkYqqur6datGzk5ORhj2LlzJ926dYu3WfvRll0iwubNm6muriY3NzekNiLi\n9jHG/N0Ys9EYs8jP68YY82djzApjzJfGmBMi0a43U6eC09lyn9Np9yuKotTV1ZGZmYkxJt6mhIUx\nhszMzKY7mFCIlM//OeDcVl4/DxjqehQCT0So3RaMGAFjxjRfAJxOuz1iRDRaUxTlQORAF3434X6O\niLh9RORDY0xOK1UuAV5wzUGdb4zpaYzpJyLrI9G+m/x8eOYZuPBCyMs7lE8+gZkz7X5FURSlmVj5\n/AcA33psV7v2tRB/Y0wh9s6A7OxsKisrg24oOTmZ2trTmT27P+PGrcGYNYRwmKhRU1MT0ueKNmpX\ncKhdwdFe7OrRowc7d+5s2m5oaGix3V4I1K66urrQz2sgK8ECeQA5wCI/r/0LOM1j+31geGvHC3WF\nb0WFSFKSSE7OTsnKstvtiY64ojCaqF3BoXa1zpIlS1ps79ixw2/dhx/eXz8qKuz+aNOaXZ54fx6R\n9rfCtxo42GN7ILAu0o24ffwnnwxgmDmz5RiAoihKoERjDHHBggUcc8wx1NXVsWvXLo488kgWLfI5\nTybqxMrtMwu43RgzAzgZ2C4R9vcDLFhgffxz58LDD3fm1FPt9oIF6vdXFKUld90Fn33WmeRk/3X6\n94fRo6FfP1i/HoYNg8mT7cMXxx0Hjz7q/3gjRozg4osv5v7772f37t1cffXVHHXUUeF9kBCJiPgb\nY14G8oAsY0w1MAlIBRCRJ4HZwPnACqAWuD4S7Xozfrx93rABGhqSWLbMir4Kv6IoodCrlxX+tWth\n0CC7HS7FxcWMGDGC9PR0/vznP4d/wBCJ1GyfsW28LsAvItFWIBx9tH3+8ks45phYtaooyoHEo4/C\nzp27W11M5Xb1PPAAPPEETJoUfmdyy5Yt1NTUsG/fPurq6uIWUrpDxvY59FBISWnkq6/ibYmiKAcq\nbuGfOROmTCFiY4iFhYU8+OCDFBQUcO+990bG2BDokOKfmgqDBtWq+CuKEjLuMUR3Tz8/v3kMMVRe\neOEFUlJSuOqqq5gwYQILFiygoqIiMgYHSYeN7XPIITV89VXXeJuhKMoBinsM0ZNwxxCvueYarrnm\nGsCuSfr4449DP1iYdMieP0Bu7i6qq2Hr1nhboiiK0v7osOI/ePAuAOI0hVZRFKVd0+HFX/3+iqIo\n+9NhxT8raw89e9rpnoqiKEpLOqz4G2Pn+GvPX1EUZX86rPiDXey1aBHYWHKKoiiKmw4r/i+/fDCd\nOsGOHXZpNmhWL0VRFDcdVvwPP3wnzz5ry199pVm9FEUJgfJyyMmBpCT7XF4ekcPm5+fz3nvvAXD/\n/fdzxx13ROS4wdBhF3kdf/w2XnwRLroI/vQnO/CrWb0URQmY8nIoLITaWrtdVWW3AQoKwjr05MmT\nKS4uZuPGjSxcuJBZs2aFaWzwdDzxLy+HoiKcuVVMXu3gmpRSXqgo4IEHVPgVRfHgrrvo/Nln+I3p\nPH8+7NnTcl9tLdx4Izz1lO/3tBXT2cXIkSMRER555BEqKytJTk7mjTfe4F//+hc7duxg7NixnH32\n2dx222106tSJvLw8CsK84HjTsdw+7it1VRVT8oCqKh6vL+TX/cp54glN6qIoShB4C39b+4Pgq6++\nYv369aSlpTVFFb300kt56qmnePLJJ3n99dd5/fXXufzyy3nqqaeicmfQsXr+RUXNt2guulDLb7YU\nccHbBU0R+vQOQFEUHn2U3Tt3+g/pnJNjXT3eOByEkxh8/fr1FBQU8Oabb3LHHXcwZ84cRo8e3fT6\nb3/7W26++Wb+85//cLQrPn1yaxlnQiQiPX9jzLnGmOXGmBXGmAk+Xh9kjHEaYxYaY740xpwfiXa9\nKcmtwpSAKXG16yo/cUpVRCLyKYqSQJSWQkZGy30ZGXZ/iNTW1nLZZZcxbdo0hg0bxgMPPEBJSQlg\n86nfe++9nHfeeRx33HEMHDiQ6upqABobG0Nu0x9h9/yNMcnAX4Gzsbl6FxhjZonIEo9q9wMzReQJ\nY8wR2MxeOeG27U3JagcllfZKbUpASuz+KhyIaFYvRVGCwO1jLypqTuVVWhrWYG9GRgbz5s1r2h45\ncmTT9mOPPcbcuXPZvn07ixYt4uabb+b222/nX//6FxdddFFYH8UXkXD7nASsEJFVAK48vZcAnuIv\nQHdXuQdRSN4O2C/Gc3Qe2JeawcR9pTyxA3r0iEqriqJ0VAoKwp7ZEyh33HFH05TPnTt30qVLF551\nz1ePApFw+wwAvvXYrnbt86QEuNqV33c28MsItLs/BQUwfToMGMCkSqBnTz65eTovU8D6iKeLVxRF\nOXCJRM/f+NjnHVBhLPCciEwzxpwCvGiMOUpEWjiyjDGFQCFAdnY2laEMqgwYAGVl/Obqq9l82EA+\nGXIkAG+//Tnff78t+ONFmJqamtA+V5RRu4JD7QqO9mJXjx492LlzZ9N2Q0NDi+32QqB21dXVhX5e\nRSSsB3AKMMdjeyIw0avOYuBgj+1VQN/Wjjt8+HAJh+pLLxXJyJClX+wRECkvD+twEcPpdMbbBJ+o\nXcGhdgVHe7FryZIlLbZ37NgRJ0taJ1C7vD+PiAjwqQSg3ZFw+ywAhhpjco0xnYArAe9JqWuBswCM\nMcOAdOCHCLTtl63Dh0NtLQO/s2nS1O2jKIrSTNjiLyL1wO3AHGApdlbPYmPMFGPMxa5qvwZuNsZ8\nAbwMXOe6QkWNbccdB0lJdJk3l86d4fvvo9maoijKgUVEFnmJyGzsQK7nvmKP8hLg1Ei0FSj1XbtC\nTg5m6sPcc8oUbn/cAceFN01LURSlo9Cxwjt40HfuXDs3d88epuRBn1pXUKYIReVTFEU5kOmw4j/4\n6aehvr7lztpau2BDURQlwemw4v/7YRt8hnooyfURq0NRFCXB6LDiP3FpNlLSHOLBXS5Z7YifUYqi\nHJCUVJZE9HjtIZlLhxX/VTfdtF9QpsbO4QVlUhQlMZn8weTIHm/yZEpLSykvL2fhwoX86U9/iujx\nA6HDiv/GUaNsqIdu3SiuhDUMYsX46TrbR1GUuOOZzGXGjBkkJyezatUqbrzxRi6//HIAKisrOf30\n07nllluisjq6w4o/YIW+tJTJlXASC1h0jAq/oiiBUVJZgplsMJNtBBt3ORIuIF/JXAYPHswzzzzT\nVMcYQ9euXamrq2PgwIFht+lNxxZ/gCFDABjKN7rKV1GUgCnJK0EmCTLJrkd1l0vySsI6rmcyly5d\nujBnzhyf9U4//XTefvttHn74YSZNmhRWm75IGPE/1KxQ8VcUJa60lszFm6QkK8+9evViTwRSR3rT\nsdI4+iInB5KTOabzChZriAdFUUJg0hmR6Xm3lsxl8+bNFBUVsXDhQqZNm8YxxxzDnDlz2LZtG7ff\nfntE2vek44t/airk5HDElm+Yqz1/RVFCIFxXTyBkZmby5JNPAjaZS7du3bjsssui1l7Hd/sADBnC\n4EZ1+yiKorhJDPEfOpQBu1ewfl1UA4kqiqIcMCSG+A8ZQsbe7TRu3ERDQ7yNURRFiT8JI/4Ag2UF\nP0Q1hYyiKMqBQWKI/9ChAAxB/f6KkuhEOY9UzAj3c0RE/I0x5xpjlhtjVhhjJvipM8YYs8QYs9gY\n81Ik2g2UP76WQ6NJarHQy+mEqVNjaYWiKPEmPT2dzZs3H/AXABFh8+bNpKenh3yMsKd6GmOSgb8C\nZwPVwAJjzCxX9i53naHYxO6nishWY0zfcNsNhuGndGKtcTBEbM/f6YQxY2DmzFhaoShKvBk4cCDV\n1dX84PL/1tXVhSWg0SIQu9LT08MK+xCJef4nAStEZBWAMWYGcAmwxKPOzcBfRWQrgIhsjEC7AZOf\nD5uOG8KQ/63g3jJYtMgKf35+LK1QFCXepKamkpub27RdWVnJ8ccfH0eLfBMLuyLh9hkAfOuxXe3a\n58mhwKHGmP8aY+YbY86NQLtBkfWjoQzlGyorhVtvVeFXFCWxMeH6vowxVwCjReQm1/Y44CQR+aVH\nnbeAfcAYYCDwb+AoEdnmdaxCoBAgOzt7+IwZM0K2q6amhq5duzbvePQt8t6cxvEHr2H1jn5MmrSE\n44/f5v8AUWI/u9oJaldwqF3BoXYFRzh25efnfyYiJ7ZZUUTCegCnAHM8ticCE73qPAlc57H9PjCi\nteMOHz5cwsHpdDaVKypE/pJ2twhIcR5S29chhV3LpKIirCbCtqs9oXYFh9oVHGpXcIRjF/CpBKDd\nkXD7LACGGmNyjTGdgCuBWV513gDyAYwxWVg30KoItB0QNX8r5xZ5HIApedB5YxWP7yuk5m/lsTJB\nURSlXRG2+ItIPXA7MAdYCswUkcXGmCnGmItd1eYAm40xSwAncI+IbA637UC5aH4RyXvrWuxL3lPL\nRfOLYmWCoihKuyIiUT1FZDYw22tfsUdZgLtdj5hTklvF5Oubt02JfZ5UWUVJPAxSFEWJMx0/pDNQ\nstpBSWUVYIVfSlwvOBzxMklRFCWuJEZ4h9JSyMhouS8jw+5XFEVJQBJD/AsKYPp09nXuRnEl7DnI\nAdOn2/2KoigJSEK4fQAoKGDdvHVM/ut43nvtS87+afd4W6QoihI3EqPn7yJtsF14vHvFd3G2RFEU\nJb4klPh3OcwGQdq3RsVfUZTEJrHE/1BXyKHq6vgaoiiKEmcSSvyTDrbin7JBe/6KoiQ2CSX+pKez\nNTmT9E3a81cUJbFJLPEHNqUPpOsO7fkripLYJJz47+g6gF67tOevKEpik3Div6vXQPrs1Z6/oiiJ\nTcKJ/94+A+jTuBGp2xNvUxRFUeJGwol/Y38713/n1+vjbImiKEr8SDjxd0/33LFE/f6KoiQuCSf+\naYfYnr+GeFAUJZGJiPgbY841xiw3xqwwxkxopd7lxhgxxrSdXDhKuFf57lutPX9FURKXsMXfGJMM\n/BU4DzgCGGuMOcJHvW7AHcDH4bYZDr1yelBDF6Rae/6KoiQukej5nwSsEJFVIrIXmAFc4qPeg8BU\noM7HazEjq4/hOwaQ8r32/BVFSVwiIf4DgG89tqtd+5owxhwPHCwib0WgvbDo2hXWmYGkbW7HPf/y\ncsjJoSTfQE6O3VYURYkgxuZWD+MAxlwBjBaRm1zb44CTROSXru0koAK4TkTWGGMqgd+IyKc+jlUI\nFAJkZ2cPnzFjRsh21dTU0LVrV5+vbTj3KfLNByx5uyzk44dKa3YB9J07l8P++EeS9+xpyjfckJbG\n8t/8ho2jRsXNrnihdgWH2hUcHdGu/Pz8z0Sk7XFVEQnrAZwCzPHYnghM9NjuAWwC1rgedcA64MTW\njjt8+HAJB6fT6fe1p/pMlH0mRaShIaw2QqE1u0RExOEQAREQSmgqi8MRX7vihNoVHGpXcHREu4BP\nJQDtjoTbZwEw1BiTa4zpBFwJzPK4uGwXkSwRyRGRHGA+cLH46PnHitpeA0iReti4MV4mtGDqVFhc\nZF09eflVmBIwJfY1d7kktyp+BiqK0uEIW/xFpB64HZgDLAVmishiY8wUY8zF4R4/GuzpY+f68137\n8PtfsL2c3N8VQlUVH+RYV4+U2Nc8y+r/VxQlUkQkgbuIzAZme+0r9lM3LxJthkNjP4+MXsOHx9cY\n4MjyIqC21TqT86CksgoKC+2OgoKo26UoSscl4Vb4Ahyx6xMAJj16abvoTZfk+nb1nLHGR+XaWigq\nipVpiqJ0UBJP/MvLOXfuPQBMyQOqXL3pOF4Ailc6fLp6zliD+v8VRYkKiSf+RUWk7vNyscShNz11\nKjidtvz4wFLqSW7x+i4yuOvfmT4vCiWrHTG0VFGUjkjCib8/F0use9MjRsCYMfDMM3Dn/KvYldSN\nWjIorgQZ5ODeXtP5RcP/Y29KRss3ZmRAaWlMbVUUpeMRkQHfA4mS1Q47cApNi6gAcMS2N52fD6+8\nAqNHw7GdltJjzzaW3T2dydNuxumEFy+FHUDnFPhL7yImVVZBWhpMn66DvYqihE3C9fwpLaWxc/vo\nTffvD/X1cNqeuQAc/suzcTrtHcHMmdCzJ6w9rYBBjWu4bkCRrXzBBTG3U1GUjkfCif/U7wr48vbp\nbKUHkyqBQYNYfNd0pn4X+9709On2+WbHe6xKGoJzdQ4LFljhHz0arrgC5s2D55+H/3Y7Fxoa4L33\nYm6noigdj4QT/xEj4OxnC7gv5Y+UVMJHD/+bvOkFjBgRWzucTnj8cejTcx9Hba4k7cKzGTPG2pef\nb+sUFMCuXbBtGxQ89iN7K/D227E1VFGUDknCiX9+vu1ZVzf0A2DyLeuZObNZcGPFJ59At25w6/Hz\nMTU1DLh2FDNnwoIFzXXmz4c+fVyzUFNS4Oyz2TPrHaY+HF4wPkVRlIQTf7BCv7f3QQBce876mAs/\nwOWXw9mbyrlv3kWU5AF33kn+unLGj2+uc9JJUFMD77wDP/wA39X2JG3zena9k9QuFqcpinLgkpDi\n73TCsu22579w9vqm+fax5NuHy3mKQtLqtjM5DxtqwmuxWX4+PPYYNDbC30eV0/tfNgT1lDzaxeI0\nRVEOXBJO/N2zaQ49rS+NGH5x2XrGjCHmF4CjXiqii3c8Hx+LzW680bp+fvZlEZ3Z3WZ9RVGUQEg4\n8XfPphl2dAo/mL7kpK3fz9ceCx4bEdhiM6cTdu6EZ/Pax+I0RVE6Bgm3yMvtU//3v2Gd9KPPOuvz\nj6Xff+1auL7SweQ2Fpu571ImTIBrS9quryiKEigJ1/N3k5UF6+lHQ/X6mLf94YdwH6VIste112ux\nmfsu5bbbbP19qe1jcZqiKAc+CSv+mZlW/FkfO/GfOhUWLuzJBx/A7B4FMOxwGlNSKa7E9uC9QjeM\nH2/vSPr0gc+HFfCnI6bDoEF2cVq3bhrqQVGUkImI+BtjzjXGLDfGrDDGTPDx+t3GmCXGmC+NMe8b\nY+Luq3D3/FM2b7ArZ2PAiBEwefIRvPMOnH467Pl+K68mjyWvWGDNmlaFfORIKF1dQMOqKkr2nALH\nHafCryhKyIQt/saYZOCvwHnAEcBYY8wRXtUWYhO2HwO8BkwNt91wcYu/aWy0k+hjQH4+3H3311RX\nQ/KOraRv+o4Trz0yoPGGM86AHTvgiy+AE0+E//0vZhetgCm3eYhL8o2uQ1CUdk4kev4nAStEZJWI\n7AVmAJd4VhARp4i45zXOBwZGoN2wcIs/EFPXT9eu9QBs+nAxAIdcclRA7zv9dPv84YdY8d+1C5Yv\nj4aJoVFebtcdVFXZdQu6DkFR2jWRmO0zAPjWY7saOLmV+jcCPgPUGGMKgUKA7OxsKisrQzaqpqam\n1ffv3ZvUJP5fzpnDlu3bQ24rGCorBwFwzfEfwkJ4er5hSEZlQO/t3/9kXn+9hh9f38BJwNIXX2TD\n6NERsaut8+WPvnPnMvjpp/n9sA1M9k5DXFvLpKevZuKvf82qm25i46hRMbMr2qhdwaF2BUdM7BKR\nsB7AFcDTHtvjgMf81L0a2/NPa+u4w4cPl3BwOp1t1hnWebUIiDz9dFhtBUpFhUinTvWSkSHS+Ivb\nZV9GN8nKbJSKisDef911IpmZIg1760UyMkR++cuI2RbI+dqPsjJrBwgl/h8Ctl5ZWWzsigFqV3Co\nXcERjl3ApxKAdkfC7VMNHOyxPRBY513JGDMKKAIuFpE9EWg3bBr62Pg+fP99TNpbsAAGDarl6KPB\nLF5EyrFHMfNVE9ACs6lToW9f2LwZln6dDCecwPaKz5gaz9GToiK7ytiFr5STTehqZEVpV0RC/BcA\nQ40xucaYTsCVwCzPCsaY44G/YYV/YwTajAjd+qSzM7VXzHz+48fDli2dGHa4wFdfwZF2sNczmJs/\nRoyAp5+25Q8/hG+zTyR18UJOOqE+uka3gr+UmP7KuhpZUdoPYYu/iNQDtwNzgKXATBFZbIyZYoy5\n2FXtD0BX4FVjzOfGmFl+DhdTMjPhh5TYzfXftg22bElj+MEbbRf+qMAGe8HOFHr1VUhKgr/8BUrn\nnEgGu8nLXhpFi1unZLXDZ29/UqXvuwBNPK8o7YeIzPMXkdkicqiIHCIipa59xSIyy1UeJSLZInKc\n63Fx60eMDU0zfmIk/ktdOj083c70CUb8Ac48Ew4/HJYsgUtOtjaX3HFM/KZVPvggGNNyX0YGJcNu\ntauPPUlP19XIitKOSNgVvmDFv7o+duK/ZIl9HrpnkS0EKf5Op51BOZZyznh/EkB8p1WmpIAI9Olj\nVx27Vyk//rh9djjsfoD6ekqevlrn/ytKOyHhxb9qXz9k/XorYlFm6VJITW2k9/pF1ufUt2/A73UH\neXvoIfgdRWQEEA46apSXW2GffhWNySnwyCNknCc4n12Ds3+BHYQuKMD57BqGZ5dBairU1+v8f0Vp\nRyS0+GdmwvcchNmzxzrko8zSpXBbz+dIeuF5So7eDLm5AYugO8jbTTfFObyzezHX2rVMyYOkhnoa\nbvo5F2wv5yc/gUsvtYPT7ovVKGcR7NvX8hg680dR4k7ChXT2JCsL/u25yrdXr6i2N/STch7acjs0\n7mVyHpRUunrB0GacHs8ZQbd8FMfwzl7TOwGS99SS9UgRtQ0FNDbaj1JTA2++CQ9PqWLybc113Res\nSZVVlKAoSrxI6J5/LEM81NbCXZuKSG8MPxvXnDNK2UV8wjv7m975+I+qSEmB/v3tqdy5E956Cyat\n8j0jKO4zfzQOkZLgJLT4N4V1hqiL//LlkXPXpF5bwM1MZ+9BrvDO3bvHLLzzrz73Lea3znOQmmqD\nz912GyQnwyOPwN11pdTG6ULlF19xiMaN0wuBklAktPhnZcGp/BuAkmfGRfWPv3Spzd4ViV7wySfD\nyxTw3OQqSjYdZaO+xSi8c/VtpdST3GJfLRl89tNSjLHj5pdfDu++a8d5H91YwO8He838ueGG+Iaj\n9uG6QkQHpJWEIqHFv8975TzGnUD0p0wuXQpFlNKY3FI4Q+kFH3KIvWv5+GPg2GNdcZ5jw5EP/BST\nnMQOulJcCWuNg9X3TWfxsQX84x/wxht2cNoY+9F69IDfripg/Jg1lMyth4ED2fzJiviEpXC5evLy\n/a9MBqC2VqelKh2ehB7w7VRSRCd/UyYj3DNduhS+OrSAGvMg3VeuZFJlvR2gLS0Nui1jbO9//nzg\numOtQG3ZAr17R9Rmn1RUkNywjzHMYk7lucgDMGUKHOlVbcwY+Mc/bMqB886DP/wBBg1K5meDhpP1\n0ZvsyjDweGifPyTcrp7aWj7Iab77aiH6XtvBDMgryoFGQvf8/Q1eRnLK5NSpdtrjkiUwbBgk1dez\n6cTRZJzXdvau1vjRj+wFZdeQY+2OGPX+d7/6T2rowkepeTzwADzxhP18nrinpebnw6hR8NprNizF\nJ3eV0/WjORhgSh6xdbH4cvW40IB0SiKS2OLvJzZNJGeijBhhe8HLl8Oww4VO367j/z4fwogRoR9z\n6lRIS7P+9c/qrfh/89oXUXOluC9giFD7yj+Zw2hu/VU6XbtakR8zpuUFwJ172M0ll8D118OUhiI6\nU9fy4DESV38X+jPW0CJERbsMSOc5MykrC7KydHBaCZuEFn9KS6lLiu5MlPx8+OMfobER1n2+kU57\nd5N/85CAUjf6Y8QImoT+319ns6dXNp8+80VYF5S22ptxcTm7e/fnsZO/46zkD9j2eDkjRtjPN3Mm\nrYaldjrtnH9/s50mRVNcXcIJvnv4b77hYPHEF8HhoLjSdx2B2Autp+CPG9c8M2nzZti8WQenOyqu\n7905OT/6v7lAgv7H4xGLZC4iIn/5cZlsS+opk/IQOfjgkBKOtMVbb9l8Jj/mP7Ywe3bYx6yoEElO\nFhk6VKQi9RzZMeT4sI7X6vkqK5P6tOakLQJ2O4BzVVEhkpVln3f1cdjP73EcAant6wjNrrbwSjaz\nX9sZGbLovjLJyhKZNk3kqfwy2YX/+vVpGTLrZ2Xh2xWu3V7lSXmIOBwiZWXxSU5SVibicFg7MjNF\nMjNb2CQS56Qp/uzLzJQ93bvvZ2tc7TNmv99osHYRYDKXuIu8v0esxP+OO0Ru6FxuT8XSpWG16Y/b\nb7eHLz/nORGQeS98HZHjHnusPe6Hp9wj0qmTyN69IR+r1fPlcPgUHXE42jzuww9Lc6aysjLZmxrc\nRSQs0fBjd3EesjbJIa9fXibXXCOSl2df6tpV5CrKZA0OKc7zLbTuC1VUxMxTBDzaDuThForFRUWR\nt6stW70Fy7NsjEzKQ3ZnZ8dWXAO1L0yhjYidbVzoA/mfeRKo+Ce22wfrQl25u7/dWLdfArKwcTrh\nmWesj37siBU0JiVz2a8c+w2ShnLclSttueyLY2Dv3qgldA9nYLyF/7+ggNRnp7Mu1bpYAJLvuTtq\nM2laSzbz89FruOy1Al54ASoroU8fG5Ji/uACTum3Zr/67vLDR9jPvHBhz8iOsXgtPNtv+qmLtgan\nX/lvKeTksLioPPJjQP5cUdJKUESx6yfSN2yI3UI670V8rdnnSSyn+LrP5dNXt5iIENMxp0CuEPF4\nxKrn//had8O7AAAgAElEQVTjIkNZbq+wL74YVpu+ePhhkR//WOSoo0Tkyiultn9/qaiw+0PF7Up5\n9llr9rTrvhQBWXxf6L2W1s5XbV9Hq73gYKiocN15p9dIDRny7YW3hGxXm/jp+VcZh4DIwIF211VX\n2fP5wAMi3buL9OghsqW77/du6e6QadNEevTY03RHE+732ZqtrZVbuwvYRYY8lV8mFRUt7QvZ1lZ6\nqAHfmUS7lx3GnVOkck4HbGegua/be8/fGHOuMWa5MWaFMWaCj9fTjDGvuF7/2BiTE4l2I0GL+D5R\n6PmPH29DHgweDKxcye7+/QNO3egP91TKa66xU/u7r1yIAK98FJ1ey3OHlrLPa0lIQ1oGc/ODGxh3\nR/p89VU46+Iu/JOL6PzWq1S+t6/p9Uj2Vp1H3IZ3n28XGUyQUkaPhro624F9+WWYONHeobhXKa+7\nvZR9qS0nAwjw6AlVXPbrHB4+5nHy85s/U8iD7YEuPDOG4krYTCa1GZl+B6fdZFDLlV8WtYi0+vOf\nN5fdM7g8z3mr599rqmwgdyae+FpIt+GmoiYbzj+/pT1Tp9rwIO6yr3qe5X9eGd6dk6+7qElPX83u\nbHsXFa59558Pi4vK2Z2dwySv3n5rdjSkZfDPH0UpFEogV4jWHkAysBIYDHQCvgCO8KpzG/Ckq3wl\n8Epbx41Vz//99+3FdV9GN5E77wyrTV80NtqL/K9+JSK9ekn1xRdH9Ph/OK5Mao1XjyyEXktr52v0\naJHvzUHSmN45rMExT///+++LPGJ+JW4ffG1fhxR2LWseHwjALr+4eoDFeUgDSF2X3lKch6zGIWMp\nk1GjrBt42jRr07RptudfWChNPeXCQpHCrmWyvZc9TqNXz7WGDLmrb5l07y4tbA6qZx2Av9c9PjHr\nZ2Xy85+LnHeefema5DKpofXe48jr7Ge+ijIZPNgOC2Vk2DvG996zdzjdu9vP6uscnHeeyKL7yqS2\nr0NGXhdErz7Au4OR19m7sOvTymTiRJFu3ezjpZdESkrsd/SHP1hb3Lbau67m8o2dy6TK+B+jidRd\n1FtXlTW1W1Fh23b/hlqzz11+66qWkwlaa68B0+p/oi2I1YAvcAowx2N7IjDRq84c4BRXOQXYBJjW\njhsr8f/iC3sWtvc/TOTyy8Nq0xfff2+PP/2hzSIg39x6a0SPv62nw+cPOthbRe/z5RbqVatEBlDt\nsn1a+O4NN2Vlsq9T24O/QYu/D0GtIUOuTirzKfj+XDctBqr9uGRW4xBjRCZMEKmvt8IZ1MWgDVfP\nntQM+fOPyuSYY5q/VhDp2dM+T3QEJny1JkPGUiZJSc3HSE0VSUqyF4RTTrFC9vvftxQyb8FqrQ23\nYP1ApvxApuuCGdhAaw3WPs/P6PlISbH2JSVZV11KirX//sFtXwBbs8+z3NYFojgPWYNDrk4qk6OO\nEklLs67cTp1EhgyxNqWktCwffri1MzVV5LvU1r9rzzauTS2T/PzvW/w+gyFQ8W95Lx8aA4BvPbar\ngZP91RGRemPMdiDTdRFowhhTCBQCZGdnU1lZGbJRNTU1Ab1/06ZOwI/5PrkXsnQpC8No0xeLF3cH\nTqDT2rcB2JqZGdbn8sZ5XBVT8pq33be6xZVV5AfRjvf5Sk3tyU9+cgQnnriFM00FCNzw4mn8ZMrn\nVFaGn/jmR7/+Nel7988LUPfrXzN/wAC/drXF0beMJ9NrJW8XanmwsYgNZ5/FCScsa9p/3309mTGj\nG8Z8izFw0kl28BdsGey2M7eKKdc3H8/zHJsPG3nooSSefHIvtbXJpKYKn3++iBkzupGcLLz00iBO\nP30TqakbAZgx42CuvNL+Xbb5Oe7INVCdfDDj9/2el+cXYIwwYMBuvvsug7y8DSxc2Iuzz97CQ3Ov\nYtstJ3PhjpnsopQu3qFKXHSWWv7c5V5mp1zBzTev4qOPsvjooyz69Kljy5ZOzJtnvb8TJwoZGQ1c\nWjuD31HE39et5YJWznUjhuJKoYpB3MfveHPBz5hCMmCgEr6hnN9RRHFlFY0YkvZzwjV/P0Pzrmb9\nvAk8OaiIyd/cwkknbeaEE7Yxb15vvviiF0OG7GTv3iTWru3CzV2e5/7aYp4ZtJYuq5qP4xmyw12e\nlAdrcHAfpbxcWcAUd+VKmspjK8vZRWHT+fN0F3mWJ1dW8WRjITcvgrXdx7Br1x66dEllxYo0MjP3\nALBiRRq9e+9BxLBsWSduSH+eSXuKGXvtWj7M8X1csO7Ibyqnk0MB7AOcMG7cGoxZQ4QlqZlArhCt\nPYArgKc9tscBj3nVWQwM9NheCWS2dtxY9fzr6uyF98tjrxbJyQmrTV+Uldnjf/cHO53047//PaLH\nbzjYEZWev4jIu+/aHter3W+QLaaXVMxtiIjNIiKT8nz30ibl0aZdrVHs57jFeYTck/J3d+W+NR+X\n3Nxr7d3b3nhcfLF/t0BbrorVOCQpSeSII+yuq6+2dynjxvl3V7Xlnhl5nbX1hnTrvvAc3P7FL2z5\n2GNFxlImu0xg7iR3bz0ry/51QOTUU0XOP9+WL7zQDu4/8IDIdZ2elyrTtvvI3UNfm+SQt64q228g\nfualgbtPdifbO4oLLmh2l51/vsgFF9jypZeKjBy5QUCk9IjA3Uet2ffAAyI3pDcfq6GNu5/iPOv6\nuiG9TO6+u9n19dOfro16zz8S4n/Aun3ct/Zdu4rMHT5epFMnqXi/MXKuDRGZMsWe5b3FtvDBO+9E\n7uAiImVlUpsUHZ//m2/aw60iRxYfdmmEDLa0NYPI7S4J5nsUkVbXJHguOAuKNnzzu8iQ3x5h/f8n\nntjcZEaGdQ+ceaZIerrd/tOIlq4Kfy6atgTfn7tqU9f+rYpXHamyIy2zhehUVFgXz+pW1jd4ljd3\nc7QQO19lT7srKkS6dNkn3buLbO/l+/vxVd5NquzKCM1FU2UccmPnMr/2uctduuxrYeui+4K4uPix\nry3B9zeW4DlG8MgjC0P+vQYq/pGY7bMAGGqMyTXGdMIO6M7yqjMLuNZVvhyocBkZV9xxd7p0gXWm\nP+zdy8+v2BLRMAmrVsGAAZC6ZgUMHEhjWlrkDg5QUMDr506nCtfc+c6dI5bYZdo0yGENuazhxeoz\nw16b4Mnc/FIa0lrOpqklncf6lQY9g8b9PTqdsOCkX+zvXHCF7AgkFIVPCgrsOfXMSeB5eGop2lXE\nG2/A119Dt25wwQVQX2//7RUVdmZRbS1cuqDIp3umuBI2dXFwV5fppIwroLzczkA66igbHuT3v7ef\n8+67W34G75ljG+683n5eP6Sxj257NjMlDwZJFU/XjaNyiuG8l8aRg3UhtjbrqJYM5l3UnLuhZ098\nljt3brZ7xgx48MFFvPEGPHlwy++9tbbS2UdGrbU1i81ksblN+/alZnCwKeOY7mt4NbXAr33uMkgL\nW/+8qYDV901nd1//4T7ass/TveXL1pFrYHdfGwr9r1ubQ6HPnUtTedmybqH/XgMlkCtEWw/gfOBr\nrDunyLVvCnCxq5wOvAqsAD4BBrd1zFi5fSoq7ODM7dkzRUA+fvrLsNr1ZuRIkdNPFzvZPy8vKitD\nX3/ddio2nDtOpH//kI7hbdebb1oXwPZUuxS+LrNfSDMPWsVjXnYjyBtcIsa0HDQN5nvs0kXkOa6R\n3aTKt2Zg0+35ojDWP3jjz111xnX2ruXGzmVSWGjrumeEXHmlyPVpZa26PYrzCKqH7w+n09l0Xs/w\n05b7MwTaA/d0cS26r0zOO6/ljCBfZe91BS2+xyDm4wdqq9s+KStruhNszT53+ZFHFvpfAxHAfHxf\n9rV2t+C+Cw3oewwRNLyDM+C6Q4Z4xN2JsFtm4ECRJ04vE0lKisoy94cfFnntNWv63Av/JALyn9fW\nB+268j5fj5y4v3vCM7ZNxDnvPNnadYAkUS/G2HhIbrt8CZ+nq6euzsbl+Zb+UpyHbKdbk+CH7Orx\nRxszdNwzltztBhIzSEA2dXWELPietPgevWz1d+FqU7AisODJ5/+xFXdaoLaGa1+bOuHnQhWUfQRv\nq4p/GATTY0xLEzkyY5UIyJJ7Ijcgu3u3jRWzN6VzSD+AQHCLTP/+IhNO/VAEZGz3t4IWO+/ztT7N\n4VOkgh1IDpS1l90h7h6ce2769deL3HrrN02C6O4Ru4U/K0vk/vtFftHL94XKfY4jsgLXTSCxWJKT\npTgP2dPNBhFry0/teVEN19b9etg+bA20hxrJgGd+/49txOBp7TxHwr6ARbaNc9lq2RXfKBhbVfzD\nIJCT5xaQsWNFuiTvFgEpzfhtxHqJy5aJrMYRdRF1X8Ay03ZKA0ZWXTc56GM4nc4mUV2xwv+sGe/Z\nOJHAe5BN8Jz73SjFxbae58KaTZvs4jMQWePnHIcSfiIgyspkd3Z2UG4Vf66iSEeT3O937x3RslOn\nwAQrwuENAhKzYGyNkH1BiWww9oUg+CHb5YWKfxCzRP5kvSXS0Ku3fHvJbRHrJc6eHTsRvewy+xm+\n63G4yCWXBP1+t3slK0vk2mtjc9Fy428q5WocAo0CIiecYP/vZ59tZ86kpNhqAwa0Pr0zWjidTv8z\niwIR1yidy6DcGN497TAFKyy72rLVT6jouNgVA/tU/MMgmJP30kv2TOweelRIwumPv/wldj3/zEx7\n2BnJV8nuPgODPob7fL3/vl1JeUdWmewjuaXdUQp25c9/WpyHdOu2Vw45pIWuSrI1Sy680F6sdnQ5\nKGbi6qZpYNXDFRArP3WbdgWKp3hFOZ59XOP5t0JHtCtQ8U/4kM4A2dn2ubZH/4gGd1u1CkpSS5GU\nlJYvRDBbmGewtNNOg5U9h5P+QzX/eX1jSMdLTbVZxx7bNJb61M7QpYud3uhwRGwKqTf+0mlOXu3g\n6qurWLXKNpuVZafk9eplA7J1f6uclfUOpo34HvE+aIQzsvnEawpoyb+T2wwcFu1zGRQFBbBmDSVO\nCSuftHJgouJPs/hv6xJ58f/00ALMYYdBp05MqoS67OyI/vE9k6VfcAHM2TwcgPVvfRbS8R580M5/\n/tP1X5G+r4Ylv3w8+uJQWrr/3PTkZBYXlPLSS4P44x/hmGPsvPd77rHPV+wt57nUQrpvW2ujOAKC\njXwZU3H1FNDnn291jj0ZGZTcVKZCq7QLVPxpFv9NnfrD999DQ0NEjrtqlSuUc20tXHYZJU5h/owZ\nEf3jeyZLueACOIRvAFi8+vygwzv/85/w3ntw4YVw59EVAIz9W35EF3f5xHsRVXo6NDQw86OrWZuU\nw93Z5YwfbxdN/fGP9vmi+UWk7mu5WMogTF7tiJ+4en+OzEzIzGxfvX1FcaHij42Jn5wM680AK/w/\n/BDysdyxvUWs+A9x7KOxai0fbTwkghb75qgvy/mLuRMgpATfzzxjn0tKsMtShw7l0f87OHorDD1x\n96BvKnMvvWRKHnTfur7pM4wfb1e4jh8fXnaxqOJ5J7BpE2zapG4VpV2i4g8kJUHfvvBtQ/jpHN2h\nBt54w6YF7LF9LUmNDfQaPjhC1vrHFBWRIV6hA2prbSKONhCXPh1/PJxwTD188AGcdVbYiWeCpqgI\ndu9uuc8zvd5tt9lnfPvUS1Y7YmaqohzIqPi7yM6GVXXhi787HscNN9jtz//PJtoddkH0e/6T/PSG\nJ7XSG3bfqSxf3o0vvoCbb4bP/vYp7NwJZ54ZdZu98dejb7qTeeKJ5tys3sRikFdROggq/i6ys2H5\nzsgkcs/Pb/bDXzHcFXD8kOiL/4QlvmfNTFjivzfsvlN5/nkHnTvDSSvKOeSX51GSh/WxRDuRtRf+\nZv60hvrUFSV4VPxdZGfD0G/tIGfJyz8PKxeu0wlz5tjylgUraUhNg/79I2Spfzo/sn+kzIa0DDo/\n4r83nJ9vJ6nMn5/JHVnlDHukkJ6yzfasq6uDGjOICF4zf3zdBfiM6Kg+dUUJChV/F6M3l1O69VYg\ntMFSN+5593l5dhzhZyeuZGVjLs4PYnCqCwpIfmY6GzvbcLQNJpnkZ9ruDX/1FYDhlm+LyPAONxzg\nmEHE8Jox4+suQP38ihI+Kv4uLvivjzjrIQife959XR3k5kKfnavoPeKQ2MyYAZz9Cziyyxq2zJtO\nsjQwv/7EVutv3Wrn9qekNPJsXjuZQeMx86ehrfwH6udXlJBQ8XfxyHGRET73vPvVqyE3R2DlSrJG\nDI7JjBn3XcfMmXDCb+xg7Wu3vt/qPP3bboNdu+Cuu75msh9/e9x61gUFLP/Nb5rnzTsccOutLbfV\nz68oIRGW+Btjehtj3jPGfON67uWjznHGmHnGmMXGmC+NMT8Lp81o4W+wNBThq6+HtWvh6IN+sPM9\nYzDYCy1X+46ZMJi1SQ4u6Px+i7sO9+wegO++g9des9M7d+xIhSlT9j9onHvWG0eNahmC4PHHNSSB\nokSAcHv+E4D3RWQo8L5r25ta4BoRORI4F3jUGNMzzHYjzg93lbILr6X5IQpfdbVdK3Zk59jN9IGW\nq33/+rih6pAzOXaLk9FnNwJW9FeubE55OGWKnd9fVQWHH76zeVC6Tx/tWStKBydc8b8EeN5Vfh64\n1LuCiHwtIt+4yuuAjUCfMNuNOGk3FHAz09ndJcsKX79+IQvf6tX2+RBj5/jbGA+xZcQIKF9/Fr3Z\nyu/GfM4jj1jRv/JKe3dw0UX246Wmunv/2+D//s9e8LRnrSgdnnDFP1tE1gO4nvu2VtkYcxLQCZvr\nt12RlQWvJBXw3M/epqQSu5goROFb5erw99/tKuTmRsTGYMjPh19dvxWAwwcO57Jf5/DkyHLy8uD1\n162fH6zPPz8fe6vyj3/A+ee3HpxMUZQOQUpbFYwxc4GDfLwU1DQYY0w/4EXgWhFp9FOnECgEyM7O\nprKyMpgmWlBTUxP0+3v0+DH/XpvGrcA3773Hdz16hNR2ZWUuSUmDqFv8X/ZkZTHv44/DsisU+s6d\ny2HT/wjYGDmTK6vo83ohv+i1lye2X09yciOXX17NM88cxMCBSzh576ewYQNLhg1jYwzsC5RYna9g\nUbuCQ+0KjpjYFUjQf38PYDnQz1XuByz3U6878D/gikCPHctkLm6OPlrk4osabaKNX/0q5LavusqV\nR+S000ROPz1su0LCT5ap1TgkNVVk7lxbraJCpLBrmdSlZdikHgcfHNWkHsHSEZNtRBO1Kzg6ol3E\nKJnLLOBaV/la4E3vCsaYTsA/gBdE5NUw24sq2dmwYaOxq3vXrAn5OEM+KWf+9zmUpPwH/ve/mIdI\nAP9xfp7Nq6JzZxvMDiB/XTmP7yskbU+tXdz27bexX9WrKErMCVf8HwLONsZ8A5zt2sYYc6Ix5mlX\nnTHASOA6Y8znrsdxYbYbFbKzYcMGwhP/8nLuXVHIQXtcwcd27YqLmN79ue+pq3d/7uCNN5pn/FBU\nRPKeOK/qVRQl5oQl/iKyWUTOEpGhructrv2fishNrnKZiKSKyHEej88jYXykcYu/DHKELP6NE9tB\niASgx198ZMfq1Ikefyltijy6YEE7jouvKEpU0RW+HmRn21Dye/vn2LgHO3YEfYzJh7QTMfXOKpWU\nBPX1TXHx89fZ5CglX/RqX6t6FUWJCSr+HrjTOW7tkWMLVcEL9r2L21GIBM/sWCkp0NjYHLRu3DhK\n8o29yLkyZzWh8XIUpcOj4u+BW/w3dM6xhRBcP/8+L3IrhSNGURHs3dtyn0hzQpTkZPZ2766rehUl\ngWhznn8i4Rb/6pQcjoWQxP/drAJeSoHnkm5gUuVeK6alpXEV05LcKiZf37zdFAPfTX09jZ07U+Lc\nHkuzFEWJI9rz98At/mt394HOnUNy+6xeDQuGXIXJyKBk2K3tIkSCd3YsT9xjEr8ftiG2RimKEle0\n5+9BH1fEoQ0bje2xh9DzX7UKjhvwAyzbBocfHlkDQ6W01E43rbWzkNwXAVPSXK5zX/kURUkItOfv\nYupU+M9/IDOzea7/jq/WMHVq4McQsT3/Ed2X2R3tRfy9Z/74GOBdddNN8bBMUZQ4oeLvwp3IvGtX\nK/7fpeaw75s1jBjR9nvdMfLds0OPTLbi/4SznYg/NM/8cQq8+OJ+CVE2jhoVbwsVRYkhKv4u3Auf\n1q2DTz6BZ97PIVM2kz+ips33ui8cr7qCV3SrXsYuMjh81MAoWx0inheCdjAmoShK7FHx9yA/H4YN\nsxmuDh/tmpcfwKCv+8Jxzz12u+bTZTQOOYz8s/T0KorSPlF18sDphBUrbPnvFTm2EOCgb34+DB9u\nyyd0WUa3Ew+LuH2KoiiRQsXfhTv5+R132O0bpuQA8PW7awJ+/7x50KfrbnptX8PqtHbk71cURfFC\nxd+FO/n5+efb7cOq30eA8i9vt1E+W4nK6b5w5ObCJUd8QxLC714/vClRuqIoSntDxd+FO/l5bi6M\npZwj/l8hBpsFi6qqVsMyL1gAr7wC69fDyT3sTJ+bHzmcBQtiZb2iKEpw6CIvL/r3h99TROpeP2GZ\nfcyMGT/eTg/dvt01zdMYTrpqKCdpKlxFUdopKv5eJCXB3/OqbI/fhTsWzqTKKkr8vG/5cvvsqFtm\n585rEnRFUdoxYYm/MaY38AqQA6wBxojIVj91uwNLgX+IyO3htBttbpvvYHKlneLpGQIBh/+wzG7x\nz/xhWftZ2asoiuKHcH3+E4D3RWQo8L5r2x8PAh+E2V5M+OcppdQGGZY5441yqnDwuz4L4b//1Ry4\niqK0a8IV/0uA513l54FLfVUyxgwHsoF3w2wvJmw6p4CbmE7jwINtCIQePVqPcV9ezk/nFDKItdZd\ntHOnJkFXFKVdE674Z4vIegDXc1/vCsaYJGAacE+YbcWM3Fx4mQIWz15LyQ9HwsiRrYdAKCoivSH+\neXsVRVECpU2fvzFmLnCQj5cCVbbbgNki8q3xjia5f1uFQCFAdnY2lZWVATaxPzU1NSG/f8uWbsBw\nZs36ir79+tFj/nzmt3IsZ24VU3wkSymurCLf633h2BVN1K7gULuCQ+0KjpjYJSIhP4DlQD9XuR+w\n3EedcmAtdkB4E7ADeKitYw8fPlzCwel0hvzejRtFQOTRR0XkD3+wGz/84Lf+3v4OWweEEprK4nBE\n1K5oonYFh9oVHGpXcIRjF/CpBKDf4bp9ZgHXusrXAm/6uLgUiMggEckBfgO8ICKtDQzHnaws6NLF\nxubn+OPtzi++8Fv/i5+1w7y9iqIorRCu+D8EnG2M+QY427WNMeZEY8zT4RoXL4yxfv/Vq4Fjj7U7\nP//cb/2KgwqYTDGAJkFXFOWAIKx5/iKyGTjLx/5Pgf1SQ4nIc8Bz4bQZK5rEPysLBg5sVfyXL4fe\n3XrCTih56hsYMiR2hiqKooSAxvbxg1v8RYDjjoOFC/3WXb4cRnX+L/TtC4ccEjsjFUVRQkTF3w+5\nuVBTA5s3Y8V/2TLYvdtn3eXL4YTd/4XTTts/P66iKEo7RMXfD7m59rlp0LehARYvbnrdnbd3yxZI\n2bSePjtXsSL71KASviuKosQLFX8/tBB/VyrHkntGNMX2d+ftfeklOJX/AvCLl04NKOG7oihKvFHx\n90NOjn3u/Ho53H8/AJPzaIrtn7+unJkzYeJEK/67SWfizOPJz4+XxYqiKIGj4u+DqVPhs8+gd2/4\n8b+KbKgGT1yhG/LXlbN8Tw7b8h4lJVnI++HV+BisKIoSJCr+PnC7dLKy4M8nVmFKmkM2uMsluVU0\n3FhI/3029n9qwx4abtRgboqiHBio+PsgP9/m812zBm78wIGUNMf0d5eLP0wmeU/LO4LkPbXsvluD\nuSmK0v5R8fdDfj7k5cEEKWVvilfohs6dmTyywecdwUNHVMXUTkVRlFDQNI5+cDptYvatFJCeDH/t\nXcQkV3YvzjiDyUuXBp3tS1EUpb2gPX8fOJ3W5//aazBgAKz6UQGDGtdwRrHA6afDO+9QkluFeL9R\ng7kpinKAoOLvgwULrM//zDPhnHPgyy/h5Zeh5m/l9kXstE8DCEaDuSmKcsChbh8fjB/fXD7nHHj2\nWejeHUbNL4K6uhZ1DULJaocdHVYURTlAUPFvg1GjbLied9+F2blVTPaRsWtSZRUl8TBOURQlRFT8\n2yArC044wYr/P79yUOJjkHd7Tx3kVRTlwEJ9/gFwzjkwbx6svLGUWq+MXbVkUH2bDvIqinJgEZb4\nG2N6G2PeM8Z843ru5afeIGPMu8aYpcaYJcaYnHDajSVTp0KfPlBfD9UjC1h933TWGgfFlbDWOFh9\n33SOLNVBXkVRDizCdftMAN4XkYeMMRNc2/f6qPcCUCoi7xljugKNYbYbM9yhHtLTbeSG2toC3pIC\nqIR9E+F32ulXFOUAJFy3zyXA867y88Cl3hWMMUcAKSLyHoCI1IhIrXe99oo71ENjo31+6y3o1Anu\nvBOeesquCVAURTnQCFf8s0VkPYDrua+POocC24wxrxtjFhpj/mCMSQ6z3ZiSn988fT89Hd55Bx59\n1F4MxozRC4CiKAceRmS/daotKxgzFzjIx0tFwPMi0tOj7lYRaeH3N8ZcDjwDHA+sBV4BZovIMz7a\nKgQKAbKzs4fPmDEjuE/jQU1NDV27dg35/Z4sXNiTyZOP4JBDali+vBsPPriY44/f1vTasmXdGDv2\n25jbFUnUruBQu4JD7QqOcOzKz8//TERObLOiiIT8AJYD/VzlfsByH3V+BFR6bI8D/trWsYcPHy7h\n4HQ6w3q/m4oKkaws++xrO152RRq1KzjUruBQu4IjHLuATyUA/Q7X7TMLuNZVvhZ400edBUAvY0wf\n1/aZwJIw240Z7lAP7gxd7jEAV5QHRVGUA5JwZ/s8BMw0xtyIdelcAWCMORG4RURuEpEGY8xvgPeN\nMQb4DHgqzHZjhmeoBzf5+Wi6RkVRDmjCEn8R2Qyc5WP/p8BNHtvvAceE05aiKIoSOXSFr6IoSgKi\n4q8oipKAqPgriqIkICr+iqIoCUibi7zihTHmByCcbOhZwKYImRNJ1K7gULuCQ+0Kjo5ol0NE+rRV\nqbS8vtYAAAUTSURBVN2Kf7gYYz6VQFa5xRi1KzjUruBQu4Ijke1St4+iKEoCouKvKIqSgHRk8Z8e\nbwP8oHYFh9oVHGpXcCSsXR3W568oiqL4pyP3/BVFURQ/HNDib4w51xiz3BizwpVG0vv1NGPMK67X\nP45F7mBjzMHGGKcrX/FiY8ydPurkGWO2G2M+dz2Ko22XR9trjDFfudr91MfrxhjzZ9c5+9IYc0IM\nbDrM41x8bozZYYy5y6tOTM6ZMebvxpiNxphFHvsCzVV9ravON8aYa33VibBdfzDGLHN9T/8wxvT0\n895Wv/Mo2FVijPnO47s63897W/3/RsGuVzxsWmOM+dzPe6N5vnzqQ1x+Y4HEfW6PDyAZWAkMBjoB\nXwBHeNW5DXjSVb4SeCUGdvUDTnCVuwFf+7ArD3grTudtDZDVyuvnA28DBpuL4eM4fK/fY+cqx/yc\nASOBE4BFHvumAhNc5QnAwz7e1xtY5Xru5Sr3irJd52BTpAI87MuuQL7zKNhVAvwmgO+51f9vpO3y\nen0aUByH8+VTH+LxGzuQe/4nAStEZJWI7AVmYHMKe+KZY/g14CxXWOmoISLrReR/rvJOYCkwIJpt\nRphLgBfEMh/oaYzpF8P2zwJWikg4C/xCRkQ+BLZ47W4zVzUwGnhPRLaIyFbgPeDcaNolIu+KSL1r\ncz4wMFLthWNXgATy/42KXS4NGAO8HKn2AqUVfYj5b+xAFv8BgGfuxGr2F9mmOq4/yXYgMybWAS43\n0/HAxz5ePsUY84Ux5m1jzJGxsgkQ4F1jzGfGps30JpDzGk2uxP+fMl7nLJBc1fE+bzdg79h80dZ3\nHg1ud7mj/u7HhRHP83U6sEFEvvHzekzOl5c+xPw3diCLv68evPfUpUDqRAVjTFfg/4C7RGSH18v/\nw7o1jgUeA96IhU0uThWRE4DzgF8YY0Z6vR7Pc9YJuBh41cfL8TxngRDP81YE1APlfqq09Z1HmieA\nQ4DjgPVYF4s3cTtfwFha7/VH/Xy1oQ9+3+ZjX8jn7EAW/2rgYI/tgcA6f3WMMSlAD0K7RQ0KY0wq\n9ostF5HXvV8XkR0iUuMqzwZSjTFZ0bbL1d461/NG4B/Y229PAjmv0eI84H8issH7hXieM2CD2/Xl\net7oo05czptr0O9CoEBcjmFvAvjOI4qIbBCRBhFpxGbt89VevM5XCnAZ8Iq/OtE+X370Iea/sQNZ\n/BcAQ40xua4e45XYnMKeeOYYvhyo8PcHiRQuf+IzwFIRecRPnYPcYw/GmJOw38PmaNrlaquLMaab\nu4wdMFzkVW0WcI2x/AjY7r4djQF+e2TxOmcuAslVPQc4xxjTy+XmOMe1L2oYY84F7gUuFpFaP3UC\n+c4jbZfnGNFP/LQXyP83GowClolIta8Xo32+WtGH2P/GojGiHasHdmbK19hZA0WufVOwfwaAdKwL\nYQXwCTA4Bjadhr0V+xL43PU4H7gFm9cY4HZgMXaGw3zgxzE6X4NdbX7hat99zjxtM8BfXef0K+DE\nGNmWgRXzHh77Yn7OsBef9cA+bE/rRuw40fvAN67n3q66JwJPe7z3BtdvbQVwfQzsWoH1Abt/Z+6Z\nbf2B2a1951G260XXb+dLrKj187bLtb3f/zeadrn2P+f+TXnUjeX58qcPMf+N6QpfRVGUBORAdvso\niqIoIaLiryiKkoCo+CuKoiQgKv6KoigJiIq/oihKAqLiryiKkoCo+CuKoiQgKv6KoigJyP8HUqeo\nD0VSDfYAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fc7f3877690>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "0.071881806500767961"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "plt.plot(t, x, '-bx')\n",
    "plt.plot(t, yr, '-ro')\n",
    "window = np.zeros(N)\n",
    "window[:15] = 1\n",
    "yr=idct(y*window, norm='ortho')\n",
    "plt.plot(t,yr,'g+')\n",
    "plt.legend(['x', '$x_{20}$', '$x_{15}$'])\n",
    "plt.grid()\n",
    "plt.show()\n",
    "sum(abs(x-yr)**2) / sum(abs(x)**2)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}