testFFTDataContainer.cc

Go to the documentation of this file.

#include <iostream>

#include <tst/Verify.h>

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

#include <utl/FFTDataContainer.h>
#include <utl/AugerException.h>
#include <utl/AugerUnits.h>

#include <cppunit/extensions/HelperMacros.h>

using namespace std;
using namespace utl;
using namespace tst;
using namespace revt;

class testFFTDataContainer : public CppUnit::TestFixture {

  static const long numtestsamples = 1999; // must be an uneven number!

  CPPUNIT_TEST_SUITE(testFFTDataContainer);
//CPPUNIT_TEST(testEqualAndCopyCtor);
  CPPUNIT_TEST(testChannelFFTs);
  CPPUNIT_TEST(testStationFFTs);
//CPPUNIT_TEST_EXCEPTION(testBinException, InvalidTraceBoundException);
  CPPUNIT_TEST_SUITE_END();

private:
  ChannelFFTDataContainer ourChannelContainer;
  StationFFTDataContainer ourStationContainer;

public:
  void
  setUp()
  {
  }

  void
  tearDown()
  {
  }

/*
  void
  testEqualAndCopyCtor()
  {
    ChannelFFTDataContainer copyFFTDataContainer(ourChannelContainer);
    CPPUNIT_ASSERT(copyFFTDataContainer == ourChannelContainer);
    CPPUNIT_ASSERT(!(copyFFTDataContainer != ourChannelContainer));
  }
*/
  void
  testChannelFFTs()
  {
   
    // put some test data into the time series
    ChannelTimeSeries& timeSeries = ourChannelContainer.GetTimeSeries();
    CPPUNIT_ASSERT(Verify<Equal>(timeSeries.GetSize(), ChannelTimeSeries::SizeType(0)));                // trace of fresh container should be empty
    for(unsigned int i=0; i<numtestsamples; i++)
    {
      timeSeries.PushBack(numtestsamples-i);
    }
    timeSeries.SetBinning(12.5*ns); // fictitious 12.5 ns spacing (80 MSPS)
    ourChannelContainer.SetNyquistZone(2); // 40 to 80 MHz

    // test reading the time trace
    const ChannelTimeSeries ourTimeSeries1 = ourChannelContainer.GetConstTimeSeries();
    CPPUNIT_ASSERT(Verify<Equal>(ourTimeSeries1.GetSize(), ChannelTimeSeries::SizeType(numtestsamples)));
    for (unsigned int i=0; i<numtestsamples; i++)
    {
      CPPUNIT_ASSERT(Verify<CloseTo>(ourTimeSeries1[i], static_cast<double>(numtestsamples-i)));
    }

    // test automatic forward transform
    const ChannelFrequencySpectrum ourFrequencySpectrum1 = ourChannelContainer.GetConstFrequencySpectrum();
    CPPUNIT_ASSERT(Verify<CloseTo>(ourChannelContainer.GetFrequencyOfBin(0),80.0*megahertz));
    CPPUNIT_ASSERT(Verify<CloseTo>(ourChannelContainer.GetFrequencyOfBin(ourFrequencySpectrum1.GetSize()-1),40.0*megahertz));
    
    // sum up energy in frequency spectrum
    double energyInFreq = 0.0;
    for (unsigned int i=0; i<ourFrequencySpectrum1.GetSize(); i++)
    {
      energyInFreq += pow(abs(ourFrequencySpectrum1[i]),2)*ourFrequencySpectrum1.GetBinning();
      if (i!=0)
        energyInFreq += pow(abs(ourFrequencySpectrum1[i]),2)*ourFrequencySpectrum1.GetBinning();  // negative frequencies
    }
    
    // access time series data without transform
    const ChannelTimeSeries ourTimeSeries2 = ourChannelContainer.GetConstTimeSeries();
    CPPUNIT_ASSERT(Verify<Equal>(ourTimeSeries2.GetSize(), ChannelTimeSeries::SizeType(numtestsamples-1)));  // automatic clipping of the trace done?
    for (unsigned int i=0; i<numtestsamples-1; i++)
    {
//    cout << ourTimeSeries2[i] << "\n";
      CPPUNIT_ASSERT(Verify<CloseTo>(ourTimeSeries2[i], static_cast<double>(numtestsamples-i)));
    }
    
    // test automatic backward transform
    ChannelFrequencySpectrum frequencySpectrum = ourChannelContainer.GetFrequencySpectrum(); // write access, invalidates time domain data
    // CAUTION: the user could now access the time series with the old reference ourTimeSeries2;
    // then the time domain data would not be updated although the frequency domain data was changed
    // in the last step; if the user asks for a new TimeSeries as done here, everything is fine
    const ChannelTimeSeries ourTimeSeries3 = ourChannelContainer.GetConstTimeSeries();
    CPPUNIT_ASSERT(Verify<CloseTo>(ourTimeSeries3.GetBinning(),12.5*ns));
    double energyInTime = 0.0;
    for (unsigned int i=0; i<numtestsamples-1; i++)
    {
//    cout << ourTimeSeries2[i] << "\n";
      CPPUNIT_ASSERT(Verify<CloseTo>(ourTimeSeries3[i], static_cast<double>(numtestsamples-i)));
          energyInTime += pow(ourTimeSeries3[i],2)*ourTimeSeries3.GetBinning(); // sum up energy in time domain
    }
    // check if equal energy in both domains
    CPPUNIT_ASSERT(Verify<CloseTo>(energyInTime,energyInFreq));
  }

