Added conclusion for sufficient -r.

libraries/base/common/python/try_round_weight.py

@@ -28,16 +28,27 @@
 #   . quantized subbands --> sigma_qq
 #   . unquantized subbands --> sigma_sq
 #   Preliminary conclusion:
-#   . for small input noise with sigma < 2 the output sigma gets disturbed
+#   . For small input noise with sigma < 2 the output sigma gets disturbed
 #     due to the weighting if the weighting is applied after the subband
 #     quantisation
-#   . increasing -N improves the results, for LOFAR subbands N = 195312
+#   . Increasing -N improves the results, for LOFAR subbands N = 195312
 #   . it may be preferred to apply the subband weights to the unquantized
 #     WPFB output.
+#   . Choosing sufficient intermediate resolution (-r) before applying
+#     weights:
+#     - Rounding noise changes the sigma of the noise anyway, as shown by
+#       sigmas_ratio_sq_input_T. Therefore choose resolution such that
+#       jumps in sigma due to rounding weighted noise, as shown by
+#       sigmas_ratio_qq_sq_T, is much smaller than sigmas_ratio_sq_input_T
+#       ==> -r 2.
+#     - For input sigma >= 1
+#       . Choose sigmas_ratio_qq_sq_T < 10%   ==> -r 3
+#       . Choose sigmas_ratio_qq_sq_T <  1%   ==> -r 6
+#       . Choose sigmas_ratio_qq_sq_T <  0.1% ==> -r 9
 # 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.
+# . For values exactly halfway between rounded decimal values, NumPy of
+#   Python3 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. (Python2 rounds half away from zero).
 # Usage:
 #   > python3 try_round_weight.py -N 195312
 
@@ -74,6 +85,7 @@ _parser = argparse.ArgumentParser(
 
         # 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
+        > python try_round_weight.py --w_lo 0.3 --w_hi 1.1 --w_step 0.001 --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')
@@ -118,7 +130,8 @@ 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_ratio_qq_sq = np.nan * np.zeros((N_weights, N_sigmas))  # w rows, s cols
+sigmas_ratio_sq_input = 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):
@@ -143,21 +156,25 @@ for s, sigma in enumerate(sigmas):
        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
+           sigmas_ratio_qq_sq[w][s] = s_qq / s_sq  # weighted rounded noise sigma / weighted noise sigma
+           sigmas_ratio_sq_input[w][s] = s_sq / sigma  # weighted noise sigma / input noise sigma
 # Transpose [w][s] to have index ranges [s][w]
-sigmas_ratio_T = sigmas_ratio.transpose()
+sigmas_ratio_qq_sq_T = sigmas_ratio_qq_sq.transpose()
+sigmas_ratio_sq_input_T = sigmas_ratio_sq_input.transpose()
 sigmas_qq_T = sigmas_qq.transpose()
 sigmas_sq_T = sigmas_sq.transpose()
 
 # Plot results
 figNr = 0
+s_colors = plt.cm.jet(np.linspace(0, 1, N_sigmas))
+w_colors = plt.cm.jet(np.linspace(0, 1, N_weights))
 
 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.plot(weights, sigmas_qq_T[s], color=s_colors[s], label='s = %4.2f' % sigma)
 plt.title("Sigma of weighted quantized noise")
 plt.xlabel("Weight")
 plt.ylabel("Sigma_qq")
@@ -171,14 +188,28 @@ for s, sigma in enumerate(sigmas):
     # 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.plot(weights, sigmas_qq_T[s] / weights, color=s_colors[s], label='s = %4.2f' % sigma)
+    plt.plot(weights, sigmas[s]*np.ones(N_weights), '--', color=s_colors[s])  # 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_sq / sigma as function of weight for different
+    # input sigma. The ratio deviation from 1 tells how much quantized
+    # weighted noise deviates from the input noise. This shows that rounding
+    # noise cause a change in sigma even when weight is 1.
+    plt.plot(weights, sigmas_ratio_sq_input_T[s] / weights, color=s_colors[s], label='s = %4.2f' % sigma)
+plt.title("Relative sigma difference of weighted quantized data / input data")
+plt.xlabel("Weight")
+plt.ylabel("Relative sigma difference")
+plt.legend(loc='upper right')
+plt.grid()
+
 figNr += 1
 plt.figure(figNr)
 for s, sigma in enumerate(sigmas):
@@ -186,8 +217,8 @@ for s, sigma in enumerate(sigmas):
     # 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.plot(weights, sigmas_ratio_qq_sq_T[s], color=s_colors[s], label='s = %4.2f' % sigma)
+plt.title("Relative sigma difference of weighting before / after quantisation")
 plt.xlabel("Weight")
 plt.ylabel("Relative sigma difference")
 plt.legend(loc='upper right')
@@ -198,8 +229,8 @@ 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.plot(sigmas, sigmas_ratio_qq_sq[w], color=w_colors[w], label='w = %4.2f' % weight)
+plt.title("Relative sigma difference of weighting before / after quantisation")
 plt.xlabel("Sigma")
 plt.ylabel("Relative sigma difference (s_qq / s_sq)")
 plt.legend(loc='upper right')
@@ -207,7 +238,7 @@ plt.grid()
 
 figNr += 1
 plt.figure(figNr)
-plt.imshow(sigmas_ratio, origin='lower', interpolation='none', aspect='auto', extent=[sigma_lo, sigma_hi, weight_lo, weight_hi])
+plt.imshow(sigmas_ratio_qq_sq, origin='lower', interpolation='none', aspect='auto', extent=[sigma_lo, sigma_hi, weight_lo, weight_hi])
 plt.colorbar()
 plt.title("Relative sigma difference of weighting after / before quantisation")
 plt.xlabel("Sigma")