{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "## ThinkDSP\n", "\n", "This notebook contains an implementation of an algorithm to generate pink noise.\n", "\n", "Copyright 2015 Allen Downey\n", "\n", "License: [Creative Commons Attribution 4.0 International](http://creativecommons.org/licenses/by/4.0/)" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "# Get thinkdsp.py\n", "\n", "import os\n", "\n", "if not os.path.exists('thinkdsp.py'):\n", " !wget https://github.com/AllenDowney/ThinkDSP/raw/master/code/thinkdsp.py" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "\n", "from thinkdsp import decorate" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Generating pink noise\n", "\n", "The Voss algorithm is described in this paper:\n", "\n", "Voss, R. F., & Clarke, J. (1978). \"1/f noise\" in music: Music from 1/f noise\". *Journal of the Acoustical Society of America* 63: 258–263. Bibcode:1978ASAJ...63..258V. doi:10.1121/1.381721.\n", "\n", "And presented by Martin Gardner in this *Scientific American* article:\n", "\n", "Gardner, M. (1978). \"Mathematical Games—White and brown music, fractal curves and one-over-f fluctuations\". *Scientific American* 238: 16–32.\n", "\n", "McCartney suggested an improvement here:\n", "\n", "http://www.firstpr.com.au/dsp/pink-noise/\n", "\n", "And Trammell proposed a stochastic version here:\n", "\n", "http://home.earthlink.net/~ltrammell/tech/pinkalg.htm" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The fundamental idea of this algorithm is to add up several sequences of random numbers that get updated at different rates. The first source should get updated at every time step; the second source every other time step, the third source ever fourth step, and so on.\n", "\n", "In the original algorithm, the updates are evenly spaced. In the stochastic version, they are randomly spaced.\n", "\n", "My implementation starts with an array with one row per timestep and one column for each of the white noise sources. Initially, the first row and the first column are random and the rest of the array is Nan." ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "jupyter": { "outputs_hidden": false } }, "outputs": [ { "data": { "text/plain": [ "array([[0.54878587, 0.80870645, 0.64066175, 0.63762133, 0.36665185],\n", " [0.58319314, nan, nan, nan, nan],\n", " [0.84357289, nan, nan, nan, nan],\n", " [0.18418847, nan, nan, nan, nan],\n", " [0.2622209 , nan, nan, nan, nan],\n", " [0.2106998 , nan, nan, nan, nan]])" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "nrows = 100\n", "ncols = 5\n", "\n", "array = np.full((nrows, ncols), np.nan)\n", "array[0, :] = np.random.random(ncols)\n", "array[:, 0] = np.random.random(nrows)\n", "array[0:6]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The next step is to choose the locations where the random sources change. If the number of rows is $n$, the number of changes in the first column is $n$, the number in the second column is $n/2$ on average, the number in the third column is $n/4$ on average, etc.\n", "\n", "So the total number of changes in the matrix is $2n$ on average; since $n$ of those are in the first column, the other $n$ are in the rest of the matrix.\n", "\n", "To place the remaining $n$ changes, we generate random columns from a geometric distribution with $p=0.5$. If we generate a value out of bounds, we set it to 0 (so the first column gets the extras)." ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "jupyter": { "outputs_hidden": false } }, "outputs": [ { "data": { "text/plain": [ "array([2, 2, 4, 1, 2, 1, 2, 1, 3, 3, 1, 1, 2, 2, 4, 1, 4, 2, 1, 4, 2, 1,\n", " 2, 1, 2, 2, 1, 4, 1, 1, 0, 1, 1, 1, 1, 1, 2, 0, 1, 4, 3, 1, 2, 1,\n", " 3, 1, 1, 4, 1, 2, 1, 1, 3, 2, 1, 0, 1, 1, 1, 2, 4, 1, 3, 1, 4, 2,\n", " 2, 2, 1, 3, 3, 1, 1, 1, 1, 2, 2, 3, 2, 0, 3, 1, 2, 1, 1, 2, 1, 3,\n", " 2, 4, 1, 1, 1, 3, 4, 1, 1, 1, 2, 1])" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "p = 0.5\n", "n = nrows\n", "cols = np.random.geometric(p, n)\n", "cols[cols >= ncols] = 0\n", "cols" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Within each column, we choose a random row from a uniform distribution. Ideally we would choose without replacement, but it is faster and easier to choose with replacement, and I doubt it matters." ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "jupyter": { "outputs_hidden": false } }, "outputs": [ { "data": { "text/plain": [ "array([79, 54, 81, 56, 89, 39, 70, 24, 47, 68, 18, 17, 84, 40, 35, 48, 5,\n", " 46, 98, 27, 36, 42, 28, 20, 22, 7, 49, 56, 59, 71, 78, 96, 64, 24,\n", " 28, 33, 77, 7, 18, 54, 3, 29, 71, 84, 54, 55, 66, 33, 13, 42, 66,\n", " 41, 64, 81, 19, 97, 66, 28, 55, 14, 73, 45, 39, 48, 85, 72, 50, 85,\n", " 19, 0, 35, 98, 0, 77, 33, 98, 34, 13, 50, 99, 73, 15, 7, 52, 55,\n", " 36, 95, 24, 45, 10, 0, 13, 42, 86, 26, 60, 9, 97, 36, 70])" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "rows = np.random.randint(nrows, size=n)\n", "rows" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we can put random values at rach of the change points." ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "jupyter": { "outputs_hidden": false } }, "outputs": [ { "data": { "text/plain": [ "array([[0.54878587, 0.51221941, 0.64066175, 0.15579406, 0.36665185],\n", " [0.58319314, nan, nan, nan, nan],\n", " [0.84357289, nan, nan, nan, nan],\n", " [0.18418847, nan, nan, 0.38864791, nan],\n", " [0.2622209 , nan, nan, nan, nan],\n", " [0.2106998 , nan, nan, nan, 0.99871071]])" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "array[rows, cols] = np.random.random(n)\n", "array[0:6]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Next we want to do a zero-order hold to fill in the NaNs. NumPy doesn't do that, but Pandas does. So I'll create a DataFrame:" ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "jupyter": { "outputs_hidden": false } }, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
01234
00.5487860.5122190.6406620.1557940.366652
10.583193NaNNaNNaNNaN
20.843573NaNNaNNaNNaN
30.184188NaNNaN0.388648NaN
40.262221NaNNaNNaNNaN
\n", "
" ], "text/plain": [ " 0 1 2 3 4\n", "0 0.548786 0.512219 0.640662 0.155794 0.366652\n", "1 0.583193 NaN NaN NaN NaN\n", "2 0.843573 NaN NaN NaN NaN\n", "3 0.184188 NaN NaN 0.388648 NaN\n", "4 0.262221 NaN NaN NaN NaN" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "import pandas as pd\n", "\n", "df = pd.DataFrame(array)\n", "df.head()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And then use `fillna` along the columns." ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "jupyter": { "outputs_hidden": false } }, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
01234
00.5487860.5122190.6406620.1557940.366652
10.5831930.5122190.6406620.1557940.366652
20.8435730.5122190.6406620.1557940.366652
30.1841880.5122190.6406620.3886480.366652
40.2622210.5122190.6406620.3886480.366652
\n", "
" ], "text/plain": [ " 0 1 2 3 4\n", "0 0.548786 0.512219 0.640662 0.155794 0.366652\n", "1 0.583193 0.512219 0.640662 0.155794 0.366652\n", "2 0.843573 0.512219 0.640662 0.155794 0.366652\n", "3 0.184188 0.512219 0.640662 0.388648 0.366652\n", "4 0.262221 0.512219 0.640662 0.388648 0.366652" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "filled = df.fillna(method='ffill', axis=0)\n", "filled.head()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Finally we add up the rows." ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "jupyter": { "outputs_hidden": false } }, "outputs": [ { "data": { "text/plain": [ "0 2.224113\n", "1 2.258520\n", "2 2.518900\n", "3 2.092369\n", "4 2.170402\n", "dtype: float64" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "total = filled.sum(axis=1)\n", "total.head()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "If we put the results into a Wave, here's what it looks like:" ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "jupyter": { "outputs_hidden": false } }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAABUSklEQVR4nO29eZQb53nm+3xAobCj926SzV0kRe3UElmKN8lLrqw4XuIk4ziOPcnk6jp2FieeTDLJuZ5xzuQmcTKJx7GvFec6HnvG2e0k3uQ1ciTHWkzJIkWKokhxEbvZC9ALGnsVqr77R9VXqCpUAVVYGtv3O4eH3Q10d6EaePHU826EUgoOh8PhDD6BXh8Ah8PhcDoDD+gcDoczJPCAzuFwOEMCD+gcDoczJPCAzuFwOEOC0KtfPD09Tffv39+rX8/hcDgDyVNPPZWhlM443dazgL5//34cP368V7+ew+FwBhJCyGW327jlwuFwOEMCD+gcDoczJPCAzuFwOEMCD+gcDoczJPCAzuFwOEMCD+gcDoczJPCAzuFwOEMCD+icjvC98xm8mM73+jA4nJGGB3ROR/iNfziJT3znxV4fBocz0vCAzukIZVlBriz3+jA4nJGGB3ROR5AUFUVJ6fVhcDgjDQ/onI4gKyoKlWqvD4PDGWl4QOd0BFmhXKFzOD2GB3RO2ygqhaLygM7h9Boe0DltIysqAKAoccuFw+klPKBz2kbSA3qhwhU6h9NLeEDntI1c1QJ6SVagqLTHR8PhjC48oHPaRlZqQbwkc5XO4fQKHtA5bcM8dID76I1QVYqP/cs5LGVLvT4UzpDCAzqnbSpVU0DnProrFzJ5/PE3XsAXn7na60PhDCk8oHPaxqzQC1yhu3IxUwQALGXLPT4SzrDCAzqnbayWC1foblzKFAAAi5vccuF0Bx7QOW1jUei8/d+Vi2taQL/KAzqnS/CAzmkbqVqrcuEK3R2m0HlA53QLHtA5bSNxy8UTLKBvFGVeDcTpCjygc9pGrvKyxWaUZQVXs2UcnI4DAK5u8sQop/PwgM5pG6uHzhW6Ey+taxUud18zBYDbLpzu0DSgE0IihJAnCSEnCCGnCSEfcrjPPYSQLCHkGf3fB7tzuJx+ROKNRU25qNstLz80DYAHdE53EDzcpwLgNZTSPCEkBOC7hJCHKKWP2+73KKX0jZ0/RE6/Y2795wrdGeaf33VwCgECXOW16Jwu0DSgU0opALbOPaT/4xOYOAYS99CbcmmtgMm4iMm4iNlkhCt0Tlfw5KETQoKEkGcArAL4JqX0CYe73a3bMg8RQm5w+TkPEEKOE0KOp9Pp1o+a01cwDz0VEVDgVS6OXMwUsH8qBgDYNc4DOqc7eArolFKFUnoMwG4AdxJCbrTd5WkA+yiltwD4MwD/5PJzPkkpvYNSesfMzEzrR83pK1hAn4iLKHGF7silTBH79QqXXeNRHtA5XcFXlQuldBPAdwDcZ/v6FqU0r3/8VQAhQsh0h46R0+ewpOh4NMQ9dAdKkoLlrTIOTGkBfX48iqvZMlQ+O57TYbxUucwQQsb1j6MAXgfgedt9dhBCiP7xnfrPXev40XL6ElnvFE1FQ9xDd+CS3vJvVuhSVcVaQTLu873zGT5Wl9M2XhT6TgAPE0JOAvg+NA/9y4SQ9xBC3qPf5ycAnCKEnADwUQBv15OpnBFAUhQEAwRJ7qE7wipcDpgCOlArXcyVZbz700/i/334xd4cIGdo8FLlchLArQ5ff9D08ccAfKyzh9Z9fvMfTuLY3nH89J17e30oA42sUISCBDFRQJEP56rjYp1CjwDQAvote8bx2ItrkBWKy3rzEYfTKiPdKfq108t46NRyrw9j4JGqKsRgAHExyBW6A5czRUwnwkiENf00ryt0Nkb30XMZAMAVHtA5bTKyAZ1SilxZxoV0vvmdOQ2RFRWiEEAsLKDEA3odF9cKODAdMz4fi4YQE4PGootHz2klvIsbJZ4o5bTFyAb0gqRApZpKKvPFxm0hKypCukKXFNXSaMTRPPT9eoULABBCjNLFl9aKuLRWxOHZBCRFxUqOd5ByWmdkA/pWSQYAUApcXuOXuu0gVbWAHhU1S4Gr9BqFShWruYrhnzN2jmnNRY/o6vxnXqblca6s80oXTuuMbkAvy8bHFzPcdmkHlhSNi0EAfK+omeOXNwAAh2YTlq/Pj0exuFnGo+fSmB+P4lVHtEa7l7iPzmmD0Q3opVrQeTFd6OGRDD6SokIUgojpST9ei17j0/92EdOJMO651toZvWs8iky+gn87v4ZXHZnG/EQUhPDEKKc9RjigmxU6D+jtICsqRLNC592iAIBzKzl852wa77p7H8JC0HIbq0XPV6p45eEZhIUgdqQiuLLhP6Bf3SzhH3+wgDNLW6gqPH8xyngZnzuUMMtl51jEc6XLc1e3cD6dx5tu2dXNQxs4WFI0pnvo3HLR+Mt/u4iwEDD8cTOsFj1AgJdfo03J2DMZw0ILHvrvfeUMvvLsEgAgEgrgniOz+MQ7b4PevM0ZIUZeoR/bM+5ZoX/quxfxwX8+1c3DGkhYUjQe1lQoT4oCa/kKPv/0In78tt2YSoTrbme16LfsGcdYLAQA2DMR8+2hl2UF3zm7ih+9eSf+x9uP4fZ9E/ja6WVUeKXRSDK6Ab2sqcibd49joyhjwzRXw410voJ8uQo+1cCKpFCEBLNCH42AvlWW8e0zK463fe6JlyBVVfyHV+x3vH3HWATJsIDXXz9nfG3PZBQruTIqVe/n77EX11CQFPzk7bvx5mPzeP112s/jy7pHk9EN6CUZMTGII3Na9cEFDyo9naugqlKufmzIeqdoTPfQR6H9P1uU8Y6/eBz/4TPH64ZqlWUFn33sEu69dgaHZpOO3x8Wgvj2f3w1HnjlQeNreyZioFRrMPLKN55bRiIsGLtKWWK6MAJ/A049IxvQc+UqUpEQDs5oAd2L7ZLOVQBoiSxODa1TlCA+Igo9W5Txzk89gVOLWwCATM56dffMlU1k8hLe8bJ9DX/ObDICIVh7Ce7VF2Bc8RjQFZXim8+t4J5rZ4ykK/sbcIU+moxsQN8qy0hFBeyeiEIIkKaJUUWlWC/oAb3MA7oZlhSNjoBC3yrL+Nm/fAJnl3N47z3XAAA2itaAvpbXPt87Gav7/kbsmdDu79VHf+bKBjJ5yWLbxMK8F2CUGemAnoyEEAoGsHcq1lShrxUqYGM2uEK3wpKiohCAGAwMtUL/qydewsmFLD7xztvw47fNA6gP6Ov65xN6stMrs8kwRCGABY8B/RunVxAKEtx7dNb4mqHQeenoSDK6Ab1URSqiPfkPTsdxoUlzEbNbAB7Q7UgKRUi3DmLh4FCvobu8VsB0QsRrr5vDeEwEAGwWZct9NvUEO7vdK4EAwe6JqKdadEopvn56GXcdnEIqUnvjiPFu3ZFmdAN6WUYqqr0QDs4kcHGt0HDSXSZfU2E84WRFVlSEBe2pFBeHe8nFwkbJKDkc158/Tgo9ERYgCv5fXnsmYp7muZxfzePSWhE/csMOy9fjbXTrSlUVv/53z+DyGm+0G1RGN6CXZEPZHJiOQ6qqxnxqJ7hCd0fz0LUmlpgYHOrW/8XNEuYntIAuBANIRoR6hV6UMRH3Z7cw9kxGPXno33hOK5dkZYqMdrp1FzaK+MLTi3jsRb49clAZyYBOKcVWuYpUtGa5AI0rXXhAd4clRQEtoA9r6z+lFIsmhQ4AEzGxTqFvFCVM+LRbGHsmYsiWZMvwOCeeuLiOozuS2DEWsXy9nXk6JX2MtMTHBwwsIxnQi5ICRaU1hT6jBfRGlS7pXAVCQFOhvMqlhqpSfdoiC+jC0Cr0TF5CparaAnoIGzaFvlFoPaCzyphmQ7rOreRw3c5U3dejodYVelnWAnlF5gF9UBnJgM7UD/PQZxJhJMNCY4WerxgT8biHXkNWtRc/84vj4eFV6MySm5+olSOOx0Rs1il02XeFC2OPEdDd7b9sScZStowjc/VNS8EAQTTUmu1V5gp94BnJgJ7TFTZT6IQQHJyJ43xDhV7GbDKMhCggxwO6gaxoiWTRpNBLQ7oBinVw7p6wK3RbQC9ImIi3brkAjRX6uZUcAODaHQnH2+Ph1na7soBeGdK/3ygwkgGdDeZiHjoAHJ5L4oWVxpbLTDKMeFjgCt2ErI9BYElRTaEP5/lZ3NSC7LwpoI/HRGwWapaLrKjIVaotWy5jsRBSEaFh6eJZPaA7KXRAt71a+BuwN+IKV+gDy2gGdN1ySZrqd6+dSyKdq7gO6crkJcwkwkhEBJ4UNSHrL/6QYPbQawrv6mZpaJY2LG6UkIwIlrrviZiIXKVqnIeNFpuKzOwci2I5675b9IXlHOJi0OLlm4mJrSl0NiWTe+iDy2gG9BKzXMwKXbt8fUFXP2YqVQXZkmwo9PyQesStUDEUOqtDD6Ig1SZS/trfPoNf+qune3Z8nWTBVuECwChPZKWL7P9WLRcAmEmGkclXXG8/u5LD4bmk67zzeLi1xHRZ/1tyD701KNVm6/RyychoBnRbUhQArt2hXb6+sFpvu7CmopmkljzNNykpGyWYMmWNRVFRAKVaxYSiUpxcyOK5pS1fI2H7lcXNksU/B2DqFtWeI+sFptBbD+jTCRHpBgH93Eoe17rYLQBaFh1lrtDb4vnlHP7Pzx7Hw2fTPTuG0QzoJWa51BT6jpQ2n/qF5XqFzmrQpxPhoa7iaAWWFDUUumk41IvpPEqyAlmhONcgPzEo2GvQgZq1smEo9E4E9HDdBEdGJl/BWkHCkR0NAroYbMlD51Uu7cHiysqWu13WbUYzoJeriIQClj2PhBAc2ZF0tFxYQJ9JhpEIh7iHbsLw0E1VLoA2HOrkQta436nFbP03DxDZkoxcpWpJiAK1wL1hKHRmubTuoc8kwyjJimNymQmORgrdnsfwSolXubQFO3/rHpbldIvRDOimtn8zR+YSeGElV7eRyBrQgzygm5AUW5ULG6ErV3FqMYtoKIhkWMCpq4Md0Gsli9aRuOMx5qFrL+KNDil0wNqdzDAqXFxKFgFWtsg7RbcbdoWz1sAu6zZNAzohJEIIeZIQcoIQcpoQ8iGH+xBCyEcJIecJIScJIbd153A7g3kwl5kjc0lsFOU6/5K9sKbitSoXvoZOQ9ITaUYdurExR8Gzi1ncOJ/C9btSxjKIQcVoKqqzXJhC15T5RkFCNBREJBREq8wktYDulBh9YSWHiVgIMw57Shla2SLvFN1u2BviWp8r9AqA11BKbwFwDMB9hJC7bPd5A4DD+r8HAHyikwfZabRtRULd19llrN3vTefLmIiFIAoBxMMCFL6GzoBZLkanqK7Qt8oyTl/N4sb5Mdw4P4YzS1s9zf63y+JGfQ06oJUIisGAocw3ijIm26hwAZoo9OUcjjSocAG0v4GkqMabrVe4h94ezObqa8uFarAIF9L/2eXpmwF8Vr/v4wDGCSE7O3uonWOr5KzQD+sB/awtMZrJSYZqSuoKNMfnuQBw99CfXciiLKu4aX4MN86nUKmqeLHJzPl+ZmGjhEgogClbsCaEYDwWMpqLNouSYcO0ynRS+x12hU6plly+tkFCFKhdJZV8+uhGp+gQVCT1gtIgBHQAIIQECSHPAFgF8E1K6RO2u8wDuGL6fEH/mv3nPEAIOU4IOZ5O9660Z6tctTQVMaYTIibjIs6tWgN6Ol8xAnoiUr+Edzlbxmcfu9S9A+5jpKq1yoUtWHjiojaC9ebdY7hx1xiAwU6MLm6WsGs86qiMzRMX14tS2wp9Kh5GgNQr9KVsGblK1bVDlBFvccmF4aHzq8+WKA+I5QJKqUIpPQZgN4A7CSE32u7idP1XZzJTSj9JKb2DUnrHzMyM74PtFFpStN5yIYTg8GyiTqGncxXjMpit+DInRj//9AI++M+ne5oM6RU1y0Wfh66XLT51eQMxMYgD0wkcnEkgEgoMdGJ0cbO+ZJExHgtZGov8biqyEwwQTMZFpPPWwNCs5Z/R6ghdo1OUB/SWMFe5NFqW0018VblQSjcBfAfAfbabFgDsMX2+G8DVdg6sW2iz0J0tF0BrMDq3kjeSnpRSbY5LwqrQzZYLU1KbpdFrOJLqOkW181OWVdywK4VggCAYILh+ZwqnBzgxurhR31TEsCj0gtRW2z9jOlHfLcpKFo/MuVe4AK0vuTA6RXlAb4mSpJ03RaVN59l3Cy9VLjOEkHH94yiA1wF43na3LwJ4l17tcheALKV0qdMH2wnKsgpZoY5li4Dmo+cqVSzpszQKkoKSrNQsl3C95cJeePbNNaOAPSkaNVV33DQ/bvp4DKevZnumXNqhJClYK0iuCn0irs1EryoqtspyWyWLjJlkuM5yObuSw1wq3PQKgOUx/FouZa7Q28I8ZbRXtosXhb4TwMOEkJMAvg/NQ/8yIeQ9hJD36Pf5KoALAM4D+AsA7+3K0XaAWtt/veUC1Cpd2OWtuQYdqAX0vENAz5Z65531CntSNBAgho9+0+7aAoYb5sdQkBRcGsB9laxk0V6DzmAz0bMlGZS2N5iL4aTQX1zN4/BsY7sFqHXr+i1dLFe5h94O5uXovUqMOkc1E5TSkwBudfj6g6aPKYD3dfbQuoMxOtdFobPL2XMrOdx77azHgK798UZRoUu21n+g1ql40/yY8TUjMXp1CwdnGlsG/UZtsYWb5RJCVaW4ojcftTOYi8EUOqUUhBBQSnEhXcBbb6urNaijVYVe89B5lUsrlGQFhACUAmv5/lXoQ8UWW27h4qGPx0TMJsM4u6xVajKVZK9ycVLooxjQDcvFFNDj4SDiekKUcXguATEYwOkBrHRZYDXorklRLYBfzGjPmU5YLtMJEZWqajzPMnkJuUoVB/T9t40wFHqLZYuyQgfSGus1JVnFXFLb8bpW6E2BxAgGdKbQ3S9Obts7gS+dvIqHnl2qKXQ9KRoNBREwraGTFbVW4TDSSdFaoVMyIuCG+TEEA7WvhYIBHN2ZxLMDGNCXNssIBgjmUhHH21kAv6DX2bdbtgjUBAR7/rF9t16ubgyF7nNERVlWwaoyeXORf8qSYlzFrfdIoTe1XIaN2rYid5/z93/8JvzCZ4/jvX/1NI7MJhEMEONFSwhBPCwYVS7mS6tscTQ99AABBJNC/7233ISoWN/6fmgmgScurm/n4XWErF7man6DMsM88wv6Ttp2G4uAWrdoJi/h4AyMfbcHPSh0lsPwo9AVlUJSVIxFQ8iWZFSqalvjC0YRVjyRDAt9nRQdKpjlkmyg0CfiIj73Cy/D666bw9mVHKbiIgKmF3PCtIbOnLgaSYWuqBb/HABu2TPuWCs9lRB7dinaDgWpaqheJ8a7oNBrAV1X6JkCRCGAXS62j5lQMABRCPjy0JndMqYLHe6j+6ckK4iGgphMiP2bFB02miVFGZFQEA++83b8wUNnjJnfjES4toaODfIKkBH10KvU4p83YjIeRllWUWwSIPuNQqVqJMOdYAr9kh50ox1QtvWWSwEHpuKuVwl2tJno3oNyyRbQeaWLf0qSgkgoiMk4D+jbxlZZhigEPF1OBgMEv/Oj19d9PW4K6Bn9Bbd3MobsCCp0WVGNfaLNYHNQ1vISYpOD89QrSorRAesEC4IlWcGOVKTh4CyvTMREBIhZoedxxEPJIiMmCm0q9O4E9KqioiQrjqM3Bp2SrCAqavN+Fjd7s+Ri9CyXUrWpOm9G0rQompUsHppNjGRAl6qqJSHaCGZF9HLWRSvkmyh0IRgwkuyd8M8BTUxMJbTSxaqi4qW1Ig7ONPfPGWZb0Av2gN4thf6p717Ej/zpI03vN4ivpZKkICYKmIyLPRsDMnoBvSy7NhV5JS5aPfRoKIidY1FjycEoISuq0SXajMmEFtDXB8xHL1YUI9HoBqs974R/zmDNRVc2Sqiq1FPJIiMWDvpKirK29VSXFfq51TyWsmXjDcSJ772YwS0f+gZ+7yvPDczIZUopSrJmuUwlwtgoSj3ZmTB6Ad1lW5EfEhEB+XItoE8nRYzHtOqAUavfdUqKumG2XMw89uIafvMfTnb82DpFvlJFvIFCB2qJ0U7UoDOmEyLSuYpR3+6nIcssOrzAukS7rdBZTqDRrJPHL2iVUH/x6EW8+9NPYmMArujYG2A0FMRUXISsUKMAYzsZuYCeK1cblix6wZwUzeS1SYxj0RBUCuRGbD2drKiek6JTeuWGPWH09dPL+NvjV3zP794uilLVGDrmBkuMdspyAbTEaCYvGdUzXkoWGTHRr0LX7suOv1tVLkZAL7m/Tp67msWh2QQ+/Lab8f2LG/ixj33XcdlHP8HOXzQUMK7SepEYHbmAvlV2Hp3rBxbQKaXI5CRMJ2oDk7JDUOlSlKp44LPHjQ7JRmgeurenUVwMQhQCdU/01ZyWQFrvU8uqUFGaKnSmzDtpucwkwkjnK3gxXcBELORrpEA83F5StGsKPd9coZ++uoUbdqXwUz+0B5/5+TuxsFHC108vd+V4OkVRP39RMWgK6Nv/JjR6Ab3UvkKPhwWoVOusYwp9XP+Zm0MwoOvscg7feG4FT13eaHpfWaGek6KEEEzFxbqk6MqW9sTvx0trqapCUlRjJK0b44ZC72BAT4YhVVWcXNj05Z8DukJvo2yxGx66olIjWeiW9FzLV7CULRuzf+46OIkdqQgeu7DW8ePpJEyhR0JBTMVrTWHbzUgFdEoptkpyw6YiL7B5LtmSjPWihJmEaNr+PvgKnXl/XuZpSz6SogAca3RXtnSF3ocBnS2J8K7QO2e5sOai55b8DzRrV6F3w3JZL0hgKaYtl4B++qo2M/+GXdqkTkII7jo4iScurPX1YnZ2/qKhIKYS3HLZFkqyAklR205cJfSa5CsbRVAKTCfDtYA+gOVWdph68rLxRvaRFAVQV9JFKcUqU+h9aLmwXEm8QR06YPbQO1vlAmjT+1pR6GVZheIxSV+WNUXeTcvF7IO7JQxZQL9+V2308t3XTCGTl3B+Ne/4Pf0Au8JhZYsAD+hdh6nn8baTorXOQAB6UlT30IcgoDP15EWh+0mKAqizXDaLsjEIqh8tF5ZYbKbQJ/XLbPsS6XZg3aIAcI2PGnSgtjnK6xq67bBc0qY3cjeFfupqFrsnopY3xrsOTgEAHu9j28VIiopa02JcDPZkhO5IBXSmANtVUUytsWUNU3HReCEMw4AulrDyEgz8JEUBLfCZlcuqSbWt96FdVVPojQP6a6+bxf/z1psM77cTTCdqz1PzKGIvxHyO0GWWQWq7FLpLQH9OT4ia2TsZw66x/vbR2Rsi60DX5rnwpGhXYRUo7ZaWJZlCX9OqQKaTYYhCADExOBQeOrvK8OLBygr15aFPJUQUJcUIIMw/B/pUoetXKc3KFiOhIN7xsr2WIW7tMhETEQwQEALsm3LeluRG3OcI3ZKs6CMxtL9lVxS6HtATYcGxyiVXlnExU8ANtjdFzUefwuMX1vvWR6+VLeoBPR7uSUf0SAX0jQ4FdEOhmywXQLNy+sVDf+zFNXz4a/bVr95gNcKekqI+FfqUrf2fBfSwEOjLskWvHno3CAS0qqD58ajvUbZ+R+iWJQURIWDYZ90K6HExiLlU2NGaPLOkrX28cT5Vd9tdB6ewXpDwwkp/+uglU9kiAEzHRW65dBtWUth2UlSvcrm8VoRomuMxFhP7RqF/8cRV/PkjF1pSNDUP3VtSVBS8q1IjYaQ/2Znlcngu0bJC//6ldWNzUqcxqlx6NB1yz2QMR3fUB7hmxB2WmTeiLKuIikEQQiAKga5UuWTyFcwktSY8p8ai01e15Sd2hQ5oiVGgf330eoXem4mLoxXQ9WA71oFOUUBTb9MJ0ZiuNx4N9c2i6Ey+AkWlvteQAWYP3VtS1JdCTzCFrgXyla0yxqIh7ByLtvQCWNgo4icffAz/9INF39/rhYJHD71bfPwdt+EP3naT7+/zq9DZLG9Au1rqloc+kwwjFQ05Wi6nr25hOqGtgLSzeyKK+fFo/wZ0m0JnM9G32yIasYAuIRoKtr2Jha2hAzT/nDEeC/WNQs946Mhzw4+H3kpSFKjNc1nZKmMuFcZkTGypbJH9nOeWtnx/rxcKRpVLb7b37BiLGJaeHwyF7rHKpawPlgK0gN6tKpeZZBipSMgxKXr66hau3zXmOH645qOv9eW8pLKsIEBqu3Wn4iIkRbXsHt4ORiygyx2ZtcHW0AGwvNjGY/3jobNA12hmhhvsxeal09BvUtReo7uyVcFsMoKJuIiNguxb0bBVgN2qUS5UqiAEHVlasZ0YCt1jt2jJEtCDLSn0R15I49Fzadfb07na3CO7h16pKji3ksONu9ztpbsOTmKjKOOF1ZzvY+s2JUm7wmFvRky4bLftMlIBfaMod6zxI2kE9NrPG4uKyJb8B6Vu0I5CNzpFm6g7SqmvaYuAtpw7FCRGUjSdq2A2FcZELARJUX1bROzxvbDSnRd5oaIgLgodWVqxnRhVLr4UuvZ3FFtU6P/9G2fxkW+dc7ytUlWQLcmYSYSRigrYKlctr5MXlvOoqtTRP2cc2zMOQBtN0W8UZcWyR5dZi9vd/j9SAT1bktpuKmI4KfSxaAhSVTW67npFUaoagdGt3tcNSqmpU7RxcK3ql76ix1kugHZ1oyWMKlBVitVcGXOpiDF4yq+iYY9vZavSlaauQqXaM7ulHfzXoas2D91/7iWTl1yfbyywMcvFnt95flmzzK7b6b6ViRUjtJIX6jZlSbFYuVM96hYdqYC+0SHLBag9ueyWC9D7AV2ZXO33+1XoRUmBolIEA6RphQSrLPGj0IFac9FGUYKsUMwlNQ8d8N/+nzO1kJ/vwqV4wcPo3H5EDAYgePgbMkomhdmKQqeUIp2vuD7fWA06q3IBrF3Vy1mtfLXREmz2htOPAd2cVAZq1uJ291aMVEDf7KDlwipdLEnRaH8M6DK3WOd8DtlnL7IdqQgqVbXhxhjms/oN6Kz9n01ZbEuhmwLIuS7UKBc8LLfoRwghiIlBH2WLCiJC61Uu+UoVUlV1zdmYAzrrRjX/7VZzFYxFQw0LFtgbTsnH0LHtoiRbt1o5vWltByMT0DUrQeqcQnfy0Ptk4uKah5kZbrAX2Y6xCIBalYcTbAaLn6QoUKvRXdHnoM+mIjVF04JCT4YFREKBhk0n3zm7iot6I5gftFnog2e5APpeUR+t/5E2FDqzVEqy4vhmYAno+sYwc/BfzZUdyxXNiMEAggFilAj2EyWb5RIXBQRIazmsdhiZgF6QFMgKNabitQtTbTNmy8UY0NVjyyVvtlx8KnT9zWinHtAbzXORFeah+w/oa3kJq3qXKCtbBID1gs83oJKMsVgIh2YTONfAcvmVv/4B/uzbzgm7Rgyq5QIAsbDgeThXWVZNCt1/lUvGclVY/zdkAX0qHjZ2+prV68qWlhxvBCEEsZC/TUzbRdmWFA0ECJIu5ZndpOkrkRCyhxDyMCHkDCHkNCHkVx3ucw8hJEsIeUb/98HuHG5jzq/mXStM2AJnFnTbJeFStqj9rt4qdPbimoj5f0KxNwAW0Bu1/8vMcvHRKQpolku+UsWV9RIATbUlIwKCAeLbc9wqV5GKhHBkNulquZRlBVvlKl5sSaEPpuUCaBuivIxvYAuOo6Je5RL03ylquSp0EBHpfBkTsRBEIWDYEebnZjqnla82IyoG+3JVYVFS6kpbWTXPduJFWlUBfIBSeh2AuwC8jxByvcP9HqWUHtP//W5Hj9IDL6zk8Lo/+Vc8ci7jeLvRJdohhb5nMobphGjpOu2XmeiZfAWpiICpRNj3JR97ke0c05JTjRV6i0lR3aZ6fnkLE7EQwkIQgQDBeDTke57LVllbWHJoLoHlrbKjZ8nU4cW0+xu+GwVpcC2XmOhNocsKhaLSWpVLyL+HnjZfFbr8Ddg4YMNy0Z+blFKjfLUZfnelbhf2pCgAx3r7btP0lUgpXaKUPq1/nANwBsB8tw/ML2wOxOU1ZxXGAnqntrK/++59+PYH7rFM14uGghCDgb5Q6NPJMFIRwXdjUbZktVwaKbxKG0lRQBvGNJeqqTKtuch/2WIqqil0wLnBaNXYNF81BrR5pVAZXMslHvam0MtV6+hXTaH7C+hWhd44oLONYey5yWbie1PoQl8GdHMOguHWEdtNfL0SCSH7AdwK4AmHm+8mhJwghDxECLnB5fsfIIQcJ4QcT6fdO8pagb2QzeNYzbBSwk4lRYVgoG4mDCEEqR68K9vJ5LXF1W4zMxrB7j/nyUNvLSk6pdtUi5slzJoCeivt/7lyFcmIgCNzWkA/59BgZJ7DfTHjvRJG1WulYwNquXhV6GXJGtBbUeiZJpVVmbxk5JuEYACJsGC8TtgbbrOkKKAp9JLch1UuTpZLxP/rr108vxIJIQkAnwfwfkqpfXDG0wD2UUpvAfBnAP7J6WdQSj9JKb2DUnrHzMxMi4fsDAvobJ2ZnU6Nzm3GeMx5QFdZVvDxh897TlK1QyZf0TryWlAI2ZKMZFgwLosbzaJoJynKmDO9iCfiIWz4TYqWZaQiIeyeiCISCuCcg0I3l3FeSHv30dkm98SAWi7xcNBTlQtrhKsp9KD/KpecZOSV7M85ZqmYNzClIrWZ6Kus2slDQI/2YVKU5SBidoUe9X+F3C6eXomEkBC0YP45SukX7LdTSrcopXn9468CCBFCpjt6pE0wAnrOOaCzTULtTlpsxnjUeUDXl05cxR99/Sy+fWa1q78fADK5CqYSIpIRwXcd+lapilQ0ZPjGjV48rXro5jVtZstlMi768tBVlSJfqSIVERAIEByaTTiOAEjnKiAEEALEV+liUX8ziw2o5RITBeMxNKJkWnAMtKbQ1woVY++pXZUWJAUlWbEG9GhNbDARZr5ac6OdpOj51ZyrJdsOkqJCpairoU9F+tBDJ9oQi08BOEMp/ROX++zQ7wdCyJ36z922OZeyouKyvj3ILaBvFGXExCDCQnfVltvExS+dXAIAvLRe7Orvr1S1ig6z5eInEbhV1jzpmIeNN5IR0P1VuaQiIQT13MOcKRE2EdM8dK/Hm5eqoLS2Nu3wbNLRQ0/nKpiKi9g7FfMV0NnVSWJALZe4GERRVppOJ6yNfq1VuUiK6muqYSYvYe9UTKu9tqlScw06w2xN+rdcWgvov/o3z+BDX3qupe9thH0WOmMsGnKty+8WXqTVywH8LIDXmMoS7yeEvIcQ8h79Pj8B4BQh5ASAjwJ4O93GCVWX1wqoqhSJsIB0zsVDL8odS4g2gg3oMrNekPBv57Xqm0stlM75gU1ZnNYtF1mhvmbLZEsyUhHB0zztVjtFAwFi/C1mbQq9qlLkPHY3sqsPlmQ7PJfAUrZcpxDZlL+D03F/Cl1im9wH03KJhQVQWkt6usHWARp16PqQLsnH0pBMTrf5HPI2RkBP1P7Wmr+s/f1WtspIhAVP5aGtVrlQSnEhXbAkbzuFfRY6gwkNp7r8buGlyuW7lFJCKb3ZVJb4VUrpg5TSB/X7fIxSegOl9BZK6V2U0u91/9BrMFX2sgOTWCtIju3qm0Wp63YLwBS61TZ46NQSFJViOiEaVxLdohbQRaOBw/wCW8qW8JFvveCqvljVSCgYgCgEGk7rY5ZL2GdSFKjZLrOWefL+5l+wS3bm9x+eZYlRq0pnc7gP6AHdq/IceIXOFrE0sd2MBcdircoF8L6GriwryOnLXpzyNs4KXTDup9Wge5v5Hg0JLVku6VwFJb0fodO4KfTa62/7fPSh6BRlAf2ug1Og1Hlk5WZJxkR8GwJ6NISCZL3M+vKJJRycieOea2dxqQsenhlWbTCVMLdY115gXzm5hI986xyubDi/sWyVZOONLy4GG87TbtVDB2qJUauHrv1er/NcjICuH++RuQSA+kqXTI4F9AQqVRVLLpVQdlgCe1CrXNhqxGYBpVKn0LX/vVoFbBTytGk0rhl21WwdNW3y0HNlS7BvhKbQq777CS7rVmc3PG3jDdHBQwf8j99oh6EI6OdW85gfjxqb0VcdbJfNotSxLtFGsCoaFpRWt8p4/OIa3njzLuyfimE1V+lqpQur6GCXv4BVobOyTregyTovAS2p1lChV7UXVagVha6/uM0v5AmfExftlsvuiRhCQYKLpjdNc4UFS9pd9Fjpkq8MdpWLvYHHDbtlEDYUujclzGyMaZfKqnS+gqDJZmPHlqtUoagUq7mKp4QoO0aV+l9iza6Mu7GvgFlW9VUu2z+gaygC+vnVPK6ZTRhPCqfSxc2i3LEu0UbcOD+GAAF+5W9+gHyliq8+uwRKgR+7eSf2TWkBpZu2C1Po00mxroEDqCWgnIJmVV+ZxRR6Iiw0VOitJkUB4OiOJI7uSFrUfW2bkbcXAAtULHAFAwR7JmN4yXR+t0pVSIqKmUQYB2f0gO6xFr3X+0TbxbjkbxJQSpL2dzRXuQDeFXrtqlB0rL3WchiipQnP7C+vbnm3XGLGxEV/tgurbml1z24jjPMn1idFge0d0DXwAV1VKV5M53FoJmE8KeyVLpRSzXLZhoB+694J/I+334qnLm/gXZ96Ap9/ehFHdyRxeC6J/dsR0HMSYmIQMVFwVGg1hV7/JGOKlwWCWDjYUKGzF7zfOnQAeN+9h/DlX36F5WtshK49B+GGXaEDwL7JGC6Zzm86rz3emWQYs8kwYmIQFzwmRgsDXrbI/v7NSleNpGioVuUCeFfBbP6+YbnYqlxWHea0sGB3dbOMkqz48ND1ZL3PShfza67TipldcTs1FgGtrYFslYEP6IubJZRlFYdmE8blu91yYZd222G5AMCP3bILH3/HrXh2MYtnF7P4sVt2AQD26pZQN2phGWuFijEwzCkpw65enBKPdsUbF4WGZYutdooCWletYHsjSIYFCAHi20NPRmpv1Pum4nhprWBcVq+aEnKEECMx6gXWNh8f0CoXJ8vNCbsH7FuhF2yWi4NCt3vkzN8/n9auluZ8WC6A/5nol03lwp1WzK4eukNRQrcZ+IDOnhCHZhMIBQOYiot1Cj27TV2iZu67cScefOftuHXvOH78Nm30zVg0hMm4aFGQnSaTrxj+tFNShin0NYegyZQLU0/NSsTaSYo6QQjR5rl4VeiVKiKhgOUNZf9UDAVJMRLjaVuNs5+AXpSqCAuBujeeQcGrQqzICgipVSuJQS0w+VHocTGIqBhEKhpCUVKM5wbAFLotoOvPsfN6Atu75aIFSWZzeOXyWgH7dUGV7fCspbJL2WI0FIQQINxD98OLq7WADmhKbNVWxbBpBPTtUeiM1143h39878uNyYUAsG8q1lWFnslJhkKPhIIQhYChEPKVqtEK7qjQS8xy0RV6uHFSVNJb/4VA5xYoT8ZEXwo9FbG+SdfyFNo5ttdAH5yO48p60ZP6zFeqA1uyCGgWihAgnhR6RKhtrG/FQ582Jilq54vZPIpKsZavV+hMNLBRDV4mLQIw9Ud4V+jZkozNooybdo8bn3cSt7JFQoilmmc7GPiAfn41j8m4aCTUZlOROoXOFN92eOjN2D8V73pS1DyjXas6qDVwMJxa7LNGGaDuoXsoWxSDASMQdILxmPd5Lqyr1cw+w9bSznE6X4EYDBiP6cBMHCr11rGrDeYaTLsFqA2La5oUtS1nEP1WuRQqRl+BYfPov3OtUIFK6xW4odD1gD7jYdIiUFPBfjx0liS/eX5MO7YO14WX9MY9pwY0rdGKe+ieOb+qJUQZs8lwXZULm0++nZaLG/umYriaLRmXaZ2kqqhYL0qYMdX7anXB1pkZoaDzIgl2vzGvCr2qtlTh0gg/81zYpEUzuye09nOzQmf+OQAcmNaeK15sl/wAj85lpDzM89G2FdVCgW+FbroqtCfi2XPOHrCZkr+YKSAsBIzPm9FKlQvr/bhptxbQ21XoCxtFvP9vfmAcA/PQnRrstBHWXKF7glKK82mtZJExmwwjk69YugE3jcFc22u5OLF/Kg5KtSdFp9koyqC0Np4WgGUNFksWXzOTaKzQTUnRsqxCcemslBW1pYRoI/zMRHeyXEQhgF3jUSMJls5VLIu8D0x5L10sSoO7rYjhZYRyyTbL22+Vy1rBZLlErb690RdhU+iJsLZzs6pSzKbCnq/yYiHt7+Gn9JBdjd2wK6UfW3sB9uHnV/FPz1zFM1c2AWgJ2mgo6PgYWhlh3Q4DHdDXChI2i7LhnwNaQK+q1BKwNnuQFHWDWQKXMp0P6BlTgwdDG1NqtVyu25ly9Km3SjKEADFUUG3iorPCkxW1YwlRBpuJ7qU930mhA9qbJks8p/UZI4yxWAhTcdFxobSiUsvvzVeUwQ/oHkYoV3QPneGnU1RRKdYLEqYNy8Va2ZHech68xewgAJjzaLcAzlUulFL8xt+fwBMXnOcBXl4rYDoRRjISQjIidECha6sTz+s7bO2WlZntnrg4sAFdVSn+4pELAIBr9eUGABybizaLMhJhoePBpxWMWvQuTF2sBXSz5RIyhgOtbFUQDQWxZzKGbEmum3nDPGmmNGoTF53VkFSlHT+nE3ERKvV2WezkoQNaeehL+mV2Jl+/2uzV187giyeuWgalqSrFL3zm+3j3p580vlasVAe2ZJGRjDTfa9mOh75ekKBSmJKiVg+dXRU6tfaz+3pNiAImD92k0AuSgr9/agFfeXbJ8XsurRWNCpdOLJ1Y2NQCOkvoliS1LiHK0HIY3ENvSKFSxS9+7in8+SMX8FN37Mbd10wZt8061KJvFqW+UOeAdpWQjAhdqXSpdYm6J0XnUmFMxUVQh6CZLVUtXiZT6G4+ejcsF/bCT3uYirflqtBj2CjKWC9IWCtIFoUOAL9531GEAgQf+tJpo179M49dwsNn03jiwrpRcjfIC6IZqUio6bS/smwNSMxD92K5rBWsV4X22vd0Tttva6/R1u6rnVsvq+cYRmORKaCz57FbXuSltaLRA9KJqhOm0NkQuLKsGE1Zdsw5rO1g4AL6lfUi3vaJ7+Gbz63gg2+8Hn/4tpuN2dpA7clhrnTZLMl9E9AJIRZLoJMYkxbjtql2xmYYbWYG68i013ubB3MBNYXuVukiVdWWukQbYbwhu2yeYpT1OdN2Dx2olS4+fXkDlNarw7lUBO9/3RE8fDaNb51ZxfnVHP7goecxFRchKarxQi1IysArdC9bc0qSNSD58dBZlyircomLQctM9EZzWthzzetgLkAb7xAWApaiAlZX7jT4riwrWN4qG1fGnVjcvLhhU+hNLBepqnalCMKJgQvoL6zksJQt49M/dyd+/hUH6hIR7PLNvEdyY5sGc3mlW7Xo9hI9wPqEWt0qYy4VwWTMeWZKtmS1MFgwa6TQQ0Jnq1ycrrCcqHW11itolqc4fnkDgHPA+Pcv34/Dswl86Eun8Wt/ewIxMYiP/8xtALSF45TSoVHo9iULUlXFXz/5knEloilMk0IXvFe52K8KmTeeMyl0+xWS+dgA701FDHvDGwvQixulumO+olub7DnR7lq4sqwgk69gIhZCJl/BRkFCSVKMZK0dexlntxm4gP7a6+bw6G/ei1cfcd5JGgkFkYwIluaibLF/FDqg+egLGyVLN10nWMtLmIyLljc58xNqZauCuWTYGCNsT4zaPWk2NtYtKSp1ISlq5EBcNk8xanNnHDz0ST2gX1oH4BzQQ8EAPvTmG7CwUcKzi1n8/o/fhDv3TyIuBnH66hYqVRVVlQ5+QHdYsvDouTT+8xeexZdPXgVQH9AJIRCDAW8K3TERX6u91hR6k4Duse2foS2/rg/oKkXdWGh2JcyeE+0q9EXdP2fx53w6j6KtSshMbYQxD+iuOF1mm5lNhi0BYaOPPHRAUwuKSo1Lt05RqFSRsClW9oRa3CyhJCuaQjemGtZbLuZzy8bGuiVFu1Hlkghr25KaWS725RZmYqKA2WQYJxeyAOCqEH/4mmm8955r8MuvOYT7btyJQIDg+l0pnFrMGgFj0C2XpK1zE6i9Wf7zM1pAL8n1G+vDgre9opm8pF0Vmp53bHmFMbrY5fyz6adzPpKigL5XVK49HnOwtG8EY1fCzHJpNynK/PN7rp0FoPnoZUlB1MVDZ7ZSdpsSowMZ0Jsxm6x1i6oqRba0PevnvLKfzeXusO1SkOrL7FjAM7dYO80dp5TqC6Jr32946K6WC+24hw5ob8grTSwXp0mLZvZPxY3xvo082v9031F84EeuNT6/YdcYnlvaMhTtwCt0h4mbGf218ei5DNbyFS0panvjEoWApyqXNX12kOWqUA+a+UpVm6ToErBZsPOTFAXqLReznWHPTb20XkQyIhiCbsxh1owfWP/IDx3QrubOreYc3xAZXgekdYrhDOipsFFz/dCpZagUOGIqbew1eya0y79OK3SnMjsWoNnMm9lkBJFQEHExaFHolaoKSVEtSdF407LFzle5sGNMN1PozEN3WSvIPNOkS4WFGzfsSqEoKTh9dQvAEAR0W6MPoNkkAaLVkH/l2SV9lov17+hdodeGwRm/U6+sqi1/dg7Yb711Hv/tLTcaV4xeiYbqPfQA0f7WdoV+aa2IfVMx4w2n3aUTixslCAGCHakIDs0mcH413zQpCnAPvS2Y5SIrKv74G2dx7VwS99+0s9eHZTCTDCMYIFjOeluF5pWCpNTN7mZPKDYzg13e2jsy7V2iQK3m122Erma5dDYpCmhvyM2Sos0UOgvofiooAG1BCQCjSWXwA3q9h5vOV7B/Oo4jcwn8/fEFAKjzgDWF3jygp22zg9jv3CrLjrtEzewaj+Kdd+3z/mB0omLQ0vrPkvkHpuN1lS4X0nmj6gkwLZ1oMcAubJSwczyCYIDg0GzSsFzcRIPXJSOdYkgDegRSVcX/9+hFXMwU8Bv/x7WW0sZeEwwQzCXDWOpwQNda1Z3XYNUsF00t2WemsCecWaGLQgBiMGBMaLTTjaQoYLXM3GjkoQO10kU3/9aNQ7MJiEIAT1zUEqqD7qE7KcRMTqvNf/OxeTy7qOUZzJ2iABAWgk0V+npBwpmlnNFSb/6dWyXZpND9/Q2awfaKMrJ6ue2+KWtAX9kqY2GjhFv3jNeOTQ+wLSv0zRJ2j2ti4fBcAstbZeSlquNgLsBseXEPvWWYZ/eRb72AO/ZN4LXXzfb4iOrZMRbB8lank6LuCv3KRhGJsGCMg52IOSt0u+KNhYMNW/+74qGnwihKCvINlmvkylUETWMK7LAkmF+FHgoGcN2OJM7qc7oHXaE7JUXZuNs36YtXgPpZ3l489G8+twxFpXVXvyl9UfpyVnt++/0bNCMa0mYMMbJ6Mv/AVMxSuvik/qZ854FJ4761tXCtBdiFjSLmJ7Rx2If1kSOU1o/OZURCQYSFAFfo7cA8u0pVxW+94WhHx7t2ip1j0e4odNsLMxIKIBQkoNTaYj0ZFy1LLlgyabfu7zO0rUUuVS5daP0HzM1F7udnqywjGRFc/7Z7W7RcAOD6XWNge4QHfdpiXNSGYFksF73yZM9kDLftHQdQH5DCQsBIKrvxlWeXsW8qhut32hU6y9sUIAoBy1VfJ7ArdNYQt3/aOhr5yYvriItBy/GNOXjolzIFPKX3LDSiUlWwmqtgtxHQa3m5Rnma7RzQNZwBXQ9cr7tuFnfsn2xy796wYyyC5Wy5YxvIVX35bcymKAkhxoo28xCkSZuH/vzSFsJCAAem45bvt794zEhdaP0HnLt97ThNWjQzFg3hPa++Bm+8eZfrfdy4cb4WAOwW1qARCBDLxM2yrCBXqRpvdG8+pm3Tsreui0IAFdk9oG8WJXzvfAZvuHFn3Zsqe76dT+cxk/A+SdErTo1FLKADtVLFJy+u47Z9E5aNU04W1Ie//jze/7c/aPp7lzbLoBSYH9cC+vxE1DhvbklR7Xe2PxDMK0MZ0A9MxfErrz2M//qmG3p9KK7sHIugKCkNL/3+/vgVnNI9zmawgf9Oni9TTHaFXpAUoyX5+eUcjswl63IN2kx0N4XeJQ895bzs24zbpEUzv/WGo7h934Tv33/jrjHj40G3XACWpNSeZ/YBbm85No+33jqP22znqZlC/+ZzK6iqFPfftMPh99US8Z22WwAteFaqtbHO2VIVqWjIsNkuZgrYLEo4u5LDyw5YBZ1TlcuFdMEYYdAIVoPOrmKDAYJr9F0MbpYL+53bNaBrKAN6IEDw668/Umcf9BM7xjQV2qjS5Xe//Bw+98RlTz+vqPvNTgHIGFNq6shjtehstPDzy1s4uqO+tDMeDho/28yV9SJylSp2jfurIfaCV8ulWYNZq1y7Q3tjY3NDBh3zgC5WecIqU8ZiIfzpvztWV1rYTKE/dGoZ8+NR3DQ/VncbExDZktzxhChgWnIhK3r/hKbQJ2IhpCICLq0V8P1LmoVy54Epy/caaxn1gE4pxaW1Akqy0nRpxuImsyVrKyWZj94ooI9xy2X42akH9KWsc2KUUs1C8frOzlS0k0XgNDNj0tT+n85VkMlLOGrzQgGtuchJoX/7zAoAbRRDpxmLhiAKAcs8HjteFHqrREJBHJ5NICY6Ly0YNJKR2vwStjy7mXIOC0FXhZ4tyXj0XBr337TDdakDozsKvdbwVpa1/olUVMunHJjWVjw+eXENYjCAm3fXv+GYA+zKVsVIsDZbTr6wUUKA1MQYABzW+1saWy7WCY///Mxi1/YK84DeI3boi6PdFLqkaJeUXt/ZWa24vcoFqJVqOSn0jaKEs8taRcd1TgpdDDrWoX/rzCoOzSbqPPdOQAjBTCLc3EPvcLLNzA/tn8Qu03LvQcbc7u40e8UJTaE7K9Zvn1mBrNRXtxi/z/R38dsF6gWmhkuSYlgnLNm5fzqOi5kCnry0gWN7xp3H9po8bXOZY7Pl5IsbJewci1psRqbQGy0TN1te6wUJH/i7E/jfj3u78vZL04BOCNlDCHmYEHKGEHKaEPKrDvchhJCPEkLOE0JOEkJu68rRDhGzyTAIgWulC7v881peVZs94hDQI/WWC+vuWytIeH5Z64q81iGgx8JCXVJ0qyzj8QtreF0X1DmjWXNRNxU6APz2/dfhrx+4q2s/fzsxL4pmbf/27k47jTz0rz67jF1jERwz1Xdbfp/p79INhc4sFy0HZQ3o+6biWNws4dRiFj90wDl/MmbytM2dpV4UOkuIMl5zdBb//Sdvwa173XM1bGsRpVpnblWleMut800eZWt4UehVAB+glF4H4C4A7yOEXG+7zxsAHNb/PQDgEx09yiEkFAxgJhF2VegsQDdbTsBgI26dttTXPPTai8tQ6HpzyEwybNlFytAUulWp/evZNKoqxeuv7159v9Oyb4aiUuQq1a556IB2Ce23Jb1fMU8/TOe1hRNhoXH1jpuHTinF4xfWcO/RWVc7ipVKAp1vKgJMa+jkeoV+YDoGSrXniN0/Z6RMExfNs1+aKvTNksU/BwAhGMDbbt/dsHExFQ1B0avQ/vkHizgyl6gr9ewUTQM6pXSJUvq0/nEOwBkA9reXNwP4LNV4HMA4IaR/eu37lJ1jESy5JP5qAd2jQq+4K/TZZBiiELBc/o5FQyBEexK7JUQBzcIpyYplUfS3zqxgKi7i2B7/FSReadQtmm/S9s+xkowIyFeqUFSKTL7iSTWHhSAqDgo9k5eQr1Qte3ztsFJJwN96Oa/EzJZL0Wa56JUuAQLXCifzCN1LmYLxxt1oObmsqFjKloymIj+wYzu1mMXxyxt4y63zXcvN+PLQCSH7AdwK4AnbTfMArpg+X0B90Ach5AFCyHFCyPF0Ou3zUIePHWMRLG06J0UNy8Vj/aqh0B2SM+942V586ZdeYUncCEGt4SOdr+Dcah7XuSgG5g2WdD9VVlQ8/PwqXnN0tqvjFGaTYWRLsuOml2aDuThW2HnKl6vI5KSm/jmgKXSpqtb1SRjjaJvkTljepjuWC0uKOil07bhu2DXm6mubk6KX1gq4ZfeYJm6K7q+15WwZKkWdQvcCu5L8X7pvbu7Q7TSeAzohJAHg8wDeTyndst/s8C11HTOU0k9SSu+glN4xM+O8oGKU2DkWbWC5aAG6UlU9jTEtNChbjImCoz8+GRPx9OUNSFXVXaHrFg4rXfz+pXVslat43fXd888B581TjEbbijj1mJcspPW2/2YYW4tsKp3t7dw/1SSg60HMy5uHX2qLoqt1Q+XGYyL2T8XwmqPudiCrOlFVistrRRycSWA8Gmqo0O016H5gb25fO7WMOw9MdrWc2tMrghASghbMP0cp/YLDXRYA7DF9vhvA1fYPb7jZMRZBrlJFriwbl6iMokmZ5spVhBONPU9m0bjNNnFiIi4aLc9HdzgrdGOErv7zv/XcKkQhgFcenvb8e1qh1i1axp5J6wtgQ1+d108z7vsZczNNpsHCCTMsoFeqqsVvv7xWRDBAmirVVCSEybjYlcYzow7dlBQ1X6197f2vavh7x6IhqBS4kNHqz/dPxTBhG1Znh81BtydFvcDebKoqxVu7lAxleKlyIQA+BeAMpfRPXO72RQDv0qtd7gKQpZQudfA4hxJWi77i4KObFzN78dELlSoEn40wzDsMBgiumXVWXDHTCF1FpfjmmWW84tC0Y3lkJzG6RR0So2zTvFMSl1MPCyjpfMXS9t8It72iF9cK2D0RbRqod09Ecc1M50taAWuVS7YkIxkWLPZfJBRskqTUnrsnFzYBaPbRpG1YnZ3z6TxCQdKSh87ebMRgAPff2N3UopdX5csB/CyAZwkhz+hf+20AewGAUvoggK8CuB/AeQBFAD/X8SMdQnakWHNRGYdmrZaHuVTQS6VLUVJ8N8KwZdHXzMRdqx6YhVOoVPHhrz2PK+sl/OZ9Rz3/jlZpNM/FaI7hAd0TLHl8Ia3ZJdNNShYBGDN67DPRL68VmtotAPBf33QDqkpn5hTZYbXlrMrFby6F+e0nrmwC0OyjibhoLJR24uxyDodmky1dcbDfd+/RGWPtXrdoGtAppd+Fs0duvg8F8L5OHdSosFNvXHGqRS+ZLBcv3aKtbKif0BW6m90C1NTQ/3r8Mr58cgnvuntfSwOv/DIVFxEMEMda9Ey+AiFALOvyOO6wgHIhrc3E9+Jrszd4s0KnlOJSpojbG9RcM7o5AycsBBAgmuhhbf9+YFcsJxayCAUJdo1HMRkTDcXuxPNLOdx9jXMZZDMmYiG86+59+Mnb9zS/c5vwV0QPYbaCU2LUPE3Oj0L3A2v/P7rTfT0fe2F++eQSXn5oCv/3G+0tCN0hECCYTojOlovDHkuOOyyAMYXuxXKpKfTa83CtoJUsNqtw6TaEEK2cVlKNSYt+YIr+uaUt7JmMIRgg+gYvrfnH/rzKFmUsb5UdCwu8Hu/vvvnGlr7XL7z1v4dEQkFMxUVHhW5ZgushoBekFhS6brlc10Chs5+5fyqGj7/jtq4kudxwq0XP5L2V3nE0EsxyyfhR6PUe+iWPFS7bQVQMoiRXdcvF3/OevQFIVRUH9McyGQ9BUlTHuUWNOqn7Da7Qe8zO8Yix2cVMyeKhN7dcihX/Cv2ug1O499oZ3L7f/RJ611gEH3j9EfzYLbswvs1VJbPJMK46vNlpCp0HdK8EAwTJsICVLW9t/4Czh866Knut0IHaTPStUrVlhQ7UVhWaO6ft9etse5VbaW8/wRV6j9mRct5cVJQUTMS0bk4v81wKUtX3dp09kzF8+ufubNhCTwjBL7/2cE9exLOpMNKOHrrkKbHHqcESo17a/gFnD/1SpuCpZHE7iIaCRpWL34CeDAtgrsqBaa0kllV8ObX/P7+cQyoiGEUM/QwP6D1m51gEyw5liyVJQVzfAeqlW7RQqdZtKxp0ZpIRrBUkVBVrYi7jsGme0ximSr12bjp56Jc8lixuBzExiGxJRklWfAf0gH7FApgUOgvoDrXozy9t4ejO1EDkbHr/lxlxdoxFsFmU64brF/RN4tpyAi8KXUFiwNel2ZlNhkFprUwR0B5npapyhe4Tv52bjh76WsEIgL0mJgpG/0YrO0tZ+SAbFTAZc57nQinFCyv5gbBbAB7Qew5rLrKr9KKkICoK2nICL1UulWrXm322G3ZurppyDMb41zhX6H5giUMvbf9AvYdOKcXlTBEHpvpjC1gkFDSsylZm+qQiIYSCxHiOTbhYLgsbJeQr1YFIiAI8oPecHS6bi0qSglgoaFkf5oaqUhRlxXGf6CDDZl6wORqAuUuUK3Q/MIXutRkrbAvoawUJuUq1jxR60Lh6aCWgj8dC2DMZMxZIpyJat6l9Jjpb/jIoCn24JN0AYjQXbdYr9J1joYZLMBjlqgJKMXQeOmuzXjQF9LS+zJd76P5gSVH/HroWNNmUxW5sqGoFc0VXK5bLL9172DLJkxCCiZiI9YJVPLEKlyNzPKBzPMBeYGw1GKMkK4iKQYSCAeNJ5UbBmIU+XAo9ERYwFg0Zy3mBmkLnAd0fTMV6zT2w9nq2qPtiRvsb7OsTyyXaZkB36vqcjNdPXDyztIXdE9G64Xn9CrdcekxcDCJA6mvNi3oZYjIiNE2KFiX3faKDzu6JqMVyyegKfVi2CW0XfpOiqUgIrz4ygz9/5AJOLWZxeY2VLPZHQG9XoTsxEaufuHh2OTcwdgvAA3rPIYRoOx9tPrmWFK1VudgXDZgxFPqQVbkA2rjSRZuHPhYNGZYAxxutLJz40393DFNxEb/4uadwYiGL+fFo35x3s3jpVECfjFsnLlaqCi5kCgOTEAV4QO8L2MB9BqVUS4qKQSQjgrGP0I3hVugxLG6WjDe0Nd5U1BJ3H5zGG27cgcOz3oPTZFzEx3/mNixtlvHIC+m+6BBlRHVLKKbbkp1gIi5akqIvrhagqBTXNhiN0W/wgN4HaKWJNVtFUlRUVarVoevqo5Htkm+wrWjQmZ+Ioigp2NDXg6V5239L7J2K4RPvvN3iPXvhtr0T+J0fvQ6ANs+nX2CWSycXhU/GRGwUtU1GAHB2RZvhMkiWy/BFgAHEXprImoxYHTqgDehiJY52mHofVssF0CpdJuMi1vKVgboEHgb+/Q/vRzBA8MMtjo/tBuyNqVN2C6ApdEWl2CrLGI+JOHEli0go0DeVPV7gCr0PSEUFy8xz8zo5ll1vVItu7BMdSstFD+h6pQuftLj9EELwrrv31y1h6SXMculkQGfjpFlz0SPn0njZgam+GHXglcE50iEmFbEmRc0B3Vjw22DJRSv7RAcFFtAXNkqQqtr8a94lymH5olaaitwwJi4WJSxulnAhXej67txOM3ySbgBJRa1JUcNyCdUUeqP2/4I0vB76WDSEuBjEwkbJUE7TSZ4UHXW6YbnUJi7KOL+qzY5/1ZGZjv387YAr9D4gFQmhICnGVEFz1QorN2uUFC1WFAQIfC2IHhQIIUalC2u+4gqdYyRFO7iG0DwT/ZFzGcylwjg8m+jYz98Ohi8CDCAs8cmqVYp6S3IsHDSy+M0UelwUBmK8ZyvM681Fa7pCn+EKfeSJdVGhp/MV/Nv5DF5xaGbgXlM8oPcBzAdkPnnJ5ImHhQBCQdJUoceGsMKFsXsiisWNIp+0yDFgVmQnS1hjYhCiEMB3z2WwWZTxqiOD5Z8D3EPvC1Km0kSgVrUSC2mq2954ZKeVbUWDxPx4FFvlKi7pA6L4pEXOZFzEZ37+Ttyxz319ol8IIZiMiXj84hoA4OWHBi+gc4XeB9QUuha0SzKrQ9dUd7N5LoWK/wXRgwSbunhiIYuwEKjb+cgZTV59ZKbjz/uJuAhKget3pgayPJYH9D7A7pPbyxCTkfpZL2YKkv8F0YMEGwh14somphPhgfM1OYMDq0V/5QDaLQAP6H1BrRtUT4qayhYBLZPf0EOXhlyh692i2ZLM57hwugqrdHnlocEqV2TwgN4H1FkuUhWRUACBgKZEk+HGW4uKleFW6NMJ0SjJ5HNcON1k51gEMTGIO/Z3zpvfToZX1g0QybAAQqwK3Tw50T4awM6wJ0UJIZifiOJCusAVOqervPeeQ3jb7buNBR+DRlOFTgj5S0LIKiHklMvt9xBCsoSQZ/R/H+z8YQ43gQBBIiyYFLpVcSeb7BUd9rJFoGa7cIXO6SYTcRFHB2hcrh0vsu5/AvgYgM82uM+jlNI3duSIRhTzPJdiXUAXjE5SwTYoiFI69AodqCVGp/imIg7HlaYKnVL6CID1bTiWkcZcmliQqoiaLRe9CoZ1kpqpVFWoFEOv0NmQLj8bdzicUaNTSdG7CSEnCCEPEUJucLsTIeQBQshxQsjxdDrdoV89HJgHdJUkBbGQVaEDzvNcWBPSsNdmG5YL7xLlcFzpREB/GsA+SuktAP4MwD+53ZFS+klK6R2U0jtmZgazLKhbaJaLOSlaC+isCibr0C3K9okO4/o5M684PI23HNuFW/aM9fpQOJy+pe2ATindopTm9Y+/CiBECBnMqvweolWy1DpFo6JHhc5G5w5x2SKgbav/yNtvNWZ4cDicetoO6ISQHURv3SOE3Kn/zLV2f+6oYU2KVq0KvcHERWPU7pBbLhwOpzlNowAh5K8B3ANgmhCyAOC/AAgBAKX0QQA/AeAXCSFVACUAb6dsRTvHM6loCPlKFapK6+vQI+6LopnlMuwKncPhNKdpQKeU/nST2z8GrayR0wapiABKgbxURUlys1waKPQh99A5HE5zeOt/n8BU+FpeQlWlFsWdbLBX1FDoQ162yOFwmsMDep/AVmktZUsAYKlDF4IBxMQgV+gcDqchPKD3CUyhr25pW3nsw7aSEcExKVrQJzMOex06h8NpDg/ofQKrNV/eKgOoD+j7puJ47MIaZH2RNKNQqYIQIBLif0oOZ9ThUaBPYD75clYL6FHbtLf/61UHcWW9hC88vWD5eqGiDPWCaA6H4x0e0PsEZrmsGArdaqG85ugsbtk9hj/7l/OQqjWVbq9Z53A4owsP6H2CodD1gB61BWlCCN7/+iNY2Cjh8yaVXpCUod5WxOFwvMMDep8gBAOIi0HXpCgA3HNkBsf2jONj/3IehUoVf/6vL+LbZ1Ywl+IDqzgcDg/ofUUqGjJZLvUBnRCCX3/9ESxulvDyP/wX/P5Dz+OHr5nGH/3ELdt9qBwOpw/h1+p9RCoSwlLW2UNnvPLwNF55eBoLGyX86U8dw71HZ7fzEDkcTh/DA3ofwXx0wFmhA5pK/8zP3QlCwCtbOByOBR7Q+whWiw7Uly2aCQR4IOdwOPVwD72PSOkKPRIK8KDN4XB8wwN6H8EUOp/LwuFwWoEH9D6CNRc1sls4HA7HDR7Q+wiWFOWdnxwOpxV4QO8japYLD+gcDsc/PKD3EYblwgM6h8NpAR7Q+wi25CLOk6IcDqcFeEDvI7hC53A47cADeh/Bk6IcDqcdeEDvI3gdOofDaQce0PsIbrlwOJx24FKwjxCFAH77/qN49RE+QZHD4fiHB/Q+44FXXdPrQ+BwOAMKt1w4HA5nSOABncPhcIYEHtA5HA5nSGga0Akhf0kIWSWEnHK5nRBCPkoIOU8IOUkIua3zh8nhcDicZnhR6P8TwH0Nbn8DgMP6vwcAfKL9w+JwOByOX5oGdErpIwDWG9zlzQA+SzUeBzBOCNnZqQPkcDgcjjc64aHPA7hi+nxB/1odhJAHCCHHCSHH0+l0B341h8PhcBidCOhOyy+p0x0ppZ+klN5BKb1jZmamA7+aw+FwOIxONBYtANhj+nw3gKvNvumpp57KEEIut/g7pwFkWvzeYYSfDyv8fFjh58PKoJ+PfW43dCKgfxHALxFC/gbAywBkKaVLzb6JUtqyRCeEHKeU3tHq9w8b/HxY4efDCj8fVob5fDQN6ISQvwZwD4BpQsgCgP8CIAQAlNIHAXwVwP0AzgMoAvi5bh0sh8PhcNxpGtAppT/d5HYK4H0dOyIOh8PhtMSgdop+stcH0Gfw82GFnw8r/HxYGdrzQTSBzeFwOJxBZ1AVOofD4XBs8IDO4XA4Q0JfBHRCyH2EkLP6gK/fcrjddQCY2/cSQiYJId8khJzT/5/YrsfTDl06F39ECHlev/8/EkLGt+nhtE03zofp9v9ICKGEkOluP45O0a3zQQj5Zf2204SQD2/HY+kEXXq9HCOEPE4IeUbvbL9zux5P21BKe/oPQBDAiwAOAhABnABwve0+9wN4CFpX6l0Anmj2vQA+DOC39I9/C8Af9vqx9vBc/AgAQf/4DwfhXHTzfOi37wHwdQCXAUz3+rH2+PlxL4BvAQjrn8/2+rH2+Hx8A8AbTN//nV4/Vq//+kGh3wngPKX0AqVUAvA30AZ+mXEbANboe98M4DP6x58B8JYuP45O0JVzQSn9BqW0qn//49C6eQeBbj03AOBPAfwnuIyp6FO6dT5+EcAfUEorAEApXd2OB9MBunU+KICU/vEYPHS+9wv9ENC9DPdyu0+j752jeseq/v8gbF7u1rkw8/PQFMsg0JXzQQh5E4BFSumJTh9wl+nW8+MIgFcSQp4ghPwrIeSHOnrU3aNb5+P9AP6IEHIFwB8D+M+dO+Tu0g8B3ctwL7f7eB4MNiB09VwQQn4HQBXA51o6uu2n4+eDEBID8DsAPtjmsfWCbj0/BAAT0CyJ3wDwd4QQp/v3G906H78I4NcopXsA/BqAT7V8hNtMPwR0L8O93O7T6HtX9Esr6P8PwmVkt84FCCHvBvBGAD9DdXNwAOjG+bgGwAEAJwghl/SvP00I2dHRI+8O3Xp+LAD4gm5LPAlAhTbAqt/p1vl4N4Av6B//PTR7ZjDotYkPTR1cgPYiY8mJG2z3+VFYExtPNvteAH8Ea1L0w71+rD08F/cBeA7ATK8fYz+cD9v3X8LgJEW79fx4D4Df1T8+As2KIL1+vD08H2cA3KN//FoAT/X6sXo+J70+AP2k3Q/gBWhZ598xPcneo39MAHxcv/1ZAHc0+l7961MAvg3gnP7/ZK8fZw/PxXn9RfqM/u/BXj/OXp4P288fmIDexeeHCOB/AzgF4GkAr+n14+zx+XgFgKegBfknANze68fp9R9v/edwOJwhoR88dA6Hw+F0AB7QORwOZ0jgAZ3D4XCGBB7QORwOZ0jgAZ3D4XCGBB7QORwOZ0jgAZ3D4XCGhP8fGAXkCZmtuo8AAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "from thinkdsp import Wave\n", "\n", "wave = Wave(total)\n", "wave.plot()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Here's the whole process in a function:" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [], "source": [ "def voss(nrows, ncols=16):\n", " \"\"\"Generates pink noise using the Voss-McCartney algorithm.\n", " \n", " nrows: number of values to generate\n", " rcols: number of random sources to add\n", " \n", " returns: NumPy array\n", " \"\"\"\n", " array = np.full((nrows, ncols), np.nan)\n", " array[0, :] = np.random.random(ncols)\n", " array[:, 0] = np.random.random(nrows)\n", " \n", " # the total number of changes is nrows\n", " n = nrows\n", " cols = np.random.geometric(0.5, n)\n", " cols[cols >= ncols] = 0\n", " rows = np.random.randint(nrows, size=n)\n", " array[rows, cols] = np.random.random(n)\n", "\n", " df = pd.DataFrame(array)\n", " df.fillna(method='ffill', axis=0, inplace=True)\n", " total = df.sum(axis=1)\n", "\n", " return total.values" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To test it I'll generate 11025 values:" ] }, { "cell_type": "code", "execution_count": 12, "metadata": { "jupyter": { "outputs_hidden": false } }, "outputs": [ { "data": { "text/plain": [ "array([5.3955119 , 5.49020141, 4.82135171, ..., 5.74383027, 6.24352434,\n", " 5.75793234])" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "ys = voss(11025)\n", "ys" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And make them into a Wave:" ] }, { "cell_type": "code", "execution_count": 13, "metadata": { "jupyter": { "outputs_hidden": false } }, "outputs": [], "source": [ "wave = Wave(ys)\n", "wave.unbias()\n", "wave.normalize()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Here's what it looks like:" ] }, { "cell_type": "code", "execution_count": 14, "metadata": { "jupyter": { "outputs_hidden": false } }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAABIyUlEQVR4nO2dd5gURfrHv+8mcmbJrEsSERDEJQuCoIKch3reHUbM5xlO72c49Mzh5PS8ZEIMZ9Y7A4qSVEAQyUjOCyywgLCw5LTsbv3+mJ7d3p7u6VQdZub9PM8+O9Oh+p3q7nqr3nrrfUkIAYZhGIYJG2lBC8AwDMMwerCCYhiGYUIJKyiGYRgmlLCCYhiGYUIJKyiGYRgmlGQELYATGjduLHJzc4MWg2EYhpHAkiVL9gohsrXbE1JB5ebmYvHixUGLwTAMw0iAiLbqbWcTH8MwDBNKWEExDMMwoYQVFMMwDBNKWEExDMMwoYQVFMMwDBNKWEExDMMwoYQVFMMwDBNKWEExDMMwoYQVFMMwjAecKivHc1PX4cjJ0qBFSVhYQTEMw3jA5z8V4pXvN+GFb9YHLUrCwgqKYRjGA0rKItnKS0rLA5YkcWEFxTAMw4QSVlAMwzAeIoIWIIFhBcUwDKOw88BxzN20V0pZczYWAYg18W3ddxRLt+2Xco1kJyHTbTAMw3jBBX+fhaMlZSgYO8J1WdNW7wYAbC8+VmX7ec9/DwBSrpHs8AiKYRhG4WhJWdAiMCpYQTEMUwUhBPo9Ox2fLikMWhRfmbFutyfllpSV41SZ9558E5YWYtaGIs+v4yesoBgmQckdMwl/92CNTWm5wM6DJzDmsxXSyw4zb80p8KTcpdsO4Jcv/ehJ2Wr++N/lGP3WworvCzbvw7Z9x+KcEcuzU9Yi7+lvZYvmGCkKioiGEdF6IsonojE6++8nomXK3yoiKiOihsq+AiJaqezjPO4MY4N/z8gPWoSkgci7stfuOuRd4Qb8dvx8DHx+pq1zXpu1GXuPlHgkkX1cO0kQUTqAlwFcAKAQwCIimiiEWBM9RgjxPIDnleMvAfBHIUSxqpjBQgg5rjMMw0iB3aOZoJExguoFIF8IsVkIUQLgYwAj4xx/JYCPJFyXYVIWIbxTHx4OJBjGFjIUVEsA21XfC5VtMRBRTQDDAHym2iwAfENES4joVqOLENGtRLSYiBYXFSXXRCCTGpSWlWPvkZNBi5HSCCHQ5sFJeHPOlph95KWNj3GEDAWld1eNuneXAPhRY97rL4ToAWA4gDuIaKDeiUKI8UKIPCFEXnZ2tjuJGSYAnvp6DfKe/g6HT5xyVc6eQycw4Dl7cwtMJUJE7oWWRFFP2/Yds+38kKjIUFCFAFqrvrcCsNPg2FHQmPeEEDuV/3sATEDEZMgwScfU1T8DAI6etLfWZseB45itch/+ctlOFO4/LlU2PcrKBUa+/CPKypNnNspDy6hvDHx+pm3nBzXbi48hd8wkzN+8T6JU3iBDQS0C0IGI2hBRFiJKaKL2ICKqB+A8AF+qttUiojrRzwAuBLBKgkwMEzqcNo4X/H0WrlO5D/tpiVq+/QCOnEiefEZJoJ9cM09RTImwzs21F58QopSI7gQwDUA6gLeEEKuJ6DZl/zjl0MsAfCOEOKo6vSmACYrtNwPAh0KIqW5lYpgwY1fBHAs4ukFpOaeL8ILDJ06hTvXMoMUINVJi8QkhJgOYrNk2TvP9bQBva7ZtBtBNhgwMw3hDMjkPeOn9GKVw/zE0qlXN9LhTZTyeM4ODxTKMTyRKc5QocjrB7m87caoMB46dQrN61S2fc+5fZ6J3m4amx3mp9mfHC3mUQDeYQx35SP6ew/hquZH/CMNYw+8RjR+jjqApLxe6cezu+OAn9Hl2uu3yFmwpNj3Gq9tYVi6wfPsB8+t7c3mpsILykaF/n427PloatBhMiBBC4LVZm7Dn8And/ZNX7orZptewfLlsB3o89S1KfQhKCkQa9Fe+z8fB4+5c5v3GSNeWGeyYvm6PZ7KQAxWhjrVnxNgpa/HCtxuciBQ6WEExjM9c9fr8iths634+jGenrMMfDDou9/5vecy2d+cVxGx7bOJqFB8twWGfPO5mbSjCc1PX44mJq325niyEgX1Lq59OnApn2g0r0cq/XWMtKrvZuPjV7zfhM5WnX2lZue/1wgqKYSQjhMDYKeuwYfdhzfbI/01FR/G40rCXKhPlR07qKxY9M1BBnEWaJWXlOGpQllP0TIonlSyxRnKHFavWyid1FvJKJ+Q2tr9OXYd7P6nsIF3+6lyc8Yi/TtasoBhGMsVHSzBu1iZc9fp8w2OsLn6162J++Stz0fmxabbOMSMV5qC07D6ob3JNBKzOUdrVjysKD9oXxiWsoBhGMtHmPJ4OWrx1f9VzJOmAHQfiR5goKS3HRwu3oTyOcFpZDIxitmVLJE6VCxxyGZIqKEI+MLMFKyjGU2as2433dOZMUoHYhiK2Ufd7idHLM/Px4Ocr8eXyHZbPCcMAKnfMJDz19Rr8d9E2LN223/wEA1bvrJqXSQihO5qdvaEIZz3+jePrWMJlvQohYszIyQavg2I8IdpDv/HtSA7Ka/vmBiiNv9hp0P1u/PcdjURTt+NMoXUsuPvjpZaisn+9YicOnyjFlb1y7AlpgDoCecHYEY7KWKNJHPjEV2vw9twCrH96mCvZnGDksGGVN+dswdOT1uKz3/fDOac18O26fsIKivGEjo9MQW6jWq7Luegfs1G4/xhWP+l/A+KU9T9HerV+jY7W/XwIp0o9dC/XtGdfLou/lu+zJYVVJtdlKSin7Dl8And+uBSvXt0jZt/bcwsABDNKdHvNlTsic0Lbi49VVVAWn7tECBDCJr4QU1YuMGnFLmmT1J//VIh5m/yJYHyqTGDjniOuy1m/+zCOBhyLLsr+oyW45+Olpp5r17y5AECsg4PebZTRSAz75w+hqSMAVZSTDNw+/+/MLcDCLcX4aOG2cNgrJVGwNxLWtNzhb9Ke1u2JSpNmWCLYs4IKMf/5cQvu+PAnfPaT9fmCePzf/5bjyjieZUx8Xp6Zjy+W7cRHC7ZZOr5Z3arhcUp8WkRrhXh6UWsCitdUJUIvXI3TZnfJVvPIEHZxqwKWK151dvWT0QJh9aLrF2dsdCyXTFhBhZg9hyN2/n2chTVUWLXhb957tMr3ePM+Ye7Yx5PND7nDcI1/fCu/wZZlGbE7goo+v58sKUTumEm6x4Ql8zMrqAA4ZbEnnYrrT8JMoo0WZJFIk+r5ew7jjg9+MnzHjF4ps98oow6em7rOMKSVW6rkC5NQXlpIHnZWUAHw/XrzcCWyOHziFCatiI3nxoQLK+3BGo2LtBOc9HmC7ifZufy9n6zApJW7sGqH3EWlRnWw78hJyx3JV77fhAc+XVFZpgzBFNTRyzcVHY1zpDXCoZ5YQSU9932yHHd8+BPy9+ivlzheUoYXp2/0Lcgo45xtxcYhjrwk6DmoZdtj1z0dPGZ9Ea1afxgpEzMdo7e/YO9RnPP0d1Xc3804earyPQtK8Vt51d+Zt7Xi88iXf8S4WZs8lMgYVlAJgJtGoHB/JLLAiVP6T+W/Z2zEC99uwCcJkP45LHjRsFgp0y+rS9AjJi16if2cOvsY/TQnP3mr0mGwEsC18jrBV+5DE1baOn759gMYO2WdR9LEhxWUh+TvOYKBz81E8dGSKtvtzi25aTDMzj2muEyHNXpzmIjGOJPVxMQLN+QVFVe0oe2CngvVu7x2wW3MOarPVn6qk98YPcdpfq4wKKuwwwt1PeRv09ZjW/Ex/PKlOZjzp/NtnfvuvAK8/oN104ERifgK7DxwHFkZaWhc2zxtdiLT9qHJqJ4Zvj7iMk2yu6BHVFYbciHiJ+r72zcbkNuopjQZXFdLiF7OoDshRoTv7Ugiio9FRk6F+49joypmlpVH4dEvK/PsuDHtVPby9Pf7nZ3VCv3GzkDe098FLYYvGJle9bBzp6x6iuoxanxirpWbuupn02OMUpWYvZPx2m8798WLBbBmJRYdDofLuBOkKCgiGkZE64kon4jG6OwfREQHiWiZ8veo1XMThdKy8opI0qPGz4usylY9OftUZr4geitOsnfKZHFBcWjWVjglaFVu56l54it5iQRXSvaIs43FH65eZ1bVMcL83BKTUFGy3thFBZUOHy/PzJdUanx6PvMd5ubv9eVasnGtoIgoHcDLAIYDOBPAlUR0ps6hPwghuit/T9o8N/SMnbIO/cfOwJ5DJzB/czEOHj9VxSygnuOR3YlaVFCM3DGTpLvWyuSKcfNw6cs/Bi2G7+w6GD/9hVfMXBd/4r6k1Hp21OlrvUt7bgU/unOmysIDIdSecl5z1RsL4u5ftcP9EgYvkDGC6gUgXwixWQhRAuBjACN9ODdUzN4YaRCiZj0talOa7GH+LGVd1Yx1xg2JmSXPj0Fd1KMw0bFTV1bXvFkpUuYI7qmv10jJjjpttbX04mbsOXQCff4yHZuKYuM3Wq1vN04HWkcmLXuPBjv633vkJI6VeJe92Muy3SBDQbUEsF31vVDZpqUvES0noilE1NnmuSCiW4loMREtLiryb6GrXdSmNPVwXo2fBr7oy61WUI99uQoz1slpWFKKoG18NjBLXBgEP23bj0tf/hEnS2NHbpNX7sLPh07gXSW6uN+Y9Rk3axa/HjlZiqkrzee8ZJH39HcY8e85vl0vLMhQUHqvrfZ2/wTgNCFENwAvAvjCxrmRjUKMF0LkCSHysrOzncrqGWa9PPUP9WoOSq8y1yvOGU9MXFOx7Z15WyvyNFXIJFGOfAlRzMOMF+7BVnSfDIcWJ4/eZz8V4qhJBHcrPDxhFZZtP4CNu+09H9a9+KqeFcULN/Mxn63AfxdvNz/QJR8v3Fbhmbhlr/sIEYmGDAVVCKC16nsrAFUSxgghDgkhjiifJwPIJKLGVs5NFKKPt7WXQfa1zQuct9mfNBubio5g6N9n+XKtRMCqSgmnk28lXi/UjPf73b4vVsyQdhI4ArEj1Pw9R5A7ZpIlT0I7jPl8JUb6MHcb1udPxjqoRQA6EFEbADsAjAJwlfoAImoGYLcQQhBRL0QU4z4AB8zOTTSMGiS14nKavyUR2H3Qm2CYYcIvj8iw5OQBYJoDywqW5tl0enhua8HKiF6bCt6IRQXFOFVaHhNMdeWOAwCAqavsx73M33MEDWpmolFA6/4OHCsJfK2bEa4VlBCilIjuBDANQDqAt4QQq4noNmX/OABXAPg9EZUCOA5glIiMqXXPdStTENgxEchud6INppNiZS+DClGb6hl+RQDYoXEq8UMtGseqM/7NBXuPIrexu+zJ0q0KtsuzdsKvx80DAPTMtZ5i3Yyhf5+FOtUzsPLxi5wV4LLu7I4e/URKJAnFbDdZs22c6vNLAF6yem4iYsfEt+fwCWwqOoJ22bUtlW3WY49e011IJDktRDKHbwl6LVlYef2HzXjmsq647b0l6N22oeFxTmvPn3WD5tItKqhMWqgdQc3Nj5jQv1jmbIYiqiRmbShCjcx09GpjXI9u0AuvJQQ8SwPiFo4kIQvlvhfsNViprnounpu6HkNekDdPE31VtMphzyHzh+4/PxYAAJ6etFaKLD9sTMwFgVZYsCXSCO055I/LsZ6y37rP+kS5k0SXjtJxKP+nrv4ZT3y1xvA4s/h5Rhy3mM7enRozP1sd3UWtoASAmZJS6Ix+ayF+89o8KWXpcULHg1JA4B/fbrBVzoLN+3xZd8kKyibb9h3TDSOzRWk4bnlvccw+v/jndxvxyBerMHllxA7+zGQ5SscO784rkFaWOgV17phJWP+zfsoQv1i67QAAYLGN9N8yTahEwEQbPXT9OSzjhviduQX48xf2Il0D3q+he3Sifau/XZGs/Ia1KgWrdToKYcQwXU6Vxv5QJ2b5346fj1+8OMfz5QysoGxQdPgkBj4/E09/HdtLjD7gxhk7/eG9+Vtx+wc/+XQ1+wgh8K/vNmKzzoJMLd2e+KbK9wlLd3glliFLtu5H7phJ+G6Ns3VjpS4m5bTPkhD2niO7V35s4mp8tFDfdVrm86v3jsQr32osuSAn+tMSREG9MWdzzDY3JtSrHaY9sQorKA1CCDw7Za1u9tIDSpSIOQ7iWnlpR59i4NoaxDtjNk+z72gJ/vHdBlxjEnpFt+wAftAHCyLhaG5+t3JkbGcu6s8TVlk6Tu/5cPvEePHIWTW3qckdMwnjZ1cmvNMzXZaVR6wS6nucv+cwuj/5Tcyxsug/doZKJneE1QtOywGdRI9uRDcKvisLVlAaDhw7hddmbcb1/1kotdxSnaRrP2yUY7feaOBGeyqELnXRF/mkSXBOPYLopcpwjJDVOSGy1xDqpohwKEo0bp9eKCIrTdz42bE9dzV/mRxZZ6WexH9zToFug2qEXQedMEbbcIJbx6QwK1dWUBqia5RKywUe+WIVvlpeafN304NX98CjXPumXCWoZdIKe2syZMTjMqujKQ7WiVSU7dGY8L15BdhukE5dTynabRCWxslRFEVv/c+nS9xHKjh04pSr1BtR1sWZ/7PSwKnnw+Idr+5T2e2QbFGFI7Lb6Lpdc+bF6N4Lq4v+sxteDcUKSoP6Vr03fyvu+mhp5b6A7qNfpq3h//rB82uoPaHs4kU9HDlZike+XG2YA0nGNfVGz1r0GqPPllSdcyOypxyFAM56/Bv8/v1wzUnG+wVqpWS37t+Ys6Xi8yc2wxC5VQZedJ7W75bvFKT3M8vKw5kXDmAFZUi82xXWm6nm79+st33OVo/tyW7xooMQHTEfOq5vStJreFrUqyFfEIs4qYPv1soLDKz36FuJjKJ+Z+IpA/Vx2rVGdvhkSaHjc4NCG60jK52bZ66BJOXfM/xJhuYnXiwCNmtb03TekPM6hi9YsR5Rt3g1Hy9yZja0apazVJbF4/zsBrrt/PxsYc2hGSNfqhqtvFpmuusyrUAUG7UkLEiJJJEslJcLrCyMLD6TYbcPO7IdQexg1B7MjJPTylMMW8PYHX5ElDhxqiym0SOQLRX9zRrzwKVCCBBRlTU+RuSOmYQRXZvrlGFDKBvH+2mpCMMszCZNSg/ZTkH74+S8Kglpe8cjKBVvzNmMG95eBAA4cSqcN0wmVpPp+cXuQycq6t93DFoovTbSj3bzD6q5zygCosId2y5fr9Bf4HvLu0uQO2aS5fnHSStjnVw++8ncnGY5qntAE70yguE6QUbad6tVdvZT3+o+5mGesGAFpWLtLtWkpK73VgSZ+Y4WFxRjg2YydL5PqTHChtUU5F5w2KCB0nt5/Xih9UJGbdx9BC/P3KRztD7qhmuWQWfE6vyUW/PqPlXvXe0Zq0XWHJSfGHmAWuH5afbnit1gd5F00LCJL2CuUKIjF4wdUbFttwR7diJi1hMMonMd1AhKj2dt5mQ6oHL8CNNC1EUFxRXBkgv3H6tStnokkyD6KcESCcbeSCvm3aDgEZSKsLwP6YkSN8UFQZly7KI33+THHNRxCaPJ2RsqR01uq1vm7Vq98xCGvDALuw4ex7l/nYkBz82s2PfpksIKJZX8b0F8dnsQlFjvPt798TLp15EFKygb+OX2mZHACspLyYNQabq3wuaPDIMydmuiM4u6kDtmku0y9x/Vd+0/oqSeSJQRlFd4EedOxqPopymeFZQBeu9GRro/b0yGnm+zASUOQgZ5SeM6wWQF9Qo9TzK7T0Eo1s0FryNtY1ZvO5MkVJEe5eUCRy3GPTxgsIZPDxneyXN8TKnDCsoAvVhxIegIx1Dq0KvLK67unWPpODdVuXBLMUa+NEeOcnYgiF2FEx1BNaqVZf9iGq7s1drReSF8dF2jHxcwOWj7kPUcrmNtzE1+HkBGADewggohYehwO8VLz6vdByPOIw9+vgLLCw9iW7HcyWk95xRdJ4k4ZRyMEy1axgjczui6igxh7F0xjAmsoNSYtB9hfMfDlobcjTnLrHq1SeJk86ROni9dJ4k4P3GUzryBlVh8XhO8BAxjH1ZQNgjjauuwjbaSraeuV7+zNhTh9IenxCzu/Nu09bouu/f8N3bRLRO+Z5exhp9vuBQFRUTDiGg9EeUT0Rid/VcT0Qrlby4RdVPtKyCilUS0jIiCy5dugb8EkEI9lTBTbn7rvglLC3UjMHy5bCdKSsuxUbPA+iWDqAB7jxiHmLHLoRPWJ8TVfLlsZ2DRErzioA3nAEYefnZCXSsoIkoH8DKA4QDOBHAlEZ2pOWwLgPOEEGcBeArAeM3+wUKI7kKIPLfyeMk6Hxa0vTVnC256J9R6Oi5WTXx2EtH5hVbyP/53edy1KG7SuWuJl29JjRvHgAc/X+n43CAwe5Lu/JBHpsmOjBFULwD5QojNQogSAB8DGKk+QAgxVwixX/k6H0ArCdf1lW/X7MbOg95HeNCbBzFi18HjmGqQ7j0R+H591cCwZsotbCahUza8CH/72jwp11y1w3knKWxu2dPjhFkqPlqC1TvddQgHPDfD/CDGNn4aMmSEOmoJQB3DvxBA7zjH3wRgiuq7APANEQkArwkhtKMrAAAR3QrgVgDIybHmymyXzUXGXmG36GTElcl9nyy37Tb9m9fmYXtxuBodOzpk3c+HMahjk4rvYTLxlUscHQHAgi3FaFrX2hqxufn+rTMJkr99s0F3OxFwyYtzXKdkD9u7kSz87r0lvl1LhoLSa5N0324iGoyIgjpXtbm/EGInETUB8C0RrRNCzI4pMKK4xgNAXl6eJ02Vdk7BjFe+z8ecjXvx4S19XF/70zgJ1oxGFjsPJHbMvtIQOJ0YRVgYN9t6UFYAWG4hrbtVrnpjgbSy1IRsABoXt8qJSQ5kKKhCAOrVg60AxMwsE9FZAN4AMFwIUeEvLITYqfzfQ0QTEDEZxigoP7C7hue5qZFIxB8u2OaFOKaEscGxU4XaORyrvQ71cduLj6F+zUzUqZ5p/cIGqDsCP209YOvckS//aHqMF7HVGCaZkTEHtQhAByJqQ0RZAEYBmKg+gIhyAHwO4FohxAbV9lpEVCf6GcCFAFZJkMkRRo3rzyZzTw9N8Hbyeeqq2Bw8yYDd9UGHdTzYBjw3E796da4skVRYkC2MPQSGSSJcj6CEEKVEdCeAaQDSAbwlhFhNRLcp+8cBeBRAIwCvKL3UUsVjrymACcq2DAAfCiGmupXJKUamNKeuvbJYVLC/yvejJ0vR+bFpAUkTHztTN9oRlFl7bxSbbMNuOSFvWN+Eg30S3fKZxEZKPighxGQAkzXbxqk+3wzgZp3zNgPopt0eFGHzEjPCyxf40IlTqCvBXOaEIJb4Jtm64qTgurcWWj42mrKeSU44koSKRHnMyzxsVR9KsLUyTjFL0malin///k+O0kwERaK043uPWJ+r08s8zCQPrKBU7A/h4lE9Bv/te8/K3n8s+c0rXy3fieH/+gFTVrqb20u0SAbJOFoMW7oZRi6c8t0CyRYiJqz41YBuUJYTbNxzBPVqVpozgx5h7Dl8Ar2eme5Z+VbzCzFMWOARlAUuf8ULLzF7vPp9ZF2On9ksnZAgVqQKtu07FrQIFbiNnGBGsgXyBThKe7LDCipB+OvUSFKyMZ+t8PQ6P+Y7T2lx6MQp7D5sffGw0zTkMpVgmBq4TIe5nlKZZFS6TCVs4kswVnncy3bDkBdmoeiw94tRZTZJRvHpgmj2ZCQ0TDVYPSU33GVjpOGHcpLNizMqU2Q4VQ+yevHHQ26+DSMPfr4S01YnbsBkJj6soJgYvl1jHGXaW6w19H6MM+woHaupMsy44T+LpJSTShQfLZEavLRHTn1pZTHuYQXFxPDK9/qJ98KG7OkHpws+9ZIaMolJehqbWcMEK6gEYseB475MCkvONGHIQc26M6s/LWyRA16eaS/yORNe7AaMZryFFVQC0X9sciVg276/qou3FcW4TxVlwEs96peS9hN2eDNHtoKqWz0DHZvWkVpmKsEKyiHabLB+4UsbY7Ml+3DBNkf1QZrZJCtu53+ZvE53DmrhlmLb149HorXlo/ueZnrMepv5zlIR2Sa+tDTC+Z2amB/I6JLyCmrngeOOwqUsKpDbIIYJu6OHhyasxPUOJvi1nVUrelGtxMpVJ/zGRkr1jxdFEkBrQxWpxUm09TWPXtI5aBGSguqZ6VLLI/DI1Q0pqaCEEPh2zW4cPVmKfmNnOFr8GtS8Q7y09LIot/FGvTd/q+/X3aXk53rF4T2IusNr10DtMsn7FWZ4cl8O2XWypJeZaJ2dMJGSCup/i7fjlncX4+25BQCA7zcUBStQAvPIF/7nl4zGRpy+1p07vHYEN2+z8ygaTDj5aGEw2a6jhM2hJ9FISQUV7SnvUHrQxUdLMGFpYZAihQohgMUFxXh99uagRYmL2+Cnk1caL/DkTi+w8KEhQYvgmgcDTh9DSLz5zDCRkgoqyiLVxPof/7s8QEnChQBwxbh5eGbyWpSWlWOWaoS5Ze9R5I6ZhI0eTLizUvCW6pn2Xne+He4hIpQno0uoT6S0gtq4R06q8GRDbTP/1/SNGP3WQvyYH0kMN1nJofTZTzs8uK70Ih3jNJBtmLHbTobpfiQyXI3OSWkFxeijbpi27I04ZUSznGYok/Fl5fITxdlxzvCKsnIBIURyNs42f1MY7kcywPXoHI5mzsTg1+hBO4Echmje7R6ajMvPbpmUCsrufU3CKvAdIv9Home1qocVhQf9vahHSBlBEdEwIlpPRPlENEZnPxHRv5X9K4ioh9VzvUC7QJSpipUXSoZ3khAC24sro0nUrZ4Z5+gIk1a4S9Nuhc+X7sDireFf5zbl7gG2jrdv4ks9FZUMoY5a1KsRtAjScK2giCgdwMsAhgM4E8CVRHSm5rDhADoof7cCeNXGuYzPnCqLNd9FFVI0GsFeCak1fti4FwOem4nNRZG5wE+WmHtSnrS5qLq8XOChCSttO3WcKgt/49ypeV1bx5fZ1FApqJ9ilh64pXOLur6b+JJp/lTGCKoXgHwhxGYhRAmAjwGM1BwzEsC7IsJ8APWJqLnFc6WTTDfQCwpUadDVNbVwSzE+V5wjPl+6AyslmRHWK+kq/j19o5Tyopw4VYbNe4/gwwXbcNv7kZQMOwwSFDKJze5DchZZyx5B3XV+B4zqmSO1TDOSyUIkYw6qJYDtqu+FAHpbOKalxXMBAER0KyKjL+Tk+HvDmQhb9lb1etwgydVcG3JIBqt2HMQvXpyD+y/qCADYVBRxj6+VJTeUTTKTSJP7D3xqPxqMHrIVVHoa4cwW9ka6bkkCK2UFMkZQetWhfbKNjrFybmSjEOOFEHlCiLzs7GybImqFSaI7GCBhbr6ik8TaILZuF/emClnpaQll4pslKRpMMjTuyfAboshQUIUAWqu+twKgzeBmdIyVc5kQQACKj2rzN8lpwWStY1Q7UERD0/EaSWe0blgjJeP71bHgqGMHoxp8c3Seq3K7tqzn6vxEQYaCWgSgAxG1IaIsAKMATNQcMxHAdYo3Xx8AB4UQuyyeywSJ0sAfPVmKOflFerskXEJOSR8sqAxcG+1FpqInmgya1q2O1g1r4olfplaU9NsHtZNantFoZkinpq7KrVvDeHYmmSxErhWUEKIUwJ0ApgFYC+B/QojVRHQbEd2mHDYZwGYA+QBeB3B7vHPdymTGoRPy5zySnTGfr8SP+eEOppqVUfk4R70OeQTljHbZtQEAo/vlBiuIz8hPtxGAskge/SRnoa4QYjIiSki9bZzqswBwh9Vzk4GuLeth5Y5wLZarlZUudw4mZI1/oryXV/XOwYcLgo2ybcYDwzoGLUJS4NV8ULxALonyHlghJUMd+WFaD6P5Xnbof1mmOS+tcMu2H/CucIeE8dnQkpmekk1DwpAqS2VS8in0wq1ZSxgfHyfzMfFeBFmKRd5cViVLt+2XVKp87LoyT7i9n6XjPrmtrxNxmAQk3ruXTDmoUjIW384DiZs5NUxIU8IeDKE+Wrjd/KCAsNN8tGlcC2fnNDDcP/3e85CljHYyJA7NkqiNCxSv6jHeG5NMty4lR1CpSlh7Vk9PWqsbXsku368vQuH+Y+YHJhBf3tk/7v522bXRumFNn6Rh7OKVk0RuI+N7HtLX3BGsoDwiFbybZf3Gk6XlmLrKOLutHX49bp6UcrzETkfBSgBdL7DSsC55eKgPkkQYcVZz364lgxqKN6BXyqJmlr7x68peOXjkF8kTzjQlFVQy9TDs4CR0Tby06D9u2utGnCrICquz62D4zbdePX9+j5Ab1a7m27US7ZXN0Yxqu7euL7V8I4X97OVd0djH++I1Kamg/CCMXjbHJIf58SP1RTKSCAspU7UTJ5toPb5/c298938DpZT5nxt6omduQyll6XFj/zaelW2XlFRQYZ2L8ZomdcLbszp8ojRoEXzDs4nzVLArJwjaDmrtahkVi5/DStO61VAwdgQevSQ8JsLUVFBBCxAQYY6tduJU6gRxDe9dqCRsMiZapzLaV1CPlr3+DY1qZcVsG9Wztc6RlZzVqjKm3/hr3cUH9ILUVFA+POthXOiYDNlCkwGvb4PdRIZe8K9R3aWWl2hPbnT8FO9ez75/sNRr1tBJJXPfRfEjgqiXJoTRGzR8rWiSoJ0kDQMZ6eF9zROth+wGr36rzHLdltVFcrTtEA/+dYmaW+OJnRPHVTx+4fqb9dbB2bH6hrGKU1JBhfFG+IHMhZyyCbFo0qlb3Zv18TLnoIxux5AzmliURZooAILpwLx9Q098cLNu/lTL+Cm2HRP+785rCwBoUDPWLBgmUlNBpVBvXU37JuGdpE2lO3LrQLkpHbR4WZdBvTpBXHZQxybo376xo3ODcFf552/P1pEjVpJhnZvhgYvOwJKHh6I+K6jUJIwNbpidJFIJdVoQL5ChRIzLsFq47CGU3OL8wz/Bu7bSMavq3IaureohPY1i1rGFsd+ekgpqQAdnvSI7hNHhN8zrb5JxVPvcFWd5VvanCRAYtnFteb3zMD+7ukS9+EIotlqmKp8N6viaPjk4o1kdj6XSJyUVVC2DMCEMI5Pf5MV38XXDaY1qGe6zMv/Tt22juPvddhiEAPJOk7eYNNEG/xVefJrtKx6/EOufHuaozEEds13JFMWushcCcQMWe0lKKig/CNP79IsEiGP22ETPEyknNPVqyI3J59zca802IFA16WGdaqnVKTRyWKlbPRPVMuxn7T1XNRcWnVey4rBidresOrM88cvO1g6UDCuoFGBYl2ZBi8C4RJvDLF4orbCYldqqIic8+6uursoKo8k8HpXroLy7GeOvc7aw1lAkg+09cxsiKyMNNXXWWXlNaiqokLzAfpHoEXD+cH77oEVIKLxUUG0thuvRPnODOjbB2MvdKSm/aJdtbD61wg39cys+O70VF5zZ1PQYK6NgvXdffZbRfFSUBQ8NwaVnt4w5zy9cKSgiakhE3xLRRuV/jKGSiFoT0UwiWktEq4nobtW+x4loBxEtU/4udiNPmFi981DQIiQ0D4/oVPH5mr6n2Tq3VEJuKbdker0oOqBOR3tFQV1xTqu4x2lHeBlpZHrOpd1buBNOEl/c0R9z/mQ/ykOzutXxqx6tcP9FHR13CqOKqW3jWCUp64ky6sDU1pmbb1q3uuo8/1WU2xHUGADThRAdAExXvmspBXCvEKITgD4A7iAidTTCfwghuit/k13KExo27jkStAiBcpaey6tf137im8CuHWVUzxxPyw9sUKy0UWYNcIv6Nap8r56ZbtrA9Wpj7LiRd5p/k/R1qmeiVQP7UR7q18zEC7/pViVXk902/Rzld2qrl6gyB1R6mvVmW88UXCU+oGp7msmILOFGUABGAnhH+fwOgEu1BwghdgkhflI+HwawFkBLl9dlnJBApk03vTXZaUWcYCZ+2+xaeNxF1Gi94rOVaPVmHnrxMPMUs3JXCsaOkJpo8ccx5+O3JkFPw4bTdDvR+tU6WaQR4elLu+CeoR0wwMbiYV0TXwK1A24VVFMhxC4googAxHUrIaJcAGcDWKDafCcRrSCit/RMhMlES02vkkksTrMRO81shDHj3kG43k3eHZ1GplWDmvjhgcGmAULjMe6aczDvwfMN90fdjb1wvDFq1FvWryHVvNSiXnXzgyRh16WbDEaoaQQ0qJWFe4aebjrScXK9sGKqoIjoOyJapfM30s6FiKg2gM8A3COEiE7QvAqgHYDuAHYBeCHO+bcS0WIiWlxUVGTn0oGhTfz1raSEZXbJVdbM9AhoLUOycHrTYBYr2qF1w5rIsGEC0lI9Mx3N6xl3pNo3qY2CsSMsTeLbxS9nHj/mUpz+lqhC057udFlA7eoZVVzUgapZDewo0CCUmemTLIQYKoToovP3JYDdRNQcAJT/e/TKIKJMRJTTB0KIz1Vl7xZClAkhygG8DqBXHDnGCyHyhBB52dnuFqxNWelPJtiwJP7q2qoevr9vEG5UeRcx9gl5Z9M1//ud/OgU4645B7cOjAQmTfb608Nuo240gnKqVOtWz0QzzYhRXZSdvkwiOklMBDBa+TwawJfaAyjyq94EsFYI8XfNPvUK0ssArHIpjyVmrk+MEZhMchvXSshwQq0a1AhNmJtkz6fVq438NOLDujTDQxdHPDJlmqbCzrhrzsEl3VrEOIs4Jd3msxcv3qO6pMa1w5tlG3CvoMYCuICINgK4QPkOImpBRFGPvP4ArgVwvo47+XNEtJKIVgAYDOCPLuXxnOddxFcLS0PrB9f2secariU6SRwmpTDUxKyVKOt8worXFr5aHi80Xffz4YrPXVrWw4tXnu3YNKedj7Nbzox7z8N7N+kbpNQRzDu3cJ/cco2HS2pcKSghxD4hxBAhRAflf7GyfacQ4mLl8xwhBAkhztK6kwshrhVCdFX2/TLqcBFmfu0ivloY2to5fxpsuh5FBm101nHYIdoDDEun+6Nb+uBXPeI7nw7v6l9IqbB0dt66Pg//ub5nxXe9tOOW8XASamT3FnjnxkiD3b+9PS/HV67ugZeuqprKYoRH9zpq5Yg18dkrp1WDmhjQQX8q5IzmlXOpF3W27uxiJMOP+XttyWaH1IwkkcK0alAT3VrX9/w6btNHRwONhmUE1bddI1MTaW0f483FDXVk4XxZqVfOP6NplQZPlklLNgQgL7ch1j01zLDhNuLirs3RqXlkpNGmcS0sf+xCXOyRgjK6LTJT5XRoUnm/iAjf/d9AvHhlbC4pLUYSrNkV0hEUExxuwrG08qERqV8zE0+NdB5gMtoAEwE1AogBlux4mQrEKWbjpxevPBuvXt2j4nt3Bx2t6pnuniWC/MC92vIBoFwIPDmyM/41qjsA4NLu8paOanVd+yZ1cEk351E8Jizd4VIiY1hBJShujCGDz2iCj27pI00WPQhUZUW9XaImjjQiX0cmZpwTJ6JBehrhzObubfpWcJsyRh3CxoxHfhHfG1WWZa61SfSGS7q1sGRG7ZFT3/a1m9Wtjil3D7B8vBcD+3svOL2Kie+6vrkY2b0lCsaOwGALkcut4tRZKghrBisoj/nDkA4Vn0NirQIA9Mz1fk3UpWe3xKMmjZsR0Qb0qt6RkEF2Fsl6iVnveVQvbyIetNWMmGv5qLRrV/NnBCurER7dLzdmW7xG+cmRnTHtjwMrzHh6eL1Gq1urerhrSIdKN/OEi9/uDaygPEY9sW5lYttyQxzy55coMqK48dw2eP+m3phwez9b59erkYmCsSNwg5toCwHQTFGsshWqzHxKrRrEmnjjdZ7MGuegOl6yXoHr+uZaN9vp/NZbBsh7RitDHVk7/u0bepoflMCwgpKIlQRiZqg9ouKhfn71Ih8Hjfo9PrdD47gZYK2QKClDLjizKd6/qXdMFBHX2NACZof+7rx2MducjnSBAO9NSB6KLi3lBUZupKxLambRBDuoozzTnxkjAkh8ygrKBDP7uxq916X4aIk8YdTXUr2cM+4b5Mk1ZGK3kx0mc6gao0ypUYgI53ZoLH1RqlEOHydERXvo4jMqtmWkJ2BTYFAReuY8Lx+nasqSCK0Z1gnDuzTDy1f1wO8HxXYigsaJU4pbEvCpTCz2HqlUUDIb3fJwdB4N0TYS6q9nNLMf0+5miWYUN0STt/mNOr+UndGOdv0OUGlqvnVgO/w45nzkndYAl8TpHXvxqEkxgRp0FvReM9lLK9TXGNKpKa7rexr+e6uLUFHKC0JEGHFWc086DJ2a13X07gUJKygDPri5N76+61xLx0ZvuuxeWo7LtUTxCFK//ceC3Vzb9lzXN9cbYWwysntLV+7zTvnnqEpFEy+MjZZfnNUCBWNHGO5vWb8GPv19vyrRBezStG51aSFz7hnawfQYJ9EPBp5eufbJTv1VYvzGZKan4cmRXSrSnYSVG/rnYuo9wQSsdgorKAP6t28sxbZsprTOO9140eDsB4yzeobdy0f7u9UOIlbcz8P8+4KQrGX9Ghil5ERyG0XC7ki+nUma9/Q0wpuj8+zJYLD9nqGnm54bzQ1ldB/0ft/eIycrPl/QyXkkds9MhR7Op13ZK+IJ299GHqlzTmuAAR2Mj3em5O3DCkqHUarkaFYeyKhX1Dkmrtt6Zd15fnsbklVi9Xk26tn6PsVj84J2g2P6icy2pH2T+I2/HmZVY7bOxW7NehFENh7nmzgbOXkyjpwsrficlkZSvSLDzjmnNUDB2BG28tF99vt+eO+m3lW2RR+rS7u3wIanh8sU0RBWUBp65jbA2F/ZW2XfsVkdzLxvEG4bWDmxaXXC1GljZ3ZedAK8QU1991mvJ8W1baT6uxXdk5CT9g6Yds9AjL/2HEvHylKMIdb9AIC3ru+J1yzUidHP0Bthuv3JjWpFOnojzopEXJBdh6kU6d0OqdEKuKB3W2u9xzaNa1V5yK7ubS2at9orrKOLhHjqyc/ebRrivNMjvdCwNEYhEUMKZp58dkhPI2SkW6udivBPLq8ZlkCzjolGW4i/W7NN22Oyd8kGtbKw8vELcY+y8N5NUkg10Y5sizhJIsOGnyZuVlAmdG7hbB5Kr0OkZ3pprJpYbVLX+iRruaaR/PqucyuiNYdBKcUzM4VAvIREHf7JFSG4Aerno061jLjzHTHnxtnnZZy8OtUzKzqh53V0lzQ1SnQhej0DS0eYCKJjwwrKI6zcyqyMNGk9p4z0tCquyE6cDKx41zklEZMlGiG7B2n1xa9YWhBi/eTkNn95Z/+Y+Q6n15hy9wA01enoaU8xc/wwI1OSCTqR3oqurSKd9eFdrKfocAsrKI+INsjxLAEdm9ZBjax0fHJbX+Q0rIk/DatcOGmWT8nMyhTdb6fXIyN5mRFVF5om0msZS1ABDGSZ+MLgH+nV89Cifg2cc1rDmIWu2lGn1YgtfuHmmXr3xl5434aCd0q77NooGDsCw7r4F1GCFZRHRN+H7q0rPfu0r+HKHQcBAD1zG2L2A4NtubXrjpDC0PIYIKMN8sq1NYgQLlWwWDdtlHBRTWxEIndxOU/584hOjs8tt2Dq7K3xPNQe2qBWFu4Y3A5/+3U3x3LY8YozQsZ7MfD0bJxrw0SaSKSOr6WHnK5ybnjo4jPwl8nrKmLPeeUuHZIwZI4IQwOppm51e/Z/v6r+gWEdq3y/fXB7dGtdv8qiUycEmQTyxzHno3pGGhrVroachjWxrfiY7nxtvOc7V4lC0bVlPSzbfkD3GCv36P6LzjA/KA5f3XUufj54wlUZTHxYQSlk16mGosMnsahgv63zJv9hADqpMoreMqAt+rZtXGGvjbbGdtddmHmKWW0kw2JNkzLB6pFmMKrrHjn18dO2A5aOlz1/l5FGuH1Q1TVy6WnkWjlFy7HLue0bo3qm+xGsetRRabK0J8+ADtmYes8AlJULvDd/q6VzvHgPGtbKQkM3Ke5RmTn6ws7OFw8HTbxgA25hBaVQdPik+UE6nKmZtyGiSuWEyhdDwNpLEvXOM1sXEU9/hdGN2O46KD8xmntrWCvcoWuc4mSu8f2b489xRO+pHeVXMU+q6/Ea+Z9rELPvjGZ1sXXfUcvXCisdm9WJG4oqEdBzSpGFqy4RETUkom+JaKPyXzeUAhEVENFKIlpGRIvtnp/IVOZ3sdb9L1UM7Bk6L3qeKptry/rx5yGcmADDqNj84Jo+RmvW9CtRtnnV71rv4GK9nRFdWtTD7wa2xYtXxganNSJePUb3xZP1tEa1cL1OckIgfGbkZEQdcccr3I7ZxwCYLoToAGC68t2IwUKI7kIIddAuO+cHjpPU42SyqFBLeXnkv9k8wZsh80Iyo8oIymHz4VV8vkT3KrTDHxyG1jIjLY3w4MWd0MokbbtsjLLwap+UIOfdkpWzc+p7fg23CmokgHeUz+8AuNTn833FSQQBu+600agC2hh657ZvXOWla1y7Ghpp7N/q/S2V+IDSE+cxurhp/vxUkNcZjDiC4OERnVC/ZqatBepOYfWUmLidg2oqhNgFAEKIXURkFOVRAPiGiASA14QQ422eDyK6FcCtAJCTk+NSbGe46b9b1W2nN62Dpy7tgos1i+E6t6iLRQXFVbbN+dP5KI0OuVQQVaZMDwvJZD7UG8nJHtt55SkoKy2GDIZ3bY7hXeO7+Mt6alJplOw3XnoUm46giOg7Ilql8zfSxnX6CyF6ABgO4A4isp2URAgxXgiRJ4TIy872zmskvgz2z6l0krB+8rV9TqtI/Rzl/os6xhxXIysddVQu0mc2j0x+hyW5n5owO0nYxQ8Xf5nx/hIZs1owrCfN5kR/5sKIH51O0xGUEGKo0T4i2k1EzZXRT3MAewzK2Kn830NEEwD0AjAbgKXzw0JUydw9pAP+NX2jpXOiN1HvPeqm8vYzIyM9Dc3r1wB03J6jNKiVZWvU1KpBDRTuP14pq4fPm4yi/W6zLbZ9lqiZlY5jJWWu5LHD9f1ycbHJ6CSVYP2UmLidg5oIYLTyeTSAL7UHEFEtIqoT/QzgQgCrrJ4fRv54gXlStShqN3Mtdhu6G5T5A1npq72aOJ7/4JCYbTJMLNV8SpIWxcr96dLSmsv2godi6ySKXs24jZb9+C87+57HyQvMnhrDexQn3QuTOLh948cCuICINgK4QPkOImpBRJOVY5oCmENEywEsBDBJCDE13vlBYOUBdtWDl9D7j66NkmeXl1SQhmb13IXiMeLhX5zpSbluaGRxrVQdm9Eq6noYlTupsPxesYZKRFw5SQgh9gGI6RoqJr2Llc+bAegGvDI630+a1a2Onw+dQLO61bHr4AnkNDR2k3Uy4oh3it3SEnn9TVVvRmdl2A1JFDY+uLk3murE0YvWx+lNa2PD7iMAgKGd4meVZSIYzu1qNnM+QO/w0vKe8sFir+0bWaQZVT7xfPtrVUt3fB0Znl/RCWFZL1s04rMfpjOvRmsywu8Y4SQkUDz6t28cN8V7tio3WKqbpBrXjiyhaGMxM7UZ7MWXmKS8gopyUedm+G1eazw8wtiM1NLBIsR4ThJ2iUZxltVw/rZnDgrGjkAtBwuQjehjkIFY3UDI9P6RFYxXbwGrURK96L28urec5Q7J5IIvi7zchnj/pt6478JY71U17OyY3KR8LL7R/XKxZe9R3D20g2k2Tid6IZ6ThF3KFA3lVW9QRqk9cryLVqU3CjWLWWiV/7uwIxYWFGP+5sq1Zmb1XLu6+evTp23DKmUy1knWFBJJgw/9qpRWUAv/PAS1q2W4ygljFRnrWqKBZGWn8KhTPQPFR0ukKD4rRcgUP9CxhwDOaF4HszYUIdtgAexb1/fE3sMl1ooT+p8Z9/AY1Tu8fFZTWkFlunTltYLMFyM6gpI9N/L+Tb0xbfXPrlMHeI3eixB0O37fhR1xQaemhskma2ZlIKdR/NeMp0ecY7VxDPo5YZyR0grKT5y+IP++8mzsVVKBlAlvFFTrhjVx84C2UsuMh5n0c8ecb70wBxWblZGGktLYEFG2ISAzPQ15uXLWG6kbW9lK6/ZB7aosymaYRIAVlMdURDN3qKF+2a1FxeeyMm8UlBvq18zEgWOnKr6bTfhbEb2FjVTavdo0xPR1NgOQWLwXRqLKjqrux918YJi77LFhxXLiTk+lSE38qFP24vMYmTcxugD2bEmRJGTwzKVd0baxNVfgxy45E1PuHuh4rkuvMerXvjHWPHmRzXLkKBjZ3ndquXgOyhocszB4vEqDA/AIyhbDNRHG3TC0U1Pce6H1kEkA0KVlPUy7ZyA6xFlLEwQT7zoXY6esxfvzt8U97gYl9Ud0Ls0uRo1RzSx7j7HDyzMhxOxWdmhSG4M6ZqOrwRwhE25YQdngFgfzNHqDhVpZ6XhjdF7sDgt0bCY/G6pbalfLQNM63oQ3MsNJD5p73alDywY18Oc4axuZcJPSJr76Ne2FzlGbplrUq45hnc1HVFpz1j9+2w2T/jDA1nUTATtNvlPDWP/2ctbF/O68doaLcANBVSFPjuwcnBwJCPc1gsOP6BwpPYJyU8FzdSJ2W+Gys1s5vqYffPb7fvjVq3MtH3+WJmWIl8+srGR7Aztko/iItbVJ9114ujQvPTOEcB/FPNUoN9BQXs6LMBqCTFjIpBbnnGY9EkTB2BForQTXtdOTjafEmtpM/+20B/3wLzrFbNNzerjz/A7o07aRlGsaob4ur4myR2Z6/CaMqzOxYQXFSMVtg1DLpsODk57yweMlttNf6CFbmXCf3z5DzmgS11zLdZrYsILygbbZtfDcFWcFLYYvWGkQiAgPDo+sy7m0ewtc2UsVdNWHLm80pQWT+KSlkW6AZw7A6z28DipJmHHvIPwmr3XQYoSKaEK+6pnpePbyrhXbvcryq8bqQmc2tyUGep6tPAflH5wPikkYnLbpWcpcQjOdhH7x+GHjXodXDA9qRXjBmU3Rsn4N3DygTXACJRHcx0hsUtqLzyp1q2dgdL/coMUINXZ7rFpHgw3PDMeEpYUYdLq9TLL5e+yb68xGaXcP6YB/Td+IvhrnCM8REU/FH+3EImSYJIYVlAVWPG4vlE5KY9Mupj7ciQv+8VNlutuv75eLt+cW6O7LMDHx9W7TEAVjR5heW1bvnHv5TCLihwmcTXyMFKKjEr8bWyOX73hp7GUF25Vte+d5E3mc1ao+gMrwWkxi4moERUQNAfwXQC6AAgC/EULs1xzTUTkmSlsAjwoh/klEjwO4BUCRsu8hIcRkNzIxwXDjuW2w6+Bx3DLQv7QdgLOwRVEF9cxlXdBNaciCxI8V+alG49rVLI2CGecM69IMU1b9jPsv6ujZNdyOoMYAmC6E6ABguvK9CkKI9UKI7kKI7gDOAXAMwATVIf+I7mfllLjUrpaBZy8/C7WrWevzOB0t9Mipb7F8Y6L64Orep1VJNBjdbjWYrCy1cnZOffyqRys8f4X3mZ0ZRhY1szLw+nV5ttLj2MWtghoJ4B3l8zsALjU5fgiATUKIrS6vy3jIH4dGoqzfMbgdCsaOwB2D20m/RnQeyG5oHy/d0KNmoQa14i/ilR1sNjM9DS/8phtyLaYtYZhUwa2CaiqE2AUAyn8zF6xRAD7SbLuTiFYQ0VtEZBhnh4huJaLFRLS4qKjI6DBGAhnpESUQHUl0aVE13l6ejXBIRlx2divc2L8N7rNpHtDqJyPzWDw1ZrTvvgtPx4Tb+6Gz5vdalYVhGLmY2mOI6DsAemG7/2znQkSUBeCXAB5UbX4VwFOIWGSeAvACgBv1zhdCjAcwHgDy8vJ4NtlDog1vdKCgrex3b+qFvYetBVs1IisjDY9e4j4Nwu/OczDnZaBZMtLTcHaOe+XLMIwcTBWUEGKo0T4i2k1EzYUQu4ioOYB4ubeHA/hJCLFbVXbFZyJ6HcDX1sRmvCStIk29UP5X3V8zKwM5jYJZoXB60zpYVFDph9M+Wz95Y9w5KMkyMQzjDW5NfBMBjFY+jwbwZZxjr4TGvKcotSiXAVjlUh5GAmkVzgKiyv8w0LBWVtAiMAzjE267wWMB/I+IbgKwDcCvAYCIWgB4QwhxsfK9JoALAPxOc/5zRNQdkQ5vgc5+JgCiI6joHFSYFJRVL8F4uJ07ClF1MAqz7x+MQydOBS0GIxlXb7sQYh8innna7TsBXKz6fgxATNwYIcS1bq7PeANVKKjYEdTQTvZCEcnmhv5t8OyUda7KkBXpmiNmh4ecRjWDFoHxAI4kwcRwUeemyEwnjOoZSYNRVl65L+hFpVkZaSgYOwJDOzV1XAZ73zFMYsCx+JgYWjWoiY3PVAyAq2S5DYt56/6LOmL3oRPo1944WZ0RrJ8YJjHgERRjyoAO2bihf67yLRwaqmOzOvjqrnMN56TObl3fX4EYhpEOKyjGEv3aRUYqYRlB6VGvRmUEiOFdmxsexyY+hkkM2MTHWCLqeh5W/bTuqWEAgDMemWp6rFvnhrDWAcMkG6ygGEu0bhjxkurVpmHAkuhTPTNdeow8M3gkxjDekpIKasa950nLCZQqnN60Dn54YDBaNfAucnGiUDMrHUBEKTIM4x0pqaDaGoTHYeITHUWFFTMX+Eu7t8AXy3a6vs51fXNxsrQcNw/gZHgM4yUpqaCY1CQjXY5PUFZGGu4Y3F5KWQzDGMNefEzqwdZdhkkIWEExDMMwoYQVFJO09GkbTo9DhmGswQqKSVo+uLkPp+dgmASGFRSTVPRp2xDjrukBAEhPo4rlBAseGhLqKBgMw8TCXnxMUvHxrX2DFoFhGEnwCIpJanKUtVuZklzMGYbxDx5BMUnN69flYeGW4ipzUexlzjCJAXcrmaSmYa0sDOvSLGgxGIZxACsohmEYJpSwgmIYhmFCiSsFRUS/JqLVRFRORHlxjhtGROuJKJ+Ixqi2NySib4loo/K/gRt5GCYegjM5MUxC4XYEtQrA5QBmGx1AROkAXgYwHMCZAK4kojOV3WMATBdCdAAwXfnOMJ5iFvWcYZhw4EpBCSHWCiHWmxzWC0C+EGKzEKIEwMcARir7RgJ4R/n8DoBL3cjDMAzDJA9+zEG1BLBd9b1Q2QYATYUQuwBA+d/EqBAiupWIFhPR4qKiIs+EZZKXOwe3xxnN6mBoJ8PHjGGYEGG6DoqIvgOg56f7ZyHElxauoWdPsT0ZIIQYD2A8AOTl5fFkAmObttm1MfWegUGLwTCMRUwVlBBiqMtrFAJorfreCkA0reluImouhNhFRM0B7HF5LYZhGCZJ8MPEtwhAByJqQ0RZAEYBmKjsmwhgtPJ5NAArIzKGYRgmBXDrZn4ZERUC6AtgEhFNU7a3IKLJACCEKAVwJ4BpANYC+J8QYrVSxFgAFxDRRgAXKN8ZhmEYBiQSMAdBXl6eWLx4cdBiMAzDMBIgoiVCiJi1tBxJgmEYhgklrKAYhmGYUMIKimEYhgklrKAYhmGYUMIKimEYhgklCenFR0RFALa6LKYxgL0SxEkWuD4q4bqoCtdHVbg+qiKjPk4TQmRrNyakgpIBES3Wc2tMVbg+KuG6qArXR1W4PqriZX2wiY9hGIYJJaygGIZhmFCSygpqfNAChAyuj0q4LqrC9VEVro+qeFYfKTsHxTAMw4SbVB5BMQzDMCGGFRTDMAwTSpJeQRHRMCJaT0T5RDRGZz8R0b+V/SuIqEcQcvqBhbq4WqmDFUQ0l4i6BSGnX5jVh+q4nkRURkRX+Cmf31ipDyIaRETLiGg1Ec3yW0Y/sfC+1COir4houVIfNwQhpx8Q0VtEtIeIVhns96YdFUIk7R+AdACbALQFkAVgOYAzNcdcDGAKIqnp+wBYELTcAdZFPwANlM/Dk7UurNaH6rgZACYDuCJouQN+PuoDWAMgR/neJGi5A66PhwD8VfmcDaAYQFbQsntUHwMB9ACwymC/J+1oso+gegHIF0JsFkKUAPgYwEjNMSMBvCsizAdQX0k/n2yY1oUQYq4QYr/ydT6AVj7L6CdWng0AuAvAZwD2+ClcAFipj6sAfC6E2AYAQohkrhMr9SEA1CEiAlAbEQVV6q+Y/iCEmI3I7zPCk3Y02RVUSwDbVd8LlW12j0kG7P7OmxDpESUrpvVBRC0BXAZgnI9yBYWV5+N0AA2I6HsiWkJE1/kmnf9YqY+XAHQCsBPASgB3CyHK/REvdHjSjma4LSDkkM42rV+9lWOSAcu/k4gGI6KgzvVUomCxUh//BPAnIURZpJOc1FipjwwA5wAYAqAGgHlENF8IscFr4QLASn1cBGAZgPMBtAPwLRH9IIQ45LFsYcSTdjTZFVQhgNaq760Q6e3YPSYZsPQ7iegsAG8AGC6E2OeTbEFgpT7yAHysKKfGAC4molIhxBe+SOgvVt+VvUKIowCOEtFsAN0AJKOCslIfNwAYKyKTMPlEtAXAGQAW+iNiqPCkHU12E98iAB2IqA0RZQEYBWCi5piJAK5TvFD6ADgohNjlt6A+YFoXRJQD4HMA1yZpr1iNaX0IIdoIIXKFELkAPgVwe5IqJ8Dau/IlgAFElEFENQH0BrDWZzn9wkp9bENkNAkiagqgI4DNvkoZHjxpR5N6BCWEKCWiOwFMQ8Qr5y0hxGoiuk3ZPw4R76yLAeQDOIZIryjpsFgXjwJoBOAVZdRQKpI0arPF+kgZrNSHEGItEU0FsAJAOYA3hBC6bseJjsXn4ykAbxPRSkRMXH8SQiRlGg4i+gjAIACNiagQwGMAMgFv21EOdcQwDMOEkmQ38TEMwzAJCisohmEYJpSwgmIYhmFCCSsohmEYJpSwgmIYhmFCCSsohvEJImqkRAJfRkQ/E9EO5fMRInolaPkYJmywmznDBAARPQ7giBDib0HLwjBhhUdQDBMwSo6lr5XPjxPRO0T0DREVENHlRPQcEa0koqlElKkcdw4RzVKCtk5L0gj8TIrDCophwkc7ACMQSWHwPoCZQoiuAI4DGKEoqRcRyU91DoC3ADwTlLAM4xVJHeqIYRKUKUKIU0oInXQAU5XtKwHkIhLzrQsi0bOhHJOM8SOZFIcVFMOEj5MAIIQoJ6JTonKiuByRd5YArBZC9A1KQIbxAzbxMUzisR5ANhH1BQAiyiSizgHLxDDSYQXFMAmGkoL8CgB/JaLliCTN6xeoUAzjAexmzjAMw4QSHkExDMMwoYQVFMMwDBNKWEExDMMwoYQVFMMwDBNKWEExDMMwoYQVFMMwDBNKWEExDMMwoeT/AVjBuquSuZ9sAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "wave.plot()\n", "decorate(xlabel='Time')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As expected, it is more random-walk-like than white noise, but more random looking than red noise.\n", "\n", "Here's what it sounds like:" ] }, { "cell_type": "code", "execution_count": 15, "metadata": { "jupyter": { "outputs_hidden": false } }, "outputs": [ { "data": { "text/html": [ "\n", " \n", " " ], "text/plain": [ "" ] }, "execution_count": 15, "metadata": {}, "output_type": "execute_result" } ], "source": [ "wave.make_audio()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And here's the power spectrum:" ] }, { "cell_type": "code", "execution_count": 16, "metadata": { "jupyter": { "outputs_hidden": false } }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAABI6ElEQVR4nO2dd3gc1fX3v2ertOrVklVsy3KTG8bCjWa6KcbEVBMggGMCCfxSCAkpQEIgkMCbQjDFwcak0CGhGUwoxrg33AuWZBvJtqze60r3/WN3dmd3Z7aXkXQ+z8OD9s7MnavRer97yj2HhBBgGIZhGK2hi/UCGIZhGEYJFiiGYRhGk7BAMQzDMJqEBYphGIbRJCxQDMMwjCZhgWIYhmE0iSHWC/BGZmamGDlyZKyXwTAMw0SQ7du31wkhstzHNS1QI0eOxLZt22K9DIZhGCaCENExpfGoCRQR6QD8DkAygG1CiJeidW+GYRhm4BFSDIqIVhBRDRHtdRufR0SHiKiMiO63Dy8AkAegF0BVKPdlGIZhBj+hJkmsBDBPPkBEegBLAVwKoATAIiIqATAOwEYhxE8A3BXifRmGYZhBTkgCJYRYC6DBbXgGgDIhRIUQogfAq7BZT1UAGu3n9IVyX4ZhGGbwE4k08zwAlbLXVfaxtwFcQkR/A7BW7WIiuoOIthHRttra2ggsj2EYhhkIRCJJghTGhBCiA8BiXxcLIZYR0UkA800m0/Swr45hGIYZEETCgqoCUCB7nQ/gRATuwzAMwwxiIiFQWwGMIaJRRGQCcAOAdwOZQAjxnhDijpSUlJAWUt3cha5eDncxDMMMREJy8RHRKwDmAsgkoioADwkhlhPR3QBWA9ADWCGE2BfgvPMBzC8uLg5oPUIIHK5pw0d7q/HR3mrsP9mCwnQL/rV4JgozLAHNxTAMw8QW0nJH3dLSUuGrkoQQArurmvHRvmqs3luNirp2EAHTC9Nw1phMvLj+KMwGHf65eCbG5SRFaeUMwzCMvxDRdiFEqce4FgVKZkEtOXz4sMfxvn6BbUcb8OHeany8rxonmrug1xHmjM7AJRNzcHHJMGQnxwEAvj7VipuXb0ZXbz9W3HoGpo9Ii/JvwzAMw3hjQAmUhNyC6rH2Y0N5HVbvq8bH+06hvr0HJoMO54zJwrxJObhwQjZSLSbFeSobOnDT8s2oaenG8zdPxzljPWoSMgzDMDFiQAmUZEEVjS5esvS/a/HR3mp8erAGrV1WJJj0OH/CMMybmIO547KQYPYvjFbT2oXvrNiKsppW/OX6abh8Sm5kfwmGYRjGLwaUQEnEDR8jcm75C9IsRlxUMgzzJuVgzuhMxBn1Qc3X3NmLxSu3Yvs3jfj9tyZj0YzCMK+YYRiGCZQBKVDDiyeKN1d/gRkj02HQhycjvrOnD3f9ezvWHKrFz+eNx11zR4dlXoZhGCY41ARKkx11iWg+ES1LoB7MGZ0ZNnECgHiTHstuLsWVU4fjDx8dxGOrDkDLIs0wDDNU0aRAhWujrhomgw5/uf403DxrBJ5fW4Gfv7Ub1r7+iNyLYRiGCQ5Nd9SNJDod4eEFE5FqMeJvn5WhpdOKvy46DWZDcPEthmEYJrxo0oKKFkSEey8ehweuKMFH+6px+8qtaOu2xnpZDMMwDDQqUFIMqrm5OSr3W3zWKPy/a6diU0UDvv3CZjS290TlvgzDMIw6mhSoSMeglLh6ej6eu2k6DpxswXXPb0R1c1fU7g0A7d1WPPL+fqw7XBfV+zIMw2gVTQpUrLioZBheum0GTjZ34epnN+BIXXtU7vtNfQcWPrMBL6w7gle2fhOVezIMw2gdFig3Zo/OwCtLZqGztw/XPrcB+05E1s24oawOVy5dh+qWLozIsKCyoSOi92MYhhkosEApMDk/Ba9/bzaMeh1ueH4TthxpCPs9hBBYuf4Ibl6xBVmJZrzzgzNxZnEmCxTDMIwdTQpUtJMklCjOTsSbd81BVrIZNy/fjM8Ongrb3N3WPtz/1h785r39OG9cNt7+/hyMzExAQZoFjR29nEnIMAwDjQpULJIklMhLjccb35uNMcMSccc/tuO/Xx0Pec6a1i7c+PfNeG1bJe45vxjLbp6OpDgjAKAgPR4A2IpiGIaBRgVKS2QkmvHKklmYPiINP3ptJ17acDTouXZXNWHB0+ux/0QLlt54Ou69eBx0OnIcL0izdf1lgWIYhmGB8oukOCNeun0GLpwwDA+9uw9//eRwwPX73tl5HNc+txE6Irx512zFdh8F6XaBauwMy7ojzYayOtS2dsd6GQzDDFJYoPwkzqjHczedjoWn5+HPn3yN3763H/39vkWqr1/gsQ8P4Iev7sRpBal49+4zMXG4susyzWJEgkkfNQvqeFNn0IVytx5twI0vbMayteVhXhXDMIyNqAkUEc0loi+J6Dkimhut+4YTg16HJ6+ZitvOHImVG47ip2/sQq+XIrPNnb1Y/NJWPP9FBW6aVYh/fXcmMhLNqucTEQrSLahqjLxA7T/RgrP+8BnWlQW+MbjH2o9fvr0HAFBW0xbupTEMwwAIUaCIaAUR1RDRXrfxeUR0iIjKiOh++7AA0AYgDkBVKPeNJTod4cErSvCTi8bi7a+O465/bUdXb5/HeeW1bfjW0vVYd7gOj35rEh65ajKMfrQNKUi3oLIh8i6+Lw/XQghg7/GWgK9dtrYch2vakJcaj/La6GxmZhhm6BGqBbUSwDz5ABHpASwFcCmAEgCLiKgEwJdCiEsB/BzAb0O8b0whIvzfBWPw8IKJ+PRgDW5ZsQUtXb2O458frMFVT69Hc2cvXl4yC9+eOcLvuQvSLKhs7Ih4j6rN9r1dFbWBWUBH6trx1GdluHxyLq6eno+qxg50Wz0FmmEYJlRCEighxFoA7rtYZwAoE0JUCCF6ALwKYIEQQvKFNQJQ93MNIG6ZPRJ/uf407DjWiEXLNqGurRvPfVGO21/aioJ0C965+0zMGJUe0JwF6fHo6OlDQwQL1vb1C2yVBCqAck5CCPzqP3tg1uvw4PwSjM5KQL8AjtUH75Isr23Dt55Zj5PNAyMxhGGY6BGJGFQegErZ6yoAeUS0kIieB/BPAE+rXUxEdxDRNiLaVltbG4HlhZcFp+Xh77eUoqymDXOfWIPHPzyIyybn4s27ZiPfnjYeCI5U8whm8u0/0YLWbivSLMaA6g2+veM4NpTX42eXjsew5DgUZSYCCNwKk7Opoh5ffdOEF9cfDXoOhmEGJ5EQKFIYE0KIt4UQ3xNCXC+EWKN2sRBiGWwuwB0mkykCyws/543Pxr++OxNpCUbcd8k4PL1oGiym4HpBOlLNI5jJt/lIPQDg6tPz0dDeg6YO39ZaQ3sPHvlgP04vTMW3ZxQCAIqyEgAgpDiUFG97ZfM3jgoa1r5+3PfGrphWdt97vBkr1x+J2f0ZhomMQFUBKJC9zgdwIpAJtFJJIhDOGJmOL392Pn5wXjGIlDTaP/LTbNUkvomgQG2qaMDIDAtmj84A4J/APPXpYbR2WfH7hZMdm4sTzAbkJMehPAQLqrKxA/FGPVq7rXh9q83wXrH+CN7YXoVPDoSvvFSgvL6tEr/74EDEY4EMw6gTCYHaCmAMEY0iIhOAGwC8G8gEWqjFFysSzAZkJJgilmre3y+w9WgDZo7KQFGW/y66rUcbMKc4E+Nzkl3Gi7ISUBGSBdWB0pFpKB2RhhXrj+BYfTv+/L/DAIBGPyy7SNHS2Yu+foGWLq6LyDCxItQ081cAbAQwjoiqiGixEMIK4G4AqwEcAPC6EGJfIPMORAsqnORHMNX8YHUrmjt7MbMoHQVp8TDoyGeiRH+/QEVtO0bbXXpybALV5tPS6O3rx7F6z/tUNnSgIN2C7549ClWNnVi0bBN0ZKuD2NTRqzBTdGjutN3bH/cnwzCRIdQsvkVCiFwhhFEIkS+EWG4fXyWEGCuEGC2EeDTQeYeyBQUABWnxqIyQBbWpwhZ/mlmUAYNeh8IMi08L6mRLFzp7+zDabnHJGZ2ViJYuK+ravH+Q//o/e3Hxn9eiXVapva3bisaOXhSkWXBRSQ4K0uNxorkL9148DqOzE2MqDpLl1BhDkWSYoY4mSx0NdQuqIN2CE02d6POjlFKgbD5Sj4L0eOSl2mJdRZmJPjP5JAFTEih/3ISbKurx2rZKdFv7cbC61TEuJYIUpMdDryP89sqJWDSjAN+ZMxKp8UY0dcZOHFrs946lm5FhhjqaFCi2oCzo7ROobukK67z9/QJbjtjiTxKjsxJwtL7DqxiW28sZjc5WcPFl2sbU3ITd1j786j97kJloy8g8cNJZucIhUPbU+vPHD8NjC6dAryOkWYxo9GMvWFdvHxYsXY/tx8LbVJJdfAwTezQpUGxBRaYv1OGaNjR29GKmbPNwUVYCeqz9OO5l31V5bTuS4gzIUqgjmJcaD7NBp2pBPf9FBcpr2/HEtVORHGfAfrlA2e8ppdbLSbWY0NJlhdVLrUPAlu24q7IJ2442ej0vUKTKII3tnlZcRW0bnlh9MKgMv+aOXlzy57X4+lSr75MZZoijSYFiCyoye6Gk/U+zipwWlOSiK69Td9GV17ahKCtRMX1epyOMykxQTFWvqG3D05+X4YopuThvXDZKhid7WFAJJj3SLEaPa1PtY76y6KR2Hw1htHS6rX3o6rUJo5IF9dG+aiz9vBy1bYG3Gjne1IlDp1qx9/jQfG8zTCBoUqCGugU1PDUeROGvJrGpoh7DU+Ice60AYJTkovOSKl5e26aYwScxOitR0YJ66N19MBt0ePCKEgDAhNxkHDzZ6nAnVjXaMviUhC/NYnMJ+ooBOQTKR5JGILR0OkVRSfikRA8l68oX0u/e3sP1CxnGF5oUqKGOyaDD8JR4VIXRghLCHn8qynARhIwEE5LjDDiiYkG1dVtxqqVbMUFCoigrAZWNneixOt1xR+va8eXhOnx/bjGyk+MA2ASqs7fPkW7+jT3FXAnJgpIsmIb2Htz98g40u2XVSQIVbDLD//afwt0v73AZkxf+rWvtwe6qJpfjbXarrr49cAvK2m97Rh3dvL+KYXyhSYEa6i4+wFZRIpyp5uW1bahr68GsItfitUSEoqxEVQvKWwafRFFWAvr6Bb5pcM7x4d5qAMD8qc7OwSW5tk2++0+2QAiByoZOhzvTnVS7BSXthdpypB7v7z6JHZWusSbJzVYfZHHdTw+cwvu7T7qIa4sse/CjfdVYsHQ9qpudCStt3TbrJxgLqt8et2pngWIYn2hSoIa6iw8If1+oTRW2LDd5Bp+Et2oQUhmjYoUMPsf19qKxZTXOOT7aexJT8lNcCuaOGZYIg45w4GQL6tt70Nnb50gIcUeKS0n7kGrsllK9myvPYUEFKVAn7MLTLBMl6WdpDULA5cuCJC7BxL2sfeziYxh/0aRAMbZEiVOtXWHrtbT5SAOGJZsxIsPTYhmdlYjqli7Fb/XlNe3Q6wiF6V4EKktKNbeJ2fGmTuyqasa8STku55kNehRnJ+LAyVaPFHN3nBaUTQRqWiSBcnWrSQLlzYLacqQBSz8vAwC8vaMK972xy3HsZJPtS0Bzp/N6KTGjUOZ+PNHk/LLQ3mMXqCDiXlIMqqOHLSiG8YUmBYpdfLZUcyHgNf3bX4QQ2FRRj5mjMhQTEqS9TEobdstr21CYboHJoP5WSYozIjvJ7LDCPrK79y6dlOtx7oTcZOw/0eI1xRwAkswG6Mjp4jtl3xPmLkSSQLV2WdGrkpL+zJoy/PVTW32/Lw/XYdWekwBsz0USHnlZJcnF19PnTCN3dfFJVSaCECi7i09yEzIMo44mBYpdfM4P7nBUNT9S147a1m7MLFJunjjK0TbDM1HCVwafhFSTD7C598bnJDkyBOVMyE1CdUsX9tgTD+QZhXJ0OkKqxeQQAcnFV+duQbV1w2Cvrq7k5uu29mFTRT16rP3otvahtcuK9p4+9Pb1o8X+M2BzJd68fDM+3HPS4eIbluzc93VSLlB2CyvQppLPrinHzcu3AOAkCYbxB00KFBPexoVSe3f5/ic5IzMSQORpQfX1Cxyt6/CaICFRlJWI8tp21LR0YduxRg/3nkRJru1Lx+p9p5CRYEKCWb1vVqrFWe5IKQbV29ePhvYex/qUYkLbjzY69jS1dlnR1m2br7mz16WLb0VtG748XIe7/r0Dp1q6YDLo8OS1U/H0jdMwJjvR5VxHDCpAgVq2ttw5h93FV17bhre2VwU0D8MMFVigNEp2khkmgy4sqeabK+qRmWh2uPLciTPqkZca75EoUdXYgZ6+fr8EanRWIpo7e/Hq1koIoezeA2wWFGCzDPNV3HsSqfFGRwyqtlVy8TktKEmsxubY5lQSjLWypoetXVa02q2fpo5enGxyWkWHa5zW4z82HoPZoENmohlXTBmOnJQ4RRdfoAKlk7lX2+0uvnl/WYt7ZTGxQGhs78H9b+1GVy+7C5nBCQuURtHpCPmpoaeaCyGw+UgDZhale22kWJSV6EhykJBcfkV+uvgAW7PBoswEjB2mLGoZiWaH66xAxb0nkWYxobG9F9a+fkfsqa7VKQpS/Gm8F4H68nCtwwXY2tXrEJfmzh4clyU+SAI1rTAVAKDXOZ9Vbkqcw8UnhJC5BQMTKPnzlyyoXnucqz+IwsBPfHwIr26txFs72AJjBicsUBomHH2hKhs6cbK5C7NGKcefJIoyE3Cktt2lvly5PW3cLwvKnmre1NGLeZNyvIrhBPt+qEIfFlSKxYjmzl7UtfVACCAl3oj69m7HGmvbbKIxdphNoNxjUI3tPdh3osXRObi1y+qIH0kuPr2OkBxnQJm9Nt5frj8Nb945Gy/eeoZjntyUeNS2daO3rx/d1n709QvoyJawoVaPr6PHin9tOuZyXKZ56HBLkujxUXNQCe72ywx2NClQnMVnIxx9oeT9n7xRlJWA9p4+nGpxutDKa9uQnmBCWoLJ533y0uIdmX5q7j0JSaDUMvgk0uxJElIG34TcJPT2ObvcShbUmGybOLpn+EnV4KcVpAKwWVDuLr6c5DikJ5gcVlFmohmlI9MxrTDNMU9uShyEsGXySdfnpsSjx9qPDpX9TH/48CB+/d+9WHOo1jHm6uJzTZLotqoLVEePFb9fdUDVlUdQ/zIAACebO/HJ/lNezwGA/SdaPJJQGCaWaFKgOIvPRkG6BU0dvWjtCrxigcSmI/VITzA5PsTVkDbbyt18al10ldDrCKMyEpCfFo9Jeclezy3x04JKsxjR0dOHKnuiiCRs0l4oSaByUuKQHGfwsKCk1HEp1lXX1uOwVJo6enG8qRO5KXFIse+5SjDpFZM2Sobb7rvjm0aHsEgbjNXiUJJYyssmyY3K9h6riwXUoyBQzR29mPbwx/jBv3dg2doKrNxw1OW4vwbUwmc24Lv/2ObzvMue+hIX/ukLx2trXz9++OpXXHmdiRmaFCjGhrOqefBuvs0VDZg5ynv8CZBttpUlSthSzH279yQeml+CP14zxee9Lp44DI9cNcml7YcSknBIH5AOgWqXEie6kRxnQJxRj4xEs4cFJW2+lZ6jPBPP5uLrQm5qPFLjbRUjspI824kAwMThKUiOM+D93Sfx7BpbJp4krmoCJcW9+l1cfM7n0i9crSYlF9/2bxrQ2NGLz+1W2J6qZsVUeh+P2yVF3hfy/WAHTrbinZ0n8JPXd/p9PcOEExYoDePoCxWkm6+qsQPHmzp9CgEA5CTHId6odwhUY3sP6tt7/EqQkJhTnIk5ozN9nmc26HHTrBEw6L2//aRSQ4drWkHkTIZwWFBt3Q5RSbMYPZIWpDJJ0l6rE7KsvaaOHlQ3d2F4SpyjMK2aQOl1hNmjM/C//afw2rZKAE7Rk6e2y8sl6ewCZZVt9nUXkjaZm0/JgnJ33X2w5yTm/XWt43WgISh/Y1aSlSjAMS4mtrBAaZhQ+0Jtlurv+Yg/Ac6+TpKLT/p/IBZUuJFabhyqbkW6xYQce1X02janBSWJSnqC2aNOn2QNZCSakGDSu5QrqqhrR09fP4b7YUEBwBkjXUVeip9JFs3RunZMe/hj/PzN3bjuuY3Q29XI2q9sQQGuiRKKJa0ULCN5jFASEB8GlPN8P/Vm4kOr0dvX7zjfV4zLOb/Ax/uqvXZnZphAiKpAEVECEW0noiuied+BSqrFiCSzwRGDCZTdVU1IMOkxzp7l5gt50dhAMvgiRYpdOI7WdyDbnswAuMagspJsopWe4GlBNXX2wKTXId6oR1Kc0cXVtf+ErXGiPAal1DFY4qZZI/DUommO1+4xqFMtXegXwGvbKrHlaIMjhV2eDKFz+5xv7/FlQfmHLxefRL8XhZrx6Ccur9u7rY71Kc3/6pZvsOMb18ry7+8+iTv+uR0vrj/i34IYxgchCRQRrSCiGiLa6zY+j4gOEVEZEd0vO/RzAK+Hcs+hBBHZU82Ds6AOnWrF2Jwkh7vJF0WZCahq7EC3tQ/ldW0w6XWqpYiigZQ92NcvkJ1khkGvQ5rF6LCUalu7HaKSnmBGY3uvixuruaMXKRYjiAhJcQZHDCrNYnTEq/y1oOKMelw5dTh+/63JAGzCbdCRQ6DcY0gbym3Zk61dVnyw+yTOf3KNh8NM7hJUEihfBOri82bYSJU6JK5+dgNu/PtmAMpCef/be7DwmQ2Kc8j3l7nz+cEal1ggw3gjVAtqJYB58gEi0gNYCuBSACUAFhFRCRFdCGA/AN/5royDYFPNhRA4VN3qt/UE2Dbr9gvgm/oOlNe0Y2SmxWecKJLIW8Fn28XDlgzRjebOXrT39CE7WRIoI3r6+l3iOk0dvQ7xSYozODbFyluA5PoRg5Jz48xCVPz+MqRabOn3ktWmVqi2vduK+9/ajYq6dsceLAm5KCgJlL+uMn9dcN4sKHfK5VVF3Ewotd9VQn6burZul/T421ZuxYKn1/u9DmZoE9KnjxBiLYAGt+EZAMqEEBVCiB4ArwJYAOA8ALMA3AhgCREp3puI7iCibUS0rba2VumUIYXUFyrQTZm1rd1o7OjFuJxABEoqGtuOigAz+CJBvFEPk10gJSHKSDChrq0H24/Z3nZT8mxbEdITbMflTQSbOnsc4pMU5xQ7ySo0G3RITzAFJFCAMwEi3WJyWHM9VuW/T1u31bE/rN2txUZNi9Pl2K3woe/LqlK649ajDarC9rfPDuNPHx/yOqc/LFtb4fJ6y5EGCCEUZbL0kU+wxC3F3d1aYxg1IvH1OA9Apex1FYA8IcSvhBA/AvAygL8LIRT/9QkhlgkhSoUQpVlZWRFY3sCiIC0enb19qAuw99DBaltqdiACJVUfP1TdimMNHQFl8EUCInKIR7Y91pSZaEZ9Wzc2VzTAqCfHhtr0BNt5dbJafU0dvUiJt7kJk+Kc+5vyUm0CNTw1HkSE6YXpWHh6HqaP8J3tKCctwYj/HTiFXZVNqlZFa7cVZrtASUVrJeSxRSUx2ntCeaP6h3tO4oPdJ52Wil0ZRt7/Aa59bqOjtYg7Sz8vx1OflXn7lRRxFx55av3H+6px3fMb8a9Nx1Sv/1JWD5FhAiESAqX0RcrxlU4IsVII8b7XCbiShAMpWyxQN98hu0CNz/G+aVaO1Ndpzdc16OsXMbegAGcmn9PFZ0J9ew82HWnA1PxUxJv0AGzljgw6wr83feO4trmz1+EmlCwok0HnsMZyU2yil2Ix4k/XneZIyvAXvY4gBPDj13Z6dfGp9dKSxxbdBaqitg1LPy93vwQAcNe/d+AHL+9wvHb/B3fwZIvXdR+ta3dxhR6s9n6+O3JjXoo3ldV4tmrxRoVCaxeGcScSAlUFoED2Oh/AiQjcZ0jgEKgAEyUOVrciK8nsyHzzl1GZCdhZ2QQgthl8EimSBWVPMc9IMKOpoxd7jze79LfKT7Pge+cW4a0dVdhQbvvG3tTR67DAku0WVJLZgFS7VZWbEloCyF3nFtt+IHV3XFuXukDJe3119fbhy8O1Dleutw7Bvlh7uBZn//Ez1dJIc59cgxv/vsnx+sevea+m7p7Ft+2Y06svFdXt8+GC/sfGo7DKRPy7L/mubMEwkRCorQDGENEoIjIBuAHAu4FMwKWOnEjxkkBTzb8+FViChERRVqLjG3KsXXyAM1FCsqAyk5yZfTNHue7vuuf8MShMt+DX/9mLrt4+dPb2OVrHSy6+xDgDku2WUl5qXEhrO2tMJi6fkgsI9cSBtm4rzAa9x7heRy5W8bK1Fbh5+Rbc+/ouvLj+iMeeKSXUqph39fajsqETR+s9OyRL7K5yeid85cHIV1Je2+ZyrbTO5k6rQ8gkkZXHTR98Zx9e2eK0br0JWllNq4uYMUOXUNPMXwGwEcA4IqoiosVCCCuAuwGsBnAAwOtCiH0BzssuPjsWkwGZiaaALKi+fmETqADiTxJS7b1hyWaXxIJYIbn4pASGDHsyhF5HmD4izeXcOKMed547GhV17dhz3PbeSYl3dfElxRkcVlVuaugp9BajHp29fS7t4eW0qlhQCSa9S0xKavfx9lfH8dv39geUcQcoV4noVUnccMeXGMpLVy1xs3ykkk7v7TqBY/W29+gHe05i3l/WeqS1t7hlMbqzp6oZh6pbceGf1uKxDw/6tXbAZn3+061yPDM4UG9n6gdCiEUq46sArAph3vcAvFdaWrok2DkGE/lploBiUMfq29Ft7Q9KoCSrSSoeG2umj0hDZWMH4ow2KyQz0SZYk/NSFAu7Sg0Rt9i7CDuz+OwWlNmA4uxEjM5KQKmbwAWDxaRHh72FvBKtXb2Ke8nMRj3QZUWS2YBWhfbvgTQhFALYe9wzjtTTF/5Ghi1uhYvlGXn77EkddW09tsK8XrIQj9V34MM9J1EyPBl3/GM7XvhOKeY/vc5hMS9fdwTbjzXivz840+eanlx9CC+sO4KsRBPm+aikzwwsQhKoSEFE8wHMLy4ujvVSNEFBugW77HEhf3AmSAQhUHZhGp0de/ceAFxbWoBrS50hzQz7xlx5/EnOmGFuAuXI4rN98CWajchMNOPTe+eGZX1xJpsF1Sv7MM5MNDmyLlu6rIpiI6XPD0uJQ6tCgkGnShsPJf78ydeKBWGvfnYjvvzZearX9fcLvzZxy89wN1L+9L+vHT/LC80CQG+/dzfdo6sO4JyxWTh0qhXv7rKFqRtlc+z08z0vxetafVhozMBDk7X4OAblSkFaPE40dfq9cfNgta246pjswAUqPy0ec0Zn4IIJwwK+NhoUpltw25kjceOMQsXjiWYD8lLjsf2YrQyPuwWVHBfe72QWowE91n50ykRIckNKyOvnSUip52qZg11u1sdUe08rJbxVKz/7j5+rHtt0pF71mBx51XVv70B3K7LX7XdQ6hrsqPfnZ7mmg9Utqq68+97c7d8kzIBBkwLFMShXCtItsPYLv0vEfH2qFSPSLY4U7EAw6HV4ecksnDcuO+Bro4FeR3ho/kSMyFC38MblJDnSqFPi3Vx8YRaoeJPtn5C8bFGKxVV0qls8BUSKS5kNOsUYlbvVZfSzXFUgHK2zuY19hW6keF5TR49qexHAtTCu0mvl20gFb33/fhvK6jDvL1/iX5u/cRnn2NPgRZMCxRaUK4H2hTpUHVyCxGBhzDBn/MyZZi65+MItULb55AkA/txDLlBmBYHqlgnUS7fP8NvCCASrDxecnAMnWzwqQrjjnmnqHoNy1xEi59gfPlJOipDv1zpiz0qUCv065vW6KmYgo0mBYlwJpC9UV28fjta3Y1wAG3QHG1J6vV5HDrGQUsuTA9yM6wuLPXlDbkFNHJ6My6fk4vGFk1Wvk2JQJhWBkmf4nTs2MhVVpNqE/vR9uvSvX2Lr0Uaf58lxt6AOnQpsQzAATHpotcKowIp1RzDy/g/wk9d2BjwnM3DQpECxi8+V4anx0BFQ5Ueq+eFTbegXwSVIDBbG2gUqNd7oSJFOiTfiT9dNxcLT88J6L8mN2iITqHiTHktvPN1r80anBeWsNyhHcvFdOz0fAHx2KQ4GX0VfQ6XaLTa2ak+1y+vKhk68urUS/vD2jioXy+nh9/fbxr86HnBVd2bgoEmBYhefK0a9Drkp8aj0Y7OuVLZmKLv4irMTQeQZC1p4er6jpl+4UBIoSXBSE9StNUmgTCoxqM7ePhh0hCeunQrA/95QgRDpzbC/eTeg7Y+qfPF1LX7y+i782y32JOGvPlU3d+G+N3YpN4dkNIkmBYrxpCA93qU0jhpfn2qFyaDDiHSLz3MHK3FGPUakWxytNiKJ5OJrkguUXXCSvMSiJBFTT5Lod+z9ihRPfvy175NCoKZVPbswEL6zYovLa3eLyd8kiQff2Ys3tlfh84M1YVkXE3lYoAYIBWn+NS48WN2KMdmJMe3jpAXuOX8MvjNnZMTv482CIiLF+BLg24LqsvYhzugcj0SSBOC9pUdeiJU2Gt32RYULX21I2rqt+N37+z0yIaWKGSs3HEVTRw/2VDVzBqDG0eSnGMegPClIt6CmtdtnhYGhnsEncfX0fCw4LbzxJiUsdoGSJwTIBWe4yoe8rxhUW5dVsVKGL+SipoaUdAPYqpG3d3u+p+48dzSm5GvTxf72V8ddXsslpsfaj2c+L8PydUc8WoBIIr+pogEzHv0U859eh/EPfIRn19iqxjd39uL+t3ajo4c3/GoFTQoUx6A8kT5UvBWNbWzvQU1r95BOkIg2Sm44o0xwpNJL7rhm8XnOcbypEznJzniZv11z/TnvrGJnVuB5T67BkTrPorI68l2jTzPIFIoIeMYuOO5ZhPLfp8cef+u29uMPHx3E54dqcPWzG/Dq1kosXsmV1rWCJgWK8cSxF8pLqrmzSeHQTTGPNhaTp5Ujt6AmqPwtpM9Ks0EHo4KLr7Khw8XF5q9W+HOeP95fHVFkMjMigLz9xz63PVIuePl9bntxq6On1cYK/ypsMJGHBWqAIPWF8pZq/vUpu0AF0WaDCQ6LQrUOuUAtmlmIqQWpOL0wVfF6s0Gn+LlZ09qNXFk7kHAaM3o/JvvOnJEDxoKSl5K6aul6x88eG4OjtSAmbLBADRCyEs0wGXReU80PVrciJd6IYclm1XOY8GI26DzEQx5Tykw0450fnInibNfq8NKHp9mgU22tIY9f+evi8we9zvc/+6wk84D/QHffgByuvWQ91n785LWdqAqwy7U/rDtch01swTnQpEBxkoQnOh0hPy3eaybfoeoWjMtJisimTkYZIkK8WxxKKStPLatSr9OpbjQdHkTHX3/+8r5cfHNG2xpBDvS30R8/OoTNFfUor23DyPs/cLQD8QepMsiObxrR7tYOZV1ZLd7+6jge+O9eAM629+7UtXVj7/HAPsNuWr4ZNyzb5PvEIYImBYqTJJQp8NIXSgiBr0+1cYJEDHB38yll5Rncir3Kq3j7ZUEFIRYThyvHv3Q6wuWT1fsmSckZA8XF540nVh/CZwds+54qatU7DLsz9bcf408fH8LCZzbgnle+Uj3vX5uO4czHP8P/9p/yOHbZX7/EFX9bF/iiGQeaFChGmYL0eNWCsVWNnWjrtnKKeQxwr++nZEHp3QXKUcVbXaByg2hJL7eeX7z1DOVzQN73iNmnUIubDSS2HWsM2hJ86rMyAM5q7gBwy4oteOSDA47Xv7ZbUTsrbXUKy2vbMPmh1TjZ3OnSzJEJDhaoAURBmgXNnb0eXU0BZ5NCTpCIPrkprkJiVLCg5GPP3TTd5Ziai0+qwA54xk/mjvNdQFbN1WtLIfd5OW6aNcL3SQOAUGM68r/P2q9rHZbY54dqPc65YdkmtHZbMefxzxzH/r62IqAOyYwTFqgBhJTJpxSHOmTP4BvLFlTUyXWLFSlVj5AsqIXT8jBvUo6Li09Jn9xdgmPdkiyyEs1YNKMA3lATISL410l3ELj4AOCTA6GVNqpr820JSX9DqSmjXNQeXXUAt8jKNVn7+v0ut9Rj7Y94UV8tEzWBIqIJRPQcEb1JRHdF676DCW99oQ5VtyIvNd7lWzcTHdwtKMUkCbsgKAmDUrkdg971vJ9fOt7ltV5HeGzhFI/rpKuevnGaqsAQaFDEl6JNeW2b6jHJTVuv0tBxy5EG9PULtHdb8cyacty2cqtfIjXpodWY8egnwS14EBCSQBHRCiKqIaK9buPziOgQEZUR0f0AIIQ4IIS4E8B1AEpDue9QpdCbBcUljmKGZEFNzkvBohmFyEz0TPM3uKV2S5JEICh0QofR7XyjXoeSXGfSg7umlY5Iw9M3TnMo1PicJO8WlBd9CmdK+2Chq7cPF/y/L1SP91qFzyzBB97Zi4kPrcaf/mcr0lvrJUb13q4TAGwVLyJV03AgEKoFtRLAPPkAEekBLAVwKYASAIuIqMR+7EoA6wB8GuJ9hyQpFiOS4gwemXw91n6U17axQMUIyYJKiTfisYWTPRIiAKdFJAnLneeOxuisBFxYMkwxSULJ0pKf577H58275uCKKcOd1xN5saAGR4ZeNPn2C5u9Hl+x/gguf8p7xt6b26tcXqslxwDAPa985Vdx6IraNuyqbHIZa+7sDTi9XauEJFBCiLUAGtyGZwAoE0JUCCF6ALwKYIH9/HeFEHMAfDuU+w5llKqaV9S1wdovOEEiRmQl2Swmb7EKycUnufOKsxPx6b1zkZ5gUrRX3GNQANAnM7WUrC7A6eLT60g9e428u/hYuzzZfiywbsJKuFdhl/6EQghH7EpOt4+q7QBw/v/7Agtk1TMA4NsvbPKZ3t7ebUWrQrKV1ohEDCoPgLxNZhWAPCKaS0RPEdHzAFapXUxEdxDRNiLaVltbq3bakKUg3bNxoSODjy2omCBZUEVZCarnSFaVkq78+frT8O2ZhUhPMDnG3GNQgKtA+ULnRYR0BHgrJsH6FB2EAPadaMaoX6xC0S9VPxIVWbn+CO57Y5fisb3HvdQjBFDV2IHTHv4Yk3/zcUD3dGf/iRbc/fKOiDa+DLyev2+U3t9CCLEGwBpfFwshlhHRSQDzTSbTdF/nDzUK0iz44utaCCEcLpxD1a0w6AijsxJ9XM1EgoxEM97+/hxHq3klJItIya0zIiMBj35rMj6WbfZ0j1kBQJ/cxadmQdnfE3odqcegOElCE+w53oRf/meP47V7EsaFf1KOeb21vQq/eW+/4rHHVjn3aH2w+yQun+K6IXvLkQZc9/xGn2vr6u2DUa9TdFdL3PPKDpTXtuNHF45BcXZkvhxHwoKqAiDPf80HcCKQCbiShDoF6RZ09fajVuZOOlTdiqKsBMXsMSY6nF6YhkQv/ZukUkf+GkG+LCgpBvX4wsl4+/tzPM7V60g12cGgD9zFN9qLdcgEx9E6V1e9tySMTw+cwrrDdTjV0oV7FSwnyUX4/NoKx9gPXt7h+PnzQzV4YvVBx3YUidtXbvVwLx6rb8f4Bz7COX/83P9fJkJEwoLaCmAMEY0CcBzADQBuDGQCIpoPYH5xcXEEljewkfpCVTZ0IjvJ5lo6WN2K00ekxXJZjA/cY1C+UKo47vJBIm0MnVGoeL0tSUJ57ksm5njN4lOyxNMsJgDtuHRSDj7cW61+MeM3gbT1WPyS9x5V976xC3++/jTV47e9uBUAXNzIAPDZwRoU/2oVKh67HD3Wfry4/ohjH9/xpk7c+PdNyEuNxxPXTsXRunbodeTYjykRyabEIQkUEb0CYC6ATCKqAvCQEGI5Ed0NYDUAPYAVQoh9Ia+UAeDcC1XV2IHpI9LQ2tWL402dPjdtMrHFEYPy8o9ZrhlKrhWriwWlMof9Mh0pW0K/u2oSirMTcVShSSFga/O+5Owi1XlvmFHIAqVB/vPVcfzHrdOwEg0K+7T6BVDT2oUZj9qSqwtlArSh3Cai00ek4f63be7It+6ag+kj0qKykTvULL5FQohcIYRRCJEvhFhuH18lhBgrhBgthHg0iHnZxadCfprrXqivT9n81tykUNsYHS4+/75uKlU/d0kz92MeJTeePMtPiRmj0r1WmeDI1cDisQ8P+D4JwMubv3H8/I1CerskTgBw9bMbXI5F0IDSZqkjbrehTrxJj8xEs6OahJTBx1XMtY0/FpQcX2nmqhaU7Gclgbpmer7tPFaaIcHzX1Rg21H3nUCebCwPrF7h3uPNjpYkkUSTAsUWlHdsqea2bzmHqluQYNK7tAdntIcjBuXl+6ZcNJQsnCunOjfiqgldXlq843r3GW6cWYg4e++qQLP45LUDmYHFNc/5ztrbfMS3iMm54m/rvFbCCBeRSJJgIkxBmgVf2cv7H6xuxdicJL+KfzKxQ+dIM/fvfCUL6sH5E1GcnYgH3tmnKnMv3joDGyvqkWoxebgBzyrOdK4nSKXhMkiMO5FMktCkBcUuPu8UpMfjRFMXrH39OHSqlStIDACkj3VvsaMkWaFfJQtKryOP3lPuZCWZHZaWPIi966GLcZmsSaHaRl1/swy9seJWLrU5lPDmFQgVTQoUu/i8U5BmQV+/wK6qZjR19HIFiQGAZLF4+/wvynTuNVLqKQU4RSdQIUlxEzY1C8rXrP4YXuePH+bnqpjBQFdv5CpJaFKg2ILyjrQP4ZMDtsoDLFDaR/pg95bFVyTbf2RU2KgLAGOH2c45e0ym4nF/8dfF98Q1U/DWXbNDuhfDBIsmBYotKO9I+xQ+sZfGGc8p5ppHEgRvMahrpucBAKbkp+Dxqz17PQG2v/VXD1yE60pD2/emtBFYiWtLCzB9RLrjNUegGHfC4RZWg5MkBiC5KXHQ6wiHa9qQlWT22B3OaA9/LKji7CQcffxyn3OlheHvTZr8asoMRIbcPijGOwa9zlFBmxMkBgZaa58edLFYbf0ajAaI5FtCkwLFMSjfSCWPOP40MHBm8cV0GQ7UXHzBrO/9e84KcTXMQGbIWVAcg/KNVDSWBWpg4IxBaUSh7ARqSCntg5qUx/9OhzJDbh8U4xvJguISRwMDnR8xqGgi7YOSb94Fgk8zX3rj6aEvimHc4CSJAcpFE4ehrLaNM/gGCvYPdo3oE8wGPT75yTnISYnHpIdWhzzf5VNykZk4C+09Vq/nleQm42B1i98VNZTITjKjJgpldhh/GWIbdRnfjM9Jxl9vmMZNCgcIJbm2LxJ3njs6xitxUpydhASTPqBrvHkEZxZl+Nyku+qHZ2P5rWcEdE93MhLNIV3PDBw0aUFxw0JmsJFqMfmVQh5tYpFdeN647JCujzPyl7Khgib/0pwkwTCxwdemy3ALWqrFe21BRvtwkgTDMIMSTvJhvMECxTBMzFDr7MswAAsUwzAywuWtef+es/DENcr1BJnBxaDYqEtEVxHR34noHSK6OFr3ZRgm+kzKS8G1fhS05QaIjDdCEigiWkFENUS01218HhEdIqIyIrofAIQQ/xVCLAFwK4DrQ7kvwzCDA39zLqQ0fUZ7aDlJYiWAefIBItIDWArgUgAlABYRUYnslF/bjzMMozV8fNjIBeWl22fg7e/PCevtzyzOUByXdwDWymZnxkYk222EJFBCiLUAGtyGZwAoE0JUCCF6ALwKYAHZ+AOAD4UQO0K5L8MwsefcsVk4vTAtpDnc09ZnFykLlHtH4EAYbq/8zww8IhGDygNQKXtdZR+7B8CFAK4hojvVLiaiO4hoGxFtq62tjcDyGIYZaDy+MPiEi3PGZoVxJYw72cmR+wIQCYFS8ioLIcRTQojpQog7hRDPqV0shFgG4LcAdphM3IiPYSLNPxfPwH2XjIv1MrySkRj8Z8HDCyaFcSWMO0Z95BJdIiFQVQDk6Tv5AE4EMgFXkmCY6HH2mCyMyLBVxxcRTRr2xN+PtqAbLAJcr3IAE4m/3FYAY4hoFBGZANwA4N1AJuCGhQwTXWKV7m1w26irFm/XWENiJkqEmmb+CoCNAMYRURURLRZCWAHcDWA1gAMAXhdC7AtkXragGCY2RCoh6/XvzcY/bp/hMT5vUo5f18stKE7iGzqEVM1cCLFIZXwVgFXBzsvVzBkmukTaQpkxKl1x3KhX/46cHGdAS5etv1QoLj4A+OEFY/DXTw+HNAcTfTTpnGULimEYOToC/vfjc/DRj872OLbpFxcoVkV/9+4zHedLMTZmYKFJgeIYFMMMfpSyvwSASyYOsx/X4dY5IwHY9kuNGZaE8TnJeNKtxl9OShyyFJoYTslPDarjdLwxsCaOTOTQpECxBcUwg5sXbzsDn907VzFr8MH5EwEABj3hofklqPj9ZS7HxwyLbIuO88bzvqlA0HKpo4jAFhTDRBcpmU4t1BPuz6DzxmWjIF3Z7SZP7CMi6Ly05Pjwh54uv1C5/cxRYZ+TCQ5NChRbUAwTXS6YMAw3zSrEb66cGOulBMSECBSRLR2pnNDBRJ+QsvgYhhkcGPU6PHLV5KjfV3IPpVmMaOzoxbTCVMcYt+JgNGlBsYuPYYYW547NwlcPXISzx0Qm/jM5j70xAxFNChS7+BhGW0TDlklLiFztzTHDklCoEvNitIsmBYphmKGB2WBL6U6Mc0YbpISMcG8ejmRRUyYyaDIGxZUkGCY8bP7lBSFXYYgk8ybl4GfzxuGW2SM9jml31Uy00KQFxS4+hgkPw5LjkJXkuYk1FD6791xs+eUFfp07zseeJb2O8P25xUg0h/Zd+ZbZI3ye494ckdE+mhQohmG0S1FWot9N6l7/3my8f89ZAc0fTAvxmxUsMF/kePkdpAoWTGxhgWIYJmKkWIyYFGQGXaQtHm/r+s2VE3nDrgZggWIYRlOEo3ROXmp8yHMYDewSjDWaTJJgGEZbDIQeTNdOz8d/dx7HF/ed55IVKOEuNyO5wrnm0aQFxRt1GYYJlCeunYrDj16G4anxSI7zbL/hzs/mjY/CqphQ0KRAcRYfwzDhxj2kZTJ4//i7dFJuBFfD+IMmBYphGCbWWeGnFaQiN8W/bEUmMrBAMQzDBMnRxy+P9RIGNZwkwTCMpohkAzwAuHBCNi6YMMzxujDdgm8aOhTP5Ty+2BI1gSKiIgC/ApAihLgmWvdlGGZgIXXZDbeLT2rf8dNLxjlawe988CKYDXpMePAjlbUwsSQkFx8RrSCiGiLa6zY+j4gOEVEZEd0PAEKICiHE4lDuxzAMAwA7HrgI2359YVDXyi20VIsJ8SZ9QNf/8jLO/osWocagVgKYJx8gIj2ApQAuBVACYBERlYR4H4ZhGAfpCSZkJoa3xqC/LDm7KCb3HYqEJFBCiLUAGtyGZwAos1tMPQBeBbDA3zmJ6A4i2kZE22pra0NZHsMwAxCtd9TlorOuRDJmGIksvjwAlbLXVQDyiCiDiJ4DMI2IfqF2sRBimRCiVAhRmpUVme6aDMMMPbzpymkFqUHPe8/5/rcF4rT1wIiEQCm9DYQQol4IcacQYrQQ4jGvE3AlCYbRJIPBdlD6xv/S7TOCnu8nF43FQ/N9RzHG5yThZ/PGBX2foUgkBKoKQIHsdT6AExG4D8MwUSYaWW2R6qjrjZR436WR1CAiGPTOj9I7zuEYVbiIhEBtBTCGiEYRkQnADQDeDWQCLnXEMEMXg86mTKkWU9TvffaYzJDnUEuiiGTsalphasTmjiUh7YMiolcAzAWQSURVAB4SQiwnorsBrAagB7BCCLEvwHm55TvDDFEK0i14eMFEzJuYE9X7hqsqhFCxMweDezTahJrFt0gIkSuEMAoh8oUQy+3jq4QQY+3xpkeDmJctKIYZwtwye6TfXXsDRU1AvPF/F4wJ+b5qBpQ/8SstE8zz9BdN1uLjJAmGYcJNKC62G84o8H2Sz/t7jv3x6in41rS8kOcerGiyFp8Q4j0A75WWli6J9VoYhhlcBLtv58XbzkCOD6tuspc28or7uki7+720gCYFimNQDMOEm1Bl4Lxx2T7PmZLvKVCZiWbUtXXj/PHZEdvUGukCu7FCky4+jkExDBNu7rtkHCwmPUZlJkT1vjfNKsTWX12IH4YhjqXGINUnbQoUwzDaZCA7o84bn439D89DgjmyjiPJZZeZaMJXD1yE/zt/DLKSzNDpVJ6efTjJbMDnP50b0bUNNDQpUJwkwTCMFpgz2rYvKsEUnKilJZjUhUmBAE4dEmhSoNjFxzCMFvj9wkn4/KdzkWIJvtJEVIhhEEpeRSPcaFKgGIZhtIDZEFjMypdMXDB+mMeYP9nvXz9yqfcTYlhhPS81PmJzs0AxDMOEEW8bV1MsRnz843Mcr/2VFZNhaH5Ua/K35hgUwzADEfLySgmh8nO4mZSXHMHZI4cmBYpjUAzDDBWCcc6t/tE5vk+ScdVpA7NahSYFimEYZuDjn00UjOU0LifJbZLBuROKBYphGCaMBFu6KJQ0h8xEs9fjGYnRb10SDligGIZhNEAoNlBuqq1G4OVTchWPEwgHfzfP6xwxTARURZMCxUkSDKMtfnnZeBRlJqBk+MAMtmuJ/DTXtOxAdWFkhkX1mEFlp+85Y7MQZ9T7mDe6JaD8QZMCxUkSDKMtpo9Ix2c/nQtLkBUVGCehPsN3fnCWx9hZxZmwmPRYfNYoxWvSE3y7+H63YFJI64oE/G5jGIaJEcF41ZSqWkh1BkMh3uTdwooFmrSgGIZhtM5UhdYasWDexByYDdoTl3DAFhTDMEyAlP/+MlXrx2y0fe+fmp8atfUMVqImUESUAOAZAD0A1ggh/h2tezMMw4QTvZey48lxRrx795kozk4MaE4R5F6maGXffe+cImw71ojtxxqjc0OE6OIjohVEVENEe93G5xHRISIqI6L77cMLAbwphFgC4MpQ7sswDKNlpuSnek2GuPr0fMfPFGGFOXtMZljm+cVlE/DmnbPDMpe/hBqDWgnAJTJHRHoASwFcCqAEwCIiKgGQD6DSflpfiPdlGIZh/OCh+SV+neePThIR9j98Ce48d3SIq/KPkARKCLEWQIPb8AwAZUKICiFED4BXASwAUAWbSIV8X4ZhGMY/kuLC28vKYjLAqI+OXzESQpEHp6UE2IQpD8DbAK4momcBvKd2MRHdQUTbiGhbbW1tBJbHMAwzeLl1zkjHz6v+72wMS46L3WJCJBJJEkrSKoQQ7QBu83WxEGIZEZ0EMN9kMk0P++oYhmEGEe6uudNHpGHlhqMAMOArf0TCgqoCUCB7nQ/gRCATcCUJhmEGEmeOzgAAXDl1eNBzhKseedDZgGG6fziJhAW1FcAYIhoF4DiAGwDcGMgERDQfwPzi4uIILI9hGCa8FGUl4ujjlwd1rRaFQSuEmmb+CoCNAMYRURURLRZCWAHcDWA1gAMAXhdC7At9qQzDMIMPf+2dX18+ATfNKnS8XnCad2vtnLFZiuPeis1qjZAsKCHEIpXxVQBWhTDvewDeKy0tXRLsHAzDMAMJX5bUd88ucnl9cUkO3tmpHj1JjffM3tv8ywuQYDZg0kOrPY7pNNhvQ5OljtjFxzDMUEPAv6rjEheVDMOiGQX48UVjXefxYpIpZfRdPiUXRZkJmKKR2oJyNLkfiZMkGIYZzAiZY09utwSyZ8lk0OGxhVOQnaScRu6vQbT0xtNx78XjXCparL3vPL/XEUk0KVAMwzBDAaUyRzsfvMjl9S8uHR/QnCKEfMDX7piF3101CYUaiVNpUqC4oy7DMEMNyTWXajFhdlGGY/xb0/KCmi+YiNLMogzcPGuEz/MyE81BzB44mhQodvExDDOUefG2M4K+1p9tUGqt4f1dw01+iFg40KRAsQXFMMxQQ+7tizPqQ7ZSIlkl3Vu7kXCiSYFiC4phmKGGquUToBYEWUgi1NtGBE0KFMMwzFAh3IaOpE/eptXglidFWKAYhmEGIwNEhLyhSYHiGBTDMIOZC8YPAwBMHJ7sqOAwMjMhLHP7Uyz2tjNHheVekUaTlSS41BHDMIOZy6fk4oIJ8xBn1AMAVtxaiin5qS7njM9Jwrqybpj1+qDuQV5MqF9cOh43zxqheVefJgWKYRhmsCOJEwCcb7eo5Dxz0+nYd7wFKZbAOuL6kyNBRChI974ZN5JZgP6iSRcfwzDMUCc5zojZozN8n6iCBvQlZFigGIZhBhPh6nyoATQpUJwkwTAMExrhNKAWnxWbpApNChRv1GUYhgmOUIrFypEL3ANXlIRlzkDRpEAxDMMwocExKIZhGEZT5KTEAwjfvqpYwmnmDMMwg4hzx2bh5SUzMXNU8BmAWoEFimEYZpAxZ3RmrJcQFqLm4iOiIiJaTkRvRuueDMMwTHBoIYbll0AR0QoiqiGivW7j84joEBGVEdH93uYQQlQIIRaHsliGYRhm6OCvi28lgKcB/EMaICI9gKUALgJQBWArEb0LQA/gMbfrbxdC1IS8WoZhGGbI4JdACSHWEtFIt+EZAMqEEBUAQESvAlgghHgMwBXBLoiI7gBwBwAUFhYGOw3DMAwTQdITTGho74noPUKJQeUBqJS9rrKPKUJEGUT0HIBpRPQLtfOEEMuEEKVCiNKsrKwQlscwDMMEi7dq6ACw9mfnYccDF0V0DaFk8SmtXnULsxCiHsCdfk1MNB/A/OLi4iCXxjAMw0SSRLMBMEf2HqFYUFUACmSv8wGcCG05DMMwDGMjFIHaCmAMEY0iIhOAGwC8G45FcS0+hmEYbXBmcew2/PqbZv4KgI0AxhFRFREtFkJYAdwNYDWAAwBeF0LsC8eiuJo5wzCMNvCjg3zE8DeLb5HK+CoAq8K6InDLd4ZhmFgzYDbqRhu2oBiGYRhNChTHoBiGYRhNChRbUAzDMLFleKqtbcfsotglSZCIZQTMB6WlpWLbtm2xXgbDMMyQpLKhA3mp8dDpIhuQIqLtQohS93Fut8EwDMMoUpBuien92cXHMAzDaBJNChQnSTAMwzCaFCiGYRiGYYFiGIZhNIkmBYpjUAzDMIwmBYpjUAzDMIwmBYphGIZhWKAYhmEYTaLpShJEVAvgmP1lCgD3oJQ/Y5kA6iKyQGWU1hSp6/0519c5asf5eQd3Lj/v8F3Pzzu6z9uf8yP1vEcIIbI8rhJCDIj/ACwLZgzAtlivM1LX+3Our3PUjvPz5ufNz3toPe9Qnme4nrf7fwPJxfdeCGPRJNT7B3K9P+f6OkftOD/v4M7l5x2+6/l5R/d5+3N+pJ+3C5p28YUDItomFIoQMpGBn3d04ecdXfh5R5eBZEEFy7JYL2CIwc87uvDzji78vKPIoLegGIZhmIHJULCgGIZhmAEICxTDMAyjSVigGIZhGE0y5ASKiBKI6CUi+jsRfTvW6xnsEFERES0nojdjvZahABFdZX9vv0NEF8d6PYMdIppARM8R0ZtEdFes1zPYGBQCRUQriKiGiPa6jc8jokNEVEZE99uHFwJ4UwixBMCVUV/sICCQ5y2EqBBCLI7NSgcHAT7v/9rf27cCuD4Gyx3wBPi8Dwgh7gRwHQBOPw8zg0KgAKwEME8+QER6AEsBXAqgBMAiIioBkA+g0n5aXxTXOJhYCf+fNxM6KxH48/61/TgTOCsRwPMmoisBrAPwaXSXOfgZFAIlhFgLoMFteAaAMvs3+B4ArwJYAKAKNpECBsnvH20CfN5MiATyvMnGHwB8KITYEe21DgYCfX8LId4VQswBwCGDMDOYP6Dz4LSUAJsw5QF4G8DVRPQsYl/GZDCh+LyJKIOIngMwjYh+EZulDUrU3t/3ALgQwDVEdGcsFjZIUXt/zyWip4joeQCrYrO0wYsh1guIIKQwJoQQ7QBui/ZihgBqz7seAH9Qhh+15/0UgKeivZghgNrzXgNgTXSXMnQYzBZUFYAC2et8ACditJahAD/v6MLPO7rw844Bg1mgtgIYQ0SjiMgE4AYA78Z4TYMZft7RhZ93dOHnHQMGhUAR0SsANgIYR0RVRLRYCGEFcDeA1QAOAHhdCLEvluscLPDzji78vKMLP2/twMViGYZhGE0yKCwohmEYZvDBAsUwDMNoEhYohmEYRpOwQDEMwzCahAWKYRiG0SQsUAzDMIwmYYFihjRE1EdEO2X/jYz1msIFEU0johfsP99KRE+7HV9DRKotIojoVSIaE+l1Mowag7kWH8P4Q6cQ4jSlA0REsO0V7I/uksLGLwE8EsL1zwL4GYAl4VkOwwQGW1AMI4OIRhLRASJ6BsAOAAVEdB8RbSWi3UT0W9m5v7I3sPuEiF4hop/axx2WCRFlEtFR+896InpCNtf37ONz7de8SUQHiejfdnEEEZ1BRBuIaBcRbSGiJCL6kohOk61jPRFNcfs9kgBMEULs8uN3vlJmQR4ioiP2Q18CuJCI+IssExP4jccMdeKJaKf95yMAfgxgHIDbhBDfJ1vb9DGw9QMiAO8S0TkA2mGrxzYNtn9HOwBs93GvxQCahRBnEJEZwHoi+th+bBqAibAVIF0P4Ewi2gLgNQDXCyG2ElEygE4AL8DWMfdHRDQWgFkIsdvtXqUA9rqNXU9EZ8leFwO2fkaw15UjotcBfGEf7yeiMgBT/fjdGCbssEAxQx0XF589BnVMCLHJPnSx/b+v7K8TYROsJAD/EUJ02K/zp3DoxQCmENE19tcp9rl6AGwRQlTZ59oJYCSAZgAnhRBbAUAI0WI//gaAB4joPgC3w9YB1p1cALVuY68JIe6W/a5r5AeJ6GewPQ95J94aAMPBAsXEABYohvGkXfYzAXhMCPG8/AQi+hEAtUKWVjjd53Fuc90jhFjtNtdcAN2yoT7Y/m2S0j2EEB1E9D/YOrpeB5u15E6n2729QkQXALgWwDluh+LsczFM1OEYFMN4ZzWA24koEQCIKI+IsgGsBfAtIoq3x3vmy645CmC6/edr3Oa6i4iM9rnGElGCl3sfBDCciM6wn58kiwe9AFtjwq1CCPf25ICt4naxP78gEY0A8AyA64QQ7mI0FgBX7WZiAltQDOMFIcTHRDQBwEZ73kIbgJuEEDuI6DUAOwEcgy2hQOJJAK8T0c0APpONvwCb626HPQmiFsBVXu7dQ0TXA/gbEcXDZslcCKBNCLGdiFoAvKhy7UEiSiGiJCFEq49f81YAGQD+Y/8dTwghLiOiYbC5/E76uJ5hIgK322CYMEBEv4FNOJ6M0v2Gw9ZqfLxaGjwR/RhAqxDihSDv8WMALUKI5UEvlGFCgF18DDPAIKJbAGwG8Csfe7SehWtsK1CaALwUwvUMExJsQTEMwzCahC0ohmEYRpOwQDEMwzCahAWKYRiG0SQsUAzDMIwmYYFiGIZhNAkLFMMwDKNJ/j8hWoZm2onEeQAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "spectrum = wave.make_spectrum()\n", "spectrum.hs[0] = 0\n", "spectrum.plot_power()\n", "decorate(xlabel='Frequency (Hz)',\n", " xscale='log', \n", " yscale='log')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The estimated slope is close to -1." ] }, { "cell_type": "code", "execution_count": 17, "metadata": { "jupyter": { "outputs_hidden": false } }, "outputs": [ { "data": { "text/plain": [ "-0.9993534088150459" ] }, "execution_count": 17, "metadata": {}, "output_type": "execute_result" } ], "source": [ "spectrum.estimate_slope().slope" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can get a better sense of the average power spectrum by generating a longer sample:" ] }, { "cell_type": "code", "execution_count": 18, "metadata": { "jupyter": { "outputs_hidden": false } }, "outputs": [ { "data": { "text/plain": [ "6553600" ] }, "execution_count": 18, "metadata": {}, "output_type": "execute_result" } ], "source": [ "seg_length = 64 * 1024\n", "iters = 100\n", "wave = Wave(voss(seg_length * iters))\n", "len(wave)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And using Barlett's method to compute the average." ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [], "source": [ "from thinkdsp import Spectrum\n", "\n", "def bartlett_method(wave, seg_length=512, win_flag=True):\n", " \"\"\"Estimates the power spectrum of a noise wave.\n", " \n", " wave: Wave\n", " seg_length: segment length\n", " \"\"\"\n", " # make a spectrogram and extract the spectrums\n", " spectro = wave.make_spectrogram(seg_length, win_flag)\n", " spectrums = spectro.spec_map.values()\n", " \n", " # extract the power array from each spectrum\n", " psds = [spectrum.power for spectrum in spectrums]\n", " \n", " # compute the root mean power (which is like an amplitude)\n", " hs = np.sqrt(sum(psds) / len(psds))\n", " fs = next(iter(spectrums)).fs\n", " \n", " # make a Spectrum with the mean amplitudes\n", " spectrum = Spectrum(hs, fs, wave.framerate)\n", " return spectrum" ] }, { "cell_type": "code", "execution_count": 20, "metadata": { "jupyter": { "outputs_hidden": false } }, "outputs": [ { "data": { "text/plain": [ "32769" ] }, "execution_count": 20, "metadata": {}, "output_type": "execute_result" } ], "source": [ "spectrum = bartlett_method(wave, seg_length=seg_length, win_flag=False)\n", "spectrum.hs[0] = 0\n", "len(spectrum)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "It's pretty close to a straight line, with some curvature at the highest frequencies." ] }, { "cell_type": "code", "execution_count": 21, "metadata": { "jupyter": { "outputs_hidden": false } }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAAAsBklEQVR4nO3dd3zV1f3H8dcnmxlG2AEChL0FwhRQqYKKWLUKauug4qKttLWlaltbF44OUX9aikq1ilonCogDlC1LpoyEJQGEhBEgAbLO74+EGCGBe3Pvzb1J3s/Hg8fj3vP9nu/95MuFT874nmPOOUREREJNWLADEBERKYkSlIiIhCQlKBERCUlKUCIiEpKUoEREJCRFBDsAgLi4OJeQkBDsMEREJAhWrlyZ7pxrcHp5SCSohIQEVqxYEewwREQkCMxsZ0nl6uITEZGQpAQlIiIhSQlKRERCUlATlJmNNLMpGRkZwQxDRERCUFATlHPuQ+fcuNjY2GCGISIiIUhdfCIiEpKUoEREJCQpQYmISEgKiQd1feGcY9a673h16Q5aN6jJwDZx9G9Tn3o1ooIdmoiI+KBCJ6iU/cf484z1LEo5QIt61Vm/+wivf/UtZtCpSW0GJsYxoE19klrVo3pUhf5RRUSqnAr5v3bmyVyemZvCiwu3US0ynL+O6swNfVuS7xxrUzNYnJLOwpR0pi3awZT524gMN3q2qMvANnEMalufbvF1iAxX76aISCizYG75bmYjgZGJiYm3JScnn/P8U915D8/8hr0ZJ/hJr3h+P6IDcTWjSzz/eHYey3ccZNHWdBalpLNhzxGcgxpR4fRtXZ8BbeozqG0c7RvVwsz8/NOJiIgnzGylc673GeXBTFCn9O7d251rsdiU/cd4cMYGFqak06lJbR66sjO9Wtbz6nMOZWazdNuBwoR1gO3pmQDE1Yyif5s4Brapz8DEOJrXq17mn0VERLxTYRNU8e68mMhw7r2kPTf0bUl4mO8tnj2Hj7MoJZ3FWw+wMCWdtKMnAWhRrzoDE+szoE3BGFb9UlpoIiLiuwqXoJxzzF7/HQ99VNCdd02veCaepTvPV845UvYfY1FKOou2HmDp1gMcPZkLQMcmtblraBtGdm8akM8WEanKSktQITlJYmtaQXfeguR0OjapzTNjetI7wbvuPG+ZGW0b1aJto1rcPLAVuXn5rNudweKtB/ho7V5+Mf1r0o+d5JaBrQIah4iIFAipBJWVXdCdN3VBQXfeX67ozA19WxARhBl3EeFh9GxRl54t6jJ2UCt+9cbX/OXDbziUmc2EH7XTpAoRkQALmQQ1a91eHv7oG/YUduf9fngHGtQKjbGfmMhwnrv+PO5/bz2T56ZwIDObv47q4pdxMBERKVlIJKjt6Znc9doqOjapzeRy6M4ri4jwMCZd3ZW6NaJ44cutHD6ewz+u7UFUhJ6nEhEJhJBIUFnZefxtZCdu7NcyKN15njIzJo7oQL0akTw6axNHjufwwo29qBEdErdRRKRSCYls0L5xwcSEUE5OxY0b3IYnrunGopR0bpj6FYcys4MdkohIpRMSGSGiAo7lXNu7Oc/f2Itv9h7h2n8t4buME8EOSUSkUvF7gjKzFmY2w8xeMrOJ/r5+KLmkc2P+c0sSezNOcPXzi9mWdizYIYmIVBoeJajCZLPfzNafVj7czDabWUqxZNQOmOmcuxXo5Od4Q07/NvV5Y1w/TuTk8ZMXlrB+d0awQxIRqRQ8bUFNA4YXLzCzcOA5YAQFiWiMmXUCvgZGm9lcYJ7/Qg1dXZrF8r87+hMTGc7oKUtZsvVAsEMSEanwPEpQzrn5wMHTipOAFOfcNudcNvAGMAq4Bfizc+5C4LLSrmlm48xshZmtSEtLK1v0IaR1g5q8fWd/msTGcNPLy5iz4btghyQiUqH5MgbVDNhV7H1qYdnHwC/N7AVgR2mVnXNTnHO9nXO9GzRo4EMYoaNJbDXeur0/nZrU5s7/ruSt5bvOXUlERErkS4Iqaeqdc86td85d45y7wzn327NewGykmU3JyKg84zZ1a0Tx2s/7MjAxjt+9s5Zxr6wgZf/RYIclIlLh+JKgUoHmxd7HA3u8uYBz7kPn3LjY2Fgfwgg9NaIjePGmPvzmR+1YvPUAF/9jPhPfWaup6CIiXvAlQS0H2ppZKzOLAkYDM7y5QGVsQZ0SFRHGLy5qy5f3DuXmAa14Z1UqQ56cx+MfbyLjeE6wwxMRCXke7QdlZtOBoUAcsI+CSRAvmtmlwD+BcOAl59wjZQnCkx11K7pdB7P4+6dbeH/1bmrHRHL3BW34Wf8EYiLDgx2aiEhQVbgNCyurDXsyeOLjzXy5JY2msTFM+FE7rjovXiuji0iVFZIJysxGAiMTExNvS05ODlocwbA4JZ1JH29ibWoGbRvW5J5h7RjRpTFhSlQiUsWEZII6pSq1oIpzzjFr3Xf8/dPNbE3LpEPjWvzqorZc0lmJSkSqjtISVEgsFltVmRmXdWvCJxOG8PToHmTn5XPna6u4dPICPl7/HaHwy4OISLCoiy+E5OU7ZqzZzeTPU9ienkmnJrWZ8KN2DOvYUFvMi0ilpS6+CiQ3L58PVu9h8txkdh7I4mf9W/KXKzorSYlIpVRagtJWsCEoIjyMq3vFM6pHUx6bvYkXF24nJjKcP4zooCQlIlWGElQIiwgP44HLOpKTl8+U+duIiQjj1xe3D3ZYIiLlIqgJqtgYVDDDCGlmxoMjO3MyJ5/Jc1OIjgzn7gt0v0Sk8gvqLL7Kuhafv4WFGY9e1ZVRPZry5JzNvLhwe7BDEhEJOHXxVRDhYcbfftKd7Nx8HvroG2Iiw7ihb8tghyUiEjB6DqoCiQgP4+nRPbmwQ0Puf289T83ZzJETWnhWRCqnoCaoyryaeaBERYTxfzecx6geTXl2XgqDn5jH819sJSs7N9ihiYj4lZ6DqsDW787gb59sZt7mNOJqRnP3BW0Yk9RCK6SLSIWipY4qoS7NYnn5liTevqM/iQ1r8JcPv2H4P+ezcufBYIcmIuIzJahKoHdCPd4Y15//3JpETp7jmheW8OisjZzIyQt2aCIiZaYEVYkMadeAORMGM7pPC6bM38blzyxkza7DwQ5LRKRMlKAqmZrRETx2VVdeuTWJzJO5XPX8Yt5dlRrssEREvKZZfJXU4MLWVL/W9fjN/9bw2lc7zzhnb8Zx0o6eDEJ0IiLnppUkKrHaMZG8eFMfLmhf8NzU1AXbAEg/dpIH3l/HoMfnMejxuTzx8SaO6nkqEQkxmmZeBWTn5jPhzdXMXLeXS7s2ZsGWdLJy8rg+qQXHTuby3te7iasZxe8u6cBPesdrxXQRKVfabqMKi4oI4+nRPYiJDOedVakM69iQiSM6ktiwJgC3DEzg4Y828rt31rJ8x0EeurKLnqUSkaBTC6oKcc6x+/Bx4utWP+NYfr7jn58nM/nzZLo3r8MLN55Hk9hqQYhSRKoaPagrmFmJyQkKVkz/9Y/a8cKNvUjZd5Th/1zAG8u+JT8/+L/AiEjVpAQlPzC8S2M+GD+I9o1qMfHddfzkX0vYeSAz2GGJSBWkBCVnSGxYkzdv78eT13Qjed9Rbn91JSdzC1alyMt3fLR2Dwczs4McpYhUdnoOSkpkZvykd3P+fm0PNn13lH98moxzjgfeX8/417/mwr99wXR1AYpIAOk5KDmrYZ0aMSapOf+av5Xx079m+rJvubFfC9o1qsUf3l3H3z7dHOwQRaSSUhefnNMDl3Wied3qzFy7l9F9mvPQqC68Oa4fo3o05d/zt7PrYFawQxSRSkgJSs6pRnQE//5Zb34/vAMPX9kFM8PMmDiiA2Fh8PjHm4IdoohUQkpQ4pH2jWtx59A2RIR//5VpEluNcYPb8NHavcxYs4fcvPwgRigilY0SlPjk9sGtSahfnV9O/5p+j81l8db0YIckIpWEEpT4pEZ0BHMmDOaFG3tRIzqc+99bT45aUiLiB0pQ4rPoiHCGd2nMny7vxPb0TKYv+7bo2Pb0TPYdORHE6ESkolKCEr+5sEND+reuzz8/S+ZQZjY70jO5fPICLnjqC15ZsoNQWPdRRCoOJSjxGzPj/ss6cuxELlc/v5i7X19FeJhxXou6/OmDDcxYsyfYIYpIBaIEJX7VpVks//15Xw5lZbNhzxEmXd2NV25NomuzWB6btYms7NxghygiFYTfE5SZnW9mL5jZVDNb7O/rS+hLalWPD38xiKk/682lXZsQFmY8eEUnvjtygsdmbVJXn4h4xKMNC83sJeByYL9zrkux8uHA00A4MNU5N8k5twBYYGZXAsv9H7JUBPF1q/9ga49eLetxy8AEXl60gx0HMomvW40b+rakSzMtcyUiJfO0BTUNGF68wMzCgeeAEUAnYIyZdSp2yvXAdD/EKJXEny7vxB9GdGDd7gzeXpnKE3M2cyInj0Up6Vp0VkTO4FGCcs7NBw6eVpwEpDjntjnnsoE3gFEAZtYCyHDOHSntmmY2zsxWmNmKtLS0skUvFYqZcfuQNqz+08XcNTSRBclpjHt1JTdM/YqfvvQVx05qfEpEvufLGFQzYFex96mFZQBjgZfPVtk5N8U519s517tBgwY+hCEV0eik5hgwf0saA9rUZ1HKAT5YvTvYYYlICPElQVkJZQ7AOfdn59w5J0hoP6iqq0lsNYZ3aUzruBq8dHMf4mpGsXLHoTPO27LvKItStHySSFXkS4JKBZoXex8PePWgi/aDqtr+eV1PZv7yfGIiw+nVsi4rdh7ikw3fserb7xPV2P8s54apX/2gTESqBo9m8ZViOdDWzFoBu4HRFEyM8JiZjQRGJiYm+hCGVFRREd//ftS7ZT3mbNjHna+tItyMHs3r0LxedXYdPA7AhDdXc9fQNsTXrc7AxLhghSwi5cijFpSZTQeWAO3NLNXMxjrncoHxwBxgI/CWc26DNx+uFpSc0iuhLgDREWF0bFKLdbszeGdVKgDX923BzgNZ/P6dddww9SvN+BOpIjxqQTnnxpRSPguY5deIpErq0jSWRrWjuWlAAncOacN3R04wcNJc8h3cM6wt0RFhvLxoBwCt75vFzQMSWLc7g+euP4/GsTHBDV5EAsKC+VR/sS6+25KTk4MWh4SGvHxHeNj3c29uf3UFOw9k8fE9gwE4mJnNeQ99+oM6wzo24t8/64VZSXN2RKQiMLOVzrneZ5SHwrIzvXv3ditWrAh2GBJisrJzycl1xFaPLCqbvW4vLy/awbIdB6kRFU5mdh7PXt+Ty7s15UROHmYF23+ISMVRWoLSYrESsqpHRfwgOQGM6NqEO4a2BuDeS9rTuWltnvh4M/n5jqRHPmPkMwuDEaqIBEBQE5Seg5KyGNKuIZOu6sropBbcMaQN3x7MYs6G7zhyIpct+44FOzwR8ZOgJijN4pOyCA8zRie1ICYynEs6N6ZBrWh+/daaouOa5SdSOaiLTyq0qIgwrk9qwfGcvKKy1vfN4paXl7E29XDwAhMRnylBSYV304AEftyzGTf2a1FUNm9zGlc8u4iVOw+Rn+9wzvHB6t0cPZETxEhFxBsag5IKr16NKP5xXQ+uT2p5xrGrn1/M1IXb+Gzjfn71xmp+Of1rMrVqukiFoDEoqTSa16sGwM8HteKVW5OKyj9cs5dXl+4EClpWfR/9nHYPzOaIWlMiIc2XtfhEQkqtmEhWPDCMetWjCCv2wO+63T9soZ/ad2p7Wibdm9ch9VAWsdUiqRXzwyntIhJcGoOSSiWuZnRRcrq+7/djUp2b1ubeS9r/4NyM4zmcyMlj0OPzGPLkF+UZpoh4IKgtKK1mLoH06I+78tuL2/PZxn1c1KEh9WtGc/RELi98uRWAZ+Yms7xwD6qDmdms3HmQ7vF1iAjX720ioUBLHUmVsnHvEUY8vaDU43WrR7Lg9xdSM1q93yLlRUsdiQAdm9RmzZ8vLvX4oawcJr6zlpy8/HKMSkRKogQlVU5stYLJENERYUQXbpr424vbFR3/aO1e2t4/m0OZ2aQeyuLml5eRkaUZfyLlTf0YUiWteGAYMZHhVI8MZ1t6JokNazKqRzNeXLidaYt3ADBjzR6+/vYQX2xO4+MNe7muT4uzX1RE/EotKKmS4mpGUzM6grAwI7FhTQCa16vO2EGtis5ZsvUAR08UTEn/aO1eAA5nZXPg2MnyD1ikCtIsPpFi4utWK3r98Ybvil4vSE7nZG4ePf76KTGRYWx6aEQwwhOpUrSShEgxZsa0W/rw92u7F5UltaoHwOWTC/aaOpGjCRQi5UHTzEVKkZ/vOJmbz4HMkwx6fN4Pjp3Xog5dm8USWy2S/yzZedaZgSJydqVNM9ckCZFShIUZ1aLCiY+qzq0DW/HSou1Fx1Z9e5hV3x4uep+RlXPG7r8i4htNkhDxwNW9mp31ePe/flK0QoWI+IcSlIgHWtavUfT6Vxe15cWbzuiNYNLsTdz33jo+Xv/dGcdExHtKUCIeKL700YQftaNR7ZgSz3v9q2+5478rWZt6mNe/+pYt+46WV4gilY7GoEQ89MHdAwmzgpXSW9avDsBDV3bhj++vP+PcK55dVPR6x6TLyidAkUpGO+qKeKh78zp0jS94JKJWTCTbH7uUn/Y7cxff0+XnO3akZwY6PJFKR89BiZSRFbamOjWpfdbzWt83i6FPfcG0YrMAReTcNAYl4qOrziuY4bfsvovOet6DH35TtM28VksXOTeNQYn4aOygVoxOauHRHlLdHvyk6LUZLPjdBew8kEW3+FhtOS9yGiUoER+Z2RnJacekyzh2Mpc7Xl3JwpT0Eus5R9EKFcM6NmJqCVPXRaoydfGJ+NHEER0YlBgHFExNn3R1V4/qfbZxXyDDEqmQtBafSIBlHM+hRlQ4iffPPut59wxry90XJBIZrt8bpWrRWnwiQXJqB99T/ndHf37ywpIzzvvnZ8mcyMnncFY2sdUj+deX27i0a2PuGNKGbvF1yilakdChFpRIOZmxZg/N6lSjV8u6pB09yb++3MrUhZ5NPf9kwmDaNaoV4AhFgqO0FpT6EkTKyRXdm9KrZV0AGtSK5oHLO/Hq2CSP6n6XcYIDx04SCr9QipQXdfGJBNH5bRt4dN7PXlpW9HrFA8NYvv0gI7o2CVRYIiHB7wnKzMKAh4DawArn3H/8/Rkilcm1veN5f/UesnM9e3i398OfAXB+2zheHds3kKGJBJVHXXxm9pKZ7Tez9aeVDzezzWaWYmYTC4tHAc2AHCDVv+GKVD5PXNOdLQ+P8HpR2QXJ6Xy0dk+AohIJPk/HoKYBw4sXmFk48BwwAugEjDGzTkB7YIlz7tfAnf4LVaTymziiAw1qRTP7V+czYVi7c54//vWvyyEqkeDwKEE55+YDB08rTgJSnHPbnHPZwBsUtJ5SgUOF5+SVdk0zG2dmK8xsRVpamveRi1RCdwxpw/L7h9GxSW1+Nawtr//8+y68kd2blljn7tdXsaiU1SpEKjJfxqCaAbuKvU8F+gJPA8+Y2fnA/NIqO+emAFOgYJq5D3GIVFoDEuNo27AmTetU46FRnflwzZldejPX7mXm2r0A/G54e2pFR/DT/gnlHKmI//mSoKyEMuecywLGenQBs5HAyMTERB/CEKncPv31EI/PfeLjzQCEh4Vx33vruP/Sjtw2uHWgQhMJKF+eg0oFmhd7Hw94NWKr/aBEvDPnnsEenXffe+sAeGTWRrKycwMZkkjA+JKglgNtzayVmUUBo4EZ/glLRErSvnEtZv/qfP5xXXeP6/T466ccz87jRE6pQ8IiIcmjpY7MbDowFIgD9gF/ds69aGaXAv8EwoGXnHOPePXh33fx3ZacnOxl6CKSMHGmx+dGhYcxOqk5dwxpQ9M61QIYlYh3SlvqSGvxiVRgU+Zv5dFZm7yu5+0zVyKBFJJr8ZnZSDObkpGREcwwRCqsn/ZL4OeDWvH+3QO9qvfYrI1c8o/5vLMylQuf+oLj2er+k9CjFpRIJXEoM5tPN+4DB797Z63X9Tf+dTjVosIDEJnI2YVkC0pE/KdujSiu7d2ca/s0P/fJJXj3a61MJqFFXXwilVDj2jFFr5vGxpzlzO9Fhun3VQktQf1G6jkokcBYet9FRa/7tq7vUZ1v9h7hmz1HWJCcRsbxnECFJuIx/cokUknN++1QHvlxF6xwzZc/Xd6JsJLWfyk0bfEOLp28gJ++uIzuf/mEXQezeHLOJvLzgz9OLVWTNiwUqaRaxdWgVVwNsk7mAbvp1LQ2iyZeSP/H5npU//wn5gFQIzqCu4ZqOTIpfxqDEqnkxg5qxYzxA+nXuj5NYr1/QPeJjzcz5Ml5fLPnSACiEymdppmLVDFb9h0lN8+xPT2Tu19f5V3dh0cQZhARrtEB8Z/Sppmri0+kimnXqBYAnZrW5r2vG1E9KpwZJWzjUWLdB2YDsPbBi6kdExmwGEVAkyREqrSpN/Vm8pieXtfr9uAnHMrMZsnWA1otXQJGY1Aiwo5JlzGsYyOv6vR86FPG/Hspf3h3XYCikqpOz0GJCAB/HtmJC9o34LcXtyO+bjXCzzYnvZj5W9ICHJlUVZokISIlcs5x67TlzNvseQJ6Y1w/Nuw5wthBrQIYmVQ2WotPRLxiZrx8SxI3D0jwuDU1espSHvroG5ZuOxDg6KQqUIISkbN68IrObH30Ut67a4DHdUZPWcrm744GMCqpCtTFJyJe8WYX3/6t63M8J4/37hqAmWetMKl6QrKLT7P4RCqeP13eyeNzl2w7wOpdh2n1h1kkTJxJ2tGTAYxMKhvN4hMRr9wyMIF7L2lfprp9HvlMSUo8pjEoEfGKmXHX0DZF7z3db+qUPo98RsLEmXyliRRyDhqDEhGfncjJY9Szi9i8z7uJETsmXRagiKQiCckxKBGpHGIiw5kzYbDXCSdh4kxe+2pngKKSik4JSkSC6v731pMwcSYX/u2LYIciIUYJSkT86sWbejP+Au83ONyWlsmxk7lsSzsWgKikItIYlIgEhHOOdbszuOLZRV7XHdWjKU+P9n6VdamYQnIMSs9BiVReZka3+DrsmHQZ8++9gPPbxnlc94PVe1iz6zC5efmcyMkLYJQSytSCEpFycSInjw9W7+b373i+PUfbhjVJ3n9Ms/0quZBsQYlI1RETGc51fVp4VSd5f8F41K3TlvPrt1YHICoJZUpQIlKu/ujFUkmnzN20n3dX7WbP4eMs2aoHfKsKJSgRKVe3DEgoev3wlV28qjtg0lzG/HspD87Y4OeoJBQpQYlIuTKD6/u24K3b+3Njv5ZlGl+atngHCRNnMmvd3gBEKKFCCUpEypWZ8eiPu5LUql5R2fTb+lGneqTX17rrtVXsPnyck7ma6VcZaRafiISMS59ewDd7j5S5/uQxPbmie1M/RiTlQbP4RCTkzRg/kC/vHVrm+r+c/jX7j54gOzfff0FJ0ChBiUjIiAgPo2X9Grx1e/8yXyPpkc9p98BszfarBJSgRCTkJLWqx//u6M8Hdw8EoGeLOl5fY8y/l/o5Kilvfk9QZjbUzBaY2QtmNtTf1xeRqqFPQj26Ny9YKum9uwaW6RoJE2dy4Jh28K2oPEpQZvaSme03s/WnlQ83s81mlmJmEwuLHXAMiAFS/RuuiFRVTWJj6B4f63W93o98RvK+o+TmaVyqovFoFp+ZDaYg6bzinOtSWBYObAF+REEiWg6MATY55/LNrBHwd+fcDee6vmbxiYgnMo7n0P0vn5S5/ke/GESXZt4nOQksn2bxOefmAwdPK04CUpxz25xz2cAbwCjn3KlfUw4B0WcJaJyZrTCzFWlpaR79ECJStcVWi/Rp4djLn1nI0RM5foxIAinCh7rNgF3F3qcCfc3sKuASoA7wbGmVnXNTgClQ0ILyIQ4RqWL+O7Yv7329myMncvj0m31e1e364Ce0iqtBXr5j7m+GEBGuuWKhypcEZSWUOefcu8C7Hl3AbCQwMjHR+903RaTqGtQ2jkFt48jPd/z8lRXM3bTfq/rb0zOBgkVoL+7cOBAhih/48qtDKtC82Pt4YI83F3DOfeicGxcbqz5hEfFeWJhx++DWRe8HtKnvVf1xr67kxYXbWbnzkL9DEz/wpQW1HGhrZq2A3cBo4Hq/RCUi4qGkVvWYOKIDo/s0JyYynA5//Nir+g999E3R61duTWJwuwb+DlHKyNNp5tOBJUB7M0s1s7HOuVxgPDAH2Ai85Zzzag18bfkuIr4yM+4Y0oY61aOIiQxn8cQLy3ytn720TAvPhhAtFisilYpzjufmpfDUJ1to06AGW9Myvb7GG+P60aVZLDERYZpEUQ5Km2Ye1ARVbJLEbcnJyUGLQ0Qqr+zcfNo9MLtMdaMjwpj326E0rVPNz1FJcSG5mrkmSYhIoEVFhJHYsGaZ6p7MzWfApLkkTJzJln1H/RyZnIvariJS6X326yE+PeALMM/Lqeziu6AmKE2SEJHytP4vl5S57mOzNxW9zs3LL3qWSgJHXXwiUmXUjPblyRoY+cxC9h89wZNzNnPBU1+w62CWnyKTkvj2tyUiUsGs+dPFHDmRw7NzU3hzxa5zVyhm3e4Mkh75vOh92rGTNK9X3d8hSiGNQYlIlRJbPZLm9aoz6equ3Dm0DWv+dDEJ9cuWZLJO6pmpQNI0cxGp8vLzHe3/OJucPO//P7y0a2O6xdfh9sGtMStpiVI5l5B8DuoUPagrIsGWn+/IOJ5Dz4c+LVP9f1zXnR/3jPdzVFVDSD4HJSISKsLCjLo1otgx6TLGJDU/d4XTTHhzDQAnc/M4nq2uP3/QJAkRkdM8dlU3pi/zbgIFQKs/zORUp5Svz12JWlAiIiWac89gPhw/yKs6xUdMVu48RH5+8IdQKjI9qCsiUoL2jWvRNT6W7Y9dWqb6Vz+/mKc/T+ZEjrr8ykqTJEREzuHn/1nBZxu921r+dJ//ZghtGpRtTcDKTrP4RET8IGHizDLX3fTQcGIiw/0YTeWgWXwiIn7w/t0Dy1y3wx8/5sstaX6MpnJTghIR8UKP5nWoEVX2VtBNLy0j43gO/126k5y8fD9GVvloJQkRkTLIyMohJz+f6/+9lC37jpXpGjcPSODBKzr7ObKKJyS7+LSauYhUVLHVI4mrGc0Hd3s3Fb24aYt3EArzAEKVHtQVEfFBtahw/nh5J95avovNZdh1t9UfZhW9TnlkBBHhGnk5RXdCRMRHYwe1Ys6EwT5f5+iJXD9EU3koQYmI+ElEWMFq5r1b1i1T/Xx19/2AnoMSEfGTU0sbhYUZb69M5bf/W1Pma/VvXZ/p4/r5K7SQFpKTJEREKpOwMCOssBV1Ta94hnVsWOZrLdl2wF9hVVhKUCIiATL1pj7MGF/2B3sTJs5k1HOLquyis1osVkQkgLrF1+H12/oCcMvABK/rr9l1mNb3zTr3iZWQxqBERMpJXr6jjQ/JpmlsDA9c3olLuzbxY1TBV9oYlJ6DEhEpJ4XDU2W2J+MEd722CoDWcTV4f/xAasdE+iGy0KQxKBGRcmJmPHlNN79ca1t6Jt0e/ISEiTPZnp7pl2uGGiUoEZFydE2veB6/uiv/d8N5frvmjVO/Yvfh4367XqhQghIRKUdmxnV9WnBp1yZlmjRRkt2HjzNw0lza3DeLyZ8XLLx9MDObrOyKvTKFJkmIiATRml2HGfXcooBcu12jmiTUr8En3+xjx6TLAvIZ/qAHdUVEQlD35nWY+cuyr4h+Nlv2HeOTbwq2qj+RkxeQzwgktaBEREJI+rGTTJq9ibdXpvr92nN/M4TWDWoCMG/zfiLCjEGJcZj5OL3QR6W1oJSgRERCUF6+48CxkyQ9+rnfr/3Li9oWjVUBLPz9BcTXre73z/GUuvhERCqQ8DCjfs3ogFy7eHICGPT4PBImzmTDnu9X9Vmbepi8IC+xFJAHdc2sBjAf+LNz7qNAfIaISGUXHmb8747+pB09WfSAbiBdNnkhMZFhNKodw84DWUXltw5sxRU9mlKvehQt6pdfS8ujLj4zewm4HNjvnOtSrHw48DQQDkx1zk0qLP8rkAls8CRBqYtPROTsFqWks3LnIf7+6ZagxjG8c2PWpB7mlxe1ZUxSC79c06cxKDMbDBwDXjmVoMwsHNgC/AhIBZYDY4CmQBwQA6QrQYmI+M+FT33BthBZOaJ7fCwfjPd9BqJPa/E55+abWcJpxUlAinNuW+EHvAGMAmoCNYBOwHEzm+Wcyy8hoHHAOIAWLfyThUVEKrvZ95xP1sk8Xl/2LfF1q/GrN1YHLZY1qYHdicKXMahmwK5i71OBvs658QBmdjMFLagzkhOAc24KMAUKWlA+xCEiUmVER4QTHRHO3RckAgXPN/3+nXVBjiowfJnFV9LE+aJE45ybdq7uPe0HJSLim+v6tGDLwyOCHUZA+JKgUoHmxd7HA3u8uYBz7kPn3LjY2FgfwhARqdqiIsLYMekyhrZvEOxQ/MqXBLUcaGtmrcwsChgNzPBPWCIi4q1ptyQFOwS/8ihBmdl0YAnQ3sxSzWyscy4XGA/MATYCbznnNnjz4eriExHxrx2TLuPiTo0Yk9SC567335YewaCljkREKinnHDPX7WV458a8sXwXX397mHdW+XeNP3+skh6Sa/GZ2UhgZGJi4m3JycnnPF9ERHyXMHGm364VyAQVkKWOPOWc+xD4sHfv3rcFMw4Rkapk+f3DcM4RExXO0q0HWPntIf715bZgh3WGoCYoEREpfw1qfb8I7cWdG3NRx0ZFCer1n/dlTWoGuw9n8d+l3wYrRCDICapYF18wwxARqdLCw4xxg1szZf422jeuxYDEOAAevrIrC5PT6de6HnszTnD+E/OK6tw+pDVXnxcf0Lg0SUJERAr2n8o8ScNaMeX+2doPSkREShUeZkFJTmejBCUiIiEpqAlKD+qKiEhpgpqgtBafiIiURl18IiISkpSgREQkJClBiYhISNIkCRERCUmaJCEiIiEpJFaSMLM0YGeALh8LlKWJ5mm9s51X2rGSyj0pK/4+Dkj3IL6yKst988c9O9txb+9RSe8Ded8C+V0L5D0rqawqfNe8Pabv2tmP+/Jdq+OcO3M7YOdcpf4DTAlkvbOdV9qxkso9KSv+HlgRavfNH/fM1/vmwfuA3bdAftcCec/Odd8q63fN22P6rgX+u3b6n6owSeLDANc723mlHSup3JOysv4sZVGWz/LHPTvb8bLco1C/Z57WC+Q9K6ks1O9bIP99lnZM37WzH/f7dy0kuvjEe2a2wpWwuKKcne6b93TPykb3zXdVoQVVWU0JdgAVlO6b93TPykb3zUdqQYmISEhSC0pEREKSEpSIiIQkJSgREQlJSlAiIhKSlKAqCTOrYWb/MbN/m9kNwY6nojCz1mb2opm9HexYKgozu7Lwe/aBmV0c7HgqAjPraGYvmNnbZnZnsOOpKJSgQpiZvWRm+81s/Wnlw81ss5mlmNnEwuKrgLedc7cBV5R7sCHEm/vmnNvmnBsbnEhDh5f37P3C79nNwHVBCDckeHnPNjrn7gCuBfRslIeUoELbNGB48QIzCweeA0YAnYAxZtYJiAd2FZ6WV44xhqJpeH7fpMA0vL9nDxQer6qm4cU9M7MrgIXA5+UbZsWlBBXCnHPzgYOnFScBKYW/+WcDbwCjgFQKkhRU8b9XL++b4N09swKPA7Odc6vKO9ZQ4e33zDk3wzk3AFAXvIeq9H9kFVQzvm8pQUFiaga8C1xtZs9TvmuCVRQl3jczq29mLwA9zewPwQktZJX2XfsFMAy4xszuCEZgIay079lQM5tsZv8CZgUntIonItgBiNeshDLnnMsEbinvYCqQ0u7bAUD/yZastHs2GZhc3sFUEKXdsy+AL8o3lIpPLaiKJxVoXux9PLAnSLFUJLpv3tM9857umR8pQVU8y4G2ZtbKzKKA0cCMIMdUEei+eU/3zHu6Z36kBBXCzGw6sARob2apZjbWOZcLjAfmABuBt5xzG4IZZ6jRffOe7pn3dM8CT6uZi4hISFILSkREQpISlIiIhCQlKBERCUlKUCIiEpKUoEREJCQpQYmISEhSgpIqx8zyzGx1sT8JwY7JX8ysp5lNLXx9s5k9e9rxL8ys1O0ezOwNM2sb6DhFPKG1+KQqOu6c61HSATMzCp4PzC/fkPzmPuBhH+o/D/wOuM0/4YiUnVpQUuWZWYKZbTSz/wNWAc3N7F4zW25ma83sL8XOvb9wM7rPzGy6mf22sLyoZWJmcWa2o/B1uJk9WexatxeWDy2s87aZbTKz1wqTI2bWx8wWm9kaM1tmZrXMbIGZ9SgWxyIz63baz1EL6OacW+PBz3xFsRbkZjPbXnhoATDMzPTLqwSdvoRSFVUzs9WFr7cDE4D2wC3OubusYBvzthTs7WPADDMbDGRSsLZaTwr+7awCVp7js8YCGc65PmYWDSwys08Kj/UEOlOwmOgiYKCZLQPeBK5zzi03s9rAcWAqBTvY3mNm7YBo59za0z6rN7D+tLLrzGxQsfeJULA3EYVrxJnZW8CXheX5ZpYCdPfgZxMJKCUoqYp+0MVXOAa10zm3tLDo4sI/Xxe+r0lBwqoFvOecyyqs58kioBcD3czsmsL3sYXXygaWOedSC6+1GkgAMoC9zrnlAM65I4XH/wf80czuBW6lYDfX0zUB0k4re9M5N77Yz/pF8YNm9jsK7kfxnXH3A01RgpIgU4ISKZBZ7LUBjznn/lX8BDO7Byht8cpcvu8yjzntWr9wzs057VpDgZPFivIo+PdoJX2Gcy7LzD6lYHfWayloLZ3u+GmffVZmdhHwE2DwaYdiCq8lElQagxI50xzgVjOrCWBmzcysITAf+LGZVSsc7xlZrM4OoFfh62tOu9adZhZZeK12ZlbjLJ+9CWhqZn0Kz69VbDxoKgUbBS53zp2+1TgUrJ6d6MkPaGYtgf8DrnXOnZ6M2gFagVuCTi0okdM45z4xs47AksJ5C8eAG51zq8zsTWA1sJOCCQWnPAW8ZWY/BeYWK59KQdfdqsJJEGnAlWf57Gwzuw54xsyqUdCSGQYcc86tNLMjwMul1N1kZrFmVss5d/QcP+bNQH3gvcKfcY9z7lIza0RBl9/ec9QXCThttyFSRmb2IAWJ46ly+rymFGwb3qG0afBmNgE46pybWsbPmAAccc69WOZARfxEXXwiFYCZ/Qz4Crj/HM9oPc8Px7a8dRj4jw/1RfxGLSgREQlJakGJiEhIUoISEZGQpAQlIiIhSQlKRERCkhKUiIiEpP8HB+9/RHc37h8AAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "spectrum.plot_power()\n", "decorate(xlabel='Frequency (Hz)',\n", " xscale='log', \n", " yscale='log')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And the slope is close to -1." ] }, { "cell_type": "code", "execution_count": 22, "metadata": { "jupyter": { "outputs_hidden": false } }, "outputs": [ { "data": { "text/plain": [ "-1.0016799964612606" ] }, "execution_count": 22, "metadata": {}, "output_type": "execute_result" } ], "source": [ "spectrum.estimate_slope().slope" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true, "jupyter": { "outputs_hidden": true } }, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.8.5" } }, "nbformat": 4, "nbformat_minor": 4 }