...

My code, when run as a 2048-point FFT, agrees with the MIT FFTW one to a max error margin of 6 parts per million.

This might be of interest: I've ported FFTW's arbitrary-precision FFT to boost::multiprecision

Awesome. I am a co-author of multiprecision and interested in your work. May I ask... Which FFTs did you concentrate on, dimensions, complex, etc.? What precisions did you use? Which multiprecision FP backend did you use? How was the performance? Where do you see the main application of your work in this area? How was your experience with multiprecision (critical suggestions are OK)?

Generally, I think it would be better if a boost::fft library would primarily be a wrapper around existing FFT libraries, with the C++ implementation only used as a fallback for multiprecision or licensing issues since it's unlikely a template implementation would catch up with the years of optimization work that went into single- and double-precision libraries like MKL and FFTW, especially FFTW's kernel generator and planner.

I completely agree. This identically mirrors the philosophy used in multiprecision.

But a reasonably fast (and open) multidimensional multiprecision FFT implementation doesn't seem to exist yet.

I know. You're preaching to the choir. Got to put FFT on the back burner. GSoC 2014? Sincerely, Chris.

...

...

This might be of interest: I've ported FFTW's arbitrary-precision FFT to boost::multiprecision

Which FFTs did you concentrate on, dimensions, complex, etc.?

only 1D complex (radix-2 Cooley-Tukey, Bluestein for non-power-of-2 sizes), using std::complex with the multiprecision type.

Just so you know - std::complex isn't guarenteed to work with user defined types (and often doesn't sadly) - we really must get around to adding a complex number adapter to the multiprecision lib.

...

How was your experience with multiprecision (critical suggestions are OK)?

Straightforward, this is the first time I've used it; there was an inconsistency between float50 and float instances because pow(value,2) occasionally seemed to return wrong results which caused some confusion thanks to template magic; I ended up using boost::math::pow<2>(value) which didn't show the same problem. I think it was a problem in my code, not a bug in boost.

If you used 1.53.0 there was a bug in the expression template code that caused pow to misbehave in some situations, basically if you wrote: a = pow(a, n); rather than: b = pow(a, n); It's fixed for the next release. Regards, John.

On 3/06/2013 7:39 AM, Mathias Gaunard wrote:

It's a pure computation function; the interface should just be a couple of functions taking pointers.

That would be a very naive and pedestrian approach. Factors ranging from the input type (real, complex, integer), dimensionality of the data, the way the data is stored, the precision required, application of the transform (lets think about multiplications for Strassen), phase rotation et al. And that's all before we realize that hey FFTs are really just a subset of transforms in general such as Wavelets and Co. So what we should be heading for is not another CRC library (which should really be part of a generalized message digest library), but rather a well thought-out and designed interface that is extensible with provisions for future support of multiple backends (fftw, ipp, cuda, opencl et al) and types (float, double, mp, et al) - and not just the first meta idea that easily scales in our minds. ;) In short - A couple functions with a few pointers would fall dismally short of the 'mark'.

On 03/06/13 14:15, Arash Partow wrote:

On 3/06/2013 7:39 AM, Mathias Gaunard wrote:

...

It's a pure computation function; the interface should just be a couple of functions taking pointers.

That would be a very naive and pedestrian approach.

Factors ranging from the input type (real, complex, integer), dimensionality of the data, the way the data is stored, the precision required, application of the transform (lets think about multiplications for Strassen), phase rotation et al.

Different function names for different variants. You may also use overloading.

So what we should be heading for is not another CRC library (which should really be part of a generalized message digest library), but rather a well thought-out and designed interface that is extensible with provisions for future support of multiple backends (fftw, ipp, cuda, opencl et al) and types (float, double, mp, et al) - and not just the first meta idea that easily scales in our minds. ;)

In short - A couple functions with a few pointers would fall dismally short of the 'mark'.

Yet it is what is needed to fulfill your requirements.

On 3.6.2013 16:48, Mario Mulansky wrote:

afaik fft methods need temporary working space. this might already be a strong enough argument to encapsulate methods into classes. I dont think "just a few functions" is enough here.

