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

#include "sopnamsp.h"
#include "fftservintf.h"

//// VOIR GRAND BLABLA EXPLICATIF A LA FIN DU FICHIER

/*!
  \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 !");
}

/* --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);
}

DECL_TEMP_SPEC  /* equivalent a template <> , pour SGI-CC en particulier */
r_8 FFTArrayChecker< r_8 >::ZeroThreshold()
{
  return(1.e-39);
}

DECL_TEMP_SPEC  /* equivalent a template <> , pour SGI-CC en particulier */
r_4 FFTArrayChecker< r_4 >::ZeroThreshold()
{
  return(1.e-19);
}

/* --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





/**********************************************************************

Memo uniquement destine aux programmeurs:    (cmv 03/05/04)
-- cf programme de tests explicatif: cmvtfft.cc

=====================================================================
=====================================================================
============== Transformees de Fourier et setNormalize ==============
=====================================================================
=====================================================================

- si setNormalize(true):  invTF{TF{S}} = S
- si setNormalize(false): invTF{TF{S}} = N * S

=====================================================================
=====================================================================
============== Transformees de Fourier de signaux REELS =============
=====================================================================
=====================================================================

-------
--- FFT d'un signal REEL S ayant un nombre pair d'elements N=2p
-------
  taille de la FFT: Nfft = N/2 + 1 = p + 1
  abscisses de la fft: | 0 | 1/N | 2/N | ..... | p/N=1/2 |
                         ^continu                 ^frequence de Nyquist

  ... Ex:  N=6 -> Nfft = 6/3+1 = 4

  le signal a N elements reels, la fft a Nfft elements complexes
    cad 2*Nfft reels = 2*(p+1) reels = 2p + 2 reels = N + 2 reels
    soit 2 reels en trop:
    ce sont les phases du continu et de la frequence de Nyquist

  relations:
    - si setNormalize(true)  : fac = N
         setNormalize(false) : fac = 1/N
    sum(i=0,N-1){S(i)^2}
                        = fac* [[ 2* sum(j=0,Nfft-1){|TF{S}(j)|^2} 
                                 - |TF{S}(0)|^2 - |TF{S}(Nfft-1)|^2 ]]
    (On ne compte pas deux fois le continu et la freq de Nyquist)
     

-------
--- FFT d'un signal REEL ayant un nombre impair d'elements N=2p+1
-------
  taille de la FFT: Nfft = N/2 + 1 = p + 1
  abscisses de la fft: | 0 | 1/N | 2/N | ..... | p/N |
                         ^continu                 
                       (la frequence de Nyquist n'y est pas)

  ... Ex:  N=7 -> Nfft = 7/3+1 = 4

  le signal a N elements reels, la fft a Nfft elements complexes
    cad 2*Nfft reels = 2*(p+1) reels = 2p + 2 reels = N + 1 reels
    soit 1 reel en trop: c'est la phase du continu

  relations:
    - si setNormalize(true)  : fac = N
         setNormalize(false) : fac = 1/N
    sum(i=0,N-1){S(i)^2}
                        = fac* [[ 2* sum(j=0,Nfft-1){|TF{S}(j)|^2} 
                                 - |TF{S}(0)|^2 ]]
    (On ne compte pas deux fois le continu)


------------
--- FFT-BACK d'un signal F=TF{S} ayant un nombre d'elements Nfft
------------
  Sback = invTF{TF{S}}

  Remarque: Nfft a la meme valeur pour N=2p et N=2p+1
            donc Nfft conduit a 2 possibilites:
                 { N = 2*(Nfft-1)  signal back avec nombre pair d'elements
                 { N = 2*Nfft-1    signal back avec nombre impair d'elements

  Pour savoir quel est la longueur N du signal TF^(-1){F} on regarde
    si F(Nfft-1) est reel ou complexe
       (la frequence de Nyquist d'un signal reel est reelle)

    - Si F(Nfft-1) reel      cad  Im{F(Nfft-1)}=0:  N = 2*(Nfft-1)
    - Si F(Nfft-1) complexe  cad  Im{F(Nfft-1)}#0:  N = 2*Nfft-1

  Si setNormalize(true):  invTF{TF{S}} = S
     setNormalize(false): invTF{TF{S}} = N * S

=========================================================================
=========================================================================
============== Transformees de Fourier de signaux COMPLEXES =============
=========================================================================
=========================================================================

-------
--- FFT d'un signal COMPLEXE S ayant un nombre d'elements N
-------
  taille de la FFT: Nfft = N
  abscisses de la fft: | 0 | 1/N | 2/N | ..... | (N-1)/N |
                         ^continu

  Frequence de Nyquist:
    si N est pair: la frequence de Nyquist est l'absicce d'un des bins
      abscisses de TF{S}: Nfft = N = 2p
      | 0 | 1/N | 2/N | ... | (N/2)/N=p/N=0.5 | ... | (N-1)/N |
                                     ^frequence de Nyquist
    si N est impair: la frequence de Nyquist N'est PAS l'absicce d'un des bins
      abscisses de TF{S}: Nfft = N = 2p+1
      | 0 | 1/N | 2/N | ... | (N/2)/N=p/N | ((N+1)/2)/N=(p+1)/N  | ... | (N-1)/N |

  ... Ex: N = 2p =6  ->  Nfft = 2p = 6
      abscisses de TF{S}: | 0 | 1/6 | 2/6 | 3/6=0.5 | 4/6 | 5/6 |
  ... Ex: N = 2p+1 = 7   ->  Nfft = 2p+1 = 7
      abscisses de TF{S}: | 0 | 1/7 | 2/7 | 3/7 | 4/7 | 5/7 | 6/7 |

  relations:
    - si setNormalize(true)  : fac = N
         setNormalize(false) : fac = 1/N
    sum(i=0,N-1){S(i)^2} = fac* [[ sum(j=0,Nfft-1){|TF{S}(j)|^2} ]]

------------
--- FFT-BACK d'un signal F=TF{S} ayant un nombre d'elements Nfft
------------
  taille du signal: N = Nfft

  Si setNormalize(true):  invTF{TF{S}} = S
     setNormalize(false): invTF{TF{S}} = N * S

**********************************************************************/