[SciPy-dev] fft
Chuck Harris
Chuck.Harris at sdl.usu.edu
Mon Feb 10 16:54:29 CST 2003
> -----Original Message-----
> From: Pearu Peterson [mailto:pearu at cens.ioc.ee]
> Sent: Monday, February 10, 2003 2:06 PM
> To: scipy-dev at scipy.net
> Subject: Re: [SciPy-dev] fft
>
>
>
> On Mon, 10 Feb 2003, Chuck Harris wrote:
>
> > I was playing around with the fft programs in fftpack and fftw and
> > noticed that the real transforms return the coefficients in
> different
> > formats. If rk,ik are the real and imaginary components of the k'th
> > Fourier coefficient, then
> >
> > fftpack: [r0,r1,i1,.....,rm]
> > where rm is absent in odd length transforms
> >
> > fftw: [r0,r1,...,i2,i1]
> > where again rm is absent in odd length transforms
>
> There is no need to use fftw package anymore. Current scipy.fftpack
> uses fftw libraries if available.
>
> Btw, scipy.fftpack uses the same convention for real transforms
> that is used in fftpack from netlib.
The main virtue of libraries is that someone else already wrote
them :o) I suspect that they chose this format so they could do
in place transforms, back in the days when memory mattered.
>
> > Besides the incompatibility of the forms, neither is of much use
> > in Numeric, where the coefficients would preferably be exposed as an
> > array of complex. Of course they do have the virtue that the number
> > of (real) components is exactly equal to the number of points
> > transformed.
>
> There are better reasons to keep real and imaginary parts separate
> than the above. For instance, in pseudo-spectral methods for
> integrating
> PDEs it is essential to minimize floating point operations that, in
> general, are the main source of numerical noise causing
> possible stability problems for integrators.
coef.real, coef.imag. Same algorithm, just with rearranged output.
>
> > However, I think the best solution would be to return
> > a complex array [r0+i0j,r1+i1j,...]. In this case i0==0 and when
> > the number of samples is even im==0 also. This means that the length
> > of the returned complex array is n//2 + 1, which may seem a bit
> > inconvenient, but the result is then easy to use in Numeric for
> > such things as convolution, amplitude, and phase.
>
> You can safely use scipy.fftpack.fft function if complex result is
> needed, just take the first half of the resulting array. Note that
> scipy.fftpack.fft takes into account if the input array is
> real and then
> a more efficient algorithm is used for calculating FT.
Sure, I wrote my first split-radix algorithm 32 years ago in
Algol 68 (blush), so I'm not unfamiliar with these things. The
virtue of the real transforms is that they *are* more efficient,
that is why it would be nice to have the output in a compatible
format.
>
> > By the way, what happened to the other fft package that was being
> > discussed?
>
> What other fft package?
I must be confused, but weren't you folks benchmarking another package?
Chuck
More information about the Scipy-dev
mailing list