diff --git a/applications/lofar2/designs/lofar2_unb2c_sdp_station/revisions/lofar2_unb2c_sdp_station_bf/tb_lofar2_unb2c_sdp_station_bf.vhd b/applications/lofar2/designs/lofar2_unb2c_sdp_station/revisions/lofar2_unb2c_sdp_station_bf/tb_lofar2_unb2c_sdp_station_bf.vhd
index f00a6c10ef0ac8ddfcfff3f27a1ea5975db8445f..c7cca711f2f77c247415e164bbe5e6f86734710a 100644
--- a/applications/lofar2/designs/lofar2_unb2c_sdp_station/revisions/lofar2_unb2c_sdp_station_bf/tb_lofar2_unb2c_sdp_station_bf.vhd
+++ b/applications/lofar2/designs/lofar2_unb2c_sdp_station/revisions/lofar2_unb2c_sdp_station_bf/tb_lofar2_unb2c_sdp_station_bf.vhd
@@ -20,39 +20,47 @@
 
 -------------------------------------------------------------------------------
 --
--- Author: R. van der Walle, E. Kooistra
+-- Author: R. van der Walle (original), E. Kooistra (updates)
 -- Purpose: Self-checking testbench for simulating lofar2_unb2c_sdp_station_bf using WG data.
 --
 -- Description:
 --   MM control actions:
 --
---   1) Enable calc mode for WG via reg_diag_wg with:
---        freq = 19.921875MHz
---        ampl = 0.5 * 2**13
+--   1) Enable calc mode for WG on signal input g_sp via reg_diag_wg with:
+--        g_subband = 102 --> WG freq = 19.921875MHz
+--        g_ampl = 1.0 --> WG ampl = 2**13
 --   
 --   2) Read current BSN from reg_bsn_scheduler_wg and write reg_bsn_scheduler_wg 
 --      to trigger start of WG at BSN.
 --     
---   3) Read subband statistics (SST)
---   
---   4) Read beamlet statistics (BST) via ram_st_bst and verify with 
---      c_exp_beamlet_power_sp_0 at c_sdp_N_sub-1 - c_subband_sp_0.
---      View sp_beamlet_power_0  in Wave window
---   5) Compare SST with BST.
---   6) Verify 10GbE output.
---   
+--   3) Read and verify subband statistics (SST)
+--
+--   4) Select subband g_subband for beamlet g_beamlet
+--
+--   5) Apply BF weight (g_bf_gain, g_bf_phase) to g_beamlet X beam and Y beam
+--
+--   6) Read and verify beamlet statistics (BST)
+--        View sp_subband_sst in Wave window
+--        View pol_beamlet_bst in Wave window
+--
+--   7) Verify 10GbE output header and output payload for g_beamlet.
 --
 -- Usage:
 --   > as 7    # default
 --   > as 12   # for detailed debugging
+--   # Manually add missing signal
+--   > add wave -position insertpoint  \
+--     sim:/tb_lofar2_unb2c_sdp_station_bf/sp_subband_ssts_arr2 \
+--     sim:/tb_lofar2_unb2c_sdp_station_bf/pol_beamlet_bsts_arr2
 --   > run -a  
---   Takes about 40 m
+--   Takes about   40 m when g_read_all_* = FALSE
+--   Takes about 1h 5 m when g_read_all_* = TRUE
 --
 -------------------------------------------------------------------------------
 LIBRARY IEEE, common_lib, unb2c_board_lib, i2c_lib, mm_lib, dp_lib, diag_lib, lofar2_sdp_lib, wpfb_lib, tech_pll_lib, tr_10GbE_lib, lofar2_unb2c_sdp_station_lib;
 USE IEEE.std_logic_1164.ALL;
 USE IEEE.numeric_std.ALL;
-USE IEEE.MATH_REAL.ALL;
+USE IEEE.math_real.ALL;
 USE common_lib.common_pkg.ALL;
 USE common_lib.common_mem_pkg.ALL;
 USE common_lib.common_field_pkg.ALL;
@@ -70,6 +78,17 @@ USE lofar2_sdp_lib.tb_sdp_pkg.ALL;
 USE tech_pll_lib.tech_pll_component_pkg.ALL;
 
 ENTITY tb_lofar2_unb2c_sdp_station_bf IS
+  GENERIC (
+    g_sp                : NATURAL := 0;     -- WG signal path index in range(S_pn = 12)
+    g_wg_ampl           : REAL := 1.0;      -- WG normalized amplitude
+    g_subband           : NATURAL := 102;   -- select g_subband at index 102 = 102/1024 * 200MHz = 19.921875 MHz
+    g_beamlet           : NATURAL := 10;    -- map g_subband to g_beamlet index in beamset in range(c_sdp_S_sub_bf = 488)
+    g_beamlet_scale     : REAL := 1.0 / 2.0**9;  -- g_beamlet output scale factor
+    g_bf_gain           : REAL := 1.0;      -- g_beamlet BF weight normalized gain
+    g_bf_phase          : REAL := 30.0;      -- g_beamlet BF weight phase rotation in degrees
+    g_read_all_SST      : BOOLEAN := FALSE;  -- when FALSE only read SST for g_subband, to save sim time
+    g_read_all_BST      : BOOLEAN := FALSE   -- when FALSE only read BST for g_beamlet, to save sim time
+  );
 END tb_lofar2_unb2c_sdp_station_bf;
 
 ARCHITECTURE tb OF tb_lofar2_unb2c_sdp_station_bf IS
@@ -95,16 +114,19 @@ ARCHITECTURE tb OF tb_lofar2_unb2c_sdp_station_bf IS
   CONSTANT c_sa_clk_period       : TIME := tech_pll_clk_644_period; -- 644MHz
 
   CONSTANT c_tb_clk_period       : TIME := 100 ps; -- use fast tb_clk to speed up M&C
+  CONSTANT c_cross_clock_domain_latency : NATURAL := 20;
 
   CONSTANT c_nof_block_per_sync  : NATURAL := 16;
   CONSTANT c_nof_clk_per_sync    : NATURAL := c_nof_block_per_sync*c_sdp_N_fft; 
   CONSTANT c_pps_period          : NATURAL := c_nof_clk_per_sync;
   CONSTANT c_wpfb_sim            : t_wpfb := func_wpfb_set_nof_block_per_sync(c_sdp_wpfb_subbands, c_nof_block_per_sync);
-  CONSTANT c_stat_data_sz        : NATURAL := c_longword_sz/c_word_sz;  -- = 2
+  CONSTANT c_stat_data_sz        : NATURAL := c_wpfb_sim.stat_data_sz;  -- = 2
 
-  CONSTANT c_percentage          : REAL := 0.05;  -- percentage that actual value may differ from expected value
-  CONSTANT c_lo_factor           : REAL := 1.0 - c_percentage;  -- lower boundary  
-  CONSTANT c_hi_factor           : REAL := 1.0 + c_percentage;  -- higher boundary
+  CONSTANT c_stat_percentage     : REAL := 0.05;  -- +-percentage margin that actual value may differ from expected value
+  CONSTANT c_stat_lo_factor      : REAL := 1.0 - c_stat_percentage;  -- lower boundary
+  CONSTANT c_stat_hi_factor      : REAL := 1.0 + c_stat_percentage;  -- higher boundary
+
+  CONSTANT c_beamlet_output_delta : INTEGER := 2;  -- +-delta margin
 
   -- header fields
   CONSTANT c_cep_eth_dst_mac     : STD_LOGIC_VECTOR(47 DOWNTO 0) := c_sdp_cep_eth_dst_mac;   -- 00074306C700 = DOP36-eth0
@@ -117,7 +139,7 @@ ARCHITECTURE tb OF tb_lofar2_unb2c_sdp_station_bf IS
 
   CONSTANT c_exp_ip_header_checksum : NATURAL := 16#5BDE#;  -- value obtained from rx_sdp_cep_header.ip.header_checksum in wave window
 
-  CONSTANT c_exp_beamlet_scale   : NATURAL := 2**15;
+  CONSTANT c_exp_beamlet_scale   : NATURAL := NATURAL(g_beamlet_scale * REAL(c_sdp_unit_beamlet_scale));  -- c_sdp_unit_beamlet_scale = 2**15;
 
   CONSTANT c_exp_sdp_info        : t_sdp_info := (
                                      TO_UVEC(601, 16),   -- station_id
@@ -135,39 +157,72 @@ ARCHITECTURE tb OF tb_lofar2_unb2c_sdp_station_bf IS
                                    );
 
   -- WG
-  CONSTANT c_full_scale_ampl      : REAL := REAL(2**(14-1)-1);  -- = full scale of WG
   CONSTANT c_bsn_start_wg         : NATURAL := c_init_bsn + 2;  -- start WG at this BSN to instead of some BSN, to avoid mismatches in exact expected data values
-  CONSTANT c_ampl_sp_0            : NATURAL := 2**(c_sdp_W_adc-1) / 2;  -- in number of lsb
+  -- .ampl
+  CONSTANT c_wg_ampl_full_scale   : NATURAL := 2**(c_sdp_W_adc-1);  -- full scale (FS) of WG, will just cause clipping of +FS to +FS-1
+  CONSTANT c_wg_ampl_lsb          : REAL := c_diag_wg_ampl_unit / REAL(c_wg_ampl_full_scale);  -- amplitude in number of LSbit resolution steps
+  CONSTANT c_wg_ampl              : NATURAL := NATURAL(g_wg_ampl * REAL(c_wg_ampl_full_scale));  -- in number of lsb
+  CONSTANT c_exp_sp_power         : REAL := REAL(c_wg_ampl**2) / 2.0;
+  CONSTANT c_exp_sp_ast           : REAL := c_exp_sp_power * REAL(c_nof_clk_per_sync);
+  -- . phase
+  CONSTANT c_subband_phase        : REAL := 0.0;  -- wanted subband phase in degrees = WG phase at sop
+  CONSTANT c_subband_freq         : REAL := REAL(g_subband) / REAL(c_sdp_N_fft);  -- normalized by fs = f_adc = 200 MHz = dp_clk rate
+  CONSTANT c_wg_latency           : INTEGER := c_diag_wg_latency - 0;  -- -0 to account for BSN scheduler start trigger latency
+  CONSTANT c_wg_phase_offset      : REAL := 360.0 * REAL(c_wg_latency) * c_subband_freq;  -- c_diag_wg_latency is in dp_clk cycles
+  CONSTANT c_wg_phase             : REAL := c_subband_phase + c_wg_phase_offset;  -- WG phase in degrees
+  -- . freq
   CONSTANT c_wg_subband_freq_unit : REAL := c_diag_wg_freq_unit/REAL(c_sdp_N_fft);  -- subband freq = Fs/1024 = 200 MSps/1024 = 195312.5 Hz sinus
-  CONSTANT c_wg_freq_offset       : REAL := 0.0/11.0; -- in freq_unit
-  CONSTANT c_subband_sp_0         : REAL := 102.0;  -- Select subband at index 102 = 102/1024 * 200MHz = 19.921875 MHz 
-  CONSTANT c_wg_ampl_lsb          : REAL := c_diag_wg_ampl_unit / c_full_scale_ampl;  -- amplitude in number of LSbit resolution steps
-  CONSTANT c_exp_wg_power_sp_0    : REAL := REAL(c_ampl_sp_0**2)/2.0 * REAL(c_nof_clk_per_sync);
 
   -- WPFB
-  CONSTANT c_wb_leakage_bin                 : NATURAL := c_wpfb_sim.nof_points / c_wpfb_sim.wb_factor;   -- = 256, leakage will occur in this bin if FIR wb_factor is reversed 
-  CONSTANT c_exp_sp_beamlet_power_ratio     : REAL := 1.0/8.0;   -- depends on internal WPFB quantization and FIR coefficients
-  CONSTANT c_exp_sp_beamlet_power_sum_ratio : REAL := c_exp_sp_beamlet_power_ratio;   -- because all sinus power is expected in one subband
-  CONSTANT c_exp_beamlet_power_sp_0         : REAL := c_exp_wg_power_sp_0 * c_exp_sp_beamlet_power_ratio;
+  CONSTANT c_pol_index                      : NATURAL := g_sp MOD c_sdp_Q_fft;
+  CONSTANT c_pfb_index                      : NATURAL := g_sp / c_sdp_Q_fft;  -- only read used WPFB unit out of range(c_sdp_P_pfb = 6)
+  CONSTANT c_subband_phase_offset           : REAL := -90.0;  -- WG with zero phase sinues yields subband with -90 degrees phase (negative Im, zero Re)
+  CONSTANT c_subband_weight_gain            : REAL := 1.0;  -- use default unit subband weights
+  CONSTANT c_subband_weight_phase           : REAL := 0.0;  -- use default unit subband weights
+  CONSTANT c_exp_subband_sp_ampl_ratio      : REAL := 7.96;  -- ~= 8 for unit FIR DC gain, depends on internal WPFB quantization and FIR coefficients
+  CONSTANT c_exp_subband_ampl               : REAL := REAL(c_wg_ampl) * c_exp_subband_sp_ampl_ratio * c_subband_weight_gain;
+  CONSTANT c_exp_subband_power              : REAL := c_exp_subband_ampl**2.0;  -- complex, so no divide by 2
+  CONSTANT c_exp_subband_sst                : REAL := c_exp_subband_power * REAL(c_nof_block_per_sync);
 
   TYPE t_real_arr IS ARRAY (INTEGER RANGE <>) OF REAL; 
-  TYPE t_slv_64_subbands_arr IS ARRAY (INTEGER RANGE <>) OF t_slv_64_arr(0 TO c_sdp_N_sub);
-  TYPE t_slv_64_beamlets_arr IS ARRAY (INTEGER RANGE <>) OF t_slv_64_arr(0 TO c_sdp_S_sub_bf);
+  TYPE t_slv_64_subbands_arr IS ARRAY (INTEGER RANGE <>) OF t_slv_64_arr(0 TO c_sdp_N_sub-1);           -- 512
+  TYPE t_slv_64_beamlets_arr IS ARRAY (INTEGER RANGE <>) OF t_slv_64_arr(0 TO c_sdp_N_beamlets_sdp-1);  -- 2*488 = 976
 
   -- BF
-  CONSTANT c_exp_beamlet_index : NATURAL := NATURAL(c_subband_sp_0) * c_sdp_N_pol_bf;
-  CONSTANT c_exp_beamlet_re    : STD_LOGIC_VECTOR(7 DOWNTO 0) := x"81"; -- = -127, derived from simulation
-  CONSTANT c_exp_beamlet_im    : STD_LOGIC_VECTOR(7 DOWNTO 0) := x"7F"; -- = +127, derived from simulation
+  -- . select
+  CONSTANT c_exp_beamlet_index        : NATURAL := g_beamlet * c_sdp_N_pol_bf;  -- in beamset 0
+  -- . Beamlet weights for selected g_sp
+  CONSTANT c_bf_weight_re             : INTEGER := INTEGER(g_bf_gain * REAL(c_sdp_unit_bf_weight) * COS(g_bf_phase * MATH_2_PI / 360.0));
+  CONSTANT c_bf_weight_im             : INTEGER := INTEGER(g_bf_gain * REAL(c_sdp_unit_bf_weight) * SIN(g_bf_phase * MATH_2_PI / 360.0));
+  -- . Beamlet internal
+  CONSTANT c_exp_beamlet_ampl         : REAL := c_exp_subband_ampl * g_bf_gain;
+  CONSTANT c_exp_beamlet_phase        : REAL := c_subband_phase_offset + c_subband_weight_phase + g_bf_phase;
+  CONSTANT c_exp_beamlet_re           : REAL := c_exp_beamlet_ampl * COS(c_exp_beamlet_phase * MATH_2_PI / 360.0);
+  CONSTANT c_exp_beamlet_im           : REAL := c_exp_beamlet_ampl * SIN(c_exp_beamlet_phase * MATH_2_PI / 360.0);
+  -- . BST
+  CONSTANT c_exp_beamlet_power        : REAL := c_exp_beamlet_ampl**2.0;  -- complex, so no divide by 2
+  CONSTANT c_exp_beamlet_bst          : REAL := c_exp_subband_sst * g_bf_gain**2.0;  -- = c_exp_beamlet_power *  REAL(c_nof_block_per_sync)
+  -- . Beamlet output
+  CONSTANT c_exp_beamlet_output_ampl  : REAL := c_exp_beamlet_ampl * g_beamlet_scale;
+  CONSTANT c_exp_beamlet_output_phase : REAL := c_exp_beamlet_phase;
+  CONSTANT c_exp_beamlet_output_re    : REAL := c_exp_beamlet_re * g_beamlet_scale;
+  CONSTANT c_exp_beamlet_output_im    : REAL := c_exp_beamlet_im * g_beamlet_scale;
 
   -- MM
   -- . Address widths of a single MM instance
+  --   . c_sdp_S_pn = 12 instances
+  CONSTANT c_addr_w_reg_diag_wg     : NATURAL := 2;
+  --   . c_sdp_N_beamsets = 2 instances
   CONSTANT c_addr_w_ram_ss_ss_wide : NATURAL := ceil_log2(c_sdp_P_pfb * c_sdp_S_sub_bf * c_sdp_Q_fft);
-  CONSTANT c_addr_w_ram_bf_weights : NATURAL := ceil_log2(c_sdp_N_pol * c_sdp_P_pfb * c_sdp_S_sub_bf * c_sdp_Q_fft);
+  CONSTANT c_addr_w_ram_bf_weights  : NATURAL := ceil_log2(c_sdp_N_pol_bf * c_sdp_P_pfb * c_sdp_S_sub_bf * c_sdp_Q_fft);
   CONSTANT c_addr_w_reg_bf_scale   : NATURAL := 1;
   CONSTANT c_addr_w_reg_hdr_dat    : NATURAL := ceil_log2(field_nof_words(c_sdp_cep_hdr_field_arr, c_word_w));
   CONSTANT c_addr_w_reg_dp_xonoff  : NATURAL := 1;
-  CONSTANT c_addr_w_ram_st_bst     : NATURAL := ceil_log2(c_sdp_S_sub_bf*c_sdp_N_pol*(c_longword_sz/c_word_sz));
+  CONSTANT c_addr_w_ram_st_bst      : NATURAL := ceil_log2(c_sdp_S_sub_bf*c_sdp_N_pol_bf*c_stat_data_sz);
   -- . Address spans of a single MM instance
+  --   . c_sdp_S_pn = 12 instances
+  CONSTANT c_mm_span_reg_diag_wg    : NATURAL := 2**c_addr_w_reg_diag_wg;
+  --   . c_sdp_N_beamsets = 2 instances
   CONSTANT c_mm_span_ram_ss_ss_wide : NATURAL := 2**c_addr_w_ram_ss_ss_wide;
   CONSTANT c_mm_span_ram_bf_weights : NATURAL := 2**c_addr_w_ram_bf_weights;
   CONSTANT c_mm_span_reg_bf_scale   : NATURAL := 2**c_addr_w_reg_bf_scale;
@@ -179,11 +234,15 @@ ARCHITECTURE tb OF tb_lofar2_unb2c_sdp_station_bf IS
   CONSTANT c_mm_file_reg_bsn_source_v2    : STRING := mmf_unb_file_prefix(c_unb_nr, c_node_nr) & "REG_BSN_SOURCE_V2";
   CONSTANT c_mm_file_reg_bsn_scheduler_wg : STRING := mmf_unb_file_prefix(c_unb_nr, c_node_nr) & "REG_BSN_SCHEDULER";
   CONSTANT c_mm_file_reg_diag_wg          : STRING := mmf_unb_file_prefix(c_unb_nr, c_node_nr) & "REG_WG";
+  CONSTANT c_mm_file_ram_equalizer_gains  : STRING := mmf_unb_file_prefix(c_unb_nr, c_node_nr) & "RAM_EQUALIZER_GAINS";
+  CONSTANT c_mm_file_reg_dp_selector      : STRING := mmf_unb_file_prefix(c_unb_nr, c_node_nr) & "REG_DP_SELECTOR";
+  CONSTANT c_mm_file_ram_st_sst           : STRING := mmf_unb_file_prefix(c_unb_nr, c_node_nr) & "RAM_ST_SST";
   CONSTANT c_mm_file_ram_st_bst           : STRING := mmf_unb_file_prefix(c_unb_nr, c_node_nr) & "RAM_ST_BST";
   CONSTANT c_mm_file_reg_dp_xonoff        : STRING := mmf_unb_file_prefix(c_unb_nr, c_node_nr) & "REG_DP_XONOFF";
-  CONSTANT c_mm_file_ram_st_sst           : STRING := mmf_unb_file_prefix(c_unb_nr, c_node_nr) & "RAM_ST_SST";
-  CONSTANT c_mm_file_reg_sdp_info         : STRING := mmf_unb_file_prefix(c_unb_nr, c_node_nr) & "REG_SDP_INFO";
+  CONSTANT c_mm_file_ram_ss_ss_wide       : STRING := mmf_unb_file_prefix(c_unb_nr, c_node_nr) & "RAM_SS_SS_WIDE";
+  CONSTANT c_mm_file_ram_bf_weights       : STRING := mmf_unb_file_prefix(c_unb_nr, c_node_nr) & "RAM_BF_WEIGHTS";
   CONSTANT c_mm_file_reg_bf_scale         : STRING := mmf_unb_file_prefix(c_unb_nr, c_node_nr) & "REG_BF_SCALE";
+  CONSTANT c_mm_file_reg_sdp_info         : STRING := mmf_unb_file_prefix(c_unb_nr, c_node_nr) & "REG_SDP_INFO";
   CONSTANT c_mm_file_reg_hdr_dat          : STRING := mmf_unb_file_prefix(c_unb_nr, c_node_nr) & "REG_HDR_DAT";  -- c_sdp_N_beamsets = 2 beamsets
 
   -- Tb
@@ -211,16 +270,32 @@ ARCHITECTURE tb OF tb_lofar2_unb2c_sdp_station_bf IS
   -- WG
   SIGNAL current_bsn_wg                 : STD_LOGIC_VECTOR(c_dp_stream_bsn_w-1 DOWNTO 0);
 
-  -- WPFB
-  SIGNAL sp_subband_powers_arr2         : t_slv_64_subbands_arr(c_sdp_N_pol-1 DOWNTO 0);   -- [sp][sub]
+  -- FSUB
+  -- . Read sp_subband_ssts_arr2 = SST for one WPFB unit that processes g_sp
+  SIGNAL sp_subband_ssts_arr2  : t_slv_64_subbands_arr(c_sdp_N_pol-1 DOWNTO 0);   -- [pol][sub], for X,Y pair of A, B
+  SIGNAL sp_subband_sst        : REAL := 0.0;
+  SIGNAL stat_data               : STD_LOGIC_VECTOR(c_longword_w-1 DOWNTO 0);
+
+  -- . Selector
+  SIGNAL sst_offload_weighted_subbands : STD_LOGIC;
+
+  -- . Subband equalizer
+  SIGNAL sp_subband_weight_re    : INTEGER := 0;
+  SIGNAL sp_subband_weight_im    : INTEGER := 0;
+  SIGNAL sp_subband_weight_gain  : REAL := 0.0;
+  SIGNAL sp_subband_weight_phase : REAL := 0.0;
+
+  -- BF
+  SIGNAL sp_subband_select       : NATURAL :=  0;
+  SIGNAL sp_subband_select_arr   : t_natural_arr(0 TO c_sdp_S_sub_bf * c_sdp_N_pol-1) := (OTHERS => 0);  -- Q_fft = N_pol = 2
+  SIGNAL sp_bf_weights_re_arr    : t_integer_arr(0 TO c_sdp_S_pn-1) := (OTHERS => 0);
+  SIGNAL sp_bf_weights_im_arr    : t_integer_arr(0 TO c_sdp_S_pn-1) := (OTHERS => 0);
+  SIGNAL sp_bf_weights_gain_arr  : t_real_arr(0 TO c_sdp_S_pn-1) := (OTHERS => 0.0);
+  SIGNAL sp_bf_weights_phase_arr : t_real_arr(0 TO c_sdp_S_pn-1) := (OTHERS => 0.0);
 
-  SIGNAL sp_beamlet_powers_arr2         : t_slv_64_beamlets_arr(c_sdp_N_beamsets*c_sdp_N_pol_bf-1 DOWNTO 0);   -- [sp][sub]
-  SIGNAL sp_beamlet_power_0             : REAL;
-  SIGNAL sp_beamlet_power_sum           : t_real_arr(c_sdp_N_beamsets*c_sdp_N_pol_bf-1 DOWNTO 0) := (OTHERS=>0.0);
-  SIGNAL sp_beamlet_power_sum_0         : REAL;
-  SIGNAL sp_beamlet_power_ratio_0       : REAL;
-  SIGNAL sp_beamlet_power_sum_ratio_0   : REAL;
-  SIGNAL sp_beamlet_power_leakage_sum_0 : REAL;
+  SIGNAL pol_beamlet_bsts_arr2 : t_slv_64_beamlets_arr(c_sdp_N_pol_bf-1 DOWNTO 0);  -- [pol_bf][blet]
+  SIGNAL pol_beamlet_bst_X_arr : t_real_arr(0 TO c_sdp_N_beamsets-1) := (OTHERS => 0.0);  -- [bset]
+  SIGNAL pol_beamlet_bst_Y_arr : t_real_arr(0 TO c_sdp_N_beamsets-1) := (OTHERS => 0.0);  -- [bset]
 
   -- 10GbE
   SIGNAL rx_beamlet_arr_re   : t_slv_8_arr(c_sdp_cep_nof_blocks_per_packet-1 DOWNTO 0);   -- [3:0]
@@ -309,7 +384,7 @@ BEGIN
     g_sim_node_nr            => c_node_nr,
     g_wpfb                   => c_wpfb_sim,
     g_bsn_nof_clk_per_sync   => c_nof_clk_per_sync,
-    g_scope_selected_subband => NATURAL(c_subband_sp_0)
+    g_scope_selected_subband => g_subband
   )
   PORT MAP (
     -- GENERAL
@@ -421,12 +496,14 @@ BEGIN
   
   p_mm_stimuli : PROCESS
     VARIABLE v_bsn                   : NATURAL;
-    VARIABLE v_sp_beamlet_power      : REAL;
-    VARIABLE v_sp_subband_power      : REAL;
-    VARIABLE v_W, v_T, v_U, v_S, v_B : NATURAL;  -- array indicies
-    VARIABLE v_re, v_im              : INTEGER;
-    VARIABLE v_re_exp, v_im_exp      : INTEGER;
-    VARIABLE v_offset                : NATURAL;
+    VARIABLE v_sp_subband_sst                       : REAL := 0.0;
+    VARIABLE v_pol_beamlet_bst                      : REAL := 0.0;
+    VARIABLE v_data_lo, v_data_hi                   : STD_LOGIC_VECTOR(c_word_w-1 DOWNTO 0);
+    VARIABLE v_stat_data                            : STD_LOGIC_VECTOR(c_longword_w-1 DOWNTO 0);
+    VARIABLE v_len, v_span, v_offset, v_addr, v_sel : NATURAL;  -- address ranges, indices
+    VARIABLE v_W, v_P, v_S, v_A, v_B, v_G           : NATURAL;  -- array indicies
+    VARIABLE v_re, v_im, v_weight                   : INTEGER;
+    VARIABLE v_re_exp, v_im_exp                     : REAL := 0.0;
   BEGIN
     -- Wait for DUT power up after reset
     WAIT FOR 1 us;
@@ -474,8 +551,15 @@ BEGIN
     ------------------------------------------------------------------------------
     FOR bset IN 0 TO c_sdp_N_beamsets-1 LOOP
       -- MM beamlet_scale
+      -- . write
       v_offset := bset * c_mm_span_reg_bf_scale;
-      mmf_mm_bus_rd(c_mm_file_reg_bf_scale,  v_offset + 0, rd_data, tb_clk); rd_beamlet_scale <= rd_data(15 DOWNTO 0);
+      mmf_mm_bus_wr(c_mm_file_reg_bf_scale, v_offset + 0, c_exp_beamlet_scale, tb_clk);
+      proc_common_wait_some_cycles(tb_clk, c_cross_clock_domain_latency);
+      proc_common_wait_some_cycles(ext_clk, c_cross_clock_domain_latency);
+
+      -- . readback
+      mmf_mm_bus_rd(c_mm_file_reg_bf_scale, v_offset + 0, rd_data, tb_clk);
+      rd_beamlet_scale <= rd_data(15 DOWNTO 0);
       proc_common_wait_some_cycles(tb_clk, 1);
       ASSERT TO_UINT(rd_beamlet_scale) = c_exp_beamlet_scale REPORT "Wrong MM read beamlet_scale for beamset " & NATURAL'IMAGE(bset) SEVERITY ERROR;
 
@@ -560,17 +644,19 @@ BEGIN
     pps_rst <= '0';
 
     ----------------------------------------------------------------------------
-    -- Enable WG
+    -- Enable and start WG
     ----------------------------------------------------------------------------
     --   0 : mode[7:0]           --> off=0, calc=1, repeat=2, single=3)
     --       nof_samples[31:16]  --> <= c_ram_wg_size=1024
     --   1 : phase[15:0]
     --   2 : freq[30:0]
     --   3 : ampl[16:0]
-    mmf_mm_bus_wr(c_mm_file_reg_diag_wg, 0, 1024*2**16 + 1, tb_clk);  -- nof_samples, mode calc
-    mmf_mm_bus_wr(c_mm_file_reg_diag_wg, 1, INTEGER(  0.0 * c_diag_wg_phase_unit), tb_clk);  -- phase offset in degrees
-    mmf_mm_bus_wr(c_mm_file_reg_diag_wg, 2, INTEGER((c_subband_sp_0+c_wg_freq_offset) * c_wg_subband_freq_unit), tb_clk);  -- freq
-    mmf_mm_bus_wr(c_mm_file_reg_diag_wg, 3, INTEGER(REAL(c_ampl_sp_0) * c_wg_ampl_lsb), tb_clk);  -- ampl
+    -- . Put wanted signal on g_sp input
+    v_offset := g_sp * c_mm_span_reg_diag_wg;
+    mmf_mm_bus_wr(c_mm_file_reg_diag_wg, v_offset + 0, 1024*2**16 + 1, tb_clk);  -- nof_samples, mode calc
+    mmf_mm_bus_wr(c_mm_file_reg_diag_wg, v_offset + 1, INTEGER(c_wg_phase * c_diag_wg_phase_unit), tb_clk);  -- phase offset in degrees
+    mmf_mm_bus_wr(c_mm_file_reg_diag_wg, v_offset + 2, INTEGER(REAL(g_subband) * c_wg_subband_freq_unit), tb_clk);  -- freq
+    mmf_mm_bus_wr(c_mm_file_reg_diag_wg, v_offset + 3, INTEGER(REAL(c_wg_ampl) * c_wg_ampl_lsb), tb_clk);  -- ampl
 
     -- Read current BSN
     mmf_mm_bus_rd(c_mm_file_reg_bsn_scheduler_wg, 0, current_bsn_wg(31 DOWNTO  0), tb_clk);
@@ -584,7 +670,117 @@ BEGIN
     mmf_mm_bus_wr(c_mm_file_reg_bsn_scheduler_wg, 0, v_bsn, tb_clk);  -- first write low then high part
     mmf_mm_bus_wr(c_mm_file_reg_bsn_scheduler_wg, 1,     0, tb_clk);  -- assume v_bsn < 2**31-1
     
+    ----------------------------------------------------------------------------
+    -- Read weighted subband selector
+    ----------------------------------------------------------------------------
+    mmf_mm_bus_rd(c_mm_file_reg_dp_selector, 0, rd_data, tb_clk);
+    sst_offload_weighted_subbands <= NOT rd_data(0);
+    proc_common_wait_some_cycles(tb_clk, 1);
+
+    ----------------------------------------------------------------------------
+    -- Subband weight
+    ----------------------------------------------------------------------------
+
+    -- . MM format: (cint16)RAM_EQUALIZER_GAINS[S_pn/Q_fft]_[Q_fft][N_sub] = [S_pn][N_sub]
+    v_addr := g_sp * c_sdp_N_sub + g_subband;
+    -- . read
+    mmf_mm_bus_rd(c_mm_file_ram_equalizer_gains, v_addr, rd_data, tb_clk);
+    v_re := unpack_complex_re(rd_data, c_sdp_W_sub_weight);
+    v_im := unpack_complex_im(rd_data, c_sdp_W_sub_weight);
+    sp_subband_weight_re <= v_re;
+    sp_subband_weight_im <= v_im;
+    sp_subband_weight_gain <= SQRT(REAL(v_re)**2.0 + REAL(v_im)**2.0) / REAL(c_sdp_unit_sub_weight);
+    sp_subband_weight_phase <= atan2(Y => REAL(v_im), X => REAL(v_re)) * 360.0 / MATH_2_PI;
+
+    -- No need to write subband weight, because default it is unit weight
+
+    ----------------------------------------------------------------------------
+    -- Subband select to map subband to beamlet
+    ----------------------------------------------------------------------------
+    -- . MM format: (uint16)RAM_SS_SS_WIDE[N_beamsets][A_pn]_[S_sub_bf][Q_fft], Q_fft = N_pol = 2
+
+    -- . write selection, only for g_beamlet to save sim time
+    v_span := true_log_pow2(c_sdp_N_pol * c_sdp_S_sub_bf);  -- = 1024
+    FOR U IN 0 TO c_sdp_N_beamsets-1 LOOP
+      -- Same selection for both beamsets
+      -- Select beamlet g_beamlet to subband g_subband
+      FOR A IN 0 TO c_sdp_A_pn-1 LOOP
+        -- Same selection to all SP
+        FOR P IN 0 TO c_sdp_N_pol-1 LOOP
+          v_addr := P + g_beamlet * c_sdp_N_pol + A * v_span + U * c_mm_span_ram_ss_ss_wide;
+          v_sel := P + g_subband * c_sdp_N_pol;
+          mmf_mm_bus_wr(c_mm_file_ram_ss_ss_wide, v_addr, v_sel, tb_clk);
+        END LOOP;
+      END LOOP;
+    END LOOP;
+
+    -- . read back selection for g_sp = c_pfb_index * c_sdp_N_pol + c_pol_index
+    v_P := c_pol_index;
+    v_A := c_pfb_index;
+    FOR U IN 0 TO c_sdp_N_beamsets-1 LOOP
+      -- Same selection for both beamsets, so fine to use only one sp_subband_select_arr()
+      FOR B IN 0 TO c_sdp_S_sub_bf-1 LOOP
+        -- Same selection for all SP, so fine to only read subband selection for g_sp
+        v_addr := v_P + B * c_sdp_N_pol + v_A * v_span + U * c_mm_span_ram_ss_ss_wide;
+        mmf_mm_bus_rd(c_mm_file_ram_ss_ss_wide, v_addr, rd_data, tb_clk);
+        v_sel := (TO_UINT(rd_data) - v_P) / c_sdp_N_pol;
+        sp_subband_select_arr(B) <= v_sel;
+        sp_subband_select <= v_sel;  -- for time series view in Wave window
+      END LOOP;
+    END LOOP;
+    proc_common_wait_some_cycles(tb_clk, 1);
+    proc_common_wait_some_cycles(ext_clk, 100);  -- delay for ease of view in Wave window
+
+    ----------------------------------------------------------------------------
+    -- Write beamlet weight for g_beamlet
+    ----------------------------------------------------------------------------
+    -- . MM format: (cint16)RAM_BF_WEIGHTS[N_beamsets][N_pol_bf][A_pn]_[N_pol][S_sub_bf]
+
+    -- . write BF weights, only for g_beamlet to save sim time
+    v_span := true_log_pow2(c_sdp_N_pol * c_sdp_S_sub_bf);  -- = 1024
+    FOR U IN 0 TO c_sdp_N_beamsets-1 LOOP
+      -- Same BF weights for both beamsets
+      FOR A IN 0 TO c_sdp_A_pn-1 LOOP
+        FOR P IN 0 TO c_sdp_N_pol-1 LOOP
+          v_S := A * c_sdp_N_pol + P;
+          IF v_S = g_sp THEN
+            -- use generic BF weight for g_sp in g_beamlet
+            v_weight := pack_complex(re => c_bf_weight_re, im => c_bf_weight_im, w => c_sdp_W_bf_weight);
+          ELSE
+            -- default set all weights to zero
+            v_weight := 0;
+          END IF;
+          v_addr := g_beamlet + A * v_span + U * c_mm_span_ram_bf_weights;
+          v_addr := v_addr + P * c_sdp_S_sub_bf;
+          mmf_mm_bus_wr(c_mm_file_ram_bf_weights, v_addr, v_weight, tb_clk);
+        END LOOP;
+      END LOOP;
+    END LOOP;
+
+    -- . read back BF weights for g_beamlet
+    FOR U IN 0 TO c_sdp_N_beamsets-1 LOOP
+      FOR A IN 0 TO c_sdp_A_pn-1 LOOP
+        FOR P IN 0 TO c_sdp_N_pol-1 LOOP
+          v_addr := g_beamlet + A * v_span + U * c_mm_span_ram_bf_weights;
+          v_addr := v_addr + P * c_sdp_S_sub_bf;
+          mmf_mm_bus_rd(c_mm_file_ram_bf_weights, v_addr, rd_data, tb_clk);
+          v_re := unpack_complex_re(rd_data, c_sdp_W_bf_weight);
+          v_im := unpack_complex_im(rd_data, c_sdp_W_bf_weight);
+          -- same BF weights for both beamsets, so fine to use only one sp_bf_weights_*_arr()
+          v_S := A * c_sdp_N_pol + P;
+          sp_bf_weights_re_arr(v_S) <= v_re;
+          sp_bf_weights_im_arr(v_S) <= v_im;
+          sp_bf_weights_gain_arr(v_S) <= SQRT(REAL(v_re)**2.0 + REAL(v_im)**2.0) / REAL(c_sdp_unit_bf_weight);
+          sp_bf_weights_phase_arr(v_S) <= atan2(Y => REAL(v_im), X => REAL(v_re)) * 360.0 / MATH_2_PI;
+        END LOOP;
+      END LOOP;
+    END LOOP;
+    proc_common_wait_some_cycles(tb_clk, 1);
+    proc_common_wait_some_cycles(ext_clk, 100);  -- delay for ease of view in Wave window
+
+    ----------------------------------------------------------------------------
     -- Wait for enough WG data and start of sync interval
+    ----------------------------------------------------------------------------
     mmf_mm_wait_until_value(c_mm_file_reg_bsn_scheduler_wg, 0,                               -- read BSN low
                             "UNSIGNED", rd_data, ">=", c_init_bsn + c_nof_block_per_sync*3,  -- this is the wait until condition
                             c_sdp_T_sub, tb_clk);
@@ -595,135 +791,186 @@ BEGIN
     ---------------------------------------------------------------------------
     -- Read subband statistics
     ---------------------------------------------------------------------------   
-    -- . the subband statistics are c_wpfb_sim.stat_data_sz = 2 word power values.
-    -- . there are c_sdp_N_sub = 512 subbands per signal path
-    -- . one complex WPFB can process two real inputs A, B, is c_sdp_Q_fft = c_sdp_N_pol = 2
+    -- . the subband statistics are c_stat_data_sz = 2 word power values.
+    -- . there are c_sdp_S_pn = 12 signal inputs A, B, C, D, E, F, G, H, I, J, K, L
+    -- . there are c_sdp_N_sub = 512 subbands per signal input (SI, = signal path, SP)
+    -- . one complex WPFB can process two real inputs A, B, so there are c_sdp_P_pfb = 6 WPFB units,
+    --   but only read for the 1 WPFB unit of the selected g_sp, to save sim time
+    -- . the outputs for A, B are time multiplexed, c_sdp_Q_fft = 2, assume that they
+    --   correspond to the c_sdp_N_pol = 2 signal polarizations
     -- . the subbands are output alternately so A0 B0 A1 B1 ... A511 B511 for input A, B
     -- . the subband statistics multiple WPFB units appear in order in the ram_st_sst address map
     -- . the subband statistics are stored first lo word 0 then hi word 1
-    
-    FOR I IN 0 TO c_sdp_N_pol*c_sdp_N_sub*c_stat_data_sz-1 LOOP  -- 2048 = 2 * 512 * 64/32
-      v_W := I MOD c_stat_data_sz;                               -- 0, 1 per statistics word
-      v_T := (I / c_stat_data_sz) MOD c_sdp_N_pol;               -- 0, 1 per pol
-      v_U := I / (c_sdp_N_pol*c_stat_data_sz*c_sdp_N_sub);       -- / 2048
-      v_S := v_T + v_U * c_sdp_N_pol;
-      v_B := (I / (c_sdp_N_pol*c_stat_data_sz)) MOD c_sdp_N_sub; -- 0:511 per dual pol
-      -- Only read sp 0, pol 0 (v_S = 0)
-      IF v_S=0 THEN
+    v_len := c_sdp_N_sub * c_sdp_N_pol * c_stat_data_sz;  -- 2048 = 512 * 2 * 64/32
+    v_span := true_log_pow2(v_len);                       -- = 2048
+    FOR I IN 0 TO v_len-1 LOOP
+      v_W := I MOD c_stat_data_sz;                   -- 0, 1 per statistics word, word index
+      v_P := (I / c_stat_data_sz) MOD c_sdp_N_pol;   -- 0, 1 per SP pol, polarization index
+      v_B := I / (c_sdp_N_pol * c_stat_data_sz);     -- subband index, range(N_sub = 512) per dual pol
+      v_addr := I + c_pfb_index * v_span;            -- MM address
+      -- Only read SST for g_subband for dual pol SP, to save sim time
+      IF g_read_all_SST = TRUE OR v_B = g_subband THEN
         IF v_W=0 THEN
           -- low part
-          mmf_mm_bus_rd(c_mm_file_ram_st_sst, I, rd_data, tb_clk);
-          sp_subband_powers_arr2(v_S)(v_B)(31 DOWNTO 0) <= rd_data;
+          mmf_mm_bus_rd(c_mm_file_ram_st_sst, v_addr, rd_data, tb_clk);
+          v_data_lo := rd_data;
         ELSE      
           -- high part
-          mmf_mm_bus_rd(c_mm_file_ram_st_sst, I, rd_data, tb_clk);
-          sp_subband_powers_arr2(v_S)(v_B)(63 DOWNTO 32) <= rd_data;
+          mmf_mm_bus_rd(c_mm_file_ram_st_sst, v_addr, rd_data, tb_clk);
+          v_data_hi := rd_data;
+          v_stat_data := v_data_hi & v_data_lo;
+
+          sp_subband_ssts_arr2(v_P)(v_B) <= v_stat_data;
+          stat_data <= v_stat_data;  -- for time series view in Wave window
         END IF;
       END IF;
     END LOOP;
+    proc_common_wait_some_cycles(tb_clk, 1);
+
+    -- Subband power of g_subband in g_sp
+    -- . For the selected g_subband in g_sp the sp_subband_sst will be close
+    --   to sp_subband_sst_sum_arr(c_pol_index), because the input is a
+    --   sinus, so most power will be in 1 subband.
+    sp_subband_sst <= TO_UREAL(sp_subband_ssts_arr2(c_pol_index)(g_subband));
+    proc_common_wait_some_cycles(tb_clk, 1);
+    proc_common_wait_some_cycles(ext_clk, 100);  -- delay for ease of view in Wave window
  
     ---------------------------------------------------------------------------
     -- Read beamlet statistics
     ---------------------------------------------------------------------------
     -- . the beamlet statistics are c_stat_data_sz = 2 word power values.
-    -- . there are c_sdp_S_sub_bf = 488 subbands per signal path
-    -- . the subbands are output alternately so A0 B0 A1 B1 ... A5487 B487 for input A, B
-    -- . the subband statistics multiple units appear in order in the ram_st_bst address map
-    -- . the subband statistics are stored first lo word 0 then hi word 1
-    FOR I IN 0 TO c_sdp_N_pol_bf*c_sdp_S_sub_bf*c_stat_data_sz-1 LOOP
-      v_W := I MOD c_stat_data_sz;
-      v_T := (I / c_stat_data_sz) MOD c_sdp_N_pol_bf;
-      v_U := I / (c_sdp_N_pol_bf*c_stat_data_sz*c_sdp_S_sub_bf);
-      v_S := v_T + v_U * c_sdp_N_pol_bf;
-      v_B := (I / (c_sdp_N_pol_bf*c_stat_data_sz)) MOD c_sdp_S_sub_bf;
-      -- Only read beamset 0, pol 0 (v_S = 0)
-      IF v_S=0 THEN
-        IF v_W=0 THEN
-          -- low part
-          --mmf_mm_bus_rd(c_mm_file_ram_st_bst, I+(c_sdp_N_pol_bf*c_sdp_N_sub*c_stat_data_sz), rd_data, tb_clk);
-          mmf_mm_bus_rd(c_mm_file_ram_st_bst, I, rd_data, tb_clk);
-          sp_beamlet_powers_arr2(v_S)(v_B)(31 DOWNTO 0) <= rd_data;
-        ELSE      
-          -- high part
-          --mmf_mm_bus_rd(c_mm_file_ram_st_bst, I+(c_sdp_N_pol_bf*c_sdp_N_sub*c_stat_data_sz), rd_data, tb_clk);
-          mmf_mm_bus_rd(c_mm_file_ram_st_bst, I, rd_data, tb_clk);
-          sp_beamlet_powers_arr2(v_S)(v_B)(63 DOWNTO 32) <= rd_data;
-  
-          -- Convert STD_LOGIC_VECTOR to REAL
-          v_sp_beamlet_power := REAL(TO_UINT(rd_data(29 DOWNTO 0) & 
-              sp_beamlet_powers_arr2(v_S)(v_B)(31 DOWNTO 30)))*2.0**30 + 
-              REAL(TO_UINT(sp_beamlet_powers_arr2(v_S)(v_B)(29 DOWNTO 0)));
-          -- sum
-          sp_beamlet_power_sum(v_S) <= sp_beamlet_power_sum(v_S) + v_sp_beamlet_power;
+    -- . there are c_sdp_S_sub_bf = 488 dual pol beamlets per beamset
+    -- . the beamlets are output alternately so X0 Y0 X1 Y1 ... X487 Y487 for polarizations X, Y
+    -- . the beamlet statistics for multiple beamsets appear in order in the ram_st_bst address map
+    -- . the beamlet statistics are stored first lo word 0 then hi word 1
+    v_len := c_sdp_S_sub_bf * c_sdp_N_pol_bf * c_stat_data_sz;  -- = 1952 = 488 * 2 * 64/32
+    v_span := true_log_pow2(v_len);                             -- = 2048
+    FOR U IN 0 TO c_sdp_N_beamsets-1 LOOP
+      FOR I IN 0 TO v_len-1 LOOP
+        v_W := I MOD c_stat_data_sz;                     -- 0, 1 per statistics word, word index
+        v_P := (I / c_stat_data_sz) MOD c_sdp_N_pol_bf;  -- 0, 1 per BF pol, polarization index
+        v_B := I / (c_sdp_N_pol_bf * c_stat_data_sz);    -- beamlet index in beamset, range(S_sub_bf = 488) per dual pol
+        v_G := v_B + U * c_sdp_S_sub_bf;                 -- global beamlet index, range(c_sdp_N_beamlets_sdp)
+        v_addr := I + U * v_span;                        -- MM address
+        --Only read BST for g_beamlet and dual pol_bf 0 and 1 and for both beamsets, to save sim time
+        IF g_read_all_BST = TRUE OR v_B = g_beamlet THEN
+          IF v_W = 0 THEN
+            -- low part
+            mmf_mm_bus_rd(c_mm_file_ram_st_bst, v_addr, rd_data, tb_clk);
+            v_data_lo := rd_data;
+          ELSE
+            -- high part
+            mmf_mm_bus_rd(c_mm_file_ram_st_bst, v_addr, rd_data, tb_clk);
+            v_data_hi := rd_data;
+            v_stat_data := v_data_hi & v_data_lo;
+
+            pol_beamlet_bsts_arr2(v_P)(v_G) <= v_stat_data;
+            stat_data <= v_stat_data;  -- for time series view in Wave window
+          END IF;
         END IF;
-      END IF;
+      END LOOP;
     END LOOP;
+    proc_common_wait_some_cycles(tb_clk, 1);
 
-    -- sp_beamlet_power_sum is the sum of all subband powers per SP, this value will be close to sp_beamlet_power
-    -- because the input is a sinus, so most power will be in 1 subband. The sp_beamlet_power_leakage_sum shows
-    -- how much power from the input sinus at a specific subband has leaked into the 511 other subbands.
-    sp_beamlet_power_0 <=
-        REAL(TO_UINT(sp_beamlet_powers_arr2(0)(INTEGER(ROUND(c_subband_sp_0)))(61 DOWNTO 30)))*2.0**30 +
-        REAL(TO_UINT(sp_beamlet_powers_arr2(0)(INTEGER(ROUND(c_subband_sp_0)))(29 DOWNTO 0)));
-
-    sp_beamlet_power_sum_0 <= sp_beamlet_power_sum(0);
-    
+    -- Beamlet power of g_beamlet X and Y, same for both beamsets
+    FOR U IN 0 TO c_sdp_N_beamsets-1 LOOP
+      v_G := g_beamlet + U * c_sdp_S_sub_bf;  -- global beamlet index, range(c_sdp_N_beamlets_sdp)
+      pol_beamlet_bst_X_arr(U) <= TO_UREAL(pol_beamlet_bsts_arr2(0)(v_G));  -- X pol beamlet
+      pol_beamlet_bst_Y_arr(U) <= TO_UREAL(pol_beamlet_bsts_arr2(1)(v_G));  -- Y pol beamlet
+    END LOOP;
     proc_common_wait_some_cycles(tb_clk, 1);
+    proc_common_wait_some_cycles(ext_clk, 100);  -- delay for ease of view in Wave window
  
     ---------------------------------------------------------------------------
-    -- Verify subband statistics
+    -- Log WG, subband and beamlet statistics
     --------------------------------------------------------------------------- 
-    FOR I IN 0 TO c_sdp_N_pol*c_sdp_S_sub_bf-1 LOOP
-      v_T := I  MOD c_sdp_N_pol;
-      v_U := I / (c_sdp_N_pol*c_sdp_S_sub_bf);
-      v_S := v_T + v_U * c_sdp_N_pol;
-      v_B := (I / c_sdp_N_pol) MOD c_sdp_S_sub_bf;
-      IF v_S=0 THEN
-        -- Convert STD_LOGIC_VECTOR to REAL
-        v_sp_beamlet_power := REAL(TO_UINT(rd_data(29 DOWNTO 0) & 
-            sp_beamlet_powers_arr2(v_S)(v_B)(31 DOWNTO 30)))*2.0**30 + 
-            REAL(TO_UINT(sp_beamlet_powers_arr2(v_S)(v_B)(29 DOWNTO 0)));
   
-        -- Convert STD_LOGIC_VECTOR to REAL
-        v_sp_subband_power := REAL(TO_UINT(rd_data(29 DOWNTO 0) & 
-            sp_subband_powers_arr2(v_S)(v_B)(31 DOWNTO 30)))*2.0**30 + 
-            REAL(TO_UINT(sp_subband_powers_arr2(v_S)(v_B)(29 DOWNTO 0)));
+    print_str("");
+    print_str("WG:");
+    print_str(". c_wg_ampl                            = " & int_to_str(c_wg_ampl));
+    print_str(". c_exp_sp_power                       = " & real_to_str(c_exp_sp_power, 20, 1));
+    print_str(". c_exp_sp_ast                         = " & real_to_str(c_exp_sp_ast, 20, 1));
   
-        -- verify if subband power and beamlet power are the same. This is expected because we only use 1 WG input and the BF weights have unit value.
-        -- the difference should not be larger than 0.5% (+/- 2^13 for low values)
-        ASSERT v_sp_beamlet_power > 0.995 * v_sp_subband_power -2.0**13 REPORT "index ("& integer'image(v_S) &","& integer'image(v_B)  &"): Subband power = "& real'image(v_sp_subband_power) &" and Beamlet power = "& real'image(v_sp_beamlet_power) &" are not equal" SEVERITY ERROR;
-        ASSERT v_sp_beamlet_power < 1.005 * v_sp_subband_power +2.0**13 REPORT "index ("& integer'image(v_S) &","& integer'image(v_B)  &"): Subband power = "& real'image(v_sp_subband_power) &" and Beamlet power = "& real'image(v_sp_beamlet_power) &" are not equal" SEVERITY ERROR;
-      END IF;
+    print_str("");
+    print_str("Subband selector:");
+    print_str(". sst_offload_weighted_subbands        = " & sl_to_str(sst_offload_weighted_subbands));
+
+    print_str("");
+    print_str("Subband weight:");
+    print_str(". sp_subband_weight_gain               = " & real_to_str(sp_subband_weight_gain, 20, 6));
+    print_str(". sp_subband_weight_phase              = " & real_to_str(sp_subband_weight_phase, 20, 6));
+
+    print_str("");
+    print_str("SST results:");
+    print_str(". c_exp_subband_ampl                   = " & int_to_str(NATURAL(c_exp_subband_ampl)));
+    print_str(". c_exp_subband_power                  = " & real_to_str(c_exp_subband_power, 20, 1));
+    print_str(". c_exp_subband_sst                    = " & real_to_str(c_exp_subband_sst, 20, 1));
+    print_str("");
+    print_str(". sp_subband_sst                       = " & real_to_str(sp_subband_sst, 20, 1));
+    print_str(". sp_subband_sst / c_exp_subband_sst   = " & real_to_str(sp_subband_sst / c_exp_subband_sst, 20, 6));
+
+    print_str("");
+    print_str("BST results:");
+    print_str(". c_exp_beamlet_ampl                   = " & int_to_str(NATURAL(c_exp_beamlet_ampl)));
+    print_str(". c_exp_beamlet_power                  = " & real_to_str(c_exp_beamlet_power, 20, 1));
+    print_str(". c_exp_beamlet_bst                    = " & real_to_str(c_exp_beamlet_bst, 20, 1));
+    print_str("");
+    FOR U IN 0 TO c_sdp_N_beamsets-1 LOOP
+      v_G := g_beamlet + U * c_sdp_S_sub_bf;  -- global beamlet index, range(c_sdp_N_beamlets_sdp)
+      print_str(". pol_beamlet_bst_X beamlet(" & INTEGER'IMAGE(v_G) & ") = " & real_to_str(pol_beamlet_bst_X_arr(U), 20, 1));
+      print_str(". pol_beamlet_bst_Y beamlet(" & INTEGER'IMAGE(v_G) & ") = " & real_to_str(pol_beamlet_bst_Y_arr(U), 20, 1));
+    END LOOP;
+    FOR U IN 0 TO c_sdp_N_beamsets-1 LOOP
+      v_G := g_beamlet + U * c_sdp_S_sub_bf;  -- global beamlet index, range(c_sdp_N_beamlets_sdp)
+      print_str(". pol_beamlet_bst_X beamlet(" & INTEGER'IMAGE(v_G) & ") / c_exp_beamlet_bst = " & real_to_str(pol_beamlet_bst_X_arr(U) / c_exp_beamlet_bst, 20, 6));
+      print_str(". pol_beamlet_bst_Y beamlet(" & INTEGER'IMAGE(v_G) & ") / c_exp_beamlet_bst = " & real_to_str(pol_beamlet_bst_Y_arr(U) / c_exp_beamlet_bst, 20, 6));
     END LOOP;
 
-    -- verify expected subband power based on WG power
-    IF sp_beamlet_power_sum_0>0.0 THEN ASSERT sp_beamlet_power_0 > c_lo_factor * c_exp_beamlet_power_sp_0 REPORT "Wrong beamlet power for SP 0" SEVERITY ERROR; END IF;
-    IF sp_beamlet_power_sum_0>0.0 THEN ASSERT sp_beamlet_power_0 < c_hi_factor * c_exp_beamlet_power_sp_0 REPORT "Wrong beamlet power for SP 0" SEVERITY ERROR; END IF;
+    print_str("");
+    print_str("Beamlet output:");
+    print_str(". rd_beamlet_scale                     = " & int_to_str(TO_UINT(rd_beamlet_scale)));
+    print_str(". c_exp_beamlet_output_ampl            = " & int_to_str(NATURAL(c_exp_beamlet_output_ampl)));
     
-    -- view c_exp_sp_beamlet_power_ratio in Wave window
-    IF sp_beamlet_power_sum_0>0.0 THEN sp_beamlet_power_ratio_0 <= sp_beamlet_power_0/sp_beamlet_power_sum_0; END IF;
+    ---------------------------------------------------------------------------
+    -- Verify SST and BST
+    ---------------------------------------------------------------------------
     
-    -- view c_exp_sp_beamlet_power_sum_ratio in Wave window
-    -- The sp_beamlet_power_sum_ratio show similar information as sp_beamlet_power_leakage_sum, because when
-    -- sp_beamlet_power_leakage_sum is small then sp_beamlet_power_sum_ratio ~= sp_beamlet_power_ratio.
-    IF sp_beamlet_power_sum_0>0.0 THEN sp_beamlet_power_sum_ratio_0 <= sp_beamlet_power_sum_0/sp_beamlet_power_0; END IF;
-
-    -- View sp_beamlet_power_leakage_sum in Wave window
-    IF sp_beamlet_power_sum_0>0.0 THEN sp_beamlet_power_leakage_sum_0 <= sp_beamlet_power_sum_0 - sp_beamlet_power_0; END IF;
+    -- verify expected subband power based on WG power
+    ASSERT sp_subband_sst > c_stat_lo_factor * c_exp_subband_sst REPORT "Wrong subband power for SP " & NATURAL'IMAGE(g_sp) SEVERITY ERROR;
+    ASSERT sp_subband_sst < c_stat_hi_factor * c_exp_subband_sst REPORT "Wrong subband power for SP " & NATURAL'IMAGE(g_sp) SEVERITY ERROR;
+
+    -- verify expected beamlet power based on WG power and BF weigths
+    --
+    -- All co and cross polarization weights are equal: w = w_xx = w_xy = w_yx = w_yy
+    -- With one g_sp either SP X or SP Y is used, so the other antenna polarization is 0
+    -- Hence the c_exp_beamlet_bst will be the same for beamlet X and Y independent of whether g_sp is an X or Y signal input:
+    --   g_sp = X --> w_xx --> beamlet X = c_exp_beamlet_bst
+    --   g_sp = Y --> w_xy --> beamlet X = c_exp_beamlet_bst
+    --   g_sp = X --> w_yx --> beamlet Y = c_exp_beamlet_bst
+    --   g_sp = Y --> w_yy --> beamlet Y = c_exp_beamlet_bst
+    --
+    FOR U IN 0 TO c_sdp_N_beamsets-1 LOOP
+      ASSERT pol_beamlet_bst_X_arr(U) > c_stat_lo_factor * c_exp_beamlet_bst REPORT "Wrong beamlet power for X in beamset " & NATURAL'IMAGE(U) SEVERITY ERROR;
+      ASSERT pol_beamlet_bst_X_arr(U) < c_stat_hi_factor * c_exp_beamlet_bst REPORT "Wrong beamlet power for X in beamset " & NATURAL'IMAGE(U) SEVERITY ERROR;
+      ASSERT pol_beamlet_bst_Y_arr(U) > c_stat_lo_factor * c_exp_beamlet_bst REPORT "Wrong beamlet power for Y in beamset " & NATURAL'IMAGE(U) SEVERITY ERROR;
+      ASSERT pol_beamlet_bst_Y_arr(U) < c_stat_hi_factor * c_exp_beamlet_bst REPORT "Wrong beamlet power for Y in beamset " & NATURAL'IMAGE(U) SEVERITY ERROR;
+    END LOOP;
 
     ---------------------------------------------------------------------------
-    -- Verify 10GbE UDP offload
+    -- Verify beamlet output in 10GbE UDP offload
     ---------------------------------------------------------------------------
-    v_re := TO_SINT(rx_beamlet_list_re(c_exp_beamlet_index));  v_re_exp := TO_SINT(c_exp_beamlet_re);
-    v_im := TO_SINT(rx_beamlet_list_im(c_exp_beamlet_index));  v_im_exp := TO_SINT(c_exp_beamlet_im);
-    ASSERT v_re = v_re_exp REPORT "Wrong 10GbE output (re) " & INTEGER'IMAGE(v_re) & " != " & INTEGER'IMAGE(v_re_exp) SEVERITY ERROR;
-    ASSERT v_im = v_im_exp REPORT "Wrong 10GbE output (im) " & INTEGER'IMAGE(v_re) & " != " & INTEGER'IMAGE(v_re_exp) SEVERITY ERROR;
+    v_re := TO_SINT(rx_beamlet_list_re(c_exp_beamlet_index)); v_re_exp := c_exp_beamlet_output_re;
+    v_im := TO_SINT(rx_beamlet_list_im(c_exp_beamlet_index)); v_im_exp := c_exp_beamlet_output_im;
+    ASSERT v_re > INTEGER(v_re_exp) - c_beamlet_output_delta REPORT "Wrong 10GbE output (re) " & INTEGER'IMAGE(v_re) & " != " & REAL'IMAGE(v_re_exp) SEVERITY ERROR;
+    ASSERT v_re < INTEGER(v_re_exp) + c_beamlet_output_delta REPORT "Wrong 10GbE output (re) " & INTEGER'IMAGE(v_re) & " != " & REAL'IMAGE(v_re_exp) SEVERITY ERROR;
+    ASSERT v_im > INTEGER(v_im_exp) - c_beamlet_output_delta REPORT "Wrong 10GbE output (im) " & INTEGER'IMAGE(v_im) & " != " & REAL'IMAGE(v_im_exp) SEVERITY ERROR;
+    ASSERT v_im < INTEGER(v_im_exp) + c_beamlet_output_delta REPORT "Wrong 10GbE output (im) " & INTEGER'IMAGE(v_im) & " != " & REAL'IMAGE(v_im_exp) SEVERITY ERROR;
 
     ---------------------------------------------------------------------------
     -- End Simulation 
     ---------------------------------------------------------------------------   
     tb_almost_end <= '1';
-    proc_common_wait_some_cycles(ext_clk, 100);
+    proc_common_wait_some_cycles(ext_clk, 100);  -- delay for ease of view in Wave window
     proc_common_stop_simulation(TRUE, ext_clk, tb_almost_end, tb_end);
     WAIT;
   END PROCESS;
@@ -852,7 +1099,7 @@ BEGIN
     FOR I IN 0 TO (c_sdp_cep_nof_blocks_per_packet * c_sdp_cep_nof_beamlets_per_block/2)-1 LOOP
       proc_common_wait_until_high(ext_clk, test_offload_sosi.valid);
       rx_beamlet_valid <= '1';
-      -- Capture rx beamlets per longword in rx_beamlet_arr
+      -- Capture rx beamlets per longword in rx_beamlet_arr, for time series view in Wave window
       rx_beamlet_arr_re(0) <= test_offload_sosi.data(55 DOWNTO 48);  -- X
       rx_beamlet_arr_im(0) <= test_offload_sosi.data(63 DOWNTO 56);
       rx_beamlet_arr_re(1) <= test_offload_sosi.data(39 DOWNTO 32);  -- Y
@@ -863,13 +1110,13 @@ BEGIN
       rx_beamlet_arr_im(3) <= test_offload_sosi.data(15 DOWNTO 8);
       IF I < c_sdp_cep_nof_beamlets_per_block/2 THEN
         -- Only capture the first beamlets block of each packet in rx_beamlet_list
-        rx_beamlet_list_re(I*4 + 0) <= test_offload_sosi.data(55 DOWNTO 48);
+        rx_beamlet_list_re(I*4 + 0) <= test_offload_sosi.data(55 DOWNTO 48);  -- X
         rx_beamlet_list_im(I*4 + 0) <= test_offload_sosi.data(63 DOWNTO 56);
-        rx_beamlet_list_re(I*4 + 1) <= test_offload_sosi.data(39 DOWNTO 32);
+        rx_beamlet_list_re(I*4 + 1) <= test_offload_sosi.data(39 DOWNTO 32);  -- Y
         rx_beamlet_list_im(I*4 + 1) <= test_offload_sosi.data(47 DOWNTO 40);
-        rx_beamlet_list_re(I*4 + 2) <= test_offload_sosi.data(23 DOWNTO 16);
+        rx_beamlet_list_re(I*4 + 2) <= test_offload_sosi.data(23 DOWNTO 16);  -- X
         rx_beamlet_list_im(I*4 + 2) <= test_offload_sosi.data(31 DOWNTO 24);
-        rx_beamlet_list_re(I*4 + 3) <= test_offload_sosi.data( 7 DOWNTO 0);
+        rx_beamlet_list_re(I*4 + 3) <= test_offload_sosi.data( 7 DOWNTO 0);   -- Y
         rx_beamlet_list_im(I*4 + 3) <= test_offload_sosi.data(15 DOWNTO 8);
       END IF;
       proc_common_wait_until_high(ext_clk, test_offload_sosi.valid);
diff --git a/applications/lofar2/designs/lofar2_unb2c_sdp_station/revisions/lofar2_unb2c_sdp_station_fsub/tb_lofar2_unb2c_sdp_station_fsub.vhd b/applications/lofar2/designs/lofar2_unb2c_sdp_station/revisions/lofar2_unb2c_sdp_station_fsub/tb_lofar2_unb2c_sdp_station_fsub.vhd
index c28b78463e1570614b135e582e07f9291283c770..ef550d46f766c8083a10c39f05cfa4133039c085 100644
--- a/applications/lofar2/designs/lofar2_unb2c_sdp_station/revisions/lofar2_unb2c_sdp_station_fsub/tb_lofar2_unb2c_sdp_station_fsub.vhd
+++ b/applications/lofar2/designs/lofar2_unb2c_sdp_station/revisions/lofar2_unb2c_sdp_station_fsub/tb_lofar2_unb2c_sdp_station_fsub.vhd
@@ -20,7 +20,7 @@
 
 -------------------------------------------------------------------------------
 --
--- Author: R. van der Walle
+-- Author: R. van der Walle (original), E. Kooistra (updates)
 -- Purpose: Self-checking testbench for simulating lofar2_unb2c_sdp_station_fsub using WG data.
 --
 -- Description:
@@ -29,44 +29,68 @@
 --   1) Enable calc mode for WG via reg_diag_wg with:
 --        freq = 19.921875MHz
 --        ampl = 0.5 * 2**13
---   
---   2) Read current BSN from reg_bsn_scheduler_wg and write reg_bsn_scheduler_wg 
+--
+--   2) Read current BSN from reg_bsn_scheduler_wg and write reg_bsn_scheduler_wg
 --      to trigger start of WG at BSN.
---     
---   3) Read subband statistics (SST) via ram_st_sst and verify with 
---      c_exp_subband_power_sp_0 at c_subband_sp_0.
---      View sp_subband_power_0  in Wave window
---   
+--
+--   3) Read subband statistics (SST) via MM and verify with exp_subband_sst at g_subband.
+--      . use weighted subbands (default selected by MM)
+--
+--   4) View in wave window
+--      . in_sosi.sop and in_data in u_si_arr(g_sp) to check that:
+--        - WG starts with zero phase sine when c_subband_phase = 0.0 degrees
+--        - WG amplitude = 8191 (is full scale - 1) when wg_ampl = 1.0
+--      . pfb_sosi_arr(c_pfb_index).im/re and fsub_sosi_arr(c_pfb_index).im/re
+--        in u_fsub in decimal radix and analog format to check that subband
+--        phase is g_subband_weight_phase phase as set by the subband weight.
+--        - Raw: pfb_sosi_arr = atan2(-65195 / 0) = -90 degrees
+--        - Weighted: fsub_sosi_arr = atan2(-56457 / 32598) = -60 degrees
+--          --> rotated expected g_subband_weight_phase = -60 - -90 = +30 degrees.
 --
 -- Usage:
 --   > as 7    # default
 --   > as 12   # for detailed debugging
---   > run -a  
+--   # Manually add missing signal
+--   > add wave -position insertpoint  \
+--     sim:/tb_lofar2_unb2c_sdp_station_fsub/sp_subband_ssts_arr2
+--   > run -a
+--   # Takes about 30 m when g_read_all_SST = FALSE
+--   # Takes about 40 m when g_read_all_SST = TRUE
 --
 -------------------------------------------------------------------------------
 LIBRARY IEEE, common_lib, unb2c_board_lib, i2c_lib, mm_lib, dp_lib, diag_lib, lofar2_sdp_lib, wpfb_lib, lofar2_unb2c_sdp_station_lib;
 USE IEEE.std_logic_1164.ALL;
 USE IEEE.numeric_std.ALL;
-USE IEEE.MATH_REAL.ALL;
+USE IEEE.math_real.ALL;
 USE common_lib.common_pkg.ALL;
-USE unb2c_board_lib.unb2c_board_pkg.ALL;
 USE common_lib.tb_common_pkg.ALL;
 USE common_lib.common_str_pkg.ALL;
 USE mm_lib.mm_file_pkg.ALL;
-USE dp_lib.dp_stream_pkg.ALL;
 USE mm_lib.mm_file_unb_pkg.ALL;
+USE dp_lib.dp_stream_pkg.ALL;
 USE diag_lib.diag_pkg.ALL;
 USE wpfb_lib.wpfb_pkg.ALL;
 USE lofar2_sdp_lib.sdp_pkg.ALL;
+USE unb2c_board_lib.unb2c_board_pkg.ALL;
 
 ENTITY tb_lofar2_unb2c_sdp_station_fsub IS
+  GENERIC (
+    g_sp                   : NATURAL := 3;    -- signal path index in range(S_pn = 12)
+    g_wg_ampl              : REAL := 1.0;     -- WG normalized amplitude
+    g_subband              : NATURAL := 102;  -- select subband at index 102 = 102/1024 * 200MHz = 19.921875 MHz
+    g_subband_weight_gain  : REAL := 1.0;     -- subband weight normalized gain
+    g_subband_weight_phase : REAL := 30.0;    -- subband weight phase rotation in degrees
+    g_read_all_SST         : BOOLEAN := TRUE  -- when FALSE only read SST for g_subband, to save sim time
+  );
 END tb_lofar2_unb2c_sdp_station_fsub;
 
 ARCHITECTURE tb OF tb_lofar2_unb2c_sdp_station_fsub IS
 
   CONSTANT c_sim             : BOOLEAN := TRUE;
   CONSTANT c_unb_nr          : NATURAL := 0; -- UniBoard 0
-  CONSTANT c_node_nr         : NATURAL := 0; 
+  CONSTANT c_node_nr         : NATURAL := 0;
+  CONSTANT c_init_bsn        : NATURAL := 17;  -- some recognizable value >= 0
+
   CONSTANT c_id              : STD_LOGIC_VECTOR(7 DOWNTO 0) := "00000000";
   CONSTANT c_version         : STD_LOGIC_VECTOR(1 DOWNTO 0) := "00";
   CONSTANT c_fw_version      : t_unb2c_board_fw_version := (1, 0);
@@ -78,63 +102,108 @@ ARCHITECTURE tb OF tb_lofar2_unb2c_sdp_station_fsub IS
   CONSTANT c_tb_clk_period       : TIME := 100 ps; -- use fast tb_clk to speed up M&C
 
   CONSTANT c_nof_block_per_sync  : NATURAL := 16;
-  CONSTANT c_nof_clk_per_sync    : NATURAL := c_nof_block_per_sync*c_sdp_N_fft - (c_sdp_N_fft/2); --15.5 block per sync
+  CONSTANT c_nof_clk_per_sync    : NATURAL := c_nof_block_per_sync * c_sdp_N_fft;
   CONSTANT c_pps_period          : NATURAL := c_nof_clk_per_sync;
   CONSTANT c_wpfb_sim            : t_wpfb := func_wpfb_set_nof_block_per_sync(c_sdp_wpfb_subbands, c_nof_block_per_sync);
-   
+  CONSTANT c_stat_data_sz        : NATURAL := c_wpfb_sim.stat_data_sz;  -- = 2
+
   CONSTANT c_percentage          : REAL := 0.05;  -- percentage that actual value may differ from expected value
-  CONSTANT c_lo_factor           : REAL := 1.0 - c_percentage;  -- lower boundary  
+  CONSTANT c_lo_factor           : REAL := 1.0 - c_percentage;  -- lower boundary
   CONSTANT c_hi_factor           : REAL := 1.0 + c_percentage;  -- higher boundary
 
   -- WG
-  CONSTANT c_full_scale_ampl      : REAL := REAL(2**(14-1)-1);  -- = full scale of WG
-  CONSTANT c_bsn_start_wg         : NATURAL := 2;  -- start WG at this BSN to instead of some BSN, to avoid mismatches in exact expected data values
-  CONSTANT c_ampl_sp_0            : NATURAL := 2**(c_sdp_W_adc-1)/2;  -- in number of lsb
+  CONSTANT c_bsn_start_wg         : NATURAL := c_init_bsn + 2;  -- start WG at this BSN to instead of some BSN, to avoid mismatches in exact expected data values
+  -- .ampl
+  CONSTANT c_wg_ampl_full_scale   : NATURAL := 2**(c_sdp_W_adc-1);  -- full scale (FS) of WG, will just cause clipping of +FS to +FS-1
+  CONSTANT c_wg_ampl_lsb          : REAL := c_diag_wg_ampl_unit / REAL(c_wg_ampl_full_scale);  -- amplitude in number of LSbit resolution steps
+  CONSTANT c_wg_ampl              : NATURAL := NATURAL(g_wg_ampl * REAL(c_wg_ampl_full_scale));  -- in number of lsb
+  CONSTANT c_exp_sp_power         : REAL := REAL(c_wg_ampl**2) / 2.0;
+  CONSTANT c_exp_sp_ast           : REAL := c_exp_sp_power * REAL(c_nof_clk_per_sync);
+  -- . phase
+  CONSTANT c_subband_phase        : REAL := 0.0;  -- wanted subband phase in degrees = WG phase at sop
+  CONSTANT c_subband_freq         : REAL := REAL(g_subband) / REAL(c_sdp_N_fft);                   -- normalized by fs = f_adc = 200 MHz is dp_clk
+  CONSTANT c_wg_latency           : INTEGER := c_diag_wg_latency - 0;  -- -0 to account for BSN scheduler start trigger latency
+  CONSTANT c_wg_phase_offset      : REAL := 360.0 * REAL(c_wg_latency) * c_subband_freq;  -- c_diag_wg_latency is in dp_clk cycles
+  CONSTANT c_wg_phase             : REAL := c_subband_phase + c_wg_phase_offset;  -- WG phase in degrees
+  -- . freq
   CONSTANT c_wg_subband_freq_unit : REAL := c_diag_wg_freq_unit/REAL(c_sdp_N_fft);  -- subband freq = Fs/1024 = 200 MSps/1024 = 195312.5 Hz sinus
-  CONSTANT c_wg_freq_offset       : REAL := 0.0/11.0; -- in freq_unit
-  CONSTANT c_subband_sp_0         : REAL := 102.0;  -- Select subband at index 102 = 102/1024 * 200MHz = 19.921875 MHz 
-  CONSTANT c_wg_ampl_lsb          : REAL := c_diag_wg_ampl_unit / c_full_scale_ampl;  -- amplitude in number of LSbit resolution steps
-  CONSTANT c_exp_wg_power_sp_0    : REAL := REAL(c_ampl_sp_0**2)/2.0 * REAL(c_nof_clk_per_sync);
-
-  -- WPFB
-  CONSTANT c_nof_pfb                        : NATURAL := 1; -- Verifying 1 of c_sdp_P_pfb = 6 pfb to speed up simulation.
-  CONSTANT c_wb_leakage_bin                 : NATURAL := c_wpfb_sim.nof_points / c_wpfb_sim.wb_factor;   -- = 256, leakage will occur in this bin if FIR wb_factor is reversed 
-  CONSTANT c_exp_sp_subband_power_ratio     : REAL := 1.0/8.0;   -- depends on internal WPFB quantization and FIR coefficients
-  CONSTANT c_exp_sp_subband_power_sum_ratio : REAL := c_exp_sp_subband_power_ratio;   -- because all sinus power is expected in one subband
-  CONSTANT c_exp_subband_power_sp_0         : REAL := c_exp_wg_power_sp_0 * c_exp_sp_subband_power_ratio;
-
-  TYPE t_real_arr IS ARRAY (INTEGER RANGE <>) OF REAL; 
-  TyPE t_slv_64_subbands_arr IS ARRAY (INTEGER RANGE <>) OF t_slv_64_arr(0 TO c_sdp_N_sub);
-
-  -- MM  
+
+  -- FSUB
+  -- . WPFB
+  CONSTANT c_pol_index                      : NATURAL := g_sp MOD c_sdp_Q_fft;
+  CONSTANT c_pfb_index                      : NATURAL := g_sp / c_sdp_Q_fft;  -- only read used WPFB unit out of range(c_sdp_P_pfb = 6)
+  CONSTANT c_exp_subband_sp_ampl_ratio      : REAL := 7.96;  -- ~= 8 for unit FIR DC gain, depends on internal WPFB quantization and FIR coefficients
+  CONSTANT c_exp_subband_ampl_raw           : REAL := REAL(c_wg_ampl) * c_exp_subband_sp_ampl_ratio;
+  CONSTANT c_exp_subband_ampl_weighted      : REAL := c_exp_subband_ampl_raw * g_subband_weight_gain;
+  CONSTANT c_exp_subband_power_raw          : REAL := c_exp_subband_ampl_raw**2.0;  -- complex, so no divide by 2
+  CONSTANT c_exp_subband_power_weighted     : REAL := c_exp_subband_ampl_weighted**2.0;  -- complex, so no divide by 2
+  CONSTANT c_exp_subband_sst_raw            : REAL := c_exp_subband_power_raw * REAL(c_nof_block_per_sync);
+  CONSTANT c_exp_subband_sst_weighted       : REAL := c_exp_subband_power_weighted * REAL(c_nof_block_per_sync);
+
+  -- . expected limit values, obtained with print_str() for g_subband = 102,
+  --   g_wg_ampl = 1.0, g_subband_weight_gain = 1.0, g_subband_weight_phase = 30.0
+  CONSTANT c_exp_subband_sst_leakage_snr_dB   : REAL := 70.0;  -- < 74.913
+  CONSTANT c_exp_subband_sst_crosstalk_snr_dB : REAL := 90.0;  -- < 96.284
+
+  TYPE t_real_arr IS ARRAY (INTEGER RANGE <>) OF REAL;
+  TyPE t_slv_64_subbands_arr IS ARRAY (INTEGER RANGE <>) OF t_slv_64_arr(0 TO c_sdp_N_sub-1);
+
+  -- . Subband weights for selected g_sp
+  CONSTANT c_subband_weight_re   : INTEGER := INTEGER(g_subband_weight_gain * REAL(c_sdp_unit_sub_weight) * COS(g_subband_weight_phase * MATH_2_PI / 360.0));
+  CONSTANT c_subband_weight_im   : INTEGER := INTEGER(g_subband_weight_gain * REAL(c_sdp_unit_sub_weight) * SIN(g_subband_weight_phase * MATH_2_PI / 360.0));
+
+  -- MM
+  -- . Address widths of a single MM instance
+  CONSTANT c_addr_w_reg_diag_wg           : NATURAL := 2;
+  -- . Address spans of a single MM instance
+  CONSTANT c_mm_span_reg_diag_wg          : NATURAL := 2**c_addr_w_reg_diag_wg;
+
   CONSTANT c_mm_file_reg_bsn_source_v2    : STRING := mmf_unb_file_prefix(c_unb_nr, c_node_nr) & "REG_BSN_SOURCE_V2";
   CONSTANT c_mm_file_reg_bsn_scheduler_wg : STRING := mmf_unb_file_prefix(c_unb_nr, c_node_nr) & "REG_BSN_SCHEDULER";
   CONSTANT c_mm_file_reg_diag_wg          : STRING := mmf_unb_file_prefix(c_unb_nr, c_node_nr) & "REG_WG";
+  CONSTANT c_mm_file_ram_equalizer_gains  : STRING := mmf_unb_file_prefix(c_unb_nr, c_node_nr) & "RAM_EQUALIZER_GAINS";
+  CONSTANT c_mm_file_reg_dp_selector      : STRING := mmf_unb_file_prefix(c_unb_nr, c_node_nr) & "REG_DP_SELECTOR";
   CONSTANT c_mm_file_ram_st_sst           : STRING := mmf_unb_file_prefix(c_unb_nr, c_node_nr) & "RAM_ST_SST";
 
   -- Tb
   SIGNAL tb_end              : STD_LOGIC := '0';
   SIGNAL sim_done            : STD_LOGIC := '0';
-  SIGNAL tb_clk              : STD_LOGIC := '0';  
+  SIGNAL tb_clk              : STD_LOGIC := '0';
   SIGNAL rd_data             : STD_LOGIC_VECTOR(c_32-1 DOWNTO 0) := (OTHERS => '0');
 
+  SIGNAL pps_rst             : STD_LOGIC := '1';
+  SIGNAL gen_pps             : STD_LOGIC := '0';
+
   -- WG
-  SIGNAL current_bsn_wg          : STD_LOGIC_VECTOR(c_dp_stream_bsn_w-1 DOWNTO 0);
-
-  -- WPFB
-  SIGNAL sp_subband_powers_arr2         : t_slv_64_subbands_arr(c_nof_pfb*c_nof_complex-1 DOWNTO 0);   -- [sp][sub]
-  SIGNAL sp_subband_power_0             : REAL;
-  SIGNAL sp_subband_power_sum           : t_real_arr(c_nof_pfb*c_nof_complex-1 DOWNTO 0) := (OTHERS=>0.0);
-  SIGNAL sp_subband_power_sum_0         : REAL;
-  SIGNAL sp_subband_power_ratio_0       : REAL;
-  SIGNAL sp_subband_power_sum_ratio_0   : REAL;
-  SIGNAL sp_subband_power_leakage_sum_0 : REAL;
-  
+  SIGNAL current_bsn_wg      : STD_LOGIC_VECTOR(c_dp_stream_bsn_w-1 DOWNTO 0);
+
+  -- FSUB
+  -- . WPFB
+  SIGNAL sp_subband_ssts_arr2            : t_slv_64_subbands_arr(c_sdp_N_pol-1 DOWNTO 0);   -- [pol][sub]
+  SIGNAL sp_subband_sst_sum_arr          : t_real_arr(c_sdp_N_pol-1 DOWNTO 0) := (OTHERS => 0.0);
+  SIGNAL sp_subband_sst                  : REAL := 0.0;
+  SIGNAL sp_subband_sst_leakage          : REAL := 0.0;
+  SIGNAL sp_subband_sst_leakage_snr_dB   : REAL := 0.0;  -- signal to noise (leakage) ratio
+  SIGNAL sp_subband_sst_crosstalk        : REAL := 0.0;
+  SIGNAL sp_subband_sst_crosstalk_snr_dB : REAL := 0.0;  -- signal to noise (crosstalk) ration
+
+  SIGNAL exp_subband_ampl    : REAL := 0.0;
+  SIGNAL exp_subband_power   : REAL := 0.0;
+  SIGNAL exp_subband_sst     : REAL := 0.0;
+  SIGNAL stat_data           : STD_LOGIC_VECTOR(c_longword_w-1 DOWNTO 0);
+
+  -- . Selector
+  SIGNAL sst_offload_weighted_subbands     : STD_LOGIC;
+
+  -- . Subband equalizer
+  SIGNAL sp_subband_weight_re    : INTEGER := 0;
+  SIGNAL sp_subband_weight_im    : INTEGER := 0;
+  SIGNAL sp_subband_weight_gain  : REAL := 0.0;
+  SIGNAL sp_subband_weight_phase : REAL := 0.0;
+
   -- DUT
   SIGNAL ext_clk             : STD_LOGIC := '0';
-  SIGNAL pps                 : STD_LOGIC := '0';
-  SIGNAL ext_pps             : STD_LOGIC := '0'; 
-  SIGNAL pps_rst             : STD_LOGIC := '0';
+  SIGNAL ext_pps             : STD_LOGIC := '0';
 
   SIGNAL WDI                 : STD_LOGIC;
   SIGNAL INTA                : STD_LOGIC;
@@ -158,7 +227,6 @@ BEGIN
   -- System setup
   ----------------------------------------------------------------------------
   ext_clk <= NOT ext_clk AFTER c_ext_clk_period/2;  -- External clock (200 MHz)
-  eth_clk(0) <= NOT eth_clk(0) AFTER c_eth_clk_period/2;  -- Ethernet ref clock (125 MHz)
   JESD204B_REFCLK <= NOT JESD204B_REFCLK AFTER c_bck_ref_clk_period/2;  -- JESD sample clock (200MHz) 
 
   INTA <= 'H';  -- pull up
@@ -167,9 +235,9 @@ BEGIN
   ------------------------------------------------------------------------------
   -- External PPS
   ------------------------------------------------------------------------------  
-  proc_common_gen_pulse(5, c_pps_period, '1', pps_rst, ext_clk, pps);
-  jesd204b_sysref <= pps;
-  ext_pps <= pps;
+  proc_common_gen_pulse(5, c_pps_period, '1', pps_rst, ext_clk, gen_pps);
+  jesd204b_sysref <= gen_pps;
+  ext_pps <= gen_pps;
 
   ------------------------------------------------------------------------------
   -- DUT
@@ -183,12 +251,12 @@ BEGIN
     g_sim_node_nr            => c_node_nr,
     g_wpfb                   => c_wpfb_sim,
     g_bsn_nof_clk_per_sync   => c_nof_clk_per_sync,
-    g_scope_selected_subband => NATURAL(c_subband_sp_0)
+    g_scope_selected_subband => g_subband
   )
   PORT MAP (
     -- GENERAL
     CLK          => ext_clk,
-    PPS          => pps,
+    PPS          => ext_pps,
     WDI          => WDI,
     INTA         => INTA,
     INTB         => INTB,
@@ -215,125 +283,231 @@ BEGIN
     JESD204B_SYNC_N => jesd204b_sync_n
   );
 
+  -- Raw or weighted subbands
+  exp_subband_ampl  <= sel_a_b(sst_offload_weighted_subbands = '0', c_exp_subband_ampl_raw, c_exp_subband_ampl_weighted);
+  exp_subband_power <= sel_a_b(sst_offload_weighted_subbands = '0', c_exp_subband_power_raw, c_exp_subband_power_weighted);
+  exp_subband_sst   <= sel_a_b(sst_offload_weighted_subbands = '0', c_exp_subband_sst_raw, c_exp_subband_sst_weighted);
+
   ------------------------------------------------------------------------------
   -- MM slave accesses via file IO
   ------------------------------------------------------------------------------
-  tb_clk  <= NOT tb_clk AFTER c_tb_clk_period/2;    -- Testbench MM clock
-  
+  tb_clk <= NOT tb_clk AFTER c_tb_clk_period/2;    -- Testbench MM clock
+
   p_mm_stimuli : PROCESS
-    VARIABLE v_bsn                   : NATURAL;
-    VARIABLE v_sp_subband_power      : REAL;
-    VARIABLE v_W, v_T, v_U, v_S, v_B : NATURAL;  -- array indicies
+    VARIABLE v_bsn                           : NATURAL;
+    VARIABLE v_data_lo, v_data_hi            : STD_LOGIC_VECTOR(c_word_w-1 DOWNTO 0);
+    VARIABLE v_stat_data                     : STD_LOGIC_VECTOR(c_longword_w-1 DOWNTO 0);
+    VARIABLE v_len, v_span, v_offset, v_addr : NATURAL;  -- address ranges, indices
+    VARIABLE v_W, v_P, v_U, v_S, v_B         : NATURAL;  -- array indicies
+    VARIABLE v_re, v_im, v_weight            : INTEGER;
+    VARIABLE v_power                         : REAL;
   BEGIN
     -- Wait for DUT power up after reset
     WAIT FOR 1 us;
 
-    -- wait for pps
-    proc_common_wait_until_hi_lo(ext_clk, ext_pps);
- 
     ----------------------------------------------------------------------------
     -- Enable BSN
     ----------------------------------------------------------------------------
