"""
Script to add a static RFI signal to the measurement set.
Input: HDF file of signal (currently using a modified version for LOW signals)
"""
import h5py
import numpy
from numpy.random import default_rng
from rascil.processing_components import simulate_rfi_block_prop
from ska_sdp_datamodels.visibility import (
create_visibility_from_ms,
export_visibility_to_ms,
)
# add correct path
SOURCE = "point" # "double" or "s3sky"
INPUT_MS = "rcal_" + SOURCE + ".ms"
OUTPUT_MS = "aa05_mid_rfi_84chans_" + SOURCE + ".ms"
RFI_HDF = "tv_transmitter_attenuation_cube_aa05_200chans.hdf5"
[docs]
def gen_rfi_for_ms():
"Generate RFI for given MS"
rng = default_rng(1805550721)
bvis_list = create_visibility_from_ms(INPUT_MS)
for bvis in bvis_list:
bvis.configuration.attrs["name"] = "MID"
# We only use the first 120 time samples and
# 4 stations for MID-AA0.5 (what input_ms contains)
# While the input RFI array contains 240 time samples and 6 antennas
with h5py.File(RFI_HDF, "r") as hdf:
coords = hdf["coordinates"][()][:, :120, :4, :]
rfi_freq_channels = hdf["frequency_channels"][()]
source_ids = hdf["source_id"][()]
rfi_signal = hdf["signal"][()][:, :120, :4, :]
source_ids = [sid.decode() for sid in source_ids]
# center the 87 channels that have signal. Originally those are the last 87 channels
# pylint: disable-next=assignment-from-no-return,no-value-for-parameter
new_rfi_signal = rfi_signal.copy()
# pylint: disable-next=unsupported-assignment-operation
new_rfi_signal[:, :, :, : (200 - 87) // 2] = (
0.0 # pylint: disable=unsupported-assignment-operation
)
ind = numpy.where(rfi_signal != 0.0)
new_rfi_signal[..., (200 - 87) // 2 : (200 - 87) // 2 + 87] = rfi_signal[
..., numpy.unique(ind[3])
]
new_rfi_signal[..., (200 - 87) // 2 + 87 :] = 0.0
# Change to MID's frequency sampling : 365 - 415 MHz
rfi_freq_channels = numpy.array(rfi_freq_channels) * 250 / 375
rfi_freq_channels = rfi_freq_channels + 2.9e8
noise_bvis_list = []
for bvis in bvis_list:
n_times = bvis.dims["time"]
new_bvis = bvis.copy(deep=True, zero=True)
new_bvis = simulate_rfi_block_prop(
new_bvis,
new_rfi_signal,
coords,
source_ids,
rfi_freq_channels,
low_beam_gain=None,
apply_primary_beam=False, # not needed for Low
)
# remove signal from half of the time samples randomly chosen
rand_times = rng.choice(range(n_times), n_times // 2, replace=False)
rand_times.sort()
new_bvis["vis"].data[rand_times, ...] = 0.0
noise_bvis_list.append(new_bvis)
# ended up with 84 channels with signal
export_visibility_to_ms(OUTPUT_MS, noise_bvis_list)
if __name__ == "__main__":
gen_rfi_for_ms()