fft.cpp File Reference

#include <cmath>
#include <ossim/matrix/include.h>
#include <ossim/matrix/newmatap.h>
#include <ossim/matrix/newmatio.h>

Go to the source code of this file.

Macros
#define	WANT_MATH
#define	REPORT   {}

Functions
void	FFTI (const ColumnVector &U, const ColumnVector &V, ColumnVector &X, ColumnVector &Y)
void	RealFFT (const ColumnVector &U, ColumnVector &X, ColumnVector &Y)
void	RealFFTI (const ColumnVector &A, const ColumnVector &B, ColumnVector &U)
void	FFT (const ColumnVector &U, const ColumnVector &V, ColumnVector &X, ColumnVector &Y)
void	DCT_II (const ColumnVector &U, ColumnVector &V)
void	DCT_II_inverse (const ColumnVector &V, ColumnVector &U)
void	DST_II (const ColumnVector &U, ColumnVector &V)
void	DST_II_inverse (const ColumnVector &V, ColumnVector &U)
void	DCT_inverse (const ColumnVector &V, ColumnVector &U)
void	DCT (const ColumnVector &U, ColumnVector &V)
void	DST_inverse (const ColumnVector &V, ColumnVector &U)
void	DST (const ColumnVector &U, ColumnVector &V)
void	FFT2 (const Matrix &U, const Matrix &V, Matrix &X, Matrix &Y)
void	FFT2I (const Matrix &U, const Matrix &V, Matrix &X, Matrix &Y)

Macro Definition Documentation

REPORT

#define REPORT   {}

Definition at line 23 of file fft.cpp.

Referenced by DCT(), DCT_II(), DCT_II_inverse(), DCT_inverse(), DST(), DST_II(), DST_II_inverse(), DST_inverse(), FFT(), FFT2(), FFT2I(), FFTI(), and RealFFT().

WANT_MATH

#define WANT_MATH

Definition at line 6 of file fft.cpp.

Function Documentation

DCT()

void DCT	(	const ColumnVector &	U,
		ColumnVector &	V
	)

Definition at line 393 of file fft.cpp.

References DCT_inverse(), GeneralMatrix::Nrows(), and REPORT.

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

DCT_II()

void DCT_II	(	const ColumnVector &	U,
		ColumnVector &	V
	)

Definition at line 249 of file fft.cpp.

References A, n, GeneralMatrix::Nrows(), RealFFT(), REPORT, ColumnVector::ReSize(), GeneralMatrix::Store(), x, and y.

{
   // 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;
   }
}

DCT_II_inverse()

void DCT_II_inverse	(	const ColumnVector &	V,
		ColumnVector &	U
	)

Definition at line 277 of file fft.cpp.

References n, GeneralMatrix::Nrows(), REPORT, GeneralMatrix::Store(), x, and y.

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

DCT_inverse()

void DCT_inverse	(	const ColumnVector &	V,
		ColumnVector &	U
	)

Definition at line 359 of file fft.cpp.

References A, ColumnVector::CleanUp(), n, GeneralMatrix::Nrows(), RealFFTI(), REPORT, ColumnVector::ReSize(), GeneralMatrix::Store(), x, and y.

Referenced by DCT().

{
   // 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;
   }
}

DST()

void DST	(	const ColumnVector &	U,
		ColumnVector &	V
	)

Definition at line 430 of file fft.cpp.

References DST_inverse(), GeneralMatrix::Nrows(), and REPORT.

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

DST_II()

void DST_II	(	const ColumnVector &	U,
		ColumnVector &	V
	)

Definition at line 304 of file fft.cpp.

References A, n, GeneralMatrix::Nrows(), RealFFT(), REPORT, ColumnVector::ReSize(), GeneralMatrix::Store(), x, and y.

{
   // 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;
   }
}

DST_II_inverse()

void DST_II_inverse	(	const ColumnVector &	V,
		ColumnVector &	U
	)

