pam2pyTOB.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "pam2pyTOB.ipynb",
      "provenance": [],
      "collapsed_sections": [],
      "mount_file_id": "10TMNYSFSV-hyzbpfh39Y0GDdU5kG4SvC",
      "authorship_tag": "ABX9TyNAlAl900b5Ee8bM5Xea3ok",
      "include_colab_link": true
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    },
    "language_info": {
      "name": "python"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/DPonirakis/pyNoiseAnalyzer/blob/main/pam2pyTOB.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "2BorlDTLE4dC"
      },
      "source": [
        "# Computes calibrated acoustic spectra from (lossless) digital audio files.\n",
        "# Adapted from PG_Func.m, by Nathan D. Merchant\n",
        "\n",
        "#WAV format             Min       Max        NumPy dtype\n",
        "#32-bit floating-point -1.0      +1.0        float32\n",
        "#32-bit integer PCM  -2147483648 +2147483647 int32\n",
        "#24-bit integer PCM  -2147483648 +2147483392 int32\n",
        "#16-bit integer PCM  -32768      +32767      int16\n",
        "# 8-bit integer PCM   0           255        uint8\n",
        "\n",
        "import numpy as np\n",
        "import scipy as sp\n",
        "import soundfile as sf\n",
        "import warnings\n",
        "warnings.filterwarnings(\"ignore\", message=\"numpy.dtype size changed\")\n",
        "warnings.filterwarnings(\"ignore\", message=\"numpy.ufunc size changed\")"
      ],
      "execution_count": 273,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "aIfx3nnJJXW3"
      },
      "source": [
        "#======================================================================\n",
        "A = []\n",
        "f = []\n",
        "t = []\n",
        "#======================================================================\n",
        "S = 0.0\n",
        "olap = 0.5"
      ],
      "execution_count": 274,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "ouydG_IQJd4I"
      },
      "source": [
        "thesignal, Fs = sf.read('/content/drive/MyDrive/SoundData/PAMguide/WhiteNoise_10s_48kHz_+-0.5.wav',dtype=np.float64) # Let's hope we always get -1 < thedata < 1\n",
        "#thesignal, Fs = sf.read('/content/drive/MyDrive/SoundData/PAMguide/Sine_10s_48kHz_+-0.5.wav',dtype=np.float64) # Let's hope we always get -1 < thedata < 1\n",
        "\n",
        "wlength = 1"
      ],
      "execution_count": 275,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "BiVCzm26gKkL"
      },
      "source": [
        ""
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "4B_MNQYRJ6eZ"
      },
      "source": [
        "samples_channels = np.squeeze( thesignal.shape )\n",
        "if samples_channels.size == 1:\n",
        "    samples  = samples_channels\n",
        "    channels = 1\n",
        "    xbit     = thesignal; thesignal = []\n",
        "else:\n",
        "    samples  = samples_channels[0]\n",
        "    channels = samples_channels[1]\n",
        "    xbit     = thesignal[0,:]; thesignal = []"
      ],
      "execution_count": 276,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "JSdJfng2KH93"
      },
      "source": [
        "xbit = np.float32( xbit ) # Single precision\n",
        "nyquist = int( Fs/2 )\n",
        "N = int( Fs*wlength )\n",
        "\n",
        "Mh = -36-120\n",
        "#Mh = -200\n",
        "Vadc = 1.414\n",
        "#Vadc = 2\n",
        "#G = 0\n",
        "G = 0\n",
        "\n",
        "lcut = 1\n",
        "hcut = nyquist\n",
        "\n",
        "winname = 'Hann'\n",
        "\n",
        "pref = 20\n",
        "#pref = 1"
      ],
      "execution_count": 277,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "Nb5HtMDN-cPY"
      },
      "source": [
        "S = Mh + G + 20*np.log10(1/Vadc)\n"
      ],
      "execution_count": 278,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "PhS8a83aLC6X"
      },
      "source": [
        "xl = len(xbit)\n",
        "\n"
      ],
      "execution_count": 279,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "8DRqs2Yyg5Tn"
      },
      "source": [
        "# Mimic MATLAB routine to generate buffer array"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "yZf0xl1bMO1R"
      },
      "source": [
        "def buffer(X, n, p=0, opt=None):\n",
        "    '''Mimic MATLAB routine to generate buffer array\n",
        "\n",
        "    MATLAB docs here: https://se.mathworks.com/help/signal/ref/buffer.html\n",
        "\n",
        "    Parameters\n",
        "    ----------\n",
        "    x: ndarray\n",
        "        Signal array\n",
        "    n: int\n",
        "        Number of data segments\n",
        "    p: int\n",
        "        Number of values to overlap\n",
        "    opt: str\n",
        "        Initial condition options. default sets the first `p` values to zero,\n",
        "        while 'nodelay' begins filling the buffer immediately.\n",
        "\n",
        "    Returns\n",
        "    -------\n",
        "    result : (n,n) ndarray\n",
        "        Buffer array created from X\n",
        "    '''\n",
        "\n",
        "    n = int(n)\n",
        "\n",
        "    import numpy as np\n",
        "\n",
        "    if opt not in [None, 'nodelay']:\n",
        "        raise ValueError('{} not implemented'.format(opt))\n",
        "\n",
        "    i = 0\n",
        "    first_iter = True\n",
        "    while i < len(X):\n",
        "        if first_iter:\n",
        "            if opt == 'nodelay':\n",
        "                # No zeros at array start\n",
        "                result = X[:n]\n",
        "                i = n\n",
        "            else:\n",
        "                # Start with `p` zeros\n",
        "                result = np.hstack([np.zeros(p), X[:n-p]])\n",
        "                i = n-p\n",
        "            # Make 2D array and pivot\n",
        "            result = np.expand_dims(result, axis=0).T\n",
        "            first_iter = False\n",
        "            continue\n",
        "\n",
        "        # Create next column, add `p` results from last col if given\n",
        "        col = X[i:i+int(np.floor(n-p))]\n",
        "        if p != 0:\n",
        "            col = np.hstack([result[:,-1][-p:], col])\n",
        "        i += n-p\n",
        "\n",
        "        # Append zeros if last row and not length `n`\n",
        "        if len(col) < n:\n",
        "            col = np.hstack([col, np.zeros(n-len(col))])\n",
        "\n",
        "        # Combine result with next row\n",
        "        result = np.hstack([result, np.expand_dims(col, axis=0).T])\n",
        "\n",
        "    return result"
      ],
      "execution_count": 280,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "kJC8HyFjMQ0m",
        "outputId": "fe016c77-2b36-4db1-9344-92434b58e771"
      },
      "source": [
        "#xbit[N:]\n",
        "np.ceil(N*olap)\n",
        "xgrid = buffer(xbit,N,int(np.ceil(N*olap)),'nodelay')\n",
        "xgrid = np.transpose(xgrid)\n",
        "\n",
        "np.shape(xgrid)\n"
      ],
      "execution_count": 281,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "(19, 48000)"
            ]
          },
          "metadata": {},
          "execution_count": 281
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "y78REJVrlCmT",
        "outputId": "e2a1f49b-c6dc-4d72-d31b-547fb905df44"
      },
      "source": [
        "M = len(xgrid[:,1])\n",
        "M"
      ],
      "execution_count": 282,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "19"
            ]
          },
          "metadata": {},
          "execution_count": 282
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "xexN8XDIX6lw"
      },
      "source": [
        "# Equation (6) - Create a Hann Window"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "3Mxb-Zh0Ufy7"
      },
      "source": [
        "#create a Hann window\n",
        "\n",
        "w = (0.5 - 0.5*np.cos(2*np.pi*np.linspace(1,N,N)/N))\n",
        "alpha = 0.5 #scaling factor"
      ],
      "execution_count": 283,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 265
        },
        "id": "TrwDO977V6jx",
        "outputId": "473738ee-9bc8-468b-fdc2-b959c5e545f2"
      },
      "source": [
        "import matplotlib.pyplot as plt\n",
        "\n",
        "plt.plot(w, linestyle = 'dotted')\n",
        "plt.show()"
      ],
      "execution_count": 284,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYEAAAD4CAYAAAAKA1qZAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3deXyU1b3H8c8v+54QErYs7CoILhhRQAUULS6Vam2rXWxtLVWr9dYu1157e2vb21ar2NparbVerVWpa0utS3EBV5Ygyir7FkQIWwiErHPuHxnSgEBCmOTMzPN9v168mOeZh+T7hJn85pzzPOeYcw4REQmmBN8BRETEHxUBEZEAUxEQEQkwFQERkQBTERARCbAkX9+4oKDA9evXz9e3FxGJSfPmzdvqnCuM1NfzVgT69etHeXm5r28vIhKTzGxdJL+euoNERAJMRUBEJMBUBEREAkxFQEQkwFQEREQCrM0iYGYPmtkWM1t0iOfNzO42s5VmtsDMRkQ+poiIdIb2tAQeAiYe5vnzgcHhP5OBe48+loiIdIU2i4Bz7nVg+2EOmQT82TWbBeSZWe9IBRTpbBu21/Dmiq0t269+sJlH3lnbsr1jTz2NTaGuDybSBSIxJlAEbGi1XRHe9zFmNtnMys2svLKyMgLfWuTILd9czQNvrG7Zvm/mKr41dX7L9nMLNnHfzH8//6Npizn3rtdbtqv2NnRNUJEu0KV3DDvn7gfuBygrK9NqNtJldtbUk5OWTEKCMXNZJb94YSmTTiqiMDuVq8b05/JTS1uOveOyE6lv9cn/0yOKOHNQQcv2Fx6YRXFeBvd96ZQuPQeRzhCJIrARKGm1XRzeJxIV3t+wk8vvn8XvvzCC8cf14HMjS7jslGK6ZaYAMKhH1n7HJyQYaQmJLdvjju3R8tg5x+fKSuiZkwZAY1OIe2es4vOnldI9K7ULzkYksiLRHTQNuDJ8ldDpQJVzblMEvq5Ih+2ua2T55moAhvbJ4XOnltC3ewYAOWnJLQXgSJkZXxrVj/OO7wXAvHU7uOvl5ZSv2xGZ4CJdzNpaY9jMHgfGAQXAZuB/gGQA59x9ZmbA72i+gqgGuMo51+bMcGVlZU4TyEln+ewf3mHr7jqmf3ssiQnWqd9r7dY99O2egZnxwsJNZKUlcebgiE3yKLIfM5vnnCuL1NdrszvIOXdFG8874JuRCiTSUVU1DWSlJZGYYNx8/nEAnV4AAPoVZALNXUX3v7Ga5IQEzhhUQPPnI5Ho5m0qaZFI2lJdyyX3vM2Vo/ryjbEDGVHarcszmBmPf/10qvY2YGbU1DeyblsNQ3rndHkWkfbStBESFwqzUrlgeC9GDezuNUdacmLLoPGUfy3nU/e8xUdVtV4ziRyOWgISsxqbQtz18nK+PLofPbLTuOXCob4j7ee68YM4sSSPXrlpvqOIHJJaAhKz1m7bw/+9tZZ/Ld7sO8pB5Wem8MkT+wCwaGMV1z/2LnvrmzynEtmfWgIScxqbQiQlJjCoRzavfGcsvXPTfUdq09JNu3i/YidVextIT0ls+x+IdBG1BCSmVFbXcdFv3+S5BR8CxEQBAPhMWQnTvz2WXrlpOOfYXdfoO5IIoCIgMSYjJZFeuWkUxODduWnJzS2A3726kk///m3NQSRRQd1BEhN27KknOy2JzNQkHrpqpO84R+WUvt3Ytqee7FS9/cQ/vQol6tU1NvGFB2YzqEcWd19xsu84R230oAJGhyek21lTT3JiApkqCOKJXnkS9VKTErliZAml3TN9R4moxqYQX/zTbAqzUnnwK6fqDmPxQkVAolZjU4iPdtVS3C2DL43q5ztOxCUlJvC1M/rTIztNBUC80cCwRK0p05dz8e/eYtvuOt9ROs0lJxczJtw1tHHnXs9pJIhUBCRqXXZKMZPPGhCIefrL125n/K9m8M8FmoVdupaKgESdreFP/gMKs7hm7EDPabrGCcV5TD5rAGMG+Z37SIJHRUCiSsWOGs6dMnO/NYCDICUpge9+4ljyMlIIhRzVtbqHQLqGioBElR7ZaVw+spQJQ3r6juLNDVPn87WHy2lstc6xSGfR1UESNRqaQqQkJfCfE4/zHcWrc4f0ZFdtQ5csiCOiloBEhWfnV/DJ377JlmrNvf+pk4u4clQ/zIxQ6PDLv4ocLRUBiQr5man07Z5BfkbHFoCPR/PW7eATv36dih01vqNIHFN3kESFsccUMvYYLc7eWkFWChmpSdQ2aA0C6TwqAuLVL1/4gH7dM7h8ZKnvKFGnb/dM/nbdaN1NLJ1K3UHiTUNTiIUbd7Jsc7XvKFHLzGhoCnH7ix/w9qqtvuNIHFJLQLxJTkzgka+eRqMGPw+roSnEC4s+AmD0wALPaSTeqAhIl3PO8dDba7l0RDG56cmk6FLIw8pISeLv148hJy3ZdxSJQ+oOki63+MNd/OyfS/n7ext9R4kZ+wrAhu01lK/d7jmNxBO1BKTLDSvK5YUbz2RQYZbvKDHn2399j+176pl+01jdTCYRoSIgXcY5x/rtNfTtnskxPbN9x4lJv/z0cFKTElUAJGLUHSRd5vmFH3H2nTOZs0bdGR01qEc2JfkZQPPSlCJHS0VAuswZgwr41tmDGVGa5ztKzJsyfTnn/+YNdtc1+o4iMU7dQdIlnHPkZiRz44TBvqPEhXHHFoJzJCeqW0iOjloC0ule/WAzVz44J66XiexqI0q7cdN5x5KalOg7isS4dhUBM5toZsvMbKWZ3XyQ50vN7DUzm29mC8zsgshHlVhVtbeBPXWNZOs694hbWFHFDY/Pp75Raw9Ix7RZBMwsEbgHOB8YClxhZkMPOOyHwBPOuZOBy4HfRzqoxK5LTi7m6WtHk5Kkhmekbd5Vy5w121i/fY/vKBKj2jMmMBJY6ZxbDWBmU4FJwJJWxzggJ/w4F/gwkiElNq2q3M36bTWMP66HJkHrJBOG9mTMoALSU9QtJB3Tno9mRcCGVtsV4X2t/Rj4oplVAM8DNxzsC5nZZDMrN7PyysrKDsSVWPLH11dz49T5VO3VermdKT0lkVDI8cYKvafkyEWqfX4F8JBzrhi4AHjEzD72tZ1z9zvnypxzZYWFmjs+3v1k0jAevfp0ctM1FtDZnn63gi/9aY7uwZAj1p7uoI1ASavt4vC+1r4GTARwzr1jZmlAAbAlEiElttTUN5KcmEBKUgLDi3N9xwmESScVkZacSFnfbr6jSIxpT0tgLjDYzPqbWQrNA7/TDjhmPXAOgJkNAdIAtU0D6o6XlvPJ376pFbG6UEpSAp88sQ8JCVqXWI5Mm0XAOdcIXA+8BCyl+SqgxWb2EzO7OHzYd4Cvm9n7wOPAV5xzeiUG1OiB3blgeG/SkjVY2dXmrNnOOVNmsnHnXt9RJEa0645h59zzNA/4tt73o1aPlwBjIhtNYtWEoT2ZMLSn7xiBVNQtncKsVPZoOglpJ00bIRHz5oqtrN66m8+PLCUpUfcE+FCUl84T14zyHUNiiN6pEjH/XPgh//fWWtQl7d/e+iYeeWetxgekTWoJSMT8/JLhbN1drzuDo8ArH2zmv/++mIGFWYwepHWJ5dBUBOSo7bsZLDc9mcLsVM9pBODC4b0pvT6DE4o1bbccnj6yyVG781/LmDBlpgYjo4iZtRQAXaorh6MiIEftc6eWcMPZg8hMVcMy2ry4aBNjfvkqm6p0yagcnN61ctSO75PL8X10Z3A0Or5PLmM0JiCHoZaAdNjctdv5+fNLtcRhFCvJz+DuK06md2667ygSpVQEpMPK1+7gufc/JClB00RHu01Ve/nDzFXoRn45kIqAdNi14wYy/aaxmh4iBry06CPu/NdyVm/V4jOyP40JyBFraArxUVUtJfkZGgyOEZ8/rS/nDOlJSX6G7ygSZdQSkCM2de4Gzr5zBis2V/uOIu2UkpTQUgD21uuSUfk3fYyTI3bukJ7sqWtkUI8s31HkCP3m5RU8M7+C6d8eqzu7BVARkA7olZvGNWMH+o4hHTCibx41DY00hkKkqCNAUHeQHIHK6jq+/df3qNhR4zuKdNCZgwv5wflDyEjR5z9ppiIg7bZw405eWbqZ+saQ7yhylBZU7ORv8w9cJVaCSB8HpN3OPq4ns/9rAukpuiQ01t03cxULKqq46ITeWvsh4FQEpF02bK+hJD9DBSBO/PiTx5OekqgCIOoOkrat2bqHcXfM4NHZ63xHkQjpkZNGdloyzjnqGnXJaJCpCEibemSn8p3zjuG8ob18R5EICoUcX/zTbH48bYnvKOKRuoOkTZmpSVw3bpDvGBJhCQnGyH7dtRBQwKkIyGFN+dcyRg0sYNTA7r6jSCe4ccJg3xHEM3UHySHtqm3gqXkVzFmz3XcU6UTOOV5ctInVlbt9RxEPVATkkHLSknn1u+OYfNYA31GkE+2oaeCmJ97nkVka+A8idQfJQW3bXUe3jBRNEx0A+ZkpPPGNURzXK9t3FPFALQH5GOcc1z76Llc9NNd3FOkiw4pySUpMIBTSojNBoyIgB/WF00q5dESR7xjShZZ8uIuz75zBwooq31GkC6k7SD7GzJh0kgpA0JTkp9MnL536Js0NFSQqArKfOWu2s3bbHi49uUhTCgRMdloyj339dN8xpIvpXS77eXZ+Bb+evpwmLUgeWLUNTbywcJPvGNJF2lUEzGyimS0zs5VmdvMhjvmsmS0xs8Vm9lhkY0pX+fklw3nq2tGkJumqoKB6fM56rn30XZZ8uMt3FOkCbXYHmVkicA9wLlABzDWzac65Ja2OGQz8ABjjnNthZj06K7B0Ducc9U0hUpMS6ZOX7juOePS5U0s4rlcOQ/vk+I4iXaA9LYGRwErn3GrnXD0wFZh0wDFfB+5xzu0AcM5tiWxM6WzvrN7GGbe9xuIPdWVI0GWkJGmakABpTxEoAja02q4I72vtGOAYM3vLzGaZ2cSDfSEzm2xm5WZWXllZ2bHE0imyU5MZ2S+fgYVaPF6aPfTWGr71+HzfMaSTRWpgOAkYDIwDrgD+aGZ5Bx7knLvfOVfmnCsrLCyM0LeWSBhenMs9XxihO4SlRW1jiJr6RmobtN5APGtPEdgIlLTaLg7va60CmOaca3DOrQGW01wUJAY8Na+Cqr0NvmNIlPnGWQN44Mun6oNBnGtPEZgLDDaz/maWAlwOTDvgmL/R3ArAzApo7h5aHcGc0klWbK7mu0++z1PzKnxHkShjZgBs3V3H+m01ntNIZ2nz6iDnXKOZXQ+8BCQCDzrnFpvZT4By59y08HPnmdkSoAn4nnNuW2cGl8gY3DOb5244Q2MBclBNIcek373Fsb2yefArp/qOI53AnKebgsrKylx5ebmX7y3NnHMtn/ZEDuXlJZvpX5ipDwpRwszmOefKIvX1dMdwgF37l3f53asrfMeQKDdhaE8VgDimIhBQDU0hMlITdWewtMv2PfX84JmFLKjY6TuKRJgmkAuo5MQEpnz2JN8xJEYkJxqvLN3M8KJcTij+2NXfEsNUBAJozdY9JCUYJfkZvqNIjMhOS+b174/X5aJxSN1BAfSL55dyye/fpkHzxssR2FcANlXt9ZxEIklFIIBunXQ8d372RJK1XoAcoRcWbmLML19l0UbNMRUv1B0UQL1z0+mdq5lC5ciNGVzAN8cP0kyzcUQfBQNk8YdVXP/Yu3xUVes7isSonLRkvnPeseRnpviOIhGiIhAgK7fsZt66HaSnaHBPjs77G3Zy74xVvmNIBKgIBMikk4p4/fvjyU1P9h1FYtwrSzfzxzdWU1WjiQdjnYpAQKzdugdAg8ESEZPHDmz+QJGhDxSxTr8RAmDZR9WMv3MGT5RvaPtgkXbISk0iKzUJ5xx767XeQCxTEQiA4m7p3DzxOM4d0tN3FIkjzjmuemgu3396ge8ochR0iWgAZKYm8Y2xA33HkDhjZpwxqIBU3UUc01QE4tw9r63kpJI8xgwq8B1F4tDVZw7wHUGOkrqD4tje+iYen7Oe11dU+o4icSwUcry4aBNrwhcfSGxREYhj6SmJvPKdsdxwtpZ7ls6zo6aeb//1fR6dtc53FOkAdQfFqd11jWQkN68XkKr/ZelE3bNSefKaUQzpneM7inSAWgJx6sfTFnPJvW8TCvlZPlSCZVhRLokJhq/laqXjVATi1FnHFHLR8N4kJGgNYeka89fv4JwpM1m3TWMDsUQdBXHq4hP7+I4gAVOUl05+Rgq79jb6jiJHQC2BOLNhew1PzN2gBWOky/XISeOpa0czvDjXdxQ5AioCceaZdzfyw78vYvueet9RJKBqG5p4c8VW3zGknVQE4sy3zhnEP284g545ab6jSEDd9fJyvvJ/c9i8S+tWxAKNCcSRppAjMcEY3DPbdxQJsKtG92fcMT30QSRGqCUQJzbu3MtZt7/G2yvVDBe/euWmMWpgd98xpJ1UBOLE3vpGBvbIom9Bpu8oIgDcO2MVt734ge8Y0gZ1B8WJQT2y+fNXR/qOIdJiw44aqmoacM5hpvtVopWKQBx47YMtjOjbTctGSlT56aRhJOpmxain7qAYt7Omnm/8ZR5T/rXMdxSR/ewrAFt317Ftd53nNHIoKgIxLi8jhWeuHa1FYyQq7a5rZPyvZvDbV1f6jiKH0K4iYGYTzWyZma00s5sPc9ynzcyZWVnkIkpbhhXl0icv3XcMkY/JSk3ivy8aypWj+vqOIofQ5piAmSUC9wDnAhXAXDOb5pxbcsBx2cCNwOzOCCofd/uLHxBy8J8Tj9XAm0Stz55a4juCHEZ7WgIjgZXOudXOuXpgKjDpIMf9FLgN0G2CXWTn3gZ21tSrAEjU+6iqllueXciWav16iDbtKQJFwIZW2xXhfS3MbARQ4pz75+G+kJlNNrNyMyuvrNSSh0fr55cM5xeXDvcdQ6RNtQ1NPPPuRt5dt8N3FDnAUV8iamYJwBTgK20d65y7H7gfoKysTKtPdND2PfXsqWukJD9DrQCJCf0KMpl9yznkpOky5mjTnpbARqB1p15xeN8+2cAwYIaZrQVOB6ZpcLjz3DtjJROmzNRMoRJT9hWAnTV63UaT9hSBucBgM+tvZinA5cC0fU8656qccwXOuX7OuX7ALOBi51x5pyQWrj5zAL+4dDj5mSm+o4gckT+/s5YzbntN9w1EkTaLgHOuEbgeeAlYCjzhnFtsZj8xs4s7O6B8XM+cNC4dUew7hsgRGz2wgCtH9SUpUbcoRQvztTB0WVmZKy9XY+FIbNy5l1++8AH/OfFYirtl+I4jIh6Y2TznXMS621WOY8iSD3fx5opKEjQYLDFu/vodPDu/wncMQRPIxZRzh/bknR+cQ1pyou8oIkflgTfXsKBiJ588oY+6hjxTEYgR67fVUNo9QwVA4sL/XDSUjNQkFYAooP+BGLC6cjfj75zBY7PX+44iEhE9ctLISk3COUdtQ5PvOIGmIhADeuWm8b1PHMt5x/f0HUUkYkIhxxcemM1Pn1vS9sHSadQdFAMyUpK4RlNFS5xJSDBG9s+nlxak90pFIMpNmb6c0wfkM3pgge8oIhH3HxOO8R0h8NQdFMWqaxt4el4Fs1Zt8x1FpNOEQo7nF25i5ZbdvqMEkopAFMtOS+bV747lmnHqCpL4tau2ge8/tYC/zFrnO0ogqTsoSlVW19E9M4XUJF0SKvEtLyOFJ68ZxTE9s31HCSS1BKKQc47Jj5TztYfn+o4i0iWG9M4hMcFobAr5jhI4KgJR6qox/flMmZblk+BY9lE14++cwdy1231HCRR1B0UhM+PiE/v4jiHSpUrzMxhYmEWCpsbqUioCUeaFhZvYVdvAZ04pIUHvBgmQ9JREHrpqpO8YgaPuoCjzjwUf8tjs9WiiUAmq2oYmHp29TuMDXUQtgShzz+dHsKOmQWsHS2C9sWIrtzy7iD656Yw/rofvOHFPRSBK7KlrxKx5iggtGylBNmFID569bjQnl3bzHSUQ1B0UJe6dsYpxv5pB1d4G31FEvDKzlgLQoC6hTqciECXOGdKDL53el9z0ZN9RRKLCi4s+4qzbX2OrFqXvVOoOihInl3ZT81eklcE9sxhelEtdo1oDnUktAc8WVlRxx0vLqKlv9B1FJKoMLMzi/ivLKMpL9x0lrqkIePbGykqmzl1PU8j5jiISlSqr65g6R6vqdRYVAc+uGzeIl28aS3aaxgJEDmbqnPX88G+L2LC9xneUuKQxAU/qG0Ns3V1Hn7x08jJ0SajIoVx95gDOH96bkvwM31HikloCnjwyax3j75jB2q17fEcRiWrpKYkM6pEFNN9PI5GlIuDJJ47vyY0TBtO3uz7diLTHg2+uYfwdM9hVq3tpIkndQZ4Ud8vgunGDfMcQiRmn9stn3Ta1nCNNLYEutvjDKm6cOp9tugFG5IgML87l1knDyNFFFBGlItDFlny4i7lrtpOUqB+9SEcs31zNr19e7jtG3NBvoi72mbISXvveOE0PIdJBr32whQffXMOmqr2+o8SFdhUBM5toZsvMbKWZ3XyQ528ysyVmtsDMXjGzvpGPGtuqaxuYv34HgBaPFzkKXxnTjxnfG0/vXN1JHAltFgEzSwTuAc4HhgJXmNnQAw6bD5Q5504AngJuj3TQWPfHN9bw6Xvf1g0vIkcpNSmR/MwUnHOsrtztO07Ma09LYCSw0jm32jlXD0wFJrU+wDn3mnNu32+3WUBxZGPGvm+cNYB7Pj9CN7yIRMgf31jNxN+8wfpt+mB1NNpziWgRsKHVdgVw2mGO/xrwwsGeMLPJwGSA0tLSdkaMbaHwnECZqUmcP7y35zQi8WPSSUUkmNEnL813lJgW0YFhM/siUAb86mDPO+fud86VOefKCgsLI/mto9YT5Ru47L632VlT7zuKSFzpmZPG1WcO0JV2R6k9P72NQEmr7eLwvv2Y2QTgFuBi55wugg/LTkumV26argYS6STz1m3nc394h2rdSdwh7ekOmgsMNrP+NP/yvxz4fOsDzOxk4A/AROfcloinjGEXntCbC09QN5BIZ0lKSKCyuo5NVbWajbcD2iwCzrlGM7seeAlIBB50zi02s58A5c65aTR3/2QBT5oZwHrn3MWdmDvqvbliK9v21HHxiX0I/0xEpBOcWJLH9JvGkpig91lHtGvuIOfc88DzB+z7UavHEyKcK+Y9Nmcdyzfv5vxhvUlJ0otTpDMlJhiNTSGmzt3AZacUk5ase3HaSxPIdZLfXjGCLdW1pCRp0EqkK7y3YSc//Nsi0pITuewUXaXeXioCEbZiczVF3dLJSEnSHY0iXaisXz5//+YYTijO9R0lpuhjagQ1NIX46sNzufYv7/qOIhJIJ5bkYWZUVtdR19jkO05MUEsggpITE7jjshNJVX+kiDeV1XWce9dMrjy9Lzedd6zvOFFPRSBC6htDpCQlcNqA7r6jiARaYXYqXz9zAOcP6+U7SkxQd1AE7NhTz7l3zeTpeRW+o4gI8M3xgxhQ2Lwu8b6pW+TgVAQiZGjvHI7tle07hoi08tPnlvDdJ9/3HSOqqTsoArplpnDvF0/xHUNEDpCbnoxz0BRyupnsEFQEjsLSTbu4+5UV/O8lw8nPTPEdR0QOcMPZg3THfhvUHXQUlm+uZuHGKkJOfY4i0WhfAVhduZubnniP2gZdNnogtQSOwqSTipg4rJeWixSJcqsq9zBjWSVrtu5hSO8c33GiiopABzwxdwPF+emMHligAiASA84d2pPTB4zTLKMHoe6gI9TQFOLBt9bw0FtrfUcRkSOwrwA8Ub6BDz7a5TlN9FBL4AglJybw1LWjNQ4gEoN21TZwx0vLOGdIT35x6XDfcaKCikA77a1v4tHZ67hqTH+yUvVjE4lFOWnJPHnNKIryNLnjPuoOaqfnF27if59fyvz1O3xHEZGj0Ld7JkmJCVTXNvDgm2twAW/V6yNtO336lGKG9snRlQUiceKZdzfy8+eXctqAfI7vE9zpp1UE2vD39zZyckk3SrtnqACIxJErR/Xl1H75DO0T7Pe1uoMOY1dtA7f+Ywm/eWWF7ygiEmFm1lIAZq3exvQlmz0n8kMtgcPISUvmqWtG0UeDSCJxyznHr19eTnVtI2cf1yNwcwypCBzEzOWVbNyxl8+fVtoyHa2IxCcz474vnoJzBK4AgLqDDuqJuRv4y6x11DeGfEcRkS6Ql5FCt8wUnHPc+o/FvLGi0nekLqOWQCuhkCMhwbjzsydS19C8UpiIBMfuukbeWbWNrNQkzhxc6DtOl1ARCHtk1jpeWLiJB79yKmnJiaRpnWCRwMlOS+bpa0eTkdL8/t9V20BOnM83pI+6YTlpSboTWETITE3CzNi+p54L736D389Y6TtSpwr0b70de+pZs20PI0q7MemkIi4+sY8WoBARoPmD4TnH9WT0wALfUTpVoFsCP3hmIVc/XM6eukYAFQARaZGUmMCPLz6ek0ryAHhs9noWVlR5ThV5gWsJ1DeGaAyFyEhJ4pYLh7CzpoFMdQOJyGHUNjTx+xkrGdkvnymfO8l3nIgK1G+/+sYQn7rnLU4uzeN/LxlOSX4GJfm+U4lItEtLTuQf159BUmJzb0HFjhp27W2MiyknAtEdVLW3AYCUpAQuPKE344/t4TmRiMSabpkpLQvT3P7iMq744yz21sf+msVxXwRmLNvC6T9/hRWbqwH45vhBTBja03MqEYllt158PHdfcTLp4UtJ//zOWj7cuddvqA5qVxEws4lmtszMVprZzQd5PtXM/hp+fraZ9Yt00PZqaArxwsJNvLdhJwDDinK5+MQ+WltURCKmW2YKY49pvpls/bYabv3HEp5fuAlo/h1U1xg7LYQ2i4CZJQL3AOcDQ4ErzGzoAYd9DdjhnBsE3AXcFumgh/Pehp3MWbO9Zfu/nl3I1DnrASjISuW2y06gV25aV0YSkYAo7Z7BzO+N4/KRpQC8+sEWyn72Miu37AZgZ00923bXRe3iNe1pCYwEVjrnVjvn6oGpwKQDjpkEPBx+/BRwjnXi9Zbfe/J9rn64vGX7Z88t4ZcvLAX+vQbwzz41rLO+vYjIfoq7ZbTcbFqUl85FJ/SmND8DgEdnr+eUn71MTXj8YFXlbm85D6Y9VwcVARtabVcApx3qGOdco5lVAd2Bra0PMrPJwGSA0tLSDkaGY3tlU13b2KT4QYkAAAWFSURBVLL9s0uGkZ+R0rI9UDN/iognw4py+cWlJ7Rsjz+2B7npyS2XoqdG2ZxkXXqJqHPufuB+gLKysg63ja4+c8B+28f1iv3LtEQkPg3tk7PfpaTF3TI8pvm49pSkjUBJq+3i8L6DHmNmSUAusC0SAUVEpPO0pwjMBQabWX8zSwEuB6YdcMw04Mvhx5cBr7poHQUREZEWbXYHhfv4rwdeAhKBB51zi83sJ0C5c24a8CfgETNbCWynuVCIiEiUa9eYgHPueeD5A/b9qNXjWuAzkY0mIiKdLbqGqUVEpEupCIiIBJiKgIhIgKkIiIgEmPm6ktPMKoF1HfznBRxwN3LABPn8g3zuEOzz17k36+ucK4zUF/ZWBI6GmZU758p85/AlyOcf5HOHYJ+/zr1zzl3dQSIiAaYiICISYLFaBO73HcCzIJ9/kM8dgn3+OvdOEJNjAiIiEhmx2hIQEZEIUBEQEQmwmCsCbS16HyvM7EEz22Jmi1rtyzez6Wa2Ivx3t/B+M7O7w+e8wMxGtPo3Xw4fv8LMvtxq/ylmtjD8b+7uzOU+j5SZlZjZa2a2xMwWm9mN4f1BOf80M5tjZu+Hz//W8P7+ZjY7nPmv4anbMbPU8PbK8PP9Wn2tH4T3LzOzT7TaH9XvEzNLNLP5ZvZceDsQ525ma8Ovy/fMrDy8z+/r3jkXM39onsp6FTAASAHeB4b6ztXBczkLGAEsarXvduDm8OObgdvCjy8AXgAMOB2YHd6fD6wO/90t/Lhb+Lk54WMt/G/P933Orc6zNzAi/DgbWA4MDdD5G5AVfpwMzA5nfQK4PLz/PuDa8OPrgPvCjy8H/hp+PDT8HkgF+offG4mx8D4BbgIeA54Lbwfi3IG1QMEB+7y+7mOtJdCeRe9jgnPudZrXXmhtEvBw+PHDwKda7f+zazYLyDOz3sAngOnOue3OuR3AdGBi+Lkc59ws1/zK+HOrr+Wdc26Tc+7d8ONqYCnN61QH5fydc27fauPJ4T8OOBt4Krz/wPPf93N5Cjgn/AlvEjDVOVfnnFsDrKT5PRLV7xMzKwYuBB4IbxsBOfdD8Pq6j7UicLBF74s8ZekMPZ1zm8KPPwJ6hh8f6rwPt7/iIPujTrh5fzLNn4YDc/7h7pD3gC00v4lXATudc43hQ1pnbjnP8PNVQHeO/OcSLX4NfB8Ihbe7E5xzd8C/zGyemU0O7/P6uu/Sheal/Zxzzszi+vpdM8sCngb+wzm3q3X3Zbyfv3OuCTjJzPKAZ4HjPEfqEmZ2EbDFOTfPzMb5zuPBGc65jWbWA5huZh+0ftLH6z7WWgLtWfQ+lm0ON+kI/70lvP9Q5324/cUH2R81zCyZ5gLwqHPumfDuwJz/Ps65ncBrwCiam/v7Ppi1ztxynuHnc4FtHPnPJRqMAS42s7U0d9WcDfyGYJw7zrmN4b+30Fz8R+L7de97oOQIB1WSaB4E6c+/B32O953rKM6nH/sPDP+K/QeIbg8/vpD9B4jmuH8PEK2heXCoW/hxvjv4ANEFvs+31Xkazf2Vvz5gf1DOvxDICz9OB94ALgKeZP/B0evCj7/J/oOjT4QfH8/+g6OraR4YjYn3CTCOfw8Mx/25A5lAdqvHbwMTfb/uvf9gOvCDvIDmq0lWAbf4znMU5/E4sAlooLnv7ms093W+AqwAXm71H2vAPeFzXgiUtfo6X6V5UGwlcFWr/WXAovC/+R3hu8Oj4Q9wBs19owuA98J/LgjQ+Z8AzA+f/yLgR+H9A8Jv4pXhX4qp4f1p4e2V4ecHtPpat4TPcRmtrgSJhfcJ+xeBuD/38Dm+H/6zeF823697TRshIhJgsTYmICIiEaQiICISYCoCIiIBpiIgIhJgKgIiIgGmIiAiEmAqAiIiAfb/nGrzBdaNxLoAAAAASUVORK5CYII=\n",
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ]
          },
          "metadata": {
            "needs_background": "light"
          }
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "PEWbON74YDvC"
      },
      "source": [
        "# Multiply Segments by Window Function (Hann)"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "VcSXMM5mWCnS"
      },
      "source": [
        "xgrid = xgrid*np.tile(w/alpha,[M,1])"
      ],
      "execution_count": 285,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 265
        },
        "id": "nfQrXVVvW8_j",
        "outputId": "095a9225-3367-4705-95ba-8eda96ede70e"
      },
      "source": [
        "plt.plot(xgrid[M-1,:], linestyle = '-')\n",
        "plt.show()"
      ],
      "execution_count": 286,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY8AAAD4CAYAAAAUymoqAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3dd3xW5f3/8dcngYACslFkheVApQgpIOJGRVCxaq21KrYqtWqXdoDiwoX21zqqXyvVVq1VatUqFawC7gESZA9lCBVkBNnDQJLP74/7QG8gjOS+k+se7+fjcT9yznXOufM+cCefnHVd5u6IiIhURE7oACIikn5UPEREpMJUPEREpMJUPEREpMJUPEREpMJqhA5QGU2aNPH8/PzQMURE0srkyZNXuXvTZLxXWhaP/Px8CgsLQ8cQEUkrZrY4We+l01YiIlJhKh4iIlJhKh4iIlJhKh4iIlJhKh4iIlJhSSkeZvYXM1tpZjP3sNzM7GEzm29m082sa9yygWY2L3oNTEYeERGpWsk68ngK6LuX5WcBHaPXIOAxADNrBNwG9AC6A7eZWcMkZRIRkSqSlOc83P09M8vfyyoDgGc81v/7BDNrYGbNgZOBse6+GsDMxhIrQs8nI5dIVdlUXAJAnVq7/wjNW7GB1Zu2ArBkzRYu6NayWrOJVIfqekiwBfBl3PySqG1P7bsxs0HEjlpo3bp11aQUKUdJaRk1cnPYWlLGYUNf54SOTXh/3ioAFg3vz5T/ruE7//cRL/2kF93aNOT0B97bafteHRqzfksJDQ+sSbODaofYBZGkS5snzN19BDACoKCgQCNYSZX783sLuXvMHAB+dHxbXiiM/Z2zvXAAPPH+Qu4aHVvngsc+4rh2jXd7nxXriznv0Q8BeOjiLvx85FQOOag2r15/PD3uGc+vzzyc607pUNW7I5JUlqyRBKPTVq+5+9HlLHsceMfdn4/mPyN2yupk4GR3/3F56+1JQUGBq3sSSaZNxSUMeXkGx7VvTN+jDqFhnTzyB49O+vdp17QOC4s2lbvsV2ccxmU986lTK5cauboRUpLPzCa7e0FS3quaikd/4HqgH7GL4w+7e/fogvlkYPvdV58C3bZfA9kTFQ9JJnen7ZAxO7WNv/EkTvv9u4ESxXw0+FQObXBA0AySWZJZPJJ1q+7zwMfA4Wa2xMyuNLNrzOyaaJUxwEJgPvBn4FqAqEjcCUyKXsP2VThEkmHd5m388h9TyR88mr98uGi35aELB8D5//cRy9Zt4fZRs1iyZnPoOCI7SdqRR3XSkYdUxuatJXy2fAPPfLyYf01ZuqO9Vo0cikvKAibbs+5tG/HJF7G/p567uge92jcJnEjSWTKPPNLmgrlIIs7+4/vMXLq+3GWpWjiAHYUD4G8fL6ZX+yaUlJaRm2OYWcBkku1UPCSjdbhpDCVl6Xd0XZ7XZy5n9PRlXPfcp1zRK5/bzz0qdCTJYrqlQzLWsnVbMqZwbHfdc58C8LcJsTF9iktKKc2wfZT0oOIhGaukNHN/qZaWOQ+Nm8fhQ//Dj/82mXVbtvHRglX73lAkSXTaSjLKkjWbmfLftazZvJVTDm8WOk6VemDc5wCMm7OCb93xJgB3nXc0l/ZsEzKWZAkVD8kYo6Z9xc+enxLXMitYllCGvjKTdz5bSetGdbj1nE6h40gGU/GQtPfo2/OpV7sGf3pnQegoKWHcnJUANK6bx9mdm9OmcZ3AiSQT6TkPSTvuzrQl6+jSqgFAlXQjkkn+fHkBp3c6OHQMSQEp94S5SHV6dsJiznv0Q/IHj2bVxuLQcVLemBnLQkeQDKTiIWmjuKSU/3tnPu/F9Wr75/cXBkyUHtyd9+cVsaBoY+gokkF02krSwupNW+l659jQMdJeiwYH8PovTuCg2jVDR5EAdNpKssr6b7apcCTJ0rVbmLDg69AxJAOoeEhKG/LydDrf/mboGBmlzJ2fPj+F/MGjeWPW8tBxJE2peEjKWrH+G57/5Mt9rygVcs2zn/LvaV8B8OO/TQ6cRtKVioekpKINxfS4Z3zoGFnhjn/P4oVJKtJSMSoekjLcnX9M+i9btpby7bvHhY6TNf764SJ+89L00DEkzSRrJMG+ZvaZmc03s8HlLH/AzKZGr8/NbG3cstK4ZaOSkUfS05gZy/ntSzM48tb/hI6Sla56ehLPTlhMcUlp6CiSBhLunsTMcoFHgdOBJcAkMxvl7rO3r+Puv4xb/6fAsXFvscXduySaQ9JbWZlTtOGb0DGy2rg5Kxk3ZyXPTljM0P6d6N1RoxbKniXjyKM7MN/dF7r7VmAkMGAv638feD4J31cyQFmZ89H8VbS7aQwvfrokdBwB5i7fwKVPTtwxZohIeZLRMWILIP5q2xKgR3krmlkboC3wVlxzbTMrBEqA4e7+yh62HQQMAmjdunUSYktom7eW0OnWN3bM72mYWAnjlldm0rhOHn2OPJi8Gro8Kjur7k/ExcCL7h5/UrVN9MTjJcCDZta+vA3dfYS7F7h7QdOmTasjq1ShktIy7nt9bugYsg/X/v1TDhv6OttKU3ecdwkjGcVjKdAqbr5l1Faei9nllJW7L42+LgTeYefrIZKBNm8tocPNr/P0xzotki666+432UUyisckoKOZtTWzPGIFYre7pszsCKAh8HFcW0MzqxVNNwGOB2bvuq1kDnff6VSVpIc1m7fh7qRjX3hSNRIuHu5eAlwPvAHMAV5w91lmNszMzo1b9WJgpO/86TsSKDSzacDbxK55qHhkqBXrv+HlT/d0UCqpru2QMQzSE+kSUa+6Ui0ufWIiH8xfte8VJeV9cW8/zCx0DKkE9aoraUeFI3MM/49udBAVD6liG77ZpmFiM8zj7y5k7eatlJU5pWXpd+ZCkkPFQ6rMZ8s3cIy6U89IM5eup91NY2h/05jQUSQQFQ9JurIyp6S0jDMffC90FKkilz45cce0OrHMTsl4wlxkhy1bS+l+9zg2FJeEjiLVpGhDMfNWbKDjwfVCR5FqpCMPSaojb/2PCkcWuvTJibwxazmffLE6dBSpJjryEJGErVhfvGNUwkXD+wdOI9VBRx6SFBuLSzhs6OuhY0gKOPeRD5ixZF3oGFLFVDwkKY6+7Q22lqjzPIHpS9ZxziMfsOGbbaGjSBXSaSuptG+2lfK7Nz6jRD2uSjlOuP9tPh16Ojk5eho9E6l4SKUdcYuGi5U9W7t5G+1uGsOr1x3Pt1o1CB1HkkynraRS/ln45b5XEgEGPPohK9driOFMo+IhFbZs3RZ+/eL00DEkjcQ/VCiZQaetZL+VljmvTf+KL1ZtCh1F0sznKzYyd/l6mtWrTaM6eaHjSBKoeMh+u+jxj5m8eE3oGJKm+j74PqDnQDJFUk5bmVlfM/vMzOab2eByll9hZkVmNjV6XRW3bKCZzYteA5ORR5Jv2pdrVTgkKQY8+mHoCJIECR95mFku8ChwOrAEmGRmo8oZEfAf7n79Lts2Am4DCgAHJkfb6rdUClGX6pJM075cS7+H3mdIvyM4oWPT0HGkkpJx5NEdmO/uC919KzASGLCf254JjHX31VHBGAv0TUImSQJ35/F3F4SOIRlo9rL1XPbkJ6FjSAKScc2jBRB/3+YSoEc5611gZicCnwO/dPcv97BtiyRkkiRoO0RjNUjV+mLVJto2qRM6hlRCdd2q+28g3907Ezu6eLqib2Bmg8ys0MwKi4qKkh5QdrZC9+VLNTjl/70TOoJUUjKKx1KgVdx8y6htB3f/2t2Lo9kngG77u23ce4xw9wJ3L2jaVOdJq9IzHy+ixz3jQ8eQLJE/eDTXPfdp6BhSQckoHpOAjmbW1szygIuBUfErmFnzuNlzgTnR9BvAGWbW0MwaAmdEbRJISWkZt746K3QMyTKjpy8jf/BoPUOURhIuHu5eAlxP7Jf+HOAFd59lZsPM7NxotZ+Z2Swzmwb8DLgi2nY1cCexAjQJGBa1STVzd371z2l0uFndqks4Oo2VPszdQ2eosIKCAi8sLAwdI6P0+cO7zF+5MXQMEX7/3W9x2pHNaHCgnkRPNjOb7O4FyXgvPWEuvPzpEhUOSRk3/nMaoCfRU506Rsxy67Zs44YXpoWOISJpRkceWWzwS9MZOUldq0tq+tO7C2hatxYXdGsZOoqUQ8Uji6lwSCob/vpcAIo2FnPNSe0Dp5Fd6bRVFnp/XhG/GDkldAyR/TL89bms3byVb7aVho4icXTkkWVenbqUn4+cGjqGSIV0GTYWgIk3ncbBB9UOnEZARx5ZR4VD0ln/h98PHUEiKh5Z4oVJX6prdUl7qzZupXCRniNOBSoeWeDdz4v4zUsac1wyw4V/+pj8waN1DSQwFY8M9/KnSxj4F42bIJmn8x1vsnKDen8ORcUjg32zrVQPAErG2lpSRve7x7N609bQUbKSikeG2vDNNv5ZqOc4JPN1vXNs6AhZSbfqZqDbR83iqY8WhY4hUm1+MXIKv7+oC7k5FjpK1lDxyCAlpWXqUl2y0itTv+KVqV/xm76Hc+3JHULHyQo6bZVBxsxcHjqCSFD3/+ez0BGyhopHBigpLSN/8Gh+9ry6HBHJHzya9+cVhY6R8ZJSPMysr5l9ZmbzzWxwOctvMLPZZjbdzMabWZu4ZaVmNjV6jdp1W9m7r9Zu0akqkV1c9uQnLCzSGDVVKeHiYWa5wKPAWUAn4Ptm1mmX1aYABe7eGXgRuD9u2RZ37xK9zkX22+yv1tNr+FuhY4ikpFN//y5L1mwOHSNjJePIozsw390XuvtWYCQwIH4Fd3/b3bf/L04A1EF/gr5au4V+6udHZK963/c2BXeNY1tpWegoGScZxaMFEP9AwZKobU+uBOLPs9Q2s0Izm2Bm5+1pIzMbFK1XWFSU3eczx85eoSMOkf20amMxHW9+nc1bS0JHySjVesHczC4FCoDfxTW3iQZkvwR40MzKHfXF3Ue4e4G7FzRt2rQa0qamt+eu5OpnCkPHEEk7nW59I3SEjJKM4rEUaBU33zJq24mZ9QFuBs519+Lt7e6+NPq6EHgHODYJmTKOu/OHNz/jh09NCh1FJG1d8NhHzFuxIXSMjJCM4jEJ6Ghmbc0sD7gY2OmuKTM7FnicWOFYGdfe0MxqRdNNgOOB2UnIlFE2FZfQdsgYHn5rfugoImlt8uI1nP7AeyxatSl0lLSXcPFw9xLgeuANYA7wgrvPMrNhZrb97qnfAXWBf+5yS+6RQKGZTQPeBoa7u4pHnLIy56jbdLgtkkwn/793KC5Rl+6JMHcPnaHCCgoKvLAw88/7T1+ylnMf+TB0DJGMNnloHxrXrRU6RrUws8nRNeaE6QnzFOTuPDRungqHSDXodtc4VqzXuCAVpeKRgtoOGcMD4z4PHUMka/S4ZzxfrtYDhRWh4pFC5i5fr3HGRQI54f63uelfM0LHSBsqHinij+Pn0fdBPTEuEtJzE/9L/uDRbCzWA4X7ovE8AispLaPnvW+xamPxvlcWkWpx9G1vcEmP1tzznWNCR0lZOvIIpKzMeeL9hXS4+XUVDpEU9NzE/3L1M4WUlqXfHanVQcUjgG2lZbS7aQx3jZ4TOoqI7MXY2Stof9MYxs9ZETpKytFpq2p266szeebjxaFjiEgFXPl07Lmy2cPO5MA8/doEHXlUm8JFq8kfPFqFQySNdbr1DV6dulvXfVlJT5hXMXen7ZAxoWOISJL96dKu9D26eegYFZLMJ8x1/FWFbnllJn+boCMNkUx0zbOfAjDttjOof0DNwGmqn4pHFfhg3ioufXJi6BgiUg2+dcebAMy4/Qzq1c6eIqLikUS3j5rFUx8tCh1DRAI45vY3OfGwpjx1xbfJybHQcaqcikeC1m3exreGvRk6hoikgPc+L6LdTWM4oWMTbjunEx2a1QsdqcqoeFTSl6s3c8L9b4eOISIp6P15q+jzh/d2zH9xbz/MMutoRMWjAlZv2krXO8eGjiEiaWb7HZc/P60jPzutI7kZcForKbfqmllf4CEgF3jC3YfvsrwW8AzQDfga+J67L4qWDQGuBEqBn7n7PofNq65bdbeWlHHt3yczbs7Kfa8sIlJBE4acxiH1a1fb90upW3XNLBd4FDgdWAJMMrNRuwwneyWwxt07mNnFwH3A98ysE7Exz48CDgXGmdlh7l7t40MWl5Qy5KUZvDxFDwCJSPXoee/43dpGDupJj7aNUv40VzJOW3UH5rv7QgAzGwkMAOKLxwDg9mj6ReARi/3LDABGunsx8IWZzY/e7+Mk5NrN4+8uYNZX63Fg2dotFC5eUxXfRkSk0i4eMWGPywZ0OZSb+h3JwQdV39HKniSjeLQAvoybXwL02NM67l5iZuuAxlH7hF22bVHeNzGzQcAggNatW1cq6MKiTUxfspYcMxau2lSp9xARCWXql2v5Zlu1n5gpV9pcMHf3EcAIiF3zqMx73Hdh593aiktKeXHyEm7+18zEAoqIJMn5XVvwqzMOp3n92il7+ioZxWMp0CpuvmXUVt46S8ysBlCf2IXz/dm2StWqkcsPerThBz3a7NReVuY89dEihr02ew9biogk5sC8XN668eRqvWieLMkoHpOAjmbWltgv/ouBS3ZZZxQwkNi1jAuBt9zdzWwU8JyZ/YHYBfOOwCdJyJSwnBzjR73b8qPebQH4emMxFz3+MQuKdLpLRCrvhR8fR/e2jULHSFjCxSO6hnE98AaxW3X/4u6zzGwYUOjuo4Angb9FF8RXEyswROu9QOzieglwXYg7rfZH47q1GH/jyUDsVFenW9/QCGMisl9e+kkvurVpGDpGUqlL9gRNXryGCx77KHQMEUkxRxxSj9d/fkJKXbNIqec8sl23Ng1ZNLw/G4tL6HLHm5ToaEQkq912Tid+eHzb0DGqnIpHktStVYP59/Rj3ZZtO7poFpHscf8Fnbno2632vWKGUPFIsvoH1GTR8P4A5A8eHTiNiFS1Ry45lv7HNE+p01PVQcWjCi0a3p/iklKu+/sUxs1ZETqOiCTRtFvPoP6B2TP4065yQgfIdLVq5PLEwAL6HHlw6CgikiSzh52Z1YUDdORRbZ4YWMDG4hK6DhvL1tKy0HFEpBIm3nQazerVyrpTVOVR8ahGdWvV4PO7z2Lm0nWc/ccPQscRkf30xOUF9OmkswfxdNoqgKNb1GfhPf1CxxCR/TDj9jNUOMqhI49AcnKMRcP7M2Hh13vtgllEwnj+6ti4GjkZMOpfVdCRR2A92zXmi3t1FCKSSl77aW+Oa99YhWMvVDxSgFnsKOSxH3QNHUUk631xbz+OblE/dIyUp+KRQs46pjlv3XhS6BgiWek7x7Zg4T39dCfVflLxSDHtmtZl3A0nho4hklXuv6AzD3yvi05TVYAumKegDs3qsWh4fz5asIpL/jwxdByRjLa9OyGpGB15pLBe7ZtwVe/M751TJIRbzu6kwpEAFY8UN/TsTvzzmuNCxxDJOFfqD7OEJFQ8zKyRmY01s3nR192GyjKzLmb2sZnNMrPpZva9uGVPmdkXZjY1enVJJE+m+nZ+I/2FJJIki4b3189TEiR65DEYGO/uHYHx0fyuNgOXu/tRQF/gQTNrELf81+7eJXpNTTBPRlt4Tz/+ckVSBgETyUoThpwWOkLGSLR4DACejqafBs7bdQV3/9zd50XTXwErgaYJft+slJNjnHqEukkQqYxxN5zIIfVrh46RMRItHge7+7Joejmw199sZtYdyAMWxDXfHZ3OesDMau1l20FmVmhmhUVFRQnGTm8L7ulHvdq6UU5kf829sy8dmtULHSOjmPvex9w2s3HAIeUsuhl42t0bxK27xt13u+4RLWsOvAMMdPcJcW3LiRWUEcACdx+2r9AFBQVeWFi4r9Uy3vvzirjsyU9CxxBJWX2ObMYTA78dOkbKMLPJ7p6Uc9/7PPJw9z7ufnQ5r1eBFVEB2F4IVu4h8EHAaODm7YUjeu9lHlMM/BXonoydyhYndGzKSz/RnVgi5RnQ5VD+fLmuEVaVRE9bjQIGRtMDgVd3XcHM8oB/Ac+4+4u7LNteeIzY9ZKZCebJOt3aNOLzu86iY7O6oaOIpIzLerbhoYuPVVcjVSjR4jEcON3M5gF9onnMrMDMnojWuQg4EbiinFty/25mM4AZQBPgrgTzZKW8GjmMveEk9c4rApzftQXDBhwVOkbGS+iqq7t/Dex275u7FwJXRdPPAs/uYftTE/n+sjP9lSXZrFf7xgw560iOaakecauDnjDPMJOH9gkdQSSI567uqcJRjVQ8MkzjurVYNLw/+Y0PDB1FpNr89Ye6o6q6qXhkqHd+fUroCCLV4vpTOnDK4c1Cx8g6Kh4ZbO6dfTmjk55Il8w1986+/OrMw0PHyEoqHhmsds1cRlxewA+Pzw8dRSTpZg87k9o1c0PHyFoqHlngtnN026JkltnDzuTAPHXRE5KKR5b44t5+HN3ioNAxRBL21o0nqXCkABWPLGFmvPbTE0LHEKm0K3rlM+WW02nXVL0ppAIVjywz986+5OXqv13Sy/NX9+T2c4+iYZ280FEkot8iWaZ2zVw+v/ssntJ98ZImbup3BMe1bxw6huxCxSNLnXx4M/X/I2nhqt7tQkeQcqh4ZLHLj8vndxd2Dh1DZI8e/v6x5OSoz7ZUpOKR5b5b0IpZd5wZOobIbj695XTO/dahoWPIHqh4CHVq1eAc/ZBKiqhdM4dFw/vTSBfHU5qKhwDwx+8fS+tG6kxRwvpWy/pMueWM0DFkPyRUPMyskZmNNbN50dc9jV9eGjcQ1Ki49rZmNtHM5pvZP6JRByWQ8TeexGM/6Bo6hmSxEZcXcECeuhxJB4keeQwGxrt7R2B8NF+eLe7eJXqdG9d+H/CAu3cA1gBXJphHElAzN4ezjmnOs1f2CB1FstCYn53AwQfVDh1D9lOixWMA8HQ0/TSxccj3SzRu+anA9nHNK7S9VJ3eHZuwaHj/0DEki/znFyfQ6VB1n5NOEi0eB7v7smh6ObCn/r9rm1mhmU0ws+0FojGw1t1LovklQIsE80gSvfSTXqEjSIYb2v9IFg3vzxGHqHCkm332LmZm44BDyll0c/yMu7uZ+R7epo27LzWzdsBbZjYDWFeRoGY2CBgE0Lp164psKpXUrU25l7BEkuaqE/QAYLra55GHu/dx96PLeb0KrDCz5gDR15V7eI+l0deFwDvAscDXQAMz217AWgJL95JjhLsXuHtB06ZNK7CLkog5w/pyWc82oWNIBpozrG/oCJKARE9bjQIGRtMDgVd3XcHMGppZrWi6CXA8MNvdHXgbuHBv20tYB+Tlcts5nULHkAxzRa983VWV5iz2O7ySG5s1Bl4AWgOLgYvcfbWZFQDXuPtVZtYLeBwoI1asHnT3J6Pt2wEjgUbAFOBSdy/e1/ctKCjwwsLCSueWinN3np34X1Zv3MoD4z4PHUfS2GM/6MpZxzQPHSMrmdlkdy9IynslUjxCUfEIK3/w6NARJE3pLr6wklk89IS5VNhzV+s5EKk49eKcWTSWo1RYr/ZNaNe0DguLNoWOImli/t1nUUODkGUU/W9KpYy/4aTQESRN/OWKAhWODKT/UakUM2Pe3WfRtXWD0FEkRXVp1YDP7urLqUfs6dlhSWcqHlJpNXNzeOknvThWBUTK8fhl3ahVQ7fjZioVD0mImfGva4/nrKPL64RAstW95x+jTg4znIqHJMVjl3Zj/t1nhY4hKeDe84/h+93VhVCmU/GQpKmRm8Mr1x0fOoYE9NDFXbj4261Cx5BqoFt1Jam6tNL1j2x1/wWdGdBFHWNnCxUPSbrtTxHPWLKOcx75IHAaqQ5v/vJEDju4XugYUo1UPKTKHNOyfugIUg0mDDmNQ+rr4ni20TUPqVJv/+rk0BGkirx140nMvbOvCkeWUvGQKtW2SR3+fHlS+mGTFHLXeUfTrmldatfUcxzZSqetpMqd3klPGGeS+y/ozEW6oyrr6chDqsWCe/rx6zMPDx1DEvSdY1uocAig4iHVJDfHuO6UDnRqflDoKFJJJ3Rswn0XdA4dQ1JEQsXDzBqZ2Vgzmxd9bVjOOqeY2dS41zdmdl607Ckz+yJuWZdE8kjqG3X98Xx6y+l01p1YaaN2zdiviScHfpu8Gvp7U2ISHYb2fmC1uw83s8FAQ3f/7V7WbwTMB1q6+2Yzewp4zd1frMj31UiCmWHz1hKOuf1NSsvSbzTLbDJyUE96tmscOoYkQSqNJDgAeDqafho4bx/rXwi87u6bE/y+kgEOzKvBgnv6hY4he3HngKPo0bZR6BiSghItHge7+7Joejmwr9tqLgae36XtbjObbmYPmFmtPW1oZoPMrNDMCouKihKILCL74/CD63HZcfmYWegokoL2edrKzMYB5fW3fTPwtLs3iFt3jbvvdt0jWtYcmA4c6u7b4tqWA3nACGCBuw/bV2idtsosC4o2ctrv3w0dQyLjbjiR1o3qUDPXVDgyTDJPW+3zOQ9377OXICvMrLm7L4sKwcq9vNVFwL+2F47ovbcftRSb2V+BX+1nbskg7ZvWZdHw/jw7YTFDX5kZOk5WO+Sg2nRopj6qZN8SPW01ChgYTQ8EXt3Lut9nl1NWUcHBYn/enAfoN0cW+0GP1hQO7cOALoeGjpK1PvjtKaEjSJpI9Anz4cALZnYlsJjY0QVmVgBc4+5XRfP5QCtg13MTfzezpoABU4FrEswjaczMaFK3Fp2aH8SrU78KHSdrDD7rCBrXyeOYlvWpkatbcWX/JFQ83P1r4LRy2guBq+LmFwG7dfTv7qcm8v0lM119QjtKypzibaU8/Nb80HEy3jUntQ8dQdKQ+raSlJMTPY0OkN+kDje8MC1wosx1Ve+2oSNImlLxkJR2fteW9O7QhBHvLeSJD74IHSdjzL2zLzlm1MzV3VRSOTrBKSmv2UG1GXp2J0Zc1i10lLSXY3Bl77bUrplLXo0c3YorlaYjD0kbZxx1CAfUzGXLttLQUdLWqOt7c3QL9SsmiVPxkLQy586+AOQPHh04SXo5v2sL7hxwNHVq6UdekkOfJElLb914Eh8t+Bp355ZXZ4WOk3Ka1qtF0YZiAN751cnkN6kTOJFkGhUPSUvtmtalXdO6ACoe5Rj7yxMZN2cltWvmqHBIldAFc9Euh8QAAAhESURBVEl7c6NTWfI/DQ7M48JuLTm7s57Wl6qh4iFpr3bNXAqH9uF3F2qUO0B3pUm1UPGQjNCkbi2+W9CKmXecuaPt/GN369QgY/U5shmTbu7D6J/15oyjyusEWyS5dM1DMkrdWjX4cPCprNm0lcZ183h5ytLQkarcx0NOpXn9A4DYhXKR6qDiIRmnRYMDaNHggNAxqtTQ/kfyn5nLObpF/R2FQ6Q6qXhIRptx+xks/nozX67ezLGtG9Lz3vGhIyXkkUuO5frnptCrfROuOqFd6DiSxfY5kmAq0kiCUlkzlqzjnEc+CB2jUj4afCqHZvgRlVStZI4kqAvmklWOaVmfRcP7s2h4f45J8W46nr+6JwA1coxFw/urcEhK0WkryVr//mlv3J1Ln5zIwOPyGfS3yaEj7XBQ7Roc174xc+/sS54GaJIUlNCn0sy+a2azzKwsGj1wT+v1NbPPzGy+mQ2Oa29rZhOj9n+YWV4ieUQqysz4+1U9OeOoQ7j25P8NitSrfeMgv7QXDe/Pwnv6MfXWM4DYMyw5Oer5VlJPokceM4Hzgcf3tIKZ5QKPAqcDS4BJZjbK3WcD9wEPuPtIM/sTcCXwWIKZRCrlN32P4McntuegA2rs6Ko82R0wxvc5BdC60YH8d/Vmrj25Pce1bwygYiFpIdFhaOcA+xoToDsw390XRuuOBAaY2RzgVOCSaL2ngdtR8ZCA6h9Yc6f5WXecydK1W1i9aSs92zVma0kZxSWlnPr7dynaUMy4G04EjFlfrePnI6cCcNJhTXn38yKa1M1j1catANSrXYNf9jmMU49oRm6OsXJDMe5Ol1YNNG64pKXquObRAvgybn4J0ANoDKx195K49j0+Emxmg4BBAK1bt66apCK7qFOrBocdXG/HfF6NHPJq5DDp5j47rdehWV1Wb9rKN9vK+MnJ7Rk9fRk92zWiVs1c6uTl7vYHVqtGB1ZLfpGqss/iYWbjgPL6O7jZ3V9NfqTyufsIYATEbtWtru8rsr9+ePz/xgPv37l5wCQiVW+fxcPd++xrnX1YCrSKm28ZtX0NNDCzGtHRx/Z2ERFJcdVxsnUS0DG6syoPuBgY5bGnE98GLozWGwhU25GMiIhUXqK36n7HzJYAxwGjzeyNqP1QMxsDEB1VXA+8AcwBXnD37aP3/Ba4wczmE7sG8mQieUREpHqoexIRkSyh7klERCQoFQ8REakwFQ8REakwFQ8REamwtLxgbmZFwOJKbt4EWJXEOOkkm/cdsnv/s3nfIbv3P37f27h702S8aVoWj0SYWWGy7jZIN9m875Dd+5/N+w7Zvf9Vte86bSUiIhWm4iEiIhWWjcVjROgAAWXzvkN273827ztk9/5Xyb5n3TUPERFJXDYeeYiISIJUPEREpMKyqniYWV8z+8zM5pvZ4NB5KsvM/mJmK81sZlxbIzMba2bzoq8No3Yzs4ejfZ5uZl3jthkYrT/PzAbGtXczsxnRNg/bPsYZrk5m1srM3jaz2WY2y8x+HrVn/P6bWW0z+8TMpkX7fkfU3tbMJkZ5/xENfYCZ1Yrm50fL8+Pea0jU/pmZnRnXntI/I2aWa2ZTzOy1aD6b9n1R9LmcamaFUVu4z727Z8ULyAUWAO2APGAa0Cl0rkruy4lAV2BmXNv9wOBoejBwXzTdD3gdMKAnMDFqbwQsjL42jKYbRss+ida1aNuzQu9z3H42B7pG0/WAz4FO2bD/UZ660XRNYGKU8wXg4qj9T8BPoulrgT9F0xcD/4imO0Wf/1pA2+jnIjcdfkaAG4DngNei+Wza90VAk13agn3us+nIozsw390XuvtWYCQwIHCmSnH394DVuzQPAJ6Opp8Gzotrf8ZjJhAbvbE5cCYw1t1Xu/saYCzQN1p2kLtP8Ngn6pm49wrO3Ze5+6fR9AZiY8S0IAv2P9qHjdFszejlwKnAi1H7rvu+/d/kReC06K/JAcBIdy929y+A+cR+PlL6Z8TMWgL9gSeieSNL9n0vgn3us6l4tAC+jJtfErVlioPdfVk0vRw4OJre037vrX1JOe0pJzoVcSyxv8CzYv+j0zZTgZXEfvAXAGs9Nuga7Jx3xz5Gy9cRG3Stov8mqeJB4DdAWTTfmOzZd4j9ofCmmU02s0FRW7DP/T7HMJf04+5uZhl9D7aZ1QVeAn7h7uvjT89m8v67eynQxcwaAP8CjggcqVqY2dnASnefbGYnh84TSG93X2pmzYCxZjY3fmF1f+6z6chjKdAqbr5l1JYpVkSHnkRfV0bte9rvvbW3LKc9ZZhZTWKF4+/u/nLUnDX7D+Dua4G3iQ0B3cDMtv8hGJ93xz5Gy+sDX1Pxf5NUcDxwrpktInZK6VTgIbJj3wFw96XR15XE/nDoTsjPfeiLQNX1InaUtZDYRbLtF8SOCp0rgf3JZ+cL5r9j5wtn90fT/dn5wtkn/r8LZ18Qu2jWMJpu5OVfOOsXen/j9tOInY99cJf2jN9/oCnQIJo+AHgfOBv4JztfNL42mr6OnS8avxBNH8XOF40XErtgnBY/I8DJ/O+CeVbsO1AHqBc3/RHQN+TnPvg/SjX/B/QjdnfOAuDm0HkS2I/ngWXANmLnJq8kdj53PDAPGBf3gTDg0WifZwAFce/zI2IXDOcDP4xrLwBmRts8QtQTQSq8gN7Ezv1OB6ZGr37ZsP9AZ2BKtO8zgVuj9nbRD/786Jdprai9djQ/P1reLu69bo727zPi7qpJh58Rdi4eWbHv0X5Oi16ztucL+blX9yQiIlJh2XTNQ0REkkTFQ0REKkzFQ0REKkzFQ0REKkzFQ0REKkzFQ0REKkzFQ0REKuz/AyOcTnGCP6wPAAAAAElFTkSuQmCC\n",
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ]
          },
          "metadata": {
            "needs_background": "light"
          }
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Kkh00LhgZJ-I"
      },
      "source": [
        "# Equation (7) - Compute DFT"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "kbc0dOE2XerG",
        "outputId": "5248ba77-50df-4dc0-e95e-834acba03e53"
      },
      "source": [
        "X = abs(sp.fft.fft(xgrid))\n",
        "\n",
        "np.shape(X)"
      ],
      "execution_count": 287,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "(19, 48000)"
            ]
          },
          "metadata": {},
          "execution_count": 287
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "NTu-d9gMZoZm"
      },
      "source": [
        "# Equation (8) - Compute Power Spectrum\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "LzFD4QiFZvXD",
        "outputId": "4e5252de-e8a4-4ef8-f444-a9996718227c"
      },
      "source": [
        "P = (X/N)**2;                      #power spectrum = square of amplitude\n",
        "\n",
        "np.shape(P)"
      ],
      "execution_count": 288,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "(19, 48000)"
            ]
          },
          "metadata": {},
          "execution_count": 288
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "CJVkfHtfZ7vQ"
      },
      "source": [
        "# Equation (9) - Compute single-sided power spectrum\n",
        "\n",
        "remove DC (0 Hz) component and \n",
        "frequencies above Nyquist frequency\n",
        "Fs/2 (index of Fs/2 = N/2+1), divide\n",
        "by noise power bandwidth"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "k3fXvmLQaGxG",
        "outputId": "f38e7028-9768-4970-bd64-48ca3900d091"
      },
      "source": [
        "Pss = 2*P[:,1:int(np.floor(N/2)+1)]\n",
        "\n",
        "np.shape(Pss)"
      ],
      "execution_count": 289,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "(19, 24000)"
            ]
          },
          "metadata": {},
          "execution_count": 289
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "hYPWfETYe0BR"
      },
      "source": [
        "# compute frequencies of DFT bins"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "M81uIj-Qe6SO",
        "outputId": "c0b99f89-d8e9-4088-8392-0ea57ef952cc"
      },
      "source": [
        "f = np.floor(Fs/2)*np.linspace(1/(N/2),1,int(N/2))\n",
        "\n",
        "f\n",
        "\n",
        "np.shape(f)"
      ],
      "execution_count": 290,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "(24000,)"
            ]
          },
          "metadata": {},
          "execution_count": 290
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "9BtQUyxKN48v",
        "outputId": "9470a29c-ff58-4da4-d98b-38c3dff54391"
      },
      "source": [
        "flow = np.where((f >= lcut) & (f<=hcut))\n",
        "f = f[flow]\n",
        "nf = len(f)\n",
        "nf"
      ],
      "execution_count": 291,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "24000"
            ]
          },
          "metadata": {},
          "execution_count": 291
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "eoaAbRidRqer"
      },
      "source": [
        "# Compute noise power bandwidth and delta(f)"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "Tg3NEW5GRuv4",
        "outputId": "7dd3d799-3bf5-4f86-b969-511a94569bb8"
      },
      "source": [
        "B = (1/N)*(np.sum((w/alpha)**2))\n",
        "delf = Fs/N\n",
        "\n",
        "print(B)"
      ],
      "execution_count": 292,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "1.5\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "aRUYC4uEcXrQ",
        "outputId": "8a3766e3-d8b3-4812-c2de-c44c34117e4a"
      },
      "source": [
        "print(flow[0])"
      ],
      "execution_count": 293,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "[    0     1     2 ... 23997 23998 23999]\n"
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "myvoa-8JSRv1"
      },
      "source": [
        "# 1/3 octave analysis\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "k6a8RCwCSWnI"
      },
      "source": [
        "if lcut < 25:\n",
        "  lcut = 25\n"
      ],
      "execution_count": 294,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "xeYt29kPi043"
      },
      "source": [
        "#(Equation 13) Build TOL frequency"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "fIaq8zDehOSa",
        "outputId": "2e78b2c7-d71d-40ca-eadf-e88ff6161a40"
      },
      "source": [
        "  lobandf = np.floor(np.log10(lcut)) #lowest power of 10 frequency for 1/3 \n",
        "  # octave band computation\n",
        "  hibandf = np.ceil(np.log10(hcut));  #highest \"\"\n",
        "  nband = 10*(hibandf-lobandf)+1;     #number of 1/3-octave bands\n",
        "  fc = np.zeros((1,int(nband)));      #initialise 1/3-octave frequency vector\n",
        "  fc[0][0] = 10**lobandf;             #lowest frequency = lowest power of 10\n",
        "\n",
        "#get the center frequencies\n",
        "  for i in range(1,int(nband)):\n",
        "    fc[0][i] = fc[0][i-1]*10**0.1\n",
        "\n",
        "fc = np.squeeze(fc)\n",
        "fc = fc[np.logical_and(fc>lcut, fc <=hcut)] #crop to cutoff freq's\n",
        "\n",
        "nfc = len(fc) #number of TOL bands\n",
        "\n",
        "print(nfc)\n"
      ],
      "execution_count": 295,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "30\n"
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "P0Xy1tQq0MOL"
      },
      "source": [
        "#(Equations 14-15) Calculate Boundary Frequencies for each band"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "m6a7kRxx0dXq",
        "outputId": "9ef960bc-624a-4106-8579-cf803739564b"
      },
      "source": [
        "fb = fc*10**-0.05 #calculate lower bounds from Fc\n",
        "fb = np.append(fb,fc[nfc-1]*10**0.05) #calculate upper bounds from Fc\n",
        "print(hcut)"
      ],
      "execution_count": 296,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "24000\n"
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "PNutx3S33BQu"
      },
      "source": [
        "#(Equation 16) Calculate 1/3-octave band levels"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "lv3qMiIH3Eyc"
      },
      "source": [
        "P13 = np.zeros((M,nfc)) #create the array to hold the TOL data\n",
        "\n",
        "for i in range(0, nfc): #loop through center frequencies\n",
        "  fall = np.logical_and(f >= fb[i], f < fb[i+1]) #index freq's falling into band\n",
        "\n",
        "  for q in range(0,M): # loop through DFT's of data segments\n",
        "    fcl = np.sum(Pss[q,fall]) # integrate over m'th band frequencies\n",
        "    P13[q,i] = fcl #store TOL value\n",
        "\n",
        "\n",
        "if np.all(10*np.log10(P13[0][:]/(pref**2))  <= -10**6):\n",
        "  lowcut = np.where(10*np.log10(P13[0][:]/(pref**2))  <= -10**6)\n",
        "\n",
        "  P13 = P13[:,lowcut[np.size(lowcut)-1]:nfc-1]\n",
        "  fb = fb[lowcut[np.size(lowcut)-1]:nfc]\n",
        "  fc = fc[lowcut[np.size(lowcut)-1]:nfc-1]\n",
        "  nfc = len(fc)\n",
        "\n",
        "a = 10*np.log10((1/B)*P13/(pref**2))-S"
      ],
      "execution_count": 297,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "tfJiRQgZTz7j"
      },
      "source": [
        "# Compute Time Vector"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "DPXGW4DGT4-q",
        "outputId": "92a0d470-f9bb-4f0b-d09f-dc6f5a26f2a7"
      },
      "source": [
        "tint = (1-olap)*N/Fs\n",
        "ttot = M*tint-tint\n",
        "t = np.linspace(0,ttot,int(ttot/tint+1))\n",
        "t"
      ],
      "execution_count": 298,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "array([0. , 0.5, 1. , 1.5, 2. , 2.5, 3. , 3.5, 4. , 4.5, 5. , 5.5, 6. ,\n",
              "       6.5, 7. , 7.5, 8. , 8.5, 9. ])"
            ]
          },
          "metadata": {},
          "execution_count": 298
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "rZ90P_J8VaRB"
      },
      "source": [
        "# Construct and output array"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "L4oHRPJXVi8X"
      },
      "source": [
        "#a = np.ndarray.flatten(a)\n",
        "#A = np.column_stack((f,np.around(a,1)))\n",
        "\n",
        "A = np.zeros((M+1,nfc+1))\n",
        "\n",
        "A[1:M+1,1:nfc+1] = a\n",
        "\n",
        "A[0,1:nfc+1] = fc\n",
        "A[1:M+1,0] = t\n",
        "f = fc"
      ],
      "execution_count": 299,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "2TApbN5CfJba"
      },
      "source": [
        "# Plot the Data"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "IO9mfytnfNu-",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 295
        },
        "outputId": "4ca85db3-297c-485e-9890-97f54bf1c56e"
      },
      "source": [
        "import matplotlib.pyplot as plt\n",
        "\n",
        "z_min, z_max = a.min(), a.max()\n",
        "\n",
        "fig, ax = plt.subplots()\n",
        "\n",
        "\n",
        "c = ax.imshow(a.transpose(),cmap ='magma', vmin = z_min, vmax = z_max,\n",
        "                 extent =[t.min(), t.max(), 0, len(f)],\n",
        "                    interpolation ='none', origin ='lower', aspect='auto')\n",
        "ax.ticklabel_format(useOffset=False)\n",
        "\n",
        "ax.set_yticks(np.arange(0,30,3))\n",
        "ax.set_yticklabels(np.round(f[(np.arange(0,30,3))],1))\n",
        "\n",
        "plt.colorbar(c,label = '1/3-Octave SPL [dB re 20 $\\mu$Pa]')\n",
        "\n",
        "#plt.gca().invert_yaxis()\n",
        "\n",
        "plt.xlabel('Time [s]')\n",
        "plt.ylabel('Frequency [Hz]')\n",
        "plt.title('WhiteNoise_10s_48kHz_+-0.5.wav')\n",
        "\n",
        "plt.show()"
      ],
      "execution_count": 300,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaAAAAEWCAYAAAAgpUMxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO2dd5hlRbW339/kQGYAh4xKEFDSCKgXL0EUgStiQFDEzEXRi1kMV1Ef7sUsZkdFwMAFCYoICB+SRAkDDDAECZKGKGly7P59f1QdONOc7t6n+nSfc2bWO89+Zu/atapq77N7r11Vq9aSbYIgCIJgpBnV7gYEQRAEqyahgIIgCIK2EAooCIIgaAuhgIIgCIK2EAooCIIgaAuhgIIgCIK2EAooCIIgaAuhgLocSZb04n7OvUPSRSPdprr6L5D0rnbV3ylIuk/Sa9rdjiDoNEIBdRiSPivpgj5pd/WTduhAZdn+je3X1sn0q6z6act9kh6XNLku7f2SLqsib/v1tk+pWl8zSNpL0qWS5ki6r8H5zfP5hZLuaJUCkHREvo/vr0sbL+knkh6T9JSkP0raqGJ5z1NOkt4t6a+taG+zSNon36+F+f5tNkDe+yQtkjQ/b2372Am6k1BAnccVwCsljQaQNBUYC+zUJ+3FOe9wMxo4ZgTqaZYFwEnAp/o5fxpwI7Au8HngTEnrDaVCSWsDnwNu7XPqGOAVwMuADYGnge8Ppa5WIWnPqh8MkqYAZwP/DawDzABOH0TsP2yvlrfXDpI3CFYgFFDncR1J4eyYj/cALgX+0SftHtsP5+PX5B7RM5J+KEmw4pe0pJqyuil/rb4tpx8oaWaW/Zukl/VpzzeAT0paq1FjJb1S0nW5J3KdpFfWnbus1lOQ9GJJl+d8T0g6vS7fNpIuzr2Hf0g6ZLCbZPta278C/tmgTVsBOwNfsr3I9lnALcCb8/ldJc2QNDf3Wr49WH2Z/wW+BzzRJ30L4M+2H7O9mPTS3q5RAZJeIuleSYdVqVDS2+p6GPMlLamqUAp4E3Cr7d/l6zgO2EHSNkMpNPdWb6k7vljSdXXHV0p6Y94/VtI9kuZJuk3SwTl9fH5Gt6+TWy/3wNYfSvuC9hEKqMOwvRS4Bnh1Tno1cCXw1z5p9b2fA4GXk77ADwFe16DcmuwO+Wv1dEk7kXoR/0nqKfwUOFfS+DrRGcBlwCf7lilpHeBPpJfyusC3gT9JWrfBpX0VuAhYG9iY3EPIw3sXA78F1gcOBX4kadsGZVRlO+CftufVpd3Ec0rhROBE22sALwLOGKxASbsC04CfNDj9C+BVkjaUNAl4B3BB30ySdgb+DHzE9mlVLsT26bUeBql39U9S72442I50n2p1LwDuoR9lmvmNpH9JukjSDv3kuRrYUtIUSWPJPUVJq0uaSLqvV+a895A+sNYEvgz8WtJU20tIvbN6xX0IcLntx5u+0qAjCAXUmVzOc8pmD9If55V90i6vy3+C7WdsP0DqLe1INY4Efmr7Gts9eb5mCbB7n3xfBD7SYAjrAOAu27+yvTy/VO8A/qNBXcuAzYANbS+2XZvjOBC4z/Yvcxk3AmcBb614DY1YDZjTJ20OsHpdW14saYrt+bavHqiwPPT5I+DDtnsbZLkLeBB4CJgLvAT4Sp88ewDnAkfYPq/Pud/nr/tnJD2T6+rbhlEkJX2Z7Z8O1N4hMNh968s7gM1Jv+ulwJ8b9ZRtLyL17F8N7EJSclcBryI9a3fZfjLn/Z3th2332j6ddG93zUX9lvSBUuPtOS3oUkIBdSZXAP+Wexjr2b4L+BtpbmgdYHtW7AE9Wre/kPQiqcJmwCf6vPw2IX1pP4vtWcB5wLF95DcE7u+Tdj/QaAL+04CAayXdKum9dW3YrU8b3gG8oOI1NGI+sEaftDWAWo/ofcBWwB152PDAQcr7EHDzAIrqh8B4Ui9wMulLvW8P6Cjgb7YvayD/Rttr1bZcX1+OJymC/xqkrbVhrNq9PI/0LNXfXyRtWj+0l0UHu28rYPuqPMS50Pb/As+QFG0jLgf2JCmhy0m96n/P27MfU0pGHjPr2ro9MCWfvhSYJGk3SZuTPrTOGex+BJ1LKKDO5O+kIYgPkL4UsT0XeDinPWz73hbU8yBwfP3Lz/akfoaHvpTrrlcuD5MUSD2bknoCK2D7UdsfsL0hacjvR0oWeQ+ShlHq27Ca7Q8O4bpuBV4oqf7LfYecju27bB9GGvL7GslAYfLzi3mWfYCDJT0q6VHglcC3JP0gn98RONn2U3mo6PvArkqT+jWOAjaV9J1mL0bJ2vEw4C22lw2W3/YJdcrsQOCvfRQcth+oMx6ofbDcSrpPtXonk4Yo+xpd9Fs16SOjEX0V0OX0UUBKFnc/Az4MrJvbOqtWpu0e0nDpYXk7r88wa9BlhALqQPKQxQzg4zw3Ng5pHujjlFu/PQa8sO74Z8BR+YtSkiZLOqDPi7vWprtJk+v1X+DnA1tJerukMUqGDduSvrpXQNJbJW2cD58mvax6c96tJL1T0ti8vVzSSwa6EEmjJE0gGWxI0gRJ43Jb7wRmAl/K6QeT5h3OyrKHS1ovD6c9k4tsNLRW492kYbUd8zaDND/x+Xz+OuAISWvmOY4PkT4S6o0V5gH7Aa+WdMJA19bnOnciKbQ32v5XVblCzgG2l/TmfG+/SOr53dGgXZtKepWkcfkef4rUU7mqn7L/BmxNGk671vat5N4vzz3Pk0nPxb9yHe8h9YDq+S3wNlIvOYbfupxQQJ3L5aQv9Pr1IFfmtFIFdBxwSh7eOMT2DFKv5gckpXA36WXbH18hvSQAyOP2BwKfAJ4kDbMd2OfFW+PlwDV5uOdc4BjbNUOB15LG9h8mDSd+jTSkNRCvBhaRlOCmeb9+HcqhpMntp4ETSL2H2gt8P+DW3JYTgUOz0m9Inl97tLYBS4G5tmvzJZ8EFpPmK/4F7A8c3KgcYF/g9ZK+Osj11TiIZLjx17ohs+cZOLSCfH/eTBrue5qkHJ6dc1Fa61Qzwlgd+HHO9xDpnr6+NpejtAj61rqyFwA3kKzslubkvwP314wIbN8GfCunPwa8lD4KzfY1JBP8DWlg6BF0F4qIqEEQBEE7iB5QEARB0BZCAQUdS7aWm99ge8cw1PWOfuqqOgE/oozkvQm6C0knKbnQmlWX9tb8zPRKmlaXvrnSYt6ZeWu0zm342hpDcEEQBCsPkl5NMqk/1fb2Oe0lJEObnwKfzPO/ZHP282r5Rpox7ai005kyZU1vtvlQlqGMEKXfDv0Zyg4nI/2dU3qN3XJPR7qdQ/n92lHnCHLDDXc+YXtIfgZf97pd/eSTfdcAN+b66+/8s+39+jtv+4qsWOrTbgeQ2vHH3z+hgBqw2eYv4JprhmuxeQvpHchyeABGtWHktbStpZReY7fc05Fu51B+v3bUOYKMHbdP38XYTfPkk3O45tpq75wxo/faRtKMuqTptqcPofotJN1I8uLxBdtXDibQKkIBBUEQtBvTjMJ9wva0wbNV4hFgU9tPStqF5BZqu7zwfdgJBRQEQdB2DMuXj3ytyXPHkrx/vaR7SG6qZgwo2CLCCi4IgqDdGLCrbS1EKaRFLc7YC4EtaRDiZLiIHlAQBEHbccvmvCSdRvK7N0XSbJIfx6dILp3WI4VMmWn7dSSPIl+RtIxkJXeU7ada0pAKhAIKgiDoBFqkgLKj3UY8z3O4U7DGs1pScQGhgBqgnh40t5pJ5AqMKjRx7C3sVo90fe2gW6zSRvq3GGkLsTGjy2WX95TJdctv2AqaM0JYaQgFFARB0HZaNwTXTYQCCoIgaDc26hl5K7h2EwooCIKgE4geUBAEQTDimO6am20RoYCCIAjaTswBBUEQBO0grOCCIAiC9mAII4QAgJ4eNKf5dUAeU3Y7NdI+oIayvqJbvDB3S32rAoX3tGv+nlpBzAEFQRAE7SHmgIIgCIJ2EQooCIIgGHEMWgUV0LANXEs6SdLjkmbVpX1D0h2SbpZ0jqS1cvrmkhZJmpm3n9TJHCbplixzoaQpOX0HSX/P5/4oaY0GbZgg6VpJN0m6VdKXh+t6gyAIyqkYiqHF4RjazXDOnJ4M9I1bfjGwve2XAXcCn607d4/tHfN2FICkMcCJwF5Z5mbgwzn/z4Fjbb+U5OX1Uw3asATY2/YOwI7AfpJ2b8nVBUEQtAqTAtJV2VYihk0B2b6CFIOiPu0i27U7eDWw8SDFKG+TJQlYA3g4n9sKuCLvXwy8uUEbbHt+Phybt5XrEyIIgpUAJyu4KttKRDvngN4LnF53vIWkG4G5wBdsX2l7maQPArcAC4C7gKNz/luBg4DfA28FNmlUSY72dz3wYuCHtq/pJ9+RwJEAm05dB+bNb5RtQIqdubtw7Hfs2DK5oTzEIz1O3S3mzaX3pfS3V+F96Zb7yRD+nrpxLmUVXYjalqdR0ueB5cBvctIjwKa2dwI+DvxW0hqSxgIfBHYCNiQNwdWG7d4LfEjS9cDqwNJGddnusb0jqbe1q6Tt+8k33fY029PWW3u1llxnEARBZXp7q20rESPeA5L0buBAYB87zajZXkKar8H29ZLuIQ2xKafdk2XPAI7NaXcAr83pWwEHDFSv7WckXUqal5o1UN4gCIKRZeUzMKjCiPaAJO0HfBp4g+2Fdenr5aEyJL0Q2BL4J/AQsK2k9XLWfYHbc7718/+jgC8Az1rO9Sm3Zmk3McvfMTxXFwRBUIhJkWOrbCsRw9YDknQasCcwRdJs4Euk4bPxwMXJpoCrs8Xbq4GvSFoG9AJH2X4ql/Nl4Ip87n7g3bmKwyTV5oPOBn6Z828I/Nz2/sBU4JSs3EYBZ9g+b7iuOQiCoIzwhNBSbB/WIPkX/eQ9Czirn3M/oUHvxvaJJBPtvukPA/vn/ZtJ80dBEASdyypqhBCeEIIgCDqBlczEugqhgBrR24PmN2+GXfwAlZpTz1swsvXByLuML72no0e3th2DsaShEebgjBvhP8HRhfUtW1Ze56hCg+qRfiEP5e9iyLjcJL+L6Z5FAUEQBCsrtXAMLViI2o8btLdmd2S9kqb1yf9ZSXdL+oek17X+4vonFFAQBEG7aa0V3Mk83w3aLOBNPOc9BgBJ2wKHAttlmR/VLJJHghiCC4IgaDuts4KzfYWkzfuk1Zav9M1+EPB/eS3mvZLuBnYF/t6SxgxCKKAgCIJOoPqc1xRJM+qOp9ueXljrRiS/nDVm57QRIRRQEARBuzHNGCE8YXva4Nk6n1BAQRAEbadtnq4fYkVHzhvntIZIWqdCmb22n6lSeSigRvT0wpwCM+xSxhTO+ZU+sL2F5tsw8m0tNacuNRteXjgOP6HQhHfR4jK5UpNhLyqTG8rLsdQMe2mhyX9pfe2mPQroXJLz52+THD5vCVw7QP6H8zbQTR4NbFql8lBAQRAE7aZmBdcC+nGD9hTwfWA94E+SZtp+ne1bs5Pn20gRCo62PVBDbs9RCwaq/8aqbQ0FFARB0G7cuiG4ftygQYoc3Sj/8cDxFYt/RYvyAKGAgiAIOoMu8AVn+9kxY0lrk4bsJtSdv6I+z2CEAgqCIOgEuigekKT3A8eQjBZmAruT1g7t3Uw57YqIupakMyXdIel2Sa+Q9FVJN0uaKemiHFYBJb6XXUXcLGnnunK+JmlW3t42QH2HSLotu6L47UhcYxAEQWVa6IpnhDgGeDlwv+29SFEHKlm+1dOuHtCJwIW23yJpHDAJuNX2fwNI+i/gi8BRwOtJ3bwtgd2AHwO7SToA2BnYkRRj6DJJF9ieW1+RpC1JcYheZfvpWiC7IAiCzsHdFmxuse3FkpA03vYdkrZutpB2hORekxSA7t0AtpcCfV0JTyZ9E0ByFXFqDt99de49TQW2Ba6wvRxYLulmki+jM/qU9QHgh7afzvU93vqrCoIgGAK1HlD3MDtHm/49KcDo06SAoU3Rjh7QFsC/gF9K2gG4HjjG9gJJxwNHAHOAvXL+jYAH6+RrriJuAr4k6VukHtReJFPCvmwFIOkqkn36cbYv7JtJ0pHAkQCbTlkd5hWslyh9gMYUjoQWrwMawoNe2tbitSCFX4XjCtcPla4DKl13VPxbFK4fKlxb4yE8MypdI1X6zCwu/C1Gtdk3cxcpINsH593jJF0KrAk87706GO2442NIQ2c/zvbkC4BjAWx/3vYmwG+ADw9UiO2LgPOBvwGnkSbAGr2txpCG7/YEDgN+ljV33/Km255me9p6a0wqvLQgCIICXHH+p81KStIESR+V9ANJ/ylpjO3LbZ+bR7Oaoh0KaDYw2/Y1+fhMkkKq5zfAm/N+v64ibB9ve0fb+5JW5t7ZT33n2l5m+96cZ8uWXEkQBEGLcK8rbW3mFGAacAtpfv5bQylsxBWQ7UeBB+smrPYBbsvGAjUOAu7I++cCR2RruN2BObYfkTRa0roAkl4GvAy4qEGVvyf1fpA0hTQk988WX1YQBMHQsKtt7WVb24fb/inwFmCPoRTWLiu4jwC/yRZw/wTeA/w8K6Ve0mTWUTnv+cD+wN3AwpwXYCxwZY5vMRc4PBskIOkrwAzb5wJ/Bl4r6TbSEN2nbD85/JcYBEFQEVM+/ziyPDvBZnt5g/hCTdEWBWR7JqkbV8+b+8lr4OgG6YtJlnCNZL7YR/7jeQuCIOg8WuiKZ5jZQdJcnnNGOrHu2LbXaKaw8IQQBEHQCXSBArLd0nDdoYAa0dOL5yxsXq7UbHSku96lptRQ3tah1FlCacSJ0pfAxHFlcssKQw6UmqeXUmpKDXjBkjLB0oWZpSFDGhrRjhwdYGBQGUm3ADfXbbcA78qOTSvTZsP3IAiCoAtd8fw78DNgEXAoMIs0V98U0QMKgiDoBDpHuQyK7aeAy/JWc3n2hWbLCQUUBEHQbuwUiblLkLSV7WfXXdq+Ky+HaYpQQEEQBG3GgLtH/wD8VNKLSE4BbibFBJolaZLtyhPooYCCIAjaTZc5I80hGJC0KbADKSrBDsBMSb22t6lSTiigIAiCTqCLFFAN2w8ADwB/rKVJWq2qfCigRvQYzyk0HR1JCm0Yvbjc3FQTCk1cS4cXxpSZtntpWYWaWPgnMa/weRlpT+ilLCr0MD0ESp/T4md0eXsVQJcNwfWL7flV84YCCoIgaDem7QqwHYQCCoIgaDfuCE/XI04ooCAIgk6gRUNwkk4CDgQet719TlsHOB3YHLgPOMT205L2BP4A3JvFz7b9lUHK34YUsWCjnPQQKeTN7c22dcQ9IeSARtdKuknSrZK+nNM/LOluSc5hE+pl9pQ0M+e/PKdtIulSSbfl9GP6qe8gSTdn+RmS/m34rzIIgqBJXHEbnJOB/fqkHQtcYntL4JJ8XOPKHFdtxwrK5zPA/5Gcj16bNwGnSTp2INlGtKMHtATY2/Z8SWOBv0q6ALgKOI+8srZGjl76I2A/2w9IWj+fWg58wvYNklYHrpd0se2+YbkvIWln54VSZwCVTASDIAhGBLfOF5ztKyRt3if5IHJcNFJQucuAzxQU/z5gO9srWKVI+jZwK3BCM4W1IyCd66wkxubNtm+0fV8DkbeTuoUPZPnH8/+P2L4h788Dbue5LmF9ffNzSAaAyVT9hgiCIBhJeituMCWP5tS2IyuUvoHtR/L+o8AGdedekUekLpC0XYVWbtggfSoFg4htmQOSNBq4Hngx8MO68NyN2AoYK+kyYHXgRNun9ilvc2AnoGE5kg4G/hdYHzhgsPa5x/TOG0Gz09KYTl2kSr2k1Cy6zKS2d16Zl2ktLPROXegJXSNtZl5o9j2Ur3MvK2zr2MK2LipcZjB6aMHVhoTB1R+9J2z3jadWvao0GlT7QW8ANssjUvuTIkhv2b80HwUukXQX8GBO25T0Lv9Is21pV0C6HmDHPLx2jqTtbc/qJ/sYYBdS6O6JwN8lXV3zQ5QXPZ0FfNT23H7qOyfX82rgq8Br+ubJXxFHAmy65uQhXV8QBEEzjIArnsckTbX9iKSpQG0k6dl3pu3zJf1I0hTbTzRsp32hpK2AXVnRCOG6/F5viraGY7D9DHApz58wq2c28GfbC/JNuYLk8oE8h3QW8BvbZ1eo7wrghX2NHPK56ban2Z42ZdKEgqsJgiAoxDQzBFfCucC78v67SJZvSHqBclxtSbuSdMKTg5S1FWk648+2z7J9te0eSQO9xxvSDiu49XLPB0kTgX2BOwYQ+QPwb5LGSJoE7Abcnm/aL4DbbX97gPpeXHeDdwbGM/gNDoIgGFHsattgSDoN+DuwtaTZkt5HMg7YNw+dvYbnjAXeQnIiehPwPeDQujnzRmX/F+md/JEsd1Dd6f9p9prbMQQ3FTglzwONAs6wfV6+sE8DLwBulnS+7ffbvl3ShSSPq73Az23PyubU7wRukTQzl/253I08CsD2T4A3A0dIWkYKnvS2gW5wEARBO2jVEJztw/o5tU+DvD8AftBE8R8AdslzRpsDZ0ra3PaJFMxmj7gCsn0zyWCgb/r3SBq4kcw3gG/0Sfsr/VxwVjy1/a8BXxtCk4MgCIaX2hBc5zOqZsVs+768kPVMSZtRoIAiJHcQBEEH0NtTbWszj0nasXaQldGBwBTgpc0WFq54giAI2o2B3jaagVfnCJITgGexvZw0zfHTZgsLBdQA95jlTzf/qaFCL/Cla0FK+6+9S8qnwKSyT7DSWTfNLXTJX3hPexYUrlkp/Ak1rnAdUOFypd4lZfdz1Pjyl2P581b2W5S2VW0cD+qWiKi2Zw9w7qpmy+tXAUn6eAX5Bbab1npBEARBPcLuih5QSxlI538KWI3kfaC/7RPD3cAgCIKVHqceUJVtZWKgIbhfVfCMGi4DgiAIhoiB3p7oAT2L7U8PJlwlTxAEQTAIBveq0tYJKHG4pC/m402zJ4WmGHTaTVKPpBNq3gRy2g3NVhQEQRD0T6s8IYwQPwJeAdQWvc4DfthsIVXsPm7N+S7KUfWg3H9zEARB0ABblbYOYTfbRwOLAWw/DYxrtpAqZtjLbX9a0tuAKyUdQVcFAiigF3oWFfzQowpvy4hPLJY/xKPGll1j77JSU/Oy+p7zNt8cpV+YowpN8Jc9XXZfRhf+DqVLBZYuKJMDGNX0aylTONzUu6SsulHj2/ta65ThtYosy+7UDMnHJ8MUD0gAtk+XdCvwW1L8hyAIgqAFdNjwWhW+B5wDrC/peJJT0y80W0gVBfT+2k52AroHKbxrEARB0BJET093eEbL9gBXkIKK7kPqpLzR9u3NljXQQtQ31e1v1uf0fIIgCILW0EU9oBxR9XzbL2XgUDqDMpDK/Y+6bXqf4wOHUqmkkyQ9LmlWXdpxkh6SNDNv++f0sZJOkXSLpNslfbZPWaMl3SjpvH7q2kzSJZJulnSZpI2H0vYgCIJWY7rOCOEGSS8faiH99oBsv6e2L+nG+uMWcDIpBsWpfdK/Y/ubfdLeCoy3/dIckO42SafZvi+fPwa4HVijn7q+CZxq+xRJewP/S4ojFARB0DF0kHKpwm7AOyTdDywgDcPZ9suaKaSqM9KWdg5tX5GDGVWte7KkMcBEYCkwFyD3Zg4Ajgf68123bd25S4Hfl7U6CIJg+OjtLgX0ulYU0mnesD+czbxnAJ/ItuVnkoweHgEmAR+z/VTO/11SFNXVByjzJuBNwInAwcDqkta1vUJYbklHAkcCbDx5MssWN/8wLFlQZm86bvzIepgeCqVen4vNm8eW2aiPHlNoLl7oDqVnadkEcun1LZ4/tkiu+Fkbionw4jIxlS5rKERD8BI/VGx1lSse2/e3opx+/2ok/VHSuZLOBV5Y269LazU/Bl4E7EhSNt/K6bsCPcCGwBbAJyS9UNKBwOO2rx+k3E8C/y7pRuDfgYdyeStge7rtabanTRk/viUXFARBUJVeq9K2MjFQD6h+LuZb/eZqEbYfq+1L+hlQMyp4O3Ch7WXA45KuAqaRwnq/IRsrTADWkPRr24f3KfdhUg8ISasBb7b9zHBfTxAEQTN02RxQSxjICOHykWyIpKm2H8mHBwM1C7kHgL2BX2Xv27sD37V9BvDZLLsn8Mm+yiefmwI8Zbs35z9pWC8kCIKgSUzXzQE9S37HPmk3P9A+0BDc9AoVD5qnH7nTgL8DW0uaLel9wNezqfXNwF7Ax3L2HwKrZS8M1wG/tH3zIOV/RdIb8uGewD8k3QlsQDJYCIIg6BzcHWbYknbPy1nOlrRTXkozC3hM0n7NljfQENwbJQ00fSiSomga24c1SP5FP3nnk0yxByrvMuCyuuMv1u2fSTJkCIIg6Fha5RJS0kmktZqP294+p60DnA5sDtwHHGL76ezV4ERgf2Ah8G7bA0U7+AHwOWBN4C/A621fLWkb4DTgwmbaOpAC+lQF+SubqSwIgiB4Pkb09LbMFc/JPH+d5bHAJbZPkHRsPv4M8Hpgy7ztRjIG222AssfYvgjSSJPtqwFs36ECE9mB5oBOabq0lYSenlHMe3riiNW3ZMnIWsOPH7+8WLZnedkfyfJCP1fjxpW1tbewvtJ2jhld9v26cEmZOfVqk8pcPi9aUFZfO1i2rMx199ixZabmpfW1ilYNr/WzzvIg0nQEwCmkEaPP5PRT8/zN1ZLW6jMf35f6B31R36qbbWunrQMKgiBYJemt/vqeImlG3fF024PNx29Qp1QeJc2HA2wEPFiXb3ZO608B7SBpLmkKZqKkeSTFI5I1clOEAgqCIGgzdlM9oCdsTyuvy1ZhwCzbLe0mDqqAJL3U9i2trDQIgiBYkd7hDTT9WG1oTdJU4PGc/hCwSV2+jXNaQyT15/IMANvfbqZRVQa8fyTpWkkfkrRmM4UHQRAE1agFpRtsK+Rc4F15/13AH+rSj1Bid2DOAPM/kNyerU5yBvBB0nDdRsBRwM7NNmrQHpDtPSRtCbwXuF7StaS1OBc3W1kQBEHwfIxY7tZYweV1lnuS5opmA18CTgDOyGsu7wcOydnPJ5lg300ywx4w6oHtL+c6rgB2tj0vHx8H/KnZtlaaA7J9l6QvkJyEfg/YKduPf8722c1WGgRBEKxIq5wL97POElL00r55DRxdUM0GpMgENZbynGFDZarMAb2MpBUPAC4G/sP2DZI2JHkzCAUUBEEwBLrQFc+pwLWSzsnHbyStP2qKKowiEVsAAB8ySURBVD2g7wM/J/V2nrX7tv1w7hWtdPT0jmLeouY9Yo8bU7ZmZfQIu51/4pnJxbLjx5Zd44TC9TwLC34HgJ7C8AFLl5cZ+ZTel+U9ZfXNmV+2Tm30qML1SkvLQo0ATBy7rEiu9O+idEHnvMXt9YLv4TVCaCm2j5d0AbBHTnqP7RubLaeKAjoAWGS7B0DSKGCC7YW2f9VshUEQBEEf3NQ6oI4gu+wZyG3PoFT5VPh/pEikNSbltCAIgqAFGNHjUZW2diJpUIVTJU+NKj2gCdkhKJCcg0qaVLWCIAiCYHC6pAf0khyxoD9EclRaiSoKaIGknWseUiXtwvN9ALUESVuTPLbWeCHwReBS4CfAaiRPru+wPbdOblPgNuA42/WB9Grnf0GyWxdwJ8nj6/y++YIgCNpFl8wBbVMhT2VnfFUU0EeB30l6mPQCfwHwtqoVNIPtf5BCciNpNGlF7jmkcAqftH25pPeSPHX/d53ot4ELBij6YzWFJenbwIdJdvFBEARtJ1nBtbsVg2P7/laWV2Uh6nU51sPWOekfOTz2cLMPcI/t+yVtBVyR0y8G/kxWQJLeCNwLLOivoDrlI9J8Vhf81EEQrEp0mRl2S6jqjPTlpEBGY4CdJWH71IFFhsyhpABHALeS3Ib/nhScbhMASauRXIrvC3xyoMIk/ZK04vc24BMNzh8JHAkwdfzqLC8w5exZVubqftzoMvfxpSxeXu6DdkHhNbKwTGxcYZiDUhYXmmGXTg4X+oRkSeFvOFpl93PZEGLVlJrEl/5djPRv0SpWxa/iQX8pSb8Cvgn8G0kRvZw0nzJsSBoHvAH4XU56L/AhSdeT/BDVVuAeB3ynynyO7fcAGwK302AI0fZ029NsT1t73MjFAgqCILBhuVVp61QkfbRZmSqfUdOAbbPLhpHi9cANth+DFG0PeC1AHo47IOfbDXiLpK8DawG9khbb/kGjQm33SPo/4NPAL4f5GoIgCCrTqoB0beTjwHebEaiigGaRDA8G8pDaag7jueE3JK1v+/G8CPYLJIs4bO9Rl+c4YH5f5ZPnfV5k++68/wbgjuG/hCAIgmqYFUONdilNa9AqCmgKcFv2gv1sHGDbb2i2sipImkya0/nPuuTDJNUc5p1Nhd6LpPOB95Oi/50iaQ3SDbqJ5EY8CIKgY+gGK7hBGJaQ3Mc1345ybC8A1u2TdiJw4iByx/U53r/u8FWtal8QBEHrUVesA+oTghueUzo1C+OmqGKGfbmkzYAtbf+/7AWhpWFZgyAIVmW6aB3Q6q0sr0o4hg+QzJPXAV5Ein73ExrEllhZ6OkVc5c07/231Hvvkp4ys9HSL6blhWaxABMLTWMnjCmTm7+0zOx7cqEH5qWF5saLFk8oklu9sJ2lv+GYcgv8YhYWmowvKfQU/q8lZV6t1y78LVpFTxcYIUiaQIp++mLgZuAk22Wu4KnmjPRo0hDWXEjB6YD1SysMgiAIVsTZG3aVrc2cQrKMvoW0rvJbQymsyqfJEttLkwEZSBrDqrlmKgiCYNjohjkg0pKcl8KzPjavHUphVRTQ5ZI+B0yUtC/wIeCPQ6k0CIIgWJEO6N1U4dlxStvLax2TUqoooGOB95G6XP8JnE+KkBoEQRC0gC5aB7SDpFokApE6JvPysW2v0UxhVazgeoGf5S0IgiBoMaY7jBBst9QCuooV3L00mPOx/cJWNiQIgmBVphuG4CS9HHjQ9qP5+AjgzaQ4bV+2/VQz5VX1BVdjAskb9TrNVBIEQRAMTBfoH4CfAq8BkPRqUly1j5DiuE0H3tJMYVWG4J7sk/Td7JX6i81U1E0YirzO9hYO4pZ2vUvlhtLRn1+4pmNR4ZqO0nU5pePppWtPFpe2s/CtUxoeYU3K1rosLrwv7WD1MWXLUuYNIUzJUEkLUVs3BCfpGOADpD/3n9n+bvaX+QHgXznb52yf32TRo+t6OW8Dpts+CzhL0sxm21llCG7nusNRpB5R+36pIAiClZBW9YAkbU9SNLuSQtdcKOm8fPo7tr85hOJHSxqTF5/uQ46hlmlaL1QRqF9otJw01ndIsxUFQRAE/dDaRaYvAa6xvRBA0uXAm1pU9mmkpTlPAIuAK3MdLwbmNFtYlSG4vZotNAiCIKhOsoKrnH2KpBl1x9NtT687ngUcL2ldkpLYH5gBPAl8OBsOzAA+YfvpptppHy/pEmAqcFFdnLhRpLmgpqgyBPfxQRr07WYrHaCuCcAVwHhS2860/SVJe5Oiso4DrgfelxdBrQn8Gtg05/+m7eeFapB0PHAEsLbt1VrV3iAIgtYgeqvPzj5hu9+o1LZvl/Q14CJgATAT6AF+DHyVpO++Shrdem+zLbV9dYO0O5stB6r5gptGip+zUd6OAnYmhcZuqWdUUryhvW3vQLKq2E/SK0n+hw61vT1wP/CunP9o4Lacf0/gWzmcd1/+SBoPDYIg6Ejsalu1svwL27vYfjXwNHCn7cds99St7Wz7O7HKHNDGwM6258GzkUf/ZPvwVjcmd+fm58OxeesBltZp2IuBzwK/IGny1XOk09WAp0jzVH3LvTq3vdVNDoIgGDKt9oRQF0V6U9L8z+6SptquRbY+mDRU11aqKKANSJYUNZbmtGFB0mjSMNuLgR+SnN2NkTTN9gySnfkmOfsPgHOBh0m9sbdl7V5S75Fki471xq3BggKTzFKz6NEa2RUACwvDP0C5WfRaY8tMY0vv6RMF4TQAxhR+oywpDI+wrFBu0uiy19UDC5qOGQbAOuPKQxX8a0lZSI1JY8qucfyoMrmSpRetpMULUc/Kc0DLgKNtPyPp+5J2JOm7+1gx6nRT5I/+dwAvtP2VrOheYLsp56RV3rKnAtdKOicfv5E0JDYs2O4BdpS0FnAOsB1wKPAdSeNJ45q14DKvI41v7k2KVXSxpCttz31+yYPWO520kIotJ0/tkjVhQRCsLLTypWN7jwZp72xhFT8iddr2Br4CzAPOAl7eTCFVrOCOl3QBULug99i+sbm2Nk/W2JcC+2W79T0AJL0W2KrWFuCEPHR3d3YbtA1DdBEeBEEwktjQ0yXeSDO72d5Z0o0Atp/uZ/59QKqOp0wC5to+EZgtaYtmK6qCpPVyzwdJE4F9gTskrZ/TxgOfIUVkBXiAHJlV0gbA1sA/h6NtQRAEw0lvxa1DWJanSwzp3U1B8wZVQJK+RHrpfzYnjSWZPg8HU4FLJd0MXAdcbPs84FOSbieFgP2j7b/k/F8FXinpFuAS4DO2n8jtftYthKSvS5oNTJI0OxtSBEEQdATJFU9XRESt8T3SFMn6eZnLX4H/abaQKnNABwM7ATcA2H5YUqvNr8ll35zr6pv+KeBTDdIfBl7bT1k71u1/Gvh061oaBEHQWjpHtwyO7d9kn6D7kPzNvdH27c2WU0UBLbVtSbWu1uRmKwmCIAgGpoN6N4OSHRScbvuHQymnigI6Q9JPgbUkfYC0cnalDk63zOKRxc2bjo4pNKeeOLpUrmxEuNT0F2BhT5nshNFl5ttPLS3zwlxqUrta4T0dVXhLVfjd+9TSMn/Ape+4BT3jCyVH/sX6VG/ZvVlvfLmp+VDploB0dawOXCTpKeB04He2H2u2kAF/qWzrfTrJsmwuaZL/i7Yvbr69QRAEQX90Uw/I9peBL0t6GSksw+WSZtt+TTPlDKiA8tDb+bZfSvJAEARBELQY011zQHU8DjxKcnS6frPCVcZFbshhWIMgCILhoKIFXKf0kiR9SNJlJOvjdYEP2H5Zs+VUGSzdDThc0n0kz6oidY6ariwIgiBojLurD7QJ8FHbTUdBradfBSRpU9sPkNzdBEEQBMNEbR1Qt2D7s4PnGpyBekC/J3nBvl/SWbbf3IoKgyAIgufTREC6jkDS2sCWwIRamu0rmiljIAVUbxP4wuaa1t2UTgg+tqTMZHjsqLInb8KosvqG4s6j1I/2E4X3ZkGh2feC5YVetClr54sml3n7XlDombzUYLfUBH8oBsrjCp/vhxaXmVOX3psFy8s8qLeKqrF+OgFJ7weOIYXrmQnsDvyd5Jy0MgM9/e5nPwiCIGghtXhAXeQL7hiS5+v7be9F8mDzTLOFDPSJsYOkuaQPiol5H54zQlij2cqCIAiCxnTTHBCw2PZiSUgab/sOSVs3W0i/Csh22VhEEARB0BxNhNvuEGbnyAW/J8Vhe5oU5K4pykNjDgOSTpL0uKRZdWk7SPq7pFsk/VHSGjl9nKRf5vSbJO3ZT5nHSXpI0sy87T9ClxMEQVCJ5Iqn2tYh3Gn7GdvHAf8N/AK4o9lCOkoBAScD+/VJ+zlwbPbGcA7PecX+AEBO3xf4lqT+ruc7tnfM2/mtb3YQBMHQ6LI5oH1rO7Yvt31ufVpVOkoBZRO+p/okbwXUTPsuBmrm4NsCf8lyj5MmwKaNQDODIAhaijF2ta2dSPpgjr+2jaSb67Z7SfHamqKjFFA/3AoclPffSlqBC3AT8AZJY3KE1l3qzvXlw/kmnZRt14MgCDqKLnHF81vgP4A/5P9r2y62D2+2sDJD+5HlvcD3JP03cC6wNKefBLwEmAHcD/wN6Gkg/2NS5FTn/7+Vy1wBSUcCRwKsPWZN5hWsIyldz/PggrKVCxtNKhLjkUVlcgDrFnrlX6P56BZA+XqexY2ehApMnVj2G85bXvYtV7rOaXThYpeFhfdzzbHlb77Zi0b2O7f03sxd1t5wCO3XLYNjew4wBzisFeV1fA/I9h22X2t7F+A04J6cvtz2x/K8zkHAWsCdDeQfs91ju5cUx2jXfuqZbnua7WmTxxS+2YMgCArotpDckk7JVnC147UlndRsOR2vgCStn/8fBXwB+Ek+nlSLzippX2C57dsayE+tOzwYmNU3TxAEQVsx9NiVtipIOkbSLEm3SvpoTltH0sWS7sr/D2U64mW2n114avtp0mLUpugoBSTpNJI7h60lzZb0PuAwSXeSTPweBn6Zs69PChVxO/AZ4J115fxcUs0g4evZVPtmYC/gYyN0OUEQBJVoZQ9I0vYkK+FdgR2AAyW9GDgWuMT2lqQwCscOocmj6hWYpHUomNLpqDkg2/2NK57YIO99pAitjcp5f93+OxvlCYIg6CRaaOD2EuAa2wsBJF0OvIlkzLVnznMKcBnp472EbwF/l/Q7knectwDHN1tIRymgIAiCVZXe6mYIUyTNqDuebnt63fEs4HhJ6wKLgP1Jxlob2H4k53kU2KC0rbZPzW2oOR99U6MpkMEIBRQEQdABNNEDesJ2v2sebd8u6WvARaQgojPpYyFs25KK+1ySJpCmcC4D7ra9uKScUEANWNID981vXm71sWVmnEsLTVseW1Rqvl3e139maVmdpWbYI+165JHCezp7Qdka9amTCs23l5XdmPUnFolRWB0A4wtnmieMLqv06cJndEIbvV/WvGG3rDz7FyT3OEj6H2A28JikqbYfycZZjzdbrqQxwP+QlrLcTxp+20TSL4HP224qckdHGSEEQRCskhh6el1pq0Kd9fCmpPmf35LWUb4rZ3kXaTFps3wDWAfYwvYutncGXkRaBvPNZguLHlAQBEGbST2glnb3z8pzQMuAo20/I+kE4IxsXXw/cEhBuQcCW7nOJ5DtuZI+SLJUPqaZwkIBBUEQdACtdPNme48GaU8C+wy96Oe31HZPyZxSDMEFQRC0GWN6K25t5jZJR/RNlHQ4BeEYogcUBEHQAXRJQLqjgbMlvRe4PqdNAyaSPM00RSigIAiCNmNguTso2k8/2H4I2E3S3sB2OflG278uKS8UUAPGjhLrT2jelPORRWWfMA/MX1Ikt8XqZa6pL3643B32hhMnFMmN7jdW4MA8OL/sj/Jl65SZ4t7wZFl9W6xeZsP7zNLB8zRi7rKydq4xrqydiwu9dkO5efOiwjrHFU4slHrRbhVu//BaZWz/hRyPTdINQCigIAiCbqXz+z/9Uqy6QwEFQRC0mWEwwx5JflYqGAooCIKg7bQ/3HYptn9UKttRZtg5ZPbjkmbVpTWMYaHE9yTdncNt79xPmbvkcAx35/xtHukNgiB4Pl1iho2kXSW9PO9vK+njkvYvKaujFBBwMrBfn7T+Yli8Htgyb0eSQm834sek2Bi1vH3LD4IgaCsGeir+ayeSvgR8D/ixpP8FfgBMBo6V9Plmy+uoITjbV0javE9yfzEsDgJOzatyr5a0Vs3RXk0wO9xbw/bV+fhU4I3ABcN4GUEQBE3SGb2bCrwF2BEYTwrpsHF2xfNN4BqajAnUUQqoH/qLYbER8GBdvtk57ZG6tI1yet88z0PSkaSeFJNGrcmtTy9vuqHrTSyzN11rXJmr6H/MXVgkt+mkQpfIwKLlZX8kd80p+3KrGoK4L48vLnu0n1rSlDPfZ1la6LZ7dOGI8GaFZt+PLiyztRo/BBvlDSaWyT5QaIK/oPAZbacZdhcZISy33QMslHSP7bkAthdJavoH67QhuAHJvZ1h+ZVsT7c9zfa08aMmDUcVQRAE/dJb8V+bWSqp9oLcpZYoaU0KLMm7oQfUXwyLh4BN6vJtnNPqeSinD5QnCIKgzRg334FoBx8jRVnFXsF1w1ieC/VQmW7oAfUXw+Jc4IhsDbc7MKd+/gcgH8+VtHu2fjuCshgYQRAEw0ZtCK4LrODeCVwv6f8kvVvSCwBsP2H7lmYL66gekKTTSAYHUyTNBr4E9BfD4nxSrPO7gYXAe+rKmWl7x3z4IZJ13USS8UEYIARB0GGYHpqfdx5pbH8QQNI2JEvkk/Pw26XAhcBVeY6oEh2lgGwf1s+p58WwyPNBR/dTzo51+zOA7VvSwCAIgmHAQG93DMEBYPsOUviF70iaCOwFvBX4Nsk7diU6SgEFQRCsqnSAgUERthcB50vawPZHmpHthjmgIAiClRx3ixXcQHy5WYHoATWgx2bO8ubXg8ybVzaGu8RlcnM0r0ju/sX3F8kBTBuzdZHcA0vL2rrZ+NWL5P45t+yebr92WbiJu+eWrR9a2FO2PuruOWWT0bf03l0kt92oFxXJAVy94LEiubEuWx9Xak320gkbDJ5pmEjrSzpauQAg6eb+TvHcGs3KhAIKgiBoO6aHso+YEWYD4HXA033SBfyt2cJCAQVBELQZ424xQjgPWM32zL4nJF3WbGGhgIIgCDqA3jY7Gq2C7fcNcO7tzZYXRghBEARtxxWXoVbrJUn6mKRbJc2SdJqkCZJOlnSvpJl523HwkoaX6AEFQRC0GQO91ddvDoikjYD/ArbNTkLPAA7Npz9l+8yWVNQCQgEFQRC0HbfaCm4MMFHSMmAS8HArC28VoYAaIGC8mnd3P6d3cVF9a40uM/2d0Lt2kdxWY9YvkgN4cumSIrkllMk9vqTMFHetMeOL5K566pkiuVJz8Qmjy0bBS6M3H7DaNkVypSEOAJ6Zv1aR3AQVhtTwgiK5Py28uEiuNTRlBTdF0oy64+m2pz9bkv1Qjs/zAMlx6EW2L5L0duB4SV8kB/e0XfaH2SJCAQVBELSZJtcBPWG7X3c3ktYmBezcAngG+J2kw4HPkmKqjQOmkwJ7fmUIzR4yYYQQBEHQdozdU2mrwGuAe23/y/Yy4GzglbYfcWIJ8Etg12G8oEqEAgqCIOgAWuiK5wFgd0mTchiafYDbczw1ctobgVnDdjEV6ZohOEn3AfOAHlJY2GmS1gFOBzYH7gMOsd13hS6SLgR2B/5q+8CRanMQBEE1jFu0Dsj2NZLOBG4AlgM3kobcLpC0HmmaeyZwVEsqHAJdo4Aye9l+ou74WOAS2ydIOjYff6aB3DdIliD/OQJtDIIgaAoDKwYYHWJ59pdI8dTq2btlFbSIbh+COwg4Je+fQupWPg/bl5B6T0EQBJ2HTY+XVdpWJrqpB2TgIkkGfprNDjeoC8P9KAXeWGtIOhI4EmCcVuPR3jlNl/GicesU1X3/0ubrAlhe2GWfueyOIjmAaaN2K5LbRGX35q/L/1Ikt0fhx97sUWXeotdYWhbz8CHKPEVvmiIhN81dzzwxeKYGjKPMHB5gicosfR/U80bTKzFtbJnH9klLX1Ukd0nzPjgb0PJ1QF1BNymgf8v27esDF0ta4S1q21k5FZEV2nSA1cas3/bA60EQrDq0egiuW+gaBWT7ofz/45LOIZkQPiZpqu1HsoXH421tZBAEQRGtM0LoJrpiDkjSZEmr1/aB15JMCM8F3pWzvQv4Q3taGARBMDTs3krbykS39IA2AM5J5uuMAX5r+0JJ1wFnSHofcD9wCICkacBRtt+fj68EtgFWkzQbeJ/tP7fhOoIgCJ6HMT2FkZG7ma5QQLb/CezQIP1J0iKrvukzgPfXHe8xrA0MgiAYIitb76YKXaGAgiAIVmrsqm52VipCATWgh17mj5rbtNwFi64qqm/JsjIz7PFj1yySG13oZRjgD/NOKJLbZ82PF8mNGVXm1fqvyy8rkhurSUVy9+uhIrk1ess8RY8dVTZ9u6DguQaY8cw5RXIAO615eJFcE96hV+DepWUezW/vubxIrlWEGXYQBEHQBhxDcEEQBMHIE+uAgiAIgjZhesMKLgiCIGgH0QMKgiAI2oAhjBCCIAiCEcfRAwqCIAjagAkz7CDTw1Ke8cNNy71izH5F9V3lPxbJla4fWn/SdkVyAMsmv6RIbsrYCUVyWj66SG5rdiqSu3Hp+UVyHlf28rhvSZkr/zkTdy+Su2f+pUVyu6/xviI5gCWUhWN4aMF1RXKrT1q3SG7n0fsWyf2Fm4rkViTMsIMgCIK2YLySBZurQiigIAiCjiB6QEEQBMGIY1gFh+C6JR7QJpIulXSbpFslHZPTj5P0kKSZedu/H/mTJD0uadbItjwIgqAarvhvZaJbekDLgU/YviEHprte0sX53Hdsf3MQ+ZOBHwCnDmMbgyAIhsCq1wPqCgVk+xHgkbw/T9LtwEZNyF8hafPhaV0QBMFQWTWt4GR3V5cuK5IrgO2BjwPvBuYCM0i9pKcHkDvP9vb9nD8SODIfbk8K+b2qMAV4ot2NGCFWpWuFuN6RYDPb6w2lAEkXktpehSdsl6356DC6SgFJWg24HDje9tmSNiA9bAa+Cky1/d5+ZDdnAAXUJ+8M29Na1vAOZ1W63lXpWiGuN+hsusIIAUDSWOAs4De2zwaw/ZjtHqe+68+AXdvZxiAIgqA6XaGAJAn4BXC77W/XpU+ty3Ywq9awWRAEQVfTFQoIeBXwTmDvPibXX5d0i6Sbgb2AjwFI2lDSsz5VJJ0G/B3YWtJsSYP5FZk+PJfRsaxK17sqXSvE9QYdTFfNAQVBEAQrD93SAwqCIAhWMkIBBUEQBG0hFFAfJO0n6R+S7pZ0bLvbM5z05+JoZUbSaEk3Sjqv3W0ZbiStJelMSXdIul3SK9rdpuFC0sfyMzxL0mmSyuJ/BCNKKKA6JI0Gfgi8HtgWOEzStu1t1bBSc3G0LbA7cPRKfr0AxwC3t7sRI8SJwIW2twF2YCW9bkkbAf8FTMvr/EYDh7a3VUEVQgGtyK7A3bb/aXsp8H/AQW1u07Bh+xHbN+T9eaQXVGUXR92GpI2BA4Cft7stw42kNYFXk5YvYHup7Wfa26phZQwwUdIYYBLQfETJYMQJBbQiGwEP1h3PZiV+IdeTPUXsBFzT3pYMK98FPs2q4fVxC+BfwC/zkOPPJU1ud6OGA9sPAd8EHiD5jJxj+6L2tiqoQiigoObi6Czgo7bntrs9w4GkA4HHbV/f7raMEGOAnYEf294JWACslHOaktYmjVRsAWwITJZ0eHtbFVQhFNCKPARsUne8cU5baWnk4mgl5VXAGyTdRxpa3VvSr9vbpGFlNjDbdq1HeyZJIa2MvAa41/a/nOJanw28ss1tCioQCmhFrgO2lLSFpHGkicxz29ymYaM/F0crI7Y/a3tj25uTfte/2F5pv5JtPwo8KGnrnLQPcFsbmzScPADsLmlSfqb3YSU1uFjZ6Ip4QCOF7eWSPgz8mWRJc5LtW9vcrOGk5uLoFkkzc9rnbJ8/gEzQPXwE+E3+mPon8J42t2dYsH2NpDOBG0iWnTcSLnm6gnDFEwRBELSFGIILgiAI2kIooCAIgqAthAIKgiAI2kIooCAIgqAthAIKgiAI2kIooCAIgqAthAIKVlokrVsXwv1RSQ/l/fmSfjQM9Z0s6V5JRw2QZ48c/mJWq+sPgm4j1gEFqwSSjgPm2/7mMNZxMnCe7TMHybd5zrf9cLUlCLqB6AEFqxyS9qwFpJN0nKRTJF0p6X5Jb5L0dUm3SLow+8pD0i6SLpd0vaQ/S5paoZ635gBpN0m6YrivKwi6jVBAQQAvAvYG3gD8GrjU9kuBRcABWQl9H3iL7V2Ak4DjK5T7ReB1tnfIZQdBUEf4ggsCuMD2Mkm3kHwAXpjTbwE2B7YGtgcuTr4uGU2KOzMYVwEnSzqD5KE5CII6QgEFASwBsN0raZmfmxjtJf2NCLjV9iuaKdT2UZJ2I0VhvV7SLrafbGXDg6CbiSG4IBicfwDrSXoFpBhKkrYbTEjSi2xfY/uLpOikmwwmEwSrEtEDCoJBsL1U0luA70lak/R3811gsFAd35C0JakHdQlw0/C2NAi6izDDDoIWEWbYQdAcMQQXBK1jDvDVwRaiAn8EnhixVgVBhxI9oCAIgqAtRA8oCIIgaAuhgIIgCIK2EAooCIIgaAuhgIIgCIK28P8BdWXYa8qo2uUAAAAASUVORK5CYII=\n",
            "text/plain": [
              "<Figure size 432x288 with 2 Axes>"
            ]
          },
          "metadata": {
            "needs_background": "light"
          }
        }
      ]
    }
  ]
}