fft.cpp

非常好用的用C编写的矩阵类,可在不同编译器下编译使用.

      inzee = !inzee; after *= now;
   }
   while (before != 1);

   if (inzee) { REPORT A.release(); X = A; B.release(); Y = B; }
}

// Trigonometric transforms
// see Charles Van Loan (1992) "Computational frameworks for the fast
// Fourier transform" published by SIAM; section 4.4.

void DCT_II(const ColumnVector& U, ColumnVector& V)
{
   // Discrete cosine transform, type II, of a real series
   Tracer trace("DCT_II");
   REPORT
   const int n = U.Nrows();                     // length of arrays
   const int n2 = n / 2; const int n4 = n * 4;
   if (n != 2 * n2)
      Throw(ProgramException("Vector length not multiple of 2", U));
   ColumnVector A(n);
   Real* a = A.Store(); Real* b = a + n; Real* u = U.Store();
   int i = n2;
   while (i--) { *a++ = *u++; *(--b) = *u++; }
   ColumnVector X, Y;
   RealFFT(A, X, Y); A.cleanup();
   V.resize(n);
   Real* x = X.Store(); Real* y = Y.Store();
   Real* v = V.Store(); Real* w = v + n;
   *v = *x;
   int k = 0; i = n2;
   while (i--)
   {
      Real c, s; cossin(++k, n4, c, s);
      Real xi = *(++x); Real yi = *(++y);
      *(++v) = xi * c + yi * s; *(--w) = xi * s - yi * c;
   }
}

void DCT_II_inverse(const ColumnVector& V, ColumnVector& U)
{
   // Inverse of discrete cosine transform, type II
   Tracer trace("DCT_II_inverse");
   REPORT
   const int n = V.Nrows();                     // length of array
   const int n2 = n / 2; const int n4 = n * 4; const int n21 = n2 + 1;
   if (n != 2 * n2)
      Throw(ProgramException("Vector length not multiple of 2", V));
   ColumnVector X(n21), Y(n21);
   Real* x = X.Store(); Real* y = Y.Store();
   Real* v = V.Store(); Real* w = v + n;
   *x = *v; *y = 0.0;
   int i = n2; int k = 0;
   while (i--)
   {
      Real c, s; cossin(++k, n4, c, s);
      Real vi = *(++v); Real wi = *(--w);
      *(++x) = vi * c + wi * s; *(++y) = vi * s - wi * c;
   }
   ColumnVector A; RealFFTI(X, Y, A);
   X.cleanup(); Y.cleanup(); U.resize(n);
   Real* a = A.Store(); Real* b = a + n; Real* u = U.Store();
   i = n2;
   while (i--) { *u++ = *a++; *u++ = *(--b); }
}

void DST_II(const ColumnVector& U, ColumnVector& V)
{
   // Discrete sine transform, type II, of a real series
   Tracer trace("DST_II");
   REPORT
   const int n = U.Nrows();                     // length of arrays
   const int n2 = n / 2; const int n4 = n * 4;
   if (n != 2 * n2)
      Throw(ProgramException("Vector length not multiple of 2", U));
   ColumnVector A(n);
   Real* a = A.Store(); Real* b = a + n; Real* u = U.Store();
   int i = n2;
   while (i--) { *a++ = *u++; *(--b) = -(*u++); }
   ColumnVector X, Y;
   RealFFT(A, X, Y); A.cleanup();
   V.resize(n);
   Real* x = X.Store(); Real* y = Y.Store();
   Real* v = V.Store(); Real* w = v + n;
   *(--w) = *x;
   int k = 0; i = n2;
   while (i--)
   {
      Real c, s; cossin(++k, n4, c, s);
      Real xi = *(++x); Real yi = *(++y);
      *v++ = xi * s - yi * c; *(--w) = xi * c + yi * s;
   }
}

void DST_II_inverse(const ColumnVector& V, ColumnVector& U)
{
   // Inverse of discrete sine transform, type II
   Tracer trace("DST_II_inverse");
   REPORT
   const int n = V.Nrows();                     // length of array
   const int n2 = n / 2; const int n4 = n * 4; const int n21 = n2 + 1;
   if (n != 2 * n2)
      Throw(ProgramException("Vector length not multiple of 2", V));
   ColumnVector X(n21), Y(n21);
   Real* x = X.Store(); Real* y = Y.Store();
   Real* v = V.Store(); Real* w = v + n;
   *x = *(--w); *y = 0.0;
   int i = n2; int k = 0;
   while (i--)
   {
      Real c, s; cossin(++k, n4, c, s);
      Real vi = *v++; Real wi = *(--w);
      *(++x) = vi * s + wi * c; *(++y) = - vi * c + wi * s;
   }
   ColumnVector A; RealFFTI(X, Y, A);
   X.cleanup(); Y.cleanup(); U.resize(n);
   Real* a = A.Store(); Real* b = a + n; Real* u = U.Store();
   i = n2;
   while (i--) { *u++ = *a++; *u++ = -(*(--b)); }
}

