[SciPy-dev] fftpack revisited
Pearu
pearu at scipy.org
Mon Aug 19 17:19:13 CDT 2002
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
More information about the Scipy-dev
mailing list