Changeset 3359 in Sophya for trunk


Ignore:
Timestamp:
Oct 23, 2007, 12:20:08 PM (18 years ago)
Author:
ansari
Message:

Ajout methode version r_4 (float) a FFTWServer (sous condition de flag) - Reza 23/10/2007

Location:
trunk/SophyaExt/IFFTW
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • trunk/SophyaExt/IFFTW/fftw3server.cc

    r3000 r3359  
    1414//! Constructor - If preserve_input==true, input arrays will not be overwritten.
    1515FFTWServer::FFTWServer(bool preserve_input)
    16   : FFTServerInterface("FFTServer using FFTW3 package")
    17   , ckR8("FFTWServer - ", true, false),
    18   _preserve_input(preserve_input) 
     16#ifdef ALSO_FFTW_FLOAT_EXTSOP
     17  : FFTServerInterface("FFTServer using FFTW3 package (r_4 and r_8)") ,
     18#else
     19  : FFTServerInterface("FFTServer using FFTW3 package (r_8 only)") ,
     20#endif
     21  ckR8("FFTWServer - ", true, false),
     22  ckR4("FFTWServer - ", true, false),
     23 _preserve_input(preserve_input) 
    1924{
    2025}
    … …  
    219224}
    220225
     226
     227#ifdef ALSO_FFTW_FLOAT_EXTSOP
     228 /* ---------------------------------------------------------------------------
     229   We define here single precision (float) version of FFTWServr methods
     230   Needs the libfftw3f.a , in addition to libfftw3.a   
     231   --------------------------------------------------------------------------- */
     232/* --Methode-- */
     233void
     234FFTWServer::FFTForward(TArray< complex<r_4> > & in, TArray< complex<r_4> > & out)
     235{
     236  int rank = ckR4.CheckResize(in, out);
     237  if (rank == 1) { // One dimensional transform
     238    fftwf_plan p = fftwf_plan_dft_1d(in.Size(), (fftwf_complex *)in.Data(),
     239                         (fftwf_complex *)out.Data(),
     240                         FFTW_FORWARD, FFTW_ESTIMATE);
     241    if (p == NULL)
     242      throw ParmError("FFTWServer::FFTForward( complex<r_4>, complex<r_4> ) Error creating fftwf_plan");
     243    fftwf_execute(p);
     244    fftwf_destroy_plan(p);
     245  }
     246  else {   // Multi dimensional
     247    if (in.NbDimensions() > MAXND_FFTW)
     248      throw ParmError("FFTWServer::FFTForward( complex<r_4>, complex<r_4> ) rank > MAXND_FFTW !");
     249    int sz[MAXND_FFTW];
     250    FillSizes4FFTW(in, sz);
     251    fftwf_plan p = fftwf_plan_dft(rank, sz,
     252                      (fftwf_complex *)in.Data(), (fftwf_complex *)out.Data(),
     253                      FFTW_FORWARD, FFTW_ESTIMATE);
     254    if (p == NULL)
     255      throw ParmError("FFTWServer::FFTForward( complex<r_4>, complex<r_4> ) Error creating fftwf_plan");
     256    fftwf_execute(p);
     257    fftwf_destroy_plan(p);
     258  } 
     259  if(this->getNormalize()) out=out/complex<r_4>((double)in.Size(),0.);
     260  return;
     261}
     262
     263/* --Methode-- */
     264void FFTWServer::FFTBackward(TArray< complex<r_4> > & in, TArray< complex<r_4> > & out)
     265{
     266  int rank = ckR4.CheckResize(in, out);
     267  if (rank == 1) { // One dimensional transform
     268    fftwf_plan p = fftwf_plan_dft_1d(in.Size(), (fftwf_complex *)in.Data(),
     269                                   (fftwf_complex *)out.Data(),
     270                                   FFTW_BACKWARD, FFTW_ESTIMATE);
     271    if (p == NULL)
     272      throw ParmError("FFTWServer::FFTBackward( complex<r_4>, complex<r_4> ) Error creating fftwf_plan");
     273    fftwf_execute(p);
     274    fftwf_destroy_plan(p);
     275  }
     276  else {   // Multi dimensional
     277    if (in.NbDimensions() > MAXND_FFTW)
     278      throw ParmError("FFTWServer::FFTBackward( complex<r_4>, complex<r_4> ) rank > MAXND_FFTW !");
     279    int sz[MAXND_FFTW];
     280    FillSizes4FFTW(in, sz);
     281    fftwf_plan p = fftwf_plan_dft(rank, sz,
     282                                (fftwf_complex *)in.Data(), (fftwf_complex *)out.Data(),
     283                                FFTW_BACKWARD, FFTW_ESTIMATE);
     284    if (p == NULL)
     285      throw ParmError("FFTWServer::FFTBackward( complex<r_4>, complex<r_4> ) Error creating fftwf_plan");
     286    fftwf_execute(p);
     287    fftwf_destroy_plan(p);
     288  } 
     289
     290  return;
     291}
     292
     293
     294void FFTWServer::FFTForward(TArray< r_4 > & in, TArray< complex<r_4> > & out)
     295{ 
     296  int rank = ckR4.CheckResize(in, out);
     297  if (rank == 1) { // One dimensional transform
     298    fftwf_plan p = fftwf_plan_dft_r2c_1d(in.Size(), in.Data(),
     299                                       (fftwf_complex *)out.Data(),
     300                                       FFTW_ESTIMATE); 
     301    if (p == NULL)
     302      throw ParmError("FFTWServer::FFTForward(r_4, complex<r_4> ) Error creating fftwf_plan");
     303    fftwf_execute(p);
     304    fftwf_destroy_plan(p);
     305
     306   
     307  }
     308  else {   // Multi dimensional
     309    if (in.NbDimensions() > MAXND_FFTW)
     310      throw ParmError("FFTWServer::FFTForward(r_4, complex<r_4> ) rank > MAXND_FFTW !");
     311    int sz[MAXND_FFTW];
     312    FillSizes4FFTW(in, sz);
     313    fftwf_plan p = fftwf_plan_dft_r2c(rank, sz, in.Data(),
     314                                      (fftwf_complex *)out.Data(),
     315                                      FFTW_ESTIMATE);
     316    if (p == NULL)
     317      throw ParmError("FFTWServer::FFTForward(r_4, complex<r_4> ) Error creating fftwf_plan");
     318    fftwf_execute(p);
     319    fftwf_destroy_plan(p);
     320  }
     321  if(this->getNormalize()) out=out/complex<r_4>((double)in.Size(),0.);
     322  return;
     323
     324}
     325
     326
     327
     328void FFTWServer::FFTBackward(TArray< complex<r_4> > & in, TArray< r_4 > & out,
     329                             bool usoutsz)
     330{
     331  // ATTENTION dans les TF (Complex->Reel), c'est la taille logique de la TF
     332  // qu'il faut indiquer lors de la creation des plans, cad taille tableau reel
     333  int rank = ckR4.CheckResize(in, out, usoutsz);
     334  bool share = (_preserve_input) ? false : true;
     335  TArray< complex<r_4> > inp(in, share);
     336  if (rank == 1) { // One dimensional transform
     337    fftwf_plan p = fftwf_plan_dft_c2r_1d(out.Size(), (fftwf_complex *)inp.Data(),
     338                             out.Data(),
     339                             FFTW_ESTIMATE); 
     340    if (p == NULL)
     341      throw ParmError("FFTWServer::FFTBackward(r_4, complex<r_4> ) Error creating fftwf_plan");
     342    fftwf_execute(p);
     343    fftwf_destroy_plan(p);
     344  }
     345  else {   // Multi dimensional
     346    if (in.NbDimensions() > MAXND_FFTW)
     347      throw ParmError("FFTWServer::FFTBackward(r_4, complex<r_4> ) rank > MAXND_FFTW !");
     348    int sz[MAXND_FFTW];
     349    FillSizes4FFTW(out, sz);
     350    fftwf_plan p = fftwf_plan_dft_c2r(rank, sz, (fftwf_complex *)inp.Data(),
     351                          out.Data(),
     352                          FFTW_ESTIMATE); 
     353    if (p == NULL)
     354      throw ParmError("FFTWServer::FFTBackward(r_4, complex<r_4> ) Error creating fftwf_plan");
     355    fftwf_execute(p);
     356    fftwf_destroy_plan(p);
     357  }
     358  return;
     359}
     360
     361
     362
     363/* --Methode-- */
     364void FFTWServer::ReShapetoReal(TArray< complex<r_4> > const & ina, TArray< r_4 > & outa,
     365                               bool usz)
     366{
     367  TVector< complex<r_4> > in(ina);
     368  TVector< r_4> out(outa);
     369  int n = in.NElts();
     370  r_4 thr = FFTArrayChecker<r_4>::ZeroThreshold();
     371  sa_size_t ncs;
     372  if (usz) {
     373    if ( (out.NElts() != 2*n-1) && (out.NElts() != 2*n-2) )
     374      throw SzMismatchError("FFTWServer::ReShapetoReal(..., true) - Wrong output array size ");
     375    ncs = out.NElts();
     376  }
     377  else {
     378    ncs = ( (in(n-1).imag() < -thr) || (in(n-1).imag() > thr) ) ?
     379                    ncs = 2*n-1 : ncs = 2*n-2;
     380
     381    if (out.NElts() != ncs) {
     382      cerr << " DEBUG-FFTWServer::ReShapetoReal() ncs= " << ncs
     383           << " out.NElts()= " << out.NElts() << endl;
     384      throw SzMismatchError("FFTWServer::ReShapetoReal() - Wrong output array size !");
     385    }
     386  }
     387  //  if (ncs == 2*n-2)  {
     388  out(0) = in(0).real();
     389  for(int k=1; k<n; k++) {
     390      out(k) = in(k).real();
     391      out(ncs-k) = in(k).imag();
     392  }
     393  if (ncs == 2*n-2)
     394    out(n-1) = in(n-1).real();
     395 
     396  //  for(int k=0; k<out.NElts(); k++) cout << "ReShapetoReal out " << k << " " << out(k) << endl;
     397}
     398
     399
     400/* --Methode-- */
     401void FFTWServer::ReShapetoCompl(TArray< r_4 > const & ina, TArray< complex<r_4> > & outa)
     402{
     403  TVector< r_4> in(ina);
     404  TVector< complex<r_4> > out(outa);
     405
     406  sa_size_t n = in.NElts();
     407  sa_size_t ncs = n/2+1;
     408  sa_size_t nc = (n%2 != 0) ? n/2+1 : n/2;
     409  if (out.NElts() != ncs)
     410    throw SzMismatchError("FFTWServer::ReShapetoCompl() - Wrong output array size !");
     411
     412  out(0) = complex<r_4> (in(0),0.);
     413  for(int k=1; k<ncs; k++)
     414    out(k) = complex<r_4>(in(k), in(n-k));
     415  if (n%2 == 0) out(ncs-1) = complex<r_4>(in(n/2), 0.);
     416
     417}
     418
     419#endif

  • trunk/SophyaExt/IFFTW/fftwserver.h

    r3000 r3359  
    2424  virtual FFTServerInterface * Clone();
    2525
    26   // Transforme unidimensionnelle , N-dimensionnel
     26  // Transforme unidimensionnelle , N-dimensionnel  (double precision - r_8)
    2727  virtual void FFTForward(TArray< complex<r_8> > & in, TArray< complex<r_8> > & out);
    2828  virtual void FFTBackward(TArray< complex<r_8> > & in, TArray< complex<r_8> > & out);
    … …  
    3131                           bool usoutsz=false);
    3232
    33 // Methodes statiques pour reordonner les donnees en 1-D
     33// Methodes statiques pour reordonner les donnees en 1-D (double precision - r_8)
    3434  static void ReShapetoReal(TArray< complex<r_8> > const & in, TArray< r_8 > & out, bool usz);
    3535  static void ReShapetoCompl(TArray< r_8 > const & in, TArray< complex<r_8> > & out);
    3636
     37#ifdef ALSO_FFTW_FLOAT_EXTSOP
     38 /* ---------------------------------------------------------------------------
     39   Declare and compile single precision (float) version of FFTWServr methods
     40   Needs the libfftw3f.a , in addition to libfftw3.a   
     41   --------------------------------------------------------------------------- */
     42
     43  // Transforme unidimensionnelle , N-dimensionnel  (single precision - r_4)
     44  virtual void FFTForward(TArray< complex<r_4> > & in, TArray< complex<r_4> > & out);
     45  virtual void FFTBackward(TArray< complex<r_4> > & in, TArray< complex<r_4> > & out);
     46  virtual void FFTForward(TArray< r_4 > & in, TArray< complex<r_4> > & out);
     47  virtual void FFTBackward(TArray< complex<r_4> > & in, TArray< r_4 > & out,
     48                           bool usoutsz=false);
     49
     50// Methodes statiques pour reordonner les donnees en 1-D (single precision - r_4)
     51  static void ReShapetoReal(TArray< complex<r_4> > const & in, TArray< r_4 > & out, bool usz);
     52  static void ReShapetoCompl(TArray< r_4 > const & in, TArray< complex<r_4> > & out);
     53 
     54#endif
     55
    3756 protected:
    38 #ifdef FFTW_V2_EXTSOP
     57
     58#ifndef FFTW_V2_EXTSOP   
     59/*  For FFTW >= V3  */
     60 FFTArrayChecker<r_8> ckR8;
     61 FFTArrayChecker<r_4> ckR4;
     62 bool _preserve_input; // if true, input arrays not overwritten
     63
     64#else
     65/*  For FFTW  V2  */
    3966 FFTWServerPlan * _p1df;
    4067 FFTWServerPlan * _p1db;
    … …  
    4673 FFTWServerPlan * _pndrf;
    4774 FFTWServerPlan * _pndrb;
    48 #endif
     75
    4976 FFTArrayChecker<r_8> ckR8;
    5077 bool _preserve_input; // if true, input arrays not overwritten
     78
     79#endif
    5180};
    5281
Note: See TracChangeset for help on using the changeset viewer.