[SciPy-dev] fftpack revisited

[eric jones]

Hey Pearu,

Here are my thoughts:

1.  Given a choice, I like f2py wrappers over hand wrappers because they
are easier to maintain.
2.  I think the axis argument is important.  It allows users to choose
how they lay out their data instead of SciPy defining how they should do
it.  The axis should default to -1 so that the default behavior operates
on contiguous memory without copying for efficiency.  This is the
default used throughout the SciPy library.

As far as whether looping over the data sets in C or Python has much
impact on speed, a simple test would be to rewrite one of Travis'
functions in Python (such as fft) and then compare the two.  If Python
is within a few percent of the wrapped fft when testing "normal sized"
2d and 3d arrays, then the benefit of having the code in C is not to
great.

So, I'd be happy to see fftpack become an f2py generated library as long
as (1) we're very sure that the new version is at least as robust as the
old version, and (2) there is very little penalty for using Python
looping.  Checking the repository, I didn't see any unit tests for the
fftpack library which will make verifying (1) more difficult.

As far as your new fftr function, it sounds like a good addition to me,
but I don't like the name that well.  It is too similar to the original
function name (rfft and fftr) and could cause confusion.  I'd go for
rfft2, except that it sounds like it is doing a 2d rfft.  No other
suggestions come immediately to mind though.

On exposing more transforms from fftpack, that sounds like a good idea
that is easily done with f2py than by hand.  I wonder though, if the
standard jpeg library doesn't have a blazing fast cosine transform
function that we might use.  Perhaps it is not worth the trouble of
adding another set of source code...

Summary: fftpack2 sounds like a good thing to pursue, but it would need
to support axis arguments if it is to become a replacement for fftpack.
Is there any other functionality you were thinking of dropping?  These
would need to be discussed also.

Thanks,
eric

> -----Original Message-----
> From: scipy-dev-admin at scipy.net [mailto:scipy-dev-admin at scipy.net] On
> Behalf Of Pearu Peterson
> Sent: Sunday, August 18, 2002 1:30 PM
> To: scipy-dev at scipy.org
> Subject: [SciPy-dev] fftpack revisited
> 
> 
> Hi scipy devels,
> 
> Due to emerged need for the fft algorithm in my work, I have returned
to
> the question whether the fftpack module provided by scipy is optimal.
> 
> Leaving aside license issues, I considered fftw again. Running some
> performance tests with fftpack and fftw, I found to my surprise
> that fftpack is actually 2-2.5 times faster (!) than fftw (could
someone
> verify this as well?). I checked fftw site for comparison and there
> fftpack was found to be about 2-3 times slower than fftw.
> 
> I don't have a good explanation for this difference, only few notes:
> 1) fftw site timings were obtained with rather old gcc-2.7.2. May be
gcc
> has become faster since then (I am using gcc-2.95.4).
> 2) scipy uses rather aggressive Fortran compiler flags, some of them
was
> used in fftw site as well. So, compiler flags should not be the
> (only) reason.
> 3) fftw tests were run on PII with 200MHz and 300MHz and Redhat 5.0, I
am
> running them on PII with 400MHz and Debian Woody. Other parameters
seem to
> be the same (including compiler flags for building fftw).
> 
> 
> Next, I needed to implement some transforms based on DFT (like
> differentiation of periodic sequences, hilbert transforms, tilbert
> transforms) that has to be very efficient and stable (they will be
called
> tens of thousands times on rather large data set).
> 
> The current interface of fftpack is rather complete and general,
> Travis O. has made remarkable job in creating it, but it turns out to
be
> not straightforward enough for my needs.
> 
> Therefore I, have created another interface to the Fortran fftpack
> library using f2py, let's call it fftpack2 for now. Currently it
exposes
> the same functions to Python as fftpack written by Travis O. but
differs
> from the following aspects
> 
> 1) y = fftpack.rfft(x) (and its inverse) is replaced with
> y = fftpack2.fftr(x) that returns real data array of length len(x):
>   y = [y(0),Re(y(1)),Im(y(1)),...]
> 
> The differences with the fftpack.rfft is that fftpack.rfft returns a
> complex array of length
>   len(x)/2     if len(x) is even
>   (len(x)-1)/2 if len(x) is odd
> and reads as follows
>   y = [y(0),Re(y(1))+1j*Im(y(1)),...]
> 
> To my opinion, fftpack2.fftr has the following advantages:
>   a) the property ifft(fft(x))=x holds what ever is the length of x.
>      With fftpack.rfft one has to specify the length when doing
inverse,
>      i.e. ifft(fft(x),len(x))=x holds, in general. I find it as a
possible
>      source of confusion and bugs.
>   b) there is no real->complex->real conversion in fftpack2.fftr that
>      makes the interface more straightforward (good for performance,
>      crucial for time critical algorithms).
> However, in certain situations fftpack.rfft has an advantage (and this
is
> the reason why I introduced different name for this function, both
rfft
> and fftr should be available).
> 
> 2) fftpack2 does not support axis keyword argument. The purpose of
axis
> argument is that FFT can be applied in one C loop to the rows of an
> array. This feature is missing in the original Fortran fftpack
library. On
> the other hand, Matlab's fft has this feature (where FFT is applied to
the
> columns of a matrix, though).
> 
> Why I have not added this feature to fftpack2? First, I am not sure
how
> often this feature will be used so that code complexity would be
> justified. In SciPy, it used only in signal/signaltools.py for hilbert
> transform where array transpose could be used instead. Second, I don't
see
> any remarkable performance gain for doing the loop in C: most of the
time
> will be spent probably in doing FFT or converting the input array to a
> contiguous one. So, the loop could be easily done in Python without
> too much of performance loss.
> 
> So, my questions are:
> Is the axis keyword argument actually used among scipy users? Should I
add
> this feature also to fftpack2 because Matlab has one? (But notice that
> Matlab uses different convention from scipy, and so Matlab users will
be
> confused anyway).
> 
> Right now I am in a position of keeping fftpack functions as simple as
> possible and add additional features only when they are explicitly
asked
> for.
> 
> 3) Finally, fftpack2 is easier to maintain due to using f2py. I have
found
> number bugs in fftpack causing either infinite loops or segfaults and
> fixing them can be a tedious task.
> Fortran fftpack provides also other transforms like sin,cos,etc that
> could be exposed to scipy. This would be really easy when using f2py.
> The alternative of manual wrapping and subsequent maintainance pain is
> not very attractive for me.
> 
> 
> In conclusion, I am a bit in the middle of deciding whether ...
> * I should proceed with fftpack2? It has an advantage of consisting
> simpler code base and also maintaining it will be easier.
> But I would like to drop some of the features (like axis argument)
from
> fftpack if there is no need for it.
> * Or should I try to fix the current fftpack bugs, extend it but
> keeping all the current features? It has a disadvantage that
> the fftpack will contain two types of C codes, one manually written
and
> other wrapped with f2py. This setup will make the fftpack module
seemingly
> more complex than necessary: I find it both from the usage as well
from
> the maintainance point of views.
> 
> So, what do you think of all that? I'll appreciate any comments that
> would help to solve my dilemma.
> 
> Thanks,
> 	Pearu
> 
> _______________________________________________
> Scipy-dev mailing list
> Scipy-dev at scipy.net
> http://www.scipy.net/mailman/listinfo/scipy-dev