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

#include <iostream.h>


/*! 
  \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.
Otherwise, the output is in the fftpack format.


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.
*/


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

{
//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);
}


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, false);
  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, false);
  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_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);
}

// 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.);
  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.);
  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.);
}