From b0b6e3d09509a621a5c97f40d7c736af8a03142b Mon Sep 17 00:00:00 2001
From: Eric Kooistra <kooistra@astron.nl>
Date: Thu, 18 Nov 2021 14:50:45 +0100
Subject: [PATCH] Try rounding.
---
libraries/base/common/python/try_round.py | 117 ++++++++++++++++++++++
1 file changed, 117 insertions(+)
create mode 100755 libraries/base/common/python/try_round.py
diff --git a/libraries/base/common/python/try_round.py b/libraries/base/common/python/try_round.py
new file mode 100755
index 0000000000..bf284b5d8c
--- /dev/null
+++ b/libraries/base/common/python/try_round.py
@@ -0,0 +1,117 @@
+#! /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()
--
GitLab