Re: [eigen] FFT for Eigen

[ Thread Index | Date Index | More lists.tuxfamily.org/eigen Archives ] 

Benoit,
Thanks for the feedback. I'm excited about contributing to Eigen. I feel it is a worthwhile project that addresses a previously unfilled niche within the scientific computing community.
Thanks for the tip about the Sparse Module. I will take a look. It looks like there is some useful stuff in Numtraits.h and MathFunctions.h. I wish the STL had some of that. For now I am just assuming the use of the standard C++ library <complex> and <vector> headers, but its good to know. 

Here's a rough stab at what I see as a chronological order of FFT "releases"

  1. one-dimensional FFTs, forward and inverse for float, double,long
     double  (i.e. minimal functionality) this is about 90% complete
  2. real-only optimization (this involves mostly cannibalizing the
     existing kissfft C code)
  3. multi-dimensional FFTs ( more cannibalizing from C version)
  4. backend for FFTW ( using this would require application to be GPL,
     or the developer to have a commercial MIT license)
  5. backend for intel Math Kernel Library (this allows use in closed
source applications that have paid their licensing fee to intel. Explicitly coding to this interface may not be strictly necessary, 
     since IMKL has a FFTW adapter interface already)
Even though it shows up last in the development order, I *really* want the ability to choose at compile time (maybe even at runtime) between kissfft, FFTW, and IMKL. Allowing the different backends makes sure everybody can get what they need based on their individual requirements regarding licensing, finances, code size, and execution speed. My hope is that this will become the only FFT interface a C++ developer will ever need to learn ;)
FFTW and IMKL both support arbitrary sizes, real FFTs, and multiple dimensions. So I don't think we would need to restrict functionality yet. I could see a need to someday include a CUDA backend (CUFFT). Their FFTs are currently limited to complex power-of two.
You are correct that the existing kissfft code is "as hardcore C as it gets". That was the aim for the original kissfft. There are a lot of embedded systems that lack FPUs, megabytes of RAM and C++ compilers. However, I don't think Eigen will get used much on those systems. The need for macros to handle different types goes away with C++ templates. My point is : Don't worry, this module will not simply be C compiled by a C++ compiler. 

-- Mark


Benoit Jacob wrote:

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

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