.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "tutorial/plot_2_simulation.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. .. rst-class:: sphx-glr-example-title .. _sphx_glr_tutorial_plot_2_simulation.py: ========== Simulation ========== Finite-wordlength analysis is a critical step in the design of algorithms. Unfortunately, is is often overlooked or oversimplified due to the complexity of the analysis. With B-ASIC, a correct analysis becomes easier. SFGs can be simulated in arbitrary precision using :class:`~b_asic.simulation.Simulation` and :class:`~b_asic.data_type.DataType`. This tutorial showcases this through a finite-wordlength analysis of a wave digital filter. .. GENERATED FROM PYTHON SOURCE LINES 18-26 Algorithm --------- We begin by designing an elliptic low-pass filter with the following specifications: * Order: 7 * Passband edge frequency: 0.3 times the Nyquist frequency * Passband ripple: 0.1 dB * Stopband attenuation: 60 dB .. GENERATED FROM PYTHON SOURCE LINES 26-36 .. code-block:: Python from scipy.signal import iirfilter from b_asic.wdf import lattice_coeffs_from_tf b, a = iirfilter(N=7, Wn=0.3, rp=0.1, rs=60, btype="low", ftype="ellip") print("Transfer function coefficients:") print("Numerator:", b) print("Denominator:", a) .. rst-class:: sphx-glr-script-out .. code-block:: none Transfer function coefficients: Numerator: [0.00698396 0.00219871 0.0142081 0.01065808 0.01065808 0.0142081 0.00219871 0.00698396] Denominator: [ 1. -3.97699108 7.99282358 -9.90310426 8.08290621 -4.30338914 1.38344285 -0.20759045] .. GENERATED FROM PYTHON SOURCE LINES 37-39 We then use the module :mod:`b_asic.wdf` to derive adaptor coefficients for a lattice structure that implements the filter. .. GENERATED FROM PYTHON SOURCE LINES 39-44 .. code-block:: Python adaptor_coeffs = lattice_coeffs_from_tf(a) print("Lattice adaptor coefficients:") for i, coeff in enumerate(adaptor_coeffs): print(f" a{i}: {coeff}") .. rst-class:: sphx-glr-script-out .. code-block:: none Lattice adaptor coefficients: a0: 0.6072682339178634 a1: -0.9233304755403088 a2: 0.5683222541662454 a3: -0.7388670919448468 a4: 0.635533951021361 a5: -0.5010756697748922 a6: 0.7804684909217103 .. GENERATED FROM PYTHON SOURCE LINES 45-46 Now, use an SFG generator to construct the SFG of the lattice wave digital filter. .. GENERATED FROM PYTHON SOURCE LINES 46-51 .. code-block:: Python from b_asic.sfg_generators.wave_digital_filters import lattice_wdf wdf_sfg = lattice_wdf(adaptor_coeffs) wdf_sfg .. image-sg:: /tutorial/images/sphx_glr_plot_2_simulation_001.png :alt: plot 2 simulation :srcset: /tutorial/images/sphx_glr_plot_2_simulation_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 52-63 Here, each adaptor is represented as a single operation with the following equations: .. math:: y_0 & = x_1 + \alpha\times\left(x_1 - x_0\right)\\ y_1 & = x_0 + \alpha\times\left(x_1 - x_0\right) where :math:`\alpha` is the adaptor coefficient. The coefficients for the adaptors are now represented as 64-bit floating point numbers. Multiplication with such a coefficient will be extremely expensive in terms of hardware resources. As such, we want to quantize the coefficients, as much as possible. But before we do that, we need to investigate the need for scaling the coefficients. .. GENERATED FROM PYTHON SOURCE LINES 65-71 Coefficient Quantization ------------------------ Coefficient quantization introduces a static error to the transfer function. In order to investigate the effects of quantization, we quantize the coefficients to 5 to 11 fractional bits, and derive the transfer functions of the resulting filters. .. GENERATED FROM PYTHON SOURCE LINES 71-94 .. code-block:: Python import matplotlib.pyplot as plt import numpy as np from scipy.signal import chirp, freqz from b_asic.quantization import quantize from b_asic.signal_generator import Impulse from b_asic.simulation import Simulation w_ref, H_ref = freqz(b, a) impulse_responses = {} sfgs = {} for n_bits in range(8, 16): q_coeffs = [quantize(c, fractional_bits=n_bits) for c in adaptor_coeffs] sfg = wdf_sfg.copy() for i in range(len(q_coeffs)): sfg.find_by_name(f"a{i}")[0].value = q_coeffs[i] sfgs[n_bits] = sfg sim = Simulation(sfg, [Impulse()]) sim.run_for(1_000) impulse_responses[n_bits] = np.array([float(v) for v in sim.results["out0"]]) .. image-sg:: /tutorial/images/sphx_glr_plot_2_simulation_002.png :alt: plot 2 simulation :srcset: /tutorial/images/sphx_glr_plot_2_simulation_002.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 95-96 Stopband attenuation plotted for different wordlengths. .. GENERATED FROM PYTHON SOURCE LINES 96-107 .. code-block:: Python from mplsignal import freqz_fir fig, ax = plt.subplots() for n_bits, h in impulse_responses.items(): freqz_fir(h, ax=ax, style="magnitude", label=f"{n_bits} bits") ax.legend() ax.set_xlim(0.3 * np.pi) ax.set_ylim(-80, -40) ax.grid(True) .. image-sg:: /tutorial/images/sphx_glr_plot_2_simulation_003.png :alt: plot 2 simulation :srcset: /tutorial/images/sphx_glr_plot_2_simulation_003.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 108-109 Passband attenuation plotted for different wordlengths. .. GENERATED FROM PYTHON SOURCE LINES 109-118 .. code-block:: Python fig, ax = plt.subplots() for n_bits, h in impulse_responses.items(): freqz_fir(h, ax=ax, style="magnitude", label=f"{n_bits} bits") ax.legend() ax.set_xlim(0, 0.32 * np.pi) ax.set_ylim(-0.15, 0.05) ax.grid(True) .. image-sg:: /tutorial/images/sphx_glr_plot_2_simulation_004.png :alt: plot 2 simulation :srcset: /tutorial/images/sphx_glr_plot_2_simulation_004.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 119-123 As a conclusion, one can see that coefficient quantization leeds to a similar but different algorithm. We can now choose a number of bits that we deem to meet our specification. In this case, let us continue with 14 bits as this meets the stopband attenuation requirement of 60 dB. .. GENERATED FROM PYTHON SOURCE LINES 123-125 .. code-block:: Python q_sfg = sfgs[14] .. GENERATED FROM PYTHON SOURCE LINES 126-132 Scaling ------- As implementation will be done in fixed-point arithmetic, thus, we need to investigate the need for scaling to maximize the dynamic range of the filter and protect against overflow. To investigate this, we can once again render the SFG, this time with the l2-norm values of all signals. .. GENERATED FROM PYTHON SOURCE LINES 132-135 .. code-block:: Python from b_asic.core_operations import LeftShift, RightShift q_sfg.sfg_digraph(signal_info="l1-norm") .. image-sg:: /tutorial/images/sphx_glr_plot_2_simulation_005.png :alt: plot 2 simulation :srcset: /tutorial/images/sphx_glr_plot_2_simulation_005.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 136-140 From the figure, we see that the l1-norm exceeds 1 at a lot of places, meaining that overflow can occur if we do not scale the input, which we want to prevent in this example. Note that one cannot scale inside loops... .. GENERATED FROM PYTHON SOURCE LINES 142-143 Inserting appropriate shifts to safe-scale. .. GENERATED FROM PYTHON SOURCE LINES 143-151 .. code-block:: Python q_sfg = q_sfg.insert_operation_before("a0", RightShift(2), 0) q_sfg = q_sfg.insert_operation_before("a1", RightShift(5), 0) q_sfg = q_sfg.insert_operation_after("a0", RightShift(2), 0) q_sfg = q_sfg.insert_operation_after("a3", LeftShift(1), 0) q_sfg = q_sfg.insert_operation_after("a5", LeftShift(2), 0) q_sfg = q_sfg.remove_operation("cmul0") q_sfg.sfg_digraph(signal_info="l1-norm") .. image-sg:: /tutorial/images/sphx_glr_plot_2_simulation_006.png :alt: plot 2 simulation :srcset: /tutorial/images/sphx_glr_plot_2_simulation_006.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 152-154 Note that the output is scaled by 1/4 now. However, it is not a problem since its just a matter of selecting the correct bits. .. GENERATED FROM PYTHON SOURCE LINES 156-162 SNR Analysis ------------ SNR is calculated as the ratio of the signal power to the noise power, where the noise is the difference between simulated outputs with and without quantization. This is calculated for a range of wordlengths along with magnitude truncation for data quantization, as it protects against limit-cycles while still being simple to implement. .. GENERATED FROM PYTHON SOURCE LINES 162-194 .. code-block:: Python from b_asic.data_type import DataType from b_asic.quantization import QuantizationMode N = 1_000 rng = np.random.default_rng(0) wl_range = range(8, 19) snr_values = [] for WL in wl_range: data = rng.integers(-2**(WL - 1), 2**(WL - 1), size=(1, N)) * 2**(-(WL - 1)) dt = DataType( (1, WL), quantization_mode=QuantizationMode.MAGNITUDE_TRUNCATION ) sim_quant = Simulation(q_sfg, data, data_type=dt) sim_quant.run_for(N) sim_ref = Simulation(q_sfg, data) sim_ref.run_for(N) out_ref = np.array([float(v) for v in sim_ref.results["out0"]]) out_quant = np.array([float(v) for v in sim_quant.results["out0"]]) P_signal = np.mean(out_ref**2) P_noise = np.mean((out_quant - out_ref) ** 2) snr_values.append(10 * np.log10(P_signal / P_noise)) fig, ax = plt.subplots() ax.plot(wl_range, snr_values, marker="o") ax.set_xlabel("Fractional bits") ax.set_ylabel("SNR, dB") ax.grid(True) .. image-sg:: /tutorial/images/sphx_glr_plot_2_simulation_007.png :alt: plot 2 simulation :srcset: /tutorial/images/sphx_glr_plot_2_simulation_007.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 195-202 As expected from theory, the SNR increases linearly with the number of bits, with a slope of approximately 6 dB per bit. Here, one should simply pick a number of bits that meets the SNR requirements posed by the application. Note that if we had more information about the input data, there may exist a scaling that yields a better SNR. Since we have no such information, we safe-scale. .. GENERATED FROM PYTHON SOURCE LINES 204-210 Conclusion ---------- In this tutorial, we have seen how B-ASIC can be used to perform finite-wordlength analysis of algorithms targeting arbitrary precision. We are now ready to move on to the next step of the design process. .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 11.570 seconds) .. _sphx_glr_download_tutorial_plot_2_simulation.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: plot_2_simulation.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_2_simulation.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: plot_2_simulation.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_