1 #include "cheetah/rfim/sum_threshold/Rfim.h" 2 #include "panda/Buffer.h" 4 #include "panda/TypeTraits.h" 5 #pragma GCC diagnostic push 6 #pragma GCC diagnostic ignored "-Wunused-parameter" 7 #include <boost/accumulators/accumulators.hpp> 8 #include <boost/accumulators/statistics.hpp> 9 #pragma GCC diagnostic pop 16 namespace sum_threshold {
19 typedef data::TimeFrequencyFlags<Cpu> FlagType;
21 inline float calcThresholdI(
float threshold1,
unsigned window,
float p)
26 return (
float) (threshold1 * std::pow(p, std::log2(window)) / window);
29 template<
typename IteratorType>
30 static void calculateNormalStatistics(
31 IteratorType&& begin, IteratorType
const&& end,
float& mean,
float& median,
float& std_dev)
33 boost::accumulators::accumulator_set<
34 double, boost::accumulators::stats<boost::accumulators::tag::variance,
35 boost::accumulators::tag::median>> acc;
36 for_each(begin, end, std::bind<void>(std::ref(acc), std::placeholders::_1));
38 mean = boost::accumulators::mean(acc);
39 std_dev = sqrtf(boost::accumulators::variance(acc));
40 median = boost::accumulators::median(acc);
43 template<
typename IteratorType>
44 static void calculateStatistics(IteratorType&& begin, IteratorType
const&& end,
float& mean,
float& median,
float& std_dev) {
46 calculateNormalStatistics(std::forward<IteratorType>(begin), std::forward<IteratorType const>(end), mean, median, std_dev);
51 template<
typename DataType>
52 static unsigned sumThreshold2DHorizontal(DataType
const & powers, FlagType &flags,
const unsigned window,
const float threshold) {
53 unsigned extra_flagged = 0;
55 for(data::DimensionIndex<data::Frequency> channel(0); channel<powers.number_of_channels(); ++channel) {
56 for (
unsigned base = 0; base + window < powers.number_of_spectra(); ++base) {
58 for(data::DimensionIndex<data::Time> time(base); time < base + window; ++time) {
59 if (flags[channel][time]) {
62 sum += powers.spectrum(time)[channel];
66 if (sum >= window * threshold) {
68 for(data::DimensionIndex<data::Time> time(base); time < base + window; ++time) {
69 if(!flags[channel][time]) {
71 flags[channel][time] =
true;
81 template<
typename DataType>
82 static unsigned sumThreshold2DVertical(DataType
const & powers
84 ,
const unsigned window
85 ,
const float threshold)
87 unsigned extra_flagged = 0;
89 for(data::DimensionIndex<data::Frequency> channel(0); channel<powers.number_of_channels(); ++channel) {
90 for (
unsigned base = 0; base + window < powers.number_of_spectra(); base++) {
93 for (data::DimensionIndex<data::Time> time(base); time < base + window; ++time) {
94 if(flags[channel][time]) {
97 sum += powers.spectrum(time)[channel];
101 if (sum >= window * threshold) {
103 for (data::DimensionIndex<data::Time> time(base); time < base + window; ++time) {
104 if(!flags[channel][time]) {
106 flags[channel][time] =
true;
113 return extra_flagged;
117 template<
typename RfimTraits>
118 Rfim<RfimTraits>::Rfim(
const Config& config)
123 template<
typename RfimTraits>
124 Rfim<RfimTraits>::~Rfim()
131 template<
typename DataAdapterType,
typename Enable=
void>
132 struct FlagTypeHelper {
133 typedef FlagType ReturnType;
134 template<
typename DataType>
135 FlagType flags(DataAdapterType&, DataType
const& data) {
136 return FlagType(data.template dimension<data::Time>(), data.template dimension<data::Frequency>());
141 template<
typename DataAdapterType>
142 struct FlagTypeHelper<DataAdapterType, typename
std::enable_if<std::is_convertible<typename panda::is_pointer_wrapper<typename DataAdapterType::ReturnType>::type, FlagType&>::value>::type> {
143 typedef FlagType& ReturnType;
144 template<
typename DataType>
145 FlagType& flags(DataAdapterType& adapter, DataType
const&) {
146 return static_cast<FlagType&
>(*adapter.data());
152 template<
typename RfimTraits>
153 template<
typename DataType>
154 void Rfim<RfimTraits>::operator()(DataType& data, DataAdapter& adapter)
156 float mean, std_dev, median;
157 calculateStatistics(data.begin(), data.end(), mean, median, std_dev);
159 if (std_dev == 0.0f) {
160 factor = _config.base_sensitivity();
162 factor = std_dev * _config.base_sensitivity();
165 unsigned extra_flagged = 0;
167 FlagTypeHelper<DataAdapter> flag_helper;
168 typename FlagTypeHelper<DataAdapter>::ReturnType flags = flag_helper.flags(adapter, data);
169 for (
unsigned window : _config.thresholding_data_sizes()) {
170 float thresholdI = median + calcThresholdI(_config.its_cutoff_threshold(), window, 1.5f) * factor;
171 extra_flagged += sumThreshold2DHorizontal(data, flags, window, thresholdI);
172 extra_flagged += sumThreshold2DVertical(data, flags, window, thresholdI);
177 auto it=data.cbegin();
178 auto flag_it =flags.cbegin();
179 for ( data::DimensionIndex<data::Time> spectrum_num(0); spectrum_num < data.number_of_spectra(); ++spectrum_num ) {
180 std::size_t bad_channel=0;
181 for ( data::DimensionIndex<data::Frequency> channel_num(0); channel_num < data.number_of_channels(); ++channel_num ) {
183 adapter.mark_bad(*it, channel_num);
190 if(bad_channel==0) adapter.mark_good(data.spectrum(spectrum_num));
Some limits and constants for FLDO.