{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "nbsphinx": "hidden"
   },
   "source": [
    "# Realization of Non-Recursive Filters\n",
    "\n",
    "*This jupyter notebook is part of a [collection of notebooks](../index.ipynb) on various topics of Digital Signal Processing. Please direct questions and suggestions to [Sascha.Spors@uni-rostock.de](mailto:Sascha.Spors@uni-rostock.de).*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Fast Convolution\n",
    "\n",
    "The straightforward convolution of two finite-length signals $x[k]$ and $h[k]$ is a numerically complex task. This has led to the development of various techniques with considerably lower complexity. The basic concept of the *fast convolution* is to exploit the correspondence between the convolution and the scalar multiplication in the frequency domain."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Convolution of Finite-Length Signals\n",
    "\n",
    "The convolution of a causal signal $x_L[k]$ of length $L$ with a causal impulse response $h_N[k]$ of length $N$ is given as\n",
    "\n",
    "\\begin{equation}\n",
    "y[k] = x_L[k] * h_N[k] = \\sum_{\\kappa = 0}^{L-1} x_L[\\kappa] \\; h_N[k - \\kappa] = \\sum_{\\kappa = 0}^{N-1} h_N[\\kappa] \\; x_L[k - \\kappa]\n",
    "\\end{equation}\n",
    "\n",
    "where $x_L[k] = 0$ for $k<0 \\wedge k \\geq L$ and $h_N[k] = 0$ for $k<0 \\wedge k \\geq N$. The resulting signal $y[k]$ is of finite length $M = N+L-1$. The computation of $y[k]$ for $k=0,1, \\dots, M-1$ requires $M \\cdot N$ multiplications and $M \\cdot (N-1)$ additions. The computational complexity of the convolution is consequently [in the order of](https://en.wikipedia.org/wiki/Big_O_notation) $\\mathcal{O}(M \\cdot N)$. Discrete-time Fourier transformation (DTFT) of above relation yields\n",
    "\n",
    "\\begin{equation}\n",
    "Y(e^{j \\Omega}) = X_L(e^{j \\Omega}) \\cdot H_N(e^{j \\Omega})\n",
    "\\end{equation}\n",
    "\n",
    "Discarding the effort of transformation, the computationally complex convolution is replaced by a scalar multiplication with respect to the frequency $\\Omega$. However, $\\Omega$ is a continuous frequency variable which limits the numerical evaluation of this scalar multiplication. In practice, the DTFT is replaced by the discrete Fourier transformation (DFT). Two aspects have to be considered before a straightforward application of the DFT\n",
    "\n",
    "1. The DFTs $X_L[\\mu]$ and $H_N[\\mu]$ are of length $L$ and $N$ respectively and cannot be multiplied straightforwardly\n",
    "    \n",
    "2. For $N = L$, the multiplication of the two spectra $X_L[\\mu]$ and $H_L[\\mu]$ would result in the [periodic/circular convolution](https://en.wikipedia.org/wiki/Circular_convolution) $x_L[k] \\circledast h_L[k]$ due to the periodicity of the DFT. Since we aim at realizing the linear convolution $x_L[k] * h_N[k]$ with the DFT, special care has to be taken to avoid cyclic effects."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Linear Convolution by Periodic Convolution\n",
    "\n",
    "The periodic convolution of the two signals $x_L[k]$ and $h_N[k]$ is defined as\n",
    "\n",
    "\\begin{equation}\n",
    "x_L[k] \\circledast h_N[k] = \\sum_{\\kappa=0}^{M-1} \\tilde{x}_M[k - \\kappa] \\; \\tilde{h}_M[\\kappa]\n",
    "\\end{equation}\n",
    "\n",
    "where the periodic continuations $\\tilde{x}_M[k]$ of $x_L[k]$ and $\\tilde{h}_M[k]$ of $h_N[k]$ with period $M$ are given as\n",
    "\n",
    "\\begin{align}\n",
    "\\tilde{x}_M[k] &= \\sum_{m = -\\infty}^{\\infty} x_L[m \\cdot M + k] \\\\\n",
    "\\tilde{h}_M[k] &= \\sum_{m = -\\infty}^{\\infty} h_N[m \\cdot M + k]\n",
    "\\end{align}\n",
    "\n",
    "The result of the circular convolution has a periodicity of $M$.\n",
    "\n",
    "To compute the linear convolution by the periodic convolution one has to take care that the result of the linear convolution fits into one period of the periodic convolution. Hence, the periodicity has to be chosen as $M \\geq N+L-1$. This can be achieved by zero-padding of $x_L[k]$ and $h_N[k]$ to a total length of $M$\n",
    "\n",
    "\\begin{align}\n",
    "x_M[k] &= \\begin{cases}\n",
    "x_L[k] & \\mathrm{for} \\; k=0, 1, \\dots, L-1 \\\\\n",
    "0 & \\mathrm{for} \\; k=L, L+1, \\dots, M-1\n",
    "\\end{cases}\n",
    "\\\\\n",
    "h_M[k] &= \\begin{cases}\n",
    "h_N[k] & \\mathrm{for} \\; k=0, 1, \\dots, N-1 \\\\\n",
    "0 & \\mathrm{for} \\; k=N, N+1, \\dots, M-1\n",
    "\\end{cases}\n",
    "\\end{align}\n",
    "\n",
    "This results in the desired equality of linear and periodic convolution\n",
    "\n",
    "\\begin{equation}\n",
    "x_L[k] * h_N[k] = x_M[k] \\circledast h_M[k]\n",
    "\\end{equation}\n",
    "\n",
    "for $k = 0,1,\\dots, M-1$ with $M = N+L-1$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Example - Linear by periodic convolution\n",
    "\n",
    "The following example computes the linear, periodic and linear by periodic convolution of a rectangular signal $x[k] = \\text{rect}_L[k]$ of length $L$ with a triangular signal $h[k] = \\Lambda_N[k]$ of length $N$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmgAAADkCAYAAAAsAtjDAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAFFlJREFUeJzt3X/Un3V93/Hnq0lQVouUxnliEknapUomIJ4M8ehOqdpjoD1grXOw0qp1jW7F407RCuuGyuo2y5kyV6amrWJlAtFZmrmsmRU83agiwWggyaKR1ZGA4g9AW/kR4nt/fK/Il5s7ub833vd1f775Ph/n3Cff67re38/1vrhObl65fqaqkCRJUjt+bKEbkCRJ0mMZ0CRJkhpjQJMkSWqMAU2SJKkxBjRJkqTGGNAkSZIaY0CTJElqjAFNkiSpMQY0SWMpyc4kZ/awnr9O8tIn8L1K8rdJ3jll/p1JTpum/oYkDyb53z9Kv5KODgY0Sc1K8qIkf5Xk/iTfSXJTkn8AUFV/v6o+s8AtzuTUqvrdQxNJfhJYBuyeWlhVLwbe0GNvkhq2eKEbkKTpJDkO+CTwz4BNwDHAPwQeWsi+fkQnA3ur6sGFbkRS2zyCJqlVPwtQVddU1cGqeqCq/mdV7YDHnnpM8rwk25N8L8nHklyX5PcODdTVvjnJju5o3HVJnjy0/OIkX+2+vyvJL4/SYJLfT3L90PTlST6d5JjDfOUU4Pau9u8k+WiSTyR5ymz/40g6uhnQJLXqy8DBJB9OclZ3evBxujD0p8BVwAnANcB0AetVwHpgNYOg9JqhZV9lcHTuqcA7gKuTLBuhx3cBP5/ktCRv6MZ/RVU9fJj6k4HbkqwGbgL2AL9SVX8zwrokTRADmqQmVdV3gRcBBfwh8M0km5M8fUrpGQwu13hvVR2oqk8An59myPdW1V1V9R3gvwHPHVrXx7plP6iq64CvAKeP0OO3gfcAHwYuAc6uqvuP8JVTGFyDdiPw9qp6R1XVTOuRNHkMaJKaVVW7q+o1VbUCeA7wDOCKKWXPAPZPCTp3TjPc14c+fx/44WnFJL+e5ItJ7ktyX7eupSO2uZ3BkbFLqmq69R5aR7pxfxl4X1X92YjjS5pABjRJY6Gq/g+D05jPmbLobmB5F4AOWTnquElOZHCE7kLgp6rqeAbXieWIXxx892TgfQyOoP3GDOWruz9fClyUZN2oPUqaPAY0SU1K8uwkFyVZ0U2vBM4HPjel9LPAQeDCJIuTnMsIpyeH/DiD06jf7NbzWh4fAqfrbzmDU6VvAP45cPIMz2U7BdhRVbcBG4A/HfE6N0kTyIAmqVXfA54P3JzkbxkEs9uBi4aLugvyXwG8DrgPuIDB4zlGehxHVe0C/gODoPcNBqcrbzrSd7pHgGwB3l1Vm6vq+8DlwDuP8LWTgR3dOq8HNgLXD99NKkmHxOtTJR1tktwMvL+qPrSAPTzIICS+t6r+9Qj1n2Jww8Pnq+ol892fpLYZ0CSNvSQ/x+CRFd8CfhV4P/DTVXX3gjYmSU+QbxKQdDR4FoO3Dfw4cAfwSsOZpHHmETRJkqTGeJOAJElSYwxokiRJjRn7a9CWLl1aq1atWug2JEmSZnTrrbd+q6qeNlPd2Ae0VatWsW3btoVuQ5IkaUZJvjZKnac4JUmSGmNAkyRJaowBTZIkqTEGNEmSpMYY0CRJkhpjQJMkSWqMAU2SJKkxBjRJkqTGGNAkSZIa01tAS/LBJPckuf0wy5PkvUn2JtmR5Hl99SZJktSSPl/1dBXwB8CfHGb5WcCa7uf5wPu6P4/otv3388J/fwNvedmzePlpyw9bd/32/Vy+dQ933fcAzzj+2CPWz6Z2vuvHdWx7meztbKmXSdnOlnqZ7diSHi9V1d/KklXAJ6vqOdMs+wDwmaq6ppveA5xZVXcfacwnLVtTy159BccuWcS/e8XJ0/4SuH77fi75xG08cODgD+cdrn42tfNdP65j28tkb2dLvUzKdrbUy2zHliZNklurat1MdS1dg7YcuHNoel83byQPHDjI5Vv3TLvs8q17HvPL4kj1s6md7/pxHdteJns7W+plUrazpV5mO7ak6bUU0EaWZEOSbUm2Dc+/674Hpq2fzfy5GGMcepmU7Wypl0nZzpZ6mZTtbKmX2Y4haXotBbT9wMqh6RXdvMepqo1VtW7qIcJnHH/stAPPZv5cjDEOvUzKdrbUy6RsZ0u9TMp2ttTLbMeQNL2WAtpm4Ne7uznPAO6f6fqzYccuWcRbXvasaZe95WXP4tgli0aqn03tfNeP69j2Mtnb2VIvk7KdLfUy27ElTa+3uziTXAOcCSxNsg94G7AEoKreD2wBzgb2At8HXjvq2MtnuEvo0Pzf+fgOHj74gyPWz6Z2vuvHdWx7meztbKmXSdnOlnqZ7diSptfrXZzz4YQTT6rvfG33SLX/+AOfBeC6179gTmvnu35cx7aX/se2l/7Htpe5GVuaFON4F6ckSZIwoEmSJDXHgCZJktQYA5okSVJjDGiSJEmNMaBJkiQ1xoAmSZLUGAOaJElSYwxokiRJjTGgSZIkNcaAJkmS1BgDmiRJUmMMaJIkSY0xoEmSJDXGgCZJktQYA5okSVJjDGiSJEmNMaBJkiQ1xoAmSZLUGAOaJElSYwxokiRJjek1oCVZn2RPkr1JLp5m+TOT3Jhke5IdSc7usz9JkqQW9BbQkiwCrgTOAtYC5ydZO6XsXwGbquo04DzgP/fVnyRJUiv6PIJ2OrC3qu6oqoeBa4Fzp9QUcFz3+anAXT32J0mS1IQ+A9py4M6h6X3dvGFvBy5Isg/YArxxuoGSbEiyLcm2AwcOzEevkiRJC6a1mwTOB66qqhXA2cBHkjyux6raWFXrqmrdkiVLem9SkiRpPvUZ0PYDK4emV3Tzhr0O2ARQVZ8Fngws7aU7SZKkRvQZ0G4B1iRZneQYBjcBbJ5S8/+AlwAkOYlBQPtmjz1KkiQtuN4CWlU9AlwIbAV2M7hbc2eSy5Kc05VdBPxmki8B1wCvqarqq0dJkqQWLO5zZVW1hcHF/8PzLh36vAt4YZ89SZIktaa1mwQkSZImngFNkiSpMQY0SZKkxhjQJEmSGmNAkyRJaowBTZIkqTEGNEmSpMYY0CRJkhpjQJMkSWqMAU2SJKkxBjRJkqTGGNAkSZIaY0CTJElqjAFNkiSpMQY0SZKkxhjQJEmSGmNAkyRJaowBTZIkqTEGNEmSpMYY0CRJkhrTa0BLsj7JniR7k1x8mJpXJdmVZGeSj/bZnyRJUgsW97WiJIuAK4FfAPYBtyTZXFW7hmrWAJcAL6yqe5P83b76kyRJakWfR9BOB/ZW1R1V9TBwLXDulJrfBK6sqnsBquqeHvuTJElqwoxH0JKcMMI4P6iq+2aoWQ7cOTS9D3j+lJqf7dZ5E7AIeHtV/fkI65ckSTpqjHKK867uJ0eoWQQ8c476WQOcCawA/jLJyVPDX5INwAaApyz7mTlYrSRJUjtGCWi7q+q0IxUk2T7COPuBlUPTK7p5w/YBN1fVAeD/Jvkyg8B2y3BRVW0ENgKccOJJNcK6JUmSxsYo16C9ACDJ701d0F34/8OaGdwCrEmyOskxwHnA5ik11zM4ekaSpQxOed4xwtiSJElHjRkDWlU92H1cnuSfHJrf3WH5F1NqjjTOI8CFwFZgN7CpqnYmuSzJOV3ZVuDbSXYBNwJvqapvz2aDJEmSxt1sHrPxemBrkr1AAR8C3jqblVXVFmDLlHmXDn0u4Le7H0mSpIk0yl2cfwJ8AdgO/BbwUeAR4OVVtXd+25MkSZo8o1yDdhWDOzhfC1wNrALuBS5I8sp560ySJGlCzXgErapuAG44NJ1kMXAScCqD55h9fN66kyRJmkCzftVTd7H/bd3P1XPekSRJ0oSb8RRnki/MRY0kSZJGM8oRtJOS7DjC8gBPnaN+JEmSJt4oAe3ZI9Qc/FEbkSRJ0sAoNwl8DSDJp4A3V9WX5r0rSZKkCTbKYzYOeStwRZIPJVk2Xw1JkiRNupEDWlV9oap+Hvgk8OdJ3pbk2PlrTZIkaTLN5ggaSQLsAd4HvBH4SpJfm4/GJEmSJtXIAS3JTcB+4D3AcuA1wJnA6Uk2zkdzkiRJk2g2D6rdAOzqXmg+7I1Jds9hT5IkSRNt5IBWVTuPsPgX56AXSZIkMctr0A6nqu6Yi3EkSZI0RwFNkiRJc8eAJkmS1BgDmiRJUmMMaJIkSY0xoEmSJDXGgCZJktSYXgNakvVJ9iTZm+TiI9T9SpJKsq7P/iRJklrQW0BLsgi4EjgLWAucn2TtNHU/AbwJuLmv3iRJklrS5xG004G9VXVHVT0MXAucO03dvwHeBTzYY2+SJEnN6DOgLQfuHJre1837oSTPA1ZW1X8/0kBJNiTZlmTbgQMH5r5TSZKkBdTMTQJJfgx4N3DRTLVVtbGq1lXVuiVLlsx/c5IkST3qM6DtB1YOTa/o5h3yE8BzgM8k+WvgDGCzNwpIkqRJ02dAuwVYk2R1kmOA84DNhxZW1f1VtbSqVlXVKuBzwDlVta3HHiVJkhZcbwGtqh4BLgS2AruBTVW1M8llSc7pqw9JkqTWLe5zZVW1BdgyZd6lh6k9s4+eJEmSWtPMTQKSJEkaMKBJkiQ1xoAmSZLUGAOaJElSYwxokiRJjTGgSZIkNcaAJkmS1BgDmiRJUmMMaJIkSY0xoEmSJDXGgCZJktQYA5okSVJjDGiSJEmNMaBJkiQ1xoAmSZLUGAOaJElSYwxokiRJjTGgSZIkNcaAJkmS1BgDmiRJUmMMaJIkSY3pNaAlWZ9kT5K9SS6eZvlvJ9mVZEeSTyc5sc/+JEmSWtBbQEuyCLgSOAtYC5yfZO2Usu3Auqo6Bfg48Pt99SdJktSKPo+gnQ7srao7quph4Frg3OGCqrqxqr7fTX4OWNFjf5IkSU3oM6AtB+4cmt7XzTuc1wH/Y7oFSTYk2ZZk24EDB+awRUmSpIXX5E0CSS4A1gGXT7e8qjZW1bqqWrdkyZJ+m5MkSZpni3tc135g5dD0im7eYyR5KfC7wM9V1UM99SZJktSMPo+g3QKsSbI6yTHAecDm4YIkpwEfAM6pqnt67E2SJKkZvQW0qnoEuBDYCuwGNlXVziSXJTmnK7sceArwsSRfTLL5MMNJkiQdtfo8xUlVbQG2TJl36dDnl/bZjyRJUouavElAkiRpkhnQJEmSGmNAkyRJaowBTZIkqTEGNEmSpMYY0CRJkhpjQJMkSWqMAU2SJKkxBjRJkqTGGNAkSZIaY0CTJElqjAFNkiSpMQY0SZKkxhjQJEmSGmNAkyRJaowBTZIkqTEGNEmSpMYY0CRJkhpjQJMkSWqMAU2SJKkxvQa0JOuT7EmyN8nF0yx/UpLruuU3J1nVZ3+SJEkt6C2gJVkEXAmcBawFzk+ydkrZ64B7q+rvAe8B3tVXf5IkSa3o8wja6cDeqrqjqh4GrgXOnVJzLvDh7vPHgZckSY89SpIkLbhUVT8rSl4JrK+qf9pN/xrw/Kq6cKjm9q5mXzf91a7mW4cbd81xx9dfnHvOSD3suvu7AKxddtyc1s53/biObS/9j20v/Y9tL4ev/frTVvLaa/7TSGNLkyLJrVW1bqa6xX00M9eSbAA2dJMPrbr6I7cvZD+aU0uBwwZyjSX36dFlVvvzN679g3lsRXPAv5/9O3GUoj4D2n5g5dD0im7edDX7kiwGngp8e+pAVbUR2AiQZNsoSVTjwf159HGfHl3cn0cX92e7+rwG7RZgTZLVSY4BzgM2T6nZDLy6+/xK4Ibq6xysJElSI3o7glZVjyS5ENgKLAI+WFU7k1wGbKuqzcAfAx9Jshf4DoMQJ0mSNFF6vQatqrYAW6bMu3To84PAP5rlsBvnoDW1w/159HGfHl3cn0cX92ejeruLU5IkSaPxVU+SJEmNGeuANtOro9S2JB9Mck/3/LtD805I8qkkX+n+/MmF7FGjS7IyyY1JdiXZmeRN3Xz36RhK8uQkn0/ypW5/vqObv7p7Fd/e7tV8xyx0rxpdkkVJtif5ZDft/mzU2Aa0EV8dpbZdBayfMu9i4NNVtQb4dDet8fAIcFFVrQXOAH6r+zvpPh1PDwEvrqpTgecC65OcweAVfO/pXsl3L4NX9Gl8vAnYPTTt/mzU2AY0Rnt1lBpWVX/J4G7dYcOv+/ow8PJem9ITVlV3V9UXus/fY/A/geW4T8dSDfxNN7mk+yngxQxexQfuz7GSZAXwi8AfddPB/dmscQ5oy4E7h6b3dfM03p5eVXd3n78OPH0hm9ETk2QVcBpwM+7TsdWdDvsicA/wKeCrwH1V9UhX4u/d8XIF8DvAD7rpn8L92axxDmg6ynUPKfY24zGT5CnAfwX+RVV9d3iZ+3S8VNXBqnougze/nA48e4Fb0hOU5JeAe6rq1oXuRaMZy3dxdkZ5dZTGzzeSLKuqu5MsY/Avd42JJEsYhLP/UlWf6Ga7T8dcVd2X5EbgBcDxSRZ3R138vTs+Xgick+Rs4MnAccB/xP3ZrHE+gjbKq6M0foZf9/Vq4M8WsBfNQnc9yx8Du6vq3UOL3KdjKMnTkhzffT4W+AUG1xXeyOBVfOD+HBtVdUlVraiqVQz+f3lDVf0q7s9mjfWDart/CVzBo6+OeucCt6RZSHINcCawFPgG8DbgemAT8Ezga8CrqmrqjQRqUJIXAf8LuI1Hr3H5lwyuQ3OfjpkkpzC4aHwRg3/Mb6qqy5L8NIObsk4AtgMXVNVDC9epZivJmcCbq+qX3J/tGuuAJkmSdDQa51OckiRJRyUDmiRJUmMMaJIkSY0xoEmSJDXGgCZJktQYA5okSVJjDGiSJEmNMaBJ0pAkL0nykYXuQ9JkM6BJ0mOdyuCJ6pK0YAxokvRYpwLbkzwpyVVJ/m33nlFJ6s3ihW5AkhpzCnAPsBX4o6q6eoH7kTSBfBenJHWSLAG+xeCl7q+vqs8ucEuSJpSnOCXpUScBtwCPAAcXuBdJE8yAJkmPOhX4K+A84ENJnr7A/UiaUAY0SXrUqcDtVfVl4K3Apu60pyT1ymvQJEmSGuMRNEmSpMYY0CRJkhpjQJMkSWqMAU2SJKkxBjRJkqTGGNAkSZIaY0CTJElqjAFNkiSpMf8fHH0W6wzgnjIAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x112c8bda0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAl8AAADiCAYAAABwdKKfAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAF0NJREFUeJzt3X+0XWV54PHv4yWUixRSIXXIBQ22NpVVfmRWloDUjqN2hQpTMyxWi1atdhyYtcSxUw2STmfVtoPSpuOvOqvLiKItNg2lmYwjTKNTcI2/Bky4SATmKiKgNyjBepXWOxAuz/xxdvAm3tycfXPOe/Y+5/tZ66ycs8/+8Zzs7Pd9st93v29kJpIkSSrjGYMOQJIkaZSYfEmSJBVk8iVJklSQyZckSVJBJl+SJEkFmXxJkiQVZPIlqZaIeHFETA06jiaKiM9ExBuPYPt/jIjn9TImSc1j8iVpQRHxQES8/ODlmfnZzFw9iJiGyUKJWmYel5n3DyomSWWYfElqhYg4atAxSFIvmHxJqiUiXhIR35r3+YGIeFtE3BUR34+IrRFxzLzvL4qIOyNiJiK+EBFnzvvuqoj4ekQ8FhH3RMS/nvfd6yPi8xHxnoj4LvCOBWIZi4jfnbePXRFxavXdiyLiS1VMX4qIF83b7jMR8UfV/h+LiE9FxEnVd/8zIq446DhfjoiLD7ffg7Z5R0RcP+/zqojIiDgqIq4GXgx8oGpq/EC1TkbEz1bvT4iIv4iIvRHxYET8XkQ8Y97fzeci4k8j4nsR8Y2I+JVuzp+kwTP5ktQLvwZcAJwGnAm8HiAi1gAfAS4HTgQ+CHwiIn6i2u7rdJKQE4A/AK6PiJPn7fcc4H7g2cDVCxz3d4BXAa8Ajgd+C/hhRDwLuAl4f3XcdwM3RcSJ87Z9NfAG4KeBo4G3Vcu3VPuk+g2nA8+ttu9mv4eVmf8R+CxwRdXUeMUCq/0Znb+X5wH/AnhdFe9+5wBTwEnAnwAfjoioE4ekwTD5ktQL78/MPZn5D8D/AM6ull8GfDAzb8vMucz8GPA4cC5AZv5Ntd1TmbkV+Brwwnn73ZOZf5aZT2bm7ALHfSPwe5k5lR1fzszvAhcCX8vMv6y23QL8X+Bfzdv2usz8arXfG+bF/N+AsyPiudXn3wC2ZebjXe73iEXEGHApsDEzH8vMB4D/Arx23moPZuaHMnMO+BhwMp0kVVLDmXxJ6oVvz3v/Q+C46v1zgbdWTY4zETEDnAqsBIiI181rkpwBfoHOnZz9vnmY455K5+7ZwVYCDx607EFg4nAxZ+ZjdO5uXVp99yrg4zX22wsnAcsOOtYh48/MH1Zvj0NS45l8SeqnbwJXZ+byea9jM3NLdWfpQ8AVwImZuRz4CjC/6Sy72P/PLLB8D53Eb77nANNdxr0FeFVEnAccA9y6hP3+E3DsvM//7KDvF/ttjwL7DjpWnfglNZjJl6TFLIuIY+a96j5x+CHg30XEOdHxzIi4MCJ+EngmnQRkL0BEvIHOna86rgX+KCKeX+3/zKr/1c3Az0XEq6sO7r8OnA58ssv93kwn8flDYGtmPjVvebf7vRP4pYh4TkScAGw86Pvv0OnP9WOqpsQbgKsj4ierRPV3gOsXWl9Su5h8SVrMzcDsvNc76mycmTuBfwt8APgecB9VZ/zMvIdOP6Yv0klEzgA+XzO+d9NJUj4F/AD4MDBe9fu6CHgr8F3gSuCizHy0y7gfB7YBLwf+at7yrvebmZ8GtgJ3Abv48QTtfcAl1dOK718gjDfTuXt2P/C5Ko6PdBO/pGaLzMPd1ZckSVKveOdLkiSpIJMvSZKkgky+JEmSCjL5kiRJKsjkS5IkqaC6Y/YUddJJJ+WqVasGHYYkSdJh7dq169HMXHG49RqdfK1atYqdO3cOOgxJkqTDioiDpx9bkM2OkiRJBZl8SZIkFWTyJUmSVJDJlyRJUkEmX5IkSQWZfEmSJBVk8iVJklSQyZckSVJBJl+SJEkFFU2+IuI/RMTdEfGViNgSEceUPL4kSdKgFZteKCImgH8PnJ6ZsxFxA3Ap8NFSMUj9tH1ymk07ptgzM8vK5eNsWLea9WsmerJ+W/dtLKP9OyUtrPTcjkcB4xGxDzgW2FP4+FLX6lZgG7ftZnbfHADTM7Ns3LYbYMFt6qzf1n0by2j/zvnbmKxJByrW7JiZ08CfAg8BDwPfz8xPlTq+VMf+SmZ6ZpbkR5XM9snpBdfftGPq6Qppv9l9c2zaMXXE67d138Yy2r8T6l9H0qgolnxFxE8BrwROA1YCz4yI1yyw3mURsTMidu7du7dUeNIB6lYye2Zm+7a8rfs2ltH+nVD/OpJGRckO9y8HvpGZezNzH7ANeNHBK2Xm5sxcm5lrV6xYUTA8jYLtk9Ocf80tnHbVTZx/zS2H/B943Upm5fLxvi1v676NZbR/J9S/jrq9PqW2K5l8PQScGxHHRkQALwPuLXh8jbg6TSB1K5kN61YzvmzsgGXjy8bYsG71Ea/f1n0by2j/Tqh3HdlEqVFSss/XbcCNwB3A7urYm0sdX6rTBFK3klm/ZoJ3XXwGR491LqmJ5eO86+IzDtmxuM76bd23sYz274R615FNlBolkZmDjuGQ1q5dmzt37hx0GBoSp111Ewv9aw/gG9dc+GPLt09Oc+WNd/HE3FNMdPmU1q9/8IsAbL38vK5iqrN+W/dtLOX33aRYur2O6l6fUhNFxK7MXHu49UoPNSENzMrl40wv0NfkUE0j69dMsOX2h4DuKyVJB+r2Oqp7fUpt5vRCGhl1mxIlleP1qVHinS+1Wp0BHPcvr9uUKKn/lnJ9OoCr2srkS621lNG2bUqUmqvO9bmU619qCpsd1Vo+HSWNLq9/tZnJl1qr7gCOkoaH17/azORLrVV3IFRJw8PrX21m8qXW8ukoaXR5/avN7HCv1vLpRWl0ef2rzUy+1Go+vSiNLq9/tZXJlxrFcXsk9YNli5rE5EuN4bg9kvrBskVNY4d7NYbj9kjqB8sWNY3JlxrDcXsk9YNli5rG5EuN4bg9kvrBskVNY/KlxnDcHkn9YNmiprHDvRrDcXsk9YNli5rG5EuN4rg9kvrBskVNYrOjJElSQSZfkiRJBZl8SZIkFWSfL/Wd03pIahvLLfWTyZf6ymk9JLWN5Zb6zWZH9ZXTekhqG8st9ZvJl/rKaT0ktY3llvrN5Et95bQektrGckv9ZvKlvnJaD0ltY7mlfrPDvfrKaT0ktY3llvqtaPIVEcuBa4FfABL4rcz8YskYVJ7TekhqG8st9VPpO1/vA/4uMy+JiKOBYwsfX5IkaaCKJV8RcQLwS8DrATLzCeCJUseXJElqgpId7k8D9gLXRcRkRFwbEc8seHxJkqSBK5l8HQX8c+DPM3MN8E/AVQevFBGXRcTOiNi5d+/eguFJkiT1X8k+X98CvpWZt1Wfb2SB5CszNwObAdauXZvlwlO3nPNMkg5kuag6iiVfmfntiPhmRKzOzCngZcA9pY6v3nDOM0k6kOWi6io9yOqbgY9HxF3A2cA7Cx9fR8g5zyTpQJaLqqvoUBOZeSewtuQx1VvOeSZJB7JcVF1OL6RanPNMkg5kuai6TL5Ui3OeSdKBLBdVl3M7qhbnPJOkA1kuqi6TL9XmnGeSdCDLRdVhs6MkSVJBJl+SJEkFmXxJkiQVZPIlSZJUkMmXJElSQSZfkiRJBZl8SZIkFeQ4XwJg++Q0m3ZMsWdmlpUOEChJfWN5K5MvsX1ymo3bdjO7bw6A6ZlZNm7bDWCBIEk9ZHkrsNlRwKYdU08XBPvN7ptj046pAUUkScPJ8lZg8iVgz8xsreWSpKWxvBWYfAlYuXy81nJJ0tJY3gq6SL4i4lldvJaXCFb9sWHdasaXjR2wbHzZGBvWrR5QRJI0nCxvBd11uN9TvWKRdcaA5/QkIhW3v5PnlTfexRNzTzHh0zeS1BeWt4Lukq97M3PNYitExGSP4tGArF8zwZbbHwJg6+XnDTgaSRpelrfqps/XeQAR8Z8P/iIixuavI0mSpMUdNvnKzP9XvZ2IiFfvXx4RPw38r4PWkSRJ0iLqDLJ6ObAjIu4DErgOeHtfopIkSRpSh02+IuIvgDuASeBNwF8BTwLrM/O+/oYnSZI0XLrp8/VROk86vgG4HlgFfA94TURc0rfIJEmShtBh73xl5i3ALfs/R8RRwAuAs4BzgBv7Fp0kSdKQqT2xdmY+CeyuXtf3PCJJkqQh1s0I93f0Yh1JkiR1d+frBRFx1yLfB3BCj+KRJEkaat0kXz/fxTpz3R6wGph1JzCdmRd1u50kSdIw6KbD/YMAEfFp4G2Z+eUjPOZbgHuB449wP1rE9slpNu2YYs/MLCudO0ySWsvyfPh0M9TEfm8H3hsR10XEyUs5WEScAlwIXLuU7dWd7ZPTbNy2m+mZWRKYnpll47bdbJ+cHnRokqQaLM+HU9fJV2bekZn/Evgk8HcR8fsRMV7zeO8FrgSeqrmdati0Y4rZfQe2BM/um2PTjqkBRSRJWgrL8+FU584XERHAFPDnwJuBr0XEa7vc9iLgkczcdZj1LouInRGxc+/evXXCU2XPzGyt5ZKkZrI8H05dJ18R8XlgGngPMAG8HngJ8MKI2NzFLs4HfjUiHgD+GnhpRPzYOGGZuTkz12bm2hUrVnQbnuZZuXzhG5KHWi5JaibL8+FU587XZcBEZv5yZv6nzPxkZt6XmW8GXny4jTNzY2aekpmrgEuBWzLzNUsLW4vZsG4148vGDlg2vmyMDetWDygiSdJSWJ4Pp65HuM/Muxf5+sIexKIe2f8UzJU33sUTc08x4dMxktRKlufDqfb0QgvJzPtrrv8Z4DO9OLYWtn7NBFtufwiArZefN+BoJElLZXk+fGp1uJckSdKRMfmSJEkqyORLkiSpIJMvSZKkgky+JEmSCjL5kiRJKsjkS5IkqSCTL0mSpIJMviRJkgrqyQj36r/tk9Ns2jHFnplZVjq9hCTpEKwvms/kqwW2T06zcdtuZvfNATA9M8vGbbsBvKAkSU+zvmgHmx1bYNOOqacvpP1m982xacfUgCKSJDWR9UU7mHy1wJ6Z2VrLJUmjyfqiHUy+WmDl8vFayyVJo8n6oh1Mvlpgw7rVjC8bO2DZ+LIxNqxbPaCIJElNZH3RDna4b4H9nSSvvPEunph7igmfXpEkLcD6oh1Mvlpi/ZoJttz+EABbLz9vwNFIkprK+qL5bHaUJEkqyORLkiSpIJMvSZKkgky+JEmSCjL5kiRJKsjkS5IkqSCTL0mSpIIc52uAtk9Os2nHFHtmZlnpQHiSpMKshwbD5GtAtk9Os3Hb7qdnn5+emWXjtt0A/sOXJPWd9dDg2Ow4IJt2TD39D36/2X1zbNoxNaCIJEmjxHpocEy+BmTPzGyt5ZIk9ZL10OAUS74i4tSIuDUi7omIuyPiLaWO3UQrl4/XWi5JUi9ZDw1OyTtfTwJvzczTgXOBN0XE6QWP3ygb1q1mfNnYAcvGl42xYd3qAUUkSRol1kODU6zDfWY+DDxcvX8sIu4FJoB7SsXQJPs7M1554108MfcUEz5lIkkqyHpocAbytGNErALWALcN4vhNsX7NBFtufwiArZefN+BoJEmjxnpoMIp3uI+I44C/BX47M3+wwPeXRcTOiNi5d+/e0uFJkiT1VdHkKyKW0Um8Pp6Z2xZaJzM3Z+bazFy7YsWKkuFJkiT1XcmnHQP4MHBvZr671HElSZKapOSdr/OB1wIvjYg7q9crCh5fkiRp4Eo+7fg5IEodT5IkqYmc27GHnKBUkjTMrOd6w+SrR5ygVJI0zKznese5HXvECUolScPMeq53TL56xAlKJUnDzHqud0y+esQJSiVJw8x6rndMvnrECUolScPMeq537HDfI05QKkkaZtZzvWPy1UNOUCpJGmbWc71hs6MkSVJBJl+SJEkF2ey4CEfylSRpaaxDD83k6xAcyVeSpKWxDl2czY6H4Ei+kiQtjXXo4ky+DsGRfCVJWhrr0MWZfB2CI/lKkrQ01qGLM/k6BEfylSRpaaxDF2eH+0NwJF9JkpbGOnRxJl+LcCRfSZKWxjr00EYu+XLcEUmSmmeU6ueRSr4cd0SSpOYZtfp5pDrcO+6IJEnNM2r180glX447IklS84xa/TxSyZfjjkiS1DyjVj+PVPLluCOSJDXPqNXPre9wX+fpCMcdkSSpeZZSP7f56chWJ19LeTrCcUckSWqeOvVz25+ObHWz46g9HSFJktpf/7c6+Rq1pyMkSVL76/+iyVdEXBARUxFxX0Rcdbj1d09/n/OvuYXtk9MLfj9qT0dIkqT69f/2yWnOv+YWTrvqpkXzilKKJV8RMQb8V+BXgNOBV0XE6Yfbbn877kJ/UaP2dIQkSapX/+/vHzY9M0uyeF5RSsk7Xy8E7svM+zPzCeCvgVd2s+Gh2nHXr5ngXRefwdFjnZ8xsXycd118Ris620mSpKWpU/83sX9YZGaZA0VcAlyQmW+sPr8WOCczrzjUNj9z3An5zrN+8enP5z7vxAXXu+fhHwBw+snHdxVLnfX7uW9jKb9vYym/b2Mpv29jKb9vYym/727X/z/3f/fp9/efMMEHz+zc9wngG9dc2NVxuhURuzJz7eHWa9xQExFxGXAZwDPGj+d1938VgJx78ol9X3hg9yBj0xE7CXh00EGoZzyfw8XzOXw8p8CyFavOiLGjju58+ipM3gp08or444t6nVc8t5uVSiZf08Cp8z6fUi07QGZuBjYDRMTOx3/4/cNmkGqHiNjZzf8I1A6ez+Hi+Rw+ntPmKtnn60vA8yPitIg4GrgU+ETB40uSJA1csTtfmflkRFwB7ADGgI9k5t2lji9JktQERft8ZebNwM01Ntncr1g0EJ7P4eL5HC6ez+HjOW2oYk87SpIkqeXTC0mSJLVNI5OvutMQqXki4iMR8UhEfGXesmdFxKcj4mvVnz81yBjVvYg4NSJujYh7IuLuiHhLtdxz2kIRcUxE3B4RX67O5x9Uy0+LiNuqsndr9XCUWiIixiJiMiI+WX32fDZU45KvpU5DpMb5KHDBQcuuAv4+M58P/H31We3wJPDWzDwdOBd4U3Vdek7b6XHgpZl5FnA2cEFEnAv8MfCezPxZ4HvAvxlgjKrvLcC98z57PhuqcckXRzANkZojM/838A8HLX4l8LHq/ceA9UWD0pJl5sOZeUf1/jE6BfwEntNWyo5/rD4uq14JvBS4sVru+WyRiDgFuBC4tvoceD4bq4nJ1wTwzXmfv1UtU/s9OzMfrt5/G3j2IIPR0kTEKmANcBue09aqmqjuBB4BPg18HZjJzCerVSx72+W9wJXAU9XnE/F8NlYTky+NgOw8Zuujti0TEccBfwv8dmb+YP53ntN2ycy5zDybzmwjLwR+fsAhaYki4iLgkczcNehY1J3Gze1Il9MQqZW+ExEnZ+bDEXEynf9xqyUiYhmdxOvjmbmtWuw5bbnMnImIW4HzgOURcVR1t8Sytz3OB341Il4BHAMcD7wPz2djNfHOl9MQDa9PAL9Zvf9N4L8PMBbVUPUf+TBwb2a+e95XntMWiogVEbG8ej8O/DKdfny3ApdUq3k+WyIzN2bmKZm5ik6deUtm/gaez8Zq5CCrVfb+Xn40DdHVAw5JNUXEFuAlwEnAd4DfB7YDNwDPAR4Efi0zD+6UrwaKiF8EPgvs5kd9Sn6XTr8vz2nLRMSZdDpgj9H5T/gNmfmHEfE8Og85PQuYBF6TmY8PLlLVFREvAd6WmRd5PpurkcmXJEnSsGpis6MkSdLQMvmSJEkqyORLkiSpIJMvSZKkgky+JEmSCjL5kiRJKsjkS5IkqSCTL0kjIyJeFhF/Oeg4JI02ky9Jo+QsOiN9S9LAmHxJGiVnAZMR8RMR8dGIeGc1b6UkFXPUoAOQpILOBB4BdgDXZub1A45H0ghybkdJIyEilgGP0pkA/PLM/OKAQ5I0omx2lDQqXgB8CXgSmBtwLJJGmMmXpFFxFvAF4FLguoh49oDjkTSiTL4kjYqzgK9k5leBtwM3VE2RklSUfb4kSZIK8s6XJElSQSZfkiRJBZl8SZIkFWTyJUmSVJDJlyRJUkEmX5IkSQWZfEmSJBVk8iVJklTQ/wezTfCxznDm1QAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x112d12898>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAm8AAADiCAYAAAAYhyCnAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAG6FJREFUeJzt3XmYZXV95/H3J92FEhEBGxGazShBIWOj6cHdtDFBJAhmhhjQuEWnNdE8cUajxmQAjcnjxCiOK+nIoqC4REGiILRiABMXmn0XZCD0Ai2boKJ243f+OKfjpfpW9+2m6t57qt6v57lPnfM7v/O733tPQX36rKkqJEmS1A2/MuoCJEmSNDjDmyRJUocY3iRJkjrE8CZJktQhhjdJkqQOMbxJkiR1iOFN6oAkL0ty7laue2ySU9vpPZP8KMm86a1wuJKcnOTdD2H9q5MsmcaSpvX9k/xrktcOsaTe935Okuu3ct1XJfnmdNck6cEMb9IMSXJzkvvbsHR7Gzi225qxqupTVXXQQ62pqv6jqrarqgce6lhd0S/oVdX+VfWvIyrpQe/fG67HQVVdWFX7Tve4SfZOUkkundS+IMnPk9w8A++5LMn1SX6R5FV9lv9aki8nuS/JHUn+frprkGaC4U2aWS+qqu2ApwKLgb/e0gGSzJ/2qqQ+hvS79qtJfqNn/qXA/5uh97oc+FPgkskLkmwDLAfOAx4L7A6MTYiWNsXwJg1BVa0CzgZ+AyDJo5KckGRNklVJ3r3hUGZ76OnfkhyX5E7g2MmHo5I8M8lFSX7Y/nxmz7LHJTm/3ZuwHFjQs2zD3o/57fxOSU5KsjrJ3UnOmOozJPkfSa5tx70myVPb9ie1h/nuaQ8HHtazzslJPpLkK+1630ny+HbZx5L8w6T3+FKS/7W5cSets9GhuvYzPiHJUuBlwFvbPaD/0i6/OcnvtNMPS/KB9jtY3U4/rF22JMnKJG9OsrbdXq+eoo7nJbmyZ355kot65i9M8uLe909yMPAO4A/b+i7vGXKv9vfgviTnJllAHz01vqPde3Rzkpf1LH9Ykn9I8h/tHuDjk2w7ad23JbkNOGlDW8/6m9q+j05yZpJ7k3wXeHy/Gic5BXhlz/wrgE8OsN4Wq6qPVNXXgZ/2WfwqYHVVvb+qflxVP62qK2aiDmm6Gd6kIUiyB3AIsOGQ0cnAeuAJwFOAg4Dec5yeBtwE7AL87aSxdgK+AnwQeDTwfuArSR7ddvk0cDFNaPsbHvyHcrJTgF8F9gceAxw3Rf1/ABxL84d2e+Aw4M4kE8C/AOe26/8Z8KkkvYfdjgTeCewI3NjzeU6jCS1p32PH9nv4zIDjblZVLQM+Bfx9e7j4RX26/RXwdOAAYBFwIA/eQ/pY4FHAQuA1wEfaWif7NrBPmsOAE8CTgd2SPLINS4uBCyfV91Xg74DPtvUt6ln8UuDV7effBnjLJj7qY2m290Ka7b2s57t6D/Dr7ed7Qtvn6Enr7gTsBSztHXSA7fARmmC0K/DH7WtzTgWOTDIvyX7AdsB3NrVCkiva8Njv9dEB3rOfpwM3Jzm7Db3/muS/bOVY0lAZ3qSZdUaSe4BvAucDf5dkF5og96b2X/xraULTkT3rra6qD1XV+qq6f9KYvwfcUFWntMtPA64DXpRkT+C/Av+7qn5WVRfQ/PHdSJJdgRcCr6+qu6tqXVWdP8XneC1NALqoGjdW1S00fwC3A95TVT+vqvOALwNH9ax7elV9t6rW0wSpA9r2C4ECntPOHwF8q6pWDzjudHkZ8K6qWltVP6AJmi/vWb6uXb6uqs4CfgRsFCLb7XQR8FzgN2kO2f0b8Kz289xQVXduQV0nVdX32nE/xy+/t6ls2Obn04T7l7TBeCnwP6vqrqq6jyYs9v6u/QI4pl138u/alNshzZ7i/w4c3f4eXwV8YoDPtRK4Hvgdmn8MnLK5FarqyVW1wxSvPx3gPfvZneZ7+CCwG8139qU0h1Olsea5NNLMenFVfa23of3X/QSwpt3pBM0/pG7t6dY7PdluwC2T2m6h2aOyG3B3Vf140rI9+oyzB3BXVd29uQ/R9v3+FLXcWlW/6FPLBrf1TP+EJgxQVZXkMzSB7AKaPU0bzjkaZNzpMvn7vKVt2+DONnhu8J+foY/zgSU0AeV84G7gt4CftfNbou/3NoV+23w3YGeaPasX9/yuBei92vgHVdXvsCJsejvsTPM35NZJywbxSZrDls+kCe+/PuB60+l+4JtVdTZAewj/r4En0QRvaWy5500avltp/pgv6Nl7sH1V7d/Tpzax/mqaQ1y99gRWAWuAHZM8YtKyqerYKckOA9bc73ym1cAeSXr/X7KhlkGcBhyRZC+aQ8Vf2Ipxf0wTUABI8thJyzf1XW54r97vc8+2bWtsCG/PbafPpwlvv8XU4W1z9Q2i3zZfDdxBE1L27/lde1R7Ec0g77+p7fADmkP/e0xaNogv0OxBvqmq/mNzndtz7X40xev4Ad9zsiuYnu9eGjrDmzRkVbWG5hyi9yXZPsmvJHl8kt8acIizgF9P8tIk85P8IbAf8OX2UOYK4J1JtknybKDfeV4b6jgb+GiSHZNMJHnuFO/5ceAtSX4zjSe0ges7NHuF3tquv6R9v88M+F1cShMwPg6cU1X3tIu2ZNzLgf2THJDk4TTn5vW6Hfi1TZRxGvDXSXZuLwo4mq2/6vDfaQ6pHgh8t6qupgmGT6PZu9jP7cDekwLS1tiwzZ8DHAp8vt1j9k/AcUkeA5BkYZIXDDjmlNuhvd3MF2kuqPnV9vy1TZ1f+Z/avYS/zYPP89xU//3bcwL7vV4/1Xrt9/Fwmr2NE0ke3vM9nwo8Pc2FI/OAN9H8Ll47SE3SKBnepNF4Bc1J6NfQHFr7Z5qTvjerPW/qUODNwJ3AW4FDq+qOtstLacLCXcAxbPpKvpfTnNN1HbCW5g9Yv/f8PM2FBp8G7gPOAHaqqp/T/DF/Ic0fvo8Cr6iq6wb5LK1P05z/9Ome9xt43Kr6HvAu4GvADTTnF/Y6AdivPbm939W076YJvFcAV9LcVmKrbgDchpJLgKvbzwDwLeCW9tzGfj7f/rwzyUa3tBjQbTS/R6tpzit8fc939TaaC0W+neRemu9poAs/BtgOb6Q5nHsbzUU4Jw1acFWtqKp+h+Kn07k0ex6fCSxrp5/bvv/1wB8Bx9N8d4cDh/VsN2lspcq9xpLUVe3esFOravdR1yJpONzzJkmS1CGGN0mSpA4Z2mHTJCfSnKeztqo23GX+s/zy3IsdgHuqaqN7GaV55t19wAPA+qpaPJSiJUmSxswww9tzaW5u+ckN4W3S8vcBP6yqd/VZdjOwuOeEbEmSpDlpaDfpraoLkuzdb1l7F/CX0Fw6LkmSpCmMyxMWngPcXlU3TLG8gHOTFPCP1TyvcLMWLFhQe++99zSVKEmSNHMuvvjiO6pq5831G5fwdhTNjTKn8uyqWtXeZHJ5kuvaZzZuJMlS2ocr77nnnqxYsWL6q5UkSZpmSQZ6xNzIrzZNMh/4b8Bnp+pTVavan2uB02nuXj5V32VVtbiqFu+882bDqyRJUqeMPLzR3Fn9uqpa2W9hkkckeeSGaeAg4Koh1idJkjQ2hhbekpxG85iYfZOsTPKadtGRTDpkmmS3JGe1s7sA30xyOfBd4CtV9dVh1S1JkjROhnm16VFTtL+qT9tq4JB2+iZg0YwWJ0mS1BHjcNhUkiRJAzK8SZIkdYjhTZIkqUMMb5IkSR1ieJMkSeoQw5skSVKHjMvjsWbElat+yLPecx5/8YJ9efFTFvbtc8alq3jvOdez+p772W2HbTfZd0v7z+TY1jK3P6ckae6ad+yxx466hhnz7vd96Nh64u9w/vd+wO47bssTd93+QcvPuHQVf/nFK7nrJz8H4L6frp+y75b2n8mxrWVuf05J0uz0zne+c82xxx67bHP95sRh0/vXPcB7z7l+o/b3nnM99697YKC+W9p/Jse2lrn9OSVJc9ucCG8Aq++5f6C26WqfybGtZW5/TknS3DZnwttuO2w7UNt0tc/k2NYytz+nJGlumxPhbduJefzFC/bdqP0vXrAv207MG6jvlvafybGtZW5/TknS3DbrL1h44pLf5+gX7df3qr0n7ro9u++4Leddt5YHqli4w7ZT9t3S/jM5trXM7c8pSZqdBr1gIVU1jHpGYqe9nlR33XLtZvv94T9+C4DPvu4ZA427Jf1ncmxrGf7Y41aLJGn2SHJxVS3eXL+hHTZNcmKStUmu6mk7NsmqJJe1r0OmWPfgJNcnuTHJ24dVsyRJ0rgZ5jlvJwMH92k/rqoOaF9nTV6YZB7wEeCFwH7AUUn2m9FKJUmSxtTQwltVXQDctRWrHgjcWFU3VdXPgc8Ah09rcZIkSR0xDlebvjHJFe1h1R37LF8I3Nozv7JtkyRJmnNGHd4+BjweOABYA7zvoQ6YZGmSFUlWrFu37qEOJ0mSNFZGGt6q6vaqeqCqfgH8E80h0slWAXv0zO/etk015rKqWlxViycmJqa3YEmSpBEbaXhLsmvP7O8DV/XpdhGwT5LHJdkGOBI4cxj1SZIkjZv5w3qjJKcBS4AFSVYCxwBLkhwAFHAz8Lq2727Ax6vqkKpan+SNwDnAPODEqrp6WHVLkiSNk6GFt6o6qk/zCVP0XQ0c0jN/FrDRbUQkSZLmmlFfsCBJkqQtYHiTJEnqEMObJElShxjeJEmSOsTwJkmS1CGGN0mSpA4xvEmSJHWI4U2SJKlDDG+SJEkdYniTJEnqEMObJElShxjeJEmSOsTwJkmS1CGGN0mSpA4ZWnhLcmKStUmu6ml7b5LrklyR5PQkO0yx7s1JrkxyWZIVw6pZkiRp3Axzz9vJwMGT2pYDv1FVTwa+B/zlJtZ/XlUdUFWLZ6g+SZKksTe08FZVFwB3TWo7t6rWt7PfBnYfVj2SJEldNE7nvP0xcPYUywo4N8nFSZYOsSZJkqSxMn/UBQAk+StgPfCpKbo8u6pWJXkMsDzJde2evH5jLQWWAmy36+NnpF5JkqRRGfmetySvAg4FXlZV1a9PVa1qf64FTgcOnGq8qlpWVYuravHExMQMVCxJkjQ6Iw1vSQ4G3gocVlU/maLPI5I8csM0cBBwVb++kiRJs90wbxVyGvAtYN8kK5O8Bvgw8EiaQ6GXJTm+7btbkrPaVXcBvpnkcuC7wFeq6qvDqluSJGmcDO2ct6o6qk/zCVP0XQ0c0k7fBCyawdIkSZI6Y+TnvEmSJGlwhjdJkqQOMbxJkiR1iOFNkiSpQwxvkiRJHWJ4kyRJ6hDDmyRJUocY3iRJkjrE8CZJktQhhjdJkqQOMbxJkiR1iOFNkiSpQwxvkiRJHWJ4kyRJ6hDDmyRJUocMNbwlOTHJ2iRX9bTtlGR5khvanztOse4r2z43JHnl8KqWJEkaH8Pe83YycPCktrcDX6+qfYCvt/MPkmQn4BjgacCBwDFThTxJkqTZbKjhraouAO6a1Hw48Il2+hPAi/us+gJgeVXdVVV3A8vZOARKkiTNeuNwztsuVbWmnb4N2KVPn4XArT3zK9u2jSRZmmRFkhXr1q2b3kolSZJGbBzC23+qqgLqIY6xrKoWV9XiiYmJaapMkiRpPIxDeLs9ya4A7c+1ffqsAvbomd+9bZMkSZpTNhve2qtBN/fa4SHUcCaw4erRVwJf6tPnHOCgJDu2Fyoc1LZJkiTNKfMH6LO6fWUTfeYBe25uoCSnAUuABUlW0lxB+h7gc0leA9wCvKTtuxh4fVW9tqruSvI3wEXtUO+qqskXPkiSJM16g4S3a6vqKZvqkOTSQd6sqo6aYtHz+/RdAby2Z/5E4MRB3keSJGm2GuSct2cAJHn35AVJ5vX2kSRJ0szabHirqp+2kwuTvHRDe5LHAF+b1EeSJEkzaJDDphu8DjgnyY00t/M4CXjbjFQlSZKkvjYb3pJ8ErgEuBR4A/BpYD3w4qq6cWbLkyRJUq9Bznk7meZK01cDpwJ7A3cDf5TkiBmrTJIkSRvZ7J63qjoPOG/DfJL5wJOARTQPiv/nGatOkiRJD7Il57wBUFXrgSvb16nTXpEkSZKmNMgTFi6Zjj6SJEl66AbZ8/akJFdsYnmAR01TPZIkSdqEQcLbEwfo88BDLUSSJEmbN8gFC7cAJFkOvKWqLp/xqiRJktTXILcK2eBtwAeSnJRk15kqSJIkSVMbOLxV1SVV9Tzgy8BXkxyTZNuZK02SJEmTbcmeN5IEuB74GPBnwA1JXj4ThUmSJGljA4e3JP8GrAKOAxYCrwKWAAcmWba1BSTZN8llPa97k7xpUp8lSX7Y0+forX0/SZKkLtuSm/QuBa6pqprU/mdJrt3aAqrqeuAAgCTzaALi6X26XlhVh27t+0iSJM0GA4e3qrp6E4t/bxpqAXg+8P0NV7hKkiTpwbbonLepVNVN0zEOcCRw2hTLnpHk8iRnJ9l/qgGSLE2yIsmKdevWTVNZkiRJ42Fawtt0SLINcBjw+T6LLwH2qqpFwIeAM6Yap6qWVdXiqlo8MTExM8VKkiSNyNiEN+CFwCVVdfvkBVV1b1X9qJ0+C5hIsmDYBUqSJI3aOIW3o5jikGmSx7a3KSHJgTR13znE2iRJksbCllxtOmOSPAL4XeB1PW2vB6iq44EjgD9Jsh64Hziyz1WvkiRJs95YhLeq+jHw6Eltx/dMfxj48LDrkiRJGjfjdNhUkiRJm2F4kyRJ6hDDmyRJUocY3iRJkjrE8CZJktQhhjdJkqQOMbxJkiR1iOFNkiSpQwxvkiRJHWJ4kyRJ6hDDmyRJUocY3iRJkjrE8CZJktQhhjdJkqQOGZvwluTmJFcmuSzJij7Lk+SDSW5MckWSp46iTkmSpFGaP+oCJnleVd0xxbIXAvu0r6cBH2t/SpIkzRljs+dtAIcDn6zGt4Edkuw66qIkSZKGaZzCWwHnJrk4ydI+yxcCt/bMr2zbHiTJ0iQrkqxYt27dDJUqSZI0GuMU3p5dVU+lOTz6hiTP3ZpBqmpZVS2uqsUTExPTW6EkSdKIjU14q6pV7c+1wOnAgZO6rAL26JnfvW2TJEmaM8YivCV5RJJHbpgGDgKumtTtTOAV7VWnTwd+WFVrhlyqJEnSSI3L1aa7AKcngaamT1fVV5O8HqCqjgfOAg4BbgR+Arx6RLVKkiSNzFiEt6q6CVjUp/34nukC3jDMuiRJksbNWBw2lSRJ0mAMb5IkSR1ieJMkSeoQw5skSVKHGN4kSZI6xPAmSZLUIYY3SZKkDjG8SZIkdYjhTZIkqUMMb5IkSR1ieJMkSeoQw5skSVKHGN4kSZI6xPAmSZLUISMPb0n2SPKNJNckuTrJn/fpsyTJD5Nc1r6OHkWtkiRJozZ/1AUA64E3V9UlSR4JXJxkeVVdM6nfhVV16AjqkyRJGhsj3/NWVWuq6pJ2+j7gWmDhaKuSJEkaTyMPb72S7A08BfhOn8XPSHJ5krOT7D/UwiRJksbE2IS3JNsBXwDeVFX3Tlp8CbBXVS0CPgScsYlxliZZkWTFunXrZq5gSZKkERiL8JZkgia4faqqvjh5eVXdW1U/aqfPAiaSLOg3VlUtq6rFVbV4YmJiRuuWJEkatpGHtyQBTgCurar3T9HnsW0/khxIU/edw6tSkiRpPIzD1abPAl4OXJnksrbtHcCeAFV1PHAE8CdJ1gP3A0dWVY2iWEmSpFEaeXirqm8C2UyfDwMfHk5FkiRJ42vkh00lSZI0OMObJElShxjeJEmSOsTwJkmS1CGGN0mSpA4xvEmSJHWI4U2SJKlDDG+SJEkdYniTJEnqEMObJElShxjeJEmSOsTwJkmS1CGGN0mSpA4xvEmSJHXIWIS3JAcnuT7JjUne3mf5w5J8tl3+nSR7D79KSZKk0Rt5eEsyD/gI8EJgP+CoJPtN6vYa4O6qegJwHPB/hlulJEnSeBh5eAMOBG6sqpuq6ufAZ4DDJ/U5HPhEO/3PwPOTZIg1SpIkjYVU1WgLSI4ADq6q17bzLweeVlVv7OlzVdtnZTv//bbPHZsae5/td6ivHX7YZmu4Zs29AOy36/YD1bwl/WdybGsZ/tjDqOW2nffg1ad9aKD+kqTZI8nFVbV4c/3mD6OYYUqyFFjazv5s71NPuWqU9WhaLQA2Gdhniz/+zIdHXcIwzJntOUe4PWcft+nw7TVIp3EIb6uAPXrmd2/b+vVZmWQ+8Cjgzn6DVdUyYBlAkhWDJFh1g9tzdnF7zi5uz9nHbTq+xuGct4uAfZI8Lsk2wJHAmZP6nAm8sp0+AjivRn28V5IkaQRGvuetqtYneSNwDjAPOLGqrk7yLmBFVZ0JnACckuRG4C6agCdJkjTnjDy8AVTVWcBZk9qO7pn+KfAHWzH0sodYmsaL23N2cXvOLm7P2cdtOqZGfrWpJEmSBjcO57xJkiRpQLMyvG3ucVsaf0lOTLK2vcffhradkixPckP7c8dR1qjBJdkjyTeSXJPk6iR/3ra7TTsoycOTfDfJ5e32fGfb/rj2EYY3to803GbUtWpwSeYluTTJl9t5t+eYmnXhbcDHbWn8nQwcPKnt7cDXq2of4OvtvLphPfDmqtoPeDrwhva/S7dpN/0M+O2qWgQcAByc5Ok0jy48rn2U4d00jzZUd/w5cG3PvNtzTM268MZgj9vSmKuqC2iuLO7V+5i0TwAvHmpR2mpVtaaqLmmn76P5A7EQt2knVeNH7exE+yrgt2keYQhuz05Jsjvwe8DH2/ng9hxbszG8LQRu7Zlf2bap+3apqjXt9G3ALqMsRlsnyd7AU4Dv4DbtrPYQ22XAWmA58H3gnqpa33bx/73d8gHgrcAv2vlH4/YcW7MxvGkOaG/S7KXSHZNkO+ALwJuq6t7eZW7TbqmqB6rqAJqn4hwIPHHEJWkrJTkUWFtVF4+6Fg1mLO7zNs0GedyWuun2JLtW1Zoku9L8i18dkWSCJrh9qqq+2Da7TTuuqu5J8g3gGcAOSea3e2v8f293PAs4LMkhwMOB7YH/i9tzbM3GPW+DPG5L3dT7mLRXAl8aYS3aAu35MycA11bV+3sWuU07KMnOSXZop7cFfpfmPMZv0DzCENyenVFVf1lVu1fV3jR/M8+rqpfh9hxbs/Imve2/Hj7ALx+39bcjLklbKMlpwBJgAXA7cAxwBvA5YE/gFuAlVTX5ogaNoSTPBi4EruSX59S8g+a8N7dpxyR5Ms0J7PNodgJ8rqreleTXaC4S2wm4FPijqvrZ6CrVlkqyBHhLVR3q9hxfszK8SZIkzVaz8bCpJEnSrGV4kyRJ6hDDmyRJUocY3iRJkjrE8CZJktQhhjdJkqQOMbxJkiR1iOFNkgaU5PlJThl1HZLmNsObJA1uEc2d5iVpZAxvkjS4RcClSR6W5OQkf9c+t1WShmb+qAuQpA55MrAWOAf4eFWdOuJ6JM1BPttUkgaQZAK4A7gFeF1VfWvEJUmaozxsKkmDeRJwEbAeeGDEtUiawwxvkjSYRcC/A0cCJyXZZcT1SJqjDG+SNJhFwFVV9T3gbcDn2kOpkjRUnvMmSZLUIe55kyRJ6hDDmyRJUocY3iRJkjrE8CZJktQhhjdJkqQOMbxJkiR1iOFNkiSpQwxvkiRJHfL/Afb4vFiB/QN0AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x112d86cf8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAl8AAADiCAYAAABwdKKfAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAHCRJREFUeJzt3XuYXHWZ4PHvSxOkASECUUmDBm9RdgAzm1URdRjUJwiMMi47oqLiZcBnhHVmNEDU8TIOFyeOqOusY0RAB4lgzGZdb/ECrNcBA41EYNtB7h0uQWlgtDUhvPvHOU2q2+ruqk7Xqaqu7+d5+knVufzOW+f2e3PO7/xOZCaSJEmqxk7tDkCSJKmXmHxJkiRVyORLkiSpQiZfkiRJFTL5kiRJqpDJlyRJUoVMvqQJIuLFETHU7jg6UURcGRFv24H5/yMinjabMZXl3hYRL5vtcnfUjuxLEXFSRPyw5ntL1l2VIuKDEXHxDsz/zYh402zGJLWDyZd61mQVdmb+IDMXtyOmuaReopaZe2TmLe2KqWqzuS/12rqrl6hl5isy8/PtikmaLSZfUoeIiJ3bHYNmj9tT0mRMvqQJIuKIiLir5vttEfHuiLg+Ih6MiEsjYtea8cdGxHURMRIRP46IQ2rGnRkRv4yIhyPixoj485pxJ0XEjyLivIj4FfDBOrH0RcR7asq4JiIOKMe9MCJ+Wsb004h4Yc18V0bEh8vyH46Ib0fEvuW4b0bEqROW87OIePV05U6YZ9yViYhYFBEZETtHxFnAi4FPlbfLPlVOkxHxjPLzXhHxhYjYHBG3R8T7ImKnmnXzw4j4aEQ8EBG3RsQrptl0/6Vcxw9ExIVj2ygifh4Rf1YT57yIuD8iltT5TUdExF3lOr+/3Pavrxn/uDKmOyLi3oj4l4jonzDvGRFxD3BhnX3pOeW2GYmIGyLilTXj9omIr0bEQxFxNfD0CbHVrrv+iPincr09WK6r/km206vK/fOhcj86qhy+sFzeryPi5oj4ywnb9rJy+zxcxrq0HHdGRKyZsIxPRMQnpyu33rqeMOy2iHhZGeN7gNeU+8/PyvGPXU2NiJ3Kfeb2iLivjHWvctzYvvimclvdHxHvrReH1A4mX1Jj/gI4CjgQOAQ4CaCswC8ATgH2AT4DfDUiHlfO90uKJGQv4EPAxRGxX025zwduAZ4EnFVnuX8LvBY4GtgTeAvw24jYG/g68MlyuR8Dvh4R+9TM+zrgzcATgV2Ad5fDV5dlUv6Gg4CnlvM3Uu60MvO9wA+AU8vbZafWmex/UKyXpwF/AryxjHfM84EhYF/gH4HPRURMsdjXA8sokpZnAe8rh38BOLFmuqOBuzNzcJJynlwucwB4E7AqIsZuHZ5blv1c4BnlNO+fMO/eFOvz5NpCI2Ie8H+Ab1Nsk9OAL9aU/c/A74D9KLbzW6b4rR8F/jPwwnJ5pwOPTpwoIp5X/v7lwHzgJcBt5egvAXcBC4HjgbMj4sia2V9ZTjMf+CrwqZr5jo6Ix5fL6KM4Pi5psNxpZea3gLOBS8v959A6k51U/v0pxT60R02MY14ELAZeCrw/Ip7TTBxSy2Smf/715B9FJfSyOsOPAO6aMN2JNd//EfiX8vOngQ9PmH8I+JNJlnkd8Kry80nAHdPEODQ2/YThbwCunjDsJ8BJ5ecrgffVjPsr4Fvl58cDvwGeWn4/C7igiXLfVn7+IHBxzXSLgAR2njhtzTRJkbT0AVuAg2rGnQJcWbNubq4Zt1s575On2JZvr/l+NPDL8vNC4GFgz/L7GuD0Sco5AngE2L1m2GXA3wFRrren14w7DLi1Zt4twK719iWKJPweYKea8avL9dgHbAWeXTPubOCHddbdTsAocGgD+/hngPPqDD8A2AY8vmbYOcBFNdv2uzXjDgJGa77/EHhj+fnlNeu6kXIvrnecTTwmmbB/1dn/vgf8Vc24xeU63Jnt++L+NeOvBk6Ybp35518Vf175khpzT83n31L8LxuKKxzvKm8jjUTECEUFtBAgIt4Y229JjgB/RHFVZcyd0yz3AIqrZxMtBG6fMOx2iisxU8acmQ9TXN06oRz3WuCLTZQ7G/YF5k1Y1qTxZ+Zvy497MLnadXk75TbIzE3Aj4D/GhHzgVew/ffW80Bm/qZOWQsoksBrarbnt8rhYzZn5u8mKXchcGdm1l6hGvvNCyiShom/oZ59gV2pv19MNNX+8+tyX5gYy5iJ+8+usb0d2yVsv3r6OrZf9Wqk3NkycV+9nWIdPqlm2GTHrdRWJl/SjrkTOCsz59f87ZaZqyPiqcBngVOBfTJzPvBziisoY7KB8p9eZ/gmisSv1lOA4QbjXg28NiIOo6jIr5hBub+hSEbGPHnC+Kl+2/0UVylql9VM/PUcMKGsTTXfP09x6/G/AT/JzKmW84SI2L1OWfdTXHH6TzXbeq/MrK3Qp/rNm4ADomzXVlP2MLCZ4orbxN9Qz/0Utyfr7RcTTbX/7D1263BCLI34MnBEROwP/Dnbk69myh23/5S3L2sT2emOjYn76lMo1uG9jfwAqZ1MvtTr5kXErjV/zT6h9lng7RHx/CjsHhHHlJXP7hQVyGaAiHgzxZWvZpwPfDginlmWf0jZ/uobwLMi4nVRNHB/DcWtoa81WO43KCquv6doV/NozfBGy70OeElEPKVs6Lxiwvh7Kdri/IHM3EZxO++siHh8maj+LTDjPqCAd0TE/mW7tfcCl9aMWwf8MfBOijZQ0/lQROwSES8GjgW+XK6jzwLnRcQTASJiICKWNRjfVRRXX06PotH/EcCfAV8q18da4IMRsVvZDq9uf1ZlHBcAHysbt/dFxGE17QxrfQ54c0S8tGygPhARz87MO4EfA+eU+/0hwFtpcP1n5maKW4AXUtx2vakc3ky5v6C4mnZM2R7ufUDtb7gXWDQhWa21GvibiDgwIvZgexuxRxr5DVI7mXyp132D4mrG2N8Hm5k5MzcAf0nR0PcB4GbKxviZeSPwTxRtpu4FDqa4/dWMj1EkKd8GHqKoTPsz81cUScG7gF9RNLg+NjPvbzDu31NU9i9j+1ULmik3M79DkeBcD1zDHyZonwCOj+Lpw0/WCeM0iqsft1C0IbqEIqmYqUso1tMtFLfa/qEm1lHgKxQPTKydppx7KLblJorbk2/PzP9XjjuDYhv/W0Q8BHyXoq3RtDJzC0Wy9QqKq1f/k6Ld1FjZp1LcFrsHuIgisZnMu4GNwE+BXwMfoc75PDOvpniI4TzgQeD/sv1q0Wsp2kZtAv4X8IHM/G4jv6V0CRP2n2bKzcwHKdoink9xZew3FA31x3y5/PdXEXFtneVfAPwr8H3gVoqrgac1Eb/UNpE53ZVdSep+EfF+4FmZeeIU0xxB0ch7/8oCk9Rz7ARQ0pxX3op8K8XTnJLUVt52lDSnRdHJ553ANzPz++2OR5K87ShJklQhr3xJkiRVyORLkiSpQh3d4H7ffffNRYsWtTsMSZKkaV1zzTX3Z+aC6abr6ORr0aJFbNiwod1hSJIkTSsiJnst2DjedpQkSaqQyZckSVKFTL4kSZIqZPIlSZJUIZMvSZKkCpl8SZIkVcjkS5IkqUImX5IkSRUy+ZIkSapQpclXRPxNRNwQET+PiNURsWuVy5ckSWq3yl4vFBEDwH8HDsrM0Yi4DDgBuKiqGKRWWjc4zMr1Q2waGWXh/H6WL1vMcUsGZmX6bi3bWHr7d0qqr+p3O+4M9EfEVmA3YFPFy5ca1mwFtmLtRka3bgNgeGSUFWs3AtSdp5npu7VsY+nt31k7j8maNF5ltx0zcxj4KHAHcDfwYGZ+u6rlS80Yq2SGR0ZJtlcy6waH606/cv3QYxXSmNGt21i5fmiHp+/Wso2lt38nNH8cSb2isuQrIp4AvAo4EFgI7B4RJ9aZ7uSI2BARGzZv3lxVeNI4zVYym0ZGWza8W8s2lt7+ndD8cST1iiob3L8MuDUzN2fmVmAt8MKJE2XmqsxcmplLFyxYUGF46gXrBoc5/NzLOfDMr3P4uZdP+j/wZiuZhfP7Wza8W8s2lt7+ndD8cdTo8Sl1uyqTrzuAF0TEbhERwEuBmypcvnpcM7dAmq1kli9bTP+8vnHD+uf1sXzZ4h2evlvLNpbe/p3Q3HHkLUr1kirbfF0FrAGuBTaWy15V1fKlZm6BNFvJHLdkgHNefTC79BWH1MD8fs559cGTNixuZvpuLdtYevt3QnPHkbco1UsiM9sdw6SWLl2aGzZsaHcYmiMOPPPr1NvbA7j13GP+YPi6wWFOX3M9W7Y9ykCDT2m95jM/AeDSUw5rKKZmpu/Wso2l+rI7KZZGj6Nmj0+pE0XENZm5dLrpqu5qQmqbhfP7Ga7T1mSyWyPHLRlg9dV3AI1XSpLGa/Q4avb4lLqZrxdSz2j2VqKk6nh8qpd45UtdrZkOHMeGN3srUVLrzeT4tANXdSuTL3WtmfS27a1EqXM1c3zO5PiXOoW3HdW1fDpK6l0e/+pmJl/qWs124Chp7vD4Vzcz+VLXarYjVElzh8e/upnJl7qWT0dJvcvjX93MBvfqWj69KPUuj391M5MvdTWfXpR6l8e/upXJlzqK/fZIagXPLeokJl/qGPbbI6kVPLeo09jgXh3DfnsktYLnFnUaky91DPvtkdQKnlvUaUy+1DHst0dSK3huUacx+VLHsN8eSa3guUWdxgb36hj22yOpFTy3qNOYfKmj2G+PpFbw3KJO4m1HSZKkCpl8SZIkVcjkS5IkqUK2+VLL+VoPSd3G85ZayeRLLeVrPSR1G89bajVvO6qlfK2HpG7jeUutZvKllvK1HpK6jecttZrJl1rK13pI6jaet9RqJl9qKV/rIanbeN5Sq9ngXi3laz0kdRvPW2q1SpOviJgPnA/8EZDAWzLzJ1XGoOr5Wg9J3cbzllqp6itfnwC+lZnHR8QuwG4VL1+SJKmtKku+ImIv4CXASQCZuQXYUtXyJUmSOkGVDe4PBDYDF0bEYEScHxG7V7h8SZKktqsy+doZ+GPg05m5BPgNcObEiSLi5IjYEBEbNm/eXGF4kiRJrVdlm6+7gLsy86ry+xrqJF+ZuQpYBbB06dKsLjw1yneeSdJ4nhfVjMqSr8y8JyLujIjFmTkEvBS4sarla3b4zjNJGs/zoppVdSerpwFfjIjrgecCZ1e8fO0g33kmSeN5XlSzKu1qIjOvA5ZWuUzNLt95JknjeV5Us3y9kJriO88kaTzPi2qWyZea4jvPJGk8z4tqlu92VFN855kkjed5Uc0y+VLTfOeZJI3neVHN8LajJElShUy+JEmSKmTyJUmSVCGTL0mSpAqZfEmSJFXI5EuSJKlCJl+SJEkVsp8vAbBucJiV64fYNDLKQjsIlKSW8Xwrky+xbnCYFWs3Mrp1GwDDI6OsWLsRwBOCJM0iz7cCbzsKWLl+6LETwZjRrdtYuX6oTRFJ0tzk+VZg8iVg08hoU8MlSTPj+VZg8iVg4fz+poZLkmbG862ggeQrIvZu4G9+FcGqNZYvW0z/vL5xw/rn9bF82eI2RSRJc5PnW0FjDe43lX8xxTR9wFNmJSJVbqyR5+lrrmfLtkcZ8OkbSWoJz7eCxpKvmzJzyVQTRMTgLMWjNjluyQCrr74DgEtPOazN0UjS3OX5Vo20+ToMICL+YeKIiOirnUaSJElTmzb5yszflR8HIuJ1Y8Mj4onAdydMI0mSpCk008nqKcD6iLgZSOBC4IyWRCVJkjRHTZt8RcQXgGuBQeAdwCXAI8BxmXlza8OTJEmaWxpp83URxZOObwYuBhYBDwAnRsTxLYtMkiRpDpr2yldmXg5cPvY9InYGngMcCjwfWNOy6CRJkuaYpl+snZmPABvLv4tnPSJJkqQ5rJEe7q+djWkkSZLU2JWv50TE9VOMD2CvWYpHkiRpTmsk+Xp2A9Nsa3SBZcesG4DhzDy20fkkSZLmgkYa3N8OEBHfAd6dmT/bwWW+E7gJ2HMHy9EU1g0Os3L9EJtGRlnou8MkqWt5Pp97GulqYswZwMcj4sKI2G8mC4uI/YFjgPNnMr8as25wmBVrNzI8MkoCwyOjrFi7kXWDw+0OTZLUBM/nc1PDyVdmXpuZfwp8DfhWRHwgIvqbXN7HgdOBR5ucT01YuX6I0a3j7wSPbt3GyvVDbYpIkjQTns/npmaufBERAQwBnwZOA/49It7Q4LzHAvdl5jXTTHdyRGyIiA2bN29uJjyVNo2MNjVcktSZPJ/PTQ0nXxHxI2AYOA8YAE4CjgCeFxGrGijicOCVEXEb8CXgyIj4g37CMnNVZi7NzKULFixoNDzVWDi//gXJyYZLkjqT5/O5qZkrXycDA5n58sz8u8z8WmbenJmnAS+ebubMXJGZ+2fmIuAE4PLMPHFmYWsqy5ctpn9e37hh/fP6WL5scZsikiTNhOfzuanhHu4z84YpRh8zC7Folow9BXP6muvZsu1RBnw6RpK6kufzuanp1wvVk5m3NDn9lcCVs7Fs1XfckgFWX30HAJeeclibo5EkzZTn87mnqQb3kiRJ2jEmX5IkSRUy+ZIkSaqQyZckSVKFTL4kSZIqZPIlSZJUIZMvSZKkCpl8SZIkVcjkS5IkqUKz0sO9Wm/d4DAr1w+xaWSUhb5eQpI0CeuLzmfy1QXWDQ6zYu1GRrduA2B4ZJQVazcCeEBJkh5jfdEdvO3YBVauH3rsQBozunUbK9cPtSkiSVInsr7oDiZfXWDTyGhTwyVJvcn6ojuYfHWBhfP7mxouSepN1hfdweSrCyxftpj+eX3jhvXP62P5ssVtikiS1ImsL7qDDe67wFgjydPXXM+WbY8y4NMrkqQ6rC+6g8lXlzhuyQCrr74DgEtPOazN0UiSOpX1RefztqMkSVKFTL4kSZIqZPIlSZJUIZMvSZKkCpl8SZIkVcjkS5IkqUImX5IkSRWyn682Wjc4zMr1Q2waGWWhHeFJkipmPdQeJl9tsm5wmBVrNz729vnhkVFWrN0I4I4vSWo566H28bZjm6xcP/TYDj9mdOs2Vq4falNEkqReYj3UPiZfbbJpZLSp4ZIkzSbrofapLPmKiAMi4oqIuDEiboiId1a17E60cH5/U8MlSZpN1kPtU+WVr0eAd2XmQcALgHdExEEVLr+jLF+2mP55feOG9c/rY/myxW2KSJLUS6yH2qeyBveZeTdwd/n54Yi4CRgAbqwqhk4y1pjx9DXXs2Xbowz4lIkkqULWQ+3TlqcdI2IRsAS4qh3L7xTHLRlg9dV3AHDpKYe1ORpJUq+xHmqPyhvcR8QewFeAv87Mh+qMPzkiNkTEhs2bN1cdniRJUktVmnxFxDyKxOuLmbm23jSZuSozl2bm0gULFlQZniRJUstV+bRjAJ8DbsrMj1W1XEmSpE5S5ZWvw4E3AEdGxHXl39EVLl+SJKntqnza8YdAVLU8SZKkTuS7HWeRLyiVJM1l1nOzw+RrlviCUknSXGY9N3t8t+Ms8QWlkqS5zHpu9ph8zRJfUCpJmsus52aPydcs8QWlkqS5zHpu9ph8zRJfUCpJmsus52aPDe5niS8olSTNZdZzs8fkaxb5glJJ0lxmPTc7vO0oSZJUIZMvSZKkCnnbcQr25CtJ0sxYh07O5GsS9uQrSdLMWIdOzduOk7AnX0mSZsY6dGomX5OwJ19JkmbGOnRqJl+TsCdfSZJmxjp0aiZfk7AnX0mSZsY6dGo2uJ+EPflKkjQz1qFTM/magj35SpI0M9ahk+u55Mt+RyRJ6jy9VD/3VPJlvyOSJHWeXqufe6rBvf2OSJLUeXqtfu6p5Mt+RyRJ6jy9Vj/3VPJlvyOSJHWeXqufeyr5st8RSZI6T6/Vz13f4L6ZpyPsd0SSpM4zk/q5m5+O7OrkayZPR9jviCRJnaeZ+rnbn47s6tuOvfZ0hCRJ6v76v6uTr157OkKSJHV//V9p8hURR0XEUETcHBFnTjf9xuEHOfzcy1k3OFx3fK89HSFJkpqv/9cNDnP4uZdz4JlfnzKvqEplyVdE9AH/DLwCOAh4bUQcNN18Y/dx662oXns6QpIkNVf/j7UPGx4ZJZk6r6hKlVe+ngfcnJm3ZOYW4EvAqxqZcbL7uMctGeCcVx/MLn3FzxiY3885rz64KxrbSZKkmWmm/u/E9mGRmdUsKOJ44KjMfFv5/Q3A8zPz1Mnmefoee+XZh77ose8veNo+dae78e6HADhovz0biqWZ6VtZtrFUX7axVF+2sVRftrFUX7axVF92o9P/2y2/euzzLXsN8JlDius+Adx67jENLadREXFNZi6dbrqO62oiIk4GTgbYqX9P3njLLwDIbY9s2frj2za2MzbtsH2B+9sdhGaN23NucXvOPW5TYN6CRQdH3867FN9+AYNXAEVeER85drbziqc2MlGVydcwcEDN9/3LYeNk5ipgFUBEbPj9bx+cNoNUd4iIDY38j0Ddwe05t7g95x63aeeqss3XT4FnRsSBEbELcALw1QqXL0mS1HaVXfnKzEci4lRgPdAHXJCZN1S1fEmSpE5QaZuvzPwG8I0mZlnVqljUFm7PucXtObe4Pecet2mHquxpR0mSJHX564UkSZK6TUcmX82+hkidJyIuiIj7IuLnNcP2jojvRMS/l/8+oZ0xqnERcUBEXBERN0bEDRHxznK427QLRcSuEXF1RPys3J4fKocfGBFXlefeS8uHo9QlIqIvIgYj4mvld7dnh+q45GumryFSx7kIOGrCsDOB72XmM4Hvld/VHR4B3pWZBwEvAN5RHpdu0+70e+DIzDwUeC5wVES8APgIcF5mPgN4AHhrG2NU894J3FTz3e3ZoTou+WIHXkOkzpGZ3wd+PWHwq4DPl58/DxxXaVCascy8OzOvLT8/THGCH8Bt2pWy8B/l13nlXwJHAmvK4W7PLhIR+wPHAOeX3wO3Z8fqxORrALiz5vtd5TB1vydl5t3l53uAJ7UzGM1MRCwClgBX4TbtWuUtquuA+4DvAL8ERjLzkXISz73d5ePA6cCj5fd9cHt2rE5MvtQDsnjM1kdtu0xE7AF8BfjrzHyodpzbtLtk5rbMfC7F20aeBzy7zSFphiLiWOC+zLym3bGoMR33bkcafA2RutK9EbFfZt4dEftR/I9bXSIi5lEkXl/MzLXlYLdpl8vMkYi4AjgMmB8RO5dXSzz3do/DgVdGxNHArsCewCdwe3asTrzy5WuI5q6vAm8qP78J+N9tjEVNKNuPfA64KTM/VjPKbdqFImJBRMwvP/cDL6dox3cFcHw5mduzS2TmiszcPzMXUdSZl2fm63F7dqyO7GS1zN4/zvbXEJ3V5pDUpIhYDRwB7AvcC3wAWAdcBjwFuB34i8yc2ChfHSgiXgT8ANjI9jYl76Fo9+U27TIRcQhFA+w+iv+EX5aZfx8RT6N4yGlvYBA4MTN/375I1ayIOAJ4d2Ye6/bsXB2ZfEmSJM1VnXjbUZIkac4y+ZIkSaqQyZckSVKFTL4kSZIqZPIlSZJUIZMvSZKkCpl8SZIkVcjkS1LPiIiXRsS/tjsOSb3N5EtSLzmUoqdvSWobky9JveRQYDAiHhcRF0XE2eV7KyWpMju3OwBJqtAhwH3AeuD8zLy4zfFI6kG+21FST4iIecD9FC8APyUzf9LmkCT1KG87SuoVzwF+CjwCbGtzLJJ6mMmXpF5xKPBj4ATgwoh4UpvjkdSjTL4k9YpDgZ9n5i+AM4DLyluRklQp23xJkiRVyCtfkiRJFTL5kiRJqpDJlyRJUoVMviRJkipk8iVJklQhky9JkqQKmXxJkiRVyORLkiSpQv8fHh2sB3PtGfgAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x112dcc6d8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import scipy.signal as sig\n",
    "\n",
    "L = 32  # length of signal x[k]\n",
    "N = 16  # length of signal h[k]\n",
    "M = 16  # periodicity of periodic convolution\n",
    "\n",
    "\n",
    "def periodic_summation(x, N):\n",
    "    \"Zero-padding to length N or periodic summation with period N.\"\n",
    "    M = len(x)\n",
    "    rows = int(np.ceil(M/N))\n",
    "    \n",
    "    if (M < int(N*rows)):\n",
    "        x = np.pad(x, (0, int(N*rows-M)), 'constant')\n",
    "    \n",
    "    x = np.reshape(x, (rows, N))\n",
    "    \n",
    "    return np.sum(x, axis=0)\n",
    "\n",
    "\n",
    "def periodic_convolve(x, y, P):\n",
    "    \"Periodic convolution of two signals x and y with period P.\"\n",
    "    x = periodic_summation(x, P)\n",
    "    h = periodic_summation(y, P)\n",
    "\n",
    "    return np.array([np.dot(np.roll(x[::-1], k+1), h) for k in range(P)], float)\n",
    "\n",
    "\n",
    "# generate signals\n",
    "x = np.ones(L)\n",
    "h = sig.triang(N)\n",
    "\n",
    "# linear convolution\n",
    "y1 = np.convolve(x, h, 'full')\n",
    "# periodic convolution\n",
    "y2 = periodic_convolve(x, h, M)\n",
    "# linear convolution via periodic convolution\n",
    "xp = np.append(x, np.zeros(N-1))\n",
    "hp = np.append(h, np.zeros(L-1))\n",
    "y3 = periodic_convolve(xp, hp, L+N-1)\n",
    "\n",
    "# plot results\n",
    "def plot_signal(x):\n",
    "    plt.figure(figsize = (10, 3))\n",
    "    plt.stem(x)\n",
    "    plt.xlabel(r'$k$')\n",
    "    plt.ylabel(r'$y[k]$')\n",
    "    plt.axis([0, N+L, 0, 1.1*x.max()])\n",
    "\n",
    "\n",
    "plot_signal(x)\n",
    "plt.title('Signal $x[k]$')\n",
    "plot_signal(y1)\n",
    "plt.title('Linear convolution')\n",
    "plot_signal(y2)\n",
    "plt.title('Periodic convolution with period M = %d' %M)\n",
    "plot_signal(y3)\n",
    "plt.title('Linear convolution by periodic convolution');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise**\n",
    "\n",
    "* Change the lengths `L`, `N` and `M` and check how the results for the different convolutions change"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The Fast Convolution\n",
    "\n",
    "Using the above derived equality of the linear and periodic convolution one can express the linear convolution $y[k] = x_L[k] * h_N[k]$ by the DFT as\n",
    "\n",
    "\\begin{equation}\n",
    "y[k] = \\text{IDFT}_M \\{ \\; \\text{DFT}_M\\{ x_M[k] \\} \\cdot \\text{DFT}_M\\{ h_M[k] \\} \\; \\}\n",
    "\\end{equation}\n",
    "\n",
    "This operation requires three DFTs of length $M$ and $M$ complex multiplications. On first sight this does not seem to be an improvement, since one DFT/IDFT requires $M^2$ complex multiplications and $M \\cdot (M-1)$ complex additions. The overall numerical complexity is hence in the order of $\\mathcal{O}(M^2)$. The DFT can be realized efficiently by the [fast Fourier transformation](https://en.wikipedia.org/wiki/Fast_Fourier_transform) (FFT), which lowers the computational complexity to $\\mathcal{O}(M \\log_2 M)$. The resulting algorithm is known as *fast convolution* due to its computational efficiency. \n",
    "\n",
    "The fast convolution algorithm is composed of the following steps\n",
    "\n",
    "1. Zero-padding of the two input signals $x_L[k]$ and $h_N[k]$ to at least a total length of $M \\geq N+L-1$\n",
    "\n",
    "2. Computation of the DFTs $X[\\mu]$ and $H[\\mu]$ using a FFT of length $M$\n",
    "\n",
    "3. Multiplication of the spectra $Y[\\mu] = X[\\mu] \\cdot H[\\mu]$\n",
    "\n",
    "4. Inverse DFT of $Y[\\mu]$ using an inverse FFT of length $M$\n",
    "\n",
    "The overall complexity depends on the particular implementation of the FFT. Many FFTs are most efficient for lengths which are a power of two. It therefore can make sense, in terms of computational complexity, to choose $M$ as a power of two instead of the shortest possible length $N+L-1$. For real valued signals $x[k] \\in \\mathbb{R}$ and $h[k] \\in \\mathbb{R}$ the computational complexity can be reduced significantly by using a real valued FFT."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Example - Fast convolution\n",
    "\n",
    "The implementation of the fast convolution algorithm is straightforward. Most implementations of the FFT include zero-padding to a given length $M$, e.g in `numpy` by `numpy.fft.fft(x, M)`. In the following example an implementation of the fast convolution is shown. For illustration the convolution of a rectangular signal $x[k] = \\text{rect}_L[k]$ of length $L$ with a triangular signal $h[k] = \\Lambda_N[k]$ of length $N$ is considered."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAsgAAAGoCAYAAABbtxOxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJzt3X2YXXV56P3v3UmgA1giBikZCAG1ESyGeFKRYj0UbKPFag4XpUKhah8lPo9arBogtop61HAZ68s5tVqqgMi7mEaqHke5glZbBUMGEyUnikiACUhQR0BHSCb388deAzvjvOw9s/dee+/5fq5rrtl7vfzWvWfNWuueNff6/SIzkSRJklTxW2UHIEmSJLUTE2RJkiSpigmyJEmSVMUEWZIkSapigixJkiRVMUGWJEmSqpggS5p1IuKkiLivge2dGBE/jIhHI2LFOPPvjogXF6/fHhGfbNS2u0lEZEQ8c5rr/lFEbGt0TJJmJxNkSaUqksfhIrl8ICIuj4gDSojhxTNo4j3AP2XmAZm5frIFM/P9mfnaGWxL/GYynZnfyMzFZcYkqXuYIEtqB3+emQcAxwFLgdUlx1OvI4Dvlx3EZCKip+wYJKlTmCBLahuZ+QDQTyVRBiAi9o2ID0bEPRHxk4j4RET0FvPmR8QXImIoIn4WEd+IiN8q5u11h7G4M/3esduMiM8AC4F/L+5inz9ebBHxuoi4s9jOjRGxoJj+I+CoqvX3newzRsS7IuLK4vWiIs5XFZ/voYj4+6plfysiLoyIH0XETyPi+og4qGr+Z4u77r+IiP+IiOeM+bwfj4gvRcQvgT8eJ5aDIuKyiNgRET+PiPVV88b9vFU/29cXZSVDEfGxqNi3eP/7VcseXPyH4OlTtTsmtq9FxGur3r86Ir5ZvP6PYvJ3i5/5X44tm4mIo4s2hiLi+xHx8jE/m49FxBcj4pGIuCUinjHZfpM0u5ggS2obEXEY8FLgzqrJFwO/RyVpfibQB7yzmPdW4D7gYOAQ4O1A1rPNzDwHuIfiLnZmfmCcuE4G1gBnAIcC24Fri/WfMWb9x+rZfuGFwGLgFOCdEXF0Mf1NwArgvwMLgJ8DH6ta7/8AzwKeDmwCrhrT7lnA+4CnAN8cZ7ufAfYDnlO08eGpPm+VlwF/ADy3WG558dnXAWdWLXcG8PXMfLDGdqeUmS8qXi4pfubXVc+PiLnAvwNfKT7Xm4CrIqK6BOOVwLuBp1L5fXtfvXFI6l4myJLawfqIeAS4F3gQuAggIgI4F/i7zPxZZj4CvJ9KcgOwi0qidURm7irqUOtKkGv0V8ClmbmpSAJXAydExKIGtf/uzBzOzO8C3wWWFNNfD/x9Zt5XbPddwOkRMQcgMy/NzEeq5i2JiAOr2v18Zv5nZu7JzF9XbzAiDqXyx8jrM/Pnxc/v63V83oszcygz7wFu5sm7/lfz5P6BSpJ+dR3tNsILgAOKGB/PzA3AF9g7cf+3zLw1M3dT+cPiuHHakTRLmSBLagcrMvMpwEnAs4H5xfSDqdzhvK34V/kQ8OViOsBaKnf/vhIRd0XEhU2KbwGVu50AZOajwE+p3M1uhAeqXv+KSnIHldrmf6v67FuBEeCQiOiJiIuL8ouHgbuLdeZXtXXvJNs8HPhZZv58nHm1fN6JYr4Z2C8iji8S3+OAf6uj3UZYANybmXuqpm2ntvglyQRZUvso7mBeDnywmPQQMAw8JzPnFV8HFg/0Udw9fWtmHgW8HHhLRJxSrPsrKsn1qN+dbNNThLaDSrIKQETsDzwNGKztk03bvcBLqz77vMz87cwcpHJn9hXAi4EDgUWj4VWtP9nnuhc4KCLmjTNv2p83M0eA66ncrT0T+EJx57/edn9J7ftvvPgPH61HLyysJX5JAhNkSe3nI8CfRMSS4g7gvwIfrnrIqy8ilhevXxYRzyxKMX5B5e7q6F3D24GzijutL6FSxzuRn1B50G4i1wCviYjjiofw3g/ckpl3T/9j1uQTwPsi4gh44oG3VxTzngI8RuUO7H5FTDXLzPup1DD/c0Q8NSLmRsRobe9MP+/VwF9SKam4ump6Pe3eDpwWEftF5WHL/2fM/Mn22S1U/kA6v/hcJwF/zjTqnSXNTibIktpKZu4EruDJB/EuoFJG8e2ilOAmKg+0QeUBtZuAR4FvAf+cmTcX886jkhQNUUnUJuufeA3wD0Upw9vGiekm4B3A54D7gWewd51ts3wUuJFKCckjwLeB44t5V1ApGxgE7ijm1escKnXc/5dK7febYeafNzNvoXIHeAGVJHx0ej3tfhh4nEoi/Gl+8wHEdwGfLvbZGWO2/ziVff9SKv+F+GfgrzPz/9b6GSTNbtGc51kkSZKkzuQdZEmSJKmKCbIkSZJUxQRZkiRJqmKCLEmSJFWZU3YAk5k/f34uWrSo7DAkSZLUBW677baHMvPgqZZr6wR50aJFbNy4sewwJEmS1AUiYvvUS1liIUmSJO3FBFmSJEmq0tISi4j4O+C1QAJbgNdk5q9bGYMkdYv1A4Os7d/GjqFhFszrZdXyxaxY2ld2WJLU8VqWIEdEH/C3wDGZORwR11MZYvTyVsUgSWVqZEK7fmCQ1eu2MLxrBIDBoWFWr9sCMKM2TbglqfUlFnOA3oiYA+wH7Gjx9iWpFKMJ7eDQMMmTCe36gcFptbe2f9sTyfGo4V0jrO3f1hbxSVIna1mCnJmDwAeBe4D7gV9k5lfGLhcR50bExojYuHPnzlaFJ0lN1eiEdsfQcF3Tp9Lo+CSpk7UsQY6IpwKvAI4EFgD7R8TZY5fLzEsyc1lmLjv44Cm7qZOkjtDohHbBvN66pk+l0fFJUidrZYnFi4EfZ+bOzNwFrAP+sIXbl6S6rB8Y5MSLN3DkhV/kxIs3zKjcoNEJ7arli+md27PXtN65Paxavnha7TU6vlGN/BlKUqu0MkG+B3hBROwXEQGcAmxt4fYlqWaNrsltdEK7Ymkfa047ln16Kqfxvnm9rDnt2Gk/VNfo+MC6Zkmdq5U1yLcANwCbqHTx9lvAJa3aviTVo9E1uY1OaEfbXLpwHscfeRD/eeHJM26r0fFZ1yypU7W0H+TMvAi4qJXblKTpaEZN7oqlfVxz6z0AXLfyhGm30yyNjs+6ZkmdypH0JGkczarJnU38GUrqVCbIkrpGIx8Ia0ZN7mzTrLpmH/qT1GwtLbGQpGZp9Mhyo+ucf8NmHh/ZQ58jy9Wt0T/DZoweKEnjMUGW1BUmeyBsuslTu9cMd4JG/gybsY8laTyWWEjqCj4Q1v3cx5JaxQRZUlfwgbDu5z6W1ComyJK6gg/VdT/3saRWsQZZUlfwobru5z6W1ComyJJKs35gkLX929gxNMyCBiQ7PlTX/Rq9jxv9OyipO5ggSyqFXXapbP4OSpqINciSSjFZl11SK/g7KGkiJsiSSmGXXSqbv4OSJmKCLKkUdtmlsvk7KGkiJsiSSmGXXSqbv4OSJuJDepJKYZddKpu/g5ImYoIsqTR2y6ay+TsoaTyWWEiSJElVTJAlSZKkKpZYSKqZo45Jk/MYkbrDlAlyRBxUQzt7MnOoAfFIalOOOiZNzmNE6h613EHeUXzFJMv0AAsbEpGktjTZqGNe/CWPEamb1JIgb83MpZMtEBEDDYpHUpty1DFpch4jUveo5SG9Wvq9sW8cqcs56pg0OY8RqXtMmSBn5q8BIuK9Y+dFRE/1MpK6l6OOSZPzGJG6Rz3dvPVFxJmjbyLi6cBNjQ9JUjtasbSPNacdyz49ldNG37xe1px2rLWVUsFjROoe9XTzthLoj4gfAQlcBlxQz8YiYh7wSeD3izb+JjO/VU8bksrjqGPS5DxGpO5QSzdvVwCbgAHgDcDVwG5gRWbeWef2Pgp8OTNPj4h9gP3qXF+SJElqqlpKLC6n0sXba4ArgUXAz4GzI+L0WjcUEQcCLwI+BZCZj9t3siRJktrNlHeQM3MDsGH0fUTMAY4GlgDHAzfUuK0jgZ3AZRGxBLgNOC8zf1lv0JIkSVKz1POQHgCZuTszt2TmlZm5qo5V5wDPAz5e9Kv8S+DCsQtFxLkRsTEiNu7cubPe8CRJkqQZqaUGeVNmPm+mywD3Afdl5i3F+xsYJ0HOzEuASwCWLVuWU8UnaWLrBwZZ27+NHUPDLJjXy6rli32iXuogHsNSOWrpxeLoiNg8yfwADpyqkcx8ICLujYjFmbkNOAW4o8Y4JdVp/cAgq9dteWLo28GhYVav2wLgBVbqAB7DUnlqSZCfXcMyI1MvAsCbgKuKHizuovLgn6QmWNu/7YkL66jhXSOs7d/mxVXqAB7DUnlqeUhve6M2lpm3A8sa1Z6kie0YGq5ruqT24jEslaeegUIAiIjrqZRV3E2lf+RNRcmEpDayYF4vg+NcSBfM6y0hGkn18hiWyjOdXizOyMy/AD4BvBD4TsOjkjRjq5Yvpnduz17Teuf2sGr54pIiklQPj2GpPNO5g/xi4FQqo+B9A3h7o4OSNHOjNYrn37CZx0f20OcT8FJH8RiWylN3ggxcCnwZ+BqwMTN/0dCIJDXMiqV9XHPrPQBct/KEkqORVC+PYakcdSfImbkwIg4D/huV4aaflZlnNj40SZIkqfVqTpAj4mTgr4Ah4HvAZuDLmflYk2KTJEmSWq6eO8iXAm8G5gLPBVYAzwGe2YS4JEmSpFLUkyBvz8z1xevPNiMYSZIkqWxTdvMWEVdExJuBb0fEW1oQkyRJklSaWvpBvpzKwCCHAOdExPaIuDEi/mdE/EVTo5MkSZJarJahpjcAG0bfR8Qc4GhgCfAHWG4hNcz6gUHW9m9jx9AwC+zzVFKDeY6RajOdbt52A1uKL0kNsn5gkNXrtjC8awSAwaFhVq+rHGZewCTNlOcYqXZ1DzUtqTnW9m974sI1anjXCGv7t5UUkaRu4jlGql1NCXJEnFV8f2Vzw5Fmrx1Dw3VNl6R6eI6RalfrHeS+iDgDOKyZwUiz2YJ5vXVNl6R6eI6RaldLN28XAQcBVwEHRcQ7mx6VNAutWr6Y3rk9e03rndvDquWLS4pIUjfxHCPVbsoEOTPfDfwUOAf4aWa+p+lRSbPQiqV9rDntWPbpqRyWffN6WXPasT48I6khPMdItau1F4v7M/PaiDizqdFIs9yKpX1cc+s9AFy38oSSo5HUbTzHSLWpqQY5M68qvl/T3HAkSZKkctnNmyRJklRlyhKLiDiohnb2ZOZQA+KRJEmSSlVLDfKO4ismWaYHWNiQiCRJkqQS1ZIgb83MpZMtEBEDDYpHkiRJKlUtNcgnAETEe8fOiIie6mUkSZKkTldLP8i/Ll72jQ45DRARTwduGrOMJEmS1NFq7QcZYCXQHxF3AglcBlzQlKgkSZKkktTSi8UVwCZgAHgDcDWwG1iRmXfWu8GiLGMjMJiZL6t3faldrB8YZG3/NnYMDbNgXi+rli92RCpJs4rnQXWrWmqQL6fSg8VrgCuBRcDPgbMj4vRpbPM8YOs01pPaxvqBQVav28Lg0DAJDA4Ns3rdFtYPDJYdmiS1hOdBdbNaapA3ZOaHM/PVmfk8YD7wFuBO4Ph6NhYRhwGnAp+cTrBSu1jbv43hXSN7TRveNcLa/m0lRSRJreV5UN2snhpkADJzN7Cl+LqyztU/ApwPPGWiBSLiXOBcgIUL7VpZ7WnH0HBd0yWp23geVDeb8g5yRGxq0DIvAx7MzNsmWy4zL8nMZZm57OCDD56qWakUC+b11jVdkrqN50F1s1pqkI+OiM2TfG2hUnYxlROBl0fE3cC1wMkRUe8daKktrFq+mN65PXtN653bw6rli0uKSJJay/OgulktJRbPrmGZkakWyMzVwGqAiDgJeFtmnl1D21LbGX1K+/wbNvP4yB76fHpb0izjeVDdbMoEOTO3tyIQqdOsWNrHNbfeA8B1Kx1MUtLs43lQ3aqWEotJRcQf1btOZn7NPpAlSZLUjmacIAN/0YA2JEmSpLZQdzdvEXEj8GMqo+vdNp02JEmSpHZV8x3kiPhoRERmvhz4EPAw8JfAEc0KTpIkSWq1ekosHgFujIj9igf3fgWcnJmnNic0SZIkqfVqLo/IzH+IiLOAr0fE48CjwIVNi0ySJEkqQc0JckScArwO+CVwKPA3memA65IkSeoq9ZRY/D3wjsw8CTgduC4iTm5KVJIkSVJJ6imxOLnq9ZaIeCnwOeAPmxGY1GjrBwZZ27+NHUPDLHDEJ0lqO56n1S6m3UVbZt5flF1IbW/9wCCr121heFdlVPTBoWFWr9sC4MlXktqA52m1kxkNFJKZw40KRGqmtf3bnjjpjhreNcLafsvoJakdeJ5WO2nESHpS29sxNP7fchNNlyS1ludptZO6EuTRh/J8OE+dZsG83rqmS5Jay/O02km9d5A/OOa71BFWLV9M79yevab1zu1h1fLFJUUkSarmeVrtZLoP6UVDo5CabPQBj/Nv2MzjI3vo8+loSWornqfVTqbdi4XUaVYs7eOaW+8B4LqVJ5QcjSRpLM/Tahc+pCdJkiRVMUGWJEmSqtSbID9afH+k0YFIkiRJ7aCuBDkzX1T9XZIkSeo2llhIkiRJVUyQJUmSpCpTdvMWEQfV0M6ezBxqQDzSE9YPDLK2fxs7hoZZYH+YkqQ6eR3RdNXSD/KO4muywUF6gIUNiUiiclJbvW4Lw7tGABgcGmb1ui0AntwkSVPyOqKZqKXEYmtmHpWZR070Bfy02YFqdlnbv+2Jk9qo4V0jrO3fVlJEkqRO4nVEM1FLgnwCQES8d+yMiOipXkZqlB1Dw3VNlySpmtcRzcSUCXJm/rp42RcRZ41Oj4inAzeNWWZCEXF4RNwcEXdExPcj4rzpBq3ut2Beb13TJUmq5nVEM1FPLxYrgddFxPMj4g+ADcAH61h/N/DWzDwGeAHwhog4po71NYusWr6Y3rk9e03rndvDquWLS4pIktRJvI5oJmrpxeIKYBMwALwBuJpKsrsiM++sdUOZeT9wf/H6kYjYCvQBd0wjbnW50Qcozr9hM4+P7KHPp48lSXXwOqKZqKUXi8uBJcBrgOcCi4DvAGdHxPcy84Z6NxoRi4ClwC3jzDsXOBdg4UI7xpjNVizt45pb7wHgupWWuUuS6uN1RNM1ZYKcmRuolFMAEBFzgKOpJM3HA3UlyBFxAPA54M2Z+fA427sEuARg2bJlWU/bkiRJ0kzVcgd5L5m5G9hSfF1Zz7oRMZdKcnxVZq6rd9uSJElSs035kF5EbGrQMgF8ikq/yh+qLTxJkiSptWq5g3x0RGyeZH4AB9bQzonAOcCWiLi9mPb2zPxSDetKkiRJLVFLgvzsGpYZmWqBzPwmkw9XrQ7nmPeSpG7ntW52qOUhve0AEfFV4G2Z+d2mR6WO45j3kqRu57Vu9qhnoJALgI9ExGURcWizAlJncsx7SVK381o3e9ScIGfmpsz8Y+ALwJcj4qKIcLxGAY55L0nqfl7rZo967iCP9kSxDfg48CbghxFxTjMCU2dxzHtJUrfzWjd71JwgR8R/AoPAh6kMEf1q4CTg+RFxSTOCU+dwzHtJUrfzWjd71DNQyLnAHZk5dnS7N0XE1gbGpA7kmPeSpG7ntW72qDlBzszvTzL71AbEog7nmPeSpG7ntW52qKsGeSKZeVcj2pEkSZLKVk+JhbqMnZ1LklQur8XtyQR5lrKzc0mSyuW1uH01pMRCncfOziVJKpfX4vZlgjxL2dm5JEnl8lrcvkyQZyk7O5ckqVxei9uXCfIsZWfnkiSVy2tx+/IhvVnKzs4lSSqX1+L2ZYI8i9nZuSRJ5fJa3J5MkDuIfSVKkqTJmCs0hglyh7CvREmSNBlzhcbxIb0OYV+JkiRpMuYKjWOC3CHsK1GSJE3GXKFxTJA7hH0lSpKkyZgrNI4JchOtHxjkxIs3cOSFX+TEizewfmBw2m3ZV6IkSZpMM3KFRuYyncSH9Jqk0YXy9pUoSZIm0+hcYTY/9GeC3CSTFcpP95fKvhIlSdJkGpkrNCOX6RQmyFUa2XeghfKSJKmTNSOX6ZR+mltagxwRL4mIbRFxZ0Rc2MptT2X03wiDQ8MkT/4bYbq1NhbKS5KkTtboXKbRuVYztSxBjoge4GPAS4FjgDMj4phWbX8qje470IfqJElSJ2t0LtNJ/TRHZrZmQxEnAO/KzOXF+9UAmblmonWWLVuWGzdubEl8R174RRJYufnzHPWLvf+SecFRT5tWmw89+hg/2vlLMpN95/Rw+EG9zD9g3xnFecf9DwNwzKG/M6N2bK8922tGm7Zne2W214w2bc/2ymyvGW22c3uNzGW+fddPn3h914F9/MtzXwFAAD+++NQZx1qLiLgtM5dNtVwra5D7gHur3t8HHD92oYg4FzgXYOHCha2JjMq/CwbHqanZd07POEvXZv4B+/LoY7sBWPS0/afdTrX99pl+PLbX/u01o03bs70y22tGm7Zne2W214w227m9RuYy+87p4bHdI78xvR3LT1t5B/l04CWZ+dri/TnA8Zn5xonWaeUd5LFdmUDl3whrTju2LYvHJUmSOkk75FrteAd5EDi86v1hxbS2MLpjOuHJSkmSpE7TSblWK+8gzwF+AJxCJTH+DnBWZn5/onVaeQdZkiRJ3a3t7iBn5u6IeCPQD/QAl06WHEuSJEllaOlAIZn5JeBLrdymJEmSVI+WDhQiSZIktbuW1SBPR0TsBLaXsOn5wEMlbFfjc3+0H/dJe3F/tB/3SXtxf7SfsvbJEZl58FQLtXWCXJaI2FhLAbdaw/3Rftwn7cX90X7cJ+3F/dF+2n2fWGIhSZIkVTFBliRJkqqYII/vkrID0F7cH+3HfdJe3B/tx33SXtwf7aet94k1yJIkSVIV7yBLkiRJVUyQJUmSpComyFUi4iURsS0i7oyIC8uORxARd0fEloi4PSI2lh3PbBQRl0bEgxHxvappB0XEVyPih8X3p5YZ42wywf54V0QMFsfJ7RHxZ2XGOJtExOERcXNE3BER34+I84rpHiMlmWSfeJyUICJ+OyJujYjvFvvj3cX0IyPiliLnui4i9ik71mrWIBciogf4AfAnwH3Ad4AzM/OOUgOb5SLibmBZZtrBe0ki4kXAo8AVmfn7xbQPAD/LzIuLPyafmpkXlBnnbDHB/ngX8GhmfrDM2GajiDgUODQzN0XEU4DbgBXAq/EYKcUk++QMPE5aLiIC2D8zH42IucA3gfOAtwDrMvPaiPgE8N3M/HiZsVbzDvKTng/cmZl3ZebjwLXAK0qOSSpdZv4H8LMxk18BfLp4/WkqFx+1wAT7QyXJzPszc1Px+hFgK9CHx0hpJtknKkFWPFq8nVt8JXAycEMxve2OERPkJ/UB91a9vw8PqHaQwFci4raIOLfsYPSEQzLz/uL1A8AhZQYjAN4YEZuLEgz/nV+CiFgELAVuwWOkLYzZJ+BxUoqI6ImI24EHga8CPwKGMnN3sUjb5VwmyGp3L8zM5wEvBd5Q/HtZbSQrdVrWapXr48AzgOOA+4F/LDec2SciDgA+B7w5Mx+unucxUo5x9onHSUkycyQzjwMOo/If+2eXHNKUTJCfNAgcXvX+sGKaSpSZg8X3B4F/o3JgqXw/Ker8Ruv9Hiw5nlktM39SXID2AP+Kx0lLFXWVnwOuysx1xWSPkRKNt088TsqXmUPAzcAJwLyImFPMarucywT5Sd8BnlU8VbkP8ErgxpJjmtUiYv/iAQsiYn/gT4HvTb6WWuRG4FXF61cBny8xlllvNBEr/A88TlqmeADpU8DWzPxQ1SyPkZJMtE88TsoREQdHxLzidS+VzhC2UkmUTy8Wa7tjxF4sqhRdvnwE6AEuzcz3lRzSrBYRR1G5awwwB7jafdJ6EXENcBIwH/gJcBGwHrgeWAhsB87ITB8ca4EJ9sdJVP5tnMDdwMqq+lc1UUS8EPgGsAXYU0x+O5WaV4+REkyyT87E46TlIuK5VB7C66FyY/b6zHxPcY2/FjgIGADOzszHyot0bybIkiRJUhVLLCRJkqQqJsiSJElSFRNkSZIkqYoJsiRJklTFBFmSJEmqYoIsSZIkVTFBliRJkqqYIEtSh4qIUyLiM2XHIUndxgRZkjrXEiojUEmSGsgEWZI61xJgICL2jYjLI+L9ERFlByVJnW5O2QFIkqbtucCDQD/wycy8suR4JKkrRGaWHYMkqU4RMRd4CNgOrMzMb5UckiR1DUssJKkzHQ18B9gNjJQciyR1FRNkSepMS4D/Al4JXBYRh5QcjyR1DRNkSepMS4DvZeYPgAuA64uyC0nSDFmDLEmSJFXxDrIkSZJUxQRZkiRJqmKCLEmSJFUxQZYkSZKqmCBLkiRJVUyQJUmSpComyJIkSVIVE2RJkiSpigmyJEmSVMUEWZIkSapigixJkiRVMUGWJEmSqpggS1IhIk6KiPsa2N6JEfHDiHg0IlaMM39xRNweEY9ExN82arvtLiJeHRHfnMH6n4iIdzQyJkmqZoIsqS1FxN0RMVwklw9ExOURcUAJMbx4Bk28B/inzDwgM9ePM/984ObMfEpm/q/pbiQivhYRr512lG1svGQ6M1+fmf+zrJgkdT8TZEnt7M8z8wDgOGApsLrkeOp1BPD9GcyXJJXABFlS28vMB4B+KokyABGxb0R8MCLuiYifFP927y3mzY+IL0TEUET8LCK+ERG/VczLiHhmVTuXR8R7x24zIj4DLAT+vbiLff54sUXE6yLizmI7N0bEgmL6j4Cjqtbfd8x6G4A/Bv6pmP97EXFqRAxExMMRcW9EvKtq+d+OiCsj4qfF5/pORBwSEe8D/qiqnX+aIM4XRsR/FeveGxGvLqYfGBFXRMTOiNgeEf9Q9bN6dUR8s/g5/zwifhwRLy3m/WVEbByzjb+LiBunanfMOouKfTKnatrXIuK1EXE08AnghOKzDY23zybaB8W8jIjXF6UuQxHxsYiI8X5GkjTKBFlS24uIw4CXAndWTb4Y+D0qSfMzgT7gncW8twL3AQcDhwBvB7KebWbmOcA9FHexM/MD48R1MrAGOAM4FNgOXFus/4wx6z82pv2TgW8Abyzm/wD4JfDXwDzgVOD/rapdfhVwIHA48DTg9cBwZv79mHbeOE6cRwD/B/jfxc/kOOD2YvamZYzHAAAaz0lEQVT/Lto9CvjvxfZfU7X68cA2YD7wAeBTRYL578DiiHhW1bJnAVfX2O6UMnNr8Tm/VXy2eeN8tgn3QZWXAX8APLdYbnk9cUiafUyQJbWz9RHxCHAv8CBwEUCRoJ0L/F1m/iwzHwHeD7yyWG8XlWTpiMzclZnfyMy6EuQa/RVwaWZuKhLg1VTudi6aTmOZ+bXM3JKZezJzM3ANleQSKp/pacAzM3MkM2/LzIdrbPos4KbMvKb4efw0M2+PiB4qP7PVmflIZt4N/CNwTtW62zPzXzNzBPg0lZ/rIZn5K+DzwJkARaL8bODGGtttlFr2wcWZOZSZ9wA3U/WfCEkajwmypHa2IjOfApxEJfmaX0w/GNgPuK34t/kQ8OViOsBaKnebvxIRd0XEhU2KbwGVO5YAZOajwE+p3M2uW0QcHxE3F2UJv6By93T0M3+GSpnJtRGxIyI+EBFza2z6cOBH40yfD8yt/gzF6+r4Hxh9USTFAKMPS15NkSBTScLXF8vU0m6j1LIPHqh6/SuejF+SxmWCLKntZebXgcuBDxaTHgKGgedk5rzi68DigT6Ku5ZvzcyjgJcDb4mIU4p1f0UluR71u5NteorQdlB50A6AiNifyl3ewdo+2W+4GrgRODwzD6RSfxsAxZ3fd2fmMcAfUikb+Osa47wXeMY40x+icmf6iKppC+uI/6vAwRFxHJVEebS8op52f1l8n2iftHofSJIJsqSO8RHgTyJiSWbuAf4V+HBEPB0gIvoiYnnx+mUR8cyiFOMXwAiwp2jnduCsiOiJiJfwZAnDeH5CpYZ2ItcAr4mI44qH8N4P3FKUFEzHU4CfZeavI+L5VO7KUnymP46IY4vyhYepJKCjn2mqOK8CXhwRZ0TEnIh4WkQcV5RNXA+8LyKeUtQqvwW4spZgM3MX8Fkqd+wPopIwU0+7mbmTSjJ7drFP/oa9k/mfAIdFxD4ThNHofSBJJsiSOkORSF3Bkw/iXUCljOLbEfEwcBOwuJj3rOL9o8C3gH/OzJuLeecBfw4MUalfHa9/4lFrgH8oyjjeNk5MNwHvAD4H3E8lsXvl2OXq8P8B7ynqrt9JJckc9bvADVSS463A16mUXQB8FDi96GniN/pTLmpv/4zKw4s/o/JHwpJi9puo3MW9C/gmlbvAl9YR89XAi4HPZubuqun1tPs6YBWV0ojnAP9VNW8Dla7wHoiIh8b5bI3eB5JENOe5FUmSJKkzeQdZkiRJqmKCLEmSJFUxQZYkSZKqmCBLkiRJVeaUHcBk5s+fn4sWLSo7DEmSJHWB22677aHMPHiq5do6QV60aBEbN24sOwxJkiR1gYjYPvVSllhIkiRJezFBliRJkqq0tMQiIv4OeC2QwBbgNZn561bGIEndYv3AIGv7t7FjaJgF83pZtXwxK5b2lR2WJHW8liXIEdEH/C1wTGYOR8T1VIYDvbxVMUhSmRqZ0K4fGGT1ui0M7xoBYHBomNXrtgDMqE0TbklqfYnFHKA3IuYA+wE7Wrx9SSrFaEI7ODRM8mRCu35gcFrtre3f9kRyPGp41whr+7e1RXyS1MlaliBn5iDwQeAe4H7gF5n5lbHLRcS5EbExIjbu3LmzVeFJUlM1OqHdMTRc1/SpNDo+SepkrSyxeCrwCuBIYAj4bEScnZlXVi+XmZcAlwAsW7YsWxWfJI3VyJKDRie0C+b1MjjOugvm9U6rvUbHN8qyDUmdqJUlFi8GfpyZOzNzF7AO+MMWbl+SatbokoOJEtfpJrSrli+md27PXtN65/awavniabXX6PjAsg1JnauVCfI9wAsiYr+ICOAUYGsLty9JNWt0yUGjE9oVS/tYc9qx7NNTOY33zetlzWnHTvvubKPjA8s2JHWulpVYZOYtEXEDsAnYDQxQlFJIUrtpdMnBaOJ6/g2beXxkD30NKDdYsbSPa269B4DrVp4w7XaaFV+zyjYkqdla2g9yZl4EXNTKbUrSdDS6xhcam9A2Q6Pja8bPUJJawZH0JHWN9QODnHjxBo688IucePGGGdW6NqPkYLZpxs+wkftYkibS0jvIktQsjR44oxklB7NNo3+GzRgcRZLGY4IsqStM9kDYdJOndi+J6ASN/Bk2Yx9L0ngssZDUFXwgrPu5jyW1igmypK7QjH581V7cx5JaxQRZUlfwobru5z6W1CrWIEvqCj5U1/3cx5JaxQRZUtfwobru5z6W1AomyJJKs35gkLX929gxNMwC7waqBP4OShqPCbKkUtinrcrm76CkifiQnqRSTNanrdQK/g5KmogJsqRS2KetyubvoKSJmCBLKoV92qps/g5KmogJsqRS2KetyubvoKSJ+JCepFLYp63K5u+gpImYIEsqjX3aqmz+DkoajyUWkiRJUhUTZEmSJKmKJRaSauaoY9LkPEak7jBlghwRB9XQzp7MHGpAPJLalKOOSZPzGJG6Ry13kHcUXzHJMj3AwoZEJKktTTbqmBd/yWNE6ia1JMhbM3PpZAtExECD4pHUphx1TJqcx4jUPWp5SK+Wfm/sG0fqco46Jk3OY0TqHlMmyJn5a4CIeO/YeRHRU72MpO7lqGPS5DxGpO5RTzdvfRFx5uibiHg6cFPjQ5LUjlYs7WPNaceyT0/ltNE3r5c1px1rbaVU8BiRukc93bytBPoj4kdAApcBF9SzsYiYB3wS+P2ijb/JzG/V04ak8jjqmDQ5jxGpO9TSzdsVwCZgAHgDcDWwG1iRmXfWub2PAl/OzNMjYh9gvzrXlyRJkpqqlhKLy6l08fYa4EpgEfBz4OyIOL3WDUXEgcCLgE8BZObj9p0sSZKkdjPlHeTM3ABsGH0fEXOAo4ElwPHADTVu60hgJ3BZRCwBbgPOy8xfVi8UEecC5wIsXGjXypIkSWqtuoeazszdwJbi68o6t/U84E2ZeUtEfBS4EHjHmPYvAS4BWLZsWdYbn6QnOeyt1Nk8hqVyTFliERGbGrEMcB9wX2beUry/gUrCLKkJRoe9HRwaJnly2Nv1A4NlhyapBh7DUnlqqUE+OiI2T/K1BZg/VSOZ+QBwb0SMdgh5CnDHDGKXNInJhr2V1P48hqXy1FJi8ewalhmZehEA3gRcVfRgcReVB/8kNYHD3kqdzWNYKk8tD+ltb9TGMvN2YFmj2pM0sQXzehkc50LqsLdSZ/AYlspTz0h6AETE9RHx2YhYGxFnVpVMSGojDnsrdTaPYak80+nF4gyAiHgG8BbgX4DfaXBckmZo9En382/YzOMje+jzCXipo3gMS+WpO0GOiBcDp1IZBe8bwNsbHZSkxnDYW6mzeQxL5ai7xAK4FNgf+Dpwa2b+orEhSZIkSeWZTonFwog4DPhvVIabflZmntn40CRJkqTWqzlBjoiTgb8ChoDvAZuBL2fmY02KTZIkSWq5eu4gXwq8GZgLPBdYATwHeGYT4pIkSZJKUU+CvD0z1xevP9uMYCRJkqSyTfmQXkRcERFvBr4dEW9pQUySJElSaWq5g3w5sAQ4BPjTiDgP+G7xtTkzvZssNcj6gUHW9m9jx9AwC+zzVFKDeY6RalPLUNMbgA2j7yNiDnA0laT5D7DcQmqI9QODrF63heFdIwAMDg2zet0WAC9gkmbMc4xUu7r7Qc7M3Zm5JTOvzMzzmxGUNBut7d/2xIVr1PCuEdb2byspIkndxHOMVLvpDBQiqQl2DA3XNV2S6uE5RqpdTQlyRJxVfH9lc8ORZq8F83rrmi5J9fAcI9Wu1jvIfRFxBnBYM4ORZrNVyxfTO7dnr2m9c3tYtXxxSRFJ6iaeY6Ta1dLN20XAQcBVwEER8c6mRyXNQiuW9rHmtGPZp6dyWPbN62XNacf68IykhvAcI9Wull4s3h0RbwPOAfoy8x+bH5Y0O61Y2sc1t94DwHUrTyg5GkndxnOMVJtaSyzuz8xrgR3NDEaSJEkqW00JcmZeVXy/prnhSJIkSeWymzdJkiSpypQ1yBFxUA3t7MnMoQbEI0mSJJVqygSZSt3xDiAmWaYHWNiQiCRJkqQS1ZIgb83MpZMtEBEDDYpHkiRJKlUtNcgnAETEe8fOiIie6mUkSZKkTjdlgpyZvy5e9o0OOQ0QEU8HbhqzjCRJktTRaimxGLUS6I+IO4EELgMuqHeDxV3njcBgZr6s3vWldrF+YJC1/dvYMTTMgnm9rFq+2BGpJM0qngfVrWrpxeIKYBMwALwBuBrYDazIzDunsc3zgK3A70xjXaktrB8YZPW6LQzvGgFgcGiY1eu2AHhxkDQreB5UN6ulBvlyKj1YvAa4ElgE/Bw4OyJOr2djEXEYcCrwybqilNrM2v5tT1wURg3vGmFt/7aSIpKk1vI8qG425R3kzNwAbBh9HxFzgKOBJcDxwA11bO8jwPnAUyZaICLOBc4FWLjQnuPUnnYMDdc1XZK6jedBdbO6R9LLzN2ZuSUzr8zMVbWuFxEvAx7MzNumaP+SzFyWmcsOPvjgesOTWmLBvN66pktSt/E8qG42ZYIcEZsasQxwIvDyiLgbuBY4OSKurGE9qe2sWr6Y3rk9e03rndvDquWLS4pIklrL86C6WS29WBwdEZsnmR/AgVM1kpmrgdUAEXES8LbMPLuWIKV2M/oAyvk3bObxkT30+fS2pFnG86C6WS0J8rNrWGZk6kWk7rJiaR/X3HoPANetdKwcSbOP50F1q1oe0tve6I1m5teArzW6XUmSJGmm6n5Ib6yI+KNGBCJJkiS1gxknyMBfNKANSZIkqS3UM9Q0ABFxI/BjKqPr3TadNiRJkqR2VfMd5Ij4aEREZr4c+BDwMPCXwBHNCk6SJElqtXpKLB4BboyI/YoH934FnJyZpzYnNEmSJKn1ai6PyMx/iIizgK9HxOPAo8CFTYtMkiRJKkHNCXJEnAK8DvglcCjwN5m5rVmBSZIkSWWo5wG7vwfekZnfjIhjgesi4i2ZuaFJsUkNtX5gkLX929gxNMwCR3ySpLbjeVrtop4Si5OrXm+JiJcCnwP+sBmBSY20fmCQ1eu2MLyrMujj4NAwq9dtAfDkK0ltwPO02sm0+0HOzPuBUxoYi9Q0a/u3PXHSHTW8a4S1/VYJSVI78DytdjKjgUIyc7hRgUjNtGNo/F/ViaZLklrL87TaSSNG0pPa3oJ5vXVNlyS1ludptZO6EuSIOLn6u9QpVi1fTO/cnr2m9c7tYdXyxSVFJEmq5nla7aTeO8gfHPNd6ggrlvax5rRj2aen8ivfN6+XNacd64MfktQmPE+rndTTzVu1aGgUUgusWNrHNbfeA8B1K08oORpJ0liep9UurEGWJEmSqpggS5IkSVVMkCVJkqQq9SbIjxbfH2l0IJIkSVI7qCtBzswXVX+XJEmSuo0lFpIkSVKV6XbzJjXd+oFB1vZvY8fQMAvm9bJq+WL7w5Qk1czriKZrygQ5Ig6qoZ09mTnUgHgkoHJSW71uC8O7RgAYHBpm9botAJ7cJElT8jqimajlDvKO4muywUF6gIUNiUgC1vZve+KkNmp41whr+7d5YpMkTcnriGailgR5a2YunWyBiBhoUDwSADuGhuuaLklSNa8jmolaHtI7ASAi3jt2RkT0VC8zmYg4PCJujog7IuL7EXFefaFqNlkwr7eu6ZIkVfM6opmYMkHOzF8XL/si4qzR6RHxdOCmMctMZjfw1sw8BngB8IaIOKb+kDUbrFq+mN65PXtN653bw6rli0uKSJLUSbyOaCbq6cViJdAfEXcCCVwGXFDrypl5P3B/8fqRiNgK9AF31BGDZonR+rDzb9jM4yN76PPpY0lSHbyOaCZq6cXiCmATMAC8Abiayt3gFZl553Q2GhGLgKXALdNZX7PDiqV9XHPrPQBct3LKKh5JkvbidUTTVUsN8uVUerB4DXAlsAj4OXB2RJxe7wYj4gDgc8CbM/PhceafGxEbI2Ljzp07621ekiRJmpEp7yBn5gZgw+j7iJgDHA0sAY4Hbqh1YxExl0pyfFVmrptge5cAlwAsW7Ysa21bkiRJaoS6R9LLzN3AluLrylrXi4gAPkWl27gP1btdSZIkqRWmLLGIiE2NWAY4ETgHODkibi++/qyG9SRJkqSWqeUO8tERsXmS+QEcOFUjmflNJh+NTx3OMe8lSd3Oa93sUEuC/OwalhmZehF1M8e8lyR1O691s0ctA4Vsz8ztwCeBeaPvx3zd1/xQ1c4mG/NekqRu4LVu9qilm7dRFwAfiYjLIuLQZgWkzuSY95Kkbue1bvaoOUHOzE2Z+cfAF4AvR8RFEeGA5gIc816S1P281s0e9dxBHu2qbRvwceBNwA8j4pxmBKbO4pj3kqRu57Vu9qg5QY6I/wQGgQ8DfcCrgZOA50fEJc0ITp1jxdI+1px2LPv0VH6l+ub1sua0Y31oQZLUNbzWzR71DBRyLnBHZo4d3e5NEbG1gTGpQznmvSSp23mtmx1qTpAz8/uTzD61AbFIkiRJpaurBnkimXlXI9qRJEmSylZPiYW6jKMBSZJULq/F7ckEeZZyNCBJksrltbh9NaTEQp3H0YAkSSqX1+L2ZYI8SzkakCRJ5fJa3L5MkGcpRwOSJKlcXovblwnyLOVoQJIklctrcfvyIb1ZarT4//wbNvP4yB76fHJWkqSW8lrcvkyQO0iju4JxNCBJksrV6Gux3cY1hglyh7ArGEmSNBlzhcaxBrlD2BWMJEmajLlC45ggdwi7gpEkSZMxV2gcE+QOYVcwkiRpMuYKjWOC3CHsCkaSJE3GXKFxfEiviRr5JKldwUiSpMk0I1eYrb1imCA3STOeJLVbNkmSNJlG5gqzuVcMSyyaxCdJJUlSJ5vNuYx3kKs08t8IPkkqSZI6WTNymU4p2WjpHeSIeElEbIuIOyPiwlZueyqj/0YYHBomefLfCOsHBqfVnk+SSpKkTtboXKbRuVYztSxBjoge4GPAS4FjgDMj4phWbX8qjf43gk+SSpKkTtboXKaTSjYiM1uzoYgTgHdl5vLi/WqAzFwz0TrLli3LjRs3tiS+Iy/8Igms3Px5jvrF3n/JvOCop02rzYcefYwf7fwlmcm+c3o4/KBe5h+w74zivOP+hwE45tDfmVE7ttee7TWjTduzvTLba0abtmd7ZbbXjDbbub1G5jLfvuunT7y+68A+/uW5rwAggB9ffOqMY61FRNyWmcumWq6VNch9wL1V7+8Djh+7UEScC5wLsHDhwtZERuXfBYPj1NTsO6dnnKVrM/+AfXn0sd0ALHra/tNup9p++0w/Httr//aa0abt2V6Z7TWjTduzvTLba0ab7dxeI3OZfef08Njukd+Y3o7lp628g3w68JLMfG3x/hzg+Mx840TrtPIO8tiuTKDyb4Q1px3blsXjkiRJnaQdcq12vIM8CBxe9f6wYlpbGN0xnfBkpSRJUqfppFyrlXeQ5wA/AE6hkhh/BzgrM78/0TqtvIMsSZKk7tZ2d5Azc3dEvBHoB3qASydLjiVJkqQytHSgkMz8EvClVm5TkiRJqodDTUuSJElVWlaDPB0RsRPYXsKm5wMPlbBdjc/90X7cJ+3F/dF+3Cftxf3RfsraJ0dk5sFTLdTWCXJZImJjLQXcag33R/txn7QX90f7cZ+0F/dH+2n3fWKJhSRJklTFBFmSJEmqYoI8vkvKDkB7cX+0H/dJe3F/tB/3SXtxf7Sftt4n1iBLkiRJVbyDLEmSJFUxQZYkSZKqmCBXiYiXRMS2iLgzIi4sOx5BRNwdEVsi4vaI2Fh2PLNRRFwaEQ9GxPeqph0UEV+NiB8W359aZoyzyQT7410RMVgcJ7dHxJ+VGeNsEhGHR8TNEXFHRHw/Is4rpnuMlGSSfeJxUoKI+O2IuDUivlvsj3cX04+MiFuKnOu6iNin7FirWYNciIge4AfAnwD3Ad8BzszMO0oNbJaLiLuBZZlpB+8liYgXAY8CV2Tm7xfTPgD8LDMvLv6YfGpmXlBmnLPFBPvjXcCjmfnBMmObjSLiUODQzNwUEU8BbgNWAK/GY6QUk+yTM/A4abmICGD/zHw0IuYC3wTOA94CrMvMayPiE8B3M/PjZcZazTvIT3o+cGdm3pWZjwPXAq8oOSapdJn5H8DPxkx+BfDp4vWnqVx81AIT7A+VJDPvz8xNxetHgK1AHx4jpZlkn6gEWfFo8XZu8ZXAycANxfS2O0ZMkJ/UB9xb9f4+PKDaQQJfiYjbIuLcsoPREw7JzPuL1w8Ah5QZjAB4Y0RsLkow/Hd+CSJiEbAUuAWPkbYwZp+Ax0kpIqInIm4HHgS+CvwIGMrM3cUibZdzmSCr3b0wM58HvBR4Q/HvZbWRrNRpWatVro8DzwCOA+4H/rHccGafiDgA+Bzw5sx8uHqex0g5xtknHiclycyRzDwOOIzKf+yfXXJIUzJBftIgcHjV+8OKaSpRZg4W3x8E/o3KgaXy/aSo8xut93uw5Hhmtcz8SXEB2gP8Kx4nLVXUVX4OuCoz1xWTPUZKNN4+8TgpX2YOATcDJwDzImJOMavtci4T5Cd9B3hW8VTlPsArgRtLjmlWi4j9iwcsiIj9gT8Fvjf5WmqRG4FXFa9fBXy+xFhmvdFErPA/8DhpmeIBpE8BWzPzQ1WzPEZKMtE+8TgpR0QcHBHzite9VDpD2EolUT69WKztjhF7sahSdPnyEaAHuDQz31dySLNaRBxF5a4xwBzgavdJ60XENcBJwHzgJ8BFwHrgemAhsB04IzN9cKwFJtgfJ1H5t3ECdwMrq+pf1UQR8ULgG8AWYE8x+e1Ual49RkowyT45E4+TlouI51J5CK+Hyo3Z6zPzPcU1/lrgIGAAODszHysv0r2ZIEuSJElVLLGQJEmSqpggS5IkSVVMkCVJkqQqJsiSJElSFRNkSZIkqYoJsiRJklTFBFmSJEmqYoIsSR0qIk6JiM+UHYckdRsTZEnqXEuojEAlSWogE2RJ6lxLgIGI2DciLo+I90dElB2UJHW6OWUHIEmatucCDwL9wCcz88qS45GkrhCZWXYMkqQ6RcRc4CFgO7AyM79VckiS1DUssZCkznQ08B1gNzBSciyS1FVMkCWpMy0B/gt4JXBZRBxScjyS1DVMkCWpMy0BvpeZPwAuAK4vyi4kSTNkDbIkSZJUxTvIkiRJUhUTZEmSJKmKCbIkSZJUxQRZkiRJqmKCLEmSJFUxQZYkSZKqmCBLkiRJVf5/7Z0Km+A9G/4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x112fa4828>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "L = 16  # length of signal x[k]\n",
    "N = 16  # length of signal h[k]\n",
    "M = N+L-1\n",
    "\n",
    "# generate signals\n",
    "x = np.ones(L)\n",
    "h = sig.triang(N)\n",
    "\n",
    "# linear convolution\n",
    "y1 = np.convolve(x, h, 'full')\n",
    "# fast convolution\n",
    "y2 = np.fft.ifft(np.fft.fft(x, M) * np.fft.fft(h, M))\n",
    "\n",
    "plt.figure(figsize=(10, 6))\n",
    "plt.subplot(211)\n",
    "plt.stem(y1)\n",
    "plt.xlabel(r'$k$')\n",
    "plt.ylabel(r'$y[k] = x_L[k] * h_N[k]$')\n",
    "plt.title('Result of linear convolution')\n",
    "\n",
    "plt.subplot(212)\n",
    "plt.stem(y1)\n",
    "plt.xlabel(r'$k$')\n",
    "plt.ylabel(r'$y[k] = x_L[k] * h_N[k]$')\n",
    "plt.title('Result of fast convolution')\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Example - Numerical complexity\n",
    "\n",
    "It was already argued that the numerical complexity of the fast convolution is considerably lower due to the usage of the FFT. The gain with respect to the convolution is evaluated in the following. In order to measure the execution times for both algorithms the `timeit` module is used. The algorithms are evaluated for the convolution of two random signals $x_L[k]$ and $h_N[k]$ of length $L=N=2^n$ for $n=0, 1, \\dots, 16$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA4sAAAJgCAYAAAAj5/aWAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJzs3XuYXWV99//3RwIhBCUKFCNRg4/YFh1ABcXziFYtWNEefLSRitVGW6Vaoxb7tNXa2lIVq9b6q1FRqigq4hFPVN14RgMGwskTjkJAEXSQYItO+P7+WGt0Z/ZMMpDZsyd73q/rmmv2ute91v6uve9s5sO91tqpKiRJkiRJ6na7QRcgSZIkSVp4DIuSJEmSpB6GRUmSJElSD8OiJEmSJKmHYVGSJEmS1MOwKEmSJEnqYViUJM0oyZokn57n5xxL8uj5fM5BStJJ8qz28by/3lNq+ack1yX54aBqGISdGXNJ7pZkS5Ld5rouSRo0w6IkzYMkf5xkQ/tH5TVJPpHkoYOua0eq6vSqesyg65itXT1o9uv1TrI6SSVZMs26NydZm+RuwDrgkKq6804812iSq3am3oVs6hirqh9U1d5VtXWQdUlSPxgWJanPkrwQeB3wz8ABwN2ANwHHDbKuHZkuWGhw+vh+/C7wcZpxeX1VXdun55Ek7WIMi5LUR0n2AV4BPLeqzqqqm6rql1X10ap6cdtnaZLXJbm6/XldkqXtutEkVyV5SZJr21nJJyY5Jsm3kvwkyd90Pd/Lk5yZ5L1JbkxyQZLDutaflOS77bpLkzypa90JSb6U5N+SXA+8vG37Yrs+7bprk/wsyaYk95k8ziT/leTHSb6f5G+T3K5rv19M8pokP03yvSS/u4OX7si2vp8meXuSPbvqfHySjUnGk3w5yaFt+ztpAs9H2xnclyQ5Lcm6dv2B7ezac9vl/9O+frfb3n7bdXdJ8oH2+L6X5C+nvObva4//xiSXJDliO2Pid5JcnuSGJG8EMuU9+GLXciV5bpJvA99u234ryTlt7d9M8uSu/suSnNK+Bze0r/sy4PNtl/H2tXlQ2/9QYBz4LeAc4C7t+ne069+f5Iftvj6f5N5dz3VM+x7dmGRzkhclWQ58oms/W5LcZZrXYKY6SfKE9jUcT3OK7m93bTfWPs9F7XbvnRwbSS5L8viuvkva9+t+O9rvlNrekeSfupZ/NVM6wxjbZta2HSsfad+f7yT5s6593aqxIkmDZliUpP56ELAn8MHt9Pl/wFHA4cBhwAOAv+1af+d2HwcCfw+8BXgacH/gYcDfJTmoq/9xwPuBOwHvBj6UZPd23XfbbfYB/gF4V5KVXds+ELiCZgb0lVPqfAzwcOBe7fZPBq5v1/1723YP4BHAnwDPmLLfbwL7Aa8C3pYkzGwN8Fjg/7TP97cASe4LnAo8G9gXeDPwkSRLq+p44AfA77WnBb4KOBcYbff5iPbYHt61/IWqumV7+23D5EeBC2neg0cBL0jy2K56nwCcAawAPgK8cbqDSrIfcFZ7PPvRvB8P2c7rAPBEmtfvkDaMnUPzvv4G8BTgTUkOafu+hmZcPJjm/X8JcEvXMa9oX5uvtMvHAGdX1X/TzDBe3a4/oV3/CeDg9rkuAE7vquttwLOr6vbAfYDPVtVNU/azd1VdPc0xTVtnknsB7wFeAOxPM+P50SR7dG37ZOBxwEHAocBkre8BntrV77HAdVV1wSz3u0MzjLGpzgCuAu4C/CHwz0mO7lo/q7EiSQuBYVGS+mtfmj9YJ7bTZw3wiqq6tqp+TBPiju9a/0vglVX1S5o/MvcDXl9VN1bVJcClNCFz0vlVdWbb/7U0QfMogKp6f1VdXVW3VNV7aWarHtC17dVV9e9VNVFV/zOlzl8Ct6eZhUpVXVZV16S5scdTgJe2NY0Bp0w5hu9X1Vva67pOA1bSBNKZvLGqrqyqn9CE1skQsBZ4c1WdV1Vbq+o04ObJ45vGucBD28D3cJqgOhnOHtGu39F+jwT2r6pXVNUvquoKmsD+lK7n+WJVfbw9vney7fvR7Rjgkq7353XAjm4m8y9V9ZP2/Xg8MFZVb2/fo28AHwD+qD3GPwWeX1Wb2+P4clXdvJ19H0sTnKZVVae27+nNwMuBw9LMlkMzHg5Jcoeq+mlVXbCD4wBgB3X+X5rwek77+rwGWEYTKie9oR3DP6EJ8Ye37e8GnpBkr3b5j2kCIrPc705Lclea8fXXVfW/VbUReCvN/zyZNNuxIkkDZ1iUpP66Htgv27/e7C7A97uWv9+2/WofXTfPmAxwP+pa/z/A3l3LV04+qKpb+PUsB0n+JL8+1XKcZkZov+m2naqqPkszC/IfwLVJ1ie5Q7v97tMcw4Fdyz/s2s/P24fdNU/VXUf363F3YN1k/e0x3JVtX6/umr8L3EQTKB4GfAy4Oslvsm1Y3N5+705zWmX3ur9h27DbHfh+Duw5w3t+F7Z9f4rtvObTvBZ3Bx44pZY1NLPP+9H8j4Hv7mB/ACRZQRP8vzzD+t2SnJzmtOWfAWPtqsnx8gc04ff7Sc6dPLV1FrZX5zb/FtrxeyUzjCWa13rvtu93gMuA32sD4xNoAuRs9zsX7gL8pKpu7Gqb8d8C2x8rkjRwhkVJ6q+v0MxQPXE7fa6mCQGT7ta23VZ3nXzQzuKsoglId6eZEXsesG9VrQAupuuaOaC2t+OqekNV3R84hOb00BcD19HMMk09hs1zcQxs+3pcSTPLuqLrZ6+qmpxBmq7+c2lOB9yjqja3y08H7ghsnMV+rwS+N2Xd7avqmNtwXNew7fuTKcc6ne5juhI4d0ote1fVn9O8D/9Lc+ru9vYx6bE0p47OdBfPP6Y5pfnRNKcYr54sG6Cqvl5Vx9Gcovoh4H3bea5u26tzm38LXa/PbMfS5KmoxwGXtgHy1u73JmCvruWpd4bd3vFdDdwpye272nb234IkDYxhUZL6qKpuoLnO8D/S3JhmryS7J/ndJJPXO70H+Nsk+7fXtP098K6deNr7J/n9drbiBTRh9avAcpo/dH8MkOQZNDOLs5LkyCQPbK9/vInmD/5b2rDxPuCVSW7fhtIX7uQxPDfJqiR3ormm871t+1uA57R1JMnyJMd2/XH+I5rrJrudSxOQJ2/y0mmXv9gVlLa3368BNyb56zQ3ZtktyX2SHHkbjuts4N5d789f0htGtudjwL2SHN+Oo93b9+W329myU4HXtjdZ2S3Jg9LcLOnHNNcudr82x7T1zOT2NGPneprw9M+TK5LskeY7IfdpT+v8Wbt/aN6DfbtOV93GDup8H3Bskke142xdW8O0s5/TOIPm2to/59ezitzK/W4EjklypyR3pvk31G26MTZ5bFe2+/yXJHumuYHQM9m5fwuSNDCGRUnqs6o6hSY8/S3NH+1X0oSVD7Vd/gnYAFwEbKK5kcg/9e5p1j5Mc43WT2muG/z9au7AeinNtYRfofmDdwT40q3Y7x1oQtVPaU6tux54dbvuRJoAeQXwRZo/1E/diWN4N/Dpdn/fpX09qmoD8Gc0p8P+FPgOv77BCcC/0ATv8SQvatvOpQk+k2HxizThZ3J5u/ttA+XjaU5l/R7NzNhbaWbbbpWqug74I+BkmtfvYG7Fe9Ce3vgYmuslr6Y5pfFfgaVtlxfRjKGvAz9p192uPfX3lcCX2tfmQTQzi5/cztP9F837vJnmutivTll/PDDWnqL6HJrTYamqy2n+B8gV7XNNd4rwTHV+k+bmTf9O8zr/Hs3NZH6x41cHquoamvH9YH79Pxi4lft9J83NjMZoxuB7p6yfbox1eyrNLOzVNDe2elk1NxCSpF1OmsslJEnDIMnLgXtW1dMGXYsWriQPoLmJ0AN22FmStGg5syhJ0uL0skEXIEla2Lz7liRJi0xVfW3QNUiSFj5PQ5UkSZIk9fA0VEmSJElSj0V3Gup+++1Xq1evHnQZ0jZuuukmli9fPugypDnluNYwclxrGDmuF5/zzz//uqraf0f9Fl1YXL16NRs2bBh0GdI2Op0Oo6Ojgy5DmlOOaw0jx7WGkeN68Uny/dn08zRUSZIkSVIPw6IkSZIkqYdhUZIkSZLUw7AoSZIkSephWJQkSZIk9TAsSpIkSZJ6GBYlSZIkST0Mi5IkSZKkHoZFSZIkSVIPw6IkSZIkqYdhUZIkSZLUw7AoSZIkSephWJQkSZIk9TAsSpIkSZJ6GBYlSZIkST0Mi5IkSZKkHoZFSZIkSVIPw6IkSZIkqYdhUZIkSZLUw7AoSZIkSephWJQkSZIk9TAsSpIkSZJ6LBl0AZKkIfSCF3DPq66C0dFBVyJJkm4jw6Ikae5t3Mje4+ODrkKSJO0ET0OVJEmSJPUwLEqSJEmSehgWJUmSJEk9DIuSJEmSpB6GRUmSJElSD8OiJEmSJKmHYVGSJEmS1MOwKEmSJEnqYViUJEmSJPUwLEqSJEmSehgWJUmSJEk9DIuSJEmSpB6GRUmSJElSD8OiJEmSJKmHYVGSJEmS1MOwKEmSJEnqkaoadA3zaunKg2vl01836DKkbawbmeCUTUsGXYY0Z85490msWl489Lh/HXQp0pzy81rDyHHdP2MnHzvoEqaV5PyqOmJH/ZxZlCRJkiT1MCxKkiRJknoYFiVJkiRJPQyLkiRJkqQehkVJkiRJUg/DoiRJkiSpR1/DYpIVSc5McnmSy5I8KMnLk2xOsrH9OabtuzrJ/3S1/2fXfj6Z5MIklyT5zyS7da07sd3/JUle1c/jkSRJkqTFot9fqPJ64JNV9YdJ9gD2Ah4L/FtVvWaa/t+tqsOnaX9yVf0sSYAzgT8CzkjySOA44LCqujnJb/TpOCRJkiRpUelbWEyyD/Bw4ASAqvoF8Ism7906VfWz9uESYA+g2uU/B06uqpvbftfuXNWSJEmSJOjvaagHAT8G3p7kG0nemmR5u+55SS5KcmqSO3Zv0/Y9N8nDuneW5FPAtcCNNLOLAPcCHpbkvHabI/t4PJIkSZK0aPTzNNQlwP2AE6vqvCSvB04C3gj8I83s4D8CpwB/ClwD3K2qrk9yf+BDSe49OatYVY9NsidwOnA0cE77HHcCjgKOBN6X5B5VVd2FJFkLrAVYse/+rBuZ6ONhS7feActwXGqorFpe7LGb41rDx89rDSPHdf90Op1Bl7BT+hkWrwKuqqrz2uUzgZOq6keTHZK8BfgYQHsq6eTppOcn+S7NzOGGyf5V9b9JPkxzneI57XOc1YbDryW5BdiPZkaTru3WA+sBlq48uE7Z1O9LNaVbZ93IBI5LDZMjbwqrlpfjWkPHz2sNI8d1/4ytGR10CTulb6ehVtUPgSuT/Gbb9Cjg0iQru7o9CbgYIMn+k3c5TXIP4GDgiiR7T26TZAlwLHB5u/2HgEe26+5Fcz3jdf06JkmSJElaLPr9vxBOBE5v74R6BfAM4A1JDqc5DXUMeHbb9+HAK5L8ErgFeE5V/STJAcBHkiylCbefAya/VuNU4NQkFwO/AJ4+9RRUSZIkSdKt19ewWFUbgSOmNB8/Q98PAB+Ypv1HNNcjTrfNL4Cn7WSZkiRJkqQp+nk3VEmSJEnSLsqwKEmSJEnqYViUJEmSJPUwLEqSJEmSehgWJUmSJEk9DIuSJEmSpB5ZbF9LeMQRR9SGDRsGXYa0jU6nw+jo6KDLkObO6Cjj4+Os2Lhx0JVIc8rPaw0jx/Xik+T8qpr6FYc9nFmUJEmSJPUwLEqSJEmSehgWJUmSJEk9DIuSJEmSpB6L7gY3S1ceXCuf/rpBlyFtY93IBKdsWjLoMqQ5c8a7T2LV8uKhx/3roEuR5pSf1xpG8zGux04+tq/7163jDW4kSZIkSbeZYVGSJEmS1MOwKEmSJEnqYViUJEmSJPUwLEqSJEmSehgWJUmSJEk9DIuSJEmSpB59C4tJ9kzytSQXJrkkyT+07acn+WaSi5OcmmT3tv3FSTa2Pxcn2ZrkTu26FUnOTHJ5ksuSPKhtPzzJV9ttNiR5QL+OR5IkSZIWk37OLN4MHF1VhwGHA49LchRwOvBbwAiwDHgWQFW9uqoOr6rDgZcC51bVT9p9vR74ZFX9FnAYcFnb/irgH9pt/r5dliRJkiTtpCX92nFVFbClXdy9/amq+vhknyRfA1ZNs/lTgfe0ffYBHg6c0O73F8AvJp8GuEP7eB/g6jk9CEmSJElapPoWFgGS7AacD9wT+I+qOq9r3e7A8cDzp2yzF/A44Hlt00HAj4G3Jzms3d/zq+om4AXAp5K8hmaW9MEz1LEWWAuwYt/9WTcyMWfHKM2FA5bhuNRQWbW82GM3x7WGj5/XGkbzMa47nU5f96/+6GtYrKqtwOFJVgAfTHKfqrq4Xf0m4PNV9YUpm/0e8KWuU1CXAPcDTqyq85K8HjgJ+Dvgz4G/qqoPJHky8Dbg0dPUsR5YD7B05cF1yqa+HrZ0q60bmcBxqWFy5E1h1fJyXGvo+HmtYTQf43pszWhf96/+mJe7oVbVOPA5mhlDkrwM2B944TTdn0J7CmrrKuCqrlnJM2nCI8DTgbPax+8HvMGNJEmSJM2Bft4Ndf92RpEky4DfAS5P8izgscBTq+qWKdvsAzwC+PBkW1X9ELgyyW+2TY8CLm0fX932Bzga+HafDkeSJEmSFpV+zjevBE5rr1u8HfC+qvpYkgng+8BXkgCcVVWvaLd5EvDp9nrEbicCpyfZA7gCeEbb/mfA65MsAf6X9rpESZIkSdLO6efdUC8C7jtN+4zPWVXvAN4xTftG4Ihp2r8I3H9n6pQkSZIk9ZqXaxYlSZIkSbsWw6IkSZIkqYdhUZIkSZLUw7AoSZIkSephWJQkSZIk9UhVDbqGeXXEEUfUhg0bBl2GtI1Op8Po6Oigy5Dmzugo4+PjrNi4cdCVSHPKz2sNI8f14pPk/Krq+baJqZxZlCRJkiT1MCxKkiRJknoYFiVJkiRJPQyLkiRJkqQeSwZdwHzbtPkGVp909qDLkLaxbmSCExyXGiJnXHE9q5YXhzuuNWT8vN61jZ187KBLkHYpzixKkiRJknoYFiVJkiRJPQyLkiRJkqQehkVJkiRJUg/DoiRJkiSph2FRkiRJktTDsChJkiRJ6tHXsJjk1CTXJrl4SvuJSS5PckmSV7VtD0iysf25MMmTuvr/Vdv34iTvSbLnlP29IcmWfh6LJEmSJC0m/Z5ZfAfwuO6GJI8EjgMOq6p7A69pV10MHFFVh7fbvDnJkiQHAn/ZrrsPsBvwlK79HQHcsc/HIUmSJEmLSl/DYlV9HvjJlOY/B06uqpvbPte2v39eVRNtnz2B6tpmCbAsyRJgL+BqgCS7Aa8GXtK3g5AkSZKkRWjJAJ7zXsDDkrwS+F/gRVX1dYAkDwROBe4OHN+Gx81JXgP8APgf4NNV9el2X88DPlJV1ySZ8QmTrAXWAqzYd3/WjUzM2FcahAOW4bjUUFm1vNhjN8e1ho+f17u2Tqcz6BIWpC1btvjaaFqDCItLgDsBRwFHAu9Lco9qnAfcO8lvA6cl+QSwjOa01YOAceD9SZ4GfBb4I2B0R09YVeuB9QBLVx5cp2waxGFLM1s3MoHjUsPkyJvCquXluNbQ8fN61za2ZnTQJSxInU6H0dHRQZehBWgQn3ZXAWdVVQFfS3ILsB/w48kOVXVZe8Oa+9CExO9V1Y8BkpwFPBj4KXBP4DvtrOJeSb5TVfec16ORJEmSpCE0iK/O+BDwSIAk9wL2AK5LclB7TSJJ7g78FjBGc/rpUUn2SpMKHwVcVlVnV9Wdq2p1Va0Gfm5QlCRJkqS50deZxSTvoTlNdL8kVwEvo7km8dT26zR+ATy9qirJQ4GTkvwSuAX4i6q6jiZInglcAEwA36A9pVSSJEmS1B99DYtV9dQZVj1tmr7vBN45w35eRhM0t/dce9/qAiVJkiRJ0xrEaaiSJEmSpAXOsChJkiRJ6mFYlCRJkiT1MCxKkiRJknoYFiVJkiRJPfp6N9SFaOTAfdhw8rGDLkPaRqfTYWzN6KDLkObOV1/N+Pg4Y37easj4eS1pMXFmUZIkSZLUw7AoSZIkSephWJQkSZIk9TAsSpIkSZJ6GBYlSZIkST0W3d1QN22+gdUnnT3oMqRtrBuZ4ATHpYbIGVdcz6rlxeGOaw0ZP693zLsgS8PDmUVJkiRJUg/DoiRJkiSph2FRkiRJktTDsChJkiRJ6mFYlCRJkiT1MCxKkiRJknoYFiVJkiRJPQYSFpP8VZJLklyc5D1J9kzyvCTfSVJJ9uvqe8ckH0xyUZKvJblP237XJJ9Lcmm7r+cP4lgkSZIkaRjNe1hMciDwl8ARVXUfYDfgKcCXgEcD35+yyd8AG6vqUOBPgNe37RPAuqo6BDgKeG6SQ+bhECRJkiRp6A3qNNQlwLIkS4C9gKur6htVNTZN30OAzwJU1eXA6iQHVNU1VXVB234jcBlw4LxUL0mSJElDbsl8P2FVbU7yGuAHwP8An66qT29nkwuB3we+kOQBwN2BVcCPJjskWQ3cFzhvuh0kWQusBVix7/6sG5nY+QOR5tABy3BcaqisWl7ssZvjWsPHz+sd63Q6gy5Bt9KWLVt83zSteQ+LSe4IHAccBIwD70/ytKp61wybnAy8PslGYBPwDWBr1/72Bj4AvKCqfjbdDqpqPbAeYOnKg+uUTfN+2NJ2rRuZwHGpYXLkTWHV8nJca+j4eb1jY2tGB12CbqVOp8Po6Oigy9ACNIhPu0cD36uqHwMkOQt4MDBtWGwD4DPavgG+B1zRLu9OExRPr6qz+l+6JEmSJC0Og7hm8QfAUUn2asPfo2iuN5xWkhVJ9mgXnwV8vqp+1m77NuCyqnpt36uWJEmSpEVk3sNiVZ0HnAlcQHNa6e2A9Un+MslVNNcjXpTkre0mvw1cnOSbwO8Ck1+R8RDgeODoJBvbn2Pm81gkSZIkaVgN5KT7qnoZ8LIpzW9of6b2/Qpwr2navwikLwVKkiRJ0iI3qK/OkCRJkiQtYIZFSZIkSVIPw6IkSZIkqYdhUZIkSZLUw7AoSZIkSeoxkLuhDtLIgfuw4eRjB12GtI1Op8PYmtFBlyHNna++mvHxccb8vNWQ8fNa0mLizKIkSZIkqYdhUZIkSZLUw7AoSZIkSephWJQkSZIk9Vh0N7jZtPkGVp909qDLkLaxbmSCExyXGiJnXHE9q5YXhzuutUB58yVJ2jFnFiVJkiRJPQyLkiRJkqQehkVJkiRJUg/DoiRJkiSph2FRkiRJktTDsChJkiRJ6mFYlCRJkiT1GFhYTLJbkm8k+Vi7fFCS85J8J8l7k+zRtt89yWeSXJSkk2RV1z7uluTTSS5LcmmS1YM5GkmSJEkaLoOcWXw+cFnX8r8C/1ZV9wR+CjyzbX8N8F9VdSjwCuBfurb5L+DVVfXbwAOAa/tetSRJkiQtAgMJi+3s4LHAW9vlAEcDZ7ZdTgOe2D4+BPhs+/hzwHHtNocAS6rqHICq2lJVP5+XA5AkSZKkIbdkQM/7OuAlwO3b5X2B8aqaaJevAg5sH18I/D7weuBJwO2T7AvcCxhPchZwEPDfwElVtXXqkyVZC6wFWLHv/qwbmZjaRRqoA5bhuNRQWbW82GM3x7UWrk6nc5u227Jly23eVlqoHNeaybyHxSSPB66tqvOTjM5ikxcBb0xyAvB5YDOwlab2hwH3BX4AvBc4AXjb1B1U1XpgPcDSlQfXKZsGlZGl6a0bmcBxqWFy5E1h1fJyXGvBGlszepu263Q6jI7etm2lhcpxrZkM4r/iDwGekOQYYE/gDjSzhiuSLGlnF1fRhEKq6mqamUWS7A38QVWNJ7kK2FhVV7TrPgQcxTRhUZIkSZJ068z7NYtV9dKqWlVVq4GnAJ+tqjU01yP+Ydvt6cCHAZLsl2SyzpcCp7aPv04TMPdvl48GLp2HQ5AkSZKkobeQvmfxr4EXJvkOzTWMkzOEo8A3k3wLOAB4JUB7beKLgM8k2QQEeMt8Fy1JkiRJw2igF5NUVQfotI+voPn6i6l9zuTXd0mduu4c4ND+VShJkiRJi9NCmlmUJEmSJC0QhkVJkiRJUg/DoiRJkiSph2FRkiRJktTDsChJkiRJ6jHQu6EOwsiB+7Dh5GMHXYa0jU6nw9ia0UGXIc2dr76a8fFxxvy8lSRpl+XMoiRJkiSph2FRkiRJktTDsChJkiRJ6mFYlCRJkiT1MCxKkiRJknosuruhbtp8A6tPOnvQZUjbWDcywQmOSw2RM664nlXLi8Md1wuOd6iVJM2WM4uSJEmSpB6GRUmSJElSD8OiJEmSJKmHYVGSJEmS1MOwKEmSJEnqYViUJEmSJPUwLEqSJEmSegwkLCa5a5LPJbk0ySVJnt+23ynJOUm+3f6+45TtjkwykeQPu9pe1e7jsiRvSJL5Ph5JkiRJGjaDmlmcANZV1SHAUcBzkxwCnAR8pqoOBj7TLgOQZDfgX4FPd7U9GHgIcChwH+BI4BHzdRCSJEmSNKwGEhar6pqquqB9fCNwGXAgcBxwWtvtNOCJXZudCHwAuLZ7V8CewB7AUmB34Ed9LV6SJEmSFoElgy4gyWrgvsB5wAFVdU276ofAAW2fA4EnAY+kmT0EoKq+kuRzwDVAgDdW1WXTPMdaYC3Ain33Z93IRL8OR7pNDliG41JDZdXyYo/dHNcLUafTGXQJu7QtW7b4GmroOK41k4GGxSR708wWvqCqftZ9uWFVVZJqF18H/HVV3dLdJ8k9gd8GVrVN5yR5WFV9oft5qmo9sB5g6cqD65RNA8/I0jbWjUzguNQwOfKmsGp5Oa4XoLE1o4MuYZfW6XQYHR0ddBnSnHJcayYD+694kt1pguLpVXVW2/yjJCur6pokK/n1KadHAGe0QXE/4JgkE8DBwFeraku7z08ADwK2CYuSJEmSpFtnUHdDDfA24LKqem3Xqo8AT28fPx34MEBVHVRVq6tqNXAm8BdV9SHgB8Ajkixpw+cjaK5/lCRJkiTthEHNLD4EOB7YlGRj2/Y3wMnA+5I8E/g+8OQd7OdM4GhgE83Nbj5ZVR/tT8mSJEmStHgMJCxW1RdpbkgznUftYNsTuh5vBZ49d5VJkiRJkmBw37MoSZIkSVrADIuSJEmSpB6GRUmSJElSD8OiJEmSJKmHYVGSJEkXS+cyAAAgAElEQVSS1GNQX50xMCMH7sOGk48ddBnSNjqdDmNrRgddhjR3vvpqxsfHGfPzVpKkXZYzi5IkSZKkHoZFSZIkSVIPw6IkSZIkqYdhUZIkSZLUY9Hd4GbT5htYfdLZgy5D2sa6kQlOcFxqiIwNugBJkrTTnFmUJEmSJPUwLEqSJEmSehgWJUmSJEk9DIuSJEmSpB6GRUmSJElSD8OiJEmSJKmHYVGSJEmS1GMgYTHJqUmuTXJxV9udkpyT5Nvt7zu27WuSXJRkU5IvJzlsyr52S/KNJB+b7+OQJEmSpGE1qJnFdwCPm9J2EvCZqjoY+Ey7DPA94BFVNQL8I7B+ynbPBy7rX6mSJEmStPgMJCxW1eeBn0xpPg44rX18GvDEtu+Xq+qnbftXgVWTGyRZBRwLvLWvBUuSJEnSIrNk0AV0OaCqrmkf/xA4YJo+zwQ+0bX8OuAlwO23t+Mka4G1ACv23Z91IxM7X600hw5YhuNSQ2X8k+Ns3bqVTqcz6FKkObVlyxbHtYaO41ozWUhh8VeqqpJUd1uSR9KExYe2y48Hrq2q85OM7mB/62lPX1268uA6ZdOCPGwtYutGJnBcapicuGIF4+PjjI6ODroUaU51Oh3HtYaO41ozWUh3Q/1RkpUA7e9rJ1ckOZTmVNPjqur6tvkhwBOSjAFnAEcnedf8lixJkiRJw2khhcWPAE9vHz8d+DBAkrsBZwHHV9W3JjtX1UuralVVrQaeAny2qp42vyVLkiRJ0nAayHlvSd4DjAL7JbkKeBlwMvC+JM8Evg88ue3+98C+wJuSAExU1RHzXrQkSZIkLSIDCYtV9dQZVj1qmr7PAp61g/11gM5OFyZJkiRJAhbWaaiSJEmSpAXCsChJkiRJ6mFYlCRJkiT1MCxKkiRJknoYFiVJkiRJPQZyN9RBGjlwHzacfOygy5C20el0GFszOugypLkz+upBVyBJknaSM4uSJEmSpB6GRUmSJElSj1sVFpPcLskd+lWMJEmSJGlh2GFYTPLuJHdIshy4GLg0yYv7X5okSZIkaVBmM7N4SFX9DHgi8AngIOD4vlYlSZIkSRqo2dwNdfcku9OExTdW1S+TVJ/r6ptNm29g9UlnD7oMaRvrRiY4wXE59Ma8E7MkSdqFzGZm8c3AGLAc+HySuwM/62dRkiRJkqTB2uHMYlW9AXhDV9P3kzyyfyVJkiRJkgZtxrCY5IU72Pa1c1yLJEmSJGmB2N7M4u3nrQpJkiRJ0oIyY1isqn+Yz0IkSZIkSQvHDq9ZTLIn8Ezg3sCek+1V9ad9rEuSJEmSNECzuRvqO4E7A48FzgVWATf2syhJkiRJ0mDNJizes6r+Dripqk4DjgUe2K+Ckowl2ZRkY5INbdsfJbkkyS1Jjujq+ztJzm/7n5/k6H7VJUmSJEmLyQ5PQwV+2f4eT3If4IfAb/SvJAAeWVXXdS1fDPw+zXc+drsO+L2qurqt7VPAgX2uTZIkSZKG3mzC4vokdwT+DvgIsDfw932taoqqugwgydT2b3QtXgIsS7K0qm6ex/IkSZIkaejsMCxW1Vvbh+cC9+hvOc1TAp9OUsCbq2r9LLf7A+CC6YJikrXAWoAV++7PupGJOStWmgsHLMNxuQh0Op1BlzBvDh8fZ+vWrYvqmLU4bNmyxXGtoeO41kxmczfUpTRBbHV3/6p6RZ9qemhVbU7yG8A5SS6vqs/voMZ7A/8KPGa69W3gXA+wdOXBdcqm2UyoSvNn3cgEjsvhN7ZmdNAlzJ8VKxgfH2d0dHTQlUhzqtPpOK41dBzXmslsbnDzYeA4YAK4qeunL6pqc/v7WuCDwAO21z/Jqrbfn1TVd/tVlyRJkiQtJrOZylhVVY/reyVAkuXA7arqxvbxY4AZZzCTrADOBk6qqi/NR42SJEmStBjMZmbxy0lG+l5J4wDgi0kuBL4GnF1Vn0zypCRXAQ8Czk7yqbb/84B7An/fftXGxvb0VUmSJEnSTpjNzOJDgROSfA+4GQhQVXXoXBdTVVcAh03T/kGaU02ntv8T8E9zXYckSZIkLXazCYu/2/cqJEmSJEkLymzC4o2zbJMkSZIkDYnZXLN4AfBj4FvAt9vHY0kuSHL/fhYnSZIkSRqM2YTFc4Bjqmq/qtqX5rTUjwF/Abypn8VJkiRJkgZjNmHxqKqavPsoVfVp4EFV9VVgad8qkyRJkiQNzGyuWbwmyV8DZ7TL/xf4UZLdgFv6VlmfjBy4DxtOPnbQZUjb6HQ6jK0ZHXQZkiRJ0q/MZmbxj4FVwIfan7u1bbsBT+5faZIkSZKkQdnhzGJVXQecOMPq78xtOZIkSZKkhWDGsJjkdVX1giQfBWrq+qp6Ql8rkyRJkiQNzPZmFt/Z/n7NfBQiSZIkSVo4ZgyLVXV++/vcybYkdwTuWlUXzUNtfbFp8w2sPunsQZchbWPdyAQnLNJxOeYNpyRJkhakHd7gJkknyR2S3Am4AHhLktf2vzRJkiRJ0qDM5m6o+1TVz4DfB/6rqh4IPLq/ZUmSJEmSBmk2YXFJkpU0X5PxsT7XI0mSJElaAGYTFl8BfAr4TlV9Pck9gG/3tyxJkiRJ0iDN5nsW3w+8v2v5CuAP+lmUJEmSJGmwZjOzKEmSJElaZAyLkiRJkqQeM4bFJM9vfz9kvopJctckn0tyaZJLump4eZLNSTa2P8d0bXNokq+0/Tcl2XO+6pUkSZKkYbW9axafAbwe+HfgfvNTDhPAuqq6IMntgfOTnNOu+7eqek135yRLgHcBx1fVhUn2BX45T7VKkiRJ0tDaXli8LMm3gbskuairPUBV1aFzXUxVXQNc0z6+McllwIHb2eQxwEVVdWG7zfVzXZMkSZIkLUYzhsWqemqSO9N8bcYT5q+kRpLVwH2B84CHAM9L8ifABprZx58C9wIqyaeA/YEzqupV0+xrLbAWYMW++7NuZGJejkGarQOWsWjHZafTGXQJ6oPDx8fZunWr76+GzpYtWxzXGjqOa80kVbXjTskeNMEM4JtV1ddTPZPsDZwLvLKqzkpyAHAdUMA/Aiur6k+TvAh4LnAk8HPgM8DfVtVnZtr30pUH18qnv66f5Uu32rqRCU7ZtMNvshlKYycfO+gS1A+jo4yPj7Ni48ZBVyLNqU6nw+jo6KDLkOaU43rxSXJ+VR2xo347vBtqkkcA3wb+A3gT8K0kD9/5Emd8vt2BDwCnV9VZAFX1o6raWlW3AG8BHtB2vwr4fFVdV1U/Bz7O/F1fKUmSJElDazZfnfFa4DFV9YiqejjwWODf+lFMkgBvAy6rqtd2ta/s6vYk4OL28aeAkSR7tTe7eQRwaT9qkyRJkqTFZDbnve1eVd+cXKiqb7Wzf/3wEOB4YFOSyXOX/gZ4apLDaU5DHQOe3dby0ySvBb7ervt4VZ3dp9okSZIkadGYTVjckOStNF9RAbCG5iYzc66qvkhzt9WpPr6dbd7Fr2uTJEmSJM2B2YTFP6e5icxftstfoLl2UZIkSZI0pHYYFqvqZprrFl+7o76SJEmSpOEwmxvcSJIkSZIWGcOiJEmSJKmHYVGSJEmS1GOH1ywmuRfwYuDu3f2r6ug+1tU3Iwfuw4aTjx10GdI2Op0OY2tGB12GJEmS9CuzuRvq+4H/BN4CbO1vOZIkSZKkhWA2YXGiqv6/vlciSZIkSVowZgyLSe7UPvxokr8APgjcPLm+qn7S59okSZIkSQOyvZnF84EC0i6/uGtdAffoV1GSJEmSpMGaMSxW1UEASfasqv/tXpdkz34XJkmSJEkanNlcs/hl4H6zaNslbNp8A6tPOnvQZUjbWDcywQnzMC7HvBOwJEmSZml71yzeGTgQWJbkvvz6dNQ7AHvNQ22SJEmSpAHZ3sziY4ETgFXAa7vabwT+po81SZIkSZIGbHvXLJ4GnJbkD6rqA/NYkyRJkiRpwGZzzeLdk7xwStsNwPlVtbEPNUmSJEmSBux2s+hzBPAcmusXDwSeDTwOeEuSl/SxNkmSJEnSgMxmZnEVcL+q2gKQ5GXA2cDDab6L8VX9K0+SJEmSNAizmVn8DeDmruVfAgdU1f9Mad9pSU5Ncm2Si7vaXp3k8iQXJflgkhVt++5JTkuyKcllSV46l7VIkiRJ0mI2m7B4OnBekpe1s4pfAt6dZDlw6RzX8w6aU1y7nQPcp6oOBb4FTIbCPwKWVtUIcH/g2UlWz3E9kiRJkrQo7fA01Kr6xySfBB7cNj2nqja0j9fMZTFV9fmpga+qPt21+FXgDydXAcuTLAGWAb8AfjaX9UiSJEnSYjWbaxYBLgA2T/ZPcreq+kHfqprZnwLvbR+fCRwHXAPsBfxVVf1kADVJkiRJ0tDZYVhMciLwMuBHwFYgNLN6h/a3tJ46/h8wQXNaLMAD2nruAtwR+EKS/66qK6bZdi2wFmDFvvuzbmRifoqWZumAZczLuOx0On1/Dgng8PFxtm7d6pjT0NmyZYvjWkPHca2ZzGZm8fnAb1bV9f0uZiZJTgAeDzyqqqpt/mPgk1X1S+DaJF+i+ZqPnrBYVeuB9QBLVx5cp2ya7YSqND/WjUwwH+NybM1o359DAmDFCsbHxxkdHR10JdKc6nQ6jmsNHce1ZjKbG9xcCdzQ70JmkuRxwEuAJ1TVz7tW/QA4uu2zHDgKuHz+K5QkSZKk4TObqYwrgE6Ss+n6qoyqeu1cF5PkPcAosF+Sq2hOf30psBQ4JwnAV6vqOcB/AG9PcgnNqbFvr6qL5romSZIkSVqMZhMWf9D+7NH+9E1VPXWa5rfN0HcLzddnSJIkSZLm2Gy+OuMfAJLsNeU0UEmSJEnSkNrhNYtJHpTkUtrrAZMcluRNfa9MkiRJkjQws7nBzeuAxwLXA1TVhcDD+1mUJEmSJGmwZhMWqaorpzRt7UMtkiRJkqQFYjY3uLkyyYOBSrI7zfcuXtbfsiRJkiRJgzSbmcXnAM8FDgQ2A4cDf9HPoiRJkiRJgzWbu6FeB6zpbkvyApprGXc5Iwfuw4aTjx10GdI2Op0OY2tGB12GJEmS9CuzumZxGi+c0yokSZIkSQvKbQ2LmdMqJEmSJEkLym0NizWnVUiSJEmSFpQZr1lMciPTh8IAy/pWkSRJkiRp4GYMi1V1+/ksZL5s2nwDq086e9BlaEiMebMkSZIkDanbehqqJEmSJGmIGRYlSZIkST0Mi5IkSZKkHoZFSZIkSVIPw6IkSZIkqYdhUZIkSZLUw7AoSZIkSeqxy4TFJCuSnJnk8iSXJXlQ17p1SSrJfoOsUZIkSZKGxZJBF3ArvB74ZFX9YZI9gL0AktwVeAzwg0EWJ0mSJEnDZJeYWUyyD/Bw4G0AVfWLqhpvV/8b8BKgBlSeJEmSJA2dXWVm8SDgx8DbkxwGnA88H3g0sLmqLkwy48ZJ1gJrAVbsuz/rRib6X7EWhU6nMyf72bJly5ztS1oIDh8fZ+vWrY5rDR0/rzWMHNeaya4SFpcA9wNOrKrzkrweeDnNbONjdrRxVa0H1gMsXXlwnbJpVzlsLXRja0bnZD+dTofR0bnZl7QgrFjB+Pi441pDx89rDSPHtWayS5yGClwFXFVV57XLZ9KEx4OAC5OMAauAC5LceTAlSpIkSdLw2CXCYlX9ELgyyW+2TY8CLqiq36iq1VW1miZQ3q/tK0mSJEnaCbvS+ZgnAqe3d0K9AnjGgOuRJEmSpKG1y4TFqtoIHLGd9avnrxpJkiRJGm67xGmokiRJkqT5ZViUJEmSJPUwLEqSJEmSehgWJUmSJEk9DIuSJEmSpB67zN1Q58rIgfuw4eRjB12GJEmSJC1ozixKkiRJknoYFiVJkiRJPQyLkiRJkqQehkVJkiRJUg/DoiRJkiSpx6K7G+qmzTew+qSzB12GtI11IxOcsJ1xOeYdfCVJkjTPnFmUJEmSJPUwLEqSJEmSehgWJUmSJEk9DIuSJEmSpB6GRUmSJElSD8OiJEmSJKnHLhEWk+yZ5GtJLkxySZJ/aNtPT/LNJBcnOTXJ7oOuVZIkSZKGwS4RFoGbgaOr6jDgcOBxSY4CTgd+CxgBlgHPGlyJkiRJkjQ8lgy6gNmoqgK2tIu7tz9VVR+f7JPka8CqAZQnSZIkSUNnV5lZJMluSTYC1wLnVNV5Xet2B44HPjmo+iRJkiRpmOwSM4sAVbUVODzJCuCDSe5TVRe3q98EfL6qvjDdtknWAmsBVuy7P+tGJualZmm2DljGdsdlp9OZv2KkOXD4+Dhbt2517GrobNmyxXGtoeO41kx2mbA4qarGk3wOeBxwcZKXAfsDz97ONuuB9QBLVx5cp2za5Q5bQ27dyATbG5dja0bnrxhpLqxYwfj4OKOjo4OuRJpTnU7Hca2h47jWTHaJ01CT7N/OKJJkGfA7wOVJngU8FnhqVd0yyBolSZIkaZjsKlNsK4HTkuxGE3DfV1UfSzIBfB/4ShKAs6rqFQOsU5IkSZKGwi4RFqvqIuC+07TvEvVLkiRJ0q5mlzgNVZIkSZI0vwyLkiRJkqQehkVJkiRJUg/DoiRJkiSph2FRkiRJktTDsChJkiRJ6rHovnpi5MB92HDysYMuQ9pGp9NhbM3ooMuQJEmSfsWZRUn6/9u7/yDNrrJO4N9nGUnEZaclRCtOlCESsICRwUXLRVk7gdVoFFhFhQoIiE6BiiwVa41bVrFVlkuKNS4bSpcaNQ4KBEhMsUBE2RVfYVdKCAkw+SEsFQZNtIxIddaxCOw0Z//oNzrpMz3TPdPdt++dz6eqq7rve+65z337qTv9zbnvDQAAHWERAACAjrAIAABAR1gEAACgc9Y94Obwvfdn71U3D11G54iH7gAAADuIlUUAAAA6wiIAAAAdYREAAICOsAgAAEBHWAQAAKAjLAIAANARFgEAAOiMIixW1XVVdV9V3b5q+yur6s+r6o6qet1Q9QEAAEzNKMJikkNJLjt+Q1VdkuQ5SZ7SWntSkl8ZoC4AAIBJGkVYbK19IMnnV21+RZKrW2tfnI+5b9sLAwAAmKhdQxdwBh6f5BlV9ctJHkjyc621j5xoYFUdSHIgSRbOOz9X7ju2fVWu02w2G7oEBnT06FE9wKTsX1rK8vKyvmZyXK+ZIn3NWsYcFncleVSSb0/yrUneUVUXtdba6oGttYNJDibJORdc3K45vPNO+8gVi0OXwIBms1kWFxeHLgM2z8JClpaW9DWT43rNFOlr1jKK21DXcE+Sm9qKDyf5cpJHD1wTAADAJIw5LL4zySVJUlWPT/LwJJ8btCIAAICJ2Hn3Y55AVV2fZDHJo6vqniSvSXJdkuvm/zuNLyV58YluQQUAAGDjRhEWW2svWOOlF25rIQAAAGeJMd+GCgAAwBYRFgEAAOgIiwAAAHSERQAAADrCIgAAAJ1RPA11M+3bszu3XH350GUAAADsaFYWAQAA6AiLAAAAdIRFAAAAOsIiAAAAnbPuATeH770/e6+6eegy2AZHPMgIAABOm5VFAAAAOsIiAAAAHWERAACAjrAIAABAR1gEAACgIywCAADQERYBAADojD4sVtWrq+qOqrq9qq6vqnOHrgkAAGDsRh0Wq2pPkp9N8rTW2pOTPCzJ84etCgAAYPxGHRbndiX5yqraleQRSf5q4HoAAABGb9fQBZyJ1tq9VfUrSf4iyReSvK+19r7V46rqQJIDSbJw3vm5ct+x7S2UQcxms6FLWLejR4+Oql44lf1LS1leXtbXTI7rNVOkr1nLqMNiVX11kuckeWySpSQ3VNULW2tvPn5ca+1gkoNJcs4FF7drDo/6tFmnI1csDl3Cus1msywuLg5dBmyehYUsLS3paybH9Zop0tesZey3oT4ryWdaa3/bWvt/SW5K8vSBawIAABi9sYfFv0jy7VX1iKqqJM9MctfANQEAAIzeqMNia+3PktyY5NYkh7NyPgcHLQoAAGACRv/hvdbaa5K8Zug6AAAApmTUK4sAAABsDWERAACAjrAIAABAR1gEAACgIywCAADQGf3TUDdq357dueXqy4cuAwAAYEezsggAAEBHWAQAAKAjLAIAANARFgEAAOgIiwAAAHTOuqehHr73/uy96uahy2AHOuIpuQAA8I+sLAIAANARFgEAAOgIiwAAAHSERQAAADrCIgAAAB1hEQAAgI6wCAAAQGcSYbGqHlZVt1XVe4auBQAAYAomERaTvCrJXUMXAQAAMBWjD4tVdWGSy5P85tC1AAAATMWuoQvYBK9P8u+TPHKtAVV1IMmBJFk47/xcue/YNpXGmMxms8GOffTo0UGPD5tt/9JSlpeX9TWT43rNFOlr1jLqsFhV35/kvtbaR6tqca1xrbWDSQ4myTkXXNyuOTzq02aLHLlicbBjz2azLC4Od3zYdAsLWVpa0tdMjus1U6SvWcvYb0P9jiTPrqojSd6W5NKqevOwJQEAAIzfqMNia+0XWmsXttb2Jnl+kve31l44cFkAAACjN+qwCAAAwNaYzIf3WmuzJLOBywAAAJgEK4sAAAB0hEUAAAA6wiIAAAAdYREAAICOsAgAAEBnMk9DXa99e3bnlqsvH7oMAACAHc3KIgAAAB1hEQAAgI6wCAAAQEdYBAAAoHPWPeDm8L33Z+9VNw9dBms44uFDAACwI1hZBAAAoCMsAgAA0BEWAQAA6AiLAAAAdIRFAAAAOsIiAAAAHWERAACAzqjDYlV9fVX9cVXdWVV3VNWrhq4JAABgCnYNXcAZOpbkytbarVX1yCQfrar/0Vq7c+jCAAAAxmzUK4uttb9urd06//7vk9yVZM+wVQEAAIzf2FcW/1FV7U3y1CR/doLXDiQ5kCQL552fK/cd29baWL/ZbDZ0CYM4evToWXvuTNP+paUsLy/raybH9Zop0tesZRJhsar+eZLfS/LvWmv/d/XrrbWDSQ4myTkXXNyuOTyJ056kI1csDl3CIGazWRYXF4cuAzbPwkKWlpb0NZPjes0U6WvWMurbUJOkqr4iK0HxLa21m4auBwAAYApGHRarqpL8VpK7Wmu/OnQ9AAAAUzHqsJjkO5K8KMmlVfWx+df3DV0UAADA2I36w3uttf+VpIauAwAAYGrGvrIIAADAFhAWAQAA6AiLAAAAdIRFAAAAOsIiAAAAnVE/DfV07NuzO7dcffnQZQAAAOxoVhYBAADoCIsAAAB0hEUAAAA6wiIAAAAdYREAAIDOWfc01MP33p+9V908dBlMyBFP1wUAYIKsLAIAANARFgEAAOgIiwAAAHSERQAAADrCIgAAAB1hEQAAgI6wCAAAQGf0YbGqrquq+6rq9qFrAQAAmIrRh8Ukh5JcNnQRAAAAUzL6sNha+0CSzw9dBwAAwJTsGrqA7VBVB5IcSJKF887PlfuODVwRUzKbzc54jqNHj27KPLBT7F9ayvLysr5mclyvmSJ9zVrOirDYWjuY5GCSnHPBxe2aw2fFabNNjlyxeMZzzGazLC6e+TywYywsZGlpSV8zOa7XTJG+Zi2jvw0VAACAzScsAgAA0Bl9WKyq65N8KMkTquqeqnrZ0DUBAACM3eg/vNdae8HQNQAAAEzN6FcWAQAA2HzCIgAAAB1hEQAAgI6wCAAAQEdYBAAAoDP6p6Fu1L49u3PL1ZcPXQYAAMCOZmURAACAjrAIAABAR1gEAACgIywCAADQERYBAADoCIsAAAB0hEUAAAA6wiIAAAAdYREAAICOsAgAAEBHWAQAAKAjLAIAANARFgEAAOgIiwAAAHSERQAAADrCIgCbb//+HH3c44auAgA4A7uGLgCACXr96/Pp2SwXDl0HAHDarCwCAADQERYBAADoCIsAAAB0hEUAAAA6wiIAAAAdYREAAICOsAgAAEBHWAQAAKAjLAIAANARFgEAAOgIiwAAAHSERQAAADrCIgAAAB1hEQAAgI6wCAAAQEdYBAAAoCMsAgAA0BEWAQAA6AiLAAAAdIRFAAAAOsIiAAAAHWERAACAjrAIAABAp1prQ9ewrarqb5N8dug61rA7yf0TPf5mzX2682x0v42MX8/YU415dJLPrfN4Y6Ovt24OfT0cfb11c2xVX2/WOH09vuOP7W+Qjeyjr09NX5/YY1pr559yVGvN1w75SnJwqsffrLlPd56N7reR8esZe6oxSW4Z8ne/lV/6euvm0NfDfenrrZtjq/p6s8bp6/Edf2x/g2xkH329fb//nXb87Tovt6HuLO+e8PE3a+7TnWej+21k/HrGDv27HdLQ577T+/pM5tDXwxn63PX1xsdv9rgpGvrct+r4Y/sbZCP76OtTG/rcd3pfn9RZdxsq7ERVdUtr7WlD1wGbSV8zRfqaKdLXrMXKIuwMB4cuALaAvmaK9DVTpK85ISuLAAAAdKwsAgAA0BEWAQAA6AiLAAAAdIRFAAAAOsIi7HBVdVFV/VZV3Th0LXAmquqrqupNVfUbVXXF0PXAZnCNZoqq6rnza/Xbq+q7h66H4QiLsIWq6rqquq+qbl+1/bKq+mRVfbqqrjrZHK21u1trL9vaSuH0bLDHfzDJja21n0zy7G0vFtZpI33tGs1YbLCv3zm/Vr88yY8OUS87g7AIW+tQksuO31BVD0vya0m+N8kTk7ygqp5YVfuq6j2rvr5m+0uGDTmUdfZ4kguT/OV82PI21ggbdSjr72sYi0PZeF//4vx1zlK7hi4Apqy19oGq2rtq87cl+XRr7e4kqaq3JXlOa+21Sb5/eyuEM7ORHk9yT1YC48fiP1ayg22wr+/c3urg9Gykr6vqriRXJ3lva+3WbS2UHcU/1rD99uSfVleSlT+g96w1uKrOq6o3JnlqVf3CVhcHm2CtHr8pyQ9V1X9L8u4hCoMzcMK+do1m5Na6Xr8yybOSPK+qXj5EYewMVhZhh2ut/V1WPjMAo9Za+4ckLx26DthMrtFMUWvt2iTXDl0Hw7OyCNvv3iRff9zPF863wVTocaZIXzNF+pqTEhZh+30kycVV9diqeniS5yd518A1wWbS40yRvmaK9CKs+YsAAAUlSURBVDUnJSzCFqqq65N8KMkTquqeqnpZa+1Ykp9J8odJ7kryjtbaHUPWCadLjzNF+pop0tecjmqtDV0DAAAAO4yVRQAAADrCIgAAAB1hEQAAgI6wCAAAQEdYBAAAoCMsAgAA0BEWAdgRquprq+qtVXV3VX20qj5UVf/2FPt8XVXduEnHf0ZV3VFVH6uqr1z12s9W1V1V9ZYNzrlQVT+1GfVtlqparKr3nGLMQ+rezPcZgPEQFgEYXFVVkncm+UBr7aLW2r9M8vwkF55sv9baX7XWnrdJZVyR5LWttf2ttS+seu2nkvyb1toVG5xzYb7v2Dyk7k1+nwEYCWERgJ3g0iRfaq298cENrbXPttbekCRVtbeqPlhVt86/nn7c9tvn37+kqm6qqj+oqv9TVa870YGq6plVdVtVHa6q66rqnKr6iSQ/kuSXVq8eVtUbk1yU5L1V9eqq+rb5qudtVfWnVfWE+bgnVdWH5yuTn6iqi5NcneQb59v+8wlq+bH52I9X1e8ed07vn2//o6r6hvn2Q1V17fyYd1fV8+bb31ZVlx8356Gqel5VnVtVvz0/z9uq6pITHP8/VtXPHffz7VW1d3Xdq97nE8673vcfgPHYNXQBAJDkSUluPcnr92VlZe+BeQi7PsnTTjBuf5KnJvlikk9W1Rtaa3/54ItVdW6SQ0me2Vr7VFX9TpJXtNZeX1XfmeQ9rbWH3G7ZWnt5VV2W5JLW2ueq6l8keUZr7VhVPSvJf0ryQ0lenuS/ttbeUlUPT/KwJFcleXJrbf/qQqvqSUl+McnT5/M+av7SG5K8qbX2pqr68STXJnnu/LULknxnkm9K8q4kNyZ5e1aC7s3z4z4zySuS/PRK+W1fVX1TkvdV1eNP8h4f7yF1zwPkg04270nffwDGxcoiADtOVf3afLXtI/NNX5HkN6rqcJIbkjxxjV3/qLV2f2vtgSR3JnnMqtefkOQzrbVPzX9+U5J/vcHydie5Yb7S9l+yEnST5ENJ/kNV/XySx5zgVtbVLk1yQ2vtc0nSWvv8fPu/SvLW+fe/m5Vw+KB3tta+3Fq7M8nXzre9N8klVXVOku/Nyq28X5jv9+b53H+e5LNJ1hsWT+Zk857q/QdgRIRFAHaCO5J8y4M/tNZ+OisrZOfPN706yd8keUpWVhQfvsY8Xzzu++VszR00v5Tkj1trT07yA0nOndf81iTPTvKFJL9fVZduwbGPP7+aH/eBJLMk35PkR7Oy0rhex/LQvwXO3cT6tur9B2CbCIsA7ATvT3JuVb3iuG2POO773Un+urX25SQvysotnqfjk0n2VtXj5j+/KMmfbHCO3UnunX//kgc3VtVFSe5urV2b5L8n+eYkf5/kkWvM8/4kP1xV5833f/A21D/NysN9kpWH7nxwHTW9PclLkzwjyR/Mt31wvn/mt4l+Q1bO/3hHMg/pVfUtSR47336yutczLwATICwCMLjWWsvK5/K+q6o+U1Ufzsotoj8/H/LrSV5cVR/Pyuf1/uE0j/NAVkLVDfNbWr+c5I0n36vzuiSvrarb8tCVsx9JcntVfSzJk5P8Tmvt75L87/mDYx7ygJvW2h1JfjnJn8zP61fnL70yyUur6hNZCbOvWkdN70vyXUn+Z2vtS/Ntv57kn83P8+1JXtJa++Kq/X4vyaOq6o4kP5PkU/Pa1qx7nfMCMAG18u8zAAAA/BMriwAAAHSERQAAADrCIgAAAB1hEQAAgI6wCAAAQEdYBAAAoCMsAgAA0Pn/uN8L8ClQIgIAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1133ebb70>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import timeit\n",
    "\n",
    "n = np.arange(17)  # lengths = 2**n to evaluate\n",
    "reps = 20  # number of repetitions for timeit\n",
    "\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 numpy.fft import rfft, irfft; \\\n",
    "            x=np.random.randn(%d); h=np.random.randn(%d)' % (length, length)\n",
    "    # direct convolution\n",
    "    tc = timeit.timeit('np.convolve(x, x, mode=\"full\")', setup=tsetup, number=reps)\n",
    "    # fast convolution\n",
    "    tf = timeit.timeit('irfft(rfft(x, %d) * rfft(h, %d))' % (2*length, 2*length), setup=tsetup, number=reps)\n",
    "    # speedup by using the fast convolution\n",
    "    gain[N] = tc/tf\n",
    "\n",
    "# show the results\n",
    "plt.figure(figsize = (15, 10))\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 fast convolution')\n",
    "plt.ylabel('Length of signals')\n",
    "plt.title('Comparison between direct/fast convolution')\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise**\n",
    "\n",
    "* When is the fast convolution more efficient/faster than a direct convolution? \n",
    "* Why is it slower below a given signal length?\n",
    "* Is the trend of the gain as expected by the numerical complexity of the FFT?\n",
    "\n",
    "Solution: The gain in execution time of a fast convolution over a direct implementation of the the convolution for different signal lengths depends heavily on the particular implementation and hardware used. The fast convolution in this example is faster for two signals having a length equal or larger than 1024 samples. Discarding the outliers and short lengths, the overall trend in the gain is approximately logarithmic as predicted above."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbsphinx": "hidden"
   },
   "source": [
    "\n",
    "**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, Digital Signal Processing - Lecture notes featuring computational examples, 2016-2017*."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}