{
"cells": [
{
"cell_type": "markdown",
"id": "efbaabb5-2ac7-40b0-bca7-1e933bb53aa1",
"metadata": {},
"source": [
"# Resynthesizing sound via spectogram\n",
"\n",
"After we took a look how we can generate sounds in Python via numpy we want to take a look how we can re-synthesize a signal using numpy and some machine learning algorithms we have learned.\n",
"This allows us to yield interesting high-quality sounds without relying on neural networks that are exhausive on resources.\n",
"\n",
"Instead of writing audio analysis algorithms from scratch we can use the library [librosa](https://librosa.org/doc/latest/index.html) which is a library build with the help of numpy."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "1f8cf66f-8078-43ea-8e84-19d1a4768da0",
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"import pandas as pd\n",
"import matplotlib.pyplot as plt\n",
"import matplotlib as mpl\n",
"import librosa\n",
"import librosa.display\n",
"from IPython.display import display, Audio\n",
"\n",
"# remove warning in plotting via librosa\n",
"import warnings\n",
"import matplotlib.cbook\n",
"warnings.filterwarnings(\"ignore\",category=matplotlib.cbook.mplDeprecation)\n",
"\n",
"np.random.seed(42)\n",
"\n",
"mpl.rcParams['figure.figsize'] = (15, 8)"
]
},
{
"cell_type": "markdown",
"id": "095ea451-553f-4a35-a3d8-91561eced228",
"metadata": {},
"source": [
"## Loading a sound file"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "316d9898-53f9-4df8-8f30-43fa72590fa6",
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
" \n",
" "
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"path = 'trumpet.flac'\n",
"data, sr = librosa.load(path, sr=None, mono=True)\n",
"display(Audio(path))"
]
},
{
"cell_type": "markdown",
"id": "30a1c925-24bd-43e9-b042-ca8789683020",
"metadata": {},
"source": [
"We start by loading a recording of a trumpet.\n",
"We set the samplerate to `None` because otherwise librosa will not use the original samplerate of the file but will downsample the signal to $22'050~\\text{Hz}$, which is handier from a computational point of view as it contains less data, but not very pleasing from an aesthetic experience.\n",
"We will also force the signal to be mono, but this could be an interesting extension to the algorithm.\n",
"\n",
"As [`librosa.load`](https://librosa.org/doc/main/generated/librosa.load.html) returns two values we can unpack this list by declaring two values on the left side.\n",
"`data` will contain the data of our audio signal in a [PCM (Pulse Code Modulation)](https://en.wikipedia.org/wiki/Pulse-code_modulation) format, which gives us for each point in time its amplitude.\n",
"The resolution of the time domain is defined by the samplerate, which is stored in variable `sr`.\n",
"The resolution of the amplitude is given by its bitsize.\n",
"16 bits allow for $2^{16} = 65'536$ discerete values and 24 bits allow for $2^{24} = 16'777'216$ discrete values of resolution.\n",
"\n",
"For this notebook we will use [trumpetsolo1_hikmet.aiff by emirdemirel](https://freesound.org/people/emirdemirel/sounds/416025/) from [freesound.org](https://freesound.org/) (converted to 16 bit FLAC) but you are welcome to exchange the file with another one.\n",
"\n",
"We can take a look at `data` which is simply an array of the amplitude floats in time."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "b060f5c5-6a87-41c6-9e5d-4a8ed30367d0",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([ 0.0000000e+00, -3.0517578e-05, 6.1035156e-05, ...,\n",
" 4.5776367e-04, 1.2207031e-04, -1.5258789e-04], dtype=float32)"
]
},
"execution_count": null,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"data"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "bc707256-1316-4472-a488-26640c61e124",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(896841,)"
]
},
"execution_count": null,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"data.shape"
]
},
{
"cell_type": "markdown",
"id": "31c45a8c-604b-409a-8211-9a4711719e51",
"metadata": {},
"source": [
"`data` is a single vector in time and we can verify the length of the signal by dividing it with our samplerate `sr`."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "7efe3ad6-c3d3-4c2b-9f9c-00755ccc8134",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Samplerate is 44100\n",
"Length of the signal in seconds: 20.336530612244896\n"
]
}
],
"source": [
"print(f\"Samplerate is {sr}\")\n",
"print(f\"Length of the signal in seconds: {data.shape[0]/sr}\")"
]
},