  void
  testStationFFTs()
  {
    // put some test data into the time series
    StationTimeSeries& timeSeries = ourStationContainer.GetTimeSeries();
    CPPUNIT_ASSERT(Verify<Equal>(timeSeries.GetSize(), StationTimeSeries::SizeType(0)));                // trace of fresh container should be empty
    double tmpvector[3]={0.0,0.0,0.0};
    for(unsigned int i=0; i<numtestsamples; i++)
    {
      tmpvector[0]=(numtestsamples-i);
      tmpvector[1]=(numtestsamples-i)*2.0;
      tmpvector[2]=(numtestsamples-i)*3.0;
      timeSeries.PushBack(Vector3D(tmpvector));
    }
    timeSeries.SetBinning(12.5);

    // test reading the time trace
    const StationTimeSeries ourTimeSeries1 = ourStationContainer.GetConstTimeSeries();
    CPPUNIT_ASSERT(Verify<Equal>(ourTimeSeries1.GetSize(), StationTimeSeries::SizeType(numtestsamples)));
    for (unsigned int i=0; i<numtestsamples; i++)
    {
      CPPUNIT_ASSERT(Verify<CloseTo>(ourTimeSeries1[i][0], static_cast<double>(numtestsamples-i)));
      CPPUNIT_ASSERT(Verify<CloseTo>(ourTimeSeries1[i][1], static_cast<double>(numtestsamples-i)*2.0));
      CPPUNIT_ASSERT(Verify<CloseTo>(ourTimeSeries1[i][2], static_cast<double>(numtestsamples-i)*3.0));
    }

    // test automatic forward transform
    const StationFrequencySpectrum ourFrequencySpectrum1 = ourStationContainer.GetConstFrequencySpectrum();
    CPPUNIT_ASSERT(Verify<CloseTo>(ourStationContainer.GetFrequencyOfBin(0),0.0*megahertz));
    CPPUNIT_ASSERT(Verify<CloseTo>(ourStationContainer.GetFrequencyOfBin(ourFrequencySpectrum1.GetSize()-1),40.0*megahertz));
    
    // sum up energy in frequency spectrum
    double energyInFreq = 0.0;
    for (unsigned int i=0; i<ourFrequencySpectrum1.GetSize(); i++)
    {
      energyInFreq += ourFrequencySpectrum1[i].GetMag2()*ourFrequencySpectrum1.GetBinning();
      if (i!=0)
        energyInFreq += ourFrequencySpectrum1[i].GetMag2()*ourFrequencySpectrum1.GetBinning();  // negative frequencies
    }
    
    // access time series data without transform 
    const StationTimeSeries ourTimeSeries2 = ourStationContainer.GetConstTimeSeries();
    CPPUNIT_ASSERT(Verify<Equal>(ourTimeSeries2.GetSize(), StationTimeSeries::SizeType(numtestsamples-1)));  // automatic clipping of the trace done?
    for (unsigned int i=0; i<numtestsamples-1; i++)
    {
//    cout << ourTimeSeries2[i][0] << std::endl;
      CPPUNIT_ASSERT(Verify<CloseTo>(ourTimeSeries2[i][0], static_cast<double>(numtestsamples-i)));
      CPPUNIT_ASSERT(Verify<CloseTo>(ourTimeSeries2[i][1], static_cast<double>(numtestsamples-i)*2.0));
      CPPUNIT_ASSERT(Verify<CloseTo>(ourTimeSeries2[i][2], static_cast<double>(numtestsamples-i)*3.0));
    }
    
    // test automatic backward transform
    StationFrequencySpectrum frequencySpectrum = ourStationContainer.GetFrequencySpectrum();             // write access, invalidates time domain data
    // CAUTION: the user could now access the time series with the old reference ourTimeSeries2;
    // then the time domain data would not be updated although the frequency domain data was changed
    // in the last step; if the user asks for a new TimeSeries as done here, everything is fine
    const StationTimeSeries ourTimeSeries3 = ourStationContainer.GetConstTimeSeries();
    CPPUNIT_ASSERT(Verify<CloseTo>(ourTimeSeries3.GetBinning(),12.5*ns));
    double energyInTime = 0.0;
    for (unsigned int i=0; i<numtestsamples-1; i++)
    {
//    cout << ourTimeSeries2[i][0] << "\n";
      CPPUNIT_ASSERT(Verify<CloseTo>(ourTimeSeries3[i][0], static_cast<double>(numtestsamples-i)));
      CPPUNIT_ASSERT(Verify<CloseTo>(ourTimeSeries3[i][1], static_cast<double>(numtestsamples-i)*2.0));
      CPPUNIT_ASSERT(Verify<CloseTo>(ourTimeSeries3[i][2], static_cast<double>(numtestsamples-i)*3.0));
          energyInTime += ourTimeSeries3[i].GetMag2()*ourTimeSeries3.GetBinning(); // sum up energy in time domain
    }

    // check if equal energy in both domains
    CPPUNIT_ASSERT(Verify<CloseTo>(energyInTime,energyInFreq));
  }

/*
  void
  testBinException()
  {
  }
*/

};

CPPUNIT_TEST_SUITE_REGISTRATION(testFFTDataContainer);

// Configure (x)emacs for this file ...
// Local Variables:
// mode: c++
// End: