Derived from apertif_arts_firmware_model.py, but in a very intermediate state, so not working.

applications/lofar2/model/lofar_station_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 LOFAR2.0 Station Digital Processor (SDP) firmware data path signal levels
+
+[1] https://support.astron.nl/confluence/display/L2M/L4+SDP+Decision%3A+LOFAR2.0+SDP+Firmware+Quantization+Model
+
+The SDP FW processing is described in [1].
+The purpose of the lofar_station_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 = 14, W_subband = 18, W_beamlet = 8).
+ 
+This lofar_station_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 subband to station beamlet output.
+- Plot results of the static and the dynamic model.
+
+The model matches the drawings of the firmware parameters and signal levels in [1].
+
+  
+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 lofar_station_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
+
+################################################################################    
+# SDP settings
+################################################################################    
+
+
+W_adc             = 14           # Word width in number of bits of the ADC
+W_adc_max         = W_adc-1      # ADC full scale
+
+N_complex         = 2
+N_fft             = 1024
+N_sub             = 512          # = N_fft / N_complex
+W_fsub_proc       = math.log(math.sqrt(N_sub), 2)    # = 4.5, processing gain of F_sub in number of bits
+
+N_int_sub         = 195312.5
+W_se_weight      = 16
+W_se_magnitude   = 1
+W_se_fraction    = 14
+c_se_unit_weight = 2**W_se_fraction
+
+W_bf_weight = 16
+W_bf_magnitude = 1
+W_bf_fraction  = 14
+c_bf_unit_weight = 2**W_sub_fraction
+
+W_beamlet_scale = 16
+W_beamlet_scale_magnitude = 1
+W_beamlet_scale_fraction  = 14
+c_beamlet_scale_unit = 2**W_beamlet_scale_fraction
+
+
+W_sst_fraction  =      4        # Extra fraction bits for SST subband input compared to W_adc
+W_sst_in        =     18        # Word width in number of bits of SST subband data input, = W_adc + W_sst_fraction
+W_sst_out       =     54        # Word width in number of bits of SST subband statistics output, = 2*W_sst_in + log2(N_int_sub)
+
+W_bf_in         =     18        # Word width in number of bits of BF subband data input
+W_bf_weight     =     16        # Word width in number of bits of the BF weights
+W_cb_max_weight = W_cb_weight-1 # = 15 = 16 - 1
+W_cb_unit_scale =     -2        # Unit weight scale factor in number of bits
+W_cb_weight_fraction = W_cb_max_weight + W_cb_unit_scale
+cbMaxWeight     = 2.0**W_cb_max_weight - 1     # = 32767, maximum positive BF weight value in full scale range
+cbUnitWeight    = 2.0**W_cb_weight_fraction    # = 8192, unit BF weight value mapped in full scale range
+
+W_bst_fraction  =      8        # Extra fraction bits for BST beamlet input compared to W_adc, = W_fsub_gain - W_cb_unit_scale + 1
+W_bst_in        =     16        # Word width in number of bits of BST beamlet data input, = W_adc + W_bst_fraction
+W_bst_out       =     56        # Word width in number of bits of BST beamlet statistics output
+W_bst_out_sz    =      2        # Number of 32b words per W_bst_out
+
+W_beamlet_sc1        = 8        # Word width in number of bits of a beamlet voltage sample at the BF to Arts SC1 interface
+W_beamlet_max_sc1    = W_beamlet_sc1-1         # = 7 = 8 - 1
+W_beamlet            = 8        # Word width in number of bits of a beamlet voltage sample at the BF to XC and Arts SC3 and SC4 interface
+W_beamlet_max        = W_beamlet-1             # = 7 = 8 - 1 or = 5 = 6 - 1
+W_beamlet_fraction   = 6        # Extra fraction bits for beamlet output compared to W_adc, = W_fsub_gain - W_cb_unit_scale - 1,
+                                # do -1 to have 1 bit extra for noise at ADC interface than at BF beamlet output interface
+
+W_channel_x          = 9        # Word width in number of bits of XC input channel voltage sample
+W_channel_x_max      = W_channel_x-1           # = 8 = 9 - 1
+W_channel_x_fraction = 9        # Extra fraction bits for XC input channel compared to W_adc, = W_beamlet_fraction + W_fchan_x_gain
+
+W_tab_weight_sc1     = 16       # Word width in number of bits of the TAB weights for SC1
+W_tab_max_weight_sc1 = W_tab_weight_sc1-1               # = 15 = 16 - 1
+tabMaxWeightSc1      = 2.0**W_tab_max_weight_sc1 - 1    # = 32767 = 2**(16-1) - 1, do -1 to avoid wrap to -32768
+
+W_int_ax             = cm.ceil_log2(M_int_ax)  # = 4, number of bits to fit integration of M_int_ax = 16 power samples
+
+W_tab_weight_sc4     = 9        # Word width in number of bits of the TAB weights for SC3 and SC4, required is >= 5 bits complex
+W_tab_max_weight_sc4 = W_tab_weight_sc4-1               # = 8 = 9 - 1
+W_tab_unit_scale_sc4 = -1                               # Unit weight scale factor in number of bits
+W_tab_weight_fraction_sc4 = W_tab_max_weight_sc4 + W_tab_unit_scale_sc4
+tabMaxWeightSc4      = 2.0**W_tab_max_weight_sc4 - 1    # = 255, maximum positive TAB SC4 weight value in full scale range
+tabUnitWeightSc4     = 2.0**W_tab_weight_fraction_sc4   # = 128, unit TAB SC4 weight value mapped in full scale range
+
+W_iab_gain           = 16       # Word width in number of bits of the IAB outputgain for SC3 and SC4
+W_iab_unit_scale     = -2       # = -2
+W_iab_gain_fraction  = W_iab_gain-1 + W_iab_unit_scale  # = 13 = 16-1 + -2
+iabMaxGain           = 2.0**(W_iab_gain-1) - 1          # = 32767, maximum positive IAB gain value in full scale range
+iabUnitGain          = 2.0**W_iab_gain_fraction         # = 8192 = 2**13
+
+                                # Determine signal scale factors relative to ADC sigma level
+                                # . apply factor 2**W_cb_unit_scale to account for unit level of the BF weights relative to full scale
+wgSstScale          = 1.0 / 2.0**W_sst_fraction
+wgBstScale          = 1.0 / 2.0**(W_bst_fraction + W_cb_unit_scale)
+wgBeamletScale      = 1.0 / 2.0**(W_beamlet_fraction + W_cb_unit_scale)
+wgChannelXcorScale  = 1.0 / 2.0**(W_channel_x_fraction + W_cb_unit_scale)
+
+noiseSstScale         = 2.0**W_fsub_gain * wgSstScale
+noiseBstScale         = 2.0**W_fsub_gain * wgBstScale
+noiseBeamletScale     = 2.0**W_fsub_gain * wgBeamletScale
+noiseChannelXcorScale = 2.0**W_fsub_gain * wgChannelXcorScale * 2**W_fchan_x_gain
+noiseSubbandXcorScale = 2.0**W_fsub_gain * wgChannelXcorScale   # for XC subband, is sum of all XC channels
+
+wgScale         = math.sqrt(2)  # Apply factor wgScale = sqrt(2) to account for WG ampl = sqrt(2) * WG sigma, because WG power = A**2 / 2
+rcScale         = math.sqrt(2)  # Apply factor rcScale = sqrt(2) to account for rms real = rms imag = rms complex / sqrt(2)
+
+
+################################################################################    
+# Model settings
+################################################################################    
+
+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
+
+################################################################################    
+# Parse arguments to derive user settings
+################################################################################    
+
+_parser = argparse.ArgumentParser('lofar_station_firmware_model')
+_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 subband periods in dynamic model')
+_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('--sub_weight_adjust', default=1.0, type=float, help='Factor to adjust unit subband weight')
+_parser.add_argument('--bf_weight_adjust', default=1.0, type=float, help='Factor to adjust unit BF weight')
+args = _parser.parse_args()
+
+if args.model=='all':
+    models = ['static', 'dynamic']
+else:
+    models = [args.model]
+
+skySigma  = args.sky_sigma
+#skySigma  = 16          # 16 for 4 bit noise
+#skySigma  = 5792        # 5792 * 2**0.5 = 8191 for full scale WG amplitude
+
+wgAmpl    = args.wg_ampl
+wgSigma   = wgAmpl / np.sqrt(2)
+
+wgEnable  = args.wg_enable
+wgSweep   = args.wg_sweep
+if wgSweep != 0:
+    wgEnable = True   # force WG
+
+quantize  = args.quantize
+
+useplot   = args.useplot
+
+subWeightAdjust     = args.sub_weight_adjust
+bfWeightAdjust     = args.bf_weight_adjust
+#subWeightAdjust     = 1.0/8    # use <= 1/8 to have no beamlet (8 bit) clipping
+
+# Time series
+if 'dynamic' in models:
+    # use args.nof_periods for number of subband periods to simulate:
+    if wgEnable:
+        nofTsub = args.nof_periods
+        #nofTsub = 10          # use >~10 for WG
+        #nofTsub = 100
+    else:
+        nofTsub = args.nof_periods
+        #nofTsub = 100         # number of subband periods to simulate, use >~ 100 for noise
+        #nofTsub = 1000        # use >= 1000 for more accuracy
+        #nofTsub = 250        # use < 1000 to speed up simulation
+
+    # also use args.nof_periods for number of Ts, Tsub periods to plot
+    # . to speed up plotting for Ts
+    # . 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
+    selSub           = 65
+
+################################################################################
+# Derived gains
+
+# Aperif BF subband CB voltage weight
+subWeight = c_sub_unit_weight * subWeightAdjust      # = 16384 = 2**14 * 1.0
+bfWeight = c_bf_unit_weight * bfWeightAdjust         # = 16384 = 2**14 * 1.0
+if quantize:
+    subWeight = np.round(subWeight)
+    bfWeight = np.round(bfWeight)
+
+
+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**W_adc_max, inputResolution)
+    S_wg = A_wg / wgScale
+    S_wg_complex = S_wg
+    S_wg_complex_phase = S_wg_complex / rcScale
+    A_wg_complex_phase = S_wg_complex
+    S_noise = np.arange(inputResolution, 2**W_adc_max / nofSigma, inputResolution)
+    S_noise_complex = S_noise
+    S_noise_complex_phase = S_noise / 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 '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
+    
+
+    ################################################################################    
+    # Dynamic model time series
+    ################################################################################    
+    
+    # Time, frequency axis
+    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
+    
+    
+    ################################################################################
+    # 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()
+