[Numpy-discussion] BSD C port of FFTPACK incl. bluestein algorithm
Dag Sverre Seljebotn
Fri Nov 18 06:18:10 CST 2011
On 11/18/2011 12:58 PM, David Cournapeau wrote:
> On Fri, Nov 18, 2011 at 11:42 AM, Robert Kern<firstname.lastname@example.org> wrote:
>> On Fri, Nov 18, 2011 at 11:19, Dag Sverre Seljebotn
>> <email@example.com> wrote:
>>> I've been in touch with Martin Reinecke, author of the libpsht code for
>>> spherical harmonic transforms, about licensing issues.
>>> libpsht itself will remain under the GPL, but he is likely to release
>>> his C port of FFTPACK under BSD in the near future, as it is based on
>>> the public domain FFTPACK.
>>> I'm grateful for this change for my own purposes (allows releasing my
>>> own competing SHT library under the BSD) -- but it could perhaps be
>>> useful for NumPy or SciPy as well, depending on how complete the port
>>> is? E.g., perhaps make numpy.fft more complete (is the
>>> numpy.fft/scipy.fftpack split simply because of the Fortran dependency?).
>> It used to be the case that scipy.fftpack allowed one to build against
>> multiple different, usually faster, FFT libraries like FFTW. I think
>> we have backed away from that since the cost of maintaining the build
>> configuration for all of those different backends was so high. It's
>> worth noting that numpy.fft is already using a C translation of
>> FFTPACK. I'm not sure what the differences are between this
>> translation and Martin's.
Here's some more info forwarded from Martin:
- only FFTs are supported (no DCTs/DSTs)
- only double precision is supported (extension to single precision might
not be much work, though)
- both complex and real FFTs are supported
- real FFTs allow various storage schemes for the (half)complex frequency
domain data (classic FFTPACK scheme, FFTW or halfcomplex scheme,
uncompressed complex storage)
- precision of transforms involving large prime factors should be slightly
better than with original FFTPACK
- Bluestein's algorithm is automatically selected if considered profitable
- small accuracy self-testing code is provided.
Fairly complete interface documentation is available in Doxygen format.
I'll prepare a source package later in the afternoon and send it around.
> Having a Bluestein transformation alone would be worthwhile, as it
> would avoid the N^2 penalty for prime sizes.
> I am wondering about precision issues, though (when I tried
> implementing bluestein transforms on top of fftpack, it gave very bad
> results numerically-wise). A comparison with fftw would be good here.
Well, there's an indirect comparison: My SHT code currently uses FFTW3,
and it manages to agree with Reinecke's SHT code to a precision of
better than ~1e-12 for full SHTs. That includes several other sources of
(That's an average over several different-sized FFTs, of which half has
n=8192 and the other half all have different size, starting from 4 and
increasing up to 8192 in steps of 4 -- meaning prime factors on the
order of 1000).
I agree, a more direct comparison with FFTW would be good.
In more detail from the README:
I replaced the iterative sine and cosine calculations in radfg() and
17 by an exact calculation, which slightly improves the transform
18 real FFTs with lengths containing large prime factors.
More information about the NumPy-Discussion