[boost] ffts

From: Jason D Schmidt (jd.schmidt_at_[hidden])
Date: 2002-09-21 18:48:25

On this list, I have discussed submitting my fft code to Boost.  Some of
the feedback I have received was:

The interface should be generic:
    should be able to do arrays whose size is a power of 2 & not a power
of 2
    should be able to transform any container generically
    should be able to do real or complex data
    if feasible, it might be useful to handle integer & floating point
data

Here is the current state of my fft code with respect to the above
points:

The prototype is:
template <const size_t cPoints, class value_type>
std::valarray<std::complex<value_type> > fft(const
std::valarray<std::complex<value_type> >& data).

When one performs an n-point fft, the first thing that has to be done is
to calculate the bit-reversal of the indices (it's the critical part of
the algorithm).  The bit-reversed indices are different for each
different fft size.  With my code, a 128-point fft is a different
function from a 256-point function (because the size is a template
parameter).  The bit-reversed indices are stored in a static valarray
object, which is initialized the first time the fft function is called,
and it's never modified again.  Consequently, if I call my 256-point fft
function 10 times, the bit-reversal part of the algorithm happens only in
the first call, and the 9 subsequent calls only have to do the second
part of the fft algorithm.  Thus, on average, my function is very fast
because it only does one part of the algorithm (all but the first time).

The drawback to this design is that the size of the fft has to be known
at compile time (because it's a template parameter), and it makes other
computations (psd, correlation, convolution, etc.) that depend on the fft
very cumbersome.  Plus, the fft size has to be a power of 2.  However, I
could pad the input data with zeros so that it has a size which is a
power of 2.

Also, the above prototype shows that my fft function can't handle any
container generically.  I implemented it for only the std::valarray
because valarray has functionality that makes it much faster than any
other standard container.  For example, inside the fft function if I want
to write code like

static const std::valarray<size_t> bit_reversed
indices(do_bit_reversal<cPoints>());
const std::valarray<value_type> bit_reversed_data[bit_reversed_indices] =
data;

there's no container-generic way to do it.  I can't assume that every
container will have an operator[] that takes a generic container of
size_t variables.  In fact, valarray is the only standard container with
that capability (along with std::slice & std::gslice capabilities).

Ok, you might say that I should write my fft implementation as if it's
operating on an STL container or a range of generic iterators.  However,
then it wouldn't work with valarray, because valarray doesn't have
iterators, and it wouldn't have the performance gain that it makes with
valarray's operator[] and slicing capabilities.

You can also see from the prototype above that my fft function takes a
complex argument and returns complex values.  I have a function that
transforms real data and returns complex values, too.  As long as both
functions have the same prototype, then my code can transform both real &
complex data.  I just always make sure to return complex values because,
in general, the return value of an fft is complex, and you can't overload
functions by return type.

As far as handling integer or complex<int> input, it's certainly possible
to pass in int or complex<int> arguments.  I would just need to do a
conversion to float/double (to perform the computation properly) or both
copy and convert, which would probably degrade the performance.  I would,
however, still want complex<float/double> as the return value, because
there's no way to tell if the fft of any given integer data will result
in floats or integers.

To conclude, I fully agree with the feedback I listed at the beginning of
this message.  I have given my reasons for not implementing the generic
interface that I/we might like.  Many times the reasons were, "I don't
know of a way to do that," or "It would degrade the performance."  I
welcome any comments that would help me find ways around roadblocks and
how to maintain good performance at the same time.

Jason Schmidt

________________________________________________________________
GET INTERNET ACCESS FROM JUNO!
Juno offers FREE or PREMIUM Internet access for less!
Join Juno today!  For your FREE software, visit:
http://dl.www.juno.com/get/web/.