[SciPy-dev] OS X fft problem is serious
Tom Loredo
loredo at astro.cornell.edu
Thu Jan 12 07:24:32 CST 2006
Hi folks,
I hate to whine, especially when I don't know how to go about
solving the problem I'm whining about, but I don't think
it's good "advertising" to base the first official release of
new SciPy on a revision with a badly broken fft exposed in
the top-level namespace on widely-used platform.
Chris Fonnesbeck has verified that the problem exists also on
his Tiger build. This is not a quirk of my Panther installation.
Something is really broken.
The ipython session logged below puts it starkly. There is
a simple sine-like input array: [0 1 0 -1 0]. fft and then
ifft with numpy and scipy. numpy gives the array back (line 4).
scipy gives nonsense (line 7).
There is a clue to what is going on. The ffts, though not
exactly equal, are essentially equal before the ifft (line 6).
But *after* taking the ifft, scipy changes the input fft
in-place (compare line 10 and line 9). I believe the
default behavior is no overwrite. Perhaps this points toward
a solution. Perhaps the problem is not with SciPy but with
FFTW2; either way, I think we need to figure it out, preferably
before announcing the new release.
I wish I could help more with this, but I'll be at a workshop
for the next ~2 weeks. I'll be glad to help when I'm free again.
-Tom
In [1]: import numpy, scipy
Overwriting fft=<function fft at 0x74d4f0> from scipy.fftpack.basic (was <function fft at 0x73c4f0> from
numpy.dft.fftpack)
Overwriting ifft=<function ifft at 0x74d530> from scipy.fftpack.basic (was <function inverse_fft at 0x73c530> from
numpy.dft.fftpack)
In [2]: a=numpy.array([0.,1.,0.,-1.,0.])
In [3]: nfa=numpy.fft(a)
In [4]: print numpy.ifft(nfa).real
[ 0.00000000e+00 1.00000000e+00 9.76996262e-16 -1.00000000e+00
-9.76996262e-16]
In [5]: sfa=scipy.fft(a)
In [6]: numpy.allclose(nfa,sfa)
Out[6]: True
In [7]: print scipy.ifft(sfa).real
[ 0. 0.5 -0.5 -0.5 0.5]
In [8]: numpy.allclose(nfa,sfa)
Out[8]: False
In [9]: print nfa
[ 0. +0.j 1.11803399-1.53884177j -1.11803399+0.36327126j
-1.11803399-0.36327126j 1.11803399+1.53884177j]
In [10]: print sfa
[ 0. +0.j 0.5+0.j -0.5-0.j -0.5+0.j 0.5-0.j]
In [11]: numpy.__version__
Out[11]: '0.9.3.1863'
In [12]: scipy.__version__
Out[12]: '0.4.4.1544'
-------------------------------------------------
This mail sent through IMP: http://horde.org/imp/
More information about the Scipy-dev
mailing list