From 4945d29a5c1dc97d0fdee6ffe8d5671eaf7b5552 Mon Sep 17 00:00:00 2001 From: Ryan Hammonds Date: Tue, 27 Aug 2024 18:33:57 -0700 Subject: [PATCH] derivations for acf->psd, the spectral ar form, and lorentzian bias --- docs/examples/derivations_ar_acf_psd.ipynb | 363 +++++++++++++++++++++ 1 file changed, 363 insertions(+) create mode 100644 docs/examples/derivations_ar_acf_psd.ipynb diff --git a/docs/examples/derivations_ar_acf_psd.ipynb b/docs/examples/derivations_ar_acf_psd.ipynb new file mode 100644 index 0000000..92b3080 --- /dev/null +++ b/docs/examples/derivations_ar_acf_psd.ipynb @@ -0,0 +1,363 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "id": "21184777-9849-4d6b-ae61-95a86aa2102b", + "metadata": {}, + "outputs": [], + "source": [ + "import matplotlib.pyplot as plt\n", + "import numpy as np" + ] + }, + { + "cell_type": "markdown", + "id": "688fdf61-a4e3-4622-8c7f-41f9dba3d2cc", + "metadata": {}, + "source": [ + "## Derivation of the AR(1) PSD from the Fourier Transform of an Exponentially Decaying ACF\n", + "\n", + "First, validate the AR(1) PSD as the spectral representation of the exponentially decaying ACF. Note that and exponentially decaying ACF, $A(t)$, can be represented as the AR(1) coefficient, $\\varphi$, which is the first lagged correlation of the ACF.\n", + "\n", + "$$\n", + "\\begin{align}\n", + "A(t) &= e^{\\frac{-t}{\\tau f_s}} \\\\\n", + "\\tau &= -\\frac{1}{\\ln(\\varphi) f_s} \\\\\n", + "A(t) &= \\varphi^t\n", + "\\end{align}\n", + "$$" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "83dfb16e-1188-4256-a8da-060aae6271c9", + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiMAAAGzCAYAAAD9pBdvAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8g+/7EAAAACXBIWXMAAA9hAAAPYQGoP6dpAABHSElEQVR4nO3deXhU5cH+8e+ZNQvZIGQPhB1klyVEULFGUHGpaIuIQnFXVJCfG1XxtVWxr2KxiqJWW2uxuFStK8obFaXsgSD7IkvCko2QheyZOb8/opEIwQSSnExyf65rrjpnnnPmztNezt0z5zxjmKZpIiIiImIRm9UBREREpG1TGRERERFLqYyIiIiIpVRGRERExFIqIyIiImIplRERERGxlMqIiIiIWEplRERERCylMiIiIiKWUhkRERERS6mMiEije+GFFzAMg8TExDrHZGVlcc8999C7d28CAgIIDAxkyJAhPPbYY+Tn59eMGz16NIZhnPCxbdu2ZvhrRKSpOawOICKtz8KFC0lISGD16tXs2rWL7t2713p9zZo1XHzxxRw9epRrr72WIUOGALB27VqefPJJvvnmG7744oua8XFxccyZM+e494mJiWnaP0REmoXKiIg0qj179rB8+XLee+89brnlFhYuXMgjjzxS83p+fj5XXHEFdrud9evX07t371r7P/7447zyyiu1toWEhHDttdc2S34RaX76mkZEGtXChQsJCwtj3LhxXHXVVSxcuLDW6y+99BIHDhzgmWeeOa6IAERGRvLQQw81V1wRaQFURkSkUS1cuJDx48fjcrmYOHEiO3fuZM2aNTWvf/jhh/j7+3PVVVfV+5gej4fc3Nxaj6NHjzZFfBGxgMqIiDSa1NRUtm3bxtVXXw3AqFGjiIuLq3V2ZOvWrfTs2ROXy1Xv427bto2OHTvWetxxxx2Nnl9ErKFrRkSk0SxcuJDIyEjOO+88AAzDYMKECfzzn/9k7ty52O12CgsLCQoKatBxExISjruORBevirQeKiMi0ig8Hg+LFi3ivPPOY8+ePTXbExMTmTt3LikpKYwZM4bg4GCKiooadOzAwECSk5MbO7KItBAqIyLSKL788ksOHTrEokWLWLRo0XGvL1y4kDFjxtC7d2/S0tKoqKho0Fc1ItJ6qYyISKNYuHAhERERzJ8//7jX3nvvPd5//30WLFjApZdeyooVK/j3v//NxIkTLUgqIi2NYZqmaXUIEfFtpaWlREZG8pvf/IZXX331uNeXL1/OyJEjWbRoEWPGjKFv374AfP311/Ts2bPW2OzsbF5++eWa23tHjx5Nbm4umzZtavo/REQsoTMjInLaPvzwQ4qKirjssstO+PqIESPo2LEjCxcuZMKECbz//vtcfPHFDBo0qNYKrOvWreNf//oXSUlJzRlfRCymMiIip23hwoX4+flxwQUXnPB1m83GuHHjWLhwIYcPHyYxMZFNmzbx1FNP8cknn/DGG29gs9no06cPDzzwgG7bFWlj9DWNiIiIWEqLnomIiIilVEZERETEUiojIiIiYimVEREREbGUyoiIiIhYSmVERERELOUT64x4vV4OHjxIUFAQhmFYHUdERETqwTRNioqKiImJwWar+/yHT5SRgwcPEh8fb3UMEREROQUZGRnExcXV+bpPlJGgoCCg+o8JDg62OI2IiIjUR2FhIfHx8TWf43XxiTLy41czwcHBKiMiIiI+5pcusdAFrCIiImIplRERERGxlMqIiIiIWMonrhkRERFpCUzTpKqqCo/HY3WUFsFut+NwOE572Q2VERERkXqoqKjg0KFDlJSUWB2lRQkICCA6OhqXy3XKx1AZERER+QVer5c9e/Zgt9uJiYnB5XK1+UU4TdOkoqKCnJwc9uzZQ48ePU66sNnJqIyIiIj8goqKCrxeL/Hx8QQEBFgdp8Xw9/fH6XSyb98+Kioq8PPzO6Xj6AJWERGRejrV/+ffmjXGnGhWRURExFIqIyIiImKpBpeRb775hksvvZSYmBgMw+CDDz74xX2+/vprzjzzTNxuN927d+fvf//7KUQVERGR5nTkyJFmeZ8Gl5Hi4mIGDhzI/Pnz6zV+z549jBs3jvPOO4+0tDRmzJjBjTfeyOeff97gsCIiItJ87r333mZ5nwbfTXPRRRdx0UUX1Xv8ggUL6NKlC3PnzgWgT58+LFu2jD//+c+MHTv2hPuUl5dTXl5e87ywsLChMevl9eV72XKwkFtHd6NLeGCTvIeIiIgv+uijj9i1axcLFizg1ltvbdL3avJrRlasWEFycnKtbWPHjmXFihV17jNnzhxCQkJqHvHx8U2S7dDKt+iZ9gRZm75ukuOLiIj4go0bN3LJJZfUPC6//HI6dOjA9ddf3+RFBJphnZHMzEwiIyNrbYuMjKSwsJDS0lL8/f2P22fWrFnMnDmz5nlhYWGTFJJkcxVDHSmszOgJjGv044uIiPiC/v378/HHH9fa9tJLLzF8+PBmef8WeTeN2+0mODi41qMpVAZ3AsDI39skxxcREfFVHTt25LnnniMjI6PJ36vJz4xERUWRlZVVa1tWVhbBwcEnPCvSnOzhXeEABB5NtzSHiIhISzN+/HjGjx/fLO/V5GdGkpKSSElJqbVtyZIlJCUlNfVb/6J2UT0ACKs4aHESERER6xiGUe9HU2hwGTl69ChpaWmkpaUB1bfupqWlkZ5efXZh1qxZTJ48uWb8rbfeyu7du7nvvvvYtm0bL7zwAm+//TZ333134/wFpyG8U28AIr05VFaU/8JoERER35Wens4111xDWFgY7du3Z9KkSTXriJimWe9HU2hwGVm7di2DBw9m8ODBAMycOZPBgwcze/ZsAA4dOlRTTAC6dOnCJ598wpIlSxg4cCBz587lr3/9a5239Tan8KhOlJtOHIaXnP3fWx1HRER8iGmalFRUWfJoaCnYtWsXQ4YMoXv37qxcuZIlS5awa9euk64j8sgjj3DfffeRmZl5ulP1ixp8zcjo0aNPOgknWl119OjRrF+/vqFv1eRsdjuZ9ig6ezM4vH87MV3PsDqSiIj4iNJKD2fMtmYBzy1/GEuAq/4f4bfffju33347jz76aM22++67r6aMLFu2jGeffZbNmzczZMgQbr75Zj755BPuuusuoqKiGj3/zzX5Bawt3ctRj/L57nJmOAfT3+owIiIijWzfvn0sWbKEZcuW1SxACuDxeGqWzRg1ahQBAQEsX76cO+64g6ysLC6++OJal100pTZfRpxRvcndvZeMI6VWRxERER/i77Sz5Q/WXHLg77TXe+yGDRto3749q1atOv44x9zVumXLFvr27Vuzz4ABA3jnnXfIzs7mggsuoGfPnqcfvA5tvox07hAAwL7DJRYnERERX2IYRoO+KrGK0+mkqKiImJgYAgIC6hwXFhbG22+/zXnnnceGDRu4/PLLyczMZPPmzcTGxjZpxha56Flz6umXz4OOfzJ2/7NWRxEREWl0iYmJBAcHM3nyZDZs2MCuXbtYvHgxM2bMqDVu3LhxvPjii0D1D+T17NmTLVu2UFFRgdfrbdKMLb/SNbG4dgYjHZ9ytMwf0+vFsLX5fiYiIq1I+/bt+fTTT7n//vs555xzME2THj16MGXKlF/ctzl+lwZURojs1AOvadDOKCUv9xDtI5r2VJSIiEhzGz58OF999ZXVMerU5k8D+PkHkmO0ByA7fZvFaURERNqeNl9GAHJdMQAcPbjT4iQiIiJtj8oIUBxQfZ91Ve5ui5OIiIi0PSojgCe0MwD2gr3WBhEREWmDVEYAZ8duAPiVNv36+yIiIlJbm7+bBsDd+yKGfeuPzRnJ8evTiYiISFNSGQHioiLIIQyKKiir9ODXgGV2RURE5PToaxogNMBJkF91L0vP07LwIiIizUlnRqj+fYFpASl09qyjcJsJkZdYHUlERKTN0JmRHwy3beMi+xqqDqy3OoqIiEibojLyg/Kg6tt7jbw9FicRERFpfrt37+bDDz+05L1VRn5g79AFAP+j6RYnERERaX6fffYZW7ZsseS9VUZ+EBDVHYDQioMWJxEREWleS5cu5eGHH+bVV19l8ODBFBcXN+v7q4z8oEN8LwCiPFl4qqosTiMiItJ8zj33XAYMGMCSJUtYv349gYGBzfr+KiM/iIjtRoVpx2VUkXNQ142IiEjbkp6eTkJCgiXvrTLyA7vDQZYtkjLTSc5BXTciIiJtx/79+4mJibHs/VVGjjEn5jn6lP+NLfYeVkcRERFpVJmZmUyaNIm///3v3H777QQHB3PJJZdQUFDAvn37iI6OtiybysgxOnSMwsTGvsNahVVEROqporjuR2VZA8aW1m/sKSgpKeH888+nrKyMkJAQXn/9dT766COWL1/OE088Qb9+/di9ezf9+/e35I4arcB6jM4dAgAtCS8iIg3wxEm+3ugxBia989Pzp7pDZR2fMZ1HwdRPfno+rz+UHD5+3P8UNDji3LlzKSsr45///CdPPPEEZ599Nueeey5jxowhLS2NkJAQUlNTG3zcxqIycow+joO84JyH/75A4AOr44iIiDSKd955h9/97nf4+/uzfv16Bg4cCFT/HEpQUJDF6VRGaokN9WOkfTWF5QGYXi+GTd9iiYjIL/j9SdanMn72K/D37jrJ2J995szYeOqZjlFaWsrGjRt57LHHAEhLS2PSpEkArFu3jltuuaVR3ud0qIwcIyqhD17TINgo4XDOQTpExlkdSUREWjpXA9bkaKqxJ1H1w9pZ5eXl5OTkcODAAQYOHMh///tf9uzZw/jx4xvlfU6Hysgx/PwDOWiLIMbMImvPJpURERHxeUFBQQwaNIi//OUvFBUV4XQ62blzJ9OmTeOJJ56wbG2RY+l7iJ/JdccDcPTAVouTiIiINI433niD4uJibrrpJqqqqpgxYwazZs3innvusToaoDMjxykJ6gJla/Hm7LQ6ioiISKPo168f69at45prriEoKIiXXnrJ6ki16MzIzxjh1T+Y5y7UkvAiItK6bNq0iQEDBlgd4zgqIz8TEN2HctNBaXml1VFEREQaTVVVFdu3b6d///5WRzmOvqb5mQ79fkWfT/+OvcrOVo8Xh119TUREfJ/D4aC8vNzqGCekT9qfiQpth8vpoNJjsv9I6S/vICIiIqdFZeRnbDaDLuHtANide9TiNCIiIq2fysgJXG//jI9dv8e9/m9WRxEREWn1VEZOIN6vhH62vThym/+XC0VERNoalZETsHfsCUC7It3eKyIiPzFN0+oILU5jzInKyAmExJ8BQMeK/RYnERGRlsDpdAJQUlJicZKW58c5+XGOToVu7T2ByC79AIggj6OFR2gXHGZxIhERsZLdbic0NJTs7GwAAgICMAzD4lTWMk2TkpISsrOzCQ0NxW63//JOdVAZOYGQsHAOE0IHCsjcs5nuA0dZHUlERCwWFRUFUFNIpFpoaGjN3JwqlZE6ZDvj6FBZQH7GFlAZERFp8wzDIDo6moiICCortUo3VH81czpnRH6kMlKH/KAefJ+bR3ZhmdVRRESkBbHb7Y3yASw/0QWsddgwcDbnV8xlsXG21VFERERaNZWROnQNDwS0CquIiEhTUxmpQ9eO1WVkT85RTK/X4jQiIiKtl8pIHTqF+fOe6xFWGNeTm5ludRwREZFWS2WkDi6ngwj7UYKNErL3bLY6joiISKulMnISh/3iATh6cKvFSURERFovlZGTKA3uCoA3d5fFSURERFovlZGTMDp0ByCgcLfFSURERFovlZGTCIztDUD7sgyLk4iIiLReKiMnEZHQH4BobyYV5VqJVUREpCmojJxEx+hObDM78Y13APszs6yOIyIi0irpt2lOwrDZ+H8dXmDzwUJeKfajq9WBREREWqFTOjMyf/58EhIS8PPzIzExkdWrV590/Lx58+jVqxf+/v7Ex8dz9913U1bmG197dOvYDoBd2VoWXkREpCk0uIy89dZbzJw5k0ceeYR169YxcOBAxo4dS3Z29gnHv/nmmzzwwAM88sgjbN26lVdffZW33nqL3//+96cdvjn0iGgHmOzT1zQiIiJNosFl5JlnnuGmm25i6tSpnHHGGSxYsICAgABee+21E45fvnw5I0eO5JprriEhIYExY8YwceLEXzyb0lIMt21lnfsWbth5p9VRREREWqUGlZGKigpSU1NJTk7+6QA2G8nJyaxYseKE+5x11lmkpqbWlI/du3fz6aefcvHFF9f5PuXl5RQWFtZ6WCU6rjPtjaPEVmXg9XgsyyEiItJaNaiM5Obm4vF4iIyMrLU9MjKSzMzME+5zzTXX8Ic//IFRo0bhdDrp1q0bo0ePPunXNHPmzCEkJKTmER8f35CYjSomoQ/lppMAo5xD+3ZYlkNERKS1avJbe7/++mueeOIJXnjhBdatW8d7773HJ598wh//+Mc695k1axYFBQU1j4wM6xYdczhdHLDHApCzO82yHCIiIq1Vg27tDQ8Px263k5VV+2LOrKwsoqKiTrjPww8/zHXXXceNN94IQP/+/SkuLubmm2/mwQcfxGY7vg+53W7cbndDojWpvMBudC3aS+mBTVZHERERaXUadGbE5XIxZMgQUlJSarZ5vV5SUlJISko64T4lJSXHFQ673Q6AaZoNzWuJyg49AXAe3m5xEhERkdanwYuezZw5kylTpjB06FCGDx/OvHnzKC4uZurUqQBMnjyZ2NhY5syZA8Cll17KM888w+DBg0lMTGTXrl08/PDDXHrppTWlpKXzi+kHeyG0WD+YJyIi0tgaXEYmTJhATk4Os2fPJjMzk0GDBrF48eKai1rT09NrnQl56KGHMAyDhx56iAMHDtCxY0cuvfRSHn/88cb7K5pYh25n8u23/djq7UZXr4nNZlgdSUREpNUwTB/4rqSwsJCQkBAKCgoIDg5u9vf3eE36zF5MRZWXpfeOpnOHwGbPICIi4mvq+/mtH8qrB7vNqFkWfkeWloUXERFpTCoj9dQzsh3BFLM/XdeNiIiINCaVkXq6qvIjvvO7if5b5lodRUREpFVRGamndpFdAd1RIyIi0thURuopvOtAAOKq0vFUVVmcRkREpPVQGamnmIQ+lJlO/IxKDu3dZnUcERGRVkNlpJ7sDgcHHNU/2KffqBEREWk8KiMNcCSw+rqR8oNbLE4iIiLSeqiMNEBlh14AOPL0GzUiIiKNRWWkAexdRvHPqvNZ4hlidRQREZFWo8G/TdOWRfY9l99+auIusHGf18Su36gRERE5bToz0gBxYQG4HTbKq7yk55VYHUdERKRVUBlpALvNoG9HB/2M3aTv0XUjIiIijUFlpIF+73mJj90P4dr8ttVRREREWgWVkQb68Y4aZ95Oi5OIiIi0DiojDeQf2w+AsOLvLU4iIiLSOqiMNFDHroMAiKvaT1VlhbVhREREWgGVkQaK7tyTYtMPt1HJ/l3fWR1HRETE56mMNJDNbifdVb0sfO7OtRanERER8X0qI6egMKQ3AFUHdWZERETkdGkF1lNwtNs4Hj/kR5mZxAirw4iIiPg4nRk5BeH9knnFcwmf5kZhmqbVcURERHyaysgp6BUVhM2Aw8UV5BSVWx1HRETEp6mMnAI/p52RHYq41LacPdvXWx1HRETEp+makVM0w3iLIa4vWbnFhGFJVscRERHxWTozcooqOlavxOrM3WJxEhEREd+mMnKKAjsNBKBj8Q6Lk4iIiPg2lZFTFNN7OABxnoOUHC2wOI2IiIjvUhk5ReFRncglFJthkrEt1eo4IiIiPktl5DQc8usGQMGedRYnERER8V0qI6fhaFgfAMzMjRYnERER8V26tfc0lPa+ipv3dcRkMIlWhxEREfFRKiOnoVOfoXzxWQkBOXa8XhObzbA6koiIiM/R1zSnIaFDIG6HjZIKD/vySqyOIyIi4pNURk6Dw27jig7p3GF/n0ObvrY6joiIiE9SGTlN423fcI/zHWw7v7A6ioiIiE9SGTlNZlR/APzztlqcRERExDepjJymkITBAESX7rQ4iYiIiG9SGTlNcX2GARBBHvk5hyxOIyIi4ntURk5Tu+D27DeiADiwbbXFaURERHyPykgjyA7oCUDRvvUWJxEREfE9KiONoLxjXwBs2ZstTiIiIuJ7tAJrYxh0DWO3d6bK05MUq7OIiIj4GJWRRtC7R2+2mxlwuIyC0kpC/J1WRxIREfEZ+pqmEYQFuujUPgCAjfsLLE4jIiLiW3RmpJFMCNtOVNHHVKzZCj3usTqOiIiIz9CZkUYyxD+TK+3fEnbgK6ujiIiI+BSVkUYS3HU4ADHFWhZeRESkIVRGGkmnfkl4TYNIDpObmW51HBEREZ+hMtJI2gWHkW6PA+DApv9anEZERMR3qIw0opyg6sXPSveusTiJiIiI71AZaUTemOpf8A3I/c7iJCIiIr5DZaQRhXVPxGsaVJYVY5qm1XFERER8gtYZaUSd+5/Fme+9Sr7Hj2+PlBL/w0JoIiIiUjedGWlEbpeb+KhIAL7TSqwiIiL1ojLSyAbEhQDw3f58a4OIiIj4iFMqI/PnzychIQE/Pz8SExNZvXr1Scfn5+czbdo0oqOjcbvd9OzZk08//fSUArd0owP38a7rfxi74S6ro4iIiPiEBl8z8tZbbzFz5kwWLFhAYmIi8+bNY+zYsWzfvp2IiIjjxldUVHDBBRcQERHBu+++S2xsLPv27SM0NLQx8rc43WI60tW2g6Iyf7weDza73epIIiIiLZphNvC2j8TERIYNG8bzzz8PgNfrJT4+njvvvJMHHnjguPELFizgqaeeYtu2bTidzlMKWVhYSEhICAUFBQQHB5/SMZpLVWUFlY/F4m9UsG/i13TuNdjqSCIiIpao7+d3g76mqaioIDU1leTk5J8OYLORnJzMihUrTrjPhx9+SFJSEtOmTSMyMpJ+/frxxBNP4PF46nyf8vJyCgsLaz18hcPpYq+rOwDZ2048JyIiIvKTBpWR3NxcPB4PkZGRtbZHRkaSmZl5wn12797Nu+++i8fj4dNPP+Xhhx9m7ty5PPbYY3W+z5w5cwgJCal5xMfHNySm5QrC+gPg3Z9qcRIREZGWr8nvpvF6vURERPDyyy8zZMgQJkyYwIMPPsiCBQvq3GfWrFkUFBTUPDIyMpo6ZqNyxA8BIPTIJouTiIiItHwNuoA1PDwcu91OVlZWre1ZWVlERUWdcJ/o6GicTif2Yy7k7NOnD5mZmVRUVOByuY7bx+1243a7GxKtRYnscxakQkLl91RWlON0+e7fIiIi0tQadGbE5XIxZMgQUlJSarZ5vV5SUlJISko64T4jR45k165deL3emm07duwgOjr6hEWkNYjr2pfdxPK1dyA79vrWWR0REZHm1uCvaWbOnMkrr7zC66+/ztatW7ntttsoLi5m6tSpAEyePJlZs2bVjL/tttvIy8tj+vTp7Nixg08++YQnnniCadOmNd5f0cIYNhuPJbzOLZUzWZmtW3tFREROpsHrjEyYMIGcnBxmz55NZmYmgwYNYvHixTUXtaanp2Oz/dRx4uPj+fzzz7n77rsZMGAAsbGxTJ8+nfvvv7/x/ooWaEjnML7cls26fUe4YVQXq+OIiIi0WA1eZ8QKvrTOyI9W7j7M1S+vYGBQIR/MuhrDppX3RUSkbanv57d+tbeJDIxpxzL3dOIqczmYPoSYhF5WRxIREWmR9H/Xm4i/n5sSZ3sADn73lcVpREREWi6VkSaU1756Kfiq9JUWJxEREWm5VEaakKtL9e3OHY+kWRtERESkBVMZaUKdBv4KgISqvRQV5FmcRkREpGVSGWlC4TGdOWhEYjdM9qZ9bXUcERGRFkllpIkdCBoAwNFd/7U4iYiISMukW3ubWHG3cfx1jYNDZb058YL5IiIibZvKSBOLHH4lv1sRQWCWnVkeLw67TkaJiIgcS5+MTaxnZBBBbgfFFR62ZRZZHUdERKTFURlpYnabQWKcmyTbZvZ+963VcURERFoclZFmcIP9M/7lepzoLa9aHUVERKTFURlpBkE9zgIgtmijxUlERERaHpWRZtBl0LlUmTaiyCFr//dWxxEREWlRVEaaQWBQKHsdXQDYv0E/miciInIslZFmcjhsEABVe5ZbG0RERKSFURlpJo6uIwGIyEu1OImIiEjLojLSTBKGjAWgi3cvedkHLE4jIiLScqiMNJMOkXHM9Z/OxeVPsDLT6jQiIiIth8pIMyrqM4EtZgLLd+dZHUVERKTFUBlpRiO6dgBgxfeHLU4iIiLScuiH8prRiK7tGW//lpH5m8g91JXw6E5WRxIREbGczow0o9AAF7f7fcGV9m/Zm7rY6jgiIiItgspIM8sNHw6A9/ulFicRERFpGVRGmpl/r9EAxOSvtTaIiIhIC6Ey0sy6DhlDlWkjzswkM32n1XFEREQspzLSzIJC2rPb2R2AjHWfW5xGRETEeiojFjgcMaL6H/Z8a20QERGRFkBlxALtep0HQGVRDqZpWpxGRETEWlpnxALdho0l8YuXyPIEsTSvhM4dAq2OJCIiYhmdGbFAQEAg8XHVC55pNVYREWnrVEYscla36qXhV36fbXESERERa+lrGoucE1VOkvMxOm3PxfRux7CpF4qISNukT0CL9OvRjTNtO4klm/QdaVbHERERsYzKiEX8/APZ4dcPgEPrPrE4jYiIiHVURixU3Kn6Ft+A9K+tDSIiImIhlRELRZ95CQC9SjdQWlxkcRoRERFrqIxYqFOvwWTSEbdRyY5Vi62OIyIiYgmVEQsZNhv72icBULZVZURERNomlRGLOc64hC89g/iiMN7qKCIiIpbQOiMW6zlqPIO/DKKqwOS63GISwrU0vIiItC06M2KxID8nQxPCAFi6I8fiNCIiIs1PZaQFGN0rgmgOk5em9UZERKTt0dc0LcAFHQu41e9OSrNclJX+Dj9/fVUjIiJth86MtABdew8im/b4GxXsXP251XFERESalcpIC2DYbOwNHQFA8RaVERERaVtURloIR68xAETnLLM4iYiISPNSGWkhuo24lCrTRmfvfg7u3W51HBERkWajMtJChISFs9PVB4CM1R9anEZERKT5qIy0IPmx5wDg2vulxUlERESaj27tbUHCR1zDnTvsLK86k28rqghw6b8eERFp/XRmpAXp3qs/60PO53CVH9/syLU6joiISLNQGWlBDMNgbN8oAL7YnGlxGhERkeah7wFamIt6heBe8QHnbt1CZcWXOF1uqyOJiIg0KZ0ZaWEGd4nkRudiEtnE9lWLrY4jIiLS5FRGWhi7w8GusLMBKN7wgbVhREREmsEplZH58+eTkJCAn58fiYmJrF69ul77LVq0CMMw+PWvf30qb9tmuPpdBkCX3K/xerwWpxEREWlaDS4jb731FjNnzuSRRx5h3bp1DBw4kLFjx5KdnX3S/fbu3cs999zD2Weffcph24reZ11KselHBHns2vCt1XFERESaVIPLyDPPPMNNN93E1KlTOeOMM1iwYAEBAQG89tprde7j8XiYNGkSjz76KF27dj2twG2Bn38g24ISATi89t8WpxEREWlaDSojFRUVpKamkpyc/NMBbDaSk5NZsWJFnfv94Q9/ICIightuuKFe71NeXk5hYWGtR1tj9hoHQEymVmMVEZHWrUFlJDc3F4/HQ2RkZK3tkZGRZGaeeF2MZcuW8eqrr/LKK6/U+33mzJlDSEhIzSM+Pr4hMVuFnmdfxRGzHRsrY9l9KMfqOCIiIk2mSe+mKSoq4rrrruOVV14hPDy83vvNmjWLgoKCmkdGRkYTpmyZgkM7MCP+He6ovIvF2wusjiMiItJkGrToWXh4OHa7naysrFrbs7KyiIqKOm78999/z969e7n00ktrtnm91XeHOBwOtm/fTrdu3Y7bz+1243Zrsa8L+sWydNcRvticxe2ju1sdR0REpEk06MyIy+ViyJAhpKSk1Gzzer2kpKSQlJR03PjevXuzceNG0tLSah6XXXYZ5513HmlpaW3y65eGuOCMSMCkaP9msrIOWR1HRESkSTR4OfiZM2cyZcoUhg4dyvDhw5k3bx7FxcVMnToVgMmTJxMbG8ucOXPw8/OjX79+tfYPDQ0FOG67HC8y2I83Ql7m7PKlrPwyh8iJD1odSUREpNE1uIxMmDCBnJwcZs+eTWZmJoMGDWLx4sU1F7Wmp6djs2lh18bi7Dwcdiwl7Pv/ACojIiLS+himaZpWh/glhYWFhISEUFBQQHBwsNVxmlXuwX2EvTQQu2FyYPJKYrv2sTqSiIhIvdT381unMFq48JjObPEbBED6N69bG0ZERKQJqIz4gLJeVwAQlf4JPnAiS0REpEFURnxAz/MmUWE66OJNZ8/m+v0ooYiIiK9QGfEBIWHhbAocAcChFYssTiMiItK4Gnw3jVijdPgdXP/FCL7PHc5XXhObzbA6koiISKPQmREfMWTkGFY5hrGvoIp16UesjiMiItJoVEZ8hJ/Tzti+1Uvuf7jhoMVpREREGo/KiA+5om8w9zoW8dv1U6iqrLA6joiISKNQGfEhI3rGMtHxNf3YxZb/fmR1HBERkUahMuJDnC43O8OTAShf/5bFaURERBqHyoiPCRk2EYC++V9ztFAXsoqIiO9TGfExvYYlk2HEEGCUs2WJlocXERHfpzLiYwybjQNdrgIgeKsWQBMREd+nMuKDuo+5iSrTRu+qrezdmmp1HBERkdOiFVh9UHhUJ5YHJbM738ORLUe4s4/ViURERE6dzoz4qNJxz/NQ1Q38bYtJeZXH6jgiIiKnTGXER53bsyORwW7yiiv4vy3ZVscRERE5ZSojPspht3HVmbGcaewg/8t5VscRERE5ZSojPmxiL4P33P/DxCMvcWjfDqvjiIiInBKVER8W16U3m10DsRkme1NesTqOiIjIKVEZ8XGl/ScB0CX9fTxVVRanERERaTiVER/XL/laCggkihw2L/uP1XFEREQaTGXEx/n5B7K148UAeFe9bHEaERGRhlMZaQVix9wJwICSVezftcniNCIiIg2jMtIKxPcYyAb/4RykAymrtDy8iIj4Fi0H30oUX/QcV7y5k4AdLsaXVRLk57Q6koiISL3ozEgrkdS/F106BnG0vIp/p+63Oo6IiEi9qYy0EoZh8LuRXXBSxZ5v/4XXo9+rERER36CvaVqR8YOiOWvxRXQrO8B3S7sx4FcTrI4kIiLyi3RmpBUJ9HORE3lO9ZPVL1kbRkREpJ5URlqZThfOwGMaDChLZd+2dVbHERER+UUqI61MTJfefBd4FgCZS561OI2IiMgvUxlphVwjpwHQP/cz8rIPWJxGRETk5FRGWqEzki5ip6MHAUY52//zlNVxRERETkplpBUybDaKhk0HoOjAVgrLKi1OJCIiUjeVkVZqUPI13Bb4Z24uu4s3VuyzOo6IiEidVEZaKZvdzpjkMQC8umwPpRVaBE1ERFomlZFW7NIBMXRqH4C9OJsvv/iP1XFEREROSGWkFXPYbTw0oJBl7ukMW/v/KC8rsTqSiIjIcVRGWrlzR19AgRFEBHls+HiB1XFERESOozLSyrn9AtjdYyoAsZsXUFVZYXEiERGR2lRG2oABl0/nCMHEmlmkffaq1XFERERqURlpAwLahbA94ToAotf/mYryMosTiYiI/ERlpI0YcNX95BJKrJnF+g/0mzUiItJyqIy0EQHtQvi+z+2Um05Wb99HcXmV1ZFEREQAcFgdQJrP4F/P4Jo9XVibH4jx3z3c8aseVkcSERHRmZG2xOV2c92FIwF4aelujhTrzhoREbGeykgbc+mAGM6IDqZ7xVZWvP201XFERERURtoam83gD8Mred/9COfvfYbM9J1WRxIRkTZOZaQNGpI4mi2u/riNStLfm211HBERaeNURtogw2bDNuZRAIYc+Yzvv1tucSIREWnLVEbaqN5Dzyc16DzshknlR/8P0+u1OpKIiLRRKiNtWNyEuZSYbnpXbiH145esjiMiIm2UykgbFhnXjQ1dbwQgYd2TFBXkWZxIRETaIpWRNu7MCQ/xna0PT1ZO4Lllh6yOIyIibZDKSBvn9gvg8G/+w7uec3ntv/vYlV1kdSQREWljTqmMzJ8/n4SEBPz8/EhMTGT16tV1jn3llVc4++yzCQsLIywsjOTk5JOOl+Z3Xp9IkvtEUOU1eeKDdbqYVUREmlWDy8hbb73FzJkzeeSRR1i3bh0DBw5k7NixZGdnn3D8119/zcSJE/nqq69YsWIF8fHxjBkzhgMHDpx2eGk8D19yBpc5V/PEgSms/+IfVscREZE2xDBN02zIDomJiQwbNoznn38eAK/XS3x8PHfeeScPPPDAL+7v8XgICwvj+eefZ/LkyfV6z8LCQkJCQigoKCA4OLghcaUBVvz1bpL2v0YOYbjuWkNI+45WRxIRER9W38/vBp0ZqaioIDU1leTk5J8OYLORnJzMihUr6nWMkpISKisrad++fZ1jysvLKSwsrPWQpjf4mj+SYcTQkSPs+MedVscREZE2okFlJDc3F4/HQ2RkZK3tkZGRZGZm1usY999/PzExMbUKzc/NmTOHkJCQmkd8fHxDYsop8gtoR/GFz+I1DYblf8Z3X71jdSQREWkDmvVumieffJJFixbx/vvv4+fnV+e4WbNmUVBQUPPIyMhoxpRtW+/EMayO/A0AkUvv19ojIiLS5BpURsLDw7Hb7WRlZdXanpWVRVRU1En3ffrpp3nyySf54osvGDBgwEnHut1ugoODaz2k+QyY/DQHjEgiOczW16dbHUdERFq5BpURl8vFkCFDSElJqdnm9XpJSUkhKSmpzv3+93//lz/+8Y8sXryYoUOHnnpaaRYB7ULIT/4zAPuy81m+I+sX9hARETl1jobuMHPmTKZMmcLQoUMZPnw48+bNo7i4mKlTpwIwefJkYmNjmTNnDgB/+tOfmD17Nm+++SYJCQk115a0a9eOdu3aNeKfIo2p78hxPJvxBn9OsxP3/mYWz+hAO3eD/+ciIiLyixp8zciECRN4+umnmT17NoMGDSItLY3FixfXXNSanp7OoUM/LSv+4osvUlFRwVVXXUV0dHTN4+mnn268v0KaxA1XXExsqD/7j5TyyPvfWR1HRERaqQavM2IFrTNinTV785j20mc87XgR/yETGfbraVZHEhERH9Ek64xI2zMsoT3P9N7GOfaN9F3/KBk7N1gdSUREWhmVEflFSZMeYbNrAAFGOeWLplJeVmJ1JBERaUVURuQX2R0OOk75B0cIorvne9a/ptt9RUSk8aiMSL1ExHZh39nVFx2PyH6btP/7l8WJRESktVAZkXobdP7VrIz4LQBdl83U9SMiItIoVEakQQZf/yzbnGeQ4w3hoQ82UlhWaXUkERHxcSoj0iBuvwA63Pg2t7j/l6WHw5ixKA2Pt8XfHS4iIi2Yyog0WMfIeJ6Zcg5uh40vt2Xz14++sjqSiIj4MJUROSUD4kJ58sr+TLZ/zvXrrmLtJ69aHUlERHyUyoicsisGx3FxfCVOw0Pf1Q+wY91SqyOJiIgPUhmR0zLsxr/wnd9Q/I0KOn54LRm7NlodSUREfIzKiJwWu8NJ19vfZZe9G2EUYl94JbmZ6VbHEhERH6IyIqetXXAYoTf9h/1GFDFmFgWvXE5RQZ7VsURExEeojEijCI+Kx7juPQ4TQjfPbt55bS7lVR6rY4mIiA9QGZFGE9u1L0eueJO53mv4Q1YSMxalUenxWh1LRERaOJURaVTdB45ixOQ/4rLb+WxTJve8uYLKinKrY4mISAumMiKNbmT3cBZcdybB9kp+u+NeNj57lQqJiIjUSWVEmsSvekfytwtdDLNt58zib9j47JUqJCIickIqI9Jkhpx9MVvPfZEK08GZxd+y6dnxVJSXWR1LRERaGJURaVIDf/XbmkIyuHgZ2/58MUcLj1gdS0REWhCVEWlyA3/1W7ad9zIlppsBZakcejaZ3MwMq2OJiEgLoTIizWLA6CvZf/nbHCGYjlWHmPn6V+w7XGx1LBERaQFURqTZ9DxzNEev/ZQH/B7imyMdGP/Ccr7bn291LBERsZjKiDSr+O79+cMd19M3JpjDxRU899ILpH7yitWxRETEQioj0uwigvxYdPMIxnf1Mtf2HEPW3MPKl+7AU1VldTQREbGAyohYIsjPyVPXX8TmmCsBGHHoDTY/fSEFeTkWJxMRkeamMiKWsTscJN3yPGuHPk2p6WJA2RqKnjubPZtXWR1NRESakcqIWG7oJTdxYPwHZNKROPMQMW+PY9Vbf8I0TaujiYhIM1AZkRah+8CRuG5fygb/4biNSrZsTOWmf6RypLjC6mgiItLEDNMH/u9nYWEhISEhFBQUEBwcbHUcaUKm18u37/6F2zYkUOxxEhXsx5+v6kNSzxiro4mISAPV9/NbZ0akRTFsNs757Qzeuv08uoYHkl1YguOfl7Pq+espLsq3Op6IiDQBlRFpkfrFhvDRnaP4fe9shtl2kJj7bwrnDmXjN+9bHU1ERBqZyoi0WIFuBzf+7gY2/up1DhoRRJND/y9/x+p5E3ULsIhIK6IyIi1e/3N+TcjMNazqeBVe02B4/qd4/3Imq999Bo/HY3U8ERE5TSoj4hMCg0JJnPYq2y9+i722eMIoxPndQq6Y/19S9x2xOp6IiJwG3U0jPqeyopzUd59i3vYwVpYnAHDNgGDuHBlBdOde1oYTEZEauptGWi2ny82Iax7i+XtvZMLQeAwD4ra8TIfXzmLlCzdxOGu/1RFFRKQBdGZEfN53GUeoWjiBM8uql5EvNv34rtO19LvqQYJC2lucTkSk7arv57fKiLQOpsmmb/+De+kf6eHZBUAhAWyOu5rel99LWEctmiYi0txURqRNMr1e1n/xDzqsfprO3gwA3vBeyJ5hj3DTOV2IDvG3OKGISNuhMiJtmtfjIe3/3sR/zfPcePQ2DtARh83ghp6lXNEvjF5DzsMwDKtjioi0aiojIoBpmny7M5f5X+1i1Z48XnI+w1j7WnY4elI44Hr6j5mC2y/A6pgiIq2SyojIz2zaf4Tid6cx+MjnuIwqAI4QxPaIi4k89wa69E20OKGISOuiMiJSh7ys/ez4bD5d9y4igrya7Snu88n81Z8Z1z+a0ACXhQlFRFoHlRGRX1BVWcnmb9/Dk/oG/Y8u57Gqa3ndMxan3eDCrn5MidhJn9ETCAwKtTqqiIhPUhkRaYDDWfv5aPMRFn13hG2ZRUywf8WfnK9QYrrZFpSIt9cl9Bh1JSFh4VZHFRHxGSojIqdoZ1YRu5a8Qv/vXyLOzKzZXmna2eY3kJKuY4k9Zypx0ZEWphQRaflURkROk+n1smvDMnLX/pvoQykk/LBuicc0OLP8JTp0jGR0zwjGxpYxsHdP/ALaWZxYRKRlqe/nt6MZM4n4FMNmo8fgc+gx+BwAMnZu4MDKf5N/6HuOVgZRkFPM7pw9XOT6HwxjD5vdvSmMSiK493l0P3O0bhkWEaknnRkROQUFpZUs35XLsm0Hmbb5amLIrvV6melkt6sXWdGj8Yy4iyGdwwgL1B06ItK26GsakWZier3s/34jB9OW4EhfRuei9YSTD8CnnuHcXjkDgK4dAnjc8TJG5BmEdkukU99EAtqFWBdcRKSJqYyIWMT0eknf+R3Zm5eyobAd/zrcje9ziulkZPGN++6acR7TIMMeR25gTyo79sXZ4zw69T2LjkFuLVUvIq2CyohIC3KkuIKtO3dgW/8G7pzviCvZSkeO1Bozv+oynqq6mhB/J0M7VHKL9194w3sTENOH8IR+RMZ1w+7QZV4i4jt0AatICxIW6OKsQf1g0J9qtuUc3MvB7WsoSU/DlbuZjKrB2PKrr0cpP7iR4a6PIO8j2FE9vtx0kmGP5oh/J7bHXklV1/Pp1D6AzmEuYkL8cLrc1vxxIiKnSWdGRFqQskoPu3OKOfT9dwRsexe/IzsIK0snxnOo5vd0AO6tvJl3PKMBOMu2iX8655BrhJHnjOKoXzSVQbHYQmJxtY/Hr9NQwmMS6BDowmbT1z8i0nx0ZkTEB/k57ZwRE8wZMaPg7FE12z1VVRxI38nhfZspydxOJ8cQkks6sO9wCV2P5GAzTCLII6IyDyq3QBFwsHrfe76+hXc95+KwGVwQuIu7zTcodoVT4d8Rb0BHbO0icIZE4hcahX9sX8I6RBLs51RxEZFmozIi4gPsDgexXfsQ27UPACOOec3rGUVuzl3kHdxNUeYeKvP2YuRn4CrNIrA8m0J7HEYJVHlN2pVk0NO5A6p2QAlwuPb7zKi4nQ+8o7DbDMb6beU+/kaJPYRyZzCVrlC8fmHgF4oREMrR2FHYw7sT4u8k2F5FiFFMQHAY/gFBGDZbs82NiPi+Uyoj8+fP56mnniIzM5OBAwfy3HPPMXz48DrHv/POOzz88MPs3buXHj168Kc//YmLL774lEOLyE9sdjvhUZ0Ij+oEjD7u9ZeBSo+XnKJyjmR2Y336QMrzD0HhIYzSXFylufhXHqFdVR5Fzg5QDh6vSWB5JgnODPBmQCXV5eUY09Nu5z/e6rM3F9jW8orrGQCqTBtHjQBK8afMFkC5PZBPQyeyK3QkAS4HnbwHGFawGFwBGK5AbK5AbO5A7O5AHH6BeDr0whEai5/Tjp9Rib9ZgtsvAD//QBxOrdUi0ho1uIy89dZbzJw5kwULFpCYmMi8efMYO3Ys27dvJyIi4rjxy5cvZ+LEicyZM4dLLrmEN998k1//+tesW7eOfv36NcofISIn57TbiAn1Jya0N/TuXee4V4HyKg/5JZUU5PRm08FRlBflUVV8GG9JHkZJHrbyApyVhQSH9+QMTzAFpZW0L6vCYxrYDROH4SWUo4RyFLyAF148cIhPM6p/5yfZlspdrtfrzPBg5fUs9CQD1dfDvOl6oua1StNOBU4qDBcVOPmH62r+z/9CXA4b3cjgpqIX8dhceGxOTJsTr82F1+bEtLvYGXYOGe3Pwmm3EezJo3/ORxh2J9icGPbqB3YnNruTo6G9KGvfG7vNhssspf3h9Rh2Bza7E9sx/2l3OCAwHDMwEofNho0qnKU52G0ODLsdu92BYXdg/+GfbQ4ndocLu83Q7dsix2jwBayJiYkMGzaM559/HgCv10t8fDx33nknDzzwwHHjJ0yYQHFxMR9//HHNthEjRjBo0CAWLFhwwvcoLy+nvLy85nlhYSHx8fG6gFWkBTO9XkqOFlBclE9pYR5lxQVUlBRQVVrIgYA+5DkiOFpeRWDeVnplfoitshh7VSk2TxlOTwlOTxlObxmvOa/m/8zhlFZ6GF61llfsf6rzPWdXTuEfnrEAJBpbecv9xzrH/qnyal70XAZAP2M3H7sfqnPsX6p+zTNVvwWgm3GAFPe9dY79a9VFPFZ1HQAx5LLc7646xy6sOp8Hq24AIMwoYo3rNrwYeLH99DAMTAw+N0bxv/abMAwDPyp4t/IOTGyYgGnY8GJgUv112BrHYF70vwWbYWAY8HzRDAxMTKqPxQ/HNDHY7uzD34Nurhn7YN5DOKkADMDA/KEkmRhkuLrxr9CbMQwDA7g1dw4B5tEfxlLr2DnOGN7pcFv1UQy4JmceQZ786rEGtfbJd4TzfsQ0AAzD4PLsBYRW5VTv+DPF9lA+iLqTH18Zk/t32lccotbIH/YrtwXwn+jpNYc5N+dfdKxI/9kRq1/02Jx8GDuzZmtS7r+JLNv9s6E/vct/Yu+BH+Zh6OGPiCndfuL/koHPYqfjsVWfxRuY9zlxJZuoq3r+X8ytlNsDAeh75EsSjqbVedyvo6+n1BkGQK/8b+hatPaE40wMlkVdx1Fn9a+Mdy9YSffCFXUed2XE1eS7o7lhVBfi2zfuz1g0yQWsFRUVpKamMmvWrJptNpuN5ORkVqw48R+6YsUKZs6cWWvb2LFj+eCDD+p8nzlz5vDoo482JJqIWMyw2QgMDiMwOAxiu9R6bXCtZ92BS+s8zmM/PKqNweu5n4ryEspLSygvK6GitJiqyjIqy0sZ745grLMDFVVezKMJrM0MwVtZhllVjllVAZ4KqKrA9JQTFzSMmwK7UOkxaVdiZ3XmxRhmFTZvFYa36od/rsRmVmEP6cZQvzAqvSYRlUfZXZiADQ9204MdDzbTW/0cD1WOdgQ5HHi9Jn5eqDDt2PFiN47//3neYz6SDLP6LFI1z3FjbVWl5JZWAOBHOZF+x1zg87NDbyyLY2fR0Z9m2O/7Ouf3YLkf647k1zzv6d5EgFF+wrElpaV8m5Nb83yuey3hRuEJx37n7cKSrKya5w+5VhJvyznh2J3eWD7OHF/z/A7XMnrZ9p9w7H4znPcP/TR2smspg2y7Tzg2z2zHtQevrHl+hXMpw+1bTji21HQx8cBvap4nO5eSaF9/wrEAEw78lh+LTJJzKYn2lXWOnXrgMkrwA2CA42tGOL6pc+y0/ReQS/VKzI86viHRsaTOsf/vwNnsN6u/gXjA8V8SHR/VOXb2/qHsMOMBmG7/L9c6/13n2Cf39yfNLOeyQTGNXkbqq0FnRg4ePEhsbCzLly8nKSmpZvt9993H0qVLWbVq1XH7uFwuXn/9dSZOnFiz7YUXXuDRRx8l65j/4R5LZ0ZExNeZXi9erxePpwpPVSVeTxUew47X7ofHa+LxVGEU5+D1evB6qvB6qr/T8no81fs5AqgIjMI0weupwp27Ga/XC5h4vR7wejFNL6bXS6U7lOKQHmCCxzQJPfgtpmliml445j8xTcrcYeS3H1z9OhCxfwmGWVk9DrNmXPXY9mR3PKv6KRC3/xNs3up/NxumSfXHxw9jnWFkRP2KHz9REg5+jLPyx4L0wzEB0zQpdwazO3rcD8+h26GPcVfmU7Phx7ZlmlTaA9kaO77mpR6ZH+NfcRjjx+MdM+cem4vv4n76rOmR+SntyjP5eS808WIadtbF/65mW/eszwkpyzhmUO2dVsffiPlDl+yem0KHkhMXIkyT1XFT8dqcAHQ7/BUdi3fUGmIcc+w1sddRaa8uAF3yviX66JYfMh7/0Zwacw3ljurPwM5HVhBbmHbiDMD6qN9S4uoAQHzBGjoVHHsWpfaxN0SOp8gdyTWJnYkN9a/zmKfCp2/tdbvduN1awElEfJdhs2G32aqvK3H7nWCEG0IC63/AuLPrP7bH+F8e86P+k+s/dsgt9R/LtJO+em6tZ3eedGztv3z6SceOqvWs7q/MAM6q9azHSccm1XrWs9HGjqj1rFcDjtvQsdfWc6w1GnT/XXh4OHa7/bgzGllZWURFRZ1wn6ioqAaNFxERkbalQWXE5XIxZMgQUlJSarZ5vV5SUlJqfW1zrKSkpFrjAZYsWVLneBEREWlbGvw1zcyZM5kyZQpDhw5l+PDhzJs3j+LiYqZOnQrA5MmTiY2NZc6cOQBMnz6dc889l7lz5zJu3DgWLVrE2rVrefnllxv3LxERERGf1OAyMmHCBHJycpg9ezaZmZkMGjSIxYsXExkZCUB6ejq2Y1ZfPOuss3jzzTd56KGH+P3vf0+PHj344IMPtMaIiIiIAPqhPBEREWki9f381g9IiIiIiKVURkRERMRSKiMiIiJiKZURERERsZTKiIiIiFhKZUREREQspTIiIiIillIZEREREUu1yF/t/bkf12UrLCy0OImIiIjU14+f27+0vqpPlJGioiIA4uPjLU4iIiIiDVVUVERISEidr/vEcvBer5eDBw8SFBSEYRiNdtzCwkLi4+PJyMjQMvNNTHPdfDTXzUvz3Xw0182nsebaNE2KioqIiYmp9bt1P+cTZ0ZsNhtxcXFNdvzg4GD9D7uZaK6bj+a6eWm+m4/muvk0xlyf7IzIj3QBq4iIiFhKZUREREQs1abLiNvt5pFHHsHtdlsdpdXTXDcfzXXz0nw3H81182nuufaJC1hFRESk9WrTZ0ZERETEeiojIiIiYimVEREREbGUyoiIiIhYSmVERERELNWmy8j8+fNJSEjAz8+PxMREVq9ebXUknzdnzhyGDRtGUFAQERER/PrXv2b79u21xpSVlTFt2jQ6dOhAu3btuPLKK8nKyrIocevw5JNPYhgGM2bMqNmmeW5cBw4c4Nprr6VDhw74+/vTv39/1q5dW/O6aZrMnj2b6Oho/P39SU5OZufOnRYm9k0ej4eHH36YLl264O/vT7du3fjjH/9Y64fWNNen5ptvvuHSSy8lJiYGwzD44IMPar1en3nNy8tj0qRJBAcHExoayg033MDRo0dPP5zZRi1atMh0uVzma6+9Zm7evNm86aabzNDQUDMrK8vqaD5t7Nix5t/+9jdz06ZNZlpamnnxxRebnTp1Mo8ePVoz5tZbbzXj4+PNlJQUc+3ateaIESPMs846y8LUvm316tVmQkKCOWDAAHP69Ok12zXPjScvL8/s3Lmz+bvf/c5ctWqVuXv3bvPzzz83d+3aVTPmySefNENCQswPPvjA3LBhg3nZZZeZXbp0MUtLSy1M7nsef/xxs0OHDubHH39s7tmzx3znnXfMdu3amc8++2zNGM31qfn000/NBx980HzvvfdMwHz//fdrvV6feb3wwgvNgQMHmitXrjS//fZbs3v37ubEiRNPO1ubLSPDhw83p02bVvPc4/GYMTEx5pw5cyxM1fpkZ2ebgLl06VLTNE0zPz/fdDqd5jvvvFMzZuvWrSZgrlixwqqYPquoqMjs0aOHuWTJEvPcc8+tKSOa58Z1//33m6NGjarzda/Xa0ZFRZlPPfVUzbb8/HzT7Xab//rXv5ojYqsxbtw48/rrr6+1bfz48eakSZNM09RcN5afl5H6zOuWLVtMwFyzZk3NmM8++8w0DMM8cODAaeVpk1/TVFRUkJqaSnJycs02m81GcnIyK1assDBZ61NQUABA+/btAUhNTaWysrLW3Pfu3ZtOnTpp7k/BtGnTGDduXK35BM1zY/vwww8ZOnQov/nNb4iIiGDw4MG88sorNa/v2bOHzMzMWvMdEhJCYmKi5ruBzjrrLFJSUtixYwcAGzZsYNmyZVx00UWA5rqp1GdeV6xYQWhoKEOHDq0Zk5ycjM1mY9WqVaf1/j7xq72NLTc3F4/HQ2RkZK3tkZGRbNu2zaJUrY/X62XGjBmMHDmSfv36AZCZmYnL5SI0NLTW2MjISDIzMy1I6bsWLVrEunXrWLNmzXGvaZ4b1+7du3nxxReZOXMmv//971mzZg133XUXLpeLKVOm1Mzpif6dovlumAceeIDCwkJ69+6N3W7H4/Hw+OOPM2nSJADNdROpz7xmZmYSERFR63WHw0H79u1Pe+7bZBmR5jFt2jQ2bdrEsmXLrI7S6mRkZDB9+nSWLFmCn5+f1XFaPa/Xy9ChQ3niiScAGDx4MJs2bWLBggVMmTLF4nSty9tvv83ChQt588036du3L2lpacyYMYOYmBjNdSvWJr+mCQ8Px263H3dnQVZWFlFRURalal3uuOMOPv74Y7766ivi4uJqtkdFRVFRUUF+fn6t8Zr7hklNTSU7O5szzzwTh8OBw+Fg6dKl/OUvf8HhcBAZGal5bkTR0dGcccYZtbb16dOH9PR0gJo51b9TTt+9997LAw88wNVXX03//v257rrruPvuu5kzZw6guW4q9ZnXqKgosrOza71eVVVFXl7eac99mywjLpeLIUOGkJKSUrPN6/WSkpJCUlKShcl8n2ma3HHHHbz//vt8+eWXdOnSpdbrQ4YMwel01pr77du3k56errlvgPPPP5+NGzeSlpZW8xg6dCiTJk2q+WfNc+MZOXLkcbeo79ixg86dOwPQpUsXoqKias13YWEhq1at0nw3UElJCTZb7Y8mu92O1+sFNNdNpT7zmpSURH5+PqmpqTVjvvzyS7xeL4mJiacX4LQuf/VhixYtMt1ut/n3v//d3LJli3nzzTeboaGhZmZmptXRfNptt91mhoSEmF9//bV56NChmkdJSUnNmFtvvdXs1KmT+eWXX5pr1641k5KSzKSkJAtTtw7H3k1jmprnxrR69WrT4XCYjz/+uLlz505z4cKFZkBAgPnPf/6zZsyTTz5phoaGmv/5z3/M7777zrz88st1u+kpmDJlihkbG1tza+97771nhoeHm/fdd1/NGM31qSkqKjLXr19vrl+/3gTMZ555xly/fr25b98+0zTrN68XXnihOXjwYHPVqlXmsmXLzB49eujW3tP13HPPmZ06dTJdLpc5fPhwc+XKlVZH8nnACR9/+9vfasaUlpaat99+uxkWFmYGBASYV1xxhXno0CHrQrcSPy8jmufG9dFHH5n9+vUz3W632bt3b/Pll1+u9brX6zUffvhhMzIy0nS73eb5559vbt++3aK0vquwsNCcPn262alTJ9PPz8/s2rWr+eCDD5rl5eU1YzTXp+arr7464b+fp0yZYppm/eb18OHD5sSJE8127dqZwcHB5tSpU82ioqLTzmaY5jHL2omIiIg0szZ5zYiIiIi0HCojIiIiYimVEREREbGUyoiIiIhYSmVERERELKUyIiIiIpZSGRERERFLqYyIiIiIpVRGRERExFIqIyIiImIplRERERGx1P8HeC0HzS9JZOAAAAAASUVORK5CYII=", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "fs = 1000\n", + "phi = 0.9\n", + "tau = -1/(np.log(phi) * fs)\n", + "lags = np.arange(0, 100)\n", + "\n", + "plt.plot(lags, np.exp(-lags / (tau * fs)), label=r\"$e^{\\frac{-t}{\\tau f_s}}$\")\n", + "plt.plot(lags, phi**lags, ls='--', label=r\"$\\varphi^t$\")\n", + "plt.legend()\n", + "plt.title(\"ACF\");" + ] + }, + { + "cell_type": "markdown", + "id": "71500018-27e6-40fc-94b5-67ac22d1ff63", + "metadata": {}, + "source": [ + "Take the discrete FT of the ACF is:\n", + "\n", + "$$\n", + "\\begin{align*}\n", + "FFT(f) &= \\int_{-\\infty}^{\\infty} A(t) e^{-i 2 \\pi f t} \\, dt \\\\\n", + "&= \\int_{-\\infty}^{\\infty} \\varphi^t e^{-i 2 \\pi f t} \\, dt \\\\\n", + "&= \\sum_{t=-\\infty}^{\\infty} \\varphi^t e^{-i 2 \\pi f \\frac{1}{f_s} t} \\quad\\quad\\quad \\text{(discrete representation)} \\\\\n", + "&= \\sum_{t=-\\infty}^{\\infty} \\left(\\varphi e^{-i 2 \\pi f \\frac{1}{f_s}}\\right)^t\n", + "\\end{align*}\n", + "$$\n", + "\n", + "The sum of this geometric series is:\n", + "\n", + "$$\n", + "\\begin{align*}\n", + "FFT(f) &= \\frac{1}{1 - \\varphi e^{-i 2 \\pi f \\frac{1}{f_s}}}\n", + "\\end{align*}\n", + "$$\n", + "\n", + "The power spectral density is the absolute value squared of the Fourier transform of the ACF.\n", + "\n", + "$$\n", + "\\begin{align*}\n", + "P(f) &= \\frac{1}{|1-\\varphi e^{-i2\\pi f \\frac{1}{f_s}}|^2}\n", + "\\end{align*}\n", + "$$\n", + "\n", + "This is the form of an AR(1) PSD." + ] + }, + { + "cell_type": "markdown", + "id": "1124b743-47ca-49d9-b91c-9a9134ed9c88", + "metadata": {}, + "source": [ + "## Derivation of a Lorentizan PSD from an AR(1) PSD" + ] + }, + { + "cell_type": "markdown", + "id": "0c74dd7a-ae20-4ebf-a606-732a0508ee73", + "metadata": {}, + "source": [ + "#### Euler's formula to remove complex exponential\n", + "\n", + "$$\n", + "\\begin{align}\n", + "P(f) &= \\frac{1}{|1-\\varphi e^{-i2\\pi f \\frac{1}{f_s}}|^2} \\\\\n", + "&= \\frac{1}{(1-\\varphi e^{-i 2 \\pi f \\frac{1}{f_s}})(1-\\varphi e^{i 2 \\pi f \\frac{1}{f_s}})} \\\\\n", + "&= \\frac{1}{1-\\varphi e^{-i 2 \\pi f \\frac{1}{f_s}} - \\varphi e^{i 2 \\pi f \\frac{1}{f_s}} + \\varphi^2} \\\\\n", + "&= \\frac{1}{1 - \\varphi (e^{-i 2 \\pi f \\frac{1}{f_s}} + e^{i 2 \\pi f \\frac{1}{f_s}}) + \\varphi^2} \\\\\n", + "&= \\frac{1}{1 - 2 \\varphi \\cos{(2 \\pi f \\frac{1}{f_s})} + \\varphi^2} \\quad\\quad\\quad \\text{(Euler's formula)}\n", + "\\end{align}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "id": "1877810e-4837-49fe-b20f-50a3fab5e9b9", + "metadata": {}, + "source": [ + "#### 1st order Taylor approximation of the cosine term\n", + "This gives bias in the high frequencies. This could be alleviate by lowering the upper frequency bound during fitting.\n", + "\n", + "$$\n", + "\\begin{align}\n", + "\\cos(2 \\pi f \\frac{1}{f_s}) &\\approxeq 1 - \\frac{(2 \\pi f \\frac{1}{f_s})^2}{2} \\\\\n", + "P(f) & \\approxeq \\frac{1}{1 - 2\\varphi (1 - \\frac{(2 \\pi f \\frac{1}{f_s})^2}{2}) + \\varphi^2}\n", + "\\end{align}\n", + "$$" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "a23a13cb-5690-47dc-a510-76bd79186373", + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkIAAAHLCAYAAAAk8PeNAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8g+/7EAAAACXBIWXMAAA9hAAAPYQGoP6dpAABu9klEQVR4nO3dd3gUVd/G8e/upneSQAq9996LFEEBFQEFEZUuNhQVK/ooKvooKgpKAAUECyo2sCtFAohILyI9BukdEtKT3fP+kZd9jAkhCQmbcn+uay+yM2dmfrs7ZO/MnJljMcYYRERERMogq6sLEBEREXEVBSEREREpsxSEREREpMxSEBIREZEyS0FIREREyiwFIRERESmzFIRERESkzFIQEhERkTJLQUhERETKLAUhEZECiI6OxmKxEB0d7epSROQyKAiJyBU1b948LBaL8+Hl5UWdOnW4//77OX78eJa2+/fvZ8SIEdSsWRMvLy/Cw8Pp3LkzEyZMyNKua9euzvVZrVYCAgKoW7cuQ4YMYcmSJXmubfjw4Vlqc3Nzo3Llytx6663s2LGjUF6/iBQvbq4uQETKphdeeIHq1auTkpLCr7/+yowZM/jhhx/Yvn07Pj4+7Nu3j9atW+Pt7c3IkSOpVq0aR48eZdOmTUyaNInnn38+y/oqVarEyy+/DEBiYiL79u3jq6++4qOPPuKWW27ho48+wt3d/ZJ1eXp6Mnv2bAAyMjKIiYlh5syZ/PTTT+zYsYPIyEgAOnfuTHJyMh4eHoX8zojIlaQgJCIu0bt3b1q1agXAnXfeSUhICG+88QZff/01gwcP5s033yQhIYEtW7ZQtWrVLMueOHEi2/oCAwO54447skx75ZVXGDt2LNOnT6datWpMmjTpknW5ubllW0+7du244YYb+P777xk9ejQAVqsVLy+vfL1mESl+dGpMRIqFq6++GoDY2FgAYmJiqFSpUrYQBFChQoU8rdNms/HWW2/RoEEDpk2bRlxcXIFqCw8PBzJD0gU59RFatWoVAwcOpEqVKnh6elK5cmUefvhhkpOTs6zv2LFjjBgxgkqVKuHp6UlERAR9+/Zl//79BapPRApOQUhEioWYmBgAQkJCAKhatSoHDx7kl19+uaz12mw2Bg8eTFJSEr/++mueljl16hSnTp3i+PHjrFmzhocffpiQkBBuuOGGXJf7/PPPSUpK4t577+Xtt9+mZ8+evP322wwdOjRLu5tvvpmFCxcyYsQIpk+fztixYzl//jwHDhwo8OsUkYLRqTERcYm4uDhOnTpFSkoKq1ev5oUXXsDb29sZNsaOHcuHH35I9+7dadasGV26dKFbt25cc801+Pj45GtbjRo1Av4XtnKTmJhI+fLls0yrWLEiixcvzjb93yZNmoS3t7fz+V133UWtWrV46qmnOHDgAFWqVOHcuXP89ttvvPbaazz66KPOtuPHj8/PSxKRQqIjQiLiEj169KB8+fLOq7L8/PxYuHAhFStWBKBhw4Zs2bKFO+64g/379zN16lT69etHWFgYs2bNyte2/Pz8ADh//vwl23p5ebFkyRKWLFnCzz//zDvvvIOfnx/XXXcde/bsyXXZf4agxMRETp06RYcOHTDGsHnzZmcbDw8PoqOjOXv2bL5eh4gUPh0REhGXiIqKok6dOri5uREWFkbdunWxWrP+bVanTh0+/PBD7HY7O3bs4LvvvuPVV1/lrrvuonr16vTo0SNP20pISADA39//km1tNlu29V533XXUrl2b8ePH8+WXX1502QMHDvDss8/yzTffZAs5F/oneXp6MmnSJB555BHCwsKcHbGHDh3q7IskIleOgpCIuESbNm2cV41dis1mo3HjxjRu3Jj27dvTrVs35s+fn+cgtH37dgBq1apVoForVapE3bp1Wbly5UXb2O12rrnmGs6cOcMTTzxBvXr18PX15fDhwwwfPhyHw+Fs+9BDD9GnTx8WLVrEzz//zDPPPMPLL7/ML7/8QvPmzQtUo4gUjE6NiUiJciE8HT16NE/t7XY7H3/8MT4+PnTq1KnA283IyHAeWcrJH3/8wZ49e5g8eTJPPPEEffv2pUePHs77Dv1bzZo1eeSRR1i8eDHbt28nLS2NyZMnF7g+ESkYBSERKZZWrVpFenp6tuk//PADAHXr1r3kOux2O2PHjmXnzp2MHTuWgICAAtWyZ88edu/eTdOmTS/axmazAWCMcU4zxjB16tQs7ZKSkkhJSckyrWbNmvj7+5Oamlqg+kSk4HRqTESKpUmTJrFx40ZuuukmmjRpAsCmTZv44IMPCA4O5qGHHsrSPi4ujo8++gjIDBsX7iwdExPDrbfeysSJE/O03YyMDOd6HA4H+/fvZ+bMmTgcjmxDe/xTvXr1qFmzJo8++iiHDx8mICCAL7/8MltfoT179tC9e3duueUWGjRogJubGwsXLuT48ePceuuteX17RKSQKAiJSLH01FNP8fHHH7NixQrmz59PUlISERER3HrrrTzzzDNUr149S/tDhw4xZMgQIPMqsYiICNq3b8+MGTO45ppr8rzd1NRU53oAAgICaN26tfNS/otxd3fn22+/ZezYsbz88st4eXnRv39/7r///ixHkipXrszgwYNZtmwZH374IW5ubtSrV4/PPvuMm2++Oc91ikjhsJh/HscVERERKUPUR0hERETKLAUhERERKbMUhERERKTMUhASERGRMktBSERERMosBSEREREps3QfoVw4HA6OHDmCv78/FovF1eWIiIhIHhhjOH/+PJGRkdkGc/43BaFcHDlyhMqVK7u6DBERESmAgwcPUqlSpVzbKAjlwt/fH8h8Iws6RpGIiIhcWfHx8VSuXNn5PZ4bBaFcXDgdFhAQoCAkIiJSwuSlW4s6S4uIiEiZpSAkIiIiZZaCkIiIiJRZCkIiIiJSZikIiYiISJmlICQiIiJlloKQiIiIlFkKQiIiIlJmlfogdPDgQbp27UqDBg1o0qQJn3/+uatLEhERkWKi1N9Z2s3NjSlTptCsWTOOHTtGy5Ytue666/D19XV1aSIiIuJipT4IRUREEBERAUB4eDihoaGcOXNGQUhERESK/6mxlStX0qdPHyIjI7FYLCxatChbm6ioKKpVq4aXlxdt27Zl3bp1Oa5r48aN2O12jSgvIiIiQAk4IpSYmEjTpk0ZOXIkN910U7b5CxYsYNy4ccycOZO2bdsyZcoUevbsye7du6lQoYKz3ZkzZxg6dCizZs266LZSU1NJTU11Po+Pjy/cF/P/0jIcjHp/PY+efDpzgsWKsVgxWMBiAawc9KrD4pA7sFosWCxw67HXcDMZmW2xYLFYMBYrWKyc9qzM72GDsfx/225HZuPhSP7/9doy12mxYrFYSfQoz58VB2CzWnCzWqh37Bs8HUlYrG5gdcNqtYHNDYvVht0zkDOR3bD+f9typzfh4UgFmw2rzR2rze3/H+7YPH0htDbuNgseNivuGYl4uFlx9/DC3cMDi9VWJO+liIjI5bAYY4yri8gri8XCwoUL6devn3Na27Ztad26NdOmTQPA4XBQuXJlHnjgAZ588kkgM+Bcc801jB49miFDhlx0/c899xzPP/98tulxcXGFOvp8cpqd+s/+xH6v2y7aZrm9KSPSn3A+3+E5Ah9Lao5t1zrqMSjtWefzDZ73EGrJOcRtc1TnxrSXnM9XeTxIZevJHNvuc0TSI+115/OfPR6nrvVQjm2PmGA6pE5zPl/o8SzNrfucz9ONjXTcSLe4cZZAbvOclhmUbFYeTplOTcd+MqyeZNi8sNs8sdu8cNi8cLj5sLLGw3i52fByt1Lz3BoC7Kexenjj5umHm48/nj6BePgE4O0XiFdwRXw93bFZLz3isIiIlE7x8fEEBgbm6fu72B8Ryk1aWhobN25k/PjxzmlWq5UePXqwZs0aAIwxDB8+nKuvvjrXEAQwfvx4xo0b53weHx9fJKfR3G0W3hzUlA0H/osxBmMc4HCAcWT+bBw4PMN4LrQBDgMOY9h6YCwWkw6OzDYWYzDGDsYQ7xHG2Aq1nG13HroVN3sSFmPg/9eHcWAxDs66V2Bw+SrYHQ7sDvjraGdOZJzB4rBjMZkPq7FjMRmctpanTWQwGQ4HdgOn4qrgnuGGBQc2MtvZjB0bds5YgvH3dCPN7iDN7sCdjKyv2WLHHTuQSppx40hcinNeqEcsda27c3yvEowXgw/c6Hw+z/1d2ti2XvS9rZHyEQ6seLvbeNFtNm35g1SrN6k2P1Ld/MnwCMTuEYDxCmJP7VH4+/oS6O1O+fQjBLjbCQiJIDA4DJtbif6vISIieVSif9ufOnUKu91OWFhYlulhYWHs2rULgNWrV7NgwQKaNGni7F/04Ycf0rhx42zr8/T0xNPTs8jrdrNZ6d+8EjQfk4+l/pPr3J5Znr2Wa9teWZ69l2vbrlme/XjRdhWAP/7x3J7Rk5S0FNLSUslISyEjPY2MtFTS01LIsNv5JqAm6XYHaRkGt+MvsjnhBPa0JExaMo70ZEx6MiY9hXS7YURYNVLSHaSm20k83pStyV642VNwdyTj+f8PH5KxGAeO/+/2lpxupxwnqWQ7BnYyH2lA0v9qvGNPJ+xknrKb6j6NvrbfMms3Fk5bAoi3BpHoXo5Uj2B+rjGegMByhPh5UiXjb0I80vAPjiQ4rBLevv65vociIlJ8legglBedOnXC4XC4uowyx+bmhs3NDy8fv0s3rnltrrM7Z3n29kXbGYeDXXZDYmoGial20k5WZlfcMdKS4klPOkdG0jkcSWchJQ6TlkSv4ErEp6QTl5yO2zlvzmX4EUQCNoshhDhCHHGQ+jekwi1rR2An8xTiVPdpdPz/0ARwlgBO28pz3jOMNN9I/qz/EOVDQogM8qaSj53QckE6wiQiUkyV6N/OoaGh2Gw2jh8/nmX68ePHCQ8Pd1FV4ioWqxUvK3i52wjxA0KaAE0u2r5DlmedAEhPSyXuzHHiTx0l8cxRUuOOk5pwlrvK1eF0QiqnE9JwOxbMkdQKlHOcw9uSRjniKWePh6QY7IkWbj94IxkcAOBN9yhusP7OCUswpz0rkuhXDRNSC+/wuoRUbUh41bq4KSSJiLhMif4N7OHhQcuWLVm2bJmzA7XD4WDZsmXcf//9BV5vVFQUUVFR2O32QqpUSgp3D09Cw6sQGl4ly/SOWZ7NBzKPQMXFneH04Rjij8eSeupvUhPOcL1vFY6cS+bIuRQiks7gbrETwUkiUk9C6hY4DewBh7HQOGMeYcFB1Aj143rrGir5Ogiq0YLKdVvg5a17XYmIFLVif9VYQkIC+/ZlXn3UvHlz3njjDbp160ZwcDBVqlRhwYIFDBs2jHfeeYc2bdowZcoUPvvsM3bt2pWt71B+5afXuUhO7BkZnDr2N2cPxxB/ZDf2k3vxjIulXPIBMuwZXJP6qrPtZx7P0+b/O41nGCuHbBU55VuH9PIN8anclPAW11MhwAuLRVfEiYjkJj/f38U+CEVHR9OtW7ds04cNG8a8efMAmDZtGq+99hrHjh2jWbNmvPXWW7Rt2/ayt60gJEXJYXdw9HwqsScTiT2VQOWtUwk9u4mKqTGU43yWtkdNMO1TpxHq50GzyuUY4LORyhUrU71JR3z8Al30CkREiqdSFYRcSUFIXME4HJw8+jdHd68n6eAW3E/t5HCKJw8nDMFhAAzrPe+jvCWODGNlv1t1Tgc1xlq5LZVbXEN4ldqufgkiIi6lIFRIFISkOElJt/PnkTj+iD1Kkw3jqZTwBxU4k63dSls7fmz4Gu1rhtK+Rgjl/Yv+lhAiIsWJglAhURCS4u74oRgO/7GStP1rCT69kZrpe5llv4FXMgYD4EMKX/u8yOnQ1vg2uo46bXvi6ent4qpFRIqWgtBl+udVY3v27FEQkhLj/LnTbNl/ghWHDKtjTlPh+Cre95jknJ9gvNnt1xp7rZ7U7NCPkLBKLqxWRKRoKAgVEh0RkpLu7OmT/LX2Oxx7llDj3GpCOOec5zAW3goYh2/rO7iuSQQVg3SkSERKBwWhQqIgJKWJw24nZtuvnN70DaFHllPLHsPVqa/zl4kEYGj4fnqXP0ONrrcTVrGGi6sVESk4BaFCoiAkpdnJI3/z034H3/1xjHX7zzDD7U162dYD8Kd7YxLrD6LBNUPw8w9ybaEiIvmkIFRIFISkrDgRn0LMj29Tbt9X1Evf4ZyeaDz5M6gbPm2H0qBtb6w2qwurFBHJGwWhQqIgJGXR8UP7iF32HhX3L6SyOQLAVkcN7vebzG1tqnJLq0qE+OmSfBEpvhSEComCkJRlxuFg94ZfiP99Ht+ersSHqVcBUM6WyozyXxB01WjqNu+CxaqjRCJSvCgIXSZdPi+SVXKane+2HeHD3/+mydEveNF9LgB7bLWJb3Evza4dgpu7h4urFBHJpCBUSHRESCS73VtWc375VJqcW4aHJQOAw5YwDtQdSdM+9+Hjq/8rIuJaCkKFREFI5OLOnjzC7m/fpO6BT5yDxJ4mkE/aLmLwVQ3Vj0hEXEZBqJAoCIlcWnLiebZ9N53Ku+awNaMK96Y/jJe7laHtq3F3h0hCggJdXaKIlDH5+f5WL0cRuSzevv60HfQEYU//iVvft2hSKZCUdAc/rFqL9c2G/DbrYeLOnHR1mSIiOdIRoVzoiJBI/hljiN59kmNfP8vg5E8AiDc+bK86hEY3P0lAYLCLKxSR0k6nxgqJgpBIwRmHna1LPiJo7etUcxwA4AwB7GnwAC37P4S7rjITkSKiIHSZdPm8SOFx2O1s+Wku5TdMdt6gcZe1Fodu/o7uDcKxWCwurlBEShsFoUKiI0IihSc9LZXNX02m9q7pfJBxDW9mDKB9jRD+c0N9GkaqQ7WIFB4FoUKiICRS+OLPnWLWrwd55/djpGU4aGHdy+MVt1P/9kkEBoW4ujwRKQUUhAqJgpBI0Tl4JolXf9zB3btH0ci6n1ME8Vfz8bTuc5eG7RCRy6LL50Wk2Ksc7MPbt7fC0uM5DlkiCeUcbTY/wZ+vdOGvnRtdXZ6IlBEKQiLiUg0796fCE5tYW/0+Uow7jdK2UenTa/n1vfGkpaa6ujwRKeUUhETE5Ty8vGk77GXiRq5mm3dbPCwZdDownTemvsrWg+dcXZ6IlGIKQiJSbIRVrUvjx35iU8tJLKUd75xpRv/pq3n5h52kpNtdXZ6IlEIKQiJSrFisVlr0uYfmj35Dn6aVcBj4aOWfbH2lBzs3rnB1eSJSyigI5SAqKooGDRrQunVrV5ciUmaF+Hny1uDmzB7aisd9vqWtfRO1vunPr++NJz093dXliUgpocvnc6HL50WKh/gzJ4iZexfNzy8HYIdbAwJue49KNeq7uDIRKY50+byIlCoBwRVoPm4hm1q+TALeNMjYQdD73fj9q2kYh8PV5YlICaYgJCIlg8VCiz73kTAiml0eDfCzJNNu29Msinqc+BSdKhORglEQEpESJbxqPeo8vpL1Ne7jmAnmv4eb0eftX9l+OM7VpYlICaQgJCIljtXNndZDX+bEsNV4BEXy9+kkbpr+Gz/99A3q9igi+aEgJCIlVpMakfww9ip61A/jWrOaXr8PYe2bt5KQEO/q0kSkhFAQEpESLdDHnVlDW3JHA3fsxkK7+J848UYn/tq11dWliUgJoCAkIiWexWKh3R0T2Nf7Y04RRA3H34R+0osNSz51dWkiUswpCIlIqVG33XXY7v2V3R4NCbAk0eLXe1j13lM47LrEXkRypiAkIqVKubDK1Hz0FzaW74fVYrjqQBSvzp7HeV1iLyI5UBASkVLHzcOLlmPeZ1PjCUQ5+jMzNox+UauJOZng6tJEpJhREMqBxhoTKR1a3DyOq+6aQkSgFzEnExkV9T3bfl/q6rJEpBjRWGO50FhjIqXDyfOpjP1wNU8ce5T6lgNsbPkyHW4c7eqyRKSIaKwxEZF/KO/vydxhrbAGhONpSafDpkdZMedJdaIWEQUhESkbvHwDaPTQN2yMuBWALgdn8PvU20hJSXZxZSLiSgpCIlJmWN3caHn3O2xu9DR2Y6FD/I/smdyT06eOu7o0EXERBSERKXOaD3icPd1nk4gXTdK3snPGHRw4neTqskTEBRSERKRMqt95AGdu+Zrdlho8lTSYm2b8xp9HNIK9SFmjICQiZVblBu0o99Bv+IbX5lRCKoPe+Z21O/9ydVkicgUpCIlImVYh0JsFd7ejbfVgWqevp/6nHVm39HNXlyUiV4iCkIiUeQFe7rw/sg0PBv1GgCWJ5qvuZtVX011dlohcAQpCIiKAl7uNxg8vZGtQd9wtdq7aNp7lH0xE95wVKd0UhERE/p/N3ZMmYz9nc8QgALr99TrRc55SGBIpxRSERET+wWK10fyud9hc/S4Auh2aTvTMh3QXapFSSkFIROTfLBaaD3uNLXUfBuDY4b957IttZCgMiZQ6bq4uoDiKiooiKioKu93u6lJExIWaDX6O1T834NkVPqRvPkyK3cGUQc1wt+lvSJHSQqPP50Kjz4sIwE/bj/LAJ5tx2DN4vuI6Btz1DF6enq4uS0QuQqPPi4gUol6NInh3aCte85jFHaffZsuUgSQnp7i6LBEpBApCIiJ50K1uBRp0u5V0Y6Nd8gq2vTWQlNRUV5clIpdJQUhEJI/qdbud2B7vkm5stE1eyeYpCkMiJZ2CkIhIPtS5agCxPWaSbmy0T17B5qm3KAyJlGAKQiIi+VTnqluIvXpGZhhKimbtW0NISddVpiIlkYKQiEgB1OkyiNhuUSQYb2afbc5dH25UGBIpgRSEREQKqE7Xwey6dTUbbM1ZueekwpBICaQgJCJyGVrVr8ncEa3xdrdxcO82lk27j/SMDFeXJSJ5pCAkInKZ2tUIYd7tDfnU40Wuj/uUVdPuwq7hOERKBAUhEZFC0LZeZc60exKAq899yS8zHsTh0I37RYo7BSERkUJSv/c9/Nn8WQCuOfUBS2c/hUYxEineFIRERApRw76PsL1B5qj11x6ZzuL3X1IYEinGFIRERApZo1ueY3vN0QD03P8aP34+y8UVicjFKAiJiBSBRne8xvZKg9nsqMWTmwKZ82usq0sSkRwoCImIFAWLhUYjp/Nbx3nE48fE73bw2fqDrq5KRP5FQUhEpKhYrdx3bSPu7lwDgG2L3uD3lT+7uCgR+ScFIRGRImSxWHiydz1errWDF93fo+6yUWzdvM7VZYnI/1MQEhEpYhaLhYG3302sR13KWc5TftFg9uzZ7eqyRAQFoRxFRUXRoEEDWrdu7epSRKSUcPMOIOK+bzliq0ik5RS2j2/m0OHDri5LpMyzGN3g4qLi4+MJDAwkLi6OgIAAV5cjIqVA/LEY0t7pQag5wx/WeoTd/yMVgoNdXZZIqZKf728dERIRuYICwmtiueNLzuNLY8cuYmfcQnxSiqvLEimzFIRERK6wkJotSLjpI1Lw4JekWoz+YBMp6XZXlyVSJikIiYi4QESTq9k/eBXz3fqzdv9ZHvx0M3YN0ipyxSkIiYi4SL269Zg1tBUeblZ+/XM/Cz56V+OSiVxhCkIiIi7UvmYIb/WvyQKPiQyKeYLFX811dUkiZYqCkIiIi/VqURu3Si2wWQydtz3BymjdfVrkSlEQEhFxNYuFuqPeZW9AO7wtaTRYPpotWze7uiqRMkFBSESkGLDY3Klx7+cc8KhFqCWOwIW3EfP3AVeXJVLqKQiJiBQTNu8AKtzzNSet5anOERLm3cLxM+dcXZZIqaYgJCJSjHgFV8Jj2Fck4EuE4yj/mfcTCakZri5LpNRSEBIRKWYCqzYh4eaPGen2MktO+HPf/E2k2x2uLkukVFIQEhEphsIbd+Wl4dfh7W5j5Z6TvPTZKt1jSKQIKAiJiBRTTSsHMe225vSyreexXQP5/st5ri5JpNRREBIRKca61w9jXPUD+FpS6frHk6xcFe3qkkRKFQUhEZFirs7wmcT6t8DPkkKtpSP5Y9duV5ckUmooCImIFHduHlS55yuOuVUi0nIay6e3cejEaVdXJVIqKAiJiJQANt9y+I/6iniLP43Yx75ZQ4lPTnV1WSIlnoKQiEgJ4RtRl7Sb3ycdN7qm/8rHs18nQ5fVi1wWBSERkRIktFF3jnd+hbmO63n1cBOe+/ZPXVYvchkUhERESphKV48mctAbGIuVj34/wNzV+11dkkiJpSAkIlIC9WwYzvje9XAnA356ktXr1ru6JJESSUFIRKSEGn1VDT6ouJCRbj8R/v0wdu4/6OqSREocBSERkRLKYrHQ6o4XOWMNoablMHHv38GJuARXlyVSoigIiYiUYO7lKuE+5DOS8aSd2cLvM8eQkm53dVkiJYaCkIhICedfvRXxvaYBcGPyIha+N0lXkonkkYKQiEgpENbuFg40HgvATUcms+jbhS6uSKRkUBASESklqvR/nr8rdCcVD778fQ9Ldxx3dUkixV6ZCEL9+/enXLlyDBgwwNWliIgUHauVqqM+YHb9OfzqaMyDn25m97Hzrq5KpFgrE0HowQcf5IMPPnB1GSIiRc/TjwcG9qJdjWAS0+w8Me8nziRoTDKRiykTQahr1674+/u7ugwRkSvC3WZlxu0t6Re4l3nJD7J45mOkZWhMMpGcFPsgtHLlSvr06UNkZCQWi4VFixZlaxMVFUW1atXw8vKibdu2rFu37soXKiJSjJTz9WB8Oy+CLIncmvA+n344Q1eSieSg2AehxMREmjZtSlRUVI7zFyxYwLhx45gwYQKbNm2iadOm9OzZkxMnTuR7W6mpqcTHx2d5iIiUVGHd7uFQ7SEA3LT/Bb7+eYmLKxIpfop9EOrduzcvvvgi/fv3z3H+G2+8wejRoxkxYgQNGjRg5syZ+Pj48N577+V7Wy+//DKBgYHOR+XKlS+3fBERl6p065scLtcaP0sKLdbcx2/bdru6JJFipdgHodykpaWxceNGevTo4ZxmtVrp0aMHa9asyff6xo8fT1xcnPNx8KDG7RGREs7mTuSdCzjlHkkVy0ncvhpBzLGzrq5KpNgo0UHo1KlT2O12wsLCskwPCwvj2LFjzuc9evRg4MCB/PDDD1SqVOmiIcnT05OAgIAsDxGRks7iG0LAyM9JtnjThj9Z+t4E4pLTXV2WSLHg5uoCroSlS5e6ugQREZfyiGhE3I3vsOTb93kjvhu/f7qZ2cNaY7NaXF2aiEuV6CNCoaGh2Gw2jh/PevfU48ePEx4e7qKqRESKp8Dmfak+ci64ebF890neWKL+QiIlOgh5eHjQsmVLli1b5pzmcDhYtmwZ7du3L/B6o6KiaNCgAa1bty6MMkVEio3GlQKZdHMTrDgwK99gybo/XF2SiEsV+1NjCQkJ7Nu3z/k8NjaWLVu2EBwcTJUqVRg3bhzDhg2jVatWtGnThilTppCYmMiIESMKvM0xY8YwZswY4uPjCQwMLIyXISJSbPRrXpFKa56h1YkvWf/9NnZGfk/9SiGuLkvEJSymmN9hKzo6mm7dumWbPmzYMObNmwfAtGnTeO211zh27BjNmjXjrbfeom3btpe97QtBKC4uTh2nRaRUyTi+i/SZ3fA2SXzhdj3dH55HOV8PV5clUijy8/1d7IOQKykIiUhplrD1G/wWZt5w8Z1yjzLq/qdxs5XoHhMiQP6+v7XHi4iUUX5Nb+RUy4cBGH5mKu9/sdDFFYlceQpCIiJlWOj1z3I8ohuelnR673iU73/f6uqSRK4oBaEc6KoxESkzrFbChr3PGa8qBHOer77/kT8Oxbm6KpErRn2EcqE+QiJSVjiO72Liok3MjQ0kItCLb+7vRHl/T1eXJVIg6iMkIiL5Yg2rx8PDBlKjvC9H41J44KP1pGU4XF2WSJFTEBIREQACvNyZNbQV7T338/LRkcz57CtXlyRS5BSERETEqWZ5P6ZWXkF163Fu3P04C1dtcXVJIkVKQUhERLKocPu7nPOuQkXLacKX3MvG2BOuLkmkyCgIiYhIVt5BBI74jBSLF+2tO9j+4WOciE9xdVUiRUJBKAe6fF5EyjpLhfrQNwqAYY5FzJvztjpPS6mky+dzocvnRaSsi1v0OIFb3uG88WZug/cYO+g6V5ckckm6fF5ERApFYJ+XOFu+NWscDZi1OZEvNx5ydUkihUpBSERELs7mTrlRX/HnVdM5jw9PLfyD7Yd152kpPRSEREQkd14BPNijLlfXq0Bqhp3X3v+Cs4lprq5KpFAoCImIyCVZrRbeHNCAWb7vMDv1Maa9/yF2h7qYSsmnICQiInkS6OtDm+rlcLfYufv4C8z4brWrSxK5bApCOdDl8yIiObBYCLxlBvEBtalgOUfrDY/w09YDrq5K5LLo8vlc6PJ5EZEcnNpHyozOeNkTed9cT4f73qF2mL+rqxJx0uXzIiJSdEJr4X7TTACGWb7nk/emEJ+S7uKiRApGQUhERPLN1vBGktqMBeChlOn855PVONR5WkogBSERESkQn54TOFvjRu6zP8I3uxOZHr3P1SWJ5JuCkIiIFIzNjXJDP6RP31sAmLxkD9G7NVK9lCwKQiIiclkGta7C4DZVqMFhln3yJgdOJ7m6JJE8UxASEZHL9txVvnzn9SzPmRlMmzuX5DS7q0sSyRMFIRERuWyeodUx9W7AZjE8dn4S//1sObo7i5QECkI50A0VRUTyyWLBp/9bJAbVpbwljhv2PM37v6rztBR/uqFiLnRDRRGRfDq1j7QZnfGwJzLTfiPNR0yhbY0QV1clZYxuqCgiIq4RWgv3/lEA3GP7hgUfvcOxuBQXFyVycfkOQunp6bi5ubF9+/aiqEdEREo4S6P+ZLS+G4Dr0hdz7/yNpGU4XFyVSM7yHYTc3d2pUqUKdruuCBARkZy59XyRM51f4HHbY2w+cI6J3+1wdUkiOSrQqbGnn36ap556ijNnzhR2PSIiUhq4eRB89YO8fmsrAD78/W++2HjIxUWJZOdWkIWmTZvGvn37iIyMpGrVqvj6+maZv2nTpkIpTkRESrar64Ux7urqeK+cyNpFq6kX/gSNKga6uiwRpwIFoX79+hVyGSIiUlrdH7wOq9sPJBlP7v6gDm+NvY1yvh6uLksE0OXzudLl8yIihcBhJ/2Dm3DfH02MI4JXKs9g5qiu2KwWV1cmpdQVuXz+3LlzzJ49m/Hjxzv7Cm3atInDhw8XdJUiIlIaWW24D5xDum84Na1H6XNgEm8u3u3qqkSAAgahbdu2UadOHSZNmsTrr7/OuXPnAPjqq68YP358YdYnIiKlgW8o7oPex2Fx40bbGuJWzWDxn8dcXZVIwYLQuHHjGD58OHv37sXLy8s5/brrrmPlypWFVpyIiJQiVdphvfYFAJ5x+5C5n33JXycTXFyUlHUFCkLr16/n7rvvzja9YsWKHDtW8hO+xhoTESki7e7DUa8PGVZPfNJPc/eHG0lMzXB1VVKGFSgIeXp6Eh8fn236nj17KF++/GUX5Wpjxoxhx44drF+/3tWliIiULhYL1n5RJA9fyh++Hdh7IoHHv9ymkerFZQoUhG688UZeeOEF0tPTAbBYLBw4cIAnnniCm2++uVALFBGRUsYrkJCqDZlxRwvcrBaWbvub2atiXV2VlFEFCkKTJ08mISGBChUqkJycTJcuXahVqxb+/v689NJLhV2jiIiUQi2rBhPVIYGVng+x8ucvWBNz2tUlSRlUoBsqBgYGsmTJEn799Ve2bdtGQkICLVq0oEePHoVdn4iIlGLX2ldhsZzjDbe3GTK/CnMf7EtEoLery5IypEA3VExJSclytVhppRsqiogUsfRkHLO6Yz3xJ+scdXmlwmt8ck8nPN1srq5MSrAiv6FiUFAQnTt35plnnuGXX34hOTm5QIWKiEgZ5+6NddCHODz8aWPdzbXHZvH8txqpXq6cAgWhpUuX0qtXL9auXcuNN95IuXLl6NSpE08//TRLliwp7BpFRKQ0C6mJtV8UAPe4fcux9Yv4bP1BFxclZcVljzWWkZHB+vXreeedd5g/fz4OhwO73V5Y9bmUTo2JiFxBPzwO697hrPGjn/0V3r6nD00qBbm6KimB8vP9XaDO0pB5z6Do6GjnIzU1lRtuuIGuXbsWdJUiIlKWXTsRc2g9u+PKcfK0D/d+tIlvH+hEsEaqlyJUoCNCFStWJDk5ma5du9K1a1e6dOlCkyZNsFhK10jCOiIkInKFpcQT5/Cib9Rq9p9OolOtUN4f2UYj1Uu+FHln6fLly5OUlMSxY8c4duwYx48fV4dpERG5fF4BBPp48M6QVni7WzkYs53XNVK9FKECBaEtW7Zw7NgxnnzySVJTU3nqqacIDQ2lQ4cOPP3004Vdo4iIlDF1gwy/VHyXbzz+w7cr1vDT9pI/jqUUT5fdWfr06dNER0fz9ddf88knn6iztIiIXL6MNJjbGw5vYIujBsMtL/LFmC7UquDn6sqkBCjyU2NfffUVY8eOpUmTJoSFhXHvvfeSkJDA5MmT2bRpU4GKLk40+ryIiIu5ecDAuRivIJpZ/2Ks/QPu/nADCRqpXgpZgY4IVahQgc6dOzs7Sjdu3LgoanM5HRESEXGx3T/BJ4MAuDvtIawNbmT67S1K3cU5UriK/PL5EydOFKgwERGRfKnbCzqMhd/e4jX3d7n+z6q8uzKIu7vUdHVlUkoU+D5CdrudRYsWsXPnTgAaNGhA3759sdk0PoyIiBSi7s/CwbUEHFzLVPcoBvxUgUYVA+lYK9TVlUkpUKA+Qvv27aN+/foMHTqUr776iq+++oohQ4bQsGFDYmJiCrtGEREpy2zuMOA9TERTfq35KA5j4YFPNnP4nG7bIpevQH2ErrvuOowxzJ8/n+DgYCDz6rE77rgDq9XK999/X+iFuoL6CImIFCPGkJLhYMDM39h+OJ6mlQJZcHd7vNx1JkKyys/3d4GCkK+vL7///nu2TtJbt26lY8eOJCQk5HeVxZKCkIhI8XPwTBLj3p7PiRQr7Vu15pWbm7i6JClmivzyeU9PT86fP59tekJCAh4eGhNGRESKTuWTK/nU+jTT3aeycH0Mn6474OqSpAQrUBC64YYbuOuuu1i7di3GGIwx/P7779xzzz3ceOONhV2jiIjI/0Q0webpT0Pr3zzj9iHPfv0nWw6ec3VVUkIVKAi99dZb1KpViw4dOuDl5YWXlxcdO3akVq1aTJ06tbBrFBER+Z+ASLjpXQwW7nBbRi/zK/d9tJHTCamurkxKoHz1EXI4HLz22mt88803pKWlUaVKFYYNG4bFYqF+/frUqlWrKGu94tRHSESkGFs2EVa9ThJe3JD6IuE1GvPByDa42Qr0N76UIkXWR+ill17iqaeews/Pj4oVK/LDDz+waNEi+vTpU+pCkIiIFHNdx0PVTviQwnSPt9gYc5TXNFK95FO+gtAHH3zA9OnT+fnnn1m0aBHffvst8+fPx+FwFFV9IiIiObO5wYA54FueepYD3GFbyjsr/uKHP466ujIpQfIVhA4cOMB1113nfN6jRw8sFgtHjhwp9MJEREQuyT8cbp4NXcfj1uFeAB77fCv7TmS/slkkJ/kKQhkZGXh5eWWZ5u7uTnp6eqEWJSIikmc1ukLXJ3msVwPa1wghMc3OXR9u5HyKvpvk0vI11pgxhuHDh+Pp6emclpKSwj333IOvr69z2ldffVV4FYqIiOSBm83K27fUZ+HbjzP5ZE8e/XwrM+9oqZHqJVf5CkLDhg3LNu2OO+4otGJEREQuR+h3oxhtX0I5j2M8+uddzFgRw31ddTGPXFy+gtDcuXOLqg4REZHL13EsxCxjgDWaX60NeP1naFwxkKtql3d1ZVJM6WYLIiJSelTvDF2eAGCS11yqc5ixn2zm0NkkFxcmxZWCkIiIlC6dH4NqV+HpSGa2z3SSkhK596NNpKTbXV2ZFEMKQjmIioqiQYMGtG7d2tWliIhIflltmZfU+5anuj2Wl7zm88fhOJ5ZtJ18DKYgZUS+htgoazTEhohICRbzC3x4E+keAXQ6/1+Om3K81L8Rt7et6urKpIgV2RAbIiIiJUbNq+HGt3C/71eG92wPwHPf/MnmA2ddXJgUJwpCIiJSerUYCkFVuKdLDXo1DCfdbrj3o02c0kj18v8UhEREpNSzWCy80fwYEwMWcSw+hfs/3kSGXeNkioKQiIiUBaf24vPF7QxJ+4y+Hhv4/a8zTPppl6urkmJAQUhEREq/0NrQ4QEAXvecRSXLCWatiuW7bRo0vKxTEBIRkbKh+7NQqTXu6ef5LPhd3Mng8S+2see4RqovyxSERESkbLC5w4D3wCuIyMQdTAn9mqQ0O3d/uJF4jVRfZikIiYhI2RFUBfpNB+D6hC8Z4Led2FOJPPLZVhwO3VavLFIQEhGRsqXe9dD2XgAer38aD5uVJTuOM2NFjIsLE1dQEBIRkbLnmhfgts+ocPOrvNC3IQCvL97Nij0nXVyYXGkKQiIiUva4eUCdngDc2qYKg1tXwhh48NPNHDyjkerLEgUhEREp2xJO8FLiBIZX2Mu5pHTu+WijRqovQxSERESkbFsThfWv5TyTPpU6Pgn8eSSepxdqpPqyQkFIRETKtq7jIawxtuQzfB46B3eLnS83HeKjtQdcXZlcAQpCIiJStrl7wcB54OFH4Im1fFp3JQAvfPsnG//WSPWlnYKQiIhIaC24YQoALfbP5uGaR0i3G+6bv5ET51NcW5sUKQUhERERgCYDofkQLBgeOPcqrULTOR6fyv0fbyZdI9WXWgpCIiIiF/R+FcrXx+rpzxt9quLn6ca62DO88qNGqi+tFIREREQu8PCB2z6Fu1dQpW5zXh/YFIA5v8by9ZbDLi5OioKCkIiIyD+Vqwae/gD0ahTOA50rA/Dkl3+w61i8CwuToqAgJCIikhNjYM10xu25g9413ElOt3PPhxuJS9ZI9aWJgpCIiEhO0pNh41ws5/5mqudMKgV6sv90EuMWbNFI9aWIgpCIiEhOPHwy7y/k5oVH7DI+b7oRDzcry3ad4M2le1xdnRQSBSEREZGLCWsIvScBELF+Eu92zRyD7O1f9vHt1iOurEwKiYKQiIhIbloMg0Y3g7HT9Y8nGdshBIDHvtjK9sNxLi5OLpeCkIiISG4slsy7TgfXgLiDPJwwhS61Q0lJdzD6gw2cPJ/q6grlMigIiYiIXIpXQGZ/IXdfLFU78NbgZtQo78vRuBTu+WgjqRl2V1coBVTqg9B3331H3bp1qV27NrNnz3Z1OSIiUlJFNIWHt0PHsQT6eDJ7aCv8vdzY+PdZnlm0HWN0JVlJVKqDUEZGBuPGjeOXX35h8+bNvPbaa5w+fdrVZYmISEnlE+z8sUaghRkD62C1wGcbDjHvt/2uq0sKrFQHoXXr1tGwYUMqVqyIn58fvXv3ZvHixa4uS0RESrrjO+DdbnTaOZGnetcDYOJ3O1i196SLC5P8KtZBaOXKlfTp04fIyEgsFguLFi3K1iYqKopq1arh5eVF27ZtWbdunXPekSNHqFixovN5xYoVOXxYY8WIiMhlSkuEMzGw/UtG+f7KzS0q4TBw/8ebiT2V6OrqJB+KdRBKTEykadOmREVF5Th/wYIFjBs3jgkTJrBp0yaaNm1Kz549OXHiRIG2l5qaSnx8fJaHiIhINpVbw9XPAGD58Qn+28lG8ypBxCWnM/qDDcSnaBiOkqJYB6HevXvz4osv0r9//xznv/HGG4wePZoRI0bQoEEDZs6ciY+PD++99x4AkZGRWY4AHT58mMjIyItu7+WXXyYwMND5qFy5cuG+IBERKT06jIWa3SEjGc+Fo3hnUH3CA7zYdyKBhz7dgl3DcJQIxToI5SYtLY2NGzfSo0cP5zSr1UqPHj1Ys2YNAG3atGH79u0cPnyYhIQEfvzxR3r27HnRdY4fP564uDjn4+DBg0X+OkREpISyWqH/O+AXBid3UWH1BN4d2hJPNyu/7DrBaz/vdnWFkgclNgidOnUKu91OWFhYlulhYWEcO3YMADc3NyZPnky3bt1o1qwZjzzyCCEhIRddp6enJwEBAVkeIiIiF+VXHm56F7DApg9ocu4XXh3QBICZK2JYuPmQa+uTS3JzdQFF7cYbb+TGG290dRkiIlJa1egKnR+FA79Dlfb0DYhg97HzTI+O4Ykv/qBKsA8tqwZfcjXiGiX2iFBoaCg2m43jx49nmX78+HHCw8NdVJWIiJRJXZ6EoV9DQAQAj15bl2sbhJFmd3DXBxs5eCbJxQXKxZTYIOTh4UHLli1ZtmyZc5rD4WDZsmW0b9/+stYdFRVFgwYNaN269eWWKSIiZYHNDaw251PriT+ZcmszGlUM4HRiGiPnrdeVZMVUsQ5CCQkJbNmyhS1btgAQGxvLli1bOHDgAADjxo1j1qxZvP/+++zcuZN7772XxMRERowYcVnbHTNmDDt27GD9+vWX+xJERKQsMQa+GwczO+ITu5TZQ1sTFuDJ3hMJjJm/iQy7w9UVyr8U6yC0YcMGmjdvTvPmzYHM4NO8eXOeffZZAAYNGsTrr7/Os88+S7NmzdiyZQs//fRTtg7UIiIiV4TFAjaPzJ8X3UM4p5kzrDXe7jZW7T3F89/u0JhkxYzF6BO5qPj4eAIDA4mLi9MVZCIikjcZqTDnGji6Fap0gGHf8vOuU9zz0UaMgQl9GjCiY3VXV1mq5ef7u1gfERIRESlx3DxhwFzw8IcDv8GKSfRsGM6Tvf43JtnyXQUbAUEKn4KQiIhIYQupCX2mZP688jX4awV3da7BoFaVcRh44JPN7DqmYZyKAwWhHOiqMRERuWyNB0DzIYCBr+7Ckp7ExH6NaF8jhITUDEbN28CJ8ymurrLMUx+hXKiPkIiIXJa0JPjoZuhwP9S7HoBzSWncNP03/jqVSLPKQXx6Vzu83G2XWJHkh/oIiYiIFAcePjDiB2cIAgjy8WDO8NYE+biz5eA5xn22BYcGaHUZBSEREZGiZLH87+e4w3B0G9VDfZl5R0s8bFZ++OMYL/2w03X1lXEKQiIiIlfC4Y0wsxN8ejskn6VdjRBeG5g5QOucX2N579dYFxdYNikIiYiIXAkhtcErAOIOwDcPgDH0bVaRx3vVBWDi9zv48Y+jLi6y7FEQEhERuRK8AjLvL2R1h53fwvrZANzbpSZ3tKuCMfDQgi1s/PuMiwstWxSEcqDL50VEpEhUbAHXPJ/5889Pw9FtWCwWnuvTkO71KpCa4eDO9zfw18kE19ZZhujy+Vzo8nkRESl0xsAnt8KenyCkFty1Ajz9SErLYPC7v7P1UBxVgn346r4OhPp5urraEkmXz4uIiBRXFgv0nQ7+kXB6H/z6BgA+Hm7MHtaaysHeHDiTxKh560lKy3BxsaWfgpCIiMiV5hsCA+ZAy+Fw1aPOyeX9PZk3og1BPu5sPRTH2E82k2F3uK7OMkBBSERExBWqdoA+UzNvuvgPNcv7MXtoKzzcrCzdeYKnF25HvViKjoKQiIiIqzkcsGEupGeOPdaqWjBv3docqwUWbDjIqz/vdnGBpZeCkIiIiKstvAu+ewiWPOOc1KtROP/t3xiAGdExzFr5l4uKK90UhHKgy+dFROSKanJr5r/r3s28x9D/u7VNFecNF1/6YSdfbDzkiupKNV0+nwtdPi8iIlfM4mfgt7fAKxDu+RWCqgBgjOGl73cy+9dYbFYL79zRkh4NwlxcbPGmy+dFRERKmu7PQsVWkBIHX4wCezoAFouFp66rz00tKmJ3GMZ8vIm1f512cbGlh4KQiIhIcWBzhwHvgWcgHFoHy19yzrJaLUy6uUmWu0/vOBLvwmJLDwUhERGR4qJcVej7dubPv8+A+P8NwupusxJ1ewtaVyvH+dQMhr63TkNxFAIFIRERkeKkQV+4+hkYtRgCIrLM8nK3MXtYa+pHBHAqIZXbZ6/l4JkkFxVaOigIiYiIFDedH4WIpjnOCvR258NRbahVwY+jcSkMnvU7R84lX+ECSw8FIRERkeLs8CbY8F6WSaF+nnx8Z1uqhfhw6Gwyt89ey4n4FBcVWLIpCImIiBRXJ/fAnGvh+0fg79+yzKoQ4MX80e2oGORN7KlEbp+9ltMJqS4qtORSEMqBbqgoIiLFQvk60OhmMA748k5IOpNldsUgbz4Z3Y7wAC/2nkhgyJx1xCWlu6jYkkk3VMyFbqgoIiIul5oA73aB0/ugTm8Y/AlYLFmaxJxMYNA7v3MqIZWmlQL5YFRbAr3dXVSw6+mGiiIiIqWFpx8MmAs2D9jzI6ydma1JzfJ+zL+zLeV83Nl6KI7bZ//OuaQ0FxRb8igIiYiIFHcRTeDa/7/B4uJn4MjmbE3qhvvz8eh2hPh6sP1wPINnqc9QXigIiYiIlARtRkO9G8CRDps+yLFJ/YgAPr2rHaF+nuw8Gs/gWb9z8rzCUG4UhEREREoCiwX6ToPrXofrJl+0We0wfxbc3Y6wAE/2HE/g1nfX6NL6XCgIiYiIlBTe5TKPDFlz//quWd6PBXe1JzLQi5iTiQx6VzddvBgFIRERkZIoNQG+fRBO7MxxdrVQXxbc3d55n6EBM35j3wmNTfZvCkIiIiIl0eL/wMZ58PkISMt5vLHKwT58dk97apT35UhcCgNn/saWg+euaJnFnYKQiIhISdTtKfCtACd3ws/jL9qsYpA3n9/dniaVAjmblM5ts35n1d6TV7DQ4k1BSEREpCTyqwA3vQtYMo8Mbf/qok1D/Dz5eHQ7OtUKJSnNzsh56/lu25ErVmpxpiAkIiJSUtXsBp0ezvz52wfh7P6LNvXzdGPO8FZc3ziCdLvhgU82M+fXWMr6ABMKQjnQWGMiIlJidHsKKrWB1Hj4YiTYLz7WmKebjbcGN2dIu6oYAxO/28EzX28nw+64ggUXLxprLBcaa0xEREqEs3/DO1eB1Q2G/wAV6uXa3BjD7FWx/PfHnRgDV9UOJer2FgR4lY7xyfLz/a0glAsFIRERKTFiV0FILQiIyPMiP/95jIc+3UJyup06YX7MGdaaysE+RVjklaFBV0VERMqa6lflKwQB9GwYzuf3tHfehbpv1GpW7ztVRAUWTwpCIiIipc2Or+GzYeC4dN+fRhUDWTSmI40qBnAmMY0hc9YyPXpfmelErSAkIiJSmiScgIX3wI5FsHpKnhaJCPTmi3s6MLBlJRwGXv1pN3d/uJH4lIt3vC4tFIRERERKE78KcN1rmT//8iIcWJunxbzcbbw6oAkv39QYD5uVxTuO0+ftX0v9nagVhEREREqbZrdD44Fg7PDlKEg+m6fFLBYLg9tU4fN7Msco+/t0EjfP+I1pv+zF7iidp8oUhEREREobiwWufwPKVYe4g/DNA5CPPj9NKwfxw9iruKFJBHaH4fXFe7j13TUcPJPzmGYlmYKQiIhIaeQVAAPeA6s77PwWNszJ1+KBPu68Pbg5b9zSFD9PN9bvP8u1b65k1sq/StUNGBWERERESquKLeCa5zN/Pvt3vhe3WCzc1KISPz54FW2rB5OcbuelH3bSN2o1fxyKK+RiXUM3VMyFbqgoIiIlnjFwcC1UaXeZqzF8vuEQL/2wk7jkdKwWGNCyEuOuqUt4oFchFVs4dGfpQqIgJCIipY7DAdaCnxA6eT6VF7/fwddbMkev93K3MvqqGtx5VQ0CvYvHEB26s7SIiIhkF3cY3r8B/viiwKso7+/J1Fub8+W9HWhZtRwp6Q7e/mUfnV75hVd+3MWJ8ymFWHDR0xGhHERFRREVFYXdbmfPnj06IiQiIqXDytcy7y3k4Q93r4CQmpe1OmMMP/95jMmL97D3RAIAHjYrvRqFM7BVJTrWDMVqtRRG5fmiU2OFRKfGRESkVLFnZB4ROrAGIprBqCXg5nHZq3U4DL/sOsH06H1sOnDOOT0i0Ise9cO4ul4F2tYIxsfD7bK3lRcKQoVEQUhEREqduEMws1PmTRbb3w89Xyq0VRtj+ONwHJ9vOMSiLYc5n5LhnGe1QK0KfjSMDCQyyIuwAC98PdzwcrdxfZP8DRZ7KQpChURBSERESqVdP8CngzN/vu0zqNOz0DeRkm5n9b5T/LLrBNG7T3L4XHKO7bzcreya2LtQt52f7+8rc4xKREREio9610Hbe2DtTFh0L9zzKwREFuomvNxtdK8fRvf6YQCciE/hj8Nx7Dp2nuPxKRyPTyE53YGbC/oQ/ZOCkIiISFl0zQvw929gT4e0oh86o0KAF90DvJzBqLhQEBIRESmL3Dxh8KfgXQ48fFxdjcsoCImIiJRVgRWzPk9PBndv19TiIrqhooiISFnncMCqyRDVBpLOuLqaK0pBSEREpKzLSIYtH8O5A/D1mMzxycoIBSEREZGyzsMXBswFmwfs/gHWvevqiq4Y9REqBHa7nfT0dFeXISJS5Nzd3bHZbK4uQ4pCRBO49iX48TFY/J/M0eojmrq6qiKnIHQZjDEcO3aMc+fOuboUEZErJigoiPDwcCwW197/RYpAm9HwVzTs/h4+H5E5Hpmnv6urKlIKQpfhQgiqUKECPj4++qUgIqWaMYakpCROnDgBQERE4Q6LIMWAxQJ9p8HMrXAmBr5/FG56x9VVFSkFoQKy2+3OEBQSEuLqckRErghv78xLq0+cOEGFChV0mqw08gmGm2fDRzdBxZaZHadL8R/6CkIFdKFPkI9P2b0JlYiUTRd+76WnpysIlVZV28ND28G39P+hr6vGLpNOh4lIWaPfe2XEP0NQSjykp7iuliKkICQiIiIXd3gTvNMZljzr6kqKhIKQiIiIXFziKTgbC+vegV3fu7qaQqcgJJctLS2NWrVq8dtvv+V5mR07dlCpUiUSExOLsLIrr1q1akyZMiXXNs899xzNmjW7IvXkxf79+7FYLGzZssXVpYhIcVTnWmh/f+bPX4+BuEOuraeQKQjlICoqigYNGtC6dWtXl1Kk1qxZg81m4/rrr88278KX44VHcHAwXbp0YdWqVdnazpw5k+rVq9OhQwfntJdeeokOHTrg4+NDUFBQtmUaNGhAu3bteOONN3Ktcd68ec4arFYrlSpVYsSIEc7LdwFWrFjB1VdfTXBwMD4+PtSuXZthw4aRlpYGQHR0dJZ1BAYG0rx5cx5//HGOHj2a17crT9avX89dd93lfG6xWFi0aFGhbqOwVa5cmaNHj9KoUaMi3c7w4cPp169fkW5DRIpI9wkQ2RySz8KXo8Ge4eqKCo2CUA7GjBnDjh07WL9+vatLKVJz5szhgQceYOXKlRw5ciTHNkuXLuXo0aOsXLmSyMhIbrjhBo4fP+6cb4xh2rRpjBo1KstyaWlpDBw4kHvvvfei2x8xYgQzZswgIyP3/1ABAQEcPXqUQ4cOMWvWLH788UeGDBkCZB5Z6tWrF61atWLlypX88ccfvP3223h4eGC327OsZ/fu3Rw5coT169fzxBNPsHTpUho1asQff/yR6/bzo3z58iXuSkKbzUZ4eDhubiXjItILAVdEriA3DxjwHnj4w4HfYOWrrq6o8Bi5qLi4OAOYuLi4bPOSk5PNjh07THJysnOaw+EwianpLnk4HI58vbbz588bPz8/s2vXLjNo0CDz0ksvZZkfGxtrALN582bntG3bthnAfP31185p69evN1ar1cTHx+e4nblz55rAwMAc56WmphpPT0+zdOnSi9aZ0/IvvfSSsVqtJikpybz55pumWrVqub7W5cuXG8CcPXs2y/SkpCRTt25d07Fjx4su27JlS/Paa685n/ft29e4ubmZ8+fPG2OMOXjwoAHM3r17jTHGVK1a1bz55pvOnwHno2rVqsYYYyZMmGCaNm1qPvjgA1O1alUTEBBgBg0adNH38IJff/3VdOnSxXh7e5ugoCBz7bXXmjNnzhhjjElJSTEPPPCAKV++vPH09DQdO3Y069atcy575swZc9ttt5nQ0FDj5eVlatWqZd577z1jTPbP+sL7tXTpUtOyZUvj7e1t2rdvb3bt2pWlnkWLFpnmzZsbT09PU716dfPcc8+Z9PT0HGufMGFClvcCMMuXLzfGGHPgwAEzcOBAExgYaMqVK2duvPFGExsb61x22LBhpm/fvubFF180ERERplq1as6aFyxYYDp16mS8vLxMq1atzO7du826detMy5Ytja+vr+nVq5c5ceJEru+r5F9Ov/+kjNj2uTETAoyZEGjMXytdXc1F5fb9/W8l40/AEiI53U6DZ392ybZ3vNATH4+8f5yfffYZ9erVo27dutxxxx089NBDjB8//qKXxSYnJ/PBBx8A4OHh4Zy+atUq6tSpg79//m/B7uHhQbNmzVi1ahXdu3fP83Le3t44HA4yMjIIDw93HrHq3Llzvrbv7e3NPffcw8MPP+y8Ody/denShejoaB599FGMMaxatYqgoCB+/fVXevXqxYoVK6hYsSK1atXKtuz69eupUKECc+fOpVevXlnutxITE8OiRYv47rvvOHv2LLfccguvvPIKL730Uo61btmyhe7duzNy5EimTp2Km5sby5cvdx71evzxx/nyyy95//33qVq1Kq+++io9e/Zk3759BAcH88wzz7Bjxw5+/PFHQkND2bdvH8nJybm+P08//TSTJ0+mfPny3HPPPYwcOZLVq1cDmZ/70KFDeeutt7jqqquIiYlxnhKcMGFCtnU9+uij7Ny5k/j4eObOnQtAcHAw6enp9OzZk/bt27Nq1Src3Nx48cUX6dWrF9u2bXPua8uWLSMgIIAlS5ZkWe+ECROYMmUKVapUYeTIkdx22234+/szdepUfHx8uOWWW3j22WeZMWNGrq9VRPKo8YDMITjOH4Py9VxdTaFQECqj5syZwx133AFAr169iIuLY8WKFXTt2jVLuw4dOmC1WklKSsIYQ8uWLbOElr///pvIyMgC1xEZGcnff/+d5/Z79+5l5syZtGrVCn9/fwYOHMjPP/9Mly5dCA8Pp127dnTv3p2hQ4cSEBBwyfXVq5f5H3n//v05BqGuXbsyZ84c7HY727dvx8PDg0GDBhEdHU2vXr2Ijo6mS5cuOa67fPnywP/GZfonh8PBvHnznAFyyJAhLFu27KJB6NVXX6VVq1ZMnz7dOa1hw4YAJCYmMmPGDObNm0fv3r0BmDVrFkuWLGHOnDk89thjHDhwgObNm9OqVSsgs1P3pbz00kvO1/bkk09y/fXXk5KSgpeXF88//zxPPvkkw4YNA6BGjRpMnDiRxx9/PMcg5Ofnh7e3N6mpqVnei48++giHw8Hs2bOdIXzu3LkEBQURHR3NtddeC4Cvry+zZ892BqP9+/cDmQGrZ8+eADz44IMMHjyYZcuW0bFjRwBGjRrFvHnzLvlaRSQfrp8MVnewlo7eNQpChcjb3caOF3q6bNt5tXv3btatW8fChQsBcHNzY9CgQcyZMydbEFqwYAH16tVj+/btPP7448ybNw93d3fn/OTkZLy8vApet7c3SUlJubaJi4vDz88Ph8NBSkoKnTp1Yvbs2UBm/5a5c+fy4osv8ssvv7B27Vr++9//MmnSJNatW3fJsZCMMcDFbxB31VVXcf78eTZv3sxvv/1Gly5d6Nq1K6+88gqQ2VH7sccey+/Lplq1almOokVERGTpAP5vW7ZsYeDAgTnOi4mJIT093fnlD5kjhLdp04adO3cCcO+993LzzTezadMmrr32Wvr165elc3tOmjRpkqU+yBxWoUqVKmzdupXVq1dnCW52u52UlBSSkpLy3E9q69at7Nu3L9sRxZSUFGJiYpzPGzdunOVIZE41hoWFOdv+c1pu76uIFICbZ9bnp/ZCaG3X1FIIFIQKkcViydfpKVeZM2cOGRkZWY7kGGPw9PRk2rRpBAYGOqdXrlyZ2rVrU7t2bTIyMujfvz/bt2/H0zPzP0JoaOhldTY+c+YMNWvWzLWNv78/mzZtwmq1EhER4Rzr6J8qVqzIkCFDGDJkCBMnTqROnTrMnDmT559/Ptd1XwgKFztCEhQURNOmTYmOjmbNmjVcc801dO7cmUGDBrFnzx727t170SNCuflnmITMfcfhcFy0fU6vOT969+7N33//zQ8//MCSJUvo3r07Y8aM4fXXX89TjReC4oUaExISeP7557npppuyLZefYJyQkEDLli2ZP39+tnkXjqhB5hGhvNb472m5va8ichky0uC7h2DbZzDq58xxyUqg0nFcS/IsIyODDz74gMmTJ7NlyxbnY+vWrURGRvLJJ59cdNkBAwbg5uaW5fRM8+bN2bVrl/PISn5t376d5s2b59rGarVSq1YtatSokadAUK5cOSIiIi55j6Lk5GTeffddOnfunOVL99+6dOnC8uXLWblyJV27diU4OJj69evz0ksvERERQZ06dS66rLu7e7ar1wqiSZMmLFu2LMd5NWvWxMPDw9l/BzLHgFq/fj0NGjRwTitfvjzDhg3jo48+YsqUKbz77rsFrqdFixbs3r2bWrVqZXtYL3K4PKcr+Vq0aMHevXupUKFCtvX8M5CLSDFkc4fU8+BIhy9GZg7DUQIpCJUxFzrnjho1ikaNGmV53HzzzcyZM+eiy1osFsaOHcsrr7ziPJ3VrVs3EhIS+PPPP7O0PXDgAFu2bOHAgQPY7XZn4EpISHC22b9/P4cPH6ZHjx4Ffj3vvPMO9957L4sXLyYmJoY///yTJ554gj///JM+ffpkaXvixAmOHTvG3r17+fTTT+nYsSOnTp26ZEfarl278vPPP+Pm5ubsU9S1a1fmz59/yaNB1apVY9myZRw7doyzZ88W+HWOHz+e9evXc99997Ft2zZ27drFjBkzOHXqFL6+vtx777089thj/PTTT+zYsYPRo0eTlJTkvK3Bs88+y9dff82+ffv4888/+e6776hfv36B63n22Wf54IMPeP755/nzzz/ZuXMnn376Kf/5z38uuky1atXYtm0bu3fv5tSpU6Snp3P77bcTGhpK3759WbVqFbGxsURHRzN27FgOHSpdN20TKXUsFrjxLQisAmf3Zx4dKuAfxa6kIFTGzJkzhx49euT41/bNN9/Mhg0b2LZt20WXHzZsGOnp6UybNg2AkJAQ+vfvn+3UxrPPPkvz5s2ZMGECCQkJNG/enObNm7NhwwZnm08++YRrr72WqlWrFvj1tGnThoSEBO655x4aNmxIly5d+P3331m0aFG2kFK3bl0iIyNp2bIlr7zyCj169GD79u1Zjprk5KqrrsLhcGRZX9euXbHb7dn6VP3b5MmTWbJkCZUrV77kka/c1KlTh8WLF7N161batGlD+/bt+frrr533/nnllVe4+eabGTJkCC1atGDfvn38/PPPlCtXDsg8GjN+/HiaNGlC586dsdlsfPrppwWup2fPnnz33XcsXryY1q1b065dO958881cP8vRo0dTt25dWrVqRfny5Vm9ejU+Pj6sXLmSKlWqcNNNN1G/fn1GjRpFSkpKnjq7i4iLeZeDAXPAYoPtX8Lmj1xdUb5ZTEHPaZQB8fHxBAYGEhcXl+2XckpKCrGxsVSvXv2yOguXBtu2beOaa64hJiYGPz+/PC2TlpZG7dq1+fjjj7N08hWR4k+//ySbVW/AsufB3QfuiobydV1aTm7f3/+mI0Jy2Zo0acKkSZOIjY3N8zIHDhzgqaeeUggSESkNOj4ENbpCehJ8MQocl9838kop/pc4SYkwfPjwfLW/0CFWRERKAasV+r8LH/aHa54Da95v6eJqCkIiIiJy+fzD4J5fS9yNFktWtSIiIlJ8/TMEnd0P5w66rJS8UhASERGRwhWzHGZ2zry/kD3d1dXkSkFIRERECldw9cx/D62D5f91bS2XoCAkIiIihatcNbhxaubPv74JMb+4tJzcKAiJiIhI4WvYH1qOAAx8dTckFM8BkBWEpESoVq0aU6ZMcXUZ+RIdHY3FYuHcuXOuLqXUWr16NY0bN8bd3Z1+/fpd8e0PHz7cJdsVKTF6vQwVGkDiCVh4NxTDQZAVhMqglStX0qdPHyIjI7FYLCxatChfy3ft2pWHHnqoSGq7koYPH47FYrno42Ij0kvxMW7cOJo1a0ZsbCzz5s274tufOnVqkW93//79WCwWtmzZUqTbESkS7t4w4D1w8848Pbal+A3BoSBUBiUmJtK0aVOioqJcXUoWdrsdRxH9tZCWlpZt2tSpUzl69KjzATB37lzn8/Xr1xdJLfmRU92uVJSfUUHExMRw9dVXU6lSJYKCgq749gMDA12y3YJKTy/eV+9IKVWhPvSeBG3vgSaDXF1NNgpCZVDv3r158cUX6d+//0XbTJ8+ndq1a+Pl5UVYWBgDBgwAMo+irFixgqlTpzqPnOzfvz/HdZw9e5ahQ4dSrlw5fHx86N27N3v37nXOnzdvHkFBQXzzzTc0aNAAT09PDhw4wIkTJ+jTpw/e3t5Ur14924CuAOfOnePOO++kfPnyBAQEcPXVV7N161bn/Oeee45mzZoxe/bsi46HFBgYSHh4uPMBEBQU5Hz++uuvU6dOHXx8fKhRowbPPPOM84tk//79WK3WLIPIAkyZMoWqVateNCx8+eWXNGzYEE9PT6pVq8bkyZOzzK9WrRoTJ05k6NChBAQEcNddd+W4np9++olOnToRFBRESEgIN9xwAzExMc75F44ifPrpp3To0AEvLy8aNWrEihUrnG0unLr7/vvvadKkCV5eXrRr147t27df8jPK7bM9efIk4eHh/Pe//7tS5LfffsPDw4Nly5bl+HpykpqaytixY6lQoQJeXl506tTJGU4vvL7Tp08zcuRILBbLRY/MpKam8sQTT1C5cmU8PT2pVasWc+bMcc5fsWIFbdq0wdPTk4iICJ588kkyMjKc87/44gsaN26Mt7c3ISEh9OjRg8TERCD7qbGuXbsyduxYHn/8cYKDgwkPD+e5557LUs+l9t1/q1498+qb5s2bY7FYsgz0O3v2bOrXr4+Xlxf16tVj+vTpznkX3qMFCxbQpUsXvLy8mD9/vrPm//73v4SFhREUFMQLL7xARkYGjz32GMHBwVSqVIm5c+fm+vmI5EvLYZlhyM3T1ZVkZ+Si4uLiDGDi4uKyzUtOTjY7duwwycnJ2RdMTbj4Iy05H22T8tb2MgBm4cKFWaatX7/e2Gw28/HHH5v9+/ebTZs2malTpxpjjDl37pxp3769GT16tDl69Kg5evSoycjIyHHdN954o6lfv75ZuXKl2bJli+nZs6epVauWSUtLM8YYM3fuXOPu7m46dOhgVq9ebXbt2mUSExNN7969TdOmTc2aNWvMhg0bTIcOHYy3t7d58803nevu0aOH6dOnj1m/fr3Zs2ePeeSRR0xISIg5ffq0McaYCRMmGF9fX9OrVy+zadMms3Xr1ny/FxMnTjSrV682sbGx5ptvvjFhYWFm0qRJzvnXXHONue+++7Kso0mTJubZZ581xhizfPlyA5izZ88aY4zZsGGDsVqt5oUXXjC7d+82c+fONd7e3mbu3LnO5atWrWoCAgLM66+/bvbt22f27duXY61ffPGF+fLLL83evXvN5s2bTZ8+fUzjxo2N3W43xhgTGxtrAFOpUiXzxRdfmB07dpg777zT+Pv7m1OnTmWpr379+mbx4sVm27Zt5oYbbjDVqlW75Gd0qc/2+++/N+7u7mb9+vUmPj7e1KhRwzz88MOX/Az+aezYsSYyMtL88MMP5s8//zTDhg0z5cqVM6dPnzYZGRnm6NGjJiAgwEyZMsUcPXrUJCUl5bieW265xVSuXNl89dVXJiYmxixdutR8+umnxhhjDh06ZHx8fMx9991ndu7caRYuXGhCQ0PNhAkTjDHGHDlyxLi5uZk33njDxMbGmm3btpmoqChz/vx5Y4wxw4YNM3379nVuq0uXLiYgIMA899xzZs+ePeb99983FovFLF682NnmUvvuv61bt84AZunSpebo0aPOdh999JGJiIgwX375pfnrr7/Ml19+aYKDg828efOy7APVqlVztjly5IgZNmyY8ff3N2PGjDG7du0yc+bMMYDp2bOneemll8yePXvMxIkTjbu7uzl48GCONeX6+0/kUjLSjfnjC2McjiLbRG7f3/+mIJSLAgehCQEXf3w0IGvbF8Mv3va967K2nVQ953aXIacg9OWXX5qAgAATHx+f4zJdunQxDz74YK7r3bNnjwHM6tWrndNOnTplvL29zWeffWaMyfySBcyWLVucbXbv3m0As27dOue0nTt3GsAZhFatWmUCAgJMSkpKlm3WrFnTvPPOO8aYzCDk7u5uTpw4kWud/5TTe/FPr732mmnZsqXz+YIFC0y5cuWcdWzcuNFYLBYTGxtrjMkehG677TZzzTXXZFnnY489Zho0aOB8XrVqVdOvX78813zByZMnDWD++OMPY8z/vgRfeeUVZ5v09HRTqVIlZ5i7UN+FUGCMMadPnzbe3t5mwYIFxpicP6O8fLbGGHPfffeZOnXqmNtuu800btw42+eVm4SEBOPu7m7mz5/vnJaWlmYiIyPNq6++6pwWGBiYJUj+24X9acmSJTnOf+qpp0zdunWN4x+/kKOiooyfn5+x2+1m48aNBjD79+/PcfmcglCnTp2ytGndurV54oknjDF523f/7cJnuXnz5mzLfPzxx1mmTZw40bRv3z7LclOmTMlWc9WqVZ2h2Rhj6tata6666irn84yMDOPr62s++eSTHGtSEJICs9uNmdcn87tr/XtFtpn8BCGdGpNsrrnmGqpWrUqNGjUYMmQI8+fPJykpKV/r2LlzJ25ubrRt29Y5LSQkhLp167Jz507nNA8PD5o0aZJtuZYtWzqn1atXL0s/jK1bt5KQkEBISAh+fn7OR2xsbJbTQ1WrVqV8+fL5qvufFixYQMeOHQkPD8fPz4///Oc/HDhwwDm/X79+2Gw2Fi5cCGSeRurWrdtFO1nv3LmTjh07ZpnWsWNH9u7di93+v5GaW7Vqdcna9u7dy+DBg6lRowYBAQHObf6zPoD27ds7f3Zzc6NVq1ZZ3v9/twkODs7zZ3Spz/b1118nIyODzz//nPnz5+PpmfdD4jExMaSnp2d5v9zd3WnTpk22+nOzZcsWbDYbXbp0yXH+zp07ad++PRaLxTmtY8eOJCQkcOjQIZo2bUr37t1p3LgxAwcOZNasWZw9ezbXbf7zvQKIiIjgxInMy4bzuu9eSmJiIjExMYwaNSrLel588cVs68lpf2rYsCHWfwyFEBYWRuPGjZ3PbTYbISEhzrpFCo3VCrWvyfz5pyfh+A7X1oMGXS0aTx25+DzLv0bkfWxfLm3/lVMf+qPgNeWDv78/mzZtIjo6msWLF/Pss8/y3HPPsX79+kLvGOrt7Z3lSygvEhISiIiIIDo6Otu8f9bn6+tb4LrWrFnD7bffzvPPP0/Pnj0JDAzk008/zdKnx8PDg6FDhzJ37lxuuukmPv74Y6ZOnVrgbean7j59+lC1alVmzZpFZGQkDoeDRo0aFUnn6oJ8RpAZZo4cOYLD4WD//v1ZvmivFG9v78ta3mazsWTJEn777TcWL17M22+/zdNPP83atWudfXf+zd3dPctzi8Xi7DOW1333UhISEgCYNWtWlkB6oeZ/yml/yqnG3OoWKVTtxsBfK2DfEvhiBIxeDh4+LitHR4SKgofvxR/uXvlo6523tkXAzc2NHj168Oqrr7Jt2zb279/PL79k3hnUw8MjyxGMnNSvX5+MjAzWrl3rnHb69Gl2795NgwYNLrpcvXr1yMjIYOPGjc5pu3fvznIvnhYtWnDs2DHc3NyoVatWlkdoaGgBX3FWv/32G1WrVuXpp5+mVatW1K5dm7///jtbuzvvvJOlS5cyffp0MjIyuOmmmy66zvr167N69eos01avXk2dOnWyfXnl5sL7+J///Ifu3btTv379ix6l+P33350/X3hf69evf9E2Z8+eZc+ePdna/Pt1XOqzTUtL44477mDQoEFMnDiRO++8M19HF2rWrImHh0eW9ys9PZ3169fnuv/8W+PGjXE4HFk6if/7taxZswZjjHPa6tWr8ff3p1KlSkBmIOjYsSPPP/88mzdvxsPDw3kUML8Ksu96eHgAZPk/FxYWRmRkJH/99Ve29VwsoIkUG1Yr9JsBfuFwchf8+oZLy9ERoTIoISGBffv+dyQqNjaWLVu2EBwcTJUqVfjuu+/466+/6Ny5M+XKleOHH37A4XBQt25dIPPKprVr17J//378/PwIDg7OcpgdoHbt2vTt25fRo0fzzjvv4O/vz5NPPknFihXp27fvRWurW7cuvXr14u6772bGjBm4ubnx0EMPZfnLvkePHrRv355+/frx6quvUqdOHY4cOcL3339P//7983Rq6VJq167NgQMH+PTTT2ndujXff/99jl9+9evXp127djzxxBOMHDky1yMQjzzyCK1bt2bixIkMGjSINWvWMG3atCxX+uRFuXLlCAkJ4d133yUiIoIDBw7w5JNP5tg2KiqK2rVrU79+fd58803Onj3LyJEjs7R54YUXCAkJISwsjKeffprQ0NBcbxKYl8/26aefJi4ujrfeegs/Pz9++OEHRo4cyXfffQfAunXrGDp0KMuWLaNixYrZtuHr68u9997rvIqpSpUqvPrqqyQlJTFq1Kg8v1fVqlVj2LBhjBw5krfeeoumTZvy999/c+LECW655Rbuu+8+pkyZwgMPPMD999/P7t27mTBhAuPGjcNqtbJ27VqWLVvGtddeS4UKFVi7di0nT57MNSjmpiD7boUKFfD29uann36iUqVKeHl5ERgYyPPPP8/YsWMJDAykV69epKamsmHDBs6ePcu4ceMKVJ/IFeNXHm56F7YtgI4PubaWIuupVIz069fPBAUFmZtvvjlfyxW4s3Qxd6GT7L8fw4YNM8Zkdujs0qWLKVeunPH29jZNmjRxdp41JrMDart27Yy3t7cBnJ2D/+3MmTNmyJAhJjAw0Hh7e5uePXuaPXv2OOfPnTvXBAYGZlvu6NGj5vrrrzeenp6mSpUq5oMPPjBVq1bNctVYfHy8eeCBB0xkZKRxd3c3lStXNrfffrs5cOCAMSazs3TTpk3z9b7wr87Sjz32mAkJCTF+fn5m0KBB5s0338yx3gtX3fyzg7cx2TtLG5N5tVeDBg2Mu7u7qVKlinnttdeyLPPv13kxS5YsMfXr1zeenp6mSZMmJjo6Okv9FzrKfvzxx6ZNmzbGw8PDNGjQwPzyyy/Z6vv2229Nw4YNjYeHh2nTpk2WK+wu9hnl9tkuX77cuLm5mVWrVjnbx8bGmoCAADN9+vQs277YvmNM5v+xBx54wISGhhpPT0/TsWPHbO/xpTpLX1jPww8/bCIiIoyHh4epVauWee+9/3XSjI6ONq1btzYeHh4mPDzcPPHEEyY9Pd0YY8yOHTtMz549Tfny5Y2np6epU6eOefvtt53L5tRZ+t8XEvTt29f5f8uYS++7OZk1a5apXLmysVqtpkuXLs7p8+fPN82aNTMeHh6mXLlypnPnzuarr74yxly8k/W/a75Y3bntiyX595+UDfnpLG0x5h/HhEup6Ohozp8/z/vvv88XX3yR5+Xi4+MJDAwkLi6OgICALPNSUlKIjY296D1qpOyYOHEin3/+Odu2bXN1KU779++nevXqbN68mWbNmuXYJjo6mm7dunH27NkSdVNAcT39/pPiLrfv738rE32Eunbtir+/v6vLkFImISGB7du3M23aNB544AFXlyMiIgXg8iCUl3GvoqKiqFatGl5eXrRt25Z169Zd+UJF/uX++++nZcuWdO3aNVu/GxERKRlc3ln6wrhXI0eOzPGKmwULFjBu3DhmzpxJ27ZtmTJlCj179mT37t1UqFABgGbNmmW5Jf4FixcvJjIyMs+1pKamkpqa6nweHx9fgFckZcW8efNcMtBnXlSrVo1LnfXu2rXrJduIiJR2Lg9CvXv3pnfv3hed/8YbbzB69GhGjBgBwMyZM/n+++957733nFfKFNaozC+//DLPP/98oaxLREREij+XnxrLTVpaGhs3bqRHjx7OaVarlR49erBmzZpC39748eOJi4tzPg4ePHjJZfQXtYiUNfq9J6WJy48I5ebUqVPY7XbCwsKyTA8LC2PXrl15Xk+PHj3YunUriYmJVKpUic8//zzLsAIXeHp65nkYgAt3YU1KSrrsu9eKiJQkF4bc+ffdqEVKomIdhArL0qVLC32dNpuNoKAg591yfXx8CjQMgYhISWGMISkpiRMnThAUFJSvO6KLFFfFOgiFhoZis9k4fvx4lunHjx8nPDzcRVX9z4UaNDChiJQlQUFBxeJ3sEhhKNZByMPDg5YtW7Js2TLnLf8dDgfLli3j/vvvL7LtRkVFERUVdcnxtCwWCxEREVSoUIH09PQiq0dEpLhwd3fXkSApVVwehC417tW4ceMYNmwYrVq1ok2bNkyZMoXExETnVWRFYcyYMYwZM8Z5Z8pLsdls+sUgIiJSArk8CG3YsIFu3bo5n18YLHDYsGHMmzePQYMGcfLkSZ599lmOHTtGs2bN+Omnn7J1oBYRERHJrzIx1lhB5WesEhERESkeNNaYiIiISB64/NRYcXbhYJmG2hARESk5Lnxv5+Wkl4JQLs6fPw9A5cqVXVyJiIiI5Nf58+cvedGT+gjlwuFwcOTIEfz9/XO8WWLr1q1Zv359gddfkOXzs0xe2+al3cXaxMfHU7lyZQ4ePFgq+lFd7mdaXLapfVP7ZnHdZnHfN/PaXvvm/xTHfdMYw/nz54mMjMRqzb0XkI4I5cJqtVKpUqWLzrfZbJe1Exdk+fwsk9e2eWl3qTYBAQGl4j/05X6mxWWb2jf/R/tm8dpmcd8389pe++b/FNd9My+3vwF1lr4sY8aMueLL52eZvLbNS7vLfa0lhSteZ1FsU/tm6aN9s+DL53eZwtrvtG+WjG3q1JhcFt1iQIor7ZtSXGnfLF50REgui6enJxMmTMDT09PVpYhkoX1Tiivtm8WLjgiJiIhImaUjQiIiIlJmKQiJiIhImaUgJCIiImWWgpCIiIiUWQpCIiIiUmYpCEmR+u6776hbty61a9dm9uzZri5HxKl///6UK1eOAQMGuLoUEaeDBw/StWtXGjRoQJMmTfj8889dXVKpp8vnpchkZGTQoEEDli9fTmBgIC1btuS3334jJCTE1aWJEB0dzfnz53n//ff54osvXF2OCABHjx7l+PHjNGvWjGPHjtGyZUv27NmDr6+vq0srtXRESIrMunXraNiwIRUrVsTPz4/evXuzePFiV5clAkDXrl3x9/d3dRkiWURERNCsWTMAwsPDCQ0N5cyZM64tqpRTEJKLWrlyJX369CEyMhKLxcKiRYuytYmKiqJatWp4eXnRtm1b1q1b55x35MgRKlas6HxesWJFDh8+fCVKl1LucvdNkaJSmPvmxo0bsdvtVK5cuYirLtsUhOSiEhMTadq0KVFRUTnOX7BgAePGjWPChAls2rSJpk2b0rNnT06cOHGFK5WyRvumFFeFtW+eOXOGoUOH8u67716Jsss2I5IHgFm4cGGWaW3atDFjxoxxPrfb7SYyMtK8/PLLxhhjVq9ebfr16+ec/+CDD5r58+dfkXql7CjIvnnB8uXLzc0333wlypQyqKD7ZkpKirnqqqvMBx98cKVKLdN0REgKJC0tjY0bN9KjRw/nNKvVSo8ePVizZg0Abdq0Yfv27Rw+fJiEhAR+/PFHevbs6aqSpYzIy74p4gp52TeNMQwfPpyrr76aIUOGuKrUMkVBSArk1KlT2O12wsLCskwPCwvj2LFjALi5uTF58mS6detGs2bNeOSRR3TFmBS5vOybAD169GDgwIH88MMPVKpUSSFJilxe9s3Vq1ezYMECFi1aRLNmzWjWrBl//PGHK8otM9xcXYCUbjfeeCM33nijq8sQyWbp0qWuLkEkm06dOuFwOFxdRpmiI0JSIKGhodhsNo4fP55l+vHjxwkPD3dRVSLaN6X40r5ZPCkISYF4eHjQsmVLli1b5pzmcDhYtmwZ7du3d2FlUtZp35TiSvtm8aRTY3JRCQkJ7Nu3z/k8NjaWLVu2EBwcTJUqVRg3bhzDhg2jVatWtGnThilTppCYmMiIESNcWLWUBdo3pbjSvlkCufqyNSm+li9fboBsj2HDhjnbvP3226ZKlSrGw8PDtGnTxvz++++uK1jKDO2bUlxp3yx5NNaYiIiIlFnqIyQiIiJlloKQiIiIlFkKQiIiIlJmKQiJiIhImaUgJCIiImWWgpCIiIiUWQpCIiIiUmYpCImIiEiZpSAkIiIiZZaCkIiIiJRZCkIicsUNHz4ci8WS7fHPwSpFRK4EjT4vIi7Rq1cv5s6dm2Va+fLlszxPS0vDw8PjSpYlImWMjgiJiEt4enoSHh6e5dG9e3fuv/9+HnroIUJDQ+nZsycA27dvp3fv3vj5+REWFsaQIUM4deqUc12JiYkMHToUPz8/IiIimDx5Ml27duWhhx5ytrFYLCxatChLDUFBQcybN8/5/ODBg9xyyy0EBQURHBxM37592b9/v3P+8OHD6devH6+//joRERGEhIQwZswY0tPTnW1SU1N54oknqFy5Mp6entSqVYs5c+ZgjKFWrVq8/vrrWWrYsmWLjoaJuJCCkIgUK++//z4eHh6sXr2amTNncu7cOa6++mqaN2/Ohg0b+Omnnzh+/Di33HKLc5nHHnuMFStW8PXXX7N48WKio6PZtGlTvrabnp5Oz5498ff3Z9WqVaxevRo/Pz969epFWlqas93y5cuJiYlh+fLlvP/++8ybNy9LmBo6dCiffPIJb731Fjt37uSdd97Bz88Pi8XCyJEjsx0Fmzt3Lp07d6ZWrVoFe8NE5PIYEZErbNiwYcZmsxlfX1/nY8CAAaZLly6mefPmWdpOnDjRXHvttVmmHTx40ABm9+7d5vz588bDw8N89tlnzvmnT5823t7e5sEHH3ROA8zChQuzrCcwMNDMnTvXGGPMhx9+aOrWrWscDodzfmpqqvH29jY///yzs+6qVauajIwMZ5uBAweaQYMGGWOM2b17twHMkiVLcnzdhw8fNjabzaxdu9YYY0xaWpoJDQ018+bNy8O7JiJFQX2ERMQlunXrxowZM5zPfX19GTx4MC1btszSbuvWrSxfvhw/P79s64iJiSE5OZm0tDTatm3rnB4cHEzdunXzVc/WrVvZt28f/v7+WaanpKQQExPjfN6wYUNsNpvzeUREBH/88QeQeZrLZrPRpUuXHLcRGRnJ9ddfz3vvvUebNm349ttvSU1NZeDAgfmqVUQKj4KQiLiEr69vjqeDfH19szxPSEigT58+TJo0KVvbiIiIPPetsVgsGGOyTPtn356EhARatmzJ/Pnzsy37z07c7u7u2dbrcDgA8Pb2vmQdd955J0OGDOHNN99k7ty5DBo0CB8fnzy9BhEpfApCIlKstWjRgi+//JJq1arh5pb9V1bNmjVxd3dn7dq1VKlSBYCzZ8+yZ8+eLEdmypcvz9GjR53P9+7dS1JSUpbtLFiwgAoVKhAQEFCgWhs3bozD4WDFihX06NEjxzbXXXcdvr6+zJgxg59++omVK1cWaFsiUjjUWVpEirUxY8Zw5swZBg8ezPr164mJieHnn39mxIgR2O12/Pz8GDVqFI899hi//PIL27dvZ/jw4VitWX+9XX311UybNo3NmzezYcMG7rnnnixHd26//XZCQ0Pp27cvq1atIjY2lujoaMaOHcuhQ4fyVGu1atUYNmwYI0eOZNGiRc51fPbZZ842NpuN4cOHM378eGrXrk379u0L540SkQJREBKRYi0yMpLVq1djt9u59tprady4MQ899BBBQUHOsPPaa69x1VVX0adPH3r06EGnTp2y9TWaPHkylStX5qqrruK2227j0UcfzXJKysfHh5UrV1KlShVuuukm6tevz6hRo0hJScnXEaIZM2YwYMAA7rvvPurVq8fo0aNJTEzM0mbUqFGkpaUxYsSIy3hnRKQwWMy/T5qLiJQCXbt2pVmzZkyZMsXVpWSzatUqunfvzsGDBwkLC3N1OSJlmvoIiYhcIampqZw8eZLnnnuOgQMHKgSJFAM6NSYicoV88sknVK1alXPnzvHqq6+6uhwRQafGREREpAzTESEREREpsxSEREREpMxSEBIREZEyS0FIREREyiwFIRERESmzFIRERESkzFIQEhERkTJLQUhERETKrP8DF6vIczmKPUcAAAAASUVORK5CYII=", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "# Bias\n", + "fs = 500.\n", + "freqs = np.linspace(.5, fs/2, 2000)\n", + "phi = 0.9\n", + "\n", + "plt.loglog(freqs, 1/(1 -2 * phi * np.cos(2*np.pi*freqs*1/fs) + phi**2), label='AR(1) PSD with cosine term')\n", + "plt.loglog(freqs, 1/(1 -2 * phi * (1 - ((2*np.pi*freqs*1/fs)**2 * 0.5)) + phi**2), label='1st order Taylor approx. of cosine term', ls='--')\n", + "plt.xlabel(\"Frequency\")\n", + "plt.ylabel(\"Power\")\n", + "plt.title(\"PSD Bias\")\n", + "plt.legend();" + ] + }, + { + "cell_type": "markdown", + "id": "9c9c3bfc-ab16-4110-a565-a17c45edba10", + "metadata": {}, + "source": [ + "#### Rearrange and factor the quadratic term\n", + "\n", + "$$\n", + "\\begin{align}\n", + "P(f) & \\approxeq \\frac{1}{1 - 2\\varphi (1 - \\frac{(2 \\pi f \\frac{1}{f_s})^2}{2}) + \\varphi^2} \\\\\n", + "&= \\frac{1}{1 - 2\\varphi + \\varphi (\\frac{2 \\pi f}{f_s})^2 + \\varphi^2} \\\\\n", + "&= \\frac{1}{\\varphi^2 - 2\\varphi + 1 + \\varphi (\\frac{2 \\pi f}{f_s})^2} \\\\\n", + "&= \\frac{1}{(1-\\varphi)^2 + \\varphi (\\frac{2 \\pi f}{f_s})^2} \\quad\\quad\\quad \\text{Factored quadratic term}\\\\\n", + "&= \\frac{\\frac{1}{\\varphi}}{\\frac{(1-\\varphi)^2}{\\varphi} + (\\frac{2 \\pi f}{f_s})^2} \\quad\\quad\\quad \\text{Multiply each term by } \\frac{1}{\\varphi} \\text{ to isolate } f\n", + "\\end{align}\n", + "$$\n", + "\n", + "Let constant scaling factors be absorbed in $b$, e.g. the numerator of where we left off. This term translates the PSD along the y-axis and allows us to multiply all terms with constants, e.g. $(\\frac{2 \\pi}{f_s})^2$:\n", + "\n", + "$$\n", + "\\begin{align}\n", + "P(f) &\\approxeq \\frac{b}{\\frac{(1-\\varphi)^2}{\\varphi} + (\\frac{2 \\pi f}{f_s})^2} \\\\\n", + "&= \\frac{b}{\\frac{(1-\\varphi)^2}{\\varphi} (\\frac{f_s}{2 \\pi})^2 + f^2} \\\\\n", + "&= \\frac{b}{(\\frac{(1-\\varphi)}{\\sqrt{\\varphi}} \\cdot \\frac{f_s}{2 \\pi})^2 + f^2} \\\\\n", + "\\end{align}\n", + "$$\n", + "\n", + "Next, $\\frac{(1-\\varphi)}{\\sqrt{\\varphi}}$ needs to be approximated using $-\\ln(x)$. This leads to the commonly used Lorentzian form. Take the first order approximations of these two forms around $\\varphi = 1$:\n", + "\n", + "$$\n", + "\\begin{align}\n", + "\\frac{(1-\\varphi)}{\\sqrt{\\varphi}} \\approxeq 0 - 1 (x+1) = -x-1 \\\\\n", + "-\\ln(\\varphi) \\approxeq 0 - 1 (x+1) = -x-1\n", + "\\end{align}\n", + "$$\n", + "\n", + "Since we took the approximations around $\\varphi = 1$, bias increases as $\\varphi \\to 0$." + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "id": "a0269b41-f272-4b44-a89d-d6a4d0508989", + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAloAAAHLCAYAAAAQiTaRAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8g+/7EAAAACXBIWXMAAA9hAAAPYQGoP6dpAACJmElEQVR4nOzdd3wUVdfA8d/uppNOGmmElpBQAoReA4Qu0hQEFVBEUaSIWNBHQV+7wqM+RGx0C1FBUHpvoZfQQidAgBQgjfRkd94/IivLJkBCkk05389nTXbumZmzmzU53HvnjkpRFAUhhBBCCFHq1KZOQAghhBCiqpJCSwghhBCijEihJYQQQghRRqTQEkIIIYQoI1JoCSGEEEKUESm0hBBCCCHKiBRaQgghhBBlRAotIYQQQogyIoWWEEIIIUQZkUJLCCGqmdDQUEJDQ02dhhDVghRaQohytWDBAlQqFSqVip07dxq1K4qCj48PKpWKRx55xAQZFs3Pz0+fu0qlws3NjU6dOvHnn38axOl0OhYtWkSbNm1wdnbGzs4Of39/Ro4cyZ49e/RxW7duNTiepaUl7u7uhIaG8tFHH3H9+vUHyuvixYsGx1GpVNjb29OsWTNmz56NVqst1fdBCPHgzEydgBCierKysuKXX36hY8eOBtu3bdvGlStXsLS0NFFm99asWTNeffVVAK5du8Z3333H4MGDmTNnDuPGjQNg4sSJhIeHM2DAAJ588knMzMw4ffo0a9asoW7durRt29bgmBMnTqRVq1ZotVquX7/Orl27mD59OrNmzeK3336jW7duD5Tb8OHD6du3LwCpqamsXr2aCRMmcOnSJT7//HN93Pr160vjrRBCPAhFCCHK0fz58xVAGTx4sOLi4qLk5eUZtI8dO1YJCQlRateurfTr169UzqnVapWsrKyHPk5hOcXFxSk1atRQ/P39FUVRlPj4eEWlUiljx4412l+n0ykJCQn651u2bFEA5ffffzeKjYqKUtzc3BRHR0fl2rVr98wrJiZGAZTPP//c6HytWrVSPD09H/g1CiFKlwwdCiFMYvjw4dy8eZMNGzbot+Xm5vLHH38wYsSIQvfJyMjg1VdfxcfHB0tLSwICAvjiiy9QFMUgTqVS8fLLL/Pzzz/TqFEjLC0tWbt2LQBXr17l2Wefxd3dHUtLSxo1asS8efNK/Do8PDwIDAwkJiYGgJiYGBRFoUOHDkaxt4cbH0RwcDBffvklKSkpzJ49u0S5qVQq3N3dMTMzHLy4e45Wbm4u7777LiEhITg4OFCjRg06derEli1bjI65ZMkSQkJCsLOzw97eniZNmvDVV1+VKD8hqgMptIQQJuHn50e7du349ddf9dvWrFlDamoqTzzxhFG8oig8+uij/Pe//6V3797MmjWLgIAAXnvtNaZMmWIUv3nzZl555RWGDRvGV199hZ+fHwkJCbRt25aNGzfy8ssv89VXX1G/fn3GjBnDl19+WaLXkZeXR2xsLDVr1gSgdu3aAPz+++9kZmaW6Ji3PfbYY1hbWz/wUF9mZiY3btzgxo0bXLhwgfDwcNauXcuoUaPuuV9aWho//vgjoaGhfPrpp8yYMYPr16/Tq1cvoqKi9HEbNmxg+PDhODk58emnn/LJJ58QGhpKZGTkw7xMIao2E/eoCSGqmdtDh/v371dmz56t2NnZKZmZmYqiKMrjjz+udO3aVVEU42G65cuXK4DywQcfGBzvscceU1QqlXLu3Dn9NkBRq9XKiRMnDGLHjBmj1KpVS7lx44bB9ieeeEJxcHDQ51GU2rVrKz179lSuX7+uXL9+XTly5IjyxBNPKIAyYcIEfdzIkSMVQHFyclIGDRqkfPHFF8rJkyeNjnevocPbgoODFScnp3vmdXvosLDHiy++qOh0OoP4Ll26KF26dNE/z8/PV3JycgxikpOTFXd3d+XZZ5/Vb5s0aZJib2+v5Ofn3zMfIcS/pEdLCGEyQ4cOJSsri5UrV3Lr1i1WrlxZ5LDh6tWr0Wg0TJw40WD7q6++iqIorFmzxmB7ly5dCAoK0j9XFIWlS5fSv39/FEXR9/zcuHGDXr16kZqayqFDh+6b8/r163F1dcXV1ZXg4GB+//13nn76aT799FN9zPz585k9ezZ16tThzz//ZOrUqQQGBtK9e3euXr1anLcIW1tbbt269UCxzz//PBs2bGDDhg0sXbqU8ePH89133xXa43cnjUaDhYUFUHDFZFJSEvn5+bRs2dLgPXF0dCQjI8NguFcIcW9y1aEQwmRcXV0JCwvjl19+ITMzE61Wy2OPPVZo7KVLl/D09MTOzs5ge2BgoL79TnXq1DF4fv36dVJSUvj+++/5/vvvCz1HYmLifXNu06YNH3zwASqVChsbGwIDA3F0dDSIUavVjB8/nvHjx3Pz5k0iIyP59ttvWbNmDU888QQ7duy473luS09PN3rNRWnQoAFhYWH654MHD0alUvHll1/y7LPP0qRJkyL3XbhwITNnzuTUqVPk5eXpt9/5Pr700kv89ttv9OnTBy8vL3r27MnQoUPp3bv3A78eIaobKbSEECY1YsQIxo4dS3x8PH369DEqWkrK2tra4LlOpwPgqaeeKnLOUtOmTe97XBcXF4Ni5n5q1qzJo48+yqOPPkpoaCjbtm3j0qVL+rlc95KXl8eZM2do3LjxA5/vbt27d2f27Nls3769yELrp59+YvTo0QwcOJDXXnsNNzc3NBoNH3/8MefPn9fHubm5ERUVxbp161izZg1r1qxh/vz5jBw5koULF5Y4RyGqMim0hBAmNWjQIF544QX27NlDREREkXG1a9dm48aN3Lp1y6CH59SpU/r2e3F1dcXOzg6tVlusQqk0tWzZkm3bthEXF/dAhdYff/xBVlYWvXr1KvE58/PzgYKesXudp27duixbtgyVSqXfPn36dKNYCwsL+vfvT//+/dHpdLz00kt89913vPPOO9SvX7/EeQpRVckcLSGESdna2jJnzhxmzJhB//79i4zr27cvWq3WaKmD//73v6hUKvr06XPP82g0GoYMGcLSpUs5fvy4UfuDrsJ+P/Hx8URHRxttz83NZdOmTajV6gcqSI4cOcLkyZNxcnJi/PjxJc7n77//BgqWiyiKRqMBMFgmY+/evezevdsg7ubNmwbP1Wq1vhcwJyenxDkKUZVJj5YQwuTut/wAQP/+/enatStvv/02Fy9eJDg4mPXr17NixQomT55MvXr17nuMTz75hC1bttCmTRvGjh1LUFAQSUlJHDp0iI0bN5KUlPTQr+XKlSu0bt2abt260b17dzw8PEhMTOTXX3/VF08uLi4G++zYsYPs7Gy0Wq1+Ttdff/2Fg4MDf/75Jx4eHg907kOHDvHTTz8BcOvWLTZt2sTSpUtp3749PXv2LHK/Rx55hGXLljFo0CD69etHTEwM3377LUFBQQY9Yc899xxJSUl069YNb29vLl26xP/+9z+aNWumnysnhLiLaS96FEJUN3cu73Avha3CfuvWLeWVV15RPD09FXNzc6VBgwbK559/brR8AaCMHz++0OMmJCQo48ePV3x8fBRzc3PFw8ND6d69u/L999/fN/cHWa0+LS1N+eqrr5RevXop3t7eirm5uWJnZ6e0a9dO+eGHHwxyvb28w+2Hubm54urqqnTu3Fn58MMPlcTExPvmpCiFL+9gZmam1K1bV3nttdeUW7duGcTfvbyDTqdTPvroI6V27dqKpaWl0rx5c2XlypXKqFGjlNq1a+vj/vjjD6Vnz56Km5ubYmFhofj6+iovvPCCEhcX90B5ClEdqRTlriWVhRBCCCFEqZA5WkIIIYQQZUQKLSGEEEKIMiKFlhBCCCFEGZFCSwghhBCijEihJYQQQghRRqTQEkIIIYQoI7JgqQnpdDquXbuGnZ2dwW0vhBBCCFFxKYrCrVu38PT0RK2+d5+VFFomdO3aNXx8fEydhhBCCCFKIDY2Fm9v73vGSKFlQrdvjBsbG4u9vb2JsxFCCCHEg0hLS8PHx8fgBvdFkULLhG4PF9rb20uhJYQQQlQyDzLtRybDCyGEEEKUESm0hBBCCCHKiBRaQgghhBBlRAotIYQQQogyIoWWEEIIIUQZkUJLCCGEEKKMSKElhBBCCFFGpNASQgghhCgjUmg9pJUrVxIQEECDBg348ccfTZ2OEEIIISoQWRn+IeTn5zNlyhS2bNmCg4MDISEhDBo0iJo1a5o6NSGEEEJUANKj9RD27dtHo0aN8PLywtbWlj59+rB+/XpTpyWEEEKICqJaF1rbt2+nf//+eHp6olKpWL58uVFMeHg4fn5+WFlZ0aZNG/bt26dvu3btGl5eXvrnXl5eXL16tTxSF0IIIUQlUK2HDjMyMggODubZZ59l8ODBRu0RERFMmTKFb7/9ljZt2vDll1/Sq1cvTp8+jZubW7HPl5OTQ05Ojv55WlraQ+VflC2nElm+9X/cUK0uMsZV6UeK5SMAOOZu4AbL7or490aZbnQnxfIxAOxzd3BT+aWISHClM6lWwwEVtrkHSNbNLTLWhbakWY8CwDrvBKn54fc4bgvSbMYCYKU9T1rOTMNA5Y7jqhqTUWM8AGb5caRn/5/R67+tpiqATNtJAGh0qWRmvHVXxL9ZOKnqkmP/asETXQ5Z6VOKiFThiA+5Dm/qt2SnTkBBV2i0PbXId3wb1T/bclImoyW/0HxtVa7gOF3/PC/5NXLJLuSoYKNyQu1U8NpVKshLmkaOko7qnzMVPAq+t1TZoXb5ELWqYBs3/49cXTKgLohTqfWx5uoaWHu+j1oFapWKvPgvyclP/OfmqmrUqn8iVWo0Kkucav+nYJsKsuLmkZ17DTVqVCqV/rhqlRq1SoNL3df0sTnxK8jPvoaZSo1arUGjtkCjMUejNsdMbYZd7SFYmFugVqvQpRxHlZOEucYcjZkF5maW/zzMMTezQu3gi7m5ORq1CjNtNmq0mJmZozEzx8zMHLVGU+j7LYQQD6taF1p9+vShT58+RbbPmjWLsWPH8swzzwDw7bffsmrVKubNm8ebb76Jp6enQQ/W1atXad26dZHH+/jjj3nvvfdK7wUUIS41m/jki0R75hYZE3TtIntTrwPQ0v4ip72Kjm2ScJFdSYkABNte5IJPTpGxwYkX2BlTEBtkE0Ns7ewiY5vdOM+OmAQAGlhdJL5OVpGxzZPOsT0mHgBfi0sk18ssOjb5HNsvxAHganaJ7AYZReeQepYd/8Taqa9DQHqRsU1vnSEy6hoA5mRjFXiryNjG6WfZffjfz4ZTw2Tyi7jLe8PMDPYf+jfWIyCRDHXhnc31s1I5fPCK/rlfg6vcNCs81i8nmWMHYvXP/etdJM6y8FjP3Ouc3nv53/zrRnPJstBQXPJ1xGw9r3/evE4k56yUQmNt83XErTmlf96y9kpO22gLjTVTFJKX9dY/b+sznxO2RX8us7Z6kU9Bkh283+OoXdGfH05P5ZbOBYDuHh9z3i4JcxQslIKHuQLmioKZouLK9de4pfbDwkxNG4tv0VqcxAw1GkWDBjVmmGGGBnOVObfsX0FnUwdzjZpa6euxz4rCQm2NhXkNrMztsLK0w8rCHitre3I822BRwxEbCzNslEysVblY1bDD2toWjVm1/lUsRJUm/3cXITc3l4MHDzJt2jT9NrVaTVhYGLt37wagdevWHD9+nKtXr+Lg4MCaNWt45513ijzmtGnTmDLl316QtLQ0fHx8Sj33tnWdUTIep/1NZ/02RTH8Q+jSOozHHYMByElVkXjDEqXwv5XUbN6FgU5NAci9ZUFi4r+9LXcft2aTDjzi3ASAvAxr4uNTCj2mAjg3bEOfjgWx+VkOxF27YhBx55Gd6regZ7vGKIA2uybx107pD2QYCQ5+TQhr0wgAXW4trl35d7j37nztvQPpHhKEAih5t7h2pfkdKRjG2rrXp1uzoIImXS7XLgXf9YruiHWpTfemgfrDXLsYDP/0aN39Nts4eRIW1FD/POFSc7RKHndTAGs7N3r2/jc28XJL8pVMfcSdx7a2caZ3rwD985uxrcjR3UKHoo/U/fPV0tKOfj380SkKOgVuXmtKfW0yyj8RivLPXoqChcqarh3qoPvn/bl13R/X3OsFcf/895+9MMOc9i28URSl4NhpvljlJvwbqyjoVAU/FzUqagS6oSigUxTMclzxzYtHC+hUSsHXf46tVYGPhxNaRUO+ToeitkSjZKItopjVqC3QoEKrU8jR5HPTrOgeLPNshSRtwT8m6taK47BtwTuFvkfy339o2J09z7Xcgl+jXdxXccg58d8fVu4/DwqKSOcdwzif3QKAPo4/oDgcx06nw1anw1qnwkqnxlwxx0Kx4Kjly6TbNsfOyoyGeUdolHUQMytH1FYOaGwcMLdxxMLWCRtHV2zd/HC0t0ejLvy1CyFMSwqtIty4cQOtVou7u7vBdnd3d06dKvgjb2ZmxsyZM+natSs6nY7XX3/9nlccWlpaYmlZRDdBKarrakvdbo8Cjz7gHt5A3weM9QF6PGCsL9CtGLEdHzC2NtD2AWP9gAUPGAuwqBixPz14aOefy+a4LCxG7IJixBYn3z+KEftXMWLXFiN2B8A/BZ0OnaIjX8lHq9OiVbTYW9ijUqlQFIW4tGbczLxOTl4WWXnZ5OZlk52XRU5+Ftm5mTQLG4AOS/K0Ok5cGoN/8hFy83PI0+WQp80lT5tHri6XHCWHlp3aotO4kpuv43JCfVzzU8lFS45KRw46lH9qn3yVCicHV7ytrcnK1ZJhlclhG+siXouCe8w1zl0pmP/ZxeV3wl0vY6vV4ZSmxSlZh5NOh5NWi5NWR2TSSA7ltsXR2py+lpH0163HXONMnlVNdNbOqGxc0Ni5Yl3Thxq+zXB1dcXOyrwY760Q4mFIofWQHn30UR599EELGiFEWVKpVGhUGjRoMMe4mFCpVHg6eOHp4FXI3saa+44uxtm/M3imKAp5ujyy8rPIys+ipnVNzNUFOZ1K8uHUzWiSM5NIyUgmNTuZW9mp3MpJ41ZuGgP79Aa1J7ey8zlw2QPyL5OuVpOuVhN718tyTTFDUSA5M49r1gcZ65GOpS4NV+153LRa3FK0uN3U4nZOy8p1z7I/tw01LDQMsTnMY8p6si1dybP1QGXvhZVbPRw86+Ph04AaNWoU47ULIYoihVYRXFxc0Gg0JCQkGGxPSEjAw8PDRFkJISoLlUqFhcYCC40FDpYOBm0NnRvS0LlhEXsaelb5gbScNJJzkknJSSEpO4mU7BSSc5JJzk5m1GPPoNLZkZyZy8JDOyHxEjlqNVfUaq6YG1ZltbKtIRcycrUk2O7iQ8ereOVfxjsrH+9b+dS+mIdDXh5W+QrPa6aTULM1Pk7WhFheoaE6FlvvRtSq15SaTk7/XPwghLgfKbSKYGFhQUhICJs2bWLgwIEA6HQ6Nm3axMsvv2za5IQQ1YZapcbRyhFHK8d7xrnaWfJhn495VzuD65nXScxMNHgkZCbwnydewowaJKRlE753PbtuXiS6kNkMdlod6svZXIlN4UhsCu6WEfQyX0OtI1pUwDVcSLCoTYZdXRTXANRB/anjW5taDlZSgAlxl2pdaKWnp3Pu3Dn985iYGKKionB2dsbX15cpU6YwatQoWrZsSevWrfnyyy/JyMjQX4UohBAVjaXGEm87b7ztvIuMqetqyytdXuOR5IFcSb/C1fSrxN6KJSY1hqvpV7mlgR9HDSc5w5LYpCz2n8+ll9oLe62OoNwcAnPyCMo9SWDqEXxu5tPjiAPnFS+cbMx50vkUIdZxWPs0x6NhG3y9fVHLRH1RjVXrQuvAgQN07dpV//z2FYGjRo1iwYIFDBs2jOvXr/Puu+8SHx9Ps2bNWLt2rdEEeSGEqGyKKsZytDlcTrtMAyc//bZr5o0xO3uONE0+e6yt2WP970R+a50Kh2QXzG6oSM7Mwy9vLV01kRAL7II4pSYx1o3IdG+JQ0AHGjTtgKNtURcCCFH1qJS7r3cX5SYtLQ0HBwdSU1Oxt7c3dTpCCFGkXG0uZ1POcurmKU4mneTkzZOcTj6NrbktW4ZuISdfx5mEW8zaNpLr2VdpmpVNx8wk2mdn46QrWBpDq6homvMjtdxcaVnbiU5u2YQEBeBR0+E+ZxeiYinO3+9q3aMlhBDiwVhoLGhUsxGNajbSb8vX5XM98zoqlQorcw1NvBy4ZH6LG/m5XLZTs9LOBRUqamNPcGY+gbcgI8eac4npnEtMZ4jFDOw3XmKveWOS3Tvg0KgHjZu3xc7awoSvVIjSJT1aJiQ9WkKIquZ65nUOJx7mcOJh9sbv5WzyWX1b45qNmR26gEOXUzgUk8jwQwPw1d402D9BceSoTXtyA/oT1HEAdVxkmQlR8RTn77cUWiYkhZYQoqq7kXWD3dd2syduD/5O/oxqVHB/04y8DEIjQgmo4Utwjg1t4q7QOvU4Vv8sp79BG8LYvFep72ZLjyB3evk70NTPQybWiwpBCq1KQgotIUR1FXk1knEbxxlsa1qzCaGWdQm+FMvuzMb8L74R+ToFL66zwfJ1dmtaktJgEI07DybAq+i7cAhR1qTQqiSk0BJCVGdx6XFsid3C5sub2Z+wH51SMGneTGXGx50/pp17d7aeTiR717cMS/xav1+yYstOy07kNh5Oh8498XCUqxhF+ZJCq5KQQksIIQrcyLrBmpg1rLywkpM3T7JuyDpq2dYC4GzSGVSJ51AfWEPNmL9x0Cbr9zuu8yPC5126depEZ39Xubm2KBdSaFUSUmgJIYSxa+nX8LT11D8ft3EckVcj6eDZgWENHqNlej43dy3C89p6chUNbXLCycAaL0drRoS48Xjb+rjZWZnwFYiqTpZ3EEIIUWndWWTl6fJQo0aFishrkURei8TP3o/RoaOpVWsmqWcO8URcHZYeusLVlEw67HiaYzucOF13NN17DSSglvwjVpiW9GiZkPRoCSHEg4m9Fcvvp3/njzN/cCvvFgA1rWoyqcUkBjUYRHaelshd2+i6ZTBqCv6sHdbVZ4fbCFr0fJoODVzlPoyi1MjQYSUhhZYQQhRPRl4Gf5z5g8XRi0nITGBGuxkM8R/yb8CNs1xfPxPHs39gruQBEK2rzZ+Oo+jY72k6+0vBJR6eFFqVhBRaQghRMnm6PNZdXEeP2j2w1FgCsO7iOuIz4hkWMAyr7DRSt/4Py0M/YqXLBGB07mukeXfllR7+dKzvIgWXKDEptCoJKbSEEKJ05Gnz6L+8P1fTr+Jm7cbkkMn0q9sPdVYKGVtnkXxiC2Gpb5GdXxDfydeGqf2bE+zjaNK8ReUkhVYlIYWWEEKUjnxdPn+f/5s5R+YQlxEHQBOXJrzR+g2CXYNBpyMxPZc5286zbO8ZVmmmsl3blOiGE3jxkXZ4yVpcohik0KokpNASQojSlaPNYXH0Yn44+gOZ+QVDhv3q9mNqy6m4WLsAkLxvCU6rXwAgTbEmXPcY5h3G8WLXhtSwlIvxxf0V5++3upxyEncIDw8nKCiIVq1amToVIYSoUiw1ljzX5DlWDlrJwPoDUaFi1YVVnEo6pY9xav0EPLOWLJcm2KuymKZZTK9dT/LSF/NYdyLehNmLqkh6tExIerSEEKJsnbhxgu1Xt/Ni8Iv6bVn5WVibWYNOh3J4Mflr/4N5Xhr5ipq52j4cqD+J6QOa4O1kY8LMRUUmPVpCCCEE0MilkUGRlZCRQO+lvVl4YiE6FahCRmE+8QD5gYMwU+nwU19nw6kb9Ji1nR93XECnk74I8XCk0BJCCFFt/HX+L5Kyk/jiwBeMXT+W+Ix4sHPHbNgCGPEbDUZ/Q2s/Z7LytMxadZgnv99JbFKmqdMWlZgMHZqQDB0KIUT5UhSFpWeX8tn+z8jKz8LOwo532r5Dnzp9DGJ+2XcZh9Uv4qPE8RYTeOqRMJ5o5SNrbwlArjqsNKTQEkII07iUdolpO6Zx7MYxAIY0GMK0NtP0i5+SehXtN+3R5KSQoVgyLW8stxoM4IvHg6lpa2nCzEVFIHO0hBBCiHuobV+bhX0WMi54HCpULD27lHnH5v0b4OCFZvxuFL/O1FDl8LXFbLqe/4yBX21hX0yS6RIXlY70aJmQ9GgJIYTp7bq2i3nH5/G/bv8ruBrxTjotbPkIdnwBQJSuHhPyJ/FEjw682KUearUMJVZHMnRYSUihJYQQFY+iKGy+vJluvt3+nZN1Zh3KsudRZadwTudJz9zP6OjvztdPNMPRxsK0CYtyJ0OHQgghRAl9c+QbJm+dzNs73yZXm1uw0b8Xqhe2o3i1JLb9h1iYm7H9zHUGhEdyJuGWaRMWFZoUWkIIIcQdXK1d0ag0/H3hb8auH0tydnJBg1NtVM9tpGvvwfz5Uge8naxRJ51ncPgONkQnmDZpUWFJoSWEEELcYWjAUL4J+wY7czsOJR5ixKoRxKbFFjT+M5QYWMueVU+4sMrqHf5Pmc1Li/cwe/NZZDaOuJsUWkIIIcRd2nu2Z3HfxXjZenEl/Qoj147kdNJpgxiHW+ewVuUySBPJ92Yzmb3+GK//cZQ8rc5EWYuKSAotIYQQohD1HOvxU9+f8Hfy50bWDcauH0tGXsa/AY0Hoxq+BMys6ao5wk8WH7Pu4CmeW3iAjJx80yUuKhQptIQQQogiuFi7ML/3fFq4teD11q9Tw7yGYUCDHjByOVg50FJ9ht8tPyD6zFmGfb+bxFvZJslZVCyyvIMJyfIOQghROegUHWrVv30TWp0WjVrzb0DCCVg8GNLjicGLx7L/g42zBz+NaUPtmjUKOaKozGR5ByGEEKIU3VlkJWQk8Njfj7H9yvZ/A9wbwbNrwN4bD+86uDg7E5uUxdDvdnMuMd0EGYuKQgotIYQQohgWRi/kXMo5Jm+ZzM6rO/9tcK4Lz67BeuTv/PRiVwLc7UhIy2HYd7s5GZdmuoSFSUmhJYQQQhTDKyGv0N23O3m6PCZtnsSua7v+bXT0BQsbXO0s+XVsG95w3oYu4wbDf9jD0SspJstZmI4UWkIIIUQxmKvN+bzz53Tz6UauLpfJWyZz7PoxozjnI9/yYuZ3RNjOIjfzFk/+sJeo2JTyT1iYlBRaQgghRDGZa8z5vMvntK3Vlqz8LF7a9BIXUi8YBvn3Bmsn/PPP8Iv9bLJzshk5dy8nrqWaJmlhElJoCSGEECVgobHgy65f0rhmY1JyUvi/3f9nGOAaAE/+AeY2NMs9xHyHH7mVncvTc/dxVu6PWG1IoSWEEEKUUA3zGnwT9g09avfgs86fGQd4t4Rhi0FtRsec7Xzl+BtJGbmM+HEvMTcyjONFlSOFlhBCCPEQnKycmBU6C1cb18ID6ofBoO8AeDT7L15z2s71Wzk8+cMerqZklWOmwhSk0BJCCCFK0eoLq/l036eGN5hu8hh0nw4qNSPb+lLXtQbXUrMZPW8fKZm5pktWlDkptEwgPDycoKAgWrVqZepUhBBClKJLaZeYtnMaP538iZ9P/mzY2PEVeH4bdl3G89OYNnjYW3E2MZ2xiw6Qnac1TcKizMkteExIbsEjhBBVz8ITC/niwBeoVWr+1+1/dPbuXGjcmUuxjJm/h9hsG3o1cuebJ0PQqFXlnK0oCbkFjxBCCGEiI4NGMqTBEHSKjje2v8GltEvGQUkx+P89hJXuc7HWKKw7kcCMv04gfR9VjxRaQgghRClSqVS83fZtWri1ID0vnVe2vkJmXqZhUH4OpF7BIWE3fzfagkoFi/dc4scdMaZJWpQZKbSEEEKIUmauNueLLl/gYu3C2eSzvLf7PcPeKreGMGA2APXP/Mj3LeMA+GjNSTafSjBFyqKMSKElhBBClAFXG1e+6PIFZiozPGp4oHDXsGDjwdB2PABhp2cwIVhBUWDir1GckQVNqwyZDG9CMhleCCGqvthbsfjY+RTeqM2DhY/C5V3oXBsySvMJOy5m4uNszYrxHXGuYVG+yYoHIpPhhRBCiAriziIrX5dvOF9LYw6PLwBbd9TXT/GD9zp8nW2ITcpi3E8Hyc3XlX/ColRJoSWEEEKUg6vpVxm9drTxfC07dxj0LdTpglXnicwd1RJbSzP2xSTxfyujTZewKBVSaAkhhBDl4EbWDY7fOM7qmNX8feFvw8Z63WDkCrD3pIG7HV8Pb6a/EvHPw1dMk7AoFVJoCSGEEOUg2DWYl5q9BMCHez7kctplwwDVv4uVdrM8w8SudQGYtuwYJ+PSyi1PUbqk0BJCCCHKyZjGY2jl0YrM/Exe3/46ebo846D1/4GFjzDJchWd/V3JztPx4k8HSc0qJFZUeFJoCSGEEOVEo9bwUcePsLew58TNE8w/Pt84yC0IAPW2T/hfNwu8HK25eDOTqb8fQaeThQIqGym0hBBCiHLkUcODaW2mATDnyBxOJ502DAgeDgF9QZeHw9oJzBnRBAuNmg3RCXy/44IJMhYPQwotIYQQopz1q9OPrj5daeDYADO1mWGjSgWPfAnWThB/lKYX5vHegEYAfLHuNIcvJ5d/wqLEZMFSE5IFS4UQovpKy03D2swac7V54QHH/oClY0BthjJ2MxO2aFl5NA4fZ2tWTeyEvVUR+4kyJwuWCiGEEBWcvYW9QZGl1WkNAxoPgcD+oMtHtXw8Hz7aEG8na2KTsnj7z+NIP0nlIIWWEEIIYULZ+dn87/D/eGbdM4bFlkoF/f4L7o2h81Qcaljz1RPN0ahV/H3kGn8clPW1KgMptIQQQggTSs9L59eTv3I48TBLTi8xbLR1hXE7odFAUKkIqe3ElB7+ALy74gTnr6eXf8KiWKTQEkIIIUzIxdqFySGTAfjf4f+RkJFgGHDHQqZkpzKuc13a16tJVp6Wib8elvshVnBSaAkhhBAm9pj/YzR1bUpGXgaf7v+08KCjv8PXzdGc/JP/DmuGo405J66l8b/NZ8s3WVEsUmgJIYQQJqZWqXm37btoVBo2XNrAtthtxkHJMZB5E9a8ibt5Nh8ObALAN1vPExWbUr4JiwcmhZYQQghRAQQ4BzCy0UgAPtz7IZl5mYYBHSaBiz9kJMKm9+jXtBaPBnui1SlM+S2KrFxtIUcVpiaFlhBCCFFBjGs6Ds8aniRnJ3Pi5gnDRjNLeOS/Bd8fmA9xR3l/QCPc7Cy5cD2DT9eeKv+ExX1JoSWEEEJUEDbmNnza+VP+GvgXrTxaGQf4dYRGgwEF1ryBo7U5nz3WFIAFuy6y69yN8k1Y3JcUWkIIIUQF0sytGbVsaxUd0ON9MLOGy7vgxJ+EBrgxoo0vAFN/P0Jadl45ZSoehBRaQgghRAV1MOEgR68fNdzo6AMdJxd8H38MgLf7BuLrbMO11Gw+WSNDiBWJFFpCCCFEBfTn2T8ZvXY07+1+j3xdvmFj+4kwdjOETQeghqUZnw4pGEL8Ze9l9ly4Wd7piiJIoSWEEEJUQF19umJvYc+Z5DMsPbPUsNHCBrxCDDa1q1eT4a0LhhCnLTtGdp5chVgRSKFlAuHh4QQFBdGqVSETHYUQQgjA0cqR8c3GA/DNkW9Izy3idjtJF2D/jwBM69sQd3tLYm5k8OVGWci0IpBCywTGjx9PdHQ0+/fvN3UqQgghKrDHAx7Hz96PpOwkFkYvNA5IuwbhbWDVVIg/hr2VOR/8s5DpDzsucPxqajlnLO4mhZYQQghRQZmrzZnYYiIAC08s5EbWXcs32HtCw36AAhtnANAjyJ1+TWuh1Sm8/sdR8rRyL0RTkkJLCCGEqMDCfMNo6tKUrPws5kTNMQ7o9g6ozeDcRrhQcOueGf0b4WhjTnRcGt9vv1DOGYs7SaElhBBCVGAqlYpXQl7B3cadZm7NjANq1oOWzxZ8v+Fd0OlwtbPknX5BAHy96SyXb2Ya7yfKhRRaQgghRAXX0qMlawavoX+9/oUHdH4dLGwhLgpOLANgcAsvOtZ3ISdfx7t/HUdRlPJLWOhJoSWEEEJUAuYa86IbbV0LbjoNsPn/QJuHSqXi/QGNsNCo2Xr6OmuPx5dPosKAFFpCCCFEJZGvy+fPs38yY9cM48a2L0HN+hA8HP5Z4LSuqy3jutQF4L2/o0nPyTfeT5QpKbSEEEKISiI+I573d7/P0rNLOZx42LDR0hbG74PQN8HcWr/5pa718XW2IT4tmy83nCnnjIUUWkIIIUQl4W3nzYD6AwAIjwo3DlBrjDZZmWt4b0AjAObvusjJuLQyzVEYkkJLCCGEqESeb/o8Zmoz9sbtZX98IQtfK0rBUg9LnoT8HAC6BrjRp7EHWp3C238eQ6eTifHlRQotIYQQohLxtPVkcP3BQEGvltHVhNpcWDEBTq2EQ4v0m9/tH0QNCw2HLqfw24HY8ky5WpNCSwghhKhkxjYdi7nanIMJB9kXv8+w0cwSOk0p+H7HTMjLBqCWgzWv9PAH4NO1p0jNzCvPlKstKbSEEEKISsajhgeP+T8GFNGr1WIk2HvDrTg4uEC/eVR7P/zdbUnOzOO/G2VifHmQQksIIYSohJ5r8hydvDoxofkE40YzS+j8asH3kV9Bfi4A5ho10/sXTIxfvOcSp+NvlVe61ZYUWkIIIUQl5Gbjxjdh39DKoxUqlco4oNmTYOsBt67B0SX6zR3qu9C7UcHE+Pf+PiErxpcxKbSEEEKIqsjMEtq/XPD9zi9Bp9U3vd0vEAszNbvO35QV48uYFFpCCCFEJXYz6yazDszi3ch3jRtDnoHaHQoWMb2Dj7MN4zoXrBj/waqTZOdpjfcVpUIKLSGEEKISu5l9k/kn5rP83HIupl40bLS0hWdWQ9OhRouZvhhaH08HK66mZPHdtgvll3A1I4WWEEIIUYn5O/kT6h2KgsK84/MeeD9rCw3T+gYCMGfbOa6mZJVVitWaFFpCCCFEJfdc0+cA+PvC38RnFDLnKjcT9syBP8YYbH6kaS1a13EmO0/HR6tPlkeq1Y4UWkIIIUQlF+waTCuPVuTr8ll4YqFxQHYKbHgXjv8Bl/fqN6tUKmb0b4RaBauOxrH3ws3yS7qakEJLCCGEqALGNC7orVp2dhm3cu9aH8veE5oOK/h+92yDpiBPe4a39gUKJsbLfRBLlxRaQgghRBXQ3rM99R3rk5mfybKzy4wD2o0v+HpqJSTFGDS90sMfW0szjl1NZXnU1XLItvqQQksIIYSoAlQqFWOajOHJwCfp5tvNOMAtEOqHgaKDvd8ZNLnYWvJS13oAfLb2NFm5stxDaZFCSwghhKgiHqn7CG+2fhMfO5/CA273ah1eDFkpBk3PdqiDl6M18WnZ/LhDlnsoLVJoCSGEENVF3a7gFgS56XBokUGTlbmGN/o0BGDOtvMkpmWbIsMqRwqtUjRo0CCcnJx47LHHTJ2KEEKIaiwqMYrJWyYTlRhl2KBSQfuJEDQQ/Doa7de/aS2a+TiSmatl5voz5ZJrVSeFVimaNGkSixYtun+gEEIIUYb+PPcnmy5vYlF0IX+Tmg2HoQvBq4VRk0ql4p1HChYx/e1gLNHX0so61SpPCq1SFBoaip2dnanTEEIIUc09Hfg0AJsubyL2Vmyx9g2p7Uy/prVQFPhwdTSKIss9PIwKUWhdvXqVp556ipo1a2JtbU2TJk04cOBAqR1/+/bt9O/fH09PT1QqFcuXLy80Ljw8HD8/P6ysrGjTpg379u0rtRyEEEKI8lLfqT4dPDugU3REnIooPOj6GVjzBlw9aNT0Zu+GWGjURJ67yZbTiWWcbdVm8kIrOTmZDh06YG5uzpo1a4iOjmbmzJk4OTkVGh8ZGUleXp7R9ujoaBISEgrdJyMjg+DgYMLDw4vMIyIigilTpjB9+nQOHTpEcHAwvXr1IjHx3w9Ys2bNaNy4sdHj2rVrxXzVQgghRNka3nA4UDCMmJVfyH0Md86Cvd/C3u+NmnycbXimgx8AH646SZ5WV5apVmkmL7Q+/fRTfHx8mD9/Pq1bt6ZOnTr07NmTevXqGcXqdDrGjx/PiBEj0Gr/XePj9OnTdOvWjYULC7ntANCnTx8++OADBg0aVGQes2bNYuzYsTzzzDMEBQXx7bffYmNjw7x5/96gMyoqiuPHjxs9PD09H+IdEEIIIUpfR6+OeNl6kZabxtqYtcYBrcYWfD2xDDJuGDW/1LU+TjbmnL+ewZJ9l8s426rL5IXWX3/9RcuWLXn88cdxc3OjefPm/PDDD4XGqtVqVq9ezeHDhxk5ciQ6nY7z58/TrVs3Bg4cyOuvv16iHHJzczl48CBhYWEG5woLC2P37t0lOua9hIeHExQURKtWrUr92EIIIQSARq1haMBQAH499avxXCvvEPBsDtpco6UeAByszXmlhz8A/914lrRs49EkcX8mL7QuXLjAnDlzaNCgAevWrePFF19k4sSJRfZOeXp6snnzZnbu3MmIESPo1q0bYWFhzJkzp8Q53LhxA61Wi7u7u8F2d3d34uMLuQt6EcLCwnj88cdZvXo13t7eRRZp48ePJzo6mv3795c4ZyGEEOJ+BtcfjL+TP/3r9UenFDL8d7tX68B80BmvBj+8tS91XWuQlJHLN1vOl3G2VZOZqRPQ6XS0bNmSjz76CIDmzZtz/Phxvv32W0aNGlXoPr6+vixevJguXbpQt25d5s6di0qlKs+0C7Vx40ZTpyCEEELoOVo5svTRpUUHNB4M69+G1MtwZh007GvQbK5RM61PIGMXHWBeZAxPtfXF28mmjLOuWkzeo1WrVi2CgoIMtgUGBnL5ctHjwQkJCTz//PP079+fzMxMXnnllYfKwcXFBY1GYzSZPiEhAQ8Pj4c6thBCCFFhmVtD84KlINhf+LSdsEA32tZ1JjdfxxfrTpdjclWDyQutDh06cPq04Q/uzJkz1K5du9D4Gzdu0L17dwIDA1m2bBmbNm0iIiKCqVOnljgHCwsLQkJC2LRpk36bTqdj06ZNtGvXrsTHFUIIISqC7Pxslp9bzt/n/zZubDUGbFzAoynojIcXVSoVb/ct6BBZHnWNI7EpZZxt1WLyocNXXnmF9u3b89FHHzF06FD27dvH999/z/ffG19uqtPp6NOnD7Vr1yYiIgIzMzOCgoLYsGED3bp1w8vLq9DerfT0dM6dO6d/HhMTQ1RUFM7Ozvj6+gIwZcoURo0aRcuWLWndujVffvklGRkZPPPMM2X34oUQQohysP7Set6JfAfPGp70rdMXjVrzb6OTH7x6CjTmRe7fxNuBwc29WHb4Kh+uPknE820rxJSdSkGpAP7++2+lcePGiqWlpdKwYUPl+++/LzJ2/fr1SlZWltH2Q4cOKbGxsYXus2XLFgUweowaNcog7n//+5/i6+urWFhYKK1bt1b27NnzUK/rflJTUxVASU1NLdPzCCGEqN6y8rKUjr92VBovaKxsvby1RMe4mpyp+L+9Wqn9xkpl7fG4Us6wcinO32+Vosja+qaSlpaGg4MDqamp2NvbmzodIYQQVdjn+z9nUfQiuvp05etuXxsH6HQQsw3MrKB24dNmPl93ivAt56njUoP1r3TGXGPyGUgmUZy/39XzHRJCCCGqmcENBgOw/cp2rmdeNw7YPRsWD4QtHxZ5jHFd6lGzhgUxNzL4Za8sYvogpNASQgghqoF6jvVo5toMraJlxfkVxgGNBgEquLgDbha+ZpadlTmT/1nE9MuNZ0jNkkVM70cKLSGEEKKaGOI/BIClZ5YaL2Dq6AP1/7lDSiErxd82vJUP9VxrkJyZxzdbzxUZJwpIoSWEEEJUEz1r98TO3A4vOy9Sc1KNA0L+WSg86hfQFt5bZaZR81bfQADmR14kNimzrNKtEqTQEkIIIaoJG3Mb1gxZw489f8TJysk4wL831HCDjEQ4U8iNqP/RraEb7evVJDdfx+eyiOk9SaElhBBCVCMOlg5FN2rModmIgu8PFn7PYShYxPStvoGoVPDXkWtEySKmRZJCSwghhKiGrmde5/iN48YNLUYWfE27Cvk5Re7f2MuBQc29APhwVTSyWlThpNASQgghqpntV7bT448evBP5jnGBVLMevLir4GFmec/jvNYrAEszNfsvJrPuRMI9Y6srKbSEEEKIaqa5W3PM1GacSznHiZsnjAPcG8ED3GKnloM1YzvVBeCTNSfJzTe+V2J1J4WWEEIIUc3YWdjRzbcbACvOFbKm1m15WZBeyOKmdxgXWg8XWwsu3szk572XSjPNKkEKLSGEEKIaGlBvAABrLq4hV5trHBD1K3zhDxtn3PM4tpZmvPLPIqZfbTori5jeRQotIYQQohpqW6strtaupOaksuPKDuMAp9qQkwbRyyH33mtlDWvpQwM3W1Iy8/hmiyxieicptIQQQohqSKPW8EjdRwAKvyWPbztw8oPcdDi18p7HkkVMiyaFlhBCCFFNPVrvUQD2xu0lM++u4kilguDhBd9H/XLfY4UGuNKhfk1ytTo+k0VM9aTQEkIIIaqp+k71+bzz56x/bD025jbGAU2HFXy9sBXSrt3zWHcuYvr3kWscvpxc+glXQlJoCSGEENVY7zq9i14t3rkO+LYHFDgacd9jNfJ0YEgLbwA+XHVSFjFFCi0hhBBC/EOr0xpvbPbP8OGRJfAAhdPUngFYmas5cCmZdSfiSznDykcKLSGEEKKa23l1JyNWjeCrQ18ZNwYNgC5vwBO/PNAiph4OVjyvX8T0VLVfxFQKLSGEEKKay8nP4diNY6y8sNK4V8vKAbq+VXBrngf0fJd6uNhacvFmJj/tqd6LmEqhJYQQQlRznb07Y29hz/Ws6xxIOPDQx7O1NOPVngWLmH69+SypmdV3EVMptIQQQohqzlxjTo/aPQBYE7Om8KAz6yHiqYIrEB/A4yHe+LsXLGI6e8vZUsq08pFCSwghhBD0rdMXgPWX1hd+S54za+Dk33Dk/lcfguEipgt3XeLyzeq5iKkUWkIIIYQgxD0EV2tXbuXeIvJqpHFAk8cLvp5aCXnZD3TMLv6udGrgQq5Wx6frTpVitpWHFFpCCCGEQKPW0MuvF1DE8KFPW7D3Krj/4dn1D3TMOxcxXXU0joOXqt8iplJoCSGEEAKAR+o+QlefrvSq08u4Ua2GxkMKvj/2+wMfM7CWPY+HFCxi+tHq6reIqRRaQgghhACgkUsjvu72Nd19uxce0OSxgq9n1kF22gMf99WeAVibazh4KZmVR+NKIdPKQwotIYQQQjwYj6bg4g/anIK5Wg/I3d6KcV0K1uH6YFU06Tn5ZZVhhSOFlhBCCCEMxKbF8t2R70jJTjFsUKmg8WPg2Rws7Yp1zBe61KV2TRsS0nL4auOZ0ku2gpNCSwghhBAGpmybwuyo2Wy4vMG4sfNUeH4rBPYv1jGtzDXMeLQRAPMiL3Iq/sGHHiszKbSEEEIIYaBPnT4ArL6w2rhRrSnxcbsGuNGrkTtancK7y09Ui4nxUmgJIYQQwkAfv4JC62DCQRIzEwsPyk6Fs4X0eN3Hu/0bYW2uYd/FJJYduvowaVYKUmiZQHh4OEFBQbRq1crUqQghhBBGatnWoqlrUxQUNl3eZByQmQSfN4CfH4db8cU6tpejNRO7NwAKlnuo6vdBlELLBMaPH090dDT79+83dSpCCCFEoXrW7gnA+ouFLE5q4wy1mgJKwW15imlMxzrUd7PlZkZulV8xvtiFVl5eHmZmZhw/frws8hFCCCFEBXD7JtMHEw5yI+uGcUDQgIKv0SuKfWwLMzX/N6AxAL/svczu8zdLnGdFV+xCy9zcHF9fX7RabVnkI4QQQogKwNPWk8Y1G2NlZsWZpEKWY7hdaF2KhPQi5nHdQ7t6NRne2heAN5cdJSu3atYVJRo6fPvtt3nrrbdISkoq7XyEEEIIUUF82vlTtg/bTnuv9saNjr7g2QIUXYmGDwGm9W2Ih70Vl25mMmvD6YfMtmIqUaE1e/Zstm/fjqenJwEBAbRo0cLgIYQQQojKz9feFyszq6ID9MOHy0t0fHsrcz4aXDCEOHdnDFGxKSU6TkVmVpKdBg4cWMppCCGEEKKiUhSFzPxMapjXMGwIGgAbp8OlXZCVAtaOxT52t4buDGruxZ+Hr/L6H0f4e0JHLM1KvlZXRaNSqsNqYRVUWloaDg4OpKamYm9vb+p0hBBCCCMHEw7y3u738LDx4Pue3xsHHImAOp3A3rPE50jOyKXHf7dxIz2X5zrW4T+PBD1ExmWvOH+/S7y8Q0pKCj/++CPTpk3Tz9U6dOgQV69W/cXHhBBCiOrC1dqVmNQY9sXvIzk72TggeNhDFVkATjUs+GRwUwB+3BnDjrPXH+p4FUmJCq2jR4/i7+/Pp59+yhdffEFKSgoAy5YtY9q0aaWZnxBCCCFMyNfel4bODdEqWrbEbimz84QFufNU24KrEF/97QhJGblldq7yVKJCa8qUKYwePZqzZ89iZfXvJLm+ffuyffv2UktOCCGEEKZ3z8VLAU6tgkUDIeqXhzrP232DqO9mS+KtHF7/42iVuBdiiQqt/fv388ILLxht9/LyIj6+eEvxCyGEEKJiu7146d64vaTmpBoHJJ6EC1vg+NKHOo+1hYavnmiGhUbNxpMJzIu8+FDHqwhKVGhZWlqSlpZmtP3MmTO4uro+dFJCCCGEqDj8HPyo71iffCWf7VcKGbm6vczDhW2QbVwfFEcjTwfe6tsQKLgXYmVfNb5Ehdajjz7K+++/T15ewY0gVSoVly9f5o033mDIkCGlmqAQQgghTK+7b3cANl/ebNzo0gBc/EGXB2eLGF4shlHt/RjU3AutTuHlXw4Rl5r10Mc0lRIVWjNnziQ9PR03NzeysrLo0qUL9evXx87Ojg8//LC0cxRCCCGEifWo3YNB9QcxxL+IDpWG/Qq+nlr10OdSqVR8NKgJgbXsuZmRywuLD5KZm//QxzWFh1pHa+fOnRw9epT09HRatGhBWFhYaeZW5ck6WkIIIaqMKwfgx+5gYQevnwczy4c+5OWbmTwavpOUzDy6N3Tju6dDMNOUeGWqUlOcv98lKrSys7MNrjYUJSOFlhBCiCpDp4NZgZAeD08uhQal0/ly8FISI37YS06+juGtffhoUBNUKlWpHLukynzBUkdHRzp37sw777zD5s2bycqqvGOnQgghhHgwiqJw/MZx5kTNQavTGjaq1QWT4ut0ATOLUjtnSG1nvh7eHLUKft0Xy0erT1aqZR9K1KO1c+dOtm/fztatW9m1axf5+fm0bNmSLl26EBoaSo8ePcoi1ypHerSEEEJUJnm6PEIjQknLTWNB7wWEuIcYBigKlFFv0y97L/PWn8cAGN3ej+n9g0zWs1XmPVodO3bkrbfeYv369aSkpLBlyxbq16/PZ599Ru/evUuUtBBCCCEqNnO1OV28uwBFXH1YhoXPiDa+fDy4CSoVLNh1kVd/O0J2nvb+O5qYWUl3PHPmDFu3btU/cnJyeOSRRwgNDS3F9IQQQghRkXTz7cbfF/5m8+XNTG05tfBepVvxkHoFvFuW6rmHt/bFTK3izWXHWHb4KjE3M/juqRDc7CvuvPESDR16eXmRlZVFaGgooaGhdOnShaZNm5p8clplI0OHQgghKpvMvEw6R3QmR5vD0keX4u/kbxhwbiP8NARqNoAJB8okh51nb/DSzwdJy87H0cac9wc0pn/TWuVWh5T50KGrqyuZmZnEx8cTHx9PQkKCTIgXQgghqgEbcxva1WoHFDF86N0aNBZw8yxcP1MmOXRs4MKKlzsSVMuelMw8Jv56mOE/7GFfTFKFmyhfokIrKiqK+Ph43nzzTXJycnjrrbdwcXGhffv2vP3226WdoxBCCCEqkG6+3YAiCi0r+4IrDwFOrSyzHOq41GDFyx2Y1L0BFho1ey4kMfS73fT473Y+Xn2SFVFXiTx3g0s3M8oshwfxUAuWAty8eZOtW7eyYsUKfv31V3Q6HVptxZ+cVhHI0KEQQojKKCk7ia6/dcXBwoGVg1dib3HX37AD82HlZPBqCWM3lXk+V1OymL35LH8evkp2ns6gbWhLbz57LLhUz1fmC5YuW7ZMPwk+OjoaZ2dnOnbsqJ+vFRxcui+oqpJCSwghRGV18uZJGjg1wExdyHV1txJg5j9zt149DXYe5ZLTrew8NkQnsP9iMucT00nJyqV/U08mdG9Qqucp80LLzc2Nzp076wurJk2alDjZ6kwKLSGEEFXW913h2iF49H/QYqSpsylVxfn7XaLlHRITE0uUmBBCCCGqFkVR0Ck6NGqNYYN/74JC68y6KldoFUeJ19HSarUsX76ckydPAhAUFMSAAQPQaDT32VMIIYQQVcG3R75l6dmlvNn6Tbr7djdsbPIYOPlB/dK552FlVaJC69y5c/Tt25erV68SEBAAwMcff4yPjw+rVq2iXr16pZqkEEIIISqelJwU4jPi2X5lu3GhVbNewaOaK9HyDhMnTqRevXrExsZy6NAhDh06xOXLl6lTpw4TJ04s7RyFEEIIUQF19u4MwPYr29EpuvtEV08l6tHatm0be/bswdnZWb+tZs2afPLJJ3To0KHUkhNCCCFExdXSvSU2ZjbcyLrByZsnaeTSyDAgOxUOzIO4I/DY/DK9F2JFVaIeLUtLS27dumW0PT09HQsLi4dOSgghhBAVn4XGgvae7QHYdmWbcYDaDLZ8DCf+hOunyzm7iqFEhdYjjzzC888/z969e1EUBUVR2LNnD+PGjePRRx8t7RyFEEIIUUHdHj4stNCyqAF1Cto5s7Ycs6o4SlRoff3119SvX5/27dtjZWWFlZUVHTp0oH79+nz11VelnaMQQgghKqhO3p0AiL4ZzfXM68YB/r0Kvp5ZV45ZVRzFmqOl0+n4/PPP+euvv8jNzWXgwIGMGjUKlUpFYGAg9evXL6s8hRBCCFEBuVi70LN2T9xs3NAqhdyCz78XrJ4KsXsgMwlsnI1jqrBiFVoffvghM2bMICwsDGtra1avXo2DgwPz5s0rq/yEEEIIUcHNDJ1ZdKOjL7g1gsQTcG4jNB1afolVAMUaOly0aBHffPMN69atY/ny5fz999/8/PPP6HRySacQQgghiqAfPqx+87SKVWhdvnyZvn376p+HhYWhUqm4du1aqScmhBBCiMojT5vH3ri9XEi5YNzo3xvMbUBjWf6JmVixCq38/HysrKwMtpmbm5OXl1eqSQkhhBCicvls/2c8t/45lpxeYtzo3RJej4FBc8o/MRMr1hwtRVEYPXo0lpb/VqTZ2dmMGzeOGjVq6LctW7as9DIUQgghRIXX3rM9S04vYfuV7UxrPQ3VnYuTqjUFj2qoWIXWqFGjjLY99dRTpZaMEEIIISqnNrXaYKY242r6VS7fukxt+9qFB6ZdA3vP8k3OhIpVaM2fP7+s8qgSBg0axNatW+nevTt//PGHqdMRQgghyo2NuQ0hbiHsjd/Lzqs7jQut/Bz4rjNcPwVTToF9LdMkWs5KtGCpKNykSZNYtGiRqdMQQgghTKKDV8H9jiOvRho3mlmCuXXB9+c3l2NWpiWFVikKDQ3Fzs7O1GkIIYQQJnG70Nofv58cbY5xQP0eBV/PbSzHrEyrQhVan3zyCSqVismTJ5fqcbdv307//v3x9PREpVKxfPnyQuPCw8Px8/PDysqKNm3asG/fvlLNQwghhKjKGjg2wM3GjWxtNgcTDhoH1A8r+Hp+M2jzyzc5E6kwhdb+/fv57rvvaNq06T3jIiMjC11OIjo6moSEhEL3ycjIIDg4mPDw8CKPGxERwZQpU5g+fTqHDh0iODiYXr16kZiYqI9p1qwZjRs3NnrIOmJCCCEEqFQqprebTsQjEbSt1dY4wCsErBwgOwWuHSr3/EyhQhRa6enpPPnkk/zwww84OTkVGafT6Rg/fjwjRoxAq/33fkqnT5+mW7duLFy4sND9+vTpwwcffMCgQYOKPPasWbMYO3YszzzzDEFBQXz77bfY2NgY3F4oKiqK48ePGz08PYt39UR4eDhBQUG0atWqWPsJIYQQFV1n784E1QxCrSqkxNCYQb1uBd9Xk+HDClFojR8/nn79+hEWFnbPOLVazerVqzl8+DAjR45Ep9Nx/vx5unXrxsCBA3n99ddLdP7c3FwOHjxocH61Wk1YWBi7d+8u0THvZfz48URHR7N///5SP7YQQghRod0ePqwmhVaxlncoC0uWLOHQoUMPXHR4enqyefNmOnXqxIgRI9i9ezdhYWHMmVPy1WZv3LiBVqvF3d3dYLu7uzunTp164OOEhYVx5MgRMjIy8Pb25vfff6ddu3YlzksIIYSojPbE7eHv83/TwbMDfev2NWys1x1ajIQGPU2TXDkzaaEVGxvLpEmT2LBhg9Gtfe7F19eXxYsX06VLF+rWrcvcuXMNV6A1kY0bq0d1LoQQQtzLkcQj/HX+LzLyMowLLfta8Oj/TJOYCZh06PDgwYMkJibSokULzMzMMDMzY9u2bXz99deYmZkZzMO6U0JCAs8//zz9+/cnMzOTV1555aHycHFxQaPRGE2mT0hIwMPD46GOLYQQQlQ3Hb06AgU9W3m66n0/ZJMWWt27d+fYsWNERUXpHy1btuTJJ58kKioKjcb4vkg3btyge/fuBAYGsmzZMjZt2kRERARTp04tcR4WFhaEhISwadMm/TadTsemTZtk6E8IIYQopsCagThbOZORl0FUYpRxgE4HVw7AjpkF31dhJh06tLOzo3HjxgbbatSoQc2aNY22Q0Hx06dPH2rXrk1ERARmZmYEBQWxYcMGunXrhpeXV6G9W+np6Zw7d07/PCYmhqioKJydnfH19QVgypQpjBo1ipYtW9K6dWu+/PJLMjIyeOaZZ0r5VQshhBBVm1qlpr1ne1ZeWEnk1Uhaedx1lb0uHxYNgNz0gjlbns1Mkmd5MPlk+OJQq9V89NFHdOrUCQsLC/324OBgNm7ciKura6H7HThwgK5du+qfT5kyBSi4SfaCBQsAGDZsGNevX+fdd98lPj6eZs2asXbtWqMJ8kIIIYS4vw5eHQoKrWuRTA6ZbNhoZgF1usDpVQVXH1bhQkulKIpi6iSqq7S0NBwcHEhNTcXe3t7U6QghhBClJik7idCIUBQUNj++GVebuzpD9s+FVVPAtx08u9Y0SZZQcf5+V6oeLSGEEEJUDs5WzjSq2YgcXQ6JmYnGhdbt9bRi90F2GlhVzQ4HKbSEEEIIUSYW9FmApcay8Ean2uBcF5IuwMWd0LBv4XGVXIVYGV4IIYQQVU+RRdZtdf+ZP31hS9knYyJSaAkhhBCiTGXnZ5ORl2HcUDe04OuVqntLOim0hBBCCFFmvj70NR2XdGTpmaXGjfW6wrPrYUzVvbOKFFpCCCGEKDMOlg7kaHPYHbfbuNHSDnzbgKbqThmXQksIIYQQZaZtrbYAHEw4SK4218TZlD8ptIQQQghRZho4NcDZypms/CyOXD9iHJCZBH9Phh+6V8nb8UihJYQQQogyo1ap9b1au68VMnxoYQtHf4OrByDxRDlnV/ak0BJCCCFEmWrn2Q6APXF7jBvNLMCvQ8H356veMg9SaAkhhBCiTN3u0Tpx8wSpOanGAbeXeaiC62lV3Wn+QgghhKgQPGp4MKTBEOo41Ck84PbCpZd2Q142mFuVX3JlTAotIYQQQpS5Ge1nFN3oFgi2HpAeD7F7oW6XcsurrMnQoRBCCCFMS6WqssOHUmgJIYQQolwkZCSw4twKEjMTjRvrdS24ybS1c/knVoZk6FAIIYQQ5eK17a9xOPEw09tN5zH/xwwbmwyF4CdMk1gZkh4tIYQQQpSLdrUKlnkodD0tddUsSarmqxJCCCFEhXN7Pa298XvRKUWsAq/Ng+SL5ZdUGZNCSwghhBDlopFLI2zMbEjNSeVM8hnjgCsH4FM/WDyo3HMrK1JoCSGEEKJcmKvNCXEPAWBv3F7jABd/yMuCpAuQeqWcsysbUmgJIYQQoty09mgNwP74/caNVvbg2bzg+5gd5ZhV2ZFCSwghhBDlpnWtgkLrYMJBtDqtcUCdTgVfY7aXY1ZlRwotIYQQQpSbAKcAvujyBSsHrUSj1hgH1Olc8PXiDlCU8k2uDEihJYQQQohyo1Fr6OXXi5rWNQsP8GkLanNIjYXkmPJNrgxIoSWEEEKIisPCBrxbFXxfBYYPpdASQgghRLnK0ebw/dHvGbdxHLnaXOOAkFHQ7R3wbV/+yZUyuQWPEEIIIcqVhdqCn0/+TFJ2EsduHNMv+aBXhW7FIz1aQgghhChXKpVKv8zDvrh9Js6mbEmhJYQQQohy18qjYB7W3vhCFi4FyLgBx5dW+vW0pNASQgghRLlrU6sNAEevHyUrP8s44OAC+ONZ2Pdd+SZWyqTQEkIIIUS587Xzxd3GnTxdHlGJUcYBdboUfL24E3RF3IC6EpBCSwghhBDl7s55WoXejsezOVjYQVYyJBwv5+xKjxRaQgghhDCJ1rVaY2NmU/gSDxozqN2u4PtKvJ6WLO8ghBBCCJPoU6cP/er2w1xtXnhAnc5wdn3B7Xjav1y+yZUSKbSEEEIIYRKWGst7B/j9c4Ppi5GgzS/o5apkZOhQCCGEECaXmZdpvNGjCVg5Qu6tSjtPq/KVhkIIIYSoMo7fOM60HdOwNrPmt/6/GTaqNTBsMTjXAwcv0yT4kKTQEkIIIYTJuNm4cTHtIipUpOWmYW9hbxhQp7NpEislMnQohBBCCJNxs3HD184XBaXw9bQqOSm0hBBCCGFSLT1aAnAg/kDhAft/hJ8fh7gj5ZhV6ZBCSwghhBAm1dL9n0IroYhC6+zGgmUeLmwrx6xKhxRaQgghhDCpEPcQAKJvRpORl2Ec4Neh4OulyHLMqnRIoSWEEEIIk/K09cSzhidaRcuRxEKGB2vfLrR2g05bvsk9JLnqUAghhBAm169uP5JzknG2djZu9GhacN/DnFRIOAG1mpZ/giUkhZYQQgghTG5ii4lFN2rMwLcNnNtYMHxYiQotGToUQgghRMV3e/jw4k7T5lFMUmgJIYQQokLI0+Vx5PoRrqZfNW706whmVgWrxVciUmiVokGDBuHk5MRjjz1m6lSEEEKISmfGrhk8tfopVpxbYdzoFQJvXoahi8o/sYcghVYpmjRpEosWVa4PgBBCCFFRBLsGA0Wsp6XWgJllOWf08KTQKkWhoaHY2dmZOg0hhBCiUrq9cOnR60fJ1eYWHZibWU4ZPTyTF1pz5syhadOm2NvbY29vT7t27VizZk2pnmP79u30798fT09PVCoVy5cvLzQuPDwcPz8/rKysaNOmDfv27SvVPIQQQghRtDoOdXC2ciZHm8PxG8eNA5IuQHhb+Lo5KEr5J1gCJi+0vL29+eSTTzh48CAHDhygW7duDBgwgBMnThQaHxkZSV5entH26OhoEhISCt0nIyOD4OBgwsPDi8wjIiKCKVOmMH36dA4dOkRwcDC9evUiMTFRH9OsWTMaN25s9Lh27VoxX7UQQggh7qZSqfSrxBc6fGjnWVBspcfDzfPlnF3JmLzQ6t+/P3379qVBgwb4+/vz4YcfYmtry549e4xidTod48ePZ8SIEWi1/64Me/r0abp168bChQsLPUefPn344IMPGDRoUJF5zJo1i7Fjx/LMM88QFBTEt99+i42NDfPmzdPHREVFcfz4caOHp6fnQ7wDQgghhLhNf9/Dwm4wbW4F3gXtXKocyzyYvNC6k1arZcmSJWRkZNCuXTujdrVazerVqzl8+DAjR45Ep9Nx/vx5unXrxsCBA3n99ddLdN7c3FwOHjxIWFiYwbnCwsLYvXt3iV9PUcLDwwkKCqJVq1alfmwhhBCiMmvpUVBIRV2PIk9nPIL173paleO+hxViZfhjx47Rrl07srOzsbW15c8//yQoKKjQWE9PTzZv3kynTp0YMWIEu3fvJiwsjDlz5pT4/Ddu3ECr1eLu7m6w3d3dnVOnTj3wccLCwjhy5AgZGRl4e3vz+++/F1owjh8/nvHjx5OWloaDg0OJ8xZCCCGqmvqO9ZkSMoXmbs1RF9Yf5NcBtlOwQryigEpV7jkWR4UotAICAoiKiiI1NZU//viDUaNGsW3btiKLLV9fXxYvXkyXLl2oW7cuc+fORVUB3uiNGzeaOgUhhBCiUlOr1DzT+JmiA7xbg9oM0q5CyiVw8iu33EqiQgwdWlhYUL9+fUJCQvj4448JDg7mq6++KjI+ISGB559/nv79+5OZmckrr7zyUOd3cXFBo9EYTaZPSEjAw8PjoY4thBBCiFJkYQOeLQq+rwTDhxWi0LqbTqcjJyen0LYbN27QvXt3AgMDWbZsGZs2bSIiIoKpU6eW+HwWFhaEhISwadMmgxw2bdpU6NCfEEIIIcpOni6PNTFr+GTfJ2h1WuOARoOg2ZMVvjcLKsDQ4bRp0+jTpw++vr7cunWLX375ha1bt7Ju3TqjWJ1OR58+fahduzYRERGYmZkRFBTEhg0b6NatG15eXoX2bqWnp3Pu3Dn985iYGKKionB2dsbX1xeAKVOmMGrUKFq2bEnr1q358ssvycjI4Jln7tF9KYQQQohSp0bNe7vfIyMvg0H1BxHgHGAY0O4l0yRWAiYvtBITExk5ciRxcXE4ODjQtGlT1q1bR48ePYxi1Wo1H330EZ06dcLCwkK/PTg4mI0bN+Lq6lroOQ4cOEDXrl31z6dMmQLAqFGjWLBgAQDDhg3j+vXrvPvuu8THx9OsWTPWrl1rNEFeCCGEEGVLo9bQzLUZkdciOZhw0LjQqkRUilJJllatgm5fdZiamoq9vb2p0xFCCCEqjO+OfMfsqNn09uvN510+Nw7QaSHhOFjYQs165Zpbcf5+V8g5WkIIIYSo3lq4F0x4P5RwiEL7hNa9Bd91hv0/lnNmxSOFlhBCCCEqnMYujTFTmZGYlci1jEJudef9z6Lfl0t/YfHSJIWWEEIIISocazNrgmoWrKd5KOGQcYDvP6sCxB2FnPRyzKx4pNASQgghRIXU3K05AGdTzho3OniBgy8oWrhayH0RKwiTX3UohBBCCFGYp4Ke4qmgp/CoUcTi4b5t4dhluLwH6oaWa24PSnq0hBBCCFEhedTwKLrIgoJCC+DSrvJJqASk0BJCCCFE5XR7ntaVA6DNM20uRZChQyGEEEJUWDuu7ODXU78S7BrMC8EvGDa6NoRu74BPa0BlkvzuR3q0hBBCCFFhJWUnsePqDnZe3WncqFZD56lQpzNoKmbfkRRaQgghhKiwbi9cevzmcbLzs02cTfFJoSWEEEKICsvb1htXa1fydfmcuHnCOCA/B6L/go3vQQW8q6AUWkIIIYSosFQqlX49rUIXLlUUWDoGds6CpAvlnN39VcwBTWFAq9WSl1cxr6YQQlRc5ubmaDQaU6chxENr4d6C9ZfWcyixkELL3Aq8QgpuxXN5d7nfYPp+pNCqwBRFIT4+npSUFFOnIoSopBwdHfHw8EClqphXZAnxIG73aB1JPIJWp0WjvusfEL5t/y20mj9lggyLJoVWBXa7yHJzc8PGxkZ+UQohHpiiKGRmZpKYmAhArVq1TJyRECXn7+RPTaua1HOsR2puKs5WzoYBvu2A/xasEF/BSKFVQWm1Wn2RVbNmTVOnI4SohKytrQFITEzEzc1NhhFFpWWmNmPz0M2oVUVMLb+9jtbNc5CeCLZu5Zrfvchk+Arq9pwsGxsbE2cihKjMbv8OkXmeorIrssgCsHYCt6CC7ytYr5YUWhWcDBcKIR6G/A4RVU1qTmrhDbfvexh3pPySeQAydCiEEEKICk+r0zJ05VDOJJ9h/ZD11LK9a95hh4kFD8fapkmwCNKjJYQQQogKT6PWYKmxBOBg4kHjACe/gkcF68WVQktUWH5+fnz55Zf3jFGpVCxfvrxc8inKg+Q5Y8YM3N3dK0S+QghRWTVzawYULPNQWUihJUpVaGgokydPNtq+YMECHB0dS/18cXFx9OnTp9SPW5pOnjzJe++9x3fffVcp8hVCiIoq2DUYgCPXiyi0Tq+FJU/Cvh/KMat7kzlaolLz8PAwdQr3df78eQAGDBjwUBOT8/LyMDc3L620hBCi0rldaJ1OPk1mXiY25nddmZ8cA6dWgjYPWo81QYbGpEerElEUhczc/HJ/KGVwk87Ro0czcOBAvvjiC2rVqkXNmjUZP3680SXot27dYvjw4dSoUQMvLy/Cw8MN2u8einvjjTfw9/fHxsaGunXr8s477xgc88iRI3Tt2hU7Ozvs7e0JCQnhwIED+vadO3fSqVMnrK2t8fHxYeLEiWRkZOjbExMT6d+/P9bW1tSpU4eff/75nq9zxowZ9O/fHwC1Wq0vtHQ6He+//z7e3t5YWlrSrFkz1q5dq9/v4sWLqFQqIiIi6NKlC1ZWVvz888/69+2jjz7C3d0dR0dH3n//ffLz83nttddwdnbG29ub+fPnP+BPQgghKg+PGh7UqlELnaLj2I1jxgE+rQu+XtkHOl35JlcE6dGqRLLytAS9u67czxv9fi9sLEr/o7JlyxZq1arFli1bOHfuHMOGDaNZs2aMHfvvv0I+//xz3nrrLd577z3WrVvHpEmT8Pf3p0ePHoUe087OjgULFuDp6cmxY8cYO3YsdnZ2vP766wA8+eSTNG/enDlz5qDRaIiKitL3Ep0/f57evXvzwQcfMG/ePK5fv87LL7/Myy+/rC9cRo8ezbVr19iyZQvm5uZMnDhRv/J2YaZOnYqfnx/PPPMMcXFx+u1fffUVM2fO5LvvvqN58+bMmzePRx99lBMnTtCgQQN93JtvvsnMmTNp3rw5VlZWbN26lc2bN+Pt7c327duJjIxkzJgx7Nq1i86dO7N3714iIiJ44YUX6NGjB97e3iX/AQkhRAXUzLUZcRlxRCVG0aZWG8NGj6ZgZg1ZyXDzLLgGmCbJO0ihJUzGycmJ2bNno9FoaNiwIf369WPTpk0GhVaHDh148803AfD39ycyMpL//ve/RRZa//nPf/Tf+/n5MXXqVJYsWaIvtC5fvsxrr71Gw4YNAQyKmo8//pgnn3xSP8esQYMGfP3113Tp0oU5c+Zw+fJl1qxZw759+2jVqhUAc+fOJTAwsMjXaGtrq5+bducw5xdffMEbb7zBE088AcCnn37Kli1b+PLLLw167SZPnszgwYMNjuns7MzXX3+NWq0mICCAzz77jMzMTN566y0Apk2bxieffMLOnTv1xxdCiKoi1CcUKzMrmrg0MW7UmBfcYPrSTojdK4WWKB5rcw3R7/cyyXnLQqNGjQxuCVKrVi2OHTPsCm7Xrp3R83td4RcREcHXX3/N+fPnSU9PJz8/H3t7e337lClTeO6551i8eDFhYWE8/vjj1KtXcKf3I0eOcPToUYPhQEVR0Ol0xMTEcObMGczMzAgJCdG3N2zYsNiT/NPS0rh27RodOnQw2N6hQweOHDGc4NmyZUuj/Rs1aoRa/e+ov7u7O40bN9Y/12g01KxZ8549bUIIUVn1rduXvnX7Fh3g0/rfQqvFyPJLrAgyR6sSUalU2FiYlfujOBO47e3tSU01XrU3JSUFBwcHg213T+xWqVToHmJMfffu3Tz55JP07duXlStXcvjwYd5++21yc3P1MTNmzODEiRP069ePzZs3ExQUxJ9//glAeno6L7zwAlFRUfrHkSNHOHv2rL4YK281atQw2lbY+1ba76UQQlRaPv8MJ8buM20e/5AeLVGqAgICWL9+vdH2Q4cO4e/vX+zj7dmzx+h5UUN1u3btonbt2rz99tv6bZcuXTKK8/f3x9/fn1deeYXhw4czf/58Bg0aRIsWLYiOjqZ+/fqFHr9hw4bk5+dz8OBB/dDh6dOnSUlJKdZrsre3x9PTk8jISLp06aLfHhkZSevWrYt1LCGEqI7ydfmcST6DpcaSeo53/UPYpzVoLMHaGfJzwczCNEn+Q3q0RKl68cUXOXPmDBMnTuTo0aOcPn2aWbNm8euvv/Lqq68W+3iRkZF89tlnnDlzhvDwcH7//XcmTZpUaGyDBg24fPkyS5Ys4fz583z99df63iqArKwsXn75ZbZu3cqlS5eIjIxk//79+sLtjTfeYNeuXbz88stERUVx9uxZVqxYwcsvvwwUFJG9e/fmhRdeYO/evRw8eJDnnnsOa2vrYr+u1157jU8//ZSIiAhOnz7Nm2++SVRUVJGvTQghxL++ifqGYSuHsSh6kXGjjTNMuwJj1pm8yALp0RKlrG7dumzfvp23336bsLAwcnNzadiwIb///ju9e/cu9vFeffVVDhw4wHvvvYe9vT2zZs2iV6/C56k9+uijvPLKK7z88svk5OTQr18/3nnnHWbMmAEUzF26efMmI0eOJCEhARcXFwYPHsx7770HQNOmTdm2bRtvv/02nTp1QlEU6tWrx7Bhw/TnmD9/Ps899xxdunTB3d2dDz74gHfeeafYr2vixImkpqby6quvkpiYSFBQEH/99ZfB5HwhhBCFu72eVlRiVOEBFaDAuk2llMUiSeKBpKWl4eDgQGpqqsGEbYDs7GxiYmKoU6cOVlZWJspQCFHZye8SURUlZyfTOaIzADuf2ImDpUPhgXnZYF76n/t7/f2+mwwdCiGEEKJScbJyws/eDyjidjwZN2BOR/i8XsEq8SYkhZYQQgghKp173vfQ2hnSrkBuOsQfLefMDEmhJYQQQohKJ9jtn0IrsZBCS60G73+u4jbxMg8yGV4IIYQQlU4z12YAHL1xlHxdPmbqu0qatuMg+Anw61j+yd1BCi0hhBBCVDr1HOsxqcWkwm/FA1CvW/kmVAQptIQQQghR6ahVap5r8pyp07gvmaMlhBBCCFFGpNASQgghRKWUnZ/N+ovrmRM1x9SpFEkKLSHuMnr0aAYOHGiSc1+8eBGVSkVUVJRJzl9WQkNDmTx5sqnTEEJUMfm6fKZum8o3R77heuZ1U6dTKJmjJcrE7t276dixI71792bVqlWmTqdYvvrqK8rjhgmjR48mJSWF5cuX67f5+PgQFxeHi4tLmZ+/PC1btgxzc3NTpyGEqGJsLWxp4NSAM8lnOHL9CGG1w0ydkhHp0RJlYu7cuUyYMIHt27dz7dq1cjlnbm5uqRzHwcEBR0fHUjlWcWk0Gjw8PDAzq1r/BnJ2dsbOzs7UaQghqqDbyzwUed9DE5NCS5S69PR0IiIiePHFF+nXrx8LFiwwaN+6dSsqlYpVq1bRtGlTrKysaNu2LcePH9fHLFiwAEdHR5YvX06DBg2wsrKiV69exMbG6mNmzJhBs2bN+PHHHw3u43b58mUGDBiAra0t9vb2DB06lISEBABOnTqFjY0Nv/zyi/44v/32G9bW1kRHRwPGQ4ehoaFMmDCByZMn4+TkhLu7Oz/88AMZGRk888wz2NnZUb9+fdasWaPfR6vVMmbMGOrUqYO1tTUBAQF89dVXBrkvXLiQFStWoFKpUKlUbN26tdChw23bttG6dWssLS2pVasWb775Jvn5+Qb5TZw4kddffx1nZ2c8PDz0N9IujpSUFF544QXc3d2xsrKicePGrFy5Ut++dOlSGjVqhKWlJX5+fsycOdNg/2+++Ub/s3J3d+exxx4zyPHOoUM/Pz8++ugjnn32Wezs7PD19eX77783OF5sbCxDhw7F0dERZ2dnBgwYwMWLF4v9uoQQVdvthUujrkeZNpEiSKFVGeVmFP3Iyy5GbNb9Y0vgt99+o2HDhgQEBPDUU08xb968QofiXnvtNWbOnMn+/ftxdXWlf//+5OX9e0+qzMxMPvzwQxYtWkRkZCQpKSk88cQTBsc4d+4cS5cuZdmyZURFRaHT6RgwYABJSUls27aNDRs2cOHCBYYNGwZAw4YN+eKLL3jppZe4fPkyV65cYdy4cXz66acEBQUV+ZoWLlyIi4sL+/btY8KECbz44os8/vjjtG/fnkOHDtGzZ0+efvppMjMzAdDpdHh7e/P7778THR3Nu+++y1tvvcVvv/0GwNSpUxk6dCi9e/cmLi6OuLg42rdvb3Teq1ev0rdvX1q1asWRI0eYM2cOc+fO5YMPPjDKr0aNGuzdu5fPPvuM999/nw0bNujbR48eTWhoaJGvT6fT0adPHyIjI/npp5+Ijo7mk08+QaPRAHDw4EGGDh3KE088wbFjx5gxYwbvvPOOvog+cOAAEydO5P333+f06dOsXbuWzp07F3k+gJkzZ9KyZUsOHz7MSy+9xIsvvsjp06cByMvLo1evXtjZ2bFjxw4iIyOxtbWld+/epdZzKYSoGm73aEXfjCZXWwF/PyjCZFJTUxVASU1NNWrLyspSoqOjlaysLOMdp9sX/fjpMcPYDzyKjp3X1zD20zrGMSXQvn175csvv1QURVHy8vIUFxcXZcuWLfr2LVu2KICyZMkS/babN28q1tbWSkREhKIoijJ//nwFUPbs2aOPOXnypAIoe/fuLXgbpk9XzM3NlcTERH3M+vXrFY1Go1y+fFm/7cSJEwqg7Nu3T7+tX79+SqdOnZTu3bsrPXv2VHQ6nb5t1KhRyoABA/TPu3TponTs2FH/PD8/X6lRo4by9NNP67fFxcUpgLJ79+4i35fx48crQ4YMKfI8iqIoMTExCqAcPnxYURRFeeutt5SAgACD/MLDwxVbW1tFq9UWmp+iKEqrVq2UN954Q//8zTffNMj3buvWrVPUarVy+vTpQttHjBih9OjRw2Dba6+9pgQFBSmKoihLly5V7O3tlbS0tEL379KlizJp0iT989q1aytPPfWU/rlOp1Pc3NyUOXPmKIqiKIsXLzZ63Tk5OYq1tbWybt26Il+HMHbP3yVCVAE6nU7p9GsnpfGCxsqRxCPlcs57/f2+m/RoiVJ1+vRp9u3bx/DhwwEwMzNj2LBhzJ071yi2Xbt2+u+dnZ0JCAjg5MmT+m1mZma0atVK/7xhw4Y4OjoaxNSuXRtXV1f985MnT+Lj44OPj49+W1BQkNF+8+bN4+jRoxw6dIgFCxagUqnu+bqaNm2q/16j0VCzZk2aNPl3NWJ3d3cAEhMT9dvCw8MJCQnB1dUVW1tbvv/+ey5fvnzP89zt5MmTtGvXziC/Dh06kJ6ezpUrVwrND6BWrVoGuXz88ccsWrSoyPNERUXh7e2Nv79/kXl06NDBYFuHDh04e/YsWq2WHj16ULt2berWrcvTTz/Nzz//rO/dK8qdOatUKjw8PPQ5HzlyhHPnzmFnZ4etrS22trY4OzuTnZ3N+fPn73lcIUT1olKpaOJa8Pv45M2T94kuf1Vrxm118dY9JperNIbPXzt3j9i76uzJx0qe0z/mzp1Lfn4+np6e+m2KomBpacns2bNxcHB46HPcqUaNGiXa78iRI2RkZKBWq4mLi6NWrVr3jL/7ijmVSmWw7XYhpNPpAFiyZAlTp05l5syZtGvXDjs7Oz7//HP27t1bonzvp7D8bufyIKytrR/q/HZ2dhw6dIitW7eyfv163n33XWbMmMH+/fuLvLDgXjmnp6cTEhLCzz//bLTfnYW1EEIATAmZwttt3qZWjXv/LjcFKbQqI4tiFBdlFVuI/Px8Fi1axMyZM+nZs6dB28CBA/n1118ZN26cftuePXvw9fUFIDk5mTNnzhAYGGhwvAMHDtC6dcEd2E+fPk1KSopBzN0CAwOJjY0lNjZW36sVHR1NSkqKfg5WUlISo0eP5u233yYuLo4nn3ySQ4cOPXSxcafIyEjat2/PSy+9pN92d0+MhYUFWq32nscJDAxk6dKlKIqiL+YiIyOxs7PD29u71PJt2rQpV65c4cyZM4X2agUGBhIZGWmwLTIyEn9/f/08LjMzM8LCwggLC2P69Ok4OjqyefNmBg8eXOx8WrRoQUREBG5ubtjb25fsRQkhqo16jvVMnUKRZOhQlJqVK1eSnJzMmDFjaNy4scFjyJAhRsOH77//Pps2beL48eOMHj0aFxcXg6v9zM3NmTBhAnv37uXgwYOMHj2atm3b6guvwoSFhdGkSRN98bRv3z5GjhxJly5daNmyJQDjxo3Dx8eH//znP8yaNQutVsvUqVNL9b1o0KABBw4cYN26dZw5c4Z33nmH/fv3G8T4+flx9OhRTp8+zY0bNwwuBLjtpZdeIjY2lgkTJnDq1ClWrFjB9OnTmTJlCmr1g//vO23aNEaOHFlke5cuXejcuTNDhgxhw4YNxMTEsGbNGtauXQvAq6++yqZNm/i///s/zpw5w8KFC5k9e7b+fVu5ciVff/01UVFRXLp0iUWLFqHT6QgICHjgHO/05JNP4uLiwoABA9ixYwcxMTFs3bqViRMnGgyZCiFERSeFlig1c+fOJSwsrNDhwSFDhnDgwAGOHj2q3/bJJ58wadIkQkJCiI+P5++//8bCwkLfbmNjwxtvvMGIESPo0KEDtra2RERE3DMHlUrFihUrcHJyonPnzoSFhVG3bl39fosWLWL16tUsXrwYMzMzatSowU8//cQPP/xgsDzDw3rhhRcYPHgww4YNo02bNty8edOgdwtg7NixBAQE0LJlS1xdXY16jAC8vLxYvXo1+/btIzg4mHHjxjFmzBj+85//FCufuLi4+84PW7p0Ka1atWL48OEEBQXx+uuv63vcWrRowW+//caSJUto3Lgx7777Lu+//z6jR48GwNHRkWXLltGtWzcCAwP59ttv+fXXX2nUqFGx8rzNxsaG7du34+vry+DBgwkMDGTMmDFkZ2dLD5cQolDLzy1n0uZJ7L6229SpGFApSjksgS0KlZaWhoODA6mpqUZ/PLKzs4mJiTFYH6qq2Lp1K127diU5ObnI+TsLFixg8uTJpKSklGtuQlQ1Vfl3iRB3ejfyXf489ydjm4xlYouJZXque/39vpv0aAkhhBCi0mvqWnAl89HrR+8TWb6k0BJCCCFEpXe70Dp24xha3b0vNCpPUmiJchcaGoqiKPe8n+DtGy4LIYQQD6KeQz1szGzIzM/kfGrFWW9PCi0hhBBCVHoatYYmLgULlx67/vDrQpYWKbSEEEIIUSXcXiH+6I2KM09LCi0hhBBCVAlNXZpiobYgT2u8LqGpyMrwQgghhKgSOnp1ZM+IPZhrzO8fXE6k0BJCCCFElVCRCqzbZOhQCCGEEFWOTtGZOgVACi0hjIwePdrgnovl6eLFi6hUKqKiokxy/rISGhrK5MmTTZ1GsWzduhWVSlXsZUZyc3OpX78+u3btKtVzzZgxg2bNmumfv/nmm0yYMKFYuQlRHey+tpvBfw3mlS2vmDoVQAotUUZ2796NRqOhX79+pk6l2L766isWLFhQ5ucprKDz8fEhLi6Oxo0bl/n5y9OyZcv4v//7vzI/j06n44033sDT0xNra2uaNm3KihUryvy8d/r222+pU6cO7du3L9XjTp06lU2bNhk8X7hwIRcuXCjV8whR2dmY23A2+SxR16OoCHcZlEJLlIm5c+cyYcIEtm/fzrVr18rlnLm5uaVyHAcHh3suplqWNBoNHh4emJlVremTzs7O2NnZlfl5fvrpJ/773/8ya9YsTp48yaxZs6hRo0aZn/c2RVGYPXs2Y8aMKfVj29raUrNmTf1zFxcXevXqxZw5c0r9XEJUZoHOgZirzUnKTuJK+hVTpyOFlih96enpRERE8OKLL9KvXz+j3qHbwySrVq2iadOmWFlZ0bZtW44fP66PWbBgAY6OjixfvpwGDRpgZWVFr169iI2N1cfcHkr58ccfDW6Ye/nyZQYMGICtrS329vYMHTqUhIQEAE6dOoWNjQ2//PKL/ji//fYb1tbWREdHA8Y9TaGhoUyYMIHJkyfj5OSEu7s7P/zwAxkZGTzzzDPY2dlRv3591qxZo99Hq9UyZswY6tSpg7W1NQEBAXz11VcGuS9cuJAVK1agUqlQqVRs3bq10KHDbdu20bp1aywtLalVqxZvvvkm+fn5BvlNnDiR119/HWdnZzw8PJgxY0axf24pKSm88MILuLu7Y2VlRePGjVm5cqW+fenSpTRq1AhLS0v8/PyYOXOmwf7ffPON/mfl7u7OY489ZpDjnUOHfn5+fPTRRzz77LPY2dnh6+vL999/b3C82NhYhg4diqOjI87OzgwYMICLFy/e8zWo1WpcXV154okn8PPzIywsjLCwsGK/F4W5/Zlct24dgYGB2Nra0rt3b+Li4vQxBw8e5Pz58wY9ubd/pkuWLKF9+/b693bbtm1G5zh48CAtW7bExsaG9u3bc/r0aX3b3UOHAP3792fJkiWl8vqEqCosNBYEOgcCFeO+h1JoVUKZeZlFPnK0OQ8cm52ffd/Ykvjtt99o2LAhAQEBPPXUU8ybN6/Q7tvXXnuNmTNnsn//flxdXenfvz95ef+ufZKZmcmHH37IokWLiIyMJCUlhSeeeMLgGOfOnWPp0qUsW7aMqKgodDodAwYMICkpiW3btrFhwwYuXLjAsGHDAGjYsCFffPEFL730EpcvX+bKlSuMGzeOTz/9lKCgoCJf08KFC3FxcWHfvn1MmDCBF198kccff5z27dtz6NAhevbsydNPP01mZsF7ptPp8Pb25vfffyc6Opp3332Xt956i99++w0oGPYZOnSo/g91XFxcoUNNV69epW/fvrRq1YojR44wZ84c5s6dywcffGCUX40aNdi7dy+fffYZ77//Phs2bNC3jx49mtDQ0CJfn06no0+fPkRGRvLTTz8RHR3NJ598gkajAQoKgKFDh/LEE09w7NgxZsyYwTvvvKMvog8cOMDEiRN5//33OX36NGvXrqVz585Fng9g5syZtGzZksOHD/PSSy/x4osv6guLvLw8evXqhZ2dHTt27CAyMlJf2Nyr57J79+6kpqbyzjvv3PPcJZWZmckXX3zB4sWL2b59O5cvX2bq1Kn69h07duDv719o791rr73Gq6++yuHDh2nXrh39+/fn5s2bBjFvv/02M2fO5MCBA5iZmfHss8/eM5/WrVtz5cqV+xagQlQ3d9730OQUYTKpqakKoKSmphq1ZWVlKdHR0UpWVpZRW+MFjYt8vLjhRYPYVj+1KjJ29JrRBrGdfu1kFFMS7du3V7788ktFURQlLy9PcXFxUbZs2aJv37JliwIoS5Ys0W+7efOmYm1trURERCiKoijz589XAGXPnj36mJMnTyqAsnfvXkVRFGX69OmKubm5kpiYqI9Zv369otFolMuXL+u3nThxQgGUffv26bf169dP6dSpk9K9e3elZ8+eik6n07eNGjVKGTBggP55ly5dlI4dO+qf5+fnKzVq1FCefvpp/ba4uDgFUHbv3l3k+zJ+/HhlyJAhRZ5HURQlJiZGAZTDhw8riqIob731lhIQEGCQX3h4uGJra6totdpC81MURWnVqpXyxhtv6J+/+eabBvnebd26dYparVZOnz5daPuIESOUHj16GGx77bXXlKCgIEVRFGXp0qWKvb29kpaWVuj+Xbp0USZNmqR/Xrt2beWpp57SP9fpdIqbm5syZ84cRVEUZfHixUavOycnR7G2tlbWrVtX6DkyMjKURo0aKWPHjlXatGmjvPrqqwb729nZKb///nuR78Hdbn9Ok5OTFUX59zN57tw5fUx4eLji7u6ufz5p0iSlW7duBse5/TP95JNP9Nvy8vIUb29v5dNPPzU418aNG/Uxq1atUgD974Dp06crwcHBBse+/Ttk69atRb6Oe/0uEaKqWnV+ldJ4QWNl+MrhZXL8e/39vpv0aIlSdfr0afbt28fw4cMBMDMzY9iwYcydO9cotl27dvrvnZ2dCQgI4OTJk/ptZmZmtGrVSv+8YcOGODo6GsTUrl0bV1dX/fOTJ0/i4+ODj4+PfltQUJDRfvPmzePo0aMcOnSIBQsWoFKp7vm6mjZtqv9eo9FQs2ZNmjRpot/m7u4OQGJion5beHg4ISEhuLq6Ymtry/fff8/ly5fveZ67nTx5knbt2hnk16FDB9LT07ly5d+5B3fmB1CrVi2DXD7++GMWLVpU5HmioqLw9vbG39+/yDw6dOhgsK1Dhw6cPXsWrVZLjx49qF27NnXr1uXpp5/m559/1vfuFeXOnFUqFR4eHvqcjxw5wrlz57Czs8PW1hZbW1ucnZ3Jzs7m/PnCbxa7YMECUlJSCA8PZ82aNWzYsIFnnnmG/Px8Ll68SHp6utFruK1Pnz768zRq1KjInG1sbKhXr57++d3vc1ZWln4I+253ft7NzMxo2bKlwWfy7vekVq1agOFn6m7W1tYA932vhahubvdonUw6aTTSU96q1ozbamLviL1FtmnUGoPnW4duLTJWrTKss9cOWftQeUHBJPj8/Hw8PT312xRFwdLSktmzZ+Pg4PDQ57hTSSc6HzlyhIyMDNRqNXFxcfo/akUxNzdcBE+lUhlsu10I6XQF67YsWbKEqVOnMnPmTNq1a4ednR2ff/45e/cW/bN7GIXldzuXB3H7D3ZJ2dnZcejQIbZu3cr69et59913mTFjBvv37y/ywoJ75Zyenk5ISAg///yz0X53FtZ3Onr0KI0aNcLc3BwnJyc2bNhAp06dGDRoEA0aNKB3795F/px//PFHsrKyCs3rfjkrdwyLu7i4cOxYyYcq7vWZKkxSUhJQ9HsiRHXlZetFsGswPnY+pOemY2ltabJcpNCqhGzMbUweW5j8/HwWLVrEzJkz6dmzp0HbwIED+fXXXxk3bpx+2549e/D19QUgOTmZM2fOEBgYaHC8AwcO0Lp1a6CgtywlJcUg5m6BgYHExsYSGxur79WKjo4mJSVFPwcrKSmJ0aNH8/bbbxMXF8eTTz7JoUOHHrrYuFNkZCTt27fnpZde0m+7uyfGwsICrVZ7z+MEBgaydOlSFEXR/+GNjIzEzs4Ob2/vUsu3adOmXLlyhTNnzhTaqxUYGEhkZKTBtsjISPz9/fXzuMzMzPSTz6dPn46joyObN29m8ODBxc6nRYsWRERE4Obmhr29/QPt4+XlxZ9//smtW7ews7PDzc2NjRs30qlTJ1auXMnBgwfvuW9paN68OXPmzDH4ed22Z88e/by1/Px8Dh48yMsvv/xQ5zt+/Djm5ub37IUTojpSqVT81PcnU6cByGR4UYpWrlxJcnIyY8aMoXHjxgaPIUOGGA0fvv/++2zatInjx48zevRoXFxcDK72Mzc3Z8KECezdu5eDBw8yevRo2rZtqy+8ChMWFkaTJk30xdO+ffsYOXIkXbp0oWXLlgCMGzcOHx8f/vOf/zBr1iy0Wq3BhObS0KBBAw4cOMC6des4c+YM77zzDvv37zeI8fPz4+jRo5w+fZobN24YXAhw20svvURsbCwTJkzg1KlTrFixgunTpzNlyhTU6gf/33fatGmMHDmyyPYuXbrQuXNnhgwZwoYNG4iJiWHNmjWsXVvQy/nqq6+yadMm/u///o8zZ86wcOFCZs+erX/fVq5cyddff01UVBSXLl1i0aJF6HQ6AgICHjjHOz355JO4uLgwYMAAduzYQUxMDFu3bmXixIkGQ6Z3GjNmDFqtlkcffZRdu3Zx+vRp1q1bR3p6OjY2NoUOX5e2rl27kp6ezokTJ4zawsPD+fPPPzl16hTjx48nOTn5vpPd72fHjh106tSpVP+RIIQoXVJoiVIzd+5cwsLCCh0eHDJkCAcOHODo0X8vtf3kk0+YNGkSISEhxMfH8/fff2NhYaFvt7Gx4Y033mDEiBF06NABW1tbIiIi7pmDSqVixYoVODk50blzZ8LCwqhbt65+v0WLFrF69WoWL16MmZkZNWrU4KeffuKHH34wWJ7hYb3wwgsMHjyYYcOG0aZNG27evGnQuwUwduxYAgICaNmyJa6urkY9RlDQ07J69Wr27dtHcHAw48aNY8yYMfznP/8pVj5xcXH3nR+2dOlSWrVqxfDhwwkKCuL111/X97i1aNGC3377jSVLltC4cWPeffdd3n//fUaPHg2Ao6Mjy5Yto1u3bgQGBvLtt9/y66+/lrinxcbGhu3bt+Pr68vgwYMJDAxkzJgxZGdnF9nD5enpyb59+3BxcWHw4ME0b96cRYsWsWjRIlatWsX333/PrFmzSpTPg6pZsyaDBg0qdMjzk08+4ZNPPiE4OJidO3fy119/4eLi8lDnW7JkCWPHjn2oYwghypZKUSrAsqnVVFpaGg4ODqSmphr98cjOziYmJsZgfaiqYuvWrXTt2pXk5OQi5+8sWLCAyZMnF/v2J0KY2tGjR+nRowfnz5/H1taWixcvUqdOHQ4fPmy0DtbDWLNmDa+++ipHjx695wK3Vfl3iRCmcq+/33eTHq1SNGjQIJycnAwWahRCVC9Nmzbl008/JSYmpkzPk5GRwfz586vcXQSEqGrk/9BSNGnSJJ599lkWLlxo6lSEECZ0e0i1LMk/6ISoHKRHqxSFhoaWy/3cKrvQ0FAURbnn/QRHjx4tw4aiSvDz80NRlFIdNhRCVB4mL7Q+/vhjWrVqpb8ce+DAgQb39yoN27dvp3///nh6eqJSqVi+fHmhceHh4fj5+WFlZUWbNm3Yt29fqeYhhBBCiOrF5IXWtm3bGD9+PHv27GHDhg3k5eXRs2dPMjIyCo2PjIws9DL46Oho/Y2D75aRkUFwcDDh4eFF5hEREcGUKVOYPn06hw4dIjg4mF69ehmsytysWTOjZQsaN27MtWvXivmqH5xcqyCEeBjyO0QI0zL5HK3b6/TctmDBAtzc3Dh48KDRTWl1Oh3jx4+nQYMGLFmyRL9Q4unTp+nWrRtTpkzh9ddfNzpHnz596NOnzz3zmDVrFmPHjuWZZ54B4Ntvv2XVqlXMmzePN998Eyi4TUl5ub1CdGZmpqyRI4Qosdu357nXivdCiLJj8kLrbqmpqUDBve/uplarWb16NZ07d2bkyJEsXryYmJgYunXrxsCBAwstsh5Ebm4uBw8eZNq0aQbnCgsLY/fu3SV7IfcQHh5OeHj4PVcF12g0ODo66nvUbGxs7ns/PiGEuE1RFDIzM0lMTMTR0VH/D1MhRPmqUIWWTqdj8uTJdOjQgcaNGxca4+npyebNm+nUqRMjRoxg9+7dhIWFMWfOnBKf98aNG2i1Wv2NgW9zd3fn1KlTD3ycsLAw/T30vL29+f333w1uJHvb+PHjGT9+vH4djqJ4eHgA976prBBC3Iujo6P+d4kQovxVqEJr/PjxHD9+nJ07d94zztfXl8WLF9OlSxfq1q3L3LlzK0Rvz8aNG0v1eCqVilq1auHm5lbovDQhhLgXc3Nz6ckSwsQqTKH18ssvs3LlSrZv337fm+UmJCTw/PPP079/f/bv388rr7zC//73vxKf28XFBY1GYzSZPiEhoUL8S1Cj0cgvSyGEEKISMvlVh4qi8PLLL/Pnn3+yefNm6tSpc8/4Gzdu0L17dwIDA1m2bBmbNm0iIiLioW4KbGFhQUhICJs2bdJv0+l0bNq0qdChPyGEEEKIB2HyHq3x48fzyy+/sGLFCuzs7IiPjwfAwcHB6Go7nU5Hnz59qF27NhEREZiZmREUFMSGDRvo1q0bXl5evPLKK0bnSE9P59y5c/rnMTExREVF4ezsjK+vLwBTpkxh1KhRtGzZktatW/Pll1+SkZGhvwpRCCGEEKK4TH5T6aLmVs2fP7/Q21hs2LCBTp06Gd0c9fDhw7i6uhY67Hj7JsZ3GzVqFAsWLNA/nz17Np9//jnx8fE0a9aMr7/+mjZt2hTvBRVDcW5KKYQQQoiKoTh/v01eaFVnqampODo6EhsbK4WWEEIIUUmkpaXh4+NDSkrKPVcPgAowdFid3bp1CwAfHx8TZyKEEEKI4rp169Z9Cy3p0TIhnU7HtWvXsLOzMxpCbdWqFfv373+o45fkGMXZ50Fj7xd3r/bb/2qoCr1+pfEzrSjnfdhjymezYpHP5sPtL5/NslNRP5uKonDr1i08PT1Rq+99XaH0aJmQWq0ucikLjUbz0P+DlOQYxdnnQWPvF/cgx7G3t6/0vzBK42daUc77sMeUz2bFIp/Nh9tfPptlpyJ/Nu/Xk3WbyZd3EIUbP368SY5RnH0eNPZ+caXxWisDU73Osjjvwx5TPpsVi3w2H25/+WyWnarw2ZShQ1GhyZWZoqKSz6aoqOSzWbFIj5ao0CwtLZk+fTqWlpamTkUIA/LZFBWVfDYrFunREkIIIYQoI9KjJYQQQghRRqTQEkIIIYQoI1JoCSGEEEKUESm0hBBCCCHKiBRaQgghhBBlRAotUWmtXLmSgIAAGjRowI8//mjqdITQGzRoEE5OTjz22GOmTkUIA7GxsYSGhhIUFETTpk35/fffTZ1SlSfLO4hKKT8/n6CgILZs2YKDgwMhISHs2rWLmjVrmjo1Idi6dSu3bt1i4cKF/PHHH6ZORwi9uLg4EhISaNasGfHx8YSEhHDmzBlq1Khh6tSqLOnREpXSvn37aNSoEV5eXtja2tKnTx/Wr19v6rSEACA0NBQ7OztTpyGEkVq1atGsWTMAPDw8cHFxISkpybRJVXFSaAmT2L59O/3798fT0xOVSsXy5cuNYsLDw/Hz88PKyoo2bdqwb98+fdu1a9fw8vLSP/fy8uLq1avlkbqo4h72sylEWSrNz+fBgwfRarX4+PiUcdbVmxRawiQyMjIIDg4mPDy80PaIiAimTJnC9OnTOXToEMHBwfTq1YvExMRyzlRUN/LZFBVZaX0+k5KSGDlyJN9//315pF29KUKYGKD8+eefBttat26tjB8/Xv9cq9Uqnp6eyscff6woiqJERkYqAwcO1LdPmjRJ+fnnn8slX1F9lOSzeduWLVuUIUOGlEeaopoq6eczOztb6dSpk7Jo0aLySrVakx4tUeHk5uZy8OBBwsLC9NvUajVhYWHs3r0bgNatW3P8+HGuXr1Keno6a9asoVevXqZKWVQTD/LZFMJUHuTzqSgKo0ePplu3bjz99NOmSrVakUJLVDg3btxAq9Xi7u5usN3d3Z34+HgAzMzMmDlzJl27dqVZs2a8+uqrcsWhKHMP8tkECAsL4/HHH2f16tV4e3tLESbKxYN8PiMjI4mIiGD58uU0a/b/7d1vSFNrHAfw7+5yops1bKYTnEKGRQiasRA0TUWtUCOcYqGbWhAZJZFE74JeVUpSgfWipkRIQmT0RiubfxgkSSkKYTkMLGL0T0jLzdxzX1w4dO6EW17Pndd9P3BePM+enfPb4WF8ec4zloKUlBSMjo4GotygsSbQBRAtVXFxMYqLiwNdBpGfx48fB7oEokVlZGTA5/MFuoygwhUtWnEMBgPUajXcbres3+12IyYmJkBVEXFu0srG+bkyMWjRiqPRaJCWloaenh6pz+fzoaenB+np6QGsjIId5yatZJyfKxMfHVJAzMzMYGJiQmpPTk5ieHgYkZGRMJlMOHnyJKxWK7Zv3w6z2Yzm5mbMzs6iuro6gFVTMODcpJWM8/N/KNA/e6Tg5HA4BAC/w2q1SmOuXLkiTCaT0Gg0wmw2i6dPnwauYAoanJu0knF+/v/wvw6JiIiIFMI9WkREREQKYdAiIiIiUgiDFhEREZFCGLSIiIiIFMKgRURERKQQBi0iIiIihTBoERERESmEQYuIiIhIIQxaRERERAph0CIiIiJSCIMWEa06NpsNKpXK7/j5z3iJiP4LawJdABGREgoLC2G322V9UVFRsrbX64VGo/kvyyKiIMMVLSJalUJDQxETEyM7cnNzcezYMdTX18NgMKCgoAAAMDY2ht27d0On0yE6OhqVlZX4+PGjdK7Z2VlUVVVBp9PBaDSiqakJ2dnZqK+vl8aoVCp0dnbKatDr9WhtbZXaU1NTKCsrg16vR2RkJEpKSvDmzRvpdZvNhn379qGxsRFGoxHr169HXV0d5ufnpTEejwenT59GXFwcQkNDkZiYiBs3bkAIgcTERDQ2NspqGB4e5moeUQAxaBFRUGlra4NGo4HT6cS1a9cwPT2NnJwcpKamYmhoCF1dXXC73SgrK5Pe09DQgL6+Pty/fx8PHz5Eb28vnj9//lvXnZ+fR0FBASIiIjAwMACn0wmdTofCwkJ4vV5pnMPhgMvlgsPhQFtbG1pbW2VhraqqCu3t7bh8+TJevnyJ69evQ6fTQaVSoaamxm8Vz263Y+fOnUhMTFzaDSOif0cQEa0yVqtVqNVqodVqpaO0tFRkZWWJ1NRU2dhz586J/Px8Wd/U1JQAIMbHx8XXr1+FRqMRHR0d0uufPn0SYWFh4sSJE1IfAHHv3j3ZedatWyfsdrsQQohbt26JpKQk4fP5pNc9Ho8ICwsT3d3dUt3x8fHix48f0hiLxSLKy8uFEEKMj48LAOLRo0eLfu53794JtVotBgcHhRBCeL1eYTAYRGtr6y/cNSJSAvdoEdGqtGvXLrS0tEhtrVaLiooKpKWlycaNjIzA4XBAp9P5ncPlcuH79+/wer3YsWOH1B8ZGYmkpKTfqmdkZAQTExOIiIiQ9c/NzcHlckntrVu3Qq1WS22j0YjR0VEAfz0GVKvVyMrKWvQasbGx2Lt3L27evAmz2YwHDx7A4/HAYrH8Vq1EtHwYtIhoVdJqtYs+LtNqtbL2zMwMioqKcP78eb+xRqPxl/c2qVQqCCFkfT/vrZqZmUFaWhpu377t996fN+mHhIT4ndfn8wEAwsLC/rGOQ4cOobKyEpcuXYLdbkd5eTnCw8N/6TMQ0fJj0CKioLZt2zbcvXsXCQkJWLPG/ytx48aNCAkJweDgIEwmEwDgy5cvePXqlWxlKSoqCu/fv5far1+/xrdv32TXuXPnDjZs2IC1a9cuqdbk5GT4fD709fUhLy9v0TF79uyBVqtFS0sLurq60N/fv6RrEdHy4GZ4IgpqdXV1+Pz5MyoqKvDs2TO4XC50d3ejuroaCwsL0Ol0qK2tRUNDA548eYKxsTHYbDb88Yf86zMnJwdXr17FixcvMDQ0hCNHjshWpw4ePAiDwYCSkhIMDAxgcnISvb29OH78ON6+fftLtSYkJMBqtaKmpgadnZ3SOTo6OqQxarUaNpsNZ86cwaZNm5Cenr48N4qIloRBi4iCWmxsLJxOJxYWFpCfn4/k5GTU19dDr9dLYerixYvIzMxEUVER8vLykJGR4bfXq6mpCXFxccjMzMSBAwdw6tQp2SO78PBw9Pf3w2QyYf/+/diyZQtqa2sxNzf3WytcLS0tKC0txdGjR7F582YcPnwYs7OzsjG1tbXwer2orq7+F3eGiJaDSvx9UwEREf2j7OxspKSkoLm5OdCl+BkYGEBubi6mpqYQHR0d6HKIghr3aBERrRIejwcfPnzA2bNnYbFYGLKIVgA+OiQiWiXa29sRHx+P6elpXLhwIdDlEBH46JCIiIhIMVzRIiIiIlIIgxYRERGRQhi0iIiIiBTCoEVERESkEAYtIiIiIoUwaBEREREphEGLiIiISCEMWkREREQK+RMFilkC2VOEgQAAAABJRU5ErkJggg==", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "# Bias\n", + "fs = 500.\n", + "freqs = np.linspace(.5, fs/2, 4000)\n", + "phi = 0.9\n", + "\n", + "# Show bias as phi -> 0\n", + "phi = 0.2\n", + "\n", + "# Unbiased form\n", + "p0 = 1/(1 -2 * phi * np.cos(2*np.pi*freqs*1/fs) + phi**2)\n", + "\n", + "# Less biased form\n", + "p1 = 1 / ((((1 - phi)/np.sqrt(phi)) * (fs / (2*np.pi)))**2 + freqs**2)\n", + "\n", + "# Biased approximation\n", + "p2 = 1 / ((-np.log(phi) * (fs/(2*np.pi)))**2 + freqs**2)\n", + "\n", + "# Normalize scaling\n", + "p0 = p0 / p0[0]\n", + "p1 = p1 / p1[0]\n", + "p2 = p2 / p2[0]\n", + "\n", + "# Plot\n", + "plt.loglog(freqs, p0, label=\"Unbiased form\")\n", + "plt.loglog(freqs, p1, ls='--', label='Approximation: cosine')\n", + "plt.loglog(freqs, p2, ls='--', label=\"Approximation: cosine & -ln(phi)\")\n", + "plt.legend()\n", + "plt.xlabel(\"Frequency\")\n", + "plt.ylabel(\"Power\")\n", + "plt.title(\"More PSD Bias\");" + ] + }, + { + "cell_type": "markdown", + "id": "b18a2abe-5d34-4691-a4c1-eaf0319c7a68", + "metadata": {}, + "source": [ + "Substituting the natural log appoximation:\n", + "\n", + "$$\n", + "\\begin{align}\n", + "P(f) &\\approxeq \\frac{b}{(\\frac{(1-\\varphi)}{\\sqrt{\\varphi}} \\cdot \\frac{f_s}{2 \\pi})^2 + f^2} \\\\\n", + "&= \\frac{b}{(-\\ln(\\varphi) \\cdot \\frac{f_s}{2 \\pi})^2 + f^2} \\\\\n", + "\\end{align}\n", + "$$\n", + "\n", + "Remember that $\\varphi$ is related to $\\tau$ through the following, where $f_s$ is sampling rate:\n", + "\n", + "\n", + "$$\n", + "\\begin{align}\n", + "\\varphi &= e^{\\frac{-1}{\\tau f_s}} \\quad\\quad\\quad \\text{(e.g. the first lagged correlation of the ACF)} \\\\\n", + "\\tau &= -\\frac{1}{\\ln(\\varphi) f_s}\n", + "\\end{align}\n", + "$$\n", + "\n", + "Now we can substitute $\\tau$ in:\n", + "\n", + "$$\n", + "\\begin{align}\n", + "P(f) &\\approxeq \\frac{b}{(-\\ln(\\varphi) \\cdot \\frac{f_s}{2 \\pi})^2 + f^2} \\\\\n", + "&= \\frac{b}{(\\frac{1}{2 \\pi \\tau})^2 + f^2} \\\\\n", + "\\end{align}\n", + "$$\n", + "\n", + "\n", + "Finally, this leads us to the approximation of the knee frequency and a commonly used PSD form:\n", + "\n", + "$$\n", + "\\begin{align}\n", + "f_k &= \\frac{1}{2 \\pi \\tau} \\\\\n", + "P(f) &\\approxeq \\frac{b}{f_k^2 + f^2} \\\\\n", + "\\end{align}\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "id": "bde0188e-0a21-4f23-a56b-02db778f8e3d", + "metadata": {}, + "source": [ + "## Summary\n", + "\n", + "1) The AR(1) PSD is the the spectral representation of an exponentially decaying ACF, confirmed using the Fourier transform and Wiener–Khinchin theorem.\n", + "\n", + "$$\n", + "\\begin{align}\n", + "A(t) &= \\varphi^t \\\\\n", + "FFT(f) &= \\sum_{t=-\\infty}^{\\infty} \\left(\\varphi e^{-i 2 \\pi f \\frac{1}{f_s}}\\right)^t \\\\\n", + "&= \\frac{1}{1 - \\varphi e^{-i 2 \\pi f \\frac{1}{f_s}}} \\\\\n", + "P(f) &= |FFT(f)|^2 \\\\\n", + "&= \\frac{1}{|1-\\varphi e^{-i2\\pi f \\frac{1}{f_s}}|^2}\n", + "\\end{align}\n", + "$$\n", + "\n", + "2) There are two sources of bias in the commonly use Lorentzian spectral form: a cosine approximation and an approximation of $\\frac{1-\\varphi}{\\sqrt{\\varphi}}$. This leads bias as $\\varphi$ becomes smaller, timescales become shorter, or as the knee frequency approaches the Nyquist frequency. Instead, it's recommended to learn $\\varphi$ and $\\tau$ the AR(1) form:\n", + "\n", + "$$\n", + "\\begin{align}\n", + "P(f) &= \\frac{1}{|1-\\varphi e^{-i2\\pi f \\frac{1}{f_s}}|^2} \\\\\n", + "&= \\frac{1}{1 - 2 \\varphi \\cos{(2 \\pi f \\frac{1}{f_s})} + \\varphi^2}\n", + "\\end{align}\n", + "$$\n", + "\n", + "where\n", + "\n", + "$$\n", + "\\begin{align}\n", + "\\tau &= -\\frac{1}{\\ln(\\varphi) f_s}\n", + "\\end{align}\n", + "$$\n", + "\n", + "The knee frequency is not recommended due these sources of bias and poor interpretations, e.g. stochastic processes to not oscillate like is suggested by a frequency." + ] + } + ], + "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.5" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +}