source: Sophya/trunk/SophyaLib/NTools/fftmserver.cc@ 3195

#include "sopnamsp.h"
#include "fftmserver.h"
#include <iostream>
#include "fftmayer.h"


/* --Methode-- */
FFTMayerServer::FFTMayerServer()
   : FFTServerInterface("FFTMayerServer using extended FFTMayer package")
  , ckR4("FFTMayerServer: ", true, true) , ckR8("FFTMayerServer: ", true, true)

{
}

/* --Methode-- */
FFTMayerServer::~FFTMayerServer()
{
}

/* --Methode-- */
FFTServerInterface * FFTMayerServer::Clone()
{
  return (new FFTMayerServer);
}


/* --Methode-- */
void FFTMayerServer::FFTForward(TArray< complex<r_8> > const & ina, TArray< complex<r_8> > & outa)
{
 ckR8.CheckResize(ina, outa);
 TVector< complex<r_8> > in(ina);
 TVector< complex<r_8> > out(outa);

 r_8 a,b,c,d;
 r_8 q,r,s,t;
 int i,j,k;
 int n = in.NElts();
 if (checkLength(n))  {
   char buff[256];
   sprintf(buff, "FFTMayerServer::FFTForward/Error Size(%d) != 2^n ! (File=%s Line=%d)",
           n, __FILE__,__LINE__);
   throw ParmError(string(buff));
 }
 TVector< r_8 > inoutre(n);
 TVector< r_8 > inoutim(n);
 inoutre(0) = in(0).real();
 inoutim(0) = in(0).imag();
 inoutre(n/2) = in(n/2).real();
 inoutim(n/2) = in(n/2).imag();
 if (n > 2) {
   inoutre(n/2+1) = in(n/2).real();
   inoutim(n/2+1) = in(n/2).imag();
 }
 for (i=1,j=n-1,k=n/2;i<k;i++,j--) {
  a = in(i).real(); b = in(j).real();  q=a+b; r=a-b;
  c = in(i).imag(); d = in(j).imag();  s=c+d; t=c-d;
  inoutre(i) = (q+t)*.5; inoutre(j) = (q-t)*.5;
  inoutim(i) = (s-r)*.5; inoutim(j) = (s+r)*.5;
 }
 fht_r8(inoutre.Data(),inoutre.NElts());
 fht_r8(inoutim.Data(),inoutim.NElts());

 // out.ReSize(in.NElts());
 r_8 fn = 1./n;
 if (getNormalize()) 
   for(k=0; k<out.NElts(); k++)
     out(k) = complex<r_8>(inoutre(k)*fn, inoutim(k)*fn);
 else 
   for(k=0; k<out.NElts(); k++)
     out(k) = complex<r_8>(inoutre(k), inoutim(k));
}

/* --Methode-- */
void FFTMayerServer::FFTBackward(TArray< complex<r_8> > const & ina, TArray< complex<r_8> > & outa)
{
 ckR8.CheckResize(ina, outa);
 TVector< complex<r_8> > in(ina);
 TVector< complex<r_8> > out(outa);

 r_8 a,b,c,d;
 r_8 q,r,s,t;
 int i,j,k;
 int n = in.NElts();
 if (checkLength(n))  {
   char buff[256];
   sprintf(buff, "FFTMayerServer::FFTBackward/Error Size(%d) != 2^n ! (File=%s Line=%d)",
           n, __FILE__,__LINE__);
   throw ParmError(string(buff));
 }
 TVector< r_8 > inoutre(n);
 TVector< r_8 > inoutim(n);
 for(k=0; k<out.NElts(); k++) {
   inoutre(k) = in(k).real();
   inoutim(k) = in(k).imag();
 }
 fht_r8(inoutre.Data(),inoutre.NElts());
 fht_r8(inoutim.Data(),inoutim.NElts());

 // out.ReSize(in.NElts());
 out(0) = complex<r_8>(inoutre(0), inoutim(0));
 out(n/2) = complex<r_8>(inoutre(n/2), inoutim(n/2));
 if (n > 2) out(n/2+1) = complex<r_8>(inoutre(n/2+1), inoutim(n/2+1));

 for (i=1,j=n-1,k=n/2;i<k;i++,j--) {
  a = inoutre(i); b = inoutre(j);  q=a+b; r=a-b;
  c = inoutim(i); d = inoutim(j);  s=c+d; t=c-d;
  out(i) = complex<r_8>( (q-t)*0.5 , (s+r)*0.5 );
  out(j) = complex<r_8>( (q+t)*0.5 , (s-r)*0.5 );
 }

}

