Re: [eigen] A complex FFT for Eigen

[ Thread Index | Date Index | More Archives ]

On Fri, Nov 28, 2008 at 1:51 AM, Gael Guennebaud <gael.guennebaud@xxxxxxxxx> wrote:


here are my thoughts.

First of all, I'm 100% OK to provide a Fourier module in Eigen.

One issue with FFTW is that it is GPL, not LGPL, so it is not possible for many companies to use FFTW (despite its world class performance).

It would still be useful to have FFT (including non-powers of 2) for Eigen under the Eigen license terms.

My issues are more about how to do it. Looking at FFTW benchmarks, it is pretty clear to me that we *must* provide support for it. So, in the same vein of the other modules, I would provide a default implementation based on Tim's lightweight source code such that we can make the dependency to FFTW optional for most use cases of this module. Of course this requires to extend a bit Tim's code to support at least: real data, dynamic sizes, 1D and 2D data (I guess all are trivial to add). IMO non pot sizes is not mandatory.

About current Tim's code, what about adding a specialization of Radix2_Decimation<> for sizes greater than, let's say 32, which would use normal recursive calls untill the size is 32. This would automatically add support for dynamic sizes, and significantly reduce the size of the generated code.

About the vectorization of complex types, this is indeed perfectly doable. However note that for std::complex<double> the size of the packet would be 1 ! So, somehow, this means that for double it is rather the std::complex's operators that should be updated to use the SSE instruction set....


On Fri, Nov 28, 2008 at 3:12 AM, Benoit Jacob <jacob.benoit.1@xxxxxxxxx> wrote:
2008/11/27 Tim Molteno <tim@xxxxxxxxxxxxxxxxxxx>:
> I'd love to vectorize this better. I'll look into this, but ran into the
> issue of addressing sub-vectors in-place. I could not work out how to do it
> in Eigen. How do I carry out operations like
> data[1:N] += data[N+1:2N]

assuming that data is a vector of length 2N :

data.start<N>() += data.end<N>();


data.start<N>() += data.segment<N>(N);

(I think that segment only got that name after beta1 so you should use
trunk and "make doc" instead of using the online documentation).

> In other words, can I create a vector object that refers into part of
> another vector object?
> I tried using MatrixBlock objects, but could not get it to work?

Yes Eigen::Block is exactly what you need here and the above-mentioned
start/end/segment methods return objects of type Block.

> I think that FFT's are inherently non-local operations.

OK, so be it, then :)

> There might be
> better ways of doing this, but as it stands this code is pretty fast for an
> FFT and extremely small and lightweight, so I'd rather vectorize the apply
> function first as this is where most of the computation is happening.

Good to know.
As I said for now Eigen doesn't vectorize operations involve complex
numbers but I don't see anything making it impossible in theory, I'll
think about that when I can.

So, from Matthias's and your answers, I gather that there is no
generalization to non-complex numbers or to non-power-of-two sizes.
What about dynamic-size: are there, in your experience, use cases
where the size N isn't known at compile-time and it's be useful to
have a FFT able to handle any size N as a runtime parameter? If yes,
that's also something that we want.

More nitpicking:
Why are fft() and ifft() methods of a class CFFT ? I'd suggest either
making them static methods, or global functions. In Eigen, we tend to
do the latter as it makes for shorter calling syntax.

Let's give us a bit more time to determine whether we want FFT inside
Eigen -- don't want to sound un-enthusiastic, it's rather that I'm
currently busy and I don't want to make the wrong decision.

Anybody reading this, opinions?


> Cheers
> Tim
> ---


Mail converted by MHonArc 2.6.19+