Merge branch 'master' into L2SDP-836

libraries/base/common/python/try_round.py

@@ -40,6 +40,10 @@ import matplotlib.pylab as plt
 Usage:
 > python try_round.py -h
 
+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.
 """
 
 import argparse

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,55 +52,156 @@ 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
+        > python try_round_weight.py --w_lo 0.2 --w_hi 2.0 --w_step 0.01 --s_lo 1 --s_hi 10 --s_step 1 -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('--weight_lo', default=0.3, type=float, help='Lowest weight')
-_parser.add_argument('--weight_hi', default=2.0, type=float, help='Highest weight')
-_parser.add_argument('--weight_step', default=0.1, type=float, help='Step weight')
+_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')
+_parser.add_argument('--w_lo', default=0.3, type=float, help='Lowest weight')
+_parser.add_argument('--w_hi', default=2.0, type=float, help='Highest weight')
+_parser.add_argument('--w_step', default=0.1, type=float, help='Step weight')
 args = _parser.parse_args()
 
-N_samples = args.N
-weight_lo = args.weight_lo
-weight_hi = args.weight_hi
-weight_step = args.weight_step
+np.random.seed(args.S)
 
-# Prepare signals
+# Prepare noise signal
+N_samples = args.N
 noise = np.random.randn(N_samples)
 noise /= np.std(noise)
 
-# Noise level range, 1 unit = 1 LSbit
-sigma_lo = 0.1
-sigma_hi = 25
-sigmas = np.arange(sigma_lo, sigma_hi, 0.1)
+# Noise levels range, 1 unit = 1 LSbit
+sigma_lo = args.s_lo
+sigma_hi = args.s_hi
+sigma_step = args.s_step
+sigmas = np.arange(sigma_lo, sigma_hi, sigma_step)
 N_sigmas = len(sigmas)
 
-# Weight range, unit weight = 1
+# Weights range, unit weight = 1
+weight_lo = args.w_lo
+weight_hi = args.w_hi
+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
-       sigma_qq = np.std(noise_q_weighted_q)
-       sigma_sq = np.std(noise_s_weighted_q)
-       if sigma_sq != 0:
-           sigmas_ratio[w][s] = sigma_qq / sigma_sq  # weighted rounded noise sigma / weighted noise sigma
+       s_qq = np.std(noise_q_weighted_q)
+       s_sq = np.std(noise_s_weighted_q)
+       sigmas_qq[w][s] = s_qq
+       sigmas_sq[w][s] = s_sq
+       if s_sq != 0:
+           sigmas_ratio[w][s] = s_qq / s_sq  # weighted rounded noise sigma / weighted noise sigma
+# Transpose [w][s] to have index ranges [s][w]
+sigmas_ratio_T = sigmas_ratio.transpose()
+sigmas_qq_T = sigmas_qq.transpose()
+sigmas_sq_T = sigmas_sq.transpose()
 
 # Plot results
 figNr = 0
 
+figNr += 1
+plt.figure(figNr)
+for s, sigma in enumerate(sigmas):
+    # Plot sigma_qq of twice quantized noise as function of weight for
+    # different input 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_qq")
+plt.legend(loc='upper right')
+plt.grid()
+
+figNr += 1
+plt.figure(figNr)
+for s, sigma in enumerate(sigmas):
+    # Plot sigma_qq of twice quantized noise as function of weight for
+    # different input sigmas.
+    # Normalize the sigma_qq by the weight, so that it can be compared with
+    # the input sigma that is shown by the horizontal sigma reference lines.
+    plt.plot(weights, sigmas_qq_T[s] / weights, label='s = %4.2f' % sigma)
+    plt.plot(weights, sigmas[s]*np.ones(N_weights))  # add sigma reference lines
+plt.title("Sigma of weighted quantized noise, normalized for weight")
+plt.xlabel("Weight")
+plt.ylabel("Sigma_qq")
+plt.legend(loc='upper right')
+plt.grid()
+
+figNr += 1
+plt.figure(figNr)
+for s, sigma in enumerate(sigmas):
+    # Plot ratio of sigma_qq / sigma_sq as function of weight for different
+    # input sigma. The ratio deviation from 1 tells how much the twice
+    # quantized noise deviates from the noise that is only quantized after
+    # the weighting.
+    plt.plot(weights, sigmas_ratio_T[s], label='s = %4.2f' % sigma)
+plt.title("Relative sigma difference of weighting after / before quantisation")
+plt.xlabel("Weight")
+plt.ylabel("Relative sigma difference")
+plt.legend(loc='upper right')
+plt.grid()
+
 figNr += 1
 plt.figure(figNr)
 for w, weight in enumerate(weights):
+    # Plot ratio of sigma_qq / sigma_sq as function of input sigma for
+    # different 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()