
{
"cell_type": "markdown",
"id": "e6f6272e-62fa-4398-9adf-c4f9c15bc404",
"metadata": {},
"source": [
"## FFT\n",
"\n",
"Before we talk about the infamous *FFT* we need to talk about Fourier analysis, especially the [Fourier series](https://en.wikipedia.org/wiki/Fourier_series).\n",
"The Fourier series allows us to deconstruct any given periodic function into a series of summed sines and cosines, so\n",
"\n",
"$$\n",
"f(t) \\approx \\sum_{k=0}^{\\infty} a_k \\cos(kt) + b_k \\sin(kt)\n",
"$$\n",
"\n",
"where $a_k, b_k \\in \\mathbb{R}$.\n",
"This may seem a bit underhelming at first but it is indeed not trivial that we can describe any periodic (*\"repeating\"*) function as a sum of sine and cosines.\n",
"As we increase the number of $k \\in \\mathbb{N}$ we increase the *resolution* and therefore come closer to the signal we want to approximate.\n",
"\n",
"But there is one problem: The audio signal we want to analyze is most of the time not periodic but aperiodic.\n",
"It turns out that this problem can be solved by the [Fourier transformation](https://en.wikipedia.org/wiki/Fourier_transform) which is generalization of the former mentioned Fourier series.\n",
"This Fourier transformation works on continous signals and returns a continous spectrum (the space of time, frequency and amplitude) which works well in a mathematical context but computers need discrete signals - in the end computers work binary and rely on discrete values.\n",
"For more information about the resolution we refer to the groundbreaking paper *On Computable Numbers, with an Application to the Entscheidungsproblem* by Alan Turing {cite}`turing1936a`.\n",
"The [difference between discrete and continious](https://en.wikipedia.org/wiki/Continuous_or_discrete_variable) is indeed not trivial as it may seems and for a proper definition of these terms rely on measure theory and [topology](https://en.wikipedia.org/wiki/Continuous_function#Continuous_functions_between_topological_spaces).\n",
"A classical example form statistics is that the height of a person is continous (as there are infinite many [infinitesimal](https://en.wikipedia.org/wiki/Infinitesimal) possibilities between $160~\\text{cm}$ and $161~\\text{cm}$) and the outcome of rolling a traditional dice is finite (1,2,3,4,5,6) is discrete.\n",
"\n",
"Thankfully there is a discrete version of the Fourier transformation, called the *Discrete Fourier transformation* (DFT) of which an algorithmic optimized version called *Fast Fourier Transformation (FFT)* exists.\n",
"So in the end the FFT allows us analyze the finite set of frequencies of sine waves which make up a discrete signal within a given time range.\n",
"*Sweeping* over the signal in time range *windows* with a step size (also called hop size) is called *Short-time Fourier transformation*.\n",
"But this *time range* is also more complicated than it may seem at first as introduces a new problem: We now have a trade off between frequency and time resolution:\n",
"The shorter the selected time range, the more precise we know when a certain frequency is audible.\n",
"At the same time if our selected time range is too small we can not pick up low frequencies.\n",
"This is famously known as the *uncertainty principle*, see [Gabor limit](https://en.wikipedia.org/wiki/Uncertainty_principle#Signal_processing) for more information on it.\n",
"A way to tackle this problem is the [discrete wavelet transformation](https://en.wikipedia.org/wiki/Discrete_wavelet_transform) which we will not cover here.\n",
"\n",
"It is also possible to convert the FFT back into the PCM format, which is called the inverse Fourier transformation.\n",
"\n",
"We now know the motivation behind FFT but we have not talked about what it actually returns.\n",
"Depending on the chosen window and hop size we have $T = \\lbrace t_0, t_1, \\dots t_n \\rbrace$ windows and depending on our sample rate $r$ we have $B = \\lbrace b_0, b_1, \\dots b_m \\rbrace$ frequency bins.\n",
"\n",
"For each tuple $(t_i, b_j)$ the FFT gives us a value $z$.\n",
"This value is covering not just the amplitude of the associated sine wave for the given tuple but also gives information about the position of the sine wave that currently occurs.\n",
"If we do not respect the position of the sine wave we can not tell `SinOsc.ar(freq: 200.0, phase: 0.5*pi)` and `SinOsc.ar(freq: 200.0, phase: 0.0)` apart and we can not re-construct the original PCM signal perfectly anymore.\n",
"\n",
"Instead of storing the necesarry tuple of phase and amplitude in 2 numbers we can use the complex numbers.\n",
"Let $x, y \\in \\mathbb{R}$.\n",
"A complex number $z := x + iy$ consisits of a real part $x$ and an imaginary part $y$, where $i := \\sqrt{-1} \\notin \\mathbb{R}$.\n",
"We can therefore imagine the complex number as a 2 dimensional plane.\n",
"The notion only $i$ introduces convenient circumstances which we will not discuss in detail here but it spawned a whole new field of mathematics called *Funktionentheorie*, see {cite}`Königsberger2013` for more information on this topic.\n",
"\n",
"Let us define $z_{i, j}$ as the associated complex number for the $i$-th time window and $j$-th frequency bin.\n",
"The amplitude of this tuple is defined by $\\text{amp}(z) := |z| = \\sqrt{x^2 + y^2}$ which is equal to the [euclidean distance](https://en.wikipedia.org/wiki/Euclidean_distance) of a 2 dimensional point $(x, y)$ to the 0-point $(0, 0)$.\n",
"The phase of our sine wave is defined by $\\text{phase}(z) := \\text{arctan2}(x, y)$ where [arctan2](https://en.wikipedia.org/wiki/Atan2) is a modified version of the arctangent.\n",
"It measures the angle between the 2 vectors $(x, y)$ and $(1, 0)$.\n",
"\n",
"A good visual representation of the FFT is availaible on the website [Seeing circles, sines and signals](https://jackschaedler.github.io/circles-sines-signals/dft_introduction.html).\n",
"\n",
"Before we continue with our examples in Python we will need to define 2 constants for this notebook which we will need to define to perform the STFT.\n",
"`WIN_LENGTH` defines is the size of our window in samples, `HOP_LENGTH` the number of samples between each window and `N_FFT` the number of frequency bins used to analyse the signal which needs to be at least half the sample size.\n",
"As we want the windows of the STFT to overlap we need to choose a smaller hop length than the window length."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "7f0bfeef-7361-4122-ac1a-b86b99436343",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(11026, 123)"
]
},
"execution_count": null,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"WIN_LENGTH = int(sr/4)\n",
"HOP_LENGTH = int(sr/6)\n",
"N_FFT = int(sr/2)\n",
"\n",
"data_fft = librosa.stft(data, n_fft=N_FFT, hop_length=HOP_LENGTH, win_length=WIN_LENGTH)\n",
"data_fft.shape"
]
},
{
"cell_type": "markdown",
"id": "1611a4fe-92c3-4b29-a29f-cef0e20db8d4",
"metadata": {},
"source": [
"The mono signal got transfered to a 2 dimensional array where the first dimension gives us the number of frequency bins and the second dimension the number of time windows used.\n",
"\n",
"We can now take a look at the raw data of the STFT."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "fd94f288-f5a6-44ff-bb14-9e4f63b9fe1b",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([[ 9.66972053e-01+0.0000000e+00j, -8.16252291e-01+0.0000000e+00j,\n",
" -3.73325378e-01+0.0000000e+00j, ...,\n",
" -6.43045232e-02+0.0000000e+00j, -1.68261409e-01+0.0000000e+00j,\n",
" 1.67719409e-01+0.0000000e+00j],\n",
" [-1.20298362e+00+1.5370498e-04j, 7.32871294e-01+1.0257830e-01j,\n",
" 2.27439582e-01+5.3965297e-02j, ...,\n",
" 3.31833810e-02-2.7371196e-02j, 1.09452739e-01-1.4467488e-01j,\n",
" -2.07813919e-01+2.3785539e-02j],\n",
" [ 1.78492057e+00-2.3307584e-04j, -5.49815536e-01-1.4575520e-01j,\n",
" 1.41225010e-01-1.4921878e-01j, ...,\n",
" 4.48968671e-02+6.2899567e-02j, 2.30876897e-02+1.9451287e-01j,\n",
" 2.71128535e-01-2.9871117e-02j],\n",
" ...,\n",
" [ 5.46663720e-03-4.7056403e-07j, 2.62198132e-03+3.2184920e-03j,\n",
" 1.10907322e-02-7.4709475e-04j, ...,\n",
" 1.42986812e-02-1.1434390e-02j, -6.79174345e-03-2.4079946e-03j,\n",
" -7.96438847e-03+1.9497715e-04j],\n",
" [-6.69784704e-03-1.8464819e-07j, -1.50201190e-03-2.0925521e-03j,\n",
" -1.30195906e-02+8.5797138e-04j, ...,\n",
" -1.55045651e-02+7.2950507e-03j, 1.60112455e-02+2.0183965e-03j,\n",
" 1.51164755e-02-1.9905735e-04j],\n",
" [ 7.11801788e-03+0.0000000e+00j, 1.18425570e-03+0.0000000e+00j,\n",
" 1.37792919e-02+0.0000000e+00j, ...,\n",
" 1.59060229e-02+0.0000000e+00j, -1.98078137e-02+0.0000000e+00j,\n",
" -1.83464885e-02+0.0000000e+00j]], dtype=complex64)"
]
},
"execution_count": null,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"data_fft"
]
},