/* --Methode-- */
void FFTMayerServer::FFTForward(TArray< complex<r_4> > const & ina, TArray< complex<r_4> > & outa)
{
 ckR4.CheckResize(ina, outa);
 TVector< complex<r_4> > in(ina);
 TVector< complex<r_4> > out(outa);

 r_4 a,b,c,d;
 r_4 q,r,s,t;
 int i,j,k;
 int n = in.NElts();
 if (checkLength(n))  {
   char buff[256];
   sprintf(buff, "FFTMayerServer::FFTForward/Error Size(%d) != 2^n ! (File=%s Line=%d)",
           n, __FILE__,__LINE__);
   throw ParmError(string(buff));
 }
 TVector< r_4 > inoutre(n);
 TVector< r_4 > inoutim(n);
 inoutre(0) = in(0).real();
 inoutim(0) = in(0).imag();
 inoutre(n/2) = in(n/2).real();
 inoutim(n/2) = in(n/2).imag();
 if (n > 2) {
   inoutre(n/2+1) = in(n/2).real();
   inoutim(n/2+1) = in(n/2).imag();
 }
 for (i=1,j=n-1,k=n/2;i<k;i++,j--) {
  a = in(i).real(); b = in(j).real();  q=a+b; r=a-b;
  c = in(i).imag(); d = in(j).imag();  s=c+d; t=c-d;
  inoutre(i) = (q+t)*.5; inoutre(j) = (q-t)*.5;
  inoutim(i) = (s-r)*.5; inoutim(j) = (s+r)*.5;
 }
 fht_r4(inoutre.Data(),inoutre.NElts());
 fht_r4(inoutim.Data(),inoutim.NElts());

 // out.ReSize(in.NElts());
 r_4 fn = 1./n;
 if (getNormalize()) 
   for(k=0; k<out.NElts(); k++)
     out(k) = complex<r_4>(inoutre(k)*fn, inoutim(k)*fn);
 else 
   for(k=0; k<out.NElts(); k++)
     out(k) = complex<r_4>(inoutre(k), inoutim(k));

}

/* --Methode-- */
void FFTMayerServer::FFTBackward(TArray< complex<r_4> > const & ina, TArray< complex<r_4> > & outa)
{
 ckR4.CheckResize(ina, outa);
 TVector< complex<r_4> > in(ina);
 TVector< complex<r_4> > out(outa);

 r_4 a,b,c,d;
 r_4 q,r,s,t;
 int i,j,k;
 int n = in.NElts();
 if (checkLength(n))  {
   char buff[256];
   sprintf(buff, "FFTMayerServer::FFTBackward/Error Size(%d) != 2^n ! (File=%s Line=%d)",
           n, __FILE__,__LINE__);
   throw ParmError(string(buff));
 }

 TVector< r_4 > inoutre(n);
 TVector< r_4 > inoutim(n);
 for(k=0; k<out.NElts(); k++) {
   inoutre(k) = in(k).real();
   inoutim(k) = in(k).imag();
 }
 fht_r4(inoutre.Data(),inoutre.NElts());
 fht_r4(inoutim.Data(),inoutim.NElts());

 // out.ReSize(in.NElts());
 out(0) = complex<r_4>(inoutre(0), inoutim(0));
 out(n/2) = complex<r_4>(inoutre(n/2), inoutim(n/2));
 if (n > 2) out(n/2+1) = complex<r_4>(inoutre(n/2+1), inoutim(n/2+1));

 for (i=1,j=n-1,k=n/2;i<k;i++,j--) {
  a = inoutre(i); b = inoutre(j);  q=a+b; r=a-b;
  c = inoutim(i); d = inoutim(j);  s=c+d; t=c-d;
  out(i) = complex<r_4>( (q-t)*0.5 , (s+r)*0.5 );
  out(j) = complex<r_4>( (q+t)*0.5 , (s-r)*0.5 );
 }

}

