[SciPy-dev] fftpack revisited

[Travis Oliphant]

On Sun, 2002-08-18 at 12:30, Pearu Peterson wrote:
> 
> 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.
> 

This is interesting, anecdotally I've noticed fftpack to be "fast
enough" for me so that I haven't been using fftw at all.

What tests did you run?

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

No idea what could be happening.  

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

Actually I should not take credit here.  This is just the Numeric
interface contributed by several people.  I only moved it over to SciPy
once ND fft's were supported. 

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

I don't have a problem conceptually with using f2py for this.  In fact,
I would like to see most if not all of the interfaces (even the ones
that I did contribute :-) written with f2py)


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

This is a long-standing ambiguity with fft's.  You always have to check
what the definition is.  What does Matlab do?  We need an interface that
reduces surprise...


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

Yes, I think this is reasonable.

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

I think the axis keyword argument is necessary.  Sure it may not always
be used, but I very often do FFT's along one direction or another and I
know several other folks who code in MATLAB do as well.  The default
should be the one that gives the most speed, but we need to allow this
feature.  It would be O.K. to not have this feature in the low-level
interface, but the command exposed to the casual user in Python, must
have it.

I'm not opposed to a new f2py-generated interface to fftpack as long as
the command most users will use has the axis keyword.  Like I said it's
O.K. if this has to be a Python function with the low-level wrapper
having no axis keyword.

I think a new interface is better than a mixed f2py + current interface.

Just my opinion.

-Travis