RdChannelBandstopFilter.cc

Go to the documentation of this file.

#include "RdChannelBandstopFilter.h"

#include <fwk/CentralConfig.h>

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

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

#include <det/Detector.h>
#include <rdet/RDetector.h>
#include <rdet/Station.h>


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


namespace RdChannelBandstopFilter {

  VModule::ResultFlag
  RdChannelBandstopFilter::Init()
  {
    Branch topBranch = CentralConfig::GetInstance()->GetTopBranch("RdChannelBandstopFilter");
    topBranch.GetChild("InfoLevel").GetData(fInfoLevel);
    topBranch.GetChild("LowerBoundsNS").GetData(fLowerBoundsNS);
    topBranch.GetChild("UpperBoundsNS").GetData(fUpperBoundsNS);
    topBranch.GetChild("LowerBoundsEW").GetData(fLowerBoundsEW);
    topBranch.GetChild("UpperBoundsEW").GetData(fUpperBoundsEW);
    topBranch.GetChild("UseOnlineBandstopFilter").GetData(fUseOnlineBandstopFilter);
    topBranch.GetChild("OnlineBandstopFilterSpectrumBaseline").GetData(fOnlineBandstopFilterSpectrumBaseline);
    topBranch.GetChild("OnlineBandstopFilterDeviationFactor").GetData(fOnlineBandstopFilterDeviationFactor);

    // check for valid bounds in xml
    if (fLowerBoundsNS.size() != fUpperBoundsNS.size()) {
      ERROR("not same amount of upper/lower bounds in NS channel");
      return eFailure;
    }

    for(unsigned int i = 0; i < fLowerBoundsNS.size(); ++i) {
      if(fLowerBoundsNS[i] > fUpperBoundsNS[i]) {
        ERROR("NoiseStartFrequency greater than NoiseStopFrequency for NS channel");
        return eFailure;
      }
    }

    if (fLowerBoundsEW.size() != fUpperBoundsEW.size()) {
      ERROR("not same amount of upper/lower bounds in EW channel");
      return eFailure;
    }

    for(unsigned int i = 0; i < fLowerBoundsEW.size(); ++i) {
      if(fLowerBoundsEW[i] > fUpperBoundsEW[i]) {
        ERROR("NoiseStartFrequency greater than NoiseStopFrequency for EW channel 2");
        return eFailure;
      }
    }

    return eSuccess;
  }


  VModule::ResultFlag
  RdChannelBandstopFilter::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;

        const rdet::RDetector& rDetector = det::Detector::GetInstance().GetRDetector();
        const rdet::Channel& detChannel = rDetector.GetStation(station.GetId()).GetChannel(channel.GetId());

        // Work on the frequency spectrum of this channel
        ChannelFFTDataContainer& myFFTDataContainer = channel.GetFFTDataContainer();
        ChannelFrequencySpectrum& freqSpectrum = myFFTDataContainer.GetFrequencySpectrum();

        // online filtering for time-varying narrow-band RFI
        if (fUseOnlineBandstopFilter) {
          Trace<double> amplitude;
          for (const auto& freq : freqSpectrum)
            amplitude.PushBack(abs(freq));

          // get design frequencies
          const double lowerfreq = detChannel.GetDesignLowerFreq();
          const double upperfreq = detChannel.GetDesignUpperFreq();

          const unsigned int idxLowFreq = myFFTDataContainer.GetBinOfFrequency(lowerfreq);
          const unsigned int idxUpFreq = myFFTDataContainer.GetBinOfFrequency(upperfreq);

          double baseline = 0;  // fOnlineBandstopFilterSpectrumBaseline == "Zero"
          if (fOnlineBandstopFilterSpectrumBaseline == "Mean")
            baseline = TraceAlgorithm::Mean(amplitude, idxLowFreq, idxUpFreq);
          else if (fOnlineBandstopFilterSpectrumBaseline == "Median")
            baseline = TraceAlgorithm::Median(amplitude, idxLowFreq, idxUpFreq, 50);

          const double deviation = TraceAlgorithm::StandardDeviation(amplitude, idxLowFreq, idxUpFreq);

          info.str("");
          info << "Channel " << channel.GetId()
               << " in Station " << channel.GetStationId()
               << " has a baseline of: " << baseline
               << " (" << fOnlineBandstopFilterSpectrumBaseline << ") "
               << " and a std. deviation of: " << deviation
               << " with a deviation factor of: " << fOnlineBandstopFilterDeviationFactor << ".";
          INFODebug(info);

          for (auto& frequency : freqSpectrum) {
            if (abs(frequency) > baseline + fOnlineBandstopFilterDeviationFactor * deviation) {
              info.str("");
              info << "Suppress frequency " << frequency / MHz
                   << " of Channel " << channel.GetId()
                   << " in Station " << channel.GetStationId()
                   << " with a deviation factor of " << abs(abs(frequency) - baseline) / deviation
                   << ".";
              INFOIntermediate(info);
              frequency = 0;
            }
          }
        }

        //fixme TH: frequencies specified in xml file are always suppressed in addition to the dynamic filter!?
        CutNoiseFrequencies(channel, detChannel);
      }
    }

    return eSuccess;
  }


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


  void
  RdChannelBandstopFilter::CutNoiseFrequencies(revt::Channel& channel, const rdet::Channel& detChannel)
  {
    // get orientation of channel
    vector<double> lowerBounds, upperBounds;
    if (CloseTo(detChannel.GetOrientationAzimuth(), fNorthSouth, fTolerance)) {
      // North-South channel
      lowerBounds = fLowerBoundsNS;
      upperBounds = fUpperBoundsNS;
    } else if (CloseTo(detChannel.GetOrientationAzimuth(), fEastWest, fTolerance)) {
      // East-West channel
      lowerBounds = fLowerBoundsEW;
      upperBounds = fUpperBoundsEW;
    } else {
      INFOIntermediate("Channel not NS or EW oriented, not cutting out frequencies");
      return;
    }

    ostringstream info;

    ChannelFFTDataContainer& fftContainer = channel.GetFFTDataContainer();
    ChannelFrequencySpectrum& freqSpectrum = fftContainer.GetFrequencySpectrum();

    for(unsigned int i = 0; i < lowerBounds.size(); ++i) {
      const int idxNoiseStart = fftContainer.GetBinOfFrequency(lowerBounds[i]);
      const int idxNoiseStop = fftContainer.GetBinOfFrequency(upperBounds[i]);

      info << "Cut out noise frequency from " << lowerBounds[i] / MHz << " MHz to "
           << upperBounds[i] / MHz << " MHz (Bin " << idxNoiseStart<< " to " << idxNoiseStop
           << ") in Channel " << channel.GetId() << endl;
      INFOIntermediate(info);

      for(int j=idxNoiseStart; j<=idxNoiseStop; ++j) {
        freqSpectrum[j] = 0;
      }
    }
  }


  bool
  RdChannelBandstopFilter::CloseTo(const double a, const double b, const double tolerance)
  {
    return fabs(a - b) < tolerance;
  }

}