RdChannelSpectrumAverager.cc

Go to the documentation of this file.

#include "RdChannelSpectrumAverager.h"

#include <cmath>
#include <string>

#include <evt/Event.h>

#include <fwk/CentralConfig.h>

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

#include <revt/REvent.h>
#include <revt/Header.h>
#include <revt/Station.h>
#include <revt/Channel.h>

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

namespace RdChannelSpectrumAverager {

  RdChannelSpectrumAverager::RdChannelSpectrumAverager()
  {
  }

  RdChannelSpectrumAverager::~RdChannelSpectrumAverager()
  {
  }

  VModule::ResultFlag
  RdChannelSpectrumAverager::Init()
  {
    // Initialize your module here. This method
    // is called once at the beginning of the run.
    // The eSuccess flag indicates the method ended
    // successfully.  For other possible return types,
    // see the VModule documentation.

    INFO("RdChannelSpectrumAverager::Init()");

    // Read in the configurations of the xml file
    // Branch topBranch =
    //   CentralConfig::GetInstance()->GetTopBranch("RdChannelSpectrumAverager");
    // topBranch.GetChild("").GetData(f);

    return eSuccess;
  }

  VModule::ResultFlag
  RdChannelSpectrumAverager::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();

    // loop through each station and each channel of the event
    for (REvent::StationIterator sIt = rEvent.StationsBegin();
         sIt != rEvent.StationsEnd(); ++sIt) {
      Station& station = *sIt;

      for (Station::ChannelIterator cIt = station.ChannelsBegin();
           cIt != station.ChannelsEnd(); ++cIt) {
        Channel& channel = *cIt;
        stringstream fMessage("");
        fMessage << "Replacing amplitude spectrum by mean ampltitude spectrum for channel "
                 << channel.GetId()
                 << " of station "
                 << channel.GetStationId()
                 << ".";
        INFO(fMessage.str());

        // construct the key for accesing the maps
        pair<int,int> key(station.GetId(), channel.GetId());

        // get the spectrum and calculate the amplitude spectrum
        revt::ChannelFrequencySpectrum& spectrum =
          channel.GetChannelFrequencySpectrum();

        Trace<double> amplitude;
        for (ChannelFrequencySpectrum::SizeType i = 0; i < spectrum.GetSize(); ++i)
          amplitude.PushBack(abs(spectrum[i]));

        // check if it is first run for this station and channel
        // (then the ID does not exist in the counter map)
        if (fCounter.find(key) == fCounter.end()) {
          fCounter[key] = 1;  // set counter to 1 for this channel and station
          fMeanAmplitude[key] = amplitude;  // current amplitude spectrum is also the mean spectrum
        } else {
          // check if size of the spectrum fits and then calculate the mean amplitude
          if (fMeanAmplitude[key].GetSize() != amplitude.GetSize()) {
            ERROR("Size of current spectrum and mean spectrum do not match!");
            return eContinueLoop;
          }

          // increase counter and calculate new mean value
          ++fCounter[key];
          fMeanAmplitude[key] = (amplitude + double(fCounter[key]-1)*fMeanAmplitude[key]) / double(fCounter[key]);

          // change amplitude of spectrum to mean amplitude
          // spectral bins with an amplitude of 0 will be set to a phase 0, otherwise the phase is not touched.
          for (ChannelFrequencySpectrum::SizeType i = 0; i < amplitude.GetSize(); ++i)
            if (amplitude[i] > 0)
              spectrum[i] *= fMeanAmplitude[key][i]/amplitude[i];
//               spectrum[i] = fMeanAmplitude[key][i];
            else
              spectrum[i] = fMeanAmplitude[key][i];
        }
      }
    }

    return eSuccess;
  }

  VModule::ResultFlag
  RdChannelSpectrumAverager::Finish()
  {
    // Put any termination or cleanup code here.
    // This method is called once at the end of the run.

    INFO("RdChannelSpectrumAverager::Finish()");
    return eSuccess;
  }
}