{ "cells": [ { "cell_type": "markdown", "id": "b36fe48a", "metadata": {}, "source": [ "# Data Workshop 2" ] }, { "cell_type": "markdown", "id": "1c953193", "metadata": {}, "source": [ "**Instructor:** Jared Brzenski jabrzenski@ucsd.edu\n", "\n", "**TAs:** Tommy Stone tstone@ucsd.edu\n", "\n", "\n", "This can be run as MATLAB or Python, depending on the environment chosen.\n", "\n", "MATLAB will be run first. If you want to skip to Python, scroll down to the **Python** header!\n", "\n", "For MATLAB, run ```pip install jupyter-matlab-proxy``` in your environmnet and activate MATLAB in the upper right corner.\n" ] }, { "cell_type": "markdown", "id": "8be50903", "metadata": {}, "source": [ "# MATLAB \n", "\n", "# Raw Data\n", "Lets say we are given the task of analyzing the hsitorical data fro mthe water gauge station at the [Prado Dam in Los Angeles](https://waterdata.usgs.gov/monitoring-location/USGS-11074000/#dataTypeId=continuous-00065-0&period=P7D).\n", "\n", "We want to download the [data](data/PradoDam.txt), clean it, and do some spectral analysis on it to see if there is anything interesting.\n", "\n" ] }, { "cell_type": "markdown", "id": "38a8b2ee", "metadata": {}, "source": [ "## Reading in the data\n", "\n", "If we do a quick check of the file, we note there is a giant header, then some columns of data.\n", "```\n", "# ---------------------------------- WARNING ----------------------------------------\n", "# Some of the data that you have obtained from this U.S. Geological Survey database\n", "# may not have received Director's approval. Any such data values are qualified\n", "# as provisional and are subject to revision. Provisional data are released on the\n", "# condition that neither the USGS nor the United States Government may be held liable\n", "# for any damages resulting from its use.\n", "#\n", "# Additional info: https://help.waterdata.usgs.gov/policies/provisional-data-statement\n", "#\n", "# File-format description: https://help.waterdata.usgs.gov/faq/about-tab-delimited-output\n", "# Automated-retrieval info: https://help.waterdata.usgs.gov/faq/automated-retrievals\n", "#\n", "# Contact: gs-w_support_nwisweb@usgs.gov\n", "# retrieved: 2020-04-29 18:30:02 EDT (caww01)\n", "#\n", "# Data for the following 1 site(s) are contained in this file\n", "# USGS 11074000 SANTA ANA R BL PRADO DAM CA\n", "# -----------------------------------------------------------------------------------\n", "#\n", "# Data provided for site 11074000\n", "# TS parameter statistic Description\n", "# 8183 00060 00003 Discharge, cubic feet per second (Mean)\n", "#\n", "# Data-value qualification codes included in this output:\n", "# A Approved for publication -- Processing and review completed.\n", "# P Provisional data subject to revision.\n", "# e Value has been estimated.\n", "# \n", "agency_cd\tsite_no\tdatetime\t8183_00060_00003\t8183_00060_00003_cd\n", "5s\t15s\t20d\t14n\t10s\n", "USGS\t11074000\t1940-09-30\t51.0\tA\n", "USGS\t11074000\t1940-10-01\t47.0\tA\n", "USGS\t11074000\t1940-10-02\t47.0\tA\n", "USGS\t11074000\t1940-10-03\t47.0\tA\n", "```\n" ] }, { "cell_type": "markdown", "id": "2e27df20-b2e9-4c36-ae28-8e3828df4cb9", "metadata": {}, "source": [ "This tells us the pertinent information about the file, where it came from, and what format the data displayed is in.\n", "\n", "In MATLAB, we can load this data in by giving the filename of the data location, and using readtable." ] }, { "cell_type": "code", "execution_count": 103, "id": "6b76be2c-0024-4690-a1d2-6db0af119162", "metadata": { "scrolled": true }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "Warning: Column headers from the file were modified to make them valid MATLAB identifiers before creating variable names for the table. The original column headers are saved in the VariableDescriptions property.\n", "Set 'VariableNamingRule' to 'preserve' to use the original column headers as table variable names." ] } ], "source": [ "% MATLAB\n", "% read in a text file\n", "\n", "filename = 'data/PradoDam.txt';\n", "\n", "% Offer a helpful hint if we cant find the file\n", "[path, name, ext] = fileparts(filename);\n", "if ext ~= '.txt'\n", " fprint(\"Wrong file extension given.\\n\");\n", " return;\n", "end\n", "\n", "% We could read in the data raw with\n", "% ff = importdata(filename);\n", "\n", "% Or, read it in as a table with\n", "f = readtable(filename);" ] }, { "cell_type": "markdown", "id": "27d0099e-e6d6-48f6-9832-fe0b06c44de8", "metadata": {}, "source": [ "WE can examine the data by viewing the ```f``` variable, and note it is in columns already. we are interested in column 3 and 4, the date and measurements." ] }, { "cell_type": "code", "execution_count": 104, "id": "da287276-a57c-4e49-93d7-2ac3ab1f7eab", "metadata": {}, "outputs": [], "source": [ "% MATLAB\n", "date = f{:,3};\n", "flow = f{:,4};\n", "% Convert the date from a string into datenum object so MATLAB can do date specific work.\n", "date = datenum( f{:,3} );" ] }, { "cell_type": "code", "execution_count": 105, "id": "13c01d63-a0b2-4707-8d23-6057d24ea7a4", "metadata": {}, "outputs": [ { "data": { "text/html": [ "
months = 29066x1 double\n", " 9\n", " 10\n", " 10\n", " 10\n", " 10\n", " 10\n", " 10\n", " 10\n", " 10\n", " 10\n", "..." ], "text/plain": [ "months = 29066x1 double\n", " 9\n", " 10\n", " 10\n", " 10\n", " 10\n", " 10\n", " 10\n", " 10\n", " 10\n", " 10\n", "..." ] }, "execution_count": 105, "metadata": {}, "output_type": "execute_result" } ], "source": [ "%MATLAB\n", "% Lets say we need a monthly average, we can do that with indexing ?!?\n", "months = month(date)\n", "\n", "for ii= 1:12\n", " indexes = find(months == ii);\n", " Total(ii) = sum( flow(indexes) );\n", " Average(ii) = Total(ii) / length(indexes);\n", "end\n", "\n", "date = datetime( date , 'ConvertFrom', 'datenum' );\n", "\n", "save ('save_data.mat', 'date', 'flow' );\n", "\n" ] }, { "cell_type": "markdown", "id": "3f84cb62-98d6-4eb6-bdab-8e798d1e685e", "metadata": {}, "source": [ "And conversely, if we are missing a newer MATLAB, or want to do this the way MATLAB is doing it under the hood, this is the equivalent to the above code. Be carefuly running it, it will clobber the variables in the above cells!" ] }, { "cell_type": "markdown", "id": "614b18bf-83a8-4f93-8b7f-60c1a1fddb63", "metadata": {}, "source": [ "## Removing NaNs\n", "If you scrolled through the data, you noticed there were some missing values, dates, etc. WE need to remove those from our data set. How can we find them efficiently?" ] }, { "cell_type": "code", "execution_count": 106, "id": "0cff4d02-bbce-4abc-8e44-8dcc9e970c55", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", "inad =\n", "\n", " 0x1 empty double column vector\n", "\n" ] } ], "source": [ "% MATLAB \n", "% clean NaNs\n", "inan = find(isnan(flow));\n", "flow(inan) = []; % This effectively removes the entry\n", "date(inan) = []; % Do not forget the dates as well\n", "\n", "% Other equivalent ways of finding nans in dates\n", "% isnat == is not a time\n", "% if date is a date string\n", "\n", "% If we did everything correct, then this should be equal to zero\n", "inad = find(isnat(date))" ] }, { "cell_type": "markdown", "id": "6cdfde0f-cfe9-44bd-aca3-98099dae3f26", "metadata": {}, "source": [ "## Interpolation\n", "We look at the data, there may be time gaps, holes, etc. But we want a nice, steady time series.\n", "\n", "**We INTERPOLATE!**" ] }, { "cell_type": "code", "execution_count": 107, "id": "c27f15f8-75a5-464b-8a23-866a53ecf7dd", "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjAAAAGkCAIAAACgjIjwAAAAB3RJTUUH6QkREDoAIaXPaAAAIABJREFUeJzt3X1wG/Wdx/GfJKuQpkA9QRmTMQz2RAi5OMUDd6c6OE6TTMCXBx6clpxRPCElHricE4WQS+lhkpK2U2C4wkBJyxDIJAEaHjppJzkewigkkwrXvZj47ibYCKdGEEMm2Idx4gfJku6PLTrxW1mSZVv7k/x+/ZGRfvrt7ncftJ9d7cZrikajAgAAo5mNLgAAACEIJACAIggkAIASCCQAgBIIJACAEggkAIASCCQAgBIIJACAEggkAIASCCQAgBIIJACAEggkAIASCCQAgBIIJACAEggkAIASCCRgwoyMjOzbt+/48eMTMrbjx4/v27dvZGRkQsYGqI9AQl75l3/5F5vNZrPZSkpKknQ7ffr0iRMnJnzqQ0NDK1eufP755ydkus8///zKlSuHhoZiLdqsXXbZZUuXLn388ceTDz5J8whMHgIJeWX9+vWXX375//7v/+7atSvWGIlEpPOMn/zkJw8//HB8h2AwqB+bNKB+PKMNGPt0YGAgyXRHG8PIyEgkEtGP8PPPPy8sLKytrS0oKNi4ceMtt9ySZHIJpyXVA6glCuSXlStXXnDBBdFotL+/XwjR2Ng4bdo0IURlZeXg4GA0Gt22bVts+//d73736KOPWq1WIcQVV1zh8/liA+7YsaOoqKi2tlZ7e88990yfPl0Icdddd2kTeuyxx7Qxz5o1680334wNuG7dumg0Wltba7FYhBCXXnqp1+vVTzcajeonHY1G16xZI4SYPn36smXLhBD9/f2xWRNCrFy5Unu9ZcsWIcTRo0e1t9Lk9NPS1wOohkBCvpECqbS09M0332xsbBRCvPbaa9Fo9L333rNarfPnzz906NBbb70lhHjwwQfPnz9/ww03OJ3O2IDTpk3btm3b0aNHtbff/e5333zzzdWrVwshjh071tzcLISor69/8803b7jhhksuueT8+fPxgeTz+f77v/+7u7u7qKjo5ptvlqb72Wef/ed//qd+0l6vVwjhdrtff/317373u0kC6fXXXxdCPPvss9pbaXLStBLWA6imYCLPtgD1LFq0aPHixUKIJ5988vz580KIa6+91mw2FxUVLVq0aOfOnUKIo0ePfvDBB2fOnPH7/bEB77jjjq1btwohzp07J4SorKzUxrNr164PPvhA+0lt9erV3//+98+ePet2u48ePXrDDTfEBv/www/feuut8+fPRyIR/XSFEAcOHNBP+tSpU0KIhoaGqqqqU6dOrVu3LvnchcPhhJOTppWwHkA1BBLynPabmNmc7HLpvHnzZs+e/Y//+I/xjRdccMFo/c1msxZIsTwQQsRf9XnhhRfq6+sbGxvvuOOOQCAw1kkXFBTE/h3Nu+++K4QoLS1NZ3Lp1wMYiEBCXjl37twnn3wSiUROnTo1c+bM0bqZzWa/3//WW29ddtllQohPPvnkn//5nz/99NPPPvtstEHefvvtN954Y8+ePUKIq6666qKLLhJC7NmzJxKJPP3009OnT58/f34skz788EMhRG1t7YUXXnjy5Mm5c+dK03U6nddee61+0lrA/OY3v+nv73/mmWf0ZQQCgZdeeundd9998skn582bp50AJZxc/LRGqwdQi9G/GQITqb6+XtuwL7nkkvgrOocOHRJC7N69W+v2wAMPaN1efPHFHTt2FBYWam83b94c/fq9CbG3N998s3ZTQ6x9x44dWiwVFRUdPHgwfkC/3z9r1iwhhNPpXL169aJFi/TT1cYgTToajWqXqQoLC7XrXtI1JCGExWKx2+1btmw5f/681p5wcvHTGq0eQCmm6FdbOTClaDdwx34WCwaDBQUFCX/ZO3fu3EUXXbRhw4Z///d/j0Qi0i9pwWDwG9/4RsJJJPxImm7CSY+MjJjN5uQ/M6YzOf08jlYqoAJ+ssMUJeVKyj31aCGRZMCEH+mvDOm7Jb96lP7kxjqPgLH4j7FACgUFBfX19f/wD/9gdCFAnuMnOwCAEjhDAgAoIZ+vITkcDqNLAKCcD5buuOrAPUZXYYyOjg6jS0gmnwNJKLn0HQ6HglVJ1C9S/QpFLhSpfoViEoo0bfJO7AhzaDEqXio/2QEAlEAgAQCUQCABAJSQz7d9K/5rKQBDmDZ5o48tMLoKYyi+V+QMCQCgBAIJAKAEAgkAoAQCCQCgBAIJAKAEAgkAoAQCCQCgBAIJAKAEAgkAoATjAykcDg8PD8e3tLa2dnV1JWlJ2QEAkHMMfvxEV1fXhg0brrvuugcffFBrcbvdDocjEAhUV1e73W59S8oORs4PACBTBgfS008/vWrVqpMnT2pvfT5fcXFxU1OTEKKmpqaurq65uTm+paSkJHmHuro6s9n40z4AwFgZHEiPPPLIoUOHYm+7u7tdLpf22m63+/1+qSUQCCTv4Pf74x8UG3ut8t8TBIDJk0PPzlbrZKK9vd1qtWqvLRZLT0+P1NLW1pa8Q09PT/wIO76SrTkAALV0xDG6lhTUCiSn0zk0NKS9DgaDNptNaqmoqEjewWazZb9sAMD4KRdIXq9XCBEKhTo7O+12u9RSXl6evIPdbjd2FgAAmTH4GpKkrKysqKiooaGht7fX4/HoW1J2MHoOAAAZUvGJseFw2GQyxd8sJ7Wk7KBR/NmIAAzBE2ONrmJUap0haSwWS/KWlB0AADlHrWtIAIApi0ACACiBQAIAKIFAAgAogUACACiBQAIAKIFAAgAogUACACiBQAIAKIFAAgAogUACACiBQAIAKIFAAgAogUACACiBQAIAKEGtQIpEIu+9995f//rX+MbW1taurq4kLfoOAICco9AD+gYHB++8887KyspPP/20oKBg+/btQgi32+1wOAKBQHV1tdvt1rfoOwAAcpFCgdTS0nLllVeuX79eCPG9731v+/btPp+vuLi4qalJCFFTU1NXV9fc3BzfUlJSInWQnmIOAMgVCgVSZWXlM888c/jw4ffff7+2tlYI0d3d7XK5tE/tdrvf75daAoGA1MHhcMSPM/ZW5cfIA8DkkfaKKlMokKxW69y5c3fu3NnX1/fQQw8JIdrb2ysqKrRPLRZLT0+P1NLW1lZVVRXfQRonOQRgiovfDSoeTgoF0pEjR/r6+vbu3fvll1+uXr36t7/9rdPpHBoa0j4NBoM2m01qcblcUgdjSgcAjJtCV1w6OzvtdrsQ4uKLL54xY8YXX3zhdDq9Xq8QIhQKaZ9KLeXl5VIHY2cBAJAxhc6Qli1bdt9997W0tPT395eWlmrpUlRU1NDQ0Nvb6/F4hBBlZWXxLdJbo+cAAJA5UzQaNbqGrwmHwyaTKf5muZQt+g4ah8PBNSQAEtMmb/SxBUZXYQzF94oKnSFpLBbLWFv0HQAAOUeha0gAgKmMQAIAKIFAAgAogUACACiBQAIAKIFAAgAogUACACiBQAIAKIFAAgAogUACACiBQAIAKIFAAgAogUACACiBQAImi2mT1+gSgFxCIAEAlKBcIIXD4ebmZr/fH2tpbW3t6uqK7yO16DsAAHKOWg/oO3To0M6dO+fOnTtr1iztEeZut9vhcAQCgerqarfbrW/RdwAweaby41Yx2RQKpMHBwaeeeur3v/997AmwPp+vuLi4qalJCFFTU1NXV9fc3BzfUlJSInXQP8gcgOIIOWgU2n0fO3asqqrK6/Xu37//7NmzQoju7m6Xy6V9arfb/X6/1BIIBKQOhlQOABg/hQLp3LlzLS0tZrN52rRpa9asCYfD7e3tVqtV+9RisfT09EgtbW1tUgdpnI6vZG0uAEApjjhG15KCQj/ZmUym5cuXL1y4UAjx9ttv/+Uvf3E6nUNDQ9qnwWDQZrNJLS6XS+ogjbOjoyNb5QOAiuJ3g4pnkkJnSE6n0+fzaa+7u7uvuOIKp9Pp9XqFEKFQqLOz0263Sy3l5eVSBwPrB5BQ8v+Pxf/WQoxCZ0gOh8PpdK5du1Y79Zk1a9asWbOKiooaGhp6e3s9Ho8QoqysLL5Femv0HGAScd0byHsKBZIQorGxMRKJCCFiN8s1NTWFw2GTyTRai74DAANx6ICMKbcTN5vNUrRYLJbkLfoOyEtT4bedqTCPwGjYjwMAlEAgAQCUQCABk4if4GCUXNz2CCQAgBIIJACAEggkAIASCCQAxsvFCx5J5NnsZA2BBAC5Jy8zj0ACACiBQAKQVeM/tFdhDElGmJfnLtlBIAEAlEAgAQCUQCABAJRAIAEAlEAgAchz3GWQK1QMpDNnzpw9ezb2trW1taurK76D1KLvAADIOWo9MVYIMTg4uHr16h/96EcrVqwQQrjdbofDEQgEqqur3W63vkXfAQCQi5QLpF/84hcLFvzt+cc+n6+4uLipqUkIUVNTU1dX19zcHN9SUlIideDRsQCQo9TafR8+fPjSSy8tKyvT3nZ3d7tcLu213W73+/1SSyAQkDpkv2YAUxzXqCaKQoHU09Pz4osvNjY2xlra29utVqv22mKx9PT0SC1tbW1SB2mcjq9MfvlAfmJvm+sccYyuJQWFfrL71a9+dc011/zxj3/8r//6r8HBwdmzZzudzqGhIe3TYDBos9mkFpfLJXWQxtnR0ZG1+gFAQfG7QcUzSaEzpFtuuWX27NlWq9VsNhcUFFitVqfT6fV6hRChUKizs9Nut0st5eXlUgeD5wGAAjiry1EKnSFdf/31sdeDg4Pf+c53hBBFRUUNDQ29vb0ej0cIUVZWFt8ivTWsdACTxrTJG31sgdFVIBsUCqSYJUuWxF43NTWFw2GTyRS7fU5q0XcAYCwiBJlRMZAkFosleYu+AwAg53BWAeD/cfUFBiKQAHwNmQSjEEgAACUQSADAeaESCCQAgBIIJACAEggkAIASCCQAgBIIJACAEggkAIASCCRAFdx5nAUsZJURSAAAJRBIAJCMaZOX86rsIJCAzLGfkrBAMB4EEgBACQQSkMM4I1EWqyYDagVSJBL585//fObMmfjG1tbWrq6uJC36DkCumOK7rSk++5Ao9MRYv9+/bdu22bNnf/jhhwsXLlyzZo0Qwu12OxyOQCBQXV3tdrv1LfoOAIBcpFAg2e32PXv2mM3mSCQyf/78NWvW+Hy+4uLipqYmIURNTU1dXV1zc3N8S0lJidTBbFbrnA8AkCaFAkkIocXJ6dOnS0tLhRDd3d0ul0v7yG63+/1+qSUQCEgdHA5H/Ahjbzs6OrIzCwDynmmTN/rYAqOrSJe0V1SZcucTw8PDW7duvf/++4UQ7e3tVqtVa7dYLD09PVJLW1ub1EEaW8dXslU+MGH011e44oIMdMQxupYU1AqkUCi0cePGu+66S4t0p9M5NDSkfRQMBm02m9RSUVEhdchgonzJAUAFCgVSOBz2eDxut7uyslJrcTqdXq9XCBEKhTo7O+12u9RSXl4udRjrREkjTLb828YUnyPFy0MSCl1DeuKJJ06ePLlr165du3YJITweT1lZWVFRUUNDQ29vr8fjEUJILfoOAMZEwcshCpaE7FAokO699957771XamxqagqHwyaTKXb7nNSi7yBh40ZuYYvNDpazghT6yW40FotFChupRd9hPDjfRx5j84bKciCQlML3OW+wKgHV5HkgsdMBgFyR54E0fkQaAGQHgQQD5EfM58dc5BwWex4jkGA8djE5jdWHiUIgAcAkIrDTRyDlIb4A0Jg2edkYkEMIpK9R4durQg3Zxx8SjTemec/mgprKKwVZQCBNUexZJCyQ7FBwOStY0pQ1JQJpYje4LGy+kz0JvoGAyibpG6r+F39KBFLGtPWn2lpUrR6kL/11l9laVnDbULCkvCQt5xxd7ATSuEzZi8YZzHWWr4skHMOkrqypuSXkh3TWnbHrNz/yJiUCCRNgsvNJZWOdkdya8Un9kSD9nWzso4xPMXNrsU+4XDl0JpAmXvZXvLI3ZY2Jykepyb/PGVeVcsDkNx+mOd0M9ubjH+eYzlCjjy2IPQlinEWO9buQcHmms/ue8DpTJmgGd6Iq+2UfzdQKpDQPE7K/LY5nDOlUm80bMcY/OQPPt/RzMdYD7SQ74uSjGudWpy84najI2lVS/cymMztJuklrarTFPp6tcbSFk34Gx09dej2ekjL7NCfkSSC1trZ2dXVlPHg2V6TD4RjrIJNdXgZHXhmMc2JHol+Mae6vEw6S5i4m5UmS/t90jH8vk06d6RzEpDm/Cc9mJuQwLs3T0IQP1hvnkdmknqKNdh6WfNgJOZTMrZTKh0Byu90HDx78+c9/vnfv3nT6T/gaSvK11B+jfbB0R5IxjPZitIOvJMeJYzo8HNMymYxxjmnqpk1ebTGO/0hTf9yaTjyP/xxurGMY62lNyskl3FrS6Z8kjWItWp8J2Z8KIa46cE/CqvTjlwpILv53Qv0YEo4w4adjSuWERjuRStKScoRilBlUnEKPMM+Mz+crLi5uamoSQtTU1NTV1cWeHivt+qXtJvrYgiSr2fT1xxunc9QTP0KT7unI+g5JxpawZZyZoS8p+bDx/RMOO9o8plPMWMW+YCl3mvE7jjSHSr4lJC9JGjz2WjqcH+vZQ2wQ7UXyHb00fqme5MVcdeAe6WsiTWK04hPOtUiUBFKAJdlmtPFIW50+BbXXSba0+BEm76AvcrQZ1I9Bn8Fj2mwSvk1H8jGM9lUd61SMYopGo0bXMC6vvvpqQUHBLbfcIoRYv379unXrYj/mSLsnMcr3U98hySApv+EJj6Ti6XcB8UPpv3gJs1P/Jc+hbS5j45/NhDtrae8/1slJO834EWrrOrafTb4ShW6PPCYJa9BvV1JPh8PR0dGRweT0o528kYynyIk12kFYrMKUh276fI11S3lgIVJtriLpruyqA/d0dHSoszATyvlA+tnPflZRUbFkyRIhxMaNG3/wgx9UVlZqH2VwtQYA8pvKgZTzP9k5nc6hoSHtdTAYtNlssY9UXu4AAEnO39TgdDq9Xq8QIhQKdXZ22u12oysCAGQi58+QysrKioqKGhoaent7PR6P0eUAADKU89eQNOFw2GQyxe6vAwDknDwJJABArrNs27bN6BryTTgcDoVCBQV/+zk0EomcOHGioKBg+vTpWsvAwEAoFAqFQpFIROvW2to6PDz87W9/W50itT4tLS2Dg4MzZsxQsMjh4eFgMBj6itlsNpvNWS4y5WLUWoaGhgoLC7UW1RZjwpZsFhmJRFpaWqxW67e+9a1Yo74AqSXLizHNIqVFrdpi1PfJ/taYHIE0wbq6uu68885Tp05VV1cLIQYGBurr67/44ot9+/YNDw9fc801QojFixefPHnS6/V+/vnn1157rdvt7u3tPXDgwJdffjlnzhxFijx06NADDzxgsVgGBgbKysoULPK1117bvXu31+v1er2//OUvKysrN27cmM0iU1Y4ODhYX18fjUaPHDly5MiRBQsWKLgY9S3ZLNLv92/YsKGnp+ell17q6+urqKgQQugLkFqyvBjTLFJa1KotRn2f7G+NqUUxoTZv3vzKK6/89Kc/1d7+4Q9/2Llzp/Z66dKl4XB4ZGTk7rvvjvX/05/+tGXLFu31TTfdFA6HVShyYGBg+fLlIyMjKhcZ6/nZZ5/dfvvt2S8yZYXvvPNOrCSXy6XmYpRajh07luUitUmEw+Gqqqpooo1Nasl+hekUGf36os7+uk6nwvg+hmyNKXEXwAR75JFHLrnkktjbYDAYe1tUVNTe3t7b29vX1/fGG2/4fD4hRHd3t8vl0jrY7Xa/369CkceOHauqqvJ6vfv37z979qyaRcY+2r1796pVq7JfZMoKKysrP/7448OHDz/99NO1tbVqLkap5aOPPspykdq9SKdPny4tLRWJvhFSSyAQyP5iTFmk+Pqizv66TqfC+D6GbI0p5fxt34q77rrr1q9fb7FYPvroo1OnTlksloKCgttuuy0cDh89evTZZ58tLS3Vzq+FEBaLpaenR4Uiz50719LSUlFRMW3atDVr1tx2220zZ85UrUitfWBg4J133tm8ebP2NzsMLFJfodVqnTt37s6dO/v6+h566KGDBw8quK6lluLi4uuvvz7LRQ4PD2/duvX+++8XQrS3t0tLSWppa2urqqrKcoUpi5Q6p+xgVIWxPq+88orhW6MeZ0iTq6Sk5LnnnrvoootWrFhRWFg4e/bswsLCFStWLFmy5Mc//vHQ0FBVVdVof2nCwCJNJtPy5csXLlx44403Xn311ZFIRMEitfY9e/b88Ic/FEn/ZodRFR45cqSvr2/v3r0vvPDC9u3b58yZo+BilFpuvfXWLBcZCoU2btx41113aX/rS78epZaKiorsL8aURUr9s781plNhfB/Dvy8JEUiTzmazLVy48MyZMzNnzowd1wshQqFQb29vdXW1Cn9pQirS6XRqvygKIbq7u+fNm6dgkUKIcDi8f//+lStXCjX+ZodUYayMiy++eMaMGVdeeaXhFeqLlFoKCgqyWWQ4HPZ4PG63O/YnKPXrUWopLy/P8mJMp0hpkCxvjelUKPVR4fuix092k27z5s3BYLC/v//RRx8VQvh8vl//+tczZ8789NNP161bJ4RQ4S9NSEVqB1Br164NBoMul8tutytYpBDi5Zdfvummmy644AKhxt/skCpctmzZfffd19LS0t/fX1paOmfOHMMr1BcptcyYMSObRT7xxBMnT57ctWvXrl27hBAej0e/HqWW7K/odIqUZLnIdCrU91Fha5TwH2MnXcK/IhEKhaxWa/I+2ZSwgEgkIr66EDpan2xKpwBji0w4dalRwcWYTkuWpSzJ8ArTqcHwItWvUEIgAQCUoEowAgCmOAIJAKAEAgkAoAQCCQCgBAIJAKAEAgkAoAQCCQCgBAIJAKAEAgkAoAQCCQCgBAIJAKAEAgkAoAQCCQCgBAIJAKAEAgkAoAQCCQCgBAIJAKAEAgkAoAQCCQCgBAIJAKAEAgkAoAQCCQCgBAIJAKAEAgkAoAQCCQCgBAIJSObAgQNnzpwZTwfJyMjIvn37jh8/Pu7SgHxDIGGKGhoa+td//dfLL7/8m9/85jXXXPPUU0/p+7S3ty9btuzhhx8ebSQpOySc7sqVK59//nmp3Waz2Wy2yy67bOnSpY8//niSMZw+ffrEiRPpTxHIFQVGFwAYw+12v/baa42NjX/3d3+3e/fuxsbGkZERj8ejfRqJRCKRyNVXX93c3FxeXh4bamRkxGw2m81/O5LTd9AGLCj42jcrEokMDQ1985vfTFLP559/brfbFy9e/Mknn2zcuPGdd97Zv39/wsF/8pOfBIPBl156aUzjB3JAFJh6/H6/EOLWW2/V3obD4dLS0lmzZvX39wshduzYUVRUVFtbq71dt26d1m3NmjVCiOnTp9fX1wsh+vv7Yx20F42NjdOmTRNCVFZWDg4OakPV1tZaLBYhxKWXXur1eqVxxgghVq5cqb3esmWLEOLo0aP6wbdt2xb78v7ud7/Td8jC0gMmCT/ZYSr6n//5HyHErbfeqr01m81///d/393dHQwGhRD33nvv3XffvWHDhvhBjhw58txzz9XX17/66qttbW0JR3vw4MH9+/c3Njb6fL7/+I//0Bo3bdp04sSJ7u7ugoKCJ554Ip3y5s+fL4T44IMP9IPffPPNVqt1/vz5hw4d0rplMH5ATfxkh6koEolILdqvcFr7HXfcsXXrViHEuXPnYh0+/PBDIcSqVasWLVp06tSpdevW6Ue7aNGixYsXCyGefPLJ8+fPxwZ86623zp8/H4lEYo3pCIfD+sGvvfZas9lcVFS0aNGicY4fUA1nSJiKysrKhBBerzfW8t57711yySUXXnihEOKCCy4YbUCr1SqEmD59epJPY1eYhBAvvPBCfX19YWFhQ0PD5ZdfnmZ57777rhCitLQ05eCZjR9QE4GEqejqq6+eN2/enj17Hn744QMHDqxdu/b9999vbGxMMkhpaakQYvfu3UeOHPnNb36T5oS086ra2trCwsKTJ08m7xwIBF566aX169c/9NBD8+bNW7RoUcLBzWaz3+9/6623Pv744zGNH1AcgYQp6tVXX12+fPm//du/LVu2bM+ePZs2bdq+fXuS/t///vfr6uqee+65W2+9VTuRSscdd9wxa9as+fPn33nnnbfffnvyzj6fb9WqVW+88caWLVtef/310QbftGnT8ePHb7zxxmPHjo1p/IDiTNFo1OgaAMNEIpGRkZFvfOMbafbXbvteu3btc889FwqFpNu7RxMMBtOfRDqDj4yMCCFiUx/n+AFFcFMDpjSz2ZzmrvzcuXMLFy68+uqrzWbzvn37br755jTTSAgxzrTQDy5NmjRCfuAnOyAtF1544T/90z9pNyw8/vjjL7/8stEVAfmGn+wAAErgDAkAoAQCCQCghHy+qcHhcBhdAgDlfLB0x1UH7jG6CmN0dHQYXUIy+RxIQsml73A4FKxKon6R6lcocqFI9SsUk1CkaZN3YkeYQ4tR8VL5yQ4AoAQCCQCghHy+7Vvxk1MAhjBt8kYfW2B0FcZQfK/IGRIAQAkEEgBACQQSAEAJBBIAQAkEEgBACQQSAEAJBBIAQAkEEgBACQQSAEAJxgdSOBweHh6Ob2ltbe3q6krSkrIDACDnGPzXvru6ujZs2HDdddc9+OCDWovb7XY4HIFAoLq62u1261tSdjByfgAAmTI4kJ5++ulVq1adPHlSe+vz+YqLi5uamoQQNTU1dXV1zc3N8S0lJSXJO9TV1ZnNxp/2AQDGyuBAeuSRRw4dOhR7293d7XK5tNd2u93v90stgUAgeQe/3x//XL7Ya5X/niAATJ4celSpWicT7e3tVqtVe22xWHp6eqSWtra25B16enriR9jxlWzNAQCopSOO0bWkoFYgOZ3OoaEh7XUwGLTZbFJLRUVF8g42my37ZQMAxk+5QPJ6vUKIUCjU2dlpt9ullvLy8uQd7Ha7sbMAAMiMwdeQJGVlZUVFRQ0NDb29vR6PR9+SsoPRcwAAyJCKT4wNh8Mmkyn+ZjkDgH6UAAASlUlEQVSpJWUHjeLPRgRgCJ4Ya3QVo1LrDEljsViSt6TsAADIOWpdQwIATFkEEgBACQQSAEAJBBIAQAkEEgBACQQSAEAJBBIAQAkEEgBACQQSAEAJBBIAQAkEEgBACQQSAEAJBBIAQAkEEgBACQQSAEAJagVSJBJ57733/vrXv8Y3tra2dnV1JWnRdwAA5ByFHtA3ODh45513VlZWfvrppwUFBdu3bxdCuN1uh8MRCASqq6vdbre+Rd8BAJCLFAqklpaWK6+8cv369UKI733ve9u3b/f5fMXFxU1NTUKImpqaurq65ubm+JaSkhKpg/QUcwBArlAokCorK5955pnDhw+///77tbW1Qoju7m6Xy6V9arfb/X6/1BIIBKQODocjfpyxtyo/Rh4AJo+0V1SZQoFktVrnzp27c+fOvr6+hx56SAjR3t5eUVGhfWqxWHp6eqSWtra2qqqq+A7SOMkhAFNc/G5Q8XBSKJCOHDnS19e3d+/eL7/8cvXq1b/97W+dTufQ0JD2aTAYtNlsUovL5ZI6GFM6AGDcFLri0tnZabfbhRAXX3zxjBkzvvjiC6fT6fV6hRChUEj7VGopLy+XOhg7CwCAjCl0hrRs2bL77ruvpaWlv7+/tLRUS5eioqKGhobe3l6PxyOEKCsri2+R3ho9BwCAzJmi0ajRNXxNOBw2mUzxN8ulbNF30DgcDq4hAZCYNnmjjy0wugpjKL5XVOgMSWOxWMbaou8AAMg5Cl1DAgBMZQQSAEAJBBIAQAkEEgBACQQSAEAJBBIAQAkEEgBACQQSAEAJBBIAQAkEEgBACQQSAEAJBBIAQAkEEgBACQQSMFlMm7xGlwDkEgIJAKAE5QIpHA43Nzf7/f5YS2tra1dXV3wfqUXfAQCQc9R6QN+hQ4d27tw5d+7cWbNmaY8wd7vdDocjEAhUV1e73W59i74DgMkzGY9bncqPcEU8hQJpcHDwqaee+v3vfx97AqzP5ysuLm5qahJC1NTU1NXVNTc3x7eUlJRIHfQPMgcA5ASFAunYsWNVVVVer/f8+fNz58612Wzd3d0ul0v71G63+/1+qSUQCEgdHA5H/Dhjb1V+jDwATB5pr6gyhc4nzp0719LSYjabp02btmbNmnA43N7ebrVatU8tFktPT4/U0tbWJnWQxtnxlazNBQAopSOO0bWkoNAZkslkWr58+cKFC4UQb7/99l/+8hen0zk0NKR9GgwGbTab1OJyuaQOhlQOIGPcHI8Yhc6QnE6nz+fTXnd3d19xxRVOp9Pr9QohQqFQZ2en3W6XWsrLy6UOBtYPICEiB2lS6AzJ4XA4nc61a9dqpz6zZs2aNWtWUVFRQ0NDb2+vx+MRQpSVlcW3SG+NngNMIm7EyhWsKWRMoUASQjQ2NkYiESFE7Ga5pqamcDhsMplGa9F3AHIXe3NMZWoFkoiLopjYXeCjteg7AAByDmcVyBlcigDyG4EEAHkoFw/gCCRgEuXiTgEwCoEEAFACgQQAUAKBBMB4/LYJQSABwIQjXzNDIAFA7snLzCOQAGSVCnvSCa9BhZnKAwQSgBzD3j9fEUgAACUQSAAwkTiByxiBBABQAoEEIM9xypIrCCQAgBJUDKQzZ86cPXs29ra1tbWrqyu+g9Si7wAAyDnKPaBvcHBw9erVP/rRj1asWCGEcLvdDocjEAhUV1e73W59i74DACAXKRdIv/jFLxYs+NsjnH0+X3FxcVNTkxCipqamrq6uubk5vqWkpETqwIPMASBHqRVIhw8fvvTSS2fPnj04OCiE6O7udrlc2kd2u93v90stgUBA6uBwOOJHGHvb0dGRpXkA8otpkzf62AKjq0DmpL2iyhQ6n+jp6XnxxRcbGxtjLe3t7VarVXttsVh6enqklra2NqmDNM6Or0x++QCgoo44RteSgkKB9Ktf/eqaa6754x//ePz48ePHj584ccLpdA4NDWmfBoNBm80mtVRUVEgdjCkdgEq4zztHKRRIt9xyy+zZs61Wq9lsLigosFqtTqfT6/UKIUKhUGdnp91ul1rKy8ulDgbPA4CJRrpMHQpdQ7r++utjrwcHB7/zne8IIYqKihoaGnp7ez0ejxCirKwsvkV6a1jpAKYwInOiKBRIMUuWLIm9bmpqCofDJpMpdvuc1KLvAMBY3AeBzKgYSBKLxZK8Rd8BAJBzOKsA8P9Mm7z8AAWjEEgAACUQSADAjQlKIJAAAEogkAAASiCQAABKIJAAAEogkAAASiCQAABKIJAAVeTBncfqz4L6FU5lBBIAQAkEEgCkwHlVdhBIAJAMaZQ1BBKQOXZVwAQikIAcRiIqi1WTAbUCKRKJ/PnPfz5z5kx8Y2tra1dXV5IWfQcgV0zx3dYUn31IFHpAn9/v37Zt2+zZsz/88MOFCxeuWbNGCOF2ux0ORyAQqK6udrvd+hZ9BwBALlIokOx2+549e8xmcyQSmT9//po1a3w+X3FxcVNTkxCipqamrq6uubk5vqWkpETqwIPMASBHKRRIQggtTk6fPl1aWiqE6O7udrlc2kd2u93v90stgUBA6uBwOOJHGHvb0dGRnVkAkPdMm7zRxxYYXUW6pL2iypQ7nxgeHt66dev9998vhGhvb7darVq7xWLp6emRWtra2qQO0tg6vpKt8oEJo7++whUXZKAjjtG1pKBWIIVCoY0bN951111apDudzqGhIe2jYDBos9mkloqKCqlDBhPlSw4AKlAokMLhsMfjcbvdlZWVWovT6fR6vUKIUCjU2dlpt9ullvLycqnDWCdKGgFjpfi3RvHykIRC15CeeOKJkydP7tq1a9euXUIIj8dTVlZWVFTU0NDQ29vr8XiEEFKLvgOgGsWvNyhYnoIlITsUCqR777333nvvlRqbmprC4bDJZIrdPie16DtI2LiRW9his4PlrCCFfrIbjcVikcJGatF3GA/O9wHAEDkQSEohrvLG1FyVU3OukSvyPJD4+gFArsjzQBo/Ig0AsoNAggHyI+bzYy5yDos9jxFIAMaFhMBEIZBgPPZoyGNs3ukjkPIQXwBoTJu8bAzIIQTS16jw7VWhhuzjD4nGG9O8Z3NBTeWVgiwgkKYo9iwSFsiUxapXx5QIpJzb4Ca74JxbIMCUMknfUPW/+FMikDKWcP0ZvlINLwAZS3/dZbaWFdw2FCwpL0nLOUcXO4GUWpJVO2UvGmcw11m+LpL9gwm2hMnonx3pVGVs5fmRNykRSJgAk51PKsuPPfJotGqz8wtS8iO/yRv5VJArh84E0sTL/opX9qasMVH5KDX59znjqlIOmPzmwzSnG+s2gUsvg8qTiD62IPYkiOQDTux0pdUav6AmdkLpDJ4yQTO4E1XZL/toplYgpXmYkP1tcVJl+eDIkGOxiZqifpc01gPtJD8VJh/VOLc6fcGj9dfvgrOwvvS1pTM7SbqlM6rkYxhrzfpJp2yPn7r0ejwlZfZpTsiTQGptbe3q6sp48GyuSIfDMdZBJuQwOePxZzaS8VSVzvd2rIsx4Y5pTPvu5PXE77zGupcf/15morItzfEkPJuZkMO4NE9D9Q/WS3PkSbqleYqWzlSSNybZ6vTLfzK+3SrLh0Byu90HDx78+c9/vnfv3nT6T/gaSvK11B+jfbB0R5IxjPZitIMv6fhLqiT9DXpMy2QyxjlW2mLMOJiTHLemE8/jD+zxDJ5O4KXctSXcWpJUG+uQJI1iLVqfCdmfCiGuOnBPwqqSzEI6j4KN/50wfiQJ38bmSP/p+I8ORzuRStKScoRilBlUnEKPMM+Mz+crLi5uamoSQtTU1NTV1cWeHivt+qXtJvrYgiSr2fT1xxunc9QTP0KT7unI+g5JxpawZZyZoS8p+bDx/RMOO9o8plPMWMW+YPr8SFhA/J5CWtH6oZJvCclLkgaPvZYO58e6cmODSLOccEcvjV+qJ3kxVx24p6OjQ+oQP4mExY821yJREkgBlmSb0cYjbXX6FNReJ9nS4keYvIO+yNFmUD8GfQaPabNJ+DYdyccw2ld1rFMxiikajRpdw7i8+uqrBQUFt9xyixBi/fr169ati/2YI+2exCjfT32HJIOk/IYnPJKKd9WBe/QnSbGh9F+8hNmp/5Ln0DaXsfHPZsKdtbT3H+vkpJ1m/Ai1dR3bzyZfiUK3Rx6ThDXotyupp8Ph6OjoyGBy+tFO3kjGU+TEGu0gLFZhytDV52usW8oDC5FqcxVJd2XawYc6CzOhnA+kn/3sZxUVFUuWLBFCbNy48Qc/+EFlZaX2UQZXawAgv6kcSDn/k53T6RwaGtJeB4NBm80W+0jl5Q4AkOT8TQ1Op9Pr9QohQqFQZ2en3W43uiIAQCZy/gyprKysqKiooaGht7fX4/EYXQ4AIEM5fw1JEw6HTSZT7P46AEDOyZNAAgDkOsu2bduMriHfhMPhUChUUPC3n0MjkciJEycKCgqmT5+utQwMDIRCoVAoFIlEtG6tra3Dw8Pf/va31SlS69PS0jI4ODhjxgwFixweHg4Gg6GvmM1ms9mc5SJTLkatZWhoqLCwUGtRbTEmbMlmkZFIpKWlxWq1futb34o16guQWrK8GNMsUlrUqi1GfZ/sb43JEUgTrKur68477zx16lR1dbUQYmBgoL6+/osvvti3b9/w8PA111wjhFi8ePHJkye9Xu/nn39+7bXXut3u3t7eAwcOfPnll3PmzFGkyEOHDj3wwAMWi2VgYKCsrEzBIl977bXdu3d7vV6v1/vLX/6ysrJy48aN2SwyZYWDg4P19fXRaPTIkSNHjhxZsGCBgotR35LNIv1+/4YNG3p6el566aW+vr6KigohhL4AqSXLizHNIqVFrdpi1PfJ/taYWhQTavPmza+88spPf/pT7e0f/vCHnTt3aq+XLl0aDodHRkbuvvvuWP8//elPW7Zs0V7fdNNN4XBYhSIHBgaWL18+MjKicpGxnp999tntt9+e/SJTVvjOO+/ESnK5XGouRqnl2LFjWS5Sm0Q4HK6qqoom2tikluxXmE6R0a8v6uyv63QqjO9jyNaYEncBTLBHHnnkkksuib0NBoOxt0VFRe3t7b29vX19fW+88YbP5xNCdHd3u1wurYPdbvf7/SoUeezYsaqqKq/Xu3///rNnz6pZZOyj3bt3r1q1KvtFpqywsrLy448/Pnz48NNPP11bW6vmYpRaPvrooywXqd2LdPr06dLSUpHoGyG1BAKB7C/GlEWKry/q7K/rdCqM72PI1phSzt/2rbjrrrtu/fr1Fovlo48+OnXqlMViKSgouO2228Lh8NGjR5999tnS0lLt/FoIYbFYenp6VCjy3LlzLS0tFRUV06ZNW7NmzW233TZz5kzVitTaBwYG3nnnnc2bN2t/s8PAIvUVWq3WuXPn7ty5s6+v76GHHjp48KCC61pqKS4uvv7667Nc5PDw8NatW++//34hRHt7u7SUpJa2traqqqosV5iySKlzyg5GVRjr88orrxi+NepxhjS5SkpKnnvuuYsuumjFihWFhYWzZ88uLCxcsWLFkiVLfvzjHw8NDVVVVY32lyYMLNJkMi1fvnzhwoU33njj1VdfHYlEFCxSa9+zZ88Pf/hDkfRvdhhV4ZEjR/r6+vbu3fvCCy9s3759zpw5Ci5GqeXWW2/NcpGhUGjjxo133XWX9re+9OtRaqmoqMj+YkxZpNQ/+1tjOhXG9zH8+5IQgTTpbDbbwoULz5w5M3PmzNhxvRAiFAr19vZWV1er8JcmpCKdTqf2i6IQoru7e968eQoWKYQIh8P79+9fuXKlUONvdkgVxsq4+OKLZ8yYceWVVxpeob5IqaWgoCCbRYbDYY/H43a7Y3+CUr8epZby8vIsL8Z0ipQGyfLWmE6FUh8Vvi96/GQ36TZv3hwMBvv7+x999FEhhM/n+/Wvfz1z5sxPP/103bp1QggV/tKEVKR2ALV27dpgMOhyuex2u4JFCiFefvnlm2666YILLhBq/M0OqcJly5bdd999LS0t/f39paWlc+bMMbxCfZFSy4wZM7JZ5BNPPHHy5Mldu3bt2rVLCOHxePTrUWrJ/opOp0hJlotMp0J9HxW2Rgn/MXbSJfwrEqFQyGq1Ju+TTQkLiEQi4qsLoaP1yaZ0CjC2yIRTlxoVXIzptGRZypIMrzCdGgwvUv0KJQQSAEAJqgQjAGCKI5AAAEogkAAASiCQAABKIJAAAEogkAAASiCQAABKIJAAAEogkAAASiCQAABKIJAAAEogkAAASiCQAABKIJAAAEogkAAASiCQAABKIJAAAEogkAAASiCQAABKIJAAAEogkAAASiCQAABKIJAAAEogkAAASiCQAABKIJAAAEogkAAASiCQAABKIJAAAEogkAAASvg/6UiavHEIO90AAAAASUVORK5CYII=" }, "execution_count": 107, "metadata": {}, "output_type": "execute_result" } ], "source": [ "% Interpolate for hourly times\n", "% using interp1\n", "\n", "firstday = datenum( date(1) ); % first time\n", "lastday = datenum( date(end) ); % last time\n", "xq = firstday : 1/24 : lastday; % make a vector with perfect hour spacing\n", "yq = interp1(datenum(date), flow, xq, 'pchip'); % Interpolate to that perfectly spaced vector\n", "close all\n", "subplot(2,1,1)\n", "xq_datetime = datetime( xq , 'ConvertFrom', 'datenum' );\n", "plot(xq_datetime,yq); title('Interpolated Data');\n", "subplot(2,1,2)\n", "plot(date, flow); title('Original Data');\n", "\n" ] }, { "cell_type": "markdown", "id": "02d9ddf9-250f-4de7-aa47-fa31c39fef0a", "metadata": {}, "source": [ "## Filtering Data\n", "\n", "OK, now lets filter this data, to see some sort of pattern in the data.\n", "\n", "This is usually the first part of spectral analysis, where we look and see if there are any time (frequency) components to our data." ] }, { "cell_type": "markdown", "id": "c15559c4-10bb-4623-8d3b-46d9bda988ff", "metadata": {}, "source": [ "Generic Example:" ] }, { "cell_type": "code", "execution_count": 108, "id": "78324a5f-e0d9-47a9-bb1f-83fa8827b63c", "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjAAAAGkCAIAAACgjIjwAAAAB3RJTUUH6QkREDoDuKye0gAAIABJREFUeJztnV+oX9WVx/dPQ+lYTaGhFGQefDEXAgMdJ0KvJQFTgim0tB0IJRSH6QxMtIVK3mZEBa0vlTCgfeiDUChISx+mdISBiCODI2IJQxPy4BChVfBBZmSgvSlJb6z+5uEkx333n3XW/r/2/n0/FHvzu+d3zt77rL2+e621z7mr9XqtAAAAgNbc0roBAAAAgFIQJAAAAEKAIAEAABABBAkAAIAIIEgAAABEAEECAAAgAggSAAAAEUCQAAAAiACCBAAAQAQQJAAAACKAIAEAABABBAkAAIAIIEgAAABEAEECAAAgAggSAAAAEUCQAAAAiACCBAAAQAQQJAAAACKAIAEAABABBAkAAIAIIEgAAABEAEECAAAgAggSAAAAEUCQAAAAiACCBAAAQAQQJAAAACKAIAEAABBBx4J09erV73//+8eOHfvyl7986dKl1s0BAACQxGq9XrduQySPPPLIF77whVOnTrVuCAAAgAz0GiG9//77v/3tb6FGAAAwDL1GSC+//PIvfvGLT37ykxcvXjx06NBTTz114MAB45itra0mbQMAALFcvny5dRO87GvdgEiuX7/+zjvv/PznP9+/f/+Pf/zj55577sknn7QPEzX0W1tbaA8B2rOItCahPTTS2qPEL9N7Tdl94hOf+PznP79//36l1KFDh957773WLQIAAJBEr4J0zz33XLx48erVq0qpt99++3Of+1zrFgEAAEii1xqSUurcuXM/+tGP7rzzzvfee+/555//7Gc/axwgMF4GAICGCPeKvdaQlFInTpw4ceJE61YAAADIQ68pOwAAAIMBQQIAACACCBIAAAARQJAAAACIAIIEAABABBAkAAAAIoAgAQAAEAEECQAAgAggSAAAAEQAQQIAACACCBIAAAARQJAAAACIAIIEAABABBAkAAAAIoAgAQAAEAEECQAAgAggSAAAAEQAQQIAACACCBIAAAARQJAAAFxWq9YtAEMDQQIAcFmvP9YkiBPIzr7WDQAA9MSsSet166aA4UCEBAAAQAQQJABAAHN4hJQdyA4ECQDAZbX6OFOHlB3IDgQJAMAFIgSKAkECAAAgAggSAAAAEUCQAAAAiACCBAAAQAQQJAAAACKAIAEAABABBAkAAIAIIEgAAABEAEECAAAgAggSAAAAEUCQAAAAiACCBAAAQAQQJAAAACKAIAEAABABBAkAAIAIIEgU+JuYAABQDQgShf53miFOAABQFAjSApMm6X+5GWwCWH8AUB8IEgAONjY43qjOAmlAkBaY5qfunsCGsJnB8cYqMROMSVEgSBS6M9oorwQ2mc1UYiYQ7KJAkCgwITcZBMfACQS7HN0L0rlz544dO9a6FWA0Njk4hhKDVvQtSO+///7PfvazK1eutG4IGI1NE6GZTVZiJhDscvQtSI8++uhjjz324Ycftm4IAIMAEaKBYBdlX+sGxPPCCy9sb2/ffffdxDFbW1vTD5cvX67SKAA2go2toPTY69kNyqdXQfrNb37z6quvPv/88/RhfB2aw3AQysb6JoLhx2RKWE19nH4Yvsv9ortB4eLUa8ru2WefvXr16pkzZ86cObO7u3vmzJnExN08tab/IjvMBxthbTZhTIzNZpvQZVCaXiOkhx566He/+9308yuvvHLy5MlV8iSYJ5hyhUpYABLMzqijISp9Q3sck0SG6TImeyt6FaRDhw7NP+/bt+++++4rfUU7RwG6Bjc0EX2z2WCjB9toRa8pO51f//rXeU/o29CJB+IGw3dDNzzjxOm+c7PZSPuhMdmbMIIg5ULPgMMKgxjJE6lM5ZB+x4RTT6Vz2pg+IA4I0g0406lfF1OU8TzRvG1MLfXIZwljjIlt8ITl99tNg1mJFSZ7dSBIN1icTmO4mBL0Oxo+p8OPkHxH9jsmTjYqf4UCUkMgSFw2zS6HXxgaK4zo/g7vrIe3BJthtgt2BwSJy6ZNy+EfKzF8jdHfTY6GffEi8lcGGI3sQJC4DO+gbYZf+xvo/d20e62j325DhDbEEhITtiAaCFIASCsDG0QPg+GLj/EoSAUgSA7gWQw2bUD4KbsNT+4Nie8+Ys1RAQiSA18kvmlrYb2uuwmpiYj+Di9CqO3r6G8Xm9g0n1CaDRKkRYsxcuVGJL5pa2Gjv9GpiV4mqtG1jbrXBJvT/QhD3TSfUIENEiS6Aml8Yh+waQaXq7+9FH6dagQ2B46hGvEi/QQ9iGCDBEmRFUi6egkLS0lNoPALuoA2VL4B97IIE8hmCRKNzwp9v+oLzsTwhYlITYTiHO1+fVO/LY+D+fYN+kVKWIRFAEHaFDirNt/DocRXOMgp/OqvKTM+9P0zDudod7Fw9m1upg8YDOcrPJxpfM4N3YQRy8VmCRLtGYc3Kc6qzbmbQyd0KERFV/NryvS9Uvazn/pvoyXE+dSa/IUzsW4Q3vJ0bP+gW4tzZHzPJuqn6mIhIoQNEiR6LtELw4EnoQ+fZocOhbShs72DLVH28WM7Yp3N6anB4sppPoCzX9c41aZZUTQbJEjEC1GGKRTRcFJn+lDQe42GwWkMvn86P+EjJ3tp4IsXU04lrY80zrTk4mg4b+jwnqQcGyRIBvOaZVq2OC2vrxlFw0mdhR4zDPSuyxlnSseHL/9DnL8hnGWK/TNxqiHtxEDyDe2UDRKk9HVf0LSUBmfC6HWj7jrog4h45kCQ01l7lwfhc52uSrjPIvKW+jFOxrAWPTAyPlGuv0nBv6H66mSMsSrEBgmSsVadXYZtZ76vz4xRqIyeG4nbHCrjjFEWiwF8v0xXHzuC6PJid8aYDnQ37RVJxJmNZEynY1WODRIktdcNOfc+Mb9ubNbq0fuovaPBF6dpRhmfCJ9gvjtF7zZcPEZ1mMB03iljcOwe2YpuG4xzs6J9BrHMDSbmgt1BTr+cCeGuXUc5NkiQUqqOdkzQuyXZmxdSeuTb/NoLhAvWjxGuuxyIXjB7x1y3Kdeaj3mGVkRbr9gedccGCZJzv6/+T/q7ylpUdo1vE4edRrePGQbDAAwHah9j667A/XKL2Llrguiqie+69Vct2W/QqIsVCWyQINkErevtw3KtKFvB31hlfKs79JJy0FeU9dC+XT5cLEeJhS8zdn/t34r1yMRKlPiK82cdW6r5chthkHFIuxeLjCxIb711Wf8nkSUPMiObxd0QouYqkbFJTFnIDBec99258Hcmde3hclZKetRpndANY/o/55Jq6KmqYaiFb0oGzVPDZoLWOtXGqrtK3siCpPwru9Ab4wzS5x84C65WNUxf8ZlOOwSdVr4zYn6oQ99x5ldkYq8emMV5AufkctaohKxa7CnpM2PCEvRjfOkT56Vrom/WEDg9DYYVpNVKHTy4pSzXqTxPn/DXRMbxwm+zsZWOOIwz8ZQ2bnpw4HQ98gla3jo764yWmCdsgp1jjLZen9d2ZsKFr1qUv1XMlJ1uHtJ8ghHCSmZf6waUYr1Wq9Vlw2n6fKiQVVshOCGRs/uLAkaPZ1CJriGhoaFvP8h6vWdAZGJUR1ZWVYwPc002j0xDDCcQ8V2l2YlzrHyzTAKIkERw8OAWXRWYYd6qaGsTlazwSUgodK1Y5lIxCz5/1BdGzpnZ/qBpIsTsncFZXNucyTr9PPzEXTVWrsqZWEYWJB1nPBRaB/J97vuuvTAR5bbSi2q+eR5xwozE9SUOvePyZ7uP9NUJ5/ytxodT8knEToQKCZSNnjZvzyIjC9K8y85ethjWw4EwYlqriMisDrReTgRZqk9+nMW5+nBqOdHrDxoJeuyEaI8zYRB6HqctiXV/dt5ygo59fXlvVV7Ro5HTEiYjC9JMlmXC4reIzHLKdRPxzRO75hx0TkNrfT6oVZagWsLQd9OlOQLC7Wa5R9L6GwRtLSnVNRDKsIK0urnLTsVqCfOwxYJh29St7hz1OCbdg/h2NMy/NX7Q29AQZi6leTt1cjVmZW10znJyIyNnWMVqxR3zavhqn07NTmm2KCvqgsF32U0/+6IB2qWGXGhhl5EQMkqjPm5zwKSf3Heh5hvwuijtGmQZNHqTmPPIOHwVfjkYE9/2A5ysrw7H2gGHYSMkpfZESPRUTLE5Yxlo/7Z5DakQRMaPnoEV8mmL64yi96LEybMPGnGq0PbTx88xuij719WdjpD4Z/P9SmxlUSAjC9KMUfOoluLPUruSic/FSChoG8UtZuiWDr/vTRwT/6IjGaoOkXIvivyQUQ4bIUjKky4PssVFY7IjBiH2VyhKsJfYzP6mp0lpFuOzOuTtXeKgBYVWcfGBEGsnsIOV0H2zoDSbIkjpcyZUveTE6SVmndPBMTOfJaIo/iBXuB2+/G3KCRMHjd4/Frew8GH3Xcj6TC920gW5UHMi1nzNp39fbIQg2evKElYyG59R8Od8sQ6+PUXRZwsqHfGPicC3AtBdIVF5LtQk4qIRp2KSHhBHB2G+z5tL0USWclHoFePWphsrYyMLkv7nJ4x8Wt4loY6++HKe3OnH64dTifUt37qSWDAm9nHxW3bZ3+5j0SqCc9FT/+b6rmjHWBn3RzjPJiQwaohhkPyhEJViqcnIgjTvslN7g4PSAZMuOYZhrVbmokm5MgmFyLg3nT7eeaHEPkZM0cqukEiNVri5/CsWihLsOVWuTBhHoXQxTbS01DcbCYwsSDNGsFK0ls5ch87WltHgFnM1FfKWReH42cWt3qX3ms9slB9xoq+02jZjokm02pbu5viwguT0DuVuj54mUn63SGeNEleUvhR5ufyk71Tl+uiDv/itX7Frko9VpBFuFPM4+LLo5YhO2amo+WIfaSdjhDOsIK3Xe2pI1VavtvE5K0a+Ukdiw+yTE/nJ9EEg7NsuJhn7miLmBjFF+fsp6jgj+yo1nSAzU1enMaquBEYkk4sSvakhV+VJ12P5a5FhBUnHV27Ni+2CfVfXj1+nPaXL6UiTSTj/YGRs5ukRFF2pfLJdzR07NanC1Tn7OEpf1/5nNXxOWaeyJMdJSxz6klR+PGQzsiDNmxrq54uJvHnextibJtSSG6rgKfhLQjrJYAdYHSGkYucMmkvgrNQ2uWvOJITxw2DY/dJTEXNuQH73hxWk1WpPyq7C5Tg3WzeO+ZNE9CjEmanjSGOF3X1GBLnyb09wbl5IGShf8FoOe23ehTvIQkTRoujIFC2jLsK/Vpb6sb0MMiJm4QwrSDKxvVKJLQa6N5eQOM7iblK2SDUfAdXuRthBc2ld9GXI5/8664sZm2d02fbIlVcGdtecDUgcBMO3rPt80XivgnTlypWnn376yJEj29vbZ8+etQ9Yr6m/h9TwPmW/NJ2mY16uWm2DU+Hw9cgnsUIiD3orQROcmbrSxu/rshHv2lELETQHXZ2zsagadqBP9DHXIPRLr4K0Xq/vvffe11577fXXXz9//vxrr71mHzOn7JwJ1jqUvpCvNqDn8Zr7a2J2OVfToa6z5vI/lGoyMKMPhdMj1xyiXEnXuCtKo+hCylY+yUPho1dB2r9//wMPPKCUuuWWW+66666dnZ3WLXJT2iZ8Jk6ETa1ctq2OoSEd0Sl6816dLkvQwpX15+bqSzVHhOybSAf6GZtEtKoEztWSz1bpQSCa7ZQfCQYZSq+CNLOzs3PhwoWjR48an6+0P2Guv0NobEQlJw0S2xYROTF3muSC1sISLfEJs53brJyk4uBbzhdtalsfTXdtcRCITAB95q2trem/0w+S6VuQdnd3v/Od7zzxxBN33HGH8Sv9DtXcbmeQZQLEnWSwbV2Lk9n4YfEr2aFHu0RjfB5KTljMpPKdalhDMn6l/8BpVVCRaT7m8uXL03+nHyTTsSDt7u6ePn36wQcfPHLkSOu2eOHYzaKziK6RiNKkcl6SWFrKCRALYXuoDdxozqetPRjCw7lT0cW/Tu9+r4L0wQcfPPzww9/61remShKNBCsk4OgNZ2Xk9PhyPDKdsguaPES+Qsg+DptqdSznpoYNp+2el/midilL/5Uzwzx/blcHF+nx7u9r3YBIXnzxxfPnz7/55puPP/64Uuro0aPPPPOM7+CG7omfT0/f7unbJ90FRvEj6GB7e6FvejfEWcZPv0G2hxK+NE5vXsS4GRdtu69nJtRKDU0aktV60J5NmxoaVo+c6D5UNyzfGmpG9zv0Vmn7VM3dE+GDnNuffN1U2tDRA9i8yzZ2wwhZdWIfY4+GUntGRuA4TBA3l8PiuBkHNN/LYN935TcJ41f6LSamxoSvhjp/a2trS3IlqdeUHQcJauR0uMYPi2tbokZin1ksibvsjC8SBWGZQ0EnZOKGwlhzGA5O5jhMxDXSTnBxVl1tU5dG+K6Wmur8iv0z57rGIHdhGL2m7BZpvixy4lu90pE4f/uNgZAqArOu41vXB/XCmRWRiS/iiYBYdEvDDnD5vSYO9qV8mydviVRhaKtCzcMe5NVKKXVZgk/wMXKEJOHxI47NZXGazjMIdMdEhOeLAmmJomeXtO7P2NFM9tWDTL9jr9YjnCyxqrNPKHMcCIj8dhdRTgrDRkjTWqB1K0yc9lTUtqRpktFxpq74AkpleTRp/dXxtY3oDo3z4JF81uJocMZK7Dg4K3+K7FRQjteQ7em/Bw9urVZyg6QBIySx9qdKLtZ8jkn4aER4z5XrsRtfLVcU9irEjpOi6wRxZ6gPJ/x1/la3Z99GBvuEdgq3Ocy6VyKEJUgYBB8DCpKS7YWdpJuIs8sCx8Hpjxbd6KLecCa5KIxtVFnO4/yVNOymEhI1Vz7mH3xSRKxR5G94MSCmgzPoMX47oc+I6WeZ9mAwpiAlslartWpguSkxTUc2p8PprDGj9JK4jeSUnY5TivQCSXQvhKfs+Dt0FqXL+bkzCTYGtugSfbT1WLhhTIz5HFL0iBs6tFKZB4eeJPaiRv95cd20eP62LKqIsVHK9jLOM3QXG+nQY2JAVIyM/VSLJ49pqlqpMjNixjB75donppRpKmrvGsWwHztlxx8Tu8t5B4GoIdk/O53D/An9Rb60NwcR0sfYUVGrOGleIYaWBCQ7ZTqtRNeTciW4pJHYI6dvyj4+esKgzozwiZC3hexiJHOC2F2uOQiGVcSlHO2hc55QGgMKkmSnTEOsaBaRbGQ+jJ5yXI9aur8lxqF+CpdpA87oMG9mRu949vBIuVprJKYiupM4AqUXpnExvW95quuNr6RkJIEle8gBBSkvNT2RbW18cZJsZD74y3l9zZgiV6HUDw7ilrH0Ro9czLFCznNaOxEiFmTO40sEizMpg+Br7Z7z8wYhoo8Sns4kGLCG5EsfL+IzshILQ3cDrDlJ7BfSyVs2cCbKi66UjYycM2lu/Nf39ZytsuyhQk3RqI7subpVVyvU8Y8b4xqB+SaoTANC94K2h48blrtatkhRt+CshNlD8XFjAo1BrNcfMEIqPUtLE5qmcEpX5KVdAYGdsCoRKxiJmj2XIyfPXG+rQPbcnd3sgJJ7iw0sRpCUZTR8yaj1eqEyr59B/3plb1sio6t7ADutzWoVWXITy4CCNCF50H3MqV47h05UC3L11Dmpss80X2t9q91Fn1vIKRMdrzMmHJUtp0nVymZOHZpYef7w4OJJFjMKuciY0bUXYfPsjusFrd9iGVaQesRWHXvlaPggI5bKOwNLOGWOwPj8i1Gb3dOeWh3Pi88dh97NjN0PlaLEsbIDoDh3ySnMcOB3p3SsnJGOVufDClKuAtL82wpOys5d+IrV85FN8jbZMZISs/AYfSQ2NVQuITR5HsAm+76p+v2it1b6zNu3nawtGUePeWeZ4+CsucpkWEEK26S0pEacwzISbTFCZmYK5YrVi9Tf3m3DT1cW3UUmh8UNhMQ4dDQ+vropp4AaxMGDW5K9xLCCJHnQCYwSvbOQoGc5Ou2mgd4XfY4ZOzXKpjXavSzKiVOVjdjReXB90vdAcyIh50Swy66+r3NouBwJaqovo8A0A8k6PawgFRp0e79ZCSMmGk97IsmmRkCnGpirxcoUcl7GzV1p7xVVUSvigEsn9Cj6u/TNtUdjwt54lnF/qQR88kxscfJ9SwmbODTDClJpZguen8nIAlG3v3Fdcl1cOcGV+YSWZ+l9B38uVq6nTPLGx83dMacv9qZT+4GkCJr33WaxR/oCpSO9WWRAQaoQp+tfKW3NrQoqQhi7+7Qr0Z2vrdP2kW0p9HTa4udj7OsxoLeSGplbzs7VXhjwTQ1Kce/TRPpEmp8cTzzPjK/xukU6o4eUyRkxDiUeVjdqSPatLOSAGnafuN3OG519BHJpScSA+Po4/1btjQaIA6LJvIE766TwjUBix8V6/TEjpAwzdn3zf6xjMztIYjlspNF9vw3Cm0APHIS8zwbONN7j6x+EbH58aTedb+dhPbcSYgkRLN5Te1Eyf7GSPRQeAfNqURfCpgaJJKmRbnYr6xM/JZZFjqsk17dZ7wEKHIT059WJvJzv0dG8RA9CogzHpex8H0a2gWi/bxD8pwodjTpbKJMIdAjz25XKrVONXan2Af0yoCDxMY1mNrvVTePTf664PiJ2uxr5Db4JcmVjcRCs03Jb4IfoS+k55nAc65hBiBsHZ5jrPObGhdbLxwfhbTYxCIrlkcOawa7hV8UYBGWZgX+JZvywfCkyEcL5ivMTJ5J1a0BBIradzI7D7YaU8tqP3xmVxlfBzvjkr36cUuQgqOURyLswXFvv8iq7tcxwQDb+QYjuePRu3VIOmh6EAnOBk7IrBCXJqt4g+BIDdt1I+V0cZ2rgz0/UxpnRWnAWtCNmHFnuXfRZKtjO7psb1pmDwDgsdDToTRwZN26Y5w9alDCOrPOXSsrm65IHQd34OHMSu9DzWN7qqcowCOrGL5OGwtnrFBuQ7PIHjJBssqmR8i6LBD7KsEiMGqnMC8Mbp/Rf1/5VezVSzcyg+INZOQahe6QOglFDGpIBBWkd9MrRIONT4Qcns1jkz0DoICxNwozbz+yUXUEizKBk2xaztZlHJnou1B0EVW5RoooMQuJ0oDtLPylvM73LTqyeDZiym2Dl6yKMz//FCs8hxZ+Qs4vBg/dVFBUzNiWIDBCtsxhfzG4Gvudv8qfs1rGLrZLTQR+EvPMijxrp3yUnUdRJ915hqfsBD19KnZ0DRkj2cxvk0XHXUEq1SVa0XxTPNFodZzitb3dlBCVXmlWXsSlD7V6u5Gx96J7ShCslfLeWQ2A+sdcjAwqSDrUfN9F6SrrjoJpK8sU4h/hmQN6mWKfPff7g3ZXs8+qXSDuXdqqiYq8PBW8QVmpNrfRLPhKQd7+72y0Idgg3rpBvOoiVrgEF6YbVrshHQwLdkHsemp6tD0+kXWb5kLnjQZ4oF3nHwXuDEhNulifK+wYH4vGsyOf5ne44MeHWzwaH1G11BJS/yRwvJk4NpOzqsVrxbn8uC5kfDvDVWsIJehpukbgZaC+HKVWWat8z1CDchPbC04Bwlibz5Uo8rr/nsuHn5kQGzj4uL01KJu5qpOxyO4SipAyFWDVSQwrSMhnvB2frRBrlHgsl9iOELZCLuQk5NQN9QLzjszb+tTJ+kIh/XaIrUIA9FCr7ZU3ZWWfPd6qm67PlJ6nVaq1WYvN1ahMFKTw2X56NBYKkj69ewnr8HYrfDuRxxymUXcp5F/pRV+0kUnTgT9b5pGiYmqJ2apGnCr0ymTPsYoW0eYKkFizGNyEXJuG0git/p7P5aHJRnOVs6dR8ApToO3dYapUQ4vBu6IiCSmAWUOVc77KLzlj6fuU9WNLSRIL5MRlQkCirdef/90jOdPPC0hRlHkNxXyrQtPi7iTK0XNIkXMZ/1zjZqi480QxnQ0ceim3xUCUWKOHrElYmU77/l5q2G1CQFoaasYwN9kE3v696WYxYWRqmGlGHeQr7EoiIDIILaTe/1g2ZhNMxSh0NgoppLdc2CiSxh2dAQfLCKJwsmpo3U9GLsZVev8ubhO425GhXX0GSA8fe75imLybuJJiBg6K3qeSew4HZJEFSYbucY1jP/7+a/5d6zuyktSgoSFICJ2FWN9Tr0qS0ZAoLl72T0bOhI8gVRDyiJ9QzCGBAQXLXkHLPQMITGabW0PIcqaoci+KF7Ln8+KDCDcEgKBHhstf1Bz/4Gy5RnvVZF7vdWjGgIHlrSHXuvlhPVLNhAjyRmzqDIKa7bkrsghMZKS4YHm+3LSd3Qh0s1iGIZEBBckRI1WzCY+JSPHLWVngnoYy+eslUOFlGrCqr3Jv+CWQOAncbf0AaP660LGVAxDCgILkjpJr3XcaayPFC69xQM1agJyo2CM5PhVLTOMUOgqq4LlFSHEIXDChIJqQ1pFihWE/Eqd/e/KyUk7aprEn8QSjaiOY4xqGmJSgRg7AH15OwBbsvwCF0xICCRLuhyDd0hTeiFcwnYct233PdOprkrhuX7KvkpYnxb5tNdMf1WyVNlaUyoCDtofAMFOuJTPzvK6tzxWpQmleySd4H1No9jrP4AFaNRcnNpmg/Np0eTYRBoEOQSseC9NJLLx0/fvyLX/ziD3/4Q+q4wtZQ841e6VRyQP5LF/VHEX+MseyALG38LUf+P0C3RDfv0WnVHpEOQRq9CtLOzs4PfvCDf/mXf3n99df/+7//+1e/+pXjIFeyuFB7zDPLmYTW5oJCgyD9cX1XEbthfNBmNBreAQmW0FAS5DgE2fQqSG+88cbhw4f379+vlPra1772y1/+0n3cav7/4t5HVLpmDxXf/eq7bks8adtCV5OwNKkfHk2ICpKYr4yqOiMQJC3RqyD94Q9/uP3226efb7vttp2dnbbtcdDi7SkZ/75ANtyBU91Xp9QvoRl4BqHAdfwDW2W8me/RabA+az4XJFRVmw/CEr0KEovqFkAHSW2o6vZZDwZWTVtVTNt6L1HrFgSV0OrFiDc/bY+ENtR1CCKe/wuhV0H61Kc+de3atennP/7xj3O01Bz6nQVVS/rNtXCmbUvaFGuEvUJGiF9qNQhC5kLbuyBkEEh6FaR77733woULu7tE40trAAAgAElEQVS7SqlXXnlle3u7dYs+pr4n8iZqJMSIZBvKruBahEduaoSCwTsMa9PWHQuRZFXpjoh4MDycfa0bEMmBAwe+/e1vf/3rX//MZz5z4MCBb3zjG61b5KH826+zvMw4Iyu1dhfVm0wGUTPQGoRymx730GJdwllqzH+duXyL3BTNW3rXiOU3vrba1ZLOap3+d+ql4XqZXWWjdxQStf1+Ba+lX7TdhqIgQcrYKmLYs18ruD2q7E2h1iUtLKFtexYtodB1vW3Y+wu7PXlbwpqAIj1/ryk74VR7+7W08Ggi6Enhglm71uERZ5NL8bJzo0FotbWh/YscrfPX3OTi7n5rh8AHglSL9fz/K+2zYruf2z5vYVO5iCJzBhYbhKB1STVL4L/EJNcskOaLQ1+Hn+gQWr03KyMQpFIsbrcrtftZhi+u/DolzqK4iSQ3fhLAygvVHwT+n8sq5Y4F+uICj0OE7WoRma9TEKSi+LbbjZqlWYD0RDkvJHSuKaUWPFFOJA+CooKkiNHI9XfKyxG6PstsEjIdggsIUi3q7H5unaVZpk5DhM9AchCyJXLFDkLm5Qd5uvJbaSKpc3e6Co8UBKk2FSxBQJbGS+5J6PDd8iWZl7NK0qS6r2aIgXwcIo8eezZ2tpoRzYIkO3ctdqXS73NIvbDncYSszyRJq98GkOOZJGNvyJ7fNd3qzcX1TFLO86/0H0WOgCr5dJpISTYfTir9kKJvO4VUc1CIkBqQY/dzF/Xbcht/+fXb5j7IjWuTi30U89lSIkaUFR8blLZVSXOBIvn+NH+jbkYGF6To2ZhxEVHpmSSRS0Iv5drV0QxMG4RFNyT37uuUeDqth37fIM1cKRuQn7t2MawgJf590rJp1uyT0PMsevN3/VJBkvUcBtFaVp2/l/BoItkT+X7RE+XMs9puxkAyPg6x3J323Q1mWEGSYHwzeYMkiX95LBSvOzV/EVbn76JwopO3jdq6RFr3mX8nKQMlI4Mi1RfX+iwVchAkb2oYVpCkQ+6uCVNTT2QgRJJTPFFAF7LO50pV3xKeKOqet6xyJ2XmAl5yn+UdpllcOf+N+PpSLHguSN5t6we77EzW6yIrCMd2O9cWI90EAwzIs19LtAl6thjRjXcPiydjGUchA3DD3me1fE8TbnX7JTO53U73xXP3Zb4tiQn/jfjO3aTL3fnY07gmkazB2AMipEYwXABrTeSfgSUmYZwpGy3Z88/w5zD4i+Jo8s7Y5RuxFC5Te9z3XqlL/JGiHRlQ3fesS3K+Tl7GZpzlQaAvJdhOIEgm5e5WxGvNnDUV8+sVY/PowZla5X6tWWJjC2wvLD1j6wxCaXJ658A9CK12ORd35ewhDfrrHh0BQapKoidivlc/O0IWhg46mYH6XyII/aO6Nu4Hj6oPQn7vHJi95Hw3898ZKvQ0yM2PlEpeW4ifCzQQJBPC5vL75ZTEXcVFcXbXk/OvvBeYgSUE2JCibG8B70SS1aI2pLjjWoNQPEJKH4Q9JxNcLPIwuCDlfTC2lDlGtLHYDGxQ8IyehLWaWmNMUjxRIzUqtT4LPW0/kswiriMhg4BNDQ3I+2Bs8ZwV4/wf96XkDHSKbnELjvBEFd3QalUqZjL+rVTUIDQiYn3G3R7G71cBMxDhr0vOBWxqaMDkvqMr/OWMMiJ3vGevEWl8hYL0QvvgjX8rxZ6HdRfF9TyUmEFo/2QSpwGDxUYzYsygPgMK0g1Pt15Pm/2zPB8qJXfsN75eHnyjSB6EPK2Ypb/JcMoYhMZwBoExAnEzwjfZnfaQYiRU83I4BO+XBPuJAQXpBqtVwHNkjm87Psx1Ixfe8Oa8yrqSD6pprFSw2G4Q5hGYbWC1qpvloJ1RgUGwbzqnv2VNZXEQFDUC1Z7DK2gY9FxQkWoknAEFab3+2HDypuwyGp93+69uhWvtZ/231qmyNWupj9l9kGMcVtZUXC8PgnFO+oDoXjTI3UUNQqhJ2De9asraf6hS5CCkXyKEBqUXwiEMGiKv1pLjtyimEnS09UzfNc7g/DALrIeuk12w9wo1X5NDsvBwCa+RuTyRffdVGX8U9jIOftVa9gI5JoseMgiJ3fdNCufnibYR+dbKtOWIZJc/ZoSUsi3KaVvTh7V99ypPQEB9l700blZQYYyAypqicaZlbH3yHRwEq83sQQg4Zw7ip9jNFgbctRBLiGyWpi6+BMl6Xf21UnuPDrKEHhlQkLKEMtW0R9pi1tfx6XMRmRxJCIkvJ4qqsuNyCX3Xm5qr2bnOM3sPTt5elAGobufRzICCRKsRPdMkB7M26cbH2TikrwqLTj/Hy+7YXww6nr7Lzl120yDYX0w3GN0vR58h0RcTtfpFC5EwZbJI0ew39M0szgOMD30/MDGEOSnh0bkaqSEFiXaai78lovJCcy++ApT8J16YAlPoydBosi+udQz1pYsEuf5AztSLCmLsPonfzy6u7XzJzOA2CHCmU/t9HoBI5SUywjMbmRhQkFJo5XYjFkflfLHTPTFnXekB1H13fEgRuGTJXjlIoagYJ9L4Ea409PtOl5EMjFmTRZ9ST+Hn4MGtcidPZ0BBSlqmeWYUcy9NxgSO77fpcb2O0WAiDnAOSKcw07a+w0rnLTnHFPJZQdLLL7RwTxj3HGumupEeFDrLSLQTiE7ZOZG21KjGgIKUkl+SUGEqaouLXUislEir8TpZbGTbXjQMfYJe4Ug46JQNeL6+27/KOFD6Rsp5m66dGBAVKw/JgIKUuMVO0LZvow1ZnZTPm/jeSmDIvE/1iy4P67jpxQXNkC5pLp9MzAbPX/gb5pGC80brb6csWjtUaeFOF2sysQwoSHqs7VzsE0Y2z8D5JL7zEJdOp+hKUP+BKSpzEkMfEEK5hZNRUYYUJ7V3+vj24BHbGbIPS4W1CH+DD59RzaMcAwqS4u0REk6hlaCx9GM+GOt0TP1OtsXXI9Hbfzm77wpR1BgWLr20my7vVkxjN3z6hlKDuCUmp4+5HpMI7e8YZacBBWnxKRMi6+0sZua6tED4NVuDfpWeZm295oPYWKU6vOMTEc32ZeRms8mYsrtxEk2KKnhbX2JA/6ccSjxZLIEBBUknos5BmODiNBZlsk6IzalBb0bxpexyOegmE8xWI99hi59wfkVTM0kVZ7cVxLicqzU01WnPzhtdM2W32PeRpGhicEFipqR8GPYn7fnQdEI9Ua6nICOoM/Hs9T7tg5jZvyxk3OtvpAFCU0y6B59P5cs9yIeTFJm7HJTeFP5MkkAGF6R0jD0OdBghX658dSDfItH+esbJRtMqI6Hv4JgICoA4exG5LSkzAkSrmArq07NO4a8b+PmDoONpNkqTBvzzEzOhwbWveOALC/paD9KzLihRaTyxwb+KBIhS2WIVwfhEly7DSCQPgh4FGhs0jMHxRU7G5wPsITLuqY7dR/3+EoOZs3l7X/GdKFGSXf7gEVJcct+XpnP6317wzRBadAmcbrdfl5QLmSbhXLMHBbu2m56zBbaADYNzl119Ix9y/4KTwQXJWUNaXMIY+XHlWkJ253l9+ztCdxUa2Uvju727pDzbw0KeJ60DrT12YoBu9mLKTkivOdDJamfKhJ8mMYYoMX+bRYqEv8tuwJSdM9W2iDNlp9jZCWfqRjiLaQrjMKWtEH0ZDLHo0kukGYkEnf2JM2XnS38JwXnjnCNA32Wn2/CNsGQ4d9w+0s7Q+gbTWMBJGBnJLn/ACImZQwi6K3R2Yl7+SLC2IJhVJfrD7npNqG8WdDWSPPntdnJa6wspBMaFFXAOhdF9O+MCfAwoSMqfc9OhnYWuQPRM66Wi6xsBO5kT2hGB3ieoSbMv5i8pmBkwgSNDxL72AQR0mk74XAiFsA3muAm0BJmMKUg+mGZhKBAzO6EfKdD+7DyVvigebHlLeAd6Iwaz70btrSP/azhWo2KUK6QbwIQmnLaxWEKTjPAa0piC5EsE01vmZhYtzOe+o4OMOqw9f4c77yWEkyuz6pMimbd+Rs8cRCjxGPAzz5zENVF20n8oPfXGYPBNDb4Mr+9D/mDYm818hU05+OrVyuq+0RcD++v6Jex6b33s6zpjGr3U7HNShPexa9fKM1yiMEZgwtiFEWQPOnJK9wTOEVBk3w2MaW4P3Xwhe5ybI9nlDxghEcNtFxtn4gonvWS6ihaiiTO0mn5E4JKrSYR3NnJi8rFz1InDJXAu2Dlq5zGKnbf0lR6J5ZqETS5I2dWGcAf6ctheyzu/MhJ0IdqA75IE+l96UZIuTkQGZnFx3ZyV9UyrzzBktj8do1/pRkvn/9fa3ighE0QsvQrSlStXnn766SNHjmxvb589e1b/FX3X6fIP31ycqjYSEZETPwtfmiADiLiDvpqBfZgc89Db3EVYXxSjuqP/PN81571bLDz7UnZ6fAZ89CpI6/X63nvvfe21115//fXz58+/9tpre3/rrhLxI3GaxUxXczfEtPtc5X2BWTuC7Ok7X75LiXFATjUiSFlv2dU7aSyGtrN42Kbii7AXhUoIb711uXUTKHoVpP379z/wwANKqVtuueWuu+7a2dlhftFZw9ThzMBFI2tuhUw/QuSvjCp3NM2HYsInmcxlCo18F6xcJXfiSHu4+GYgIS7kz1Aibzl/aKT4fPnMWcLX649/Nmg+MqghlWVnZ+fChQtHjx7VP7QNyBczqb1rH5muhElc44nZxT8tEXM0n4ETTieb63YTBQnJ3fdBuOb5bMxrtYJf+MzSWn31RgSjEkZGOD1t+7506dK7776rlDpx4sStt96qlNrd3f37v//706dPHzlyZD7MMA7aCPTFUbq5OOP6atgX9a3+nEVd4xNiTOx9RPM429mqJh7Zt/4gGrmYxrE/0Q1Mv/XGwRI0KdHCOSGgczTkoM8F/qx3TgR6dixSenB840/HgkLoKUK6fv36tWvXrl27Nv1zd3f39OnTDz74oK5GKny4M66X297pOENP353hXBJmj0JSWC9VqqNPm/2cYuHEVXLGgV5FKZ59+lKXomybz3R3hKfseoqQdD744IPTp0+fOnXq+PHjvmOYEVLeAZAzJxUjQlqMDDgRkhGPGleXECw64znl8TL6J754yDjSuC4/u1WN9AjJiDB857ErVRIMQP+wocOTED4ePLh1+bLcfQ37WjcgkhdffPH8+fNvvvnm448/rpQ6evToM88807pRy1QwxyAZKCHGzZ2RcsnA7ImYXdZbbgeRhF8zPm/ugGYSmxGUp2rueX1q1BYJbRBOrxESh8oRkr2KbIIuSMwakv2ruAhJeVbi9b2DfX6jkYsR0vy5Uo6+2GGT8i/AhfiglPiA+cXmOjRDBOvNHV7zCOmtty43HwQfgwvSBGGIeQ10UY0aRkiGV6VFWrm8tnJpPF311U/YFiKvqGIFib/iaWgPziuG2jw/KjI+UQLCAr7l16GOByA+bz4CBD1taoggMXUefTkfpcMC40LzJ0HdJwbNSFtxaOKM9JbbN4U5LL5MnTGxo9WoHHYHF2PBjBdVe0ebn94s2hgJqyIVMv6SZaMcgwtSK5oYk3Phn+WE/JDCPoPzsNLjYzvEUCnVW253OcinLAaj2bFDk8qRAVOGa4pExFoqO0GFzOyNXN3cZSdEm50MLki2Myp9ObWUsnP+nL0NxlWCuj8Pmv7POJqna3TsoVitWHfK90muJhU6oZ04rTkXdPGubwNEQNwwa9c26InLl1RmcEGqeQ8iVtDlIOpAi1+kexHq15xLwgpp9PkHupxWxzAifpWO7+4HdTm6hU5JqIZudXKWRDNFm+Q7efMAkcPgmxoiKrGJV1TsTQ2FVo62G3JWm4MSbsTOBXpTg1FE0dtWbVY4W2X8ttyq2XlR/erZx8F5uZTe6bcsYmVjt6SaDRhtjutIdvSJY5hilgFZPI9klz94hFQZn+uZM2DVIoN0nKGSEVssrrl80UmFcTAiM7s7dfJXKdFJHNlDkwjNNkbb98Vyo2GfWU4CY7Y6Z0616HXlA0EqC+GP8rrC2WvkWn5yzNcZG3HO0yqHM1E5m+SLjXy/ynXFeRlev4pDjGqrfJEcd6xvbUiMkHxJdecx2NSwWehiUAen2jVfBoamqtuKU53L2WNSwSnk0t28qlZuTWYjs2riG8mIRsb1S44220CQ8qPL0hybV97aVAfbVdmhUmWRXqSJb7K7b6weCl0u/fxZMl0+/1vaKkoXCKPxNayEVcwTUNogOIEgZWNN/uHUCdtxpxAx0zK6Y5+r0j+0FctuQGL5PeL4hgvnCk4h7yUyVqEm9DEPnQ5duNRFnCOgsm5qIP4pHAhSA3Lt84k4Q0Z3vCi6tmLRva6c6qy5cHbuV1QFFHG6v61KR0FUyC3LHAciZZdrnUr8UzgQpJw4U1h29FDNFdLl9BJtiMgPrNO2IMqfb4tZr4w3wkiQSl4dM/e/0F+hkTkOlff7ier7IgMKUisPZdiZESjwkxJBBrT4wEH99JRzstF7/3xJjBSI8+jrhgrDYtxQ2x5S2mBXHcQ6IKJhi2a8eIagE7bFd+ulRXJNGFCQjPta3xEb/zTsz7cLYD4g+qLMeKhoHsPpMuglYeKifvGKdGMqa7bRsOz1HplOjd7vtzj+9gFiRZeJc6ePkc6NJtGltGVAQVL+R8+q4TMCO4SyS/rRmSvOdyunC4IuFDcbUyTcuG79/EmJcFCgs3bGyvYBYwRAiyzuBtpkxhSkRUobNGfd7ZyEiWV/+7pGPLQJRl8/7okm1BEbX/T9UxSLSVql3IuzGcm9i2AxdZmYJFgMNCUzpiBNXvjgwS1f0FAhrRe97o5ojC05c8Aks65rUCJG4Z9Twl6saCOR72WMgfVFqM5fTdiZOlrIc9XnKmMnSxZxJhWco9QqVxTKgILETNM3T+tFYK8rmVNdCMTiIOUuEGUJTpO60OwZwxdnqTqURm+kvmxi+kq9/Gn8sHhF4ffUmL+LuU3i6/aHwq3CyYCCJMQEg9bd6alke1OfHHPUV7g+Tcq4qUGPDhfHoa210L5DiCVnxxffdJ1risMZzOXqr67KohwCwYCCxGS+TyVuFX/dbSwb6cYsFuRlrvftTuUd8PlWGp/MVxcLvZDnr5El91H5n3ywC5zMUIlzRSV+IzVnQxNnOUUXFLuYCDMbKkjMtF40QSJkHL9Y85wPsy1YrM0FVXSYZ5uPN0ZP7CDYcBbFkiM8DsZqyZkqj9MhQ4z1zJ5kL0wkDPRmc1puV9eUZ/u4ZGHW2VBBamum6zXlmp1GtnhkLwbnjGZsfFkd35HKUuvBMLrZy+1WrrUI/+ZOOD2sfmbnFcXim60r12Pazr4w9dtYqHVRW9pQQWqOvk5UbBOROdM4jScEmE44GB/6fJDvizLhrzkWj5HsXHTi7o4vR+20mV6GQjF2x9HiEWQ/xvnfeuuy5JkCQRKEkYJYTCivtbfANTQyzvyZ/UhiO7somSxiLEeIpL/9q9DwomuMUhBnsUKfSjJEiL8ht1tBkFqh+2j9v/KnjRNaF+nYiFM+6XRYCIgymLPX440AB856yx4uXwAtx5x8NWxfROj8MNRdrG4+ndm8+wQQpAbo+5INaCOzj+9iN5GNbzHo64URYNl79vrqvo1zQJzOiAitOmVRJ2gtseOnxQVQ23GzEwah1mvsCgnNmUs2GwhSA4ysi2++2ckZ2xELMTJjA9UihMtwOlw6ilr3866gmVC3mL41QCy6otgWbvyTiBh6MQBfF6KXlZwv9mInECQR6DPK8MKc7Hlb8uoiHTDZm/SMLzZf/3LQQ+T5f4pXnO80Jqbx3bvQe2qvYIy0Z4rfLwGnSObE6AgRTCvxCq0DQZIFXfw3PhRiZ9ECEDphnA5LlH9h4ot9Ca9kL1bk625GaBm23bqtZHLSCTP66iquSc4UgqFSfZkKBEkEi2s33WHpX+ka54QhHIdvfPTz9KhPaq9EEX1UYpxpXvR7R4fIzjIqx60LHLdoEaKN3Mh79wUEqT1OF6znc5TI6RTNYkHIlpbFtZ7A9W8WiKRlj+7GiV6cd6bsfAkD/RjnOUcaJQN7DeeMlrobBAhScSLqjYvRD6fYIJa1tXPBEJ4IaRlDhOwKmeFNhtTd7D0acpQIbMmZ51F3uzEhSMWxfQrnK8bBzrxWR3ZmYPiLxY7oC95OZZiDc52hf0gEiMYPm0y/8yIa38Klu6GAINVAjwmCTGS9d0ds1yKUgi5g+twb0v/SRRTiKwOUFUEcRlDYrxlAkKRj+JoxkuOJvdiQhMxcUVOj93Rm7AiYJqXjw5gHBKkG+g6FFLMbIzk+Ri/qYCu3b9/zAMsUFfhgw2CBclDoP+puFwhScYxESqgL9u2I7Zf0Xujx4sDYZuPbFT2SwIfuZBnJGPip11F3u0CQipNiH8PYWQT0I1kTQ44Jv4MbVU7jHDz2aBhEV6YlA0ESzTB2FoHPv/S4r4HfSEKGnQmZ8bYz8LekRjyvJhzfXXYm6IYEggTkQkSHfbkevoIaGTl7CdxFf6Oxk5DKP2KDjQlxl+3Smp7D771upANBAqAGQeUB42fiK2OUsmeM/toONzRm6gimrOpqhJQdAPUg6tVdu55cDBk5OTvCkfMhR2PTgCABoRhP4fS+p4ivoPqDR/z83kgkPp02MPqG2yFTdvtaNwAAN/pkS9w63xy+gupq1F03c6EnoyZG2tsdgb2TBSk7AEAkQc/WDONcUjA2NfS1hyU7vk0944EICYgGzgjMjJSbSkefEcPMDkRIQC491opAFoxbbz9pNDaLuwp7T2L76F6Qzp07d+zYsdatAEUYZpqBIPSi/WaW0/TtPMolS6MOSN+C9P777//sZz+7cuVK64YAALJBiJCuUmOzIbGgQd81pEcfffSxxx775je/6Ttga2tr+uHy5cu1GgUAKMIwlZJQjEeDQwdhdoPy6ViQXnjhhe3t7bvvvps4BjoEwDBslBrpT+ClqJHa6waFi1NPgnTp0qV3331XKXXixIl33nnn1Vdfff7551s3CgAAMrOxsWBPgnT9+vVr165NPz/77LNXr149c+aMUmp3d/fMmTNnz5699dZbmzYQAAAyYOznNqKlgbWqJ0E6fPjw4cOHp58feuih3/3ud9PPr7zyysmTJ1ebVv4DAAyN/fcYB5aiiZ4ESefQoUPzz/v27bvvvvsaNgaAFDbB0YAInE9ijU3f274nfv3rX7duAgDx8P+8AthMNufd9iMIEgC9g9cjAR8b9b4SCBIAAMhleBHSgSAB0J7NyckAQABBAqAxG5WTAYAAggRAYyBCAExAkAAAAIgAgtQxqDcAAEYCgtQxeH4FADASEKS+wfMrAIBhgCABAAAQAQSpb/D8Ctg0YOoDA0HqGDy/AjYQlE4HBoLUMRAhsJmgdDoqECQAAAAigCABADoDpdNRgSABAHoCpdOBgSABAHoCIjQwECQAAAAigCABAAAQAQQJAACACCBIAAAARABBAgAAIAIIEgAAABFAkAAAAIgAggQAAEAEECQAAAAigCABAAAQAQQJAACACCBIAAAARABBAgAAIAIIEgAAABFAkAAAAIgAggQAAEAEECQAAAAigCABAAAQAQQJAACACCBIAAAARABBAgAAIAIIEgAAABFAkAAAAIgAggQAAEAEECQAAAAigCABAAAQAQQJAACACCBIAAAARABBAgAAIAIIEgAAABFAkOqxtbXVugl7QHtopLVHyWsS2kMjrT3y6ViQrl69+v3vf//YsWNf/vKXL1261Lo5AAAAktjXugHx/NM//dMXvvCFxx9/vHVDAAAAZKDXCOn999//7W9/e+rUqdYNAQAAkIfVer1u3YYYXn755V/84hef/OQnL168eOjQoaeeeurAgQPGMUjgAgCAweXLl1s3wUtPKbtLly69++67SqkTJ05cv379nXfe+fnPf75///4f//jHzz333JNPPmkcL3ncAQAAGPSUsrt+/fq1a9euXbumlPrEJz7x+c9/fv/+/UqpQ4cOvffee61bBwAAIImeIqTDhw8fPnx4+vmee+7553/+56tXr952221vv/325z73ubZtAwAAkEivNSSl1Llz5370ox/deeed77333vPPP//Zz362dYsAAADE07EgAQAAGImeakgAAAAGBoIEAABABBAkAAAAIoAgAQAAEMGYgvTSSy8dP378i1/84g9/+MMKl7ty5crTTz995MiR7e3ts2fPEm2wPyza1HPnzh07dkxIe+yX4bZt0k9+8pOvfvWrf/3Xf/2P//iPu7u7Ddvz4osvHj9+PPSK5Rqmt0eIbRtDNNHQvI32NLdtoz1ybDuY9XD8/ve/v//++3//+9+v1+uHH374jTfeqHDFc+fOrdfrDz/88OTJk//5n//pbIP9YdGm/u///u/f/M3fHD582Hnp+u353ve+99Of/nT+Z9smvf3221/5ylf+9Kc/rdfrp5566l//9V9btefJJ5984okn/vIv/3L6J/OK5Rpmt6e5bRtNmmho3nZ72tq20R45th3BgBHSG2+8cfjw4eklDl/72td++ctflr7i/v37H3jgAaXULbfcctddd+3s7DjbYH9YtKmPPvroY4899uGHHzovXbk99stw2zbpjjvu+Oijj/70pz8ppXZ3d//sz/6sVXtOnTqlv/WKecVyDTPaI8G2jSZNNDRvoz3NbdtojxzbjmBAQfrDH/5w++23Tz/fdtttOzs71S69s7Nz4cKFo0ePOttgf1iuqS+88ML29vbdd989/bN5ey5evPjnf/7nZ86cuf/++7/73e/+3//9X9smHThw4KGHHvrmN7/5D//wD3feeefx48dbtWe+RxPMK5ZrmNGemYa2bTeprXkb7Wlu20Z75Nh2BAMKUit2d3e/853vPPHEE3fccUfblvzmN7959dVX/+7v/q5tM3Sml+E++eST//Ef//FXf/VXzz33XNv2fPTRRy+//PL999//F3/xF//+7//+P//zP23bIxw5tq3kmTdsOyM9vcuOyac+9UMWBLAAAAJRSURBVKnpBaxKqT/+8Y+z5hdld3f39OnTDz744JEjR3xtsD8s1NRnn3326tWrZ86cmRp25syZ733ve//1X//Vqj3Kehnur371q+3t7YZD9G//9m+f/vSnH3nkEaXUXXfddfbs2S996UsN2zPDvGLNhomybSXPvGHbGRkwQrr33nsvXLgw7S155ZVXtre3S1/xgw8+ePjhh7/1rW9N2XZfG+wPCzX1oYce+u53v3vy5MmTJ0/eeuutJ0+evP322xu2Ryl1zz33XLx48erVq0qp6WW4bYdovV5PSfbpZ+ela7ZnhnnFag2TZttKnnnDtjMyYIR04MCBb3/721//+tc/85nPHDhw4Bvf+EbpK7744ovnz59/8803p7+nfvTo0WeeecZug7NhJZp66NCh+ed9+/bdd999zgtVa890rUceeeTUqVPzy3CdV6/WpAceeOCll17627/9209/+tPvvvvus88+27Y9M8wrVmuYNNtW8swbtp2RYV+u+tFHH63X61tvvVVaG+wPqzW1eXuY16rWJGntkd8wse1p3iRpt0xae5gMK0gAAAD6YsAaEgAAgB6BIAEAABABBAkAAIAIIEgAAABEAEECAAAgAggSAAAAEUCQAAAAiACCBAAAQAQQJAAAACKAIAEAABABBAkAAIAIIEgAAABEAEECAAAgAggSAAAAEUCQAAAAiACCBAAAQAQQJAAAACKAIAEAABABBAkAAIAIIEgAAABEAEECAAAgAggSAAAAEUCQAAAAiACCBAAAQAQQJAAAACKAIAEAABDB/wP28uQlbBnnuAAAAABJRU5ErkJggg==" }, "execution_count": 108, "metadata": {}, "output_type": "execute_result" } ], "source": [ "%% Filters\n", "% deciding what filter to use on your data. This is personal, decide what\n", "% is the best option. This uses the Butterworth filter \"butter\"\n", "\n", "t = 0 : 1 : 5*3600; % Time, maybe when measurements were taken\n", "y = sin(t*2*pi/3600) + randn(size(t)); % Nice cyclic data, with noise added\n", "\n", "% Setup some values for our filters\n", "T = t(2) - t(1); % 1 second\n", "samplingf = 1/T; % 1 Hz\n", "cutoff = 1/(60); % 1/60 Hz, 1/1 min cutoof frequencies higher than this\n", "\n", "% First have to setup filter, with parameters A,B\n", "[B A] = butter(2, cutoff/(0.5*samplingf), 'low');\n", "% Filter the y data, save is as fy\n", "fy = filtfilt(B, A, y);\n", "\n", "%% Plot the results\n", "% Original noisy data\n", "figure\n", "plot(t,y, 'b-x');\n", "hold on;\n", "% Filtered data, (red line)\n", "plot(t,fy,'r-o');\n", "hold on;\n", "% Original, non-noisy function (yellow line)\n", "plot(t, sin(t*2*pi/3600), 'y');\n" ] }, { "cell_type": "code", "execution_count": 109, "id": "e7815538-b035-4b5a-95ac-ae420e961c30", "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjAAAAGkCAIAAACgjIjwAAAAB3RJTUUH6QkREDoGyMZqXQAAIABJREFUeJzt3X10VOWdwPHnZhIoYqEIYSObWqGM40TjgnCOUxSCyqIpRl1hKxsGRLTZnrJAtLKsVoTVPT1bPT3WrWDrAd9QES0retxVsRs0ZWPKLkFwmyaE0DhKNIvJSkNIMsPM3T9ucrlz752XTOblmcn3czg698kzd37PfXl+89x7515FVVUBAECm5WU6AAAAhCAhAQAkQUICAEiBhAQAkAIJCQAgBRISAEAKJCQAgBRISAAAKZCQAABSICEBAKRAQgIASIGEBACQAgkJACAFEhIAQAokJACAFEhIwDlvvfVWR0fH2bNnd+3adfDgwaTP/+DBg7t27Tp79mzS5wzkABISRq5Cgx07djQ1NVVUVPz0pz/t6+tbunTps88+K4Q4ceLERx99lKxPfPbZZ5cuXdrX12cN48ILL7zpppt+/vOfR59DcuMBpJKf6QCAjPnyyy9dLteCBQuEENOmTbv00kvr6+tLS0tDoZBe54EHHvD7/Tt37tRLQqHQ2bNnR40aZZpbKBQKhUL5+fm2dc6ePZuXZ//978svv3Q6nQsXLvzss8/uueee999/f8+ePcbZ9vX1nXfeeVHiMVYAspgKjFRCiBUrVuiT3d3dQojVq1frLzZv3qzvKa+88oqqqo899lhBQYEQ4qKLLqqrq9Pf9dRTTxUVFS1evNi2zqpVq4QQY8eOraioEEJ0d3ebwli6dKn2esOGDUKI2tpabXLx4sUOh0MIMWnSpJqaGms8pgrpWGpAynDIDiPasWPHtm/fvn379k8//dT611tuuaWgoGD+/Pnvvffe/PnzDx48uH79+vvvv7+np+eiiy6666679Jr33nvvD37wg3Xr1lnr7Nu375lnnvF6vb/+9a99Pl/0eObPny+EOHr0qDb5ox/96KOPPmpvb8/Pz3/iiSdM8VgrJGmpAJnBITuMaAcPHtR6/29961sej8f01xkzZuTl5RUVFWmH9d566y0hRG1t7dGjRzs6OlpaWvSay5Yt27RpkxBi+/btpjrHjx8XQlRVVc2dO/f48eOrV6+OGVUwGNReHDt2bO/evT09PaFQqKenxxSPtcKwlweQSSQkjGi33377888/r70+ffp0PG+ZN2/e9OnTv/vd7xoLR48eHamO3+8XQmjnlrT/RvHhhx8KIaZNmyaEeOmll1asWLFmzZply5bZDq1iVgCyC4fsgGjy8vJaWlr27t376aefzpgxQwjx2WefLVy4sLS0tLCw0FrfWkfLLr/85S/37t379NNP236Kz+fbuXPn2rVrH3744Xnz5mkDoGPHjgkhFi9ePGHChMbGRms8thWALJbpk1hAxohYFzWoqvrggw9qe8rLL7+squpTTz01YcIErWT9+vXGd+nzsdZZuXKlEGLChAlr1qwRdhc1CCEcDofT6dywYUNPT49W3tLSMmXKFCGE2+1euXLlggULTPHYVgCyl6IO7g8AbGm/YzUebfP7/fn5+ZEu47ato132Hf0tkeZjusTcFI+1ApClSEgAAClwDgkAIAUSEgBACiQkAIAUcvl3SC6XK9MhAIBEmpubMx1CNLmckISUS9/lckkYlYn8QcofociGIOWPUCQ7SEUZeJHEy7myaDFKHiqH7AAAUiAhAQCkQEICMLKoQoldCZmQyz+MlfxoKYD0UxShCkUROdzzRSN5r8gICQAgBRISAEAKJCQAgBRISAAAKZCQAABSICEBAKRAQgIASIGEBACQAgkJACAFEhIAQAokJACAFDKfkILBYH9/v7GkoaGhra0tSknMCgCArJPhB/S1tbWtW7du1qxZDz30kFbi9XpdLpfP5ysrK/N6vdaSmBUy2R4AQKIynJC2bt26fPnyxsZGbbKurq64uHjjxo1CiPLy8srKyvr6emPJ1KlTo1eorKzMy8v8sA8AMFQZ7rsfffTR8ePH65Pt7e0ej0d77XQ6W1paTCU+ny96hZaWljSGDwBIGrkGE01NTQUFBdprh8PR2dlpKjl8+HD0Cp2dncYZugalqwUAIBeXQaZjiSHDh+xM3G53X1+f9trv9xcWFppKPB5P9AqFhYXGGcr8KCoASANjNyh5TpJrhOR2u2tqaoQQgUCgtbXV6XSaSkpLS6NXcDqdmW0CACAxco2QSkpKioqKqqqqurq6qqurrSUxK2S6BQCABCkSPlk+GAwqimK8WM5UErOCRvKnxwNIP0URqlAUIWHPlw6S94pyjZA0DocjeknMCgCArCPXOSQAwIhFQgIASIGEBACQAgkJACAFEhIAQAokJACAFEhIAAApkJAAAFIgIQEApEBCAgBIgYQEAJACCQkAIAUSEgBACiQkAIAU5EpIoVDo0KFDf/zjH42FDQ0NbW1tUUqsFQAAWUei5yH19vbeeeedc+bM+fzzz/Pz8x955BEhhNfrdblcPp+vrKzM6/VaS6wVAADZSKKEdODAgYsvvnjt2rVCiO985zuPPPJIXV1dcXHxxo0bhRDl5eWVlZX19fXGkqlTp5oqmB4aCwDIFhJ133PmzPn000/37du3devWxYsXCyHa29s9Ho/2V6fT2dLSYirx+XymChmJHAAwfBKNkAoKCq6++urt27efOnXq4YcfFkI0NTXNnDlT+6vD4ejs7DSVHD58eO7cucYKpnm6XC7thcyPkQeA1NG7QflJlJA++OCDU6dOvfjii3/6059Wrlz5q1/9yu129/X1aX/1+/2FhYWmEo/HY6pgmid5CMAIZ+wGJU9OEh2ya21tdTqdQohx48ZNnDjxq6++crvdNTU1QohAIKD91VRSWlpqqpDZJgAAEibRCKmiouK+++47cOBAd3f3tGnTtOxSVFRUVVXV1dVVXV0thCgpKTGWmCYz3QIAQOIUVVUzHUOYYDCoKIrxYrmYJdYKGpfLxSE7AEaKIlShKEK2ni9NJO8VJRohaRwOx1BLrBUAAFlHonNIAICRjIQEAJACCQkAIAUSEgBACiQkAIAUSEgAACmQkAAAUiAhAQCkQEICAEiBhAQAkAIJCQAgBRISAEAKJCQgJRRFCEVRFKEomQ4FyBLS3e0bwMhhzNYj83kQMJJuhBQMBuvr61taWvSShoaGtrY2Yx1TibUCACDryDVCeu+997Zv33711VdPmTJFe2Ks1+t1uVw+n6+srMzr9VpLrBUApMgIf7odUk2ihNTb2/vkk0/+67/+q/7Avbq6uuLi4o0bNwohysvLKysr6+vrjSVTp041VbA+NxaA5LQkl+kokHkSdd/79++fO3duTU3Nnj17Tp48KYRob2/3eDzaX51OZ0tLi6nE5/OZKmQkcgDA8EmUkE6fPn3gwIG8vLwxY8asWrUqGAw2NTUVFBRof3U4HJ2dnaaSw4cPmyqY5ukalLZWAIBUXAaZjiUGiQ7ZKYpy8803X3/99UKI3/zmN//1X//ldrv7+vq0v/r9/sLCQlOJx+MxVTDNs7m5OV3hA4CMjN2g5DlJohGS2+2uq6vTXre3t1900UVut7umpkYIEQgEWltbnU6nqaS0tNRUIYPxI3W0X/Pwmx752a+pqKtNFaxUDJBohORyudxu9/e//31t6DNlypQpU6YUFRVVVVV1dXVVV1cLIUpKSowlpslMtwApNBLOe+v9tvzXsClKxCBHwppCiiiyXb8ZCoWEEMaL5YLBoKIoUUqsFTQul4tDdrlBv9pYZENnrUkgZvkTUsxG2VeInL60+kKItF1KPsKvXJe8V5RohKSx5hX9KvBIJdYKQJZieIGkkP/LjS2JziEB0XGyAchtJCQAyEHZ+AWOhASkUDZ2CgnjGkgMk3TnkAAge3H/8uFghAQAkAIJCUDmjahjm4iEQ3YAkGQpvXw/h48KMkICgOyTk2NKEhKA9Br+1XjDnkMSe/OBG/cpChcZDh8JCUBccqzDzbHm5AYSEgAkU04eTEsPEhIAQAokJACpoj8YieNjiAcJCUCuIx9mCRkTUkdHx8mTJ/XJhoaGtrY2YwVTibUCACDrSPfD2N7e3pUrV951111LliwRQni9XpfL5fP5ysrKvF6vtcRaAQCQjaRLSD/5yU+uu+467XVdXV1xcfHGjRuFEOXl5ZWVlfX19caSqVOnmipYn+8HABoO3UlOru573759kyZNKikp0Sbb29s9Ho/22ul0trS0mEp8Pp+pQvpjBpKCvhKQKCF1dna+/PLLa9as0UuampoKCgq01w6Ho7Oz01Ry+PBhUwXTPF2DUh8+kKNIlVnOZZDpWGKQ6JDd448/fvnll7/55ptHjhzp7e2dPn262+3u6+vT/ur3+wsLC00lHo/HVME0z+bm5rTFDwASMnaDkuckiUZIt9566/Tp0wsKCvLy8vLz8wsKCtxud01NjRAiEAi0trY6nU5TSWlpqalChtsAIOkSGKIxqstOEo2QZs+erb/u7e297LLLhBBFRUVVVVVdXV3V1dVCiJKSEmOJaTJjoQMYwbhXULIoqvTP0wgGg4qiGC+fM5VYK2hcLheH7HKDogzs84qQf4MdMNSYte/02nN05GqjomhP3TFGKOyexKNXEIOttm/U4AyFYSkJ2zla6kcsif4WSxNE8rYl04pOw7Yaz6aihWGtIHmvKNEIKRKHwxG9xFoBADRa/oyZwiADic4hAcgs7bk+qXu0D7e2Q3QkJABSIEuBhATA3sga0IyIRsqOhAQAkAIJCQAgBRISgBGBY3LyIyEBAKRAQgKAnJK9Y0ESEgBACiQkQBY5cEu0pDdBv+5c/9Y/zCvRc2Ah5zASEoCRbgT93EpuJCQAuSMlP+Y1DtCQSiQkYBjop4DkISEBiRhZt9XJaaxBeciVkEKh0O9+97uOjg5jYUNDQ1tbW5QSawVgpMjy3jSHLzHI4aaljkTPQ2ppadm8efP06dOPHTt2/fXXr1q1Sgjh9XpdLpfP5ysrK/N6vdYSawUgmwzjQT086ScNBh5HyHJOC4kSktPp3LFjR15eXigUmj9//qpVq+rq6oqLizdu3CiEKC8vr6ysrK+vN5ZMnTrVVMH63FgAQFaQq/vW0smJEyemTZsmhGhvb/d4PNqfnE5nS0uLqcTn85kqZCJqACNMfEdKs/x4agZINELS9Pf3b9q06f777xdCNDU1zZw5Uyt3OBydnZ2mksOHD8+dO9dYwTQ3l8ulvZD5MfIAkDp6Nyg/uRJSIBC455577r77bm0Jut3uvr4+7U9+v7+wsNBU4vF4TBVMMyQPIUspilAHz1+cO3XBeQwMnbEblDw5SXTILhgMVldXe73eOXPmaCVut7umpkYIEQgEWltbnU6nqaS0tNRUYUifqCgDV+8yskZWS/MGLPn1Y5KHhygkGiE98cQTjY2Nzz333HPPPSeEqK6uLikpKSoqqqqq6urqqq6uFkKYSqwVANmoQlEUVd6BjYSjLglDQlrIvKMMCAaDiqIYL58zlVgraFwu19GjzdpxD9tWakdFhBCRKkASEq4pbVASz9YlIoetH5dT1bD6QlXjOWRnXyeBlgy+2bScTW3UJyO1yNrkSHM493bD4C6sMDwkYxuNC82+NeERGucgTOMnVRWxVqX9thd1ocVc78MUz7YXqYLL5ZL5RIZEI6RIHA5H9BJrhYRF39aBTBnsYhg8IJdJdA5Jfto5J044IdVSuI3JtPnKFAukkOMJidObiIjuEEPHPQxTKpcT0tGjyThUynYHwA59Q9LlckKCnHLmyCfj74xgsecwElI0jM3ThEWc1Vh9SBISEjJMz/rp7NZG8lcNY6vjXw4jc1klB8subiSkXCP5ATFrb4gUsd6IxPpCNlwyMMKRkMKc+0Vb+P6Qtt1D8nSSOqYTA9m+HIYZ+ZBOk6TznIp1NaU5f2TvJoF4kJBGomzv7lNiGIsjsbcOc/mz+pJlhFwlMTBiltuISEjJ3eBUkfKuPKUfEXaLGgCSSen3xeT8GCZlRkRCSpjWcUfrvtN8Ll7/0KTVSm0M2fM5Q5P4kCjujiaBLw1yDnxH+Lef9B3SDP+MLF3sIyUhxdwgEsgsST9FHP98ZOt0kn6dVhJ2p/CPG1i/Efrr5PQa0q2VdBjqmjLWT3WnGf9WGU8k6e/iw3qY8E3X5pxrrPZmRYoaKQkp6aKt+1R2TPbzTm9XmED3bX2LdfcYuGVyMsac6TzTroctogavlyfxeOyQupjEPjTKQYKY3V/0Cta/RmmO/qeYTY5U0/6NimJKkNG/ZWbdV45zG2eWhD7yElLcKyZKRRm+awxhA0vxWN6mL07Z1p++JW/omYzXT9u2LLGoVKHE370aSyIt3UhdcJTePPZB6bgDi+cgRKQPUoRqejyE1kzTW0xfaAYmFfsKcUZlG6eWqKIv56EmaeNuYh3xxBwDCbtVbyqJMySZc1OOJKSGhoa2trbY9eJbFQn3evFeBTvsk0DGwUTMfcD4Krnbov3OMOzhSfSvw/Y7mN2vbSLN1vTesF5POTfcsS5cLYVEyUkDOcZSw9qFWfvWBJKTuSV6NcuRSdPHWVthuyFFyR9DDc+aLONpb8w6xjVlG388M4kUajzrJWYKsd0wIqU0a4KxZqDowQ9p3CmhLHgeUkxer9flcvl8vrKyMq/XG6maou9GimJ64pHx+V3DZPwlk+3TxsTgU22E4YFj0WhdpOH5ZqYPMlUzVlCEqjdMsXuOzkChMRpz9EMTthyHPpOhZjHjJ6jC/sGsg8s8LBnY9ggDk4oiwjcG+44j1gYT1yGjOAy0y3xkyVwypGAGSiJ3XoMfEcdMDC+iVIgnyOh1wtagUK2ZZqgLYeChf8aH6Q09e8X8lOgV9I1TNfzVvBeHr+t4vrtkr6xPSHV1dcXFxRs3bhRClJeXV1ZWao+Ota428+aYQN8X+ZCIopi7sMFkMNDHmbb+wSFOWKHNzBXDZBw5U9+arR2ENUhDuwZehMcW7VOMb4/Cur/ZVoju3NBH7zLOhSqMjbK20fYL41BzSfQ6WlRRatpWiLMzjdRLWntn2wqmrS6K+CO0f/RqeMwifHxjDsPynNbosQ18fVTC5hbWalU17p42rbatoIZtWjYfGt6E6CkwsXwWbd+Pb1KIsH3AmMDiX/vyyPpDdu3t7R6PR3vtdDpbWlpMFWzXinEsb9p/ok8a52Cav6lO9MlIhYpQ4xxVGA+72wZm+ydrSMYcaZ00lkepEKk5Q/rEKHNIrFHWmG03hqEOjqN0xKa+MtKcrZ172KSq6v9UNayOtnmErXpVNc7HGIM1NuMbz02qqvGvpprWEi0041/DnqouVL2C6RPPtSvikjU3X/t3ySUu42fpfzF+uqoOvNG6uPQK+qQehXHxanMwThoXkb5IjQvN1OpIy9y4OiKteqs4K8T8JqGXaGFfcolL5icOZ/0IqampaebMmdprh8PR2dkp9LMYQhVCXHKJSzlq2OZMh31VIcLHHqoa/qVbFUKoprdYJtVIf42T/q6BbcXwFTJKTUUJ251ijgb07VJfDpEq6JPWEu2Fbdq2nYOx/jAraIs6ZsyRZqiqQhtRaXMwHOsb/I9l2zBvKtEHjgMjzcFoRdib9aOLlnkOfnlXB7o962z1Q5HquRL988I+V39tfNe5qsJc0+VyCdEcuYey/kW1/7P+SeEFsReX6SOGFIpdXUOh5ZBjxPqmdthEZo323LpTxSXnFqNxyzS92zAuVAe2OuO+bNze9G1Vn7RuPzF7nsH6A3UuucTlctksAalk/QjJ7Xb39fVpr/1+f2FhoRj8jqV9F2hubjZ9MzN+/TJNWktiviXmDE3/tKgifaiR7duNNWN+rjB8nVQHv/EZG2X8zqtXEGIgSP0rnv4RphLbCsLuO6xpDolVEOJco/Qveqb6QrX5Hq2/XW+1aVHbLvOYW4LxQ401zhUa82KEbUlvlDDVC48qZolljrZzstRMVFjMUauZGptjhtou06IwbmlC2GwYpgrWwig9gLV+c7PUt2kQQtidBc4qjY2NW7Zs2bJlSyAQqKioeOedd/Q/ZcH3AQBIL5nTUtYfsispKSkqKqqqqurq6qqurjb+SeblDgAwyfoRkiYYDCqKol1fBwDIRjmSkAAA2Y4hBQBACiQkAIAUHJs3b850DLkmGAwGAoH8/IELRkKh0EcffZSfnz927Fit5MyZM4FAIBAIhEIhrVpDQ0N/f/83vvENeYLU6hw4cKC3t3fixIkSBtnf3+/3+wOD8vLy8vLy0hxkzMWolfT19U2YMEErkW0x2pakM8hQKHTgwIGCgoLzzz9fL7QGYCpJ82KMM0jTopZtMVrrpH9rjI6ElGRtbW133nnn8ePHy8rKhBBnzpxZsWLFV199tWvXrv7+/ssvv1wIsXDhwsbGxpqami+//HLGjBler7erq+utt97605/+dMUVV0gS5Hvvvffggw86HI4zZ86UlJRIGOTu3btfeOGFmpqampqaf/7nf54zZ84999yTziBjRtjb27tixQpVVT/44IMPPvjguuuuk3AxWkvSGWRLS8u6des6Ozt37tx56tQp7Ufu1gBMJWlejHEGaVrUsi1Ga530b42xqUiq9evXv/baa//4j/+oTb7xxhvbt2/XXt90003BYPDs2bM/+MEP9Pr/+Z//uWHDBu31jTfeGAwGZQjyzJkzN99889mzZ2UOUq/5xRdf3H777ekPMmaE77//vh6Sx+ORczGaSvbv35/mILWPCAaDc+fOVe02NlNJ+iOMJ0g1fFGnf13HE6GxTka2xpg4h5Rkjz766Pjx4/VJv9+vTxYVFTU1NXV1dZ06deqdd96pq6sTcdyLLyNB7t+/f+7cuTU1NXv27Dl58qScQep/euGFF5YvX57+IGNGOGfOnE8//XTfvn1bt25dvHixnIvRVPLJJ5+kOUjt1xonTpyYNm2asNsjTCU+ny/9izFmkCJ8Uad/XccTobFORrbGmLL+h7GSmzVr1tq1ax0OxyeffHL8+HGHw5Gfn3/bbbcFg8Ha2tpt27ZNmzbNei++jAd5+vTpAwcOzJw5c8yYMatWrbrtttsmT54sW5Ba+ZkzZ95///3169f/0z/9U2aXpDXCgoKCq6++evv27adOnXr44Yf/7d/+TcJ1bSopLi6ePXt2moPs7+/ftGnT/fffL+zuTmkqOXz48Ny5c9McYcwgTZVjVshUhHqd1157LeNboxUjpNSaOnXqM8888/Wvf33JkiUTJkyYPn36hAkTlixZsmjRon/4h3/o6+ubO3eu9V58GQ9SUZSbb775+uuvv+GGGy699NJQKCRhkFr5jh07vve974kIdzXMbIQffPDBqVOnXnzxxZdeeumRRx654oorJFyMppK/+qu/SnOQgUDgnnvuufvuu7V7fVnXo6lk5syZ6V+MMYM01U//1hhPhMY6Gd9fbJGQUq6wsPD666/v6OiYPHmy/r1eCBEIBLq6usrKympqarTJ1tZWp9MpQ5But1s7oiiEaG9vnzdvnoRBCiGCweCePXuWLl0qhHC73RkP0hShHsa4ceMmTpx48cUXZzxCa5Cmkvz8/HQGGQwGq6urvV7vnDlztBLrejSVlJaWpnkxxhOk6S1p3hrjidBUR4b9xYpDdim3fv16v9/f3d392GOPCSHq6uq2bNkyefLkzz//fPXq1UKISPfiy2CQ2heo73//+36/3+PxOJ1OCYMUQrz66qs33njj6NGjRdS7GmYqwoqKivvuu+/AgQPd3d3Tpk274oorMh6hNUhTycSJE9MZ5BNPPNHY2Pjcc88999xzQojq6mrrejSVpH9FxxOkSZqDjCdCax0ZtkYTbh2Ucrb32QsEAgUFBdHrpJNtAKFQSAyeCI1UJ53iCSCzQdp+uqlQwsUYT0maxQwp4xHGE0PGg5Q/QhMSEgBACrIkRgDACEdCAgBIgYQEAJACCQkAIAUSEgBACiQkAIAUSEgAACmQkAAAUiAhAQCkQEICAEiBhAQAkAIJCQAgBRISAEAKJCQAgBRISAAAKZCQAABSICEBAKRAQgIASIGEBACQAgkJACAFEhIAQAokJACAFEhIAAApkJAAAFIgIQEApEBCAgBIgYQEAJACCQkAIAUSEmDvrbfe6ujoGM57z549u2vXroMHDyY3MM3Bgwd37dp19uzZVMwcyAgSErLM3/3d3xWGS8WnNDU1VVRU/PSnP42nsjGYHTt26O/t6+tbunTps88+q1U7ceLERx99lKwIn3322aVLl/b19ZnCuPDCC2+66aaf//zn0d+e3GCApMjPdADA0Fx11VXai+7u7hdeeMHj8WiToVDo7Nmzo0aNMla2LdTKQ6FQfn6+EOLs2bN5eXl5eWFfzi699NL6+vrS0lLbt5h8+eWXLpdrwYIFQohp06bp7w2FQsZqDzzwgN/v37lzZ/yx2dbRArYNw+l0Lly48LPPPrvnnnvef//9PXv2GGfb19d33nnn2QZjrQBkgApkpxUrVowdO7a1tVVV1Z/97GdjxowRQkyZMuXdd9/VKlgLu7u7hRCrV68ePXq0EKK8vPzuu+92OBwOh+NnP/uZceZ6Te3FmjVrtFnNmTOnt7fXFIkQYsWKFZHeu3r1alVVN2/erO90r7zyymOPPVZQUCCEuOiii+rq6vR3PfXUU0VFRYsXL1ZV1VpHVdVVq1YJIcaOHVtRUSGE6O7uNoaxdOlS7fWGDRuEELW1tdrk4sWLHQ6HEGLSpEk1NTWmYKwVkrSKgKEhISErvfrqq0KIbdu2qapaX1+vpYR33333mmuuGT9+fE9Pj22h1uk7nc633377lltu0Xrwd9991+l0jh8/3jh/U1KZNm3au+++u2bNGiHE7t27TcFoiWrbtm3btm3z+Xy2CenQoUMFBQXz589/77339u7dK4R46KGHenp6rrnmGrfbrX/imDFjNm/eXFtb+9///d/WOjU1NUIIr9f79ttv/8Vf/EWUhPT222/ry0dV1bq6uo8//ri9vb2oqOiWW24xBvPFF19YK6RorQHRccgO2aejo+OHP/zhokWL7rrrLiHE//zP/wghVq5BlkBPAAAXzklEQVRcee211548edLr9dbW1p44ccJaeM011wghFixYcOONN548efKNN95Yvnz5woULZ82a9frrr0f5xAULFixcuFAI8Ytf/KKnp8da4eDBg0ePHhVCfOtb39KPIhrNmDEjLy+vqKhowYIF27dvF0LU1tYePXq0o6OjpaVFr7Zs2bJNmzYJIWzrHD9+XAhRVVU1d+7c48ePr169OvqCCgaD2otjx47t3bu3p6cnFAr19PQYg7GtEH22QIqQkJB97rjjjry8PK3L1umdrxBCP3ljW2g8AWM9hWNLO3Rme+ZGc/vttz///PPa69OnT8czz3nz5k2fPv273/2usVA7lhi9jnZ6yfZslu7DDz8UQkybNk0I8dJLL61YsWLNmjXLli3z+XzWyjErAOlBQkKWefLJJ999990VK1a8//77Wsns2bOFEDt27AiFQlu3bh07duz8+fO1wYSp0HSVQZrl5eW1tLTs3bv3wgsvFEJ89tlnP/zhDz///PMvvvjCWnnGjBnWOlqC+eUvf9nd3f30009b3+Xz+Xbu3Pnhhx/+4he/mDdvnjYAOnbsmBBi8eLFX/va1xobG6+++mpjMG6327YCkH5c9o0so51HeeGFF5YO+va3v/3UU0+9/vrrN9xww/Hjx1999dXzzjvv8ssvtxZmNvIf/ehHBw8evOGGG06dOqXFVlRUNHPmTK1FJrNmzbLWufbaa1euXKm1fc6cOdZ31dXVLV++/J133tmwYYN2GkkIsWzZsilTpsyfP//OO++8/fbbTcHs37/ftgKQfoqqqpmOAUgOv99vPQRnW5gp2u9Y9aNtfr8/Pz8/ypFA2zq216nHZF0O1mDkWVAYmUhIAAApcMgOACAFEhIAQAokJACAFEhIAAAp5PLvkFwuV6ZDwMhy9GizKhRFqEKISy5h84N0mpubMx1CNLmckIRh6btcLsnXRHLR3oxQlHOvUxePJI1NDxqbdR8xHNIdsgsGg/39/frkmUHGwoaGhra2tkiTAIBsJNcIqa2tbd26dbNmzXrooYe0kkWLFmn3UJkxY8Ydd9whhPB6vS6Xy+fzlZWVeb1e02QmowcADINcCWnr1q3Lly9vbGzUJoPB4KWXXvr444/rFerq6oqLizdu3CiEKC8vnzp1qnGysrIy0s/XZR6lpgLtzWE0NleNqMbakuuQ3aOPPjp+/Hh9squr69SpU++8805dXZ1W0t7ert/b3+l0+nw+46TxNv4AgOwiV0Iyyc/Pv+2224LBYG1trfagzKamJu1BAEIIh8Nx+PBh42RnZ6dpDq5B6QwbI5wqFBF+gQOQQS6DTMcSg1yH7EwmTJiwZMkSIcSiRYsqKyvPnDnjdrv7+vq0v/r9fo/HY5wsLCw0zYEhMIARztgNSp6TpB4h6QKBQFdX1+jRo91ut3Yf/kAg0NraWlpaapx0Op2ZjhQAkCCpR0h1dXVbtmyZPHny559/vnr1aofDUVJSUlRUVFVV1dXVVV1dbZrMdLwAgMRlweMnAoGAfqJIEwwGFUXRL6gzTeok/wkYco+iDJ5AEqoQQvp9CyOO5L2i1CMkjSkbCSEcDkeUSQBANsqOc0iA/LisDhgmEhIAQAokJCDJFKFqZ5IADAkJCQAgBRISAEAKJCQAgBRISAAAKZCQAABSICEBAKRAQgIASIGEBACQAgkJACAFEhIAQAokJCBVuN0qMCQkJACAFEhIQPJxf1UgASQkAIAUSEgAACmQkAAAUiAhAUmjCkURaqajALIVCQkAIAUSEgBACiQkAIAUSEgAACmQkAAAUiAhAUmmcp0dkBASEgBACiQkAIAUSEgAACmQkAAAUiAhAQCkQEICAEiBhAQAkAIJCQAgBRISAEAK0iWkYDDY399vLGloaGhra4tSYq0AAMg6+ZkOIExbW9u6detmzZr10EMPaSVer9flcvl8vrKyMq/Xay2xVgAAZCO5EtLWrVuXL1/e2NioTdbV1RUXF2/cuFEIUV5eXllZWV9fbyyZOnWqqUJennRjPgBAPOTqvh999NHx48frk+3t7R6PR3vtdDpbWlpMJT6fz1QhzQEDAJJFrhGSSVNT08yZM7XXDoejs7PTVHL48OG5c+caK5jm4HK5tBfNzc1pCRkA5KJ3g/KTOiG53e6+vj7ttd/vLywsNJV4PB5TBdMcyENID0XJdARABMZuUPLkJNchOxO3211TUyOECAQCra2tTqfTVFJaWmqqkOGIAQCJknqEVFJSUlRUVFVV1dXVVV1dbS2xVgAAZClFlf7xlsFgUFEU4+VzphJrBY3L5eKQHdJDO2SnCkURqqr9b/A1IA/Je0WpR0gah8MRvcRaAQCQdaQ+hwQAGDlISAAAKZCQAABSICEBKcTvk4D4kZAAAFIgIUEmDCiAEYyEBKSEIlRVkF+BISAhAQCkQEICAEiBhAQAkAIJCQAgBRISAEAKJCQAgBRISEAy8bwJIGEkJACAFEhIAAApkJAAAFIgIQEApEBCAgBIgYQEJB/X2gEJICEBAKRAQoJksvaRSKpQFMHICEgcCQkAIAUSEgBACiQkAIAUSEgAACmQkAAAUiAhQRb65XVZe50dgGEhIQEApEBCAgBIgYQEAJACCQkAIAUSEgBACiQkAIAUSEgAACmQkCAL7pYNjHCyJ6Qzg/r7+/XChoaGtrY2YzVrCQAgu+RnOoAYFi1aNGPGDCHEjBkz7rjjDiGE1+t1uVw+n6+srMzr9dqWINspCg9dBUYcqRNSMBi89NJLH3/8cb2krq6uuLh448aNQojy8vLKysr6+npTSV6e7MM+jDgkWCAOUiekrq6uU6dOvfPOO+PGjZszZ44Qor293ePxaH91Op0tLS3WEpfLpc9Bf93c3Jze2DGCcPM9yMzYJUpO6sFEfn7+bbfdFgwGa2trV61aJYRoamoqKCjQ/upwODo7O60lxjk0D0pz5ICOdIXMajbIdCwxSD1CmjBhwpIlS4QQixYtqqysPHPmjNvt7uvr0/7q9/sLCwutJRkLFwAwDFKPkHSBQKCrq2v06NFut7umpkYraW1tdTqd1pJMBwsASITUI6S6urotW7ZMnjz5888/X716tcPhKCkpKSoqqqqq6urqqq6uFkJYS5CN+AkSAEWV/uKfQCCgnyXSBINBRVGMV9NZS4QQLpdL/mOm0CjKwA9j9Z/HSr9hnqOdJRqIXz1XeK4tg03LplYhF0neK0o9QtKYspEQwuFwxCwBAGSX7DiHBGSjgVERgPiQkIDUIicBcSIhAenDb5KAKEhIQBKYrmgw0Y/dMVoCoiAhAWnCle1AdCQkyCVnLgSIOFqyNo4DeYAQgoSErEGvDeQ6EhIwIpHgIR8SErKf9H1rlMN3A7FL3wQgDUhIQKqE5SHDhDH7qNzGDxhEQgJSzJiKDJdsmAZFjJEAEhIQ2fCzhCEbWQ/c5cb1hECykJCADNAzXUaO1+knrhiWQSokJGQBKc78J+nTtaN22j89Gw08d4P0gJGNhATZGR8slO2MR+1sH/t07rq7lMml5YkcQ0JCNrG5EEBRrOVJ/KyYc07goxWhGrORnpPOjZwYKmFEIiFBUvF3yqm421DYE2yVaCddhvrRqnruX9gnClX/p+ek9AyYAEmQkJCtIvakQ+1iE+6SB/PT8I+AWfOTGEyKqRsw5cxtA5EzSEiQmtYpx1kzPSIdyot0O4b42aYl/UP1f0CuIiFh5NL794THVCm6OkC/Ek8/gqf/KbkDJgZJkEp+pgMAEmftTAfyg3b2ZyhMt/MRhjRj22UbC6M8mi8BqjqYJu1iGEwhycyCQ19aQEqQkCCdgbP6EfrcKL2n6UqE6J2s8fSP/nF6v68/AVZRwuZi/vGQmpKrp9XwC8JNMQhSCHIUCQly0ccHA5MRMlM8dzoYSDkR/m6as5aKrHML++VQki5hGCprDAkLW7aGRU2Gy33ZsI5JSMg8ayc7OMoJ24FMI6eYiWEgxyj2vz9VByuYhkExxlWGgVRm9+7o6XYI84k6HgXSiYSE7GAaOVnLzSd1DEfhFKGeO5qnht2qYEipJWxIkaF8pKVPfTw3zLSktehc+icr5Sj93OpRub95kJAgKeMFYNbxjZF1iGOuMNh9a92u+c1Dca4Hz9x+rbc3sbRkPlYZfuAu/vkgW2TRTwVISJCX7UXJtid7VJssM1BmzEDmg34J9bwy9Nd6XhThaSn+2Kz3iTA8qMnmICeQBiQkDDCd7k7unIc0w7AhjmUoYjwQF89sI12VkO0drha/fsBNPzIZpXKsWZ07OSdIS7kl/Bi1vCuVhDRyRRnIW4+JDWe2NgfKzOd+ol3bFl54bsiUQF+pn4DJmX7WmEtsfpWlX6Eex0GbSBkusbWfM0s425m+aLpclzRnLpiYSEgji+3JlWj1B39mGvYuu77GePMC43ttP+jcaY/BnwFFmbOJqn+PT8hAVsstpiHOQKFQoqSomLOKOeoSwn70KcKPHw7nyw3iYXsXqyw6b2REQsp9UZJQ7O7JkkisMzTVN16WHelTTH8yvAWJs1z6YbM447zC23TdRNhsB79ARMpV+ok62+8ocaa36FJ3eDmDYl5REn2/O1fN8G0vrE42LCgSUu6I/p0oUnqIdTsD85fueCIJO45kdy8D274SyWV7aaLtUwEjvV2/biKs3HJXi7BPGfyGYfu9J8r2E0+ustzeSSuMeIlKYmfUkiueBGysE+debP6UCEcjwsuzQC4npKNHm5MybjXtDLbD4Xg28TQMoqPv8OeqJXQtlm33FOUtUT4libceQJwSO+UmLNekaD9XirQxGDpZ8/ee6NtP9O86pgsuhF2SszlPGSFlRvpBmxh6DhvS1ms+3qDYZ5Eooi1AwylAY2F2XcqfuwkpGV/CTUeojAcibMfIMZl+rRljJ7Td4Yf+loE3DntxpGibzpZdZcSyu6Q+9lqzVoh+MCqem24Iy/cq47Xv1p3LdngX6byaaTePfwePclTT9oNsI4n/ThnxXDAZf33Z5G5CUuNNEtHmYTmDYtwxoh+7sBV2Aj++O99Y5xAj5qza/gAR+7hxxPG9aVgQ9gtou4OTkc6rhdWxJLNIe6LtT+JE5F3beIwuypHz7BrTJFeOPA+poaGhra3NVHjJJS518FnRxtdx/hPaNyPDP2OJMF5TG14tyj9TZSEifnSk2cYTthDC5XKleplLZUS1d6Q11nYLNzLtsNoL28rx7+bW2Zr3RO14id0uGWnXjh6GGFyzkZo5EuTCCMnr9bpcLp/PV1ZW5vV6kzXbDG4TI3ZzBBKW2F6TwLFHof9aK8J72X8TlvUJqa6urri4eOPGjUKI8vLyysrKvLwcGfYBkBMpJ0Wyvu9ub2/3eDzaa6fT2dLSktl4AACJyfoRUlNT08yZM7XXDoejs7PT+Ffj0fYRdeRd0N6cRmNz1YhqrFXWJyS3293X16e99vv9hYWF+p+am2W+aRMAIEzWH7Jzu901NTVCiEAg0Nra6nQ6Mx0RACARWT9CKikpKSoqqqqq6urqqq6uznQ4AIAE5chTi4PBoKIoXF8HANkrRxISACDbZd+Qor+/37b8zJkzL730kqkwGAwa64dCoUOHDp08edJUraOjQy8MBoP19fXt7e3JCzlxWjBHjhxJ4I16w7Oo1ZFWbnTxrGVjYyOVpFl3d7f19iLxSGCrtk6mX6r33EglGZG2PTdSSZZybN68OdMxDEEgEKiqqrr11lutfzp9+vTzzz+/aNEivaStre3OO+88fvx4WVmZEOLMmTMrVqz46quvdu3a1d/ff/nll2vVent7/+Zv/uaCCy4oKSlpaGjYsGFDf3//7t2729rarrrqqvS0y9a+ffvWrFnz1Vdfffzxx1u3bp09e/YFF1xgrTN16lRTobHhWdTqKCs3injWsrGx2rusJem3Y8eOf/mXf1m6dOmQ3pXAVm2dTL9U77mRSjIibXtupJLsla0XNZw8efKzzz7TfoH0+9//ftKkSaNGjTLV2bp16/LlyxsbG7XJ3/zmNwsXLly1apUQoqKi4q//+q+1c04/+clPrrvuOq3OlVde+eKLL+qv165dm57mWPX29m7atOn111+fOHGiEKKxsfHee+994403uru7a2trg8Hg7NmzOzo6nnjiiTNnzsyePfvP/uzP9PcaG55drda0trb+/ve/HzVq1LXXXjt69OhAIPCHP/whGAx+8sknTqfzsssuM1aOZy0bG6uxlqTfv//7v1955ZWHDh2aOXOmsY2lpaXf/va3tZLRo0c3NzfffPPN+rsS2KqtkxmUoj1XI0Mz07nnamRodbJk3yE7TUtLywsvvKC93r17t+3Q+NFHHx0/frw+6ff79cmioqKmpiYhxL59+yZNmmT8ZhEKhfr7+59//vmKiooUNiCW/fv3X3XVVdo2LYQoKSkJhUJCiMrKyt7eXkVRenp6vvrqq7y8vFGjRpmu5jA2PLtardm/f//YsWOPHz/+yCOPCCFOnz79t3/7t4cOHfr6179+3333mW7GEXMtWxtrLUm/Q4cOlZSUVFRU7N69Wwy28ciRI2PGjFm7dm1ra+vp06d//OMfv/baa6b+OoGtWob26lK350rSzDTvuZK0OlmydYSUgFmzZq1du9bhcHzyySfHjx/Xbuvw8ssv/+pXv3r77bf1agcOHHj55ZdPnTq1cuXKzAUr/H7/6NGjjSUXXnjhK6+8MmfOnCVLlmglF1988c6dO//yL/8yynyyq9WaO+64449//OOoUaM++ugjreSb3/ym9lWxubk5+q/NTO294IILHn/8cWNjbZuffrt3777ttttmzpz50EMPaecMvvnNb95xxx1CiNOnT7/xxht33nmnEOLBBx+MPp+Y61eS9g5HPNuwPM1M554rT6uTJRcSUiAQiKfa1KlTn3nmmSNHjixZsuS3v/3t9OnTN23adPnll7/55ptHjhzp7e2dPn36jBkzPB6Px+MJhUKLFi265pprCgoKUh2/rXHjxnV3dxtLjh079r3vfW+oV0VmV6uFEC0tLT/+8Y9vvPFG42F3/QYc48ePDwaDUd5uau9rr71mauyvf/1ra/NT2ySL3t7e2tra//3f/xVCOByOPXv2LFy4cNKkSdpfzz///K6uLiHERRddFHNWMdfv7373u4y3N5Ik7rkyrFZNOvdceVqdLNl6yG7ixImnT5/WXp84cSLOdxUWFl5//fUdHR2TJ092OBy33nrr9OnTCwoK8vLy8vPzCwoKtMG1JrMXxF999dVHjhxpbW3VJv/jP/7j4osvPv/88/fv36/XycvL6+npiTmrLGq1EGL37t1er3fVqlWmc0XxM7b3O9/5jqmx1uYnN/54vPnmm8uWLXv66aeffvrpLVu2aEftDh8+rP21qanJ7XbHP7fo67eioiLj7TVK0Z4rw2rVpHPPlafVyZKtIySXyxUMBlevXp2fn3/++efH+a7169f7/f7u7u7HHntMCDF79mz9T729vZdddtm2bdv2798/adKktra25cuXZ3AF5+XlPfnkk/fee29xcbHf7w8Gg4899tjEiRPfe++9lStXjhs3bvHixdo1SOvWrVu5cqV+h1mrLGq1EKKiouLv//7vf/vb344bNy4/P5Ht09he/VC+GGyssaa1JD1ee+21p556Snv953/+52PGjAmFQoWFhatXrx41alRHR8e2bdvivwI++vqdN2+ecTIj7TVK0Z5rrJzZZqZzzzVWlmHlDl+W/TC2paXl8ccf37p1qzYZDAYdDkf8b4/nhg6hUCgYDMrzXcMaszag0UsCgUD0aLOl1frKDYVCqqoOac0aZeNtO/7v//7vwQcf3LJly1A3aZEl7U3Dniub9Oy5OSabRkibNm36wx/+8MADD+glQ91146mfl5cn1RZgjdkUXswskhWtNq7cYUaScCaTQQLBy9/e9Oy5sknPnptjsmyEBOSkQCDw8ccfX3nllZkOBMgkEhIAQAoSHZsCAIxkJCQAgBRISAAAKZCQAABSICEBAKRAQgIASIGEBACQAgkJACAFEhIAQAokJACAFEhIAAApkJAAAFIgIQEApEBCAgBIgYQEAJACCQkAIAUSEgBACiQkAIAUSEgAACmQkAAAUiAhAQCkQEICAEiBhAQAkAIJCQAgBRISAEAKJCQAgBRISAAAKZCQAABS+H/LNsVAZdLZBQAAAABJRU5ErkJggg==" }, "execution_count": 109, "metadata": {}, "output_type": "execute_result" } ], "source": [ "%%\n", "figure\n", "cutoff = 1/7; % (less than weekly )\n", "T = 1/( datenum(date(2)) - datenum( date(1) ) );\n", "[B A] = butter(2, cutoff/(0.3*T), 'low');\n", "\n", "% Aplly the filt forwards and backwards, to try to prevent phase shifting\n", "fx = filtfilt(B,A,flow);\n", "\n", "\n", "tiledlayout(2,1)\n", "nexttile\n", "plot(date,fx, 'b-', 'LineWidth',3); hold on;\n", "plot(date, flow, 'r');\n", "title('Filtered Data');\n", "\n", "nexttile\n", "plot(date,fx, 'b-', 'LineWidth',3); hold on;\n", "plot(date, flow, 'r');\n", "title('Zoom in Filtered Data')\n", "xlim([date(1000) date(1500)]);" ] }, { "cell_type": "markdown", "id": "77162886-bb13-4333-b927-1fdd459b2208", "metadata": {}, "source": [ "## Spectral Analysis - Fourier Space\n", "\n", "Now, we can take out filtered data, and analyze it in Fourier (frequency) space" ] }, { "cell_type": "code", "execution_count": 110, "id": "6ace5ea1-0843-45a1-9334-b5fd91eaafe8", "metadata": {}, "outputs": [ { "data": { "text/html": [ "
dim = 1" ], "text/plain": [ "dim = 1" ] }, "execution_count": 110, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjAAAAGkCAIAAACgjIjwAAAAB3RJTUUH6QkREDoIL35HWgAAFu9JREFUeJzt3X9oVff9+PGjSZo1aoaNWNj2h6VmV4QVW+Jq5iJol5nRT1Ep7RakFt1q1DFWWQvrVgWjdPhr69ptn7KW/YCWtbAyiR+krih1waaIVeeG9q6t2xDthjjqVRKvNt7PH/f7vd98Y0ztvTHnfU4ej7/OzT2R1/Hknqfn3JPruEKhEAFA3MbHPQAARJEgARAIQQIgCIIEQBAECYAgCBIAQRAkAIIgSAAEQZAACIIgARAEQQIgCIIEQBAECYAgCBIAQRCkj5HP53fu3Ll169Z///vfcc8CkGaC9DEeffTRT33qUx0dHTfffHPcswCkWXXcAwTtzTff/MIXvtDa2hr3IADp5wxpOB988MGJEyeeeOKJFStWHD16NO5xANJsTAepq6tr4NnP7t27W1tb586d++yzzxa/8tFHH335y1/+0Y9+tHnz5l//+tcxjQkwJozdIHV2dr799ttnz54tPszlcps3b3711Vf3799//Pjxt956K4qiqVOnnjt3LoqiS5cuVVVVxTkuQNqlPEinTp3q7e0tPcxms6Xl9vb2DRs2lB729PQ0NTXV19dHUbRo0aIdO3ZEUTR//vy//vWv3//+93/4wx92dHSM4uAAY07Kb2qYMmXKc88998gjj9TV1XV3d3/00UeZTKb4VGNj48A1L1y4MHHixOJyXV1dLpcrLm/durW/v9/pEcCNlvIg1dbWrlq16rnnnstkMrW1tfPnzy/jD1EjgFGQ8kt2URTV1tbOmDFj7969d9999zCrTZgwoa+vr7h88eLF0tkSAKMj/UHat29fdXX1xo0bn3/++YHvJw0ye/bsw4cP5/P5KIr27NnT3Nw8ijMCkPYgvfvuu1euXLnnnnuK1+5eeumla63Z0NCwfPnyxYsXL126tLe3d8mSJaM5JwDjCoVC3DME5MqVK4VCwZtGAKMvwUHatm3b73//+9LD5ubmn/zkJzHOA0AlEhykgX7xi1/U1NQ88sgjcQ8CQJnSEKT+/v6vfOUrXV1dkyZNGvj10q8cAVA08PMBQpOG30N6+eWXv/rVrw6qUVHIf/Vly2QytitZ0rpptitxAv9nehrusnvxxReXLVsW9xQAVCTxQXrttdc+//nPf/azn417EAAqkvgg/fKXv/zmN78Z9xQAVCrZQTp48OBNN910xx13xD3IqErr1e20bleU3k2zXYysZAepqanp5ZdfjnsKAEZAsoMEQGoIEgBBECQAgiBIAARBkAAIgiABEARBAiAIggRAEAQJgCAIEgBBECQAgiBIAARBkAAIQpqD9Lf/+u+4RwDgeqU5SAAkiCABEARBAiAIggRAEAQJgCAIEgBBECQAgiBIAARBkAAIgiABEARBAiAIyQ5Sb2/vxo0bFyxY8LWvfe3o0aNxjwNA+arjHqAiTzzxxJw5c9atWxf3IABUKsFnSGfOnDlx4kR7e3vcgwAwAhJ8hnTkyJHPfe5za9euPXLkyMyZMzs7OxsaGgatk8lkigvZbHbUBwSIX+kwGL4EB+nSpUv/+Mc/Xnnllfr6+l/96lfPPPPMhg0bBq2jQ8AYN/AwGHicEnzJ7qabbpo1a1Z9fX0URTNnzvzggw/ingiA8iU4SHfdddeRI0d6e3ujKPr73/9+6623xj0RAOVL8CW7hoaG7373u+3t7Z/5zGc++OCD559/Pu6JAChfgoMURVFbW1tbW1vcUwAwAhJ8yQ6ANBEkAIIgSAAEQZAACIIgARAEQQIgCIIEQBAECYAgCBIAQRAkAIIgSAAEQZAACIIgARAEQQIgCIIEQBAECYAgCBIAQRAkAIIgSAAEQZAACIIgARAEQQIgCIIEQBAECYAgCBIAQRAkAIJQHfcAFbnvvvvOnDlTXN67d29dXV288wBQtmQH6fTp02+//XbcUwAwApIdpJqamuFXyGQyxYVsNnvjxwEITukwGL4EBymfz/f3969cubKqqmrWrFkdHR1Xr6NDwBg38DAYeJwSHKTa2to333yzpqYml8s9+uijNTU1K1asiHsoAMqU7Lvsipfs6uvrH3zwwePHj8c9DgDlS3CQ+vr6SssHDx6cPn16jMMAUKEEX7J75513HnvssUwmc/r06dtuu+3xxx+PeyIAypfgIN1555179uyJewoARkaCL9kBkCaCBEAQBAmAIAgSAEEQJACCIEgABEGQAAiCIAEQBEECIAiCBEAQBAmAIAgSAEEQJACCIEgABEGQAAiCIAEQBEECIAiCBEAQBAmAIAgSAEEQJACCIEgABEGQAAiCIAEQBEECIAhpCNJrr722YMGCuKcAoCKJD9KZM2d+97vfnT9/Pu5BAKhIddwDVOoHP/jBk08++fWvf33IZzOZTHEhm82O4lAAoSgdBsOX7CC9+OKLzc3NjY2N11pBh4AxbuBhMPA4JfiS3fvvv79v374VK1bEPQgAIyDBQfrpT3/a29u7du3atWvX5vP5tWvX9vf3xz0UAGVK8CW7VatWffjhh8XlPXv2PPDAA+PGjYt3JADKluAgzZw5s7RcXV39pS99KcZhAKhQgi/ZDXTo0KG4RwCgIikJEgBJJ0gABEGQAAiCIAEQBEECIAiCBEAQBAmAIAgSAEEQJACCIEgABEGQAAiCIAEQBEECIAiCBEAQBAmAIAgSAEEQJACCIEgABEGQAAiCIAEQBEECIAiCBEAQBAmAIAgSAEEQJACCkOAgnTx5cu3atS0tLXPnzu3s7Ix7HAAqkuAgTZw4cfny5d3d3fv37z9x4sTrr78e90QAlK867gHKN3ny5MmTJ5eW+/v7450HgEokOEhRFJ06daqnp+fgwYNTp05ta2u7eoVMJlNcyGazozsaQBBKh8HwJTtI1dXVN9988/Tp03ft2pXNZq/+e9chYIwbeBgMPE7JDtKtt9567733RlE0ZcqUF154YevWrXFPBECZEnxTw+XLl0vL77333tSpU2McBoAKJfgM6Y033ti+ffvtt99+5syZhoaGp556Ku6JAChfgoPU2tra2toa9xQAjIwEX7IDIE0ECYAgCBIAQRAkAIIgSAAEQZAACIIgARAEQQIgCIIEQBAECYAgCBIAQRAkAIIgSAAEQZAACIIgARAEQQIgCIIEQBAECYAgCBIAQRAkAIIgSAAEQZAACIIgARAEQQIgCIIEQBAECYAgJDhI58+f37RpU0tLS3Nz87Zt2+IeB4CKJDhIhUJh9uzZ3d3d+/fvP3DgQHd3d9wTAVC+6rgHKF99ff3ChQujKBo/fvy0adNyudzV62QymeJCNpsd1eEAwlA6DIYvwUEqyeVyhw8fXrdu3dVP6RAwxg08DAYepwRfsivK5/Nr1qxZv379pEmT4p4FgPIlO0j5fL6jo+Ohhx5qaWmJexYAKpLgS3aXL19evXr10qVLW1tb454FgEolOEhdXV0HDhw4duxY8d2jefPmbdmyJe6hAChTgoN0//3333///XFPAcDISPZ7SACkhiABEARBAiAIggRAEAQJgCAIEgBBECQAgiBIAARBkAAIgiABEARBAiAIggRAEAQJgCAIEgBBECQAgiBIAARBkAAIgiABEARBAiAIggRAEAQJgCCkPEjjvrc37hEAuC4pDxIASSFIAARBkAAIQuKD1NXV1draGvcUAFSqOu4BKtLZ2dnf33/27Nm4BwGgUskOUnt7e2Nj486dO4dZJ5PJRFGUzWZHayiAgBSPgYmQ7CA1NjZ+7DpSBIxlA4+Bgccp8e8hAZAOggRAEAQJgCAIEgBBSEOQDh06FPcIAFQqDUECIAUECYAgpD9I/gcKgERIf5AASARBAiAIggRAEAQJgCAIEgBBECQAgiBIAARBkAAIgiABEARBAiAIggRAEMZKkAZ+op1PtwMI0JgI0rUKpEwA4RgTQQIgfGMoSM6HAEI2hoI00Ljv7R3UJ7kCiNcYDVKJDgEEojruAUbVkPm5niaN+97ewvYFN2AiAP6PsX6GNIxiqIa5uOfsCmAECdL/M/wbS6U+Xb2mMgFUTpAGK8bmk17cG/K0SagArp8gVWTIU6VPmrRBMRs+aR8bORWEYJXx+r3WdZrr+dMSZ1yhUIh7hvLt3r1727Ztvb293/jGN77zne8MejZZe6uwfcHAgYv3UJS+UrqlYtBGlb5r0ApX34LhvgwYNUO+3K71Ur36RT1wheiqo8Hw3zW8TCaTzWaveztGW4KDlMvlFi9evGPHjvr6+jVr1ixbtmzOnDkDV0hWkOJS+lkv/oiLFgxv0IslkONMOoKU4Nu+e3p6mpqa6uvroyhatGjRjh07BgWJ6zHMBcMQlP0vQRLqWtcAghLmiyUFEhykCxcuTJw4sbhcV1eXy+XinYcboYx30UgTu7tymUwm7hGuV4KDxEgZ8vrDkFfwXNMj9a71Y3/1Ozphnr1dbeA1usDjlOAgTZgwoa+vr7h88eLF0tnSmHI9V7GvJyGldQaufPU3qhGpd60f+0FfH/KLV7v6tTlk2yqYN1USHKTZs2c//fTT+Xy+trZ2z549zc3NcU90TQPvhRv476+rvzjwPpzrP0HRCQjT9bxgP/ZlPszxIUEnatcjwXfZRVH0yiuv/OY3v7nlllsaGhqeeeaZQc+Ozk4a8odjRArh+hgQjeihIPC77JIdpOGNZpAAwhd4kBJ8yW5kDfl7qdFQv5pa+roUAYwgQfr/DPl+5pBnWmoEMLLG1mfZDX//zLUaU9i+QH4AbrSxFaSiIW9xBiBeYzFIkRQBhGeMBgmA0IyhIA3/y2gAxGusBGmY/CgTQAjGSpBK5AcgTGMuSACESZAACMKYCJLLdADhGxNBAiB8ggRAENIfJNfrABIh/UECIBEECYAgCBIAQRAkAIIgSAAEQZAACIIgARAEQQIgCIIEQBAECYAgJD5IXV1dra2tcU8BQKWq4x6gIp2dnf39/WfPno17EAAqlewgtbe3NzY27ty5M+5BAKhUsoPU2Nj4setkMpkoirLZ7I0fByA4xWNgIiQsSEePHj158mQURW1tbVVVVdfzLVIEjGUDj4GBxylhQbp06VJfX1/cUwAw8hIWpKampqamprinAGDkJf62bwDSIQ1BOnToUNwjAFCpNAQJgBQQJACCIEgABEGQAAiCIAEQBEECIAiCBEAQBAmAIAgSAEEQJACCIEgABEGQAAiCIAEQBEECIAiCBEAQBAmAIAgSAEEQJACCIEgABEGQAAiCIAEQBEECIAiCBEAQBAmAIAgSAEEQJACCkOAgnT9/ftOmTS0tLc3Nzdu2bYt7HAAqkuAgFQqF2bNnd3d379+//8CBA93d3XFPBED5quMeoHz19fULFy6Momj8+PHTpk3L5XKDVvj8/6z+23/9dyaTiaIom83GMCJA3IrHwERIcJBKcrnc4cOH161bN+SzUgSMZQOPgYHHKWFBOnr06MmTJ6Moamtrq6qqiqIon8+vWbNm/fr1kyZNins6AMqXsCBdunSpr6+v9DCfz3d0dDz00EMtLS0xTgVA5RIWpKampqampuLy5cuXV69evXTp0tbW1ninAqByCQvSQF1dXQcOHDh27Fjx3aN58+Zt2bJl0DqF7QviGA2AT2xcoVCIe4YbJZPJuKMBoCTwo2KCfw8JgDQRJACCIEgABEGQAAiCIAEQBEECIAiCBEAQBAmAIAgSAEEQJACCIEgABEGQAAiCIAEQBEECIAiCBEAQBAmAIAgSAEEQJACCIEgABEGQAAiCIAEQBEECIAiCBEAQBAmAIAhS8mQymbhHuCHSul1RejfNdjGyEhykkydPrl27tqWlZe7cuZ2dnXGPA0BFEhykiRMnLl++vLu7e//+/SdOnHj99dfjngiA8lXHPUD5Jk+ePHny5NJyf39/vPMAUIlxhUIh7hnKd+rUqZ6enoMHD376059+4oknBj3rQjDAINlsNu4RrilhZ0hHjx49efJkFEVtbW1VVVXV1dU333zz9OnTd+3alc1mBxUo5L93AAZJWJAuXbrU19dXenjrrbfee++9URRNmTLlhRde2Lp1a3yjAVCRhAWpqampqampuHz58uWampri8nvvvTd16tT45gKgUgkL0kBvvPHG9u3bb7/99jNnzjQ0NDz11FNxTwRA+ZJ9UwMAqZHg30MCIE0ECYAgCBIAQUhnkHbv3t3a2jp37txnn3027lkqMsyG3HfffXP+r97e3ljGG1ldXV2tra1xTzECrrUhqdll58+f37RpU0tLS3Nz87Zt2+IepyLDb0tqdlliPvmzkDrnzp2bP3/+uXPnCoXC6tWre3p64p6oTMNvyF133RXTXDfEhg0b1q9ff+edd8Y9SKWG2ZDU7LJz58699tprhUKhv7//gQce+NOf/hT3ROUbfltSs8v+85///PnPfy4uP/zww3/84x/jnedaUniG1NPT09TUVF9fH0XRokWLduzYEfdEZRp+Q0q/g5UO7e3tGzZsiHuKETDMhqRml9XX1y9cuDCKovHjx0+bNi2Xy8U9UfmG35bU7LLJkyffcccdpeVgP/kzwb+HdC0XLlyYOHFicbmuri65r5ZhNiSfz/f3969cubKqqmrWrFkdHR0xzThiGhsb4x5hZFxrQ9K3y6IoyuVyhw8fXrduXdyDjICrtyVlu6z0yZ9Tp05ta2uLe5yhpTBIY0Ftbe2bb75ZU1OTy+UeffTRmpqaFStWxD0Uw0nfLsvn82vWrFm/fv2kSZPinqVSQ25LynbZ8J/8GYgUXrKbMGFC6fPuLl68WDrJSJzhN6R4MaG+vv7BBx88fvx4DPPxCaVpl+Xz+Y6OjoceeqilpSXuWSo1zLakaZcVP/nzW9/61rJly1544YW4xxlaCoM0e/bsw4cP5/P5KIr27NnT3Nwc90RlGnJD9u3bF0XRwE+YPXjw4PTp0+Mako+Vvl12+fLl1atXL126tPjuS6INuS2p3GWl5ZA/+TOFl+waGhqWL1++ePHiW265paGhYcmSJXFPVKarN+T9999fuXLlkSNH3nnnncceeyyTyZw+ffq22257/PHH4x6WoaVyl3V1dR04cODYsWPFd1zmzZu3ZcuWuIcq09Xb0tHRkb5dlpRP/kztZ9lduXKlUChUVVXFPUilBm3IwM84v3z5clVV1fjxKTzNTRO7LHHSusvC35bUBgmAZAk3lQCMKYIEQBAECYAgCBIAQRAkAIIgSAAEQZAACIIgAaTcH/7wh7hHuC6CBJBmv/3tb3/2s5/FPcV1ESSA1Hr//fdPnjw5Y8aMuAe5LoIEkE5Xrlz58Y9/nKCPhRUkgETq6upqbW0tPdy9e3dra+vcuXOfffbZ4leefvrpL37xi3/5y18+/PDDQ4cOxTTmJyBIAMnT2dn59ttvnz17tvgwl8tt3rz51Vdf3b9///Hjx996660oihoaGv75z3/u2rXrX//6165du2Kd97qk8P9DAkiHU6dOTZ48ua6urvhw4H893t7e3tjYuHPnzuLDnp6epqam+vr6KIoWLVq0Y8eOOXPmPPzww8Vnv/3tbz/55JOjPv4n5gwJIFBTpkx5/vnne3t7oyjq7u4+ffp06anGxsaBa164cGHixInF5bq6ulwuN/DZn//85zd+2BHgDAkgULW1tatWrXruuecymUxtbe38+fPjnujGcoYEEK7a2toZM2bs3bv37rvvHma1CRMm9PX1FZcvXrxYOltKFkECCNe+ffuqq6s3btxYunY3pNmzZx8+fDifz0dRtGfPnubm5lGcccQIEkCg3n333StXrtxzzz3Fa3cvvfTStdZsaGhYvnz54sWLly5d2tvbu2TJktGcc6SMKxQKcc8AwAi4cuVKoVCoqqqKe5AyCRIAQXDJDoAgCBIAQRAkAIIgSAAEQZAACIIgARCE/wVSbFeBpVj7QAAAAABJRU5ErkJggg==" }, "execution_count": 110, "metadata": {}, "output_type": "execute_result" }, { "name": "stderr", "output_type": "stream", "text": [ "Warning: Imaginary parts of complex X and/or Y arguments ignored." ] } ], "source": [ "% Make date back into a number\n", "d = datenum(date);\n", "d = d - d(1) + 1;\n", "\n", "L=length(date);\n", "n=2^nextpow2(L);\n", "dim=1\n", "\n", "% Take the Fast-Fourier-Transform to put in phase space\n", "orig=fft(fx,n,dim);\n", "\n", "figure\n", "plot(d, orig(1:length(d)));" ] }, { "cell_type": "code", "execution_count": 111, "id": "fc022240-97a9-4438-9ffa-843ae3db2884", "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjAAAAGkCAIAAACgjIjwAAAAB3RJTUUH6QkREDoQPBLfDAAAGY5JREFUeJzt3U9oHOf5B/BXloWprGyIBek1pk4WeihukMCqUKAFkUIxtg82mJAUh2LV7sUG51JqF4dQ6pAeWh/bQw+99GAQPpQ4YGhrbBnjIHqJotCeSnVoSSFrY3sx6v4O22z100rr1f6ZeWb38zkYabU788y7M+9335l3xyO1Wi0BQN525V0AAKQkkAAIQiABEIJAAiAEgQRACAIJgBAEEgAhCCQAQhBIAIQgkAAIQSABEIJAAiAEgQRACAIJgBAEEgAhCCQAQhBIAIQgkAAIQSABEIJAAiAEgQRACAIJgBAEEgAhCCQAQhBIAIQgkAAIQSABEEK4QHrw4MF77703Nzc3MzPzwQcf1B88fPjwoS89evSo/uCNGzfm5+dnZ2evXr2aX70A9Ea4QKrVatPT07du3bp9+/a9e/du3bqVUlpbW7v7pfHx8ZRSpVK5cuXKtWvXbt++vbKycvfu3bwLB6Aru/MuYLNSqfT666+nlHbt2vXSSy9VKpWU0tjY2KanLS0tTU1NlUqllNKRI0cWFxcPHTq06TnlcjmTkgEKY3V1Ne8SthUukBoqlcry8vLFixer1er6+vrp06dHR0cPHjy4sLCQUnr48OHExET9mePj4/XcapZv05fL5SEvIEINuRcQoQYFRKgh9wJS+I/pQQOpWq2ePXv20qVLzz33XErpzp07Y2NjlUrl3LlzY2Njb7/9dt4FAtBj4a4hpZSq1erCwsKbb745NzdXf6R+yq5UKp04cWJlZSWltHfv3sePH9f/+uTJk8ZoCYCCChdIT58+PXPmzBtvvFG/kpRSagRPSun+/fsHDhxIKU1PTy8vL1er1ZTSzZs3Z2ZmcqkWgF4ZqdVqedfw/1y7du2nP/1pY8Tz2muvnTx58sKFC+VyeW1tbf/+/T//+c/37NmTUvr973//29/+dt++fZOTk7/61a+aFxXhjC1AHMF7xXCB1EPBmx4gY8F7xXCn7AAYTgIJgBAEEgAhCCQAQhBIAIQgkAAIQSABEIJAAiAEgQRACAIJgBAEEgAhCCQAQhBIAIQgkAAIQSABEIJAAiAEgQRACAIJgBAEEgAhCCQAQhBIAIQgkAAIQSABEIJAAiAEgQRACAIJgBAEEgAhCCQAQhBIAIQgkAAIQSABEIJAAiCEQQ6kzz5bzbsEANo1yIEEQIEIJABCEEgAhCCQAAhBIAEQgkACIASBBEAIAgmAEAQSACEIJABCEEgAhCCQAAhBIAEQgkACIASBBEAIAgmAEAQSACEIJABCEEgAhCCQAAhBIAEQgkACIASBBEAIAgmAEAQSACEIJABCEEgAhBAukB48ePDee+/Nzc3NzMx88MEH9Qdv3LgxPz8/Ozt79erVxjO3fBCAggoXSLVabXp6+tatW7dv3753796tW7cqlcqVK1euXbt2+/btlZWVu3fvppS2fBCA4tqddwGblUql119/PaW0a9eul156qVKpLC0tTU1NlUqllNKRI0cWFxcPHTq05YPNSyuXy/UfVldXM9wIgCga3WB84QKpoVKpLC8vX7x48aOPPpqYmKg/OD4+XqlUUkoPHz5sfrCZHAKG3MZuMHg4hTtlV1etVs+ePXvp0qXnnnsu71oAyELEQKpWqwsLC2+++ebc3FxKae/evY8fP67/6cmTJ/WB0ZYPAlBc4QLp6dOnZ86ceeONN+pXklJK09PTy8vL1Wo1pXTz5s2ZmZntHgSguMJdQ7p+/fq9e/c++eSTixcvppRee+21999//9SpU0ePHt23b9/k5OSxY8dSSpOTk80PbvLKK+WRkdVaLetNAKADI7XB7bDL5fJnnwkkgP8ql8uRp3qFO2UHwHASSACEIJAACEEgARCCQAIgBIEEQAgCCYAQBBIAIQgkAEIQSACEIJAACEEgARCCQAIgBIEEQAgCCYAQBBIAIQgkAEIQSACEIJAACEEgARCCQAIgBIEEQAgCCYAQBBIAIQgkAEIQSACEIJAACEEgARCCQAIgBIEEQAgCCYAQBBIAIQgkAEIQSACEIJAACEEgARCCQAIgBIEEQAgCCYAQBBIAIQgkAEIQSACEIJAACEEgARCCQAIgBIEEQAgCCYAQBBIAIQgkAEIQSACEIJAACEEgARCCQAIgBIEEQAgCCYAQBBIAIQgkAEIQSACEIJAACCFoIF2/fn1+fr7x6+HDhw996dGjR/UHb9y4MT8/Pzs7e/Xq1ZzKBKBnduddwBbefffd9fX1zz//vPHI2traxx9/vPE5lUrlypUri4uLpVLp7Nmzd+/ePXToUOaVRjEykmq1vIsA6E7EEdLJkycvX7688ZGxsbFNz1laWpqamiqVSimlI0eOLC4uZlcfAH0QcYT08ssvb/y1Wq2ur6+fPn16dHT04MGDCwsLKaWHDx9OTEzUnzA+Pl6pVLZbWrlcTimtrq72s2SAoOp9YCFEDKRN9uzZc+fOnbGxsUqlcu7cubGxsbfffrv9l4siYJht7AODh1PEU3bN6qfsSqXSiRMnVlZWUkp79+59/Phx/a9PnjxpjJYAKKgCBFIjeFJK9+/fP3DgQEppenp6eXm5Wq2mlG7evDkzM5NbfQD0QgFO2X366acXLlwol8tra2v79+9/5513UkqTk5OnTp06evTovn37Jicnjx07lneZAHRlpDa484XL5fJnn60O7vb9j2nfQDvK5XLky+oFOGUHwDAQSACEIJAACEEgARCCQAIgBIEEQAgCCYAQBBIAIQgkAEIQSACEIJAACEEgARCCQAIgBIEEQAgCCYAQBBIAIQgkAEIQSACEIJAACEEgARCCQAIgBIEEQAgCCYAQBBIAIQgkAEIQSACEIJAACEEgARCCQAIgBIEEQAgCCYAQBBIAIQgkAEIQSACEIJAACEEgARCCQAIgBIEEQAgCCYAQBBIAIQgkAEIQSACEMOCBVKulkZG8iwCgDQMeSAAUhUACIASBBEAIAgmAEAQSACEIJABCEEgAhCCQAAhBIAEQgkACIASBBEAIAmlAuGUfUHQCCYAQBBIAIQgkAEIQSACEIJAACCFoIF2/fn1+fr7x640bN+bn52dnZ69evdr6QQAKKmIgvfvuux9//PHnn39e/7VSqVy5cuXatWu3b99eWVm5e/fudg8CUFwRA+nkyZOXL19u/Lq0tDQ1NVUqlVJKR44cWVxc3O7BZuVyufEvwBAqb5B3Lc+wO+8CtvDyyy9v/PXhw4cTExP1n8fHxyuVynYPNltdXR0ZSaurq/2sFyCujR1g8EyKOEICYAgVIJD27t37+PHj+s9PnjypD4y2fBCA4ipAIE1PTy8vL1er1ZTSzZs3Z2ZmtnsQgOKKeA1pk8nJyVOnTh09enTfvn2Tk5PHjh3b7kEAimukVqvlXUO/lMvl+qSGwd3E/6rf6nvgNxPoUr1XzLuKbRXglB0Aw0AgARCCQAIgBIEEQAgCCYAQBBIAIQgkAEIQSACEIJAACEEgARCCQAIgBIEEQAgCCYAQBBIAIQgkAEIQSACEIJAACEEgARCCQAIghOEKpJGRvCsAYBvDFUgAhCWQAAhBIAEQgkACIASBBEAIAgmAEAQSACEIJABCEEgAhCCQAAhBIAEQgkACIASBBEAIAgmAEAQSACEIJABCEEgAhCCQAAhBIAEQgkACIASBBEAIAgmAEAQSACEIJABCEEgAhCCQAAhBIAEQgkACIASBBEAIAgmAEAQSACEIJABCEEgAhCCQAAhBIAEQgkDq3MhI3hW0FLy8jg3qdgECCYAQBj+QajWfqQEKYPADafDIV2AgFSaQDh8+fOhLjx49SinduHFjfn5+dnb26tWreVcHQLd2511Au9bW1j7++OPGr5VK5cqVK4uLi6VS6ezZs3fv3j106FCO5bVjZCTVankXARBVYUZIY2NjG39dWlqampoqlUoppSNHjiwuLra5HOe7AGIqxgipWq2ur6+fPn16dHT04MGDCwsLDx8+nJiYqP91fHy8Uqls+cJyuZxSSmm1+xqMb4Ai+rIbLIBiBNKePXvu3LkzNjZWqVTOnTs3Njb2/PPPt/PC1dXVZFQEDLF6N1gXPJwKdsquVCqdOHFiZWVl7969jx8/rv/pyZMnjdESAAVVjEBqZE9K6f79+wcOHJienl5eXq5WqymlmzdvzszM5FcdAD1QjFN2n3766YULF8rl8tra2v79+9955509e/acOnXq6NGj+/btm5ycPHbsWC6FubCUMWdfYYCN1Aa3Qy2Xy41rSPWtrHdnnW1xc/bsdGm9Sq82KxnIsOzmHQQavWJMxThlNwD699F+IIMHGEICqRUniAAyI5AACEEgsQVDQyB7AmlnnW//ngww5AQSACEIpA51MLfNfxUI0IJAyoKZ2QDPJJB6o3noYzAEsCMCid4TxkAHBBJDTXZCHAJpMHXTz+qjgVwMeyDt6NrPTnvq/k2rkxnA4Bn2QMqFOAFoJpBCE13A8BBInej594oED4BA6rG+/r9HORKZQL8JpP/S4QLkSyA9g6DKTMdNPajv0aBuF2xHIP3PMBz/w7CN7dMaEIpAaks7sxjy7d36eivx1kve9Fe9PNAZgbSFUF1qD7+o2xOhGof2eeOIb4gCyQGZGf/zE9CBoQik9vvHzP7jIv11v2lhqCvQsTAUgdTadnHVTjjVX7vpmQV6+wmuV/uSfZJCEEj5K25n0ai8uJsAxCGQ+qLF6ErfHUqXl7u8m9F4RwptqANpU2z06epRnHDK/lgtYu/wzJqLuFFQCEMdSPnqsl+r51yczjHIt6AGRmbzawbPoO4Sw0Ag/VebB/9OM2DTYgN2MY7efsiyVQf1fr7NotVDzwmkEBqz9Rq2m/jXQz1cWvtBm02foueCumIdCwJpBzIe32zKp36vvSc7brH2/mY9rN8dlXqr5w3oHQlIIA2yvh5y2y28EZzNT2gxBNxuae1sQp+upbU/u0HXBj0hkAqmJ0Ol1lGx00X1fOi2XUfvOn+z9t+45tiWo2EN7VsjkP6fLPu7bqZR9K/Oeqcfav5eGuLjcxPtQAeiHc4tDG8gtf9xu+fP7EmcdLyQbL4HOqhDmU0tEOFQ76yG3MseYPm2baHf2eENpJ5rdMGF64uf+UWfTV8fbnOPbz/Fe7KcblaRsWj19EQGk2IixH8GhmEbtyOQYtm4L/Yj2Oqn49Jw7/Rd6rhbjHkPiPZvhN/Ny9sXqhEcJhkbokAKOHDpxye+FofWpjFcO2tvPKG3pWbfNW/6mtcAG5LN7FK0G510aWA2ZIgCqTNbxlg3d8DL7L9C39HchJ1uRfPVsjbXteXT+vqNq40L7HJFLWa6h7o1Qw9nJ/Zwdj48k0BKqevBU+uXBxyZZSyX0UlzSGR8O4mifwDf0Vye7tfVQxmP5rt8/paG9ib0QxpIffpGS6/mNXRZXo7f18lgvZ0dq+0U1pjyvuWfOtPD7iz3Xib3AprtqKQO3ovmnaHFl7t3tFI3RtnSkAZSKC36ysbFnmzKyH7/jv9d1+2GWa2TY9NGZbaZGbfnlo0weL3kJjvNpO2emUFDxT++NhFIzxb5He2+tvaX0L9h33aDksgt37FuLqdn8KFhpxOvN31m6njOXq9GDG1+YuhgyR2Pdfr0lg1k8A9XIAXv4PK6T8R2Eze6j4TuO9DWBWx3UqVXLVkv/pl9dJfLb/yc7/7Z8aW+FpM28+00M1h7i2MnbXhPW1eSy2WqmIYrkNqXy0X4tOF7Qlv+dWB0MDWxg8+h3WTqM1/YwZ2fet5l9zzDtltai6trHdvp5Z/m2OvmAtLGpTWPWYvS4+/oA18hNkogbWtgMqAfO2LPb0vR14mOm+y0Qfq08DbPL+1o1Zt67Z2ef9u4qNYr6lgjAHqrhyfonvn8jTHWvTZP4bZ4Q9tZQlF6M4FUJO1/Km/zJZnFQFGOhxaaO6BnblS0re7JJah2tjrjD+PtpG/qdZDsVDbNkvuJ3y4NXSAV/dshrfV7d+zJwttZSJ8GXh2PCVovticvzHh2XMbd1sYLVNtNQcylM93p8HdH32Zr/0RC82WnxkKaF/tMxe3lhi6Qoon8cSb72rrskto8Vjv4Uw+1/42oLvVvIl+b/V1nz9nYNXdZSZunLttZxU4Hx+0vuR09SZdCpJRA4hmG6rxcs03TmnOx3ZWDZ+prH9TbAU3r6WodLG3T6GRHpbYzHWbj8jPu6HsyAzYmgbS1sG92c2G9/cSd44aHbfN29KpX2m45/e6AthsEdNn3tWiW1i02AB1uiymLGayloARSRNnM5Y2zwHzl9XmzeR5zbxee16X7Fn/q4bmvHb1l+e6xzWvf+I7vdKTV2UWsohiWQOrHKeDB06dmaf+cSfvzA9M2H+pzP5HVmRYbu2kzW3zuzn2vzvEsVmp77+3sHMNO9XCZ9ZZsPdYcGLvzLoDhtdMPuUX86mLqT+/cnFKtr+23fuaWMwu2/LnNNW564XZvdJvvfjdTG3baWWeT691cG+t6vatdLaLPhmWE1EOD9HlkkLQ/xXZHT+umkgg6KKZXU8h68sy8zj1uqmHTD9msdLuJ8ts9UqBPadsRSFvI/QCgHf2b4hVwB+hrScOzvc+8xNXima1Lav9ccTatvd1aXnmlnMXqOyWQ8hGwC4gv48+nO11j/6YnZCCbgVQ3s/V6vswequ3kf2duvCTFKD4UgUQhbXckR5tP1f0z+/HyPi0qG3nNisxG/MF6XwmkPiqXcx4dNwoIuFtnVlKv3oVuCt7ReZJ2JhS0+WDDxkbIeGeory6bY6H9Rujg5e2svfUgqbMCdjpML7QCB9KNGzfm5+dnZ2evXr2ady1RFH13HAB5TSLoTKgdpijfK2rhmYUN27XAnSpqIFUqlStXrly7du327dsrKyt3795t/7UD8LYNhn5cUchsCX3V7/JyGST17/k9FHnXKlAzdqyogbS0tDQ1NVUqlVJKR44cWVxczLui/MXc/4LP6mEwxNz52amifjH24cOHExMT9Z/Hx8crlcqWT4tzFWdoC4hQQ+4FRKhBARFq2FEBr7yS8q43a0UNpHasrob+TjIAGxX1lN3evXsfP35c//nJkyeN0RIABVXUQJqenl5eXq5WqymlmzdvzszM5F0RAF0p6im7ycnJU6dOHT16dN++fZOTk8eOHcu7IgC6MlIzPQWAAIp6yg6AASOQAAhBIAEQgkACIITBDKS87rt6+PDhQ1969OhRlpVcv359fn6+8euW6+1rMZsKaG6Kvhbw4MGD9957b25ubmZm5oMPPmixuj7VsGUBWTbC3//+9/Pnz8/Nzc3Ozr777rst1tW/d2HLGjLeE+o+/PDD73znOy3WlcGBubGGjBuhzY4o4v2pawPniy+++Pa3v/3FF1/UarUzZ84sLS1ltupXX301l0ouX7586dKlb37zmy3W29diNhVQa2qKfhfwxRdffPjhh7VabX19/fjx43/+858zboTmAmrZNsK///3vv/zlL/Wfv//973/00UfZ7wbNNdQy3xNqtdo///nPt956a2pqart1ZXBgbqyhlnkjtNMR5dhPtjCAI6Qc77s6NjaWSyUnT568fPly6/X2tZhNBaSmpuh3AaVS6fXXX08p7dq166WXXqpUKhk3QnMBKdtGeOGFF77xjW80fl5fX89+N2iuIWW+J6SUfvzjH//kJz+prz37RmiuIWXeCO10RDHvT13UL8a20OZ9V3uuWq2ur6+fPn16dHT04MGDCwsLmVXy8ssvb/x1y/X2tZhNBTQ3Rb8LaKhUKsvLyxcvXvzoo48yboRNBWTfCP/4xz+Wlpbu37//4osvfve737127Vr2LbCphuwb4Xe/+93MzExjh8z+WGiuIeNGaLMjyqufbG0AAykve/bsuXPnztjYWKVSOXfu3NjY2PPPP593Ufloboq33347g/VWq9WzZ89eunTpueeey2B1zywg40bYvXv3V77ylQMHDvzhD3/I687Cm2ool8tZNsLf/va3P/3pT7/+9a/7t4oOasj4cCh0RzSAp+xyvO9qfaRcKpVOnDixsrKSVyVbrjfjYjY1RQYFVKvVhYWFN998c25ubrvV9bWGTQWkzBvhq1/96ve+970f/OAHb7311m9+85tcdoNNNaRsG+GXv/zlo0ePzp8/f/78+Wq1ev78+VdffTXjRmiuYX19PeM9oZ2OKOb9qQcwkPK672rj3U0p3b9//8CBA3lVsuV6syymuSn6XcDTp0/PnDnzxhtv1C/kbLe6/tXQXEDGjfD06dPGz3/9619ffPHF7HeD5hoyboQf/vCHP/rRj44fP378+PHR0dHjx49PTExk3AjNNWTcCG12RDHvTz2Ap+zyuu/qp59+euHChXK5vLa2tn///nfeeWfPnj25VLJlC2TZLM1N0e8Crl+/fu/evU8++eTixYsppddee+3999/PshGaCzh58mSWjfDHP/7xF7/4xde+9rV//etfk5OTP/vZz1544YWMd4PmGjLeE77+9a83ft69e/e3vvWtlFLGjdBcw/LycpaN0GZHFPP+1AN7c9X//Oc/tVptdHQ04/U+ffp0dHR0167/DT3zqmTL9WZZTHNTZFzAdqsb4EZoc/fraws015D7npD7bpCGck/owMAGEgDFMoDXkAAoIoEEQAgCCYAQBBIAIQgkAEIQSACEIJAACEEgARCCQAIgBIEEQAgCCYAQBBIAIQgkAEIQSACEIJAACEEgARCCQAIgBIEEQAgCCYAQBBIAIQgkAEIQSACEIJAACEEgARCCQAIgBIEEQAj/B1+pfPAPOaCbAAAAAElFTkSuQmCC" }, "execution_count": 111, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjAAAAGkCAIAAACgjIjwAAAAB3RJTUUH6QkREDoQPBLfDAAAIABJREFUeJzs3XlYE9f+P/BDAgW0YlFRxOUqNUYQFVyuFFFRULHUBbFXxYDWrUV+KF7LtdZa9XazLm3tgrWtluva1qK4i9SoqIgLaLVSlKoUBaUIvSp7SPL749zOd8xGCCRzEt6vx8dnMjk585kZmA/nzMkZO7VaTQAAAIQmEjoAAAAAQpCQAACAEUhIAADABCQkAABgAhISAAAwAQkJAACYgIQEAABMQEICAAAmICEBAAATkJAAAIAJSEgAAMAEJCQAAGACEhIAADABCQkAAJiAhAQAAExAQoLm5eDBg8XFxY0poKGuru7777/PyspqdGgNoFKpUlJSduzYcfHiRUtuF8CskJD0OnHixNq1a48ePSp0IGCU6urqf/3rX126dGnRooWPj8/nn3+uXSY3N3fcuHEffvihvkrqLaBzu1OnTv3222/5K6dNm+b2tNLSUuPrrNerr74aHh4eExNz+PDhJqwWQFj2QgfAqLy8vDNnzixfvlzoQMBYMpksOTk5Li5u0KBBW7dujYuLq6uri4+Pp++qVCqVStWrV6/MzMw+ffpwn6qrqxOJRCLR//4y0y5AP2hv/9Rvikqlqq6ubtGihb5ghg4d2rZt259++unGjRuxsbGEEEdHR+0K9dWjvdG6ujpCCLdm//79gYGBp0+f5srX1dU988wzhuvhv9TYcQAmqEGXpKSkpKSkVatWXbx4UehYoH55eXmEkPDwcPpSqVR6enp6eHg8efKEELJx40Z3d/eIiAj6MjY2lhabNWsWIaRly5bR0dGEkCdPnnAF6EJcXJyzszMhJCAgoKqqin4qIiJCLBYTQtq1ayeXyzXq5KPV0mWNSLTr4cpobDQ/P9/Ly4v+toaFhanV6iVLlnC/v+np6evXr6flPTw8UlNTtbfFRUiT4tixY+fMmSMWi8Vi8fr16817YgAaovkmpLq6uurqav6arKysO3fu0OW33nrrxx9/VKvVMTExDx8+tHx40CB79+4lhGzdupVbM3XqVEII7ShzdnZeuXJleno6P3mcPHmSEBIdHX3kyJF+/frpTEienp6pqalxcXGEkOTkZFpzRkbGtWvXioqK3N3dJ0yY0KCExEWiXQ9XRmOjmzZtIoScPHny+vXrNG9dvny5VatW/v7+aWlpx48fp3uRmpoaGBjYunXriooKjW3RlxKJ5MiRIxMmTCCETJ06NTU1VSKRtG7d2sxnBqABmmmXXX5+/sKFCwcMGPD222/TNTKZTCqVFhQUDB8+XCaT9ejRo2PHjoQQDw+PsrKytm3bChov1EOlUmmsoZ1RdP306dNXrFhBCCkvL+cK/Pbbb4SQqKiokJCQ27dv0441DSEhIaNHjyaEfPbZZxUVFdwHjx07VlFRoVKpuJVG4iIxUI/GRv38/BwdHWmzZtmyZYQQX19fR0fHtm3bhoSEbN68mRAyc+bMESNGlJSUyGSy9PT0wMBA7b0OCQkJDQ0tKSnZt29fVFTU6NGjBwwYQBM5ACOaaQ9yYmJiVFQU9zIjI6Nz587Lly//+uuvd+zYoVKpwsPDd+/evWLFisePH0skEgFDBWN4e3sTQuRyObfm8uXLrVu3dnJyIrz7N9ocHBwIIS1btjTwLv9Gy44dO6Kjo11dXefNm9elS5eGxslFYqAejY0OGjQoJydn5syZiYmJ/v7+OqtVKpXcMpebNfaavxfad5sAWGDjLaRLly7Z29v7+voSQlJTU729vekv/5o1a9LS0rhiRUVF3K+6RCLJy8uTSqUff/yxSqXCXV+r0KtXr2HDhm3btq1Xr169e/fet2/fr7/++tZbbxn4iKenJyFk69athJAvv/zSyA3RdlVERISTk1NOTs6QIUNMC9j4erKyssrLy9etW1dQUHDs2DGNdwcPHkwI2bZtm0qlSkxMbNmyZVBQkHZ7EcAq2PjVduDAgQ8ePLhy5UpqamqHDh30/Umbm5tL/ywlhIjFYm6ELrKRFfnxxx/Hjx+/bNmycePGbdu2bfHixe+8846B8iNGjIiMjNyyZUt4eDhtSBlj+vTpHh4eQUFBr7zyypQpU0yO1vh6bt68GRwc3LJly2PHjq1bt07jXR8fn40bN+7du3fMmDG3b9/+4YcfDIz9A2CcnVqtFjoGs/voo486duw4bdo0/sq0tLRz587Re0jJycmEkIiICEJIbGxsfHw8uumslL4B0PrQ0c9z587dsmWLQqHQGN6tT21tbZP0ehlZD90pe3t7A38hNVVIAAKy/RbA0aNHR44c2a5duytXrugr4+XlRW8/KBSKW7duIRtZL5FIZOR1uby8fPDgwbNnz549e/auXbsmTJhgZDYiTXcPxsh66E4Zbq8jG4ENsPF7SNnZ2e7u7vQeUlpa2t27d3X22nl7e7u7u8+bN6+srIz7KiXYNicnp2nTpv3888+EkE8++WTmzJlCRwTQ3DWLLjsjKZVKOzs73DcCABAEEhIAADABrQEAAGCCjdxDUiqVFy9ebNu2LX88gt1iec+DMQJGBQDAmhs3bggdgl62kJDS0tI2b948ZMgQDw8PjQFyLB96A6RSKSK3MERueYjc8qRSqdAhGGL1Camqqurzzz/fs2cPnTgZAACslNUnpDNnzgwdOlQul1dUVAwZMsTNzY3/LvfngJX+OQMA0EiMt4r4rD4hlZeXX7hwwc/Pz9nZedasWSkpKfymkpXmISsNmyByISByy7OuyPnRMp6crH6UnZ2d3fjx44ODg8eMGdOrV6+LFy8KHREAAJjC6hOSl5dXRkYGXS4qKuratauw8QAAgGmsvstOKpV6eXnNnTu3trbW39/fw8ND6IgAAMAUVp+QCCFxcXH0ATCY9QcAwHrZQkIiSEUAANYP13EAAGACEhIAADABCQkAAJiAhAQAAExAQgIAACYgIQEAABOQkAAAgAm2k5CKi4tLSkqEjgIAAExkIwmpqqpq5syZp06dEjoQAAAwkY0kpPfff3/kyJFCRwEAAKazhamDTpw40a5dux49elRVVWm8hQf0AUAzx/gzkPisPiGVlpbu3Llz06ZNR44c0X4XeQgAmjkrekCf1Sekjz/+2MfHZ//+/VevXq2qqurRo4evr6/QQQEAQINZfUKaOHFicXExIUQkEtnb2zs4OAgdEQAAmMLqE9LAgQO55aqqqt69ewsYDAAAmMzqExInLCxM6BAAAMB0NjLsGwAArB0SEgAAMAEJCQAAmICEBAAATEBCAgAAJiAhAQAAE5CQAACACUhIAADABFtISCqV6vz583QCIQAAsFJWP1NDXl7eypUre/To8dtvvwUHB8+aNUvoiAAAwBRWn5AkEsm2bdtEIpFKpQoKCkJCAgCwUlafkAghIpGIEFJYWOjp6anxFh7QBwDNHOPPQOKzhYRECKmpqVmxYsXSpUs11iMPAUAzZ0UP6LOFQQ0KhWLRokVz5sxh/FgDAIABVp+QlEplfHy8TCYLCAgQOhYAADCd1SekDRs25OTkJCUlzZs3b968eTk5OUJHBAAAprBTq9VCx2Audovl6vUjhY4CAIAVUqmU5TvrVt9CAgAA22DjCclusVzoEAAAwCg2npAAAMBaICEBAAATkJAAAIAJtp+QcBsJAMAq2H5CAgAAq4CEBAAATLCRhJSdnZ2fn6/vXbvFcn7HHbeM3jwAAHbYQkKSyWSHDh167733tm/fbqAYTT/c/xovLRIpAADoZfVTB2VkZOzfv3/16tWEkLFjxx46dIg+HolYJM2o14/ktkKnKeLmK2r8xEWY+gig2TLTrz/jUwdZfUL68ccf7e3tJ06cSAhZsGBBbGws9xAKtHsAwLYZk7Q0nsvDckKy+gf05ebm+vn50WWxWFxaWipsPCAIflOVWyNUMMChf+bz/yd//aWofcrATKzoAX1Wn5C8vLyqq6vpcm1trZubm4UD4C58TdhZBybAMWcQPSn8/3UuAFC2kJC++OKLiIgIhUJx69YtiURioLC+P8qa5BcDv2YAAI1h9QnJ29vb3d193rx5ZWVl8fHx+orxkwQSBgAAg6w+IRFCli9frlQq7ezsuPF1BiAbAQCwyRYSEiFELBbrewsZCADAKtjCF2MBAMAGICEBAAATkJAAAIAJSEgAAMAEJCQAAGACEhIAADABCQkAAJhg9QlJpVKdP3++uLhY6EAAAKBRrDsh5eXlRUVFHT58+J///OeWLVuEDqfJMD4jrwGI3PIQueVZb+SMs+6ZGiQSybZt20QikUqlCgoKmjVrltARAQCAiay7hUQIofPXFRYWenp6Ch0LAACYzipbSEqlUqVSiUQiOoVdTU3NihUrli5dqlGs58GYmy9tpI1rlh+SCABgPlbUwWiVjzDfsWNHRkZGr1694uLiFArFwoULZTJZQECARjGpVHrzpY09D8YIEiQAAINY/uvcKhMSR6lULliwYPr06drZCAAArIt1J6SPPvrowIED3FNi4+Pjvb29hQ0JAABMY90JCQAAbIbVj7IDAADbgIQEAABMEK9cuVLoGMwiOzu7pqbmueeeEzoQQ/QFqbG+srJSoVAoFAqVSmVvz9ZIfX27oFQqFQoFa9ESowO2umOuUqkuXLjg4ODw7LPPChWYPkYGbHXHXKlUXrhwQSQStWrVSqjA9DEyYNaOuW0mJJlMVlZWdvDgwcePH/ft21focHTTF6T2+tGjR+fk5Mjl8ocPH/r6+goXsiZ9u5Cfn//KK6/cvn17+PDhAoanzfiAreuY5+XlLVy4sLS0dNeuXY8ePfLz8xM2SD7jA7auY56dnb1kyZKamprk5OT8/PzBgwcLGySf8QEzd8zVNufs2bNLliyhy6GhoUqlUth4dNIXpPb6urq61157TZgoDTJwnBMSEnbv3r1q1SqBQtPN+ICt8ZjTZaVSOXToUGGC08X4gK3xmHP8/PwsG5QhxgfM4DG3wXtIRUVF/v7+dFkikeTl5Qkbj076gtReX1ZW9ujRo6NHj2ZkZAgTqx4GjvOaNWtat24tUFx6GR+wNR5zNufQMj5gazzmKpWqpqbmP//5z7hx4wSKTgfjA2bwmNtgQsrNzXVwcKDLYrG4tLRU2Hh00hek9np7e/tJkyYplcr09HSmZo+1iuPMZ3zAVnrM9c2hJSDjA7bGY37hwoWEhAS5XB4UFCRMcLoYHzCDx1z4u1hNzsvLq7q6mi7X1ta6ubkJG49O+oLUXu/q6jp58mRCSFhYWGRkZGVlZYsWLQSJWYNVHGc+4wO2xmOuUCgWLVo0Z84cpiYuMz5gazzm/v7+/v7+KpUqLCwsMDCQSwPCMj5gBo+5DbaQvLy85HI5IUShUNy6dYubx4Ep2kHevXtX53ruIwqFoqyszNHRUaiYNejbBWaZELC1HHOlUhkfH69zRkdhmRCwtRxzlUrFlVGzNL2ACQGzc8xtsIXk7e3t7u4+b968srKy+Ph4ocPRTSPIzMxM+r928BkZGV988UX79u3v378fGxtLJzhngb5dEDouvYwP2OqO+YYNG3JycpKSkpKSkghLc2gZH/B///tf6zrmW7ZsOXPmTLt27fLz86OiohhpHpGGBMzgz7nNTh2kVCrt7OzojVNm8YOsqanh/kLRDl6hULDzE8+nbxeYZXzAOOZNxVaPuUqlUiqV1h4wU8fcZhMSAABYF6YbEAAA0HwgIQEAABOQkAAAgAlISAAAwAQkJAAAYAISEgAAMAEJCQAAmICEBAAATEBCAgAAJiAhAQAAE5CQAACACUhIAADABCQkAABgAhISAAAwAQkJAACYgIQEAABMQEICAAAmICEBAAATkJAAAIAJSEgAAMAEJCQAAGACEhIAADABCQkAAJiAhAQAAExAQgIAACYgIQEAABOQkAAAgAlISAAAwAQkJAAAYAISEgAAMAEJCQAAmICEBM3dwYMHi4uLG/qWgFQqVUpKyo4dOy5evCh0LABNCQkJbISbm5ubm1uHDh3Gjx+/Z88eIz+Vm5s7bty4Dz/8kL4sLCy8cuWKzreMN23aNLenlZaWNrQSA1599dXw8PCYmJjDhw83YbUAgrNTq9VCxwDQBOzs7CQSSUhIiFwuv3HjxrVr13x8fAghKpWqrq7umWee4Reuq6sjhNjb2xNCzp8/36dPnxYtWhBCZsyYUVtbu2vXLlqM/5Y2nTUTQhITE3Nycn766acbN27ExsYSQlavXv3ss8/Sj6hUKrpdlUpVXV2tXTm/jHa0hJAOHTr07Nnz9OnThsPQqIf/sq6uTiQSiUT4exQYowawCYSQqVOnqtXqffv2EUJ27typVqvXr1/v7OxMCPHw8EhNTVWr1fn5+V5eXvSHPyws7MmTJ4SQ2NhYtVq9cuVK7vfiu+++4976+9//3q1bN7qVWbNmOTg4PHz4cO3atQ4ODoSQrl27ZmRkaMcTHR3N/X7RqjZu3Oju7h4REaFWqyMiIsRiMSGkXbt2crmcKxMXF0cDDggIqKqq0ohWrVYvWbKECzI9PV17BzW2xe2Fo6MjIWTs2LFz5swRi8VisXj9+vXmPy0ADYCEBDaCEBIeHn7p0iWaCc6cOZOZmUkIiY6OTk1NDQwMbN26dUVFxaZNmwghJ0+evH79ulwu5yeky5cvOzg4BAUFpaWlPXjwgHtr48aNhJALFy4olcrWrVvTrRBC3n777YqKisDAQC8vL+14tBOSs7PzypUr09PT1Wp1RkbGtWvXioqK3N3dJ0yYwJXx9PRMTU2Ni4sjhCQnJ2tES4Ns1aqVv79/Wlra8ePHtXdQY1v0pUQiOXLkyIQJE2jaTk1NlUgkrVu3ttzpATACEhLYCH67nyaYb775hhBCr+Pbt28nhBw5cuTChQuOjo7Ozs5xcXH8rEMrcXR0pM0s9V8ZIjY29smTJw4ODosWLTpy5AghJDU1ldYcFBQ0depUiUSis6dBOyHNmTOHe3fr1q0ymSw8PLx9+/YhISFcmXnz5qnV6tTUVELI1q1bNaKln23Xrh1tLencQY1t0ZcxMTF0o4SQQ4cOqdXqqVOnOjo6NvE5AGgcdCKD7QgKCpLL5VVVVZ9//jm3UqlUcssqlWrQoEE5OTkzZ85MTEz09/c3ptpnn312ypQpe/fuTU5O9vDwGD16NF0/bNiwF198cfny5fRCXy/aaUYI2bFjR3R0tKur67x587p06cIvQ7sBubs7xkSrsYMa26L4t4u07zYBMMK+/iIAVsLd3X3EiBHcy8GDBxNCtm3bplKpEhMTW7ZsGRQUlJWVVV5evm7duoKCgmPHjmnUIBKJ8vLyjh075uXl5erqyq2fMWPG9u3bjx07NnPmTEKIr68vIeTevXvz58+/f//+gwcPGhTnb7/9RgiJiIhwcnLKyckZMmSIvpKGo9W5g1xOArA6aCGBzfLx8dm4cePevXvHjBlz+/btH374oUWLFjdv3gwODm7ZsuWxY8fWrVun8ZHFixdnZWWNGTPmzJkz/PUhISFdu3YtKCiYMWMGIWTAgAG0Znd3dz8/P7lc3qDApk+f7uHhERQU9Morr0yZMsVAScPR6tzBBkUCwBQM+wbbV1tby++nouOk7e3tdY571hhjXW/N+uppaFT6GI62QVUBMA4JCQAAmIAuOwAAYAISEgAAMAEJCQAAmICEBAAATLDl7yHZLZb3PBgjdBQAAAy5ceOG0CHoZcsJibB96E0mlUqxX9bFVncN+2V1pFKp0CEYgi47AABgAusJSaVSnT9/XuOpndnZ2fn5+YbXAACAdWG6yy4vL2/lypU9evT47bffgoODZ82aRQiRyWRSqbSgoGD48OEymUznGttmq50JtrpfxHZ3DfsFTYvphCSRSLZt2yYSiVQqVVBQ0KxZszIyMjp37rx8+XJCyNixYyMjIzMzMzXW4DmYAADWiPVrN80uhYWFnp6ehJCioiJuEn6JRJKXl6e9RqhQAQCgMZhuIVE1NTUrVqxYunQpISQ3N9fPz4+uF4vFpaWl2mv4n+WGlKANDgDNE+Mj6/hYT0gKhWLRokVz5syhx9TLy6u6upq+VVtb6+bmpr2G/3HkIQBo5viXQcaTE9NddkqlMj4+XiaTBQQE0DVeXl702TMKheLWrVsSiUR7jZARAwCAqZhuIW3YsCEnJycpKSkpKYkQEh8f7+3t7e7uPm/evLKysvj4eEKI9hoAALBGVvk8JKVSaWdnxx9Np72GEGK3WK5eP9Li0QEAMIrxSSiYbiHpIxaL610DAADWhel7SAAA0HwgIQEAABOQkAAAgAlISAAAwAQkJAAAYAISEgAAMAEJCQAAmICEBAAATEBCAgAAJiAhAQAAE5CQAACACUhIAADABCQkAABgAhISAAAwAQkJAACYgIQEAABMQEICAAAmICEBAAATkJAAAIAJSEgAAMAEJCQAAGACEhIAADABCQkAAJhgBQlJqVTW1NTw12RnZ+fn5xteAwAA1sVe6ADqkZ+fv3DhwgEDBrz99tt0jUwmk0qlBQUFw4cPl8lkOtcAAIDVYT0hJSYmRkVF5eTk0JcZGRmdO3devnw5IWTs2LGRkZGZmZkaa0QiK2j2AQCABtav3WvWrGndujX3sqioyN/fny5LJJK8vDztNQJECQAAjcZ6C0lDbm6un58fXRaLxaWlpdpr+OWlUilduHHjhiXjBABgBHcZZJ+VJSQvL6/q6mq6XFtb6+bmpr2GXx55CACaOf5lkPHkxHqXnQYvLy+5XE4IUSgUt27dkkgk2muEjhEAAExhZS0kb29vd3f3efPmlZWVxcfH61wDAADWyE6tVgsdQ4MplUo7Ozv+aDrtNYQQu8Vy9fqRFo8OAIBRUqmU5RsZVtZCosRicb1rAADAuljZPSQAALBVSEgAAMAEG09IdovlQocAAABGsfGEBAAA1gIJCQAAmICEBAAATEBCAgAAJiAhAQAAE5CQAACACUhIAADABCQkAABgAhISAAAwAQkJAACYgIQEAABMQEICAAAmICEBAAATkJAAAIAJtp+Q8AQKAACrYPsJCQAArAISEgAAMKFZJCT02gEAsK9ZJCQAAGBfc0lIdovltJ3ELRC0nAAAWGIjCSk7Ozs/P7/eYvxUpJGfrCg5SaVSoUMwC1vdL2K7u4b9gqZlL3QATUAmk0ml0oKCguHDh8tkMhNqaGhOUq8fyRWmy+r1I03YLgAAcKw+IWVkZHTu3Hn58uWEkLFjx0ZGRopEZm/28VOXyQ0sflYz8JKmOv5bPeuLDdkRgE31/npqXAS4NRrXB2228Vtvp1arhY6hUX788Ud7e/uJEycSQhYsWBAbG8s1t62oFw4AwGTGZyOpVHrjxg2zBtMYVt9Cys3N9fPzo8tisbi0tJR7q+fBGIGCAgCwHOlBoSNoIlafkLy8vKqrq+lybW2tm5sb9xbLfwgAAIAGqx9l5+XlJZfLCSEKheLWrVsSiUToiAAAwBRW30Ly9vZ2d3efN29eWVlZfHy80OEAAICJrH5QA6VUKu3s7Cwwvg4AAMzERhISAABYO1a67CorK7/55psDBw7QEQotW7aMiIiYOXOmg4ODaRVmZ2e3adOmW7duTRmlQHTuS2VlJV0Qi8WOjo4ChNXUlEplXV2dDe+LLZ0ylUp18eLFbt26dejQQehYmoC+3bGlU6ZUKi9evNi1a1cPDw+hY9FLvHLlSqFjIJ999tnChQsHDhz4xhtvREdHT5s2LTQ0NDMzMy4urk2bNj4+Pg2tUCaTlZWVHTx48PHjx3379jVHzBajb19Gjx6dk5Mjl8sfPnzo6+srYIRNIj8//5VXXrl9+/bw4cOFjqWx9O2LzZyyvLy8hQsXlpaW7tq169GjR9z3LqyUgd2xmVOWnZ29ZMmSmpqa5OTk/Pz8wYMHCx2RHmoG7Nu3T99be/fubWhtZ8+eXbJkCV0ODQ1VKpWmRyY0fftSV1f32muvCRdX00tISNi9e/eqVauEDqQJ6NwXGztl9EdRqVQOHTpU6FiagM7dsbFTxvHz8xM6BL2YGAUwfvx4fW/RKRgapKioyN/fny5LJJK8vDzTIxOavn0pKyt79OjR0aNHMzIyhIuuKa1Zs6Z169ZCR9E0dO6LjZ0yOoCosLDQ09NT6FiagM7dsbFTplKpampq/vOf/4wbN07oWPRiIiERQiorKz/99NNRo0YNHTp06NChoaGhX3/9tUKhMKGq3Nxc7s6TxtwNVkffvtjb20+aNEmpVKanp8+aNUu4AMFYtnfKampqVqxYsXTpUqEDaRrau2Njp+zChQsJCQlyuTwoKEjoWPRiYlDDZ599lpSUNHv27KSkpBYtWhBCysvLt2/fPmDAgGXLlk2ZMqVBtRmYu8Hq6NsXV1fXyZMnE0LCwsIiIyMrKyvpcQNm2dgpUygUixYtmjNnjm08qUHn7tjYKfP39/f391epVGFhYYGBgSaPFzMrJlpIf/vb37KysubPn9+pUydXV1dXV9cuXbosXbr06tWrJoxssaW5G7T35e7du/wCCoWirKzM2sf/2DbbO2VKpTI+Pl4mkwUEBAgdSxPQ3h3bO2UqlYpbVjP8VR8mWkjcPaTk5ORt27YlJiZ6eHjk5OR4e3ubcA/JluZu0NiXzMxM+n9GRsYXX3zRvn37+/fvx8bGisVioSMF3WzylG3YsCEnJycpKSkpKYkQEh8f7+3tLXRQptPYnenTpy9ZssTGTtmWLVvOnDnTrl27/Pz8qKgoNptHhKkvxsbExDz33HO///77unXrPDw8hgwZcvbsWZNrs6W5G/j7UlNTw/2lplAomP3BAg5OmdWxyVOmUqmUSiXj+8LQ9fr69esffPAB99LJyamqqsrk2sRisW1kI/L0vvD7DRj/2QIKp8zq2OQpE4lE7O8LQ5fsiooK/ss///zT2dlZqGAAAMDCmLiHRK1evdrX17eqqiomJubWrVv0qeQAANBMMHQPiRBSU1NTUFBACOnUqZNVj7AEAICGYishAQBAsyV8l11eXt7LL7+svV4kEmVnZ1s+HgAAEAQTLST7mhz5AAAgAElEQVSlUqlzvVUP/AcAYMSpU6fy8/NnzJghdCD1EL6FRJ5OPI8fPyaEuLi4CBcOAIDtKC0tPX36NDcJGcsYGvZ96dIlX1/fESNGjBgxwsfHJzU1VeiIAACsXmJiYkxMjNBRGIWJFhL16quvpqSk0OeilpaWDhs27Pr160IHBQDANO3HE/OfMZ2cnDx27Fhr6XNiqIXUsmVL7indbdu27dKlS01NjaARAQAwLT8/f9KkSR9++CG3RiaTHTp06L333tu+fTshpKio6IcffnjzzTfPnz/P/jAxhlpIPj4+hYWFnTp1IoSUlpY+++yzVj29LgCAuSUmJkZFReXk5NCXGRkZnTt3prMKjB07NjIyMi4ujhCiUChWrVrVv39/IWM1AkMJ6ffffx85ciRNQjU1NS4uLvQh9r169fruu++Ejg4AQDCXLl2yt7enl8TU1FRvb+8uXboQQtasWZOWlsYV037GNH3Ck4ODw7vvvitE4A3DUELav3+/0CEAALBo4MCBR48evXLlSnFxcYcOHWg20pabm+vn50eXrfF52QwlJLFYfOvWLf59I6t+yAoAQBMKDQ396KOPOnbsOGbMGH1lrP152QwlpNDQUDs7O1dXV/rSzs5ux44dwoYEAMCIo0ePjhw5sqSk5MqVK7TvTpuXl9cXX3wRERFhpc/LZighPXr06Ny5c0JHAQDAnOzsbHd3d5qH0tLS7t69q7PXztqfl83E1EHUyJEj5XK50FEAAFg3631eNkMJqbS0dNSoUZ06daLH0c7OLiUlReigAADAQhjqsgsNDX3zzTcHDhwodCAAACAAhhKSs7Pz5MmThY4CAACEwVAno1QqvXv3rtBRAACAMBhqIRUVFYWEhDg6OtJ7SHhAHwBAs8LQoAbtx/ThAX0AAM0HQwmJEFJYWKhQKLiX3OTfAABg8xjqspPJZPfv3793716PHj1KS0vDwsLonLUAANAcMDSooaCg4Pjx44MGDdq1a1dmZuaxY8eEjggAACyHoYREOw9ff/31J0+eEEJcXFyqqqqEDgoAACyEoS47iURSXFzs6ekZGBjYp0+fwsJCZ2dnoYMCAAALYWtQA1VYWHj//n2pVNqqVSuhYwEAAAthMSEBAEAzxNA9pE8//ZQufPbZZyNGjDh+/Liw8QAAgCUxlJC+//57QkhJScmOHTv27NmzbNkyoSMCU9TV1X3//fdZWVkNeqt5stgBOXjwYHFxsW1sBWwYQwnpmWeeIYSsWbNm48aNrq6ubdu2xSg7wfn4+Lzyyit0+Y033nBzc7t69SohpLy83M3N7Z133tH+SHV19dSpU7/99tsGvcVxc3Nzc3Pr0KHD+PHj9+zZ0xQ7UY9169Z16dKlRYsWQ4YM+emnnxpTVWFh4ZUrV4wvr/OA0CPQsWPHl1566ZNPPmlMPFRubu64ceM+/PDDBkXoxvPGG2/oLMOvjb8V0zT06IHtYSghOTo6Llq06Pjx435+foSQx48fY5Sd4Pr160dbroSQn3766eHDhydPniSEyOXyhw8f0jNFCFGpVLW1tQbqqa6uVqlUdNnNzU2lUtXV1eks+fDhQ1dX14iIiJs3b0ZERPzyyy9GbkKDvvIamz59+nRCQsLgwYO/+OKLjh07tmvXTl9JfXXyd+3NN9/UviLz61GpVJWVlYYj546Avb39okWLJk6c2NCdIoTU1dVxa3r16pWZmfnuu+/qi9BAGFOmTJkyZcqAAQO0q9Wojb8V7ajq6uq4o8S9q3EoNGJr6BkHW6Bmyc8//1xbW8stCxsMqNXqrVu3EkKuXbtGvxzWrVu38PBwtVq9ZMkSsVisUCjUavXatWsdHBwIIV27ds3IyKAlY2NjaQ3r169v2bIlIcTZ2fnXX38lhERGRtLxk4GBgVVVVRpbJIRMnTpVrVbv27ePELJz505aCf3rxMPDIzU1Va1Wh4SESCQStVq9d+9eQsiYMWPUavXGjRsJIb/++qtGSGq1mka1ceNGd3f3iIgIbnM7d+4khGzdupVbQ0vGxMTQsOfMmUPXa9fJ37WCgoKVK1dyv1bfffed9hYjIiLo9Izt2rWTy+XctrhjpXEE6HEmhKSnp+sMgH48Li6OHpyAgICqqqr8/HwvLy8aRlhYGH8rGhH6+/u7uroqlUq1Wh0ZGeng4PDw4UOdYajVau1qde5vbGwst+Do6EgIGTt27Jw5c8RisVgsXr9+Pa1N+1Bo1Ka9s9AcMJGQQkNDv/7664qKCqEDAU35+fmEkG+++SY5OVksFm/ZssXZ2VmpVA4bNmzYsGFqtfrSpUuEkLfffruioiIwMNDLy4t/kc3MzKTXtdTU1C1bttC3PD09U1NT4+LiCCHJyckaWySEhIeHX7p0KTo6mhBy5swZWkl0dHRqampgYGDr1q0rKirWr19PCLl3797ChQs9PT0dHR2VSmVERISnp6d2SOq/LsrOzs4rV66k13fq4cOH7u7uhJBZs2bdu3ePK9mvX7/U1NSZM2fSGLTr1Ng1tVp9+fJlBweHoKCgtLS0Bw8eaG8xIyPj2rVrRUVF7u7uEyZMUBuRkI4cOUKPv4Gd0jiemzZtIoScPHny+vXrGmlPI8JvvvmGEJKamqpQKJydnWlI/DA8PDzCwsLCwsIyMjK0q9W5v1xCkkgkR44cmTBhAneUJBJJ69ataeXah4JfG52lRWNnoTlg4ouxe/bsSU5ODg4ObtGixfvvvz948GChI4L/+dvf/ta1a1e5XP7ss88OHjz4hRdeqKqqOnv27Llz5+gf77TTPz09/ebNm8XFxXl5efyP0w632bNnh4SEEELKy8sJISEhIaNHjyaEfPbZZxUVFdob3bt3L233xMbGDhkyZPPmzYSQmTNnjhgxoqSkRCaTpaenv/jii4sXLz558uTJkydXrlw5e/bss2fPHj58eM6cOQZCmj59+ooVK/jbatu27ZUrV95+++3Nmzfv3r07IyODTukbEBBAg0xKSrp58ybtbuLXqbFrhBBfX1+RSOTu7s7fWf4Wf/vtt2PHjlVUVKhUKp07ro9SqTSwUxrH08/Pz9HRkbZLNEYGaUQYFRW1cOHC3bt3V1dXV1VVzZs3T2O7Tk5OHTt2JIQ4OztrV6tzf/lRhYaGlpSU7Nu3LyoqavTo0QMGDKCnVeeh4NdGz7i+HyqwYUwkJGdnZ5lMJpPJiouL165dO3v27IEDB7733nudOnUSOjQgQUFB586dc3V1DQkJ6dWrl7u7+1dffaVQKIYOHcqVGTZsWI8ePV588UWNz2o/UoQQQrti6FOv9G3x7bfffuGFF5ycnHRWpVKpevXq5eHhIZfLf/nll8GDB4eEhGzdurWqqio0NPT+/fv6QqKdSBo6dOiwadOmmJgYPz+/zz77jLa9+EQiEU1I/DqNHHHDbXHHjh3R0dFxcXHTp08vKCgw5rOEkHPnzhFCPD09f//9d307pXE8Bw0alJOTs27dusTExAMHDty5c0df5c8888ycOXO2bt1aXV3dvn370NBQjQJ///vfv/76a+6lkdVS/PNLxytxjDwU+n6owIYxkZA4HTp0WLdu3bp1627cuDF79uyjR48KHRGQF198cceOHSKRiI6pCwkJ+f777x0cHOjfxb6+voSQe/fuzZ8///79+w8ePOB/ViqVEkK+/PJLQsjNmzf/8Y9/GLNFd3f3ESNGcC9pi3nbtm0qlSoxMbFly5ZBQUE0sGPHjrm4uPTs2bNv376JiYkODg6hoaGXL182EJKGlJSUy5cvDxo0qKioiBDi4eFB1//0009Hjx7dtm0bIaRnz570ksqvk9624XZt1qxZTk5OIpEoLy/v2LFjXl5erq6u/A399ttvhJCIiAgnJ6ecnJwhQ4YYiKqgoGDXrl3nzp377LPPhg0bFhISQoeGG7NTWVlZ5eXl69atKygo0J6hmB9hly5dZs2atWHDhmPHjkVFRWn/iVBQUECHtPTo0YMQol2tgf01QN+h4GqjzTIjzyDYFKH7DIF19Ert4OBA735/9913hJCRI0dyBegwffrjlJCQoHFfJCEhgd6+bt++fUlJCfdWWloaeXo0AUWevpfObYKOg3B3dz906BBdSft/Xn75ZW553LhxOkNS67lbo1arDx06RFOLWCwODw9/8uQJLTlhwgQ6YIH7iHad/F2jozPeeustWmDnzp0aW8zLy6PZzsvLa+bMmSEhIfqiojWIxWKJRLJkyRLu3qrhneKO586dO2lUDg4OGzZs0CjGj5BW6+/vTwi5du2a9ongxMbGalerb3/5m6ODYtLS0tRq9dSpUx0dHfUdCo3atHcWmgPhpw7Ky8ujt441iESi06dPWzwcMFFtba29vb3Ojjg6/Fej38a0TTSoEgMhGShZXl7eqlWrhQsXfvTRRyqVyt7eXl9JomvX6EBnjU+ZvAuGQ9WHRqWvmEaEvr6+9vb2dNCEYTqrNby/Bug8FBq1GX8GwTYI32VHR+MQQhITE6VSaVBQkFgszsnJwfOQrIuB66xIJGp8NjK8iUaW13llFIlE2pdCjZLau2b40tz442BMDYYPOBfhp59+unXr1p9//vnQoUPGbFpntSakIkpnhBq1NcmPDVgRJv70cHFxcXFxOXz4cHh4uKurq4uLi7+//+nTp/mPMwewDHt7++jo6OYw1LNly5b9+vWjQxaFjgWAEKZm+x4wYAB/Uq8XXnjh7NmzaK0DADQTwnfZcRYsWNC7d29PT0+xWHz37t1p06YhGwEANB8MtZAIIZWVlYWFhYSQNm3atG3bVuhwAADActhKSIWFhfz7RvQ78yazWyzveTCmsTEBANiQGzduCB2CXgx12clksvv379+7d69Hjx6lpaVhYWHLly9vZJ3CHnqpVNrMA2AhBsEDYCEGBMBCDIIHQP76rjqzGLpJU1BQcPz48UGDBu3atSszMxPDvgEAmhWGEhLtPHz99dfpN71dXFzwgD4AgOaDoS47iURSXFzs6ekZGBjYp0+fwsJCa39An+Bz8QkeAAsxCB4ACzEgABZiEDwA9jGUkJYtW9ahQwdCyJEjR+7fv6+zr7OmpoY/67OjoyOdXys7O7tNmzaNHAQBAAACYighTZ8+nT70rFOnTvoePJGSkkLLEEIuXLjw9ddfe3t7y2QyqVRaUFAwfPhwmUxmuYgBAKDpMJSQpk2bduDAgXHjxhkoM2XKlClTphBCiouLFy5c6O3tnZGR0blzZzoeb+zYsZGRkfg6LQCANWLo2i2Xy19//XUfHx9fX19fX9/+/fsbKLx169aoqChCSFFREZ0/nxAikUjwcEkAACvFUAtpz549RpasrKw8efJkQkICISQ3N9fPz4+uF4vFpaWl/JJSqfTmSxt7HowR5HbivXv3LL9RpgJgIQbBA2AhBgTAQgxCBaD9LGBmMZSQ6PAEY2zbto179qiXl1d1dTVdrq2tdXNz45e8ceOG3WK5gF9G6969u1CbZiQAFmIQPAAWYkAALMQgSAD8CyC+GFu/6dOn02eJaigsLIyMjNRYqVQqU1JSpk6dSl96eXnJ5XJCiEKhuHXrlkQiMXe0AABgDky0kD755JOIiAhCSK9evWJiYggha9euvXPnznPPPafdj/fDDz+EhoY6OjrSl97e3u7u7vPmzSsrK4uPj7dw5AAA0FSYSEhubm7p6eklJSX3799PTk52cnJatmxZ+/btNfrfqGnTpmmsWb58uVKptLOzw/g6AADrxURCGjt2bERERGRkZN++ffv27WtCDcbffwIAADYx0aTYs2ePk5NTcHBwcHDw+fPnm7Bmu8XyJqwNAADMh4mE5OzsLJPJzp07t3Pnzt27d/v4+MycOZM+qQ8AAJoJJhISp0OHDuvWrfvll1+WLl06e/ZsocMBAADLYSUhFRUVFRcXE0KePHmyYMGCnJycpv0qK/ruAAAYx0RC+v7770eNGvXyyy+/++67Q4cOtbOz+/zzz9944w195ZVKZWZmJn+WoOzs7Pz8fH3lkY0AANjHxCi7TZs2ZWZmtmrVasKECevXrw8ODiaEDBgwYPXq1dqF09LSNm/ePGTIEA8PD/o1WMz2DQBgA5hISFVVVa1atSKEtGrV6vnnn6crn3nmGZ0lP//88z179nDjvA3M9t3zYMzNlzZaYgcAAKDRmEhIhJCamhq6oFAouGVtZ86cGTp0qFwur6ioGDJkiJubm/Zs3/oma7pz506Th21Ys53MkakYBA+AhRgQAAsxYHLVejGRkJycnAYPHkyXX375Zbqgc9qF8vLyCxcu+Pn5OTs7z5o1KyUlxfBs3+r1I8lf95A8P7/DrbGY5jmZI2sxCB4ACzEgABZiwOSqhjGRkE6cOGFkSTs7u/Hjx9ObTD/99NPFixcNz/ZNqdePNGFcg91iuYWzFwBAc8bEKDvjeXl5ZWRk0OWioqKuXbs2dLZvu8XyBiUnjNADALAM4VtIeXl5M2fO1F4vEolOnz6tsVIqlXp5ec2dO7e2ttbf39/Dw8PDw8OY2b41Gkn81g9d5tZwLwmyEQCABQmfkCQSyZEjRwghiYmJUqk0KChILBbn5OQcO3ZMZ/m4uDiVSkV4N5lMm+2bJhuNTjmdeYj/Ep14AABmwkSXnYuLi4uLy+HDh8PDw11dXV1cXPz9/U+fPq1QKHSWF4lEGrlHLBab9uwJfgZCewgAQEBMJCSqoqKC/7K8vLxpHyqhXj+y8e0bpC4AADNhKCEtWLCgd+/e48aNmzhx4oABAyIiIszxwD0jc5KB7MVlI33JiVtDB5obk73QUAMAEP4eEmfGjBkvv/wyfepEmzZt2rZta6YNNWgUuDGFNYZFEK2kZeCDBPelAAAIIYy0kKZPn15SUkIIadGihUQikUgkNBsVFhZGRkZqFK78C39CB8OTq5qMpgqdCUNfptG5XmOsuUZLSLthhKYSADRDTLSQPvnkk4iICEJIr169YmJiCCFr1669c+fOc889t2fPHo3CYWFhvr6+hBBfX98ZM2YQS02uaripZGTyqLflhCQEAM0WEwnJzc0tPT29pKTk/v37ycnJTk5Oy5Yta9++vfa0C0qlslevXh9//DG3xsDkqkbSzjTG96EZGCBeb+EmofH1qSavHwDAYpjosqPc3Nz69u27atWqpUuX9u7dW+ckQGVlZY8ePTp69Ch/vgaNyVVN2DR3Kdd3TTfQd2e4TpOTBNdrhzYTADQTTLSQjGdvbz9p0iSlUpmenv7NN99s2bLF8OSq3EyCGs+fvf3/uhNCPD+/c/v/defPAk6XDc8Lfvv/dafD5wwUoJWkT3SgVdX7EcM04qFh0/8JIXaL5dwWuZVU084urFG5kZrtFMtMxYAAWIgBs33Xy8oSkqur6+TJkwkhYWFhkZGRlZWVhidX5U9zq029vrv2Mn+l/o/8b+Jwnc0X/oS+3bt3J+QO/d9AtYbRCnmdcv+rkEtydIH+rzGdcL2zCxvT6ffXW3f4tWmPLTQcv4AED4CFGBAACzFgtm/DGOqyS0lJMfAkJA0KhaKsrMzR0bGhk6uaQ709fvrKG8ncvXbGf1PK+PJcscY0DQGgWWEoIZ0/f37QoEGhoaEG7gNlZGRMnz590aJFUVFRsbGxYrHY29ubTq46bdo0A5OrNi3tjGIgJxl+y7TkZHi8X71D+PR9mdfAl3ybVkMHJQJAc8BQQvrggw+uXr365ZdfLlu2zNfXd+nSpfTLSXwBAQE7duxYs2bNd999N27cOLpy+fLlGzdu/OGHHyzfVcrPKA1NLfzyGp/VHhDR0EuzMcP/6l3ZoPyk82tV9X77ysj4jYHsBWDtGEpIVLdu3VasWNGxY8f09PRJkyZNmjRJu4yDg4PGGpMnVzWZmcZYN2G1dovlxs9dZLge4ws04Vd6jUxd2p8y7YMAIDiGElJJScm//vUvHx+fDRs2bN68+ezZs6dPn46Ojk5ISBA6NAsx0GYymb75IAxnjia/lBvT08gtGG66aQyI15eBNGYUBAD2MZSQoqKiXnjhhatXr3711VceHh505ZAhQzRGclspI7/qZKAAqS9PGJnDGnn/xsAjoxrZz1Zvuqr3Hpjh+gGAcQwlJDs7uxdffJHreQsLCyOEuLm5bdmyRdC4mpLOnNHkvX+sTdlg+IaW4U62Bs2FYcwmhMVIGABsYigh/fnnn46OjtzL6urqqqoqAeMxq8ZM+mDgs02Sioy/aJr77lST19BUmQm3qQDMgaGEpFQqlUol9/LPP/90dnbWV7i4uJg/Bs9Ms30LTmNcOGtNH+ti5OgME3r/6u2xbNB9LCQ5aLYYSkgrV67s27fvhAkTJk6c6OvrO3/+fH0lq6qqZs6ceerUKfpSJpMdOnTovffe2759u6WCZQLyk8nqHUPBL2n8oHmN9frutCHlAOjEUEIKCwu7dOnSunXr1q5de+LEiTlz5ugr+f77748c+b9rMTfb99dff71jxw6VSmWpeC3K5NzDeNJqTHiNn7iWv8bw4HUjc4mBplJDv62sLwbtgJHewGYwlJAIIc7OzvQBfa6urvrKnDhxol27dt7e3vRlk8z2bUWaZLqHJmdaAIY/1aA5L/hlzH00jEkA9WYLA+lQu6QxaakxjO/MbPLNmVyhRhuX+6dvQ43fIlgAQ5OrpqSkrF69mpspVSQSZWdna5QpLS3duXPnpk2bjhw5QteYNtu3ZZhvcl862/edO0/NVs7NLE7XE0K4MuSv2xjcdN1N++0c7UnNjZkTXWMWc+4j3N7x13PluQWNAvxIjGRgUJ8xA9nptvSV5NboC4nOzq5dIfdZ/gf5W9HeFjfRO/c/4R038r+Z2p+ar11j7naNqeK1Z3bn7w6/fsL7oTKA+13Q/rI2fzeNnFGeVqJ9AInW0bZbLOeOAxcDLWPa7PWNgdm+68VQQlq5cuXhw4c7dOhgoMzHH3/s4+Ozf//+q1evVlVV9ejRozGzfVtAdzNO7nvn6crv6NoWf+Wdp5sOmhOWG34kLjcjuPZ6u8Xy7lqTmtNlrk6NBe1IuI+Qp6dIf3qCc/5c7He6d+/Of8mr5P9eGt6pRjIyqXfXM9e79nrPz+n+PjWJu/Gb6P70HPC8lUTjXW7udrun5nHXfJfWo3EANeaYJ0/PRq91cp8KUue5+GuvCRengUq4YkTPgaW18TfEzYLv+fkd9XrtY2JRlt8iwWzfpnFxcfHw8BDzaJeZOHFijx49HBwcRCKRvb29g4MDC7N9C0Ln9Hca+H8AWvLpgjq3aExt+j5ioPuOe2l4c4b3yKy9fA0atmda+jSmx89Av5++d43sozOmWL3Jtd6voGmvMWE8pHYBY7YrLNbiMSuGElLv3r3rvQM0cODAsLCwsLAwPz+/fv369e7dW5DZvm2ACfnJmDLa1/16k2ITJgONBGbM5izzVWVBNOi6bMLHDXxE445Ovf2f2h2S2netmuQ+loEcZkwnrQmbM7JCk5OrjWEoIRUUFLz00kt9+/b19fX19fXt37+/gcJhYWH0SX1E0Nm+bYCRbQv+cmPmNdf3cWNyiVlx29XXLGM/SzF48WqqkMy6aw0apWJ8nTobmhppr97xLA3aqLnrsQCGEtL+/ftzcnIuX76clZWVlZV18eJF4z9r+dm+bYDhAXv6EoM5rssCph9jiulcMLmkzk8ZX9Lkx51YTJNf/sxxPTXcX2p4WKNpTU/t/sZG7pc5mnSCY+giLhaLU1JSIiIiiouLxWKx4OMRmhUTrm7W0pLQHknV+DtJGq1GY0aoG2iBGd861K6EfYxcMRvTxDG+TpM/VW+EBpKcOVp4QmEoIcXExFy6dKlFixb05dy5c4WNBwywogsi0XOhN7Kpoa9Lk66vNzNRJgwvNqanVLuM9nLjae+pzTD5Mm1a+6ZBIzsaNPLCmCRnFTmJoYR0/fr1Dz74gHvp5ORkw5Or2gabvEhR5hjsYI7eNo1cZZnWqr4BLM0QvxPP3Ld8dPYoNknN7GAoIVVUVPBfGp5cFYBZZrp7VG8BAxnC5G+AajcBtbfS5DmJnSSnPfDPklvU+a7x4/HYTz/aGEpIq1ev9vX1zcrKiomJ8fHxWbJkic5iSqUyMzOzqKiIv9JWZ/sGm2e4O7Gp6tdY1tcF19CxjtphG/5gg9KYWZt3TY6Fqz8LMTQSQzM1jBo16vz58wUFBYSQTp06cTeT+LKzsz/66KN+/frdvHmzT58+CxYsIITIZDKpVFpQUDB8+HCZTGbpuAH0aNC1r8kzU1N1MDa0sM57XfzJFyzcAFKbc7YOpjZK/spJ/CNsXVmKoYRECHF0dDQ81UL//v25Z0z0799/wYIF3GzfhJCxY8dGRkZi/DcwRb1+ZINm2ON/sMmD0bkJc1yzuOAbWrkxLSd9dXK7I1Q+0MnywbCz7w3FUEIKDg7m30ays7M7d+6cdjGVSqVQKL777rtx48YRXbN98ydrstXJVa0lABZiEDwAQWLQmLuWBqAxDy9/QeMjpmVQffiT4WpPqquzMH8uXe0paPVN3av9WaK1Lxq7bPx0wA2iMcGxdgGTN2dazD0Pxtx8aaMJm7MwhhLSgQMHuOWsrCz+S74LFy7s3Lnz0aNHM2fOJPXN9i34l5kEmUuRqQBYiEHwAISO4U7nzp35s992f3qmWvX67nb/myGX8As01da7PzXN61NvaTfRni72f3Pscu2e7k/Pn8vRalo9NbMwv37enLz/9xZ3cAy3wIxsfPA33V3XJLA6V/K3om9b/A8aqETDjRs3rKLZxFDvVgueoUOHXr58mf9Ec46/v/+nn3767bffrlmzRqFQGJ7tGwCMGaFg1q8xkfqG+RnzhV/DX8MypryBjzTJrTvT9qLeCnWuZ2cgYtNiKCFpKC8v1378K3+NWq0mhDTb2b4BTGBdFzJjpq5o/NgNw4MPjazfcKozMk7DGUjfF7TrrdCKTjpDXXb8e0i1tbUhISEODg4aZbZs2XLmzJl27drl5+dHRUU5ODhws32XlZVhtm+AepnQsPcwqr8AABMGSURBVGBHY9ouDb0xZkwHnQkDFtRazwYz0EGnsd5AxtL+uM7kKj3YoGAtjaGExL9pJBaLHR0dtcvMmTNn1qxZSqWSn6uWL1+uVCrt7Owwvg7ASjU001gycap1PcfSyBzAZ3haRX2tmYb28rE2yLBBGLqC8+8h6cxGlEgk0m45sTnbt+CPwxA8ABZiEDwAFmKw7QCMvGSbKQbjuw25AMz3PWhzV25uDLWQQkNDnzx5or1eKpVu2bLF8vEAgFmZaUIKM1Vo4TZZIwtbVx7iMJSQwsLCOnfuHBQUZG9vf+3atf3797/55ptCBwUAzZdQs0uYxiqCNIyhbq7du3eHh4e7urq2atUqICAgKyurZcuWLi4uLi4uQocGAPA/bF732Yyqoezo4GkWDBgwICsri3vp7++fkZHRmDtD/CkbAACays2XNvY8GCN0FCYSfLoAAxjqsnv99dd79+7t6ekpEol+//33V155pZHjFFg+7gBg3dbj8tL0GGohEUIqKysLCwsJIe3atXN1dRU6HAAAsByG7iERQo4cOZKQkNCyZUtXV9ecnByhwwEAAMthKCHFxMRcunSJewzS3LlzhY0HAAAsiaGEdP369Q8++IB76eTkVFVVJWA8AABgSQwlJP7DkAghf/75p7Ozs8m1CfJQc30bVSqVNTU1AsagUqnOnz9fXFwsVAA6Hzxv4Rio4uLikpISQQKo/IsFfhgM/ChmZmbm5eWZOwCdMdTU1FTy6JzO36wBEEJUKtXly5eb9oFPpsVggZ9DDZa8EJlGvHLlSqFj+B9PT8+IiIiCgoLz58+/9957b775po+Pj2lVyWSysrKygwcPPn78uG/fvk0bZ0M3mp+f/8orr9y+fXv48OGCxJCXl7dw4cLS0tJdu3Y9evSIe3yUxQLIzs5esmRJTU1NcnJyfn7+4MGDzReAvhioqqqqadOmtWnTxtvb2/IBjB49OicnRy6XP3z40NfX1/IBpKWlvfXWW2KxuLKy0qxHQF8MycnJW7dulcvlcrl89erVAQEB5ntejM4AqqqqoqOj1Wr1qVOnTp06NXKkeb+7ozOGysrK6Ojo//73v99//31NTY3JV7mGsuSFyHRqllRXV9+8efPmzZsVFRUmV3L27NklS5bQ5dDQUKVS2UTRmbjRhISE3bt3r1q1SsAY6LJSqRw6dKggAXD8/PzMF0C9Mbz11ltr1qzZvXu35QOoq6t77bXXzLfdegOorKwcP358XV2dgDFwHjx4MGXKFMsHcPLkSW69v7+/+QIwEMO+ffs2b95Ml1966SXLXKDUFrwQNQZDXXa3bt1ydHSUSCQSiYQb2mAC7YeaN1GAJm50zZo1rVu3FjYG+o2uwsJCT09PQQJQqVQ1NTX/+c9/6IPnBYnhxIkT7dq1M3fLQF8AZWVljx49Onr0aEZGhiABnDlzZujQoXK5PCUlxdydRfX+Dm7dujUqKsryAQQEBNy9e/fEiROJiYkRERHmC8BADLW1tdwFwd3dPTc316xhcCx2IWoMhhLS9OnTm6Se3Nxcbjpw7Yeam4kgG21QDDU1NStWrFi6dKkgAVy4cCEhIUEulwcFBZkvAAMxlJaW7ty5My4uzqxbNxCAvb39pEmTlEplenr6rFmzLB9AeXn5hQsXRCKRs7MzfYaL5WOgKisrT548GRYWZvkAHBwchgwZsnnz5iNHjgQHB5svAAMxDBgwICkpKSUlZcOGDbdv3xaLxWYNw7owNFPDtGnTDhw40Pi/oAV5qDkLT1I3EINCoVi0aNGcOXPMOp2SgQD8/f39/f1VKlVYWFhgYKD2A0TMHcPHH3/s4+Ozf//+q1evVlVV9ejRw0x3cfQF4OrqOnnyZEJIWFhYZGRkZWVlY7oBTAjAzs5u/Pjx9Cr8008/Xbx4kfv73WIxUNu2bfvHP/5hpk0bDuDUqVOPHj3avn3748ePZ86cuWnTJvP9quqLoXv37lu2bLl69erkyZNPnz7do0cPMwVgjRhqIcnl8tdff93Hx8fX19fX17d///6m1SPIQ821N3r37l0LbNeYGJRKZXx8vEwmCwgIECQA7QfPWz6GiRMn9ujRw8HBQSQS2dvbmzUjGv5JUCgUZWVlBp74ZaYAvLy8uN7CoqKirl27mikAAzEQQpRKZUpKytSpU823dQMBcBcEFxeXtm3b/ve//7V8DIQQNze34ODg4uLi9u3bo4XEx1ALac+ePU1SjyAPNdfYaGZmJv3fMls3HMOGDRtycnKSkpKSkpIIIfHx8Wa6j6IvAO0Hz5tj64ZjGDhwIFemqqqqd+/eFg4gIyPjiy++aN++/f3792NjY813GdIXgFQq9fLymjt3bm1trb+/v4eHh5kCMBADIeSHH34IDQ01Xz42HMC4ceNef/31CxcuPHnyxNPT06x/rRo4CAkJCbW1tU+ePFm7dq35ArBGws9ld+PGjT/++GPo0KHz5s376quvmqpaQR5qzt9oTU2NuX/r2IxBZwAqlUrjwfOWj8GS9AWgUCgscxD0BUBbq5b5vWD2LFjy4qAzBkGuTlZB+BaSvb39gQMHevXq9dtvv2nc/Gzbtq3J1QrSEOZvVJBsxEIMOgMQiUSW/PVj8yAQQiyWkvUFgLNALHtx0BkDuun0ET4hPf/8861atYqIiCguLg4PD+fWi8XiEydOCBgYAABYkvBddpzY2NgvvvhC6CgAAEAYTCSksWPHRkREREZGmmkgLAAAsI+Ju2p79uxxcnIKDg4ODg4+f/680OEAAIAAmGghcYqLi9euXXv06NGBAwe+9957nTp1EjoiAACwECZaSJwOHTqsW7ful19+Wbp06ezZs4UOBwSjUqlSU1O5ZRNqUCqVaG2bwJLHjX+WAQhTCenWrVvcslQqPXr0qIDBgE4jR47cvn079/Lf//63maYKVSqVdNq9mpqaQYMGmVBDTU3Nxo0b+WtGjhzp/5eUlJSmCdT8rl+/npCQYLHNaRy3kpKSpppkUht3lrXt3bt306ZNZtouMIuhhGS+n3toKpWVlV999RX3nL3y8nKFQmHWLTo6OjbVhBf//e9/T58+nZmZmZmZOXHixCap0wJWrFjx2muvCbX1lJSUESNGWH674eHhKSkpFnuiIzBC+O8hcZpqclUwqwULFrz11ltbtmzhr0xNTf3444+dnJyee+65tWvXurm5KRSKpUuXSqXSzZs3f/DBB4GBgQsWLGjdunVBQcGff/755ptvbtq06dGjR25ubps2bXJwcNi0adORI0dEItHzzz/Pn09FoVCMGTNGLpefOHHinXfeIYTU1dUNGzbs3XffTU1N3bhxo5OTU7du3d577z36ZcPk5ORvvvmmTZs2Q4YMMbwjGhGOGDFCu8K9e/d+9dVXbdq0CQ4OvnLlyqefflpZWTllypQDBw4QQiorKydNmkSb8hqfValUCxYseO65537//fd79+6tWLGCzmp69OjR9evXV1ZWvvDCCw8ePFixYgWdvWbt2rVdu3adMmWKRpA3btygx4S+LCwsdHV15Qaj3rhxw5jZcjX29ObNm9qHWt9x279//5YtW0pLSxMSElQq1YMHD1577bWJEydqnC+FQmH4/NICbdu2vXnz5uPHj99//33+ZJU6T+W4ceNSUlLmz59f7w6C7RD0aUxPGT9+fM+ePXv37t2vX79+/fqZ+0luYAJ6UqKjo/ft26dWqxMSEk6ePPnHH38EBQWVlZWp1eqjR4/OmDFDrVbX1tb27t1779699IO1tbU9e/bMzc1Vq9VJSUkjR46srKxUq9Vz5849duyYWq1+9OgRLRkZGXnx4sXa2lq6LW6BUiqVEyZMyM3N/eOPP0JCQuiDHFetWrVt2za1Wl1QUODv70+rSkpKopHwg09PTz979uzZs2fLyso0ItSu8N69e4GBgY8fP1ar1Zs3b6a1VVRUBAYG0o9UVFQMHjxY52fp/mZnZ6vV6tzc3JdeeokWCwgIePjwId2R7du3v//++3Q5MDCQrtewd+9e/hPVqqurP/nkE7qh9PR0uVxO19fV1b3zzjtxcXHx8fH0/4KCAu5TGnuqcagNHLfr16/T5X379n3wwQfcKdB5vgyfX1rg+vXrtFr6oEh6cnWeSrVaffLkyfnz52sfE7BhDLWQmmpyVTC3d999d/r06UOHDqUvs7KyBgwY4OrqSggZM2bM4sWL6XqRSMTvGXN2dqZ/znt4ePj5+Tk7OxNC2rdvX1tbSwgpLS09fvz4tWvX/vjjj4qKCn2bXr9+/UsvvcTdYly3bh0hJD8/v7y8nBBy/fr1IUOGuLi4EEJGjhypPdPHjz/+SOfOmTt3rkQi4UeYlZWlUWG7du0GDx7cqlUrQsioUaPS09P1RaX92SlTpjg7O9OnxUskksLCQlrshRdeoBNiiUSi8PDw0NDQJUuWZGZm9unTR+dEWZcvX+Y/JsPR0fG111778ssvpVKpo6Mj15kml8v/+c9/nj59etSoUWlpaWPGjNGoh7+n2oda33Hbs2cPfYrdgAEDEhMTCSEvv/wyba5pV2L4/NICdFZfb2/vqqqqqqoqe3t7nUePlvf29ubfV4bmgKGEJBaLk5OTt23blpiY6OHhkZOTY+6He4JpunTpEh0d/e9//5ubm8vOzo5fgI6Lo5cbDn8WNY0Z1a5evbp06dLXX389Njb24cOH+rZ7/vz5K1eu7Nixg7708PAICQmhy/yn/tCFZ599VruGdevWcbPJKRQKjQg1Krx165bh2mpqavR9lr+P/J3lH6gWLVoMGzYsLS3t+PHj2p11lPakZ46Ojr169Tp+/PjKlSu5laNGjSKElJSUXL161cnJSbsebk/1HWrtPVWpVMePH6eDDjp16nT48OFDhw7Fx8dPnjzZz89PuxID51fnfql5XzjReSpVKpXGzxXYPIYGNcTExFy6dInrH587d66w8YABc+bMuXPnDv3btk+fPllZWfTqfP78+Z49ezZ0Bs9ff/114MCBI0aMcHFxuXnzps4yjx8/fvfddz/66CP6slevXnfv3vXx8QkICAgICKB3Yp5//vns7Gz6LNRTp041KAbtCnXW5uDg8OTJE5px6dNu9AWjTSqVZmVlVVZWcmumTZt2+PDhrKwsrrmpoW/fvtnZ2fw1p06dsre3f+edd77++mt+VadOnfrb3/5WUlLCPRROJ52HWueenjhxYtiwYVxGFIlE48aN+/e//52ZmWnM+dJWVVVFWzz5+flisZj7Tdd39HJzcz09PY2sHGwDQy2k69evp6enR0ZG0pdOTk5VVVW04Q8MWrt27fjx4wkhnTp1io+P/8c//tGpU6eioiITHvHy97///auvvlq4cGFlZWWfPn10lklMTCwvL1++fDkhxNfXd/78+XFxcRERET179iwsLHzjjTf8/f0lEsmECRMmTpzYtWvXLl26NCiGbt26aVfI1daxY0dazMHBYdKkSRMmTOjatWuHDh3on/Danx0wYID2Jrp37z5jxoxJkyZ5eno+++yza9as6d27d2FhYWhoqL4U3rt3b/oUKyovL0+lUtHxEa+99trWrVu5v9u2b9++adOm7Ozsw4cPa3fZcXQeap3H7ccff+QqP3Xq1MaNGzt27Pj777/HxcV169at3vOlrWXLlp9++qlarb5+/fq7777Lrdd55AkhV65cMdNTfYFdQt/E+j/9+/dXq9XTpk0rLCxU/3X/HKyFUqmsra1tTA2mfby2tpbeZudHUldX15gw+BUqlUqlUllRUcEfIlFXV6dzE9rBaNM4UBEREbdv3zZQPjo6+pdffqk3bC4eY/Zd56E2fNy0z2+Dzhc3OKWurk7fIeIfvbq6ujFjxugc6AE2jKEW0urVq319fauqqmJiYm7dukX/FgZr0fgnHpn2rCDtTzVtGDpr0/c8G2N2gTtQOTk5mZmZLi4u3bt3N1B+1apVqamp9T7ilgvJmGft6IzT8HHTPr+mnS8D4fErzMjImD9/fmOeiAbWiKGENGrUqPPnzxcUFBBCOnXqhJm/gR1isbjeLzY11J49exwdHT/55BPDxbp16/bqq6827aYtTyQShYWFGV9e3001sG0MTa566tQpPz8/OvYUAACaG4ZG2R04cGDMmDH+/v6xsbH5+flChwMAABbFUAuJqqysLCwsXL58+c8///zrr78KHQ4AAFgIQ/eQCCH5+fnvv//+pUuX2rdv/+233wodDgAAWA5DCal///6+vr7Lly83POgIAABsEkP3kJKTk5VK5eTJk6dOnXrnzh2hwwEAAIti7h4SIaSwsHDDhg0HDhzAPSQAgOaDoS67oqKi5OTkvXv3Pnz4sEuXLriHBADQrDDUQgoLC4uIiAgLC3Nzc2vkl+0BAMDqMJSQyP9v7w5tAAahIAy/sAJMwQhICJOxE0vgGQDLAAhMBUnDBC1p/88999wlZ04kxthaU0oZY3LOa4oGAPAHB1V2IYSUknNORGqt3vtSyttPAQAeclAzNudcaSQi1lqt9T6ABgD4toMCaYyxn733e5AUAPB5F7D30lz7moc4AAAAAElFTkSuQmCC" }, "execution_count": 111, "metadata": {}, "output_type": "execute_result" }, { "name": "stderr", "output_type": "stream", "text": [ "Warning: Imaginary parts of complex X and/or Y arguments ignored." ] } ], "source": [ "P2 = abs(orig / L);\n", "\n", "P1 = P2(1:n/2 + 1);\n", "P1(2:end-1) = 2*P1(2:end-1);\n", "d = datenum(date);\n", "d = d - d(1) + 1;\n", "figure\n", "plot(d(1:500), P1(1:500), 'b-');\n", "\n", "\n", "figure\n", "tiledlayout(3,1)\n", "nexttile\n", "% Original Transformed Data\n", "plot(d/86400, orig(1:length(d))); title('Original Transform');\n", "nexttile\n", "% All Positive\n", "plot(d,P2(1:length(d))); title('Positive Transform');\n", "nexttile\n", "% MATLAB has built-in functions or other stuff\n", "pwelch(flow)" ] }, { "cell_type": "code", "execution_count": 112, "id": "2c7acb18-129d-4f66-891d-c7c3e862d835", "metadata": {}, "outputs": [], "source": [ "% MATLAB - Not necessary to run, also SUPER SLOW!!!!!!\n", "%raw_data = importdata(filename);\n", "%[nr nc] = size(raw_data);\n", "%date = zeros(nr,1);\n", "%Flow = zeros(nr,1);\n", "\n", "% Scan the rows for a string, string, string, float, and a string\n", "\n", "%for ii = 1:nr\n", " % row = textscan( raw_data{ii}, '%s%s%s%f%s');\n", " % agency = row{1};\n", " % % Check if the first string says USGS, we know we in the data!\n", "% if strcmp(agency, 'USGS')\n", "% date(ii) = datenum(row{3});\n", "% flow(ii) = row{4};\n", "% end\n", "%end\n", "\n", "%[nr nc] = size(f);\n", "%date = zeros( nr, 1);\n", "%Q = zeros( nr, 1);\n", "\n", "% Import the date and time and plot\n", "\n", "%date = f{:,3} ;\n", "%flow = f{:,4};\n", "\n", "%figure\n", "%plot(date, flow)\n", "%xtickformat('dd-MMM-yyyy');" ] }, { "cell_type": "markdown", "id": "bfbe07af-7e11-4be2-9def-94726e4df699", "metadata": {}, "source": [ "# Python\n", "\n", "# Raw Data\n", "Lets say we are given the task of analyzing the hsitorical data fro mthe water gauge station at the [Prado Dam in Los Angeles](https://waterdata.usgs.gov/monitoring-location/USGS-11074000/#dataTypeId=continuous-00065-0&period=P7D).\n", "\n", "We want to download the [data](data/PradoDam.txt), clean it, and do some spectral analysis on it to see if there is anything interesting.\n", "\n" ] }, { "cell_type": "markdown", "id": "4812aee1-75bd-471d-9cbe-01a881ae97e1", "metadata": {}, "source": [ "## Reading in the data\n", "\n", "If we do a quick check of the file, we note there is a giant header, then some columns of data.\n", "```\n", "# ---------------------------------- WARNING ----------------------------------------\n", "# Some of the data that you have obtained from this U.S. Geological Survey database\n", "# may not have received Director's approval. Any such data values are qualified\n", "# as provisional and are subject to revision. Provisional data are released on the\n", "# condition that neither the USGS nor the United States Government may be held liable\n", "# for any damages resulting from its use.\n", "#\n", "# Additional info: https://help.waterdata.usgs.gov/policies/provisional-data-statement\n", "#\n", "# File-format description: https://help.waterdata.usgs.gov/faq/about-tab-delimited-output\n", "# Automated-retrieval info: https://help.waterdata.usgs.gov/faq/automated-retrievals\n", "#\n", "# Contact: gs-w_support_nwisweb@usgs.gov\n", "# retrieved: 2020-04-29 18:30:02 EDT (caww01)\n", "#\n", "# Data for the following 1 site(s) are contained in this file\n", "# USGS 11074000 SANTA ANA R BL PRADO DAM CA\n", "# -----------------------------------------------------------------------------------\n", "#\n", "# Data provided for site 11074000\n", "# TS parameter statistic Description\n", "# 8183 00060 00003 Discharge, cubic feet per second (Mean)\n", "#\n", "# Data-value qualification codes included in this output:\n", "# A Approved for publication -- Processing and review completed.\n", "# P Provisional data subject to revision.\n", "# e Value has been estimated.\n", "# \n", "agency_cd\tsite_no\tdatetime\t8183_00060_00003\t8183_00060_00003_cd\n", "5s\t15s\t20d\t14n\t10s\n", "USGS\t11074000\t1940-09-30\t51.0\tA\n", "USGS\t11074000\t1940-10-01\t47.0\tA\n", "USGS\t11074000\t1940-10-02\t47.0\tA\n", "USGS\t11074000\t1940-10-03\t47.0\tA\n", "```\n" ] }, { "cell_type": "code", "execution_count": 3, "id": "a436a73c-82e4-4a58-8733-d40df2a1720b", "metadata": { "scrolled": true }, "outputs": [], "source": [ "# Python\n", "import os\n", "import pandas as pd\n", "import numpy as np\n", "import scipy.io\n", "from scipy.interpolate import PchipInterpolator\n", "import matplotlib.pyplot as plt" ] }, { "cell_type": "code", "execution_count": 20, "id": "fe54ec2e-1a18-4fbc-acd0-66869a1686e8", "metadata": { "scrolled": true }, "outputs": [ { "data": { "text/html": [ "
\n", " | agency_cd | \n", "site_no | \n", "datetime | \n", "8183_00060_00003 | \n", "8183_00060_00003_cd | \n", "
---|---|---|---|---|---|
0 | \n", "USGS | \n", "11074000 | \n", "1940-09-30 | \n", "51.0 | \n", "A | \n", "
1 | \n", "USGS | \n", "11074000 | \n", "1940-10-01 | \n", "47.0 | \n", "A | \n", "
2 | \n", "USGS | \n", "11074000 | \n", "1940-10-02 | \n", "47.0 | \n", "A | \n", "
3 | \n", "USGS | \n", "11074000 | \n", "1940-10-03 | \n", "47.0 | \n", "A | \n", "
4 | \n", "USGS | \n", "11074000 | \n", "1940-10-04 | \n", "55.0 | \n", "A | \n", "
... | \n", "... | \n", "... | \n", "... | \n", "... | \n", "... | \n", "
29061 | \n", "USGS | \n", "11074000 | \n", "2020-04-24 | \n", "140 | \n", "P | \n", "
29062 | \n", "USGS | \n", "11074000 | \n", "2020-04-25 | \n", "124 | \n", "P | \n", "
29063 | \n", "USGS | \n", "11074000 | \n", "2020-04-26 | \n", "120 | \n", "P | \n", "
29064 | \n", "USGS | \n", "11074000 | \n", "2020-04-27 | \n", "137 | \n", "P | \n", "
29065 | \n", "USGS | \n", "11074000 | \n", "2020-04-28 | \n", "156 | \n", "P | \n", "
29066 rows × 5 columns
\n", "