[boost] Request to contribute boost::FFT

[Nathan Bliss]

Dear Boost Community Members,
I am writing to ask if I could contribute a C++ FFT implementation to the boost library. I noticed that there is a boost::CRC so I thought it would be a good addition to have boost::FFT as well. MIT have an existing open-source FFTW library, but it is GPL-licensed which is much more restrictive than the Boost license. I should imagine a commercial FFTW license would be very expensive.
I have working FFT code which I have tested on Ubuntu Linux GNU gcc and also Visual Studio Express 2012 (MS VC++) for Windows. 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.
It is implemented using templates in a .hpp file and is very easy to use:
------------------------------------------------------------------------------------------------------------[fft.hpp]
template<class FLOAT_TYPE, int FFT_SIZE, int NUM_THREADS>class FFT;///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////[Invocation example]
typedef double My_type;FFT<My_type, 2048, 4> user_fft;user_fft.initialise_FFT();get_input_samples(*user_fft.input_value_array, user_fft.m_fft_size);user_fft.execute_FFT();print_output_csv(user_fft.output_value_array, 2048);------------------------------------------------------------------------------------------------------------
It's uniqueness compared to other FFT implementations is that it fully utilises the boost::thread library, to the extent that the user can give the number of threads they want to use as a parameter to the class template (above case is for 4 parallel threads).
It is structured and organised so that users could customise/optimise it to specific processor architectures.
I've also tried to develop it in the spirit of Boost in that all class members which users should not access are private, only making public what is necessary to the API. My code is a decimation-in-time radix-2 FFT, and users could in theory use the existing API as a basis to extend it to more complex implementations such as the Reverse/Inverse FFT, bit-reversal of inputs/outputs and multi-dimensional FFTs.
I look forward to your reply.
Kind regards,Nathan Bliss


_______________________________________________
Unsubscribe & other changes: http://lists.boost.org/mailman/listinfo.cgi/boost

[Andy Jost]

Do you know how the performance compares to FFTW or, say, Intel MKL?

[Jason Roehm]

On 05/28/2013 06:12 PM, Nathan Bliss wrote:
> Dear Boost Community Members,
> I am writing to ask if I could contribute a C++ FFT implementation to the boost library. I noticed that there is a boost::CRC so I thought it would be a good addition to have boost::FFT as well. MIT have an existing open-source FFTW library, but it is GPL-licensed which is much more restrictive than the Boost license. I should imagine a commercial FFTW license would be very expensive.

I certainly don't speak for the Boost community at large, but in my
opinion, an FFT library would be a great addition to Boost. While for
performance-critical applications, I doubt you would be able to approach
the speed of a specially-tuned implementation like Intel's Math Kernel
Library, a portable, parallel-capable, easy-to-use, freely licensed
library that provides decent throughput would be useful. Some thoughts
that come to mind:

