[eigen] A complex FFT for Eigen

```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.

*/

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

#include <Eigen/Core>
USING_PART_OF_NAMESPACE_EIGEN

*/
template<unsigned N, typename T=double>
{
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>
{
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 };

/*!\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_ */

```

 Mail converted by MHonArc 2.6.19+ http://listengine.tuxfamily.org/