FFTDataContainerAlgorithm.h

Go to the documentation of this file.

#ifndef _utl_FFTDataContainerAlgorithm_h_
#define _utl_FFTDataContainerAlgorithm_h_

#include <utl/fft++.h>
#include <cmath>
#include <complex>
#include <boost/rational.hpp>
#include <utl/AugerException.h>
#include <utl/FFTDataContainer.h>
#include <utl/MathConstants.h>
#include <utl/Test.h>
#include <utl/TimeInterval.h>


class testFFTDataContainerAlgorithm;  // forward declaration needed for declaration of friend class


namespace utl {

  class FFTDataContainerAlgorithm {
  public:

    template<template<typename X> class C, typename T, typename F>
    static
    void
    ShiftTimeSeries(FFTDataContainer<C, T, F>& container, const TimeInterval& shift)
    {
      // Get the frequency spectrum
      C<F>& spectrum = container.GetFrequencySpectrum();

      // multiply phase gradient to each value of the spectrum
      for (unsigned int i = 0; i < spectrum.GetSize(); ++i)
        spectrum[i] *=
          std::exp(-std::complex<double>(0,1) * kTwoPi * container.GetFrequencyOfBin(i) * shift.GetInterval());
    }

    template<template<typename X> class C, typename T, typename F>
    static
    void
    UpsampleTimeSeries(FFTDataContainer<C, T, F>& container, const unsigned int upsamplingFactor, const bool removeOffset)
    {
      // Get the frequency spectrum
      C<F>& spectrum = container.GetFrequencySpectrum();

      // remove constant offset from trace (bin 1)
      // store offset from trace (first bin) to later readd it (if wanted)
      const F offset = spectrum[0];
      spectrum[0] *= 0;

      // construct upsampled trace:
      // start with a new spectrum of zeros
      // then fill in the values of the old spectrum at the correct position
      // the bin size is reduced by the upsamplingFactor so the new
      // spectrum has a size of (originalSize-1)*upsamplingFactor+1
      // as there is only one first and one last value
      C<F> upsampledSpectrum((spectrum.GetSize()-1)*upsamplingFactor+1, spectrum.GetBinning(), 0);

      // as the values must be filled in in the proper order there are two cases
      // as the order of the frequencies depends on the Nyquist zone
      if (container.GetNyquistZone() == 1) {
        for (unsigned int i = 0; i < spectrum.GetSize(); ++i)
          upsampledSpectrum[i] = spectrum[i];
      } else if (container.GetNyquistZone() == 2) {           // WARNING: Check missing if data are sorted as assumed
        const int lastentry = 2*(spectrum.GetSize()-1);
        for (int i = spectrum.GetSize()-1; i >= 0; --i)
          upsampledSpectrum[lastentry-i] = std::conj(spectrum[i]);  // complex conjugate avoids time inversion
        upsampledSpectrum[lastentry] *= 0.5;
        upsampledSpectrum[spectrum.GetSize()-1] *= 0.5;
      } else {
        throw utl::DoesNotComputeException("Upsampling for Nyquist zones greater than 2 not implemented!");
      }

      // replace the original spectrum by the new one (after zero padding)
      spectrum.Swap(upsampledSpectrum);

      // Add old mean to retain offset, if substraction is not requested
      if (!removeOffset)
        spectrum[0] = offset;

      // the new spectrum is allways in the 1st Nyquist zone
      container.SetNyquistZone(1);
    }

