[eigen] A complex FFT for Eigen

[Tim Molteno <tim@xxxxxxxxxxxxxxxxxxx>]

Hi,

Eigen is great. Thanks to everyone who has done so much work on it.
I've written a complex fft algorithm for eigen. It uses template 
metaprogramming, and is fairly snappy.

It defines a templated class CFFT, and here is how to use it.

   #define LOG_LENGTH 13
   #define N 1<<LOG_LENGTH

   CFFT<LOG_LENGTH, double> gfft;
   Matrix<fft_complex, LENGTH, 1> cdata;
   cdata.setRandom();

   gfft.fft(cdata);    // Forward transform
   gfft.ifft(cdata);    // Reverse transform

I did run into one issue, which is how to create an eigen object that 
refers into a vector nicely (there is a TODO in the code at the point 
where this would be nice to have). I ended up using pointers into the 
arrays :(

Anyway, I have licensed it under the GPL 3 at the moment, please let me 
know if you'd like to include this in Eigen. I can arrange testharness 
code, and would welcome suggestions for improvements.

Cheers

Tim Molteno
Physics Department
University of Otago




#ifndef _fft_complex_h_
#define _fft_complex_h_
/*
	A FFT and Inverse FFT C++ class library for complex arrays

	Author: Tim Molteno, tim@xxxxxxxxxxxxxxxxxxx

	Based on the article "A Simple and Efficient FFT Implementation in C++" by Volodymyr Myrnyy
	with just a simple Inverse FFT modification. This class uses Eigen http://eigen.tuxfamily.org
	as the Container class for doing its complex arrays.

	Copyright Tim Molteno, 2008.
	Licensed under the GPL v3.
*/


#include <cmath>
#include <complex>
#include <algorithm>

#include <Eigen/Core>
USING_PART_OF_NAMESPACE_EIGEN

/*!\brief Radix-2 Decimation
*/
template<unsigned N, typename T=double>
class Radix2_Decimation
{
	Radix2_Decimation<N/2,T> next;
	enum { N2 = N>>1 };
public:
	void apply(std::complex<T>* data, int iSign)
	{
		//TODO Modify the following to take a MatrixBlock object (reference inside a Vector)
		//	Problem is that I'm not quite sure how to do this so at the moment I'm using
		//	a pointer to the data. Oh well.
		next.apply(data, iSign);
		next.apply(data + N2, iSign);

		T wtemp = iSign*std::sin(M_PI/N);
		T wpi = -iSign*std::sin(2*M_PI/N);

		const std::complex<T> wp(-2.0*wtemp*wtemp, wpi);

		std::complex<T> w(1.0,0.0);

		for (unsigned i=0; i<N2; i++)
		{
			std::complex<T> temp(data[i+N2]*w);

			data[i+N2] = data[i]-temp;

			data[i] += temp;
		
			w += w*wp;
		}
	}
}; 


/*!\brief N=1 case for decimation.
*/
template<typename T>
class Radix2_Decimation<1,T>
{
public:
	void apply(std::complex<T>* data, int iSign) { }
};


/*!\brief Create a templated FFT/Inverse FFT object

	// How to use this FFT.....

	#define LOG_LENGTH 13
	#define N 1<<LOG_LENGTH

	CFFT<LOG_LENGTH, double> gfft;
	Matrix<fft_complex, LENGTH, 1> cdata;
	cdata.setRandom();

	gfft.fft(cdata);	// Forward transform
	gfft.ifft(cdata);	// Reverse transform
*/
template<unsigned P,typename T=double>
class CFFT
{
	enum { N = 1<<P };
	Radix2_Decimation<N,T> decimation;

	/*!\brief Cooley-Tukey factorization 
	*/
	void factorize(Matrix< std::complex<T>, N, 1>& data)
	{
		int j=1;
		for (int i=1; i<N; i++)
		{
			if (j>i)
			{
				std::swap(data[j-1], data[i-1]);
			}
			int m = N/2;
			while (m>=2 && j>m)
			{
				j -= m;
				m >>= 1;
			}
			j += m;
		}
	}

public:
	/*!\brief Replaces the complex array data[1..N] by its DFT */
	void fft(Matrix< std::complex<T>, N, 1>& data)
	{
		factorize(data);
		decimation.apply(data.data(),1);
	}

	/*!\brief Replaces complex array data[1..N] by its inverse DFT */
	void ifft(Matrix< std::complex<T>, N, 1>& data)
	{
		factorize(data);
		decimation.apply(data.data(),-1);

		/** Scale all results by 1/n */
		T scale = static_cast<T>(1)/N;
		data *= scale;
	}
};

#endif /* _fft_complex_h_ */