{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Simulation of a coherent optical homodyne (intradyne) transmission using polarization division multiplexed signals and data-aided equalization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "INFO:skcomm:\n",
      "============\n",
      " SNR: 10.0\n",
      "INFO:skcomm:header lengths:\n",
      "              header length: 3072 symbols, \n",
      "              overall frame length: 32768 symbols, \n",
      "              overhead: 9.38%\n",
      "INFO:skcomm:header lengths:\n",
      "              header length: 3072 symbols, \n",
      "              overall frame length: 32768 symbols, \n",
      "              overhead: 9.38%\n",
      "INFO:skcomm.rx:frame number: 0\n",
      "INFO:skcomm.rx:Frame 0: est. shift: [150 200] samples\n",
      "INFO:skcomm.rx:Frame 0: est. CFO (coarse): [301.53956657 301.3172634 ] MHz, mean coarse CFO: 301.428 MHz\n",
      "INFO:skcomm.rx:Frame 0: est. CFO (fine): [-1.03118483 -0.63198944] MHz, mean fine CFO: -0.832 MHz\n",
      "INFO:skcomm.rx:Frame 0: est. CFO (coarse+fine): [300.50838174 300.68527396] MHz, mean fine+coarse CFO: 300.597 MHz\n",
      "INFO:skcomm.rx:Frame 0: estimated SNR (per Rx aperture): [9.96353123 9.84739076] dB\n",
      "INFO:skcomm.rx:frame number: 1\n",
      "INFO:skcomm.rx:Frame 1: est. shift: [0 0] samples\n",
      "INFO:skcomm.rx:Frame 1: est. CFO (coarse): [301.62661329 300.72609994] MHz, mean coarse CFO: 301.176 MHz\n",
      "INFO:skcomm.rx:Frame 1: est. CFO (fine): [-1.21501075 -1.10536297] MHz, mean fine CFO: -1.160 MHz\n",
      "INFO:skcomm.rx:Frame 1: est. CFO (coarse+fine): [300.41160254 299.62073696] MHz, mean fine+coarse CFO: 300.016 MHz\n",
      "INFO:skcomm.rx:Frame 1: estimated SNR (per Rx aperture): [10.07732052  9.93294119] dB\n",
      "INFO:skcomm.rx:frame number: 2\n",
      "INFO:skcomm.rx:Frame 2: est. shift: [0 0] samples\n",
      "INFO:skcomm.rx:Frame 2: est. CFO (coarse): [300.89932729 302.05770271] MHz, mean coarse CFO: 301.479 MHz\n",
      "INFO:skcomm.rx:Frame 2: est. CFO (fine): [-0.51684178 -0.89005712] MHz, mean fine CFO: -0.703 MHz\n",
      "INFO:skcomm.rx:Frame 2: est. CFO (coarse+fine): [300.38248552 301.1676456 ] MHz, mean fine+coarse CFO: 300.775 MHz\n",
      "INFO:skcomm.rx:Frame 2: estimated SNR (per Rx aperture): [9.90358006 9.8422421 ] dB\n",
      "INFO:skcomm.rx:frame number: 3\n",
      "INFO:skcomm.rx:Frame 3: est. shift: [ 0 -1] samples\n",
      "INFO:skcomm.rx:Frame 3: est. CFO (coarse): [297.55038354 297.98818076] MHz, mean coarse CFO: 297.769 MHz\n",
      "INFO:skcomm.rx:Frame 3: est. CFO (fine): [1.24312414 1.30400562] MHz, mean fine CFO: 1.274 MHz\n",
      "INFO:skcomm.rx:Frame 3: est. CFO (coarse+fine): [298.79350768 299.29218638] MHz, mean fine+coarse CFO: 299.043 MHz\n",
      "INFO:skcomm.rx:Frame 3: estimated SNR (per Rx aperture): [9.98468497 9.96824584] dB\n",
      "INFO:skcomm.rx:frame number: 4\n",
      "INFO:skcomm.rx:Frame 4: est. shift: [-1  1] samples\n",
      "INFO:skcomm:EVM X: 31.75%, EVM Y: 31.69%\n",
      "INFO:skcomm:BER X = 1.15e-03 (201 error(s)), BER Y = 1.31e-03 (228 error(s)), mean BER = 1.23e-03 (429 error(s))\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(0.0, 13.540131656211258)"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAi8AAAGlCAYAAAAyFxZnAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/H5lhTAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA90UlEQVR4nO3de3RU5b3/8c/kQkKABCGSkJAIVNEikhRCUrReUC5GG0RrbaXViEovhi5tVrWw2ooctbRySiOSHrSKqVaEnz0HtNJq0tQQtVgsGI+KgByRopDQgGRyIckws39/pJmQG8mEyZPZm/drrSw7e/bMPPNplC/7+T77cVmWZQkAAMAmwgZ6AAAAAIGgeAEAALZC8QIAAGyF4gUAANgKxQsAALAVihcAAGArFC8AAMBWKF4AAICtULwAAABboXgBAAC2QvECAABsJSSLl5dfflnnn3++zjvvPD355JMDPRwAABBCXKG2MeOJEyc0ceJEvfbaa4qLi9PUqVP1t7/9TSNHjhzooQEAgBAQcldetm3bpgsvvFDJyckaOnSosrOzVVxcPNDDAgAAISLoxUt5eblycnKUlJQkl8ulTZs2dTqnsLBQY8eOVXR0tLKysrRt2zb/cwcPHlRycrL/cXJysj777LNgDxMAANhU0IuX+vp6paWlqbCwsMvnN2zYoPz8fC1dulQ7duxQWlqa5syZo8OHDwd7KAAAwIEigv2G2dnZys7O7vb5lStXauHChVqwYIEkac2aNdq8ebPWrl2rxYsXKykpqd2Vls8++0yZmZndvl9TU5Oampr8j30+n44ePaqRI0fK5XIF4RsBAID+ZlmWamtrlZSUpLCwHq6tWP1IkrVx40b/46amJis8PLzdMcuyrFtvvdWaO3euZVmW5fF4rHPPPdf69NNPrdraWmvChAlWdXV1t5+xdOlSSxI//PDDDz/88OOAnwMHDvRYXwT9ysupVFdXy+v1KiEhod3xhIQE7dq1S5IUERGhX/3qV5oxY4Z8Pp/uu+++U640WrJkifLz8/2Pa2pqlJqaqj179mjYsGH980Vs5MSJE3r99dd16aWXKiLC6P/dZxRyNoOczSBnc8i6TW1trSZMmNCrP7tDMqm5c+dq7ty5vTo3KipKUVFRKiwsVGFhobxeryTpH//4h2JiYvpzmLYRExOj7du3D/QwHI+czSBnM8jZHLJu0dDQIEm9avkwWrzEx8crPDxcVVVV7Y5XVVUpMTHxtN47Ly9PeXl5crvdiouL0+zZsxUbG3ta7+kEHo9HJSUlmjVrliIjIwd6OI5FzmaQsxnkbA5Zt3G73b0+12jxMmjQIE2dOlWlpaWaN2+epJYG29LSUi1atCionxUZGXnG/yKcjDzMIGczyNkMcjaHrBXQ9w968VJXV6e9e/f6H+/bt08VFRUaMWKEUlNTlZ+fr9zcXGVkZCgzM1MFBQWqr6/3rz4KFo/HI4/HE9T3tKPWDMiif5GzGeRsBjmbQ9ZtAskg6NsDlJWVacaMGZ2O5+bmqqioSJK0evVqrVixQpWVlUpPT9eqVauUlZV1Wp97cs/Lnj17tG7dOnpeAACnFBYW1vOyXASFz+eTz+fr9vmGhgbNnz9fNTU1PbZ9hNzeRqerteelurqanhcxn2oKOZtBzmacCTl7PB5VVVXp+PHjAzoOy7LU2Nio6OjoM+LeZIMHD1ZCQkKXv1dut1vx8fG9Kl5CcrURAAD9xefz6Z///KfCw8OVlJSkyMjIASscLMtSfX29hgwZ4ujixbIseTwe/etf/9I///lPjRs37rSueDnmygvTRgCA3oiIiFBiYqLGjBmjqKiogR7OGaWpqUmffvqpDh065L+1SSumjZg28jsTLv+GAnI2g5zNcHrOjY2NOnDggH+D4IFk/fuW+MOGDXP0lZdWjY2N+uSTT5SSktIpe6aNxLKzjsjDDHI2g5zNcGrOXq9XLpcrJJp1WxtYW8fjdGFhYXK5XF3+bg3oUulQwVLpFizDM4OczSBnM5yes8fjkWVZPa5+MaF18qN1PE7n8/n8/S/h4eHtnhvQpdIDhZ4XAEBvtPa8pKSkaNCgQX1+H6/P0o4DblXXNyt+yCBNSYlVeFj/Tf3cddddev755/2PzzrrLE2ZMkUPPPCAJk2a1O54V5588kl97Wtf0xtvvKGcnBz/8ZEjR+pLX/qSHnjgAV144YXdfv5TTz2l//iP/9Cbb76pMWPG+I/fd999+utf/6ry8vIe/9xtbm7WgQMHVFlZqRMnTrR7jp4Xel78nD53HSrI2QxyNsPpOQej5+WV9yv1Hy9/qEp3o/9YYmy07v/qF3X1pN5vdxNIz8uCBQtUVVWltWvXSpIqKyv1s5/9TO+9954++eQT/3nh4eF66qmndPXVV7d7/fDhwxUdHa2ysjJdddVV+vDDDxUbG6uDBw/qxz/+sT766CPt2bOn24LOsiz/e7766quSpNLSUmVnZ6usrEwXX3xxj9+XnpceOHWutq/IwwxyNoOczXBqzqfb8/LK+4eUt+4ddfybf5W7UXnr3tF/fXuKrp40ulfvFUjPi8vlUnR0tJKSkiRJSUlJWrJkiS699FIdOXJEZ599tv/cESNG+M/rqPVzEhMTNXz4cCUlJemee+7R3LlztWfPHk2ePLnbMaxdu1aTJk3SE088ofnz5+vOO+9Ufn6+vvKVr/Tq+war58X53UEAAASJ12dp2R93dipcJPmPLfvjTnl9/T+pUVdXp9///vc699xzNXLkyD6/T01NjdavXy9JPU6jpaSkqKCgQPfee6++/e1va+jQoXrwwQf7/Nl95dgrLzTstnB6412oIGczyNkMp+d8Og27f//4iA7VNHb7vCXpUE2j/v5xtb48vueCIpCGXcuy9PLLL2vo0KGSpPr6eo0ePVovvfSSJLV7/c0339ypIfb9999Xamqq/7zWvpX6+npJUk5OjiZMmNDjOHJzc/XEE0/oj3/8o7Zu3arIyMhe5xishl3HFC8nN+xKUnFxMQ27JykpKRnoIZwRyNkMcjbDqTm3NuzW1dWpubk5oNfuP3ys1+dNjO/9NEhtbW2P53g8Hl166aX61a9+JUk6duyYnnrqKV1zzTX6y1/+otTUVP+5Dz/8sK644op2rx86dKjcbrcaGhokSX/60580ePBg/eMf/9DKlSv1yCOPyO129ziO9957Tzt27FBMTIz+8pe/6IILLuj192xubtbx48dVXl7eZcNubzmmeMnLy1NeXp6/YXf27Nk07Mr5jXehgpzNIGcznJ5za8Pu0KFDA27YPWdU764OnDNqeK/+DAqkYTcyMlKxsbFKT0/3H7v00kt11llnacOGDe2mb8aOHdvuvJO1/sV+0qRJGj58uKZOnara2lp95zvfUVlZ2SnH0NzcrEWLFmn+/Pm67LLLdNddd+lrX/uazj///B6/q9SS/eDBg3XZZZd12bDbW44pXjpyaqNZX5GHGeRsBjmb4dScT6dhN2t8vEbHRauyprHLvheXpMS4aGWNj1dYL5ZNB9qw29V5YWFhamxsbHf8VN+t9fjJ5yxatEi/+MUv9OKLL+r666/vdgwPPfSQjh49qoKCAsXFxWnjxo2644479MYbb/QqSxp2AQAwLDzMpaU5EyW1FCona328NGdiv93vpampSZWVlaqsrNSHH36oH/zgB6qrq2t33xapZUqp9bzWn9belq7ExMRo4cKFWrp0qbq7g8rbb7+tX/7yl3rqqacUFxcnSXr88ce1e/du/frXvw7el+wFihcAAAJw9aTR+q9vT1FiXPtpj8S46ICWSffFK6+8otGjR2v06NHKysrS22+/rRdeeKFTf8uCBQv857X+PPbYY6d870WLFunDDz/UCy+80Om5pqYm5ebmasGCBZo9e7b/eOv7/vSnP9Xu3buD8h17w7HTRqw2auH0VQOhgpzNIGcznJ5zMLYHmD0xQVddMEpvf3JUh2ubNGpYlKaNHaHwMFdA7xnIaqO1a9f6b1DX0cmv7bhbc8fzLrvsMv85J78uOTlZTU1NnY5LLVM677//fpfPffOb39Q3v/nNLp/r6vPZHuAkbA8AAOiNYG0PgMCxPUA32B6gPaevGggV5GwGOZvh9JyDsT1AsASy2sgJ2B6gB07tku8r8jCDnM0gZzOcmvPpbg8QTIGsNnICVhsBAIAzEsULAACwFYoXAABgKxQvAADAVhzbsMt9Xlo4/X4NoYKczSBnM5yeczDu8xIsgdznxQm4z0sH3OcFANAb3Odl4HCfl25wn5f2nH6/hlBBzmaQsxlOz5n7vAwc7vPSA6fen6CvyMMMcjaDnM1was5Bu8+Lzyvt/5tUVyUNTZDOuVgKC+/5dSe/RQD3ebntttv0u9/9zv94xIgRmjZtmh555BFNnjzZf7y7Iuj555/XN7/5TZWVlWnGjBn+4/Hx8Zo2bZp++ctf6qKLLurytc8++6y+973v6d1339W5557rP37w4EFdeOGFevDBB7Vo0aIevy/3eQEAYKDsfEkqmCT97qvSf9/R8s+CSS3H+9HVV1+tQ4cO6dChQyotLVVERIS++tWvdjrv6aef9p/X+jNv3rx25+zevVuHDh3Sq6++qqamJl177bVqbm7u8nNvueUWzZkzR7fddlu73pyFCxdq6tSpysvLC+r37AnFCwAAgdj5kvT/bpXcB9sfdx9qOd6PBUxUVJQSExOVmJio9PR0LV68WAcOHNC//vWvducNHz7cf17rT8dpmlGjRikxMVFTpkzRPffcowMHDmjXrl3dfvbjjz+uPXv2aOXKlZKkoqIivfnmm3r66aeNT3k5dtoIAICg83mlV34sqat2UUuSS3plsXTBtQFPIQWqrq5Ov//973Xuuedq5MiRfX6fmpoarV+/XpJO2cB89tln64knntDNN9+stLQ0/fCHP9Sjjz6qlJSUPn92X1G8AADQW/v/1vmKSzuW5P6s5bxxlwb9419++WUNHTpUklRfX6/Ro0fr5Zdf7tQvc/PNN3dairxz506lpqb6H48ZM8b/PpI0d+5cXXDBBaf8/Hnz5ummm27S1VdfrZycHOXm5p72d+oLpo0AAOituqrgnhegGTNmqKKiQhUVFdq2bZvmzJmj7Oxs7d+/v915v/71r/3ntf4kJSW1O+f111/X9u3bVVRUpAkTJmjNmjW9GsPPfvYz+Xw+/fSnPw3a9woUV14AAOitoQnBPS9AQ4YMabfa58knn1RcXJx++9vf6qGHHvIfT0xMbHdeV8aNG6fhw4fr/PPP1+HDh/WNb3xD5eXlPY4hIiKi3T8HAldeAADorXMulmKTJHXXoOqSYpNbzjOgdYn18ePHT+t98vLy9P7772vjxo1BGln/cuyVF7YHaOH023yHCnI2g5zNcHrOp7c9gEua8wu5XsiV5JLrpMZd698FjTVnect5vXjvQLYHsCxLjY2NOniwpefm888/V2Fhoerq6nTttde2e/3Ro0f957UaNmyYhgwZ4j/v5O8fHR2tO++8U0uXLtXcuXNPuXqoq9f3VrC2B3BM8XLy9gCSVFxczPYAJykpKRnoIZwRyNkMcjbDqTm3bg9QV1fX7X1NTin5ckV+9b80uGyZXHWH/IetoYk6fsVSeZIvl9zugN6ytra2x3M8Ho9effVVJScnS2opRs477zwVFRVpypQpcp/0mXfccUen199///364Q9/qIaGBv9nntzoe+utt+rXv/61nnnmGV1//fXdjqOurk5SS6OvO8Dv2dzcrOPHj6u8vLzL7QF6i+0BHM7pt/kOFeRsBjmb4fScg7Y9gM8r/XOrVFcpDU2UUqcHvDya7QHasD2AnHtb674iDzPI2QxyNsOpOQdte4CwMGn8Zac1lkC2B3ACtgcAAABnJIoXAABgKxQvAADAViheAACArVC8AAAAW6F4AQAAtkLxAgAAbIXiBQAA2ArFCwAANlVWViaXy6Vjx44N9FCMCsni5frrr9dZZ52lG2+8caCHAgBAyLjiiit0zz33DPQwBlxIFi933323nnnmmYEeBgAA6EKfNrQMopAsXq644goNGzZsoIcBAEDIuO2227RlyxY9+uijcrlccrlc+uSTTyRJ27dvV0ZGhmJiYnTxxRdr9+7d7V774osvasqUKYqOjtb48eO1bNmydrs6//Of/9R1112noUOHKjY2VjfddJOqqqr8zz/wwANKT0/Xk08+qXHjxik6OlrPPPOMRo4cqaampnafNW/ePN1yyy39F4T6ULyUl5crJydHSUlJcrlc2rRpU6dzCgsL/bt1ZmVladu2bcEYKwAA/cKyLDV4Goz/WJbV6zE++uijmj59uhYuXKhDhw7p0KFDSklJkST95Cc/0a9+9Sv94x//UEREhG6//Xb/615//XXdeuutuvvuu7Vz5049/vjjKioq0sMPPyypZXPI6667TkePHtWWLVtUUlKijz/+WN/4xjfaff7evXv13//93/qf//kfVVRU6Otf/7q8Xq9eeukl/zmHDx/W5s2b231+fwh4V+n6+nqlpaXp9ttv1w033NDp+Q0bNig/P19r1qxRVlaWCgoKNGfOHO3evVujRo2SJKWnp7er+FoVFxcrKSmpD18DAIC+O37iuLLWZRn/3K3f3Nrrc+Pi4jRo0CDFxMQoMTFRkrRr1y5J0sMPP6zLL79ckrR48WJde+21amxsVHR0tJYtW6bFixcrNzdXkjR+/Hg9+OCDuu+++7R06VKVlpbqvffe0759+/zF0DPPPKMLL7xQb7/9tqZNmyapZaromWee0dlnn+0f0/z58/X000/r61//uiTp97//vVJTU3XFFVecXjA9CLh4yc7OVnZ2drfPr1y5UgsXLtSCBQskSWvWrNHmzZu1du1aLV68WJJUUVHRt9F2oampqd0lK7fbLUnyeDzyeDxB+xy7as2ALPoXOZtBzmY4PWePxyPLsuTz+eTz+STJ/8+B0jqeQM9t/eekSZP8/zshIUGSVFlZqdTUVL377rt68803/VdaJMnr9aqxsVF1dXXauXOnUlJSlJyc7H+PCy64QMOHD9cHH3ygqVOnyrIsnXPOORo5cmS7cd5xxx3KysrSgQMHlJycrKKiIuXm5sqyrC6vKvl8PlmWJY/Ho/Dw8HbPBfL7FnDxcirNzc3avn27lixZ4j8WFhammTNnauvW3leXgVi+fLmWLVvW6XhxcbFiYmL65TPtqKSkZKCHcEYgZzPI2Qyn5hwREaHExETV1dX5G08ty1LxtcXGx+Jp8Mjlcqm2trZX5584cULNzc3+v6g3NDRIkhobGzsdc7vdcrvdqqur0+LFi5WTk9Pp/Zqbm9XY2Cifz+d/fSvLsvzv29TUpOjo6E7nfOELX9CkSZP029/+VldeeaU++OADrVu3rtN5J3/e8ePHVV5e3mkGpnXcvRHU4qW6ulper9df9bVKSEjwX9rqjZkzZ+rdd99VfX29xowZoxdeeEHTp0/v8twlS5YoPz/f/9jtdislJUWzZ89WbGxs376Ig3g8HpWUlGjWrFmKjIwc6OE4FjmbQc5mOD3nxsZGHThwQEOHDlV0dLT/eJzijI/FsizV1tZq2LBhcrlcPZ4/ePBghYeH+/98a/1L+rBhw/zHhgwZIkn+5tspU6Zo//79Sk9P7/I909PT9dlnn6mmpsY/bbRz507V1NRoypQpio2NVVRUVLvPPdnChQu1atUqHTlyRFdddZUmTpzY7fgbGxs1ePBgXXbZZe2yl9RtwdOVoBYvwfKXv/yl1+dGRUUpKiqqH0cDAEBoGDt2rLZt26ZPPvlEQ4cO7dVU009/+lPNnTtXqamp+trXvqawsDC9++67+uCDD/Tggw9q5syZuuiii3TLLbdo5cqVOnHihBYtWqTLL79cGRkZPb7//Pnzdd999+nJJ59UUVFREL5lz4JavMTHxys8PLzd8ipJqqqq8jcX9ZfCwkIVFhbK6/VKYtqoI6de/g015GwGOZvh1Jy7mjYaaL2dNvrud7+ru+66S5MmTdLx48dVWFjof31YWMsC4vr6eklSXV2d3G63pk+frvXr1+uRRx7RI488ooiICE2YMEG33HKL/2rHM888ox//+Me64oorFBYWpquuukq//OUv/c83NTXJ6/V2eXXE5XIpJydHxcXFuvLKK095BSVY00YuK5B1Wl0MeOPGjZo3b57/WFZWljIzM/XYY49JamnOSU1N1aJFi/wNu/3J7XYrLi5O1dXVTBvJ+Zd/QwU5m0HOZjg959Zpo9ZbegykQKeNQtWsWbM0ceJEPfroo6c8r7GxUZ988olSUlK6nDaKj49XTU1Nj39+B3zlpa6uTnv37vU/3rdvnyoqKjRixAilpqYqPz9fubm5ysjIUGZmpgoKClRfX+9ffQQAAJzh888/V1lZmcrKyrR69WpjnxvwlZeysjLNmDGj0/Hc3Fz/XNfq1au1YsUKVVZWKj09XatWrVJWVv+unz952mjPnj1at24d00YAgE5ap41SUlI0aNCggR6OrU2ePFnHjh3Tvffeqx/84Ac9nt/c3KwDBw6osrKyy2mj+fPn9+rKy2lNG4Uipo3ac/rl31BBzmaQsxlOz5lpo4EzYNNGdhEZGenIf+n6ijzMIGczyNkMp+bs9XrlcrkUFhbmb3IdKK2rhVrH43RhYWFyuVxd/m4F8rvm2OKFO+y2cPqdMkMFOZtBzmY4Peeu7rA7UFonPwK5w66dBesOu46ZNqLnBQDQG609L2PGjOE+YYY1NTXp008/peelI3pe2nP63HWoIGczyNkMp+fs9Xr18ccf6+yzz9bIkSMHdCxnWs/LkSNH9K9//Uvjx4/vdOWFnhc5d662r8jDDHI2g5zNcGrOkZGROuuss1RdXa2wsDDFxMQMWOHg8/nU3NyspqYmR/e8WJalhoYGVVdX66yzzuqyUZqeF9Hz0srpc9ehgpzNIGczzoScR44cKa/X2+mO8Ka1bn4YHR19Rlx5iY2N1ciRI7v83aLnhZ4XAEAvuFyuTtMX6B9er1enKjnoeaHnxc/pc9ehgpzNIGczyNkcsm5Dz4ucO1fbV+RhBjmbQc5mkLM5ZB1Yz4tzu4MAAIAjOfbKCw27Lc6ExrtQQM5mkLMZ5GwOWbehYZeGXQAAbIWGXRp2/WgGM4OczSBnM8jZHLJuQ8OuaH7qiDzMIGczyNkMcjaHrGnYBQAADkbxAgAAbMWx00asNmpBJ7sZ5GwGOZtBzuaQdRtWG7HaCAAAW2G1EauN/OhkN4OczSBnM8jZHLJuw2oj0bndEXmYQc5mkLMZ5GwOWbPaCAAAOBjFCwAAsBWKFwAAYCuO7XlhqXQLluGZQc5mkLMZ5GwOWbdhqTRLpQEAsBWWSrNU2o9leGaQsxnkbAY5m0PWbVgqLZaddUQeZpCzGeRsBjmbQ9YslQYAAA5G8QIAAGyF4gUAANgKxQsAALAVihcAAGArFC8AAMBWKF4AAICtOPY+L2wP0IJbT5tBzmaQsxnkbA5Zt2F7ALYHAADAVtgegO0B/Lj1tBnkbAY5m0HO5pB1G7YHELda7og8zCBnM8jZDHI2h6zZHgAAADgYxQsAALAVihcAAGArFC8AAMBWKF4AAICtULwAAABboXgBAAC2QvECAABsheIFAADYCsULAACwlZArXg4cOKArrrhCEydO1OTJk/XCCy8M9JAAAEAICbm9jSIiIlRQUKD09HRVVlZq6tSpuuaaazRkyJCBHhoAAAgBIVe8jB49WqNHj5YkJSYmKj4+XkePHqV4AQAAkvowbVReXq6cnBwlJSXJ5XJp06ZNnc4pLCzU2LFjFR0draysLG3btq1Pg9u+fbu8Xq9SUlL69HoAAOA8ARcv9fX1SktLU2FhYZfPb9iwQfn5+Vq6dKl27NihtLQ0zZkzR4cPH/afk56erkmTJnX6OXjwoP+co0eP6tZbb9UTTzzRh68FAACcKuBpo+zsbGVnZ3f7/MqVK7Vw4UItWLBAkrRmzRpt3rxZa9eu1eLFiyVJFRUVp/yMpqYmzZs3T4sXL9bFF1/c47lNTU3+x263W5Lk8Xjk8Xh685UcrTUDsuhf5GwGOZtBzuaQdZtAMghqz0tzc7O2b9+uJUuW+I+FhYVp5syZ2rp1a6/ew7Is3Xbbbbryyit1yy239Hj+8uXLtWzZsk7Hi4uLFRMT0/vBO1xJSclAD+GMQM5mkLMZ5GwOWUsNDQ29PjeoxUt1dbW8Xq8SEhLaHU9ISNCuXbt69R5vvvmmNmzYoMmTJ/v7aZ599llddNFFXZ6/ZMkS5efn+x+73W6lpKRo9uzZio2N7dsXcRCPx6OSkhLNmjVLkZGRAz0cxyJnM8jZDHI2h6zbtM6c9EbIrTb6yle+Ip/P1+vzo6KiFBUV1Y8jAgAAoSSoxUt8fLzCw8NVVVXV7nhVVZUSExOD+VGdFBYWqrCwUF6vVxLTRh1xSdIMcjaDnM0gZ3PIegCnjQYNGqSpU6eqtLRU8+bNkyT5fD6VlpZq0aJFwfyoTvLy8pSXlye32624uDimjf6NS5JmkLMZ5GwGOZtD1m36ddqorq5Oe/fu9T/et2+fKioqNGLECKWmpio/P1+5ubnKyMhQZmamCgoKVF9f7199BAAAcDpclmVZgbygrKxMM2bM6HQ8NzdXRUVFkqTVq1drxYoVqqysVHp6ulatWqWsrKygDLg7J08b7dmzR+vWrWPaCAAAm2hoaND8+fNVU1PT48xJwMVLqGudNqqurmbaSFySNIWczSBnM8jZHLJu43a7FR8f36viJeR2lQYAADgVx1x5YdoIAAD7YtqIaSM/LkmaQc5mkLMZ5GwOWbcJZNoo5G5SFyyRkZFn/C/CycjDDHI2g5zNIGdzyFoBfX96XgAAgK049soLu0q3YMdSM8jZDHI2g5zNIes2gWTgmJ4XGnYBALAvGnZp2PWjGcwMcjaDnM0gZ3PIug0Nu6L5qSPyMIOczSBnM8jZHLIOrGHXscULPS8tmE81g5zNIGczyNkcsm5Dzws9LwAA2Ao9L/S8+DGfagY5m0HOZpCzOWTdhp4XMX/YEXmYQc5mkLMZ5GwOWXOTOgAA4GAULwAAwFYcO23EaqMWdLKbQc5mkLMZ5GwOWbdhtRGrjQAAsBVWG7HayI9OdjPI2QxyNoOczSHrNqw2Ep3bHZGHGeRsBjmbQc7mkDWrjQAAgINRvAAAAFuheAEAALZC8QIAAGzFsQ273OelBfcQMIOczSBnM8jZHLJuw31euM8LAAC2wn1euM+LH/cQMIOczSBnM8jZHLJuw31exJr5jsjDDHI2g5zNIGdzyJr7vAAAAAejeAEAALZC8QIAAGyF4gUAANgKxQsAALAVihcAAGArFC8AAMBWHHufF7YHaMGtp80gZzPI2QxyNoes27A9ANsDAABgK2wPwPYAftx62gxyNoOczSBnc8i6DdsDiFstd0QeZpCzGeRsBjmbQ9ZsDwAAAByM4gUAANgKxQsAALAVihcA9uDzyrX/DSUf3SrX/jckn3egRwRggDi2YReAg+x8SXrlx4pwH1SGJO3/Lyk2Sbr6l9LEuQM9OgCGceUFQGjb+ZL0/26V3AfbH3cfajm+86WBGReAAUPxAiB0+bzSKz+W1NXtqP597JXFTCEBZxiKFwCha//fOl9xaceS3J+1nAfgjEHxAiB01VUF9zwAjhByxcuxY8eUkZGh9PR0TZo0Sb/97W8HekgABsrQhOCeB8ARQm610bBhw1ReXq6YmBjV19dr0qRJuuGGGzRy5MiBHhoAw7wp01WtkTrbOqIwV+fnfZZ02DVSZ6dMV7j54QEYICF35SU8PNy/G3RTU5Msy5LD9o4E0Evb9tfo/uZbJLUUKidrfby0+RZt219jeGQABlLAxUt5eblycnKUlJQkl8ulTZs2dTqnsLBQY8eOVXR0tLKysrRt27aAPuPYsWNKS0vTmDFjdO+99yo+Pj7QYQJwgMO1jXrVl6nve+5RpUa0e65SI/V9zz161Zepw7WNAzRCAAMh4Gmj+vp6paWl6fbbb9cNN9zQ6fkNGzYoPz9fa9asUVZWlgoKCjRnzhzt3r1bo0aNkiSlp6frxIkTnV5bXFyspKQkDR8+XO+++66qqqp0ww036MYbb1RCAnPawJlm1LBoSdKrvkyVNGUoM2yXRumYDmu4tvkukO/ff/9qPQ/AmSHg4iU7O1vZ2dndPr9y5UotXLhQCxYskCStWbNGmzdv1tq1a7V48WJJUkVFRa8+KyEhQWlpaXr99dd14403dnlOU1OTmpqa/I/dbrckyePxyOPx9OpznKw1A7LoX+TcP740ZpgSY6NU5W6ST2F6yzex3fMuSYlxUfrSmGFkH0T8PptD1m0CySCoDbvNzc3avn27lixZ4j8WFhammTNnauvWrb16j6qqKsXExGjYsGGqqalReXm5vv/973d7/vLly7Vs2bJOx4uLi/29M5BKSkoGeghnBHIOvmsSXVrrbp3hPrlr15IlKTuhQa++8ucBGJnz8ftsDllLDQ0NvT43qMVLdXW1vF5vpymehIQE7dq1q1fvsX//fn3nO9/xN+r+4Ac/0EUXXdTt+UuWLFF+fr7/sdvtVkpKimbPnq3Y2Ni+fREH8Xg8Kikp0axZsxQZGTnQw3Escu4/10ia8kGVHvrTLlW6266yjo6L1k+yL9CcC5lSDjZ+n80h6zatMye9EXJLpTMzM3s9rSRJUVFRioqK6r8BARhwcy5M0MwvjtJb//cv/XXrdl05faq+/IWzFd7V+mkAjhfU4iU+Pl7h4eGqqmp/t8uqqiolJiYG86M6KSwsVGFhobzelj1OmDZqj0uSZpBz/5saL9V89A+9+tFAj8T5+H02h6wHcNpo0KBBmjp1qkpLSzVv3jxJks/nU2lpqRYtWhTMj+okLy9PeXl5crvdiouLY9ro37gkaQY5m0HOZpCzOWTdpl+njerq6rR3717/43379qmiokIjRoxQamqq8vPzlZubq4yMDGVmZqqgoED19fX+1UcAAACnw2UFePvasrIyzZgxo9Px3NxcFRUVSZJWr16tFStWqLKyUunp6Vq1apWysrKCMuDunDxttGfPHq1bt45pIwAAbKKhoUHz589XTU1NjzMnARcvoa512qi6upppI3FJ0hRyNoOczSBnc8i6jdvtVnx8fK+Kl5BbbRQskZGRZ/wvwsnIwwxyNoOczSBnc8haAX1/xxYv3GG3BXdvNIOczSBnM8jZHLJuE0gGjpk2oucFAAD7oueFnhc/5lPNIGczyNkMcjaHrNvQ8yLmDzsiDzPI2QxyNoOczSFrel4k0fPSivlUM8jZDHI2g5zNIes29LzQ8wIAgK3Q80LPix/zqWaQsxnkbAY5m0PWbeh5EfOHHZGHGeRsBjmbQc7mkHVgPS9h/TgOAACAoHPslRcadlvQDGYGOZtBzmaQszlk3YaGXRp2AQCwFRp2adj1oxnMDHI2g5zNIGdzyLoNDbui+akj8jCDnM0gZzPI2RyypmEXAAA4GMULAACwFcdOG7HaqAWd7GaQsxnkbAY5m0PWbVhtxGojAABshdVGrDbyo5PdDHI2g5zNIGdzyLoNq41E53ZH5GEGOZtBzmaQszlkzWojAADgYBQvAADAViheAACArVC8AAAAW3Fswy73eWnBPQTMIGczyNkMcjaHrNtwnxfu8wIAgK1wnxfu8+LHPQTMIGczyNkMcjaHrNtwnxexZr4j8jCDnM0gZzPI2Ryy5j4vAADAwSheAACArVC8AAAAW6F4AQAAtkLxAgAAbIXiBQAA2ArFCwAAsBXH3ueF7QFacOtpM8jZDHI2g5zNIes2bA/A9gAAANgK2wOwPYAft542g5zNIGczyNkcsm7D9gDiVssdkYcZ5GwGOZtBzuaQNdsDAAAAB6N4AQAAtkLxAgAAbIXiBQAA2ArFCwAAsBWKFwAAYCsULwAAwFYoXgAAgK1QvAAAAFsJ2eKloaFB55xzjn70ox8N9FAAAEAICdni5eGHH9aXv/zlgR4GAAAIMSFZvHz00UfatWuXsrOzB3ooAAAgxARcvJSXlysnJ0dJSUlyuVzatGlTp3MKCws1duxYRUdHKysrS9u2bQvoM370ox9p+fLlgQ4NAACcAQIuXurr65WWlqbCwsIun9+wYYPy8/O1dOlS7dixQ2lpaZozZ44OHz7sPyc9PV2TJk3q9HPw4EG9+OKLmjBhgiZMmND3bwUAABwrItAXZGdnn3I6Z+XKlVq4cKEWLFggSVqzZo02b96stWvXavHixZKkioqKbl//1ltvaf369XrhhRdUV1cnj8ej2NhY3X///V2e39TUpKamJv9jt9stSfJ4PPJ4PIF+PcdpzYAs+hc5m0HOZpCzOWTdJpAMXJZlWX39IJfLpY0bN2revHmSpObmZsXExOgPf/iD/5gk5ebm6tixY3rxxRcDev+ioiK9//77+s///M9uz3nggQe0bNmyTsfXrVunmJiYgD4PAAAMjIaGBs2fP181NTWKjY095bkBX3k5lerqanm9XiUkJLQ7npCQoF27dgXzo/yWLFmi/Px8/2O3262UlBTNnj27xy9/JvB4PCopKdGsWbMUGRk50MNxLHI2g5zNIGdzyLpN68xJbwS1eAm22267rcdzoqKiFBUV1f+DAQAAISGoxUt8fLzCw8NVVVXV7nhVVZUSExOD+VGdFBYWqrCwUF6vV5JUXFzMtNFJSkpKBnoIZwRyNoOczSBnc8i6Zdqot4JavAwaNEhTp05VaWmpv+fF5/OptLRUixYtCuZHdZKXl6e8vDy53W7FxcUxbfRvXJI0g5zNIGczyNkcsm7Tr9NGdXV12rt3r//xvn37VFFRoREjRig1NVX5+fnKzc1VRkaGMjMzVVBQoPr6ev/qIwAAgNMR8GqjsrIyzZgxo9Px3NxcFRUVSZJWr16tFStWqLKyUunp6Vq1apWysrKCMuDunDxttGfPHlYbAQBgI4GsNjqtpdKhqHXaqLq6mmkjcUnSFHI2g5zNIGdzyLqN2+1WfHy8+aXSoSQyMvKM/0U4GXmYQc5mkLMZ5GwOWSug7+/Y4oU77Lbg7o1mkLMZ5GwGOZtD1m2M3WE3lNDzAgCAfdHzQs+LH/OpZpCzGeRsBjmbQ9Zt6HkR84cdkYcZ5GwGOZtBzuaQNT0vkuh5acV8qhnkbAY5m0HO5pB1G3pe6HkBAMBW6Hmh58WP+VQzyNkMcjaDnM0h6zb0vIj5w47IwwxyNoOczSBnc8g6sJ6XsH4cBwAAQNA59soLDbstaAYzg5zNIGczyNkcsm5Dwy4NuwAA2AoNuzTs+tEMZgY5m0HOZpCzOWTdhoZd0fzUEXmYQc5mkLMZ5GwOWdOwCwAAHIziBQAA2Ipjp41YbdSCTnYzyNkMcjaDnM0h6zasNmK1EQAAtsJqI1Yb+dHJbgY5m0HOZpCzOWTdhtVGonO7I/Iwg5zNIGczyNkcsma1EQAAcDCKFwAAYCsULwAAwFYc2/PCUukWLMMzg5zNIGczyNkcsm7DUmmWSgMAYCsslWaptB/L8MwgZzPI2QxyNoes27BUWiw764g8zCBnM8jZDHI2h6xZKg0AAByM4gUAANgKxQsAALAVihcAAGArFC8AAMBWKF4AAICtULwAAABbcex9XtgeoAW3njaDnM0gZzPI2RyybsP2AGwPAACArbA9ANsD+HHraTPI2QxyNoOczSHrNmwPIG613BF5mEHOZpCzGeRsDlmzPQAAAHAwihcAAGArFC8AAMBWKF4AAICtULwAAABboXgBAAC2QvECAABsheIFAADYCsULAACwFYoXAABgKyG5PcDYsWMVGxursLAwnXXWWXrttdcGekgAACBEhGTxIkl/+9vfNHTo0IEeBgAACDFMGwEAAFsJuHgpLy9XTk6OkpKS5HK5tGnTpk7nFBYWauzYsYqOjlZWVpa2bdsW0Ge4XC5dfvnlmjZtmp577rlAhwgAABws4Gmj+vp6paWl6fbbb9cNN9zQ6fkNGzYoPz9fa9asUVZWlgoKCjRnzhzt3r1bo0aNkiSlp6frxIkTnV5bXFyspKQkvfHGG0pOTtahQ4c0c+ZMXXTRRZo8eXIfvh4AAHCagIuX7OxsZWdnd/v8ypUrtXDhQi1YsECStGbNGm3evFlr167V4sWLJUkVFRWn/Izk5GRJ0ujRo3XNNddox44d3RYvTU1Nampq8j92u92SJI/HI4/H0+vv5VStGZBF/yJnM8jZDHI2h6zbBJJBUBt2m5ubtX37di1ZssR/LCwsTDNnztTWrVt79R719fXy+XwaNmyY6urq9Ne//lU33XRTt+cvX75cy5Yt63S8uLhYMTExgX8JhyopKRnoIZwRyNkMcjaDnM0ha6mhoaHX5wa1eKmurpbX61VCQkK74wkJCdq1a1ev3qOqqkrXX3+9JMnr9WrhwoWaNm1at+cvWbJE+fn5/sdut1spKSmaMWOGYmNj+/AtnOXEiRN67bXXNGPGDEVEhOziMtsjZzPI2QxyNoes27TOnPRGyCU1fvx4vfvuu70+PyoqSlFRUSosLFRhYaG8Xq8k6bXXXuPKy0m4V44Z5GwGOZtBzuaQ9QBeeYmPj1d4eLiqqqraHa+qqlJiYmIwP6qTvLw85eXlye12Ky4uTrNnz+bKi1rmEEtKSjRr1ixFRkYO9HAci5zNIGczyNkcsm4zYFdeBg0apKlTp6q0tFTz5s2TJPl8PpWWlmrRokXB/KgeRUZGnvG/CCcjDzPI2QxyNoOczSFrBfT9Ay5e6urqtHfvXv/jffv2qaKiQiNGjFBqaqry8/OVm5urjIwMZWZmqqCgQPX19f7VR6aw2qgFnexmkLMZ5GwGOZtD1m0CycBlWZYVyJuXlZVpxowZnY7n5uaqqKhIkrR69WqtWLFClZWVSk9P16pVq5SVlRXIxwTs5J6XPXv2aN26dfS8AABgEw0NDZo/f75qamp6bPsIuHgJda09L9XV1fS8iPlUU8jZDHI2g5zNIes2brdb8fHxvSpe2NsIAADYimOuvDBtBACAfTFtxLSRH5ckzSBnM8jZDHI2h6zbBDJtFHI3qQsWlp21Rx5mkLMZ5GwGOZtD1v28VNouWCrdgmV4ZpCzGeRsBjmbQ9Zt+nWpdKii5wUAAPui54WeFz/mU80gZzPI2QxyNoes29DzIuYPOyIPM8jZDHI2g5zNIevAel64zwsAALAVx155oWG3Bc1gZpCzGeRsBjmbQ9ZtaNilYRcAAFuhYZeGXT+awcwgZzPI2QxyNoes29CwK5qfOiIPM8jZDHI2g5zNIWsadgEAgINRvAAAAFtx7LQRq41a0MluBjmbQc5mkLM5ZN2G1UasNgIAwFZYbcRqIz862c0gZzPI2QxyNoes27DaSHRud0QeZpCzGeRsBjmbQ9asNgIAAA5G8QIAAGyF4gUAANgKxQsAALAVxzbscp+XFtxDwAxyNoOczSBnc8i6Dfd54T4vAADYCvd54T4vftxDwAxyNoOczSBnc8i6Dfd5EWvmOyIPM8jZDHI2g5zNIWvu8wIAAByM4gUAANgKxQsAALAVihcAAGArFC8AAMBWKF4AAICtULwAAABbcex9XtgeoAW3njaDnM0gZzPI2RyybsP2AGwPAACArbA9ANsD+HHraTPI2QxyNoOczSHrNmwPIG613BF5mEHOZpCzGeRsDlmzPQAAAHAwihcAAGArFC8AAMBWKF4AAICtULwAAABboXgBAAC2QvECAABsheIFAADYCsULAACwlZAsXvbt26cZM2Zo4sSJuuiii1RfXz/QQwIAACEiJLcHuO222/TQQw/p0ksv1dGjRxUVFTXQQwIAACEi5IqXDz74QJGRkbr00kslSSNGjBjgEQEAgFAS8LRReXm5cnJylJSUJJfLpU2bNnU6p7CwUGPHjlV0dLSysrK0bdu2Xr//Rx99pKFDhyonJ0dTpkzRz3/+80CHCAAAHCzgKy/19fVKS0vT7bffrhtuuKHT8xs2bFB+fr7WrFmjrKwsFRQUaM6cOdq9e7dGjRolSUpPT9eJEyc6vba4uFgnTpzQ66+/roqKCo0aNUpXX321pk2bplmzZvXh6wEAAKcJuHjJzs5WdnZ2t8+vXLlSCxcu1IIFCyRJa9as0ebNm7V27VotXrxYklRRUdHt65OTk5WRkaGUlBRJ0jXXXKOKiopui5empiY1NTX5H9fU1EiSjh49Ko/HE9B3cyKPx6OGhgYdOXLkjN9uvT+RsxnkbAY5m0PWbWprayVJlmX1eG5Qe16am5u1fft2LVmyxH8sLCxMM2fO1NatW3v1HtOmTdPhw4f1+eefKy4uTuXl5frud7/b7fnLly/XsmXLOh0fN25c4F8AAAAMqNraWsXFxZ3ynKAWL9XV1fJ6vUpISGh3PCEhQbt27erVe0REROjnP/+5LrvsMlmWpdmzZ+urX/1qt+cvWbJE+fn5/sc+n09Hjx7VyJEj5XK5+vZFHMTtdislJUUHDhxQbGzsQA/HscjZDHI2g5zNIes2lmWptrZWSUlJPZ4bcquNpJ6npk4WFRXVaSn18OHD+2FU9hYbG3vG/4thAjmbQc5mkLM5ZN2ipysurYJ6k7r4+HiFh4erqqqq3fGqqiolJiYG86MAAMAZKqjFy6BBgzR16lSVlpb6j/l8PpWWlmr69OnB/CgAAHCGCnjaqK6uTnv37vU/3rdvnyoqKjRixAilpqYqPz9fubm5ysjIUGZmpgoKClRfX+9ffQSzoqKitHTpUu5S3M/I2QxyNoOczSHrvnFZvVmTdJKysjLNmDGj0/Hc3FwVFRVJklavXq0VK1aosrJS6enpWrVqlbKysoIyYAAAcGYLuHgBAAAYSCG5qzQAAEB3KF4AAICtULwAAABboXixudraWt1zzz0655xzNHjwYF188cV6++23T/mapqYm/eQnP9E555yjqKgojR07VmvXrjU0YnvqS87PPfec0tLSFBMTo9GjR+v222/XkSNHDI049PW0Q71lWbr//vs1evRoDR48WDNnztRHH33U4/uezq72TtUfWS9fvlzTpk3TsGHDNGrUKM2bN0+7d+/ux28R+vrrd7rVL37xC7lcLt1zzz3BHbgNUbzY3J133qmSkhI9++yzeu+99zR79mzNnDlTn332Wbevuemmm1RaWqqnnnpKu3fv1vPPP6/zzz/f4KjtJ9Cc33zzTd16662644479MEHH+iFF17Qtm3btHDhQsMjD12tO9QXFhZ2+fwjjzyiVatWac2aNfr73/+uIUOGaM6cOWpsbOz2PVt3tV+6dKl27NihtLQ0zZkzR4cPH+6vr2EL/ZH1li1blJeXp7feekslJSXyeDyaPXu26uvr++trhLz+yLnV22+/rccff1yTJ08O9rDtyYJtNTQ0WOHh4dbLL7/c7viUKVOsn/zkJ12+5s9//rMVFxdnHTlyxMQQHaEvOa9YscIaP358u2OrVq2ykpOT+22cdibJ2rhxo/+xz+ezEhMTrRUrVviPHTt2zIqKirKef/75bt8nMzPTysvL8z/2er1WUlKStXz58n4Ztx0FK+uODh8+bEmytmzZEszh2lYwc66trbXOO+88q6SkxLr88sutu+++u59GbR9cebGxEydOyOv1Kjo6ut3xwYMH64033ujyNS+99JIyMjL0yCOPKDk5WRMmTNCPfvQjHT9+3MSQbakvOU+fPl0HDhzQn/70J1mWpaqqKv3hD3/QNddcY2LItrdv3z5VVlZq5syZ/mNxcXHKysrqdof61l3tT35NoLvan4n6knVXampqJEkjRowI+hid4HRyzsvL07XXXtvutWe6kNyYEb0zbNgwTZ8+XQ8++KC++MUvKiEhQc8//7y2bt2qc889t8vXfPzxx3rjjTcUHR2tjRs3qrq6WnfddZeOHDmip59+2vA3sIe+5HzJJZfoueee0ze+8Q01NjbqxIkTysnJ6fZyMtqrrKyUpC53qG99rqNg7Gp/JupL1h35fD7dc889uuSSSzRp0qSgj9EJ+prz+vXrtWPHjh577M40XHmxuWeffVaWZSk5OVlRUVFatWqVbr75ZoWFdf1/rc/nk8vl0nPPPafMzExdc801WrlypX73u99x9eUUAs15586duvvuu3X//fdr+/bteuWVV/TJJ5/oe9/7nuGRA/0vLy9P77//vtavXz/QQ3GUAwcO6O6779Zzzz3X6crvmY7ixea+8IUvaMuWLaqrq9OBAwe0bds2eTwejR8/vsvzR48ereTk5Hbbjn/xi1+UZVn69NNPTQ3bdgLNefny5brkkkt07733avLkyZozZ45+85vfaO3atTp06JDh0dtP6y70gexQz672fdOXrE+2aNEivfzyy3rttdc0ZsyYfhmjE/Ql5+3bt+vw4cOaMmWKIiIiFBERoS1btmjVqlWKiIiQ1+vt93GHKooXhxgyZIhGjx6tzz//XK+++qquu+66Ls+75JJLdPDgQdXV1fmP7dmzR2FhYfyHpxd6m3NDQ0OnqzLh4eGSWpZL4tTGjRunxMTEdjvUu91u/f3vf+92h3p2te+bvmQttfweL1q0SBs3btRf//pXjRs3zsRwbasvOV911VV67733VFFR4f/JyMjQt771LVVUVPj/m3JGGshuYZy+V155xfrzn/9sffzxx1ZxcbGVlpZmZWVlWc3NzZZlWdbixYutW265xX9+bW2tNWbMGOvGG2+0PvjgA2vLli3WeeedZ915550D9RVsIdCcn376aSsiIsL6zW9+Y/3f//2f9cYbb1gZGRlWZmbmQH2FkFNbW2u988471jvvvGNJslauXGm988471v79+y3Lsqxf/OIX1vDhw60XX3zR+t///V/ruuuus8aNG2cdP37c/x5XXnml9dhjj/kfr1+/3oqKirKKioqsnTt3Wt/5znes4cOHW5WVlca/Xyjpj6y///3vW3FxcVZZWZl16NAh/09DQ4Px7xcq+iPnjlht1ILixeY2bNhgjR8/3ho0aJCVmJho5eXlWceOHfM/n5uba11++eXtXvPhhx9aM2fOtAYPHmyNGTPGys/PP6P/g9Mbfcl51apV1sSJE63Bgwdbo0ePtr71rW9Zn376qeGRh67XXnvNktTpJzc317KslqWlP/vZz6yEhAQrKirKuuqqq6zdu3e3e49zzjnHWrp0abtjjz32mJWammoNGjTIyszMtN566y1D3yh09UfWXb2fJOvpp58298VCTH/9Tp+M4qUFu0oDAABboecFAADYCsULAACwFYoXAABgKxQvAADAViheAACArVC8AAAAW6F4AQAAtkLxAgAAbIXiBQAA2ArFCwAAsBWKFwAAYCsULwAAwFb+P1I5r0MZriaNAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAh8AAAGdCAYAAACyzRGfAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/H5lhTAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAyK0lEQVR4nO3dd3xUdb7/8fdk0oEEqUkgCaDgBdEgAhILZQ1EVmkWdMGVYrlKUFiuiqiUAIpYuIAILOpSlPbbu1IWr0AEAVF6EQXEIC1gQicB0iaZ8/uDm9GYBBKY+aa9no8HD53v+Z5zPueTPJg3p8zYLMuyBAAAYIhXaRcAAAAqF8IHAAAwivABAACMInwAAACjCB8AAMAowgcAADCK8AEAAIwifAAAAKO8S7uAP3I6nfr1119VrVo12Wy20i4HAAAUg2VZunDhgsLCwuTldeVzG2UufPz6668KDw8v7TIAAMA1SEpKUv369a84p8yFj2rVqkm6XHxQUFApV1P6HA6HVq1apc6dO8vHx6e0y6mw6LMZ9Nkcem0Gff5NWlqawsPDXe/jV1LmwkfepZagoCDChy7/YgcGBiooKKjS/2J7En02gz6bQ6/NoM8FFeeWCW44BQAARhE+AACAUYQPAABgVJm75wMAUPFYlqWcnBzl5uaWdilu5XA45O3trczMzAp3bIXx8fGR3W6/7u0QPgAAHpWdna3k5GSlp6eXdiluZ1mWQkJClJSUVCk+m8pms6l+/fqqWrXqdW2H8AEA8Bin06lDhw7JbrcrLCxMvr6+FepN2ul06uLFi6patepVP1irvLMsS6dOndKxY8fUuHHj6zoDQvgAAHhMdna2nE6nwsPDFRgYWNrluJ3T6VR2drb8/f0rfPiQpNq1a+vw4cNyOBzXFT4qfqcAAKWuMrwxVwbuOmvFbwMAADCK8AEAAIwifAAAyrxcp6WNv5zR0l3HtfGXM8p1WqVdklvNnj1b1atXd9v22rVrp/nz55donRkzZqhr165uq+FKCB8AgDJtxY/JumfCGv3lo00avHCX/vLRJt0zYY1W/Jjs0f0mJSVpwIABrqd0IiMjNXjwYJ05c8Y1Z+nSperZs6eio6N111136dChQ27Z9+zZs2Wz2Vx/qlatqjvuuEOff/75VdddtmyZTpw4occff1zS5W+Lv+GGGzRlypR88zZv3iwfHx+tWrVKkjRgwADt2LFD33zzjVuO4UoIHwCAMmvFj8l6/rMdSk7NzDeekpqp5z/b4bEAcvDgQbVq1UqJiYlasGCBDhw4oBkzZmj16tWKjo7W2bNnJUldunTR4sWLtXHjRjVr1kwrVqxwWw1BQUFKTk5WcnKydu7cqdjYWPXq1Uv79++/4npTpkxR//79XTf5hoWF6YMPPtDw4cOVmJgoScrIyFDfvn319NNPq3PnzpIkX19f9e7du0BI8QTCBwCgTMp1Wor/914VdoElbyz+33s9cgkmLi5Ovr6+WrVqldq3b6+IiAh16dJFX331lY4fP67XX39d0uU3bEn64osvdOzYMfXv379Y2589e7YiIiIUGBionj175jubksdmsykkJEQhISFq3Lixxo0bJy8vL+3evbvI7Z46dUpr1qwpcPnkiSeeUGxsrPr16yen06nhw4fL4XDo3XffzTeva9euWrZsmTIyMop1HNeK8AEAKJO2HDpb4IzH71mSklMzteXQWbfu9+zZs1q5cqUGDhyogICAfMtCQkLUp08fLVq0SJZlyel06r333tOSJUu0ZMkS+fv7X3X7mzdv1lNPPaVBgwZp165d6tixo8aNG3fFdXJzczVnzhxJUsuWLYuct2HDBgUGBqpp06YFls2YMUOJiYnq06ePpk6dqlmzZhX4pNJWrVopJydHmzdvvupxXA8+ZAwAUCadvFB08LiWecWVmJgoy7IKfQOXpKZNm+rcuXM6deqUPvvsM73//vuKiopShw4d1KdPH73wwgtX3P7kyZN1//3365VXXpEkNWnSRN99912BSzapqamucJCRkSEfHx/NnDlTN954Y5HbPnLkiOrWrVvo56rUqVNHY8eO1XPPPafnn39e7dq1KzAnMDBQwcHBOnLkyBWP4XoRPgAAZVKdalc/i1CSeSVlWVe+nOPr66shQ4ZowIABCgoKKvYHqe3bt089e/bMNxYdHV0gfFSrVk07duyQJKWnp+urr77Sc889p5o1axb5VEpGRkaRZ19yc3M1e/ZsBQYGatOmTcrJyZG3d8EYEBAQ4PHv4eGyCwCgTGrTsIZCg/1V1Gdq2iSFBvurTcMabt3vTTfdJJvNpn379hW6fN++fapdu7ZbH40tjJeXl2666SbddNNNuu222zR06FB16NBBEyZMKHKdWrVq6dy5c4Uue++993Tw4EFt27ZNx44d01tvvVXovLNnz6p27dpuOYaiED4AAGWS3cumUV2bSVKBAJL3elTXZrJ7ufeL6mrWrKlOnTpp2rRpBW68TElJ0bx589SvX79r3n7Tpk0L3FOxadOmYq1rt9uveDPo7bffrpSUlAIBZM+ePRo1apSmT5+upk2bavr06Ro3blyBm1d/+eUXZWZm6vbbby/m0VwbwgcAoMy6v3mopj/RUiHB+S8lhAT7a/oTLXV/81CP7Hfq1KnKyspSbGys1q9fr6SkJK1YsUKdOnVSkyZNNHLkyGve9osvvqgVK1bovffeU2JioqZOnVroI7qWZSklJUUpKSk6dOiQZs6cqZUrV6p79+5Fbvv2229XrVq19O2337rGcnJy1LdvXz300EN66KGHJEkPP/ywHn74YfXr1085OTmuud98840aNWp0xftK3IHwAQAo0+5vHqoNw/6kBc+01eTHW2jBM221YdifPBY8JKlx48baunWrGjVqpF69eikyMlJdunRRkyZN9O233xZ4SqQk2rZtq48++kiTJ09WVFSUVq1apTfeeKPAvLS0NIWGhio0NFRNmzbV+++/rzFjxrge8y2M3W5X//79NW/ePNfYW2+9pePHj2vq1Kn55n744YdKTk7Od/llwYIFeuaZZ6752IrNKmNSU1MtSVZqamppl1ImZGdnW0uWLLGys7NLu5QKjT6bQZ/NKSu9zsjIsPbu3WtlZGSUah3uMHLkSKtq1arWxo0bXWO5ubnWuXPnrNzc3FKsLL/k5GSrRo0a1uHDh0u03o8//mjVqVPHOn/+fJFzrvTzLMn7N2c+AAAohvj4eE2ZMkWbNm2S0+ks7XKKFBISok8++URHjx4t0XrJycmaO3eugoODPVTZb3jUFgCAYirOJ5h26dKlyO9Hee211/Taa6+5u6wCevToUeJ1YmJi3F9IEQgfAAC40ccff1zkEyk1arj3seDyivABAIAb1atXr7RLKPO45wMAABhF+AAAAEYRPgAAgFGEDwAAYBThAwAAGEX4AACUfc5c6dA30g//c/m/ztzSrsitZs+e7dZvyW3Xrp3mz59f7PmvvvqqXnjhBbft/2oIHwCAsm3vMmlSc2nOg9K/nrr830nNL497UFJSkgYMGKCwsDD5+voqMjJSgwcP1pkzZ1xzli5dqp49eyo6Olp33XWXDh065JZ9z549WzabzfWnatWquuOOO/T5559fdd1ly5bpxIkTevzxx3X69GmFhITk+/6WPL169VLbtm2Vm5url156SXPmzNHBgwfdUv/VED4AAGXX3mXS/3tSSvs1/3ha8uVxDwWQgwcPqlWrVkpMTNSCBQt04MABzZgxQ6tXr1Z0dLTOnj0r6fKnmS5evFgbN25Us2bNCv122msVFBSk5ORkJScna+fOnYqNjVWvXr20f//+K643ZcoU9e/fX15eXqpVq5Zmzpyp+Ph4/fDDD645//znP7V8+XLNmTNHdrtdtWrVUmxsrKZPn+62+q+E8AEAKJucudKKYZKsQhb+39iKVz1yCSYuLk6+vr5atWqV2rdvr4iICHXp0kVfffWVjh8/7vpmWV9fX0nSF198oWPHjhXr49ely2c2IiIiFBgYqJ49e+Y7m5LHZrMpJCREISEhaty4scaNGycvLy/t3r27yO2eOnVKa9asUdeuXV1j3bp1U+/evdW3b185HA6dOnVKcXFxevvtt3XzzTe75nXt2lULFy4sVv3Xq8ThY/369eratavCwsJks9m0ZMkS1zKHw6Fhw4bp1ltvVZUqVRQWFqYnn3xSv/76a9EbBACgMEe+K3jGIx9LSjt+eZ4bnT17VitXrtTAgQMVEBCQb1lISIj69OmjRYsWybIsOZ1Ovffee1qyZImWLFkif3//q25/8+bNeuqppzRo0CDt2rVLHTt21Lhx4664Tm5urubMmSNJatmyZZHzNmzYoMDAQDVt2jTf+OTJk3XmzBmNHTtWAwcOVPPmzQvc49GmTRsdO3ZMhw8fvuoxXK8Sf7z6pUuXFBUVpQEDBuihhx7Ktyw9PV07duzQiBEjFBUVpXPnzmnw4MHq1q2btm3b5raiAQCVwMUT7p1XTImJibIsq8AbeJ6mTZvq3LlzOnXqlD777DO9//77ioqKUocOHdSnT5+r3rg5efJk3X///XrllVckSU2aNNF3331X4JJNamqqqlatKknKyMiQj4+PZs6cqRtvvLHIbR85ckR169aVl1f+cwtBQUGaNWuWOnfurCpVqmj37t2y2Wz55oSFhbm20aBBgysew/Uqcfjo0qWLunTpUuiy4OBgJSQk5BubOnWq2rRpo6NHjyoiIuLaqgQAVD5V67p3XglZVmGXe37j6+urIUOGaMCAAQoKCirwhl+Uffv2qWfPnvnGoqOjC4SPatWqaceOHZIu/+P+q6++0nPPPaeaNWvmu6zyexkZGUWeffnTn/6ktm3bqkWLFoqMjCywPO8sT3p6erGO43p4/IvlUlNTZbPZinyEKCsrS1lZWa7XaWlpki5fwnE4HJ4ur8zL6wG98Cz6bAZ9Nqes9NrhcLguTzidzpKtHN5WtqAwKS1ZtkLu+7Bkk4LCZIW3lUq67Sto1KiRbDab9u7dq+7duxdYvnfvXtWuXVtBQUGugJJ3jMX1x/l528kbczqd8vLyUqNGjVxzmjdvrpUrV2rChAl64IEHCt1ujRo1dO7cuSJr8fb2lt1uL3T56dOnJUk1a9Yscn2n0ynLsuRwOGS32/MtK8nvmkfDR2ZmpoYNG6a//OUvCgoKKnTO+PHjFR8fX2B81apVCgwM9GR55cofzyjBM+izGfTZnNLutbe3t0JCQnTx4kVlZ2eXeH2fdiMVuPx5WbLlCyCWLl8ySG83Qo6Ll9xWryT5+PioY8eOmjZtmgYMGJDvvo8TJ05o/vz5euqpp1z/WJakCxcuFHv7N954o7799tt863/zzTeyLMs1lpmZme91HsuydPHixQLjeZo0aaKUlBQdPXq00H/05+TkKDs7u9D1t2zZIh8fH4WHhxe5/ezsbGVkZGj9+vXKycnJt6wkZ0w8Fj4cDod69eoly7Ku+OjO8OHDNXToUNfrtLQ0hYeHq3PnzkUGlsrE4XAoISFBnTp1ko+PT2mXU2HRZzPoszllpdeZmZlKSkpS1apVi3UzZgEtH5MVECjbylfz33waFCYrdrwCmnZVQNFrX7Np06bpnnvu0WOPPaYxY8aoYcOG2rNnj4YNG6YmTZpo3Lhxqlq1qizL0oULF1StWrUC91AUZejQobr33nv10UcfqVu3blq1apXWrFkjm83met/L61XeG3pGRoYSEhK0Zs0ajRgxosj3x3vuuUe1atXS7t279eCDDxZY7u3tLV9f30LX37Fjh+69917VrVv0ZazMzEwFBASoXbt2BX6eRQWWwngkfOQFjyNHjmjNmjVXDBF+fn7y8/MrMO7j48NfTr9DP8ygz2bQZ3NKu9e5ubmy2Wzy8vIq9j0RBdzSXWr64OWnWi6ekKrWlS3yLtm87Fdf9xrdfPPN2rp1q0aPHq3HH39cJ0+elGVZeuihh/Tpp5+6zsznXZ7IO8biuOuuu/TRRx9p1KhRGjVqlGJiYvTGG29o7Nixrm14eXkpLS1N9erVk3T5vTIyMlJjxozRsGHDityXl5eX+vfvrwULFqhbt26Fzimq1kWLFmn06NFXPA4vLy/ZbLZCf69K8nvm9vCRFzwSExP19ddfq2bNmu7eBQCgsvGySw3vNbrLBg0aaPbs2a7Xo0aN0sSJE7V79261bdv2urY9YMAADRgwIN/Yf/3Xf7n+v1+/furXr981bftvf/ubbrnlFh05cqTAjaVr164tdJ0vv/xSXl5eeuSRR65pnyVV4vBx8eJFHThwwPX60KFD2rVrl2rUqKHQ0FA98sgj2rFjh5YvX67c3FylpKRIunwTTN6HsQAAUN7Ex8erQYMG2rRpk9q0aXPtZ3I8LCQkRJ988omOHj1a6FMthbl06ZJmzZolb2+PP4ci6RrCx7Zt29SxY0fX67z7Nfr27avRo0dr2bLLH3XbokWLfOt9/fXX6tChw7VXCgBAKSvOJ5h26dJF33zzTaHLXnvtNb322mvuLquAHj16lGi+qTMeeUocPjp06HDFZ5+v9lw0AAAV2ccff6yMjIxCl9WoUcNwNWWTmfMrAABUEnk3iaJoZfOCFQAAqLAIHwAAwCjCBwAAMIrwAQAAjCJ8AAAAowgfAAAUw9q1a2Wz2XT+/PnSLqXcI3wAAFCIDh06aMiQIaVdRoVE+AAAoAzLzs4u7RLcjvABADDGsiylO9JL5U9JPoG7X79+WrdunSZPniybzSabzabDhw9LkrZv365WrVopMDBQ99xzjxITE/Otu3TpUrVs2VL+/v5q1KiR4uPjlZOT41p+9OhRde/eXVWrVlVQUJB69eqlEydOuJaPHj1aLVq00Mcff6yGDRvK399fc+fOVc2aNZWVlZVvXz169NBf//rXa/hJlC4+4RQAYExGTobunH9nqex7c+/NCvQJLNbcyZMn6+eff1bz5s01ZswYSdKePXskSa+//rref/991a5dW88995wGDRqkjRs3SpK++eYbPfnkk5oyZYruvfde/fLLL3r22WclXf5WXKfT6Qoe69atU05OjuLi4vTYY4/l+8bZAwcO6F//+pc+//xz2e12NW7cWC+++KKWLVumRx99VJJ08uRJffHFF1q1apW7WmQM4QMAgD8IDg6Wr6+vAgMDFRISIkn66aefJElvvvmm2rdvL0l65ZVX1LVrV2VmZiowMFDx8fF69dVX1bdvX0lSo0aNNHbsWL3yyisaNWqUVq9erR9++EGHDh1SeHi4JGnu3Lm65ZZbtHXrVrVu3VrS5Ustc+fOVe3atV019e7dW7NmzXKFj88++0wRERHl8ktbCR8AAGMCvAO0uffmUtu3O9x2222u/w8NDZV0+SxEgwYN9P333+vbb7/Vm2++6ZqTm5urzMxMpaena9++fQoPD3cFD0lq1qyZqlevrn379rnCR2RkZL7gIUnPPPOMWrdurePHj6tevXqaPXu2+vXrJ5vN5pbjMonwAQAwxmazFfvSR1nl4+Pj+v+8N36n0ylJunjxouLj4/XQQw8VWM/f37/Y+6hSpUqBsdtvv11RUVGaO3euOnfurD179uiLL74oafllAuEDAIBC+Pr6Kjc3t0TrtGzZUvv379dNN91U6PKmTZsqKSlJSUlJrrMfe/fu1fnz59WsWbOrbv/pp5/WpEmTdPz4ccXExOQ7g1KeED4AAChEgwYNtHnzZh0+fFhVq1Z1nd24kpEjR+rBBx9URESEHnnkEXl5een777/Xjz/+qHHjxikmJka33nqr+vTpo0mTJiknJ0cDBw5U+/bt1apVq6tuv3fv3nrppZf00Ucfae7cue44zFLBo7YAABTipZdekt1uV7NmzVS7dm0dPXr0quvExsZq+fLlWrVqlVq3bq22bdvqv//7vxUZGSnp8mWapUuX6oYbblC7du0UExOjRo0aadGiRcWqKTg4WA8//LCqVq2qHj16XM/hlSrOfAAAUIgmTZq4HqHN069fv3yvW7RooXPnzikoKMg1Fhsbq9jY2CK3GxERoaVLlxa5fPTo0Ro9enSRy48fP64+ffrIz8/vygdQhhE+AAAoB86dO6e1a9dq7dq1mjZtWmmXc10IHwAAlAO33367zp07pwkTJujmm28u7XKuC+EDAIByIO/j3SsCbjgFAABGET4AAB5Xki91Q9nlrp8j4QMA4DF5nwaanp5eypXAHbKzsyVJdrv9urbDPR8AAI+x2+2qXr26Tp48KUkKDAwsl99FUhSn06ns7GxlZmbKy6ti/3ve6XTq1KlTCgwMlLf39cUHwgcAwKPyvhU2L4BUJJZlKSMjQwEBARUqVBXFy8tLERER132shA8AgEfZbDaFhoaqTp06cjgcpV2OWzkcDq1fv17t2rXL94VzFZWvr69bzvAQPgAARtjt9uu+V6CssdvtysnJkb+/f6UIH+5SsS9QAQCAMofwAQAAjCJ8AAAAowgfAADAKMIHAAAwivABAACMInwAAACjCB8AAMAowgcAADCK8AEAAIwifAAAAKMIHwAAwCjCBwAAMIrwAQAAjCJ8AAAAo0ocPtavX6+uXbsqLCxMNptNS5YsybfcsiyNHDlSoaGhCggIUExMjBITE91VLwAAKOdKHD4uXbqkqKgoffjhh4Uuf+eddzRlyhTNmDFDmzdvVpUqVRQbG6vMzMzrLhYAAJR/3iVdoUuXLurSpUuhyyzL0qRJk/TGG2+oe/fukqS5c+eqbt26WrJkiR5//PHrqxYAAJR7JQ4fV3Lo0CGlpKQoJibGNRYcHKw777xTGzduLDR8ZGVlKSsry/U6LS1NkuRwOORwONxZXrmU1wN64Vn02Qz6bA69NoM+/6YkPXBr+EhJSZEk1a1bN9943bp1Xcv+aPz48YqPjy8wvmrVKgUGBrqzvHItISGhtEuoFOizGfTZHHptBn2W0tPTiz3XreHjWgwfPlxDhw51vU5LS1N4eLg6d+6soKCgUqysbHA4HEpISFCnTp3k4+NT2uVUWPTZDPpsDr02gz7/Ju/KRXG4NXyEhIRIkk6cOKHQ0FDX+IkTJ9SiRYtC1/Hz85Ofn1+BcR8fn0r/g/w9+mEGfTaDPptDr82gzyrR8bv1cz4aNmyokJAQrV692jWWlpamzZs3Kzo62p27AgAA5VSJz3xcvHhRBw4ccL0+dOiQdu3apRo1aigiIkJDhgzRuHHj1LhxYzVs2FAjRoxQWFiYevTo4c66AZQzuU5Lmw+d1fbTNtU8dFbRN9WR3ctW2mUBKAUlDh/btm1Tx44dXa/z7tfo27evZs+erVdeeUWXLl3Ss88+q/Pnz+uee+7RihUr5O/v776qAZQrK35MVvy/9yo5NVOSXXMTtyk02F+jujbT/c1Dr7o+gIqlxOGjQ4cOsiyryOU2m01jxozRmDFjrqswABXDih+T9fxnO/THvzVSUjP1/Gc7NP2JlgQQoJLhu10AeEyu01L8v/fKkuQlp9p67VU3r+/U1muvbHJKkuL/vVe5zqL/QQOg4in1R20BVFxbDp1VcmqmYr22aJTPXIXZzrqW/WrVULzjSa1MbaMth84q+saapVgpAJM48wHAY05euBw8pvtMUojO5lsWorOa7jNJsV5bdPIC3/0EVCaEDwAeU6eKj0b5zJUk/fHBlrzXo3w+VZ0qlfvzEYDKhssuADymjf0n2W1ni1zuZZPCdEZ17T9JqmOuMAClijMfADzGfumkW+cBqBgIHwA8p2rdq88pyTwAFQLhA4DnRN4lBYVJKuqTTG1SUL3L8wBUGoQPAJ7jZZfun/B/L/4YQP7v9f1vX54HoNIgfADwrGbdpF5zpaA/fIppUNjl8WbdSqcuAKWGp10AeF6zbtJ/PKCcg+u165uVanFvrLwbteOMB1BJceYDgBledlmR9+h4jWhZkfcQPIBKjPABAACMInwAAACjCB8AAMAowgcAADCK8AEAAIwifAAAAKMIHwAAwCjCBwAAMIrwAQAAjCJ8AAAAowgfAADAKMIHAAAwivABAACMInwAAACjCB8AAMAowgcAADCK8AEAAIwifAAAAKMIHwAAwCjCBwAAMIrwAQAAjCJ8AAAAowgfAADAKMIHAAAwivABAACMInwAAACjCB8AAMAowgcAADCK8AEAAIwifAAAAKMIHwAAwCjCBwAAMMrt4SM3N1cjRoxQw4YNFRAQoBtvvFFjx46VZVnu3hUAACiHvN29wQkTJmj69OmaM2eObrnlFm3btk39+/dXcHCwXnzxRXfvDgAAlDNuDx/fffedunfvrgceeECS1KBBAy1YsEBbtmxx964AAEA55Pbwcdddd2nmzJn6+eef1aRJE33//ffasGGDJk6cWOj8rKwsZWVluV6npaVJkhwOhxwOh7vLK3fyekAvPIs+m0GfzaHXZtDn35SkBzbLzTdjOJ1Ovfbaa3rnnXdkt9uVm5urN998U8OHDy90/ujRoxUfH19gfP78+QoMDHRnaQAAwEPS09PVu3dvpaamKigo6Ipz3R4+Fi5cqJdfflnvvvuubrnlFu3atUtDhgzRxIkT1bdv3wLzCzvzER4ertOnT1+1+MrA4XAoISFBnTp1ko+PT2mXU2HRZzPoszn02gz6/Ju0tDTVqlWrWOHD7ZddXn75Zb366qt6/PHHJUm33nqrjhw5ovHjxxcaPvz8/OTn51dg3MfHp9L/IH+PfphBn82gz+bQazPos0p0/G5/1DY9PV1eXvk3a7fb5XQ63b0rAABQDrn9zEfXrl315ptvKiIiQrfccot27typiRMnasCAAe7eFQAAKIfcHj4++OADjRgxQgMHDtTJkycVFham//zP/9TIkSPdvSsAAFAOuT18VKtWTZMmTdKkSZPcvWkAAFAB8N0uAADAKMIHAAAwivABAACMInwAAACjCB8AAMAowgcAADCK8AEAAIwifAAAAKMIHwAAwCjCBwAAMIrwAQAAjCJ8AAAAowgfAADAKMIHAAAwivABAACMInwAAACjCB8AAMAowgcAADCK8AEAAIwifAAAAKMIHwAAwCjCBwAAMIrwAQAAjCJ8AAAAowgfAADAKMIHAAAwivABAACMInwAAACjCB8AAMAowgcAADCK8AEAAIwifAAAAKMIHwAAwCjCBwAAMIrwAQAAjCJ8AAAAowgfAADAKMIHAAAwivABAACMInwAAACjCB8AAMAowgcAADCK8AEAAIzySPg4fvy4nnjiCdWsWVMBAQG69dZbtW3bNk/sCgAAlDPe7t7guXPndPfdd6tjx4768ssvVbt2bSUmJuqGG25w964AAEA55PbwMWHCBIWHh2vWrFmusYYNG7p7NwAAoJxye/hYtmyZYmNj9eijj2rdunWqV6+eBg4cqGeeeabQ+VlZWcrKynK9TktLkyQ5HA45HA53l1fu5PWAXngWfTaDPptDr82gz78pSQ9slmVZ7ty5v7+/JGno0KF69NFHtXXrVg0ePFgzZsxQ3759C8wfPXq04uPjC4zPnz9fgYGB7iwNAAB4SHp6unr37q3U1FQFBQVdca7bw4evr69atWql7777zjX24osvauvWrdq4cWOB+YWd+QgPD9fp06evWnxl4HA4lJCQoE6dOsnHx6e0y6mw6LMZ9Nkcem0Gff5NWlqaatWqVazw4fbLLqGhoWrWrFm+saZNm+pf//pXofP9/Pzk5+dXYNzHx6fS/yB/j36YQZ/NoM/m0Gsz6LNKdPxuf9T27rvv1v79+/ON/fzzz4qMjHT3rgAAQDnk9vDxt7/9TZs2bdJbb72lAwcOaP78+Zo5c6bi4uLcvSsAAFAOuT18tG7dWosXL9aCBQvUvHlzjR07VpMmTVKfPn3cvSsAAFAOuf2eD0l68MEH9eCDD3pi0wAAoJzju10AAIBRhA8AAGAU4QMAABhF+AAAAEYRPgAAgFGEDwAAYBThAwAAGEX4AAAARhE+AACAUYQPAABgFOEDAAAYRfgAAABGET4AAIBRhA8AAGAU4QMAABhF+AAAAEYRPgAAgFGEDwAAYBThAwAAGEX4AAAARhE+AACAUYQPAABgFOEDAAAYRfgAAABGET4AAIBRhA8AAGAU4QMAABhF+AAAAEYRPgAAgFGEDwAAYBThAwAAGEX4AAAARhE+AACAUYQPAABgFOEDAAAYRfgAAABGET4AAIBRhA8AAGAU4QMAABhF+AAAAEYRPgAAgFGEDwAAYBThAwAAGEX4AAAARnk8fLz99tuy2WwaMmSIp3cFAADKAY+Gj61bt+rvf/+7brvtNk/uBgAAlCMeCx8XL15Unz599NFHH+mGG27w1G4AAEA54+2pDcfFxemBBx5QTEyMxo0bV+S8rKwsZWVluV6npaVJkhwOhxwOh6fKKzfyekAvPIs+m0GfzaHXZtDn35SkBx4JHwsXLtSOHTu0devWq84dP3684uPjC4yvWrVKgYGBniivXEpISCjtEioF+mwGfTaHXptBn6X09PRiz7VZlmW5c+dJSUlq1aqVEhISXPd6dOjQQS1atNCkSZMKzC/szEd4eLhOnz6toKAgd5ZWLjkcDiUkJKhTp07y8fEp7XIqLPpsBn02h16bQZ9/k5aWplq1aik1NfWq799uP/Oxfft2nTx5Ui1btnSN5ebmav369Zo6daqysrJkt9tdy/z8/OTn51dgOz4+PpX+B/l79MMM+mwGfTaHXptBn1Wi43d7+Ljvvvv0ww8/5Bvr37+//uM//kPDhg3LFzwAAEDl4/bwUa1aNTVv3jzfWJUqVVSzZs0C4wAAoPLhE04BAIBRHnvU9vfWrl1rYjcAAKAc4MwHAAAwivABAACMInwAAACjCB8AAMAowgcAADCK8AEAAIwifAAAAKMIHwAAwCjCBwAAMIrwAQAAjCJ8AAAAowgfAADAKMIHAAAwivABAACMInwAAACjCB8AAMAowgcAADCK8AEAAIwifAAAAKMIHwAAwCjCBwAAMIrwAQAAjCJ8AAAAowgfAADAKMIHAAAwivABAACMInwAAACjCB8AAMAowgcAADCK8AEAAIwifAAAAKMIHwAAwCjCBwAAMIrwAQAAjCJ8AAAAowgfAADAKMIHAAAwivABAACMInwAAACjCB8AAMAowgcAADCK8AEAAIwifAAAAKMIHwAAwCi3h4/x48erdevWqlatmurUqaMePXpo//797t4NAAAop9wePtatW6e4uDht2rRJCQkJcjgc6ty5sy5duuTuXQEAgHLI290bXLFiRb7Xs2fPVp06dbR9+3a1a9fO3bsDAADljNvDxx+lpqZKkmrUqFHo8qysLGVlZblep6WlSZIcDoccDoenyyvz8npALzyLPptBn82h12bQ59+UpAc2y7IsTxXidDrVrVs3nT9/Xhs2bCh0zujRoxUfH19gfP78+QoMDPRUaQAAwI3S09PVu3dvpaamKigo6IpzPRo+nn/+eX355ZfasGGD6tevX+icws58hIeH6/Tp01ctvjJwOBxKSEhQp06d5OPjU9rlVFj02Qz6bA69NoM+/yYtLU21atUqVvjw2GWXQYMGafny5Vq/fn2RwUOS/Pz85OfnV2Dcx8en0v8gf49+mEGfzaDP5tBrM+izSnT8bg8flmXphRde0OLFi7V27Vo1bNjQ3bsAAADlmNvDR1xcnObPn6+lS5eqWrVqSklJkSQFBwcrICDA3bsDAADljNs/52P69OlKTU1Vhw4dFBoa6vqzaNEid+8KAACUQx657AIAAFAUvtsFAAAYRfgAAABGET4AAIBRhA8AAGAU4QMAABhF+AAAAEYRPgAAgFGEDwAAYBThAwAAGEX4AAAARhE+AACAUYQPAABgFOEDAAAYRfgAAABGET4AAIBRhA8AAGAU4QMAABhF+AAAAEYRPgAAgFGEDwAAYBThAwAAGEX4AAAARhE+AACAUYQPAABgFOEDAAAYRfgAAABGET4AAIBRhA8AAGAU4QMAABhF+AAAAEYRPgAAgFGEDwAAYBThAwAAGEX4AAAARhE+AACAUYQPAABgFOEDAAAYRfgAAABGET4AAIBRhA8AAGAU4QMAABhF+AAAAEYRPgAAgFEeCx8ffvihGjRoIH9/f915553asmWLp3YFAADKEY+Ej0WLFmno0KEaNWqUduzYoaioKMXGxurkyZOe2B0AAChHPBI+Jk6cqGeeeUb9+/dXs2bNNGPGDAUGBuof//iHJ3YHAADKEW93bzA7O1vbt2/X8OHDXWNeXl6KiYnRxo0bC8zPyspSVlaW63Vqaqok6ezZs3I4HO4ur9xxOBxKT0/XmTNn5OPjU9rlVFj02Qz6bA69NoM+/+bChQuSJMuyrjrX7eHj9OnTys3NVd26dfON161bVz/99FOB+ePHj1d8fHyB8YYNG7q7NAAA4GEXLlxQcHDwFee4PXyU1PDhwzV06FDXa6fTqbNnz6pmzZqy2WylWFnZkJaWpvDwcCUlJSkoKKi0y6mw6LMZ9Nkcem0Gff6NZVm6cOGCwsLCrjrX7eGjVq1astvtOnHiRL7xEydOKCQkpMB8Pz8/+fn55RurXr26u8sq94KCgir9L7YJ9NkM+mwOvTaDPl92tTMeedx+w6mvr6/uuOMOrV692jXmdDq1evVqRUdHu3t3AACgnPHIZZehQ4eqb9++atWqldq0aaNJkybp0qVL6t+/vyd2BwAAyhGPhI/HHntMp06d0siRI5WSkqIWLVpoxYoVBW5CxdX5+flp1KhRBS5Nwb3osxn02Rx6bQZ9vjY2qzjPxAAAALgJ3+0CAACMInwAAACjCB8AAMAowgcAADCK8FHKLly4oCFDhigyMlIBAQG66667tHXr1iuuk5WVpddff12RkZHy8/NTgwYN+NK+q7iWPs+bN09RUVEKDAxUaGioBgwYoDNnzhiquHxYv369unbtqrCwMNlsNi1ZsiTfcsuyNHLkSIWGhiogIEAxMTFKTEy86nY//PBDNWjQQP7+/rrzzju1ZcsWDx1B+eCJPo8fP16tW7dWtWrVVKdOHfXo0UP79+/34FGUfZ76fc7z9ttvy2azaciQIe4tvBwifJSyp59+WgkJCfr000/1ww8/qHPnzoqJidHx48eLXKdXr15avXq1PvnkE+3fv18LFizQzTffbLDq8qekff7222/15JNP6qmnntKePXv0z3/+U1u2bNEzzzxjuPKy7dKlS4qKitKHH35Y6PJ33nlHU6ZM0YwZM7R582ZVqVJFsbGxyszMLHKbixYt0tChQzVq1Cjt2LFDUVFRio2N1cmTJz11GGWeJ/q8bt06xcXFadOmTUpISJDD4VDnzp116dIlTx1GmeeJPufZunWr/v73v+u2225zd9nlk4VSk56ebtntdmv58uX5xlu2bGm9/vrrha7z5ZdfWsHBwdaZM2dMlFghXEuf3333XatRo0b5xqZMmWLVq1fPY3WWd5KsxYsXu147nU4rJCTEevfdd11j58+ft/z8/KwFCxYUuZ02bdpYcXFxrte5ublWWFiYNX78eI/UXd64q89/dPLkSUuStW7dOneWW265s88XLlywGjdubCUkJFjt27e3Bg8e7KGqyw/OfJSinJwc5ebmyt/fP994QECANmzYUOg6y5YtU6tWrfTOO++oXr16atKkiV566SVlZGSYKLlcupY+R0dHKykpSf/7v/8ry7J04sQJ/c///I/+/Oc/myi5Qjh06JBSUlIUExPjGgsODtadd96pjRs3FrpOdna2tm/fnm8dLy8vxcTEFLlOZXctfS5MamqqJKlGjRpur7EiuJ4+x8XF6YEHHsi3bmVX6t9qW5lVq1ZN0dHRGjt2rJo2baq6detqwYIF2rhxo2666aZC1zl48KA2bNggf39/LV68WKdPn9bAgQN15swZzZo1y/ARlA/X0ue7775b8+bN02OPPabMzEzl5OSoa9euRZ6ORUEpKSmSVOCTjevWreta9kenT59Wbm5uoev89NNPnim0nLuWPv+R0+nUkCFDdPfdd6t58+Zur7EiuNY+L1y4UDt27LjqPWaVDWc+Stmnn34qy7JUr149+fn5acqUKfrLX/4iL6/CfzROp1M2m03z5s1TmzZt9Oc//1kTJ07UnDlzOPtxBSXt8969ezV48GCNHDlS27dv14oVK3T48GE999xzhisHPC8uLk4//vijFi5cWNqlVChJSUkaPHiw5s2bV+DMa2VH+ChlN954o9atW6eLFy8qKSlJW7ZskcPhUKNGjQqdHxoaqnr16uX72uKmTZvKsiwdO3bMVNnlTkn7PH78eN199916+eWXddtttyk2NlbTpk3TP/7xDyUnJxuuvnwKCQmRJJ04cSLf+IkTJ1zL/qhWrVqy2+0lWqeyu5Y+/96gQYO0fPlyff3116pfv75HaqwIrqXP27dv18mTJ9WyZUt5e3vL29tb69at05QpU+Tt7a3c3FyP111WET7KiCpVqig0NFTnzp3TypUr1b1790Ln3X333fr111918eJF19jPP/8sLy8v/uIohuL2OT09vcBZEbvdLuny43a4uoYNGyokJESrV692jaWlpWnz5s2Kjo4udB1fX1/dcccd+dZxOp1avXp1ketUdtfSZ+ny7/GgQYO0ePFirVmzRg0bNjRRbrl1LX2+77779MMPP2jXrl2uP61atVKfPn20a9cu198plVJp3u0Ky1qxYoX15ZdfWgcPHrRWrVplRUVFWXfeeaeVnZ1tWZZlvfrqq9Zf//pX1/wLFy5Y9evXtx555BFrz5491rp166zGjRtbTz/9dGkdQrlQ0j7PmjXL8vb2tqZNm2b98ssv1oYNG6xWrVpZbdq0Ka1DKJMuXLhg7dy509q5c6clyZo4caK1c+dO68iRI5ZlWdbbb79tVa9e3Vq6dKm1e/duq3v37lbDhg2tjIwM1zb+9Kc/WR988IHr9cKFCy0/Pz9r9uzZ1t69e61nn33Wql69upWSkmL8+MoKT/T5+eeft4KDg621a9daycnJrj/p6enGj6+s8ESf/4inXS4jfJSyRYsWWY0aNbJ8fX2tkJAQKy4uzjp//rxred++fa327dvnW2ffvn1WTEyMFRAQYNWvX98aOnRopf4Loziupc9TpkyxmjVrZgUEBFihoaFWnz59rGPHjhmuvGz7+uuvLUkF/vTt29eyrMuPJ44YMcKqW7eu5efnZ913333W/v37820jMjLSGjVqVL6xDz74wIqIiLB8fX2tNm3aWJs2bTJ0RGWTJ/pc2PYkWbNmzTJ3YGWMp36ff4/wcZnNsjiHDAAAzOGeDwAAYBThAwAAGEX4AAAARhE+AACAUYQPAABgFOEDAAAYRfgAAABGET4AAIBRhA8AAGAU4QMAABhF+AAAAEYRPgAAgFH/H8LHAkjtX5JEAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import logging\n",
    "\n",
    "import numpy as np\n",
    "import scipy.signal as ssignal\n",
    "import matplotlib.pyplot as plt\n",
    "#import scipy.interpolate as sinterp\n",
    "\n",
    "import skcomm as skc\n",
    "\n",
    "logging.basicConfig()\n",
    "skc_log = logging.getLogger('skcomm')\n",
    "skc_log.setLevel(logging.INFO)\n",
    "\n",
    "plt.ion()\n",
    "\n",
    "############### PARAMS ########################\n",
    "symbol_rate = 32e9\n",
    "sample_rate_tx = 2*symbol_rate #16e9#12.8e9\n",
    "sample_rate_rx = 2*symbol_rate\n",
    "mod_order = 4\n",
    "roll_off = 0.1\n",
    "\n",
    "linewidth = 200e3 * 1\n",
    "f_offset  = 300e6 * 1\n",
    "cfo_comp  = True\n",
    "\n",
    "snr_seed_X = 123\n",
    "snr_seed_Y = 456\n",
    "SNR = 10\n",
    "\n",
    "# dbg_plt = False\n",
    "dbg_plt = False\n",
    "dbg_dim = 0   # starts with 0,  # -1 means last dim\n",
    "dbg_frame = 0 # starts with 0; -1 means last frame\n",
    "\n",
    "snrs = np.arange(0.0,1.0,1.0) + SNR #np.arange(-5,15,1)#np.arange(2,13,2.0) #np.arange(10,11,1.0)#np.arange(0,13,1.0) #np.arange(5,6,1.0)#np.arange(100,101,1)#\n",
    "bersX = np.ones_like(snrs,dtype=np.float64)\n",
    "bersY = np.ones_like(snrs,dtype=np.float64)\n",
    "\n",
    "for snr_idx, snr in enumerate(snrs):\n",
    "\n",
    "    skc_log.info(f'\\n============\\n SNR: {snr}')\n",
    "    \n",
    "    ############## TRANSMITTER ############################################################################################################################################\n",
    "    # contruct signal\n",
    "    sig_tx = skc.signal.Signal(n_dims=2)\n",
    "    sig_tx.symbol_rate = symbol_rate\n",
    "\n",
    "    # generate bits    \n",
    "    sig_tx.generate_bits(n_bits=int(2**15*np.log2(mod_order)), seed=[1975,1979])\n",
    "    # sig_tx.bits=[np.tile(True,int(2**15*np.log2(mod_order))),np.tile(True,int(2**15*np.log2(mod_order)))]\n",
    "    # set constellation (modulation format)\n",
    "    sig_tx.generate_constellation(format='QAM', order=mod_order)\n",
    "\n",
    "    # create symbols\n",
    "    sig_tx.mapper()\n",
    "    \n",
    "    sig_tx.samples = sig_tx.symbols[:]\n",
    "    sig_tx.sample_rate = sig_tx.symbol_rate[:]\n",
    "    ce_seed = [165,165] #[833,833] # equal seeds for PAPR optimization\n",
    "    pilot_offset = [0,1]\n",
    "    \n",
    "    for dim in np.arange(sig_tx.n_dims):\n",
    "        # insert frame sync header\n",
    "        sig_tx = skc.tx.insert_frame_sync_seq(sig_tx, dim=dim, n_symb=32, fs_header_len=1024, \n",
    "                                              mean_power=1*np.mean(np.abs(sig_tx.constellation[dim])**2), \n",
    "                                              seed=10+dim)\n",
    "\n",
    "        # # insert CE header\n",
    "        n_fft = 128; n_pilots=64; n_repetitions=16\n",
    "        # n_fft =  8; n_pilots= 4; n_repetitions = 8\n",
    "        sig_tx = skc.tx.insert_ch_est_seq(sig_tx, dim=dim, offset=sig_tx.info[dim]['fs_header']['length'], \n",
    "                                          n_fft=n_fft, n_pilots=n_pilots, pilot_offset=pilot_offset[dim], n_repetitions=n_repetitions, \n",
    "                                          mean_power=1*np.mean(np.abs(sig_tx.constellation[dim])**2), seed=ce_seed[dim])\n",
    "        \n",
    "        # insert CPE pilots\n",
    "        sig_tx = skc.tx.insert_cpe_pilots(sig_tx, dim=dim, offset=sig_tx.info[dim]['fs_header']['length']+sig_tx.info[dim]['ce_header']['length'],\n",
    "                                        pilot_distance=100, mean_power=1*np.max(np.abs(sig_tx.constellation[dim])**2), \n",
    "                                        seed=12+dim)    \n",
    "        \n",
    "        skc_log.info(f\"\"\"header lengths:\n",
    "              header length: {sig_tx.info[dim]['fs_header']['length']+sig_tx.info[dim]['ce_header']['length']:d} symbols, \n",
    "              overall frame length: {sig_tx.samples[dim].size:d} symbols, \n",
    "              overhead: {(sig_tx.info[dim]['fs_header']['length']+sig_tx.info[dim]['ce_header']['length'])/sig_tx.samples[dim].size*100:.2f}%\"\"\")\n",
    "\n",
    "        # resample to AWG (Tx) sample rate\n",
    "        upsam = sample_rate_tx/symbol_rate\n",
    "        sr = symbol_rate * upsam\n",
    "        sig_tx.samples[dim] = skc.tx.pulseshaper(sig_tx.samples[dim], upsampling=upsam, pulseshape='rrc', roll_off=roll_off)\n",
    "        sig_tx.sample_rate[dim] = sample_rate_tx\n",
    "\n",
    "        # emulate DAC (quantization and normalization)\n",
    "        n_quantize = None\n",
    "        if n_quantize:\n",
    "            maxVal = np.max(np.abs(np.concatenate((np.real(sig_tx.samples[dim]), np.imag(sig_tx.samples[dim])))))\n",
    "            sig_tx.samples[dim] = sig_tx.samples[dim]/maxVal\n",
    "            # scale i and Q independently (\"Load and Norm\")???\n",
    "            maxVal_r = np.max(np.abs(np.real(sig_tx.samples[dim])))\n",
    "            maxVal_i = np.max(np.abs(np.imag(sig_tx.samples[dim])))\n",
    "            maxVal = np.max([maxVal_r,maxVal_i])\n",
    "            sig_tx.samples[dim] = sig_tx.samples[dim].real/maxVal + 1j*sig_tx.samples[dim].imag/maxVal\n",
    "\n",
    "            sig_tx.samples[dim] = sig_tx.samples[dim] * 2**(n_quantize-1)\n",
    "            sig_tx.samples[dim] = np.round(sig_tx.samples[dim])\n",
    "            sig_tx.samples[dim] = sig_tx.samples[dim] / 2**(n_quantize-1)\n",
    "\n",
    "        # repeat signal (emulate Scope)\n",
    "        repetitions = 4.2 #4.2 # TODO: non-cyclic shift\n",
    "        rem_samples = int((repetitions%1)*sig_tx.samples[dim].size)\n",
    "        sig_tx.samples[dim] = np.concatenate((np.tile(sig_tx.samples[dim], int(repetitions//1)), sig_tx.samples[dim][:rem_samples]), axis=0)\n",
    "\n",
    "    if dbg_plt:\n",
    "        skc.visualizer.plot_signal(np.abs(sig_tx.samples[dbg_dim]), sample_rate=sig_tx.sample_rate[0], fNum=205,tit=f'Signal before channel (abs) (dim = {dbg_dim})')\n",
    "        sig_tx.plot_spectrum(fNum=206,dimension=dbg_dim, tit=f'Spectrum before channel (dim = {dbg_dim})')\n",
    "\n",
    "    ############## CHANNEL ############################################################################################################################################\n",
    "    # b, a = ssignal.bessel(5, 6e9, btype='low', analog=False, output='ba', norm='phase', fs=sample_rate_tx)\n",
    "    # signal = ssignal.lfilter(b,a,signal)\n",
    "    sig_tx.samples[0] = skc.filters.filter_arbitrary(sig_tx.samples[0], FILTER=np.asarray([[0, 6.4e9, 10e9],[0, 0*-1, 0*-5],[0,0* -90,0*-180]]).T,sample_rate=sample_rate_tx)\n",
    "    sig_tx.samples[1] = skc.filters.filter_arbitrary(sig_tx.samples[1], FILTER=np.asarray([[0, 6.4e9, 10e9],[0, 0*-1, 0*-5],[0,0*-180,-0*370]]).T,sample_rate=sample_rate_tx)\n",
    "    \n",
    "    # (pol rotation): linear unitary combination of X and Y pol\n",
    "    rng = np.random.default_rng(seed=None)\n",
    "    alpha = 0 * np.pi/180 #1*0.5*np.pi/2 #0*rng.random()*2*np.pi*1 #4*22.5\n",
    "    phi   = 0 * rng.random()*2*np.pi*1\n",
    "\n",
    "    tmp = sig_tx.samples[0]\n",
    "    sig_tx.samples[0] =  np.cos(alpha) * sig_tx.samples[0]*np.exp(-1j*phi) + np.sin(alpha)*sig_tx.samples[1]\n",
    "    sig_tx.samples[1] = -np.sin(alpha)*tmp + np.cos(alpha)*sig_tx.samples[1]*np.exp(1j*phi)\n",
    "\n",
    "    # add noise\n",
    "    sig_tx.samples[0] = skc.channel.set_snr(sig_tx.samples[0],snr_dB=snr,sps=upsam,seed=snr_seed_X)\n",
    "    sig_tx.samples[1] = skc.channel.set_snr(sig_tx.samples[1],snr_dB=snr,sps=upsam,seed=snr_seed_Y)\n",
    "    ## noise as CW tones on frequency sub-grid of repeated pilot tones => doesn't work :(\n",
    "    # rng = np.random.default_rng(seed=123)\n",
    "    # noise_0 = np.fft.ifft(sig_tx.samples[0].std()/10**(snr/20) * np.exp(1j*2*np.pi*np.random.uniform(low=0,high=1,size=sig_tx.samples[0].shape))) * np.sqrt(sig_tx.samples[0].size * sig_tx.sample_rate[0]/sig_tx.symbol_rate[0])\n",
    "    # sig_tx.samples[0] += noise_0\n",
    "    # rng = np.random.default_rng(seed=345)\n",
    "    # noise_1 = np.fft.ifft(sig_tx.samples[1].std()/10**(snr/20) * np.exp(1j*2*np.pi*np.random.uniform(low=0,high=1,size=sig_tx.samples[1].shape))) * np.sqrt(sig_tx.samples[1].size * sig_tx.sample_rate[1]/sig_tx.symbol_rate[1])\n",
    "    # sig_tx.samples[1] += noise_1\n",
    "        \n",
    "    # add CFO\n",
    "    if f_offset != 0:\n",
    "        sig_tx.samples[0] = skc.channel.add_frequency_offset(sig_tx.samples[0],sample_rate=sr,f_offset=f_offset)\n",
    "        sig_tx.samples[1] = skc.channel.add_frequency_offset(sig_tx.samples[1],sample_rate=sr,f_offset=f_offset)\n",
    "\n",
    "    if linewidth != 0.0:\n",
    "        #seed = np.random.randint(low=1, high=1000)\n",
    "        seed = 789\n",
    "        pn_res = skc.channel.add_phase_noise(sig_tx.samples[0],s_rate=sr, linewidth=linewidth, seed=seed)\n",
    "        sig_tx.samples[0] = pn_res['samples']\n",
    "        # pn_signature = pn_res['phaseAcc']\n",
    "        pn_res = skc.channel.add_phase_noise(sig_tx.samples[1],s_rate=sr, linewidth=linewidth, seed=seed)\n",
    "        sig_tx.samples[1] = pn_res['samples']\n",
    "        # pn_signature = pn_res['phaseAcc']\n",
    "  \n",
    "    if dbg_plt:\n",
    "        sig_tx.plot_spectrum(fNum=211,dimension=dbg_dim, tit=f'Spectrum after channel (dim = {dbg_dim})')\n",
    "\n",
    "    # add delay\n",
    "    channel_delay_0 = 150\n",
    "    sig_tx.samples[0] = np.roll(sig_tx.samples[0], channel_delay_0)\n",
    "    channel_delay_1 = 200\n",
    "    sig_tx.samples[1] = np.roll(sig_tx.samples[1], channel_delay_1)\n",
    "\n",
    "    ############## RECEIVER ############################################################################################################################################\n",
    "    sig_rx = sig_tx.copy()\n",
    "\n",
    "    # # normalize to mean constellation\n",
    "    # sig_rx.samples[0] = sig_rx.samples[0] * np.sqrt(np.mean(np.abs(sig_rx.constellation[0])**2)/np.mean(np.abs(sig_rx.samples[0])**2))\n",
    "    # sig_rx.samples[1] = sig_rx.samples[1] * np.sqrt(np.mean(np.abs(sig_rx.constellation[1])**2)/np.mean(np.abs(sig_rx.samples[1])**2))\n",
    "\n",
    "    # resample to 2 samples per symbol (Rx sample rate)\n",
    "    upsam = sample_rate_rx/sample_rate_tx\n",
    "    sr = sr * upsam\n",
    "    for dim in np.arange(sig_rx.n_dims):\n",
    "        sig_rx.samples[dim] = ssignal.resample(sig_rx.samples[dim],int(sig_rx.samples[dim].size*upsam))\n",
    "        sig_rx.sample_rate[dim] = sample_rate_rx\n",
    "\n",
    "    if dbg_plt:\n",
    "        sig_rx.plot_spectrum(fNum=251, dimension=dbg_dim, tit=f'Spectrum after Rx resampling (dimension {dbg_dim})')\n",
    "        \n",
    "    frame_len = sig_rx.info[0]['fs_header']['length'] + sig_rx.info[0]['ce_header']['length'] + sig_rx.symbols[0].size + sig_rx.info[0]['cpe_pilot_info']['cpe_pilots'].size\n",
    "\n",
    "    results_da_EQ = skc.rx.da_equalizer(sig_tx, sig_rx, cfo_comp=cfo_comp, decimate=True, dbg_plt=dbg_plt, dbg_frame=dbg_frame, dbg_dim=dbg_dim)\n",
    "    sig_rx = results_da_EQ['sig_rx']\n",
    "\n",
    "    # exit()\n",
    "\n",
    "    # pilot_based - CPE\n",
    "    sig_rx, _ = skc.rx.pilot_based_cpe_wrapper(sig_rx, pilot_average=5, tap_len=128, test_phases=91, remove_filter_effects=True, debug=False, fNum=4000)\n",
    "\n",
    "    if dbg_plt:\n",
    "        modMax = np.max(np.concatenate((np.abs(sig_rx.constellation[0].real), np.abs(sig_rx.constellation[0].imag))))\n",
    "        # skc.visualizer.plot_constellation(sig_rx.samples[dbg_dim],fNum=530,hist=True,axMax=modMax*1.5,tit=f'received equalized constellation after CPE (dimension {dbg_dim})')\n",
    "        skc.visualizer.plot_constellation(sig_rx.samples[0],fNum=530,hist=True,axMax=modMax*1.5,tit=f'received equalized constellation after CPE (X Pol)')\n",
    "        skc.visualizer.plot_constellation(sig_rx.samples[1],fNum=531,hist=True,axMax=modMax*1.5,tit=f'received equalized constellation after CPE (Y Pol)')\n",
    "\n",
    "    # delay and phase ambiguity estimation and compensation --> only necessary if signal is cropped after pilot based CPE in line 213 (sig_rx.samples[dim] = sig_rx.samples[dim][filter_len:-filter_len])\n",
    "    # sig_rx = skc.rx.symbol_sequence_sync(sig_rx, dimension=-1)[\"sig\"]\n",
    "\n",
    "    evm = skc.utils.calc_evm(sig_rx)\n",
    "    skc_log.info(f'EVM X: {evm[0]*100:.2f}%, EVM Y: {evm[1]*100:.2f}%')\n",
    "\n",
    "    # decision and demapper\n",
    "    sig_rx.decision()\n",
    "    sig_rx.demapper()\n",
    "\n",
    "    # BER counting\n",
    "    ber_res_X = skc.rx.count_errors(sig_rx.bits[0], sig_rx.samples[0])\n",
    "    ber_res_Y = skc.rx.count_errors(sig_rx.bits[1], sig_rx.samples[1])\n",
    "    skc_log.info(f\"BER X = {ber_res_X['ber']:.2e} ({ber_res_X['err_idx'].sum():d} error(s)), BER Y = {ber_res_Y['ber']:.2e} ({ber_res_Y['err_idx'].sum():d} error(s)), mean BER = {(ber_res_X['ber']+ber_res_Y['ber'])/2:.2e} ({ber_res_X['err_idx'].sum()+ber_res_Y['err_idx'].sum():d} error(s))\")\n",
    "\n",
    "    if dbg_plt:\n",
    "        plt.figure(654448)\n",
    "        if dbg_dim == 0:\n",
    "            plt.plot(ber_res_X['err_idx'][::int(np.log2(mod_order))],label='BER error idx X')\n",
    "        elif dbg_dim==1:\n",
    "            plt.plot(ber_res_Y['err_idx'][::int(np.log2(mod_order))],label='BER error idx Y')\n",
    "        plt.axvline(sig_tx.symbols[dbg_dim].size,ls='--')\n",
    "        plt.axvline(sig_tx.symbols[dbg_dim].size*2,ls='--')\n",
    "        plt.legend()\n",
    "\n",
    "    monitor_num=0\n",
    "    # skc.visualizer.place_figures(monitor_num=monitor_num)\n",
    "\n",
    "    bersX[snr_idx] = ber_res_X['ber']\n",
    "    bersY[snr_idx] = ber_res_Y['ber']\n",
    "\n",
    "\n",
    "plt.figure(600)\n",
    "plt.semilogy(snrs, bersX, 'o', label='BER X')\n",
    "plt.semilogy(snrs, bersY, 'o', label='BER Y')\n",
    "plt.semilogy(snrs, skc.utils.ber_awgn(value=snrs, modulation_format='QAM', modulation_order=mod_order, type='SNR', \n",
    "                                      symbol_rate=symbol_rate, PDM=False), label='theory')\n",
    "plt.legend()\n",
    "plt.grid(visible=True,which='both')\n",
    "plt.ylim([1e-6, 1])\n",
    "\n",
    "plt.figure(601)\n",
    "plt.plot(snrs,skc.utils.ber2q(bersX), 'o', label='Q²_dB (X)')\n",
    "plt.plot(snrs,skc.utils.ber2q(bersY), 'o', label='Q²_dB (Y)')\n",
    "plt.plot(snrs,skc.utils.ber2q(skc.utils.ber_awgn(value=snrs, modulation_format='QAM', modulation_order=mod_order, type='SNR', \n",
    "                                      symbol_rate=symbol_rate, PDM=False)), label='theory')\n",
    "plt.legend()\n",
    "plt.grid(visible=True,which='both')\n",
    "plt.ylim([0, skc.utils.ber2q(1e-6)])"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "noelle311",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}