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

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

#include <iostream>


/*! 
  \class SOPHYA::FFTPackServer
  \ingroup NTools
  An implementation of FFTServerInterface based on fftpack, for 
  one dimensional arrays.

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.

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.

  \code
  #include "fftpserver.h"
  // ...
  TVector<r_8> in(32);
  TVector< complex<r_8> > out;
  in = RandomSequence();
  FFTPackServer ffts;
  ffts.setNormalize(true);  // To have normalized transforms
  cout << " FFTServer info string= " << ffts.getInfo() << endl;
  cout << "in= " << in << endl;
  cout << " Calling ffts.FFTForward(in, out) : " << endl;
  ffts.FFTForward(in, out);
  cout << "out= " << out << endl;
  \endcode
*/


//! Constructor - If preserve_input==true, input arrays will not be overwritten.
FFTPackServer::FFTPackServer(bool preserve_input)
   : FFTServerInterface("FFTPackServer using extended FFTPack (C-version) package")
  , ckR4("FFTPackServer: ", true, true) , ckR8("FFTPackServer: ", true, true),
  _preserve_input(preserve_input)
{
//the working array and its size for the different
//possible numerical types
  sz_rfft = 0;     
  ws_rfft = NULL;  
  sz_dfft = 0;     
  ws_dfft = NULL;  
  sz_cfft = 0;
  ws_cfft = NULL;
  sz_cdfft = 0;
  ws_cdfft = NULL;
}

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

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


void FFTPackServer::FFTForward(TArray< complex<r_8> > & in, TArray< complex<r_8> > & out)
{
  ckR8.CheckResize(in, out);
  out = in;
  fftf(out.Size(), out.Data());
  if (getNormalize()) out *= (1./(r_8)(in.Size()));
}

void FFTPackServer::FFTBackward(TArray< complex<r_8> > & in, TArray< complex<r_8> > & out)
{
  ckR8.CheckResize(in, out);
  out = in;
  fftb(out.Size(), out.Data());
}


void FFTPackServer::FFTForward(TArray< complex<r_4> > & in, TArray< complex<r_4> > & out)
{
  ckR4.CheckResize(in, out);  
  out = in;
  fftf(out.Size(), out.Data());
  if (getNormalize()) out *= (1./(r_4)(in.Size()));
}

void FFTPackServer::FFTBackward(TArray< complex<r_4> > & in, TArray< complex<r_4> > & out)
{
  ckR4.CheckResize(in, out);  
  out = in;
  fftb(out.Size(), out.Data());
}

void FFTPackServer::FFTForward(TArray< r_4 > & in, TArray< complex<r_4> > & out)
{
  ckR4.CheckResize(in, out);
  bool share = (_preserve_input) ? false : true;
  TArray< r_4 > inout(in, share);
  fftf(inout.Size(), inout.Data());
  ReShapetoCompl(inout, out);
  if (getNormalize()) out *= complex<r_4>((1./(r_4)(in.Size())), 0.);
}

void FFTPackServer::FFTBackward(TArray< complex<r_4> > & in, TArray< r_4 > & out,
                                bool usoutsz)
{
  ckR4.CheckResize(in, out, usoutsz);    
  ReShapetoReal(in, out, usoutsz);
  fftb(out.Size(), out.Data());
}


void FFTPackServer::FFTForward(TArray< r_8 > & in, TArray< complex<r_8> > & out)
{
  ckR8.CheckResize(in, out);  
  bool share = (_preserve_input) ? false : true;
  TArray< r_8 > inout(in, share);
  fftf(inout.Size(), inout.Data());
  ReShapetoCompl(inout, out);
  if (getNormalize()) out *= complex<r_8>((1./(r_8)(in.Size())), 0.);
}

void FFTPackServer::FFTBackward(TArray< complex<r_8> > & in, TArray< r_8 > & out,
                                bool usoutsz)
{
  ckR8.CheckResize(in, out, usoutsz);  
  ReShapetoReal(in, out, usoutsz);
  fftb(out.Size(), out.Data());
}


template <class T>
void FFTPack_ReShapetoReal(TArray< complex<T> > const & ina, TArray< T >  & outa, bool usz)
{
  TVector< complex<T> > in(ina);
  TVector< T > out(outa);
  sa_size_t n = in.NElts();
  sa_size_t ncs;
  if (usz) {
    if ( (out.NElts() != 2*n-1) && (out.NElts() != 2*n-2) )
      throw SzMismatchError("FFTPack_ReShapetoReal(..., true) - Wrong output array size ");
    ncs = out.NElts();
  }
  else {
    T thr = FFTArrayChecker<T>::ZeroThreshold();
    ncs = ( (in(n-1).imag() < -thr) || (in(n-1).imag() > thr) ) ? 2*n-1 : 2*n-2;
    
    if (out.NElts() != ncs) {
      cerr << "DEBUG-FFTPack_ReShapetoReal() ncs = " << ncs 
           << " out.NElts()= " << out.NElts() << endl;
      throw SzMismatchError("FFTPack_ReShapetoReal() - Wrong output array size !");
    }
  }

  sa_size_t 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(); }

  return;
}

template <class T>
void FFTPack_ReShapetoCompl(TArray< T > const & ina, TArray< complex<T> > & outa)
{
  TVector< T > in(ina);
  TVector< complex<T> > out(outa);
  sa_size_t n = in.NElts();
  sa_size_t ncs = n/2+1;
  sa_size_t nc = (n%2 != 0) ? n/2+1 : n/2;
  if (out.NElts() != ncs) {
    cerr << "DBG-ReShapetoCompl() ncs=" << ncs 
         << " out.NElts()= " << out.NElts() << endl;
    throw SzMismatchError("FFTPack_ReShapetoCompl() - Wrong output array size !");
  }
  out(0) = complex<T> (in(0),0.);
  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<T>(in(n-1), 0.);

  return;
}

void FFTPackServer::ReShapetoReal(TArray< complex<r_8> > const & in, TArray< r_8 >  & out, 
                                  bool usz)
{
  FFTPack_ReShapetoReal<r_8>(in, out, usz);
}

void FFTPackServer::ReShapetoCompl(TArray< r_8 > const & in, TArray< complex<r_8> > & out)
{
  FFTPack_ReShapetoCompl<r_8>(in, out);
}

void FFTPackServer::ReShapetoReal(TArray< complex<r_4> > const & in, TArray< r_4 >  & out,
                                  bool usz)
{
  FFTPack_ReShapetoReal<r_4>(in, out, usz);
}

void FFTPackServer::ReShapetoCompl(TArray< r_4 > const & in, TArray< complex<r_4> > & out)
{
  FFTPack_ReShapetoCompl<r_4>(in, out);
}

void FFTPackServer::checkint_rfft(int_4 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_4 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_4 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_4 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_4 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_4 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_4 l, complex<r_4>* inout)
{
  checkint_cfft(l);
  cfftf_(&l, (r_4 *)(inout), ws_cfft);
}

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

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

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

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

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