{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Compare Estimators\n\nIn this tutorial we investigate different estimator implementation techniques\nand compare their performance.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import timeit\nimport matplotlib.pyplot as plt\nimport numpy as np\nimport cbadc"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Analog System\n\nWe will commit to a leap-frog control-bounded analog system throughtout\nthis tutorial.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Determine system parameters\nN = 6\nM = N\nbeta = 6250\n# Set control period\nT = 1.0 / (2.0 * beta)\n# Adjust the feedback to achieve a bandwidth corresponding to OSR.\nOSR = 128\nomega_3dB = 2 * np.pi / (T * OSR)\n\n# Instantiate analog system.\nbeta_vec = beta * np.ones(N)\nrho_vec = - omega_3dB ** 2 / beta * np.ones(N)\nGamma = np.diag(-beta_vec)\nanalog_system = cbadc.analog_system.LeapFrog(beta_vec, rho_vec, Gamma)\n\nprint(analog_system, \"\\n\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Analog Signal\n\nWe will also need an analog signal for conversion.\nIn this tutorial we will use a Sinusodial signal.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Set the peak amplitude.\namplitude = 1.0\n# Choose the sinusodial frequency via an oversampling ratio (OSR).\nfrequency = 1.0 / (T * OSR * (1 << 0))\n\n# We also specify a phase an offset these are hovewer optional.\nphase = 0.0\noffset = 0.0\n\n# Instantiate the analog signal\nanalog_signal = cbadc.analog_signal.Sinusodial(amplitude, frequency, phase, offset)\n\nprint(analog_signal)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Simulating\n\nEach estimator will require an independent stream of control signals.\nTherefore, we will next instantiate several digital controls and simulators.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Set simulation precision parameters\natol = 1e-6\nrtol = 1e-12\nmax_step = T / 10.\n\n# Instantiate digital controls. We will need four of them as we will compare\n# four different estimators.\ndigital_control1 = cbadc.digital_control.DigitalControl(T, M)\ndigital_control2 = cbadc.digital_control.DigitalControl(T, M)\ndigital_control3 = cbadc.digital_control.DigitalControl(T, M)\ndigital_control4 = cbadc.digital_control.DigitalControl(T, M)\nprint(digital_control1)\n\n# Instantiate simulators.\nsimulator1 = cbadc.simulator.StateSpaceSimulator(\n    analog_system,\n    digital_control1,\n    [analog_signal],\n    atol=atol,\n    rtol=rtol,\n    max_step=max_step\n)\nsimulator2 = cbadc.simulator.StateSpaceSimulator(\n    analog_system,\n    digital_control2,\n    [analog_signal],\n    atol=atol,\n    rtol=rtol,\n    max_step=max_step\n)\nsimulator3 = cbadc.simulator.StateSpaceSimulator(\n    analog_system,\n    digital_control3,\n    [analog_signal],\n    atol=atol,\n    rtol=rtol,\n    max_step=max_step\n)\nsimulator4 = cbadc.simulator.StateSpaceSimulator(\n    analog_system,\n    digital_control4,\n    [analog_signal],\n    atol=atol,\n    rtol=rtol,\n    max_step=max_step\n)\nprint(simulator1)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Default, Quadratic Complexity, Estimator\n\nNext we instantiate the quadratic and default estimator\n:py:class:`cbadc.digital_estimator.DigitalEstimator`. Note that during its\nconstruction, the corresponding filter coefficients of the system will be\ncomputed. Therefore, this procedure could be computationally intense for a\nanalog system with a large analog state order or equivalently for large\nnumber of independent digital controls.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Set the bandwidth of the estimator\nG_at_omega = np.linalg.norm(\n    analog_system.transfer_function_matrix(np.array([omega_3dB])))\neta2 = G_at_omega**2\nprint(f\"eta2 = {eta2}, {10 * np.log10(eta2)} [dB]\")\n\n# Set the batch size\nK1 = 1 << 14\nK2 = 1 << 14\n\n# Instantiate the digital estimator (this is where the filter coefficients are\n# computed).\ndigital_estimator_batch = cbadc.digital_estimator.DigitalEstimator(\n    analog_system, digital_control1, eta2, K1, K2)\ndigital_estimator_batch(simulator1)\n\nprint(digital_estimator_batch, \"\\n\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Visualize Estimator's Transfer Function (Same for Both)\n\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Logspace frequencies\nfrequencies = np.logspace(-3, 0, 100)\nomega = 4 * np.pi * beta * frequencies\n\n# Compute NTF\nntf = digital_estimator_batch.noise_transfer_function(omega)\nntf_dB = 20 * np.log10(np.abs(ntf))\n\n# Compute STF\nstf = digital_estimator_batch.signal_transfer_function(omega)\nstf_dB = 20 * np.log10(np.abs(stf.flatten()))\n\n# Signal attenuation at the input signal frequency\nstf_at_omega = digital_estimator_batch.signal_transfer_function(\n    np.array([2 * np.pi * frequency]))[0]\n\n# Plot\nplt.figure()\nplt.semilogx(frequencies, stf_dB, label='$STF(\\omega)$')\nfor n in range(N):\n    plt.semilogx(frequencies, ntf_dB[0, n, :], label=f\"$|NTF_{n+1}(\\omega)|$\")\nplt.semilogx(frequencies, 20 * np.log10(np.linalg.norm(\n    ntf[0, :, :], axis=0)), '--', label=\"$ || NTF(\\omega) ||_2 $\")\n\n# Add labels and legends to figure\nplt.legend()\nplt.grid(which='both')\nplt.title(\"Signal and noise transfer functions\")\nplt.xlabel(\"$\\omega / (4 \\pi \\\\beta ) $\")\nplt.ylabel(\"dB\")\nplt.xlim((frequencies[1], frequencies[-1]))\nplt.gcf().tight_layout()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## FIR Filter Estimator\n\nSimilarly as for the previous estimator the\n:py:class:`cbadc.digital_estimator.FIRFilter` is initalized. Additionally,\nwe visualize the decay of the $\\|\\cdot\\|_2$ norm of the corresponding\nfilter coefficients. This is an aid to determine if the lookahead and\nlookback sizes L1 and L2 are set sufficiently large.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Determine lookback\nL1 = K2\n# Determine lookahead\nL2 = K2\ndigital_estimator_fir = cbadc.digital_estimator.FIRFilter(\n    analog_system, digital_control2, eta2, L1, L2)\n\nprint(digital_estimator_fir, \"\\n\")\n\ndigital_estimator_fir(simulator2)\n\n# Next visualize the decay of the resulting filter coefficients.\nh_index = np.arange(-L1, L2)\n\nimpulse_response = np.abs(np.array(digital_estimator_fir.h[0, :, :])) ** 2\nimpulse_response_dB = 10 * np.log10(impulse_response)\n\nfig, ax = plt.subplots(2)\nfor index in range(N):\n    ax[0].plot(h_index, impulse_response[:, index],\n               label=f\"$h_{index + 1}[k]$\")\n    ax[1].plot(h_index, impulse_response_dB[:, index],\n               label=f\"$h_{index + 1}[k]$\")\nax[0].legend()\nfig.suptitle(f\"For $\\eta^2 = {10 * np.log10(eta2)}$ [dB]\")\nax[1].set_xlabel(\"filter taps k\")\nax[0].set_ylabel(\"$| h_\\ell [k]|^2_2$\")\nax[1].set_ylabel(\"$| h_\\ell [k]|^2_2$ [dB]\")\nax[0].set_xlim((-50, 50))\nax[0].grid(which='both')\nax[1].set_xlim((-50, 500))\nax[1].set_ylim((-200, 0))\nax[1].grid(which='both')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## IIR Filter Estimator\n\nThe IIR filter is closely related to the FIR filter with the exception\nof an moving average computation.\nSee :py:class:`cbadc.digital_estimator.IIRFilter` for more information.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Determine lookahead\nL2 = K2\n\ndigital_estimator_iir = cbadc.digital_estimator.IIRFilter(\n    analog_system, digital_control3, eta2, L2)\n\nprint(digital_estimator_iir, \"\\n\")\n\ndigital_estimator_iir(simulator3)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Parallel Estimator\n\nNext we instantiate the parallel estimator\n:py:class:`cbadc.digital_estimator.ParallelEstimator`. The parallel estimator\nresembles the default estimator but diagonalizes the filter coefficients\nresulting in a more computationally more efficient filter that can be\nparallelized into independent filter operations.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Instantiate the digital estimator (this is where the filter coefficients are\n# computed).\ndigital_estimator_parallel = cbadc.digital_estimator.ParallelEstimator(\n    analog_system, digital_control4, eta2, K1, K2)\n\ndigital_estimator_parallel(simulator4)\nprint(digital_estimator_parallel, \"\\n\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Estimating (Filtering)\n\nNext we execute all simulation and estimation tasks by iterating over the\nestimators. Note that since no stop criteria is set for either the analog\nsignal, the simulator, or the digital estimator this iteration could\npotentially continue until the default stop criteria of 2^63 iterations.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Set simulation length\nsize = K2 << 4\nu_hat_batch = np.zeros(size)\nu_hat_fir = np.zeros(size)\nu_hat_iir = np.zeros(size)\nu_hat_parallel = np.zeros(size)\nfor index in range(size):\n    u_hat_batch[index] = next(digital_estimator_batch)\n    u_hat_fir[index] = next(digital_estimator_fir)\n    u_hat_iir[index] = next(digital_estimator_iir)\n    u_hat_parallel[index] = next(digital_estimator_parallel)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Visualizing Results\n\nFinally, we summarize the comparision by visualizing the resulting estimate\nin both time and frequency domain.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "t = np.arange(size)\n# compensate the built in L1 delay of FIR filter.\nt_fir = np.arange(-L1 + 1, size - L1 + 1)\nt_iir = np.arange(-L1 + 1, size - L1 + 1)\nu = np.zeros_like(u_hat_batch)\nfor index, tt in enumerate(t):\n    u[index] = analog_signal.evaluate(tt * T)\nplt.plot(t, u_hat_batch, label=\"$\\hat{u}(t)$ Batch\")\nplt.plot(t_fir, u_hat_fir, label=\"$\\hat{u}(t)$ FIR\")\nplt.plot(t_iir, u_hat_iir, label=\"$\\hat{u}(t)$ IIR\")\nplt.plot(t, u_hat_parallel, label=\"$\\hat{u}(t)$ Parallel\")\nplt.plot(t, stf_at_omega * u, label=\"$\\mathrm{STF}(2 \\pi f_u) * u(t)$\")\nplt.xlabel('$t / T$')\nplt.legend()\nplt.title(\"Estimated input signal\")\nplt.grid(which='both')\nplt.xlim((-100, 500))\nplt.tight_layout()\n\nplt.figure()\nplt.plot(t, u_hat_batch, label=\"$\\hat{u}(t)$ Batch\")\nplt.plot(t_fir, u_hat_fir, label=\"$\\hat{u}(t)$ FIR\")\nplt.plot(t_iir, u_hat_iir, label=\"$\\hat{u}(t)$ IIR\")\nplt.plot(t, u_hat_parallel, label=\"$\\hat{u}(t)$ Parallel\")\nplt.plot(t, stf_at_omega * u, label=\"$\\mathrm{STF}(2 \\pi f_u) * u(t)$\")\nplt.xlabel('$t / T$')\nplt.legend()\nplt.title(\"Estimated input signal\")\nplt.grid(which='both')\nplt.xlim((t_fir[-1] - 50, t_fir[-1]))\nplt.tight_layout()\n\nplt.figure()\nplt.plot(t, u_hat_batch, label=\"$\\hat{u}(t)$ Batch\")\nplt.plot(t_fir, u_hat_fir, label=\"$\\hat{u}(t)$ FIR\")\nplt.plot(t_iir, u_hat_iir, label=\"$\\hat{u}(t)$ IIR\")\nplt.plot(t, u_hat_parallel, label=\"$\\hat{u}(t)$ Parallel\")\nplt.plot(t, stf_at_omega * u, label=\"$\\mathrm{STF}(2 \\pi f_u) * u(t)$\")\nplt.xlabel('$t / T$')\nplt.legend()\nplt.title(\"Estimated input signal\")\nplt.grid(which='both')\n# plt.xlim((t_fir[0], t[-1]))\nplt.xlim(((1 << 14) - 100, (1 << 14) + 100))\nplt.tight_layout()\n\nbatch_error = stf_at_omega * u - u_hat_batch\nfir_error = stf_at_omega * u[:(u.size - L1 + 1)] - u_hat_fir[(L1 - 1):]\niir_error = stf_at_omega * u[:(u.size - L1 + 1)] - u_hat_iir[(L1 - 1):]\nparallel_error = stf_at_omega * u - u_hat_parallel\nplt.figure()\nplt.plot(t, batch_error,\n         label=\"$|\\mathrm{STF}(2 \\pi f_u) * u(t) - \\hat{u}(t)|$ Batch\")\nplt.plot(t[:(u.size - L1 + 1)], fir_error,\n         label=\"$|\\mathrm{STF}(2 \\pi f_u) * u(t) - \\hat{u}(t)|$ FIR\")\nplt.plot(t[:(u.size - L1 + 1)], iir_error,\n         label=\"$|\\mathrm{STF}(2 \\pi f_u) * u(t) - \\hat{u}(t)|$ IIR\")\nplt.plot(t, parallel_error,\n         label=\"$|\\mathrm{STF}(2 \\pi f_u) * u(t) - \\hat{u}(t)|$ Parallel\")\nplt.xlabel('$t / T$')\nplt.xlim(((1 << 14) - 100, (1 << 14) + 100))\nplt.ylim((-0.00001, 0.00001))\nplt.legend()\nplt.title(\"Estimation error\")\nplt.grid(which='both')\nplt.tight_layout()\n\n\nprint(f\"Average Batch Error: {np.linalg.norm(batch_error) / batch_error.size}\")\nprint(f\"Average FIR Error: {np.linalg.norm(fir_error) / fir_error.size}\")\nprint(f\"Average IIR Error: {np.linalg.norm(iir_error) / iir_error.size}\")\nprint(\n    f\"\"\"Average Parallel Error: { np.linalg.norm(parallel_error)/\n    parallel_error.size}\"\"\")\n\nplt.figure()\nu_hat_batch_clipped = u_hat_batch[(K1 + K2):-K2]\nu_hat_fir_clipped = u_hat_fir[(L1 + L2):]\nu_hat_iir_clipped = u_hat_iir[(K1 + K2):-K2]\nu_hat_parallel_clipped = u_hat_parallel[(K1 + K2):-K2]\nu_clipped = stf_at_omega * u\nf_batch, psd_batch = cbadc.utilities.compute_power_spectral_density(\n    u_hat_batch_clipped)\nf_fir, psd_fir = cbadc.utilities.compute_power_spectral_density(\n    u_hat_fir_clipped)\nf_iir, psd_iir = cbadc.utilities.compute_power_spectral_density(\n    u_hat_iir_clipped)\nf_parallel, psd_parallel = cbadc.utilities.compute_power_spectral_density(\n    u_hat_parallel_clipped)\nf_ref, psd_ref = cbadc.utilities.compute_power_spectral_density(u_clipped)\nplt.semilogx(f_ref, 10 * np.log10(psd_ref),\n             label=\"$\\mathrm{STF}(2 \\pi f_u) * U(f)$\")\nplt.semilogx(f_batch, 10 * np.log10(psd_batch), label=\"$\\hat{U}(f)$ Batch\")\nplt.semilogx(f_fir, 10 * np.log10(psd_fir), label=\"$\\hat{U}(f)$ FIR\")\nplt.semilogx(f_iir, 10 * np.log10(psd_iir), label=\"$\\hat{U}(f)$ IIR\")\nplt.semilogx(f_parallel, 10 * np.log10(psd_parallel),\n             label=\"$\\hat{U}(f)$ Parallel\")\nplt.legend()\nplt.ylim((-200, 50))\nplt.xlim((f_fir[1], f_fir[-1]))\nplt.xlabel('frequency [Hz]')\nplt.ylabel('$ \\mathrm{V}^2 \\, / \\, (1 \\mathrm{Hz})$')\nplt.grid(which='both')\nplt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Compute Time\n\nCompare the execution time of each estimator\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "def dummy_input_control_signal():\n    while True:\n        yield np.zeros(M, dtype=np.int8)\n\n\ndef iterate_number_of_times(iterator, number_of_times):\n    for _ in range(number_of_times):\n        _ = next(iterator)\n\n\ndigital_estimator_batch = cbadc.digital_estimator.DigitalEstimator(\n    analog_system,\n    digital_control1,\n    eta2,\n    K1,\n    K2)\ndigital_estimator_fir = cbadc.digital_estimator.FIRFilter(\n    analog_system,\n    digital_control2,\n    eta2,\n    L1,\n    L2)\ndigital_estimator_parallel = cbadc.digital_estimator.ParallelEstimator(\n    analog_system,\n    digital_control4,\n    eta2,\n    K1,\n    K2)\ndigital_estimator_iir = cbadc.digital_estimator.IIRFilter(\n    analog_system,\n    digital_control3,\n    eta2,\n    L2)\n\ndigital_estimator_batch(dummy_input_control_signal())\ndigital_estimator_fir(dummy_input_control_signal())\ndigital_estimator_parallel(dummy_input_control_signal())\ndigital_estimator_iir(dummy_input_control_signal())\n\nlength = 1 << 14\nrepetitions = 10\n\nprint(\"Digital Estimator:\")\nprint(timeit.timeit(lambda: iterate_number_of_times(\n    digital_estimator_batch, length), number=repetitions), 'sec \\n')\n\nprint(\"FIR Estimator:\")\nprint(timeit.timeit(lambda: iterate_number_of_times(\n    digital_estimator_fir, length), number=repetitions), 'sec \\n')\n\nprint(\"IIR Estimator:\")\nprint(timeit.timeit(lambda: iterate_number_of_times(\n    digital_estimator_iir, length), number=repetitions), 'sec \\n')\n\nprint(\"Parallel Estimator:\")\nprint(timeit.timeit(lambda: iterate_number_of_times(\n    digital_estimator_parallel, length), number=repetitions), 'sec \\n')"
      ]
    }
  ],
  "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.8.5"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}