[SciPy-dev] fftpack revisited

Pearu Peterson pearu at cens.ioc.ee
Sun Aug 18 13:30:07 CDT 2002

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

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.


More information about the Scipy-dev mailing list