FFTDataContainer.h

Go to the documentation of this file.

#ifndef _utl_FFTDataContainer_h_
#define _utl_FFTDataContainer_h_

#include <complex>
#include <utl/Trace.h>
#include <utl/ShadowPtr.h>
#include <utl/AugerException.h>

#include <utl/fft++.h>


namespace utl {

  template<template<typename X> class C, typename T, typename F>
  class FFTDataContainer {

  public:
    enum Domain { eTime, eFrequency, eBoth };

    FFTDataContainer() = default;

    unsigned int GetNyquistZone() const { return fNyquistZone; }

    void
    SetNyquistZone(const unsigned int zone)
    {
      if (zone < 1)
        throw OutOfBoundException("Requested illegal Nyquist zone of zero.");
      fNyquistZone = zone;
    }

    double
    GetFrequencyOfBin(const typename C<F>::SizeType bin)
      const
    {
      if (fUpToDateDomain == eTime)
        UpdateFrequencySpectrum();  // do the fft
      const double binning = fFrequencySpectrum.GetBinning();
      const typename C<F>::SizeType highestbin = fFrequencySpectrum.GetSize() - 1;
      if (bin > highestbin)
        throw OutOfBoundException("Requested frequency for a bin that is not in the FrequencySpectrum.");
      const double bandwidth = double(highestbin) * binning;
      const unsigned int evenedzone = (fNyquistZone / 2) * 2;
      const double startingfrequency = evenedzone * bandwidth;
      // we are in an even Nyquist zone, i.e., the frequencies count downwards
      const int sign = (fNyquistZone == evenedzone) ? -1 : 1;
      return startingfrequency + sign * double(bin) * binning;
    }

    int
    GetBinOfFrequency(const double frequency)
      const
    {
      if (fUpToDateDomain == eTime)
        UpdateFrequencySpectrum(); // do the fft
      const double binning = fFrequencySpectrum.GetBinning();
      const typename C<F>::SizeType highestbin = fFrequencySpectrum.GetSize() - 1;
      const double bandwidth = double(highestbin) * binning;
      if ((fNyquistZone - 1)*bandwidth > frequency || fNyquistZone*bandwidth < frequency)
        throw OutOfBoundException("Requested bin for a frequency that is not in the FrequencySpectrum.");
      const unsigned int evenedzone = (fNyquistZone / 2) * 2;
      const double startingfrequency = evenedzone * bandwidth;
      // we are in an even Nyquist zone, i.e., the frequencies count downwards
      const int sign = (fNyquistZone == evenedzone) ? -1 : 1;
      const double timeBinning = fTimeSeries.GetBinning();
      const int timeTraceLength = fTimeSeries.GetSize();
      const double bin = (frequency - startingfrequency) * sign *  timeBinning * timeTraceLength;
      return std::round(bin);
    }

    const C<T>&
    GetConstTimeSeries()
      const
    {
      if (fUpToDateDomain == eFrequency)
        UpdateTimeSeries();  // do the fft^-1
      return fTimeSeries;
    }

    C<T>&
    GetTimeSeries()
    {
      if (fUpToDateDomain == eFrequency)
        UpdateTimeSeries(); // do the fft^-1
      fUpToDateDomain = eTime;  // write access is possible, so have to assume data will change
      return fTimeSeries;
    }

    const C<F>&
    GetConstFrequencySpectrum()
      const
    {
      if (fUpToDateDomain == eTime)
        UpdateFrequencySpectrum();  // do the fft
      return fFrequencySpectrum;
    }

    C<F>&
    GetFrequencySpectrum()
    {
      if (fUpToDateDomain == eTime)
        UpdateFrequencySpectrum();  // do the fft
      fUpToDateDomain = eFrequency;  // write access is possible, so have to assume data will change
      return fFrequencySpectrum;
    }

    void
    Clear()
    {
      fTimeSeries.Clear();
      fFrequencySpectrum.Clear();
      fUpToDateDomain = eBoth;  // no data there
    }

  private:
    mutable Domain fUpToDateDomain = eBoth;

    unsigned int fNyquistZone = 1;  // data corresponds to which Nyquist zone

    mutable C<T> fTimeSeries;
    mutable C<F> fFrequencySpectrum;

    void UpdateTimeSeries() const;
    void UpdateFrequencySpectrum() const;

  };

}


#endif