source: Sophya/trunk/SophyaLib/NTools/fftpserver.cc@ 710

#include "fftpserver.h"
#include "fftpackc.h"

#include <iostream.h>


/*! \class FFTPackServer
 \brief An FFT server based on fftpack

A class that calculates Fourier transforms forwards and backwards.

The class calls the c library ``fftpack'', which is accessible and documented
at http://www.netlib.org/fftpack/.  However, the class functions do not
necessarily correspond with the equivalent fftpack function.  For example,
fftpack "forward" transformations are in fact inverse fourier transformations.
Otherwise, the output is in the fftpack format.

Complex input must be entered using the <complex> template.  Otherwise,
all input is assumed to be real.

Due to the way that fftpack manages
its work arrays, an object can run faster if the length of the input arrays
does not change.  For example, if you need to do a series of FFT's
of differing length, it may be more efficient to create an fftserver object
for each length.
*/

/*! \fn virtual void FFTServer::fftf(int l, float* inout)
  \param l length of array
  \param inout input array /output forward FFT (original array destroyed)
*/
/*! \fn virtual void FFTServer::fftb(int l, float* inout)
  \param l length of array
  \param inout input array /output backward FFT (original array destroyed)
*/
/*! \fn virtual void FFTServer::fftf(int l, double* inout)
  \param l length of array
  \param inout input array /output forward FFT (original array destroyed)
  \param inout input/output array (original array destroyed)
*/
/*! \fn virtual void FFTServer::fftb(int l, double* inout)
  \param l length of array
  \param inout input array /output backward FFT(original array destroyed)
*/
/*!\fn  virtual void FFTServer::fftf(int l, complex<float>* inout)
  \param l length of array
  \param inout input array /output forward FFT (original array destroyed)
*/
/*! \fn virtual void FFTServer::fftb(int l, complex<float>* inout)
  \param l length of array
  \param inout input array /output backward FFT (original array destroyed)
*/
/*! \fn virtual void FFTServer::fftf(int l, complex<double>* inout)
  \param l length of array
  \param inout input array /output forward FFT (original array destroyed)
*/
/*! \fn virtual void FFTServer::fftb(int l, complex<double>* inout)
  \param l length of array
  \param inout input array /output backward FFT(original array destroyed)
*/
/*!\fn  virtual void FFTServer::fftf(Vector& in, Vector& out)
  \param in input array
  \param out forward FFT
*/
/*! \fn virtual void FFTServer::fftb(Vector& in, Vector& out)
  \param in input array
  \param out backward FFT
*/

FFTPackServer::FFTPackServer()
   : FFTServerInterface("FFTServer using FFTPack (C-version) package")

{
  sz_rfft = 0;     //the working array and its size for the different
  ws_rfft = NULL;  //possible numerical types
  sz_cfft = 0;
  ws_cfft = NULL;
  sz_cdfft = 0;
  ws_cdfft = NULL;
}

FFTPackServer::~FFTPackServer()
{
if (ws_rfft) delete[] ws_rfft;
if (ws_cfft) delete[] ws_cfft;
if (ws_cdfft) delete[] ws_cdfft;
}

FFTServerInterface * FFTPackServer::Clone()
{
  return (new FFTPackServer);
}

void FFTPackServer::FFTForward(TVector< complex<r_8> > const & in, TVector< complex<r_8> > & out)
{
  out = in;
  fftf(out.NElts(), out.Data());
}

void FFTPackServer::FFTBackward(TVector< complex<r_8> > const & in, TVector< complex<r_8> > & out)
{
  out = in;
  fftb(out.NElts(), out.Data());
}

void FFTPackServer::FFTForward(TVector< complex<r_4> > const & in, TVector< complex<r_4> > & out)
{
  out = in;
  fftf(out.NElts(), out.Data());
}

void FFTPackServer::FFTBackward(TVector< complex<r_4> > const & in, TVector< complex<r_4> > & out)
{
  out = in;
  fftb(out.NElts(), out.Data());
}

void FFTPackServer::FFTForward(TVector< r_4 > const & in, TVector< complex<r_4> > & out)
{
  TVector< r_4 > inout(in);
  fftf(inout.NElts(), inout.Data());
  ReShapetoCompl(inout, out);
}

void FFTPackServer::FFTBackward(TVector< complex<r_4> > const & in, TVector< r_4 > & out)
{
  ReShapetoReal(in, out);
  fftb(out.NElts(), out.Data());
}

void FFTPackServer::FFTForward(TVector< r_8 > const & in, TVector< complex<r_8> > & out)
{
  TVector< r_8 > inout(in);
  fftf(inout.NElts(), inout.Data());
  ReShapetoCompl(inout, out);
}

