Try rounding.

libraries/base/common/python/try_round.py

+#! /usr/bin/env python3
+###############################################################################
+#
+# Copyright (C) 2021
+# ASTRON (Netherlands Institute for Radio Astronomy) <http://www.astron.nl/>
+# P.O.Box 2, 7990 AA Dwingeloo, The Netherlands
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program.  If not, see <http://www.gnu.org/licenses/>.
+#
+###############################################################################
+
+# Author: Eric Kooistra
+# Date: may 2021
+# Purpose:
+#   Simulate linearity of rounding
+# Description:
+#   Usage:
+#   > python3 try_round.py -N 1024 -ampl 12 --useplot
+
+
+import argparse
+import numpy as np
+import matplotlib
+matplotlib.use('tkagg')  # to make X11 forwarding work
+import matplotlib.pylab as plt
+
+"""Try histogram
+
+Usage:
+> python try_histogram.py -h
+
+"""
+
+import argparse
+import numpy as np
+import matplotlib.pylab as plt
+import common as cm
+
+figNr = 1
+
+# Parse arguments to derive user parameters
+_parser = argparse.ArgumentParser('try_histogram')
+_parser.add_argument('-N', default=1024, type=int, help='Number of points of FFT')
+_parser.add_argument('-w', default=8, type=int, help='Total number of bits')
+_parser.add_argument('-f', default=0, type=int, help='Number of bits of fraction that gets rounded')
+_parser.add_argument('--useplot', action='store_true', dest='useplot', default=False, help='Default without plotting, else with plotting')
+args = _parser.parse_args()
+
+N = args.N
+width = args.w
+ampl = 2**(width-1)
+fraction = args.f
+scale = 2**fraction
+useplot = args.useplot
+
+# Sinus
+t = 2 * np.pi * np.arange(0, N) / N
+a_adc = np.round(ampl*np.sin(t))
+a_away = [cm.int_round(inp = x, r = fraction, direction="HALF_AWAY") for x in a_adc]
+a_even = [cm.int_round(inp = x, r = fraction, direction="HALF_EVEN") for x in a_adc]
+a_adc = a_adc / scale
+
+# Spectrum
+# . use rfft for real input instead of fft, to only calculate the >=0 frequencies of the FFT
+f = np.arange(0, N/2 + 1)
+A_adc = np.power(np.abs(np.fft.rfft(a_adc)), 2)
+A_away = np.power(np.abs(np.fft.rfft(a_away)), 2)
+A_even = np.power(np.abs(np.fft.rfft(a_even)), 2)
+
+A_adc_cw = A_adc[1]  # bin 1 with CW
+A_away_cw = A_away[1]  # bin 1 with CW
+A_even_cw = A_even[1]  # bin 1 with CW
+A_adc_spurious = A_adc[0] + np.sum(A_adc[2:])  # skip bin 1 with CW
+A_away_spurious = A_away[0] + np.sum(A_away[2:])  # skip bin 1 with CW
+A_even_spurious = A_even[0] + np.sum(A_even[2:])  # skip bin 1 with CW
+SNR_adc = 10*np.log10(A_adc_cw / A_adc_spurious)
+SNR_adc_scale = SNR_adc - fraction * 20 * np.log10(2)
+SNR_away = 10*np.log10(A_away_cw / A_away_spurious)
+SNR_even = 10*np.log10(A_even_cw / A_even_spurious)
+
+print('SNR_adc       = %7.3f dB' % SNR_adc)
+print('SNR_adc_scale = %7.3f dB' % SNR_adc_scale)
+print('SNR_away      = %7.3f dB' % SNR_away)
+print('SNR_even      = %7.3f dB' % SNR_even)
+
+if useplot:
+    # Plot sinus
+    plt.figure(figNr)
+    plt.plot(t, a_adc, '-r', t, a_away, '.-g', t, a_even, '.b')
+    plt.xlabel('Time')
+    plt.ylabel('Voltage')
+    plt.title('Sinus (%d values)' % N)
+    plt.grid(True)
+    figNr += 1
+
+    # Plot spectrum
+    plt.figure(figNr)
+    plt.plot(f, A_adc, '-r', f, A_away, '.-g', f, A_even, '.b')
+    plt.xlabel('Frequency')
+    plt.ylabel('Power')
+    plt.title('Frequency (%d bins)' % (N/2))
+    plt.grid(True)
+    figNr += 1
+
+    # Show plots
+    plt.show()