source: Sophya/trunk/SophyaLib/NTools/fftservintf.cc@ 1652

#include "fftservintf.h"


/*!
  \class SOPHYA::FFTServerInterface
  \ingroup NTools
  Defines the interface for FFT (Fast Fourier Transform) operations.
  Definitions : 
    - Sampling period \b T
    - Sampling frequency \b fs=1/T
    - Total number of samples \b N
    - Frequency step in Fourier space \b =fs/N=1/(N*T)
    - Component frequencies
        - k=0      ->  0
        - k=1      ->  1/(N*T)
        - k        ->  k/(N*T)
        - k=N/2    ->  1/(2*T)   (Nyquist frequency)
        - k>N/2    ->  k/(N*T)   (or negative frequency -(N-k)/(N*T))

  For a sampling period T=1, the computed Fourier components correspond to :
  \verbatim
  0  1/N  2/N  ... 1/2  1/2+1/N  1/2+2/N ... 1-2/N  1-1/N
  0  1/N  2/N  ... 1/2                   ...  -2/N   -1/N
  \endverbatim

  For complex one-dimensional transforms:
  \f[
  out(i) = F_{norm} \Sigma_{j} \ e^{-2 \pi \sqrt{-1} \ i \  j} \ {\rm (forward)}
  \f]
  \f[
  out(i) = F_{norm} \Sigma_{j} \ e^{2 \pi \sqrt{-1} \ i \  j} \ {\rm (backward)}
  \f]
  i,j= 0..N-1 , where N is the input or the output array size.

  For complex multi-dimensional transforms:
  \f[
  out(i1,i2,...,id) = F_{norm} \Sigma_{j1} \Sigma_{j2} ... \Sigma_{jd} \ 
  e^{-2 \pi \sqrt{-1} \ i1 \ j1} ... e^{-2 \pi \sqrt{-1} \ id \ jd} \ {\rm (forward)}
  \f]
  \f[
  out(i1,i2,...,id) = F_{norm} \Sigma_{j1} \Sigma_{j2} ... \Sigma_{jd} \ 
  e^{2 \pi \sqrt{-1} \ i1 \ j1} ... e^{2 \pi \sqrt{-1} \ id \ jd} \ {\rm (backward)}
  \f]

  For real forward transforms, the input array is real, and
  the output array complex, with Fourier components up to k=N/2.
  For real backward transforms, the input array is complex and
  the output array is real. 
*/

/* --Methode-- */
FFTServerInterface::FFTServerInterface(string info)
{
  _info = info;
  _fgnorm = true;
}

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

// ----------------- Transforme pour les double -------------------

/* --Methode-- */
//! Forward Fourier transform for double precision complex data 
/*!
  \param in : Input complex array
  \param out : Output complex array
 */
void FFTServerInterface::FFTForward(TArray< complex<r_8> > const &, TArray< complex<r_8> > &)
{
  throw NotAvailableOperation("FFTServer::FFTForward(TArray...) Unsupported operation !");
}

/* --Methode-- */
//! Backward (inverse) Fourier transform for double precision complex data 
/*!
  \param in : Input complex array
  \param out : Output complex array
 */
void FFTServerInterface::FFTBackward(TArray< complex<r_8> > const &, TArray< complex<r_8> > &)
{
  throw NotAvailableOperation("FFTServer::FFTBackward(TArray...) Unsupported operation !");
}

/* --Methode-- */
//! Forward Fourier transform for double precision real input data
/*!
  \param in : Input real array
  \param out : Output complex array
 */
void FFTServerInterface::FFTForward(TArray< r_8 > const &, TArray< complex<r_8> > &)
{
  throw NotAvailableOperation("FFTServer::FFTForward(TArray...) Unsupported operation !");
}

/* --Methode-- */
//! Backward (inverse) Fourier transform for double precision real output data 
/*!
  \param in : Input complex array
  \param out : Output real array
  \param usoutsz : if true, use the output array size for computing the inverse FFT.
 */
