Skip to content

Instantly share code, notes, and snippets.

@danibene

danibene/filtering_lowfreq_engzeemod2012.ipynb

Created July 22, 2022 05:55
Show Gist options
  • Save danibene/12e95a7cdf881d5fc1a8eab742ffa0f2 to your computer and use it in GitHub Desktop.
Save danibene/12e95a7cdf881d5fc1a8eab742ffa0f2 to your computer and use it in GitHub Desktop.
Download ZIP
filtering_lowfreq_engzeemod2012.ipynb
Display the source blob
Display the rendered blob
Raw
filtering_lowfreq_engzeemod2012.ipynb
{
"nbformat": 4,
"nbformat_minor": 0,
"metadata": {
"colab": {
"name": "filtering_lowfreq_engzeemod2012.ipynb",
"provenance": [],
"authorship_tag": "ABX9TyPXVJCxyz+Iw9/9+WHcC7ui",
"include_colab_link": true
},
"kernelspec": {
"name": "python3",
"display_name": "Python 3"
},
"language_info": {
"name": "python"
}
},
"cells": [
{
"cell_type": "markdown",
"metadata": {
"id": "view-in-github",
"colab_type": "text"
},
"source": [
"<a href=\"https://colab.research.google.com/gist/danibene/12e95a7cdf881d5fc1a8eab742ffa0f2/filtering_lowfreq_engzeemod2012.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "C23RXuWzy_w0",
"outputId": "2c13a593-82e0-47a3-a8d5-cca558b294d5"
},
"outputs": [
{
"output_type": "stream",
"name": "stdout",
"text": [
"Looking in indexes: https://pypi.org/simple, https://us-python.pkg.dev/colab-wheels/public/simple/\n",
"Requirement already satisfied: neurokit2 in /usr/local/lib/python3.7/dist-packages (0.2.0)\n",
"Requirement already satisfied: scipy in /usr/local/lib/python3.7/dist-packages (from neurokit2) (1.7.3)\n",
"Requirement already satisfied: numpy in /usr/local/lib/python3.7/dist-packages (from neurokit2) (1.21.6)\n",
"Requirement already satisfied: matplotlib in /usr/local/lib/python3.7/dist-packages (from neurokit2) (3.2.2)\n",
"Requirement already satisfied: scikit-learn>=1.0.0 in /usr/local/lib/python3.7/dist-packages (from neurokit2) (1.0.2)\n",
"Requirement already satisfied: pandas in /usr/local/lib/python3.7/dist-packages (from neurokit2) (1.3.5)\n",
"Requirement already satisfied: joblib>=0.11 in /usr/local/lib/python3.7/dist-packages (from scikit-learn>=1.0.0->neurokit2) (1.1.0)\n",
"Requirement already satisfied: threadpoolctl>=2.0.0 in /usr/local/lib/python3.7/dist-packages (from scikit-learn>=1.0.0->neurokit2) (3.1.0)\n",
"Requirement already satisfied: kiwisolver>=1.0.1 in /usr/local/lib/python3.7/dist-packages (from matplotlib->neurokit2) (1.4.4)\n",
"Requirement already satisfied: cycler>=0.10 in /usr/local/lib/python3.7/dist-packages (from matplotlib->neurokit2) (0.11.0)\n",
"Requirement already satisfied: python-dateutil>=2.1 in /usr/local/lib/python3.7/dist-packages (from matplotlib->neurokit2) (2.8.2)\n",
"Requirement already satisfied: pyparsing!=2.0.4,!=2.1.2,!=2.1.6,>=2.0.1 in /usr/local/lib/python3.7/dist-packages (from matplotlib->neurokit2) (3.0.9)\n",
"Requirement already satisfied: typing-extensions in /usr/local/lib/python3.7/dist-packages (from kiwisolver>=1.0.1->matplotlib->neurokit2) (4.1.1)\n",
"Requirement already satisfied: six>=1.5 in /usr/local/lib/python3.7/dist-packages (from python-dateutil>=2.1->matplotlib->neurokit2) (1.15.0)\n",
"Requirement already satisfied: pytz>=2017.3 in /usr/local/lib/python3.7/dist-packages (from pandas->neurokit2) (2022.1)\n"
]
}
],
"source": [
"!pip install neurokit2"
]
},
{
"cell_type": "code",
"source": [
"import neurokit2 as nk\n",
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"import scipy\n"
],
"metadata": {
"id": "w73hineB3I84"
},
"execution_count": 2,
"outputs": []
},
{
"cell_type": "code",
"source": [
"sampling_rate = 100\n",
"ecg_signal = nk.ecg.ecg_simulate(sampling_rate=sampling_rate)\n",
"clean_ecg_signal_current = nk.ecg_clean(ecg_signal, sampling_rate=sampling_rate, method=\"engzeemod2012\")"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/",
"height": 360
},
"id": "5GbUw5uUdcD5",
"outputId": "15a292f3-42f3-406f-ae64-c2b946927300"
},
"execution_count": 3,
"outputs": [
{
"output_type": "error",
"ename": "ValueError",
"evalue": "ignored",
"traceback": [
"\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
"\u001b[0;31mValueError\u001b[0m Traceback (most recent call last)",
"\u001b[0;32m<ipython-input-3-09d1e0ad3246>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m 1\u001b[0m \u001b[0msampling_rate\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;36m100\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 2\u001b[0m \u001b[0mecg_signal\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnk\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mecg\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mecg_simulate\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0msampling_rate\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0msampling_rate\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 3\u001b[0;31m \u001b[0mclean_ecg_signal_current\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnk\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mecg_clean\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mecg_signal\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0msampling_rate\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0msampling_rate\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmethod\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m\"engzeemod2012\"\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
"\u001b[0;32m/usr/local/lib/python3.7/dist-packages/neurokit2/ecg/ecg_clean.py\u001b[0m in \u001b[0;36mecg_clean\u001b[0;34m(ecg_signal, sampling_rate, method, **kwargs)\u001b[0m\n\u001b[1;32m 108\u001b[0m \u001b[0mclean\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0m_ecg_clean_elgendi\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mecg_signal\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0msampling_rate\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 109\u001b[0m \u001b[0;32melif\u001b[0m \u001b[0mmethod\u001b[0m \u001b[0;32min\u001b[0m \u001b[0;34m[\u001b[0m\u001b[0;34m\"engzee\"\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m\"engzee2012\"\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m\"engzeemod\"\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m\"engzeemod2012\"\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 110\u001b[0;31m \u001b[0mclean\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0m_ecg_clean_engzee\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mecg_signal\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0msampling_rate\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 111\u001b[0m elif method in [\n\u001b[1;32m 112\u001b[0m \u001b[0;34m\"christov\"\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;32m/usr/local/lib/python3.7/dist-packages/neurokit2/ecg/ecg_clean.py\u001b[0m in \u001b[0;36m_ecg_clean_engzee\u001b[0;34m(ecg_signal, sampling_rate)\u001b[0m\n\u001b[1;32m 260\u001b[0m \u001b[0mf2\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;36m52\u001b[0m \u001b[0;34m/\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0;36m0.5\u001b[0m \u001b[0;34m*\u001b[0m \u001b[0msampling_rate\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 261\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 262\u001b[0;31m \u001b[0msos\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mscipy\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msignal\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mbutter\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m4\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m[\u001b[0m\u001b[0mf1\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mf2\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mbtype\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m\"bandstop\"\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0moutput\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m\"sos\"\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 263\u001b[0m \u001b[0mzi_coeff\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mscipy\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msignal\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msosfilt_zi\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0msos\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 264\u001b[0m \u001b[0mzi\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mzi_coeff\u001b[0m \u001b[0;34m*\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mmean\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mecg_signal\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;32m/usr/local/lib/python3.7/dist-packages/scipy/signal/filter_design.py\u001b[0m in \u001b[0;36mbutter\u001b[0;34m(N, Wn, btype, analog, output, fs)\u001b[0m\n\u001b[1;32m 2955\u001b[0m \"\"\"\n\u001b[1;32m 2956\u001b[0m return iirfilter(N, Wn, btype=btype, analog=analog,\n\u001b[0;32m-> 2957\u001b[0;31m output=output, ftype='butter', fs=fs)\n\u001b[0m\u001b[1;32m 2958\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 2959\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;32m/usr/local/lib/python3.7/dist-packages/scipy/signal/filter_design.py\u001b[0m in \u001b[0;36miirfilter\u001b[0;34m(N, Wn, rp, rs, btype, analog, ftype, output, fs)\u001b[0m\n\u001b[1;32m 2422\u001b[0m raise ValueError(\"Digital filter critical frequencies \"\n\u001b[1;32m 2423\u001b[0m \"must be 0 < Wn < fs/2 (fs={} -> fs/2={})\".format(fs, fs/2))\n\u001b[0;32m-> 2424\u001b[0;31m raise ValueError(\"Digital filter critical frequencies \"\n\u001b[0m\u001b[1;32m 2425\u001b[0m \"must be 0 < Wn < 1\")\n\u001b[1;32m 2426\u001b[0m \u001b[0mfs\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;36m2.0\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;31mValueError\u001b[0m: Digital filter critical frequencies must be 0 < Wn < 1"
]
}
]
},
{
"cell_type": "code",
"source": [
"sampling_rate = 150\n",
"ecg_signal = nk.ecg.ecg_simulate(sampling_rate=sampling_rate)\n",
"clean_ecg_signal_current = nk.ecg_clean(ecg_signal, sampling_rate=sampling_rate, method=\"engzeemod2012\")"
],
"metadata": {
"id": "NRarXl3h28TR"
},
"execution_count": 4,
"outputs": []
},
{
"cell_type": "code",
"source": [
"plt.plot(ecg_signal)\n",
"plt.plot(clean_ecg_signal_current)"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/",
"height": 285
},
"id": "RdpxXqn93PAZ",
"outputId": "f072336f-1869-4f38-aed1-6bfaba2cbaed"
},
"execution_count": 5,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"[<matplotlib.lines.Line2D at 0x7fb0fcc27690>]"
]
},
"metadata": {},
"execution_count": 5
},
{
"output_type": "display_data",
"data": {
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
],
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAD6CAYAAACs/ECRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nOy9ebQkV3kn+PsiIvftvVf1alFpKaEFAQYDLoRt5K3dYNz2gHva7ga3bexjBveZ9niO7TltfOYc7MFtN9ge77iBNmrb092WbdqAphGDAS8sAqQCCQm0lqSSVPtbc4nMjIzlzh/33ogbmRmZWS8y7nt6L3/nSJUvl/giM27cb/t930eMMSywwAILLHBwYez2CSywwAILLLC7WCiCBRZYYIEDjoUiWGCBBRY44FgoggUWWGCBA46FIlhggQUWOOBYKIIFFlhggQOOuSgCIrqTiK4Q0dcTXv/XRPQQET1MRPcS0Tcrr50Vzz9IRKfncT4LLLDAAgvMDppHHQERfSeADoA/Z4x905jXvx3Ao4yxLSL6fgC/yhh7rXjtLIBTjLH1WeUdPnyYnTx5MvV5L7DAAgscJHzlK19ZZ4ytDj9vzePgjLHPEtHJCa/fq/z5JQDXppF38uRJnD69cB4WWGCBBa4GRPTsuOd3I0fw0wA+ofzNAPwtEX2FiN6xC+ezwAILLHCgMRePYFYQ0feAK4I7lKfvYIydJ6IjAD5FRI8xxj475rPvAPAOALj++uu1nO8CCyywwEGANo+AiF4B4E8AvJkxtiGfZ4ydF/9eAfARALeP+zxj7IOMsVOMsVOrqyMhrgUWWGCBBXYILYqAiK4H8DcAfpwx9oTyfIWIavIxgDcAGMs8WmCBBRZYIBvMJTRERH8B4LsBHCaicwB+BUAOABhj7wfwLgCHAPwxEQGAxxg7BeAogI+I5ywA/40x9v/N45wWWGCBBRaYDfNiDb11yutvB/D2Mc8/DeCbRz+xwAILLLCALiwqixdYYIEFDjgWiuCFjof+GnDau30WCyywwAsYC0XwAobz/APA37wdrb/6X3f7VBZYYIEJYL4Hb/2Z3T6NRCwUwQsYT56/AgC48NyZXT6TBRZYYBIe/n/+D1h/9EpsXnhqt09lLBaK4AUM2SbKQLC7J7LAruNzT67hcqu/26exQAKKz38BALB+6fldPpPxWCiCFzB8EADAoF0+kQUS8bEHz+Pur13IXM6Pf+g+/OAffj5zOQvsDIzETcr2ptG2UARZ4d8fBT7+i5mK8BhfXIT0HWQXyAb/+10P4uf+4quZy3kVPQm3PXMD3wU0gwmjjQX+Lp/JeCwUQVbw+sD9f5KpiEVoaO/jN60P4BP5d2Yu5yOFX8EnCr+cuZwFdgYfJn+wRz0CrU3nDgoYk/o/WwTCI1ho872Lf2n9ozZZx2lTm6wFrhIyNLTwCA4O1ppdLXJ84REsQkMLLLC3Ecitdo96BAtFkAE8t6dFjs/45VuEhhZYYGfotJu4fOG5zOUwWiiCAwc20KQIxL9Zh6G2z9yHcx9/b8ZSFtgp5jFu9qDi4u9+F45+8OWZywmDxYvQ0MEBeY4WOSzg1kXWoaGl//J6XHv/b2QqQzeutPr7hncfLPTAjnFLoKfaNxDJYrZHPYIDlSx+eq0DL2C49WgtUzkB06P1Zf6JFqGhq8Y73/NbMMDwJ7/xrt0+ldSQBsECexeyjmCv0kcPlCJ49/94BJv2AHf/7B3T35wGmkw0Q4QESFdogLFI+7zAcWf+t8WjF74iCBaKYM9jkSzeQ7AMgq9hkw40XWwmQkK6PAIWeFrkLHB10LXeAAAP/Ffg4Q/rk7dPsNeTxQfKIzBIjyLQlrzTFIKSCHwfppnTKnOBGaBzc/mY6HT78h/WJ1MHMvZ2mbS596j3dqA8AlOTR6Brg5YhIV11BIHvapGzL5GhcRDs0bjzCwoZr22pYvZqXv9AKQLDIPgarHWmKUcgGQikySL0vWxDQ1v2AL/9ycf1KGvdyHDd6Vpv+xnMH2QrQPacy1bKjnGgFIFJhEBLjkCTIhBupp6GFtl7BL/2sQfw+D/ehb9/7EqmcnYFGSrrfekRuD2gs6ZNnO9l7RHwPUHX3nC1OFCKwNLkEWiL2WpeVIGfrUfwxisfwn/K/w4qF7+YqZzdQJb8cbZn7cydw/ng64HfvlmbPM/N1iOgsPvo3rxWc1EERHQnEV0hoq8nvE5E9AdEdIaIHiKiVyuvvY2InhT/vW0e55MEwyD4vgaPQFNCSHdxSpBxaGjF5xZgrns5Uzm7gSy5/kzDmgb0VjAX1h7WJgsAgkBP/muv1hHMyyP4UwBvnPD69wO4Rfz3DgD/EQCIaAXArwB4LYDbAfwKES3P6ZxGYNJ+8whkaChb+KLLaZAxfdQzigAAw9fTokMnsrTamS768N40ZucC39WkCPao9zYXRcAY+yyAST1w3wzgzxnHlwAsEdFxAN8H4FOMsU3G2BaAT2GyQkkFwyD4Gu4ZbTeMJoUji2GyzhF4RgGAvhYdOpGlJairsnivxrfnAW1FeXv0N9SVIzgBQB3WeU48l/R8JjANTYt5ny2qQMY3/WzdWqkIDG//eQRZkhS0KYI9Gt+eB/zMf0Px2+3RgrIXTLKYiN5BRKeJ6PTa2s7YBJZh7KvK4mhRZcsakqykIGOKHZPFavuwXoFlWFuiLSe1j3taLTwCPTgP4Drl72vFc0nPj4Ax9kHG2CnG2KnV1dUdnYSuymJdWl9XsjhUBBnfLOEQp715r6RClolWbTmCPVoVOw/4OmLGwIFXBHcD+AnBHvpWAE3G2EUAnwTwBiJaFkniN4jnMoFpQJNHsL9YHJIDnTX1zZQKZ4/eLGmQ6W+nyyPYjxpaQFcSd1+3oSaivwDw3QAOE9E5cCZQDgAYY+8HcA+AfwbgDIAugJ8Sr20S0a8BuF8c6t2MscwGr+qqLNam9TUvqqzjqJFHsD86nKrIcgPQFdbYl4VrAtq8nT1q5MxFETDG3jrldQbg3ya8dieAO+dxHtNw29Y/4M+MPwdns2aHLOPBQ5I0yeHIeiOgsAxf4/fS1Fo7Sy9Hl2e4V6mPc0HGRlXoVe9nj+CFgn/+xDt5MCwIACPDqJgurR+2mNCDzEMDu2Et6ZqxkCl9VI/hsb89gmzXXri096hH8IJhDc0VGVcR6nMzNVcW71FrJhU0z47I5NjaelvtkqLWIkbX2t6b99CBVARZdxrUl3jSlSyWAvfmIk4FXYogQ1aKLsNjV8Ia+4x4sVcT7gdSEWTdM0dbPxF5Y2oaH6mv14zOZPEL3+LUZnjsikew35Tc3jSmDqQiyLz3uO4WE9pGFu/NRZwG2uLrmc4j2G+bpSp0f9RIRMnihUewZxBk3Ht8r17stMjaItwN0qi2a5Vl99F9ZzWrQnXJ1BRm3aN7w4FUBCzrlrO66KPSmtW1g2quZNYiax9sovp6De1fRbCvvaoZcCAVQda9bPRVxrLYP1nL2auLOA10USIzZVxp722lEdoUtbb4qh45V4kDpQjOL50CALCMQ0P7rfvoLonTAm0J0AzXhC7DYze6j+oqztRm5OxRY+pAKYJvnPhXAADGsmUN6eY+M22soaxvyt3YaPZBHYFulppGaAvZaGpDvVfzhwdKETCDF1Jn7RHo6gZJum/MPbqI02A/bDS6NpfdaAaob+Pcb3KuDgdKEUi+fea8631WBBPJ06R4NHk4gMZCokx7De3fRKeuHM5+7QYwKw6UIpDbS9Y3v77EkybPI0xK701rJg32w0azH/IcSdh3XUH36D10oBQBSHzdrG8cbTFb3R7B3lzE6fDCryzed8yamMx9wOqCakwtPIJdR9TmOFtoI1doGlUZicu4DXWYUMtUTAz6Eq0v/NDQbnQf3W8ewV41pg6WIhD/Zp30Im3zCPS0oY5CahkL2gVkOlRe+cEyLSjTFtbQIyYmct8Veu3Nm+hAKYIwNJTxRd+3xSm6FJxO9zlTS115nOmoSl3hk93wCHTVEeiaR7AIDe06KGQNZQt9jDfZdG6/KB79Sekswx2BJo9A24zsXSko2x/Ei6iV+8Ij2H3ImbhZu5uaLCeSoaF914Za86jKjKDum1kqAm31JLtBH91vifCFR7D7IE3ZYn2DaUKBmUI740HjzZKlR6Cug33hEexGZfEB36B14WApAuESZL649mlxSta/m2zRq5NZkW1XUOV7ZBlW0cUa2oWwBu2Dyu+4oEVoaPehq7JYW2/z/aUIFEma5CDbZLGyuWTZdkTbJrYLBWWL0JAezEURENEbiehxIjpDRO8c8/rvEtGD4r8niGhbec1XXrt7Hucz4TwBZK+UtVHetFPRMmZbhQ/03SzZhmwURZApTVWXR7AboSFdkvQ0ndur9FEr7QGIyATwPgCvB3AOwP1EdDdj7BH5HsbYzyvv/98AvEo5RI8x9sq05zHTuYbnk3WyWK9HoM1O11R9qXOKU5ZKW80/ZFohuw+mrCVhv9BH93qblnl4BLcDOMMYe5oxNgBwF4A3T3j/WwH8xRzkXj1IT45gv/U2Dxdx1vRBTTmCXSn0yjRFoCt8okXMkMz9oeQkw28/h4ZOAHhe+fuceG4ERHQDgBsB/J3ydJGIThPRl4joh5KEENE7xPtOr62t7ehEKSwo29HHrwJ6Kz31eQSZx9Ti/2YlRhOtk/lKaCjT76TL8NiFgjJN3y3rMGtE8N6/HsHV4C0APszifvINjLFTAH4UwO8R0U3jPsgY+yBj7BRj7NTq6uqOhOvqPqpP62u+MbMODYXH1+gRZBm7hx5FoG1CWVyDapGpq4gte0bc/vcIzgO4Tvn7WvHcOLwFQ2Ehxth58e/TAP4B8fzBXEGGnspifZW+esRoU6CapjjFkp5ZbtBquGEftJiI/VaZtuaI9ebITE4MC/poatwP4BYiupGI8uCb/Qj7h4huA7AM4IvKc8tEVBCPDwN4HYBHhj87b2Qdw9cV14ySxboGuWSdZOfHp4zlxJK4miaHZRriUJdblpXSMWWmp0eTtjnJmSeL9Xi7O0Vq1hBjzCOinwXwSQAmgDsZY98goncDOM0Yk0rhLQDuYvFd8iUAPkBEAbhSeo/KNpo3aL8mi3WtrX1iNanhhmwLynTRRxWPgAUAmZnIiSloFoDf7vNHXN1oCg1lbHwYezw0lFoRAABj7B4A9ww9966hv391zOfuBfDyeZzDTBDJ4uxzBPur8MogGbLJWJCmJnqxDVpTQVmWiVamKXbPNIWGAsZCFcP8fdbxdh+Hhl4w0NYBcJ/RRyNoShZnrQg0JXFVazbThK6mmPpuhIZ0GTtZbwl7PTR0sBTBPhter63rXChOV/Vl1jkCNWSjp6Asy3Yg8euiqaeRpiS7vjqCjD0P+T0WHsHuQ1uLCc2FXqTLytgnIbXdCNlkuqExPR5OoCk0pCqzfUMfPQCVxS8YhF2o98mGFrFs9kv1pbxZNHoEmip+M63+1SQn3gk0S3aSwurSVSyXeUGZHkbcTnGgFIH8urr48JlD94SyzENqmnIEcX6iHjmaePeZJr810VRjrK590tJ9N1qsXw0OlCLQ5xHotWIMXWX42qovM74+u9AeOtu8lK7kt6YcQSzUtb9CQ7pbx8+KA6YIwhrZbAVpKyjTbF1oShbvm4KyWAhKz2zkLMewMk2hIVWZ6ct/Hew21AdKEQD7K1msqxI3FKcroZZ1ryFtIyQV1lCmPec0WdC6urb6ehRbTOYiWXxwYMjmo/skR6CbNZS5W6spjhqbU6wtdp9liwk9MXX1O2T7dfRcnyGhmR6edBM7rhIHShFQmCzW0zNHCMtcDum6VzQli7NWOEzTxhmjdWbqtakDcPYDTVUtxNsnyWJNDRV3igOlCBjp6T6qi12hK6YeStMWGsoYKjd9Hwyvj9sder5Ppht0rGvr/lrbi2TxHoAR0YYylUPDTcCygmbrQldoKOvvpYuVwjRZ0Los9VhuJUPFFsRCavuDGq0r/7VTHChFIJE9fTTxj7lCbsz7hj4KXaEhPWye+KasK3afpcKJfqsgw5YM2lpmqFKyLpbU3A7manGgFIGxGy0mtJTi7xOKnaabJT5CUk9hVJbrgDQlpZmupnO70WtI1zyCRY5gDyBUBLo4w9nKkovL2Ce9hnS1mIgr6he+pc40hbqgKWSjLZkfg55k8SJHsAcQFZRlDE0bgH53U1NlccbQRevU5RnqmsEcq2DOtHBN7TWk53fTZuQsQkO7D4KeCWUxyynTtsCixYS2Vr2aPCmdHkGms4Q11RHE6KN6aJ3aktKZfh9V6KKg7OBADK/P9OYHYhc728pIeez95RFkTSNlgR5WV6Ar6alrQ1PWcpaDdgJfVdQZXp9YLmKRLD4wMDTVEcSbgGWfJNTFv98vVD59cwJ0Fa7pZ/Nkuxb0J6Wzp49KRtxCEew6ojICXeyXbOmJUgHoSxZnnVDT1IZakyWorVVCLNSdoeEBPYpNF2so5rFlPmsjFJSpnJ3iQCmCsOlc1mJilmCWcjSHhrT1Y9HXYiJb1pAey1b1CDJtq61p49TleSw8gghzUQRE9EYiepyIzhDRO8e8/pNEtEZED4r/3q689jYielL897Z5nM+E8+QPNM4jyNZCy54+Gj9/XX1fMj68Jktdl2WrHpqyzH9pIkHomimts7ld5BHsTUVgpT0AEZkA3gfg9QDOAbifiO5mjD0y9Na/ZIz97NBnVwD8CoBT4L/QV8Rnt9KeV8LJAsg20QUMudAZhoYiZLvJkPpHhojoo/vDI4jHbPS0mMhybau/W6a9hjTROgNtrK6DUVB2O4AzjLGnGWMDAHcBePOMn/0+AJ9ijG2Kzf9TAN44h3MaC0NTHUG8kEhHiwlN1ln2MTUA2bvP2nIEvq6uoLtQF5FlEz1dXU41tSMHont0P88sPgHgeeXvc+K5YfwLInqIiD5MRNdd5WdBRO8gotNEdHptbW1HJxrOJ9NY3ZetrOxZQ6qFScjWu9E2mGYX+PDjPIJut4PTd78frtNLLSkSk2VX0N0oxNOTw8m6xfpenUMgoStZ/P8COMkYewW41f9nV3sAxtgHGWOnGGOnVldXd3YWu5EjGLMALjz7BB7+7EdTi5lmOd9316/joY+/P5WM+E2Z6lCzSJNCsxUT6LFsMaXX0P13vQenvvpL+Prdv5dKjC4PRxctOt6jKTMxsd8t66r2aS0mnnr0AaxdOpfpOUzCPBTBeQDXKX9fK54LwRjbYIw54s8/AfAts352njDEiDK9rKHRC7/9p2/By//ubVg/fyaVmChZPCqj2e7g9sd+E6+4/5dSyYidf8Y3i7EbBWWZsmwmyylvfh0A4K0/lVZS9GicR8AYAneQUga0kSAY9Cg2nS0mjPB7jMrp9rq46S+/G/n3vzbTc5iEeSiC+wHcQkQ3ElEewFsA3K2+gYiOK3++CcCj4vEnAbyBiJaJaBnAG8RzmSLzdspTeg29lPEb/+KjX0onKCwoG8XzzzwWPnbs7R2LiA163yedIOMhGz2FUeNyRdX+JQBAoZPS9pmyWX75L9+Dzq/fiO0r6SzOeI49w8I1ZWZxpr2GNFWYx4WOfp+zT3CDoIEOAs/Vcx5DSK0IGGMegJ8F38AfBfBXjLFvENG7iehN4m0/R0TfIKKvAfg5AD8pPrsJ4NfAlcn9AN4tnssERHJUZVYSJJItp62tiBDV33geaUATElC9S5G3sfb8kzuWoY0LD40tJnR1BZ3CU697G/zf/oWUgibH7l/y6B+gjg7Onv5ESjmamsFp8ghi7V8yTxYn18i0Lkb35/Z6yrWwQ6SmjwIAY+weAPcMPfcu5fEvA/jlhM/eCeDOeZzHNISkIZ2VxUMvXTn3JJbla610F33Sxultng0fNy+ewbW3vWZHMgKtVpPMEWRMH41Z55qSkUMbzcD1cZhtAgTU/J17bPzgkxWbD5PLvLxzgwAY9giyzK2osXtN7KTMk8WhoBEEyr3aWjuHlWM3ZHou43CgKoujNtSZVyxFD4dcaKd5OXycsy+mFQQAMMd8H9eOPI9Ba2csK2C4aZ4eZkXmJF/Vss2UBhnJGU4SNrfXUCAPA2ahxjopLdJkr83uOaijAwAw7MtIA9L2u+0ddtKzD/4D3H4ntSya4BGwdnRd7PXMUqQTcSAVgS4+PDB6www6Tf4vs1Dor6eSMql+IOhGVqbf2di5EOX8dTErslY48UuSZWGU+jgup9viEdBLxlHkyYPn2CnkKAVlQx5Bc+MCTOInku9d2bkMDG/Q2eUI4u2uMxMzdYLck48/ghs++maced+PpJYV3qtj5AS9Zvh4sL1QBJkjqiPQmCwe2mi8LrfULxlHkfdSWhoTBp+YgxY2aBkOywHdnaddggky5g1dvYagq5AoFlOPo9/m62C7cAwA0Gnu/BpNYr84m1H4sTxIZ3hMWm9zxS5MkBtn5Fx6+DMAgFvaKUkdUL3d0d/NGLSwTjxg7HYyS5FOxMFSBDJZnHloKLk5V9Djlvp24ThKQTpFEFtUQwvZctvomTVsowqzn2KTiYW29glrSD3+mJj6ubOP42vvfQPOP/lgZnKcDlcEvTKvn7S3d75JxwyPoXUgQ5HP4SiqXrpchK4CLH2VxZOPbaw/DgCwUU4ta5K3m/c6aFuH4DAL6DdHXteBA6UItBWUxXoNxWWxfgsA0C0dRzlIEQ4YkjP8nQpeB45ZQdtowHLS0EcnV8fOE/oqi5M9NgB46jN/im/ufRmXP5W20Ct54xzY/IYP6ryMptdOEb6bMNFrIKjDm7ljKLOU621CCGqeYLrYPFOaAlpdHkproAPfdUZevxpEBWWjcgp+B65VRYcqMJxWKjk7xYFSBFNDQ1vPAq20Cdz4xR6WRf0mbBThFVZQZXaqhW5McKGLfhsDq4aeWYXltXcsI96hcceHmQlhaEjjKNFxv7/ZfBYAYHVSroUJ3S09kcOxDp0EADhpFMGEAkYpxy5egwrrpmsfrakAS1+bcPXYo3IKTuSlNdd2XoPBGFPyeaO/WymwMbCqsKkCc7BQBJlDVhaPxdf/O/D7rwB+57Y5SEpWBOagBRsVoNRAjvyUjIRktkiZ2fDydbhWBYUUuQh9nTojRZ19HcFkSmyjx2/6SsqY+qTfTiYIq0dfBABwOyka7k6o+PWFHK92AiYx9LtpQg+aqsw1VTAHU5LFVXcDHuN7RvPKzmt+AqY2nRuVUw668PM19IwKLHehCDJHGBkaQ31rPfix+cmJ8brj7ArLbaNrVEDFBgCg09y5JZiUIwgChiqzEeRr8HNVFPydhwSmFUUxxnD6M3+N7bU5eFK6ZjBP2QCWPBES8NKEazAxuRr0uZe2dM1NAAC/u3NFMJFu2W/BZwSzwYv77ea85GRZ8auLoDBZTiPYxvMmD911mzs3ChhjMGh8aMgPGGqwEeTr6Js1FFJ472lwwBSBnFA2etEv4PAcJSVbNHmvg75ZhVXhLIFeaz4hAfVxx3FRQxcoNuDnayixbgoRk2+WR77yWZz63Ntx/j+/bccyFAEAsg8NxZlQcUXNGAtzNytsG8zfecn/pKE+5DThsBwah6/hr/bSJHInWNBOCx2Ukavy9dZNEYKiSQpnjmATPOp5Qs1zDDPVGGOoMRvNAr8+nr1zBRpMKGDs9AeoogcUahjkqij66WsWdoIDpQjCUZVjFteaulemXXwT2jIU/A4GVg1WWSqCnTN64s3mIjmtThsF8kDFJbB8DdUUiiCeFBzT3O7xzwEAVrvpGugBkyl2c0Us9BB/qdN3UYONDsowiKG9maIIa4ISNQZt2FRGqZCHzQqAs3NLMKY4h3IAppCTF+vNae9c4Uybs9Hv9/HYfX+bum0Hm7Lm5gaW3EfL7tookot+lbO6vG4KwoVKUx1a253WNgxiMIoNeLn6HAgkO8OBUgSTCosHTrRZ9rvp4nTqjTnMrigHHbi5GgrVJQCAk4I3nBQaku6/WV4CFerIkwenv0NlMM1Nb/J4eh7pWBWA+rvpTBbHr0/b7iBPPjasIwAAO4WinhRKMWWIkAg2lUGDFCGBCbRO022jZ1SQFx6BkyIXQVMs9a/8+S/htnt+BI9+/iM7lsEPPnnNMcbw9S98HH6aIjzEldlwHUFnm1931rgeQEqPbcL36Yp6EqPcQJBvoIqFR5A5ZB3BOCvj2zajhqkbl9M2fkqOQVdE7L5UWwEAuHaa5N14Ob02X8S5cgMo1gEAndbONoBpRTeWaJOxhA68VInIyRS7eWJSa+2u+J1aIiTQTaEIMOG3y7sd9I0Kl0FlWKkUQXKvobzbQd+soCgUQRrLNu5tjFYWV9e/BgBoP/XlnctA/DuMUzj3f+oufNOnfhRf/W+/mkpOMMFz74p7qLB0nBdl9nduHKr9uoY9AqmYrXIDrFhDES6CQdpBRVePg6UIjIQWE4Ef41h3t9P1ZEliPQR+gCrrIig0UBQegdebcfNkDP1e3KqnBMvWEYs4X1uBWeKKoLdDRRBM6SlTcaK2Bc2NdAnjmeoI5qAk2AS2lVSiAxESSOOxTaJb5v0OBhZXBH2jCstLY90mezgFv4OBWUOpzg0PP4UiiMXux7XV9oQn2klnSE1jqrnPfBEAUFx/OJWcSaGhvsil5CtLaFMZ5OzcyFHvoSRFkK8sg4TRZrezGdk+CQdKERhJBWU+H9rxGE4CAJxWup4slNBG17ZbyJEPFBuoihtTskem4Ysf+wCK7z2Oy89GcwYoIUcgi4iK1RXuFQDodXa2AUyjWtb9LdisCABopYmnY3JjLgD40h++DY//hztSK4NJzcbC2Q1LnC3ippjlMGlSXcm34Vo1LtOsIJ+GLTJheH0xsOHmqqhIRZCmcnVCboUxhsM+b25Y7F7auQwA09g8Vpv34ykO0rVjiG3Qw91hxf2Sry7Dpmoqfr/KHByRY0dyjGK6ezUNDpQiIJEsHtlmhCI4U3w5AMBr77xbJ4AhAy2SJsMzRqmOSqUCl5nAjJWESw9/CABw/sFPh89RQmjIE4urVFtGviJzETtVBJPpo5Wgg8sWD6PY2+l+t0l1BL4f4Fs3PooXD76O9fNPpJIziT4qN/6CKPRKwxZhE3j3JWbDy3FF4FpVFILZcjjNrTU8+euvxaOf/XD4nPp7DYegyiIUWa1U4DALbEbDY1NQiYEAACAASURBVBzix47L2dzaRIO4V1MbpG1uN7k9dFXMb1j2UsqZkJQeiOterK2gZ1SRc3f+u00KQfkinFquLcMsCUWQIqG/UxwoRZDYFFxQBL0qbwAWdNJuaOM3T0kVtcrLsCwTNkozs0UC0VPe33w2kpNA55OJrUpjBYUqX1yDnYYEYnHuoU3T81FDF60SD6O47bQNs1jsHxXnnlWGdzz3SEoxyRu0VAS1ozcCiHeGvHoxyQqHhwi5IvByVZRmZIuc/eqncYv7GOjzv6MKUsSo4RsZiqzDNAg2SjBS5CKGj62iK3olucxExU9ZFDWlzmNFKICVYBvMSzOCM1mO7N5brq+gb1Wvit//jc99DA9+4s7o0BP6dfk9/luV6svIVfi92l94BNmCooqy2POyj0ihxhs/sRTdOvnxx7MrZMdJS1jpNpVhurOxBCqBWIi96NxiISj1xhTuv1VeRknUK7g7TORO4kA3W03kyEe/dgN/tZuuACuaWTyGpnr+8fBxZyNtDDpZEQTixlw5dh0GzASbF1tE2QD6zgBV6gEFHhP28zWUZ6X4NnlYxPSU9yd4oE6vA4uCUE6XyjBTWLY0oV5B9kq6ZBzl8xVSIHZ9xlBRa8yGxwxO703RrG/SBh1WfjdWMLBqKF0Fv/9ln/kJvPLLPw/me0JOdOyR1vEi91CqLqMoFIGbJqG/QxxIRTBsZNhdnqUvl8qwUQKl4HQDyTS70N0UDI6+UYY5mG2ByY3CVJJW6gAXVQ45LXgwgVwJ5TqXFeyQ9TCpoKwjGDXGMqfYIUV1rIpxoSE1Vt/fShmDntDCgAklWq0vo43Zm4AFfoAv/95b8dDH3z/+2IrMTkt8F5EcZPkaquiFG8ckUIsrgkC9dROuUXh9ZMjBKMN0U1SZT6DD9kU4Yzt/DAVy4Ts7r11JUqAAMBi4qFIPlwxO7+2mKMicVMHMnCY8ZqBYrsPP1Wb22Hq9fvi4LUKlk+4hclpwYYJyZRSq6Yy2NDhQisDA+GSxLdg4hUIR3bScbgy3BY6sDrmZFQV1tG9UkJuhDxBjDBXBasoNVEUwPiltDlroUAUgQkUoAtabbUN74oHP4+xjUfvlSS0mbNFDv7J0BB1WnDnfkYRoruuoIlBpj0Z3fgNWhmPq5LQQgECFOmxj9iZgly+exWu378Er7v+l8cdWZEpKqowJS4XgdKevOznVLm51j9+gZT2JJQgDjlFJNQNjUo5ATsTrirbanVSW+iTjg8vZyvEwbipFMGEdGE4HNpUBIvj5GiqYTRFsrUWDZdrr/FpNqmA2Bm3e5poIpRqPFPgp65h2ggOlCCL66HAVYaQIrsZKFwcb00N8fOxRUvfKDa4IBlZlpj5AzmCACvHwVdGLFknc84gWWM5ro2tU+eNCCQNmzbRJM8Zw80d/ECfv+q7ouQmMh76kqVaXYVMJdDW/2wSMTRaLEE2TVZDvpwtBTUqAG4M2uigBhoGeUUFuxiZg68+NJrCT5PQkd7zEb3yzyHMFsxSvmaKleJ2p60BNrkbXq9+R9SRczsCqIJ+i71RsUx4K2chwhmyrbafpzTOhcM1uyTkOnKDQT0G1nLwOWugRn0PACnUU4YJ504smVcXUEZ5rrBZn6P3WgBf8AYiMtl1oRX2gFIH8usPbzKDNF22xWIRjVq6K033uM+8H3nM9upeUweAJpfhR2EEUk1kVFGZwOdtKY7qSH1mNSRXMea+DvlAEAETl6vRNumn3w+ZYA2GVxO/DYYod32hKtRX0qAxjxnxHEsKWGeOqSYVHc9k8FvOKdoQJg0+sQQtdsQE4ZhW5GddC+9JT0TFliCepzkMkA2VyMAzdzJAktMRvXEUPzO0LMeoGHT2UcgoiJ+VaVRRnZCeNR/IGLQkK1uGTANJZ6pjA5pF1Hp5QOKkmek0I2Viejb7B14Hk9/dmUDpqolcOBZqkcHJeJ1IEpQoGzExVvLZTzEURENEbiehxIjpDRO8c8/ovENEjRPQQEX2GiG5QXvOJ6EHx393Dn50nwsLi4eKe87wSko5/M6fyXYXVtH7/XwMAnvjK30dylNdji6Df5LOKi/zC+1YVpRluTOni91g+ShojOTRU9NoYWJEi6M2YlF47/0z4ePOKcHEnuOky1FVpHELfqISbVFqMHefnNNFBCV2zno5zj+GQwPAGwJsCAsDArKA4Y2yYKXMs5CaYFBuW3HEZE5YWe38G1pXKXpGJUkpIfof1JDVRVZyrzBzrHotJuRWRXC2tcrbVIA2DbEIuQtZ5WCs8L+Wl6tqavEHnPRt9k9+nhijK7M5QlDlQWnh4sh5JDQ2NK/gT96ppGuigDMzJs74apFYERGQCeB+A7wfwUgBvJaKXDr3tAQCnGGOvAPBhAL+pvNZjjL1S/PemtOcz8Vxl07mh563ms7jIVlA+fD08qzrzzQ8AJuPWn7OlMFkSNgBy2miL2D0ABPk6ypiuCKSLf8U4EhtmE3MzlS8l+xmFnzfLsGaIDffWnwsfd7f5Ig6YGhoaYluJcEC1cUgURc3HI0ia62pTGW6uikLaDo0TNrSCogj83OxrgSlWXEtWWCcUvslWHCWxQeeFZzCYod1Iwe/wlgeI2iCoUBucyVkEJUkYyNVQRor2BRMqfslpoYcCSstHAUTe4s7EJFOW5UZbPsLbdwdpKqUn0FQLgQ1XbNAyhNebQbmpjB+Z15pUWVz0IzmAYHalzFHuBPPwCG4HcIYx9jRjbADgLgBvVt/AGPt7xkJ+3JcAXDsHuVcNI4k+6jlwWA61ogU/X72qts0hrawdWYRGQmWk5bbQpUr0UqGKEgZTudDSxW/mjyJHftiLJE4fjTbsMuvCz9fDvwdmBfkZQhxOJ3Ln+60xjIcEil2uvCTCXGnCDtHxxykCy+2gR1V4uRrKKWc9T0oSFoNOeGN6udlpnaQklbtbosI6wVKX+Y6KIA3IBoSDGdgi5cDGeXYIQNQ2hBKuUUiBFKFIVqihADcMKU3C+WefxrP/18vw+Of++9jvMKxAjUELNsqoNPi5pZmvEL8/h3MR/DstHzkhegClqPOIxe6H14ENT6yD0GObIXTnK9eQxLnFjKlxA6RUo82YzWibN+ahCE4AUMf3nBPPJeGnAXxC+btIRKeJ6EtE9ENzOJ9EhMni4Q3NdeHBRK1oIchXUWGzW01lEbMv9lQmi3Izqklct42eGWl/GXvsdiYvZunid0WCTFqC42iqrh/wQReFSBG4MyalVWaOLA6L3/DD1LcmHOSAXPGqPalxCCuLx1jSBa+NvlVBUGjwkYtpMGGjKQddeHlxYxZqfC3M0FZZteKi+cPjLU7pPVQEaaAsQkT+tOI130MZPWxY3OoOO4kmFHqxfhMBI5RFUSGJAjZnBlbKua99Gjewc/A+/wfKd0gOE5puB12jilqdK4JUhXix7zM+F1FtrFx1D6B+z4bnRkbXJFonN6b4vZqXinqGdiMyD9hkFRgil8USPALOBuzCV+7VvlmZOS81T2hNFhPRjwE4BeC3lKdvYIydAvCjAH6PiG5K+Ow7hMI4vba2s8pfkl93aJ8JfBc+WbBMA1Soo0DuTG2b+SQwrr3VvifxZnDR46LXhqMoAkMqgilJKBmC8Wtcv9qhJTiqCNp2jzOMRN8SAGKTnv59VCvOs8VmFiTP3TUHLdjCwwny1dmLohIQeVKjikDOYEahjhIcBG6KitIJzeAqzIafE/z+Qh0GMQx60111y+uggxKAKHeSSL11WnCZCTPP3x/WekzZPCWbpF04JuTI6zV+QzOcFjpUBhm8Kp1mXG8AMNjgFewmxhddDRtTebeNvllBuVRElxXSUYkneKGhEq2vwKYqrBkVQeAH2HrPy/Ho7/5PysHGGzmMMVRZF0xY6kVB65ypc2u/iQEzsWmsRB1lE4b4OK6PKrpAIfIIBuYcQp87wDwUwXkA1yl/Xyuei4GI/imA/xPAmxhjIQ+LMXZe/Ps0gH8A8KpxQhhjH2SMnWKMnVpdXd3ZmSa0mAi8AQKDx10NQeULi34mwHbEdCHwbpLK2SrnrSRxh2L3YWfQaYpAbBDWMv+ZpSUYty74DWuLRCWJuCYwe+WqWkUb2KMewbBbm3fb6Al2EssL6zlFQzhBWBqtvoRo0pargwTDpjOn9tDq4/7AQxVdMHFjyrVgz9D7Je91sGHydSnDA0nzIgynHdZ5AEC11oDPaGofIBkidKo8shp6EAmxe1VRA4Al1lt/hu9jNHm+yAiiCW2TRrAW/A4cqybmK5RS1eKwCYoa/RYGzIJVKPMZvzNazxcvnsNx2sDLu19SvsP4XES310ORXEBc/7II4U312CBzWRV0jVpIPY63oVYKC9tNmMRABcVoy1VR9NOGWK8e81AE9wO4hYhuJKI8gLcAiLF/iOhVAD4ArgSuKM8vE1FBPD4M4HUAUjaSSUbUfHTIyvBdMMMCEG3O3RluFru1GdIt1RFzlBBLLQfxeGBeFvpMSRJKK6hwiJOt+kIRGAgQMBHuEu6nVCpWKVpcEJWr/pSW0tTnzBxeHNYWx01ux1Dw22FiFUVuPacZ6jNpZrGc4yCLsNIMjEkKcbQ7fCiNtJxDWucMScKib6Od5yGbsIo7wbK13HZIUQWAnGWig3IszzAOcr61ITqjylxDUiV7zm3DNqL1JtdEf4akdMnhjKSKn/DeoUtUDjrh2u5SJd18BSR3vDUG3MsBOL13VgbZ+nrUGdfrc+WRFBqSip+EIoiKMqf/btaAX9u+pVj2ym+lKhxptBmluCJIM1p2p0itCBhjHoCfBfBJAI8C+CvG2DeI6N1EJFlAvwWgCuCvh2iiLwFwmoi+BuDvAbyHMZaZIjCM8aEhBC4gPAIr5HRPd59lE7kuK8RoeXHLOVpsVWYjULS/ZIs4U2KPkpGRq/G5ymq7BV9eQrFhD/czAgAUayiQi05nsvUkOfQ2lWGITWlS7/6i38FA3Pym7KWeYjh6UvdR3w9QBZ/jIEd8zkK1TELSLGGpROWNGbbwnsEoKDMbXnEZPZYH9WUNxnjLNucqClSgS9MbwsmK5NKh67gHIavFE5oPFoZCkdLwmKUBodxga4GinCbUX5SZDV/kVrilvvPwRvzQQ/ReNyrAupqureqAodaWsEUT1oE0AmVb6EqpNHO4K+d10DOrcJX+REnJYrmuZD0JID1r/YrAmsdBGGP3ALhn6Ll3KY//acLn7gXw8nmcwywIJ1WOzCPwACsPQNmcpyRwgYi1ccVcxbX+Bb6wiBBbvMIKdwd9lGgApiSGZKHPtMlepmBkSLqhHGZjIIAPAzn44YYtaXuF6kr4eZmL6LS30KjHNyAVlisqkgM/rDuYVFlcCmy0c5zPLbnWs/ZS/+qn74Jz7iF820/+RnT8BPpop91EgwJQsY5cZXYGBwC0mpu48MdvBvv2n8NLvutHIL5U9IYxk91MEVaLPLbpsiqsh8u5Go/JCyWqeobq46LfgmPVYp/vGhVYU/jjUvnVlg6hjShRSgk5goJvY7sYEfRkJ1rXnr6hSQ+3DhvMG4CsfFwRqISIgKHGuggEU80xqyilYb5MYCfxAizJ6qqi1JstNDRQ2ol3tq5g5fiNiR6BNKZyZf59TIO4xzZDD7KCqOHx8nWUbXkPqeykUTl5xWhjhRry5MEf9MIckg4cqMpiOapy2CGgwAVMrgjk5jwLQ0BaVq3cMVgUIBAzVGNJXCGttT3qBoaFPlP6AFkDPt9WJhXVGLT0CGRlsSdi+8X64ejzM4a78n4HjlkVrRWkWzspocbDNYBiPc84E/fVn/8ZfNvZ96G1FtVfJHkEtgyJlJZCqqU7Iz3xsfs+jduch1D6/HuUrzHeUu+H06L4d5HKdFoTMMf1RNKvzov3hGWf1KSt5LcxyNXjx5jBipYhxEp9BR2qREylhBxBOejAU2jEcipeMAPlUh2ibkvvK4F2a3c7KJALEha0a1VQTNHKYtJ41ILXgSNonX6+Fvbgmga13qDfHKVGq2vO6Ua0aImuMVtRZjGwMbDqCAp1VCFqfoQcjxlQ76FogNRydB6F2RP688SBUgRGmCweijsGLsjkoSHZ+MmdoUmbbMPQL/PYcJTAZErsnsvqCEWghmyqNckWmSwrJ4qcKrUlBIzCnAGBRV0oxWJjgu1TakSKQG7S06zoom9jYFUxMCvIiRuZJUxxcrwAdRGuAWbPdwDxTXHt+WjiWlJBmcwHWOWlaMTnjLmI/gZnNqutOZJmCctiJTnovVhrCFmTv5Ntd8TkuTpXonJDVzYA9TtVWQd+vhE7hjNDRbssVqrUl0UcXv4G4xUOD0VGiiBsajbD2i4zGy3GY/FhvUJCG2oZDpSJfC9XQ2nWDdoP8IUP/TucOf234XOUYHwAQCnohJPdWL6GCvozdW1VGVmOpPcmsLpk6DWvhGy4cTTdIygFNvxcFazQgIUAQb8dFpQFMGJECFkVXaxH3rvMUc5CUJgnDpQiCD0C1YAa2LiRPY9BkW+c5au4WWSXwKA6zO8Poti92Axk24G8ErKpVGvwGU0dTsPL0GuoFPLooAhyIkUg5cgbMxCzFGrLR8LPz1q5WhLJ7EGuFoYGYqMqlY2g1WnHrMDQUp9BEdhdxdrciAhmSXUEsrI6X10KvaJpXpSEJebnespST6qNkJZ/SSiCSjhOdPJ36rZkcrEeT2DKawIKH/t+gBqzERSXYsfwrOntLORmVls+hJ5ZDSu5iTFhbUYynYGDKvXC6wMA1Wqd5xamrLeBywcOrZt8DYXJ8pgjpSZXI0UNXN18hSeffASve/4DOPw/fko5dvIsYb5GI4ICMFtdBBRF4LZlC5DxnqFcW6VqdI36ZnWIGTgK6SX7hUaoFLvtzdDz4PeqElIT3nulEbEgw4T+QhFkCKEIVCvQufAI8uTjwup3AgAq0kqfofGTL95jLvM4bE9uCCyy1OX6kmGHoqIILMEWmaYIyqLa1TAItsIuMRAocoQ13duCjSKsfDH8fFGGu6ZYtuXAhp+rxXogJcU3u6IFtaSplofyF5MgvSMg3jQsKUcQ9eZZCb2oWRgcAFDscbaIGuqQytlnFPMOJD2wLIqiKiKmPo3WKT0ts1THwKoqYZFA/D/aANqdFgrkAaXl2DF4xfRkRcB6LfRZDtVSGYPYxqSsAyGztS1nEUSbWSlvoTPDvI12exsmMbREvULkSY5nkDlSEYh1djWWevsC79q6hIRNdqTOoxt6OTL3NQuDTGVk+V1JjR7vGUojUK5pABjM4LE5gwGq1OcMOvG7d1ub4T0UwIChfp++aJch1hsAWCGBZBEayg4Ut5oAoGfzm0Jm7q1iFT6jmRJDQV8Ukx3mCVOnE1X8Boi3s5Dcf3VxAZwtYk5pdVxlHfgiBMOTivzciCmsIfGdDGcLHYonIqOCmOTNMwiYYObUhvqvjw8N9YaswHIYf56uQNX4Z9BVFYH8N37zS1e9XFtGsVhEn+VmHvFpierOBmuDiZGk0hL0YcQZN0IRVBoiNFTIo8MiDywJMpmcKzV4WERu6CwKCch10BHDSoxy3CMIZrGiBXXSMAiDnDI1izGFPcZlSqqpGorkHP8yjClJabmx9ivc05WKOIkWLesbimIanrTU+zMQLmjr6ZHnYjUKap2HIzZaoQjknIVZ4unmoIUmary7pyxCS8p/iTUsvU8A8KwaSlNam7SF8qViHZb4LXrtjfD7BENEEqPfhI0iT8QLRDlKvcNpDqgiUKhifd53pVgQFrS4WWaaUiY25PwSr/gNb5iYhcYvvBwoX106FDtEz6hMnBrFAh5KkIpAjUHTGI8g72zDNuOJSOnl+BM26U63gwJ5nENfqKEMB8z3EnsNRa2UZRhFWIMzhGzCGC0A9KKbWCqAYUWgzo8lIsHMmVERiASfQQy2aKRHoatuxuLe5DThM4KpVHryjXNKKEWEJnKVBgJViY4JCUhPyqqsxI7BCjUUyUXgJve8NwftsFeVn6uhLOLwaq6IhaHIaFaEit4M4yp74rNBXRSuSes0IZQiK5wLYp0ZYa3H9FbU1In4/VJRJ23QciiNDHfJZK4zQxgl57bQNatx9k/SqEqnDZ8RcsWIYefPQOu0hUIyS0thZ1mnsxV+Hb7elHoSZxvtIaOtVJ0tLzVvHDBFMNp0TiqCUrEQPjdrb30adGCjiLKIJYcl6IzBl0pHuIVMJIaqjbgicIzyxKKYVpvPBZYhGEeNDSMaWSiTukWvCScXT0QWKtNDHDLhZxQbocXV7zQTp0W5IU2Vn1cuX0CP5WfiWjtK0troi3YMjIWJtBFF0JeWumi2dhUFS2pct9OU3oeyQY+0ZKgARnRbzLIWZG6hUGmAFepCibpR3oYiOX2ZK6rEN2gjZIskhzkspf4gKNRRAe+DpOaK5G8nN/NiNb7e+kZ5ai+bkD65wgsYZcgsiaYqw4Eyp2JdDZVYMU5CxZHQkmF4spukEs9C780LwoWtsK3UVaZ6u8agHU4nk2CFGkpwAN9FEqSXnKssoSCqkb3OFiByHsPJ4pzbCgdISVxNQn+eOPCKoN/nFlipGHF2Z+0AaLod9KgcKgI/5PezEY/A7K2jh3w4jUqiZzVQnBAakj3nTVFIpcagCQHfZIBwsVX8Ftx8POxA+Qp8GKAJA13kBmSWG2FNgN3ejFMGlY1gHOPBpjKMGZgVaqteOXErYAirtIcVASltBYCr69BY9G3YjCt5uZEwNWSjfCdzEO8OC3APLD9FEUTJxeUwLNJtb4ebSwAj9EJkV9eqkswHABIeQmcreQxn3ot617NCHQYY7z/EonUgPQJXdJKtDBseCiMsCdLCrx7hiiBssR1zCJQw0VBI7Wq6dapht86m9A5UQZHVHlbNC8MmYpBNl1P02nCsGrrK1Lk4ESKSabrtcDpZiIIcLJ98DzkK66xUi4zDILbelDoPr4WeFffeZTuLnc4Y3ykOmCIIJ9OET/Udrgiq5cgjcGZs22y6fIpRrbbEE48hrVMN2XBZuf46tmkpZmUAwKCwglqQvJDbW1wRSBdfpeZxhcMbigWMIfADrLJNDASdNfrehDaqsCZY6044PrERWlzd1lZ4w4esFIGgKxNdEU1VzV9MgmyidxkrISUvKSkNAMagyS11gf5VTA4rBl1cMTgrI+rWOd4jKHgtdM24oh5Y0zdOuVGWao0wLNJtqWwRM3yv2+TtyhtH4p3YCw1BQd64lCinoPauD9tfbA3RiPn3YW2+qVYOHY9/H3N60lOGJZYOHR2aRZ3QrbPPQ2oFcU6yQGowQ02JGqbqymHvMS80eig9SelNDRdYTkIpsOFatZhHLb31gFHMI7BcGz0jrgjChn0TEtMDcX6l6nLY+jvobYcKJ4AZ8whKXgvuEI24Uq7AmXG07DxxMBWBynhweGioUoo8glmnlOV8GwOjjHzORBvl2MUb5vcXnE20zeWRY/ilQ1gKmiPsCImeKH6RBWJqEY2sLAa4wmltraFMDnxBZ1XRNuooDJJvTLmI89WlMHHu2M3w/FWrFlCtwMji7Bmz9X6RdRNb1pGwH0sQ8zyG2goM2rAVF3rWEZ8AUEYX2zkxLMWOQncAVwQUs9A6I60fXLM6tTgqENe9XF0O6X88gTnqebD2ZfiMUFuJb9Bl8XdvO1kRVFk7rNswFYUDNpojgH0FPiNUGnHPg08pmxzrDucl1FfQQRmGpCsnJIvJEaEUEVKThZLTCvEA3m7DYbzBgQxJxa+/WoAlchHCE6iEXVunb5oV0QLDsWrhvR0jDShKLj/UmgOYrUGk1xH9mZZWUauUuCfab0YKh+IeQTVoIyjEvXfODEzXtG8nOJCKIEYfFR6Bqgi8GSdT5f0uXItbDl2qwFRaCwxbaGV3E/3CysgxWPkwcuQnLjC3yTeG8ooY8VCoIw8PzO1xj4CEtRkEWLvwDD+vQ9eNHKdjLaHsTlAESpVjIYy9boXW2fCmaThN9FkuVgYfs7YmQSRlO/kj4SbLJgzvyHltOGbkEfi56kwjFwN3gDIc9Ep8k/XDHI5C64xV4o5W/MZYQAkwnDYGzIKRLyotMLbCdcaVqAgRdq9gmxogM97dpXZIMHSalzEOjutimbUQlLl3I3su9dqbYz0Co7uObWrAsOJyghmamqmtnrtGRbHax2/Q5qAZ63IaxbmnK4K818FlEj20hJeZ1PpBhoBknUe1UoPHjNh0uHGIZnQ04Oaq0WAj9foo36fsNeEMhVfldZ1UOS+95NryERRzJtqogJxWrKBMxgO6zgDLaAHlQyPH2Y0pZQdSEaiLyxlwRVAoRBSuIFdFeYbhNMWgG7rqXWVm7zgWRyPYglc8PHIMq8Zv7Ob6hZHXAMBvc0VQW+WKQC7I9vYmTARh2IExhtaVs/y9R06OHKdnLSV3koTCoa8tx3sgqQlPBabTRJuGrOcZFSg5fOykl482WbWCebgNddHrwFFCNn5uttkHnU58joPfHeMRjHSHjSuCIF+dOk6UBh1uESOaQzzobEdhGlBofRb7a2haowbBodVjCBgh6IzPEWyvX4JFAajGLfy80gplnOGR76+jNcYDpWId5Skcf+pH8xLUdiPxwrVIcVtu1P8HiCavzaIIikEHWxb/TiFTJqGVRRB6oSIpLbu2TgmjdJQZHSrbiiWsgxprwc3Hf7so3DUhH9HbgstM5EWPIlsYh3I4VUBGWD2/deUCLAqA+qj33jEbyA8WBWXZYUxoSDJpyIwUASvUUEEXrj+evSBRVXq58LyCOJZqqTPGtT9rgVVG5ygU6vwmaG+ODwkErcsYMAs10TIiL1ggW+sXRbJYyvHRFzOHDx2/ceQ4Tn4Z9WC6Iqg2DkUWXbcZKrLh0JA1hvHgWdOLogBOg7SpDF+Z2TycI1BDD3yOQySLFeqosN7Q4JdRyHiutXxC9PuPd+v0hzyC4ZYMgGTnTN44+RhN7hmVYgnMUW+q4m7Azo1agcVCAduoguz1sTKaazy3kKvzMJcMv/CNiSmkAS6zNNhEYDdnfQAAIABJREFUNzeqcMwZGD281TOfl+CYFRS8yCOI6mOUc/fiTfSq1RpcZs5U61EJbNjCYwskqy0hBCUTqFLRALNZz7IugkpLfMId+rzRJBsN2fh+gAbrIBgq+JOMqEmFXqazxY0jkQe0jRpyg1ZYR8CUZHFnnbc+sZaOjxynm1tBZTCdejtPHEhFoLJffuCCGMWnKAIq1FAhB51uMqebF2DZ4Wxgx6qGCeZhV33j0nlYFMCoHxs5TnFZVG9ujw8JmPZlbBnLIBF/za/wsE/7ynMjHoG59ihsFLF0ZHQktFtcQZ21k8cu2utwmYlibSWy6PqtoYRnnPrWH0qsBoXZWujmxEAbVqiHM5vVkZ6EeMqkwuxY8zQq1GBRgF53SttmwS8vVJZi3Tpl9W2gfCfH6aFMDthQ6wc53rE/oWOnNWiGuYWQStyLvCmm0EeX/A0MSuMHKzWNJeT64xWBvclbcRSX+HppLHPDYDg0JDfOur8FtziqcMwKf257I5mdZLptdEWydGDVwjyO6nkw5R4q+e0YZdk0DR7nnpbwZHxUY1A5CpeZ4ZjHpDoCw2mix/IwlKr5rlGBOSUcqTLi5DQwrxc3cqTMZqvJW6eU4kq0scx/N2eCAs05TXSU+Q99s4qC11HICVGyuCeuZ2llNIzrFA+h7i8qi7PDmNCQC95sDuVxbZuTL7rd66FK/XDjcJXKQwIDE5YTYwFaF58EAOQOv2jkOI0jYurY5vMjrwHAcv85rOejEdBVEfbprz8Lk+Kex5Hm1/BU/raR+DMAGJXDsChAtz3e0jB762hSHSBCtVbnIYB+K8awibVS9trhLAIJVqijRAP4U8ZIyj75kolht7cw3AlSJo/DsYEKuyL83JQpcv2w6G0JXapEhWGxil/+uCVoulSKszhk/YbdTB6PWvaa6Obi/YlYjyfaAyZ9gQDNZhNHsAXWODn2OM38EVT64z3D7uWnAACHrr0ZAHD4EFcmbmcjliwGCzAYuFhlG/Aqo4aH9Ci6W+NDkYBs9cyvrRtrIBeMzL8AeBJbVdQAYnz9JAz6NnLkwyjV4yGehCJG02miMxSO7JuVqXmpqK30cng9u62t2DqQxmFb0HfNajyMu7TMf+/ATrbUi4NNdHORIcF7drUR+LKyOCoo62+c48c9ev3IcfzyESyx5kwtOuaFA6kI1A1t01jB5yuvj73NLE9vBSubyKEoG23Vw7g1sShkwxgLFcGh6148cpyjx6+HzQpgG6Ol9p7n47h7Dt16pEBOXH8jPGZgsM4Tw3ID8Hot3OA+ja2VsZM+UVjiIaj1S+fGvp53ttA2+XeRFh0G7dDCZIjTXstBG94Q9S2aiTu590tBdDmVgz+6rc2wjTYgFSmH3e2iRIOQnw8Iyw7TWwvIeG6xusT544P4RqPyultb4wu9crXJoTuAT/EaiORiuVxGi5Vh2FcgQylMJAkvP/soP+bRW8Yep1e5Doe9i+OFbD6NAbNw6DhfC6VSCWtYQq59nv9eFHkEl88/jTz5YCsnRw4jPYokDxQQrZ5NOYta8fKYUsAoLlAQMNRZZ8ST6hnlMGeWhHZT9kNqwFZDPAk5gtygGWOPAYBrVqbO+JU01kJ1OWJbtTfjNGIBWyiCXC3uTeULBWyiBqM73mMDgCVvDb1CRN12Rf+okD5KZqjX3O0LCBjh0BjvnapHYBJDZzvZa5s3DpgikJtZtLgM5sJQen0AESOjP6FEXlZumqJ4hhVqqDLef1yt+AVjCC4/BpeZOHrdrSPHMU0DF81rUGifHXntuTMPoUE2rOPfFD5XLBRwwTiO4vbj/PBC4Ww8fi9MYlh68R1jz7d0mOcN2pdHFQ7AWU3SqgUQ3pgha0ixZhhjWGZN+KX4zRK1Fpi8QZcCG26uFhYG9dpbI8li6RHInjmqpZ4P22pPluP2IiYUrz2Ikp78O0Veji08gkI9HhIoLQsLejNhgwbQYC34Rf45IsKGcQi57mWASc+QN7fbOPswAGD5upeMPQ6t3IgldLC9Mep9lJpP4ZJ5PObtXTGPodJ9XngEslaB4dmnuMI5PMbwqIq6AreVrAiK4voAXBEU4AKeIwwcuWWI3kl2m8/3LcUVgWMq7bgTIOdMWOVlXiQoa0rU2pWhOo/eUDjSzdWmzviVtNNyfSXsjdVvb8UKC+U6cFqSrj0avts2VpDvj/cMfVHD41WimD8VG6igA1I8ApksLraexhVjFUYuP3Is6bVtr42Mfs8MB0sRIN4SGABM5sIauhjSCnTbyRo5nGYlLEgqNmASQ99uxit+EaC+/Q08Z90Qi22q2C5ei+X+qKV+4XN/DgA4ceoHY89fKd+Ma3pn+NGJbwy5C/cBAG599feMlbF0zU0AAGftmfGv+xsYFNWagAosN1IETHGfW60mKuQgKMfdZznZqz/FUq8ynluRTIy+EhqSiVV5iaRSsZQmbbE6hwmQbcIr9SUMrKi1dizchXib8Eo9/p0iWuf4tdDr9VGnLpiiFFu5VZSdKwALOGOIC4X73Gk4yOHaW1899ljFa14KALjyxJdjz7ueh5t6D2Nt6RWx59ula7HsnIsVMIIFsM+eBgBcd+uod9gQ3yeYsLYrQRu+UAThsJmuSEoPsZNkwSOV4wrUUWf2JkB6jrnyEm/zrPTQAkYrcYteG84Yem95yuwDLyx+PBQWZg7srbH0UekpVVeOjBynZa2g7Iw3DtfXLqJALtCIFEG+eggmGAai3TSDGfrVh7tP43JxNFQMAKUVfo06a+O99yxwQBWB2AgCBov5sHKF2HuKogjH7yS7gQNR+FMQWf8obr0lQhuRq359/3Fcqd6WeCyn8SIc8y+Fk5FYEOBL//Xd+PZz/xkP1L4bR66PW3b9ldtwDPzcZAjqJucRPG9eh1JjNEEIAIeP3wCH5cC2zo681rFtHGUb8BoR26hr1lF0t6H2SZH3ZHuDW8dmLV7BHG7Qk7jWfoAKumCFehiGGdjbcetMUQR9pX+LhJzoNG32gaQbVusrPNYtqarh1hxZnJ7g71cPxel8S6v8b6893oJeuyzYH7Vo4+iXjmDZvRJ2oWXE6YnLWw/hXP4mGEPrTeLoS78TASO0n/x89CRj+Npf/ToaZCP/kjfG3u+u3IojbAOloBOGhsAYlq/chwvmCVhLJzCMWqWELVYFEmiqnudjibVCJS8VQa+1KXoayZyUUKBNvg6toZCaZ9VQnFK4JhOvheoyBlYV+XGFXjF672g4kuWrPHSVUJAJICx+rDZWlLWjrDmFGu01uRU+jnnXyx9G3RuvCC4/zb29yrHI688f4vF/d/1pIYd71f1+D9f659BfGfXYAGBFeIzdi4+NfT0LHDhFoN787b6LHDzk8nGPoCTi6Uig8gHAoCX4/cLVlskle/MiiEUx28HaM2igAzoxPnYPAPlbvgc58vHM/Z9Ev9/DF//wp/CtT/7feKB6B279mf8y8v7c0WgBSUVQxADnS6OhJ4liPodLtAqzOZqUvnj2URjEYK1Gsetu8Qjq7np0U5IVWmqdTUlljFtNBdXaSoDsdU/FRjQNrrsdteqVikCyeYR3UVA2GknRlBWwiXCaGDAL+WJFVGQr7BcmQjZSKXT49WwcjiuCeq2GDisB9viQQPMiv8llK3IAGDRuxGFswXQ7wiMgMH+Am70zaB1OXgfHjx3D0+ZJNJ7/DMAYmptX8I0/+J9x6onfwQPF2/Gy7/nR2PsL13IP4Sb/mXAdBL6HFztfx+WVU2NlGAahaSzB6I7/Ptvbm3xeQpWHRmQ+ptfmBXJhrkg20RPkg8JQczs+nGaypS5j9+X6ElwrKhKkMWweALEuvNGP0ECOfHjOBKXTb8IHwSjUUAzHvW6DlNYPUqbRvoA2yuFUPxVeeRVLbGus0mk/9xAA4Mgt0fWtHeeJfUMYX1IRPPHg55EjH6VrXzFyHAA4fs31aLEy2NoTyd9pzjhwiiAAhRvadlcqgnjIprq0yvuPdJNzBEy41vXD3Ooqi42gfeVsjM6XuyIWyK23Jx7rxbe/HjYrovGP78Jzv3UHvn3ro7j/xI/jVb/wMVSqtZH3y4QfECkCAPDH0AVVbJeuRcMeDQ2tPXovAODQi6KF6ZWP4lCwCQScuRAoLXS7W1wRlJfjHOiS8KS8TnKyWMaFzVIjatbXbYbVlzI0JFMGA9k8bSlSOuHs5t5kVorptNAWfPio9sCHTOIGiHrMmPYVdFCCVYr/3kSELWMJZsJa6K7x2o36sciCzB/lFl29/RQYcUVwvP8USjSAcX3yOgCAi7f+GG52n8Clf38byr9/G168+Q/45PF/g5f8wj0jVcLLx29Sz5T/s/44GmSDrv+2RBnNwjHUeuPjzy1R2Cg9nFxZDfeN0lTl/F91NCrAC/GqUyx1dRCQOtUs8g4jRd3v91GlXkjOCCEbJE7oAWT2t9BGFTCM2OApSVBgSh1BsXsJm+Z4em+ucRwFuGhujnpTuXNfwjZqOKRck+M3cIOtbD8n5PB7qPPVD2PALNz8un8+Vo5hGrhgXYdye3w+LwvMRREQ0RuJ6HEiOkNE7xzzeoGI/lK8/mUiOqm89svi+ceJ6PvmcT6TwECh9h+c/SIsCnjVoYJ8PodtVGH1kxWB2b2CJqugUuZ86yWxAJyN5wWLg2/QR7pnEDDCiYS4MADUq1U88LJ3YsVfQ93fxtdf94d4zf/yRzBMc+z7y8qmyBRFIIuFktA99HJc7z+HXifid29tbWHpsb/AFTqE625TNqn6NciTB7L5oueMB8mB5opg5Wg89NA4xENFfieZahnG/CuNsDEX67eUxlxGLFnsC7peZTnaaCpyCM4Unnqsm2hxCQYxOHZLid1HPWby/XVs02glLgC0zBWUnfGsIXfzLADg8Iko3rt0w8sAAMcGz4YewTHGf5PyiW8aOYaKb/sXP4d7j/0Y1mgVnz/yVjz/I/fg+37mvSjmcyPvrSxFv4k0CFY6PHe0enPyeutUTuK4f37sJm2La5sXTfBkTN2xt+IVzGLjHIgWKEurQ+yXQh0WBXCdZK8gUEI2TLKTGIsK/iiq85BdeI1y/BqZSrfXJOQGTXQM/r5GtYwuK3DZY2jEtcFltPOj+QEAKB7hFv7ls4/Enr9y8Xm8rP0FPLXynWG9DwDUl1fRRhmrA650A7JggOH42hfwROmbY5PJhrFdfRGu6T8VG9Jjt5u47+4PTC2k3AlGCedXCSIyAbwPwOsBnANwPxHdzRhTf62fBrDFGLuZiN4C4L0A/hURvRTAWwC8DMA1AD5NRLeyeH/YuULtL1N67CMAgMNrXx55X9NowOpP6A3fW8cWLaEhmEirx66BwywE28/HknfXsou4YBzFNeXJm/Qd//Ln0bZ/BqulEo4ZNPG9ZWWWqkwWA1Ef+CQUb74D5rkP4al7P4ybv+MtuP8jv49bnvhPeCk28MCrfwNHlEWcO3wj8CRgbkXxTUPmVkRYbHl1OIzSQI/lwSZwrWWSvVhdQalYFDMMmjGaqpojYD05xyGy0sxcHl1WmDo8KOdGLBPJaOo0N0I2D6MoDFF01tEZU4kLAM3y9bit/cWxr+XXvoHLtIqjygCYEy96GVxmYplafBSpgiOr461NCSuXx7f/m/dNfI9EVVEEMhR5rf8cpyUKZTQO3spNqGz0YW88j4oS0gKAjqhOXzrCn5eN3WRNBp+tEDF7WPsSAkZoDK0FmVvoNDexXIxTPsNz7m3DYwaK5XqoOAa9thKyidg8neYGVjGai8gpoaskFL0memJYU8EycZnKYP1tsGrkEchq9hV/DWdL4/N5S9fxZH7n/KPAt3BSRr/fw+U7fxQvhodj3//v4h8gwrp1HDd6Twk5JgxiuCE4h/sOvyHxfAGAnfwOLH3tE3j6a59FaeUEnr3nd/DSyx/D7ejiiWtfjFtf/d0TP3+1mIdHcDuAM4yxpxljAwB3AXjz0HveDODPxOMPA/heIiLx/F2MMYcx9gyAM+J4mUFuNACwBb5J2G/83ZH32dYS8oNkRVDsxzeOQi6Hi8ZR5JrPxnjdALcoZ0GtUoY5RQkAQL4cNfhSPQLKV8a9PcTLXvcmnMNR3HTvO2G/9zZ8xxPvQadwFM/90N/gVW/6t7H3rtzALdelDl/EahVuvv0c1mgFlCvFPmMYhG2qw5igQGUiuVhbDqeNGU5U8RyQISqLuSzqbaOLwgjjqkulqZPDCkr/fsk66rY2hjwCjpq7Aac4fpN2l27CCtsO+/RLBH6AE/Y3cKkcT/pVyyWcIx42Y0BMTr0xFNpIgWqlBodxT0GugxIG6FIRxfL4zRcAatfwje65xx8cec1f4x7FERHWOLrKreNua1MYOFyOGlLbpjrIjHssYSuLCVRi09lGS7RkoDHU40Dp0SST0vmhXISckOdMCEeWvBb6CtvINuow+xFrSE4O27h8DofRRHBofK7t2Mnb4LAc/Is8MRz4AU7/8U/j5e5DeOy1/wEnbnnlyGc65UhBSq/NIIagOlrsp+LF3/EjaLEyDn3sX2P1ztfi1KW/xJO11+LxH/gwbn3Vd0387E4wD0VwAoCagTwnnhv7HsaYB6AJ4NCMn50rVNZQr9+Hxww0bv7Wkfc5+WWU3GRWStndRH+oidxW4TrUes+NKALHSr4pd4JSJVrUao6ArPFsFIl8zkL7B/4jnsndgqcqr8Zj3/unuPmd9+L6V37vyHuvvfE2DJiJazxuIapueq13Duu50WZZANAxGsg507ucyvxAl3h3y2h4B7eapEdgOVsj4/wAOeJzMj2xFHTgiorXqHvkNi+MErF7YgEGjoNj7Arc2mhxDwBYIjl/STBDJL722Y/gBK4guPWNI5+xRU1GAANMmUFB+dHvslMYBqElQl+qZ+jS5HXwom/i+YPmk/eOvJZrnsUaVsIxjfWlFfRZDqx1gXtS4rvIZH7BWcP2GEPHUn/vBFiDFmxRwWyq7bvDZLFSidviiqA4VOhVEp7ioJUcjqwEbbhKC4xm7ghqgt4LyCRugLXHuddXOPmasccpFos4k78Vy+unwRjDvR/6RdzR+jgeuOHteMU/e8fYz3TL0ZpiRnSNzOrkfN7K4SN46rv+EM8WX4LT1/4E1t9+P77lFz+KF7/m9SMzTeaB1KEhXSCidwB4BwBcf/1oWfasYEqy2On34FIOJXNUH3rFQ6jboxYTwKspl4MtXCrHY4lu4ySOXXoANpViG/RgzoqgWFAscWVxJdESVbzkNd8LvGZ04x+VUcAzxjW4kXE9zciEIXbnVfcinl167djPda3GxHbXfjeKCwNA36gg57bDjUWG1GQiLzdojsxglp/LT5jsBgCVoIN1oQiirqCbAIIwdk8I8PxTD+Mm8oCjLx17nBO3vBK4H7jyxH247uXfCcYYTn/6r3HTF34R54zjeOkbfmrkM65ZBlz5ffiN20MepTHtP9KgQ1Wssu2YIvCM0SIlFUurx3HGuBGHzv9d7HnGGI61H8bl8s2QvhEZJi6ax1FqPwtSeltJD646WIedH+2qm5fzhCdQifNuMwzdyUH0/Xa89QOGktKNw3HKckN05U0qkPM8H4fYFi4qDR/t4jGcbD6OC0G84td++j74jHDNS0YNQ4nta74Hr3v2j/DQb74Bd/Tuw+mVH8S3vO23Et/vVSO7VvXeS6XyuLfH8Kp/8sPAP/nhqe+bB+bhEZz//9s78zi5qjrRf3/33tp67066OwnZN2KCEGKzix8XkOUjoA4qqDO44cdRnjIqIwzv+RyXGbeZccbnPAcX4IHOc2XgOfh4LqAzikBYE8CQhUD2pZPupNPpparO++OeW3Wqq6o73V2nuqDO9/PppOrWrbqnTt17f+e3A2blpPl6W8l9RCQAWoHeE3wvAEqpW5RSPUqpns4J7KzjoYxIkZHh46Sl2AkHMNSyhHaOMlKiBtCOXTtpkuPE2wtXxUHnMhpkmE76MKc2E0z8o0+GIMifUKbm4U2gEUyW/nj+ostqjaCvv48uDqHai+OsAYYT7TSmy2tSkc0/pVd2w0ET8cwASmdfqjEtF5PpIwwFxYJgMGgdt6x2WIN+EKVt1VGk0eixsA9zJAhQsGvTYwCctLJ0yOXCZaewm04Sm+9l4+9/zoa/eQ1n/O5aBr0mvHf9mESqWNCng1BYZ8TPhVwOUtnzAEKBCJD1TlwjANi37EpWjG7i4e9+kg2/vYsnHvgJv7/zcyxSuxlaXLhQ6E8tou34i2FY9BhncWvmUEmTWkpHto2Mk8GcSB/JJYjFc0mCh8k18TH+zfaH0UxtXYVF2tpndTGi/FxHtrH0HtxHUkZRLfkbste2IOwFoH1MUVhnbN8TvOgvpLOjvCl39eUf4wU5iVMGH+GhOVez7sO3FTiIx+IZ/pECQdBQ+XNhOlRCEDwCrBCRJSISJ3T+3jNmn3uAa/TjK4Ffq9D4dw9wlY4qWgKsAB6uwJjKovBCu2M2ywVH76ZZlTEvtIdRIP37S8TdP/kLANpOPr9ge6IrH4evJmGymQ5myYET0Qgmg9lIJwp92/nMHwBIzX9lyfeMJmeVjbUGiB0/GFYCjYU2/6gHc2HIILlqpA2ZI4zEiu3qQ/F2mscRBP1Hj5KQ0ZyTOJWrkR9myJoaQfbFhxklYM6y0t/J8z22LXo7pww9yin/7yrmjb7I46s/xdwbH2festJRQJlA36AN01BUqrqSpP1wHiejEQD0XPkJHk+dzZkvfotX/vo9rH3gfZy39e/YHKzglZddV7BvpmMpc7N7Udl0QeLa0aP9dHOIbNvios9vnxPesNN95cskNGQGyMSj/sPh7zNyLOzjkFHR76N9Rcf26vDewkVBLPA5JG1lawAd3hOGSyeMKp+peWF4b/JwmLCVJcyRWTC0iYMt5Z3sEAqeOTc+Rv91z3LWh76JF5ReSEZ4pt/Oy98TGhrG9+dVm2nrqUqptIhcB9wH+MB3lVJPi8hngfVKqXuA7wB3iMgW4BChsEDv90PgGSANfMRmxBDk8wjUYC/jWdpa2sMTs/dQL+Z6R2WzeBt+xDGSzD+lUBA0dOdDCE1BMLaWUSURL38ierHSJSymitlGLysBgqJ/038AsPC015Z8j2qZT8OBYQb7DxSEuUbEhw/S57UTWcozsWYajh1jSEUagS7Wlw2jONqyffSmildow/EOWlV5QTDQu5fZhFVXAZp1xnV2sB9RmfwNWinm9z3MlobTeEWs/I36nD/9LA/+WxcxUay5+FpObxzf3BcJgtCUYlEQePo3N863tDfxgiCRSLH2hp+z4/lnOXpgB77vk2xsYenKdfhjb24dy0jsSDM7vY9BbaZTSrHv+adpBmLdxc7V9raOsJhimZV6NqtoV33s06U5zB4Y6DBVs6FP/PgB+rwOSs36Eb+deJny3Ud1s6amrrw5uXN5DzwIc4+FCVtKfGLpQTo4wvbO5SU/xySRSJLoLO4jUAo/kb/hm+bieLzy58J0qIjBUil1L3DvmG2fNh4PAW8r894vAF+oxDhOhKg2fH//YdqAgUR3yZNrbpcuM7H1N3DuZQAM736a0e9cwpmZfh5ffh2nj6lR1NRihLZ5psmm8oJgWMXC2iaGRuCXKGA1HVRs7EmsaNrzIC/4C1k0u7SzOD5rEWyF3l1bSgqC5HAvx4K8gEk3dNF+uI+d6dHwOFHCUjbLkYEB2mQA1VTiomuYRZIRMkMD+CXCE4/rMst+1MglmaRPNeEf20smG91oPGLpAeZnX+Tx7kvGnQs/CDjnyuvH3cdE6ZWgEi9nGjoRk81kyWiNQBmrzcwJCAIA8TwWLFsDy8ZfBcusMEfmJPbznIS/nZCld/tGlgNtC4rf73nCIa+DYLB0KYsj/YdpkyEyukhbFFOfPZ7vk23m/DQPl4/vH4x10DxaOmR5aPdGAOYtOy23bd7ikxlQKeZLGAadFZ8udRAE4h1T9z+WIkgaJiDDfBck7FkJpkIdZhZ7eCrNzn2h8+n5nv9acr/OeYsBWLP1WwCk7/8iiVvOpSnTz3OzL2TtOz9X9J7G5nxkghLTZFPZlTpAOrph+vnPDsoUtZsyxgpZiU+g0iwf2si+jtJRFQDN3aHvoG/P1pKvN6UPMZzIOxf9jgX4oji6L4xOiuq+ZFWWQ3tfCPcpUTMnavF5uEyLz8FD4fao17OIsN/vJjmwg3yGrNA5uhtPFPG5pR3FU0YLgkClc4Ig7Y1vRpgKkSBAfF02A7J+ZW8yCaMSZ67KqYIjLzxFFmHRitImtSNBB6nh0tE8vXu2A+T8bC0tbQypGDKwTx8ncrKHGkFXZi8DDaWjuoYSnbSVqQGUOPgMe6Qr19ENwPN9dgf5cyorfljqnLCbXSUJyvgIgnG0z5mg7gRB2CUow5FnHwBgwZzSqwy/uYtDng7x+kwrwW/+FghNJCs//MOSDqJUQ3PuYswYdlobGkEuesPwP1ReEORXM1nxmZPdH2ZhzytfLydqnDJ88IWi17KZDN3ZA4wYkRRNnYsAOKLLYyujqNlR7Z9JdRRfnAmd+dp/sHTG7/HD4fZZ3Xnb8JHkPFqHdyMqE/qKREgxVDCOihEPbwBxNZIL9zPt+JUiq53SyvNz/pWMX9nzoLElr8GZzvxFvb/l+dQavDJmjuPxTppHS5ts+vaF50dzZ7gCDwKf/d5sgoHdOs8DbboLs99n0Y9qK/0bZdoWM4s+jo/JLh4aGmLlwHr2tBbH9x8zcoDyJbwhXmLRMR1ihmlIDI1gbFmbmabuBEFafPxsmnM3fxmAtmXl89eaZKhom/ff9heYfUzE8xiUSFU3LnoLgiDSCMTQNsZWUZ0uUmDfDPAkXJ3FWssnw8yePZejKoX0bil6be/u7aRkpCDiqEOXxx7pDW8M0Y0mm1Uc112cWrqK1fXGjlAQDBwqrRGo/h1kxmS8DjcvpjuzFz87QsYw2QDMmle6JPBU8fTcxRnJd6uzoBEoLQgwopOyfmXPt5ZW46apf5/M4RdYyQscXXJp2fcNNc37lAHwAAAZm0lEQVSnM7OfbImOdcP7Qvu82TOhP9ZFw9Bena/g5Zz5e7eGYdypOaUTvVJzwiCNXduezm1T2SyP3foXtMkAyXVXFR8/XizcIOy3XEniqTKCIOEEwYySIWDxUNi4Y198YUGLyrGke64t3HDdowWe/1IcJ/yBTY1AKqyqg6ERGJ/tV3iV4cUN+6aRxJJsKY4bz73H99geX0Hr4Y1Fr+3ZEiZktc3Pp/B36ho9Xn940zejUkZ1xEl7d/FKsLUr3DbaW7rFZ0P/NvZ4cwoitvy5pxCXDO3HthVkFmeV0FzG5zFVIkEQI03kLM7aEASxvCCINAJVYY2gxdQI9Hk360gYcdP5ivNLvgcgmLOauGTYu/3potdk/7McJUWr8dseS53E7NE9uqZRFN6rOLpFF0VcVbrp0uwloWmqf/uTuW0Pf/fjnLvv+6yfdQWrz39r0Xsi/1daFS4ImhorG9YZN4oYKmNB6FswF0+HOhQEPp26HWDvnNeMu2/DxZ/JP3ntX8HsiSMKhiKNwFiVyQQhZlMhUmfNkNFYvLICxwx9M8skNDaNX9Oov+NUFo5uY3RMaeChLb8JE3bWnJfbFmto5SgNNBwPV/aRHTWrMsQObeIwLaRaioV1x9xFoU25b3vJMXQMbuNww+KCbW1Lw0JsS9NbC8I6h4mV7PM8HYJEaBqKkc5n41oQBOgcFc8QOKrC4cpBLMYxpT9Tf5cl2e0AtC0uXUoZoHVR+NqBrcWJmZ19T7ArsbzAxJqevYoO+kkMH8x3dkMR7HyIPXQyd0FprW3u8rUcVSnUjrBm2O9+8FXO2nkrD7Vfxqs+clvpTFxt9sziFbyeSlXWdm+W+igQ0BYjCadC/QkCw2HT0DhBLK8IfHILXL8RXvupE/r84ZKCwIJGoL+H+dmVtjsGhmnIvFiaJ1g1JZedR1wybHnwZ7lt/X2HOHnXT3kueSpNY5rnHAq6mJUOnYS5VXpW0X30GV5Mrix5ISdiMXZ5c4j3F5fV3v3CFpaoHRyfU5ggtmjFqYwqH1/CxKh8NE/lb9CxpGFWi76TDUGg7fOSGcn18lJB5VebgxL+5pGgjpHhsGqmsbl87aQFK9eSUcLwzkJB8Mi9t7E8+zwDiwsLrzXo3JT2Y1tQuTLhWeYPbGBPy6lImdIKsSDgucYeVh78BQ/987Wc9+zn2JA6g1f9+XfKJ3vFzHPYKAFS4bkrEATmfcCClWA61J8gMCJmG08kqaOpE9oWTLyfZlRCASAFzuLK/+iR49EMGY0nKruaCVKmRmAKgvHn7ZTXvJVeWuHBb6CyGbY8/Qh7vn4xbeoIwRv/umj//lgX8wijS6IMWTU8wIL0C/S1lY5IAehPzqd5sNg0tP3evwNg4auvLtieSiU5LKE2E9q6dfMgLAgCI9s4163OgiCIzHehINC/UVD5iJRcDoRhTx/0xj8PmpqaeTa+hu4996OyWXZs28QfvnolZzz8Mbb6y1hz2ccK9m/TXcHaM73hdxEhNbSPLg6RnT9+LcrU6z9BQg1z1v4fsr7tYlZd/38IxtGQI203LukCHwEV9q8kDV9Agd/Qt7AomAZ1KAjyGsFEK9spfb4WAMr4oW1oBL0r3g5Ae0feXh+vsGmonEYQTJAMk0wmefbkj/CK4Sfo/exSlv/oAk7K7OCZV/8jK15V3FN5INEddsUiv+Ls2/Zo2MlsnM5uw23LmJfZXWCC2nDXV0LbcPulzFlSHN/e7+leBnhkxZ5GYDoJc4lEFi5+0YLAy47iEyblVXpVCzAUaQRePkx1bCP5UvSv+BMWZV5g89+eQ/ft53D60QdYv/B9zL/hP0k1FXYBa5kVRuy0czTnI+hIh3kIqa7xzbKre17Hvnffz+Y3/zs91/+A2ASLIi+Rv/YjTWqYWMULuvkFdcyMxxYKx02Hl0zRuUphmobi664eZ8+pkUvvN27+NsJHT7v6r6Hv/fgH8yUyKp2vYJo3TPX5ROyb573jBn5/V4r4tl+wefZq1lx2PafOLp2NmYkXV1Pt3xKWslh6amkHIYC36Fzie+5k8xO/YcVZl7D+p1+j56nP80jibF5x7XdLvmfYb4Js5JTW8f1S+d8nYZRCyK04LQiCKGooUKPERSflW3BEjvhJPW9B7sY54k+8kDrzLdfx0J7Hmdv/KI91X8nyKz5Fz0nLSu7b1j6LYRXoRUH428xSfSCckDN/4Yry/oqx+CUWOaPEsGqwqbGbv0kdCoLwK/epRtraSieoTIesFgRiqJhehTN+wwMItC/CO2Q0uK6ww9Nc1VKgPk98uYjnce6ffBT46MQHSphJN+F3aO59iv10MH9R6ZsGwKqzLuH4g3EG/uOf+cOG+zhzx208merhlI/eRaqh9Iow4ydhlAIfQbnCg9Mh2Wh+pyj5r/LnQS401fh9xvaJqASjnpmv4OGTOSGfRywW56yP3QHARDm7QeBzUFqZQ69uJerlQpbbOiud6FUcCGHjPCighgVB3ZmG0to0ZJqIKkkkCMxiVDZDxbwTaGQzVZIFVTVNjaDC0UkJI8RO21E7hnfRmxjfN9Pa3sFjc97G6QO/5eydt/Jo+0Wc/NG7x63smNEraLM8tA2NwHQSRv4cseEjyIYx+mb5Cs+CIIiqqZoaQfYEittNlqO6paQac2tqnjV+I5fJEjMyfsWiibAALaxHa3D9XXsjskx0UZomoop+fk4jMHoJ29AINN44JXCniykIohDIEQLilbajJk0ziu71LIfZkhq/vy/AOR/8Ok/89nU0ts7mjHXF/oexZP3ohublVtEZCzfohgZDuGmBM7aLVyVI6+9zKJ43u9koaZLx81FDymJexIjfEFah1o2DAI6RpLHCi494qnoawX8u/i90je7KCxwCC+EJ06PuNIIIWxpBdFMJCorOWSxDbUmgASQbihPKRi2cwn5DsUbQKMMcT0zcd8Lzfda+7m2sOAEhAKCM+PF8DaDK/z6+7/Hl0bfz/VNvyxVRsxE7Pq/ncj4+8iH8N9ycP7aFypbZwAwfjRLXKv99Rv2oamteUI9Q+eOYiV4RNjRDgFe/5/OsvPbWnEaQrsH1d+2NyDK51YwtjSDqS6qbawAENk1DFrWNeGDOkc2Ye0PzMELshhLjt/Ob2sEijUByWo6V+H7gL78QFix85Iu3AIV+o0qxtKuZv/+bLxVs8yocRgz5DGblBfkMZgumoVGteZjlu22cc6mGYmexdR9BVAvK4uJtqtSfRhCZAyzJQC+qq29cJEHcomnIQqhghJnAo3KrGRsx96ZpyGiwkmgvtfv00OGWWeXlk9cs3NBMRJ8TNvtSmEwU3jsVIk1KRMhgz/mdiUU9mP2cL8LGtVqQ6JUzEVr+fbQAcBpBDZBzdFmSyr7SdfUtNowxiVsuXrXvDV8Li2U9/gBgZ9WUaDTKdxsagViYNz+Ku9e/E0DWgmmoAN1cRaolCBKVz4+Jqql62ZFcprSN7NhMLDxO1tAIbJxzqYbiQIiMJdNQRNQzoiCBrUaovRFZx15JYIDRVWGBq/aTjXo6FkvOppJ2BUH3+e+l67w/zfsILKyakoYgMCs02vCteDpUNciO5CudWljZmog2E76UNYKoiJ5kRvKmIQvfR2lBEJrutBZqwXRX6EexayLMEWXN1+Btt/40gujit6QRnHfpO+l/7dtY2pA/qWyahhKxKtkbI/XZRqhlU2mNoFyd++ng69VyjNF86QfLdV+iLlvVEgSxZOU1gqgkg2RGcqYhGz4PEkZDH22ZtLJoM2oNRb6ijO36PzpqTNVgPkHtiSbrRAk49mRga0PhysKms7hcIS5b2FDTG4zObmY7P6/C/RUgny0tKpObu0pX6xyLr8LyGTbi+0sRs2Aa8vS8+dnhfEN5G/OWa+gznKuwa8V2b45dquMrinKLsjV42629EVkmHylSPWWo0uWhZwSLDrWUuYL1zUS8yt844zlBoHImjkoXGhtLoP0RXpW6UsUtaARZvQDwyRJoU5cNQRD1nw7UaP5atRjNM0KQ7yBnWRD4+jwzCzjWCnVnGsprBNUL4bKycjL447lfpemk1VS+YEYem0lEBaWCjcY/41WPnCoJnSTnkYl8uNZ/H08Xg6t046By2BAEMtQPhLWaouJ2Nkxdfq6PQ950ZyuaZ/AjT4Va5x03AnaioEyixNKXnUYgIh0i8gsR2az/L4r3E5G1IvKgiDwtIk+JyDuM124TkedF5An9V9xctNLMgEYwUVez6bLqjdcy32j2YoPIjGI9xM74XWwkRiUaIkGQxdMmG9u14SOB41sM9TWxYYKS5a9jR7aTras/nNcILJjuombvcUZzfjwrDX2Ahs5FJNvm5I9jWRBEpqFadBZPd0Q3Ar9SSq0AfqWfj2UQ+DOl1BrgYuBrImJ2s7hBKbVW/xW3MqowcRXWZlk9uN72oV5W5JzslkMtfWOVaSPaKqFNQx4KyeoQUgs3NJOoSmeginv3WsFC2ZEzT1nFgQ88wlsueqNR7tqeDydGxnDm243miQpRWhcEL2Nn8RXA7frx7cCbx+6glHpOKbVZP94N7Acmrh1gibbMoZk69EscO43Ri45itMcMLGTIplrD/g2NDOVyPmz0lDZZvjbs6zu/20KmtEH/xd9gaOXl1j5/3cJ2PE8IJNQIbNTQipvlRqKqrZa10ChhzUZnN5ORrP1AlakyXUHQrZTaox/vBbrH21lEzgTiwFZj8xe0yegfRMS6V3XEctLIy5VcL2HbkRVGxEvMQrJcoEuPH0gszEfzWHbmN138GXj3TwgWn2v1OK1nv5vkO++wegwTG6aukoLA8uIjWqFbCYc1SA/2AdCbqGxJ7UowoSAQkV+KyMYSf1eY+ymlFPnE3VKfMxe4A3ivUlEVLm4CVgFnAB1A2cbAIvJBEVkvIusPHDgw8Tcr/0FTf+8k2b3qveydd2HVjmeTXHx/hXsejMUsSV3p1psA+AED776PuR+5lyDn9LQc1ukHsPwCu8eYAWwI0ESBIKiO7T7K87AdPdZ92kV8M/0mBi74stXjTIUJr2qlVNkzWET2ichcpdQefaPfX2a/FuDfgZuVUn8wPjvSJoZF5Fbgk+OM4xbgFoCenp6yAmciqmmdm3fV16p4NLtEgkBsO9kNR+dE7QanStPys4F8ORAb+Qr1gG9h3lINRre6XCkLy1q89ubbLOkOsLi7jQ99/ntWjzFVpvvN7wGu0Y+vAe4eu4OIxIG7gP+llPrxmNfm6v+F0L+wcZrjmZCc9HdMChXFcluqkzKs9OcbEVYJy3WUYrlELycIpoKNRMlUY6gRPJRdVbU8j6gEiFgWBLXMdL/5F4ELRWQzcIF+joj0iMi39T5vB14DvKdEmOj3RGQDsAGYDXx+muM5AULpf//KmyfYz2GSjQSAJdPa8+/4Nb88/euIKQjidqNFfEJBYKM2Tz3g2zANxWN8Zcm3UVf/gIyqjkYQLQ5rsRhctZiWnq+U6gXeUGL7euAD+vGdwJ1l3v/66Rx/KuR+9HhxYwrHeCj9r52LZdXqU1m1+lS2/v6nuW2JwO6FGWhBYMPEUQ8EFuZNRLjhmrcB8Lu7tQXYdo0mXSacGuwTUC3qUARqe6Blp+fLDS+3arLrZTF9EIHFfsyQNw05jWBqBJZNd1FvD9vlu73RY4A9n9RLgboTBKIdQ74TBJMjEgS2TxljVWa7oF5UA6ix0UL9/jogbrmInlRJEOzyw5DioTnrrB6nlqk7QRBpBIFfv2rgVMimw6xY2wJULJfjMBlQoWmjZe6Kqh3z5USQapp4p2mQ0wgsJ/yd867/zldW3Mmprzrf6nFqmboVBH7gBMFkSAz3ApBOzbZ7oCpGbrx44bf4fONNJJsttMSsA+Kpxol3mgb5Fp92BcHCzhZueNdlxC37pGqZurOPRM7iwJmGJsWu5e9k18HfwJq3Wj2O2O4SZXDJ+WdxyflnVe14Lzf8uF1B4FksbucopO7uhpLTCOruq0+Lt150AU+ueZQzFlpePddxLPdLDttOXK0R+FXq7FbP1N1VFwmCwAmCSeF5wum2hQAgTlOreXoXXVqV4wQ6qkssdFxzFFJ/V52KBEH1TBCOE6eazmLH1Jh1zZ2QGbV+nKgGgB9zgsA2dScIstlQ3WxpqE6TEMckqcESvY4xeL71ZktArpyJ5wI7rFN3pqHZTaG9sbOlfpNHahmnETgiPH0ujI6mZ3gkL3/qThB0pMIVpzMN1ShOEDg0TTq/Y1aH3YY+jjo0DZGN+tQ6QVCLuNIfjogl7/4fqOcuY9bJdhv6OOpQIyCj+8barnHumBIqbjdb1fESIt6AnPKWmR5FXVB/gkA7i50gqE1mzZ4z00NwOOqO+tPDo7A3JwhqklQyjvJijKz7AC6f1OGoDvUnCLJOENQ68umDTgg4HFWk/kxDGecsdjgcDpP6EwROI3A4HI4C6k8QBDqj2GkEDofDAdSjj+B9/xeeuw8s1zh3OByOlwr1Jwg6Tw7/HA6HwwFM0zQkIh0i8gsR2az/L1mnWEQyIvKE/rvH2L5ERB4SkS0i8gMRcYZ7h8PhqDLT9RHcCPxKKbUC+JV+XorjSqm1+u9yY/uXgH9QSi0HDgPvn+Z4HA6HwzFJpisIrgBu149vB958om8UEQFeD/x4Ku93OBwOR2WYriDoVkrt0Y/3At1l9kuKyHoR+YOIRDf7WUCfUiqqMbsTOGma43E4HA7HJJnQWSwivwRKFYC52XyilFIiosp8zCKl1C4RWQr8WkQ2AP2TGaiIfBD4IMDChQsn81aHw+FwjMOEgkApdUG510Rkn4jMVUrtEZG5wP4yn7FL/79NRB4ATgd+ArSJSKC1gvnArnHGcQtwC0BPT085geNwOByOSTJd09A9wDX68TXA3WN3EJF2EUnox7OB84BnlFIKuB+4crz3OxwOh8Mu0xUEXwQuFJHNwAX6OSLSIyLf1vu8AlgvIk8S3vi/qJR6Rr/2KeDjIrKF0GfwnWmOx+FwOByTRMKF+UsLETkAvDDFt88GDlZwODao9THW+vig9sdY6+MDN8ZKUGvjW6SU6hy78SUpCKaDiKxXSvXM9DjGo9bHWOvjg9ofY62PD9wYK0Gtjy+i/orOORwOh6MAJwgcDoejzqlHQXDLTA/gBKj1Mdb6+KD2x1jr4wM3xkpQ6+MD6tBH4HA4HI5C6lEjcDgcDodBXQkCEblYRDbpstflKqXaHsMCEblfRJ4RkadF5GN6e8mS3hLyT3rMT4nIuiqN0xeRx0XkZ/p5yZLhIpLQz7fo1xdXaXxtIvJjEfmjiDwrIufU4Bz+hf6NN4rIv4pIcqbnUUS+KyL7RWSjsW3S8yYi1+j9N4vINaWOVcHxfUX/zk+JyF0i0ma8dpMe3yYRucjYbu1aLzVG47VPiIjSybMzModTQilVF3+AD2wFlgJx4Elg9QyMYy6wTj9uBp4DVgNfBm7U228EvqQfXwr8HBDgbOChKo3z48D3gZ/p5z8ErtKPvwn8uX78YeCb+vFVwA+qNL7bgQ/ox3GgrZbmkLCA4vNAypi/98z0PAKvAdYBG41tk5o3oAPYpv9v14/bLY7vjUCgH3/JGN9qfR0ngCX6+vZtX+ulxqi3LwDuI8xxmj1Tczil7zRTB676F4VzgPuM5zcBN9XAuO4GLgQ2AXP1trnAJv34X4Crjf1z+1kc03zC/hKvB36mT+KDxsWYm0t94p+jHwd6P7E8vlZ9k5Ux22tpDk8CdugLPdDzeFEtzCOweMyNdlLzBlwN/IuxvWC/So9vzGtvAb6nHxdcw9EcVuNaLzVGwpL6pwHbyQuCGZnDyf7Vk2koujAjZrzstVb/TwceonxJ75kY99eAvwSy+vl4JcNz49Ov9+v9bbIEOADcqs1X3xaRRmpoDlVYaPGrwIvAHsJ5eZTamseIyc7bTF5L7yNcYTPOOKo+PhG5AtillHpyzEs1M8bxqCdBUFOISBNhBdbrlVJHzNdUuESYkXAuEXkTsF8p9ehMHP8ECQhV8/+plDodOMaY7ngzOYcQFlskbNy0BJgHNAIXz9R4TpSZnrfxEJGbgTTwvZkei4mINAB/BXx6pscyVepJEOwitOFFjFv22iYiEiMUAt9TSv1Ub94nYSlvpLCkd7XHfR5wuYhsB/43oXnoH9Elw0uMITc+/Xor0GtxfBCunnYqpR7Sz39MKBhqZQ4hLML4vFLqgFJqFPgp4dzW0jxGTHbeqj6fIvIe4E3Au7SwqqXxLSMU+E/q62Y+8JiIzKmhMY5LPQmCR4AVOmojTuiQu6fagxARIayy+qxS6u+Nl8qV9L4H+DMdfXA20G+o8RVHKXWTUmq+Umox4Rz9Win1LsqXDDfHfaXe3+qKUim1F9ghIifrTW8AnqFG5lDzInC2iDTo3zwaY83Mo8Fk5+0+4I0SlphvJ3Tm3mdrcCJyMaGp8nKl1OCYcV+lI66WACuAh6nyta6U2qCU6lJKLdbXzU7CgJC91MgcTshMOSdm4o/Qg/8cYUTBzTM0hlcTqt5PAU/ov0sJ7cG/AjYDvwQ69P4CfEOPeQPQU8WxvpZ81NBSwotsC/AjIKG3J/XzLfr1pVUa21pgvZ7HfyOMvKipOQT+GvgjsBG4gzC6ZUbnEfhXQp/FKOEN6/1TmTdCW/0W/fdey+PbQmhPj66Xbxr736zHtwm4xNhu7VovNcYxr28n7yyu+hxO5c9lFjscDkedU0+mIYfD4XCUwAkCh8PhqHOcIHA4HI46xwkCh8PhqHOcIHA4HI46xwkCh8PhqHOcIHA4HI46xwkCh8PhqHP+P4zYzfm7qhaLAAAAAElFTkSuQmCC\n"
},
"metadata": {
"needs_background": "light"
}
}
]
},
{
"cell_type": "code",
"source": [
"f1 = 52 / (0.5 * sampling_rate)\n",
"f2 = 48 / (0.5 * sampling_rate)\n",
"sos = scipy.signal.butter(4, [f1, f2], btype=\"bandstop\", output=\"sos\")\n",
"zi_coeff = scipy.signal.sosfilt_zi(sos)\n",
"zi = zi_coeff * np.mean(ecg_signal)\n",
"clean_ecg_signal_scipy = scipy.signal.sosfilt(sos, ecg_signal, zi=zi)[0]"
],
"metadata": {
"id": "aGLu3_Op54ay"
},
"execution_count": 6,
"outputs": []
},
{
"cell_type": "code",
"source": [
"plt.plot(ecg_signal)\n",
"plt.plot(clean_ecg_signal_scipy)"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/",
"height": 294
},
"id": "Hju_Gc_QcL68",
"outputId": "847b9402-62d4-48bc-e10f-dbe7159e5995"
},
"execution_count": 7,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"[<matplotlib.lines.Line2D at 0x7fb0fc720310>]"
]
},
"metadata": {},
"execution_count": 7
},
{
"output_type": "display_data",
"data": {
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
],
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAEDCAYAAAAoWo9tAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAARNUlEQVR4nO3df5BdZX3H8fcXloSfJUAiRpIYoIilVAV3VLRVqggp48BYtINjFfwxse3YarWj0HTU1nGmiqPW0YopUm0ngopQKf6IYG2xrVISGyAQIlEElgHZiAQUK8Z8+8c9i5d1k/1xzzn3PMv7NXMn5z7n5Hm++yT3s2efe+7ZyEwkSeXaa9gFSJIGY5BLUuEMckkqnEEuSYUzyCWpcAa5JBVuaEEeERdHxH0RsXkGx54bEeMRsal6vL5v3zkRcVv1OKfZqiWpe2JY15FHxPOBHwP/lJnHT3PsucBoZr5xUvuhwAZgFEhgI/DMzPxRI0VLUgcN7Yw8M68F7u9vi4ijI+IrEbExIr4REU+dppvTgKsz8/4qvK8GVjVUsiR1UtfWyNcCf5qZzwT+Avj7vn1nRcSNEXFZRCyv2o4A7uo7Zqxqk6THjZFhFzAhIg4Engt8LiImmhdWf/4rcElm/iwi3gB8Cnhh+1VKUvd0Jsjp/XTwQGY+Y/KOzPxh39OLgPdV23cDJ/ftWwb8e0P1SVIndWZpJTMfBG6PiJcDRM/Tq+2lfYeeAWypttcDp0bEIRFxCHBq1SZJjxtDOyOPiEvonU0vjogx4J3AK4GPRcRfAfsAlwI3AH8WEWcAO+m9QXouQGbeHxHvBq6vuv2bzHzMG6iSNN8N7fJDSVI9OrO0Ikmam6EsrSxevDhXrlw5jKElqVgbN27cnplLJrcPJchXrlzJhg0bhjG0JBUrIu6Yqt2lFUkqnEEuSYUzyCWpcAa5JBXOIJekwhnkklQ4g1ySCmeQS1Ibtt8Gt1/bSNdduo2tJM1fHxnt/fmuHbV37Rm5JBXOIJekwhnkklQ4g1ySCmeQS1LhDHJJKpxBLkmFM8glqXAGuSQVziCXpMIZ5JJUOINckgpnkEtS4QxySSqcQS5JhTPIJalwBrkkFW7gII+I5RHx9Yi4JSJujog31VGYJGlm6vhVbzuBt2bmtyPiIGBjRFydmbfU0LckaRoDn5Fn5j2Z+e1q+yFgC3DEoP1Kkmam1jXyiFgJnABcN8W+1RGxISI2jI+P1zmsJD2u1RbkEXEg8HngzZn54OT9mbk2M0czc3TJkiV1DStJ3bdrV6Pd1xLkEbEPvRBfl5mX19GnJM0bm9Y12n0dV60E8AlgS2Z+YPCSJGme+fEPGu2+jjPy5wGvAl4YEZuqx+k19CtJ80Nmo90PfPlhZv4nEDXUIknzVLNB7ic7JalpDZ+RG+SS1DiDXJLK5hm5JJXOIJeksnlGLkmlM8glqWyekUtS6QxySSqbZ+SSVDqDXJLK5hm5JJXOIJeksnlGLkmF++8PN9q9QS5JhTPIJalwBrkkFc4gl6TCGeSSVDiDXJIKZ5BLUuEMckkqnEEuSYWrJcgj4uKIuC8iNtfRnyRp5uo6I/8ksKqmviRJs1BLkGfmtcD9dfQlSZqd1tbII2J1RGyIiA3j4+NtDStJw/Wd9Y0P0VqQZ+bazBzNzNElS5a0NawkDdcVb2h8CK9akaQm/eLnjQ9hkEtSk37xSOND1HX54SXAN4FjI2IsIl5XR7+SVLwWgnykjk4y8xV19CNJmj2XViSpcAa5JBXOIJekwhnkklQ4g1ySmrJjrJVhDHJJasrYhlaGMcglqSk/aOfO3ga5JDXl2gtaGcYgl6TCGeSS1IT/eN/U7Q/cVftQBrkkNeG/Pjx1+/attQ9lkEtS3Xb9Ah55qLXhDHJJqtv4ra0OZ5BLUt3u/NYedkbtwxnkklS3L76l1eEMckmq0//t2PP+8Ixckrrta+9ufUiDXJLqdP0/tD6kQS5JdXn4/hkc5NKKJHXXxauGMqxBLkl12DE2s09t+manJHXUB39zaEMb5JI0qCFcqdKvliCPiFURsTUitkXEeXX0KUlFuOky+Mb7Z/EX6l9aGRm0g4jYG/go8GJgDLg+Iq7MzFsG7VuSOmt8K3z0WcOuAqghyIFnAdsy83sAEXEpcCZQe5BvuW49D915U93dDk9mK8ME7YxDS+O09vW09O/Tpvk3d82PEyRL7t/IwT+5nf1/eg8ju342UH/bf/IIi2uqbUIdQX4E0H+n9DHg2ZMPiojVwGqAFStWzGmgB6+/lGdvv3xOf1eSumDsRz/tZJDPSGauBdYCjI6Ozunb6FP/8P2MP/yuOsvavQYuERqmbO3raWmc+fb1tDYO82/umv56Mtl7xx3s9bMH2OvBu1kw9k0W3n4Ne/30h3Pq7rilv1ZzgfUE+d3A8r7ny6q22h286DBYdFgTXUvS7j3h8L4nr3nsvgfugg8dP+OuFozsXU9Nfeq4auV64JiIODIiFgBnA1fW0K8kdd+i5bDm3qGWMHCQZ+ZO4I3AemAL8NnMvHnQfiWpGPvsB6e8a2bHdvWTnZn5pcx8SmYenZnvqaNPSSrK8948tKH9ZKck1SECfusPZnJg7UMb5JJUlxNfPZRhDXJJqsuTnzuUYQ1ySarLXjO4tLCrb3ZKkirHnt76kAa5JNXp5Z+c5gDPyCWp20YWtj6kQS5JhTPIJalui5+y+32+2SlJBVhxUqvDGeSSVLc9nnV7Ri5J3bfvolaHM8glqW4nt/s76A1ySarbPvvtfp9vdkqSJjPIJalVnpFLkiYxyCWpCS94e2tDGeSS1ITjz5q63Tc7JUmTGeSS1IjdnXl7Ri5JmsQgl6TCDRTkEfHyiLg5InZFxGhdRUlS8fY/dOr2+ldWBj4j3wz8PnBtDbVI0vxxwOLWhhoZ5C9n5haAaOByGkkq3oKD4JGHJjX6ZqckaZJpz8gj4hrgiVPsWpOZX5jpQBGxGlgNsGLFihkXKEnFamm1Ytogz8xT6hgoM9cCawFGR0ezjj4lqTh+slOSNNmglx++NCLGgJOAL0bE+nrKkqT5qmNn5Jl5RWYuy8yFmXl4Zp5WV2GSVLyVv9PKMC6tSFJTzrqolWEMcklqyoL9f7XNNzslSZMZ5JJUOINckgpnkEtS4QxySSqcQS5JhTPIJalVXn4oSZrEIJekJj3zNY0PYZBLUpMOfELjQxjkktSo5n+5hEEuSY2a9Ht0vNeKJBUu6/8FaQa5JBXOIJekNrm0IkmazCCXpMIZ5JJUOINckgpnkEtSkw5e1vgQBrkkNemEVzU+hEEuSU36lcsNO3b5YURcEBG3RsSNEXFFRCyqqzBJ0swMekZ+NXB8Zj4N+A5w/uAlSZJmY6Agz8yvZubO6um3gOZX9SVJj1HnGvlrgS/vbmdErI6IDRGxYXx8vMZhJenxbWS6AyLiGuCJU+xak5lfqI5ZA+wE1u2un8xcC6wFGB0drf/2X5JUggbutTJtkGfmKXvaHxHnAi8BXpTZwP0ZJWk+aSAmpw3yPYmIVcDbgBdk5sP1lCRJmo1B18g/AhwEXB0RmyLiwhpqkqT5axhLK3uSmb9eVyGSpLnxk52SVDiDXJIKZ5BLUuEMcklq2pLfaLR7g1ySmvbarzTavUEuSU2L/qjt2G1sJUnDZ5BLUuEMckkqnEEuSYUzyCWpcX13PGzgXisGuSS1qYHb2BrkklQ4g1yS2uTSiiQVqOFfnmaQS1LhDHJJKpxBLkmNc2lFkuYR3+yUJE1ikEtS4QxySWravosa7d4gl6SmRcDiYxvrfqAgj4h3R8SNEbEpIr4aEU+qqzBJ0swMekZ+QWY+LTOfAVwFvKOGmiRpHmruEsSBgjwzH+x7egBNXywpSaVr4F4rI4N2EBHvAV4N7AB+dw/HrQZWA6xYsWLQYSVJlWnPyCPimojYPMXjTIDMXJOZy4F1wBt3109mrs3M0cwcXbJkSX1fgSQ9zk17Rp6Zp8ywr3XAl4B3DlSRJGlWBr1q5Zi+p2cCtw5WjiTNUw3eynbQNfK/jYhjgV3AHcAfDV6SJGk2BgryzDyrrkIkSXPjJzslqVXe/VCSNIlBLkmFM8glqXAGuSS1oqP3WpEkzVID91oxyCWpTQ18MMggl6TCGeSS1CaXViSpUPvsV2108H7kkqQZOPvTsOkSOOzo2rs2yCWpDYtWwMlvb6Rrl1YkqXAGuSQVziCXpMIZ5JJUOINckgpnkEtS4QxySSqcQS5JhYts4E5c0w4aMQ7cMce/vhjYXmM5Teh6jV2vD6yxDl2vD7pfY9fqe3JmLpncOJQgH0REbMjM0WHXsSddr7Hr9YE11qHr9UH3a+x6fRNcWpGkwhnkklS4EoN87bALmIGu19j1+sAa69D1+qD7NXa9PqDANXJJ0mOVeEYuSepjkEtS4YoK8ohYFRFbI2JbRJw3pBqWR8TXI+KWiLg5It5UtR8aEVdHxG3Vn4dU7RERH65qvjEiTmypzr0j4n8j4qrq+ZERcV1Vx2ciYkHVvrB6vq3av7Kl+hZFxGURcWtEbImIkzo4h39e/RtvjohLImLfYc9jRFwcEfdFxOa+tlnPW0ScUx1/W0Sc03B9F1T/zjdGxBURsahv3/lVfVsj4rS+9sZe61PV2LfvrRGREbG4et76HM5JZhbxAPYGvgscBSwAbgCOG0IdS4ETq+2DgO8AxwHvA86r2s8D3lttnw58md4v6nsOcF1Ldb4F+DRwVfX8s8DZ1faFwB9X238CXFhtnw18pqX6PgW8vtpeACzq0hwCRwC3A/v1zd+5w55H4PnAicDmvrZZzRtwKPC96s9Dqu1DGqzvVGCk2n5vX33HVa/jhcCR1et776Zf61PVWLUvB9bT+7Di4mHN4Zy+pmENPIfJPwlY3/f8fOD8DtT1BeDFwFZgadW2FNhabX8ceEXf8Y8e12BNy4CvAS8Erqr+E27vezE9OpfVf9yTqu2R6rhouL6Dq5CMSe1dmsMjgLuqF+pINY+ndWEegZWTgnJW8wa8Avh4X/tjjqu7vkn7Xgqsq7Yf8xqemMM2XutT1QhcBjwd+D6/DPKhzOFsHyUtrUy8sCaMVW1DU/34fAJwHXB4Zt5T7boXOLzaHkbdHwLeBuyqnh8GPJCZO6eo4dH6qv07quObdCQwDvxjtfxzUUQcQIfmMDPvBt4P3AncQ29eNtKteZww23kb5mvptfTOcNlDHa3XFxFnAndn5g2TdnWmxj0pKcg7JSIOBD4PvDkzH+zfl71v0UO5rjMiXgLcl5kbhzH+DI3Q+9H2Y5l5AvATeksCjxrmHAJU68xn0vum8yTgAGDVsOqZqWHP255ExBpgJ7Bu2LX0i4j9gb8E3jHsWuaqpCC/m94a1oRlVVvrImIfeiG+LjMvr5p/EBFLq/1Lgfuq9rbrfh5wRkR8H7iU3vLK3wGLImJkihoera/afzDwwwbrg97Zy1hmXlc9v4xesHdlDgFOAW7PzPHM/DlwOb257dI8TpjtvLU+nxFxLvAS4JXVN5su1Xc0vW/YN1Svm2XAtyPiiR2qcY9KCvLrgWOqqwYW0HtD6cq2i4iIAD4BbMnMD/TtuhKYeOf6HHpr5xPtr67e/X4OsKPvx+DaZeb5mbksM1fSm6N/y8xXAl8HXrab+ibqfll1fKNndJl5L3BXRBxbNb0IuIWOzGHlTuA5EbF/9W8+UWNn5rHPbOdtPXBqRBxS/eRxatXWiIhYRW+p74zMfHhS3WdXV/wcCRwD/A8tv9Yz86bMfEJmrqxeN2P0Lmi4l47M4bSGtTg/xzcoTqd3lch3gTVDquG36f3oeiOwqXqcTm899GvAbcA1wKHV8QF8tKr5JmC0xVpP5pdXrRxF70WyDfgcsLBq37d6vq3af1RLtT0D2FDN47/Qe+e/U3MI/DVwK7AZ+Gd6V1cMdR6BS+it2f+cXuC8bi7zRm+telv1eE3D9W2jt5488Xq5sO/4NVV9W4Hf62tv7LU+VY2T9n+fX77Z2foczuXhR/QlqXAlLa1IkqZgkEtS4QxySSqcQS5JhTPIJalwBrkkFc4gl6TC/T8QnrfvAHSkNgAAAABJRU5ErkJggg==\n"
},
"metadata": {
"needs_background": "light"
}
}
]
},
{
"cell_type": "markdown",
"source": [
"Seems like scipy does not play well with lowcut being higher than highcut, maybe the critical frequencies should be sorted before passing them to scipy.signal.butter"
],
"metadata": {
"id": "mLcNzEvydEtC"
}
},
{
"cell_type": "code",
"source": [
"clean_ecg_signal_nk = nk.signal_filter(ecg_signal, lowcut=52, highcut=48, sampling_rate=sampling_rate)"
],
"metadata": {
"id": "WmFU14BBbmD-"
},
"execution_count": 8,
"outputs": []
},
{
"cell_type": "code",
"source": [
"plt.plot(ecg_signal)\n",
"plt.plot(clean_ecg_signal_nk)"
],
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/",
"height": 294
},
"id": "sc0eyauhcUpt",
"outputId": "1fa2cfe0-29e0-489b-cfdf-9568344f125e"
},
"execution_count": 9,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"[<matplotlib.lines.Line2D at 0x7fb0fc6916d0>]"
]
},
"metadata": {},
"execution_count": 9
},
{
"output_type": "display_data",
"data": {
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
],
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAEDCAYAAAA2k7/eAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAATK0lEQVR4nO3de7CcdX3H8feXhCRcWggkQkiAhMugqVXEKKB4qSIGakm1OhOGUbAg44Xaap0OFEut05niZeqtVEyRio4FlIJGiqXIpXhFguWOkQAKSYMEKHEUBQPf/rHPOSyHPUnO2T37PHl+79fMznkuv31+3/NL9rPP+e2zu5GZSJLab7u6C5AkDYeBL0mFMPAlqRAGviQVwsCXpEIY+JJUiMYHfkScFxEPRsRtW9H2ExFxU3X7SUQ82rXvoxFxe0TcGRGfjoiY2solqVkaH/jAF4ClW9MwM9+XmQdn5sHAZ4BLACLiZcDLgRcAzwdeArxqSqqVpIZqfOBn5nXAI93bImL/iPjPiLgxIr4dEc/tcdfjgAtGDgPMAmYAM4HtgZ9PYdmS1DiND/xxrAD+LDNfDHwA+OfunRGxL7AIuBogM78PXAOsr25XZOadQ61Ykmo2ve4CJioidgZeBny1axp+5phmy4GLM/PJ6j4HAM8DFlT7r4yIV2Tmt4dQsiQ1wjYX+HT+Knm0mqcfz3LgPV3rbwR+kJm/BIiIbwKHAwa+pGJsc1M6mfkL4N6IeAtAdLxwZH81nz8b+H7X3e4DXhUR0yNiezov2DqlI6kojQ/8iLiATngfFBFrI+Ik4HjgpIi4GbgdWNZ1l+XAhfnMjwG9GLgbuBW4Gbg5M78xlF9Akhoi/HhkSSpD48/wJUmD0dgXbefMmZMLFy6suwxJ2qbceOOND2Xm3F77Ghv4CxcuZNWqVXWXIUnblIj42Xj7nNKRpEIY+JJUCANfkgph4EtSIQx8SSqEgS9JhTDwJakQ7Qz8e6+Dh9bUXYUkNUpj33jVl/P/qPPzQxvrrUOSGqSdZ/iSpGcx8CWpEAa+JBXCwJekQhj4klQIA1+SCmHgS1IhBhL4EXFeRDwYEbeNsz8i4tMRsSYibomIQwbRryRp6w3qDP8LwNLN7D8aOLC6nQJ8dkD9SpK20kACPzOvAx7ZTJNlwBez4wfArhExbxB9b9b6m6e8C0naVgxrDn8+cH/X+tpq2zNExCkRsSoiVm3YsKH/Xu+5tv9jSFJLNOpF28xckZlLMnPJ3Lk9v3RdkjRJwwr8dcDeXesLqm2SpCEZVuCvBN5WXa1zGLAxM9cPqW9JEgP6eOSIuAB4NTAnItYCfwtsD5CZ5wCXA8cAa4DHgLcPot+tqGw43UjSNmAggZ+Zx21hfwLvGURfkqTJadSLtpKkqdPuwA+ndCRpRLsDX5I0ysCXpEIY+JJUiHYH/n99sO4KJKkx2h34kqRRBr4kFcLAl6RCGPiSVAgDX5IKYeBLUiEMfEkqhIEvSYUw8CWpEAa+JBXCwJekQhj4klSI9gf+dR+vuwJJaoT2B/4N59ZdgSQ1QvsDX5IEGPiSVAwDX5IKYeBLUiEMfEkqhIEvSYUw8CWpEAa+JBXCwJekQrQ/8DPrrkCSGqH9gf/LB+quQJIaof2BD/CL9XVXIEm1KyPw88m6K5Ck2hUS+M7jS1IZgS9JMvAlqRQGviQVwsCXpEIMJPAjYmlErI6INRFxWo/9J0bEhoi4qbqdPIh+JUlbb3q/B4iIacDZwOuAtcANEbEyM+8Y0/SizDy13/4mx6t0JGkQZ/gvBdZk5j2Z+QRwIbBsAMeVJA3QIAJ/PnB/1/raattYfxIRt0TExRGxd68DRcQpEbEqIlZt2LBhAKVJkkYM60XbbwALM/MFwJXA+b0aZeaKzFySmUvmzp07uZ56vclq5XsndyxJapFBBP46oPuMfUG1bVRmPpyZj1er5wIvHkC/W++ea4banSQ10SAC/wbgwIhYFBEzgOXAyu4GETGva/VY4M4B9NubH6MgST31fZVOZm6KiFOBK4BpwHmZeXtEfBhYlZkrgfdGxLHAJuAR4MR++5UkTUzfgQ+QmZcDl4/ZdmbX8unA6YPoayuqGU43krSN8Z22klSI9gW+c/iS1FP7Al+S1FMLA98zfEnqpYWBL0nqpX2BP94c/qYnhluHJDVM+wJ/vCmd735yuGVIUsO0MPDH8etH665AkmrVvsAf97JMX8yVVLb2Bf54vD5fUuFaGPie4UtSLy0M/HF4hi+pcO0LfINdknpqX+CPyycCSWVrYeCPE+z3/Pdwy5Ckhmlh4I/jodV1VyBJtWpf4DuHL0k9tS/wJUk9tTDwPcOXpF5aGPiSpF7aF/jO4UtST+0LfElSTy0M/M2c4d/yleGVIUkN08LA34xL3lF3BZJUm/YFvnP4ktRT+wJfktRTCwPfM3xJ6qV9gb+lKZ2H7x5OHZLUMO0L/C35zCFw9d/XXYUkDV15gQ9w3cfg0fvqrkKShqrMwAf47BF1VyBJQ9W+wN/ayzIf3wiPPTK1tUhSg7Qv8CfiwuPrrkCShqaFgT+ByzLv+97UlSFJDdPCwJ+g736q7gokaSjaF/jTZkys/ZVnwqbHp6YWSWqQ9gX+rN+d+H0ufefg65CkhhlI4EfE0ohYHRFrIuK0HvtnRsRF1f7rI2LhIPodmNsvgY1r665CkqZU34EfEdOAs4GjgcXAcRGxeEyzk4D/y8wDgE8AH+m334H7xO/BtWfBj77om7IktdL0ARzjpcCazLwHICIuBJYBd3S1WQZ8qFq+GPiniIjMwX+W8ZNPJdMme+dr/+FZm36z4Ah+fcAf8uQOu5MzdiZjGk/tOAeeepKcPgu2mw4kGdOe+R6ACIj2zZhJmnozZsxizl77Dvy4gwj8+cD9XetrgUPHa5OZmyJiI7A78NAA+n+GRx97gt0HeLxZa7/DrLXfGeARJWnzVk8/iDkf/OHAjzuIwB+YiDgFOAVgn332mdQxdpo5Nb/S49vvwmM77MWvdpzPht1fwqZpOzDtqSfYNH0Hku2IfAoCkiBISAieIokpqUdSe83cZe6UHHcQ6bgO2LtrfUG1rVebtRExHdgFeHjsgTJzBbACYMmSJZOa7pm1/aQndJ72zu/Cns9/xqaZ1W02nV9QkrY1g5hkvgE4MCIWRcQMYDmwckyblcAJ1fKbgaunYv6+b0e8Hz608VlhL0lt0PcZfjUnfypwBTANOC8zb4+IDwOrMnMl8HngSxGxBniEzpNCs7z7enjOc+uuQpKmzEAmvDPzcuDyMdvO7Fr+DfCWQfQ1UIteCUvPgucs7lxVI0kt1qgXbYfubSsNeknFKPtCccNeUkHKDfz9X1N3BZI0VOUG/oyd6q5Akoaq3MB/zd/UXYEkDVW5gT/3oLorkKShKjfwJakwZQb+tJl1VyBJQ1dm4L/1krorkKShKzPw/QRLSQUqM/Cnz6q7AkkaujIDf/4hdVcgSUNXZuD7kQqSClRm4EtSgcoL/J33rLsCSapFeYF/8pV1VyBJtSgv8Hed3JejS9K2rrzAl6RCGfiSVAgDX5IKUVbgv+QddVcgSbUpK/DJuguQpNoUFviSVK6yAv/Qd9VdgSTVpqzAn3NA3RVIUm3KCnxJKpiBL0mFKCfwl55VdwWSVKtyAv/Ao+quQJJqVU7gS1Lhygn82YvqrkCSalVG4M9/MWxXxq8qSeOZXncBU+4VH4DD3l13FZJUu/af9i5eBjvtXncVklS79gf+du3/I0aStkYBgT+t7gokqREKCHzP8CUJigh8z/AlCfoM/IjYLSKujIi7qp+zx2n3ZETcVN1W9tPnhHmGL0lA/2f4pwFXZeaBwFXVei+/zsyDq9uxffY5MQa+JAH9B/4y4Pxq+Xzgj/s83uBtv0PdFUhSI/Qb+Htk5vpq+QFgj3HazYqIVRHxg4gY90khIk6p2q3asGFDn6WN9LzLYI4jSdu4Lc53RMS3gD177DqjeyUzMyLG+5bwfTNzXUTsB1wdEbdm5t1jG2XmCmAFwJIlS/r/xnHDXpJGbTHwM/PI8fZFxM8jYl5mro+IecCD4xxjXfXznoi4FngR8KzAH7yY+i4kaRvR75TOSuCEavkE4OtjG0TE7IiYWS3PAV4O3NFnv5KkCeo38M8CXhcRdwFHVutExJKIOLdq8zxgVUTcDFwDnJWZwwn88Axfkkb0dc1iZj4MvLbH9lXAydXy94Df76efyTPwJWlEu99p6xm+JI1qd+B7hi9Jo9od+J7hS9Kodge+Z/iSNKrdge8ZviSNanfge4YvSaNaHviSpBHtDvxo968nSRPR7kR0Dl+SRrU78CVJo1oe+J7hS9KIdge+UzqSNKrdge8ZviSNanfgm/eSNKrdgW/iS9Kodgf+m/6l7gokqTHaHfj7HFp3BZLUGO0OfEnSKANfkgph4EtSIQx8SSqEgS9JhTDwJakQBr4kFcLAl6RCGPiSVAgDX5IKYeBLUiEMfEkqhIEvSYUw8CWpEAa+JBXCwJekQhj4klQIA1+SCmHgS1IhDHxJKoSBL0mF6CvwI+ItEXF7RDwVEUs2025pRKyOiDURcVo/fUqSJqffM/zbgDcB143XICKmAWcDRwOLgeMiYnGf/UqSJmh6P3fOzDsBImJzzV4KrMnMe6q2FwLLgDv66VuSNDF9Bf5Wmg/c37W+Fji0V8OIOAU4BWCfffaZfI9v/Ro89vDk7y9JLbTFwI+IbwF79th1RmZ+fZDFZOYKYAXAkiVLctIH2v8PBlWSJLXGFgM/M4/ss491wN5d6wuqbZKkIRrGZZk3AAdGxKKImAEsB1YOoV9JUpd+L8t8Y0SsBQ4H/iMirqi27xURlwNk5ibgVOAK4E7gK5l5e39lS5Imqt+rdC4FLu2x/X+BY7rWLwcu76cvSVJ/fKetJBXCwJekQhj4klQIA1+SChGZk39/01SKiA3Az/o4xBzgoQGVMxWaXh80v8am1wfWOAhNrw+aVeO+mTm3147GBn6/ImJVZo77CZ51a3p90Pwam14fWOMgNL0+2DZqBKd0JKkYBr4kFaLNgb+i7gK2oOn1QfNrbHp9YI2D0PT6YNuosb1z+JKkZ2rzGb4kqYuBL0mFaF3gN+UL0yNi74i4JiLuqL7o/c+r7btFxJURcVf1c3a1PSLi01Xdt0TEIUOqc1pE/E9EXFatL4qI66s6Lqo+0pqImFmtr6n2LxxSfbtGxMUR8eOIuDMiDm/SGEbE+6p/39si4oKImFX3GEbEeRHxYETc1rVtwmMWESdU7e+KiBOGUOPHqn/nWyLi0ojYtWvf6VWNqyPi9V3bp+Tx3qu+rn1/GREZEXOq9VrGcFIyszU3YBpwN7AfMAO4GVhcUy3zgEOq5d8BfkLnS9w/CpxWbT8N+Ei1fAzwTSCAw4Drh1Tn+4F/Ay6r1r8CLK+WzwHeVS2/GzinWl4OXDSk+s4HTq6WZwC7NmUM6Xx9573ADl1jd2LdYwi8EjgEuK1r24TGDNgNuKf6Obtanj3FNR4FTK+WP9JV4+LqsTwTWFQ9xqdN5eO9V33V9r3pfNT7z4A5dY7hpH6vOjsf+C/T+Vz+K7rWTwdOr7uuqpavA68DVgPzqm3zgNXV8ueA47raj7abwpoWAFcBrwEuq/7DPtT1oBsdz+o/+eHV8vSqXUxxfbtUgRpjtjdiDHn6+5p3q8bkMuD1TRhDYOGYMJ3QmAHHAZ/r2v6MdlNR45h9bwS+XC0/43E8Mo5T/XjvVR9wMfBC4Kc8Hfi1jeFEb22b0un1henza6plVPWn+4uA64E9MnN9tesBYI9quY7aPwn8FfBUtb478Gh2vrRmbA2j9VX7N1btp9IiYAPwr9W007kRsRMNGcPMXAd8HLgPWE9nTG6kWWM4YqJjVvdj6U/pnDWzmVqGWmNELAPWZebNY3Y1or6t0bbAb5yI2Bn4d+AvMvMX3fuy87Rfy3WxEfEG4MHMvLGO/rfSdDp/Vn82M18E/IrOdMSomsdwNrCMzhPTXsBOwNI6apmIOsdsa0TEGcAm4Mt11zIiInYE/ho4s+5a+tG2wG/UF6ZHxPZ0wv7LmXlJtfnnETGv2j8PeLDaPuzaXw4cGxE/BS6kM63zKWDXiBj5JrTuGkbrq/bvAjw8hfVB54xobWZeX61fTOcJoCljeCRwb2ZuyMzfApfQGdcmjeGIiY5ZLY+liDgReANwfPXE1JQa96fzxH5z9ZhZAPwoIvZsSH1bpW2B35gvTI+IAD4P3JmZ/9i1ayUw8mr9CXTm9ke2v616xf8wYGPXn+ADl5mnZ+aCzFxIZ5yuzszjgWuAN49T30jdb67aT+lZYmY+ANwfEQdVm14L3EFDxpDOVM5hEbFj9e89Ul9jxrDLRMfsCuCoiJhd/SVzVLVtykTEUjpTjMdm5mNjal9eXeW0CDgQ+CFDfLxn5q2Z+ZzMXFg9ZtbSuSjjARo0hltU5wsIU3Gj84r5T+i8en9GjXUcQefP5luAm6rbMXTmbK8C7gK+BexWtQ/g7KruW4ElQ6z11Tx9lc5+dB5Ma4CvAjOr7bOq9TXV/v2GVNvBwKpqHL9G52qHxowh8HfAj4HbgC/RuZKk1jEELqDzmsJv6QTTSZMZMzrz6Guq29uHUOMaOnPeI4+Xc7ran1HVuBo4umv7lDzee9U3Zv9PefpF21rGcDI3P1pBkgrRtikdSdI4DHxJKoSBL0mFMPAlqRAGviQVwsCXpEIY+JJUiP8HeWd/DzVeKoAAAAAASUVORK5CYII=\n"
},
"metadata": {
"needs_background": "light"
}
}
]
},
{
"cell_type": "markdown",
"source": [
"Accordingly, the default neurokit signal_filter does not like lowcut being higher than highcut :("
],
"metadata": {
"id": "6t0TfSaWduVh"
}
},
{
"cell_type": "code",
"source": [
""
],
"metadata": {
"id": "ynnPzsNUciUY"
},
"execution_count": null,
"outputs": []
}
]
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment