Put data path processing and quantization model related scripts and drawings here in model/.

applications/lofar2/model/apertif_arts_firmware_model.py

+###############################################################################
+#
+# Copyright (C) 2018
+# ASTRON (Netherlands Institute for Radio Astronomy) <http://www.astron.nl/>
+# P.O.Box 2, 7990 AA Dwingeloo, The Netherlands
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program.  If not, see <http://www.gnu.org/licenses/>.
+#
+###############################################################################
+
+# Author: Eric Kooistra, 28 Oct 2018,
+
+"""Model the Apertif BF, XC / Arts TAB, IAB firmware data path signal levels
+
+The Apertif / Arts data path processing is described in:
+
+  ASTRON_SP_062_Detailed_Design_Arts_FPGA_Beamformer.pdf
+
+The purpose of the apertif_arts_firmware_model is to understand signal levels
+in the data path processing and to prove that the internal requantization in
+the firmware implementation is does not increase the system noise while
+providing the maximum dynamic range, given data widths at the external
+interfaces (i.e. W_adc = 8, W_beamlet = 6, W_vis = 32, W_volt = 8 and
+W_pow = 8).
+ 
+This apertif_arts_firmware_model.py model consists of the following sections:
+
+- User parameters
+- Firmware parameters
+- A static model to determine the expected internal signal levels as function
+  of the WG amplitude or the ADC noise sigma level.
+- A dynamic model using FFT, BF and ()**2 to model the data path signal
+  processing for a WG input signal or a system noise input signal. The FIR
+  part of the F_sub and F_chan is not modelled.
+  . WG input for test purposes to verify the linearity and dynamic range
+    as function of the WG amplitude and the number of dishes.
+  . Noise input for operational purposes to verify that the system noise 
+    maintains represented by suffcient LSbits throughout the data path
+    processing chain from ADC via beamlet and fine channel to correlator
+    visibility, Arts SC beamlet TAB, Arts SC4 power TAB and IAB.
+- Plot results of the static and the dynamic model.
+
+The model matches the drawings of the firmware parameters and signal levels in:
+
+  apertif_arts_firmware_model.vsd
+  apertif_arts_firmware_model.pdf 
+
+More information and a log an plots are also reported in:
+
+  ASTRON_RP_010888_Apertif_Arts_firmware_model.pdf
+  
+Parameters:
+* wgEnable
+  When True then a WG input is used for the dynamic model, else a noise input.
+* wgSweep
+  The WG makes a sweep through wgSweep number of subbands if wgSweep != 0. The
+  wgSweep can be < 1 and negative. If wgSweep != 0 then an image plot is made
+  of the subbands spectrum and the channels spectrum. The expected signal
+  levels that are reported by the dynamic model are only relevant when wgSweep
+  = 0.
+* quantize
+  The quantize parameter provides rounding of internal signal levels, weights
+  and gain factors. For the dynamic model clipping is not supported, because
+  the static model already shows when the signals clip. In the static model the
+  clipping can be controlled via the adjust parameter.
+
+Remark:
+- power = rms**2
+- mean powers = (rms voltages)**2
+- rms**2 = std**2 + mean**2 = std**2 when mean = 0
+- rms = std = ampl / sqrt(2) when mean is 0
+  --> pi_apr.wgScale = sqrt(2)
+  
+- power complex = power real + power imag
+                = (rms real)**2 + (rms imag)**2
+                = (ampl real/sqrt(2))**2 + (ampl imag/sqrt(2))**2
+- power real = power imag = power complex / 2, when signal is random or periodic
+- rms real = rms imag = rms complex / sqrt(2), when signal is random or periodic
+  --> pi_apr.rcScale = sqrt(2) = sqrt(nof_complex)
+- rms complex = sqrt(power complex)
+- ampl real = ampl imag = sqrt(power complex / 2) * sqrt(2)
+                        = sqrt(power complex)
+                        = rms complex
+                        = (rms real) * sqrt(2)
+                        = (rms imag) * sqrt(2)
+
+A = 1.0
+N = 100000
+realSamples = A * np.random.randn(N)
+realPowers = realSamples * realSamples
+x = rms(realSamples)**2
+y = np.mean(realPowers)
+z = np.std(realPowers)
+x # = A * 1.0
+y # = A * 1.0
+z # = A * sqrt(2)
+complexSamples =  A * np.random.randn(N) +  A * np.random.randn(N) * 1j
+complexPowers = complexSamples * np.conjugate(complexSamples)
+complexPowers = complexPowers.real
+x  = rms(complexSamples)**2
+xr = rms(complexSamples.real)**2
+xi = rms(complexSamples.imag)**2
+y = np.mean(complexPowers)
+z = np.std(complexPowers)
+x # = A * 2.0
+xr
+xi
+y # = A * 2.0
+z # = A * 2.0
+
+Use prefix rmsComplexY for rms of complex signal Y
+Use prefix rmsY for rms of real signal Y or of Y.real or of Y.imag, so
+    for complex signal Y then rmsY = rmsComplexY / pi_apr.rcScale,
+    for real signal Y then rmsY = rmsY.
+
+Usage:
+# in iPython:
+import os
+filepath='E:\\svn_radiohdl\\trunk\\applications\\apertif\\matlab'
+os.chdir(filepath)
+run apertif_arts_firmware_model.py
+
+# on command line
+python apertif_arts_firmware_model.py
+
+"""
+
+import argparse
+import math
+import numpy as np
+import common as cm
+import test_plot
+import pi_apertif_system as pi_apr
+
+################################################################################    
+# General
+################################################################################    
+
+def rms(arr):
+    """Root mean square of values in arr
+    Use this root-mean-square rms() instead of std() for CW signals that are at
+    center of bin and result in static phasor with std() = 0. For noise
+    signals with zero mean rms() is equal to std(). This rms() also works for
+    complex input thanks to using np.abs().
+    """
+    return np.sqrt(np.mean(np.abs(arr)**2.0))
+
+nofSigma = 4    # Number of sigma to fit in full scale
+
+sc1_use_qua = True       # use SC1 TAB bf quantized and clipped output
+sc1_use_qua = False      # use SC1 TAB bf raw LSbits and wrapped MSbits output (as implemented in arts_unb1_sc1)
+
+################################################################################    
+# Parse arguments to derive user settings
+################################################################################    
+
+_parser = argparse.ArgumentParser('apertif_arts_firmware_model')
+_parser.add_argument('--app', default='all', type=str, help='Application to model: dish, apertif, arts_sc1, arts_sc4 or all')
+_parser.add_argument('--model', default='all', type=str, help='Model: static, dynamic or all')
+_parser.add_argument('--sky_sigma', default=4.0, type=float, help='Signal sigma level at ADC input')
+_parser.add_argument('--wg_enable', default=False, action='store_true', help='Model coherent ADC input using WG sinus or model incoherent ADC input using Gaussian white noise')
+_parser.add_argument('--wg_ampl', default=4.0 * 2**0.5, type=float, help='WG amplitude at ADC input')
+_parser.add_argument('--wg_sweep', default=0.0, type=float, help='Number of subbands fsub to sweep with WG in dynamic model')
+_parser.add_argument('--nof_periods', default=250, type=int, help='Number of channel periods in dynamic model')
+_parser.add_argument('--nof_dishes', default=1, type=int, help='Number of dishes')
+_parser.add_argument('--quantize', default=False, action='store_true', help='Use fixed point rounding or no rounding of weights, gains and internal signals')
+_parser.add_argument('--useplot', default=False, action='store_true', help='Use plots or skip plotting')
+_parser.add_argument('--cb_weight_adjust', default=1.0, type=float, help='Factor to adjust unit CB weight')
+_parser.add_argument('--tab_weight_adjust_sc1', default=1.0, type=float, help='Factor to adjust unit TAB weight for SC1')
+_parser.add_argument('--tab_weight_adjust_sc4', default=1.0, type=float, help='Factor to adjust unit TAB weight for SC4')
+_parser.add_argument('--iab_gain_adjust', default=1.0, type=float, help='Factor to adjust unit IAB gain for SC4')
+args = _parser.parse_args()
+
+if args.app=='all':
+    apps = ['dish', 'apertif', 'arts_sc1', 'arts_sc4']
+else:
+    apps = [args.app]
+if args.model=='all':
+    models = ['static', 'dynamic']
+else:
+    models = [args.model]
+
+skySigma  = args.sky_sigma
+#skySigma  = 4
+#skySigma  = 89.8        # 89.8 * 2**0.5 = 127 for full scale WG amplitude
+#skySigma  = 26.1        # 26.1 * 2**0.5 =  37
+#skySigma  = 11.3        # 11.3 * 2**0.5 =  16 for full scale 8 bit beamlet amplitude
+
+wgAmpl    = args.wg_ampl
+wgSigma   = wgAmpl / np.sqrt(2)
+
+wgEnable  = args.wg_enable
+wgSweep   = args.wg_sweep
+#wgSweep   = 5.0          # number of fsub to sweep over nofTchanx
+#wgSweep   = -5.0         # number of fsub to sweep over nofTchanx
+#wgSweep   = 2.0          # number of fsub to sweep over nofTchanx
+#wgSweep   = 1.0/32       # number of fsub to sweep over nofTchanx
+if wgSweep != 0:
+    wgEnable = True   # force WG
+    apps = ['dish']   # only support sweep for dish
+
+quantize  = args.quantize
+
+useplot   = args.useplot
+
+cbWeightAdjust     = args.cb_weight_adjust
+tabWeightAdjustSc1 = args.tab_weight_adjust_sc1
+tabWeightAdjustSc4 = args.tab_weight_adjust_sc4
+iabGainAdjust      = args.iab_gain_adjust
+#cbWeightAdjust     = 1.0/8    # use <= 1/8 to have no beamlet (8 bit) clipping
+#cbWeightAdjust     = 1.0/32   # use <= 1/32 to have no fine channel_x (9 bit) clipping
+#tabWeightAdjustSc1 = 1.0      # = 1.0 to use unit weight
+#tabWeightAdjustSc4 = 1.0/4    # use <= 1/4 and cbWeightAdjust = 1/32 to have no TAB Stokes I (8 bit) clipping
+#iabGainAdjust      = 1.0/16   # use <= 1/16 and cbWeightAdjust = 1/32 to have no IAB Stokes I (8 bit) clipping
+
+# Number of dishes
+# - not used for Apertif correlator, because it only models auto visibilities
+# - only used for Arts array beamformers:
+#   . use  1 to verify the expected TAB, IAB output signal level
+#   . use 12 is all dishes for Arts SC1 or for IAB
+#   . use 10 is equidistant dishes for Arts SC4
+#   . use  8 is equidistant and fixed dishes for Arts SC3
+# The actual number of dishes does not alter the TAB or IAB output level, 
+# because the number of dishes is compensated via the TAB weights in case of
+# SC1 when sc1_use_qua = True and SC4 or via the IAB gain in case of the IAB.
+# When sc1_use_qua = False then the TAB output does increase with the number
+# of dishes, because then the TAB weigth is fixed at 1. Still for SC1 the
+# maximum gain with N_dish = 12 dishes is log2(sqrt(12)) = 1.8 bits, so with
+# 6 bit beamlet input and 8 bit voltage TAB data this gain still fits, and
+# therefore it is fine to use sc1_use_qua = False for SC1.
+# Typical values:
+# . nofDishes = 16    # model only
+# . nofDishes = 12    # RT2-9,a-d, all for Apertif XC and Arts SC1
+# . nofDishes = 10    # RT2-9,a,b at equidistant for Arts SC4
+# . nofDishes = 8     # RT2-9 fixed at equidistant for commensal Arts SC3
+# . nofDishes = 1     # model default, because tabWeightSc1, tabWeightSc4 and
+#                     # iabGainReg nomalize for nofDishes. Therefore use 
+#                     # nofDishes = 1 to speed up simulation, except when modeling
+#                     # Arts SC1 with fixed tabWeightSc1 = 1 (when sc1_use_qua =
+#                     # False).
+nofDishes = args.nof_dishes
+
+# Time series
+if 'dynamic' in models:
+    # use args.nof_periods for number of channel periods to simulate:
+    if wgEnable:
+        nofTchanx = args.nof_periods
+        #nofTchanx = 10          # use >~10 for WG
+        #nofTchanx = 100
+    else:
+        nofTchanx = args.nof_periods
+        #nofTchanx = 100         # number of channel periods to simulate, use >~ 100 for noise
+        #nofTchanx = 1000        # use >= 1000 for more accuracy
+        #nofTchanx = 250        # use < 1000 to speed up simulation
+    nofTStokes = nofTchanx
+
+    # also use args.nof_periods for number of Ts, Tsub, Tchanx periods to plot
+    # . to speed up plotting for Ts and Tsub
+    # . to have same density of samples per plot
+    nof_periods = args.nof_periods
+
+# Frequency selection
+if 'dynamic' in models:
+    # Select frequency for WG, SST, BST, XC
+    selSub           = 65
+    selChanx         = 33          # N_chan_x / 2 = center of subband
+    selChanxFraction = 1 * 1.0/nofTchanx   # yields 1 rotation in fine channel
+    #selChanxFraction = 0.0         # center channel, yields constant in fine channel
+    if wgSweep:
+        selChanx         = 32      # N_chan_x / 2 = center of subband
+        selChanxFraction = 0.0     # center channel, yields constant in fine channel
+
+################################################################################
+# Derived gains
+
+# Aperif BF subband CB voltage weight
+cbWeight = pi_apr.cbUnitWeight * cbWeightAdjust         # = 8192 = 2**13 * 1.0
+if quantize:
+    cbWeight = np.round(cbWeight)
+cbGain = cbWeight/pi_apr.cbMaxWeight   # normalized gain relative to maximum gain 1.0
+
+# Arts SC1 beamlet TAB voltage weight
+if sc1_use_qua:
+    # Intended implementation using quantized and clipped BF data with 
+    # TAB weights that are scaled by 1/nofDishes for coherent WG input
+    # and by 1/sqrt(nofDishes) for incoherent noise input.
+    W_tab_unit_scale_sc1 = -2                                   # = -2
+    W_tab_weight_fraction_sc1 = pi_apr.W_tab_max_weight_sc1 + W_tab_unit_scale_sc1  # = 
+    tabUnitWeightSc1 = 2.0**W_tab_weight_fraction_sc1           # = 2**13 = 8192
+    if wgEnable:
+        tabWeightSc1 = tabUnitWeightSc1 / nofDishes             # coherent TAB input, like WG
+    else:
+        tabWeightSc1 = tabUnitWeightSc1 / math.sqrt(nofDishes)  # incoherent TAB input, like sky noise
+    tabWeightSc1 *= tabWeightAdjustSc1
+else:
+    # Current implementation using LSbits and wrapped MSbits of raw BF data
+    # and unit weight = 1 independent of nofDishes. The TAB weight
+    # tabWeightSc1 = 1 is the smallest weight value, so effectively the range
+    # of W_tab_weight_sc1 is not used.
+    W_tab_unit_scale_sc1 = -pi_apr.W_tab_max_weight_sc1         # = -15
+    W_tab_weight_fraction_sc1 = pi_apr.W_tab_max_weight_sc1 + W_tab_unit_scale_sc1 # = 0 = 15 + -15
+    tabUnitWeightSc1 = 1.0                                      # = 1 = 2.0**W_tab_weight_fraction_sc1 = 2**0
+    tabWeightSc1 = tabUnitWeightSc1                             # = 1 independent of nofDishes
+if quantize:
+    tabWeightSc1 = np.round(tabWeightSc1)
+tabGainSc1 = tabWeightSc1/pi_apr.tabMaxWeightSc1  # Normalized gain relative to maximum gain 1.0
+
+# Arts SC4 channel TAB voltage weight
+if wgEnable:
+    tabWeightSc4 = pi_apr.tabUnitWeightSc4 / nofDishes             # coherent TAB input, like WG
+else:
+    tabWeightSc4 = pi_apr.tabUnitWeightSc4 / math.sqrt(nofDishes)  # incoherent TAB input, like sky noise
+tabWeightSc4 = tabWeightSc4 * tabWeightAdjustSc4
+if quantize:
+    tabWeightSc4 = np.round(tabWeightSc4)
+tabGainSc4 = tabWeightSc4/pi_apr.tabMaxWeightSc4    # Normalized gain relative to maximum gain 1.0
+
+# Arts SC4 power IAB power gain
+iabGainReg = pi_apr.iabUnitGain * iabGainAdjust / nofDishes     # IAB input powers add coherently
+if quantize:
+    iabGainReg = np.round(iabGainReg)
+iabGain = iabGainReg/pi_apr.iabMaxGain    # Normalized gain relative to maximum gain 1.0
+
+
+################################################################################
+# Apertif / Arts settings
+# . BF
+N_fft         = pi_apr.N_FFT            # = 1024
+N_sub         = pi_apr.N_sub            # = 512
+N_int_sub     = pi_apr.N_int_x          # = 800000
+
+# . XC, SC4
+N_chan_x      = pi_apr.N_chan_x         # = 64
+N_chan_bf     = pi_apr.N_chan_bf        # = 4
+N_int_chan_x  = pi_apr.N_int_chan_x     # = 12500 = N_int_x / N_chan_x
+
+if 'static' in models:
+    ################################################################################    
+    # Static model of signal level as function of:
+    # . WG input amplitude level
+    # . Noise input sigma level
+    ################################################################################    
+    
+    # Input axis
+    inputResolution = 0.05
+    
+    # Dish
+    A_wg = np.arange(inputResolution, 2**pi_apr.W_adc_max, inputResolution)    # 0:127
+    S_wg = A_wg / pi_apr.wgScale
+    S_wg_complex = S_wg
+    S_wg_complex_phase = S_wg_complex / pi_apr.rcScale
+    A_wg_complex_phase = S_wg_complex
+    S_noise = np.arange(inputResolution, 2**pi_apr.W_adc_max / nofSigma, inputResolution)  # 0:31
+    S_noise_complex = S_noise
+    S_noise_complex_phase = S_noise / pi_apr.rcScale
+    
+    SST_wg_inputs = S_wg_complex / pi_apr.wgSstScale
+    SST_wg_values = SST_wg_inputs**2 * N_int_sub
+    SST_wg_dBs = 10*np.log10(SST_wg_values)
+    SST_noise_inputs = S_noise_complex / pi_apr.noiseSstScale
+    SST_noise_values = SST_noise_inputs**2 * N_int_sub
+    SST_noise_dBs = 10*np.log10(SST_noise_values)
+    
+    BST_wg_inputs = S_wg_complex * cbWeightAdjust / pi_apr.wgBstScale
+    BST_wg_values = BST_wg_inputs**2 * N_int_sub
+    BST_wg_dBs = 10*np.log10(BST_wg_values)
+    BST_noise_inputs = S_noise_complex * cbWeightAdjust / pi_apr.noiseBstScale
+    BST_noise_values = BST_noise_inputs**2 * N_int_sub
+    BST_noise_dBs = 10*np.log10(BST_noise_values)
+    
+    clip = 2.0**pi_apr.W_beamlet_max - 1
+    beamlet_wg_complex_phase_values = S_wg_complex_phase * cbWeightAdjust / pi_apr.wgBeamletScale
+    beamlet_wg_complex_phase_values = np.clip(beamlet_wg_complex_phase_values, -clip, clip)
+    beamlet_noise_complex_phase_values = S_noise_complex_phase * cbWeightAdjust / pi_apr.noiseBeamletScale
+    beamlet_noise_complex_phase_values = np.clip(beamlet_noise_complex_phase_values, -clip, clip)
+    
+    # Central
+    clip = 2.0**pi_apr.W_channel_x_max - 1
+    channel_x_wg_complex_phase_values = S_wg_complex_phase * cbWeightAdjust / pi_apr.wgChannelXcorScale
+    channel_x_wg_complex_phase_values = np.clip(channel_x_wg_complex_phase_values, -clip, clip)
+    channel_x_wg_complex_values = channel_x_wg_complex_phase_values * pi_apr.rcScale
+    channel_x_noise_complex_phase_values = S_noise_complex_phase * cbWeightAdjust / pi_apr.noiseChannelXcorScale
+    channel_x_noise_complex_phase_values = np.clip(channel_x_noise_complex_phase_values, -clip, clip)
+    channel_x_noise_complex_values = channel_x_noise_complex_phase_values * pi_apr.rcScale
+    
+    if 'apertif' in apps:
+        # Correlator
+        XC_wg_auto_visibilities = channel_x_wg_complex_values**2 * N_int_chan_x
+        XC_wg_auto_visibilities_dB = 10*np.log10(XC_wg_auto_visibilities)
+        XC_noise_auto_visibilities = channel_x_noise_complex_values**2 * N_int_chan_x
+        XC_noise_auto_visibilities_dB = 10*np.log10(XC_noise_auto_visibilities)
+    
+    if 'arts_sc1' in apps:
+        # Arts SC1 beamlet voltage TAB
+        if sc1_use_qua:
+            beamlet_tab_wg_complex_phase_values_sc1 = beamlet_wg_complex_phase_values * tabWeightAdjustSc1
+            beamlet_tab_noise_complex_phase_values_sc1 = beamlet_noise_complex_phase_values * tabWeightAdjustSc1
+        else:
+            beamlet_tab_wg_complex_phase_values_sc1 = beamlet_wg_complex_phase_values * nofDishes
+            beamlet_tab_noise_complex_phase_values_sc1 = beamlet_noise_complex_phase_values * nofDishes
+        clip = 2.0**pi_apr.W_volt_max - 1
+        beamlet_tab_wg_complex_phase_values_sc1 = np.clip(beamlet_tab_wg_complex_phase_values_sc1, -clip, clip)
+        beamlet_tab_noise_complex_phase_values_sc1 = np.clip(beamlet_tab_noise_complex_phase_values_sc1, -clip, clip)
+    
+    if 'arts_sc4' in apps:
+        # Arts SC4 channel power TAB
+        clip = 2.0**pi_apr.W_volt - 1
+        power_tab_xx_wg_complex_values_sc4 = (channel_x_wg_complex_values * tabWeightAdjustSc4)**2    # Stokes I = XX* + 0
+        power_tab_xx_wg_complex_values_sc4 = np.clip(power_tab_xx_wg_complex_values_sc4, 0, clip)  # Stokes I is unsigned >= 0
+        power_tab_xx_noise_values_sc4 = (channel_x_noise_complex_values * tabWeightAdjustSc4)**2    # Stokes I = XX* + 0
+        power_tab_xx_noise_values_sc4 = np.clip(power_tab_xx_noise_values_sc4, 0, clip)  # Stokes I is unsigned >= 0
+        
+        # Arts SC4 channel IAB
+        clip = 2.0**pi_apr.W_pow - 1
+        iab_xx_wg_complex_values_sc4 = channel_x_wg_complex_values**2 * iabGainAdjust   # Stokes I = XX* + 0
+        iab_xx_wg_complex_values_sc4 = np.clip(iab_xx_wg_complex_values_sc4, 0, clip) # Stokes I is unsigned >= 0
+        iab_xx_noise_values_sc4 = channel_x_noise_complex_values**2 * iabGainAdjust   # Stokes I = XX* + 0
+        iab_xx_noise_values_sc4 = np.clip(iab_xx_noise_values_sc4, 0, clip) # Stokes I is unsigned >= 0
+
+
+if 'dynamic' in models:
+    ################################################################################    
+    # Dynamic model parameters
+    ################################################################################    
+    
+    ################################################################################
+    # Apertif / Arts settings
+    # . BF
+    fs            = 800e6
+    fsub          = fs / N_fft
+    Ts            = 1/fs
+    Tsub          = 1/fsub
+
+    # . XC, SC4
+    fchan_x       = fsub / N_chan_x
+    Tchanx        = 1/fchan_x
+    
+    # . SC4
+    Tstokesx      = Tchanx
+    M_int_ax      = pi_apr.M_int_ax         # = 16
+    W_int_ax      = pi_apr.W_int_ax         # = 4 = log2(16)
+    #M_int_ax      = 1
+    #W_int_ax      = cm.ceil_log2(M_int_ax)
+    
+    ################################################################################
+    # Dish CB
+
+    # Number of bits for internal product
+    W_cb_product = pi_apr.W_cb_in + pi_apr.W_cb_weight - 1  # = 31 = 16 + 16 -1, do -1 to skip double sign bit of product
+    
+    # Number of bits for fraction < 1
+    W_beamlet_fraction = pi_apr.W_beamlet_fraction    # number of bits for fraction < 1
+    
+    if 'arts_sc1' in apps:
+        ################################################################################    
+        # Arts SC1 beamlet TAB
+        
+        # Number of bits for internal product
+        W_tab_product_sc1 = pi_apr.W_beamlet + pi_apr.W_tab_weight_sc1 - 1  # = 23 = 8 + 16 -1, do -1 to skip double sign bit of product
+        
+        # Number of bits for fraction < 1
+        W_tab_fraction_sc1 = W_beamlet_fraction - W_tab_unit_scale_sc1  # = 8  = 6 - -2  when sc1_use_qua = True
+                                                                        # = 21 = 6 - -15 when sc1_use_qua = False
+    
+    if 'arts_sc4' in apps:
+        ################################################################################    
+        # Arts SC4 channel TAB
+        
+        # Fine channel TAB        
+        # Number of bits for internal product
+        W_tab_voltage_product_sc4 = pi_apr.W_channel_x + pi_apr.W_tab_weight_sc4 - 1  # = 17 = 9 + 9 -1, do -1 to skip double sign bit of product
+        W_tab_power_product_sc4 = 2*W_tab_voltage_product_sc4 - 1  # = 33 = 2*17 -1, do -1 to skip double sign bit of product
+        
+        # For SC4 it is possible to use W_channel_x_fraction > 9, because the F_chan is internal
+        # For SC3 it W_channel_x_fraction = 9, because the F_chan is in Apertif XC
+        # Number of bits for fraction < 1
+        W_tab_voltage_fraction_sc4 = pi_apr.W_channel_x_fraction - pi_apr.W_tab_unit_scale_sc4  # = 10 = 9 - -1
+        
+        # Arts channel TAB
+        # Number of bits for fraction < 1
+        W_tab_power_fraction_sc4 = 2*W_tab_voltage_fraction_sc4             # = 20 = 2*10
+        W_tab_power_int_fraction_sc4 = W_tab_power_fraction_sc4 - W_int_ax  # = 16 = 20 - 4
+                                                                                                                 
+        ################################################################################    
+        # Arts SC4 channel IAB
+                
+        # Number of bits for internal product
+        W_iab_power_product = 2*pi_apr.W_channel_x_fraction - 1  # = 17 = 2*9 -1, do -1 to skip double sign bit of product
+        W_iab_gain_product = W_iab_power_product + pi_apr.W_iab_gain - 1  # = 32 = 17 + 16 - 1, do -1 to skip double sign bit of product
+        
+        # Number of bits for fraction < 1
+        W_iab_power_fraction = 2*pi_apr.W_channel_x_fraction              # = 18 = 2*9
+        W_iab_gain_fraction = W_iab_power_fraction - W_int_ax - pi_apr.W_iab_unit_scale    # = 16 = 18 - 4 - -2
+
+    ################################################################################    
+    # Dynamic model time series
+    ################################################################################    
+    
+    # Time, frequency axis
+    nofTsub = nofTchanx * N_chan_x
+    nofTs = nofTsub * N_fft
+    t = np.arange(nofTs) * Ts
+    
+    ################################################################################
+    # Dish Compound Beamformer
+    
+    # ADC input
+    if wgEnable:
+        wgFreq = (selSub + (selChanx + selChanxFraction - N_chan_x/2.0) / N_chan_x) * fsub
+        chirpFreq = 0
+        if wgSweep:
+            chirpFreq = wgSweep * fsub * np.arange(nofTs)/nofTs/2.0
+        spSamples = wgAmpl * np.sin(2*np.pi*(wgFreq + chirpFreq) * t)
+    else:
+        spSamples = np.random.randn(nofTs)
+        spSamples *= skySigma / rms(spSamples)
+    
+    if quantize:
+        # The unit level is at 1 lsb of the ADC input, so the W_adc = 8 bit input
+        # has no fraction and can thus directly be rounded to model quantization.
+        spSamples = np.round(spSamples)
+    rmsSp = rms(spSamples)  # real signal
+    amplSp = rmsSp * pi_apr.wgScale
+    
+    # Fsub subbands full bandwidth
+    spBlocks = spSamples.reshape((nofTsub, N_fft))  # [nofTsub][N_fft]
+    fftSamples = np.fft.fft(spBlocks) / N_fft       # [nofTsub][N_fft]
+    fSubSamples = fftSamples[:, 0:N_sub]            # [nofTsub][N_sub], real input, so only keep positive frequencies
+    
+    # SST for selSub
+    sstSamples = fSubSamples[:, selSub].copy()
+    if quantize:
+        sstSamples *= 2.0**pi_apr.W_sst_fraction  # fixed point fraction
+        sstSamples = np.round(sstSamples)
+        sstSamples /= 2.0**pi_apr.W_sst_fraction  # back to normalized fraction
+    rmsComplexSST = rms(sstSamples)       # complex signal
+    rmsSST = rmsComplexSST / pi_apr.rcScale
+    amplSST = rmsSST * pi_apr.wgScale
+    SSTpower = N_int_sub * rmsComplexSST**2
+    
+    # BF samples full bandwidth
+    bfSamples = fSubSamples * cbGain
+    
+    # BST for selSub
+    bstSamples = bfSamples[:, selSub].copy()
+    if quantize:
+        bstSamples *= 2.0**pi_apr.W_bst_fraction  # fixed point fraction
+        bstSamples = np.round(bstSamples)
+        bstSamples /= 2.0**pi_apr.W_bst_fraction  # back to normalized fraction
+    rmsComplexBST = rms(bstSamples)    # complex signal
+    rmsBST = rmsComplexBST / pi_apr.rcScale
+    amplBST = rmsBST * pi_apr.wgScale
+    BSTpower = N_int_sub * rmsComplexBST**2
+    
+    # CB beamlet for selSub
+    beamletSamples = bfSamples[:, selSub].copy()
+    if quantize:
+        beamletSamples *= 2.0**pi_apr.W_beamlet_fraction  # fixed point fraction
+        beamletSamples = np.round(beamletSamples)
+        beamletSamples /= 2.0**pi_apr.W_beamlet_fraction  # back to normalized fraction
+    rmsComplexBeamlet = rms(beamletSamples)    # complex signal
+    rmsBeamlet = rmsComplexBeamlet / pi_apr.rcScale
+    amplBeamlet = rmsBeamlet * pi_apr.wgScale
+    
+    ################################################################################
+    # Beamlet channelizer
+    # . Fchan output
+    transposeSamples = bfSamples.transpose()                                # [nofTsub][N_sub] --> [N_sub][nofTsub]
+    transposeBlocks = transposeSamples.reshape(N_sub, nofTchanx, N_chan_x)  # [N_sub][nofTchanx][N_chan_x]
+    fSubChanSamples = np.fft.fft(transposeBlocks) / N_chan_x                # [N_sub][nofTchanx][N_chan_x], complex input and output
+    fSubChanSamples = np.fft.fftshift(fSubChanSamples, 2)
+    fSubChanSamples = fSubChanSamples.transpose(1, 0, 2)                    # [nofTchanx][N_sub][N_chan_x]
+    
+    # . CB fine channel for selSub, selChanx
+    fineChannelSamples = fSubChanSamples[:, selSub, selChanx].copy()
+    if quantize:
+        fineChannelSamples *= 2.0**pi_apr.W_channel_x_fraction  # fixed point fraction
+        fineChannelSamples = np.round(fineChannelSamples)
+        fineChannelSamples /= 2.0**pi_apr.W_channel_x_fraction  # back to normalized fraction
+    rmsComplexFineChannel = rms(fineChannelSamples)    # complex signal
+    rmsFineChannel = rmsComplexFineChannel / pi_apr.rcScale
+    amplFineChannel = rmsFineChannel * pi_apr.wgScale
+    
+    if 'apertif' in apps:
+        ################################################################################
+        # Aperif XC channel auto-visibility
+        XCpower = N_int_chan_x * rmsComplexFineChannel**2
+    
+    if 'arts_sc1' in apps:
+        ################################################################################
+        # Arts SC1 beamlet voltage TAB
+        if wgEnable:
+            # Coherent sum using beamlets from original ADC input
+            beamletTabSamples = nofDishes * beamletSamples.copy()
+        else:
+            # Incoherent sum using copies with the same rms beamlet as for the original ADC input
+            beamletTabSamples = np.zeros(nofTsub) + np.zeros(nofTsub) * 1j
+            for tel in range(nofDishes):
+                complexSamples = np.random.randn(nofTsub) + np.random.randn(nofTsub) * 1j
+                complexSamples *= rmsComplexBeamlet / rms(complexSamples)
+                beamletTabSamples += complexSamples
+        beamletTabSamples *= tabGainSc1
+        if quantize:
+            beamletTabSamples *= 2.0**W_tab_fraction_sc1  # fixed point fraction
+            beamletTabSamples = np.round(beamletTabSamples)
+            beamletTabSamples /= 2.0**W_tab_fraction_sc1  # back to normalized fraction
+        rmsComplexBeamletTab = rms(beamletTabSamples)    # complex signal
+        rmsBeamletTab = rmsComplexBeamletTab / pi_apr.rcScale
+        amplBeamletTab = rmsBeamletTab * pi_apr.wgScale
+    
+    if 'arts_sc4' in apps:
+        ################################################################################
+        # Arts SC4 Power TAB for single polarization Stokes I = XX* + 0
+        # . Fine channel voltage TAB
+        if wgEnable:
+            # Coherent sum using channel from original ADC input
+            fineChannelTabSamples = nofDishes * fineChannelSamples.copy()
+        else:
+            # Incoherent sum using copies with the same rms channel as for the original ADC input
+            fineChannelTabSamples = np.zeros(nofTchanx) + np.zeros(nofTchanx) * 1j
+            for tel in range(nofDishes):
+                complexSamples = np.random.randn(nofTchanx) + np.random.randn(nofTchanx) * 1j
+                complexSamples *= rmsComplexFineChannel / rms(complexSamples)
+                fineChannelTabSamples += complexSamples    # voltage sum
+        fineChannelTabSamples *= tabGainSc4
+        if quantize:
+            fineChannelTabSamples *= 2.0**W_tab_voltage_fraction_sc4  # fixed point fraction
+            fineChannelTabSamples = np.round(fineChannelTabSamples)
+            fineChannelTabSamples /= 2.0**W_tab_voltage_fraction_sc4  # back to normalized fraction
+        rmsComplexFineChannelTab = rms(fineChannelTabSamples)    # complex signal
+        rmsFineChannelTab = rmsComplexFineChannelTab / pi_apr.rcScale
+        amplFineChannelTab = rmsFineChannelTab * pi_apr.wgScale
+        
+        # . Fine channel power TAB
+        fineChannelTabPowers = fineChannelTabSamples * np.conjugate(fineChannelTabSamples)  # powers
+        fineChannelTabPowers = fineChannelTabPowers.real           # keep real, imag = 0
+        meanFineChannelTabPower = np.mean(fineChannelTabPowers)    # = rmsComplexFineChannelTab**2, use mean for power data and rms for voltage data
+        stdFineChannelTabPower = np.std(fineChannelTabPowers)
+        
+        # . Integrated channels power TAB
+        if wgEnable:
+            # The WG sinus only occurs in one fine channel if it is close to the center fine channel frequency, so
+            # the other fine channels in N_chan_bf fine channels carry almost zero and do not contribute much
+            integratedTabPowers = fineChannelTabPowers.copy()
+        else:
+            # Integrate using copies with the same rmsComplexFineChannelTab as for the orginal ADC input
+            integratedTabPowers = np.zeros(nofTchanx)
+            for ci in range(M_int_ax):
+                complexSamples = np.random.randn(nofTchanx) + np.random.randn(nofTchanx) * 1j
+                complexSamples *= rmsComplexFineChannelTab / rms(complexSamples)
+                powerSamples = complexSamples * np.conjugate(complexSamples)
+                powerSamples = powerSamples.real           # keep real, imag = 0
+                integratedTabPowers += powerSamples
+        if quantize:
+            integratedTabPowers *= 2.0**W_tab_power_int_fraction_sc4  # fixed point fraction
+            integratedTabPowers = np.round(integratedTabPowers)
+            integratedTabPowers /= 2.0**W_tab_power_int_fraction_sc4  # back to normalized fraction
+        meanIntegratedTabPower = np.mean(integratedTabPowers)    # use mean for power data and rms for voltage data
+        stdIntegratedTabPower = np.std(integratedTabPowers)
+        
+        ################################################################################
+        # Arts SC4 IAB for single polarization Stokes I = XX* + 0
+        
+        # . Fine channel powers
+        fineChannelPowers = fineChannelSamples * np.conjugate(fineChannelSamples)
+        fineChannelPowers = fineChannelPowers.real           # keep real, imag = 0
+        meanFineChannelPower = np.mean(fineChannelPowers)    # = rmsComplexFineChannel**2, use mean for power data and rms for voltage data
+        
+        # . Integrated channel powers
+        if wgEnable:
+            # The WG sinus only occurs in one fine channel if it is close to the center fine channel frequency, so
+            # the other fine channels in N_chan_bf fine channels carry almost zero and do not contribute much
+            integratedChannelPowers = fineChannelPowers
+        else:
+            # Integrate using copies with the same rmsComplexFineChannel as for the orginal ADC input
+            integratedChannelPowers = np.zeros(nofTchanx)
+            for ci in range(M_int_ax):
+                complexSamples = np.random.randn(nofTchanx) + np.random.randn(nofTchanx) * 1j
+                complexSamples *= rmsComplexFineChannel / rms(complexSamples)
+                powerSamples = complexSamples * np.conjugate(complexSamples)
+                powerSamples = powerSamples.real           # keep real, imag = 0
+                integratedChannelPowers += powerSamples
+        meanIntegratedChannelPower = np.mean(integratedChannelPowers)    # use mean for power data and rms for voltage data
+        
+        # . Sum of powers is coherent sum
+        if wgEnable:
+            # Sum of powers using integrated channel from original ADC input
+            integratedIabPowers = nofDishes * integratedChannelPowers
+        else:
+            # Sum of powers using copies with the same rmsComplexFineChannel as for the original ADC input
+            integratedIabPowers = np.zeros(nofTchanx)
+            for tel in range(nofDishes):
+                telPowerSamples = np.zeros(nofTchanx)
+                for ci in range(M_int_ax):
+                    complexSamples = np.random.randn(nofTchanx) + np.random.randn(nofTchanx) * 1j
+                    complexSamples *= rmsComplexFineChannel / rms(complexSamples)
+                    powerSamples = complexSamples * np.conjugate(complexSamples)
+                    powerSamples = powerSamples.real      # keep real, imag = 0
+                    telPowerSamples += powerSamples       # powers sum integration
+                integratedIabPowers += telPowerSamples    # powers sum IAB
+        meanIntegratedIabPower = np.mean(integratedIabPowers)    # use mean for power data and rms for voltage data
+        
+        # . Output gain
+        outputIabPowers = integratedIabPowers * iabGain
+        if quantize:
+            outputIabPowers *= 2.0**W_iab_gain_fraction    # fixed point fraction
+            outputIabPowers = np.round(outputIabPowers)
+            outputIabPowers /= 2.0**W_iab_gain_fraction    # back to normalized fraction
+        meanOutputIabPower = np.mean(outputIabPowers)    # use mean for power data and rms for voltage data
+        stdOutputIabPower = np.std(outputIabPowers)
+    
+    ################################################################################
+    # Dynamic model: Normalized floating point fraction to fixed point fraction
+    ################################################################################    
+    
+    # SST
+    rmsComplexSST *= 2.0**pi_apr.W_sst_fraction
+    rmsSST *= 2.0**pi_apr.W_sst_fraction
+    amplSST *= 2.0**pi_apr.W_sst_fraction
+    SSTpower *= 2.0**(2*pi_apr.W_sst_fraction)
+    
+    # BST
+    rmsComplexBST *= 2.0**pi_apr.W_bst_fraction
+    rmsBST *= 2.0**pi_apr.W_bst_fraction
+    amplBST *= 2.0**pi_apr.W_bst_fraction
+    BSTpower *= 2.0**(2*pi_apr.W_bst_fraction)
+    
+    # CB beamlet for selSub
+    beamletSamples *= 2.0**pi_apr.W_beamlet_fraction
+    rmsComplexBeamlet *= 2.0**pi_apr.W_beamlet_fraction
+    rmsBeamlet *= 2.0**pi_apr.W_beamlet_fraction
+    amplBeamlet *= 2.0**pi_apr.W_beamlet_fraction
+    
+    # CB fine channel for selSub, selChanx
+    fineChannelSamples *= 2.0**pi_apr.W_channel_x_fraction
+    rmsComplexFineChannel *= 2.0**pi_apr.W_channel_x_fraction
+    rmsFineChannel *= 2.0**pi_apr.W_channel_x_fraction
+    amplFineChannel *= 2.0**pi_apr.W_channel_x_fraction
+    
+    if 'apertif' in apps:
+        # Aperif XC channel auto-visibility
+        XCpower *= 2.0**(2*pi_apr.W_channel_x_fraction)
+    
+    if 'arts_sc1' in apps:
+        # Arts SC1 beamlet voltage TAB
+        beamletTabSamples *= 2.0**W_tab_fraction_sc1
+        rmsComplexBeamletTab *= 2.0**W_tab_fraction_sc1
+        rmsBeamletTab *= 2.0**W_tab_fraction_sc1
+        amplBeamletTab *= 2.0**W_tab_fraction_sc1
+    
+    if 'arts_sc4' in apps:
+        # Arts SC4 Power TAB for single polarization Stokes I = XX* + 0
+        # . Fine channel voltage TAB
+        fineChannelTabSamples *= 2.0**W_tab_voltage_fraction_sc4
+        rmsComplexFineChannelTab *= 2.0**W_tab_voltage_fraction_sc4
+        rmsFineChannelTab *= 2.0**W_tab_voltage_fraction_sc4
+        amplFineChannelTab *= 2.0**W_tab_voltage_fraction_sc4
+        
+        # . Fine channel power TAB
+        fineChannelTabPowers *= 2.0**W_tab_power_fraction_sc4
+        meanFineChannelTabPower *= 2.0**W_tab_power_fraction_sc4
+        stdFineChannelTabPower *= 2.0**W_tab_power_fraction_sc4
+        
+        # . Integrated channels power TAB
+        integratedTabPowers *= 2.0**W_tab_power_int_fraction_sc4
+        meanIntegratedTabPower *= 2.0**W_tab_power_int_fraction_sc4
+        stdIntegratedTabPower *= 2.0**W_tab_power_int_fraction_sc4
+        
+        # Arts SC4 IAB for single polarization Stokes I = XX* + 0
+        # . Fine channel powers
+        fineChannelPowers *= 2.0**W_iab_power_fraction
+        meanFineChannelPower *= 2.0**W_iab_power_fraction
+        
+        # . Integrated channel powers
+        integratedChannelPowers *= 2.0**W_iab_power_fraction
+        meanIntegratedChannelPower *= 2.0**W_iab_power_fraction
+        
+        # . Sum of powers is coherent sum
+        integratedIabPowers *= 2.0**W_iab_power_fraction
+        meanIntegratedIabPower *= 2.0**W_iab_power_fraction
+        
+        # . Output gain
+        outputIabPowers *= 2.0**W_iab_gain_fraction
+        meanOutputIabPower *= 2.0**W_iab_gain_fraction
+        stdOutputIabPower *= 2.0**W_iab_gain_fraction
+
+################################################################################
+# Logging
+################################################################################    
+print '--------------------------------------------------------------------'
+print '-- User settings'
+print '--------------------------------------------------------------------'
+if wgEnable:
+    print 'WG sinus input:'
+    print '. wgSigma              = %f' % wgSigma
+    print '. wgAmpl               = %f [fsub]' % wgAmpl
+    print '. wgSweep              = %f [fsub]' % wgSweep
+else:
+    print 'ADC noise input:'
+    print '. skySigma             = %f' % skySigma
+print '. quantize             = %s' % quantize
+print '. nofDishes            = %d   (= %.1f bit)' % (nofDishes, math.log(nofDishes, 2))
+print ''
+
+if 'dynamic' in models:
+    print 'Dynamic model:'
+    print '. selSub               = %d' % selSub
+    print '. selChanx             = %d' % selChanx
+    print '. selChanxFraction     = %d' % selChanxFraction
+    print '. nofTchanx            = %d' % nofTchanx
+    print '. nofTsub              = %d' % nofTsub
+    print '. nofTs                = %d' % nofTs
+    print ''
+
+print '--------------------------------------------------------------------'
+print '-- Design settings'
+print '--------------------------------------------------------------------'
+print 'Apertif BF:'
+print '. N_fft                = %d' % N_fft
+print '. N_sub                = %d' % N_sub
+print '. N_int_sub            = %d' % N_int_sub
+print ''
+if 'apertif' in apps:
+    print 'Apertif XC:'
+    print '. N_chan_x             = %d' % N_chan_x
+    print '. N_chan_bf            = %d' % N_chan_bf
+    print '. N_int_chan_x         = %d' % N_int_chan_x
+    print '. cbWeightAdjust       = %f' % cbWeightAdjust
+    print '. cbMaxWeight          = %d' % pi_apr.cbMaxWeight
+    print '. cbUnitWeight         = %d' % pi_apr.cbUnitWeight
+    print '. cbWeight             = %d' % cbWeight
+    print '. cbGain               = %f' % cbGain
+    print ''
+if 'arts_sc1' in apps:
+    print 'Arts SC1:'
+    print '. sc1_use_qua          = %s' % sc1_use_qua
+    print '. tabWeightAdjustSc1   = %f' % tabWeightAdjustSc1
+    print '. tabUnitWeightSc1     = %d' % tabUnitWeightSc1
+    print '. tabWeightSc1         = %d' % tabWeightSc1
+    print '. tabGainSc1           = %.8f' % tabGainSc1
+    print ''
+if 'arts_sc4' in apps:
+    print 'Arts SC4:'
+    print '. tabWeightAdjustSc4   = %f' % tabWeightAdjustSc4
+    print '. iabGainAdjust        = %f' % iabGainAdjust
+    print '. iabMaxGain           = %d' % pi_apr.iabMaxGain
+    print '. iabUnitGain          = %d' % pi_apr.iabUnitGain
+    print '. iabGainReg           = %d' % iabGainReg
+    print '. iabGain              = %.8f' % iabGain
+    print ''
+
+if 'dynamic' in models:
+    print '--------------------------------------------------------------------'
+    print '-- Dynamic model settings'
+    print '--------------------------------------------------------------------'
+    print 'Design:'
+    print '. fs                           = %f [Hz]' % fs
+    print '. Ts                           = %.2f [ns]' % (Ts * 1e9)
+    print '. fsub                         = %f [Hz]' % fsub
+    print '. Tsub                         = %.2f [us]' % (Tsub * 1e6)
+    print '. W_fsub_gain                  = %d [bit]' % pi_apr.W_fsub_gain
+    if 'apertif' in apps:
+        print '. fchan_x                      = %f [Hz]' % fchan_x
+        print '. Tchanx                       = %.2f [us]' % (Tchanx * 1e6)
+        print '. W_fchan_x_gain               = %d [bit]' % pi_apr.W_fchan_x_gain
+    if 'arts_sc4' in apps:
+        print '. Tstokesx                     = %.2f [us]' % (Tstokesx * 1e6)
+    print ''
+
+    print 'Dish:'
+    print '. W_adc                        = %d [bit]' % pi_apr.W_adc
+    print '. W_adc_max                    = %d [bit], adc max = %d' % (pi_apr.W_adc_max, 2**pi_apr.W_adc_max-1)
+    print '. W_sst_in                     = %d [bit]' % pi_apr.W_sst_in
+    print '. W_sst_fraction               = %d [bit]' % pi_apr.W_sst_fraction
+    print '. W_cb_weight                  = %d [bit]' % pi_apr.W_cb_weight
+    print '. W_cb_unit_scale              = %d [bit]' % pi_apr.W_cb_unit_scale
+    print '. W_cb_weight_fraction         = %d [bit]' % pi_apr.W_cb_weight_fraction
+    print '. W_cb_in                      = %d [bit]' % pi_apr.W_cb_in
+    print '. W_cb_product                 = %d [bit]' % W_cb_product
+    print '. W_bst_in                     = %d [bit]' % pi_apr.W_bst_in
+    print '. W_bst_fraction               = %d [bit]' % pi_apr.W_bst_fraction
+    print '. W_beamlet_sc1                = %d [bit]' % pi_apr.W_beamlet_sc1
+    print '. W_beamlet                    = %d [bit]' % pi_apr.W_beamlet
+    print '. W_beamlet_fraction           = %d [bit]' % pi_apr.W_beamlet_fraction
+    print ''
+
+    if 'apertif' in apps:
+        print 'Apertif XC:'
+        print '. W_channel_x                  = %d [bit]' % pi_apr.W_channel_x
+        print '. W_channel_x_fraction         = %d [bit]' % pi_apr.W_channel_x_fraction
+        print ''
+        
+    if 'arts_sc1' in apps:
+        print 'Arts SC1 beamlet TAB:'
+        print '. W_tab_weight_sc1             = %d [bit]' % pi_apr.W_tab_weight_sc1
+        print '. W_tab_unit_scale_sc1         = %d [bit]' % W_tab_unit_scale_sc1
+        print '. W_tab_weight_fraction_sc1    = %d [bit]' % W_tab_weight_fraction_sc1
+        print '. W_beamlet_sc1                = %d [bit]' % pi_apr.W_beamlet_sc1
+        print '. W_tab_product_sc1            = %d [bit]' % W_tab_product_sc1
+        print '. W_tab_fraction_sc1           = %d [bit]' % W_tab_fraction_sc1
+        print ''
+    
+    if 'arts_sc4' in apps:
+        print 'Arts SC4 TAB:'
+        print '. W_channel_x                  = %d [bit]' % pi_apr.W_channel_x
+        print '. W_channel_x_fraction         = %d [bit]' % pi_apr.W_channel_x_fraction
+        print '. W_tab_weight_sc4             = %d [bit]' % pi_apr.W_tab_weight_sc4
+        print '. W_tab_unit_scale_sc4         = %d [bit]' % pi_apr.W_tab_unit_scale_sc4
+        print '. W_tab_weight_fraction_sc4    = %d [bit]' % pi_apr.W_tab_weight_fraction_sc4
+        print '. tabUnitWeightSc4             = %d' % pi_apr.tabUnitWeightSc4
+        print '. tabWeightSc4                 = %d' % tabWeightSc4
+        print '. tabGainSc4                   = %.8f' % tabGainSc4
+        print '. W_beamlet                    = %d [bit]' % pi_apr.W_beamlet
+        print '. W_tab_voltage_product_sc4    = %d [bit]' % W_tab_voltage_product_sc4
+        print '. W_tab_voltage_fraction_sc4   = %d [bit]' % W_tab_voltage_fraction_sc4
+        print '. W_tab_power_product_sc4      = %d [bit]' % W_tab_power_product_sc4
+        print '. W_tab_power_fraction_sc4     = %d [bit]' % W_tab_power_fraction_sc4
+        print '. W_tab_power_int_fraction_sc4 = %d [bit]' % W_tab_power_int_fraction_sc4
+        print '. W_int_ax                     = %d [bit]' % W_int_ax
+        print '. M_int_ax                     = %d' % M_int_ax
+        print ''
+        
+        print 'Arts SC4 IAB:'
+        print '. W_channel_x                  = %d [bit]' % pi_apr.W_channel_x
+        print '. W_channel_x_fraction         = %d [bit]' % pi_apr.W_channel_x_fraction
+        print '. W_int_ax                     = %d [bit]' % W_int_ax
+        print '. M_int_ax                     = %d' % M_int_ax
+        print '. W_iab_gain                   = %d [bit]' % pi_apr.W_iab_gain
+        print '. W_iab_unit_scale             = %d [bit]' % pi_apr.W_iab_unit_scale
+        print '. W_iab_gain_fraction          = %d [bit]' % pi_apr.W_iab_gain_fraction
+        print '. W_iab_power_product          = %d [bit]' % W_iab_power_product
+        print '. W_iab_power_fraction         = %d [bit]' % W_iab_power_fraction
+        print '. W_iab_gain_product           = %d [bit]' % W_iab_gain_product
+        print '. W_iab_gain_fraction          = %d [bit]' % W_iab_gain_fraction
+        print ''
+
+if 'dynamic' in models and wgSweep == 0:
+    print '--------------------------------------------------------------------'
+    print '-- Dynamic model results'
+    print '--------------------------------------------------------------------'
+    if wgEnable:
+        print 'Apertif BF internal levels (excluding sign bit) for coherent WG sinus input:'
+        print '. ADC input rms                = %20.2f               (= %4.1f [bit])' % (rmsSp, math.log(rmsSp, 2))
+        print '. ADC input ampl               = %20.2f               (= %4.1f [bit])' % (amplSp, math.log(amplSp, 2))
+        print '. SST input rms complex        = %20.2f               (= %4.1f [bit])' % (rmsComplexSST, math.log(rmsComplexSST, 2))
+        print '. SST input ampl               = %20.2f               (= %4.1f [bit])' % (amplSST, math.log(amplSST, 2))
+        print '. BST input rms complex        = %20.2f               (= %4.1f [bit])' % (rmsComplexBST, math.log(rmsComplexBST, 2))
+        print '. BST input ampl               = %20.2f               (= %4.1f [bit])' % (amplBST, math.log(amplBST, 2))
+        print '. CB beamlet rms complex       = %20.2f               (= %4.1f [bit])' % (rmsComplexBeamlet, math.log(rmsComplexBeamlet, 2))
+        print '. CB beamlet ampl              = %20.2f               (= %4.1f [bit])' % (amplBeamlet, math.log(amplBeamlet, 2))
+    else:
+        print 'Apertif BF internal levels (excluding sign bit) for incoherent ADC noise input:'
+        print '. ADC input rms                = %20.2f               (= %4.1f [bit])' % (rmsSp, math.log(rmsSp, 2))
+        print '. SST input rms complex        = %20.2f               (= %4.1f [bit])' % (rmsComplexSST, math.log(rmsComplexSST, 2))
+        print '. SST input rms                = %20.2f               (= %4.1f [bit])' % (rmsSST, math.log(rmsSST, 2))
+        print '. BST input rms complex        = %20.2f               (= %4.1f [bit])' % (rmsComplexBST, math.log(rmsComplexBST, 2))
+        print '. BST input rms                = %20.2f               (= %4.1f [bit])' % (rmsBST, math.log(rmsBST, 2))
+        print '. CB beamlet rms complex       = %20.2f               (= %4.1f [bit])' % (rmsComplexBeamlet, math.log(rmsComplexBeamlet, 2))
+        print '. CB beamlet rms               = %20.2f               (= %4.1f [bit])' % (rmsBeamlet, math.log(rmsBeamlet, 2))
+    print '. SST power                    = %20.2f  = %7.2f dB (= %4.1f [bit])' % (SSTpower, 10*math.log(SSTpower, 10), math.log(SSTpower, 2))
+    print '. BST power                    = %20.2f  = %7.2f dB (= %4.1f [bit])' % (BSTpower, 10*math.log(BSTpower, 10), math.log(BSTpower, 2))
+    print ''
+    
+    if 'apertif' in apps:
+        if wgEnable:
+            print 'Apertif XC output for coherent WG sinus input:'
+            print '. Fine channel rms complex     = %20.2f               (= %4.1f [bit])' % (rmsComplexFineChannel, math.log(rmsComplexFineChannel, 2))
+            print '. Fine channel ampl            = %20.2f               (= %4.1f [bit])' % (amplFineChannel, math.log(amplFineChannel, 2))
+        else:
+            print 'Apertif XC output for incoherent ADC noise input:'
+            print '. Fine channel rms complex     = %20.2f               (= %4.1f [bit])' % (rmsComplexFineChannel, math.log(rmsComplexFineChannel, 2))
+            print '. Fine channel rms             = %20.2f               (= %4.1f [bit])' % (rmsFineChannel, math.log(rmsFineChannel, 2))
+        print '. XC auto power                = %20.2f  = %7.2f dB (= %4.1f [bit])' % (XCpower, 10*math.log(XCpower, 10), math.log(XCpower, 2))
+        print ''
+    
+    if 'arts_sc1' in apps:
+        if wgEnable:
+            print 'Arts SC1 beamlet TAB internal levels (excluding sign bit) for coherent WG sinus input:'
+            print '. Output TAB rms complex       = %20.2f               (= %4.1f [bit])' % (rmsComplexBeamletTab, math.log(rmsComplexBeamletTab, 2))
+            print '. Output TAB ampl              = %20.2f               (= %4.1f [bit])' % (amplBeamletTab, math.log(amplBeamletTab, 2))
+        else:
+            print 'Arts SC1 beamlet TAB internal levels (excluding sign bit) for incoherent ADC noise input:'
+            print '. Output TAB rms complex       = %20.2f               (= %4.1f [bit])' % (rmsComplexBeamletTab, math.log(rmsComplexBeamletTab, 2))
+            print '. Output TAB rms               = %20.2f               (= %4.1f [bit])' % (rmsBeamletTab, math.log(rmsBeamletTab, 2))
+        print ''
+    
+    if 'arts_sc4' in apps:
+        if wgEnable:
+            print 'Arts SC4 TAB internal levels (excluding sign bit) for coherent WG sinus input:'
+            print '. Fine channel rms complex     = %20.2f               (= %4.1f [bit])' % (rmsComplexFineChannel, math.log(rmsComplexFineChannel, 2))
+            print '. Fine channel ampl            = %20.2f               (= %4.1f [bit])' % (amplFineChannel, math.log(amplFineChannel, 2))
+            print '. Fine channel TAB rms complex = %20.2f               (= %4.1f [bit])' % (rmsComplexFineChannelTab, math.log(rmsComplexFineChannelTab, 2))
+            print '. Fine channel TAB ampl        = %20.2f               (= %4.1f [bit])' % (amplFineChannelTab, math.log(amplFineChannelTab, 2))
+            print '. Fine channel TAB power       = %20.2f               (= %4.1f [bit])' % (meanFineChannelTabPower, math.log(meanFineChannelTabPower, 2))
+            print '. Output TAB power             = %20.2f               (= %4.1f [bit])' % (meanIntegratedTabPower, math.log(meanIntegratedTabPower, 2))
+        else:
+            print 'Arts SC4 TAB internal levels (excluding sign bit) for incoherent ADC noise input:'
+            print '. Fine channel rms complex     = %20.2f               (= %4.1f [bit])' % (rmsComplexFineChannel, math.log(rmsComplexFineChannel, 2))
+            print '. Fine channel rms             = %20.2f               (= %4.1f [bit])' % (rmsFineChannel, math.log(rmsFineChannel, 2))
+            print '. Fine channel TAB rms complex = %20.2f               (= %4.1f [bit])' % (rmsComplexFineChannelTab, math.log(rmsComplexFineChannelTab, 2))
+            print '. Fine channel TAB rms         = %20.2f               (= %4.1f [bit])' % (rmsFineChannelTab, math.log(rmsFineChannelTab, 2))
+            print '. Fine channel TAB power       = %20.2f               (= %4.1f [bit])' % (meanFineChannelTabPower, math.log(meanFineChannelTabPower, 2))
+            print '. Fine channel TAB power std   = %20.2f               (= %4.1f [bit])' % (stdFineChannelTabPower, math.log(stdFineChannelTabPower, 2))
+            print '. Output TAB power             = %20.2f               (= %4.1f [bit])' % (meanIntegratedTabPower, math.log(meanIntegratedTabPower, 2))
+            print '. Output TAB power std         = %20.2f               (= %4.1f [bit])' % (stdIntegratedTabPower, math.log(stdIntegratedTabPower, 2))
+        print ''
+        
+        if wgEnable:
+            print 'Arts SC4 IAB internal levels (excluding sign bit) for coherent WG sinus input:'
+            print '. Fine channel rms complex     = %20.2f               (= %4.1f [bit])' % (rmsComplexFineChannel, math.log(rmsComplexFineChannel, 2))
+            print '. Fine channel ampl            = %20.2f               (= %4.1f [bit])' % (amplFineChannel, math.log(amplFineChannel, 2))
+            print '. Fine channel power           = %20.2f               (= %4.1f [bit])' % (meanFineChannelPower, math.log(meanFineChannelPower, 2))
+            print '. Integrated channel power     = %20.2f               (= %4.1f [bit])' % (meanIntegratedChannelPower, math.log(meanIntegratedChannelPower, 2))
+            print '. Integrated IAB power         = %20.2f               (= %4.1f [bit])' % (meanIntegratedIabPower, math.log(meanIntegratedIabPower, 2))
+            print '. Output IAB power             = %20.2f               (= %4.1f [bit])' % (meanOutputIabPower, math.log(meanOutputIabPower, 2))
+        else:
+            print 'Arts SC4 IAB internal levels (excluding sign bit) for incoherent ADC noise input:'
+            print '. Fine channel rms complex     = %20.2f               (= %4.1f [bit])' % (rmsComplexFineChannel, math.log(rmsComplexFineChannel, 2))
+            print '. Fine channel rms             = %20.2f               (= %4.1f [bit])' % (rmsFineChannel, math.log(rmsFineChannel, 2))
+            print '. Fine channel power           = %20.2f               (= %4.1f [bit])' % (meanFineChannelPower, math.log(meanFineChannelPower, 2))
+            print '. Integrated channel power     = %20.2f               (= %4.1f [bit])' % (meanIntegratedChannelPower, math.log(meanIntegratedChannelPower, 2))
+            print '. Integrated IAB power         = %20.2f               (= %4.1f [bit])' % (meanIntegratedIabPower, math.log(meanIntegratedIabPower, 2))
+            print '. Output IAB power             = %20.2f               (= %4.1f [bit])' % (meanOutputIabPower, math.log(meanOutputIabPower, 2))
+            print '. Output IAB power std         = %20.2f               (= %4.1f [bit])' % (stdOutputIabPower, math.log(stdOutputIabPower, 2))
+        print ''
+
+################################################################################
+# Plotting
+################################################################################    
+if useplot:
+    tplot = test_plot.Testplot()
+    tplot.close_plots()
+    
+    if 'static' in models:
+        ################################################################################
+        # Static model:
+        # . WG, coherent input
+        # . Noise, incoherent input
+        ################################################################################
+        
+        if wgEnable:
+            # WG, coherent input
+            # Plot internal levels as function of WG amplitude
+            Lx = A_wg
+            if 'apertif' in apps or 'arts_sc4' in apps:
+                A = []
+                A.append(A_wg)
+                A.append(beamlet_wg_complex_phase_values * pi_apr.wgScale)
+                A.append(channel_x_wg_complex_phase_values * pi_apr.wgScale)
+                Alegend = ['adc', 'beamlet', 'channel_x']
+                lineFormats = ['-', '-', '-']
+                Title = 'Amplitudes for cbWeight=%d' % cbWeight
+                Ylim = None
+                tplot.plot_two_dimensional_list(3, A, Alegend=Alegend, Lx=Lx, representation='', offset=0, Title=Title, Xlabel='WG amplitude', Ylabel='Level', Xlim=None, Ylim=Ylim, lineFormats=lineFormats)
+                tplot.save_figure('plots/apertif_arts_firmware_model_static_wg_amplitudes.png')   # Save current figure
+
+            if 'arts_sc1' in apps:
+                A = []
+                A.append(A_wg)
+                A.append(beamlet_wg_complex_phase_values * pi_apr.wgScale)
+                A.append(beamlet_tab_wg_complex_phase_values_sc1 * pi_apr.wgScale)
+                Alegend = ['adc', 'beamlet', 'tab sc1']
+                lineFormats = ['-', '-', '*']
+                Title = 'Amplitudes for cbWeight=%d, tabWeightSc1=%d, nofDishes=%d' % (cbWeight, tabWeightSc1, nofDishes)
+                Ylim = None
+                tplot.plot_two_dimensional_list(3, A, Alegend=Alegend, Lx=Lx, representation='', offset=0, Title=Title, Xlabel='WG amplitude', Ylabel='Level', Xlim=None, Ylim=Ylim, lineFormats=lineFormats)
+                tplot.save_figure('plots/apertif_arts_firmware_model_static_wg_amplitudes_sc1.png')   # Save current figure
+        
+            if 'arts_sc4' in apps:
+                A = []
+                A.append(power_tab_xx_wg_complex_values_sc4)
+                A.append(iab_xx_wg_complex_values_sc4)
+                Alegend = ['TAB XX* power', 'IAB XX* power']
+                lineFormats = ['-', '*']
+                Title = 'Powers for cbWeight=%d, tabWeightSc4=%d, iabGainReg=%d, nofDishes=%d' % (cbWeight, tabWeightSc4, iabGainReg, nofDishes)
+                Ylim = [0, 1.1 * 2.0**pi_apr.W_pow]
+                tplot.plot_two_dimensional_list(3, A, Alegend=Alegend, Lx=Lx, representation='', offset=0, Title=Title, Xlabel='WG amplitude', Ylabel='Level', Xlim=None, Ylim=Ylim, lineFormats=lineFormats)
+                tplot.save_figure('plots/apertif_arts_firmware_model_static_wg_powers_sc4.png')   # Save current figure
+            
+            # Plot power statistics as function of WG amplitude
+            A = []
+            A.append(SST_wg_dBs)
+            A.append(BST_wg_dBs)
+            if 'apertif' in apps:
+                A.append(XC_wg_auto_visibilities_dB)
+                Alegend = ['SST', 'BST', 'XC']
+            else:
+                Alegend = ['SST', 'BST']
+            Title = 'Auto correlation power statistics for CB weight = %d' % cbWeight
+            Ylim = None
+            tplot.plot_two_dimensional_list(3, A, Alegend=Alegend, Lx=Lx, representation='', offset=0, Title=Title, Xlabel='WG amplitude', Ylabel='Level [dB]', Xlim=None, Ylim=Ylim, lineFormats=None)
+            tplot.save_figure('plots/apertif_arts_firmware_model_static_wg_auto_powers.png')   # Save current figure
+        
+        else:
+            # Noise, incoherent input
+            # Plot internal levels as function of WG amplitude
+            Lx = S_noise
+            if 'apertif' in apps or 'arts_sc4' in apps:
+                A = []
+                A.append(S_noise)
+                A.append(beamlet_noise_complex_phase_values)
+                A.append(channel_x_noise_complex_phase_values)
+                Alegend = ['adc', 'beamlet', 'channel_x']
+                lineFormats = ['-', '*', '-']
+                Title = 'Sigmas for cbWeight=%d' % cbWeight
+                Ylim = None
+                tplot.plot_two_dimensional_list(3, A, Alegend=Alegend, Lx=Lx, representation='', offset=0, Title=Title, Xlabel='ADC noise', Ylabel='Level', Xlim=None, Ylim=Ylim, lineFormats=lineFormats)
+                tplot.save_figure('plots/apertif_arts_firmware_model_static_noise_sigmas.png')   # Save current figure
+        
+            if 'arts_sc1' in apps:
+                A = []
+                A.append(S_noise)
+                A.append(beamlet_noise_complex_phase_values)
+                A.append(beamlet_tab_noise_complex_phase_values_sc1)
+                Alegend = ['adc', 'beamlet', 'tab sc1']
+                lineFormats = ['-', '-', '*']
+                Title = 'Sigmas for cbWeight=%d, tabWeightSc1=%d, nofDishes=%d' % (cbWeight, tabWeightSc1, nofDishes)
+                Ylim = None
+                tplot.plot_two_dimensional_list(3, A, Alegend=Alegend, Lx=Lx, representation='', offset=0, Title=Title, Xlabel='ADC noise', Ylabel='Level', Xlim=None, Ylim=Ylim, lineFormats=lineFormats)
+                tplot.save_figure('plots/apertif_arts_firmware_model_static_noise_sigmas_sc1.png')   # Save current figure
+
+            if 'arts_sc4' in apps:
+                A = []
+                A.append(power_tab_xx_noise_values_sc4)
+                A.append(iab_xx_noise_values_sc4)
+                Alegend = ['TAB XX* power', 'IAB XX* power']
+                lineFormats = ['-', '*']
+                Title = 'Powers for cbWeight=%d, tabWeightSc4=%d, iabGainReg=%d, nofDishes=%d' % (cbWeight, tabWeightSc4, iabGainReg, nofDishes)
+                Ylim = [0, 1.1 * 2.0**pi_apr.W_pow]
+                tplot.plot_two_dimensional_list(3, A, Alegend=Alegend, Lx=Lx, representation='', offset=0, Title=Title, Xlabel='ADC noise', Ylabel='Level', Xlim=None, Ylim=Ylim, lineFormats=lineFormats)
+                tplot.save_figure('plots/apertif_arts_firmware_model_static_noise_powers_sc4.png')   # Save current figure
+            
+            # Plot power statistics as function of WG amplitude
+            A = []
+            A.append(SST_noise_dBs)
+            A.append(BST_noise_dBs)
+            A.append(XC_noise_auto_visibilities_dB)
+            Alegend = ['SST', 'BST', 'XC']
+            Title = 'Auto correlation power statistics for CB weight = %d' % cbWeight
+            Ylim = None
+            tplot.plot_two_dimensional_list(3, A, Alegend=Alegend, Lx=Lx, representation='', offset=0, Title=Title, Xlabel='ADC noise', Ylabel='Level [dB]', Xlim=None, Ylim=Ylim, lineFormats=None)
+            tplot.save_figure('plots/apertif_arts_firmware_model_static_noise_auto_powers.png')   # Save current figure
+
+        
+    if 'dynamic' in models and wgSweep == 0:
+        ###############################################################################
+        # Dynamic model: Processing
+        ###############################################################################
+    
+        # Plot dish ADC
+        L = spSamples[0:nof_periods]  # do not plot all to save time
+        if wgEnable:
+            Title = 'Dish ADC input (sigma = %.2f, ampl = %.2f)' % (rmsSp, wgAmpl)
+            ylim = 5 * math.ceil(0.2 * 1.1 * amplSp)
+        else:
+            Title = 'ADC input (sigma = %.2f)' % rmsSp
+            ylim = 5 * math.ceil(0.2 * 1.1 * rmsSp * nofSigma)
+        Ylim = [-ylim, ylim]
+        tplot.plot_one_dimensional_list(3, L, Lx=None, representation='', offset=0, Title=Title, Xlabel='Time [Tsub]', Ylabel='Voltage', Xlim=None, Ylim=Ylim)
+        tplot.save_figure('plots/apertif_arts_firmware_model_dynamic_adc.png')   # Save current figure
+        
+        # Plot dish output beamlets for selSub
+        A = []
+        A.append(beamletSamples.real[0:nof_periods])
+        A.append(beamletSamples.imag[0:nof_periods])
+        Alegend = ['real', 'imag']
+        if wgEnable:
+            Title = 'Dish CB output samples (ampl = %.2f)' % amplBeamlet
+            ylim = 5 * math.ceil(0.2 * 1.1 * amplBeamlet)
+        else:
+            Title = 'Dish CB output samples (sigma = %.2f)' % rmsBeamlet
+            ylim = 5 * math.ceil(0.2 * 1.1 * rmsBeamlet * nofSigma)
+        Ylim = [-ylim, ylim]
+        tplot.plot_two_dimensional_list(3, A, Alegend=Alegend, Lx=None, representation='', offset=0, Title=Title, Xlabel='time [Tsub]', Ylabel='Level', Xlim=None, Ylim=Ylim, lineFormats=None)
+        tplot.save_figure('plots/apertif_arts_firmware_model_dynamic_beamlet.png')   # Save current figure
+        
+        if 'apertif' in apps or 'arts_sc4' in apps:
+            # Plot correlator fine channels for selSub, selChanx
+            A = []
+            A.append(fineChannelSamples.real[0:nof_periods])
+            A.append(fineChannelSamples.imag[0:nof_periods])
+            Alegend = ['real', 'imag']
+            if wgEnable:
+                Title = 'Central fine channel samples (ampl = %.2f)' % amplFineChannel
+                ylim = 5 * math.ceil(0.2 * 1.1 * amplFineChannel)
+            else:
+                Title = 'Central fine channel samples (sigma = %.2f)' % rmsFineChannel
+                ylim = 5 * math.ceil(0.2 * 1.1 * rmsFineChannel * nofSigma)
+            Ylim = [-ylim, ylim]
+            tplot.plot_two_dimensional_list(3, A, Alegend=Alegend, Lx=None, representation='', offset=0, Title=Title, Xlabel='time [Tsub]', Ylabel='Level', Xlim=None, Ylim=Ylim, lineFormats=None)
+            tplot.save_figure('plots/apertif_arts_firmware_model_dynamic_fine_channel.png')   # Save current figure
+        
+        if 'arts_sc1' in apps:
+            # Plot Arts SC1 output TAB samples for selSub and nofDishes
+            # . Independent of nofDishes when tabWeightSc1 is scaled by nofDishes or sqrt(nofDishes)
+            A = []
+            A.append(beamletTabSamples.real[0:nof_periods])
+            A.append(beamletTabSamples.imag[0:nof_periods])
+            Alegend = ['real', 'imag']
+            if wgEnable:
+                Title = 'Arts SC1 output TAB samples (ampl = %.2f, nofDishes = %d)' % (amplBeamletTab, nofDishes)
+                ylim = 5 * math.ceil(0.2 * 1.1 * amplBeamletTab)
+            else:
+                Title = 'Arts SC1 output TAB samples (sigma = %.2f, nofDishes = %d)' % (rmsBeamletTab, nofDishes)
+                ylim = 5 * math.ceil(0.2 * 1.1 * rmsBeamletTab * nofSigma)
+            Ylim = [-ylim, ylim]
+            tplot.plot_two_dimensional_list(3, A, Alegend=Alegend, Lx=None, representation='', offset=0, Title=Title, Xlabel='time [Tsub]', Ylabel='Level', Xlim=None, Ylim=Ylim, lineFormats=None)
+            tplot.save_figure('plots/apertif_arts_firmware_model_dynamic_beamlet_tab.png')   # Save current figure
+    
+        if 'arts_sc4' in apps:
+            # Plot Arts SC4 output TAB powers for selSub, selChanx and nofDishes
+            # . Independent of nofDishes, because tabWeightSc4 is scaled by nofDishes or sqrt(nofDishes)
+            A = []
+            A.append(integratedTabPowers[0:nof_periods])
+            Alegend = None
+            Title = 'Arts SC4 output TAB powers (mean = %.2f = %.1f [bit], nofDishes = %d)' % (meanIntegratedTabPower, math.log(meanIntegratedTabPower, 2), nofDishes)
+            if wgEnable:
+                ylim = 5 * math.ceil(0.2 * 1.1 * meanIntegratedTabPower)
+            else:
+                ylim = 5 * math.ceil(0.2 * 1.1 * meanIntegratedTabPower * nofSigma)
+            Ylim = [0, ylim]
+            tplot.plot_two_dimensional_list(3, A, Alegend=Alegend, Lx=None, representation='', offset=0, Title=Title, Xlabel='time [Tchanx]', Ylabel='Level', Xlim=None, Ylim=Ylim, lineFormats=None)
+            tplot.save_figure('plots/apertif_arts_firmware_model_dynamic_power_tab.png')   # Save current figure
+            
+            # Plot Arts SC4 output IAB powers for selSub, selChanx and nofDishes
+            # . Independent of nofDishes, because tabWeightSc4 is scaled by nofDishes or sqrt(nofDishes)
+            A = []
+            A.append(outputIabPowers)
+            Alegend = None
+            Title = 'Arts SC4 output IAB powers (mean = %.2f = %.1f [bit], nofDishes = %d)' % (meanOutputIabPower, math.log(meanOutputIabPower, 2), nofDishes)
+            if wgEnable:
+                ylim = 5 * math.ceil(0.2 * 1.1 * meanOutputIabPower)
+            else:
+                ylim = 5 * math.ceil(0.2 * 1.1 * meanOutputIabPower * nofSigma)
+            Ylim = [0, ylim]
+            tplot.plot_two_dimensional_list(3, A, Alegend=Alegend, Lx=None, representation='', offset=0, Title=Title, Xlabel='time [Tchanx]', Ylabel='Level', Xlim=None, Ylim=Ylim, lineFormats=None)
+            tplot.save_figure('plots/apertif_arts_firmware_model_dynamic_power_iab.png')   # Save current figure
+        
+        
+    if 'dynamic' in models and wgSweep != 0:
+        ###############################################################################
+        # Dynamic model: WG frequency sweep
+        ###############################################################################
+
+        # Plot subband spectrum abs(fSubSamples[nofTsub][N_sub])
+        A = np.abs(fSubSamples)
+        tplot.image_two_dimensional_list(3, A, transpose=True, grid=True, Title='Subband magnitudes', Xlabel='time [Tsub]', Ylabel='freq [fsub]')
+        tplot.save_figure('plots/apertif_arts_firmware_model_wg_sweep_subband_spectrum.png')   # Save current figure
+        
+        # Plot fine channels spectrum abs(fSubChanSamples[N_sub][nofTchanx][N_chan_x])
+        A = np.abs(fSubChanSamples.reshape(nofTchanx, N_sub*N_chan_x))   # [nofTchanx][N_sub*N_chan_x]
+        chanLo = 0
+        chanHi = N_sub*N_chan_x - 1
+        nofSubBorder = 0
+        nofSubBorder = 1 + int(abs(wgSweep))
+        if nofSubBorder>0:
+            if wgSweep>0:
+                chanLo = N_chan_x * (selSub - nofSubBorder)
+                chanHi = N_chan_x * (selSub + nofSubBorder + wgSweep + 1) - 1
+            else:
+                chanLo = N_chan_x * (selSub - nofSubBorder + wgSweep)
+                chanHi = N_chan_x * (selSub + nofSubBorder + 1) - 1
+            A = A[:, chanLo:chanHi+1]
+        extent = [0, nofTchanx-1, 1.0 * chanLo / N_chan_x - 0.5, 1.0 * chanHi / N_chan_x - 0.5]
+        #print chanLo, chanHi, extent
+        #extent = None
+        tplot.image_two_dimensional_list(3, A, transpose=True, grid=True, Title='Channel magnitudes', Xlabel='time [Tchanx]', Ylabel='freq [fsub]', extent=extent)
+        tplot.save_figure('plots/apertif_arts_firmware_model_wg_sweep_fine_channel_spectrum.png')   # Save current figure
+    
+    tplot.show_plots()
+