    template<template<typename X> class C, typename T, typename F>
    static
    void
    zeroPad(FFTDataContainer<C, T, F>& container, const unsigned int newSize, const bool removeOffset)
    {
      // Get the frequency spectrum
      C<F>& spectrum = container.GetFrequencySpectrum();

      // remove constant offset from trace (bin 1)
      // store offset from trace (first bin) to later readd it (if wanted)
      const F offset = spectrum[0];
      spectrum[0] *= 0;

      // construct upsampled trace:
      // start with a new spectrum of zeros
      // then fill in the values of the old spectrum at the correct position
      // the bin size is reduced by the upsamplingFactor so the new
      // spectrum has a size of (originalSize-1)*upsamplingFactor+1
      // as there is only one first and one last value
      C<F> upsampledSpectrum(newSize, spectrum.GetBinning(), 0);

      if (container.GetNyquistZone() == 1) {
        for (unsigned int i = 0; i < spectrum.GetSize(); ++i)
          upsampledSpectrum[i] = spectrum[i];
      } else {
        throw utl::DoesNotComputeException("Zero Padding for Nyquist zones greater than 1 not implemented!");
      }

      // replace the original spectrum by the new one (after zero padding)
      spectrum.Swap(upsampledSpectrum);

      // Add old mean to retain offset, if substraction is not requested
      if (!removeOffset)
        spectrum[0] = offset;

      // the new spectrum is allways in the 1st Nyquist zone
      container.SetNyquistZone(1);
    }

    template<template<typename X> class C, typename T, typename F>
    static
    void
    ResampleTimeSeries(FFTDataContainer<C, T, F>& container, const double resampleBinning)
    {
      if (resampleBinning <= 0)
        throw utl::DoesNotComputeException("Requested illegal resampling factor, aborting ...");
      const long maxUpsamplingFactor = 500;
      // first we have to convert the time bases to rationals
      double origBinning = container.GetConstTimeSeries().GetBinning();
      double targetBinning = resampleBinning;
      const boost::rational<long> origBinningRational = ConvertDecimalDoubleToRational(origBinning, false);
      const boost::rational<long> targetBinningRational = ConvertDecimalDoubleToRational(targetBinning, false);
      const boost::rational<long> resamplingFactor = origBinningRational / targetBinningRational;
      if (resamplingFactor.numerator() > maxUpsamplingFactor)
        throw utl::DoesNotComputeException("Requested upsampling factor is too large, aborting ...");
      if (resamplingFactor.numerator() != 1)
        UpsampleTimeSeries(container, resamplingFactor.numerator(), false);
      if (resamplingFactor.denominator() != 1)
        ThinOutTimeSeries(container, resamplingFactor.denominator());
    }

    template<template<typename X> class C, typename T, typename F>
    static
    void
    ThinOutTimeSeries(FFTDataContainer<C, T, F>& container, const long thinOutFactor)
    {
      // Get the time series
      C<T>& timeseries = container.GetTimeSeries();
      C<T> thinnedtimeseries(static_cast<unsigned int>(timeseries.GetSize()/thinOutFactor), timeseries.GetBinning()*thinOutFactor);

      // fill the thinned out time series
      for (unsigned int i = 0; i < thinnedtimeseries.GetSize(); ++i)
        thinnedtimeseries[i] = timeseries[i * thinOutFactor];

      // overwrite original time series
      timeseries.Swap(thinnedtimeseries);
    }

    template<template<typename X> class C, typename T, typename F>
    static void HilbertEnvelope(FFTDataContainer<C, T, F>& container);

    template<template<typename X> class C, typename T, typename F>
    static
    void
    HilbertTransform(FFTDataContainer<C, T, F>& container);


  private:
    static
    boost::rational<long>
    ConvertDecimalDoubleToRational(double fpn, const bool secondcall)
    {
      const double abortprecision = 1e-5;  // limit number of significant digits
      const double orgfpn = fpn;
      long denominator = 1;                       // starting value for denominator is 1
      double fpn_int;                             // double to store the integer part of fpn
      while (std::modf(fpn, &fpn_int)/fpn > abortprecision) {
        if (Test<CloseTo>(fpn - fpn_int, 1.0))  // floating point rounding error, we are finished
          break;
        fpn *= 10.0;
        denominator *= 10;
      }

      if (secondcall)
        return boost::rational<long>(denominator, long(fpn+0.5));

      // if this is the first call, try the conversion again with the inverse
      const boost::rational<long> firstresult =
        boost::rational<long>(long(fpn+0.5), denominator);
      const boost::rational<long> secondresult =
        ConvertDecimalDoubleToRational(1.0/orgfpn, true);

      // return the result with the smaller denominator
      return (firstresult.denominator() < secondresult.denominator()) ? firstresult : secondresult;
    }

    friend class ::testFFTDataContainerAlgorithm;  // so that we can test private functions

  };

} // utl

#endif