Changeset 717 in Sophya for trunk/SophyaLib/NTools/fftpserver.cc

Timestamp:

Feb 5, 2000, 6:22:55 PM (26 years ago)

Author:

ansari

Message:

Introduction FFTMayer, debug de FFTPack - Reza 5/2/2000

File:

• trunk/SophyaLib/NTools/fftpserver.cc (modified) (13 diffs)

trunk/SophyaLib/NTools/fftpserver.cc

@@ -26 +26 @@
 */
 
-/*! \fn virtual void FFTServer::fftf(int l, float* inout)
+/*! \fn virtual void FFTServer::fftf(int l, r_4* inout)
   \param l length of array
   \param inout input array /output forward FFT (original array destroyed)
 */
-/*! \fn virtual void FFTServer::fftb(int l, float* inout)
+/*! \fn virtual void FFTServer::fftb(int l, r_4* inout)
   \param l length of array
   \param inout input array /output backward FFT (original array destroyed)
 */
-/*! \fn virtual void FFTServer::fftf(int l, double* inout)
+/*! \fn virtual void FFTServer::fftf(int l, r_8* inout)
   \param l length of array
   \param inout input array /output forward FFT (original array destroyed)
   \param inout input/output array (original array destroyed)
 */
-/*! \fn virtual void FFTServer::fftb(int l, double* inout)
+/*! \fn virtual void FFTServer::fftb(int l, r_8* inout)
   \param l length of array
   \param inout input array /output backward FFT(original array destroyed)
 */
-/*!\fn  virtual void FFTServer::fftf(int l, complex<float>* inout)
+/*!\fn  virtual void FFTServer::fftf(int l, complex<r_4>* inout)
   \param l length of array
   \param inout input array /output forward FFT (original array destroyed)
 */
-/*! \fn virtual void FFTServer::fftb(int l, complex<float>* inout)
+/*! \fn virtual void FFTServer::fftb(int l, complex<r_4>* inout)
   \param l length of array
   \param inout input array /output backward FFT (original array destroyed)
 */
-/*! \fn virtual void FFTServer::fftf(int l, complex<double>* inout)
+/*! \fn virtual void FFTServer::fftf(int l, complex<r_8>* inout)
   \param l length of array
   \param inout input array /output forward FFT (original array destroyed)
 */
-/*! \fn virtual void FFTServer::fftb(int l, complex<double>* inout)
+/*! \fn virtual void FFTServer::fftb(int l, complex<r_8>* inout)
   \param l length of array
   \param inout input array /output backward FFT(original array destroyed)
@@ -69 +69 @@
 
 FFTPackServer::FFTPackServer()
-   : FFTServerInterface("FFTServer using FFTPack (C-version) package")
+   : FFTServerInterface("FFTPackServer using extended FFTPack (C-version) package")
 
 {
@@ -92 +92 @@
 }
 
+
 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()));
 }
 
@@ -104 +106 @@
 }
 
+
+
 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()));
 }
 
@@ -121 +126 @@
   fftf(inout.NElts(), inout.Data());
   ReShapetoCompl(inout, out);
+  if (getNormalize()) out *= (1./(r_4)(in.NElts()));
 }
 
@@ -128 +134 @@
   fftb(out.NElts(), out.Data());
 }
+
 
 void FFTPackServer::FFTForward(TVector< r_8 > const & in, TVector< complex<r_8> > & out)
@@ -134 +141 @@
   fftf(inout.NElts(), inout.Data());
   ReShapetoCompl(inout, out);
+  if (getNormalize()) out *= (1./(r_8)(in.NElts()));
 }
 
@@ -142 +150 @@
 }
 
+ 
 
 void FFTPackServer::checkint_rfft(int l)
@@ -149 +158 @@
   if (ws_rfft) delete[] ws_rfft;  //a transform
   sz_rfft = l;
-  ws_rfft = new float[2*l+15];
+  ws_rfft = new r_4[2*l+15];
   rffti_(&l, ws_rfft);
 }
@@ -159 +168 @@
   if (ws_cfft) delete[] ws_cfft;
   sz_cfft = l;
-  ws_cfft = new float[4*l+15];
+  ws_cfft = new r_4[4*l+15];
   cffti_(&l, ws_cfft);
 }
@@ -169 +178 @@
   if (ws_dfft) delete[] ws_dfft;
   sz_dfft = l;
-  ws_dfft = new double[2*l+15];
+  ws_dfft = new r_8[2*l+15];
   dffti_(&l, ws_dfft);
 }
@@ -179 +188 @@
   if (ws_cdfft) delete[] ws_cdfft;
   sz_cdfft = l;
-  ws_cdfft = new double[4*l+15];
+  ws_cdfft = new r_8[4*l+15];
   cdffti_(&l, ws_cdfft);
 }
@@ -186 +195 @@
    return inverse transformations */
 
-void FFTPackServer::fftf(int l, float* inout)
+void FFTPackServer::fftf(int 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 l, double* inout)
+  //  for (int k= 2;k<=(l+1)/2;k++) inout[2*k-2]=-inout[2*k-2];
+}
+
+void FFTPackServer::fftf(int 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 l, complex<float>* inout)
+  //  for (int k= 2;k<=(l+1)/2;k++) inout[2*k-2]=-inout[2*k-2];
+}
+
+void FFTPackServer::fftf(int l, complex<r_4>* inout)
 {
   checkint_cfft(l);
-  cfftf_(&l, (float *)(inout), ws_cfft);
-
-  /*
-  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, (r_4 *)(inout), ws_cfft);
+}
+
+void FFTPackServer::fftf(int l, complex<r_8>* inout)
+{
+  checkint_cdfft(l);
+  cdfftf_(&l, (r_8*)(inout), ws_cdfft);
+}
+
+void FFTPackServer::fftb(int l, r_4* inout)
+{
+  checkint_rfft(l);
+  rfftb_(&l, inout, ws_rfft);
+}
+
+void FFTPackServer::fftb(int l, r_8* inout)
+{
+  checkint_dfft(l);
+  dfftb_(&l, inout, ws_dfft);
+}
+
+void FFTPackServer::fftb(int l, complex<r_4>* inout)
+{
+  checkint_cfft(l);
+  cfftb_(&l, (r_4 *)(inout), ws_cfft);
+}
+
+void FFTPackServer::fftb(int 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();
   }
-  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 FFTPackServer::fftf(int l, complex<double>* inout)
-{
-  checkint_cdfft(l);
-  cdfftf_(&l, (double*)(inout), ws_cdfft);
-  /*
-  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();
+  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();
   }
-  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 FFTPackServer::fftb(int l, float* inout)
-{
-  checkint_rfft(l);
-  rfftf_(&l, inout, ws_rfft);
-}
-
-void FFTPackServer::fftb(int l, double* inout)
-{
-  checkint_dfft(l);
-  dfftf_(&l, inout, ws_dfft);
-}
-
-void FFTPackServer::fftb(int l, complex<float>* inout)
-{
-  checkint_cfft(l);
-  cfftb_(&l, (float *)(inout), ws_cfft);
-  /*
-  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 FFTPackServer::fftb(int l, complex<double>* inout)
-{
-  checkint_cdfft(l);
-  cfftb_(&l, (float *)(inout), ws_cfft);
-  /*
-  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;
-  */
-}
-
+  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.);
+  int k;
+  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.);
+  int k;
+  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.);
+}