nyanpasu64/fir_filter.ipynb

## fir_filter.ipynb
{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import control"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "from numbers import Number\n",
    "from functools import singledispatchmethod\n",
    "\n",
    "FIR_LEN = 8\n",
    "FIR_SPAN = FIR_LEN - 1\n",
    "\n",
    "class Polynomial(dict[int, complex]):\n",
    "    def _normalize(self):\n",
    "        for k in list(self.keys()):\n",
    "            if not self[k]:\n",
    "                del self[k]\n",
    "        return self\n",
    "    \n",
    "    @singledispatchmethod\n",
    "    def __mul__(self, _):\n",
    "        return NotImplemented\n",
    "    \n",
    "    @__mul__.register\n",
    "    def _(self, factor: Number):\n",
    "        return Polynomial({k: v * factor for k, v in self.items()})._normalize()\n",
    "    \n",
    "    def to_control(self, size: int, min_tap: int) -> np.ndarray:\n",
    "        out = np.zeros(size)\n",
    "        for k, v in self.items():\n",
    "            out[(size - 1) - (k - min_tap)] = v\n",
    "        return out\n",
    "\n",
    "@Polynomial.__mul__.register\n",
    "def _(self, other: Polynomial):\n",
    "    \"\"\"python was a mistake.\n",
    "    https://stackoverflow.com/a/24064102\"\"\"\n",
    "    out = Polynomial()\n",
    "    for k1, v1 in self.items():\n",
    "        for k2, v2 in other.items():\n",
    "            out[k1 + k2] = out.get(k1 + k2, 0) + v1 * v2\n",
    "    return out._normalize()\n",
    "\n",
    "def make_transfer(num: Polynomial, den: Polynomial):\n",
    "    all_keys = list(num.keys()) + list(den.keys())\n",
    "    min_tap = min(all_keys)\n",
    "    max_tap = max(all_keys)\n",
    "    \n",
    "    size = max_tap - min_tap + 1\n",
    "\n",
    "    out_num = num.to_control(size, min_tap)\n",
    "    out_den = den.to_control(size, min_tap)\n",
    "    return control.tf([[out_num]], [[out_den]], True)\n",
    "    \n",
    "def SIGMA(fir) -> Polynomial:\n",
    "    assert len(fir) == FIR_LEN\n",
    "    poly = Polynomial()\n",
    "    for i in range(FIR_LEN):\n",
    "        amplitude = fir[FIR_SPAN - i] / 128\n",
    "        if amplitude:\n",
    "            poly[-i - 1] = amplitude\n",
    "    return poly\n",
    "\n",
    "def H(fir, feedback=127) -> control.TransferFunction:\n",
    "    num = SIGMA(fir)\n",
    "    den = num * (-feedback / 128)\n",
    "    den[0] = 1\n",
    "    return make_transfer(num, den)\n",
    "\n",
    "def H_den_only(fir, feedback=127) -> control.TransferFunction:\n",
    "    \"\"\"Hide all zeros, so they don't zoom out the plot.\"\"\"\n",
    "    num = Polynomial({0: 1})\n",
    "    den = SIGMA(fir) * (-feedback / 128)\n",
    "    den[0] = 1\n",
    "    return make_transfer(num, den)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "assert Polynomial({-1: 1, 0: 2}) * Polynomial({-1: 1, 0: 2}) == Polynomial({-2: 1, -1: 4, 0: 4})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$\\frac{-2 z - 1}{z^2 + 2 z + 1}$$"
      ],
      "text/plain": [
       "\n",
       "  -2 z - 1\n",
       "-------------\n",
       "z^2 + 2 z + 1"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "H([0, 0, 0, 0, 0, 0, -128, -256], 128)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$$\\frac{z^2}{z^2 - 1}$$"
      ],
      "text/plain": [
       "\n",
       "  z^2\n",
       "-------\n",
       "z^2 - 1"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "H_den_only([0, 0, 0, 0, 0, 0, 128, 0], 128)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# TODO https://stackoverflow.com/questions/39894896/set-points-outside-plot-to-upper-limit\n",
    "\n",
    "def pzmap(sys, ax, fig):\n",
    "    from numpy import real, imag\n",
    "    if not isinstance(sys, control.lti.LTI):\n",
    "        raise TypeError('Argument ``sys``: must be a linear system.')\n",
    "\n",
    "    poles = sys.pole()\n",
    "    zeros = sys.zero()\n",
    "\n",
    "    # Plot the locations of the poles and zeros\n",
    "    if len(poles) > 0:\n",
    "        ax.scatter(real(poles), imag(poles), s=50, marker='x', facecolors='k')\n",
    "    if len(zeros) > 0:\n",
    "        ax.scatter(real(zeros), imag(zeros), s=50, marker='o',\n",
    "                    facecolors='none', edgecolors='k')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Direct FIR control\n",
    "\n",
    "Drag the sliders to control the FIR filter coefficients directly. The poles and zeros of the echo section (not the master audio) are shown.\n",
    "\n",
    "All poles must be within the unit circle for the filter to be stable. Poles on the unit circle are marginally stable. \"Poles around the circle\" is generally fine on SNES for low input amplitudes, but may lead to unstable low-level oscillation after input is stopped. Note that coefficients of 1.0 are unachievable (127 is only 127/128 of full-scale)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "0c6b1ad0c9ca4c56a8c00dfa9ee659f3",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(IntSlider(value=0, continuous_update=False, description='fir0', max=127, min=-128), IntS…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from ipywidgets import interact, interactive, fixed, interact_manual\n",
    "import ipywidgets as widgets\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "def move_axis_to_origin(ax):\n",
    "    # set the x-spine (see below for more info on `set_position`)\n",
    "    ax.spines['left'].set_position('zero')\n",
    "\n",
    "    # turn off the right spine/ticks\n",
    "    ax.spines['right'].set_color('none')\n",
    "    ax.yaxis.tick_left()\n",
    "\n",
    "    # set the y-spine\n",
    "    ax.spines['bottom'].set_position('zero')\n",
    "\n",
    "    # turn off the top spine/ticks\n",
    "    ax.spines['top'].set_color('none')\n",
    "    ax.xaxis.tick_bottom()\n",
    "\n",
    "def plot_fir(fir, feedback):\n",
    "    sys = H(fir, feedback)\n",
    "\n",
    "    fig = plt.gcf()\n",
    "    ax = plt.gca()\n",
    "    \n",
    "    ax.axis(\"scaled\")\n",
    "    ax.set_xlim(-1.5, 1.5)\n",
    "    ax.set_ylim(-1.5, 1.5)\n",
    "    move_axis_to_origin(ax)\n",
    "    \n",
    "    unit_circle = plt.Circle((0, 0), 1, fill=False)\n",
    "    ax.add_artist(unit_circle)\n",
    "    \n",
    "    pzmap(sys, ax, fig)\n",
    "    \n",
    "    plt.show()\n",
    "\n",
    "def fir_callback(**kwargs):\n",
    "    # don't put feedback as a non-kwarg, it reorders it\n",
    "    plot_fir([kwargs[f\"fir{i}\"] for i in range(FIR_LEN)], kwargs[\"FEEDBACK\"])\n",
    "\n",
    "def fir_slider(default=0):\n",
    "    return widgets.IntSlider(min=-128, max=127, step=1, value=default, continuous_update=False)\n",
    "\n",
    "interact(fir_callback, **{f\"fir{i}\": fir_slider() for i in range(FIR_LEN)}, FEEDBACK=fir_slider(127))\n",
    "None"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Parametric four-pole FIR filter\n",
    "\n",
    "The two-pole resonant filter requires FIR coefficient values up to -256 to 256, which are unachievable on SNES hardware. However, adding two poles on the left-hand side brings all coefficient values in-bounds within \\[-128, 127)!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "c07a807b555c4bca81d5cf0833c6c888",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(FloatSlider(value=1.5, continuous_update=False, description='coeff', max=2.0, min=0.3, s…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "def clamp(n, smallest, largest):\n",
    "    return max(smallest, min(n, largest))\n",
    "\n",
    "def four_pole_fir(coeff, resonance):\n",
    "    fixed_poles = Polynomial({0: 1, -1: 1, -2: 0.3})\n",
    "    sweep_poles = Polynomial({0: -1, -1: coeff, -2: -resonance})\n",
    "    all = fixed_poles * sweep_poles\n",
    "    \n",
    "    assert all.get(0, 0) == -1\n",
    "    def get_coeff(i):\n",
    "        x = all.get(i, 0)\n",
    "        x = int(x * 128)\n",
    "        return clamp(x, -128, 127)  \n",
    "    \n",
    "    fir = [0, 0, 0, 0] + [get_coeff(i) for i in range(-4, 0)]\n",
    "    return fir\n",
    "\n",
    "def four_pole_callback(*, coeff, resonance):\n",
    "    fir = four_pole_fir(coeff, resonance)\n",
    "    plot_fir(fir, 127)\n",
    "    print(fir)\n",
    "\n",
    "interact(\n",
    "    four_pole_callback,\n",
    "    coeff=widgets.FloatSlider(min=0.3, max=2, step=0.01, value=1.5, continuous_update=False),\n",
    "    resonance=widgets.FloatSlider(min=0.5, max=1, step=0.01, value=1, continuous_update=False),\n",
    ")\n",
    "None"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "At a given resonance level, when sweeping the frequency, what shape do the curves take on? Can we generate them through linear DAW sweeps? (yes.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXkAAAD4CAYAAAAJmJb0AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAAAycklEQVR4nO3deXxU9b3/8dd3si8TIAFCSEiClSUssgUU2ddJtNbWrW5QXEq9Xq3a1h+4tL3q73d/9nd729tWQOm9va2Pa9VW605mIFGsCtlQUAuyhYCBhJCQZIYAySzf3x9ngAABAplklvN5Ph59ZDLn5JxPTvHNl+/5zPcorTVCCCEikyXYBQghhOg5EvJCCBHBJOSFECKCScgLIUQEk5AXQogIFh3sAjrq37+/zs3NDXYZQggRVjZt2tSgtR7Q2baQCvnc3FwqKyuDXYYQQoQVpdTec22T6RohhIhgEvJCCBHBJOSFECKChdScfGfcbjc1NTUcP3482KUERHx8PFlZWcTExAS7FCGECYR8yNfU1GC1WsnNzUUpFexyukVrTWNjIzU1NQwdOjTY5QghTCDkp2uOHz9OWlpa2Ac8gFKKtLS0iPlXiRAi9IV8yAMREfAnRNLvIoQIfSE/XSOEEJHM09CAc+1aolPTSCmwBfz4YTGSDwW//e1vycvLo1+/fjz77LPn3beiooLo6Ghee+21XqpOCBFOPIcP0/TKq+z93hJ2zpzFwaefwbVuXY+cS0byXbRy5UqKi4vJysrqdLvH4yE6Ohqv18uyZctYuHBhL1cohAhlnqYmXMXFuIqKaC0rB6+X2Nxc+t/3A1IKC4kbNqxHzish3wX33XcfVVVVFBYWcvfdd7N7926ee+45lixZQnx8PJ999hnTpk3jV7/6Fb/73e+48cYbqaioCHbZQogg8zY34yopwVlkp3XjRvB6icnJJu379xrBPnx4j9+nC6uQf+qdf7D1gDOgxxw1OIWfXzf6vPs8//zz2O12PvjgA959993TttXU1LBhwwaioqLYv38/b7zxBh988IGEvBAm5Wtrw7mmCGfRGlo3bASPh5ghQ0i7+26sBTbiR43q1QaMsAr5UHTzzTcTFRUFwMMPP8wvfvELLBa51SGE2fja22l+7TUaX1iN5+BBYgYPJvV7i0kpvIb40b0b7B2FVchfaMQdDElJSSdfV1ZWcuuttwLQ0NDAmjVriI6O5tvf/naQqhNC9DTd3k7z3/5Gwwur8dTWkjBpEoP/77+SOHVqSLRMh1XIh7o9e/acfL1kyRK++c1vSsALEaF8bW20vPUWDc8/j+dALQnjx5Pxv58h6eqrQyLcT5CQF0KILvK1t9P6ySe47HZcJe/jO3KE+HFXkPHU0yRNnxZS4X6C0loHu4aT8vPz9ZkPDdm2bRt5eXlBqqhnROLvJESk0u3ttJaW4lxThKukBJ/LhSUlBev8+fT55rUhMS2jlNqktc7vbJuM5IUQ4gza7aa1rBxn0RpcxSX4WlqwWK1Y580jpbCApKlTUbGxwS6zSyTkhRAC0B4PR8vLcRbZca1bh7e5GUtSEtb587DaCkiaPg1LmAR7RxcV8kqpPwDfBOq11mP876UCrwK5QDVwi9a6SRn/fvkNcA1wFFiitf40cKULIUT3nWx9XP17PHV1WBITSZ471xixT5+OJS4u2CV2y8WO5P8IPAe82OG95UCJ1vpZpdRy//fLgEJgmP9/VwKr/F+FECLozmp9nDiR9MceI3nWTCzx8cEuL2AuKuS11n9XSuWe8fb1wGz/6z8B6zFC/nrgRW3c2S1VSvVVSmVorWu7VbEQQnSDdrtpfuONkG99DJRAzMmndwjuOiDd/zoT+LrDfjX+904LeaXUUmApQHZ2dgDKEUKI02mtafvqK5xrimh57108B2qJv+IKMp56iqTp0yMy3E8I6Ofv/aP2i+rJ1Fqv1lrna63zBwwYEMhyAqorSw2vX7+ePn36MH78eMaPH8/TTz/dy1UKIU7QWnN8+w7q/+M/qCooZM93bqDxD38gLncoQ154ntxXXyF5xoyIDngIzEj+4IlpGKVUBlDvf38/MKTDfln+98JSV5YaBpgxY8ZZi5gJIXpP286dOIvsOO122quqwGIh6aorSb3nbqwLFhDdr1+wS+xVgQj5t4HvAc/6v77V4f0HlFKvYNxwbQnX+fiuLjX8rW99K9ilCmFKbVVVOIuKcNnttO3cBRYLiZMnk7p4kRHsaWnBLjFoLraF8mWMm6z9lVI1wM8xwv0vSql7gL3ALf7d12C0T+7CaKG8q9vVFi2Hui+6fZjTDBoLhed/0lNXlxpev349GzduZNy4cQwePJhf/vKXjB4deouqCREJTrQ+Nr/yKm07doBSJObnk/7TJ0lZuJDoEJ7+7U0X211z2zk2zetkXw3886UUFU46LjU8ceJE9u7dS3JyMmvWrOHb3/42O3fuDHKFQkSWM1sf46+4gvQnnsC6cCEx6QODXV7ICa9PvF5gxB0MHZcaTklJOfn6mmuu4f7776ehoYH+/fsHozQhIorZWh8DJbxCPsTV1dWRnp6OUory8nJ8Ph9pJp4LFCIQPIcO4bQ7OPynP+GuqQn5VR9DjYR8AL322musWrWK6OhoEhISeOWVV+QPoRCXwNPYiGvdOpxrijhaUQFaEz92LIN++iRJM2fKf1cXQZYaDoJI/J2E6C5PU5MR7EVFHC0rB5+P2KFDSSksJKWwgLhhw4JdYsiSpYaFECHJ29yMq7gYZ5Gd1tJS8HqJyckmben3SSm8hrjhw2TU3k0S8kKIXuVtacFV8j5OexGtGzaCx0PMkCGk3XOPMWIfOdJ0wV5RV0G0JZoJAycE/NgS8kKIHud1uXCVlOAqsnNkwwZwu4nJzCTtriVYCwqIHzXKdMEOsOngJlZuXkl5XTkzMmewcv7KgJ9DQl4I0SN8bW241q41pmI++gjtdhM9OIPURYtIKbARP3asKYMdYHP9ZlZsXkFpbSlp8Wksm7yMm4bf1CPnkpAXQgSUr72d5r/+lcYXVuOpryd60CD63X47KYUFxI8bZ9pg11qz5dAWnt/yPJ8c+ITU+FR+kv8TbhlxCwnRCT12Xgl5IURAnPwk6vMv4KmrI2HSJAY/+39JvOoqlCWgC96GDa0125u246h24Kh28LXra/rF9eNHk37Ed0d8l8SYxB6vQUK+i37729+yatUq6urqWLZsGcuXL+90v/Xr1/Pwww/jdrvp378/H374YS9XKkTv0u3tNL/xJg0v+D+JOmECg//1/5A4daopR+1aa3Y07cBR7WDt3rXsde4lSkVxVcZV3Dv2XgpyC3ol3E+QkO+iriw1fOTIEe6//37sdjvZ2dnU19d3uq8Q4U673bSWlhoPvS4uxud0mv6TqLuaduHYa4zY97TswaIsTBk0hSWjlzAvex794oOzxLGEfBd0danhyy+/nBtuuOHkE64GDpTFkkTk0B4PrWVluOx2XGvX4W1pwZKcjHXeXFKu+xZJ08y3hkxVc9XJqZjdLbuxKAv56fncMfIO5ufMJy0h+MuahFXI/6L8F3x1+KuAHnNk6kiWTVl23n26utTwiWma2bNn43K5eOihh1i8eHFA6xWiN2mPh6OVlTjXFOFatw5vUxOWxESS580jpbCApGnTsMTFBbvMXlXdUo292o6j2sGu5l0oFJPSJ/HEyCeYnzOf/gmhtSBhWIV8KOq41LDH42HTpk2UlJRw7Ngxpk6dylVXXcXw4cODXKUQXae9Xo5u2mQ8hGPtOryNjajERKxz5mAtsJE8YwaW+Phgl9mr9jn3nRyxb2/aDsDEgRNZPmU5C3MWMiAxdNeuD6uQv9CIOxg6LjWclZVFWloaSUlJJCUlMXPmTLZs2SIhL8LCydbH3/8nnro6VEICybNnkVJQSPLMGVgSeq7NL1SV15azcstKNh3cBMD4AeNZNnkZC3IWkJ6UHuTquiasQj7UXX/99TzwwAN4PB7a29spKyvjkUceCXZZQpxXZ62P6cv+F8mzZmFJ7L0ukFBSUVfBys0rqTxYycCEgfxo0o8oyC0gIzkj2KVdtG6HvFJqBPBqh7cuA34G9AW+Dxzyv/+41npNd88XyvLy8igoKOCKK67AYrFw7733MmbMmGCXJUSnpPXxbJ8e/JSVm1dSVldG/4T+LJ+ynJuG30RcVPjedwjoUsNKqShgP8aDu+8Cjmitf9nVn5elhoXoWVpr2nbuxGW30/LmW7gPHCBh3Dj6P/igKbtjAA4fP0zJvhLeq3qPTQc3kRafxj1j7+Hm4TcTHx0e9x56c6nhecBurfVeM/5hESJUte3ejbPIjrOoiPbdu8FiIXHKFAb9y89JmjHDdOHe0tZCyb4S7HvslNeV49VeclJyemWZgd4W6JC/FXi5w/cPKKUWA5XAj7XWTWf+gFJqKbAUONlfLoTovraqPbgcdpxrimjbuROUInHyZFLvvAPrggVEm+zZwy1tLXzw9Qc4qh2UHijFoz0MsQ7h7jF3Y8u1Mbzf8Ij8yy5gIa+UigW+BTzmf2sV8Ayg/V//Hbj7zJ/TWq8GVoMxXROoeoQwo/a9e40Ru91O21dfgVIkTJxI+pNPYl24gBiTfUDP1e46GewbDmzA4/OQmZzJ4tGLseXayEvNi8hg7yiQI/lC4FOt9UGAE18BlFK/B9491w8KIS6dr72dltdfp+mvf6Vt6zYAEiZMIP3xx7DabMSkh0erXyBV1FXw4tYX+WT/J7h9bjKSMrgz705suTZGp42O+GDvKJAhfxsdpmqUUhla61r/t98BvgzguYQwvZOtjy+sxlNbS/zo0QxctoyUAhsxGeHX6hcIlXWVrNyykoq6CtLi07h15K0U5BYwtr95164PSMgrpZKABcAPOrz9/5RS4zGma6rP2CaEuETa7ab5jTdoeN7f+jh+PBn/+xmSrjZndwzAZ/WfsWLzCspqjdbHEw/hCJfumJ4UkJDXWrcCaWe8tygQxw4VXVlq+N/+7d946aWXAGOJg23btnHo0CFSU1N7u1wRgTyHD+O02zn8X3/AvX8/8VdcQcZTT5E0fbopw93tdbOxdiMvbXuJDQc2kBqfyqP5j3LziJsjqjumu+QTr13UlaWGH330UR599FEA3nnnHX79619LwItu8TQ14SouxlVkp7WsDLxe4seMYdDPfkrSzJmmC3e3z015bTn2ajsl+0pwtbvoF9ePH0/6MbeMuKVX12kPFxLyXdDVpYZ/9atfnfyZl19+mdtuuy2IVYtw5W1pwVVcgtNup3XjRvB4iMnOJu3ee0kpLCBuxAhThbvH56G8rpy11Wsp3ldMS1sLyTHJzM2eiy3XxtSMqcRExQS7zJAVViFf96//Stu2wC41HJc3kkGPP37efbq61PAJR48exW6389xzzwW0VhG5vE4nrpL3cdqLaN2wEdxuYrKySLtrCdaCAuJHjTJdsG86uMkYse8toamticToROZkz8GWY2Na5jRio2KDXWZYCKuQD0Udlxo+4Z133mHatGkyVSPOy9fejsvhwLmmiNaPP0a73UQPziB10SLjoddjxpgq2AE212/m3ap3Wbd3HYePHyYhOoHZWbOx5RrBLjdSL15YhfyFRtzB0HGp4RNeeeUVmaoR52S0Pr5Bwwsv4KmtJXrQIPrdcYcR7FdcYbpgh9NbHxOiE5iZNRNbro3pmdPlJmo3hVXIh4OWlhY+/PBD/ud//ifYpYgQc6L1sfH5F4yFwcaPJ+OZZ0i6eirKYgl2eUHRWevjDcNukBuoASQhH2BvvPEGCxcu7HSEL8xJt7fT8s47NKx6HndNDfFXXMGgp54y7QOvtdZsObSFlZtXsrF2o7Q+9rCALjXcXbLUsIgU2u2mtbTUeDZqSQk+p5P4MWMY8OADpmx91Fqzq3kX9mo7a6vXUu2sJjU+lbtG3yWtjwHQm0sNC2Fa2uOhtbQMp72II+uK8ba0YLFasc6bR8o3v2nK9dp3N+8++WzUqpYqLMrC5PTJLB69mGuHXivh3gsk5IXoBu3xcLSiAmeRHdfatXibm7EkJZE8by4pBYUkTZ+GJdZcrX57WvacDPZdzbtQKPIH5XP7yNuZlzOP/gnmWuI42MIi5LXWETMCCqXpMXFptNfL0cpNOIvW4Fq7Du/hw6jERKxz5pBSWEDSjBlY4sL3cXGXYp9zH45qB/ZqOzuadqBQTBg4geVTlrMwZyEDEgcEu0TTCvmQj4+Pp7GxkbS0tLAPeq01jY2NxMdLr284Ornq4+rVeA7UohISsM6ZjbWggOSZM7GY8P/XyrpKVm1ZRXldOQDjB4xn+ZTlzM+eT3qS+ZY4DkUhH/JZWVnU1NRw6NChC+8cBuLj48+5/o0ITZ2t+pj+6KMkz5qFJdGcc8pntj4+MukRrhl6DYOSBgW7NHGGkA/5mJgYhg4dGuwyhAlpt5uWt94yWh/37yd+3BVkPPW0aVsfwfhEqrQ+hpeQD3khepPWmrYdO3Hai3C+/Y4R7iZe9RHg8PHDFO8tpmhPEZUHK2XVxzAjIS8E0LZz58lno7ZXVYHFQuKUKaQ/8QTJc2abLtybjzdTsq8Ee7WdiroKvNpLbkouj0x6hFtH3CrhHkYk5IVptVVV4SwqwmW307ZzFyhF4uTJpC66E+vChUSnpV34IBGkpa2F9/e9j6PaQWltKV7tJduazd1j7saWa2N4v+Gm+8suEgQs5JVS1YAL8AIerXW+UioVeBXIxXgE4C1a66ZAnVOIi9VeXY3TbsdZZKdt+3ZQioRJE0n/6ZOkLFxI9ABztfo52518sO8DHNUONtZuxOPzkJmcyZLRS7Dl2hiZOlKCPcwFeiQ/R2vd0OH75UCJ1vpZpdRy//fLAnxOIc7rxKqPTa++Stu2bQAkTJxI+uOPY7XZiEkfGOQKe19lXSV/2vonPtn/CW6fm8zkTBaNWoQt18aoVHOtXR/penq65npgtv/1n4D1SMiLXnJm62P86NGkP7bcCPZB5mz169j6mBafxu0jb8eWa2NMf/OtXW8WgQx5DaxVSmngBa31aiBda13r314HnPXpCKXUUmApQHZ2dgDLEWal3W5a3n6bhpWr5IHXftL6aF6BDPnpWuv9SqmBwDql1GnP6dNaa/9fAJzx/mpgNRirUAawHmEynsZGnA4Hh//7j7i//tr0rY/t3nY2HtjIy9tf5pP9n5AanyqtjyYUsJDXWu/3f61XSr0BTAEOKqUytNa1SqkMoD5Q5xMCwNPUhGvdOpxFRRwtKwefj/hRo0hfudKUrY9ur5vS2lLs1XY+2PcBLreLfnH9eHjiw9w28jYJdxMKSMgrpZIAi9ba5X+9EHgaeBv4HvCs/+tbgTifMDdvczOu4mKcRXZaS0vB6yUmJ5u0pd8npbCQuOHmavVz+9yU15bjqHZQsq8EZ7sTa4zVeOh1ro2pGVOJiYoJdpkiSAI1kk8H3vD/hxUN/FlrbVdKVQB/UUrdA+wFbgnQ+YTJeFtacJW8j9NeROuGjeDxEDNkCGn33ENKYQFxI83V6ufxeaioqzgZ7M1tzSTFJDF3yFwj2AdPJTbKXEsci84FJOS11lXAuE7ebwTmBeIcwnx0eztOhwPne2s48skn4HYTk5lJ2pLvYS0oJH60+Vr9Ntdv5p3d71C8r5jDxw+TGJ3I7CGzseXamJY5jbgocy1xLC5MPvEqQo52u2l+800aVz2P+8ABogdnkLpoESmFBcSPMWer3+b6zazYvILS2lISohOYlTWLgtwCpmVOIz7afEsci66TkBch42Tr42kPvP4XU7c+nnjg9YYDG0iNT+Un+T/h5uE3yw1U0WUS8iLodHs7Le+toWHVKtz79hmtjz990rStj1prPm/4nOe3PM/H+z+mX1w/fjTpR3x3xHcl3MVFk5AXQaHb22ktLTWejVpcjM/lIm5UHlkmbX3UWrOzeSeOagdrq9dS7aymT1wfaX0U3SYhL3qNdrtpLS3DaS/CVVyCr6UFi9WKdd48Uq691pQP49jdvBt7tR1HtYM9LXuwKAuTB01m8ejFXDP0GpJikoJdoghzEvKiR2mPh6Pl5caIfd06vM3NWJKSsM6fh7WggKRp07DEmqvVr6ql6uSIfVfzLhSK/EH53Jl3J/Oy55GWYK4ljkXPkpAXAae9Xo5WVBjBvnYt3qYmLImJJM+dS0qBjaQZM7DEmavVr7qlGke1A8deBzubdqJQTEyfyGNTHmNh7kL6J/QPdokiQknIi4A52fr4/Au49+9HJSRgnTMba0EByTNnYok3X6tfx9ZHgAkDJ7B8ynIW5CxgYKL5ljgWvU9CXnTbWa2PY8cy8Cc/Jnn2bCwJ5lzl8MzWx4cnPsy1l13LoCRzLnEsgkdCXlwy7fHQ8s67p1ofR48m/cknSJ41y3Q3UE/4suFLVmxeIa2PImRIyIuLorWmbccOnEVFON99D3dNjb/1cQXJc+aYMtwbjzWe9tDrPnF9eGjiQ9w+8nYJdxF0EvKiS9p27sRZZMdpt9NeVQUWC4lXTiF9+TKS580zXbg3HW86Ldh92kduSi4PTXyI20beJq2PImRIyItzatu92x/sRbTv2m0E++TJpC5ejHXBfKLTzNXq19LWQsm+EhzVDspqy/BqL9nWbO4Zcw8FQwsY1neY6f6yE6FPQl6cpm3PHlx2O84iO207doBSJE6aRL+f/ZSUBQuIHjAg2CX2qpa2Fj74+gPs1XbKDpTh0R6ykrNYMnoJBUMLGNFvhAS7CGkS8uJk62PTS3+m7SvjqY0JkyaR/sQTWBcuJCbdfK1+m+s3819f/BcfH/gYj89DZnImi0YvwpZrY1Sq+ZY4FuFLQt7Ezmp9HDWK9McfM4J9kDlb/c5sfbxj5B3Ycm2M6W/OJY5F+JOQNyFpfTybtD6KSNXtkFdKDQFexHgEoAZWa61/o5T6F+D7wCH/ro9rrdd093zi0nkaGnCuXUvTn16kfe9eU6/6CNDubWfjgY38Zcdf+HvN36X1UUSkQIzkPcCPtdafKqWswCal1Dr/tl9rrX8ZgHOIS+Q5fBjXumKcRUUcLS8Hn4+4vDyynvudKVsf3V43G2s34qh28MG+D3C5XfSJ68MPJ/yQ2/Nul9ZHEXG6HfJa61qg1v/apZTaBmR297ji0nmamnAVF+MqstNaVgZeL7E5OaT9YCkpBYXEDTdXq5/b56a8thx7tZ2SfSW42l1YY6zMzTYeen1VxlXERMUEu0whekRA5+SVUrnABKAMmAY8oJRaDFRijPabOvmZpcBSgOzs7ECWYyrelhZcxSU47XZaN24Ej4eY7GzS7rmHlGsKiRthrlY/j89DRV0FjmoHxfuKaWlrISkmiblDjGCfOngqsVHmWuJYmJPSWgfmQEolAx8C/0dr/TelVDrQgDFP/wyQobW++3zHyM/P15WVlQGpxwy0243Tbsf57nsc2bAB3G5isrJIKSzAWlBA/Cjztfp9fuhz3tz1JsV7i2lqayIxOpE52XOw5di4OvNq4qLMtcSxMAel1CatdX5n2wIykldKxQCvAy9prf8GoLU+2GH774F3A3EucXbrY/TgDFIXLSKlsID4MeZs9evY+pgQncDsrNnYcm1My5xGfLT5ljgW4oRAdNco4L+AbVrrX3V4P8M/Xw/wHeDL7p7L7KT18WzS+ijE+QViJD8NWAR8oZTa7H/vceA2pdR4jOmaauAHATiXKfna23EVFdGwcpW0PmKshPlFwxes/nw1H9Z8KK2PQpxHILprPgY6Sxrpie8G3d7OkQ0bcBXZcZWU4DtyhLiRI03b+qi1ZkfTDuMRetUO9rn2kRKbIq2PQlyAfOI1hGi3m9bSUuPZqMXF+JxOLFYr1gULSLmmkKRp01AWS7DL7DVaa3Y17zoZ7NXOaizKwpRBU7hrzF3Ycm1YY63BLlOIkCYhH2Ta46G1rAxnURFH1hXjbWnBkpyMdd48rIUFJF99NSrWXK1+Vc1VOKod2KvtVLVUYVEWJqdPZtGoRczPmU9qfGqwSxQibEjIB4H2ejlaUYFzTRGudevwNjVhSUoied5cUgoKSJo+HYvJgr26pRp7tR1HtYNdzbtQKCalT+K2kbcxP2c+/RP6B7tEIcKShHwvOrP1USUmYp0zh5TCApJmzMASZ74e7o6tjwrFhIETeGzKYyzIWcCARHOtXS9ET5CQ7wXa46Hl3XdpWHmq9XHgT35M8uzZWOLN2cN9ZuvjQxMf4rrLriM9KT3YpQkRUSTke5D2enGuWUPDcyuk9dFva+NWVm5eKa2PQvQSCfkA01rTtmMHzqIinO+twf3116ZufQRoPNZI8d5i7NV2Kg9WSuujEL1IQj5Aju/YcfLZqO179kBUFElXXsnAR3+Cdf58U7U+AjQdb6J4XzGOagcVdRX4tI+hfYby0MSHuHXErSTHJge7RCFMQUK+G9p278ZZZMdpL6J9126wWEicMoXU730P68IFRKeaq9Wv+Xgz73/9Po5qB2W1ZXi1l5yUHO4dey+2XBvD+ppriWMhQoGE/EVqq9qD016Eq8hO286doBSJkybR72c/JWXhQqL7m6vVr6Wthff3vY9jr4OyA2V4tIch1iEsGb2EgqEFjOhnriWOhQg1EvJdoD0eWt56m8Mvvkjb9u2gFAmTJpL+5JNYFy4gZuDAYJfY6z4/9Dm//+L3fLz/Yzw+D5nJmSwevRhbro281DwJdiFChIT8eZzZ+hiXl0f6449htdmISTdnq98/Gv7Bis0r+Gj/R/SN68sdI++gYGgBo9NGS7ALEYIk5DuhvV6c771Hw4qVRutjXh5ZK1eQPGeOaYOss9bH20beJt0xQoQ4CfkOPA0NONeupel/XqK9qoq4ESPI/N1vje4YE4Z7u7edDQc28PrO11n/9XpSYlN4cMKD3D7ydumOESJMmD7kPYcP41q7DqfdztHycvD5iBs5kszf/AbrAvO1Prq9bjbWbsRR7eD9fe9zxH2EPnF9uH/8/dyZd6es+ihEmDFlyHuamnAVF+MqKqK1rBy8XmJzc+l/3w+wFhQQN8xcrX5un5uy2jIc1Q5K9pXgandhjbUyP2c+tlwbV2ZcSYwlJthlCiEugWlC3tvcjKukBGeRndaNG8HrJSY7m7R77yWlsIC4EeZq9fP4PJTXluPYawR7S1sLyTHJzM2eiy3XxtSMqcRESbALEe56POSVUgXAb4Ao4D+11s/29DlP8DqduErex1m0htYNG8HjIWbIENLuvtsI9jxztfp5fB4qD1biqHZQvLeY5rZmkmKSmDNkDrZcG1cPvprYKHMtcSxEpOvRkFdKRQErgAVADVChlHpba721p87pPXKEI/4R+5FPPgG3m5jBg0n93mJSCq8hfvQoUwW71+fl0/pPse+xU7yvmMPHD5MQncDsIbOx5dqYnjmduCjzLXEshFn09Eh+CrBLa10FoJR6BbgeCGjI+1pbcb3/AU67ndaPPkK3txOdkUHqnXeSUlhA/Nixpgp2n/bx6cFPcVQ7WLd3HY3HG0mITmBm1kwKcguYnjmd+GhzLnEshNn0dMhnAl93+L4GuLLjDkqppcBSgOzs7Es6yfHtOzjw6KNEp6fT77ZbsRYUkDBunKk6Y3zax5ZDW3BUO1hbvZZDxw4RHxXPjKwZ2HJtzMicIcv5CmFCQb/xqrVeDawGyM/P15dyjITx48j5859JGG++YP/80OdGsO9dS/3RemItsSeDfVbWLAl2IUyup0N+PzCkw/dZ/vcCSlksJE6cEOjDhiStNV82fImj2oFjr4O61jpiLDFMy5zGjyb9iFlZs+SDSkKIk3o65CuAYUqpoRjhfitwew+fM+Jordl6eOvJqZj9R/YTbYlm2uBpPDjhQeYMmSMfUhJCdKpHQ15r7VFKPQA4MFoo/6C1/kdPnjPSVNRVsGLzCjYd3ES0iuaqwVfxT+P+iTnZc0iJTQl2eUKIENfjc/Ja6zXAmp4+T6SprKtk5ZaVVNRVMCBhAMsmL+O6b1xHn7g+wS5NCBFGgn7jVZzus/rPWLF5BWW1ZfRP6M+yycu4afhN0vIohLgkEvIhoMZVw9q9a3FUO9jauJXU+FQezX+Um0fcTEJ0QrDLE0KEMQn5IKk9Unsy2L9o+AKAsf3HsmzyMm4YdoO0PgohAkJCvhfVtdaxbu86HNUOthzaAsCotFE8MukRFuYsJMuaFeQKhRCRRkK+h9UfrT8Z7J/VfwZAXmoeD018CFuOjSEpQy5wBCGEuHQS8j2g4VjDyWD/9OCnaDTD+w3nwQkPYsu1kZOSE+wShRAmISEfII3HGinZV4K92k5lXSUazeV9L+efxv0TtqE2LutzWbBLFEKYkIR8N2itKa0t5b+//G/K6srwaR9D+wzlvnH3sTBnIZf3uzzYJQohTE5C/hKV15azYvMKPq3/lIGJA7l37L3Ycm0M62uuRwcKIUKbhPxFqqirYOXmlVQerGRgwkAev/Jxbhx2ozxRSQgRkiTku+Co+ygf7f+Iv27/K2V1xidRl09Zzk3Db5KnKgkhQpqE/Dkc8xzj4/0f46h28Peav3PMc4wBCQN4NP9RbhlxiywzIIQICxLyHbR5204G+/qv13PMc4zU+FSuu+w6bLk2JqVPIsoSFewyhRCiy0wf8u3edjYc2IC92s76r9fT6m6lb1xfrr3sWmy5NvLT84m2mP4yCSHClCnTS2tNRV0Fb+1+i/f3vc8R9xFSYlOw5dqw5diYnDGZGEtMsMsUQohuM13Id2x9tMZamZ8zH1uujSszrpRgF0JEHNOE/Jmtj49NeYwbh98o3TFCiIjWrZBXSv0bcB3QDuwG7tJaNyulcoFtwHb/rqVa6/u6c65L4dM+Nh3cxAtbXpDWRyGEKXV3JL8OeMz/LNdfAI8By/zbdmutx3fz+BdNa80XDV8YD73eu5a61rqTD+GQ1kchhNl0K+S11ms7fFsK3NS9ci65DrY2bsVR7cBR7eBA6wGiLdFMHzydH074IfNz5ssTloQQphTIOfm7gVc7fD9UKfUZ4ASe1Fp/FMBznaa8rpx7195LtIpm6uCp3D/+fuZkzyElNqWnTimEEGHhgiGvlCoGBnWy6Qmt9Vv+fZ4APMBL/m21QLbWulEpNQl4Uyk1Wmvt7OT4S4GlANnZ2Zf0S0xMn8gz055hzpA59Inrc0nHEEKISKS01t07gFJLgB8A87TWR8+xz3rgJ1rryvMdKz8/X1dWnncXIYQQZ1BKbdJa53e2zdLNAxcA/wv4VseAV0oNUEpF+V9fBgwDqrpzLiGEEBevu3PyzwFxwDr/GuonWiVnAk8rpdyAD7hPa324m+cSQghxkbrbXdPpo4+01q8Dr3fn2EIIIbqvW9M1QgghQpuEvBBCRDAJeSGEiGARs0DZU+/8g60HzmrDF0KIsDBqcAo/v250wI8rI3khhIhgETOS74m/AYUQItzJSF4IISKYhLwQQkQwCXkhhIhgEvJCCBHBJOSFECKCRUx3jRBChB1POzTugkPbICEVvjEn4KeQkBdCiJ7mdcPhKqjfBoe+OvW1cRf4PMY+I66VkBdCiJDgPg4NO8C5v/PtnuPQsAvqtxph3rATfG7/RgX9cmFgHoy4xvg6YCT0H94jpUrICyHEuXjajNF2xxF4/TZo2gPad+Gf75tjhPiwhaeHeWxiz9fuJyEvhDCP9qNweLcxF34m7YOWr8+YTtkN2mtsV1GQehmkj4IxN8LAkdA3Fyyd9K9YoqHfUIhL7tFfpysk5IUQkcd9zJhOqf/KuKl54mvTXuBCz7VWkDoUBuRB3nUwcJR/BD4MouN6o/qAkpAXQoS+Nhcca+582/GW00ffZ06nWGIg7XIYPAHG3Q4DhkNMUufHsqYb0ykxCT3yawRDt0JeKfUvwPeBQ/63Htdar/Fvewy4B/ACP9RaO7pzLiGECbS3+oO6wwi8fhs4ay78syrKCPNBY2DszcZ0yoA8SPsGRMX0fO0hKhAj+V9rrX/Z8Q2l1CjgVmA0MBgoVkoN1/rE5JYQwhQ87XQ6PeJtN+a7O97MPLQNmved2icq1hhV50w1pkuSBwLq7GPFJhrb0y4Py+mUntZT0zXXA69orduAPUqpXcAUYGMPnU8IEUzHnXBou3/0ve3U1Imr9sI/a4kx5rsz82HCIiOwB+YZNy6jZEa5uwJxBR9QSi0GKoEfa62bgEygtMM+Nf73zqKUWgosBcjOzg5AOUKIHtN25PQwPzG10nE6JToBBoyAy2Yb3SiWqLOPo6JO9YqnXmbq6ZSedsGQV0oVA4M62fQEsAp4BuPfY88A/w7cfTEFaK1XA6sB8vPzL3TbWwjRG9pb/WH+1elh3tJxOiXOuImZc7Ux/32iC6VvTudthSIoLhjyWuv5XTmQUur3wLv+b/cDQzpszvK/J4QIhuMtZ7cTNlad+kh9R9oLR+o5OZceFQtpw2DIZJi42Bh9D8wzRuKdjdJFSOlud02G1vrEpNt3gC/9r98G/qyU+hXGjddhQHl3ziWE6II2lzECP3ONlI4fv49JNKZTsq86943KPkNOdaekXiZz42Gsu//P/T+l1HiMv/KrgR8AaK3/oZT6C7AV8AD/LJ01QlyCo4dPdZ7Uf2UEdpvrHPs2Gp/YPCE6wZhOyZ1xKrAHjoQ+2TKdYiJK69CZBs/Pz9eVlZXBLkOI3nes6ezplPqvoLX+1D6xViOkE1I7P0Z8yqnOlIF5/rlxmU4xA6XUJq11fmfb5N9gQvSEI4dOrUBYv+0CqxXuPL3VMCbJmE4ZtuD00E7JBNVJn7gQ5yEhL0R3tDae3U54aJsxdXJCfF/om935qNoSDUNnnQryASON+XCZThEBIiEvRFd0ZTolzj9dMvLaU+2EA/MgOV1G4CJoJOSF+WhttAieCOrmvZ2vDe51G8vS1n8FR+pOvR+bbAT48IWnbmYOyIOUwRLmIuRIyIvI1tpw9gMfDm0zRuYnxCR13iKoLEYv+DfmntGdMkTCXIQNCXkRPrQ2blDWbzN6wd1HO9sJXAdPhfrRhlOb4voYIZ33rVPz3zKdIiKchLwIPVrDkYNnf6Cn/itoa7nwz59oNRxReHqYWzMkzIXpSMiLnufr8Fi1I/Wd79PeCg3bT93U7DidkpBqhPTYm06F9oCREN+n82NFxUiYC+EnIS8CR2ujH/zMLpRD26H9yIV/Pr6PMe896tunj8CTBkhoC3GJJOTFxdMaXHVntBP6w7zNeWq/pIHGtMn4O07duOyTSacPfoiOkzAXogdIyItT0yktNXT6FB9Pm/GpzI6hfrzD3HhimhHgV9xyavQ9IA+S0nrtVxBCdE5C3ky0NoL8zIceH9oO7tYL/3xCPyO8x9x4en948oCer10IcUkk5MOZz2uMwD3tnWzUxrbT5se3Q3uHFQyT042R90T/I9fOtT64JRpSv2E8Y1OmU4QIKxLy4cDnMz6VeeYIvGGHscDVhST2N6ZQxt16+hopiedYzVAIETEk5IPB6zl3t0mb8/TRd/1WI8w7fvAnJdMI6aEzjafZxyZ1fqzkdH93Sv/A/w5CiLAgId+TfF44vOfsRa0ad4K3symWM1gzjDCftKTDDc0R5+4PF0KIM0jIB4LPB83Vp7cT1n9ljMC9baf265tt3KgcNt8YZXfWShibCP1H+B8O0a+3fgMhRITq7jNeXwVG+L/tCzRrrccrpXKBbcB2/7ZSrfV93TlXSPD5jKfVn/Vhnx3gOXZqv5QsI6Qvm3WqnXDACIhLDl7tQghT6lbIa62/e+K1UurfgY4Li+zWWo/vzvF7nNcDTXvOXiOl41N6OnIfPz3MT0yn5N/ln04Z5Z9OSemd+oUQ4gICMl2jlFLALcDcQBwv4HxeaKo+e7Grs6ZTcoyRd87VxjKzZ4qKhbTLT3WnJPTtrd9ACCEuSaDm5GcAB7XWOzu8N1Qp9RngBJ7UWn/U2Q8qpZYCSwGys7Mv7ewH/wGv3d35tpO95B1aDftkG9Mp35jTYcGrEefuUhFCiDB1wZBXShUDgzrZ9ITW+i3/69uAlztsqwWytdaNSqlJwJtKqdFaa+eZB9FarwZWA+Tn53fymfouiI43QvpchttO706Js17SaYQQItxcMOS11vPPt10pFQ3cAEzq8DNtQJv/9Sal1G5gOFDZrWrPJe0bcMuLPXJoIYQIZ4F4JPx84Cutdc2JN5RSA5RSUf7XlwHDgKoAnEsIIcRFCMSc/K2cPlUDMBN4WinlBnzAfVrrwwE4lxBCiIvQ7ZDXWi/p5L3Xgde7e2whhBDdE4jpGiGEECFKQl4IISKYhLwQQkQwCXkhhIhgEvJCCBHBlNaX9iHTnqCUOgTs7cYh+gMNASonEsj1OJ1cj7PJNTlduF6PHK11pw9bDqmQ7y6lVKXWOj/YdYQKuR6nk+txNrkmp4vE6yHTNUIIEcEk5IUQIoJFWsivDnYBIUaux+nkepxNrsnpIu56RNScvBBCiNNF2kheCCFEBxLyQggRwcIu5JVSf1BK1SulvjzHdqWU+q1SapdS6nOl1MTerrE3deF63OG/Dl8opTYopcb1do297ULXpMN+k5VSHqXUTb1VWzB05XoopWYrpTYrpf6hlPqwN+vrbV34b6aPUuodpdQW//W4q7drDKSwC3ngj0DBebYXYjykZBjGs2NX9UJNwfRHzn899gCztNZjgWeIwBtLnfgj578m+B9q8wtgbW8UFGR/5DzXQynVF1gJfEtrPRq4uXfKCpo/cv4/H/8MbNVajwNmA/+ulIrthbp6RNiFvNb678D5HkByPfCiNpQCfZVSGb1TXe+70PXQWm/QWjf5vy0FsnqlsCDqwp8RgAcxnnlQ3/MVBVcXrsftwN+01vv8+0f0NenC9dCAVSmlgGT/vp7eqK0nhF3Id0Em8HWH72v87wm4BygKdhHBppTKBL5D5P8rr6uGA/2UUuuVUpuUUouDXVCQPQfkAQeAL4CHtNa+4JZ06QLx+D8RBpRSczBCfnqwawkB/wEs01r7jMGa6UUDk4B5QAKwUSlVqrXeEdyygsYGbAbmAt8A1imlPtJaO4Na1SWKxJDfDwzp8H2W/z3TUkpdAfwnUKi1bgx2PSEgH3jFH/D9gWuUUh6t9ZtBrSp4aoBGrXUr0KqU+jswDjBryN8FPKuNDxHtUkrtAUYC5cEt69JE4nTN28Bif5fNVUCL1ro22EUFi1IqG/gbsMjEI7PTaK2Haq1ztda5wGvA/SYOeIC3gOlKqWilVCJwJbAtyDUF0z6Mf9WglEoHRgBVQa2oG8JuJK+Uehnjjnd/pVQN8HMgBkBr/TywBrgG2AUcxfhbOWJ14Xr8DEgDVvpHrp5IW2XvTF24JqZyoeuhtd6mlLIDnwM+4D+11udtPw1nXfjz8QzwR6XUF4DCmNoLx+WHAVnWQAghIlokTtcIIYTwk5AXQogIJiEvhBARTEJeCCEimIS8EEJEMAl5IYSIYBLyQggRwf4/bSY58uKQIREAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0, 0, 0, 0, -34, -76, -25, 0]\n",
      "[0, 0, 0, 0, -34, -42, 88, 113]\n"
     ]
    }
   ],
   "source": [
    "def f():\n",
    "    # Phendrana Drifts has a filter sweep from 0.75 KHz to 4.5 KHz.\n",
    "    # This corresponds approximately to this region in the parameter space.\n",
    "    coeffs = np.linspace(1.0, 1.89)\n",
    "    \n",
    "    # There's a second (quieter) filter sweep from 4.5 KHz to just over 8 KHz.\n",
    "    # I could try making a 4-6 pole filter for this,\n",
    "    # but I'm afraid all the poles/zeros will distort the volume balance.\n",
    "    resonance = 0.9\n",
    "    \n",
    "    fir_coeff_amplitude = np.array([four_pole_fir(coeff, resonance) for coeff in coeffs]).T\n",
    "    \n",
    "    for fir_index in range(4, 8):\n",
    "        coeff_amplitude = fir_coeff_amplitude[fir_index]\n",
    "        plt.plot(coeffs, coeff_amplitude, label=f\"fir{fir_index}\")\n",
    "    plt.legend()\n",
    "    plt.show()\n",
    "    \n",
    "    print(four_pole_fir(coeffs[0], resonance))\n",
    "    print(four_pole_fir(coeffs[-1], resonance))\n",
    "f()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.9.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}