Re: [eigen] Eigen/FFT

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

I forgot a .eval() here:

template<typename SrcDerived, typename DstDerived>
inline void foo(const MatrixBase<SrcDerived>& src, MatrixBase<DstDerived> *dst)
{
   foo_impl(src.eval(), dst);
}

template<typename SrcMatrixType, typename DstDerived>
void foo(const SrcMatrixType& src, MatrixBase<DstDerived> *dst)
{
   foo_impl<typename SrcDerived::PlainMatrixType, DstDerived>(src, dst);
}

2009/6/23 Benoit Jacob <jacob.benoit.1@xxxxxxxxx>:
> (it seems that GMail now defaults to HTML mode. Is it worth that i
> click everytime on "text-only" or will it still send a text-only
> version alongside the HTML version if I leave HTML mode enabled ?)
>
> OK, finally here's a review of FFT!
>
> First of all thanks a lot for the good work and the wiki page, it's very useful.
>
> 2009/6/16 Mark Borgerding <mark@xxxxxxxxxxxxxx>
>> I would like to add fwd & inv functions for the Eigen Matrix types into the FFT interface, but I could use a little guidance on how to do it since I am still an Eigen novice.
>
> your function could have a prototype like this:
>
> template<typename SrcDerived, typename DstDerived>
> void foo(const MatrixBase<SrcDerived>& src, MatrixBase<DstDerived> *dst);
>
> now if each coefficient in src is going to be used many times, it is
> probably a good idea to evaluate the expression src into a matrix, if
> it wasn't already a plain matrix. To do so, use
> typename SrcDerived::PlainMatrixType
>
> The ultimate refinement is to make this foo() a thin wrapper around a
> function foo_impl() that is template only in the "plain matrix type"
> so that instantiating foo on many different Src types doesn't result
> in the whole code being instantiated everytime...
>
> template<typename SrcDerived, typename DstDerived>
> inline void foo(const MatrixBase<SrcDerived>& src, MatrixBase<DstDerived> *dst)
> {
>    foo_impl(src, dst);
> }
>
> template<typename SrcMatrixType, typename DstDerived>
> void foo(const SrcMatrixType& src, MatrixBase<DstDerived> *dst)
> {
>    foo_impl<typename SrcDerived::PlainMatrixType, DstDerived>(src, dst);
> }
>
> Of course that's not always needed, depends how your function is going
> to be used in practice, if it's big, etc.
>
>>
>> ### Design choices:
>> I should mention that in the FFT module, I've made a couple of design decisions that stray away from the behavior of most FFT libraries.  The intent is to facilitate generic programming and ease migrating code from Matlab/octave.   I've based these decisions on my experience both using FFTs and answering questions about KISSFFT.  I think the default behavior of Eigen/FFT should favor correctness and generality over speed.
>
> Agree!
>
>> 1) Scaling:
>> Other libraries (FFTW,IMKL,KISSFFT)  do not perform scaling, so there is a constant gain incurred after the forward&inverse transforms , so IFFT(FFT(x)) = Kx;  this is done to avoid a vector-by-value multiply.  The downside is that algorithms that worked correctly in Matlab/octave don't behave the same way once implemented in C++.
>
> I don't have the experience, but since you have, I believe that your
> solution can be more convenient in practice!
>
> This forces a retraversal of the array though so i wonder about the
> performance cost? First, O(N) is not that small compared to O(N log N)
> ; second, i'm concerned about the re-traversal of the array, which can
> be bad for caches on large data?
>
> By the way, you have room for improvement in the implementation of the
> scaling: by letting Eigen do it, you would have the vectorization. So
> this:
>
>        Scalar s = 1./nfft;
>        for (int k=0;k<nfft;++k)
>          dst[k] *= s;
>
> becomes:
>
>        Matrix<Scalar,Dynamic,1>::Map(dst, nfft) /= nfft;
>
> Likewise, in KissFFT, you can optimize scale().
>
>> 2) Real FFT half-spectrum
>> Other libraries use only half the frequency spectrum (plus one extra sample for the Nyquist bin) for a real FFT, the other half is the conjugate-symmetric of the first half.  This saves them a copy and some memory.  The downside is the caller needs to have special logic for the number of bins in complex vs real.
>>
>> How Eigen/FFT differs: The full spectrum is returned from the forward transform.  This facilitates generic template programming by obviating separate specializations for real vs complex.  On the inverse transform, only half the spectrum is actually used if the output type is real.
>
> Again, i don't have the experience to express an opinion here.
>
> Now for more comments:
>
> * The size "nfft" is a runtime variable, which is great, but (honest
> naive question) isn't there also a use case for tiny sizes that are
> known at compile time? If yes, we can help making it work the Eigen
> way, the Eigen-ish approach is to add a NfftAtCompileTime parameter
> that can take the value Dynamic in which case the runtime variable
> nfft is used, otherwise it plays the role of the nfft variable. Of
> course that is a bit of work as there are template specializations to
> do to take advantage of the fixed size when possible. Also I
> understand that this probably would only affect the KissFFT backend.
> But it would then beat hands down all others, fixed size
> specializations mean 10x speed improvements.
>
> * The backend selection should IMO work like in the Sparse module:
> require the user to define a symbol to tell which backend to use,
> don't use whatever was included before Eigen. This way, there are no
> surprises, and it becomes possible to choose your backend for
> Eigen/FFT regardless of what you use in other parts of your program.
> (One thing though, Sparse uses #include"..." but that should be
> #include<....>, will fix this)
>
> * I'm uncomfortable about:
> typename _Traits=DEFAULT_FFT_IMPL<_Scalar>
> because 1) it's strange to have 'traits' as shorthand for an "impl",
> why not call it impl instead of traits, and 2) I don't think that
> there's a justification for using a macro here (a typedef would do)
> and 3) i'd rather have DEFAULT_FFT_IMPL be an integer (an enum) rather
> than a class name, and then the actual impl class could take that
> integer as template parameter and you would have partial template
> specializations for each possible value. It amounts to the same, but
> somehow feels cleaner to me.
>
> * We have a loose rule in Eigen that a "src" argument is a const
> reference, and a "dst" argument is a pointer. This helps the user
> catch a programming error where he mixes src and dst. Though at some
> places we have dst as a reference. What we don't have AFAIK is src as
> a pointer. Anyway you choose if it's a good rule in the context of
> your FFT module.
>
> That's all I can see for now ;)
>
> Feel free to push FFT to eigen/eigen2 when you want, since you put it
> in unsupported there's no reason not to import it. Do you want to get
> write access and push it, or should we pull it?
>
> As for when to move it to Eigen, I think that besides finishing the
> TODOs the main factors are if you think that you can commit to be the
> maintainer for this module, and if there are people starting to use
> it.
>
> Cheers,
> Benoit
>