void FFTServerInterface::FFTBackward(TArray< complex<r_8> > const &, TArray< r_8 > &, bool)
{
  throw NotAvailableOperation("FFTServer::FFTBackward(TArray...) Unsupported operation !");
}


// ----------------- Transforme pour les float -------------------

/* --Methode-- */
//! Forward Fourier transform for complex data 
/*!
  \param in : Input complex array
  \param out : Output complex array
 */
void FFTServerInterface::FFTForward(TArray< complex<r_4> > const &, TArray< complex<r_4> > &)
{
  throw NotAvailableOperation("FFTServer::FFTForward(TArray r_4 ... ) Unsupported operation !");
}

/* --Methode-- */
//! Backward (inverse) Fourier transform for complex data 
/*!
  \param in : Input complex array
  \param out : Output complex array
 */
void FFTServerInterface::FFTBackward(TArray< complex<r_4> > const &, TArray< complex<r_4> > &)
{
  throw NotAvailableOperation("FFTServer::FFTBackward(TArray r_4 ... ) Unsupported operation !");
}

/* --Methode-- */
//! Forward Fourier transform for real input data
/*!
  \param in : Input real array
  \param out : Output complex array
 */
void FFTServerInterface::FFTForward(TArray< r_4 > const &, TArray< complex<r_4> > &)
{
  throw NotAvailableOperation("FFTServer::FFTForward(TArray r_4 ... ) Unsupported operation !");
}

/* --Methode-- */
//! Backward (inverse) Fourier transform for real output data 
/*!
  \param in : Input complex array
  \param out : Output real array
  \param usoutsz : if true, use the output array size for computing the inverse FFT.
 */
void FFTServerInterface::FFTBackward(TArray< complex<r_4> > const &, TArray< r_4 > &, bool)
{
  throw NotAvailableOperation("FFTServer::FFTBackward(TArray r_4 ... ) Unsupported operation !");
}


// ----------------- Transformation de normalisation pour les energies -------------------

/* --Methode-- */
//! Compute the transform to be applied to "fftx=FFT(x)" so that "x" and "FFT(x)" have the same energy
/*!
\return The factor to be applied to the FFT energy such that we get the same energy as for the x.
\verbatim
  fftx is computed by: FFTForward(x,fftx)
  Energy of x    : Ex    = sum{|x(i)|^2}      i=0,x.Size()-1
  Energy of fftx : Efftx = sum{|fftx(i)^2|}   i=0,fftx.Size()-1
        ( usually x.Size() != fftx.Size() )
  -------------------------------------------------------------------
  |  TransfEnergyFFT return A and B such that : Ex = A * Efftx + B  |
  |                  and Norm such that : Ex = Norm * Efftx         |
  -------------------------------------------------------------------
  To normalize the fftx vector apply : "fftx *= sqrt(Norm)"
\endverbatim
*/
r_8 FFTServerInterface::TransfEnergyFFT
     (TVector< complex<r_8> > const& x, TVector< complex<r_8> > const& fftx, r_8& A, r_8& B)
{
  B=0.;
  if(getNormalize()) A = x.Size(); else A = 1./x.Size();
  r_8 norm = A;
  return norm;
}

//! Compute the transform to be applied to "fftx=FFT(x)" so that "x" and "FFT(x)" have the same energy
r_8 FFTServerInterface::TransfEnergyFFT
     (TVector< complex<r_4> > const & x, TVector< complex<r_4> > const & fftx, r_8& A, r_8& B)
{
  B=0.;
  if(getNormalize()) A = x.Size(); else A = 1./x.Size();
  r_8 norm = A;
  return norm;
}

