Merge branch 'L2SDP-846' into 'master'

Resolve L2SDP-846

Closes L2SDP-846

See merge request desp/hdl!305

libraries/base/common/python/try_round_weight.py

@@ -28,16 +28,29 @@
 #   . 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
+#       . Hence the disturbance (jumps) on the sigma is about proportional
+#         to 1/2**r
 # 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
 
@@ -46,7 +59,7 @@ import textwrap
 
 import numpy as np
 import matplotlib
-matplotlib.use('tkagg')
+matplotlib.use('tkagg')  # to make X11 forwarding work
 import matplotlib.pyplot as plt
 
 import common as cm
@@ -72,8 +85,15 @@ _parser = argparse.ArgumentParser(
         # 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
+        # Use -r = 6 to see effect of applying weights with more (intermediate)
+        # 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
+
+        # Reproduce plots/
+        > 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 0 --noplot --save
+        > 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 4 --noplot --save
+        > 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 --noplot --save
         \n""")),
     formatter_class=argparse.RawTextHelpFormatter)
 _parser.add_argument('-S', default=0, type=int, help='Random number seed')
@@ -86,6 +106,8 @@ _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')
+_parser.add_argument('--noplot', action="store_true", help='Do not show plots')
+_parser.add_argument('--save', action="store_true", help='Do save plots')
 args = _parser.parse_args()
 
 np.random.seed(args.S)
@@ -118,7 +140,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,26 +166,33 @@ 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
+dpi = 254  # 10 dots per mm
+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")
 plt.legend(loc='upper right')
 plt.grid()
+if args.save:
+    plt.savefig('plots/try_round_weight_r%d_w_sigma_qq.jpg' % resolution, dpi=dpi)
 
 figNr += 1
 plt.figure(figNr)
@@ -171,13 +201,31 @@ 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()
+if args.save:
+    plt.savefig('plots/try_round_weight_r%d_w_sigma_qq_normalized.jpg' % resolution, dpi=dpi)
+
+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()
+if args.save:
+    plt.savefig('plots/try_round_weight_r%d_w_sigmas_ratio_sq_input.jpg' % resolution, dpi=dpi)
 
 figNr += 1
 plt.figure(figNr)
@@ -186,32 +234,39 @@ 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')
 plt.grid()
+if args.save:
+    plt.savefig('plots/try_round_weight_r%d_w_sigmas_ratio_qq_sq.jpg' % resolution, dpi=dpi)
 
 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.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')
 plt.grid()
+if args.save:
+    plt.savefig('plots/try_round_weight_r%d_s_sigmas_ratio_qq_sq.jpg' % resolution, dpi=dpi)
 
 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")
 plt.ylabel("Weight")
 plt.grid()
+if args.save:
+    plt.savefig('plots/try_round_weight_r%d_sw_sigmas_ratio_qq_sq.jpg' % resolution, dpi=dpi)
 
-plt.show()
+if not args.noplot:
+    plt.show()

libraries/base/common/python/try_wrap.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: jan 2023
+# Purpose:
+#   Try wrapping in a summator
+# Description:
+#   Usage:
+#   > python3 try_wrap.py -N 10
+#
+
+import argparse
+import textwrap
+
+import numpy as np
+import matplotlib
+matplotlib.use('tkagg')  # to make X11 forwarding work
+import matplotlib.pyplot as plt
+
+import random
+
+import common as cm
+
+# Parse arguments to derive user parameters
+_parser = argparse.ArgumentParser(
+    description="".join(textwrap.dedent("""\
+        The int_wrap function keeps w LSbits and remove MSbits from an integer,
+        so that the result is in range range -2**(w-1) to + 2**(w-1)-1. This
+        try_wrap.py show that wrap(sum) = wrap(sum(wrap)), so intermediate
+        wrapping in a summator does not change the final sum.
+
+        The wrap function is distributive, similar to modulo [1]:
+
+            (a + b) mod n = [(a mod n) + (b mod n)] mod n
+                 ab mod n = [(a mod n)(b mod n)] mod n
+
+        > python try_wrap.py -N 10 --w_inp 10 --w_sum 4
+
+        References:
+        [1] https://en.wikipedia.org/wiki/Modulo_operation
+        \n""")),
+    formatter_class=argparse.RawTextHelpFormatter)
+_parser.add_argument('-N', default=1000, type=int, help='Number of input samples')
+_parser.add_argument('--w_inp', default=10, type=int, help='Number bits of input samples')
+_parser.add_argument('--w_sum', default=4, type=int, help='Number bits of output sum')
+args = _parser.parse_args()
+
+N_samples = args.N
+W_input = args.w_inp
+W_sum = args.w_sum
+
+lo = -1 * 2**(W_input - 1)
+hi = 2**(W_input - 1) - 1
+x = [random.randint(lo, hi) for i in range(N_samples)]
+x_wrap = [cm.int_wrap(i, W_sum) for i in x]
+x_wrap_sum_wrap = cm.int_wrap(cm.add_list_elements(x_wrap), W_sum)
+x_wrap_sum      = cm.int_wrap(cm.add_list_elements(x), W_sum)
+
+print(x)
+print()
+print(x_wrap)
+print()
+
+if x_wrap_sum_wrap == x_wrap_sum:
+    print('OK')
+else:
+    print('Error:')
+    print(x_wrap_sum_wrap)
+    print(x_wrap_sum)
+

libraries/dsp/wpfb/src/vhdl/wpfb_pkg.vhd

@@ -82,6 +82,9 @@ package wpfb_pkg is
   -- LOFAR2 subband filter
   -----------------------------------------------------------------------------
 
+  -- Use guard_w = 1, instead of 2 to avoid overflow in first FFT stage,
+  -- because fil_backoff_w = 1 already provides sufficient FFT input margin.
+
   -- Fsub settings:
   -- . Settings used on LTS and DTS until at least March 2022
   constant c_wpfb_lofar2_subbands_lts_2021 : t_wpfb := (1, 1024, 0, 6,
@@ -97,7 +100,7 @@ package wpfb_pkg is
                                                        true, false, true, 23, 18, 1, 24, 1, true, 54, 2, 195313,
                                                        c_fft_pipeline, c_fft_pipeline, c_fil_ppf_pipeline);
 
-  -- . Settings used on DTS with fft_out_dat_w = 19b, to preserve FFT processing gain of 4.5 bits
+  -- . Settings used in tb_tb_verify_pfb_wg with fft_out_dat_w = 19b, to preserve FFT processing gain of 5 bits
   --   - use stage_dat_w = 25 --> fil_out_dat_w = fft_in_dat_w = 24
   --   - with fft_out_dat_w = 19 --> stat_data_w = 2*19 + 18 = 56 b
   constant c_wpfb_lofar2_subbands_dts_19b : t_wpfb := (1, 1024, 0, 6,
@@ -105,14 +108,26 @@ package wpfb_pkg is
                                                        true, false, true, 24, 19, 1, 25, 1, true, 56, 2, 195313,
                                                        c_fft_pipeline, c_fft_pipeline, c_fil_ppf_pipeline);
 
+  -- . Settings for L2TS with fft_out_gain_w = 2b, to have W_fsub_gain = W_fft_proc = 5b
+  constant c_wpfb_lofar2_subbands_l2ts_18b : t_wpfb := (1, 1024, 0, 6,
+                                                       16, 1, 14, 23, 16,
+                                                       true, false, true, 23, 18, 2, 24, 1, true, 54, 2, 195313,
+                                                       c_fft_pipeline, c_fft_pipeline, c_fil_ppf_pipeline);
+
   constant c_wpfb_lofar2_subbands : t_wpfb := c_wpfb_lofar2_subbands_dts_18b;
 
   -- The FFT output has more bits to be able to preserve the sensitivity of
   -- the processing gain of the FFT. The FFT has a processing gain of
-  -- sqrt(N_sub = N_fft / 2 = 512), so 4.5 bits. Therefore choose
-  -- fft_out_dat_w = fil_in_dat_w + 5 = 14 + 5 = 19b. Using fft_out_gain_w =
-  -- 1 compensates for the fil_backoff_w = 1 of the FIR filter. The
-  -- func_wpfb_subband_scale_w then thus returns 19 + 1 - (14 + 1) = 5 bits.
+  -- W_fft_proc = sqrt(N_fft = 1024), so 5 bits. Therefore choose
+  -- fft_out_dat_w = fil_in_dat_w + 5 = 14 + 5 = 19b. Using fft_out_gain_w
+  -- = 1 compensates for the fil_backoff_w = 1 of the FIR filter.
+  -- However, instead keep fft_out_dat_w = 18b to fit a 18x19 multiplier in
+  -- the SST. To preserve the sensitivity increase fft_out_gain_w by 1 at the
+  -- expense of loosing factor 2 (1 bit) in subband dynamic range. Therefore
+  -- fft_out_gain_w = 2 and the func_wpfb_subband_scale_w() then thus returns
+  -- (fft_out_dat_w + fft_out_gain_w) - (fil_in_dat_w + fil_backoff_w) =
+  -- (18 + 2) - (14 + 1) = 5 bits = W_fft_proc, to preserve the subband
+  -- sensitivity.
   function func_wpfb_subband_scale_w(wpfb : t_wpfb) return natural;
 
   -- The WPFB subband gain is the expected factor between subband amplitude
@@ -122,11 +137,12 @@ package wpfb_pkg is
   -- . DC gain of the FIR filter (= fir_filter_dc_gain ~= 1.0),
   -- . the FFT gain for a real input (= c_fft_real_input_gain_sine = 0.5) and
   -- . the extra bits to preserve the sensitivity of the FFT processing gain
-  --   (derived from wpfb).
+  --   W_fft_proc = 5b (derived from wpfb by func_wpfb_subband_scale_w()).
   -- For example:
   -- . func_wpfb_subband_gain() ~=  8 for c_wpfb_lofar2_subbands_lts_2021 and
   --                                  for c_wpfb_lofar2_subbands_dts_18b
   -- . func_wpfb_subband_gain() ~= 16 for c_wpfb_lofar2_subbands_dts_19b
+  --                                  for c_wpfb_lofar2_subbands_l2ts_18b
   function func_wpfb_subband_gain(wpfb : t_wpfb; fir_filter_dc_gain : real) return real;
 
   -- The expected WPFB SST level for subband amplitude A_sub and an integration