Many many time ago, when FFTW isn't yet wroted...

dcomplex.cpp

//---------------------------------------------------------------------------
#ifndef DcomplexCPP
#define DcomplexCPP
//---------------------------------------------------------------------------
template <class T>
class dcomplex
{
    public:
	T re, im;

	dcomplex (const T& re_arg=0, const T& imag_arg=0)
	{
		re = re_arg;
		im = imag_arg;
	}

	dcomplex<T>& operator=  (const T&);
	dcomplex<T>& operator+= (const T&);
	dcomplex<T>& operator-= (const T&);
	dcomplex<T>& operator*= (const T&);
	dcomplex<T>& operator/= (const T&);

        dcomplex<T>& operator=  (const dcomplex<T>&);
        dcomplex<T>& operator+= (const dcomplex<T>&);
        dcomplex<T>& operator-= (const dcomplex<T>&);
        dcomplex<T>& operator*= (const dcomplex<T>&);
        dcomplex<T>& operator/= (const dcomplex<T>&);
};

//
// Temporarily turn off the warnings under the Borland compiler that
// say 'Functions containing ... cannot be inlined'
//
//#if defined(__BORLANDC__)
//#pragma option -w-inl
//#endif
//! @todo Check this

template <class T>
inline dcomplex<T>& dcomplex<T>::operator= (const T& rhs)
{
	re = rhs; im = 0.0;
	return *this;
}

template <class T>
inline dcomplex<T>& dcomplex<T>::operator+= (const T& rhs)
{
	re += rhs;
	return *this;
}

template <class T>
inline dcomplex<T>& dcomplex<T>::operator-= (const T& rhs) 
{
	re -= rhs;
	return *this;
}

template <class T>
inline dcomplex<T>& dcomplex<T>::operator*= (const T& rhs) 
{     
	im *= rhs; 
	re *= rhs;
	return *this;
}

template <class T>
inline dcomplex<T>& dcomplex<T>::operator/= (const T& rhs)
{
	re /= rhs;
	im /= rhs;
	return *this;
}

template <class T>
inline dcomplex<T>& dcomplex<T>::operator= (const dcomplex<T>& rhs)
{
        re = rhs.re;
        im = rhs.im;
        return *this;
}

template <class T>
inline dcomplex<T>& dcomplex<T>::operator+= (const dcomplex<T>& rhs)
{
        re += rhs.re;
        im += rhs.im;
        return *this;
}

template <class T>
inline dcomplex<T>& dcomplex<T>::operator-= (const dcomplex<T>& rhs)
{
        re -= rhs.re;
        im -= rhs.im;
        return *this;
}

template <class T>
inline dcomplex<T>& dcomplex<T>::operator*= (const dcomplex<T>& rhs)
{
        T tmp = re*rhs.re - im*rhs.im;
        im = im*rhs.re + re*rhs.im;
        re = tmp;
        return *this;
}

template <class T>
inline dcomplex<T>& dcomplex<T>::operator/= (const dcomplex<T>& rhs)
{
        T norm = rhs.re*rhs.re + rhs.im*rhs.im;
        T tmp = (re*rhs.re + im*rhs.im)/norm;
        im = (im*rhs.re - re*rhs.im)/norm;
        re = tmp;
        return *this;
}

template <class T>
inline dcomplex<T> operator+ (const dcomplex<T>& lhs, const dcomplex<T>& rhs)
{
	dcomplex<T> tmp = lhs;
	return tmp += rhs;
}

template <class T>
inline dcomplex<T> operator+ (const dcomplex<T>& lhs, const T& rhs)
{
	return dcomplex<T>(rhs+lhs.re, lhs.im);
}

template <class T>
inline dcomplex<T> operator+ (const T& lhs, const dcomplex<T>& rhs)
{
	return dcomplex<T>(lhs+rhs.re, rhs.im);
}

template <class T>
inline dcomplex<T> operator- (const dcomplex<T>& lhs, const dcomplex<T>& rhs)
{
	dcomplex<T> tmp = lhs;
	return tmp -= rhs;
}

template <class T>
inline dcomplex<T> operator- (const dcomplex<T>& lhs, const T& rhs)
{
	return dcomplex<T>(lhs.re-rhs, lhs.im);
}

template <class T>
inline dcomplex<T> operator- (const T& lhs, const dcomplex<T>& rhs)
{
	return dcomplex<T>(lhs-rhs.re, -rhs.im);
}

template <class T>
inline dcomplex<T> operator* (const dcomplex<T>& lhs, const dcomplex<T>& rhs)
{
	dcomplex<T> tmp = lhs;
	return tmp *= rhs;
}

template <class T>
inline dcomplex<T> operator* (const dcomplex<T>& lhs, const T& rhs)
{
	return dcomplex<T>(rhs*lhs.re, rhs*lhs.im);
}

template <class T>
inline dcomplex<T> operator* (const T& lhs, const dcomplex<T>& rhs)
{
	return dcomplex<T>(lhs*rhs.re, lhs*rhs.im);
}

template <class T>
inline dcomplex<T> operator/ (const dcomplex<T>& lhs, const dcomplex<T>& rhs)
{
	dcomplex<T> tmp = lhs;
	return tmp /= rhs;
}

template <class T>
inline dcomplex<T> operator/ (const dcomplex<T>& lhs, const T& rhs)
{
	return dcomplex<T>(lhs.re/rhs, lhs.im/rhs);
}

template <class T>
inline dcomplex<T> operator/ (const T& lhs, const dcomplex<T>& rhs)
{
	register T denom = rhs.re*rhs.re + rhs.im*rhs.im;
	return dcomplex<T>(lhs*rhs.re/denom,(-lhs*rhs.im)/denom);
}

template <class T>
inline dcomplex<T> operator+ (const dcomplex<T>& lhs)
{
	return lhs;
}

template <class T>
inline dcomplex<T> operator- (const dcomplex<T>& lhs)
{
	return dcomplex<T>(-lhs.re, -lhs.im);
}

template <class T>
inline bool operator== (const dcomplex<T>& lhs, const dcomplex<T>& rhs)
{
	return lhs.re == rhs.re && lhs.im == rhs.im;
}

template <class T>
inline bool operator== (const T& lhs, const dcomplex<T>& rhs)
{
	return lhs == rhs.re && rhs.im == 0;
}

template <class T>
inline bool operator== (const dcomplex<T>& lhs, const T& rhs)
{
	return lhs.re == rhs && lhs.im == 0;
}

template <class T>
inline bool operator!= (const dcomplex<T>& lhs, const dcomplex<T>& rhs)
{
	return lhs.re != rhs.re || lhs.im != rhs.im;
}

template <class T>
inline bool operator!= (const T& lhs, const dcomplex<T>& rhs)
{
	return lhs != rhs.re || rhs.im != 0;
} 

template <class T>
inline bool operator!= (const dcomplex<T>& lhs, const T& rhs)
{
	return lhs.re != rhs || lhs.im != 0;
}

template <class T>
inline T real (const dcomplex<T>& a)
{
	return a.re;
}

template <class T>
inline T imag (const dcomplex<T>& a)
{
	return a.im;
}

template <class T>
inline T norm (const dcomplex<T>& a)
{
	return a.re*a.re + a.im*a.im;
}

template <class T>
inline T abs (const dcomplex<T>& a)
{
	return sqrt(norm(a));
}

template <class T>
inline T arg (const dcomplex<T>& a)
{
	return a == dcomplex<T>(0,0) ? T(0) : atan2(a.im, a.re);
}

template <class T>
dcomplex<T> conj (const dcomplex<T>& a)
{
	return dcomplex<T>(a.re, -a.im);
}

template <class T>
inline dcomplex<T> polar (const T& r, const T& theta)
{
	return dcomplex<T>(r*cos(theta), r*sin(theta));
}

//
// transcendentals
//

//
// dcomplex<T> cosine of dcomplex<T> number a
//      cos (a) = cos u * cosh v - i * sin u * sinh v
//
template <class T>
inline dcomplex<T> cos (const dcomplex<T>& a)
{
	return dcomplex<T>(cos(a.re)*cosH(a.im), -sin(a.re)*sinH(a.im));
}

