cfftw.hh

Go to the documentation of this file.

// #Id: cfftw.hh,v 1.5 2001/05/01 19:20:37 zweck Exp #

// #####################################################################
//
//    cfftw.cc
//
//    Methods module for the class of objects used as
//    an interface to the FFTW toolkit.
//
//    Copyright:
//       Optical Fiber Communications Laboratory (OFCL)
//       Computer Science & Electrical Engineering Department (CSEE)
//       University of Maryland Baltimore County (UMBC)
//
// #####################################################################
//    Author: Ronald Holzloehner  <holzloehner@umbc.edu>
// #####################################################################


//#define _USING_FFTW3_


#ifdef _USING_FFTW3_

//## You should include the next line if you are using an old GNU C++ compiler 
//#define _OLD_CPP_COMPILER_

#ifndef _CFFTW3_HH_
#define _CFFTW3_HH_



#include <iostream.h>
#include <complex.h>
//#include <string.h>
#include <string>
   using namespace std;
#include <malloc.h>
#include <new>

//#include <fftw.h>
#include <fftw3.h>

typedef complex<double> cplx; 


//## These lines must be commented for the old GNU C++ compiler in malkhut
#ifndef _OLD_CPP_COMPILER_
#ifndef _OVERLOADED_OPERATORS_
#define _OVERLOADED_OPERATORS_
//## Added by Ivan Lima to make the source compatible with g++
//## Most of the new code consists of overloading operators
using namespace std;

inline cplx operator+( cplx a, double b);
inline cplx operator+( cplx a, double b) 
{
   return a + cplx(b);
}  

inline cplx operator+( double b, cplx a);
inline cplx operator+( double b, cplx a) 
{
   return a + cplx(b);
}   

inline cplx operator*( cplx a, double b);
inline cplx operator*( cplx a, double b) 
{
   return a * cplx(b);
}   

inline cplx operator*( double b, cplx a);
inline cplx operator*( double b, cplx a) 
{
   return a * cplx(b);
}   

inline cplx operator/( cplx a, double b);
inline cplx operator/( cplx a, double b) 
{
   return a / cplx(b);
}   

inline double log (int a);
inline double log(int a)
{
   return log( (double) a );
}

inline double sqrt(int a);
inline double sqrt(int a)
{
   return sqrt( (double) a );
}
#endif   // _OVERLOADED_OPERATORS_
#endif   // _OLD_CPP_COMPILER_
//## End of the extra lines that must be added for the new GNU C++ compiler


#ifdef GLOBAL_NEW_MEMORY
extern volatile unsigned long GLOBAL_NEW_MEMORY_VAL;
#endif //GLOBAL_NEW_MEMORY

// overloaded "new" and "new[]" operator
// purpose: allocate memory aligned to 16-byte boundaries
inline void* operator new[](size_t t) throw (std::bad_alloc);
inline void* operator new[](size_t t) throw (std::bad_alloc)
{
  void *p = memalign(16, t);
  if(!p) {
    #ifdef _OLD_CPP_COMPILER_  
    throw bad_alloc();
    #endif  // _OLD_CPP_COMPILER_
  }    
  if(((unsigned int)(p))%16)
  {
    cerr << "\n\noperator new[]: bad pointer alignment (cfftw.hh)\n\n";
    exit(-1);
  }

#ifdef GLOBAL_NEW_MEMORY
  GLOBAL_NEW_MEMORY_VAL += (unsigned long)(t);
#endif //GLOBAL_NEW_MEMORY
#ifdef DEBUG_NEW
  cout << "\nAllocating " << t << " bytes! (new[])" << endl;
#ifdef GLOBAL_NEW_MEMORY
  cout << "Total allocated memory: " << GLOBAL_NEW_MEMORY_VAL << " bytes" << endl;
#endif //GLOBAL_NEW_MEMORY
#endif // DEBUG_NEW
  return(p);
}
inline void* operator new(size_t t) throw (std::bad_alloc);
inline void* operator new(size_t t) throw (std::bad_alloc)
{
  void *p = memalign(16, t);
  if(!p) {
    #ifdef _OLD_CPP_COMPILER_    
    throw bad_alloc();
    #endif  // _OLD_CPP_COMPILER_    
  }    
  if(((unsigned int)(p))%16)
  {
    cerr << "\n\noperator new: bad pointer alignment (ron_base.hh)\n\n";
    exit(-1);
  }
#ifdef GLOBAL_NEW_MEMORY
  GLOBAL_NEW_MEMORY_VAL += (unsigned long)(t);
#endif //GLOBAL_NEW_MEMORY
#ifdef DEBUG_NEW
  cout << "\nAllocating " << t << " bytes! (new)" << endl;
#ifdef GLOBAL_NEW_MEMORY
  cout << "Total allocated memory: " << GLOBAL_NEW_MEMORY_VAL << " bytes" << endl;
#endif //GLOBAL_NEW_MEMORY
#endif // DEBUG_NEW
  return(p);
}

