RdChannelBeaconSuppressor.cc

Go to the documentation of this file.

#include "RdChannelBeaconSuppressor.h"

#include <set>

#include <fwk/CentralConfig.h>

#include <utl/ErrorLogger.h>
#include <utl/Reader.h>
#include <utl/config.h>
#include <utl/TimeStamp.h>
#include <utl/UTCDateTime.h>

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

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



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

namespace RdChannelBeaconSuppressor {

  VModule::ResultFlag
  RdChannelBeaconSuppressor::Init()
  {
    INFO("RdChannelBeaconSuppressor::Init()");

    // Read in the configurations of the xml file
    const Branch topBranch =
      CentralConfig::GetInstance()->GetTopBranch("RdChannelBeaconSuppressor");
    topBranch.GetChild("InfoLevel").GetData(fInfoLevel);
    
    topBranch.GetChild("Frequencies").GetData(fXMLFrequencies);
    topBranch.GetChild("GetBeaconFrequenciesFromDatabase").GetData(fGetBeaconFrequenciesFromDatabase);
    topBranch.GetChild("SurpressionFactor").GetData(fSurpressionFactor);
    topBranch.GetChild("NeighbouringBinsToBeSurpressed").GetData(fNeighbouringBinsToBeSurpressed);
    topBranch.GetChild("NeighbouringBandwidthToBeSurpressed").GetData(fNeighbouringBandwidthToBeSurpressed);
    topBranch.GetChild("CompensateTotalSpectrumPower").GetData(fCompensateTotalSpectrumPower);
    return eSuccess;
  }


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

    REvent& rEvent = event.GetREvent();

    // copy frequencies provided in the XML file to a local vector
    std::vector<double> frequencies = fXMLFrequencies;
    
    // add beacon frequencies of database to those specified in the xml file
    if (fGetBeaconFrequenciesFromDatabase) {
      det::Detector& detector = det::Detector::GetInstance();
      const rdet::RDetector& RDetector = detector.GetRDetector();
      vector<double> BeaconFrequencies = RDetector.GetBeaconFrequencies();
      for (auto it = BeaconFrequencies.begin(); it != BeaconFrequencies.end(); ++it)
        frequencies.push_back(*it);
    }  

    // loop through stations and for every station through every channel
    for (auto sIt = rEvent.StationsBegin(); sIt != rEvent.StationsEnd(); ++sIt) {
      Station& station = *sIt;

      std::set<int> frequencyBins;
    
      for (auto cIt = station.ChannelsBegin(); cIt != station.ChannelsEnd(); ++cIt) {
        Channel& channel = *cIt;
        if (!channel.IsActive())
          continue;
        stringstream message;
        message << "Suppressing sine waves of channel " << channel.GetId()
                << " of station " << channel.GetStationId() << ".";
        INFOIntermediate(message.str());

        auto& spectrum = channel.GetChannelFrequencySpectrum();
            
        //fixme TH: the following should be refactored to use FFTDataContainer::GetBinOfFrequency() method
        
        // loop through frequencies and find corresponding frequency bin in the spectrum 
        for (unsigned int freqN=0; freqN < frequencies.size(); ++freqN) {
          // check if frequency is inside the range (take care that the order of the spectrum is unknown)
          if(  (frequencies[freqN] <
                min(channel.GetFrequencyOfBin(0),channel.GetFrequencyOfBin(spectrum.GetSize()-1)))
            || (frequencies[freqN] >
                max(channel.GetFrequencyOfBin(0),channel.GetFrequencyOfBin(spectrum.GetSize()-1))) ) {
            ERROR("Frequency must be inside of the spectrum!");
            return eFailure;
          }
          
          // find nearest frequency bin:
          double bin = (frequencies[freqN]-channel.GetFrequencyOfBin(0)) * (spectrum.GetSize()-1)
                    / (channel.GetFrequencyOfBin(spectrum.GetSize()-1) - channel.GetFrequencyOfBin(0));
                    
          unsigned int intBin = round(bin);                    
                    
          // consistency check, whether bin is approximately integer
          if (abs( double(intBin)/bin -1.) > 1e-3) {
            double nearestIntFreq = channel.GetFrequencyOfBin(intBin);
            stringstream message;
            message << "Frequency does not correspond to approxemately integer bin of spectrum! The nearest integer frequency is "
                    << nearestIntFreq / MHz  << " MHz (bin #" << intBin << ")" ;
            INFODebug(message.str());
          }
          
          // store frequencies of the closest spectrum bin
          frequencyBins.insert(intBin);
          
          // also surpress surounding bins if specified:
          // calculate number of bins from bandwidth if specified
          if (fNeighbouringBandwidthToBeSurpressed > 0)         
            fNeighbouringBinsToBeSurpressed = round(fNeighbouringBandwidthToBeSurpressed/spectrum.GetBinning());

          // add neighbouring bins, but do not allow for bins outside of the spectrum
          for (unsigned int i = 1; i <= fNeighbouringBinsToBeSurpressed; ++i) {
            if ( intBin >= i)
              frequencyBins.insert(intBin-i);
            if ( (intBin+i) < spectrum.GetSize())
              frequencyBins.insert(intBin+i);
          } 
        }
        
        // count how many bins are in the design bandwidth (for later compensation of surpression)
        unsigned int numberOfbinsInDesignBandwidth = 0;
        double lowerDesignFreq = 0;
        double upperDesignFreq = 0;
        // for performance: only fill values, when compensation will be done later
        if (fCompensateTotalSpectrumPower) {
          det::Detector& detector = det::Detector::GetInstance();
          const rdet::Channel& detChannel 
            = detector.GetRDetector().GetStation(station.GetId()).GetChannel(channel.GetId());
          lowerDesignFreq = detChannel.GetDesignLowerFreq();
          upperDesignFreq = detChannel.GetDesignUpperFreq();          }
        
        for (auto it = frequencyBins.begin(); it != frequencyBins.end(); ++it) {
          message.str("");
          message << "Suppressing spectral bin " << *it << " corresponding to frequency "
                    << channel.GetFrequencyOfBin(*it)/megahertz << " MHz.";
          INFODebug(message.str());

          spectrum[*it] *= fSurpressionFactor; // surpress spectral bin
          if (fCompensateTotalSpectrumPower)
            if( (channel.GetFrequencyOfBin(*it) >= lowerDesignFreq) && (channel.GetFrequencyOfBin(*it) <= upperDesignFreq) )
              ++numberOfbinsInDesignBandwidth;
        }
        
        // calculate compensation: the power of a possible signal in the design bandwidth should be preserved
        // For this, a flat frequency spectrum of the signal is assumend, which corresponds to a first order correction
        // the correction compensates for the fraction of surpressed bins within the design bandwidth
        if (fCompensateTotalSpectrumPower) {
           // get maximum frequency of spectrum (depending on Nyquist zone could be the first or last bin of the spectrum)
          double maxFrequency = 0;
          if (maxFrequency<channel.GetFrequencyOfBin(0))
            maxFrequency = channel.GetFrequencyOfBin(0);
          if (maxFrequency<channel.GetFrequencyOfBin(spectrum.GetSize()-1))
            maxFrequency = channel.GetFrequencyOfBin(spectrum.GetSize()-1);
          
          const double relativeBandWidth = (upperDesignFreq-lowerDesignFreq)/maxFrequency;         
          // total power should be constant, but compensation is applied on amplitude spectrum
          double compensationFactor = sqrt(1.+double(numberOfbinsInDesignBandwidth)/double(spectrum.GetSize())/relativeBandWidth);
            
          message.str("");
          message << "Multiplying each bin in the spectrum by " << compensationFactor
                  << " such that the total signal power within the design bandwidth is preserved.";
          INFOIntermediate(message.str());
          
          // apply compensation factor
          for (ChannelFrequencySpectrum::SizeType i = 0; i < spectrum.GetSize(); ++i)
            spectrum[i]*=compensationFactor;         
        }
      }
    }

    return eSuccess;
  }


  VModule::ResultFlag
  RdChannelBeaconSuppressor::Finish()
  {
    INFOFinal("RdChannelBeaconSuppressor::Finish()");

    return eSuccess;
  }

}