synaps/linalg/FFT.h

Go to the documentation of this file.
00001 /*********************************************************************
00002 *      This file is part of the source code of SYNAPS kernel.        *
00003 *   Author(s): B. Mourrain, GALAAD, INRIA                            *
00004 **********************************************************************
00005 History: 
00006 $Id: FFT.h,v 1.1 2005/07/11 08:33:33 mourrain Exp $
00007 **********************************************************************/
00008 #ifndef SYNAPS_LINALG_FFT_H
00009 #define SYNAPS_LINALG_FFT_H
00010 #include <complex>
00011 #include <synaps/init.h>
00012 
00013 __BEGIN_NAMESPACE_SYNAPS
00014 //====================================================================
00018 //----------------------------------------------------------------------
00020 template<class T> inline void Root_Unity(std::complex<T> &c, int n) {
00021   c = std::complex<T>(cos(8.0*atan(1.0)/n), sin(8.0*atan(1.0)/n));
00022 }
00023 //----------------------------------------------------------------------
00028 inline unsigned int BitReverse(int k, int n)
00029 {
00030   int s, j = 0, m = k;
00031   double l =  std::log((double)n)/std::log((double)2);
00032   for (int q = 0; q < l; q++)
00033     {
00034       s = (int) std::floor(m/2.);
00035       j = 2*j + (m - 2*s);
00036       m = s;
00037     }
00038   return j;
00039 }
00040 //----------------------------------------------------------------------
00042 template <class T>
00043 inline void Permut(T *Tab, int N)
00044 {
00045   int j = 0;
00046   for (int k = 0; k < N; k++)
00047     {
00048       j = BitReverse(k, N);
00049       if (j > k)
00050         std::swap(Tab[j], Tab[k]);
00051     }
00052 }
00053 //-------------------------------------------------------------------
00061 template <class C, class U>
00062 void FFT(int n, const std::complex<U> & w,  std::complex<C>  *W, int inverse)
00063 {
00064   switch (n) {
00065   case 1:
00066     break;
00067   case 2: 
00068     {
00069       std::complex<C> t(W[0]);
00070       W[0] = t + W[1];
00071       W[1] = t - W[1];
00072     }
00073     break;
00074   case 4:
00075     {
00076        std::complex<C> x0(W[0]+W[1]), x1(W[0]-W[1]), y0(W[2]+W[3]), y1(W[2]-W[3]);
00077       if (inverse) {
00078         W[0] = (x0 + y0);
00079         W[1] = (x1 - std::complex<C>(-y1.imag(), y1.real()));
00080         W[2] = (x0 - y0);
00081         W[3] =(x1 + std::complex<C>(-y1.imag(), y1.real()));}
00082       else {
00083         W[0] = x0 + y0;
00084         W[1] = x1 + std::complex<C>(-y1.imag(), y1.real());
00085         W[2] = x0 - y0;
00086         W[3] = x1 - std::complex<C>(-y1.imag(), y1.real()); }
00087       }
00088     break;
00089   default:
00090   { 
00091     /* Since the array w has been permutated before the call to FFT,
00092        elements with even index in the initial array, are now the
00093        first n/2 elements and elements with odd index the last n/2.
00094        So the recursive process takes the initial array divided in two parts.
00095     */
00096 
00097     FFT(n/2,pow(w, 2), W, inverse);
00098     FFT(n/2, pow(w, 2), W+n/2,inverse);
00099     
00100     for (int k = 0; k< n/2; k++)
00101       {
00102         std::complex<C> z = pow(w,k)*W[k+n/2] ;
00103         W[k+n/2] = W[k] - z;
00104         W[k] = W[k] + z;
00105       }
00106   }
00107   }
00108 }
00109 //-----------------------------------------------------------------------------
00110 __END_NAMESPACE_SYNAPS
00111 #endif // SYNAPS_LINALG_FFT_H
00112 

SYNAPS DOCUMENTATION
logo