Cheetah - SKA - PSS - Prototype Time Domain Search Pipeline
RfiExcision.cpp
1 /*
2  * The MIT License (MIT)
3  *
4  * Copyright (c) 2016 The SKA organisation
5  *
6  * Permission is hereby granted, free of charge, to any person obtaining a copy
7  * of this software and associated documentation files (the "Software"), to deal
8  * in the Software without restriction, including without limitation the rights
9  * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
10  * copies of the Software, and to permit persons to whom the Software is
11  * furnished to do so, subject to the following conditions:
12  *
13  * The above copyright notice and this permission notice shall be included in all
14  * copies or substantial portions of the Software.
15  *
16  * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
17  * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
18  * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
19  * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
20  * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
21  * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
22  * SOFTWARE.
23  */
24 #include "cheetah/sps/detail/RfiExcision.h"
25 #include <algorithm>
26 
27 namespace ska {
28 namespace cheetah {
29 namespace sps {
30 
31 template<typename TimeFrequencyType, typename RfiFlagDataType>
32 RfiExcision<TimeFrequencyType, RfiFlagDataType>::RfiExcision( RfiFlagDataType& data, size_t offset, float ideal_rms, float threshold, bool active)
33  : _flag_it(data.rfi_flags().begin() + offset)
34  , _tf_it(data.tf_data().begin() + offset)
35  , _wrapindex(0)
36  , _ideal_rms(ideal_rms)
37  , _threshold(threshold)
38  , _active(active)
39  , _tf_end_it(data.tf_data().cend())
40 {
41  assert(data.rfi_flags().data_size() == data.tf_data().data_size());
42  _noise.resize(8192); // A vector to save pre calculated random numbers.
43  // The size is fixed such that it s not too big to waste resources or too small to reduce the 'randomness' of the replacement values
44 
45  // we take the median of the means for each channel as our default replacement value
46  typedef typename RfiFlagDataType::Statistics Statistics;
47  auto channel_stats = data.channel_stats();
48  auto subband_nchans = channel_stats.size() / 8; // Dividing the band into 8 subbands so as to choose the best median from the middle of the band
49  std::nth_element(channel_stats.begin() + 2*subband_nchans, channel_stats.begin() + (std::size_t)(2.5)*subband_nchans, channel_stats.begin() + 3*subband_nchans - 1,
50  [](Statistics const& s1, Statistics const& s2)
51  {
52  return s1.mean < s2.mean;
53  } );
54  _replacement_value = channel_stats[(std::size_t)((2.5)*subband_nchans)].mean;
55 
56  /* Setting up the random generator */
57  static std::random_device rd;
58  std::seed_seq seed { rd(), rd(), rd(), rd(), rd() };
59  auto e2 = std::mt19937( seed );
60  auto distribution = std::normal_distribution<float>(_replacement_value,_ideal_rms);
61  std::generate(_noise.begin(),_noise.end(), std::bind(distribution, e2));
62  _value = _noise[0];
63 }
64 
65 template<typename TimeFrequencyType, typename RfiFlagDataType>
66 RfiExcision<TimeFrequencyType, RfiFlagDataType>::~RfiExcision()
67 {
68 }
69 
70 template<typename TimeFrequencyType, typename RfiFlagDataType>
71 bool RfiExcision<TimeFrequencyType, RfiFlagDataType>::operator!=( RfiExcision<TimeFrequencyType, RfiFlagDataType> const& iterator) const
72 {
73  return this->_tf_it != iterator._tf_it;
74 }
75 
76 template<typename TimeFrequencyType, typename RfiFlagDataType>
77 bool RfiExcision<TimeFrequencyType, RfiFlagDataType>::operator==( RfiExcision<TimeFrequencyType, RfiFlagDataType> const& iterator) const
78 {
79  return this->_tf_it == iterator._tf_it;
80 }
81 
82 template<typename TimeFrequencyType, typename RfiFlagDataType>
83 RfiExcision<TimeFrequencyType, RfiFlagDataType>& RfiExcision<TimeFrequencyType, RfiFlagDataType>::operator++()
84 {
85  ++_tf_it;
86  ++_flag_it;
87  ++_wrapindex;
88  if(_tf_it!=_tf_end_it)
89  {
90  bool do_replacement;
91  if (_active)
92  do_replacement = (*(_flag_it)) || (std::fabs((float) *(_tf_it) - _replacement_value) >= _threshold *_ideal_rms);
93  else
94  do_replacement = (*(_flag_it));
95 
96  if (do_replacement)
97  {
98  if (_wrapindex >= _noise.size())
99  {
100  _wrapindex = 0;
101  }
102  _value = static_cast<NumericalRep>(_noise[_wrapindex]);
103  }
104  else
105  _value = *(_tf_it);
106 
107  }
108  return *this;
109 }
110 
111 template<typename TimeFrequencyType, typename RfiFlagDataType>
112 typename RfiExcision<TimeFrequencyType, RfiFlagDataType>::const_reference RfiExcision<TimeFrequencyType, RfiFlagDataType>::operator*() const
113 {
114  return _value;
115 }
116 
117 } // namespace rfim
118 } // namespace cheetah
119 } // namespace ska
Some limits and constants for FLDO.
Definition: Brdz.h:35