{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Cacophony for the whole family\n",
"==============================\n",
"\n",
"Allen Downey\n",
"\n",
"This is an example that demonstrates some of the features in the *Think DSP* library.\n",
"\n",
"It is inspired by the performance of a grade school band I witnessed recently. My goal is to simulate the sound of a beginner band.\n"
]
},
{
"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": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"First, a function that translates from a MIDI number to a frequency:"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"def midi_to_freq(midi_num):\n",
" \"\"\"Converts MIDI note number to frequency.\n",
"\n",
" midi_num: int MIDI note number\n",
" \n",
" returns: float frequency in Hz\n",
" \"\"\"\n",
" x = (midi_num - 69) / 12.0\n",
" freq = 440.0 * 2**x\n",
" return freq"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now here's a randomized version that simulates three kinds of errors: poor tuning, playing the wrong note, and [popping an overtone](https://en.wikipedia.org/wiki/Overtone). Notice that it is possible to make all three errors."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"import random\n",
"\n",
"def random_freq(midi_num):\n",
"\n",
" # simulate poor tuning by adding gaussian noise to the MIDI number\n",
" midi_num += random.gauss(0, 0.5)\n",
" \n",
" # one kid out of 10 plays the wrong note\n",
" if random.random() < 0.1:\n",
" midi_num += random.randint(-5, 5)\n",
" \n",
" freq = midi_to_freq(midi_num)\n",
" \n",
" # and one kid in 10 pops an overtone\n",
" if random.random() < 0.1:\n",
" freq *= random.randint(2, 5)\n",
"\n",
" return freq"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This function takes a MIDI number and duration and makes a Wave:"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"from thinkdsp import SawtoothSignal\n",
"\n",
"def make_note(midi_num, duration, framerate=22050):\n",
" \"\"\"Make a MIDI note with the given duration.\n",
"\n",
" midi_num: int MIDI note number\n",
" duration: float seconds\n",
" sig_cons: Signal constructor function\n",
" framerate: int frames per second\n",
"\n",
" returns: Wave\n",
" \"\"\"\n",
" freq = random_freq(midi_num)\n",
" signal = SawtoothSignal(freq)\n",
" wave = signal.make_wave(duration, framerate=framerate)\n",
" wave.apodize()\n",
" return wave"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let's test make_note. MIDI number 60 is middle C."
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
" \n",
" "
],
"text/plain": [
""
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"note = make_note(60, 1.0)\n",
"note.make_audio()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Sounds good.\n",
"\n",
"Now we can make 10 notes and play them at the same time. Since Wave provides `__add__`, we can use `sum`:"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [],
"source": [
"def make_ensemble(midi_num, duration):\n",
" notes = [make_note(midi_num, duration) for i in range(10)]\n",
" ensemble = sum(notes)\n",
" ensemble.make_audio()\n",
" return ensemble"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"And we can test it with a middle C:"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
" \n",
" "
],
"text/plain": [
""
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"c = make_ensemble(60, 1.0)\n",
"c.make_audio()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Good, sounds like angry bees.\n",
"\n",
"And now, a rousing chorus of that old crowd favorite, _Hot Cross Buns_.\n",
"\n",
"Wave provides `__or__`, which concatenates notes, so we can use `reduce`:"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
" \n",
" "
],
"text/plain": [
""
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"from thinkdsp import Wave\n",
"from functools import reduce\n",
"\n",
"midi_nums = [64, 62, 60, 64, 62, 60, 60, 60, 60, 60, 62, 62, 62, 62, 64, 62, 60]\n",
"durations = [0.5, 0.5, 1.0, 0.5, 0.5, 1.0, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.5, 0.5, 1.0]\n",
"\n",
"waves = [make_ensemble(midi_num, duration) for midi_num, duration in zip(midi_nums, durations)]\n",
"wave = reduce(Wave.__or__, waves)\n",
"wave.make_audio()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"And that brings a tear of pride to any parent's eye.\n",
"\n",
"On a more serious note, this example tells us something about how the ear interprets complex sounds with many tones and harmonics.\n",
"\n",
"Let's take a look at the spectrum of that middle C:"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3deXxV1bn/8c8DyKTIrCKDoFKpWq2YOqC3ivTiUKdabbVWaWuvHbxea3vb63D92Vu11dbetrbWoYpS660DdaBOiPOMBBABmYJMYUqYhwAhyfP7Y68TTsJJckLOsE/yfb9e55W91157n7Wyk/Octfbaa5u7IyIiEjft8l0AERGRVBSgREQklhSgREQklhSgREQklhSgREQkljrkuwC51qdPHx88eHC+iyEiIsHUqVPXuHvf+ultLkANHjyY4uLifBdDREQCM1uSKl1dfCIiEksKUCIiEksKUCIiEksKUCIiEksKUCIiEksKUCIiEksKUCIiEksKUCIiEksKUBl075sLWbq2It/FEBFpFRSgMqR88w5uf3Eul42dnO+iiIi0CgpQGZJ4MnFFZXWeSyIi0jooQImISCwpQImISCwpQImISCwpQGWI57sAIiKtjAJUhlm+CyAi0kooQImISCwpQImISCwpQImISCwpQImISCwpQGWIaxifiEhGKUBlmGkYn4hIRihAiYhILClAiYhILClAiYhILClAiYhILClAZYhrNj4RkYxSgMow02x8IiIZoQAlIiKxpAAlIiKxpAAlIiKxpACVYRosISKSGVkLUGY21szKzGxWUlovM5tkZgvCz54h3czsLjMrMbOPzWx40j5jQv4FZjYmKf1YM5sZ9rnLTJMMiYi0JtlsQT0MnFEv7TrgVXcfCrwa1gHOBIaG15XAPRAFNOBm4HjgOODmRFALea5M2q/+e+WFRvGJiGRG1gKUu78FrKuXfB4wLiyPA85PSv+rRz4AephZP+B0YJK7r3P39cAk4IywbV93f9/dHfhr0rFERKQVyPU1qP3dfSVA+LlfSO8PLEvKVxrSGksvTZGekpldaWbFZlZcXl7e4kqIiEj2xWWQRKp+Md+D9JTc/X53L3L3or59++5hEUVEJJdyHaBWh+45ws+ykF4KDEzKNwBY0UT6gBTpIiLSSuQ6QE0AEiPxxgDPJqVfHkbznQBsDF2AE4HRZtYzDI4YDUwM2zab2Qlh9N7lScfKCz1RV0Qkszpk68Bm9nfgVKCPmZUSjca7HXjCzK4AlgIXhewvAGcBJUAF8G0Ad19nZrcAU0K+X7h7YuDFD4hGCnYBXgyvvNNgdxGRzMhagHL3SxrYNCpFXgeuauA4Y4GxKdKLgSNbUkYREYmvuAySEBERqUMBSkREYkkBSkREYkkBKkM0iE9EJLMUoDJMg/hERDJDAUpERGJJAUpERGJJASpDxhdHc9du3LYzzyUREWkdFKAy5KH3FgGwtbI6zyUREWkdFKAyRHPxiYhklgKUiIjEkgKUiIjEkgJUhmgWcxGRzFKAyhDFJxGRzFKAEhGRWFKAEhGRWFKAEhGRWFKAEhGRWFKAEhGRWFKAEhGRWFKAEhGRWFKAEhGRWFKAEhGRWFKAEhGRWFKAEhGRWFKAEhGRWFKAyhA9r1BEJLMUoEREJJbyEqDM7Fozm21ms8zs72bW2cyGmNlkM1tgZo+bWceQt1NYLwnbBycd5/qQPs/MTs9HXUREJDtyHqDMrD/wH0CRux8JtAcuBu4AfufuQ4H1wBVhlyuA9e5+KPC7kA8zOzzsdwRwBvBnM2ufy7ok0/OgREQyK19dfB2ALmbWAegKrAROA8aH7eOA88PyeWGdsH2UmVlIf8zdd7j7IqAEOC5H5RcRkSzLeYBy9+XAncBSosC0EZgKbHD3qpCtFOgflvsDy8K+VSF/7+T0FPuIiEiBy0cXX0+i1s8Q4EBgb+DMFFkTA+NS9Z55I+mp3vNKMys2s+Ly8vLmF1pERHIuH118XwIWuXu5u+8EngJGAD1Clx/AAGBFWC4FBgKE7d2BdcnpKfapw93vd/cidy/q27dvpusjIiJZkI8AtRQ4wcy6hmtJo4BPgNeBC0OeMcCzYXlCWCdsf83dPaRfHEb5DQGGAh/mqA4iIpJlHZrOklnuPtnMxgPTgCpgOnA/8DzwmJndGtIeDLs8CDxiZiVELaeLw3Fmm9kTRMGtCrjK3atzWpkkUawVEZFMyXmAAnD3m4Gb6yV/SopReO6+HbiogePcBtyW8QKKiEjeaSaJDFH7SUQksxSgMkRz8YmIZJYClIiIxJIClIiIxJICVIboGpSISGYpQGWIrkGJiGSWApSIiMSSApSIiMSSApSIiMSSApSIiMSSApSIiMSSApSIiMSSAlSG6D4oEZHMUoDKEN0HJSKSWQpQIiISSwpQIiISS00GqPBo9pvM7C9hfaiZnZ39oomISFuWTgvqIWAHcGJYLwVuzVqJCpQGSYiIZFY6AeoQd/81sBPA3behz2MREcmydAJUpZl1IQxUM7NDiFpUIiIiWdMhjTw3Ay8BA83sUeAk4FvZLJSIiEiTAcrdJ5nZNOAEoq69a9x9TdZLVmBMnZ4iIhnVYIAys+H1klaGn4PMbJC7T8tesQqP605dEZGMaqwF9dvwszNQBMwgakEdBUwGTs5u0UREpC1rcJCEu49095HAEmC4uxe5+7HAMUBJrgpYKGrUhBIRyah0RvENc/eZiRV3nwV8PntFKkzrK3bmuwgiIq1KOqP45pjZA8DfiIaafxOYk9VSiYhIm5dOgPo28APgmrD+FnBP1kokIiJCesPMtwO/Cy8REZGcSGey2EVm9mn9V0ve1Mx6mNl4M5trZnPM7EQz62Vmk8xsQfjZM+Q1M7vLzErM7OPk4e9mNibkX2BmY1pSJhERiZd0uviKkpY7AxcBvVr4vn8AXnL3C82sI9AVuAF41d1vN7PrgOuA/wLOBIaG1/FE3YvHm1kvolkuioiujU01swnuvr6FZRMRkRhosgXl7muTXsvd/ffAaXv6hma2L/BF4MFw/Ep33wCcB4wL2cYB54fl84C/euQDoIeZ9QNOBya5+7oQlCYBZ+xpuUREJF6abEHVm1GiHVGLpVsL3vNgoBx4yMyOBqYSDcDY391XArj7SjPbL+TvDyxL2r80pDWUnqoOVwJXAgwaNKgFRRcRkVxJp4vvt0nLVcAi4GstfM/hwNXuPtnM/kDUndeQVLPceSPpuye63w/cD1BUVKQ7akVECkA6AeoKd68zKMLMhrTgPUuBUnefHNbHEwWo1WbWL7Se+gFlSfkHJu0/AFgR0k+tl/5GC8olIiIxks5MEuPTTEuLu68ClpnZYSFpFPAJMAFIjMQbAzwblicAl4fRfCcAG0NX4ERgtJn1DCP+Roc0ERFpBRqbzXwYcATQ3cwuSNq0L9Fovpa4Gng0jOD7lOhm4HbAE2Z2BbCUaLQgwAvAWUTz/1WEvLj7OjO7BZgS8v3C3de1sFwiIhITjXXxHQacDfQAzklK3wz8W0ve1N0/ou7w9YRRKfI6cFUDxxkLjG1JWUREJJ4aDFDu/izwrJmd6O7v57BMIiIijXbx/czdfw18w8wuqb/d3f8jqyUTEZE2rbEuvsSM5cW5KIiIiEiyxrr4/hl+jmsoj4iISLY01sX3Txq48RXA3c/NSolERERovIvvzpyVQkREpJ7GuvjeTCyH+5WGEbWo5rl7ZQ7KJiIibVg6k8V+GbgXWEg0/90QM/ueu7+Y7cKJiEjblc5UR78FRrr7qe5+CjASPV232Rat2crg657n/YVr810UEZGCkE6AKnP3kqT1T9k1kauk6YNPo8D07EfL81wSEZHCkM5s5rPN7AXgCaJrUBcBUxLz87n7U1ksn4iItFHpBKjOwGrglLBeTvTI93OIApYClIiIZFyTAcrdv52LgoiIiCRLZxTfEKLHYwxOzq8bdUVEJJvS6eJ7BngQ+CdQk93iiIiIRNIJUNvd/a6sl0RERCRJOgHqD2Z2M/AysCOR6O7TslYqERFp89IJUJ8DLgNOY1cXn4d1ERGRrEgnQH0FOFjz77WMNzgvvIiIpJLOTBIzgB7ZLkhbYZbvEoiIFIZ0WlD7A3PNbAq7rkG5u5+XvWKJiEhbl06Aujlp2YCTgUuyU5zWT119IiLpabKLLzwXaiPwZeBhYBTR4zekGdS1JyLSPI098v0zwMVEraW1wOOAufvIHJVNRETasMa6+OYCbwPnJB63YWbX5qRUIiLS5jXWxfdVYBXwupn9xcxGEV2DEhERyboGA5S7P+3uXweGAW8A1wL7m9k9ZjY6R+UTEZE2Kp1BElvd/VF3PxsYAHwEXJf1komISJuWzo26tdx9nbvf5+6a5qiZNLxcRKR5mhWgMsnM2pvZdDN7LqwPMbPJZrbAzB43s44hvVNYLwnbBycd4/qQPs/MTs9PTZpHw81FRNKTtwAFXAPMSVq/A/iduw8F1gNXhPQrgPXufijwu5APMzucaBj8EcAZwJ/NrH2Oyi4iIlmWlwBlZgOIbvx9IKwb0ezo40OWccD5Yfm8sE7YPirkPw94zN13uPsioAQ4Ljc1EBGRbMtXC+r3wM/Y9fiO3sAGd68K66VA/7DcH1gGELZvDPlr01PsU4eZXWlmxWZWXF5ensl6iIhIluQ8QJnZ2UCZu09NTk6R1ZvY1tg+dRPd73f3Incv6tu3b7PKKyIi+ZHOZLGZdhJwrpmdBXQG9iVqUfUwsw6hlTQAWBHylwIDgVIz6wB0B9YlpSck7xNbGs0nIpKenLeg3P16dx/g7oOJBjm85u6XAq8DF4ZsY4Bnw/KEsE7Y/pq7e0i/OIzyGwIMBT7MUTWaTaP3RESaJx8tqIb8F/CYmd0KTAceDOkPAo+YWQlRy+liAHefbWZPAJ8AVcBV7l6d+2KLiEg25DVAufsbRNMo4e6fkmIUnrtvBy5qYP/bgNuyV0IREcmXfN4H1abo2pOISPMoQOWYrkWJiKRHAUpERGJJAUpERGJJAUpERGJJAUpERGJJAUpERGJJAUpERGJJAUpERGJJASrHdMOuiEh6FKBERCSWFKByxMOjqjSThIhIehSgcuTGp2fluwgiIgVFASoHNm7bme8iiIgUHAWoHPjdpPn5LoKISMFRgMqB6hoN3RMRaS4FqBybtmRDvosgIlIQFKByIHnk3rzVm/NXEBGRAqIAJSIisaQAJSIisaQAJSIisaQAJSIisaQAFQNlm7bzm4lzqdFwdBGRWgpQOdDU9Hs/eXIGd7++kOIl63NSHhGRQqAAFQM7dtYAUKNncYiI1FKAihFNdC4isosCVIx8/f4PWL5hW76LISISCwpQMfDh4nW1y7OXb8xjSURE4iPnAcrMBprZ62Y2x8xmm9k1Ib2XmU0yswXhZ8+QbmZ2l5mVmNnHZjY86VhjQv4FZjYm13XJBl2FEhGJ5KMFVQX8xN0/C5wAXGVmhwPXAa+6+1Dg1bAOcCYwNLyuBO6BKKABNwPHA8cBNyeCWtyYHqMrItJsOQ9Q7r7S3aeF5c3AHKA/cB4wLmQbB5wfls8D/uqRD4AeZtYPOB2Y5O7r3H09MAk4I4dVyYpHJy/NdxFERGIhr9egzGwwcAwwGdjf3VdCFMSA/UK2/sCypN1KQ1pD6ane50ozKzaz4vLy8kxWIePemh/v8omI5EreApSZ7QP8A/iRu29qLGuKNG8kffdE9/vdvcjdi/r27dv8woqISM7lJUCZ2V5EwelRd38qJK8OXXeEn2UhvRQYmLT7AGBFI+kiItIK5GMUnwEPAnPc/X+TNk0AEiPxxgDPJqVfHkbznQBsDF2AE4HRZtYzDI4YHdJERKQV6JCH9zwJuAyYaWYfhbQbgNuBJ8zsCmApcFHY9gJwFlACVADfBnD3dWZ2CzAl5PuFu++6oUhERApazgOUu79Dw7P6jEqR34GrGjjWWGBs5konIiJxoZkkREQklhSgciAT9+m6O3e/XkLZ5u0tP5iISAHIxzWoNqelT9F4ongZZZu2c+fL83ljXhlPfn9EZgomIhJjClAF4GfjP65dnrJYDzUUkbZBXXw5oKn4RESaTwFKRERiSQFKRERiSQEqB6wFD3P/+YTZDW7btH0nY99ZxODrnmfwdc/z8uxVe/w+IiJxowAVcw+/t7jBbTc/O5tfPPdJ7fo1j33UYF4RkUKjAFXANlRU1lnftrM6TyUREck8BagCls6TetdvrWTrjqoclEZEJLN0H1QrVVPjXD72Q94pWQPA9WcOo3/PLpx91IF5LpmISHoUoFqpzduraoMTwK9enAvAFz/Tl30775WvYomIpE1dfDmQrRt19+Swv35pboPbdlbX7HlhREQyTAGqlVm2rqLR7RWVqQdSDLvpRYbe+CIfLtIjtUQkHhSgCliqltk3H5wcNjawU9LEtW/MK6u9h2r7zqj19P7CtVRV11C6vvFAJyKSbQpQrUxTM6fXJGUYl+IeK8e59fk5nHzH65Rv3lGbvmTtViqr1AUoIrmjANWGNRTL3l5QDsCNT8/k3ZI1XHTve5zymze4/qmZuSuciLR5ClAFbfd+PG8w7CS2N8591/1VL3+ymksfmFz7iI9/TCtlQ0UlE2ev4pnpyzn5jtfUqhKRrNEw8xzI5dM2Ej14TY0cPOP3bzF31eZmH/++tz7lnjcW1q7PX70ZdzjnT+/QzuDTX30ZgPcWruEbf5nMny8dzsjD9uOz/+8lABb+8izat9PzR5Jt3r4TgG4a/i9ShwJUDrTwgboN2pPh64kA1lhwauywFfVmpTj7j+/ULtd4NIpwYK+ufOMv0WCNHz46rU7+Q254gd9edDQ/eXIGZx/VjzOP7MdV/xflmfnz0Sk/pHdUVbOz2tmnU+v8c/3cz18GYPHtX85zSUTiRV18bUxNE6MoWhpMP1y0jm0NDGVP+O9nZgHw3Mcra4MTRB/U2yqr+eeMFfzgb1M59TevA3DOH9/hyJsnRuVzpyrpfq1tldWs3bKDF2eu1H1cIq1M6/xK2oY1NYovnQDUWMusqf1/8uQMbnn+k0bzNDap7T+mldYGsIT5q7cAsGVHFXe8OJdHPljCtJv+lV57d6ztOgQY3LsrI4ftR7/unbnyi4fUOcbT00uZumQ9t57/uSZqIJmyYsM2Rtz+Gs9dfTJH9u+e7+JIAVILKkP6de8MQJ99Ou62ramgsaf26EpOGqMkGnt+VTp12VCxs3llSjJhxooGt3133BQe+WAJAMNvmcRT00rrbF+8toKH3l3ML1+Yy7MfLQeia2GPfLCEax+fwd8+WEpVdQ13TpzHpnDdZ/vO6toW2fIN25iyeB3vL1xb57iL1mzl7tdL2FBRSUnZ5pTD82V3r88rA+DRyUvyXBIpVGpBZciBPbqwcuN29t+3827bmhpZl0txKks6pi5ZX7v8wad1Z7n48RMzGtzvmsc+4l+G9q29FpbwjQcm8+GidUxetJYnvz+CYTe9xL8M7cPow/fnpmd3PRzynkuHM2Xxesa+u4gDu3dmxcbt/GbivNrtY0YMBqCquoYrxhXz5vxyvnnCIL41Ygg9u+5F7306sW5rJV07tqfzXu1TlnFjxU66d01/YMTaLTvo0bVjwQwySXzRac4XNHdn0ierGTlsP/Zqr+/PbZ0CVIZU1UT/hamug2SrBZWKhzdr6D2berqvk725A/fEV+95b4/3fWXO6t3SElM5TVm8nkmfRNvfXrCGtxesqZPv4fcWMznkXbFx+27HcXfMjP9+ZhZvzo/uG/vbB0v52wdL6da5A4N7783M5Rtr8196/CCqqp1h/boxfFBPPly0jttemFPnmK/NXc3Iw/ar8xiVl2evYkifvSlesp7rn5rJhccO4I6vHkXp+gpWb9pB9y57cVDvrmzctpOqGqd/jy61+9bUOI9NWcZXj+3Pmi2VdNmrPb323r2Fn6pupeu3MbBX1ybzNiZRjeb8/b85v5wrH5nK1acdyk9GH9ai98+VtVt2sHRdBccM6pnvorQ6ClAZUl0TBaaq6t3/G72R/9CqFlzYTxVICqt9lEIGK9DUc7CmLV3f4LamPlSrapyhN76Qctvm7VV1ghPAo5OXNn5A4DsPF/PHS47hnKMP5OPSDdz58nzeCsEvYfzUUnp23Yu/vL2oNi3RwgP41QWf45LjBjFt6Xr+a/zHLCjbwg1PRzdYJ24DqK5xfvjoVCbO3hXArxk1lO+dcjBdO3bg6enL+fETMxjSZ2+OPagnN5z1WYbfMon7LjuW0484AICSsi0c2KMzXTum/ghJ/t0+XryMm845fLdRmD+fMJuH31vM0z8cwbbKak48pDdrtkQP4Vy+fhtL11awf/dOdGzfjqXrKjio9961+9bUOC/MWslZR/ajXb0W5YaKSrp0bE+nDqlbro0p27ydPnt32u2YjTn3T++yfMO2vIzCnLV8I1t3VHH8wb1327atsprtO6vpmcaXkoSa8EW7OfXPJgWoDEkEpsoUAefApG+19W3dkdmn4CZunG0oKDbVOmrqg7k6y83BpkYZNkdLusKa6gq9782FjW7fU1f/fToLyrYwcdYq5q1OfSvAi7NW1VlPbuFd/9RMjh/Siwv+vHvLM3z2sGLDtjrBCeAPry7g+ZkreeXHp/DRsg1AdO1t0ZqtjJ8aXev702slnH7EAbw5v5wxYz+ss//++3bir985nsMO6MZ3x03hlTllDDugW+32p6cv5xvHDaJ9O+OM37/FuZ8/kIfDtbyvhLKOGrYfZxwZBcAXZq3kqenRdcQRh/TmvYVrGf/9E1mzpZL/fmYWa7YkpuGazjGDevD0D0/i3jcXcnt4rMywA7pxwfD+rK/YybbKaiqra7j6tEP56/tLuOeNhfTv0YX27Yw3f3oqa7ZU8t7CNbjDjx7/iK8XDWTMiMGYwapN2zmkzz48Nb2Ug3p35dyj+9O+nfHa3NV03qs9Gyt2snzDNgDKN++gb7dObN6+kyVrKzi83760a2f8+qW5LF1XwZ++MXy3c7Jx205G/OpV7r+8iJMO7VNn287qGp7/eCXdu+7FyMP2221fd6+9zSNVcPzKn99l7qrNKbdVVddQWV2z2xeMEbe/xo6qaqb/v9G77ZPwwsyVHHZANw7pu0+DeTJFAaoZZq/YyKzlGzlqQA8+22/fOtsSXXypWlA9uzb8DWZL5Z4/7bYsaa68hO1hhNyezvDw3sI1jW7f3sQQ8pbKZABs10Q0bmxrU8X4dM3W5hcoTXe9uqDR7U2V7YJGukXXb62s/Vutr6QsGi3Z0O8t8eWhfnACWL1pB/87aR73XVbEK3OiwRHJ99rd9Mwstmyv4vunHMzcVZuZ+9K83Y7x6tyy2g/7xOTFAO+FQStzVm3mpnojPAGmL40CaiI4Jd77ly/UfbTM/yW1YhPvM3XJem557hNmlO5q8T5evIzHi5el/B3c+PQsLjp2AOPe333gxxdue4XvnXIw/5i6PCmA7lK8+FVWbdrOBcP789S05Vx92qH88bUSAC59YDL/Ofoz/HbSfJ76wQiOGdST256fUxvEn7nqJM6/+106tDMm3zCKh95dXGc07ODrngega8f2VFRW88NTD6n9/X/t3vc548gDGHFob+6cOI+d1V7bLQ0w+39O56ZnZlHtzqpN0ZedR95fzND9u3H7i3P56rEDeGnWSi44ZgCvzS3j+ZkrAZhy45eYv3ozRw/skbV7FK2x7qdCYGZnAH8A2gMPuPvtjeUvKiry4uLiPXqvoltfqf3D++e/n0yPrnvRv0cX3phfxnce3nXMOy86mkVrtnDJcYMAGPvOYsa+u6tL5tNfnkW7dkZVdQ3Tlm7ga/e9X+d9Fv3qLKYt3UCvvTsy8s43divHTWcfzi3PNTyU+/3rT2Ptlso6N9Em7NetEzedfThX/316s+qecNyQXll9JMdn++3LnJWbMnKsX37lc7XdW6l89+QhPPDOopTbBvTsQun6bQ3uW3RQT4qXNNxFGGdXjTyEu19P3QL85gmD2LStqtHRlI054sB9mb0iM+dPCkdLuzfNbKq7F+2WXsgByszaA/OBfwVKgSnAJe7e4Kd3SwJU4luKiIjskq0AVejjOI8DStz9U3evBB4DzstzmURE2pTtjdx83xKFHqD6A8mdxaUhrQ4zu9LMis2suLy8vP7mjHnrpyM5qHfDQ3NPrncRNNnMn4/m+6cc0uD2fxnap84Q4lwb3Ei9RKRt29LEiNk9VehdfBcBp7v7d8P6ZcBx7n51Q/u0pItPREQyr7V28ZUCA5PWBwB7dnVXRERipdAD1BRgqJkNMbOOwMXAhDyXSUREMqCg74Ny9yoz+3dgItEw87HuPruJ3UREpAAUdIACcPcXgNRzzoiISMEq9C4+ERFppRSgREQklhSgREQklhSgREQklhSgREQklgp6Jok9YWblwO5z5aevD9D4MylaD9W1dVJdW6dCrutB7t63fmKbC1AtZWbFqabkaI1U19ZJdW2dWmNd1cUnIiKxpAAlIiKxpADVfPfnuwA5pLq2Tqpr69Tq6qprUCIiEktqQYmISCwpQImISCwpQKXJzM4ws3lmVmJm1+W7PHvCzAaa2etmNsfMZpvZNSG9l5lNMrMF4WfPkG5mdleo88dmNjzpWGNC/gVmNiZfdWqKmbU3s+lm9lxYH2Jmk0O5Hw/PEcPMOoX1krB9cNIxrg/p88zs9PzUpHFm1sPMxpvZ3HB+T2yt59XMrg1/v7PM7O9m1rk1nVczG2tmZWY2KyktY+fSzI41s5lhn7vMzHJbw2Zwd72aeBE9a2ohcDDQEZgBHJ7vcu1BPfoBw8NyN2A+cDjwa+C6kH4dcEdYPgt4ETDgBGBySO8FfBp+9gzLPfNdvwbq/GPg/4DnwvoTwMVh+V7gB2H5h8C9Yfli4PGwfHg4352AIeHvoH2+65WinuOA74bljkCP1nhegf7AIqBL0vn8Vms6r8AXgeHArKS0jJ1L4EPgxLDPi8CZ+a5zg7+LfBegEF7hZE5MWr8euD7f5cpAvZ4F/hWYB/QLaf2AeWH5PuCSpPzzwvZLgPuS0uvki8sLGAC8CpwGPBf+IdcAHeqfV6KHXp4YljuEfFb/XCfni8sL2Dd8aFu99FZ3XkOAWhY+eDuE83p6awgJ8c8AAAYDSURBVDuvwOB6ASoj5zJsm5uUXidf3F7q4ktP4p8ioTSkFazQ1XEMMBnY391XAoSf+4VsDdW7UH4fvwd+BtSE9d7ABnevCuvJ5a6tU9i+MeQvhLoeDJQDD4XuzAfMbG9a4Xl19+XAncBSYCXReZpK6zyvyTJ1LvuH5frpsaQAlZ5UfbQFOz7fzPYB/gH8yN03NZY1RZo3kh4bZnY2UObuU5OTU2T1JrbFvq5ELYPhwD3ufgywlagbqCEFW9dw7eU8om65A4G9gTNTZG0N5zUdza1fQdVbASo9pcDApPUBwIo8laVFzGwvouD0qLs/FZJXm1m/sL0fUBbSG6p3Ifw+TgLONbPFwGNE3Xy/B3qYWYeQJ7nctXUK27sD6yiMupYCpe4+OayPJwpYrfG8fglY5O7l7r4TeAoYQes8r8kydS5Lw3L99FhSgErPFGBoGCnUkehi64Q8l6nZwmidB4E57v6/SZsmAIlRPmOIrk0l0i8PI4VOADaG7oWJwGgz6xm+0Y4OabHh7te7+wB3H0x0vl5z90uB14ELQ7b6dU38Di4M+T2kXxxGgw0BhhJdZI4Nd18FLDOzw0LSKOATWuF5JeraO8HMuoa/50RdW915rScj5zJs22xmJ4Tf3+VJx4qffF8EK5QX0WiZ+USjfW7Md3n2sA4nEzXnPwY+Cq+ziPrkXwUWhJ+9Qn4D7g51ngkUJR3rO0BJeH0733Vrot6nsmsU38FEH0QlwJNAp5DeOayXhO0HJ+1/Y/gdzCOmI56AzwPF4dw+QzRyq1WeV+B/gLnALOARopF4rea8An8nur62k6jFc0UmzyVQFH53C4E/UW9wTZxemupIRERiSV18IiISSwpQIiISSwpQIiISSwpQIiISSwpQIiISSwpQImkws2oz+yjpNTjfZcokMzvGzB4Iy98ysz/V2/6GmRU1sv9jZjY02+WUtqVD01lEBNjm7p9vaKOZdfBdc8EVohuAW1uw/z1E8x7+W2aKI6IWlMgeCy2NJ83sn8DLIe2nZjYlPJvnf5Ly3hieO/RKeIbRf4b02paJmfUJUzMlnmP1m6RjfS+knxr2STz76dHE83zM7Atm9p6ZzTCzD82sm5m9bWafTyrHu2Z2VL16dAOOcvcZadT53KRW5DwzWxQ2vQ18KWm6IZEW0x+TSHq6mNlHYXmRu38lLJ9I9OG+zsxGE02ZcxzRHf4TzOyLRJO3Xkw0e3wHYBrRDNyNuYJo2povmFkn4F0zezlsOwY4gmgOtXeBk8zsQ+Bx4OvuPsXM9gW2AQ8QPS/pR2b2GaIZFj6u916JmQWSfd3MTk5aPxTA3ScQpvkysyeAN0N6jZmVAEenUTeRtChAiaSnoS6+Se6+LiyPDq/pYX0fooDVDXja3SsAzCydeRxHA0eZWWJ+ue7hWJXAh+5eGo71EdGzgzYCK919CoCHWerN7EngJjP7KdHUNw+neK9+RI/rSPa4u/97YsXM3kjeaGY/I/qd3J2UXEY0w7gClGSEApRIy2xNWjbgV+5+X3IGM/sRDT/SoIpdXe2d6x3ranevM1mrmZ0K7EhKqib6P7ZU7+HuFWY2iegRFV8jai3Vt63eezfKzEYBFxE9+TVZ53AskYzQNSiRzJkIfMei521hZv3NbD/gLeArZtYlXO85J2mfxcCxYfnCesf6gUWPR8HMPmPRQwgbMhc40My+EPJ3S7oe9ABwFzAlqbWXbA6hC68pZnYQ8Gfga+5ePxh9BpidznFE0qEWlEiGuPvLZvZZ4P0wbmEL8E13n2ZmjxPNHr+EaEBBwp3AE2Z2GfBaUvoDRF1308IgiHLg/Ebeu9LMvg780cy6ELVkvgRscfepZrYJeKiBfeeaWXcz6+bum5uo5reIZtZ+OtRxhbufZWb7E3X5rWxif5G0aTZzkRwzs58TBY47c/R+BwJvAMPcvaaBPNcCm939gT18j2uBTe7+4B4XVKQedfGJtGJmdjkwmegZZimDU3APda9tNdcGYFwL9hfZjVpQIiISS2pBiYhILClAiYhILClAiYhILClAiYhILClAiYhILP1/qgw6Doxx0oUAAAAASUVORK5CYII=\n",
"text/plain": [
"