//
// dcomplex<T> hyperbolic cosine of dcomplex<T> number a
//      cosh (a) = cosh u * cosv + i * sinh u * sin v
//
template <class T>
inline dcomplex<T> cosh (const dcomplex<T>& a)
{
	return dcomplex<T>(cosH(a.re)*cos(a.im), sinH(a.re)*sin(a.im));
}

//
// dcomplex<T> exponential of  dcomplex<T> number a
//      exp (a) = exp(u) * (cos v + i * sin v)
//
template <class T>
inline dcomplex<T> exp (const dcomplex<T>& a)
{
	register T e = exp(a.re);
	return dcomplex<T>(e*cos(a.im), e*sin(a.im));
}

//
// dcomplex<T> natural log of dcomplex<T> number a
//      log(a) = log(r) + i * theta
//
template <class T>
inline dcomplex<T> log (const dcomplex<T>& a)
{
	return dcomplex<T>(log(abs(a)), arg(a));
}

//
// For all the power functions:
//
//   0**0 == 1
//   0**x == 0 for x != 0
//

//
// dcomplex<T> number a raised to an integer power n
//
// a**n = r**n * (cos(n theta) + i sin (n theta))
//
template <class T>
inline dcomplex<T> pow (const dcomplex<T>& a, int n)
{
	if( a == dcomplex<T>(0, 0) )
	{
		if( n == 0 ) 
			return dcomplex<T>(1, 0);
		else
			return dcomplex<T>(0, 0);
	}

	if( a.im == 0 )
	{
		if( a.re < 0 )
			return pow(a, dcomplex<T>(n, 0));
		else
			return dcomplex<T>(pow(a.re, n), 0);
	}

	register T r  = pow(T(abs(a)), T(n));
	register T th = T(n) * arg(a);

	return dcomplex<T>(r*cos(th), r*sin(th));
}

//
// dcomplex<T> number a raised to a real power s
//
// a**s = exp(s * log(a))
//
template <class T>
inline dcomplex<T> pow (const dcomplex<T>& a, T s)
{
	if( a == dcomplex<T>(0,0) )
	{
		if( s == T(0) )
			return dcomplex<T>(1, 0);
		else
			return dcomplex<T>(0, 0);
	}
	
	if( a.im == 0 )
	{
		if( a.re < 0 )
			return pow(a, dcomplex<T>(s,0));
		else
			return dcomplex<T>(pow(a.re,s), 0);
	}
	
	return exp(s*log(a));
}   

//
// real number s raised to a dcomplex<T> power a
//
//  s**a = exp(a * log (s))
//
template <class T>
inline dcomplex<T> pow (T s, const dcomplex<T>& a)
{
	if( s == T(0) )
	{
		if( a == dcomplex<T>(0, 0) )
			return dcomplex<T>(1, 0);
		else
			return dcomplex<T>(0, 0);
	}
	
	if( s < 0 )
		return pow(dcomplex<T>(s,0), a);

	if( a.im == 0 )
		return dcomplex<T>(pow(s, a.re), 0);

	return dcomplex<T>(exp(a*(T)log(s)));
}

//
// dcomplex<T> number a1 raised to a dcomplex<T> power a2
//
// a1**a2 = rho * (cos(phi) + i sin(phi))
//      rho = r1 **u2   *  exp (-v2* theta1)
//      phi = v2 * log(r1) + u2 * theta1
//
template <class T>
inline dcomplex<T> pow (const dcomplex<T>& a1, const dcomplex<T>& a2)
{
	if( a1 == dcomplex<T>(0, 0) )
	{
		if( a2 == dcomplex<T>(0, 0) )
			return dcomplex<T>(1, 0);
		else
			return dcomplex<T>(0, 0);
	}

	T r1  = abs(a1);
	T u2  = real(a2);
	T v2  = imag(a2);
	T th1 = arg(a1);
	T rho = pow(r1, u2)*exp(-v2*th1);
	T phi = v2*log(r1) + u2*th1;

	return dcomplex<T>(rho*cos(phi), rho*sin(phi));
}    
  
