[SciPy-dev] fftpack revisited

Pearu Peterson pearu at cens.ioc.ee
Sun Aug 18 17:19:17 CDT 2002


On Sun, 18 Aug 2002, eric jones wrote:

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

Ok then. Now that I look it more closely, adding axis support is pretty
much paste-copy from the current fftpack.

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

(1) I'll create the testing suite while proceeding with fftpack2
(2) Actually, I can create this loop also in C while still using
f2py (sometimes f2py strikes me with its flexibility ;).

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

How about
  rrfft,irrfft
? Read it as 'real real fft' -- the second 'real' for the fact that input
is assumed to be real and the first 'real' for the fact that the returned
array is really real ;-)

Hmm, currently fftpack provides the following set of fft functions:
 fft,ifft
 rfft,irfft
 hfft,ihfft
 fftn,ifftn
 rfftn,irfftn
 fft2,ifft2
 rfft2,irfft2
These would be the additions:
 rrfft,irrfft
 rrfftn,irrfftn
 rhfft,irhfft
 rrfft2,irrfft2 
+ various cos/sin transforms. Is this pattern acceptable?

In addition, I would like to introduce the following functions to
fftpack, all use fft extensively:

y = diff(x, k=1, period=2*pi) - k-th derivative (integral) of a periodic
    sequence. Exact for trigonometric polynomials.
y = hilbert(x) - Hilbert transform (and its inverse as ihilbert)
y = tilbert(x,a=1) - Tilbert transform (and its inverse as itilbert)
y = difftilbert(x,k=1,a=1,period=2*pi) - optimized diff(tilbert(..),..).

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

I'll keep the axis arguments.

There is only one thing in the interface that I would like to discuss. For
example, fftpack.fft has three arguments:
  fft(a, n=None, axis=-1)
The second argument n stands for the size of the Fourier transform
(matrix). It seems to me that this argument could be dropped as this
information already is contained in a.shape and axis:
  n = a.shape[axis]
In fftpack._raw_fft this is a special case and, in general, the following
rules are applied: 
*  If n>a.shape[axis] then a is padded with zeros.
*  If n<a.shape[axis] then a is trunctated. 
Both cases create a new array that is non-contiguous (except in
special cases). So, in total there will be two copies made of input data
before passing it to Fortran. 
I wonder if we could move this functionality out from the fft functions?
I think that the padding functionality is already covered by the
fftpack.zeropad function and truncating can be done by slicing.

This brings me to a general note: fft and relatives are fundamental tools
that should perform specific tasks. When using them I'll assume that they
perform their tasks in an optimal way. As an user, I am not interested in
all kinds of additional bells and whistles that these tools could
carry. They just make me wonder what other magic is done in these
functions, that is, whether they can be optimized away for my specific
problems.

Pearu




More information about the Scipy-dev mailing list