source: Sophya/trunk/Poubelle/DPC:FitsIOServer/NTools/fftserver.cc@ 2306

#include "fftserver.h"
#include <iostream.h>

extern "C" {
  int rffti_(int *n, float *wsave);             //fftpack functions
  int rfftf_(int *n, float *r__, float *wsave);
  int rfftb_(int *n, float *r__, float *wsave);
  int dffti_(int *n, double *wsave);
  int dfftf_(int *n, double *r__, double *wsave);
  int dfftb_(int *n, double *r__, double *wsave);
  int cffti_(int *n, float *wsave);
  int cfftf_(int *n, float *c__, float *wsave);
  int cfftb_(int *n, float *c__, float *wsave);
  int cdffti_(int *n, double *wsave);
  int cdfftf_(int *n, double *c__, double *wsave);
  int cdfftb_(int *n, double *c__, double *wsave);
           }

FFTServer::FFTServer()
{
  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;
}

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

void FFTServer::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 FFTServer::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 FFTServer::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 FFTServer::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 FFTServer::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 FFTServer::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 FFTServer::fftf(int l, complex<float>* inout)
{
  checkint_cfft(l);
  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 FFTServer::fftf(int l, complex<double>* inout)
{
  checkint_cdfft(l);
  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 FFTServer::fftb(int l, float* inout)
{
  checkint_rfft(l);
  rfftf_(&l, inout, ws_rfft);
}

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

void FFTServer::fftb(int l, complex<float>* inout)
{
  checkint_cfft(l);
  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 FFTServer::fftb(int l, complex<double>* inout)
{
  checkint_cdfft(l);
  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;
}

void FFTServer::fftf(Vector& in, Vector& out)
{
  int l = in.NElts();
/* ----- Si c'etait un Vector<float> on aurait ecrit comme ca 
   Pour le moment Vector est double, on passe donc par un tableau 
   intermediare
// La transformee sur le tableau de flaot se fait en place,
// on utilise donc out comme in-out
  out = in;
  fftf_float(l, out.Data() );
  ------------------------------------------------------------- */
  float * inout = new float[l];
  int i;
  for(i=0; i<l; i++) inout[i] = in(i);
  fftf(l, inout);
  out.Realloc(l);
  for(i=0; i<l; i++) out(i) = inout[i];
}

void FFTServer::fftb(Vector& in, Vector& out)
{
  int l = in.NElts();
/* ----- Si c'etait un Vector<float> on aurait ecrit comme ca 
   Pour le moment Vector est double, on passe donc par un tableau 
   intermediare
// La transformee sur le tableau de flaot se fait en place,
// on utilise donc out comme in-out
  out = in;
  fftf_float(l, out.Data() );
  ------------------------------------------------------------- */
  float * inout = new float[l];
  int i;
  for(i=0; i<l; i++) inout[i] = in(i);
  fftb(l, inout);
  out.Realloc(l);
  for(i=0; i<l; i++) out(i) = inout[i];
}