
{
"cell_type": "markdown",
"id": "9c031872-2729-4b3b-bc01-a04f8055e51a",
"metadata": {},
"source": [
"We can take a look at the amplitudes (also called magnitude) of each bin for the first window.\n",
"We use pandas to plot the distribution of the values for us."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "8e3c3553",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Amplitude values: [0.96697205 1.2029836 1.7849206 ... 0.00546664 0.00669785 0.00711802]\n"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAA4IAAAHgCAYAAADwqtMZAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Z1A+gAAAACXBIWXMAAAsTAAALEwEAmpwYAAAiaUlEQVR4nO3df7Rmd10f+vfHhKAr6EEbrstChglMpEa8RTo30PpjYa/axHgSir+I1gs0y7lxFSjL1hqLrVTb3lBL15UrrWta01AvJsQf2IREArdtzOoVagKmkBDREGJJSg0/DyAKhnz6x3mGnIxzzjzzZPZ5zjPf12uts+bZ3+fXe86efZ7znr33d1d3BwAAgHF80bIDAAAAsLsUQQAAgMEoggAAAINRBAEAAAajCAIAAAxGEQQAABjM6csO8FicddZZvX///mXHAAAAWIp3vvOdH+nuJ53o81a6CO7fvz+33377smMAAAAsRVX9wSLPc2goAADAYBRBAACAwSiCAAAAg1EEAQAABqMIAgAADEYRBAAAGIwiCAAAMBhFEAAAYDCKIAAAwGAUQQAAgMEoggAAAINRBAEAAAajCAIAAAxGEQQAABiMIggAADAYRRAAAGAwiiAAAMBg9lQRrKozq+r2qvrOZWcBAAA4VZ0+5YtX1VVJvjPJg939zC3jFyT52SSnJfk33X3l7K4fS3LdlJl2w/4rbtzx/vuuvGiXkgAAAPxZU+8RvDrJBVsHquq0JK9LcmGS85JcWlXnVdW3JXlvkgcnzgQAADC0SfcIdvetVbX/qOHzk9zT3fcmSVVdm+SSJE9IcmY2y+EfV9VN3f3wlPkAAABGNGkR3MaTk3xwy/L9SZ7T3S9Nkqp6cZKPbFcCq+pQkkNJsm/fvmmTAgAAnIL21GQxSdLdV3f3m3e4/3B3H+zug0960pN2MxoAAMApYRlF8IEkZ29ZfspsDAAAgF2wjCJ4W5Jzq+qcqjojyQuTXL+EHAAAAEOatAhW1TVJ3p7kGVV1f1Vd1t0PJXlpkpuT3J3kuu6+a8ocAAAAPGLqWUMv3Wb8piQ3Lfq6VbWeZP3AgQOLvgQAAMCw9txkMfPo7hu6+9Da2tqyowAAAKyclSyCAAAALE4RBAAAGIwiCAAAMJiVLIJVtV5Vhzc2NpYdBQAAYOWsZBE0WQwAAMDiVrIIAgAAsDhFEAAAYDCKIAAAwGBWsgiaLAYAAGBxK1kETRYDAACwuJUsggAAACxOEQQAABiMIggAADAYRRAAAGAwiiAAAMBgVrIIunwEAADA4layCLp8BAAAwOJWsggCAACwOEUQAABgMIogAADAYBRBAACAwSiCAAAAg1EEAQAABrOSRdB1BAEAABa3kkXQdQQBAAAWt5JFEAAAgMUpggAAAINRBAEAAAajCAIAAAxGEQQAABiMIggAADAYRRAAAGAwK1kEXVAeAABgcStZBF1QHgAAYHErWQQBAABYnCIIAAAwGEUQAABgMIogAADAYBRBAACAwSiCAAAAg1EEAQAABqMIAgAADEYRBAAAGIwiCAAAMJiVLIJVtV5Vhzc2NpYdBQAAYOWsZBHs7hu6+9Da2tqyowAAAKyclSyCAAAALE4RBAAAGIwiCAAAMBhFEAAAYDCKIAAAwGAUQQAAgMEoggAAAINRBAEAAAajCAIAAAxGEQQAABiMIggAADAYRRAAAGAwiiAAAMBgFEEAAIDBKIIAAACDUQQBAAAGs5JFsKrWq+rwxsbGsqMAAACsnJUsgt19Q3cfWltbW3YUAACAlbOSRRAAAIDFKYIAAACDUQQBAAAGowgCAAAMRhEEAAAYjCIIAAAwGEUQAABgMIogAADAYBRBAACAwSiCAAAAg1EEAQAABqMIAgAADEYRBAAAGIwiCAAAMBhFEAAAYDCKIAAAwGAUQQAAgMEoggAAAINRBAEAAAajCAIAAAxGEQQAABiMIggAADAYRRAAAGAwe6YIVtXXVNXPV9WvVNUPLzsPAADAqWrSIlhVV1XVg1V151HjF1TV+6rqnqq6Ikm6++7uvjzJ9yb5hilzAQAAjGzqPYJXJ7lg60BVnZbkdUkuTHJekkur6rzZfRcnuTHJTRPnAgAAGNakRbC7b03ysaOGz09yT3ff292fS3Jtkktmj7++uy9M8gNT5gIAABjZ6Ut4zycn+eCW5fuTPKeqnpfkBUkenx32CFbVoSSHkmTfvn2ThQQAADhVLaMIHlN335LkljkedzjJ4SQ5ePBgT5sKAADg1LOMWUMfSHL2luWnzMYAAADYBcsogrclObeqzqmqM5K8MMn1S8gBAAAwpKkvH3FNkrcneUZV3V9Vl3X3Q0lemuTmJHcnua6775oyBwAAAI+Y9BzB7r50m/Gb8hguEVFV60nWDxw4sOhLAAAADGsZh4Y+Zt19Q3cfWltbW3YUAACAlbOSRRAAAIDFKYIAAACDUQQBAAAGs5JFsKrWq+rwxsbGsqMAAACsnJUsgiaLAQAAWNxKFkEAAAAWpwgCAAAMRhEEAAAYzEoWQZPFAAAALG4li6DJYgAAABa3kkUQAACAxSmCAAAAg1EEAQAABqMIAgAADEYRBAAAGMxKFkGXjwAAAFjcShZBl48AAABY3EoWQQAAABanCAIAAAxGEQQAABiMIggAADAYRRAAAGAwiiAAAMBgVrIIuo4gAADA4layCLqOIAAAwOJWsggCAACwOEUQAABgMIogAADAYBRBAACAwSiCAAAAg1EEAQAABqMIAgAADGYli6ALygMAACxuJYugC8oDAAAsbiWLIAAAAItTBAEAAAajCAIAAAxGEQQAABiMIggAADAYRRAAAGAwiiAAAMBgFEEAAIDBKIIAAACDOW4RrKqnV9XjZ7efV1Uvr6onTp4MAACAScyzR/BXk3y+qg4kOZzk7CS/NGmq46iq9ao6vLGxscwYAAAAK2meIvhwdz+U5K8n+X+6+0eTfNW0sXbW3Td096G1tbVlxgAAAFhJ8xTBP62qS5O8KMmbZ2OPmy4SAAAAU5qnCL4kyV9O8k+6+wNVdU6SX5w2FgAAAFM5/XgP6O73VtWPJdk3W/5AkldPHQwAAIBpzDNr6HqSO5K8Zbb8rKq6fuJcAAAATGSeQ0NfleT8JJ9Iku6+I8nTJksEAADApOaaLKa7j75Ow8NThAEAAGB6xz1HMMldVfX9SU6rqnOTvDzJb00bCwAAgKnMs0fwZUm+Nslns3kh+Y0kr5gwEwAAABOaZ9bQzyR55ewLAACAFTfPrKFvq6onbln+8qq6edJUAAAATGaeQ0PP6u5PHFno7o8n+V8mSwQAAMCk5imCD1fVviMLVfXUJD1dJAAAAKY0z6yhr0zyn6vqN5NUkm9KcmjSVAAAAExmnsli3lJVz07y3NnQK7r7I9PGAgAAYCrz7BFMkscn+djs8edVVbr71uliAQAAMJXjFsGqenWS70tyV5KHZ8OdRBEEAABYQfPsEXx+kmd092cnzjK3qlpPsn7gwIFlRwEAAFg588waem+Sx00d5ER09w3dfWhtbW3ZUQAAAFbOPHsEP5Pkjqr6D0m+sFewu18+WSoAAAAmM08RvH72BQAAwClgnstHvL6qviTJvu5+3y5kAgAAYELHPUdwNjHLHUneMlt+VlXZQwgAALCi5pks5lVJzk/yiSTp7juSPG2yRAAAAExqniL4p929cdTYw8d8JAAAAHvePJPF3FVV35/ktKo6N8nLk/zWtLEAAACYyjx7BF+W5GuzeemIX0qykeQVE2YCAABgQjvuEayq05Lc2N3fkuSVuxMJAACAKe24R7C7P5/k4apa26U8AAAATGyecwQ/neQ9VfW2JH90ZLC7Xz5ZKgAAACYzTxH8tdkXAAAAp4DjFsHufv1uBAEAAGB3HLcIVtUHkvTR493tovIAAAAraJ5DQw9uuf3FSb4nyVdMEwcAAICpHfc6gt390S1fD3T3/53koumjAQAAMIV5Dg199pbFL8rmHsJ59iQCAACwB81T6F6z5fZDST6Q5HuniQMAAMDU5pk19Ft2IwgAAAC747jnCFbVP62qJ25Z/vKq+seTpgIAAGAyxy2CSS7s7k8cWejujyf5jskSAQAAMKl5iuBpVfX4IwtV9SVJHr/D4wEAANjD5pks5g1J/kNV/dvZ8kuSvH66SAAAAExpnsliXl1V/zXJt86Gfrq7b54iTFU9P5vXKPyyJL/Q3W+d4n0AAABGNs9kMeckuaW7/253/90kt1bV/nnfoKquqqoHq+rOo8YvqKr3VdU9VXVFknT3r3f3DyW5PMn3ndDfBAAAgLnMc47gLyd5eMvy52dj87o6yQVbB6rqtCSvS3JhkvOSXFpV5215yE/M7gcAAOAkm6cInt7dnzuyMLt9xrxv0N23JvnYUcPnJ7mnu++dvd61SS6pTa9O8hvd/a553wMAAID5zVMEP1xVFx9ZqKpLknzkMb7vk5N8cMvy/bOxl2XzXMTvrqrLj/XEqjpUVbdX1e0f/vCHH2MMAACA8cwza+jlSd5QVUcO1fxgkh+cIkx3vzbJa4/zmMNJDifJwYMHe4ocAAAAp7J5Zg19f5LnVtUTZsufPgnv+0CSs7csP2U2BgAAwMTmmTV0rar+RZJbktxSVa+pqrXH+L63JTm3qs6pqjOSvDDJ9Y/xNQEAAJjDPOcIXpXkU0m+d/b1yST/dsdnbFFV1yR5e5JnVNX9VXVZdz+U5KVJbk5yd5LruvuuEw0PAADAiZvnHMGnd/d3bVn+R1V1x7xv0N2XbjN+U5Kb5n2drapqPcn6gQMHFnk6AADA0ObZI/jHVfWNRxaq6huS/PF0kY6vu2/o7kNra4/1CFUAAIDxzDtr6L/bcl7gx5O8aLpIAAAATGmeWUP/a5K/WFVfNlv+5OSpAAAAmMw8ewSTKIAAAACninnOEdxzqmq9qg5vbGwsOwoAAMDK2bYIVtX3zP48Z/fizMdkMQAAAIvbaY/gj8/+/NXdCAIAAMDu2OkcwY9W1VuTnFNV1x99Z3dfPF0sAAAAprJTEbwoybOT/GKS1+xOHAAAAKa2bRHs7s8leUdV/ZXu/nBVPWE2/uldS7eNqlpPsn7gwIFlRwEAAFg588wa+pVV9TtJ7kry3qp6Z1U9c+JcOzJZDAAAwOLmKYKHk/xIdz+1u/cl+TuzMQAAAFbQPEXwzO7+T0cWuvuWJGdOlggAAIBJ7TRZzBH3VtU/yOakMUnyN5LcO10kAAAApjTPHsG/meRJSX4tm9cUPGs2BgAAwAo67h7B7v54kpfvQhYAAAB2wTx7BPecqlqvqsMbGxvLjgIAALByVrIIunwEAADA4layCAIAALC4454jWFXnJHlZkv1bH9/dF08XCwAAgKnMc/mIX0/yC0luSPLwpGkAAACY3DxF8E+6+7WTJwEAAGBXzFMEf7aqfjLJW5N89shgd79rslSnuP1X3LjtffddedEuJgEAAEY0TxH8uiQ/mOSv5pFDQ3u2DAAAwIqZpwh+T5Kndffnpg4zr6paT7J+4MCBZUcBAABYOfNcPuLOJE+cOMcJcR1BAACAxc2zR/CJSX63qm7Lo88RdPkIAACAFTRPEfzJyVMAAACwa45bBLv7N3cjCAAAALvjuEWwqj6VzVlCk+SMJI9L8kfd/WVTBgMAAGAa8+wR/NIjt6uqklyS5LlThgIAAGA688wa+gW96deT/LVp4gAAADC1eQ4NfcGWxS9KcjDJn0yWCAAAgEnNM2vo+pbbDyW5L5uHhy6NC8oDAAAsbp5zBF+yG0FORHffkOSGgwcP/tCyswAAAKyabYtgVf3DHZ7X3f3TE+QBAABgYjvtEfyjY4ydmeSyJH8uiSIIAACwgrYtgt39miO3q+pLk/ztJC9Jcm2S12z3PAAAAPa2Hc8RrKqvSPIjSX4gyeuTPLu7P74bwQAAAJjGTucI/kySFyQ5nOTruvvTu5YKAACAyex0Qfm/k+TPJ/mJJP+9qj45+/pUVX1yd+IBAABwsu10juBOJREAAIAVpewBAAAMRhEEAAAYjCIIAAAwmJUsglW1XlWHNzY2lh0FAABg5axkEezuG7r70Nra2rKjAAAArJyVLIIAAAAsThEEAAAYjCIIAAAwGEUQAABgMIogAADAYBRBAACAwSiCAAAAg1EEAQAABqMIAgAADEYRBAAAGIwiCAAAMBhFEAAAYDCKIAAAwGAUQQAAgMGcvuwAPNr+K27c8f77rrxol5IAAACnKnsEAQAABrOSRbCq1qvq8MbGxrKjAAAArJyVLILdfUN3H1pbW1t2FAAAgJWzkkUQAACAxSmCAAAAg1EEAQAABqMIAgAADEYRBAAAGIwiCAAAMBhFEAAAYDCKIAAAwGAUQQAAgMEoggAAAINRBAEAAAajCAIAAAxGEQQAABjM6csOwInZf8WNO95/35UX7VISAABgVdkjCAAAMBhFEAAAYDCKIAAAwGAUQQAAgMEoggAAAINRBAEAAAajCAIAAAxGEQQAABiMIggAADCY05cd4IiqelqSVyZZ6+7vXnYeTsz+K27c8f77rrxol5IAAADHM+kewaq6qqoerKo7jxq/oKreV1X3VNUVSdLd93b3ZVPmAQAAYPpDQ69OcsHWgao6LcnrklyY5Lwkl1bVeRPnAAAAYGbSItjdtyb52FHD5ye5Z7YH8HNJrk1yyZQ5AAAAeMQyzhF8cpIPblm+P8lzqurPJfknSb6+qn68u/+vYz25qg4lOZQk+/btmzrrUHY6z885fgAAcOrYM5PFdPdHk1w+x+MOJzmcJAcPHuypcwEAAJxqlnH5iAeSnL1l+SmzMQAAAHbBMorgbUnOrapzquqMJC9Mcv0ScgAAAAxp6stHXJPk7UmeUVX3V9Vl3f1QkpcmuTnJ3Umu6+67pswBAADAIyY9R7C7L91m/KYkNy36ulW1nmT9wIEDi77EkI530XcAAGAMyzg09DHr7hu6+9Da2tqyowAAAKyclSyCAAAALE4RBAAAGIwiCAAAMJiVLIJVtV5Vhzc2NpYdBQAAYOWsZBE0WQwAAMDiVrIIAgAAsDhFEAAAYDCKIAAAwGBOX3aARVTVepL1AwcOLDvKnrP/ihuXHQEAANjjVnKPoMliAAAAFreSRRAAAIDFKYIAAACDUQQBAAAGowgCAAAMRhEEAAAYjMtHMBeXpQAAgFPHSu4RdPkIAACAxa1kEQQAAGBxiiAAAMBgFEEAAIDBKIIAAACDUQQBAAAGowgCAAAMxnUEGdpO10e878qLdjEJAADsnpXcI+g6ggAAAItbySIIAADA4hRBAACAwSiCAAAAg1EEAQAABqMIAgAADEYRBAAAGIwiCAAAMJiVLIJVtV5Vhzc2NpYdBQAAYOWsZBF0QXkAAIDFrWQRBAAAYHGKIAAAwGAUQQAAgMEoggAAAINRBAEAAAajCAIAAAxGEQQAABiMIggAADAYRRAAAGAwiiAAAMBgVrIIVtV6VR3e2NhYdhQAAICVs5JFsLtv6O5Da2try44CAACwclayCAIAALA4RRAAAGAwiiAAAMBgFEEAAIDBKIIAAACDUQQBAAAGowgCAAAMRhEEAAAYjCIIAAAwGEUQAABgMIogAADAYBRBAACAwSiCAAAAg1EEAQAABqMIAgAADEYRBAAAGMzpyw6wiKpaT7J+4MCBZUdhTvuvuHHb++678qJdTLI37PT9SFbze3Iq/p04cbZ1AFgNK7lHsLtv6O5Da2try44CAACwclayCAIAALA4RRAAAGAwiiAAAMBgFEEAAIDBKIIAAACDUQQBAAAGowgCAAAMRhEEAAAYjCIIAAAwGEUQAABgMIogAADAYBRBAACAwSiCAAAAg1EEAQAABqMIAgAADEYRBAAAGIwiCAAAMBhFEAAAYDCKIAAAwGAUQQAAgMEoggAAAINRBAEAAAajCAIAAAzm9GUHOKKqzkzyL5N8Lskt3f2GJUcCAAA4JU26R7CqrqqqB6vqzqPGL6iq91XVPVV1xWz4BUl+pbt/KMnFU+YCAAAY2dSHhl6d5IKtA1V1WpLXJbkwyXlJLq2q85I8JckHZw/7/MS5AAAAhjVpEezuW5N87Kjh85Pc0933dvfnklyb5JIk92ezDE6eCwAAYGTLOEfwyXlkz1+yWQCfk+S1SX6uqi5KcsN2T66qQ0kOJcm+ffsmjMlu2X/FjZO99n1XXrTwc6fMNaWdcj+W78cyHW9dTLmel/VvaMp1dSr+G3ksHuu2vle/Z3v13x+nvil/rq6qx/Jzd5nfz9E+L0b7t7tnJovp7j9K8pI5Hnc4yeEkOXjwYE+dCwAA4FSzjEMwH0hy9pblp8zGAAAA2AXLKIK3JTm3qs6pqjOSvDDJ9UvIAQAAMKSpLx9xTZK3J3lGVd1fVZd190NJXprk5iR3J7muu++aMgcAAACPmPQcwe6+dJvxm5LctOjrVtV6kvUDBw4s+hIAAADDWsnLNHT3Dd19aG1tbdlRAAAAVs5KFkEAAAAWpwgCAAAMRhEEAAAYzEoWwapar6rDGxsby44CAACwclayCJosBgAAYHErWQQBAABYnCIIAAAwGEUQAABgMCtZBE0WAwAAsLiVLIImiwEAAFjcShZBAAAAFqcIAgAADEYRBAAAGEx197IzLKyqPpzkD5ad4xjOSvKRZYfgC6yPvcO62Fusj73F+tg7rIu9xfrYW6yPvePIunhqdz/pRJ+80kVwr6qq27v74LJzsMn62Dusi73F+thbrI+9w7rYW6yPvcX62Dse67pwaCgAAMBgFEEAAIDBKILTOLzsADyK9bF3WBd7i/Wxt1gfe4d1sbdYH3uL9bF3PKZ14RxBAACAwdgjCAAAMBhF8CSrqguq6n1VdU9VXbHsPCOpqrOr6j9V1Xur6q6q+tuz8VdV1QNVdcfs6zuWnXUUVXVfVb1n9n2/fTb2FVX1tqr6/dmfX77snKe6qnrGln//d1TVJ6vqFbaN3VNVV1XVg1V155axY24Ltem1s8+Rd1fVs5eX/NS0zfr4mar63dn3/E1V9cTZ+P6q+uMt28nPLy34KWqb9bHtz6eq+vHZ9vG+qvpry0l9atpmXbxxy3q4r6rumI3bNia2w++2J+Xzw6GhJ1FVnZbk95J8W5L7k9yW5NLufu9Sgw2iqr4qyVd197uq6kuTvDPJ85N8b5JPd/c/X2a+EVXVfUkOdvdHtoz9syQf6+4rZ/9Z8uXd/WPLyjia2c+pB5I8J8lLYtvYFVX1zUk+neTfdfczZ2PH3BZmv/C+LMl3ZHM9/Wx3P2dZ2U9F26yPb0/yH7v7oap6dZLM1sf+JG8+8jhOvm3Wx6tyjJ9PVXVekmuSnJ/kzyf5/5J8dXd/fldDn6KOtS6Ouv81STa6+6dsG9Pb4XfbF+ckfH7YI3hynZ/knu6+t7s/l+TaJJcsOdMwuvtD3f2u2e1PJbk7yZOXm4pjuCTJ62e3X5/NH2jsnv89yfu7+w+WHWQk3X1rko8dNbzdtnBJNn8J6+5+R5Inzn4Z4CQ51vro7rd290OzxXckecquBxvUNtvHdi5Jcm13f7a7P5Dknmz+/sVJsNO6qKrK5n+uX7OroQa2w++2J+XzQxE8uZ6c5INblu+PIrIUs/+l+vok/2U29NLZLvKrHIq4qzrJW6vqnVV1aDb2ld39odnt/5HkK5cTbVgvzKM/xG0by7PdtuCzZPn+ZpLf2LJ8TlX9TlX9ZlV907JCDehYP59sH8vzTUn+sLt/f8uYbWOXHPW77Un5/FAEOeVU1ROS/GqSV3T3J5P8qyRPT/KsJB9K8prlpRvON3b3s5NcmORvzQ45+YLePDbd8em7pKrOSHJxkl+eDdk29gjbwt5RVa9M8lCSN8yGPpRkX3d/fZIfSfJLVfVly8o3ED+f9p5L8+j/SLRt7JJj/G77BY/l80MRPLkeSHL2luWnzMbYJVX1uGxuKG/o7l9Lku7+w+7+fHc/nORfxyEku6a7H5j9+WCSN2Xze/+HRw5TmP354PISDufCJO/q7j9MbBt7wHbbgs+SJamqFyf5ziQ/MPvlKrNDED86u/3OJO9P8tVLCzmIHX4+2T6WoKpOT/KCJG88Mmbb2B3H+t02J+nzQxE8uW5Lcm5VnTP7n/cXJrl+yZmGMTt2/ReS3N3d/2LL+NZjo/96kjuPfi4nX1WdOTuxOVV1ZpJvz+b3/vokL5o97EVJ/v1yEg7pUf+ba9tYuu22heuT/B+z2d+em82JGT50rBfg5KmqC5L8vSQXd/dntow/aTbJUqrqaUnOTXLvclKOY4efT9cneWFVPb6qzsnm+vjt3c43oG9N8rvdff+RAdvG9Lb73TYn6fPj9AkyD2s209hLk9yc5LQkV3X3XUuONZJvSPKDSd5zZGrjJH8/yaVV9axs7ja/L8n/uYxwA/rKJG/a/BmW05P8Une/papuS3JdVV2W5A+yeeI5E5uV8W/Lo//9/zPbxu6oqmuSPC/JWVV1f5KfTHJljr0t3JTNGd/uSfKZbM7uykm0zfr48SSPT/K22c+td3T35Um+OclPVdWfJnk4yeXdPe/EJsxhm/XxvGP9fOruu6rquiTvzeYhvH/LjKEnz7HWRXf/Qv7s+eWJbWM3bPe77Un5/HD5CAAAgME4NBQAAGAwiiAAAMBgFEEAAIDBKIIAAACDUQQBAAAGowgCsOdV1fOrqqvqL5zE13xeVb15dvviqrpiy3udt8Dr3VJVB09WvmO8/tVV9d1TvT4AY1EEAVgFlyb5z7M/T7ruvr67r5wtPj/JCRdBAFgliiAAe1pVPSHJNya5LJsXNT4y/ryq+s2q+vdVdW9VXVlVP1BVv11V76mqp88ed3VV/XxV3V5Vv1dV33mM93hxVf1cVf2VJBcn+ZmquqOqnr51T19VnVVV981uf0lVXVtVd1fVm5J8yZbX+/aqentVvauqfnn2d9j6fn+hqn57y/L+qnrP7PY/rKrbqurOqjpcs6ubH/X8+6rqrNntg1V1y+z2mVV11ex78DtVdcli33UATnWKIAB73SVJ3tLdv5fko1X1l7bc9xeTXJ7ka5L8YJKv7u7zk/ybJC/b8rj9Sc5PclGSn6+qLz7WG3X3byW5PsmPdvezuvv9O+T64SSf6e6vSfKTSf5SslkWk/xEkm/t7mcnuT3Jjxz1Pr+b5IyqOmc29H1J3ji7/XPd/b919zOzWS7/THHdwSuT/MfZ9+BbsllozzyB5wMwCEUQgL3u0iTXzm5fm0cfHnpbd3+ouz+b5P1J3jobf082y98R13X3w939+0nuTXIyzjX85iT/b5J097uTvHs2/txsHlr6/1fVHUlelOSpx3j+ddksgMmji+C3VNV/me0h/KtJvvYEMn17kitm73tLki9Osu8Eng/AIE5fdgAA2E5VfUU2y9DXVVUnOS1JV9WPzh7y2S0Pf3jL8sN59GdcH/XSRy/v5KE88h+nx9yTeHTsJG/r7uOdz/jGJL9cVb+WpLv792d7Kv9lkoPd/cGqetU277ldpkryXd39vjlyAjAwewQB2Mu+O8kvdvdTu3t/d5+d5ANJvukEX+d7quqLZucNPi3JTkXpU0m+dMvyfZkd9jnLc8StSb4/SarqmUn+19n4O5J8Q1UdmN13ZlV99dFvMjvs9PNJ/kEe2Rt4pNR9ZHZe4XazhG7N9F1bxm9O8rIj5xVW1ddv95cEYGyKIAB72aVJ3nTU2K/mxGcP/W9JfjvJbyS5vLv/ZIfHXpvkR2eTrTw9yT9P8sNV9TtJztryuH+V5AlVdXeSn0ryziTp7g8neXGSa6rq3Unenu0PRX1jkr+RzcNE092fSPKvk9yZzVJ32zbP+0dJfraqbs9mmTzip5M8Lsm7q+qu2TIA/BnVfSJHxwDAaqmqq5O8ubt/ZdlZAGCvsEcQAABgMPYIAgAADMYeQQAAgMEoggAAAINRBAEAAAajCAIAAAxGEQQAABiMIggAADCY/wk25p00+szSSwAAAABJRU5ErkJggg==\n",
"text/plain": [
"