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

#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, r_4* inout)
  \param l length of array
  \param inout input array /output forward FFT (original array destroyed)
*/
/*! \fn virtual void FFTServer::fftb(int l, r_4* inout)
  \param l length of array
  \param inout input array /output backward FFT (original array destroyed)
*/
/*! \fn virtual void FFTServer::fftf(int l, r_8* 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, r_8* inout)
  \param l length of array
  \param inout input array /output backward FFT(original array destroyed)
*/
/*!\fn  virtual void FFTServer::fftf(int l, complex<r_4>* inout)
  \param l length of array
  \param inout input array /output forward FFT (original array destroyed)
*/
/*! \fn virtual void FFTServer::fftb(int l, complex<r_4>* inout)
  \param l length of array
  \param inout input array /output backward FFT (original array destroyed)
*/
/*! \fn virtual void FFTServer::fftf(int l, complex<r_8>* inout)
  \param l length of array
  \param inout input array /output forward FFT (original array destroyed)
*/
/*! \fn virtual void FFTServer::fftb(int l, complex<r_8>* 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("FFTPackServer using extended 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());
  if (getNormalize()) out *= (1./(r_8)(in.NElts()));
}

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());
  if (getNormalize()) out *= (1./(r_4)(in.NElts()));
}

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);
  if (getNormalize()) out *= (1./(r_4)(in.NElts()));
}

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);
  if (getNormalize()) out *= (1./(r_8)(in.NElts()));
}

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 r_4[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 r_4[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 r_8[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 r_8[4*l+15];
  cdffti_(&l, ws_cdfft);
}

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

void FFTPackServer::fftf(int l, r_4* 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, r_8* 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<r_4>* inout)
{
  checkint_cfft(l);
  cfftf_(&l, (r_4 *)(inout), ws_cfft);
}

void FFTPackServer::fftf(int l, complex<r_8>* inout)
{
  checkint_cdfft(l);
  cdfftf_(&l, (r_8*)(inout), ws_cdfft);
}

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

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

void FFTPackServer::fftb(int l, complex<r_4>* inout)
{
  checkint_cfft(l);
  cfftb_(&l, (r_4 *)(inout), ws_cfft);
}

void FFTPackServer::fftb(int l, complex<r_8>* inout)
{
  checkint_cdfft(l);
  cdfftb_(&l, (r_8 *)(inout), ws_cdfft);
}

// Methodes pour reordonner les donnees 

/* --Methode-- */
void FFTPackServer::ReShapetoReal( TVector< complex<r_8> > const & in, TVector< r_8 >  & out)
{
  int n = in.NElts();
  int ncs = (fabs(in(n-1).imag()) > 1.e-12) ? ncs = 2*n-1 : ncs = n*2-2; 
  out.ReSize(ncs);
  int k;
  out(0) = in(0).real();
  for(k=1;k<n-1;k++) {
    out(2*k-1) = in(k).real();
    out(2*k) = in(k).imag();
  }
  if (ncs == n*2-2)  out(ncs-1) = in(n-1).real();
  else { out(ncs-2) = in(n-1).real();  out(ncs-1) = in(n-1).imag(); }
}

/* --Methode-- */
void FFTPackServer::ReShapetoReal( TVector< complex<r_4> > const & in, TVector< r_4 >  & out)
{
  int n = in.NElts();
  int ncs = (fabs(in(n-1).imag()) > 1.e-12) ? ncs = 2*n-1 : ncs = n*2-2; 
  out.ReSize(ncs);
  int k;
  out(0) = in(0).real();
  for(k=1;k<n-1;k++) {
    out(2*k-1) = in(k).real();
    out(2*k) = in(k).imag();
  }
  if (ncs == n*2-2)  out(ncs-1) = in(n-1).real();
  else { out(ncs-2) = in(n-1).real();  out(ncs-1) = in(n-1).imag(); }
}


/* --Methode-- */
void FFTPackServer::ReShapetoCompl(TVector< r_8 > const & in, TVector< complex<r_8> > & out)
{
  uint_4 n =  in.NElts();
  uint_4 ncs = n/2+1;
  uint_4 nc = (n%2 != 0) ? n/2+1 : n/2;
  out.ReSize(ncs);
  out(0) = complex<r_8> (in(0),0.);
  int k;
  for(int k=1;k<nc;k++) 
    out(k) =  complex<r_4> (in(2*k-1), in(2*k));
  if (n%2 == 0) out(ncs-1) = complex<r_8>(in(n-1), 0.);

}

/* --Methode-- */
void FFTPackServer::ReShapetoCompl(TVector< r_4 > const & in, TVector< complex<r_4> > & out)
{
  uint_4 n =  in.NElts();
  uint_4 ncs = n/2+1;
  uint_4 nc = (n%2 != 0) ? n/2+1 : n/2;
  out.ReSize(ncs);
  out(0) = complex<r_4> (in(0),0.);
  int k;
  for(int k=1;k<nc;k++) 
    out(k) =  complex<r_4> (in(2*k-1), in(2*k));
  if (n%2 == 0) out(ncs-1) = complex<r_4>(in(n-1), 0.);
}