Updated md files and checked that MATLAB and Pyhton scripts indeed run ok.

.gitignore

+*log
+__pycache__
+*pyc
+*#
\ No newline at end of file

LICENSE.md

+# apertif_matlab license
+
+The apertif_matlab package is open source and  available under the Apache Licence 2.0. The following text fragment shall be placed in comments at the top of each (new) source code file:
+
+```
+  Copyright [year] [name of copyright owner]
+
+  Licensed under the Apache License, Version 2.0 (the "License");
+  you may not use this file except in compliance with the License.
+  You may obtain a copy of the License at
+
+     http://www.apache.org/licenses/LICENSE-2.0
+
+  Unless required by applicable law or agreed to in writing, software
+  distributed under the License is distributed on an "AS IS" BASIS,
+  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+  See the License for the specific language governing permissions and
+  limitations under the License.
+```
+#### Replace ```<year>``` by the current year.
+
+#### For [ASTRON] replace ```Copyright [year] [name of copyright owner]``` by:
+```
+  Copyright <year> Stichting Nederlandse Wetenschappelijk Onderzoek Instituten,
+  ASTRON Netherlands Institute for Radio Astronomy
+```
+
+
+[ASTRON]:https://www.astron.nl
+
+
+

NOTICE.md

+Citation Notice version 1.0
+
+```
+This Citation Notice is part of the apertif_matlab software suite.
+
+Parties that use ASTRON Software resulting in papers and/or publications are requested to
+refer to the DOI(s) that correspond(s) to the version(s) of the ASTRON Software used:
+
+10.5281/zenodo.4748331
+
+Parties that use ASTRON Software for purposes that do not result in publications (e.g.
+commercial parties) are asked to inform ASTRON about their use of ASTRON Software, by
+sending an email to including the DOIs mentioned above in the message.
+```
+

README.md

+# Apertif MATLAB
+
+The apertif_matlab package contains MATLAB scripts and Python scripts that model the subband filterbank and the channel filterbank that are used the Apertif radio telescope of [ASTRON] [1].
+
+Both filterbanks are critically sampled, polyphase filterbanks (PFB). The subband filterbank is a wideband polyphase filterbank (WPFB), meaning that the input sample rate of 800 MSps is a wideband factor 4 higher than the processing rate of 200 MHz. Both filterbanks are implemented in VHDL and run on FPGAs of the UniBoard [2].
+
+For more information on the the apertif_matlab package see:
+
+* The matlab/apertif_matlab_readme.txt briefly decribes each script
+* The matlab/apertif_arts_firmware_model.pdf describes the quantization model of the subband filterbank and channel filterbank. The Python script apertif_arts_firmware_model.py runs the model.
+
+## References
+
+[1] W.A. van Cappellen et al., "Apertif, Phased Array Feeds for the Westerbork Synthesis Radio Telescope", unpublished, 2021
+
+[2] "The application of UniBoard as a beam former for APERTIF”, A. Gunst, A. Szomoru, G. Schoonderbeek, E. Kooistra, D. van der Schuur, H. Pepping, Experimental Astronomy, Volume 37, Issue 1, pp 55-67, February 2014, doi:1
+0.1007/s10686-013-9366-x.
+
+
+If you use (part of) the **apertif_matlab** package please attribute the use as indicated this citation NOTICE:
+
+### [NOTICE](NOTICE.md)
+
+The apertif_matlab package is Open Source and available under the following LICENSE:
+
+### [LICENSE](LICENSE.md)
+
+The apertif_matlab package was developed at [ASTRON] and used in several projects, but others have contributed and are welcome to contribute as well. The CREDITS lists the contributers of apertif_matlab:
+
+
+[ASTRON]:https://www.astron.nl

matlab/apertif_arts_firmware_model.py

@@ -141,7 +141,7 @@ import math
 import numpy as np
 import common as cm
 import test_plot
-import pi_apertif_system as pi_apr
+import pi_apertif_system_constants as pi_apr
 
 ################################################################################    
 # General

matlab/pi_apertif_system_constants.py

+###############################################################################
+#
+# Copyright (C) 2017
+# 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, 12 May 2017,
+
+"""Apertif system constants
+
+"""
+
+
+###############################################################################
+# System imports
+import math
+import common as cm
+
+
+###############################################################################
+# Apertif constants (from ASTRON-SP-062)
+# - The relation between constants is explained in the comment, but not used
+#   in the constant definition. The disadvantage of this is the dependencies
+#   between constants need to be maintained manually, however since they are
+#   unlikely to change this is no problem. The advantage is that the constant
+#   values are immediately clear in the code.
+# - Use bm6 to select between 6 bit beamlet mode or 8 bit beamlet mode
+bm6 = True                      # Init for 6 bit mode when True, else for 8 bit mode
+#bm6 = False
+
+N_complex       =      2        # Two part of a complex number, the real and imaginary part
+N_pol           =      2        # Number of polarizations, X and Y
+N_page          =      2        # Number of pages, e.g. 2 pages in a dual page memory
+N_Stokes        =      4        # Number of power values in the Stokes vector [I, Q, U, V]
+N_dish          =     12        # Number of WSRT dishes in Apertif
+K_dish          =      8        # Number of WSRT dishes in that are used for Arts SC3
+N_tp            =     24        # Number of single polarization telescope paths of the Apertif BF, = N_pol * N_dish
+K_tp            =     16        # Number of single polarization telescope paths for Arts SC3, = N_pol * K_dish
+N_ant_x         =     61        # Number of antennas in the frontend FPA X polarization of the Apertif BF
+N_ant_y         =     61        # Number of antennas in the frontend FPA Y polarization of the Apertif BF
+S_ant           =     64        # Number of ADC signal paths in the frontend FPA of the Apertif BF
+S_BN            =      4        # Number of ADC signal paths per BN in the frontend FPA of the Apertif BF
+N_BN            =     16        # Number of BN per subrack in the Apertif BF, equals the number of FN and = N_band
+N_FN            =     16        # Number of FN per subrack in the Apertif BF, equals the number of BN and = N_band
+T_sample        =      1.25e-9  # s, Digitizer sample period, = 1/ f_sample
+f_sample        =    800e6      # Hz Digitizer sample frequency of the ADC at the Apertif BF frontend
+f_clk           =    200e6      # Hz Data processing clock rate in the FPGA, = f_sample / P_sub
+P_sub           =      4        # Wideband parallelization factor for the subband filterbank, = f_sample / f_clk
+f_low           =     50e6      # Hz, Lower edge frequency of the 400 MHz RF band
+RF_BW           =    400e6      # Hz, Sampled RF bandwidth, = f_sample/2
+CB_BW           =    300e6      # Hz, Full bandwidth of the CB and also of the TAB and IAB (SR-0.2)
+N_FFT           =   1024        # FFT size of the real input FFT in the Apertif BF subband polyphase filterbank
+T_sub           =      1.28e-6  # s, Subband period, = 1/f_sub
+f_sub           = 781250.0      # Hz Subband rate and beamlet rate in Apertif BF, f_sub = B_sub = f_sample / N_FFT
+B_sub           = 781250.0      # Hz Subband bandwidth in Apertif BF, = beamlet bandwidth
+B_max           =   3000.0      # m Maximum baseline between two dishes in the WSRT array
+f_chan_x        =  12207.03125  # Hz Channel rate in Apertif = f_sub / N_chan_x
+f_chan_a        =  24414.0625   # Hz Channel rate in Arts = f_sub / N_chan_a
+t_chan_bf       =      5.12e-6  # s, Channel period, = 1/f_chan_bf
+f_chan_bf       = 195312.5      # Hz Channel rate in Arts, f_chan_bf = f_sub / N_chan_bf = B_chan_bf
+B_chan_bf       = 195312.5      # Hz Channel bandwidth within a beamlet, for SC3 and SC4, = B_sub/N_chan_bf
+N_sub           =    512        # Number of subbands that covers RF_BW=400MHz, = N_FFT/N_complex
+M_sub           =    128        # Number subbands per sp per row in R_sub, = N_sub/P_sub
+Q_sp            =      2        # Interleave N_complex = 2 SP FFT outputs for complex FFT with two real inputs
+N_wpfb_units    =      2        # Number off WPFB units per BN, = S_BN/Q_sp
+N_sel           =    384        # Number of selected subbands to cover CB_BW=300 MHz
+N_band          =     16        # Number of bands in the Apertif BF to process the full CB_BW
+N_bf_units      =     64        # Number of bf units per polarization (is per subrack), = N_band * P_BF
+N_bb            =     24        # Number of subbands per band or per FN in the Apertif BF, = N_sel/N_band
+N_bb_per_unit   =      6        # Number of subbands per BF unit, = N_bb/P_BF
+N_bb_pol        =     48        # Number of subbands per dual polarization band, = N_bb*N_pol
+N_work          =     40        # Number of workstations in the GPU cluster of the Arts PL, = K_CB = 40
+N_CB            =     37        # Required number of compound beams
+K_CB_max        =     42        # Maximum number of beamlets per subband, = N_slot/N_bb = 1024/24 = 42.6
+K_CB            =     40        # Implemented number of beamlets per subband (>= N_CB) in W_beamlet = 6 bit mode
+K_CB_complete   =     24        # Number of complete CB with full CB_BW = 300 MHz in W_beamlet = 8 bit mode
+K_CB_fraction   =      8        # Number of fractional CB with 200 MHz in W_beamlet = 8 bit mode
+K_CB_total      =     32        # Number of 24 complete and 8 fractional CB in W_beamlet = 8 bit mode
+PN_CB           =      5        # Implemented number of CB per PN in W_beamlet = 6 bit mode, = K_CB / nof_pn
+PN_CB_complete  =      3        # Number of complete CB per PN with full CB_BW = 300 MHz in W_beamlet = 8 bit mode
+PN_CB_fraction  =      1        # Number of fractional CB per PN with 200 MHz in W_beamlet = 8 bit mode
+PN_CB_total     =      4        # Number of 3 complete and 1 fractional CB per PN in W_beamlet = 8 bit mode
+G_nof_beamlets  =     11/16.0   # Ratio of number of beamlets in W_beamlet = 8 bit mode compared to 6 bit mode
+G_load          =     44/45.0   # Ratio of data load in W_beamlet = 8 bit mode compared to 6 bit mode
+P_BF            =      4        # Number of parallel BF units per FN in the Apertif BF
+N_clk           =    256        # Number of DP clock cycles per subband period, = N_FFT/P_sub
+N_blk           = cm.sel_a_b(bm6, 240, 176)      # Number of valid DP clock cycles per subband period in the Apertif BF, <= N_clk
+N_slot          =   1024        # Maximum number of compound beamlet slots per FN output, = P_BF * N_clk
+N_beamlet       = cm.sel_a_b(bm6, 960, 704)      # Number of compound beamlet slots per FN output, maximum Nslot = 1024, actual P_BF * N_blk = 960, required N_CB * N_bb = 888
+N_pre_transpose =     16        # Number of beamlet blocks in the pre transpose for efficient DDR3 access
+Q_interleave    =      2        # Additional beamlet output interleave factor, = nof_pn/P_BF
+Q_chan          =      2        # Narrowband serialization factor for the fine channel filterbank, = Q_interleave
+Q_blk           =      1        # Number of blocks per payload, >= 1
+M_blk           = cm.sel_a_b(bm6, 120, 88)       # Default block size per subband period in the Apertif X and Arts, = N_blk/Q_interleave
+K_blk           =     64        # block size per fine channel period in the Arts BF for SC4
+K_blk           =    128        # block size per fine channel period in the Arts BF for SC3, Apertif X
+N_gr            =     12        # Required number of TAB grating lobe patterns to cover the full CB (SR-0.41)
+N_VLBI          =     12        # Required number of TABs in the central CB for VLBI, choose = N_gr (SR-0.23)
+K_TAB           =     12        # Implemented number of TABs per beamlet for SC4 (>= N_gr)
+C_TAB           =      9        # Implemented number of TABs per beamlet for SC3
+N_TAB           =    480        # Number of TABs, = K_CB*K_TAB
+N_IAB           =     40        # Number of IABs, = K_CB
+N_link          =    384        # Number of physical 10G output links of the Apertif BF, = N_PN
+N_PN            =    384        # Total number of parallel processing nodes in the Apertif BF, = Ntp * N_band
+M_PN            =    128        # Total number of parallel processing nodes in the Apertif X, Arts SC1,4 using Muni=16 UniBoards
+M_PN2           =     16        # Total number of parallel processing nodes in the Arts SC2,SC3 using Muni2=4 UniBoard2s
+N_uni           =      4        # Number UniBoards per Apertif BF subrack, so for one telescope path (= one dish, one polarization), = N_FN / nof_fn = N_BN / nof_bn
+M_uni           =     16        # Number UniBoards in the Arts BF and in Apertif X, = N_band
+M_uni2          =      4        # Number UniBoard2s in the Arts BF for commensal modes, = N_band / nof_pn2
+N_taps_sub      =     16        # Number of taps in the FIR filter section of the subband filterbank in the Apertif BF
+N_taps_chan     =      8        # Number of taps in the FIR filter section of the channel filterbank in the Apertif X
+N_chan_x        =     64        # Number of fine channels per beamlet for Arts SC3 as separated by Apertif X
+N_chan_a        =     32        # Number of fine channels per beamlet for Arts SC4
+N_chan_bf       =      4        # Required number of output channels per beamlet, for TAB and IAB beamformers in Arts SC3 and SC4
+N_vis           =    300        # Number of visibilities, N_tp * (N_tp + 1) / 2
+t_int_x         =      1.024    # s, Integration interval of the Apertif X, = N_int_x / f_sub
+t_sync          =      1.024    # = t_int_x
+N_int_x         = 800000        # Number of beamlet time samples per integration interval in the Apertif X
+N_int_chan_x    =  12500        # Number of fine channel time samples per integration interval in the Arts SC3
+N_int_chan_a    =  25000        # Number of fine channel time samples per integration interval in the Arts SC4
+M_int_a         =      8        # Number of fine channel IQUV power values that are integrated in Arts SC4
+M_int_ax        =     16        # Number of fine channel IQUV power values that are integrated in Arts SC3
+T_chan_x        =     81.92e-6  # s, Period of a fine channel for SC3, = N_chan_x * T_sub
+T_chan_a        =     40.96e-6  # s, Period of a fine channel for SC4, = N_chan_a * T_sub
+T_Stokes        =     50e-6     # s, Minimum required sample period for the Stokes power values
+f_Stokes        =     20e3      # Hz, Minimum required sample frequency for the Stokes power values, = 1/T_Stokes
+nof_bn          =      4        # Number of back node FPGAs (BN) per UniBoard
+nof_fn          =      4        # Number of front node FPGAs (FN) per UniBoard
+nof_pn          =      8        # Number of processing node FPGAs per UniBoard, = nof_fn + nof_bn
+nof_pn2         =      4        # number of processing node FPGAs per UniBoard2
+nof_10g         =      3        # Number of external 10G links per FPGA node on UniBoard
+nof_5g_mesh     =      3        # Number of 5G links between a pair of FN and BN via the on board UniBoard mesh
+nof_10g2        =     24        # Number of external 10G links per FPGA node on UniBoard2
+nof_10g2_ring   =     12        # Number of 10G links between PN2 via the on board UniBoard2 ring
+byte_w          =      8        # Number of bits in a byte or an octet
+word_sz         =      4        # Number of bytes per 32 bit long word
+longword_sz     =      8        # Number of bytes per 64 bit long word
+bram_sz         =   1024        # Number of bytes in a BRAM, 1024 for M9K in Stratix IV on UniBoard
+bram2_sz        =   2048        # Number of bytes in a BRAM, 2048 for M20K in Arria10 on UniBoard2
+
+W_adc           =      8        # Word width in number of bits of the ADC
+W_adc_max       = W_adc-1       # = 7 = 8 - 1
+W_fsub_gain     = math.log(math.sqrt(N_FFT), 2)    # = 5, incoherent noise processing gain of F_sub in number of bits
+W_fchan_x_gain  = math.log(math.sqrt(N_chan_x), 2) # = 3, incoherent noise processing gain of F_chan_x in number of bits
+W_sst_fraction  =      8        # Extra fraction bits for SST subband input compared to W_adc, = W_fsub_gain + 3
+W_sst_in        =     16        # Word width in number of bits of SST subband data input, = W_adc + W_sst_fraction
+W_sst_out       =     56        # Word width in number of bits of SST subband statistics output
+W_sst_out_sz    =      2        # Number of 32b words per W_sst_out
+
+W_cb_in         =     16        # Word width in number of bits of BF subband data input
+W_cb_weight     =     16        # Word width in number of bits of the Apertif 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)
+
+W_FIR_sub       =     16        # Word width in number of bits of a subband filterbank FIR filter coefficient
+W_FIR_chan      =      9        # Word width in number of bits of a channel filterbank FIR filter coefficient
+W_vis           =     32        # Word width in number of bits of the visibility statistics
+W_volt          =      8        # Word width in number of bits of an Arts BF output voltage sample
+W_volt_max      = W_volt-1      # 7 = 8 - 1
+W_pow           =      8        # Word width in number of bits of an Arts BF output power sample
+W_pow_max       = W_pow-1       # 7 = 8 - 1
+P_rcmult        =      2        # Number of real multiplications per real * complex multiplication (9 bit)
+P_cmult         =      4        # Number of real multiplications per complex * complex multiplication
+W_rcmult        =      9        # bit, Maximum operand width for a hard IP real * complex multiplier
+W_cmult         =     18        # bit, Maximum operand width for a hard IP complex * complex multiplier
+W_DSP           =     18        # Default data width, matches data width of multipliers and block RAM in the FPGA
+
+nof_samples_per_sync = f_sample * t_sync       # = 800e6 * 1.024 = 819200000
+
+MAC_beamlets_scheme = '300MHz_40CB'    # Default R_beam beamlet mapping used by Apertif MAC
+