Add -f and -r to investigate effect of resolution.

libraries/base/common/python/try_round_weight.py

@@ -34,10 +34,15 @@
 #   . increasing -N improves the results, for LOFAR subbands N = 195312
 #   . it may be preferred to apply the subband weights to the unquantized
 #     WPFB output.
-#   Usage:
+# Note:
+# . For values exactly halfway between rounded decimal values, NumPy rounds to
+#   the nearest even value. Thus 1.5 and 2.5 round to 2.0, -0.5 and 0.5 round
+#   to 0.0, etc.
+# Usage:
 #   > python3 try_round_weight.py -N 195312
 
 import argparse
+import textwrap
 
 import numpy as np
 import matplotlib
@@ -47,9 +52,33 @@ import matplotlib.pyplot as plt
 import common as cm
 
 # Parse arguments to derive user parameters
-_parser = argparse.ArgumentParser('try_round_weight')
+_parser = argparse.ArgumentParser(
+    description="".join(textwrap.dedent("""\
+        Model effect of applying a weight after or before rounding:
+
+        . weights in range(w_lo, w_hi, w_step)
+        . sigma of noise in range(s_lo, s_hi, s_step)
+        . use N = 195312 to model number of subband periods in 1 s of LOFAR SST
+        . use different seed S to check whether result in plots depends on seed
+        . use r = 0 for default integer rounding, use r > 0 to have extra
+          LSbits resolution before rounding
+        . use f = 0 for full double resolution before rounding, use f > 0 and
+          r = 0 to have noise with a fraction of f LSbits before rounding
+
+        # Get an overview
+        > python try_round_weight.py --w_lo 0.2 --w_hi 3.0 --w_step 0.01 --s_lo 0.5 --s_hi 10 --s_step 0.2 -N 195312 -S 1
+
+        # Zoom in at w = 0.75
+        > python try_round_weight.py --w_lo 0.7 --w_hi 0.8 --w_step 0.0001 --s_lo 1 --s_hi 10 --s_step 1 -N 195312 -S 0
+
+        # Use -r = 6 to see effect of having more resolution before rounding
+        > python try_round_weight.py --w_lo 0.7 --w_hi 0.8 --w_step 0.0001 --s_lo 1 --s_hi 10 --s_step 1 -N 195312 -S 0 -r 6
+        \n""")),
+    formatter_class=argparse.RawTextHelpFormatter)
 _parser.add_argument('-S', default=0, type=int, help='Random number seed')
 _parser.add_argument('-N', default=1000, type=int, help='Number of input samples')
+_parser.add_argument('-f', default=0, type=int, help='Number of LSbits in fraction, use default 0 for double accuracy')
+_parser.add_argument('-r', default=0, type=int, help='Number of LSbits extra resolution before rounding')
 _parser.add_argument('--s_lo', default=0.1, type=float, help='Lowest sigma')
 _parser.add_argument('--s_hi', default=25.0, type=float, help='Highest sigma')
 _parser.add_argument('--s_step', default=0.1, type=float, help='Step sigma')
@@ -63,7 +92,6 @@ np.random.seed(args.S)
 # Prepare noise signal
 N_samples = args.N
 noise = np.random.randn(N_samples)
-noise = np.random.randn(N_samples)
 noise /= np.std(noise)
 
 # Noise levels range, 1 unit = 1 LSbit
@@ -80,13 +108,32 @@ weight_step = args.w_step
 weights = np.arange(weight_lo, weight_hi, weight_step)
 N_weights = len(weights)
 
+# Fraction
+fraction = args.f
+fraction_factor = 2**fraction
+
+# Resolution
+resolution = args.r
+resolution_factor = 2**resolution
+
 # Determine weighted rounded noise sigma / weighted noise sigma for range of weights and input noise sigmas
 sigmas_ratio = np.nan * np.zeros((N_weights, N_sigmas))  # w rows, s cols
 sigmas_qq = np.zeros((N_weights, N_sigmas))
 sigmas_sq = np.zeros((N_weights, N_sigmas))
 for s, sigma in enumerate(sigmas):
    noise_s = noise * sigma
-   noise_q = np.round(noise_s)
+   if resolution == 0:
+       if fraction == 0:
+           # Default use full double resolution of fraction of noise_s
+           noise_q = np.round(noise_s)
+       else:
+           # First use round() to get noise_f with a fraction of fraction number of LSbits
+           noise_f = np.round(noise_s * fraction_factor) / fraction_factor
+           # Then use round to round the fraction of fraction number of LSbits in noise_f
+           noise_q = np.round(noise_f)
+   else:
+       noise_q = np.round(noise_s * resolution_factor) / resolution_factor
+
    for w, weight in enumerate(weights):
        noise_q_weighted_q = np.round(noise_q * weight)  # apply weight to rounded noise
        noise_s_weighted_q = np.round(noise_s * weight)  # apply weight to original noise
@@ -110,7 +157,7 @@ for s, sigma in enumerate(sigmas):
     plt.plot(weights, sigmas_qq_T[s], label='s = %4.2f' % sigma)
 plt.title("Sigma of weighted quantized noise")
 plt.xlabel("Weight")
-plt.ylabel("Sigma")
+plt.ylabel("Sigma_qq")
 plt.legend(loc='upper right')
 plt.grid()
 
@@ -130,7 +177,7 @@ for w, weight in enumerate(weights):
     plt.plot(sigmas, sigmas_ratio[w], label='w = %4.2f' % weight)
 plt.title("Relative sigma difference of weighting after / before quantisation")
 plt.xlabel("Sigma")
-plt.ylabel("Relative sigma difference")
+plt.ylabel("Relative sigma difference (s_qq / s_sq)")
 plt.legend(loc='upper right')
 plt.grid()