/* --Methode-- */
void FFTMayerServer::FFTForward(TArray< r_4 > const & ina, TArray< complex<r_4> > & outa)
{
 ckR4.CheckResize(ina, outa);
 TVector< r_4 > in(ina);
 TVector< complex<r_4> > out(outa);

 r_4 a,b;
 int i,j,k;
 int n = in.NElts();
 if (checkLength(n))  {
   char buff[256];
   sprintf(buff, "FFTMayerServer::FFTForward/Error Size(%d) != 2^n ! (File=%s Line=%d)",
           n, __FILE__,__LINE__);
   throw ParmError(string(buff));
 }

 TVector< r_4 > inout(in);
 fht_r4(inout.Data(),inout.NElts());
 for (i=1,j=n-1,k=n/2;i<k;i++,j--) {
  a = inout(i);  
  b = inout(j);
  inout(j) = (a-b)*0.5;
  inout(i) = (a+b)*0.5;
  out(i) = complex<r_4>( (a+b)*0.5 , (a-b)*0.5 );
 }

  
 // out.ReSize(n/2+1);
 out(0) = complex<r_4>(inout(0), 0.);

 if (n%2 == 0) { // n pair
   for (i=1,j=n/2+1; i<n/2; i++,j++) 
     out(i) = complex<r_4>(inout(i), inout(j));

   out(n/2) = complex<r_4>(inout(n/2), 0.);
 }
 else { // n impair
   for (i=1,j=n/2+2; i<n/2; i++,j++) 
     out(i) = complex<r_4>(inout(i), inout(j));

   out(n/2) = complex<r_4>(inout(n/2), inout(n/2+1));
 }

 if (getNormalize()) out *= (1./(r_4)(in.NElts()));
}

/* --Methode-- */
void FFTMayerServer::FFTBackward(TArray< complex<r_4> > const & ina, TArray< r_4 > & outa,
                                 bool usoutsz)
{
 ckR4.CheckResize(ina, outa);
 TVector< complex<r_4> > in(ina);
 TVector< r_4 > out(outa);

 r_4 a,b;
 int i,j,k;
 int n = in.NElts();
 int nc = (fabs(in(n-1).imag()) < 1.e-9) ? nc = n*2-2 : nc = n*2-1;
 if (checkLength(nc))  {
   char buff[256];
   sprintf(buff, "FFTMayerServer::FFTBackward/Error Size(%d) != 2^n ! (File=%s Line=%d)",
           nc, __FILE__,__LINE__);
   throw ParmError(string(buff));
 }

 // out.ReSize(nc);
 out(0) = in(0).real();
 if (nc%2 == 0) { // nc pair
   for (i=1,j=nc/2+1; i<nc/2; i++,j++) {
     out(i) = in(i).real();
     out(j) = in(i).imag();
   }
   out(nc/2) = in(nc/2).real();
 }
 else { // nc impair
   for (i=1,j=nc/2+2; i<nc/2; i++,j++) {
     out(i) = in(i).real();
     out(j) = in(i).imag();
   }
   out(nc/2) = in(nc/2).real();
   out(nc/2+1) = in(nc/2).imag();
 }

 for (i=1,j=nc-1,k=nc/2;i<k;i++,j--) {
  a = out(i);
  b = out(j);
  out(j) = (a-b);
  out(i) = (a+b);
 }
 fht_r4(out.Data(),out.NElts());
}


/* --Methode-- */
void FFTMayerServer::FFTForward(TArray< r_8 > const & ina, TArray< complex<r_8> > & outa)
{
 ckR8.CheckResize(ina, outa);
 TVector< r_8 > in(ina);
 TVector< complex<r_8> > out(outa);

 r_8 a,b;
 int i,j,k;
 int n = in.NElts();
 if (checkLength(n))  {
   char buff[256];
   sprintf(buff, "FFTMayerServer::FFTForward/Error Size(%d) != 2^n ! (File=%s Line=%d)",
           n, __FILE__,__LINE__);
   throw ParmError(string(buff));
 }

 TVector< r_8 > inout(in);
 fht_r8(inout.Data(),inout.NElts());
 for (i=1,j=n-1,k=n/2;i<k;i++,j--) {
  a = inout(i);  
  b = inout(j);
  inout(j) = (a-b)*0.5;
  inout(i) = (a+b)*0.5;
  out(i) = complex<r_8>( (a+b)*0.5 , (a-b)*0.5 );
 }

  
 // out.ReSize(n/2+1);
 out(0) = complex<r_8>(inout(0), 0.);

 if (n%2 == 0) { // n pair
   for (i=1,j=n/2+1; i<n/2; i++,j++) 
     out(i) = complex<r_8>(inout(i), inout(j));

   out(n/2) = complex<r_8>(inout(n/2), 0.);
 }
 else { // n impair
   for (i=1,j=n/2+2; i<n/2; i++,j++) 
     out(i) = complex<r_8>(inout(i), inout(j));

   out(n/2) = complex<r_8>(inout(n/2), inout(n/2+1));
 }
 if (getNormalize()) out *= (1./(r_4)(in.NElts()));
}

/* --Methode-- */
void FFTMayerServer::FFTBackward(TArray< complex<r_8> > const & ina, TArray< r_8 > & outa,
                                 bool usoutsz)
{
 ckR8.CheckResize(ina, outa);
 TVector< r_8 > out(outa);
 TVector< complex<r_8> > in(ina);

 r_8 a,b;
 int i,j,k;
 int n = in.NElts();
 int nc = (fabs(in(n-1).imag()) < 1.e-18) ? nc = n*2-2 : nc = n*2-1;
 if (checkLength(nc))  {
   char buff[256];
   sprintf(buff, "FFTMayerServer::FFTBackward/Error Size(%d) != 2^n ! (File=%s Line=%d)",
           nc, __FILE__,__LINE__);
   throw ParmError(string(buff));
 }

 // out.ReSize(nc);
 out(0) = in(0).real();
 if (nc%2 == 0) { // nc pair
   for (i=1,j=nc/2+1; i<nc/2; i++,j++) {
     out(i) = in(i).real();
     out(j) = in(i).imag();
   }
   out(nc/2) = in(nc/2).real();
 }
 else { // nc impair
   for (i=1,j=nc/2+2; i<nc/2; i++,j++) {
     out(i) = in(i).real();
     out(j) = in(i).imag();
   }
   out(nc/2) = in(nc/2).real();
   out(nc/2+1) = in(nc/2).imag();
 }

 for (i=1,j=nc-1,k=nc/2;i<k;i++,j--) {
  a = out(i);
  b = out(j);
  out(j) = (a-b);
  out(i) = (a+b);
 }
 fht_r8(out.Data(),out.NElts());
}



/* --Methode-- */
bool FFTMayerServer::checkLength(int n)
{
  if (n < 2) return(true);
  int nc = n;
  while (nc > 1) {
    if (nc%2 != 0)  return(true);
    nc /= 2;
  }
  return(false);
}