RdChannelMedianFilter.cc

Go to the documentation of this file.

#include "RdChannelMedianFilter.h"

#include <fwk/CentralConfig.h>

#include <utl/ErrorLogger.h>
#include <utl/Reader.h>
#include <utl/config.h>
#include <utl/TraceAlgorithm.h>

#include <evt/Event.h>
#include <revt/REvent.h>
#include <revt/Station.h>
#include <revt/Channel.h>


using namespace fwk;
using namespace utl;
using namespace std;
using namespace revt;


namespace RdChannelMedianFilter {

  VModule::ResultFlag
  RdChannelMedianFilter::Init()
  {
    Branch topBranch = CentralConfig::GetInstance()->GetTopBranch("RdChannelMedianFilter");
    topBranch.GetChild("InfoLevel").GetData(fInfoLevel);
    topBranch.GetChild("FilterBandwidth").GetData(fFilterBandwidth);
    topBranch.GetChild("FilterSizeBins").GetData(fFilterSizeBins);
    topBranch.GetChild("MedianSortAlgorithmLimit").GetData(fMedianSortAlgorithmLimit);

    return eSuccess;
  }


  VModule::ResultFlag
  RdChannelMedianFilter::Run(evt::Event& event)
  {
    // Check if there are events at all
    if(!event.HasREvent()) {
      WARNING("No radio event found!");
      return eContinueLoop;
    }

    REvent& rEvent = event.GetREvent();
    ostringstream info;

    // loop through stations and for every station through every channel
    for (auto& station : rEvent.StationsRange()) {
      for (auto& channel : station.ChannelsRange()) {

        if (!channel.IsActive())
          continue;

        info.str("");
        info << "Supressing RFI of channel " << channel.GetId()
             << " of station " << channel.GetStationId() << ".";
        INFODebug(info);

        // Work on the frequency spectrum of this channel
        revt::ChannelFrequencySpectrum& spectrum = channel.GetChannelFrequencySpectrum();
        MedianFilter(spectrum);
      }
    }

    return eSuccess;
  }


  VModule::ResultFlag
  RdChannelMedianFilter::Finish()
  {
    return eSuccess;
  }


  void
  RdChannelMedianFilter::MedianFilter(revt::ChannelFrequencySpectrum& spectrum)
  {
    // get the amplitude
    Trace<double> amplitude;
    for (const auto& val : spectrum)
      amplitude.PushBack(abs(val));

    // check if filter bandwidth was provided instead of a filter size in bins,
    // and calculate size in bins from bandwidth
    // use provided size in bins as default, and modify if bandwidth > 0
    unsigned int filterSizeBins = fFilterSizeBins;
    if (fFilterBandwidth > 0) {
      filterSizeBins = trunc(fFilterBandwidth/spectrum.GetBinning()/2.) * 2 + 1;

      ostringstream info;
      info << "Use median filter size of " << filterSizeBins
           << " bins, corresponding to a bandwidth of "
           << filterSizeBins*spectrum.GetBinning()/megahertz << " MHz.";
      INFOIntermediate(info);
    }

    // get the filter size and check if it is small enough
    if (filterSizeBins > amplitude.GetSize()) {
      ostringstream info;
      info << "Size of median filter (" << fFilterSizeBins
           << " samples) should not be larger than the size of the spectrum ("
           << amplitude.GetSize() << " samples)";
      WARNING(info);
      filterSizeBins = amplitude.GetSize();
    }

    // index for looping through the spectrum
    ChannelFrequencySpectrum::SizeType i(0);

    // use constant median for the first part of the trace
    double median = TraceAlgorithm::Median(amplitude, 0, filterSizeBins-1, fMedianSortAlgorithmLimit);

    // adjust amplitude of the spectrum to the median amplitude
    // if the amplitude of a frequency bin is 0 it means that the
    // phase is abitrary and must be defined while setting the amplitude
    // to the median. At the moment phase = 0 is used,
    // but an interpolation might be better.
    // For real spectra there are hopefully no frequencies where the
    // amplitude is exactly 0, so the problem should not be severe.
    for (i = 0; i < ceil(filterSizeBins/2.0); ++i)
      if (abs(spectrum[i]) > 0)
        spectrum[i] *= median/abs(spectrum[i]);
      else spectrum[i] = median;

    // use gliding median for the middle part
    for(;i < (spectrum.GetSize()-ceil(filterSizeBins/2.0)); ++i)
      if (abs(spectrum[i]) > 0)
        spectrum[i] *=
          TraceAlgorithm::Median(amplitude, i-trunc(filterSizeBins/2.0),
                                 i+ceil(filterSizeBins/2.0)-1, fMedianSortAlgorithmLimit)/abs(spectrum[i]);
      else
        spectrum[i] =
          TraceAlgorithm::Median(amplitude, i-trunc(filterSizeBins/2.0),
                                 i+ceil(filterSizeBins/2.0)-1, fMedianSortAlgorithmLimit);

    // use constant median for the end of the spectrum
    median = TraceAlgorithm::Median(amplitude, spectrum.GetSize()-filterSizeBins,
                                    spectrum.GetSize()-1, fMedianSortAlgorithmLimit);
    for(;i < spectrum.GetSize(); ++i)
      if (abs(spectrum[i]) > 0)
        spectrum[i] *= median/abs(spectrum[i]);
      else spectrum[i] = median;
  }

}