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