How to reconstruct filter parameters

itu_r_bs_1770_4_kfilter.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "b9568aa6",
   "metadata": {},
   "source": [
    "# How to reconstruct filter parameters\n",
    "\n",
    "The [RECOMMENDATION ITU-R BS.1770-4](https://www.itu.int/rec/R-REC-BS.1770) defines a pre-filter as two biquads.\n",
    "The coefficients for stage1 and stage2 at a samplerate of 48kHz are given on pages 4-5. We translate those tables to Python as:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "5b1f8379",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "stage1 = np.array([\n",
    "    1.53512485958697, # b_0\n",
    "    -2.69169618940638, # b_1\n",
    "    1.19839281085285, # b_2\n",
    "    1.0,\n",
    "    -1.69065929318241, # a_1\n",
    "    0.73248077421585, # a_2\n",
    "])\n",
    "\n",
    "stage2 = np.array([\n",
    "    1.0, # b_0\n",
    "    -2.0, # b_1\n",
    "    1.0, # b_2\n",
    "    1.0,\n",
    "    -1.99004745483398, # a_1\n",
    "    0.99007225036621, # a_2\n",
    "])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c8570625",
   "metadata": {},
   "source": [
    "For other samplerates then `samplerate = 48000` we need to construct filters with a similar frequency response."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cc2d9a1f",
   "metadata": {},
   "source": [
    "## First attempt"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "60379ab0",
   "metadata": {},
   "source": [
    "How to calculate the coefficients of digital filters is given in\n",
    "\"DAFX: Digital Audio Effects\" (Udo Zölzer 2011, 2. edition).\n",
    "Table 2.3 on page 64 defines coefficients for second-order high-frequency boost shelving filters (`stage1`) and\n",
    "Table 2.2 on page 50 defines coefficients for second-order highpass filters (`stage2`).\n",
    "\n",
    "We invent two functions to calculate the second-order-sections:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "9c048003",
   "metadata": {},
   "outputs": [],
   "source": [
    "def hf_boost_shelv(K, V_0):\n",
    "    \"\"\"Calculate coefficients for a second-order high-frequency boost shelving filter.\"\"\"\n",
    "    denominator = (1 + np.sqrt(2)*K + K**2)\n",
    "    b_0 = (V_0 + np.sqrt(2*V_0)*K + K**2) / denominator\n",
    "    b_1 = 2*(K**2 - V_0) / denominator\n",
    "    b_2 = (V_0 - np.sqrt(2*V_0)*K + K**2) / denominator     \n",
    "    \n",
    "    a_1 = 2*(K**2 - 1) / denominator\n",
    "    a_2 = (1 - np.sqrt(2)*K + K**2) / denominator\n",
    "\n",
    "    return np.array([b_0, b_1, b_2, 1.0, a_1, a_2])\n",
    "\n",
    "\n",
    "def highpass(K, Q):\n",
    "    \"\"\"Calculate coefficients for a second-order highpass filter.\"\"\"\n",
    "    denominator = ((K**2) * Q + K + Q)\n",
    "    b_0 = Q / denominator\n",
    "    b_1 = -2*Q / denominator\n",
    "    b_2 = b_0\n",
    "    \n",
    "    a_1 = (2*Q * (K**2 - 1)) / denominator\n",
    "    a_2 = ((K**2) * Q - K + Q) / denominator\n",
    "\n",
    "    return np.array([b_0, b_1, b_2, 1.0, a_1, a_2])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "179085ff",
   "metadata": {},
   "source": [
    "`K` is defined as `K = np.tan(np.pi * f_c / samplerate)`, thus"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "ce33abec",
   "metadata": {},
   "outputs": [],
   "source": [
    "def cutoff_freq(K, samplerate):\n",
    "    \"\"\"Return the cutoff frequency in Hertz.\"\"\"\n",
    "    return (samplerate * np.arctan(K))/np.pi"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f857eceb",
   "metadata": {},
   "source": [
    "`V_0` is defined as `V_0 = np.power(10.0, G / 20.0)`, where `G` is a gain we are not interested in."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "25985191",
   "metadata": {},
   "source": [
    "Now we can try to find values for `K`, `Q` and `V_0`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "51da6f3c",
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.optimize import minimize\n",
    "\n",
    "def solve_stage1(samplerate):\n",
    "    \"\"\"Minimize the difference between target biquad and calculated biquad.\"\"\"\n",
    "    def helper1(z):\n",
    "        K = z[0]\n",
    "        V_0 = z[1]\n",
    "        sos = hf_boost_shelv(K, V_0)\n",
    "        \n",
    "        result = np.sum(np.square(stage1 - sos))\n",
    "        return result\n",
    "    \n",
    "    guess = np.array([1.0, 1.0], dtype=np.float64)\n",
    "    result = minimize(helper1, guess, tol=1e-11)\n",
    "    return result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "a907cdaa",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "      fun: 1.2547183412375365e-09\n",
      " hess_inv: array([[0.03148047, 0.0707926 ],\n",
      "       [0.0707926 , 0.27269465]])\n",
      "      jac: array([ 7.18314297e-14, -5.01126918e-14])\n",
      "  message: 'Optimization terminated successfully.'\n",
      "     nfev: 36\n",
      "      nit: 9\n",
      "     njev: 12\n",
      "   status: 0\n",
      "  success: True\n",
      "        x: array([0.11051784, 1.58485266])\n",
      "freq = 1681.7632251028442 Hz\n",
      "V_0 = 3.9997778685513232 dB\n"
     ]
    }
   ],
   "source": [
    "samplerate = 48000\n",
    "\n",
    "result = solve_stage1(samplerate)\n",
    "print(result)\n",
    "\n",
    "K1, V_0 = result.x\n",
    "\n",
    "print(f'freq = {cutoff_freq(K1, samplerate)} Hz')\n",
    "print(f'V_0 = {20.0 * np.log10(np.abs(V_0))} dB')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "140e6974",
   "metadata": {},
   "source": [
    "Let's look at the coefficient errors:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "757afe31",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[-1.74123774e-05  2.54301517e-06  2.51807174e-05  0.00000000e+00\n",
      "  1.65097513e-05 -6.19839611e-06]\n"
     ]
    }
   ],
   "source": [
    "err = stage1 - hf_boost_shelv(*result.x)\n",
    "print(err)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eefc3ddc",
   "metadata": {},
   "source": [
    "Okay.\n",
    "\n",
    "Now we repeat that for the second stage."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "38d47933",
   "metadata": {},
   "outputs": [],
   "source": [
    "def solve_stage2(samplerate):\n",
    "    \"\"\"Minimize the difference between target biquad and calculated biquad.\"\"\"\n",
    "    def helper2(z):\n",
    "        K = z[0]\n",
    "        Q = z[1]\n",
    "        sos = highpass(K, Q)\n",
    "        result = np.sum(np.square(stage2 - sos))\n",
    "        return result\n",
    "    \n",
    "    guess = np.array([0.002, 0.25], dtype=np.float64)\n",
    "    result = minimize(helper2, guess, tol=1e-09)\n",
    "    return result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "c6c1af49",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "      fun: 8.469157205235147e-05\n",
      " hess_inv: array([[0.00227113, 0.01264941],\n",
      "       [0.01264941, 1.00309113]])\n",
      "      jac: array([9.96333256e-08, 9.05220077e-09])\n",
      "  message: 'Desired error not necessarily achieved due to precision loss.'\n",
      "     nfev: 159\n",
      "      nit: 5\n",
      "     njev: 49\n",
      "   status: 2\n",
      "  success: False\n",
      "        x: array([0.00070414, 0.24727469])\n",
      "freq = 10.758482410249737 Hz\n",
      "Q = 0.2472746919148863\n"
     ]
    }
   ],
   "source": [
    "result = solve_stage2(samplerate)\n",
    "print(result)\n",
    "\n",
    "K2, Q = result.x\n",
    "print(f'freq = {cutoff_freq(K2, samplerate)} Hz')\n",
    "print(f'Q = {Q}')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "95887f38",
   "metadata": {},
   "source": [
    "The errors:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "deb8ab24",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 0.00284001 -0.00568003  0.00284001  0.          0.00427153 -0.00424871]\n"
     ]
    }
   ],
   "source": [
    "err = stage2 - highpass(*result.x)\n",
    "print(err)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a474ffba",
   "metadata": {},
   "source": [
    "Now let's have a look at the transfer functions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "a95cfb1b",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/tmp/ipykernel_28555/4003190015.py:16: RuntimeWarning: divide by zero encountered in log10\n",
      "  ax1.plot(w, 20 * np.log10(abs(h)), label='EBU')\n",
      "/tmp/ipykernel_28555/4003190015.py:17: RuntimeWarning: divide by zero encountered in log10\n",
      "  ax1.plot(w2, 20 * np.log10(abs(h2)), label='jkbd')\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYcAAAEaCAYAAAD65pvjAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAA330lEQVR4nO3dd5xU5fX48c/ZXmFhl74sIFJFpCzYgiBBY0eNsQQLNtREExNb1JivmpiYxCTKz4oFjTVERcCKqNjpRTrSWeqysL3vnN8f9y4ZtzG7O7MzO3ver9e8dm577pm5cM88z3Pvc0VVMcYYY7xFBDsAY4wxoceSgzHGmFosORhjjKnFkoMxxphaLDkYY4ypxZKDMcaYWiw5GABE5GkRuc/f69azvYrI0fUs6yIiX4hIgYj8Q0TuEZHn3GW93W2jmrpvH2L7wf4DtR9jQl3A/pOZ0CEi24AuQCVQBawF/g1MU1UPgKre6Gt53uuKyDjgFVVN91O4U4ADQDs9wk04IjLf3fdzftp3o/ZvTDizmkPbca6qJgO9gIeBu4DngxtSnXoBa1vixCwikY3dfyBrLaGonu/ItAWqaq8wfwHbgAk15o0GPMAQd/pF4E9ey+8E9gC7gesABY72XhdIBErccgrdV3e37G+BXLeMx4EYr7IPl1UjpheBCqDcLWsCcD9O7QCgt7ttFPAQTi2o1F33cXedgcDHwEFgA3BxjfKfAt4Hiur4Turb/5vAK0C++120x0mse4Bd7ncR6ZYRCTyCU/vYAvyyOua6joX353OnTwC+cb+7lcA4r2XzgT8CXwMFwFwgzWv5j7y23QlMBkYB+6r37673U2BFPf9Wan1H7jF9C8gGtgK/qvHvaIn73ewD/lnjWE3B+Te0B7jNa7tY4FF32W73fay7bByQBdwG7He3vdpr27Nwar8F7vd/u9eyc4AV7nfwDTA02P//Wusr6AHYqwUOch3JwZ2/A7jJff8ibnIAzgD2AscACcDL1JEc3PfjgKwa5Y50T3JR7kliHXCr1/I6k0PNst3p+6kjObjT84HrvNZNdE+KV7v7HoFzkj7Gq+w84GScWnOcj/uvAM53t4kH3gGecffXGVgE3OCufyOwHugJdAQ+w8fkAPQActyTXwRwmjvdyevzbgb6u3HMBx52l2XgnCwvA6KBVGCYu2wtcKbXPmfidaKu4/N7f0cJwFLgD0AMcBRO0vuJu/63wBXu+yTghBrH6nX3ezoWJ7lMcJc/CCxwv79OOCfyP3r9m6p014l2v49ioIO7fA8wxn3fARjhvh+Bk0yOx0nSV7nfd2yw/w+2xpc1K7Vtu3FOYDVdDExX1TWqWgw80JhCVXWpqi5Q1UpV3YZzIh3b7GiP7Bxgm6pOd/e9DOcX70Ve68xS1a9V1aOqpT6W+62qvqNO/0w74EycZFekqvuBfwGXuuteDDyqqjtV9SDwl0bEfznwvqq+78b3Mc6v8rO81pmuqhtVtQSYAQxz508C5qnq66paoao5qrrCXfaSWzYi0hH4CfBaA3Ec/o5wTuqdVPVBVS1X1S3As16ftwI4WkTSVLVQVRfUKOsB93taBUzHSV7V8T6oqvtVNRvn39gVXttVuMsrVPV9nJrcAK9lg0Wknaoeco8zwPXAM6q6UFWrVPUloAznh4ppJEsObVsPnOaXmrrj/AKvtrOOdeolIv1F5F0R2Ssi+cCfgbSmh+mzXsDxIpJb/cI5CXX1WqdRn6WObXrh/Jrd47WPZ3B+AUPt7257I/bTC/hZjfh/BHTzWmev1/tinF/r4NRUNtdT7ivAuSKShJO8vlTVPQ3EUfPzdq8R0z04FzgAXItTk1kvIotF5JwGytqO8/3g/t1ezzKAHFWt9Jr2/qw/xUmY20XkcxE50SvW22rE2rNGucZHbapzzfyPiIzCSQ5f1bF4D+B99VHPBoqqq+P2KWA5cJmqFojIrfzw17u/1Nz3TuBzVT2tEds0dj87cX6NptU4eVXbww+/r4way4twmmqq1UxcL6vq9U2IcSdO+38tqrpLRL4FLsD5df7UEcqq+Xm3qmq/esr+HrhMRCKAC4E3RSTVa5WeOM1s4HwXu933u3FO5mvqWNZwcKqLgYkiEg3cjFOD6unG+pCqPuRLOaZhVnNoY0Sknfvr7g2ctu5Vdaw2A7haRAaJSAJOe3N99gGpItLea14yTgdloYgMBG7yU/h17fsor+l3gf4icoWIRLuvUSIyyF87dH9xzwX+4X6XESLSV0Sqm81mAL8SkXQR6QD8rkYRK4BL3dgy+WHSrP6F/xMRiRSROBEZJyK+XCb8KjBBRC4WkSgRSRWRYV7L/41zkcGxOH0OvloE5IvIXSIS78Y1xP1xgYhcLiKd3CaoXHebKq/t7xORBBE5Bqcv6D/u/NeB34tIJxFJw/k39sqRghGRGBGZJCLtVbUC599Z9f6eBW4UkePFkSgiZ4tIciM+r3FZcmg75ohIAc6vq3uBf+L8Z61FVT8ApuJ0pm7C6XQE5xdzzXXX4/xH3+JW5bsDtwM/x+kgfZb/nRD87THgIhE5JCJTVbUAOB2nPXw3ThPMX3GujPGnK3E6Z9cCh3CuZqpu+nkW+AjnSqNlwNs1tr0P6Otu9wBebf+quhOYiNNsk41zrO7Ah/+nqroDp6nlNpymwhXAcV6rzMT5pT5TVYt8/aCqWgWci9O3sRWng/85nCu2wLl4YY2IFOIcj0tr9OV8jvNv6BPgEVWd687/E05/ynfAKpzv6k8+hnUFsM1tsrwRtz9FVZfg9Ds8jvP9bsK5Yss0gajafT6mYe4v79U4V33U1ZRi6iEivXFOqtHB/u5EZDPOVVXzWmBfvQmRz22axmoOpk4icoFbhe+A8+t7jv0nb71E5Kc4fQmfBjsW0zqEbHIQkRQReVNE1ovIOq8rEkzLuAGnaWMzTptuoPoNTIC5w4w8BfzS7Rsw5ohCtllJRF7CueTuORGJARJUNTfIYRljTJsQkslBRNrhdOgdpaEYoDHGhLlQbVY6CqdJY7qILBeR50QkMdhBGWNMWxGqNYdMnHFXTlbVhSLyGJCvqvd5rTMFZ1AvEhMTRw4cODA4wRpjTCu1dOnSA6raqa5loZocugILVLW3Oz0G+J2qnl3X+pmZmbpkyZIWjNAYY1o/EVmqqpl1LQvJZiVV3QvsFJHqgbZ+jHPDkTHGmBYQymMr3QK86l6ptIV67uY1xhjjfyGbHNzhhuus7hhjjAmskE0OzVVRUUFWVhalpb4O2d/6xMXFkZ6eTnR0dLBDMcaEmbBNDllZWSQnJ9O7d29EJNjh+J2qkpOTQ1ZWFn369Al2OMaYMBOSHdL+UFpaSmpqalgmBgARITU1NaxrRsaY4Anb5ACEbWKoFu6fzxgTPGGdHIItMjKSYcOGHX49/PDDAIwbN44BAwYwbNgwBg0axLRp0w5vk5SU9IMyXnzxRW6++eYWjdsYY8K2zyEUxMfHs2LFijqXvfrqq2RmZnLw4EH69u3L5MmTiYmJadkAjTGmHlZzCLLCwkISExOJjIwMdijGGHNYm6g5PDBnDWt35/u1zMHd2/F/5x7T4DolJSUMGzbs8PTdd9/NJZdcAsCkSZOIjY3l+++/59FHH7XkYIwJKW0iOQSLL81K2dnZnHTSSZxxxhn06tWrznWt49kY09LaRHI40i/8YOrUqRMjRoxg4cKF9OrVi/j4eMrLyw/3Pxw8eJC0tLQgR2mMaWuszyHIiouLWb58OX379gVg7NixvPLKK4DTLDVjxgxOPfXUYIZojGmD2kTNIVhq9jmcccYZhy9nnTRpEvHx8ZSVlTF58mRGjhwJwGOPPcYNN9zA1KlTUVWuvPJKTjnllGCEb4xpwyw5BFBVVVWd8+fPn1/vNj169ODdd98NUETGGOMba1YyxhhTiyUHY4wxtYRschCRSBFZLiLWxmKMMS0sZJMD8GtgXbCDMMaYtigkO6RFJB04G3gI+G2QwzHGhDBVpazSQ1FZJcXlVRSXV1FWWUVFlVJR5aHS/eu8lEqPh/JKD5UepdKjqCoej+JRULc8j7rTCh5116lrGudvHUGhqohTIgCiHH6v1LVR9bZHni01VlKtv7wGFjUoJJMD8ChwJ5Ac5DiMMS2srLKKvXml7MotYdehErILyygqyKc8P5vyokN4ig+hpXlEV+QTV1lATFUxsZQT575ipfp9hTtdQSIeIqki6vDfKiLwONPyv/mReIjEg6Du638n4uppaiz733KIkCaeiYPkngaWhVxyEJFzgP2qulRExjWw3hRgCkBGRkbLBNcEJ510En/+85955JFHal2iev/995OUlMTtt9/eYBmTJ0/mnHPO4aKLLgpkqMa0qLySCtbuzuf7fflk7c6ieM9GonO3kFq2k26SQxcOMUxy6SKHaCfFdRciQBRUSRRVkXFURcaikXF4ouLQyDg0Kg6i2kNEFEREIZHRREREIpFRSKQ7z30REQmRUYhEIBIBIs7QNe5fgR/Ox1n2v7/UMc/rb3W8+Gk4HH8Mq/NA/eeekEsOwMnAeSJyFhAHtBORV1T1cu+VVHUaMA0gMzMzZNP1N9980+B9Dca0FXvzSvni+2w2f7+eih2LSStYx7GyhXMittFRCg+v54mKpCS+C1VJXYhsN4zYDj2gfTdITIO4FIhr77zi3fcxyURGRmFDVzZFK0oOqno3cDeAW3O4vWZiaE2SkpJ+UGNYvHgxU6ZM4a233gJg5cqVjB8/np07d3LnnXdy/fXXo6rccsstfPrpp/Tp06fB9kRjQpWqsmZ3Ph+v2Ezh2o/pnbeIkyNWc3HEXgCqoiMpThlIZI+JaPfBSFo/SD2aiJQMEiOjgxy9CbnkEBAf/A72rvJvmV2PhTMfbtQm33zzDbfccguzZs063BT23XffsWDBAoqKihg+fDhnn302CxYsYMOGDaxatYp9+/YxePBgrrnmGv/Gb0yA7Msv5Z1F35O9dBYjC+dzU8QK4qSCith4ynqchGfgLURknEhkl2NIjo4LdrimHiGdHFR1PjA/yGH4xbp165gyZQpz586le/fuh+dPnDiR+Ph44uPjOfXUU1m0aBFffPEFl112GZGRkXTv3p3x48cHMXJjfLMqK485n8yn26bXuTTic9pLMSUJacjgK2DohURnnEC01QhajZBODn7TyF/4gdCtWzdKS0tZvnz5D5JDzWc1VE/bMxxMa/FdVi5vzp7NmD3TuSdyGVVRUZQefRacdD3xvU52OnpNqxPKN8GFlZSUFN577z3uueeeH3RQz5o1i9LSUnJycpg/fz6jRo3ilFNO4Y033qCqqoo9e/bw2WefBS9wY+qxPaeIPz0/g5xnzuPB/bcwJnYzpWN+R+Rt60ic9DL0OcUSQyvWNmoOQeRdA+jSpQtz5szhzDPP5IUXXgBg9OjRnH322ezYsYP77ruP7t27c8EFF/Dpp59y7LHH0r9/f8aOHRus8I2ppayyin/PW0biN3/lbplHeWw7yk6+j7iTboBYuzUpXFhyCKCcnBw6duzIuHHjGDduHODck7FmzRoAjj/++Dq3ExEef/zxlgrTGJ+t2HGID157lJtKniU5ooTSYdeS+JPfQ3yHYIdm/MySQ4Ds3r2bcePGHfEGN2NagyqPMv2jhfT69l7ujlhCbqeRRF78BImdBwU7NBMglhwCpHv37mzcuDHYYRjTbPvyS3lu+rPcdPBhkiPLKBn3ACljbrH+hDBnycEYU6+lWw+w9OW7ubvqvxS070f0Fa8Q3WlAsMMyLSCsk4OqhvUloXbntAmkt79ZTdqHv2BKxEryBlxE+4v+H8QkBDss00LCNjnExcWRk5NDampqWCYIVSUnJ4e4OLvD1PiXx6M8PftTTlt2C30i9lF8+iO0P/E6/wz0ZlqNsE0O6enpZGVlkZ2dHexQAiYuLo709PRgh2HCSEWVh8df/g+Xb72LpGgP8vOZJPQ9JdhhmSAI2+QQHR1Nnz59gh2GMa1GUVklTz33FL/c/yAVcWnEXTsT6Tww2GGZIAnb5GCM8V1ecQWPPz2VO/IeojBlAB2vnwVJnYMdlgkiSw7GtHG5xeU88eS/uLPgrxSlDqHj9bOdZyWYNs2SgzFt2KGicp588h/cVfg3itKOI+X6Wc4DdEybZ8nBmDbqYFE5zzz5d+4qfISizsNpf90sGxvJHGbJwZg2KKewjGef/Ct3Fv2Twi6ZtL/2HYhNCnZYJoSE5JDdItJTRD4TkXUiskZEfh3smIwJFwcKy3j+ib9wR9E/Kew62q0xWGIwPxSqNYdK4DZVXSYiycBSEflYVdcGOzBjWrMDhWVMf+Ihbi+eSn63E0m55i2769nUKSRrDqq6R1WXue8LgHVAj+BGZUzrdqCwjH8//iC3FU8lv/vJpFz7tiUGU6+QTA7eRKQ3MBxYWGP+FBFZIiJLwvkuaGP84UBhGa88fj+/LX2cvB6nkHLNmxAdH+ywTAgL6eQgIknAW8CtqprvvUxVp6lqpqpmdurUKTgBGtMKHCgs4/XH7+PW0ic51GMcHa6eYYnBHFGo9jkgItE4ieFVVX072PEY0xodKCxjxuP3cEvpcxxM/zEdJ78OUbHBDsu0AiGZHMQZRvV5YJ2q/jPY8RjTGh0oLOOdx+/kF6UvktPzdFKvehWiYoIdlmklQrVZ6WTgCmC8iKxwX2cFOyhjWovduSXMnnor15W+SHbvc0id/JolBtMoIVlzUNWvABs83pgm2JpdyOdP38o1Vf/lQN8L6DTpeXukp2m0kEwOxpimWbfrECue+wWT9X0ODriUtEuetMRgmsSSgzFh4tuNuyh47RouYwG5x11Px4l/g4hQbTk2oc6SgzFhYOa3a+j+wbWcHrGO3B/9HykTfhvskEwrZ8nBmFbM41FemP0xY5fd6jzv+dynSRl5WbDDMmHAkoMxrdShonJeeGka1+97iMjoaPj52yT0HRvssEyYsORgTCu0fNsBFr1yH7+peJ3cdv1JvmYG0qF3sMMyYcSSgzGtSEWVh1c++prBC+7ghoh15PY9j46XPg0xicEOzYQZSw7GtBIb9uQz59WpXFfwJPGRVRSfMZWU0VeC2C1Bxv8sORgT4vJLK3j53c847rsHuT1iFbmpxxE7aTqk9g12aCaMWXIwJkSVV3qYtWAdBZ88wrWeORAVTdG4v5DyoxvsxjYTcJYcjAkxZZVVzFy4kQOfPcWkirfoIIXkHj2RlIl/hXbdgh2eaSMsORgTInblljD7y6XELX+eC6o+IkWKyOl+CnruH0npPizY4Zk2xpKDMUGUV1LBvNU7yVo0m4F73+XaiGVEiYeDGaejp/2W1IwTgh2iaaMsORjTglSVTfsLWbB+O3mrP6LLvi8YL0tJlQKK4jpQNuRaYn50I2nW2WyCzJKDMQFUUFrBmt35rN+WRdHmb4nbs5jBFWu5JGIDMVJFaXQSJb3Go8dPIvHoCRBp/yVNaAjZf4kicgbwGBAJPKeqDwc5JGPqVF7pYX9BKTsOFrNjTza5ezbh2b+BhNyNdCndQn/JYrTsJUKUKiLI6zCA0r7XE3PcucT1PJ64yOhgfwRjagnJ5CAikcATwGlAFrBYRGar6trgRmbaAlUlv7SSvOIKckvKyS0qpyg/h9K8/VTkZ1OWv5+qgmwoziGm9AAdK/fTQ7IZJAc4SQoPl+MhgtyknlSmDaUk40oS+55MZI+RdIxNCuKnM8Y3IZkcgNHAJlXdAiAibwATAUsO5sg8VXjKCiksyKOoII+iwjxKi/IoK8qnvDifypICqkoL8JQVQlkhVBQRUV5IVGUx0VXFxHhKSKSERCklg1IGUUq0VNW5q7KIeAqTu1KWmE5pyknkd+5DUpc+RHQaQERafzpGx7XwhzfGP+pNDiIy1Yft81X1936Mp1oPYKfXdBZwvPcKIjIFmAKQkZERgBBMi6ssg9I8KMl1/pbmQWkulOZSVnCQkoKDVBTnU1FSQFVZIZQVQUUhkZXFRFcWE+0pJU5LiKOcCKCd+2pIGTGUSDzlEfGURyZQFZeIJyoNYhIpj03GE5dESUI7opI7EduuM/EpnYlO7gQJaZCQSmxMArEB/2KMaXkN1RwmAn84wva/AwKRHOoaLEZ/MKE6DZgGkJmZqXWsb0JBZRnk74K8XVC4D4qyoXA/FO2HogNQuB8t2o8WHSCisrTeYmIB1WgqiKdC4ygmjiLiKJN4KqO64IlORKMTIToBiU1CYpOIjEsiOj6ZmPgkYhPaE5/UnoTkdiQmdyAxqR0RscnERkbZyd2YOjSUHP6lqi81tLGIdPBzPNWygJ5e0+nA7gDtyzSHx+Oc/HO+h5zNkLsdcndCXpbzKtxbexOJpDAyhUOSwl5PO3aV9+GAHkeeJpJHIoUkIPEdiEnqQHy7VJJSUklJSaNjSjtSk2JJTYyhY2IMfRJjiIu2YSSMCYR6k4OqPnqkjX1Zp4kWA/1EpA+wC7gU+HmA9mV8oeqc+Peucl7ZGyBnk5MQKkv+t15kLNo+nbLEHuxLPYnt7VNZX9KelXmJbChO5IC2J49EEmNj6J2WQO/UROeVlsjw1ATSO8TTKSmWqEh79rExwdRQn0MccAlwCJgD3AmMATYDf1TVA4EKSlUrReRm4COcS1lfUNU1gdqfqUPeLtjxLexa6iaE75w+AACJgJRekNYP+oylqmNfNnm6sTC/I1/tiWR5Vh7Zu8sAiImMoF+XJAYOaMclXZMZ2C2ZAV2T6ZQUi9hQ08aELFGtu7leRGYAFUAi0AFYjZMkfgQMU9VzWirII8nMzNQlS5YEO4zWSxUOfA/bvoAdC2HHAsjb4SyLioMux0DXodBtKHQ9Du08kM25ylffZ/PVphwWbMmhsKwSgN6pCYzI6MDwjBSGZ3RgQNdkoq0WYExIEpGlqppZ17KG+hwGq+oQEYkCslS1+uG0H4rISr9HaVpWaT5s/Rw2feK8qpNBUlfIOB5O/AVknABdjoXIKDweZWVWLh+u2stHqxexLacYgIyOCZx7XHfG9Evj+D4dSU2y7l1jwkFDyaEcDjfx1OwMrvuibxPainJg/buwdpaTGDyVEJMMR42FH90KfU+FDn1+8GSxtbvzeXNpFu+v2sPe/FKiIoQT+6Zy7ZijGNuvExmpCcH7PMaYgGkoOaS79zqI13vc6R4Bj8z4R1kBrHkHVr8JW78ErYIOveHEX0K/n0DP0VBj+IaDReW8s3wXby7NYu2efKIjhXEDOnPnkAH8eGAX2ifYcA/GhLuGksMdXu9rNuhbA38oU4WdC2HZy7BmJlQUQce+Tu1g8ESn/6COzuANewuY/vVWZi7fRVmlh6Hp7Xlw4jGcO7Q7HRJjWv5zGGOCpqFLWRu8x8GEoIoSWPkGLHwastdDTBIMuRBGXAnpo+pMCKrK15tyeOrzTXy9KYe46AguHJHOVSf1YmDXI91fbIwJVw1dyjqHGncle1PV8wISkWm8gn2w+FlY/DyUHHRqBhOfgMHnQz2DvFUnhUfnbWTJ9kN0aRfLnWcM4LJRGVZLMMY02Kz0iPv3QqAr8Io7fRmwLYAxGV8V7IWvHoWl051hKgac5Vxl1OvkOmsJ1ZbvOMSf31/H4m2H6NY+jj+eP4SLM9OJjbK7jY0xjoaalT4HEJE/quopXovmiMgXAY/M1K9wv5MUljwPVRVw3GUw5rdwhKeH7c0r5a8frmfm8l2kJcXy4MRjuGRUT0sKxphafBmyu5OIHOU1fHYfoFNgwzJ1qiiFBU/Cl/+AimIYeimccvsRk0JFlYfnvtzK1E++p0qVX4zryy9OPZqk2FAdsd0YE2y+nB1+A8wXkS3udG/cobJNC1GFdbNh7n3O+EYDzoLTHnSGrziC1bvyuOut71izO5/TB3fhvnMG07Oj3ZtgjGnYEZODqn4oIv2Age6s9apaFtiwzGEHt8K7v4Etn0HnwXDFO87NakdQUeXhsXnf89Tnm+mQEMNTk0Zw5rHdAh+vMSYsNHS10ghVXQbgJoNaQ2Z4r2P8rKoSFj4Fnz4EEVFw5t8h8xqfHkC/82Axt7y+nBU7c/npiHTuO2cQKQl2BZIxxncNnWmmi8g46n7wTrXngeH+DMjgDIc98wbYvRz6nwlnPwLt033adPbK3dz79ioQeOLnIzh7qNUWjDGN11ByaA8speHkkO3fcNo4VVjyAnx0L8QkwEXT4ZgLGrwstVpllYeH3l/H9K+3MSIjhccuHW59C8aYJmvoUtbeLRiHKcqB2TfDhveh73g4/ylI7urTprnF5dz82nK+2nSAq0/uzT1nDbJhso0xzRJy1zKKyN+Bc3FGhd0MXK2quUENKtB2LoYZV0LxAfjJX+D4GyHCt5P75uxCrn1xMbtyS/jbT4dy8aieR97IGGOOIBR/Xn4MDFHVocBG4O4gxxNYS6bD9DOdkVGv+8S5w9nHxLBiZy4XPfUNBaWVvDHlBEsMxhi/Cbmag6rO9ZpcAFwUrFgCqrIM3r8dlv0b+v4YfvocJHT0efMvNmZz4ytLSU2K4eVrjqd3WmIAgzXGtDVH/IkqjstF5A/udIaIjA58aABcA3xQT1xTRGSJiCzJzm5l/eLFB+HfE53EMOZ2mPTfRiWG977bw7UvLaZXaiJv3XiSJQZjjN/5UnN4EvAA44EHgQLgLWBUU3cqIvNwBvOr6V5VneWucy9QCbxaVxmqOg2YBs4zpJsaS4s7uAVe/Rnk7nSuRhpyYaM2f++7PfzqjeWMyEjh+cmjaBdnD94xxvifL8nheFUdISLLAVT1kIg0644qVZ3Q0HIRuQo4B/ixqraeE/+RZC2B1y5xnsZ25SzodWKjNv9wtZMYhvdM4cWrR5NoYyMZYwLEl7NLhYhE4j7bQUQ64dQkAkJEzgDuAsaqanGg9tPivp8H/7kckrvApLcg7ehGbT53zV5ufm05x6W358VrLDEYYwLLl8tipgIzgc4i8hDwFfDnAMb0OJAMfCwiK0Tk6QDuq2Wsexdev9RJCNfOa3RiWLT1IDe/vpxjerTnpWtG22iqxpiA82XgvVdFZCnwY5y7pc9X1XWBCkhVG3fmDHWr3oS3p0D34XD5mxDfoVGbb9xXwHUvLSY9JZ7pk0eRbH0MxpgW0NDAe96Xz+wHXvdepqoHAxlYWFj+Csy62Xky28/fgNjkRm2+J6+Eq15YRGx0JC9dM5qO9vhOY0wLaajmsBSnn0GADOCQ+z4F2AH0CXRwrdrK/ziJoe+pcMmrzlhJjVBcXsk1Ly6hoLSS/9xwgo2TZIxpUfX2OahqH1U9CvgIOFdV01Q1FecqordbKsBWae0seOdG6DMGLn2t0YlBVbn9vyvZsDefJyaN4Jju7QMUqDHG1M2XDulRqvp+9YSqfgCMDVxIrdzGufDmtZA+Ci59HaLjG13E459u4v1Ve7n7zEGM7W9PZDXGtDxfLns5ICK/B17BaWa6HMgJaFSt1dYvnctVuxzj3PUcm9ToIuau2cs/Pt7IBcN7cN0Ya7kzxgSHLzWHy4BOOJezvgN0ducZb/vWwBs/h4594IqZENf4pqAdOcXcNmMlQ9Pb85cLj0V8eI6DMcYEgi+Xsh4Eft0CsbReebucITFiEuHytxo1TlK18koPN7++DBF4ctII4qIjAxCoMcb45ojJQUQ+w7072puqjg9IRK1NaZ6TGErz4ZoPfX6cZ01//XA932Xl8fTlI0jvYFcmGWOCy5c+h9u93scBP8UZEM9Uljt9DAc2ODWGrkOaVMy8tft4/qutXHliL84YYs98NsYEny/NSktrzPpaRD4PUDythyq8fxts/QLOfxqOGtekYvYXlHLHmysZ3K0d95w1yL8xGmNME/nSrOTdgB4BjKTu4bbblkXPus9juA2GNa1/XlW55+3VFJVXMfWyYdbPYIwJGb40K3nfKV0JbAWuDWRQIW/L5/Dh76D/mXDq75tczNvLdjFv3T5+f/Ygju7cuKE1jDEmkHxJDoNUtdR7hojEBiie0HdwK/z3KkjrBxdO8/l5zzXtySvh/jlrGNW7A1efbPczGGNCiy9ntm/qmPetvwNpFcoKnXsZVJ1hMeLaNakYVeWut1ZRWaU88rPjiIyw+xmMMaGloVFZuwI9gHgRGY7TrATQDmh711qqwru/gf3rnCuTUvs2uaiZy3fxxcZsHjjvGHql2vOfjTGhp6FmpZ8Ak4F04J9e8wuAewIYEwAicjvwd6CTqh4I9P6OaNlLsGoGnHovHP3jJheTW1zOQ++tY3hGClec0MuPARpjjP/UmxxU9SXgJRH5qaq+1YIxISI9gdNwhgYPvj3fwft3Qt/xMOb2I6/fgIc/WE9uSQWvXHAsEdacZIwJUQ01K12uqq8AvUXktzWXq+o/69jMX/4F3AnMCuA+fFOa53RAJ6TChc82uQMaYMm2g7yxeCdTTjmKQd2a1l9hjDEtoaFmperG8MYPLdoMInIesEtVVzY08JyITAGmAGRkZAQmGFWYfQsc2g6T34PEtCYXVVHl4d6Zq+nePo5f/7ifH4M0xhj/a6hZ6Rn37wP+3qmIzKPuG+nuxenPOP1IZajqNGAaQGZmZq2xn/xi6XTnwT0THoBeJzarqH9/u50N+wqYdsVIEmN9uYLYGGOCx5c7pDsB1wO9vddX1WuaulNVnVDPvo7Fefxoda0hHVgmIqNVdW9T99ck2Rvhw3vgqFPhpF81q6iDReU8Nm8jY/qlcdrgLn4K0BhjAseXn7CzgC+BeUBVIINR1VU4z4sAQES2AZktfrVSZTm8fZ3zFLfzn2pWPwPAo/M2UlRexX3nDLZnNBhjWgVfkkOCqt4V8EhCyfw/w56VcMkr0K55o6Ru3FfAqwt3MOn4DPp3sSEyjDGtgy8/id8VkbMCHkkdVLV3i9catn0FXz0Kw6+AQec2qyhV5Y/vriUxJpJbJ/T3T3zGGNMCfEkOv8ZJECUiki8iBSKSH+jAgqIkF2be6Dzq84yHm13c/I3ZfPn9AX49oT8dE2OaH58xxrQQX57n0HbaQj66B/J3w7VzIbZ5V/B6PMrfPtxAr9QEuxPaGNPq+HK10og6ZucB21U1fJ4It3EurHjVuQM6PbPZxc35bjfr9uTz2KXDiIlqXoe2Mca0NF86pJ8ERgCr3OljgZVAqojcqKpzAxVciynNgzm/hk6DYOydzS6uosrDvz7eyMCuyZw7tLsfAjTGmJbly0/abcBwVR2pqiOBYcBqYALwt8CF1oI+uhcK98L5T0BU8x9V8d8lWWzLKeaOnwyw8ZOMMa2SL8lhoKquqZ5Q1bU4yWJL4MJqQZvmwfKX4eRfQ4+RzS6utKKKxz7ZyMheHRg/sPORNzDGmBDkS7PSBhF5CnjDnb4E2Og+Da4iYJG1hNI8mP0rSBsAY3/nlyJf/nY7+/LLmHrpcLvhzRjTavmSHCYDvwBuxXngz1fA7TiJ4dRABdYiPv4DFOyBaz+G6LhmF1dSXsUzX2xmTL80jj8q1Q8BGmNMcPhyKWsJ8A/3VVOh3yNqKdu/haUvwok3++XqJIDXFu3gQGG5jbpqjGn1fLmUtR/wF2AwcPjntaoeFcC4AquyHN69Fdr3hFP981C70ooqpn2xmROO6khm745+KdMYY4LFlw7p6cBTQCVOM9K/gZcDGVTAffv/IHs9nPUIxPjnGc7/XZrFvvwyfjXeag3GmNbPl+QQr6qfAKKq21X1fmB8YMMKoINb4PO/OeMmDTjDL0WWV3p4ev5mRvbqwIl9ra/BGNP6+dIhXSoiEcD3InIzsAuvYbVbFVV473aIiIYz/XeLxszlWezKLeGhC4bYFUrGmLDgS83hViAB+BUwErgCuCqAMQXO6rdg8ycw/vfQzj93Lld5lCfnb2ZoenvG9u/klzKNMSbYfLlaabH7thC4OrDhBFBJLnx4N3QbBqOv91uxH63Zy/acYp6+fITVGowxYaPe5CAisxvaUFXP8384h/d9C3AzTif4e6ra/AGPPvszFB+ASTMgIrLZxYHzvIZnvthC79QEThtc1yOxjTGmdWqo5nAisBN4HViIcwNcwInIqcBEYKiqlolI8/s39q2Bxc/ByKuh+/BmF1dtyfZDrNyZyx8nHkOkjaFkjAkjDSWHrsBpwGXAz4H3gNe9x1kKkJuAh1W1DEBV9zerNFV4/06Ia+f0NfjRs19soUNCNBeN7OnXco0xJtjq7ZBW1SpV/VBVrwJOADYB890mn0DqD4wRkYUi8rmIjGpWaWvehu1fwfj7IMF/N6dtyS7k43X7uOKEXsTH+KeZyhhjQkWDHdLu4Hpn49QeegNTgbebu1MRmYdTM6npXjemDjgJaRQwQ0SOUlWtUcYUYApARkZG3TsqL4K590HXoTBycnPD/oHnv9pKdGQEV5zY26/lGmNMKGioQ/olYAjwAfCAqq72105VdUID+70JeNtNBotExAOkAdk1ypgGTAPIzMzUWgUBfPlPyN8FF73gt05ogJzCMt5cmsWFw3vQKbn5z38wxphQ01DN4QqgCKeZ51del2kKoKraLkAxvYNzB/Z8EekPxAAHGl3KwS3wzVQYeglknODXAF9buIOySg/Xjenj13KNMSZU1JscVDVYDz5+AXhBRFYD5cBVNZuUfPLhPRAZAxMe8GtwFVUeXl24g1P6d+Lozsl+LdsYY0KFL8NntChVLQcub1Yhm+bBxg+cxNCum38Cc328dh9780v50/lD/FquMcaEkmDVDgLHU+V0QnfoAyfc5PfiX/pmG+kd4jnVHgFqjAlj4Zcclr8M+9fCaQ9AlH87i9fvzWfh1oNccUIvu+nNGBPWwis5lBXApw9BxokwyP+je7z87XZioyK4ONNuejPGhLeQ63Nolq8fg6L9cNkb4OdB8PJKKnh72S7OO647HRJj/Fq2McaEmvCpOeRlwTePw5CLIH2k34t/a2kWJRVVXHVSb7+XbYwxoSZ8ksMnfwT1wIT/83vRHo/y8oLtjMhIYUiP9n4v3xhjQk14JIeKYvjuDefqpJR6htJohgVbc9h6oIjLT+jl97KNMSYUhUdyyN8NCakw5rcBKf6NRTtpFxfFWcf6954JY4wJVeGRHMoKYNzdEOf/Jp9DReV8uHovFwzvQVy0jb5qjGkbwiM5xCQ5D/IJgJnLd1Fe5eHS0f5vrjLGmFAVHskhrR9E+v+qXFXljcU7OC69PYO6BWqcQWOMCT3hkRwCZPnOXDbuK7RagzGmzbHk0IA3Fu0gISaSc4/rHuxQjDGmRVlyqEdBaQVzVu7h3KHdSYoNrxvJjTHmSCw51GPOyj2UVFRxyWgbR8kY0/aEXHIQkWEiskBEVojIEhEZHYw4ZizZyYAuyQzvmRKM3RtjTFCFXHIA/obzzOphwB/c6Ra1ObuQFTtzuWhkOuLnAfyMMaY1CMXkoED1daPtgd0tHcDby7KIEJg4zDqijTFtUyj2tN4KfCQij+Akr5NacucejzJz2S7G9OtE53ZxLblrY4wJGUFJDiIyD+hax6J7gR8Dv1HVt0TkYuB5YEIdZUwBpgBkZPjvPoQFW3PYnVfKXWcO9FuZxhjT2gQlOahqrZN9NRH5N/Brd/K/wHP1lDENmAaQmZmp/ort7WW7SIqN4vTBdeUuY4xpG0Kxz2E3MNZ9Px74vqV2XFxeyQer9nD2sd2Ij7FB9owxbVco9jlcDzwmIlFAKW7TUUuYu2YfReVVXDiiR0vt0hhjQlLIJQdV/Qrw/3M+ffDWsizSO8QzqnfHYOzeGGNCRig2KwXF3rxSvtp0gAuH9yAiwu5tMMa0bZYcXO+s2IUqXDAiPdihGGNM0FlycM1asZthPVPok5YY7FCMMSboLDkAm/YXsG5PPufZ0NzGGANYcgBg9so9RAicM7RbsEMxxpiQ0OaTg6oyZ+VuTjgq1YbLMMYYV5tPDmt257P1QJE97c0YY7y0+eQwe+VuoiKEM4fYcBnGGFOtTScHj0d5d+VuTunfiZSEmGCHY4wxIaNNJ4elOw6xO6/UrlIyxpga2nRymLNyN7FREUwY3CXYoRhjTEhps8mhssrD+6v2MGFQF5JiQ26IKWOMCao2mxy+3ZLDgcJyu0rJGGPq0GaTw5yVu0mKjWLcgE7BDsUYY0JOm0wOFVUe5q7dx2mDuxAXbQ/1McaYmoKSHETkZyKyRkQ8IpJZY9ndIrJJRDaIyE8Csf8FW3LILa7gDLu3wRhj6hSsntjVwIXAM94zRWQwcClwDNAdmCci/VW1yp87/2D1XhJiIhnb35qUjDGmLkGpOajqOlXdUMeiicAbqlqmqluBTcBof+67yqPMXbOX8QM7W5OSMcbUI9T6HHoAO72ms9x5frN420EOFJZz5hAbgdUYY+oTsGYlEZkH1NWof6+qzqpvszrmaT3lTwGmAGRkZPgc1wer9hAXHWFXKRljTAMClhxUdUITNssCenpNpwO76yl/GjANIDMzs84EUpPHo3y4Zi9j+3ci0W58M8aYeoVas9Js4FIRiRWRPkA/YJG/Cl++8xD78ss461hrUjLGmIYE61LWC0QkCzgReE9EPgJQ1TXADGAt8CHwS39eqfT+qr3EREYwfmBnfxVpjDFhKShtK6o6E5hZz7KHgIcCsE8+XL2XMf3SSI6L9nfxxhgTVkKtWSlgvsvKY1duid34ZowxPmgzyeGD1XuJihBOs+G5jTHmiNpEclBVPli9hxP7ptoT34wxxgdtIjls3FfI9pxia1IyxhgftYnk8PHavQBMGGRNSsYY44s2khz2cVzPFLq0iwt2KMYY0yqEfXLYl1/Kyqw8TreOaGOM8VnYJ4eP1+4DsKuUjDGmEdpEcuiVmkC/zknBDsUYY1qNsE4OhWWVfLs5h9MHd0GkrgFfjTHG1CWsk8PnG7Ipr/Jw2mC7hNUYYxojrJPD3LV76ZgYw8heHYIdijHGtCphmxwqqjx8tn4/4wd2JjLCmpSMMaYxwjY5LNp6kPzSSrtKyRhjmiBsk8PHa/cRGxXBmH5pwQ7FGGNanbBMDqrKx2v3MaZfGgkx9jhQY4xprGA9Ce5nIrJGRDwikuk1/zQRWSoiq9y/45tS/to9+ezKLbEmJWOMaaJg/axeDVwIPFNj/gHgXFXdLSJDgI+AHo0tfN7a/YjA+IGWHIwxpimC9ZjQdUCtG9NUdbnX5BogTkRiVbWsMeV/un4fw3qm0Ck5ttmxGmNMWxTKDfI/BZbXlxhEZAowxZ0sFZE1tdb5ZZP22x7IC/B2R1q3oeWNXVbXvDScWlowNfV79mdZvm7ny3p2zFqmrGAfs8bMbw3HrFe9S1Q1IC9gHk7zUc3XRK915gOZdWx7DLAZ6Ovjvqb5Me4mldWY7Y60bkPLG7usnnlLAnXcA/09B+OY+bKeHbO2ccwaM7+1H7OA1RxUdUJTthORdGAmcKWqbvZxszlN2Zefy2rMdkdat6HljV3mz+/Gn1rTMfNlPTtmLVNWsI9ZY+cHW5PjEje7BIWIzAduV9Ul7nQK8DnwoKq+FbTAwpyILFHVzCOvaUKFHbPWp7Ufs2BdynqBiGQBJwLvichH7qKbgaOB+0RkhfvqHIwYw9y0YAdgGs2OWevTqo9ZUGsOxhhjQlNY3iFtjDGmeSw5GGOMqcWSgzHGmFosORhE5HwReVZEZonI6cGOxxyZiAwSkadF5E0RuSnY8ZgjE5FEd8y4c4Idiy8sOYQpEXlBRPaLyOoa888QkQ0isklEfgegqu+o6vXAZOCSIIRraPQxW6eqNwIXA632csnWrDHHy3UXMKNlo2w6Sw7h60XgDO8ZIhIJPAGcCQwGLhORwV6r/N5dboLjRRpxzETkPOAr4JOWDdO4XsTH4yUiE4C1wL6WDrKpLDmEKVX9AjhYY/ZoYJOqblHVcuANYKI4/gp8oKrLWjpW42jMMXPXn62qJwGTWjZSA40+XqcCJwA/B64XkZA/94bywHvG/3oAO72ms4DjgVuACUB7ETlaVZ8ORnCmTnUeMxEZhzPsfSzwfsuHZepR5/FS1ZsBRGQycEBVPUGIrVEsObQtUsc8VdWpwNSWDsb4pL5jNh9n4EoTWuo8XoffqL7YcqE0T8hXbYxfZQE9vabTgd1BisX4xo5Z6xI2x8uSQ9uyGOgnIn1EJAa4FJgd5JhMw+yYtS5hc7wsOYQpEXkd+BYYICJZInKtqlbiDG74EbAOmKGqtR6SZILDjlnrEu7HywbeM8YYU4vVHIwxxtRiycEYY0wtlhyMMcbUYsnBGGNMLZYcjDHG1GLJwRhjTC2WHExYEZEqEVnh9eod7Jj8QUQmi0i2iDznTo8TkXdrrPOiiFzUQBl/F5G9InJ7oOM1rZ+NrWTCTYmqDqtrgYgIzr09IT/oWT3+Uz2AW1Oo6h0iUuTPgEz4spqDCWsi0ltE1onIk8AyoKeI3CEii0XkOxF5wGvde92HtMwTkderf2GLyHwRyXTfp4nINvd9pPtrvLqsG9z549xt3hSR9SLyqpuYEJFRIvKNiKwUkUUikiwiX4rIMK84vhaRoc34zJleNadVImJ3uppGs5qDCTfxIrLCfb8V+A0wALhaVX/hPga1H864+wLMFpFTgCKccXCG4/y/WAYsPcK+rgXyVHWUiMQCX4vIXHfZcOAYnEHXvgZOFpFFwH+AS1R1sYi0A0qA53CewneriPQHYlX1Ox8+6xivzwqQAbyrqkuAYeA0JQEf+lCWMT9gycGEmx80K7l9DttVdYE763T3tdydTsJJFsnATFUtdrfzZbC004GhXu387d2yyoFFqprllrUC6A3kAXtUdTGAqua7y/8L3CcidwDX4DxhzBdfqurh5xGLyA+2E5GLgRFunMY0iiUH0xZ4t7ML8BdVfcZ7BRG5Fa9x92uo5H9NsHE1yrpFVT+qUdY4oMxrVhXO/zWpax+qWiwiH+M8Mcwvz4QWkWOAB4BTVLWqueWZtsf6HExb8xFwjYgkAYhIDxHpDHwBXCAi8SKSDJzrtc02YKT7/qIaZd0kItFuWf1FJLGBfa8HuovIKHf9ZBGp/oH2HM4Dlxaras1HTzaKiLTHeTzllaqa3ZyyTNtlNQfTpqjqXBEZBHzr9hEXAper6jIR+Q+wAtgOfOm12SPADBG5AvjUa/5zOM1Fy9wO52zg/Ab2XS4ilwD/T0TicfobJgCFqrpURPKB6X74mOcDvYBn3c9IfVdwGVMfG7LbmDqIyP04J+1HWmh/3XEe+zmwrktt3WcPZzbnUla3nPtpwc9lWi9rVjImyETkSmAhcG8D92CUAGdW3wTXxP38HbicH/bBGFMnqzkYY4ypxWoOxhhjarHkYIwxphZLDsYYY2qx5GCMMaYWSw7GGGNqseRgjDGmlv8PjCVy9JdmLPEAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import scipy\n",
    "\n",
    "target_sos = np.vstack( (np.array([stage1]), np.array([stage2])) )\n",
    "w, h = scipy.signal.sosfreqz(target_sos, worN=2**20, fs=samplerate)\n",
    "\n",
    "calculated_sos = np.vstack(\n",
    "    (np.array(hf_boost_shelv(K1, V_0)),\n",
    "     np.array(highpass(K2, Q)))\n",
    ")\n",
    "w2, h2 = scipy.signal.sosfreqz(calculated_sos, worN=2**20, fs=samplerate)\n",
    "\n",
    "fig, ax1 = plt.subplots()\n",
    "ax1.set_title('Digital filter frequency response')\n",
    "\n",
    "ax1.plot(w, 20 * np.log10(abs(h)), label='EBU')\n",
    "ax1.plot(w2, 20 * np.log10(abs(h2)), label='jkbd')\n",
    "\n",
    "ax1.set_xlabel('Frequency [Hz]')\n",
    "ax1.set_xscale('log')\n",
    "ax1.set_ylabel('Magnitude [dB]')\n",
    "\n",
    "plt.xlim(20, samplerate/2.0)\n",
    "plt.ylim(-12, 6)\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0ff1320e",
   "metadata": {},
   "source": [
    "Is there something wrong with highpass filter?\n",
    "\n",
    "Looks like the optimisation of the second filter is weak. The result depends heavily on the initial guess."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5510b831",
   "metadata": {},
   "source": [
    "## Second try"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "a114569c",
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.optimize import fsolve\n",
    "\n",
    "def solve_stage2(samplerate):\n",
    "    \"\"\"Minimize the difference between target biquad and calculated biquad.\"\"\"\n",
    "    def helper(z):\n",
    "        K = z[0]\n",
    "        Q = z[1]\n",
    "        sos = highpass(K, Q)\n",
    "        result = stage2[4:] - sos[4:]\n",
    "        return result\n",
    "    \n",
    "    guess = np.array([0.002, 0.25], dtype=np.float64)\n",
    "    result = fsolve(helper, guess)\n",
    "    return result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "28c61d50",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.00249597 0.50032704]\n",
      "freq = 38.135470876002174 Hz\n",
      "Q = 0.5003270373223665\n"
     ]
    }
   ],
   "source": [
    "result = solve_stage2(samplerate)\n",
    "print(result)\n",
    "\n",
    "K2, Q = result\n",
    "print(f'freq = {cutoff_freq(K2, samplerate)} Hz')\n",
    "print(f'Q = {Q}')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a2517242",
   "metadata": {},
   "source": [
    "Frequency response, again:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "b3f94c6d",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/tmp/ipykernel_28555/1576023010.py:13: RuntimeWarning: divide by zero encountered in log10\n",
      "  ax1.plot(w, 20 * np.log10(abs(h)), label='EBU')\n",
      "/tmp/ipykernel_28555/1576023010.py:14: RuntimeWarning: divide by zero encountered in log10\n",
      "  ax1.plot(w2, 20 * np.log10(abs(h2)), label='jkbd')\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYcAAAEaCAYAAAD65pvjAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAAzPElEQVR4nO3deXwV5bnA8d9zshMSlhCWEAIooIh4UYIbioi4b3WpiivqFW2L1V61brXVWlvvrbXWpSpu4FKtCoq7iIoLoAICsrkAIoSwBpKQPTnnuX/MQA/ZOAk5mbM838/nfHJme9/nzOTMc+admXdEVTHGGGOC+bwOwBhjTOSx5GCMMaYBSw7GGGMasORgjDGmAUsOxhhjGrDkYIwxpgFLDgYAEXlMRO5o63mbWF5FZEAT03qIyKciskNE/iYit4nIk+60fu6yia2tO4TYdqs/XPUYE+nC9iUzkUNE1gA9gDrADywHngUmqWoAQFWvCbW84HlFZDTwvKrmtlG4E4CtQKbu4SYcEZnl1v1kG9XdovqNiWV25BA/TlfVDKAvcC9wM/CUtyE1qi+wvD12zCKS0NL6w3nUEomaWEcmHqiqvWL8BawBxtYbdygQAA50hycDfwqa/ltgA1AI/DegwIDgeYF0oNItp8x95bhlzwWK3TIeBpKDyt5VVr2YJgO1QI1b1ljgTpyjA4B+7rKJwD04R0FV7rwPu/PsD3wAbAO+A86rV/6jwDtAeSPrpKn6XwWeB0rdddEJJ7FuANa76yLBLSMBuA/n6GM18KudMTe2LYI/nzt8ODDHXXeLgdFB02YBdwOzgR3ADKBb0PSjgpZdB4wHRgCbdtbvzncOsKiJ/5UG68jdplOBLcCPwK/r/R/Nd9fNJuD+ettqAs7/0AbghqDlUoAH3GmF7vsUd9pooAC4AdjsLnt50LKn4Bz97nDX/41B004DFrnrYA5wkNffv2h9eR6AvdphIzeSHNzxa4FfuO8n4yYH4CRgIzAE6AA8RyPJwX0/GiioV+5wdyeX6O4kVgDXB01vNDnUL9sdvpNGkoM7PAv476B5092d4uVu3Yfg7KSHBJVdAozEOWpODbH+WuBn7jJpwOvA42593YGvgKvd+a8BvgX6AF2BjwkxOQC9gSJ35+cDjneHs4M+7ypgkBvHLOBed1oezs5yHJAEZAHD3GnLgZOD6nyNoB11I58/eB11ABYAvweSgX1wkt6J7vxzgUvc9x2Bw+ttqxfd9TQUJ7mMdaf/EfjCXX/ZODvyu4P+p+rceZLc9VEBdHGnbwCOdt93AQ5x3x+Ck0wOw0nSl7nrO8Xr72A0vqxZKb4V4uzA6jsPeEZVl6lqBXBXSwpV1QWq+oWq1qnqGpwd6TF7He2enQasUdVn3Lq/xvnFe27QPNNVdbaqBlS1KsRy56rq6+qcn8kETsZJduWquhn4O3CBO+95wAOquk5VtwF/aUH8FwPvqOo7bnwf4PwqPyVonmdU9XtVrQReBoa54y8CZqrqi6paq6pFqrrInTbFLRsR6QqcCPyrmTh2rSOcnXq2qv5RVWtUdTXwRNDnrQUGiEg3VS1T1S/qlXWXu56WAM/gJK+d8f5RVTer6hac/7FLgpardafXquo7OEdy+wVNO0BEMlV1u7udAa4CHlfVL1XVr6pTgGqcHyqmhSw5xLfeOM0v9eXg/ALfaV0j8zRJRAaJyFsislFESoE/A91aH2bI+gKHiUjxzhfOTqhn0Dwt+iyNLNMX59fshqA6Hsf5BQwN191PLainL/DzevEfBfQKmmdj0PsKnF/r4ByprGqi3OeB00WkI07y+kxVNzQTR/3Pm1MvpttwLnAAuBLnSOZbEZknIqc1U9ZPOOsH9+9PTUwDKFLVuqDh4M96Dk7C/ElEPhGRI4JivaFerH3qlWtCFFcn18x/iMgInOTweSOTNwDBVx/1aaaoxk7cPgosBMap6g4RuZ7df723lfp1rwM+UdXjW7BMS+tZh/NrtFu9nddOG9h9feXVm16O01SzU/3E9ZyqXtWKGNfhtP83oKrrRWQucBbOr/NH91BW/c/7o6oObKLsH4BxIuIDzgZeFZGsoFn64DSzgbMuCt33hTg782WNTGs+ONV5wJkikgRMxDmC6uPGeo+q3hNKOaZ5duQQZ0Qk0/119xJOW/eSRmZ7GbhcRAaLSAec9uambAKyRKRT0LgMnBOUZSKyP/CLNgq/sbr3CRp+CxgkIpeISJL7GiEig9uqQvcX9wzgb+669InIviKys9nsZeDXIpIrIl2AW+oVsQi4wI0tn92T5s5f+CeKSIKIpIrIaBEJ5TLhF4CxInKeiCSKSJaIDAua/izORQZDcc45hOoroFREbhaRNDeuA90fF4jIxSKS7TZBFbvL+IOWv0NEOojIEJxzQf92x78I/E5EskWkG87/2PN7CkZEkkXkIhHppKq1OP9nO+t7ArhGRA4TR7qInCoiGS34vMZlySF+vCkiO3B+Xd0O3I/zZW1AVd8FHsQ5mboS56QjOL+Y68/7Lc4XfbV7KJ8D3AhciHOC9An+s0Noa/8AzhWR7SLyoKruAE7AaQ8vxGmC+V+cK2Pa0qU4J2eXA9txrmba2fTzBPA+zpVGXwPT6i17B7Cvu9xdBLX9q+o64EycZpstONvqJkL4nqrqWpymlhtwmgoXAf8VNMtrOL/UX1PV8lA/qKr6gdNxzm38iHOC/0mcK7bAuXhhmYiU4WyPC+qdy/kE53/oQ+A+VZ3hjv8TzvmUb4AlOOvqTyGGdQmwxm2yvAb3fIqqzsc57/AwzvpdiXPFlmkFUbX7fEzz3F/eS3Gu+misKcU0QUT64exUk7xedyKyCueqqpntUFc/IuRzm9axIwfTKBE5yz2E74Lz6/tN+5JHLxE5B+dcwkdex2KiQ8QmBxHpLCKvisi3IrIi6IoE0z6uxmnaWIXTphuu8wYmzNxuRh4FfuWeGzBmjyK2WUlEpuBccvekiCQDHVS12OOwjDEmLkRkchCRTJwTevtoJAZojDExLlKblfbBadJ4RkQWisiTIpLudVDGGBMvIvXIIR+n35WRqvqliPwDKFXVO4LmmYDTqRfp6enD999/f2+CNcaYKLVgwYKtqprd2LRITQ49gS9UtZ87fDRwi6qe2tj8+fn5On/+/HaM0Bhjop+ILFDV/MamRWSzkqpuBNaJyM6Oto7DueHIGGNMO4jkvpWuBV5wr1RaTRN38xpjjGl7EZsc3O6GGz3cMcYYE14Rmxz2Vm1tLQUFBVRVhdplf/RJTU0lNzeXpKQkr0MxxsSYmE0OBQUFZGRk0K9fP0TE63DanKpSVFREQUEB/fv39zocY0yMicgT0m2hqqqKrKysmEwMACJCVlZWTB8ZGWO8E7PJAYjZxLBTrH8+Y4x3Yjo5eC0hIYFhw4btet17770AjB49mv32249hw4YxePBgJk2atGuZjh077lbG5MmTmThxYrvGbYwxMXvOIRKkpaWxaNGiRqe98MIL5Ofns23bNvbdd1/Gjx9PcnJy+wZojDFNsCMHj5WVlZGenk5CQoLXoRhjzC5xceRw15vLWF5Y2qZlHpCTyR9OH9LsPJWVlQwbNmzX8K233sr5558PwEUXXURKSgo//PADDzzwgCUHY0xEiYvk4JVQmpW2bNnCkUceyUknnUTfvn0bnddOPBtj2ltcJIc9/cL3UnZ2Nocccghffvklffv2JS0tjZqaml3nH7Zt20a3bt08jtIYE2/snIPHKioqWLhwIfvuuy8AxxxzDM8//zzgNEu9/PLLHHvssV6GaIyJQ3Fx5OCV+uccTjrppF2Xs1500UWkpaVRXV3N+PHjGT58OAD/+Mc/uPrqq3nwwQdRVS699FJGjRrlRfjGmDhmySGM/H5/o+NnzZrV5DK9e/fmrbfeClNExhgTGmtWMsYY04AlB2OMMQ1EbHIQkQQRWSgi1sZijDHtLGKTA3AdsMLrIIwxJh5F5AlpEckFTgXuAf7H43CMMRFMVamuC1BeXUdFjZ+KGj/VdX5q/UqtP0Cd+9d5KXWBADV1AeoCSl1AUVUCASWgoG55AXWHFQLqztPYMM6yTcXVtEBzH6iJ8c0t0ory9iAikwPwAPBbIMPjOIwx7ay6zs/GkirWF1eyfnslW8qq2VFWTnXpVvxlW9GKInxV2/HVlpFYV0mSv4JUqkmnijSqSZcqUqglET9J1JGEn1SpI5E6kvGTSJ0zXvwkEMBHAJ+TFvC5Lwn6K7vGB+oNKz5p3Y43UtzWzLSISw4ichqwWVUXiMjoZuabAEwAyMvLa5/gWuHII4/kz3/+M/fdd1+DS1TvvPNOOnbsyI033thsGePHj+e0007j3HPPDWeoxrSrkspalheW8sPmHaxdv56aDStILF5NZvUGcigiR7YyXIrIlhIypLLxQgRIhAA+6hI74E9II5DovEhIAl8SmpCG+JIgMRkSkpCEJCQhGV9iMvh8iC8RQcDnA/EhIoj4nGkizjifD8EXNN4H7jQQ532TATY1qblucZqY1mxPOq2o666m9z0RlxyAkcAZInIKkApkisjzqnpx8EyqOgmYBJCfnx+x6XvOnDnN3tdgTLzYWFLFpz9s4fuVq6j96Ut67VjCgbKaE33r6SHFu+bTRKEqJZu6jN74Oo8gpUsOdMyGDl2hQ5bzN60rpGRAckdI7oAvMZVk64OsFaIoOajqrcCtAO6Rw431E0M06dix425HDPPmzWPChAlMnToVgMWLFzNmzBjWrVvHb3/7W6666ipUlWuvvZaPPvqI/v3776Ht0pjIpKosKyxlxqIfKV7xEQOK5zDK9w3n+TYB4E9KpKLL/iT2PBHNGYJ0HwzdBiCZuaQl2rNNvBZxySEs3r0FNi5p2zJ7DoWT723RInPmzOHaa69l+vTpu5rCvvnmG7744gvKy8s5+OCDOfXUU/niiy/47rvvWLJkCZs2beKAAw7giiuuaNv4jQmTTaVVTJu3ms3zp3NY2Yf8yreIFKmlNjmVqj5HERg4EV/eYST0+i8yklK9Dtc0IaKTg6rOAmZ5HEabWLFiBRMmTGDGjBnk5OTsGn/mmWeSlpZGWloaxx57LF999RWffvop48aNIyEhgZycHMaMGeNh5MaEZklBCa9/9Bm5PzzHON9ndJZyKtOy4MDLYMipJPUdSVJiitdhmhBFdHJoMy38hR8OvXr1oqqqioULF+6WHOo/q2HnsD3DwUSLbwqKmfrGdI7c8Cy3+xagCT6qBp0Oh15CWv/RkBAfu5lYE8k3wcWUzp078/bbb3PbbbftdoJ6+vTpVFVVUVRUxKxZsxgxYgSjRo3ipZdewu/3s2HDBj7++GPvAjemCT8VlXP3U6+y4fFzuGvzrxmVupLakb8h4X+WkX7hFBgw1hJDFLMtF2bBRwA9evTgzTff5OSTT+bpp58G4NBDD+XUU09l7dq13HHHHeTk5HDWWWfx0UcfMXToUAYNGsQxxxzjVfjGNFBd52fyh4vJmP1nbpeZ1CZ3oPrIW0k76lfOFUQmJlhyCKOioiK6du3K6NGjGT16NODck7Fs2TIADjvssEaXExEefvjh9grTmJAtWrudt//1IBMqnyTLV0blsCtJP+F25/JSE1MsOYRJYWEho0eP3uMNbsZEA39AeXrGfPrOuYXbffMp7XYQvnMfJr3Xf3kdmgkTSw5hkpOTw/fff+91GMbstU2lVUx65ikmbPs/shLKqRx9F5lHXwu+BK9DM2FkycEY06QFa4qY/+yt3BH4N6UZ+5B40Zsk9jrI67BMO4jp5KCqMX1JqN05bcJp6twVZL57LVf75lEy6Fw6nfsQJHfwOizTTmI2OaSmplJUVERWVlZMJghVpaioiNRUu8PUtK1AQHls+keMXXgt+/o2UDHmT3Q6euIeOoozsSZmk0Nubi4FBQVs2bLF61DCJjU1ldzcXK/DMDGk1h/gwedf5dLVN5KeFIBx0+gw4FivwzIeiNnkkJSURP/+/b0Ow5ioUV5dx4NPPcXETX+A1EzSrnwD6b6/12EZj8RscjDGhK6kopbHH72PG0r/SkVmfzpf9SZk5ux5QROzLDkYE+eKK2p4+p9/5sYdD1CSPZwuV74KaV28Dst4zJKDMXFse3kNkx/5E9eXP0hxzyPoeuVUuyLJAJYcjIlb28prmPLI3VxX/hDFvUbS9cpXISnN67BMhLDkYEwcKiqr5vlH7uI3lY+wLWcUXa94BezBOyZIRHbZLSJ9RORjEVkhIstE5DqvYzImVmwtq+ZfD/+B6yofYVvOaEsMplGReuRQB9ygql+LSAawQEQ+UNXlXgdmTDTbWlbNiw/fwbVVj7Ot93F0vfxFsKezmUZE5JGDqm5Q1a/d9zuAFUBvb6MyJrrtlhhyj6fr5S9ZYjBNisjkEExE+gEHA1/WGz9BROaLyPxYvgvamLbQIDGM/xckJnsdlolgEZ0cRKQjMBW4XlVLg6ep6iRVzVfV/OzsbG8CNCYKWGIwrRGxyUFEknASwwuqOs3reIyJRk5i+L2bGMZaYjAhi8jkIE43qk8BK1T1fq/jMSYa/ScxPOYmhhctMZiQRWRyAEYClwBjRGSR+zrF66CMiRaFxZW8+NDtbmI4zhKDabGIvJRVVT8HrPN4Y1rhxy1lfPD4TVxb9y+2551A10tfsMRgWiwik4MxpnVWFJYw/4mJTNA32D7wHLpcMAkS7GtuWs7+a4yJEXN+2MjGF37BJXxE8dAr6HLW38AXqS3HJtJZcjAmBkybu4Ksd6/mbN9iSkZcT+dT7rTHepq9YsnBmCgWCCiT3vqEY+ZPZJBvPRUn3k+nI670OiwTAyw5GBOltpfX8M9nn+eqjXeRkViHjptKh4FjvA7LxAhLDsZEoYU/bWP2c3dyc+1zVKTnkjr+ZaT7YK/DMjHEkoMxUaTWH2DyBwvoO+cWJvrmU9z/ZDpf8DikdvI6NBNjLDkYEyW+27iDV194jAmlD9LFV0HlsXfTedS1duLZhIUlB2MiXGlVLU+/M5f+C+/l9oTZlHYeTOK4J0nseaDXoZkYZsnBmAhVUxdg2ler2TLzfv7bP5XURD8Vh99I5thbICHJ6/BMjLPkYEyEqa7zM23ej/z04VNcVPMKfXxbKO17PIln/i+JWft6HZ6JE5YcjIkQ64srmTp7OVULXuBC/3TGyVZKux2EnvI4mQOO8zo8E2csORjjoZLKWj5YtpEl8z9hv/VTudI3m3SppjT7YPTEx8gcMNZOOBtPWHIwph2pKis3l/H5D5v5aelceq5/n5PlC871baY2MYWawWfByAlk9h7udagmzllyMCaMdlTVsqywlCXriin8cRmJ6+ZwUM1CzvAtI0t2EEhIoKz3SPSQ35F0wBkkpXXxOmRjgAhODiJyEvAPIAF4UlXv9TgkYxpVUxdg844q1m6rYPWWcgo2baVq43f4tq0ku/x7hspqzvf9SKZUAFDZIZtA/xNh/7H4Bp5AZnqWx5/AmIYiMjmISALwCHA8UADME5E3VHW5t5GZeKCqlFbVUVJRS3FlDcUVtRRX1lJcUcP2shp2FG+idvt6KC0kqWIjHWu20IPt9JVNHOfbSC/Ztqssf1Ii5V32J6nPuZCXD30OIy17PzuPYCJeRCYH4FBgpaquBhCRl4AzAUsOJiSBgLKjuo7SylpKKmspraqltLKW0so6531FDVXlJdRWFBMo34ZWbsNXVUxidQkptcVkUkZnyukiO+gk5exHGV2kjE6UkSJ1u9WliUJVcldqM/qQ2P04Aj33x9dtAGQNIKHbQDITUzxaC8a0XpPJQUQeDGH5UlX9XRvGs1NvYF3QcAFwWPAMIjIBmACQl5cXhhBMJNlRVcvWshq2lVe7f51XSUkxtaWb8ZdthYoiEqpLSKwtIaWujAzK6UQ5mVJBJuX0lnIGU7FrOEG08coSoc6XSk1yJ/wpnSG1M5Len4T0LBI7doXMXpCZ47wyeiEZPUlLSCKtXdeIMeHV3JHDmcDv97D8LUA4kkNjx9y7fZNVdRIwCSA/P7+Jb7mJBrX+AAXbK1lTVM6G4io2FpdTVrSBuhKn6Sa1YgNd/VvpJiV0pZSeUsoQKaUrO+gg1Y0X6u7ga5My8adkEkjphKT1xJfWiaT0LvjSu0BaZ6fDutTO0KErpHXZ9UpMSovYw2pj2kNz//9/V9UpzS0sIuG6tKIA6BM0nAsUhqku0062ldfw7cZSvt+4g8JNm6navBLf9tV0LF9HXzbQ17eJUVJED7aTJP7/LCjgT0qiOjULf2oWpOeRkNGN5Mzu0DEb0rMhvRt0yHJ27qmdILUTiYkptoM3ppWa/O6o6gN7WjiUeVppHjBQRPoD64ELgAvDVJdpY6pKwfZKFq4r5tt1m9lRsIzELcvoXb2a/WUtp/oKyJaS/yyQCJWp3anr1JeErkNJ6NoHOvWGzN5u801vEjpk0cGeh2xMu2nunEMqcD6wHXgT+C1wNLAKuFtVt4YrKFWtE5GJwPs4l7I+rarLwlWf2Tu1/gDfFBQz78ciCld+Q+KG+QyoXsEhvpWcKgW72vbrklOp7DKIxF6noD32Q7L2ha77QNf+pCWne/wpjDHBmjvqfhaoBdKBG4ClwMPAUcBk4LRwBqaq7wDvhLMO0zqqyqot5Xz+/WZWr1hAh4LPGBH4hnG+7+jkXstfnZpJTa/h0O8C6DUUug8hsWt/MnwJHkdvjAlFc8nhAFU9UEQSgQJVPcYd/56ILG6H2EwECQSUxQXFfPjNanYseZeDKuZysm8pPaQYBMo69SNxn3Og/xGQO4KUrAGkWDOQMVGrueRQA7uaeOqfDPY3Mr+JQcsLS3nry+VULX2LI2rmcK3vG1KklqrUzvj7jYbBx8E+x9Kxc589lmWMiR7NJYdc914HCXqPO9w77JEZz2wrr2H6gp/48au3OLTkPa7zLSBFaqno2BMOuAKGnkFq3hFgTUTGxKzmksNNQe/n15tWf9jEgO827mDqrK/osvw5zpaP6SHFVKV2RoeOh+Hj6JBziHX7YEycaO5S1mbvcTCxQVWZvbKI9z94i/wN/+Ym35ck+AJU9D0ODh9P6sATITHZ6zCNMe2suUtZ36TeXcnBVPWMsERk2sWupPDOVE7eOpm7E5ZTnZyOf9hVJI28ho5d+3sdojHGQ801K93n/j0b6Ak87w6PA9aEMSYTZgvXbmfa669w8pZnuDthORUdulE36h5S8i+DlAyvwzPGRIDmmpU+ARCRu1V1VNCkN0Xk07BHZtrcxpIqnnpjJvnf/527E+ZTmZZF3TH30OHQKyHJuo0zxvxHKF3PZIvIPkHdZ/cHssMblmlLtf4AUz5egnz6V26SdyApmeqjbydt5ERI7uB1eMaYCBRKcvgNMEtEVrvD/XC7yjaRb+n6El598Smu3vEQvXzbKB98Pumn/BEyenodmjEmgu0xOajqeyIyENjfHfWtqjbRT7KJFLX+AJPe/YrcL//InQmz2dFpIJz3Cum5+V6HZoyJAs1drXSIqn4N4CaDBl1mBM9jIse6bRU8PeUJfll8H10SKqgceTMZx95ol6QaY0LW3JHDMyIymsYfvLPTU8DBbRmQ2TtvLlxD0eu38wd5i9JOg0i8aAqJPQ7wOixjTJRpLjl0AhbQfHLY0rbhmNaq8wd4+LWPOPabmzjdt5odQ8eTeca9dhWSMaZVmruUtV87xmH2QnFFDQ89M4Vfbr6L9CSl7uznyBhi9ygaY1ov4p6iKCJ/BU7H6RV2FXC5qhZ7GlQEW7WljNee+BO3VE+ismMeqZe/At0Geh2WMSbKRWKH+x8AB6rqQcD3wK0exxOxFq3dxuxHruHGmkepyD2KzImzLDEYY9pExCUHVZ2hqnXu4BdArpfxRKrPvi1k7VOXcilvUjr0cjpd+RqkdfY6LGNMjNhjchDHxSLye3c4T0QODX9oAFwBvNtEXBNEZL6IzN+yJb7Oi7/39Sr8/xrHGfIZZUfeTObZf7dnKxhj2lQoRw7/BI7A6XAPYAfwyN5UKiIzRWRpI68zg+a5HagDXmisDFWdpKr5qpqfnR0/vXm89/Uqurx+IaN831B50t/oeMJt9owFY0ybC+WE9GGqeoiILARQ1e0isld3U6nq2Oami8hlwGnAcaraZLfh8WbGotV0fv1i8n3fU3Pm46QdfJ7XIRljYlQoyaFWRBJwn+0gItlAIFwBichJwM3AMapaEa56os3MxT/ScdrFjPB9S+0Zj5JqicEYE0ahNCs9CLwGdBeRe4DPgT+HMaaHgQzgAxFZJCKPhbGuqDBv1SaSp17G4b7l1Jz2MKmHXOB1SMaYGBdKx3sviMgC4Dicu6V/pqorwhWQqg4IV9nR6PuNpWx47r85w7eY8hPvJz3/Iq9DMsbEgeY63usaNLgZeDF4mqpuC2dgBjaUVDL3ieu4jE8pOewmOh1xpdchGWPiRHNHDgtwzjMIkAdsd993BtYC9pDhMKqoqWPqY3cx0T+NbftfSNeTbvc6JGNMHGnynIOq9lfVfYD3gdNVtZuqZuFcRTStvQKMR6rKE1Oe4ZqKx9naewxdf/6QXa5qjGlXoZyQHqGq7+wcUNV3gWPCF5J59u1ZXFbwB0o79qfbpc9CQsR1gWWMiXGh7HW2isjvgOdxmpkuBorCGlUc+3DRSg7/aiJJiT46XPEqpGR4HZIxJg6FcuQwDsjGuZz1daA7/7lb2rShtVvLCbz2Cwb4Ckm6YAqStY/XIRlj4lQol7JuA65rh1jiWk1dgHefuYur5SuKj/4DnQcd53VIxpg4tsfkICIf494dHUxVx4Qlojg1Zep0xpc9xeac0XQf8xuvwzHGxLlQzjncGPQ+FTgHp0M800Y+XryKsctupiqlC90vftquTDLGeC6UZqUF9UbNFpFPwhRP3Nm8o4rK135Nnm8LdRe8CelZXodkjDEhNSsF3yntA4YDPcMWURxRVV579iGu5nOKRtxA1r5HeR2SMcYAoTUrBd8pXQf8CFg/Dm3g7TmLOG/zA2zudCDdT7rN63CMMWaXUJLDYFWtCh4hIilhiidubCiuIP2DG0mXahIuespudDPGRJRQ7nOY08i4uW0dSDxRVd549u8cy3zKjrqNhB77ex2SMcbsprleWXsCvYE0ETkYp1kJIBPo0A6xxax35y5iXNHDbOpyMD3G2C0kxpjI01xbxonAeCAXuD9o/A4g7A3kInIj8FcgW1W3hru+9lJcUUPCB7eRKrV0vOhJ8CV4HZIxxjTQZHJQ1SnAFBE5R1WntmNMiEgf4HicrsFjytSXp3ClzmFL/g1kZ9tzjYwxkam5ZqWLVfV5oJ+I/E/96ap6fyOLtZW/A78Fpoexjna3YGUhY1f/H0Ud+pJ90s1eh2OMMU1qrlkp3f3bsT0C2UlEzgDWq+piaeZOYRGZAEwAyMvLa6foWq/WH+D7V37PcN9mKs+eDol2wZcxJnI116z0uPv3rrauVERm0viNdLfjnM84YU9lqOokYBJAfn5+g76fIs30mbM4p2oa6/udRe9Bo70OxxhjmhXKHdLZwFVAv+D5VfWK1laqqmObqGsozuNHdx415AJfi8ihqrqxtfV5bVt5Dd3n/pE6Xwo5P/8/r8Mxxpg9CuXOq+nAZ8BMwB/OYFR1Cc7zIgAQkTVAfrRfrfTWtGe5lIVsOfwOOnTsvucFjDHGY6Ekhw6qamdPW+n7wm0csfJ+ilJzyT7u116HY4wxIQnlDum3ROSUsEfSCFXtF81HDarK3JfvY6CsJ/mUv0BistchGWNMSEJJDtfhJIhKESkVkR0iUhruwGLB50t+4IztU1jf9TAyDjrd63CMMSZkoTzPwZ5w3wqBgLL57T+TKRV0POc+e4CPMSaqhHK10iGNjC4BflJVeyJcEz748mtOq3qL9XlnkNf7IK/DMcaYFgnlhPQ/gUOAJe7wUGAxkCUi16jqjHAFF61q/QGqP7wXnyi5Z7X5bSLGGBN2oZxzWAMcrKrDVXU4MAxYCowF7KL9Rrz3yWxOqZ3JxgEX4Ovaz+twjDGmxUJJDvur6rKdA6q6HCdZrA5fWNGrqtZPyud/oU6SyD3z916HY4wxrRJKs9J3IvIo8JI7fD7wvfs0uNqwRRal3p4xg3MCs1k/9Jf0zujhdTjGGNMqoRw5jAdWAtcDvwFWu+NqgWPDFFdUqqzx023+fZRLR3qfYvcNGmOiVyiXslYCf3Nf9ZW1eURR7L0PZ3CWzmf9wb8hPa2z1+EYY0yrhXIp60DgL8ABQOrO8aq6TxjjijpVtX46ffUAFdKB3idc73U4xhizV0JpVnoGeBSow2lGehZ4LpxBRaMZs2YxRr+gaMjlYEcNxpgoF0pySFPVDwFR1Z9U9U5gTHjDii41dQFS5z5ApaSSe/INXodjjDF7LZSrlapExAf8ICITgfUEdatt4IPPZ3OS/zPWD76SvPQsr8Mxxpi9FsqRw/VAB+DXwHDgEuCyMMYUVfwBxff5/dRJMn1O/a3X4RhjTJsI5Wqlee7bMuDy8IYTfT6Zt5CxtZ+wfuCF9LP7GowxMaLJ5CAibzS3oKqe0fbh7Kr7WmAizknwt1U1In+Sqyolsx7CJ2pHDcaYmNLckcMRwDrgReBLoF36nBaRY4EzgYNUtVpEIvb8xsIffmJsxbus7XUi/bvkeR2OMca0meaSQ0/geGAccCHwNvBicD9LYfIL4F5VrQZQ1c1hrq/VVr33Tw6RSpJOvsnrUIwxpk01eUJaVf2q+p6qXgYcjtOFxiy3ySecBgFHi8iXIvKJiIwIc32tsnrjdkYWvcK6zOGk9h3udTjGGNOmmj0h7XaudyrO0UM/4EFg2t5WKiIzcY5M6rvdjakLTkIaAbwsIvuoqtYrYwIwASAvr/2bdOa/8zTnyTZKjnuo3es2xphwa+6E9BTgQOBd4C5VXdpWlarq2Gbq/QUwzU0GX4lIAOgGbKlXxiRgEkB+fr42KCiMinZUMeSnZ9mU2pceQ09pz6qNMaZdNHefwyU4TTzXAXNEpNR97RCR0jDG9DruHdgiMghIBraGsb4W+3TGNIbIGuSIX4EvlFtFjDEmujR55KCqXu31ngaeFpGlQA1wWf0mJS/V+gN0XfoMpb5OdB9p9wIaY2JTKN1ntCtVrQEu9jqOpnw+bwGjAvNYe8DVZCal7nkBY4yJQtYm0kIln00CEfJOmOh1KMYYEzaWHFrgu4LNjCp7h7XdjiGhSx+vwzHGmLCx5NACS96fTFcpo9uYcN/qYYwx3rLkEKKSihr2X/sim1L6kjHYHmdhjIltlhxC9OnH73GgrMaffxVIu3QzZYwxnrHkEIJAQElZ+BQVkkbOqPFeh2OMMWFnySEE81f8wDG1n7Oh31mQkuF1OMYYE3aWHEJQMGsyKVJH7thfeh2KMca0C0sOe7C9rJqDNr/Oug5DSOk91OtwjDGmXVhy2IPZs95hgKwnYcR4r0Mxxph2Y8mhGapK0uLnqJQ0co680OtwjDGm3VhyaMbiVWs5uuZz1ueeAikdvQ7HGGPajSWHZqz+aAodpJqcMdd4HYoxxrQrSw5N2FFVy37rp1GYui8d+kXkk0qNMSZsLDk04fNPP2SI/Ih/2KV2R7QxJu5EXHIQkWEi8oWILBKR+SJyqBdx6NfPUk0yuaPsgT7GmPgTcckB+D+cZ1YPA37vDrerVRuLGFn5Met6HId06NLe1RtjjOciMTkokOm+7wQUtncASz/6N52kgm4jx7d31cYYExEi7jGhwPXA+yJyH07yOrI9Kw8ElK4rp7E9oStdDjy+Pas2xpiI4UlyEJGZQM9GJt0OHAf8RlWnish5wFPA2EbKmABMAMjLy2uz2Oav+J7D/V+zZuBldPEltFm5xhgTTTxJDqraYGe/k4g8C1znDr4CPNlEGZOASQD5+fnaVrEVfvY8h4qfPqOvaKsijTEm6kTiOYdC4Bj3/Rjgh/aquKKmjgEb3mJ96kBSc62TPWNM/IrEcw5XAf8QkUSgCrfpqD3Mnfs5x8lq1gy9o72qNMaYiBRxyUFVPweGe1F3xfx/UYePvFGXelG9McZEjEhsVvLExu3ljCidwU9djsSX0d3rcIwxxlOWHFzzZ71OT9lOx8Mu8ToUY4zxnCUHV+qKVyiXdHoM/5nXoRhjjOcsOQCrCjdzWPUXrO91PCSleh2OMcZ4zpIDsPyTaWRIJdlH2NPejDEGLDmgqmSsnE6JrzNdDjjO63CMMSYixH1yWLFmPYfXzWNTn5MhIeKu7DXGGE/EfXL44bOXSZVaeo28yOtQjDEmYsR1cggElOw1b7I1oTsZA0Z6HY4xxkSMuE4Oi75fzQj/Yrb3Pw18cb0qjDFmN3G9R1w3+yWSxE/vo+3GN2OMCRa3yaHOHyCn4B02JvWhQ97BXodjjDERJW6Tw4KlyxkeWEbZgDNBxOtwjDEmosRtctg090V8ovQZdbHXoRhjTMSJy+RQ6w+Qt/ED1qcMIKXXYK/DMcaYiONJchCRn4vIMhEJiEh+vWm3ishKEflORE4MR/1fL13GML6jYuBp4SjeGGOinle3BC8FzgYeDx4pIgcAFwBDgBxgpogMUlV/W1a++ctXAcgbeUFbFmuMMTHDkyMHVV2hqt81MulM4CVVrVbVH4GVwKFtWbc/oOQUzqAwuZ81KRljTBMi7ZxDb2Bd0HCBO67NLFrxHQfrCsr2OaUtizXGmJgStmYlEZkJ9Gxk0u2qOr2pxRoZp02UPwGYAJCXlxdyXIVfvMpwUXJHjgt5GWOMiTdhSw6qOrYVixUAfYKGc4HCJsqfBEwCyM/PbzSB1BcIKN0L3mdTUm965A5tRXjGGBMfIq1Z6Q3gAhFJEZH+wEDgq7Yq/JsfVjM8sJSS/qfYjW/GGNMMry5lPUtECoAjgLdF5H0AVV0GvAwsB94DftWWVyqtnfsqiRKg95F2lZIxxjTHk0tZVfU14LUmpt0D3BOGOum29n22JPYku+/wti7eGGNiSqQ1K4XN0lVryfcvYlveSdakZIwxexA3yWHN3Gkki5+cI8/3OhRjjIl4cZEcVJXOa96hKKEbGfsc7nU4xhgT8eIiOfxQsJkRdQvZknuCPfHNGGNCEBd7ypVzXydVauk+4hyvQzHGmKgQF8khddX77JCOdB082utQjDEmKsR8cthUXMbBVV9SmD0KErzqhNYYY6JLzCeHxXPep4uUkTHsDK9DMcaYqBHzySGw4m1qSKTXIad6HYoxxkSNmE4OZVW1HFD6GWs7jUBSM70OxxhjokZMJ4cF8+aQJ5tJGGyPAzXGmJaI6eSwY7Hz2Ii8I872OBJjjIkuMZscav0B+m6dxdq0wSR0yvE6HGOMiSoxmxwWLVvBUFZRte9JXodijDFRJ2aTw6Z5To/gfY441+NIjDEm+sRkclBVstbPZFNib9JyhngdjjHGRB2vngT3cxFZJiIBEckPGn+8iCwQkSXu3zGtKf/bn9Yz3P8N2/uMtWc3GGNMK3h15LAUOBv4tN74rcDpqjoUuAx4rjWF/zj3DZLFT89DraM9Y4xpDa8eE7oCQOr9qlfVhUGDy4BUEUlR1eqWlJ+2Zgalkknn/Y7a61iNMSYeRXJPdOcAC5tKDCIyAZjgDlaJyLIGM93Zqo/XCSgJ83J7mre56S2d1ti4bjhHaV5q7Xpuy7JCXS6U+WybtU9ZXm+zloyPhm3Wt8kpqhqWFzATp/mo/uvMoHlmAfmNLDsEWAXsG2Jdk9ow7laV1ZLl9jRvc9NbOq2JcfPDtd3DvZ692GahzGfbLD62WUvGR/s2C9uRg6qObc1yIpILvAZcqqqrQlzszdbU1cZltWS5Pc3b3PSWTmvLddOWommbhTKfbbP2KcvrbdbS8V5rdVziZhdPiMgs4EZVne8OdwY+Af6oqlM9CyzGich8Vc3f85wmUtg2iz7Rvs28upT1LBEpAI4A3haR991JE4EBwB0issh9dfcixhg3yesATIvZNos+Ub3NPD1yMMYYE5li8g5pY4wxe8eSgzHGmAYsORhjjGnAkoNBRH4mIk+IyHQROcHreMyeichgEXlMRF4VkV94HY/ZMxFJd/uMi4pHU1pyiFEi8rSIbBaRpfXGnyQi34nIShG5BUBVX1fVq4DxwPkehGto8TZboarXAOcBUXu5ZDRryfZy3Qy83L5Rtp4lh9g1GdjtSUcikgA8ApwMHACME5EDgmb5nTvdeGMyLdhmInIG8DnwYfuGaVyTCXF7ichYYDmwqb2DbC1LDjFKVT8FttUbfSiwUlVXq2oN8BJwpjj+F3hXVb9u71iNoyXbzJ3/DVU9EriofSM10OLtdSxwOHAhcJWIRPy+N5I73jNtrzewLmi4ADgMuBYYC3QSkQGq+pgXwZlGNbrNRGQ0Trf3KcA77R+WaUKj20tVJwKIyHhgq6oGPIitRSw5xJfGnnykqvog8GB7B2NC0tQ2m4XTcaWJLI1ur11vVCe3Xyh7J+IPbUybKgD6BA3nAoUexWJCY9ssusTM9rLkEF/mAQNFpL+IJAMXAG94HJNpnm2z6BIz28uSQ4wSkReBucB+IlIgIleqah1O54bvAyuAl1W14UOSjCdsm0WXWN9e1vGeMcaYBuzIwRhjTAOWHIwxxjRgycEYY0wDlhyMMcY0YMnBGGNMA5YcjDHGNGDJwcQUEfGLyKKgVz+vY2oLIjJeRLaIyJPu8GgReavePJNF5NxmyviriGwUkRvDHa+Jfta3kok1lao6rLEJIiI49/ZEfKdnTfj3zg7cWkNVbxKR8rYMyMQuO3IwMU1E+onIChH5J/A10EdEbhKReSLyjYjcFTTv7e5DWmaKyIs7f2GLyCwRyXffdxORNe77BPfX+M6yrnbHj3aXeVVEvhWRF9zEhIiMEJE5IrJYRL4SkQwR+UxEhgXFMVtEDtqLz5wfdOS0RETsTlfTYnbkYGJNmogsct//CPwG2A+4XFV/6T4GdSBOv/sCvCEio4BynH5wDsb5XnwNLNhDXVcCJao6QkRSgNkiMsOddjAwBKfTtdnASBH5Cvg3cL6qzhORTKASeBLnKXzXi8ggIEVVvwnhsx4d9FkB8oC3VHU+MAycpiTgvRDKMmY3lhxMrNmtWck95/CTqn7hjjrBfS10hzviJIsM4DVVrXCXC6WztBOAg4La+Tu5ZdUAX6lqgVvWIqAfUAJsUNV5AKpa6k5/BbhDRG4CrsB5wlgoPlPVXc8jFpHdlhOR84BD3DiNaRFLDiYeBLezC/AXVX08eAYRuZ6gfvfrqeM/TbCp9cq6VlXfr1fWaKA6aJQf57smjdWhqhUi8gHOE8Pa5JnQIjIEuAsYpar+vS3PxB8752DizfvAFSLSEUBEeotId+BT4CwRSRORDOD0oGXWAMPd9+fWK+sXIpLkljVIRNKbqftbIEdERrjzZ4jIzh9oT+I8cGmeqtZ/9GSLiEgnnMdTXqqqW/amLBO/7MjBxBVVnSEig4G57jniMuBiVf1aRP4NLAJ+Aj4LWuw+4GURuQT4KGj8kzjNRV+7J5y3AD9rpu4aETkfeEhE0nDON4wFylR1gYiUAs+0wcf8GdAXeML9jDR1BZcxTbEuu41phIjcibPTvq+d6svBeezn/o1daus+ezh/by5ldcu5k3b8XCZ6WbOSMR4TkUuBL4Hbm7kHoxI4eedNcK2s56/Axex+DsaYRtmRgzHGmAbsyMEYY0wDlhyMMcY0YMnBGGNMA5YcjDHGNGDJwRhjTAOWHIwxxjTw/9BP8+yoHTgcAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "target_sos = np.vstack( (np.array([stage1]), np.array([stage2])) )\n",
    "w, h = scipy.signal.sosfreqz(target_sos, worN=2**20, fs=samplerate)\n",
    "\n",
    "calculated_sos = np.vstack(\n",
    "    (np.array(hf_boost_shelv(K1, V_0)),\n",
    "     np.array(highpass(K2, Q)))\n",
    ")\n",
    "w2, h2 = scipy.signal.sosfreqz(calculated_sos, worN=2**20, fs=samplerate)\n",
    "\n",
    "fig, ax1 = plt.subplots()\n",
    "ax1.set_title('Digital filter frequency response')\n",
    "\n",
    "ax1.plot(w, 20 * np.log10(abs(h)), label='EBU')\n",
    "ax1.plot(w2, 20 * np.log10(abs(h2)), label='jkbd')\n",
    "\n",
    "ax1.set_xlabel('Frequency [Hz]')\n",
    "ax1.set_xscale('log')\n",
    "ax1.set_ylabel('Magnitude [dB]')\n",
    "\n",
    "plt.xlim(20, samplerate/2.0)\n",
    "plt.ylim(-12, 6)\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ae6c75dd",
   "metadata": {},
   "source": [
    "What is the maximum difference beween the spectra?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "2acb641b",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.04334313168283899 dB\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/tmp/ipykernel_28555/3750599922.py:1: RuntimeWarning: divide by zero encountered in log10\n",
      "  diffs = np.abs(20 * np.log10(abs(h)) - 20 * np.log10(abs(h2)))\n",
      "/tmp/ipykernel_28555/3750599922.py:1: RuntimeWarning: invalid value encountered in subtract\n",
      "  diffs = np.abs(20 * np.log10(abs(h)) - 20 * np.log10(abs(h2)))\n"
     ]
    }
   ],
   "source": [
    "diffs = np.abs(20 * np.log10(abs(h)) - 20 * np.log10(abs(h2)))\n",
    "m = np.max(diffs[np.logical_not(np.isnan(diffs))])\n",
    "print(f'{m} dB')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "475e7e00",
   "metadata": {},
   "source": [
    "## Is it correct for other samplerates?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "1aebcd19",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/tmp/ipykernel_28555/1822753161.py:23: RuntimeWarning: divide by zero encountered in log10\n",
      "  ax1.plot(w, 20 * np.log10(abs(h)))\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYcAAAEaCAYAAAD65pvjAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAAtSklEQVR4nO3dd5xU9b3/8ddnK7BLZwEpSxGQYsUVNWhiFKOxRk2xxha5MTE3Jppq9Jdmyk29uZqi3qhRY7kRYi8Qu6JSRLqISF1YOgss2z+/P85ZMm4ZZpedPTM77+fjMY+d077fz8zZmc+c7/me7zF3R0REJFZW1AGIiEjqUXIQEZEmlBxERKQJJQcREWlCyUFERJpQchARkSaUHAQAM/uzmd3c3uu2sL2b2agWlg0ws1fMbJeZ/cbMvm9md4XLhofb5rS17gRi+0j9yapHJNUl7UMmqcPMVgEDgFqgDlgC/A24w93rAdz9y4mWF7uumZ0E3O/uQ9op3KnAFqCH7+ciHDN7Kaz7rnaqu1X1i3RmOnLIHGe7e3dgGPAL4DvA/0YbUrOGAUs64ovZzLJbW38yj1pSUQvvkWQCd9ejkz+AVcCURvMmAfXAoeH0PcBPY5Z/G9gAlAJfAhwYFbsuUADsDcvZHT4GhWXPAnaEZdwG5MWUva+sRjHdA9QA1WFZU4AfEhwdAAwPt80BbiU4CqoM170tXGcsMAPYBrwHfL5R+X8Cngb2NPOetFT/P4D7gfLwvehJkFg3AOvD9yI7LCMb+DXB0cdK4KsNMTe3L2JfXzh9HPBG+N69C5wUs+wl4CfA68Au4HmgX8zyE2K2XQtcARwDlDXUH653ATC/hf+VJu9RuE8fBTYDHwL/2ej/aE743pQBv220r6YS/A9tAG6I2S4f+H24rDR8nh8uOwlYB9wAbAq3vTJm2zMIjn53he//jTHLzgLmh+/BG8DhUX/+0vUReQB6dMBObiY5hPPXANeGz+8hTA7A6cBGYALQDbiPZpJD+PwkYF2jco8Ov+Rywi+JpcD1McubTQ6Nyw6nf0gzySGcfgn4Usy6BeGX4pVh3RMJvqQnxJS9E5hMcNTcJcH6a4DPhNt0Bf4J/CWsrz/wNvAf4fpfBpYBQ4E+wIskmByAwcDW8MsvCzg1nC6Keb0fAGPCOF4CfhEuKyb4srwIyAX6AkeGy5YAn46pczoxX9TNvP7Y96gbMBe4BcgDRhIkvdPC9WcBl4XPC4HjGu2rB8P36TCC5DIlXP5j4M3w/Ssi+CL/Scz/VG24Tm74flQAvcPlG4ATw+e9gYnh84kEyeRYgiR9efh+50f9GUzHh5qVMlspwRdYY58H7nb3xe5eAfyoNYW6+1x3f9Pda919FcEX6ScOONr9OwtY5e53h3XPI/jF+9mYdR5z99fdvd7dKxMsd5a7/9OD8zM9gE8TJLs97r4J+B1wYbju54Hfu/tad98G/LwV8V8KPO3uT4fxzSD4VX5GzDp3u/tyd98LPAIcGc6/BJjp7g+6e427b3X3+eGye8OyMbM+wGnA3+PEse89IvhSL3L3H7t7tbuvBO6Meb01wCgz6+fuu939zUZl/Sh8nxYCdxMkr4Z4f+zum9x9M8H/2GUx29WEy2vc/WmCI7lDYpaNN7Me7r493M8A1wB/cfe33L3O3e8Fqgh+qEgrKTlktsEEzS+NDSL4Bd5gbTPrtMjMxpjZk2a20czKgZ8B/doeZsKGAcea2Y6GB8GX0MCYdVr1WprZZhjBr9kNMXX8heAXMDR971a3op5hwOcaxX8CcFDMOhtjnlcQ/FqH4EjlgxbKvR8428wKCZLXq+6+IU4cjV/voEYxfZ+ggwPA1QRHMsvMbLaZnRWnrNUE7w/h39UtLAPY6u61MdOxr/UCgoS52sxeNrPjY2K9oVGsQxuVKwnKqJNr8m9mdgxBcnitmcUbgNjeR0PjFNXcids/Ae8AF7n7LjO7no/+em8vjeteC7zs7qe2YpvW1rOW4Ndov0ZfXg028NH3q7jR8j0ETTUNGieu+9z9mjbEuJag/b8Jd19vZrOA8wh+nf9pP2U1fr0fuvvoFsp+H7jIzLKA84F/mFnfmFWGEjSzQfBelIbPSwm+zBc3syx+cO6zgXPNLBe4juAIamgY663ufmsi5Uh8OnLIMGbWI/x19xBBW/fCZlZ7BLjSzMaZWTeC9uaWlAF9zaxnzLzuBCcod5vZWODadgq/ubpHxkw/CYwxs8vMLDd8HGNm49qrwvAX9/PAb8L3MsvMDjazhmazR4D/NLMhZtYb+G6jIuYDF4axlfDRpNnwC/80M8s2sy5mdpKZJdJN+AFgipl93sxyzKyvmR0Zs/xvBJ0MDiM455Cot4FyM/uOmXUN4zo0/HGBmV1qZkVhE9SOcJu6mO1vNrNuZjaB4FzQw+H8B4EfmFmRmfUj+B+7f3/BmFmemV1iZj3dvYbg/6yhvjuBL5vZsRYoMLMzzax7K16vhJQcMscTZraL4NfVTcBvCT6sTbj7M8AfCE6mriA46QjBL+bG6y4j+KCvDA/lBwE3AhcTnCC9k39/IbS3/wY+a2bbzewP7r4L+BRBe3gpQRPMLwl6xrSnLxKcnF0CbCfozdTQ9HMn8BxBT6N5wLRG294MHBxu9yNi2v7dfS1wLkGzzWaCffUtEvicuvsagqaWGwiaCucDR8SsMp3gl/p0d9+T6At19zrgbIJzGx8SnOC/i6DHFgSdFxab2W6C/XFho3M5LxP8D/0L+LW7Px/O/ynB+ZQFwEKC9+qnCYZ1GbAqbLL8MuH5FHefQ3De4TaC93cFQY8taQNz13U+El/4y3sRQa+P5ppSpAVmNpzgSzU36vfOzD4g6FU1swPqGk6KvG5pGx05SLPM7LzwEL43wa/vJ/QhT19mdgHBuYQXoo5F0kPKJgcz62Vm/zCzZWa2NKZHgnSM/yBo2viAoE03WecNJMnCYUb+BHw1PDcgsl8p26xkZvcSdLm7y8zygG7uviPisEREMkJKJgcz60FwQm+kp2KAIiKdXKo2K40kaNK428zeMbO7zKwg6qBERDJFqh45lBCMuzLZ3d8ys/8Gyt395ph1phIM6kVBQcHRY8eOjSZYEZE0NXfu3C3uXtTcslRNDgOBN919eDh9IvBddz+zufVLSkp8zpw5HRihiEj6M7O57l7S3LKUbFZy943AWjNrGGjrFIILjkREpAOk8thKXwMeCHsqraSFq3lFRKT9pWxyCIcbbvZwR0REkislm5VERCRaSg4iItKEkoOIiDSh5CAiIk0oOYiISBNKDiIi0oSSg4iINKHkICIiTSg5iIhIE0oOIiLShJKDiIg0oeQgIiJNKDmIiEgTSg4iItKEkoOIiDSRssnBzLLN7B0zezLqWEREMk3KJgfg68DSqIMQEclEKXknODMbApwJ3Ap8M+JwRERa5O4tzI+zTVvKi7tNvLribdmylEwOwO+BbwPdI45DRFJAVW0dOypq2Lanmu0V1WzfU8OeqloqqmvZU13H3uo69lTXhn/rqKqpo7beqamrDx9ObV091eHfhnl19U69O/UOEPytd8dj/rp/dL43Wq+zSrnkYGZnAZvcfa6ZnRRnvanAVIDi4uKOCU5EkmJnRQ3vb9rFyi17WL99L6U79lK6cy+lOyrZvKuK3VW1cbfPMijIy6FrXjbd8rLpkptNbnYWOdlGbnYWXXKzyO2SQ05WFnk5Rk5WFrnZWWRnQXaWAUaWQZYZFvPXGuZnNZpuWG6GAWbNxxUsbWFZy4ta3CruNvEWtuC6X8Ypr6VDmKiY2c+By4BaoAvQA5jm7pe2tE1JSYnPmTOngyIUkQOxeVcV89ZsZ96a7Sxav5P3y3azaVfVvuVm0L97PoN6dWVQr670755Pn2559C7Io09BHr265dK7Wx6F+TkU5OfQLS+b/JysNn05Zjozm+vuJc0tS7kjB3f/HvA9gPDI4cZ4iUFEUtve6jpmrdzCi8s288r7m1m9tQKA3Gxj7MAenDi6iDEDChkzoDsjiwo4qGdX8nJSua9MZki55CAi6a+qto4Xl21i+jvrefG9zVTX1tM1N5vJo/pyybHFHD2sNxMG9aRLbnbUoUoLUjo5uPtLwEsRhyEiCVq9dQ/3vLGKafPWs3NvDf0K87l4UjGnjOvPpBF9yM9RMkgXKZ0cRCQ9zF+7gz++uIIZS8vINuOMww7igqOHMPngvuRkq4koHSk5iEibvbdxF795/j2eX1JG7265fPWkUVx2/DAG9OgSdWhygJQcRKTVdu6t4VfPLeOBt9ZQmJfDDaeO4coTRlCYr6+UzkJ7UkQS5u48/m4pP3lyCdv2VHP58cP5+imj6V2QF3Vo0s6UHEQkIdv2VPPdRxfw/JIyjhjai3uunMShg3tGHZYkiZKDiOzXq+9v5oZH3mVHRQ03nTGOq04YEV5ZLJ2VkoOItKi+3vmfF1bwu5nLGdW/kLuvPIYJg3S0kAmUHESkWburavnmw/N5fkkZ508czK2fOYyuebpOIVMoOYhIE2u3VXDVPbNZuWUPt5w1nisnD9fYRRlGyUFEPmLR+p1ccfdsaurq+dtVk5g8ql/UIUkElBxEZJ83Vmxh6n1z6dElh4emHs+o/rqlSqZSchARAJ5cUMo3Hp7PyH6F3HvVJAb21FXOmUzJQUR4dO46bvzHuxwzrA93frGEnt1yow5JIqbkIJLhHpmzlu88uoDJB/fjzi+WqEeSAEoOIhnt4dlr+O60hZwwKkgMur+CNNBYuiIZ6u9vreE7jy7k46OLlBikiZRMDmY21MxeNLOlZrbYzL4edUwincl9b67m+9MX8slDivjLZUcrMUgTqdqsVAvc4O7zzKw7MNfMZrj7kqgDE0l3f5u1ilseW8yUcf25/ZKJujubNCsljxzcfYO7zwuf7wKWAoOjjUok/TUkhlPHD+CPlxytxCAtSsnkEMvMhgNHAW81mj/VzOaY2ZzNmzdHEptIOolNDLdfPJG8nJT/+EuEUvq/w8wKgUeB6929PHaZu9/h7iXuXlJUVBRNgCJpQolBWitl/0PMLJcgMTzg7tOijkckXd237xyDEoMkLiX/SywY/vF/gaXu/tuo4xFJV/fNWsXNYWL44yVKDJK4VP1PmQxcBpxsZvPDxxlRByWSTu55/cMwMfRXYpBWS8murO7+GqDB40XawN257YUV/GbGcj41fgC3qSlJ2iAlk4OItI2787Onl3Lnqx9y/sTB/NcFh5OTrcQgrafkINJJ1NbVc9P0RTw8Zy1XfGw4t5w1nqwsHYBL2yg5iHQCuypruO7v7/Dy8s187eRRfPPUMbqtpxwQJQeRNLd+x16uvmc272/azc/PP4yLJhVHHZJ0AkoOImls9qptfOWBeVTW1HHvlZM4YbTu9yztQ8lBJA25O3e+upJfPvseQ3t35e9fOpbRA3S/Z2k/Sg4iaWbbnmq+++gCnl9SxqcPHcgvP3s4Pbrotp7SvpQcRNLIc4s3ctP0hezcW8MPzhzH1SeM0IlnSQolB5E0sKm8klufXspj80sZf1AP7rv6WMYd1CPqsKQTU3IQSWFVtXX89bVV3PbC+9TUOV8/ZTTXnTyKXF3YJkmm5CCSgqpr65k2bx23vbiCddv3MmXcAH5w5jiG9yuIOjTJEEoOIilkV2UN0+at545XVrJ+x16OGNKTn513GB8fo3uWSMdSchCJmLuzaH05f397DY/NX09FdR0Ti3tx63mH8okxRTrhLJFQchCJQH29s7i0nKcXbeCpBRtYs62C/JwszjliEJceN4wjhvaKOkTJcEoOIh3A3Vm1tYK3Vm7ltRVbeOODrWzbU012lvGxg/vy1U8ezOkTDqJnN12vIKkhZZODmZ0O/DeQDdzl7r+IOCSRhFRU17Jy8x5WbtnDktJyFq7fwYJ1O9lVWQtA/+75nDSmiMmj+vHJsf3pU5AXccQiTaVkcjCzbOB24FRgHTDbzB539yXRRiaZzt3ZXlHDxp2VlJVXsrG8ct/z1Vsr+HDLHjaWV+5bPzfbGDuwB2cfMYjDB/fk6GG9GdW/UOcRJOWlZHIAJgEr3H0lgJk9BJwLKDlIu3B3dlfVUl5Zy46KanZW1LBjbw07KmrYXlHNzr017KioZntFTbjs38+r6+o/UpYZ9C3IZ0jvrnxsVF8OLipkRL8CRvQrYGRRAfk52RG9SpG2azE5mNkfEti+3N1/0I7xNBgMrI2ZXgccG7uCmU0FpgIUF2uI4kxVUV3L1t3VbNtTzbaKasr31rBzbw3le2sor6xlZ0UN5ZXBI5hfG0zvraHeWy63S24Wvbrm0atbLj275jKyX2HwvFsu/bt34aCeXRjQowsDe3ahf/d8XZQmnU68I4dzgVv2s/13gWQkh+aOuT/yUXb3O4A7AEpKSuJ8zCXd1Nc7W/ZUsXFnJRt2VrJhx142lFeyZVc12/ZUsXVP9b6EsLemrsVyuuRm0bNrLj265NKja/ClPqoohx5dc2Pm59Czay69ugWJoCEhdMnVr33JbPGSw+/c/d54G5tZ73aOp8E6YGjM9BCgNEl1SQR2VdawemsFq7bu2ddWv3rrHkp3BO33tY1+1udlZ9GvMI8+hXn0LchnVFEhfQqC6X4F+fQpyKN3QfDF3vClr+YckbZrMTm4++/3t3Ei67TRbGC0mY0A1gMXAhcnqS5JosqaOlZs2s3SDeUs27iLZRvLeW/jbrbsrvrIegN65DOsTwHHjujDwJ5Bs83Anl3Dv13o0y1P90MW6UDxzjl0Ab4AbAeeAL4NnAh8APzE3bckKyh3rzWz64DnCLqy/tXdFyerPmkf9fXOyi27mbd6B++s3c681Tt4f9OufW37XXKzOGRAd04eW8SIfoWM6NeNYX0LGNa3G93yUrVvhEhmiveJ/BtQAxQANwCLgNuAE4B7gLOSGZi7Pw08ncw65MC4O8vLdvPaii28vmILc1Ztozzsy9+zay5HFffitAkDGHdQDw4Z2J1hfQvI1q9/kbQQLzmMd/dDzSwHWOfunwjnP2tm73ZAbJKC9lTV8tJ7m/nX0jJeW7GFTbuC5qGR/Qo48/BBTCzuxVHFvRnZr0DNQCJpLF5yqIZ9TTyNTwa33EVEOp2dFTXMWFrGs4s28sr7m6murad3t1xOGF3ECaP6csLoIgb36hp1mCLSjuIlhyHhtQ4W85xwenDSI5NI1dbV8+qKLfxj7jpmLCmjuraeQT27cMmxxZw2YSDHDO+jJiKRTixecvhWzPM5jZY1npZOYuPOSu5/czWPzFnLpl1V9O6Wy8WTijnvqMEcPqSnhn0QyRDxurLGvcZBOpf5a3fw19c+5OmFG6hz5+RD+vO5kqGcPLY/eTm6+lck08TryvoEja5KjuXu5yQlIulQb63cyu9nvs+slVvpnp/D5R8bzuXHD6e4b7eoQxORCMVrVvp1+Pd8YCBwfzh9EbAqiTFJB3j7w238bsZyZq3cSlH3fG4+azxfOGYohfm63kBE4jcrvQxgZj9x94/HLHrCzF5JemSSFKu27OFnTy/l+SVl9CsMksIlxxZrLCER+YhEfiYWmdnImOGzRwC623ma2VVZw20vrOCvr39IbnYW3zrtEK6aPIKueUoKItJUIsnhG8BLZrYynB5OOFS2pIcXlpXx/WmL2FheyWePHsK3TzuE/j26RB2WiKSw/SYHd3/WzEYDY8NZy9y9Kt42khq27q7ix08u4bH5pRwyoDt/unQiRxUnayBdEelM4vVWmuju8wDCZNBkyIzYdSS1vPTeJm78v3fZubeGb0wZw7UnHawuqSKSsHhHDneb2Uk0f+OdBv8LHNWeAcmBqa6t51fPLePOVz9k7MDuPPCl4zhkYPeowxKRNBMvOfQE5hI/OWxu33DkQKzbXsFXHpjHgnU7uey4Ydx05jj1QhKRNonXlXV4B8YhB+itlVu59oF51NTV8+dLj+b0QwdGHZKIpLGUu+LJzH4FnE0wKuwHwJXuviPSoFLc399awy2PLaK4bzfu/GIJBxcVRh2SiKS5VDxDOQM41N0PB5YD34s4npRVX+/c+tQSvj99IZNH9WP6VyYrMYhIu0i5Iwd3fz5m8k3gs1HFkspq6ur59j8WMP2d9Vx+/DBuOXuChtAWkXaz3yMHC1xqZreE08VmNin5oQFwFfBMC3FNNbM5ZjZn8+bMOi9eUV3LNX+bw/R31nPDqWP44TlKDCLSvhJpVvojcDzBgHsAu4DbD6RSM5tpZouaeZwbs85NQC3wQHNluPsd7l7i7iVFRZkzmkdFdS1X3D2bV5Zv5mfnHcbXThmteyyISLtLpFnpWHefaGbvALj7djPLO5BK3X1KvOVmdjlwFnCKu7c4bHim2Vtdx1X3zGbOqm38/sKjOOeIQVGHJCKdVCLJocbMsgnv7WBmRUB9sgIys9OB7wCfcPeKZNWTbvZW13H1vbODoba/cKQSg4gkVSLNSn8ApgP9zexW4DXgZ0mM6TagOzDDzOab2Z+TWFdaqKmr58v3z2XWyq385vNHcO6RuoW3iCRXIgPvPWBmc4FTCK6W/oy7L01WQO4+KlllpyN35zuPLuDl5Zv5xfmHcd5RQ6IOSUQyQLyB9/rETG4CHoxd5u7bkhmYBH713HtMm7eeb0wZw4WTiqMOR0QyRLwjh7kE5xkMKAa2h897AWuAEckOLtPdN2sVf3zpAy6aVMx/nqIDKhHpOC2ec3D3Ee4+EngOONvd+7l7X4JeRNM6KsBM9fqKLfzwiSVMGdefn5w7Qd1VRaRDJXJC+hh3f7phwt2fAT6RvJBk9dY9fOWBeRxcVMDvLzyKnOxUHOVERDqzRLqybjGzHwD3EzQzXQpsTWpUGWx3VXD1M8CdXyyhMD/lRjgRkQyQyE/Si4Aigu6s/wT68++rpaUduTs3PvIuKzbt5vaLJzKsb0HUIYlIhkqkK+s24OsdEEvGu/eNVTy7eCM3nTGOE0b3izocEclg+00OZvYi4dXRsdz95KRElKEWrd/Jz55exilj+/OlE9URTESilUiD9o0xz7sAFxAMiCftZHdVLdf9fR59CvL41eeOUM8kEYlcIs1KcxvNet3MXk5SPBnppukLWbOtgoemHk+fggMa01BEpF0k0qwUe6V0FnA0oBsUt5Mn3i3lsfmlfGPKGCaN6LP/DUREOkAizUqxV0rXAh8CVyczqEyxaVclNz+2iCOG9uKrnzw46nBERPZJJDmMc/fK2Blmlp+keDKGu/P9aYuoqK7jN587XBe6iUhKSeQb6Y1m5s1q70AyzbR565m5tIxvn3YIo/p3jzocEZGPiDcq60BgMNDVzI4iaFYC6AF064DYOq1N5ZX88InFHDO8N1dOVrdVEUk98ZqVTgOuAIYAv42Zvwv4fhJjAsDMbgR+BRS5+5Zk19eRfvzkEqpq6/mvzx5Bdpa6rYpI6mkxObj7vcC9ZnaBuz/agTFhZkOBUwmGBu9UXl6+mScXbOCbp45hRD8NjyEiqSles9Kl7n4/MNzMvtl4ubv/tpnN2svvgG8DjyWxjg5XWVPHzf9cxMiiAv7jEyOjDkdEpEXxmpUaftYWdkQgDczsHGC9u78b70phM5sKTAUoLk6PO6T9zwvvs2ZbBQ9ecxz5OdlRhyMi0qJ4zUp/Cf/+qL0rNbOZNH8h3U0E5zM+tb8y3P0O4A6AkpKSJmM/pZoVm3ZzxysruWDiEI4/uG/U4YiIxJXIFdJFwDXA8Nj13f2qtlbq7lNaqOswgtuPNhw1DAHmmdkkd9/Y1vpSwa1PLaFLTjbfO2Ns1KGIiOxXIhfBPQa8CswE6pIZjLsvJLhfBABmtgooSffeSi+9t4kX39vMTWeMo1+hrh8UkdSXSHLo5u7fSXoknVRNXT0/fWopw/t24/KPDY86HBGRhCRyhfSTZnZG0iNphrsPT/ejhr+/tYYVm3Zz05njycvREBkikh4S+bb6OkGC2Gtm5Wa2y8zKkx1YZ7CjoprfzVzO5FF9mTKu//43EBFJEYncz0ED/7TR7S+uoHxvDT84c7xu4CMiaSWR3koTm5m9E1jt7rojXAs27NzLvbNWc95RQxh3UI+owxERaZVETkj/EZgILAynDwPeBfqa2Zfd/flkBZfO/vCvFbg7108ZHXUoIiKtlsg5h1XAUe5+tLsfDRwJLAKmAP+VvNDS14db9vDInLVcPKmYoX00gK2IpJ9EksNYd1/cMOHuSwiSxcrkhZXefjtjOXnZWVx3so4aRCQ9JdKs9J6Z/Ql4KJz+ArA8vBtcTdIiS1NLSst54t1SvvrJgynqrgveRCQ9JXLkcAWwArge+AawMpxXA3wySXGlrd/OWE6PLjlM/bjuCS0i6SuRrqx7gd+Ej8Z2t3tEaWxx6U5mLi3jm6eOoWfX3KjDERFps0S6so4Gfg6MB7o0zHd33ZCgkdteWEH3/BwNkyEiaS+RZqW7gT8BtQTNSH8D7ktmUOloedkunlm0kSsmD9dRg4ikvUSSQ1d3/xdg7r7a3X8InJzcsNLP7S+uoFteNldNHhF1KCIiByyR3kqVZpYFvG9m1wHriRlWW4LrGp54t5RrThxJ74K8qMMRETlgiRw5XA90A/4TOBq4DLg8iTGlndtfXEFeThZfOlGnYUSkc0ikt9Ls8Olu4MrkhpN+Snfs5Z/vrOfS44bpugYR6TRaTA5m9ni8Dd39nPYPZ1/dXwOuIzgJ/pS7fztZdR2ou1//EAeu+biOGkSk84h35HA8sBZ4EHgL6JAxp83sk8C5wOHuXmVmKXt+o7yyhgffXsuZhx3E4F5dow5HRKTdxEsOA4FTgYuAi4GngAdjx1lKkmuBX7h7FYC7b0pyfW328Ntr2V1VyzU61yAinUyLJ6Tdvc7dn3X3y4HjCIbQeCls8kmmMcCJZvaWmb1sZsckub42qamr56+vf8hxI/tw2JCeUYcjItKu4p6QDgfXO5Pg6GE48Adg2oFWamYzCY5MGrspjKk3QUI6BnjEzEa6uzcqYyowFaC4uPhAQ2q1pxZsYMPOSm4979AOr1tEJNninZC+FzgUeAb4kbsvaq9K3X1KnHqvBaaFyeBtM6sH+gGbG5VxB3AHQElJiTcpKIncnTtfXcmo/oWcNCZlT4mIiLRZvOscLiNo4vk68IaZlYePXWZWnsSY/kl4BbaZjQHygC1JrK/VZn2wlcWl5XzphBFkZene0CLS+bR45ODuiVwglwx/Bf5qZouAauDyxk1KUbv7jVX0KcjjM0cNjjoUEZGkSGT4jA7l7tXApVHH0ZK12yr419Iyrj3pYLrkZkcdjohIUkR1dJC2HnhrDQCXHDss4khERJJHyaEVKmvqeHj2Gk4dP4BBuuhNRDoxJYdWeHLBBrZX1HD58cOjDkVEJKmUHBLk7tz7xipG9S/k+IP7Rh2OiEhSKTkkaP7aHSxcv5PLjx+GmbqvikjnpuSQoPtmraYwP4fzJg6JOhQRkaRTckjAtj3VPLlgAxdMHExhfsr1/hURaXdKDgmYNm8d1XX1XKzuqyKSIZQc9sPdeWj2Wo4q7sUhA7tHHY6ISIdQctiPeWu2s2LTbi48ZmjUoYiIdBglh/148O21FORlc9bhg6IORUSkwyg5xFFeWcNTCzZwzpGDKNCJaBHJIEoOcTw+v5S9NXVceEzH30xIRCRKSg5xPDR7DWMHdudw3QZURDKMkkMLFq3fyaL15Vw0qVhXRItIxkm55GBmR5rZm2Y238zmmNmkKOJ4ePZa8nOy+MyRuqGPiGSelEsOwH8R3LP6SOCWcLpDVdXW8dj89Zx+6EB6dsvt6OpFRCKXisnBgR7h855AaUcH8MLSTZRX1nKBxlESkQyViv0zrweeM7NfEySvj3V0AI/OW0//7vlMHtWvo6sWEUkJkSQHM5sJDGxm0U3AKcA33P1RM/s88L/AlGbKmApMBSgubr+uplt3V/HSe5u4+oQRZGfpRLSIZKZIkoO7N/myb2BmfwO+Hk7+H3BXC2XcAdwBUFJS4u0V2+PvllJb75yvJiURyWCpeM6hFPhE+Pxk4P2OrHzavPVMGNRDg+yJSEZLxXMO1wD/bWY5QCVh01FHWF62i4Xrd3LLWeM7qkoRkZSUcsnB3V8Djo6i7mnz1pOdZZxzpAbZE5HMlorNSpGoq3emv7OOk8YU0a8wP+pwREQipeQQeuODLZSVV+lEtIgISg77TJ+3nu5dcjhlXP+oQxERiZySA1BZU8fzS8r49KED6ZKbHXU4IiKRU3IAXly2id1VtZxzhAbZExEBJQcAnlhQSr/CPI4b2SfqUEREUkLGJ4ddlTX8a+kmzjzsIHKyM/7tEBEBlByYubSMqtp6zj5C1zaIiDTI+OTw+PxSBvfqysTi3lGHIiKSMjI6OWzfU82r72/hrMMPIksjsIqI7JPRyeGZRRuprXc1KYmINJLRyeGJd0sZWVTAhEE99r+yiEgGydjkUFZeyZsfbuXswwdhpiYlEZFYGZscnlywAXfUpCQi0oyMTQ5PL9zAuIN6MKp/YdShiIiknEiSg5l9zswWm1m9mZU0WvY9M1thZu+Z2WnJqH/jzkrmrt7OmYc1dxtrERGJ6mY/i4Dzgb/EzjSz8cCFwARgEDDTzMa4e117Vv7sog0AnH7oQe1ZrIhIpxHJkYO7L3X395pZdC7wkLtXufuHwApgUnvX/8yijYwZUKgmJRGRFqTaOYfBwNqY6XXhvHazeVcVb6/apqMGEZE4ktasZGYzgeYa9W9y98da2qyZed5C+VOBqQDFxcUJx/X8ko24wxk63yAi0qKkJQd3n9KGzdYBQ2OmhwClLZR/B3AHQElJSbMJpDnPLNzIiH4FHDKgexvCExHJDKnWrPQ4cKGZ5ZvZCGA08HZ7Fb59TzWzVm7l04cO1IVvIiJxRNWV9TwzWwccDzxlZs8BuPti4BFgCfAs8NX27Kk0Y0kZdfXOGYfpfIOISDyRdGV19+nA9BaW3Qrcmox6n1m0gSG9u2osJRGR/Ui1ZqWk2bm3htdWbOGMww5Sk5KIyH5kTHJ4YVkZNXXO6Yeql5KIyP5kTHJ4euFGDurZhSOH9Io6FBGRlJcRyaGiupZXlm/mtAkDdcc3EZEEZERyeGX5Zqpq6/nUhAFRhyIikhYyIjk8v6SMnl1zmTS8T9ShiIikhU6fHGrr6nlh2SZOHtufnOxO/3JFRNpFp/+2nLN6Ozsqajh1vJqUREQS1emTw4wlZeRlZ/HxMUVRhyIikjY6dXJwd2YsKeNjo/pSmB/VfY1ERNJPp04Oy8t2s2ZbhZqURERaqVMnhxlLNgIwZZySg4hIa3Ty5FDGEUN7MaBHl6hDERFJK502OZSVV/Luup18Sk1KIiKt1mmTw4wlZQBKDiIibdCpk8Pwvt0Y1b8w6lBERNJOVHeC+5yZLTazejMriZl/qpnNNbOF4d+T21L+rsoa3vhgC6eOH6B7N4iItEFUnf8XAecDf2k0fwtwtruXmtmhwHPA4NYW/vLyzdTUOaeO170bRETaIqrbhC4Fmvyqd/d3YiYXA13MLN/dq1pT/r+WbqJ3t1yOHtb7gGMVEclEqXzZ8AXAOy0lBjObCkwNJyvNbHHjdXL+X5vq7QnsTPJ2+1s33vLWLmtuXj+Co7QotfV9bs+yEt0ukfW0zzqmrKj3WWvmp8M+G9biEndPygOYSdB81Phxbsw6LwElzWw7AfgAODjBuu5ox7jbVFZrttvfuvGWt3ZZC/PmJGu/J/t9jmKfJbKe9llm7LPWzE/3fZa0Iwd3n9KW7cxsCDAd+KK7f5DgZk+0pa52Lqs12+1v3XjLW7usPd+b9pRO+yyR9bTPOqasqPdZa+dHrc1xWZhdImFmLwE3uvuccLoX8DLwY3d/NLLAOjkzm+PuJftfU1KF9ln6Sfd9FlVX1vPMbB1wPPCUmT0XLroOGAXcbGbzw0f/KGLs5O6IOgBpNe2z9JPW+yzSIwcREUlNnfYKaRERaTslBxERaULJQUREmlByEMzsM2Z2p5k9Zmafijoe2T8zG2dmfzazf5jZtVHHI/tnZgXhmHFnRR1LIpQcOikz+6uZbTKzRY3mn25m75nZCjP7LoC7/9PdrwGuAL4QQbhCq/fZUnf/MvB5IG27S6az1uyv0HeARzo2yrZTcui87gFOj51hZtnA7cCngfHARWY2PmaVH4TLJRr30Ip9ZmbnAK8B/+rYMCV0DwnuLzObAiwByjo6yLZScuik3P0VYFuj2ZOAFe6+0t2rgYeAcy3wS+AZd5/X0bFKoDX7LFz/cXf/GHBJx0Yq0Or99UngOOBi4BozS/nv3lQeeE/a32Bgbcz0OuBY4GvAFKCnmY1y9z9HEZw0q9l9ZmYnEQx7nw883fFhSQua3V/ufh2AmV0BbHH3+ghiaxUlh8zS3J2P3N3/APyho4ORhLS0z14iGLhSUkuz+2vfE/d7Oi6UA5PyhzbSrtYBQ2OmhwClEcUiidE+Sy+dZn8pOWSW2cBoMxthZnnAhcDjEcck8WmfpZdOs7+UHDopM3sQmAUcYmbrzOxqd68lGNzwOWAp8Ii7N7lJkkRD+yy9dPb9pYH3RESkCR05iIhIE0oOIiLShJKDiIg0oeQgIiJNKDmIiEgTSg4iItKEkoN0KmZWZ2bzYx7Do46pPZjZFWa22czuCqdPMrMnG61zj5l9Nk4ZvzKzjWZ2Y7LjlfSnsZWks9nr7kc2t8DMjODanpQf9KwFDzcM4NYW7v4tM9vTngFJ56UjB+nUzGy4mS01sz8C84ChZvYtM5ttZgvM7Ecx694U3qRlppk92PAL28xeMrOS8Hk/M1sVPs8Of403lPUf4fyTwm3+YWbLzOyBMDFhZseY2Rtm9q6ZvW1m3c3sVTM7MiaO183s8AN4zSUxR04LzUxXukqr6chBOpuuZjY/fP4h8A3gEOBKd/9KeBvU0QTj7hvwuJl9HNhDMA7OUQSfi3nA3P3UdTWw092PMbN84HUzez5cdhQwgWDQtdeByWb2NvAw8AV3n21mPYC9wF0Ed+G73szGAPnuviCB13pizGsFKAaedPc5wJEQNCUBzyZQlshHKDlIZ/ORZqXwnMNqd38znPWp8PFOOF1IkCy6A9PdvSLcLpHB0j4FHB7Tzt8zLKsaeNvd14VlzQeGAzuBDe4+G8Ddy8Pl/wfcbGbfAq4iuMNYIl519333Izazj2xnZp8HJoZxirSKkoNkgth2dgN+7u5/iV3BzK4nZtz9Rmr5dxNsl0Zlfc3dn2tU1klAVcysOoLPmjVXh7tXmNkMgjuGtcs9oc1sAvAj4OPuXneg5Unm0TkHyTTPAVeZWSGAmQ02s/7AK8B5ZtbVzLoDZ8dsswo4Onz+2UZlXWtmuWFZY8ysIE7dy4BBZnZMuH53M2v4gXYXwQ2XZrt741tPtoqZ9SS4PeUX3X3zgZQlmUtHDpJR3P15MxsHzArPEe8GLnX3eWb2MDAfWA28GrPZr4FHzOwy4IWY+XcRNBfNC084bwY+E6fuajP7AvA/ZtaV4HzDFGC3u881s3Lg7nZ4mZ8BhgF3hq+RlnpwibREQ3aLNMPMfkjwpf3rDqpvEMFtP8c219U2vPdwyYF0ZQ3L+SEd+LokfalZSSRiZvZF4C3gpjjXYOwFPt1wEVwb6/kVcCkfPQcj0iwdOYiISBM6chARkSaUHEREpAklBxERaULJQUREmlByEBGRJpQcRESkif8PUXPwxwU2f0MAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "samplerate = 48000\n",
    "\n",
    "def decibel2factor(decibel):\n",
    "    return np.power(10, (decibel/20.0))\n",
    "\n",
    "def freq2k(freq, samplerate):\n",
    "    return np.tan(np.pi * freq / samplerate)\n",
    "\n",
    "K1 = freq2k(1681.7632251028442, samplerate)\n",
    "V_0 = decibel2factor(3.9997778685513232)\n",
    "K2 = freq2k(38.135470876002174, samplerate)\n",
    "Q = 0.5003270373223665\n",
    " \n",
    "sos = np.vstack(\n",
    "    (np.array(hf_boost_shelv(K1, V_0)),\n",
    "     np.array(highpass(K2, Q)))\n",
    ")\n",
    "w, h = scipy.signal.sosfreqz(sos, worN=2**20, fs=samplerate)\n",
    "\n",
    "fig, ax1 = plt.subplots()\n",
    "ax1.set_title('Digital filter frequency response')\n",
    "\n",
    "ax1.plot(w, 20 * np.log10(abs(h)))\n",
    "\n",
    "ax1.set_xlabel('Frequency [Hz]')\n",
    "ax1.set_xscale('log')\n",
    "ax1.set_ylabel('Magnitude [dB]')\n",
    "\n",
    "plt.xlim(20, samplerate/2.0)\n",
    "plt.ylim(-12, 6)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2edea3a6",
   "metadata": {},
   "source": [
    "Looks okay."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c271b8ea",
   "metadata": {},
   "source": [
    "The magnitude at 997Hz is:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "9e6afe43",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(array([ 997., 1000.]), array([1.07762056, 1.07845159]))\n",
      "1.0776205611408702\n",
      "1.0784515947925883\n"
     ]
    }
   ],
   "source": [
    "w2, h2 = scipy.signal.sosfreqz(sos, np.array([997.0, 1000.0]),\n",
    "                               fs=samplerate)\n",
    "print((w2, np.abs(h2)))\n",
    "for mag in np.abs(h2):\n",
    "    print(mag)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3a353372",
   "metadata": {},
   "source": [
    "So we need to normalize with the factor"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "505a11c2",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.9273671710547968"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "1/1.0783215442731167"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a6e449a6",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.10.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}