inline void operator delete[](void *p) throw() ;
inline void operator delete[](void *p) throw()
{
#ifdef DEBUG_NEW
  cout << "Freeing memory (delete[])" << endl;
#endif // DEBUG_NEW
  free(p);
}
inline void operator delete(void *p) throw() ;
inline void operator delete(void *p) throw()
{
#ifdef DEBUG_NEW
  cout << "Freeing memory (delete)" << endl;
#endif // DEBUG_NEW
  free(p);
}



inline cplx* NewCplx(int N)
{

  cplx * fftw_data = (cplx*) fftw_malloc(sizeof(fftw_complex)*N); 

  if (!fftw_data)
    {
    cerr << "NewCplx: null pointer \n";
    exit(-1);
    }

  cplx * ptr = fftw_data;
  cplx * end_ptr = fftw_data + N;


  for( ; ptr < end_ptr; ptr++)
    *ptr = 0.0;

  return fftw_data;  
};

// ##############################


inline double* NewDouble(int N)
{
  double * start_ptr = (double*) memalign(16,N*sizeof(double));

  if((!start_ptr) || ((unsigned int)(start_ptr) % 16))
   {
     cerr << "NewDouble: null or unaligned pointer\n"; 
     exit(-1);
   }

  double * ptr = start_ptr;
  double * end_ptr = start_ptr + N;

  for( ; ptr < end_ptr; ptr++)
    *ptr = 0.0;

  return start_ptr;  

};

// ##############################


inline int* NewInt(int N)
{
  int * start_ptr = (int*) memalign(16,N*sizeof(int));

  if((!start_ptr) || ((unsigned int)(start_ptr) % 16))
   {
     cerr << "NewInt: null or unaligned pointer\n"; 
     exit(-1);
   }

  int * ptr = start_ptr;
  int * end_ptr = start_ptr + N;

  for( ; ptr < end_ptr; ptr++)
    *ptr = 0;

  return start_ptr;  

};


inline void FreeCplx(cplx * ptr) { fftw_free((fftw_complex*)ptr); };


inline void FreeDouble(double * ptr) { free(ptr); };


inline void FreeInt(int * ptr) { free(ptr); };


// ##### This class was created by Ronald Holzloehner


class cfftw
{
    int N, N_half;
    fftw_plan  forw, backw;
    cplx  *shift_work;
    int shift_len1, shift_len2;
    
  public:


    double T;


    double F;

 
    double  delta_f;
    double  delta_t;
    fftw_complex *fftw_in, *fftw_out; 
    cplx * work;


    cfftw(int N, double T);
    
    ~cfftw();

    void FFT(cplx *in, cplx *out);

    void IFFT(cplx *in, cplx *out);
    void fftshift(cplx *in);
    void ifftshift(cplx *in);

    // gives frequency at index i after FFT() ran
    // frequencies go like  0 1 2 3 -4 -3 -2 -1  (N=8)
    // input in [0, N-1], output in [-N/2, N/2-1]


    inline double freq(int i) {return delta_f 
                *( (i < N_half)? (i):(i-N) ); };

    // gives index of frequency with index k after FFT() ran
    // input in [-N/2, N/2-1], output in [0, N-1]


    inline int index(int k) {return ( (k >= 0)? (k):(k+N) ); };

    inline int give_N(void) {return N;};


};

//## These lines must be commented for the new GNU C++ compiler
// ##### Additional function used by the methods of cfftw
#ifdef _OLD_CPP_COMPILER_
template <class T> inline T max(T,T);
template <class T> inline T max(T a, T bb)
{
  return  ((a >= bb)? a:bb);
}
inline cplx max(cplx &a, cplx &bb);
#endif // _OLD_CPP_COMPILER_
//## End of the lines that must be included for the old GNU C++ compiler of malkhut


inline cplx max(cplx &a, cplx &bb)
{
  return (norm(a) >= norm(bb))? a:bb;
}


// ##### Example of use of cfftw
/*****************************  All these lines are commented
int main(void)     
{
cfftw *fft; 
int N;
cout "\n Enter N: " << endl;
cin >> N;
double T = 1.0;
fft = new cfftw(N, T, FFTW_MEASURE);
int i;
cplx *time = NewCplx(N);
cplx *freq = NewCplx(N);
for(i=0; i<N; i++)
  in[i] = i;  // put in something  
fft->FFT(time, freq);
fft->fftshift(freq);
// do something...
fft->ifftshift(freq);
fft->IFFT(freq, time);
}
*************************************/ // end of comment

#endif /* _CFFTW3_HH_ */ 


#endif /* _USING_FFTW3_ */


