FFTDataContainer.cc

Go to the documentation of this file.

#include <utl/FFTDataContainer.h>
#include <utl/ErrorLogger.h>

using namespace std;


// the specialisations of the template functions have to be defined in namespace utl

namespace utl {

  template<>
  void
  FFTDataContainer<Trace, double, complex<double>>::UpdateTimeSeries()
    const
  {
    const unsigned int np = fFrequencySpectrum.GetSize();  // number of points in frequency domain
    const unsigned int n = (np - 1) * 2;  // number of points in time domain

    if (!np) {  // empty frequency spectrum
      fUpToDateDomain = eBoth;  // now both domains are up to date
      WARNING("Warning: trying to do FFT on empty channel frequency spectrum!");
      return;
    }

    double* const f = FFTWdouble(n);  // allocate memory through fftw, making sure it is correctly aligned
    complex<double>* const g = FFTWComplex(np);

    fftwpp::crfft1d Backward(n, g, f);  // "plan" fft^-1 through fftw

    // fill fftw array, normalising on the way - maybe implement a more efficient way through copying ("reverse adopt") in Trace
    for (unsigned int i = 0; i < np; ++i)
      g[i] = fFrequencySpectrum[i] * fFrequencySpectrum.GetBinning();
    // should this explicitly "try" for exceptions?

    Backward.fft(g, f);

    fTimeSeries.Adopt(f, n);  // fill result into our data structure
    fTimeSeries.SetBinning(1. / (fFrequencySpectrum.GetBinning() * n));

    FFTWdelete(f/*, n*/);
    FFTWdelete(g/*, np*/);

    fUpToDateDomain = eBoth;  // now both domains are up to date
  }


  template<>
  void
  FFTDataContainer<Trace, double, complex<double>>::UpdateFrequencySpectrum()
    const
  {
    // check for an uneven number of samples and clip one sample if this is the case
    unsigned int nop = fTimeSeries.GetSize();  // number of points in time domain
    if ((nop / 2)*2 != nop) {
      fTimeSeries.PopBack();
      nop = fTimeSeries.GetSize();
    }

    if (!nop) {  // empty frequency spectrum
      fUpToDateDomain = eBoth;  // now both domains are up to date
      WARNING("Warning: trying to do FFT on empty channel time series!");
      return;
    }

    const unsigned int n = fTimeSeries.GetSize();  // number of points in time domain
    const unsigned int np = n/2 + 1;  // number of points in frequency domain

    double* const f = FFTWdouble(n);  // allocate memory through fftw, making sure it is correctly aligned
    complex<double>* const g = FFTWComplex(np);

    fftwpp::rcfft1d Forward(n, f, g);  // "plan" fft through fftw

    // fill fftw array, normalising on the way - maybe implement a more efficient way through copying ("reverse adopt") in Trace
    for (unsigned int i = 0; i < n; ++i)
      f[i] = fTimeSeries[i] * fTimeSeries.GetBinning();
    // should this explicitly "try" for exceptions?

    Forward.fft(f, g);

    fFrequencySpectrum.Adopt(g, np);        // fill result into our data structure
    fFrequencySpectrum.SetBinning(0.5 / (fTimeSeries.GetBinning() * (np - 1)));

    FFTWdelete(f/*, n*/);
    FFTWdelete(g/*, np*/);

    fUpToDateDomain = eBoth;  // now both domains are up to date
  }


  template<>
  void
  FFTDataContainer<Trace, Vector3D, Vector3C>::UpdateTimeSeries()
    const
  {
    const unsigned int np = fFrequencySpectrum.GetSize();  // number of points in frequency domain
    const unsigned int n = (np - 1) * 2; // number of points in time domain

    if (!np) {  // empty frequency spectrum
      fUpToDateDomain = eBoth;  // now both domains are up to date
      WARNING("Warning: trying to do FFT on empty station frequency spectrum!");
      return;
    }

    double* const f = FFTWdouble(n);  // allocate memory through fftw, making sure it is correctly aligned
    complex<double>* const g = FFTWComplex(np);
    Vector3D* const resultTimeSeries = new Vector3D[n];  // allocate temporary array of Vector3D to store the result time series

    fftwpp::crfft1d Backward(n, g, f);  // "plan" fft^-1 through fftw

    // now do the three sub-ffts
    for (int j = 0; j < 3; ++j) {
      // fill fftw array, normalising on the way - maybe implement a more efficient way through copying ("reverse adopt") in Trace
      for (unsigned int i = 0; i < np; ++i)
        g[i] = fFrequencySpectrum[i][j] * fFrequencySpectrum.GetBinning();
      // should this explicitly "try" for exceptions?
      Backward.fft(g, f);
      for (unsigned int i = 0; i < n; ++i)
        resultTimeSeries[i][j] = f[i];  // fill result into our data structure
    }
    fTimeSeries.Adopt(resultTimeSeries, n);
    fTimeSeries.SetBinning(1. / (fFrequencySpectrum.GetBinning() * n));

    FFTWdelete(f/*, n*/);
    FFTWdelete(g/*, np*/);
    delete[] resultTimeSeries;

    fUpToDateDomain = eBoth;  // now both domains are up to date
  }


  // ORIGINAL SERIAL FUNCTION
  template<>
  void
  FFTDataContainer<Trace, Vector3D, Vector3C>::UpdateFrequencySpectrum()
    const
  {
    // check for an uneven number of samples and clip one sample if this is the case
    unsigned int nop = fTimeSeries.GetSize();  // number of points in time domain
    if ((nop / 2)*2 != nop) {
      fTimeSeries.PopBack();
      nop = fTimeSeries.GetSize();
    }

    if (!nop) {  // empty frequency spectrum
      fUpToDateDomain = eBoth;  // now both domains are up to date
      WARNING("Warning: trying to do FFT on empty station time series!");
      return;
    }

    const unsigned int n = fTimeSeries.GetSize();  // number of points in time domain
    const unsigned int np = n/2 + 1;  // number of points in frequency domain

    double* const f = FFTWdouble(n);  // allocate memory through fftw, making sure it is correctly aligned
    complex<double>* const g = FFTWComplex(np);
    Vector3C* const resultSpectrum = new Vector3C[np];  // allocate temporary array of Vector3C to store the result spectrum

    fftwpp::rcfft1d Forward(n, f, g);  // "plan" fft through fftw

    // now do the three sub-ffts
    for (int j = 0; j < 3; ++j) {
      // fill fftw array, normalising on the way - maybe implement a more efficient way through copying ("reverse adopt") in Trace
      for (unsigned int i = 0; i < n; ++i)
        f[i] = fTimeSeries[i][j] * fTimeSeries.GetBinning();
      // should this explicitly "try" for exceptions?
      Forward.fft(f, g);
      for (unsigned int i = 0; i < np; ++i)
        resultSpectrum[i][j] = g[i];  // fill result into our data structure
    }
    fFrequencySpectrum.Adopt(resultSpectrum, np);
    fFrequencySpectrum.SetBinning(0.5 / (fTimeSeries.GetBinning() * (np - 1)));

    FFTWdelete(f/*, n*/);
    FFTWdelete(g/*, np*/);
    delete[] resultSpectrum;

    fUpToDateDomain = eBoth;  // now both domains are up to date
  }

}