FFTDataContainerAlgorithm.cc

Go to the documentation of this file.

#include <utl/FFTDataContainerAlgorithm.h>
#include <utl/Math.h>


using namespace std;

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

  template<>
  void
  FFTDataContainerAlgorithm::HilbertTransform(FFTDataContainer<Trace, double, complex<double> >& container)
  {
    // Get the time series
    Trace<double>& timeseries = container.GetTimeSeries();
    const unsigned int n = timeseries.GetSize();  // number of points in time and frequency domain

    if (!n) {                            // empty frequency spectrum
      WARNING("Warning: trying to do a Hilbert transform on an empty time series!");
      return;
    }

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

    fftwpp::fft1d Forward(n, -1, fftin, fftout);               // plan fft through fftw
    fftwpp::fft1d Backward(n, 1, fftout, fftin);               // plan fft^-1 through fftw

    for (unsigned int i = 0; i < n; ++i) {
      fftin[i].real(timeseries[i]);
      fftin[i].imag(0.);
    }
    Forward.fft(fftin, fftout);

    for (unsigned int i = 0; i < n; ++i) {
      double tmp;
      if (i < (n / 2)) {
        tmp = imag(fftout[i]);
        fftout[i].imag(-real(fftout[i]));
      } else {
        tmp = -imag(fftout[i]);
        fftout[i].imag(real(fftout[i]));
      }
      fftout[i].real(tmp);
    }

    Backward.fft(fftout, fftin);

    for (unsigned int i = 0; i < n; ++i)
      timeseries[i] = real(fftin[i])/n;

    FFTWdelete(fftin/*, n*/);
    FFTWdelete(fftout/*, n*/);
  }

  template<>
  void
  FFTDataContainerAlgorithm::HilbertTransform(FFTDataContainer<Trace, Vector3D, Vector3C>& container)
  {
    // Get the time series
    Trace<Vector3D>& timeseries = container.GetTimeSeries();
    const unsigned int n = timeseries.GetSize();  // number of points in time and frequency domain

    if (!n) {                            // empty frequency spectrum
      WARNING("Warning: trying to do a Hilbert envelope on an empty time series!");
      return;
    }

    std::complex<double>* const fftin = FFTWComplex(n);   // allocate memory through fftw, making sure it is correctly aligned
    std::complex<double>* const fftout = FFTWComplex(n);

    fftwpp::fft1d Forward(n, -1, fftin, fftout);               // fft through fftw
    fftwpp::fft1d Backward(n, 1, fftout, fftin);               // fft^-1 through fftw

    for (unsigned int j = 0; j < 3; ++j) {

      for (unsigned int i = 0; i < n; ++i) {
        fftin[i].real(timeseries[i][j]);
        fftin[i].imag(0.);
      }
      Forward.fft(fftin, fftout);

      for (unsigned int i = 0; i < n; ++i) {
        double tmp;
        if (i < (n / 2)) {
          tmp = imag(fftout[i]);
          fftout[i].imag(-real(fftout[i]));
        } else {
          tmp = -imag(fftout[i]);
          fftout[i].imag(real(fftout[i]));
        }
        fftout[i].real(tmp);
      }

      Backward.fft(fftout, fftin);

      for (unsigned int i = 0; i < n; ++i)
        timeseries[i][j] = real(fftin[i])/n;

    }

    FFTWdelete(fftin/*, n*/);
    FFTWdelete(fftout/*, n*/);
  }

  template<>
  void
  FFTDataContainerAlgorithm::HilbertEnvelope(FFTDataContainer<Trace, double, complex<double> >& container)
  {
    FFTDataContainer<Trace, double, complex<double> > transformedContainter = container;

    HilbertTransform(transformedContainter);
    Trace<double>& timeseries =container.GetTimeSeries();
    Trace<double>& transformedTimeseries =transformedContainter.GetTimeSeries();
    const unsigned int n = timeseries.GetSize();
    for (unsigned int i = 0; i < n; ++i)
      timeseries[i] = sqrt(Sqr(timeseries[i]) + Sqr(transformedTimeseries[i]));
  }

  template<>
  void
  FFTDataContainerAlgorithm::HilbertEnvelope(FFTDataContainer<Trace, Vector3D, Vector3C>& container)
  {

    FFTDataContainer<Trace, Vector3D, Vector3C> transformedContainter = container;

    HilbertTransform(transformedContainter);

    Trace<Vector3D>& timeseries = container.GetTimeSeries();
    const unsigned int n = timeseries.GetSize();
    Trace<Vector3D>& transformedTimeseries =transformedContainter.GetTimeSeries();
    // could also be done on GPU
    for (unsigned int j = 0; j < 3; ++j) {
      for (unsigned int i = 0; i < n; ++i) {
        timeseries[i][j] = sqrt(Sqr(timeseries[i][j]) + Sqr(transformedTimeseries[i][j]));
      }
    }
  }

}