- Do you have any benchmark information available? A comparison to
contemporary competitors like FFTW, Intel's MKL on x86 platforms, or
kissfft (http://kissfft.sourceforge.net/) would be useful.

- The 6e-6 relative error specification that you quoted (when compared
to FFTW) seems a bit large to me for double precision. Have you
characterized what causes the difference?

- You might already support some of these, and I've not been involved in
Boost development myself, but from observing the process, I would guess
that an FFT library would need to provide a lot of flexibility in order
to be accepted. Some of the features that I would consider important
would be:

- The ability to compute both forward and inverse FFTs
- Optimized forms for real inputs (i.e. only computing the half
spectrum due to the conjugate symmetry in the result) or real outputs
(i.e. calculating the inverse FFT on a half spectrum that is assumed to
be conjugate symmetric)
- An interface to compute multiple transforms in a batch (which can
improve throughput, especially for smaller transforms), maybe even
customizable input and output strides and distances (like FFTW)
- Multiple backend implementations to allow for non-power-of-2
transform sizes

- One other idea that is somewhat pie-in-the-sky right now: FFT
algorithms are a great fit for SIMD implementations, so perhaps support
for the still-in-development Boost.SIMD library
(http://nt2.metascale.fr/doc/html/the_boost_simd_library.html) would be
a future possibility.

Jason

[Karsten Ahnert]

On 05/29/2013 12:12 AM, Nathan Bliss wrote:
> Dear Boost Community Members,
> I am writing to ask if I could contribute a C++ FFT implementation to the boost library. I noticed that there is a boost::CRC so I thought it would be a good addition to have boost::FFT as well. MIT have an existing open-source FFTW library, but it is GPL-licensed which is much more restrictive than the Boost license. I should imagine a commercial FFTW license would be very expensive.
> I have working FFT code which I have tested on Ubuntu Linux GNU gcc and also Visual Studio Express 2012 (MS VC++) for Windows. 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.

I would really like to see an FFT implementation in boost. It should of
course focus on performance. I also think that interoperability with
different vector and storage types would be really great. It prevents
that the user has to convert its data types in different formats (which
can be really painful if you use different libraries, for example
odeint, ublas, fftw, ...).

[Paul A. Bristow]

> -----Original Message-----
> From: Boost [mailto:boost-...@lists.boost.org] On Behalf Of Nathan Bliss
> Sent: Tuesday, May 28, 2013 11:12 PM
> To: bo...@lists.boost.org
> Subject: [boost] Request to contribute boost::FFT
>

> I am writing to ask if I could contribute a C++ FFT implementation to the boost library.

Definitely, a good templated C++ FFT would be very welcome.

Would/does it work with Boost.Multiprecision to give much higher precision ? (at a snail's pace of
course).

Paul

---
Paul A. Bristow,
Prizet Farmhouse, Kendal LA8 8AB UK
+44 1539 561830 07714330204
pbri...@hetp.u-net.com

[Mathias Gaunard]

You may want to give a look to the FFT functions bundled with NT2
courtesy of Domagoj Saric.

They also generate a FFT for a given compile-time size and use
Boost.SIMD for vectorization (but no threads). Unfortunately the code is
not so generic that it can work with arbitrary vector sizes, so it
limits portability somewhat.

<https://github.com/MetaScale/nt2/blob/master/modules/core/signal/include/nt2/signal/static_fft.hpp>

[Christopher Kormanyos]

> Dear Boost Community Members,
> I am writing to ask if I could contribute a C++ FFT implementation
> to the boost library. I noticed that there is a boost::CRC so I thought
> it would be a good addition to have boost::FFT as well. MIT have
> an existing open-source FFTW library, but it is GPL-licensed
> which is much more restrictive than the Boost license.

Yes, in my opinion, Boost definitely needs FFT support --- even if
only having a radix-2 decimation in frequency to start out.

So it sounds pretty exciting.

<snip>

> I have working FFT code which I have tested on Ubuntu Linux GNU
> gcc and also Visual Studio Express 2012 (MS VC++) for Windows.

Is the code available for public viewing and benchmark?

> 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.
> It is implemented using templates in a .hpp file and is very easy to use:

<snip>

> It's uniqueness compared to other FFT implementations is that it fully
> utilises the boost::thread library, to the extent that the user can give the
> number of threads they want to use as a parameter to the class template
> (above case is for 4 parallel threads).

Just out of interest, can it also use std::thread? Your VS Express 2012
and GCC > 4.7 compilers should support the C++ thread support library.

> It is structured and organised so that users could customise/optimise
> it to specific processor architectures.

Everyone always wants to do this. But I strongly prefer parallel
distribution over these kinds of specializations.

<snip>

> My code is a decimation-in-time radix-2 FFT, and users could
> in theory use the existing API as a basis to extend it to more
> complex implementations such as the Reverse/Inverse FFT,
> bit-reversal of inputs/outputs and multi-dimensional FFTs.

I did some work on radix-2 decimation using template
metaprogramming. Perhaps our research overlaps here.

As Paul mentioned, this is work in progress for Boost.Multiprecision
for ultra high precision multiplication routines. It would really be great
to pull a Boost FFT off the rack instead of rolling one solely for
Boost.Multiprecision. For this application, though, we would also
need the inverse transformation. And the convolution of the
half-complex arrays is nice to have.

> I look forward to your reply.
> Kind regards, Nathan Bliss
                       

Thanks for your heads-up.

Sincerely, Chris.

[Pascal Germroth]

Hi,

> 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
<https://github.com/neapel/vexcl-fft-accuracy/blob/master/multi_precision_fft.hpp>,
if you want to compare the numerical performance of your implementation
with FFTW as in <http://www.fftw.org/accuracy/Pentium4-3.60GHz-icc/>.
It's more meaningful than comparing two low-precision results.

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.
But a reasonably fast (and open) multidimensional multiprecision FFT
implementation doesn't seem to exist yet.


Cheers,

--
Pascal Germroth

[Christopher Kormanyos]

>> 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.

[Pascal Germroth]

Hi,

>>> 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
>

> 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.

> What precisions did you use?
> Which multiprecision FP backend did you use?

static_mpfr_float_50 to check float/double FFTs (and float1000 to check
float50 has the same error behaviour, around ~1e-50).

> How was the performance?

Very good, especially compared to FFTW's original implementation.

> Where do you see the main application of your work in this area?

I use it to verify an OpenCL-FFT implementation when debugging, so
performance isn't that important.

> 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.


Cheers,

--
Pascal Germroth

[Arash Partow]

On 29/05/2013 8:12 AM, Nathan Bliss wrote:
> Dear Boost Community Members,
> I am writing to ask if I could contribute a C++ FFT implementation to
> the boost library.

An FFT library would be an awesome addition to Boost.

Some commenters have mentioned performance of your library when
compared to others. I believe that initially performance should be a
secondary concern over the primary objective which would be 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).

[Mathias Gaunard]

On 02/06/13 23:19, Arash Partow wrote:
> On 29/05/2013 8:12 AM, Nathan Bliss wrote:
> > Dear Boost Community Members,
> > I am writing to ask if I could contribute a C++ FFT implementation to
> > the boost library.
>
> An FFT library would be an awesome addition to Boost.
>
> Some commenters have mentioned performance of your library when
> compared to others. I believe that initially performance should be a
> secondary concern over the primary objective which would be 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).

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