void FFTPackServer::FFTBackward(TVector< complex<r_8> > const & in, TVector< r_8 > & out)
{
  ReShapetoReal(in, out);
  fftb(out.NElts(), out.Data());
}


void FFTPackServer::checkint_rfft(int l)
{
  if (sz_rfft == l) return;       //checkint functions check and reallocate
                                  //memory for the work arrays when performing
  if (ws_rfft) delete[] ws_rfft;  //a transform
  sz_rfft = l;
  ws_rfft = new float[2*l+15];
  rffti_(&l, ws_rfft);
}

void FFTPackServer::checkint_cfft(int l)
{
  if (sz_cfft == l) return;

  if (ws_cfft) delete[] ws_cfft;
  sz_cfft = l;
  ws_cfft = new float[4*l+15];
  cffti_(&l, ws_cfft);
}

void FFTPackServer::checkint_dfft(int l)
{
  if (sz_dfft == l) return;

  if (ws_dfft) delete[] ws_dfft;
  sz_dfft = l;
  ws_dfft = new double[2*l+15];
  dffti_(&l, ws_dfft);
}

void FFTPackServer::checkint_cdfft(int l)
{
  if (sz_cdfft == l) return;

  if (ws_cdfft) delete[] ws_cdfft;
  sz_cdfft = l;
  ws_cdfft = new double[4*l+15];
  cdffti_(&l, ws_cdfft);
}

/* In general forward transformations are resorted since fftpack functions
   return inverse transformations */

void FFTPackServer::fftf(int l, float* inout)
{
  checkint_rfft(l);
  rfftf_(&l, inout, ws_rfft);
  for (int k= 2;k<=(l+1)/2;k++) inout[2*k-2]=-inout[2*k-2];
}

void FFTPackServer::fftf(int l, double* inout)
{
  checkint_dfft(l);
  dfftf_(&l, inout, ws_dfft);
  for (int k= 2;k<=(l+1)/2;k++) inout[2*k-2]=-inout[2*k-2];
}

void FFTPackServer::fftf(int l, complex<float>* inout)
{
  checkint_cfft(l);
  cfftf_(&l, (float *)(inout), ws_cfft);

  /*
  float* foo = new float[2*l];
  int i;
  for (i=0;i<l;i++){
    foo[2*i]=inout[i].real();
    foo[2*i+1]=inout[i].imag();
  }
  cfftf_(&l, foo, ws_cfft);
  inout[0]=complex<float> (foo[0],foo[1]);
  for (i=1;i<l;i++) inout[l-i]= complex<float> (foo[2*i], foo[2*i+1]);
  delete[] foo;
  */
}

void FFTPackServer::fftf(int l, complex<double>* inout)
{
  checkint_cdfft(l);
  cdfftf_(&l, (double*)(inout), ws_cdfft);
  /*
  double* foo=new double[2*l];
  int i;
  for (i=0;i<l;i++){
    foo[2*i]=inout[i].real();
    foo[2*i+1]=inout[i].imag();
  }
  cdfftf_(&l, foo, ws_cdfft);
  inout[0]=complex<double> (foo[0],foo[1]);
  for (i=1;i<l;i++) {
    inout[l-i]= complex<double> (foo[2*i],foo[2*i+1]);
  }
  delete[] foo;
  */
}

void FFTPackServer::fftb(int l, float* inout)
{
  checkint_rfft(l);
  rfftf_(&l, inout, ws_rfft);
}

void FFTPackServer::fftb(int l, double* inout)
{
  checkint_dfft(l);
  dfftf_(&l, inout, ws_dfft);
}

void FFTPackServer::fftb(int l, complex<float>* inout)
{
  checkint_cfft(l);
  cfftb_(&l, (float *)(inout), ws_cfft);
  /*
  float* foo = new float[2*l];
  int i;
  for (i=0;i<l;i++){
    foo[2*i]=inout[i].real();
    foo[2*i+1]=inout[i].imag();
  }
  cfftf_(&l, foo, ws_cfft);
  for (i=0;i<l;i++) inout[i]=complex<float> (foo[2*i],foo[2*i+1]);
  delete[] foo;
  */
}

void FFTPackServer::fftb(int l, complex<double>* inout)
{
  checkint_cdfft(l);
  cfftb_(&l, (float *)(inout), ws_cfft);
  /*
  double* foo = new double[2*l];
  int i;
  for (i=0;i<l;i++){
    foo[2*i]=inout[i].real();
    foo[2*i+1]=inout[i].imag();
  }
  cdfftf_(&l, foo, ws_cdfft);
  for (i=0;i<l;i++) inout[i]=complex<double> (foo[2*i],foo[2*i+1]);
  delete[] foo;
  */
}