Re: [eigen] A complex FFT for Eigen

["Benoit Jacob" <jacob.benoit.1@xxxxxxxxx>]

Hi,

Thank you for your interest and your code, it's at least interesting to read.

Eigen is a linear algebra library, so FFT is beyond its basic scope,
but we do have a couple of modules doing non-linear-algebra stuff
(Eigen/Regression and Eigen/Geometry) so FFT is not ruled out at that
level.

I lack knowledge on FFT actually to make a decision whether or not we
want it inside Eigen:
- is it usual to restrict it to complex numbers and power-of-two fixed
sizes ? Wouldn't it be desirable to allow reals, to allow any
non-pow-of-2 size, to allow dynamic-size? If any is "yes" then we need
it before inclusion in Eigen
- your code doesn't leverage Eigen as much as it could to enable high
performance. For example, in apply(), replacing the for loop by single
operations on vectors would let Eigen handle loop unrolling and
vectorization. Well, except that Eigen doesn't vectorize std::complex
currently, so I wonder if your code should be adapted to that or if
this means that we should consider vectorizing std::complex. (So far
we haven't seen many uses for std::complex anyway so we didn't
optimize it as much as we maybe could).
- the factorize() function here is not vectorizable at all because it
doesn't do any contiguous access. Aren't there other, more
hardware-friendly algorithms?

Side note:
data /= n; is optimized by Eigen as *= 1/n.


> 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 :(

Not sure if that's your question, but, for lack of template typedefs
in C++98, there is no nicer way to refer to general N-vectors than
Matrix<Scalar, N, 1>. Eigen doesn't have a Vector class, vectors are
matrices. Note there are typedefs like Vector3f, VectorXf etc. We're
looking forward to C++0x ;)

Cheers,
Benoit



2008/11/27 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
>
>
>
>
>

---