applications/lofar2/model/pfb_pk/PFB.py

+#!/usr/bin/env python
+# coding: utf-8
+
+# In[1]:
+
+
+import numpy as np
+import matplotlib.pyplot as plt
+
+
+# In[8]:
+
+
+#Calculate the weight of the filter bank
+def GetW(Nfft,P=4):
+    M=Nfft*P
+    W=np.arange(0,M,1)
+    flt=np.zeros([M])
+    #P=P-1
+    flt[0:P-2]=1;
+    flt[P-2]=0.5;
+    flt[-P+3:]=1;
+    flt[-P+2]=0.5;
+    W=np.fft.ifft(flt).real #Imag part should be zero
+    W=np.roll(W,M//2)
+    return W*Nfft
+
+#Apply the PFB to data
+def fftw(A,W,P=4):
+  Nfft=A.shape[-1]
+  Dx=A[:,:-P,:]*W[:Nfft]
+  for i in range(P-1):
+    S=(i+1)*Nfft
+    Dx+=A[:,i+1:-P+i+1,:]*W[S:S+Nfft]
+  return (np.fft.fft(Dx,axis=2))
+
+
+# In[11]:
+
+
+P=4
+Nfft=1024; #FFT length
+W=GetW(Nfft,P=P)
+plt.plot(W);
+plt.title("PFB coefficients");
+
+
+# In[22]:
+
+
+#Load Data
+Nc=3; #Number of channels
+
+A=np.array(np.load("../LabTest8/D6b.npy")) #Raw ADC data
+A-=65536*(A>65536//2)
+print(A.shape)
+[n1,n2,n3]=A.shape
+
+A=A.reshape([n1,n3//Nfft,Nfft])
+print(A.shape) #Channels*seconds, blocks, fft length
+
+
+# In[ ]:
+
+
+
+
+
+# In[23]:
+
+
+#Do FFT
+F=fftw(A,W,P=P)
+[n1,n2,n3]=F.shape
+F=F.reshape([n1//Nc,Nc,n2,n3])
+F=F.swapaxes(0,1)
+F=F.reshape([Nc,n1//Nc*n2,n3])
+print(F.shape) #channels, blocks, ffts
+
+
+# In[26]:
+
+
+#plot power
+FP=np.sum((F*F.conj()).real,axis=1)[:,:Nfft//2]
+FP=10*np.log10(FP)
+freq=np.arange(0,Nfft//2,1)/(Nfft//2)*100
+plt.figure(dpi=300)
+plt.plot(freq,FP.T,linewidth=0.5);
+plt.xlabel("Frequency (MHz)")
+plt.ylabel("Power (dB)")
+
+
+# In[ ]:
+
+
+
+

applications/lofar2/model/pfb_pk/readme.txt

+2021-01-20 Paulus Kruger script for offline analysis of LTS raw data and comparison with SST data.
+Use Jupyter notebook: https://jupyter.org/
+LTS results: https://support.astron.nl/confluence/display/RCU/2020-11-23+Stati+measurement+vs+raw+data