void DCT_inverse(const ColumnVector& V, ColumnVector& U)
{
   // Inverse of discrete cosine transform, type I
   Tracer trace("DCT_inverse");
   REPORT
   const int n = V.Nrows()-1;                     // length of transform
   const int n2 = n / 2; const int n21 = n2 + 1;
   if (n != 2 * n2)
      Throw(ProgramException("Vector length not multiple of 2", V));
   ColumnVector X(n21), Y(n21);
   Real* x = X.Store(); Real* y = Y.Store(); Real* v = V.Store();
   Real vi = *v++; *x++ = vi; *y++ = 0.0;
   Real sum1 = vi / 2.0; Real sum2 = sum1; vi = *v++;
   int i = n2-1;
   while (i--)
   {
      Real vi2 = *v++; sum1 += vi2 + vi; sum2 += vi2 - vi;
      *x++ = vi2; vi2 = *v++; *y++ = vi - vi2; vi = vi2;
   }
   sum1 += vi; sum2 -= vi;
   vi = *v; *x = vi; *y = 0.0; vi /= 2.0; sum1 += vi; sum2 += vi;
   ColumnVector A; RealFFTI(X, Y, A);
   X.cleanup(); Y.cleanup(); U.resize(n+1);
   Real* a = A.Store(); Real* b = a + n; Real* u = U.Store(); v = u + n;
   i = n2; int k = 0; *u++ = sum1 / n2; *v-- = sum2 / n2;
   while (i--)
   {
      Real s = sin(1.5707963267948966192 * (++k) / n2);
      Real ai = *(++a); Real bi = *(--b);
      Real bz = (ai - bi) / 4 / s; Real az = (ai + bi) / 2;
      *u++ = az - bz; *v-- = az + bz;
   }
}

void DCT(const ColumnVector& U, ColumnVector& V)
{
   // Discrete cosine transform, type I
   Tracer trace("DCT");
   REPORT
   DCT_inverse(U, V);
   V *= (V.Nrows()-1)/2;
}

void DST_inverse(const ColumnVector& V, ColumnVector& U)
{
   // Inverse of discrete sine transform, type I
   Tracer trace("DST_inverse");
   REPORT
   const int n = V.Nrows()-1;                     // length of transform
   const int n2 = n / 2; const int n21 = n2 + 1;
   if (n != 2 * n2)
      Throw(ProgramException("Vector length not multiple of 2", V));
   ColumnVector X(n21), Y(n21);
   Real* x = X.Store(); Real* y = Y.Store(); Real* v = V.Store();
   Real vi = *(++v); *x++ = 2 * vi; *y++ = 0.0;
   int i = n2-1;
   while (i--) { *y++ = *(++v); Real vi2 = *(++v); *x++ = vi2 - vi; vi = vi2; }
   *x = -2 * vi; *y = 0.0;
   ColumnVector A; RealFFTI(X, Y, A);
   X.cleanup(); Y.cleanup(); U.resize(n+1);
   Real* a = A.Store(); Real* b = a + n; Real* u = U.Store(); v = u + n;
   i = n2; int k = 0; *u++ = 0.0; *v-- = 0.0;
   while (i--)
   {
      Real s = sin(1.5707963267948966192 * (++k) / n2);
      Real ai = *(++a); Real bi = *(--b);
      Real az = (ai + bi) / 4 / s; Real bz = (ai - bi) / 2;
      *u++ = az - bz; *v-- = az + bz;
   }
}

void DST(const ColumnVector& U, ColumnVector& V)
{
   // Discrete sine transform, type I
   Tracer trace("DST");
   REPORT
   DST_inverse(U, V);
   V *= (V.Nrows()-1)/2;
}

// Two dimensional FFT
void FFT2(const Matrix& U, const Matrix& V, Matrix& X, Matrix& Y)
{
   Tracer trace("FFT2");
   REPORT
   int m = U.Nrows(); int n = U.Ncols();
   if (m != V.Nrows() || n != V.Ncols() || m == 0 || n == 0)
      Throw(ProgramException("Matrix dimensions unequal or zero", U, V));
   X = U; Y = V;
   int i; ColumnVector CVR; ColumnVector CVI;
   for (i = 1; i <= m; ++i)
   {
      FFT(X.Row(i).t(), Y.Row(i).t(), CVR, CVI);
      X.Row(i) = CVR.t(); Y.Row(i) = CVI.t();
   }
   for (i = 1; i <= n; ++i)
   {
      FFT(X.Column(i), Y.Column(i), CVR, CVI);
      X.Column(i) = CVR; Y.Column(i) = CVI;
   }
}

void FFT2I(const Matrix& U, const Matrix& V, Matrix& X, Matrix& Y)
{
   // Inverse transform
   Tracer trace("FFT2I");
   REPORT
   FFT2(U,-V,X,Y);
   const Real n = X.Nrows() * X.Ncols(); X /= n; Y /= (-n);
}


#ifdef use_namespace
}
#endif


///@}