[Numpy-discussion] BSD C port of FFTPACK incl. bluestein algorithm

[Dag Sverre Seljebotn]

On 11/18/2011 12:58 PM, David Cournapeau wrote:
> On Fri, Nov 18, 2011 at 11:42 AM, Robert Kern<robert.kern at gmail.com>  wrote:
>> On Fri, Nov 18, 2011 at 11:19, Dag Sverre Seljebotn
>> <d.s.seljebotn at astro.uio.no>  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.

Best regards,
   Martin
"""

>
> 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 
errors though.

(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 
radbg()
17	by an exact calculation, which slightly improves the transform 
accuracy for
18	real FFTs with lengths containing large prime factors.

Dag Sverre