//
// dcomplex<T> sine of dcomplex<T> number a
//      sin (a) = sin u * cosh v + i * cos u * sinh v
//
template <class T>
inline dcomplex<T> sin (const dcomplex<T>& a)
{
	return dcomplex<T>(sin(a.re)*cosH(a.im), cos(a.re)*sinH(a.im));
}

//
// dcomplex<T> hyperbolic sine of dcomplex<T> number a
//      sinh (a) = sinh u cos v + i cosh u sin v
//
template <class T>
inline dcomplex<T> sinh (const dcomplex<T>& a)
{
	return dcomplex<T>(sinH(a.re)*cos(a.im), cosH(a.re)*sin(a.im));
}

//
// dcomplex<T> square root of dcomplex<T> number a
//      sqrt(a) = sqrt(r) * ( cos (theta/2) + i sin (theta/2) )
//
template <class T>
inline dcomplex<T> sqrt (const dcomplex<T>& a)
{
	register T r  = sqrt(abs(a));
	register T th = arg(a)/2.;
	return dcomplex<T>(r*cos(th), r*sin(th));
}

template <class T>
inline dcomplex<T> tan (const dcomplex<T>& a)
{
	return sin(a)/cos(a);
}

template <class T>
inline dcomplex<T> tanh (const dcomplex<T>& a)
{
	return sinh(a)/cosh(a);
}
//---------------------------------------------------------------------------
#endif
//---------------------------------------------------------------------------

fft.cpp

//---------------------------------------------------------------------------
#ifndef FftCPP
#define FftCPP
//---------------------------------------------------------------------------
#include "dcomplex.cpp"
#include <stdlib.h>
#include <math.h>
//---------------------------------------------------------------------------
/*
! 1D Fast Fourier Transform.
! Replaces XOK by its discrete Fourier transform,
! if IC is input as 1 or replaces XOK by its inverse discrete
! Fourier transform, if IC is input as -1. XOK is a complex array
! of length N. N  must be an integer power of 2. This is checked for.
!
! In: complex xok(n) - complex field or its spectrum
!     integer n - array dimension  - integer power of two
!     integer ic - direct (ic=+1) or inverse (ic=-1) transform
!
! Out: xok(n) - complex field or its spectrum
!      n -  unchanged
!      ic - unchanged
!
! Calls to external FUNCTIONS and SUBROUTINES:
*/
template <class ddouble>
int FFT1D(dcomplex <ddouble> *xok, int n, int ic);
//---------------------------------------------------------------------------
/*
! 2D FFT. Replaces AM by its 2D discrete Fourier transform,
! if IC is input as 1 or replaces AM by its inverse discrete
! Fourier transform, if IC is input as -1. AM is a complex array
! (N*N2). Both N and N2 must be an integer power of 2.
!
! NOTE  Internal static complex array XOK is used for temporal
! storage of data. Parameter JJ defines the number of ellements of
! static array. Both N and M in this case must be less or equal to JJ.
!
! In: complex AM(N,M) - complex field or its spectrum
!     integer N - number of rows  - integer power of two <= JJ
!     integer M - number of columns  - power of two <= JJ
!     integer IC - direct (IC=+1) or inverse (IC=-1) transform
!
! Out:AM(N,M) - complex field or its spectrum
!     N -  unchanged
!     IC - unchanged
!
! Calls to external FUNCTIONS and SUBROUTINES:
!     CALL fort1(XOK,N,IC)
*/
template <class ddouble>
int FFT2D(dcomplex <ddouble> **AM,int N,int M,int IC);
//---------------------------------------------------------------------------
#endif
//---------------------------------------------------------------------------

fft.h

