{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## ThinkDSP\n",
"\n",
"This notebook contains code examples from Chapter 6: Discrete Cosine Transform\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",
"PI2 = np.pi * 2"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Synthesis\n",
"\n",
"The simplest way to synthesize a mixture of sinusoids is to add up sinusoid signals and evaluate the sum."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"from thinkdsp import CosSignal, SumSignal\n",
"\n",
"def synthesize1(amps, fs, ts):\n",
" components = [CosSignal(freq, amp)\n",
" for amp, freq in zip(amps, fs)]\n",
" signal = SumSignal(*components)\n",
"\n",
" ys = signal.evaluate(ts)\n",
" return ys"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Here's an example that's a mixture of 4 components."
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
" \n",
" "
],
"text/plain": [
""
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"from thinkdsp import Wave\n",
"\n",
"amps = np.array([0.6, 0.25, 0.1, 0.05])\n",
"fs = [100, 200, 300, 400]\n",
"framerate = 11025\n",
"\n",
"ts = np.linspace(0, 1, framerate, endpoint=False)\n",
"ys = synthesize1(amps, fs, ts)\n",
"wave = Wave(ys, ts, framerate)\n",
"wave.apodize()\n",
"wave.make_audio()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can express the same process using matrix multiplication."
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"def synthesize2(amps, fs, ts):\n",
" args = np.outer(ts, fs)\n",
" M = np.cos(PI2 * args)\n",
" ys = np.dot(M, amps)\n",
" return ys"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"And it should sound the same."
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
" \n",
" "
],
"text/plain": [
""
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"ys = synthesize2(amps, fs, ts)\n",
"wave = Wave(ys, framerate)\n",
"wave.apodize()\n",
"wave.make_audio()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"And we can confirm that the differences are small."
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"1.2789769243681803e-13"
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"ys1 = synthesize1(amps, fs, ts)\n",
"ys2 = synthesize2(amps, fs, ts)\n",
"np.max(np.abs(ys1 - ys2))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Analysis\n",
"\n",
"The simplest way to analyze a signal---that is, find the amplitude for each component---is to create the same matrix we used for synthesis and then solve the system of linear equations."
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [],
"source": [
"def analyze1(ys, fs, ts):\n",
" args = np.outer(ts, fs)\n",
" M = np.cos(PI2 * args)\n",
" amps = np.linalg.solve(M, ys)\n",
" return amps"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Using the first 4 values from the wave array, we can recover the amplitudes."
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([0.6 , 0.25, 0.1 , 0.05])"
]
},
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"n = len(fs)\n",
"amps2 = analyze1(ys[:n], fs, ts[:n])\n",
"amps2"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"What we have so far is a simple version of a discrete cosine tranform (DCT), but it is not an efficient implementation because the matrix we get is not orthogonal."
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [],
"source": [
"# suppress scientific notation for small numbers\n",
"np.set_printoptions(precision=3, suppress=True)"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([[ 1. , 1. , 1. , 1. ],\n",
" [ 1. , 0.707, 0. , -0.707],\n",
" [ 1. , 0. , -1. , -0. ],\n",
" [ 1. , -0.707, -0. , 0.707]])"
]
},
"execution_count": 12,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"def test1():\n",
" N = 4.0\n",
" time_unit = 0.001\n",
" ts = np.arange(N) / N * time_unit\n",
" max_freq = N / time_unit / 2\n",
" fs = np.arange(N) / N * max_freq\n",
" args = np.outer(ts, fs)\n",
" M = np.cos(PI2 * args)\n",
" return M\n",
"\n",
"M = test1()\n",
"M"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"To check whether a matrix is orthogonal, we can compute $M^T M$, which should be the identity matrix:"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([[ 4., 1., -0., 1.],\n",
" [ 1., 2., 1., 0.],\n",
" [-0., 1., 2., 1.],\n",
" [ 1., 0., 1., 2.]])"
]
},
"execution_count": 13,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"M.transpose().dot(M)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"But it's not, which means that this choice of M is not orthogonal.\n",
"\n",
"Solving a linear system with a general matrix (that is, one that does not have nice properties like orthogonality) takes time proportional to $N^3$. With an orthogonal matrix, we can get that down to $N^2$. Here's how:"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([[ 0.981, 0.831, 0.556, 0.195],\n",
" [ 0.831, -0.195, -0.981, -0.556],\n",
" [ 0.556, -0.981, 0.195, 0.831],\n",
" [ 0.195, -0.556, 0.831, -0.981]])"
]
},
"execution_count": 14,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"def test2():\n",
" N = 4.0\n",
" ts = (0.5 + np.arange(N)) / N\n",
" fs = (0.5 + np.arange(N)) / 2\n",
" args = np.outer(ts, fs)\n",
" M = np.cos(PI2 * args)\n",
" return M\n",
" \n",
"M = test2()\n",
"M"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now $M^T M$ is $2I$ (approximately), so M is orthogonal except for a factor of two."
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([[ 2., -0., 0., 0.],\n",
" [-0., 2., -0., -0.],\n",
" [ 0., -0., 2., -0.],\n",
" [ 0., -0., -0., 2.]])"
]
},
"execution_count": 15,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"M.transpose().dot(M)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"And that means we can solve the analysis problem using matrix multiplication."
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [],
"source": [
"def analyze2(ys, fs, ts):\n",
" args = np.outer(ts, fs)\n",
" M = np.cos(PI2 * args)\n",
" amps = M.dot(ys) / 2\n",
" return amps"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"It works:"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([0.6 , 0.25, 0.1 , 0.05])"
]
},
"execution_count": 17,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"n = len(fs)\n",
"amps2 = analyze1(ys[:n], fs, ts[:n])\n",
"amps2"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### DCT\n",
"\n",
"What we've implemented is DCT-IV, which is one of several versions of DCT using orthogonal matrices."
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {},
"outputs": [],
"source": [
"def dct_iv(ys):\n",
" N = len(ys)\n",
" ts = (0.5 + np.arange(N)) / N\n",
" fs = (0.5 + np.arange(N)) / 2\n",
" args = np.outer(ts, fs)\n",
" M = np.cos(PI2 * args)\n",
" amps = np.dot(M, ys) / 2\n",
" return amps"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can check that it works:"
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"8.326672684688674e-17"
]
},
"execution_count": 19,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"amps = np.array([0.6, 0.25, 0.1, 0.05])\n",
"N = 4.0\n",
"ts = (0.5 + np.arange(N)) / N\n",
"fs = (0.5 + np.arange(N)) / 2\n",
"ys = synthesize2(amps, fs, ts)\n",
"\n",
"amps2 = dct_iv(ys)\n",
"np.max(np.abs(amps - amps2))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"DCT and inverse DCT are the same thing except for a factor of 2."
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {},
"outputs": [],
"source": [
"def inverse_dct_iv(amps):\n",
" return dct_iv(amps) * 2"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"And it works:"
]
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"8.326672684688674e-17"
]
},
"execution_count": 21,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"amps = [0.6, 0.25, 0.1, 0.05]\n",
"ys = inverse_dct_iv(amps)\n",
"amps2 = dct_iv(ys)\n",
"np.max(np.abs(amps - amps2))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Dct\n",
"\n",
"`thinkdsp` provides a `Dct` class that encapsulates the DCT in the same way the Spectrum class encapsulates the FFT."
]
},
{
"cell_type": "code",
"execution_count": 22,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
" \n",
" "
],
"text/plain": [
""
]
},
"execution_count": 22,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"from thinkdsp import TriangleSignal\n",
"\n",
"signal = TriangleSignal(freq=400)\n",
"wave = signal.make_wave(duration=1.0, framerate=10000)\n",
"wave.make_audio()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"To make a Dct object, you can invoke `make_dct` on a Wave."
]
},
{
"cell_type": "code",
"execution_count": 23,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAAAadElEQVR4nO3de7hddZ3f8feHkxCC3DUgJmjiGAeBqkhEFLU+g5WoM+K0XtBHQWUGa9WqrfqAtl6ep4x2aqdT6kiHwRGsVERHK/p4ozjqlEEggCjILcgtkiHHCxAuBpJ8+8dah2zONSRnn73OOe/Xk/3stX/r9lu/c7I/Z/32b6+VqkKSpK7ZZdAVkCRpPAaUJKmTDChJUicZUJKkTjKgJEmdtGDQFeiXJzzhCbV8+fJBV0OSNIUrrrjiV1W1ZHT5nA2o5cuXs2bNmkFXQ5I0hSS3jVduF58kqZMMKElSJxlQkqROMqAkSZ1kQEmSOsmAkiR1kgElSeokA0qS1El9Dagk70tybZJrknwxyW5J9ktyYZKb2ud9e5Y/NcnaJDckOban/IgkP2vnnZ4k/az3TPrrH97M7b9+YNDVkKTO6VtAJVkK/FtgVVUdBgwBxwOnABdV1UrgovY1SQ5p5x8KrAY+k2So3dwZwMnAyvaxul/1nkm/vf8hPvHt63njWT8edFUkqXP63cW3AFicZAGwO3AncBxwTjv/HODV7fRxwHlVtamqbgHWAkcmORDYq6ouqeb2v5/vWWdW29rezfiBh7YMuCaS1D19C6iq+iXwKeB2YD1wT1V9Dzigqta3y6wH9m9XWQrc0bOJdW3Z0nZ6dPkYSU5OsibJmuHh4ek8HEnSDOtnF9++NGdFK4AnAY9L8qbJVhmnrCYpH1tYdWZVraqqVUuWjLkwriRpFulnF99LgVuqariqHga+CrwAuKvttqN93tAuvw44qGf9ZTRdguva6dHlkqQ5rJ8BdTtwVJLd21F3xwDXARcAJ7bLnAh8vZ2+ADg+yaIkK2gGQ1zWdgNuTHJUu50TetaRJM1RfbsfVFVdmuQrwJXAZuAq4ExgD+D8JCfRhNhr2+WvTXI+8PN2+XdW1cjogXcAZwOLgW+3D0nSHNbXGxZW1UeBj44q3kRzNjXe8qcBp41TvgY4bNorKEnqLK8kIUnqJANKktRJBpQkqZMMKElSJxlQkqROMqAkSZ1kQEmSOsmAkiR1kgElSeokA0qS1EkGlCSpkwwoSVInGVADNO5dFyVJgAHVCePdMliS5jsDSpLUSQaUJKmTDChJUicZUJKkTjKgOsDRfJI0lgE1QI7ek6SJGVCSpE4yoCRJnWRASZI6yYCSJHWSASVJ6iQDSpLUSQaUJKmTDChJUicZUJKkTjKgJEmdZEBJkjrJgJIkdZIBNUBexVySJmZAdYBXNZeksQwoSVInGVCSpE4yoCRJnWRASZI6yYCSJHWSASVJ6qS+BlSSfZJ8Jcn1Sa5L8vwk+yW5MMlN7fO+PcufmmRtkhuSHNtTfkSSn7XzTk/iyGxJmuP6fQb134HvVNXBwLOA64BTgIuqaiVwUfuaJIcAxwOHAquBzyQZardzBnAysLJ9rO5zvSVJA9a3gEqyF/Bi4LMAVfVQVd0NHAec0y52DvDqdvo44Lyq2lRVtwBrgSOTHAjsVVWXVFUBn+9ZR5I0R/XzDOqpwDDwuSRXJTkryeOAA6pqPUD7vH+7/FLgjp7117VlS9vp0eVjJDk5yZoka4aHh6f3aCRJM6qfAbUAeA5wRlUdDtxP2503gfE+V6pJyscWVp1ZVauqatWSJUsea30lSR3Sz4BaB6yrqkvb11+hCay72m472ucNPcsf1LP+MuDOtnzZOOVzhheNlaSx+hZQVfVPwB1Jfr8tOgb4OXABcGJbdiLw9Xb6AuD4JIuSrKAZDHFZ2w24MclR7ei9E3rWmdUciihJE1vQ5+2/Gzg3ya7AL4C30oTi+UlOAm4HXgtQVdcmOZ8mxDYD76yqLe123gGcDSwGvt0+JElzWF8Dqqp+AqwaZ9YxEyx/GnDaOOVrgMOmtXKSpE7zShKSpE4yoCRJnWRASZI6yYAaIIeXS9LEDKgOcLi5JI1lQEmSOsmAkiR1kgElSeokA0qS1EkGlCSpkwwoSVInGVCSpE4yoCRJnWRASZI6yYCSJHWSASVJ6iQDSpLUSQaUJKmTDChJUicZUJKkTjKgJEmdZEBJkjrJgJIkdZIB1QE16ApIUgcZUJKkTjKgOiCDroAkdZABJUnqJANKktRJBpQkqZMMKElSJxlQkqROMqAkSZ1kQEmSOsmAkiR1kgElSeokA0qS1EkGlCSpkwwoSVInGVCSpE4yoCRJndT3gEoylOSqJN9sX++X5MIkN7XP+/Yse2qStUluSHJsT/kRSX7Wzjs9iXeokKQ5bibOoN4DXNfz+hTgoqpaCVzUvibJIcDxwKHAauAzSYbadc4ATgZWto/VM1BvSdIATRpQSf5sZzaeZBnwSuCsnuLjgHPa6XOAV/eUn1dVm6rqFmAtcGSSA4G9quqSqirg8z3rSJLmqKnOoHb2TOUvgQ8CW3vKDqiq9QDt8/5t+VLgjp7l1rVlS9vp0eVjJDk5yZoka4aHh3ey6pKkQZoqoIaS7Nt+bjTmMdmKSf4Q2FBVV2xnXcb7XKkmKR9bWHVmVa2qqlVLlizZzt0OTo17FJIkgAVTzD8YuIKJQ+Kpk6x7NPCqJK8AdgP2SvIF4K4kB1bV+rb7bkO7/DrgoJ71lwF3tuXLximfMxzyIUljTXUG9fOqempVrRjnMVk4UVWnVtWyqlpOM/jh+1X1JuAC4MR2sROBr7fTFwDHJ1mUZAXNYIjL2m7AjUmOakfvndCzjiRpjprqDKofPgmcn+Qk4HbgtQBVdW2S84GfA5uBd1bVlnaddwBnA4uBb7cPSdIcNlVA/U2SJVX1qBEHSfYH7q2q323PTqrqB8AP2ulfA8dMsNxpwGnjlK8BDtuefUmS5oapuvieDbxonPJ/Afy3aa+NJEmtqQLqhVX11dGFVXUu8OL+VGn+cTSfJI01VUBNNr7M6/jtJEfvSdLEpgqZDUmOHF2Y5LmA34SVJPXNVIMkPkAz4u5smu9DAayiGep9fB/rJUma5yY9g6qqy4Dn0XT1vaV9BHheVV3a78pJkuavKb8HVVV3AR9NsqR9bdeeJKnvprqaeZJ8LMkwcD1wQ5LhJB+ZmepJkuarqQZJvJfmmnpHVtXjq2o/mi6/o5O8r9+VkyTNX1MF1AnAG9r7MwFQVb8A3tTOkySpL6YKqIVV9avRhe3nUAv7UyVJkqYOqId2cJ4kSTtlqlF8z0py7zjlobnHkyRJfTFpQFXV0ExVRJKkXl5PT5LUSQbUAHkVc0mamAHVAV7VXJLGMqAkSZ1kQEmSOsmAkiR1kgElSeokA0qS1EkGlCSpkwwoSVInGVCSpE4yoCRJnWRASZI6yYDqAK/JJ0ljGVAD5DX4JGliBpQkqZMMKElSJxlQkqROMqAkSZ1kQEmSOsmAkiR1kgE1QH7/SZImZkB1gN+HkqSxDChJUicZUJKkTjKgJEmd1LeASnJQkr9Pcl2Sa5O8py3fL8mFSW5qn/ftWefUJGuT3JDk2J7yI5L8rJ13euKnNpI01/XzDGoz8O+r6hnAUcA7kxwCnAJcVFUrgYva17TzjgcOBVYDn0ky1G7rDOBkYGX7WN3HekuSOqBvAVVV66vqynZ6I3AdsBQ4DjinXewc4NXt9HHAeVW1qapuAdYCRyY5ENirqi6pqgI+37OOJGmOmpHPoJIsBw4HLgUOqKr10IQYsH+72FLgjp7V1rVlS9vp0eXj7efkJGuSrBkeHp7WY5Akzay+B1SSPYC/A95bVfdOtug4ZTVJ+djCqjOralVVrVqyZMljr6wkqTP6GlBJFtKE07lV9dW2+K622472eUNbvg44qGf1ZcCdbfmyccolSXNYP0fxBfgscF1V/UXPrAuAE9vpE4Gv95Qfn2RRkhU0gyEua7sBNyY5qt3mCT3rSJLmqAV93PbRwJuBnyX5SVv2IeCTwPlJTgJuB14LUFXXJjkf+DnNCMB3VtWWdr13AGcDi4Fvtw9J0hzWt4Cqqv/H+J8fARwzwTqnAaeNU74GOGz6aidJ6jqvJCFJ6iQDSpLUSQaUJKmTDChJUicZUJKkTjKgOsBbv0vSWAaUJKmTDChJUicZUB3g7RclaSwDSpLUSQaUJKmTDChJUicZUJKkTjKgJEmdZEBJkjrJgJIkdZIBJUnqJANKktRJBpQkqZMMKElSJxlQkqROMqA07S65+df87uEtg66GpFnOgNK0unn4Pt7wNz/mP/yfawZdFUmznAGlaXXvgw8DcNOG+wZcE0mznQGlaeXd6yVNFwNKktRJBpQkqZMMqAGqOdwh5l3sJe0sA2qAqs2nxLdzSRrNgBqgkYDaxXySpDEMqAHa2iZU5lCHWM3dXktJM8yAGqCR93LPoCRpLANqgLZubc+g/AxKksYwoCRJnWRADdAjgyT8KUjSGL41DtBcHCQx8snaXOu1XH/Pg5QjQKQZZUANkIMkZodrfnkPz//E9/nCj28bdFWkecWAGqBHzqDm2unGHHPzcHNl9ktv+c2AayLNLwbUAG27ksRg6zGdtmxtnofm0kG1dpmDxyR1mQE1QBt/19w76bf3PzTgmkyfzVubhFowNHfezEfOdOdaV+xP193t52rqtFkTUElWJ7khydokpwy6PtPh3EtvB+C3Dzw84JpMn81bmje8BXNoaGKbuXPqDOo716znVZ++mK9e+ctBV2VaPfjQlkFXQdNoVryLJBkC/gp4OXAI8IYkhwy2VjvvK1esm7F93fPAw/xiuP93ub2nvaPuTFyp/Z4HHubOux/s+34efLh509u8tf/HdN+mzY+0YT/deFfzu3Drr+/v+742b9n6yJfS++kbV9/JMz7yHW68a2Pf97X+ngdn5Of0i+H7eP+Xr2bzSN95H/3oxmHu+M0Dfd/PY5HZcIqf5PnAx6rq2Pb1qQBV9YmJ1lm1alWtWbNmh/b3wxuH+c41//TIL8VIC22t4oFNWx55o6pq3oa3VnHTXffxq/s2sWnz9P0iHbliP55+wB4s2GUXkm3D0Uf/Id/7IxwvGDZvKf5XOwLtrUcvHzN/e38Fen9XtlRRBVureQO6ccN9XH3H3WPWeeayvXnakj3Ye/eFhGx33UOmDLnPXXwrAG87esWky+7o8W0tuO3X93Px2l8/Un7Mwfvz5Mfv/qhj2ZljGG17j2ls3bdnmeaYHtq8lS+tueOR8n/5nKXss3hXkrHbKWqHjqO3Tmf/460AnPTCFY90l+6IiVatKrZU8YUf3/5I2duOXvGo4xk5jpHpndl/VXHOJc3/p5NeuIKq/v1BNvL7kMBbX7BizM9jR382o1Vt+zm99ejlj2m7T9xrN97+z39vp/af5IqqWjWmfJYE1GuA1VX1J+3rNwPPq6p3jVruZOBkgCc/+clH3Hbbjg0LPvNHN/Nfv3cjC4eaE8yRNtpa2/6anklDu6QZdND8Y6j9MGTkvb53FOCj3v+zbf27227EPXdbMO4+QhPEI+v3To+n97dm69bi/im6VnYdakJ2e+v+qP1PUJGNv9sMwB6LFjSLTFLhkeMb/XrSYyzYuGnzmPJdh3Zhl12aLr/t+hkU27dDxjmmMRuc2HjHOJ4t4/y8Fi8cav8I2s5jge06ngD3tsf0uF2HSE+b7czX/0avWmxrO2h+RguHRv6gy8RVfgzH0rv4pMe0ndvaHr3HtPuuQ490Mz+q2o/xd2yqfY3sZ3t/9gc/cU++/K9fsGM7Hdn8BAE1/rtV94zX5GOStarOBM6E5gxqR3d24guW86cvemrfh39v2Vr83oe+BcDzVuzHaX/8z1i6z2IWDmXbG/k01qGq+n5Mt/7qfl7yqR8AsPa0lzO0S/q2z02bt/C7h7ey9+KFfdn+iCtu+w3/6oxLSOAXf/aKvrbh8MZNPPjQFp78+N37tg9oupff/+WrWX3oE/mfbz6ir/u65pf3sGVr8ayD9unrfj538S18/Bs/5+OvOpQTX7C8r/v6xtV38rhFQ/zBwQf0fT/v/uJVfOe9L+LgJ+7V13198CtX87T99+DkF+/c2dB0mi0BtQ44qOf1MuDOfu1s0YKhfm16Ql96+/P7vo+Z+L5V70CCBUP9/Yhz0YKhGflZbXq46bZ93or9+t6GS/Zc1Nftj3h2GxZvfN6T+76vw5bu3fd9ALzpqKew28IhXrfqoKkX3kl/9Kwn9X0fI/s55hn7s/uu/X+r/vPXPKvv+3isZktAXQ6sTLIC+CVwPPDGwVZJ4xnptz5ov8UDrsn0WXnAngC85QUrBlyT6fO0/ffg1k++ctDVmFYLh3bhDUf2P3Bn2kyEU1fNiiOvqs1J3gV8FxgC/raqrh1wtTSObXcJnjtDspfsuWjOvZlLs8GsCCiAqvoW8K1B10OT22f35vOglx924IBrImm2mzUBpdlhn9135eqPvGzC0YKStL18F9G023v3/o6qkzQ/zIorSUiS5h8DqgPm0HgCSZo2BpQkqZMMqA7wBEqSxjKgBmg2XAdRkgbFgOoAb/kuSWMZUAPk+ZMkTcyA6gDPnyRpLL+oO0ALdglHP+3xvO3ouXMRUkmaLgbUACXh3D85atDVkKROsotPktRJBpQkqZMMKElSJxlQkqROMqAkSZ1kQEmSOsmAkiR1kgElSeokA0qS1EmZq7d8SDIM3LYTm3gC8Ktpqs5sZjtsY1s0bIeG7dCYjnZ4SlUtGV04ZwNqZyVZU1WrBl2PQbMdtrEtGrZDw3Zo9LMd7OKTJHWSASVJ6iQDamJnDroCHWE7bGNbNGyHhu3Q6Fs7+BmUJKmTPIOSJHWSASVJ6iQDahxJVie5IcnaJKcMuj7TLcnfJtmQ5Jqesv2SXJjkpvZ53555p7ZtcUOSY3vKj0jys3be6Uky08eyM5IclOTvk1yX5Nok72nL51VbJNktyWVJrm7b4eNt+bxqhxFJhpJcleSb7et51w5Jbm3r/5Mka9qymW+HqvLR8wCGgJuBpwK7AlcDhwy6XtN8jC8GngNc01P258Ap7fQpwH9upw9p22ARsKJtm6F23mXA84EA3wZePuhje4ztcCDwnHZ6T+DG9njnVVu0dd6jnV4IXAocNd/aoac9/h3wv4Fvtq/nXTsAtwJPGFU24+3gGdRYRwJrq+oXVfUQcB5w3IDrNK2q6kfAb0YVHwec006fA7y6p/y8qtpUVbcAa4EjkxwI7FVVl1Tzm/j5nnVmhapaX1VXttMbgeuApcyztqjGfe3Lhe2jmGftAJBkGfBK4Kye4nnXDhOY8XYwoMZaCtzR83pdWzbXHVBV66F54wb2b8snao+l7fTo8lkpyXLgcJqzh3nXFm231k+ADcCFVTUv2wH4S+CDwNaesvnYDgV8L8kVSU5uy2a8HRbsQMXnuvH6SOfzWPyJ2mPOtFOSPYC/A95bVfdO0k0+Z9uiqrYAz06yD/C1JIdNsvicbIckfwhsqKorkrxke1YZp2zWt0Pr6Kq6M8n+wIVJrp9k2b61g2dQY60DDup5vQy4c0B1mUl3tafktM8b2vKJ2mNdOz26fFZJspAmnM6tqq+2xfOyLQCq6m7gB8Bq5l87HA28KsmtNF37f5DkC8y/dqCq7myfNwBfo/noY8bbwYAa63JgZZIVSXYFjgcuGHCdZsIFwInt9InA13vKj0+yKMkKYCVwWXuKvzHJUe3InBN61pkV2np/Friuqv6iZ9a8aoskS9ozJ5IsBl4KXM88a4eqOrWqllXVcpr/99+vqjcxz9ohyeOS7DkyDbwMuIZBtMOgR4t08QG8gmZE183Ahwddnz4c3xeB9cDDNH/lnAQ8HrgIuKl93q9n+Q+3bXEDPaNwgFXtL+7NwKdpr0wyWx7AC2m6HH4K/KR9vGK+tQXwTOCqth2uAT7Sls+rdhjVJi9h2yi+edUONCOYr24f1468Bw6iHbzUkSSpk+zikyR1kgElSeokA0qS1EkGlCSpkwwoSVInGVDSBJJsaa/mPPJYPug6TZckhyc5q51+S5JPj5r/gySrJln/vCQr+11PzW9e6kia2INV9ezxZrRfPExVbR1v/izwIeA/7cT6Z9Bcs+5Pp6c60lieQUnbKcnyNPeO+gxwJXBQkg8kuTzJT9PeR6ld9sPtvXH+b5IvJnl/W/7ImUmSJ7SX1Rm5WOt/6dnW29vyl7TrfCXJ9UnOHbmnTpLnJvnHNPdxuizJnkn+Icmze+pxcZJnjjqOPYFnVtXV23HMr+o5g7whyS3trH8AXprEP3LVN/5ySRNb3F7hG+AW4H3A7wNvrap/k+RlNJd1OZLmwpgXJHkxcD/NpXIOp/k/diVwxRT7Ogm4p6qem2QRcHGS77XzDgcOpbmO2cXA0UkuA74EvL6qLk+yF/AgzW0i3gK8N8nTgUVV9dNR+xr5dn+v1yd5Yc/rpwFU1QW0l/pKcj7ww7Z8a5K1wLO249ikHWJASRN7VBdf+xnUbVX147boZe3jqvb1HjSBtSfwtap6oF1ve67l+DLgmUle077eu93WQzTXNVvXbusnwHLgHmB9VV0OUFX3tvO/DPzHJB8A3gacPc6+DgSGR5V9qare1XOsP+idmeSDNO3xVz3FG4AnYUCpTwwo6bG5v2c6wCeq6q97F0jyXia+rcBmtnWt7zZqW++uqu+O2tZLgE09RVto/t9mvH1U1QNJLqS5idzraM6WRntw1L4nleQY4LU0d2LutVu7Lakv/AxK2nHfBd6W5n5SJFma5v45PwL+OMni9vOeP+pZ51bgiHb6NaO29Y40t/8gydPbK0lP5HrgSUme2y6/Z8/nQWcBpwOXV9XoOydDc+fgp23PASZ5CvAZ4HVVNTqMnk5zMVGpLzyDknZQVX0vyTOAS9pxC/cBb6qqK5N8iebq6LfRDCgY8Sng/CRvBr7fU34WTdfdle0giGEmuT12VT2U5PXA/2hvkfEgzW0y7qvmhnv3Ap+bYN3rk+ydZM9qbnU/mbfQXMX6a+0x3llVr0hyAE2X3/op1pd2mFczl/osycdoguNTM7S/J9HcdPDgiYbBJ3kfsLGqztrBfbwPuLeqPrvDFZWmYBefNIckOQG4lOYePpN9R+sMHv3Z1mN1N3DOTqwvTckzKElSJ3kGJUnqJANKktRJBpQkqZMMKElSJxlQkqRO+v+QkRv4Eq5yywAAAABJRU5ErkJggg==\n",
"text/plain": [
"