Actually, no, you don't need any temporary space for FFT/IFFT. There might be some implementations that require this, but by all means that doesn't cover all FFT implementations. At least with the basic radix-2 implementations there's even no need for separate input and output buffers; All you need is one (complex) buffer of size N, where N is your FFT size and everything happens in place. (If you're doing some contiguous/real life signal processing, obviously you can't have a buffer of an unlimited size. Then you need to split your input buffer into blocks and this is where overlap add/save etc. techniques step in; And these need temporary buffers between blocks. However, the temporary buffers needed for these are "outside" the FFT/IFFT calculation.) It's not that I would need a Boost.FFT, but I completely agree with the "few just won't cut it" ideology. Mostly because there are so many free, tested and ready implementations around that work on the POD/trivial complex types. 100-200 lines worth of good ole copy-paste and you're ready to go. However, what the most - if all - lack is the ability to use actually complex (not complex as in complex number) types that would reveal and enable one to test how many bits are used or required to implement a FFT of a certain precision. That might not be a good use case for a rainy day project, but if you're implementing FFT related things (say, a digital filter; FFT, filter, IFFT) on a embedded system you don't just start directly writing code for that particular device. What you do is first verify (protype) your plans with eg. Matlab, then write a C/C++ version that runs on your desktop and only after that port the code for the device. As it could be, that you don't know until you have your PoC running that what kind of embedded platform you do need for your implementation to be feasible. Do you need a DSP chip, 16bit or 32bit, floating or fixed point or does a multi purpose micro controller cut it. And that's the point where you'd like stick your custom complex number class in; That tells you at what kind of level (as in how many bits) precision it operates and what kind of operations it makes. This is a bit what Matlab's Fixed Point toolbox does (I have used it a bit) but still I find a bit lacky for real life problems. I'd say that in the end you'll write your FFT completely by hand as that's where all the processing time counts. But! You don't want to start by doing yet another FFT implementation that is off-by-one bug etc. ridden for the first few days. As I'm not the one doing the implementation, please do regard this as a comment coming from the back row - easier said than done, that is. Just liked to share some thoughts on FFT matters that I've encountered during the years. -- Pekka

In short - A couple functions with a few pointers would fall dismally short of the 'mark'.

I agree with this. We need to define an interface that can accommodate the various forms of usage and implementation of the DFT. Even if the first implementation of that interface is limited in terms of performance. So far everyone has identified a number of things that may belong in a DFT interface. Maybe we need to work on a list. However I don't know if the original poster is up for this kind of work. He proposed a simple implementation of FFT that covered a single usage. Some of the things to support so far in the interface include: * Abstraction of data it operates on (std::complex, fixed point, ...) * Abstraction of dimensionality * Abstraction of memory layout of data (Is it strided, are the real/imag interleaved, if doing multiple dimensions how are they arranged in memory). Using iterators is probably sufficient here * Size of the transform (supporting non power of 2 is important) * Does it use temporary storage or not On this point I think this is a must. If for example you want to support using FFTW as a backend implementation, you want to be able to store data across transform calls. * Is the transform to happen in-place or not * Do we allow definition of phase rotations as part of the interface or require them as pre/post processing steps I am not sold on this as part of the transform itself. * Do we want to support (probably separate interface) real (vs complex) DFT transforms? * How will we support backend customizations? Allowing a single interface to be backed by say FFTW is probably a very good thing in this situation. Not just FFTW but many embedded platforms provide their own custom implementation of the FFT that a user may wish to utilize, but wrap it in a compatible API. I dont agree that we should have one "transform" that is specialized for DFT, DCT, ... But whatever interface we define is probably worth thinking how it applies to other similar transforms. I.e. The usage could be almost identical. I have been looking at this as I am working on a concept for a library called Boost.GAL (Same as GIL but for audio). In doing so I feel that GIL has been made too specific working on 2D data for images only, as opposed to a DSP library that can work in N dimensions. I.e. If we define a DFT or DCT, having that same algo work both for a Graphics or Audio library and any other kind of DSP application is probably a good thing.