Add try_round_weight.py.

libraries/base/common/python/try_round_weight.py

+#! /usr/bin/env python3
+###############################################################################
+#
+# Copyright 2022
+# ASTRON (Netherlands Institute for Radio Astronomy) <http://www.astron.nl/>
+# P.O.Box 2, 7990 AA Dwingeloo, The Netherlands
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+# http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
+###############################################################################
+
+# Author: Eric Kooistra
+# Date: sep 2022
+# Purpose:
+#   Try applying a weight after or before rounding.
+# Description:
+#   Model to investigate usage of subband weights on:
+#   . quantized subbands --> sigma_qq
+#   . unquantized subbands --> sigma_sq
+#   Preliminary conclusion:
+#   . 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
+#   . it may be preferred to apply the subband weights to the unquantized
+#     WPFB output.
+#   Usage:
+#   > python3 try_round_weight.py -N 195312
+
+import argparse
+
+import numpy as np
+import matplotlib
+matplotlib.use('tkagg')
+import matplotlib.pyplot as plt
+
+import common as cm
+
+# Parse arguments to derive user parameters
+_parser = argparse.ArgumentParser('try_round_weight')
+_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')
+args = _parser.parse_args()
+
+N_samples = args.N
+weight_lo = args.weight_lo
+weight_hi = args.weight_hi
+weight_step = args.weight_step
+
+# Prepare signals
+noise = np.random.randn(N_samples)
+noise /= np.std(noise)
+
+# Noise level range, 1 unit = 1 LSbit
+sigma_lo = 0.1
+sigma_hi = 5
+sigmas = np.arange(sigma_lo, sigma_hi, 0.1)
+N_sigmas = len(sigmas)
+
+# Weight range, unit weight = 1
+weights = np.arange(weight_lo, weight_hi, weight_step)
+N_weights = len(weights)
+
+# 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
+for s, sigma in enumerate(sigmas):
+   noise_s = noise * sigma
+   noise_q = np.round(noise_s)
+   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
+
+# Plot results
+figNr = 0
+
+figNr += 1
+plt.figure(figNr)
+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.legend(loc='upper right')
+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.colorbar()
+plt.title("Relative sigma difference of weighting after / before quantisation")
+plt.xlabel("Sigma")
+plt.ylabel("Weight")
+plt.grid()
+
+plt.show()