#ifndef _USING_FFTW3_


// #####################################################################
//
//    cfftw.cc
//
//    Methods module for the class of objects used as
//    an interface to the FFTW toolkit.
//
//    Copyright:
//       Optical Fiber Communications Laboratory (OFCL)
//       Computer Science & Electrical Engineering Department (CSEE)
//       University of Maryland Baltimore County (UMBC)
//
// #####################################################################
//    Author: Ronald Holzloehner  <holzloehner@umbc.edu>
// #####################################################################

//## You should include the next line if you are using an old GNU C++ compiler 
//#define _OLD_CPP_COMPILER_

#ifndef _CFFTW_HH_
#define _CFFTW_HH_

#include <iostream.h>
#include <complex.h>
//#include <string.h>
#include <string>
   using namespace std;
#include <malloc.h>
#include <new>
#include <fftw.h>


typedef complex<double> cplx; 


//## These lines must be commented for the old GNU C++ compiler in malkhut
#ifndef _OLD_CPP_COMPILER_
#ifndef _OVERLOADED_OPERATORS_
#define _OVERLOADED_OPERATORS_
//## Added by Ivan Lima to make the source compatible with g++
//## Most of the new code consists of overloading operators
using namespace std;

inline cplx operator+( cplx a, double b);
inline cplx operator+( cplx a, double b) 
{
   return a + cplx(b);
}  

inline cplx operator+( double b, cplx a);
inline cplx operator+( double b, cplx a) 
{
   return a + cplx(b);
}   

inline cplx operator*( cplx a, double b);
inline cplx operator*( cplx a, double b) 
{
   return a * cplx(b);
}   

inline cplx operator*( double b, cplx a);
inline cplx operator*( double b, cplx a) 
{
   return a * cplx(b);
}   

inline cplx operator/( cplx a, double b);
inline cplx operator/( cplx a, double b) 
{
   return a / cplx(b);
}   

inline double log (int a);
inline double log(int a)
{
   return log( (double) a );
}

inline double sqrt(int a);
inline double sqrt(int a)
{
   return sqrt( (double) a );
}
#endif   // _OVERLOADED_OPERATORS_
#endif   // _OLD_CPP_COMPILER_
//## End of the extra lines that must be added for the new GNU C++ compiler


#ifdef GLOBAL_NEW_MEMORY
extern volatile unsigned long GLOBAL_NEW_MEMORY_VAL;
#endif //GLOBAL_NEW_MEMORY

// overloaded "new" and "new[]" operator
// purpose: allocate memory aligned to 16-byte boundaries
inline void* operator new[](size_t t) throw (std::bad_alloc);
inline void* operator new[](size_t t) throw (std::bad_alloc)
{
  void *p = memalign(16, t);
  if(!p) {
    #ifdef _OLD_CPP_COMPILER_  
    throw bad_alloc();
    #endif  // _OLD_CPP_COMPILER_
  }    
  if(((unsigned int)(p))%16)
  {
    cerr << "\n\noperator new[]: bad pointer alignment (cfftw.hh)\n\n";
    exit(-1);
  }

#ifdef GLOBAL_NEW_MEMORY
  GLOBAL_NEW_MEMORY_VAL += (unsigned long)(t);
#endif //GLOBAL_NEW_MEMORY
#ifdef DEBUG_NEW
  cout << "\nAllocating " << t << " bytes! (new[])" << endl;
#ifdef GLOBAL_NEW_MEMORY
  cout << "Total allocated memory: " << GLOBAL_NEW_MEMORY_VAL << " bytes" << endl;
#endif //GLOBAL_NEW_MEMORY
#endif // DEBUG_NEW
  return(p);
}
inline void* operator new(size_t t) throw (std::bad_alloc);
inline void* operator new(size_t t) throw (std::bad_alloc)
{
  void *p = memalign(16, t);
  if(!p) {
    #ifdef _OLD_CPP_COMPILER_    
    throw bad_alloc();
    #endif  // _OLD_CPP_COMPILER_    
  }    
  if(((unsigned int)(p))%16)
  {
    cerr << "\n\noperator new: bad pointer alignment (ron_base.hh)\n\n";
    exit(-1);
  }
#ifdef GLOBAL_NEW_MEMORY
  GLOBAL_NEW_MEMORY_VAL += (unsigned long)(t);
#endif //GLOBAL_NEW_MEMORY
#ifdef DEBUG_NEW
  cout << "\nAllocating " << t << " bytes! (new)" << endl;
#ifdef GLOBAL_NEW_MEMORY
  cout << "Total allocated memory: " << GLOBAL_NEW_MEMORY_VAL << " bytes" << endl;
#endif //GLOBAL_NEW_MEMORY
#endif // DEBUG_NEW
  return(p);
}