//---------------------------------------------------------------------------
#include "fft.h"
//---------------------------------------------------------------------------
/*
! 1D Fast Fourier Transform.
! Replaces XOK by its discrete Fourier transform,
! if IC is input as 1 or replaces XOK by its inverse discrete
! Fourier transform, if IC is input as -1. XOK is a complex array
! of length N. N  must be an integer power of 2. This is checked for.
!
! In: complex xok(n) - complex field or its spectrum
!     integer n - array dimension  - integer power of two
!     integer ic - direct (ic=+1) or inverse (ic=-1) transform
!
! Out: xok(n) - complex field or its spectrum
!      n -  unchanged
!      ic - unchanged
!
! Calls to external FUNCTIONS and SUBROUTINES:
*/
template <class ddouble>
int FFT1D(dcomplex <ddouble> *xok, int n, int ic)
{
	ddouble pi = 3.14159265358979;
	dcomplex <ddouble> u, w, t;
	ddouble si, co, pi1;
	int mm, m, i, nv2, nm1, le, le1, j, k, ip, l;
	
	//! calculate degree
	mm = n;
	m = 1;
	for( i = 1; i <= n; i++ )
	{
		mm = mm / 2;
		if( mm < 2 ) break;
		m = m + 1;
	}
	
	//! must be integer of 2
	if( (int)pow(2, m) != n ) return(-1);

	nv2 = n/2;
	nm1 = n - 1;
	j = 1;

	for( i = 1; i <= nm1; i++ )
	{
		if( i < j )
		{
			t = xok[j-1];
			xok[j-1] = xok[i-1];
			xok[i-1] = t;
		}
		k = nv2;

		gt6: if( k < j )
		{
			j = j - k;
			k = k/2;
			goto gt6;
		}
		j = j + k;
	}

	for( l = 1; l <= m; l++ )
	{
		u.re = 1.0; u.im = 0.0;
		le = pow(2.0, l);
		le1 = le/2.0;
		pi1 = pi/(ddouble)(le1);
		si = -ic*sin(pi1);
		co = cos(pi1);
		w.re = co; w.im = si;
		
		for( j = 1;j<=le1;j++)
		{
			for( i = j; i <= n; i += le )
			{
				ip = i + le1;
				t = xok[ip-1]*u;
				xok[ip-1] = xok[i-1]-t;
				xok[i-1] = xok[i-1]+t;
			}
			
			u = u*w;
		}
	}
	
	// Delenie na n
	if( ic > 0)
	{
		for( i = 1; i <= n; i++ )
		{
			xok[i-1].re = xok[i-1].re/(ddouble)n;
			xok[i-1].im = xok[i-1].im/(ddouble)n;
		}
	}
	
	return 0;
}
//---------------------------------------------------------------------------
/*
! 2D FFT. Replaces AM by its 2D discrete Fourier transform,
! if IC is input as 1 or replaces AM by its inverse discrete
! Fourier transform, if IC is input as -1. AM is a complex array
! (N*N2). Both N and N2 must be an integer power of 2.
!
! NOTE  Internal static complex array XOK is used for temporal
! storage of data. Parameter JJ defines the number of ellements of
! static array. Both N and M in this case must be less or equal to JJ.
!
! In: complex AM(N,M) - complex field or its spectrum
!     integer N - number of rows  - integer power of two <= JJ
!     integer M - number of columns  - power of two <= JJ
!     integer IC - direct (IC=+1) or inverse (IC=-1) transform
!
! Out:AM(N,M) - complex field or its spectrum
!     N -  unchanged
!     IC - unchanged
!
! Calls to external FUNCTIONS and SUBROUTINES:
!     CALL fort1(XOK,N,IC)
*/
template <class ddouble>
int FFT2D(dcomplex <ddouble> **AM,int N,int M,int IC)
{
	int L, K;
	dcomplex <ddouble> *XOK;
	if( ( XOK = (dcomplex <ddouble> *)malloc(M*sizeof(dcomplex <ddouble>)) ) == NULL ) return 1;

	for( L = 1; L <= M; L++ )
	{
		FFT1D( AM[L-1], N, IC);
	}
	
	for( K = 1; K <= N; K++ )
	{
		for( L = 1; L <= N2; L++ )
		{
			XOK[L-1] = AM[L-1][K-1];
		}
		
		FFT1D(XOK, M, IC);
		
		for( L = 1; L <= M; L++ )
		{
			AM[L-1][K-1] = XOK[L-1];
		}
	}
	
	free(XOK);
	return 0;
}
//---------------------------------------------------------------------------