{
 "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": [
    "## Linear Mean\n",
    "\n",
    "In the following we aim at finding a relation between the linear mean $\\mu_x[k]$ of the input signal $x[k]$ and the linear mean $\\mu_y[k]$ of the output signal $y[k] = \\mathcal{H} \\{ x[k] \\}$ of a linear time-invariant (LTI) system."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Non-Stationary Input Signal\n",
    "\n",
    "Let's first impose no restrictions in terms of stationarity to the input signal. The [linear mean](../random_signals/ensemble_averages.ipynb#Linear-mean) of the output signal is then given as\n",
    "\n",
    "\\begin{equation}\n",
    "\\mu_y[k] = E\\{ y[k] \\} = E\\{ x[k] * h[k] \\}\n",
    "\\end{equation}\n",
    "\n",
    "where $h[k]$ denotes the impulse response of the system. Since the convolution and the ensemble average are linear operations, and $h[k]$ is a deterministic signal, this can be rewritten as\n",
    "\n",
    "\\begin{equation}\n",
    "\\mu_y[k] = \\mu_x[k] * h[k]\n",
    "\\end{equation}\n",
    "\n",
    "The linear mean of the output signal $\\mu_y[k]$ is given as the convolution of the linear mean of the input signal $\\mu_x[k]$ with the impulse response $h[k]$ of the system."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Example\n",
    "\n",
    "The linear mean $\\mu_y[k]$ of the output of an LTI system with given impulse response $h[k]$ and non-stationary random input signal $x[k]$ is computed. The estimated linear means $\\hat{\\mu}_x[k]$ and $\\hat{\\mu}_y[k]$ of the input and output signals are plotted."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA2AAAAFACAYAAADaoR+7AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAzO0lEQVR4nO3de1yUdf7//+cIMogCnjkkKpknJPGUiuWhWklb+Wha6boadrDMrIy1g22l7raaZqWbWR9ty8wyP/3U0rU1rYRqxQIVNXP92H40SCEUFdASBK7vH/6YdWQ4DDDXMMPjfrvN7ea85zq8rus97wufc11zjcUwDEMAAAAAAJdr5O4CAAAAAKChIIABAAAAgEkIYAAAAABgEgIYAAAAAJiEAAYAAAAAJiGAAQAAAIBJCGAAAAAAYBICGAAAAACYhAAGAAAAACYhgAEAAACASQhgAADUI//+97/VtGlTDRkyRIZhuLscAEAdI4ABAFBPGIahe+65R4899pjOnz+vZcuWubskAEAdsxh8vAYAQL2wdOlSbd68Wdu2bdPRo0c1ZMgQffXVV7r66qvdXRoAoI5wBgyAW61atUoWi6XCR1JSUrWXtXPnTs2dO1dnz551uI5jx47Vae21qakuzJ07VxaLpcrpHG2/O/cJKvboo4/qs88+U6NGjdSpUycdP37cqfA1ZcoU29iJjo4u93peXp4aNWqkJUuWOJz/o48+sht/aWlpNd2UctatW6cePXqoSZMmslgsSk9PdzhdfXxv1vU4rm/bWN/qAbwdAQxAvfD2228rJSWl3KNPnz7VXsbOnTs1b968cv9J+u1vf6uUlBSFhYXVcdU1r8nd3LlP4FqhoaFKSUnR+++/X+61tLQ0GYah6667zuG8Q4cOVUpKip555pk6renkyZOaPHmyOnXqpK1btyolJUVdunRxOG19fG/W9Tiuj9sIwDy+7i4AACQpOjpa/fr1c8my27RpozZt2rhk2Z6qvu+TX375RQEBAe4uwyNZrVYNHDjQ4WtpaWny9fWt8IONFi1aaODAgfrXv/5VpzX97//+ry5evKhJkyZp6NChlU5b39+bdaEhbCOAinEGDIBHOHnypO6//35FRETIarWqTZs2uv766/XZZ59JunQ53uOPPy5JioyMtLuE0dHlNWWX7+3fv1933HGHgoOD1bJlSyUmJqq4uFiHDx/WiBEjFBgYqI4dO2rRokV29fzwww+6++671blzZwUEBOiqq65SfHy8Dhw4YLeOimoqc+TIEU2cOFFt27aV1WpV9+7d9dprr5Xb/i1btqhXr16yWq2KjIzU4sWLa7U/r9wnZfvj4MGD+t3vfqfg4GCFhITonnvuUV5eXrn5q1N3dfbR5eves2ePbr/9drVo0UKdOnWqsPba9p0rt6G6+8+RyMhIJSQklGu/6aabqgwt1ZWamqro6Gg1adJE0qWbfrzyyivy9/fXU089pZKSEqeX+fXXX+vmm29WYGCgAgICNGjQIG3ZssX2+pQpU3TDDTdIksaPHy+LxaJhw4ZVuLzKxmt19m/ZtHv37tXYsWMVFBSk4OBgTZo0SSdPnrSrq2PHjuXWf+WlvdUZx1eq6nhV0SV/H3/8sXr27Cmr1aqrr75aS5cudVhPdfdFdd+/AMzFGTAA9UJJSYmKi4vt2iwWi3x8fCRJkydP1p49e/SXv/xFXbp00dmzZ7Vnzx7l5uZKku677z6dPn1ar776qjZs2GC7tCcqKqrS7zXceeedmjRpkh544AFt375dixYt0sWLF/XZZ59p+vTpmjVrlt5//309+eSTuuaaazR27FhJ0okTJ9SqVSu98MILatOmjU6fPq133nlHAwYM0N69e9W1a9dKa5Kk77//XoMGDVL79u310ksvKTQ0VJ9++qkeeeQRnTp1SnPmzJEkff755xo9erRiY2P1wQcfqKSkRIsWLdLPP/9cdx3w/xs3bpzGjx+ve++9VwcOHNDs2bMlSW+99ZZtmurWXZ19dLmxY8dqwoQJmjZtms6fP19lrTXtO1duQ3X2nyO5ubk6duyYEhMT7doNw9CePXt0zz33VLk/qiMtLU1xcXGSpFOnTmnKlCnatWuXNm7cqJEjRzq9vOTkZA0fPlw9e/bU3/72N1mtVi1fvlzx8fFau3atxo8fr2effVb9+/fXQw89pPnz5+vGG29UUFBQjep3Zv/edtttuvPOOzVt2jQdPHhQzz77rL7//nt98803aty4cbXXWdU4dqSq45UjW7du1dixYzVkyBCtW7dOxcXFWrx4cYXjvDr7wtn3LwCTGADgRm+//bYhyeHDx8fHNl2zZs2MmTNnVrqsF1980ZBkHD161OE6Lm+fM2eOIcl46aWX7Kbt1auXIcnYsGGDre3ixYtGmzZtjLFjx1a47uLiYqOoqMjo3Lmz8dhjj1VZk2EYxi233GK0a9fOyMvLs2ufMWOG4e/vb5w+fdowDMMYMGCAER4ebvz666+2afLz842WLVsa1TmMO9r+K9vK9seiRYvs5p0+fbrh7+9vlJaWOl33lSraR2Xrfu6556rclsunr03fuWobqrP/HNm6dashydi5c6dd++HDhw1Jxpo1ayqdv0xCQoLRoUMHh6+dPHnSkGSsXLnSSE5ONq666irj+uuvNzIzM8tNW/b+SE1NrXR9AwcONNq2bWsUFBTY2oqLi43o6GijXbt2tu3esWOHIcn48MMPq9yGysZrdfZv2bSX949hGMZ7771nty8r2ldl81+usnHsSFXHK0fbeN111xkRERFGYWGhra2goMBo1aqVXT21ea9V9P51VA8A1+ESRAD1wurVq5Wammr3+Oabb2yv9+/fX6tWrdLzzz+vXbt26eLFi3Wy3lGjRtk97969uywWi93ZAF9fX11zzTX68ccfbW3FxcWaP3++oqKi5OfnJ19fX/n5+enIkSM6dOhQleu9cOGCPv/8c912220KCAhQcXGx7XHrrbfqwoUL2rVrl86fP6/U1FSNHTtW/v7+tvkDAwMVHx9fB3vA3n/913/ZPe/Zs6cuXLignJwcp+qWnN9H48aNc6rWmvadK7ehqv1XkdTUVPn6+qpXr1527bt375Yk9e3bt+odUoXU1FRJ0vbt23XzzTdr4sSJSkpKUrt27Wq0vPPnz+ubb77R7bffrmbNmtnafXx8NHnyZP300086fPhwreu+nDP79/e//73d8zvvvFO+vr7asWNHndbkiLPHq/PnzystLU1jxoyRn5+frb1Zs2YVjvPq7IvaHqcAuAYBDEC90L17d/Xr18/ucfl/OtetW6eEhAS9+eabio2NVcuWLXXXXXcpOzu7Vutt2bKl3XM/Pz8FBATYhZ2y9gsXLtieJyYm6tlnn9WYMWO0efNmffPNN0pNTVVMTIx+/fXXKtebm5ur4uJivfrqq2rcuLHd49Zbb5V06RKxM2fOqLS0VKGhoeWW4aittlq1amX33Gq1SpJtm6pbt+T8PnL2jnA17TtXbkNV+68iaWlpioqKsn036/L2Zs2a2d0x8KGHHtK9994rSSotLdXo0aP1yCOPVLr8smX5+/vrk08+0ZAhQ7Ro0SL5+tb8mwhnzpyRYRgO+y08PFySKr3kriac2b9Xjg9fX1+1atWqzmtyxNnjVdm+DAkJKfeaozapevuitscpAK7Bd8AAeITWrVtryZIlWrJkiTIyMrRp0yY99dRTysnJ0datW02vZ82aNbrrrrs0f/58u/ZTp06pefPmVc7fokUL25mChx56yOE0kZGR8vf3l8Vicfgft9qGz5qobt2S8/uoOr9pVhdcuQ01lZaWpuHDh5drT0pKUu/evdWo0X8+L509e7a6d++uOXPmaNmyZSopKdErr7xSrXX07t1bc+bM0ahRozR79mwtWLCgxjW3aNFCjRo1UlZWVrnXTpw4IenSuHWX7OxsXXXVVbbnxcXFys3NtQUXf39/FRYWlpuvLHzXhrPHqxYtWshisTj8vldtxrlZ718AziGAAfA47du314wZM/T555/rn//8p629umcb6oLFYrGtr8yWLVt0/PhxXXPNNVXWFBAQoBtvvFF79+5Vz5497S47ulL//v21YcMGvfjii7azOwUFBdq8eXNdbU61OVN3dfeR2erbNmRnZ+v48ePlzkYlJydrz549mjlzpl17u3btdNddd9kuQfv6669tN6upTGpqqsaNG6dbbrlFK1eu1N1336127dpVGEKr0rRpUw0YMEAbNmzQ4sWLbWfvSktLtWbNGrVr167C3/oyw3vvvWd3Fv1//ud/VFxcbLsDY8eOHZWTk6Off/7ZdpapqKhIn376abll1ebYUtHx6nJNmzZVv3799NFHH2nx4sW29+S5c+f097//3el1lqmvYxBo6AhgAOqF7777rtxdECWpU6dO8vPz04033qiJEyeqW7duCgwMVGpqqu2uYWWuvfZaSdLSpUuVkJCgxo0bu+wuX6NGjdKqVavUrVs39ezZU7t379aLL75Y7vs0FdUUGBiopUuX6oYbbtDgwYP14IMPqmPHjiooKNAPP/ygzZs364svvpAk/fnPf9aIESM0fPhw/eEPf1BJSYkWLlyopk2b6vTp0y7ZvspUt+7q7iN3qE/bUPbdrA8//FBRUVG65pprlJ6ebrslfk5Ojr777jtFR0fb5unVq5eWL1+unTt32n3/qiJZWVnKysqyBZIpU6bop59+0iOPPKLQ0FCnv39XZsGCBRo+fLhuvPFGzZo1S35+flq+fLm+++47rV271rSzmo5s2LBBvr6+Gj58uO0uiDExMbrzzjslXbol/nPPPacJEybo8ccf14ULF/TXv/7V4a34KxvHV8rLy6vW8epKf/rTn/Tb3/5Wt9xyix599FGVlJToxRdfVLNmzWo8zuvzGAQaMgIYgHrh7rvvdti+cuVKTZ48WQMGDNC7776rY8eO6eLFi2rfvr2efPJJPfHEE7Zphw0bptmzZ+udd97RypUrVVpa6rIv3C9dulSNGzfWggULdO7cOfXp00cbNmzQM888YzddRTUNGzZMUVFR2rNnj/785z/rmWeeUU5Ojpo3b67OnTvbvoskScOHD9dHH32kZ555RuPHj1doaKimT5+uX3/9VfPmzXPJ9lWmunVXdx+5Q33ahrIfR37zzTf1+OOPKzs7WwMHDtSmTZv0+9//Xjt27NCMGTNs03/77beaN2+e7rjjDr3zzjuKjY2tch1lIe/yM0LPPPOMMjMzNWnSJLVt21aDBw92uvahQ4fqiy++0Jw5czRlyhSVlpYqJiZGmzZtKneTFLNt2LBBc+fO1euvvy6LxaL4+HgtWbLEdnYpMjJSH3/8sZ5++mndfvvtCgsLU2Jiok6ePFluXFU2jq/k7+9frePVlUaMGKH169frueeesxvnJ06c0LvvvlujfVCfxyDQkFkMwzDcXQQAAA3VrbfequzsbO3Zs6fKaTMyMnTDDTdo9erV6tq1q7p27aoDBw6oQ4cOki6d2UpKStIPP/xg9zt61WUYhkpKSrR69Wrde++9Sk1NVb9+/Wq0Xe4yd+5czZs3TydPnnTrd9DqwsWLF9WrVy9dddVV2rZtm7vLAVBHuAsiAAButHv3bvXv37/K6QoKCjRq1CjNmTNHw4YNU1hYmCZNmqTnn3/ebroff/xRjRs3VkxMjNO1fPzxx2rcuLHtLosw17333qsPPvhAycnJWrduneLi4nTo0KFKz5wB8DxcgggAgJtkZGQoJyenWgEsMDBQ+/fvt2tbvny53fO5c+faLle88pb21TFs2DDb5YrSpUs1YZ6CggLNmjVLJ0+eVOPGjdWnTx998skn+s1vfuPu0gDUIS5BBAAAAACTcAkiAAAAAJiEAAYAAAAAJiGAAQAAAIBJCGAAAAAAYBLuglhDpaWlOnHihAIDA2WxWNxdDgAAAAA3MQxDBQUFCg8PV6NGlZ/jIoDV0IkTJxQREeHuMgAAAADUE5mZmWrXrl2l0xDAaigwMFDSpZ0cFBTk5moAAAAAuEt+fr4iIiJsGaEyBLAaKrvsMCgoiAAGAAAAoFpfTeImHAAAAABgEgIYAAAAAJiEAAYAAAAAJiGAAQAAAIBJCGAAAAAAYBICGAAAAACYhAAGAAAAACYhgAEAAACASQhgAAAAAGASAhgAAAAAmMTX3QUAAOCpSkoNfXv0tHIKLqhtoL/6R7aUTyOLu8sCANRjBDAAAGpg63dZmrf5e2XlXbC1hQX7a058lEZEh7mxMgBAfeYVlyB++eWXio+PV3h4uCwWiz766KMq50lOTlbfvn3l7++vq6++Wm+88YbrCwUAeIWt32XpwTV77MKXJGXnXdCDa/Zo63dZbqoMAFDfeUUAO3/+vGJiYrRs2bJqTX/06FHdeuutGjx4sPbu3aunn35ajzzyiNavX+/iSgEAnq6k1NC8zd/LcPBaWdu8zd+rpNTRFACAhs4rLkEcOXKkRo4cWe3p33jjDbVv315LliyRJHXv3l1paWlavHixxo0b56IqAQDe4Nujp8ud+bqcISkr74K+PXpasZ1amVcYAMAjeMUZMGelpKQoLi7Oru2WW25RWlqaLl686HCewsJC5efn2z0AAA1PTkHF4asm0wEAGpYGGcCys7MVEhJi1xYSEqLi4mKdOnXK4TwLFixQcHCw7REREWFGqQCAeqZtoH+dTgcAaFgaZACTJIvF/jbBhmE4bC8ze/Zs5eXl2R6ZmZkurxEAUP/0j2ypsGB/VXSzeYsu3Q2xf2RLM8sCAHiIBhnAQkNDlZ2dbdeWk5MjX19ftWrl+Hp9q9WqoKAguwcAoOHxaWTRnPgoSSoXwsqez4mP4vfAAAAONcgAFhsbq+3bt9u1bdu2Tf369VPjxo3dVBUAwFOMiA7T65P6qG2Q1a49NNhfr0/qw++AAQAq5BUB7Ny5c0pPT1d6erqkS7eZT09PV0ZGhqRLlw/eddddtumnTZumH3/8UYmJiTp06JDeeust/e1vf9OsWbPcUT4AwAONiA7TZ4lDbc9X3X2dvn7yJsIXAKBSXnEb+rS0NN14442254mJiZKkhIQErVq1SllZWbYwJkmRkZH65JNP9Nhjj+m1115TeHi4/vrXv3ILegCAUy6/zLB/ZEsuOwQAVMkrAtiwYcNsN9FwZNWqVeXahg4dqj179riwKgAAAACw5xWXIAIAAACAJyCAAQAAAIBJCGAAAAAAYBICGAAAAACYhAAGAAAAACYhgAEAAACASQhgAAAAAGASAhgAAAAAmIQABgAAAAAmIYABAAAAgEkIYAAAAABgEgIYAAAAAJiEAAYAAAAAJiGAAQAAAIBJCGAAAAAAYBICGAAAAACYhAAGAAAAACYhgAEAAACASQhgAAAAAGASAhgAAAAAmIQABgAAAAAmIYABAAAAgEkIYAAAAABgEgIYAAAAAJiEAAYAAAAAJiGAAQAAAIBJCGAAAAAAYBICGAAAAACYhAAGAAAAACYhgAEAAACASQhgAAAAAGASAhgAAAAAmIQABgAAAAAmIYABAAAAgEkIYAAAAABgEgIYAAAAAJiEAAYAAAAAJvGaALZ8+XJFRkbK399fffv21VdffVXhtElJSbJYLOUe//rXv0ysGAAAAEBD4xUBbN26dZo5c6b++Mc/au/evRo8eLBGjhypjIyMSuc7fPiwsrKybI/OnTubVDEAAACAhsgrAtjLL7+se++9V/fdd5+6d++uJUuWKCIiQq+//nql87Vt21ahoaG2h4+Pj0kVAwAAAGiIPD6AFRUVaffu3YqLi7Nrj4uL086dOyudt3fv3goLC9PNN9+sHTt2VDptYWGh8vPz7R4AAAAA4AyPD2CnTp1SSUmJQkJC7NpDQkKUnZ3tcJ6wsDCtWLFC69ev14YNG9S1a1fdfPPN+vLLLytcz4IFCxQcHGx7RERE1Ol2AAAAAPB+vu4uoK5YLBa754ZhlGsr07VrV3Xt2tX2PDY2VpmZmVq8eLGGDBnicJ7Zs2crMTHR9jw/P58QBgAAAMApHn8GrHXr1vLx8Sl3tisnJ6fcWbHKDBw4UEeOHKnwdavVqqCgILsHAAAAADjD4wOYn5+f+vbtq+3bt9u1b9++XYMGDar2cvbu3auwsLC6Lg8AAAAAbLziEsTExERNnjxZ/fr1U2xsrFasWKGMjAxNmzZN0qXLB48fP67Vq1dLkpYsWaKOHTuqR48eKioq0po1a7R+/XqtX7/enZsBAAAAwMt5RQAbP368cnNz9ac//UlZWVmKjo7WJ598og4dOkiSsrKy7H4TrKioSLNmzdLx48fVpEkT9ejRQ1u2bNGtt97qrk0AAAAA0ABYDMMw3F2EJ8rPz1dwcLDy8vL4PhgANFC/FBUr6rlPJUnf/+kWBfh5xeeaAAAnOZMNPP47YAAAAADgKQhgAAAAAGASAhgAAAAAmIQABgAAAAAmIYABAAAAgEkIYAAAAABgEgIYAAAAAJiEAAYAAAAAJiGAAQAAAIBJCGAAAAAAYBICGAAAAACYhAAGAAAAACYhgAEAAACASQhgAAAAAGASAhgAAAAAmIQABgAAAAAmIYABAAAAgEkIYAAAAABgEgIYAAAAAJiEAAYAAAAAJiGAAQAAAIBJCGAAAAAAYBICGAAAAACYhAAGAAAAACYhgAEAAACASQhgAAAAAGASAhgAAAAAmIQABgAAAAAmIYABAAAAgEl8nZ1h06ZNTq9k+PDhatKkidPzAQAAAIA3cTqAjRkzxqnpLRaLjhw5oquvvtrZVQEAAACAV6nRJYjZ2dkqLS2t1iMgIKCuawYAAAAAj+R0AEtISHDqcsJJkyYpKCjI2dUAAAAAgNdx+hLEt99+2/bvn3/+WSEhIZVO//rrrztfFQAAAAB4oVrdBXHcuHEqLi52+FpF7QDQUJSUGkr5d64+Tj+ulH/nqqTUcHdJAADAzZw+A3a5Fi1a6OGHHy53lis3N1fjxo1TUlJSbRYPAB5r63dZmrf5e2XlXbC1hQX7a058lEZEh7mxMgANSUmpoW+PnlZOwQW1DfRX/8iW8mlkcXdZDYqr+4A+9jy1CmDvvvuu+vfvrzfffFP33XefJOnQoUMaNWqUevToUScF1nff/t9p3dgzkDc6nObpB2RPP+C7sv6t32XpwTV7dOX5ruy8C3pwzR69PqlPnYQwT+8DVI1x7H6evI+85YMg+sB9yzdLQzvWWQzDqNU1MQcOHNDQoUP1j3/8Q2fOnNGECRN0//33a+HChbJYzNuw5cuX68UXX1RWVpZ69OihJUuWaPDgwRVOn5ycrMTERB08eFDh4eF64oknNG3atGqvLz8/X8HBwYqY+T+6qm3LOn2jm/EmqW9vRGd5ev2efkD29AO+K+svKTV0w8Iv7JZ9OYuk0GB/ff3kTbV6z3p6H0jecaz7pahYUc99Kkn6/k+3KMCvVp9r2vGGcezpfezJ+6iiD4LKllxXHwRJ9EFFXN0HZvaxK3nDsU76TzbIy8ur8gaETgew0aNHq1evXurdu7d69eqljh07au3atXr44Yd14cIFvfbaa0pISKjVBjhr3bp1mjx5spYvX67rr79e//3f/60333xT33//vdq3b19u+qNHjyo6OlpTp07VAw88oH/+85+aPn261q5dq3HjxlVrnZcHMB/rpVvt18Ub3Yw3iScfzCTvqN+TD8hmHfA99Q9iyr9z9buVu6qcbu3UgYrt1KpG6/D0PpC851jnqgDmDePY0/vYk/eRWR8ESfRBRVzdB2b2cdn6PPFvspkh1aUBbNasWUpPT9e+ffuUm5ur5s2bKyYmRvv379ftt9+u6dOnKyoqSo0bN67VRjhjwIAB6tOnj9130bp3764xY8ZowYIF5aZ/8skntWnTJh06dMjWNm3aNO3bt08pKSnVWmfZTu404135WANkkRQS5K/PEofW+A25/ftsPfpBeoVvkqUTeml4VGiNlm32OuZ/8i9l5//noBAa5K+nb+3W4OsvKTX0m5eT7ZZ9udq+jzx9+WVc1Qdm1P/3/Sf0+P+3v8rpXry9p0b1DHd6+Z7eB2XL9oZjnXQpgPV9/jNJ0u5nflMnAcwbxrGn97Gn76Nvjp7WlLe/rXK6VXf314DIljVah0QfVMbVfWBWH0ue+ze5ouUX+vhJFkudh1SXBrDL/fTTT0pPT7d7HD16VL6+vurWrZv27dtX00VXW1FRkQICAvThhx/qtttus7U/+uijSk9PV3Jycrl5hgwZot69e2vp0qW2to0bN+rOO+/UL7/84jA8FhYWqrCw0PY8Pz9fERER+vaazmrm41PHWwUAAACgro0Z9RcV+lptz2tzRcrlnAlgtfqorl27dmrXrp1GjRplazt37pz27t2r/fur/gS4Lpw6dUolJSXlfo8sJCRE2dnZDufJzs52OH1xcbFOnTqlsLDypyIXLFigefPm1V3hAAAAANwqp8DxGThXcjqA7d+/X9HR0WrUyPFPiDVr1kyDBw+23QDj4MGD6tq1q3x96+6LyY5cecMPwzAqvQmIo+kdtZeZPXu2EhMTbc/LzoD9bsRztu+ASTU/1evqy5bMWIerT4d7ev1mrIPlu3f5Zcoua5Fkd2kLl7V4x7HO1ehj96/D0/dR2aVXP+dfKHd5neQZl1vTB+5dvuT5x6KKll/o42f3vG2gv9PLri2nU1Hv3r2VnZ2tNm3aVGv62NhYpaen6+qrr3a6uOpo3bq1fHx8yp3tysnJKXeWq0xoaKjD6X19fdWqleNTkFarVVartVx7ka9VjXyttutI+3e/So1q8EZv06aF3enQyqZrFBBQ5XTuWEfOxTPVWn7ORUuNlu/p9UtS/+5N1LJVsLLzKj5g1uZ95OnLd3UfuLr+Mrf0u1pL/JtU+MXuW2rxhV9P7wNvONa5mqePY2/oY0/fR40kzR7bWw+u2SPJ8QdBs8f2VuNmTZ1e9uW10QcVc3UfmNHHnv43udrLr+V35GrC6QBmGIaeffZZBVRzRxcVFTldlDP8/PzUt29fbd++3e47YNu3b9fo0aMdzhMbG6vNmzfbtW3btk39+vWr0c1Dyt4Sc+KjavwpQ//IlgoL9nfpm8TV66juJwg1/aTB0+uXJJ9GFs2Jj9KDa/bIIscHzNq8jzx9+a7uA1fXf7kR0WEaHhVa53eN8vQ+8IZjnat5+jj2hj72hn00IjpMr0/qU+6DoNA6ukMhfVA1V/eBq5fv6X+Tzfyb7yynb8IxbNgwp3/f6/3333f4vaq6UnYb+jfeeEOxsbFasWKFVq5cqYMHD6pDhw6aPXu2jh8/rtWrV0v6z23oH3jgAU2dOlUpKSmaNm1ajW9DX1e/A1Z2q0zJ8Zukrm636qp1lN0StaqDWW3uNuPp9Zfx9N+8cPVte13dB97wG1qe3Aeefqwzi6eO47Jle0Mfe/o+ksz5WRWJPqiMp/7IsLf8TfaK3wGrr5YvX65FixYpKytL0dHReuWVVzRkyBBJ0pQpU3Ts2DElJSXZpk9OTtZjjz1m+yHmJ598skY/xLx971Hd2LMDv41z2bI5GFePpx6QXb18b/mDaAZP7gNPP9aZxVPHseQ9fezp+8jV6APv5i1/k834m29aAMvMzFRERERNZ/dozuxkZ5nxJvH0g5mn14/K0Qfu5+nj2Mx1oGL0cdU8vX7J87fB0+t3Nf4mV4/LAtjatWv1u9/9zva8UaNGatGihWJiYhQTE6NevXopJiZGhYWFeu2112yX/HkjVwYwb+DpBzNPr98b0AfuRx8AACT+HlRHnQew7OxsTZ8+Xc2bN9dbb71laz927JjtB5j37t2rPXv26MSJE5KkoKAgnTlzppabUn8RwAAAAABILvgh5hUrVqi4uNgufElSx44d1bFjR40ZM8bWlpKSooSEBC1cuND5ygEAAADAizn+NeUrPPLII2revHm17hAYGxurpUuX6vnnn691cQAAAADgTap1Bqx58+ZavXp1ud/OunjxosPfzercubMOHjxYNxUCAAAAgJdw6oeY4+Pj7Z43bdpUUVFR6t27t3r16qXevXsrPDxcr776quLi4uq0UAAAAADwdLW6Df3XX3+tffv2ad++fUpPT9fBgwf166+/SpLi4uLUt29f9ezZUz179lT37t3rrOj6gJtwAAAAAJDc+EPMpaWlOnz4sO3OiGXhLCcnRyUlJXW1mnqBAAYAAABAcmMAq8jPP/+skJAQV6/GVAQwAAAAAJJz2aBad0GsLW8LXwAAAABQE6YEMAAAAAAAAQwAAAAATEMAAwAAAACTEMAAAAAAwCQEMAAAAAAwCQEMAAAAAExCAAMAAAAAkxDAAAAAAMAkBDAAAAAAMAkBDAAAAABMQgADAAAAAJMQwAAAAADAJAQwAAAAADAJAQwAAAAATEIAAwAAAACTEMAAAAAAwCQEMAAAAAAwCQEMAAAAAExCAAMAAAAAkxDAAAAAAMAkBDAAAAAAMAkBDAAAAABMQgADAAAAAJMQwAAAAADAJAQwAAAAADAJAQwAAAAATEIAAwAAAACTEMAAAAAAwCQeH8DOnDmjyZMnKzg4WMHBwZo8ebLOnj1b6TxTpkyRxWKxewwcONCcggEAAAA0WL7uLqC2Jk6cqJ9++klbt26VJN1///2aPHmyNm/eXOl8I0aM0Ntvv2177ufn59I6AQAAAMCjA9ihQ4e0detW7dq1SwMGDJAkrVy5UrGxsTp8+LC6du1a4bxWq1WhoaHVXldhYaEKCwttz/Pz82teOAAAAIAGyaMvQUxJSVFwcLAtfEnSwIEDFRwcrJ07d1Y6b1JSktq2basuXbpo6tSpysnJqXT6BQsW2C5zDA4OVkRERJ1sAwAAAICGw6MDWHZ2ttq2bVuuvW3btsrOzq5wvpEjR+q9997TF198oZdeekmpqam66aab7M5wXWn27NnKy8uzPTIzM+tkGwAAAAA0HPXyEsS5c+dq3rx5lU6TmpoqSbJYLOVeMwzDYXuZ8ePH2/4dHR2tfv36qUOHDtqyZYvGjh3rcB6r1Sqr1Vqd8gEAAADAoXoZwGbMmKEJEyZUOk3Hjh21f/9+/fzzz+VeO3nypEJCQqq9vrCwMHXo0EFHjhxxulYAAAAAqK56GcBat26t1q1bVzldbGys8vLy9O2336p///6SpG+++UZ5eXkaNGhQtdeXm5urzMxMhYWF1bhmAAAAAKiKR38HrHv37hoxYoSmTp2qXbt2adeuXZo6dapGjRpldwfEbt26aePGjZKkc+fOadasWUpJSdGxY8eUlJSk+Ph4tW7dWrfddpu7NgUAAABAA+DRAUyS3nvvPV177bWKi4tTXFycevbsqXfffddumsOHDysvL0+S5OPjowMHDmj06NHq0qWLEhIS1KVLF6WkpCgwMNAdmwAAAACggbAYhmG4uwhPlJ+fr+DgYOXl5SkoKMjd5QAAAABwE2eygcefAQMAAAAAT0EAAwAAAACTEMAAAAAAwCQEMAAAAAAwCQEMAAAAAExCAAMAAAAAkxDAAAAAAMAkBDAAAAAAMAkBDAAAAABMQgADAAAAAJMQwAAAAADAJAQwAAAAADAJAQwAAAAATEIAAwAAAACTEMAAAAAAwCQEMAAAAAAwCQEMAAAAAExCAAMAAAAAkxDAAAAAAMAkBDAAAAAAMAkBDAAAAABMQgADAAAAAJMQwAAAAADAJAQwAAAAADAJAQwAAAAATEIAAwAAAACTEMAAAAAAwCQEMAAAAAAwCQEMAAAAAExCAAMAAAAAkxDAAAAAAMAkBDAAAAAAMAkBDAAAAABMQgADAAAAAJMQwAAAAADAJAQwAAAAADAJAQwAAAAATOLxAewvf/mLBg0apICAADVv3rxa8xiGoblz5yo8PFxNmjTRsGHDdPDgQdcWCgAAAKDB8/gAVlRUpDvuuEMPPvhgtedZtGiRXn75ZS1btkypqakKDQ3V8OHDVVBQ4MJKAQAAADR0FsMwDHcXURdWrVqlmTNn6uzZs5VOZxiGwsPDNXPmTD355JOSpMLCQoWEhGjhwoV64IEHHM5XWFiowsJC2/P8/HxFREQoLy9PQUFBdbYdAAAAADxLfn6+goODq5UNPP4MmLOOHj2q7OxsxcXF2dqsVquGDh2qnTt3VjjfggULFBwcbHtERESYUS4AAAAAL9LgAlh2drYkKSQkxK49JCTE9pojs2fPVl5enu2RmZnp0joBAAAAeJ96GcDmzp0ri8VS6SMtLa1W67BYLHbPDcMo13Y5q9WqoKAguwcAAAAAOMPX3QU4MmPGDE2YMKHSaTp27FijZYeGhkq6dCYsLCzM1p6Tk1PurBgAAAAA1KV6GcBat26t1q1bu2TZkZGRCg0N1fbt29W7d29Jl+6kmJycrIULF7pknQAAAAAg1dNLEJ2RkZGh9PR0ZWRkqKSkROnp6UpPT9e5c+ds03Tr1k0bN26UdOnSw5kzZ2r+/PnauHGjvvvuO02ZMkUBAQGaOHGiuzYDAAAAQANQL8+AOeO5557TO++8Y3tedlZrx44dGjZsmCTp8OHDysvLs03zxBNP6Ndff9X06dN15swZDRgwQNu2bVNgYKCptQMAAABoWLzmd8DM5sy9/gEAAAB4L34HDAAAAADqIQIYAAAAAJiEAAYAAAAAJiGAAQAAAIBJCGAAAAAAYBICGAAAAACYhAAGAAAAACYhgAEAAACASQhgAAAAAGASAhgAAAAAmIQABgAAAAAm8XV3AZ7KMAxJUn5+vpsrAQAAAOBOZZmgLCNUhgBWQ7m5uZKkiIgIN1cCAAAAoD4oKChQcHBwpdMQwGqoZcuWkqSMjIwqdzI8U35+viIiIpSZmamgoCB3lwMXoI8bBvrZ+9HH3o8+9n6e3seGYaigoEDh4eFVTksAq6FGjS59fS44ONgj3ySovqCgIPrYy9HHDQP97P3oY+9HH3s/T+7j6p6U4SYcAAAAAGASAhgAAAAAmIQAVkNWq1Vz5syR1Wp1dylwEfrY+9HHDQP97P3oY+9HH3u/htTHFqM690oEAAAAANQaZ8AAAAAAwCQEMAAAAAAwCQEMAAAAAExCAAMAAAAAkxDAamj58uWKjIyUv7+/+vbtq6+++srdJaGOzJ07VxaLxe4RGhrq7rJQC19++aXi4+MVHh4ui8Wijz76yO51wzA0d+5chYeHq0mTJho2bJgOHjzonmJRI1X18ZQpU8qN64EDB7qnWNTIggULdN111ykwMFBt27bVmDFjdPjwYbtpGMuerTp9zFj2bK+//rp69uxp+7Hl2NhY/eMf/7C93lDGMAGsBtatW6eZM2fqj3/8o/bu3avBgwdr5MiRysjIcHdpqCM9evRQVlaW7XHgwAF3l4RaOH/+vGJiYrRs2TKHry9atEgvv/yyli1bptTUVIWGhmr48OEqKCgwuVLUVFV9LEkjRoywG9effPKJiRWitpKTk/XQQw9p165d2r59u4qLixUXF6fz58/bpmEse7bq9LHEWPZk7dq10wsvvKC0tDSlpaXppptu0ujRo20hq8GMYQNO69+/vzFt2jS7tm7duhlPPfWUmypCXZozZ44RExPj7jLgIpKMjRs32p6XlpYaoaGhxgsvvGBru3DhghEcHGy88cYbbqgQtXVlHxuGYSQkJBijR492Sz1wjZycHEOSkZycbBgGY9kbXdnHhsFY9kYtWrQw3nzzzQY1hjkD5qSioiLt3r1bcXFxdu1xcXHauXOnm6pCXTty5IjCw8MVGRmpCRMm6P/+7//cXRJc5OjRo8rOzrYb01arVUOHDmVMe5mkpCS1bdtWXbp00dSpU5WTk+PuklALeXl5kqSWLVtKYix7oyv7uAxj2TuUlJTogw8+0Pnz5xUbG9ugxjABzEmnTp1SSUmJQkJC7NpDQkKUnZ3tpqpQlwYMGKDVq1fr008/1cqVK5Wdna1BgwYpNzfX3aXBBcrGLWPau40cOVLvvfeevvjiC7300ktKTU3VTTfdpMLCQneXhhowDEOJiYm64YYbFB0dLYmx7G0c9bHEWPYGBw4cULNmzWS1WjVt2jRt3LhRUVFRDWoM+7q7AE9lsVjsnhuGUa4NnmnkyJG2f1977bWKjY1Vp06d9M477ygxMdGNlcGVGNPebfz48bZ/R0dHq1+/furQoYO2bNmisWPHurEy1MSMGTO0f/9+ff311+VeYyx7h4r6mLHs+bp27ar09HSdPXtW69evV0JCgpKTk22vN4QxzBkwJ7Vu3Vo+Pj7lknhOTk65xA7v0LRpU1177bU6cuSIu0uBC5Td4ZIx3bCEhYWpQ4cOjGsP9PDDD2vTpk3asWOH2rVrZ2tnLHuPivrYEcay5/Hz89M111yjfv36acGCBYqJidHSpUsb1BgmgDnJz89Pffv21fbt2+3at2/frkGDBrmpKrhSYWGhDh06pLCwMHeXAheIjIxUaGio3ZguKipScnIyY9qL5ebmKjMzk3HtQQzD0IwZM7RhwwZ98cUXioyMtHudsez5qupjRxjLns8wDBUWFjaoMcwliDWQmJioyZMnq1+/foqNjdWKFSuUkZGhadOmubs01IFZs2YpPj5e7du3V05Ojp5//nnl5+crISHB3aWhhs6dO6cffvjB9vzo0aNKT09Xy5Yt1b59e82cOVPz589X586d1blzZ82fP18BAQGaOHGiG6uGMyrr45YtW2ru3LkaN26cwsLCdOzYMT399NNq3bq1brvtNjdWDWc89NBDev/99/Xxxx8rMDDQ9il5cHCwmjRpIovFwlj2cFX18blz5xjLHu7pp5/WyJEjFRERoYKCAn3wwQdKSkrS1q1bG9YYdtv9Fz3ca6+9ZnTo0MHw8/Mz+vTpY3eLVHi28ePHG2FhYUbjxo2N8PBwY+zYscbBgwfdXRZqYceOHYakco+EhATDMC7dvnrOnDlGaGioYbVajSFDhhgHDhxwb9FwSmV9/MsvvxhxcXFGmzZtjMaNGxvt27c3EhISjIyMDHeXDSc46l9Jxttvv22bhrHs2arqY8ay57vnnnts/39u06aNcfPNNxvbtm2zvd5QxrDFMAzDzMAHAAAAAA0V3wEDAAAAAJMQwAAAAADAJAQwAAAAADAJAQwAAAAATEIAAwAAAACTEMAAAAAAwCQEMAAAAAAwCQEMAAAAAExCAAMAAAAAkxDAAACooT/84Q+Kj493dxkAAA9CAAMAoIbS09PVq1cvd5cBAPAgBDAAAGpo37596t27t7vLAAB4EAIYAAA1kJmZqdzcXNsZsLNnzyo+Pl6DBg1SVlaWe4sDANRbBDAAAGogPT1dwcHBioyM1IEDB3TdddcpLCxMSUlJCgsLc3d5AIB6igAGAEANpKenKyYmRmvXrtWQIUM0a9YsrVixQn5+fu4uDQBQj1kMwzDcXQQAAJ5m3Lhx2rFjhyTp73//uwYNGuTmigAAnoAzYAAA1EB6errGjRunCxcu6OzZs+4uBwDgITgDBgCAkwoKChQcHKzdu3dr3759evTRR7Vz50716NHD3aUBAOo5X3cXAACAp0lPT5ePj4+ioqLUu3dvHTx4UPHx8fr222/VunVrd5cHAKjHuAQRAAAn7du3T926dZPVapUkLVy4UFFRURo7dqyKiorcXB0AoD7jEkQAAAAAMAlnwAAAAADAJAQwAAAAADAJAQwAAAAATEIAAwAAAACTEMAAAAAAwCQEMAAAAAAwCQEMAAAAAExCAAMAAAAAkxDAAAAAAMAkBDAAAAAAMAkBDAAAAABM8v8AMwPx40toGLYAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 1000x300 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA2AAAAFBCAYAAAAR/cweAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA5lklEQVR4nO3de3gU9b3H8c8mIQmBZCFAbhAgKLcQAkiMBLmpBUHJUaEVtNBQrRaRKqZUi60CPRYErYUWwYOtIsWqRwGVahFUiLbcEiDIrRTbKIgJgQBJQJOQZM4fnGxZcttNZmc3yfv1PPs87G/n8t3f/mbYT2Z2xmYYhiEAAAAAgMf5ebsAAAAAAGgpCGAAAAAAYBECGAAAAABYhAAGAAAAABYhgAEAAACARQhgAAAAAGARAhgAAAAAWIQABgAAAAAWIYABAAAAgEUIYAAAAABgEQIYAAAAAFiEAAYAgBf861//Ups2bTRixAgZhuHtcgAAFiGAAQBgMcMwdM899+iRRx7RhQsXtGzZMm+XBACwiM3gz24AAFhq6dKl2rBhgzZt2qScnByNGDFCn376qXr06OHt0gAAHsYRMACWWLVqlWw2W62PrVu3urysbdu2ad68eTp37lyN6/jiiy9Mrb0xNZlh3rx5stls9U5X0/v3Zp+gdg8//LA+/PBD+fn56aqrrtKJEyfcCl/Tpk1zbDsJCQnVXi8sLJSfn5+WLFlS4/xvv/220/aXlZXV0LdSzRtvvKF+/fqpdevWstlsys7ONm3Z9fHkdmjF+nxte/W1eoDmggAGwFIvv/yytm/fXu1xzTXXuLyMbdu2af78+dW+9Nx6663avn27oqOjTa664TV5mzf7BJ4VFRWl7du3689//nO117KysmQYhq699toa5x05cqS2b9+uX/7yl6bWdOrUKU2dOlVXXXWVNm7cqO3bt6tXr16mrqMuVm+HZq+P7RVoGQK8XQCAliUhIUFJSUkeWXanTp3UqVMnjyy7qfL1Pvnmm28UEhLi7TKapKCgIA0ZMqTG17KyshQQEFDrHzbat2+vIUOG6B//+IepNf3zn//UxYsXNWXKFI0cOdLUZbcEvr69AjAHR8AA+JRTp07p/vvvV2xsrIKCgtSpUyddf/31+vDDDyVdOh3vZz/7mSQpLi7O6RTGmk6XqTp977PPPtP3vvc92e12hYeHKz09XeXl5Tpy5IjGjh2r0NBQde/eXYsXL3aq5/PPP9cPf/hD9ezZUyEhIercubNSU1O1f/9+p3XUVlOVo0eP6u6771ZERISCgoLUt29fPf/889Xe/3vvvaeBAwcqKChIcXFxevbZZxvVn1f2SVV/HDx4UHfddZfsdrsiIyN1zz33qLCwsNr8rtTtSh9dvu49e/bou9/9rtq3b6+rrrqq1tob+9l58j242n81iYuLU1paWrX2G2+80bTQkpmZqYSEBLVu3VrSpYt+/Pa3v1VwcLB+/vOfq6Kiwu1l/u1vf9NNN92k0NBQhYSEaOjQoXrvvfccr0+bNk3Dhg2TJE2aNEk2m02jRo1q1DIvX3b37t2rtV9+em5922HVtHv37tWECRMUFhYmu92uKVOm6NSpU26ty5X11aS+/Vttp/y98847SkxMVFBQkHr06KGlS5fWeGqyq+PT1fEOwDM4AgbAUhUVFSovL3dqs9ls8vf3lyRNnTpVe/bs0a9//Wv16tVL586d0549e1RQUCBJ+tGPfqQzZ87o97//vdatW+c4VSc+Pr7O3ynceeedmjJlin784x9r8+bNWrx4sS5evKgPP/xQM2bM0OzZs/XnP/9Zjz32mK6++mpNmDBBkvT111+rQ4cOevrpp9WpUyedOXNGr7zyiq677jrt3btXvXv3rrMmSTp06JCGDh2qrl276je/+Y2ioqL0wQcf6KGHHtLp06c1d+5cSdJHH32k2267TSkpKXr99ddVUVGhxYsX6+TJk+Z9AP9v4sSJmjRpku69917t379fc+bMkSS99NJLjmlcrduVPrrchAkTNHnyZE2fPl0XLlyot9aGfnaefA+u9F9NCgoK9MUXXyg9Pd2p3TAM7dmzR/fcc0+9/eGKrKwsjRkzRpJ0+vRpTZs2TTt27ND69es1btw4t5eXkZGh0aNHKzExUX/84x8VFBSk5cuXKzU1Va+99pomTZqkJ554QsnJyXrwwQe1YMEC3XDDDQoLC2vUMt1R33ZY5Y477tCdd96p6dOn6+DBg3riiSd06NAh7dy5U61atTJ9fZerb/9Wk40bN2rChAkaMWKE3njjDZWXl+vZZ5+tc79Q3/h0d7wDMJkBABZ4+eWXDUk1Pvz9/R3TtW3b1pg1a1ady3rmmWcMSUZOTk6N67i8fe7cuYYk4ze/+Y3TtAMHDjQkGevWrXO0Xbx40ejUqZMxYcKEWtddXl5ulJWVGT179jQeeeSRemsyDMO4+eabjS5duhiFhYVO7TNnzjSCg4ONM2fOGIZhGNddd50RExNjfPvtt45pioqKjPDwcMOV3XVN7//Ktqr+WLx4sdO8M2bMMIKDg43Kykq3675SbX1Ute4nn3yy3vdy+fSN+ew89R5c6b+abNy40ZBkbNu2zan9yJEjhiRjzZo1dc5fJS0tzejWrVuNr506dcqQZLz44otGRkaG0blzZ+P66683jh8/Xm3aqvGRmZlZ5/qGDBliREREGMXFxY628vJyIyEhwejSpYvjfW/ZssWQZLz55pv1vgdXl1nX+636PKrUtR1WTXv552kYhvHqq6869b2r66pvfTWpb/9W0zZ87bXXGrGxsUZpaamjrbi42OjQoUO1eho6Pmsb7zXVA6DxOAURgKVWr16tzMxMp8fOnTsdrycnJ2vVqlV66qmntGPHDl28eNGU9Y4fP97ped++fWWz2ZyOBgQEBOjqq6/Wl19+6WgrLy/XggULFB8fr8DAQAUEBCgwMFBHjx7V4cOH611vSUmJPvroI91xxx0KCQlReXm543HLLbeopKREO3bs0IULF5SZmakJEyYoODjYMX9oaKhSU1NN6AFn//Vf/+X0PDExUSUlJcrPz3erbsn9Ppo4caJbtTb0s/Pke6iv/2qTmZmpgIAADRw40Kl99+7dkqTBgwfX3yH1yMzMlCRt3rxZN910k+6++25t3bpVXbp0adDyLly4oJ07d+q73/2u2rZt62j39/fX1KlT9dVXX+nIkSNeX6arvv/97zs9v/POOxUQEKAtW7Z4ZH2Xc3f/duHCBWVlZen2229XYGCgo71t27Z17hfqG5+N3a8BaBwCGABL9e3bV0lJSU6Py790vvHGG0pLS9Mf/vAHpaSkKDw8XD/4wQ+Ul5fXqPWGh4c7PQ8MDFRISIhT2KlqLykpcTxPT0/XE088odtvv10bNmzQzp07lZmZqQEDBujbb7+td70FBQUqLy/X73//e7Vq1crpccstt0i6dIrY2bNnVVlZqaioqGrLqKmtsTp06OD0PCgoSJIc78nVuiX3+8jdK7w19LPz5Huor/9qk5WVpfj4eMdvsy5vb9u2reOKgefOnZO/v7+Ki4sd07z11luO31jVt47g4GC9//77GjFihBYvXqyAgIb/4uDs2bMyDKPGzy0mJkaS6jyFzqpluurK7SkgIEAdOnTw2Pou5+7+raqfIiMjq71WU1uV+sZnY/drABqH34AB8CkdO3bUkiVLtGTJEh07dkzvvvuufv7znys/P18bN260vJ41a9boBz/4gRYsWODUfvr0abVr167e+du3b+/4q/6DDz5Y4zRxcXEKDg6WzWar8YtYY8NnQ7hat+R+H7lyTzMzePI9NFRWVpZGjx5drX3r1q0aNGiQ/Pwu/V20Xbt26tGjh/bt26dhw4apoqJCTz75pFasWOHSOgYNGqS5c+dq/PjxmjNnjhYuXNjgmtu3by8/Pz/l5uZWe+3rr7+WdGm79eQyg4ODVVpaWm3aqgDtjry8PHXu3NnxvLy8XAUFBY7QYua6ruTu/q19+/ay2Ww1/t6rMfsFq8Y7gJpxBAyAz+ratatmzpyp0aNHa8+ePY52V482mMFmsznWV+W9997TiRMnnNpqqykkJEQ33HCD9u7dq8TExGpH/5KSktShQwe1adNGycnJWrdundNRnOLiYm3YsMFD7652rtYtud5HVvO195CXl6cTJ05UOxqVkZGhPXv2VDv9cPDgwY6bGK9evVqxsbEuXSUxMzNTgwcP1s0336wXX3xRTz/9dI1X3HRVmzZtdN1112ndunVO47uyslJr1qxRly5d3L7Xl7vL7N69u/Lz852CSFlZmT744AOn5bqyb3j11Vednv/v//6vysvLHVdsdHVdrq6vNrXt3y7Xpk0bJSUl6e2331ZZWZmj/fz58/rLX/7i9jqr+Oo2C7QUHAEDYKkDBw5UuwqiJF111VUKDAzUDTfcoLvvvlt9+vRRaGioMjMzHVcBq9K/f39J0tKlS5WWlqZWrVp57Kpd48eP16pVq9SnTx8lJiZq9+7deuaZZ6r9nqa2mkJDQ7V06VINGzZMw4cP1wMPPKDu3buruLhYn3/+uTZs2KCPP/5YkvTf//3fGjt2rEaPHq2f/vSnqqio0KJFi9SmTRudOXPGI++vLq7W7WofeYMvvYeq32a9+eabio+P19VXX63s7GxHOMrPz9eBAweUkJAg6T8BrLS0VPPnz9dbb71V7zpyc3OVm5vrCHPTpk3TV199pYceekhRUVFu//6uysKFCzV69GjdcMMNmj17tgIDA7V8+XIdOHBAr732WoOOarqzzEmTJunJJ5/U5MmT9bOf/UwlJSX63e9+V+1y+nVth1XWrVungIAAjR492nEVxAEDBujOO+90a12urq9KYWGhS/u3K/3qV7/SrbfeqptvvlkPP/ywKioq9Mwzz6ht27YN3i/48jYLtAhevggIgBairqsg6v+v2FZSUmJMnz7dSExMNMLCwozWrVsbvXv3NubOnWtcuHDBaXlz5swxYmJiDD8/P0OSsWXLljqvgnjq1Cmn+dPS0ow2bdpUq3PkyJFGv379HM/Pnj1r3HvvvUZERIQREhJiDBs2zPj000+NkSNHGiNHjqy3pio5OTnGPffcY3Tu3Nlo1aqV0alTJ2Po0KHGU0895bSMd99910hMTDQCAwONrl27Gk8//XSNV1+rq49duQrilf1R29XOXKnb1T6qbd21aexnZ+V7cOVqcU8++aQREBBgvPnmm0b37t2N4OBgY9SoUcbOnTuNq6++2oiOjna6OuJHH31kDB482FiyZEm1qzvWdqW+d955x5BkfPbZZ07t999/vxEcHGx88sknNdZd31UQDcMwPv30U+PGG2802rRpY7Ru3doYMmSIsWHDBqdp3LkKoqvLrPL+++8bAwcONFq3bm306NHDWLZsWY3bRm3bYdW0u3fvNlJTU422bdsaoaGhxl133WWcPHmyQeuqa31XcmX/Vts4Wr9+vdG/f3+n/cJDDz1ktG/f3mk6V8enq+OdqyACnmEzDMPwfMwDAKBlu+WWW5SXl1fr6WZXOnfunGJiYmS32/XRRx853V9q2rRp2rp1qz7//HOn++i5yjAMVVRUaPXq1br33nuVmZmppKQkt5bR1MybN0/z58/XqVOn3P7Nmq+5ePGiBg4cqM6dO2vTpk3eLgeAmzgFEQAAC+zevVt33HGHy9O3a9dOnTt3VkpKSo039/3yyy/VqlUr9evXTwcOHHCrlnfeecetWuBd9957r0aPHq3o6Gjl5eXphRde0OHDh7V06VJvlwagAQhgAAB42LFjx5Sfn6/k5GSX5zl//rzOnz+v+fPnV3tt3rx5mjlzpiRVu6S9K0aNGuX4TZqkGgMefEdxcbFmz56tU6dOqVWrVrrmmmv0/vvv6zvf+Y63SwPQAJyCCACAD3rkkUckSb/97W+9XAkAwExchh4AAB+SnZ0tu92ugwcP6qmnnvJ2OQAAk3EEDAAAAAAswhEwAAAAALAIAQwAAAAALMJVEBuosrJSX3/9tUJDQ2Wz2bxdDgAAAAAvMQxDxcXFiomJkZ9f3ce4CGAN9PXXXys2NtbbZQAAAADwEcePH1eXLl3qnIYA1kChoaGSLnVyWFiYl6sBAAAA4C1FRUWKjY11ZIS6EMAaqOq0w7CwMAIYAAAAAJd+msRFOAAAAADAIgQwAAAAALAIAQwAAAAALEIAAwAAAACLEMAAAAAAwCIEMAAAAACwCAEMAAAAACxCAAMAAAAAi3AjZgBoISoqDe3KOaP84hJFhAYrOS5c/n713zASAACYhwAGAC3AxgO5mr/hkHILSxxt0fZgzU2N19iEaC9WBgBAy8IpiADQzG08kKsH1uxxCl+SlFdYogfW7NHGA7leqgwAgJaHAAYAzVhFpaH5Gw7JqOG1qrb5Gw6porKmKQAAgNmaRQD75JNPlJqaqpiYGNlsNr399tv1zpORkaHBgwcrODhYPXr00AsvvOD5QgHAYrtyzlQ78nU5Q1JuYYl25ZyxrigAAFqwZhHALly4oAEDBmjZsmUuTZ+Tk6NbbrlFw4cP1969e/X444/roYce0tq1az1cKQBYK7+49vDVkOkAAEDjNIuLcIwbN07jxo1zefoXXnhBXbt21ZIlSyRJffv2VVZWlp599llNnDjRQ1UCgPUiQoNNnQ4AADROszgC5q7t27drzJgxTm0333yzsrKydPHixRrnKS0tVVFRkdMDAHxdcly4ou3Bqu1i8zZduhpicly4lWUBANBitcgAlpeXp8jISKe2yMhIlZeX6/Tp0zXOs3DhQtntdscjNjbWilIBoFH8/WyamxovSdVCWNXzuanx3A8MAACLtMgAJkk2m/OXDcMwamyvMmfOHBUWFjoex48f93iNAGCGsQnRWjHlGkWEBTm1R9mDtWLKNdwHDAAACzWL34C5KyoqSnl5eU5t+fn5CggIUIcOHWqcJygoSEFBQTW+BgC+bmxCtK6/uqP6z9skSVr1w2s1vGcnjnwBAGCxFnkELCUlRZs3b3Zq27Rpk5KSktSqVSsvVQUAnnV52EqOCyd8AQDgBc0igJ0/f17Z2dnKzs6WdOky89nZ2Tp27JikS6cP/uAHP3BMP336dH355ZdKT0/X4cOH9dJLL+mPf/yjZs+e7Y3yAQAAALQQzeIUxKysLN1www2O5+np6ZKktLQ0rVq1Srm5uY4wJklxcXF6//339cgjj+j5559XTEyMfve733EJegAAAAAe1SwC2KhRoxwX0ajJqlWrqrWNHDlSe/bs8WBVAAAAAOCsWZyCCAAAAABNAQEMAAAAACxCAAMAAAAAixDAAAAAAMAiBDAAAAAAsAgBDAAAAAAsQgADAAAAAIsQwAAAAADAIgQwAAAAALAIAQwAAAAALEIAAwAAAACLEMAAAAAAwCIEMAAAAACwCAEMAAAAACxCAAMAAAAAixDAAAAAAMAiBDAAAAAAsAgBDAAAAAAsQgADAAAAAIsQwAAAAADAIgQwAAAAALBIgLcLAAA0PRWVhnblnFF+cYkiQoOVHBcufz+bt8tCM8DYAtDcEcAAAG7ZeCBX8zccUm5hiaMt2h6suanxGpsQ7cXK0NQxtgC0BJyCCABw2cYDuXpgzR6nL8iSlFdYogfW7NHGA7leqgxNHWMLQEtBAAMAuKSi0tD8DYdk1PBaVdv8DYdUUVnTFEDtGFsAWhICGADAJbtyzlQ7OnE5Q1JuYYl25Zyxrig0C4wtAC0JAQwA4JL84tq/IDdkOqAKYwtAS0IAAwC4JCI02NTpgCqMLQAtCQEMAOCS5LhwRduDVdsFwW26dMW65LhwK8tCM8DYAtCSEMAAAC7x97Npbmq8JFX7olz1fG5qPPdsgtsYWwBaEgIYAMBlYxOitWLKNYoIC3Jqj7IHa8WUa7hXExqMsQWgpeBGzAAAt4xNiNb1V3dU/3mbJEmrfnithvfsxNEJNBpjC0BLwBEwAIDbLv9CnBwXzhdkmIaxBaC5I4ABAAAAgEUIYAAAAABgEQIYAAAAAFiEAAYAAAAAFmk2AWz58uWKi4tTcHCwBg8erE8//bTWabdu3SqbzVbt8Y9//MPCigEAAAC0NM0igL3xxhuaNWuWfvGLX2jv3r0aPny4xo0bp2PHjtU535EjR5Sbm+t49OzZ06KKAQAAALREzSKAPffcc7r33nv1ox/9SH379tWSJUsUGxurFStW1DlfRESEoqKiHA9/f3+LKgYAAADQEjX5AFZWVqbdu3drzJgxTu1jxozRtm3b6px30KBBio6O1k033aQtW7bUOW1paamKioqcHgAAAADgjiYfwE6fPq2KigpFRkY6tUdGRiovL6/GeaKjo7Vy5UqtXbtW69atU+/evXXTTTfpk08+qXU9CxculN1udzxiY2NNfR8AAAAAmr8AbxdgFpvN5vTcMIxqbVV69+6t3r17O56npKTo+PHjevbZZzVixIga55kzZ47S09Mdz4uKighhAAAAANzS5I+AdezYUf7+/tWOduXn51c7KlaXIUOG6OjRo7W+HhQUpLCwMKcHAAAAALijyQewwMBADR48WJs3b3Zq37x5s4YOHerycvbu3avo6GizywMAAAAAh2ZxCmJ6erqmTp2qpKQkpaSkaOXKlTp27JimT58u6dLpgydOnNDq1aslSUuWLFH37t3Vr18/lZWVac2aNVq7dq3Wrl3rzbcBAAAAoJlrFgFs0qRJKigo0K9+9Svl5uYqISFB77//vrp16yZJys3NdbonWFlZmWbPnq0TJ06odevW6tevn9577z3dcsst3noLAAAAAFqAZhHAJGnGjBmaMWNGja+tWrXK6fmjjz6qRx991IKqAAAAAOA/mvxvwAAAAACgqSCAAQAAAIBFCGAAAAAAYBECGAAAAABYhAAGAAAAABYhgAEAAACARQhgAAAAAGARAhgAAAAAWIQABgAAAAAWIYABAAAAgEUIYAAAAABgEQIYAAAAAFiEAAYAAAAAFiGAAQAAAIBFCGAAAAAAYBECGAAAAABYJMDbBQAA0JxVVBralXNG+cUliggNVnJcuPz9bN4uCwDgJQQwAAA8ZOOBXM3fcEi5hSWOtmh7sOamxmtsQrQXKwMAeAunIAIA4AEbD+TqgTV7nMKXJOUVluiBNXu08UCulyoDAHgTAQwAAJNVVBqav+GQjBpeq2qbv+GQKiprmgIA0JwRwAAAMNmunDPVjnxdzpCUW1iiXTlnrCsKAOATCGAAAJgsv7j28NWQ6QAAzQcBDAAAk0WEBps6HQCg+SCAAQBgsuS4cEXbg1XbxeZtunQ1xOS4cCvLAgD4AAIYAAAm8/ezaW5qvCRVC2FVz+emxnM/MABogQhgAAB4wNiEaK2Yco0iwoKc2qPswVox5RruAwYALRQ3YgYAwEPGJkTr+qs7qv+8TZKkVT+8VsN7duLIFwC0YBwBAwDAgy4PW8lx4YQvAGjhCGAAAAAAYBG3T0F899133V7J6NGj1bp1a7fnAwAAAIDmxO0Advvtt7s1vc1m09GjR9WjRw93VwUAAAAAzUqDTkHMy8tTZWWlS4+QkBCzawYAAACAJsntAJaWlubW6YRTpkxRWFiYu6sBAAAAgGbH7VMQX375Zce/T548qcjIyDqnX7FihftVAQAAAEAz1KirIE6cOFHl5eU1vlZbO4CWraLS0PZ/Feid7BPa/q8CVVQa3i4JAADAMo26EXP79u31k5/8pNpRroKCAk2cOFFbt25tzOKBJqOi0tCunDPKLy5RRGgw9/qpxcYDuZq/4ZByC0scbdH2YM1NjdfYhGi3l0e/A7CKmfsb9l2uo6/cQ381DY0KYH/605+UnJysP/zhD/rRj34kSTp8+LDGjx+vfv36mVKgr9v17zO6ITGUwd2CmRkqmvOOc+OBXD2wZo+uPN6VV1iiB9bs0Yop17jVX2aHOeByzXlb9ITm3l9m7m/Yd7mOP9q5pyX0V3P5Q0ijAli7du20du1ajRw5Uv3799fZs2c1efJk3X///Vq0aJFZNbpk+fLleuaZZ5Sbm6t+/fppyZIlGj58eK3TZ2RkKD09XQcPHlRMTIweffRRTZ8+3e313vNKpjpH/MsnBrcvbii+zIz+MjNUNOf/lCsqDc3fcKhaP0mSIckmaf6GQxodH+XSZ2B2mPNlbNfWa87boif4cn/54n7ezH2Xr36H8LV+r1qer4YT+sv6mry933I7gN12220aOHCgBg0apIEDB6p///56/vnndeutt6qkpETPP/+80tLSPFFrrd544w3NmjVLy5cv1/XXX6//+Z//0bhx43To0CF17dq12vQ5OTm65ZZbdN9992nNmjX6+9//rhkzZqhTp06aOHGi2+v3hcHtixuKJ5ZlFjP6y8xQ4YlA4Uuf4a6cM059fSVDUm5hiXblnFHKVR3qrcXMMOfLvP0fhFV8aR/h6+Hel/pK8u3+8rX9vCf+EOWL3yF8rd+ravLVcEJ/Wf/HC1/Yb7kdwHr27Km///3vWr58uQoKCtSuXTsNGDBAhmHo+9//vgYOHKiLFy+qVatWnqi3Rs8995zuvfdex2mQS5Ys0QcffKAVK1Zo4cKF1aZ/4YUX1LVrVy1ZskSS1LdvX2VlZenZZ591O4AFlpfK399fNkkL1+3VTd3DXBrcmw/ladbr2TIkBV3WfragVLNWbdfSyQM1Oj7KpRrMXFbV8ha8/w/lFf1nQ4kKC9bjt/RxazlmL0u6tJPJ+vKsThWXqFNosJK6tXf7y4dZ/bUr54zOFBQ6LeNKZwpKtevwCV0XF17rNBWVhhau26vA8tIaX3d3bEm+9xmeOnVWQbW8vyunq4yu+zYXZvX7lcwYW2Yye7uWpMqycsfnUPnNN6osb/hJEGYty+x9RGN4YluUmmdfSb7dX762nzd7Wb76HcIX+93scUp/Wd9fZtZU27JK/QNl2GyW/RHXZhhGgy9B9tVXXyk7O9vpkZOTo4CAAPXp00f79u0zs9YalZWVKSQkRG+++abuuOMOR/vDDz+s7OxsZWRkVJtnxIgRGjRokJYuXepoW79+ve6880598803NYbH0tJSlZb+58MqKipSbGysdl3dU239/U1+VwAAAACscPv4X6s04D8R8bX7htR7Rs6VioqKZLfbVVhYWO89kBv1G7AuXbqoS5cuGj9+vKPt/Pnz2rt3rz777LPGLNplp0+fVkVFRbX7kUVGRiovL6/GefLy8mqcvry8XKdPn1Z0dPXDjgsXLtT8+fPNKxwAAACAz8kvrv1nE2ZwO4B99tlnSkhIkJ9fzbcQa9u2rYYPH+64AMbBgwfVu3dvBQQ0KuvVy2ZzPkxoGEa1tvqmr6m9ypw5c5Senu54XnUE7K6xT8o/KMTRvuqHyfUe3v3LZ1/rZ2/VH1Cf+W6ixifGWLasnTlnNO3lXfUuy5X3aOayKioNfee5DKfTbi5nkxQZFqwP00fWe7jYzP6qqutkUUmN5127Wpev9ruZy5IunYbw8OvZkuTUX1U94+qpFmb1++XLMmNsVdl8KE9PvXdY+cX/OWLu7mliZo5TX2X2+JKkb8rKNfipDyVJu3/5HYUEuvf/jidqMoOn6mqu/eWL+3kzl+Wr3yF8td/NHKf0139Y2V9WfLcp9Q90eh4RGlzv+hrD7VQ0aNAg5eXlqVOnTi5Nn5KSouzsbPXo0cPt4lzRsWNH+fv7VzvalZ+fX+0oV5WoqKgapw8ICFCHDjUfbgwKClJQUPWza8sCguQXECSbpCh7sJL7dpZfPYO7U6f2Toc565rOLySk3mnMWlb+xbMuLSv/os3SZe38V4G+/MaQ6ljel98YyjpZUu/hYjP7y0/SnAmD9MCaPZJqDhVzJgxSq7Zt6lxOct/WCu9gV15h7TtOV8eWr36GknRzUg8tCW5d6w9xb3bxB69m9btk7tiS/v+HvW8dvlTTZcs89o2h6W8d1ooprV36Ya+Z49RXmT2+JMkvoNyxTL+QEPm5GSjM3BbN5Im+kppvf/nift7MZfnqdwhf7Xczxyn95Z3+MrMml5fl4T8auR3ADMPQE088oRAXd/JlZWVuF+WOwMBADR48WJs3b3b6DdjmzZt122231ThPSkqKNmzY4NS2adMmJSUlNejiIVUf9dzUeJf+Sp4cF65oe7ApH76Zy3I17bsynZnLcvUwsCvTmdlfkjQ2IVorplxTLVREuXF1H38/m+amxuuBNXtkU807TlfHlq9+hlXGJkRrdHxUo6/iZka/S+aOLTOvQmX2OPVFnhhfjWXmtmgmX+wryXf7yxf382Yuy1e/Q/hqv5s5TumvS6zuLzNr8pX9ltsBbMSIETpy5IjL06ekpKh167qvatZY6enpmjp1qpKSkpSSkqKVK1fq2LFjjvt6zZkzRydOnNDq1aslSdOnT9eyZcuUnp6u++67T9u3b9cf//hHvfbaaw1avzcHty9uKGYvy8wvH57Y8MwIFWbtOH31M7ycv5/N7R+21sSMfjdzbJl5qX1f+Q/Ck3w1ZJr5ZdssvtpXkm/2l6/u581alq9+h/DlfvfFcEJ/XeLqe/S1P4Q0VqOuguhLli9frsWLFys3N1cJCQn67W9/qxEjRkiSpk2bpi+++EJbt251TJ+RkaFHHnnEcSPmxx57zK0bMVdd6WTz3hzdkNitQTthX7uPRNVy6jqE3ZD7LDR2WRWVhoYt+rjeLx9/e+xGlz8HX72/kpk3Y5R85zP0VWaOrXeyTzh+41aXpZMH6raBnV2qz1fHqVnMHl/flJUr/skPJEmHfnWz279pupyv3m9LMm9bbM79JbWM7cfXvkOYvSyzmTVO6S/3mPkefekep1dy5yqIjQpgx48fV2xsbENnb9Lc6eS6+OJA8sUdiye+fPjiFwaz+OJn6KvMGlvb/1Wgu17cUe907l7atjmPU8nc8WVmoPBFZm+Lzb2/pOa//fjidwizl+Wr6C/3tIT36LEA9tprr+muu+5yPPfz81P79u01YMAADRgwQAMHDtSAAQNUWlqq559/3nHKX3NkVgDzVb64Y2nuQcBsvvgZ+iozxpYnjtS2FGaNLwKFe1pCfwGAVUy/D1heXp5mzJihdu3aOQWwf//7344bMO/du1dvvfWWvv76a0lqlqGkJTHrdzpmLsvMc/BbAl/8DH2VGWOrJfxuy1Oa+/gyE30FAE2fSwFs5cqVKi8v10svveTU3r17d3Xv3l233367o2379u1KS0vTokWLTC0UkPjyAc8xY2z5wg97AQCAb3MpgD300EN66KGHNHHiRK1du7bOaVNSUrR06VL98pe/dLosPAC0BBypBQAAdXEpgLVr106rV6+udu+sixcv1njfrJ49e+rgwYPmVAgATQxHagEAQG3c+sVtamqq0/M2bdooPj5egwYN0sCBAzVo0CDFxMTo97//vcaMGWNqoQAAAADQ1DXqkkcff/yx9u3bp3379unVV1/V448/rm+//VaSNGbMGP3iF79QYmKiEhMT1bdvX1MKBgAAAICmqlEBbNiwYRo2bJjjeWVlpY4cOeK4MuLu3bv10ksvKT8/XxUVFY0uFgAAAACaMlNv+uHn56e+ffuqb9++TperP3nypJmrAQAAAIAmyc+KlURGRlqxGgAAAADwaZYEMAAAAAAAAQwAAAAALEMAAwCgBaqoNBz/3pVzxuk5AMBzCGAAgGaBQOG6jQdy9Z3nMhzPp72cqWGLPtbGA7lerAoAWgYCGACgySNQuG7jgVw9sGaPThaVOrXnFZbogTV76DMA8DACGACgSSNQuK6i0tD8DYdU07HBqrb5Gw5x9BAAPIgABgBosggU7tmVc0a5hSW1vm5Iyi0s0a6cM9YVBQAtDAEMANBkESjck19ce181ZDoAgPsIYACAJotA4Z6I0GBTpwMAuI8ABgBosggU7kmOC1e0PVi2Wl63SYq2Bys5LtzKsgCgRSGAAQCaLAKFe/z9bJqbGi9J1fqs6vnc1Hj5+9XWowCAxiKAAQCaLAKF+8YmRGvFlGsUZXc+KhhlD9aKKddobEK0lyoDgJbBZhgGl4ZqgKKiItntdhUWFiosLMzb5QBAi7bxQK7mbzjkdEGOaHuw5qbGEyhqUVFpaFfOGeUXlygi9NJRQoIqADSMO9mAANZABDAA8C0ECgCAt7iTDQIsqgkAAI/y97Mp5aoO3i4DAIA68RswAAAAALAIAQwAAAAALEIAAwAAAACLEMAAAAAAwCIEMAAAAACwCAEMAAAAACxCAAMAAAAAixDAAAAAAMAiBDAAAAAAsAgBDAAAAAAsQgADAAAAAIsQwAAAAADAIk0+gJ09e1ZTp06V3W6X3W7X1KlTde7cuTrnmTZtmmw2m9NjyJAh1hQMAAAAoMUK8HYBjXX33Xfrq6++0saNGyVJ999/v6ZOnaoNGzbUOd/YsWP18ssvO54HBgZ6tE4AAAAAaNIB7PDhw9q4caN27Nih6667TpL04osvKiUlRUeOHFHv3r1rnTcoKEhRUVFWlQoAAAAATfsUxO3bt8tutzvClyQNGTJEdrtd27Ztq3PerVu3KiIiQr169dJ9992n/Pz8OqcvLS1VUVGR0wMAAAAA3NGkA1heXp4iIiKqtUdERCgvL6/W+caNG6dXX31VH3/8sX7zm98oMzNTN954o0pLS2udZ+HChY7fmdntdsXGxpryHgAAAAC0HD4ZwObNm1ftIhlXPrKysiRJNput2vyGYdTYXmXSpEm69dZblZCQoNTUVP31r3/VP//5T7333nu1zjNnzhwVFhY6HsePH2/8GwUAAADQovjkb8BmzpypyZMn1zlN9+7d9dlnn+nkyZPVXjt16pQiIyNdXl90dLS6deumo0eP1jpNUFCQgoKCXF4mAAAAAFzJJwNYx44d1bFjx3qnS0lJUWFhoXbt2qXk5GRJ0s6dO1VYWKihQ4e6vL6CggIdP35c0dHRDa4ZAAAAAOrjk6cguqpv374aO3as7rvvPu3YsUM7duzQfffdp/HjxztdAbFPnz5av369JOn8+fOaPXu2tm/fri+++EJbt25VamqqOnbsqDvuuMNbbwUAAABAC9CkA5gkvfrqq+rfv7/GjBmjMWPGKDExUX/605+cpjly5IgKCwslSf7+/tq/f79uu+029erVS2lpaerVq5e2b9+u0NBQb7wFAAAAAC2EzTAMw9tFNEVFRUWy2+0qLCxUWFiYt8sBAAAA4CXuZIMmfwQMAAAAAJoKAhgAAAAAWIQABgAAAAAWIYABAAAAgEUIYAAAAABgEQIYAAAAAFiEAAYAAAAAFiGAAQAAAIBFCGAAAAAAYBECGAAAAABYhAAGAAAAABYhgAEAAACARQhgAAAAAGARAhgAAAAAWIQABgAAAAAWIYABAAAAgEUIYAAAAABgEQIYAAAAAFiEAAYAAAAAFiGAAQAAAIBFCGAAAAAAYBECGAAAAABYhAAGAAAAABYhgAEAAACARQhgAAAAAGARAhgAAAAAWIQABgAAAAAWIYABAAAAgEUIYAAAAABgEQIYAAAAAFiEAAYAAAAAFiGAAQAAAIBFCGAAAAAAYBECGAAAAABYhAAGAAAAABYhgAEAAACARQhgAAAAAGARAhgAAAAAWKTJB7Bf//rXGjp0qEJCQtSuXTuX5jEMQ/PmzVNMTIxat26tUaNG6eDBg54tFAAAAECL1+QDWFlZmb73ve/pgQcecHmexYsX67nnntOyZcuUmZmpqKgojR49WsXFxR6sFAAAAEBLZzMMw/B2EWZYtWqVZs2apXPnztU5nWEYiomJ0axZs/TYY49JkkpLSxUZGalFixbpxz/+cY3zlZaWqrS01PG8qKhIsbGxKiwsVFhYmGnvAwAAAEDTUlRUJLvd7lI2aPJHwNyVk5OjvLw8jRkzxtEWFBSkkSNHatu2bbXOt3DhQtntdscjNjbWinIBAAAANCMtLoDl5eVJkiIjI53aIyMjHa/VZM6cOSosLHQ8jh8/7tE6AQAAADQ/PhnA5s2bJ5vNVucjKyurUeuw2WxOzw3DqNZ2uaCgIIWFhTk9AAAAAMAdAd4uoCYzZ87U5MmT65yme/fuDVp2VFSUpEtHwqKjox3t+fn51Y6KAQAAAICZfDKAdezYUR07dvTIsuPi4hQVFaXNmzdr0KBBki5dSTEjI0OLFi3yyDoBAAAAQPLRUxDdcezYMWVnZ+vYsWOqqKhQdna2srOzdf78ecc0ffr00fr16yVdOvVw1qxZWrBggdavX68DBw5o2rRpCgkJ0d133+2ttwEAAACgBfDJI2DuePLJJ/XKK684nlcd1dqyZYtGjRolSTpy5IgKCwsd0zz66KP69ttvNWPGDJ09e1bXXXedNm3apNDQUEtrBwAAANCyNJv7gFnNnWv9AwAAAGi+uA8YAAAAAPggAhgAAAAAWIQABgAAAAAWIYABAAAAgEUIYAAAAABgEQIYAAAAAFiEAAYAAAAAFiGAAQAAAIBFCGAAAAAAYBECGAAAAABYhAAGAAAAABYJ8HYBTZVhGJKkoqIiL1cCAAAAwJuqMkFVRqgLAayBCgoKJEmxsbFergQAAACALyguLpbdbq9zGgJYA4WHh0uSjh07Vm8nw1xFRUWKjY3V8ePHFRYW5u1yWgz63Xvoe++g372HvvcO+t176HvvMLPfDcNQcXGxYmJi6p2WANZAfn6Xfj5nt9vZULwkLCyMvvcC+t176HvvoN+9h773Dvrde+h77zCr3109KMNFOAAAAADAIgQwAAAAALAIAayBgoKCNHfuXAUFBXm7lBaHvvcO+t176HvvoN+9h773Dvrde+h77/BWv9sMV66VCAAAAABoNI6AAQAAAIBFCGAAAAAAYBECGAAAAABYhAAGAAAAABYhgDXQ8uXLFRcXp+DgYA0ePFiffvqpt0tq1ubNmyebzeb0iIqK8nZZzdInn3yi1NRUxcTEyGaz6e2333Z63TAMzZs3TzExMWrdurVGjRqlgwcPeqfYZqS+fp82bVq1bWDIkCHeKbYZWbhwoa699lqFhoYqIiJCt99+u44cOeI0DWPeM1zpe8a9Z6xYsUKJiYmOm8+mpKTor3/9q+N1xrxn1NfvjHdrLFy4UDabTbNmzXK0WT3mCWAN8MYbb2jWrFn6xS9+ob1792r48OEaN26cjh075u3SmrV+/fopNzfX8di/f7+3S2qWLly4oAEDBmjZsmU1vr548WI999xzWrZsmTIzMxUVFaXRo0eruLjY4kqbl/r6XZLGjh3rtA28//77FlbYPGVkZOjBBx/Ujh07tHnzZpWXl2vMmDG6cOGCYxrGvGe40vcS494TunTpoqefflpZWVnKysrSjTfeqNtuu83xhZMx7xn19bvEePe0zMxMrVy5UomJiU7tlo95A25LTk42pk+f7tTWp08f4+c//7mXKmr+5s6dawwYMMDbZbQ4koz169c7nldWVhpRUVHG008/7WgrKSkx7Ha78cILL3ihwubpyn43DMNIS0szbrvtNq/U05Lk5+cbkoyMjAzDMBjzVrqy7w2DcW+l9u3bG3/4wx8Y8xar6nfDYLx7WnFxsdGzZ09j8+bNxsiRI42HH37YMAzv7Oc5AuamsrIy7d69W2PGjHFqHzNmjLZt2+alqlqGo0ePKiYmRnFxcZo8ebL+/e9/e7ukFicnJ0d5eXlO4z8oKEgjR45k/Ftg69atioiIUK9evXTfffcpPz/f2yU1O4WFhZKk8PBwSYx5K13Z91UY955VUVGh119/XRcuXFBKSgpj3iJX9nsVxrvnPPjgg7r11lv1ne98x6ndG2M+wCNLbcZOnz6tiooKRUZGOrVHRkYqLy/PS1U1f9ddd51Wr16tXr166eTJk3rqqac0dOhQHTx4UB06dPB2eS1G1Rivafx/+eWX3iipxRg3bpy+973vqVu3bsrJydETTzyhG2+8Ubt371ZQUJC3y2sWDMNQenq6hg0bpoSEBEmMeavU1PcS496T9u/fr5SUFJWUlKht27Zav3694uPjHV84GfOeUVu/S4x3T3r99de1Z88eZWZmVnvNG/t5AlgD2Ww2p+eGYVRrg3nGjRvn+Hf//v2VkpKiq666Sq+88orS09O9WFnLxPi33qRJkxz/TkhIUFJSkrp166b33ntPEyZM8GJlzcfMmTP12Wef6W9/+1u11xjznlVb3zPuPad3797Kzs7WuXPntHbtWqWlpSkjI8PxOmPeM2rr9/j4eMa7hxw/flwPP/ywNm3apODg4Fqns3LMcwqimzp27Ch/f/9qR7vy8/OrJWd4Tps2bdS/f38dPXrU26W0KFVXnmT8e190dLS6devGNmCSn/zkJ3r33Xe1ZcsWdenSxdHOmPe82vq+Jox78wQGBurqq69WUlKSFi5cqAEDBmjp0qWMeQ+rrd9rwng3x+7du5Wfn6/BgwcrICBAAQEBysjI0O9+9zsFBAQ4xrWVY54A5qbAwEANHjxYmzdvdmrfvHmzhg4d6qWqWp7S0lIdPnxY0dHR3i6lRYmLi1NUVJTT+C8rK1NGRgbj32IFBQU6fvw420AjGYahmTNnat26dfr4448VFxfn9Dpj3nPq6/uaMO49xzAMlZaWMuYtVtXvNWG8m+Omm27S/v37lZ2d7XgkJSXp+9//vrKzs9WjRw/LxzynIDZAenq6pk6dqqSkJKWkpGjlypU6duyYpk+f7u3Smq3Zs2crNTVVXbt2VX5+vp566ikVFRUpLS3N26U1O+fPn9fnn3/ueJ6Tk6Ps7GyFh4era9eumjVrlhYsWKCePXuqZ8+eWrBggUJCQnT33Xd7seqmr65+Dw8P17x58zRx4kRFR0friy++0OOPP66OHTvqjjvu8GLVTd+DDz6oP//5z3rnnXcUGhrq+Auo3W5X69atHfeKYcybr76+P3/+POPeQx5//HGNGzdOsbGxKi4u1uuvv66tW7dq48aNjHkPqqvfGe+eExoa6vTbUunSmVQdOnRwtFs+5j1ybcUW4Pnnnze6detmBAYGGtdcc43TZXNhvkmTJhnR0dFGq1atjJiYGGPChAnGwYMHvV1Ws7RlyxZDUrVHWlqaYRiXLtc6d+5cIyoqyggKCjJGjBhh7N+/37tFNwN19fs333xjjBkzxujUqZPRqlUro2vXrkZaWppx7Ngxb5fd5NXU55KMl19+2TENY94z6ut7xr3n3HPPPY7vMJ06dTJuuukmY9OmTY7XGfOeUVe/M96tdfll6A3D+jFvMwzD8Ey0AwAAAABcjt+AAQAAAIBFCGAAAAAAYBECGAAAAABYhAAGAAAAABYhgAEAAACARQhgAAAAAGARAhgAAAAAWIQABgAAAAAWIYABAAAAgEUIYAAANNBPf/pTpaamersMAEATQgADAKCBsrOzNXDgQG+XAQBoQghgAAA00L59+zRo0CBvlwEAaEIIYAAANMDx48dVUFDgOAJ27tw5paamaujQocrNzfVucQAAn0UAAwCgAbKzs2W32xUXF6f9+/fr2muvVXR0tLZu3aro6GhvlwcA8FEEMAAAGiA7O1sDBgzQa6+9phEjRmj27NlauXKlAgMDvV0aAMCH2QzDMLxdBAAATc3EiRO1ZcsWSdJf/vIXDR061MsVAQCaAo6AAQDQANnZ2Zo4caJKSkp07tw5b5cDAGgiOAIGAICbiouLZbfbtXv3bu3bt08PP/ywtm3bpn79+nm7NACAjwvwdgEAADQ12dnZ8vf3V3x8vAYNGqSDBw8qNTVVu3btUseOHb1dHgDAh3EKIgAAbtq3b5/69OmjoKAgSdKiRYsUHx+vCRMmqKyszMvVAQB8GacgAgAAAIBFOAIGAAAAABYhgAEAAACARQhgAAAAAGARAhgAAAAAWIQABgAAAAAWIYABAAAAgEUIYAAAAABgEQIYAAAAAFiEAAYAAAAAFiGAAQAAAIBFCGAAAAAAYJH/A0Hd7VybrH/GAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 1000x300 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "L = 32  # number of random samples\n",
    "N = 10000  # number of sample functions\n",
    "\n",
    "# generate input signal (white Gaussian noise)\n",
    "np.random.seed(2)\n",
    "x = np.random.normal(size=(N, L))\n",
    "x[:, L // 2] += 1\n",
    "# generate output signal\n",
    "h = 2 * np.fft.irfft([1, 1, 1, 0, 0, 0])\n",
    "y = np.asarray([np.convolve(x[n, :], h, mode=\"full\") for n in range(N)])\n",
    "\n",
    "\n",
    "def estimate_plot_linear_mean(x):\n",
    "    \"\"\"Estimate and plot linear mean.\"\"\"\n",
    "    # estimate linear mean by ensemble average\n",
    "    mu = 1 / N * np.sum(x, 0)\n",
    "    # plot linear mean\n",
    "    plt.stem(mu)\n",
    "    plt.xlabel(r\"$k$\")\n",
    "    plt.ylabel(r\"$\\hat{\\mu}[k]$\")\n",
    "    plt.axis([0, x.shape[1], -1.2, 1.2])\n",
    "\n",
    "\n",
    "plt.figure(figsize=(10, 3))\n",
    "plt.title(r\"Estimated linear mean $\\hat{\\mu}_x[k]$ of input signal\")\n",
    "estimate_plot_linear_mean(x)\n",
    "\n",
    "plt.figure(figsize=(10, 3))\n",
    "plt.title(r\"Estimated linear mean $\\hat{\\mu}_y[k]$ of output signal\")\n",
    "estimate_plot_linear_mean(y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise**\n",
    "\n",
    "* Can you estimate the impulse response $h[k]$ of the system from above plots of $\\hat{\\mu}_x[k]$ and $\\hat{\\mu}_y[k]$?\n",
    "* You can check your results by plotting the impulse response $h[k]$, for instance with the command `plt.stem(h)`.\n",
    "\n",
    "Solution: Inspecting above plot, the linear mean of the input signal can be approximated as $\\mu_x[k] = \\delta[k]$. The linear mean of the output is then given as $\\mu_y[k] = \\delta[k] * h[k] = h[k]$. It follows that the impulse response of the LTI system can be estimated from the linear mean $\\mu_y[k]$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Stationary Input Signal\n",
    "\n",
    "For a (wide-sense) stationary process, the linear mean of the input signal $\\mu_x[k] = \\mu_x$ does not depend on the time index $k$. For a (wide-sense) stationary input signal, also the output signal of the system is (wide-sense) stationary. Using the result for the non-stationary case above yields\n",
    "\n",
    "\\begin{equation}\n",
    "\\begin{split}\n",
    "\\mu_y &= \\mu_x * h[k] \\\\\n",
    "&= \\sum_{\\kappa = -\\infty}^{\\infty}\\mu_x[k-\\kappa]h[\\kappa] \\\\\n",
    "&= \\mu_x \\cdot \\sum_{\\kappa = -\\infty}^{\\infty}h[\\kappa] \\\\\n",
    "&= \\mu_x \\cdot \\sum_{\\kappa = -\\infty}^{\\infty}h[\\kappa]\\cdot\\mathrm{e}^{-\\mathrm{j}\\Omega\\kappa} \\hspace{5mm} \\text{for}\\,\\,\\Omega=0 \\\\\n",
    "&= \\mu_x \\cdot H(\\mathrm{e}^{\\,\\mathrm{j}\\,\\Omega}) \\big\\vert_{\\Omega = 0}\n",
    "\\end{split}\n",
    "\\end{equation}\n",
    "\n",
    "where $H(\\mathrm{e}^{\\,\\mathrm{j}\\,\\Omega}) = \\mathcal{F}_* \\{ h[k] \\}$ denotes the discrete time Fourier transformation (DTFT) of the impulse response. The linear mean of a (wide-sense) stationary input signal is weighted by the transmission characteristics for the constant (i.e. DC, $\\Omega = 0$) component of the LTI system. This implies that for a system which just attenuates the input signal $y[k] = A \\cdot x[k]$, e.g. an ideal amplifier, the linear mean at the output is given as $\\mu_y = A \\cdot \\mu_x$. Furthermore, if the input signal is zero-mean $\\mu_x = 0$, the output signal is also zero-mean $\\mu_y = 0$."
   ]
  },
  {
   "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."
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.9.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}