{
"cells": [
{
"cell_type": "markdown",
"metadata": {
"nbsphinx": "hidden"
},
"source": [
"# The Discrete Fourier Transform\n",
"\n",
"*This Jupyter notebook is part of a [collection of notebooks](../index.ipynb) in the bachelors module Signals and Systems, Comunications Engineering, Universität Rostock. Please direct questions and suggestions to [Sascha.Spors@uni-rostock.de](mailto:Sascha.Spors@uni-rostock.de).*"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Fast Fourier Transform\n",
"\n",
"The discrete Fourier transformation (DFT) can be implemented computationally very efficiently by the [fast Fourier transform (FFT)](https://en.wikipedia.org/wiki/Fast_Fourier_transform). Various algorithms have been developed for the FFT resulting in various levels of computational efficiency for a wide range of DFT lengths. The concept of the so called [radix-2 Cooley–Tukey algorithm](https://en.wikipedia.org/wiki/Cooley%E2%80%93Tukey_FFT_algorithm) is introduced in the following."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Radix-2 Decimation-in-Time Algorithm\n",
"\n",
"Let's first consider the straightforward implementation of the DFT $X[\\mu] = \\text{DFT}_N \\{ x[k] \\}$ by its definition\n",
"\n",
"\\begin{equation}\n",
"X[\\mu] = \\sum_{k=0}^{N-1} x[k] \\, w_N^{\\mu k}\n",
"\\end{equation}\n",
"\n",
"where $w_N = e^{-j \\frac{2 \\pi}{N}}$. Evaluation of the definition for $\\mu = 0,1,\\dots,N-1$ requires $N^2$ complex multiplications and $N \\cdot (N-1)$ complex additions. The numerical complexity of the DFT scales consequently [on the order of](https://en.wikipedia.org/wiki/Big_O_notation) $\\mathcal{O} (N^2)$.\n",
"\n",
"The basic idea of the radix-2 decimation-in-time (DIT) algorithm is to decompose the computation of the DFT into two summations: one over the even indexes $k$ of the signal $x[k]$ and one over the odd indexes. Splitting the definition of the DFT for an even lengths $N$ and rearranging terms yields\n",
"\n",
"\\begin{align}\n",
"X[\\mu] &= \\sum_{\\kappa = 0}^{\\frac{N}{2} - 1} x[2 \\kappa] \\, w_N^{\\mu 2 \\kappa} + \n",
"\\sum_{\\kappa = 0}^{\\frac{N}{2} - 1} x[2 \\kappa + 1] \\, w_N^{\\mu (2 \\kappa + 1)} \\\\\n",
"&= \\sum_{\\kappa = 0}^{\\frac{N}{2} - 1} x[2 \\kappa] \\, w_{\\frac{N}{2}}^{\\mu \\kappa} +\n",
"w_N^{\\mu} \\sum_{\\kappa = 0}^{\\frac{N}{2} - 1} \\, x[2 \\kappa + 1] w_{\\frac{N}{2}}^{\\mu \\kappa} \\\\\n",
"&= X_1[\\mu] + w_N^{\\mu} \\cdot X_2[\\mu]\n",
"\\end{align}\n",
"\n",
"It follows from the last equality that the DFT can be decomposed into two DFTs $X_1[\\mu]$ and $X_2[\\mu]$ of length $\\frac{N}{2}$ operating on the even and odd indexes of the signal $x[k]$. The decomposed DFT requires $2 \\cdot (\\frac{N}{2})^2 + N$ complex multiplications and $2 \\cdot \\frac{N}{2} \\cdot (\\frac{N}{2} -1) + N$ complex additions. For a length $N = 2^w$ with $w \\in \\mathbb{N}$ which is a power of two this principle can be applied recursively till a DFT of length $2$ is reached. This comprises then the radix-2 DIT algorithm. It requires $\\frac{N}{2} \\log_2 N$ complex multiplications and $N \\log_2 N$ complex additions. The numerical complexity of the FFT algorithm scales consequently on the order of $\\mathcal{O} (N \\log_2 N)$. The notation DIT is due to the fact that the decomposition is performed with respect to the (time-domain) signal $x[k]$ and not its spectrum $X[\\mu]$.\n",
"\n",
"The derivation of the FFT following above principles for the case $N=8$ is illustrated using signal flow diagrams. The first decomposition results in\n",
"\n",
"\n",
"\n",
"where $\\underset{a}{\\oplus}$ denotes the weigthed summation $g[k] + a \\cdot h[k]$ of two signals whereby the weighted signal $h[k]$ is denoted by the arrow. The same decomposition is now applied to each of the DFTs of length $\\frac{N}{2} = 4$ resulting in two DFTs of length $\\frac{N}{4} = 2$\n",
"\n",
"\n",
"\n",
"where $w_\\frac{N}{2}^0 = w_N^0$, $w_\\frac{N}{2}^1 = w_N^2$, $w_\\frac{N}{2}^2 = w_N^4$ und $w_\\frac{N}{2}^3 = w_N^6$. The resulting DFTs of length $2$ can be realized by\n",
"\n",
"\n",
"\n",
"where $w_2^0 = 1$ und $w_2^1 = -1$. Combining the decompositions yields the overall flow diagram of the FFT for $N=8$\n",
"\n",
"\n",
"\n",
"Further optimizations can be applied by noting that various common terms exist and that a sign reversal requires to swap only one bit in common representations of floating point numbers."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Benchmark\n",
"\n",
"The radix-2 DIT algorithm presented above can only be applied to lengths $N = 2^w$ which are are power of two. Similar and other principles can be applied to derive efficient algorithms for various other cases. A wide variety of implemented FFTs is available for many hardware platforms. Their computational efficiency depends heavily on the particular algorithm and hardware used. In the following the performance of the [`numpy.fft`](http://docs.scipy.org/doc/numpy/reference/generated/numpy.fft.fft.html#numpy.fft.fft) function is evaluated by comparing the execution times of a [DFT realized by matrix/vector multiplication](definition.ipynb#Matrix/Vector-Representation) with the FFT algorithm. Note that the execution of the following cell may some time."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAA40AAAJsCAYAAABZO0aqAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3XeYJHd16P3vEQMIkMSKNDKStXvFJYMsuESTmmBy8MVg\nMpJxvCZj4ALXtvTACwbuS7Rf+xJ1SSJjMoigLTBgkaQREhI44JFE2CVpBLIwSNZ5/6gabc9sd+/M\n9HT/qqu/n+fpZ/pUVVefnvlt7Zyp36mKzESSJEmSpEEOKp2AJEmSJKm9LBolSZIkSUNZNEqSJEmS\nhrJolCRJkiQNZdEoSZIkSRrKolGSJEmSNJRFoyQVFhGPjYhPTmC/B0fERyJiJSLevd37n4aI+HhE\nPGFK73VORNx9Gu815P13RsQVEeH/zS0TESdHxAtL5yFJpfgfkyRtUkQsR8SlEfGziPh+8wvlNTf4\n2v0Kg8w8JTPvP4FUHwFcHzg8Mx81gf1vq4g4MSLe2r8sMx+YmW+bwHvtVwRk5q0y8/Pb/V6btKWb\nJ89CwRkR94iI/2z+3fwsIi6IiHdHxO3WbXdFRPy82ebnEfHT5g8rq8su7dvPzyPiZ6U+U5Pv8RFx\neX8+EfHaZt3/jYhfrlv3u33Pf9Z8lkv7lj2m5OeRpEFa+5+LJLVYAg/KzMOA44DbAM/f4GujeX1M\nKLd+O4F/yswtFSKaKdMcV+P4XmYe1vzbuRPwLeAfIuKefdskcGyz3aGZeZ3mDyuHNq97QN9+VpeV\n9qX+fDLzac3yBF62bt17+p4fBpxPfTxZXfbOch9DkgazaJSkrQmAzPwhcCp18ViviHhgRJwRERdH\nxPkRcWLf6z7XfF1pzircsTlT8Q99r//NiPhKRFwUEV+OiDsPTSLiZhGxu9n27Ih4SLP8JOAvgUc3\n7/N7A14bEfG8iPiXiPhRRLwrInY06343Iv41Ig5p4gdExA8i4rp97/upiPhJRJwXEY/s2+/BEfGK\n5ozsRRHx+Yi4enOm6cJ1OfxbRNwrIu4HvAB4VHPG5cxm/e6IeFJfvn/e7HdPcxbnsGbd6pm2Jzbf\n8x9GxAuGfM/+EHgc8Nzme/Oh/lya5ydGxHsi4m3NNmdFxI2b79fe5j3u07fPwyLijVGfeb4wIl4U\nEdGsu1FEVFFPE/5hRIwqCgL4/Yj4XvN41kZ+Xuw/ru7UfJ9u07z28c3352ZN/PsR8fcb2C/Nvr7Y\n/CzPjIh79K3bHREvjIgvNO/7yYi4zojPd6XM/H5mngi8EXjZuu/B2MVvRLw66rOZF0fEVyPirn3r\nToz6LOdbmrzPjojb9q2/TUR8vXntu4CDx81nVKq0v9iXNOcsGiVpDBFxFPWZj3/uW3wJ8ITMvDbw\nIOBPIuKhzbrVnrnDmrMKX27ibPZ3OPBR4NXAdYFXAR9rlq9/7wXgI8AnqaehPg14R0TcODNPAl4C\nvKt5n5MHpP904KHA3YAbAhcBfwuQme8BvgS8tikC3gg8KTN/EvVU3E8BbweuBzwG+NuIuHmz31dQ\nn329E3Ad4LnAFf2fc73MPLXJ993NGZfbDNjs94AnAvcAjgEOBf5m3TZ3AW4M3Af4y4i46YD3egPw\nDuDlzffmYYNyAh4MvAXYASxR/3EgqL9XLwJe37ftW4FfNXndBvgt4A+adS8CTs3MHcBRwF8Peb9V\nPeBGwP2A560Wsoz4ebH/uDodqJp90bzmX6m/d6vbVwfab0QcST0eX5iZhwPPBt4fzR8PGo8Bjqce\ng1dvttmMDwC3jYhrbPJ1B/IV4FjgcOAU4L0RcbW+9Q9pll+b+t/R/wcQEVcF/p76Z38d4L3A72xz\nbpI0UywaJWlrPhh1L9UFwF7gpNUVmfn5zPxm8/wc4F3s+2V91bAzCw+inlJ6SmZekZnvop7C95AB\n294JuFZmviwzL8/M3dS/4G+0J+qPgP+VmT/IzMuAFwKPiH19cU8B7k1dXHwoMz/RLH8w8G+Z+das\nLQHvb14b1MXd0zJzT7P+9Gb/43os8MrMPD8zL6WeEvzovnwTOCkzf5WZ3wDOAn5jjPf7h8z8TGZe\nQV04XA94aWb+J/XPdGdzhnERuD/wzMz8j8z8MXXR/+hmP5c12x7Z5PalA7zvSc1+zgFOZt/Pc9TP\na3U89Y+rz7Nv3N0N+Ku++B7sKxpH7fdxwMeaop7M/CzwNeCBfe9zcmb+a2b+EngPfWfdN+j7Td47\n+pad0ZzZ/GlEvHqT+6PJ9ZTMXGn+Hb2KuqDt/yPCFzLz1Gb69tuoC0yAOwMLmfnazPzPzHw/8NUD\nvN2dm1xXc75D37rn9K374VY+iySVZtEoSVvzsKYf6R7AzagLCgAi4g4RcVozFXEF+OP+9QdwQ+oe\np37nA0cO2fbCDW47yE7g75tfaH8KnEtd4CwCZObF1MXSLYFXrnvdnVZfFxEXURd0i9Sf82DgOxvM\nYTPWf2/OBxZW823s7Xt+KXDIGO/Xv69fAD/u6w/9BXWhcwhwNHBV4Ad934//Q33mDeA51P/ffqWZ\nBrnfVOE+CXy3Lz6f+nPD6J/XoDO4nwPu1hS1BwHvBu4aETupz0ietYH97gR+d93P+i7AEX3vs6fv\n+Va+50c2+a/0LbtNZh7e9DM+Y5P7AyAi/iwizm2KtYuAw1j773B93gc3hfKvAd9bt7v1/ybX+8cm\n19Wcv9K37n/3rbvBVj6LJJVm0ShJW7Pa0/gP1NPYXtG37hTgg8CRzZTE17HvDNCBLkrzfWDXumVH\ns/8vsavb/voGtx3kAuABzS+0q7/UXiszfwAQEccBTwLeydoplRcC1brXHZaZTwF+TF1Q3WjA+/07\ncOVVZiPiKuwrrGBj35udffFO6uJm7+DNR9rOiwNdCPwHcN2+78eOzDwW6r7XzPyjzDwS+BPqqbzH\njNhf/8/0aOrPDaN/Xvt9nsz8V+qfxdOAz2fmv1MXSn8EfKFv01H7vRB467p1h2bm/978t2mohwNn\nZOYv+paN1ePX9C8+F3hEk/PhwM82uN8fsP8fXo4eJx9JmnUWjZI0vlcDvxURq9PbDgEuyszLmmlq\nj+3b9kfU/X2DiiqAjwM3johHR8RVIuJRwM2pp52u92Xg3yPiuRGxEBE96qmjG7364uuAl0TE0QAR\ncf3V3suIOJh6yt7zqAvHG0bE/2he91HgJs3FVRYi4qoRcbuIuGlzJu5k4JUR8WsRcVBzIZWrAv9E\nfTbnAU0/5p8D/T1me4FdzRTXQd4JPDMidkV9gZ4XU/dsrvZLbqbQ2Evdfzi2zNxD3eP5qog4NGrH\nRHPPx4h4RNMbCPXZtCuA/xyyuwD+IiKuERG3pJ7q+65m3dCfF8PH1eeopxmvXiinWhcfaL9vBx4S\nEfdtfpYHR31BoxuyNVf+jCLihlFfJOpJbPzqwxt1KPUfFH4SEVeLiL9slm0kt38ELo+Ipzb/Bh8O\n3GHE6ySp8ywaJWnz1pzVaXrY3gL8RbPoycCLIuJi6sLo3X3b/oK62PnigN4nMvOn1IXfs6nP2j2b\n+nL8P90vibr/7KHU/WU/pr4ozBMy85/XbzvEa4APAZ9qcv0S+345fglwQWa+PjN/BTyh+Uw3ysxL\ngPtS9+x9v3m8lLpnjCbns6n7wH7SrDsoM38G/CnwJuopmD9n7VTM91L/4v6TiPja6sfsW/9m6kL2\n89QXdbmU+iwaA7YdFPd7E3DL5mfwgQ1sP0j/9k+kLoDPBX5K/VlWp3DeHvhy0wP7Qep+z2HTHZO6\noPsX4NPUF+v5bLNu6M9rxLj6HPUfMT4/JD7Qfr8LPIz6yrY/op6m+Wz2/f6w2e/Zr0Vzz0LqC9Xc\nErhH32fcyj4HOZX6AlH/BPwb9VhZP5V7vYQr/109nLpg/ynwSOqe3a3YyGfxljiSWi/S23dJkiRJ\nkobwTKMkSZIkaSiLRkmSJEnSUBaNkiRJkqShFkonMGkRYdOmJEmSpLmWmVu+ndFcnGnMzJl/nHji\niZ14z+3Y51b2sZnXbHTbA213oPXHH3/81H+mk3iUGJuTet9x9zkrY/NA2xzfkePmdvxM2/KeJcbm\nZl/nsXO6P9M2va/HzvrRlbG5HT/Ttrxn28fmRrffjm3GNRdFYxf0er1OvOd27HMr+9jMaza67YG2\nK/EzK6HU52zj+JyVsbnZ951lHjvHe73Hzsnx2Dne6z12TpbHzq2/frOv2a5xN+mfWedvuRER2fXP\nqNl00kkncdJJJ5VOQ9rPSRGc5HFTLeWxU23l2FSbRQTp9FRp9szLXys1e3qlE5BG8NiptnJsqsss\nGiVJkiRJQ1k0SpIkSZKGsqdRkrRWBHjclCSpM+xplCRJkiRNjEWjVEhVVaVTkAaqSicgjeCxU23l\n2FSXWTRKkiRJkoayp1GStJY9jZIkdYo9jZIkSZKkibFolAqx90FtVZVOQBrBY6fayrGpLrNolCRJ\nkiQNZU+jJGktexolSeoUexolSZIkSRNj0SgVYu+D2qoqnYA0gsdOtZVjU11m0ShJkiRJGsqeRknS\nWvY0SpLUKfY0SpIkSZImxqJRKsTeB7VVVToBaQSPnWorx6a6zKJRkiRJkjSUPY2SpLXsaZQkqVPs\naZQkSZIkTcxC6QSmIWLLRbUkzZ3dwD09bkqSNDWLizvZs2e5dBpDzUXRCE6zUhtVQK9wDtIggcdN\ntVeFx061U4VjU1u1d2+7/1g7Fz2N/vIjSRuXBOFxU5KkKQomWZfZ0yhJkiRJmhiLRqmYqnQC0kBV\n6QSkkarSCUhDVKUTkCbGolGSJEmSNNTUi8aIeGZEnBMR34iId0TE1SPiyRHxzxHxnxFxnb5td0TE\nByLirIg4PSJu0Sw/KiJOi4hzI+LsiHjatD+HNL5e6QSkgXqlE5BG6pVOQBqiVzoBaWKmWjRGxA2B\npwK3zcxjqa/e+ijgC8C9gfPXveQFwJmZ+RvA8cBrm+WXA8/KzFsAdwaeHBE3m8JHkCRJkqS5UmJ6\n6lWAa0XEAnBN4PuZeVZmXkB9nfd+twA+C5CZ3wZ2RcT1M3NPZi41yy8BzgOOnNonkLZFVToBaaCq\ndALSSFXpBKQhqtIJSBMz1aIxM78PvAK4APgesJKZnxnxkrOAhwNExB2Ao4Gj+jeIiF3AccCXtz9j\nSZIkSZpvC9N8s4jYATwM2AlcDLwvIh6bmacMeclLgddExBnA2cCZ1FNTV/d3CPA+4OnNGcchTgB2\nNc93UNeYvSaumq/GxtOOey3Lx9i4WhO1JR9jY2Pj2Yk5wHpj41FxE1V13Ov1thwvLS2xsrICwPLy\nMuOKSd5Ecr83i3gEcL/M/MMmfgJwx8x8ShN/B7hdZv50yOv/Dbh1Zl7STG/9KPCJzHzNiPdMvEm1\nJG1YEoTHTUmSpiiYZF0WEWTm+lbADTtoO5PZgAuAO0XEwRER1Be/Oa9vfdDX1xgR146IqzbP/xD4\nXN8ZxTcD544qGKV2q0onIA1UlU5AGqkqnYA0RFU6AWliplo0ZuZXqKeTnkndrxjA6yPiqRFxIfXF\nbM6KiNc3L7k58M2IOBe4H/B0gIi4C/A44F4RcWZEnBER95/mZ5EkSZKkeTDV6aklOD1VkjbH6amS\nJE2b01MlSZIkSTPKolEqpiqdgDRQVToBaaSqdALSEFXpBKSJsWiUJEmSJA1lT6MkaQ17GiVJmjZ7\nGiVJkiRJM8qiUSqmKp2ANFBVOgFppKp0AtIQVekEpImxaJQkSZIkDWVPoyRpDXsaJUmaNnsaJUmS\nJEkzyqJRKqYqnYA0UFU6AWmkqnQC0hBV6QSkibFolCRJkiQNZU+jJGkNexolSZo2exolSZIkSTNq\noXQC07HlolqS5k4FeNyUJGl6Fhd3lk5hpLkoGrs+BVezqaoqer1e6TSk/VQx2Sky0jg8dqqtHJvq\nsrnoaez6Z5SkbRUBHjclSeoMexolSZIkSRNj0SgVUlVV6RSkgarSCUgjeOxUWzk21WUWjZIkSZKk\noexplCStZU+jJEmdMm5P41xcPTXCS8dL0kYlHjclSbNrcXEne/Ysl06jU+aiaKx/BZLapgJ6hXOQ\n9lcReNxUe1V47FQ7VTg222HvXv/wud3saZQkSZIkDTUXPY3+xVySNi4JwuOmJGlmBV2vcTbL+zRK\nkiRJkibGolEqpiqdgDRQVToBaaSqdALSEFXpBKSJKVI0RsRBEXFGRHy4iXdFxOkR8e2IeGdELDTL\nj46Iz0TEWRFxWkTcsG8fvx4Rp0bEuRFxTkQcXeKzSJIkSVKXFelpjIhnAv8NOCwzHxoR7wbel5nv\njYi/A5Yy83UR8R7gw5n59ojoAU/KzCc2+9gNvCgzT4uIawJXZOZ/DHgvexolaRPsaZQkzTZ7Gteb\nuZ7GiDgKeCDwxr7F9wLe3zx/C/DbzfNbAKcBZGYFPKzZx82Bq2Tm6rpLBxWMkiRJkqTxlJie+irg\nOTSn/yLiusBFmXlFs/67wJHN8yXgd5rtHg4cEhGHAzcBLo6I90fE1yPiZeGdqDVzqtIJSANVpROQ\nRqpKJyANUZVOQJqYhWm+WUQ8CNibmUvNdFOAaB79Vs8nPwf4m4g4Afg88D3gcuq87wocB1wIvAc4\nATh58DufAOxqnu9oXrb69lXz1djY2Nh4rbbkY2zcH3OA9cbGpeKlluUzzzFUVUWv17vyOTBX8dLS\nEisrKwAsLy8zrqn2NEbES4DHUxd+1wAOBT4I3Bc4IjOviIg7ASdm5gPWvfZawHmZeXRE3BH4q8y8\nV7Pu8cAdM/OpA97TnkZJ2gR7GiVJs82exvVmqqcxM1+QmUdn5jHAo4HTMvPxwG7gkc1mxwMfgnrq\nat+00+cDb26efxU4vJnaCnVP5LnT+AySJEmSNE+mWjSO8DzgWRHxT8B1gDc1y3vAtyPiW8ANgBcD\nNP2PzwZOi4izmm3fMNWMpbFVpROQBqpKJyCNVJVOQBqiKp2ANDFT7Wnsl5mfAz7XPP834I4Dtnk/\n+66qun7dZ4HfmGSOkiRJkjTvityncZrsaZSkzbGnUZI02+xpXG+meholSZIkSbPFolEqpiqdgDRQ\nVToBaaSqdALSEFXpBKSJsWiUJEmSJA1lT6MkaQ17GiVJs82exvXsaZQkSZIkTYxFo1RMVToBaaCq\ndALSSFXpBKQhqtIJSBNj0ShJkiRJGsqeRknSGvY0SpJmmz2N69nTKEmSJEmaGItGqZiqdALSQFXp\nBKSRqtIJSENUpROQJsaiUZIkSZI0lD2NkqQ17GmUJM02exrXG7encWE7k2mvLX9/JGlOedyUJM2m\nxcWdpVPonLkoGv1Lg9qoqip6vV7pNKT9VOFfaNVeHjvVVo5NdZk9jZIkSZKkoeaip7Hrn1GStlUE\neNyUJKkzvE+jJEmSJGliLBqlQqqqKp2CNFBVOgFpBI+daivHprrMolGSJEmSNNRcXD01wkvHS9JG\nJR43JakLFhd3smfPcuk01AFzcSEcvEm1JG1YEoTHTUnqAG+hpJoXwpFmVlU6AWmgqnQC0khV6QSk\nIarSCUgTY9EoSZIkSRrK6amSpDWcnipJXeH0VNWcnipJkiRJmpipFo0RcVREnBYR50bE2RHxtGb5\n4RHxqYj4dkScGhHXXve620fE5RHx8L5lL4uIcyLimxHx6ml+Dml7VKUTkAaqSicgjVSVTkAaoiqd\ngDQx0z7TeDnwrMy8BXBn4MkRcTPgecBnMvOmwGnA81dfEBEHAS8FPtm37M7Ab2bmrYBbAXeIiLtP\n72NIkiRJ0nyYatGYmXsyc6l5fglwHnAU8DDgLc1mbwF+u+9lTwXeB/ywf1fAwRFxMHAN6vtN7p1s\n9tJ265VOQBqoVzoBaaRe6QSkIXqlE5AmplhPY0TsAo4DTgcWM3Mv1IUlcP1mmyOpC8j/A1zZuJmZ\np1PPAfgB8D3g1Mz89vSylyRJkqT5sFDiTSPiEOqzh0/PzEvqK5wO9Crgf2ZmRgQ0hWNE3Ai4GXDD\nZtlnIuLUzPzC4N2cAOxqnu+grlV7TVw1X42Npx2vPm9LPsbG1ZqlbcnH2HhtvLqsLfkYG6/GS8Az\nWpTPPlVVx71ez3hO4qWlJVZWVgBYXl5mXFO/5UZELAAfBT6Rma9plp0H9DJzb0QcAezOzJtHxHdW\nXwZcD/h34I+AmwBXz8wXN6//C+AXmfn/Dng/b7mhlqrYd3CX2mM3wT09bqq1Kjx2qp0q2jc2veWG\narN4y403A+euFoyND1OfDgQ4HvgQQGYe0zz+C/WZyT/NzA8DFwD3iIirRMRVgXtQ90dKM6RXOgFp\noF7pBKSReqUTkIbolU5AmpipTk+NiLsAjwPOjogzqU8BvgB4GfCeiHgSdUH4yAEv7/8zyfuAewFn\nA1dQn7X82CRzlyRJkqR5NPXpqdPm9FS1V4V/lVQbOT1V7VbhsVPtVNG+sen0VNVmcXqqJEmSJGlG\neKZRkrRGEoTHTUnqAM80quaZRkmSJEnSxFg0SsVUpROQBqpKJyCNVJVOQBqiKp2ANDEWjZIkSZKk\noexplCStYU+jJHWFPY2q2dMoSZIkSZoYi0apmKp0AtJAVekEpJGq0glIQ1SlE5AmxqJRkiRJkjSU\nPY2SpDXsaZSkrrCnUTV7GiVJkiRJE2PRKBVTlU5AGqgqnYA0UlU6AWmIqnQC0sRYNEqSJEmShloo\nncB0bHn6riTNnR7gcVOSZt/i4s7SKagj5qJotAFYkjYhvHCCJEnax+mpUiFVVZVOQRqoKp2ANILH\nTrWVY1NdZtEoSZIkSRpqLu7T2PXPKEnbKgI8bkqS1Bnep1GSJEmSNDEWjVIh9j6orarSCUgjeOxU\nWzk21WVzcfXUCC8dL0kbtRuPm5I0zOLiTvbsWS6dhjRVc9HTCN3+jJK0nZIgPG5K0hDelkizx55G\nSZIkSdLEWDRKxVSlE5AGqkonII1UlU5AGsieRnWZRaMkSZIkaSh7GiVJa9jTKEmj2NOo2TNTPY0R\n8aaI2BsR3+hbdnhEfCoivh0Rp0bEtZvlj42IsyJiKSK+EBG3XrevgyLijIj48DQ/gyRJkiTNk2lP\nTz0ZuN+6Zc8DPpOZNwVOA57fLP8OcPfMPA74f4A3rHvd04FzJ5irNGFV6QSkgarSCUgjVaUTkAay\np1FdNtWiMTO/AFy0bvHDgLc0z98C/Haz7emZeXGz/HTgyNUXRMRRwAOBN040YUmSJEmac224EM4N\nMnMvQGbuAa4/YJs/AD7RF78KeA42K2qm9UonIA3UK52ANFKvdALSQL1er3QK0sQslE7gQCLinsDv\nAXdt4gcBezNzKSJ6wAYaOk8AdjXPdwDHse8/nar5amxsbGy8VlvyMTY2Nm5XvDoVdbVQNDZuW7y0\ntMTKygoAy8vLjGvqV0+NiJ3ARzLz2CY+D+hl5t6IOALYnZk3b9YdC7wfuH9m/muz7CXA44HLgWsA\nhwIfyMwnDnk/r56qlqrY95+R1B67Ce7pcVOtVeGxU2UNvnpqVVWebVRrzdTVUxvB2rODH6Y+FQhw\nPPAhgIg4mrpgfMJqwQiQmS/IzKMz8xjg0cBpwwpGSZIkSdJ4pnqmMSJOof7z4HWBvcCJwAeB9wK/\nDlwAPDIzVyLiDcDDgfOpi8zLMvMO6/Z3D+DPMvOhI97TM42StAnep1GSRvE+jZo9455pnPr01Gmz\naJSkzbFolKRRLBo1e2ZxeqokYF9TvdQuVekEpJGq0glIA3mfRnWZRaMkSZIkaSinp0qS1nB6qiSN\n4vRUzR6np0qSJEmSJsaiUSqmKp2ANFBVOgFppKp0AtJA9jSqyywaJUmSJElD2dMoSVrDnkZJGsWe\nRs0eexolSZIkSRNj0SgVU5VOQBqoKp2ANFJVOgFpIHsa1WUWjZIkSZKkoexplCStYU+jJI1iT6Nm\njz2NkiRJkqSJWSidwHRsuaiWpLlTAR43JWmwxcWdA5dXVUWv15tuMtKUzEXR6BQCtZH/uaitqnDq\nldrLY6ckTd9c9DR2/TNK0raKAI+bkiR1hj2NkiRJkqSJsWiUCvF+TmqrqnQC0ggeO9VWjk11mUWj\nJEmSJGkoexolSWvZ0yhJUqeM29M4F1dPjfDS8ZK0UYnHTXXb4uJO9uxZLp2GJM2MuSga61+BpLap\ngF7hHKT9VQQeN9VeFeMeO/fu9Y8i2n7eDkZdZk+jJEmSJGmouehp9C/mkrRxSRAeN9VpQdd//5Gk\nft6nUZIkSZI0MRaNUjFV6QSkgarSCUgjVaUTkAbyPo3qMotGSZIkSdJQrelpjIhl4GLgCuCyzLxD\nRDwCOAm4OXD7zDyj2fY+wEuBqwK/Ap6bmbuH7NeeRknaBHsa1X32NEqaL126T+MVQC8zL+pbdjbw\n34HXrdv2R8CDM3NPRNwSOBU4ajppSpIkSdL8aNP01GBdPpn57cz852Zd//KzMnNP8/ybwNUj4qpT\ny1TaFlXpBKSBqtIJSCNVpROQBrKnUV3WpqIxgVMj4qsR8YcbfVEzhfXMzLxscqlJkiRJ0nxq0/TU\n32ymm14f+HREnJeZXxj1gmZq6l8BvzV61ycAu5rnO4DjgF4TV81XY+Npx72W5WNsXK2J2pKPsfEk\n4tWzQr2esfH2xavako/x/MZLS0usrKwAsLy8zLhacyGcfhFxIvDzzHxlE+8G/mz1QjjNsqOAzwLH\nZ+bpI/blhXAkaRO8EI66zwvhSJov414I56DtTGarIuKaEXFI8/xawH2Bc9Zv1rf9tYGPAs8bVTBK\n7VaVTkAaqCqdgDRSVToBaaD1ZxulLmlF0QgsAl+IiDOB04GPZOanIuK3I+JC4E7ARyPiE832TwFu\nBPxFRJwZEWdExPXKpC5JkiRJ3dXK6anbyempkrQ5Tk9V9zk9VdJ86cT0VEmSJElSO1k0SsVUpROQ\nBqpKJyCNVJVOQBrInkZ1mUWjJEmSJGkoexolSWvY06jus6dR0nyxp1GSJEmSNDEWjVIxVekEpIGq\n0glII1WlE5AGsqdRXWbRKEmSJEkayp5GSdIa9jSq++xplDRf7GmUJEmSJE2MRaNUTFU6AWmgqnQC\n0khV6QSkgexpVJdZNEqSJEmShrKnUZK0hj2N6j57GiXNl3F7Ghe2M5n22vL3R5LmlMdNddfi4s7S\nKUjSTJmHMoLCAAAgAElEQVSLotG/JqqNqqqi1+uVTkPaTxWehVF7eexUWzk21WX2NEqSJEmShpqL\nnsauf0ZJ2lYR4HFTkqTO8D6NkiRJkqSJsWiUCvF+TmqrqnQC0ggeO9VWjk11mUWjJEmSJGmouehp\nLJ2DJM2SxBtuaLDFxZ3s2bNcOg1J0iaN29M4J0Vjtz+jJG2nJAiPmxrI27FI0izyQjjSzKpKJyAN\nVJVOQBrBvjG1lWNTXWbRKEmSJEkayumpkqQ1nJ6q4ZyeKkmzyOmpkiRJkqSJaUXRGBFHRcRpEXFu\nRJwdEU9tlp8YEd+NiDOax/37XnNsRHwpIs6JiLMi4mrlPoG0FVXpBKSBqtIJSCPYN6a2cmyqyxZK\nJ9C4HHhWZi5FxCHA1yPi0826V2bmK/s3joirAG8DHpeZ50TE4cBl001ZkiRJkrqvlT2NEfFB4K+B\nuwKXZOYr1q1/APCYzHziBvZlT6MkbYI9jRrOnkZJmkWd62mMiF3AccCXm0VPjoiliHhjRFy7WXaT\nZttPRsTXIuI5089UkiRJkrqvVUVjMzX1fcDTM/MS4G+BG2XmccAeYHWa6gJwF+AxwN2A/x4R9yyQ\nsjSGqnQC0kBV6QSkEewbU1s5NtVlbelpJCIWqAvGt2XmhwAy80d9m7wB+Ejz/LvA5zLzoua1Hwdu\nC+wevPcTgF3N8x3UJzJ7TVw1X42NjY2N12pLPsbtipuo+QW51+tNNS79/sbGw+KlpaVW5WM83/HS\n0hIrKysALC8vM67W9DRGxFuBH2fms/qWHZGZe5rnzwRun5mPjYgdwGeoex4vBz5BfcGcTwzYrz2N\nkrQJ9jRqOHsaJWkWjdvT2IozjRFxF+BxwNkRcSZ1lfcC4LERcRxwBbAM/DFAZq5ExCuBrzXrPjao\nYJQkSZIkjac1ZxonxTONaq+KfdO+pPbYTXBPj5saqPyZxqqqrpyCJbWJY1Nt1rmrp0qSJEmS2sMz\njZKkNexp1HDlzzRKkjbPM42SJEmSpImxaJSKqUonIA1UlU5AGmH10vJS2zg21WUWjZIkSZKkoexp\nlCStYU+jhrOnUZJmkT2NkiRJkqSJsWiUiqlKJyANVJVOQBrBvjG1lWNTXWbRKEmSJEkayp5GSdIa\n9jRqOHsaJWkW2dMoSZIkSZoYi0apmKp0AtJAVekEpBHsG1NbOTbVZRaNkiRJkqSh5qSnUZK0UQls\nuelBnba4uJM9e5ZLpyFJ2qRxexoXtjOZtup6YSxJ2yq82IkkSdrH6alSIfY+qK2q0glII3jsVFs5\nNtVlFo2SJEmSpKHmoqex659RkrZVBHjclCSpM7xPoyRJkiRpYiwapULsfVBbVaUTkEbw2Km2cmyq\ny+bi6qkRXjxekvp56wRJkrRRc9HTWN91TJK0z4jbatjTKElSp9jTKEmSJEmaGItGqZiqdALSQFXp\nBKQR7BtTWzk21WUWjZIkSZKkoexplKS5ZE+jJEnzojM9jRHxpojYGxHf6Fv28og4LyKWIuL9EXFY\ns3whIv5vRHwjIr4ZEc8rl7kkSZIkdVdrikbgZOB+65Z9CrhlZh4H/DPw/Gb5I4GrZeaxwO2AP46I\no6eWqbQtqtIJSANVpROQRrBvTG3l2FSXtaZozMwvABetW/aZzLyiCU8HjlpdBVwrIq4CXBP4JfCz\naeUqSZIkSfOiVT2NEbET+EhzBnH9ug8D78rMUyJiAXgbcG/gGsAzM/ONQ/ZpT6Mk7ceeRkmS5kVn\nehpHiYj/BVyWmac0i+4AXA4cARwDPDsidpXJTpIkSZK6a6F0AgcSEccDDwTu1bf4scAnm6mrP4qI\nL1L3Ni4P3ssJwK7m+Q7gOKDXxFXz1dh42vHq87bkYzxv8Wr/Ta+3Nl7dYth6Y+OS8eqytuRjbLwa\nLy0t8YxnPKM1+RjPd7y0tMTKygoAy8vLjKtt01N3UU9PvXUT3x94BXD3zPxJ33bPBW6amb8fEdcC\nvgI8KjPPGbBPp6eqpSr2/TIvTdvw6alVBL0W/d8g9auq6spfjKQ2cWyqzcadntqaojEiTqH+Dfq6\nwF7gROAFwNWA1YLx9Mz806ZQPBm4RbP8zZn5yiH7tWiUpP3Y0yhJ0rzoTNE4KRaNkjSIRaMkSfNi\nLi6EI3VTVToBaaCqdALSCKu9O1LbODbVZRaNkiRJkqShnJ4qSXPJ6amSJM0Lp6dKkiRJkibGolEq\npiqdgDRQVToBaQT7xtRWjk11mUWjJEmSJGkoexolaS7Z0yhJ0rywp1GSJEmSNDEWjVIxVekEpIGq\n0glII9g3prZybKrLLBolSZIkSUMdsKcxIh6fmW+fUj7bzp5GSRrEnkZJkubFuD2NCxvY5pkR8Uvg\nZ8DXM/PHW30zSZIkSdJs2cj01Kdm5nuBLwK3jIhHRsSjIuIpEfGbE85P6rCqdALSQFXpBKQR7BtT\nWzk21WUHPNOYmV9qvl4SEXuA+wCPBM4Fzp9settly2diJamTFhd3lk5BkiTNiI30NB4BPBp4PHVz\n4FuAd2bmTyaf3vgiIrt+L0pJ2lb2NEqS1CnT6Gm8EPgQcHxmfnOrbyRJkiRJmj0b6Wl8DvA3wO0j\n4okR8fiIuEdEXCsiHjHh/KTOsvdBbVWVTkAawWOn2sqxqS7bSE/jq9cva6as3hN4HvC+CeQlSZIk\nSWqBA/Y0jnxxxD0zc/c25rPt7GmUpE2yp1GSpE4Zt6dxrKJxFlg0StImWTRKktQp07gQzsyL8JYb\nkrbf4uJO9uxZLp3GtquAXuEcpGGqqqLX65VOQ9qPY1NdNhdFY32nEKltKvzVfLbt3esfpCRJUvdt\neHpqRFwd+B1gF33FZma+cCKZbZOISItGSZMRdHL6u9NTJUnqlGlOT/0QcDHwdeCXW31DSZIkSdLs\n2Mh9GlcdlZmPysyXZ+YrVh8Ty0zqvKp0AtJAVekEpBG8F57ayrGpLttM0filiLj1xDKRJEmSJLXO\nAXsaI+Js6qbABeDGwHeop6cGkJl57MSTjLg28EbgVsAVwJMy88vNumcDLweul5k/HfBaexolTYg9\njZIkqf2m0dP44K3ufBu9Bvh4Zj4yIhaAawJExFHAfYDzSyYnSZIkSV11wOmpmXl+Zp4P/Onq8/5l\nk04wIg4F7paZJzf5XJ6ZP2tWvwp4zqRzkCajKp2ANFBVOgFpBPvG1FaOTXXZZnoaf2vAsgdsVyIj\nHAP8OCJOjogzIuL1EXHNiHgIcGFmnj2FHCRJkiRpLh1wempE/A/qM4o3iohv9K06FPjSpBLrswDc\nFnhyZn4tIl4FnATcnbWF7Ig5uidQ314SYAdwHPtuql41X42Npx33WpaP8ebj+i/LvV7vyufA7Md9\nn60V+RgbGxvPSLyqLfkYz2+8tLTEysoKAMvLy4xrIxfCuTZwOPBXwPP6Vv180IVntltELAL/mJnH\nNPFdqYvGWwGXUheLRwHfA+6QmT9c93ovhCNpQrwQjiRJar9xL4Rz0IE2yMyLM3MZ+BrwO32PEyLi\n9yPiuK2++UZk5l7gwoi4SbPo3sDXM/OIzDwmM/8L8F3gNusLRqndqtIJSANVpROQRlh/RkdqC8em\numwjV09d9d+A2wEfaeIHA98A/iQi3puZL9/u5Po8DXhHRFyV+pYfv7dufTJyeqokSZIkaSsOOD31\nyg0jPg88MDMvaeJDgI8B96c+83eLiWU5BqenSpocp6dKkqT2m/j01D43AH7VF18GLGbmL4BfbjUB\nSZIkSVJ7baZofAdwekScGBEnAl8ETomIawHnTiQ7qdOq0glIA1WlE5BGsG9MbeXYVJdtuKcxM18U\nER8H7krdP/gnmfm1ZvXjJpGcJEmSJKmsDfc0zip7GiVNjj2NkiSp/cbtadzwmcaIuDr1rTZ29b8u\nM1+41TeXJEmSJLXbZnoaPwQ8DLgc+Pe+h6QtqUonIA1UlU5AGsG+MbWVY1Ndtpn7NB6VmfefWCaS\nJEmSpNbZzH0aXw/8dWaePdmUtpc9jZImx55GSZLUfuP2NG6maDwXuDHwHer7MgaQmXnsVt98Giwa\nJU2ORaMkSWq/cYvGzfQ0PgD4r8B9gYcAD26+StqSqnQC0kBV6QSkEewbU1s5NtVlmykaLwDuBhyf\nmedTn75bnEhWkiRJkqRW2Mz01L8DrgDulZk3j4jDgU9l5u0nmeC4nJ4qaXKcnipJktpvavdpBO6Y\nmbeNiDMBMvOiiLjaVt9YkiRJktR+m5meellEXIXmtF1EXJ/6zOMMCB8+fPjY9sfi4k66qCqdgDSC\nfWNqK8emumwzZxpfC/w9cIOIeDHwCODPJ5LVNuvk9DHNvKqq6PV6pdOQJEmSRtpwTyNARNwMuDf1\nn9k/C9wvM189ody2RUSkRaMkbYI9jZIkdcrU7tM45M0vyMyjt7yDKbBolKRNsmiUJKlTpnmfxoHv\nP+brpbll74PaqiqdgDSCx061lWNTXTZu0eifoiVJkiSpww44PTUifs7g4jCAa2TmZi6mM3VOT5Wk\nTXJ6qiRJnVK0p3EWRES3P6A0JYuLO9mzZ7l0GpoGi0ZJkjqldE/jjEgfPlr42N2CHDb+2Lv3fDQf\nqtIJSCPYN6a2cmyqy+akaJQkSZIkbcWcTE/t9meUpiPo+vFCDaenSpLUKU5PlSRJkiRNjEWjVExV\nOgFpoKp0AtII9o2prRyb6rLWF40RcfWI+HJEnBkRZ0fEic3yt0fEtyLiGxHxxoi4SulcJUmSJKlr\nZqKnMSKumZmXNoXhF4GnAdfJzE82608BPpeZrxvwWnsapW1hT+PcsKdRkqROmYuexsy8tHl6dWCh\nXlQXjI2vAEdNPTFJkiRJ6riZKBoj4qCIOBPYA3w6M7/at24BeALwyWGvl9qpKp2ANFBVOgFpBPvG\n1FaOTXXZQukENiIzrwBuExGHAR+MiFtk5rnN6r+lnpr6xeF7OAHY1TzfARwH9Jq4ar4aGxsfKF79\nD7HXM+5yvKot+Rgb98er2pKPsfFqvLS01Kp8jOc7XlpaYmVlBYDl5WXGNRM9jf0i4i+BSzLzlc1F\ncX4jMx8+Ynt7GqVtYU/j3LCnUZKkTul8T2NEXC8irt08vwZwH+BbEfEHwH2Bx5TMT5IkSZK6rPVF\nI/BrwO6IWAK+DJyamR8H/g64AXB6RJwREX9eMklp86rSCUgDVaUTkEZYnYYltY1jU13W+p7GzDwb\nuO2A5VctkI4kSZIkzZWZ62ncLHsape1iT+PcsKdRkqRO6XxPoyRJkiSpHItGqZiqdALSQFXpBKQR\n7BtTWzk21WUWjZIkSZKkoexplLRB9jTODXsaJUnqFHsaJUmSJEkTY9EoFVOVTkAaqCqdgDSCfWNq\nK8emusyiUZIkSZI0lD2NkjbInsa5YU+jJEmdYk+jJEmSJGliLBqlYqrSCUgDVaUTkEawb0xt5dhU\nl1k0SpIkSZKGsqdR0gbZ0zg37GmUJKlTxu1pXNjOZNpry98fSY3FxZ2lU5AkSVIBc1E0enZEbVRV\nFb1er3Qa0n4qoFc4B2kYj51qK8emusyeRkmSJEnSUHPR09j1zyhJ28qeRkmSOsX7NEqSJEmSJsai\nUSrE+zmprarSCUgjeOxUWzk21WUWjZIkSZKkoeaip7F0DtK0LC7uZM+e5dJpaNbZ0yhJUqeM29M4\nJ0Vjtz+jtE94ixmNz6JRkqRO8UI40syqSicgDVSVTkAawb4xtZVjU11m0ShJkiRJGsrpqVKnOD1V\n28DpqZIkdYrTUyVJkiRJE9P6ojEi3hQReyPiG+uWPzUivhURZ0fES0vlJ21dVToBaaCqdALSCPaN\nqa0cm+qyhdIJbMDJwF8Db11dEBE94CHArTLz8oi4XqHcJEmSJKnTZqKnMSJ2Ah/JzGOb+N3A6zLz\ntA281p5GzRF7GrUN7GmUJKlT5rWn8SbA3SPi9IjYHRG3K52QJEmSJHXRLExPHWQB2JGZd4qI2wPv\nAY4ZvvkJwK7m+Q7gOKDXxFXz1dh42vHq8+3cf91T0ev1rnwOGBtvKoZmhLYkH2Pj/nh1WVvyMTZe\njZeWlnjGM57RmnyM5zteWlpiZWUFgOXlZcY1q9NTPw68NDM/38T/AtwxM38y4LVOT1VLVdS/mm8n\np6dqfFUEPceRWqqqqit/MZLaxLGpNht3euqsFI27qIvGWzfxHwFHZuaJEXET4NOZuXPIay0aNUcs\nGrUN7GmUJKlTxi0aWz89NSJOoT4dc92IuAA4EXgzcHJEnA38EnhiuQwlSZIkqbtm4kzjODzTqPaq\ncHqq2sjpqWozpwCqrRybarN5vXqqJEmSJGkKPNModYpnGrUN7GmUJKlTPNMoSZIkSZoYi0apmKp0\nAtJAVekEpBFW70cmtY1jU11m0ShJkiRJGsqeRqlT7GnUNrCnUZKkTrGnUZIkSZI0MRaNUjFV6QSk\ngarSCUgj2DemtnJsqsssGiVJkiRJQ9nTKHWKPY3aBvY0SpLUKfY0SpIkSZImxqJRKqYqnYA0UFU6\nAWkE+8bUVo5NdZlFoyRJkiRpqDnpaZTmw+LiTvbsWS6dhmadPY2SJHXKuD2NC9uZTFt1vTCWJEmS\npElxeqpUiL0PaquqdALSCB471VaOTXWZRaMkSZIkaai56Gns+meUpG1lT6MkSZ3ifRolSZIkSRNj\n0SgVYu+D2qoqnYA0gsdOtZVjU102F1dPjdjymVhp5njbDUmSJG2nuehphG5/Rmmt8DYzGo89jZIk\ndYo9jZIkSZKkibFolIqpSicgDVSVTkAawb4xtZVjU11m0ShJkiRJGsqeRqlz7GnUmOxplCSpU+a6\npzEinhkR50TENyLiHRFxtdI5SZIkSVKXzGzRGBE3BJ4K3DYzj6W+fcijy2YlbUZVOgFpoKp0AtII\n9o2prRyb6rJZv0/jVYBrRcQVwDWB7xfOR5IkSZI6ZaZ7GiPiacCLgUuBT2XmEwZsY0+j5ow9jRqT\nPY2SJHXKuD2NM3umMSJ2AA8DdgIXA++LiMdm5in7b30CsKt5vgM4Dug1cdV8NTbuUtxEzVSZXq9n\nbLzxGNqVj7GxsbGxsfGm4qWlJVZWVgBYXl5mXDN7pjEiHgHcLzP/sImfANwxM5+ybjvPNKqlKvb9\ner6dPNOo8VQR9BxDaqmqqq78xUhqE8em2myer556AXCniDg4IgK4N3Be4ZwkSZIkqVNm9kwjQESc\nSH3F1MuAM4E/yMzL1m3jmUbNGc80akz2NEqS1Cnjnmmc6aJxIywaNX8sGjUmi0ZJkjplnqenSjOu\nKp2ANFBVOgFphNULPkht49hUl1k0SpIkSZKGcnqq1DlOT9WYnJ4qSVKnOD1VkiRJkjQxFo1SMVXp\nBKSBqtIJSCPYN6a2cmyqyywaJUmSJElD2dModY49jRqTPY2SJHWKPY2SJEmSpImxaJSKqUonIA1U\nlU5AGsG+MbWVY1NdZtEoSZIkSRrKnkapc+xp1JjsaZQkqVPsaZQkSZIkTYxFo1RMVToBaaCqdALS\nCPaNqa0cm+qyhdIJTMeWz8RKM2dxcWfpFCRJktQhc9HT2PXPKEnbyp5GSZI6xZ5GSZIkSdLEWDRK\nhdj7oLaqSicgjeCxU23l2FSXWTRKkiRJkoayp1GStJY9jZIkdYo9jZIkSZKkiZmLojEifPiYuccR\nR+wq/U9Hc6oqnYA0gn1jaivHprpsTu7T6DQrtVEF9Iau3bvX+4tKkiSpvLnoabRo1GwKuv7vUy0V\n9jRKktQlEfY0SpIkSZImxKJRKqYqnYA0UFU6AWkE+8bUVo5NdZlFoyRJkiRpqJnvaYyIg4CvAd/N\nzIcOWG9Po2aUPY0qxJ5GSZI6xZ5GeDpwbukkJEmSJKmLZrpojIijgAcCbyydi7R5VekEpIGq0glI\nI9g3prZybKrLZrpoBF4FPAfnn0qSJEnSRCyUTmCrIuJBwN7MXIqIHjBiju4JwK7m+Q7gOPbdVL1q\nvhobTzvuHXD71b9a9nrGxlOMoV35GBsbG89IvKot+RjPb7y0tMTKygoAy8vLjGtmL4QTES8BHg9c\nDlwDOBT4QGY+cd12XghHM8oL4agQL4QjSVKnzO2FcDLzBZl5dGYeAzwaOG19wSi1W1U6AWmgqnQC\n0gjrz+hIbeHYVJfNbNEoSZIkSZq8mZ2eulFOT9XscnqqCnF6qiRJnTK301MlSZIkSZNn0SgVU5VO\nQBqoKp2ANIJ9Y2orx6a6zKJRkiRJkjSUPY1Sa9nTqELsaZQkqVPsaZQkSZIkTYxFo1RMVToBaaCq\ndALSCPaNqa0cm+oyi0ZJkiRJ0lD2NEqtZU+jCrGnUZKkTrGnUZIkSZI0MRaNUjFV6QSkgarSCUgj\n2DemtnJsqsssGiVJkiRJQ9nTKLWWPY0qxJ5GSZI6ZdyexoXtTKa9tvz9kYpZXNxZOgVJkiRpPqan\nZqYPH6177N69e+T6PXuWS//T0ZyqSicgjWDfmNrKsakum4uiUZIkSZK0NXPR09j1zyhJ28qeRkmS\nOsX7NEqSJEmSJsaiUSrE3ge1VVU6AWkEj51qK8emusyiUZIkSZI01Fz0NJbOQd2yuLjTK5uq2+xp\nlCSpU8btaZyTorHbn1HTFnT9343mnEWjJEmd4oVwpBll74PaqiqdgDSCx061lWNTXWbRKEmSJEka\nyump0qY5PVUd5/RUSZI6xempkiRJkqSJsWiUCrH3QW1VlU5AGsFjp9rKsakum9miMSKOiojTIuLc\niDg7Ip5WOidJkiRJ6pqZ7WmMiCOAIzJzKSIOAb4OPCwzv7VuO3satc3saVTH2dMoSVKnzG1PY2bu\nycyl5vklwHnAkWWzkiRJkqRumdmisV9E7AKOA75cNhNp4+x9UFtVpROQRvDYqbZybKrLFkonMK5m\naur7gKc3ZxwHOAHY1TzfQV1f9pq4ar4aG28mbqLmP4her2ds3Jl4VVvyMTbuj1e1JR9j49V4aWmp\nVfkYz3e8tLTEysoKAMvLy4xrZnsaASJiAfgo8InMfM2Qbexp1Dazp1EdZ0+jJEmdMrc9jY03A+cO\nKxglSZIkSeOZ2aIxIu4CPA64V0ScGRFn/P/t3W3IZOdZB/D/ZZYUFbRYMFRjt0qthWINgi8fIlmt\n0oLYCJVSpa2tolZpISooiOAHP6gIgmliqrbEUExDfANjKxS0k9oWIZJsq9VWjay21eSLWUQrsS63\nH5554rPJc2b3eTnPfc+Z3w+WnXNmzplrkmvvnWvn/J+pqlf3rguu1/6lBDCaVe8CYANrJ6PSmyzZ\n1mYaW2sfSXJD7zoAAACWbKszjddDppHTJ9PIwsk0AsCi7HqmEQAAgBkZGqET2QdGtepdAGxg7WRU\nepMlMzQCAAAwSaYRjkymkYWTaQSARZFpBAAAYDaGRuhE9oFRrXoXABtYOxmV3mTJDI0AAABMkmmE\nI5NpZOFkGgFgUWQaAQAAmI2hETqRfWBUq94FwAbWTkalN1kyQyMAAACTdiTTCKfnppvO54knLvUu\nA+Yj0wgAi3LSTOO50yxmVEsfjAEAAObi8lToRPaBUa16FwAbWDsZld5kyQyNAAAATNqJTOPSXyPA\nqZJpBIBF8T2NAAAAzMbQCJ3IPjCqVe8CYANrJ6PSmyzZTvz01KpjfxILx+arOQAAWIKdyDQmy36N\njKp83QvbSaYRABZFphEAAIDZGBqhE9kHRrXqXQBsYO1kVHqTJTM0AgAAMEmmEWYj08iWkmkEgEXZ\n2UxjVb27qp6sqo/3rgUAAGCptnZoTHJvklf1LgKOS/aBUa16FwAbWDsZld5kybZ2aGytfTjJU73r\nAAAAWLKtzjRW1fkkD7XWXrHhMTKNdCLTyJaSaQSARdnZTCMAAADzO9e7gLPx5iQvXt9+fpJbklxY\nb6/Wv9u2ffrb+/mGCxeeu30w+3DY/bZt99pO9jp4lHps2z64vb9vlHps297fvnjxYu64445h6rG9\n29sXL17M5cuXkySXLl3KSW375akvzt7lqV+/4TEuT6WTzZenrlarZ/5ww0hWVbmwxX83sGzWTkal\nNxnZSS9P3dqhsaruz94/hr8gyZNJfqG1du8hjzM00olMI1tKphEAFmVnh8brZWikH0MjW8rQCACL\n4gfhwJbav/4cRrPqXQBsYO1kVHqTJTM0AgAAMMnlqTAbl6eypVyeCgCL4vJUAAAAZmNohE5kHxjV\nqncBsIG1k1HpTZbM0AgAAMAkmUaYjUwjW0qmEQAWRaYRAACA2RgaoRPZB0a16l0AbGDtZFR6kyUz\nNAIAADBJphFmI9PIlpJpBIBFkWkEAABgNoZG6ET2gVGtehcAG1g7GZXeZMnO9S7gbBz7k1g4tptu\nOt+7BAAAOLGdyDQu/TUCnCqZRgBYFJlGAAAAZmNohE5kHxjVqncBsIG1k1HpTZbM0AgAAMAkmUYA\nribTCACLItMIAADAbAyN0InsA6Na9S4ANrB2Miq9yZIZGgEAAJgk0wjA1WQaAWBRZBoBAACYjaER\nOpF9YFSr3gXABtZORqU3WTJDIwAAAJNkGgG4mkwjACyKTCMAAACzMTRCJ7IPjGrVuwDYwNrJqPQm\nS2ZoBAAAYJJMIwBXk2kEgEWRaQQAAGA2hkboRPaBUa0++MHeJcAkayej0pssmaERAACASTKNAAAA\nCybTCAAAwGwMjdCJ7AOj0puMTH8yKr3JkhkaAQAAmCTTCAAAsGAyjQAAAMzG0AidyD4wKr3JyPQn\no9KbLJmhEQAAgEkyjQAAAAsm0wgAAMBsDI3QiewDo9KbjEx/Miq9yZIZGgEAAJgk0wgAALBgMo0A\nAADMxtAIncg+MCq9ycj0J6PSmyyZoREAAIBJMo0AAAALJtMIAADAbAyN0InsA6PSm4xMfzIqvcmS\nGRoBAACYJNMIAACwYDKNAAAAzMbQCJ3IPjAqvcnI9Cej0pssmaERAACASTKNAAAACybTCAAAwGwM\njdCJ7AOj0puMTH8yKr3JkhkaoZOLFy/2LgEOpTcZmf5kVHqTJTM0QieXL1/uXQIcSm8yMv3JqPQm\nSxDsny4AAAYESURBVGZo3BI9LnmY4zlP45zHOcdRjrnex17rcbtymUqv1zlif25Lbx71ebeZtfNk\nx1s752PtPNnx1s55WTuPf/xRjzmtvpv7/5mhcUv4w3uyc4z4xufSpUvX9Tyj88bnZMeP+MZnKb2Z\nWDtPery1cz7WzpMdb+2cl7Xz+McvdWjcia/c6F0DAABATyf5yo3FD40AAAAcn8tTAQAAmGRoBAAA\nYJKhEQAAgEmGRgAAACYZGgEAAJi0k0NjVd1eVb9VVe+tqu/qXQ/sq6qvrqp3VdWDvWuBg6rqi6rq\nd6rqN6vqB3rXA/usm4zMe05GVVUvq6p7qurBqnrrNR+/y1+5UVXPT/KrrbUf6V0LHFRVD7bWXte7\nDthXVW9I8lRr7X1V9UBr7fW9a4KDrJuMzHtORlVVleS+1tqbNj1uqz9prKp3V9WTVfXxZ+1/dVV9\nsqr+vqp+dsMpfj7J3fNWyS46hd6EWR2jR29O8un17StnVig7x/rJyE7Qn95zMqvj9GZVfU+SP0ny\n/mudf6uHxiT3JnnVwR1V9QVJ7lrvf3mS76+ql63ve2NV/VpVfUVV/XKS97fWLp510eyE4/bmC/cf\nfpbFspOO1KPZGxhv3n/oWRXJTjpqbz7zsLMpjx135P70npMzcuTebK091Fr77iRvuNbJt3pobK19\nOMlTz9r9zUn+obX2z621zyd5IMnt68e/p7X2U0lem+SVSb6vqn70LGtmN5ygN5+uqnuS3OJf0pnT\nUXs0yR9lb828O8lDZ1cpu+aovVlVX2bd5Kwcoz/fHu85OQPH6M3bqurXq+qdSd53rfOfO+2CB/CV\n+f9LqJLkM9n7D/aM1to7krzjLIuCXF9v/nuSHz/LouCAyR5trX0uyQ/1KAqyuTetm/S2qT+956Sn\nTb35cJKHr/dEW/1J44TDLk/Z3Z/2w0j0JqPTo4xKbzIy/cmoTq03lzg0fibJiw5s35zkXzvVAgfp\nTUanRxmV3mRk+pNRnVpvLmForFw9RT+S5CVVdb6qbkzy+iR/3KUydp3eZHR6lFHpTUamPxnVbL25\n1UNjVd2f5KNJXlpV/1JVb2mtXUny9iQfSPKJJA+01v6uZ53sHr3J6PQoo9KbjEx/Mqq5e7Nac8k1\nAAAAh9vqTxoBAACYl6ERAACASYZGAAAAJhkaAQAAmGRoBAAAYJKhEQAAgEmGRgAAACYZGgHYOVX1\n5VX1u1X1j1X1SFV9pKpuv8YxL6yqB0/p+W+tqr+pqker6nnPuu/Kev9j699fVFW3VdXlA/s/UFU/\nt779WFX97/q+R6vqbadRIwDsq9Za7xoA4ExV1UeT3Nta++319lcleU1r7e4zev57kvxFa+3+Q+77\nj9balzxr321Jfrq19pqJ8z3nGAA4LT5pBGCnVNV3JHl6f2BMktbap/cHxqo6X1Ufqqq/Wv/61gP7\n/3p9+wer6g+q6k+r6lNV9SsTz/XK9ad/H6uqd1XVjVX1w0lel+QXq+o9hx02VfoJXjYAHNu53gUA\nwBl7eZJHN9z/ZJLvbK39T1W9JMl7k3zT+r6Dl+d8Q5Jbknw+yaeq6s7W2mf371xfdnpvkm9vrT1e\nVfcleWtr7c6qujXJQ621Pzzk+b+wqh7N3pD4T6211673f9t6f5L8Xmvtl470qgHgmAyNAOy0qror\nya3Z+/TxW5LcmOSuqrolyZUkXztx6J+11v5zfY6/TXI+yWcP3P912Rv6Hl9v35fkJ5LceY2SPtda\n+8ZD9n9o6vJUAJiToRGAXfOJJPuf3qW19raqekGSR9a7fjLJE621V1TVDUn+e+I8Tx+4fSXP/Tu1\n4pJSABZAphGAndJa+/Mkz6uqHzuw+4sP3P7SJP+2vv2mJDcc86k+meR8VX3NevuNSR6+juOOM2ga\nTgGYjaERgF30vUkuVNXjVfWX2cse/sz6vt9I8uaqeizJS5P813Wc7zk/iry19nSStyT5/ar6WPY+\njXzn1OM3nes4zw8Ap8VXbgAAADDJJ40AAABMMjQCAAAwydAIAADAJEMjAAAAkwyNAAAATDI0AgAA\nMMnQCAAAwKT/AzqNsKd+31GvAAAAAElFTkSuQmCC\n",
"text/plain": [
"<matplotlib.figure.Figure at 0x111552f98>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"%matplotlib inline\n",
"import matplotlib.pyplot as plt\n",
"import numpy as np\n",
"import timeit\n",
"\n",
"n = np.arange(14) # lengths = 2**n to evaluate\n",
"reps = 100 # number of repetitions per measurement\n",
"\n",
"# measure execution times\n",
"gain = np.zeros(len(n))\n",
"for N in n:\n",
" length = 2**N\n",
" # setup environment for timeit\n",
" tsetup = 'import numpy as np; from scipy.linalg import dft; \\\n",
" x=np.random.randn(%d)+1j*np.random.randn(%d); F = dft(%d)' % (length, length, length)\n",
" # DFT\n",
" tc = timeit.timeit('np.matmul(F, x)', setup=tsetup, number=reps)\n",
" # FFT\n",
" tf = timeit.timeit('np.fft.fft(x)', setup=tsetup, number=reps)\n",
" # gain by using the FFT\n",
" gain[N] = tc/tf\n",
"\n",
"# show the results\n",
"plt.figure(figsize = (15, 10))\n",
"plt.barh(n-.5, gain, log=True)\n",
"plt.plot([1, 1], [-1, n[-1]+1], 'r-')\n",
"plt.yticks(n, 2**n)\n",
"plt.xlabel('Gain of FFT')\n",
"plt.ylabel('Length $N$')\n",
"plt.title('Ratio of execution times between DFT and FFT')\n",
"plt.grid()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Exercise**\n",
"\n",
"* For which lengths $N$ is the FFT algorithm faster than the direct computation of the DFT? \n",
"* Why is it slower below a given length $N$?\n",
"* Does the trend of the gain follow the expected numerical complexity of the radix-2 FFT algorithm?"
]
},
{
"cell_type": "markdown",
"metadata": {
"nbsphinx": "hidden"
},
"source": [
"**Copyright**\n",
"\n",
"The notebooks are provided as [Open Educational Resource](https://de.wikipedia.org/wiki/Open_Educational_Resources). Feel free to use the notebooks for your own educational purposes. The text is licensed under [Creative Commons Attribution 4.0](https://creativecommons.org/licenses/by/4.0/), the code of the IPython examples under the [MIT license](https://opensource.org/licenses/MIT). Please attribute the work as follows: *Lecture Notes on Signals and Systems* by Sascha Spors."
]
}
],
"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.1"
},
"nbsphinx": {
"execute": "never"
}
},
"nbformat": 4,
"nbformat_minor": 0
}