random_signals/power_spectral_densities.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "nbsphinx": "hidden"
   },
   "source": [
    "# Random Signals\n",
    "\n",
    "*This jupyter notebook is part of a [collection of notebooks](../index.ipynb) on various topics of Digital Signal Processing. Please direct questions and suggestions to [Sascha.Spors@uni-rostock.de](mailto:Sascha.Spors@uni-rostock.de).*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Power Spectral Density\n",
    "\n",
    "The (auto-) [power spectral density](https://en.wikipedia.org/wiki/Spectral_density#Power_spectral_density) (PSD) is defined as the Fourier transformation of the [auto-correlation function](correlation_functions.ipynb) (ACF)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Definition\n",
    "\n",
    "For a continuous-amplitude real-valued wide-sense stationary (WSS) random signal $x[k]$ the PSD is given as\n",
    "\n",
    "\\begin{equation}\n",
    "\\Phi_{xx}(\\mathrm{e}^{\\,\\mathrm{j}\\,\\Omega}) = \\mathcal{F}_* \\{ \\varphi_{xx}[\\kappa] \\}\n",
    "\\end{equation}\n",
    "\n",
    "where $\\mathcal{F}_* \\{ \\cdot \\}$ denotes the [discrete-time Fourier transformation](https://en.wikipedia.org/wiki/Discrete-time_Fourier_transform) (DTFT) and $\\varphi_{xx}[\\kappa]$ the ACF of $x[k]$. Note, the DTFT is performed with respect to $\\kappa$. The ACF of a random signal of finite length $N$ can be expressed by way of a linear convolution\n",
    "\n",
    "\\begin{equation}\n",
    "\\varphi_{xx}[\\kappa] = \\frac{1}{N} \\cdot x_N[k] * x_N[-k]\n",
    "\\end{equation}\n",
    "\n",
    "Taking the DTFT of the left- and right-hand side results in\n",
    "\n",
    "\\begin{equation}\n",
    "\\Phi_{xx}(\\mathrm{e}^{\\,\\mathrm{j}\\,\\Omega}) = \\frac{1}{N} \\, X_N(\\mathrm{e}^{\\,\\mathrm{j}\\,\\Omega})\\, X_N(\\mathrm{e}^{-\\,\\mathrm{j}\\,\\Omega}) = \n",
    "\\frac{1}{N} \\, | X_N(\\mathrm{e}^{\\,\\mathrm{j}\\,\\Omega}) |^2\n",
    "\\end{equation}\n",
    "\n",
    "The last equality results from the definition of the magnitude and the symmetry of the DTFT for real-valued signals. The spectrum $X_N(\\mathrm{e}^{\\,\\mathrm{j}\\,\\Omega})$ quantifies the amplitude density of the signal $x_N[k]$. It can be concluded from above result that the PSD quantifies the squared amplitude or power density of a random signal. This explains the term power spectral density."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Properties\n",
    "\n",
    "The properties of the PSD can be deduced from the properties of the ACF and the DTFT as\n",
    "\n",
    "1. From the link between the PSD $\\Phi_{xx}(\\mathrm{e}^{\\,\\mathrm{j}\\,\\Omega})$ and the spectrum $X_N(\\mathrm{e}^{\\,\\mathrm{j}\\,\\Omega})$ derived above it can be concluded that the PSD real valued\n",
    "\n",
    "    $$\\Phi_{xx}(\\mathrm{e}^{\\,\\mathrm{j}\\,\\Omega}) \\in \\mathbb{R}$$\n",
    "\n",
    "2. From the even symmetry $\\varphi_{xx}[\\kappa] = \\varphi_{xx}[-\\kappa]$ of the ACF it follows that\n",
    "\n",
    "    $$ \\Phi_{xx}(\\mathrm{e}^{\\,\\mathrm{j} \\, \\Omega}) = \\Phi_{xx}(\\mathrm{e}^{\\,-\\mathrm{j}\\, \\Omega}) $$\n",
    "\n",
    "3. The PSD of an uncorrelated random signal is given as\n",
    "\n",
    "    $$ \\Phi_{xx}(\\mathrm{e}^{\\,\\mathrm{j} \\, \\Omega}) = \\sigma_x^2 + \\mu_x^2 \\cdot {\\bot \\!\\! \\bot \\!\\! \\bot}\\left( \\frac{\\Omega}{2 \\pi} \\right) $$\n",
    "    \n",
    "    which can be deduced from the [ACF of an uncorrelated signal](correlation_functions.ipynb#Properties).\n",
    "\n",
    "4. The quadratic mean of a random signal is given as\n",
    "\n",
    "    $$ E\\{ x[k]^2 \\} = \\varphi_{xx}[0] = \\frac{1}{2\\pi} \\int\\limits_{-\\pi}^{\\pi} \\Phi_{xx}(\\mathrm{e}^{\\,\\mathrm{j}\\, \\Omega}) \\,\\mathrm{d} \\Omega $$\n",
    "\n",
    "    The last relation can be found by expressing the ACF by the inverse DTFT."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example - Power Spectral Density of a Speech Signal\n",
    "\n",
    "In this example the PSD $\\Phi_{xx}(\\mathrm{e}^{\\,\\mathrm{j} \\,\\Omega})$ of a speech signal $x[k]$ is estimated by applying a discrete Fourier transformation (DFT) to its ACF. For a better interpretation of the PSD, the frequency axis $f = \\frac{\\Omega}{2 \\pi} \\cdot f_s$ has been chosen for illustration, where $f_s$ denotes the sampling frequency of the signal. \n",
    "\n",
    "In Python the ACF is stored in a vector with indexes $0, 1, ..., 2N -1$ where the indexes correspond to the lags $\\kappa = -N+1,-N+2,....,N-1$. When computing the discrete Fourier transform (DFT) of the ACF numerically by the fast Fourier transform (FFT) one has to take this shift into account. For instance, by multiplying the DFT $\\Phi_{xx}[\\mu]$ by $e^{j \\mu \\frac{2 \\pi}{2N - 1} (N-1)}$ where $N$ denotes the length of the signal $N$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAnUAAAHyCAYAAACNqjipAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvFvnyVgAAIABJREFUeJzs3Xl8XHW9//H3h7Iqu1SURYuC1YoCsul1C1dUcIHrdQNXfurlqhevy3WpGyKKCii4AEIRZN8qO91ooWmBbrR039M16ZI0bdMkzZ58fn/MSTpJJsnMZGbOmTOv5+PBg5lzzpzzmW+mmXe+53y/x9xdAAAAKG77hV0AAAAAho9QBwAAEAOEOgAAgBgg1AEAAMQAoQ4AACAGCHUAAAAxQKgDYs7M3m9mq8OuIxUzKzOzqrDrKDVmttHMzk9z28vM7MU813Ormf0yn8cASgGhDoio4Iu32cwak/67KY3XuZmd3P3c3V9w99F5qvEuM/ttPvaN/sxsVPDz3T/sWnLJ3b/p7r+RCPrAcMTqFwMQQ59092lhFwHJzPZ3946w6xhKsdQJIPfoqQOKkJmdbGYzzGyPmdWa2cPB8pnBJouDnr3P9+35CHoAf2RmS8xsr5ndYWbHmtkkM2sws2lmdlTS9uPNbHtwrJlm9vZg+eWSvijpx8Gxng6WH2dmj5rZDjPbYGb/m7SvQ4Levd1mtkLS2UO8Tzez/zWz9cH7vN7M9gvW7WdmvzCzTWZWY2b3mNkRwbq7zez/gsfHB/v5n+D5m81sV9J+PmFmi8yszsxmmdk7+7TVT8xsiaS9fXvILOHG4Pj1ZrbUzE4N1t0VnFacGrTrDDN7Y9Jr3xqs22Vmq83sc33a6U/Be9tjZi+a2SGSun++dUGbvyc4PfpSUMdOSVcF7/F5M9sZtNv9ZnbkYG2ddOzXmNlTwfuZJ+nNfdYPVvddZnazmU0I3vNcM3tzmm31WzN7taRJko6zfb3Tx5lZk5m9Juk47wo+Xwek856AUkGoA4rTbyQ9K+koSSdI+pskufsHgvWnufuh7v7wAK//tKQPS3qLpE8q8UX6M0kjlfi98L9J206SdIqk10p6RdL9wbHGBY+vC471ySAoPS1psaTjJX1I0vfM7KPBvn6lREh4s6SPSvpqGu/1U5LOkvQuSRdL+lqw/LLgv/MkvUnSoZK6T0/PkFQWPP6gpPWSPpD0/AV37zKzMyTdKem/Jb1G0m2SnjKzg5KOf6mkj0s6MkUP2EeC/b5F0hGSPidpZ9L6LyrxszpG0iIFbReEl6mSHlCiXS+RdIuZjQle90dJZ0r6N0lHS/qxpK6k93Bk0Oazg+fnBu/xWEnXSDJJv5d0nKS3STpR0lVKz82SWiS9Xom27m7vdOpWsOzXSnw2K4J60mkrufteSRdK2hq8v0Pdfauk8mD7bl+W9JC7t6f5noCSEMtQZ2Z3Bn8NLktj2zea2XOW6LUoN7MTClEjkKYngh6k7v/+K1jeLumNko5z9xZ3z/RC9r+5e7W7b5H0gqS57r7Q3VskPS7pjO4N3f1Od29w91YlgsFp3T1iKZwtaaS7X+3ube6+XtLtSnzRS4kv5mvcfZe7V0r6axq1Xhtsv1nSn5UIWVIiMN3g7uvdvVHSTyVdEvSmzZD0viBkfkDSdZLeG7zug8F6Sbpc0m3uPtfdO939bkmtkt6ddPy/unuluzenqK1d0mGS3irJ3H2lu29LWj/B3WcGbfdzSe8xsxMlfULSRnf/p7t3uPtCSY9K+mxQ89ckfdfdtwR1zQr2MZCt7v63YF/N7l7h7lPdvdXdd0i6IXjfgzKzEUoE/ivdfa+7L5N0d9ImA9adtM3j7j4vCMD3Szo9zbYazN2SvpRU46WS7k3ztUDJiGWok3SXpAvS3PaPku5x93dKulqJv26BqPgPdz8y6b/bg+U/VqI3Zp6ZLTezrw2yj1Sqkx43p3h+qJT4AjWzP5jZOjOrl7Qx2OaYAfb7RiVOnfUEUSV6AI8N1h8nqTJp+01p1Np3++OS9rWpz7r9JR3r7usk7VUiULxf0jOStprZaPUOdW+U9H996j0x6Rh9j9+Luz+vRO/gzZJqzGycmR2e6rVB8NwV7PuNks7tc9wvSnqdEm17sKR1gzdLL71qtMTp9IfMbEvwc7tPA//Mko1Uog0H+hkNVne37UmPmxR8ltJoq8E8KWmMmZ2kRA/zHnefl+ZrgZIRy1Dn7jOV+OXZI7jGZLKZLTCzF8zsrcGqMZKeDx5PV+L0DhBp7r7d3f/L3Y9T4tThLZY04jWHvqDEv4nzlThlNipYbt2l9Nm+UtKGPkH0MHf/WLB+mxKhqdsb0qih7/Zbg8dblQgZyes6tC+gzpD0GUkHBj2SM5Q43XuUEqdCu+u9pk+9r3L3B5P22/c99uLuf3X3M5X4XfIWST9KVbuZHarEqdStwXFn9Dnuoe7+LUm1Spz+7HUt2xC19F3+u2DZO9z9cCV6uazfq/rboUQbDvQzGqzuIQ3RVgO9FwU9yI8E7+PLopcOSCmWoW4A4yR9J/iF8kNJtwTLF0v6z+DxpyQdlnxBLhBFZvbZpEsFdivxRdgVPK9W4hqzXDhMidOROyW9SomwkKzvseZJarDE4IJDgp6+U82se0DEI5J+amZHBfV/J40afhRsf6Kk70rqvk7wQUnfN7OTgsD0O0kPJ133NkPSFdo3uKA8eP6iu3cGy26X9E0zOze4kP/VZvZxMzssjbpkZmcHrz1AiZ7BFu37OUjSx8zsfWZ2oBLX1s0JTjs/I+ktZvZlMzsg+O9sM3ubu3cpcZ3fDcEggRGWGBBxkBKhq0tD/3wPk9QoaY+ZHa/U4amfoF0eU2KwxauCa+WSr3scsO6h9p1GW3WrlvSaFKf471HiGsqLRKgDUiqJUBf8wv83SePNbJESF0O/Plj9Q0kfNLOFSpyW2SKpM+WOgMJ72nrPU/d4sPxsSXPNrFHSU0pcf7U+WHeVpLuD02OfS7HPTNyjxOm3LZJWSJrTZ/0dSpwWqzOzJ4JQ8AklTntuUKLX6R9K9PJJiQvoNwXrnlV6X85PSlqgRO/ahOCYUiL43KtEaNugREhIDokzlAg33aHuRSWCafdzuft8Sf+lxGnB3Upc2H9ZGjV1O1yJYLg7eF87JV2ftP4BJQaH7FJi4MOXguM2KDFw4BIleu62S7pWUvcAjR9KWirp5eC110raz92blBh48FLQ5snX/iX7tRIDS/Yo0WaPZfCerlDilOl2JS5l+Wf3ijTqHsxQbdV9jFVKBPb1wXs8Llj+khIh8BV3T+e0PVByzH3QMwtFy8xGSXrG3U8NrttY7e6vH+I1h0pa5e4MlgAiwMxc0inuXhF2LZkys7skVbn7L8KuJS7M7HlJD7j7P8KuBYiikuipc/d6SRvM7LNSz3xJpwWPjwlGm0mJ0XN3hlQmAGAAwSn8d2nf6XcAfcQy1JnZg5JmSxptZlVm9nUlRmh93cwWS1qufQMiyiStNrM12jfHEwAgIszsbknTJH0vOAUMIIXYnn4FAAAoJbHsqQMAACg1hDoAAIAY2H/oTYrLkUce6SefnI85WIvb3r179epXvzrsMiKHdumPNkmNdkmNdkmNdumPNkltwYIFte4+Mhf7il2oO/bYYzV//vywy4ic8vJylZWVhV1G5NAu/dEmqdEuqdEuqdEu/dEmqZlZzuZd5PQrAABADBDqAAAAYoBQBwAAEAOEOgAAgBgg1AEAAMQAoQ4AACAGCHUAAAAxQKgDAACIAUIdAABADBDqAAAAYoBQBwAAEAOEOgAAgBgg1AEAAMQAoQ4AACAGCHUAAAAxQKgDAACIAUIdAABADBDqAAAAYoBQBwAAEAOEOgAAgBgILdSZ2Z1mVmNmy4bY7mwz6zCzzxSqNgAAgGITZk/dXZIuGGwDMxsh6VpJzxaiIAAAgGIVWqhz95mSdg2x2XckPSqpJv8VAQAAFK/IXlNnZsdL+pSkv4ddCwAAQNSZu4d3cLNRkp5x91NTrBsv6U/uPsfM7gq2+9cA+7lc0uWSNHLkyDMfeeSRvNVcrBobG3XooYeGXUbk0C790Sap0S6p0S6p0S790SapnXfeeQvc/axc7CvKoW6DJAueHiOpSdLl7v7EYPscPXq0r169OseVFr/y8nKVlZWFXUbk0C790Sap0S6p0S6p0S790SapmVnOQt3+udhJPrj7Sd2Pk3rqBg10AAAApSq0UGdmD0oqk3SMmVVJ+pWkAyTJ3W8Nqy4AAIBiFFqoc/dLM9j2sjyWAgAAUPQiO/oVAAAA6SPUAQAAxAChDgAAIAYIdQAAADFAqAMAAIgBQh0AAEAMEOoAAABigFA3TI2tHara3RR2GQAAoMQR6obps7fO1vuunR52GQAAoMQR6oZp5bb6sEsAAAAg1AEAAMQBoQ4AACAGCHUAAAAxQKgDAACIAUIdAABADBDqAAAAYoBQBwAAEAOEOgAAgBgg1AEAAMQAoQ4AACAGCHUAAAAxQKgDAACIAUIdAABADBDqAAAAYoBQBwAAEAOEOgAAgBgg1AEAAMQAoS6Hllbt0X/c/JJa2jvDLgUAAJQYQl0OXfX0ci2qrNOyLXvCLgUAAJQYQh0AAEAMEOoAAABigFAHAAAQA4Q6AACAGCDUAQAAxAChDgAAIAYIdQAAADFAqAMAAIgBQh0AAEAMEOoAAABigFAHAAAQA4Q6AACAGCDUAQAAxAChDgAAIAYIdTly+8z1YZcAAABKGKEuR66ZuDLsEgAAQAkj1AEAAMQAoQ4AACAGCHV54GEXAAAASg6hLocs7AIAAEDJCi3UmdmdZlZjZssGWP9FM1tiZkvNbJaZnVboGgEAAIpFmD11d0m6YJD1GyR90N3fIek3ksYVoigAAIBitH9YB3b3mWY2apD1s5KezpF0Qr5rAgAAKFbFck3d1yVNCrsIAACAqDL38MZqBj11z7j7qYNsc56kWyS9z913DrDN5ZIul6SRI0ee+cgjj+S+2AFcNnlvz+NTjtxPa+u69LNzD9ZbjhpRsBrS0djYqEMPPTTsMiKHdumPNkmNdkmNdkmNdumPNkntvPPOW+DuZ+ViX6Gdfk2Hmb1T0j8kXThQoJMkdx+n4Jq70aNHe1lZWWEKlKTJE3oeHnHEEVLdbp1xxhk6e9TRhashDeXl5SpouxQJ2qU/2iQ12iU12iU12qU/2iT/Inv61czeIOkxSV929zVh1wMAABBlofXUmdmDksokHWNmVZJ+JekASXL3WyVdKek1km4xM0nqyFX3JAAAQNyEOfr10iHWf0PSNwpUDgAAQFGL7OlXAAAApI9QBwAAEAOEOgAAgBgg1AEAAMQAoQ4AACAGCHUAAAAxQKjLofBuuAYAAEodoS4PLOwCAABAySHUAQAAxAChLg84DZsb7q7vPrRQs9bVhl0KAACRR6jLIU675lZHl+vJRVv1lTvmhV0KAACRR6gDAACIAUIdAABADBDqAAAAYoBQBwAAEAOEOgAAgBgg1OXBpKXb9albXpI7k5uUgjnrd2pva0fYZQAAShyhLg/ufGmDFm6uC7sMFEB1fYsuGTdHP3hkUdilAABKHKEOGIamtk5J0urtDSFXAgAodYQ6AACAGCDUAQAAxAChDgAAIAYIdQAAADFAqAMAAIgBQh0AAEAMEOqAYWCCaQBAVBDqgBwws7BLAACUOEIdAABADBDqAAAAYoBQl0OcgQMAAGEh1AEAAMQAoQ4AACAGCHUAAAAxQKgrUu6u8tU1zJMWMlofABAVhLoiNX5BlS7758t6ZH4lwS4CGCMDAAgboa5Iba1rliT95NGlOumnEwl2AACUOEIdAABADBDqEFl0PgIAkD5CHSKPSZ0BABgaoS4mdje1h10CAAAIEaEuJr70j7lhlwAAAEJEqIuJ1dUNYZeQN8nX1rV3dumPU1arsbUjvIKScN0fACAqCHUx0dkVv3SR6lq6xxdu0U3TK/THKasLX9BguO4PABAyQl0Ohd1rM2nptnALKID2zi5JUmtHV8iVAAAQLYS6PCp0yFuwaXdhD1ggHV2unY2tYZcBAECkEepQFM787bSwSwAAINIIdQAAADFAqAMAAIgBQh0AAEAMEOpyiNtZ5d8DczdLktoY/QoAQC+hhTozu9PMasxs2QDrzcz+amYVZrbEzN5V6BoRPcu31kuS6lu4LRoAAMnC7Km7S9IFg6y/UNIpwX+XS/p7AWoCAAAoSqGFOnefKWnXIJtcLOkeT5gj6Ugze31hqsuNr/5zXtglAACAEhHla+qOl1SZ9LwqWFY0XlhbG3YJAACgROwfdgG5YGaXK3GKViNHjlR5eXkodezZs6ffsnzVsnFjW79llVWVKi+vSbl9Y2NjaO2SrY4+97NNrr+2tjYn72e47bK1MTFgo6mpqejadyDF+FkpBNolNdolNdqlP9ok/6Ic6rZIOjHp+QnBsn7cfZykcZI0evRoLysry3txPSZP6Hl4xBFHSLt736orX7UsbF8jVazttezEE05UWdmYlNuXl5fnrZZ8ae/skp6d1PO8rKysp72POeYYlZWdNexjDLddKmoapRdn6FWvelXRte9AivGzUgi0S2q0S2q0S3+0Sf5F+fTrU5K+EoyCfbekPe4e/zvWAwAAZCG0njoze1BSmaRjzKxK0q8kHSBJ7n6rpImSPiapQlKTpP8XTqUAAADRF1qoc/dLh1jvkv6nQOUgBDdPr9C5Jx2ts0YdHXYpw+BDbwIAQAFE+Zo6xNz1U1ZLkjb+4eMhVzJ83EwEABC2KF9Th0GU+i3JnA4yAAB6IdTFCDkHAIDSRaiLsTXVDXpyUcpZYAAAQMxwTV2R+vO0tUNu85EbZ0qSLj69qG7EAQAAskBPHSKL6+YAAEgfoQ4AACAGCHXAMNCbCACICkIdkANW6nPMAABCR6jLIQt5CtpSihVkKAAAeiPUAQAAxAChDgAAIAYIdQAAADFAqMuheRt3hV0CAAAoUYS6GFlfuzfsEnKKwRAAAKSPUBcjz6+qCbuEksM0dQCAqCDUATlApyIAIGyEOgAAgBgg1BWRBZt2ybkvFQAASIFQVySeW1mtT/99tu6dsynsUiIhn9l21942rd7ekL8DAACQB4S6IlG5q0mStHJbfciVREX+Ut1H/zxTH/3zzLztHwCAfCDUFYkdja2SpAfnVYZcSfztaGgNuwQAADJGqCsSDS0dYZcAAAAijFAHAAAQA4Q6AACAGCDUAQAAxAChrkgwPR0AABgMoQ4AACAGCHUAAAAxQKgDAACIAUIdAABADBDqgGFgAAsAICoIdUWo+z6wcVdMgcks7AoAAKWOUFckkkPD+6+bHl4hAAAgkgh1BbCzsVWdXUXU7VQU6BoDACAZoS7P6lvadeZvp+maCSvDLiVWON0JAEBvhLo8q29ulyRNWb495EoAAECcEepQlIppEAUAAIVAqAMAAIgBQh0AAEAMEOqAYXBxHhgAEA2EugKpaWhRRU1j1q/nGrJoM6ZYAQCEjFBXIO2drvNvmBF2GQAAIKYIdQAAADFAqAMAAIgBQl2R4A4KAABgMIS6PLMcpTEyHQAAGAyhDkWK4cAAACQj1BUJIkw0MdUMACAqQg11ZnaBma02swozG5ti/RvMbLqZLTSzJWb2sTDqHI75G3eFXQIKgGseAQBhCy3UmdkISTdLulDSGEmXmtmYPpv9QtIj7n6GpEsk3VLYKofvuw8tCrsEFMCq7Q16Ye2OsMsAAJSwMHvqzpFU4e7r3b1N0kOSLu6zjUs6PHh8hKStBawPyMjfy9eFXQIAoITtP9wdmNmrJbW4e2eGLz1eUmXS8ypJ5/bZ5ipJz5rZdyS9WtL52dYJAAAQZxmHOjPbT4lToV+UdLakVkkHmVmtpAmSbnP3ihzVd6mku9z9T2b2Hkn3mtmp7t7Vp6bLJV0uSSNHjlR5eXmODp972da2ZUtr1vsvLy9XY2NjZNtloLo6unzA7Wprd+bk/QzWLunsv7Jh30dx9+7dkW3jTET5sxIm2iU12iU12qU/2iT/sumpmy5pmqSfSlrWHbDM7GhJ50m61swed/f7htjPFkknJj0/IViW7OuSLpAkd59tZgdLOkZSTfJG7j5O0jhJGj16tJeVlWXxtrI0eUJGm2db2/N7lkmbN6W1/7qmNh128AHS5Ik9y8rLy7M+dt4EbTdQXW0dXdKzk3qel5WV9bympu2AnLyflO0yRF3JVm6rl156QZJ05JFHqazs3cOuKWyR/KxEAO2SGu2SGu3SH22Sf9mEuvPdvV2SzOwwSQ2S5O67JD0q6VEzOyCN/bws6RQzO0mJMHeJpC/02WazpA9JusvM3ibpYEkleTV6ulNnNLd16vSrp+rL735jfgsKWXVDej2XAACUioxCnZn9m6SxZlanxCCLIyV9ou923aFvMO7eYWZXSJoiaYSkO919uZldLWm+uz8l6f8k3W5m31di0MRl7swMNpimtg5J0oSl20KuBAAAFFKmPXXflvQpd+80s2OUOA2bNXefKGlin2VXJj1eIem9wzlGXDAPWvQ5U0QDAEKU6ZQmDd2jXN29VtJ3c18SUolzphs/v1Lb97QMuV1dU1sBqgEAoDhlGurMzD5kZseY2beVmFsOGJYf/WuJvnLn3H7L+/Z8La7a0/M4ziEXAIBsZBrqvivpHZJ+rsQcc7/JeUUoSbWN0eiF45JNAECxyvSaup9JOkjSi5KWuPvTuS8JCM9zK2uG3mgA5EEAQJgyCnXu/iszG6nEpMOXBdORrHX3a/JSHXqQF3rLV3vs5ro9AECRynieOnffocSI1YmSFFxbhzyjF4hTowAADCbre7+a2SmSxkoaetgikAP3zB76jhqFRs4EAERFpgMlkt0r6V+S3i9JZnaqmd2Tk6qAFJZU1YVdwqDIdwCAMA0n1O3n7pMkdc9bt0zSqTmpCiWppb1T3394karrw+v8bW7vDO3YAAAMx3BC3dZgoIRLiQnsJB2Sk6pQkqYs367HF27RbyesDK2GW8vXhXZsAACGI+tr6iR9T9I/JL3OzP6fpAskLctJVcgapwCHp62zK/sX0/gAgBBlHercfaOZXSDpPySdJmmGpDtzVRiGhzsuAABQWjIOdWZmHswt4e4dSgyW+NdA2yA3jJQGAAAGkc01ddPN7Dtm9obkhWZ2oJn9u5ndLemruSkPSBZusm3v7NLEpduYLw8AEEnZhLoLlBjx+qCZbTWzFWa2QdJaSZdK+rO735XDGgFJUm1ja6jH/8u0tfr2/a/o+VX7biW2edfensfORXUAgBBlc0eJFkm3SLrFzA6QdIykZkmNwelYIO/C6C3buqdZkrS7qb1n2Tfve6XgdQAAkMpwpjSRu7e7+zZ3r5N0q5kdIklm9oGcVAcUXPaneI3hKQCAEA1nSpO+fiXpTjPrkLRI0swc7hsZivuJwK4IvkFOvwIAwjSsnro+fiNptRJ54pEc7hcZOOu308IuocgRzAAAxSmXPXU/dvdaM3u1pL9I+kYO913yGHAJAAAGk8ueut2S5O57Jf13DvcLAACAIeQy1I0zs1cFj9+bw/1G1prqhrBLQITQmwoACFMuT79eKemOUhoo8ZEbY/8WEZi2olqrtxPiAQDRlctQ1z1Q4k1ioASKVuppSb5xz/wC1wEAQGayDnVmZpJOUeJbcI0YKJFXpXTv15Xb6lXX1KZDDhxR8GOXUjsDAOIlq2vqzGyMpOWSHpP0qKSVkl4rMVAiKoo5m1TUNOo/b5kVdhkZ45I6AECYhgx1ZvZ2M7u/z+J/SvqZu5/q7qdK+omknm3cvTO3ZaIUWFI32fravYNsCQAA+krn9Os0Se/pfmJm35J0vKTDzewr3YslvcbMvu3ut+S+TAAAAAwmnVD3EUnXSPpi8PxMSQdIOrvPdgcG64CcCGOKkOEc05nTBAAQoiFDnbsv1b5AJ3f/hpmtk/QXd6+QJDN7k6RPuPvX81YpAAAABpTt6NcfSJphZo8Hz/9D0hW5KQnITm1jq6rrW/T2447Ieh97WztyWBEAAIWTVahz9yfNbKkSp2Yl6c/dvXZAWD58wwztbmrXxj98POt9NLczxgcAUJyynqfO3ddLujWHtQDDsrupPdTjc0UdACBMubz3K/KIa/ABAMBgCHWIrGeWbAu7hJQY5QoAiCJCHSLrh+MXh11CRsh6AIAwEeqKBPckjQ7jhwEAiCBCHQAAQAwQ6opEpn1DO/e25aWOKHl68dawSwAAIDIIdYiMTAcgfOfBhXmqBACA4kOoAwAAiAFCXZEohYGV7Z2l8C4BAMgPQh0iY21NQ9glZKTv6eJFlXUaNXaCNu3cG1JFAIBSRqgDsvTjfy1JuXzltuIKpwCAeCDUAUkymYJu/IKq/BUCAECG9g+7AEDKfMqWfDH1vn5xR0OrqutbwioHAIC0EeqKRNxvQRXVt/fRP8/UrhKY8w8AUPw4/Qok6RsuCXQAgGJBqAMAAIgBQh2QJO6nuQEA8RVqqDOzC8xstZlVmNnYAbb5nJmtMLPlZvZAoWsEAAAoBqENlDCzEZJulvRhSVWSXjazp9x9RdI2p0j6qaT3uvtuM3ttONWGL5OpNgAAQOkJs6fuHEkV7r7e3dskPSTp4j7b/Jekm919tyS5e02Ba4wEd9c9szeFXQYAAIiwMKc0OV5SZdLzKknn9tnmLZJkZi9JGiHpKnef3HdHZna5pMslaeTIkSovL89HvTmRTW1tw7wnanl5uRobGyPdLu1tbarcXDn0hn2kek+ZvM/B2mWg5atWrVR5Q8WA+1y+fJkOrl2Vdg1RE/XPSlhol9Rol9Rol/5ok/yL+jx1+0s6RVKZpBMkzTSzd7h7XfJG7j5O0jhJGj16tJeVlRWmuskTMn5JNrW1tHdKU/tl2YyOWV5entWx8yqp/Q448ECd+IYTpA3rM9pFr/cU7C+T99mvXZJqKisrS/kzHj36rSo768QBf/5r247SR085Wacef0TadURJJD8rEUC7pEa7pEa79Eeb5F+Yp1+3SDox6fkJwbJkVZKecvd2d98gaY0SIS902c5ftm5HY44rQdRMXr5dn/jbi2GXAQAoMWGGupclnWJmJ5nZgZIukfRUn22eUKKXTmZ2jBKnYzPrysmTrXXNWb3ui7f/+N50AAAgAElEQVTPzXElxWdRZZ3e+4fney3btbdNjS0dIVUEAEDxCy3UuXuHpCskTZG0UtIj7r7czK42s4uCzaZI2mlmKyRNl/Qjd98ZTsW50dbZFXYJobtx6hptSRGK75+7OYRqAACIh1CvqXP3iZIm9ll2ZdJjl/SD4D8gEoz5ZQAAEcQdJQAAAGKAUIeCK5aOrqrdTWGXAABA2gh1iAXPw01b33ft9JzvEwCAfCHUZYkbvwMAgCgh1KGoubsmLt3GqGIAQMkj1BVYkVxOlle5bIMZa3bo2/e/ohumrulZ1tLeqc07uR4OAFBaCHUoanVN7ZKkbXUtPcu+ed8CfeB6rocDAJQWQl2B7czy9mJIX/nqHXnd/08fW6K2Dk73AgCihVAHZKi90/Xsiu1hlwEAQC+EOiALxtWRAICIIdRlycWcJgAAIDoIdSg47p0KAEDuEeqALNBTCwCIGkJdlrijBAAAiBJCHQqOk68AAOQeoQ4AACAGCHUAAAAxQKjLEpfUxc8rm3envS3XVAIAooZQBwT+85ZZYZcAAEDWCHUAAAAxQKhDweVj7mHOhgIASh2hDrHgBb7IjRAJAIgaQh1iIdWtxwod9AAACBOhDiFg+mEAAHKNUJcleoFKG7EUABA1hDoAAIAYINQBWaCfFgAQNYQ6AACAGCDUAQAAxAChLkucfgMAAFFCqEPB5eOOEoXG6GcAQNQQ6rLEdzoAAIgSQh0AAEAMEOpQMBU1jTrt18+qur4l7FIAAIgdQh0K5r45m7SnuV1LqvbkfN9c4wYAKHWEOhS1OAy6AAAgFwh1KGp00AEAkECoQywYXXYAgBJHqMsaXURRkuqaOnrxAAClhFCHgqEzDQCA/CHUlZimtg61d3aFXQYAAMgxQl0JWLZl3xQiY66coq/eOS/EagAAQD4Q6krAJ/72oqoa9vXOzVq3U60dnSFWBAAAco1QVwRyccF/fVvvnSzNwwTAxYxT0gCAYkeoyxIjK+Plt8+syGh7fv4AgKgh1GWpkN/pxThq9OGXN+vLd8zttcwU3Tcyd8OusEsAAGBY9g+7AAxubXWDfvHEspzvN9+h9CePLk1xzNwftaG1I+f7BACgGNFTF3G/fHIZvUiD+GUeAi8AAMWIUAdk4XsPLwq7BAAAegk11JnZBWa22swqzGzsINt92szczM4qZH1xEt2r2XKDcQsAgFIXWqgzsxGSbpZ0oaQxki41szEptjtM0nclze27DpmpaWjpeRzG6M0oD5QAAKDYhdlTd46kCndf7+5tkh6SdHGK7X4j6VpJLSnWhaYYp7Q455rnwi4BAADkSZih7nhJlUnPq4JlPczsXZJOdPcJhSysFBTjNCkAAGBgkZ3SxMz2k3SDpMvS2PZySZdL0siRI1VeXp7X2iTprmWtWb82k/rq6pqzPk6y5uZmJV9Zt3DhQu3dOCIn+x5M8nutqsq+zYayo6am/7FnlGu/IdJrY2OjysvL1djYlPOaCvE5zIfuNkFvtEtqtEtqtEt/tEn+hRnqtkg6Men5CcGybodJOlVSuSW+mF8n6Skzu8jd5yfvyN3HSRonSaNHj/aysrI8lp3wnelTsn5tJvXdtmaOtGtn1sfqdsghhyj5DPbpp5+hc046etj7HdDkROdq8nud2bBC2rQhL4d77WtfK23f1mvZBz9YphH7DR7qysvLVVZWpkMXzZQaG3JaUyE+h/nQ3SbojXZJjXZJjXbpjzbJvzBPv74s6RQzO8nMDpR0iaSnule6+x53P8bdR7n7KElzJPULdEhPFE63RqEGAADiKrRQ5+4dkq6QNEXSSkmPuPtyM7vazC4Kq650FeE4iVgb7s9jS45OcwMAEJZQr6lz94mSJvZZduUA25YVoiaUpoYWbjcGAChu3FEiS16gOU1ydb/Uvmc+C1V/72MW/JAAAJQMQl2WyCcAACBKCHVZKvZeJwth1AIDJQAAyB9CXZZydVp0KHPW7yrIcQAAQHEj1GWp2Hrqbl/ae+LfMK6pi5LKXU369z+Wq66lK+xSAADICUJdiahtLu0Q19fdszZqfe1ezd7WGXYpAADkBKEuS0SkiOEHAgAocYS6bBV5iAijfMZJAACQP4Q6xMKEpduG3iiFzq4iT+cAAAQIdVkq1OhX5F5FTYOeX1UjSXp0bXvI1QAAkBuh3iasmBX74NFSPhV6/g0zex4X+Y8RAIAe9NRlKZ9hYEdDqxpa8tuDRJgBACBeCHURdPY10/SB66aHXUbOFfqOEu6uhpZ2VdQ0FvbAAACEgFAXUbub4netVxinrC+9fY7Ov2FG4Q8MAECBEeqyVOp3ZCgWy7bUh3bsUWMnaEtdc2jHBwCUFkJdloo90pFJC2Pu+p1hlwAAKBGEOhTE1rpm1TVH85Ty/zzwStglAAAwbExpkiV6ujLzb394PuwSJEnlq2v6LZuwJLuJiwEAiBJ66kpUoUeiRsVl/3w57BIAAMgLQl2JoqcRAIB4IdQBAADEAKEOAAAgBgh1EVbX1BZ2Cb1U7mrS9x5aqLaOrrBLAQAAfRDqIuz0q6fmbd+exUx7P31sqZ5YtFVzNww891pLe6eWb93T8/zFtbUlPVFzV+m+dQBAgRHqkFNjH12ij//1xZ7nX7pjru6buznEisJ154sbwi4BAFAiCHUl6m/PVeRlv69sruu3rHJXU16OVQzW1zaGXQIAoEQQ6krU7ALevqq2obVgxwIAoFQR6kIwc82OsEsoqNXVDWGXAABA7BHqQjB+QVXYJQxLsYx7eGFtbdglSJKWVNXp/dc9r/qWaN77Fsin9s4ubajdG3YZQEkg1CFt2d5aLKxbkl0zcWU4B+7jxqlrVLmrWfM37gq7FKDgrnpquc77Y7l2cBkGkHf7h10Aom/ltnot2LS753ltY2a/nE0leqNZlfZ7BzbvbNL9wej3+pZ2jTzsoJArAuKNnroQFNu8bRf+5QX94ollPaddf/DI4gG3DatXDkD0fPjGGWGXAJQUQh1yqsjyKoA8auXuM0BBEeqQd6Xce5d85w4CLwAgnwh1QJ5ZKadaAEDBEOpCUKwdNsWWTaJQLgMlAACFQqhDTqUKfmHFmr7h+flV1aHUASAaf2QBcUeoC0Gp/XJbXLUnlONW1Oy77+oLa3foa3fND6UOAMV7hgIoJoQ6lIQv3zEvlOM2t3dq6ZZEqGWgBAAgnwh1IYjTd/u0FdX60J/K1dHJ1AUDYSZ9AEAhcEeJiHhmyVa94ehXqbOruCLf2MeWqLaxTbubmC0eAIAwEeoioKOzS1c8sDC048+qqNW5b3qNRuyX+dV+tY1teagIAABkitOvERBm39xLFbX6wj/m6ubpFSFWURqKbUoYAEBxoacuDBE6w1pd3yJJ2lC7N6PX1TS0aHFlOKNaixUDJVDK+JsGyD966iIgKl/201fXaNzMdfI0Cmpo6dDEpdsKUBWAOIjIrzkg1uipC0ME/2Rt7ejU//vny5Kk0a87XB98y8iQKwIAAJmgpw6SpI7OfX9HN7d1DLl9Z5fr8YVb8lkSAADIAKEuDN73afGdmHh0QVXYJQAAgCSEuhBMWhada9EaWhK9cr1HZg59fri1I/VkwxE8swwAQEkg1IUgSvML/+qp5ZIkGyCObd7ZlNH+IvTWIoe2QSnjDz4g/0INdWZ2gZmtNrMKMxubYv0PzGyFmS0xs+fM7I1h1FkKkk8BJ/fafeD66SFUAwAAMhVaqDOzEZJulnShpDGSLjWzMX02WyjpLHd/p6R/SbqusFUWRhSmNEmugb+oAQAoPmH21J0jqcLd17t7m6SHJF2cvIG7T3f37vN/cySdUOAaS5Jx6wMAORaBv12B2AtznrrjJVUmPa+SdO4g239d0qRUK8zsckmXS9LIkSNVXl6eoxLzJ7nG9pAusvvNfVN7Hs+pqOl5vGzZUh1Qs7Lf9rt37e55XFnVe/TrrFmzdMRBpubm5jxUGg8DtWvUNDY2FsW/oUKjXVJLt13mzZ2nykNL5zJuPi/90Sb5VxSTD5vZlySdJemDqda7+zhJ4yRp9OjRXlZWlv+iJk8Y1suTa2xp75SenTzMgjJ3x7K2nsf1bfuC5TtOfYfKxhybeJL0Po86+ihpZ60kafb23kH0nWeeo9cdfrBe9fJMqSmzwRWl4tTkdo2w8vJyFeTfUJGhXVIbtF2Sfn+cc+45evPIQwtTVATweemPNsm/MEPdFkknJj0/IVjWi5mdL+nnkj7o7q0Fqg1DaGrr7PX8vD+Wa+RhB+nVB44IqaLo21DbKCn6oQ4AUJzC7At/WdIpZnaSmR0o6RJJTyVvYGZnSLpN0kXuXpNiH6HoitKcJHmQ7SV1OxrI3IP53cRVYZcAhIYrdYH8Cy3UuXuHpCskTZG0UtIj7r7czK42s4uCza6XdKik8Wa2yMyeGmB3BfXgy5vDLiGvBgp1S6r2DPnaeMddAACiK9Rr6tx9oqSJfZZdmfT4/IIXlYbq+uH3SH3m77M0/pvvieRI04EmIt7T3D7ka7uiMD8LAAAlqHSGIkXM/E271d4ZvwBUuYvRrwAAhIFQFwGR69yKXuchgCIXtV9zQBwR6kIUwTOvAACgSBHqQhTVTBfVugAUr7aOLrm7OmM+ewAQpqKYfDiuojhIQtpXV2NrR8iVAIiLJxZtUeWuJk1cul0b//DxsMsBYomeuhBFM9Ltq2v7npZQ6wBQXDo6uwZdP3Hp9gJVApQmQl0WohrGco3JhAGk68lFW3Tyzydp/Y7GsEsBShahLgu5uiJkW300e8LMpOVb9+jS2+eEXQqAIjF5WaIXbtX2hpTrB5r/cri272nR4wur8rJvoNgQ6kL03j88L0nyCA72r6jhr+182Mt1ioip7qmZHpi7WfM37hpy+/qWdrV2dA653VC++I85+v7Di7kGGBChDimYLHpz58XEN+9bEHYJQF5MXp7oqXuxolafuXX2kNu/86pn9fnbhn82oCa4ww93swEIdVkplWvqkHvzN+4OuwQgMhZV1mnqimo9v6p62Psi0wFMaRIJd8/aFHYJvZhF85RwHDS3D/90ExAl7q4tdUPfHnCgGZz+6575kpT9NCf8lQ30INRlYNqKarV1dqklB9eBJLt28qqc7m+4TPzVCyChpqFFm3Y26exRR6dcP2lDu/7flOcz2uee5vZclNYbv7MAQl0mvhH8RRl7RqgDkPDJv72o6vrWXj1pdU1tOuKQA2RmWrlr8LnpUjnt18+mtd2MNTv01Tvn6ZVfflhHv/rAlNvQUQfswzV1AIABVdf3nq9y/Y5GnX71VH3m1tm6b87Al47MXb+z1/Nsbg92+8z1kqR3/WZqxq8FShE9dejHZJzJAJDSxp17JUkLNu3Wgk279c5jRqTc7vPjeo9sHRcENAD5Q08dUmJ6AKD47GxsVeWuprDLyJkXK2rT3pbBXQChLnRRnAl9oFFqyK26pjYtrdqjMVdOjtUXMcLz7t8/p/dfNz1v++/o7NKuvXkY5JBCfUt6xzF+YQE9CHUh+/7Di8MuoR+TGElWAKdfPVWfvOlFNbV1asLSbWGXgxD93yOLdc/sjSnXzVizI+07L7R3DvwPt6KmUX+ZtjaL6vb55ZPL9cPxvX9nLanNzzQ9bR2ZDcDg5AJAqEMKy7fWcyojj878zVR958GFYZeBCHn0lSpd+eTyfsuXVNXpq3fO0zUTVg76+nkbdunDN8wYdJtLxs3RjdPWaE9T9j1tk5YV7o+PrjQHVtBRB+xDqEM/v5+0UvfN2Rx2GbG1c2+bnl68tdcyvpeQyu4ggG2o3Tvodlc/s1xrh7hfc1uK3r7G1g5dfPNLqqhpyL7IPOmk6w3IGKEO/bR3upZu2RN2GSWF3gak4nkINsm98C+s2aHFlXX645Q1OT/OcDW0dPR6/uiCKnV1ubbWNeu6yavS7skDSglTmgARYPTVIQSZxKLJBTz1KkkfuXFmr+f/N36xWju69MTCLZq3cZc++vbX6bQTj+xZT8QD6KkDIoGeOiRr7ejUzdMr1DHIwIdkTa3ZDVbo7ghM5/P3zfteUd0wrsfLhd1NbWrrTAygqK5vkcSlC0AyQh0QAUzLgGT/eGGDrp+yWnfP3ihp8M/HN+9doPVDXHM31D6i9PGr2t2kxtaOlOuST0dffu+CAdcBpYrTr0AEROg7FRGwNwg1re1DT+sxefn2jPadnH2iOMr9fddO11tfd1jKdV3e+3Zj3b11krRtT4tec+hBA+63pqFFSyr36Pwxx+auWCBi6KkDgIgpVNTqOf06wJ8Vt4d0a69V21OPxr1h6ppeg7jO/d1zPSOEv3H3/EH3ecm4OfrGPfOzugctUCwIdUAEROn0F0rQAJ+/ayYOPj9elDS1pT5l221jcIqa07SIM0IdEAE7GlrDLgER0pM7rNf/hmVP88CDHDbv7H+bulkZ3Hc1Fzo6M7uDRF/pRrXqhtZhTcAMRBmhDoiAW8rXhV0CQrClrlmrttf3W959rZv1PM+P7v2mmpfy/nmFnYD8v+9dMGjwHErfee0G8t4/PK/3/OG5rI8DRBkDJQAgJO/9w/OpV2Qw1UimFlbu1r+/NXqDBZ5bVaPTfv1s3vZvZj1doE1t+6aAqW1sVU19q8Ycd3jejg0UCj11QAlo6+jSrHWFPZ2G7HX3oHX1DGTobU9Tu66fsiqrU5bJp1oHu76sM8058orVuJmJ3vEP3zBDH/vrCyFXA+QGoS5Nu/e2hV0CkLVrJqzQF26fq+Vbuf1bMZm3YVe/Zf/+x3KddvWzunn6Ok1ZXj2s/c9csy/ov7B2R691mU6VEnV9R73eOiMxsnc319chRgh1aZq6Yni/PIEwra5OTBExnGuWkFt1Tan/ULzrpQ3614KqAV+XPNHwxp1DTzo8kGeXb9ejr+w7TuWu5qz3BSAauKYuTR3MbYQiNmd9/x4fhOPP09boz9PWDrj+qqdXpL2vtdWp53MbyoJNu/vdkSFO0+q8sHaHTj/xSB128AEDbsPUJogjeurSFMWZ1xEvn711lmobczu1SVeX68M3zOh5/ruJK/kyy9LKnZ166y8n6azfTst6Hy3tnYMGuqG0dfS+hu6JRVsH3X5LXbOm9TnL4JI+/fdZ/bbtznT1Le266fnsawzb9j0t+vId8/Tpv8/qud1YqgmHXdLc9TsLXB2QX4S6NG3e1X8eJyCXXt64W2f9dpq6ctgr3NbZpbU1jT3Pl22pT+s+oaVi0869GjV2gpZU1fVbd9uMdRo/v7Ln+bUvt6ilvWtYwfvnjy/L+rWStCzDayIv+tuL+sY9ve+00DcYdlu+tV4dnV36zdMr9Mdn12RdY9i6JyFeU92oS8fN0Z7mdn3sL6kHQvzlueINr0AqhLo03TYjnNvloPS0dHQOvVEBrNvRmDLsxMlzK2skSY+9sqXfut9PWqUf/WtJTo+3qHL3oOsHGs1qlrh36T2zNmZ0vJ3BAK/k3tm/PV+Rctt752zSn6et7TXdR7G56KYXe96zlJh/79nl23uuKU1W19Su+hauMUW8EOqAiHGXVm9vUGtHp55ctGXYM+2n2n86PvSnGbroppdyeuyo6Qoao+/1ZMkhaNTYCfrWfQsGXD+QDbV7NW7mOk1Zvl3vv+55tXd2aW/r4IFp7GNLUy7fvbdN51zz3JCnWwfy1OJ9r+s+JZnKTdMrcnP7ipAsqdrTq3dV0qDBfNmW/hM/p1JT36KnF2fX9kAhMVACiJgJS7bpx4/u+yLaUtesb5edHGJF8dUd6kb0SXWTlm0f9HlHl+uAEYOnn8/dNrvX7d9eWLtD2+tbBn3NQKNeF1dlPhVNcvDclOI2YAMp4kwnSXpk/sAjh4eyanu9vvfQIo3/5nt6DbL44j/mam1Noz70ttfqVQfytYno4tMJRExyoJOyuy9sV5dr4rJtev0Rh/Rbl+qaqtaOTtU2tun4I/tv39eCTbt01KsO1JtGHppxXbm0uLJOrz38oJTvMV3dly/O2dD7gvmdQ8xL2R18Glra9Y6rEndB+MK5b9DvPvWOnm36/ty+dlfva9vyLbkzsT2D3t794jQMNgO/eGKpqutbtWp7g2at26mPvv11Peu21CWme2ESBEQdoQ6IOMui7+T+eZv1yydSX5T/xKIt/W6J9KPxS/TU4q265OwTVbmlVSeM2XcN0sk/m6jHv/1eveOEIyRJn/77bEnSvJ9/SK897OCMaxuO5rZOtXV26YhDDtDFN7+kA0fspzXXXJj1/rp76rpPwy3bske3lFfoba8b/JZRt81cr+unrNbpJx7Zs+yBuZs1edl23f6VMzPqGcuHpraOXqcW2zIIdc8sKc3TjPfN2ayPjEncPm2gs+udna411Q16y7GH5ey4bR1d6uxyHXLgiJztE6WLa+qAiFtb06D/vnf+gKMWU6kZ5DTf/I2956x7+OXNPddcPfRypV7a2tEz276UONX4yZtelCTtTboe65nF29KuJxfaOrr0weun97o/aFtn14DhdSh7mtt13eTVPc8/f9tsfe622Zq4dHu/nru+rp+SeN2iyt4DSXbtbdOn/z5bP3hkcVY15cpnb028l26ZDPQq5d6oZ3umf0ndCH+YvFIfuXFmr1utDccVD7yit/xikt525eSc7A8g1AER98LaWk1ZXq1bylOPWkxl5tqB7/P6yuZ9QWT8/Er95NH+F+cP1Df4rftfSbuG4Xhi4RYt7XMd2RlXP6uaFKei752zqdfzH41frF8/vTzlflvaO3sm7H38ld7XXs3dsKtn5GexT+W3fGt6AwCQWn1zR68BSt2fh/kbE6OXa/fmZj7JZ5YU9g8jxB+hDigS6Uxa6+4aNXaCFlcOPhVJ91QOA40MHCjTLNyUekoOd9f01TVZTWx8zYQVGjV2gh6at7ln2fceXqRP3vSitu1p1qade3Xlk8u0N2mqjVFjJ/Tax/VTVmlPc7u6ulzjF1Tpny9t1IV/eUFfv+tlLdycqHnX3jb974ML9eEbZ2r7nhY9vrD/NCbdZq1jUtpS9uNHl+jkn0+SJP3jhfVqbk989pLnfBw/v1Kjxk4YdDQxUGhcUwcUub9MW6vqhhb97lPv0O0vpHea7WePLdVnzjxhwPV7U3xR/Wj8YjUkLX9uVbW+9r6TJEkPv1ypsY8t1fWfeac+e9aJGdV/+wsbJCWm8xi/oEo/+PBbeta95/fP67QTjhhy9OfN09fp5unr9D/nvbln2cpt9Vq5rV7PrarRA984V1/4x9yedZffO19LshhRitLyt+fW6k9T+0/E7L7vD6I7X9ygScu2a0dDi1748b8PeW1cV5frTT+bmJd6AXrq0jDQjbcRT1df/PawSxjQnuZ2/fKJZb0mBb5x2ho9MHezOrtcv5u4Kq397G5q02X/fHnA9X2n8JCk8X2m23ipYl9v1tZgdGCq036fv222Ji5NfZqp7+2bFmzarS8mhS8ps+k8bp6+LuXyL/TZJ4EO6UgV6Pp6ctEWrdxWr9rGNs1cu2PQbTs6uwYMdN13wgCGg1CXhvbOIr/ABhn5yntGhV3CgE779bO6d84mXXTTS6ppaOl1mvX0q58d5JW9deVoPuNRYyfoc7fN1l+DuxTcNWujPnLjDL2YdE3f3A279O0BrsWbsrx/eASiLvkyg3U79t1277/vTUxS3dLeqcsm79Wvnlwmd1ddU5teXFurlwY5rT/myinavbdN77/uec0OtntuZbW++9DClPeuBVLh9GsaVm7jouNSceCI4vk755xrnuv1vKEl/b/0Z+fwRubzNvQeTbumulFfumOunrrivQPekmpp1R7taGzRBu5DiyL0mVtnD7iuua1Te4Net7tnb9LdszcNuG1fZ/xmqiTpOw8u1PxfnK+v352Y23D06w5jAnKkJdRvMDO7wMxWm1mFmY1Nsf4gM3s4WD/XzEYNtc+mDtemnf2/KL7/8CL9buLKrOrse9oJ8fX5szO7HgwDu+imlzS1Z4oIaUlVnW4pr1BnMEXK1+6ar/vnpP+FBxSDt105edgDbWobW3sNBrpu8mo1cJ9apCG0UGdmIyTdLOlCSWMkXWpmY/ps9nVJu939ZEk3Srp2qP3WNLnK/lguSZq2olqjfzFJU1dU6/GFWzRu5vpec1yla8Igk3Ge/7bXZrw/oBRddNNLum7yar056ZqirXsGv20WUIz+98GFOd/nO656Vi3do3CrG7R6e4M6OruK+lq8zTub1NjaoZqGlp4R+anUNLSoanfmcwPuaW4f8t7Zq7c36NJxc/SpW15SQ0u7pq6oVkt7p9xdVz+9QrPW1eqkn07QqLETerV1TX2LRo1NLG/tGPyezoVk2UxBkJMDm71H0lXu/tHg+U8lyd1/n7TNlGCb2Wa2v6Ttkkb6IEUf9PpT/PVf/bMOP3h/1WdwOipbj37r3/Tpv8/K+3Ek6bgjDtavLz5VN0xdwynhPPlW2Zv1kwve2m/KjHQcuP9+GU0QDACF8p43vabXZRfHHHqQahsT8+297fWHq6Ozq9eULZk47OD95S7tZ9Jxwa0GV21PzAd52olHavSxh+rk1x6a9kCu4Xjvya/pNYirGGy69hML3P2sXOwrzFD3GUkXuPs3gudflnSuu1+RtM2yYJuq4Pm6YJsBZ1btDnVx1H0DcffELP+lxiw/k8IeMMJ6DYbp+xwAgHzJZaiLxUAJM7tc0uXB09ZN134iu/sGxdsxkga+zUDpol36o01So11So11So136o01SG52rHYUZ6rZISr4q/YRgWaptqoLTr0dI6tev6u7jJI2TJDObn6vEGye0S2q0S3+0SWq0S2q0S2q0S3+0SWpmNj9X+wpz9OvLkk4xs5PM7EBJl0h6qs82T0n6avD4M5KeH+x6OgAAgFIVWk+du3eY2RWSpkgaIelOd19uZldLmu/uT0m6Q9K9ZlYhaZcSwQ8AAAB9hHpNnbtPlDSxz7Irkx63SPpshrsdl4PS4oh2SY126Y82SY12SY12SY126Y82SS1n7RLa6A6gHGwAAAZMSURBVFcAAADkTvHcEwkAAAADilWoG+q2Y3FlZiea2XQzW2Fmy83su8Hyq8xsi5ktCv77WNJrfhq002oz+2h41eeXmW00s6XB+58fLDvazKaa2drg/0cFy83M/hq0yxIze1e41eeHmY1O+kwsMrN6M/teKX5ezOxOM6sJ5sTsXpbx58PMvhpsv9bMvprqWMVigDa53sxWBe/7cTM7Mlg+ysyakz4ztya95szg315F0G4WxvvJlQHaJeN/M3H7nhqgXR5OapONZrYoWF4Sn5dBvpPz/7vF3WPxnxKDLdZJepOkAyUtljQm7LoK9N5fL+ldwePDJK1R4tZrV0n6YYrtxwTtc5Ckk4J2GxH2+8hT22yUdEyfZddJGhs8Hivp2uDxxyRNkmSS3i1pbtj1F6B9Rihxp5Y3luLnRdIHJL1L0rJsPx+Sjpa0Pvj/UcHjo8J+bzluk49I2j94fG1Sm4xK3q7PfuYF7WRBu10Y9nvLQ7tk9G8mjt9Tqdqlz/o/SbqylD4vg3wn5/13S5x66s6RVOHu6929TdJDki4OuaaCcPdt7v5K8LhB0kpJxw/ykoslPeTure6+QVKFEu1XKi6WdHfw+G5J/5G0/B5PmCPpSDN7fRgFFtCHJK1z902DbBPbz4u7z1RiZH2yTD8fH5U01d13uftuSVMlXZD/6vMjVZu4+7Pu3n3fxTlKzCs6oKBdDnf3OZ74drpH+9qxKA3wWRnIQP9mYvc9NVi7BL1tn5P04GD7iNvnZZDv5Lz/bolTqDteUmXS8yoNHmxiycxGSTpD0txg0RVBd+6d3V29Kq22cknPmtkCS9x5RJKOdfdtwePtko4NHpdSu3S7RL1/4Zb650XK/PNRau3zNSV6FbqdZGYLzWyGmb0/WHa8Eu3QLc5tksm/mVL7rLxfUrW7r01aVlKflz7fyXn/3RKnUFfyzOxQSY9K+p6710v6u6Q3Szpd0jYlusFLzfvc/V2SLpT0P2b2geSVwV+FJTkE3BKTfl8kaXywiM9LH6X8+UjFzH4uqUPS/cGibZLe4O5nSPqBpAfM7PCw6gsB/2YGd6l6/9FYUp+XFN/JPfL1uyVOoS6d247FlpkdoMSH5353f0yS3L3a3TvdvUvS7dp3yqxk2srdtwT/r5H0uBJtUN19WjX4f02wecm0S+BCSa+4e7XE5yVJpp+PkmgfM7tM0ickfTH4QlJwenFn8HiBEteLvUWJ9598ijaWbZLFv5mS+KxIkiVu7fmfkh7uXlZKn5dU38kqwO+WOIW6dG47FkvBdQt3SFrp7jckLU++HuxTkrpHJz0l6RIzO8jMTpJ0ihIXqcaKmb3azA7rfqzExd7L1Pv2c1+V9GTw+ClJXwlGIr1b0p6krvI46vVXdKl/XpJk+vmYIukjZnZUcPrtI8Gy2DCzCyT9WNJF7t6UtHykmf3/9u5eNYooDAPw+6GdiKCVlRALCxWbFIKWVoKVvSA23kPatFYqWFgI3oG1glegiVHwJ4qdlY2NRYqxmIkGye5qYDPJ2eeBgWWYXc58nJl5OTNn58jweSl93/gy1OVHVV0ezk+38qeOzdjDMbNI16lrSd53Xff7tuqi9JdJ1+Tsx7llnjNA9ntJP4PkY/r0vzJ2e/Zxv6+mH8Z9k2RtWK4neZpkY1j/LMnpHd9ZGer0IYd4ltGMuiyln122nuTddp9IcirJiySfkjxPcnJYX0keDnXZSLI89j7MsTbHknxPcmLHuoXrL+lD7bckW+mfV7mzl/6R/jmzzWG5PfZ+zaEmm+mf7dk+vzwatr05HFtrSV4lubHjd5bTh5zPSR5k+LP7w7pMqMt/HzOtXad2q8uw/kmSu39tuxD9JZOvyXM/t3ijBABAA1q6/QoAsLCEOgCABgh1AAANEOoAABog1AEANECoAwBogFAHANAAoQ7gH1TVvapar6r7Y7cFYDdHx24AwEFXVWeTXOm67tLYbQGYxEgdwBRVdS7JyyRnqur18B5hgAPHa8IAZqiq1SRfu657PHZbACYxUgcw28Uk62M3AmAaoQ5gtvNJ3o7dCIBphDqAKarqeJKtrut+jt0WgGmEOoDpLsQoHXAImCgBANAAI3UAAA0Q6gAAGiDUAQA0QKgDAGiAUAcA0AChDgCgAUIdAEADhDoAgAb8AhlpRGljN0IQAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 720x576 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.io import wavfile\n",
    "\n",
    "# read audio file \n",
    "fs, x = wavfile.read('../data/speech_8k.wav')\n",
    "x = np.asarray(x, dtype=float)\n",
    "N = len(x)\n",
    "\n",
    "# compute ACF\n",
    "acf = 1/len(x) * np.correlate(x, x, mode='full')\n",
    "# compute PSD\n",
    "psd = np.fft.fft(acf)\n",
    "psd = psd * np.exp(1j*np.arange(2*N-1)*2*np.pi*(N-1)/(2*N-1))\n",
    "f = np.fft.fftfreq(2*N-1, d=1/fs)\n",
    "\n",
    "# plot PSD\n",
    "plt.figure(figsize = (10, 8))\n",
    "plt.plot(f, np.real(psd))\n",
    "plt.title('Estimated power spectral density')\n",
    "plt.ylabel(r'$\\hat{\\Phi}_{xx}(e^{j \\Omega})$')\n",
    "plt.xlabel(r'$f$')\n",
    "plt.axis([0, 2000, 0, 1.1*max(np.abs(psd))]);\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise**\n",
    "\n",
    "* What does the PSD tell you about the spectral contents of a speech signal?\n",
    "\n",
    "Solution: It can be concluded from the shown PSD that the main power of a speech signal is contained in the frequency range below 500 Hz. The speech signal exhibits furthermore a harmonic structure with a dominant fundamental frequency and a number of harmonics."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Cross-Power Spectral Density\n",
    "\n",
    "The cross-power spectral density is defined as the Fourier transformation of the [cross-correlation function](correlation_functions.ipynb#Cross-Correlation-Function) (CCF)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Definition\n",
    "\n",
    "For two continuous-amplitude real-valued wide-sense stationary (WSS) random signals $x[k]$ and $y[k]$ the cross-power spectral density is given as\n",
    "\n",
    "\\begin{equation}\n",
    "\\Phi_{xy}(\\mathrm{e}^{\\,\\mathrm{j} \\, \\Omega}) = \\mathcal{F}_* \\{ \\varphi_{xy}[\\kappa] \\}\n",
    "\\end{equation}\n",
    "\n",
    "where $\\varphi_{xy}[\\kappa]$ denotes the CCF of $x[k]$ and $y[k]$. Note, the DTFT is performed with respect to $\\kappa$. The CCF of two random signals of finite lengths $N$ and $M$ can be expressed by way of a linear convolution\n",
    "\n",
    "\\begin{equation}\n",
    "\\varphi_{xy}[\\kappa] = \\frac{1}{N} \\cdot x_N[k] * y_M[-k]\n",
    "\\end{equation}\n",
    "\n",
    "Taking the DTFT of the left- and right-hand side results in\n",
    "\n",
    "\\begin{equation}\n",
    "\\Phi_{xy}(\\mathrm{e}^{\\,\\mathrm{j}\\,\\Omega}) = \\frac{1}{N} \\, X_N(\\mathrm{e}^{\\,\\mathrm{j}\\,\\Omega})\\, Y_M(\\mathrm{e}^{-\\,\\mathrm{j}\\,\\Omega})\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Properties\n",
    "\n",
    "1. The symmetries of $\\Phi_{xy}(\\mathrm{e}^{\\,\\mathrm{j}\\, \\Omega})$ can be derived from the symmetries of the CCF and the DTFT as\n",
    "\n",
    "    $$ \\underbrace {\\Phi_{xy}(\\mathrm{e}^{\\,\\mathrm{j}\\, \\Omega}) = \\Phi_{xy}^*(\\mathrm{e}^{-\\,\\mathrm{j}\\, \\Omega})}_{\\varphi_{xy}[\\kappa] \\in \\mathbb{R}} = \n",
    "\\underbrace {\\Phi_{yx}(\\mathrm{e}^{\\,- \\mathrm{j}\\, \\Omega}) = \\Phi_{yx}^*(\\mathrm{e}^{\\,\\mathrm{j}\\, \\Omega})}_{\\varphi_{yx}[-\\kappa] \\in \\mathbb{R}}$$\n",
    "\n",
    "    from which can be concluded that $|\\Phi_{xy}(\\mathrm{e}^{\\,\\mathrm{j}\\, \\Omega})| = |\\Phi_{yx}(\\mathrm{e}^{\\,\\mathrm{j}\\, \\Omega})|$\n",
    "\n",
    "2. The cross PSD of two uncorrelated random signals is given as\n",
    "\n",
    "    $$ \\Phi_{xy}(\\mathrm{e}^{\\,\\mathrm{j} \\, \\Omega}) = \\mu_x^2 \\mu_y^2 \\cdot {\\bot \\!\\! \\bot \\!\\! \\bot}\\left( \\frac{\\Omega}{2 \\pi} \\right) $$\n",
    "    \n",
    "    which can be deduced from the CCF of an uncorrelated signal."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example - Cross Power Spectral Density\n",
    "\n",
    "The following example estimates and plots the cross PSD $\\Phi_{xy}(\\mathrm{e}^{\\,\\mathrm{j}\\, \\Omega})$ of two random signals $x_N[k]$ and $y_M[k]$ of finite lengths $N$ and $M = 2 N$. The estimated cross PSD is calculated from the valid part of the CCF $\\varphi_{xy}[\\kappa]$ by an DFT in order to exclude boundary effects."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAnYAAAEZCAYAAAAaHmjrAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvFvnyVgAAIABJREFUeJzt3XuUHWWV/vHvk4RAIECCYA+ESEACCuiAtMCgaEcQkFEJM46Cyk00OqKjS5SbODAgAw4qynIEozCAKBFBIDIghkuL/oYAiUTCLSTcJE0kQNJAIELo7N8f9Z5Y6Zy+Jae7Tqqez1q1uuqtt6r2rjqne3ddzlFEYGZmZmbrv2FFB2BmZmZmjeHCzszMzKwkXNiZmZmZlYQLOzMzM7OScGFnZmZmVhIu7MzMzMxKwoWd2XpO0n6S5hUdRz2S2iQtLDoOK46kMyRdMYD+IWnHQYynad8vZo3gws6sIJKekLRc0rLc8IN+LLfaH76I+H1E7DxIMV4q6ZuDsW5rfpLaJX266Dgaqfv7Jb0PDygyJrNGGlF0AGYV96GIuKXoIJqdpBER8XrRcQyl9SHn9SFGs6rxGTuzJiRpR0m/k/SCpOck/SK135G6/Cmd4ftY98ud6QzE1yTdJ+llSRdLapF0k6SXJN0iaWyu/y8l/SVt6w5Ju6b2KcAngBPTtn6d2reRdI2kZyU9Lunfcusalc7yLZX0IPDOPvLcVdIMSUskPSPp1NR+hqSrJV0h6UXgGEkbSvqepKfT8D1JG6b+W0q6QVJnWtfvJQ1L806S1JFynydp/x5iuVTSRSmel9L+3y43f19J96T9dI+kfVP7JElzc/1mSLonN/17SZP7se/WyLlOjIdIejDF1yHpq6m9TdJCSaem18sTkj6RW25DSd+W9Oe0ny+SNCo3/1BJcyS9KOlRSQdLOhvYD/hB/mxyOmN8vKT5wPzU9n1JT6XlZ0var7fj3i2nr0lalI7pp7rN6zHuXM4nSFqc1nFsf/dVGv8p8Cbg1ynHEyX9r6QvdovjPkmH9Tcns0JFhAcPHgoYgCeAA3qYdyXwdbJ/vjYC3p2bF8COuek2YGG39c4EWoBxwGLgj8AeaV23Aafn+n8K2BTYEPgeMCc371Lgm7npYcBs4N+BkcAOwGPAQWn+ucDvgS2A8cD9+di65bgpsAg4IcW1KbB3mncGsAKYnLY5Cjgz5fVGYCvg/4CzUv9zgIuADdKwHyBgZ+ApYJvUbwLw5h7iuRR4CXhP2hffB/6Q5m0BLAWOJLvScUSafkOK7a/AlmnbzwAdKZ9RwPLUr699t0bOdWJcBOyXxscC78i9Bl4Hvptify/wMrBzmn8+MD3lsSnwa+CcNG8v4AXg/Wm744C3pHntwKe7xRDAjLSuUantkynHEel4/gXYKJfXFT3s84PT/toN2AT4ObnXdx9x13I+M+33Q4BXgLH92Ffd3y8H5KY/CtyVm/574HlgZNG/Mzx46M9QeAAePFR1SH9QlgGdueEzad7lwFRg2zrL9aew+0Ru+hrgwtz0F4HreohpTFr/5mn6UlYv7PYG/txtmVOA/0njjwEH5+ZNoefC7gjg3h7mnQHc0a3tUeCQ3PRBwBNp/Ezg+vx+Se07khW2BwAb9HE8LgWm5aZHA11kBeqRwN3d+t8JHJPGfw/8E7AP8FvgKrKiZRJwXz/33Ro514nxz8Bngc26tbeRFTmb5NquAr5BVuC+TK6gBf4BeDyN/wg4v4fttVO/sHtfH3EuBf4+l1dPhd0lwLm56Z3S+nfsR9xtZEXziNz8xcA+/dhXvRV2G6X4J6bpbwM/7Ov97MFDswy+FGtWrMkRMSY3/Di1n0j2h+1uSQ90v0TVD8/kxpfXmR4NIGm4pHPT5bcXyf7IQXb2qZ7tgG3SJc9OSZ3AqWRnBwG2ITtDVvNkLzGOJyvWevJUt+ltuq3vydQGcB6wAPitpMcknQwQEQuAL5MVF4slTZO0DT1btc2IWAYsSdvovu3a9sel8d+RFQzvSePtZGfN3pumoe99Vy/n7v6Z7MzUk+lS8T/k5i2NiJe7xbcN2dnNjYHZue3+JrVD38ehntXilPRVSQ+ly9SdwOb0/BrK6+310lfcAM/H6vf4vUJ6bdP7vupRRPwV+AXwSWWX848AftqfZc2agQs7syYUEX+JiM9ExDZkZx1+qMH5CIiPA4eSndHanOxSJWRFJWRnT/KeIjtjki9GN42IQ9L8RWSFQs2betn2U2SXI3vSfdtPkxVH+XU/DRARL0XECRGxA/Bh4Cu1e+ki4ucR8e60bADf6mWbq2KXNJrsEuDTdbZd235HGu9e2P2ONQu7vvZdvZxXExH3RMShZJejryM7K1czVtIm3eJ7GniOrJjfNbfdzSOiVgA9Bby5p0321Z7upzuR7BLm2IgYQ3ZpVz0sm9fb66WvuHvVx76qm0vOZWT3l+4PvBIRd/Znm2bNwIWdWROS9C+Stk2TS8n++KxM08/Qe0E0EJsCr5LdQ7Qx8J/d5nff1t3AS8oeSBiVzvjtJqn2kMRVwCmSxqb4v0jPbgC2lvTldJP8ppL27qX/lcBpkraStCXZvWpXAEj6oLIHTkRWVHQBKyXtLOl9yh6y+CtZobCyh/UDHCLp3ZJGAmcBMyPiKeBGYCdJH5c0QtLHgF1SDpDd77cz2f1qd0fEA2SF4N5A7YGXvvZdrySNlPQJSZtHxArgxTq5/Efqtx/wQeCXEbES+DFwvqQ3pnWNk3RQWuZi4FhJ+0salua9Jc3rz2ttU7LLwM8CIyT9O7BZf3Iie70cI2kXSRsDp9dm9CPuHvVzX9WskWMq5FYC38Fn62w948LOrFi1p/Fqw7Wp/Z3AXZKWkd08/qWIeCzNOwO4LF2e+ug6bv9ysstfHcCDZA8n5F0M7JK2dV1EdJEVDLsDj5OdVfkJ2dk+gP9I63uc7F6zHv8oRsRLZDfsf4jsZvv5ZPek9eSbwCzgPmAu2QMhtc/YmwjcQnbP4p1k90TdTvYgwbkpzr+Qnb05pZdt/JysuFgC7En2UAAR8XzK+wSyIvhE4IMR8Vya/3KK54GIeC2t607gyYhYnPr0te/640jgiXTZ/HNkZ5Vq/kL2T8DTwM+Az0XEw2neSWSXqmemZW8hK0SJiLuBY8keVHiB7Axj7ezk94GPKHvK+YIeYrqZ7BLpI2TH/q/0fUmZtO2byB7YuS3Fd1u3Lj3G3Q+97au8c8j+YeisPTmbXA68jfTPg9n6QhG9nvk3M6sESZeS3VR/WtGxDJSkNrIHFLbtq6/1j6SjgCnpMr7ZesNn7MzMzHLSZeHPkz2ZbrZecWFnZmaWpHv4niW79+7nBYdjNmC+FGtmZmZWEj5jZ2ZmZlYSI4oOoChbbrllTJgwYVC38fLLL7PJJpv03bGkqpx/lXOHaudf5dyh2vlXOXeodv5Dkfvs2bOfi4it+upX2cJuwoQJzJo1a1C30d7eTltb26Buo5lVOf8q5w7Vzr/KuUO1869y7lDt/Icid0m9fZPPKr4Ua2ZmZlYSLuzMzMzMSsKFnZmZmVlJuLAzMzMzKwkXdmZmZmYlUdmnYs3MGum6ezs47+Z5dHQuZ9zM2/jaQTszeY9xRYdlZhXjws7MbB1dd28Hp/xqLstXdAHQ0bmcU341F8DFnZkNqUIvxUoaL+l2SQ9KekDSl1L7FpJmSJqffo5N7ZJ0gaQFku6T9I7cuo5O/edLOrqonMyses67ed6qoq5m+Youzrt5XkERmVlVFX2P3evACRGxC7APcLykXYCTgVsjYiJwa5oG+AAwMQ1TgAshKwSB04G9gb2A02vFoJnZYHu6c/mA2s3MBkuhhV1ELIqIP6bxl4CHgHHAocBlqdtlwOQ0fihweWRmAmMkbQ0cBMyIiCURsRSYARw8hKmYWYVtM2bUgNrNzAaLIqLoGACQNAG4A9gN+HNEjEntApZGxBhJNwDnRsQf0rxbgZOANmCjiPhmav8GsDwivt1tG1PIzvTR0tKy57Rp0wY1p2XLljF69OhB3UYzq3L+Vc4dqpf//z29gkvvf43XVv6tbeQwOGa3key7zQbFBVaAqh37vCrnDtXOfyhynzRp0uyIaO2rX1M8PCFpNHAN8OWIeDGr5TIREZIaUn1GxFRgKkBra2sM9ve6Vfl786Da+Vc5d6he/m3ALvd2cOLV9/Fa10rGjRlV2adiq3bs86qcO1Q7/2bKveh77JC0AVlR97OI+FVqfiZdYiX9XJzaO4DxucW3TW09tZuZDYnJe4xjjzeNYeexw/h/J7+vkkWdmRWv6KdiBVwMPBQR383Nmg7Unmw9Grg+135Uejp2H+CFiFgE3AwcKGlsemjiwNRmZmZmVhlFX4p9F3AkMFfSnNR2KnAucJWk44AngY+meTcChwALgFeAYwEiYomks4B7Ur8zI2LJ0KRgZmZm1hwKLezSQxDqYfb+dfoHcHwP67oEuKRx0ZmZmZmtXwq/x87MzMzMGsOFnZmZmVlJuLAzMzMzKwkXdmZmZmYl4cLOzMzMrCRc2JmZmZmVhAs7MzMzs5JwYWdmZmZWEi7szMzMzErChZ2ZmZlZSbiwMzMzMysJF3ZmZmZmJeHCzszMzKwkXNiZmZmZlYQLOzMzM7OSKLSwk3SJpMWS7s+1/ULSnDQ8IWlOap8gaXlu3kW5ZfaUNFfSAkkXSFIR+ZiZmZkVaUTB278U+AFwea0hIj5WG5f0HeCFXP9HI2L3Ouu5EPgMcBdwI3AwcNMgxGtmZmbWtAo9YxcRdwBL6s1LZ90+ClzZ2zokbQ1sFhEzIyLIisTJjY7VzMzMrNk18z12+wHPRMT8XNv2ku6V9DtJ+6W2ccDCXJ+Fqc3MzMysUoq+FNubI1j9bN0i4E0R8bykPYHrJO06kBVKmgJMAWhpaaG9vb1Rsda1bNmyQd9GM6ty/lXOHaqbf2fncrq6uiqZe01Vjz1UO3eodv7NlHtTFnaSRgD/BOxZa4uIV4FX0/hsSY8COwEdwLa5xbdNbWuIiKnAVIDW1tZoa2sbjPBXaW9vZ7C30cyqnH+Vc4fq5n/hvDvp7OysZO41VT32UO3codr5N1PuzXop9gDg4YhYdYlV0laShqfxHYCJwGMRsQh4UdI+6b68o4DriwjazMzMrEhFf9zJlcCdwM6SFko6Ls06nDUfmngPcF/6+JOrgc9FRO3Bi88DPwEWAI/iJ2LNzMysggq9FBsRR/TQfkydtmuAa3roPwvYraHBmZmZma1nmvVSrJmZmZkNkAs7MzMzs5JwYWdmZmZWEi7szMzMzErChZ2ZmZlZSbiwMzMzMysJF3ZmZmZmJeHCzszMzKwkXNiZmZmZlYQLOzMzM7OScGFnZmZmVhIu7MzMzMxKwoWdmZmZWUm4sDMzMzMrCRd2ZmZmZiXhws7MzMysJAot7CRdImmxpPtzbWdI6pA0Jw2H5OadImmBpHmSDsq1H5zaFkg6eajzMDMzM2sGRZ+xuxQ4uE77+RGxexpuBJC0C3A4sGta5oeShksaDvw38AFgF+CI1NfMzMysUkYUufGIuEPShH52PxSYFhGvAo9LWgDsleYtiIjHACRNS30fbHC4ZmZmZk2t0MKuF1+QdBQwCzghIpYC44CZuT4LUxvAU93a9663UklTgCkALS0ttLe3Nzjs1S1btmzQt9HMqpx/lXOH6ubf2bmcrq6uSuZeU9VjD9XOHaqdfzPl3oyF3YXAWUCkn98BPtWIFUfEVGAqQGtra7S1tTVitT1qb29nsLfRzKqcf5Vzh+rmf+G8O+ns7Kxk7jVVPfZQ7dyh2vk3U+5NV9hFxDO1cUk/Bm5Ikx3A+FzXbVMbvbSbmZmZVUbRD0+sQdLWucnDgNoTs9OBwyVtKGl7YCJwN3APMFHS9pJGkj1gMX0oYzYzMzNrBoWesZN0JdAGbClpIXA60CZpd7JLsU8AnwWIiAckXUX2UMTrwPER0ZXW8wXgZmA4cElEPDDEqZiZmZkVruinYo+o03xxL/3PBs6u034jcGMDQzMzMzNb7zTdpVgzMzMzWzsu7MzMzMxKwoWdmZmZWUm4sDMzMzMrCRd2ZmZmZiXhws7MzMysJFzYmZmZmZWECzszMzOzknBhZ2ZmZlYSLuzMzMzMSsKFnZmZmVlJuLAzMzMzKwkXdmZmZmYl4cLOzMzMrCRc2JmZmZmVRKGFnaRLJC2WdH+u7TxJD0u6T9K1ksak9gmSlkuak4aLcsvsKWmupAWSLpCkIvIxMzMzK1LRZ+wuBQ7u1jYD2C0i3g48ApySm/doROyehs/l2i8EPgNMTEP3dZqZmZmVXqGFXUTcASzp1vbbiHg9Tc4Etu1tHZK2BjaLiJkREcDlwOTBiNfMzMysmSmrhQoMQJoA3BARu9WZ92vgFxFxRer3ANlZvBeB0yLi95JagXMj4oC0zH7ASRHxwTrrmwJMAWhpadlz2rRpg5JTzbJlyxg9evSgbqOZVTn/KucO1c3/nLuW09XVxWn7Vi/3mqoee6h27lDt/Ici90mTJs2OiNa++o0Y1CjWgaSvA68DP0tNi4A3RcTzkvYErpO060DWGRFTgakAra2t0dbW1sCI19Te3s5gb6OZVTn/KucO1c3/wnl30tnZWcnca6p67KHauUO182+m3JuysJN0DPBBYP90eZWIeBV4NY3PlvQosBPQweqXa7dNbWZmZmaVUvTDE2uQdDBwIvDhiHgl176VpOFpfAeyhyQei4hFwIuS9klPwx4FXF9A6GZmZmaFKvSMnaQrgTZgS0kLgdPJnoLdEJiRPrVkZnoC9j3AmZJWACuBz0VE7cGLz5M9YTsKuCkNZmZmZpVSaGEXEUfUab64h77XANf0MG8WsMbDF2ZmZmZV0nSXYs3MzMxs7axzYSdpk9q9b2ZmZmZWnAEXdpKGSfq4pP+VtBh4GFgk6cH0dWA7Nj5MMzMzM+vL2pyxux14M9lDDn8XEeMj4o3Au8m+KeJbkj7ZwBjNzMzMrB/W5uGJAyJiBYCkTYGXANITqtcA10jaoHEhmpmZmVl/DKiwk7QvcLKkTrKzfWPIPkh4NbXCz8zMzMyGzkDP2H0eOCwiuiRtCdwyCDGZmZmZ2VoY6D12L0VEF0BEPAd8qfEhmZmZmdnaGGhhJ0n7S9pS0ueBaYMRlJmZmZkN3EALuy8BbwO+DjwFnNXwiMzMzMxsrQz0HrtTyb7H9Q/AfRHx68aHZGZmZmZrY0CFXUScLmkr4J3AMZK2B+ZHxNmDEp2ZmZmZ9duAP8cuIp4FbkwD6V47MzMzMyvY2nxAMQCSJgInA39tXDhmZmZmtrbW5ivFan4KXA3sByBpN0mXNyQqMzMzMxuwdSnshkXETUDtc+3uB3ZrSFRmZmZmNmDrUtg9nR6eCMg+4A4YNZAVSLpE0mJJ9+fatpA0Q9L89HNsbf2SLpC0QNJ9kt6RW+bo1H++pKPXISczMzOz9da6FHZfBn4M/J2kY8k+rPj+3hdZw6XAwd3aTgZujYiJwK1pGuADwMQ0TAEuhKwQBE4H9gb2Ak6vFYNmZmZmVbLWhV1EPEFWlP0bsAPwO+DIAa7jDmBJt+ZDgcvS+GXA5Fz75ZGZCYyRtDVwEDAjIpZExFJgBmsWi2ZmZmalp4gY2AKSoo+F+tMn13cCcENE7JamOyNiTG09wNKIGCPpBuDciPhDmncrcBLQBmwUEd9M7d8AlkfEt+tsawrZ2T5aWlr2nDZtcL8RbdmyZYwePXpQt9HMqpx/lXOH6uZ/zl3L6erq4rR9q5d7TVWPPVQ7d6h2/kOR+6RJk2ZHRGtf/dbm405ul3QNcH1E/LnWKGkk8G7gaOB2ssus6yQiQtLAKs/e1zcVmArQ2toabW1tjVp1Xe3t7Qz2NppZlfOvcu5Q3fwvnHcnnZ2dlcy9pqrHHqqdO1Q7/2bKfW0uxR5M9iTslZKelvSgpMeB+cARwPci4tJ1iOmZdImV9HNxau8Axuf6bZvaemo3MzMzq5QBF3YR8deI+GFEvAvYDtgf2APYOSI+ExH3rmNM08nO+pF+Xp9rPyo9HbsP8EJELAJuBg6UNDY9NHFgajMzMzOrlLX+5gmAiFgBLAKQdLGkaRExo7/LS7qS7B65LSUtJHu69VzgKknHAU8CH03dbwQOARYArwDHphiWSDoLuCf1OzMiuj+QYWZmZlZ661TY5UXEcZKOk3QBWXH1XD+WOaKHWfvX6RvA8T2s5xLgkoHEa2ZmZlY26/I5dquRdBCwPbAj8GNJk/tYxMzMzMwaqGGFHbA1cElEHBIRhwGTGrhuMzMzM+tDwy7FAj+NiK7c9DcauG4zMzMz60Mjz9j9SNLGAJLeExEvNnDdZmZmZtaHRp6xOx24WNLrwBzgjgau28zMzMz60MgzdmcB84AArmrges3MzMysH9b6jF36HteJgIBHgBMj4jlJmwDfBz7dmBDNzMzMrD/WqrCTtAtwNbAyt55/Ap6LiJclfbZB8ZmZmZlZP/V5KVbSrpJ+1q35f4BTI2K3iNgNOAlY1afb07FmZmZmNgT6c8buFuAfahOS/hUYB2wm6ahaM/AGSZ+PiB82PkwzMzMz60t/CrsDgbOBT6TpPYENgHd26zcyzTMzMzOzAvRZ2EXEXP5W1BERn5b0KPD9iFgAIGkH4IMRcdygRWpmZmZmvVrbp2K/AvxO0rVpejLwhcaEZGZmZmZrY60Ku4i4XtJcssu0AN+rnb0zMzMzs2Ks9efYRcRjwEUNjMXMzMzM1kEjv3miYSTtLGlObnhR0pclnSGpI9d+SG6ZUyQtkDRP0kFFxm9mZmZWhEZ+V2zDRMQ8YHcAScOBDuBa4Fjg/Ij4dr5/+sDkw4FdgW2AWyTt5M/TMzMzsyppyjN23ewPPBoRT/bS51BgWkS8GhGPAwuAvYYkOjMzM7MmoYgoOoZeSboE+GNE/EDSGcAxwIvALOCEiFgq6QfAzIi4Ii1zMXBTRFzdbV1TgCkALS0te06bNm1QY1+2bBmjR48e1G00syrnX+Xcobr5n3PXcrq6ujht3+rlXlPVYw/Vzh2qnf9Q5D5p0qTZEdHaV7+mvBRbI2kk8GHglNR0IXAWEOnnd4BP9Xd9ETEVmArQ2toabW1tjQx3De3t7Qz2NppZlfOvcu5Q3fwvnHcnnZ2dlcy9pqrHHqqdO1Q7/2bKvdkvxX6A7GzdMwAR8UxEdEXESuDH/O1yawcwPrfctqnNzMzMrDKavbA7AriyNiFp69y8w4D70/h04HBJG0raHpgI3D1kUZqZmZk1gaa9FCtpE+D9wGdzzf8laXeyS7FP1OZFxAOSrgIeBF4HjvcTsWZmZlY1TVvYRcTLwBu6tR3ZS/+zgbMHOy4zMzOzZtXsl2LNzMzMrJ9c2JmZmZmVhAs7MzMzs5JwYWdmZmZWEi7szMzMzErChZ2ZmZlZSbiwMzMzMysJF3ZmZmZmJeHCzszMzKwkXNiZmZmZlYQLOzMzM7OScGFnZmZmVhIu7MzMzMxKwoWdmZmZWUm4sDMzMzMrCRd2ZmZmZiXRtIWdpCckzZU0R9Ks1LaFpBmS5qefY1O7JF0gaYGk+yS9o9jozczMzIZe0xZ2yaSI2D0iWtP0ycCtETERuDVNA3wAmJiGKcCFQx6pmZmZWcGavbDr7lDgsjR+GTA51355ZGYCYyRtXUSAZmZmZkVRRBQdQ12SHgeWAgH8KCKmSuqMiDFpvoClETFG0g3AuRHxhzTvVuCkiJjVbZ1TyM7o0dLSsue0adMGNYdly5YxevToQd1GM6ty/lXOHaqb/zl3Laerq4vT9q1e7jVVPfZQ7dyh2vkPRe6TJk2anbuC2aMRgxrFunl3RHRIeiMwQ9LD+ZkREZIGVJVGxFRgKkBra2u0tbU1LNh62tvbGextNLMq51/l3KG6+V847046OzsrmXtNVY89VDt3qHb+zZR7016KjYiO9HMxcC2wF/BM7RJr+rk4de8AxucW3za1mZmZmVVGUxZ2kjaRtGltHDgQuB+YDhyduh0NXJ/GpwNHpadj9wFeiIhFQxy2mZmZWaGa9VJsC3BtdhsdI4CfR8RvJN0DXCXpOOBJ4KOp/43AIcAC4BXg2KEP2czMzKxYTVnYRcRjwN/XaX8e2L9OewDHD0FoZmZmZk2rKS/FmpmZmdnAubAzMzMzKwkXdmZmZmYl4cLOzMzMrCRc2JmZmZmVhAs7MzMzs5JwYWdmZmZWEi7szMzMzErChZ2ZmZlZSbiwMzMzMysJF3ZmZmZmJeHCzszMzKwkXNiZmZmZlYQLOzMzM7OScGFnZmZmVhJNWdhJGi/pdkkPSnpA0pdS+xmSOiTNScMhuWVOkbRA0jxJBxUXvZmZmVkxRhQdQA9eB06IiD9K2hSYLWlGmnd+RHw731nSLsDhwK7ANsAtknaKiK4hjdrMzMysQE15xi4iFkXEH9P4S8BDwLheFjkUmBYRr0bE48ACYK/Bj9TMzMyseSgiio6hV5ImAHcAuwFfAY4BXgRmkZ3VWyrpB8DMiLgiLXMxcFNEXN1tXVOAKQAtLS17Tps2bVBjX7ZsGaNHjx7UbTSzKudf5dyhuvmfc9dyurq6OG3f6uVeU9VjD9XOHaqd/1DkPmnSpNkR0dpXv2a9FAuApNHANcCXI+JFSRcCZwGRfn4H+FR/1xcRU4GpAK2trdHW1tbwmPPa29sZ7G00syrnX+Xcobr5XzjvTjo7OyuZe01Vjz1UO3eodv7NlHtTXooFkLQBWVH3s4j4FUBEPBMRXRGxEvgxf7vc2gGMzy2+bWozMzMzq4ymLOwkCbgYeCgivptr3zrX7TDg/jQ+HThc0oaStgcmAncPVbxmZmZmzaBZL8W+CzgSmCtpTmo7FThC0u5kl2KfAD4LEBEPSLoKeJDsidrj/USsmZmZVU1TFnYR8QdAdWbd2MsyZwNnD1pQZmZmZk2uKS/FmpmZmdnAubAzMzMzKwkXdmZmZmYl4cLOzMzMrCRc2JmZmZmVhAs7MzMzs5JwYWdmZmZWEi7szMzMzErChZ2ZmZlZSbiwMzMzMysJF3ZmZmZmJeHCzszMzKwkXNiZmZmZlYQLOzMzM7OScGFnZraOrru3g93/47fc9fgS5i1dyR5n/pbr7u0oOiwzqyBFRNExNIykg4HvA8OBn0T4Hb71AAALqklEQVTEuT31bW1tjVmzZg1aLNfd28HXr5nDy69n02M33oB/fPvW3P7wszzduZxtxoziawftzOQ9xnHdvR2cd/O8NdrrrbPWb/NRGyBB5ysrVi0D1F3Pdfd2cMb0B+hcvmJVLKd/aNe6/YHV+ubVW27zURvw2utdvLJi5Wp9Ju8xjvb2djo3n8h5N8+jo3M5wyW6IhhXJ956+fS0D+rFl99uvf3VfZ3X3dvBKb+6j+UpbgABAavF131bG2+Q/S/UPd9ZTy7hyrueoiuC4RI7jYUXV25ER+fyVcsOlzhi7/F8c/LbesxnmGBlZDFMestWXDN74WoxjhwuNtlwBJ2vrOj3Psvvi47O5avy7K7ePsw77bq5/Gzmn9dYtpZX63Zb9JhL/njn90mt38f3ftNq+6Xevukpvu7HedJbtuL2h5+t+5rrKbe+nHbd3FXHt7+GDxMjh2vV8RvIa7ReO/T8vtt4g2FsuMFwlr6yYrWca/uip/dq/ndS7fW09JUVq71GejqOPf2++s+fz+Cq+bHG62DUBsN6fS2//OoKcrMBVjtuff3+675fe/o91t2YXN79fb3k30/5Zf7xTV2c+vH3r9Y3/9rJ/w6od4xnPblktfdYfh/V+x22rr9b+3qdre3vldrrsd7vqfxrsrd4BnqM89vsz9+//Guq+++k/hg5XLzWFavW+S87ao1j32iSZkdEa5/9ylLYSRoOPAK8H1gI3AMcEREP1us/mIXddfd28JVfzGFl313NzMysBDYcMYxv/fPb1/ofyL70t7Ar06XYvYAFEfFYRLwGTAMOLSKQM6Y/4KLOzMysQl59fSUn/PJPhd+GMaLQrTfWOOCp3PRCYO8iAqmd8v3sfdezwwu+z8bMzKzMHtt8HD96+6F0rQzOu3neoJ21648yFXZ9kjQFmALQ0tJCe3t7sQGZmZlZqXR0Li+0vihTYdcBjM9Nb5vaVomIqcBUyO6xa2trG5RAxt7xW5a+soIfvb2QK8FmZmZWkHFjRjFY9UV/lOkeu3uAiZK2lzQSOByYXkQgp39oV4apiC2bmZlZUYYP06qni4tSmjN2EfG6pC8AN5N93MklEfFAEbHUrq3nP+7EzMzMymuTkcM5+7C3FXp/HZSosAOIiBuBG4uOA7LibswL8ws9HVu09vb2yuZf5dyh2vlXOXeodv5Vzh2qnX8z5V6mS7FmZmZmlebCzszMzKwkXNiZmZmZlYQLOzMzM7OScGFnZmZmVhKKiKJjKISkZ4EnB3kzWwLPDfI2mlmV869y7lDt/KucO1Q7/yrnDtXOfyhy3y4ituqrU2ULu6EgaVZEtBYdR1GqnH+Vc4dq51/l3KHa+Vc5d6h2/s2Uuy/FmpmZmZWECzszMzOzknBhN7imFh1Awaqcf5Vzh2rnX+Xcodr5Vzl3qHb+TZO777EzMzMzKwmfsTMzMzMrCRd2ZmZmZiXhwm4dSfoXSQ9IWimpx0edJR0saZ6kBZJOzrVvL+mu1P4LSSOHJvLGkLSFpBmS5qefY+v0mSRpTm74q6TJad6lkh7Pzdt96LNYO/3JPfXryuU3PddehWO/u6Q703vkPkkfy81b7459T+/j3PwN07FckI7thNy8U1L7PEkHDWXcjdCP3L8i6cF0nG+VtF1uXt33wPqkH/kfI+nZXJ6fzs07Or1P5ks6emgjX3f9yP38XN6PSOrMzVuvj72kSyQtlnR/D/Ml6YK0b+6T9I7cvGKOe0R4WIcBeCuwM9AOtPbQZzjwKLADMBL4E7BLmncVcHgavwj416JzGmD+/wWcnMZPBr7VR/8tgCXAxmn6UuAjRecxmLkDy3poL/2xB3YCJqbxbYBFwJj18dj39j7O9fk8cFEaPxz4RRrfJfXfENg+rWd40Tk1OPdJuff1v9ZyT9N13wPry9DP/I8BflBn2S2Ax9LPsWl8bNE5NTL3bv2/CFxSomP/HuAdwP09zD8EuAkQsA9wV9HH3Wfs1lFEPBQR8/rothewICIei4jXgGnAoZIEvA+4OvW7DJg8eNEOikPJ4ob+xf8R4KaIeGVQoxoaA819laoc+4h4JCLmp/GngcVAn5+c3qTqvo+79cnvk6uB/dOxPhSYFhGvRsTjwIK0vvVFn7lHxO259/VMYNshjnEw9efY9+QgYEZELImIpcAM4OBBinMwDDT3I4ArhySyIRARd5CdjOjJocDlkZkJjJG0NQUedxd2Q2Mc8FRuemFqewPQGRGvd2tfn7RExKI0/hegpY/+h7Pmm/7sdAr7fEkbNjzCwdPf3DeSNEvSzNolaCp47CXtRfYf/6O55vXp2Pf0Pq7bJx3bF8iOdX+WbWYDjf84srMYNfXeA+uT/ub/z+n1fLWk8QNctln1O/50+X174LZc8/p+7PvS0/4p7LiPGIqNrO8k3QL8XZ1ZX4+I64c6nqHWW/75iYgIST1+fk76L+ZtwM255lPIioKRZJ8DdBJw5rrG3CgNyn27iOiQtANwm6S5ZH/wm16Dj/1PgaMjYmVqbupjb2tH0ieBVuC9ueY13gMR8Wj9Nay3fg1cGRGvSvos2Znb9xUc01A7HLg6IrpybVU49k3FhV0/RMQB67iKDmB8bnrb1PY82WnbEem/+1p7U+ktf0nPSNo6IhalP96Le1nVR4FrI2JFbt21Mz6vSvof4KsNCbpBGpF7RHSkn49Jagf2AK6hIsde0mbA/5L9IzQzt+6mPvZ19PQ+rtdnoaQRwOZk7/P+LNvM+hW/pAPIiv73RsSrtfYe3gPr0x/3PvOPiOdzkz8huwe1tmxbt2XbGx7h4BnIa/dw4Ph8QwmOfV962j+FHXdfih0a9wATlT0FOZLsxT89sjssbye77wzgaGB9OwM4nSxu6Dv+Ne69SAVB7Z6zyUDdJ4+aVJ+5Sxpbu8QoaUvgXcCDVTn26fV+Ldk9KFd3m7e+Hfu67+NuffL75CPAbelYTwcOV/bU7PbARODuIYq7EfrMXdIewI+AD0fE4lx73ffAkEXeGP3Jf+vc5IeBh9L4zcCBaT+MBQ5k9asWza4/r3skvYXsIYE7c21lOPZ9mQ4clZ6O3Qd4If3TWtxxH4onNMo8AIeRXTt/FXgGuDm1bwPcmOt3CPAI2X8qX8+170D2C34B8Etgw6JzGmD+bwBuBeYDtwBbpPZW4Ce5fhPI/oMZ1m3524C5ZH/UrwBGF51TI3MH9k35/Sn9PK5Kxx74JLACmJMbdl9fj3299zHZ5eMPp/GN0rFckI7tDrllv56Wmwd8oOhcBiH3W9LvwNpxnp7ae3wPrE9DP/I/B3gg5Xk78Jbcsp9Kr4kFwLFF59Lo3NP0GcC53ZZb74892cmIRen32EKy+0c/B3wuzRfw32nfzCX36RhFHXd/pZiZmZlZSfhSrJmZmVlJuLAzMzMzKwkXdmZmZmYl4cLOzMzMrCRc2JmZmZmVhAs7MzMzs5JwYWdmZmZWEi7szMwGSNIESTdKmifpEUmndJt/kaR3FRWfmVWXCzszswGQNIzsu34vioidgbcBrZKm5LrtA8yst7yZ2WByYWdmNjAHAU9ExHSAyL7s/gvAVwEkvRV4JCK6JL1Z0rOSnpA0R9ISSY9K2qy48M2szFzYmZkNzFvJvvtylci+9Huz9CXpHwB+k9ofBf4AHBkRuwP3AZMj4sWhDdnMqsKFnZnZwHQBo/MNkgRsDLxOdkbvN7nZuwL3p/G3AvOGIEYzqygXdmZmA9MOHJKKuZr3A38ENgLGRMTTAJJGARtFxFJJ44HnIuK1oQ7YzKrDhZ2Z2QBExJ+Ae4EzASS1AN8FTgUmAbfnuu8CPJTG35obNzMbFCOKDsDMbH0i6WSgFfikpNvJHpzYDvghsBD4dq57/jLscuAdkt4SEQ8PYchmViGKiKJjMDMrBUl/BPaOiBVFx2Jm1eTCzszMzKwkfI+dmZmZWUm4sDMzMzMrCRd2ZmZmZiXhws7MzMysJFzYmZmZmZWECzszMzOzknBhZ2ZmZlYS/x+idL7nkWKIxQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 720x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "N = 1024  # length of random signal x\n",
    "\n",
    "# generate two random signals\n",
    "np.random.seed(2)\n",
    "x = 2 + np.random.normal(size=N)\n",
    "y = 1 + np.random.normal(size=2*N)\n",
    "\n",
    "# compute cross PSD via CCF\n",
    "acf = 1/N * np.correlate(x, y, mode='valid')\n",
    "psd = np.fft.fft(acf)\n",
    "psd = psd * np.exp(1j*np.arange(N+1)*2*np.pi*(N-1)/(2*N-1))\n",
    "\n",
    "# plot results\n",
    "f = np.fft.fftfreq(len(psd), d=1/2)\n",
    "plt.figure(figsize = (10, 4))\n",
    "plt.stem(f, np.real(psd))\n",
    "plt.title('Estimated cross power spectral density')\n",
    "plt.ylabel(r'$\\hat{\\Phi}_{xy}(e^{j \\Omega})$')\n",
    "plt.xlabel(r'$\\Omega/ \\pi$')\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise**\n",
    "\n",
    "* What does the cross PSD $\\Phi_{xy}(\\mathrm{e}^{\\,\\mathrm{j} \\, \\Omega})$ tell you about the statistical properties of the two random signals?\n",
    "\n",
    "Solution: The cross PSD $\\Phi_{xy}(\\mathrm{e}^{\\,\\mathrm{j} \\, \\Omega})$ is essential only non-zero for $\\Omega=0$. It hence can be concluded that the two random signals are not mean-free and uncorrelated to each other."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbsphinx": "hidden"
   },
   "source": [
    "**Copyright**\n",
    "\n",
    "This notebook is provided as [Open Educational Resource](https://en.wikipedia.org/wiki/Open_educational_resources). Feel free to use the notebook for your own purposes. The text is licensed under [Creative Commons Attribution 4.0](https://creativecommons.org/licenses/by/4.0/), the code of the IPython examples under the [MIT license](https://opensource.org/licenses/MIT). Please attribute the work as follows: *Sascha Spors, Digital Signal Processing - Lecture notes featuring computational examples, 2016-2018*."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}