[SciPy-dev] fftpack revisited

[eric jones]

By the way, interesting timing results on fftpack and fftw.  I'm glad to
see the penalty paid for using fftpack is either non existent or at
least not obvious. Any chance you can get the standard LAPACK to be as
fast as ATLAS? :-)  That sure would save a lot of wailing and gnashing
of teeth.

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