hilbert_transform_inline.hpp

/* hilbert_transform_inline.hpp */
/* $Id: hilbert_transform_inline.hpp 21022 2011-07-11 11:45:03Z urij $ */
#ifndef __hilbert_transform_inline_hpp
#define __hilbert_transform_inline_hpp

template <typename FFT_T>
bool Hilbert::transform(FFT_T &fft, Trace &trace, int direction)
{
  // check direction
  if (direction != -1 && direction != 1)
    {
      printf ("S2Corr::hilbertTransform: warning: incorrect direction (%d), fixed to 1.\n",
          direction);
      direction = 1;
    }
  
  // FFT with rotated frequency scale, where f=0 at the center
  fft.setRotateFlag(true); 
  
#warning proper t0 calc should be fixed in CDFT2D_FFTW
  double t0 = trace.header().t0;
  
  // calc spectrum
  Trace trace_im(trace.header(), trace.size(), 0.);
  fft.setDirection(1);
  //fft.setNormType(CDFT2D::INORM);
  if (!fft.calc(trace, trace_im))
    {
      printf ("S2Corr::hilbertTransform: error: can't calc direct FFT.\n");
      return false;
    }
  
  // calc Hilbert transform in spectrum domain
  for (size_t i=0; i<trace.size(); i++)
    {
      double w = trace.header().t0 + i*trace.header().timestep;
      double sig_w = (w>0)?1:(-1);
      if (w == 0.)
    sig_w = 0.;
      
      double re = direction*sig_w*trace_im[i];
      trace_im[i] = -direction*sig_w*trace[i];
      trace[i] = re;
    }
  
  // inverse Fourier
  fft.setDirection(-1);
  if (!fft.calc(trace, trace_im))
    {
      printf ("S2Corr::hilbertTransform: error: can't calc inverse FFT.\n");
      return false;
    }
  trace.header().t0 = t0;

  return true;
}


#endif /* hilbert_transform_inline.hpp */