diff --git a/libraries/base/common/python/try_round.py b/libraries/base/common/python/try_round.py
index 7de2e47c1fd87eee6e4b59f31ea7ccb8345c5ed6..0959b31a6677248904744e3273e4a6cae7ade637 100755
--- a/libraries/base/common/python/try_round.py
+++ b/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
diff --git a/libraries/base/common/python/try_round_weight.py b/libraries/base/common/python/try_round_weight.py
index e2a6986247d7e14a2a23ec22b1c187f66f17ec45..7b8a3ec60a04224192f8610714d77899e66a34e3 100644
--- a/libraries/base/common/python/try_round_weight.py
+++ b/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()