Amplitude and phase reconstruction simulation

sine_wave_reconstruction.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from numpy.random import randn\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "        Amplitude: 2.69926\n",
      "        Phase shift: 2.66615\n",
      "        \n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAD8CAYAAABjAo9vAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJztvXeUZFd5r/3syrmrc+7piRrNKGskRJBIQpKJBmOSjc2HjcDXYHD4fO2rBU6Xe+3PYWEbbIwBG2MsggkmySCBAAklRiNNzrmnc6qqrhz298cJXd1TPdOh4qn9rNVruqtOnbPn1Dm/8+53v0FIKVEoFAqFdbDVegAKhUKhKC9K2BUKhcJiKGFXKBQKi6GEXaFQKCyGEnaFQqGwGErYFQqFwmIoYVcoFAqLoYRdoVAoLIYSdoVCobAYjloctKOjQw4PD9fi0AqFQtGwPPvss9NSys6rbVcTYR8eHmbv3r21OLRCoVA0LEKI86vZTrliFAqFwmIoYVcoFAqLoYRdoVAoLIYSdoVCobAYStgVCoXCYihhVygUCouhhF2hUCgshhJ2haIOefjIBGORZK2HoWhQlLArFHVGJlfgvZ/fy78/tapcFEUF2Htulq8/N1LrYawbJewKRZ0xvZCmIGE+ka31UJqWf3niHH/+0LFaD2PdKGFXKOqMyVgagFgqV+ORNC/RZLahz78SdoWizpgyhV1Z7LUimsqRyOTJ5Qu1Hsq6UMKuUNQZU8pirzmxpPZQXUjnSOfyRBrMLaaEXaGoMyZjKUAJey2J6MIeS+X4xKOnee3HH6vxiNaGEnaFogaksvkVXS3KFVNbpJRE9XMfTWU5PxPn4mySeLpxHrRK2BWKGvDnDx3jlz/9dMn31OJpbUllC2TzEtC+AyM66dJ84+QVKGFXKGrA6akFLswmSr5nWOwLmRyFgqzmsBRgWuugC7vulhmZK/191SMbFnYhhEcI8YwQYr8Q4rAQ4k/KMTCFwspMxdJEUzmkvFy4DWGXUhN3RXWJJouFPUskkQHg0lxzWexp4BVSyhuBm4D7hBB3lGG/CoVlmV7IkC9I4pn8ktellEzF0rT6nIByx9SCSHIFi72ZXDFSY0H/06n/qPljg5HK5jk4Eqn1MJqCfEEyG9es8mLrUPs7RyZfYEtnoOT7ispT7IqJJLOm0Bdb7M+en2Munqn62FZLWXzsQgi7EOJ5YBJ4WEp52aqQEOJ+IcReIcTeqampchxWUUa++fwor/v44+w9N2u+VijIkq4CxcaYS2QwXOfRZZEvRqjj1k4/oCz2WhBNLp7z0fkkxi1gLJ5Gklne+k9P8pnHz9ZieKuiLMIupcxLKW8CBoDbhRDXldjmU1LKPVLKPZ2dneU4rKKMjOqVBD/2yEkAEpkcb/jET/mTbx2p5bAsyfRC2vy9WERg0b9uWOwq5LH6GA9bu01wUV8wdTlspsX+1JkZcgXJaB27Zhzl3JmUcl4I8ShwH3ConPtWVJaZBW1a+fipab59YJRHjkxw8FIEv9te45FZj+nY4hR+uatlShf9LR3KYq8VxnfSE/JwcVYT72t7guwfiZDO5fnpqWkAJvTZVT1SjqiYTiFEWP/dC7wKaNyyaE3KbDxDf9hLT8jD+//jOb7x/Chuh01VGKwASyz25a6YqPbe1i5lsdeKSDKL12mnPeAyrfJdfS0AjM2neNwQ9mh6xX3UmnJY7L3A54QQdrQHxZellN8uw34VGySayvLzH/8p/+dN13PHlvYrbjsTT9MX9vDJX76VZ8/PIdGaPTx+cnrJdsfHY+zoDiCEqODIrc1SV8xS4b40n8TnstMf9mrv6xZ7Ll/gwCVtcXtXbwiPU82kKkU0mSPkdRD0OMjpiyG7+0IA7D0/x5mpOG6HjcmohS12KeUBKeXNUsobpJTXSSn/tBwDU2ycgyMRzkzH2Xdh7rL3pJT8y0/P8qEvPke+IJlZyNDmd9EecHPP7h7u3d1Du9/FXGLRbfDjE1Pc+7Gf8OSZmWr+NyzH1EIah017MEaW+diPjkW5pieI22HDaRemK+ar+0Z40z88wZv+4Qn+7Ntq3aOSRFNZQh4nQbfTfO26fs1i//sfamtQ9+zuIZrKkVwWrlovqMxTC3NkNArARORyy+K3vvg8f/KtI3zj+VEmoilm4xnaA+4l27T4nKRzBVJZ7eL9yt6LABwdi1V45Nbk1GSMAyPzTMcydAbd+F32Ja4YKSXHJ2Ls7AkhhCDocZqumOcuzNPidbKzJ8jpqYWVDqEoA9FUlpDXSdCz6NDY2RPk9s1tZHIFXrS1nTu3dwCLUUz1RlkXTxX1xZExXdiX+QInoym+tX+UGwda2D8S4eJsgrlEhna/a8l2rT7t77lEBr/bwcNHJgA4O62EZT18+BuHOTsdZ0dPkA79IVrsipmIpplPZLm2NwhAyOMwLfbDo1Gu6w/R7nfz/MX56g++iYgks3QG3AQ9msUecDvwOO18+b0vNLd57KQWsj0ZS7Op3V+TcV4JZbFbjEQmxwcefI5z0/FFiz2WIpsv8M7PPM0Tp6Y5rL/+ltsGAe0BUJDQtkzYw17twp6LZ3no4BjpXIEWr5MzU/Eq/o+sgZSSQ6MRxqMp9p2foyPgIuRxLrHYj45r38vOHs2fa1js2XyB4+Mxdve10Bf2MhZJqhoyFUTzsS9a7C1e52XbdAU9AEzUqZ9dWewW4+BIhG/t1yJaTulT9olIiguzCR47OU1nwM1mPZTuVdd288DXD5kZp8tdMWHdYp9PZvjOwXG2dPi5aTCsfOzr4OJs0rS+F9I5OgJuFtK5JXHsx3QX1zU9msUe1C32U5MLZPIFdveFiCazZPOS6YU0XSFP9f8jFuYTj57Sa/hkaSkS9rDvcmHvDmn3SqnIGCkl49EUvS3eyg74CiiL3WIYFQO/um+EfEGypcPPZCzNWd3KfvLMDAcvRdjc4acr5CHsc5rRFpe5YvzaBT2fyHJuOs51/S1s6fQzFkmRUMWp1sShUe0cB92aWHQE3ZdZ7MfGo/SHvaaFGPQ4iKay5gxrd18L/a2aWIzMJ3nPv+3lw99Q6SLlYDKa4m8fOcm/PnGO+YS2eBrSXTGlhL3F68RVFBlz6FKEd37maVLZPD8+McWL//yHXFyhemc1UMJuMYyLyUiDfvnOLnIFaUbGjEW0OFwjfKs/7DUX4y53xWh/z8QzjEWS9IW9Zkbk2WnljlkLh0cj2G2Cd7xgCICOgJsW71JhPz4eM611MFwxOQ6PRvA67Wzu8NOnh0FemEnw4xNTfPvAqHLLlIFPP36WXKGA26FJohHuCIv3QTFCCLpDbrN2/g+OTvLYyWkuzCY4PRWnIKnpIrcSdotxYTZBf9jLcLuPgNvBnk2tADx9dhY9wo5EJm+Gb/WFveZDoD2wTNh1S+X05ALZvKQ/7DHdOMrPvjYOj0bZ3hXg1df3AtoDNeR1mq6YeFpzuSwVdocu7FGu7Q1itwlT2H9ycopMrsBcImsukivWRySR5d+fOs/rb+zjzbcOABDyOAkYPvYSFjtofnbDx35uRrsfJqIpsyzE6Hzt/O9K2C3GhdkEm9p9/OkbruOB11xLT4vmhz0wMs/uvhY6g5pv8Do9k85IhIHFKBgDj9OOx2kzF2H7wl5T2JXFvjYOXYqyqy/EjYNhvvGbL+ZVu7r1qJcshYLkC0+fJ1eQ3Lu7x/xM0ONkIZ3jmbOz7BluA9Djqx384OikuZ2RCalYH4dGIyQyeX5xzyDvuXMLHQEX27uDZlRMuMTiKWh+9uXCPhlNmyGQtawlo4TdYlyYTTLU5uOuHZ28/fYhU9izeclwh58X6hmoxa4Y0HyGTvvll0Orz2VahH1hLx6nlhV5RsVSr5qJaIrphTS79YfpTYNh7DZByOukIDVX1z8/dpYXb2vnpsGw+bnbhlu5tjfEA6++lt++e4f5el/YSySZxeeys60rYNYuUawPQ4h7WzwMd/j52QN3c+um1isunoJmsRslIM7phs5kLF1ksddO2FVUjIVIZHJML6QZbPOZr3UE3Aih+dyH233cu7uH3X0hWnV/ujG1X+6GMQj7XIzpCU59+ir/5g6/sthXSS5f4IGvH0QIuGNL25L3jMW5zzx+lqlYmr99201L3r9zeycPffDySqh9YQ/HJ2Ls6g1xw0CYLzx9nlQ2r8oMrJInTk8TdDu5fkB70BpCbMxmjXIZfS1edvYEuWmwteR+els8xNI5PQ9EWyuZjKVMsa9lj1RlsVsIoxLdUJGwO+02MxlmuN3Pdf0tvPelW833jSiL5RExBsY01O+yE/JqdkBf2MN4ncbv1ht/8/AJHjk6yZ+8frdpsRsY5/PzT57j+v4WczZ1NYyH8XX9LbxgSxvpXIGjys++aj78jUN87JET5t+T0TQep42Ae6md63XZ+e8P3cXtm9uW7wLAnF19Wc/IBt1i12sBGaWwa4ESdgthhDoWCzssxtwOd/gu+0xfWHPVLI+IMTBCHvvCXtOS6Q55mIqlyatojKvy4xNTvGhrO7/ywuHL3jMs9ngmzzteMLTqwmqGsO/uC9GlW5mqCufqkFIyFkkRSy+G604tpOkKetZc2O6moTBuh80U9na/i0tzSWbjGew2wXgkVbN7RAm7hVhR2PUsueESqc8dfjdux6JVv5wW71KXDUBXyKP5hhfqt2xpvTC9kGagtXSiSkifDflcdl53Y9+q97lVDzm9eShs7iOiWuitilg6RyKTJ14s7LG06YZZC26HnVs3tTIRTSME3LqplWN69vA13UEzkawWKGFvcCaiKTNJ4uJsgqDbcdliz1C7jza/q6RVbrMJ/u7tN/NrL9lccv9GU+ViYe8Orpx1p1ikoFfNXOmhaVjsr7+x7zI3wJW4Z1c33/vQXWzrCprJTErYV4dxrySKqjJOxtJ0rvAdXQ3DfdYb8jDY5iOVLQCaNQ+187MrYW9wPvjF5/jAg88BcGY6zlC777Ip5QdfuZ0vv/eFK041793dYyYeLcd4SPSHF9PXu0P1XSejXogks+QK8rJSDQYDrV5+65Xb+cArt69pvzabMOPdDWFXTa9Xx3hEM0aWW+xdoXUK+1ZN2Ic7/KbLExb977WKjFFRMQ3OiYkFEpkcuXyBw5civGJn12XbhH0us+7LWjE+t8RiN4S9TkuW1gvGNLxjhYgjm03wO6/aUfK91eK02/C57MpiXyXGor8h7Olc3qzmuB5uGAgTdDvY0R00C4NB7YVdWewNTCSZZTaeIZUt8NjJaWbiGW4YaLn6B9eAIUrFiUwdARdCKFfM1TCiI9YrGqsl5HESSWZJZfP83+8eZSGt6vishDHLTGTzFAryslDHteJy2Pja/3gRH7p7u7mQDVpIcNDtqFn2qRL2BubCzGKRoS88fQGA6wfCK22+Lu7c3slf/eKN3Da8GPLl0EMo67k1WD1gNAjvWKdorJYWrybse8/N8U8/OaMSlq7AuJ6TISUks3lT2NfrigHY3h0k7HOZ+2jzu3DabfSGPcoVo1g7RhozwA+PTeCwCXYW1RopB067zayfUUxxOrWiNIuumMoLezSVNY83X9TOULGU4vyLeCa3aLEHNl4C2SijbFjuYZ+rZmGoymJvYM7rwn5df4iChB3dwaplH3YHPcoVcxWmF9LYbWLFWiPlIuR1EkkuipSKaV+ZYmMkkc6b1RnX64opJuh24HHazH0ZM6lasGFhF0IMCiEeFUIcEUIcFkJ8sBwDU1ydczMJukNu001Sbv/6legKeeq232O9MB3T2g3abGtLfFkrIa+DaHLRYp9Twr4i45GUGcK7kNYehkKsXFJjLQgh2NkTMmfNtRT2crhicsDvSin3CSGCwLNCiIellKqVeoU5PxNnuN1vCvr1VRT27pCb6YUM2XyhZPEwhWaxrxTqWE5avE6iyWyRxa5cMaXI5QtML6S5eaiVZ8/PkchoFnubz1W2a/hL770Dux5W3NAWu5RyTEq5T/89BhwF+je631IcHYvyw2MTldh1Q3JuJsFwu5+7tndy147OkqGOlcIIeTTERHE50wvpFUMdy0mL10ksnTPDT+eUsJdkeiFDQcIWvfR0XLfYy+GGMXA77Dj0h0TY6ySZzZPJFcq2/9VSVlNLCDEM3Aw8Xc79GvzH0xf43S/vr8SuGw7jotzU4aM94Obf3n17VXssLvZ8VO6Y5ey7MMfx8RjTC5mKhzrCYpKS0fxEuWJKYyycGsl4cb0aaqUWt40GHbWw2ssWFSOECABfBT4kpbys1JwQ4n7gfoChoaF1HcPnthNP56++YRNwXg91LFX/pRoYyRhjkRQ312QE9csffPUABalb7BUOdYTF0gRGeWXlirmcXL7A1/eNALCtSxP2RDrPbDzDpvbLi+OVg+JyD+WcFayGsljsQggnmqh/QUr5tVLbSCk/JaXcI6Xc09l5eY3p1RBwOcjkCzWZ2tQbKxX8qhZbOv0IAScmYjU5fj0zEU1zanKBdK5QNVdMMSoq5nJ+9yv7+dyT53n77YPcNqzVV19I55hZSK9Y2XSj1LJAWzmiYgTwGeColPJvNj6klfHphZKSGWW1G4kPxRmh1cTncrC53a/qgC/DSFE3qHQMOyztyel32ZlPZJFSlVQu5rGT0/z8TX383zfdgM+l6chcIkM8k1+xF8FGKa7j83tf2c/vfPn5ihynFOWw2F8MvBN4hRDief3n1WXY72X4XVqM9kJGpUyPzifxOu0rtu2qBtf2hlQj5WUY2abGTV0VYS+y2Ld1BcjkC0uqFzY7+YJkPpExZ7cuhw2X3cZFfdbb5q+Qj73IYt93YY595+cqcpxSbNjHLqV8HKhsoK6OX7fYE6oWBmORFL3htTcHKCfX9gb5zsExYqms2fi32TGihD7wim08fGTC7C1bSUJF5357d5D9IxHmEhnzfml2osksBYnZDhLA77ab7sxKuWKKhX1sXmu6USjIiuc1QINlnvrdmsUeV9YIl+aTNXPDGOzSRevYuPKzGxhJQnuG2/jSe19YtTh2g+36wqDysy8yE9dmUcUC7nM5uDinuTMrtQ5ifC8XZhNa2GO+YBaGqzSNJey6byyuLHZG55P0tmy8vsVGuLZXE3bDz/6j45P81feO13JINWej1QLXg8epuRYAtndrwq5i2RcxzkWxsAfcDvO7qpTFbpRUPl5k+IzMVacoWGMJu1sJO0Ampz35+2pssfeEPIR9To6MasL+nQNj/MtPz9Z0TLXGEItKLciVQghByOvA7bAx2Kr5kVUs+yKzusXeWtSTwOderKnUXiEfO2hW+7Elwp64wtblo6GE3ecyXDHNLewT0RRSQl8VE5JKIYRgV2/ItNhjqRzxjFbnulmZXkgT8jiqVozNIOR10hl0m41RIspiN5kr4YoxWhE6bNpDsVK0eJ1L+p5Wq1VeQwl7wLTYm9vHboQ61tpiBxhs9ZmJMbG0ZiU284N3qkpJScsJm8Ku+XWVxb5IaR+79uBt9bsqGoBgxLILASGPo2qumIZaNjfi2BNNLBwAoxHt4ugN19bHDloMdTSliUg0qX0vC+lc00bJTMeqU0ZgOf/vvTsBza8bcDuUj72IuXgGn8u+ZBZlrNdV2mVmLKB2Bd10BT1cUsJ+OT79i1loeotds5Br7YoB7cJNZQukc3liusA38xrI1ELajBaqJkZTZdAakKuomEVmE5kl/nVYXK+r1MKpgSHsvS1eekIeTk0tVPR4Bg3lirHZBD6Xvenj2Efnk7T6nHhd1fXjlqI4bTqa0r6XWKp5v5/pWLomFnsxrT6XstiLmItnLhNwY/G00uGoi8LuYaDVy8hcoipZwQ0l7KDFnzazDxc0Ya8H/zosTZtetNibc0aVyuaJpXNVL/i0nLDPqXzsRcwmskuSk0CrOwXVc8X0tHjob/WSyhbMKJ1K0nDCHlAVHjk3k6h5cpJByKPdIBPRNNm8ZokspJtTVBb7Z9beYlcVHheZjadpW1Z6w1clV4yxmN3b4jHv2WosoDacsPtcjqZePJ2MpTg7HefWTa21HgqwaJEUx+c26xqIkVVYa4u9VfnYlzAXz15WD8aoO1VNH/tAq492v6sqrsqGWjwFLeRxoYl97E+fmQXgji3tV9myOiwK+6IVspBqTlExLPZqFP66EmGfi2gqS74gsVehLkk9k87lWUjnaPMvtdiNxdNKu2IMl+nmDj+7+kI8++FXVfR4Bo1nsbvtTV257qkzMwTcjqoUl1oNpYS9WWv5GNUC+1tr6yYL+5xIWZs64PWGMXNZ7mM3wxBDlX0I79nUysO/fRfX9VevHzE0oLD7XY6mDqd76swMtw23mn0Va02ohCumWaNizk7HafE6aa1hKWVYTJ1XkTGL5QTaloU7vmhrOx9/x83cMlRZl6YQgu3dwYoeoxQN54rxN/Hi6WQsxempOG/ZM1jroZgYhY4uzhZZ7E364D07HWdzh7+mpZRhccGuWRdQZxbSfP25S6RzBbM42nJfusNu47U39NVieFWh4YS9mcMdnz2nFeq/fXNbjUeylBav02wU7HLYmnYN5Nx0vC7WPkyLPZ7luQtzjMwled2N1hWx5Xz68bP8449OL3mt0ouk9UZ9zOfXQMCtuWKasfXXed2HazTjrRdavJpPF7SwrmYU9mQmz2gkxeaO2jQXL8YQ9vlklk/++DQf/q9DNR5RdXni1DR7NrWy78Ov4m23DbK5w183eR/VovEsdredgoR0rlD1Cnq1ZmQuQYvXWXd1WIwOPnaboDPgZqEJfeznZuIAbO6svbCH/YuumNH5FPOJLPOJjFn50cpEklkOXorwgVdsp83v4s9/4YZaD6kmNJzF3szNNi7NJRmoccRFKYwF1KDHQcDTfK6ybL7A2Wld2OvAYg+6HdhtgrlEhjG9YNz5merUAa81z5ydpSC1xdFmpvGEvYlL947M1b4dXimM0LGQx6nlGTSRxf7wkQlu/tOHeeToBADD7bUXdiEEYa+TiWiaab25tjGjsDpPnJ7G7bBx01C41kOpKWURdiHEZ4UQk0KIijvz/E3abENKychckgG9Q0490VJssTdZAtnpqQUW0jm+tu8S3SF33TSQDvucHBuPmn+fm24Oi/3J0zPcNtyG29FcbtrllMti/1fgvjLt64r4m7Qm+1wiSzKbr0tXzGUWexMJe7QoCage3DAGrT4XJ8YXS8Q2g8WeyuY5Nh6rm3IbtaQs5oWU8idCiOFy7Otq+N3NWZPdSACqdVZjKYzWYkGPA7/bQSKTb5p09kgyS4vXyTXdQe7c3lnr4ZiEfS4y+QKgRSo1g7AbbeeGO+pvVltt6mPeuAZ8+uJps9VkNzqv1LPFHvQ4CerVHuOZnBktY2UiySxtfhdfft8Laz2UJRRnv96xpZ0fn5iq4Wiqw4h5jyhhr9riqRDifiHEXiHE3qmp9V9kRt/TZpruQ9FFG66/i9Z0xXgdRYvbzfH9RFM5MyqonjBqo7T6nOzsCTIbz1i+doxRq2dQCXv1hF1K+Skp5R4p5Z7OzvVPWY0mtM0iHAaX5pME3Y6KdlRfL8UWu/ngbZLImEgya9akrycW64B7GdZ9/+ct7o4ZmUvistvoqnHZ5Hqg4cIdW7xOXHYbY5FUrYdSVUbmEvS3emteh6QUIXPx1NF0M6qo7mOvN8JezWLvC3vMEEwj1t6qXNTvEVsTrO1cjXKFOz4IPAlcI4QYEUL8Wjn2WwqH3caWTj/HJ2KVOkRdMlKnyUmgLc6FfU52dAcJeJpP2OvSFVNksW/u8ONx2njuwnyNR1VZ6vkeqTZlEXYp5dullL1SSqeUckBK+Zly7HcldvYEOT6+KOxz8QwPH5mo5CFrSiZX4Mx0nC2d9VUjxiDocfL8R+7hrh2dZmbw8fEY3z88XuORVRYppRkVU28Y5QN6wx5cDht7NrXx1JmZGo+qMjx/cZ5kJs/IbEIJu07DuWIArukJMRZJEdGL6H/0u0e5//N7zWbKVuPERIxMrsD1VS7Wvx6MqJj//Z2jvO/fnyVpwaYbqWyeo2NREpk8uYKsS2HvafEAi5mwL9jcxvGJmOVK+cZSWd78j0/wZ985wkw8oyJidBpU2DXL9cRkjOmFNN98fhQpMdOnrcahSxEAbhiof2EPFGVeFqSWmWk1/uHRU7z+44+bdVjqMaxzc4ef/3zfC7l3dw8AL9jSjpRaLRUrMRFNkStIvvyzi0B9hgPXggYVdq0t3LHxGF985oKZiDGjNxO2GgcuRQh5HAy11b81EvY5efeLN/Pnb7oe0GYbVkJKybcPjpHNSw6Pain79WixA+wZbjOTxG4cbMHtsPHUGasJu3bP5wpa3ejBBrhHqkFDCntfi4eg28Hec7P825Pn6dWnnVa12A+ORLh+oKUuI2KWI4TgI6/bxS/cOoDTLjgxYS2L/eTkAmemtOiSo2PaQ6tehb0Yt8POLUOtPHF62lK9DCb0Bi9GGLSy2DUaUtiFEOzoCfJfz48yG8/wR6/bBcBM3HoWezqX59h4lOv7G6tandNuY0tHgJMWs9gfOri4IHx0TLPY6zG3oBSvvbGXY+MxPv/U+VoPpWwYFvtvvnwbWzr8dAZUDDs0qLAD7NAbxD7wmmt5xc5uAKZj1rPYj4/HyOZlQ/jXl7O9O8DJSWtZ7A8dGuOWoTAuu40jY/XtilnO228b4uXXdPK/v310SeXHRmYimiLodvCbL9/GD373pQ0xq60GDSvs737xMH/8ul2860XDuBw2WrxOS1rsB/WF00aIiFnOju4gF+cSlomMSWRyHBuP8fJruugNe5iKaddbPS6elsJmE/z1W27CZoMv/2yk1sMpC5OxFF0hzUpXor5Iwwr79u4g73rxZvPLbA+4mLGgj/3YWIyg29GQvsMd3QGkhFMWsdqNej2bOvxLGp7UY4LSSrT5Xdw0GOZn56yxiDoRTdMd8tR6GHVHwwr7cjr8bqYtGBVzbDzKNT3BhrRGtuvuMqtExiwWmfKazZGNNnSNxG3DbRwejVgiO3gimlLCXgLLCHt7wMVM3FoWu5SSY+MxdvYGaz2UdbGpzYfTLjhlkVj2C4awt/lMi72RrHWD24bbKEh47sJcrYeyIaSUTEbTpitGsYhlhL0jYD2LfTSSIpbKsVOP2280HHYbA60+LlikkfLF2SRep512v8tIL69ZAAAgAElEQVRseNKIwn7LplZsAn7W4MlK84ksmXyB7qCy2JdjGWFvD7iYT2TJ6slKVuCYHnVxbYNa7ABDbT7T0m10Ls4lGGzTKmwaFntLg4Q6FhNwO9jVF+KZBvezT8S0GHZlsV+OhYRd+3LnLOSOOaYXOjNCOxuRoTafZeqAX5xNmNm/piumQSJilnP7cDvPXZhvaD+7EcOufOyXYxlh79A7xlgp+/TYeIyBVi/BBhUPgE3tPqKpnFmwrVGRUuplYTVh7w1rYtIoMezLec0NvaRzBb57YKzWQ1k3RtapcsVcjnWEXe+aYpVY9kJBcvhSpGH96wZG7Y7zs41ttc8lsiykc+b/x+2wc+/ubl6wpb3GI1sftwyF2dLp58t7L9Z6KOtmMqpcMSthGWFvNy32xhf2fEHyP796gDPTcV6+c/1tBOuBTe26sDf4AmpxqKPBP71zD2++daBWQ9oQQgh+8dZB9p6f40yDRi1NxtKEPA48Tnuth1J3WEfYdR+7FZKUvrL3Il95doTfeuV23nH7UK2HsyGMxsKNvoB6cW4x1NEq/MIt/QA8dKgxG6LMLGToULVhStJ4S/orEPI4cNltTFnAYj8+ESPgdvA7r9pR66FsGL/bQUfA3fAhjxdntaxTKwl7V8hDwO1oWGNoJp6mTZ+pK5ZiGYtdCMFQu4/TFkhfH4+kzA44VmCozWsJi73N71rSSMQK+N124g0aGTMbzyhhXwHLCDvA7r6Q2fygkRmLpMwa81ZgU7u/8YV9NrHEv24V/G4HC5nGFfb2gBL2UlhO2MciKWYbPJZ9PJKix0KxuYNtPkYjyYau8nhxNsGAhdwwBgG3oyEt9kJBMpfIKot9Bcoi7EKI+4QQx4UQp4QQf1COfa6H3X1aadsjDWy15/IFJmPWsthv3dSKlPDU2ZlaD2Vd5AuSS/NJcyHYSvhdjSnskWSWfEHS5leLp6XYsLALIezAJ4CfA3YBbxdC7NroftfD7j4t5vvwaKQWhy8LUwtpChJ6Wqwz7X/B5jY8Ths/OjZZ66Gsi4loimxeNkTP2bXidztYSDfeTMoo+NeuLPaSlMNivx04JaU8I6XMAF8E3lCG/a6ZsM9Ff9jb0H72sYiWdGEli93jtPPirR08enyqIfttmjHsbdZ52BoEGnTx1HC3KldMacoh7P1AcfraiP7aEoQQ9wsh9goh9k5NTZXhsKXZ1RdqaIt9XBd2K0XFALzsmk4uzCY4O914GagX9QYblnTFNKiPfVbPMFfCXpqqLZ5KKT8lpdwjpdzT2Vm5bMrdfSHOTMd5+MhEQ1qHVrTYAV52TRcAjx6v3EO9UlyYTSAEZnMNKxFwOxqyEJjpilFRMSUph7BfAgaL/h7QX6sJb751gG2dAd7zb3v59GNnazWMdTMeSeJx2hq2uNRKGM0pDo7M13ooa2ZkNkFvyIPLYakgMkCz2NO5ArkGK3c9u6BcMVeiHFfqz4DtQojNQggX8Dbgm2XY77oYaPXx3Q/eyVCbj+cuNl6HGC2G3duQrfCuxuYOf4O6YhKWyjgtxq8nXMUbbAF1Jp4h4Hbgdqg6MaXYsLBLKXPA+4HvAUeBL0spD290vxvBabfR5ncRSy1OMaWUDdEKzGox7MUMd/g4Ox1vOBfZxdmkZYU94NaEsVGSlC7MJNh3YU5lnV6FsswtpZTflVLukFJulVJ+tBz73ChBj4NokbA/eXqGN/7DExwbr++IGatlnRYz3O4nmsox10C12TO5AhOxFAMWzDqFYou9MYT9rx8+zq9+5hnGoykl7FfAek5DnZDHyUJqUUAmY9oqej1npRYKkomoterEFLO5ww/QUO6Y+WQGKRerh1oNQ9gbZQF1LJIils6x99ysimG/ApYV9oDbscQVE9VFvp7T2qfjaXIFaVmLvRGF3biGQh5rFf8yCDSYxT6lG2gFqRZOr4RlhT3oWSbsSV3Ys/Ur7Isx7Nac9g+2+bDbBOcaUtitFaVk4Hc1lrAbXZMA2lSo44pYWNidJLN5snoYl3GDJurYYrdqDLuB025joNXL2QZqbm0YBEGLW+yNUFYgns4Rz+TN0iHKFbMyFhb2pZaI4YpJNYTFbk1hB20B9eholNd//HE+/sOTtR7OipyaXCCXL5gGQSM3FL8Sfj0qphEsdmOd7G23DfLCLe3cvrkx+81WA2uaIUBAF/ZYKkfY5yKa1C7cevaxj0VSuOw22nzWtUQ2d/j58Qkt+7ReLa7ZeIZ7P/YT/uYtN5rXi1Ut9kZaPDXcMMMdfh68/44aj6a+sebVyuJil2GpG//WsytmPJKku8WNzWa95CSD7d0BAFp9Tiai9dnGcDaeIV+QjMwlcdm1Sa1Vhd3tsOGwiYay2DuD1oxQKifWvFpZnDobU2kjpr2eXTFG1qmVefOtA9zQH+Y/njnPw0fqs4xvQk/WmYtn8LnsCLG4yGg1hBB66d7GEfauoHVdleXC8j72BV3QY40QFRO1bnKSgdth5/qBFjqDHmbi6bqsUWLM6uYSWaKpHEG3w9KzqEYpBDYVS+O0C1p91lzvKCcWFnbdYk8brpj6joqRUjJmsSbWV6I75EZKmF6ov4Qx02JPZIimspZdODVolIbWk7EUnQG3JesolRvLCrsRxrXoiqlvi30ukSWTK9Br0ToxyzGm0xNFccn1wqLFniGWylnWv26g1WSvz/uimKlYms4muT82imWFPVgUFZPK5snktCl/vUbFjEW0Zg5WTU5aTndIWwCrS2HXRW4uniGWylo2OcmgUVwxk9E0XWrhdFVYVtg9Tjsuu41YKrckA7VehX3c4slJy+nWLS9jQayeMFwxs/EmsdgbpKH1ZCylhH2VWFbYwSgrkDXdMFCfrph8QZr1U5pF2Nv9LoRYmiJeL8T1h380lWM+kSVksaYny2mE9niZXIG5RFZFxKwSS5siRr0YIy3c5bDVncUeS2V55V//mMlYGo/TZtkqgstx2G10BNx1abEXXyNjkaTlLfaA2173rpipBT3UMdQc98dGsfQVG9AtdsMV0x1y153F/t2DY0zG0vzOq3Zw5/YO7BYOq1tOd8hdlz72eFHTiYK0bnKSgd/tIJ7JI6Ws24iTAxe1lopbOwM1HkljYOkrNuh2aha77orpDno4N5Oo8aiW8tVnL7Glw88HXrGtbm+qStEV9JhrC/VEYlmEiNXDHQMeB/mCJJnN46vTRKxHj08S9Di4ZShc66E0BJb3sS+kc2admO6Qp64yTy/MJHjm3Cy/cOtA04k6aBZ7Pbli9l2YI1+QJLLLhb0+xa5cDLdrdfLPTNVn1U0pJT86PsVdOzpx2C0tWWXD0mcp6HHqUTGaxd6lu2LqpefmN/dfQgh44839tR5KTSjOPr00n+QDDz63ZKG7mpyeWuBN//AEjxydIJHOLclutLrFvkOv33NiIlbjkZTmyFiUyVial1/TVeuhNAwWF3YHUT0qxm4TdATc5AuSTJ2ksZ+aXKA/7KUv3Byx68sxsk/HIin+9pETfGv/KM+er03D8UtzWh7BZDRFIpOnv6jHqVW7Jxlsavfjsts4MbFQ66GU5EfHtWqgL93RWeORNA4bEnYhxC8KIQ4LIQpCiD3lGlS5MFwxkWSWoMeB16nVnk5l6kPYZ+IZOpokCqYUL9jchsMm+F9fP8jX9l0CFgW22hguoblElkQmR5vfjdthVHa0tsXutNvY0umvW4t93/k5rukOqqqOa2CjFvsh4E3AT8owlrIT9DiQEsYjaUIeJ16XJuyJbH2Edk0vZOho4vZe27qC/M49O3js5DQAdpvg0nxthH3KFPYMiUwev8tOq14X3+oWO8CO7iDHx+tT2GcTGRXmuEY2JOxSyqNSyuPlGky5MSytkbkEIa8Dny7s9RLLPr2QbmqLHeC9d23ltTf08hsv20pf2FNDi12LzokksiQyWnRIq94IxOoWO8A1PUEuzSfrMp49ksjSYvEksXJjaR+7sSh0bDxG0O3Eo7ti6iGWvVCQzMYztDexxQ6alf7xd9zC795zDf1hb51Y7Dl8Lru5gGr1qBjQLHaAk3XojplLZMzZk2J1XFXYhRCPCCEOlfh5w1oOJIS4XwixVwixd2pqav0jXgO3bmrj07+yh7DPyeZOv+ljrweLPZLMki9I2v3NbbEX0x/21YWPPZ7J43PbafW7sNuEOdOzMtfowl5vfvZCQRJJZgmrGuxr4qqmiJTy7nIcSEr5KeBTAHv27KlavOHdu7p5+n+9EoCDIxGgPiz2mbgmJM1usRfT3+plIpYikyvgclR3MjmtC/tMPE0mV8DndDAQ9tIdbI763wOtXrxOe91FxsRSOQoSwspiXxPWn2Oide0BFl0xdWCxGw0mmt3HXsxAq1cPf0yySU+aqRaGxW5kwvrddt79kmF++Y5NVR1HrbDZBAOtXkbm6iszez6p3Sdh5WNfExsNd3yjEGIEeCHwHSHE98ozrMpgLp7Wg8WuhP0yBvR4/mq7YxKZHAvpHG6HjWxem0x6XXaCHieDbb6qjqWW9LfWbo1jJeYSWsJaq18J+1rYaFTM16WUA1JKt5SyW0p5b7kGVgm8NYiKOTYeLZnpqlwxl2MkBY1UWVyMhdPiAlNWbV59JfrD3pqtcazEfEIzgFq86j5ZC5aOilmOsXharb6nR0aj3Pexx9h3YTGb8jsHxvjYIyeYjqURArXaX0Rvixchqm+xG8JuRFEBTbFgupz+Vq+ZoPWWTz7JX/z3sVoPiXnDYleLp2uiuYS9yq6YCT02eqqo0NV/PnuRTzx6ipH5JG0+V1OV6b0aLoeNrqC76u4Aw7++oydovlavVQ4rSb/uCjs9Gedn52d59Nhk1cdwfibOK//6R2Y5Z8NiV4una6OphN1lt2ETVK3C44JeB764UfBYJEU2L/nJiSnlhilBLdwBpsXeVSTs7ia02HVh/9HxSaTUQh+r3Vnp6FiU01NxDl3SItgMH7tKUFobTSXsQgi8TnvVXDHGTZEoatwwpkddaOUE1MLpcgZafVW32Kdiaew2webOxUicZnXFAPxAt9QLElNgq4VhBBmzqEgyS8jjUDPbNdJUwg7gdTmq5oox0rONHpqJjFaQzKBZ2uCthf5WL6PzSfKF6pVWnoyl6Ai4aPcvzqCacfG0K+jBYRPsH5nHadeEdP/IfFXHYHSvKs4EbvWrme1aaUJhr17f0wXTYteOZ1jrRlGpdnXBXkZ/2EuuIM3aLdVgKpamM+gm5HFi5CJ5m9Bit9sEvWEPUsK1vSEGWr3sv1hdi924Z4zvfz6RVTHs66DphN3ndFRN2E1XjP7v2Lx2sd53XQ9AU1d2XAnDHVBNP/tkLE1X0IPNJkxfbjNa7LDoZ9/ZE+TGwTDPX6yyxW4Ie1Sz2OcTGbVwug6aTtg9LnvNXDGjEU2sfv7mfuw20VTJL6vFTFKqop99KpamU3eLtfpcCAEeZ9PdGgBm05drekLcNBDm0nyyqg3Hl/vY51WdmHXRdFdvu99VtQt1Ib3oW4dFi/3WTa384Hdeymuu763KOBoJM0mpShZ7viCZXkibTRxavE58TntT1IcpxUCRxf7ynVrHogefuVC14xvGkOljj2eUK2YdNJ2wb+30c2Y6XpXFuQW9f6dhhYxHk3QEXLgddoY7/Koxbwl8LgdtflfVLPbZeIaCxGzk0Opz4m1SNwzADQNhgh4Hu/tCbOsKcve1XXzuiXNVc18mihZPc/kC0VROuWLWQdMpy7auAJlcoSrFjgxBT+odm0bnU/S2NGd/07VQzVh2Y5HOcMVs7Qww2Na839Err+3i+Y/cY4rpe1+6lblElq88e7EqxzdmuZl8gQuz2j2qXDFrpymFHbRG0pXG9LGnjaiYJL0tnooft9HpD1evyqAx5Tcs9v/5czt58D13VOXY9YgQYknM+J5NrezuC/Gt/aNVOX5xQpRRQliV3Vg7zSfsnVp2YTWFvdjHbixOKVbGqDJYqnhauTEW6ToD2gPXabeZ5Z0VmtDfsaWdAyMRMrnKN4GPp3NmXRgjIqdNhQWvmaYT9hafk46Am9NTlRf2eJHFHktliaVz9CiL/ar0h72ksgVm45mKH8uw2I3FU8Xl3DLUSjpX4OhYtOLHWkjnGO7QMoC/um8El93GLZtaK35cq9F0wg6wrctfFYs9VmSxG5E4yhVzdYzImItV8LNPxdIE3Y6mTEhaLbdsCgMsqVJaKeLpHJt1YZ+KpXnxtnYC7uZdzF4vTSrsAU5NLlR0qp/NF8ypayKTVx2T1sD1/S24HTb++vvHKVQ4emkqlqYzpL6TK9Hb4qW3xcO+C5VPVopn8nQG3fj1B+09u3sqfkwr0pzC3hkgmsoxtZC++sbrxHDDtPqcpHMFc8qv/IVXpy/s5Y9et5vHTk7zTz85U9Z9SymXzNYmYykzIkaxMrcMtbLv/BynJhcYi1RmJmUYQwGXg66QByG0KB3F2mlKYd+ud2R/7sI8qWyeTzx6irky+3NjesnerqDmejESbpSwr4633z7Ii7a2lz3M7u9+cIq7/+bH7NcX5ow6MYorc/OQloV699/8mN//zwMVOYZhDPndDrZ2+rljc7t5/yjWRlM6r24bbmOg1csnHj3FoUsR/v6Hp7DbBO976dayHcOoUtcVcnN8ImaG76nQrdUhhODa3hDPXZhHSlmWTNAfHpvgYz84AcCz5+e4cTBs1olRXJl7dvXwrf2jzCezXJytTCiqEUUWcDv42NtursgxmoWmtNhdDhu/9YrtHBiJ8PFHTwHwg6MTZT3GwjKL/eJckqDHgcvRlKd8XfS2eEhm80ST5Wn28FffO8H2rgAdATcHL0WIp3MkdJ+u4soMtfv4r/e/hFfu7F7SEaycGPkePredgNuhFk03wIZURgjxl0KIY0KIA0KIrwshwuUaWKV54y39bGr3EXA7ePvtgzyr+w/f8c9PcaAMNagN68NIfBmZTSg3zBoxsnRHy+DTlVJyfibOi7d1cNNgCwdG5heTk5Swr5qukJt4Jl+RzkoLRa4YxcbYqPn4MHCdlPIG4ATwhxsfUnVw2m184ddfwNd+40W89bYhChJ+5TNP88TpGZ46M7Ph/RvWhyEaI/NJJexrxIj5H49svGhbJJklnsnTH/ZyfX+YM9Nxnjk3C6CSxtaAsdBcCavdSORTlvrG2ZCwSym/L6U0Ht1PAQMbH1L1GGj1sb07yA39LXQEXIzqAjIR3fhFu5DWCoB1hzRxyuQKqrHGGukLa+euHBa7sXg90OrlhoEWpIS/eOgY/WEvtw2rBJjVYritNhJR9r7PP8uX916+KG4unjZxEbZyUU6H77uBh8q4v6phswneeHM/u3pD9Ie9GyrrOxfP8Ef/dch8OBRP89XC6droDLixiUWLPV+QfPQ7R9a1eLco7D6u628BYCae4ZfuGFJVNteAKezrtNiTmTz/fXic7x0av+w9owCYstg3zlXPoBDiEaBUlsADUsr/0rd5AMgBX7jCfu4H7gcYGhpa12AryQOv2YWUkrd+6imze8t6ePjoBJ978ryZPVcccdGmOiatCYfdRnfIw6hex/7cTJx/fuws3SEPv37nljXtyygD3B/20up30dviYWYhw1v3DJZ93Fama4PCblRsPD4RM1/L5Ar89PS0abH73CoLeKNcVdillHdf6X0hxLuA1wKvlFdI5ZRSfgr4FMCePXuq16l4DQgh6A55OLhs8TSdy/Mn3zrCm27uZ89w2xX3cWxMu2DPTsfxOG0EPYunuE1Z7Gumt8VjJsQYtWPG1uFzvzSXxOeymyVg3/WiYfJSqobia6TV58JuE+vuSXtuJg5oM6iFdI6A28FDh8b44Bef5z49y1RZ7BtnQ2dQCHEf8PvAS6WU1amzWmG6g24eiaaXxE7/8TcP8+AzF/E57VcV9uJCSQG3Y4n1oRZP105vi9c8pzN6WYb1LKaOzCXoD3vN7/S9ZcxZaCZsNkFHwLVui/28LuwAJydi3DzUypkp7bWfnprGbhO4VUjwhtnoGfw4EAQeFkI8L4T4ZBnGVFO6Qm6S2bxZwOtb+0d58Bltoedq1QallBwbjxLSrfSA24HLbsOh17duV66YNdPb4mE0opXwXbTY176Yemk+aRYXU2yMzqB73cJ+biZh1ns/obtjLurJe7F0Dr+redsSlpONRsVsk1IOSilv0n/eV66B1QojimVSX0D96alp2vwuru9vYfoqwj4RTTOXyPJLd2wCtHhcIQQ+vaCRWjxdOz0tHlLZApFkltm4JibrsdgvzScZUMJeFjoD7nVHxZybjnNdfwsep81spDEyu/igVm6Y8qDmPMswFjuNqJaLcwmG2nx0Bt2msKzE0XHNZfCyHZ1c0x00hdxIuGj3K3/uWjFizEfnU8zoD9aJWHpNPWsX0jnmE1n6w76KjLHZ6Ap6lgQYzMUzJSulPnFqmqje99fg/EyCLR1+tncFl1jsRtMmnxL2sqCEfRndeqaosTg0MpdksM1Hm99l+nhXwlg43dkb4hO/dAt/+obdAKbFrqJi1o5Rv34skjRdMfmCZHqVFmO+ILkwo031lSumPHQG3czEM+QLkrPTcW776CP85OT0km3GIkne8emn+ezjZ83XUtk8o5Ekm9p97OgOcnw8RjqXZzya4sXbOgCVdVoulLAvoyu0aLHnC5JRfQrfHnAxU8Iy+dwT5/jJiSlAWzjtD3tp8TrZ1hVgS6fWX9Xn0nztftXMYc0YYnxpPrlkjWO1kTFv+scneN3HH9f2pTJMy0Jn0E2+IJlLZHjy9Ay5guTwaGTJNk+e1rK3jSqaoC1gSwnD7X6u6QkwGUtz6FIEKeG+63pw2gUBFepYFpSwLyPgduB32ZmIphiPpsjmJYOtPtr9LjK5glnPArT4249+9yh/9u0j5AuSZ8/PcW1v6LJ9+lx22vwutSi0DjoDbjxOGxdmEkwvZEwLfjyS5OhYlHQuv+JnY6ks+y/Oc8NAC2+8uZ/dfZd/N4q1YyQpTUbT7D2vlWVYnjRmlOU4eCliGkPnprVtNrX7uE2PLvuiHpiwvSvIS7Z1sL0rWPn/QBOg5j0l6A5pPsQR/WIdbPMyGdWegbPxDEGPFgt9bDxKJlfg5OQCf/eDk1yaT/KHr9552f46gm6y+co3ArYiQggGW31cmE0wG09zfX8LY5EU+y7M85v/8RzvetEwH37trpKfPTauucY+8IptvGJndzWHbWmMGj6nphZ49rzWLu/8zHJhn8VuE0wvZBiPpuht8Zox7MPtfkJeJ60+J9/cPwpo99hn33WbMn7KhLLYS9AVcjMRTZk9NwdafaZ/fLrIz250UXfYBH/7g5N0Bd3cW6KV1x+/bjcff8ctVRi5NRlsM4Q9w9auAC67jQefuUC+IPmPpy+s2CTlyKi2mL2rt6Waw7U8Nw6E6Q97+eSPTnN+JoEQixmloLnNLswmeM31vQAcGNHcNOdnEoQ8DsI+J3ab4M7tnaRzBVx2G91BjxL1MqKEvQTdIQ8TsRQXZ7WLti/soUOPaCn28z5/cZ6OgJuf0y/gt98+hLNE3ZHOoFtVENwAQ20+Tk8tkM1LOvxuelo8xFI5ekJavfbPPXmu5OcOj0Zo97vMBXFFebDbBO94wRBH9MSxF21tZ3Q+ac5Kn9L96+9+yWbsNsFBXdjPzcTZ3OE3BfylOzoBbR3FZlOiXk6UsJfg2t4QF2eTPHZyiu6gB7fDblrsM0XRGM9fnOemwRbe/eJhru0N8UsvqL8aOFZgsM1HNq/5adv8LtMV8M4XbuLua7v4/JPnS37uyFiUXX0hZQlWgLfeNojTLnA5bLz6+l4KUivbAPD9I+O0+13c0N/Cju4gBy4tWuyb2v3mPu7ShX2wTYWhlhsl7CV40y39OGyCfRfmGWzTLG2j5K4RSx1JZDkzFeemwTA3D7Xy0AfvNCNqFOVlsChMsT3gMhdQX3N9LzcPtTITz5DKLi6ifu6Jc/z4xBQnxhfYVWIxW7FxOgJu3nnHMK++rodtevTXhdkEl+aTPHxkgl/cM4jNJrihv4WDI/NkcgVG5hIMty+KeGfQzVv3DJo1YhTlQy2elqAr6OHua7v578PjDLZqF6LHacfvspux7Pv1QmE3Dapa3pVmqEgM2v1ufv6mfrpDHoY7/GaRtVgqh8dpZzyS4o++eRibgIKEXSoSpmJ85HXaorVR4uHCbMKMhvnlO7TZ601DYb609yKPnZyiIFlisQP8xZtvqOKImwdlsa/A23W3SnEaelvAZWaffmv/KB6njRsH1cJcpTEerqB9By/f2cX/evW1AEXCrmU4fv+IVufbaKunQhwrT3fQg8th4/BolC/+7CKv2tXNgP6d7dmkGT5f3TcCwHCHcrtUA2Wxr8Cd2zp4/8u38fqb+s3X2v1axt3ofJJvPH+JX3rBJjP0UVE5/G4H7X4tQWx5F6qgWzv/Mb15+H8fGmdrp58H33MHjx6fZKvuJlBUDptNMNjq5cFnLmC3CX7jZdvM97Z2BmjxOnnkyCRwucWuqAzKYl8Bm03we/dew7auRWFo18sKfObxsxQk/Pqdm2s4wuZioM2Hz2XH41yamWhY7AvpHHPxDE+fneW+63roCnl4621DauG0SgzpC6B/+HM7uWlwsae9zSa4dVMrmXyBgP6AVlQeZbGvgfaAiydOz3BqaoE33NhnTjcVlWd7V4BEUdavgTFjiqWy/PDYJPmCLJlLoKgsb799iN19LfzaSy43dm7d1MoPj02yqd2nHrRVQgn7Gmjza7Xat3T4V8x2VFSGB159LfFMKWHXLuFoKsfFWa1K4HV9at2j2tyzu4d7VnigGn72YeWGqRpK2NfAzUNhtnUF+Oyv3karmlJWlVa/q+Q5D3kWfexziQxhn0slu9QZNw6G8bvsXNOj6sBUCyXsa+De3T1qml9nBIqiYubiWbOnqaJ+8DjtfO+376JD9ZetGkrYFQ2N3aZ1qDIsdtUwvD5R61HVRUXFKBqeoMdBLJVlNp5RLun1LLYAAAc8SURBVDKFAiXsCgsQ9DiJpbT2d63KFaNQKGFXND6axZ5jNqEsdoUCNijsQog/E0IcEEI8L4T4vhCir1wDUyhWS9DjZCKaIpMrmA3EFYpmZqMW+19KKW+QUt4EfBv4SBnGpFCsiaDHYTZ6UIunCsUGhV1KGS360w/IlbZVKCpFyOMgndOaPKhwR4WiDOGOQoiPAr8CRICXX2G7+4H7AYaGVEMKRfkoLsTWpnzsCsXVLXYhxCNCiEMlft4AIKV8QEo5CHwBeP9K+5FSfkpKuUdKuaezs7N8/wNF0xN0L9onavFUoViFxS6lvHuV+/oC8F3gjzY0IoVijRj1YgC1eKpQsPGomO1Ff74BOLax4SgUa8dwxQgBLV7lY1coNupj/3MhxDVAATgPvG/jQ1Io1oZhsbd4ndhVATCFYmPCLqX8hXINRKFYL0YhMBXqqFBoqMxTRcNjlO5VoY4KhYYSdkXDY7hiVKijQqGhhF3R8ARNi10Ju0IBStgVFsCw2FVlR4VCQzXaUDQ8TruNB159LXfu6Kj1UBSKukAJu8ISvOeuLbUegkJRNyhXjEKhUFgMJewKhUJhMZSwKxQKhcVQwq5QKBQWQwm7QqFQWAwl7AqFQmExlLArFAqFxVDCrlAoFBZDSFn9/tNCiCm0+u3roQOYLuNwKkG9j1GNb+PU+xjV+DZOPY5xk5Tyqr1FayLsG0EIsVdKuafW47gS9T5GNb6NU+9jVOPbOI0wxpVQrhiFQqGwGErYFQqFwmI0orB/qtYDWAX1PkY1vo1T72NU49s4jTDGkjScj12hUCgUV6YRLXaFQqFQXIG6FXYhxH1CiONCiFNCiD8o8b5bCPEl/f2nhRDDVRzboBDiUSHEESHEYSHEB0ts8zIhREQI8bz+85Fqja9oDOeEEAf14+8t8b4QQvydfg4PCCFuqeLYrik6N88LIaJCiA8t26bq51AI8VkhxKQQ4lDRa21CiIeFECf1f1tX+Oyv6tucFEL8ahXH95dCiGP6d/h1IUR4hc9e8Xqo4Pj+WAhxqeh7fPUKn73iPV/hMX6paHznhBDPr/DZip/DsiClrLsfwA6cBrYALmA/sGvZNv8D+KT++9uAL1VxfL3ALfrvQeBEifG9DPh2jc/jOaDjCu+/GngIEMAdwNM1/L7H0WJ0a3oOgbuAW4BDRa/9f8Af6L//AfAXJT7XBpzR/23Vf2+t0vjuARz6739RanyruR4qOL4/Bn5vFdfAFe/5So5x2ft/DXykVuewHD/1arHfDpySUp6RUmaALwJvWLbNG4DP6b//J/BKIYSoxuCklGNSyn367zHgKNBfjWOXmTcA/yY1ngLCQojeGozjlcBpKeV6k9bKhpTyJ8DsspeLr7XPAT9f4qP3Ag9LKWellHPAw8B91RiflPL7Usqc/udTwEC5j7taVjh/q2E193xZuNIYdQ15C/BgJY5dLepV2PuBi0V/j3C5cJrb6Bd1BGivyuiK0F1ANwNPl3j7hUKI/UKIh4QQu6s6MA0JfF8I8awQ4v4S76/mPFeDt7HyjVTrcwjQLaUc038fB7pLbFMv5/LdaLOwUlzteqgk79ddRZ9dwZVVL+fvTmBCSnlyhfdreQ5XTb0Ke0MghAgAXwU+JKWMLnt7H5pr4Ubg74FvVHt8wEuklLcAPwf8phDirhqM4YoIIVzA64GvlHi7Hs7hEqQ2H6/LUDIhxANADvjCCpvU6nr4R2ArcBMwhubqqFfezpWt9bq/p6B+hf0SMFj094D+WslthBAOoAWYqcrotGM60UT9C1LKry1/X0oZlVIu6L9/F3AKITqqNT79uJf0fyeBr6NNd4tZzXmuND8H7JNSTix/ox7Ooc6E4aLS/50ssU1Nz6UQ4l3Aa4Ff0h8+l7GK66EiSCknpJR5KWUB+OcVjlvza1HXkTcBX1ppm1qdw7VSr8L+M2C7EGKzbtG9Dfjmsm2+CRiRB28GfrjSBV1udD/cZ4CjUsq/WWGbHsPnL4S4He1cV/PB4xdCBI3f0RbYDi3b7JvAr+jRMXcAkSKXQ7VY0UKq9Tksovha+1Xgv0ps8z3gHiFEq+5quEd/reIIIe4Dfh94vZQyscI2q7keKjW+4nWbN65w3NXc85XmbuCYlHKk1Ju1PIdrptartyv9oEVsnEBbKX9Af+1P0S5eAA/a9P0U8AywpYpjewnadPwA8Lz+82rgfcD79G3eDxxGW91/CnhRlc/fFv3Y+/VxGOeweIwC+IR+jg8Ce6o8Rj+aULcUvVbTc4j2kBkDsmh+3l9DW7v5AXASeARo07fdA3y66LPv1q/HU8D/U8XxnULzTxvXohEt1gd890rXQ5XG93n9+jqAJta9y8en/33ZPV+tMeqv/6tx7RVtW/VzWI4flXmqUCgUFqNeXTEKhUKhWCdK2BUKhcJiKGFXKBQKi6GEXaFQKCyGEnaFQqGwGErYFQqFwmIoYVcoFAqLoYRdoVAoLMb/D3VwaYrFul8pAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fe879d09e50>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "def log_normal(mean, log_err):\n",
    "    return np.exp(np.log(mean) + randn()*log_err)\n",
    "\n",
    "class Simulation:\n",
    "    def __init__(self):\n",
    "        self.amplitude = log_normal(5.0, 0.5)\n",
    "        self.phase = np.random.rand() * np.pi * 2\n",
    "        self.dc_shift = np.random.randn()*2\n",
    "        self.error = 0.3\n",
    "    \n",
    "    def simulate_measurement(self, n_periods = 3, n_samples = 200):\n",
    "        t = np.linspace(0, 2*np.pi*n_periods, num=n_samples)\n",
    "        noise = np.random.randn(*t.shape)*self.error\n",
    "        return (t, np.sin(t + self.phase)*self.amplitude + self.dc_shift + noise)\n",
    "    \n",
    "    def __str__(self):\n",
    "        return \"\"\"\n",
    "        Amplitude: %g\n",
    "        Phase shift: %g\n",
    "        \"\"\" % (self.amplitude, self.phase)\n",
    "\n",
    "simulation = Simulation()\n",
    "print(simulation)\n",
    "t, v = simulation.simulate_measurement()\n",
    "plt.plot(t, v);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "regrsssion with dc_component\n",
      "(2.7395462833861544, 2.6564203019480406)\n",
      "regrsssion without dc_component\n",
      "(2.740523808052972, 2.6557441598346117)\n",
      "fast dot product (vectorized)\n",
      "(2.732820345831146, 2.6516032891165473)\n",
      "fast dot product (loop)\n",
      "(2.732820345831147, 2.651603289116549)\n",
      "actual\n",
      "\n",
      "        Amplitude: 2.69926\n",
      "        Phase shift: 2.66615\n",
      "        \n"
     ]
    }
   ],
   "source": [
    "def sin_and_cos_coeff_to_amplitude_and_shift(sin_coeff, cos_coeff):\n",
    "    amplitude = np.sqrt(sin_coeff**2 + cos_coeff**2)\n",
    "    phase_shift = (np.arctan2(cos_coeff, sin_coeff) + 2*np.pi) % (2*np.pi)\n",
    "    return amplitude, phase_shift\n",
    "\n",
    "def reconstruction_regression(t, v, dc_component=True):\n",
    "    sine = np.sin(t)\n",
    "    cosine = np.cos(t)\n",
    "    \n",
    "    regressors = [sine[:,np.newaxis], cosine[:,np.newaxis]]\n",
    "    if dc_component:\n",
    "        regressors.append(np.ones((len(t),1)))\n",
    "    \n",
    "    # solve least squares problem / linear regression\n",
    "    M = np.hstack(regressors)\n",
    "    inv_MTM = np.linalg.inv(np.dot(M.T, M))\n",
    "    coeffs = np.dot(inv_MTM, np.dot(M.T, v))\n",
    "    \n",
    "    sin_coeff, cos_coeff = coeffs[:2]\n",
    "    \n",
    "    return sin_and_cos_coeff_to_amplitude_and_shift(sin_coeff, cos_coeff)\n",
    "\n",
    "def reconstruction_fast_vectorized(t, v):\n",
    "    sine = np.sin(t)\n",
    "    cosine = np.cos(t)\n",
    "    \n",
    "    # = 1 / np.sum(sine**2)\n",
    "    norm_constant = 1.0 / (len(t)/2) \n",
    "    sin_coeff = np.sum(sine*v) * norm_constant\n",
    "    cos_coeff = np.sum(cosine*v) * norm_constant\n",
    "    \n",
    "    return sin_and_cos_coeff_to_amplitude_and_shift(sin_coeff, cos_coeff)\n",
    "\n",
    "def reconstruction_fast_loop(t, v):\n",
    "    n = len(t)\n",
    "    norm_constant = 1.0 / (n/2)\n",
    "    sin_coeff = 0\n",
    "    cos_coeff = 0\n",
    "    for i in range(n):\n",
    "        sin_coeff += v[i]*np.sin(t[i])\n",
    "        cos_coeff += v[i]*np.cos(t[i])\n",
    "    \n",
    "    sin_coeff *= norm_constant\n",
    "    cos_coeff *= norm_constant\n",
    "    return sin_and_cos_coeff_to_amplitude_and_shift(sin_coeff, cos_coeff)\n",
    "\n",
    "print('regrsssion with dc_component')\n",
    "print(reconstruction_regression(t, v, dc_component=True))\n",
    "\n",
    "print('regrsssion without dc_component')\n",
    "print(reconstruction_regression(t, v, dc_component=False))\n",
    "\n",
    "print('fast dot product (vectorized)')\n",
    "print(reconstruction_fast_vectorized(t, v))\n",
    "\n",
    "print('fast dot product (loop)')\n",
    "print(reconstruction_fast_loop(t, v))\n",
    "\n",
    "print('actual')\n",
    "print(simulation)\n",
    "    "
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}