1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
|
--[[
This is a strategy for ATCA wideband L-band, data version 2021-03-30
Author: André Offringa
It's based on 3 observations provided by Björn Adebahr and Ancla Müller. It
works well on both the target as calibrator observations. Baseline CA4 x CA5
has so much RFI that it probably should be flagged in full -- This strategy
flags 70% of that baseline, but the other 30% seems also bad.
Changes compared to the generic strategy:
- Don't flag on Stokes I, only on Q, U and V. The bandpass structure has too many
features.
- Base theshold is set to 1.2. This can be further tweaked to control the level of flagged
data vs. the level of residual RFI.
- frequency resize factor is set to 3, which makes the background fit somewhat more stable.
The large bandwidth vs. the "resolved" RFI features make this necessary.
]]
aoflagger.require_min_version("3.0")
function execute(input)
--
-- Generic settings
--
-- What polarizations to flag? Default: input:get_polarizations() (=all that are in the input data)
-- Other options are e.g.:
-- { 'XY', 'YX' } to flag only XY and YX, or
-- { 'I', 'Q' } to flag only on Stokes I and Q
local flag_polarizations
if #input:get_polarizations() == 4 then
flag_polarizations = { "Q", "U", "V" }
else
flag_polarizations = input:get_polarizations()
end
local base_threshold = 1.6 -- lower means more sensitive detection
-- How to flag complex values, options are: phase, amplitude, real, imaginary, complex
-- May have multiple values to perform detection multiple times
local flag_representations = { "amplitude" }
local iteration_count = 6 -- how many iterations to perform?
local threshold_factor_step = 1.5 -- How much to increase the sensitivity each iteration?
-- If the following variable is true, the strategy will consider existing flags
-- as bad data. It will exclude flagged data from detection, and make sure that any existing
-- flags on input will be flagged on output. If set to false, existing flags are ignored.
local use_input_flags = true
local frequency_resize_factor = 3.0 -- Amount of "extra" smoothing in frequency direction
local transient_threshold_factor = 1.0 -- decreasing this value makes detection of transient RFI more aggressive
--
-- End of generic settings
--
local inpPolarizations = input:get_polarizations()
if not use_input_flags then
input:clear_mask()
end
-- For collecting statistics. Note that this is done after clear_mask(),
-- so that the statistics ignore any flags in the input data.
local copy_of_input = input:copy()
for ipol, polarization in ipairs(flag_polarizations) do
local pol_data = input:convert_to_polarization(polarization)
local original_data
for _, representation in ipairs(flag_representations) do
data = pol_data:convert_to_complex(representation)
original_data = data:copy()
for i = 1, iteration_count - 1 do
local threshold_factor = threshold_factor_step ^ (iteration_count - i)
local sumthr_level = threshold_factor * base_threshold
if use_input_flags then
aoflagger.sumthreshold_masked(
data,
original_data,
sumthr_level,
sumthr_level * transient_threshold_factor,
true,
true
)
else
aoflagger.sumthreshold(data, sumthr_level, sumthr_level * transient_threshold_factor, true, true)
end
-- Do timestep & channel flagging
local chdata = data:copy()
aoflagger.threshold_timestep_rms(data, 3.5)
aoflagger.threshold_channel_rms(chdata, 3.0 * threshold_factor, true)
data:join_mask(chdata)
-- High pass filtering steps
data:set_visibilities(original_data)
if use_input_flags then
data:join_mask(original_data)
end
local resized_data = aoflagger.downsample(data, 3, frequency_resize_factor, true)
aoflagger.low_pass_filter(resized_data, 51, 51, 2.5, 5.0)
aoflagger.upsample(resized_data, data, 3, frequency_resize_factor)
-- In case this script is run from inside rfigui, calling
-- the following visualize function will add the current result
-- to the list of displayable visualizations.
-- If the script is not running inside rfigui, the call is ignored.
aoflagger.visualize(data, "Fit #" .. i, i - 1)
local tmp = original_data - data
tmp:set_mask(data)
data = tmp
aoflagger.visualize(data, "Residual #" .. i, i + iteration_count)
aoflagger.set_progress((ipol - 1) * iteration_count + i, #flag_polarizations * iteration_count)
end -- end of iterations
if use_input_flags then
aoflagger.sumthreshold_masked(
data,
original_data,
base_threshold,
base_threshold * transient_threshold_factor,
true,
true
)
else
aoflagger.sumthreshold(data, base_threshold, base_threshold * transient_threshold_factor, true, true)
end
end -- end of complex representation iteration
if use_input_flags then
data:join_mask(original_data)
end
-- Helper function used below
function contains(arr, val)
for _, v in ipairs(arr) do
if v == val then
return true
end
end
return false
end
if contains(inpPolarizations, polarization) then
if input:is_complex() then
data = data:convert_to_complex("complex")
end
input:set_polarization_data(polarization, data)
else
input:join_mask(data)
end
aoflagger.visualize(data, "Residual #" .. iteration_count, 2 * iteration_count)
aoflagger.set_progress(ipol, #flag_polarizations)
end -- end of polarization iterations
if use_input_flags then
aoflagger.scale_invariant_rank_operator_masked(input, copy_of_input, 0.2, 0.2)
else
aoflagger.scale_invariant_rank_operator(input, 0.2, 0.2)
end
aoflagger.threshold_timestep_rms(input, 4.0)
if input:is_complex() and input:has_metadata() then
-- This command will calculate a few statistics like flag% and stddev over
-- time, frequency and baseline and write those to the MS. These can be
-- visualized with aoqplot.
aoflagger.collect_statistics(input, copy_of_input)
end
end
|