//! Compute the transform to be applied to "fftx=FFT(x)" so that "x" and "FFT(x)" have the same energy
r_8 FFTServerInterface::TransfEnergyFFT
     (TVector< r_8 > const & x, TVector< complex<r_8> > const & fftx, r_8& A, r_8& B)
{
  A= 2.;
  B=  - abs(fftx(0)*fftx(0));
  if(x.Size()%2 == 0) B -= abs(fftx(fftx.Size()-1)*fftx(fftx.Size()-1));
  if(getNormalize()) {A *= x.Size(); B *= x.Size();}
     else            {A /= x.Size(); B /= x.Size();}
  r_8 norm = 0.;
  for(int_4 i=0;i<fftx.Size();i++) norm += abs(fftx(i)*fftx(i));
  if(norm>0.) norm = (A*norm+B)/norm;
  return norm;
}

//! Compute the transform to be applied to "fftx=FFT(x)" so that "x" and "FFT(x)" have the same energy
r_8 FFTServerInterface::TransfEnergyFFT
     (TVector< r_4 > const & x, TVector< complex<r_4> > const & fftx, r_8& A, r_8& B)
{
  A= 2.;
  B=  - abs(fftx(0)*fftx(0));
  if(x.Size()%2 == 0) B -= abs(fftx(fftx.Size()-1)*fftx(fftx.Size()-1));
  if(getNormalize()) {A *= x.Size(); B *= x.Size();}
     else            {A /= x.Size(); B /= x.Size();}
  r_8 norm = 0.;
  for(int_4 i=0;i<fftx.Size();i++) norm += abs(fftx(i)*fftx(i));
  if(norm>0.) norm = (A*norm+B)/norm;
  return norm;
}


/* --Methode-- */
/*!
  \class SOPHYA::FFTArrayChecker
  \ingroup NTools
  Service class for checking array size and resizing output arrays,
  to be used by FFTServer classes
*/

template <class T>
FFTArrayChecker<T>::FFTArrayChecker(string msg, bool checkpack, bool onedonly)
{
  _msg = msg + " FFTArrayChecker::";
  _checkpack = checkpack;
  _onedonly = onedonly;
}

/* --Methode-- */
template <class T>
FFTArrayChecker<T>::~FFTArrayChecker()
{
}

template <class T>
T FFTArrayChecker<T>::ZeroThreshold()
{
  return(0);
}

r_8 FFTArrayChecker< r_8 >::ZeroThreshold()
{
  return(1.e-18);
}

r_4 FFTArrayChecker< r_4 >::ZeroThreshold()
{
  return(1.e-9);
}

/* --Methode-- */
template <class T>
int FFTArrayChecker<T>::CheckResize(TArray< complex<T> > const & in, TArray< complex<T> > & out)
{
  int k;
  string msg;
  if (in.Size() < 1) {
    msg = _msg + "CheckResize(complex in, complex out) - Unallocated input array !";
    throw(SzMismatchError(msg));
  }
  if (_checkpack) 
    if ( !in.IsPacked() ) {
      msg = _msg + "CheckResize(complex in, complex out) - Not packed input array !";
      throw(SzMismatchError(msg));
    }
  int ndg1 = 0;
  for(k=0; k<in.NbDimensions(); k++) 
    if (in.Size(k) > 1)  ndg1++;
  if (_onedonly) 
    if (ndg1 > 1) {
      msg = _msg + "CheckResize(complex in, complex out) - Only 1-D array accepted !";
      throw(SzMismatchError(msg));
    }
  out.ReSize(in);
  //  sa_size_t sz[BASEARRAY_MAXNDIMS];
  //  for(k=0; k<in.NbDimensions(); k++) 
  //    sz[k] = in.Size(k);
  //  out.ReSize(in.NbDimensions(), sz);

  return(ndg1);
}

