{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "nbsphinx": "hidden"
   },
   "source": [
    "# Quantization of Signals\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 Error of a Linear Uniform Quantizer\n",
    "\n",
    "As illustrated in the [preceding section](linear_uniform_characteristic.ipynb), quantization results in two different types of distortions. Overload distortions are a consequence of exceeding the minimum/maximum amplitude of the quantizer. Granular distortions are a consequence of the quantization process when no clipping occurs. Various measures are used to quantify the distortions of a quantizer. We limit ourselves to the signal-to-noise ratio as commonly used measure."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "### Signal-to-Noise Ratio\n",
    "\n",
    "A quantizer can be evaluated by its [signal-to-noise ratio](https://en.wikipedia.org/wiki/Signal-to-noise_ratio) (SNR), which is defined as the power of the continuous amplitude signal $x[k]$ divided by the power of the quantization error $e[k]$. Under the assumption that both signals are drawn from a zero-mean wide-sense stationary (WSS) process, the average SNR is given as\n",
    "\n",
    "\\begin{equation}\n",
    "SNR = 10 \\cdot \\log_{10} \\left( \\frac{\\sigma_x^2}{\\sigma_e^2} \\right) \\quad \\text{ in dB}\n",
    "\\end{equation}\n",
    "\n",
    "where $\\sigma_x^2$ and $\\sigma_e^2$ denote the variances of the signals $x[k]$ and $e[k]$, respectively. The SNR quantifies the average impact of the distortions introduced by quantization. The statistical properties of the signal $x[k]$ and the quantization error $e[k]$ are required in order to evaluate the SNR of a quantizer. First, a statistical model for the quantization error is introduced."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Model for the Quantization Error\n",
    "\n",
    "In order to derive the statistical properties of the quantization error, the probability density functions (PDFs) of the quantized signal $x_\\text{Q}[k]$ and the error $e[k]$, as well as its bivariate PDFs have to be derived. The underlying calculus is quite tedious due to the nonlinear nature of quantization. Please refer to [[Widrow](../index.ipynb#Literature)] for a detailed treatment. The resulting model is summarized in the following. We focus on the non-clipping case $x_\\text{min} \\leq x[k] < x_\\text{max}$ first, hence on granular distortions. Here the quantization error is in general bounded $|e[k]| < \\frac{Q}{2}$.\n",
    "\n",
    "Under the assumption that the input signal has a wide dynamic range compared to the quantization step size $Q$, the quantization error $e[k]$ can be approximated by the following statistical model\n",
    "\n",
    "1. The quantization error $e[k]$ is not correlated with the input signal $x[k]$\n",
    "\n",
    "2. The quantization error is [white](../random_signals/white_noise.ipynb)\n",
    "\n",
    "    $$ \\Phi_{ee}(\\mathrm{e}^{\\,\\mathrm{j}\\,\\Omega}) = \\sigma_e^2 $$\n",
    "\n",
    "3. The probability density function (PDF) of the quantization error is given by the zero-mean [uniform distribution](../random_signals/important_distributions.ipynb#Uniform-Distribution)\n",
    "\n",
    "    $$ p_e(\\theta) = \\frac{1}{Q} \\cdot \\text{rect} \\left( \\frac{\\theta}{Q} \\right) $$\n",
    "\n",
    "The variance of the quantization error is then [derived from its PDF](../random_signals/important_distributions.ipynb#Uniform-Distribution) as\n",
    "\n",
    "\\begin{equation}\n",
    "\\sigma_e^2 = \\frac{Q^2}{12}\n",
    "\\end{equation}\n",
    "\n",
    "Let's assume that the quantization index is represented as binary or [fixed-point number](https://en.wikipedia.org/wiki/Fixed-point_arithmetic) with $w$-bits. The common notation for the mid-tread quantizer is that $x_\\text{min}$ can be represented exactly. Half of the $2^w$ quantization indexes is used for the negative signal values, the other half for the positive ones including zero. The quantization step is then given as\n",
    "\n",
    "\\begin{equation}\n",
    "Q = \\frac{ |x_\\text{min}|}{2^{w-1}} = \\frac{ x_\\text{max}}{2^{w-1} - 1}\n",
    "\\end{equation}\n",
    "\n",
    "where $x_\\text{max} = |x_\\text{min}| - Q$. Introducing the quantization step, the variance of the quantization error can be expressed by the word length $w$ as\n",
    "\n",
    "\\begin{equation}\n",
    "\\sigma_e^2 = \\frac{x^2_\\text{max}}{3 \\cdot 2^{2w}}\n",
    "\\end{equation}\n",
    "\n",
    "The average power of the quantization error quarters per additional bit spend. Introducing the variance into the definition of the SNR yields\n",
    "\n",
    "\\begin{equation}\n",
    "\\begin{split}\n",
    "SNR &= 10 \\cdot \\log_{10} \\left( \\frac{3 \\sigma_x^2}{x^2_\\text{max}}  \\right) + 10 \\cdot \\log_{10} \\left( 2^{2w} \\right) \\\\\n",
    "& \\approx 10 \\cdot \\log_{10} \\left( \\frac{3 \\sigma_x^2}{x^2_\\text{max}}  \\right) + 6.02 w \\quad \\text{in dB}\n",
    "\\end{split}\n",
    "\\end{equation}\n",
    "\n",
    "It now can be concluded that the SNR decays approximately by 6 dB per additional bit spend. This is often referred to as the 6 dB per bit rule of thumb for linear uniform quantization. Note, this holds only under the assumptions stated above."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Uniformly Distributed Signal\n",
    "\n",
    "A statistical model for the input signal $x[k]$ is required in order to calculate the average SNR of a linear uniform quantizer. For a signal that conforms to a zero-mean uniform distribution and under the assumption $x_\\text{max} \\gg Q$ its PDF is given as\n",
    "\n",
    "\\begin{equation}\n",
    "p_x(\\theta) = \\frac{1}{2 x_\\text{max}} \\text{rect}\\left( \\frac{\\theta}{2 x_\\text{max}} \\right)\n",
    "\\end{equation}\n",
    "\n",
    "Hence, all amplitudes between $-x_\\text{max}$ and $x_\\text{max}$ occur with the same probability. The variance of the signal is then calculated to\n",
    "\n",
    "\\begin{equation}\n",
    "\\sigma_x^2 = \\frac{4 x_\\text{max}^2}{12}\n",
    "\\end{equation}\n",
    "\n",
    "Introducing $\\sigma_x^2$ and $\\sigma_e^2$ into the definition of the SNR yields\n",
    "\n",
    "\\begin{equation}\n",
    "SNR = 10 \\cdot \\log_{10} \\left( 2^{2 w} \\right) \\approx 6.02 \\, w \\quad \\text{in dB}\n",
    "\\end{equation}\n",
    "\n",
    "The word length $w$ and resulting SNRs for some typical digital signal representations are\n",
    "\n",
    "| | $w$ | SNR |\n",
    "|----|:----:|:----:|\n",
    "| Compact Disc (CD) | 16 bit | 96 dB |\n",
    "| Digital Video Disc (DVD) | 24 bit | 144 dB |\n",
    "| Video Signals | 8 bit | 48 dB |\n",
    "\n",
    "Note that the SNR values hold only if the continuous amplitude signal conforms reasonably well to a uniform PDF and if it uses the full amplitude range of the quantizer. If the latter is not the case this can be considered by introducing the level $0 < A \\leq 1$ into above considerations, such that $x_\\text{min} \\leq \\frac{x[k]}{A} < x_\\text{max}$. The resulting variance is given as\n",
    "\n",
    "\\begin{equation}\n",
    "\\sigma_x^2 = \\frac{4 x_\\text{max}^2 A^2}{12}\n",
    "\\end{equation}\n",
    "\n",
    "introduced into the definition of the SNR yields\n",
    "\n",
    "\\begin{equation}\n",
    "SNR = 10 \\cdot \\log_{10} \\left( 2^{2 w} \\right) + 20 \\cdot \\log_{10} ( A ) \\approx 6.02 \\, w + 20 \\cdot \\log_{10} ( A ) \\quad \\text{in dB}\n",
    "\\end{equation}\n",
    "\n",
    "From this it can be concluded that a level of -6 dB is equivalent to a loss of one bit in terms of SNR of the quantized signal."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Example - Quantization of a uniformly distributed signal \n",
    "\n",
    "In this example the linear uniform quantization of a random signal drawn from a uniform distribution is evaluated. The amplitude range of the quantizer is $x_\\text{min} = -1$ and $x_\\text{max} = 1 - Q$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "SNR = 48.090272 in dB\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAoAAAAEYCAYAAADMEEeQAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzs3Xd4VGX2wPHvSYcQICEJNfQOCkKk2bAja3ctWLG7lt1Vd1fdddVVd39u013bsq4iVuwFxYZYEUGK9C41hEBIgARCIOX8/rg3cDOZSWbSJuV8nmeezNz6zp07b868VVQVY4wxxhjTfESEOwHGGGOMMaZ+WQBojDHGGNPMWABojDHGGNPMWABojDHGGNPMWABojDHGGNPMWABojDHGGNPMWADoQ0SOE5HV4U6HPyIyVkQyKlmvItI7wLrLROSzuktdwyciLUTkAxHZIyJvhjs9lamrz0tEfi8iz9b2cY0JVWPOaxsSETlGRNaKyF4ROTfc6amMiHwsIlfVwXGXi8jY2j5uU9dkAkAR2Sgi+90vQdnjySD2Kxc0qeq3qtqvjtI4RUQerotjV0VVX1HV06raLpxprAc/B9oD7VT1wnAnpoyIdHfvw6iyZcF+XlUct8I/MVX9i6peV5PjmubN8tpD72Wf+963isijIhLprjtWRGa7PzRzReQ7ETnaXTdRREo8122DiDwvIn1rkJwHgSdVtZWqvlcb7682iMgDIvKyd5mqnqGqL9TwuBU+W1UdpKpf1eS4zVGTCQBdZ7lfgrLHreFOkDnMG+CESTdgjaoWhzkdzZq/+yDUe0McTS3/akwsr4UhqtoKOBm4FLheRFoDHwJPAElAZ+BPwAHPft+7+7UBTgH2AwtEZHA109ENWF7NfU0tqaV8rX7/R6pqk3gAG4FTAqzrDXwN7AF2Aq+7y78BFNgH7AUuBsYCGT7H/S2wxN3uOZxSpI+BfOBzINGz/ZtAlnuub4BB7vIbgCLgoHuuD9zlnYC3gWxgA/BLz7FaAFOAXcAKNx0ZlVwDBW4C1rr7PAWIu24iMMt9LsBjwA43nUuAwZWkcQDwFbAbJ6M523POdsAHQB4wD3i47DyeNN3ipmmDu+zfwBZ3nwXAcZ7tH3Cv4cvu9V0K9AXucdO7BTitkmvgN604mfBB9/3tBa71s2+l19t9L709r6cAD7vPE3Ey/mx3/w+BLp5tvwIeAr5z39dnQLK7brN77L3uY7TP5/U7z7q97nuY4q67GljpHnM9cKO7PB7nH0upZ79O7vV92ZOus93rtNtN4wCfe/83OPfHHuB1IK6Sa3+Nm5ZdwKdAtyruA3/LxuDcR3vcv2N8ruGf3Wu43/tZ2MPyWuo/r+3tk5YngXRgdyX7TcSTP3qWfwi8Vcl+1wPrgFxgGtDJXf4Tznd8v/teY/3sexSw0L2GrwOvcTjfqpAe73sDfgb8iJNXbwEe8GzX3d32Kpw8bCfwB3fdOMrnt4vd5V8B17nPF1M+X1NgbDU/24249yQQC/wLyHQf/yq7Lrj3HHAnzv+TbcDVlVz3Njj34TZgK87/t0jPtfsO539prrvO37II4F5gk3vOF4E2PtfwWvcaflOv3+VwZya19kYqz5SmAn9wP4g44Fh/N7v3BvE57hycjKiz+wEuxPlSxQJfAPd7tr8GSPDchIs866bgfvHc1xE4AdB9QAzQE+ef+Onu+keAb3F+SaYBy6g6U/oQaAt0xcnoxvl+0YHT3fO2xQkGBwAdA6QxGifj+b2bxpNwMpJ+7vrX3EdLYCBOJuEbAM5w30MLd9nlOIFjFM4XMQs3sMAJUArdNEa5X5YN7ucXjZMRbgjw/qtK6wN4gh8/+1d6vf3cK4eulft+LnCvQwJOBvaeZ9uvcDLrvjj/bL4CHnHXdXePHeXZ/tDn5ZPGNJxMbbz7+mdAL/dzPAEoAIb5u5d9r4Gbln3Aqe61+517/WI89/4POP84k3CCu5sCXLtz3X0HuJ/bvcDsKu6Dcsvcv7uAK9xjTHBft/Ncw83AIHd9dLjzneb4wPLacu8FJ9/Lwvkn3hrIAV4AzsATsLrbTsT/9/oaYHuAc52EE1wNc9/rE3gChSo+jxicwON2nO/4z3ECqGADwLHAEe71OxLYDpzrruvubvs/9/s7BKekc4C7/gF88ls8AaDP8huAVUDrUD9b32uAUyU+B0gFUoDZwEOe91PsbhMNjMfJMxMDXL/3gP/i/KBOxckPb/Rcu2LgNpz8qEWAZdfg5I09gVbAO8BLPtfwRfccLer1u1yfJ6vTN+LcAHtxSjLKHte7614EnsFTIuPvZvfcIL6Z0mWe128D//G8vg3PP3qfY7d1j18W7Ze7cYGRwGaffe4Bnnefr8cN4DxfkqoyJW+G+wZwt+dmLQsATwLWAKOACJ9j+KbxOJzMLcKzbCrOlzsSJzPp51nnrwTwpCo+u1041Sm4x53hWXeW+7mW/epKcI/Z1s9xAqbVc+zKAsBKr7efe6XctfI51lBgl+f1V8C9ntc3A5+4z7sTRACIk5ksAO6q5D28B/zK373sew2APwJveNZF4PzKHeu59y/3rP8bMCnAeT/GU6rqHqsAtxTQ333guwwn8PvBZ5vvgYmea/hgZfeSPer+geW1Ze8lDyfv+gm3pMddN8A9fwZOMDANaO+uq/C9dpePA4oCnOs54G+e161w8t3unusWKAA8HucHo3iWzSbIANDP8f4FPOY+7+5u663p+AG4xH3+AEEEgMCxOMF+3+p8tr7XwP08xnvWnQ5s9Nxz+ymf1+4ARvk5b3ucgLaFZ9kE4EvPtfO9p/wtmwnc7Hndz/38ojzXsGddfmcDPZpaG5pzVbWt5/E/d/nvcEpIfnB7C10T4nG3e57v9/O6FYCIRIrIIyLyk4jk4dyUAMkBjtsN6CQiu8seOKVX7d31nXBK1MpsCiKtWZ7nBWVp81LVL3CqK54CtovIM27bFX86AVtUtdQnHZ1xfl1F+aTR+9zvMhG5U0RWuo2kd+MUs3uvke/13amqJZ7X+HtfVaQ1GNW53gCISEsR+a+IbHI/+2+AtmUNw11VfjZVeA5Yrap/9Zz3DBGZ4zY2343zizbQ/earE5736F63LZS/XsGmuRvwb899nIvznfMeq6p7o1x6XL6fn79jmPpnea1T0p6oqr1U9d6yfEdVV6rqRFXtgtO0phNO4FSZzjjfGX98v6d7cUoZg8nXOgFb1Y08XKHkayNF5EsRyRaRPThNjHyvcbXzNRFJwymouEpV17jLQv1sffnmI5vcZWVytHw78EBp7oZTSrjNc8/8F6cksEyV/+8CpCeKw/deoOPUuaYWAPqlqlmqer2qdgJuBJ6WAMOl1NClwDk4DXvb4ET34GSI4ET6XltwqjO9GWmCqo5312/DqY4o07W2Eqqqj6vqcJzqtL44bV78pTETSPNpcN8Vp6QoG+cXbhfPOm96D52u7ImIHAfcBVyEU+zeFqedh/jZL1SVpTUYVV3vApwq3jIdPM/vxPllN1JVW+P88obg3pfvNa9ARO52j3+tZ1ksTinJP3BKGNoCHxH4fvOViZPJlR1PcN5/sNfLawtO1Yj3Xm6hqrM92/hLj3dZufS4fD+/Kq+VCR/La8tT1VU4JVZVdfA4D6cK2h/f72k8TpOTYL6n24DO7ne7jPe97cOTp4mIN08DeBWnBDNNVdsAkwg+r670uyoiLXBqLP6lqh97VoX62fryzUe6ustCtQWnBDDZc8+0VtVBnm2qytMCpaeY8j9uwpKvNYsAUEQuFJGyIGUXzsUuK1HajlM3XxsScG6YHJwv1V981vue6wcgT0TuEmeMukgRGVw2ZADOL6N7RCTRTf9ttZFIETna/WUXjZMBFBL4esx1t/mdiESLM9bSWcBrbqncO8ADbglYf+DKKk6fgHPzZwNRInIfTruZ2hAwrUHuX9X1XgRc6n5O43Da3JVJwCmh2C0iScD9IaQ7G6cht9/7UETOAH6JU+qy37MqBqeNTDZQ7G7nHTpmO9BORNoEOO8bwM9E5GT3XrgT5/6dHWD7ykzCuXaD3DS3EZFQh9r5COgrIpeKSJSIXIzTvurDaqTHhEFzz2tFpL9bw9HFfZ2GU204x8+2kSLSQ0SewKma/FOAw74KXC0iQ90ffX8B5qrqxiCS9D1OfvtL9zt1PjDCs34xMMg9dhxOta1XApCrqoUiMgInOAvWdqB7Jb31JwOrVPVvfs4Zymfraypwr4ikiEgyTrvPlyvZ3i9V3YbTWe+fItJaRCJEpJeInFDVvn7Sc7v7WbfCeT+vawMYjaKpBYAfSPmxqd51lx8NzBWRvTi/Zn6lqhvcdQ8AL7hFvBfV8Pwv4hTvbsXpSeb7pX8OGOie6z03gDoLp73YBpyGvs/i/OoBJ0PY5K77DHiphukr0xqn4e4u9/g5OKVI/tJ4EKen6Blu+p4GrnR/2QLc6qY3y03fVMoPeeDrU5z2YmvccxdSS8XfQaS1KlVd71/hfF67gctwfr2W+RdOG72dOJ/7JyGkuwC3d6t73Uf5bHIxTnX7Ss+9PUlV83ECwzdwPstLce7vsuOuwvk81rvH9VaDoKqrcTrkPOGm+yyc4T0OBpt2z7HeBf4KvOZW2yzD+RxCOUYOcCZOIJqDU514pqruDDU9ps5ZXutfPk57w7kiss9N1zKce7rMaPf65OG0iWsNHK2qS/0dUFVn4rTXfRunRK8XcEkwiXG/y+fjtE3bhZOXvONZvwanQ8TnOL3xZ/kc4mbgQRHJxwmk3gjmvK6ywfZzRGShn/WXAOf53EfHEeJn6+e4DwPzcXqTL8XpSFTdMSGvxPmhvQLn+r0FdAzxGJNx7qdvcO6vQmqpMKemyoYIMaZWiMhfgQ6qelW401JTbgniy25bHmOMafREZApOB5d7w50WE15NrQTQ1DO3yuNIcYzAaaP2blX7GWOMMSZ8whIAishkEdkhIssCrL9MRJa4j9kiMqS+02iCloBTpbAPp3rgn8D7YU2RMQ2ciKSJ07typTi9ZX/lZxsRkcdFZJ2bFw7zrLtKnPlf10odzK1qjGn6wlIFLCLH44wj9aKqVugdJSJjgJWqustt2P6Aqo6s73QaY0xdEJGOOIOvLxSRBJzxHc9V1RWebcbjtBUaj9Ou7N+qOtLtZDQfZ9YJdfcdrqq76vt9GGMar7CUAKrqNwQe8whVne3JzOZQfpgRY4xp1FR1m6oudJ/n48yy4juu2zk4P5JVVefgjCvZEWdg2xmqmuvmkzNwBhI2xpig1e/Ew9VzLU6v0QpE5AacEduJj48f3r9///pMlzEmTBYsWLBTVVPCnY7aICLdcaY7m+uzqjPle8hnuMsCLfc9ruWPxjRDweaPDToAFJETcQLAY/2tV9VncKYdIj09XefPn1+PqTPGhIuIBD2bQUPmjgv2NvBrVc3zXe1nF61kefkFlj8a0ywFmz822F7AInIkzjhN57jjgxljTJMhzuDbbwOvqOo7fjbJoPzsFF1wZhUItNwYY4LWIANAEemK07P0irL5AY0xpqkQEcEZ0Halqj4aYLNpwJVub+BRwB53doJPgdPcWSsScWZ/+bReEm6MaTLCUgUsIlNxpr5JFpEMnGmzogFUdRLOiOPtcOaRBChW1fRwpNUYY+rAMcAVwFIRWeQu+z3uPK1uPvgRTg/gdTjzUF/trssVkYeAee5+D6pqwE51xhjjT1gCQFWdUMX664Dr6ik5xhhTr1R1Fv7b8nm3UeCWAOsm40wxZYwx1dIgq4CNMcYYY0zdsQDQGGOMMaaZsQDQGGOMMaaZsQDQGGOMMaaZsQDQGGOMMaaZsQDQGGOMMaaZsQDQGGOMMaaZsQDQGGOMMaaZsQDQGGOMMaaZsQDQGGOMMaaZsQDQGGOMMaaZsQDQGGOMMaaZsQDQGGOMMaaZsQDQGGOMMaaZsQDQGGOMMaaZsQDQGGOMMaaZsQDQGGOMMaaZsQDQGGOMMaaZiQp3AowxprkRkcnAmcAOVR3sZ/1vgcvcl1HAACBFVXNFZCOQD5QAxaqaXj+pNsY0JVYCaIwx9W8KMC7QSlX9u6oOVdWhwD3A16qa69nkRHe9BX/GmGqxANAYY+qZqn4D5Fa5oWMCMLUOk2OMaYYsADTGmAZKRFrilBS+7VmswGciskBEbqhk3xtEZL6IzM/Ozq7rpBpjGhlrA2jK6X739JD32fjIz+ogJcYY4CzgO5/q32NUNVNEUoEZIrLKLVEsR1WfAZ4BSE9P1/pJrjGmsbAAsBGx4MyYZucSfKp/VTXT/btDRN4FRgAVAkBjjKmMVQEbY0wDJCJtgBOA9z3L4kUkoew5cBqwLDwpNMY0ZlYCaIwx9UxEpgJjgWQRyQDuB6IBVHWSu9l5wGequs+za3vgXREBJ/9+VVU/qa90G2OaDgsAjTGmnqnqhCC2mYIzXIx32XpgSN2kyhjTnIQlAAxiEFQB/g2MBwqAiaq6sH5TaRoSa/9ojDHG1J5wlQBOAZ4EXgyw/gygj/sYCfzH/WuagOoEc8YYY4ypPWEJAFX1GxHpXskm5wAvqqoCc0SkrYh0VNVt9ZJAY4IUajBrpZLGGGMagobaBrAzsMXzOsNdVi4AdAdBvQGga9eu9ZY4Y6rLqrKNMcY0BA11GBjxs6zCQKaq+oyqpqtqekpKSj0kyxhjjDGm8WuoJYAZQJrndRcgM0xpMVWwNn3GGGNM49JQSwCnAVeKYxSwx9r/GWOMMcbUjnANA1PVIKgf4QwBsw5nGJirw5FOY4wxxpimKFy9gCsdBNXt/XtLPSXHGGOMMaZZaahVwMYYY4wxpo401E4gxtSYjdFnjDHG+GcBYBNnPXSDZ9fKGGNMc2FVwMYYY4wxzYwFgMYYY4wxzYwFgMYYY4wxzYwFgMYYY4wxzYwFgMYYY4wxzYwFgMYYY4wxzYwFgMYYU89EZLKI7BCRZQHWjxWRPSKyyH3c51k3TkRWi8g6Ebm7/lJtjGlKLAA0xpj6NwUYV8U236rqUPfxIICIRAJPAWcAA4EJIjKwTlNqjGmSLAA0xph6pqrfALnV2HUEsE5V16vqQeA14JxaTZwxplmwANAYYxqm0SKyWEQ+FpFB7rLOwBbPNhnuMmOMCYlNBWeMMQ3PQqCbqu4VkfHAe0AfQPxsq/4OICI3ADcAdO3ata7SaYxppKwE0BhjGhhVzVPVve7zj4BoEUnGKfFL82zaBcgMcIxnVDVdVdNTUlLqPM3GmMbFAkBjjGlgRKSDiIj7fAROXp0DzAP6iEgPEYkBLgGmhS+lxpjGyqqAjTGmnonIVGAskCwiGcD9QDSAqk4Cfg78QkSKgf3AJaqqQLGI3Ap8CkQCk1V1eRjegjGmkbMA0Bhj6pmqTqhi/ZPAkwHWfQR8VBfpMsY0H1YFbIwxxhjTzFgAaIwxxhjTzFgAaIwxxhjTzFgAaIwxxhjTzFgAaIwxxhjTzFgAaIwxxhjTzFgAaIwxxhjTzFgAaIwxxhjTzFgAaIwxxhjTzIQlABSRcSKyWkTWicjdftZ3FZEvReRHEVkiIuPDkU5jjDHGmKao3gNAEYkEngLOAAYCE0RkoM9m9wJvqOpROJOdP12/qTTGGGOMabrCUQI4AlinqutV9SDwGnCOzzYKtHaftwEy6zF9xhhjjDFNWjgCwM7AFs/rDHeZ1wPA5SKSgTPp+W31kzRjjDGm6XlpziaueG5uuJNhGpBwBIDiZ5n6vJ4ATFHVLsB44CURqZBWEblBROaLyPzs7Ow6SKoxxhjTuO3IL+SRj1by7dqdFBaVhDs5poEIRwCYAaR5XnehYhXvtcAbAKr6PRAHJPseSFWfUdV0VU1PSUmpo+QaY4wxjddjM9aw76AT+GXnHwhzakxDEY4AcB7QR0R6iEgMTiePaT7bbAZOBhCRATgBoBXxGWOMMSFYlZXH6/O20L9DAgA7LAA0rnoPAFW1GLgV+BRYidPbd7mIPCgiZ7ub3QlcLyKLganARFX1rSY2xhhjTACqyp+nryQhLpr7znIG28jOLwxzqkxDERWOk6rqRzidO7zL7vM8XwEcU9/pMsYYY5qKr9Zk8+3anfzxzIH0SbUSQFOezQRijDH1TEQmi8gOEVkWYP1l7iD4S0RktogM8azbKCJLRWSRiMyvv1SbxqS4pJQ/T19J93YtuWJUN9rFxxAZIWzPsxJA47AA0Bhj6t8UYFwl6zcAJ6jqkcBDwDM+609U1aGqml5H6TON3LTFmazbsZd7xg8gJiqCiAghuVUMO/KsBNA4wlIFbIwxzZmqfiMi3StZP9vzcg7OaAnGBG1Jxh5axUZx2sD2h5a1bx1nVcDmECsBNMaYhu1a4GPPawU+E5EFInJDmNJkGritu/fTJbEFIoeH3k1NiLUA0BxiAaAxxjRQInIiTgB4l2fxMao6DGc+9VtE5PgA+9pA+c1Yxq79dG7botyylIS4BtsL+KfsveFOQrNjAaAxxjRAInIk8CxwjqrmlC1X1Uz37w7gXZz51SuwgfKbt627CuiSWD4ATE2IJWffQYpLSsOUKv+WZOzm5H9+zUdLt4U7Kc2KBYDGGFMJETlVRP4nIkPd13Ve7SoiXYF3gCtUdY1nebyIJJQ9B04D/PYkNs1XXmEReYXFdPYNAFvHogo79x4MU8r8W7xlNwAvfr8xrOlobqwTiDHGVO5m4GrgXhFJAoZ6V4rIU0C8qk4UkdNU9bOqDigiU4GxQLKIZAD3A9EAqjoJuA9oBzzttuEqdnv8tgfedZdFAa+q6ie18i5Nk7F1134AuiS2LLc8NSEOcOYG7tAmrt7TFciqrHwA5qzPZd2OvfRObRXmFDUPFgAaY0zlslV1N/AbEXkEONpn/UFgu/v8JKDKAFBVJ1Sx/jrgOj/L1wNDKu5hzGEZbgDo2wYwNSEWoMENBbM6K5/eqa3YlLOPV+Zu4v6zBoU7Sc2CVQEbY0zlppc9UdW7gRd91hcAbUQkGuhanwkzxp+tuwoAKrYBbO0EgNsbUEcQVWX19nxG9Uxi3OCOvL0gg/0HS8KdrGbBAkBjjKmEqr5f9lxE7lfVJ3w2uR/4CXgKeLU+02YaloYyZX3Grv3ERUeQFB9Tbnlyq1hEGlYJYOaeQvILi+nXoTWXjexKXmExHy7JDHeymgWrAjbGmODdLyItgSRgIfCaqu4Cng5vsky4lZQq5z/9HWN6J3PXuP5hTYszBmDLcmMAAkRHRtAuPqZBjQW4OisPgP4dEkjvlkivlHhembuZC9PTwpyyps9KAI0xJngKFAKfAmlAuXl6TfP1/qKtLM7Y0yCGMvE3BmCZhjYW4OosZ/y/vqkJiAiXjezGoi27WbZ1T5hT1vRZAGiMMcFbpar3q+pbqvp74BzgsXAnyoRXUUkp/565lgiBTTkFZO7eH9b0lM0C4k9Dmw1kdVYeHdvE0aZlNAAXDOtCXHQEr/6wOcwpa/qCDgBFpJOIXCEiN4rICXWZKGOMaaB2isjwshfuGH02ynIz987CDDblFHD7KX0BmLM+p4o96k7BwWJy9x2sMAZgmdSE2AbVBnBVVj79OiQcet2mZTRnHdmJ937cSn5hURhT1vQFFQCKyGnAAmA8MBr4l4isFpExdZk4Y4xpYH4JvCwiL4vIXSLyCrAh3Iky4XOwuJTHZ65jSFpbbj6xN21bRvP9T+ELAAONAVgmtXUsO/ceoLQ0/B1WikpK+Sl7b7kAEODSkV0pOFjC+4usM0hdCrYE8GHgOFWdoKoTVfUoYCLwjIiMFJHBdZZCY4xpIFR1Mc5A0FPdRV8ClY7pZ5q2N+ZvYevu/dxxal8iI4SRPZKYsyF8AWCgMQDLpCbEUVyq5BaEfzaQDTv3UVSi9PcJAIemtWVAx9a8Mndzg+lZ3RQFGwDGqOo67wJV/R44H3iJw5mhMcY0OSIyWtwulap6QFWnq+pfVfVZVd0X7vSZ8CgsKuHJL9aR3i2R4/skAzC6Zzu25O4nwx2Lr75luO0P0yqpAoaGMRRM2Qwg/dq3LrdcRLh0ZFdWbstjkTtNnKl9wQaAhSJSoZ2L2/6lBDi5VlNljDENy1XAAhF5TUQmikiHcCfIhN/LczaRlVfIHaf1PTTkyqhe7QDCVg2csauAmMgIklvF+l3fkAaDXp2VR2SE0Cs1vsK6c4d2omVMJK/Otc4gdSXYAPDvwHsi0sm7UESSgQOquqPWU2aMMQ2Eqt6kqsOAB4BEYIqIfC8ifxGR40UkMrwpNPVpRWYe170wj4enr+S4PsmM6ZV8aF3f1ASS4mOYsz43LGnbums/nRNbEBEhfteXzQec3QBKAFdn5dMzOZ7YqIpfn4S4aM4e0okPlmSyZ39onUEmz9rA2U/OoqiktLaS2iQFNRC0qr4tIrHA9yKyAFgMxAAX4bQPNMaYJk9VVwGrROQZIBI4FrgQeBRID2faTN3blLOPv326mulLttE6Lorfnt6PiWO6l9smIkIY1TOJOetzUNUKgzHXtcrGAARIKasCbgAlgKuy8hma1jbg+stGduO1eVt478etXOVznSvzxvwtrMrKZ/qSbZx7VOdaSGnTFPQwMKr6KjAA+BBoAxQBl6rqC3WUNmOMaRBEJEJELhWR6SKyHVgFrMOpHSnEOoI0eRt37uOC/8zmy1U7uPXE3nx710nccmJv4mMrlqOM6tmOrbv3syU38HiAe/YXMeGZOby/aGutprOyMQAB4qIjadMiOuxjAe49UEzGrv0VOoB4HdGlDUd0bsOrIXQGydhVcKht4X+/WV/vnUgOFJdw11tLWO2moSELdhiY0SIiqlqgqpNV9Q5VfUBV59V1Ao0xpgH4EugF3AN0VNU0VU0FjgPmAP8nIpeHM4EmdHsPFPNT9t4qt9uRV8gVk+dSqjDt1mP5zen9aNMiOuD2o3s67QArGw/wqS/X8f36HH7z5mJ+2FA71cWFRSVk5x+otAQQGsZYgGu2OwFS3/aBA0BwhoRZvT2fhZt3BXXcmSudFmm/GNuLldvy+HbtzpolNETfrtnJ6/O38NCHK+r1vNURbAmgNYA2xjRnp6jqQ6q6RFUPNSxS1VxVfVtVfw68Hsb01ZkNO/fx/qKt1S5J+XL1Dla58702JGu25/Ozx7/l9Me+Yd7GwAFYXmERVz0/j5y9B3l+4tH0Tm3gpSdCAAAgAElEQVRV5bF7p7YiuVUM3wcIADfl7OP57zbwsyM7kpbYkhtfms+mnJp3Ji+bgaRLUhUBYOvYsFcBl5WQ9e/QutLtzh7SiVaxUbwSZGeQmat20DM5nl+f0ofUhFie+WZ9jdMaik+WZwEwa91O5oZxQPBgBBUAWgNoY0xzpqqHWqGLiN8iC+82Tcl97y/jV68t4pZXF7LvQHFI+076+ieufn4et736Y4Maz+3T5Vmc99R3FBwsoWPbOG5+ZSE78ioGRIVFJVz/wnzW7chn0uXDGVJJezUvEWFkz3Z8/1OO3/f9yMeriI6M4P4zB/LcxKNR4Jop80Lu7ODr8BiA/geBLpOaEBf2KuDVWfm0jImstLoaID42inOGdmL6km3s2lf52IV7DxQz56ccTh6QSmxUJNcc24NZ63bW27zCRSWlfL5yO2cM7kBKQiyPzlhTq/f9jrxC3pi3pdZ6Roc0F7CqrlLVx1R1HHASMAunAfTcWkmNMcY0MCIyRkSmiciLIvIyzWjc02179jNr3U6GpLXlk2VZXPCf2WzOqXp8O1Xln5+t5pGPV9E7tRVrd+xl1rr6rYrzp7RUeXTGGm58aQG9U1vxwa3H8uyVR7O3sJibX1nIweLDvUYzdhVw1eQfmLshl39cOITj+4Y249/onu3Iyitkk8/1mrs+h4+XZfGLE3qR2jqOHsnxTLp8OJtzC7j11YU16rm61S0BDDQNXJmy+YDDGZSvysqjb/uEgL2Vva4Y3Y0DxaW8PGdTpdvNWpvNwZJSTurfHnCqj1vFRtVbKeAPG3LZXVDEuUd15paxvZi7IbfGwwHtPVDMozPWcOYT3zLiLzP53dtL+P27S9m4s+YlxiEFgGVEJB44qKofqeptqmq934wxTdXNwHmqeiXwa6BLmNNTb979cSuq8PglQ5ly9Qi27SnkrCdnMauSdlWqyoMfruCJL9ZxcXoaH9x6LMmtYnj+u411ksZQAqYpszfy+My1XDCsC6/fOJoObeLo1yGBv/38SOZv2sVfPlqJqvLm/C2M+9e3LNu6h8cuHsI5Q0PvSTrKbQc4fem2Q9OulZYqD09fScc2cVx3XM9y2/75vCP4du3OGpXuZOwqICpCaJ/gfwzAMikJsRwsLq1xiWN1qSqrs/Ir7QDi1b9Da8b2S2HK7I0UFpUE3G7myh20josivXsiAK3jorl0ZFemL93Glty6H5j7k2VZtIiO5Pg+KVwyoisd28TVuBTw0c/W8MQXa4mLiuS3p/fjhWtGAPDhkppPkxdsJxBvD7gdwGogS0SWi8jfRaRPKCcVkXHuXMLrROTuANtcJCIr3HO8GsrxjTGmFuWragmAqu4EflXTA4rIZBHZISLLAqwXEXnczSOXiMgwz7qrRGSt+7iqpmkJRFV5e0EGR3dPpFu7eI7vm8K0W4+hQ+s4Jj7/g9/eq4VFJdz55mKe/24jVx/TnUcuOIIWMZFcOrIbX6zawYYqSi32HSjmpH98FXTP2Kk/bGbInz5jSUbVs0WoKq/N28zQtLb848IjiYs+3HLprCGduPbYHkyZvZHznp7Nb99awsCOrfnk18dz3lHVi/d7pcQzoGNr/v7pak7851c8/dU6Jn+3gaVb9/C7cf1oEVO+5dRF6WkMTWvLC99vrHbAsHXXfjq0iSMqsvJ/7amtnbEAw1UNvD3vALsKiirMAVyZm07oRc6+g7y1IMPv+tJS5cvVOxjbL5Voz/u/+pjuCPD4zLUBfyxszilg596aXYvSUuXT5VmM7ZdCi5hI4qIjueXE3szftItvqtkRpbiklGmLt3L6wA689Ysx3HJib07om0J6t0Q+WLytRumF4EsAvT3gOqhqF1VN4XAPuEeC7QHnthd8CjgDGAhMEJGBPtv0cc91jKoOwvnVbYwx4SAicrKIJIvIzcBrtXDMKcC4StafAfRxHzcA/3ETkgTcD4wERgD3i0hiLaSngsUZe/gpex8XDDscAHVrF89bvxjNsG6J/Pr1RbzkqZLbklvAhZO+552FW/n1KX2478yBh8bAu3xUV6IjhRdmb6z0nDNWbGf9zn089eW6KoOgDxZn8vt3l1JwsIR3FlYdMC7PzGPN9r1cMLyL37H57j6jP6N6JrEiM48/jB/A1BtGkZZUeVu6yogI7948hscuHkL71nH87ZPVPDx9JUd2acM5Q/yXKF41phvrs/dVu7o8Y1flQ8CUCfd0cO+5Af6xvZOr2PKwkT2SGJLWlv99u56S0or3xqKM3ezce5CTB6SWW96xTQsuPjqNNxdkcOxfv+DxmWvZkV/IvgPFvDF/CxdOms3xf/+Sa6fUbFCTH7fsZkf+AcYNPtxH9qL0NDq3bVHtUsBv1+5k596DnDes/P1y1pBOrN6eX+OhZoINAKvqAXcBwfeAGwGsU9X1qnoQJzM9x2eb64GnVHWXex6bacQYEy6/Ao4A/gBsAR6q6QFV9RugsrE/zgFeVMccoK2IdAROB2a4ee8uYAaVB5LV9vaCDGKjIhh/ZMdyyxPionnxmhGc1C+VP763jCe/WMs3a7I568lZbNy5j/9dmc6vT+lbLshKTYjjzCM78eb8LeQVBq52fG/RVkRgzfa9zK1kaJQvV+3g9tcXcXS3JI7vm8Iny7IOVbMG8u6PW4mOFM7yeT9loiMjeOGaEcy6+0SuP74nkUG0TatKXHQk5x3VhTduHM0Xd57Anaf25dGLhgRs9zb+iI60i4/hhdmVt3ULZOvu/VV2AAFof6gE8HDHl/d+3MrkWRuqdd5QlJYqr87dzIgeSfSpYggYLxHhpuN7simngE/dnrZeM1duJzJCGNs3tcK6B88ZzHNXpdOvQ2senbGGYx75ghF//pzfvbWEnL0HOWVAexZn7GF5ZvU7i3y6PIvoSOHE/ofPHxMVwS9P7s3iLbv5fGXoYczbCzNo2zKaE/uVf0/jj+hIhNS8GjjYALCziPxNRN4RkWdF5FYR6ebdIIQecJ1xMtEyGe4yr75AXxH5TkTmiEidZHDGGBOE3wPtcTq9LVHVp+vhnIHyyWDyzxo7UFzCtMWZnD6oA63jKo53FxcdyaQrhnPeUZ35x2druHLyD7RPiGPabcdy6sD2fo95zTE92HewhDfn+6/Cy9l7gG/X7uSq0d1p0yKal773HwTNXZ/DTS8voH/HBJ6dmM55R3UiK6+QRZVUAxeXlPL+okxO6p9K25YxAbeLjYo8NFVabeuZ0orbTu5D79TAQU9sVCQTRnRl5qrtIbdZO1hcSlZeYWglgG4V8Ferd3DHG4t48MMVlQ6JUxu+XbeTzbkFXD6qW9Ub+zhtUAd6JMfz369/qlCiNnPlDtK7JdKmZcX7NTJCOHlAe168ZgRf3HkCV47uztlDO/PWTaOZeecJ/OPCI4mJiuCNeVsq7BsMVeWTZVkc0zu5wvfl/GFd6Jkcz98+WUVxCO1V8wqL+GzFds4e0omYqPKhWkpCLGN6JfPB4swatS8MNgB8H6fd31PAqcAQ4BsRecqdIi4U/n76+L6DKJyqj7E4I+w/KyIV+t+LyA0iMl9E5mdnZ4eYDGOMqZqq3o8z1dt+YKKITBGRP9TxaQPlk8Hkn84BapA/zly5gz37i7hgeOD2b9GREfzzwiHcemJvLhvZlXdvGUOP5PiA2x/RpQ3p3RJ5YfZGv1V405duo6RUuWREGheld+HT5Vlk7Sk/NMu6HXu57oX5dElswQtXj6B1XDQnD2hPdKTw8dLAbaK+XbeTnXsPcP6wht9/57JRXYkQqbLHq6+sPYWoVt0DGJyhVeJjItmRd4ANO/dx29Qf6ds+gc5tW3Dvu8vqdA7dl+dsol18DKcP8v9DoTKREcL1x/VkccaecmMsls3+ccqAqo/ZM6UVfzxzIP93/hGkd09CRGjbMoZxgzrw7o9bK+1kEsjKbflszi1g3KCKQyRHR0bwu3H9WLtjb8D2i/58vHQbB4tLOS/AVHZnDenIxpwClm2t/hibwQaAkar6nKrOBHJV9XqcNoEbgWdCPGcGkOZ53QXwLcfMAN5X1SJV3YATfFboaKKqz6hquqqmp6SE1kXfGGOCparZ7qgHf1LViUBw0xJUX6B8Mpj8E6hZ/vj2ggzat46tso1WRITwm9P78efzjqBlTNVTy199TA825xbwxaqK1WHvL8qkf4cE+ndozeWjulGiyqs/HO4RW1hUwq2vLiQ6KoKXrh1Ju1ZO2UPruGiO65PCR0uzApaGvLNwq9+qtIaoY5sWnD6oPa/N28L+g8EHIxm7nBLDYEoAwekIsmHnXq5/cT5REcL/rkzn/rMGsnp7fpVtNasrc/d+Zq7czkVHpxEbVb3hg88f1pnkVrH8/dPVPD5zLXe8vohr3PZ7vu3/QnHJ0WnkFRbzybKK1ctV+WR5FhECpwQo/T59UAeGdW3LY5+voeBgcGNpvr1wKz2T4wPOlXz6oA5ERwof1KAaONgA8HMRudV9rgCqWqyqfwdGh3jOeUAfEekhIjHAJcA0n23eA04EEJFknCrh+h3O2xhjfIhIHxF5DhhUx6eaBlzp9gYeBexR1W3Ap8BpIpLodv44zV1Wa7LzD/DVmmzOPapzrbSD8zp9UHvSklrwfx+tLPePcEtuAQs27eLsoZ0Ap7PJ2L4pTP1h86Gx+f70wQpWZeXz6EVD6OQz1dkZgzuwdfd+lvoZ8DevsIjPlmdx1pEVq9IaqitHd2fP/iI+WHz4n/v+gyWVDmj8k9vDuksQbQDBqUb8cnU2G3bu46nLhpGW1JJTB7bn5P6pPDZjDdv2BJ7HuLpe+2EzClw6omu1jxEXHcn1x/Xgx827eXTGGuasz6FdfCy/PLkPPVOqnqUlkFE925GW1ILXQ6wGdqp/t3F09ySSW/mvEBURfj9+ANvzDgTVznJLbgE/bMjl/GGd/XZYAmjbMobj+6Tw4eLMKtu/BhLst+EOoI2IzAc6uVULl4vIU0BIoxyqajFwK06mtRJ4Q1WXi8iDInK2u9mnQI6IrMDpgfxbVW3Yc6oYY5qDl4C3cEZAQEQGi8iLoR5ERKYC3wP9RCRDRK4VkZtE5CZ3k49wfvSuA/6HMxYhqpqL0wllnvt40F1Wa6YtzqSkVPl5HVSXRkVG8NcLjmRDzj4enr6y3DnBmfarzJWju5Odf4BPl2cxbXEmU3/YzE0n9GKsn1K8Uwe2JypC+GhpxdKbT5ZmcaC4lPOH1XpTyTozskcS/donMGX2Rj5dnsVtU39k+MMzOPOJWXyyrGJVt6ry2g+b6ZUSH3wJoNsO8N6fDWBML6ekV0R44OxBFJcqD3+4srLdK1iRmVdpu8WiklJem7eFE/qm1Kh3NcD1x/Xkq9+MZdVD45h9z8lMvWEUd5zat0bHjIgQLk5P4/v1OSFNy/fPz9awZvveKseKTO+exKkD2zPp6/XkVDHkzHs/Or2kqzrmWUM6kbmnMOh5kn0FOxVcqar+GTgeZ0iCDsBwYBnOcAUhcatS+qpqL/e4qOp9qjrNfa6qeoeqDlTVI1S1NoZdMMaYmopQ1Y+BsnEBlwGDQz2Iqk5Q1Y6qGu0Oq/Wcqk5S1UnuelXVW9w88ghVne/Zd7Kq9nYfz9fWGyvzweJMBnVqHVIPzVCM6ZXMDcf15NW5m/l8xXZUlfd+3MrR3RPpkng4MDihbwpdk1ry5BfruOftJQzvlsidp/n/J9+2ZQxjeifz8bJtFaqB316YUWlVWkMkIlw1pjsrtuVx40sLmLU2m3OGdqZbu5Y8PrPiEDnfr89heWYe1x3XM6iZNQCuGNWNP4wfwMQx3cstT0tqyW0n9Wb60m18vSa4tqM79x7g3Ke/4/i/f8l1L8xn1tqdFdL4+Yrt7Mg/wOUjQ+/84SsiQuieHF9uLMfa8PPhaUQIvDE/uFLASV//xJNfrmPCiDQmjEircvu7xvVnf1EJT3yxLuA2qso7P25lZI+kKgPlUwa2JzYqolxJcSiCHQha3IQVqOo0VX1QVW9X1f+o6m7vNsYY04RlikgP3KYwbr4XXJFLI7Alt4BFW3Zz5pGdqt64Bu44rS8DO7bmrreX8O3anazdUbEEJSJCuGJUN1ZvzycqMoLHJxxVboBfX+MHd2BTTgErth1uFL85p4C5G3I576jAVWkN1QXDO3PXuP68eM0IfvjDKfzf+Udw20l9WLEtj5k+Q4o8++0G2sXHBOww4M/Inu24/viefq/L9cf3pGdyPPe9vyyoThFT5zpV9VeO6sbCzbu4/Lm5nPLo19zzzlKemLmWN+dv4X/frqdz2xblhklpaDq0iWNsv1TenJ9RZY/dl+ds4pGPV3HWkE48fO4RQd1fvVNbcVF6Gi/P2RSwlPHHLbvZsHNfUCXWrWKjOLFfKjNWbK9yW3+CHghaRG4TkXIV9yISIyInicgLQJ2NSG+MMQ3Er3GqZDuIyNU445j6nc2jMZru9qQ9M8BYebUlNiqSf18ylL0Hig91Qhh/RMVzXpSexnF9knl8wlF0blt5nH3aoA5ERggfL82iuKSUl77fyNlPzSImKqLCQLqNQWxUJL8Y24vj+6YcCnzPGdqJtKQWPPHF2kMlbOt25PPFqh1cObp7rZWIxUZF8vB5g9mUU8DjM9dWum1RSSkvz93EcX2S+dM5g5l990n888IhtIuP5bPlWfxzxhp++9YSFm7ezaUju9Z6u9LadvHRaezIP8BXqwOXfr77YwZ/fH8ZJ/dP5dGLhoT0nn59Sh8iIoT/fPWT3/X/+2Y9CbFRnOHn++BP79RWbK/mvM5Vd9tyjAOuAaa6v3534/zqjQA+Ax5T1UUhn90YYxoRVd3ojkt6Ls5wWF8Dk8Obqtrz4ZJMhqS1rXEbrWD0aZ/APWf054EPVnBy/1SS4iuOz9emZTQvXTsyqOMlxccwqmcS7yzMYMaK7azens/onu2476yB5aqWG7PoyAhuHtube95Zyjdrd3JC3xSem7WB2KgILh9V/Y4V/ozplczPh3fhmW/Wc/bQTvTv0Nrvdp8sy2J73gH+7/wjAKeTxgXDuxwaQqiwqISsPYXk7DvAkV0afjX8Sf1TSUmI5fEv1nJc3+QKvZUXbMrlt28uYVSPdjx12bBKS6X9ad86jovT03ht3mZ+eXKfch2aVmfl8/GyLG49sbff8Tf9SYyPoaRUydtf7HcMxMoE2wawUFWfVtVjgG7AycBRQC9Vvd6CP2NMU+Zt4uKOgPCWqv7RzRcLfbdpjDbu3MeyrXkBZ8qoC1eN6c5vTuvL7TVswF9m/BEdydxTyP6iEiZdPpxXrx/JgI7+A5fG6oJhXejUJo4nZq5l594DvL1wKxcM73JoWJza9IfxA2jdIprfv7M0YE/TF2ZvpFu7ln5n4AAnIOyeHM/wbkkhB0vhEB0ZwUPnDGZJxh7uf395uXXZ+Qe4+ZWFdGrbgklXDK92ieuNJ/REFZ75pvzgJk9+uY6WMZFce2yPoI+VFO8EfbkFB0NOR8ifhjs23za37d8kEWkBICLHh3x2Y4xpHJp8M5iyaaX8VcXWFRHh1pP6MLhzm1o53sXpaUyemM5ntx/PuMEdGl27v2DEREVw09hezN+0i9tfX8TB4tKQAoZQJMbH8IfxA1i4eTeveMZkLLNs6x7mb9rFFaO6Bd35pDEYN7gDt5zYi9fmbeHVuc77Li4p5bapC9ldUMR/Lh9GmxahlbZ5dUlsyfnDOjP1h81ku7OxrNuxlw+XZHLl6O4k+ikNDyTRndkmd189BIA+7gcmi8hLwNE1PJYxxjRU43B6/k4VkUwRWSEiG4C1OLMVPaaqU8KZwJr6cMk20rslVhhjrzGJiozgpP7ta713aENzUXoaqQmxfLt2J6cMSKVXDca/q8r5wzpzTO92/O3jVWzPKz8zy5TZG2kRHcmF6VX3gG1s7ji1Hyf0TeH+actYsGkX//hsDXPW5/Ln845gUKea/2D5xdjeFJWU8uwspxTw6S/XERcVyXXHhRbMlzWd2BWGAPAhnFk6FHijhscyxpgGKVAzGFXt1hSawazbkc+qrPw67/xhakdcdCQ3ndALcMbEq0siwsPnHsGBklIue3buoXEic/YeYNriTC4Y3rlGpWENVWSE8PglR9GxTQuumTKPSV//xKUju/LzSqZHDEWP5HjOPLITL3+/iUVbdvPeoq1cNrJrwMGkAykLAOulCtjH71T1AeAXOKWBxhjTpHmbwYjIkyIyBUBETgtz0qrtg8XbEKnf6l9TMxPHdOfjXx3HyJ7t6vxcPZLj+c9lw1BVfjn1R07+51fc/c5SDhaXctXo7nV+/nBp0zKa/14xnIPFpRzZpQ33nzWwVo9/y4m92XewhKsm/0B0ZAQ3HB96MF9vJYAicoaIzBWR1SLyBtAbQFX3ATeGfHZjjGncijg8TeVJ4UxIdakqHy7JZGSPJFJbx4U7OSZIERFSrx1cTh7Qnhm3n8Cky4eTEBfNjBXbObZ3cp0NGN5QDOjYms/vPIHXbhhV7fmLA+nXIYHTBrZnz/4iJozoWq3vX4voSGKjIqpVAljpMDAiMhD4vape7i56GrgcWIEzE8g/ROQpVZ2qqsHPWm2MMU1DAc40mdFA7Y7DUU9WZeXzU/Y+rj6mbjoSmKYjIkIYN7gDpw9qz8LNu0hrIsPrVKWqMShr4jen9yO/sJhfjO1Vrf1FhKT4GHL31nIACMwERnteb1fV79znn4vI98BcYGrIZzbGmMbvfpzpMZ8CXg1zWqpl+pJtRAicMbhDuJNiGgkRYXi3pHAno0no2z6BqTeMqtExElvGsKsO2gCeBvzZ83qjiDwsImV9lIuA/JDPaowxTYA7JuDTqnqDqn4Y7vSESlWZvnQbY3ol18k4csaYupcUH1P7w8Co6lJVvcy7CDgf2CIis4B1wFci0ifkMxtjTCMnIk815k4gK7fls2HnPn5mvX+NabQS42PYVVAU8n7BTgUHgKpOABCROGAwzlRIQ4BnRaSnqja9wYCMMSawg0DZTOwn4UyN2WhMX5pJZIRw+iCr/jWmsUpqGV2tEsCQAsAy7tRH892HMcY0V422E4iqMn3JNsb0aud3Hl5jTOOQGB/Dnv1FFJeUEhXCdHsNf2I+Y4xpIMTRV0T6uXP/3g/8RCPsBLI8M4+NOQX8zMb+M6ZRa+f+gNu9P7Rq4GqVABpjTHPjDov1FlAKCBAJnK+qT4c1YdU0fek2q/41pglI9AwGHcpMIlYCaIwxPkRkoIi87LP4eZxxUQer6iDgLuCV+k9dzZVV/x7TOzmkieeNMQ1PUkvnO5wTYjtACwCNMaaimcC9ZS9E5BdAZ6C1iFwpIlcCbYF2InJzdU4gIuPcWZXWicjdftY/JiKL3McaEdntWVfiWTct1HMv25rH5twCzrTqX2MavcRqTgdnVcDGGFNR2RioZcNgDQeigaN9totx14VERCJx2g2eCmQA80RkmqquKNtGVW/3bH8bcJTnEPtVdWio5y0zfek2oiKE0wa1r+4hjDENRFknrlCng7MA0BhjfKjqUg4Hf6jqdSLyE/BvVV0HICI9gTNV9dpqnGIEsE5V17vHeg04B2eaTX8m4HQ4qTFn8OdMjumdTNuWVv1rTGPXtmU0EHoJoFUBG2NMcO4AvhaRJ0XkCeAb4PYq9gmkM7DF8zrDXVaBiHQDegBfeBbHich8EZkjIueGcuKlW/ewJXe/Df5sTBMRGxVJq9gocvdZL2BjjKl1qvq+iCzFqR4GT2lgNYi/UwTY9hLgLVUt8SzrqqqZbinkFyKyVFV/KncCkRtw5imma9fDQxQuz8wDYEyvdtVMujGmoUmMjw55PmALAI0xJkhule2kWjhUBuCdOakLkBlg20uAW3zSkVmWHhH5Cqd94E8+2zwDPAOQnp5+KLjML3RKCaz615imI6ll6PMBWxWwMcbUv3lAHxHpISIxOEFehd68ItIPSAS+9yxLFJFY93kycAyB2w5WkF9YTIRAfExkDd+CMaahSIq3ANAYYxo8VS0GbgU+BVYCb6jqchF5UETO9mw6AXhNVb3VwwOA+SKyGPgSeMTbe7gqefuLaBUbhTORiTGmKUisRgBoVcDGGBMGqvoR8JHPsvt8Xj/gZ7/ZwBHVPW9+YTEJcdHV3d0Y0wAltYwJuQ1gWEoAqxoA1bPdz0VERSS9PtNnjDFNVV5hMQlx9tvfmKYkMT6GgoMlFBaVVL2xq94DQM8AqGcAA4EJ7hybvtslAL8E5tZvCo0xpunKLyyitZUAGtOklA0GHUopYDhKAA8NgKqqB4GyAVB9PQT8DSisz8QZY0xTlm8lgMY0OYlur/5Q2gGGIwCscgBUETkKSFPVD+szYcYY09TlHyiyANCYJuZQCWAIg0GHIwCsdABUEYkAHgPurPJAIje4o+HPz87OrsUkGmNM02SdQIxpepLine90KPMBhyMArGoA1ARgMPCViGwERgHT/HUEUdVnVDVdVdNTUlLqMMnGGNP4qapVARvTBB2qAt57IOh9whEAVjoAqqruUdVkVe2uqt2BOcDZqjo/DGk1xpgmY39RCSWlaiWAxjQxbVvGIAK5BQ24CjiEAVCNMcbUovzCYgArATSmiYmMENq2iGZXCJ1AwpILBDMAqmf52PpIkzHGNHVl8wBbAGhM05MYH9Pg2wAaY4wJgzy3BNDGATSm6UlqGRNSCaAFgMYY00xYFbAxTVeo8wFbAGiMMc3E4SpgKwE0pqkJdT5gCwCNMaaZsBJAY5ouKwE0xhjjl3UCMabpSoqPpqhEq97QZQGgMcY0E/mFxYhAfIwFgMY0NWWDQQfLAkBjjGkm8guLaRUbRUSEvxk5jTGNWdl8wMGyANAYY5qJvMIiGwLGmCbKAkBjjDF+2TzAxjRdFgAaY4zxK7+wyAJAY5qoRAsAjTGm4RORcSKyWkTWicjdftZPFJFsEVnkPq7zrLtKRNa6j6uCPWd+YbFVARvTRCXERhEVQvte+ylojDH1TEQigaeAU4EMYJ6ITFPVFT6bvq6qt3Ob9LAAABEcSURBVPrsmwTcD6QDCixw991V1XnzC4vpk2rZvjFNkYiEVApoJYDGGFP/RgDrVHW9qh4EXgPOCXLf04EZqprrBn0zgHHB7OhUAVsJoDFNVVIIQ8FYAGiMMfWvM7DF8zrDXebrAhFZIiJviUhaKPuKyA0iMl9E5mdnZ6Oq1gnEmCYuMT74H3gWABpjTP3z11DHdwj/D4Duqnok8DnwQgj7oqrPqGq6qqanpKRQWFRKcalaCaAxTVgoPYEtADTGmPqXAaR5XncBMr0bqGqOqh5wX/4PGB7svv7k2TRwxjR5l43sFvS2FgAaY0z9mwf0EZEeIhIDXAJM824gIh09L88GVrrPPwVOE5FEEUkETnOXVcrmATam6Tumd3LQ21pOYIwx9UxVi0XkVpzALRKYrKrLReRBYL6qTgN+KSJnA8VALjDR3TdXRB7CCSIBHlTV3KrOmVdYDGDDwBhjAAsAjTEmLFT1I+Ajn2X3eZ7fA9wTYN/JwORQzpfvBoBWAmiMAasCNsaYZuFwFbCVABpjLAA0xphmwUoAjTFeFgAaY0wzYJ1AjDFeFgAaY0wzkF9YjAjEx1gAaIyxANAYY5qF/MJiWsVGERHCZPHGmKbLAkBjjGkG8gqLbAgYY8whFgAaY0wzYPMAG2O8LAA0xphmIL+wyAJAY8whFgAaY0wz4JQAWhWwMcYRlgBQRMaJyGoRWScid/tZf4eIrBCRJSIyU0SCn93YGGNMBVYFbIzxqvcAUEQigaeAM+D/27v/YDvK+o7j708SkghpSiARMSAkBZEgGpgMoHQKaORH/kjiSDF0GKPSyYDSqe04UygdpmWmU9AZGR11IFqqVUeEWNprSaTRJO0fJZEMBEKI+YlT0kQDKlgKDST32z/2Obg5Oefec++eu7v3nM9r5k727D57nu+zu2fzPc/u2Yd5wPWS5jUVexJYEBHvAVYBnys3SjOz3uJLwGaWV0UP4EXA7ojYGxGvAw8AS/IFImJ9RLyaXm4ETis5RjOznuJLwGaWV0UCOBt4Pvd6X5rXzo3AmlYLJK2QtFnS5hdeeKGLIZqZ9Y7BgMOD4R5AM3tTFQlgq6eQRsuC0g3AAuDzrZZHxMqIWBARC2bNmtXFEM3MesfgYHaKdQ+gmTVU8XVwH3B67vVpwP7mQpIWArcDl0XEoZJiMzPrOUciSwCnuwfQzJIqegAfB86WNEfSZGAZMJAvIOkC4D5gcUQcrCBGM7OeceTNHkAngGaWKT0BjIjDwC3Ao8B24MGI2CbpTkmLU7HPA9OAhyRtkTTQ5u3MzGwYvgRsZs0q+ToYEauB1U3z7shNLyw9KDOzHtW4BOweQDNr8EggZmYVKPJAfElH0tWRjq6QNC4BT3cPoJkl/jpoZlay3APxP0T2w7jHJQ1ExLO5Yo0H4r8q6WayB+J/NC17LSLmd1rfoHsAzayJewDNzMpX6gPxjwwGEpww2QmgmWWcAJqZla/oA/Gnpofgb5S0tNUK+Qflv/K/rzFtyiQmTGj1GFYz60f+OmhmVr7RPBD/stzsd0TEfklzgXWStkbEnqPeLGIlsBLglN+bF77/z8zy3ANoZla+kT4Qf3H+gfgRsT/9uxfYAFwwVGWDg77/z8yO5gTQzKx8o34gvqQZkqak6ZnApUD+xyPHOOJxgM2sic8IZmYli4jDkhoPxJ8I3N94ID6wOSIGOPqB+AD/FRGLgXOB+yQNkn2Jv6vp18PHOBLhh0Cb2VGcAJqZVWC0D8SPiP8Ezh9JXYPuATSzJr4EbGbW47IeQCeAZvZbTgDNzHpc1gPoS8Bm9ltOAM3MelzgXwGb2dGcAJqZ9QH3AJpZnhNAM7M+MN09gGaW4wTQzKwP+BKwmeU5ATQz6wO+BGxmeU4Azcx63FuOm8jJJ0yuOgwzqxEngGZmPe6st05j7qxpVYdhZjXiBNDMzMyszzgBNDMzM+szTgDNzMzM+owTQDMzM7M+4wTQzMzMrM84ATQzMzPrM04AzczMzPqME0AzMzOzPuME0MzMzKzPVJIASrpa0g5JuyXd2mL5FEnfS8s3STqz/CjNzMZOkfOgpNvS/B2SriozbjPrDaUngJImAl8BrgHmAddLmtdU7Ebg1xFxFnAPcHe5UZqZjZ0i58FUbhlwHnA18NX0fmZmHauiB/AiYHdE7I2I14EHgCVNZZYA30zTq4APSlKJMZqZjaUi58ElwAMRcSgingN2p/czM+vYpArqnA08n3u9D7i4XZmIOCzpZeBk4MV8IUkrgBXp5SFJz4xJxJ2ZSVN8fVJ31fX3fNvVvv+759s+hHMqqrdbipwHZwMbm9ad3VxBzc6Pnaj6eB6O4yvG8RUzkvjO6KRQFQlgq568GEUZImIlsBJA0uaIWFA8vNGpsn633W3vt/olba6i3i4qch4cd+fHTtQ9RsdXjOMrZiziq+IS8D7g9Nzr04D97cpImgT8LvCrUqIzMxt7Rc6DnaxrZjakKhLAx4GzJc2RNJnsZuaBpjIDwPI0fS2wLiKO+YZrZjZOFTkPDgDL0q+E5wBnAz8pKW4z6xGlXwJO97LcAjwKTATuj4htku4ENkfEAPD3wLck7Sb7xrusg7deOWZBd6bK+t32/qzfbR+nipwHU7kHgWeBw8CnI+LIMFWOh+1V9xgdXzGOr5iuxyd3rJmZmZn1F48EYmZmZtZnnACamZmZ9ZlxlQBK+kNJ2yQNSmr7c+h2QyylG643SdqVhliaPML6T5K0Nq2/VtKMFmWukLQl9/d/kpamZd+Q9Fxu2fxu1p3KHcm9/0Bufhltny/psbSPnpb00dyyEbe9yqGyOqj7zyU9m9r5Y0ln5Ja13Addrv/jkl7I1fPHuWXL037aJWl587pdqPueXL07Jb2UW1ao7ZLul3RQbZ5Zp8yXUmxPS7owt6xQu3tBkc9MTeJre1yXFN+oj7+axHe5pJdz2++OEmM7XdJ6SdvT/wF/2qJMZduvw/gq236p/qmSfiLpqRTj37Qo073PcESMmz/gXLIHwG4AFrQpMxHYA8wFJgNPAfPSsgeBZWn6XuDmEdb/OeDWNH0rcPcw5U8iu3n7+PT6G8C1o2x7R3UDr7SZP+ZtB94JnJ2m3w4cAE4cTduH2o+5Mp8C7k3Ty4Dvpel5qfwUYE56n4ldrvuK3H69uVH3UPugy/V/HPhym2Nub/p3Rpqe0c26m8r/CdkPGLrV9j8ALgSeabN8EbCG7Fl4lwCbutHuXvgr8pmpUXwtj+sSYxzV8Vej+C4H/rWibXcqcGGa/h1gZ4v9W9n26zC+yrZfql/AtDR9HLAJuKSpTNc+w+OqBzAitkfEjmGKtRxiSZKAD5ANqQTZEEtLRxhCfmimTta/FlgTEa+OsJ5u1P2mstoeETsjYlea3g8cBGaNsJ6GKofKGrbuiFif268byZ7F1i2dtL2dq4C1EfGriPg1sJZsvNixqvt64LsjeP8hRcR/MPQzP5cA/xiZjcCJkk6leLt7Qd2H2SxyXJeiwPFXig7iq0xEHIiIJ9L0/wDbOXaEmsq2X4fxVSptl1fSy+PSX/Mvdbv2GR5XCWCHWg2xNJtsCKWXIuJw0/yROCUiDkB2MAFvHab8Mo79z/FvU9f3PZKmjEHdUyVtlrRR6dIzFbRd0kVk3/L35GaPpO3t9mPLMqlt+aGyhlu3aN15N5J9q21otQ9GotP6P5K25ypJjQcDl9Z2ZZe95wDrcrOLtn208RVtdy8o8pkpQ5Hjui7Gw3H2vnQJcY2k86oIIF2WvICsByuvFttviPig4u0naaKkLWQdKGsjou02LPoZrmIouCFJ+hHwthaLbo+If+nkLVrM63j4pKHq76Du/PucCpxP9pyvhtuAn5MlRiuBvwDu7HLd74iI/ZLmAuskbQV+06LcWLf9W8DyiBhMs4dse6u36SDmQvu6YN1ZQekGYAFwWW72MfsgIva0Wr9A/T8AvhsRhyTdRPaN8AMjib1A3Q3LgFVx9DPoirZ9tPEVbXcv6Nowm2OkyHFdF3U/zp4AzoiIVyQtAv6Z7EHhpZE0Dfg+8JmIaP6/p/LtN0x8lW+/dD6dL+lE4GFJ746I/D2fXduGtUsAI2JhwbdoN0zSi2TdzZNS1txy+KSh6pf0C0mnRsSBlOQcHCKO64CHI+KN3HsfSJOHJP0D8Nlu150uvRIReyVtIPuW831Karuk6cAjwF+lLv6O2t7CSIbK2qfuDpXV0fqSFpIlx5dFxKHG/Db7YCRJ0LD1R8Qvcy+/BtydW/fypnU3dLPunGXAp5viKtr20cZXtN29oMhnpgxFjuu6qPUwfPmEJiJWS/qqpJkR8WIZ9Us6juz/m+9ExD+1KFLp9hsuvqq3X1MsL6Vz6NVAPgHs2me4Fy8BtxxiKSICWE92Xx5kQyx10qOYlx+aabj1j7k3qnGvQ7pev5Sjd2rhuiXNaFxalTQTuBR4tqy2p+39MNk9Hg81LRtp26scKmvYuiVdANwHLI6Ig7n5LffBCOrutP78fTOLye5ngazH+coUxwzgSo7uhS5cd6r/HLIfWzyWm9eNtg9nAPiYMpcAL6cvF0Xb3QvqPsxmkeO6Ltodf7Ug6W3pHNu4DWcC8Muh1+pa3SIbvWZ7RHyhTbHKtl8n8VW5/VKds1LPH5LeAiwEftpUrHuf4ajo1y6j+QM+TJb9HgJ+ATya5r8dWJ0rt4jsFz57yC4dN+bPJUsEdgMPAVNGWP/JwI+BXenfk9L8BcDXc+XOBP4bmNC0/jpgK1ny823Sr326VTfw/vT+T6V/byyz7cANwBvAltzf/NG2vdV+JLtsvDhNT01t2Z3aNje37u1pvR3ANaM41oar+0fpGGy0c2C4fdDl+v8O2JbqWQ+8K7fuJ9M22Q18ott1p9d/DdzVtF7htpN9aTqQjqN9ZPdX3gTclJYL+EqKbSu5pwEUbXcv/BX5zNQkvrbHdUnxjfr4q0l8t+S230bg/SXG9vtklyKfzp0XF9Vl+3UYX2XbL9X/HuDJFOMzwB1p/ph8hj0UnJmZmVmf6cVLwGZmZmY2BCeAZmZmZn3GCaCZmZlZn3ECaGZmZtZnnACamZmZ9RkngGZmZmZ9xgmgmZmZWZ9xAmi1lwbH/qKkbZK2pnFmG8vuk3Rpmv6IpE3KBvLeLOmq6qI2M6uGpDMlrZa0Q9JOSbdVHZPVjxNAGw9uA/ZGxHnAl4BP5ZZdDGyU9Edk4wsviYj3kg3F901Jp5UerZlZRSRNIBvv9t6IOAc4H1ggaUW1kVndOAG0WpN0AvDhiPhimvUccFZadi7ZsFJTgbuA6yLi5wARsQvYAHyw7JjNzCp0FfCziBgAiIhDZEOcfbbSqKx2JlUdgNkwFgKnS9qSXp9ENg4vwDXAD8kGlX8iIp5vWvcQcHwpUZqZ1cO5ZGPZvikiDkiaLmlyRLxeUVxWM+4BtLqbTzYg9vyImA/8G9kg3pB90/0h8G6aTnjJe4GflhKlmVk9HAGm5WdIEtmX4cOVRGS15ATQ6m4G8CqApEnAlcAPJB0PnBgR+4HfAJPzK0l6H3AC8O/lhmtmVqkNwKKU9DV8iOwqyWA1IVkdOQG0utsJXJKm/wx4JCKeA64A1qf5jwDXSZoFIOmdwNeBT/qEZ2b9JCKeAp4E7gSQdArwBeAvq4zL6kcRUXUMZm1JmgGsAWYCjwErIuI1SV8GVkXEhlTuE8BngCnAROBjEfFYNVGbmVVD0q3AcuBdZD+Cu4WsB3APsDQiflZddFYnTgBtXJL0BHBxRLzRNP8UYB1wfUQ8XUlwZmZmNecE0MzMzKzP+B5AMzMzsz7jBNDMzMyszzgBNDMzM+szTgDNzMzM+owTQDMzM7M+4wTQzMzMrM84ATQzMzPrM/8PMoTC/clEVXkAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 648x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAooAAAGDCAYAAACskzHZAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3XucXVV98P/PlyGhU1SGm5cMRLCk+YlijaaAxVZ+eElorcnDg0/RKqFikVaftj9tNKn+1AelhObXYq1opUrFK1BKQ34VjdRIW61cYsNDBE2JgCQTFBAGEUYCyXr+2GvIyck+M2cm57LPmc/79ZrXnLP22muvvc/tu/e67EgpIUmSJNXbr9sVkCRJUjUZKEqSJKmUgaIkSZJKGShKkiSplIGiJEmSShkoSpIkqZSBojTDRUSKiGOmue7vRsTXWl2nXPYfRMSPI+JnEXFoO7bRYLt/FhGfakO5Z0XEN1tdblW067i1UkT8ekRs7sB2jsqfq/07sK0PRsTn270dzVwGippxIuKNEbEhByD3RsRXIuLl3a5X1ZX9+KWUvpBSek0btjUL+CvgNSmlp6WUftLqbeTtnBwR22rTUkp/nlJ6azu2N129EGR26riVvWbNSin9e0ppfqvrJPUzA0XNKBHxTuAjwJ8DzwLmAh8HljTI3/YrAq1SVtdeqn+dZwG/ANzW7YpIM0Grvj8iYqA1NVJVGChqxoiIg4DzgLenlK5OKT2aUnoipfT/p5SW5zwfjIirIuLzEfFT4KyIOCAiPhIR2/PfRyLigJz/sIj454gYjYgHI+LfI2K/vOw9ETESEY9ExOaIeOUEdXt5RPxHLmdrRJw1XueI+GxE3B8RP4yI99WUf1ZEfCsiLoqIB4EPlqXlvG+JiO9FxEMRsS4intugHr8VERsj4qe5Hh+sWfxv+f9ovhr7svorXRHxaxFxc0Q8nP//Ws2y6yPiQ7l+j0TE1yLisJI6/DKwuWZb68uuZuby3lpzLL4ZEf9f3se7IuLUmryHRMTf59fvoYhYExEHAl8B5uT9+VlEzKlvyouI10XEbfm1uT4inl+z7O6I+NOIuDXv8xUR8QuNXudilfibnPf7te+J/Fp/Ooqr3CMR8eGIGMjb+1vgZbmOoxFxdP4//l74VETcV1PW5yPiTyYqtyZvw/dGPubnRsQdefnFERENduyp41bzei2LiHsi4oGIeG9d3qvy8XokIv4zIn6lbrvH1Dz/TK536WtWUpffjIjbc9kjEfGnOX2Pq5ER8ZL8fn8kIv4h1+fDtXkj4l0RcV8+fr9Xs+5En5UJRcSKiPhB3u7tEfHfapZN9l4+OiL+Na97HbDXZ6huW6+NiFvy++U/IuJFNcvujuJ76lbg0YjYv0Ha8/N7fzSKz8Lr6l6bT0TEtRHxKPB/N3sc1CNSSv75NyP+gMXAk8D+E+T5IPAEsJTiRGqQIri8AXgmcDjwH8CHcv4LKH7EZ+W/XwcCmA9sBebkfEcBv9Rgm3OBR4A35DIOBV6cl30WuAZ4ei7jv4Cz87Kz8v78T2D/XNeytKXAFuD5Oe19wH/UbD8Bx+THJwPH5X1/EfBjYGnNPqTa45e39838+BDgIeDNeTtvyM8PzcuvB34A/HKu1/XAqgbHZI9tNdj29cBba+rxBPD7wADwB8B2IPLyLwNXAAfnY/yKmv3dVvIe+Hx+/MvAo8Cr83rvzsdydl5+N3ATMCfv//eAcxvs0/hr8//ksn4HeBg4JC9fA3wSOJDivXYT8Lb641xT3j3AS/PjzcCdwPNrli1ootxm3hv/DAxRvE/vBxZP8NkZP27jr9ff5df6V4DHa+r3wfx6nZ6PxZ8CdwGz6t+T+flngA83es1K6nIv8Ov58cHAS+rXBWYDPwT+ONfhNGBH3XaepPj8zwJ+E3gMOHg6n5W6+r2e4j2zX34fPAo8p8n38rcpumUcAPwGxXfH5xts5yXAfcAJuaxlFO/ZA2rev7cARwKDZWl537cAf5aP2Sl5m/NrXpuHgZPy/vxCt7/r/WvtX9cr4J9/nfoDfhf40SR5Pgj8W13aD4DfrHm+CLg7Pz6PIpA7pm6dY/IX9KvIP34TbHMl8E8l6QMUP67H1qS9Dbg+Pz4LuKdunbK0r5CDy/x8v/yD99z8fI8f5bp1PwJclB/v9ePHnoHim4Gb6tb/NnBWfnw98L6aZX8IfLXBdvfYVoNtX8+egeKWmmW/mPM/G3gOsIv8A1+3nZOZOFD8f4Er647dCHByfn438Kaa5X8B/G2DfTqLmh/8nHZTPm7Pyq/1YM2yNwDfqD/ONcs/B7wz7+PmvO1zgaOB0VzXycpt5r3x8prlVwIrJvjs1AeKR9Tt6xk1eW+o225tcLevgeI9FJ+VZzR6vSmCrJG61+ObddsZY8/33H3AidP5rExS31uAJU28l+dSBK8H1iz/Io0DxU+QT2pr0jaz+0TpbuAtdcv3SKM4+f0RsF9N2peAD9a8Np9tZj/9680/m541k/wEOCwm73ezte75HIorD+N+mNMAVlOcbX8tIu6MiBUAKaUtwJ9Q/CDeFxGXjzeR1TSZ/Swi5lKcuf+gpB6HsfuqR+22hyeoa1nac4G/zs1Go8CDFFc9h+tXjIgTIuIbUTR1P0wReEzYtFWj/jiV1fdHNY8fA57WZNnNeKrslNJj+eHTKI7vgymlh6ZR5h77lFLaRXF8p7tPIyn/umbj76XnUly5ubfmdfokxRXARv6VIpj5DYpuAdcDr8h//57rOlm5zbw39uU1m2jdp96nua7b2P252lf/neIK4A9zM+3LSvLMYe/Xo/6z85OU0pM1z5/ah335rETEmTXNwaPAC+vWbfRengM8lFJ6tCZv/Weu1nOBd41vJ2/rSPY8zpN9h8wBtubXqHabk30PqU8YKGom+Tbwc4rmtomkuufbKb5wx83NaaSUHkkpvSul9Dzgt4F3Ru53llL6Ykrp5XndBFyY059W83cPxZfsL5XU4wGKJqj6bY9MUNeytK0UTY1DNX+DKaX/KFn3i8Ba4MiU0kEUzerjfdLKtlWr/jiV1Xe6xn8Yf7Em7dlNrrsVOCQihkqWTWmfcv+8I5n+Pg3X9fEbfy9tpbjyd1jNa/SMlNILJqjnv1Jc7Tk5P/4mRfPfK/Jzmih3Ku+NVjty/EEUfS2PIH+uKAKyRq/1ZK8ZKaWbU0pLKALiNRRXQuvdy96vx5El+RqZ6LPSUBR9QP8OeAdFt4wh4LvNrJvrfHDuqzlu7gT5twLn172+v5hS+lJNnsm+Q7YDR+bXqHabk30PqU8YKGrGSCk9DLwfuDgilkbEL0bErIg4NSL+YoJVvwS8LyIOj2LwxfuB8U77r42IY/KPzU+BncDOiJgfEadEMejl5xRNWDsblP8F4FUR8T9yx/FDI+LFKaWdFD9w50fE0/MPzDvHtz0FfwusjIgX5DofFBGvb5D36RRX334eEccDb6xZdj9FE+7zGqx7LfDLUUw/tH9E/A5wLEUft32SUrqf4ofpTVEM8HgL5cF12br3UjSxfjwiDs6v+W/kxT8GDo1ioFOZK4HfiohXRjFlz7soAq/pBlLPBP4o1+H1FH0Dr811/BrwlxHxjIjYLyJ+KSJeUVPPIyJids1+3UHxvnoTRXeJn+Z8/50cKDZR7lTeG6320og4LV/h/xOK43pDXnYL8Mb8Wi+mCH7HTfiaRcTsKOb3PCil9AS7P5f1vp3T35Hfr0uA46dQ/4k+KxM5kCKwuj/X9/corihOKqX0Q2AD8L/yfr6c4gS1kb8Dzs1XPyMiDoxiEM7Tm6wrwI0UJ2rvzu/bk/M2L59CGephBoqaUVJKf0URbL2P4ot6K8WZ/ZoJVvswxZfzrcAm4D9zGsA84F+An1H88Hw8pXQ9RUfzVRRXBX9EESD8WYM63UPRTPYuiqa/Wyg6/0MxKOVRioEK36S4inHpFPf5nyiuZl4exUju7wKnNsj+h8B5EfEIRUD81JWY3AR2PvCt3Ix1Yt12fgK8Nu/HTygGfrw2pfTAVOo7gd8HlueyX8DUgrU3U1yd/T5FP7M/yXX+PsWJwJ15n/Zo+kwpbaYIxP6G4rX8beC3U0o7prkPN1K8Zx6gOJanp91zRJ5J0dXgdopBQFdR9K8EWE8xVdCPIqL2eP4rRfPoPTXPA9hYk6dhuVN8b7TaNRQDOcYHQJ2WAzsoBpj8NkVfy9+l5vM52WuWvRm4O+/TuRSv4R7ya3gacHbezpsoTmoeb7L+DT8rE0kp3Q78JcX3xY8pBsR8q8ltQhGQnkDxXfEBigFvjba1geJz8zGK47yFog9k0/Jxeh3F++IBiunEzsyvg2aA8VFUkiR1RBRTyRyTUtorgOumiLiRYjDS33e7LlJVeEVRkjQjRcQrIuLZuel5GcU0N1/tdr2kKunVuzZIkrSv5lM0GT+NYuaB03O/TkmZTc+SJEkqZdOzJEmSShkoSpIkqZR9FFvgsMMOS0cddVS3qyFJkjSp73znOw+klA5vJm+lAsU8sepfU9zj9lMppVV1yw+gmDPqpRRzqf1OSunuvGwlxXxYO4E/Simtm6jMiDiaYsLQQyjmxXtz7dxoEXE68A/Ar+a5qBo66qij2LBhwiySJEmVEBET3fpxD5Vpeo6IAeBiikk9jwXeEBHH1mU7m+I+l8cAF5FviZbznUExCe9iijswDExS5oUUN3CfRzER6dk1dXk68EcUk+NKkiTNSJUJFClunbQlpXRnvrJ3ObCkLs8S4LL8+CrglfnWaUuAy1NKj6eU7qKYff74RmXmdU7JZZDLrL3/74eAv6C49ZokSdKMVKVAcZjidmrjtuW00jwppSeBh4FDJ1i3UfqhwGguY49tRcQCipu8T3h/2og4JyI2RMSG+++/v9l9lCRJ6hlVChSjJK1+ksdGeVqSHhH7UTRpv2uCehaZU7okpbQwpbTw8MOb6g8qSZLUU6oUKG4Djqx5fgSwvVGeiNgfOIjixuiN1m2U/gAwlMuoTX868ELg+oi4GzgRWBsRC/dx3yRJknpOlQLFm4F5EXF0RMymGJyyti7PWmBZfnw6sD4Vt5ZZC5wREQfk0czzgJsalZnX+UYug1zmNSmlh1NKh6WUjkopHQXcALxuslHPkiRJ/agy0+OklJ6MiHcA6yimsrk0pXRbRJwHbEgprQU+DXwuIrZQXEk8I697W0RcCdwOPAm8PaW0E6CszLzJ9wCXR8SHgY25bEmSJGXe67kFFi5cmJxHUZIk9YKI+E5KqaludVVqepYkSVKFVKbpWZIkaaZas3GE1es2s310jDlDgyxfNJ+lC+pnCew8A0VJkqQuWrNxhJVXb2LsiZ0AjIyOsfLqTQBdDxZtepYkSeqi1es2PxUkjht7Yier123uUo12M1CUJEnqou2jY1NK7yQDRUmSpC6aMzQ4pfROMlCUJEnqouWL5jM4a2CPtMFZAyxfNL9LNdrNwSySJEldND5g5d1X3cqOnbsYdtSzJEmSxi1dMMyXbroHgCve9rIu12Y3m54lSZJUykBRkiRJpQwUJUmSVMpAUZIkSaUMFCVJklTKQFGSJEmlDBQlSZJUykBRkiRJpQwUJUmSVMpAUZIkSaUMFCVJklTKQFGSJEmlDBQlSZJUykBRkiRJpQwUJUmSVMpAUZIkSaUMFCVJklTKQFGSJEmlDBQlSZJUykBRkiRJpQwUJUmSVMpAUZIkSaUMFCVJklTKQFGSJEmlDBQlSZJUykBRkiRJpQwUJUmSVKpSgWJELI6IzRGxJSJWlCw/ICKuyMtvjIijapatzOmbI2LRZGVGxNG5jDtymbNz+rkRsSkibomIb0bEse3da0mSpGqqTKAYEQPAxcCpwLHAG0qCtLOBh1JKxwAXARfmdY8FzgBeACwGPh4RA5OUeSFwUUppHvBQLhvgiyml41JKLwb+AvirtuywJElSxVUmUASOB7aklO5MKe0ALgeW1OVZAlyWH18FvDIiIqdfnlJ6PKV0F7All1daZl7nlFwGucylACmln9Zs70AgtXg/JUmSesL+3a5AjWFga83zbcAJjfKklJ6MiIeBQ3P6DXXrDufHZWUeCoymlJ4syU9EvB14JzCbIqCUJEmacaoUKEZJWv3VvEZ5GqWXXTGdKH/xIKWLgYsj4o3A+4Ble1U24hzgHIC5c+eWFNd/1mwcYfW6zWwfHWPO0CDLF81n6YLhyVeUJEk9qUpNz9uAI2ueHwFsb5QnIvYHDgIenGDdRukPAEO5jEbbgqKpemlZZVNKl6SUFqaUFh5++OGT7lyvW7NxhJVXb2JkdIwEjIyOsfLqTazZONLtqkmSpDapUqB4MzAvj0aeTTE4ZW1dnrXsvrp3OrA+pZRy+hl5VPTRwDzgpkZl5nW+kcsgl3kNQETMq9nebwF3tHg/e9LqdZsZe2LnHmljT+xk9brNXaqRJElqt8o0Pec+h+8A1gEDwKUppdsi4jxgQ0ppLfBp4HMRsYXiSuIZed3bIuJK4HbgSeDtKaWdAGVl5k2+B7g8Ij4MbMxlA7wjIl4FPEExGnqvZueZaPvo2JTSJUlS76tMoAiQUroWuLYu7f01j38OvL7BuucD5zdTZk6/k2JUdH36H0+54jPAnKFBRkqCwjlDg12ojSRJ6oQqNT2rwpYvms/grIE90gZnDbB80fwu1UiSJLVbpa4oqrrGRze/+6pb2bFzF8OOepYkqe8ZKKppSxcM86Wb7gHgire9rMu1kSRJ7WbTsyRJkkoZKEqSJKmUgaIkSZJKGShKkiSplIGiJEmSShkoSpIkqZSBoiRJkkoZKEqSJKmUgaIkSZJKGShKkiSplIGiJEmSShkoSpIkqZSBoiRJkkoZKEqSJKmUgaIkSZJKGShKkiSplIGiJEmSSu3f7QpIkiS1w5qNI6xet5nto2PMGRpk+aL5LF0w3O1q9RQDRUmS1HfWbBxh5dWbGHtiJwAjo2OsvHoTgMHiFNj0LEmS+s7qdZufChLHjT2xk9XrNnepRr3JQFGSJPWd7aNjU0pXOQNFSZLUd+YMDU4pXeUMFCVJUt9Zvmg+g7MG9kgbnDXA8kXzu1Sj3uRgFkmS1HfGB6y8+6pb2bFzF8OOep4WA0VJktSXli4Y5ks33QPAFW97WZdr05tsepYkSVIpA0VJkiSVMlCUJElSKQNFSZIklTJQlCRJUilHPUvqqDUbR1i9bjPbR8eY43QVklRpBoqSOmbNxhFWXr3pqfuvjoyOsfLqTQAGi1If8YSwf9j0LKljVq/b/FSQOG7siZ2sXre5SzWS1GrjJ4Qjo2Mkdp8Qrtk40u2qaRoMFCV1zPbRsSmlS+o9nhD2l0oFihGxOCI2R8SWiFhRsvyAiLgiL78xIo6qWbYyp2+OiEWTlRkRR+cy7shlzs7p74yI2yPi1oj4ekQ8t717Lc0cc4YGp5Quqfd4QthfKhMoRsQAcDFwKnAs8IaIOLYu29nAQymlY4CLgAvzuscCZwAvABYDH4+IgUnKvBC4KKU0D3golw2wEViYUnoRcBXwF+3YX2kmWr5oPoOzBvZIG5w1wPJF87tUI0mt5glhf6lMoAgcD2xJKd2ZUtoBXA4sqcuzBLgsP74KeGVERE6/PKX0eErpLmBLLq+0zLzOKbkMcplLAVJK30gpPZbTbwCOaMO+SjPS0gXDXHDaccweKL56hocGueC04+zkLvURTwj7S5VGPQ8DW2uebwNOaJQnpfRkRDwMHJrTb6hbd/yXp6zMQ4HRlNKTJflrnQ18payyEXEOcA7A3LlzJ9ovSTWWLhjmSzfdA8AVb3tZl2sjqdXGT/zefdWt7Ni5i+EZMOq5n0d5VylQjJK01GSeRullV0wnyr97QxFvAhYCryjJS0rpEuASgIULF9bXU5KkGWsmnRD2+7RfVWp63gYcWfP8CGB7ozwRsT9wEPDgBOs2Sn8AGMpl7LWtiHgV8F7gdSmlx/dpryRJUt/q91HeVQoUbwbm5dHIsykGp6yty7MWWJYfnw6sTymlnH5GHhV9NDAPuKlRmXmdb+QyyGVeAxARC4BPUgSJ97VpXyVJUh/o91HelQkUc3/BdwDrgO8BV6aUbouI8yLidTnbp4FDI2IL8E5gRV73NuBK4Hbgq8DbU0o7G5WZy3oP8M5c1qG5bIDVwNOAf4iIWyKiPliVJEkC+n+Ud5X6KJJSuha4ti7t/TWPfw68vsG65wPnN1NmTr+TYlR0ffqrplxxSZI0Iy1fNH+PPorQX6O8KxUoSpIk9ZJ+H+VtoChJkrQP+nmUd2X6KEqSJKlaDBQlSZJUyqZnSZL0lH6+y4imzkBRkiQB/X+XEU2dTc+SJAno/7uMaOoMFCVJEtD/dxnR1BkoSpIkoP/vMqKpM1CUJElAcZeRwVkDe6T1011GNHUOZpEkSUD/32VEU2egKEmSntLPdxnR1Nn0LEmSpFIGipIkSSploChJkqRSBoqSJEkqZaAoSZKkUgaKkiRJKmWgKEmSpFIGipIkSSploChJkqRS3plFarM1G0dYvW4z20fHmOPtsCRJPcRAUWqjNRtHWHn1Jsae2AnAyOgYK6/eBGCwKEmqPJuepTZavW7zU0HiuLEndrJ63eYu1UiSpOYZKEpttH10bErpkiRViYGi1EZzhganlC5JUpUYKEpttHzRfAZnDeyRNjhrgOWL5nepRpLUG9ZsHOGkVes5esWXOWnVetZsHOl2lWYkB7NIbTQ+YOXdV93Kjp27GHbUsyRNyoGA1WGgKLXZ0gXDfOmmewC44m0v63JtJKn6JhoIaKDYWTY9S5KkSnEgYHUYKEqSpEpxIGB1GCiqL9kJWpJ6lwMBq8M+iuo7doKWpN7mQMDqMFBU37ETtCT1PgcCVoNNz+o7doKWJKk1DBTVd+wELUlSa1QqUIyIxRGxOSK2RMSKkuUHRMQVefmNEXFUzbKVOX1zRCyarMyIODqXcUcuc3ZO/42I+M+IeDIiTm/vHqsd7AQtVZuDzaTeUZlAMSIGgIuBU4FjgTdExLF12c4GHkopHQNcBFyY1z0WOAN4AbAY+HhEDExS5oXARSmlecBDuWyAe4CzgC+2Yz/VfksXDHPBaccxe6B4ew8PDXLBacft1T/RHyup88YHm42MjpHYPdjMz59UTVUazHI8sCWldCdARFwOLAFur8mzBPhgfnwV8LGIiJx+eUrpceCuiNiSy6OszIj4HnAK8Mac57Jc7idSSnfnvLvasI/qkMk6QTsyWp20ZuMIq9dtZvvoGHNm+OhNB5tJvaUyVxSBYWBrzfNtOa00T0rpSeBh4NAJ1m2UfigwmstotC31sYl+rKRW6rcraPt6Jd7BZlJvqVKgGCVpqck8rUpvWkScExEbImLD/fffP5VVVQH+WKlT+umkpBVBr4PNpN5SpUBxG3BkzfMjgO2N8kTE/sBBwIMTrNso/QFgKJfRaFsTSildklJamFJaePjhh09lVVWAP1bqlH46KWlF0OtgM6m3VClQvBmYl0cjz6YYnLK2Ls9aYFl+fDqwPqWUcvoZeVT00cA84KZGZeZ1vpHLIJd5TRv3TRXjj5U6pZ9OSloR9DY72ExSNVQmUMz9Bd8BrAO+B1yZUrotIs6LiNflbJ8GDs2DVd4JrMjr3gZcSTHw5avA21NKOxuVmct6D/DOXNahuWwi4lcjYhvweuCTETGeX33EHyt1Sj+dlLQq6F26YJgFc4c44ehD+NaKU/zcSRVWpVHPpJSuBa6tS3t/zeOfUwRwZeueD5zfTJk5/U52j4yuTb+Zoilafc7bQ6kT+umetcsXzd9jtgDo3aBXUnMqFShKUj/ql5OSfgp6JTXHQFGS1LR+CXolNacyfRQlSZJULZNeUYyIQ5ooZ1dKabQF9ZEkdZF3kZFUq5mm5+35r2yS6nEDwNyW1EiS1BXe2lJSvWYCxe+llBZMlCEiNraoPpKkLvE+zJLqNRMoNtNb2R7N6iibx6TW66e7yEhqjUkHs+S5C4mIJfXLImK/2jxSJ7TifrOS9tZPd5GR1BpTGfX8tog4ASAiBiLiLcD321MtqbFW3G9W0t766S4yklpjKvMovgG4JiK+DPwBsAk4sy21kiZg89jemmmKt7lek3FCbUn1phIongS8F/gCcFZK6fq21EiaxJyhQUZKgsKZ2jzWzEhVR7OqWU6oLanWVJqeTwc+AhwGfDIiPhYR57anWlJjNo/tqZmmeJvrJUnT0fQVxZTSWwAiIoB5wHH5T+oom8f21ExTvM31kqTpmPK9nlNKCfiv/PePLa+R1ASbx3Zrpine5npJ0nRM2vQcEf/ZijyS2qOZpnib6yVJ09HMFcXnR8StEywP4KAW1UfSFDXTFG9zvSRpOpoJFP+vJvLsnDyLpHZppine5nr1G6d8ktpv0kAxpfTDTlREkqRmOeWT1BlTmR7nKRFxcUR8Jj9+TUtrJEnSJFo15dOajSOctGo9R6/4MietWu+tQKU60woUgR3AnfnxKS2qiyRJTWnFlE/eN16a3HQDxceAgyJiFjC3hfWRJGlSjaZ2msqUT/02Eb1XR9UO0w0UPwD8APg4xS39JEnqmFZM+dSqieirEKB5dVTtMt1A8W+Bz6SUfh94pIX1kSRpUksXDHPBaccxe6D4GRseGuSC046b0kCWVlyVrEqA1m9XR1Ud+3JF8dMR8TngV1tYH0mSmrJ0wTAL5g5xwtGH8K0Vp0x5tHMrrkpWJUDzNp1qlynfwi/7ELAZeB5wZeuqI0lSZ7RiIvqqBGidvE2n81fOLFMKFCPiV4HZwHtSSvdHxIHAXwNvbUflJElqp32diL4q91Ffvmj+HvNKQntu0+n8ldPTy8F1003PEbEG+ATwp8A3I+JSIAFva1PdJEmqtKrcR70VfTabUZWm9l5SlX6s0zWVK4ovBH4X2JhS2hERbwE+kVJa1p6qSTNDL59pSjNdle6j3onbdFalqb2XTBRc98J3/VQCxVXAe4BfiYjHgE3AKyLi1cAtKaX721FBqZ/ZjCP1vpl0H/WqNLX3kl4Prqcy6vm3gStTSr8EvJxiipwB4Azgq22om9T3bMaR1Euq0tTeS1oxDVM3TSVQfCvwpoj43xQTbf8V8NWU0tkppZe2pXZSn+v1M01JM0un+kL2k14Prptues5Ny6+NiDnAccDRaGfpAAAYcklEQVTDKaUb2lYzaQawGUdSr+m3pvZ29xOvUj/W6ZjyPIoppe3A9jbURZpxOjWlhSRpb53qJ97LwfV0J9yW1AK9fqYptYuzAagTen1EcicYKEpd1stnmlI7OBuAOsV+4pOb7r2eJUlqC2cDUKf0+ojkTjBQlCRVild51Cm9PiK5EyoVKEbE4ojYHBFbImJFyfIDIuKKvPzGiDiqZtnKnL45IhZNVmZEHJ3LuCOXOXuybXTLmo0jnLRqPUev+DInrVrfM7f9kaTp8CqPOsXpfiZXmUAxIgaAi4FTgWOBN0TEsXXZzgYeSikdA1wEXJjXPZZi4u8XAIuBj0fEwCRlXghclFKaBzyUy264jW7p9XtESmqNmXTC6FUeddLSBcMsmDvECUcfwrdWnGKQWCdSSt2uAwAR8TLggymlRfn5SoCU0gU1edblPN+OiP2BHwGHAytq847ny6vtVSbF7QjvB56dUnqydtuNtpEmOFALFy5MGzZsaMlxqHfSqvWMjI7xtluv4XkP7/5hOGD/ARbMHXrq+d0/eRSAow49sLScB372OHc98Cg7dyUO2H+AIw8Z5LCnHbBXvsnKuf3enwJw7HOe0bDOk5Ux2fJm80xWl2bqOlmeKtW1Fa9NVV6/VtSjU3Xt1HtgojIe+Nnj3PnAo+zatftraL/9gucdduAen+N+er8+8LPH+cH9j5LSxN9ZrfgMd6KMZsqp0merKt+vnXptOrGdZuvxo8OP5Pe+9DcN87RCRHwnpbSwmbxVGvU8DGyteb4NOKFRnhzgPQwcmtNvqFt3/JSgrMxDgdGU0pMl+Rtt44HaikTEOcA5AHPnzp3Kfk5Joz45jz+5Z0fvx3bsLM0He//IPP7kTu58oHjT13/xTlQOwC/OHphweTNlTLa82TyT1aWZuk6Wp0p1bcVrU5XXrxX1aCZPVcqAfXuvbX1wbI8gEWDXrsTWB8f2+Az30/v1sKcdwM8eL76iJ/pxbsVnuBNlNFNOlT5bVfl+7dRr04ntNFuPQw7c+4Som6oUKEZJWv1VvEZ5GqWXNa1PlL/ZepBSugS4BIoriiXrtMT4nTs++aIle6QPDw3yOytOeer5uz/5baB8epU35quS9YaHBvlWTRmTlQPw3CbqPFkZky1vNs9kdWmmrpPlqVJdW/HaVOX1a0U9mslTlTJg395rJ6/48t5fQhRfVnet+q2O17VK79dWfIY7UUYz5XTis7Vm40hT87ZW5fu1U69NJ7bTqvd8p1WmjyLFVb0ja54fwd53gHkqT24WPgh4cIJ1G6U/AAzlMuq31WgbXdGKvjqOIJR6m4M71Arjfd537NwF2OddzalSoHgzMC+PRp5NMThlbV2etcCy/Ph0YH3uO7gWOCOPWD4amAfc1KjMvM43chnkMq+ZZBtdMT4ia3hokGB6I7L8kZF6m4M71ArOT6npqEzTc+4P+A5gHTAAXJpSui0izgM2pJTWAp8GPhcRWyiu8p2R170tIq4EbgeeBN6eUtoJUFZm3uR7gMsj4sPAxlw2jbbRTUsXDO/TKCzvJyz1tvHPv7e0076wdUnTUZlAESCldC1wbV3a+2se/xx4fYN1zwfOb6bMnH4ncHxJesNt9Cp/ZKTet68njNJ4n/eydKmRSgWKah9/ZCRpZrN1SdNhoChJ0gxg65Kmw0BRkmaINRtH2HjPKDt27uKkVesNEmYgW5c0VVUa9SxJahOnRpE0HQaKkjQDODWKpOkwUJSkGcCpUSRNh4GitA/G+3zdeNeDnLRqvc14qiwn3pc0HQaK0jTZ50u9xLu7SJoOA0VpmuzzpV7SituBSpp5nB5Hmib7fKnXODWKpKnyiqI0Tfb5kiT1OwNFaZrs8yVJ6nc2PUvT5O2wJEn9zkBR2gf2+ZIk9TObniU1zXkjJWlmMVCU1BTnjZSkmcdAUVJTnDdSkmYeA0VJTXHeSEmaeQwUJTXFeSMlaeYxUJTUlKrNG+nAGklqP6fHkdSUKs0b2WhgTW09JUn7zkBRUtOqMm/kRANrqlA/SeoXNj1L6jkOrJGkzjBQlNRzHFijmch+ueoGA0VJPadqA2ukdnPCe3WLgaLUgGfv1bV0wTAXnHYcw0ODBDA8NMgFpx1n/0T1rMm+b5zwXt3iYBaphKNqp2/8B2/Hzl2ctGp920ZGV2VgjbSvmvm+sV+uusUrilIJz96nx+Yxaeqa+b6xX666xUBRKuHZ+/QYYEtT18z3jf1y1S0Gin3AvnSt59n79BhgS1PXzPeN/XLVLfZR7HH2pWuP5Yvms/LqTXtcHfPsfXJzhgYZKQkKDbClxpr9vrFfrrrBK4o9zqa+9vDsfXp6rXnMq/GqAr9vVGVeUexxNvW1j2fvU1el+0FPxqvx6qTJZgPw+0ZVZaDY42zqU9X0yg+e94tWp3hSol5m03OP67WmPqkqvBqvTrGLkHqZgWKPs2+LND2ObFeneFKiXlaJQDEiDomI6yLijvz/4Ab5luU8d0TEspr0l0bEpojYEhEfjYiYqNwofDTnvzUiXlJT1lcjYjQi/rnd+90qSxcM860Vp3DXqt/iWytOMUiUmuDVeHWKJyXqZZUIFIEVwNdTSvOAr+fne4iIQ4APACcAxwMfqAkoPwGcA8zLf4snKffUmrzn5PXHrQbe3LI9k1RJXo1Xp3hSol5WlcEsS4CT8+PLgOuB99TlWQRcl1J6ECAirgMWR8T1wDNSSt/O6Z8FlgJfmaDcJcBnU0oJuCEihiLiOSmle1NKX4+Ik5HU93pl4I16Wy/NBiDVq0qg+KyU0r0AKaV7I+KZJXmGga01z7fltOH8uD59onIblXXvvu6IJEn1PClRr+pYoBgR/wI8u2TRe5stoiQtTZA+nbKaFhHnUDRbM3fu3KmsKk3ZZHOwSZLUDh0LFFNKr2q0LCJ+PN70GxHPAe4rybaN3c3IAEdQNCVvy49r07fnx43K3QYc2WCdpqSULgEuAVi4cOGUgkxpKpyDTZLULVUZzLIWGB/FvAy4piTPOuA1EXFwHsTyGmBdblp+JCJOzKOdz6xZv1G5a4Ez8+jnE4GHx5uotW+8JVrrOQebJKlbqhIorgJeHRF3AK/Oz4mIhRHxKYA8iOVDwM3577zxgS3AHwCfArYAP6AYyNKwXOBa4M6c/++APxyvSET8O/APwCsjYltELGrLHvehRle+aoNFA8mpcw42SVK3VGIwS0rpJ8ArS9I3AG+teX4pcGmDfC+cQrkJeHuDuvz6VOqu3Sa7JZpNqNPjbRolSd1SlSuK6gOTXfmyCXV6nINNkvZmC1VnGCiqZSa7+4BNqNPjxNDT44+I1L+a6eqk1qhE07P6w/JF81l59aY9rhrWXvmyCXX6nINtauzmIPW3ybo6qXW8oqiWmezKl02o6hS7OUj9zRaqzvGKolpqoitf3saq+vplYm9/RKT+ZgtV5xgoqqNsQq2ufmqu7eSPSL8E11Ivmayrk1rHpmdJQH8113aqm4Md6qXucJBf53hFURLQX821nermYId6qXtsoeoMA0VJQP/1+enEj0g/BdeSVMamZ0mAo9KnY7K5QyWp1xkoSgLs8zMdBteS+p1Nz5KeYp+fqXHKJ0n9zkBRkvaBwbWkfmbTsyRJkkoZKEqSJKmUgaIkSZJKGShKkiSplIGiJKmlxu9/feNdD3LSqvU9fUvDftoXaToMFCVJLdNP97/up32RpstAUZLUMhPd/7rX9NO+aPpm+lVlA0VJUsv00/2v+2lfND1eVTZQlCS1UD/d/7qf9kXT41VlA0VJUgv10/2v+2lfND1eVfYWfpKkFuqn+1/3075oeuYMDTJSEhTOpKvKBoqSKme88/iOnbs4adV6f5x7TD/d/7qf9kVTt3zRfFZevWmP5ueZdlXZQFFSpTTqPA74gy2po7yqbKAoqWIm6jw+k76cJVXDTL+q7GAWSZVi53FJqg4DRUmV4pQkklQdBoqSKsUpSSSpOuyjKKlS7DwuSdVhoCipcmZ653FJqgqbniVJklTKQFGSJEmlDBQlSZJUykBRkiRJpQwUJUmSVKoSgWJEHBIR10XEHfn/wQ3yLct57oiIZTXpL42ITRGxJSI+GhExUblR+GjOf2tEvCSnvzgivh0Rt+X03+nE/kuSJFVRJQJFYAXw9ZTSPODr+fkeIuIQ4APACcDxwAdqAspPAOcA8/Lf4knKPbUm7zl5fYDHgDNTSi/IZXwkIoZauJ+SesiajSNsvGeUG+96kJNWrWfNxpFuV0mSOqoqgeIS4LL8+DJgaUmeRcB1KaUHU0oPAdcBiyPiOcAzUkrfTikl4LM16zcqdwnw2VS4ARiKiOeklP4rpXQHQEppO3AfcHhL91RST1izcYSVV29ix85dAIyMjrHy6k0Gi5JmlKoEis9KKd0LkP8/syTPMLC15vm2nDacH9enT1Ruo7KeEhHHA7OBH5RVOCLOiYgNEbHh/vvvn3QHJfWW1es2M/bEzj3Sxp7Yyep1m7tUI0nqvI7dmSUi/gV4dsmi9zZbRElamiB9OmUVC4urlJ8DlqWUdpUVkFK6BLgEYOHChZNtT1KP2T46NqV0SepHHQsUU0qvarQsIn6cm37vzUHafSXZtgEn1zw/Arg+px9Rl749P25U7jbgyLJ1IuIZwJeB9+VmaUkz0JyhQUZKgsI5Q4Mt39Z4X8gdO3dx0qr13ttaUmVUpel5LTA+inkZcE1JnnXAayLi4DyI5TXAutyk/EhEnJhHO59Zs36jctcCZ+bRzycCD+dgcjbwTxT9F/+hxfsoqYcsXzSfwVkDe6QNzhpg+aL5Ld2OfSElVVlVAsVVwKsj4g7g1fk5EbEwIj4FkFJ6EPgQcHP+Oy+nAfwB8ClgC0Wfwq9MVC5wLXBnzv93wB/m9P8B/AZwVkTckv9e3J5dllRlSxcMc8FpxzE8NEgAw0ODXHDacS2/0mdfSElV1rGm54mklH4CvLIkfQPw1prnlwKXNsj3wimUm4C3l6R/Hvj8FKsvqU8tXTDc9iZg+0JKqrKqXFGUpBmpUZ/HdvSFlKSpMlCUpC7qVF9ISZqOSjQ9S9JMNd60vXrdZraPjjFnaNBRz5Iqw0BRkrqsE30hJWk6bHqWJElSKQNFSVJfGp/I/Ma7HuSkVeudm1KaBgNFSVLfcSJzqTUMFCVJfceJzKXWMFCUJPUdJzKXWsNAUZL6gP3x9uRE5lJrGChKUo+zP97enMhcag0DRUnqcfbH29vSBcNccNpxDA8NEsDw0CAXnHac81VKU+SE23rKeNPVjp27OGnVeu8OIfUI++OVcyJzad95RVFAtZqu7GslTY398SS1i4GigOo0XVUpYJV6hf3xJLWLgaKA6jRdVSVglXqJ/fEktYt9FAUUTVQjJUFhp5uuqhKwSr3G/niS2sErigKq03RlXytJkqrDQFFAdZquqhKwSpIkm55VowpNV+PbX71uM9tHx5gzNOg0PZIkdYmBoiqnCgGrJEmy6VmSJEkNGChKkiSplIGiJEmSShkoSpIkqZSBoiRJ6qg1G0fYeM8oN971ICetWu9tWivMQFGSJHXMmo0jrLx6Ezt27gJgZHSMlVdvMlisKANFSZLUMavXbWbsiZ17pI09sZPV6zZ3qUaaiIGiJEnqmO2jY1NKV3cZKEqSpI6ZMzQ4pXR1l4GiJEnqmOWL5jM4a2CPtMFZAyxfNL9LNdJEvIWfJEnqmPFbtK5et5nto2PMGRpk+aL53rq1ogwUJUlSRy1dMGxg2CNsepYkSVIpA0VJkiSVMlCUJElSqUoEihFxSERcFxF35P8HN8i3LOe5IyKW1aS/NCI2RcSWiPhoRMRE5Ubhozn/rRHxkpz+3Ij4TkTcEhG3RcS5ndh/SZKkKqpEoAisAL6eUpoHfD0/30NEHAJ8ADgBOB74QE1A+QngHGBe/ls8Sbmn1uQ9J68PcC/waymlF+ftrIiIOS3cT0mSpJ5RlUBxCXBZfnwZsLQkzyLgupTSgymlh4DrgMUR8RzgGSmlb6eUEvDZmvUblbsE+Gwq3AAMRcRzUko7UkqP5zwHUJ3jI0mS1HFVCYSelVK6FyD/f2ZJnmFga83zbTltOD+uT5+o3EZlERFHRsStefmFKaXt+7BfkiRJPatj8yhGxL8Azy5Z9N5miyhJSxOkT6csUkpbgRflJuc1EXFVSunHexUQcQ5FszVz586dZHOSJEm9p2OBYkrpVY2WRcSPc9Pvvbkp+b6SbNuAk2ueHwFcn9OPqEsfvwrYqNxtwJEN1hmv7/aIuA34deCqkv25BLgEYOHChZMFppIkST2nKk3Pa4HxUczLgGtK8qwDXhMRB+dBLK8B1uUm5Uci4sQ82vnMmvUblbsWODOPfj4ReDgHk0dExCBA3sZJwOaW7qkkSVKPiGL8R5crEXEocCUwF7gHeH1K6cGIWAicm1J6a873FuDP8mrnp5T+PqcvBD4DDAJfAf5nSilNUG4AH6MYHf0Y8HsppQ0R8WrgL9ndpP2xfOVwsvrfD/ywBYdiMocBD3RgOzONx7U9PK7t4XFtD49re3hc22Nfj+tzU0qHN5OxEoGimhMRG1JKC7tdj37jcW0Pj2t7eFzbw+PaHh7X9ujkca1K07MkSZIqxkBRkiRJpQwUe8uk/SU1LR7X9vC4tofHtT08ru3hcW2Pjh1X+yhKkiSplFcUJUmSVMpAseIi4kMRcWtE3BIRX8t3jCHPAfnRiNiSl7+k23XtJRGxOiK+n4/dP0XEUM2ylfm4bo6IRd2sZ6+JiNdHxG0RsStPW1W7zOO6DyJicT52WyJiRbfr08si4tKIuC8ivluTdkhEXBcRd+T/B3ezjr0m3/72GxHxvfwd8Mc53eO6DyLiFyLipoj43/m4/q+cfnRE3JiP6xURMbtddTBQrL7VKaUXpZReDPwz8P6cfiowL/+dA3yiS/XrVdcBL0wpvQj4L2AlQEQcC5wBvIBins2PR8RA12rZe74LnAb8W22ix3Xf5GN1McXn/ljgDfmYano+Q/E+rLUC+HpKaR7w9fxczXsSeFdK6fnAicDb83vU47pvHgdOSSn9CvBiYHG+UciFwEX5uD4EnN2uChgoVlxK6ac1Tw9k932slwCfTYUbgKF8m0I1IaX0tZTSk/npDey+DeQS4PKU0uMppbuALcDx3ahjL0opfS+lVHY3I4/rvjke2JJSujOltAO4nOKYahpSSv8GPFiXvAS4LD++DFja0Ur1uJTSvSml/8yPHwG+Bwzjcd0n+Tf+Z/nprPyXgFPYfXvhth5XA8UeEBHnR8RW4HfZfUVxGNhak21bTtPUvYXijj7gcW0Xj+u+8fi137PyLWHJ/5/Z5fr0rIg4ClgA3IjHdZ9FxEBE3ALcR9Ea9gNgtOZiR1u/DwwUKyAi/iUivlvytwQgpfTelNKRwBeAd4yvVlKUQ9hrTHZcc573UjSZfGE8qaQoj2uNZo5r2WolaR7X5nn81BMi4mnAPwJ/UtcipmlKKe3M3c+OoGhdeH5ZtnZtf/92FazmpZRe1WTWLwJfBj5AcQZxZM2yI4DtLa5aT5vsuEbEMuC1wCvT7nmiPK6TmML7tZbHdd94/NrvxxHxnJTSvbkbz33drlCviYhZFEHiF1JKV+dkj2uLpJRGI+J6ij6gQxGxf76q2NbvA68oVlxEzKt5+jrg+/nxWuDMPPr5RODh8cv7mlxELAbeA7wupfRYzaK1wBkRcUBEHE0xWOimbtSxz3hc983NwLw80nE2xcCgtV2uU79ZCyzLj5cB13SxLj0nIgL4NPC9lNJf1SzyuO6DiDh8fFaOiBgEXkXR//MbwOk5W1uPqxNuV1xE/CMwH9gF/BA4N6U0kj+UH6MYufcY8HsppQ3dq2lviYgtwAHAT3LSDSmlc/Oy91L0W3ySovnkK+WlqF5E/Dfgb4DDgVHglpTSorzM47oPIuI3gY8AA8ClKaXzu1ylnhURXwJOBg4DfkzRSrMGuBKYC9wDvD6lVD/gRQ1ExMuBfwc2UfxeAfwZRT9Fj+s0RcSLKAarDFBc3LsypXReRDyPYlDbIcBG4E0ppcfbUgcDRUmSJJWx6VmSJEmlDBQlSZJUykBRkiRJpQwUJUmSVMpAUZIkSaUMFCVJklTKQFGSJEmlDBQlqYIi4u8j4rURMRQRX8mTmUtSRxkoSlI1HUdxd5trgA+llP6py/WRNAN5ZxZJqpiI2A94hOIWkxenlC7scpUkzVBeUZSk6pkHbAfOAs6NiFndrY6kmcpAUZKq5zjgupTSeuC7wJldro+kGcpAUZKq5ziKABHgz4GVEbF/F+sjaYayj6IkSZJKeUVRkiRJpQwUJUmSVMpAUZIkSaUMFCVJklTKQFGSJEmlDBQlSZJUykBRkiRJpQwUJUmSVOr/AMnLBttqNCuuAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 720x432 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",
    "w = 8  # wordlength of the quantized signal\n",
    "xmin = -1  # mimimum amplitude of input signal\n",
    "N = 8192  # number of samples\n",
    "K = 30  # maximum lag for cross-correlation\n",
    "\n",
    "\n",
    "def uniform_midtread_quantizer(x, Q):\n",
    "    # limiter\n",
    "    x = np.copy(x)\n",
    "    idx = np.where(x <= -1)\n",
    "    x[idx] = -1\n",
    "    idx = np.where(x > 1 - Q)\n",
    "    x[idx] = 1 - Q\n",
    "    # linear uniform quantization\n",
    "    xQ = Q * np.floor(x/Q + 1/2)\n",
    "    \n",
    "    return xQ\n",
    "\n",
    "def analyze_quantizer(x, e):\n",
    "    # estimated PDF of error signal\n",
    "    pe, bins = np.histogram(e, bins=20, density=True, range=(-Q, Q))\n",
    "    # estimate cross-correlation between input and error\n",
    "    ccf = 1/len(x) * np.correlate(x, e, mode='full')\n",
    "    # estimate PSD of error signal\n",
    "    nf, Pee = sig.welch(e, nperseg=128)\n",
    "    # estimate SNR\n",
    "    SNR = 10*np.log10((np.var(x)/np.var(e)))\n",
    "    print('SNR = %f in dB' %SNR)\n",
    "\n",
    "    # plot statistical properties of error signal\n",
    "    plt.figure(figsize=(9,4))\n",
    "\n",
    "    plt.subplot(121)\n",
    "    plt.bar(bins[:-1]/Q, pe*Q, width = 2/len(pe))\n",
    "    plt.title('Estimated histogram of quantization error')\n",
    "    plt.xlabel(r'$\\theta / Q$')\n",
    "    plt.ylabel(r'$\\hat{p}_x(\\theta) / Q$')\n",
    "    plt.axis([-1, 1, 0, 1.2])\n",
    "\n",
    "    plt.subplot(122)\n",
    "    plt.plot(nf*2*np.pi, Pee*6/Q**2)\n",
    "    plt.title('Estimated PSD of quantization error')\n",
    "    plt.xlabel(r'$\\Omega$')\n",
    "    plt.ylabel(r'$\\hat{\\Phi}_{ee}(e^{j \\Omega}) / \\sigma_e^2$')\n",
    "    plt.axis([0, np.pi, 0, 2])\n",
    "    plt.tight_layout()\n",
    "    \n",
    "    plt.figure(figsize=(10,6))\n",
    "    ccf = ccf[N-K-1:N+K-1]\n",
    "    kappa = np.arange(-len(ccf)//2,len(ccf)//2)\n",
    "    plt.stem(kappa, ccf)\n",
    "    plt.title('Cross-correlation function between input signal and error')\n",
    "    plt.xlabel(r'$\\kappa$')\n",
    "    plt.ylabel(r'$\\varphi_{xe}[\\kappa]$')\n",
    "\n",
    "\n",
    "# quantization step\n",
    "Q = 1/(2**(w-1))\n",
    "# compute input signal\n",
    "np.random.seed(1)\n",
    "x = np.random.uniform(size=N, low=xmin, high=(-xmin-Q))\n",
    "# quantize signal\n",
    "xQ = uniform_midtread_quantizer(x, Q)\n",
    "e = xQ - x\n",
    "# analyze quantizer\n",
    "analyze_quantizer(x, e)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise**\n",
    "\n",
    "* Change the number of bits `w` and check if the derived SNR holds\n",
    "* How does the SNR change if you lower the magnitude of the minimum amplitude `xmin` of the input signal?\n",
    "* What happens if you chose the magnitude of the minimum amplitude `xmin` in the range of the quantization step? Why?\n",
    "\n",
    "Solution: The numerically computed SNR conforms well to the theoretic result derived above. Lowering the magnitude of the minimum amplitude results in a lower SNR as predicted above. The input signal $x[k]$ is correlated to the quantization error $e[k]$ if the magnitude of the minimum amplitude is lowered such that it is close to the quantization step. Here the assumptions made for the statistical model of the quantization error do not hold."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Harmonic Signal\n",
    "\n",
    "For a harmonic input signal $x[k] = x_\\text{max} \\cdot \\cos[\\Omega_0 k]$ the variance $\\sigma_x^2$ is given by its squared [root mean square](https://en.wikipedia.org/wiki/Root_mean_square) (RMS) value\n",
    "\n",
    "\\begin{equation}\n",
    "\\sigma_x^2 = \\frac{x_\\text{max}^2}{2}\n",
    "\\end{equation}\n",
    "\n",
    "Introducing this into the definition of the SNR together with the variance $\\sigma_e^2$ of the quantization error yields\n",
    "\n",
    "\\begin{equation}\n",
    "SNR = 10 \\cdot \\log_{10} \\left(2^{2 w} \\cdot \\frac{3}{2} \\right) \\approx 6.02 \\, w + 1.76 \\quad \\text{in dB}\n",
    "\\end{equation}\n",
    "\n",
    "The gain of 1.76 dB with respect to the case of a uniformly distributed input signal is due to the fact that the amplitude distribution of a harmonic signal is not uniform\n",
    "\n",
    "\\begin{equation}\n",
    "p_x(\\theta) = \\frac{1}{\\pi \\sqrt{1 - (\\frac{\\theta}{x_\\text{max}})^2}}\n",
    "\\end{equation}\n",
    "\n",
    "for $|\\theta| < x_\\text{max}$. High amplitudes are more likely to occur. The relative power of the quantization error is lower for higher amplitudes which results in an increase of the average SNR."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Normally Distributed Signal\n",
    "\n",
    "So far, we did not consider clipping of the input signal $x[k]$, e.g. by ensuring that its minimum/maximum values do not exceed the limits of the quantizer. However, this cannot always be ensured for practical signals. Moreover, many practical signals cannot be modeled as a uniform distribution. For instance, a [normally distributed](../random_signals/important_distributions.ipynb#Normal-Distribution) random signal exceeds a given maximum value with non-zero probability. Hence, clipping will occur for such an input signal. Clipping results in overload distortions whose amplitude can be much higher that $\\frac{Q}{2}$. For the overall average SNR both granular and overload distortions have to be included.\n",
    "\n",
    "The root mean square (RMS) of the normal distributed input signal is given by its standard deviation $\\sigma_x$. The RMS level $A$ of the input signal normalized to the maximum level of the quantizer as\n",
    "\n",
    "\\begin{equation}\n",
    "A = \\frac{\\sigma_x}{x_\\text{max}}\n",
    "\\end{equation}\n",
    "\n",
    "The probability that clipping occurs can be derived from the [cumulative distribution function](../random_signals/important_distributions.ipynb#Normal-Distribution) (CDF) of the normal distribution as\n",
    "\n",
    "\\begin{equation}\n",
    "\\Pr \\{ |x[k]| > x_\\text{max} \\} = 1 + \\text{erf} \\left( \\frac{-1}{\\sqrt{2} A} \\right)\n",
    "\\end{equation}\n",
    "\n",
    "where $x_\\text{max} = - x_\\text{min}$ was assumed. For a normally distributed signal with a given probability that clipping occurs $\\Pr \\{ |x[k]| > x_\\text{max} \\} = 10^{-5}$ the SNR can be approximately calculated to [[Zölzer](../index.ipynb#Literature)]\n",
    "\n",
    "\\begin{equation}\n",
    "SNR \\approx 6.02 \\, w - 8.5 \\quad \\text{in dB}\n",
    "\\end{equation}\n",
    "\n",
    "The reduction of the SNR by 8.5 dB results from the fact that small signal values are more likely to occur for a normally distributed signal. The relative quantization error for small signals is higher, which results in a lower average SNR. Overload distortions due to clipping result in a further reduction of the average SNR."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Example - Quantization of a normal distributed signal\n",
    "\n",
    "The following example evaluates the SNR of a linear uniform quantizer with $w=8$ for a normally distributed signal $x[k]$. The SNR is computed and plotted for various RMS levels, the probabilities for clipping are shown additionally."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Maximum SNR = 40.854 dB for A = -11.7 dB with clipping probability 1.2e-04\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfIAAAEOCAYAAABywZWyAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xd8VFX6x/HPkx4IvYTQO9JLgohIEzsqFmyL3V3U1dVdu2tDXXft2LvYVo2uDUQQUQKIIr2GTui9Qyip5/fHDPuLbAgBZnJnJt/36zWvzNy5c+/z5EKeOefee4455xAREZHwFOV1ACIiInLsVMhFRETCmAq5iIhIGFMhFxERCWMq5CIiImFMhVxERCSMqZCLiIiEMRVyERGRMKZCLiIiEsZUyEVERMJYjNcBlEbNmjVd48aNA7a9vXv3UrFixYBtz0vKJfRESh6gXEJVpOQSKXlA4HOZMWPGVudcrdKsGxaFvHHjxkyfPj1g2xs/fjx9+vQJ2Pa8pFxCT6TkAcolVEVKLpGSBwQ+FzNbVdp11bUuIiISxlTIRUREwpgKuYiISBhTIRcREQljKuQiIiJhTIVcREQkjKmQi4iUgdz8Qhas3+11GBKBVMhFRMrAfV/O5ZyXfmbN9n1ehyIRRoVcRCTAsnPy+XbOerbsyeEfIxfw9PeL+GrWOgDGZG7kQF4Bs1bv4PxXJrFrf57H0Uq4C/rIbmYWDUwH1jnnzjWzJkA6UB2YCVzlnMsNdhwiIsG2fW8uP2Ru5MPJq1iwYTex0Uahg4JCR6WEGLJz8vl4ympeH7+cbXt9f/YmLtnCeR3rehy5hLOyGKL1dmAhUNn/+ilgqHMu3czeAG4AXi+DOEREAi47J5/v5q5n1bZ9pE9bw/a9udStksA13RuxZFM2D/RvTe3K8VROiOU/09fw0PDM331+1uqdKuRyXIJayM2sPtAfeAK4w8wMOBX4g3+VD4AhqJCLSJjZl5vPh5NX8eaE5ezYl0d0lNGuXhXev64r7etVwffn7vcGdWvEpGVb2bQ7h9lrdgIwbeX2sg5dIkywW+QvAPcAlfyvawA7nXP5/tdrgXpBjkFEJGDyCgr54NeVvDFhOVuzc+nVsha392tBpwZViY763+JdVFSU8eZVaQDMWLWDCYs380rGMpZvyaZ+tUTiY6LLIgWJMOacC86Gzc4FznHO/dnM+gB3AdcBk51zzf3rNABGOefaF/P5wcBggOTk5NT09PSAxZadnU1SUlLAtucl5RJ6IiUPUC6H2ryvkDfm5JC1q5C2NaK4oHkcLaode/HN2lnAY78dAKB73Whu7JBQqs9FynGJlDwg8Ln07dt3hnMurVQrO+eC8gD+ha/FvRLYCOwDPga2AjH+dboDY460rdTUVBdIGRkZAd2el5RL6ImUPJxTLget3rbXXf/eVNfo3pGu3cPfu5Fz1gckpoKCQtfo3pH/fezPzS/V5yLluERKHs4FPhdguitlvQ3a7WfOufudc/Wdc42By4FxzrlBQAYw0L/aNcDwYMUgInI8DuQV8OKPSznt+QlMztrGbf1a8P3fetG/Q0pAth8VZTw+oC09W9QE4KeFmwOyXSlfyuKq9UPdC6Sb2T+AWcC7HsQgInJYzjl+XLiZx0Zmsmb7fs7tkMID/VuTUiUx4Pu6qntj/tCtESc/+RNfz1obsC8JUn6USSF3zo0HxvufZwEnlsV+RUSO1sqte3n020wyFm+hRe0kPvljN05uXjOo+4yOMgZ0qsewSSvYlp1DjaT4oO5PIotGdhMRAfbnFvDsmMWcMXQi01bu4MH+rRl1e8+gF/GDLuxcj/xCx3fzNpTJ/iRyeNG1LiISMpxzfD9/I//4biHrdu7nws71uP/sE6hduXRXkAdK65TKnFCnEl/OWMtVJzUq9j50keKoRS4i5dayzdlcPWwqN388k0oJMXx+Y3eGXtapzIv4QYO6NWTO2l3c++Xcg3f/iByRWuQiUu5k5+Tz8k9LeXfSChLjonn0/LYM6taQmGhv2zZXntSIjbsP8GrGcjo1qMYfujX0NB4JDyrkIlJuOOcYMWc9/xy1kE27c7g0rT73nHUCNUPk4jIz487TWzElazuvT1jGpWn1yc7Jp2qFOK9DkxCmQi4i5cLaPYVc/tZvTFmxnXb1KvP6lal0aVjN67D+R1SUcVX3RtyePpvmD4wmJsqY+sBpVK+oYi7FUyEXkYi2a38eL/64lPd/3U/lxHyeuLAdl3dteMRx0b10Zts6pDaqxoxVO8gvdMxZu5O+rWp7HZaEKBVyEYlIhYWOz6ev4Zkxi9m+L5fe9WIYel0fqoVByzYhNpovbz6ZPQfyaD/kB+at3aVCLoelQi4iEWfe2l08+M085qzdRdfG1fjgvBPZunRWWBTxoiolxNK0ZkXmrdvldSgSwlTIRSRi7NqXx7M/LObfU1ZRo2I8L1zWiQGd6mJmjF/qdXTHpl29KvyWtQ3nnO4tl2KpkItI2CssdHw5cy1Pjl7Ejn25XNO9MXec0ZLKCbFeh3bcejSvwYg561m8aQ8n1KnsdTgSglTIRSSsZa7fxcPDM5mxagddGlblwxtOpG3dKl6HFTB9/OfGMxZtUSGXYqmQi0hY2pqdw3M/LCF92mqqV4jjmYEduLhLfaJC+Gr0Y5FcOYHWKZX5eekWbu7TzOtwJASpkItIWDnYjf7YyAXszy3g+h5NuK1fC6okhn83+uF0bVyNL2espaDQhfRtc+INFXIRCRszV+/g4eHzmb9uN2mNqvHkxR1oXjvJ67CCrkvDanw4eRWLN+6hTV11r8vvqZCLSMjblp3D098v5rPpa0iuHM/QyzoyoGO9iOtGP5yDI9D97bPZfDr4JI+jkVCjQi4iIaug0PHp1NU8M2Yxe3PyGdyrKbf1a0FSfPn609WgeiLdmlRnyortfD1rHTpTLkWVr/8NIhI2Zq/ZyUPfzGfeul2c1LQ6jw1oR8vkSl6H5Qkz47Mbu9Pz6XFMydpGM02KJkWokItISNm+N5dnxiwifdoaaiXF8+LlnTi/Y10NhgJ0a1KDHxdu4vIGkXthnxw9FXIRCQkFhY7Ppq3h6TGL2HMgnz+e4rsavVIEDOoSKN2b1uCLGWvJ2hnNqV4HIyFDhVxEPDdnzU4eHj6fOWt30a2Jrxu9VZ3y2Y1ekjPaJlNxeDQZa/L5o9fBSMgIWiE3swRgIhDv388XzrlHzOx9oDdwcBaAa51zs4MVh4iErg279vPw8EzGLthErUrqRj+SSgmxnN+pHl9OX617yuW/gtkizwFOdc5lm1ksMMnMRvvfu9s590UQ9y0iIWxfbj5vT1zBmxOX4xz89bQW3HBKE3Wjl0LH+lX4dCqs37mfBtUreB2OhICgFXLnnAOy/S9j/Q8XrP2JSHgYPW8DQ77NZNPuHM5pX4f7zmpNwxoqSKXVpGZFALK27lUhFwCigrlxM4s2s9nAZmCsc26K/60nzGyumQ01s/hgxiAioWHz7gPc+NF0bv54JjWT4vnipu68NihVRfwoNanlK+QrtmQfYU0pL8zXcA7yTsyqAl8DfwG2ARuBOOAtYLlz7rFiPjMYGAyQnJycmp6eHrB4srOzSUqKjGEdlUvoiZQ8IDC55BQ4flqdx7fL88gvhAuax3JW49gyP78bKcfFOcfNP+6lR71YrmoT3u2gSDkmEPhc+vbtO8M5l1aadcvkqnXn3E4zGw+c5Zx71r84x8zeA+46zGfewlfoSUtLc3369AlYPOPHjyeQ2/OScgk9kZIHHH8u4xZt4p/fLmDVtjx6tqjJYwPa/bdruKxF0nGpM3k0O60SvXt3D+sLAyPpmHiZS9C61s2slr8ljpklAqcBi8wsxb/MgAuA+cGKQUS8sXrbPq57byrXvz+d6Cjjkz9146MbunlWxCNNt5QYZqzawRcz1nodioSAYLbIU4APzCwa3xeGz51zI81snJnVAgyYDdwUxBhEpAzlFxTyydTVPP39YgAe7N+aq7s3Ji4mqJfjlDtnNo5h9s54vpixlkvSGngdjngsmFetzwU6F7NcAxKJRKCpK7bz8PD5LNq4hx7Na/DUxR2oX00XsgVDlBl9WtVi2C8r2JebT4U4je1Vnulrsogcl427DnDbp7O49M3J7DmQz+uDuvDvG7qpiAdZj+Y1yStwTFyyxetQxGP6GicixyQnv4Bhk1by8ril5Bc6bju1OTf3aU5iXLTXoZULJzapTqMaFbj7P3Pp2KAqKVUSvQ5JPKIWuYgctYzFmznrhZ956vtFnNysJj/+rTd3nNFKRbwMJcRG8+ZVqezJyWfcos1ehyMeUotcREpt1ba9PD5yAT8u3EzTmhV5/7qu9GlV2+uwyq1WyZVIrhzP5OXbGNStkdfhiEdUyEXkiPbnFvDa+GW8OTGLmCjjvrNP4PoeTXQ1usfMjO5NazBp2VYKCx1RmkSlXFIhF5HDcs7x3dwNPPHdAtbvOsCATnW5/+zW1KmS4HVo4tenVW2+mb2eaSu3061pDa/DEQ+okItIsZZs2sPT0w6wcPtMWqdU5oXLO3Nik+pehyWHOKNtMhXjovly5loV8nJKhVxEfmdfbj5Dxy5h2C8rSYh2PD6gLVec2JCYaHWjh6IKcTGc0z6FUfM28uj57XTBYTmk/5ki8l8Zizdz5gsTefvnFVyaVp8ne1bgqu6NVcRD3MWp9cnOyWdM5kavQxEP6H+niLBi616uf38a1703jdjoKD4bfBL/uqgDleJ08VQ4OLFxdepVTWTk3A1ehyIeUNe6SDm2a18er45fxnu/rCA+JpoHzmnNNSdrbPRwExVl9GpZk5FzN1BQ6Mp8iljxlgq5SDlUWOj4dJpvcpPdB/IY2KU+d5/VitqVdDV6uDqpaQ0+nbqGBet3075+Fa/DkTKkQi5SzizauJsHv57P9FU76N60Bg+f14bWKZW9DkuOU3f/FeuTs7aqkJczKuQi5UR2Tj7/GLmA9GlrqFohlmcGdmBgan3M1A0bCWpXTqBprYr8lrWdwb2aeR2OlCEVcpEI55xj+Oz1/HPUQrZk5zC4V1Nu7t2MahXjvA5NAqx70xoMn72e/IJC3WlQjqiQi0SwBet3M2REJlNXbqdD/Sq8eVUqnRtW8zosCZLuzWrw8ZTVzFy9U4P3lCMq5CIRaPeBPJ7/YQkfTl5JlcRYnryoPZemNdBY3BGud8ta1KgYxzNjFvH5jd112qScUCEXiSCFhY5Ppq7mhR+Xsm1vDld2a8SdZ7SkagV1o5cHlRJiufOMVvz963lMWraVni1qeR2SlAGdRBGJEEs27eHyt37jwW/m07RWRYbf0oPHL2inIl7OXJxaj5pJ8bzz8wqvQ5Eyoha5SJjbuS+Xp8cs5tOpq0mKj+HZSzpycZd66lYtp+Jjorm8awNeG7+MHXtzdVFjOaAWuUiYKix0fD59Dac+N4HPpq3h2pMbM+HuvrqlTDitTTKFDiYs2eJ1KFIGgtYiN7MEYCIQ79/PF865R8ysCZAOVAdmAlc553KDFYdIJFq4YTcPfjOfGat2kNaoGo9f0E6Dush/dahXhZpJcYxbtJkLOtfzOhwJsmB2recApzrnss0sFphkZqOBO4Chzrl0M3sDuAF4PYhxiESM7XtzGTp2CZ9MXU2VxFieHtiBgV3q62p0+Z2oKKNPq9qMXbBJ95SXA0E7us4n2/8y1v9wwKnAF/7lHwAXBCsGkUiRm1/IOz9n0fuZDD6ZuppB3Roy7s7euqVMDuvUE2qza38eM1fv9DoUCbKgXuxmZtHADKA58CqwHNjpnMv3r7IWUL+PyGE45/hx4Wb+OWohK7bupXfLWjzYvzUtkit5HZqEuJ4tahIbbXw9a50Gh4lw5pwL/k7MqgJfAw8D7znnmvuXNwBGOefaF/OZwcBggOTk5NT09PSAxZOdnU1SUlLAtucl5RJ6ApXHzgOFvDM/l/lbC6hb0bj8hDg61CrbG00i5ZhA+czlowU5jFudzxOnJFI3KfS618vjMSmtvn37znDOpZVqZedcmTyAR4C7ga1AjH9Zd2DMkT6bmprqAikjIyOg2/OScgk9x5tHXn6BGzYpy7V75HvX6sFRbtikLJebXxCY4I5SpBwT58pnLlv2HHBN7hvpnhuzKLgBHaPyeExKC5juSllfg/YVzcxq+VvimFkicBqwEMgABvpXuwYYHqwYRMLN5OXb6P/SJB79dgGdGlRl1G09ua5HE2J1sZIcg5pJ8XRqUJUJS7d6HYoEUTD76VKAD/znyaOAz51zI81sAZBuZv8AZgHvBjEGkbCwYdd+/jlqEd/OWU+9qom8eVUqZ7RJ1v3gctx6tazFiz8tZVt2DjWS4r0OR4IgaIXcOTcX6FzM8izgxGDtVySc5OQXMGzSSl4et5T8Qsft/VpwU+9mJMZFex2aRIgz29bhhR+XMnz2eq4/pYnX4UgQaIhWEY9MWLKFR0dkkrV1L6e3Seah/m1oWKOC12FJhGmdUpmODaqSPm011/VorF6eCKQTbyJlbM32fQz+cDrXDJuKA967ritvX52mIi5Bc2lafZZsymbhhj1ehyJBoBa5SBk5kFfA6+OX88aE5USZcc9ZrbjhlCbEx6gbXYLr7HYpPDw8k5Fz19OmrobyjTQq5CJB5pzjhwWbeHzkAtbu2M+5HVJ4oH9rUqokeh2alBPVK8ZxcrMajJy7gbvPbKXu9QijrnWRIFq+JZtr3pvGjR/NoEJcNJ/8qRuv/KGLiriUufM61GX19n3MW7fL61AkwEpskZtZb2CHc26umV0K9MI3zOprzrmcsghQJBztzcnn5XHLeHdSFgkx0Tx0bhuu7t5I94OLZ85sW4cHvpnHd3M30KF+Va/DkQA6bCE3s1eBDkC8mS0BkoDvgZOBYcCgMolQJIw455i8Pp97nxvPpt05XNylPvedfQK1Kun+XfFWlQqxdGtSg3GLNnP/Oa29DkcCqKQWeV/nXBv/vOLrgNrOuQIzexOYWzbhiYSPhRt288iITKauyKF9vSq8NiiV1EbVvA5L5L96t6zFE6MWsm7nfupV1emdSFFSP98BAOfcAWCVc67A/9oBeWUQm0hY2LUvj0eGz6f/Sz+zdNMerm0bxze39FARl5DTp1UtAP4zfY3HkUggldQir21mdwBW5Dn+17WCHplIiCssdHw+fQ1Pj1nMzn25XHlSI+44vSWzp/5KtOYIlxDUIrkS/Tuk8PK4ZfRvn6LpcCNESS3yt4FK+M6NH3x+8PU7wQ9NJHSt3bGPQe9M4b6v5tGsVkVG/qUnjw1oR9UKcV6HJlKixwe0IzrK+HDyKq9DkQA5bIvcOfdoWQYiEg6yc/J5+aelvP/rSmKijCcvas9lXRvovlwJG9UrxnFuhxS+mrmWB/q3JiFWAxKFu5KuWn+ppA86524LfDgiock5x4g563niu4Vs3pPDRZ3rcccZLalfTcOqSvjxFfJ1zFy1g5Ob1/Q6HDlOJZ0jn+H/2QNoA3zmf31JkfdEIl7m+l0MGZHJtJU7aF+vCm9elUrnhrqQTcJX18bViY4yflm+VYU8ApTUtf4BgJldi+9WtDz/6zeAH8okOhEPZefk8+yYxXw4eSVVK8Txr4vac2laA13IJmGvUkIsHetX4ZtZ6zm3Q11ap2j89XBWmmGm6uK7yO2gJP8ykYjknOP7+Rs496Wf+WDySq48qREZd/bhihMbqohLxLjj9Fbsy83nkRGZXocix6k0k6Y8Ccwyswz/697AkKBFJOKhRRt389ToRWQs3kLTmhX5bHB3TmxS3euwRALulBY1ubRrA4ZNWsG+3HwqxGkOrXB1xCPnnHvPzEYD3fyL7nPObQxuWCJla/eBPJ7/YQkfTl5JxfgYHuzfmut6NFELXCJaj2Y1eXNCFlNXbKdPq9pehyPHqFRfwfyFe3iQYxEpcwWFjg8nr+TVjGVs25vLoG4NueuMVrofXMqFro2rUyEumn//tkqFPIypL0XKrekrt/PYyAXMXbuLk5vV4L6zT9CsUFKuJMZFc3u/Fvxr9CImLNlC75YatDMcqZBLubNlTw4vj1vKJ1NWk1w5gaGXdeSCTvU0qIuUS9f1aMIHv67k1XHLVMjDVKkKuZlFA8lF13fOrQ5WUCLBkFdQyAe/ruTFH5dyIL+AS9Lqc/85ramcEOt1aCKeiYuJ4o89m/LYyAVkrt9F27pVvA5JjtIRbz8zs78Am4CxwHf+x8hSfK6BmWWY2UIzyzSz2/3Lh5jZOjOb7X+cc5w5iBzRr8u3cs6LP/OP7xbSpVE1vv9rL/51UQcVcRHggs71iI4yRs7d4HUocgxK0yK/HWjlnNt2lNvOB+50zs00s0rADDMb639vqHPu2aPcnshR27znAE98t5Dhs9fToHoib1+dxmmta6sbXaSI6hXjOLlZDb6euY4bezXVxZ5hpjQDwqwBdh3thp1zG5xzM/3P9wALgXpHux2RY1FY6Ph4yir6PTeB0fM2clu/Foz9W29Ob5OsIi5SjL+e1pLte3N58Jv5XociR8mccyWvYPYu0Apfl3rOweXOuedLvROzxsBEoB1wB3AtsBuYjq/VvqOYzwwGBgMkJyenpqenl3Z3R5SdnU1SUlLAtucl5fK/1uwp5IPMHJbtLKR19SiubhNPSlJpvrMGho5JaFIuR/ZBZg6T1+fzSr8KxJTBGAo6JofXt2/fGc65tNKsW5qu9dX+R5z/cVTMLAn4Evirc263mb0OPA44/8/ngOsP/Zxz7i3gLYC0tDTXp0+fo931YY0fP55Abs9LyuX/rdm+jydHL2L0/A1UrRDHc5e05qIuZX81uo5JaFIuR7avxgYyPp5JtWYdSW0U/BENdUwCozQjux3zvORmFouviH/snPvKv71NRd5/m1JcOCdSkgN5BbwxYTmvj19OdJRxY+9mDO7ZlGoVdZ5P5Gh0b1oDM/h06hq6NKym01BhoqT5yF9wzv3VzL7F13r+Hefc+SVt2Hz/At4FFhbthjezFOfcwUsjLwR0QkaO2Y8LNvHYyAWs3r6P/h1SeLB/a1KqJHodlkhYqlYxjj/1bMpbE7Po2aImAzrpsqZwUFKL/CP/z2O9urwHcBUwz8xm+5f9HbjCzDrh+3KwErjxGLcv5diqbXt59NsFjFu0mea1k/j4j93ooXmVRY7bfWedwFcz1/LTws0q5GGipPnIZ/h/TjiWDTvnJgHF9cuMOpbtiQDszy3g9fHLeGNiFrFRxgPntObaHo2JjS67i9lEIllUlNGrZS3GLdpMQaHTxEFhQEO0SlhwzvHDgk08PnIBa3fsZ0Cnuvz9nNYkV07wOjSRiNO3VW2+mrmOiUu20PcETaYS6lTIJeSt2LqXISMymbBkC62SK5E++CROalrD67BEItaZbetQr2oiL/60lD6taumitxB3TIXczBo551YFOhiRovbl5vNqxjLenriCuJgoHjq3DVd3b6RudJEgi4uJ4pa+zfn71/OYsGSLpjgNcSX+RTSz7mY20Mxq+193MLNPgEllEp2US845Rs/bwGnPTeDVjOWc2yGFcXf15oZTmqiIi5SRgan1qVc1kaE/LuVIA4eJtw77V9HMngGGARcD35nZI/gmTpkCtCib8KS8Wb4lm6uHTeXmj2dSOTGW/9zUnecv60TtSjoXLlKW4mKiuPXU5sxZs5PxS7Z4HY6UoKSu9f5AZ+fcATOrBqwHOjjnlpZNaFKe7M3J5/PFuYwdO5GEmGiGnNeGK09qRIxa4CKeubhLfV4Zt4wXflxKn5Y6Vx6qSvorud85dwDAPxb6YhVxCTTnHCPnrue05ycwakUeAzrVY9xdfbi2RxMVcRGPqVUeHkpqkTczsxFFXjcu+vpII7uJHMmyzXt4ZEQmvyzbRpuUylx/Avzpwo5ehyUiRVzcpT4v/riUd39eQV9d9BaSSirkAw55/VwwA5HyIzsnn5d+WsqwSSuoEBfN4wPa8odujfh54jGNPSQiQRQXE8XVJzfi6e8Xs3TTHlokV/I6JDlESSO76a+qBJRzjhFz1vPPUQvZtDuHy9IacM9ZraiRFO91aCJSgktSG/DsmMV8O2c9d5zRyutw5BAlTZqSQTGTpfg551y/4IQkkWjxxj08MmI+v2Vtp329KrxxZSqdG1bzOiwRKYValeI5qWkNRsxZz/WnNKFqBc0sGEpK6lq/q5hlJwH3AJuDE45Emj0H8njhx6W8/+tKKiXE8MSF7bi8a0ON3ywSZq48qRF//ngml7/1G6Nv76kr2EPIESdNATCz3sBDQDxwk3NudBnEJmHMOcc3s9fxz1GL2Jqdw+VdG3LPma00R7hImDqnfQpDzmvDkG8XsGRTNq3q6Fx5qChxiFYzOxNfAT8APOGcyyiTqCSsLdywm0eGZzJ15XY61q/CO1en0bFBVa/DEpHjdE77FIZ8u4AxmRtVyENISefIpwG1gGeAyf5lXQ6+75ybGfToJKzs2p/H0LFL+Oi3VVROiOHJi9pzaVoDotSNLhIRaldO4KSm1Xn/15UM6tZQF6qGiJJa5HuBbGAgvmFai/41dsCpQYxLwkhhoeOrWet4cvRCtu3NZVC3htx1RitdECMSgR4b0I7+L/3Mc2OX8M8L23sdjlDyOfI+ZRiHhKnM9bt4eHgmM1btoHPDqrx/3Ym0q1fF67BEJEhaJlfikrQGfDF9Lbed2oI6VTQPgtdKmjSlq5nVKfL6ajMbbmYvmVn1sglPQtWufXk8PHw+5708iZVb9/L0wA58edPJKuIi5cDNvZtR4Bxv/5zldShCyWOtvwnkAphZL+BJ4ENgF/BW8EOTUFRY6Ph82hpOfW48//5tFVed1Ihxd/bRuXCRcqRB9QoM6FSXj6esYlt2jtfhlHslFfJo59x2//PLgLecc1865x4Cmgc/NAk189bu4qLXf+WeL+fSpGZFvv3LKTw6oB1VKsR6HZqIlLE/92lOTn4hw35Z4XUo5V5JF7tFm1mMcy4f6AcMLuXnJMJs2ZPD82MXkz5tDTUqxvPcJR25qEs9DQghUo41r53EOe1SGDZpJf3b16VN3cpeh1RuldQi/xSYYGbDgf3AzwBm1hxf93qJzKyBmWWY2UIzyzSz2/3Lq5vZWDNb6v+pcTpDVE5+AW9MWE7fZ8fzn+lrue7kJoxtWnUHAAAZy0lEQVS7qzcXp9ZXERcRHjmvDVUSY7nnyzk4d7gRvSXYSrpq/Qkz+wlIAX5w/3+UooC/lGLb+cCdzrmZZlYJmGFmY4FrgZ+cc0+a2X3AfcC9x5OEBN60ldv5+1fzWLo5m34n1Obv/VvTrFaS12GJSAipXTmBW05tzkPfzGfu2l0a+MkjJbXIcc795pz72jm3t8iyJaUZDMY5t+Hges65PcBCoB6+6VE/8K/2AXDBsQYvgbdh135u+mgGl7wxmZ3783jv2q68e21XFXERKdYFnepSMS6aVzKWeR1KuVUm57rNrDHQGZgCJDvnNoCv2JuZZqoPAfkFhbz/60qGjl1CfqHj7jNbcX2PJiTGRXsdmoiEsEoJsfylXwueHL2IcYs2ceoJyV6HVO5YsM9rmFkSMAHfWO1fmdlO51zVIu/vcM79z3lyMxuM/wK75OTk1PT09IDFlJ2dTVJSZLQwjzcX5xzzthbwxZI8Vu8ppEPNaK5sE0ftCiV21gRFpByXSMkDlEuoCrVc8gsdD07aT1QUPH5yYqlnNwy1PI5HoHPp27fvDOdcWqlWds4F7QHEAmOAO4osWwyk+J+nAIuPtJ3U1FQXSBkZGQHdnpeOJ5dNu/e7P30wzTW6d6Tr/s8f3Xdz17vCwsLABXeUIuW4REoezimXUBWKuYyet941unek+3TKqlJ/JhTzOFaBzgWY7kpZa4PWtW6+y5rfBRY6554v8tYI4Bp8A8xcAwwPVgxSvJz8Aj6ftobnxi5hX24B9599Atef0oTY6LJvhYtIZDizbR26NKzK82OXcF7HulSM113KZSWYf7l7AFcBp5rZbP/jHHwF/HQzWwqc7n8tZeSXZVu54NVfeWh4Js1rJTHqtp7c2LuZiriIHBcz44H+bdi8J4dXdeFbmQraVybn3CR+P2NaUf2CtV8p3tbsHIaMyGTk3A0kV47n7avTOK11bd0PLiIBk9qoGhd1qcebE7Po17o2qY00LUdZUDMswuUVFJI+dTWnPT+BHzI3cefpLZl4T19Ob5OsIi4iATfk/LYkV4rnX6MWeR1KuaGTGBHs+/kbGTIik427D5DWqBr/uqg9LZIreR2WiESwygmxXNejCU+MWsjCDbtpnaKhW4NNLfIItG7nfv788Qxu+vcMqleMY9i1afznpu4q4iJSJgam1qdiXDSPjMikoFBDtwabWuQR5EBeAe/8nMWrGcspdL5BXQb3aqoL2USkTFWrGMeQ89ty9xdz+XLGWi7t2sDrkCKaCnmEmLhkCw8Nn8+qbfs4u10dHujfmvrVKngdloiUUwNT6/Ph5FW8nLGU8zvVJSFWo0QGi5pqYW77gUL+8uksrh42lWgz/n1DN16/MlVFXEQ8ZWbcd/YJrNm+n3+NWuh1OBFNLfIwlZtfyNs/ZzF04n6iLIe/ntaCm/s0Iz5G33pFJDT0aF6TP57ShHcmreCUFrU4vY3GYQ8GFfIwNG7RJv4xciFZW/eSlhzN0Gt60aC6WuAiEnruPqsVk7O2cc8Xcxh9ey/qVEnwOqSIo671MLJy616uf38a178/HQzeu7Yrt3ZOUBEXkZAVHxPNy1d0Zn9eAf8arS72YFAhDwMH8goYOnYJZwydyJSsbTxwTmu+v70XfU/QDLAiEvqa1krimu6N+XbOepZt3uN1OBFHhTyEOefIWLSZs1/8mRd/WsqZ7eqQcVcf/tSrKXExOnQiEj7+1KspFeNjeOibzIMzYUqAqBqEqGWb93DZW79x3fvTcM7x0Q0n8vIVnaldWeeXRCT81EyK5/6zWzM5axuvjNOkKoGki91CzM59ubwybhkfTl5FxfhoHh/Qlsu6NlQLXETC3hUnNmDqim0M/XEJPVvW8jqciKFCHiLyCwr5etY6nhmzmK3ZOQxMrc9dZ7aidiW1wEUkMpgZj13Qjt+ytnPvF3O5u6O62ANBhTwErN62j799PpsZq3bQOqUyw67tSrt6VbwOS0Qk4ConxPKPC9rxxw+n811WLKed6nVE4U+F3EMH8gp4e2IWr45fRmx0FC9c1okBnepqelERiWintUnmvI51GTF3PX/etEcTOh0nnXj1gHOOEXPW0++5CTw3dgl9W9Xm+7/24oLO9VTERaRceOS8NiTGwD1fziW/oNDrcMKaCnkZy1y/i0vfnMxtn86iSmIsn/zRNzZ6vaqJXocmIlJmaibFM6h1PLNW7+Rvn8/RLWnHQV3rZWRrdg7PjlnMZ9PXUK1CHE9e1J5L0hoQHaUWuIiUT93rxlAxuRHPj11C//YpnNWujtchhSUV8iArKHR8PGUVz4xZzIG8Aq7v0YTbTm1BlQqxXocmIuK5P/dpxog563l85AK6NKqqO3WOgbrWg2jW6h0MeHUSDw/PpEP9Koy+vRcPndtGRVxExC8mOornL+3I9r253PrJLAoK1cV+tIJWyM1smJltNrP5RZYNMbN1Zjbb/zgnWPv30o69udz/1Twuev1XtuzJ4eUrOvPvG7rRvHaS16GJiIScDvWr8tiAtkxdsZ1hk1Z4HU7YCWbX+vvAK8CHhywf6px7Noj79UxhoeM/M9bw5OhF7D6Qzw09mvDX01uSFK8zGCIiJRmYWp8fFmzimTGLSWtcjc4Nq3kdUtgIWovcOTcR2B6s7YeaeWt3MfCNX7n3y3k0r53Ed7edwoPntlERFxEpBTPjyYvaU6dKAte+N43New54HVLY8OIc+a1mNtff9R72X7k27T7AbZ/O4rxXJrFq2z6eu6Qjn9/YnRPqVPY6NBGRsFIjKZ73ruvKvtx8nhq9WLeklZIF8xdlZo2Bkc65dv7XycBWwAGPAynOuesP89nBwGCA5OTk1PT09IDFlZ2dTVLS8Z2vds4xcV0+6YtyyS+Es5rEck6TWBJjyvZ2skDkEioiJZdIyQOUS6iKlFwOl8cXS3IZmZXHeU1jubhlnAeRHb1AH5O+ffvOcM6llWpl51zQHkBjYP7RvnfoIzU11QVSRkbGcX1+StY2d8Grk1yje0e6S9741WVtyQ5MYMfgeHMJJZGSS6Tk4ZxyCVWRksvh8igoKHT3fjHHNbp3pBs5Z33ZBnWMAn1MgOmulLW2TE/gmlmKc26D/+WFwPyS1g81m/cc4MlRi/hq1jpSqiTw1MXtuSS1AVEa1EVEJGCioozHBrRj6eZs7v5iDs1rJ9GqjsZjP5ygFXIz+xToA9Q0s7XAI0AfM+uEr2t9JXBjsPYfSPkFhXz02yqe/2EJOfmF3Nq3Obf0bU5iXLTXoYmIRKS4mCheG9SFc1+exE3/nsE3t/SgSqLG4ChO0Aq5c+6KYha/G6z9Bcv0ldt58Jv5LNq4h54tavLo+W1pWiv8z02JiIS65MoJvDaoC1e89Rt3fDabt65O07DWxdDIboexZU8Od34+h4FvTGb3/jzeuLILH15/ooq4iEgZ6tq4Oo+c14afFm3mqe8XeR1OSNJNzofILyjk4ymrefYH39jof+7TjFtPbU6FOP2qRES8cFX3xizdnM1bE7NoUTuJS9IaeB1SSFF1KmLGqu089E0mCzbs5pTmNXl0QFuaqQUuIuK5h89tw7LN2TwyIpPODavSvLYufjtIXevA4o17+Muns7j49cls35vLq3/owkc3nKgiLiISImKio3j2ko4kxEZz8euTWb4l2+uQQka5LuQFhY53fs7ivJcnkbFoMzf1bsZPd/amf4cUzHRBhYhIKKlbNZGv/3wyAHd8Npud+3I9jig0lNtC/uuyrZz5wkT+8d1CerWsxYS7+3Df2SdQUWOji4iErEY1KvLUxe1ZsGE3F7/+K7v25XkdkufKXSHfvOcAb8w5wB/emUJeQSGvDerC21enUiMp3uvQRESkFM5ql8KH13dj9fZ93PLJTPILCr0OyVPlrpC/NSGL6RsLuK1fC8b8tRfntFc3uohIuOnerAZPXNCeScu2cu+X8ygsLL8TrJS7fuS/9GtBc9vI5ae39DoUERE5Dpd2bcD6Xft54celxEYb/7ywfbkcMrvcFfIqibHUqVjuOiJERCLS7f1akFdQyKsZy4mJNh4f0K7c9bKWu0IuIiKRw8y464xW5Bc43pyYRWJsNA/0b+N1WGVKhVxERMKamXHf2SdwIK+At39eQcMaFbnqpEZeh1VmVMhFRCTsmRkPn9eWNTv2M2REJomx0QxMre91WGVCJ4tFRCQiREcZL13RmRMbV+eu/8zhvV9WeB1SmVAhFxGRiJEUH8NHN5zIGW2SefTbBfxn+hqvQwo6FXIREYkoMdFRvHRFZ05pXpN7vpzL+xHeMlchFxGRiJMQG80716Rxeutkhny7gOd/WIxzkTlojAq5iIhEpITYaF4b1IXL0hrw0rhlPPjN/IgcAU5XrYuISMSKiY7iyYvbU7ViLG9OyKKg0PHQuW0iaoKsyMlERESkGGbGfWedgGG8MWE5met38+ngk0iKkGKurnUREYl4BweNeeuqVBZs2M1Fr/3Cup37vQ4rIIJWyM1smJltNrP5RZZVN7OxZrbU/7NasPYvIiJyqDPa1uH967qyYdcBrn53Ciu37vU6pOMWzBb5+8BZhyy7D/jJOdcC+Mn/WkREpMz0bFGLd65OY8ueHAa8+guZ63d5HdJxCVohd85NBLYfsngA8IH/+QfABcHav4iIyOF0a1qDb/9yComx0Vzw6i988OtKr0M6ZmV9jjzZObcBwP+zdhnvX0REBIBGNSoy4tYe9GpRi0dGZPLmhOVhea+5BTNoM2sMjHTOtfO/3umcq1rk/R3OuWLPk5vZYGAwQHJycmp6enrA4srOziYpKSlg2/OScgk9kZIHKJdQFSm5hEoeeYWOt+bmMG1jAT3rxXBN2zhioo5uTvNA59K3b98Zzrm0Uq3snAvaA2gMzC/yejGQ4n+eAiwuzXZSU1NdIGVkZAR0e15SLqEnUvJwTrmEqkjJJZTyKCgodM//sNg1unekO+/ln938dTuP6vOBzgWY7kpZa8u6a30EcI3/+TXA8DLev4iIyP+IijL+dnpLXh/UhQ27DjDw9cl8P3+j12GVSjBvP/sUmAy0MrO1ZnYD8CRwupktBU73vxYREQkJZ7dP4bvbTqFVnUrc/PEMXs1YRm5+oddhlShow9o45644zFv9grVPERGR41W7UgLpg0/irv/M4Zkxi/ly5lreuiqN5rW9P59fHI3sJiIicoiE2GhevqIz716Txu79eVz46i+MnrfB67CKpUIuIiJSDDOjX+tkht96Ck1qVeTmj2dyx2ez2X0gz+vQfkeFXEREpAT1qiby5c0nc1u/Fgyfs56zhk7k12VbvQ7rv1TIRUREjiA2Ooo7Tm/JlzefTEJsNH94ZwqPfpvJgbwCr0NTIRcRESmtTg2q8t1tPbmmeyPe+2Ul5748iXlrvR2rXYVcRETkKCTGRfPogHZ8dMOJZB/I58LXfmH4slwKCr0Z3lWFXERE5Bj0bFGLMX/tRf8OKWTtKuQoR3UNGBVyERGRY1SlQiwvXt6ZWzrFY+ZNJVchFxEROU5x0R41x1EhFxERCWsq5CIiImFMhVxERCSMqZCLiIiEMRVyERGRMKZCLiIiEsZUyEVERMKYOefNkHJHw8y2AKsCuMmaQOhMXXN8lEvoiZQ8QLmEqkjJJVLygMDn0sg5V6s0K4ZFIQ80M5vunEvzOo5AUC6hJ1LyAOUSqiIll0jJA7zNRV3rIiIiYUyFXEREJIyV10L+ltcBBJByCT2Rkgcol1AVKblESh7gYS7l8hy5iIhIpCivLXIREZGIUC4LuZndZWbOzGr6X5uZvWRmy8xsrpl18TrGIzGzx/2xzjazH8ysrn95HzPb5V8+28we9jrWkpSQRzgek2fMbJE/3q/NrKp/eWMz21/kmLzhdaxHcrhc/O/d7z8ui83sTC/jPBIzu8TMMs2s0MzSiiwPx2NSbC7+98LmmBzKzIaY2boix+Icr2M6WmZ2lv93v8zM7ivzAJxz5eoBNADG4LsvvaZ/2TnAaMCAk4ApXsdZijwqF3l+G/CG/3kfYKTX8QUgj3A8JmcAMf7nTwFP+Z83BuZ7HV+AcmkDzAHigSbAciDa63hLyKM10AoYD6QVWR6Ox+RwuYTVMSkmryHAXV7HcRzxR/t/502BOP+xaFOWMZTHFvlQ4B6g6MUBA4APnc9vQFUzS/EkulJyzu0u8rIiv88nbJSQRzgekx+cc/n+l78B9b2M53iUkMsAIN05l+OcWwEsA070IsbScM4tdM4t9jqOQCghl7A6JhHoRGCZcy7LOZcLpOM7JmWmXBVyMzsfWOecm3PIW/WANUVer/UvC2lm9oSZrQEGAUW70Lub2RwzG21mbT0Kr9QOk0dYHpMirsfXo3BQEzObZWYTzKynV0Edo6K5hPtxKSqcj0lRkXBMbvWfxhlmZtW8DuYoef77jynLnZUFM/sRqFPMWw8Af8fXZfg/Hytmmect3JJycc4Nd849ADxgZvcDtwKPADPxDe2X7T/X9A3QosyCLsYx5hGWx8S/zgNAPvCx/70NQEPn3DYzSwW+MbO2h/RGlLljzCXkjktp8ihG2B6T4j5WzDLP/68UdYS/y68Dj+OL+XHgOXxfHsOF57//iCvkzrnTiltuZu3xnT+aY2bg6yqcaWYn4vsG1aDI6vWB9UEO9YgOl0sxPgG+Ax4p+ofIOTfKzF4zs5rOOc/GMz6WPAjTY2Jm1wDnAv2c/wSacy4HyPE/n2Fmy4GWwPQgh1uiY8mFEDwuR/Hvq+hnwvKYHEbIHZNDlTYvM3sbGBnkcALN899/uelad87Nc87Vds41ds41xvfL7+Kc2wiMAK72Xyl9ErDLObfBy3iPxMyKtrLPBxb5l9cx/zcV/5eUKGBb2UdYOofLg/A8JmcB9wLnO+f2FVley8yi/c+b4ushyfImytI5XC74jsvlZhZvZk3w5TLVixiPRzgekxKE9TE55NqXC4H5XsVyjKYBLcysiZnFAZfjOyZlJuJa5MdoFL6rpJcB+4DrvA2nVJ40s1ZAIb4r8G/yLx8I3Gxm+cB+4PIiralQdLg8wvGYvILvyuGx/u9SvznnbgJ6AY/5j0kBcJNzbrt3YZZKsbk45zLN7HNgAb4u91uccwUexlkiM7sQeBmoBXxnZrOdc2cShsfkcLmE2zEpxtNm1glfd/RK4EZvwzk6zrl8M7sV391Q0cAw51xmWcagkd1ERETCWLnpWhcREYlEKuQiIiJhTIVcREQkjKmQi4iIhDEVchERkTCmQi4iIhLGVMhFRETCmAq5yHEwswL/HMrzzexb+/85yJ2ZfVRkvRgz22JmI4sse8A/v/TB+di7HWYf2UGI+7i3aWZvmlmPQMRTwj5+Pc7PDzGzu/zPDx6rOWY208xODkyUIt5SIRc5Pvudc52cc+2A7cAt/uV7gXZmluh/fTqw7uCHzKw7vnHMuzjnOgCn8fsZlMJBN3xTnAaNcy6QxfbgseoI3A/8K4DbFvGMCrlI4Ezm99MXjgb6+59fAXxa5L0UYKt/8g6cc1udcyVOtGBmV5rZVH+r8s0iY4U/ZWZ/LrLeEDO783DrH2EfHc1sopktMLNCf8/Co8Ws1xpYcnAoUDPLMLPT/c//YWYvHWlfpXGw58DMGpvZQjN729+L8UORL0mHfuYBM1vsn3Gr1WE2XRnYEYgYRbymQi4SAP4i2Y/fT5aQjm8yiwSgAzClyHs/AA3MbIl/hrreR9h+a+AyoIdzrhO+8cEHFdnPZUVWvxTfLF6HW/9w+0gAPgPucs61AZ4AngWGFLP62cD3RV4/gm8q2kFAZ+BvJe3rGLUAXnXOtQV2AhcfuoL5piS93B/DRUDXIm8n+r/ULALewTdlpkjY06QpIscn0cxmA42BGcDYg2845+aaWWN8rfFRRT/kny8+FegJ9AU+M7P7nHPvH2Y//YBUYJp/EpNEYLN/W7PMrLaZ1cU3ocYOoP3h1i/BacBM59zBmbPmAmcdZtKdMykykY1zbqJ/1r07gD4lTdphZucDO5xzPx8hnkOtcM7N9j+fge93fqiewNcHZ2wzs6JfrPb7v9QcPLXxoZm1C/FJhUSOSIVc5Pjsd851MrMq+OZRvgUo2q08Al+rtg9Qo+gH/cVuPDDezOYB1wDvH2Y/BnzgnLv/MO9/gW/muzr4WuhHWr847YB5RV53AWaaWS3gaeAh4DF8re2qRU8FmFl7/v90wR7/smvxfUnZD2wAYv37GAUU+t/vjW/60ELn3BNHiC+nyPMCfF9OinPEwuycm2xmNfF98TnSFxyRkKaudZEAcM7tAm4D7jKz2CJvDQMec84VLZCYWSv7/VzsnfBN43o4PwEDzay2//PVzaxRkffT8XUpD8RX1I+0fnG24TsFgJm1xNc1ne6c2wKsBp7z53gKkFEklxTgY2AAsNfMziyyzTEHp3J1zj2Ib57sikXe/9459zi+Ah8IE4ELzSzRzCoB5xW3kpmdgG/KyW0B2q+IZ9QiFwkQfxf3HHwF9eCytcCLxayeBLzsv10tH9+864NL2PYCM3sQ+MHMooA8fK3/Vf73M/2Fa51zbgOwoaT1D+NT4Hwzmw9sBa5wzm0zsySgKZDvPyVwNr4vC5hZBeAr4E7n3EIzexx4Ct/czAC7/T+3+H/m4pvn/KC9/p9WQlyl5pybaWafAbPx5Vq0+/7gaZCD+7smzObtFimW5iMXkcMysxjgLeBRfBfRTQOeB7o55/KO8Nlr8XW1jzSzdOfc5f57umOA+UDNQ98PZi4ikUqFXEREJIzpHLmIiEgYUyEXEREJYyrkIiIiYUyFXEREJIypkIuIiIQxFXIREZEwpkIuIiISxlTIRUREwpgKuYiISBj7P0nLgjgv++uzAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 576x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "from scipy.special import erf\n",
    "\n",
    "w = 8  # wordlength of the quantizer\n",
    "A = np.logspace(-2, 0, num=500)  # RMS levels\n",
    "N = int(1e6)  # number of samples\n",
    "np.random.seed(1)\n",
    "\n",
    "def compute_SNR(a):\n",
    "    # compute input signal\n",
    "    x = np.random.normal(size=N, scale=a)\n",
    "    # quantize signal\n",
    "    xQ = uniform_midtread_quantizer(x, Q)\n",
    "    e = xQ - x\n",
    "    # compute SNR\n",
    "    SNR = 10*np.log10((np.var(x)/np.var(e)))\n",
    "    \n",
    "    return SNR\n",
    "\n",
    "\n",
    "def plot_SNR(A, SNR):\n",
    "    # plot results\n",
    "    plt.figure(figsize=(8,4))\n",
    "    plt.plot(20*np.log10(A), SNR)\n",
    "    plt.xlabel(r'RMS level $\\sigma_x / x_\\mathrm{min}$ in dB')\n",
    "    plt.ylabel('SNR in dB')\n",
    "    plt.grid()\n",
    "\n",
    "\n",
    "# quantization step\n",
    "Q = 1/(2**(w-1))\n",
    "# compute SNR for given RMS levels\n",
    "SNR = [compute_SNR(a) for a in A]\n",
    "# plot results\n",
    "plot_SNR(A, SNR)\n",
    "# find maximum SNR\n",
    "Amax = A[np.argmax(SNR)]\n",
    "Pc = 1 + erf(-1/(np.sqrt(2)*Amax))\n",
    "print(r'Maximum SNR = {0:2.3f} dB for A = {1:2.1f} dB with clipping probability {2:2.1e}'\n",
    "      .format(np.array(SNR).max(), 20*np.log10(Amax), Pc))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise**\n",
    "\n",
    "* Can you explain the overall shape of the SNR?\n",
    "* For which RMS level and probability of clipping is the SNR optimal?\n",
    "* Change the wordlength `w` of the quantizer. How does the SNR change?\n",
    "\n",
    "Solution: The SNR is low for low RMS levels of the input signal since the relative level of the quantization error is high. The SNR increases with increasing level until the clipping errors become dominant which make the SNR decay after its maximum. The SNR is optimal for $A \\approx -12$ dB which is equivalent to $\\Pr \\{ |x[k]| > x_\\text{max} \\} \\approx 10^{-4}$. Increasing the wordlength by one bit increases the SNR approximately by 6 dB."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Laplace Distributed Signal\n",
    "\n",
    "The [Laplace distribution](../random_signals/important_distributions.ipynb#Laplace-Distribution) is a commonly applied model for speech and music signals. As for the normal distribution, clipping will occur with a non-zero probability. The probability that clipping occurs can be derived from the [cumulative distribution function](../random_signals/important_distributions.ipynb#Laplace-Distribution) (CDF) of the normal distribution as\n",
    "\n",
    "\\begin{equation}\n",
    "\\Pr \\{ |x[k]| > x_\\text{max} \\} = e^{- \\frac{\\sqrt{2}}{A}}\n",
    "\\end{equation}\n",
    "\n",
    "The SNR for a Laplace distributed signal is in general lower compared to a normal distributed signal. The reason for this is, that the Laplace distribution features low signal values with a higher and large values with a lower probability in comparison to the normal distribution. The relative quantization error for small signals is higher, which results in a lower average SNR. The probability of overload distortions is also higher compared to the normal distribution."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Example - Quantization of a Laplace distributed signal\n",
    "\n",
    "The following example evaluates the SNR of a linear uniform quantizer with $w=8$ for a Laplace distributed signal $x[k]$. The SNR is computed and plotted for various RMS levels."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Maximum SNR = 35.581 dB for A = -16.6 dB with clipping probability 7.1e-05\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfIAAAEOCAYAAABywZWyAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xd8FNX6x/HPk0YCoRNCJ/ReQxURsGJDxd47eu1dr167WK967Yq9EitiQRQRpCgovYXeey8hpJ/fH7to5BdCgN3M7ub7fr32ld3Z2ZnnyUCePWdmzjHnHCIiIhKeorwOQERERA6dCrmIiEgYUyEXEREJYyrkIiIiYUyFXEREJIypkIuIiIQxFXIREZEwpkIuIiISxlTIRUREwpgKuYiISBiL8TqAkqhRo4ZLSUkJ2PZ2795NhQoVArY9LymX0BMpeYByCVWRkkuk5AGBz2Xq1KmbnXNJJVk3LAp5SkoKU6ZMCdj2xo4dS9++fQO2PS8pl9ATKXmAcglVkZJLpOQBgc/FzFaUdF11rYuIiIQxFXIREZEwpkIuIiISxlTIRUREwpgKuYiISBhTIRcREQljKuQiErJ2Z+exPTOHhRt2kZWb73U4IiEpaPeRm1k8MA4o59/PF865B83sPaAPsMO/6mXOuRnBikNEwk9+gWN3Th5XvTeFP5ZvBeC6vk24q39LjyMTCT3BHBAmGzjaOZdhZrHABDP7wf/enc65L4K4bxEJU1syshnw8kTWbN/zj+XfzFzLvHU7+d+5HalSPs6j6ERCT9AKuXPOARn+l7H+hwvW/kQk/HwxdTWLNu6iZ+PqDBm3lPwCx+RlvhZ469qV2JSRzbld6vPymMWs3raH1dv20PGRUQzsXJfnzunocfQioSGoQ7SaWTQwFWgKvOKcm2xm/wIGm9kDwGjgHudcdjDjEJHQsmprJmMWbOSB4XMBeOPXpQDERBngK+Ijbu5NfoEjOspYvDGDkXPX//X5r6at4azUenRvVJ3svHwmLt5Cw+rlaZ5csfSTEfGY+RrOQd6JWRVgGHAjsAVYD8QBQ4AlzrlHivjMIGAQQHJycmpaWlrA4snIyCAxMTFg2/OScgk9kZIHBD4X5xyzN+fz3FTfd/eW1aJoXjWab5bk0rtuDGc3j2NNRgG1KhhV4/++FndoejY/rsijQ1I013csx81jMulZO4Zq8cbPK/PYnu37O/Ze//1PWqHjEnoiJQ8IfC79+vWb6pzrUpJ1S2XSFOfcdjMbC/R3zv3XvzjbzN4F7tjPZ4bgK/R06dLFBXIweg3UH5oiJZdIyQMCl8umXdms2pbJA8PnMGdN5l/LbzmpEz0aVyPri1ncc2JLGicV/YdwcfRSflyRTtP6tTn+mA40S59Abnwsawsc27O3/LVem9SeJFUsF9RcQkGk5BIpeYC3uQTzqvUkINdfxBOAY4GnzKy2c26dmRlwOjAnWDGIiPemrtjKeUMmkZv/d+/fye1qUz0xjmNa1SQ2Ooohl5So4UF8rK+VnlK9AtNWbmNHZu4/3v9l/gbO7dogcMGLhIFgtshrA+/7z5NHAZ85574zs1/8Rd6AGcC1QYxBRDyUk1fAW+OXUT4uhtuPb07F+Bh6Nq5BzYrliPKfDy+JU9rX4aNJK7jyyEYApFQvzzcz1/71foW4aHLzHUs27WbUvA2s3pbJ5b0aBTwfkVAUzKvWZwGdilh+dLD2KSLecs7x5bQ1fD9rLU1rJpL2xyp2ZedxSc+GXNIz5ZC3W6tyPGPv7PfX64bV/z4X3q9FErUqJzBl+VaWbd7NkHG+C+dUyKWsKJVz5CIS+XbsyeWmodP5deEmalYsx9iFm2hfrwpdGlb9qyUdKA2rlwegckIs717eDYCrP5jCii27/1onKzef+NjogO5XJBSpkIvIYSkocHw8eQUvj1nM1t05PDygDRf3aIj5e87NSt6FXlLt6lXmsiNS/vEFIaV6eUbN2/DX69RHR/HIaW05M7VewPcvEkpUyEXkkK3amsnD387j5/QNdG9UjVcu6EyXlGpB32+5mGgeGtDmH8tSavzz1rPdOfnc/vlMFXKJeCrkInLQcvMLeHXMEl76ZREOeOCU1lzeKyUore+SaluncpHL09ftpFXtSqUcjUjp0exnInJQCgocNw2dzvM/L+Tk9rUZf1c/rjiykadFHKBD/Spc1KMB5eP+eV78l/kbPYpIpHSokItIiQ2bvprG947ghznrue245rxwXifqVEnwOqy/PHZ6O2Y9ePxfr1vVrsSERZs9jEgk+NS1LiIHNHftDm77dCYLNuwCfAO6XNunicdRFS0mOoq7+regbpUE5q7dyXsTl7MnR3OZS+RSIReR/covcDwxIp0h45eSlFiOO45vzkU9Gob8NKLX9W0KQPm4GIaMW8q8dTs9jkgkeFTIRaRIizbsYvDkLJbuWMr53epz1wktqVohtAv4vtrW9V3kNnftDhrgG/O9RmKc5+fzRQJJhVxE/mHTrmzuHTabn9M3UCEGXjivI6d1rOt1WIekVqV4qlWIY+6anRTE5nPlE6MZcnEqx7RK9jo0kYBRIReRv4xbuInbPpvJrqxcbuzXlKZuLQPCtIiDbzCaNnUqMXvNDlaQS36BY9nm3Qf+oEgYUSEXEbbuzuH+4XP4ftY6micn8vFV3WlRqyJjx67zOrTD1i2lGs/9vJC9nembM3I8jUck0FTIRcq4z6es4qmR89mZlcdNRzflX32bkhAXOWOUn9WlHv8bvQicI7FcDJszsr0OSSSgVMhFyqitu3N4+Nu5DJ+xlq4pVXnw1Da0rVv06GjhrHblBC4/IoX1a1ezPDuBLSrkEmFUyEXKoO9nreOB4XPYsSeXW49tzg1HNyX6IOYHDzf/OaU1Y8du5L1l5dS1LhFHhVykDNm4K4sHvp7LyLnraVe3Mh9d1b1MjUNeI7EcC9bv4tuZa3HAgA51vA5J5LCpkIuUAc45hk1fw8PfzmNPbj5392/J1b0bERNdtkZprp4Yx7odWdw4dDqgQi6RQYVcJMKt3b6H+4bNZsyCTaQ2rMrTZ7WnSVKi12F5onqhAW0qxevPn0QG/UsWiVDOOYb+sYrHR6STX+B48NTWXNIzJaLPhR9IQpzvT16V8rHs3JNLQYHDAet27KFe1fLeBidyiFTIRSLQyi2Z3PPVLH5bsoUjmlTnyYHtaVBdhers1HokJcaxZnsWj343j+17crnj85n8Mn8jsx46nkrxsV6HKHLQVMhFIkh+geP935bzzI8LiI4ynhjYjvO61tfY4n7xsdH0b1ub4TPWADBm/sa/5itPX7uT7o2rexmeyCFRIReJEIs3ZnD3l7OYumIb/VokMfiMdiE1V3goqV6hHADP/Ljgr2Xp61TIJTypkIuEubz8At4cv4znf15IQmw0z53TgTM61VUrvBjV/Be9rd+ZxRW9GjF8xhpNdSphS4VcJIzNX7+TOz+fxew1O+jfphaPnN6GmhXjvQ4r5FVP/Pvq9f5ta7Fwwy7S1+3yMCKRQ6dCLhKG8gscQ8Yt5blRC6icEMurF3bmpHa1vQ4rbFQt/3ch71C/Mq1qV+T931eQl19Q5u6tl/AXtEJuZvHAOKCcfz9fOOceNLNGQBpQDZgGXOyc05iJIiW0cksmt38+gz+Xb+PEtrUYfEa7v7qKpWTiYv4u1uViomlVuxI5eQUs3byb5skVPYxM5OAFs0WeDRztnMsws1hggpn9ANwGPO+cSzOz14ErgdeCGIdIRMgvcPyxbCuDPpwCDp0LP0z3n9Ka5sm+gXFa1/ENU5u+bqcKuYSdoBVy55wDMvwvY/0PBxwNXOBf/j7wECrkIsVKX7eTWz+dwfz1u6hbJYG0QT2oX033hR+OK49s9NfzJkmJxEYbz/60kMY1EmlXrzLZeflEmRGrrnYJcUH9F2pm0WY2A9gIjAKWANudc3n+VVYDdYMZg0g4yy9wvDZ2CQNensDmjBweO70tw2/opSIeYLHRUXRrVI2VWzM59eUJ/LZkM72eHMPJL473OjSRAzJfwznIOzGrAgwDHgDedc419S+vD4xwzrUr4jODgEEAycnJqWlpaQGLJyMjg8TEyBhrWrmEnkDlsX53AW/PzmbR9gK6JEdzaZtyVIwr3W70SDkmcOBc8gsc63c77pu4h4HNYvlqUS4A7/WvUFohllikHJdIyQMCn0u/fv2mOue6lGTdUrlq3Tm33czGAj2AKmYW42+V1wPW7uczQ4AhAF26dHF9+/YNWDxjx44lkNvzknIJPYebR05eAf/7eSFDJi4lIS6a589tx+kdvTkXHinHBEqey0uzRzNjWzTgK+R9+vQJuesQIuW4REoe4G0uwbxqPQnI9RfxBOBY4ClgDHAWvivXLwWGBysGkXAza/V27vpiFvPX7+Ls1HrccUILkivpvvDS1Cw5kfGLNv/1emdWHpUTNAa7hK5gtshrA++bWTS+c/GfOee+M7N5QJqZPQZMB94OYgwiYSEvv4APJ61g8PfpVE+M461LunBs62SvwyqT6u4zrO36HVkq5BLSgnnV+iygUxHLlwLdgrVfkXAzf/1O7hs2h6krttGneRIvnteJyuVVOLzSo3F10v5cxT0ntuTJH+azbsceWtTSLWkSujSym4hHdmXlMvj7dNL+XEVCbDT/O7cjp3WsE3LnY8ua0zrWoU/zJDJz83nyh/ms35HldUgixVIhF/HAhEWbueuLmazfmcU1fRpzzVFNNDpbiDAzqlaIIzG/ADNYq0IuIU6FXKQU7cnJ56mR83nvt+U0SarAl/86gk4NqnodlhQhNjqKulUSWLZ5t9ehiBRLhVyklMxbu5Ob0qazeGMGl/dK4e7+LYmPjfY6LClGy1oVma/pTSXEqZCLBFlufgHvTFjGsz8tpEr5WD68shu9myV5HZaUQItaFRmzYBPZefmUi9GXLglNGkRYJIgmL93CyS+O54kf5tO3RRIjbzlKRTyMtKxVifwCR4v/jOSPZVu9DkekSCrkIkGwaVc2t306g3OHTGJ3dj5vXtKFIZd00QVtYaaNf1Y0gCvf/5MdmbkeRiNSNHWtiwTYwm353PnieLZn5nBDv6Zc368pCXHqlg1HjZMS+eyansREGwNf/Y1P/ljJv/o28ToskX9QIRcJkIzsPB76Zi5fTM2ibpUEvr3xSFrWqnTgD0pI69aoGgBHNKnOx5NXcG2fxrrXX0KKutZFAuDbmWs58YVxfDVtNac0juWnW49SEY8wp3Wsw+pte5i7VlexS2hRIRc5DDsyc7n10xncOHQ6ieViGXp1D85qHkeFcursijTHta5FlMEPc9Z5HYrIP+ivjcgh+mX+Bv791Ww2Z+Rw8zHNuOmYZkRHGWNXeh2ZBEO1CnEc1TyJV8YsYcLiLTxzVnuaJ2sMdvGeWuQiBykjO4+7v5jFFe9NoUpCHF9f14tbj2tOdJTOm0a6i3s0BGDmqu3cN2w2zjmPIxJRi1zkoExcvJl7vprFmm17uK5vE24+tpkGCilD+raoyZ0ntGDllkw+nbKKlVszaVi9gtdhSRmnQi5SAqu2ZvLY9/P4ce4GGlYvz2fX9KRLSjWvw5JSFh1lXN+vKXPX7uDTKauYuXqHCrl4ToVcpBjOOdL+XMXD384lyow7T2jBlUc20hjpZVzz5IrEx0Yxc9V2BnSo43U4UsapkIvsx+KNGTz87VzGL9rMkU1r8MzZ7aldOcHrsCQExEZH0aZOZb6ctprLjkihfrXyXockZZgudhMpwpdTV3PqSxOYvWYHD5zSmg+u6KYiLv9wUY8GZObkc/eXs7wORco4tchFCtmwM4snRqTz9Yy1dGtUjRfP60StyvFehyUh6IxO9Zi/bhfvTFxGZk4e5eP051S8oRa5CL5z4e9NXEa//47l+9nruOmYZnxyVXcVcSnWkc1qkJvvmKyZ0cRD+gopZd7mjGz+M2wOI+eup2+LJB4e0EZXIkuJdE2pRoW4aN7/bTl9mydpDHbxhFrkUmY55/hsyiqOfe5XRs/fwH0nteLdy7qqiEuJxcdGc/vxLRi7YBOPj0jXADHiCbXIpUxavS2TOz+fxe9Lt9A1pSqPn9GOZhpuUw7BpUeksGzzbt4cv4wq5eO4vl9Tr0OSMkaFXMqUvPwChoxfypvjlpKX73j8jHac17U+URpeVQ5RdJTxyGlt2JqZw/OjFnJ6p7rUraI7HKT0qGtdyozV2zK5+O0/eHrkAtrWrczXN/Tigu4NVMTlsJkZ953UCjN449clXocjZUzQCrmZ1TezMWaWbmZzzexm//KHzGyNmc3wP04KVgwi4DsX/tGkFZzw/Dhmrd7O02e158Mru9MkKdHr0CSC1KmSwIAOdfly6mp2Z+d5HY6UIcHsWs8DbnfOTTOzisBUMxvlf+9559x/g7hvEcA3RvrdX87ityVb6NW0Ok8ObK9RuCRozu9Wny+nreb72es4p0t9r8ORMiJohdw5tw5Y53++y8zSgbrB2p9IYQUFjo8nr+CJH+YTZcbjZ7Tj/G71dXuQBFVqw6o0TqrAp3+uUiGXUmOlcbuEmaUA44C2wG3AZcBOYAq+Vvu2Ij4zCBgEkJycnJqWlhaweDIyMkhMjIxuVeXy/23KLOCdOdmkby2gbfVoLm8bR/WE0rscRMckNJVWLj8sy+XTBTkM7pVA3YrB+XcXKcclUvKAwOfSr1+/qc65LiVa2TkX1AeQCEwFBvpfJwPR+M7PDwbeOdA2UlNTXSCNGTMmoNvzknL5W3Zuvvto0nLX+v4fXJsHRrqhk1e4goKCwAR3EHRMQlNp5bJpV5Zrff8P7tw3fnP5+cH59xcpxyVS8nAu8LkAU1wJ62xQmylmFgt8CXzsnPvK/8Vhg3Mu3zlXALwJdAtmDFI2zFu7kxNfGMd9w+bQrl5lRt7Sm/O6NVBXupS6GonlePDUNkxaupVnflrgdThSBgTtHLn5/oK+DaQ7554rtLy2850/BzgDmBOsGCTyOef4aPJKHv1uHlUSYnnzki4c26qmCrh46uwu9ZixejuvjV1Cs5qJDOxcz+uQJIIF86r1XsDFwGwzm+Ffdi9wvpl1BBywHLgmiDFIBJu3dif3DpvNjFXb6dM8iefO6UD1xHJehyWCmfHwgDYs3pDBI9/N46R2tYmPjfY6LIlQxRZyM+sDbHPOzTKzc4CjgCXAq8657OI+65ybABTVLBpxqMGKAOTkFfDSL4t4bewSqpSP5emz2nNW53oa2EVCSmx0FLcd35zzhkxi2PQ1nN+tgdchSYTabyE3s1eA9kA5M1uI76K1kcARwDvAhaUSoUghc9bs4I7PZzJ//S4Gdq7L/Se3pmqFOK/DEilS90bV6FCvMo+PSKdn4+qk1NCEPBJ4xV3s1s851xtfK/xE4Ezn3OvApfgKvEipyckr8I1j/cpEtuzO4e1Lu/DcOR1VxCWkmRkvX9CZ3PwC3p6wzOtwJEIVV8izAJxzWcAK51y+/7UDckshNhEApq/cxqkvTeCF0Ys4pX1tRt16FMe0SvY6LJESqV+tPP1a1GTk3PXkF2iaUwm84s6R1zSz2/Cd5977HP/rpKBHJmVedl4+b09Yxn9/XEDNivG8fWkXFXAJSye2q80Pc9bzc/oGTmhTy+twJMIUV8jfBCoW8RzgraBFJALMXr2DQR9OYd2OLE5uX5snB7ajYnys12GJHJL+bWrRPDmRB4bPoU/zJF3BLgG130LunHu4NAMRAcjKzefJH+bzwe/LqVUpnncv70rf5km6L1zCWlxMFPec2JIr3pvCb0s2c3RL9SxJ4BR31fqLxX3QOXdT4MORsmzO5nzueWYs63dmcUnPhtx6bHNdzCYRo1fTGpjBFe9N4fEz2nFBd92OJoFR3MVuU/2PeKAzsMj/6AjkBz80KSs27szi0e/m8eyULCrGx/DpoB48clpbFXGJKOViohnQoQ4A9w6brTnLJWCK61p/H8DMLsN3K1qu//XrwE+lEp1EvHELN3HrpzPYlplD73oxvD6oF+XjgjngoIh3Hju9LUe3rMnNaTP4ZuZaDRIjAVGSSVPq8M8L3RL9y0QO2bbdOdz+2UwueecPqpSP5cdbjuKKtuVUxCWiVYyPZUCHOjSoVp7R6Ru8DkciREn+aj4JTDezMf7XfYCHghaRRLzZq3dw1Qd/siUjhxv6NeWGo5sSHxvNmnSvIxMJPjPjyGY1+GbGWnLzC4iNDuoklFIGHPBfkHPuXaA7MMz/6Lm3213kYGzPzOGOz2dy6ssTAPj6+l7ccUIL3YojZU7vpjXIyM6j2X0/sHb7Hq/DkTBXoq+Czrn1zrnh/sf6YAclkcU5x3ez1tL7qTF8NW011/Zpwsibj6Jt3cpehybiiT4tkji+te8WtLfGa+hWOTzq05Gg2rY7h0EfTuWGT6bTNDmRETf35p4TW+qKdCnTysfFMOSSLgzsVJe0P1eyPTPH65AkjKmQS9B8NW01Rz09hrELNnLvSS35dFBPWtaq5HVYIiHjmj5NyMzJ592Jy70ORcJYiS4RNrNoILnw+s65lcEKSsLb1t05DP4+nS+nraZ7o2o8fFobFXCRIrSoVZH+bWrxwuhFfDx5JS+e35EjmtTwOiwJMwcs5GZ2I/AgsAEo8C92aCpT2Ydzji+mrubxEensysrjur5NuOXY5sTFqONHZH+eO7cDNUbE8dGklbz/23IVcjloJWmR3wy0cM5tCXYwEr6WbMrgvmGzmbR0K10aVuXxge1onlzxwB8UKePKx8Xw2OntiIuO5qNJK9iRmUvl8pogSEquJE2lVcCOYAci4WlvK/zEF8Yzb+1OnhjYjs+u6akiLnKQTu9Uh5z8AkbMWed1KBJmStIiXwqMNbPvgey9C51zzwUtKgkLc9fu4IHhc5m6YhvdG1XjpQs6UbNivNdhiYSldnUr07hGBb6evkZDt8pBKUmLfCUwCojDN1Tr3oeUUc450v5YyRmv/sbKrZk8ObAdH1/VXUVc5DCYGWem1mPysq0M/WMlzjmvQ5IwccAWueYll8IWb9zFnV/MYvrK7fRuVoP/nduR6onlvA5LJCJc3bsxExdv5t9fzWb1tkzuPKGl1yFJGChuPvL/OeduMbNv8V2l/g/OuQFBjUxCSkGB473flvPUyPlUKBfDU2e246zU+kRHmdehiUSMuJgo3r+iGzcNnc67E5dzTZ8mVIrXhW9SvOJa5B/6f/63NAKR0DVv7U4Gj5jHxMVbOKZlTZ48sz1JFdUKFwmG2Ogo/tW3CT/MWc9nf67iqt6NvQ5JQlxx85FP9f/89VA2bGb1gQ+AWvjuPx/inHvBzKoBnwIpwHLgHOfctkPZhwRXVm4+z/60gLcmLCMxLoYnB7bj3K71MVMrXCSY2terQs/G1Xnh50V0TalGh/pVvA5JQlgwR+rIA253zrUCegDXm1lr4B5gtHOuGTDa/1pCzPSV2zj5xfG8OX4ZF3RrwIS7j+a8bg1UxEVKye3HN2dXdh4DX/uNJZsyvA5HQljQCrlzbp1zbpr/+S4gHagLnAbsnQb1feD0YMUgBy87L5+nRs7nzNd+Y09OPh9e2Y3BZ7TTABUipaxLSjV+uLk3+QW+sRpE9ueQCrmZNTzI9VOATsBkINk5tw58xR6oeSgxSOD9Mn8DJ70wntfGLuGs1HqMvPUoejdL8joskTKrVe1KHN2yJp9PWcXu7Dyvw5EQZcXdq2hmPfG1osc55zaaWXt8XeG9nXP1S7QDs0TgV2Cwc+4rM9vunKtS6P1tzrmqRXxuEDAIIDk5OTUtLe1g8ipWRkYGiYmJAduelwKRS26B4/MFOfy0Io/aFYzzWsbRIalE8+kEVKQcl0jJA5RLKFi8LZ/HJmdxdIMYLmoVR5RZ2Oayr0jJAwKfS79+/aY657qUaGXnXJEP4Bl83eFDgT/5e+KUm4H4/X1un23EAj8CtxVatgCo7X9eG1hwoO2kpqa6QBozZkxAt+elw83l53nrXe+nfnEN7/7OPfTNHJeVmxeYwA5BpByXSMnDOeUSKu7/erZrePd37umR6c658M6lsEjJw7nA5wJMcSWos865Ym8/Oxno5JzLMrOqwFqgvXNuUUm+IJjvqqi3gXT3z+FcvwEuBZ70/xxeku1JYO3MyuV/oxbxzsRlNK2ZyIdXdlM3ukiIenhAG7Zk5PDBbyu4pk8Tr8OREFNcId/jnMsCcM5tM7MFJS3ifr2Ai4HZZjbDv+xefAX8MzO7Et/wr2cfQtxyGOas2cH1n0xjxZZMLurRgAdOaaOpRkVCmJnxr75NGDFnHZe98wfXNNfwrfK34gp5EzP7ptDrlMKv3QFGdnPOTQD2d6/SMSUPUQJlT04+z41awNsTllGtQjm+uLYnXVKqeR2WiJRA27qVefn8zlz/yTSGRcVwwrFeRyShorhCfto+r58NZiASXJOXbuHuL2exfEsmF3RvwN39W1I5QbeUiYSTk9vXZsLiBnz650rmrNlB27qVvQ5JQkBxI7sd0ohuElpWb8vkmR8XMHzGWhpUK88nV3fniCY1vA5LRA7RPf1b8v2MVfz7q9l8fX0vzXcgxU6aMoYiJkvxc845dY+HMOccn/yxkse/T8cB1/Zpwk3HNKV8XOnfViYigVO5fCwXtYrj1Zk7eGPcEq7r29TrkMRjxf1Vv6OIZT2Au4CNwQlHAmH1tkzu+XI2ExZvplfT6jw5sD31q5X3OiwRCZCutaI52dXm6ZELqFe1PAM61PE6JPHQASdNATCzPsD9QDngWufcD6UQmxwk5xxD/1jF4yPScc7x2OltubC7xkcXiTRmxrNnd2DTzmxu/2wGLWtVpHlyRa/DEo8Ue8+RmZ1gZhPwFfHBzrneKuKhadXWTC5++w/uHTab9vUqM/KWo7ioR0MVcZEIFR8bzesXpxITFcXrvy7xOhzxUHHnyP8EkvCN8Pa7f1nnve87/4Qo4q0C53hnwjKe+XEBUQaPnt6Wi9QKFykTqlWI47xu9fng9xVc0K2Bbicto4o7R74byADOAs7kn/eEO+DoIMYlJbBqayaPT85i8fZ59GuRxGNntKNulQSvwxKRUnTLsc0ZM38jNw6dzs+39aFCOV3QWtYUd468bynGIQchN7+At8Yv4/Vfl5CbW8AL53VkQIc6aoWLlEGVE2JFoTBMAAAZzElEQVR59pyOnPnab5w75HeeOrM9bero/vKyZL/nyM2sq5nVKvT6EjMbbmYvmpn6bzyycMMuznh1Ik+NnE/nBlV4oGcCp3WsqyIuUoalNqzKY6e3ZeH6DP7z9Zy9E1RJGVHcxW5vADkAZnYUvjHSPwB2AEOCH5oUticnn1fHLuaUFyewbnsWr1/UmXcv70atChojXUTgoh4Nuf/U1kxfuZ3R6bpDuCwp7mRKtHNuq//5ucAQ59yXwJeFJkGRUjBx8WZu/XQGG3dlc2LbWjx2eluqJ5bzOiwRCTHndqnPR7+v4P7hc+jeuBoV4zUMc1lQXHMu2sz2FvpjgF8KvaerKUrB7uw87v96Dhe+NZmK8TF8OqgHr17YWUVcRIoUFxPFU2e1Z8POLG7/bCZrtu/xOiQpBcUV5KHAr2a2GdgDjAcws6b4utcliGav3sF1n0xl9bY9XHlkI+44vgUJcdFehyUiIa5j/SrcdEwz/vfzIlZv28OIm3t7HZIEWXFXrQ82s9FAbeAn9/fVE1HAjaURXFmUm1/AkHFLeWH0IpISy/HZNT3pqntDReQg3HJscyrExTB4RDqzVm+nfb0qXockQVTslVLOuUnOuWHOud2Fli3UYDDBMWv1dk59aQLP/LiAY1vVZPgNvVTEReSQDOxcl5goY8DLExk1b4PX4UgQ6ZLnELAnJ5/HR6Rz+isT2bo7hzcuTuXVC1OpoXPhInKIqieW4/0rulG/WgL3DpvNwg27vA5JgkSF3GO/Ld5M/xfGMWTcUs7tWp9Rt/XhhDa1DvxBEZED6NW0Bm9c1IX8AseZr/3Gjsxcr0OSIFAh98iOPbnc/cUsLnhrMgZ8cnV3nhjYnsoJul1ERAKndZ1KfHhlN3Zl5fHR5BVehyNBoNvIPDByzjruHz6XrbtzuKZPY249tjnxsboiXUSCo02dyvRrkcQLPy+iRXJFjm2d7HVIEkBqkZeijbuy+NdHU7n2o2kkJZZj+PW9+PeJrVTERSTonj+3I61qV+Taj6by7cy1XocjAaQWeSlwzvH5lNU89v08svIKuKt/C67u3ZjYaH2PEpHSUaV8HB9d1Z0r35vCzWnTqVU5XnfFRAhVkiCbv34nF741mbu+nEXL2pUYeXNvruvbVEVcREpdxfhY3r28K7UrJ3D/13PIzsv3OiQJAFWTIMnJK+CRb+dx0gvjmbt2J4PPaEva1T1onJTodWgiUoZVKBfDQwPaMH/9Lm4aOp31O7K8DkkOk7rWg2Di4s08/eMCZq7azkU9GnDH8S2oUj7O67BERAA4rnUyd/Vvwf9GLWLDzqkMu+4ITYUcxoJWyM3sHeAUYKNzrq1/2UPA1cAm/2r3OudGBCuG0paRncfg7+cx9I9VVC0fy2sXdubEdrW9DktE5P+5rm9TKsXH8p+v5/D7ki0c0bSG1yHJIQpm1/p7QP8ilj/vnOvof0RMEf99yRaOf+5X0v5cxbV9mjDp3mNUxEUkpJ2VWo+6VRK4+6tZbM/M8TocOURBK+TOuXHA1gOuGOa2ZGTz6HfzuPjtycTHRvPFtT2558SWlIvRLWUiEtriY6N58fxObNiRzYVvTWb1tkyvQ5JD4MXFbjeY2Swze8fMqnqw/4D5ce56+v53LO9OXMbAznUZfkMvUhvqdg4RCR+pDavyxiWprNiSyYkvjGd0uiZYCTf29+ykQdi4WQrwXaFz5MnAZsABjwK1nXNX7Oezg4BBAMnJyalpaWkBiysjI4PExEO/enxnjuOrhTmMXZ1Ho8pRXN2uHHUSvbkB4HBzCSWRkkuk5AHKJVQFI5dNmQW8OD2bHdkF/LdPeeKig3/xm47J/vXr12+qc65LiVZ2zgXtAaQAcw72vX0fqampLpDGjBlzSJ8rKChwX05d5do+ONI1/vf37rHv5rqs3LyAxnawDjWXUBQpuURKHs4pl1AVrFwmLtrkGt79nWt+3wg3c9W2oOyjMB2T/QOmuBLW2lJtRppZ4au/zgDmlOb+D8f2zBxuGDqd2z6bSctaFfnxlt7cd3JrnQsXkYjRs0l1LujegOy8Al7+ZbHX4UgJBfP2s6FAX6CGma0GHgT6mllHfF3ry4FrgrX/QMnLL2DI+KW88stisvMKuPOEFlzbpwnRUbrnUkQii5nx+BntfLfPjl3Ch5NWcGG3BkTp711IC1ohd86dX8Tit4O1v2BYtnk3t302g+krt3NCm2RuPa45LWtV8josEZGgurp3Y6av3M79X8/h9yWbeeWCzhowJoRpiNYiOOf48PflnPTCeJZu2s2L53fijYu7qIiLSJlQpXwcH1/VnZuPacaI2eu57N0/WblFt6aFKhXyfSzemMEl7/zB/cPn0rVRNX685SgGdKjjdVgiIqXKzLjh6KZ0rF+FXxdu4oXRi7wOSfZDY6377cjMZeifK3l65HziY6MZfEZbLujWQN1JIlJmxUZH8fX1vbhv2Gw+nryS6olx3HF8C+Ji1AYMJSrkwOj0Ddzx+Uy2ZeZybKtknjqzHdUTy3kdlohISLi8VyNGp29kyLilbN6VzXPndvQ6JCmkTBfyXVm5PDdqIe9OXE7r2pV485IupDasqla4iEghTWsmMuneY3h+1EJeGL2IPi2SOK1jXa/DEr8yWcidc3w5dTWPfDePHXtyueyIFO45sSXxsbonXERkf248uinjF23iri9mUSkhln4tanodklAGL3bblZXLG7Oyuf3zmTRPTuTbG47koQFtVMRFRA4gJjqKty7tSuOkRK54708+nLTC65CEMljIn/1pIZPX5XPH8c1JG9STdvUqex2SiEjYqFYhjq/+dQRHt6jJ/V/PYeSc9V6HVOaVuUJ+67HNubd7PDcc3Uyjs4mIHIKEuGhevagzHepX4aa06Xzw+/K9c2iIB8pcIa9cPpZmVdWNLiJyOMrFRPPOpV04okl1Hhg+l6d/XOB1SGVWmSvkIiISGNUTy/HOpV05O7Ueb/y6hN+XbPE6pDJJhVxERA5ZVJTxn1Na06hGBS579w/em7jM65DKHBVyERE5LJUTYvn82iPo1bQGD307jzfHLdU581KkQi4iIoetWoU43rykCye0SWbwiHQeGD7X65DKDBVyEREJiOgo47ULU7miVyM+nLSC//64gPwCtcyDrUyO7CYiIsERFWXce1JLdmbl8vKYxcREG7cc29zrsCKaCrmIiARUTHQUz5zVnoICx/9+XkRsdBTXHNWYmGh1AgeDfqsiIhJwZsbjA9txcvvaPPPjAk55aQJz1+7wOqyIpEIuIiJBER8bzcvnd+L1izqzLTOHC9+azPz1O70OK+KokIuISNCYGf3b1ubza46gXEwU/f83ngEvTyArN9/r0CKGCrmIiARdg+rlSRvUkxPaJDNr9Q4e/nYeebqiPSBUyEVEpFQ0qlGBNy7uwpVHNmLoHyt5bWY2BSrmh02FXEREStX9p7TmPye3YuqGfE5/dSIbdmZ5HVJYUyEXEZFSd+WRjbiqXRyLNmQw6MOpLNmU4XVIYStohdzM3jGzjWY2p9CyamY2yswW+X9WDdb+RUQkdJkZR9aN5blzOrB4wy5OeXECI+es8zqssBTMFvl7QP99lt0DjHbONQNG+1+LiEgZdWK72oy5oy8talXk2o+m8dTI+WTm5HkdVlgJWiF3zo0Dtu6z+DTgff/z94HTg7V/EREJDzUrxZM2qAfnd6vPa2OX0PupMUxYtNnrsMJGaZ8jT3bOrQPw/6xZyvsXEZEQFB8bzRMD2/Plv46gemIcV7z3J5OWbvE6rLBgwZwz1sxSgO+cc239r7c756oUen+bc67I8+RmNggYBJCcnJyalpYWsLgyMjJITEwM2Pa8pFxCT6TkAcolVEVKLvvLIyPH8dikPWzc4zipUSwDmsQSF20eRFhygT4m/fr1m+qc61KilZ1zQXsAKcCcQq8XALX9z2sDC0qyndTUVBdIY8aMCej2vKRcQk+k5OGccglVkZJLcXls353jbv9shmt493fu6P+OcVOWby29wA5BoI8JMMWVsNaWdtf6N8Cl/ueXAsNLef8iIhIGKpeP5b9nd+D9K7qRlVvAeUN+Z8z8jV6HFZKCefvZUOB3oIWZrTazK4EngePMbBFwnP+1iIhIkfo0T2LEzb1pnlyRqz6YwkujF2mc9n0EbT5y59z5+3nrmGDtU0REIk/lhFjSBvXgni9n8+yohXw2dRVvX9qV5skVvQ4tJGhkNxERCXkV42N55cLOfHJVd7JzCzjnjd/57M9VGqsdFXIREQkjRzStwafX9KRJUiJ3fTmLM177jVmrt3sdlqdUyEVEJKw0qlGBL67tyXPndGDt9j2c9spE/v3VbPbklM1z5yrkIiISdsyMgZ3r8cvtfbiyVyPS/lzJeUN+J33dzr23O5cZKuQiIhK2KsbH8p9TWjPk4i4s2pjBiS+M5/bPZ5apYq5CLiIiYe+41sn8eMtRnNulPl9NW0PXwT/zx7J9p/uITCrkIiISEepXK88TA9vx7NkdqJQQy0VvTeal0YvIySvwOrSgUiEXEZGIERVlnJlaj8+v6clxrZN5dtRCBrw8gRGz15GdF5kXw6mQi4hIxKmeWI5XLuzMkItT2ZWVx3UfT6PvM2NZsinD69ACToVcREQi1vFtajHurn68e1lXcvMLOPv13/l8yqqIap2rkIuISESLjjL6tazJp9f0pG6VBO78Yha9nhzDm+OWkpcf/ufPVchFRKRMaJKUyDc39OLDK7vRqnZFBo9I5+w3fmf55t1eh3ZYVMhFRKTMMDN6N0viwyu789L5nVi8MYPjnv+VB4fPYeOuLK/DOyQq5CIiUiad2qEOP9/Wh7O71OfjySvp8/RYnvlxftgN9apCLiIiZVZypXgeP6MdP9/Wh+NaJ/PKmCUc8+xYPpy0ImwuiFMhFxGRMi+lRgVePL8TaYN6ULtKAvd/PYe+z4RHQVchFxER8evRuDpfXNuTj67sTp3CBf335SE7QpwKuYiISCFmxpHNavyzoA+fy6kvTWDUvA0UFITWhCwxXgcgIiISivYW9F5NqzM6fSMPfjOXqz+YQuMaFbjiyEac2bkeCXHRXoepFrmIiEhxzIxjWyfz6519efH8TiTGx/Cfr+dwxJOjeXH0InbsyfU0PrXIRURESiAmOooBHepwavvaTFmxjTd+XcJzoxby5vilHF3XOOooR1SUlX5cpb5HERGRMGZmdE2pRteUasxZs4OXf1nM8o2bPCnioK51ERGRQ9a2bmVevziVGzuV8ywGFXIREZHDFONRaxxUyEVERMKaJ+fIzWw5sAvIB/Kcc128iENERCTceXmxWz/n3GYP9y8iIhL21LUuIiISxrwq5A74ycymmtkgj2IQEREJe+Zc6Y8Za2Z1nHNrzawmMAq40Tk3bp91BgGDAJKTk1PT0tICtv+MjAwSExMDtj0vKZfQEyl5gHIJVZGSS6TkAYHPpV+/flNLfP2Yc87TB/AQcEdx66SmprpAGjNmTEC35yXlEnoiJQ/nlEuoipRcIiUP5wKfCzDFlbCOlnqL3MwqAFHOuV3+56OAR5xzI4v5zCZgRQDDqAFEyoV2yiX0REoeoFxCVaTkEil5QOBzaeicSyrJil5ctZ4MDDOzvfv/pLgiDlDSZErKzKa4CLnlTbmEnkjJA5RLqIqUXCIlD/A2l1Iv5M65pUCH0t6viIhIJNLtZyIiImGsrBbyIV4HEEDKJfRESh6gXEJVpOQSKXmAh7l4cvuZiIiIBEZZbZGLiIhEhDJZyM3sDjNzZlbD/9rM7EUzW2xms8yss9cxHoiZPeqPdYaZ/WRmdfzL+5rZDv/yGWb2gNexFqeYPMLxmDxjZvP98Q4zsyr+5SlmtqfQMXnd61gPZH+5+N/7t/+4LDCzE7yM80DM7Gwzm2tmBWbWpdDycDwmRebify9sjsm+zOwhM1tT6Fic5HVMB8vM+vt/94vN7J5SD6CkN5xHygOoD/yI7770Gv5lJwE/AAb0ACZ7HWcJ8qhU6PlNwOv+532B77yOLwB5hOMxOR6I8T9/CnjK/zwFmON1fAHKpTUwEygHNAKWANFex1tMHq2AFsBYoEuh5eF4TPaXS1gdkyLyeogDDAoWyg8g2v87bwzE+Y9F69KMoSy2yJ8H7sI33vtepwEfOJ9JQBUzq+1JdCXknNtZ6GUF/plP2Cgmj3A8Jj855/L8LycB9byM53AUk8tpQJpzLts5twxYDHTzIsaScM6lO+cWeB1HIBSTS1gdkwjUDVjsnFvqnMsB0vAdk1JTpgq5mQ0A1jjnZu7zVl1gVaHXq/3LQpqZDTazVcCFQOEu9J5mNtPMfjCzNh6FV2L7ySMsj0khV+DrUdirkZlNN7Nfzay3V0EdosK5hPtxKSycj0lhkXBMbvCfxnnHzKp6HcxB8vz37+V85EFhZj8DtYp46z7gXnxdhv/vY0Us87yFW1wuzrnhzrn7gPvM7N/ADcCDwDR8Q/tl+M81fQ00K7Wgi3CIeYTlMfGvcx+QB3zsf28d0MA5t8XMUoGvzazNPr0Rpe4Qcwm541KSPIoQtsekqI8Vsczz/yuFHeDv8mvAo/hifhR4Ft+Xx3Dh+e8/4gq5c+7YopabWTt8549mmm942HrANDPrhu8bVP1Cq9cD1gY51APaXy5F+AT4Hniw8B8i59wIM3vVzGo45zwbz/hQ8iBMj4mZXQqcAhzj/CfQnHPZQLb/+VQzWwI0B6YEOdxiHUouhOBxOYh/X4U/E5bHZD9C7pjsq6R5mdmbwHdBDifQPP/9l5mudefcbOdcTedcinMuBd8vv7Nzbj3wDXCJ/0rpHsAO59w6L+M9EDMr3MoeAMz3L69l/m8q/i8pUcCW0o+wZPaXB+F5TPoDdwMDnHOZhZYnmVm0/3ljfD0kS72JsmT2lwu+43KemZUzs0b4cvnDixgPRzgek2KE9THZ59qXM4A5XsVyiP4EmplZIzOLA87Dd0xKTcS1yA/RCHxXSS8GMoHLvQ2nRJ40sxZAAb4r8K/1Lz8L+JeZ5QF7gPMKtaZC0f7yCMdj8jK+K4dH+b9LTXLOXQscBTziPyb5wLXOua3ehVkiRebinJtrZp8B8/B1uV/vnMv3MM5imdkZwEtAEvC9mc1wzp1AGB6T/eUSbsekCE+bWUd83dHLgWu8DefgOOfyzOwGfHdDRQPvOOfmlmYMGtlNREQkjJWZrnUREZFIpEIuIiISxlTIRUREwpgKuYiISBhTIRcREQljKuQiIiJhTIVcREQkjKmQixwGM8v3z6E8x8y+tb/nIHdm9mGh9WLMbJOZfVdo2X3++aX3zsfefT/7yAhC3Ie9TTN7w8x6BSKeYvbx22F+/iEzu8P/fO+xmmlm08zsiMBEKeItFXKRw7PHOdfROdcW2Apc71++G2hrZgn+18cBa/Z+yMx64hvHvLNzrj1wLP+cQSkcdMc3xWnQOOcCWWz3HqsOwL+BJwK4bRHPqJCLBM7v/HP6wh+Ak/3PzweGFnqvNrDZP3kHzrnNzrliJ1ows4vM7A9/q/KNQmOFP2Vm1xVa7yEzu31/6x9gHx3MbJyZzTOzAn/PwsNFrNcKWLh3KFAzG2Nmx/mfP2ZmLx5oXyWxt+fAzFLMLN3M3vT3YvxU6EvSvp+5z8wW+GfcarGfTVcCtgUiRhGvqZCLBIC/SB7DPydLSMM3mUU80B6YXOi9n4D6ZrbQP0NdnwNsvxVwLtDLOdcR3/jgFxbaz7mFVj8H3yxe+1t/f/uIBz4F7nDOtQYGA/8FHipi9ROBkYVeP4hvKtoLgU7ArcXt6xA1A15xzrUBtgNn7ruC+aYkPc8fw0Cga6G3E/xfauYDb+GbMlMk7GnSFJHDk2BmM4AUYCowau8bzrlZZpaCrzU+ovCH/PPFpwK9gX7Ap2Z2j3Puvf3s5xggFfjTP4lJArDRv63pZlbTzOrgm1BjG9Buf+sX41hgmnNu78xZs4D++5l05wQKTWTjnBvnn3XvNqBvcZN2mNkAYJtzbvwB4tnXMufcDP/zqfh+5/vqDQzbO2ObmRX+YrXH/6Vm76mND8ysbYhPKiRyQCrkIodnj3Ouo5lVxjeP8vVA4W7lb/C1avsC1Qt/0F/sxgJjzWw2cCnw3n72Y8D7zrl/7+f9L/DNfFcLXwv9QOsXpS0wu9DrzsA0M0sCngbuBx7B19quUvhUgJm14+/TBbv8yy7D9yVlD7AOiPXvYwRQ4H+/D77pQwucc4MPEF92oef5+L6cFOWAhdk597uZ1cD3xedAX3BEQpq61kUCwDm3A7gJuMPMYgu99Q7wiHOucIHEzFrYP+di74hvGtf9GQ2cZWY1/Z+vZmYNC72fhq9L+Sx8Rf1A6xdlC75TAJhZc3xd02nOuU3ASuBZf45HAmMK5VIb+Bg4DdhtZicU2uaPe6dydc79B9882RUKvT/SOfcovgIfCOOAM8wswcwqAqcWtZKZtcQ35eSWAO1XxDNqkYsEiL+Leya+grp32WrghSJWTwRe8t+ulodv3vVBxWx7npn9B/jJzKKAXHyt/xX+9+f6C9ca59w6YF1x6+/HUGCAmc0BNgPnO+e2mFki0BjI858SOBHflwXMrDzwFXC7cy7dzB4FnsI3NzPATv/PTf6fOfjmOd9rt/+nFRNXiTnnppnZp8AMfLkW7r7fexpk7/4uDbN5u0WKpPnIRWS/zCwGGAI8jO8iuj+B54DuzrncA3z2Mnxd7d+ZWZpz7jz/Pd0xwBygxr7vBzMXkUilQi4iIhLGdI5cREQkjKmQi4iIhDEVchERkTCmQi4iIhLGVMhFRETCmAq5iIhIGFMhFxERCWMq5CIiImFMhVxERCSM/R/SmMAeB4gyTwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 576x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "w = 8  # wordlength of the quantizer\n",
    "A = np.logspace(-2, 0, num=500)  # relative RMS levels\n",
    "N = int(1e6)  # number of samples\n",
    "np.random.seed(1)\n",
    "\n",
    "def compute_SNR(a):\n",
    "    # compute input signal\n",
    "    x = np.random.laplace(size=N, scale=a/np.sqrt(2))\n",
    "    # quantize signal\n",
    "    xQ = uniform_midtread_quantizer(x, Q)\n",
    "    e = xQ - x\n",
    "    # compute SNR\n",
    "    SNR = 10*np.log10((np.var(x)/np.var(e)))\n",
    "    \n",
    "    return SNR\n",
    "\n",
    "\n",
    "# quantization step\n",
    "Q = 1/(2**(w-1))\n",
    "# compute SNR for given RMS levels\n",
    "SNR = [compute_SNR(a) for a in A]\n",
    "# plot results\n",
    "plot_SNR(A, SNR)\n",
    "# find maximum SNR\n",
    "Amax = A[np.argmax(SNR)]\n",
    "Pc = np.exp(-np.sqrt(2)/Amax)\n",
    "print(r'Maximum SNR = {0:2.3f} dB for A = {1:2.1f} dB with clipping probability {2:2.1e}'\n",
    "      .format(np.array(SNR).max(), 20*np.log10(Amax), Pc))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise**\n",
    "\n",
    "* Compare the SNR for the Laplace distributed signal to the case of a normally distributed signal. What is different?\n",
    "\n",
    "Solution: The overall SNR is lower compared to the case of a normally distributed signal. Its maximum is also at lower RMS levels. Both can be explained by the properties of the Laplace distribution discussed above."
   ]
  },
  {
   "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": {
  "anaconda-cloud": {},
  "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
}