{ "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": [ "## Auto-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 that the DTFT is performed with respect to $\\kappa$. The ACF of a wide-sense ergodic 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[\\kappa] * x_N[-\\kappa].\n", "\\end{equation}\n", "\n", "Taking the DTFT of the left- and right-hand side yields\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 that the PSD quantifies the squared amplitude or power density of a random signal. This is reflected by 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 is 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}[\\kappa=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 via the inverse DTFT of $\\Phi_{xx}$ and considering that $\\mathrm{e}^{\\mathrm{j} \\Omega \\kappa} = 1$ when evaluating the integral for $\\kappa=0$." ] }, { "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 of length $N$ 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. The speech signal constitutes a recording of the vowel 'o' spoken from a German male, loaded into variable `x`.\n", "\n", "In Python the ACF is stored in a vector with indices $0, 1, \\dots, 2N - 2$ corresponding to the lags $\\kappa = (0, 1, \\dots, 2N - 2)^\\mathrm{T} - (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 $\\mathrm{e}^{\\mathrm{j} \\mu \\frac{2 \\pi}{2N - 1} (N-1)}$." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAA1YAAAGJCAYAAACJq6K4AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAABGN0lEQVR4nO3dd3hUZfrG8XtSSEIIvSXUCChIL4JIF0GquKxYsIDusiIosqxiWwUsRNRlcRcF28K6SpGfoICoBIVIC9JCtYDSBQKEJJBAmCTv7w82swxJYJIzLZnv57py6Zw5c84zw5PJuee85x2bMcYIAAAAAFBsQb4uAAAAAABKOoIVAAAAAFhEsAIAAAAAiwhWAAAAAGARwQoAAAAALCJYAQAAAIBFBCsAAAAAsIhgBQAAAAAWEawAAAAAwCKCFYASa/bs2bLZbIX+rFq1yuVt/fbbb5o4caKSkpLy3Tdx4kTZbDb3FV4Eu3fv1sSJE7V//363b9uXzwuuuVJfukPe71Bx+8tfeqig5zFnzhxNmzbNZzUBCDwhvi4AAKyaNWuWGjdunG/59ddf7/I2fvvtN02aNEn169dXq1atnO774x//qD59+lgts1h2796tSZMmqXv37qpfv75PaoDvXKkv8T/9+/fX+vXrFR0d7Vg2Z84c7dy5U2PHjvVdYQACCsEKQInXrFkztWvXzmPbr127tmrXru2x7cPzcnJylJ2drbCwMF+X4lGZmZkqW7asr8vwumrVqqlatWq+LgNAgGMoIICAsGDBAnXo0EEVKlRQ2bJldc011+ihhx6SJK1atUo33HCDJOnBBx90DCWcOHGipIKHO9WvX18DBgzQ0qVL1bp1a0VERKhJkyZaunSppItDk5o0aaLIyEi1b99emzZtcnr8pk2bdPfdd6t+/fqKiIhQ/fr1dc899+jAgQOOdWbPnq0hQ4ZIknr06OGoa/bs2Y51VqxYoZ49e6p8+fIqW7asOnXqpG+++Sbf8//iiy/UqlUrhYWFKTY2Vm+88YbLr1337t3VrFkzrV69WjfeeKMiIiJUq1YtPf/888rJyXFaNyUlRaNGjVKtWrVUpkwZXXPNNXruueeUlZXlWGfIkCFq2rSp0+MGDhwom82mBQsWOJZt2bJFNptNS5YscSw7duyYHn74YdWuXVtlypRRbGysJk2apOzsbMc6+/fvl81m02uvvaaXX35ZsbGxCgsL08qVKwt9jlfqD+lij9hsNn300UcaN26catasqYiICHXr1k1bt27Nt71NmzbptttuU+XKlRUeHq7WrVvrk08+ybfekSNH9Kc//Ul16tRRmTJlFBMTozvuuEPHjx+/al8OHz5c5cqV044dO9S7d29FRUWpZ8+ekqT4+HgNGjRItWvXVnh4uBo2bKiHH35YJ0+eLPQ1uBpXe8gYo7ffflutWrVSRESEKlWqpDvuuEO//vqr03p5fbVx40Z16dLF8bq/+uqrys3NdayXm5url19+Wdddd50iIiJUsWJFtWjRQm+++aZjncuHAnbv3l1ffPGFDhw44DQ82BijRo0a6dZbb81X99mzZ1WhQgWNHj262K8RgABnAKCEmjVrlpFkEhMTjd1ud/rJzs52rLdu3Tpjs9nM3XffbZYtW2a+/fZbM2vWLHP//fcbY4xJS0tzbOuvf/2rWb9+vVm/fr05dOiQMcaYCRMmmMvfLuvVq2dq165tmjVrZubOnWuWLVtmOnToYEJDQ80LL7xgOnXqZBYuXGgWLVpkrr32WlOjRg2TmZnpePyCBQvMCy+8YBYtWmQSEhLMvHnzTLdu3Uy1atXMiRMnjDHGJCcnm8mTJxtJ5q233nLUlZycbIwx5j//+Y+x2Wzm9ttvNwsXLjRLliwxAwYMMMHBwWbFihWOfa1YscIEBwebzp07m4ULF5oFCxaYG264wdStWzff8ypIt27dTJUqVUxMTIz5xz/+Yb7++mszZswYI8mMHj3asd65c+dMixYtTGRkpHnjjTfM8uXLzfPPP29CQkJMv379HOvNnDnTSDK//fabMcYYu91uoqKiTEREhBkxYoRjvSlTppiQkBCTnp5ujDHm6NGjpk6dOqZevXrmnXfeMStWrDAvvfSSCQsLM8OHD3c8bt++fUaSqVWrlunRo4f5v//7P7N8+XKzb9++Ap/f1frDGGNWrlxpJJk6deqYQYMGmSVLlpiPPvrINGzY0JQvX9788ssvjnW//fZbU6ZMGdOlSxczf/5889VXX5nhw4cbSWbWrFmO9Q4fPmyio6NN1apVzdSpU82KFSvM/PnzzUMPPWR++OGHq/blsGHDTGhoqKlfv76Ji4sz33zzjfn666+NMcbMmDHDxMXFmcWLF5uEhATz73//27Rs2dJcd9115sKFC44a8rZf2GuTpyg9NGLECBMaGmr+8pe/mK+++srMmTPHNG7c2NSoUcMcO3bMsV5eXzVq1MjMnDnTxMfHm1GjRhlJ5t///rdjvbi4OBMcHGwmTJhgvvnmG/PVV1+ZadOmmYkTJxb6PHbt2mU6depkatas6Xjd1q9fb4wx5s033zQ2m838/PPPTnW/9dZbRpLZtWvXFV8LACgMwQpAiZV3MFXQT3BwsGO9N954w0gyqamphW5r48aN+Q588xQWrCIiIszhw4cdy5KSkowkEx0dbTIyMhzLP/vsMyPJLF68uND9Z2dnm7Nnz5rIyEjz5ptvOpYvWLDASDIrV650Wj8jI8NUrlzZDBw40Gl5Tk6OadmypWnfvr1jWYcOHUxMTIw5d+6cY1l6erqpXLmyy8FKkvn888+dlo8YMcIEBQWZAwcOGGP+F5g++eQTp/WmTJliJJnly5cbY4zZu3evkWQ+/PBDY4wxa9asMZLM+PHjTWxsrONxvXr1MjfddJPj9sMPP2zKlSvn2F+evH/fvAPivGDVoEEDpxBRGFf6Iy9YtWnTxuTm5jqW79+/34SGhpo//vGPjmWNGzc2rVu3Nna73WkbAwYMMNHR0SYnJ8cYY8xDDz1kQkNDze7duwvd75X6ctiwYUaS+de//nXF55ebm2vsdrs5cOBAvn9HV4OVqz20fv16I8n87W9/c3r8oUOHTEREhBk/frxjWV5fbdiwwWnd66+/3tx6662O2wMGDDCtWrW6Yn0FPY/+/fubevXq5Vs3PT3dREVFmccffzzffnv06HHF/QDAlQTEUMDvvvtOAwcOVExMjGw2mz777LMiPf78+fMaPny4mjdvrpCQEN1+++0FrpeQkKC2bdsqPDxc11xzjWbOnGm9eABX9eGHH2rjxo1OPxs2bHDcnzec6s4779Qnn3yiI0eOuGW/rVq1Uq1atRy3mzRpIuniMKRLr3PJW37pML+zZ8/qqaeeUsOGDRUSEqKQkBCVK1dOGRkZ+uGHH66673Xr1iklJUXDhg1Tdna24yc3N1d9+vTRxo0blZGRoYyMDG3cuFGDBw9WeHi44/FRUVEaOHCgy881KipKt912m9OyoUOHKjc3V999950k6dtvv1VkZKTuuOMOp/WGDx8uSY4hig0aNFD9+vW1YsUKSReHrTVv3lz33Xef9u3bp19++UVZWVlas2aNbrnlFsd2li5dqh49eigmJsbpOfft21fSxffgS912220KDQ296nMrSn8MHTrUaVhovXr1dNNNNzmGGe7du1c//vij7r33XklyqrNfv346evSofvrpJ0nSl19+qR49ejj6o7h+//vf51uWnJyskSNHqk6dOgoJCVFoaKjq1asnSS7116WK0kNLly6VzWbTfffd5/Tca9asqZYtW+abqbNmzZpq376907IWLVo4/a60b99e27Zt06hRo/T1118rPT29SPVfLioqSg8++KBmz56tjIwMSRd7d/fu3Xr00UctbRtAYAuIYJWRkaGWLVtq+vTpxXp8Tk6OIiIiNGbMGKc/8pfat2+f+vXrpy5dumjr1q169tlnNWbMGH366adWSgfggiZNmqhdu3ZOP23btnXc37VrV3322WfKzs7WAw88oNq1a6tZs2aaO3eupf1WrlzZ6XaZMmWuuPz8+fOOZUOHDtX06dP1xz/+UV9//bW+//57bdy4UdWqVdO5c+euuu/jx49Lku644w6FhoY6/UyZMkXGGKWkpOj06dPKzc1VzZo1822joGWFqVGjRqGPP3XqlOO/NWvWzHc9WvXq1RUSEuJYT5J69uzpCForVqxQr1691Lx5c9WoUUMrVqzQ2rVrde7cOaf33OPHj2vJkiX5nm/e9VqXXz906QxxV1KU/ijsdcx7bnn/Lk888US+OkeNGuVU54kTJyxPilK2bFmVL1/eaVlubq569+6thQsXavz48frmm2/0/fffKzExUZJc6q9LFaWHjh8/LmOMatSoke/5JyYm5vs3qlKlSr5thoWFOdX4zDPP6I033lBiYqL69u2rKlWqqGfPnvmuWyyKxx57TGfOnNHHH38sSZo+fbpq166tQYMGFXubABAQswL27dvX8YlmQS5cuKC//vWv+vjjj5WamqpmzZppypQp6t69uyQpMjJSM2bMkCStXbtWqamp+bYxc+ZM1a1b1/GdGU2aNNGmTZv0xhtvFPhpIgDvGjRokAYNGqSsrCwlJiYqLi5OQ4cOVf369dWxY0ev1pKWlqalS5dqwoQJevrppx3Ls7KylJKS4tI2qlatKkn65z//qRtvvLHAdWrUqCG73S6bzaZjx47lu7+gZYXJCwwFPT7v4LhKlSrasGGDjDFO4So5OVnZ2dmOmqWLweqDDz7Q999/rw0bNuivf/2rJOnmm29WfHy8Dhw4oHLlyjk9t6pVq6pFixZ65ZVXCqwxJibG6XZRvl/J1f4o7HXMew3ynuMzzzyjwYMHF7iv6667TtLFmewOHz7sco0FKeg57ty5U9u2bdPs2bM1bNgwx/K9e/cWax+VKlVyuYeqVq0qm82m1atXFzgDY3FmZQwJCdG4ceM0btw4paamasWKFXr22Wd166236tChQ8WaBbFhw4bq27ev3nrrLfXt21eLFy/WpEmTFBwcXORtAUCegDhjdTUPPvig1q5dq3nz5mn79u0aMmSI+vTpoz179ri8jfXr16t3795Oy2699VZt2rRJdrvd3SUDKKawsDB169ZNU6ZMkSTHjG55B3xF/TS/OPJmJ7v8IPP999/PN8teYXV16tRJFStW1O7du/Odrcv7KVOmjGNWwoULFzqdMTtz5ozTbHtXc+bMGS1evNhp2Zw5cxQUFKSuXbtKuhiWzp49m2+49Ycffui4P0/Pnj1ls9n0/PPPO23jlltu0cqVKxUfH6+uXbs6DeUbMGCAdu7cqQYNGhT4fC8PVsVRWH/kmTt3rowxjtsHDhzQunXrHB/EXXfddWrUqJG2bdtW6L9LVFSUpIsf+q1cudIxNLCweqSi9WVe2Lq8v9555x2Xt3GpovTQgAEDZIzRkSNHCnzuzZs3L1YNeSpWrKg77rhDo0ePVkpKyhW/2PjyM1+Xe/zxx7V9+3YNGzZMwcHBGjFihKXaACAgzlhdyS+//KK5c+fq8OHDjj/KTzzxhL766ivNmjVLkydPdmk7x44dyzdUpkaNGsrOztbJkyddHpICoOh27tzpNN12ngYNGqhatWp64YUXdPjwYfXs2VO1a9dWamqq3nzzTYWGhqpbt26OdSMiIvTxxx+rSZMmKleunGJiYtxysH658uXLq2vXrnr99ddVtWpV1a9fXwkJCfrggw9UsWJFp3WbNWsmSXr33XcVFRWl8PBwxcbGqkqVKvrnP/+pYcOGKSUlRXfccYeqV6+uEydOaNu2bTpx4oTjTPtLL72kPn36qFevXvrLX/6inJwcTZkyRZGRkS6fIatSpYoeeeQRHTx4UNdee62WLVum9957T4888ojq1q0rSXrggQf01ltvadiwYdq/f7+aN2+uNWvWaPLkyerXr5/TsL7q1aurWbNmWr58uXr06OE463DLLbcoJSVFKSkpmjp1qlMNL774ouLj43XTTTdpzJgxuu6663T+/Hnt379fy5Yt08yZM4s1tM6V/siTnJys3/3udxoxYoTS0tI0YcIEhYeH65lnnnGs884776hv37669dZbNXz4cNWqVUspKSn64YcftGXLFseU8i+++KK+/PJLde3aVc8++6yaN2+u1NRUffXVVxo3bpwaN25crL7Me9zTTz8tY4wqV66sJUuWKD4+vsivTR5Xe6hTp07605/+pAcffFCbNm1S165dFRkZqaNHj2rNmjVq3ry5HnnkkSLte+DAgY7vqqtWrZoOHDigadOmqV69emrUqFGhj2vevLkWLlyoGTNmqG3btgoKCnL6vrtevXrp+uuv18qVK3XfffepevXqRX9hAOBSPpw4wyckmUWLFjluf/LJJ0aSiYyMdPoJCQkxd955Z77HDxs2zAwaNCjf8kaNGpnJkyc7Lcub6ero0aPufhoAzJVnBZRk3nvvPWOMMUuXLjV9+/Y1tWrVMmXKlDHVq1c3/fr1M6tXr3ba3ty5c03jxo1NaGiokWQmTJhgjCl8VsD+/fvnq0mXTUFuzP9mqXv99dcdyw4fPmx+//vfm0qVKpmoqCjTp08fs3PnTlOvXj0zbNgwp8dPmzbNxMbGmuDg4HwzxCUkJJj+/fubypUrm9DQUFOrVi3Tv39/s2DBAqdtLF682LRo0cKUKVPG1K1b17z66qsFPq+CdOvWzTRt2tSsWrXKtGvXzoSFhZno6Gjz7LPP5pv57tSpU2bkyJEmOjrahISEmHr16plnnnnGnD9/Pt92//znPxtJ5pVXXnFa3qhRIyPJbN++Pd9jTpw4YcaMGWNiY2NNaGioqVy5smnbtq157rnnzNmzZwt9va/Elf7ImxXwP//5jxkzZoypVq2aCQsLM126dDGbNm3Kt81t27aZO++801SvXt2EhoaamjVrmptvvtnMnDnTab1Dhw6Zhx56yNSsWdOEhoaamJgYc+edd5rjx4871imsL4cNG2YiIyMLfE67d+82vXr1MlFRUaZSpUpmyJAh5uDBg06PN8b1WQGNKVoP/etf/zIdOnQwkZGRJiIiwjRo0MA88MADTq9VXl9dbtiwYU6z+f3tb38zN910k6latapj33/4wx/M/v37r/g8UlJSzB133GEqVqxobDZbgXVOnDjR8ZUNAGCVzZhLxjQEAJvNpkWLFjlm9ps/f77uvfde7dq1K9/Y6nLlyuW7MHf48OFKTU3NN9Sla9euat26tdMXFi5atEh33nmnMjMzXZqZCgD8Uffu3XXy5Ent3LnT16X4zKpVq9SjRw8tWLAg36yHKLnatWsnm82mjRs3+roUAKVAwA8FbN26tXJycpScnKwuXboUezsdO3bMN9Z8+fLlateuHaEKAAA/kZ6erp07d2rp0qXavHmzFi1a5OuSAJQSARGszp496zQb0r59+5SUlKTKlSvr2muv1b333qsHHnhAf/vb39S6dWudPHlS3377rZo3b65+/fpJknbv3q0LFy4oJSVFZ86cUVJSkqSL32MjSSNHjtT06dM1btw4jRgxQuvXr9cHH3xgeTpnAADgPlu2bFGPHj1UpUoVTZgwodDvpgSAogqIoYB5QzguN2zYMM2ePVt2u10vv/yyPvzwQx05ckRVqlRRx44dNWnSJMcMRvXr13f6wsI8l758CQkJ+vOf/6xdu3YpJiZGTz31lEaOHOm5JwYAAADALwREsAIAAAAAT+J7rAAAAADAIoIVAAAAAFhUqievyM3N1W+//aaoqCjHN9EDAAAACDzGGJ05c0YxMTEKCnL/+aVSHax+++031alTx9dlAAAAAPAThw4dUu3atd2+3VIdrKKioiRdnF69cuXKPq4GpZndbtfy5cvVu3dvvrcMHkWvwVvoNXgLvQZvSUlJUWxsrCMjuFupDlZ5w/+ioqJUvnx5H1eD0sxut6ts2bIqX748fxTgUfQavIVeg7fQa/AWu90uSR67RIjJKwAAAADAIoIVAAAAAFhEsAIAAAAAiwhWAAAAAGARwQoAAAAALCJYAQAAAIBFBCsAAAAAsIhgBQAAAAAWEawAAAAAwCKCFQAAAABYRLACAAAAAIsIVgAAAABgEcEKAAAAACwiWAEAAACARQQrAAAAALCIYAUAAAAAFhGsAAAAAMAighUAAAAAWESwAgAAAACLCFYAAAAAYBHBCgAAAAAsIlgBAAAAgEUEKwAAAACwiGAFuNGF7FxtOXhaObnG16UAAADAiwhWgBuNX7hTg99ep7/H/+zrUgAAAOBFBCvAjb7YcUyS9O53v/q4EgAAAHgTwQoAAAAALCJYAQAAAIBFBCsAAAAAsIhgBQAAAAAWEawAAAAAwCKCFQAAAABYRLACAAAAAIsIVgAAAABgEcEKAAAAACwiWAEAAACARX4brLKzs/XXv/5VsbGxioiI0DXXXKMXX3xRubm5vi4NAAAAAJyE+LqAwkyZMkUzZ87Uv//9bzVt2lSbNm3Sgw8+qAoVKujxxx/3dXkAAAAA4OC3wWr9+vUaNGiQ+vfvL0mqX7++5s6dq02bNvm4MgAAAABw5rfBqnPnzpo5c6Z+/vlnXXvttdq2bZvWrFmjadOmFfqYrKwsZWVlOW6np6dLkux2u+x2u6dLRgC7vL+MDD0Hj8jrK/oLnkavwVvoNXiLp3vMb4PVU089pbS0NDVu3FjBwcHKycnRK6+8onvuuafQx8TFxWnSpEn5lq9cuVJly5b1ZLmAk9zcXC1btszXZaAUi4+P93UJCBD0GryFXoOnZWZmenT7fhus5s+fr48++khz5sxR06ZNlZSUpLFjxyomJkbDhg0r8DHPPPOMxo0b57idnp6uOnXqqEePHqpSpYq3SkcAstvtTn8QgoKC1K/frT6sCKVVXq/16tVLoaGhvi4HpRi9Bm+h1+Atp06d8uj2/TZYPfnkk3r66ad19913S5KaN2+uAwcOKC4urtBgFRYWprCwsHzLQ0ND+UWFV9lko+fgUbyvwVvoNXgLvQZP83R/+e1065mZmQoKci4vODiY6dZRIhgZX5cAAAAAL/LbM1YDBw7UK6+8orp166pp06baunWrpk6dqoceesjXpQEAAACAE78NVv/85z/1/PPPa9SoUUpOTlZMTIwefvhhvfDCC74uDQAAAACc+G2wioqK0rRp0644vTrgr2yy+boEAAAAeJHfXmMFlGRcYwUAABBYCFYAAAAAYBHBCgAAAAAsIlgBHsA1VgAAAIGFYAV4ANdYAQAABBaCFQAAAABYRLACAAAAAIsIVoAHcI0VAABAYCFYAR7ANVYAAACBhWAFAAAAABYRrAAPsOcY/d/mw74uAwAAAF5CsAI85IkF23xdAgC41S8nzmrgP9fo613HfF0KAPgdghUAAHDJn+cnaceRND38n82+LgUA/A7BCnCT+b/w6wSgdEs/Z/d1CQDgtzgSBNxkXTK/TgAAAIGKI0EAAAAAsIhgBQAAAAAWEawAAAAAwCKCFQAAAABYRLACAAAAAIsIVgAAAABgEcEKAAAAACwiWAEAAACARQQrAAAAALCIYAUAAAAAFhGsAAAAAMAighUAAAAAWESwAgAAAACLCFYAAAAAYBHBCgAAAAAsIlgBAAAAgEUEKwAAAACwiGAFAAAAABYRrAAAgEuMrwsAAD9GsAIAAAAAiwhWAADAJTZfFwAAfoxgBQAAAAAWEawAAAAAwCKCFQAAAABYRLACAAAAAIsIVgAAAABgEcEKcIOfj5/xdQkAAADwIYIVYNF5e476T1/v6zIAAADgQwQrwKL0c3ZflwAAAAAfI1gBAAAAgEUEKwAAAACwiGAFAAAAABYRrACLjK8LAAAAgM8RrAAAAADAIoIVAAAAAFhEsAIAAAAAiwhWgEU2XxcAAAAAnyNYARYxeQUAAAAIVgAAoMhmr92nnFw+WgKAPAQrAABQZBOX7NYnmw75ugwA8Bt+HayOHDmi++67T1WqVFHZsmXVqlUrbd682ddlAQAQkC4/P7X7t3Sf1AEA/ijE1wUU5vTp0+rUqZN69OihL7/8UtWrV9cvv/yiihUr+ro0AAAAAHDit8FqypQpqlOnjmbNmuVYVr9+fd8VBABAgLt8FlQb06ICgIPfBqvFixfr1ltv1ZAhQ5SQkKBatWpp1KhRGjFiRKGPycrKUlZWluN2evrFIQp2u112u93jNSMwXam36Du4W15P0VvwtIJ6zVw2FjAnJ5dehGW8r8FbPN1jNmMuf5v0D+Hh4ZKkcePGaciQIfr+++81duxYvfPOO3rggQcKfMzEiRM1adKkfMvnzJmjsmXLerReBK60C9ILmwv+jOLNjtlergYAPOelrcE6ef5/p6k618jVkGtyfVgRALguMzNTQ4cOVVpamsqXL+/27fttsCpTpozatWundevWOZaNGTNGGzdu1Pr16wt8TEFnrOrUqaOjR4+qSpUqHq8Zgel4+nl1fv27Au/b81JvL1eD0s5utys+Pl69evVSaGior8tBKVZQr93y9zU6kJLpWOe+DnU0YUATX5WIUoL3NXjLqVOnFB0d7bFg5bdDAaOjo3X99dc7LWvSpIk+/fTTQh8TFhamsLCwfMtDQ0P5RYXHhIbmXOE++g6ewfsavOXSXrv8mqqgoCD6EG7D+xo8zdP95bfTrXfq1Ek//fST07Kff/5Z9erV81FFAAAAAFAwvw1Wf/7zn5WYmKjJkydr7969mjNnjt59912NHj3a16UBAAAAgBO/DVY33HCDFi1apLlz56pZs2Z66aWXNG3aNN17772+Lg0AAAAAnPjtNVaSNGDAAA0YMMDXZQAAAADAFfntGSsAAAAAKCkIVgAAoFhsV18FAAIGwQoAALjEL7/4EgD8BMEKAAAAACwiWAEAAJcw9A8ACkewAiwyjI0BAAAIeAQrAABQLDYb57AAIA/BCgAAAAAsIlgBAAAAgEUEK8AiRsIACFSZF7L18tLd2rQ/xdelAIDPhfi6AKCkY/IKAIHqk02HJUnvr9mn/a/293E1AOBbnLECAAAu4XMkACicW89Y2e12HTt2TJmZmapWrZoqV67szs0DfslwqAEAABDwLJ+xOnv2rN555x11795dFSpUUP369XX99derWrVqqlevnkaMGKGNGze6o1YAAOBDXFIKAIWzFKz+/ve/q379+nrvvfd08803a+HChUpKStJPP/2k9evXa8KECcrOzlavXr3Up08f7dmzx111AwAAAIDfsDQUcN26dVq5cqWaN29e4P3t27fXQw89pJkzZ+qDDz5QQkKCGjVqZGWXgN9h8goAAABYClYLFixwab2wsDCNGjXKyq4AAAAAwG8xKyBgESesAAAAUKxgtWnTJg0dOlRt27ZVhw4dNH78eB07dkzdunVzd30AAAAA4PeKPBRw7ty5+vvf/66XX35ZzZs315kzZ7R06VJ169ZNaWlpnqgRAAAAAPxakYPVK6+8otWrV6tSpUqSpOjoaI0bN079+/fX0KFD3V4g4O8Ms1cAAAAEvCIPBczJyXGEqktdd911+vjjj91SFAAAAACUJEUOVrVq1dKqVavyLZ86daqaNm3qjpqAEoUTVgAAACjyUMCZM2dqyJAhatGihVq0aOG4xqpBgwZ8RxUAAACAgFTkM1YNGzbU5s2bddddd8lms6lSpUp67733NH/+fD344IOeqBEAAPgBTtADQOGKfMZqypQpeuqpp1S7dm316tVLoaGhjvueeuoptxYHAAAAACVBkYNV586dJUkTJ07Ujz/+qJCQEDVr1kzNmzdX8+bNdcMNN6hGjRpuLxQAAPiWzdcFAIAfK3Kw6tSpkyRp4cKFkqTMzEzt3LlTO3bs0IoVKzRhwgT169dPL730knsrBfwUk1cAAACgyMHqcmXLllX79u3Vvn17x7J27doRrAAAAAAEjCJPXlGYtLQ0/elPf1LDhg2Vmpqqo0ePumvTgF8zXM4NAAAQ8NwWrEaNGqUdO3botdde09GjR3Xu3DlJ0tixY/Xmm2+6azcAAAAA4HfcFqy+/PJLvf322xo8eLCCgv632T59+ug///mPu3YDAAAAAH7HbcFKksqVK5dvWaNGjbR371537gbwK0xeAQAAALcFq379+mnOnDn5lp89e1Y2GxO0AgBQ0vE5EgAUzvKsgHni4uLUrl07SZIxRjabTefOndOLL76oNm3auGs3gN/hQAMAAABuC1Z16tTR2rVr9cgjjygzM1Pt27fXmTNnVL58eS1btsxduwEAAD7C+BMAKJzbgpUkNWzYUPHx8Tp48KC2bdum0NBQdejQQZUqVXLnbgC/YrjICgAAIOC5NVjlqVu3rurWreuJTQMAAB/hYyQAKJylySsOHjxYpPWPHDliZXcAAAAA4JcsBasbbrhBI0aM0Pfff1/oOmlpaXrvvffUrFkzLVy40MruAL/EJ7gAAgXXWAFA4SwNBfzhhx80efJk9enTR6GhoWrXrp1iYmIUHh6u06dPa/fu3dq1a5fatWun119/XX379nVX3QAAAADgNyydsapcubLeeOMN/fbbb5oxY4auvfZanTx5Unv27JEk3Xvvvdq8ebPWrl1LqEKpxdwVAAAAcMvkFeHh4Ro8eLAGDx7sjs0BAAAAQIli6YzVlbzzzjue2jQAAAAA+BWPBav169frscceU25uriTpp59+0v333++p3QE+xFhAAACAQOeR77GSpNmzZ2vq1Knq16+fKlasqD179ujJJ5/01O4AAAAAwGc8Fqy2bNmitWvX6vjx4/r555+1cuVK1atXz1O7A3yGySsABAre7gCgcB4bCjhq1Cj94Q9/0NatWzVv3jwNGjRIa9eu9dTuAAAAAMBnPBashg8frn79+kmS2rdvry+++ELjx4/31O4An+ETXACB4kpfEPx50hEZTuEDCGAeC1aJiYlOk1ecPXtW0dHRntodAADwocfnJWnZjmO+LgMAfMZjwWr27NmKjY1Vv379dPfdd2vo0KEaMmSIp3YHAAB8bOvB074uAQB8hskrAIsY+QIAAAAmrwAAAAAAizx2xioxMdHx/3mTVwwZMkTr1q3z1C4BnzBMXwEAABDw3H7GKjU1VadP5x9jXatWLX377bfu3h0AAAAA+JzbgtX69evVunVrValSRVWrVlXLli3znZ0KDw8v9vbj4uJks9k0duxYi5UCAIDi4Pw8ABSuyMFq7ty5+ZYdOHBAvXr1UkhIiOLi4vTqq68qLCxMvXr10r59+ywXuXHjRr377rtq0aKF5W0B7sbkFQAAAHA5WB07dkyDBw9WfHx8vvsmT56sm266SYmJiRo/fryefPJJJSYmqmvXrnrllVcsFXj27Fnde++9eu+991SpUiVL2wIAAMV3pS8IBoBA5/LkFe+++66ys7P1r3/9y2n5wYMHtWrVKo0fP15Hjhxxuu+uu+7SSy+9pEOHDqlOnTrFKnD06NHq37+/brnlFr388stXXDcrK0tZWVmO2+np6ZIku90uu91erP0DV2O3Z1/hPvoO7pXXU/QWPK2gXrvaGfrc3Fx6E0XG+xq8xdM9ZjPGtYFMqampGjNmjDIyMvTpp586lgcFBclmK/wzLGOMbDabcnJyilzcvHnz9Morr2jjxo0KDw9X9+7d1apVK02bNq3A9SdOnKhJkyblWz5nzhyVLVu2yPsHXHEkQ3pte8GfUbzZsfDQBQAlzUtbgnUyq/C/+d2jc/W7+rlerAgAXJeZmamhQ4cqLS1N5cuXd/v2XT5jVbFiRX344YdasmSJ0/ItW7bod7/7nZ544gl17tzZ6b5169ZpypQp+vzzz4tc2KFDh/T4449r+fLlLk968cwzz2jcuHGO2+np6apTp4569OihKlWqFLkGwBW7j6ZL2xMLvK9fv35ergalnd1uV3x8vHr16qXQ0FBfl4NSrKBe+9tPq6Wsc4U+5prYWPXre523SkQpwfsavOXUqVMe3X6Rv8dq4MCBTrdbtWqlLl26aNWqVRo9erTTfXFxceratatatmxZ5MI2b96s5ORktW3b1rEsJydH3333naZPn66srCwFBwc7PSYsLExhYWH5thUaGsovKjwmJKTwXyP6Dp7C+xq85dJeC7rCCBVJCg4Ooi9RbLyvwdM83V9u+YLgp556Sq1bt9agQYN03333yWaz6eOPP9YXX3yhrVu3FmubPXv21I4dO5yWPfjgg2rcuLGeeuqpfKEK8JXVe076ugQA8AvMkgogkLklWDVt2lTz5s3TyJEjHUMFq1Spoo8//lhNmzYt1jajoqLUrFkzp2WRkZGqUqVKvuWAL7365Y++LgEAAAA+5pZgJUmDBw/WwIEDtWPHDhlj1KJFC07nAgAQQK4yUhAASjW3BSvp4rjFNm3auHOTTlatWuWxbQMAAABAcbn8BcEAAAAAgIIRrAAAAADAIoIVAABwCZP+AUDhCFYAAAAAYBHBCgAAuIRJ/wCgcAQrAAAAALCIYAUAAFzCNVYAUDiCFQAAAABYRLACAAAu4RorACgcwQoAAAAALCJYAQAAt7DZOKcFIHARrAAAgEuuNnmFMUxvASBwEawAAAAAwCKCFeBBfHoLoDS52kA/hgICCGQEKwAA4BanMy74ugQA8BmCFeBBnLACUJpc7S1twebDXqkDAPwRwQoAALhNbi6fKAEITAQrAADgEleuoMomWAEIUAQrwIM4vAAQaLJzc31dAgD4BMEKAAC4xJUPi+w5fKQEIDARrAAPYrp1AIEmO4czVgACE8EKAAC4xJVrrHL4QAlAgCJYAR7E4QWAgMMbH4AARbACAABuQ64CEKgIVoAHMSIGQGnCWxoAFI5gBQAA3IYPlAAEKoIVAABwiSuTVxjOawEIUAQrwIM4wAAAAAgMBCsAAOASVz4qYigggEBFsAI8iAMMAIFm8rIffF0CAPgEwQoAALjElWuslm4/6vE6AMAfEawAAAAAwCKCFQAAAABYRLACgBIsIytbh1IyfV0GAgSXjQJA4QhWgAcxeQU87cbJ36jLayv18/Ezvi4FAICARrACgBLqSOo5ncnKliSt/DHZx9UgELgyeQUABCqCFeBBfEEwPOnD9fsd/5/D6VEAAHyKYAUApQC5Ct5AmwFA4QhWgAdxsAtvycml2QAA8CWCFQCUAgQreAPXWAFA4QhWgAdxqAtvMZweBQDApwhWAFAK/OPbvfpq51Ffl4FSjvgOAIUjWAFAKTHyoy2+LgEAgIBFsAI8iOFZAEoTrrECgMIRrAAAwFVt3J+i/acyfV0GAPgtghXgQZyvAlBaDJm53tclAIBfI1gBAAAAgEUEK8CDuMQKAAAgMBCsAAAAAMAighUAAAAAWESwAjyJoYAAAAABgWAFAAAAABYRrAAPMpyyAgAACAgEKwAAAACwiGAFeBDTrQMAAAQGghUAACjU0bTzOnchx9dlAIDf89tgFRcXpxtuuEFRUVGqXr26br/9dv3000++LgsoEk5YASjJTp6Xur7xnbq89q2vSwEAv+e3wSohIUGjR49WYmKi4uPjlZ2drd69eysjI8PXpQEAEBB2n7ZJkk6eveDjSgDA/4X4uoDCfPXVV063Z82aperVq2vz5s3q2rVrgY/JyspSVlaW43Z6erokyW63y263e65YoBAXe8/m6zJQiuS9l9ntduXm5BZ6P2CV3W5X/g5z/bGAqy59XwM8ydM95rfB6nJpaWmSpMqVKxe6TlxcnCZNmpRv+cqVK1W2bFmP1YZAV/iv0YoVK1Qu1IulIGDEx8fr1wNBunzgwbJly3xTEEqlXFO8D4boQxRHfHy8r0tAKZeZmenR7duM8f95y4wxGjRokE6fPq3Vq1cXul5BZ6zq1Kmjo0ePqkqVKt4oFQGo0fPLC70v8enuqhJZxovVoLSz2+2Kj49Xr1699Pdv9+m9Nfud7v9xUi8FB3GWFNbZ7XY9NXuFlhwMLvJj97zU2wMVobS69H0tNJRPI+E5p06dUnR0tNLS0lS+fHm3b79EnLF69NFHtX37dq1Zs+aK64WFhSksLCzf8tDQUH5R4RMhISH0HjwiNDRUQcH5L5Pt+sZ3+vaJ7ioXViLe3uHncov50SvveygOjtfgaZ7uL7+dvCLPY489psWLF2vlypWqXbu2r8sBAL+WfCZLy7Yf9XUZKCX8fkgLAPgRv/1I0xijxx57TIsWLdKqVasUGxvr65KAIvP/gbYojQyHw3CT4l5jBQCByG+D1ejRozVnzhx9/vnnioqK0rFjxyRJFSpUUEREhI+rAwD/RaCHuxR3KCAABCK/HQo4Y8YMpaWlqXv37oqOjnb8zJ8/39elAYBf41gY7lLc6dYBIBD57RmrEjBZIXBVDMmCpxhj9E7Cr74uA6Ucf4oBwHV+e8YKAFC41XtPFXofB8NwF4YCAoDrCFaAJ3FQAg85nVn4t8dzphRWGWP0+vKftT6ZySsAwFV+OxQQAOAsOydXX+w4ppPnpPjNhwtdjzNWsGrTgdN6d/V+SQQrAHAVwQrwII5v4U7/STygSUt26+Jb92lfl4NS7HTGBV+XAAAlDkMBAaCESPj5hEvrEehhFT0EAEVHsAI8iCFZcCeX+4nGg0W0EAAUHcEKsKhWRb6wGt7hcq7yaBUIDHQRABQVwQoASgi+3w/eQqsBQNERrAAPYtpruJOrB7scFMMqvr8KAIqOYAUAJYSrQZ0zW7CKD4UAoOgIVoBFVzqI5fgW7uTyGSvPloEAwHsXABQdwQoASggOduEtuTQbABQZwQqw6EqHHxyawJ1cHwro4UIAAEA+BCsAKCEYCghvIZwDQNERrACLOACBt7j+/cA0JayxOhTw/g826MCpDDdVAwAlA8EKsOhKw7M4wIVb0U7wEqtvXav3nNSYeUluqQUASgqCFQCUEEyBjZLkWNo5X5cAAF5FsAIsutInu5ywgjvxBcHwFne0EH0IINAQrACghHD5GivObMEid0y3ThcCCDQEK8CiKx088F0wcCdXr9mj7WCZG3qIPgQQaAhWgAdxYAF3crWdcuk7WOSes540IoDAQrACLLriNVbeKwMBwNWgzplSWOWOcE7ABxBoCFaAB3GAC3dytZtyOKKFRe546+LrJgAEGoIVYNmVvsfKi2Wg9HOxoQj0sModQwHpQgCBhmAFWHTlY1gOLeA+Ll9jxRkrWHQoxfp3UJHvAQQaghXgQbPX7fd1CShFXD1QzeGIFhbNTPjF8jYYCggg0BCsAIuudOjwUeJBr9WB0s/V4Vmr95xUauYFD1cDXBmxCkCgIVgBQAnh6gmA7YfTNHD6Gs8WA1wFJ6wABBqCFWARw13gLUVpNXdcIwNYwXsjgEBDsAKAEoLDVJQk9CuAQEOwAizi4AHewhkAlCRM+w8g0BCsAIs4doC30GsoSehXAIGGYAUAJYQ7vrQV8Ba6FUCgIVgBFjE8C95Cq8Eb3PaeRr8CCDAEKwAoIThOhTfk5Lqn0zjDCiDQEKwAizh0gLdwdhTekO2uYEW7AggwBCsAKCE4ToU3uO+MFQAEFoIVYBVHD/AWeg1e4L4zVjQsgMBCsAKAEoLDVHhDdk6uW7bjpnyGAJR85rzeWrlXJ85k+boUoEhCfF0AUNJx7ABv4QwAvMFdQwGB4vrD7E3acSRN3/6YrE8fucnX5QAu44wVYBEHu/AWOg3e4K6hgEBx7TiSJknafOC0jysBioZgBQAlBBke3pCdQ6MBQHEQrACLOASBt/C9QPCG7Fz3XGMFuGpq/B69/2OQcjlbihKOYAUAJQRnrOAN7rzGasm235R2zu627aF0mvHdPu04HaSNDP1DCUewAiziYBfeQq/BG9x5jdVjc7fqkY82u217KH0uvU45t4A3Oa5jRklCsAKAEqKggw7A3dx9jdW6X065dXsoXS4N8kE2W777T5y9OOX6iTNZhCz4PaZbByziuhd4i91N3y8EXAnXWMGbfj2R4fj/kCCbDqVkOt3/3c8ntXrPCX2e9Jvuu7GuXr69ubdLBFzGGSvAzQa3ruXrElBKZWVzwAvP43us4E23TvvO8f/BQTZ1eW2l0/1PLNimz5N+kyR9lHiwwG1kZGVzLR/8AsEKsOjSkQnVo8I09a5WPqsFpVtRz1gxbAbFYWe6dfiIK6E+6VCqvth+1HHbGKPWL8ar5aTlOnchx5PlAVdFsAKKKSfX6Pcz1jmdRcg/Olz6ZOMh7xWFUu1CEc9Y/d/mw/rxWLqyGUKIq1i796TW//daKM5YwVdcmTjl9rfWavScLdp+OFWSdCEnVxf++x538LJhhIC3EayAYvrlxFmXvhV+/KfbvVANAkFRj3ef/L/t6jNttR6ds9UzBaFUOJuVrXvf36B73ktUWqZd73z3i9v3QViDK4rSJz8fPytJOn/hfx8cMcEPfI1gBRRTUEGnpwpaBvjYV7uO+boE+KFDKZnadzJDmVnZjmUtX1yu1XtOun1fDZ5d5vZtouS7fHhzUYLVi0t2KSs7Rxv3pziWFfWsPuBuzAoIFFNIUP7PJRrXjPJBJQgEXC8Fd8rNNbp12nfKvJCjz0d38nU5CCBns7KVdDBVN15TWefsztdEHT+T5fJ20s9n67q/fuW0bP6mQ2pZp6I7ygSKhTNWQDGcu5Cj5z/f6bSsQ7Vcvfq7pgWun5x+3htloRR7dtHOq68EuOicPUeZ/73Qf/2vRfueqVHdGxRrn/tOZlx9JZR6w/71ve77YIP+tXZfvskmnlm0y9K252woeNZAwFsIVkAxvL/613zDZe5ukKuq5cIKXH9v8llvlIVS5uCpTD3/2U4dPJWpud9bO2BYt9f9w7tQcn1/yfCpjEuGArpifJ/GGtOzkdOy8uEhevLW6674uDeW/1Sk/aB0yrs2+ZNNh5nFD6WO3wert99+W7GxsQoPD1fbtm21evVqX5eEADdnw0H9Lf7nfMsLvObqv4a+v0EzE9x/QThKt6cXbtd/Eg/o/n9tsLytoe9v0KNztuj/Nh92Q2UoiX44mq6TZy8OtXp75V7H8n9+u7ewhxRqXK9rNaJLrOP2tgm9NbR93Ss+xv7f61+2HUpVm5fiXZr8ByVT5oVsx/VOJ89m6f4PNuT7cKdu5bI6cdb1oX+uqv/0F1q954Tbtwu4wq+D1fz58zV27Fg999xz2rp1q7p06aK+ffvq4EFO9cJ3nl2046rrtKxdId+yV7/8UTuPpCk184InykIpsWjrYe08kiZJWvff6a8PnHLPFMJLtx/VEwu26cdj6W7ZHkqOQymZ6vvmarV7eYVyc41OZRTvfSi6Qrjj/6PCQx3/b7PZFBZ65UOK5buPq/7TX2jQW2uVknFBv5+xTt/+eFwZWdnam3xWucwcWCrsP5mh61/4Wve+nyh7Tq76vrlaq/ec1ND3N+hvl5y13HLwtIbMXO+RGu7/4Hvd9/4GbT1IeIczT1+vbDN+fEV0hw4d1KZNG82YMcOxrEmTJrr99tsVFxd31cenp6erQoUK+tN7qxRerrwuf6aXP/HLX4p8L0y+x1+2/uX3myKsm6+WK++8KM/latu+/P7cXCMjo4jQEIUE2RQSbFOQzSZbEWa8s9pVVpvSalsn/pqisJAg1a1cVunn7dr129UPRB+6Nkdt2rRRSEiw9iaf1RvL85/VckW7epUUFX5xXpmyYf+bXybIZrviWbGrPeUr3X211+uqr+YVVri89/Pv+yqb9uS2LTz2aq/KlR5vs108EA0JsumrXcdkzMVPb9PO2ZV2zn61HeuWWrlacST/QWyT6PL64WjRQ1OL2hVks9m07VCqQoNtsucYlQkOUo/G1Yr8u5/HVtwpMou1r2Lu6r9PzOa4XcwNFaCgf/+r/c25tC7bVerJNdLxtPPa9VuaOjaoqvDQIJ235+pURpYOpZyTJF1TLVK/JJ8tdoiqHhWm5EsmE3i46zVqHB2liNBgSdL0lXu188jFfpt5XxtJ0siPthRrX64KshX9awcud0uTGkrNvKDNB08rNChInRtVVa652PNlQoIc/042XZy+++J7r80xlXfev8+VuFpi3j7ytmn0v76wXeU9P+91uLx/r/be46jR5K817+5LX4NLt5trjIyRsnNzlZNrFBxk03l7rs5dyFGmPVvnLuTolxP+eS1ds1rldeZ8tqpEllGOkcqFBSssJFiRYSFO/54228XZCY0uPv+8f/tLfy9L4wTAhbWNK+9lF9fLv7TAbRayo4L+nhe47wLryb/wx2NnHB9KhoUEOb5zNDcrU4em3am0tDSVL1++4GIs8NtgdeHCBZUtW1YLFizQ7373O8fyxx9/XElJSUpISMj3mKysLGVl/e+PQFpamurWrataj8xWUFhZr9QNAAAAwP/kZmXqyIzhSk1NVYUK+UcXWeW3062fPHlSOTk5qlGjhtPyGjVq6Nixgr+TJS4uTpMmTcq3/MiM4Z4oEQAAAEAJc+rUqcAKVnlsl42FMP89HVuQZ555RuPGjXPcTk1NVb169XTw4EGPvHhAnvT0dNWpU0eHDh3yyKllIA+9Bm+h1+At9Bq8JW80W+XKlT2yfb8NVlWrVlVwcHC+s1PJycn5zmLlCQsLU1hY/umuK1SowC8qvKJ8+fL0GryCXoO30GvwFnoN3hIU5Jn5+/x2VsAyZcqobdu2io+Pd1oeHx+vm266yUdVAQAAAEB+fnvGSpLGjRun+++/X+3atVPHjh317rvv6uDBgxo5cqSvSwMAAAAAB78OVnfddZdOnTqlF198UUePHlWzZs20bNky1atXz6XHh4WFacKECQUODwTciV6Dt9Br8BZ6Dd5Cr8FbPN1rfjvdOgAAAACUFH57jRUAAAAAlBQEKwAAAACwiGAFAAAAABYRrAAAAADAolIbrN5++23FxsYqPDxcbdu21erVq31dEkqY7777TgMHDlRMTIxsNps+++wzp/uNMZo4caJiYmIUERGh7t27a9euXU7rZGVl6bHHHlPVqlUVGRmp2267TYcPH/bis4C/i4uL0w033KCoqChVr15dt99+u3766Sendeg1uMuMGTPUokULxxexduzYUV9++aXjfnoNnhAXFyebzaaxY8c6ltFrcJeJEyfKZrM5/dSsWdNxvzd7rVQGq/nz52vs2LF67rnntHXrVnXp0kV9+/bVwYMHfV0aSpCMjAy1bNlS06dPL/D+1157TVOnTtX06dO1ceNG1axZU7169dKZM2cc64wdO1aLFi3SvHnztGbNGp09e1YDBgxQTk6Ot54G/FxCQoJGjx6txMRExcfHKzs7W71791ZGRoZjHXoN7lK7dm29+uqr2rRpkzZt2qSbb75ZgwYNchxk0Gtwt40bN+rdd99VixYtnJbTa3Cnpk2b6ujRo46fHTt2OO7zaq+ZUqh9+/Zm5MiRTssaN25snn76aR9VhJJOklm0aJHjdm5urqlZs6Z59dVXHcvOnz9vKlSoYGbOnGmMMSY1NdWEhoaaefPmOdY5cuSICQoKMl999ZXXakfJkpycbCSZhIQEYwy9Bs+rVKmSef/99+k1uN2ZM2dMo0aNTHx8vOnWrZt5/PHHjTG8r8G9JkyYYFq2bFngfd7utVJ3xurChQvavHmzevfu7bS8d+/eWrdunY+qQmmzb98+HTt2zKnPwsLC1K1bN0efbd68WXa73WmdmJgYNWvWjF5EodLS0iRJlStXlkSvwXNycnI0b948ZWRkqGPHjvQa3G706NHq37+/brnlFqfl9Brcbc+ePYqJiVFsbKzuvvtu/frrr5K832shbngufuXkyZPKyclRjRo1nJbXqFFDx44d81FVKG3yeqmgPjtw4IBjnTJlyqhSpUr51qEXURBjjMaNG6fOnTurWbNmkug1uN+OHTvUsWNHnT9/XuXKldOiRYt0/fXXOw4g6DW4w7x587RlyxZt3Lgx3328r8GdOnTooA8//FDXXnutjh8/rpdfflk33XSTdu3a5fVeK3XBKo/NZnO6bYzJtwywqjh9Ri+iMI8++qi2b9+uNWvW5LuPXoO7XHfddUpKSlJqaqo+/fRTDRs2TAkJCY776TVYdejQIT3++ONavny5wsPDC12PXoM79O3b1/H/zZs3V8eOHdWgQQP9+9//1o033ijJe71W6oYCVq1aVcHBwfkSZnJycr60ChRX3mwzV+qzmjVr6sKFCzp9+nSh6wB5HnvsMS1evFgrV65U7dq1HcvpNbhbmTJl1LBhQ7Vr105xcXFq2bKl3nzzTXoNbrN582YlJyerbdu2CgkJUUhIiBISEvSPf/xDISEhjl6h1+AJkZGRat68ufbs2eP197VSF6zKlCmjtm3bKj4+3ml5fHy8brrpJh9VhdImNjZWNWvWdOqzCxcuKCEhwdFnbdu2VWhoqNM6R48e1c6dO+lFOBhj9Oijj2rhwoX69ttvFRsb63Q/vQZPM8YoKyuLXoPb9OzZUzt27FBSUpLjp127drr33nuVlJSka665hl6Dx2RlZemHH35QdHS099/XijTVRQkxb948Exoaaj744AOze/duM3bsWBMZGWn279/v69JQgpw5c8Zs3brVbN261UgyU6dONVu3bjUHDhwwxhjz6quvmgoVKpiFCxeaHTt2mHvuucdER0eb9PR0xzZGjhxpateubVasWGG2bNlibr75ZtOyZUuTnZ3tq6cFP/PII4+YChUqmFWrVpmjR486fjIzMx3r0Gtwl2eeecZ89913Zt++fWb79u3m2WefNUFBQWb58uXGGHoNnnPprIDG0Gtwn7/85S9m1apV5tdffzWJiYlmwIABJioqynHc781eK5XByhhj3nrrLVOvXj1TpkwZ06ZNG8fUxYCrVq5caSTl+xk2bJgx5uIUnhMmTDA1a9Y0YWFhpmvXrmbHjh1O2zh37px59NFHTeXKlU1ERIQZMGCAOXjwoA+eDfxVQT0mycyaNcuxDr0Gd3nooYccfxurVatmevbs6QhVxtBr8JzLgxW9Bne56667THR0tAkNDTUxMTFm8ODBZteuXY77vdlrNmOMKfa5NgAAAABA6bvGCgAAAAC8jWAFAAAAABYRrAAAAADAIoIVAAAAAFhEsAIAAAAAiwhWAAAAAGARwQoAAAAALCJYAQAAAIBFBCsAAAAAsIhgBQAAAAAWEawAACXWW2+9pfr16yskJERPPvmkr8sBAAQwghUAoETauXOnxo4dq7feekuHDh3SpEmTNHz4cD399NOOdbp27ao//OEP+R779ttvq2zZssrJyfFmyQCAUoxgBQAokRYvXqy2bduqf//+io6OVnh4uL744gsNGjRIkmSMUVJSktq2bZvvsZs3b1bLli0VHBzs7bIBAKUUwQoAUOI0aNBAzz33nDZs2CCbzab7779fa9euVVBQkDp06CBJ2rNnj86cOVNosMpbPnnyZNlstnw/U6dO9epzAgCUbAQrAECJs379el1zzTV6/fXXdfToUb399ttavHixBg4cqKCgi3/aNm/erODgYLVo0cLpsVlZWdq1a5cjWD322GM6evSo4+eRRx5RvXr1dOedd3r9eQEASq4QXxcAAEBRlStXTvv371fnzp1Vs2ZNSReHBr7xxhuOdbZs2aKcnByVLVu2wG3kBauoqChFRUVJkiZNmqRly5YpISFBtWvX9vCzAACUJjZjjPF1EQAAFEViYqI6deqk9PR0RUZG6ocfflC7du108uRJRURESJJuvvlmVaxYUS+88ILTYxcsWKCpU6fqzJkzCgn53+eLkyZN0qxZs5SQkKB69ep59fkAAEo+hgICAEqcpKQkNWzYUJGRkZIunq3q1auXI1RJ0tatW9W9e3e1atXK6SclJUUtWrQgVAEA3IpgBQAocZKSktSyZUvH7c8//1y33Xab4/avv/6q1NRUtWnTJt9jt2zZ4jShBaEKAOAOBCsAQImTlJSkVq1aSZKSk5O1ceNGDRgwwHH/5s2bFRQU5FgnT3Z2trZv3+4IVi+//LKmT5+u+fPnKywsTMeOHdOxY8eUlZXlracCACglCFYAgBIlNzdXO3bscJyxWrJkiTp06KDq1as71tmyZYsaNWqkcuXKOT12165dOn/+vNq0aSNjjF5//XWdPHlSN954o6Kjox0/SUlJ3nxKAIBSgMkrAAAl2m233abOnTtr/Pjxvi4FABDAOGMFACjROnfurHvuucfXZQAAAhxnrAAAAADAIs5YAQAAAIBFBCsAAAAAsIhgBQAAAAAWEawAAAAAwCKCFQAAAABYRLACAAAAAIsIVgAAAABgEcEKAAAAACwiWAEAAACARf8PLbjL2hRGEvwAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "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/vocal_o_8k.wav\")\n", "x = np.asarray(x, dtype=float)\n", "N = len(x)\n", "\n", "# compute ACF\n", "acf = 1 / N * 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, 4))\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 / Hz$\")\n", "plt.axis([0, 500, 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 average spectral contents of a speech signal?\n", "\n", "Solution: The speech signal exhibits a harmonic structure with the dominant fundamental frequency $f_0 \\approx 100$ Hz and a number of harmonics $f_n \\approx n \\cdot f_0$ for $n > 0$. This due to the fact that [vowels](https://en.wikipedia.org/wiki/Vowel) generate random signals which are in good approximation periodic. To generate vowels, the sound produced by the periodically vibrating [vocal folds](https://en.wikipedia.org/wiki/Vocal_cords) is filtered by the resonance volumes and articulators above the voice box. The spectrum of periodic signals is a line spectrum." ] }, { "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 defined 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 again, that the DTFT is performed with respect to $\\kappa$. The CCF of two wide-sense ergodic random signals of finite length $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[\\kappa] * y_M[-\\kappa].\n", "\\end{equation}\n", "\n", "Note, the $\\frac{1}{N}$ normalization corresponds to the length of the signal $x$. If $N \\neq M$, care should be taken on the interpretation of this normalization. In case of $N=M$ the $\\frac{1}{N}$-averaging yields a [biased estimator](https://en.wikipedia.org/wiki/Bias_of_an_estimator) of the CCF, which consistently should be denoted with $\\hat{\\varphi}_{xy,\\mathrm{biased}}[\\kappa]$.\n", "\n", "\n", "Taking the DTFT of the left- and right-hand side from above cross-correlation 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 $|\\Phi_{xy}(\\mathrm{e}^{\\,\\mathrm{j}\\, \\Omega})| = |\\Phi_{yx}(\\mathrm{e}^{\\,\\mathrm{j}\\, \\Omega})|$ can be concluded.\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 \\mu_y \\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 = 64$ and $M = 512$." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAA2UAAAGHCAYAAAAwZfrzAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAABbjklEQVR4nO3deXxU1f3/8fckTCYJhIEEskEIEQGRTQWBIMomAZQdl0qNoKiUtQjWutQSqkIr/YF+gyJaBWUR2woqLpEgi0tAFkVlKUVFFkkAIQsQSEI4vz/ojJlsZOUy4fV8PHg8mHvP3Ps593Pvnfnk3jnXZowxAgAAAABYwsfqAAAAAADgckZRBgAAAAAWoigDAAAAAAtRlAEAAACAhSjKAAAAAMBCFGUAAAAAYCGKMgAAAACwEEUZAAAAAFiIogwAAAAALERRBqBcFi5cKJvN5vGvYcOG6tGjh95///0i7W02mxISEi5+oGXw008/yWazaeHChRdtnTt37lRCQoJ++umnIvNGjRqlpk2bXrRYCkpJSVFCQoIyMjIsWf+FvPXWW2rdurUCAgJks9m0bds2q0NCDVTa8VkVEhISZLPZKvx+K88RBRXXjxdffPGinkuBmoaiDECFLFiwQBs2bFBKSopefvll+fr6auDAgVq5cqVHuw0bNuj++++3KMpLz86dOzV9+vRiv/Q9+eSTWrFixcUPSueLsunTp1+SRdnRo0cVHx+vZs2aKSkpSRs2bFCLFi2sDgs1UGnHJ351//33a8OGDR7TKMqAyqlldQAAvFObNm3UsWNH9+t+/fqpfv36evPNNzVw4ED39C5dulgRnldq1qyZ1SFUuezsbAUGBlZqGf/973+Vl5enu+++W927d6+iyDzl5eXJZrOpVi0+FqvC5bI9q2L/9kaNGzdW48aNrQ4DqFG4UgagSvj7+8vPz092u91jeuHbF48ePapx48bp6quvVp06dRQaGqpevXrps88+K7LMefPmqX379qpTp46CgoJ01VVX6fHHH/dok5aWpjFjxqhx48by8/NTTEyMpk+frrNnz3q0O3TokO644w4FBQXJ6XTqzjvvVFpaWpn7V9b1lBbzwoULdfvtt0uSevbs6b790/XX5eJuTbLZbJowYYIWLFigli1bKiAgQB07dtTGjRtljNGsWbMUExOjOnXqqFevXvr+++893p+cnKzBgwercePG8vf315VXXqkxY8bol19+cbdJSEjQH/7wB0lSTEyMO65169ZJks6dO6dnn31WV111lRwOh0JDQ3XPPffo4MGDHuvq0aOH2rRpo08//VRdu3ZVYGCg7rvvvlK363vvvafY2FgFBgYqKChIffr08fgL/KhRo9StWzdJ0p133imbzaYePXqUusyff/5ZDz74oKKiouTn56fIyEjddtttOnz4sCRp3bp1stlsWrRokaZOnapGjRrJ4XC4t91rr72m9u3by9/fX8HBwRo6dKh27drlsY4ff/xRv/nNbxQZGSmHw6GwsDD17t3b47bKNWvWqEePHgoJCVFAQICaNGmi4cOHKzs7u9T4mzZtqgEDBmjFihVq166d/P39dcUVV+j//u//irTdv3+/7r77boWGhsrhcKhVq1b6f//v/+ncuXPuNtdff71uvfVWj/e1bdtWNptNmzdvdk9bvny5bDabvvvuO/e0PXv2aMSIER7Lf+GFFzyWdaHtWZwLHduu26STk5N17733Kjg4WLVr19bAgQP1448/Flne6tWr1bt3b9WtW1eBgYG64YYb9MknnxRp95///Ed33XWXwsLC5HA41KRJE91zzz3Kycm54PFZ2v791ltvKS4uThEREQoICFCrVq306KOP6tSpUyVugwtZuHChWrZs6d7ub7zxRrHtcnNz9fTTT7uPz4YNG+ree+/V0aNHPdq59qukpCRdd911CggI0FVXXaXXXnvNo112drYefvhhxcTEuI+Bjh076s0333S3KXz7YtOmTbVjxw6tX7/evd2aNm2qkydPql69ehozZkyRuH/66Sf5+vpq1qxZFd5GQI1iAKAcFixYYCSZjRs3mry8PJObm2sOHDhgJk2aZHx8fExSUpJHe0lm2rRp7tf/+c9/zNixY82yZcvMunXrzPvvv29Gjx5tfHx8zNq1a93t3nzzTSPJTJw40axatcqsXr3avPTSS2bSpEnuNqmpqSYqKspER0eb+fPnm9WrV5unnnrKOBwOM2rUKHe77Oxs06pVK+N0Ok1iYqL5+OOPzaRJk0yTJk2MJLNgwYJS+1zW9Vwo5iNHjpgZM2YYSeaFF14wGzZsMBs2bDBHjhwxxhgzcuRIEx0dXWT7RUdHm65du5rly5ebFStWmBYtWpjg4GDz0EMPmcGDB5v333/fLFmyxISFhZl27dqZc+fOud8/b948M3PmTPPee++Z9evXm9dff920b9/etGzZ0uTm5hpjjDlw4ICZOHGikWSWL1/ujiszM9MYY8yDDz5oJJkJEyaYpKQk89JLL5mGDRuaqKgoc/ToUfe6unfvboKDg01UVJRJTEw0a9euNevXry9xuy5ZssRIMnFxceadd94xb731lunQoYPx8/Mzn332mTHGmO+//9688MILRpKZMWOG2bBhg9mxY0eJyzx48KCJiIgwDRo0MLNnzzarV682b731lrnvvvvMrl27jDHGrF271kgyjRo1Mrfddpt57733zPvvv2+OHTvmzs9dd91lPvjgA/PGG2+YK664wjidTvPf//7XvZ6WLVuaK6+80ixatMisX7/evP3222bq1KnufXjv3r3G39/f9OnTx7zzzjtm3bp1ZsmSJSY+Pt6kp6eXGL8xxkRHR5tGjRqZJk2amNdee818+OGH5re//a2RZGbNmuVud+TIEdOoUSPTsGFD89JLL5mkpCQzYcIEI8mMHTvW3e7RRx81derUcec7LS3NSDIBAQHmmWeecbcbO3asCQsLc7/esWOHcTqdpm3btuaNN94wq1atMlOnTjU+Pj4mISHB3a607VmcshzbrvNMVFSUue+++8xHH31kXn75ZRMaGmqioqI8tuGiRYuMzWYzQ4YMMcuXLzcrV640AwYMML6+vmb16tXudtu2bTN16tQxTZs2NS+99JL55JNPzOLFi80dd9xhsrKyLnh8lrZ/P/XUU2bOnDnmgw8+MOvWrTMvvfSSiYmJMT179vTo+7Rp00xZvnq5+j948GCzcuVKs3jxYnPllVe6z0Mu+fn5pl+/fqZ27dpm+vTpJjk52fzjH/8wjRo1MldffbXJzs52t42OjjaNGzc2V199tXnjjTfMxx9/bG6//XYjyeM4HTNmjAkMDDSzZ882a9euNe+//77561//ahITE0vsx1dffWWuuOIKc+2117q321dffWWMMeahhx4ytWvXNhkZGR59/MMf/mD8/f3NL7/8csHtAVwOKMoAlIvry0Lhfw6Hw7z44otF2hcuygo7e/asycvLM7179zZDhw51T58wYYKpV69eqbGMGTPG1KlTx+zbt89j+t///ncjyf3lfd68eUaSeffddz3aPfDAA2Uqysq6nrLE/K9//ctI8ihAXUoqysLDw83Jkyfd09555x0jyVxzzTUeBdhzzz1nJJlvv/222HWfO3fO5OXlmX379hXZHrNmzTKSzN69ez3es2vXLiPJjBs3zmP6l19+aSSZxx9/3D2te/fuRpL55JNPSt0Gxpz/MhkZGWnatm1r8vPz3dNPnDhhQkNDTdeuXd3TXF/6//Wvf11wuffdd5+x2+1m586dJbZxLe+mm27ymJ6enm4CAgLMLbfc4jF9//79xuFwmBEjRhhjjPnll1+MJPPcc8+VuI5///vfRpLZtm3bBWMuLDo62thstiLv7dOnj6lbt645deqUMeZ8sSXJfPnllx7txo4da2w2m9m9e7cxxpjVq1cbSebTTz81xhizePFiExQUZMaNG+dRNDRv3tzdR2OM6du3r2ncuLG7OHeZMGGC8ff3N8ePHzfGlLw9S1KW48R1nil4TjDGmC+++MJIMk8//bQxxphTp06Z4OBgM3DgQI92+fn5pn379qZTp07uab169TL16tVzF1nFKe34LOv+7TrO1q9fbySZb775xj2vLEWZ69i47rrrPI7vn376ydjtdo9zhKvAffvttz2WsXnzZiPJ45wcHR1t/P39Pc5jp0+fNsHBwWbMmDHuaW3atDFDhgwpNcbi+tG6dWvTvXv3Im1/+OEH4+PjY+bMmeOx3pCQEHPvvfeWuh7gcsLtiwAq5I033tDmzZu1efNmffTRRxo5cqTGjx+vuXPnXvC9L730kq677jr5+/urVq1astvt+uSTTzxuEevUqZMyMjJ011136d133/W43c7l/fffV8+ePRUZGamzZ8+6//Xv31+StH79eknS2rVrFRQUpEGDBnm8f8SIEWXqa1nXU5aYK6Jnz56qXbu2+3WrVq0kSf379/e4hcg1fd++fe5pR44c0e9+9ztFRUW5t3V0dLQkFbklrzhr166VdP42woI6deqkVq1aFblFrH79+urVq9cFl7t7924dOnRI8fHx8vH59aOoTp06Gj58uDZu3HjB2/yK89FHH6lnz57ubVGa4cOHe7zesGGDTp8+XaSvUVFR6tWrl7uvwcHBatasmWbNmqXZs2fr66+/9rhdUJKuueYa+fn56cEHH9Trr79e7C13pWndurXat2/vMW3EiBHKysrSV199Jen87ZFXX321OnXq5NFu1KhRMsZozZo1kqQbbrhB/v7+Wr16taTzt7T26NFD/fr1U0pKirKzs3XgwAHt2bNHN998syTpzJkz+uSTTzR06FAFBgZ67Pe33HKLzpw5o40bN5a6PUtSnuPkt7/9rcfrrl27Kjo62r1fpqSk6Pjx4xo5cqRHjOfOnVO/fv20efNmnTp1StnZ2Vq/fr3uuOMONWzYsExxFqek/fvHH3/UiBEjFB4eLl9fX9ntdvfvH8tynBXkOjZGjBjhcXxHR0era9euHm3ff/991atXTwMHDvTo/zXXXKPw8HD3Lcgu11xzjZo0aeJ+7e/vrxYtWnicMzp16qSPPvpIjz76qNatW6fTp0+XK/7CrrjiCg0YMEAvvviijDGSpKVLl+rYsWOaMGFCpZYN1CQUZQAqpFWrVurYsaM6duyofv36af78+YqLi9MjjzxS6gh+s2fP1tixY9W5c2e9/fbb2rhxozZv3qx+/fp5fPjHx8frtdde0759+zR8+HCFhoaqc+fOSk5Odrc5fPiwVq5cKbvd7vGvdevWkuT+snfs2DGFhYUViSU8PLxMfS3resoSc0UEBwd7vPbz8yt1+pkzZySd/y1YXFycli9frkceeUSffPKJNm3a5P4yXZYvW8eOHZMkRUREFJkXGRnpnu9SXLuKLPfcuXNKT08v07IKOnr0aJkHICi87rL21Waz6ZNPPlHfvn317LPP6rrrrlPDhg01adIknThxQtL5QVtWr16t0NBQjR8/Xs2aNVOzZs30/PPPlym24vZN1zRXHMeOHSsx1oLt/P39dcMNN7iLsk8++UR9+vRRjx49lJ+fr88++8y9j7qKsmPHjuns2bNKTEwsst/fcsstklSkmCpr7stznJS0HVx9c/1O8LbbbisS59/+9jcZY3T8+HGlp6crPz+/0oNTFNfHkydP6sYbb9SXX36pp59+WuvWrdPmzZu1fPlySWU7zgpy9a20fcDl8OHDysjIcP+et+C/tLS0IjkKCQkpskyHw+ER4//93//pj3/8o9555x317NlTwcHBGjJkiPbs2VOufhT0+9//Xnv27HHn+IUXXlBsbKyuu+66Ci8TqGlq9rBIAC6qdu3a6eOPP9Z///vfIn+9d1m8eLF69OihefPmeUx3fZkt6N5779W9996rU6dO6dNPP9W0adM0YMAA/fe//1V0dLQaNGigdu3a6Zlnnil2Xa4vpyEhIdq0aVOR+WUd6KOs6ylLzBfT9u3b9c0332jhwoUaOXKke3ppAzAU5voSl5qaWuQL7aFDh9SgQQOPaWV9BlPB5RZ26NAh+fj4qH79+mWO06Vhw4ZFBiApSeFYLxRTwb5GR0fr1VdflXR+dMh//vOfSkhIUG5url566SVJ0o033qgbb7xR+fn52rJlixITEzV58mSFhYXpN7/5TamxFbdvuqa54gwJCSkxVkke8fbu3Vt//vOftWnTJh08eFB9+vRRUFCQrr/+eiUnJ+vQoUNq0aKFoqKiJJ2/IuTr66v4+HiNHz++2BhjYmI8Xpfn+VtlPU5K2g5XXnmlRx8TExNLHOk1LCxM+fn58vX1LfO+UZLi+rhmzRodOnRI69at8xgdtKKPl3Dlt7R9wKVBgwYKCQlRUlJSscsKCgoq9/pr166t6dOna/r06Tp8+LD7qtnAgQP1n//8p9zLk6RevXqpTZs2mjt3rurUqaOvvvpKixcvrtCygJqKK2UAqoxr5LnSbg+y2WxyOBwe07799tsiz7wpqHbt2urfv7+eeOIJ5ebmaseOHZKkAQMGaPv27WrWrJn7ql3Bf65iqWfPnjpx4oTee+89j+UuXbq0TP0q63rKErOr75W9JagsXF8gC2/v+fPnF2lbUlyuW7UKf4HavHmzdu3apd69e1cotpYtW6pRo0ZaunSp+5YmSTp16pTefvtt94iM5dW/f3+tXbtWu3fvLvd7Y2NjFRAQUKSvBw8e1Jo1a0rsa4sWLfSnP/1Jbdu2dd9aWJCvr686d+7sHrWwuDaF7dixQ998843HtKVLlyooKMh9daF3797auXNnkeW98cYbstls6tmzp3vazTffrLNnz+rJJ59U48aNddVVV7mnr169WmvWrHFfJZOkwMBA9ezZU19//bXatWtX7H5f3FWX8irpOHFZsmSJx+uUlBTt27fPPQLnDTfcoHr16mnnzp3FxtixY0f5+fkpICBA3bt317/+9a9Sb5esyPFZnuOsLFq2bKmIiAi9+eabHsfGvn37lJKS4tF2wIABOnbsmPLz84vte8uWLSsUg0tYWJhGjRqlu+66S7t37y71luLCV9wKmzRpkj744AM99thjCgsLc490CeA8rpQBqJDt27e7h4M/duyYli9fruTkZA0dOrTIX9ALGjBggJ566ilNmzZN3bt31+7du/WXv/xFMTExHsPLP/DAAwoICNANN9ygiIgIpaWlaebMmXI6nbr++uslSX/5y1+UnJysrl27atKkSWrZsqXOnDmjn376SR9++KFeeuklNW7cWPfcc4/mzJmje+65R88884yaN2+uDz/8UB9//HGZ+lrW9ZQl5jZt2kiSXn75ZQUFBcnf318xMTFV8gW3sKuuukrNmjXTo48+KmOMgoODtXLlymJvE2vbtq0k6fnnn9fIkSNlt9vVsmVLtWzZUg8++KASExPl4+Oj/v3766efftKTTz6pqKgoPfTQQxWKzcfHR88++6x++9vfasCAARozZoxycnI0a9YsZWRk6K9//WuFlvuXv/xFH330kW666SY9/vjjatu2rTIyMpSUlKQpU6a4i5Hi1KtXT08++aQef/xx3XPPPbrrrrt07NgxTZ8+Xf7+/po2bZqk839EmDBhgm6//XY1b95cfn5+WrNmjb799ls9+uijks7/bnLNmjW69dZb1aRJE505c8Y99HjB4qckkZGRGjRokBISEhQREaHFixcrOTlZf/vb39zF6kMPPaQ33nhDt956q/7yl78oOjpaH3zwgV588UWNHTvW4wHbHTp0UP369bVq1Srde++97uk333yznnrqqWLjev7559WtWzfdeOONGjt2rJo2baoTJ07o+++/18qVK92/WSuvshwnLlu2bNH999+v22+/XQcOHNATTzyhRo0aady4cZLO/wYxMTFRI0eO1PHjx3XbbbcpNDRUR48e1TfffKOjR4+6r8rPnj1b3bp1U+fOnfXoo4/qyiuv1OHDh/Xee+9p/vz5CgoKqtDx2bVrV9WvX1+/+93vNG3aNNntdi1ZsqRIUV1WPj4+euqpp3T//fdr6NCheuCBB5SRkaGEhIQity/+5je/0ZIlS3TLLbfo97//vTp16iS73a6DBw9q7dq1Gjx4sIYOHVqu9Xfu3FkDBgxQu3btVL9+fe3atUuLFi264B9K2rZtq2XLlumtt97SFVdcIX9/f/d5RZLuvvtuPfbYY/r000/1pz/9yX27NYD/sXKUEQDep7jRF51Op7nmmmvM7NmzzZkzZzzaq9Doizk5Oebhhx82jRo1Mv7+/ua6664z77zzTpGRB19//XXTs2dPExYWZvz8/ExkZKS54447iowsePToUTNp0iQTExNj7Ha7CQ4ONh06dDBPPPGEx4iFBw8eNMOHDzd16tQxQUFBZvjw4SYlJaVMoy+WdT1ljfm5554zMTExxtfX12P9JY2+OH78eI9pe/fuLTI8ujHFj1K4c+dO06dPHxMUFGTq169vbr/9drN///5iR8V87LHHTGRkpPHx8fEYgS4/P9/87W9/My1atDB2u900aNDA3H333ebAgQMe7+/evbtp3br1BbdlQe+8847p3Lmz8ff3N7Vr1za9e/c2X3zxxQX7VZoDBw6Y++67z4SHhxu73e7Ow+HDh8u0vH/84x+mXbt2xs/PzzidTjN48GCPYfgPHz5sRo0aZa666ipTu3ZtU6dOHdOuXTszZ84cc/bsWWOMMRs2bDBDhw410dHRxuFwmJCQENO9e3fz3nvvXTD+6Ohoc+utt5p///vfpnXr1sbPz880bdrUzJ49u0jbffv2mREjRpiQkBBjt9tNy5YtzaxZszxGtHQZOnSokWSWLFninpabm2tq165tfHx8ih2qf+/evea+++4zjRo1Mna73TRs2NB07drVPfphWbZnYWU5TlznmVWrVpn4+HhTr14998iYe/bsKbLM9evXm1tvvdUEBwcbu91uGjVqZG699dYiMe3cudPcfvvtJiQkxPj5+ZkmTZqYUaNGeZy3Sjo+S9u/U1JSTGxsrAkMDDQNGzY0999/v/nqq6+KnF/KOiS+Mef3w+bNmxs/Pz/TokUL89prrxV7jsjLyzN///vfTfv27Y2/v7+pU6eOueqqq8yYMWM8tpVrvyqse/fuHqMmPvroo6Zjx46mfv36xuFwmCuuuMI89NBDHkPXF9ePn376ycTFxZmgoCD3ozwKGzVqlKlVq5Y5ePBgmbYBcDmxGVPg2jgAALBU06ZN1aZNG73//vtWh2KZhQsX6t5779XmzZvVsWNHq8NBFcjNzVXTpk3VrVs3/fOf/7Q6HOCSw+2LAAAAqBZHjx7V7t27tWDBAh0+fNh9iy8ATxRlAAAAqBYffPCB7r33XkVEROjFF19kGHygBNy+CAAAAAAWYkh8AAAAALAQRRkAAAAAWIiiDAAAAAAsxEAfVejcuXM6dOiQgoKCZLPZrA4HAAAAgEWMMTpx4oQiIyPl41P6tTCKsip06NAhRUVFWR0GAAAAgEvEgQMH1Lhx41LbUJRVoaCgIEnnN3zdunUtjqZy8vLytGrVKsXFxclut1sdDiqAHHo38ufdyJ/3I4fejfx5t5qSv6ysLEVFRblrhNJQlFUh1y2LdevWrRFFWWBgoOrWrevVB8PljBx6N/Ln3cif9yOH3o38ebealr+y/KyJgT4AAAAAwEIUZQAAAABgIYoyAAAAALAQRRkAAAAAWIiiDAAAAAAsRFEGAAAAABaiKAMAoArlnzP6cu9xbf3Fpi/3Hlf+OWN1SACASxzPKQMAoIokbU/V9JU7lZp5RpKv3tizRRFOf00beLX6tYmwOjwAwCWKK2UAAFSBpO2pGrv4q/8VZL9KyzyjsYu/UtL2VIsiAwBc6ijKAACopPxzRtNX7lRxNyq6pk1fuZNbGQEAxaIoAwCgkjbtPV7kCllBRlJq5hlt2nv84gUFAPAaFGUAAFTSkRMlF2QVaQcAuLxQlAEAUEmhQf5V2g4AcHmhKAMAoJI6xQQrwukvWwnzbZIinP7qFBN8McMCAHgJijIAACrJ18emaQOvlqQihZnr9bSBV8vXp6SyDQBwOaMoAwCgCvRrE6F5d1+n0LoOj+nhTn/Nu/s6nlMGACgRD48GAKCK9GsToRuubKC2CaskSf+Iv1Y9W0VwhQwAUCqulAEAUIUKFmDXN61PQQYAuCCKMgAAAACwEEUZAAAAAFiIogwAAAAALERRBgAAAAAWuuSLsnnz5qldu3aqW7eu6tatq9jYWH300Ufu+cYYJSQkKDIyUgEBAerRo4d27NjhsYycnBxNnDhRDRo0UO3atTVo0CAdPHjQo016erri4+PldDrldDoVHx+vjIyMi9FFAAAAAJexS74oa9y4sf76179qy5Yt2rJli3r16qXBgwe7C69nn31Ws2fP1ty5c7V582aFh4erT58+OnHihHsZkydP1ooVK7Rs2TJ9/vnnOnnypAYMGKD8/Hx3mxEjRmjbtm1KSkpSUlKStm3bpvj4+IveXwAAAACXl0v+OWUDBw70eP3MM89o3rx52rhxo66++mo999xzeuKJJzRs2DBJ0uuvv66wsDAtXbpUY8aMUWZmpl599VUtWrRIN998syRp8eLFioqK0urVq9W3b1/t2rVLSUlJ2rhxozp37ixJeuWVVxQbG6vdu3erZcuWF7fTAAAAAC4bl3xRVlB+fr7+9a9/6dSpU4qNjdXevXuVlpamuLg4dxuHw6Hu3bsrJSVFY8aM0datW5WXl+fRJjIyUm3atFFKSor69u2rDRs2yOl0ugsySerSpYucTqdSUlJKLMpycnKUk5Pjfp2VlSVJysvLU15eXlV3/6Jyxe/t/bickUPvRv68V17eWY//k0PvxDHo3cifd6sp+StP/F5RlH333XeKjY3VmTNnVKdOHa1YsUJXX321UlJSJElhYWEe7cPCwrRv3z5JUlpamvz8/FS/fv0ibdLS0txtQkNDi6w3NDTU3aY4M2fO1PTp04tMX7VqlQIDA8vXyUtUcnKy1SGgksihdyN/3icnX3J9vK5Zs0YOX0vDQSVxDHo38ufdvD1/2dnZZW7rFUVZy5YttW3bNmVkZOjtt9/WyJEjtX79evd8m83m0d4YU2RaYYXbFNf+Qst57LHHNGXKFPfrrKwsRUVFKS4uTnXr1r1gvy5leXl5Sk5OVp8+fWS3260OBxVADr0b+fNe2bln9cimNZKkXr16yVnb3+KIUBEcg96N/Hm3mpI/1110ZeEVRZmfn5+uvPJKSVLHjh21efNmPf/88/rjH/8o6fyVroiICHf7I0eOuK+ehYeHKzc3V+np6R5Xy44cOaKuXbu62xw+fLjIeo8ePVrkKlxBDodDDoejyHS73e7VO1BBNakvlyty6N3In/exm1//mGe31yJ/Xo5j0LuRP+/m7fkrT+yX/OiLxTHGKCcnRzExMQoPD/e4tJmbm6v169e7C64OHTrIbrd7tElNTdX27dvdbWJjY5WZmalNmza523z55ZfKzMx0twEAAACA6nDJXyl7/PHH1b9/f0VFRenEiRNatmyZ1q1bp6SkJNlsNk2ePFkzZsxQ8+bN1bx5c82YMUOBgYEaMWKEJMnpdGr06NGaOnWqQkJCFBwcrIcfflht27Z1j8bYqlUr9evXTw888IDmz58vSXrwwQc1YMAARl4EAAAAUK0u+aLs8OHDio+PV2pqqpxOp9q1a6ekpCT16dNHkvTII4/o9OnTGjdunNLT09W5c2etWrVKQUFB7mXMmTNHtWrV0h133KHTp0+rd+/eWrhwoXx9f/319ZIlSzRp0iT3KI2DBg3S3LlzL25nAQAAAFx2Lvmi7NVXXy11vs1mU0JCghISEkps4+/vr8TERCUmJpbYJjg4WIsXL65omAAAAABQIV75mzIAAAAAqCkoygAAAADAQhRlAAAAAGAhijIAAAAAsBBFGQAAAABYiKIMAAAAACxEUQYAAAAAFqIoAwAAAAALUZQBAAAAgIUoygAAAADAQhRlAAAAAGAhijIAAAAAsBBFGQAAAABYiKIMAAAAACxEUQYAAAAAFqIoAwAAAAALUZQBAAAAgIUoygAAAADAQhRlAAAAAGAhijIAAAAAsBBFGQAAAABYiKIMAAAAACxEUQYAAAAAFqIoAwAAAAALUZQBAAAAgIUoygAAAADAQhRlAAAAAGAhijIAAAAAsBBFGQAAAABYiKIMAAAAACxEUQYAAAAAFqIoAwAAAAALUZQBAAAAgIUu+aJs5syZuv766xUUFKTQ0FANGTJEu3fv9mgzatQo2Ww2j39dunTxaJOTk6OJEyeqQYMGql27tgYNGqSDBw96tElPT1d8fLycTqecTqfi4+OVkZFR3V0EAAAAcBm75Iuy9evXa/z48dq4caOSk5N19uxZxcXF6dSpUx7t+vXrp9TUVPe/Dz/80GP+5MmTtWLFCi1btkyff/65Tp48qQEDBig/P9/dZsSIEdq2bZuSkpKUlJSkbdu2KT4+/qL0EwAAAMDlqZbVAVxIUlKSx+sFCxYoNDRUW7du1U033eSe7nA4FB4eXuwyMjMz9eqrr2rRokW6+eabJUmLFy9WVFSUVq9erb59+2rXrl1KSkrSxo0b1blzZ0nSK6+8otjYWO3evVstW7asph4CAAAAuJxd8kVZYZmZmZKk4OBgj+nr1q1TaGio6tWrp+7du+uZZ55RaGioJGnr1q3Ky8tTXFycu31kZKTatGmjlJQU9e3bVxs2bJDT6XQXZJLUpUsXOZ1OpaSkFFuU5eTkKCcnx/06KytLkpSXl6e8vLyq67QFXPF7ez8uZ+TQu5E/75WXd9bj/+TQO3EMejfy591qSv7KE79XFWXGGE2ZMkXdunVTmzZt3NP79++v22+/XdHR0dq7d6+efPJJ9erVS1u3bpXD4VBaWpr8/PxUv359j+WFhYUpLS1NkpSWluYu4goKDQ11tyls5syZmj59epHpq1atUmBgYGW6eslITk62OgRUEjn0buTP++TkS66P1zVr1sjha2k4qCSOQe9G/rybt+cvOzu7zG29qiibMGGCvv32W33++ece0++88073/9u0aaOOHTsqOjpaH3zwgYYNG1bi8owxstls7tcF/19Sm4Iee+wxTZkyxf06KytLUVFRiouLU926dcvcr0tRXl6ekpOT1adPH9ntdqvDQQWQQ+9G/rxXdu5ZPbJpjSSpV69ectb2tzgiVATHoHcjf96tpuTPdRddWXhNUTZx4kS99957+vTTT9W4ceNS20ZERCg6Olp79uyRJIWHhys3N1fp6ekeV8uOHDmirl27utscPny4yLKOHj2qsLCwYtfjcDjkcDiKTLfb7V69AxVUk/pyuSKH3o38eR+7+fUPeXZ7LfLn5TgGvRv5827enr/yxH7Jj75ojNGECRO0fPlyrVmzRjExMRd8z7Fjx3TgwAFFRERIkjp06CC73e5xCTQ1NVXbt293F2WxsbHKzMzUpk2b3G2+/PJLZWZmutsAAAAAQFW75K+UjR8/XkuXLtW7776roKAg9++7nE6nAgICdPLkSSUkJGj48OGKiIjQTz/9pMcff1wNGjTQ0KFD3W1Hjx6tqVOnKiQkRMHBwXr44YfVtm1b92iMrVq1Ur9+/fTAAw9o/vz5kqQHH3xQAwYMYORFAAAAANXmki/K5s2bJ0nq0aOHx/QFCxZo1KhR8vX11Xfffac33nhDGRkZioiIUM+ePfXWW28pKCjI3X7OnDmqVauW7rjjDp0+fVq9e/fWwoUL5ev76y+wlyxZokmTJrlHaRw0aJDmzp1b/Z0EAAAAcNm65IsyY0yp8wMCAvTxxx9fcDn+/v5KTExUYmJiiW2Cg4O1ePHicscIAAAAABV1yf+mDAAAAABqMooyAAAAALAQRRkAAAAAWIiiDAAAAAAsRFEGAAAAABaiKAMAAAAAC1GUAQAAAICFKMoAAAAAwEIUZQAAAABgIYoyAAAAALAQRRkAAAAAWIiiDAAAAAAsRFEGAAAAABaiKAMAAAAAC1GUAQAAAICFKMoAAAAAwEIUZQAAAABgIYoyAAAAALAQRRkAAAAAWIiiDAAAAAAsRFEGAAAAABaiKAMAAAAAC1GUAQAAAICFKMoAAAAAwEIUZQAAAABgIYoyAAAAALAQRRkAAAAAWIiiDAAAAAAsRFEGAAAAABaiKAMAAAAAC1GUAQAAAICFKMoAAAAAwEIUZQAAAABgoUu+KJs5c6auv/56BQUFKTQ0VEOGDNHu3bs92hhjlJCQoMjISAUEBKhHjx7asWOHR5ucnBxNnDhRDRo0UO3atTVo0CAdPHjQo016erri4+PldDrldDoVHx+vjIyM6u4iAAAAgMvYJV+UrV+/XuPHj9fGjRuVnJyss2fPKi4uTqdOnXK3efbZZzV79mzNnTtXmzdvVnh4uPr06aMTJ06420yePFkrVqzQsmXL9Pnnn+vkyZMaMGCA8vPz3W1GjBihbdu2KSkpSUlJSdq2bZvi4+Mvan8BAAAAXF5qWR3AhSQlJXm8XrBggUJDQ7V161bddNNNMsboueee0xNPPKFhw4ZJkl5//XWFhYVp6dKlGjNmjDIzM/Xqq69q0aJFuvnmmyVJixcvVlRUlFavXq2+fftq165dSkpK0saNG9W5c2dJ0iuvvKLY2Fjt3r1bLVu2vLgdBwAAAHBZuOSLssIyMzMlScHBwZKkvXv3Ki0tTXFxce42DodD3bt3V0pKisaMGaOtW7cqLy/Po01kZKTatGmjlJQU9e3bVxs2bJDT6XQXZJLUpUsXOZ1OpaSkFFuU5eTkKCcnx/06KytLkpSXl6e8vLyq7fhF5orf2/txOSOH3o38ea+8vLMe/yeH3olj0LuRP+9WU/JXnvi9qigzxmjKlCnq1q2b2rRpI0lKS0uTJIWFhXm0DQsL0759+9xt/Pz8VL9+/SJtXO9PS0tTaGhokXWGhoa62xQ2c+ZMTZ8+vcj0VatWKTAwsJy9uzQlJydbHQIqiRx6N/LnfXLyJdfH65o1a+TwtTQcVBLHoHcjf97N2/OXnZ1d5rZeVZRNmDBB3377rT7//PMi82w2m8drY0yRaYUVblNc+9KW89hjj2nKlCnu11lZWYqKilJcXJzq1q1b6rovdXl5eUpOTlafPn1kt9utDgcVQA69G/nzXtm5Z/XIpjWSpF69eslZ29/iiFARHIPejfx5t5qSP9dddGXhNUXZxIkT9d577+nTTz9V48aN3dPDw8Mlnb/SFRER4Z5+5MgR99Wz8PBw5ebmKj093eNq2ZEjR9S1a1d3m8OHDxdZ79GjR4tchXNxOBxyOBxFptvtdq/egQqqSX25XJFD70b+vI/d/PqHPLu9FvnzchyD3o38eTdvz195Yr/kR180xmjChAlavny51qxZo5iYGI/5MTExCg8P97i8mZubq/Xr17sLrg4dOshut3u0SU1N1fbt291tYmNjlZmZqU2bNrnbfPnll8rMzHS3AQAAAICqdslfKRs/fryWLl2qd999V0FBQe7fdzmdTgUEBMhms2ny5MmaMWOGmjdvrubNm2vGjBkKDAzUiBEj3G1Hjx6tqVOnKiQkRMHBwXr44YfVtm1b92iMrVq1Ur9+/fTAAw9o/vz5kqQHH3xQAwYMYORFAAAAANXmki/K5s2bJ0nq0aOHx/QFCxZo1KhRkqRHHnlEp0+f1rhx45Senq7OnTtr1apVCgoKcrefM2eOatWqpTvuuEOnT59W7969tXDhQvn6/voL7CVLlmjSpEnuURoHDRqkuXPnVm8HAQAAAFzWLvmizBhzwTY2m00JCQlKSEgosY2/v78SExOVmJhYYpvg4GAtXry4ImECAAAAQIVc8r8pAwAAAICarFJXymJiYi447HxxJk+erEmTJlVm1QAAAABQI1SqKFu4cGGF3te0adPKrBYAAAAAaoxKFWXdu3evqjgAAAAA4LLEb8oAAAAAwEJVOvpiXl6e0tLSlJ2drYYNGyo4OLgqFw8AAAAANU6lr5SdPHlS8+fPV48ePeR0OtW0aVNdffXVatiwoaKjo/XAAw9o8+bNVRErAAAAANQ4lSrK5syZo6ZNm+qVV15Rr169tHz5cm3btk27d+/Whg0bNG3aNJ09e1Z9+vRRv379tGfPnqqKGwAAAABqhErdvpiSkqK1a9eqbdu2xc7v1KmT7rvvPr300kt69dVXtX79ejVv3rwyqwQAAACAGqVSRdm//vWvMrVzOBwaN25cZVYFAAAAADUSoy8CAAAAgIUqVJRt2bJFI0aMUIcOHdS5c2c98sgjSktL47llAAAAAFBO5b598c0339ScOXP09NNPq23btjpx4oTef/99de/eXZmZmdURIwAAAADUWOUuyp555hl99tlnql+/viQpIiJCU6ZM0a233qoRI0ZUeYAAAAAAUJOV+/bF/Px8d0FWUMuWLbVkyZIqCQoAAAAALhflLsoaNWqkdevWFZk+e/ZstW7duipiAgAAAIDLRrlvX5w/f75uu+02tWvXTu3atVNWVpY++OADNWvWjGeQAQAAAEA5lftK2b///W9t3bpVLVu21Llz5xQSEqJ//OMfeuutt3TvvfdWR4wAAAAAUGOV+0pZt27d5OPjoy1btmjXrl2qVauWNm7cqLZt26pt27Y6fPiwwsLCqiNWAAAAAKhxyl2U3XDDDZKk5cuXS5Kys7O1fft2fffdd1q9erWmTZumW265RU899VTVRgoAAAAANVC5i7LCAgMD1alTJ3Xq1Mk9rWPHjhRlAAAAAFAG5f5NWUkyMzP14IMP6sorr1RGRoZSU1OratEAAAAAUGNVWVE2btw4fffdd3r22WeVmpqq06dPS5ImT56s559/vqpWAwAAAAA1SpUVZR999JFefPFFDRs2TD4+vy62X79+WrRoUVWtBgAAAABqlCoryiSpTp06RaY1b95c33//fVWuBgAAAABqjCorym655RYtXbq0yPSTJ0/KZrNV1WoAAAAAoEap9OiLLjNnzlTHjh0lScYY2Ww2nT59Wn/5y1903XXXVdVqAAAAAKBGqbKiLCoqSl988YXGjh2r7OxsderUSSdOnFDdunX14YcfVtVqAAAAAKBGqbKiTJKuvPJKJScna//+/frmm29kt9vVuXNn1a9fvypXAwAAAAA1RpUWZS5NmjRRkyZNqmPRAAAAAFCjVGqgj/3795er/c8//1yZ1QEAAABAjVOpouz666/XAw88oE2bNpXYJjMzU6+88oratGmj5cuXV2Z1AAAAAFDjVOr2xV27dmnGjBnq16+f7Ha7OnbsqMjISPn7+ys9PV07d+7Ujh071LFjR82aNUv9+/evqrgBAAAAoEao1JWy4OBg/f3vf9ehQ4c0b948tWjRQr/88ov27NkjSfrtb3+rrVu36osvvqhUQfbpp59q4MCBioyMlM1m0zvvvOMxf9SoUbLZbB7/unTp4tEmJydHEydOVIMGDVS7dm0NGjRIBw8e9GiTnp6u+Ph4OZ1OOZ1OxcfHKyMjo8JxAwAAAMCFVMlAH/7+/ho2bJiGDRtWFYsr4tSpU2rfvr3uvfdeDR8+vNg2/fr104IFC9yv/fz8POZPnjxZK1eu1LJlyxQSEqKpU6dqwIAB2rp1q3x9fSVJI0aM0MGDB5WUlCRJevDBBxUfH6+VK1dWS78AAAAAoFpGX5SkiRMnatq0aWrQoEGll9W/f/8LXmlzOBwKDw8vdl5mZqZeffVVLVq0SDfffLMkafHixYqKitLq1avVt29f7dq1S0lJSdq4caM6d+4sSXrllVcUGxur3bt3q2XLlpXuBwAAAAAUVm1FWY8ePdS3b18NHTpUU6dOVUBAQHWtSpK0bt06hYaGql69eurevbueeeYZhYaGSpK2bt2qvLw8xcXFudtHRkaqTZs2SklJUd++fbVhwwY5nU53QSZJXbp0kdPpVEpKSrFFWU5OjnJyctyvs7KyJEl5eXnKy8urrq5eFK74vb0flzNy6N3In/fKyzvr8X9y6J04Br0b+fNuNSV/5Ym/2oqy4cOHa/DgwZo3b566dOmiiRMnavTo0bLZbFW+rv79++v2229XdHS09u7dqyeffFK9evXS1q1b5XA4lJaWJj8/vyIPsQ4LC1NaWpokKS0tzV3EFRQaGupuU9jMmTM1ffr0ItNXrVqlwMDAKuiZ9ZKTk60OAZVEDr0b+fM+OfmS6+N1zZo1cvhaGg4qiWPQu5E/7+bt+cvOzi5z22oryiSpVq1aGjZsmJxOp6ZMmaI5c+Zo1qxZuuWWW6p0PXfeeaf7/23atFHHjh0VHR2tDz74oNTfuRljPIrE4grGwm0KeuyxxzRlyhT366ysLEVFRSkuLk5169atSFcuGXl5eUpOTlafPn1kt9utDgcVQA69G/nzXtm5Z/XIpjWSpF69eslZ29/iiFARHIPejfx5t5qSP9dddGVRbUVZ//79tWvXLjVu3FidOnVSYmKiWrRooRdeeEHJycmaM2dOda1aERERio6Odo8CGR4ertzcXKWnp3tcLTty5Ii6du3qbnP48OEiyzp69KjCwsKKXY/D4ZDD4Sgy3W63e/UOVFBN6svlihx6N/Lnfezm1z/k2e21yJ+X4xj0buTPu3l7/soTe7UVZTNmzFC7du3cIxu6vPbaa7rqqquqtSg7duyYDhw4oIiICElShw4dZLfblZycrDvuuEOSlJqaqu3bt+vZZ5+VJMXGxiozM1ObNm1Sp06dJElffvmlMjMz3YUbAAAAAFS1aivKNm3apGuvvbbYeR9++GG5lnXy5El9//337td79+7Vtm3bFBwcrODgYCUkJGj48OGKiIjQTz/9pMcff1wNGjTQ0KFDJUlOp1OjR4/W1KlTFRISouDgYD388MNq27atezTGVq1aqV+/fnrggQc0f/58SeeHxB8wYAAjLwIAAACoNpV6eHRpNmzYoIkTJ+rcuXOSpN27dys+Pl6SdMUVV5RrWVu2bNG1117rLvKmTJmia6+9Vn/+85/l6+ur7777ToMHD1aLFi00cuRItWjRQhs2bFBQUJB7GXPmzNGQIUN0xx136IYbblBgYKBWrlzpcSVvyZIlatu2reLi4hQXF6d27dpp0aJFld0UAAAAAFCiartStnDhQs2ePVu33HKL6tWrpz179ugPf/hDhZbVo0cPGWNKnP/xxx9fcBn+/v5KTExUYmJiiW2Cg4O1ePHiCsUIAAAAABVRbUXZV199pS+++EKHDx/Wf//7X61du1bR0dHVtToAAAAA8ErVdvviuHHjNHr0aH399ddatmyZBg8erC+++KK6VgcAAAAAXqnarpRt3LjR/f9OnTrpgw8+0O23366UlJTqWiUAAAAAeJ0qv1KWkZGh9PT0ItMbNWqkNWvWVPXqAAAAAMCrlbsoe/PNN4udvmHDBl177bUKCQlRgwYN1L59+yJXxfz9/SsWJQAAAADUUGUuytLS0jRs2DAlJycXmbdv3z716dNHtWrV0syZM/XXv/5VDodDffr00d69e6s0YAAAAACoScr8m7KXX35ZZ8+e1WuvvVZk3owZM9S1a1d99NFH7ud+TZ06VbfeequeeeYZ/eMf/6i6iAEAAACgBinzlbJJkyapXr16Gj58uMf0/fv3a926dbrzzjv1888/a//+/dq/f78OHjyoO++8U2vXrtWBAweqPHAAAAAAqAnKfKWsXr16euONN7Ry5UqP6U2bNpXNZtODDz5Y7PuMMWratKny8/MrFykAAAAA1EDlHhJ/4MCBHq+/+uorDR06VA8//LC6devmMS8lJUV/+9vf9O6771YuSgAAAACooSr9nLJrrrlGN954o9atW6fx48d7zJs5c6ZuuukmtW/fvrKrAQAAAIAaqUoeHv3HP/5R1157rQYPHqy7775bNptNS5Ys0QcffKCvv/66KlYBAAAAADVSlRRlrVu31rJly/S73/3O/ZuzkJAQLVmyRK1bt66KVQAAAABAjVQlRZkkDRs2TAMHDtR3330nY4zatWsnu91eVYsHAAAAgBqpyooySbLb7bruuuuqcpEAAAAAUKNVqiiLiYmRzWYr9/smT56sSZMmVWbVAAAAAFAjVKooW7hwYYXe17Rp08qsFgAAAABqjEoVZd27d6+qOAAAAADgsuRjdQAAAAAAcDmjKAMAAAAAC1GUAQAAAICFKMoAAAAAwEIUZQAAAABgIYoyAAAAALAQRRkAAAAAWIiiDAAAAAAsRFEGAAAAABaiKAMAAAAAC1GUAQAAAICFKMoAAAAAwEIUZQAAAABgIYoyAAAAALAQRRkAAAAAWMgrirJPP/1UAwcOVGRkpGw2m9555x2P+cYYJSQkKDIyUgEBAerRo4d27Njh0SYnJ0cTJ05UgwYNVLt2bQ0aNEgHDx70aJOenq74+Hg5nU45nU7Fx8crIyOjmnsHAAAA4HLmFUXZqVOn1L59e82dO7fY+c8++6xmz56tuXPnavPmzQoPD1efPn104sQJd5vJkydrxYoVWrZsmT7//HOdPHlSAwYMUH5+vrvNiBEjtG3bNiUlJSkpKUnbtm1TfHx8tfcPAAAAwOWrltUBlEX//v3Vv3//YucZY/Tcc8/piSee0LBhwyRJr7/+usLCwrR06VKNGTNGmZmZevXVV7Vo0SLdfPPNkqTFixcrKipKq1evVt++fbVr1y4lJSVp48aN6ty5syTplVdeUWxsrHbv3q2WLVtenM4CAAAAuKx4RVFWmr179yotLU1xcXHuaQ6HQ927d1dKSorGjBmjrVu3Ki8vz6NNZGSk2rRpo5SUFPXt21cbNmyQ0+l0F2SS1KVLFzmdTqWkpBRblOXk5CgnJ8f9OisrS5KUl5envLy86ujuReOK39v7cTkjh96N/HmvvLyzHv8nh96JY9C7kT/vVlPyV574vb4oS0tLkySFhYV5TA8LC9O+ffvcbfz8/FS/fv0ibVzvT0tLU2hoaJHlh4aGutsUNnPmTE2fPr3I9FWrVikwMLD8nbkEJScnWx0CKokcejfy531y8iXXx+uaNWvk8LU0HFQSx6B3I3/ezdvzl52dXea2Xl+UudhsNo/Xxpgi0wor3Ka49qUt57HHHtOUKVPcr7OyshQVFaW4uDjVrVu3POFfcvLy8pScnKw+ffrIbrdbHQ4qgBx6N/LnvbJzz+qRTWskSb169ZKztr/FEaEiOAa9G/nzbjUlf6676MrC64uy8PBwSeevdEVERLinHzlyxH31LDw8XLm5uUpPT/e4WnbkyBF17drV3ebw4cNFln/06NEiV+FcHA6HHA5Hkel2u92rd6CCalJfLlfk0LuRP+9jN7/+Ic9ur0X+vBzHoHcjf97N2/NXnti9YvTF0sTExCg8PNzj8mZubq7Wr1/vLrg6dOggu93u0SY1NVXbt293t4mNjVVmZqY2bdrkbvPll18qMzPT3QYAAAAAqppXXCk7efKkvv/+e/frvXv3atu2bQoODlaTJk00efJkzZgxQ82bN1fz5s01Y8YMBQYGasSIEZIkp9Op0aNHa+rUqQoJCVFwcLAefvhhtW3b1j0aY6tWrdSvXz898MADmj9/viTpwQcf1IABAxh5EQAAAEC18YqibMuWLerZs6f7tet3XCNHjtTChQv1yCOP6PTp0xo3bpzS09PVuXNnrVq1SkFBQe73zJkzR7Vq1dIdd9yh06dPq3fv3lq4cKF8fX/9BfaSJUs0adIk9yiNgwYNKvHZaAAAAABQFbyiKOvRo4eMMSXOt9lsSkhIUEJCQolt/P39lZiYqMTExBLbBAcHa/HixZUJFQAAAADKxet/UwYAAAAA3oyiDAAAAAAsRFEGAAAAABaiKAMAAAAAC1GUAQAAAICFKMoAAAAAwEIUZQAAAABgIYoyAAAAALAQRRkAAAAAWIiiDAAAAAAsRFEGAAAAABaiKAMAAAAAC1GUAQAAAICFKMoAAAAAwEIUZQAAAABgIYoyAAAAALAQRRkAAAAAWIiiDAAAAAAsRFEGAAAAABaiKAMAAAAAC1GUAQAAAICFKMoAAAAAwEIUZQAAAABgIYoyAAAAALAQRRkAAAAAWIiiDAAAAAAsRFEGAAAAABaiKAMAAAAAC1GUAQAAAICFKMoAAAAAwEIUZQAAAABgIYoyAAAAALAQRRkAAAAAWIiiDAAAAAAsVCOKsoSEBNlsNo9/4eHh7vnGGCUkJCgyMlIBAQHq0aOHduzY4bGMnJwcTZw4UQ0aNFDt2rU1aNAgHTx48GJ3BQAAAMBlpkYUZZLUunVrpaamuv9999137nnPPvusZs+erblz52rz5s0KDw9Xnz59dOLECXebyZMna8WKFVq2bJk+//xznTx5UgMGDFB+fr4V3QEAAABwmahldQBVpVatWh5Xx1yMMXruuef0xBNPaNiwYZKk119/XWFhYVq6dKnGjBmjzMxMvfrqq1q0aJFuvvlmSdLixYsVFRWl1atXq2/fvhe1LwAAAAAuHzWmKNuzZ48iIyPlcDjUuXNnzZgxQ1dccYX27t2rtLQ0xcXFuds6HA51795dKSkpGjNmjLZu3aq8vDyPNpGRkWrTpo1SUlJKLMpycnKUk5Pjfp2VlSVJysvLU15eXjX19OJwxe/t/bickUPvRv68V17eWY//k0PvxDHo3cifd6sp+StP/DWiKOvcubPeeOMNtWjRQocPH9bTTz+trl27aseOHUpLS5MkhYWFebwnLCxM+/btkySlpaXJz89P9evXL9LG9f7izJw5U9OnTy8yfdWqVQoMDKxsty4JycnJVoeASiKH3o38eZ+cfMn18bpmzRo5fC0NB5XEMejdyJ938/b8ZWdnl7ltjSjK+vfv7/5/27ZtFRsbq2bNmun1119Xly5dJEk2m83jPcaYItMKu1Cbxx57TFOmTHG/zsrKUlRUlOLi4lS3bt2KdOWSkZeXp+TkZPXp00d2u93qcFAB5NC7kT/vlZ17Vo9sWiNJ6tWrl5y1/S2OCBXBMejdyJ93qyn5c91FVxY1oigrrHbt2mrbtq327NmjIUOGSDp/NSwiIsLd5siRI+6rZ+Hh4crNzVV6errH1bIjR46oa9euJa7H4XDI4XAUmW632716ByqoJvXlckUOvRv58z528+sf8+z2WuTPy3EMejfy5928PX/lib3GjL5YUE5Ojnbt2qWIiAjFxMQoPDzc4/Jnbm6u1q9f7y64OnToILvd7tEmNTVV27dvL7UoAwAAAIDKqhFXyh5++GENHDhQTZo00ZEjR/T0008rKytLI0eOlM1m0+TJkzVjxgw1b95czZs314wZMxQYGKgRI0ZIkpxOp0aPHq2pU6cqJCREwcHBevjhh9W2bVv3aIwAAAAAUB1qRFF28OBB3XXXXfrll1/UsGFDdenSRRs3blR0dLQk6ZFHHtHp06c1btw4paenq3Pnzlq1apWCgoLcy5gzZ45q1aqlO+64Q6dPn1bv3r21cOFC+fryC20AAAAA1adGFGXLli0rdb7NZlNCQoISEhJKbOPv76/ExEQlJiZWcXQAAAAAULIa+ZsyAAAAAPAWNeJKGQAAVso/Z7Rp73GlZZ5WWtYZ9/TENT/I38+u2GYh6nJFiHx9Sn8UCwDg8kRRBgBAJSRtT9X0lTuVmnmmyLx/fLFPkjR37feqF2jXX4e1Vb82EUXaAQAub9y+CABABSVtT9XYxV8VW5AVlpGdp98t/kpJ21MvQmQAAG9CUQYAQAXknzOavnKnTDnfN33lTuWfK++7AAA1GUUZAAAVsGnv8TJdISssNfOMNu09Xg0RAQC8FUUZAAAVcORE+QuyqngvAKDmoSgDAKACQoP8LXkvAKDmoSgDAKACOsUEK8JZ/uIqwumvTjHB1RARAMBbUZQBAFABvj42TRt4tcr75LFpA6/meWUAAA8UZQAAVFC/NhF6YcR1qh9ov2DbOg5fvTjiWp5TBgAogqIMAIAKStqeqqc+2Kn07LwLtj2Zk6+nPtjFc8oAAEVQlAEAUAHleXC0S1rmGY3lAdIAgEIoygAAKKeKPjja1Z4HSAMACqIoAwCgnCr64GjpfGHGA6QBAAVRlAEAUE5V8fBnHiANAHChKAMAoJyq4uHPPEAaAOBCUQYAQDl1iglWvTIMg18SH5uUfiq3CiMCAHgzijIAAMopeWeaMsowDH5Jzhlp/FJGYQQAnEdRBgBAObhGXqwsI0ZhBACcR1EGAEA5VGbkxcIYhREAIFGUAQBQLqt3plXp8pKreHkAAO9DUQYAQBnlnzNase3nKl3mu9sOcQsjAFzmKMoAACijTXuP6/ipig/wUZxjp3K5hREALnMUZQAAlFF1PfCZB0kDwOWNogwAgDKqrgc+8yBpALi8UZQBAFBG6ady5GOr2mXWC7SrU0xw1S4UAOBVKMoAACiDpO2pGr/0a1X1mBwZ2XmMwAgAlzmKMgAALsD1wOjqGCPRJh4iDQCXu1pWBwAAwKUi/5zRxh+O6YsfjupQxhk1qh+grs0aSFKVPTC6MKNfHyLdKSa42PV3uSJEvlV93yQA4JJBUQYAuOzlnzOau+Z7zf/0B2Xn5nvMe2HtDwr08632GJJ3pmnskq3KyPYccv+FtT+otsNXs4a30y3tIqs9DgDAxUdRVkPlnzPatPe4jpw4owa1HZJN+uVkjhrUduicMdrw4y/6Of20bDabIur5q16AnzJO5yo144wi6vnL6ailA0dsSv1irzLP5OvQ/9o2qh+gLjEh8vGx6UjWGR0/lat6gX46fipHGafzZJNNsc1Cquyvuq5+pGWervC6Cm6L0CB/dYoJLnNsJW3H8i6nLMt25eXQ/3IQHOhQgyCHwut6rqssuT2UcUbhdR3yzbCp7zkjn2K2gSRt/OGYNvz4i/S/bXl902Bt3Zde5m1VUh8K71tZZ8qWr5Jy5bp64Yq1c0ywxz4YXOf8duoQXb9M8Re8GlIwVtc2D61zvi9f7j3m3jZVfaWiuL5KKhBXtrKP+qjeD8fUrUWYxzzXPlLati147PxyMkfHs3N1KP20JHn0N7i2nzKyz2/D0Drn96eCx7ZrXuHtW9rxUPhqU+FzTMHziGvfLbgfF7fc0pZZsF/FXVkq7r0Fc/3l3mN6+dMfdebsuRLzVbhQqw6vffFTifNO5eRr3NKv1WPLft3YPNSdr4LHfOG+FzxuzhnJGWB3b68LHZ+Fz70l7QMFj3lJRY6lks5fF1puwf2w8L7p2p9Ti9kPStoHXOsu7bh3xSqdP9Y+23NYW/b4aFfyHt3YIlRdrgiRpAp/nlxI4fOc65y8ee/xIjlMLcNV1MLnGNc2Lvh5mpFd9DgveMyVta8XOuZL2ydK2x6FP+sKn/Nd+Squnz+nn9RnP9u0K/m/8vHxVf1AP/f6S/qsKO1zufB5q7L5Lyk/lf2uUtbP7sJX5Kuyb2Vdl+t4zjidJ1No/w6v65At3aa6e37R5v0ZRc4xrnOXMSpTbr2BzRjDTexVJCsrS06nU5mZmapbt65lcSRtT9X0lTur7VabsrD7SL2uCtV10cHFfgkobZp0/qA7nZevlB+O6cSZsxVe1+HMHH28M02ncn79QlXH4aMbmjWQv93Xva7i4rrQ+gPtPmodWVeN6geWq18lxVUa17psNpt2pGaV+X2SVMtH8vXxUU6BL5z+tWw6e046e4HfsAT6+ahf6zCFOQOK9Our/ela/9+jOpNX8hfZ4pSUr+K2d6Cfj9o2qqttB7I84i+JTfL4zY+/3UfdmzfwWFdq5hl9czCzTMsryM/XpvaNncXmWyp5PyrrvmX3tcmY4nNSy+d87y6ULz9faUC7CJ3KPVeh3FxI4e1bkGtb13bY9cF3qeXeviXxt/voqrA6+s/hk+XqjytfNputQvn2Zn6+0rVN6pX5uCnI7iO1jqyrH37JLvbcV9o+UJKC+0ZJ572KLPdCXPtAZL2AMh/3FzoO7b4+Ol1gPyx8jpTKdy4oeD5d858jyssv/1ZwHfcFz9PFnWPKuo1dx1zhfaC486lrXeU935T2+Smdv523LJ91dl+bfGy2cu/nhbeF67Nm+6ET5fp8ddSyqV3jeurYtH6lvwMUF1NZ9q3itlVpcZU1XyXluzyffxXZNyqr8HaMcPpr2sCr1a9NxEWLoaDy1AYUZcV48cUXNWvWLKWmpqp169Z67rnndOONN17wfZdCUZa0PVVjF39VLT9GBwAAALyJTdK8u6+zpDArT23A6IuFvPXWW5o8ebKeeOIJff3117rxxhvVv39/7d+/3+rQLqg6RwcDAAAAvI2Rd4xwS1FWyOzZszV69Gjdf//9atWqlZ577jlFRUVp3rx5Vod2QZv2Hrf0lkUAAADgUuMa4fZSxkAfBeTm5mrr1q169NFHPabHxcUpJSWlSPucnBzl5OS4X2dlZUmS8vLylJeXV6R9dUvNOHXR1wkAAABc6lIzTikv7+L+vKg89QBFWQG//PKL8vPzFRYW5jE9LCxMaWlpRdrPnDlT06dPLzJ91apVCgwMrLY4S/Jjpk1S9Q/bDADwFOBrdDrfJqfd6C8df/3B/YLdPtp23EfDm+brk0M+ysiVzv/CAQBwMf24Y5s+PPj1RV1ndnZ2mdtSlBXDZvP8wDTGFJkmSY899pimTJnifp2VlaWoqCjFxcVZMtBH/jmjJX9bp/Tsi3+VDgAuZzOHtVPPqxrKJpsC/HyVl5en5ORkLfhdD/n41pLd10c9/3NUE5Z9Y3WoAHDZCa5t14Q7+1z04fFdd9GVBUVZAQ0aNJCvr2+Rq2JHjhwpcvVMkhwOhxwOR5Hpdrtddru92uIsiV3SM0PaaNzSi/tXAAC4nI25KUZDOjQpdl6dAH/358GAaxqrVi1fPbr8uyIPiAYAVJ+nB7eVv8Pvoq+3PPUAA30U4Ofnpw4dOig5OdljenJysrp27WpRVOVzS7tIjbkpxuowAMBr+dt99MCNMYpw+pfaro7DVy+OuFaP3XJ1mZfdr02Etv6pjx66uYUC/Uq+3bx+oF0P3Bij4NoX/0sEANQkY26K0S3trHlOWXnwnLJC3nrrLcXHx+ull15SbGysXn75Zb3yyivasWOHoqOjS33vpfCcMpcPv03VI29/q5M5ZX/wcXEPAy74oOXSHuZY1ocRV4WLua6Lzc/3/GX13Ao8QPRSUXjfOpyZU6UPE65KDl+b2v3vgdBWPOSyJLV8zt9GXdKDZF0P4C7rg5odtWzq0aKhxwPTy/pw9soq+LDwij5wvLDiHmoq6YIPCK7lI/nYbCUeX7UdvnqgW4wm9m4hXx+b8s8Zbdp7XGmZp/XLyRwdz85VasYZNaofoK7NGqjLFSEl3gqTl5enDz/8ULfcckuJfynNP2e08Ydj+uKHo/r5fw9bLbxsVwyrdqTq31/9XGy+/HxtOmeMStsNinvAcFUdn65jPrJewEV9YGxFHw5f1uO+lo907pzRuUvkN4AFHw5f0QdNXwwF81JVx3xF+cjIXsu3Uvt3TfhcLk1FH8Bd1etyPSz7VO65IvtM4XNMSQ/gLqyOw1fPDm+nW9pFVkt/yoKHR1fSiy++qGeffVapqalq06aN5syZo5tuuumC77uUijLp1w/8DT/+onNGqh/opwZBDoXX9VenmOBiv0y4vgCkZpzSjzu2acKd/Twu9xb8knL8VK6C6/y6PEnuLxiHMs4U+3T3skyTfn06fHCgQ8G1/ZSRXfF1ZZ3Jk/lf/4Nr++n4qfNfri70JPqC6z9+KkcZp/Nkk02d/7f+DT/+4v4iVd5+FY6rYF6k8483OHLijBrUduicMR7ralQ/QF1iQiRJX+49Jv0vJh8fm345maPQIH91iK6vTT8eVdKnX6rxlVfpRM65ItugYH98fGw6knVGx0/lql7gr/ONkZwB9hL7lfq/bR8c6Ch13yr45bO0fBXe3q71Z53xjPWXkzlqUNsh2eQRd0b2r/G7vkAXt66SvlQXPGYkm2Kbhej6psHavPe4xxfnC+W2LPt8wb4W92U/Ly9Pc99KUtOr2ys9+6w7X7HNQjziLnycu/JVluKhuOM5tI7ndnVNO2eMvtx7zONc4pq+4cdfPPLqylfhWAuus+D+XZblXmgfK25fK67QKXx8ySb3cVPacsurLEVZeZXl/FtwPygtD4WXW/j4LHzuLWkfKO8xX/i8V9q+Vdx+WJZ8Fd4HJJXpuHfF6mrXIaquPvroI4Vc1Vmb9mW44yvuHFWZc0HB86lr3a7zTsFzUVmO+7J8fhU837v2oQ7R9d3nucL7wIXOpxc6Pos75kv7/CyYr8KfdYXPya58udq5zhmN6geoU3Q9Hf/Pl+rfv7++Pnii2HOO6zgxhf5f3s/lyn4HcB1rpX3+lrZvFd4GF4qrcL5c/Sv8B6iK9Ks86yp8Hiu4f2/4/ohWffal4m7srNgrQ0v9DCuc2+KOGatQlFnkUivKKqM6vlDg4iKH3o38eTfy5/3IoXcjf96tpuSvPLUBvykDAAAAAAtRlAEAAACAhSjKAAAAAMBCFGUAAAAAYCGKMgAAAACwEEUZAAAAAFioltUB1CSupwtkZWVZHEnl5eXlKTs7W1lZWV49FOnljBx6N/Ln3cif9yOH3o38ebeakj9XTVCWJ5BRlFWhEydOSJKioqIsjgQAAADApeDEiRNyOp2ltuHh0VXo3LlzOnTokIKCgmSzWfsE8crKyspSVFSUDhw44PUPwr5ckUPvRv68G/nzfuTQu5E/71ZT8meM0YkTJxQZGSkfn9J/NcaVsirk4+Ojxo0bWx1Glapbt65XHwwgh96O/Hk38uf9yKF3I3/erSbk70JXyFwY6AMAAAAALERRBgAAAAAWoihDsRwOh6ZNmyaHw2F1KKggcujdyJ93I3/ejxx6N/Ln3S7H/DHQBwAAAABYiCtlAAAAAGAhijIAAAAAsBBFGQAAAABYiKIMAAAAACxEUYYyGTRokJo0aSJ/f39FREQoPj5ehw4dsjoslMFPP/2k0aNHKyYmRgEBAWrWrJmmTZum3Nxcq0NDGT3zzDPq2rWrAgMDVa9ePavDQRm8+OKLiomJkb+/vzp06KDPPvvM6pBQRp9++qkGDhyoyMhI2Ww2vfPOO1aHhHKYOXOmrr/+egUFBSk0NFRDhgzR7t27rQ4LZTRv3jy1a9fO/dDo2NhYffTRR1aHdVFQlKFMevbsqX/+85/avXu33n77bf3www+67bbbrA4LZfCf//xH586d0/z587Vjxw7NmTNHL730kh5//HGrQ0MZ5ebm6vbbb9fYsWOtDgVl8NZbb2ny5Ml64okn9PXXX+vGG29U//79tX//fqtDQxmcOnVK7du319y5c60OBRWwfv16jR8/Xhs3blRycrLOnj2ruLg4nTp1yurQUAaNGzfWX//6V23ZskVbtmxRr169NHjwYO3YscPq0KodQ+KjQt577z0NGTJEOTk5stvtVoeDcpo1a5bmzZunH3/80epQUA4LFy7U5MmTlZGRYXUoKEXnzp113XXXad68ee5prVq10pAhQzRz5kwLI0N52Ww2rVixQkOGDLE6FFTQ0aNHFRoaqvXr1+umm26yOhxUQHBwsGbNmqXRo0dbHUq14koZyu348eNasmSJunbtSkHmpTIzMxUcHGx1GECNk5ubq61btyouLs5jelxcnFJSUiyKCrh8ZWZmShKfeV4oPz9fy5Yt06lTpxQbG2t1ONWOogxl9sc//lG1a9dWSEiI9u/fr3fffdfqkFABP/zwgxITE/W73/3O6lCAGueXX35Rfn6+wsLCPKaHhYUpLS3NoqiAy5MxRlOmTFG3bt3Upk0bq8NBGX333XeqU6eOHA6Hfve732nFihW6+uqrrQ6r2lGUXcYSEhJks9lK/bdlyxZ3+z/84Q/6+uuvtWrVKvn6+uqee+4Rd79ap7z5k6RDhw6pX79+uv3223X//fdbFDmkiuUP3sNms3m8NsYUmQagek2YMEHffvut3nzzTatDQTm0bNlS27Zt08aNGzV27FiNHDlSO3futDqsalfL6gBgnQkTJug3v/lNqW2aNm3q/n+DBg3UoEEDtWjRQq1atVJUVJQ2btx4WVxSvhSVN3+HDh1Sz549FRsbq5dffrmao8OFlDd/8A4NGjSQr69vkatiR44cKXL1DED1mThxot577z19+umnaty4sdXhoBz8/Px05ZVXSpI6duyozZs36/nnn9f8+fMtjqx6UZRdxlxFVkW4rpDl5ORUZUgoh/Lk7+eff1bPnj3VoUMHLViwQD4+XCS3WmWOP1y6/Pz81KFDByUnJ2vo0KHu6cnJyRo8eLCFkQGXB2OMJk6cqBUrVmjdunWKiYmxOiRUkjHmsvi+SVGGC9q0aZM2bdqkbt26qX79+vrxxx/15z//Wc2aNeMqmRc4dOiQevTooSZNmujvf/+7jh496p4XHh5uYWQoq/379+v48ePav3+/8vPztW3bNknSlVdeqTp16lgbHIqYMmWK4uPj1bFjR/eV6f379/M7Ti9x8uRJff/99+7Xe/fu1bZt2xQcHKwmTZpYGBnKYvz48Vq6dKneffddBQUFua9aO51OBQQEWBwdLuTxxx9X//79FRUVpRMnTmjZsmVat26dkpKSrA6t2jEkPi7ou+++0+9//3t98803OnXqlCIiItSvXz/96U9/UqNGjawODxewcOFC3XvvvcXO4/D3DqNGjdLrr79eZPratWvVo0ePix8QLujFF1/Us88+q9TUVLVp00Zz5sxhOG4vsW7dOvXs2bPI9JEjR2rhwoUXPyCUS0m/3VywYIFGjRp1cYNBuY0ePVqffPKJUlNT5XQ61a5dO/3xj39Unz59rA6t2lGUAQAAAICF+GEJAAAAAFiIogwAAAAALERRBgAAAAAWoigDAAAAAAtRlAEAAACAhSjKAAAAAMBCFGUAAAAAYCGKMgAAAACwEEUZAAAAAFiIogwAgCry7bffatiwYQoJCZG/v79at26tWbNm6ezZs1aHBgC4hFGUAQBQBdavX68uXbooICBA7777rr755hs98sgj+vvf/65hw4bp3LlzVocIALhE2YwxxuogAADwZvn5+WrevLm6du2qxYsXe8zbuXOnrrnmGs2bN0+jR4+2KEIAwKWMogwAgErasGGDunbtqm3btql9+/ZF5g8ZMkTZ2dlatWqVBdEBAC513L4IAEAl7d27V5LUvHnzYue3aNFC+/btu5ghAQC8CEUZAACVVLduXUnS8ePHi52fnp7ubgMAQGEUZQAAVFJsbKzsdrtWrlxZZF5+fr5WrVqlbt26WRAZAMAbUJQBAFBJISEhmjRpkp5++mkdOnTIY96cOXN07NgxPfTQQxZFBwC41FGUAQBQSSdPntSkSZMUExOjnj176quvvpIkzZo1S48//rgSExPl5+en/Px8iyMFAFyKGH0RAIBKSkhI0PTp092vR44cqYULF8pms3m027t3r5o2bXqRowMAXOooygAAAADAQty+CAAAAAAWoigDAAAAAAtRlAEAAACAhSjKAAAAAMBCFGUAAAAAYCGKMgAAAACwEEUZAAAAAFiIogwAAAAALERRBgAAAAAWoigDAAAAAAtRlAEAAACAhf4/KZ0iWpKj/lcAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "N = 64 # length of x\n", "M = 512 # length of y\n", "\n", "# generate two uncorrelated random signals\n", "np.random.seed(1)\n", "x = 2 + np.random.normal(size=N)\n", "y = 3 + np.random.normal(size=M)\n", "N = len(x)\n", "M = len(y)\n", "\n", "# compute cross PSD via CCF\n", "acf = 1 / N * np.correlate(x, y, mode=\"full\")\n", "psd = np.fft.fft(acf)\n", "psd = psd * np.exp(1j * np.arange(N + M - 1) * 2 * np.pi * (M - 1) / (2 * M - 1))\n", "psd = np.fft.fftshift(psd)\n", "Om = 2 * np.pi * np.arange(0, N + M - 1) / (N + M - 1)\n", "Om = Om - np.pi\n", "\n", "# plot results\n", "plt.figure(figsize=(10, 4))\n", "plt.stem(Om, np.abs(psd), basefmt=\"C0:\")\n", "plt.title(\"Biased estimator of cross power spectral density\")\n", "plt.ylabel(r\"$|\\hat{\\Phi}_{xy}(e^{j \\Omega})|$\")\n", "plt.xlabel(r\"$\\Omega$\")\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*." ] } ], "metadata": { "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.11.0" } }, "nbformat": 4, "nbformat_minor": 1 }