Re: [eigen] FFT for Eigen

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

Hi Mark,

First of all, FFT is an area where I'm very incompetent so I'm hoping
that others will correct me if I'm saying something wrong.

For Eigen we wanted a FFT module that
1) has a generalistic API that covers many use cases
-->KISSFFT seems to have that since it allows dynamic non-power-of-2
size while having small fixed-size specializations, allows
real/complex, etc.
2) provides a reasonably small implementation that still has reasonable speed
--> I checked, KISSFFT has less than 500 lines of C
3) allows to use the most efficient libraries as optional backends.
It's OK if certain backends have stricter requirements than the
general API -- e.g. if a certain backend is power-of-two-only.
--> that's what you seem to be offering at the end of your email:
> Also, my hope is to have #ifdef protected code that would use the kissfft
> library by default, but could also use FFTW or IMKL if it was so compiled
> (since they are much faster).
I'd suggest having a look at how Gael is doing the selectable-backend
thing in the Sparse module, for the solvers.

So I'd say go ahead, this is very interesting. In the special case of
small power-of-2 sizes, it'd be great if you could compare the speed
with Tim's code.

Organization-wise: the next step is that you get yourself an account
at bitbucket.org, fork eigen/eigen2, go in the unsupported/ directory
and create your FFT module there. Once you have something, tell us to
have a look at your fork and we can then merge it.

Code-wise: your code is as hardcore C as it gets. Some of that needs
to be c++-ified and eigen-ified, some of that can stay that way. What
absolutely must be eigen-ified is the macros you use to handle
different numeric types: in Eigen we already have infrastructure for
that, see NumTraits.h and MathFunctions.h in Eigen/src/Core/, so
normally your FFT shouldn't bother at all anymore about fixed-point
etc, but now if you find that our infrastructure is missing something
that you need, then tell us and we can add it.

Cheers,
Benoit

2009/5/17 Mark Borgerding <mark@xxxxxxxxxxxxxx>:
> I've been thinking about porting/adapting my kissfft library for use in C++
> (  http://sourceforge.net/projects/kissfft )
>
> I am aware that there was a previous offering for an Eigen FFT library, but
> it seemed to suffer from a couple of a serious drawbacks. Namely, it could
> only perform power-of-2 transforms and the power of two was required at
> compile time.
>
> A little about KISSFFT, since that is the library this would be based on.
>
> arbitrary FFT sizes with optimized radix butterflies for factors 2,3,4,5
> arbitrary dimensions
> optimization for real-only case
> small code footprint
> reasonably good efficiency (within order of magnitude speed of FFTW and
> IMKL).
> can be compiled for fixed point or a variety of floating point (AFAIK,
> kissfft is the only free  fixed point FFT , period)
> easy going license (currently revised BSD style, for Eigen, I would be
> willing to use LGPL/GPL)
>
> My thought for Eigen was to make a template fft class that would be
> specialized for: float, double, & long double, eventually maybe others
> (fixed point, extended precision classes?)
> Much of the behavior would be relegated to a traits class which would be
> responsible for creating twiddle factors, deciding factorization, scaling,
> etc.  This should allow easy extension by users without requiring them to
> hack the library itself.
>
>
> I've done a little bit of coding and little bit of thinking.  The interface
> is shaping up to something like:
>
> template <typename T_Data,typename T_traits=kissfft_utils::traits<T_Data> >
> class kissfft
> {
>     public:
>         typedef T_traits traits_type;
>         typedef T_Data scalar_type;
>         typedef std::complex<scalar_type> cpx_type;
>
>         kissfft(int nfft, const traits_type & traits = traits_type() );
>
>        void fwd( cpx_type * dst, const cpx_type * src); // forward
> complex-to-complex transform
>        void fwd( cpx_type * dst, const scalar_type * src);  // forward
> real-to-complex transform
>        void inv( cpx_type *dst, const cpx_type * src);// inverse
> complex-to-complex transform
>        void inv(scalar_type * dst, const cpx_type *src; // inverse
> complex-to-real transform
> }
>
> Also, my hope is to have #ifdef protected code that would use the kissfft
> library by default, but could also use FFTW or IMKL if it was so compiled
> (since they are much faster).
>
> -- Mark
>
>