00001
00002
00003
00004
00005
00006
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
00092
00093
00094
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