/* --Methode-- */
template <class T>
int FFTArrayChecker<T>::CheckResize(TArray< T > const & in, TArray< complex<T> > & out)
{
  int k;
  string msg;
  if (in.Size() < 1) {
    msg = _msg + "CheckResize(real in, complex out) - Unallocated input array !";
    throw(SzMismatchError(msg));
  }
  if (_checkpack) 
    if ( !in.IsPacked() ) {
      msg = _msg + "CheckResize(real in, complex out) - Not packed input array !";
      throw(SzMismatchError(msg));
    }
  int ndg1 = 0;
  for(k=0; k<in.NbDimensions(); k++) 
    if (in.Size(k) > 1)  ndg1++;
  if (_onedonly) 
    if (ndg1 > 1) {
      msg = _msg + "CheckResize(real in, complex out) - Only 1-D array accepted !";
      throw(SzMismatchError(msg));
    }
  sa_size_t sz[BASEARRAY_MAXNDIMS];
  // 
  if (ndg1 > 1) {
    sz[0] = in.Size(0)/2+1; 
    for(k=1; k<in.NbDimensions(); k++) 
      sz[k] = in.Size(k); 
  }
  else {
    for(k=0; k<BASEARRAY_MAXNDIMS; k++)  sz[k] = 1; 
    sz[in.MaxSizeKA()] = in.Size(in.MaxSizeKA())/2+1;
    //    sz[k] = in.Size(k)/2+1; 
    //    sz[k] = (in.Size(k)%2 != 0) ? in.Size(k)/2+1 : in.Size(k)/2;
  }
  out.ReSize(in.NbDimensions(), sz);

  return(ndg1);
}

/* --Methode-- */
template <class T>
int FFTArrayChecker<T>::CheckResize(TArray< complex<T> > const & in, TArray< T > & out, 
                                    bool usoutsz)
{
  int k;
  string msg;
  if (in.Size() < 1) {
    msg = _msg + "CheckResize(complex in, real out) - Unallocated input array !";
    throw(SzMismatchError(msg));
  }
  if (_checkpack) 
    if ( !in.IsPacked() ) {
      msg = _msg + "CheckResize(complex in, real out) - Not packed input array !";
      throw(SzMismatchError(msg));
    }
  int ndg1 = 0;
  for(k=0; k<in.NbDimensions(); k++) 
    if (in.Size(k) > 1)  ndg1++;
  if (_onedonly) 
    if (ndg1 > 1) {
      msg = _msg + "CheckResize(complex in, real out) - Only 1-D array accepted !";
      throw(SzMismatchError(msg));
    }
  if (usoutsz) { // We have to use output array size 
    bool fgerr = false;
    if (ndg1 > 1) {
      if (in.Size(0) != out.Size(0)/2+1) fgerr = true;
    }      
    else {
      if (in.Size(in.MaxSizeKA()) != out.Size(in.MaxSizeKA())/2+1) fgerr = true;
    }
    if (fgerr) {
        msg = _msg + "CheckResize(complex in, real out) - Incompatible in-out sizes !";
        throw(SzMismatchError(msg));
    }
  }
  else {  // We have to resize the output array 
    sa_size_t sz[BASEARRAY_MAXNDIMS];
    if (ndg1 > 1) {
      sz[0] = 2*in.Size(0)-1;
      for(k=1; k<in.NbDimensions(); k++) 
        sz[k] = in.Size(k);
    //      sz[k] = in.Size(k)*2-1;
    }
    else {
      for(k=0; k<BASEARRAY_MAXNDIMS; k++)  sz[k] = 1; 
      T thr = ZeroThreshold();
      sa_size_t n = in.Size(in.MaxSizeKA());
      sa_size_t ncs = ( (in[n-1].imag() < -thr) || (in[n-1].imag() > thr) )
                      ? 2*n-1 : 2*n-2;
      sz[in.MaxSizeKA()] = ncs;
    }
  out.ReSize(in.NbDimensions(), sz);
  }

  return(ndg1);

}


#ifdef __CXX_PRAGMA_TEMPLATES__
#pragma define_template FFTArrayChecker<r_4>
#pragma define_template FFTArrayChecker<r_8>
#endif

#if defined(ANSI_TEMPLATES) || defined(GNU_TEMPLATES)
template class FFTArrayChecker<r_4>;
template class FFTArrayChecker<r_8>;
#endif