Definition at line 332 of file fft.cpp.

References n, GeneralMatrix::Nrows(), REPORT, GeneralMatrix::Store(), x, and y.

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

DST_inverse()

void DST_inverse	(	const ColumnVector &	V,
		ColumnVector &	U
	)

Definition at line 402 of file fft.cpp.

References A, ColumnVector::CleanUp(), n, GeneralMatrix::Nrows(), RealFFTI(), REPORT, ColumnVector::ReSize(), GeneralMatrix::Store(), x, and y.

Referenced by DST().

{
   // 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;
   }
}

FFT()

void FFT	(	const ColumnVector &	U,
		const ColumnVector &	V,
		ColumnVector &	X,
		ColumnVector &	Y
	)

Definition at line 197 of file fft.cpp.

References FFT_Controller::ar_1d_ft(), FFT_Controller::CanFactor(), n, GeneralMatrix::Nrows(), FFT_Controller::OnlyOldFFT, REPORT, and GeneralMatrix::Store().

Referenced by FFT2(), and FFTI().

{
   // from Carl de Boor (1980), Siam J Sci Stat Comput, 1 173-8
   // but first try Sande and Gentleman
   Tracer trace("FFT");
   REPORT
   const int n = U.Nrows();                     // length of arrays
   if (n != V.Nrows() || n == 0)
      Throw(ProgramException("Vector lengths unequal or zero", U, V));
   if (n == 1) { REPORT X = U; Y = V; return; }

   // see if we can use the newfft routine
   if (!FFT_Controller::OnlyOldFFT && FFT_Controller::CanFactor(n))
   {
      REPORT
      X = U; Y = V;
      if ( FFT_Controller::ar_1d_ft(n,X.Store(),Y.Store()) ) return;
   }

   ColumnVector B = V;
   ColumnVector A = U;
   X.ReSize(n); Y.ReSize(n);
   const int nextmx = 8;
   int prime[8] = { 2,3,5,7,11,13,17,19 };
   int after = 1; int before = n; int next = 0; bool inzee = true;
   int now = 0; int b1;             // initialised to keep gnu happy

   do
   {
      for (;;)
      {
     if (next < nextmx) { REPORT now = prime[next]; }
     b1 = before / now;  if (b1 * now == before) { REPORT break; }
     next++; now += 2;
      }
      before = b1;

      if (inzee) { REPORT fftstep(A, B, X, Y, after, now, before); }
      else { REPORT fftstep(X, Y, A, B, after, now, before); }

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

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

FFT2()

void FFT2	(	const Matrix &	U,
		const Matrix &	V,
		Matrix &	X,
		Matrix &	Y
	)

Definition at line 440 of file fft.cpp.

References BaseMatrix::Column(), FFT(), n, GeneralMatrix::Ncols(), GeneralMatrix::Nrows(), REPORT, BaseMatrix::Row(), and BaseMatrix::t().

Referenced by FFT2I(), and ossimFftFilter::runFft().

{
   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;
   }
}

FFT2I()

void FFT2I	(	const Matrix &	U,
		const Matrix &	V,
		Matrix &	X,
		Matrix &	Y
	)

Definition at line 461 of file fft.cpp.

References FFT2(), n, GeneralMatrix::Ncols(), GeneralMatrix::Nrows(), and REPORT.

Referenced by ossimFftFilter::runFft().

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

FFTI()

void FFTI	(	const ColumnVector &	U,
		const ColumnVector &	V,
		ColumnVector &	X,
		ColumnVector &	Y
	)

Definition at line 116 of file fft.cpp.

References FFT(), n, GeneralMatrix::Nrows(), and REPORT.

{
   // Inverse transform
   Tracer trace("FFTI");
   REPORT
   FFT(U,-V,X,Y);
   const Real n = X.Nrows(); X /= n; Y /= (-n);
}

RealFFT()

void RealFFT	(	const ColumnVector &	U,
		ColumnVector &	X,
		ColumnVector &	Y
	)

Definition at line 126 of file fft.cpp.

References A, n, GeneralMatrix::Nrows(), and REPORT.

Referenced by DCT_II(), and DST_II().

{
   // Fourier transform of a real series
   Tracer trace("RealFFT");
   REPORT
   const int n = U.Nrows();                     // length of arrays
   const int n2 = n / 2;
   if (n != 2 * n2)
      Throw(ProgramException("Vector length not multiple of 2", U));
   ColumnVector A(n2), B(n2);
   Real* a = A.Store(); Real* b = B.Store(); Real* u = U.Store(); int i = n2;
   while (i--) { *a++ = *u++; *b++ = *u++; }
   FFT(A,B,A,B);
   int n21 = n2 + 1;
   X.ReSize(n21); Y.ReSize(n21);
   i = n2 - 1;
   a = A.Store(); b = B.Store();              // first els of A and B
   Real* an = a + i; Real* bn = b + i;        // last els of A and B
   Real* x = X.Store(); Real* y = Y.Store();  // first els of X and Y
   Real* xn = x + n2; Real* yn = y + n2;      // last els of X and Y

   *x++ = *a + *b; *y++ = 0.0;                // first complex element
   *xn-- = *a++ - *b++; *yn-- = 0.0;          // last complex element

   int j = -1; i = n2/2;
   while (i--)
   {
      Real c,s; cossin(j--,n,c,s);
      Real am = *a - *an; Real ap = *a++ + *an--;
      Real bm = *b - *bn; Real bp = *b++ + *bn--;
      Real samcbp = s * am + c * bp; Real sbpcam = s * bp - c * am;
      *x++  =  0.5 * ( ap + samcbp); *y++  =  0.5 * ( bm + sbpcam);
      *xn-- =  0.5 * ( ap - samcbp); *yn-- =  0.5 * (-bm + sbpcam);
   }
}

RealFFTI()

void RealFFTI	(	const ColumnVector &	A,
		const ColumnVector &	B,
		ColumnVector &	U
	)

Definition at line 162 of file fft.cpp.

Referenced by DCT_inverse(), and DST_inverse().

{
   // inverse of a Fourier transform of a real series
   Tracer trace("RealFFTI");
   REPORT
   const int n21 = A.Nrows();                     // length of arrays
   if (n21 != B.Nrows() || n21 == 0)
      Throw(ProgramException("Vector lengths unequal or zero", A, B));
   const int n2 = n21 - 1;  const int n = 2 * n2;  int i = n2 - 1;

   ColumnVector X(n2), Y(n2);
   Real* a = A.Store(); Real* b = B.Store();  // first els of A and B
   Real* an = a + n2;   Real* bn = b + n2;    // last els of A and B
   Real* x = X.Store(); Real* y = Y.Store();  // first els of X and Y
   Real* xn = x + i;    Real* yn = y + i;     // last els of X and Y

   Real hn = 0.5 / n2;
   *x++  = hn * (*a + *an);  *y++  = - hn * (*a - *an);
   a++; an--; b++; bn--;
   int j = -1;  i = n2/2;
   while (i--)
   {
      Real c,s; cossin(j--,n,c,s);
      Real am = *a - *an; Real ap = *a++ + *an--;
      Real bm = *b - *bn; Real bp = *b++ + *bn--;
      Real samcbp = s * am - c * bp; Real sbpcam = s * bp + c * am;
      *x++  =  hn * ( ap + samcbp); *y++  =  - hn * ( bm + sbpcam);
      *xn-- =  hn * ( ap - samcbp); *yn-- =  - hn * (-bm + sbpcam);
   }
   FFT(X,Y,X,Y);             // have done inverting elsewhere
   U.ReSize(n); i = n2;
   x = X.Store(); y = Y.Store(); Real* u = U.Store();
   while (i--) { *u++ = *x++; *u++ = - *y++; }
}