Skip to content

Instantly share code, notes, and snippets.

@philpem
Created October 19, 2018 14:36
Show Gist options
  • Save philpem/b24bfb98f1fd39e856ea794a3f9f36e6 to your computer and use it in GitHub Desktop.
Save philpem/b24bfb98f1fd39e856ea794a3f9f36e6 to your computer and use it in GitHub Desktop.
Root raised cosine filter using SciPy
Display the source blob
Display the rendered blob
Raw
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Root raised-cosine filter design\n",
"\n",
"From CML FX919B datasheet: https://www.cmlmicro.com/wp-content/uploads/2017/06/FX919B_ds.pdf\n",
"\n",
"## First calculate and plot the desired frequency response"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"data": {
"application/vnd.jupyter.widget-view+json": {
"model_id": "2b7e03f4ded84fcebb6d20b7f890fdd0",
"version_major": 2,
"version_minor": 0
},
"text/plain": [
"interactive(children=(FloatSlider(value=1.0, description='symrate', max=3.0, min=-1.0), FloatSlider(value=0.2,…"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"from __future__ import print_function\n",
"\n",
"%matplotlib inline\n",
"\n",
"from ipywidgets import interact, interactive, fixed, interact_manual\n",
"import ipywidgets as widgets\n",
"\n",
"import matplotlib\n",
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"\n",
"\n",
"def rrc(symrate, b, f):\n",
" T = 1.0/symrate\n",
"\n",
" if f < (1.-b)/(2.0*T):\n",
" return 1.\n",
" elif f <= (1+b)/(2.0*T):\n",
" return np.sqrt(0.5 * (1.0 - np.sin(np.pi * T * (f - 0.5/T)/b)))\n",
" else:\n",
" return 0.\n",
"\n",
"\n",
"@interact(symrate=1.0, b=0.2)\n",
"def g(symrate, b):\n",
" f = np.linspace(0.0, 1.0, 100)\n",
" r = np.array([rrc(symrate, b, x) for x in f])\n",
" plt.plot(f, r)\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Next calculate a FIR filter to model the desired filter response"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAY0AAAEWCAYAAACaBstRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvqOYd8AAAIABJREFUeJzt3Xl8XHW5+PHPM1mbpUmaNKFtmqQbLS3QNoFSkKUVkIJaUFlaoeL1cnu5iqi4gCD8uOjPFfWicr36Q69eCpRFWcQCsrQCXgtt032l1Eya7kmTNPs2z++PmUmHNMtMMjNnJnner1fonDnfOedhtme+57uJqmKMMcYEw+V0AMYYY+KHJQ1jjDFBs6RhjDEmaJY0jDHGBM2ShjHGmKBZ0jDGGBM0SxrGGGOCZknDmB5EpEJEWkSkMeBvvIiUiIiKSKKv3O9EpN23/7iIvCoiM/o57v0isiJ6/yfGhJ8lDWN693FVzQj4O9hHuR+qagYwATgA/CZ6IRoTfZY0jAkDVW0BngLm9LZfRBYBdwM3+Gomm333/5OI7BSRBhHZJyL/GvCYBSJSJSJ3i0i1rwZ0Y8D+q0Rkh++xB0Tka5H9vzQGEp0OwJjhQETSgaXA3t72q+rLIvJdYKqq3hSw6yjwMWAfcDHwkoisU9Vy3/7TgDy8NZn5wCoRWa+qu/HWaq5X1bdEJAeYFIn/N2MCWU3DmN49JyJ1vr/n+in3NRGpAxqAC4FloZxEVf+squ+r11+BvwAX9Sh2r6q2+fb/Gbjed38HMFNERqtqbUCiMSZiLGkY07trVDXb93dNP+UeVNVsoARoAaaHchIRuVJE1voa0uuAq/DWLPxqVbUpYNsNjPfd/pSvvFtE/ioi54dybmMGw5KGMWGgqpXAl4CHRGRUX8UCN0QkBfgD8CBQ4Es+qwAJKJbju/TlVwQc9J1znapeDeQDz+FtUzEmoixpGBMmqvoq3i/05X0UOQKUiIj/c5cMpADHgE4RuRL4SC+P+3cRSRaRi/C2fzzt275RRLJUtQM4AXjC+f9jTG8saRgTXj8CvuGrRfT0tO/fGhEpV9UG4Ha8NYRa4NPACz0ec9i37yDwGHCrqu7y7VsGVIjICeBW4EaMiTCxRZiMiU0isgBYoaqFTsdijJ/VNIwxxgTNkoYxxpig2eUpY4wxQbOahjHGmKANu2lE8vLytKSkZNCPb2pqIj09feCCUWZxhcbiCo3FFZrhGNeGDRuqVXXsgAVVdVj9lZWV6VCsXr16SI+PFIsrNBZXaCyu0AzHuID1GsR3rF2eMsYYEzRLGsYYY4JmScMYY0zQLGkYY4wJmiUNY4wxQbOkYYwxJmiWNIwxxgTNkoYxUbJm91EONdqSFya+WdIwJgo6ujzcumIDz7zX7nQoxgyJJQ1jomDnoRO0dnh4v86D2iShJo5Z0jAmCsrdtQDUtSkH6locjsaYwbOkYUwUlFfWkeiS7tvGxCtLGsZEQXllLQtn5JOccLLWYUw8sqRhTIQdbWilqraF8yaNYXKWi42VljRM/HI0aYjIIhHZLSJ7ReSuXvZ/VkSOicgm398tTsRpzFCUu72Xo+YW5TA1O4HtB0/Q2tHlcFTGDI5jSUNEEoCHgSuBmcBSEZnZS9EnVXWO7++RqAZpTBhsrKwlOcHFmRNGMyXbRadH2VJV73RYxgyKkzWNecBeVd2nqu3ASuBqB+MxJiI2uGuZNWE0KYkJTMlOALxtHMbEI3Gqz7iIXAssUtVbfNvLgPNU9baAMp8FvgccA/YAX1HV/b0cazmwHKCgoKBs5cqVg46rsbGRjIyMQT8+Uiyu0MRKXJ0e5dbXmrl0YiJLz0ihsbGRb5e7mJDh4vbSVKfD6xYrz1dPFldohhLXwoULN6jqOQMWDGZ5v0j8AdcCjwRsLwN+0aNMLpDiu/2vwBsDHdeWe40ui6t/GytrtfjOF/XFzQdV1RvXV1Zu1LJvv6oej8fh6E6KleerJ4srNMN9udcDwMSA7ULffd1UtUZV23ybjwBlUYrNmLDwd68tLc7uvm9ucQ7VjW1U1dogPxN/nEwa64BpIjJJRJKBJcALgQVEZFzA5mJgZxTjM2bIyitrGZeVyrisUd33lRZld+8zJt44ljRUtRO4DXgFbzJ4SlW3i8gDIrLYV+x2EdkuIpuB24HPOhOtMYOzsbKO0qKcD9w3vSCTtOQEG+Rn4lKikydX1VXAqh733Rdw+5vAN6MdlzHhcOREKwfqWvinD5V84P7EBBezC7NtOhETl2xEuDER4q9JlBXnnLKvtDibHYdO0NzeGe2wjBkSSxrGREh5ZS3JiS5mjc86ZV9pUQ5dNsjPxCFLGsZESHllHWdNyCI58dSP2VxfO4c1hpt4Y0nDmAho6+xi64H67p5SPY1JT2ZyXnr3vFTGxAtLGsZEwPaDJ2jv9JzScyrQ3KIcNlbW2kp+Jq5Y0jAmAk4O6us7aZQWZ1PT1E7l8eZohWXMkFnSMCYCNlbWMSF7FAWj+55fqtTaNUwcsqRhTASUV9Yyt4/2DL/TCzLJSEm0dg0TVyxpGBNmh+pbOFTf2m97BkCCS5g9MctqGiauWNIwJsz8NYf+2jP8Soty2HW4wQb5mbhhScOYMCuvrCUl0cXMcaMHLOsf5Ld5vw3yM/HBkoYxYVZeWcvZhb0P6utprs14a+KMJQ1jwqits4vtB04M2J7hl52WzOSx6TbjrYkbljSMCaNtB07Q3uXpniYkGKVFOWzcX2eD/ExcsKRhTBhtrDx1pb6BlBXncLypnYoaG+RnYp8lDWPCaIO7lsKcUeRn9j2or6fuQX52icrEAUsaxoSJqlJeWRt0e4bftPwMMlMSrTHcxAVLGsaEycH6Vo6caOtzZtu+uFzCnCJbyc/EB0saxoRJMJMU9mVuUQ67D5+gsc0G+ZnYZknDmDApr6wlNcnFGUEM6uuptCgbj8KW/VbbMLHNkoYxYVJeWcfZhdkkJYT+sZo70Vs72WCN4SbGWdIwJgxaO7rYcbA+5EZwv6y0JKbmZ1hjuIl5ljSMCYNtB+rp6NKQG8EDlRZl2yA/E/McTRoiskhEdovIXhG5q59ynxIRFZFzohmfMcEqrxx8I7hfaVEOdc0d7KtuCldYxoSdY0lDRBKAh4ErgZnAUhGZ2Uu5TOBLwDvRjdCY4JW76ygak0ZeRsqgj+FPODbIz8QyJ2sa84C9qrpPVduBlcDVvZT7NvADoDWawRkTLFVlQ2XtkC5NAUwdm0FmaqKN1zAxTZy6fioi1wKLVPUW3/Yy4DxVvS2gTClwj6p+SkTWAF9T1fW9HGs5sBygoKCgbOXKlYOOq7GxkYyMjEE/PlIsrtBEM65jzR6+/mYLN52RzGXFSUOK68H1rdS1evjOhWnhDrNf9jqGZjjGtXDhwg2qOnATgKo68gdcCzwSsL0M+EXAtgtYA5T4ttcA5wx03LKyMh2K1atXD+nxkWJxhSaacT23sUqL73xRt1bVDVh2oLh++upuLbnrRT3R0h6m6IJjr2NohmNcwHoN4rvbyctTB4CJAduFvvv8MoEzgTUiUgHMB16wxnATazZW1jEqKYEZp2UO+VilRTmoYiv5mZjlZNJYB0wTkUkikgwsAV7w71TVelXNU9USVS0B1gKLtZfLU8Y4yb9SX+IgBvX1NKcoGxEb5Gdil2NJQ1U7gduAV4CdwFOqul1EHhCRxU7FZUwovIP6TlA2hK62gUanJjHNBvmZGJbo5MlVdRWwqsd99/VRdkE0YjImFFuq6un06KBHgvemtCiHVVsP4fEoLpeE7bjGhIONCDdmCPw1grlD7G4bqLQohxOtneyrbgzbMY0JF0saxgxBubuWktw0cocwqK8n/1Kx5W4br2FijyUNYwZJVSmvrAvrpSmAyXkZjE61lfxMbLKkYcwgVdW2UN3YxtwwNYL7uVzC3KIcSxomJlnSMGaQ/N1ihzp9SG9Ki3J472gjJ1o7wn5sY4bCkoYxg1ReWUtacgLTC4Y+qK+n0uJsVGGTzUNlYowlDWMGqbyyltmF2WEZ1NfTnIneQX52icrEGksaxgxCc3snOw81dPd0CrfM1CSmF2TayHATcyxpGDMIW6rq6fJo2EaC92ZuUQ6b9tfh8dhKfiZ2WNIwZhC6B/VNjFzSKC3KpqG1k73HbJCfiR2WNIwZhHJ3HZPz0slJT47YOWwlPxOLLGkYEyJVZWNlLXPDPKivp8l56WSnJVljuIkpljSMCVHl8WZqmtoj1gjuJyLMnZhty7+amGJJw5gQ+X/5h3v6kN6UFuWw92gj9c02yM/EBksaxoSo3F1HRkoip0dgUF9P/naNjfvtEpWJDZY0jAnRBnctsydmkRCFtS5mT8zGJdglKhMzLGkYE4Kmtk52HT4RlUtTQHeNxnpQmVhhScOYEGyuqsOj0WnP8Cst9g7y67JBfiYGWNIwJgQbfZeJwrlS30DKinJobOvkvaMNUTunMX2xpGFMCMrdtUwZm052WuQG9fV0cpCftWsY51nSMCZIqsrG/eFfqW8gJblpjElPtkF+JiZY0jAmSFW1LRxvamdOFC9NgXeQ35yJ2WypspqGcZ4lDWOCVFHTBMCUsRlRP/fkvHTcNc02461xnKNJQ0QWichuEdkrInf1sv9WEdkqIptE5G0RmelEnMYAVNQ0A1CSmx71cxfnpdPW6eFIQ2vUz21MIMeShogkAA8DVwIzgaW9JIXHVfUsVZ0D/BD4SZTDNKabu7qJ1CQX+ZkpUT93SW4aABXVzVE/tzGBnKxpzAP2quo+VW0HVgJXBxZQ1RMBm+mA1c2NYypqmikek44rCiPBe/LXbty+S2TGOEVUnfkeFpFrgUWqeotvexlwnqre1qPcF4A7gGTgw6r6Xi/HWg4sBygoKChbuXLloONqbGwkIyP616wHYnGFJhJx3f12M6elubi9NHXQxxhsXF0eZfmrzVxRksT108Pf3XckvY7hMBzjWrhw4QZVPWfAgqrqyB9wLfBIwPYy4Bf9lP808PuBjltWVqZDsXr16iE9PlIsrtCEO66uLo9Ou2eV/t8/7xjScYYS18IfrdZbH10/pPP3ZaS8juEyHOMC1msQ391OXp46AEwM2C703deXlcA1EY3ImD4cPtFKe6eHYl/bghOKc9O6G+ONcYqTSWMdME1EJolIMrAEeCGwgIhMC9j8KHDKpSljosHf3daJnlN+xbnpuGua/DVvYxyR6NSJVbVTRG4DXgESgN+q6nYReQBvNekF4DYRuQzoAGqBm52K14xsbt8vfCdrGiW5aTS3d3GssY38zMG3qxgzFI4lDQBVXQWs6nHffQG3vxT1oIzpRUVNE8kJLsZljXIshuI8fw+qZksaxjE2ItyYILirm5k4ZlRUFl7qi//SWEW1dbs1zrGkYUwQKmqaHG3PAJiQ7U1abmsMNw6ypGHMAFQVd00zxQ4njeREFxOyR3U3yhvjBEsaxgzgWEMbLR1dlOQ51wjuV5ybZjUN4yhLGsYMoKK755SzNQ3wtmtUWLdb4yBLGsYM4OQYjdioaTS0dlLX3OF0KGaEsqRhzADcNU0kuoQJ2c51t/Xr7kFl7RrGIZY0jBmAu6aZwpxRJCY4/3Hxt6tYu4ZxivOfAmNiXCz0nPIrzElDxGoaxjkDJg0RKRCR34jIS77tmSLyz5EPzRjnqapvjIbz7RkAqUkJjM8aZTUN45hgahq/wzs/1Hjf9h7gy5EKyJhYUtvcQUNrZ8zUNMA/263VNIwzgkkaear6FOAB70SDQFdEozImRvi/nJ2cqLAnG6thnBRM0mgSkVx8S62KyHygPqJRGRMj3N1JI5ZqGukcb2qnvsW63ZroC2aW2zvwrnMxRUT+BozFu+qeMcNeRXUzIjBxjPPdbf387SuVNc2cVZjlcDRmpBkwaahquYhcAkwHBNitqvYTx4wI7pomxmeNIiUxwelQuhUHjNWwpGGibcCkISKf6XFXqYigqv8ToZiMiRkVNc0xMedUIH/7itsaw40Dgrk8dW7A7VTgUqAcsKRhhj13TRNXnjXO6TA+IC05kfzMFFsv3DgimMtTXwzcFpFsYGXEIjImRtQ3d1Db3BEzYzQClfjWCzcm2gYzIrwJmBTuQIyJNe7jsddzys87VsNqGib6gmnT+BO+7rZ4k8xM4KlIBmVMLPB/KTu9Yl9vSvLSObahiqa2TtJTgrnKbEx4BPNuezDgdifgVtWqCMVjTMxw+9biLhoTe5enTjaGNzNz/GiHozEjSTBtGn+NRiDGxJqKmmZOG53KqOTY6W7r56/9uGuaLGmYqOozaYhIAycvS31gF6Cqau9UM6y5a5piavqQQEX+msZxa9cw0dVnQ7iqZqrq6F7+MsOVMERkkYjsFpG9InJXL/vvEJEdIrJFRF4XkeJwnNeYYFTUNMdkewbA6NQkctOTrQeVibqge0+JSL6IFPn/hnpiEUkAHgauxNu4vlREZvYothE4R1XPBp4BfjjU8xoTjMa2Tqob2yiOsYF9gYpz06iotpqGia5g1tNYLCLvAf8A/gpUAC+F4dzzgL2quk9V2/GO/bg6sICqrlZV/6diLVAYhvMaMyB397rgsVnTABurYZwhqr01WwQUENkMfBh4TVXnishC4CZVHdJCTCJyLbBIVW/xbS8DzlPV2/oo/wvgsKp+p5d9y4HlAAUFBWUrVw5+7GFjYyMZGRmDfnykWFyhGWpc6w538vCmNv79glSKR4evITycz9fze9t5dm8Hv748jeQEiZm4wsniCs1Q4lq4cOEGVT1nwIKq2u8fsN7372bA5b890OOCOO61wCMB28uAX/RR9ia8NY2UgY5bVlamQ7F69eohPT5SLK7QDDWuh1e/p8V3vqgNrR3hCcgnnM/XcxurtPjOF3XP4RNDPtZwfR0jZTjG5f+uH+gvmHEadSKSAbwFPCYiR/GOCh+qA8DEgO1C330fICKXAfcAl6hqWxjOa8yA3NXN5GWkkBHDA+dOznbbzLSCTIejMSNFn20aIvKwiFyIt52hGe8Sry8D7wMfD8O51wHTRGSSiCQDS/Cu2xEYw1zgV8BiVT0ahnMaExT38dhZF7wvJTbbrXFAfz+j9gA/AsbhnTbkCVX9fbhOrKqdInIb3vXHE4Dfqup2EXkAbzXpBd/5M4CnRQSgUlUXhysGY/rirmnmgil5TofRr+y0ZLJGJdl64Saq+kwaqvoQ8JBvbMQS4LciMgp4HFipqnuGenJVXQWs6nHffQG3LxvqOYwJVWtHF4fqW2O+pgHe2oatF26iacAut6rqVtUfqOpcYCnwCWBnxCMzxiGVvlHWRXGQNIpy062mYaIqmHEaiSLycRF5DO/4jN3AJyMemTEOqaiO/TEafiW5aRyobaG90+N0KGaE6G/uqcvx1iyuAt7FO/huuarazxozrLljeEr0nopz0/EoVNU2M3ls7I0bMMNPfzWNbwL/C5yhqotV9XFLGGYkqKhpIjstiay0JKdDGVBJwBTpxkRDfw3hH45mIMbECndNc0yu1tebk2M17PeciY7BLPdqzLBWURP7YzT88jKSSU9OsJqGiRpLGsYEaOvs4mBdS9zUNESEYutBZaLIkoYxAapqW/AocVPTACjJs7EaJnosaRgTwD8lR7zUNMAb6/7jzXR2WbdbE3mWNIwJ4F/UKK5qGrlpdHqUQ/WtTodiRgBLGsYEcNc0kZmSyJj0ZKdDCZr1oDLRZEnDmAAVNc0U56XhmyAzLpQETJFuTKRZ0jAmgLumKa7aMwDyM1NITXLhrraahok8SxrG+HR0eaiqbYmr9gwAl0soHpNuNQ0TFZY0jPE5WNdCp0fjrqYBUJybZosxmaiwpGGMT0UcTVTYU0leOu7jzXg86nQoZpizpGGMj/+XerxdngJvTaO908PhE9bt1kSWJQ1jfCqqmxmVlMDYzBSnQwlZiXW7NVFiScMYH2/PqfjqbutXbFOkmyixpGGMj/t4c1y2ZwCMyxpFcoLLahom4ixpGAN0eZTKmubuX+zxJsElFI4ZhbvaahomsixpGAMcPtFKe5cnLrvb+pXYFOkmChxNGiKySER2i8heEbmrl/0Xi0i5iHSKyLVOxGhGBv9o6njsOeXnHavRjKp1uzWR41jSEJEE4GHgSmAmsFREZvYoVgl8Fng8utGZkcY/RqM4L75rGi0dXRxraHM6FDOMOVnTmAfsVdV9qtoOrASuDiygqhWqugWwhQJMRLlrmkhOdDFudKrToQyavz3GphMxkeRk0pgA7A/YrvLdZ0zUVdQ0UTQmDZcr/rrb+tlYDRMNiU4HEA4ishxYDlBQUMCaNWsGfazGxsYhPT5SLK7QhBrXdnczeaNcEf9/ieTz1elRXAJvlu8kv/H9mIlrKCyu0EQlLlV15A84H3glYPubwDf7KPs74NpgjltWVqZDsXr16iE9PlIsrtCEEpfH49Hp31qlD/xpe+QC8on083XxD9/QLzy2IeTHDYfXMZqGY1zAeg3iO9bJy1PrgGkiMklEkoElwAsOxmNGqKMNbbR2eOK655RfcW66jQo3EeVY0lDVTuA24BVgJ/CUqm4XkQdEZDGAiJwrIlXAdcCvRGS7U/Ga4avC1902nsdo+JXkplFR02Tdbk3EONqmoaqrgFU97rsv4PY6oDDacZmRxR3HU6L3VJybTkNrJ7XNHXG1zrmJHzYi3Ix4FTVNJLqE8dnx293Wr6S72631oDKRYUnDjHjummYmjkkjMSH+Pw7+S2y2ip+JlPj/lBgzRBW+KdGHg4ljRiHiXRvEmEiwpGFGNFXFXRO/U6L3lJKYwPisUVbTMBFjScOMaDVN7TS2dQ6bmgZASV6aTSViIsaShhnRTq4LPjxqGuAfq2E1DRMZljTMiOa/9j+sahq5adQ2d1Df3OF0KGYYsqRhRjR3TRMugcKc4ZM0untQHbfahgk/SxpmRKuoaWZCziiSE4fPR+HkbLfWrmHCb/h8UowZBHdN07BqzwAoGuOtNflXIzQmnCxpmBHNfby5+0t2uBiVnEDB6BSraZiIsKRhRqy65nbqmjuGXU0DrAeViZxhsQiTib7Wji7ueXYbLbXtlJzZREkcrq3tn6hwOPWc8ivJTeONXcecDmNQVJW/76vhmfVVpLZ0sMDpgMwHWNIwg/Lvf9rOH8qrcAmsenANH5qay9J5RXxk5mlx06jsn9QvHhPeQIpz06lurKKxrZOMlPj4mFc3tvGHDVWsXLeff1R7J5Hs9CgLdxzh8pkFTodnfOLj3WRiypPrKnni3f18fsEUpnGQA8kTeeLd/dz2+EZy05O5tqyQG86dyOSxGU6H2i9/TWO4tWnAyR5U7pomZo3Pcjiavnk83lrF4+9W8pfth+noUuaVjOH2S6fy4ekFXP0fr3HHk5t44YsXMmkYJvd4ZEnDhGRrVT33Pr+dC6fm8dWPTOetNw/ziQXT+PyCqby1t5on3qnkN2//g1+9uY/zJ+eyZN5EFp15GimJCU6HfoqKmibGZaWSmhR7sQ2V/5JbZU1zTCaNow2tPLOhipXv7qfyeDPZaUl85vwSls6byNT8zO5yX5iTwv9d38m/rdjAHz9/AWnJ9pXlNHsFTNBqm9q5dcUG8tKTeWjJHBJc0r3P5RIuOX0sl5w+9gNfCF9auYmctCQ+VVrIknlFTM2PndqHu6Z5WLZnwMmkEUs9qDwe5e291TzxbiWv7jhCp0eZP3kMX/3I6Vwx67Rek/fYNBcPLZnLZ//7Xb75x638xw1zEJFejm6ixZKGCUqXR/nSk5s41tDGU7eeT25GSp9l8zNT+fyCqdx68RT+9/0anni3kt/9bwWPvP0P5k0aw6fnFbHozN6/JKLJXdPEZWcMz2vlmalJ5GUkx0QPqqMnWnl6QxUr11Wy/3gLY9KT+dyFk7jh3IlMCeIS5iWnj+WOy07nx6/uobQoh5svKIl80KZPljRMUB56bQ9v7jnGdz9xFnMmZgf1GJdLuHBaHhdOy+NYQxt/KK9i5buVfPnJTWT/KYm7rzyD68+dGOHIe9fQ2kF1Y/uwWBe8L8W56Y6u4OfxKN/4wxae3XiALo9ywZRcvnHFDD4yqyDky5VfWDiVTfvr+PaLOzhzwmjKisdEKGozkPjo5mIc9frOI/zsjb1cV1bI0nmD+5Ifm5nCrZdM4Y2vLuDxW86jJDedb/95B83tnWGONjgn1wUfnpenwHuJyu3g5ak33zvGMxuquP6cQlZ/bQGP/8t8Pj57/KDat1wu4Sc3zGFCzig+/1g5RxtaIxCxCYYlDdOviuomvvzkJmaNH823rzlzyNeTXS7hgql5fOujZ9DQ2skLmw6GKdLQnByjMXxrGiW56Ryqb6W1o8uR869YW0leRjL/vvjMsPR8yhqVxC9vLKO+pYPbHt9IR5cnDFGaUFnScMiJ1g7KK2t5ev1+vvfSTh55ax/tnbH1IWhp7+LWFRtwifBfN5WFtQ2irDiHGadl8uhaN6oatuMGy3/ZZrg2hENAD6rj0a9tHKhr4Y1dR7jh3IlhHbczc/xovvfJs3j3H8f54cu7wnbccDlc38p3V+3kZ6+/x6qth9hzpIG2TmeSdqRYm0YEqSpHTrSx92gj7x9rZO/Rxu7bRxvausslJQgdXcqzGw/w4+tnM+O00Q5G7aWq3P3sVnYfaeC/P3suE8M8lkFEuGl+Md96bhub9tcxtygnrMcfiLumibGZKaTHycC3weie7ba6idMLMgcoHV5PvFMJwNJ5RWE/9ifmFrKxso7/99Y/mDMxh4+ePS7s5wiVqvL8poPc9/w2mtu76PSc/CGU4BKKxqQxZWwGU/LTmTo2g6n5GUzJz2B0apKDUQ+Oo58YEVkEPAQkAI+o6vd77E8B/gcoA2qAG1S1ItpxBuNgXQtbD9R3J4f3jzby/rEmGttOXrPPTElkSn4GF58+lim+N87U/Awm5ozi9V1HuefZrSz++d/4yuWns/ziyR/o0hptK9a6eXbjAe64/HQWTM+PyDmumTuB77+0i0fXuqOeNCpqmod1ewYEDvCLbk2jvdPDynWVfHhGQcTWKfnWR2ey7UA9X39mM9NPy/jA2I5oq2ls41vPbeOlbYcpLcrmx9fPoWB0CvuONZ3yY/Gve47S0XUyoeRnpngTSMD3wdmFWWTGcDJaiCmNAAAV6ElEQVRxLGmISALwMHA5UAWsE5EXVHVHQLF/BmpVdaqILAF+ANwQ/WhP5fEom6vqeG3nEV7feZRdhxu69xWM9r4RPlk6wftG8L0hxmam9NkmcMWs0zinOIdvPbeNH7y8i1d3HObH189xZBTsBnctD7y4gw/PyOe2hVMjdp6MlEQ+MXcCT67fz70fnUlOenLEztWTu6aJi6aNjdr5nJCVlkR2WlLUe1C9vP0w1Y3t3DQ//LUMv+REF/95Yxkf+/lbLH90A89/4UOOfNH+Zfth7n52KydaOrlz0YwP/Ng7c0IWZ0744MDKzi4PlcebfUmkyZtQjjXy7MYD3T8wkxKE+ZNzuXRGPpeeURD2Wv5QOVnTmAfsVdV9ACKyErgaCEwaVwP3+24/A/xCRESduAgONLd38tZ71by+8whv7DpGdWMbCS6hrDiHu6+awbklY4ZU5czNSOE/byzlhc0Hufe5bVz50Jt888ozWDa/OMz/J3071tDG5x/bwLisUfz0+jm4IlzbuWl+MY+udfP0hv0sv3hKRM/l19zeyZETbcO+pgH+2W6jW9NYsdZN0Zg0Lo5wUj4tK5WfLy3lpt+8w9ef3sIvbyqN2sC/pg7ljqc28cfyA8wcN5oVtwR3WTkxwcXksRmnTLGjqhxtaGP34Qbe3lvNazuPcP+fdnD/n3YwvSCTS8/wJpA5E7MdvQIBIA59/yIi1wKLVPUW3/Yy4DxVvS2gzDZfmSrf9vu+MtU9jrUcWA5QUFBQtnLlykHH1djYSEbGyRf0eKuHTUe72HS0ix3Hu+j0wKhEOCsvgTn5iZydl0BGcvhfxNpWD7/d1s7W6i5m5rpYMqmLorzIjqbu8ig/Wt/K+3Ue7p2fStHogRu+ez5fg/Hdd1qobVV+cPEoXGH60PcX1/4GD/f+rYV/m53CeeOi+7spHM9XKP5rcyt76zw8eEn/CTJccVU1ePjW31q4fnoSV00aes0xmLhe+kcHT+5uD9s5B7KtuotHtrRwokP42OQkFk9JIjECX+SHm3zfPcc62VPrwaOQmQyzxyYyNz+BWbkJpCZ+8LxDeR0XLly4QVXPGajcsGgFVNVfA78GOOecc3TBggWDPtYbq1czZuocXttxhNd2HmXHoZO9bD5zfgGXnZHPuZPGkJQQ+Y5n11yhPPHufr7z5x18b5PwwDVTuLasMGK/pr63aie7ju/jx9fN5lNlhUE9Zs2aNQzl+QY4kXOQ25/YiGv8rLC1n/QX18vbDsPfNnDVRedyVmF052UKx/MVivKOPbz7xnucf+FF/Y6PCFdc9z63jeTE/dx1/QLGhOFyYzBxXXKJ0vB4Oc9sO8w1F5VywdS8IZ+3N83tnXxv1S4eXe9mXLqL3/3LBUEPdB2sJb5/65s7WLPnKK/tPMqa3Ud5+0AbyQkuzp+Sy2W+Wsj47FFReX85mTQOAIEjxQp99/VWpkpEEoEsvA3iYXesoY2fvLqblza3UPfK33AJlBblcOeiGVx2Rj5T8zOiPueNiPDp84q4cGoe//LIX/n6M1t4ZfsRvvvJM8nPTA3ruV7aeohfvbmPZfOLg04Y4bJo1mnkZSSzYm1lxBrdA/mn1igaAZenSnLT8ChU1bYENWXHUDS2dfLsxgN87OxxYUkYwRIRfnjtbHYfbuCLT2zkxdsvZFzWqLCeY33Fcb769GYqjzfzzxdO4rxRRyKeMAJlpSVx9ZwJXD1nAh1dHtZVHOf1nUd5fecR7n1+O/c+v50zxo1mZkY7kf5N4mTSWAdME5FJeJPDEuDTPcq8ANwM/B24FngjUu0Z6SkJvLrjCFOzXSy9+EwWzsiP6hu/P0W5adw5L5V9icX88JXdXPHTN/nONWeFravh3qONfO3pzcwtyubej80MyzFDkZzo4oZzJ/LLNe9zoK6FCdnh/cD3VFHTzJj0ZLJGxW4PlXApDpgiPdJJ4zlfY+5NUWyD88tISeRXy8q4+hd/499WlPPkv84Py8zKrR1d/PS1Pfz6zX1MyB7FE/8yn/mTc1mz5mgYoh6cpAQXF0zJ44Ip3kGy7x9r4nVfh5x99ZEf6+XY4D5V7QRuA14BdgJPqep2EXlARBb7iv0GyBWRvcAdwF2RiictOZF37r6M2+am8qmywphJGH4uEW65aDKrbr+QiWPS+MLj5XzxiY3UNbcP6biNbZ3cumIDqUkJ/OeNpY4toOTvz+/v3x9J7pqmYT2oL5C/sb+iOrKN4arKirVuZo0fzdwo/gIPNDU/kx9dN5tN++v4zos7h3y8bQfqWfyLt/nVX/ex5NwiXv7yxcyfnBuGSMNHRJian8G/XjKFp249ny/O7Xsi0XBxtE1DVVcBq3rcd1/A7VbgumjF43SvhGBMzc/kj/92Af+55n1+9vp7rN1Xw91XzaDAd7lKAVVQvBUy7226R12r7z/+/Svf3c++Y42suOW8sFfpQ1GYk8aHZxSwcl0lt186LaLJy13TzLkl0R0X4pQx6clkpCRGfLbbDe5adh1u4PufPMvRqcuvOmscyy+ezK/f3MdpWamcMS4TwRePgEB3fN7bdO8X334E3tl3nIdX72VMejL//U/nsjAKl03DIVwdSfozLBrCR5rEBBe3XzqND8/I56tPbeYrT24e0vG+eeUMLpgSmcbDUNw0v4jXdh7h5e2HWTx7fETO0dbZxcH6Fopzo9tu4xQRoTg3LeLrajy61k1maiKL50TmdQvFN66Yztaqen70yu4hHeeaOeO5f/EsstNi66qD0yxpxLEzJ2Txwhc/xJaqero86v89hYic/NWE9xcUPX5N+X9tZaQkxszCSBdPG0vRmDRWrHVHLGnsP96CKpTkjYzLU+AdGb7j0ImIHb+6sY1VWw9x43nFMbGyXmKCi//553nsPHQCj36wln2yRVS7a+H47lfV7u3RqUnMHO/8dD6xyPlX2AxJSmIC55YMj7UFXC7hpvlFfHfVLnYfbmD6aeGfGmLbgXpgeM9u21NxbhqvbD9MbVN7REbdP7V+Px1d6kgDeF+SElycXehM28pwZ7PcmphyXZl3VtQVa91hP3Z1Yxvf+fNOpuVncGYMrpsdKR+fPR4RuPMPW8I+o3CXR3lsbSXnT86NmRqriSxLGiam5KQn87Gzx31gLp5wUFW+/vRmTrR28PNPz3Wsl5gTzhg3mjsXzeAvO47wWJh7p63ZfZQDdS0sOz92ahkmskbOJ8fEjZvmF9PY1slzG3uO9Ry8//5bBat3H+Oeq86Iianno+1zH5rExaeP5dsv7mDPkYaBHxCkFWvd5GemcPnM4bnWujmVJQ0Tc+ZOzGbW+NGsCNMCTdsP1vP9l3Zx6Yx8PjNCfxG7XMKD151NZmoitz+xMSyr+VXWNLNmzzGWzCuKyrQ6JjbYK21ijoiwbH4xuw43sMFdO6RjNbd3cvsTG8lOS+KH157t6BgCp+VnpvKj62az63AD31s19MFvj73rxiUy6HXjTXyypGFi0uI548lMTeTRITaIf/vFneyrbuIn188hNyPyo2Vj3cLp+XzuQ5P4/d/dvL7zyKCP09rRxdPrq7jsjHxHB4Wa6LOkYWJSWnIinyot5KWth6lubBv4Ab14edshnni3kuUXT+bCac4PXowVd145nZnjRvP1Z7Zw9ETroI7x0rZDHG9qZ9n8kvAGZ2KeJQ0Ts26aX0x7l4en1u8P+bE1LR7u/MNWzi7M4quXT49AdPErJTGBny2dS0t7F3c8tRmPJ/R2o0f/7mZyXjoXTImtuZhM5FnSMDFran4G50/O5bG1lXSF8MXW5VF+vaWNzi4PP1sysrrXBmtqfgb/5+MzeXtvNb9+a19Ij91+sJ7yyjo+fV5RxFd2NLHHPk0mpi07v5gDdS38dU/wU1E/vHovu2s9PHD1mZQ4sMZ6vLjh3IlcddZpPPjKbvbVB9+basXaSlKTXFxXZg3gI5ElDRPTLp9ZQH5mCo/+PbgG8Q3u4zz0+nvMH5fAJ0snRDi6+CYifO8TZ1MwOpX/2twW1GDKE60dPLfxAItnjycrbfivR2JOZUnDxLSkBBdL5hWxZs8x9h/vf6bW+pYObn9iE+OzU7l5VsqI7l4brKy0JP5jyRyONSv3Pb9twPLPlh+gpaMrpuaZMtFlScPEvKXzJuIS6XcKDFXlnme3cvhEKw8tmcuoREsYwTq3ZAyLpyTxx/IDPL+p71H4qsqja93MLsyyyQBHMEsaJuaNyxrF5WcU8NT6/bR19n7t/ekNVby45RB3XH46pUUjY4GlcFo8JYlzinO459ltVPax9sbafcfZe7TRahkjnCUNExduml/M8aZ2Xtp6+JR9+441cv8L25k/eQy3XjLFgejiX4JL+I8lcxCB21dupKPr1LWmV7zjJmtUEh+P0FonJj5Y0jBx4YIpuUzOSz9lhHh7p4cvrdxEcqKLn94wJy6W7I1VhTlpfO+TZ7Fpfx0PvfbeB/YdPdHKK9sOc11ZIalJCQ5FaGKBJQ0TF1wu4dPnFbHBXcuOgydXoXvwL7vZeqCeH3zqbJvOIgw+dvZ4rj+nkIfX7OXv79d0379y3X46PcqNdmlqxLOkYeLGdWUTSU1yseIdb23jzT3H+PWb+7jxvCKumHWaw9ENH/cvnsWk3HS+8uQmapva6ezy8MS7lVw0LY9JNu5lxLOkYeJGVloSi2eP57mNB6iobuKOpzYzLT+Db310ptOhDStpyYn8bOlcaprauPMPW3ht51EO1bdaA7gBLGmYOHPT/GKa27v45C//t3sVvlHJdo093M6ckNW92t/dz25lXFYql87IdzosEwMcSRoiMkZEXhWR93z/9tpHUkReFpE6EXkx2jGa2HR2YTazC7M43tQ+Ylfhixb/an/Hm9pZOq+IRFtoyeBcTeMu4HVVnQa87tvuzY+AZVGLysSFb19zJncumjFiV+GLFpdL+Mn1s/n8gincfH6J0+GYGJHo0HmvBhb4bv8eWAPc2bOQqr4uIgt63m9GtrMLs21EcpTkZaTwjUUznA7DxBAJxxrMIZ9UpE5Vs323Baj1b/dSdgHwNVX9WD/HWw4sBygoKChbuXLloGNrbGwkIyNj0I+PFIsrNBZXaCyu0AzHuBYuXLhBVc8ZsKCqRuQPeA3Y1svf1UBdj7K1/RxnAfBisOctKyvToVi9evWQHh8pFldoLK7QWFyhGY5xAes1iO/YiF2eUtXL+tonIkdEZJyqHhKRcUDwiyUYY4xxjFMN4S8AN/tu3ww871AcxhhjQuBU0vg+cLmIvAdc5ttGRM4RkUf8hUTkLeBp4FIRqRKRKxyJ1hhjDOBQ7ylVrQEu7eX+9cAtAdsXRTMuY4wx/bPROsYYY4JmScMYY0zQHBmnEUkicgxwD1iwb3lAdZjCCSeLKzQWV2gsrtAMx7iKVXXsQIWGXdIYKhFZr8EMcIkyiys0FldoLK7QjOS47PKUMcaYoFnSMMYYEzRLGqf6tdMB9MHiCo3FFRqLKzQjNi5r0zDGGBM0q2kYY4wJmiUNY4wxQRuRSUNEFonIbhHZKyKnrBooIiki8qRv/zsiUhKFmCaKyGoR2SEi20XkS72UWSAi9SKyyfd3X6TjCjh3hYhs9Z13fS/7RUR+5nvOtohIaRRimh7wXGwSkRMi8uUeZaLynInIb0XkqIhsC7gv2GWNb/aVeU9Ebu6tTJjj+pGI7PK9Ts+KSF9r2fT7mkcgrvtF5EDAa3VVH4/t9/MbgbieDIipQkQ29fHYSD5fvX4/OPIeC2b+9OH0ByQA7wOTgWRgMzCzR5nPA//lu70EeDIKcY0DSn23M4E9vcS1gBDWFglzfBVAXj/7rwJeAgSYD7zjwOt6GO8Apag/Z8DFQCmwLeC+HwJ3+W7fBfygl8eNAfb5/s3x3c6JcFwfARJ9t3/QW1zBvOYRiOt+vAuuDfQ69/v5DXdcPfb/GLjPgeer1+8HJ95jI7GmMQ/Yq6r7VLUdWIl3YahAV+NdhhbgGbyz7Eokg1LVQ6pa7rvdAOwEJkTynGF2NfA/6rUWyPatlRItlwLvq+pQZgMYNFV9Ezje4+7A99HvgWt6eegVwKuqelxVa4FXgUWRjEtV/6Kqnb7NtUBhuM43lLiCFMznNyJx+b4DrgeeCNf5gtXP90PU32MjMWlMAPYHbFdx6pdzdxnfh6seyI1KdIDvcthc4J1edp8vIptF5CURmRWtmAAF/iIiG8S7vG5PwTyvkbSEvj/MTj1nBap6yHf7MFDQSxmnn7fP4a0h9mag1zwSbvNdNvttH5danHy+LgKOqOp7feyPyvPV4/sh6u+xkZg0YpqIZAB/AL6sqid67C7He/llNvBz4LkohnahqpYCVwJfEJGLo3jufolIMrAY79orPTn5nHVT73WCmOrfLiL3AJ3AY30UifZr/ktgCjAHOIT3UlAsWUr/tYyIP1/9fT9E6z02EpPGAWBiwHah775ey4hIIpAF1EQ6MBFJwvuGeExV/9hzv6qeUNVG3+1VQJKI5EU6Lt/5Dvj+PQo8i/cyQaBgntdIuRIoV9UjPXc4+ZwBR/yX6KTvZY0ded5E5LPAx4AbfV82pwjiNQ8rVT2iql2q6gH+Xx/nc+r5SgQ+CTzZV5lIP199fD9E/T02EpPGOmCaiEzy/UJdgnf52UCBy9FeC7zR1wcrXHzXS38D7FTVn/RR5jR/24qIzMP7+kUjmaWLSKb/Nt6G1G09ir0AfEa85gP1AdXmSOvzF6BTz5lPMMsavwJ8RERyfJdjPuK7L2JEZBHwDWCxqjb3USaY1zzccQW2gX2ij/MF8/mNhMuAXapa1dvOSD9f/Xw/RP89FomW/lj/w9vTZw/eXhj3+O57AO+HCCAV76WOvcC7wOQoxHQh3qrlFmCT7+8q4FbgVl+Z24DteHuMrAUuiNLzNdl3zs2+8/ufs8DYBHjY95xuBc6JUmzpeJNAVsB9UX/O8CatQ0AH3mvG/4y3Hex14D3gNWCMr+w5wCMBj/2c7722F/inKMS1F+81bv/7zN9TcDywqr/XPMJxPep772zB+2U4rmdcvu1TPr+RjMt3/+/876mAstF8vvr6foj6e8ymETHGGBO0kXh5yhhjzCBZ0jDGGBM0SxrGGGOCZknDGGNM0CxpGGOMCVqi0wEYM1yIiL/7I8BpQBdwzLc9T71zJRkT16zLrTERICL3A42q+qDTsRgTTnZ5ypgoEJE/+Say2y4it/juSxSROvGuQ7Ldtx5C1CbGNGYwLGkYEx03q2oZcC5wR8AMrlnA31R1FvB34F6nAjQmGJY0jImOr4jIZryJoRDvbK7gnWXWPzvvCrzTRRgTs6wh3JgIE5HL8K4IN19VW0Tkbbzzm/XGGhlNTLOahjGRlwUc9yWMWXgvUfn5p9wG+DTwdrSDMyYUljSMibw/A2kisgP4Dh9ckbEeuEhEtuO9NPUdB+IzJmjW5dYYh/gW9qlW1WynYzEmWFbTMMYYEzSraRhjjAma1TSMMcYEzZKGMcaYoFnSMMYYEzRLGsYYY4JmScMYY0zQ/j/mQGB47Kj1IgAAAABJRU5ErkJggg==\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEWCAYAAACJ0YulAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvqOYd8AAAIABJREFUeJzt3Xl8HfV57/HPc452S5YX2fIir9gGjFcsL0AgJhDikAan2YBc0qQhkLaXbG3pzW3SNCXdktzQNhduG7InbTCENInbQIAQFCeAjW3wgnd53+XdlrxoOc/9Y8bioNjWsazR6Jzzfb9e8zozc2Z5fpJ9vpr5nZkxd0dERAQgEXcBIiLSeygURESknUJBRETaKRRERKSdQkFERNopFEREpJ1CQURE2ikUJBJmts3MTplZY9owLO664mBmbmZN4c9gt5k9aGbJuOsSOReFgkTpne5enjbs6biAmRXEUVgMprp7OfBm4HbgIzHXI3JOCgXpUWY2OvzL+W4z2wH8Kpw/x8xeNLOjZrbSzOamrTPGzH5tZifM7Fkze8jM/j18b66Z7eqwj21mdnM4njCzz5jZZjM7ZGaPm9mADrV8yMx2mNlBM/ts2naSZvaX4bonzGy5mY0ws4fN7Ksd9rnQzD7dWfvdvR54AZiWtm6lmX3LzPaGRxJ/e/ZIwszGhW0/Ftb3WNp6bmafMLMt4XtfMbNEWrs/Z2bbzazBzL5vZpUZtnuWmS0zs+Nmtt/MHkx777y/J8kR7q5BQ7cPwDbg5nPMHw048H2gD1AKDAcOAbcS/KHy1nB6ULjOS8CDQDFwA3AC+PfwvbnArvPtG/gksBioCdf/OvBoh1q+EdYxFTgDXBm+fz+wGrgcsPD9gcAsYA+QCJerAk4C1ef5WTgwLhy/AtgLfDrt/Z+EdfUBBgMvAx8L33sU+Gz4cykB3tRhu88DA4CRwEbgo+F7HwHqgbFAOfCfwA8ybPdLwAfD8XJgTjh+wd+ThtwYYi9AQ24O4QdzI3A0HH4azj/7gTQ2bdn/dfYDK23e08CHwg+7VqBP2ns/vIhQWAfclPbeUKAFKEirpSbt/ZeBO8LxDcD887RvHfDWcPw+4MkL/CwcOA40heOPAsXhe9XhB3Jp2vJ3As+H498HHkmvscN256VN/wnwXDj+HPAnae9dfhHtXgT8DVDVYX/n/T3F/e9NQ/cNOn0kUXqXu/cLh3d1eG9n2vgo4H3hKYmjZnYUeBPBB/gw4Ii7N6Utv/0iahgF/CRtu+uANoIP47P2pY2fJPjrGGAEsPk82/0ecFc4fhfwg07quDrc7u3AbIKjgrP1FQJ702r8OsERA8BfEBylvGxma8ysY19E+s9xO8HPi/B1e4f3Csis3XcDE4D1ZrbUzH4vrdbz/Z4kR+RLJ5/0Pum3591J8BfoPR0XMrNRQH8z65MWDCPT1m8CytKWTwKDOmz7I+7+wjm2PbqTGncClwGvneO9fwdeM7OpwJXATzvZFh78af24mc0HPg98KtzHGYK/ylvPsc4+4J6w3jcBvzSzRR70TUAQXGvC8ZEEp7UIX0elbersEdd+glNpF6pzE3Bn2D/xbuAJMxvIBX5Pkjt0pCC9wb8D7zSzt4WduyVhB3KNu28HlgF/Y2ZF4QfjO9PW3QiUmNk7zKwQ+BxB38FZ/wb8XRgumNmg8EM5E98Evmhm4y0wJfxwxN13AUsJjhB+7O6nLqK9/wjcY2ZD3H0v8AzwVTPrG3YQX2Zmbw7rfZ+Znf0QP0IQhqm0bd1vZv3NbARB/8nZjuhHgU+HnfTlwN8Dj50reDoys7vMbJC7pwhO/RHu87y/p4tou/RyCgWJnbvvBOYDfwkcIPiL9H5e//f5AYJTLoeBvyY4z3523WME59K/CewmOHJI/zbSvwALgWfM7ARBp/PsDEt7EHic4EP7OPAtgo7Zs74HTKbzU0dv4O6rCc7b3x/O+gOgCFhL8MH/BK+fkpkJLDGzxrAdn3T3LWmb+xmwHFgB/DysEeDbYV2LgK3AaeDjGZY4D1gT7vNfCPoaTmXwe5IcYMERrUj2MLMvEHyb567Olo24jhsI/noe5TH8RzIzB8annUoSuWRKeJEuCE9VfRL4ZhyBIBIVhYLIRTKzKwnOtQ8F/jnmckS6lU4fiYhIOx0piIhIu6y7TqGqqspHjx7dpXWbmpro06dP5wvmELU5P6jN+eFS2rx8+fKD7j6os+WyLhRGjx7NsmXLurRuXV0dc+fO7d6Cejm1OT+ozfnhUtpsZhndCUCnj0REpJ1CQURE2ikURESknUJBRETaKRRERKRdZKFgZt8OHwN4rtsOE9518mtmVm9mq8zs6qhqERGRzET5ldTvAg+RdkfLDt4OjA+H2cC/kvndK+PhDs1NcPJQOBwOXw8Gr6ePQaoNPBUODslCKCyDorLwtRz6DIQ+g6CsKnwdAIlk/G07fQyaDkLTgaBNTQfSpg9BWwtYAsyC14JSKK6Akr7haz8or4aKaigfAn2q4m+XiFyUyELB3Rd18hCT+cD3w5uJLTazfmY2NLy/fLdbuu0w/7nxDK+eWUMy1UrSm0mmWihqa6Kk9fjvDKUtRyhtPUZpy1FKW48Gry1HKfDmbq/NMU4V9ONkYX9OFg3gZOEAmgoHcLJwYDCvcAAniwbQVDiQU4X9aUsUdbpN81ZKWk9waPtW9h48ELblCGUtRyhrOUxp61HKWg6H00cobTlCsvNb7V+UFElOFvanqajq9TaF7Xt96EdropSWZAmtiWJaE8WAYd5GwttIeGv7a9JbSaRagldvIektJFLB/JZkCacKB9LaZxDbdrZwfOUe3nLFYMqLs+5SHJFYxfk/ZjhvfJTgrnDe74SCmd0L3AtQXV1NXV3dRe9s8JIH+OeTr5DYc2n3ejrthRyiL4e9giNewWGC1yNewXH60EISx0iRwIFCWinjDCU0U2ZnKOckA+wEA+04AzjBQDvGAGukrPUIZa1H4NSWTms444WcIRiaKaTZkxRaG4W0UkgrRbRSbqdfX6Ehs7ad8FIOeV8OU8Ehr+SQ9+UQfTnofTniFbRQgEF760qsmQpOUm6nqOAU/ayRwRxlkB1lsB1hgDVS3nKQ8paDXfpZd9X21GDWbB7Nl3wyBwbOZMrIwUwcmNtHLI2NjV36f5HN1OZoZMWfUe7+CMHDy6mtrfUuXdG3519howMGBcWQLA5O7RRXQGk/KO3/xqGsCsoGBqd2yga2DyVFZQwnSK9u09b6+umaxoa01wZoPBC8Nh0Ixw9QTAvFtLy+vp1rowal/ThJKWUDhgVt6jPo9VNX7UNVeBqriorCUioInureLVqbg9pP7A/b0qFNjeFpqZaT0HoaWk4H42aQKAgGSwSvyaLg95UsfH08cXa8AM404icPwol9jKKBUTRwKy+TOvZtnllRy7oJH+EP77idksLcDAdd3ZsfeqLNcYbCboLny55VE86Lxvt/QN1vfsPcG2+KbBddliyAiiHB0JlUKvgAbTsDreHQ1pL2YRl+YBaVQyLBy3H+xykogsqaYOgBBtDWwtKnfsjMoQ4bfoHV/5J5LIXNS/nNV/6DsR/4Z4aPHt8j9YhkozhDYSFwn5ktIOhgPhZVfwIQfEBZDvyVmEgEndavP6te0iULaSofBTPmwowPYyf2w8uP0PbiQ1zf/FuOf/fN7LjlXxh57fvirlSkV4ryK6mPAi8Bl5vZLjO728z+yMz+KFzkSWALUA98g+A5uyLdq6Iabvorkp9YzqnRN9OXJkY+81FOPPuPwTeuROQNovz20Z2dvO/A/4xq/yJvUFlD6YeeYP8vvsKgxX9PxQv/QHPrCYrm/W3QhyEigK5olnxiRvXb/4J11/0TLZ6kaMlDpBZ9Ne6qRHoVhYLknatu+UNenP4lUm4knv8irHo87pJEeg2FguSlG+Z/lMeqgrOXqYWfgIb1MVck0jsoFCQvmRk3fvBz/MxvINF6Cv/Rh6D5ZNxlicROoSB5a0i/Uk7f8hU2pYZjB9ZD3T/EXZJI7BQKktfed83lfGfQX9CG4S89BHtWxF2SSKwUCpLXEgnjnjvfx/fb5mGegv/6ZHDVuEieUihI3htT1YcdUz/NPh8Ae1fA6h/FXZJIbBQKIsDH3jqVf0q9P5h47gFoORVvQSIxUSiIAEMqS6iYdRdrU6Pg+C5Y8vW4SxKJhUJBJPRHN07gQf9AMPHi/9VXVCUvKRREQlXlxUy4dj4rUmOD51u88r24SxLpcQoFkTQfe/M4vsm7g4kXvhY8r0IkjygURNJUlhUy8Or5rPeRcGKPvokkeUehINLBh64byzdabg0mlnxdz12QvKJQEOlg7KByjo97J0foC/tWwc6X4y5JpMcoFETO4YNvupwfts4NJl7W11MlfygURM7h+vFVvNBvPm0k8LU/gxP74y5JpEcoFETOwcx4x/Uzea5tOpZqhVWPxV2SSI9QKIicx7un1/DfybcEEyt+qA5nyQsKBZHzKC1K0nfKOzjkfeHAOtjzatwliUROoSByAe+pHc1P264LJlb8R7zFiPQAhYLIBUwb0Y/FfecFE6uf0BXOkvMUCiIXYGbMmH19cPfU00dh86/iLkkkUgoFkU78/vTh/Dw1J5hY85N4ixGJmEJBpBPVfUs4OPLtAPj6J6HldMwViURHoSCSgRvmzGF1ajTWfAI2Pxd3OSKRUSiIZODmiYN5LhF+C0mnkCSHKRREMlBckOTU+HcC4Bue0jOcJWcpFEQydE3tjPAUUiNs+XXc5YhEQqEgkqHrxlXxm8SsYGLjU/EWIxKRSEPBzOaZ2QYzqzezz5zj/ZFm9ryZvWpmq8zs1ijrEbkUhckEzZfdAoBveBpSqZgrEul+kYWCmSWBh4G3AxOBO81sYofFPgc87u7TgTuA/xdVPSLdYfqsN7PHB2CNe2HvirjLEel2UR4pzALq3X2LuzcDC4D5HZZxoG84XgnsibAekUt27bgqfmszgomNv4i3GJEImEd0O2Azey8wz90/Gk5/EJjt7velLTMUeAboD/QBbnb35efY1r3AvQDV1dUzFixY0KWaGhsbKS8v79K62Upt7n7LXlnMnx//B473GcsrM/8psv1cDP2e88OltPnGG29c7u61nS1X0KWtd587ge+6+1fN7BrgB2Y2yd3fcLLW3R8BHgGora31uXPndmlndXV1dHXdbKU2d79E9TiaFjxI36YtzJ0+DiprIttXpvR7zg890eYoTx/tBkakTdeE89LdDTwO4O4vASVAVYQ1iVyyayYMY4lNCSY2PRNvMSLdLMpQWAqMN7MxZlZE0JG8sMMyO4CbAMzsSoJQOBBhTSKXrDCZ4NCQGwBI1T8fczUi3SuyUHD3VuA+4GlgHcG3jNaY2QNmdlu42J8B95jZSuBR4MMeVSeHSDeqmhbcIC+1pQ7aWuMtRqQbRdqn4O5PAk92mPf5tPG1wHVR1iAShdpp09j21BBGN++DPa/AiFlxlyTSLXRFs0gXVJQUsqkiDAI9eEdyiEJBpIuS424C4PT6Z2OuRKT7KBREuujyOW+nxZMU7X8VTh2NuxyRbqFQEOmi4UOqWVdwBQlvg62L4i5HpFsoFEQuwbFh1wNwZsMvY65EpHsoFEQuweCpbwOgpb4u3kJEuolCQeQSjJ92PU2UUN60HY7rfo6S/RQKIpcgUVDIjvKpAKS2/ibmakQunUJB5BKlRr0JgMNrdL2CZD+Fgsglqpke9CsU7PhtzJWIXDqFgsglqhwzgyYro9/pXXBsV9zliFwShYLIpUoWsK9yOgAnN/465mJELo1CQaQbFF4W3Er70JrnYq5E5NIoFES6wdDptwDQZ8+LMVcicmkUCiLdoHDYVE4m+jCgeS9+ZHvc5Yh0mUJBpDskkhweOAOAhtX6aqpkL4WCSDcpnxD0KxzZoIvYJHspFES6Sb/Lg1Aob1gecyUiXadQEOkuw6bRYkXUtGyjpfFQ3NWIdIlCQaS7FBRzYsBkALavUL+CZCeFgkg3Kh13HQBH1+uWF5KdFAoi3aj0suDmeOUNS2OuRKRrFAoi3WnELADGnNnI6VMnYy5G5OIpFES6U2l/GivHU2wtbFqhr6ZK9lEoiHSzojHXAnBk/aKYKxG5eAoFkW5WNCbobC7dq34FyT4KBZHuNnIOAJedWUPj6eaYixG5OAoFke7WbyRnSqsZYI2sWbUs7mpELopCQaS7mVEwOjhaOLhW/QqSXRQKIhFIjpwNQMFe3QdJskukoWBm88xsg5nVm9lnzrPM+81srZmtMbMfRlmPSI+pmQnA6FPrOH66JeZiRDIXWSiYWRJ4GHg7MBG408wmdlhmPPC/gevc/SrgU1HVI9Kjhk4llShivO1ixaYdcVcjkrEojxRmAfXuvsXdm4EFwPwOy9wDPOzuRwDcvSHCekR6TkExDJlMwpw9a1+IuxqRjBVEuO3hwM606V3A7A7LTAAwsxeAJPAFd/9Fxw2Z2b3AvQDV1dXU1dV1qaDGxsYur5ut1Ob4jLNh1LCcU5t+Q13dyEj31Vva3JPU5mhEGQqZ7n88MBeoARaZ2WR3P5q+kLs/AjwCUFtb63Pnzu3Szurq6ujqutlKbY7RwIPw4/9idPMmZl97PaVFych21Wva3IPU5mhEefpoNzAibbomnJduF7DQ3VvcfSuwkSAkRLJf2Nk81Tbx6o7DMRcjkpkoQ2EpMN7MxphZEXAHsLDDMj8lOErAzKoITidtibAmkZ7TbySpPoMZYI1sWLcq7mpEMhJZKLh7K3Af8DSwDnjc3deY2QNmdlu42NPAITNbCzwP3O/ueo6h5AYzEuHRwskti2MuRiQzkfYpuPuTwJMd5n0+bdyBPw0HkdxTUwsbfk6/Qytobk1RVKDrRaV3079QkSiFD92ZwkZW7z4WczEinVMoiERp2HTcElxpO3ilvuP3LER6H4WCSJSK+mDVV1FgKQ5sXBJ3NSKdUiiIRC3sbC7e9wptKY+5GJELUyiIRC0MhYmpDWxqOBFzMSIXlvG3j8xsODAqfR13183iRToThsLViU08u/UwVwzpG3NBIueXUSiY2ZeA24G1QFs42wGFgkhnBlyGl/Sj+vRRNm/eANeMjrsikfPK9EjhXcDl7n4mymJEclIigdXUQv0vad3xMvC2uCsSOa9M+xS2AIVRFiKS04bXAjDi5Doajp+OuRiR88v0SOEksMLMngPajxbc/RORVCWSa8J+hemJTSzbfoRbJw+NuSCRc8s0FBbyuzezE5FMDb8agMm2la9uPaBQkF4ro1Bw9+9FXYhITisbAAPHUXKonkNbXgGmxF2RyDldMBTM7HF3f7+ZrSb4ttEbuLv+ZYtkangtHKqn4uAKTjW3RfrQHZGu6uxI4ZPh6+9FXYhIzquphVULmGKbWLnrKHPGDoy7IpHfccFQcPe94ev2nilHJIfVBN9Ammab+cX2IwoF6ZUy+kqqmc0xs6Vm1mhmzWbWZmbHoy5OJKdUT4KCEi5L7GXd5m1xVyNyTplep/AQcCewCSgFPgo8HFVRIjkpWQhDpwHQtms5Kd0cT3qhjG+I5+71QNLd29z9O8C86MoSyVHhKaQJLRvYcrAp5mJEflfGF6+ZWRGw0sy+DOxFd1gVuXhhKExP1PPqjiOMG1wec0Eib5TpB/sHw2X/J9AE1ADviaookZwV3u5iejIIBZHeprPrFOYDNe7+cDj9a2AwwTULLwH1kVcokksqa6B8CJWN+2jYthZdxCa9TWdHCn/BG29vUQzMAOYCfxxRTSK5y6z9FFLfQytoPNMac0Eib9RZKBS5+8606d+6+2F33wH0ibAukdzVfr1CPat2HY25GJE36iwU+qdPuPt9aZODur8ckTwQ9itMS9Tz6g6FgvQunYXCEjO7p+NMM/sY8HI0JYnkuGHTwRJMTOzgtW37465G5A06+0rqp4GfmtkHgFfCeTMI+hbeFWVhIjmruBwGT6Rg/2uc2fkK7tdiZnFXJQJ0fu+jBuBaM3sLcFU4++fu/qvIKxPJZcNnwP7XGHNmHTsPn2LkwLK4KxIBMn+ewq8ABYFId6mZCa98L7iIbecRhYL0GroqWSQO4TeQrk5sVmez9CoKBZE4VF0OxX0ZZgfZtm1z3NWItFMoiMQhkQi+hQSU7H+V0y1tMRckEog0FMxsnpltMLN6M/vMBZZ7j5m5mdVGWY9IrxKeQppqm1izR48nkd4hslAwsyTBMxfeDkwE7jSziedYroLgsZ9LoqpFpFeqmQkET2JbuVP9CtI7RHmkMAuod/ct7t4MLADmn2O5LwJfAk5HWItI7xNe2Tw1uYWVOw7FXIxIINPnKXTFcCD9vkm7gNnpC5jZ1cAId/+5md1/vg2Z2b3AvQDV1dXU1dV1qaDGxsYur5ut1ObebXZJNWWn93Ng01Lq6hq7vJ1sanN3UZujEWUoXJCZJYAHgQ93tqy7PwI8AlBbW+tz587t0j7r6uro6rrZSm3u5Q6+CV77MSOb65k682769ynq0mayqs3dRG2ORpSnj3YDI9Kma8J5Z1UAk4A6M9sGzAEWqrNZ8krYrzDd6lmpO6ZKLxBlKCwFxpvZmPBRnneQ9mwGdz/m7lXuPtrdRwOLgdvcfVmENYn0Lml3TF2hzmbpBSILBXdvBe4DngbWAY+7+xoze8DMbotqvyJZZegUSBYxPrGbDdt3d768SMQi7VNw9yeBJzvM+/x5lp0bZS0ivVJBMQyZTGL3cnz3K7jP1R1TJVa6olkkbuEppMvOrGfXkVMxFyP5TqEgErezF7GpX0F6AYWCSNxqZgBwdaKeFTuOxFyM5DuFgkjc+o+BsoEMtOPs3b4h7mokzykUROJm1t6vUNLwKi1tqZgLknymUBDpDcJ+hUmpjWzcfyLmYiSfKRREeoOwX2Faop5Vu47FXIzkM4WCSG8wfAaOMSmxjTU7GuKuRvKYQkGkNyipxKomUEQrTTtWxF2N5DGFgkhvET6Jrf/hVXo8p8RGoSDSW4ShMMU2sXavHs8p8VAoiPQW4ddSp1s9q3Rls8REoSDSWwyeiBeWMSrRwObt2+OuRvKUQkGkt0gWYMOmA5DaqceKSDwUCiK9SdivMOTEahrPtMZcjOQjhYJIbxL2K0y1zazWRWwSA4WCSG+Sdhvt1bsOx1yM5COFgkhv0ncoVI6gr53iwJaVcVcjeUihINLbjLwGgNI9i2MuRPKRQkGktxl1LQATTq/mcFNzzMVIvlEoiPQ2YSjMTKxn9S5dxCY9S6Eg0ttUTSBVOpBqO8rOzWvirkbyjEJBpLcxIzE6OFpo2/pCzMVIvlEoiPRGI4NQqDq8POZCJN8oFER6o7BfYVLLaxxqPBNzMZJPFAoivdGQybQWljMq0cCGTRvirkbyiEJBpDdKJGHEbACOr18UczGSTxQKIr1UwZjrACjeuyTmSiSfKBREeqtRQSiMPKFnNkvPUSiI9FbDptOaKOYy38HB/bvjrkbyRKShYGbzzGyDmdWb2WfO8f6fmtlaM1tlZs+Z2ago6xHJKgXFNFUHt9Lev/KZmIuRfBFZKJhZEngYeDswEbjTzCZ2WOxVoNbdpwBPAF+Oqh6RbFRy+VsASG3+dcyVSL6I8khhFlDv7lvcvRlYAMxPX8Ddn3f3k+HkYqAmwnpEsk7x+CAUhhzSHVOlZxREuO3hwM606V3A7Assfzfw1LneMLN7gXsBqqurqaur61JBjY2NXV43W6nNWc7bqLU+DGrdy+KnFnC6dMg5F8upNmdIbY5GlKGQMTO7C6gF3nyu9939EeARgNraWp87d26X9lNXV0dX181WanP227p2DuUHnmNixQn6vumOcy6Ta23OhNocjShPH+0GRqRN14Tz3sDMbgY+C9zm7rqeX6QDGxv8rXRqw69irkTyQZShsBQYb2ZjzKwIuANYmL6AmU0Hvk4QCA0R1iKStQZPvQWAir0vQioVczWS6yILBXdvBe4DngbWAY+7+xoze8DMbgsX+wpQDvzIzFaY2cLzbE4kb5UNvYIGG0hZ61Fo0PMVJFqR9im4+5PAkx3mfT5t/OYo9y+SE8zYVlHL4ONPw+bnYcjkuCuSHKYrmkWyQNOIGwBo3vBszJVIrlMoiGSBiqvmkXKjYNdLcOZE3OVIDlMoiGSBy8eO4lUfRyLVAlt0dbNER6EgkgUqSgpZUTIrmNj0dLzFSE5TKIhkiUND5wYjm54F91hrkdylUBDJEv3GXs0+7w8n9sK+1XGXIzlKoSCSJSbV9OP5tmnBxCbdSluioVAQyRKThldSlwpDYaP6FSQaCgWRLNG3pJAd/WbRbEWwaykc3xt3SZKDFAoiWWRszRAW2zTAYf1/x12O5CCFgkgWmTy8kp+cnhFMrP1ZvMVITlIoiGSRycMreS51NalEIWx/AZoOxl2S5BiFgkgWmTSskuP0YWe/WeApnUKSbqdQEMkilWWFjBhQym8Krglm6BSSdDOFgkiWmTSskscbp4Alg/sgNR6IuyTJIQoFkSwzaXglqw4X0DL2JvA2eO2JuEuSHKJQEMkyk4ZXArBl+DuDGSsfjbEayTUKBZEsM2lYXwBesJlQXAl7V9KncXvMVUmuUCiIZJmB5cUMqyxhxb7TMOn3Aaje/3zMVUmuUCiIZKGrhlfy2p5jMPVOAKr310Frc7xFSU5QKIhkoUnDKtl6sInGwTNg0BUUNx+B9f8Vd1mSAxQKIllock1f3GHt3hMw655g5svfjLcoyQkKBZEsNGlY8A2k13Yfgym305oshR0vwr7XYq5Msp1CQSQLDe5bwqCK4qBfobiCfUPeEryx5N/iLUyynkJBJEtNHl4ZHCkAu4e/AzBYuQCO7Y63MMlqCgWRLDVpWF/qGxo51dzGqbLhMOndkGqBF78Wd2mSxRQKIlnqquGVpBzW7TsezLj+z4LX5d+FE/tjq0uym0JBJEtNqQk6mxdvORTMqL4Krvg9aD0Ni74cY2WSzQriLkBEumZoZSnXXjaQ77ywjb+bkwxmvuVzsOEpWPYdmHkPDL7i3Csf3QGv/SdsXQRHtsLp41BQDJU1MPhKGHMDjL0Rygb0XIOkV1AoiGSxj79lPHd+YzG/3lXELRB8oM/4MCz7Fjz9l3DXj8Hs9RWO74VnPw+v/Ti4w2pHx3fDziXBKSgYgzMPAAAL9klEQVRLwvhbYPpdMOFtkCzsmUZJrBQKIllsztgBzBo9gCe3HOHzrW0UFyThxr+E1U/A5udg9Y9gyvuDhdf+DH72cThzDBIFcNV74PJbYchkKO0PLafgyDbY82qw7vYXYeNTwdB3OMz+I5jxISipvHBR7nBiH+xdCQ1rgiBqanj9NhwFRdBnEJRXw4CxMHgiDBwXzJfYRRoKZjYP+BcgCXzT3f+xw/vFwPeBGcAh4HZ33xZlTSK5xMz4xE3juetbS/jRsl3cNWcU9KmCt/0tLPw4PPnnMPp6WPUY/PKvg5UmzINbvwL9Rv7uBvuPgrFvhjd9Knh4z6rH4JXvwcGN8Oxfwa+/HATD7I+9cf22VthaBxufDoajF3nX1kQhDJsGI6+BUdfCiNnnP3XlDoe3UHVgMSxaBgc3wYm9cPJQ8MzqtubgUaU4FFUEgVfaDypHwIAx0H90EEKDr4TC0s5ra26Cozvh2C44vivoxG8+EcxvPgmWCI6ikkVQXAHlg18Pvf6jg0BNXET3bcspaGyApgPB0NwErWeg7Qz9jjQCczPfVhdEFgpmlgQeBt4K7AKWmtlCd1+bttjdwBF3H2dmdwBfAm6PqiaRXHTduIFcVpngX+s2c/vMERQmEzD9g7Duv2DTM/Dg2X4Fg7c+ANd+/I2nlM6nfBBcex/M+ROo/2XwVddtv4GXHoLF/y84ypj2P4LTTSsfhca0bzwVV8LQKTBkStBPUT4YCsuC91pOBh/eJ/YGH+gNa4MjlF1Lg+HsV2oHXQkjZwdBUVkDu5fDjsXBcOowkwDWdNKG08eCD/JzsSRUTQiOlIZMhkGXBx/kR7bB/jVBXfvXBH0ul6KgBPqPCY6KzobSgDFQVhWE56F6OLQ5fK0Pwu08qofcBHzi0urprNwItz0LqHf3LQBmtgCYD6SHwnzgC+H4E8BDZmbu7hHWJZJTzIz54wp5cPkpbvrqrykuCP4qHZD6AP9qLzPAjwLw5ZJP8OySqbBkURf2Ugzcz7g+7+a9Z37KDa0vULj+v2H9f7cvsSMxnEUF1/FyQS0bkuNJHU7C4XNtqw8wCLgSCK7ELqtoYmLbeia3rmVS2xquaNtE0YF1cGBd0L/RwSHrzyZGsKtwDDsSNRxIDOKIVXLMKjlDUXvolflJyr2RSj9OdaqBYal9DEvtZWRqJyNSu0me3cfqx8/b8hYK2J8YzAGr4kCiioM2kCYr45SVcoZiAApooZBW+ngT/VNH6e9HGeiHGZbay4DWo8E+DqzL6CfdTAFHrV8wJCppoowWC/YAE3hHRlvpOovq89fM3gvMc/ePhtMfBGa7+31py7wWLrMrnN4cLnOww7buBe4FqK6unrFgwYIu1dTY2Eh5eXmX1s1WanN+OHGikaf2FNFwMvWG+RNaN3L3qW/zbPHN/LLo5m7bX7/UUd7a/Cy1LcvZnhzFc0U3siF5eWZHIBko8BbGtm3lirb1XNm6gQF+mM3JsaxPXs76gitosEG0tbWRLOj637VF3syIth2MSW1ndNs2hqb2MSB1mIbEYLYnR7AjMZLtyZHsTQyj1bq+n1I/yZDUfoa07WNIah/VqQaqUw309WMcSAxiT2IoexND2ZMcxp7EUI5Y//P+HGdXtTJzRNf+bd94443L3b220wXdPZIBeC9BP8LZ6Q8CD3VY5jWgJm16M1B1oe3OmDHDu+r555/v8rrZSm3OD2pzfriUNgPLPIPP7igvXtsNjEibrgnnnXMZMysAKgk6nEVEJAZRhsJSYLyZjTGzIuAOYGGHZRYCHwrH3wv8Kkw0ERGJQWQdze7eamb3AU8TfCX12+6+xsweIDiMWQh8C/iBmdUTdEndEVU9IiLSuUivU3D3J4EnO8z7fNr4aeB9UdYgIiKZ0w3xRESknUJBRETaKRRERKSdQkFERNpFdkVzVMzsAHCRd9tqVwUc7HSp3KI25we1OT9cSptHufugzhbKulC4FGa2zDO5zDuHqM35QW3ODz3RZp0+EhGRdgoFERFpl2+h8EjcBcRAbc4PanN+iLzNedWnICIiF5ZvRwoiInIBCgUREWmXk6FgZvPMbIOZ1ZvZZ87xfrGZPRa+v8TMRvd8ld0rgzb/qZmtNbNVZvacmY2Ko87u1Fmb05Z7j5m5mWX91xczabOZvT/8Xa8xsx/2dI3dLYN/2yPN7HkzezX8931rHHV2FzP7tpk1hE+mPNf7ZmZfC38eq8zs6m4tIJMn8WTTQHCb7s3AWKAIWAlM7LDMnwD/Fo7fATwWd9090OYbgbJw/I/zoc3hchXAImAxUBt33T3wex4PvAr0D6cHx113D7T5EeCPw/GJwLa4677ENt8AXA28dp73bwWeAgyYAyzpzv3n4pHCLKDe3be4ezOwAJjfYZn5wPfC8SeAm8y66eGy8ei0ze7+vLufDCcXEzwJL5tl8nsG+CLwJeB0TxYXkUzafA/wsLsfAXD3hh6usbtl0mYH+objlcCeHqyv27n7IoLny5zPfOD7HlgM9DOzod21/1wMheHAzrTpXeG8cy7j7q3AMWBgj1QXjUzanO5ugr80slmnbQ4Pq0e4+897srAIZfJ7ngBMMLMXzGyxmc3rseqikUmbvwDcZWa7CJ7f8vGeKS02F/v//aJE+pAd6X3M7C6gFnhz3LVEycwSwIPAh2MupacVEJxCmktwNLjIzCa7+9FYq4rWncB33f2rZnYNwdMcJ7l7Ku7CslEuHinsBkakTdeE8865jJkVEBxyHuqR6qKRSZsxs5uBzwK3ufuZHqotKp21uQKYBNSZ2TaCc68Ls7yzOZPf8y5gobu3uPtWYCNBSGSrTNp8N/A4gLu/BJQQ3DguV2X0/72rcjEUlgLjzWyMmRURdCQv7LDMQuBD4fh7gV952IOTpTpts5lNB75OEAjZfp4ZOmmzux9z9yp3H+3uown6UW5z92XxlNstMvm3/VOCowTMrIrgdNKWniyym2XS5h3ATQBmdiVBKBzo0Sp71kLgD8JvIc0Bjrn73u7aeM6dPnL3VjO7D3ia4JsL33b3NWb2ALDM3RcC3yI4xKwn6NC5I76KL12Gbf4KUA78KOxT3+Hut8VW9CXKsM05JcM2Pw3cYmZrgTbgfnfP2qPgDNv8Z8A3zOzTBJ3OH87mP/LM7FGCYK8K+0n+GigEcPd/I+g3uRWoB04Cf9it+8/in52IiHSzXDx9JCIiXaRQEBGRdgoFERFpp1AQEZF2CgUREWmXc19JFTkXM2sDVqfNepe7b4upHJFeS19JlbxgZo3uXn6B9wvC+2CJ5DWdPpK8ZWYfNrOFZvYr4Llw3v1mtjS8T/3fpC37WTPbaGa/NbNHzezPw/l1Z2+dYWZV4S01MLOkmX0lbVsfC+fPDdd5wszWm9l/nL1Dr5nNNLMXzWylmb1sZhVmtsjMpqXV8Vszm9pTPyPJPzp9JPmi1MxWhONb3f33w/GrgSnuftjMbiG4T9AsgnvVLzSzG4AmgqvepxH8n3kFWN7J/u4muP3ATDMrBl4ws2fC96YDVxHc4vkF4Dozexl4DLjd3ZeaWV/gFMHV9x8GPmVmE4ASd195ST8JkQtQKEi+OOXu084x/1l3P3vv+lvC4dVwupwgJCqAn5x9HoWZZXILjVuAKWb23nC6MtxWM/Cyu+8Kt7UCGE1w+/a97r4UwN2Ph+//CPgrM7sf+Ajw3UwbLNIVCgXJd01p4wb8g7t/PX0BM/vUBdZv5fXTsCUdtvVxd3+6w7bmAul3qG3jAv8P3f2kmT1L8GCV9wMzLlCLyCVTn4LI654GPmJm5QBmNtzMBhM8zvNdZlZqZhXAO9PW2cbrH9Tv7bCtPzazwnBbE8yszwX2vQEYamYzw+Urwtu6A3wT+Bqw9OwT1USioiMFkZC7PxPeevmlsO+3EbjL3V8xs8cIng/cQHA757P+D/C4md0LpD/h7ZsEp4VeCTuSDwDvusC+m83sduD/mlkpQX/CzUCjuy83s+PAd7qpqSLnpa+kilwkM/sCwYf1/+mh/Q0D6oAr9DQxiZpOH4n0Ymb2B8AS4LMKBOkJOlIQEZF2OlIQEZF2CgUREWmnUBARkXYKBRERaadQEBGRdv8fr2wleJnRt9oAAAAASUVORK5CYII=\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"import scipy\n",
"import scipy.signal\n",
"\n",
"# Calculate the desired RRC filter response\n",
"b = 0.2\n",
"f = np.linspace(0.0, 1.0, 100)\n",
"r = np.array([rrc(1.0, b, x) for x in f])\n",
"\n",
"# Calculate the FIR taps required for the desired response\n",
"taps = scipy.signal.firls(numtaps=21, bands=f, desired=r)\n",
"\n",
"# Plot the tap values\n",
"plt.plot(taps)\n",
"plt.xlabel('Tap')\n",
"plt.ylabel('Value')\n",
"plt.title('FIR taps')\n",
"plt.grid(True)\n",
"\n",
"# Create a second figure and plot the actual frequency response against the desired frequency response\n",
"plt.figure(2)\n",
"w, h = scipy.signal.freqz(taps, worN=8000)\n",
"plt.plot(f, r)\n",
"plt.plot((w/np.pi)*1.0, np.absolute(h), linewidth=2)\n",
"plt.xlabel('Frequency')\n",
"plt.ylabel('Gain')\n",
"plt.title('Frequency Response')\n",
"plt.ylim(-0.05, 1.05)\n",
"plt.grid(True)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Calculate the passband and stopband ripple"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Passband ripple: 0.902 %\n",
"Stopband ripple: 1.213 %\n"
]
}
],
"source": [
"onequart = len(h)//4\n",
"threequart = onequart * 3\n",
"\n",
"# passband ripple\n",
"pbr = max(np.absolute(h[:onequart])) - min(np.absolute(h[:onequart]))\n",
"\n",
"# stopband ripple\n",
"sbr = max(np.absolute(h[threequart:])) - min(np.absolute(h[threequart:]))\n",
"\n",
"\n",
"print(\"Passband ripple: %0.3f %%\" % (pbr * 100.))\n",
"print(\"Stopband ripple: %0.3f %%\" % (sbr * 100.))"
]
}
],
"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.5.2"
}
},
"nbformat": 4,
"nbformat_minor": 2
}
@philpem
Copy link
Author

philpem commented Oct 19, 2018

Refer to https://scipy-cookbook.readthedocs.io/items/FIRFilter.html for info on Scipy's FIR filter design tools

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment