[SciPy-dev] fft

[Chuck Harris]

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