[eigen] Re: Eigen::FFT using Matrix

[Mark Borgerding <mark@xxxxxxxxxxxxxx>]

Matt,
Thanks for your interest and questions.  I hope you do not mind; I'm replying to your questions on the mailing list.  There may be other interested parties. 

== Stability and Accuracy
As far as the stability and accuracy of FFT goes.  It's very solid.  There are two backends: one based on kissfft and one that uses fftw. Both have been stable for years with thousands of users.
Last year, Pavel Holoborodko reported that `FFT<mpreal>` (only possible with the kissfft backend) was giving him 300 digits of accuracy.
http://listengine.tuxfamily.org/lists.tuxfamily.org/eigen/2012/04/msg00021.html
(The issue he noted with non-powers of two was resolved shortly thereafter)

FFTW, may not be specializable to arbitrary types, but it does have speed competitive with the best libs out there, at least on a general CPU. (hint: `configure` it to use simd)

The IMKL backend *will* be even faster than FFTW on intel machines.  (TODO)

== Matrix Args

Regarding the Matrix arguments to fwd.  You are correct that it does not work with non-vector types.  I hope to address that at some point in the future.
Nit: In your example, I think maybe you meant to put something other than "1". Matrix<float,Dynamic,1> is VectorXf, which works fine.

== Roadmap

I hope to refactor the FFT interface -- to make it more like the DFT unitary linear transform:
    e.g. in Matlab,octave  we can model the DFT as a linear operator

octave:1> DFT=fft(eye(16))/sqrt(16);   <---- This is what I'd like the Eigen:FFT module to be more like
octave:2> x=randn(16,1)+j*randn(16,1);
octave:3> norm(DFT'*DFT - eye(16))     <--- orthonormal
ans =  1.2585e-16
octave:4> norm(x) - norm(DFT*x)            <--- preserves lengths
ans = 0

== TODO:
- Scale the forward,inverse FFTs to be orthonormal.  The forward,inverse FFTs will be scaled by 1/sqrt(nfft) , sqrt(nfft) with respect to matlab
- use return value proxies, this necessitates making different function names for the real vs complex fft to disambiguate the inverse case.
  - e.g. crossSpec += _dft.fft(x).cwiseProduct(  _dft.rfft(y) );
- allow FFTs of matrices to be accomplished 1 column vector at a time, like Matlab does.
- add 2d FFTs
- IMKL backend

== Interface Changes

The TODO list will obviously change the interface (compile-time errors) and behavior(latent errors). The class may be in `unsupported/Eigen`, but the polite thing to do is have *frozen* version of the current interface for existing users. This should allow existing users to get their existing code working with little or no effort.  Later they can cycle back around and switch over to the newer class.
- One option is to rename the current FFT -> FFTLegacy.  When the interface changes cause code to not compile, users can make trivial changes.
- Another option is use a new name for the new code e.g. `FastFourierTransform`.  The old `FFT` would be deprecated and eventually removed.    <-- I'm leaning this way



I'd welcome your opinions (or coding help, for that matter)

-- Mark


On 03/07/2013 06:40 PM, Matt Flax wrote:

Hi Mark,

Hows it going ?

I noticed that you implemented the FFT code for the Eigen system. It is fantastic that you have done that... thanks.

I have a few questions I was hoping that you could answer for me.
First of all, I would like to use Eigen::Matrix  data types with the fft. Is that possible ?

Something like :
Matrix<float, Dynamic, 1> b;
Matrix<complex<float>, Dynamic, 1> B;

FFT.fwd(B,b);

However that doesn't work.

From the code, I think I can use the following :
FFT.fwd(B.data(), b.data(), b.size());

However I also notice that there is a definition for the 'fwd' method using Matrix types ... can you please give me an example on how to use it ?

Also, is FFTW and Eigen optimised for various platforms ?
Can you give me a feel for how stable/accurate it is ?

thanks
Matt