random_signals_LTI_systems/correlation_functions.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "nbsphinx": "hidden"
   },
   "source": [
    "# Random Signals and LTI-Systems\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": [
    "## Auto-Correlation Function\n",
    "\n",
    "The auto-correlation function (ACF) $\\varphi_{yy}[\\kappa]$ of the output signal of an LTI system $y[k] = \\mathcal{H} \\{ x[k] \\}$ is derived. It is assumed that the input signal is a wide-sense stationary (WSS) real-valued random process and that the LTI system has a real-valued impulse repsonse $h[k] \\in \\mathbb{R}$. \n",
    "\n",
    "Introducing the output relation $y[k] = h[k] * x[k]$ of an LTI system into the definition of the ACF and rearranging terms yields\n",
    "\n",
    "\\begin{align}\n",
    "\\varphi_{yy}[\\kappa] &= E \\{ y[k+\\kappa] \\cdot y[k] \\} \\\\\n",
    "&= E \\left\\{  \\sum_{\\mu = -\\infty}^{\\infty} h[\\mu] \\; x[k+\\kappa-\\mu] \\cdot \n",
    "\\sum_{\\nu = -\\infty}^{\\infty} h[\\nu] \\; x[k-\\nu] \\right\\} \\\\\n",
    "&= \\underbrace{h[\\kappa] * h[-\\kappa]}_{\\varphi_{hh}[\\kappa]} * \\varphi_{xx}[\\kappa]\n",
    "\\end{align}\n",
    "\n",
    "where the ACF $\\varphi_{hh}[\\kappa]$ of the deterministic impulse response $h[k]$ is commonly termed as *filter ACF*. This is related to the [link between ACF and convolution](../random_signals/correlation_functions.ipynb#Definition). The relation above is known as the *Wiener-Lee theorem*. It states that the ACF of the output $\\varphi_{yy}[\\kappa]$ of an LTI system is given by the convolution of the input signal's ACF $\\varphi_{xx}[\\kappa]$ with the filter ACF $\\varphi_{hh}[\\kappa]$. For a system which just attenuates the input signal $y[k] = A \\cdot x[k]$ with $A \\in \\mathbb{R}$, the ACF at the output is given as $\\varphi_{yy}[\\kappa] = A^2 \\cdot \\varphi_{xx}[\\kappa]$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example - System Response to White Noise\n",
    "\n",
    "Let's assume that the wide-sense ergodic input signal $x[k]$ of an LTI system with impulse response $h[k] = \\text{rect}_N[k]$ is normal distributed white noise. Introducing $\\varphi_{xx}[\\kappa] = N_0\\, \\delta[\\kappa]$ and $h[k]$ into the Wiener-Lee theorem yields\n",
    "\n",
    "\\begin{equation}\n",
    "\\varphi_{yy}[\\kappa] = N_0 \\cdot \\varphi_{hh}[\\kappa] = N_0 \\cdot (\\text{rect}_N[\\kappa] * \\text{rect}_N[-\\kappa])\n",
    "\\end{equation}\n",
    "\n",
    "The example is evaluated numerically for $N_0 = 1$ and $N=5$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmQAAAGFCAYAAABT4e8GAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvFvnyVgAAIABJREFUeJzt3X2UXXV97/H3l0mAEdRpIKVkQEXFqdSA8aqVpbeOWBt8JGKvtVofWr3R29LqrQ0ltba2twptbntd1odK60NvLyKoIbW1NmJxltXFg2KQIDgVAaETFAEPkDBAMvneP86ecJjMwz7JmfObmfN+rTUrmXN+Z+/v/M5v7/M5e//2OZGZSJIkqZxDShcgSZLU6wxkkiRJhRnIJEmSCjOQSZIkFWYgkyRJKsxAJkmSVJiBTJIkqTADmSRJUmEGMmmBi4jvRMRw6TqmExGfjIg/K13HXCJiKCKuiYj7IuJ3StczH7o1TiLiloj4xQN4XEbEroh475Tbb4uINdO0vywiHoiIrx1MvdJiYSCT5kn1wjUeETtbfj5Y4zGPeLHLzJ/LzJF5qq/tF9YDWM9IRPwkIg6b5r7XRsQ3q765PSK+GBHPm1Lj1D5cdQBlnA18JTMfnZkfOJi/Z6pO9uPBLGu+xkmHnZKZ75r8JSJ+CjgWuGFqw8w8DXhbF2uTijKQSfPr5Zl5ZMvPWaUL6qaIeALwX4EEXjHlvt8F3g+8DzgGeBzwYeCMKYuZ2oc7DqCUxwPfOYDHaX6tBm7MzAdKFyKVZiCTCoiI34+IseoU2mhEvDAi/oFmKPmn6kjQ2VXbfUdNqv9viIhrq9M/H4uIY6ojS/dFxJerow6T6zknIr5f3Xd9RLyyun2mda2KiM9FxI8j4uapp/ciYk1EfKta3kXA4XP8qW8ArgA+CbyxZTmPBf4U+K3M3JyZuzJzd2b+U2ZuOID+fGp1JK5Rnbp7Rct9lwEvAD5Y/a1PafPxGRFPbvl932na6fqxeo42Vv39k4j4REQcPtfyZnpOpql1v7FT3f6Io2sR8YyI2Fa1+0xEXNRS9y0R8XvVOLqnuq+1xmnHzRzPwZERMRERx7bc9rTqyOejZ3jYycB1VdtHRcSnImJzRBw51/qkpcZAJnVZRAwBZwHPysxHA2uBWzLz9cCtPHxE6C9mWMSrgBcBTwFeDnwR+ANgJc1tujVEfZ/mEarHAn8C/L+IOHa6dUXEIcA/Ad8GBoEXAu+IiLVV3YcCW4B/AFYAn6lqmc0bgAuqn7URcUx1+6k0w9wlczx+ThGxvKr7S8BPA78NXFD18+Spr38Hzqr+1v9o5/GzmeU5ex3N5/VJNJ+nPzyIZbXWOu3YmabdoTT79pM0n6sLgamh6tXA6cAJNIPRm1rum3bczFH/TuC7wDNabj4PeF9m3jfDw1YD2yPiBODrwCjwqmpZUk8xkEnza0t11GXy578DE8BhwEkRsTwzb8nM77exzL/OzB9l5hjNoHFlZm6rTvtcAuybIJ2Zn8nMHZm5NzMvAr4HPHuG5T4LWJmZf5qZD2XmTcDfAq+p7n8OsBx4f3U067PAN2YqMppzwR4PXJyZV9N8kX9tdfdRwJ2ZuafG39vah1umuf85wJHAeVXdlwH/DPxqjWV34vHT+WBm3paZdwPvPchltao7dp4DLAM+UD1Xm4GrprT5QDU27qYZSJ8+eUeb46bVN6gCWUT8AnAS8NFZ2p9Mcw7ZV4D3ZOafZGbWWI+05BjIpPm1LjMHWn7+NjNvBN4BvAe4IyI+He1NVP9Ry//Hp/l93+meiHhDNK8ubEREA3gacPQMy308sKo1QNI88jZ5VGsVMDblBfMHs9T5RuBLmXln9funePi05V3A0RGxbJbHT2rtw3XT3L8KuC0z906pa7DGsjvx+OncNmVZB3Ihwn7aGDvTPVe3TWnzw5b/38+Bj5tW+wIZ8BfAuzPzoekaRkRUy30l8JHM/Mcay5eWLAOZVEBmfiozJ48gJfDnk3d1ah0R8XiaR7jOAo7KzAGa83VihnXdBtw8JUA+OjNfUt1/OzBYvZBOetwM6+6neUrs+RHxw4j4IfA/gVMi4hTgcuBBYLqA1a4dwPHVKdfWusY69Pj7gUe13PczUx4/3XN2/JRltV6IMNvy5nz+Zxk7raZ7ro6fpt1+aoyb2XwDeEZEvIrmKelPzdL2hOrfXwTeGRHPrFOftFQZyKQui+ZnYp0WzY+BeIDmUa3JozM/Ap7YoVUdQfMF+8fVen+d5hGJSVPXdRVwXzVpvD8i+qpJ2c+q7r8c2AP8TkQsj4gzmfk01jqap9dOonkq7OnAU2meYn1DZt4D/BHwoYhYV03oXh4RL46ImebOzeRKmiHn7GoZwzTn1n26Q4+/Bnht1R+nA8+f8vjpnrPfiojjImIF8C7gopb7ZlverM//HGOn1eU0+/+siFgWEWdQ75QjzD1uZvNtmgHzL4GNc5x+PBm4NjO3A+uBS+aapyYtZQYyaX5NXjE3+XMJzTlA5wF30jxt9NPAxqr9ucAfVqeKfu9gVpyZ19N8Ybyc5gv9apoTpyc9Yl2ZOQG8jGZ4urmq7+9oTuymOvV0Js3J33cDvwJsnmH1bwQ+kZm3ZuYPJ3+ADwKvi4hlmfmXwO/SnPD+Y5pH6M6ieeFAO3/nQzQD1Iurmj9MM/R9t0OPf3t1f4PmZP2p9U33nH2K5kUCN9GcO9f64bmzLW+u53+2sTP1bzoTeHO1nl+jOS/uwdn6onrsXONmtsc+CGyneZHKF+dovhq4tnrcFuB8mvMF57pyV1qSwvmTktQ5EXEL8JbM/HLpWlpFxJXA32TmJ+ZxHYcCNwKvzswrptz3AM1A+IHMfHeNZV1K8+KEqzLzhfNRr7SQ1JlQK0laZCLi+TQ/RuJOmkfiTgb+dZ5X+8fA16eGMYDMbOvIV2a+qGNVSYuAgUySlqYh4GKac8JuAn45M2+fjxVFxDNofnTFtez/eWeSavCUpSRJUmFO6pckSSrMQCZJklTYophDNjAwkE9+8pPnbqiO2bVrF0cccUTpMnqKfd599nn32efdZ59339VXX31nZq5s5zGLIpAdc8wxfPOb3yxdRk8ZGRlheHi4dBk9xT7vPvu8++zz7rPPuy8iZvtauWl5ylKSJKkwA5kkSVJhBjJJkqTCDGSSJEmFGcgkSZIKM5BJkiQVZiCTJEkqzEAmSZJUmIFMkiSpMAOZJElSYQYySZKkwgxkkiRJhRnIJEmSCjOQSZIkFWYgkyRJKsxAJkmSVJiBTJIkqTADmSRJUmEGMkmSpMIMZJIkSYUZyCRJkgozkEmSJBVmIJMkSSpsWakVR8QtwH3ABLAnM59ZqhZJkqSSigWyygsy887CNUjqMVu2jbFp6yhjjXEGr7iMDWuHWLdmsHRZknpY6UAmSV21ZdsYGzdvZ3z3BABjjXE2bt4OYCiTVExkZpkVR9wM/ARI4KOZef6U+9cD6wFWrlz5Xy6++OLuF9nDdu7cyZFHHlm6jJ5in3fHO0fu564H9t/vHXV48JfDjypQUW9xnHeffd59L3jBC65udypWyUA2mJljEfHTwKXAb2fmV6drOzQ0lKOjo90tsMeNjIwwPDxcuoyeYp93xwnnfIHp9noB3HzeS7tdTs9xnHeffd59EdF2ICt2lWVmjlX/3gFcAjy7VC2Seseqgf62bpekbigSyCLiiIh49OT/gV8CritRi6TesmHtEP3L+x5xW//yPjasHSpUkSSVm9R/DHBJREzW8KnM/NdCtUjqIZMT98/+7LU8NLGXwYF+r7KUVFyRQJaZNwGnlFi3JK1bM8iFV91Ko9Fg6++fVrocSfKT+iVJkkozkEmSJBVmIJMkSSrMQCZJklSYgUySJKkwA5kkSVJhBjJJkqTCDGSSJEmFGcgkSZIKM5BJkiQVZiCTJEkqzEAmSZJUmIFMkiSpMAOZJElSYQYySZKkwgxkkiRJhRnIJEmSCjOQSZIkFWYgkyRJKsxAJkmSVJiBTJIkqTADmSRJUmEGMkmSpMIMZJIkSYUZyCRJkgozkEmSJBVmIJMkSSrMQCZJklSYgUySJKkwA5kkSVJhBjJJkqTCDGSSJEmFGcgkSZIKM5BJkiQVZiCTJEkqzEAmSZJUmIFMkiSpMAOZJElSYQYySZKkwgxkkiRJhRnIJEmSCjOQSZIkFWYgkyRJKsxAJkmSVJiBTJIkqTADmSRJUmEGMkmSpMIMZJIkSYUZyCRJkgpbVmrFEdEHfBMYy8yXlapD0tKyZdsYm7aOsqMxzqqBfjasHWLdmsHiy5Kk2RQLZMDbgRuAxxSsQdISsmXbGBs3b2d89wQAY41xNm7eDtB2kOrksiRpLkVOWUbEccBLgb8rsX5JS9OmraP7AtSk8d0TbNo6WnRZkjSXUnPI3g+cDewttH5JS9COxnhbt3drWZI0l66fsoyIlwF3ZObVETE8S7v1wHqAlStXMjIy0p0CBcDOnTvt8y6zzw/eisODux7IaW+f2reNxjgTExMz9nk7y1J9jvPus88Xh8jcf4czryuMOBd4PbAHOJzmHLLNmflrMz1maGgoR0c9TdBNIyMjDA8Ply6jp9jnB2/qvC+A/uV9nHvm6v3mff3KRy+n0Wiw9fdffNDLUn2O8+6zz7svIq7OzGe285iun7LMzI2ZeVxmPgF4DXDZbGFMkupat2aQc89czaF9zV3b4ED/AQeoTi5LkuZS8ipLSeq4dWsGufCqWwG46K2nLphlSdJsigayzBwBRkrWIEmSVJqf1C9JklSYgUySJKkwA5kkSVJhBjJJkqTCDGSSJEmFGcgkSZIKM5BJkiQVZiCTJEkqzEAmSZJUmIFMkiSpMAOZJElSYQYySZKkwgxkkiRJhRnIJEmSCjOQSZIkFWYgkyRJKsxAJkmSVJiBTJIkqTADmSRJUmEGMkmSpMIMZJIkSYUZyCRJkgozkEmSJBVmIJMkSSrMQCZJklSYgUySJKkwA5kkSVJhBjJJkqTCDGSSJEmFGcgkSZIKM5BJkiQVZiCTJEkqzEAmSZJUmIFMkiSpMAOZJElSYQYySZKkwgxkkiRJhRnIJEmSCjOQSZIkFWYgkyRJKsxAJkmSVJiBTJIkqTADmSRJUmEGMkmSpMIMZJIkSYUZyCRJkgozkEmSJBVmIJMkSSrMQCZJklTYshIrjYjDga8Ch1U1fDYz/7hELZIWjy3bxti0dZQdjXFWDfSzYe0Q69YMWpOkRa9IIAMeBE7LzJ0RsRz4WkR8MTOvKFSPpAVuy7YxNm7ezvjuCQDGGuNs3LwdoFgAWog1SVqcipyyzKad1a/Lq58sUYukxWHT1tF9wWfS+O4JNm0dLVTRwqxJ0uJUbA5ZRPRFxDXAHcClmXllqVokLXw7GuNt3d4NC7EmSYtTqVOWZOYE8PSIGAAuiYinZeZ1k/dHxHpgPcDKlSsZGRkpU2iP2rlzp33eZfb57FYcHtz1wP4H0lccHvv1W6MKRLP1Z6MxzsTExJx9Ptuy2qlJTY7z7rPPF4digWxSZjYi4ivA6cB1LbefD5wPMDQ0lMPDw2UK7FEjIyPY591ln8/u3Y995HwtgP7lfbz7jNUMT5mv9ZHRywEYHj51xuV9ZPRyGo3GnH0+27LaqUlNjvPus88XhyKnLCNiZXVkjIjoB14EfLdELZIWh3VrBjn3zNUc2tfcbQ0O9HPumauLTp5fiDVJWpxKHSE7Fvj7iOijGQovzsx/LlSLpEVi3ZpBLrzqVgAueuvMR7+6aSHWJGnxqR3IImJFjWZ7M7MxV6PMvBZYU3fdkiRJS1k7R8h2VD8xS5s+4HEHVZEkSVKPaSeQ3ZCZsx7ViohtB1mPJElSz2lnUn+dyRFOoJAkSWpT7UCWmQ8ARMQZU++LiENa20iSJKm+A/nYi/UR8fOw79P2fwM/skKSJOmAHcjHXrwW+MeI+ALwP4DtwBs6WpUkSVIPOZBA9jzgXcAFwJsyc6SjFUmSJPWYAwlkrwJWA0cDH42IS4HrMvNvOlqZJElSj2g7kGXmbwBERAAn0gxnqztclyRJUs844K9OyswE/qP6+VzHKpIkSeoxta+yjIhvdaKNJEmSHqmdI2RPjYhrZ7k/gMceZD2SJEk9p51A9rM12kwcaCGSJEm9qnYgy8wfTHd7RByWmQ92riRJkqTeciCf1D/VhyPiRR1YjiRJUk866ECWmW8GHhcRH4iIoztQkyRJUk856EAWEWuBE4AnA38bEesOuipJkqQe0olTlk8APpaZL8nMVwIv6MAyJUmSekYnAtmzaR4dm/TuDixTkiSpZ3R8Dllm3tuBuiRJknqGc8gkSZIKO6BAFhHPiojnRsRjgZ8BPu4cMkmSpAPT9peLR8QW4DjgNuAk4OvAZ1qaOIdMkiSpDXMeIYuIn4uIC1puehpwFvArmTkEfA34yOSdziGTJElqT50jZF8GTm35/TzgbOCUiLgf2A48PyJ+CbgmM+/ofJmSJElLV505ZL8EvLfl95cDF2fmk4DnAX8D9AGvBr7Y8QolSZKWuDmPkGXmduB1LTe9BfhERGwErgOGgH/NzLfMT4mSJElLW9uT+jPzx8DLImKQ5nyyezLzio5XJkmS1CPaDmSTMnMMGOtgLZIkST2pE1+dJEmSpINgIJMkSSrMQCZJklSYgUySJKkwA5kkSVJhBjJJkqTCDGSSJEmFGcgkSZIKM5BJkiQVZiCTJEkqzEAmSZJUmIFMkiSpMAOZJElSYQYySZKkwgxkkiRJhRnIJEmSCjOQSZIkFWYgkyRJKmxZ6QIkacu2MTZtHWVHY5xVA/1sWDvEujWDpcvqmKX+90k6eAYySUVt2TbGxs3bGd89AcBYY5yNm7cDLInQstT/PkmdUeSUZUQcHxFfiYjrI+I7EfH2EnVIKm/T1tF9YWXS+O4JNm0dLVRRZy31v09SZ5Q6QrYHeGdmfisiHg1cHRGXZub1heqRVMiOxnhbty82S/3vk9QZRY6QZebtmfmt6v/3ATcAHruXetCqgf62bl9slvrfJ6kzil9lGRFPANYAV5atRFIJG9YO0b+87xG39S/vY8PaoUIVddZS//skdUbRSf0RcSTwOeAdmXnvlPvWA+sBVq5cycjISPcL7GE7d+60z7usV/t8AHj9U/v4+PYJ9iQcdXjwqqf0MXDP9xgZ+d5+7RvVqb7Z+qpum4mJiTn7/GDX1+7ft9T16jgvyT5fHIoFsohYTjOMXZCZm6fen5nnA+cDDA0N5fDwcHcL7HEjIyPY593Vy30+DHz7o5cDcNFbT5217UdGm+2Gh2duV7dNo9GYs887sb5h6v99S10vj/NS7PPFodRVlgF8DLghM/+qRA2SJEkLRak5ZM8FXg+cFhHXVD8vKVSLJElSUUVOWWbm14AosW5JkqSFpvhVlpIkSb3OQCZJklSYgUySJKkwA5kkSVJhBjJJkqTCDGSSJEmFGcgkSZIKM5BJkiQVZiCTJEkqzEAmSZJUmIFMkiSpMAOZJElSYQYySZKkwgxkkiRJhRnIJEmSCjOQSZIkFWYgkyRJKsxAJkmSVJiBTJIkqTADmSRJUmEGMkmSpMIMZJIkSYUZyCRJkgozkEmSJBVmIJMkSSrMQCZJklSYgUySJKkwA5kkSVJhBjJJkqTCDGSSJEmFGcgkSZIKM5BJkiQVZiCTJEkqzEAmSZJUmIFMkiSpMAOZJElSYQYySZKkwgxkkiRJhRnIJEmSCjOQSZIkFWYgkyRJKsxAJkmSVJiBTJIkqTADmSRJUmEGMkmSpMIMZJIkSYUZyCRJkgpbVroASUvblm1jbNo6yo7GOKsG+tmwdoh1awZLl7Xg2E9SbzOQSZo3W7aNsXHzdsZ3TwAw1hhn4+btAIaNFvaTpCKnLCPi4xFxR0RcV2L9krpj09bRfSFj0vjuCTZtHS1U0cJkP0kqNYfsk8DphdYtqUt2NMbbur1X2U+SigSyzPwqcHeJdUvqnlUD/W3d3qvsJ0leZSlp3mxYO0T/8r5H3Na/vI8Na4cKVbQw2U+SFuyk/ohYD6wHWLlyJSMjI2UL6jE7d+60z7tsKfb5APD6p/bx8e0T7Ek46vDgVU/pY+Ce7zEy8r1HtG1Up+fm6oM67eq2mZiY6Or6ZmrTTj8tdktxnC909vnisGADWWaeD5wPMDQ0lMPDw2UL6jEjIyPY5921VPt8GPj2Ry8H4KK3njpju4+MNtsMD8/cpm67um0ajcacfd7J9c3WZph6/bTYLdVxvpDZ54uDpywlSZIKK/WxFxcClwNDEfGfEfHmEnVIkiQtBEVOWWbmr5ZYryRJ0kLkKUtJkqTCDGSSJEmFGcgkSZIKM5BJkiQVZiCTJEkqzEAmSZJUmIFMkiSpMAOZJElSYQYySZKkwgxkkiRJhRnIJEmSCjOQSZIkFWYgkyRJKsxAJkmSVJiBTJIkqTADmSRJUmEGMkmSpMIMZJIkSYUZyCRJkgozkEmSJBVmIJMkSSrMQCZJklSYgUySJKkwA5kkSVJhBjJJkqTCDGSSJEmFGcgkSZIKM5BJkiQVZiCTJEkqzEAmSZJUmIFMkiSpMAOZJElSYQYySZKkwgxkkiRJhRnIJEmSCltWugB1z5ZtY2zaOsqOxjirBvrZsHaIdWsGp20z1hhn8IrLpm0jSVr46u7P67w2aP4ZyHrElm1jbNy8nfHdEwCMNcbZuHk7wL4Nr04bSdLCV3d/7n5/4fCUZY/YtHV03wY3aXz3BJu2jrbVZtKWbWM897zLOOGcL/Dc8y5jy7ax+SlckrSfufbBdffn7ez3Nb88QrZEzHXIeUdjfNrHtd5ep83kunxHJUll1NkH192ft7Pf97Tm/PII2RIwuXGONcZJHt44W98xrRron/axrbfXaQP131F5FE2S2teJo1919+d12tV5jdHBM5AtAXU2zg1rh+hf3veINv3L+9iwdqitNlDvHZUbsCS1r86+s84+uO7+vE47T2t2h4FsgatzlKnOxrluzSDnnrmaQ/uaT/ngQD/nnrn6EYec67SBeu+o3IAlqX2dOvpVd39ep107pzU9K3LgDGQLWN2jTHUPTa9bM8iaxw3w8yes4OvnnDbt+f/JNkM/dciMbeq8o6q7AUuSHtbJo1919uet7WZ6bfC0ZncYyBawukeZ6m6cnVLnHVXdkOg7Kkm9os7+rpNHvzrF05rd4VWWBXXiykh4+Kqasz97LQ9N7GWwC1fArFszyIVX3QrARW89db/7N6wdesRVQLD/BuzVmpJ6Rd39XZ195+RjZtsHd1Kd1xiv1jx4BrJC6mycqwb6GZtmkE/3DqqbG2cddTbg2d5RuYFKWkrq7u9KvMGuY67XmDqvV74Jn52BrJA6G2fdd0oL1VwbsO+oFq66fe5z033tfAWaz0v3dOqMByy8N9h11Hm98k347Axk86QTG+dCfafUKZ1+R+WLUD1z9ZNfubJwdfIr0Nxe6uvENtPOGY/FqJOnNaE3x6eT+udBpz6oFepdGblYdXKiaN0rfHr9IoI6/eRXrixcnfoKNK+Ia6qzP+jUNtPti69K6MTVmtC747NYIIuI0yNiNCJujIhzStUxH9w46+nk59/4ItTUiU/47vRXrqhzOvUVaH5vbf39Qae2mW5fGbkQ1X3d69U3e0VOWUZEH/Ah4EXAfwLfiIjPZ+b1JerpNE9H1teJiaJw8C9CS+FUTqe+365uny/1UzALUZ0+r9OmF763dq7tuO7+oJPbzGKcG9ZJdV/3evXNXqkjZM8GbszMmzLzIeDTwBmFauk4T0d2Tt13VHX6vNNf+bTQjhx06hO+O/mVK+qsTn0FWqe/t7bb5tr2OvX1Q9DZbUb1Xvfqjs+lJjKz+yuN+GXg9Mx8S/X764Gfz8yzpmv/lMcM5KVnvGLWZd5y1y4AnnDUEfPa5s6dD3LznbuY2JsctqyP41f0c/SRh+3X5qY7d7F378N9e8ghwROPPmK/ttfffi8AJx37mBlrqtOmk8u6/vZ72bNnDycfv6Jr65utzZ07H+T7P95F5sH1+bZbGzy455EvLgCHLetjzeMGarepu77JdnXGy1xt6rS74qa7pu0/gOc88ai2656rz9tp5zjvXJu628NsbeqOgzpjanJ53RrnC3lb79S2ULfdUh3nndy/Qr3X/rrtbrlrF41VJ3DG3/3vWZcVEVdn5jNnbTTFgr3KMiLWA+sBnnTEY2g0GrO2v/f+vQA0+nbPW5t7Hkp+uGsvkxn2wT0TfP/HO9l1//089tDY124ZcMyjgh/tSvYmLD8kOPpRwbI94zSmvPvqy2p9s/x9ddp0cll9uZdDDsmurm+2NsuAgcMAgmMeBUzTj3X6fMXhyQ93Qet7kIjm7ZPrnm4HPXl7a30/aOx9xM4CYO/e5Ad37mLZnub66oyXumOqTrvlhwS79+7/Bmv5IbGv9rpjs06ft9POcd65NnW3h9na1B0HdcZUt8d5nW2vznZcZ3/QTl91cluo226pjvM6fV53TEG91/667e69fy+3330vIyMjsy7rQJQKZGPA8S2/H1fdtk9mng+cD3DYsSfmbz73HbPO5XnfRy8HZj4vv2Xb2Jznredq89zzLpt2nsDgQD9fP+e0uf7mRWVkZIRThodLl9Fx080r+YWW5/g3az7H6875AtMdWw7g5vNeCtQbL3XHVJ12N28b4z3TfA7QuWeu5hRPh09rqY7zTqkzpro9zutse3W347n2B0vFUhzndcdUndf+dtpNZo23DXd+DmCpQPYN4MSIOIFmEHsN8NrZHnAwk0kn5xM8NLF3xmXVadOrEw2XknVrBmcdP3U/jLdTE6c7eUXj5N+1GC9I0MJUZ0x1e5zX2fba+foht4/FqZ05wbO9rrfTbr4VmdSfmXuAs4CtwA3AxZn5nbkeN9tl2dtubXDlzXcf8KX+nZoQrcVt8tL0wYF+gpkvTe/UxOm6Y6qdC0W+fs5p3HzeS71QRB0x15jq9jivs+3V3Y61eNUZK53+TMW5ssbBKvY5ZJn5L5n5lMx8Uma+t+7jZrose2qybfdqmjptvJKmN9QJNXV2+J26Iq6ddlK3dXuc1w1bvjlZ2uqMlU4ema2TNQ7Wgp3UP5N2Lstu5ysr6rTxlJBazXW6o854aW0z1hifce6CY08LVYlx7qlG1Rkgv26eAAAF8UlEQVQrnfxMxW58D+eiCmTTvaOqe2RrrvkEzjnQfKgzXibbjIyMMDzLxFvHnhYqx7lK6NSc4DrtujGHvMjnkLWr71GPzb4jVzw0sfPusb3j997det/ylU9YHX3LDp36mJzY89DuH9+yffL3Q/ofs6LvyBWD0bfs0JzYM+2y6rTpIUcDd5YuosfY591nn3effd59PdvndV/X52pXN2u0GMrMR7dT66I4QrZ3/N6rJ+6/p60PWNPBiYhvtvuhdjo49nn32efdZ593n33efRHxzXYfU2xSvyRJkpoMZJIkSYUtlkB2fukCepB93n32effZ591nn3effd59bff5opjUL0mStJQtliNkkiRJS9aCDmQR8b8i4tqIuCYivhQRq6rbIyI+EBE3Vvc/o3StS0VEbIqI71b9eklEDLTct7Hq89GIWFuyzqUkIv5bRHwnIvZGxDOn3Gefz5OIOL3q1xsj4pzS9SxFEfHxiLgjIq5ruW1FRFwaEd+r/v2pkjUuNRFxfER8JSKur/Yrb69ut9/nSUQcHhFXRcS3qz7/k+r2EyLiymofc1FE7PexGa0WdCADNmXmyZn5dOCfgT+qbn8xcGL1sx74SKH6lqJLgadl5snAfwAbASLiJJpfAv9zwOnAhyOib8alqB3XAWcCX2290T6fP1U/fojmvuQk4Fer/lZnfZLm2G11DvBvmXki8G/V7+qcPcA7M/Mk4DnAb1Vj236fPw8Cp2XmKcDTgdMj4jnAnwP/JzOfDPwEePNsC1nQgSwz72359QhgcsLbGcD/zaYrgIGIOLbrBS5Bmfml6svfAa4Ajqv+fwbw6cx8MDNvBm4Enl2ixqUmM2/IzNFp7rLP58+zgRsz86bMfAj4NM3+Vgdl5leBqR/CeQbw99X//x5Y19WilrjMvD0zv1X9/z7gBmAQ+33eVFlkZ/Xr8uongdOAz1a3z9nnCzqQAUTEeyPiNuB1PHyEbBC4raXZf1a3qbN+A/hi9X/7vPvs8/lj35ZzTGbeXv3/h8AxJYtZyiLiCcAa4Ers93kVEX0RcQ1wB80zTd8HGi0HOObcxxQPZBHx5Yi4bpqfMwAy812ZeTxwAXBW2WqXhrn6vGrzLpqHvi8oV+nSUafPpV6Tzcv8vdR/HkTEkcDngHdMOdtkv8+DzJyoplcdR/MI/M+2u4ziX52Umb9Ys+kFwL8AfwyMAce33HdcdZtqmKvPI+JNwMuAF+bDn4tinx+ENsZ5K/t8/ti35fwoIo7NzNurqSZ3lC5oqYmI5TTD2AWZubm62X7vgsxsRMRXgFNpTqdaVh0lm3MfU/wI2Wwi4sSWX88Avlv9//PAG6qrLZ8D3NNyKFYHISJOB84GXpGZ97fc9XngNRFxWEScQPOCiqtK1NhD7PP58w3gxOoqqENpXjzx+cI19YrPA2+s/v9G4B8L1rLkREQAHwNuyMy/arnLfp8nEbFy8hMJIqIfeBHNuXtfAX65ajZnny/oD4aNiM8BQ8Be4AfA2zJzrBpwH6R59c79wK9nZttf5Kn9RcSNwGHAXdVNV2Tm26r73kVzXtkemofBvzj9UtSOiHgl8NfASqABXJOZa6v77PN5EhEvAd4P9AEfz8z3Fi5pyYmIC4Fh4GjgRzTPcGwBLgYeR3O//urMnDrxXwcoIp4H/DuwneZrJ8Af0JxHZr/Pg4g4meak/T6aB7ouzsw/jYgn0rxgaAWwDfi1zHxwxuUs5EAmSZLUCxb0KUtJkqReYCCTJEkqzEAmSZJUmIFMkiSpMAOZJElSYQYySZKkwgxkkiRJhRnIJPWciPhERLw8IgYi4l+qD+eVpGIMZJJ60WrgJzS/yuTPMvOSwvVI6nF+Ur+knhIRhwD30fx6sA9l5p8XLkmSPEImqeecCOwA3gS8LSKWly1HkgxkknrPauDSzLwMuA54Q+F6JMlAJqnnrKYZxADeB2yMiGUF65Ek55BJkiSV5hEySZKkwgxkkiRJhRnIJEmSCjOQSZIkFWYgkyRJKsxAJkmSVJiBTJIkqTADmSRJUmH/H1/oZKJ9RwAzAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 720x432 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "L = 10000  # number of samples\n",
    "K = 30  # limit for lags in ACF\n",
    "\n",
    "# generate input signal (white Gaussian noise)\n",
    "np.random.seed(2)\n",
    "x = np.random.normal(size=L)\n",
    "# compute system response\n",
    "y = np.convolve(x, [1, 1, 1, 1, 1], mode='full')\n",
    "\n",
    "# compute and truncate ACF\n",
    "acf = 1/len(y) * np.correlate(y, y, mode='full')\n",
    "acf = acf[len(y)-K-1:len(y)+K-1]\n",
    "kappa = np.arange(-K, K)\n",
    "\n",
    "# plot ACF\n",
    "plt.figure(figsize = (10, 6))\n",
    "plt.stem(kappa, acf)\n",
    "plt.title('Estimated ACF of output signal $y[k]$')\n",
    "plt.ylabel(r'$\\hat{\\varphi}_{yy}[\\kappa]$')\n",
    "plt.xlabel(r'$\\kappa$')\n",
    "plt.axis([-K, K, 1.2*min(acf), 1.1*max(acf)]);\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise**\n",
    "\n",
    "* Derive the theoretic result for $\\varphi_{yy}[\\kappa]$ by calculating the filter-ACF $\\varphi_{hh}[\\kappa]$.\n",
    "* Why is the estimated ACF $\\hat{\\varphi}_{yy}[\\kappa]$ of the output signal not exactly equal to its theoretic result $\\varphi_{yy}[\\kappa]$?\n",
    "* Change the number of samples `L` and rerun the example. What changes?\n",
    "\n",
    "Solution: The filter-ACF is given by $\\varphi_{hh}[\\kappa] = \\text{rect}_N[\\kappa] * \\text{rect}_N[-\\kappa]$. The convolution of two rectangular signals $\\text{rect}_N[\\kappa]$ results in a triangular signal. Taking the time reversal into account yields\n",
    "\n",
    "\\begin{equation}\n",
    "\\varphi_{hh}[\\kappa] = \\begin{cases} \n",
    "N - |\\kappa| & \\text{for } -N < \\kappa \\leq N \\\\\n",
    "0 & \\text{otherwise}\n",
    "\\end{cases}\n",
    "\\end{equation}\n",
    "\n",
    "for even $N$. The estimated ACF $\\hat{\\varphi}_{yy}[\\kappa]$ differs from its theoretic value due to the statistical uncertainties when using random signals of finite length. Increasing its length `L` lowers the statistical uncertainties."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Cross-Correlation Function\n",
    "\n",
    "The cross-correlation functions (CCFs) $\\varphi_{xy}[\\kappa]$ and $\\varphi_{yx}[\\kappa]$ between the in- and output signal of an LTI system $y[k] = \\mathcal{H} \\{ x[k] \\}$ are derived. As for the ACF it is assumed that the input signal originates from a wide-sense stationary real-valued random process and that the LTI system's impulse response is real-valued, i.e. $h[k] \\in \\mathbb{R}$.\n",
    "\n",
    "Introducing the convolution into the definition of the CCF and rearranging the terms yields\n",
    "\n",
    "\\begin{align}\n",
    "\\varphi_{xy}[\\kappa] &= E \\{ x[k+\\kappa] \\cdot y[k] \\} \\\\\n",
    "&= E \\left\\{ x[k+\\kappa] \\cdot \\sum_{\\mu = -\\infty}^{\\infty} h[\\mu] \\; x[k-\\mu] \\right\\} \\\\\n",
    "&= \\sum_{\\mu = -\\infty}^{\\infty} h[\\mu] \\cdot E \\{ x[k+\\kappa] \\cdot x[k-\\mu] \\} \\\\\n",
    "&= h[-\\kappa] * \\varphi_{xx}[\\kappa]\n",
    "\\end{align}\n",
    "\n",
    "The CCF $\\varphi_{xy}[\\kappa]$ between in- and output is given as the time-reversed impulse response of the system convolved with the ACF of the input signal. \n",
    "\n",
    "The CCF between out- and input is yielded by taking the symmetry relations of the CCF and ACF into account\n",
    "\n",
    "\\begin{equation}\n",
    "\\varphi_{yx}[\\kappa] = \\varphi_{xy}[-\\kappa] = h[\\kappa] * \\varphi_{xx}[\\kappa]\n",
    "\\end{equation}\n",
    "\n",
    "The CCF $\\varphi_{yx}[\\kappa]$ between out- and input is given as the impulse response of the system convolved with the ACF of the input signal. \n",
    "\n",
    "For a system which just attenuates the input signal $y[k] = A \\cdot x[k]$, the CCFs between input and output are given as $\\varphi_{xy}[\\kappa] = A \\cdot \\varphi_{xx}[\\kappa]$ and $\\varphi_{yx}[\\kappa] = A \\cdot \\varphi_{xx}[\\kappa]$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## System Identification by Cross-Correlation\n",
    "\n",
    "The process of determining the impulse response or transfer function of a system is referred to as *system identification*. The CCFs of an LTI system play an important role in the estimation of the impulse response $h[k]$ of an unknown system. This is illustrated in the following.\n",
    "\n",
    "The basic idea is to use a specific measurement signal as input signal to the system. Let's assume that the unknown LTI system is excited by [white noise](../random_signals/white_noise.ipynb). The ACF of the wide-sense stationary input signal $x[k]$ is then given as $\\varphi_{xx}[\\kappa] = N_0 \\cdot \\delta[\\kappa]$. According to the relation derived above, the CCF between out- and input for this special choice of the input signal becomes\n",
    "\n",
    "\\begin{equation}\n",
    "\\varphi_{yx}[\\kappa] = h[\\kappa] * N_0 \\cdot \\delta[\\kappa] = N_0 \\cdot h[\\kappa]\n",
    "\\end{equation}\n",
    "\n",
    "For white noise as input signal $x[k]$, the impulse response of an LTI system can be estimated by estimating the CCF between its out- and input signals. Using noise as measurement signal instead of a Dirac impulse is beneficial since its [crest factor](https://en.wikipedia.org/wiki/Crest_factor) is limited."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example\n",
    "\n",
    "The application of the CCF to the identification of a system is demonstrated. The system is excited by wide-sense ergodic normal distributed white noise with $N_0 = 1$. The ACF of the in- and output, as well as the CCF between out- and input is estimated and plotted."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAnMAAADgCAYAAABo3eDKAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvFvnyVgAAIABJREFUeJzt3Xm0JHV99/H3Z2YAUUTUGQ3MgKCikYiKjigHY3jcQKJi1Ci4QVxGE01MTFTQRIkrhnOMJsGFqElcIhIXnER8cI+JjyijoAiKjoDAgGyCrLLN9/mj62LPne57u+/te7uL+36dc8/tqvp11bf2b/+qflWpKiRJktROy8YdgCRJkubOZE6SJKnFTOYkSZJazGROkiSpxUzmJEmSWsxkTpIkqcVM5iQNLMnvJjln3HH0kuSAJBeNO47ZpONfklyV5Ds9hj8vyRfHEdugkuyW5LokyxdhWpXk/gs9HanNTOakJSDJ+UlubE7AU3//NMD3tjiRVtX/VNUDFyjGf03y1oUYd9c0kuTcJGf3GX5gkm8kuTbJ5Un+O8nTmmFHJLlt2GXYw2OAJwJrqmrf6QOr6uNV9aQ5jHcoSXZv1u+KYb9bVRdU1Q5VddtCxCZpOEPvxJJa66lV9eVxBzFmjwXuBaxI8siqOm1qQJJnAR8GXg08FbgW+F3g+cD6pti3quox84zhPsD5VXX9PMcjSYA1c9KSl+T+TQ3Ur5JckeSTTf9vNEW+39RCPWf6pcymxu81SX6Q5PokH0py7yRfaGq3vpzk7l3l/yPJL5ppfSPJ7zT91wHPA17bTOs/m/67JPl0U0t2XpI/6xrX9k1t3lVNTdsjB5jdw4HPASc3n6fGFeBdwFuq6oNV9auq2lxV/11VL53DMt0lyfokv0yyMclLm/4vBj4I7NfM59/2+O4RSf63q7uSvDzJT5NcneS4Jt6pst9M8k/NMv1xksd3fff8JE/o6j46yceazqn1e3UTy349Ytk3yYYk1yS5NMm7mv5b1Ool2aOrRvPLTYwfm1b28CQXNNvYG6ZN41vNvF3SzMu2wy5zaSkzmZP0FuCLwN2BNcA/AlTVY5vhD20uqX2yz/efSeey4QPo1Gh9AXg9sIrOMebPusp+AdiTTu3Y94CPN9M6vvn8d820nppkGfCfwPeB1cDjgT9PcmAzrjcB92v+DqQrOeslyZ2BZzXT+ThwaFfS8EBgV+BTM41jCCcAFwG7NNN8e5LHVdWHgJfTqeHboareNOD4nkInWX0I8Gw68zvlUcDPgJV0lslnktxjgHFOrd+dmli+1aPMe4D3VNWOdJbziX3G9e/Ad4B7AkcDL+hR5jF0lvPjgTcmeVDT/zbgL5r492uG/8kA8UtqmMxJS8dJTe3H1N9UjdMtdC797VJVv66q/51hHL38Y1VdWlWbgP8Bvl1Vp1fVr4HPAvtMFayqD1fVtVV1E52T/kOT3K3PeB8JrKqqN1fVzVV1LvDPwKHN8GcDb6uqX1bVhcA/zBLnM4Cb6CSunwe2AX6/GXbP5v8ls4zj0dOW4aOnF0iyK7A/8LpmeZ5BpzbuhbOMeybHVNXVVXUB8DXgYV3DLgPeXVW3NAn3OV3zNV+3APdPsrKqrquqU6cXSLIbnXX1xmY9/S+/uSzd7W+r6saq+j6dBP2hAFX13ao6tapurarzgQ8Avzei+KUlwWROWjqeXlU7df39c9P/tUCA7yQ5K8mLhhzvpV2fb+zRvQNAkuVJjknysyTXAOc3ZVb2Ge99gF26kyc6NX73bobvAlzYVf7ns8R5OHBikzT8Gvg0v6nNu7L5v/Ms4zh12jLcKrlp4vplVV07LbbVs4x7Jr/o+nwDzTJtbKqqmjatXeYxrW4vplPj+uMkpyV5So8yU/N7Q1e/C3uU6zkPSR6Q5L+ay+/XAG+n/zYhqQeTOWmJq6pfVNVLq2oX4GXAe7Mwj4J4LnAI8ATgbsDuTf9MhTKt/IXAedOSp7tW1cHN8EvoXBqdslu/CSdZAzwOeH6TNPyCzuXPg5OspFObdSGdS8bzdTFwjyR3nRbbphGMu5fVU/fQdU3r4ubz9cCdu4b9Vtfn6ct7K1X106o6jM5l8XcCn0pyl2nFLqEzv93T2ZXBvQ/4MbBnczn39fxmm5A0AJM5aYlL8odNsgNwFZ2T/Oam+1LgviOa1F3pXOa8kk6C8fZpw6dP6zvAtUle1zR2WJ7kwUmmGjqcCByV5O5N/H86w7RfAPyEzj1bD2v+HkDnvrbDmpqtVwN/k+SPkuyYZFmSxyQ5fpiZbC75/j/gHUnulOQhdGq4PjbzN+fsXsCfJdkmyR8CD6LTwAPgDDr3Bm6TZC2dBHbK5XTWc9/1m+T5SVZV1Wbg6qb35u4yVfVzYANwdJJtm4YUTx0i/rsC1wDXJflt4I+H+K4kTOakpeQ/s+Uz0j7b9H8k8O0k19G51+lVzf1p0Lmv7d+ay5zPnuf0P0LnEuAm4Gxg+iXKDwF7NdM6qXmG2VPoJF7nAVfQufds6h67v23Gdx6d++A+OsO0Dwfe29RC3v4HvL8ZRlV9CngO8CI6NVuXAm+l0/p1WIfRqXm8mM59g29awMfCfJtOo5IrgLcBz6qqqcvGf0On4cJVdJbXv099qbks+jbgm/3u/wMOAs5qto33AIdW1Y09yj2PTuOFK+kss0/SSdwH8Vd0am2vpXNPZL+GNpL6yJa3WkiS2iLJEcBLRvDsu5FK5/E2Px6ita6kebBmTpI0L0kemeR+zaXpg+jcG3nSuOOSlgrfACFJmq/fAj5D5xEvFwF/XFWnjzckaenwMqskSVKLeZlVkiSpxUzmJEmSWmxJ3TO3cuXK2n333ccdhiRJ0qy++93vXlFVq2Yrt6SSud13350NGzaMOwxJkqRZJZntNYWAl1klSZJazWROkiSpxUzmJEmSWsxkTpIkqcVM5iRJklrMZE6SJKnFTOYkSZJabCKTuSQfTnJZkh/2GZ4k/5BkY5IfJHn4YscoSSedvon9j/kqexz5efY/5qucdPqmcYckaQmayGQO+FfgoBmGPxnYs/lbB7xvEWKSpNuddPomjvrMmWy6+kYK2HT1jRz1mTNN6CQtuolM5qrqG8AvZyhyCPCR6jgV2CnJzosTnSTBsaecw4233LZFvxtvuY1jTzlnTBFJWqomMpkbwGrgwq7ui5p+W0myLsmGJBsuv/zyRQlO0h3fxVffOFR/SVoobU3mBlZVx1fV2qpau2rVrO+qlaSB7LLT9kP1l6SF0tZkbhOwa1f3mqafJC2K1xz4QLbfZvkW/bbfZjmvOfCBY4pI0lLV1mRuPfDCplXro4FfVdUl4w5K0tLx9H1W845n7M22yzuH0dU7bc87nrE3T9+n5x0fkrRgVow7gF6SfAI4AFiZ5CLgTcA2AFX1fuBk4GBgI3AD8EfjiVTSUvb0fVbzie9cAMAnX7bfmKORtFRNZDJXVYfNMryAVyxSOJIkSROrrZdZJUmShMmcJElSq5nMSZIktZjJnCRJUouZzEmSJLWYyZwkSVKLmcxJkiS1mMmcJElSi5nMSZIktZjJnCRJUouZzEmSJLWYyZwkSVKLmcxJkiS1mMmcJElSi5nMSZIktZjJnCRJUotNbDKX5KAk5yTZmOTIHsN3S/K1JKcn+UGSg8cRpyRJ0jhNZDKXZDlwHPBkYC/gsCR7TSv218CJVbUPcCjw3sWNUpIkafwmMpkD9gU2VtW5VXUzcAJwyLQyBezYfL4bcPEixidJkjQRVow7gD5WAxd2dV8EPGpamaOBLyb5U+AuwBMWJzRJkqTJMak1c4M4DPjXqloDHAx8NMlW85NkXZINSTZcfvnlix6kJEnSQprUZG4TsGtX95qmX7cXAycCVNW3gDsBK6ePqKqOr6q1VbV21apVCxSuJEnSeExqMncasGeSPZJsS6eBw/ppZS4AHg+Q5EF0kjmr3iRJ0pIykclcVd0KvBI4BfgRnVarZyV5c5KnNcX+Enhpku8DnwCOqKoaT8SSJEnjMakNIKiqk4GTp/V7Y9fns4H9FzsuSZKkSTKRNXOSJEkajMmcJElSi5nMSZIktZjJnCRJUouZzEmSJLWYyZwkSVKLmcxJkiS1mMmcJElSi5nMSZIktZjJnCRJUouZzEmSJLWYyZwkSVKLmcxJkiS1mMmcJElSi5nMSZIktZjJnCRJUotNZDKX5KAk5yTZmOTIPmWeneTsJGcl+ffFjlGSJGkSrBh3ANMlWQ4cBzwRuAg4Lcn6qjq7q8yewFHA/lV1VZJ7jSdaSZKk8ZrEmrl9gY1VdW5V3QycABwyrcxLgeOq6iqAqrpskWOUJEmaCJOYzK0GLuzqvqjp1+0BwAOSfDPJqUkOWrToJEmSJsjEXWYd0ApgT+AAYA3wjSR7V9XV0wsmWQesA9htt90WM0ZJkqQFN4k1c5uAXbu61zT9ul0ErK+qW6rqPOAndJK7rVTV8VW1tqrWrlq1akECliRJGpdJTOZOA/ZMskeSbYFDgfXTypxEp1aOJCvpXHY9dzGDlCRJmgQTl8xV1a3AK4FTgB8BJ1bVWUnenORpTbFTgCuTnA18DXhNVV05noglSZLGZyLvmauqk4GTp/V7Y9fnAl7d/EmSJC1ZAyVzSe4xQLHNvRogSJIkaeEMWjN3cfOXGcosB2wuKkmStIgGTeZ+VFX7zFQgyekjiEeSJElDGLQBxH4jKiNJkqQRGiiZq6pfAySZ/lotkizrLiNJkqTFM+yjSdYleRRAkuVJXgT8ePRhSZIkaRDDPprkucDnknwe+GPgTOCFI49KkiRJAxk2mXsM8Abg48ARVfX1kUckSZKkgQ2bzD0T2BtYCXwgyZeAH1bV+0cemSRJkmY1VDJXVS8CSBI6L7bfu/mTJEnSGMzpdV7N67R+0vx9eqQRSZIkaWADtWZN8r1RlJEkSdJoDVoz96AkP5hheIC7jSAeSZIkDWHQZO63Byhz23wCkSRJ0vAGSuaq6ue9+idZXlUmcZIkSWMy7Bsgpjs+yZ0Bkjx2BPFIkiRpCPNN5t4IfCjJR4FHjiCe2yU5KMk5STYmOXKGcs9MUknWjnL6kiRJbTDfZO4twDlAASfOP5yOJMuB44AnA3sBhyXZq0e5uwKvAr49qmlLkiS1yXyTuddW1dF03tP6pvmHc7t9gY1VdW5V3QycABzSo9xbgHcCvx7htCVJklpjvsncMUm2r6rrgU+MIqDGauDCru6Lmn63S/JwYNeq+vwIpytJktQqc3oDRJc3AR9OcitwBvCV+Yc0uyTLgHcBRwxQdh2wDmC33XZb2MAkSZIW2UTeMwdsAnbt6l7T9JtyV+DBwNeTnA88GljfqxFEVR1fVWurau2qVatGGKIkSdL4DV0zl+SRwLbAD+ncM3dFkrsA7wFeMqK4TgP2TLIHnSTuUOC5UwOr6lfAyq6Yvg78VVVtGNH0JUmSWmGoZC7JSXRqyS6k08r0m0leWVXXJ3nZqIKqqluTvBI4BVgOfLiqzkryZmBDVa0f1bQkSZLabMZkLsnvAK+vquc1vR4MPB/4XlXdnORFwPuAw0f9JoiqOhk4eVq/N/Ype8Aopy1JktQWs9XMfRnYr6v7GOC1wEOT3ACcCfxekicBZ1TVZQsTpiRJknqZrQHEk4C3dXU/FTixqu4HPAZ4P53LoM8GvrAgEUqSJKmvGWvmqupM4HldvV4C/EuSo+g0gHgg8H+ralQNHyRJkjSEoRpAVNXlwFOSrKZz/9yvqurUBYlMkiRJs5rTQ4OrahNbPvdNkiRJYzDfhwZLkiRpjEzmJEmSWsxkTpIkqcVM5iRJklrMZE6SJKnFTOYkSZJazGROkiSpxUzmJEmSWsxkTpIkqcVM5iRJklrMZE6SJKnFJjKZS3JQknOSbExyZI/hr05ydpIfJPlKkvuMI05JkqRxm7hkLsly4DjgycBewGFJ9ppW7HRgbVU9BPgU8HeLG6UkSdJkmLhkDtgX2FhV51bVzcAJwCHdBarqa1V1Q9N5KrBmkWOUJEmaCJOYzK0GLuzqvqjp18+LgS8saESSJEkTasW4A5iPJM8H1gK/N0OZdcA6gN12222RIpMkSVock1gztwnYtat7TdNvC0meALwBeFpV3dRvZFV1fFWtraq1q1atGnmwkiRJ4zSJydxpwJ5J9kiyLXAosL67QJJ9gA/QSeQuG0OMkiRJE2HikrmquhV4JXAK8CPgxKo6K8mbkzytKXYssAPwH0nOSLK+z+gkSZLu0CbynrmqOhk4eVq/N3Z9fsKiByVJkjSBJq5mTpIkSYMzmZMkSWoxkzlJkqQWM5mTJElqMZM5SZKkFjOZkyRJajGTOUmSpBYzmZMkSWoxkzlJkqQWM5mTJElqMZM5SZKkFjOZkyRJajGTOUmSpBZbMe4ApKXspNM3cewp53Dx1Teyy07b85oDH8jT91k97rCkOzz3Pd2RmMxJY3LS6Zs46jNncuMttwGw6eobOeozZ94+3BONNBrTE7f/89ur+PR3N/Xc99zP1EYmc9P4a+03hlkWS3G5zXeejz3lnNtPJlNuvOU2jl5/FjfdunmrE82Gn/+Sr/348q2m1yuOqfEPUnam/gsx3/0sxW1osS3kMp7vNgTz32Z7jQPY6kfTx0+9gJoW04233Maxp5wz1PJYitvssPO8FJdRPwu5LFI1fZOeDEkOAt4DLAc+WFXHTBu+HfAR4BHAlcBzqur8mca53c571tpXfWDgnR5g+22W845n7N1zgU/CCbCf+U5veq0R9F8WM5WFhalhGsXyHOaE0uu7w2wrvexx5Oe3OqHMJLBF+e23Wc4zH7F6ixoGgG2WBQK33Fazlp2p/7DretD57rfcF2obWqhtZWocz/nAtwD45Mv2W5TpzXUc02ujYLTLeD7Hi1Fss/3GcadtlnHVDbcMNB8Bzjvm90c+z1Pl5/tDaqZYFuM4O9M2tFDHyUk+dw5btteyeOYjVvf8kT4lyXerau1scU9kMpdkOfAT4InARcBpwGFVdXZXmT8BHlJVL09yKPAHVfWcmca73c571s6Hv3vonX71TtvzzSMft0W/USU7i73DQe8kdnq/Y085h01X37jVtHbafhvust2Kgct21zDNNt9z/aU9Nd7Zdorp0xr0hNIr3v2P+WrPee61rfTTbxzDWJ5w24D7cL+y/foPs64Hne9++0K/fa/fNtRvXY8igRn2BNadzA37w6Zfv4UYx/QfA7Mt42GOT/225WG2oV6G3WbnaxT77zDni2F+SE2NZ77Hw0ETkF4x99uGVnet27meR0a1/87HQp7X+20vvX6kd4+j7cncfsDRVXVg030UQFW9o6vMKU2ZbyVZAfwCWFUzzNBUMjcXq3fafuCT2vSN+oabb+2bJPbaAYa9NDGfg3a/BGb65b9+hik7ZdCD3bBJd7+dAgZPVvvFO3099fvuTL/sBz1QDVOTsJhmW9fT95FhDmjDGrSWctgEZphxTG3H3cncTEnNoPveMIntKGqj+umX2PQ6Dv3FJ88YqJZ5LseLhTLM8aLXttyvZj3A3z/nYQOdA4b5IQW9k7ZhjofDJI/D7qvT1+2w55FR7b+D/qCHrbflYc7Vw57XB91HpsYxte+1PZl7FnBQVb2k6X4B8KiqemVXmR82ZS5qun/WlLmi33jvt8Pd6u0Pfcy841u2LGze3H+5zTZ8prLLloVVO2zHFdfdxG2bi+1WLGfXe2wPwLlXXL9V2WUJt962eY5z0l8SBt02hik7ZfmybDF/F/7yRm66dfQH+RXLl7G5aqvlNuj6GfY7U/Nz3hXXD7T+Vu2wHZddexNVM5ftZxTraSHW9bJl4b4r7wKwxbJYiHU8bGyjHMej73tPzr7kGgD22nlHTj33ynlNf9Lc/147bLH+drrzNlx+3U3zOg4t9jbb7xgw6L7X75jc75jVa3rzNapjfb/ltmL5Mqpqi/nbeNl18x7vQp9HZtPvODTscXZqXPM5r8+2/s6922o+8JBDgC0rBQZN5u7wDSCSrAPWAdz3LjvOWHbQnXDz5ppx4x10hfcqu3lzcek1v769+6Zbb+PcK65nWZ+ym4e662pwVTXwxtur7Gwb721N2an5G+WBr1uv6c+0/noZdJ0uWxZ2uvM2W8zPbOvv6htu4d47bgfA7ve8yxbDBzmJrtphu636JwHYYv76le3Xv59Bt4vNm4vzr7xhi/1ppkSu3wl3mBPYKE4Ew4xjuxXLAbjztsu36LeQCeuozbSMVyxfttW2fOk1W8/b5s3FsuWZ1/FiFNtsv3Hsfs87A1ufzFfusB3NV27f906/4OqBj8mdGDZvFdvUdwYx6HFoVMf6ftPqXv9T87di+bKB9r2Z1vuw55FR63cc6ndM7mcU5/Vh9pFddtp+oGl1m9RkbhOwa1f3mqZfrzIXNZdZ70anIcQWqup44HjoXGZ93e/+yYz3Ri2HgS6nTX1nPtXKC3m5oVcV+zCXXkZRrQzw1wNeAh7mHpheVez9xjuT6etgLpecp19afMMQl2+n4u51WfY+dFr2dOt1eeuAIVr19Svbq/8w63rYy6b9Lm9N3/dGtQ0Nsy8MM453PGNv7rPPau7T1f/00zdx9BCXwnrpdwlpFOPodRkKtl7Gc2k4MOilxZmOF/PZZvuN4xHNZbbp+xOwxboDOGCIRknzvZzW79LisIY5Hs73ONtvGxrFbUej2H8XwqjO6732kX632kxty8OY1MusK+g0gHg8naTtNOC5VXVWV5lXAHt3NYB4RlU9e6bxztaadZh7fEZxw+cw9231M+xBez43/c+lIcf0A+4wyfFMsQHzuv9smBPKMDf9D9tCdZgbrhfTMOt6LvfBDXJ/XXcsg6zrmVo+DrIvDDuO+d7bOsz2PapxDBPzXO/x6R7nQjb8WgjDbMv9fogN0xik103//ZKdmRqqwGDHw2GSx17Jxyha9w/bGGTQ/XcuP+h7GbaR33wbic3WIKXV98wBJDkYeDedR5N8uKreluTNwIaqWp/kTsBHgX2AXwKHVtW5M41z7dq1tWHDhqHimEvz8/ls1HNpcQbzO2jP9/vDHJiHTY5HMW+jaBo/38RmFC0GF9t8Wr0N2zp8VLGN4rECi/kcPRiuRd4oxjGouba+my3eSd3eYfgWnP1OzvM55ozqUU+Dbt8z1aAOs58u9r433wZlwxyTx/n0itYncwthLskcLOxDEsfdFHuxjePX+mI+x2kcz9ybBAuVSGt85vpcrLab77PVeo1jHM9GHGZad5T9dCGfYTmu58qazPUw12RusbXt1+ywnL+lw2XRbq6/jjv6cnD+JpfJXA9tSeYkSZIGTeaWLUYwkiRJWhgmc5IkSS1mMidJktRiJnOSJEktZjInSZLUYiZzkiRJLWYyJ0mS1GImc5IkSS1mMidJktRiJnOSJEktZjInSZLUYiZzkiRJLWYyJ0mS1GImc5IkSS02cclcknsk+VKSnzb/796jzMOSfCvJWUl+kOQ544hVkiRp3CYumQOOBL5SVXsCX2m6p7sBeGFV/Q5wEPDuJDstYoySJEkTYRKTuUOAf2s+/xvw9OkFquonVfXT5vPFwGXAqkWLUJIkaUJMYjJ376q6pPn8C+DeMxVOsi+wLfCzhQ5MkiRp0qwYx0STfBn4rR6D3tDdUVWVpGYYz87AR4HDq2pznzLrgHVN53VJzplb1K2wErhi3EFoTlx37eb6ay/XXbvd0dfffQYplKq+udJYNMnWAVV1SZOsfb2qHtij3I7A14G3V9WnFjnMiZRkQ1WtHXccGp7rrt1cf+3lums311/HJF5mXQ8c3nw+HPjc9AJJtgU+C3zERE6SJC1lk5jMHQM8MclPgSc03SRZm+SDTZlnA48FjkhyRvP3sPGEK0mSND5juWduJlV1JfD4Hv03AC9pPn8M+Ngih9YGx487AM2Z667dXH/t5bprN9cfE3jPnCRJkgY3iZdZJUmSNCCTuTuIJH+ZpJKsbLqT5B+SbGxeefbwcceorSU5NsmPm3X02e43mSQ5qll/5yQ5cJxxqrckBzXrZ2OSXm+r0QRJsmuSryU5u3kd5Kua/rO+RlKTIcnyJKcn+a+me48k3272wU82DSSXHJO5O4AkuwJPAi7o6v1kYM/mbx3wvjGEptl9CXhwVT0E+AlwFECSvYBDgalX1r03yfKxRamtNOvjODr72l7AYc160+S6FfjLqtoLeDTwimadDfIaSU2GVwE/6up+J/D3VXV/4CrgxWOJasxM5u4Y/h54LdB9A+QhdB7dUlV1KrBT89w+TZCq+mJV3dp0ngqsaT4fApxQVTdV1XnARmDfccSovvYFNlbVuVV1M3ACnfWmCVVVl1TV95rP19JJClYzwGskNX5J1gC/D3yw6Q7wOGDqEWVLdt2ZzLVckkOATVX1/WmDVgMXdnVf1PTT5HoR8IXms+tv8rmOWizJ7sA+wLcZ8jWSGpt306m4mHrj0z2Bq7t+EC/ZfXDiHk2irc3y+rPX07nEqgk10/qrqs81Zd5A5xLQxxczNmkpSrID8Gngz6vqmk4FT8dsr5HUeCR5CnBZVX03yQHjjmfSmMy1QFU9oVf/JHsDewDfbw5Ga4DvJdkX2ATs2lV8TdNPi6zf+puS5AjgKcDj6zfPCnL9TT7XUQsl2YZOIvfxqvpM0/vSJDt3vUbysvFFqD72B56W5GDgTsCOwHvo3EK0oqmdW7L7oJdZW6yqzqyqe1XV7lW1O50q5odX1S/ovBbthU2r1kcDv+q6jKAJkeQgOpcNnlZVN3QNWg8cmmS7JHvQacjynXHEqL5OA/ZsWtNtS6fByvoxx6QZNPdYfQj4UVW9q2vQrK+R1HhV1VFVtaY51x0KfLWqngd8DXhWU2zJrjtr5u64TgYOpnPj/A3AH403HPXxT8B2wJea2tVTq+rlVXVWkhOBs+lcfn1FVd02xjg1TVXdmuSVwCnAcuDDVXXWmMPSzPYHXgCcmeSMpt/r6bw28sQkLwZ+TueVkWqH1wEnJHkrcDqdZH3J8Q0QkiRJLeZlVkmSpBYzmZMkSWoxkzlJkqQWM5mTJElqMZM5SZKkFjOZkyRJajGTOUmSpBYzmZOkOUryL0memmSnJCcn+YNxxyRp6TGZk6S52xu4is4rhN5aVZ8dczySliDfACFSD0BDAAAAk0lEQVRJc5BkGXAtcCVwXFW9c8whSVqirJmTpLnZE7gYOAJ4eZJtxhuOpKXKZE6S5mZv4EtV9VXgh8ALxxyPpCXKZE6S5mZvOkkcwNuBo5KsGGM8kpYo75mTJElqMWvmJEmSWsxkTpIkqcVM5iRJklrMZE6SJKnFTOYkSZJazGROkiSpxUzmJEmSWsxkTpIkqcX+Pw8QIyxNqlZwAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 720x216 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmsAAADgCAYAAABcpCHWAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvFvnyVgAAIABJREFUeJzt3X+cXHV97/HXO5slLEZZIFtNNoHQmm4vGjUaqVzsvVZrQ6lAClSwVkCxgVaueq831NhqkVrBm3tta/1VBCpaLgQhxthCUyy0Sq/8CAQIAWMj8iMbfoQfSwgs+bH53D/OmTCZnJmd2Z2dc2bn/Xw85pGd7/nmzOf82D2f+f44RxGBmZmZmRXTlLwDMDMzM7PqnKyZmZmZFZiTNTMzM7MCc7JmZmZmVmBO1szMzMwKzMmamZmZWYE5WTPrUJJ+TdLGvOPIIukdkjbnHcdolPg7Sc9Kuj3veCaKpE9JurQFn3OWpFsm+nPM2o2TNbM2I+khScOStpe9vlzH/wtJry29j4gfRcTABMX4TUmfm4h1l32GJD0o6f4qyxdJ+qGk5yVtlfRvkk5Ml50laaTRfZjh7cC7gdkRcfQ4Nicr/gsk/X0R1hcRn4+IDzcrFjNrzNS8AzCzMTkhIn6QdxA5+y/ALwBTJb01Iu4oLZB0KnA58D+AE4DngV8Dfh9YnVb7cUS8fZwxHAE8FBEvjHM9ZmZVuWXNbBKR9Nq0Bek5SU9JWpGW/zCtck/ainRaZVdj2mK3VNK9kl6QdJmkV0u6IW2d+oGkQ8rqf0fS4+ln/VDS69LyJcD7gfPTz/p+Wj5L0nVpK9fPJX20bF09aWvcs2lL2Vvr2Nwzge8B16c/l9Yl4IvAn0fEpRHxXETsiYh/i4g/GMM+nSVptaRnJG2S9Adp+dnApcAx6XZ+NuP/TpH0p5IelvSkpG9JOjhdtl9Xb3oMfkPSccCngNPSdd+TLv9XSRdJul3SNknfk3ToWNeXEe8fSxpMj/dGSe9Ky/dplZN0RrpNT0v6dOlzyupek27r85I2SFpY9n8/Keln6bL7Jf1Oo8fErNM4WTObXP4c+GfgEGA28DcAEfFf0uVvjIjpEbGiyv8/haRb75dJWqRuILnI95H8vfhoWd0bgHkkrVt3AVemn3VJ+vP/Sj/rBElTgO8D9wD9wLuAj0talK7rz4BfSl+LKEu+skg6CDg1/ZwrgdMlHZAuHgDmANfWWkcDrgY2A7PSz/y8pHdGxGXAuSQtdNMj4s8y/u9Z6evXgV8EpgOjdrdGxD8BnwdWpOt+Y9niM4APATOB3cCXxrk+ACQNAOcBb42IV5Ich4cy6h0FfJUkIZ8JHExyTMudSLLfeklaMsu3+WckrZwHA58F/l7SzNG2wayTOVkza0+rJA2VvUotRrtIuuZmRcRLEdHoYO2/iYgnImIQ+BFwW0Ssi4iXgO8CC0oVI+LyiHg+InYAFwBvLLUaZXgr0BcRF0bEzoh4EPgGcHq6/L3AX0TEMxHxKKMnICcDO0gS038EuoHfTpcdlv772CjreFvFPnxbZQVJc4BjgT9O9+fdJK1pZ4yy7pL3A1+MiAcjYjuwjCSxHM8QlG9HxH1p1+ungfdK6hrH+kpGgGnAUZK6I+KhiPhZRr1Tge9HxC0RsRP4DFD5kOlbIuL6iBgBvg3sTQ4j4jsRsSVt7VwB/AfQ1PF+ZpONkzWz9rQ4InrLXt9Iy88HBNyedj99qMH1PlH283DG++kAkrokXZx2Z23j5RaYGVXWewQwqzw5Immxe3W6fBbwaFn9h0eJ80zgmojYnSaS1/Fya9zT6b+jtdbcWrEPb82oMwt4JiKer4itsiWpmlnsuy0Pk4wVfnV29bpU7qduqu/3ukXEJuDjJIn3k5KuljQro+o+xyoiXuTlfV7yeNnPLwIHlhLUtAv17rLz4PXNiN9sMnOyZjaJRMTjEfEHETELOAf4qspmgDbR7wEnAb9B0p01Ny1XKZSK+o8CP69Ijl4ZEcenyx8j6bosObzaB0uaDbwT+P10zNzjJK09x0uaAWxMP++UMW/dy7YAh0p6ZUVsgw38/yMq/u9ukiT4BeCg0oK0dayvrG7lPiyp3E+7gKfGsb6XK0T833TSxRFp/S9kVHuMpIu99Dk9vNyaWZOkI0haVM8DDouIXuA+Xj5vzCyDkzWzSUTS76bJDMCzJBfcPen7J0jGTTXDK0m6IZ8mSRA+X7G88rNuB55PB7D3pC1zr5dUmkhwDbBM0iFp/P+txmd/APgpydi0N6WvXyYZV/a+iAiSWaCflvRBSa9KB/q/XdIljWxk2iX7/4CLJB0o6Q3A2UC9t8C4Cvjvko6UNJ2Xx43tTrfhQEm/Lakb+FOSbsiSJ4C56Xi/cr8v6ah03N6FwLVpd+NY1wckY9YkvVPSNOAlkpbUPRlVrwVOkPSf03GCF1B/svUKknNya/qZHyRpWTOzGpysmbWn72vfe4R9Ny1/K3CbpO0kA7s/lo4Pg+SiekXa/fTecX7+t0i64AaB+4HKLsTLSMY+DUlalSYT7yFJrH5O0hJ0KUmrHCQDzR9Ol/0zyTinas4Evpq2Iu59AV9PlxER1wKnkQzE30KSqHyOZPZoo95H0nK4hWTc3p81cNuUy9Nt+SHJtr1EmohGxHPAH5Hsh0GSlrHy2ZzfSf99WtJdZeXfBr5J0tV4IOmkj3Gsr2QacDHJsXmcZOLIsspKEbEh3YarSVrZtgNPkiTvNUXE/cD/AX5MckzmA/8+2v8z63RKvoSamVnRSfpX4O8jYsKfJlCvtMVwCJgXET/POx6zycgta2Zm1hBJJ0g6SNIrgP8NrCfjNh9m1hxO1szMrFEnkXQLbyG5197p4W4aswnjblAzMzOzAnPLmpmZmVmBOVkzMzMzK7DxPPKkcGbMmBFz587NOwwzMzOzUd15551PRUTfaPUmVbI2d+5c1q5dm3cYZmZmZqOSNNqj9QB3g5qZmZkVWi7JWvrYltsl3ZM+bPqzGXXOkrQ1feDv3ZI+nEesZmZmZnnKqxt0B/DOiNiePsPuFkk3RETlI2tWRMR5OcRnZmZmVgi5JGvpzRO3p2+705dv+GZmZmZWIbcxa5K6JN1N8gDgGyPitoxqp0i6V9K1kuZUWc8SSWslrd26deuExmxmZmbWarklaxExEhFvAmYDR0t6fUWV7wNzI+INwI3AFVXWc0lELIyIhX19o85+NTMzM2sruc8GjYgh4GbguIrypyNiR/r2UuAtrY7NzMzMLG95zQbtk9Sb/twDvBv4SUWdmWVvTwQeaF2EZmZmZsWQ12zQmcAVkrpIEsZrIuIfJF0IrI2I1cBHJZ0I7AaeAc7KKVYzMzOz3CiZmDk5LFy4MPwEAzMzm0xWrRtk+ZqNbBkaZlZvD0sXDbB4QX/eYVkTSLozIhaOVm9SPW7KzMxsMlm1bpBlK9czvGsEgMGhYZatXA/ghK2D5D7BwMzMzLItX7Nxb6JWMrxrhOVrNuYUkeXByZqZmVlBbRkabqjcJicna2ZmZgU1q7enoXKbnJysmZmZFdTSRQP0dHftU9bT3cXSRQM5RWR58AQDMzOzgipNIjj/2nvZObKHfs8G7UhO1szMzAqg2i06Fi/o56rbHwFgxTnH1Kxrk5OTNTMzs5w1cosO386j83jMmpmZWc4auUWHb+fReZysmZmZ5ayRW3T4dh6dx8mamZlZzhq5RYdv59F5nKyZmZnlrJFbdPh2Hp0nl2RN0oGSbpd0j6QNkj6bUWeapBWSNkm6TdLc1kdqZmY28RYv6Oeik+dzQFdyWe7v7eGik+dnThhopK5NDnnNBt0BvDMitkvqBm6RdENE3FpW52zg2Yh4raTTgS8Ap+URrJmZ2UTLukVHM+pa+8ulZS0S29O33ekrKqqdBFyR/nwt8C5JalGIZmZmZoWQ25g1SV2S7gaeBG6MiNsqqvQDjwJExG7gOeCwjPUskbRW0tqtW7dOdNhmZmZmLZVbshYRIxHxJmA2cLSk149xPZdExMKIWNjX19fcIM3MzMxylvts0IgYAm4GjqtYNAjMAZA0FTgYeLq10ZmZmZnlK6/ZoH2SetOfe4B3Az+pqLYaODP9+VTgpoioHNdmZmZmNqnlNRt0JnCFpC6ShPGaiPgHSRcCayNiNXAZ8G1Jm4BngNNzitXMzMwsN7kkaxFxL7Ago/wzZT+/BPxuK+MyMzMzK5rcx6yZmZmZWXVO1szMzMwKzMmamZmZWYE5WTMzMzMrsLxmg5qZmXWsVesGWb5mI1uGhpnV28PSRQNNeRD7RK3X8uVkzczMrIVWrRtk2cr1DO8aAWBwaJhlK9dP6HqdsLU3d4OamZm10PI1G/cmVCXDu0ZYvmZjIddr+XOyZmZm1kJbhoYbKs97vZY/J2tmZmYtNKu3p6HyvNdr+XOyZmZm1kJLFw3Q0921T1lPdxdLFw0Ucr2WPydrZmZmLbR4QT8XnTyfA7qSS3B/bw8XnTx/3JMAJmq9lj/PBjUzM2uxxQv6uer2RwBYcc4xhV+v5SuXljVJcyTdLOl+SRskfSyjzjskPSfp7vT1max1mZmZmU1mebWs7QY+ERF3SXolcKekGyPi/op6P4qI9+QQn5mZmVkh5NKyFhGPRcRd6c/PAw8A7lQ3MzMzq5D7BANJc4EFwG0Zi4+RdI+kGyS9rqWBmZmZmRVArhMMJE0HrgM+HhHbKhbfBRwREdslHQ+sAuZlrGMJsATg8MMPn+CIzczMzFort5Y1Sd0kidqVEbGycnlEbIuI7enP1wPdkmZk1LskIhZGxMK+vr4Jj9vMzMyslfKaDSrgMuCBiPhilTqvSesh6WiSWJ9uXZRmZmZm+curG/RY4APAekl3p2WfAg4HiIivA6cCfyhpNzAMnB4RkUewZmZmZnnJJVmLiFsAjVLny8CXWxORmZmZWTHlPhvUzMzMzKpzsmZmZmZWYE7WzMzMzArMyZqZmZlZgTlZMzMzMyswJ2tmZmZmBeZkzczMzKzAnKyZmZmZFZiTNTMzM7MCy+txU2ZmZpPeqnWDLF+zkS1Dw8zq7WHpogEWL+jvuBhsfJysmZmZTYBV6wZZtnI9w7tGABgcGmbZyvUALUuWihCDjZ+7Qc3MzCbA8jUb9yZJJcO7Rli+ZmNHxWDjl0uyJmmOpJsl3S9pg6SPZdSRpC9J2iTpXklvziNWMzOzsdgyNNxQ+WSNwcYvr5a13cAnIuIo4G3ARyQdVVHnt4B56WsJ8LXWhmhmZjZ2s3p7GiqfrDHY+NWdrEk6tI5Xbz3riojHIuKu9OfngQeAys7zk4BvReJWoFfSzHrjNTMzy9PSRQP0dHftU9bT3cXSRQMdFYONXyMTDLakL9Wo0wUc3kgAkuYCC4DbKhb1A4+Wvd+clj3WyPrNzMzyUBrAf/6197JzZA/9OczELEIMNn6NJGsPRMSCWhUkrWvkwyVNB64DPh4R2xr5v2XrWELSTcrhhzeUJ5qZmU2oxQv6uer2RwBYcc4xHRuDjU8jY9bqOcJ1nwWSukkStSsjYmVGlUFgTtn72WnZPiLikohYGBEL+/r66v14MzMzs7ZQd7IWES8BSDqpcpmkKeV1RiNJwGUkrXVfrFJtNXBGOiv0bcBzEeEuUDMzM+soY7kp7hJJj0fEbZK6gDOBTwK/3MA6jgU+AKyXdHda9inS8W4R8XXgeuB4YBPwIvDBMcRqZmZm1tbGkqz9HvA9Sf8I/CGwHjijkRVExC3UnqhARATwkTHEZ2ZmZjZpjCVZezvwJ8CVwFkR8a9NjcjMzMzM9hpLsnYKMB+YAfytpBuB+9KuSzMzMzNrooaTtYj4EOydJDCPJHGb3+S4zMzMzIyxtawBe8eU/TR9Xde0iMzMzMxsr0YeN3VXM+qYmZmZWf0aaVn7T5LurbFcwMHjjMfMzMzMyjSSrP1KHXVGxhqImZmZme2v7mQtIh7OKpc0LSJ2NC8kMzMzMytp5Nmg1XxV0rubsB4zMzMzqzDuZC0izgYOl/QlSTOaEJOZmZmZpcadrElaBBwJvBb4hqTF447KzMzMzIDmdIPOBS6LiOMj4neAX2/COs3MzMyM5iRrR5O0qpV8erT/IOlySU9Kuq/K8ndIek7S3enrM02I08zMzKztNH3MWkRsq+O/fRM4bpQ6P4qIN6WvC8cbp5mZmVk7ymXMWkT8EHhmvJ9tZmZmNtmN6dmgkt4KHADcB7wGuDwi/jRd9tfAqibEdoyke4AtwP+MiA1NWKeZmVnTrVo3yPI1G9kyNMys3h6WLhpg8YL+vMOqqR1j7lQNJ2uSVgGzgUeBo4B/B75TVmXUMWt1uAs4IiK2SzqeJPmbVyWeJcASgMMPP7wJH21mZla/VesGWbZyPcO7kof4DA4Ns2zleoDCJj/tGHMnG7UbVNLrJF1ZVvR64DzgtIgYAG4BvlZaWOeYtZoiYltEbE9/vh7ornYPt4i4JCIWRsTCvr6+8X60mZlZQ5av2bg36SkZ3jXC8jUbc4podO0Ycyerp2XtB8AxZe8vBs4H3ijpRWA98F8l/SZwd0Q8Od6gJL0GeCIiQtLRJEnl0+Ndr5mZWbNtGRpuqLwI2jHmTlbPBIPfBP6i7P0JwDUR8UvA24GvA13Ae4Eb6vlQSVcBPwYGJG2WdLakcyWdm1Y5FbgvHbP2JeD0iIi6tsjMzKyFZvX2NFReBO0YcycbtWUtItYD7y8r+jDwd5KWkUwwGAD+KSI+XO+HRsT7Rln+ZeDL9a7PzMwsL0sXDewz/gugp7uLpYsGcoyqtnaMuZM1PMEgIrYC75HUTzJ+7bmIuLXpkZmZmbWB0oD886+9l50je+hvg5mV7RhzJxvTrTsAImIQGGxiLGZmZm1p8YJ+rrr9EQBWnHPMKLWLoR1j7lTNeNyUmZmZmU0QJ2tmZmZmBeZkzczMzKzAnKyZmZmZFZiTNTMzM7MCc7JmZmZmVmBO1szMzMwKzMmamZmZWYE5WTMzMzMrMCdrZmZmZgXmZM3MzMyswHJJ1iRdLulJSfdVWS5JX5K0SdK9kt7c6hjNzMzMiiCvlrVvAsfVWP5bwLz0tQT4WgtiMjMzMyucXJK1iPgh8EyNKicB34rErUCvpJmtic7MzMysOIo6Zq0feLTs/ea0bD+SlkhaK2nt1q1bWxKcmZmZWasUNVmrW0RcEhELI2JhX19f3uGYmZmZNdXUvAOoYhCYU/Z+dlpmZmaWq1XrBlm+ZiNbhoaZ1dvD0kUDLF6Q2fnTlib79rWjoiZrq4HzJF0N/CrwXEQ8lnNMZmbW4VatG2TZyvUM7xoBYHBomGUr1+ccVfPU2j4nbPnJ69YdVwE/BgYkbZZ0tqRzJZ2bVrkeeBDYBHwD+KM84jQzMyu3fM3GvYlMyfCuEZav2ZhTRM012bevXeXSshYR7xtleQAfaVE4ZmZmddkyNFy1fPYhPS2OpvlqbZ/lp+0nGJiZmbXKrN7shKxaebuZ7NvXrpysmZmZ1WnpogF6urv2Kevp7mLpooGcImquyb597crJmpmZWZ0WL+jnopPnc0BXcvns7+3hopPnT5rB95N9+9pVUWeDmpmZFdLiBf1cdfsjAKw455ico2m+yb597cgta2ZmZmYF5mTNzMzMrMCcrJmZmZkVmJM1MzMzswJzsmZmZmZWYE7WzMzMzArMyZqZmZlZgTlZMzMzMyuw3JI1ScdJ2ihpk6RPZiw/S9JWSXenrw/nEaeZmZlZnnJ5goGkLuArwLuBzcAdklZHxP0VVVdExHktD9DMzMysIPJ63NTRwKaIeBBA0tXASUBlsmZmZpaLVesGWb5mI1uGhpnV28PSRQMd+4xM74t85dUN2g88WvZ+c1pW6RRJ90q6VtKcrBVJWiJpraS1W7dunYhYzcysw6xaN8iylesZHBomgMGhYZatXM+qdYN5h9Zy3hf5K/IEg+8DcyPiDcCNwBVZlSLikohYGBEL+/r6WhqgmZlNTsvXbGR418g+ZcO7Rli+ZmNOEeXH+yJ/eSVrg0B5S9nstGyviHg6Inakby8F3tKi2MzMrMNtGRpuqHwy877IX17J2h3APElHSjoAOB1YXV5B0syytycCD7QwPjMz62CzensaKp/MvC/yl0uyFhG7gfOANSRJ2DURsUHShZJOTKt9VNIGSfcAHwXOyiNWMzPrPEsXDdDT3bVPWU93F0sXDeQUUX68L/KX12xQIuJ64PqKss+U/bwMWNbquMzMbH+dNhuwtG3nX3svO0f20N8B21xNJ++Lopz3uSVrZmZWPFkXJ4BlK9fvHWRemg249uFnuPknW3O/kE2UxQv6uer2RwBYcc4xOUeTr8m+L4p+3jtZMzPrUJUXqF//lT6uu3Nwv4vTgd1TMmcDXnnrI0T6vlS3pAitEWaVGknKinTeO1kzM+tApXtnlV+gyi9CJcO7Rva7YJVk1b1g9QZ27N6z34UPcMJmuco652slZUU6752smZlNclmtCVn3zqq8CI3F0PCu/crK78lV1Ba3ooxNajdF3m+Vsb24c3dDSVkjqp33F6ze0JT942TNzGwSq9aa0MgFqrene59WAwDRWHJX+blFanGrto+stlr7rYjHtFHNOO+HhnftTeTGs3+K/AQDMzNrwKp1gxx78U0c+cl/5NiLb9rbspDVmtAlZa6jsrSnu4sLTnwdF508nwO6kktGf28P73/b4Zm3czjkoO7M9XZJhb0Lvu/QPzZF3m9ZsVXT29OdeS4347yvNNb945Y1M7NJoNEWtJEIerq79lne093FKW/p55o7NmfeoqFyNuDCIw7d73YOwH6fW/k55QaHhjn24pty7UardYf+2Yf4xq/VFOXJBlldsfXGUErKoPqtScZ63lczlv3jZM3MrA3VOx6nS2Ik9u+4KV1ksi5Q//HEdmD0WzTUup1D5XqXr9mY2RUlXu6iyqsbbVZvT2ZsvkN/bUXYb9W+pPQe1M2zL+4/jqy3p5sXd47UlZRVU+95/+LO3ZkxjGX/OFkzM2szjYzHqdaCVrpITcS9s6qtt7LlIWv8Tx6TEZYuGshsDVy6aGDvdtj+au23iVLvl5RpU6dknvcXnPi6CbtfXOV5X/l7WophLPvHyZpZk1S7f09RZ0pZe6h3Jmc1tVrQWinrLvjVksyJnIxQa/Zi1j5yslZdtf0GTEjXdiNfUp4b3sVfnvamXI9pM5/84GStA1T749RoeSeqNwGD/W+quPQ794Bg10jsLfNNQ62Wem9SW2+iNtEtaI2qjOPYi2/KvOBWm4ww3tsgjDZ7sQj7qN2M1po01kR7vF9SZvX2FOKYNisGJ2uTTL1/7Nc+/EzD5ZP5sTL13tW6WgKWdVPFXXv2HydU6+aJk30fV9NIi2Q7fsEYT8Jf7Sa11cah1RqPU0TVutGqXZCr3QYBsvdnPRf7UrdrkfdTOxlthuhYfxfG8iVlMsktWZN0HPDXQBdwaURcXLF8GvAt4C3A08BpEfFQq+NsJ43ckfyq2x7d7499rfJ2fKzMeC6SjSZgjdyzqtrNE9txHzdivAlxs75gjDexayRhbGT7ss63avdzqjYObSLH40yEat1E1SYjVKr25afaPq72e9rq2YuTWbV9mdW13cjvwmT6kjIWuSRrkrqArwDvBjYDd0haHRH3l1U7G3g2Il4r6XTgC8Bptda7fvA5jr34pkl5oCqN947kWSd8rfJGHq/R6haiibhINuuu1o2YTI/uGc8zJ6slxM34gjHexK7RluqJSvhrjUNrtzFW9U5GqCbry0+1fVztYu9Zn81TbYZoVtd2o78Lk+VLylgoqlycJ/RDpWOACyJiUfp+GUBEXFRWZ01a58eSpgKPA31RI+Bfmn5wfP6Nb2fKFNE3fRpPbd/ByJ5g2tQu5hzaw4zp0wB46OkXAJh72Cv2+f9Z5UWoW1n+1PYdPPjUC+wpO9GnTNE+70cjiaxdWa18PKZMEb84I9menz/1wj7HJKtsRnrs6invPaibrdt37LcvpkjsHtnT1O2oZWrXFPZE7BOH0puOlu/PZsVW2idZ+6iV52zWcQL2Oz+Lot7zu/Q3pPLcGu96m6H0+zRj+jTuf2wbAEfNfNXe5VlljZZPVN1G1vHU9h38bOsLRCTn1khEU36nK/9WNnN/Fnnftyq2Zlyfqin9jSk/L0p/94q4L8rLHu+bwwev+pv9tknSnRGxsOaGk183aD/waNn7zcCvVqsTEbslPQccBjxVXknSEmAJwC++Itk5e/YET2x7aW+dHbtHePCpF/a+f2LbDiKCZ1/YtU8S9+LO/bP5rLJW131q+459Yh6pSAog2ebxXogavUDVa8+e4KGnX9wnmdmxe4SfbU2OSSnm0nF6/qXd+8RQq/yJbfvvrz17gj1NeMphIwnY3MMOAupLRmH8yUxpn2Ttoyefb835XflHuRTDFDXnD3OlZnzBqLfenj2xdz82c721ZJ1vo33xPOiArv3Wk1XWaPlE1W1kHTOmT2P7jt1A7S+pjXz5qfUlp5HYmlG31Z/XqthK+7JyHz/6zDA7dtfXelztd6F0rMrPi4najmau46ADujj0FdMy69crr5a1U4HjIuLD6fsPAL8aEeeV1bkvrbM5ff+ztM5TWesEmDZzXsw886+qfm7Wc756uru46OT5LF7Qz2l/+2Ng32bUrLJq5RNRN+s+LbVUuyN5I2N3Ruv6Ka33wO4pmTf8G69qXRXVyser1nkB9Q2IbbRLcrz7uNq+qLxvVfl2ZHWZZZ2Hq9YNZk7FH+u4olqy9n33FO3TXV3ajlPe0p+5j7LKqz2/b6LOoVpjaerdvmrnW5G7u1utniEPUHsfe3/mI+ta5t+F4resDQJzyt7PTsuy6mxOu0EPJploMGbVBnaXZqmse2SInSN7chv3tmrd4H4xjOV+SvWe4IsX9GcuyypfeMShdf2RbPQht1kaHU+XpZGLZOlRI9X2W7V9NB7j2ce1ZsvVGve2M219KB9nVXm+lT6vvG6jg7WryUokq+37rLLFC/oz91FWebXktxmJXdZ2VEskG90+KPZYxLxV+5sFE/OlypqntO/9uzA2ebWsTQV+CryLJCm7A/i9iNhQVucjwPyIODedYHByRLy31npLLWtjSRgafUZes1vWqt3puJEWtTy+NU5EK1yjLWu1WpMmwx+HL1IrAAAF7UlEQVTwapNJxtuqlbXfmnGcqrVUVmvhnSjjbTmu1ZI3EbNMzazz1NuylkuyBiDpeOCvSG7dcXlE/IWkC4G1EbFa0oHAt4EFwDPA6RHxYK11Tps5LxZ+7G8L1aUE2YlZvV1KtS6Ir5g2tZAXhvF2VTTS3ZVHElAEWcl9M1o1G5X1JacduzDa8f5tZtb+Cp+sTYSFCxfG2rVrgfoThkZar6B2q0FWK1xlYlYtkawVQ7ULYrtdNCb7zU9bLe+xhY12u5uZ2b46PlmrppVdSs0YH+MLojWi3i8p1c7DRgfE+zw0Mxs7J2sNmKgupUZnnk2WFjQrnnpb4SbTWD8zs6JzstYgdylZp3GXsplZvjoyWZO0FXi4Weub0vOqQ7umH9qvrqkHxMjunSPbnxkEmPqqviOQpoy6giBpohu1XuzZvW3rw3uGtz0zSs0ZVNwU2NqGj1178/FrXz527W2yH78jIqJvtEqTKlmb7CStrScDt+LxsWtvPn7ty8euvfn4JUZvHTIzMzOz3DhZMzMzMyswJ2vt5ZK8A7Ax87Frbz5+7cvHrr35+OExa2ZmZmaF5pY1MzMzswJzstYmJH1CUkiakb6XpC9J2iTpXklvzjtG25+k5ZJ+kh6j70rqLVu2LD1+GyUtyjNOyybpuPT4bJL0ybzjsdokzZF0s6T7JW2Q9LG0/FBJN0r6j/TfQ/KO1bJJ6pK0TtI/pO+PlHRb+ju4QtIBeceYBydrbUDSHOA3gUfKin8LmJe+lgBfyyE0G92NwOsj4g3AT4FlAJKOAk4HXgccB3xVUlduUdp+0uPxFZLftaOA96XHzYprN/CJiDgKeBvwkfSYfRL4l4iYB/xL+t6K6WPAA2XvvwD8ZUS8FngWODuXqHLmZK09/CVwPvs+Aesk4FuRuBXolTQzl+isqoj454jYnb69FZid/nwScHVE7IiInwObgKPziNGqOhrYFBEPRsRO4GqS42YFFRGPRcRd6c/Pk1z0+0mO2xVptSuAxflEaLVImg38NnBp+l7AO4Fr0yode+ycrBWcpJOAwYi4p2JRP/Bo2fvNaZkV14eAG9KfffyKz8eojUmaCywAbgNeHRGPpYseB16dU1hW21+RNEzsSd8fBgyVfeHt2N/BqXkHYCDpB8BrMhb9CfApki5QK6haxy8ivpfW+ROSLporWxmbWSeSNB24Dvh4RGxLGmgSERGSfBuEgpH0HuDJiLhT0jvyjqdonKwVQET8Rla5pPnAkcA96R+b2cBdko4GBoE5ZdVnp2XWYtWOX4mks4D3AO+Kl++V4+NXfD5GbUhSN0midmVErEyLn5A0MyIeS4eLPJlfhFbFscCJko4HDgReBfw1yRCfqWnrWsf+DrobtMAiYn1E/EJEzI2IuSRNwG+OiMeB1cAZ6azQtwHPlTXzW0FIOo6kWf/EiHixbNFq4HRJ0yQdSTJR5PY8YrSq7gDmpbPRDiCZELI655ishnSM02XAAxHxxbJFq4Ez05/PBL7X6tistohYFhGz02vd6cBNEfF+4Gbg1LRaxx47t6y1r+uB40kGpr8IfDDfcKyKLwPTgBvT1tFbI+LciNgg6RrgfpLu0Y9ExEiOcVqFiNgt6TxgDdAFXB4RG3IOy2o7FvgAsF7S3WnZp4CLgWsknQ08DLw3p/iscX8MXC3pc8A6kmS84/gJBmZmZmYF5m5QMzMzswJzsmZmZmZWYE7WzMzMzArMyZqZmZlZgTlZMzMzMyswJ2tmZmZmBeZkzczMzKzAnKyZmVUh6e8knSCpV9L1kn4n75jMrPM4WTMzq24+8CzJI24+FxHfzTkeM+tAfoKBmVkGSVOA54Gnga9ExBdyDsnMOpRb1szMss0DtgBnAedK6s43HDPrVE7WzMyyzQdujIibgPuAM3KOx8w6lJM1M7Ns80mSNIDPA8skTc0xHjPrUB6zZmZmZlZgblkzMzMzKzAna2ZmZmYF5mTNzMzMrMCcrJmZmZkVmJM1MzMzswJzsmZmZmZWYE7WzMzMzArMyZqZmZlZgf1/+GEtXXGM2WIAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 720x216 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAnUAAADgCAYAAABlw5CNAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvFvnyVgAAIABJREFUeJzt3XmcXHWd7//Xp6qrk86ekBDSHSBAAiQQQmMFRO74U1yCXAcyjAsooI4j6k8dHTUOmbnXH9cNNHPH0REZGVdAWQ0hDnjjEma8Oga6QgIhgZAQIKQ7Cdk6a6fXz++Pc6qt7q7qrkrX3u/n45FHur7nW+d8qs6pOp/6LueYuyMiIiIilS1S6gBEREREZPiU1ImIiIhUASV1IiIiIlVASZ2IiIhIFVBSJyIiIlIFlNSJiIiIVAEldSJVxsz+zMw2lzqOdMzsTWa2o0jb+rGZfaUY28qGmf3SzD5Qgu1+0Mx+X+ztikjxKakTKRNm9rKZtZnZkZR/38nieW5ms5OP3f3/uvs5BYqxrBKlE1WKRMfd3+HuPynmNkVkZKkpdQAi0sefu/tvSh2EgJlF3b271HFUAjOrcfeuUschMtKppU6kApjZbDP7TzM7aGZ7zez+sPx3YZWnw5a99/bv4gxbAJeY2TNmdtTMfmBm08PuwMNm9hszm5xS/0Ez2xVu63dmdl5YfhPwfuAL4bZ+EZbXm9nPzWyPmb1kZn+Tsq66sHXvgJltAhYO8Tq/ZWavmtkhM1trZn+WsuwWM3vAzO4K495oZvGU5Y1m9lS47H5gdIZtzAX+Fbg0fB2tYfmPzewOM3vMzI4Cbzaz/zCzv055bp8WPjM718x+bWb7zWyzmb1nkNfWu65wPX8ws2+aWauZbTOzN4Tlr5rZa6ldtWFs/xpu63B4LJweLpsVttbWpNtWvxgs3OZr4Xu8wczOD5eNMrN/NLPtZrY73F5dhteSGv8+4Jaw/K/M7Llwf69KiXGw7WZ8beHyN5hZU3g8NpnZG/q9zi+HsRw2s1+Z2dRw2Wgzu8fM9oXvcZOZTQ+XTbTgc7DTzJrN7CtmFs2070QqhZI6kcrwZeBXwGRgJvAvAO7+xnD5Ancf5+73Z3j+XwJvA84G/hz4JfD3wDSC74G/San7S2AOcDLwFPDTcFt3hn9/I9zWn5tZBPgF8DTQALwF+IyZLQrX9f8BZ4X/FgFDjSlrAi4EpgA/Ax40s9Tk7CrgPmASsBL4DoCZ1QIrgLvD5z4YvuYB3P054GPAH8PXMSll8fuArwLjgUG7Z81sLPDrMM6TgWuB75rZvCFeY9IlwDPASeE67iNIemcD1wPfMbNxKfXfT3AcTAXWE+6XHL0deCPBcTAReA+wL1x2W1h+YRhDA/DFIeLfBkwHvmpmVxMcU9cQHFf/F7g3i+1mfG1mNgV4FPg2wfv0T8CjZnZSynPfB3yIYB/UAp8Pyz8QbuvU8LkfA9rCZT8GusLX2RjGNyAJFqk0SupEysuKsFUh+e8jYXkncDpQ7+7H3T3X8WD/4u673b2Z4GT7hLuvc/fjwMMEJzYA3P2H7n7Y3dsJWmAWmNnEDOtdCExz9y+5e4e7bwP+jSDBgeDk/VV33+/urxKcnDNy93vcfZ+7d7n7/wZGAanjA3/v7o+F3aJ3AwvC8tcDMeCf3b3T3R8iSBBz9Yi7/8Hde8L3ZjDvBF529x+F8a4Dfg68O8ttvRQ+txu4nyD5+JK7t7v7r4AOgqQj6VF3/124X/6BoKXx1JxeXXAcjQfOBczdn3P3nWZmwE3A34b76jDwNf60H9Npcfd/CV97G0HSdGu4zq7w+ReGrW5pt5vFa/vvwBZ3vzvczr3A8wQ/TJJ+5O4vhDE8QJCUJl/rScBsd+9297XufihsrbsS+Iy7H3X314BvDvFaRSqCkjqR8rLY3Sel/Pu3sPwLgAFPWtDt+Fc5rnd3yt9taR6Pg2AcmZndZmYvmtkh4OWwztQM6z0dqE9NRAlaa6aHy+uBV1PqvzJYkGb2+bD77mC4ron9tr0r5e9jwOiw27EeaHZ3z3ZbGbw6dJVepwOX9Hvt7wdOyfL5/fcB7p52v/SPzd2PAPsJXnfW3H01Qevm7cBrZnanmU0gaFkbA6xNeS3/JyzPpP97dTrwrZTn7yc4ZhsG2e5Qr62egfvxFYJWxKT+x0TyPbsbWAXcZ2YtZvYNM4uFccaAnSmxfo+gpU+koimpE6kA7r7L3T/i7vXARwm6+WYP9bwT8D7gauCtBAnVrLDckqH0q/8qQYtTaiI63t2vDJfvJGiBSjot04YtGD/3BYLWvclht+jBlG0PZifQELY4DbmtNK8jU/lRgmQnKTVhexX4z36vfZy7fzyLeE9E7/sYdstOAVrCGBkkzj7c/dvu/jpgHkF36BJgL0ESeV7Ka5no7uMyrYf0x8JH+70fde7+X4Nsd6jX1kKQhKU6DWgeJK7k6+x09//l7vOANxC0rN4YxtkOTE2Jc4K7nzfUOkXKnZI6kQpgZu82s5nhwwMEJ9Se8PFu4Mw8bWo8wQlvH0GS8LV+y/tv60ngsJn9nQWTIqJmdr6ZJSdEPAAsNbPJYfyfGmLbXcAeoMbMvghMGKR+qj+Gz/0bM4uZ2TXAxYPU3w3MDMfiDWY9cI2ZjQmT6A+nLPt34GwzuyHcZszMFlowEaMQrjSz/xbG/GVgjbu/6u57CJKc68P3/68IxjAOEMZ3SdhidRQ4DvS4ew9Bt/k3zezksG5DytjIbPwrwb5OTqyZaGbvHmy7Q7024DGC9/h9ZlZjZu8lSAr/fahgzOzNZjY/nABxiKA7tifs9v0V8L/NbIKZRczsLDP7f3J4rSJlSUmdSHn5hfW9Tt3DYflC4AkzO0IwQeDT4fg1CMa9/STsSso4+zJLdxF0bzUDm4A1/Zb/AJgXbmtFOB7snQTjmF4iaPH5PkErH8D/Ctf3EsGJ9O5Btr2KoMvvhfA5x8myO9TdOwgG6H+QoOvuvcDyQZ6yGtgI7DKzvYPU+ybB2LbdwE9ImZwQjjt7O8FYrBaCbsCvE4wDLISfEUw82Q+8jmAyRdJHCFq+9gHnAf+VYR0TCJK3AwTv8T5gWbjs74CtwJqw6/039B3POCh3f5jg9d8XPv9Z4B1ZbDfja3P3fQTH1+fC53wBeKe7D7bPkk4BHiJI6J4D/pM/HX83Ekyq2BTG9BAwI9vXKlKurO8QFBERKTdm9mNgh7v/j1LHkm/V/NpEik0tdSIiIiJVoKyTOjP7oQUXq3w2w3Izs2+b2VYLLqx6UbFjFBERESkHZd39amZvBI4Ad7n7+WmWX0kw8PpKggthfsvdLylulCIiIiKlV9Ytde7+O4KBs5lcTZDwubuvASaZmQa7ioiIyIhT1kldFhroOztuB30vStnLzG4ys0T476aiRCciIiJSJDVDV6kOHty38k6AqVOnejwe/16JQxIREREZ0tq1a/e6+2B3eAEqP6lrpu/V6meSxZXGZ82aRSKRKFhQIiIiIvliZlnd9rDSu19XAjeGs2BfDxzsd5NoERERkRGhrFvqzOxe4E3AVDPbQXDF8RiAu/8rwS1kriS4Cvox4EOliVRERESktMo6qXP364ZY7sAnihSOiIiISNmq9O5XEREREUFJnYiIiEhVUFInIiIiUgWU1ImIiIhUgbKeKCEiItVrxbpmlq3aTEtrG/WT6liy6ByAAWWLG9PeKEhE+lFSJyIiRbdiXTNLl2+grbMbgObWNpY8+DQYdHZ7b9nS5RsAlNiJZEHdryIiUnTLVm3uTeiSOnu8N6FLauvsZtmqzcUMTaRiqaVORESKrqW1jR6OsXvUUrqtdcDyWM8sTu64BcNoaW0rQYQilUdJnYiIFF39pDq2HlpPR+RF6rovIeoTe5d1WgvHo2vpZh81TKV+Ul0JIxWpHErqRESk6JYsOoebHr4HgJM6/pYo44hFDAyO9DzHrujn6YhsYXx0eu8EChEZnMbUiYhI0S1ubOCsmbuI+QyijKNhUh3L3r2AZe9awNjIWeBRase8xK3XzNckCZEsqaVORERKouXYRk6pO49LZkzh/o9e2lt+75Pb2bf3TM47fa8SOpEcqKVORESK7rWjr7H94HamxM5Nu3xy7FwSLQncPe1yERlISZ2IiBRdoiUBwJTY3LTLp8Tmsr9tPy+1vlTMsEQqmpI6EREpuqbmJgxjUs3ZaZcnk71k8iciQ1NSJyIiRZfYmeDcqecSi4xNu3xCzZnURmtpam4qcmQilUtJnYiIFJW7k2hJsLBhYcY6UYtx4SkXktipljqRbCmpExGRomo+3MyuI7uIz4gPWi8+I87alrX0eE+RIhOpbGWd1JnZFWa22cy2mtnNaZafZmaPm9k6M3vGzK4sRZwiIpK95Di5eP0QSV19nMMdh3lh3wvFCEuk4pVtUmdmUeB24B3APOA6M5vXr9r/AB5w90bgWuC7xY1SRERy1dTcRNSiXHjKhYPWS3bParKESHbKNqkDLga2uvs2d+8A7gOu7lfHgQnh3xOBliLGJyIiJyCxM8H5J59PXWzwe7qeO/VcxsTGaLKESJbKOalrAF5NebwjLEt1C3C9me0AHgM+lWllZnaTmSXMLLFnz558xyoiIlnonSRRn3mSRFJNpIaLZlykyRIiWSrnpC4b1wE/dveZwJXA3WaW9jW5+53uHnf3+LRp04oapIiIBF5qfYn9bfuHHE+XFJ8RZ93OdXT1dBU4MpHKV85JXTNwasrjmWFZqg8DDwC4+x+B0cDUokQnIiI5y3aSRFK8Pk5bVxub9mwqZFgiVaGck7omYI6ZnWFmtQQTIVb2q7MdeAuAmc0lSOrUtyoiUqaampuojdYyf/r8rOprsoRI9so2qXP3LuCTwCrgOYJZrhvN7EtmdlVY7XPAR8zsaeBe4IOuuz+LiJStxM4EC6YvoDZam1X92VNmM2HUBE2WEMlCTakDGIy7P0YwASK17Ispf28CLit2XCIikrse72Fty1puuOCGrJ8TsQjx+rgmS4hkoWxb6kREpLq8sO8FDncczno8XVJ8Rpyndz1Ne1d7gSITqQ5K6kREpChynSSRFK+P09nTyYbXNhQiLJGqoaRORESKoqm5iTGxMcydNjen52myhEh2lNSJiEhRJHYmaDylkZpIbsO5T594OifVnaSkTmQISupERKSgVqxr5g23/po/bl/LtpZTWLGu/yVHB2dmnDZ+Pveuf5wzbn6Uy25bnfM6REaCsp79KiIilW3FumaWLt/Awa6t+Oh2OtpmsXR5MDZucWP/Oz9mXsf2XTM4Zr+jm+M0t5LzOkRGArXUiYhIwSxbtZm2zm7aI1sAGNUzh7bObpat2pzTOiKdZ4H10BF5CSDndYiMBErqRESkYFpa2wDoiGzBfAw1Xt+nPNt11PbMCdZjLwxYt4gElNSJiEjB1E+qA6AjspXantlYeNpJlme7jhpOIupT6IhsHbBuEQkoqRMRkYJZsugcRsW66bCXGBW2ttXFoixZdE5O66iLRantmdPbjZvrOkRGAiV1IiJSMIsbG/jI5TVgXdT2zKFhUh23XjM/pwkOixsbuPWa+dT5HLqsmekTe3Jeh8hIoKROREQKatz47QBc3BDnDzdffkLJ2OLGBs6b1gjmfOO6CUroRNJQUiciIgXV1NxErU1kbHTGsNYzOXZusL6WpnyEJVJ1lNSJiEhBJXYmmBI7FzMb1npGRyYzJnKK7iwhkoGSOhERKZhjncfY+NpGJsdyu99rJlNi56qlTiQDJXUiIlIw63etp9u7mRJ2nQ7X5Nhcth3Yxv62/XlZn0g1UVInIiIFk+wqnZK3lrpgPWtb1uZlfSLVREmdiIgUTFNLEzPGzaAuOi0v65scO6d3vSLSV1kndWZ2hZltNrOtZnZzhjrvMbNNZrbRzH5W7BhFRCSzREuCeH08b+urjYxnzpQ5miwhkkbZJnVmFgVuB94BzAOuM7N5/erMAZYCl7n7ecBnih6oiIikdaj9EJv3bmZh/cK8rjdeH1dLnUgaZZvUARcDW919m7t3APcBV/er8xHgdnc/AODurxU5RhERyeCpnU/heF5b6gAW1i9kx6Ed7DqyK6/rFal05ZzUNQCvpjzeEZalOhs428z+YGZrzOyKokUnIiKDSnaR5jupS65PkyVE+irnpC4bNcAc4E3AdcC/mdmkdBXN7CYzS5hZYs+ePUUMUURkZGpqaeL0iaczbWx+JkkkNc5oJGIRdcGK9FPOSV0zcGrK45lhWaodwEp373T3l4AXCJK8Adz9TnePu3t82rT8fsGIiMhA+Z4kkTSudhxzp87VZAmRfso5qWsC5pjZGWZWC1wLrOxXZwVBKx1mNpWgO3ZbMYMUEZGB9rftZ9uBbXmfJJEUr4+TaEng7gVZv0glKtukzt27gE8Cq4DngAfcfaOZfcnMrgqrrQL2mdkm4HFgibvvK03EIiKSVKjxdEkL6xey++hudhzaUZD1i1SimlIHMBh3fwx4rF/ZF1P+duCz4T8RESkTyaTudfWvK8j6k8lioiXBqRNPHaK2yMiQU0udmU3J4l/aiQoiIjJyNLU0MWfKHCaNLswpYcEpC6iJ1GiyhEiKXFvqWsJ/NkidKHDaCUckIiIVL9GS4I2nv7Fg6x9dM5r5J8/XZAmRFLkmdc+5e+NgFcxs3TDiERGRCrfryC52HNpBfEZhxtMlxevjPLTpIdwds8HaGkRGhlwnSlyapzoiIlKlCj1JImlh/UIOHD/AtgO66IEI5JjUuftxADP7Sv9l4b1ae+uIiMjIlGhJELEIjTMG7dgZttTJEiJy4pc0aTCz65IPzOxk4Df5CUlERCrRinXNXHbbar6x+jFGcRq/2XiwoNt7sWUyRoyb7nuAy25bzYp1/a9PLzKynOglTT4KrDKzFwEHfgT8Xd6iEhGRirJiXTNLl2/gWGcX7aO3UNcZZ+nyDQAsbux/2+78bO9/rnieWORM2iMv0NzaVtDtiVSCXC9pcpeZfYZg3NwngDuBO4DF7v5oAeITEZEKsGzVZto6u+m2vfRYK6N6ZtPW2c2yVZsLur1RPbPpiLyI01PQ7YlUgly7X39McDmTDwH3ALOAA8D1ZvauvEYmIiIVo6W1DYAO2wJAbc+cPuWF2l5tzxzc2uiy5oJuT6QS5NT96u6rgdXJx2ZWA8wFFgCXAA/lNToREakI9ZPqaG5toz2yBTxKrZ/RW17I7SWTx/bIFmLdpxZseyKVYFj3fnX3Lnff4O73uPuSfAUlIiKVZcmic6iLRemIbKXWZ2HUUheLsmTROQXdXsxnYj6ajsiWgm5PpBLkOqbuqXzUERGR6rK4sYGv/cX5dES2UNszh4ZJddx6zfyCTVpY3NjArdfMZ1Q0Rm3PWRB7saDbE6kEuc5+nWtmzwyy3ICJw4hHREQq1AWz2umxI5wzZQF/+NvLC769xY0N3PvkduzQ+WzveIR3Lphe8G2KlLNck7pzs6jTfSKBiIhIZWtqaQJgSiybU0X+TInN5YVj97PxtY0sOGVBUbctUk5ynSjxCoCZXQ68H2gFngWeAZ519/a8RygiIhUh0ZIgQi0Ta84q6nYnx+b2bl9JnYxkJzpR4ofAL4A1wJnAF4GN+QpKREQqT1NLE5Nic4jYiV7X/sSMizYwcdTE3pZCkZHqRD95r7j7ivDvB/MVjIiIVKbunm6e2vkUp8SuKPq2zSLE6+O6B6yMeCd6R4k1ZvbZAsUkIiIV5oV9L3Ck40jRx9MlxevjPLP7Gdq7NApIRq4TvaPEdOAGM3vFzFaa2ZfN7N35Ds7MrjCzzWa21cxuHqTeX5qZm1k83zGIiMjQkl2fk2tKk9QtrF9IZ08nz+we7AINItUtn3eUWEgeu2LNLArcDrwN2AE0mdlKd9/Ur9544NPAE/natoiI5CbRkmBsbCzja04vyfbj9fHeOBY2LCxJDCKlls87SnwhX0GFLga2uvs2d+8A7gOuTlPvy8DXgeN53r6IiGSpqaWJi2ZcRMSiJdn+aRNPY+qYqZosISPasJK6AmsAXk15vCMs62VmFwGnuvujxQxMRET+pLO7k/W71ve2lpWCmbGwfqEmS8iIVs5J3aDMLAL8E/C5LOvfZGYJM0vs2bOnsMGJiIwgm/Zs4njXcRbWl7bbM14fZ+OejRzrPFbSOERKpZyTumbg1JTHM8OypPHA+cB/mNnLwOuBlZkmS7j7ne4ed/f4tGnTChSyiMjIk+zyLGVLHQSTJXq8h3U715U0DpFSyUtSZ2YzzGxUPtaVogmYY2ZnmFktcC2wMrnQ3Q+6+1R3n+XuswguhHyVu6vtXUSkiBItCSaOmsjsKbNLGkfqZAmRkShfLXV3A8+b2T/maX24exfwSWAV8BzwgLtvNLMvmdlV+dqOiIgMT1NLE/H6OGZW0jhmjJ9Bw/gGTZaQESunS5qY2ZQMi94DRAla1ia5e+uwIwPc/THgsX5lX8xQ90352KaIiGTveNdxNuzewGcvLY/r0evOEjKS5XqbsJbwX/+fYx7+bwTJ3WnDjEtERCrAM7ufobOns+STJJLi9XEe2fwIB48fZOLoiaUOR6Sock3qnnP3xsEqmJlGqIqIjBDJVrFST5JISiaXT+18ijef8eYSRyNSXLmOqbs0T3VERKQKJFoSTB0zldMmlkcHzevqXwdosoSMTDklde5+HMDMvtJ/WXhbr946IiJS/ZpamlhYv7DkkySSpo6ZyhmTztBkCRmRTnT2a4OZXZd8YGYnA7/JT0giIlIJjnYcZdOeTWXT9ZqkyRIyUp1oUvdR4CYzu9jMFgKrgbxdzkRERMrbinXNvP4fv0+P9/Dgf9WyYl3z0E8qkjqfw0utL3HazT/jsttWl1VsIoWU6yVN7gKeAtYBnwB+BnQBi919a/7DExGRcrNiXTNLl29gd89GqIXDh09j6fINpQ4LCGJ7/NkJEIX2yFaaWyf2xra4sWGIZ4tUtlxb6n5McNmSDwH3ALOAA8D1ZvauvEYmIiJladmqzbR1dtMR2ULUp1DDSbR1drNs1eZShxbE0HEGAB2RoK2hXGITKbScWurcfTVBVysAZlYDzAUWAJcAD+U1OhERKTstrW0AtEe2UNszp0/5zMl1pQqrN4YIY6npmUl7ZEufcpFqN6zbhLl7l7tvcPd73H1JvoISEZHyVT+pjh6O0hVp7pPU1U8qbUKXGsOontl0pCR15RCbSKHllNSZ2VP5qCMiIpVryaJzoHYbAKPCpK4uFg3KS2zJonOoi0Wp7ZlDt+2ji/1lE5tIoeV6R4m5ZvbMIMsN0H1ZRESq2OLGBh7ecpi7noPantk0TKpjyaJzWNzYwL1Pbi95bACf+vmzHADGj3+FW6+8XJMkZETINak7N4s63ScSiIiIVI4228LY6AzecMYZ3P/R8rqR0OLGBu5ecxHNr0V49xs6lNDJiJHrRIlXAMzscuD9QCvwLPAM8Ky7t+c9QhERKTtNLU1Mjs0tdRgZ1UTqmFBzhi5CLCPKoGPqzOw8M/tpmkU/BH4BrAHOBL4IbMx/eCIiUm72HtvLy60vM6Umm86b0pkcO5dESwJ3L3UoIkUxVEvdb4B07eqvuPuK8O8H8xuSiIiUs7UtawGYUsYtdQBTYufy8qFH2X5wO6dPOr3U4YgU3FCzX98OfDX5wMzuMrPPAGvM7LMFjUxERMpSU0sTAJNi5T2jNJl0qgtWRopBk7rwGnTvTyn6McEM1+nADWb2ipmtNLMvm9m7CxiniIiUiURLgrNPOpvayLhShzKoiTWziUVivUmoSLXL5x0lFqKuWBGRqtfU0sSbZ72Zzj2ljmRwUavlgukXqKVORox83lHiC/kKKsnMrjCzzWa21cxuTrP8s2a2ycyeMbPfmpkGTYiIFFDL4RZaDrcQr4+XOpSsxOvjmiwhI8awkrpCMrMocDvwDmAecJ2ZzetXbR0Qd/cLCO47+43iRikiMrIkJ0ksrF9Y4kiyE6+Pc7D9IFv3by11KCIFV7ZJHXAxsNXdt7l7B3AfcHVqBXd/3N2PhQ/XADOLHKOIyIjS1NJExCJceMqFpQ4lK8nkU12wMhKUc1LXALya8nhHWJbJh4FfZlpoZjeZWcLMEnv2lPlAEBGRMpVoSXDetPMYWzu21KFkZd60eYyuGa3JEjIilHNSlzUzux6IA8sy1XH3O9097u7xadOmFS84EZEq4e4kWhIVM54OIBaNceEpF6qlTkaEck7qmoFTUx7PDMv6MLO3Av8AXKXblImIFM72g9vZc2xPRSV1EHTBPrXzKbp7dGtyqW7lnNQ1AXPM7AwzqwWuBVamVjCzRuB7BAndayWIUURkxEi2dlXKJImkeH2co51HeX7v86UORaSgyjapc/cu4JPAKuA54AF332hmXzKzq8Jqy4BxwINmtt7MVmZYnYiIDFNTSxOxSIwLpl9Q6lByoskSMlLkdPHhYnP3x4DH+pV9MeXvtxY9KBGRESrRkuCC6RcwqmZUqUPJydknnc242nE0tTTxgQs/UOpwRAqmbFvqRESkfFTiJImkaCTKRTMuUkudVD0ldSIiMqSt+7dysP1gxY2nS1pYv5D1u9bT2d1Z6lBECkZJnYiIDCnZylWJLXUQxN3e3c6zrz1b6lBECkZJnYiIDCnRkmB0zWjmTet/t8bKoMkSMhIoqRMRkYxWrGvmsttWc/sfVhHrPpNHn6nMq0edOflMxsYm8j8fe4Qzbn6Uy25bzYp1Ay59KlLRlNSJiEhaK9Y1s3T5Bna0HqEj8iJ0nsXS5RsqMhl6ZH0L3n4mB7qex4Hm1raKfS0imSipExGRtJat2kxbZzed1ozbcWp75tDW2c2yVZtLHVrOlq3aTE33bDrsZZwOgIp9LSKZKKkTEZG0WlrbAOiIBInPqJ45fcorSUtrWxC/ddMeebFPuUi1KOuLD4uIFNqKdc0sW7WZltY26ifVsWTROQADyhY3NpQ40uKrn1RHc2sbR6KPE+05mRpv6C2vNPWT6tjeOh88xtHo44zumdtbPtKkO+YXNzZkLJfKoaROREaETMnb0uUbaOsMbvTe3NrGkgefBoPObu8tW7p8Q+96RtJJb8mic/js8sdojz7DpM4bMSLUxaK9710lWbLoHJYu72Bs959xNPo4kzs/yNjY+IrLhE0LAAAQh0lEQVR8Lbnof9y/+dxp/Hxtc59jfunyDSRe2Z+2PGkkHfeVTEmdiFS95ID//ies0bFIb1lSZ48PeH5bZze3rNxIe1dP2pNetZ7gFjc28INnn+ClF2sY1/V2Gir4hJ6M+VM/fydHbTWx8f/FrVd+viJfS7bSHfc/XbOd/kd4W2c39z7xKt3uA8pH4nFfyZTUiUjV6d86cayja0Dy1tbZPaBsMK1tA+9EkDrQvhpbMo51HuP3LT/n1Lo3cemMM7n/o5eWOqRhWdzYwM+euIRf75vD+JMe5+oL/6nUIeVNupbo5ESXVAN/sgT6J3RJmY77W1ZurMpjvtJpooSIVJVk60Rza1vvpSsOHCvcraGSLRep26uWS2Xc/+z9tB5vZfaYvyh1KHljZpw15i94evfTrNmxptTh5EW6Yz75OFtRs5y22drWWZXHfKVTS52IVKxsWycymVQX69O1BBCLWJ8xdQB1sSijY5G0yWHULG0r4LJVmyu+5eKOxB3MmzaPqZELSx1KXp02+u28cPwO7kjcwaWnVnbrI5D2mG/r7CZqlrYFzujbYlcXi/KXr2voM6YuWZ7puO+v2lutK4WSOhGpSJnGyWWb0NXFotxy1XkAfOGhZ+jo7ukdM5aprP/662LRjNtrbm3jsttWV+zJbW3LWppamvj2Fd/m9+tza8Upd7HIGG644AZ+sO4HfHPRNzlpzEmlDikn/X/MZGqR63YfcIwmE7gHmnb0Ob4XNzYQP31KVsd9Jv0/gxp/V3xK6kSk7GXbIjdY68SkuhjHOroHnMgA7n1yO0CfMWPpymBgsrds1ea0J1WD3vJKPLndkbiDMbEx3LjgRn6/flOpw8m7j8U/xncT3+VH63/E59/w+VKHk7V0P2b6t7wlJY/R/sfs4sYGtuw+AvQ9vhc3NmR13B/r6Mqp1Vrj74pnRCZ1G5oPctltq3VgZSmX63jpml+Sb7m2yGVqnbjlqvMynrCylemk1z+edCfZSuqeaj3eyr3P3sv7zn8fE0dPLHU4BTF/+nwuO/Uyvrf2e3z20s8SsfIbYp7L5Id0XarJYyvfx33/z2Rye5k+k61tnb0TLirxB04lGZFJHfS9Ns/jz+8p6y/YfMvlwpOQ/XW80l3naCRf80vJ8InJduZqpha5wVonkiemfEruk9TtZeoOy9Q9VW7fQ3c/fTfHOo/x8YUfL1kMxfDx+Me5/uHr+e223/K2s95W0liyvZ5cpsTJgdpoJG1LdL6lO+YHa7Xur5J+4ORboS/8XNZJnZldAXwLiALfd/fb+i0fBdwFvA7YB7zX3V/Odv1tnd19rtlTjQnHcC88mct1vNJd52gkXPOr2pPhXH8EDKduuvctk0wtcvlqnchF/+1ddtvqtLFn6p7K9D2Uy37N1376xv95nsTxZYyLnsv2XdO5aMYJvCEV4l3z3sX/++jfcO3PvsSEox15Pb5z3XfZXk9usB8zMycHd8coxTGflI/xd5D+u6EQ30PF/O48kQs/5xpf2SZ1ZhYFbgfeBuwAmsxspbunDu74MHDA3Web2bXA14H35rKddB+ackk4hvslDQNPkLleeDKX63hlus5ROpV6za9sP6SlSIYh/1+EMPAYyjURzaVuuvctk8Fa5EotuHtB9t1TuXTVZirL13460L2ezlGvMqH9032OrWr0yw17qTl+OfttOXXspbl1al6O78FaX/NxPblMP2YK0RKdi3QteCcy/q7/d1w+vlvysZ+Kff49kRn0ZZvUARcDW919G4CZ3QdcDaQmdVcDt4R/PwR8x8zMPYfsIo1cL7aYj18Kw8noczlB5nrhyVxk+gWZi3Lunsrl13Sxk+FCfRGmO4ZyTURzqZvLzNVStMhlcrj9MPs6gq+mJ5uj1J8MH3pzF99ZvZXO7h6mjR/N9a8/jXvWbGfP4eNZrXPbQfjM8g20d3WDBY8/9VCwn7q6vbfsM8s3UFtjHO7qCgZWhdq74EdPPkePe99y50+DsNLUPRR7mIiPZUz3n9HWHZxYki1A1WbZqs2M7lgEox/iUOznjO160wm/b/3Lf/jEpuC7IWU/PfLcNB5/bk+ffZq6j7ORPJb6H1v1Jzf3HoOd3XFi0dgw350Tk4/xd/3l47slU+KUrpV8sASwmOffliy6svuzYeY/BWNm7wKucPe/Dh/fAFzi7p9MqfNsWGdH+PjFsM7ewdZ91riJ/rUF/23YMUYixrRxo9hzpJ2elAPJwos4pr63mepmKs/EzCjEPsu03ppohB73E359mepGzOjq7sk6jv4iEePMqWOZOm4UL+87CsCsk8b2Lk9Xlqk8XdneI+28tPco3T3OqJoop06p49X9bcGXb54Vap+Wy/byoSYaobvHcf/T/pg6bhQAm3YeAmDejAl9npOuvFB1nxy9m/fM/NUJvrryc+XLs7jxhXm9jyfUBQlCvt+3fKxjONtbs20fALc1NrF+2h6qxdpt7+akntG9j0u9n/YeaefFPUf7fH4L9X2aD5m+I4tx/t02sYHvXXA1EPRG/OHmy5N11rp7fKh1lXNLXV6Z2U3ATQBnjp3AqJook8bE0iZZmRKO/np6nNcOtw/Yyel2eqa6mcozKcQBNVjCOeukMQADEpx0ZVPHjWL86Jqs6gJs23t0wPaySWwheN9e3R/8itl9KHj/Dhzt7I3jWEf6L4v+5XuPtA94fv/Y2ru6B8Q6lGInw7nIxzGUy5dbLnXTvW/J4/BIexcwMFEfUxtNu6505YWqe3bHJG57/g0ATJ/wp5Pp7kPHB5QdbOtk96Hj9LgTi0YYO6qGg22dfd6jfJ08TmQ/mRtzD0zpLR9VEy3Y+5aPdQxne6NqorR3dfOJZxewdWJrb3mhju9c1mNmTKyLcaits/dYmTpuFBPDBDvdsZUsG1fTt5Wu1Ptp6rhRaT+/6c4BuXzHFWo/ZapXzPNvsjciV+XcUncpcIu7LwofLwVw91tT6qwK6/zRzGqAXcC0obpf4/G4JxIJILtB7uUkly7OdFfLT154MpcxBIWSaWxJLre2yeXCmivWNfcZ69G/izv5/MHuHJDt1dlvvWY+UPgxbkPFnE6+jqH+712mOzHkUjfT+1YO4+QKabifhULtp+Q+qdb3P1PX4HDft0zXjRtqxraO+/Tfcfn4bsnHfirl+TfblrpyTupqgBeAtwDNQBPwPnffmFLnE8B8d/9YOFHiGnd/z1DrTk3qMkl3WYVcTvTp5HqgZHsrl2o6Qab7ks31g5ft+5ZpvYPJlEQWc7xfMb8IBzuGdMmWwkr3WSjFfqr2fVKI2d2ZfjCOxMT5RBTqu2W4+6mU59+KT+oAzOxK4J8JLmnyQ3f/qpl9CUi4+0ozGw3cDTQC+4FrkxMrBpNNUtdfoX7RDVWebUYP1XOCzPaDl0tLaj4mcZT7r2klWdVH+6lyFfp6ZJIfue6nUn0mqyKpK5QTSeqgcNfr0gd9aPnoqs1WpmZz/ZoWEZFSUFI3iBNN6qS8FKqrthK7rUVEpHpp9qtUvWSSNdwxEpnGwymJExGRSqKkTira4saGAclX/PQpaVvZMpWLiIhUA3W/ioiIiJSxbLtfI8UIRkREREQKS0mdiIiISBVQUiciIiJSBZTUiYiIiFQBJXUiIiIiVUBJnYiIiEgVUFInIiIiUgWU1ImIiIhUASV1IiIiIlVASZ2IiIhIFVBSJyIiIlIFlNSJiIiIVAEldSIiIiJVQEmdiIiISBUoy6TOzKaY2a/NbEv4/+Q0dS40sz+a2UYze8bM3luKWEVERETKQVkmdcDNwG/dfQ7w2/Bxf8eAG939POAK4J/NbFIRYxQREREpG+Wa1F0N/CT8+yfA4v4V3P0Fd98S/t0CvAZMK1qEIiIiImWkXJO66e6+M/x7FzB9sMpmdjFQC7xY6MBEREREylFNqTZsZr8BTkmz6B9SH7i7m5kPsp4ZwN3AB9y9Z5B6NwE3hQ+PmNnm3KOuCFOBvaUOQk6Y9l9l0/6rXNp3la3a99/p2VQy94z5UsmECdeb3H1nmLT9h7ufk6beBOA/gK+5+0NFDrMsmVnC3eOljkNOjPZfZdP+q1zad5VN+y9Qrt2vK4EPhH9/AHikfwUzqwUeBu5SQiciIiIjXbkmdbcBbzOzLcBbw8eYWdzMvh/WeQ/wRuCDZrY+/HdhacIVERERKa2SjakbjLvvA96SpjwB/HX49z3APUUOrRLcWeoAZFi0/yqb9l/l0r6rbNp/lOmYOhERERHJTbl2v4qIiIhIDpTUVRkz+5yZuZlNDR+bmX3bzLaGt1O7qNQxykBmtszMng/30cOpd0cxs6Xh/ttsZotKGaekZ2ZXhPtnq5mluwOOlBEzO9XMHjezTeGtJj8dlg95i0opD2YWNbN1Zvbv4eMzzOyJ8DN4fziZcsRRUldFzOxU4O3A9pTidwBzwn83AXeUIDQZ2q+B8939AuAFYCmAmc0DrgWSt8P7rplFSxalDBDuj9sJPmvzgOvC/Sblqwv4nLvPA14PfCLcZ9ncolLKw6eB51Iefx34prvPBg4AHy5JVCWmpK66fBP4ApA6UPJqgsu+uLuvASaF1/6TMuLuv3L3rvDhGmBm+PfVwH3u3u7uLwFbgYtLEaNkdDGw1d23uXsHcB/BfpMy5e473f2p8O/DBMlBA1ncolJKz8xmAv8d+H742IDLgeTlzUbsvlNSVyXM7Gqg2d2f7reoAXg15fGOsEzK118Bvwz/1v4rf9pHFczMZgGNwBPkeItKKZl/JmjASN5F6iSgNeWH8Yj9DJblJU0kvSFurfb3BF2vUqYG23/u/khY5x8IuoZ+WszYREYiMxsH/Bz4jLsfChp8AkPdolJKw8zeCbzm7mvN7E2ljqfcKKmrIO7+1nTlZjYfOAN4OvxSmgk8ZWYXA83AqSnVZ4ZlUmSZ9l+SmX0QeCfwFv/TtYa0/8qf9lEFMrMYQUL3U3dfHhbvNrMZKbeofK10EUoGlwFXmdmVwGhgAvAtgqFFNWFr3Yj9DKr7tQq4+wZ3P9ndZ7n7LIKm54vcfRfBLdduDGfBvh44mNK9IGXCzK4g6E64yt2PpSxaCVxrZqPM7AyCCS9PliJGyagJmBPOvqslmNiyssQxySDCMVg/AJ5z939KWTTkLSqltNx9qbvPDM911wKr3f39wOPAu8JqI3bfqaWu+j0GXEkwwP4Y8KHShiMZfAcYBfw6bG1d4+4fc/eNZvYAsImgW/YT7t5dwjilH3fvMrNPAquAKPBDd99Y4rBkcJcBNwAbzGx9WPb3BLekfMDMPgy8QnA7SqkMfwfcZ2ZfAdYRJO0jju4oISIiIlIF1P0qIiIiUgWU1ImIiIhUASV1IiIiIlVASZ2IiIhIFVBSJyIiIlIFlNSJiIiIVAEldSIiIiJVQEmdiMgwmdmPzOzPzWySmT1mZn9R6phEZORRUiciMnzzgQMEtyb6irs/XOJ4RGQE0h0lRESGwcwiwGFgH3C7u3+9xCGJyAilljoRkeGZA7QAHwQ+Zmax0oYjIiOVkjoRkeGZD/za3VcDzwI3ljgeERmhlNSJiAzPfIJkDuBrwFIzqylhPCIyQmlMnYiIiEgVUEudiIiISBVQUiciIiJSBZTUiYiIiFQBJXUiIiIiVUBJnYiIiEgVUFInIiIiUgWU1ImIiIhUASV1IiIiIlXg/we/qLsoUamivAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 720x216 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import scipy.signal as sig\n",
    "\n",
    "\n",
    "N = 10000  # number of samples for input signal\n",
    "K = 50  # limit for lags in ACF\n",
    "\n",
    "# generate input signal\n",
    "np.random.seed(5)\n",
    "x = np.random.normal(size=N)  # normally distributed (zero-mean, unit-variance) white noise\n",
    "# impulse response of the system\n",
    "h = np.concatenate((np.zeros(10), sig.triang(10), np.zeros(10)))\n",
    "# output signal by convolution\n",
    "y = np.convolve(h, x, mode='full')\n",
    "\n",
    "# compute correlation functions\n",
    "acfx = 1/len(x) * np.correlate(x, x, mode='full')\n",
    "acfy = 1/len(y) * np.correlate(y, y, mode='full')\n",
    "ccfyx = 1/len(y) * np.correlate(y, x, mode='full')\n",
    "\n",
    "def plot_correlation_function(cf):\n",
    "    cf = cf[N-K-1:N+K-1]\n",
    "    kappa = np.arange(-len(cf)//2,len(cf)//2)\n",
    "    plt.stem(kappa, cf)\n",
    "    plt.xlabel(r'$\\kappa$')\n",
    "    plt.axis([-K, K, -0.2, 1.1*max(cf)])\n",
    "\n",
    "# plot ACFs and CCF\n",
    "plt.rc('figure', figsize=(10, 3))\n",
    "plt.figure()\n",
    "plot_correlation_function(acfx)\n",
    "plt.title('Estimated ACF of input signal')\n",
    "plt.ylabel(r'$\\hat{\\varphi}_{xx}[\\kappa]$')\n",
    "\n",
    "plt.figure()\n",
    "plot_correlation_function(acfy)\n",
    "plt.title('Estimated ACF of output signal')\n",
    "plt.ylabel(r'$\\hat{\\varphi}_{yy}[\\kappa]$')\n",
    "\n",
    "plt.figure()\n",
    "plot_correlation_function(ccfyx)\n",
    "plt.plot(np.arange(len(h)), h, 'g-')\n",
    "plt.title('Estimated and true impulse response')\n",
    "plt.ylabel(r'$\\hat{h}[k]$, $h[k]$');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise**\n",
    "\n",
    "* Why is the estimated CCF $\\hat{\\varphi}_{yx}[k]$ not exactly equal to the true impulse response $h[k]$ of the system?\n",
    "* What changes if you change the number of samples `N` of the input signal?\n",
    "\n",
    "Solution: The derived relations for system identification hold for the case of a wide-sense ergodic input signal of infinite duration. Since we can only numerically simulate signals of finite duration, the observed deviations are a result of the resulting statistical uncertainties. Increasing the length `N` of the input signal improves the estimate of the impulse response."
   ]
  },
  {
   "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.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}