[John Maddock]

>>> 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.

[Arash Partow]

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'.

[Mathias Gaunard]

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.

[Karlsen, Naja (E W EMEA ECC GOV)]

Dear Mathias Gaunard,

I really don't know what you have answered here ;o)

Is this really an answer to MY question??

Can you please help me find out the Export control classification number, ECCN, for Boost?

If you don't have the ECCN, then you can maybe inform me of the content of encryption algorithms?

Thank you in advance.

Med venlig hilsen / Best regards,
Naja Karlsen

Siemens Wind Power A/S
E W EMEA ECC GOV
Borupvej 16
7330 Brande, Denmark
Tel.: +45 99 42-8070
Mobil: +45 30 37-3663
mailto:naja.k...@siemens.com
http://www.siemens.com/wind

[Mario Mulansky]

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.

I would really like to see a generic fft libraries, that can also run on
multi-precision types.

Mario

[Pekka Seppänen]

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

[Mathias Gaunard]

On 03/06/13 16:55, Pekka Seppänen wrote:
>
> 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.

You can only test if precision is sufficient if you use exactly the same
code in your prototype and in your optimized version though.

[Christopher Kormanyos]

>>> 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.

We got close a while ago with the extended complex code
from Matthieu Schaller. I don't know if this is the right ticket.
But it seemed like the right direction last year. As I recall,
Matthieu's code needs to be reviewed and probably needs
to be better documented.

That life-long TODO list keeps getting longer and longer! :-)

Sincerely, Chris.

[Brendon Costa]

> 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.

[Christopher Kormanyos]

+1

This is a great summary that shows how multifaceted the
concept of FFT within Boost actually is --- or would be.

Taking on FFT in Boost would be challenging. The best FFT
implementations tend to be proprietary. They are also written in
language forms that lack the portability and clarity of interface
to which Boost strives.

At the same, the existing implementations are, insofar as
performance is concerned, hard to
beat.

Nonetheless, a portable FFT with Boost licensing would
be a very valuable resource for the community.

If there is ever a coalition of interested developers on this
topic, I might be in on it --- with whatever minimal contribution
I could offer.

Got it on the back burner.

Sincerely, Chris.