From 0ae8fb557eb50b4e17b380609f65c1d07dc12834 Mon Sep 17 00:00:00 2001
From: Eric Kooistra <kooistra@astron.nl>
Date: Mon, 1 Feb 2021 17:14:42 +0100
Subject: [PATCH] Derived from apertif_arts_firmware_model.py, but in a very
intermediate state, so not working.
---
.../model/lofar_station_firmware_model.py | 1035 +++++++++++++++++
1 file changed, 1035 insertions(+)
create mode 100644 applications/lofar2/model/lofar_station_firmware_model.py
diff --git a/applications/lofar2/model/lofar_station_firmware_model.py b/applications/lofar2/model/lofar_station_firmware_model.py
new file mode 100644
index 0000000000..11d1ecea59
--- /dev/null
+++ b/applications/lofar2/model/lofar_station_firmware_model.py
@@ -0,0 +1,1035 @@
+###############################################################################
+#
+# 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()
+
--
GitLab