Skip to content

Commit 99c006a

Browse files
timcornwelltimcornwell
committed
Add RFI simulation
1 parent 86dde15 commit 99c006a

2 files changed

Lines changed: 260 additions & 0 deletions

File tree

‎processing_components/simulation/rfi.py‎

Lines changed: 180 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,180 @@
1+
import numpy
2+
from astropy import constants
3+
import astropy.units as u
4+
from astropy.coordinates import SkyCoord
5+
6+
from processing_library.util.array_functions import average_chunks2
7+
from processing_library.util.coordinate_support import xyz_to_uvw, skycoord_to_lmn
8+
9+
def simulate_DTV(frequency, times, power=50e3, timevariable=False, frequency_variable=False):
10+
""" Calculate DTV sqrt(power) as a function of time and frequency
11+
12+
:param frequency: (sample frequencies)
13+
:param times: sample times (s)
14+
:param power: DTV emitted power W
15+
:return: Complex array [ntimes, nchan]
16+
"""
17+
nchan = len(frequency)
18+
ntimes = len(times)
19+
shape = [ntimes, nchan]
20+
bchan = nchan // 4
21+
echan = 3 * nchan // 4
22+
amp = power / (max(frequency) - min(frequency))
23+
signal = numpy.zeros(shape, dtype='complex')
24+
if timevariable:
25+
if frequency_variable:
26+
sshape = [ntimes, nchan // 2]
27+
signal[:, bchan:echan] += numpy.random.normal(0.0, numpy.sqrt(amp/2.0), sshape) \
28+
+ 1j * numpy.random.normal(0.0, numpy.sqrt(amp/2.0), sshape)
29+
else:
30+
sshape = [ntimes]
31+
signal[:, bchan:echan] += numpy.random.normal(0.0, numpy.sqrt(amp/2.0), sshape) \
32+
+ 1j * numpy.random.normal(0.0, numpy.sqrt(amp/2.0), sshape)
33+
else:
34+
if frequency_variable:
35+
sshape = [nchan // 2]
36+
signal[:, bchan:echan] += (numpy.random.normal(0.0, numpy.sqrt(amp/2.0), sshape)
37+
+ 1j * numpy.random.normal(0.0, numpy.sqrt(amp/2.0), sshape))[numpy.newaxis, ...]
38+
else:
39+
signal[:, bchan:echan] = amp
40+
41+
return signal
42+
43+
44+
def add_noise(visibility, bandwidth, int_time):
45+
"""Determine noise rms per visibility
46+
47+
:returns: visibility with noise added
48+
"""
49+
# The specified sensitivity (effective area / T_sys) is roughly 610 m ^ 2 / K in the range 160 - 200MHz
50+
# sigma_vis = 2 k T_sys / (area * sqrt(tb)) = 2 k 512 / (610 * sqrt(tb)
51+
sens = 610
52+
k_b = 1.38064852e-23
53+
bt = bandwidth * int_time
54+
sigma = 2 * 1e26 * k_b / ((sens / 512) * (numpy.sqrt(bt)))
55+
print("RMS noise per sample = %g Jy" % sigma)
56+
sshape = visibility.shape
57+
visibility += numpy.random.normal(0.0, sigma, sshape) + 1j * numpy.random.normal(0.0, sigma, sshape)
58+
return visibility
59+
60+
61+
def create_propagators(config, interferer, frequency, attenuation=1e-9):
62+
""" Create a set of propagators
63+
64+
:return: Complex array [nants, ntimes]
65+
"""
66+
nchannels = len(frequency)
67+
nants = len(config.data['names'])
68+
interferer_xyz = [interferer.geocentric[0].value, interferer.geocentric[1].value, interferer.geocentric[2].value]
69+
propagators = numpy.zeros([nants, nchannels], dtype='complex')
70+
for iant, ant_xyz in enumerate(config.xyz):
71+
vec = ant_xyz - interferer_xyz
72+
# This ignores the Earth!
73+
r = numpy.sqrt(vec[0] ** 2 + vec[1] ** 2 + vec[2] ** 2)
74+
k = 2.0 * numpy.pi * frequency / constants.c.value
75+
propagators[iant, :] = numpy.exp(- 1.0j * k * r) / r
76+
return propagators * attenuation
77+
78+
79+
def calculate_rfi_at_station(propagators, emitter):
80+
""" Calculate the rfi at each station
81+
82+
:param propagators:
83+
:param emitter:
84+
:return: Complex array [nants, ntimes, nchannels]
85+
"""
86+
rfi_at_station = emitter[numpy.newaxis, :, :] * propagators[:, numpy.newaxis, :]
87+
rfi_at_station[numpy.abs(rfi_at_station)<1e-15] = 0.
88+
return rfi_at_station
89+
90+
91+
def calculate_station_fringe_rotation(ants_xyz, times, frequency, phasecentre, pole):
92+
"""Calculate fringe rotation vectors
93+
94+
:param ants_xyz:
95+
:param times:
96+
:param frequency:
97+
:param phasecentre:
98+
:param pole:
99+
:return: phasor, uvw
100+
"""
101+
102+
uvw = xyz_to_uvw(ants_xyz, times, phasecentre.dec.rad)
103+
nants, nuvw = uvw.shape
104+
ntimes = len(times)
105+
uvw = uvw.reshape([nants, ntimes, 3])
106+
lmn = skycoord_to_lmn(phasecentre, pole)
107+
delay = numpy.dot(uvw, lmn)
108+
nchan = len(frequency)
109+
phase = numpy.zeros([nants, ntimes, nchan])
110+
for ant in range(nants):
111+
for chan in range(nchan):
112+
phase[ant, :, chan] = delay[ant] * frequency[chan] / constants.c.value
113+
phase[...] -= phase[0, :, :][numpy.newaxis,...]
114+
return numpy.exp(2.0 * numpy.pi * 1j * phase), uvw
115+
116+
117+
def calculate_station_correlation_rfi(fringe_rotation, rfi_at_station):
118+
""" Form the correlation from the rfi at the station and the fringe rotation
119+
120+
:param fringe_rotation:
121+
:param rfi_at_station:
122+
:return: Correlation(nant, nants, ntimes, nchan] in Jy
123+
"""
124+
phased_rotated_rfi_at_station = fringe_rotation * rfi_at_station
125+
nants, ntimes, nchan = fringe_rotation.shape
126+
correlation = numpy.zeros([nants, nants, ntimes, nchan], dtype='complex')
127+
128+
for time in range(ntimes):
129+
for chan in range(nchan):
130+
correlation[..., time, chan] = numpy.outer(phased_rotated_rfi_at_station[..., time, chan],
131+
numpy.conjugate(phased_rotated_rfi_at_station[..., time, chan]))
132+
return correlation * 1e26
133+
134+
135+
def calculate_averaged_correlation(correlation, channel_width, time_width):
136+
""" Average the correlation in time and frequency
137+
138+
:param correlation: Correlation(nant, nants, ntimes, nchan]
139+
:param channel_width: Number of channels to average
140+
:param time_width: Number of integrations to average
141+
:return:
142+
"""
143+
wts = numpy.ones(correlation.shape, dtype='float')
144+
return average_chunks2(correlation, wts, (channel_width, time_width))[0]
145+
146+
147+
def simulate_rfi_block(config, times, frequency, phasecentre, emitter_location, emitter_power=5e4,
148+
attenuation=1.0):
149+
""" Simulate RFI block
150+
151+
:param config: ARL telescope Configuration
152+
:param times: observation times (hour angles)
153+
:param frequency: frequencies
154+
:param phasecentre:
155+
:param emitter_location: EarthLocation of emitter
156+
:param emitter_power: Power of emitter
157+
:param attenuation: Attenuation to be applied to signal
158+
:return:
159+
"""
160+
# Calculate the power spectral density of the DTV station: Watts/Hz
161+
emitter = simulate_DTV(frequency, times, power=emitter_power, timevariable=False)
162+
163+
# Calculate the propagators for signals from Perth to the stations in low
164+
# These are fixed in time but vary with frequency. The ad hoc attenuation
165+
# is set to produce signal roughly equal to noise at LOW
166+
propagators = create_propagators(config, emitter_location, frequency=frequency,
167+
attenuation=attenuation)
168+
# Now calculate the RFI at the stations, based on the emitter and the propagators
169+
rfi_at_station = calculate_rfi_at_station(propagators, emitter)
170+
171+
# Station fringe rotation: shape [nants, ntimes, nchan] complex phasor to be applied to
172+
# reference to the pole.
173+
pole = SkyCoord(ra=+0.0 * u.deg, dec=-90.0 * u.deg, frame='icrs', equinox='J2000')
174+
fringe_rotation, uvw = calculate_station_fringe_rotation(config.xyz, times, frequency, phasecentre, pole)
175+
176+
# Calculate the rfi correlationm using the fringe rotation and the rfi at the station
177+
# [nants, nants, ntimes, nchan]
178+
correlation = calculate_station_correlation_rfi(fringe_rotation, rfi_at_station)
179+
180+
return correlation, uvw

‎tests/processing_components/test_rfi.py‎

Lines changed: 80 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,80 @@
1+
""" Unit tests for RFI simulation
2+
3+
"""
4+
5+
import logging
6+
import unittest
7+
8+
import astropy.units as u
9+
import numpy
10+
import numpy.testing
11+
12+
from astropy.coordinates import SkyCoord, EarthLocation
13+
14+
from data_models.polarisation import PolarisationFrame
15+
from processing_components.simulation.configurations import create_named_configuration
16+
from processing_components.simulation.rfi import create_propagators, calculate_averaged_correlation, \
17+
calculate_rfi_at_station, calculate_station_correlation_rfi, calculate_station_fringe_rotation, \
18+
simulate_DTV, add_noise
19+
from processing_components.simulation.testing_support import create_test_image
20+
from processing_components.skycomponent.operations import create_skycomponent
21+
from processing_components.visibility.base import create_blockvisibility
22+
from processing_library.image.operations import create_image
23+
24+
log = logging.getLogger(__name__)
25+
26+
27+
class TestRFISim(unittest.TestCase):
28+
def setUp(self):
29+
pass
30+
31+
def test_rfi(self):
32+
33+
sample_freq = 3e4
34+
nchannels = 1000
35+
frequency = 170.5e6 + numpy.arange(nchannels) * sample_freq
36+
37+
ntimes = 100
38+
integration_time = 0.5
39+
times = numpy.arange(ntimes) * integration_time
40+
41+
phasecentre = SkyCoord(ra=+30.0 * u.deg, dec=-45.0 * u.deg, frame='icrs', equinox='J2000')
42+
pole = SkyCoord(ra=+0.0 * u.deg, dec=-90.0 * u.deg, frame='icrs', equinox='J2000')
43+
44+
# Perth from Google for the moment
45+
perth = EarthLocation(lon="115.8605", lat="-31.9505", height=0.0)
46+
47+
rmax = 1000.0
48+
low = create_named_configuration('LOWR3', rmax=rmax)
49+
antskip = 33
50+
low.data = low.data[::antskip]
51+
nants = len(low.names)
52+
53+
# Calculate the power spectral density of the DTV station: Watts/Hz
54+
emitter = simulate_DTV(frequency, times, power=50e3, timevariable=False)
55+
numpy.testing.assert_almost_equal(numpy.max(numpy.abs(emitter)), 0.00166834)
56+
assert emitter.shape == (ntimes, nchannels)
57+
58+
# Calculate the propagators for signals from Perth to the stations in low
59+
# These are fixed in time but vary with frequency. The ad hoc attenuation
60+
# is set to produce signal roughly equal to noise at LOW
61+
attenuation = 1.0
62+
propagators = create_propagators(low, perth, frequency=frequency, attenuation=attenuation)
63+
assert propagators.shape == (nants, nchannels), propagators.shape
64+
65+
# Now calculate the RFI at the stations, based on the emitter and the propagators
66+
rfi_at_station = calculate_rfi_at_station(propagators, emitter)
67+
assert rfi_at_station.shape == (nants, ntimes, nchannels), rfi_at_station.shape
68+
69+
# Station fringe rotation: shape [nants, ntimes, nchan] complex phasor to be applied to
70+
# reference to the pole.
71+
fringe_rotation, uvw = calculate_station_fringe_rotation(low.xyz, times, frequency, phasecentre, pole)
72+
numpy.testing.assert_allclose(numpy.abs(fringe_rotation), 1.0)
73+
74+
# Calculate the rfi correlationm using the fringe rotation and the rfi at the station
75+
# [nants, nants, ntimes, nchan]
76+
correlation = calculate_station_correlation_rfi(fringe_rotation, rfi_at_station)
77+
assert correlation.shape == (nants, nants, ntimes, nchannels), correlation.shape
78+
79+
if __name__ == '__main__':
80+
unittest.main()

0 commit comments

Comments
 (0)