[SciPy-dev] fftpack revisited

[Pearu]

On 19 Aug 2002, Travis Oliphant wrote:

> What tests did you run?

Attached are test codes. 
Usage:
1) Unpack fftw2.tgz and cd fftw2
2) To build, run 
 python ./setup_fftw.py build --build-platlib=.
3) To test, run
 python bench_fft.py

Here are the results that I get:

n= 1000
Testing  <fortran object at 0x819af48>  ...  0.6153 seconds
Testing  <function fft at 0x835ea84>  ...  0.2199 seconds
Testing  <fortran object at 0x819af48>  ...  0.5865 seconds
Testing  <function fft at 0x835ea84>  ...  0.2179 seconds
n= 1024
Testing  <fortran object at 0x819af48>  ...  0.5279 seconds
Testing  <function fft at 0x835ea84>  ...  0.2158 seconds
Testing  <fortran object at 0x819af48>  ...  0.5033 seconds
Testing  <function fft at 0x835ea84>  ...  0.2146 seconds

where
  <fortran object at 0x819af48> refers to fftw routine
  <function fft at 0x835ea84> refers to scipy.fftpack.fft function

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

In matlab:
++++++++++++++++++++++++
>> help fft

 FFT Discrete Fourier transform.
    FFT(X) is the discrete Fourier transform (DFT) of vector X.  For
    matrices, the FFT operation is applied to each column. For N-D
    arrays, the FFT operation operates on the first non-singleton
    dimension.
 
    FFT(X,N) is the N-point FFT, padded with zeros if X has less
    than N points and truncated if it has more.
 
    FFT(X,[],DIM) or FFT(X,N,DIM) applies the FFT operation across the
    dimension DIM.
<snip>    

>> help fftn

 FFTN N-dimensional discrete Fourier Transform.
    FFTN(X) returns the N-dimensional discrete Fourier transform
    of the N-D array X.  If X is a vector, the output will have
    the same orientation.
 
    FFTN(X,SIZ) pads X so that its size vector is SIZ before
    performing the transform.  If any element of SIZ is smaller
    than the corresponding dimension of X, then X will be cropped
    in that dimension.
++++++++++++++++++++++++++++++++++

In addition, matlab has functions IFFT, FFT2, IFFT2, IFFTN.
For matlab functions, ifft(fft(x))=x always holds.

Also I think the current fftpack provides too many fft functions (see the 
list in my previous mail) compared to what Matlab provides. Optimization 
decisions which Fortran fftpack  functions to call for given input data 
can be partly made in fft routines. I'll see if some of the routines can 
be dropped if their functionality can be implemented in some more generic 
fft routine. Basically, here follows what we have:

* fft(x),ifft(x) should accept arbitrary type sequences (int,real,complex)
and they always return complex arrays of length len(x). Also 
fft(ifft(x))=x must hold. Internally, the most optimal fftpack routine is 
chosen based on the type of x.

* rfft(x) accepts only real sequences. Raises an exception if x is 
complex. Returns a complex array of length len(x)/2. irfft(x) accepts any 
sequences, returns a real array of length 2*len(x) (I have assumed that 
len(x) is even). irfft(rfft(x),len(x)) = x holds.

* rrfft(x),irrfft(x) accept only real sequences. Both return real arrays 
of length len(x). rrfft(irrfft(x))=x holds.

* hfft(x),ihfft(x) -- opposite to rfft/irfft functions.

And then there are 2D and N-D fft routines.

Remarks:

1) What shall we decide about the signatures of fft functions?
Matlab fft has basically the following signature

  fft(x[,n[,axis]]) -> y

Shall use the same? Initially I was thinking of

  fft(x[,axis[,n]]) -> y

with the possibility of dropping the n argument. Though Matlab functions 
do truncating and zero-padding if x.shape[axis]!=n, I find it dangerous: 
if one uses incorrect n, then fft function still succeeds. As a result, 
the calculations will be incorrect (users may not even notice that) or 
some exception is raise later, far from the fft call. Such cases are 
difficult to debug and therefore I would prefer if fft takes input with 
the correct size, truncating and zero-padding should be done explicitely 
by the user, scipy can only provide the corresponding functions.
What do you think?

2) Notice that fft(x)[:len(x)/2] == rfft(x). Do we acctually need the 
rfft(x) function exposed in scipy?

Pearu
-------------- next part --------------
A non-text attachment was scrubbed...
Name: fftw2.tgz
Type: application/x-gzip
Size: 2668 bytes
Desc: 
Url : http://projects.scipy.org/pipermail/scipy-dev/attachments/20020819/86bea1f4/attachment.tgz