{
"cells": [
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
""
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# In Search of the Fourth Wave\n",
"\n",
"### Allen Downey\n",
"\n",
"Olin College\n",
"\n",
"[DSP Online Conference](https://www.dsponlineconference.com/)"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "skip"
}
},
"source": [
"Copyright 2020 Allen B. Downey\n",
"\n",
"License: [Creative Commons Attribution 4.0 International](http://creativecommons.org/licenses/by/4.0/)"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"Run this notebook\n",
"\n",
"[tinyurl.com/mysterywave](https://tinyurl.com/mysterywave)"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"When I was working on [*Think DSP*](https://greenteapress.com/thinkdsp), I encountered a small mystery. \n",
"\n",
"As you might know:"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"* A sawtooth wave contains harmonics at integer multiples of the fundamental frequency, and\n",
"\n",
"* Their amplitudes drop off in proportion to $1/f$. "
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"* A square wave contains only odd multiples of the fundamental, but \n",
"\n",
"* They also drop off like $1/f$. "
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"* A triangle wave also contains only odd multiples, \n",
"\n",
"* But they drop off like $1/f^2$."
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"Which suggests that there's a simple waveform that \n",
"\n",
"* Contains all integer multiples (like a sawtooth) and \n",
"\n",
"* Drops off like $1/f^2$ (like a triangle wave). \n",
"\n",
"Let's find out what it is."
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"In this talk, I'll \n",
"\n",
"* Suggest four ways we can find the mysterious fourth wave.\n",
"\n",
"* Demonstrate using tools from SciPy, NumPy and Pandas, and\n",
" "
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"And a tour of DSP and the topics in *Think DSP*."
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"I'm a professor of Computer Science at Olin College.\n",
"\n",
"Author of *Think Python*, *Think Bayes*, and *Think DSP*.\n",
"\n",
"And *Probably Overthinking It*, a blog about Bayesian probability and statistics.\n",
"\n",
"Web page: [allendowney.com](http://allendowney.com)\n",
"\n",
"Twitter: @allendowney"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"This talk is a Jupyter notebook.\n",
"\n",
"[You can read it here](https://nbviewer.jupyter.org/github/AllenDowney/ThinkDSP/blob/master/code/fourth_wave.ipynb).\n",
"\n",
"[And run it here](https://colab.research.google.com/github/AllenDowney/ThinkDSP/blob/master/code/fourth_wave.ipynb)."
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"Here are the libraries we'll use."
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"slideshow": {
"slide_type": "-"
}
},
"outputs": [],
"source": [
"import numpy as np\n",
"import pandas as pd\n",
"import matplotlib.pyplot as plt"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "skip"
}
},
"source": [
"And a convenience function for decorating figures."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"slideshow": {
"slide_type": "skip"
}
},
"outputs": [],
"source": [
"def decorate(**options):\n",
" \"\"\"Decorate the current axes.\n",
"\n",
" Call decorate with keyword arguments like\n",
"\n",
" decorate(title='Title',\n",
" xlabel='x',\n",
" ylabel='y')\n",
"\n",
" The keyword arguments can be any of the axis properties\n",
"\n",
" https://matplotlib.org/api/axes_api.html\n",
"\n",
" In addition, you can use `legend=False` to suppress the legend.\n",
"\n",
" And you can use `loc` to indicate the location of the legend\n",
" (the default value is 'best')\n",
" \"\"\"\n",
" plt.gca().set(**options)\n",
" plt.tight_layout()"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"slideshow": {
"slide_type": "skip"
}
},
"outputs": [],
"source": [
"def legend(**options):\n",
" \"\"\"Draws a legend only if there is at least one labeled item.\n",
"\n",
" options are passed to plt.legend()\n",
" https://matplotlib.org/api/_as_gen/matplotlib.plt.legend.html\n",
"\n",
" \"\"\"\n",
" ax = plt.gca()\n",
" handles, labels = ax.get_legend_handles_labels()\n",
" if handles:\n",
" ax.legend(handles, labels, **options)"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "skip"
}
},
"source": [
"## Basic waveforms"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"We'll start with the basic waveforms:\n",
"\n",
"* sawtooth, \n",
"\n",
"* triangle, and \n",
"\n",
"* square."
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"Sampled at CD audio frame rate."
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"slideshow": {
"slide_type": "-"
}
},
"outputs": [],
"source": [
"framerate = 44100 # samples per second\n",
"dt = 1 / framerate # seconds per sample"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"At equally-spaced times from 0 to `duration`."
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {
"slideshow": {
"slide_type": "-"
}
},
"outputs": [],
"source": [
"duration = 0.005 # seconds\n",
"ts = np.arange(0, duration, dt) # seconds"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"We'll work with signals at 1000 Hz. \n",
"\n",
"The number of complete cycles is $f t$."
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [],
"source": [
"freq = 1000 # cycles per second (Hz)\n",
"cycles = freq * ts # cycles"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"In 0.005 seconds, a 1000 Hz signal completes 5 cycles."
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"4.988662131519274"
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"np.max(cycles)"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {
"slideshow": {
"slide_type": "skip"
}
},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAAApFUlEQVR4nO3deXhU9dn/8fdNIKyBsK8J+74KAdcqqBRZ3FqsW11qq7XW+ljbWpfaWhArtIpY5UHq8qjValtrLQEjogJqXUFcSMIewg4BAoHsk/v3x4y/phggQE5mknxe18WVzJkz53ufL5N8cs59ZsbcHRERkVhTL9oFiIiIVEQBJSIiMUkBJSIiMUkBJSIiMUkBJSIiMUkBJSIiMUkBJXKMzCzLzM6Ndh2VYWavmdk10a5D5HgooKRWM7MrzOwTMztgZtsiv7DPiHZdQTCze83sz+WXuft4d38mWjWJnAgFlNRaZnYb8DBwP9AeSAZmAxdGsSwRqSQFlNRKZtYCmAL82N3/4e4H3b3E3ee5+y/MrIOZ5ZtZ63KPGWFmu8ysQeT29WaWYWZ5ZpZuZsMrGKeemd1hZuvMbLeZ/dXMWkXua2Rmf44szzWzj82sfQXbuMPM/n7Isllm9kjk+2vNbH2kjg1mdmUF2zgPuAu4NHK0+Flk+WIz+0G57bxnZjMj9aw3s9MiyzeZ2c7ypwPNrKGZ/cHMss1sh5nNMbPGx/P/IXI8FFBSW50KNAJeqehOd98OLAa+U27xd4EX3b3EzC4B7gWuBpoDFwC7K9jULcBFwFlAJ2Av8FjkvmuAFkAS0Bq4ESioYBt/ASaYWXMAM4uL1PWCmTUFHgHGu3sCcBqwooL9SSN8pPiSuzdz96EV7TdwMvB5pJ4XgBeBkUCvyP4/ambNIutOB/oAwyL3dwZ+fZjtilQ5BZTUVq2BHHcvPcI6zxD+pfxVKFwOPBe57wfADHf/2MPWuvvGCrbxQ+Bud9/s7kWEQ22ymdUHSiJ19HL3kLsvc/f9h24gst3lhIMO4Gwg390/iNwuAwaZWWN33+buKys7CRXY4O5Pu3sIeIlweE5x9yJ3XwgUA73MzIDrgZ+6+x53zyMcgJedwNgix0QBJbXVbqBNJCgO51VggJn1AMYC+9z9o8h9ScC6SozTFXglcsosF8gAQoR7Xs8BrwMvmtlWM5vx1enDCrxAOCABrojcxt0PApcSPvraZmbzzaxfJeo6nB3lvi+IjHHosmZAW6AJsKzcvqVFlotUCwWU1FbvA4X856jka9y9EPgrcCVwFf85egLYBPSsxDibCJ9+Syz3r5G7b4n0vH7r7gMIn5qbRPiUYUX+Bow2sy7AxUQCKlLn6+4+FugIZAJ/OtwuVaLeysohHFYDy+1XC3dvdrQHilQVBZTUSu6+j3C/5DEzu8jMmphZAzMbb2Yzyq36LHAt4R5T+Uu0nwB+Hrlwwsysl5l1rWCoOcC0r+4zs7ZmdmHk+zFmNjhy+nA/4VN+ocPUu4twT+xpwqfhMiLbaG9mF0R6UUXAgcNtg/DRUTczO+Gfa3cvIxyEM82sXaSWzmY27kS3LVJZCiiptdz9IeA24FfALsJHOzcD/yy3znuEezzL3T2r3PK/AdMIH8nkRR7TqoJhZgH/AhaaWR7wAeELEQA6AH8nHE4ZwBL+OwQP9QJwLuWOngj/jP4M2ArsIXwxxk2HefzfIl93m9nyI4xTWb8E1gIfmNl+YBHQtwq2K1Ippg8slLrOzN4CXnD3J6Jdi4j8hwJK6jQzGwm8ASRFrlQTkRihU3xSZ5nZM4RPW92qcBKJPTqCEhGRmKQjKBERiUlHehFjtWvTpo1369Yt2mWIiEg1WrZsWY67f+1F4DEVUN26deOTTz6JdhkiIlKNzKyitxHTKT4REYlNCigREYlJCigREYlJCigREYlJCigREYlJgV7FZ2ZZhN9oMwSUuntKkOOJiEjtUR2XmY9x95xqGEdERGoRneITEZFjlr07n5eXbQ50jKCPoJzw5+Q48Li7zz10BTO7AbgBIDk5OeByRETkRBwoKuWxt9fy5DsbaNSgHt8c2J6ERg0CGSvogDrd3bdGPpHzDTPLdPel5VeIhNZcgJSUFL1zrYhIDAqVOS8v28yM11eRc6CIb53UmdvP6xdYOEHAAeXuWyNfd5rZK8AoYOmRHyUiIrHkow17+O28lazcup/hyYk8cU0Kw5ISAx83sIAys6ZAPXfPi3z/TWBKUOOJiEjV2rQnnwdey2T+F9vo2KIRsy4bxgVDO2Fm1TJ+kEdQ7YFXIjtSn/BHaqcFOJ6IiFSBg0WlzF68lj+9s4F6Bree25sfntmTxvFx1VpHYAHl7uuBoUFtX0REqlZZmfPy8nCfaVdeERcN68Tt5/WjU2LjqNQTUx+3ISIi0fFJ1h6mpKbz+eZ9DE1K5PGrRjA8uWVUa1JAiYjUYVtyC3jgtUzmfbaV9s0b8tB3hnLRsM7Uq1c9faYjUUCJiNRB+cWlzFm8jseXrgfglrN7cePonjSJj51YiJ1KREQkcGVlzj9XbGF6WiY79hdx/tBO3DG+H52j1Gc6EgWUiEgdsWzjXqakpvPZplyGdGnBY1cMJ6Vbq2iXdVgKKBGRWm5rbgHT0zJ5dcVW2iU05A+XDOVbJ8VGn+lIFFAiIrVUQXGIOUvW8fjSdZQ53DymFz8a3ZOmDWvGr/6aUaWIiFSau/Pqiq1MT8tk275CJg7pyB3n9SOpVZNol3ZMFFAiIrXIp9nhPtOn2bkM6tycWZedxKjusdtnOhIFlIhILbB9XyHT0zJ55dMttE1oyIzJQ5g8vEvM95mORAElIlKDFRSHmLt0PXOWrCPkzk2je3LTmF40qyF9piOp+XsgIlIHuTvzPt/GAwsy2LqvkPGDOnDXhP41rs90JAooEZEa5rNNuUxJTWfZxr0M6Nichy4dxik9Wke7rCqngBIRqSF27C9kRtoqXl6+mTbN4nngW4O5JCWJuBrcZzoSBZSISIwrLAnxxDvrmb14HaUh54dn9eDmMb0C/bj1WKCAEhGJUe7O/C+28bsFmWzJLWDcwPbcNaE/XVs3jXZp1UIBJSISg77cso8p89L5KGsP/Tok8MIPTua0Xm2iXVa1UkCJiMSQnXmF/D5tFX9fvplWTeK5/+LBXDqy9vaZjkQBJSISAwpLQjz57gZmv72W4lAZPzijOz85pzfNa3mf6UgUUCIiUeTupH25nftfy2DTngLO7d+euyf2p3ubutFnOhIFlIhIlHy5ZR9TU9P5cMMe+rZP4M/fP5kzetetPtORKKBERKrZrrwiHly4ipc+2URi4wZMvWgQl49Mon5cvWiXFlMUUCIi1aSoNMTT72Xx6FtrKSwJcd3p3bnlnN60aFx3+0xHooASEQmYu/P6yh3cvyCD7D35nNOvHXdP7E+Pts2iXVpMU0CJiAQoY9t+psxL5/31u+ndrhnPXjeKM/u0jXZZNYICSkQkADkHinhw4Wpe+jib5o0bMOXCgVwxKll9pmOggBIRqULFpWU88+8sHnlzDfklIa4+tRu3ntubxCbx0S6txlFAiYhUAXdnUcZOps1PJ2t3PqP7tuVXE/vTq11CtEursRRQIiInKHP7fu5LzeDdtTn0bNuUp783kjF920W7rBpPASUicpx2Hyhi5qLVvPBhNgmNGvCb8wfw3VO60kB9piqhgBIROUbFpWU8+34Ws95cQ35xiKtO6cqt5/ahZVP1mapS4AFlZnHAJ8AWd58U9HgiIkFxd97K3Mm0+RmszznIN3q34Z5JA+jTXn2mIFTHEdT/ABlA82oYS0QkEKt35DE1NZ131uTQo01Tnro2hTF922FW9z4Go7oEGlBm1gWYCEwDbgtyLBGRIOw9WMzMRat5/sNsmsbHcc+kAVx1Slfi66vPFLSgj6AeBm4HDnv8a2Y3ADcAJCcnB1yOiEjllITKeO79jTy8aDUHikq58uSu/HRsH1qpz1RtAgsoM5sE7HT3ZWY2+nDruftcYC5ASkqKB1WPiEhlvZ25k6nz01m/6yBn9Ar3mfp2UJ+pugV5BHU6cIGZTQAaAc3N7M/u/t0AxxQROW5rd+YxNTWDJat30b1NU564OoVz+qvPFC2BBZS73wncCRA5gvq5wklEYlFufjEPL1rDcx9spEl8HL+a2J+rT+2mPlOU6XVQIlJnlYbKeP7DbGYuWs3+ghIuG5XMz8b2oXWzhtEuTaimgHL3xcDi6hhLRKQylqzexX2p6azZeYDTerbmnkkD6N9Rr4aJJTqCEpE6Zd2uA0ybn8FbmTvp2roJj181gm8OaK8+UwxSQIlInbAvv4RZb67h2fezaNQgjjvH9+Pa07vRsH5ctEuTw1BAiUitVhoq4y8fZfPQG6vJLSjhspFJ3Da2L20T1GeKdQooEam13l2Tw9TUdFbtyOPk7q349fkDGNipRbTLkkpSQIlIrbMh5yDT5qezKGMnSa0aM+e7wxk3sIP6TDWMAkpEao19BSX88c01PPN+FvFx9bj9vL5cd3p3GjVQn6kmUkCJSI0XKnNe/DibBxeuZm9+MZeM6MLPx/WlXUKjaJcmJ0ABJSI12r/X5jAlNZ3M7XmM6hbuMw3qrD5TbaCAEpEaKSvnIPcvyGBh+g46JzbmsSuGM2Gw+ky1iQJKRGqUvMISHn1rLU+/l0X9OOMX4/ry/TPUZ6qNFFAiUiOEypy/fbKJPyxcRc6BYiaP6MLt4/rSrrn6TLWVAkpEYt4H63czZV466dv2k9K1JU9dO5IhXRKjXZYETAElIjEre3c+9y/IIG3ldjonNuaPl5/EpCEd1WeqIxRQIhJzDhSV8tjba3nynQ3E1TNuG9uHG87soT5THaOAEpGYESpzXl62mRmvryLnQBHfOqkzt5/Xjw4t1GeqixRQIhITPtqwh9/OW8nKrfsZnpzIE9ekMCwpMdplSRQpoEQkqjbtyeeB1zKZ/8U2OrZoxKzLhnHB0E7qM4kCSkSi42BRKbMXr+VP72ygnsGt5/bmh2f2pHG8+kwSpoASkWpVVua8vDzcZ9qVV8RFwzpx+3n96JTYONqlSYxRQIlItfk4aw9T5qXzxZZ9DE1K5PGrRjA8uWW0y5IYpYASkcBtyS3gdwsySP18Gx2aN2LmpUO5cGhn6tVTn0kOTwElIoHJLy5lzuJ1PL50PQC3nNObG8/qQZN4/eqRo9OzRESqXFmZ888VW5ielsmO/UVcMLQTvxzfj87qM8kxUECJSJVatnEvU1LT+WxTLkO6tGD2lcMZ0bVVtMuSGkgBJSJVYmtuAdPTMnl1xVbaJTTkwUuGcvFJ6jPJ8VNAicgJKSgOMWfJOh5fuo4yh5vH9OJHo3vStKF+vciJ0TNIRI6Lu/Pqiq1MT8tk275CJg7pyJ3j+9GlZZNolya1hAJKRI7Zp9nhPtOn2bkM6tycWZedxKju6jNJ1VJAiUilbd9XyPS0TF75dAttExoyY/IQJg/voj6TBEIBJSJHVVAcYu7S9cxZso6QOzeN7slNY3rRTH0mCZCeXSJyWO7OvM+38cCCDLbuK2T8oA7cNaE/Sa3UZ5LgBRZQZtYIWAo0jIzzd3f/TVDjiUjV+mxTLlNS01m2cS8DOjbnoUuHcUqP1tEuS+qQII+gioCz3f2AmTUA3jWz19z9gwDHFJETtGN/ITPSVvHy8s20aRbP9G8PZvKIJOLUZ5JqFlhAubsDByI3G0T+eVDjiciJKSwJ8cQ765m9eB2lIeeHZ/Xg5jG9SGjUINqlSR0VaA/KzOKAZUAv4DF3/7CCdW4AbgBITk4OshwRqYC7M/+LbfxuQSZbcgsYN7A9d03oT9fWTaNdmtRxgQaUu4eAYWaWCLxiZoPc/ctD1pkLzAVISUnREZZINfpyyz6mzEvno6w99OuQwAvXn8xpPdtEuywRoJqu4nP3XDNbDJwHfHmU1UUkYDvzCvl92ir+vnwzrZrEc//Fg7l0pPpMEluCvIqvLVASCafGwLnA9KDGE5GjKywJ8eS7G5j99lqKQ2Vc/40e3Hx2L5qrzyQxKMgjqI7AM5E+VD3gr+6eGuB4InIY7k7al9u5/7UMNu0pYOyAcJ+pexv1mSR2BXkV3+fASUFtX0Qq58st+5iams6HG/bQt30Cz//gZE7vpT6TxD69k4RILbUrr4gHF67ipU82kdi4AVMvGsTlI5OoH1cv2qWJVEqlAsrMZgD3AQVAGjAUuNXd/xxgbSJyHIpKQzz9XhaPvrWWwpIQ153enVvO6U2LxuozSc1S2SOob7r77WZ2MbAZuAR4G1BAicQId+f1lTu4f0EG2XvyOadfO+6e2J8ebZtFuzSR41LZgPrqT68JwF/cfY+ZLkcViRUZ2/YzZV4676/fTe92zXj2ulGc2adttMsSOSGVDah5ZpZJ+BTfTZFLyAuDK0tEKiPnQBEPLlzNSx9n07xxA6ZcOJArRiWrzyS1QqUCyt3vMLPpwH53D5lZPnBhsKWJyOEUl5bxzL+zeOTNNRSUhLjmtG78zzm9SWwSH+3SRKpMZS+SaAL8GEgm/L55nYC+gF7XJFKN3J1FGTuZNj+drN35jOnblrsnDqBXO/WZpPap7Cm+pwm/6etpkdubgb+hgBKpNpnb93Nfagbvrs2hZ9um/N/3RjK6b7tolyUSmMoGVE93v9TMLgdw9wLTVRIi1WL3gSJmLlrNCx9mk9CoAfeeP4ArT+lKA/WZpJarbEAVR95PzwHMrCfhDyQUkYAUl5bx7PtZzHpzDfnFIa4+NdxnatlUfSapGyobUL8h/ALdJDN7HjgduDaookTqMnfnrcydTJufwfqcg5zZpy33TOxP7/YJ0S5NpFpV9iq+N8xsOXAKYMD/uHtOoJWJ1EGrd+QxNTWdd9bk0KNtU56+diSj+7ZFZ9SlLjpiQJnZ8EMWbYt8TTazZHdfHkxZInXL3oPFzFy0muc/zKZpfBy/njSAq05Vn0nqtqMdQT14hPscOLsKaxGpc0pCZTz3/kYeXrSaA0WlXHlyV346tg+t1GcSOXJAufuY6ipEpK55O3MnU+ens37XQc7o1YZ7Jg2gbwf1mUS+UtkX6v4YeN7dcyO3WwKXu/vsAGsTqZXW7sxjamoGS1bvonubpjxxdQrn9G+nPpPIISp7Fd/17v7YVzfcfa+ZXQ8ooEQqKTe/mIcXreG5DzbSJD6OX03sz9WndiO+vvpMIhWpbEDVMzNz969eBxUH6CS5SCWUhsp4/sNsZi5azf6CEi4blczPxvahdbOG0S5NJKZVNqAWAn81szmEL464kfDrokTkCJas3sV9qems2XmA03q25p5JA+jfsXm0yxKpESobULcDPwR+RPh1UAuBJ4IqSqSmW7frANPmZ/BW5k66tm7C3KtGMHZAe/WZRI5BZQNqAvC4u/9vkMWI1HT78kuY9eYann0/i0YN4rhzfD+uPb0bDevHRbs0kRqnsgF1GTDLzF4Gnnb3jABrEqlxSkNl/OWjbB56YzW5BSVcNjKJ28b2pW2C+kwix6uyb3X0XTNrDlwOPG1mTvgjOP7i7nlBFigS695dk8PU1HRW7cjjlB6tuGfSAAZ2ahHtskRqvMoeQeHu+yNHUI2BW4GLgV+Y2SPu/seA6hOJWRtyDjJtfjqLMnaS1Koxc747nHEDO6jPJFJFKvtC3fOB64CewHPAKHffGfmk3QxAASV1xr6CEv745hqeeT+L+Lh6/PK8fnzv9G40aqA+k0hVOtqbxfYCOgCXADPdfWlk+TfMLMHd15nZddVQp0jUhcqcFz/O5sGFq9mbX8x3RiTxs3F9aJfQKNqlidRKRzuCehi4y92vPmR5QeS+8939zQDqEokp/16bw5TUdDK35zGqWyt+ff4ABnVWn0kkSEcLqG7u/vmhC939EzPrFkxJIrEjK+cg9y/IYGH6DjonNuaxK4YzYbD6TCLV4WgBdaRzF42rshCRWJJXWMKjb63l6feyqB9n/GJcX75/Rnf1mUSq0dEC6mMzu97d/1R+oZl9H1gWXFki0REqc/72ySb+sHAVOQeKmTyiC7eP60u75uoziVS3owXUrcArZnYl/wmkFMJvFHtxgHWJVLsP1u9myrx00rftJ6VrS566diRDuiRGuyyROutoH1i4AzjNzMYAgyKL57v7W0fbsJklAc8SvgqwDJjr7rNOsF6RKpe9O5/7F2SQtnI7nRMb88fLT2LSkI7qM4lEWWXfSeJt4O1j3HYp8DN3X25mCcAyM3vD3dOPtUiRIBwoKuWxt9fy5DsbiKtn/GxsH64/s4f6TCIxotLvJHGs3H0bsC3yfZ6ZZQCdAQWURFWozHl52WZmvL6KnANFfGt4Z24f148OLdRnEoklgQVUeZFL0k8CPqzgvhuAGwCSk5Oroxypwz7asIffzlvJyq37GZ6cyBPXpDAsKTHaZYlIBQIPKDNrBrwM3Oru+w+9393nAnMBUlJSPOh6pG7atCefB17LZP4X2+jYohGzLhvGBUM7qc8kEsMCDSgza0A4nJ53938EOZZIRQ4WlTJ78Vr+9M4G6hn89Nw+3HBmDxrHq88kEusCCygL/2n6JJDh7g8FNY5IRcrKnJeXh/tMu/KKuGhYJ345vh8dW+j15SI1RZBHUKcDVwFfmNmKyLK73H1BgGOK8HHWHqbMS+eLLfsYlpTI41eNYHhyy2iXJSLHKMir+N4FdIJfqs2W3AJ+tyCD1M+30aF5I2ZeOpQLh3amXj09DUVqomq5ik8kSPnFpcxZvI7Hl64H4JZzenPjWT1oEq+nt0hNpp9gqbHKypx/rtjC9LRMduwv4oKh4T5T50T1mURqAwWU1EjLNu5lSmo6n23KZUiXFsy+cjgjuraKdlkiUoUUUFKjbM0tYHpaJq+u2Eq7hIY8eMlQLj5JfSaR2kgBJTVCQXGIOUvW8fjSdZQ53DymFz8a3ZOmDfUUFqmt9NMtMc3deXXFVqanZbJtXyETh3TkzvH96NKySbRLE5GAKaAkZn2aHe4zfZqdy+DOLXjk8pMY2U19JpG6QgElMWf7vkKmp2XyyqdbaJvQkBmThzB5eBf1mUTqGAWUxIyC4hBzl65nzpJ1hNy5aXRPbhrTi2bqM4nUSfrJl6hzd+Z9vo0HFmSwdV8hEwZ34M7x/UlqpT6TSF2mgJKo+mxTLlNS01m2cS8DOjbnoUuHcUqP1tEuS0RigAJKomLH/nCf6R/Lt9CmWTzTvz2YySOSiFOfSUQiFFBSrQpLQjzxznpmL15Haci58aye/HhMTxIaNYh2aSISYxRQUi3cnflfbON3CzLZklvAuIHtuWtCf7q2bhrt0kQkRimgJHBfbtnHlHnpfJS1h34dEnjh+pM5rWebaJclIjFOASWB2ZlXyO/TVvH35Ztp1SSe+y8ezKUj1WcSkcpRQEmVKywJ8eS7G5j99lqKQ2Vc/40e3Hx2L5qrzyQix0ABJVXG3Un7cjv3v5bBpj0FjB0Q7jN1b6M+k4gcOwWUVIkvt+xjamo6H27YQ9/2CTz/g5M5vZf6TCJy/BRQckJ25RXx4MJVvPTJJlo2iee+iwZx2cgk6sfVi3ZpIlLDKaDkuBSVhnj6vSwefWsthSUhvn96d35yTm9aNFafSUSqhgJKjom78/rKHdy/IIPsPfmc278dd03oT4+2zaJdmojUMgooqbT0rfuZmprO++t306d9M569bhRn9mkb7bJEpJZSQMlR5Rwo4sGFq3np42xaNG7A1AsHcvmoZPWZRCRQCig5rOLSMp75dxaPvLmGgpIQ15zWjVvP6UOLJuoziUjwFFDyNe7OooydTJufTtbufMb0bcvdEwfQq536TCJSfRRQ8l8yt+/nvtQM3l2bQ8+2Tfm/741kdN920S5LROogBZQAsPtAETMXreaFD7NJaNSAe88fwJWndKWB+kwiEiUKqDquuLSMZ9/PYtaba8gvDnH1qd34n3N607JpfLRLE5E6TgFVR7k7b2XuZNr8DNbnHOTMPm25Z2J/erdPiHZpIiKAAqpOWr0jj6mp6byzJocebZvy9LUjGd23LWb6GAwRiR2BBZSZPQVMAna6+6CgxpHK23uwmJmLVvP8h9k0jY/j15MGcNWp6jOJSGwK8gjq/4BHgWcDHEMqoSRUxnPvb+ThRas5UFTKlSd35adj+9BKfSYRiWGBBZS7LzWzbkFtXyrn7cydTJ2fzvpdB/lG7zb8auIA+nZQn0lEYl/Ue1BmdgNwA0BycnKUq6k91u7MY2pqBktW76J7m6Y8cXUK5/Rvpz6TiNQYUQ8od58LzAVISUnxKJdT4+XmF/PwojU898FGmsTH8auJ/bn61G7E11efSURqlqgHlFSN0lAZz3+YzcxFq9lfUMLlo5K5bWwfWjdrGO3SRESOiwKqFliyehf3paazZucBTuvZmnsmDaB/x+bRLktE5IQEeZn5X4DRQBsz2wz8xt2fDGq8umjdrgNMm5/BW5k76dq6CXOvGsHYAe3VZxKRWiHIq/guD2rbdd2+/BJmvbmGZ9/PonGDOO6a0I9rTutGw/px0S5NRKTK6BRfDVIaKuMvH2Xz0BuryS0o4bKRSdw2ti9tE9RnEpHaRwFVQ7y7Joepqems2pHHKT1acc+kAQzs1CLaZYmIBEYBFeM25Bxk2vx0FmXsJKlVY+Z8dzjjBnZQn0lEaj0FVIzaV1DCH99cwzPvZxEfV49fnteP753ejUYN1GcSkbpBARVjQmXOix9n8+DC1ezNL+Y7I5L42bg+tEtoFO3SRESqlQIqhvx7bQ5TUtPJ3J7HqG6t+PX5AxjUWX0mEambFFAxICvnIPcvyGBh+g66tGzM7CuHM36Q+kwiUrcpoKIor7CER99ay1PvbaBBXD1+Ma4v3z+ju/pMIiIooKIiVOb87ZNN/GHhKnYfLGby8C78Ylxf2jVXn0lE5CsKqGr2wfrdTJmXTvq2/aR0bcnT145icBf1mUREDqWAqibZu/O5f0EGaSu30zmxMY9ecRITB3dUn0lE5DAUUAE7UFTKY2+v5cl3NhBXz/jZ2D5cf2YP9ZlERI5CARWQUJnz8rLNzHh9FTkHivjW8M788rx+tFefSUSkUhRQAfhowx5+O28lK7fuZ3hyIk9ck8KwpMRolyUiUqMooKrQpj35PPBaJvO/2EbHFo2YddkwLhjaSX0mEZHjoICqAgeLSpm9eC1/emcD9Qx+em4fbjizB43j1WcSETleCqgTUFbmvLw83GfalVfERcM68cvx/ejYonG0SxMRqfEUUMfp46w9TJmXzhdb9jEsKZHHrxrB8OSW0S5LRKTWUEAdo817w32m1M+30aF5Ix6+NNxnqldPfSYRkaqkgKqk/OJS/nfxOuYuXQ/ALef05sazetAkXlMoIhIE/XY9irIy558rtjA9LZMd+4u4YGi4z9Q5UX0mEZEgKaCOYNnGvUxJTeezTbkM7dKC2VcOZ0TXVtEuS0SkTlBAVWBrbgHT0zJ5dcVW2iU05MFLhnLxSZ3VZxIRqUYKqHIKikPMWbKOx5euwx1+cnYvbjyrJ00bappERKqbfvMC7s6rK7YyPS2TbfsKmTSkI3eM70eXlk2iXZqISJ1V5wPq0+xwn+nT7FwGd27BI5efxMhu6jOJiERbnQ2o7fsKmZ6WySufbqFtQkN+P3kI3x7eRX0mEZEYUecCqqA4xNyl65mzZB0hd24a3ZObxvSimfpMIiIxpc78VnZ35n2+jQcWZLB1XyETBnfgzvH9SWqlPpOISCyqEwH12aZcpqSms2zjXgZ0bM5Dlw7jlB6to12WiIgcQa0OqB37w32mfyzfQptm8Uz/9mAmj0giTn0mEZGYF2hAmdl5wCwgDnjC3R8IcryvFJaEeOKd9cxevI7SkHPjWT358ZieJDRqUB3Di4hIFQgsoMwsDngMGAtsBj42s3+5e3pQY7o787/Yxu8WZLIlt4DzBnbgzgn96Nq6aVBDiohIQII8ghoFrHX39QBm9iJwIRBIQLk71zz9MUtX76J/x+b8/pIhnNazTRBDiYhINQgyoDoDm8rd3gycfOhKZnYDcANAcnLycQ9mZpzVpy3jB3XgOynqM4mI1HRBBlRFCeFfW+A+F5gLkJKS8rX7j8X3z+h+Ig8XEZEYUi/AbW8Gksrd7gJsDXA8ERGpRYIMqI+B3mbW3czigcuAfwU4noiI1CKBneJz91Izuxl4nfBl5k+5+8qgxhMRkdol0NdBufsCYEGQY4iISO0U5Ck+ERGR46aAEhGRmKSAEhGRmKSAEhGRmGTuJ/Ta2CplZruAjSe4mTZAThWUU9toXiqmefk6zUnFNC8Vq4p56erubQ9dGFMBVRXM7BN3T4l2HbFG81IxzcvXaU4qpnmpWJDzolN8IiISkxRQIiISk2pjQM2NdgExSvNSMc3L12lOKqZ5qVhg81LrelAiIlI71MYjKBERqQUUUCIiEpNiOqDM7DwzW2Vma83sjgruNzN7JHL/52Y2/GiPNbNWZvaGma2JfG1ZXftTVQKal0vMbKWZlZlZjbyUNqB5+b2ZZUbWf8XMEqtpd6pMQPMyNbLuCjNbaGadqmt/qkIQc1Lu/p+bmZtZm6D3o6oF9Fy518y2RJ4rK8xsQqULcveY/Ef4IzrWAT2AeOAzYMAh60wAXiP86b2nAB8e7bHADOCOyPd3ANOjva8xMi/9gb7AYiAl2vsZQ/PyTaB+5Pvper78/3lpXu7xtwBzor2v0Z6TyP1JhD9iaCPQJtr7GgvzAtwL/Px4aorlI6hRwFp3X+/uxcCLwIWHrHMh8KyHfQAkmlnHozz2QuCZyPfPABcFvB9VLZB5cfcMd19VfbtR5YKal4XuXhp5/AeEPxm6JglqXvaXe3xToCZdbRXU7xaAmcDt1Kz5+EqQ83JcYjmgOgObyt3eHFlWmXWO9Nj27r4NIPK1XRXWXB2Cmpearjrm5TrCfz3WJIHNi5lNM7NNwJXAr6uw5qAFMidmdgGwxd0/q+qCq0mQP0M3R04JPnUsbZVYDiirYNmhf5Ucbp3KPLam0rxULNB5MbO7gVLg+eOqLnoCmxd3v9vdkwjPyc3HXWH1q/I5MbMmwN3UrKA+VFDPlf8FegLDgG3Ag5UtKJYDajPh87lf6QJsreQ6R3rsjsghKZGvO6uw5uoQ1LzUdIHNi5ldA0wCrvTISfUapDqeLy8A3z7hSqtPEHPSE+gOfGZmWZHly82sQ5VWHqxAnivuvsPdQ+5eBvyJ8OnAyol2Y+4IDbv6wHrC/+lfNd0GHrLORP67YffR0R4L/J7/vkhiRrT3NRbmpdxjF1MzL5II6vlyHpAOtI32PsbYvPQu9/ifAH+P9r5Ge04OeXwWNe8iiaCeKx3LPf6nwIuVrinak3KUCZsArCZ8dcjdkWU3AjdGvjfgscj9X5T/xVrRYyPLWwNvAmsiX1tFez9jZF4uJvxXUBGwA3g92vsZI/OylvC59RWRfzXmarWA5+Vl4Evgc2Ae0Dna+xntOTlk+1nUsIAK8LnyXGTdz4F/US6wjvZPb3UkIiIxKZZ7UCIiUocpoEREJCYpoEREJCYpoEREJCYpoEREJCYpoESOg5m1LvfuzNvLvVvzATObHdCYt5rZ1Ue4f5KZ/TaIsUWiQZeZi5wgM7sXOODufwhwjPrAcmC4/+fNaw9dxyLrnO7u+UHVIlJddAQlUoXMbLSZpUa+v9fMnol8XlKWmX3LzGaY2RdmlmZmDSLrjTCzJWa2zMxe/+qtuA5xNrD8q3Ays1vMLD3yBpwvAnj4r83FhN+WSaTGU0CJBKsn4beHuRD4M/C2uw8GCoCJkZD6IzDZ3UcATwHTKtjO6cCycrfvAE5y9yGEX+n/lU+Ab1T5XohEQf1oFyBSy73m7iVm9gXhD3VLiyz/AuhG+EMiBwFvhM/QEUf4HZ8P1RHIKHf7c+B5M/sn8M9yy3cCNerTbUUORwElEqwiAHcvM7MS/0/Tt4zwz58BK9391KNspwBoVO72ROBM4ALgHjMbGDn91yiyrkiNp1N8ItG1CmhrZqcCmFkDMxtYwXoZQK/IOvWAJHd/m/CntyYCzSLr9SH8Jq4iNZ4CSiSKPPzx2JOB6Wb2GeF3TD+tglVfI3zEBOHTgH+OnDb8FJjp7rmR+8YA84OsWaS66DJzkRrCzF4Bbnf3NYe5vz3wgrufU72ViQRDASVSQ5hZX6C9uy89zP0jgRJ3X1GthYkERAElIiIxST0oERGJSQooERGJSQooERGJSQooERGJSQooERGJSf8PEj4pHt7AD3YAAAAASUVORK5CYII=\n",
"text/plain": [
"