diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000000000000000000000000000000000000..ce5c1499ba06e98ca35ddbb143bd180b8ed1fa22 --- /dev/null +++ b/.gitignore @@ -0,0 +1,4 @@ +*log +__pycache__ +*pyc +*# \ No newline at end of file diff --git a/LICENSE.md b/LICENSE.md new file mode 100644 index 0000000000000000000000000000000000000000..5dc22d844f4c601d2fbe5b083f764ade05af57bc --- /dev/null +++ b/LICENSE.md @@ -0,0 +1,32 @@ +# 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 + + + diff --git a/NOTICE.md b/NOTICE.md new file mode 100755 index 0000000000000000000000000000000000000000..0263080bdfbe827c21959b9346f7c0b875c52976 --- /dev/null +++ b/NOTICE.md @@ -0,0 +1,15 @@ +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. +``` + diff --git a/README.md b/README.md new file mode 100644 index 0000000000000000000000000000000000000000..d9b201893dc6a023ce86b14072a2908a0cc12622 --- /dev/null +++ b/README.md @@ -0,0 +1,31 @@ +# 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 diff --git a/matlab/ASTRON_RP_010888_Apertif_Arts_firmware_model.pdf b/matlab/ASTRON_RP_010888_Apertif_Arts_firmware_model.pdf new file mode 100755 index 0000000000000000000000000000000000000000..799ce70e26e928a9622cb46e54f0c070f213a0e4 Binary files /dev/null and b/matlab/ASTRON_RP_010888_Apertif_Arts_firmware_model.pdf differ diff --git a/matlab/apertif_arts_firmware_model.py b/matlab/apertif_arts_firmware_model.py index 1aae4143d01720dc1d28d0d342bbaca81f7c98ef..706cd1d86f6ec60f68644e4f87b795b0a31e09dc 100644 --- a/matlab/apertif_arts_firmware_model.py +++ b/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 diff --git a/matlab/pi_apertif_system_constants.py b/matlab/pi_apertif_system_constants.py new file mode 100755 index 0000000000000000000000000000000000000000..cc34a8fb1d9de34ce21c289d2a0ca8cca464cd3f --- /dev/null +++ b/matlab/pi_apertif_system_constants.py @@ -0,0 +1,240 @@ +############################################################################### +# +# 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 + diff --git a/readme.md b/readme.md deleted file mode 100644 index e69de29bb2d1d6434b8b29ae775ad8c2e48c5391..0000000000000000000000000000000000000000