{
"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 as representative."
]
},
{
"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 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 powers of two. Similar and other principles can be applied to derive efficient algorithms for 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 take some time."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"data": {
"application/pdf": "JVBERi0xLjQKJazcIKu6CjEgMCBvYmoKPDwgL1BhZ2VzIDIgMCBSIC9UeXBlIC9DYXRhbG9nID4+CmVuZG9iago4IDAgb2JqCjw8IC9FeHRHU3RhdGUgNCAwIFIgL0ZvbnQgMyAwIFIgL1BhdHRlcm4gNSAwIFIKL1Byb2NTZXQgWyAvUERGIC9UZXh0IC9JbWFnZUIgL0ltYWdlQyAvSW1hZ2VJIF0gL1NoYWRpbmcgNiAwIFIKL1hPYmplY3QgNyAwIFIgPj4KZW5kb2JqCjEwIDAgb2JqCjw8IC9Bbm5vdHMgWyBdIC9Db250ZW50cyA5IDAgUgovR3JvdXAgPDwgL0NTIC9EZXZpY2VSR0IgL1MgL1RyYW5zcGFyZW5jeSAvVHlwZSAvR3JvdXAgPj4KL01lZGlhQm94IFsgMCAwIDQyMi40NTc1IDMyNy4zMDg4NzUgXSAvUGFyZW50IDIgMCBSIC9SZXNvdXJjZXMgOCAwIFIKL1R5cGUgL1BhZ2UgPj4KZW5kb2JqCjkgMCBvYmoKPDwgL0ZpbHRlciAvRmxhdGVEZWNvZGUgL0xlbmd0aCAxMSAwIFIgPj4Kc3RyZWFtCniczVpNb9w2EL3rV/CYHkKTw+9jjNQGiqJAUgM9FD2kzsZxIDtI4jb9+X2UtKsZRrtJhNa7ARx7H2dGmkfqcYZaq951Vt0oo97h57P6Xf2B36+VVZf4uekMPt11nkj7kAI+9PMHR0k7kzP+7GEoPr7tujfd2TOE+ASny64LXkdXnbLOyVYjhLVWpyDAnoHG65LdiO7cOThc5IP6MrQLSVtSFKO2+P1xo35T98oAsyHmSDUT7eP0L3VGJxNyJheoqI83+y1VY9l1T7OlDNQ5cl5Fo6P3ccguGdyOIUoM7RmanA4+mQGVUfjIkOMLdVJZpqhTojFLa4suNubC0J6huWhAS1nykVPMspAuxsxZJsybDQzuOWxN0Nm4pUTF0Clmak3RhvI21aytS5GjPUPJ6uLsYp586CTzpKDhyPJMJta5m/Fe4LhnE8pirnzoJHP1VlOepMgSFAXWkcM9h3HnuN/FTPnQSWaKC3u7lSNPOkF5XOZ4L3BIsLOLkiSGTjLXYWvYilKwdSPOQ047vBd4jtr7ZVXiQyeZa90g4laVMulcfHSF4z3HyZAOaVGZxNAp5krYIlKZlIms08Yn6zjccxgSFcuiLomhk8y0bhI06RKhZsC8OJc43gvceZ1pUZnE0EnmOmwSkzI5TF82FkLE8V7gwegSFpVJDJ1krsM2MSmTw+UxocYFjvcCj0mbvKhMYugkc63bhJmUyZWgA/mUHMd7gWenyS4rEx/67lzPnlFtsEj9hGYNDZzOcwuHdGMxtkSTskc3d/Z88/ft9ebl5bm6/jRUrjF750RDxlDeZ3W/Dnc1Xc2gXfziarXJqyNfDxy027Z63Tl6z8/dB/xv1FODqGBitDMOGgZRwKJX13fd+VXtRHVJMUQK6up1d3YBa6Ou3nRP7A/q6l0XtYuE2441tqkmT8w04EM2YAFlnsKzRsMfU4g0RBgNf7zqXjwS/WQslrhPPguaOLx2Ar4h9OEpQIMzGKaCJhRaJafARJQZ/8MU2EefAlQOPgcbo+SJwaun4OuhD04BoZobDEOAzaM9BfTYU+AibiSHkKVYcHjtFHxD6INT4JD9aGhQYjzaFLh5CsAcbilSPWDI8lmewTifmyGHVKfHpkzCmqHSPCdNxeZGKBgqzbGfYSMHg8KcodLcQnBAkVhAdxxtzbHpJ2yYsbGf4cahnhOWkG1zAQa3DgGTWGKgxmGGGwdfS5dMRRLE4cYhZGz51MTfgY1xrJU6UWzun8GNQy2Kc10+0oHBjQNqC0PRoqARDgxuHbL2KcToG4cZbhyK03i4yTUpM7h1SHg6jImlcZhh6UCEKioDkSxxuHFwRVf6rEyaw40DKmtTKDZKweHWIUGTUNE0GyyDG4fotceE5iYHBrcOBZ+jaWaaw41D8hpLuFlKDG3Ns0YcnxuOGCwdnInIrdimquNw44BuDvVtLfSFA4MbB2e1L04mPIOtcdEuWyzhxnyGGwcfhk7RyHXK4cYhWI0m1Be5TjncOkBFSjIxNw4z3DhEtB/RFOObDWqGGwc8hS6jB5T2MyrNvcEmgQ6jUTsON48mdpvzKzVvZ2SGLVAVpLzdqZ5cvrq9V+/fqIuLq7pVqbpVPUqtsL1C7Q0h70m8P5rB76sSloJiiq1rYy7QUw9SR15QUxBWRKXHHoeU2kSSa0iZwXWk8KBbUkTMg6SkAnFLAyl0HFIK1CC0bxpncB0pPOiWFBHzICn1SCCMD5I/Dim12kFn17DC0HW0iLBbXmTUg8TUqq1utZWZfCRmsNlUi4aZGV3JDA+7Y0ZEXWDGoe0ZsrWUNIqKUVzikZjxSSfXai5DVzLDw+6YEVEPMuOdLm6UXXckibG1nYyt8DJ0JTM87I4ZEfUgM7Ee04zaG4+lM9D/Ulr1ZehKZnjYHTMi6lIlgw4slNEUeZZRgS0dS2lKDfiFBs/oSm542B03IuphbgrV4nzcssORtIYsafKtCjN0HTci7JYbGfUgN2TQAodRhoM9ktoQ7ge/W25mdCU3POyOGxF1iRtwByt0mGkqfQ0dSWvI4y5Mq8IMXckLD7vjRUTdy4uvb5Cm6hcN9ZF4qd87c60GM3QlLzzsjhcRdS8vAT2nmwpgU46lMfWbaqHVX4au5IWH3fEiou7lJZEucSp/ocNbXj509Sz6aT2VRi1NqHmyLjQo1nwWvfQ+5ufhHDnUl+iUy2AwnkRvppNoG8j5eoMjfL89uXZJmN9scT8deI/wwwCj/oIulszivN0TR00OFvduTdk5nF3Q9pZ/+U/ePlgd+LdYh++vvtz/8rO+VY/TtjCjhJ0w58Xtpl5nXAZm7/ri5y2L31tFyIWvvt7t+eorrL/5y7Pcdo5xILKp+UzL24rFfbNbp1TnB5SkQONEWLONUef25auH2/f1TGjzz+b6L/x9rx5u7zaf1J+bh8+bzb16fnGlXt2/FmdG3b+4wSd4CmVuZHN0cmVhbQplbmRvYmoKMTEgMCBvYmoKMTgzMgplbmRvYmoKMTYgMCBvYmoKPDwgL0ZpbHRlciAvRmxhdGVEZWNvZGUgL0xlbmd0aCA4MCA+PgpzdHJlYW0KeJw1jLsNwEAIQ3umYAT+hH2iVHf7t0EkV2A/2cJpFxJKtqRwX+HNwFGDG5RkaIGTIpM1hfifnZdus3pigVmMiznGNCz6JTl2phc8L8EIF+IKZW5kc3RyZWFtCmVuZG9iagoxNCAwIG9iago8PCAvQmFzZUZvbnQgL0RlamFWdVNhbnMtT2JsaXF1ZSAvQ2hhclByb2NzIDE1IDAgUgovRW5jb2RpbmcgPDwgL0RpZmZlcmVuY2VzIFsgNzggL04gXSAvVHlwZSAvRW5jb2RpbmcgPj4gL0ZpcnN0Q2hhciAwCi9Gb250QkJveCBbIC0xMDE2IC0zNTEgMTY2MCAxMDY4IF0gL0ZvbnREZXNjcmlwdG9yIDEzIDAgUgovRm9udE1hdHJpeCBbIDAuMDAxIDAgMCAwLjAwMSAwIDAgXSAvTGFzdENoYXIgMjU1IC9OYW1lIC9EZWphVnVTYW5zLU9ibGlxdWUKL1N1YnR5cGUgL1R5cGUzIC9UeXBlIC9Gb250IC9XaWR0aHMgMTIgMCBSID4+CmVuZG9iagoxMyAwIG9iago8PCAvQXNjZW50IDkyOSAvQ2FwSGVpZ2h0IDAgL0Rlc2NlbnQgLTIzNiAvRmxhZ3MgOTYKL0ZvbnRCQm94IFsgLTEwMTYgLTM1MSAxNjYwIDEwNjggXSAvRm9udE5hbWUgL0RlamFWdVNhbnMtT2JsaXF1ZQovSXRhbGljQW5nbGUgMCAvTWF4V2lkdGggMTM1MCAvU3RlbVYgMCAvVHlwZSAvRm9udERlc2NyaXB0b3IgL1hIZWlnaHQgMCA+PgplbmRvYmoKMTIgMCBvYmoKWyA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMAo2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDMxOCA0MDEgNDYwIDgzOCA2MzYKOTUwIDc4MCAyNzUgMzkwIDM5MCA1MDAgODM4IDMxOCAzNjEgMzE4IDMzNyA2MzYgNjM2IDYzNiA2MzYgNjM2IDYzNiA2MzYgNjM2CjYzNiA2MzYgMzM3IDMzNyA4MzggODM4IDgzOCA1MzEgMTAwMCA2ODQgNjg2IDY5OCA3NzAgNjMyIDU3NSA3NzUgNzUyIDI5NQoyOTUgNjU2IDU1NyA4NjMgNzQ4IDc4NyA2MDMgNzg3IDY5NSA2MzUgNjExIDczMiA2ODQgOTg5IDY4NSA2MTEgNjg1IDM5MCAzMzcKMzkwIDgzOCA1MDAgNTAwIDYxMyA2MzUgNTUwIDYzNSA2MTUgMzUyIDYzNSA2MzQgMjc4IDI3OCA1NzkgMjc4IDk3NCA2MzQgNjEyCjYzNSA2MzUgNDExIDUyMSAzOTIgNjM0IDU5MiA4MTggNTkyIDU5MiA1MjUgNjM2IDMzNyA2MzYgODM4IDYwMCA2MzYgNjAwIDMxOAozNTIgNTE4IDEwMDAgNTAwIDUwMCA1MDAgMTM1MCA2MzUgNDAwIDEwNzAgNjAwIDY4NSA2MDAgNjAwIDMxOCAzMTggNTE4IDUxOAo1OTAgNTAwIDEwMDAgNTAwIDEwMDAgNTIxIDQwMCAxMDI4IDYwMCA1MjUgNjExIDMxOCA0MDEgNjM2IDYzNiA2MzYgNjM2IDMzNwo1MDAgNTAwIDEwMDAgNDcxIDYxNyA4MzggMzYxIDEwMDAgNTAwIDUwMCA4MzggNDAxIDQwMSA1MDAgNjM2IDYzNiAzMTggNTAwCjQwMSA0NzEgNjE3IDk2OSA5NjkgOTY5IDUzMSA2ODQgNjg0IDY4NCA2ODQgNjg0IDY4NCA5NzQgNjk4IDYzMiA2MzIgNjMyIDYzMgoyOTUgMjk1IDI5NSAyOTUgNzc1IDc0OCA3ODcgNzg3IDc4NyA3ODcgNzg3IDgzOCA3ODcgNzMyIDczMiA3MzIgNzMyIDYxMSA2MDgKNjMwIDYxMyA2MTMgNjEzIDYxMyA2MTMgNjEzIDk5NSA1NTAgNjE1IDYxNSA2MTUgNjE1IDI3OCAyNzggMjc4IDI3OCA2MTIgNjM0CjYxMiA2MTIgNjEyIDYxMiA2MTIgODM4IDYxMiA2MzQgNjM0IDYzNCA2MzQgNTkyIDYzNSA1OTIgXQplbmRvYmoKMTUgMCBvYmoKPDwgL04gMTYgMCBSID4+CmVuZG9iagoyMSAwIG9iago8PCAvRmlsdGVyIC9GbGF0ZURlY29kZSAvTGVuZ3RoIDE2MCA+PgpzdHJlYW0KeJw9kEsSwyAMQ/ecQkfA+H+edLpK7r+tDZ1ssBiE9MB9YiKjFieCr8SHBqXDJPBsFYR7MNkRcoTkBE2GsoMkcQ0NBqXCpmOZ78mmddJKrLzRftl3NGaddIotRYd2If/n9SLco+Aa6xk8D2AxyNpKpeyZMFplpq7yqOi1H9PhPQ9Eq8Xl9Qau8NpHN6koKkvq/kR3NNj+kbf7Ht8fmWU4JAplbmRzdHJlYW0KZW5kb2JqCjIyIDAgb2JqCjw8IC9GaWx0ZXIgL0ZsYXRlRGVjb2RlIC9MZW5ndGggNzQgPj4Kc3RyZWFtCnicMzU3VTBQsLQAEqaG5grmRpYKKYZcQD6IlcsFE8sBs8xMzIAsQ0tklomxIZBlYmGGxDI2sYDKIlgGQBpsTQ7M9ByuNAADcRiTCmVuZHN0cmVhbQplbmRvYmoKMjMgMCBvYmoKPDwgL0ZpbHRlciAvRmxhdGVEZWNvZGUgL0xlbmd0aCAyNDQgPj4Kc3RyZWFtCnicTVFJbsQwDLv7FfzAAJasxXlPip7a/19LOhhMD4YYWeISdycmsvCyhboWOhxfNvJK2Az8HrTmxM+IFf/RNiKtfFBtgUzERJHQRd1o3CPd8CpE+5EKXqneY81H3K00b+nYxf7eB9OaR6qsCvGQY3NkI2ldE0XH99B6zw3RKYME+tyEHBClOXoVkv7aD9e10ezW2syeqA4emRLKJ81qaE6nmCGzoR63qVjJKNyoMiruUxlpPcjbOMsATo4Tymg92bGaiPJTn1xCXkzECbvs7FiITSxsHNJ+VPrE8vOtN+NvprWWQsYFidAUl97PeI/vP91YW7QKZW5kc3RyZWFtCmVuZG9iagoyNCAwIG9iago8PCAvRmlsdGVyIC9GbGF0ZURlY29kZSAvTGVuZ3RoIDU5ID4+CnN0cmVhbQp4nDM1NVcwULC0ABKmpkYK5kaWCimGXEA+iJXLZWhpDmblgFkWxkAGSBmcYQCkwZpzYHpyuNIAqeEQWgplbmRzdHJlYW0KZW5kb2JqCjI1IDAgb2JqCjw8IC9GaWx0ZXIgL0ZsYXRlRGVjb2RlIC9MZW5ndGggMjI3ID4+CnN0cmVhbQp4nEWQS44DIRBD95zCR6D+cJ6OsurcfzsuOtFssCUo1zO5AxN78chMlG68ZLg7zBWf4Rkwc/hKmGzETOhOXCOUrhThVJ8IjsvevOmgiXtEzqOeBVnVzg1qAWeS5oLtgi7njBU3zsmtRuXN9KPXEL5pdx/XeYf2SOPew1S+zjnVzruKCGkLWdW0vpBsFMkOaz8qTdvOyxCx4GwaVugc3gi7V3cnSxh+v/IwJRM/D936UXxdN6PrFGcnVyZrz3noSelf9cqjD8VxKegXse3MJPdfp1OSqVN7Z+9p/ae4x/sPkG5WOQplbmRzdHJlYW0KZW5kb2JqCjI2IDAgb2JqCjw8IC9GaWx0ZXIgL0ZsYXRlRGVjb2RlIC9MZW5ndGggNjQgPj4Kc3RyZWFtCnicMzM0VDBQ0DUCEmaGJgrmRpYKKYZcQD6IlcsFE8sBs8xMzIAsY1NTJJYBkDYyNYPTEBmgAXAGRH8aAClPFE4KZW5kc3RyZWFtCmVuZG9iagoyNyAwIG9iago8PCAvRmlsdGVyIC9GbGF0ZURlY29kZSAvTGVuZ3RoIDMwNCA+PgpzdHJlYW0KeJw9kjuSwzAMQ3udghfIjPiT5PNkJ5X3/u0+MslWgEmJACgvdZmypjwgaSYJ/9Hh4WI75XfYns3MwLVELxPLKc+hK8TcRfmymY26sjrFqsMwnVv0qJyLhk2TmucqSxm3C57DtYnnln3EDzc0qAd1jUvCDd3VaFkKzXB1/zu9R9l3NTwXm1Tq1BePF1EV5vkhT6KH6UrifDwoIVx7MEYWEuRT0UCOs1yt8l5C9g63GrLCQWpJ57MnPNh1ek8ubhfNEA9kuVT4TlHs7dAzvuxKCT0StuFY7n07mrHpGps47H7vRtbKjK5oIX7IVyfrJWDcUyZFEmROtlhui9We7qEopnOGcxkg6tmKhlLmYlerfww7bywv2SzIlMwLMkanTZ44eMh+jZr0eZXneP0BbPNzOwplbmRzdHJlYW0KZW5kb2JqCjI4IDAgb2JqCjw8IC9GaWx0ZXIgL0ZsYXRlRGVjb2RlIC9MZW5ndGggMjM3ID4+CnN0cmVhbQp4nEVRSXIEIQy79yv0ganCK/CeTs2p8/9rLDNJThZgazFpgYEteIkh1sDMgS+5fE3oNHw3MtvwOtkecE+4LtyXy4JnwpbAV1SXd70vXdlIfXeHqn5mZHuzSM2QlZU69UI0JtghET0jMslWLHODpCmtUuW+KFuALuqVtk47jZKgIxThb5Qj4ekVSnZNbBqr1DqgoQjLti6IOpkkonZhcWrxliEin3VjNcf4i04idsfj/qww61EkktJnB91xJqNNll0DObl5qrBWKjmIPl7RxoTqdKqBY7zXtvQTaeC59l/hBz59/48Y+rneP8buXCIKZW5kc3RyZWFtCmVuZG9iagoyOSAwIG9iago8PCAvRmlsdGVyIC9GbGF0ZURlY29kZSAvTGVuZ3RoIDIzMCA+PgpzdHJlYW0KeJw1UUluwzAMvOsV84EA4i6/x0FP7f+vHdIJYGBoS5zNERsbEXiJwc9B5MZb1oya+JvJXfG7PBUeCbeCJ1EEXoZ72QkubxiX/TjMfPBeWjmTGk8yIBfZ9PBEyGCXQOjA7BrUYZtpJ/qGhM+OSDUbWU5fS9BLqxAoT9l+pwtKtK3qz+2zLrTta0842e2pJ5VPIJ5bsgKXjVdMFmMZ9ETlLsX0QaqzhZ6E8qJ8DrL5qCESXaKcgScGB6NAO7Dntp+JV4WgdXWfto2hGikdT/82NDVJIuQTJZzZ0rhb+P6ee/38A6ZUU58KZW5kc3RyZWFtCmVuZG9iagozMCAwIG9iago8PCAvRmlsdGVyIC9GbGF0ZURlY29kZSAvTGVuZ3RoIDIyNyA+PgpzdHJlYW0KeJw1TzuyAyEM6zmFLpAZjG1gz7OZVC/3b59ksg0S/kjy9ERHJl7myAis2fG2FhmIGfgWU/GvPe3DhOo9uIcI5eJCmGEknDXruJun48W/XeUz1sG7Db5ilhcEtjCT9ZXFmct2wVgaJ3FOshtj10RsY13r6RTWEUwoAyGd7TAlyBwVKX2yo4w5Ok7kiediqsUuv+9hfcGmMaLCHFcFT9BkUJY97yagHRf039WN30k0i14CMpFgYZ0k5s5ZTvjVa0fHUYsiMSekGeQyEdKcrmIKoQnFOjsKKhUFl+pzyt0+/2hdW00KZW5kc3RyZWFtCmVuZG9iagozMSAwIG9iago8PCAvRmlsdGVyIC9GbGF0ZURlY29kZSAvTGVuZ3RoIDI0NSA+PgpzdHJlYW0KeJxFULuNQzEM6z0FFwhg/Sx7nndIldu/PUpGcIUhWj+SWhKYiMBLDLGUb+JHRkE9C78XheIzxM8XhUHOhKRAnPUZEJl4htpGbuh2cM68wzOMOQIXxVpwptOZ9lzY5JwHJxDObZTxjEK6SVQVcVSfcUzxqrLPjdeBpbVss9OR7CGNhEtJJSaXflMq/7QpWyro2kUTsEjkgZNNNOEsP0OSYsyglFH3MLWO9HGykUd10MnZnDktmdnup+1MfA9YJplR5Smd5zI+J6nzXE597rMd0eSipVX7nP3ekZbyIrXbodXpVyVRmY3Vp5C4PP+Mn/H+A46gWT4KZW5kc3RyZWFtCmVuZG9iagozMiAwIG9iago8PCAvRmlsdGVyIC9GbGF0ZURlY29kZSAvTGVuZ3RoIDM5MiA+PgpzdHJlYW0KeJw9UktuBTEI288puECl8E1ynqne7t1/W5vMVKoKLwO2MZSXDKklP+qSiDNMfvVyXeJR8r1samfmIe4uNqb4WHJfuobYctGaYrFPHMkvyLRUWKFW3aND8YUoEw8ALeCBBeG+HP/xF6jB17CFcsN7ZAJgStRuQMZD0RlIWUERYfuRFeikUK9s4e8oIFfUrIWhdGKIDZYAKb6rDYmYqNmgh4SVkqod0vGMpPBbwV2JYVBbW9sEeGbQENnekY0RM+3RGXFZEWs/PemjUTK1URkPTWd88d0yUvPRFeik0sjdykNnz0InYCTmSZjncCPhnttBCzH0ca+WT2z3mClWkfAFO8oBA7393pKNz3vgLIxc2+xMJ/DRaaccE62+HmL9gz9sS5tcxyuHRRSovCgIftdBE3F8WMX3ZKNEd7QB1iMT1WglEAwSws7tMPJ4xnnZ3hW05vREaKNEHtSOET0ossXlnBWwp/yszbEcng8me2+0j5TMzKiEFdR2eqi2z2Md1Hee+/r8AS4AoRkKZW5kc3RyZWFtCmVuZG9iagozMyAwIG9iago8PCAvRmlsdGVyIC9GbGF0ZURlY29kZSAvTGVuZ3RoIDEzMyA+PgpzdHJlYW0KeJxNj0ESwzAIA+9+hZ6AsQHznnR6Sv5/LZA27gXtjICRhjAIPGIM6zAlvHr74VWkS3A2jvklGUU8CGoL3BdUBUdjip342N2h7KXi6RRNi+sRc9O0pHQ3USptvZ3I+MB9n94fVbYknYIeW+qELtEk8kUCc9hUMM/qxktLj6ft2d4fZj4z1wplbmRzdHJlYW0KZW5kb2JqCjM0IDAgb2JqCjw8IC9GaWx0ZXIgL0ZsYXRlRGVjb2RlIC9MZW5ndGggMjQ3ID4+CnN0cmVhbQp4nE1Ru21EMQzr3xRc4ADra3meC1Jd9m9DyQiQwiChLymnJRb2xksM4QdbD77kkVVDfx4/MewzLD3J5NQ/5rnJVBS+FaqbmFAXYuH9aAS8FnQvIivKB9+PZQxzzvfgoxCXYCY0YKxvSSYX1bwzZMKJoY7DQZtUGHdNFCyuFc0zyO1WN7I6syBseCUT4sYARATZF5DNYKOMsZWQxXIeqAqSBVpg1+kbUYuCK5TWCXSi1sS6zOCr5/Z2N0Mv8uCounh9DOtLsMLopXssfK5CH8z0TDt3SSO98KYTEWYPBVKZnZGVOj1ifbdA/59lK/j7yc/z/QsVKFwqCmVuZHN0cmVhbQplbmRvYmoKMzUgMCBvYmoKPDwgL0ZpbHRlciAvRmxhdGVEZWNvZGUgL0xlbmd0aCA5MCA+PgpzdHJlYW0KeJxNjUESwCAIA++8Ik9QRND/dHrS/1+r1A69wE4CiRZFgvQ1aksw7rgyFWtQKZiUl8BVMFwL2u6iyv4ySUydhtN7twODsvFxg9JJ+/ZxegCr/XoG3Q/SHCJYCmVuZHN0cmVhbQplbmRvYmoKMzYgMCBvYmoKPDwgL0ZpbHRlciAvRmxhdGVEZWNvZGUgL0xlbmd0aCAzMzggPj4Kc3RyZWFtCnicRVJLcsUwCNvnFFwgM+Zn4/O8Tlfp/beVcDrdPPQMCAkyPWVIptw2lmSE5BzypVdkiNWQn0aORMQQ3ymhwK7yubyWxFzIbolK8aEdP5elNzLNrtCqt0enNotGNSsj5yBDhHpW6MzuUdtkw+t2Iek6UxaHcCz/QwWylHXKKZQEbUHf2CPobxY8EdwGs+Zys7lMbvW/7lsLntc6W7FtB0AJlnPeYAYAxMMJ2gDE3NreFikoH1W6iknCrfJcJztQttCqdLw3gBkHGDlgw5KtDtdobwDDPg/0okbF9hWgqCwg/s7ZZsHeMclIsCfmBk49cTrFkXBJOMYCQIqt4hS68R3Y4i8Xroia8Al1OmVNvMKe2uLHQpMI71JxAvAiG25dHUW1bE/nCbQ/KpIzYqQexNEJkdSSzhEUlwb10Br7uIkZr43E5p6+3T/COZ/r+xcWuIPgCmVuZHN0cmVhbQplbmRvYmoKMzcgMCBvYmoKPDwgL0ZpbHRlciAvRmxhdGVEZWNvZGUgL0xlbmd0aCAxNjMgPj4Kc3RyZWFtCnicRZC5dQQxDENzVYESeIA66hk/R7P9pwtpvN5A+niEeIg9CcNyXcWF0Q0/3rbMNLyOMtyN9WXG+KixQE7QBxgiE1ejSfXtijNU6eHVYq6jolwvOiISzJLjq0AjfDqyx0Nb25l+Oq9/7CHvE/8qKuduYQEuqu5A+VIf8dSP2VHqmqGPKitrHmravwi7IpS2fVxOZZy6ewe0wmcrV/t9A6jnOoAKZW5kc3RyZWFtCmVuZG9iagozOCAwIG9iago8PCAvRmlsdGVyIC9GbGF0ZURlY29kZSAvTGVuZ3RoIDY4ID4+CnN0cmVhbQp4nDMyt1AwULA0ARKGFiYK5mYGCimGXEC+qYm5Qi4XSAzEygGzDIC0JZyCiFtCNEGUglgQpWYmZhBJOAMilwYAybQV5QplbmRzdHJlYW0KZW5kb2JqCjM5IDAgb2JqCjw8IC9GaWx0ZXIgL0ZsYXRlRGVjb2RlIC9MZW5ndGggMjU1ID4+CnN0cmVhbQp4nEWRS5IDIAhE956CI4D85DyZmlVy/+00mEw2dpeo/YRKI6YSLOcUeTD9yPLNZLbptRyrnY0CiiIUzOQq9FiB1Z0p4sy1RLX1sTJy3Okdg+IN566cVLK4UcY6qjoVOKbnyvqq7vy4LMq+I4cyBWzWOQ42cOW2YYwTo81Wd4f7RJCnk6mj4naQbPiDk8a+ytUVuE42++olGAeCfqEJTPJNoHWGQOPmKXpyCfbxcbvzQLC3vAmkbAjkyBCMDkG7Tq5/cev83v86w53n2gxXjnfxO0xru+MvMcmKuYBF7hTU8z0XresMHe/JmWNy031D51ywy91Bps/8H+v3D1CKZogKZW5kc3RyZWFtCmVuZG9iago0MCAwIG9iago8PCAvRmlsdGVyIC9GbGF0ZURlY29kZSAvTGVuZ3RoIDE2MSA+PgpzdHJlYW0KeJxFkEsSwyAMQ/ecQkfwRwZ8nnS6Su+/rSFNs4CnsUAGdycEqbUFE9EFL21Lugs+WwnOxnjoNm41EuQEdYBWpONolFJ9ucVplXTxaDZzKwutEx1mDnqUoxmgEDoV3u2i5HKm7s75R3D1X/VHse6czcTAZOUOhGb1Ke58mx1RXd1kf9JjbtZrfxX2qrC0rKXlhNvOXTOgBO6pHO39BalzOoQKZW5kc3RyZWFtCmVuZG9iago0MSAwIG9iago8PCAvRmlsdGVyIC9GbGF0ZURlY29kZSAvTGVuZ3RoIDMyMCA+PgpzdHJlYW0KeJw1UbtxxTAM6zUFF/Cd+JU0j3Ovytu/DUA7FWEaBECqvGRKuVzqklWywuRHh+oUTfk+YKb8DvWQ4+ge2SG6U9aWexgIy8Q8pY5YTZZ7uAWBLwxNibmF8/cI6CsGozATgbrF3z9AsyQwaXDwU5BrrVpiiQ48LBZYsyvMrRopVMhVfDs2uQcFcnGz0KccmhS33ILwZYhkR2qxr8tlKfK79QkYhBXmiE8UiYXngQ5mIvEnA2J79tliV1cvqhEZ1kmHB1IE0mxuEjA0RbLqgxvYV8c1P09H2cHJQb+Kwfg2OJkvSXlfBaEQjxf+Ds/ZyLGSQyQU8n21wIgjbIARoU/tIxBlIDRF9+6ZUj4mVYrvAEYhHH2qVzK8F5HZaobN/xld2SoKBlVZH59GcCaDSTjzZKMK01K107/73OPzB2NjeoAKZW5kc3RyZWFtCmVuZG9iago0MiAwIG9iago8PCAvRmlsdGVyIC9GbGF0ZURlY29kZSAvTGVuZ3RoIDIxNCA+PgpzdHJlYW0KeJw9ULsRQzEI6z0FC+TOfO03z8uly/5tJJykQjZCEpSaTMmUhzrKkqwpTx0+S2KHvIflbmQ2JSpFL5OwJffQCvF9ieYU993VlrNDNJdoOX4LMyqqGx3TSzaacCoTuqDcwzP6DW10A1aHHrFbINCkYNe2IHLHDxgMwZkTiyIMSk0G/61y91Lc7z0cb6KIlHTwrvnl9MvPLbxOPY5Eur35imtxpjoKRHBGavKKdGHFsshDpNUENT0Da7UArt56+TdoR3QZgOwTieM0pRxD/9a4x+sDh4pS9AplbmRzdHJlYW0KZW5kb2JqCjQzIDAgb2JqCjw8IC9GaWx0ZXIgL0ZsYXRlRGVjb2RlIC9MZW5ndGggODAgPj4Kc3RyZWFtCnicRYy7DcAwCER7pmAEfiZmnyiVs38bIErccE+6e7g6EjJT3mGGhwSeDCyGU/EGmaNgNbhGUo2d7KOwbl91geZ6U6v19wcqT3Z2cT3Nyxn0CmVuZHN0cmVhbQplbmRvYmoKNDQgMCBvYmoKPDwgL0ZpbHRlciAvRmxhdGVEZWNvZGUgL0xlbmd0aCAzMzIgPj4Kc3RyZWFtCnicLVI5jiQxDMv9Cn5gAOvy8Z4eTNT7/3RJVQUFqmzLPORyw0QlfiyQ21Fr4tdGZqDC8K+rzIXvSNvIOohryEVcyZbCZ0Qs5DHEPMSC79v4GR75rMzJswfGL9n3GVbsqQnLQsaLM7TDKo7DKsixYOsiqnt4U6TDqSTY44v/PsVzF4IWviNowC/556sjeL6kRdo9Ztu0Ww+WaUeVFJaD7WnOy+RL6yxXx+P5INneFTtCaleAojB3xnkujjJtZURrYWeDpMbF9ubYj6UEXejGZaQ4AvmZKsIDSprMbKIg/sjpIacyEKau6Uont1EVd+rJXLO5vJ1JMlv3RYrNFM7rwpn1d5gyq807eZYTpU5F+Bl7tgQNnePq2WuZhUa3OcErJXw2dnpy8r2aWQ/JqUhIFdO6Ck6jyBRL2Jb4moqa0tTL8N+X9xl//wEz4nwBCmVuZHN0cmVhbQplbmRvYmoKNDUgMCBvYmoKPDwgL0ZpbHRlciAvRmxhdGVEZWNvZGUgL0xlbmd0aCAzMTcgPj4Kc3RyZWFtCnicNVJLckMxCNu/U3CBzpi/fZ50smruv62EJyuwLUBCLi9Z0kt+1CXbpcPkVx/3JbFCPo/tmsxSxfcWsxTPLa9HzxG3LQoEURM9+DInFSLUz9ToOnhhlz4DrxBOKRZ4B5MABq/hX3iUToPAOxsy3hGTkRoQJMGaS4tNSJQ9Sfwr5fWklTR0fiYrc/l7cqkUaqPJCBUgWLnYB6QrKR4kEz2JSLJyvTdWiN6QV5LHZyUmGRDdJrFNtMDj3JW0hJmYQgXmWIDVdLO6+hxMWOOwhPEqYRbVg02eNamEZrSOY2TDePfCTImFhsMSUJt9lQmql4/T3AkjpkdNdu3Csls27yFEo/kzLJTBxygkAYdOYyQK0rCAEYE5vbCKveYLORbAiGWdmiwMbWglu3qOhcDQnLOlYcbXntfz/gdFW3ujCmVuZHN0cmVhbQplbmRvYmoKNDYgMCBvYmoKPDwgL0ZpbHRlciAvRmxhdGVEZWNvZGUgL0xlbmd0aCAxNyA+PgpzdHJlYW0KeJwzNrRQMIDDFEMuABqUAuwKZW5kc3RyZWFtCmVuZG9iago0NyAwIG9iago8PCAvRmlsdGVyIC9GbGF0ZURlY29kZSAvTGVuZ3RoIDEzMSA+PgpzdHJlYW0KeJxFj8sNBCEMQ+9U4RLyGT6ph9We2P6v6zCaQUL4QSI78TAIrPPyNtDF8NGiwzf+NtWrY5UsH7p6UlYP6ZCHvPIVUGkwUcSFWUwdQ2HOmMrIljK3G+G2TYOsbJVUrYN2PAYPtqdlqwh+qW1h6izxDMJVXrjHDT+QS613vVW+f0JTMJcKZW5kc3RyZWFtCmVuZG9iago0OCAwIG9iago8PCAvRmlsdGVyIC9GbGF0ZURlY29kZSAvTGVuZ3RoIDMzOCA+PgpzdHJlYW0KeJw1Ujmu3UAM630KXSCAds2c5wWpfu7fhpRfCkO0VoqajhaVafllIVUtky6/7UltiRvy98kKiROSVyXapQyRUPk8hVS/Z8u8vtacESBLlQqTk5LHJQv+DJfeLhznY2s/jyN3PXpgVYyEEgHLFBOja1k6u8Oajfw8pgE/4hFyrli3HGMVSA26cdoV70PzecgaIGaYlooKXVaJFn5B8aBHrX33WFRYINHtHElwjI1QkYB2gdpIDDmzFruoL/pZlJgJdO2LIu6iwBJJzJxiXTr6Dz50LKi/NuPLr45K+kgra0zad6NJacwik66XRW83b309uEDzLsp/Xs0gQVPWKGl80KqdYyiaGWWFdxyaDDTHHIfMEzyHMxKU9H0ofl9LJrookT8ODaF/Xx6jjJwGbwFz0Z+2igMX8dlhrxxghdLFmuR9QCoTemD6/9f4ef78Axy2gFQKZW5kc3RyZWFtCmVuZG9iago0OSAwIG9iago8PCAvRmlsdGVyIC9GbGF0ZURlY29kZSAvTGVuZ3RoIDI0OCA+PgpzdHJlYW0KeJwtUTmSA0EIy+cVekJz0++xy5H3/+kKygGDhkMgOi1xUMZPEJYr3vLIVbTh75kYwXfBod/KdRsWORAVSNIYVE2oXbwevQd2HGYC86Q1LIMZ6wM/Ywo3enF4TMbZ7XUZNQR712tPZlAyKxdxycQFU3XYyJnDT6aMC+1czw3IuRHWZRikm5XGjIQjTSFSSKHqJqkzQZAEo6tRo40cxX7pyyOdYVUjagz7XEvb13MTzho0OxarPDmlR1ecy8nFCysH/bzNwEVUGqs8EBJwv9tD/Zzs5Dfe0rmzxfT4XnOyvDAVWPHmtRuQTbX4Ny/i+D3j6/n8A6ilWxYKZW5kc3RyZWFtCmVuZG9iago1MCAwIG9iago8PCAvRmlsdGVyIC9GbGF0ZURlY29kZSAvTGVuZ3RoIDE3MSA+PgpzdHJlYW0KeJxNkE0OQiEQg/ecohcwofMDj/NoXOn9t3bw+eKC9EshQ6fDAx1H4kZHhs7oeLDJMQ68CzImXo3zn4zrJI4J6hVtwbq0O+7NLDEnLBMjYGuU3JtHFPjhmAtBguzywxcYRKRrmG81n3WTfn67013UpXX30yMKnMiOUAwbcAXY0z0O3BLO75omv1QpGZs4lA9UF5Gy2QmFqKVil1NVaIziVj3vi17t+QHB9jv7CmVuZHN0cmVhbQplbmRvYmoKNTEgMCBvYmoKPDwgL0ZpbHRlciAvRmxhdGVEZWNvZGUgL0xlbmd0aCA4OCA+PgpzdHJlYW0KeJw1jLsRwDAIQ3tPwQgGi4/3yaVK9m+D7dCApHf3goM6QfK4GymcLm7ZV3obj5OeJgCx9ExD7d9gRdWLWhQtX25j0GIqvj/6JCCWdfJeOPSQEt4fxRcdewplbmRzdHJlYW0KZW5kb2JqCjUyIDAgb2JqCjw8IC9GaWx0ZXIgL0ZsYXRlRGVjb2RlIC9MZW5ndGggODcgPj4Kc3RyZWFtCnicNU25EcAwCOuZghHMo9jsk0vl7N8G7LhBOn0glBtr5AGC4Z1vIfimLxmEdQhPKrslOmyhhrMKkonhVzZ4Va6K9rWSiexspjHYoGX60c63Sc8Hpd4bmAplbmRzdHJlYW0KZW5kb2JqCjUzIDAgb2JqCjw8IC9GaWx0ZXIgL0ZsYXRlRGVjb2RlIC9MZW5ndGggMjEwID4+CnN0cmVhbQp4nDVQyw1DMQi7ZwoWqBQCgWSeVr11/2tt0DthEf9CWMiUCHmpyc4p6Us+OkwPti6/sSILrXUl7MqaIJ4r76GZsrHR2OJgcBomXoAWN2DoaY0aNXThgqYulUKBxSXwmXx1e+i+Txl4ahlydgQRQ8lgCWq6Fk1YtDyfkE4B4v9+w+4t5KGS88qeG/kbnO3wO7Nu4SdqdiLRchUy1LM0xxgIE0UePHlFpnDis9Z31TQS1GYLTpYBrk4/jA4AYCJeWYDsrkQ5S9KOpZ9vvMf3D0AAU7QKZW5kc3RyZWFtCmVuZG9iagoxOSAwIG9iago8PCAvQmFzZUZvbnQgL0RlamFWdVNhbnMgL0NoYXJQcm9jcyAyMCAwIFIKL0VuY29kaW5nIDw8Ci9EaWZmZXJlbmNlcyBbIDMyIC9zcGFjZSA0OCAvemVybyAvb25lIC90d28gL3RocmVlIC9mb3VyIC9maXZlIC9zaXggNTYgL2VpZ2h0IC9uaW5lIDY4Ci9EIDcwIC9GIC9HIDc2IC9MIDgyIC9SIDg0IC9UIDk3IC9hIC9iIC9jIC9kIC9lIC9mIC9nIC9oIC9pIDEwOSAvbSAvbiAvbwoxMTUgL3MgL3QgL3UgMTE5IC93IC94IF0KL1R5cGUgL0VuY29kaW5nID4+Ci9GaXJzdENoYXIgMCAvRm9udEJCb3ggWyAtMTAyMSAtNDYzIDE3OTQgMTIzMyBdIC9Gb250RGVzY3JpcHRvciAxOCAwIFIKL0ZvbnRNYXRyaXggWyAwLjAwMSAwIDAgMC4wMDEgMCAwIF0gL0xhc3RDaGFyIDI1NSAvTmFtZSAvRGVqYVZ1U2FucwovU3VidHlwZSAvVHlwZTMgL1R5cGUgL0ZvbnQgL1dpZHRocyAxNyAwIFIgPj4KZW5kb2JqCjE4IDAgb2JqCjw8IC9Bc2NlbnQgOTI5IC9DYXBIZWlnaHQgMCAvRGVzY2VudCAtMjM2IC9GbGFncyAzMgovRm9udEJCb3ggWyAtMTAyMSAtNDYzIDE3OTQgMTIzMyBdIC9Gb250TmFtZSAvRGVqYVZ1U2FucyAvSXRhbGljQW5nbGUgMAovTWF4V2lkdGggMTM0MiAvU3RlbVYgMCAvVHlwZSAvRm9udERlc2NyaXB0b3IgL1hIZWlnaHQgMCA+PgplbmRvYmoKMTcgMCBvYmoKWyA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMAo2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDYwMCA2MDAgNjAwIDMxOCA0MDEgNDYwIDgzOCA2MzYKOTUwIDc4MCAyNzUgMzkwIDM5MCA1MDAgODM4IDMxOCAzNjEgMzE4IDMzNyA2MzYgNjM2IDYzNiA2MzYgNjM2IDYzNiA2MzYgNjM2CjYzNiA2MzYgMzM3IDMzNyA4MzggODM4IDgzOCA1MzEgMTAwMCA2ODQgNjg2IDY5OCA3NzAgNjMyIDU3NSA3NzUgNzUyIDI5NQoyOTUgNjU2IDU1NyA4NjMgNzQ4IDc4NyA2MDMgNzg3IDY5NSA2MzUgNjExIDczMiA2ODQgOTg5IDY4NSA2MTEgNjg1IDM5MCAzMzcKMzkwIDgzOCA1MDAgNTAwIDYxMyA2MzUgNTUwIDYzNSA2MTUgMzUyIDYzNSA2MzQgMjc4IDI3OCA1NzkgMjc4IDk3NCA2MzQgNjEyCjYzNSA2MzUgNDExIDUyMSAzOTIgNjM0IDU5MiA4MTggNTkyIDU5MiA1MjUgNjM2IDMzNyA2MzYgODM4IDYwMCA2MzYgNjAwIDMxOAozNTIgNTE4IDEwMDAgNTAwIDUwMCA1MDAgMTM0MiA2MzUgNDAwIDEwNzAgNjAwIDY4NSA2MDAgNjAwIDMxOCAzMTggNTE4IDUxOAo1OTAgNTAwIDEwMDAgNTAwIDEwMDAgNTIxIDQwMCAxMDIzIDYwMCA1MjUgNjExIDMxOCA0MDEgNjM2IDYzNiA2MzYgNjM2IDMzNwo1MDAgNTAwIDEwMDAgNDcxIDYxMiA4MzggMzYxIDEwMDAgNTAwIDUwMCA4MzggNDAxIDQwMSA1MDAgNjM2IDYzNiAzMTggNTAwCjQwMSA0NzEgNjEyIDk2OSA5NjkgOTY5IDUzMSA2ODQgNjg0IDY4NCA2ODQgNjg0IDY4NCA5NzQgNjk4IDYzMiA2MzIgNjMyIDYzMgoyOTUgMjk1IDI5NSAyOTUgNzc1IDc0OCA3ODcgNzg3IDc4NyA3ODcgNzg3IDgzOCA3ODcgNzMyIDczMiA3MzIgNzMyIDYxMSA2MDUKNjMwIDYxMyA2MTMgNjEzIDYxMyA2MTMgNjEzIDk4MiA1NTAgNjE1IDYxNSA2MTUgNjE1IDI3OCAyNzggMjc4IDI3OCA2MTIgNjM0CjYxMiA2MTIgNjEyIDYxMiA2MTIgODM4IDYxMiA2MzQgNjM0IDYzNCA2MzQgNTkyIDYzNSA1OTIgXQplbmRvYmoKMjAgMCBvYmoKPDwgL0QgMjEgMCBSIC9GIDIyIDAgUiAvRyAyMyAwIFIgL0wgMjQgMCBSIC9SIDI1IDAgUiAvVCAyNiAwIFIgL2EgMjcgMCBSCi9iIDI4IDAgUiAvYyAyOSAwIFIgL2QgMzAgMCBSIC9lIDMxIDAgUiAvZWlnaHQgMzIgMCBSIC9mIDMzIDAgUgovZml2ZSAzNCAwIFIgL2ZvdXIgMzUgMCBSIC9nIDM2IDAgUiAvaCAzNyAwIFIgL2kgMzggMCBSIC9tIDM5IDAgUiAvbiA0MCAwIFIKL25pbmUgNDEgMCBSIC9vIDQyIDAgUiAvb25lIDQzIDAgUiAvcyA0NCAwIFIgL3NpeCA0NSAwIFIgL3NwYWNlIDQ2IDAgUgovdCA0NyAwIFIgL3RocmVlIDQ4IDAgUiAvdHdvIDQ5IDAgUiAvdSA1MCAwIFIgL3cgNTEgMCBSIC94IDUyIDAgUgovemVybyA1MyAwIFIgPj4KZW5kb2JqCjMgMCBvYmoKPDwgL0YxIDE5IDAgUiAvRjIgMTQgMCBSID4+CmVuZG9iago0IDAgb2JqCjw8IC9BMSA8PCAvQ0EgMCAvVHlwZSAvRXh0R1N0YXRlIC9jYSAxID4+Ci9BMiA8PCAvQ0EgMSAvVHlwZSAvRXh0R1N0YXRlIC9jYSAxID4+ID4+CmVuZG9iago1IDAgb2JqCjw8ID4+CmVuZG9iago2IDAgb2JqCjw8ID4+CmVuZG9iago3IDAgb2JqCjw8ID4+CmVuZG9iagoyIDAgb2JqCjw8IC9Db3VudCAxIC9LaWRzIFsgMTAgMCBSIF0gL1R5cGUgL1BhZ2VzID4+CmVuZG9iago1NCAwIG9iago8PCAvQ3JlYXRpb25EYXRlIChEOjIwMTkwNjI3MTQ0MzI1KzAyJzAwJykKL0NyZWF0b3IgKG1hdHBsb3RsaWIgMy4wLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZykKL1Byb2R1Y2VyIChtYXRwbG90bGliIHBkZiBiYWNrZW5kIDMuMC4zKSA+PgplbmRvYmoKeHJlZgowIDU1CjAwMDAwMDAwMDAgNjU1MzUgZiAKMDAwMDAwMDAxNiAwMDAwMCBuIAowMDAwMDE1MjU3IDAwMDAwIG4gCjAwMDAwMTUwNTIgMDAwMDAgbiAKMDAwMDAxNTA5NSAwMDAwMCBuIAowMDAwMDE1MTk0IDAwMDAwIG4gCjAwMDAwMTUyMTUgMDAwMDAgbiAKMDAwMDAxNTIzNiAwMDAwMCBuIAowMDAwMDAwMDY1IDAwMDAwIG4gCjAwMDAwMDAzOTcgMDAwMDAgbiAKMDAwMDAwMDIwOCAwMDAwMCBuIAowMDAwMDAyMzA0IDAwMDAwIG4gCjAwMDAwMDMwMDAgMDAwMDAgbiAKMDAwMDAwMjc5MiAwMDAwMCBuIAowMDAwMDAyNDc3IDAwMDAwIG4gCjAwMDAwMDQwNTMgMDAwMDAgbiAKMDAwMDAwMjMyNSAwMDAwMCBuIAowMDAwMDEzNjE3IDAwMDAwIG4gCjAwMDAwMTM0MTcgMDAwMDAgbiAKMDAwMDAxMjk1NiAwMDAwMCBuIAowMDAwMDE0NjcwIDAwMDAwIG4gCjAwMDAwMDQwODUgMDAwMDAgbiAKMDAwMDAwNDMxOCAwMDAwMCBuIAowMDAwMDA0NDY0IDAwMDAwIG4gCjAwMDAwMDQ3ODEgMDAwMDAgbiAKMDAwMDAwNDkxMiAwMDAwMCBuIAowMDAwMDA1MjEyIDAwMDAwIG4gCjAwMDAwMDUzNDggMDAwMDAgbiAKMDAwMDAwNTcyNSAwMDAwMCBuIAowMDAwMDA2MDM1IDAwMDAwIG4gCjAwMDAwMDYzMzggMDAwMDAgbiAKMDAwMDAwNjYzOCAwMDAwMCBuIAowMDAwMDA2OTU2IDAwMDAwIG4gCjAwMDAwMDc0MjEgMDAwMDAgbiAKMDAwMDAwNzYyNyAwMDAwMCBuIAowMDAwMDA3OTQ3IDAwMDAwIG4gCjAwMDAwMDgxMDkgMDAwMDAgbiAKMDAwMDAwODUyMCAwMDAwMCBuIAowMDAwMDA4NzU2IDAwMDAwIG4gCjAwMDAwMDg4OTYgMDAwMDAgbiAKMDAwMDAwOTIyNCAwMDAwMCBuIAowMDAwMDA5NDU4IDAwMDAwIG4gCjAwMDAwMDk4NTEgMDAwMDAgbiAKMDAwMDAxMDEzOCAwMDAwMCBuIAowMDAwMDEwMjkwIDAwMDAwIG4gCjAwMDAwMTA2OTUgMDAwMDAgbiAKMDAwMDAxMTA4NSAwMDAwMCBuIAowMDAwMDExMTc0IDAwMDAwIG4gCjAwMDAwMTEzNzggMDAwMDAgbiAKMDAwMDAxMTc4OSAwMDAwMCBuIAowMDAwMDEyMTEwIDAwMDAwIG4gCjAwMDAwMTIzNTQgMDAwMDAgbiAKMDAwMDAxMjUxNCAwMDAwMCBuIAowMDAwMDEyNjczIDAwMDAwIG4gCjAwMDAwMTUzMTcgMDAwMDAgbiAKdHJhaWxlcgo8PCAvSW5mbyA1NCAwIFIgL1Jvb3QgMSAwIFIgL1NpemUgNTUgPj4Kc3RhcnR4cmVmCjE1NDcxCiUlRU9GCg==\n",
"image/svg+xml": [
"\n",
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"import matplotlib.pyplot as plt\n",
"import numpy as np\n",
"import timeit\n",
"%matplotlib inline\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.barh(n, 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",
"This notebook is provided as [Open Educational Resource](https://en.wikipedia.org/wiki/Open_educational_resources). Feel free to use the notebook for your own 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: *Sascha Spors, Continuous- and Discrete-Time Signals and Systems - Theory and Computational Examples*."
]
}
],
"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.7.3"
}
},
"nbformat": 4,
"nbformat_minor": 1
}