inline void operator delete[](void *p) throw() ;
inline void operator delete[](void *p) throw()
{
#ifdef DEBUG_NEW
  cout << "Freeing memory (delete[])" << endl;
#endif // DEBUG_NEW
  free(p);
}
inline void operator delete(void *p) throw() ;
inline void operator delete(void *p) throw()
{
#ifdef DEBUG_NEW
  cout << "Freeing memory (delete)" << endl;
#endif // DEBUG_NEW
  free(p);
}







inline cplx* NewCplx(int N)
{
  cplx * start_ptr = (cplx*) memalign(16,N*sizeof(cplx));

  if (!start_ptr)
    {
    cerr << "NewCplx: null pointer \n";
    exit(-1);
    }

  if((unsigned int)(start_ptr) % 16)
   {
     cerr << "NewCplx: Unaligned pointer\n"; 
     exit(-1);
   }

  cplx * ptr = start_ptr;
  cplx * end_ptr = start_ptr + N;

  for( ; ptr < end_ptr; ptr++)
    *ptr = 0.0;

  return start_ptr;  
};

// ##############################


inline double* NewDouble(int N)
{
  double * start_ptr = (double*) memalign(16,N*sizeof(double));

  if((!start_ptr) || ((unsigned int)(start_ptr) % 16))
   {
     cerr << "NewDouble: null or unaligned pointer\n"; 
     exit(-1);
   }

  double * ptr = start_ptr;
  double * end_ptr = start_ptr + N;

  for( ; ptr < end_ptr; ptr++)
    *ptr = 0.0;

  return start_ptr;  

};

// ##############################


inline int* NewInt(int N)
{
  int * start_ptr = (int*) memalign(16,N*sizeof(int));

  if((!start_ptr) || ((unsigned int)(start_ptr) % 16))
   {
     cerr << "NewInt: null or unaligned pointer\n"; 
     exit(-1);
   }

  int * ptr = start_ptr;
  int * end_ptr = start_ptr + N;

  for( ; ptr < end_ptr; ptr++)
    *ptr = 0;

  return start_ptr;  

};


inline void FreeCplx(cplx * ptr) { free(ptr); };


inline void FreeDouble(double * ptr) { free(ptr); };


inline void FreeInt(int * ptr) { free(ptr); };


// ##### This class was created by Ronald Holzloehner
class cfftw
{
    int N, N_half;
    fftw_plan  forw, backw;
    cplx  *shift_work;
    int shift_len1, shift_len2;
    
  public:


    double T;


    double F;

 
    double  delta_f;
    double  delta_t;
    cplx *work;


    cfftw(int N, double T);
    
    ~cfftw();

    void FFT(cplx *in, cplx *out);

    void IFFT(cplx *in, cplx *out);
    void fftshift(cplx *in);
    void ifftshift(cplx *in);


    void timewindow(cplx *in, double T_half, double T_001);

    // gives frequency at index i after FFT() ran
    // frequencies go like  0 1 2 3 -4 -3 -2 -1  (N=8)
    // input in [0, N-1], output in [-N/2, N/2-1]


    inline double freq(int i) {return delta_f 
                *( (i < N_half)? (i):(i-N) ); };

    // gives index of frequency with index k after FFT() ran
    // input in [-N/2, N/2-1], output in [0, N-1]


    inline int index(int k) {return ( (k >= 0)? (k):(k+N) ); };

    inline int give_N(void) {return N;};

    // return true if we are using FFTW-GEL, false if not

    inline bool is_GEL(void) {return(strstr(fftw_version, "GEL")!=0);};

};

//## These lines must be commented for the new GNU C++ compiler
// ##### Additional function used by the methods of cfftw
#ifdef _OLD_CPP_COMPILER_
template <class T> inline T max(T,T);
template <class T> inline T max(T a, T bb)
{
  return  ((a >= bb)? a:bb);
}
inline cplx max(cplx &a, cplx &bb);
#endif // _OLD_CPP_COMPILER_
//## End of the lines that must be included for the old GNU C++ compiler of malkhut


inline cplx max(cplx &a, cplx &bb)
{
  return (norm(a) >= norm(bb))? a:bb;
}


// ##### Example of use of cfftw
/*****************************  All these lines are commented
int main(void)     
{
cfftw *fft; 
int N;
cout "\n Enter N: " << endl;
cin >> N;
double T = 1.0;
fft = new cfftw(N, T, FFTW_ESTIMATE);
int i;
cplx *time = new cplx[N];
cplx *freq = new cplx[N];
for(i=0; i<N; i++)
  in[i] = i;  // put in something  
fft->FFT(time, freq);
fft->fftshift(freq);
// do something...
fft->ifftshift(freq);
fft->IFFT(freq, time);
}
*************************************/ // end of comment

#endif /* _CFFTW_HH_ */ 

#endif /* ifndef _USING_FFTW3_ */