{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "nbsphinx": "hidden"
   },
   "source": [
    "# Realization of 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": [
    "## Quantization of Variables and Operations\n",
    "\n",
    "As for [non-recursive filters](../nonrecursive_filters/quantization_effects.ipynb#Quantization-Effects), the practical realization of recursive filters may suffer from the quantization of variables and algebraic operations. The effects of [coefficient quantization](quantization_of_coefficients.ipynb) were already discussed. This section takes a look at the quantization of variables.  We limit the investigations to the recursive part of a second-order section (SOS), since any recursive filter of order $N \\geq 2$ can be [decomposed into SOSs](cascaded_structures.ipynb).\n",
    "\n",
    "The computation of the output signal $y[k] = \\mathcal{H}\\{ x[k] \\}$ by a difference equation involves a number of multiplications and additions. As discussed already for [non-recursive filters](../nonrecursive_filters/quantization_effects.ipynb#Quantization-of-Signals-and-Operations), multiplying two numbers in a binary representation (e.g. [two's complement](https://en.wikipedia.org/wiki/Two's_complement) or [floating point](https://en.wikipedia.org/wiki/Floating_point)) requires requantization of the result to keep the word length constant. The addition of two numbers may fall outside the maximum/minimum values of the representation and may suffer from clipping.\n",
    "\n",
    "The resulting round-off and clipping errors depend on the number and sequence of algebraic operations. These depend on the structure used for implementation of the SOSs. For ease of illustration we limit our discussion to the [direct form I and II](direct_forms.ipynb). Similar insights can be achieved in a similar manner for other structures."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Analysis of Round-Off Errors\n",
    "\n",
    "Round-off errors are a consequence of reducing the word length after a multiplication. In order to investigate the influence of these errors on a recursive filter, the statistical model for [round-off errors in multipliers](../nonrecursive_filters/quantization_effects.ipynb#Model-for-round-off-errors-in-multipliers) as introduced for non-recursive filters is used. We furthermore neglect clipping. \n",
    "\n",
    "The difference equation for the recursive part of a SOS realized in direct form I or II is given as\n",
    "\n",
    "\\begin{equation}\n",
    "y[k] =  x[k] - a_1 \\, y[k-1] - a_2 \\, y[k-2]\n",
    "\\end{equation}\n",
    "\n",
    "where $a_0 = 1$, $a_1$ and $a_2$ denote the coefficients of the recursive part. Introducing the requantization after the multipliers into the difference equation yields the output signal $y_Q[k]$\n",
    "\n",
    "\\begin{equation}\n",
    "y_Q[k] = x[k] - \\mathcal{Q} \\{ a_{1} \\, y[k-1] \\} - \\mathcal{Q} \\{ a_{2} \\, y[k-2] \\}\n",
    "\\end{equation}\n",
    "\n",
    "where $\\mathcal{Q} \\{ \\cdot \\}$ denotes the requantizer. Requantization is a non-linear process which results in a requantization error. If the value to be requantized is much larger that the quantization step $Q$, the average statistical properties of this error can be modeled as additive uncorrelated white noise. Introducing the error into above difference equation gives\n",
    "\n",
    "\\begin{equation}\n",
    "y_Q[k] = x[k] - a_1 \\, y[k-1] - e_1[k] - a_2 \\, y[k-2] - e_2[k]\n",
    "\\end{equation}\n",
    "\n",
    "where the two white noise sources $e_1[k]$ and $e_2[k]$ are assumed to be uncorrelated to each other. This difference equation can be split into a set of two difference equations\n",
    "\n",
    "\\begin{align}\n",
    "y_Q[k] &= y[k] + e[k] \\\\\n",
    "y[k] &=  x[k] - a_1 \\, y[k-1] - a_2 \\, y[k-2] \\\\\n",
    "e[k] &= - e_1[k] - e_2[k] - a_1 \\, e[k-1] - a_2 \\, e[k-2]\n",
    "\\end{align}\n",
    "\n",
    "The first difference equation computes the desired output signal $y[k]$ as a result of the input signal $x[k]$. The second one the additive error $e[k]$ due to requantization as a result of the requantization error $- (e_1[k] + e_2[k])$ injected into the recursive filter.\n",
    "The power spectral density (PSD) $\\Phi_{ee}(\\mathrm{e}^{\\,\\mathrm{j}\\,\\Omega})$ of the error $e[k]$ is then given as\n",
    "\n",
    "\\begin{equation}\n",
    "\\Phi_{ee}(\\mathrm{e}^{\\,\\mathrm{j}\\,\\Omega}) = | H(\\mathrm{e}^{\\,\\mathrm{j}\\,\\Omega})|^2 \\cdot (\\Phi_{e_1 e_1}(\\mathrm{e}^{\\,\\mathrm{j}\\,\\Omega}) + \\Phi_{e_2 e_2}(\\mathrm{e}^{\\,\\mathrm{j}\\,\\Omega}))\n",
    "\\end{equation}\n",
    "\n",
    "According to the model for the requantization errors, their PSDs are given as $\\Phi_{e_1 e_1}(\\mathrm{e}^{\\,\\mathrm{j}\\,\\Omega}) = \\Phi_{e_2 e_2}(\\mathrm{e}^{\\,\\mathrm{j}\\,\\Omega}) = \\frac{Q^2}{12}$. Introducing this together with the transfer function of the SOS yields\n",
    "\n",
    "\\begin{equation}\n",
    "\\Phi_{ee}(\\mathrm{e}^{\\,\\mathrm{j}\\,\\Omega}) = \\left| \\frac{1}{1 + a_1 \\, \\mathrm{e}^{\\,-\\mathrm{j}\\,\\Omega} + a_2 \\, \\mathrm{e}^{\\,-\\mathrm{j}\\,2\\,\\Omega}} \\right|^2 \\cdot \\frac{Q^2}{6}\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Example - Round-off error of a SOS\n",
    "\n",
    "The following example evaluates the error $e[k] = y_Q[k] - y[k]$ for a SOS which only consists of a recursive part. The desired system response $y[k]$ is computed numerically by floating point operations with double precision, $y_Q[k]$ is computed by applying a uniform midtread quantizer after the multiplications. The system is excited by uniformly distributed white noise. Besides the PSD $\\Phi_{ee}(\\mathrm{e}^{\\,\\mathrm{j}\\,\\Omega})$, the signal-to-noise ratio (SNR) $10 \\cdot \\log_{10} \\left( \\frac{\\sigma_y^2}{\\sigma_e^2} \\right)$ in dB of the filter is evaluated."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "SNR due to requantization: 44.902056 dB\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmwAAAEYCAYAAADhzqpvAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzs3Xd4FWXax/Hvk04KCYGQhN4jxVAFUZqAjcUNi9iwIWIv6+vuWtZd21pwi2XX3nXtWEBdCxZio9fQW6jSQ0sgPc/7x0wwxCQkISdzkvw+13WukzMzZ+ae57Q7TxtjrUVERERE/FeA1wGIiIiISMWUsImIiIj4OSVsIiIiIn5OCZuIiIiIn1PCJiIiIuLnlLCJiIiI+DklbOL3jDGDjTGrvY6jLMaYYcaYrV7H4QvGmI3GmJFex1EWY8wDxpg9xpgdNbS/JGPMImNMpjHm5prYZyWP28YYk2WMCaytY7rHjTfGfO+e779q89i+5qsy9cfvIX+MSXxHCZv4jPuDn+1+eRbfnqzE86wxplPxY2vtD9baJB/F+Kox5gFf7NvdvzXGHHLP/WdjzKPFPyTGmEHGmJnGmAPGmL3GmJ+MMSe56yYYYwpLlNsGY8wrxpguvoq1Knxdbsc4dmvgD0A3a21CDe32NiDVWhtlrf13De3zV0onwdbazdbaSGttoa+OWY6rgT1AY2vtH2r52DXKV2Vam99D1eWPMYnvKGETXzvH/fIsvt3odUAe6GmtjQRGAOOBq4wxjYFPgf8AsUBL4D4gt8TzZrnPiwZGAtnAAmNMj9oMvjqMMUE+3H1bIMNau6uqTzSOsr732gLLjzuyuqMtsMJWcub02q4BFJEyWGt1080nN2AjMLKcdZ2A74ADOP/pv+su/x6wwCEgC7gAGAZsLbXfPwFp7nYvAfHA50Am8DXQpMT2U4Ad7rG+B7q7y68G8oE891ifuMtbAB8Au4ENwM0l9tUIeBXYB6xw49haQRlYoFOpWJ4E+gH7K3jeBODHMpZ/CrxfznOauOt3u/F9CrQqsT4V+Bvwk1tO04FmJdZfCmwCMoC7ynv9Kii3jcDt7uuSCwQBdwDr3eOtAH5X+hyBf7rxbgDOLrU+3X3uBuBifklci9xjv+puezIwE9gPLAGGlTrvB93zzi75erjrvwUKgRx3n13c50wq7/VwX9drgbVu7E8BpsT6q4CVJc67D/BfN+5s9zi3Ae3cfQWVeO99DOwF1gFXldjnvcB7wOvufpcD/Sp4D50CzMN5388DTnGXv1rq9SvrNX4VeAb4DOczNhIIdV+rzcBO4FmgUYnn/AnYDmwDJlLivV+J8nwC2AIcBBYAgytz3scqU2Cgu7z4lgNsdJ/bH5iF857ZjvO5DKnC91BX97z2uzH9tlT5PQX8z415DtCxnNepON7L3bLdA9xVYn0o8Lhbrtvcv0PddaVjuh342T3mamCEuzyAXz6LGW55xtbG74BuNXfzPADd6u+NihO2t3GSggAgDBhUYl3pJKf0l9JGYDZOktYS2AUsBHq7X27fAveU2H4iEFXii29xiXWvAg+UeByA84NxNxACdMBJGs50108GfsCpFWsNLKOSCRvQDSdxvBJo7H5xvgacTYkE0912AmUnbBOBneUcqylwLhDunu8UYGqJ9anuF3YXnMQzFZhcIrYsYIhbTo8CBRW8fkeVW4nXZbFbLo3cZefhJCEBOD96h4DEEueYj5PcBALX4fwgGSAC58c7yd02kV8S7dLvh5ZuWY5yj3O6+ziuxHlvBrrj/IgHl3E+qRydUJR+fNTr4b6unwIxQBucJPmsEuf8M3CSey6dgLZlfSb4dcL2HfA0zmeil7vf4h/de3ESjlFueT0MzC7n9YnFSSQvdc/5Ivdx0/JevzJe3wPAqfzyGX0cJ5mMxXl/fQI87G5/Fk4S18N97d6iagnbJTjv3yCc5u4dQFhlzvtYZVpiebAbR3HMfXES/SD3OSuBWyrzPeTuax3wZ5zvieE4SVJSifLbi5MUBgFvAu+UU9bF8b6A87nsifMPT1d3/f0433fNgTicf0z+VkZMSThJb4sS++3o/n2Lu49WOJ/v54C3q/qdrpu3NzWJiq9NNcbsL3G7yl2ej9Ms08Jam2Ot/bGK+/2PtXantfZnnARqjrV2kbU2F/gIJ3kDwFr7srU20113L9DTGBNdzn5Pwvmhv99am2etTcf5Ir3QXX8+8KC1dq+1dgtQmf5OC40x+3B+4F4EXrHWHgQG8csX9W5jzMfGmPhj7Gsbzg/mr1hrM6y1H1hrD1trM3FqlYaW2uwVa+0aa202zn/Zvdzl44BPrbXfu+X0V5yai6r6t7V2i7t/rLVTrLXbrLVF1tp3cWqk+pfYfpO19gXr9Dd6DScxKy6DIqCHMaaRtXa7tba8JstLgM+stZ+5x/kKmI/zA1/sVWvtcmttgbU2vxrnVZbJ1tr91trNwAx+KctJwN+ttfOsY521dtOxdub2zRsE3O5+JhbjvF8uLbHZj+55FuLULvUsZ3e/AdZaa//rnvPbwCrgnCqc3zRr7U/W2iKcBOIq4P/c934m8BBHfy5esdYus9YewvmcVZq19g33/Vtgrf0XTlJRsm9WZc+7Iv/G+YfhLveYC6y1s91jbsRJYkp/XspzMhCJ8x7Is9Z+i5PAX1Rimw+ttXOttQU4CVuvMvZT0n3W2mxr7RKcWuLic7wYuN9au8tauxun68SlZTy/EKfcuhljgq21G62169111+DU2m0t8T04zsddF6SGKWETXxtjrY0pcXvBXX4bTu3DXGPMcmPMxCrud2eJv7PLeBwJTt8bY8xkY8x6Y8xBnP/GAZqVs9+2QIuSSSbOf9HFSUQLnP9iix3zhxjoY61tYq3taK39i/sDiLV2pbV2grW2FU7NRAucWoyKtMT5z/1XjDHhxpjnjDGb3HP9Hogp1f+o5KjKw7jlVPq83B/djEqcW2klywZjzGXGmMUlyrIHR5f9kXistYfdPyPd41+A0+y43RjzP2PMCeUcsy1wXqnXbBBO8ldmXDWkvLJsjVOTWVUtgOJkqNgmnNe8vGOGlfOj24JfvzdL7+tYSpZZHE7N7YISZfyFu7z4eFX9XBxhjPmDMWalOwBnP06/zTLfJ1R83uXt/xqc2qjxxZ8/Y0wXY8ynxpgd7uflIcr/XiitBbCleF+uY71WkVSsos9myfLc5C47irV2HU5N2r3ALmPMO8aY4u3aAh+VeO1W4iR4x/oHUfyIEjbxhLV2h7X2KmttC5z//p4uOSKrBo0HUnD64ETjNBOAkyyCU8NV0hZgQ6kkM8paW1xbsx3nB7lYm5oI0lq7CqcZ5VgDCn6HU6NYlj/g1EoMsNY2xmnehF/OtSJHnZcxJhyniarckI+13BjTFqf28EacprgYnCbkysSDtfZLa+3pOInXKndfZdkC/LfUaxZhrZ1ciXjLcwgnQSlWldGoW4CO5ayrKI5tQKwxJqrEsjY4zatVtQ3nR7qkqu6rZKx7cP4R6l6ijKOtMygGjv25KLc8jTGDcfpenY/TNSAGpzm2Uu8TjvHauvv/G5BirT1QYtUzOO+rzu7n5c9VOOY2oHWpASzVfa0qc6ySr2Ubd9mvWGvfstYOcre3wCPuqi04/UNLfkbCrNNCIXWEEjbxhDHmPGNMK/fhPpwvl+Jh+Dtx+o7VhCic5pwMnB+Mh0qtL32sucBBY8ztxphGbg1dj+LpNnCaEe80xjRx47+pOkEZY05waxVauY9b4zSnzC5j20BjTHtjzH9wagnuq+Bcs4H9xphY4J4qhPQ+MNqdaiQEp99MRd8PlXmNInBe190AxpgrOHZCirttvDHmt8aYCJzXL4tf3h+lvQGcY4w50y2rMOPMj9eqnO0rYzEw1q217ITT77CyXgT+aIzp645K7eQmr1BBublN7DOBh91zSHaP+2Y14v8M6GKMGW+MCTLGXIDTT/HTauwLtybpBeAxY0xzAGNMS2PMme4m7wETjDHd3GS/9HuvovKMwukvuRsIMsbcjdPHs7LKLVP3c/UucJm1dk2p1VE4/SSz3Nrb6yq7X5xBBIeA24wxwcaYYTjNze9UIe7Kehv4izEmzhjTDKd/7RulNzLOXILDjTGhOH3+svnlM/Ms8GDx+9DdV4oPYhUfUsImvvaJOXoeto/c5ScBc4wxWTgdmX9vrd3grrsXeM2tvj//OI//Ok4Tws84o/VKJ0Qv4fT52G+Mmer2kTkHp7/JBpyahRdxaufASZY2ueum4/SnqY5MYABOGRxy41qGU0tWbKBbPgdxOks3Bk6y1i4tZ5+P43Ra3uPu74vKBuP2D7sBp7P4dpwkuqIJgY8qt3L2uQL4F85IvJ3AiTgjNSsjAKcstuE0AQ8Fri/nOFtwalH/jPOjvwVnxOLxfL89hjOKcidO37pKJ03W2ik4/Qffwnmdp/JLv8OHcX589xtj/ljG0y/CqQXehtMX8x63T16VWGszgNE4ZZiB0wVhtLV2T1X3VcLtOB3tZ7tNiF/j9jOz1n6O8/771t3m21LPrag8v8QZ4b0G57OVQ9WasCsq0xE4tXnvl/gOKu4L+UecGvhMnGT03VLPvZdyvoestXnAb3EGDO3BGShymVtTXtMewOmTmQYsxRlgVdYciKE4g6L24DSvNsf5TIAzCvdjYLoxJhPn+2GAD2IVHzLWVrWlQEREpGLGGIvT3LjO61hE6gPVsImIiIj4uVpN2IwxLxtjdhljlpVYFmuM+coYs9a9b+IuN8aYfxtj1hlj0owxfWozVhERERF/Uds1bK/iTLBY0h3AN9bazsA37mNw+gZ0dm9X44zoERGROsBaa9QcKlJzajVhs9Z+z6/nkErB6YSKez+mxPLX3YknZ+PMJ5WIiIiISAPjD7Mcx1trtwNYa7cXDxnHmYCw5Eihre6y7aV3YIy5GqcWjrCwsL5t2tTI1FhSjqKiIgIC1P3RV1S+vqcy9j2Vse+pjH2vNsp4zZo1e6y1ccfazh8StvKUNYFhmUNarbXPA88DJCUl2dWrV/syrgYvNTWVYcOGeR1GvaXy9T2Vse+pjH1PZex7tVHGxphKXRnEH1LzncVNne79Lnf5Vo6eObsV5czuLCIiIlKf+UPC9jFwufv35cC0Essvc0eLngwcKG46FREREWlIarVJ1BjzNs6ldZoZY7biXL5kMvCeMeZKYDNwnrv5Z8AonFmzDwNX1GasIiIiIv6iVhM2a+1F5awaUca2FudSOSIiIuKR/Px8tm7dSk5Ojteh1Lro6GhWrlxZI/sKCwujVatWBAcHV+v5/jzoQERERDy2detWoqKiaNeuHcaUNR6w/srMzCQqKuq492OtJSMjg61bt9K+fftq7cMf+rCJiIiIn8rJyaFp06YNLlmrScYYmjZtely1lErYREREpEJK1o7f8ZahEjYRERERP6eETURERMTPKWETERGROmnp0qUkJCSwbNkyr0PxOSVsIiIiUic99NBDzJw5k4ceesjrUHxO03qIiIhInfT2228D8NZbb3kcie+phk1ERESkHHfffTc5OTlcc801XHvttYwfP56ffvqp1uNQwiYiIiJ+b9iwYWzcuPGoZVdeeSVPP/10tfaXnZ3N0KFDKSwsLHeb3bt3k5eXxzPPPMO4ceN49tlnef3115k8eTIAeXl5DBkyhIKCggrjrAlK2ERERKROmTZtGh06dODjjz/mgQce4KSTTmL58uVV2sfLL7/M2LFjCQwMLHebBQsW0LdvXzZu3EjPnj0BCAoKIjo6mqKiIkJCQhgxYgTvvvvucZ1PZShhExERkTpj/fr1XHfddXzxxRdccMEFPPTQQ9x5552MGzeuwtqy0t58801SUlKOPN6wYQMpKSn069eP/v37s3r1ahYuXEi/fv0YMWIEU6ZMAWDPnj1YawkIcFKoMWPG8Oabb9bsSZZBCZuIiIjUGdOnT2fMmDF06dLlyLKxY8cSEBDA2rVrK7WPvLw80tPTadeuHeBc4H7SpEk8+uijzJ8/n3vvvZfJkyezadMm2rdvz5YtW7j//vuZMWMGI0aMYN++fRw4cACAHj16MG/evBo/z9I0SlREREQq5ZZbYPHimt1nr17w+ONVe05Zl3my1mKtrdTz9+zZQ0xMzJHHU6dOZfny5Zx77rkAFBQUsHz5cu666y4AbrjhBm644QYAlixZctS+AgMDCQkJqbELxZdHCZuIiIjUGSNHjuRvf/sbt95665Fl06ZNo6Cg4Khat4o0atToqAuxL1myhAcffJArr7zyqO0yMzMrtb/c3FzCwsIqtW11KWETERGRSqlqTZgvdO7cmSeffJIRI0aQnZ3NtGnTiI2N5YMPPiAwMJD8/HzuueceDh8+TF5eHk888cRRj59++mmaNGlCYWEhOTk5hIWFkZiYyJdffskVV1xBQEAAS5cupUePHpWKJyMjg7i4OIKDg3163urDJlIHzZ0LKSmQn+91JCIitW/s2LFs2LCBUaNGcdttt7FkyRJOPPFEAJ5//nmys7OJiYkhKyvrV4+LnXHGGfz4448ATJw4kaKiIrp27UqvXr145JFHymx2LcuMGTMYNWpUzZ9kKaphE6mD3noLPv4YNm+Gjh29jkZEpPYZY3jllVd+tXzRokU89dRThIaGAjBp0qSjHhe78cYbefTRRxk5ciSNGjXi/fffr1Ycb731Fg8//HC1nlsVSthE6qDiAUnbtythExEpKSUlhQkTJtC6dWuGDx/+q8dnnXUWAL179+a0006jsLCwwrnYKpKXl8eYMWNISkqqyVMokxI2kTqmoAAWLXL+3rHD21hERGrLhAkTjhrZWZ5zzjmHc84551fLyjJx4sTjiikkJITLLrvsqGWVjbOqlLCJ1DErVkB2tvP39u3exiIiUlsmTJjgdQiV4qs4NehApI6ZP/+Xv1XDJiLSMKiGTaSOmTcPGjeGiAjVsImINBRK2ETqmPnzoV8/OHBANWwiIg2FmkRF6pDcXFiyxEnYEhNVwyYi0lAoYROpQ9LSnMlyTzoJEhJUwyYi0lAoYROpQ4oHHJx0klPDtmsXFBZ6G5OIiPieEjaROmTePGjWDNq0cWrYiopg926voxIREV9TwiZSh8yf79SuGePUsIGaRUVEGgIlbCJ1xKFDsHy5k7CBU8MGGnggIlJdd999Nzk5OVxzzTVce+21jB8/np9++umYz5s6dSpXXXUVKSkpTJ8+vRYiVcImUmcsXuw0gfbr5zxWDZuINCTDhg1j48aNADz33HNcf/31R63v3r07q1atIjs7m6FDh1J4jA6+u3btIi8vj2eeeYZx48bx7LPP8vrrrzN58mTAuU7oWWedRUFBwa+OP2bMGF544QVeffVV3n333Zo90XIoYROpI4ov+F6csKmGTUQaqrS0NHr37n3kcU5ODps3b6Zz5868/PLLjB079pgXdJ83bx59+/Zl48aN9OzZE4CgoCCio6MpKioiJCSEoUOHVpiQPfDAA9xwww01c1LHoIRNpI6YPx9atvylZi0sDGJiVMMmIg3P0qVL6dOnz1GPu3TpQmBgIG+++SYpKSlH1m3YsIGUlBT69etH//79Wb16NeAkbP369WPEiBFMmTIFgD179mCtJSDASY9Gjx7Nm2+++avjW2u5/fbbOfvss4+Kw5d0pQOROmLevF/6rxVLSFANm4jUnlu+uIXFOxbX6D57JfTi8bMer9Jzli9fztixYzHGAJCVlcXo0aPJy8sjPT2ddu3aAZCfn8+kSZN4/vnn6dixI5999hmTJ0/mlVdeYcOGDbRv357PPvuM+++/n27dunHLLbfQsmVLDhw4QHR0NN26dWNecfNGCf/5z3/4+uuvOXDgAOvWrePaa6897nI4FiVsInXA/v2wZg1cfvnRyzV5rog0NFu2bCEuLo5Vq1YdWXbjjTfSoUMH9uzZQ0xMzJHlU6dOZfny5Zx77rkAFBQUMHjwYJ566ik6deoEwA033HCkWXPJkiVHHSswMJCQkBAyMzOPWn7zzTdz8803++T8yuM3CZsx5v+ASYAFlgJXAInAO0AssBC41Fqb51mQIh5ZuNC5L+6/ViwxEebOrf14RKRhqmpNmC+kpaXRvXv3o5atWLGClJQUGjVqRE5OzpHlS5Ys4cEHH+TKK6+s9vFyc3MJCwur9vNril/0YTPGtARuBvpZa3sAgcCFwCPAY9bazsA+oPolLlKHlR5wUEw1bCLS0CxdupRu3bodtWz58uUkJyfTpEkTCgsLjyRtiYmJfPnllxQVFR15rrW20sfKyMggLi6O4ODgmjuBavKLhM0VBDQyxgQB4cB2YDjwvrv+NWCMR7GJeGr+fOjYEWJjj16emOjMz1aqtl5EpN4qnbDt3bsXay3x8fEAnHHGGfz4448ATJw4kaKiIrp27UqvXr145JFHjvR7q4wffviBUaNG1ewJVJNfNIlaa382xvwT2AxkA9OBBcB+a22Bu9lWoGVZzzfGXA1cDRAXF0dqaqrPY27IsrKyVMY+VFb5zp17Em3bHiY1dflRy/ftiwe6Mm3aHFq1yq69IOs4vYd9T2Xse7VVxtHR0b/qw+WFwsJCsrKyePbZZwGOxBQcHMz69euPPL7iiit48sknGTBgAACvvPLKUfupyrm899573HfffWRmZh45/vGURU5OTrVfM79I2IwxTYAUoD2wH5gCnF3GpmXWY1prnweeB0hKSrLDhg3zTaACQGpqKipj3yldvkVFsHMnXHBBxK/KvaAAHn4Y2rQZwJAhtRtnXab3sO+pjH2vtsp45cqVREVF+fw4xxIYGEhkZOQxYxk0aBBr1qwhPDz8mHOxVSQvL4/Ro0cfmbajssevSFhY2FHzx1WFXyRswEhgg7V2N4Ax5kPgFCDGGBPk1rK1ArZ5GKOIJ7Zvh9xc6NDh1+uKJ89VPzYRqe8mTJhw1AjQikycOPG4jxcSEsL48eOrdXxf8Jc+bJuBk40x4cZpXB4BrABmAOPcbS4HpnkUn4hn0tOd+7IStuJJdDUXm4jUd14nTF4f3y8SNmvtHJzBBQtxpvQIwGnivB241RizDmgKvORZkCIeqShhi42F4GDVsImI1Hf+0iSKtfYe4J5Si9OB/h6EI+I30tMhIADatv31OmN0tQMRkYbAL2rYRKR86enQpo1Tk1YWzcUmIr5WlbnLpGzHW4ZK2ET83Pr1ZTeHFktMVA2biPhOWFgYGRkZStqOg7WWjIyM47pigt80iYpI2dLT4Zxzyl+fkACzZ9dePCLSsLRq1YqtW7eye/dur0OpdTk5OTV2WaqwsDBatWpV7ecrYRPxY4cOOXOwHauGbfduZ062IH2iRaSGBQcH0759e6/D8ERqamq1502raWoSFfFjGzY49xUlbAkJYC3s2lU7MYmISO1Twibixyqa0qOY5mITEan/lLCJ+LHihK1jR+d++vrpPDrrUTYf2HxkG13tQESk/lPCJuLH1q+H6Gho0gTu+PoOznzjTP4w/Q90e6obK3avAFTDJiLSEChhE/Fj6elOc+jMLT/xyE+PMLHXRJZfv5zw4HDOn3I++YX5xMc726qGTUSk/lLCJuLH0tOhfQfLLV/eQpvoNjxx9hN0i+vGc6OfY/nu5UxZMYXQUOcSVaphExGpv5SwifipoiJnlGh4p/nM3zaf20+9nciQSABSTkghqWkSj81+DGutrnYgIlLPKWET8VPbt0NuLqQ3eY6I4AguSb7kyLoAE8DNA25m/rb5LNy+UFc7EBGp55Swifip9HQgMJfF+e9xfvfzaRza+Kj153c/n0ATyIcrP1QNm4hIPaeETcRPpacD7WdwuDCTc7ue+6v1zcKbMaTtED5a9dGRGjZd6k9EpH5Swibip9avB06YRkRwBCM6jChzm7Fdx7Jyz0oCmq8iJwcOHqzdGEVEpHYoYRPxU+vTLYFdP+HMTmcSFlT2xYdHdxkNwI7wrwD1YxMRqa+UsIn4qeU7V1EY8TNndTyr3G3axbSjfUx71ttvAfVjExGpr5Swifip9KIZAAxvP7zC7U5rdxpLD6aCKVQNm4hIPaWETcQPHToEmc2+JYY2dGhSwZXfcRK6g/n7IWEJmzbVUoAiIlKrlLCJ+KH16UXQLpXkxqdhjKlw22HthgEQ2/t7FiyoheBERKTWKWET8UPfr1gN4RkMaTvkmNu2bNyS1o1bE5E0m7lzayE4ERGpdUrYRPzQTxvnAHB28smV2n5g64Fkxcxi82bYudOXkYmIiBeUsIn4oWX75kBuYwZ0OKFS25/c8mT22c0QtY1583wcnIiI1DolbCJ+aHPRHKIO9CcwoHIf0YGtBwIQ0EbNoiIi9ZESNhE/k1OYw8FGabRkQKWf0zuhN0EBQTRLnq+ETUSkHlLCJuJnVh1cAwGFdG1c+YQtNCiU7nHdCWm7kLlzdU1REZH6RgmbiJ+Zv30tACe3rnzCBtA7sTcHwhexb591rkMqIiL1hhI2ET8zd9ta2Neeswc3r9Lzeif0JrNoF0RtV7OoiEg9o4RNxM9sLFhG5IEB9OhRtef1SewDcKRZVERE6g8lbCJ+ZPnm7eQ3+pkBrQZwjAsc/ErP+J4YDAm9FylhExGpZ5SwifiR5z5xJlEbP6R/lZ8bFRpFp9hOhLRZxMKFkJ9f09GJiIhXlLCJ+JEvFi8CazhvSHK1nl888CA3F5YureHgRETEM0rYRPzE/v2w7tBiovI6EBUaWa199Enow+78jdBor654ICJSjyhhE/ETn34KNn4RXaI6VXsfvRN7A9C4y2L1YxMRqUf8JmEzxsQYY943xqwyxqw0xgw0xsQaY74yxqx175t4HaeIr7w9dS/EbKJ3i7bV3kfvBCdha9FXAw9EROoTv0nYgCeAL6y1JwA9gZXAHcA31trOwDfuY5F6JysLvlq2GOC4atjiIuJoGdWS4DaLWL4cMjNrKkIREfGSXyRsxpjGwBDgJQBrbZ61dj+QArzmbvYaMMabCEV867PPID92EQCdIzsf1756JvQks1Ea1sLChTURnYiIeC3I6wBcHYDdwCvGmJ7AAuD3QLy1djuAtXa7MabMqd+NMVcDVwPExcWRmppaK0E3VFlZWSrjGvb0090IabuAxiHNCMoLOq7yjcmJYUv2lxCYx9tvb8FiCyOFAAAgAElEQVTaLTUXaD2h97DvqYx9T2Xse/5Uxv6SsAUBfYCbrLVzjDFPUIXmT2vt88DzAElJSXbYsGE+CVIcqampqIxrTnY2zJ8PEbekcXK7k4mMjDyu8t3RbAdvbXmLlr1WkpHRk2HDOtZcsPWE3sO+pzL2PZWx7/lTGftFkyiwFdhqrZ3jPn4fJ4HbaYxJBHDvd3kUn4jPfPklHMrN5kDwKnrF9zru/SXHO3O4JfZKIy3tuHcnIiJ+wC8SNmvtDmCLMSbJXTQCWAF8DFzuLrscmOZBeCI+NX06hLdfShGFR6blOB5dmnYhNDCUgMQ0NmyAwsIaCFJERDzlFwmb6ybgTWNMGtALeAiYDJxujFkLnO4+FqlXZs2CNgOcAQfF03Icj6CAILo3705WxBLy82GLurCJiNR5/tKHDWvtYqBfGatG1HYsIrUlKwvS0qDXbxcRHRpNu5h2bGLTce83OT6ZaXs+B2D9emjX7rh3KSIiHvKnGjaRBmf+fCgqgqyoRfRK6IUxpkb2m9w8mX35OyFiJ+vW1cguRUTEQ0rYRDw0axYQUMDmnLQaaQ4t1jOhJwDBrdKUsImI1ANK2EQ8NGsWtO+3hpzCnBoZcFDsxOYnAtCkaxrr19fYbkVExCOVTtiMMScYY0YYYyJLLT+r5sMSqf+sdRK21v1rbsBBsbiIOBIjEwlprRo2EZH6oFIJmzHmZpwpNW4ClhljUkqsfsgXgYnUd+vXw549ENp2ESGBIZzQ7IQa3X/PhJ7kxjg1bNbW6K5FRKSWVbaG7Sqgr7V2DDAM+Ksx5vfuuprpJS3SwMya5dwfCF9EcnwywYHBNbr/5ObJ7AtcweGcfHbsqNFdi4hILatswhZorc0CsNZuxEnazjbGPIoSNpFqmT0bIqMsazMX1WhzaLHk+GQKyINmq9UsKiJSx1U2YdthjDlyzRw3eRsNNANO9EVgIvXdrFnQc8hm9uXs80nCVjxSlHgNPBARqesqm7BdBhzVqGKtLbDWXgYMqfGoROq5Q4ecCXNb9HEHHNTgCNFiSU2TCA4IxiQuUQ2biEgdV6mEzVq71b3eZ1nrfqrZkETqv3nznGt8BrdeRIAJOHLB9poUHBhMt7huhLVVDZuISF133POwGWOuqIlARBqS4gEHe0MXkdQ0ifDgcJ8cp2dCT4qaaWoPEZG6riYmzr2vBvYh0qDMmgVJSbB0z0KfNIcWS26eTG7oNtb+vMdnxxAREd+r1MXfjTFp5a0C4msuHJH6z1pnhOiIc3bzTubP9Eno47NjFTe1HghLY+/e4cTG+uxQIiLiQ5VK2HCSsjOBfaWWG2BmjUYkUs+lp8Pu3RDfaxHs9c2Ag2IlR4quWzec/v19digREfGhyjaJfgpEWms3lbptBFJ9Fp1IPVTcf80kOiNEeyX0qmDr49M8ojlNQ+M1tYeISB1X2YTtRaDM0aDW2vE1F45I/TdrFkRFwTa7iLbRbYlt5Nt2yl6JyRCvqT1EROqyyiZslwMLjDHvGGMmGGMSfBmUSH02ezb07w+Ldy7yaXNosd6JPaH5ctauL/D5sURExDcqOw/btdbaPsC9QBPgVWPMLGPMQ8aYIcaYQF8GKVJfWAtr1kCXHlmszVjrkysclJYcnwxBuSzbvtbnxxIREd+o0rQe1tpV1trHrLVnAcOBH4HzgDm+CE6kvtm3D7KyIKjlEiy29hI2IP3QEp8fS0REfKPa87BZa7OttZ9Za2+y1varyaBE6quNG537w9G+uyRVaV3juhJAEAfC0sjK8vnhRETEB46ZsBljTjfGvFB88XdjzNW+D0ukftq0ybnfHbSIuPA4Wka19PkxQwJDaBXaVSNFRUTqsMrUsF0P/Am4xBgzHPDdHAQi9Vxxwpae7VzhwBhTK8ft3swZKaqETUSkbqpMwrbbWrvfWvtH4AzgJB/HJFJvbdoEjRpnszJjKSe1qL2P0sD2PSF6K2lr99baMUVEpOZUJmH7X/Ef1to7gNd9F45I/bZpEzRPXkShLazVhK1/W2fgwYKtS2vtmCIiUnOOmbBZa6cBGGNaGGMuBfKMMUN9HplIPbRpE4R3mgfASS1rL2ErHim6en95lwUWERF/VqlRosaYM4AFwChgIPC4MWa1MeYUXwYnUt9s2gQF8fNoEdWCFlEtau24CZEJhBbEsa1QU3uIiNRFlb34+wPAYGvtkYvbGGMGAi8YY64EDllrl/kiQJH6IisLMjLAhM/j1FpsDgUwxtAiMJkNkWnk5kJoaK0eXkREjlNl52ELKZmsAVhrZwFjgf8Cb9d0YCL1zaZNQNh+9tg1tdp/rVhSTDI0X8bK1YW1fmwRETk+lU3YcowxcaUXWmvXAIXAiBqNSqQe2rQJaDEfgP4t+9f68Yd37wnB2XyUqqvAi4jUNZVN2P4BTDXGHNXpxhjTDMi11u6q8chE6hknYXMGHPRrUfsXBxnR3Rl48FWaBh6IiNQ1lerDZq39wBgTCswyxiwAlgAhwPk4/dtE5Bg2bQLTei4dYzvRpFGTWj9+t+ZdMTaQtJ1pWHsetTRnr4iI1IBKX0vUWvsW0BX4FIgG8oHx1trXfBSbSL2yaRMEtJrnSf81gLCgMBKDT+BQ1BLWrPEkBBERqabKjhIFwFp7GHjZR7GI1Gtrdm6h8ISfGdBygGcx9GmZzLY9M/nuO0hK8iwMERGpokrXsInI8dlQ8BMAp7Y51bMYTumYDDGbmP79Ac9iEBGRqlPCJlILcnNhX+RMggmnZ3xPz+LoleAcO3XlEqz1LAwREakiv0rYjDGBxphFxphP3cftjTFzjDFrjTHvGmNCvI5RpDq2bAFa/0TH0P4EBwZ7FkffFn0ByAidx/r1noUhIiJV5FcJG/B7YGWJx48Aj1lrOwP7gCs9iUrkOK1Kz4KEJfRt7l1zKEDziOa0jGgHLefy3XeehiIiIlXgNwmbMaYV8BvgRfexAYYD77ubvAaM8SY6kePz3bq5EFDI8M7eJmwAp7btT0AbJWwiInVJlUaJ+tjjwG1AlPu4KbDfWlvgPt4KtCzricaYq4GrAeLi4khNTfVtpA1cVlaWyriKvl79P4iFmKzCY5adr8u3aU5Tihpv5PP3tzBjxvoGOR+b3sO+pzL2PZWx7/lTGftFwmaMGQ3sstYuMMYMK15cxqZldpO21j4PPA+QlJRkhw0bVtZmUkNSU1NRGVfN7s8fIWhvd8aePfqY2/q6fAM3BfJM+jPsCV1Ku3ajaN/eZ4fyW3oP+57K2PdUxr7nT2XsL02ipwK/NcZsBN7BaQp9HIgxxhQnla2Abd6EJ1J9RbaIHSGzaHr4FK9DAaBPYh8CCFA/NhGROsQvEjZr7Z3W2lbW2nbAhcC31tqLgRnAOHezy4FpHoUoUm3Ldi2jMOgA7YO8778GEBESQY/mPQhup4RNRKSu8IuErQK3A7caY9bh9Gl7yeN4RKrsm/QZAPRpMszbQEro37I/ptVcUr/TZGwiInWB3yVs1tpUa+1o9+90a21/a20na+151tpcr+MTqarPV30LezuQ3Lat16Ec0b9lf/KCMti4fwObN3sdjYiIHIvfJWwi9UlhUSGztn0HG4bjR/ka/Vv2d/5oOUfNoiIidYASNhEfWrxjMVkFB2DjaX6VsHVv3p2I4AhCOs5k5kyvoxERkWNRwibiQ99u+Nb5Y8NptGnjbSwlBQUEMbD1QEI6/8CCBV5HIyIix6KETcSHZmycQUx+V+IaJRIR4XU0RxvcZjBZEWksWb2fvDyvoxERkYooYRPxkfzCfL7f9D1RGf7VHFpscJvBYCx5zWexfLnX0YiISEWUsIn4yPxt8zmUf4jCtf414KDYgFYDCDJB0FbNoiIi/k4Jm4iPfLHuCwJMAHsX+mcNW3hwOH1b9CWwww/Mn+91NCIiUhElbCI+8tm6z+gXP5CcfbF+mbCB0yxalDiXeYtyvA5FREQqoIRNxAd2Zu1k/rb5dA8ZBUC7dt7GU55BbQZhA/JI2zNfAw9ERPyYEjYRH/hi3RcALHpvFAkJMHKkxwGV49Q2zvVNC1r8wLJlHgcjIiLlUsIm4gOfrfuM2OBEFn/Rk7/+FcLDvY6obM3Cm9E5uju0n6F+bCIifkwJm0gNKygq4Mt1X2LWjaJ9e8OkSV5HVLFRJ5wObX9gzgL1YxMR8VdK2ERq2KwtsziQe4CM2aO47z4ICfE6ooqd0fF0CMrh+40/eh2KiIiUQwmbSA2btuoTKAzmhJCRjB/vdTTHNrTtUAJsMOl8RW6u19GIiEhZlLCJ1CBrLa/Nex/SR/LwvY0JDPQ6omOLCInghIhTKGr/lQYeiIj4KSVsIjVo9qZF7CncQPvs80hJ8Tqayju7y+mQuIgZc3d7HYqIiJRBCZtIDfrL21OgMIh/XZWCMV5HU3nn9T0dgM9XfeNxJCIiUhYlbCI1ZOtWy4yd79M0czi/OyvW63CqpF+LvgTlN2Fx5ldehyIiImVQwiZSQybckYZtso7/O3Oc16FUWWBAIJ0CRrI39nMOZxd5HY6IiJSihE2kBnz6KXyzfQoBBHL14DFeh1Mtp7f5LURt570fNYOuiIi/UcImcpwOHYLrbygiuM9bnNZuOHERcV6HVC0TB4+CokDeXTLN61BERKQUJWwix+m++2CL+YH8yA1M7DPB63CqrWeXWIJ+HsLs/UrYRET8jRI2keOQlgaPPgqdz3uVxqGNGXNC3WwOBTAGOhaksD94Oev2rvM6HBERKUEJm0g1FRXBNddATPMstjWZwgXdLyA82E+v8l5JZ7d3Jo97a4Fq2URE/IkSNpFqeuEFmD0bxt71AYfyDzGh1wSvQzpuF57VDnb05O1FSthERPyJEjaRatixA26/HU47DVaHv0zn2M4MbDXQ67COW79+0Gjj71h1+Ee2ZW7zOhwREXEpYROphltvhexsuOXhpXy/+Xsm9ZmEqUuXNihHYCAMa3YRGMvbS9/xOhwREXEpYROpounT4e234c474dNd/6FRUCMm9ZnkdVg1ZtxpXeDnfrw0702vQxEREZcSNpEqyM6G66+Hzp3h6t/v5Y20N7j4xIuJbVS3LkVVkTPOAJZezMr9C1m1Z5XX4YiICErYRKrkwQdh/Xp49ll4c8VLZBdkc9OAm7wOq0a1agVd8i4EG8CbaaplExHxB0rYRCpp3Tr4+9/hkktgyLACnpr3FMPaDSM5Ptnr0Grc6GEJmA0jeTPtLay1XocjItLgKWETqaQ33oCCAnjkEXhr6VtsOrCJWwbc4nVYPnHmmWAXX8KGA+l8v+l7r8MREWnwlLCJVNKHH8KgQRCfUMiDPzxIcnwyv036rddh+cTgwRCafi6hRTE8M/8Zr8MREWnwlLCJVML69bB0KfzudzBlxRTWZKzhL4P/Ui+m8ihLo0Yw9JRwwtdO4MOVH7Iza6fXIYmINGh+kbAZY1obY2YYY1YaY5YbY37vLo81xnxljFnr3jfxOlZpmD76yLlPGVPEA98/QLe4bpzb7Vxvg/KxM8+EfdOvJb8on5cWveR1OCIiDZpfJGxAAfAHa21X4GTgBmNMN+AO4BtrbWfgG/exSK376CPo1QvmZL3L8t3LuWvwXQQYf/n4+MaZZwIZSZwQMpznFjxHYVGh1yGJiDRYfvGLY63dbq1d6P6dCawEWgIpwGvuZq8BY7yJUBqyHTtg1iwYPSaHO7+5k14Jvbiwx4Veh+Vz3bpBy5YQvfY6Nh/YzCdrPvE6JBGRBsv425B9Y0w74HugB7DZWhtTYt0+a+2vmkWNMVcDVwPExcX1fe+992on2AYqKyuLyMhIr8OoNR9/nMhjjyVx3uM3MWX/k/wr+V/0adLHZ8fzp/L9+9+T+O6HJjS+qz2xIbE82fvJetFvz5/KuL5SGfueytj3aqOMTzvttAXW2n7H2i7Ip1FUkTEmEvgAuMVae7CyPwzW2ueB5wGSkpLssGHDfBajQGpqKg2pjCdPhvbdd/Nl9uuM7jKaW393q0+P50/lu3s3fP45jIu4m9f3Xo9pbxjWbpjXYR03fyrj+kpl7HsqY9/zpzL2iyZRAGNMME6y9qa19kN38U5jTKK7PhHY5VV80jDt3w/ffAORY+7gUN4hHhn5iNch1aqUFBgyBN654wpiQ+J5+MeHvQ5JRKRB8ouEzThVaS8BK621j5ZY9TFwufv35cC02o5NGrb//Q8KWs1gafDL/OmUP9EtrpvXIdWqkBD44ANoGR9G3vf/x/T101mwbYHXYYmINDh+kbABpwKXAsONMYvd2yhgMnC6MWYtcLr7WKTWvD81h8Ax19CxSUfuHnq31+F4olkz+OQTYN51BObFcNfXDbMcRES85Bd92Ky1PwLldVgbUZuxiBTLzoZPs+6lMGYtz47+ikbBjbwOyTPdu8M7rzVm9MN/5suQ2/gmfQYjOpzmdVgiIg2Gv9Swifidf304g4IBf+fs5pMY2WGk1+F47je/gYfH3AQHWnP+87exeUuR1yGJiDQYSthEAGvhySfhootg6FBo3y2Dvy68lIB9nXnrsse9Ds9v3P6HMMZEP8DeRvPplPIed94JBw54HZWISP2nhE0EZyToTTfBzJlgKSJ/1EQConbxxNC3iYmI8Do8v2EMvP/Xi+napCfBo/7E5Ecz6dgRXnjB68hEROo3JWzS4FkL993nzOq/Zg0Mv+9+fo76mEfP+gc3jvXdBLl1VWBAIC/97hkOB29l/At3c+KJcPXV8Nprx36uiIhUjxI2afBmzIAff4Q77oD/pX/Ifd/dx+U9L+fmATd7HZrfGth6INf2vZZ3Nv6bya8tYPhwJ2mbOdPryERE6iclbNKgFdeutWgByaPmcOlHl9K/ZX+eHf1svbgEky89PPJh4sLjuPbzSbz5Th6tW8PYsbBli9eRiYjUP0rYpEFLTYXvv4cJf1rJ794fRUJkAtMunEZYUJjXofm9mLAYnh39LIt3LObRxX/h44/h8GHn6giHDnkdnYhI/aKETRq0++6DuKT1vFZ0BsEBwUy/ZDoJkQleh1VnjDlhDNf0vYZ/zPwH28O+4Z13YPFiuOIKKNKsHyIiNUYJmzRY330H3y1bQ974oWQXHOaLS76gY2xHr8Oqcx4981FOaHYCl350KX2G7OCRR2DKFOjXDz77zGl2FhGR46OETRqsP/0zjYCJQwlplMeMy2fQK6GX1yHVSeHB4bw77l0O5B5g7LtjuemWXP77X9i/35lsd9Agp+lZRESqTwmbNEj//Gg6804cRGRkAKkTUkmOT/Y6pDotOT6Z18a8xqyts7j+s+u4+GLLqlXwzDOwcSOcdhr86U9eRykiUncpYZMGxVrLU3Of4bbFvyEoqz3zr5lDt7huXodVL4zrNo67h9zNK4tf4W/f/42QELj2Wli3Di67DB57DFav9jpKEZG6SQmbNBiH8w8zYdoEbvz8euz603mq7w90jm/ldVj1yj3D7uHynpdzT+o9PDn3SQAaNYJ//hPCw+GuuzwOUESkjlLCJg3C/G3z6fd8P/675L9ELbiXvqs/ZdKljb0Oq94JMAG8+NsXSUlK4abPb+KlhS8BEBcHf/wjfPABzJnjcZAiInWQEjap1/IL87k39V5OfvFkDuQe4GK+IPOTe3js0QAC9O73iaCAIN4Z9w5ndTqLSZ9M4onZTwBw663QvDncfrtGjoqIVJV+sqTemrllJgNeHMB9393HhT0uZPqYZXz49zMYNw4GD/Y6uvotLCiMqRdM5dyu53LLl7fw12//SnhEEX/9qzOdyhdfeB2hiEjdooRN6p1tmdu49KNLOfXlU9l1aBfvn/c+b4x9g3/c34SCAnjkEa8jbBhCg0J5Z9w7XNn7Sh744QEu+uAiLrniMB06wJ13amJdEZGqUMImdVp6Olx+OXTrBt/MyuDP3/yZLv/pwnvL3+OuwXex6sZVnNvtXBYsgNdeg1tugQ4dvI664QgKCOKFc17gH6f/gynLp3DaG6dy3d2rWLIE3n7b6+hEROqOIK8DEKmOLVvggQfg5ZchMDKDoMGPM/KTJzChWZzf/XweHP7gkasWzJgBV17pdHz/8589DrwBMsbwx1P+SLe4blz20WXck9GXNin/5o47JzJwoFECLSJSCaphkzpl/Xq48Ubo1Ale/mQFSX+4BvOHVhzq+wCx+87EPp1G7/Xv0KFJRw4cgGuugeHDISAApk6F6Givz6DhGtV5FGnXpTGw1UA2957EzsHn03vwDj791OvIRET8nxI28XvWwk8/wdix0KnrYZ6d+QbNbh1BwTXdWR/1Opf2vJRl1y1j66NTuHB4D+64A8aNg+7d4cUXnekk0tLglFO8PhNpEdWC6ZdOZ/KIyZgTPibriiTO+duT3HlXIYWFXkcnIuK/1CQqfufnn2HFCli1ClauhNlzClmU8RMhJ71OyJ/fI89kEhrTngd6P8DVfa8mLiLuyHPfegu6dIH774cePeDDD6F/fw9PRn4lwARw+6DbGdt1LNd9egPfBN/E5O0v88X5/+DBK0fQty/Ex3sdpYiIf1HCJn4jNxeuuMLtjB6YC+2/JbjnR3DmxxC6k+DgCM7vfj4Tek1gUJtBBJhfVxAbA/fdB+PHQ/v2EBJS++chldO5aWe+uuxL3l/xPtd8dCuLE0fym/dOg2sfoEXhKfTtC5dcAuef73WkIiLeU8ImfuHgQctZF69l1s6v6fzXr/g55BsOF2YSGhLJqM6j+N0Jv2N0l9FEhkRWan9JST4OWGqEMYbzup/HOUnn8O+fnueRkIfY2/5Uig4PZ+6s/+OTC0exfn0Ad9zhJOMiIg2VEjbxRGFRISt2r2DmlpnMWDeLjxZ/S16/LQDkRbfl4o4XkpKUwogOIwgLCvM4WvG1sKAwbht6MzcMvJJn5j/D47Mf5+cR5xA1uAt/nnojm3ZdzNP/itXVKUSkwVLCJj6xdy/84x/www9gseSHb+Rw4yUcjl5EXvNZZITNIYeDAATmxGE3D+a6kX/m1jEj6dikI0bVKQ1SREgEfzzlj/x+wO95f8X7PDb7MeaNupnnCv7I57eN4T9XTuCsLiMICVRbt4g0LErYpMZYa9m4ew8PP7+G/365jJzoNBqftIRDkWkUBmW6GwUQuOdECpdeDFsGwpZTiCjqwP8+NQwa5G384j+CA4O56MSLuOjEi1i0fTE3vvQKM/PeJOW99wg3MQxtcQ4T+o9ldNczCA8O9zpcERGfU8ImZdq/H7Kzj17WuDE0Ci8iIzeD2Vtns3H/RtZmrGXVnjUs2bqG9fvXkMN+Z+OREBnUmOQWyfSMv4ye8T1Jjk+mR/MeRIREsGcPrF4N69bBoEHQsWPtn6PUDb0Te/HTX57g+Zf+wa1PfcmhNh/y+aGP+fzn/2KmhNMi9zROano6KSeOZMwp3YiJUe2siNQ/StiE/fvh669hweI85i7fxdINO9idvQMid0DjLRC9GWI2OffRWyAwD2b/8nxzoA12TxfIGE/H2M5cd14Xzh3SjbbRbctt2mzWzLmdemotnaTUeVdfGcJVE89hw4ZzmLcwn4+XfM9PGR+xNfQrpub8j6nz4IpvE4jYcxqdQk9mYOuBnNWrJyefFKJpQkSkzlPC5pHsbFi40Pk7NNSZfiI01BkJV1T0yy0urmpzUm3dCrNnw779RezJPMiuzL3sz91LdGIGca33EhG3l8z8vezO2suSdbtZuXUHe7J3YCN2QHgG9MK5uQyG6IAWxAa2IbLwJOz+c9mbHkfOzhM4uKUtPVt3ZNCAcAaOhoEDoXXrGi0mkaMY41wLtkOHYC4YNwIYAcD8dZt4c+Y3fLvpK9aG/cCS4LdZUgjPzgmFj/rSJnAAKQN6cemZySQndPX2JEREqqHBJ2wHDkBGhnPpopK3wMCKH4eEUOURaxkZ8OmnMG0afPklHD5cYqUpgqAcCMp27oOzIfgQAY0yGXpGJqPHZhIdl0lmXiaZuUff7zmYSfrPmWzPyCQr/6CTeIXtg4AiZ98hQIZ7K5YbBYebEpSbSIfYLiR3GEKPdgm0ikkgIdK5xUfEkxiV+KsO3qmpqQwbNqzqhS3iI/06taVfp4nARAC2HtzKt2tm81nabGYGzWZzwTP8Z2sO/3kJjA0kIagtg3b258SEbnRu2pnOsZ3pFNuJ6DBdu0xE/FO9S9hyCnOYs3UORbYIi3XurT3yODOriKXLLIsWF7FkiWXdemc9pgiMe497H5gPAfll3gcEFdA4Jp/ImHyiGucTHuXcwsLzCYvIJySsgINZ+ezam8Pegznsz8omM8dJyILaZNPoDznEhGaTb3PILcqmwOaVeT5FwAxgxuyjlwcRSlBhFDY3ityDUZDbmMiQOHq06EjXNk1JiI4lvnEsCTGxNIuIJXdfLNvSY9m0simrlzShaUwwF18MI0dCUL17F0hD16pxKy7rN47L+o0DICevgJenruXFT5ayaFsa2+PTmJIxiykx7xz1vLjwuCMJXJvoNrRu3JrW0a1pFtyGZiGtaR4dRVhY1f9ZExE5XvXup3rz4c2c/NLJx96wEXCye6uGImC/e6Mw2LkdDIZ9wVAUDEVBzrKCMIIIIzKsEe2bNaJlfCyJzRoRFhRGoyD3PviXxyX/Dg8OJyo0isLsKKa9F8Xbr0ZxcE8U5EZRUBRMbHPnepmDB8MFF0C3bscIemD1zlWkrgsLCeL687ty/fld2bbtfF56aRl5eT2Yv+Aw89avJ6NoHTRdy8HEtaxuu5ZlUV+TabdjTdHRO8qJhgOtMVktCbfxdG8bz+A+zenZMZ7mEc2Jj3Tu48LjCA4M9uZkRaReqncJW8tGLXlh/At8/bXh0X8F0LqVoV27AKIiDZGRATSOMnQ7IYATTzSENwrAYAgwARjj3pd4HBwQTHBgcJn3QQFBBAcGE2gCMcZgrdN5f+tW5/b/7d1rjFxlHcfx7/+c2dmd7V66vReUghEIQqEqqRBj4o2k5QXESAK8EDEajJeIiYaoMRqJL9QXGDUKwZyGtusAAAr4SURBVEtEYwADWKtiFCNU30i5yE1RUySklWYL7badbWd355zz98Vzpp1OZ9uB6ezMdH+f5ul5zjnPmXn677Ob/8y5PLt3w9q1sGEDrFnT/lPaN18E3/o8/O534fUuvBBWrTo1MRNZTM44A971rlcJZ/WHgfW89NJ6tm0jlC3wwgswurTKmy5+mTXn72T8rJ3MFHfyaryTvUM7ObDyfxyoPs92m2T7s7Pw7PHvMz44zkRpgomhUAbSCWx2gmGbYCRexpJogpJNMBxNMMgYg4xQZJSV4yNcdN4o5795gMHBhY2NiPSu0y5hW1JYws6HN3PbJ+GKK2DLFhhegMc0mcHERCjr13fmPcbG4PrrO/PaIovZunVwww2hAJTLMDIygNk6YN28x01OOnfdU+bnv5rkuf/ugZFJCuN7WPWmSZaduZdZm+LF6Smerk6RFf8Zri0tTUFh9sQdegZIikRJSOJK8Qhjg6MsXTLC8tERlhRGibNhPBmCuRLZXIlstkQyM0S1UqJ6uEQyUyKZHSKplEgqJYaLJZaPD7FqosSqiRKrlw+xZmWR1SsHWLUyYvnycONTHIfLJAqFozdCiUj39XzCZmabgO8AMfAjd//Gidrv31/k4x+HK6+E+++HIc1qJCKv0ehoa+1WrzZuuXmMW24e48UXz2X7dnj8cXjiCXjit+Fat0suOVrOOy/csDTnFQ5nU5STKQ6lU1TSMpWsTCWd5pUDZXbtmeblvWVemZ5m73SZgzPTTFmZl4rTUJyEYhkGDoebkwoz4fpaCL8lR/JyMrUbkf4FpAVIi+FyjrR4pFhWpBgXGRoYYKhYZHiwyGChyNBAkVIxbI+8iCehpNUBkmpMMltgal+FoV8/hKcFRkoFRkdixkYLLB0tMFCIia1ARIHIYmIK+XpMZAXifHtU216rN2w3YiIijAiIyJKIpBqRJjFJNWJsNGLtmogz1kasWR0xXIryMykRnkVkaUQhiinEEXEUHVlGkbJU6T09nbCZWQx8H7gC2AU8ZmZb3f2f8x2zZ88gH/gA3HNP+MUoIrIQzjknlGuvDevuYdn8G6pSXs5o+fX37QsPmn7hhfAt2NgYjI+HMjySMLhkhniwwmw6QyWpUKlWqCQVZpKZI/VKNaxPz1bYd3CGA+Uq+6fnOHhojvLhOarpHNVsjmpWJfE5Ds2E7dPlOV6dmSP1KsSHId4fnsfYWKIkLyksq6tDuPD3QF4WyiSw43Uem0XgjSUOSyKsbrvRvO5u4Ia7ceSPWX4JzvzLWt0sHAPNl3NzVYp/HMSz8B7uFurZ0ToYUWTEeYmi5q+LW/7cTMPq6uR9x8J2aNauxhpq+d927D7ydz6mXr/bre7nxo490hqOm/f9Gvtz7PtHkYVvk2OI4/B+SQLVKiRVo1oNj9aanq5Q2vIA7uFn2jMjy+u4NTxNwogjiGJCrGOOWa/VB4t25FFel695D63q6YQN2AjscPf/ApjZPcDVwLwJ2+howr33woCu9xWRLjrVpxKXLYONG0M5XoHWv1p7/aanw6OQaqVcDmcxxsbCt5K1MjQE27aFx//U7tIvH0rYPZmyezJhLklIPSX1hIyELK+H9ZQsr6eeklFXb7IdnIyMzDPiOIMoJSpkRFHGocMZU/sz9k1l7D+QcagStkdxhsUZFqXhaD++5K9KRqinlh7Zdsw+y9tbvu4ZGSkWZZg5FjlmTuZOljlp5qRpWGa1ZXZ0mWShbZo5YSbmsAxPL6hbNyfLBohjwntEjllWVw9Lau+dHO1D7fj5XreV7Uf7lrO6OuRtm9XJX3eefW28Jhz9oHTssU3aNm2XtzWwQl5fVveznLe1Zu9Xi4aHfe6EDykZkMz/b7rv7hKt6vWE7UxgZ936LuAdjY3M7Cbgpnx1tli05xagb4vZCuDVbnfiNKb4dp5i3HmKceetSBXjjqhLq1Z4R2P8dTjRhbJ1ej1ha/YZtTEdxt3vBO4EMLPH3f3STndsMVOMO0vx7TzFuPMU485TjDuvl2Lc649/3AXUT3b0BuDlLvVFREREpCt6PWF7DDjXzM4xsyJwHbC1y30SERERWVA9fUrU3RMz+zTwB8IN6z9x93+c5LA7O9+zRU8x7izFt/MU485TjDtPMe68nomxuR93SZiIiIiI9JBePyUqIiIisugpYRMRERHpcX2bsJnZJjP7t5ntMLMvNNk/aGb35vsfNbOzF76X/auF+N5oZq+Y2VN5+Vg3+tnPzOwnZrbHrPlzAy34bv5/8IyZvW2h+9jPWojvu83sQN0Y/spC97HfmdkbzexhM3vezP5hZjc3aaNx3IYWY6yx3AYzGzKz7Wb2dB7jrzVp0/Wcoi8TtropqzYDbwGuN7O3NDT7KDDl7m8Gvg18c2F72b9ajC/Ave6+IS8/WtBOnh5+Cmw6wf7NwLl5uQm4fQH6dDr5KSeOL8Bf68bwrQvQp9NNAnzO3S8ALgM+1eR3hcZxe1qJMWgst2MWeK+7XwJsADaZ2WUNbbqeU/RlwkbdlFXuPgfUpqyqdzVwV16/D3if2ameLOa01Up8pU3u/hdg3wmaXA38zIO/AUvNbO3C9K7/tRBfaZO773b3J/N6GXieMENNPY3jNrQYY2lDPjan89WBvDTekdn1nKJfE7ZmU1Y1DuAjbdw9IUw5vHxBetf/WokvwAfzUxz3mdkbm+yX9rT6/yCv3+X5aZDfm9mF3e5MP8tPEb0VeLRhl8bxKXKCGIPGclvMLDazp4A9wEPuPu847lZO0a8JWytTVrU0rZU01UrsfgOc7e4XA3/i6CcPOXU0hjvrSWBdfhrke8CWLvenb5nZCHA/8Fl3P9i4u8khGsev0UlirLHcJndP3X0DYUaljWZ2UUOTro/jfk3YWpmy6kgbMysA4+j0SKtOGl933+vus/nqD4G3L1DfFhNNzdZB7n6wdhrE3R8EBsxsRZe71XfMbICQSPzC3R9o0kTjuE0ni7HG8qnj7vuBRzj++teu5xT9mrC1MmXVVuDDef0a4M+upwS36qTxbbgG5SrCdRVyam0FbsjvsrsMOODuu7vdqdOFma2pXYNiZhsJvw/3drdX/SWP34+B5939tnmaaRy3oZUYayy3x8xWmtnSvF4C3g/8q6FZ13OKnp6aaj7zTVllZrcCj7v7VsIA/7mZ7SBkwdd1r8f9pcX4fsbMriLcwbQPuLFrHe5TZnY38G5ghZntAr5KuNgVd78DeBC4EtgBHAY+0p2e9qcW4nsN8AkzS4AKcJ0+1L1m7wQ+BDybX/8D8CXgLNA4PkVaibHGcnvWAnflT0iIgF+6+297LafQ1FQiIiIiPa5fT4mKiIiILBpK2ERERER6nBI2ERERkR6nhE1ERESkxylhExEREelxSthEREREepwSNhEREZEep4RNRGQeZna2mT1oZv82s/+Y2Re73ScRWZyUsImINGFmEWH+xjvc/XxgPXCpmd3U3Z6JyGKkmQ5ERJows83Ax9z9g3Xb1gLb3P287vVMRBYjfcMmItLcBcDT9RvyScvHzKzYnS6JyGKlhE1EpLkUGKnfYGYGDANJV3okIouWEjYRkeYeAa7Mk7SaK4An3T3rTpdEZLFSwiYi0oS7Pw38HbgVwMxWA7cBX+pmv0RkcVLCJiLShJl9AbgU+LKZvRe4HVgH/MDMzu5i10RkEdJdoiIiIiI9Tt+wiYiIiPQ4JWwiIiIiPU4Jm4iIiEiPU8ImIiIi0uOUsImIiIj0OCVsIiIiIj1OCZuIiIhIj/s/6ljOH+B4QWAAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 720x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "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",
    "\n",
    "N = 8192  # length of signals\n",
    "w = 8  # wordlength for requantization of multiplications\n",
    "\n",
    "\n",
    "def uniform_midtread_quantizer(x):\n",
    "    # linear uniform quantization\n",
    "    xQ = Q * np.floor(x/Q + 1/2)\n",
    "    \n",
    "    return xQ\n",
    "\n",
    "\n",
    "def no_quantizer(x):\n",
    "    \n",
    "    return x\n",
    "\n",
    "\n",
    "def sos_df1(x, a, requantize=None):\n",
    "    y = np.zeros(len(x)+2)  # initial value appended\n",
    "    for k in range(len(x)):\n",
    "        y[k] = x[k] - requantize(a[1]*y[k-1])  - requantize(a[2]*y[k-2])\n",
    "        \n",
    "    return y[0:-2]\n",
    "\n",
    "\n",
    "# cofficients of the SOS\n",
    "p = 0.90*np.array([np.exp(1j*np.pi/3), np.exp(-1j*np.pi/3)])\n",
    "a = np.poly(p)\n",
    "# quantization step\n",
    "Q = 1/(2**(w-1))\n",
    "\n",
    "# compute input signal\n",
    "x = np.random.uniform(low=-1, high=1, size=N)\n",
    "# compute output signals w and w/o requantization\n",
    "yQ = sos_df1(x, a, requantize=uniform_midtread_quantizer)\n",
    "y = sos_df1(x, a, requantize=no_quantizer)\n",
    "# compute requantization error\n",
    "e = yQ-y\n",
    "# Signal-to-noise ratio\n",
    "SNR = 10*np.log10(np.var(y)/np.var(e))\n",
    "print('SNR due to requantization: %f dB'%SNR)\n",
    "\n",
    "# estimate PSD of requantization error\n",
    "nf, Pxx = sig.welch(e, window='hamming', nperseg=256, noverlap=128)\n",
    "Pxx = .5*Pxx  # due to normalization in scipy.signal\n",
    "Om = 2*np.pi*nf\n",
    "# compute frequency response of system\n",
    "w, H = sig.freqz([1,0,0], a)\n",
    "\n",
    "\n",
    "# plot results\n",
    "plt.figure(figsize=(10,4))\n",
    "plt.plot(Om, Pxx/Q**2 * 12, 'b', label=r'$|\\hat{\\Phi}_{ee}(e^{j \\Omega})|$')\n",
    "plt.plot(w, np.abs(H)**2 * 2, 'g', label=r'$|H(e^{j \\Omega})|^2$')\n",
    "plt.title('Estimated PSD and transfer function of requantization noise')\n",
    "plt.xlabel(r'$\\Omega$')\n",
    "plt.ylabel(r'$Q^2/12$')\n",
    "plt.axis([0, np.pi, 0, 100])\n",
    "plt.legend()\n",
    "plt.grid();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Small Limit Cycles\n",
    "\n",
    "Besides the requantization noise, recursive filters may be subject to periodic oscillations present at the output. These undesired oscillations are termed *limit cycles*. Small limit cycles emerge from the additive round-off noise due to requantization after a multiplication. The feedback in a recursive filter leads to a feedback of the requantization noise. This may lead to a periodic output signal with an amplitude range of some quantization steps $Q$, even after the input signal is zero. The presence, amplitude and frequency of small limit cycles depends on the location of poles and the structure of the filter. A detailed treatment of this phenomenon is beyond the scope of this notebook and can be found in the literature."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Example - Small limit cycles of a SOS\n",
    "\n",
    "The following example illustrates small limit cycles for the system investigated in the previous example. The input signal is uniformly distributed white noise till time-index $k=256$ and zero for the remainder."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAm4AAADiCAYAAAD3RUuZAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3XecnWWZ//HPdaYmk0YIBEgCCRA6CWCkKOBQxIAFV10XXAvqiiK6uuIPKa66K6irrLqsbVERsICIKE3pjNICJAFCKEmGdEIgnUySKWfO9fvjPGfmzJlTnlPmlJnv+/WaV2aedt/nzhPm4rqbuTsiIiIiUv0ila6AiIiIiISjwE1ERESkRihwExEREakRCtxEREREaoQCNxEREZEaocBNREREpEYocBORqmZm55nZIwXee7CZPW1m283sX0tdt2phZs+bWWsZymkzs38Z6nJEJDMFbiKSFzNbaWanV7oeIV0MtLn7WHe/upQPLnU7FPM8dz/c3dtKVRcRqV4K3ERkONsPeL7SlRARKRUFbiJSMmb2LjN7xsy2mtljZjYrOH6Jmd2Scu3/mNnVwffjzeyXZvaqmb1iZleYWV3IMt8TdBVuDbryDg2OPwicAvzIzDrM7KA09+5jZreb2WYzazezTyWdu87Mrkj6udXM1gbf/xrYF7gjePbFZjbdzNzMzjezdcFnuajQ56Wp6yQzuzP4nJvN7GEziwTn+rJ1ZjbKzK43sy1m9mJQt7VJz1lpZl82s0Vmts3Mfm9mzcG53YIyNgT332lmU8P8PYhIeShwE5GSMLNjgGuBTwO7A/8H3G5mTcCNwFlmNi64tg74IPC74PbrgShwIHA0cAaQcyxVEIzdCHwR2AP4C/Hgp9HdTwUeBj7n7mPcfWmaR9wIrAX2AT4AfMvMTstVrrt/BFgNvDt49neTTp8CzAw+wyVhuj9zPC/hoqCuewCTgcuAdHsWfh2YDuwPvB34cJprPgjMBWYAs4DzguMR4FfEM5X7AruAH+Wqv4iUjwI3ESmVTwH/5+5PuHuvu18PdAHHu/sqYCHw3uDaU4Gd7j7PzCYDZwJfdPcd7v468APgnBBl/hNwl7vf5+49wFXAKOAtuW40s2nAicBX3L3T3Z8BfgF8JJ8PncZ/BJ/jOeJB0LlFPi+hB9gb2M/de9z9YU+/2fQHgW+5+xZ3XwukG9t3tbuvc/fNwB3AUQDuvsnd/+juO919O3Al8LYS1V9ESkCBm4iUyn7ARUFX3lYz2wpMI57Ngnh2LRHEfIj+bNt+QAPwatJ9/wfsGaLMfYBViR/cPQasAaaEvHdzEKAkrAp5bzZrUp63T6YL8/Q9oB2418yWm9klGa7bJ6UOa9Jcsz7p+53AGAAzG21m/2dmq8zsDeDvwISw3dYiMvQUuIlIqawBrnT3CUlfo939xuD8H4DWYMzUP9AfuK0hnpmblHTfOHc/PESZ64gHfgCYmREPFl8Jee9EMxubdGzfpHt3AKOTzu2Vcn+6bBdB+cnPW1fk8+In3be7+0Xuvj/wbuBLGbp1XwWSx6VNS3NNJhcBBwPHufs44OTguOXxDBEZQgrcRKQQDWbWnPRVD/wc+IyZHWdxLWb2zkRg5O4bgDbi3Ycr3P3F4PirwL3Af5vZODOLmNkBZhami+5m4J1mdpqZNRAPPLqAx3Ld6O5rguu+HXyGWcAngd8GlzxDfFzeRDPbi/g4umSvER9Hlurfg8zV4cDHgd8X+Tygb+LHgUFw+gbQG3yluhm4NJhoMAX4XKZnpjGW+Li2rWY2kfh4ORGpIgrcRKQQfyH+Cz7x9Q13n098nNuPgC3Eu/XOS7nvd8Dp9GfbEj4KNAIvBPfeQnw8V1buvoT44Pv/BTYSz0S92927Q36Oc4kP5F8H/An4urvfF5z7NfAssJJ4YPn7lHu/DXw16N79ctLxvxH/7A8AV7n7vUU+L2EmcD/QATwO/CTD2m3/SXwSw4rg+luIB7Nh/JD4GMGNwDzg7pD3iUiZWPqxrSIikg8zm048WGpw92hla9PPzC4AznF3TTIQGQaUcRMRGUbMbG8ze2vQ5Xww8e7jP1W6XiJSGvWVroCIiJRUI/FZuTOArcBNwE8qWiMRKRl1lYqIiIjUCHWVioiIiNQIBW4iIiIiNWLYjnGbNGmST58+fUjL2LFjBy0tLUNaxnCltiuO2q9warvCqe2Ko/Yr3EhouwULFmx09z1yXTdsA7fp06czf/78IS2jra2N1tbWIS1juFLbFUftVzi1XeHUdsVR+xVuJLSdma3KfZW6SkVERERqhgI3ERERkRqhwE1ERESkRihwExEREakRCtzKrKMrytUPLCPaG6t0VURERKTGKHArsyvufIHv37eUB156vdJVERERkRqjwK3M1m7ZBUBzQ12FayIiIiK1RoFbmW3Z2Q3AmKZhu4SeiIiIDBEFbmW2dWcPAO5ekfJ7Y84vHl7Ozu5o2vOLX9nG9Evu4unVW8pcMxEREclFgVsJdUV7aX99e8bzy17bzitb412l0VhlArfbn32FK+56kf99sD3t+cdf3gTAnYteLWe1REREJAQFbiX0739ezOnf/zsbO7rSnn/7D/7e932sQoHbuq2dWcsfP6oBgG27egp6frQ3xgvr3iiscml0R2MVaysREZFqo8CthB5fHs9WdXSm74ZM1uvOru7eIanH3YvXc9p/t7FlR/egc4mAbPzoeIB23aMrWL6ho+/8mOb6vute397JX5/LL/P2rb+8xFlXP8yqTTsKrf4AB331r3zp5mdyXvfjh9p5OelziIiIDEcK3ErIsNDXti3ZwKFfu5v5KzeXvB6f+c0CXt6wg+UbdwwaS7c1mBzRXF9HdzTGN+54gX/82eN95xOXb9vZwznXzOOC3y6ksyd8gJkIXreHCF7D+vMz67Ke39Xdy/fuWcJp//23kpUpIiJSjRS4DYEw49d++cgKAK59dEVJyozFnP+443lWb9rZd+z9P32MWxe+MvC6oGpd0Rhd0XhAtjkI5rqivSxZH+/m3NkTZfmGeNYsn/F4XUGQV18XPojNJOwixdGYFjMWEZGRQYHbEOiOhg8kOnvCXbtg1WaWvpZ54sMLr77Brx5dyYW/Wzjg+ENL+hf6jcWc17d3BeX20hXUM5Flu/KuF7k6mLSw+JX+cWr57PKQeGZPtPhxaV0h21Fxm4iIjBQK3ErIiQcr3UGgc9OTq3msfWPWe+oi4TJT7//p45yRNLkhVcTiz+lJCbImj2vu+/7qB5fx96UbAOiM9g4KjF56NX1g+I3bnw/dXZq4rru3+PF7YcvsrdDSKiIiIuWmwG0IfP22xezsjnLJrc/xoV88kfXahhBdits7c8/wjAR/k6kxTLQ3xoJV8TXZ2pZs6Ds+f+UWbnpy9YBrp04clfbZf35mHX9cuDZnHaA/SxY2W5bJ9Evu4st/eDbUtTEFbiIiMkIocBsCz67dxtdue37AsYtvSR+E1EVy/xWs3Lgz5zWJjFtq9un6x1fx/p8+BgzMYC1YtWXQWm7J2blU0d70wdG2XT2850ePcNsz8bF0iXFzie7iNZt38v37lha04PBDSYFmNlouRERERgoFbiVy56J1rNm8q+/nWxb0Z6g6e3q5eX76jFVDxDj2yvv54f1LMz57R4ZdDtLJln3KNfauIUu3baZTf1u6gUVrt3Hz/DVAf8YvUdbnb3yaqx9YxtLXwi/VkW+Qp65SEREZKRS4lcjnfvd0xnPZuvzMjNe3d/HD+5f1ZcZSJQdcmYKaREYsMRM0nVxjxrIlrsz6I7d7nl9P++vxQKwnqFt9kDlMjNlLjPPrDR6az5IiPRmye5n0KuMmIiIjhHY6L4OnV2/NeC4xoQHoG4uWKjlw29XTy+jGwX9tYYKXXOPOwmauPv3rBQA8esmpPPDSa0B/pq+hLkJXNNZX56b6SKiyB9Yzv4kNmlUqIiIjhTJuZZAtqAoTcHUnzRT92d+Wc/8Lrw26JsxaZrmyXtnqku7et37nQf7y3HqgP3BLrN/WHY3x7JqtzA+C0VzB2NotO7n92XV996ba1NHFkd+4hzueHbwY786e0i32KyIiUs1qJnAzs7lmtsTM2s3skkrXJx/ZFrDNNe5s4eotfPa3/WuzXf3AMv7lhvkDrtm2sydU92JnlrIeeun1rIHbzmB7rkxruiXuTXSZdvfGOPvHj/ad78qxXt0//uxx/vXGp4nFPG127tVtnWzvjPL5Gwd3Sc/94cNZny0iIjJc1ETgZmZ1wI+BM4HDgHPN7LDK1iou2hvLmU3KlunK1YX4n3e8kPX8lh3dzP7Pe7nq3iVZr3P3rIHZx697Kuv5XcFnyBT8Je5tSMq4pbs/k1e3dcbv642lbZNMky4Kma0qIiJSq2pljNuxQLu7Lwcws5uAs4HsUU0ZnPvzeTy7ZlvWazq6MnflpQZ91z6ygq27evjiaTOJRIzdgs3gM9kSbFf15Irse56GiW+yBW6L1m7l4WUb+Mgvn0x7PnFroqs0NfjamWNmbF3E6I05tz+zjubGutB129jRnfW5IiIiw0mtBG5TgDVJP68Fjku9yMzOB84HmDx5Mm1tbUNaqe3bO3hqZXF7cm7YNHBCwn/eGY9FW7av4eCJdXQFW1SlSny2DTvDjcx/KERbrFn7SsZzj7Zv4tH2TRnPL1i1hcuuu4/uYLHgl9qXDzj/7PNLmLyj/1hHR0ffZ0jOBl78x0Vpnz9/QX93cfLf6/KtAwPfof47rxbJ7Sf5UdsVTm1XHLVf4dR2/WolcEsXHQ1Kwbj7NcA1AHPmzPHW1tYhrdTd9z8E5F4cN5uN3fXA4KzRrNlHccIBu/Nwxws8tm7wRvSrm6bz0ROms2bzTvj7Q7kL2vswYH7WS57bUlzP+e9e6mbSmCagizteHrjbw5T9ZtDaemDfz21tbbS2tuLunPk/2ceo1UeMWbOPgifmAZD899qycjPMe7zv54OPPo69x6ffAWI4SbSf5E9tVzi1XXHUfoVT2/WriTFuxDNs05J+ngoMnl5YZp0l2Eh98470XX2JZUIydZV+7bbnWbVpB8+sybzUSLJPXp89aAPYlKEu+djYkT5DmLqHasKS17bz0vr0e6QmxNwzLlWSumvCCd9+kIeWvB6ipiIiIrWnVgK3p4CZZjbDzBqBc4DbK1wnOovfRz2npvrB470S3va9trSzLKtRpi2zdnTlbsSYw4d+nn7P13RD3/61RtpEREQkXzURuLl7FPgccA/wInCzuz+f/a6hV2zGbWxT7p7qSmygPmlMY8mfmWlJlEyZuLDSzSrd3ql13UREZHiqlTFuuPtfgL9Uuh7JdhUZH4QJySqxm9NQxIqZ1n/LtUzIAXu08HKWbby025WIiIwkNZFxq1adee6pmSprNs0Tf+RXxv97x8FF1Kg0Tj5oDz510owBxzJl3DpyZMcmjM6c/fvz06/w4V+m70IVEREZjhS4FaHYHrlscVsizsk3+9WSZg20fBWbxJoyYRRvOXDSgGPRWIzemBPtjfG7J1azpTOegcvVrZnYtD6dHz/UnvHcsyEnbYiIiNQSBW5F2FXkGLds2bT/eWApHV3RvlmTf7zghFDPbAkxbi4fl555CI31+b0mEYNRDQMDyGiv848/e4wDL/8rl/3pOf736fjs046unnSP6FOfJnBLjGuLWOagLnm7LRERkeFCgVsROnqKDNyy3P7Uyi1cdc+SvtDu0L3HhXpmKdYwu+ETxxZ1f8SM0SmZv2jMWbi6Pwv2Rnf8kyX2QM0kXcYtsVhvlrhNRERkWFLgVqC1W3Zy67Ls2aJccoV9HV3RvnFwDXXh/qqOnDKea8+bw33/djKzp44vqF5HTBnP+SfvX9C9EM+4DQrcUiYnuMN/3f1S1okHkD5we2jJBnp6YznXfxMRERluFLgVaNlrHYOOffPsw/N7SI7IzUjaAzTLWK9kdXXGqYdMZubksYxpLr7btJCcopkxqnFg2T0pkxM2dTo/bXuZO57Nvo5yus/9qRvmc+vCtQXUTEREpLbVzHIg1SZdUGR59t3lmjH6hwX9wUnYZycHOpZ2p7AC5Bm9RcwYnTLGrbfAGbiZJifsytHFCrBk/XYO3mtsQeWKiIhUI2XcCtQTHbwuWepg+YktmZeyuP9LJ/eNcTsxZQZmMZIDnVyx3o8/dAxfe9dhJSs7IWIwatAYt3AL7R6wR8uAn+sj6V/RMGHgO37491BlioiI1AoFbgXqDsZsfed9R/Ydq4vA3V88Kee9s6eO58A9xyZNPCg+KzR999HxOuSR9cs2bK6YXJ0ZNNVHOPfY/u1lc01CSOjsGRjgRUJ2EYuIiIwE6iotUGLvzeTMkplxyF65Z38mArbEshZ1GbJK+fjJP7+J8aMbBgQ6ubpXI2ZDMjPTPV72t983ixufXAPAjpCB287ugeu6ZRrb9+SKzaGe1xvzrGvBiYiI1BJl3AqU2GMzeb2yfLJdkP/Eg2xamuqYMmHgUiC5nloXsZwL/JZq+6vOLIFbc0P/a5i6BVamtdr+unh9qHK37Spu5q+IiEg1UeBWoERXaXLGLTVxlm4D9HRKkRFK94xccWTYcvPddivd1dn2JG2q72/D1K7SYoPazTu6i7pfRESkmihwK1BPoqs0KeOWmh0KG+6UIuOWaRB/NukCt+9+YFbRdUmOV3/6z8fQWB+hM2vg1l/3E/bffWAd68JP+Ejn9O//LeMG9yIiIrVGgVuBEsFAc5rA7cGL3saTl50W+lmlGIBfX5f/M9J17X5wTjChoIgqJWfozjxyb046cFL2jFtSV+kvPjaHBy96W9/PqUHtZ1sPyLs+m3cq6yYiIsODJicUqCddV2kQCO2/x5is9x6Qcr40Gbc0XaU57skWMI5rbgBgTFP+m9an9hC/tH571s3kG5Omt7Y01Q9ov9SsYFND/vXZsqOHPcc2532fiIhItVHGrUDdabtKc983Z7/d+NY/HDnw2PSJAHz1nYcWXJ/0Y9yyVyjbGLdPnbQ///6uwzj32H2LnqDwytZdWc9n2yw+NSvYnOeG9wAbO7ryvkdERKQaKXArUCLjNnrA5ISUMW5pAp6TD9pj0OK0k8c1sfI77+SDb542+IaQ0gZuOe7Jdr6xPsInT5xBfcg9UpOFnZSR0J1lDFrqGLdGBW4iIjKCKXArULYxbtmkuyIxsWBccwM//+icguqTruxc1RmKNdxg8KSMSWOasl7f1ZM5cEvtAm7IEEgmLymS6gs3PcPqTTuz1kFERKQWKHArUHevYwycERkmOZUuWEqeEPr2wybnXZc3T99tQACZVFrOe8c0lX6YY2rC7Z4cu0l0RTNPXEhdnDhT4Pbwxaeyx9jMAeIdi7JvZi8iIlILFLgVqKc3Rp0NHEeWOqYsXZdhunFn+S7cm+ptB+2R9vjEloac977/TVP59Mn7F1V+qtR133bPlXFLs+9rwuCMW//PHzthv77vc40v3KL13EREZBgoKHAzs0mWa+T7MNcTjZE63Cq1uzJ9Fmyw1DXYPn/qgXnVJdNfxdfefTgnzcy2gb1RFzE+87bsS2zkOzch38kM2QK31LF7yRm3d87ap+/7XN3UWohXRESGg5yBm5kdb2ZtZnarmR1tZouBxcBrZjZ36KtYnXp6Y4O6RlMzZzeefzxnH7XPgGO5ukoBLjrj4FJUkTFN9fzLSZmzaYm6NBQw4L+UemNOQ53xwTlTB53LFrglZ99yBW6bFLiJiMgwEGaA04+Ay4DxwIPAme4+z8wOAW4E7h7C+lWt5oY6xjcODBZSu+sO2GMMl555KLc9k318VbFbXmWLWcI8eVRDfJ/TL739oKLqkRA24dZQZ307UCy78qy01ySC4SOmjOOGTxzHyk07+s4lzzC1HLHnwtVbiPbGCpolKyIiUi3C/Bard/d73f0PwHp3nwfg7i8NbdWq26VnHcq3Tho94Fi6LsuWlAVsLU0olS5btOCrp/PwxaeEqku6ZyaECaLqIsajl5zK+980OONViLBdpYlFfrNJBLXN9XVMbGkcsFhv8ve5Mm7bO6N87fbnw1VMRESkSoUJ3JIHIKWupFrk0qzDS7rMWUvjwKRmuvgi3a4Hu49pYreQ+3JW32jDcK/FuFG5A7dEd2jiMyZv7dUwIHDLXd49i9eHqpeIiEi1ChO4zTazN8xsOzDLzLYn/XxkrptHknTBQ+qivOnii0xdpSXYCSvrYrj7TRyd8VzYZyRLBFdhM25jm3P31CfaL5HNTA7WkrtKI2a8Z/bA8YSpNM5NRERqXc7fnO6e/+aQI1SYzeLTZccyzQrN1v13xXuPoDsaY/6qzXzouH1D1zFh5Xfemfc92XzzvUfQ1dPLFXe9mPZ8fcSIxgZGdM31uV+tRDYy0RKNAyYnJI1xM7jsrEPpivbym3mr86y9iIhIbcgZuJnZl7Kdd/fvl646tS3czgnh02jZnjd5XDNvP2wynzhxRtZnDFVf9o8/dAzRWIwv3PQMAB85fj9ufDIeMKXLuDU31NHRNXCj+dTZtOkk2ivRFsldpakZt7qIMT5H92tPbyzjIr4iIiLVLsxvsLHB1xzgAmBK8PUZ4LChq1rtyZRwGzD7Memar8w9JOvOBdkSeJUc1vbXL5zEO2ftzdlHTRlwPFGn1AV4If2WVIku4tMP3TNnmX1LlxQ4OSHh4796KtR1IiIi1ShMV+l/AJjZvcAx7r49+PkbwB+GtHY1JlPwcM8XT+aUq9oGHb+g9QAuaM28+G22ZUJCT0gYgpTboXuP6/v+1s++pS/Dlm2M23lvmc5V9y4dcKwuEmHFt9MvA9LHBj67IZJpjFu4uj/SvjHchSIiIlUonz6jfYHk0d3dwPRSVcTMvmdmL5nZIjP7k5lNSDp3qZm1m9kSM3tHqcostUyB24xJLZx7bP7j0LJtThE2w1QKR04Zn/HcMfvuxpv22w3o79ZMFyteeMqBfPmMgevEJbYMC7MJR+LzNtT3X5sc2CaesUeO7bVERERqWT6B26+BJ83sG2b2deAJ4PoS1uU+4Ah3nwUsBS4FMLPDgHOAw4G5wE/MrConTIQas1VEwDUqeQutkI9J122Zrxs+cRyXn3Vo7guzZNzSBWipG8iHkbo9WKqPnDCdfzu9NAsJi4iIVJvQvznd/Urg48AWYCvwcXf/dqkqEizymxi9Pg9IrAZ7NnCTu3e5+wqgHTi2VOWWUpjN4ovJkz3ylf4Fecs5xm386AbmTN8t53X51imfOQJ9Gbe67KXURYz3HJV9WZBHlqm7VEREalOYLa/6uPtCYOEQ1SXZJ4DfB99PIR7IJawNjg1iZucD5wNMnjyZtra2IawidHR0DCjjqaee4pUx6aORV9d1AfByeztt0VV5l5UoZ/dmY1On89yiRfBq7r++Ra9H0x7Pt22Wb+3Nee9La3sAWL9+PW1tWwadf3n5wHXUNm3cmLMeS16Kb9CxZfPmQdcm/5z8/fodmTetB/jwL5/gurktWa+pdqnvnoSntiuc2q44ar/Cqe365RW4FcvM7gf2SnPqcne/LbjmciAK/DZxW5rr0/b/ufs1wDUAc+bM8dbW1mKrnFVbWxutra1w910AnHD8ccyYlD4guGfzc7B2NTNnHkjrW7Mv4TFA8OzEZ2me9yB07mL27NmcfNAeOW/vffE1WDh/0PF822bCmq0w79Gs925csBYWP8vkyZNpbT1q0PlFvctgWf8EhcmT96S19Zj0BQaf+5BDDoHFi9h994m0th474Fxy2yfXacXGHfBwW9bPM9TvxlDre/ckb2q7wqntiqP2K5zarl9ZAzd3Pz3beTP7GPAu4DTvX65/LTAt6bKpQPZd28ussS5Cd28s1MzGYsa4Qf8uBuWcnADhN6sHaMmwxEksqLtZfBxcPp+g1J83FvNQCyaLiIhUk9CjjMysycw+ZGaXmdnXEl+lqoiZzQW+ArzH3XcmnbodOCcofwYwE3iyVOWWQmJZilAL8JYoVgj7nLDbT5WivLlH7MXFcw/mK2cekvb8P7152oAZtvkEsfm0234TR3PusdOyXrN1V0/4B4qIiFSJfKb13UZ8okAU2JH0VSo/Ir7Q731m9oyZ/QzA3Z8HbgZeAO4GLnT33syPKb/EgPlyZnDCljRhdO6N3MOVl7vEuojx2dYDMy4qvPf4UTz05Vam7jYKyG8v1nyCvEjE+Pb7ZmW95phv3he+cBERkSqRT1fpVHefO1QVcfcDs5y7ErhyqMouVmI1/zChRclCu5APmjN9Ij//6ByWvrad792zpPDiShiTJoLAMI88bsbuNNVH+Mzb9i9dBQLd0diARXxFRESqXT6/tR4zsyOHrCY1LBG49cZC9EuWKALKZ8/Ttx82mdGN1bP0XSLTFqZreeKYRpZccSZv2m9iyeuxTd2lIiJSY/IJ3E4EFgS7Fywys+fMbNFQVayWNAVZm65o9mUooHQZt3x7ZeuL3Fi9lBm3RMAWbseE0pWbasGqwUuWiIiIVLN8fpufSXxiwBnAu4nP/nz3UFSq1nz3A7OYs99u7DtxdM5rSzc5Ib8HNRQZAeWT4cv5rJT9R7MZytmzn/nNAp5bu23Ini8iIlJqoce4uXv+q8aOEHOmT+SWC95S1jLzjWeKzbiVUiIYKzSW/NV5bw7V9bvoG2cw6xv3Zr1mxaYdHDk1816sIiIi1SRn4GZmj7j7iWa2nYEL3xrg7j5uyGo3DJVseY48r8+1VVTO8kraVRo8M8SnSJdxO+WQPUOVM64594zazp6qmqAsIiKSVc7Azd1PDP4cO/TVGb5K3eOX7/P22724LZ5KGrhFEmPcyltuOmHGJYqIiFSL6uk/G+ZKlWlLPCbfMW5HTZvAtefNKbjcko5xS/wZanLC0EZuXcq4iYhIDVHgVmYlm5xQwD0H71V4r/ZQxE+htgjLcf7Y6cUtE6KMm4iI1JKy7lUqpVPInqfFxF6ljNv6s4Yhys1xzY3nHx9u/bwMlHETEZFaosCtRhUSSBWTNStlxi3RbRxub9fs19RFjLoiljp5cf123L2gQFhERKTcwswq3Tfks7a6+xtF1kdCKiTOKG6cWukCm1gQuZUjVBrbVM/2rmjG8/e98Bo/fqidz506swy1ERERKU6YjNv1xHu3sv2edeA64IYS1GlYK91yIAVmh30yAAATH0lEQVR0lVZZxq0cWa55l53G4V+/J+s1973wmgI3ERGpCWGWAzmlHBWR/BSWcasOfRm3MlSopSn3/5ts7Oge+oqIiIiUQJiu0jDT9mLuvrUE9Rm2Cg1S7vz8iTQ3DN4loKDArYhIqaSTE/IY41YOGzu6Kl0FERGRUMJ0la4LvrL9lq0Dwo6FkzwcMSX9dkzl7yotXZDllG+MWxhaEkRERGpFmMDtRXc/OtsFZvZ0ieojIZW7q3RIMm5ZZoPeduFbeXz5phKWmtl1H39zWcoREREpVpjA7YQSXTOilWpSQkLZu0pLGLklll3L9sjZ0yYwe9qE0hWaRWOd1qEWEZHakPM3lrt3ApjZFannzKwu+RrJrdgAqJjxYcVl3Iagq7RKxrhly/yJiIhUk3xSDVPM7NzED2a2J3B/6askYQyPBXhL98xiXHXPEv709NpKV0NERCSnfAK3TwPnm9mxZvZm4EHgqqGp1vBVbJdpIoAq/wK8pROLlW85EIBv/cORWc/PX7WFy/+0uDyVERERKUKY5UBuABYCTwMXAr8DosB73b19aKs3fJQ+SCnv7ISSZtyCP8u1HMiHjtuXhau3cMuC9Fm1E/bfnafXbClLXURERIoRJuN2fXDdJ4gHbdOBLcCHzewDQ1c1yaawyQmlr0chyrnlVUK2spobIkR7Szx7REREZAiE2TnhAeCBxM9mVg8cBswGjgduGbLaySAVm5xQynXcEjFSGSPJbGFZc0Md0Zhrs3kREal6OTNuZrYw+Wd3j7r7Inf/tbt/Od01MvQKm5xQJTsnBH9Wy+SExM4UPcq6iYhIlQuzjtuhZrYoy3kD0i/vL0Om7AvwlnRWaTxAqpYtr5ob4v//Eo3FaMxrvo6IiEh5hQncDglxTW+xFRnuSr4Ab7m3vCphzq0SY9yyaaqPZ9yuumcpY5rraaqPDGir3l6npzdGY31kwJpv7tDV00tdJEJD/cBP09kTw4CmhoGBYE/U6Y3FaGqoG1BGLOZ0R2M01EWoq8tdxvLl3TzdsxSnP/AsqoxojIhBY/3AZ3VHY8ScQW0S7XWivYPLSNS3vi5CfV24NumOxojFfNCzhqrdV6zo5gXac7dJhjK6emI44dskUUZqm2Rr90QZqX+3+bZ7ooxStfuS9m6ejS7L2CYleRd7eolELHSblPN9z/Quhn3fly/v5pno0rK+76Vq9yH/b0COdl+5sv/fbbn/O3PucdPYc2wz1SLMGLdVqcfMbCbxjeVfHpJaSU7lXg5kSNZxK2Nf6WfedgCPtm+kN+Z9gePGjm6gv6v02kdXlK0+JbF0WaVrULuWLql0DWrXy0srXYPapn+3havQv9unVm7mN/9yXEXKTidMxm0AM/smMCn4fi/gU+6+sdQVG25K3StY7lmlpaz+Pxw9hZ+0vcxZR+5dwqdmd+CeY3j80tMGHJt+yV0AjAoCNxERkVTbO3sqXYUB8g7cgAnufgGAmU0Hvmdm29z9i6WsmGRXyESDooLHEkZuMyeP5bq5LcyY1FK6hxYhtetFRESkWhXyG6tvPJu7rwRiwK9KVSEJp6BZpVUzqqy6NCvjJiIimVTJRLqEQgK3R8zs22a2j5ntAzS7+7OlqpCZfdnM3MwS3bFmZlebWbuZLTKzY0pVVi0rf1dpdb24paSMm4iI1Iq8u0rd/RYzWwd8M7j/ylJVxsymAW8HVicdPhOYGXwdB/w0+FPyVC3LgVQbZdxERKRWhE41mNn9ZjYbwN0fc/dPuvvH3P2FEtbnB8DFDFzo/mzgBo+bB0wws/KNaq8ynnUPgOyqZQHeapNYDkRERKTa5dNHdDHwAzP71VAETmb2HuCVNN2uU4A1ST+vDY7VpFIt51bQOm7FlDeMU27qKhURkVoRuqvU3RcCp5rZ+4G7zexW4LvuvivsM8zsfmCvNKcuBy4Dzkh3W7rqZHj++cD5AJMnT6atrS1s1QrS0dERuox167oAWLZ0KW2dha8X1tUVX3vssccfY2JzfgGHJ60CnG/bdHQXfm/a5+XRdkPtxcXPVboKIiJSpba/8UbV/L6CPMe4WTztsoT4OLMrgE+Z2aXu/usw97v76RmeeyQwA3g2yOxMBRaa2bHEM2zTki6fCqzL8PxrgGsA5syZ462trWGqVbC2tjbClnHvludgzWpmHnQQrcfvV3CZTY/dD11dnHDCCew9flT+D7gnvnZZvm2zdWc3PHhfQfemk0/bDZm7421x7Jyj4anHK1sXERGpSmPHjaO19a2VrkaffMa4PQK8Qnwc2hTgPKAVONbMrimmEu7+nLvv6e7T3X068WDtGHdfD9wOfDSYXXo8sM3dXy2mvEoqtsOxUrM7h/Os0vqIukpFRKQ25JNx+wzwvPugXTc/b2YvlrBOqf4CnAW0AzuBjw9hWVWvmMkJRRm+cRsNdQrcREQkvWfXbOXlDR0csMeYSlcFCBG4mdm+wbdvANMyDFI/y8zGufsbpahUkHVLfO/AhaV47nAynDNg5fKnz76Fp1ZuZv89Wnjf0VN4cuVmIL70yZimBuojxq6eXjp7+tacZnRjHY31EaK9TkdXtO94Y12E0U11xGKwoztKbyweYNdFjDFN9ZjBzq5euntjfWW0NNZTX2d0R2Ps7O4vY1RDHU0N8TJ2dEf79nZtqIvQ0lSHO2zYup3GpvimxxGLlxGJwM7uXrqjsb5njWmKl9ETjT8roak+wqjGOnpj0NHV01dGfcRoaYr/Z2FHV5Ro8DkiFj9eF4Fd3b10JZXR0lhPQ33mMmIx6OiK9u0Rm1zGzu5eepLaZExTA3WR+EbRmdo9uU2ytXtLY9AmSe0OYNEuxo8bU1C7p7ZJvu3e3FBHc0OE3lj8WcltMqa5Pm0ZiXZPbZNEuw96F+sjjE7T7ol3MV27J97FXO2+a1cnzc3Nebd7oe/76Mb4jO9M7Z76LqZrd/f4BKTE+57a7i3Bv8+Ozvze90LafePW7TQE/26T26SrJ8auIXrfE22S2u6Z3sUw73up2j3d322m952eTiaMH5t3u5vBjq7073vYdj948lgaqqhnJkzG7XrikwGyRQoOXAfcUII6SRUajpNKj953N47edzcAvv9PR1W4NvmpijGCNSredidVuho1Se9dcdR+hdO/2345Azd3P6UcFZHymDSmiQtPOSDv+4Zh3CYiIlJzCtlkXopQoRFqfeZ/Ne3E3pyG8zpuIiIitaJ6Om2HuVKFPd953yxm7jmG3cc0luiJ4ShsExERqTxl3GrMKYfsySmH7Fn2cpVwExERqTxl3CQUzWIVERGpPAVuZVLpsW3FUsZNRESk8hS4lZniHxERESmUAjcRERGRGqHArcxqtctUXaUiIiKVp8CtTGo97tHkBBERkcpT4CahKOMmIiJSeQrcJBTFbSIiIpWnwE1C0ZZXIiIilafATUJR2CYiIlJ5CtzKpFZnkyYo4SYiIlJ5CtxEREREaoQCNwlFY9xEREQqT4FbmSjsERERkWIpcBMRERGpEQrcRERERGqEArcymTGpBYC9xzVXuCYiIiJSq+orXYGR4hNvncFhe4/jLQdOqnRVREREpEYp41YmkYgpaBMREZGiKHATERERqREK3ERERERqhAI3ERERkRqhwE1ERESkRihwExEREakRCtxEREREaoQCNxEREZEaUVWBm5l93syWmNnzZvbdpOOXmll7cO4dlayjiIiISKVUzc4JZnYKcDYwy927zGzP4PhhwDnA4cA+wP1mdpC791autiIiIiLlV00ZtwuA77h7F4C7vx4cPxu4yd273H0F0A4cW6E6ioiIiFSMuXul6wCAmT0D3AbMBTqBL7v7U2b2I2Ceu/8muO6XwF/d/ZY0zzgfOB9g8uTJb7rpppuGtM4dHR2MGTNmSMuoJufdvQOA6+a2FP2skdZ2pab2K5zarnBqu+Ko/Qo3EtrulFNOWeDuc3JdV9auUjO7H9grzanLg7rsBhwPvBm42cz2ByzN9WmjTXe/BrgGYM6cOd7a2lqCWmfW1tbGUJdRVe6+C6Akn3nEtV2Jqf0Kp7YrnNquOGq/wqnt+pU1cHP30zOdM7MLgFs9ngJ80sxiwCRgLTAt6dKpwLohraiIiIhIFaqmMW5/Bk4FMLODgEZgI3A7cI6ZNZnZDGAm8GTFaikiIiJSIVUzqxS4FrjWzBYD3cDHguzb82Z2M/ACEAUu1IxSERERGYmqJnBz927gwxnOXQlcWd4aiYiIiFSXauoqFREREZEsFLiJiIiI1AgFbiIiIiI1QoGbiIiISI1Q4CZ5OfmgPSpdBRERkRGramaVSvVbesWZ1EXSbWQhIiIi5aDATUJrrFeCVkREpJL0m1hERESkRihwExEREakRCtxEREREaoQCNxEREZEaocBNREREpEaYu1e6DkPCzDYAq4a4mEnAxiEuY7hS2xVH7Vc4tV3h1HbFUfsVbiS03X7unnOx1GEbuJWDmc139zmVrkctUtsVR+1XOLVd4dR2xVH7FU5t109dpSIiIiI1QoGbiIiISI1Q4FacaypdgRqmtiuO2q9warvCqe2Ko/YrnNouoDFuIiIiIjVCGTcRERGRGqHArUBmNtfMlphZu5ldUun6VBszm2ZmD5nZi2b2vJl9ITg+0czuM7NlwZ+7BcfNzK4O2nORmR1T2U9QeWZWZ2ZPm9mdwc8zzOyJoO1+b2aNwfGm4Of24Pz0Sta70sxsgpndYmYvBe/fCXrvwjOzfwv+zS42sxvNrFnvXnpmdq2ZvW5mi5OO5f2umdnHguuXmdnHKvFZKiFD+30v+Le7yMz+ZGYTks5dGrTfEjN7R9LxEfX7WIFbAcysDvgxcCZwGHCumR1W2VpVnShwkbsfChwPXBi00SXAA+4+E3gg+BnibTkz+Dof+Gn5q1x1vgC8mPTzfwE/CNpuC/DJ4PgngS3ufiDwg+C6kex/gLvd/RBgNvE21HsXgplNAf4VmOPuRwB1wDno3cvkOmBuyrG83jUzmwh8HTgOOBb4eiLYGwGuY3D73Qcc4e6zgKXApQDB749zgMODe34S/M/tiPt9rMCtMMcC7e6+3N27gZuAsytcp6ri7q+6+8Lg++3Ef3lOId5O1weXXQ+8N/j+bOAGj5sHTDCzvctc7aphZlOBdwK/CH424FTgluCS1LZLtOktwGnB9SOOmY0DTgZ+CeDu3e6+Fb13+agHRplZPTAaeBW9e2m5+9+BzSmH833X3gHc5+6b3X0L8cAlNZgZltK1n7vf6+7R4Md5wNTg+7OBm9y9y91XAO3EfxePuN/HCtwKMwVYk/Tz2uCYpBF0nxwNPAFMdvdXIR7cAXsGl6lNB/ohcDEQC37eHdia9B+05Pbpa7vg/Lbg+pFof2AD8Kugm/kXZtaC3rtQ3P0V4CpgNfGAbRuwAL17+cj3XdM7mNkngL8G36v9AgrcCpPu/yg1PTcNMxsD/BH4oru/ke3SNMdGZJua2buA1919QfLhNJd6iHMjTT1wDPBTdz8a2EF/V1U6arskQRfd2cAMYB+ghXgXVCq9e/nL1FZqwzTM7HLiQ25+mziU5rIR2X4K3AqzFpiW9PNUYF2F6lK1zKyBeND2W3e/NTj8WqIrKvjz9eC42rTfW4H3mNlK4mn/U4ln4CYE3VcwsH362i44P57B3TcjxVpgrbs/Efx8C/FATu9dOKcDK9x9g7v3ALcCb0HvXj7yfdf0DqYIJmi8C/hn71+zTO0XUOBWmKeAmcFMq0biAyZvr3CdqkowzuWXwIvu/v2kU7cDiVlTHwNuSzr+0WDm1fHAtkR3w0jj7pe6+1R3n0783XrQ3f8ZeAj4QHBZatsl2vQDwfXD+v84M3H39cAaMzs4OHQa8AJ678JaDRxvZqODf8OJ9tO7F16+79o9wBlmtluQ8TwjODYimdlc4CvAe9x9Z9Kp24FzgpnMM4hP8niSkfj72N31VcAXcBbxGS8vA5dXuj7V9gWcSDxdvQh4Jvg6i/j4lweAZcGfE4PrjfjMoJeB54jPaqv456j0F9AK3Bl8vz/x/1C1A38AmoLjzcHP7cH5/Std7wq32VHA/ODd+zOwm967vNrvP4CXgMXAr4EmvXsZ2+pG4mMBe4hnfj5ZyLtGfCxXe/D18Up/rgq3XzvxMWuJ3xs/S7r+8qD9lgBnJh0fUb+PtXOCiIiISI1QV6mIiIhIjVDgJiIiIlIjFLiJiIiI1AgFbiIiIiI1QoGbiIiISI1Q4CYiIiJSIxS4iYiIiNQIBW4iIiGY2elm9utK10NERjYFbiIi4cwGnq50JURkZFPgJiISzmzg6WCvxOvM7FvBfp4iImVTX+kKiIjUiNnA68Q3AP+Fu/+mwvURkRFIe5WKiORgZg3ARmAV8Gl3f7zCVRKREUpdpSIiuR0GPAVEgd4K10VERjAFbiIiuc0GHgPOAX5lZpMrXB8RGaEUuImI5DYbWOzuS4GvADcH3aciImWlMW4iIiIiNUIZNxEREZEaocBNREREpEYocBMRERGpEQrcRERERGqEAjcRERGRGqHATURERKRGKHATERERqREK3ERERERqxP8HVQ28gXZm4I8AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 720x216 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmgAAADiCAYAAAD6WzveAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAHahJREFUeJzt3Xm4XHWZ4PHvSwgQ16hEIZfNVp6rNMPSjSitM31V2qBta3TQgVZaXDoy3Wr7DEaMaOMu3fEZHTdG2gXcEAdCsHGJLN6m0YAKQQLCBUYWSbCR5YqBOxiSd/6oc6HqblV1U+fWqVvfz/PUc6vO+Z1fvef9nVN5c5aqyEwkSZJUHTt1OwBJkiQ1skCTJEmqGAs0SZKkirFAkyRJqhgLNEmSpIqxQJMkSaoYCzRJlRYR10XE0By8z3BEvGWaeRERX4mI+yLip2XHUpaIeF1E/LDbcUhqzgJNmoci4viI2BgRD0bEbyLitIhY3Mbyt0bEkR2MZ9b9ZeYfZ+Zwp2KZpRcAfwHslZmHdzmWWcvMb2TmS8p+n4jYLyIyInYu+72k+coCTZpnIuJE4J+AlcATgecB+wIXRsQu3Yyth+0L3JqZD7S7YJlFigWQNH9ZoEnzSEQ8Afgg8PbM/EFmbs3MW4HXUisyXl+0OyMiPlK33FBE3FE8/xqwD/CvEbElIt5dd0RkRURsjog7i0KQ2fQ3Rdy7R8QFETEaEfdGxL9HxE7FvEeOvkXEoog4szjVeH0R2x11/dwaEe+KiGsi4ncRcXZE7FbMe1LxHr8tlr8gIvZqIadvBr4IHFHE/8Fi+t9GxM1FvN+JiKV1y2RE/H1E3ATcNEWfny36Gn88HBEfKOYtjYhzizhviYh31C33gYg4JyK+HhH3A8dHxK4R8aliXDYXz3edZl2Oj4jLJsR5QkTcVOTkcxERdW1/HBGfKXJ5Q0S8eEKuj6x7/YGI+Hrx8tLi72ixfkc0y7OkRhZo0vzyZ8BuwJr6iZm5Bfg+tdN0M8rM44Dbgb/KzMdl5j/XzX4hsD/wEuA9rZy2bNLfuBOBO4AlwNOA9wJT/Q7dKcB+wB8V6/L6Kdq8FjgKeDpwEHB8MX0n4CvUCtV9gDHgsy3E/yXgBGB9Ef8pEfEi4OPFe+0J3AZ8a8Kiy4HnAgdM0efbir4eR+306X3A+UVR+q/AL4AB4MXAOyNiWd3irwTOARYD3wBOpnaU9BDgYOBw4H3N1qvOy4HnFMu+Fqh/r+cCvwJ2p5b7NRHx5Bb6/C/F38XFeq5vIx5JWKBJ883uwN2Z+fAU8+4s5u+ID2bmA5m5kVqxc+wO9jduK7VCZ9/iqN+/59Q/FPxa4GOZeV9m3gF8eoo2n87MzZl5L7Vi5xCAzLwnM8/NzAcz8/fAR4E/n2W8rwO+nJlXZeZDwCpqR9j2q2vz8cy8NzPHpuskIpYAa6kd8dxArVBakpkfysw/ZOavgH8BjqlbbH1mrs3M7UXfrwM+lJl3ZeZvqR1BPa6NdTk1M0cz83bgRxT5KtwFfKoYk7OBEeAv2+hb0ixZoEnzy93A7tNcm7RnMX9H/Lru+W3A0ukatmk1cDPww4j4VUS8Z5p2SyfE8Osp2vym7vmDwOMAIuIxEfGFiLitOD14KbA4IhbMIt6l1NYfeOQI5T3UjnrNFNsjImIhtSNh38zM8aNv+wJLi1O9oxExSu1o4tNm6LchFtoflynzVdg0oVDu5JhLmoEFmjS/rAceAl5dPzEiHgu8FLi4mPQA8Ji6JntM6Geqo1cAe9c93wfYvIP91WZm/j4zT8zMPwL+Cvgf9dc71bkTqL9ubO8p2kznRGAQeG5mPoFHT8NFG32M20ytmKp1UMvvU4BNdW1mXGfgM8DvaTwd+WvglsxcXPd4fGa+bIZ+G2KhcVx21MD4NWlT9D3TmDdbd0lNWKBJ80hm/o7aKa7PRMRREbGwOO32f6hd4/W1ounVwMsi4skRsQfwzgld/Qe167wmen9xJOqPgTcCZ+9gfwBExMsj4plFMXA/sK14TPRtYFVxwf8A8Lbp+pzC46lddzZaXEd1ShvLTvRN4I0RcUhxQf7HgCuKGzKaioi3Uju9+teZub1u1k+B+yPipOKGiAURcWBEPGeG7s4C3hcRSyJid+Afga/P0L4dTwXeUWxHrwGeDXyvmHc1cEwx7zDg6LrlfgtsZ4YxlzQzCzRpnikuwn8v8Alqxc4V1I7MvLi4XgpqhdovgFuBH/JooTXu49T+0R+NiHfVTf83aqciLwY+kZnjX3o62/7G7Q9cBGyhdhTw89N899mHqBWatxTtz6F2xLAVnwIWUTvNeznwgxaXmyQzLwbeD5xL7ajeM2i8TqyZY6kVL5vr7uR8b2Zuo3YE8RBq63g3tTtInzhDXx8Bfg5cA2wEriqmdcIV1MbmbmrX7B2dmfcU895Pbb3vo/afgm+OL5SZDxbtf1yM+fM6FI/UN2Lq63Al6VHFUbhbgIXT3IDQFRHx34FjMnO2F/trGhFxPPCWzHxBt2OR+pFH0CT1jIjYMyKeHxE7RcQgtevKzut2XJLUaZX7FuqofankpcCu1OI7JzN35FoRSfPHLsAXqH3H2Si17x77fFcjkqQSVO4UZ3GR8GMzc0txG/plwD9k5uVdDk2SJGlOVO4IWvGdO1uKlwuLR7WqSEmSpBJV8hq04tbyq6l9i/WFmXlFt2OSJEmaK5U7ggZQ3Gp+SEQsBs6LiAMz89rx+RGxAlgBsNtuu/3pPvvs06VI1Snbt29np50q+f8FtcFxnB8cx/nBcaymG2+88e7MXNKsXeWuQZsoIk4BHsjMT0w1f3BwMEdGRuY4KnXa8PAwQ0ND3Q5DO8hxnB8cx/nBcaymiLgyMw9r1q5ypXXxbdiLi+eLgCOBG7oblSRJ0typ4inOPYEzix8w3gn4dmZe0OWYJEmS5kzlCrTMvAY4tNtxSJIkdUvlTnFKkiT1Ows0SZKkirFAkyRJqhgLNEmSpIqxQJMkSaoYCzRJkqSKsUCTJEmqGAs0SZKkirFAkyRJqhgLNEmSpIqxQJMkSaoYCzRJkqSKsUCTJEmqGAs0SZKkirFAkyRJqhgLNEmSpIqxQJMkSaoYCzRJkqSKqVSBFhF7R8SPIuL6iLguIv6h2zFJkiTNtZ27HcAEDwMnZuZVEfF44MqIuDAzf7kjna7dsInV60bYPDrG0sWLWLlskOWHDnQm4h5lThqZj0bmYzJz0sh8NDIfk5mTRuP52GWPZ/5pK+0rVaBl5p3AncXz30fE9cAAMOsCbe2GTaxas5GxrdsA2DQ6xqo1GwH6dkMxJ43MRyPzMZk5aWQ+GpmPycxJo4n5aEWlTnHWi4j9gEOBK3akn9XrRiYlZGzrNlavG9mRbnuaOWlkPhqZj8nMSSPz0ch8TGZOGk2Vj2YqdQRtXEQ8DjgXeGdm3j/F/BXACoAlS5YwPDw8bV+bRsemnT7TcvNZFXOyZcuWrr13FfPRTTuSj26OY5n6bRtpNo79lo9mqpoPP1erY7p8zKRyBVpELKRWnH0jM9dM1SYzTwdOBxgcHMyhoaFp+xu4/JIpEzOweBEzLTefVTEnw8PDXXvvKuajm3YkH90cxzL12zbSbBz7LR/NVDUffq5Wx3T5mEmlTnFGRABfAq7PzP/ZiT5XLhtk0cIFDdMWLVzAymWDnei+J5mTRuajkfmYzJw0Mh+NzMdk5qTRVPlopmpH0J4PHAdsjIiri2nvzczvzbbD8YsR333ONfxh23YGvJPEnExgPhqZj8nMSSPz0ch8TGZOGtXno1WVKtAy8zIgOt3v8kMHOOuntwNw9luP6HT3PcmcNDIfjczHZOakkfloZD4mMyeNxvNx029uvrKV9pU6xSlJkiQLNEmSpMqxQJMkSaoYCzRJkqSKsUCTJEmqGAs0SZKkirFAkyRJqhgLNEmSpIqxQJMkSaoYCzRJkqSKsUCTJEmqGAs0SZKkirFAkyRJqhgLNEmSpIqxQJMkSaoYCzRJkqSKsUCTJEmqGAs0SZKkirFAkyRJqpidux3ARBHxZeDlwF2ZeWCz9rfev53nn3oJK5cNsvzQgY7EsHbDJlavG2Hz6BhLFy/qaN9lKTPmXswHlBe3+Zi+702jYwxc3tn9sSzuM5OVNY69ng+3kRrzMVmZcVeuQAPOAD4LfLXVBTaNjrFqzUaAHU7M2g2bWLVmI2Nbt3W877KUGXMv5gPKi9t8zG3fZTEfk7nPNHIbaWQ+Jis77sqd4szMS4F7211ubOs2Vq8b2eH3X71u5JFkd7rvspQZcy/mA8qL23zMbd9lMR+Tuc80chtpZD4mKzvuKh5BayoiVgArAHbZ45mPTN80Osbw8PCUy4yOjgFMO7++j+mmN1u2W2Ybcys5mat8bNmypaP9zSbuKuWj08rKx2z77rb5sM90mvtMo/mwz3Tyc9V9ZrId2UZa0ZMFWmaeDpwOsOue++f49IHFixgaGppymdNG1gMwNHTEjH0PXH7JlEmfqe9um23MreRkrvIxPDzc0f5mE3eV8tFpZeVjtn1323zYZzrNfabRfNhnOvm56j4z2Y5sI62o3CnO2Vq0cAErlw3ucD8rlw2yaOGCUvouS5kx92I+oLy4zcfc9l0W8zGZ+0wjt5FG5mOysuPuySNoEw108M6J8T7efc41/GHb9o72XZYyY+7FfEB5cZuPue27LOZjMveZRm4jjczHZGXHXbkCLSLOAoaA3SPiDuCUzPzSdO13WwA/fs+LOhrD8kMHOOuntwNw9ltnPlRdFWXG3Iv5gPLiNh/T9z06Osq6kzq7P5bFfWayssax1/MBbiNgPqZSZtyVK9Ay89huxyBJktRN8+YaNEmSpPnCAk2SJKliLNAkSZIqxgJNkiSpYizQJEmSKqbpXZwRsRR4MfAY4IbM/LfSo5IkSepjMx5Bi4iXAFcCLwOOAD4VESMR8WdzEZwkSVI/anYE7SPAf87Mm8cnRMQRwL9ExJuBBzLz2jIDlCRJ6jfNCrRd6oszgMxcHxGvBi4AHgL+U1nBSZIk9aNmNwn8v4hYMnFiZt4IbKN2bZokSZI6qFmBthpYW9wo8IiI2B14KDPvKi0ySZKkPjXjKc7MPDcidgXWR8SVwC+AXYDXUrs+TZIkSR3W7C7OI4CzgGdTu+bsicBW4K8z88zyw5MkSeo/zW4SeAPwOeBG4AfAP2fmb0qPSpIkqY81O8V5AkBEPAt4KXBGRDwR+BG1gu3Hmbmt9CglSZL6SEs/9ZSZN2TmJzPzKOBFwGXAa4ArygxOkiSpH814BC0injzNrMuLR0TE4swc7XhkkiRJfarZNWibi0dMmJ7F3wAWAPt0OC5JkqS+1axAuz4zD52pQURs6GA8RMRRwP+iVvh9MTNP7WT/kiRJVdfsGrQjWuijlTYtiYgF1O4afSlwAHBsRBzQqf4lSZJ6QWRm80YRH8nM902YtqDTd3AW37v2gcxcVrxeBZCZH59umf0fvzgvWv6Kpn3/8s77AThgzye0FEs77W+95wEA9nvKY1vqu5327bQtcx3b7bvddXzooYcYXDrdJY+z6xfKW8ey89Fq23bbl7mOv7zzfh5++GEO2rvz4+g+Mzdtob1x7MV8tNu+7H2m1fazWccyPlerss9U5TME2l/Hl118/pWZeVizts1OcY4biIhjM/MsgIh4KnA28MIWl2/VAPDrutd3AM+d2CgiVgArAJ7x2MczOtr8HoUFuR2gpbbttr//waLtgq0t9d1O+3balrmO7fbd7jpmZim5Lmsdy85Hq23bbV/mOi7I7ey0Uznj6D4zN22hvXHsxXy0277sfabV9rNZxzI+V6uyz1TlMwRmt46taPUI2i7AOuAkajcIfAU4KTO/2/I7tRJMxGuAZZn5luL1ccDhmfn26ZYZHBzMkZGRTobRtv/2hfUAnP3W1s72ttO+3b6rot11HB0dZd1JL+1ov1VS5phXKSfDw8MMDQ01bec+M1lZ6zibfLQ6jmXql32mVbNZx/n8udrLnyERseNH0CLiq8BVwAbg74FvAg8DyzPz5k4EOsEdwN51r/eidhepJElS32h2k8CZRZs3USvO9gPuA14fEUeXEM/PgP0j4unFUbtjgO+U8D6SJEmV1eynni4GLh5/HRE7U7u78mDgecA5nQwmMx+OiLdRO526APhyZl7XyfeQJEmqulZvEgBqBRRwTfH4WhkBZeb3gO+V0bckSVIvmPEUZ0Rc1ayDVtpIkiSpdc2OoD07Iq6ZYX4AT+xgPJIkSX2vWYH2rCmm7QpsB8a/IKSjX1YrSZLU75rdJHBb/euI+DDwFGpHzvYA/jYz7y4vPEmSpP7T1k0CwOLM/DuAiNgPWB0Rv8vMd3Y6MEmSpH7V7HvQJnrkdGZm3krtVOcZHYxHkiSp77VboF0WER+PiKURsRTYLTOvLiMwSZKkftXu96CdExGbgQ8Xy360lKgkSZL6WEsFWkRcBJyYmb/IzJ8APyk3LEmSpP7V6inOdwOfjIivRMSeZQYkSZLU71oq0DLzqsx8EXAB8IOIOCUiFpUbmiRJUn9q+SaBiAhgBDgNeDtwU0QcV1ZgkiRJ/aqlAi0iLgM2AZ8EBoDjgSHg8Ig4vazgJEmS+lGrd3GeAFyXmTlh+tsj4voOxyRJktTXWirQMvPaGWb/ZYdikSRJEu1/Ue0kmfmrTgQiSZKkmh0u0CRJktRZFmiSJEkVU5kCLSJeExHXRcT2iDis2/G0au2GTWy4fZQrbrmX5596CWs3bOqJvstUVtzmY277Lov5mMx9ppHbSCPz0ahf8lGZAg24Fng1cGm3A2nV2g2bWLVmI3/Yth2ATaNjrFqzsSMDWmbfZSorbvMxt32XxXxM5j7TyG2kkflo1E/5qEyBlpnXZ+ZIt+Nox+p1I4xt3dYwbWzrNlav2/HVKLPvMpUVt/mY277LYj4mc59p5DbSyHw06qd8tPo9aJUSESuAFQBLlixheHi4K3FsGh2bdvpMMY0Wy83UZrZ9d9ts4h4dHWPbtm3mo9DK9jHbvsu2ZcuWUsbRfaZR2floNo5l6rd9ppkd2Wfm4+dqP32GzGmBFhEXAXtMMevkzDy/1X4y83TgdIDBwcEcGhrqTIBtGrj8kikHdGDxImaK6bSR9QAMDR3R8b67bTZxnzayntHR0RnXq9/yATNvH7Ptu2zDw8OljKP7TKOy89FsHMvUb/tMMzuyz8zHz9V++gyZ01OcmXlkZh44xaPl4qxKVi4bZNHCBQ3TFi1cwMplg5Xuu0xlxW0+5rbvspiPydxnGrmNNDIfjfopHz15irMqlh86ANTOW28eHWPp4kWsXDb4yPSq9l2msuI2H3Pbd1nMx2TuM43cRhqZj0b9lI/KFGgR8SrgM8AS4LsRcXVmLutyWE0tP3SgtMErs+8ylRW3+ZjbvstiPiZzn2nkNtLIfDTql3xUpkDLzPOA87odhyRJUrdV5ms2JEmSVGOBJkmSVDEWaJIkSRVjgSZJklQxFmiSJEkVY4EmSZJUMRZokiRJFWOBJkmSVDEWaJIkSRVjgSZJklQxFmiSJEkVY4EmSZJUMRZokiRJFWOBJkmSVDEWaJIkSRVjgSZJklQxFmiSJEkVY4EmSZJUMZUp0CJidUTcEBHXRMR5EbG42zFJkiR1Q2UKNOBC4MDMPAi4EVjV5XhKsXbDJjbcPsoVt9zL80+9hLUbNnU7pK4az8fIfdvNB24fUzEnjcxHI/MxmZ+rjXp1G6lMgZaZP8zMh4uXlwN7dTOeMqzdsIlVazbyh23bAdg0OsaqNRt7ZmPpNPPRyHxMZk4amY9G5mMyc9Kol/NRmQJtgjcB3+92EJ22et0IY1u3NUwb27qN1etGuhRRd5mPRuZjMnPSyHw0Mh+TmZNGvZyPnefyzSLiImCPKWadnJnnF21OBh4GvjFDPyuAFQBLlixheHi488GWYNPo2LTTe2UdOsl8NJoP+diyZUtHY50POemkucpHp8exLG4fk5mTRr2cjzkt0DLzyJnmR8QbgJcDL87MnKGf04HTAQYHB3NoaKiTYZZm4PJLptxYBhYvolfWoZPMR6P5kI/h4eGOxjofctJJc5WPTo9jWdw+JjMnjXo5H5U5xRkRRwEnAa/IzAe7HU8ZVi4bZNHCBQ3TFi1cwMplg12KqLvMRyPzMZk5aWQ+GpmPycxJo17Ox5weQWvis8CuwIURAXB5Zp7Q3ZA6a/mhA0DtnPjm0TGWLl7EymWDj0zvN/X52DQ6xoD5ANw+6pmTRuajkfmYzM/VRr28jcQMZxJ7wuDgYI6MVP9iP82sV06paGaO4/zgOM4PjmM1RcSVmXlYs3aVOcUpSZKkGgs0SZKkirFAkyRJqhgLNEmSpIqxQJMkSaoYCzRJkqSKsUCTJEmqGAs0SZKkirFAkyRJqhgLNEmSpIqxQJMkSaoYCzRJkqSKsUCTJEmqGAs0SZKkirFAkyRJqhgLNEmSpIqxQJMkSaoYCzRJkqSKqVSBFhEfjohrIuLqiPhhRCztdkySJElzrVIFGrA6Mw/KzEOAC4B/7HZAkiRJc61SBVpm3l/38rFAdisWSZKkbtm52wFMFBEfBf4G+B3wwi6HI0mSNOcic24PUkXERcAeU8w6OTPPr2u3CtgtM0+Zoo8VwIri5YHAtWXEqjm1O3B3t4PQDnMc5wfHcX5wHKtp38xc0qzRnBdorYqIfYHvZuaBTdr9PDMPm6OwVBLHcX5wHOcHx3F+cBx7W6WuQYuI/etevgK4oVuxSJIkdUvVrkE7NSIGge3AbcAJXY5HkiRpzlWqQMvM/zqLxU7veCDqBsdxfnAc5wfHcX5wHHtYZa9BkyRJ6leVugZNkiRJFS3QIuLLEXFXRFxbN+3JEXFhRNxU/H1SMT0i4tMRcXPxM1F/UrfMG4r2N0XEG7qxLv2qzTF8VkSsj4iHIuJdE/o5KiJGivF9z1yvR79rcxxfV+yD10TETyLi4LplHMcuanMcX1n3k3s/j4gX1C3jZ2oXtTOOdfOfExHbIuLoummOYw+oZIEGnAEcNWHae4CLM3N/4OLiNcBLgf2LxwrgNKhttMApwHOBw4FTJm64KtUZtD6G9wLvAD5R3zgiFgCfozbGBwDHRsQBJcasyc6g9XG8BfjzzDwI+DDF9S+OYyWcQevjeDFwcPGTe28Cvgh+plbEGbQ+juP73j8B6+qmOY49opIFWmZeSu0f7XqvBM4snp8JLK+b/tWsuRxYHBF7AsuACzPz3sy8D7iQyRu2StLOGGbmXZn5M2DrhPaHAzdn5q8y8w/At4o+NEfaHMefFPsawOXAXsVzx7HL2hzHLfnoxcn1P7nnZ2qXtflvI8DbgXOBu+qmOY49opIF2jSelpl3AhR/n1pMHwB+XdfujmLadNPVPdON4XQcw2pqZRzfDHy/eO44VtO04xgRr4qIG4DvUjuKBo5jVU05jhExALwK+N8T2juOPaKXCrTpxBTTcobp6h2OYQ+KiBdSK9BOGp80RTPHscIy87zMfBa1ozEfLiY7jr3lU8BJmbltwnTHsUf0UoH2H8WpS4q/44ds7wD2rmu3F7B5hunqnunGcDqOYTVNO44RcRC1a5ZemZn3FJMdx2pquj8Wp9SeERG74zhW1XTjeBjwrYi4FTga+HxELMdx7Bm9VKB9Bxi/2+QNwPl10/+muJvzecDvisO864CXRMSTigsgX0LdhZLqiunGcDo/A/aPiKdHxC7AMUUf6q4pxzEi9gHWAMdl5o117R3HappuHJ8ZEVE8/xNgF+Ae/EytqinHMTOfnpn7ZeZ+wDnA32XmWhzH3pGZlXsAZwF3Urto/A5qp0ueQu0OlZuKv08u2ga1O8T+L7AROKyunzcBNxePN3Z7vfrp0eYY7lG0uR8YLZ4/oZj3MuDGYnxP7vZ69dujzXH8InAfcHXx+HldP45j74zjScB1xRiuB15Q14+fqT0yjhOWOwM42nHsrYe/JCBJklQxvXSKU5IkqS9YoEmSJFWMBZokSVLFWKBJkiRVjAWaJElSxVigSZIkVYwFmiRJUsVYoElSnYg4MiK+1u04JPU3CzRJanQwsKHbQUjqbxZoktToYGBDROwaEWdExMfGf5tSkubKzt0OQJIq5mDgLmo/IP3FzPx6l+OR1If8LU5JKkTEQuBu4DbgrZm5vsshSepTnuKUpEcdAPwMeBjY1uVYJPUxCzRJetTBwE+AY4CvRMTTuhyPpD5lgSZJjzoYuDYzbwROAr5dnPaUpDnlNWiSJEkV4xE0SZKkirFAkyRJqhgLNEmSpIqxQJMkSaoYCzRJkqSKsUCTJEmqGAs0SZKkirFAkyRJqpj/DzjsDdhUDXVAAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 720x216 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# compute input signal\n",
    "x = np.random.uniform(low=-1, high=1, size=256)\n",
    "x = np.concatenate((x, np.zeros(1024)))\n",
    "# compute output signal\n",
    "yQ = sos_df1(x, a, requantize=uniform_midtread_quantizer)\n",
    "\n",
    "# plot results\n",
    "np.seterr(divide='ignore')\n",
    "plt.figure(figsize=(10, 3))\n",
    "plt.plot(20*np.log10(np.abs(yQ)))\n",
    "plt.title('Level of output signal')\n",
    "plt.xlabel(r'$k$')\n",
    "plt.ylabel(r'$|y_Q[k]|$ in dB')\n",
    "plt.grid()\n",
    "\n",
    "plt.figure(figsize=(10, 3))\n",
    "k = np.arange(1000, 1050)\n",
    "plt.stem(k, yQ[k]/Q)\n",
    "plt.title('Output signal for zero input')\n",
    "plt.xlabel(r'$k$')\n",
    "plt.ylabel(r'$y_Q[k] / Q$ ')\n",
    "plt.axis([k[0], k[-1], -3, 3])\n",
    "plt.grid();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise**\n",
    "\n",
    "* Estimate the period of the small limit cycles. How is it related to the poles of the system?\n",
    "* What amplitude range is spanned?\n",
    "\n",
    "Solution: The period of the small limit cycles can be estimated from the second illustration as $P = 6$. The normalized frequency of a harmonic exponential signal with the same periodicity is given as $\\Omega_0 = \\frac{2 \\pi}{P} = \\frac{\\pi}{3}$. The poles of the system can be extracted from the code of the first example as $z_{\\infty 0,1} = 0.9 \\cdot e^{\\pm j \\frac{\\pi}{3}}$. The periodicity of the small limit cycles is hence linked to the normalized frequency of the poles. The amplitude range spanned by the small limit cycles is $\\pm 2 Q$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Large Limit Cycles\n",
    "\n",
    "Large limit cycles are periodic oscillations of a recursive filter due to overflows in the multiplications/additions. As for small limit cycles, large limit cycles may be present even after the input signal is zero. Their level is typically in the range of the minimum/maximum value of the requantizer. Large limit cycles should therefore be avoided in a practical implementation. The presence of large limit cycles depends on the scaling of input signal and coefficients, as well as the strategy used to cope for clipping. Amongst others, they can be avoided by proper scaling of the coefficients to prevent overflow. Again, a detailed treatment of this phenomenon is beyond the scope of this notebook and can be found in the literature."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Example - Large limit cycles of a SOS\n",
    "\n",
    "The following example illustrates large limit cycles for the system investigated in the first example. In order to trigger large limit cycles, the coefficients of the filter have been doubled. The input signal is uniformly distributed white noise till time-index $k=256$ and zero for the remainder."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAm4AAADiCAYAAAD3RUuZAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3XmYXGWZ9/Hv3Wv2PemsZCEbBAhgGyAgdALKqqjji/COC7jEdZwZcRBkHJ1R0FFfGfcxKouoIKKMDCCy2RAIIYQlISxJOiEhnYQsZO0snV7u94861V1VXVV9qrq6q6r797muurrOOc8556mnT7ruPKu5OyIiIiJS+ErynQERERERCUeBm4iIiEiRUOAmIiIiUiQUuImIiIgUCQVuIiIiIkVCgZuIiIhIkVDgJiIFzcyuNLMnszx3lpm9YGYHzOwLuc5boTCzl82spgfuU2tmn+ju+4hIagrcRCQjZrbRzM7Ldz5CugaodffB7v7DXF441+XQleu5+xx3r81VXkSkcClwE5HebDLwcr4zISKSKwrcRCRnzOwSM3vRzPaa2VIzOynYf62Z3Z2Q9gdm9sPg/VAz+5WZbTOzLWb2TTMrDXnP9wRNhXuDprzjgv2PAQuAH5tZg5nNTHLueDO718x2m1mdmX0y5titZvbNmO0aM6sP3t8OHAP8b3Dta8xsipm5mS0ys63BZ7k62+slyesoM7sv+Jy7zWyJmZUEx9pq68ysv5ndZmZ7zOzVIG/1MdfZaGZfMrNVZrbPzH5vZv2CY8ODe+wMzr/PzCaG+T2ISM9Q4CYiOWFmpwI3A58CRgI/B+41s0rgDuAiMxsSpC0FLgN+F5x+G9AMTAdOAd4FdNqXKgjG7gD+CRgNPEAk+Klw94XAEuDz7j7I3dcmucQdQD0wHvgAcKOZndvZfd39w8AbwLuDa38n5vACYEbwGa4N0/zZyfWirg7yOhqoAr4CJFuz8GvAFGAa8E7gQ0nSXAZcAEwFTgKuDPaXALcQqak8BjgM/Liz/ItIz1HgJiK58kng5+7+jLu3uPttQCNwurtvAp4H3hukXQgccvdlZlYFXAj8k7sfdPcdwE3A5SHu+UHgfnd/2N2bgO8B/YH5nZ1oZpOAs4Avu/sRd38R+CXw4Uw+dBL/HnyOl4gEQVd08XpRTcA4YLK7N7n7Ek++2PRlwI3uvsfd64Fkfft+6O5b3X038L/AyQDu/pa7/9HdD7n7AeAG4Jwc5V9EckCBm4jkymTg6qApb6+Z7QUmEanNgkjtWjSI+b+017ZNBsqBbTHn/RwYE+Ke44FN0Q13bwU2AxNCnrs7CFCiNoU8N53NCdcbnyphhr4L1AEPmdkGM7s2RbrxCXnYnCTNmzHvDwGDAMxsgJn93Mw2mdl+4AlgWNhmaxHpfgrcRCRXNgM3uPuwmNcAd78jOP4HoCboM/U+2gO3zURq5kbFnDfE3eeEuOdWIoEfAGZmRILFLSHPHWFmg2P2HRNz7kFgQMyxsQnnJ6vtIrh/7PW2dvF6kYPuB9z9anefBrwb+GKKZt1tQGy/tElJ0qRyNTALOM3dhwBnB/stg2uISDdS4CYi2Sg3s34xrzLgF8Cnzew0ixhoZhdHAyN33wnUEmk+fN3dXw32bwMeAv6fmQ0xsxIzO9bMwjTR3QVcbGbnmlk5kcCjEVja2YnuvjlI963gM5wEfBz4bZDkRSL98kaY2Vgi/ehibSfSjyzRV4OaqznAVcDvu3g9oG3gx/QgON0PtASvRHcB1wUDDSYAn091zSQGE+nXttfMRhDpLyciBUSBm4hk4wEiX/DR19fdfQWRfm4/BvYQada7MuG83wHn0V7bFvURoAJ4JTj3biL9udJy9zVEOt//CNhFpCbq3e5+NOTnuIJIR/6twD3A19z94eDY7cBKYCORwPL3Ced+C/jXoHn3SzH7Hyfy2R8FvufuD3XxelEzgEeABuBp4Kcp5m77DyKDGF4P0t9NJJgN47+I9BHcBSwDHgx5noj0EEvet1VERDJhZlOIBEvl7t6c39y0M7PPAJe7uwYZiPQCqnETEelFzGycmZ0ZNDnPItJ8fE++8yUiuVGW7wyIiEhOVRAZlTsV2AvcCfw0rzkSkZxRU6mIiIhIkVBTqYiIiEiRUOAmIiIiUiR6bR+3UaNG+ZQpU7r1HgcPHmTgwIHdeo/eSmXXNSq/7Knssqey6xqVX/b6Qtk999xzu9x9dGfpem3gNmXKFFasWNGt96itraWmpqZb79Fbqey6RuWXPZVd9lR2XaPyy15fKDsz29R5qgJrKjWzm81sh5mtjtk3wsweNrN1wc/h+cyjiIiISL4UVOAG3ApckLDvWuBRd59BZCbyVAsri4iIiPRqBRW4ufsTwO6E3ZcCtwXvbwPe26OZEhERESkQBRW4pVAVLEIdXYx6TJ7zIyIiIpIXBTcBb7De333ufkKwvdfdh8Uc3+PuSfu5mdkiYBFAVVXV2+68885uzWtDQwODBg3q1nv0lA17W9h12Jk3rmfGqzQ0NPDS/n6cOKqUQRWWNu3hZqd/Wfo0fU1vevZ6msoueyq7rlH5Za8vlN2CBQuec/fqztIVw6jS7WY2zt23mdk4YEeqhO6+GFgMUF1d7d09AqW7Rrnc+tTrPPzqdn77idOzvsaDq7cxfcxgpo8J96Bfee39AFxzxXlZ3zMTf3jgMX6+6jBnTh+Z9nPW7WjgvO8/znf+7iQue/skGhqbaTjSzNih/Xokn4WqL4yw6i4qu+yp7LpG5Zc9lV27YmgqvRf4aPD+o8Cf85iXHvH1/32Fp+reSnl89ZZ9TLn2fl7dtp+jza20tHasNf30b57nvO8/3p3Z7JKm1sjPbfuOpE23bvsBAB59bTsAl/74SU7/1qPdmjcREZFCVVCBm5ndATwNzDKzejP7OPBt4J1mtg54Z7BdcNyd2jU7aE0SRGUiGqik85fV2wB45JXtzPzXv3DF4mUp07a0OkeaWrqUp0SHjjZz6U+eYvWWfW37PnLzcj5xW/h586KlVGIdm0DX72xgbVAOienW7zyYVZ5FRER6g4IK3Nz9Cncf5+7l7j7R3X/l7m+5+7nuPiP4mTjqtCDct2obV97yLLcvSz9/3rZ9h3njrUNt24+v3cnfXmtv/X3nTU90eq+mlkg4U1ISCWaWb0xdJP945wvM/uqDnV4zrD0Hj3Lpj59i5ea9fPsvr7Xtf2LtTh55dXvo60S7VpYk6bp27v97nHcF5dAaJEwW4ImIiPQ1BRW4FbOtew8DUL/nUNp0Z3zrMc7+7t/atj9683KuuvXZ0Pc5fLSFxU9sAKA0WdRDpPYv6r5V29Jeb/fBo9z85Otx+6794yp+uWRD0vSX/OhJ1u1oaLv/jv1HCDPAJbHWL3qGkT4gi1ZgvrH7EFOCfngiIiJ9VTEMTihouxoaGdKvvD0QiakZam5pZc+hJkYPruz0Op//3fO81XA0bt/T699iZtUgRg6q5MHV2xgxsJIpIwe0HS9NUQuVrM9bKgu+V8u+w01x++58djMAn3jHtLj9dTsa2BIEqAAb3zrIvBsf5d8uOT7tPVZs3M0H/vtpfv2xeZw9czSr6vfy1aci1+msIi0aFL4U0ywrIiLSV6nGrQvcnepvPsIX73qxrekvNg751/9ZzdtveCRUH7P7Vm3j6Q3xAxKu+MUy/u8vnuGuZzfz6d88z2U/f5rYkKwkRY1bcwaBW2LQls7+I/Fp39gdqV1cun5X275jv/JAh/OeeT3SlPv0hrfYceAIi379XNuxxCbQbfsOx22nqsxraGwOnW8REZHeQoFbF0Tjo/tWbcPpGLk98FKkmfKDi5dlFCDFWrP9ANf8cVXbdmwgU5qitqqppTWre0HygOhIUwubd3dsAo7ePjb4aml1frlkA2vePBC3L5Jf46pbnuXN/e0jSRNr3N7z46fa3l95y/K2Pm6JTvjaX9l5oLHTzyMiItKbKHDrgtjwqL3GrT0SiTabrty8l3tXbs3JPf/nxS1t71P1cWtu6RjshJ1oeeH3auO263YcYPZXH+Qd3/lbytqvsoQI8pv3v8olP1rStt0+wIC4plaID9weX7szLhirXbOzvVyTfNTn39gTuX6r88Xfv9i2LSIi0lspcOuCZC2S0QDjwdXb4mrZhvYvz8k9Y0dyJjaVtrQ6j6/dSVNrxxq3sAtk7IgJnHbsP8Ivl8QOXIi/SDQwTTbisykmeIxOkVJSYh2mSykx4819R3jjrUN89OblHa6TblTpp26PNLnuO9zEn17Ywvt/upQXFLyJiEgvpsEJXRAbg0QDkmh48enfPB+XNleBW6zEwQn//fh6vvvXNdz4vhM7pG11p6STEZyJ5t0YP9Ftqq5zqWr+Es8rNesQQBqknVC3JaFcEx1sbI6rjfvbazs45ZikK6KJiIgUPdW4dUFc4NY2L1nyEKOiNPdFnVjjtnFXZHLarQnNkQAtCRHTjQ+8ytnf+VuHdOncvaI+6f50gdvBxub2WrMS69BnzToZVno06K+Xqlyvvmtl3O+hs+uJiIgUM9W4dUFsDBIdnJAqhgnbxywTD73cPuHtzgON/OG5SGCVrKn03he3cteKzRw/bghfveT4trngMvH7FZvjtqMfNdW0JC9u3st7f/IUx44eGElvHWvtOouzvvvXNWnTPfjym1z29olt25qoV0REejMFbl0QG4S0xWVm3PLU6x3S5j5sI26lgnteaK8Na2rueLd/uTsyMvXZjXt415yxObl/NEZKHJwQtap+L9C+TFWJWYeav84CrQNHOp/242O3ti+11UmrrYiISFFT4NYFsc1+0RGNq+r3UrtmZ9q0AK+9uT/p6M9sbd3bPsVGc5Iat1iZzPOWTtvSWymCr8Rmy0gft/h7P7cp3GCCVNOCJEo1t52IiEhvoMCtC2LDoyXrIpPQHmpMPtluYqx0wX8tSZouW7HBU2eBWeLIzq4qSxEsJe4uKbGMVnWI1ZTDIFdERKRYFc3gBDO7wMzWmFmdmV2b7/xA8lGWqWq7wtYYZatuZ0N7HoIO/alaIXNV4xaVqpYrsSauJEkft1xTHzcREenNiiJwM7NS4CfAhcDxwBVmln6BzB6QPHBLEZl0c8DyxNr25tnmTpows631SiXV4ITEeK6zaUNykpeieKJFRESyUyxNpfOAOnffAGBmdwKXAq/kM1OHkwwCaGzKT41brKYgMEsVoCWbLqQrSlMMTkjs49YTU3UYxvqdDSyt28WC2WN4ZsNuJo0YgLuz8a2DzD92FE+s28kZ00byyrb9DKgoZdzQ/jy7cTfnHlfFk+t2Mmf8UHYeaGT/kSZOmjiMJ+t2UTNzNCs27aZqSD8qy0pZ8+YB3jFjFEvW7WLe1OGs3d5AeWkJk0cOYNmGtzj3uCqeqtvFzKrB7D/cxK6GRt42eThPrNvFOTNG88LmPYwYWMHQ/uWsqt9HzazRLFm3i1OPGc7Gtw7i7kwfM5in1+9i4XFVPL3+LaaOGkBjcyv1ew7DoVZ+s2wTZ00fxUtb9jG4XxmjB1fy/KY9LDyuiiVrd3LixKFs23uEw00tzBk/hKeCMln++m4mDOtPSYmxfkcDZ06PlMnp00by2rYD9CsvYcLw/ix/PVImT63bxXHjhrDrYCP7DjVx8qRhLAnK5LlNexgzuJL+FaW8uu0AZ88cxRNrI2WybnsDJSXGtFEDeTook6V1u5g+ZhAHG1vYvv8I1VOGs2TdLs4OymT4gAqGD6jgxfq9LJg1mifX7eLkY4bxxluHaG51ZlZFyiT6u508cgBNLc7m3Yc449iRPLFuJ/OPHcXqLfsY1K+MqsH9eG7TbhYGv9sTJgzlxR3N7F+5lRMnDI373Y4b2p+yEmPdjgbOSiiTyrISJo0YwDOvt/9uZ48dzO6DR9lz6CinTBoeVyajB1cysLKMl7fu45yZkd9t9eThrN/ZABjTxwxk6fr2Mjl2zCAOH21h274jzJsygifW7eQdM0axsn4fQ/uXM3JgBS9s3svC2WNYsnYncycNo37PYY42tzJ73OC45/2YkQNobU3/vK/Y2F4mc8YPZceBIxw40hx53tftpGbWGFZs2s3YIf2pKCth7fYDnDV9FLWbm5i4o4G12w9k/LyPHFjJ4H5lvLSl/Xl/2+ThbNh5EHCmjxnE0+vf6vC8b9lzmNOmRn63Z00fxaot+xjSr4xRgyp54Y325/2kicPYuvdw0ud94vABAGzY2f68nzFtJK9m+bz3qyjltSTPe2mJMTXJ897Q2MyO/Y20HGrl9mWbOvwNWFm/L+553/TWIVpbnRkxz/uyDbuZkuR5P/PYyN+AQf3KGBPzNyD6vL+57wgHj7ZwQlAm0d9t9HmvS/I3IJvnfUBFKa9s29/heTfr+Dcg1fP+4ua9DB9QwYgkz/vz25s5uGpb0ue9pdXZFPO8zz92FC9v3cfAijKqhsT/DejwvMf8DYg878ba7ZG/AdFn9IU39nBZ9aSC6UNt3TFNRa6Z2QeAC9z9E8H2h4HT3P3zqc6prq72FStWpDrcZS9v3cfFP3yyw/7JIwew6a2O63r+4iPVnHfcGKZe13ER9ly76MSxPPDSm91+n6jP1hzLT2vXd9j/L+fPapvOA+DG953IV+55qcfyJSIikgt3feoM5k0d0a33MLPn3L26s3TFUuOWLMztEHGa2SJgEUBVVRW1tbXdlqEl9ckXjd/fkLw265O/XsFXTuvXbfmJtfXNjqNau9PGTW8k3R8btAGsW7smaToREZFCNWlwCYc2raJ2U75zElEsgVs9MClmeyLQYdV2d18MLIZIjVtNTU23Zeit5+ph9coO+720DEge1N3zRgVwJOmxXCobMATouTU7J06aBK93PqHvzFmz4OXOa9xKuzD6VEREJJdGDhtCTc2Z+c5Gm2Lpyv0sMMPMpppZBXA5cG8+M5Sqo32qPm4Ag/v1TJy8bV/3B4exwgZZYfv5KWgTERFJrihq3Ny92cw+D/wVKAVudveX85mnVP3sG5uTz+MG3bPQfDJbcjz4oDPhA7duzoiIiEgvVxSBG4C7PwB0f8/+kFJNtZEuOOmpwK2nha1JK4aBMCIiIrEK7ZurWJpKe4V+5b2zuJONok2mqys2PH3dQv7l/FlduoaIiEgx652RRA/obD3QZFpCnpKs/9yMMYMyvl9PqCgr4fG14UaxdrWpdNzQ/lQN6ZmRuSIiIpB8Wot8UuCWpaYkk+92ek7IyK0nJ+vtqsoMliqI/Vwp5uztVIHMfygiIpIXCtyydDRs9VnsOc3hzuksblswa3Sn1ygPIqNUC8DnSkVZ+EcobDx63nFjUh7TWqQiItKXKXDLUtjas1jpRpxmYnC/zgc59CsrBbp/fdBMlgCJrXFLd9bJk4ZR+6WapMcUt4mISF+mwC1L2QRuYWvckvEU71OpLI8EbhUFtOp6XB+3NAFYujVNVeMmIiI9qdA6LxXOt3qRaWrJ/FfZGDJwm3/syA77wkylMXpwZdv76AjWsqDJNFcBXGINXiZTfIStcTNLXbPW3TWIIiIihUyBW5ayqT0LW0s3YmAFs8cOTnl86sgBSff3D2rZAPoF78uCgK0yg75o6SQGTs0ZDBX1kIFbiRmWIoXiNhER6csUuGWpp5tKo/7tkuP5P9WTkh6LrdGK1riVB5FOZcIccmNiaucAzpnZ+YAH6DjYoSWDmsewTaUlaWrc0jWjioiI5FqhfesocMtSdoMTUp8zcmBF2nOjMc8ZSZpR29LEBEZtgxNSNJUmhlvvP3UCY0PMkdaVGrew05ykqm0D9XETEZG+TYFblrLp4xatcfvcgmM7HPvwGZPTnxzcrrTEUtZGHW5qH7UabSotL4n8ihOn7YiuL3rhCWNZeu1CLj15Qqf5BygvTX6dMGKTZt/HLfTtREREeh19DWYpm3ncojVuwwe01671Ly9l47cv5uRJw0Jdo8QsZXPhnkNH294nDk4oS4h4orVfJSXG+GH9Q36CSOC48dsXt22HXUGixML3cbM0n1FNpSIi0pM0qjQFM/s/ZvaymbWaWXXCsevMrM7M1pjZ+fnKY6wb3nsCwyozCyK27D0MtNdazT92JK9+4wIgvgk03UNiljroiW2JrAyaSsuCGrfShIAnWlOWadNjYh+3sBVuJWahm0rTDUBQ2CYiIn1ZwQRuwGrg/cATsTvN7HjgcmAOcAHwUzMr7Xh6zzIz+sXkYt6UEfzg8pNDnRsN3BIDrceuPgeA95+SvtkyNtb66d+fyoYbL+qQJjoYIbqCQmLftGjglunSU9lMx/H9y+YGgVv7vnTxYmRUaXKpgr8LTxibcb5ERESKTcEEbu7+qruvSXLoUuBOd29099eBOmBez+YuuYHl7eHFzLGDGBJiRQNoD6Y8oW5t2uhBbPz2xZx7XFWHc+IHZLbf96ITx7WtXhC7FFbidCBlpckHFXS1s3+YQO6cmaMxg9aQ1XPp+rilaqH+2YfeFuraIiIimSi0lp6yfGcghAnAspjt+mBfB2a2CFgEUFVVRW1tbbdm7MqZLdz4gnG4GbZu2corzTtCnbdubSQ+3bNnb8o8Hjx4OG770KFDADy7fDmVMTV90fMXv3MAZSUH+VsQ+u56cysA+/bti1zvwIG467UEEdD27dvbrvGeKc7iVe1pjI7NtvV7DsfleUQl7DxMWk8vXQreylOvbGrf6U6qfw5169YxeN+GDvtra2tZtb056Tnd/bsWEZG+af/+/QX1HdOjgZuZPQIka9O63t3/nOq0JPuSVt24+2JgMUB1dbXX1NRkk83QamtrefzLp/O+nyzl+svmsfNAI6xY1ul5c0+YA6ueZ+jQYdTUnJE0zeCVS+DA/rbt/v37w6FDzJs3j0GVZVD7KACJn/Gb/Tfxs9r1TJk8DjZtYPCQIbB3LyOGD2X9vj1t6aIVVxPGj6WmZm7kWkC/UWv44WN1lBjcfOXbufKWZ+OuP2ZwZeSeD94PwDtmj+dPL2yhevJwVmzaQzLvOOssyp96jNW72oOudIMMZs6cyfzjq9o+Y1RNTQ2HXtoGLzzf4ZzYPImIiOTK4CFDqKk5M9/ZaNOjgZu7n5fFafVA7IyzE4GtuclR140Z3I+nrl0IwN5DTaHOKc+0Yxnt03lYutEJwIdOn8yHTp/Mt/7yKtDejy6xSXRo/3L2Hmpi7NDkI0q/cO4MjhnRcYWGP35mftz2De87kbNnjuZtk4fzju/8Lem1rAQaGpPXlAEdgr5M+7h9pqbj9CoiIiK9UcH0cUvjXuByM6s0s6nADGB5nvOUVNj1QCcMjwRLb5syPPS1f/GRaj5bcyxTRg4I1S8tOn9bdAqOxD5u5x8/lps+OJd/WDg95TUSQ6Tjxg1hUkIw17+ilPeeMqGtn10yyfIbu2vO+CEJ6VNeKum8cdecPyv1CSIiIr1IwfRxM7P3AT8CRgP3m9mL7n6+u79sZncBrwDNwOfcvSXdtfIlMThKZfqYQTzyxXOYkmLN0WQmjxzINRfMBsJ1lIyOXI0OQigt6RhUvu+UiWmvkckC8umCrVTHll9/LkeOtvLLJ+P7s5kR9yEvf/skRg6KzH2XrMZNc7uJiEhfUTCBm7vfA9yT4tgNwA09m6PMJa4qkDJdSQnTxwxKmyZdLBImUCkvi6SJji6dkDDJbuKI1qiBlZFHYkBFKaMGRdYzvfjEcdz/0ra4dN+/bG7c2qvpagGT1rgRaWZOxohfZP66i45jaP/IiN3mLFasEBERyVahVQ1kFbiZ2SjgLc+kSqYPCNNUOn5ov7TNimGEOTs6Ue4pk4bx4dMnc/6csdyx/I2246l+c1edOZUWd66cP5WKshJe+8YF1O1o6BC4vf/U+Nq6dIFbskOdfYZUl9NapSIi0pd1GmmY2elmVmtmfzKzU8xsNZHJcreb2QXdn8XiEaap9JNnT+vyfcLELtGm0RZ33nvKBPpXhJuzuKKshM/WTG8bDBGtsYP0wVb6ptLMI7fYw7Gnv3vueM5LMs+diIhIdyi0GqowbXs/Bm4E7gAeAz7h7mOBs4FvdWPeik6YptLENUOzkXrMZbtoDJlqEfhMHsSJwWCKD58xOWWabJpK04ltDo5NW1FWwtfefXwnZ4uIiPROYZpKy9z9IQAz+w93Xwbg7q+pU3i8xHU8kynvYjMpEKqttDRhcALAuKH92LbvCJC6qTSZYQMq4haWTyZ94JZsb/hySHzO9NiJiEhfFab6J3aRocQ58gutBjGvwvS/ykmNW4jAJRpEtsR05n/6unP57gdOAlIPTsg6T2k+VmcBfrIg0lK8B/VzExGRvitMFDHXzPab2QHgJDM7ELN9Yjfnr6gMHVDeaTNe2Ml3044qDXF+aXCBloSo6MSJQwFYOHtMqHyElWkwlTZ5mrVKOz1XREQkhwrtK6fTwM3dS919iLsPdvey4Gd0O9yq6n3IVWdOZdzQ5NNcAJQlmU8tU2GaqKOLvyf2cZs9dgiv/scFXHLS+C7nI+5+mQZunR6P6eOWkDhVkPi7T5yWUR5ERESKTad93Mzsi+mOu/v3c5ed3uXaC2dz+rSRvPcnT7XtK81BH7dQNW4pAjcg9AjTTOS8Fsxi34br4zZ/+qgcZ0JERPq6QusTFmZwwuDg5yzg7USWoAJ4N/BEd2Sqt7j4xHEdlojKZp3SROGmA0kduHWHjJtKOztuyd9ncy8REZHeotPAzd3/HcDMHgJOdfcDwfbXgT90a+6KXLL4IuzghHRTfoSZDiQ6OKG5tbWTlLmRaUXiR+ZUpD3eSRc4ERGRPimTDlfHAEdjto8CU3Kam14mWbNo2OlAuto5v73GLdTtuqyzWrA/f+5MzjsuMiDiKxfN5pQx2a+2pho3ERHpqzIJ3G4HlpvZ183sa8AzwG25yoiZfdfMXjOzVWZ2j5kNizl2nZnVmdkaMzs/V/fsLtGBnMlqxnIxHUgY7YFbz0RuncVScycNY9royPqsYeaQi5uAV02lIiKSJ4X2jRM6iggWer8K2APsBa5y91yunPAwcIK7nwSsBa4DMLPjgcuBOcAFwE/NLPe967tB8qbSnu3j1txDfdxiA61Hrz4neZrgZ2KWks0pFz+Pm6U+KCIi0odkVP3j7s+7+w+C1wu5zIi7P+TuzcHmMiC6ivmlwJ3u3ujurwN1wLzST3v4AAAVHUlEQVRc3rsn9S8PF3N+6V2z2t6PGVwZdyxMH7e3TR7OMSMG8M/vnJlZBnPg2KBmLZUwk/+mC05zsfiEiIhIGMU4qjQfPgb8Png/gUggF1Uf7CtY6QKTUYMqUx6LdfbM0W3vn7hmQdyxMDVug/uVdzivuw2oKOUL585IeTxaK+dOiEXm1VQqIiKSqEcDNzN7BBib5ND17v7nIM31QDPw2+hpSdInjYzMbBGwCKCqqora2tquZjmthoaGpPdobIyM4Xh66VKG9Yuv1Fz17NKM53Jb9tSSuO3YKT66+zNm4qcL+4FvprZ2c9z+aB43vxEplw0bNjC56mjb/i1bGuPSr3ntNYbtq2vbfuLxx+PKrLEl/tdfSGUgIiK9y4H9+wvqeyZ04GZmlcDfERlJ2naeu/9H2Gu4+3md3OOjwCXAue5tXdjrgUkxySYCW1NcfzGwGKC6utpramrCZi0rtbW1JLtH5dJHoLGR+fPnM2ZIP3jw/rZj5y4MXwtW8tf7aXU63KOl1eGhB4COxwrFZbtWcteKeqA9j8uPvAavr2fq1KkMKtnStv/hPS/B5jfazp09ezZnzx0PDz/Ydn5s4HakqaXtWOz1gbiyFhER6arBQ4ZQU3NmvrPRJpM+bn8m0t+sGTgY88oJM7sA+DLwHnc/FHPoXuByM6s0s6nADGB5ru7bHW764MmcPm0EIwamn6usM09+eSF/+uz8DvuLoaHwOx+Y22FftIUzzKjSuPNSXEdERKS7FdpXTiZNpRPd/YJuywn8GKgEHg76Qi1z90+7+8tmdhfwCpGg8XPu3tKN+eiy+ceOYv6xHZdfShaEpTN+WH/GD+vfYX+xBC5VQyrZvr+9GTTaby1M3JZ+cEKRFICIiBS9Yh6csNTMTnT3l7ojI+4+Pc2xG4AbuuO+PenUY4bn5DphFpkvBA/98znsP9zUtl2SQY1busEJsZvvnju+CzkUEREpLpkEbmcBV5rZ60Ajke9PD+ZdkzTu+ez80KNJe5Oh/csZ2r+8fUcQgbWGiNzi1yqNj9xia9x+dMUpXcukiIhIEckkcLuw23LRy52So5q2YhcNt0I1laY7Fhy87sLZXcyRiIhIcQkduLn7pu7MiPR+75pTxQ8eXce7jq9i17qkA4PbpGsONjM2fvviXGdPRESk4HU6qtTMngx+HjCz/TGvA2a2v/uzKL3FnPFD2fjtizlhwtC4/YXW8VNERCSq0HqVd1rj5u5nBT8Hd392RCIK7R+KiIj0TYVWuZDRWqUiPaVIBs6KiIj0KAVuUpCKZcoTERGRnqTATURERKRIKHCTgqPaNhERkeQ6HZxgZseEvNZed9coU8lYpmuXioiI9JRCq0oIM4/bbUQGVaTLuwO3Ar/OQZ6kj/pg9STufr6ed8zouM6riIhIPhRa3UKY6UAW9ERGRE6cOJT//IBWUBMREUklTFPpiBDXaXX3vV3JiJl9A7gUaAV2AFe6+1aLdHj6AXARcCjY/3xX7iUiIiJSjMI0lW4NXumaSkuBsH3hUvmuu38VwMy+APwb8Gkia6TOCF6nAT8LfkovccyIAQCMGVyZ55yIiIgUtjCB26vufkq6BGb2QlczkjCwYSDtzcqXAr92dweWmdkwMxvn7tu6ek8pDIvOnsbscYOpmTk631kREREpaGECtzNylKZTZnYD8BFgHxDtWzcB2ByTrD7Yp8CtlygtMRbMGpPvbIiIiBS8MIMTjgCY2Tfd/V9jj5lZqbu3RNN0xsweAcYmOXS9u//Z3a8Hrjez64DPA18jeRNt0kEeZrYIWARQVVVFbW1tmGxlraGhodvv0Zl83z9bmZRdsX5GEREpfgf27y+o76EwNW5RE8zsCne/A8DMxgC/p71mrFPufl7IpL8D7icSuNUDk2KOTSTS5y7Z9RcDiwGqq6u9pqYmbNayUltbS3ffI6UH7wfI3/27KFTZZfMZg3NERERyYfCQIdTUnJnvbLTJZOWETwGLzGyemb0deAz4Xq4yYmYzYjbfA7wWvL8X+IhFnA7sU/82ERER6YvCTAfya+B54AXgc0Rqw5qB97p7XQ7z8m0zm0VkOpBNREaUAjxAZCqQOiLTgVyVw3uKiIiIFI2wKyfMBT4W/JwCPAt8yMxWu/vduciIu/9div1OJGAUERER6dPCDE54FHg0um1mZcDxRIK404GcBG4iIiIikl6nfdzMLG6VAndvdvdV7n67u38pWRoRERERyb0wTaXHmdmqNMcNGJqj/IiIiIhICmECt9kh0rR0NSMiIiIikl6YPm6bEvcFU3e0uvv6bsmViIiIiHSQyQS8AJjZN4BRwfuxwCfdfVeuMyYiIiIi8TIO3IBh7v4ZADObAnzXzPa5+z/lMmMiIiIiEi+TlROi2vqzuftGIhPm3pKrDImIiIhIctkEbk+a2bfMbLyZjQf6ufvKXGdMREREROJl3FTq7neb2VbgG8H5N+Q8VyIiIiLSQejAzcweAa5295XuvhRY2n3ZEhEREZFEmTSVXgPcZGa3mNm47sqQiIiIiCQXOnBz9+fdfSFwH/CgmX3NzPrnOkNm9iUzczOLTjliZvZDM6szs1Vmdmqu7ykiIiJSDDLq42ZmBqwBfgZ8E/ikmV3n7rfnIjNmNgl4J/BGzO4LgRnB67Tg3qfl4n7F7LGrz2FQv2xmcxEREZFiFbrGzcyeBLYANwETgCuBGmCemS3OUX5uItIk6zH7LgV+7RHLgGFqqoVpowcxZnC/fGdDREREelAmVTafBl52d0/Y/w9m9mpXM2Jm7wG2uPvKSMVemwnA5pjt+mDftiTXWAQsAqiqqqK2trar2UqroaGh2+/RW2VSdipjERHJlwP79xfU91CngZuZHRO83Q9MSgiqoi4ysyHuvr+Taz0CjE1y6HrgK8C7kp2WZF9i8BjZ6b4YWAxQXV3tNTU16bLTZbW1tXT3PXqrUGX34P0AmZVxcI6IiEguDB4yhJqaM/OdjTZhatxuIxIoJY3YAg7cCvw63YXc/bxk+83sRGAqEK1tmwg8b2bziNSwTYpJPhHYGiLfIiIiIr1Kp4Gbuy/o7ky4+0vAmOi2mW0Eqt19l5ndC3zezO4kMihhn7t3aCYVERER6e2KYVjiA8BFQB1wCLgqv9kRERGRvmLaqIH5zkKcggzc3H1KzHsHPpe/3EgxWXLNAjbvOYQ7HGxspl95KWWlMa38DoeOtlBaYlSWxw+qbmxupbnFGVhRGtcxoKXVOXy0hQEVZZTEnBK9R2V5KeWl8T0JDh9twQz6lZfG7T/a3EpTizOgohRLco/+FaWUlsTn9+DRFspLjYqy9puvfHEls+acgDv0r4i/R1OL09jUwsDKsrh7tLbCoaPJy+Tg0RbKSo3KsoQyaWqlpTWS39gyaW5xjjRlUCZBuZeUdCyTbMr90NEWKspKOtzjcFML0LFMjja3crS5lUGVZaxcuZK5J8+NK5Ow5Q5wpKklUu7l8fnNptxTPos5LPfDTcmfxWzK/ZnnXuRtp8xN+rxDxzJJ9bx3Vu7JnsUjTS20ttKhTAq13JM97ytfXMnsE05MW+6JZRJ93pM9i23lnuJ5T1UmxVjuL61q/3ebzfMefRbDPu+xfwNOmzqSQlKQgZtItiaNGMCkEQPynY1ud3RzKTWzq/KdjaJ0tL6U+ceOync2ilJTfRnnzByd72wUraObS6mZNabzhNJB8xb9u43KZMkrEREREckjBW4iIiIiRUKBm4iIiEiRUOAmIiIiUiQUuImIiIgUCQVuIiIiIkVCgZuIiIhIkVDgJiIiIlIkFLiJiIiIFAkFbiIiIiJFomACNzP7upltMbMXg9dFMceuM7M6M1tjZufnM58iIiIi+VJoa5Xe5O7fi91hZscDlwNzgPHAI2Y2091b8pFBERERkXwpmBq3NC4F7nT3Rnd/HagD5uU5TyIiIiI9rtACt8+b2Sozu9nMhgf7JgCbY9LUB/tERERE+hRz9567mdkjwNgkh64HlgG7AAe+AYxz94+Z2U+Ap939N8E1fgU84O5/THL9RcAigKqqqrfdeeed3fNBAg0NDQwaNKhb79FbhSm7Kx88CMCtFwzsiSwVFT172VPZZU9l1zUqv+z1hbJbsGDBc+5e3Vm6Hu3j5u7nhUlnZr8A7gs264FJMYcnAltTXH8xsBigurraa2pqss5rGLW1tXT3PXqrUGX34P0AKuMk9OxlT2WXPZVd16j8sqeya1cwTaVmNi5m833A6uD9vcDlZlZpZlOBGcDyns6fiIiISL4V0qjS75jZyUSaSjcCnwJw95fN7C7gFaAZ+JxGlIqIiEhfVDCBm7t/OM2xG4AbejA7IiIiIgWnYJpKRURERCQ9BW4iIiIiRUKBm4iIiEiRUOAmIiIiUiQUuImIiIgUCQVuIiIiIkVCgZuIiIhIkVDgJiIiIlIkFLiJiIiIFAkFbiIiIiJFQoGbiIiISJFQ4CYiIiJSJAoqcDOzfzCzNWb2spl9J2b/dWZWFxw7P595FBEREcmXsnxnIMrMFgCXAie5e6OZjQn2Hw9cDswBxgOPmNlMd2/JX25FREREel4h1bh9Bvi2uzcCuPuOYP+lwJ3u3ujurwN1wLw85VFEREQkbwopcJsJvMPMnjGzx83s7cH+CcDmmHT1wT4RERGRPqVHm0rN7BFgbJJD1wd5GQ6cDrwduMvMpgGWJL2nuP4iYBFAVVUVtbW1Och1ag0NDd1+j94qk7JTGXekZy97Krvsqey6RuWXPZVdux4N3Nz9vFTHzOwzwJ/c3YHlZtYKjCJSwzYpJulEYGuK6y8GFgNUV1d7TU1NjnKeXG1tLd19j94qVNk9eD+AyjgJPXvZU9llT2XXNSq/7Kns2hVSU+n/AAsBzGwmUAHsAu4FLjezSjObCswAluctlyIiIiJ5UjCjSoGbgZvNbDVwFPhoUPv2spndBbwCNAOf04hSERER6YsKJnBz96PAh1IcuwG4oWdzJPk2tH85rZ60O6OIiEifVDCBm0ii5/41ZZdIERGRPkmBmxSsstJC6oIpIiKSf/pmFBERESkSCtxEREREioQCNxEREZEiocBNREREpEgocBMREREpEua9dJ4sM9sJbOrm24wisrqDZE5l1zUqv+yp7LKnsusalV/2+kLZTXb30Z0l6rWBW08wsxXuXp3vfBQjlV3XqPyyp7LLnsqua1R+2VPZtVNTqYiIiEiRUOAmIiIiUiQUuHXN4nxnoIip7LpG5Zc9lV32VHZdo/LLnsouoD5uIiIiIkVCNW4iIiIiRUKBW5bM7AIzW2NmdWZ2bb7zU2jMbJKZ/c3MXjWzl83sH4P9I8zsYTNbF/wcHuw3M/thUJ6rzOzU/H6C/DOzUjN7wczuC7anmtkzQdn93swqgv2VwXZdcHxKPvOdb2Y2zMzuNrPXgufvDD134ZnZPwf/Zleb2R1m1k/PXnJmdrOZ7TCz1TH7Mn7WzOyjQfp1ZvbRfHyWfEhRft8N/u2uMrN7zGxYzLHrgvJbY2bnx+zvU9/HCtyyYGalwE+AC4HjgSvM7Pj85qrgNANXu/txwOnA54IyuhZ41N1nAI8G2xApyxnBaxHws57PcsH5R+DVmO3/BG4Kym4P8PFg/8eBPe4+HbgpSNeX/QB40N1nA3OJlKGeuxDMbALwBaDa3U8ASoHL0bOXyq3ABQn7MnrWzGwE8DXgNGAe8LVosNcH3ErH8nsYOMHdTwLWAtcBBN8flwNzgnN+Gvznts99Hytwy848oM7dN7j7UeBO4NI856mguPs2d38+eH+AyJfnBCLldFuQ7DbgvcH7S4Ffe8QyYJiZjevhbBcMM5sIXAz8Mtg2YCFwd5AkseyiZXo3cG6Qvs8xsyHA2cCvANz9qLvvRc9dJsqA/mZWBgwAtqFnLyl3fwLYnbA702ftfOBhd9/t7nuIBC6JwUyvlKz83P0hd28ONpcBE4P3lwJ3unuju78O1BH5Lu5z38cK3LIzAdgcs10f7JMkguaTU4BngCp33waR4A4YEyRTmcb7L+AaoDXYHgnsjfmDFls+bWUXHN8XpO+LpgE7gVuCZuZfmtlA9NyF4u5bgO8BbxAJ2PYBz6FnLxOZPmt6BlP7GPCX4L3KL6DALTvJ/kep4blJmNkg4I/AP7n7/nRJk+zrk2VqZpcAO9z9udjdSZJ6iGN9TRlwKvAzdz8FOEh7U1UyKrsYQRPdpcBUYDwwkEgTVCI9e5lLVVYqwyTM7HoiXW5+G92VJFmfLD8FbtmpBybFbE8EtuYpLwXLzMqJBG2/dfc/Bbu3R5uigp87gv0q03ZnAu8xs41Eqv0XEqmBGxY0X0F8+bSVXXB8KB2bb/qKeqDe3Z8Jtu8mEsjpuQvnPOB1d9/p7k3An4D56NnLRKbPmp7BBMEAjUuAv/f2OctUfgEFbtl5FpgRjLSqINJh8t4856mgBP1cfgW86u7fjzl0LxAdNfVR4M8x+z8SjLw6HdgXbW7oa9z9Onef6O5TiDxbj7n73wN/Az4QJEssu2iZfiBI36v/x5mKu78JbDazWcGuc4FX0HMX1hvA6WY2IPg3HC0/PXvhZfqs/RV4l5kND2o83xXs65PM7ALgy8B73P1QzKF7gcuDkcxTiQzyWE5f/D52d72yeAEXERnxsh64Pt/5KbQXcBaR6upVwIvB6yIi/V8eBdYFP0cE6Y3IyKD1wEtERrXl/XPk+wXUAPcF76cR+UNVB/wBqAz29wu264Lj0/Kd7zyX2cnAiuDZ+x9guJ67jMrv34HXgNXA7UClnr2UZXUHkb6ATURqfj6ezbNGpC9XXfC6Kt+fK8/lV0ekz1r0e+O/Y9JfH5TfGuDCmP196vtYKyeIiIiIFAk1lYqIiIgUCQVuIiIiIkVCgZuIiIhIkVDgJiIiIlIkFLiJiIiIFAkFbiIiIiJFQoGbiIiISJFQ4CYiEoKZnWdmt+c7HyLStylwExEJZy7wQr4zISJ9mwI3EZFw5gIvBGsl3mpmNwbreYqI9JiyfGdARKRIzAV2EFkA/Jfu/ps850dE+iCtVSoi0gkzKwd2AZuAT7n703nOkoj0UWoqFRHp3PHAs0Az0JLnvIhIH6bATUSkc3OBpcDlwC1mVpXn/IhIH6XATUSkc3OB1e6+FvgycFfQfCoi0qPUx01ERESkSKjGTURERKRIKHATERERKRIK3ERERESKhAI3ERERkSKhwE1ERESkSChwExERESkSCtxEREREioQCNxEREZEi8f8BH71bvXRsLt8AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 720x216 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAnIAAADiCAYAAADK1yr/AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAHcxJREFUeJzt3Xu4XHV56PHvSwiSajVgMJANGKw8KBYBT0RQn55db+ClEj1owUvBS9FT0cfnKAJaL/UGlj6t1SqVUhS1XnqQRlQ0KnQfqwQlEq5qJHKRJFQksMHAFkLynj9mbZhsZmbP7MxlrZnv53nm2TO/9Zvf7/fuN7Py7rXWzERmIkmSpOrZadALkCRJ0txYyEmSJFWUhZwkSVJFWchJkiRVlIWcJElSRVnISZIkVZSFnKTKi4jrImK8D/NMRMQbm2yLiPhsRNwZET/p9Vp6JSJeHRHfHfQ6JLXHQk4aURFxQkRcExH3RsR/R8RZEbGwg+ffFBHP6+J65jxeZj4lMye6tZY5ejbwfGDvzDxswGuZs8z8t8x8Qa/niYilEZERsXOv55KGmYWcNIIi4h3Ax4CTgccAhwOPB74XEbsMcm0V9njgpsy8p9Mn9rKYsVCShpuFnDRiIuLRwN8Ab83M72Tmlsy8CXgltWLkNUW/z0XEh+ueNx4R64v7XwD2Bb4REZsj4l11R1hOjIiNEXFrUTAyl/EarHtRRHwzIiYj4o6I+K+I2KnY9uDRvIhYEBHnFac4f16sbX3dODdFxDsj4uqIuCsivhoRuxbbdivm+G3x/G9GxN5t/E7fAJwDHFGs/2+K9r+MiHXFei+MiCV1z8mIeEtEXA9c32DMfyrGmr49EBEfKLYtiYivFeu8MSLeVve8D0TE+RHxxYi4GzghIh4RER8v8rKxuP+IJrGcEBE/nLHON0fE9cXv5FMREXV9fxQRnyx+l7+IiOfO+F0/r+7xByLii8XDHxQ/J4v4jpjt9yzp4SzkpNHzTGBX4IL6xszcDHyb2unBljLztcCvgT/LzEdl5t/Wbf5TYH/gBcCp7ZwunWW8ae8A1gN7AIuBdwONvmPw/cBS4AlFLK9p0OeVwFHAfsBTgROK9p2Az1IraPcFpoB/amP9/wq8GVhVrP/9EfEc4PRirr2Am4GvzHjqcuAZwIENxjypGOtR1E7b3gl8vShevwFcBYwBzwXeHhFH1j39aOB8YCHwb8B7qB11PQQ4GDgM+OvZ4qrzEuDpxXNfCdTP9QzgBmARtd/9BRGxextj/knxc2ER56oO1iOpYCEnjZ5FwO2Z+UCDbbcW23fE32TmPZl5DbWi6LgdHG/aFmoF0eOLo4j/lY2/LPqVwEcz887MXA98okGfT2Tmxsy8g1pRdAhAZm7KzK9l5r2Z+TvgI8D/nON6Xw2cm5lXZOZ9wGnUjtgtretzembekZlTzQaJiD2AFdSOoK6hVlDtkZkfzMz7M/MG4F+AY+uetiozV2TmtmLsVwMfzMzbMvO31I7IvraDWM7IzMnM/DXwnxS/r8JtwMeLnHwVWAu8uIOxJe0ACzlp9NwOLGpy7dRexfYdcUvd/ZuBJc06duhMYB3w3Yi4ISJObdJvyYw13NKgz3/X3b8XeBRARPxBRHwmIm4uTkv+AFgYEfPmsN4l1OIHHjziuYnaUbRWa3tQRMyndmTtS5k5fTTv8cCS4hTzZERMUjs6ubjFuNuthc7z0vD3Vdgwo6DuZs4lzcJCTho9q4D7gJfXN0bEI4EXAhcXTfcAf1DXZc8Z4zQ6GgawT939fYGNOzhebWPm7zLzHZn5BODPgP9Tfz1WnVuB+uva9mnQp5l3AAcAz8jMR/PQ6b/oYIxpG6kVXbUBar/fxwIb6vq0jBn4JPA7tj8NegtwY2YurLv9YWa+qMW4262F7fOyo8amr5lrMHarnM8Wu6Q2WMhJIyYz76J2au2TEXFURMwvTvf9X2rXoH2h6Hol8KKI2D0i9gTePmOo31C7Dm2m9xZHtp4CvA746g6OB0BEvCQinlgUDXcDW4vbTP8OnFa8cWEMOKnZmA38IbXr4iaL67ze38FzZ/oS8LqIOKR4Y8FHgR8XbyyZVUS8idpp3Vdl5ra6TT8B7o6IU4o3dsyLiD+OiKe3GO7LwF9HxB4RsQh4H/DFFv078TjgbcW/o1cATwYuKrZdCRxbbFsGHFP3vN8C22iRc0mzs5CTRlDxZoJ3A39HrSj6MbUjPc8trueCWkF3FXAT8F0eKsimnU6tOJiMiHfWtf8/aqdALwb+LjOnP1x2ruNN2x/4PrCZ2lHFTzf57LgPUitIbyz6n0/tCGQ7Pg4soHZ6+TLgO20+72Ey82LgvcDXqB0l/CO2v45tNsdRK3I21r1z9d2ZuZXaEclDqMV4O7V3zD6mxVgfBlYDVwPXAFcUbd3wY2q5uZ3aNYXHZOamYtt7qcV9J7U/Hr40/aTMvLfo/6Mi54d3aT3SSInG1wpLUmeKo3o3AvObvJFiICLifwPHZuZc37SgJiLiBOCNmfnsQa9FGlUekZM0VCJir4h4VkTsFBEHULvu7T8GvS5J6gU/8VvSsNkF+Ay1z4ibpPbZbZ8e6IokqUc8tSpJklRRnlqVJEmqKAs5SZKkihqZa+QWLVqUS5cu7ekc99xzD4985CN7OofmxtyUm/kpL3NTbuan3HYkPz/96U9vz8w9Zus3MoXc0qVLWb16dU/nmJiYYHx8vKdzaG7MTbmZn/IyN+VmfsptR/ITETfP3stTq5IkSZVlISdJklRRI3NqVe1bsWYDZ65cy8bJKZYsXMDJRx7A8kPHBr2snjNu4zbu4WXcxj2scZeukIuIc4GXALdl5h832B7APwIvAu4FTsjMK/q7yuG1Ys0GTrvgGqa21L6LfMPkFKddcA3A0L4IwLiN27iNe/gY92jEXcZTq58Djmqx/YXUvqB5f+BE4Kw+rKmlFWs28KwzLuGE79zDs864hBVrNgx6SXN25sq1D/7jnza1ZStnrlz7sL7Tce936reMu6KM+yHGvb1R3a/B8OTcuB8yzHGX7ohcZv6g+PLtZo4GPp+1r6S4LCIWRsRemXlrXxY4w7BV/hsnp9pqN27jNm7jrop244bhit24Z28fhrhLV8i1YQy4pe7x+qLtYYVcRJxI7agdixcvZmJiouuL+dDEvUxt2f5rzqa2bOVDX7+KhXdd3/X5em33XYNNv3/417btvmts9/urWtybN29umf9hjXs2ZYl7tvx0W1ni7jfjfnj7zH93VYu91WtnmONupUxx92PfVsVCLhq0NfzC2Mw8GzgbYNmyZdmLz9q54zvfatz++6zkZ/u89zHb/3UCsGD+PN579EGM1/11UrW4Z/ssn2GNezZlibvfn4VVlrj7zbhbxw3Vi73Va2eY426lTHH3Y99WxmvkZrMe2Kfu8d7AxgGthSULF3TUXnbLDx3j9JcfxC7zav80xhYu4PSXH/SwQ8zGTVvtZWfcxm3c2xum2I17NOKuYiF3IfAXUXM4cNegro8DOPnIA1gwf952bQvmz+PkIw8Y0Ip23PJDxzh034U8Y7/d+dGpz2n4j9+4H2Lc1WTcxt3IsMVu3MMfd+lOrUbEl4FxYFFErAfeD8wHyMx/Bi6i9tEj66h9/MjrBrPSmul/HO86/2ru37qNsSH/vJppxm3cxj28RjVuGN3Yjbu6cZeukMvM42bZnsBb+rSctiw/dIwv/+TXTE5OsvKU5wx6OX0zHTfAV990xIBX0z/GbdyjYFT3a2DOwbirpIqnViVJkoSFnCRJUmVZyEmSJFWUhZwkSVJFWchJkiRVlIWcJElSRVnISZIkVZSFnCRJUkVZyEmSJFWUhZwkSVJFWchJkiRVlIWcJElSRVnISZIkVZSFnCRJUkVZyEmSJFWUhZwkSVJFWchJkiRVlIWcJElSRVnISZIkVZSFnCRJUkVZyEmSJFWUhZwkSVJFWchJkiRVlIWcJElSRVnISZIkVZSFnCRJUkVZyEmSJFVU6Qq5iDgqItZGxLqIOLXB9hMi4rcRcWVxe+Mg1ilJkjRoOw96AfUiYh7wKeD5wHrg8oi4MDN/NqPrVzPzpL4vUJIkqUTKdkTuMGBdZt6QmfcDXwGOHvCaJEmSSqlshdwYcEvd4/VF20z/KyKujojzI2Kf/ixNkiSpXEp1ahWIBm054/E3gC9n5n0R8WbgPOA5DQeLOBE4EWDx4sVMTEx0canbm5ycYuvWrT2do58mJ6cAZo2n3X6Dtnnz5rbWOGxxt2vQcbebn24bdNyD0knco7hf67TvILXz2hnGuNtRhrj7sW8rWyG3Hqg/wrY3sLG+Q2Zuqnv4L8DHmg2WmWcDZwMsW7Ysx8fHu7bQmc5au4rJyUl6OUc/nbV2FQDj40d0pd+gTUxMtJWbYYu7XYOOu938dNug4x6UTuIexf1ap30HqZ3XzjDG3Y4yxN2PfVvZTq1eDuwfEftFxC7AscCF9R0iYq+6hy8Fft7H9UmSJJVGqY7IZeYDEXESsBKYB5ybmddFxAeB1Zl5IfC2iHgp8ABwB3DCwBYsSZI0QKUq5AAy8yLgohlt76u7fxpwWr/XJUmSVDZlO7UqSZKkNlnISZIkVZSFnCRJUkVZyEmSJFWUhZwkSVJFWchJkiRVlIWcJElSRVnISZIkVZSFnCRJUkVZyEmSJFWUhZwkSVJFWchJkiRVlIWcJElSRVnISZIkVZSFnCRJUkVZyEmSJFWUhZwkSVJFWchJkiRVlIWcJElSRVnISZIkVZSFnCRJUkXt3GpjROzexhjbMnOyS+uRJElSm1oWcsDG4hYt+swD9u3aiiRJktSW2Qq5n2fmoa06RMSaLq5HkiRJbZrtGrkj2hijnT6SJEnqspaFXGb+HiAiPjxzW0TMq+8jSZKk/mr3XatjEXHc9IOIeBzw/d4sSZIkSe2Y7Rq5aW8CVkbEr4AEPguc0rNVSZIkaVazffzI54ErgDXAW4AvAQ8AyzNzXS8WFBFHAf9I7d2w52TmGTO2PwL4PPA/gE3An2fmTb1YS7etWLOBM1euZePkFEsWLuDkIw9g+aFjzt2HuTdMTjF22SUjE/eg53du53bu/sw9avu2UZ27ldmOyJ0HHAy8vvi5FLgceE1EXJuZ53dzMcV1d58Cng+sBy6PiAsz82d13d4A3JmZT4yIY4GPAX/ezXX0woo1GzjtgmuY2rIVgA2TU5x2wTUAPf+H4NyjNfeg53du53Zu53bu/pntzQ4XZ+bfZ+bxmXkIsAh4B/Ar4PAerOcwYF1m3pCZ9wNfAY6e0edoagUmwPnAcyOi1efclcKZK9c++A9g2tSWrZy5cq1zO/dQze/czu3czu3c/ROZOeg1PCgijgGOysw3Fo9fCzwjM0+q63Nt0Wd98fhXRZ/bW429bNmyXL16dc/W/tnj3sqiW2/kqfs0/jKMy27Y1PS5hz/hsQ3bb9p0DwBLH/vIlnPP1m8uc//s1rsBOHCvR7ece7Z+xt3+3N2Key7zdyvuucx906Z7uO+++zhgSesvkhnGuKH8+X7ggQea7tfmMne7cbfTt1dxt9PXuNubexjjvuExY3zmqbVjTAHceMaLm447MTHB+Pj4rPM3EhE/zcxls/Wb7Rq5KzLzaTvapwONjqzNrDTb6VPrGHEicCLA4sWLmZiY2KHFtbJL3s8uOyWTk42/rWz+TsGWbQ9f5vydoulz7r53GwCT87a0nHu2fnOZe14WYzbZ3m4/425/7m7FPZf5uxX3XOa++95tZDZ/7dT3g+GKG8qf751a7NfmMne7cbfTt1dxt9PXuNube9jj3n3XaFlbbN68uae1B8x+jdyTI+LqFtsDeEwX17Me2Kfu8d7UviKsUZ/1EbFzMf8djQbLzLOBs6F2RG6uVXFbxseZmJjg4CZz3LhmAx+oO78OsGD+PE5/+UEc3OT8+kc/swqAr76p9Wcuz9ZvLnN3i3H3P+65zt8tnc790c+sYnJykpWnvLDluMMYN5Q/3632a3OZu9242+k7jPlup69xlyfu9x59EOMt5t6RI3Ltmq2Qe1KDtkcA24Dpknlrgz5zdTmwf0TsB2wAjgVeNaPPhcDxwCrgGOCSLNP54SamL4YcxDtenLt4Z9eIxD3o+Z3buZ27f3OP0r5tVOeeTctCLjNvrn8cER8CHkvtSNyewF/Odm1aJzLzgYg4CVhJ7eNHzs3M6yLig8DqzLwQ+FfgCxGxjtqRuGO7NX+vLT90bGBJH/W5+/FXUbO5B6UMv3fndm7n7u3co7ZvG9W5W2n3A4GnLczMvwKIiKXAmRFxV2a+vVsLysyLgItmtL2v7v7vgVd0az5JkqSqavcruqY9eBq1+BDebcDnurgeSZIktanTQu6HEXF6RCyJiCXArpl5ZS8WJkmSpNY6OrWamedHxEbgQ8VzP9KTVUmSJGlWbRVyEfF94B2ZeVVmXgpc2ttlSZIkaTbtnlp9F/APEfHZiNirlwuSJElSe9oq5DLzisx8DvBN4DsR8f6IWNDbpUmSJKmVtt/sUHwx/VrgLOCtwPXFd6FKkiRpANoq5CLih9S+aeEfgDHgBGAcOCwizu7V4iRJktRcu+9afTNwXYOvwnprRPy8y2uSJElSG9oq5DLz2habX9yltUiSJKkDnX4g8MNk5g3dWIgkSZI6s8OFnCRJkgbDQk6SJKmiLOQkSZIqykJOkiSpoizkJEmSKspCTpIkqaIs5CRJkirKQk6SJKmiLOQkSZIqykJOkiSpoizkJEmSKspCTpIkqaIs5CRJkirKQk6SJKmiLOQkSZIqykJOkiSpoizkJEmSKqo0hVxE7B4R34uI64ufuzXptzUirixuF/Z7nZIkSWVRmkIOOBW4ODP3By4uHjcylZmHFLeX9m95kiRJ5VKmQu5o4Lzi/nnA8gGuRZIkqfR2HvQC6izOzFsBMvPWiHhck367RsRq4AHgjMxc0WzAiDgROBFg8eLFTExMdHnJ29u8eXNX55icnAKYdcx2+1VFL+Ludm56YZTzvXXr1pGMG8qf70Ht1zrtW3a9irvs+7ZRzfe0fuSnr4VcRHwf2LPBpvd0MMy+mbkxIp4AXBIR12Tmrxp1zMyzgbMBli1bluPj450uuSMTExN0c46z1q4CYHz8iK70q4pexN3t3PTCKOd7cnJy1vwMY9xQ/nwPar/Wad+y61XcZd+3jWq+p/UjP30t5DLzec22RcRvImKv4mjcXsBtTcbYWPy8ISImgEOBhoWcJEnSMCvTNXIXAscX948Hvj6zQ0TsFhGPKO4vAp4F/KxvK5QkSSqRMhVyZwDPj4jrgecXj4mIZRFxTtHnycDqiLgK+E9q18hZyEmSpJFUmjc7ZOYm4LkN2lcDbyzuXwoc1OelSZIklVKZjshJkiSpAxZykiRJFWUhJ0mSVFEWcpIkSRVlISdJklRRFnKSJEkVZSEnSZJUURZykiRJFWUhJ0mSVFEWcpIkSRVlISdJklRRFnKSJEkVZSEnSZJUURZykiRJFWUhJ0mSVFEWcpIkSRVlISdJklRRFnKSJEkVZSEnSZJUURZykiRJFWUhJ0mSVFEWcpIkSRVlISdJklRRFnKSJEkVZSEnSZJUURZykiRJFWUhV1Ir1mxgza8n+fGNd/CsMy5hxZoNg15SXxj3aMa99s5tIxn3qObbuI1b3VOaQi4iXhER10XEtohY1qLfURGxNiLWRcSp/Vxjv6xYs4HTLriG+7duA2DD5BSnXXDN0L8IjNu4jXt4Gbdxj0Lcg1CaQg64Fng58INmHSJiHvAp4IXAgcBxEXFgf5bXP2euXMvUlq3btU1t2cqZK9c+rO8w/cVj3A8x7u0Z92jFDcMTu3E/ZBTiHoSdB72AaZn5c4CIaNXtMGBdZt5Q9P0KcDTws54vsI82Tk611d7sLx6A5YeO9XaRPWDcrduN27hHIW4YrtiNe/b2YYp7EMp0RK4dY8AtdY/XF21DZcnCBW21d/oXT9kZd+t2436IcVdPu3HDcMVu3LO3D1Pcg9DXI3IR8X1gzwab3pOZX29niAZt2WK+E4ETARYvXszExEQ7y5yzzZs3d2WOF++7lc/dDfdve6htl51q7fXjb2jyF8+Gyamex9oLvYy7W7npBfP9UJtxTzzYVpa4+71fg/LE3g29jrus+7ZRzfdM/chPXwu5zHzeDg6xHtin7vHewMYW850NnA2wbNmyHB8f38HpW5uYmKAbc4wDB67ZwJkr17JxcoolCxdw8pEHPOwQ89hllzR8AYwtXNCVdfTbOL2Lu1u56YVxzPeGySnGRjDuKuS73/s1KE/s3TBOb+Mu675tnNHM90z9yE9prpFr0+XA/hGxH7ABOBZ41WCX1BvLDx2b9dqAk488gNMuuGa7Q9IL5s/j5CMP6PXyesa4mxvmuFvt7IY57lZGNW4YvtiNu7Vhi7vfSlPIRcTLgE8CewDfiogrM/PIiFgCnJOZL8rMByLiJGAlMA84NzOvG+CyB2r6BdLOXzzDxLiN27iH26jGbtyjFXe3RGbTS8yGyrJly3L16tU9naOsh7hlbsrO/JSXuSk381NuO5KfiPhpZjb9XN1pVXvXqiRJkgoWcpIkSRVlISdJklRRI3ONXET8Fri5x9MsAm7v8RyaG3NTbuanvMxNuZmfctuR/Dw+M/eYrdPIFHL9EBGr27kwUf1nbsrN/JSXuSk381Nu/ciPp1YlSZIqykJOkiSpoizkuuvsQS9ATZmbcjM/5WVuys38lFvP8+M1cpIkSRXlETlJkqSKspBrISLOjYjbIuLaurbdI+J7EXF98XO3oj0i4hMRsS4iro6Ip9U95/ii//URcfwgYhlGHebnSRGxKiLui4h3zhjnqIhYW+Tu1H7HMYw6zM2ri9fM1RFxaUQcXPccc9MDHebn6CI3V0bE6oh4dt1z3Lf1QCf5qdv+9IjYGhHH1LWZny7r8LUzHhF3Fa+dKyPifXXP6d6+LTO9NbkBfwI8Dbi2ru1vgVOL+6cCHyvuvwj4NhDA4cCPi/bdgRuKn7sV93cbdGzDcOswP48Dng58BHhnXf95wK+AJwC7AFcBBw46tqrfOszNM6dfE8AL61475qYc+XkUD12G81TgF8V9920lyE/xeB5wCXARcIz5KUdugHHgmw3G6Oq+zSNyLWTmD4A7ZjQfDZxX3D8PWF7X/vmsuQxYGBF7AUcC38vMOzLzTuB7wFG9X/3w6yQ/mXlbZl4ObJnR/zBgXWbekJn3A18pxtAO6DA3lxavDYDLgL2L++amRzrMz+Ys/vcBHglM33ff1iMd/t8D8Fbga8BtdW3mpwfmkJtGurpvs5Dr3OLMvBWg+Pm4on0MuKWu3/qirVm7eqNZfpoxP/3TTm7eQO3INpibfmuan4h4WUT8AvgW8Pqi2fz0V8P8RMQY8DLgn2f0Nz/902rfdkREXBUR346IpxRtXc2NhVz3RIO2bNGucjA/JRERf0qtkDtluqlBN3MzAJn5H5n5JGpHGj5UNJufcvg4cEpmbp3Rbn4G7wpqX7N1MPBJYEXR3tXcWMh17jfFKVOKn9OHstcD+9T12xvY2KJdvdEsP82Yn/5pmpuIeCpwDnB0Zm4qms1Nf8362ilOK/1RRCzC/PRbs/wsA74SETcBxwCfjojlmJ9+apibzLw7MzcX9y8C5vfitWMh17kLgel3/xwPfL2u/S+Kd68eDtxVHGJdCbwgInYr3snygqJNvdEsP81cDuwfEftFxC7AscUY6r6GuYmIfYELgNdm5i/r+pub/mqWnydGRBT3n0bt4uxNuG/rt4b5ycz9MnNpZi4Fzgf+KjNXYH76qdlrZ8+6185h1GquTXR73zbod4CU+QZ8GbiV2gXy66md9nkscDFwffFz96JvAJ+i9k6Ua4BldeO8HlhX3F436LiG5dZhfvYs+twNTBb3H11sexHwyyJ37xl0XMNw6zA35wB3AlcWt9V145ibwefnFOC6IjergGfXjeO+bcD5mfG8z1G8a9X8DD43wEnFa+cqam/kembdOF3bt/nNDpIkSRXlqVVJkqSKspCTJEmqKAs5SZKkirKQkyRJqigLOUmSpIqykJMkSaooCzlJkqSKspCTpA5FxPMi4guDXockWchJUucOBtYMehGSZCEnSZ07GFgTEY+IiM9FxEenv1NRkvpp50EvQJIq6GDgNmpfQn5OZn5xwOuRNKL8rlVJ6kBEzAduB24G3pSZqwa8JEkjzFOrktSZA4HLgQeArQNei6QRZyEnSZ05GLgUOBb4bEQsHvB6JI0wCzlJ6szBwLWZ+UvgFODfi9OtktR3XiMnSZJUUR6RkyRJqigLOUmSpIqykJMkSaooCzlJkqSKspCTJEmqKAs5SZKkirKQkyRJqigLOUmSpIr6/8eUyWTiP6nfAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 720x216 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "def uniform_midtread_quantizer(x, xmin=1):\n",
    "    # limiter\n",
    "    x = np.copy(x)\n",
    "    if x <= -xmin:\n",
    "        x = -1\n",
    "    if x > xmin - Q:\n",
    "        x = 1 - Q\n",
    "    # linear uniform quantization\n",
    "    xQ = Q * np.floor(x/Q + 1/2)\n",
    "    \n",
    "    return xQ\n",
    "\n",
    "\n",
    "# compute input signal\n",
    "x = np.random.uniform(low=-1, high=1, size=256)\n",
    "x = np.concatenate((x, np.zeros(1024)))\n",
    "# compute output signal\n",
    "yQ = sos_df1(x, 2*a, requantize=uniform_midtread_quantizer)\n",
    "\n",
    "# plot results\n",
    "plt.figure(figsize=(10, 3))\n",
    "plt.plot(20*np.log10(np.abs(yQ)))\n",
    "plt.title('Level of output signal')\n",
    "plt.xlabel(r'$k$')\n",
    "plt.ylabel(r'$|y_Q[k]|$ in dB')\n",
    "plt.grid()\n",
    "\n",
    "plt.figure(figsize=(10, 3))\n",
    "k = np.arange(1000, 1050)\n",
    "plt.stem(k, yQ[k])\n",
    "plt.title('Output signal for zero input')\n",
    "plt.xlabel(r'$k$')\n",
    "plt.ylabel(r'$y_Q[k]$ ')\n",
    "plt.grid();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise**\n",
    "\n",
    "* Determine the period of the large limit cycles. How is it related to the poles of the system?\n",
    "\n",
    "Solution: The period of the large limit cycles can be estimated from the second illustration as $P = 6$. The normalized frequency of a harmonic exponential signal with the same periodicity is given as $\\Omega_0 = \\frac{2 \\pi}{P} = \\frac{\\pi}{3}$. The poles of the system can be extracted from the code of the first example as $z_{\\infty 0,1} = 0.9 \\cdot e^{\\pm j \\frac{\\pi}{3}}$. The periodicity of the large limit cycles is hence linked to the normalized frequency of the poles."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbsphinx": "hidden"
   },
   "source": [
    "**Copyright**\n",
    "\n",
    "This notebook is provided as [Open Educational Resource](https://en.wikipedia.org/wiki/Open_educational_resources). Feel free to use the notebook for your own purposes. The text is licensed under [Creative Commons Attribution 4.0](https://creativecommons.org/licenses/by/4.0/), the code of the IPython examples under the [MIT license](https://opensource.org/licenses/MIT). Please attribute the work as follows: *Sascha Spors, Digital Signal Processing - Lecture notes featuring computational examples, 2016-2018*."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}