-    mmf_mm_bus_wr(c_mm_file_reg_bsn_source_v2, 3,                    0, tb_clk);
-    mmf_mm_bus_wr(c_mm_file_reg_bsn_source_v2, 2,                    0, tb_clk);  -- Init BSN = 0
-    mmf_mm_bus_wr(c_mm_file_reg_bsn_source_v2, 1,   c_nof_clk_per_sync, tb_clk);  -- nof_block_per_sync
-    mmf_mm_bus_wr(c_mm_file_reg_bsn_source_v2, 0,         16#00000003#, tb_clk);  -- Enable BSN at PPS
-    
+    mmf_mm_bus_wr(c_mm_file_reg_bsn_source_v2, 2,         c_init_bsn, tb_clk);  -- Init BSN
+    mmf_mm_bus_wr(c_mm_file_reg_bsn_source_v2, 3,                  0, tb_clk);  -- Write high part activates the init BSN
+    mmf_mm_bus_wr(c_mm_file_reg_bsn_source_v2, 1, c_nof_clk_per_sync, tb_clk);  -- nof_block_per_sync
+    mmf_mm_bus_wr(c_mm_file_reg_bsn_source_v2, 0,       16#00000003#, tb_clk);  -- Enable BSN at PPS
+
+    -- Release PPS pulser, to get first PPS now and to start BSN source
+    WAIT FOR 1 us;
+    pps_rst <= '0';
+
+    ----------------------------------------------------------------------------
+    -- Read weighted subband selector
+    ----------------------------------------------------------------------------
+    mmf_mm_bus_rd(c_mm_file_reg_dp_selector, 0, rd_data, tb_clk);
+    proc_common_wait_some_cycles(tb_clk, 1);
+    sst_offload_weighted_subbands <= NOT rd_data(0);
+
     ----------------------------------------------------------------------------
-    -- Enable WG
+    -- Enable and start WG
     ----------------------------------------------------------------------------
     --   0 : mode[7:0]           --> off=0, calc=1, repeat=2, single=3)
     --       nof_samples[31:16]  --> <= c_ram_wg_size=1024
     --   1 : phase[15:0]
     --   2 : freq[30:0]
     --   3 : ampl[16:0]
-    mmf_mm_bus_wr(c_mm_file_reg_diag_wg, 0, 1024*2**16 + 1, tb_clk);  -- nof_samples, mode calc
-    mmf_mm_bus_wr(c_mm_file_reg_diag_wg, 1, INTEGER(  0.0 * c_diag_wg_phase_unit), tb_clk);  -- phase offset in degrees
-    mmf_mm_bus_wr(c_mm_file_reg_diag_wg, 2, INTEGER((c_subband_sp_0+c_wg_freq_offset) * c_wg_subband_freq_unit), tb_clk);  -- freq
-    mmf_mm_bus_wr(c_mm_file_reg_diag_wg, 3, INTEGER(REAL(c_ampl_sp_0) * c_wg_ampl_lsb), tb_clk);  -- ampl
+    v_offset := g_sp * c_mm_span_reg_diag_wg;
+    mmf_mm_bus_wr(c_mm_file_reg_diag_wg, v_offset + 0, 1024*2**16 + 1, tb_clk);  -- nof_samples, mode calc
+    mmf_mm_bus_wr(c_mm_file_reg_diag_wg, v_offset + 1, INTEGER(c_wg_phase * c_diag_wg_phase_unit), tb_clk);  -- phase offset in degrees
+    mmf_mm_bus_wr(c_mm_file_reg_diag_wg, v_offset + 2, INTEGER(REAL(g_subband) * c_wg_subband_freq_unit), tb_clk);  -- freq
+    mmf_mm_bus_wr(c_mm_file_reg_diag_wg, v_offset + 3, INTEGER(REAL(c_wg_ampl) * c_wg_ampl_lsb), tb_clk);  -- ampl
 
     -- Read current BSN
     mmf_mm_bus_rd(c_mm_file_reg_bsn_scheduler_wg, 0, current_bsn_wg(31 DOWNTO  0), tb_clk);
     mmf_mm_bus_rd(c_mm_file_reg_bsn_scheduler_wg, 1, current_bsn_wg(63 DOWNTO 32), tb_clk);
     proc_common_wait_some_cycles(tb_clk, 1);
-    
+
     -- Write scheduler BSN to trigger start of WG at next block
     v_bsn := TO_UINT(current_bsn_wg) + 2;
     ASSERT v_bsn <= c_bsn_start_wg REPORT "Too late to start WG: " & int_to_str(v_bsn) & " > " & int_to_str(c_bsn_start_wg) SEVERITY ERROR;
     mmf_mm_bus_wr(c_mm_file_reg_bsn_scheduler_wg, 0, c_bsn_start_wg, tb_clk);  -- first write low then high part
-    mmf_mm_bus_wr(c_mm_file_reg_bsn_scheduler_wg, 1,     0, tb_clk);  -- assume v_bsn < 2**31-1
+    mmf_mm_bus_wr(c_mm_file_reg_bsn_scheduler_wg, 1,              0, tb_clk);  -- assume v_bsn < 2**31-1
 
+    ----------------------------------------------------------------------------
+    -- Write subband weight for selected g_sp and g_subband
+    ----------------------------------------------------------------------------
+    -- . MM format: (cint16)RAM_EQUALIZER_GAINS[S_pn/Q_fft]_[Q_fft][N_sub] = [S_pn][N_sub]
+    v_addr := g_sp * c_sdp_N_sub + g_subband;
+    -- . read
+    mmf_mm_bus_rd(c_mm_file_ram_equalizer_gains, v_addr, rd_data, tb_clk);
+    v_re := unpack_complex_re(rd_data, c_sdp_W_sub_weight);
+    v_im := unpack_complex_im(rd_data, c_sdp_W_sub_weight);
+    sp_subband_weight_re <= v_re;
+    sp_subband_weight_im <= v_im;
+    sp_subband_weight_gain <= SQRT(REAL(v_re)**2.0 + REAL(v_im)**2.0) / REAL(c_sdp_unit_sub_weight);
+    sp_subband_weight_phase <= atan2(Y => REAL(v_im), X => REAL(v_re)) * 360.0 / MATH_2_PI;
+    -- . write
+    v_weight := pack_complex(re => c_subband_weight_re, im => c_subband_weight_im, w => c_sdp_W_sub_weight);  -- c_sdp_W_sub_weight = 16 bit
+    mmf_mm_bus_wr(c_mm_file_ram_equalizer_gains, v_addr, v_weight, tb_clk);
+    -- . read back
+    mmf_mm_bus_rd(c_mm_file_ram_equalizer_gains, v_addr, rd_data, tb_clk);
+    v_re := unpack_complex_re(rd_data, c_sdp_W_sub_weight);
+    v_im := unpack_complex_im(rd_data, c_sdp_W_sub_weight);
+    sp_subband_weight_re <= v_re;
+    sp_subband_weight_im <= v_im;
+    sp_subband_weight_gain <= SQRT(REAL(v_re)**2.0 + REAL(v_im)**2.0) / REAL(c_sdp_unit_sub_weight);
+    sp_subband_weight_phase <= atan2(Y => REAL(v_im), X => REAL(v_re)) * 360.0 / MATH_2_PI;
+
+    ----------------------------------------------------------------------------
     -- Wait for enough WG data and start of sync interval
-    
-    mmf_mm_wait_until_value(c_mm_file_reg_bsn_scheduler_wg, 0,                   -- read BSN low
-                            "UNSIGNED", rd_data, ">=", c_nof_block_per_sync * 3,   -- this is the wait until condition
+    ----------------------------------------------------------------------------
+    mmf_mm_wait_until_value(c_mm_file_reg_bsn_scheduler_wg, 0,                                 -- read BSN low
+                            "UNSIGNED", rd_data, ">=", c_init_bsn + c_nof_block_per_sync * 3,  -- this is the wait until condition
                             c_sdp_T_sub, tb_clk);
 
     ---------------------------------------------------------------------------
     -- Read subband statistics
     ---------------------------------------------------------------------------
-    -- . the subband statistics are c_wpfb_sim.stat_data_sz = 2 word power values.
-    -- . there are c_sdp_N_sub = 512 subbands per signal path
-    -- . one complex WPFB can process two real inputs A, B
+    -- . the subband statistics are c_stat_data_sz = 2 word power values.
+    -- . there are c_sdp_S_pn = 12 signal inputs A, B, C, D, E, F, G, H, I, J, K, L
+    -- . there are c_sdp_N_sub = 512 subbands per signal input (SI, = signal path, SP)
+    -- . one complex WPFB can process two real inputs A, B, so there are c_sdp_P_pfb = 6 WPFB units,
+    --   but only read for the 1 WPFB unit of the selected g_sp, to save sim time
+    -- . the outputs for A, B are time multiplexed, c_sdp_Q_fft = 2, assume that they
+    --   correspond to the c_sdp_N_pol = 2 signal polarizations
     -- . the subbands are output alternately so A0 B0 A1 B1 ... A511 B511 for input A, B
     -- . the subband statistics multiple WPFB units appear in order in the ram_st_sst address map
     -- . the subband statistics are stored first lo word 0 then hi word 1
-    
-    FOR I IN 0 TO c_nof_pfb*c_nof_complex*c_sdp_N_sub*c_wpfb_sim.stat_data_sz-1 LOOP
-      v_W := I MOD c_wpfb_sim.stat_data_sz;
-      v_T := (I / c_wpfb_sim.stat_data_sz) MOD c_nof_complex;
-      v_U := I / (c_nof_complex*c_wpfb_sim.stat_data_sz*c_sdp_N_sub);
-      v_S := v_T + v_U * c_nof_complex;
-      v_B := (I / (c_nof_complex*c_wpfb_sim.stat_data_sz)) MOD c_sdp_N_sub;
-      IF v_W=0 THEN
-        -- low part
-        mmf_mm_bus_rd(c_mm_file_ram_st_sst, I, rd_data, tb_clk);
-        sp_subband_powers_arr2(v_S)(v_B)(31 DOWNTO 0) <= rd_data;
-      ELSE      
-        -- high part
-        mmf_mm_bus_rd(c_mm_file_ram_st_sst, I, rd_data, tb_clk);
-        sp_subband_powers_arr2(v_S)(v_B)(63 DOWNTO 32) <= rd_data;
-
-        -- Convert STD_LOGIC_VECTOR to REAL
-        v_sp_subband_power := REAL(TO_UINT(rd_data(29 DOWNTO 0) & 
-            sp_subband_powers_arr2(v_S)(v_B)(31 DOWNTO 30)))*2.0**30 + 
-            REAL(TO_UINT(sp_subband_powers_arr2(v_S)(v_B)(29 DOWNTO 0)));
-        -- sum
-        sp_subband_power_sum(v_S) <= sp_subband_power_sum(v_S) + v_sp_subband_power;
+    v_len := c_sdp_N_sub * c_sdp_N_pol * c_stat_data_sz;  -- 2048 = 512 * 2 * 64/32
+    v_span := true_log_pow2(v_len);                       -- = 2048
+    FOR I IN 0 TO v_len-1 LOOP
+      v_W := I MOD c_stat_data_sz;                   -- 0, 1 per statistics word, word index
+      v_P := (I / c_stat_data_sz) MOD c_sdp_N_pol;   -- 0, 1 per SP pol, polarization index
+      v_B := I / (c_sdp_N_pol * c_stat_data_sz);     -- subband index, range(N_sub = 512) per dual pol
+      v_addr := I + c_pfb_index * v_span;            -- MM address for WPFB unit of selected g_sp
+      -- Only read SST for g_subband for dual pol SP, to save sim time
+      IF g_read_all_SST = TRUE OR v_B = g_subband THEN
+        IF v_W = 0 THEN
+          -- low part
+          mmf_mm_bus_rd(c_mm_file_ram_st_sst, v_addr, rd_data, tb_clk);
+          v_data_lo := rd_data;
+        ELSE
+          -- high part
+          mmf_mm_bus_rd(c_mm_file_ram_st_sst, v_addr, rd_data, tb_clk);
+          v_data_hi := rd_data;
+          v_stat_data := v_data_hi & v_data_lo;
+
+          sp_subband_ssts_arr2(v_P)(v_B) <= v_stat_data;
+          stat_data <= v_stat_data;  -- for time series view in Wave window
+
+          -- sum of all subband powers per pol
+          sp_subband_sst_sum_arr(v_P) <= sp_subband_sst_sum_arr(v_P) + TO_UREAL(v_stat_data);
+        END IF;
       END IF;
     END LOOP;
+    proc_common_wait_some_cycles(tb_clk, 1);
 
-    -- sp_subband_power_sum is the sum of all subband powers per SP, this value will be close to sp_subband_power
-    -- because the input is a sinus, so most power will be in 1 subband. The sp_subband_power_leakage_sum shows
-    -- how much power from the input sinus at a specific subband has leaked into the 511 other subbands.
-    sp_subband_power_0 <= REAL(TO_UINT(sp_subband_powers_arr2(0)(INTEGER(ROUND(c_subband_sp_0)))(61 DOWNTO 30)))*2.0**30 + 
-        REAL(TO_UINT(sp_subband_powers_arr2(0)(INTEGER(ROUND(c_subband_sp_0)))(29 DOWNTO 0)));
-
-    sp_subband_power_sum_0 <= sp_subband_power_sum(0);
-    
+    -- Subband power of g_subband in g_sp
+    -- . For the selected g_subband in g_sp the sp_subband_sst will be close
+    --   to sp_subband_sst_sum_arr(c_pol_index), because the input is a
+    --   sinus, so most power will be in 1 subband.
+    sp_subband_sst <= TO_UREAL(sp_subband_ssts_arr2(c_pol_index)(g_subband));
     proc_common_wait_some_cycles(tb_clk, 1);
 
+    -- The sp_subband_sst_leakage shows how much power from the input sinus at a specific
+    -- subband has leaked into the N_sub-1 = 511 other subbands. The power ratio yields an
+    -- indication of the SNR, although that also depends on the SNR of the WG sinus.
+    v_power := sp_subband_sst_sum_arr(c_pol_index) - sp_subband_sst;
+    sp_subband_sst_leakage <= v_power;
+    IF sp_subband_sst > c_eps AND v_power > c_eps THEN
+      sp_subband_sst_leakage_snr_dB <= 10.0 * LOG10(sp_subband_sst / v_power);
+    ELSIF g_read_all_SST THEN
+      REPORT "Wrong, zero leakage is unexpected for SP-" & NATURAL'IMAGE(g_sp) SEVERITY ERROR;
+    END IF;
+
+    -- The sp_subband_sst_crosstalk shows how much power from one WPFB input cross talks
+    -- into the other output, due to quantization cross talk in the complex FFT. The power
+    -- ration indicates the suppression, provided that the other input was zero.
+    v_power := sp_subband_sst_sum_arr(not_int(c_pol_index));
+    sp_subband_sst_crosstalk <= v_power;
+    IF sp_subband_sst > c_eps AND v_power > c_eps THEN
+      sp_subband_sst_crosstalk_snr_dB <= 10.0 * LOG10(sp_subband_sst / v_power);
+    ELSIF g_read_all_SST THEN
+      REPORT "Wrong, zero crosstalk is unexpected for SP-" & NATURAL'IMAGE(g_sp) SEVERITY ERROR;
+    END IF;
+
+    proc_common_wait_some_cycles(tb_clk, 10);
+
+    ---------------------------------------------------------------------------
+    -- Log WG, subband statistics
+    ---------------------------------------------------------------------------
+    print_str("");
+    print_str("WG:");
+    print_str(". c_wg_ampl                            = " & int_to_str(c_wg_ampl));
+    print_str(". c_exp_sp_power                       = " & real_to_str(c_exp_sp_power, 20, 1));
+    print_str(". c_exp_sp_ast                         = " & real_to_str(c_exp_sp_ast, 20, 1));
+
+    print_str("");
+    print_str("Subband selector:");
+    print_str(". sst_offload_weighted_subbands        = " & sl_to_str(sst_offload_weighted_subbands));
+
+    print_str("");
+    print_str("Subband weight:");
+    print_str(". sp_subband_weight_gain               = " & real_to_str(sp_subband_weight_gain, 20, 6));
+    print_str(". sp_subband_weight_phase              = " & real_to_str(sp_subband_weight_phase, 20, 6));
+
+    print_str("");
+    print_str("SST results:");
+    print_str(". exp_subband_ampl                     = " & int_to_str(NATURAL(exp_subband_ampl)));
+    print_str(". exp_subband_power                    = " & real_to_str(exp_subband_power, 20, 1));
+    print_str(". exp_subband_sst                      = " & real_to_str(exp_subband_sst, 20, 1));
+    print_str("");
+    print_str(". sp_subband_sst                       = " & real_to_str(sp_subband_sst, 20, 1));
+    print_str(". sp_subband_sst / exp_subband_sst     = " & real_to_str(sp_subband_sst / exp_subband_sst, 20, 6));
+
+    IF g_read_all_SST THEN
+      -- Log WPFB details, these are allready verified in tb of wpfb_unit_dev.vhd, so here
+      -- quality indicators like leakage and crosstalk are also reported out of interest.
+      print_str("");
+      print_str("SST quality indicators");
+      print_str(". sp_subband_sst_leakage               = " & real_to_str(sp_subband_sst_leakage, 20, 0));
+      print_str(". sp_subband_sst_leakage_snr_dB        = " & real_to_str(sp_subband_sst_leakage_snr_dB, 20, 3));
+      print_str(". sp_subband_sst_crosstalk             = " & real_to_str(sp_subband_sst_crosstalk, 20, 0));
+      print_str(". sp_subband_sst_crosstalk_snr_db      = " & real_to_str(sp_subband_sst_crosstalk_snr_db, 20, 3));
+    END IF;
+
+    ---------------------------------------------------------------------------
+    -- Verify SST
     ---------------------------------------------------------------------------
-    -- Verify subband statistics
-    ---------------------------------------------------------------------------  
     -- verify expected subband power based on WG power
-    IF sp_subband_power_sum_0>0.0 THEN ASSERT sp_subband_power_0 > c_lo_factor * c_exp_subband_power_sp_0 REPORT "Wrong subband power for SP 0" SEVERITY ERROR; END IF;
-    IF sp_subband_power_sum_0>0.0 THEN ASSERT sp_subband_power_0 < c_hi_factor * c_exp_subband_power_sp_0 REPORT "Wrong subband power for SP 0" SEVERITY ERROR; END IF;
-    
-    -- view c_exp_sp_subband_power_ratio in Wave window
-    IF sp_subband_power_sum_0>0.0 THEN sp_subband_power_ratio_0 <= sp_subband_power_0/sp_subband_power_sum_0; END IF;
-    
-    -- view c_exp_sp_subband_power_sum_ratio in Wave window
-    -- The sp_subband_power_sum_ratio show similar information as sp_subband_power_leakage_sum, because when
-    -- sp_subband_power_leakage_sum is small then sp_subband_power_sum_ratio ~= sp_subband_power_ratio.
-    IF sp_subband_power_sum_0>0.0 THEN sp_subband_power_sum_ratio_0 <= sp_subband_power_sum_0/sp_subband_power_0; END IF;
-
-    -- View sp_subband_power_leakage_sum in Wave window
-    IF sp_subband_power_sum_0>0.0 THEN sp_subband_power_leakage_sum_0 <= sp_subband_power_sum_0 - sp_subband_power_0; END IF;
+    ASSERT sp_subband_sst > c_lo_factor * exp_subband_sst REPORT "Wrong subband power for SP " & NATURAL'IMAGE(g_sp) SEVERITY ERROR;
+    ASSERT sp_subband_sst < c_hi_factor * exp_subband_sst REPORT "Wrong subband power for SP " & NATURAL'IMAGE(g_sp) SEVERITY ERROR;
 
+    IF g_read_all_SST THEN
+      -- Verify expected SNR quality measures
+      ASSERT sp_subband_sst_leakage_snr_dB   > c_exp_subband_sst_leakage_snr_dB   REPORT "Wrong to much leakage for SP " & NATURAL'IMAGE(g_sp) SEVERITY ERROR;
+      ASSERT sp_subband_sst_crosstalk_snr_dB > c_exp_subband_sst_crosstalk_snr_dB REPORT "Wrong to much crosstalk for SP " & NATURAL'IMAGE(g_sp) SEVERITY ERROR;
+    END IF;
+
+    ---------------------------------------------------------------------------
+    -- End Simulation
     ---------------------------------------------------------------------------
-    -- End Simulation 
-    ---------------------------------------------------------------------------  
     sim_done <= '1';
     proc_common_wait_some_cycles(ext_clk, 100);
     proc_common_stop_simulation(TRUE, ext_clk, sim_done, tb_end);