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

Charles R Harris charlesr.harris@gmail....
Fri Nov 18 09:08:37 CST 2011


On Fri, Nov 18, 2011 at 5:18 AM, Dag Sverre Seljebotn <
d.s.seljebotn@astro.uio.no> wrote:

> On 11/18/2011 12:58 PM, David Cournapeau wrote:
> > On Fri, Nov 18, 2011 at 11:42 AM, Robert Kern<robert.kern@gmail.com>
>  wrote:
> >> On Fri, Nov 18, 2011 at 11:19, Dag Sverre Seljebotn
> >> <d.s.seljebotn@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.
>
>
The current numpy FFT code could use a review in any case, IIRC, there were
some possible integer related issues in the arguments. Having clean
documented C code would probably be an improvement and also help with
maintenance. I suspect there were iterative computations of sin/cos that
were spoiling the precision of the single precision versions (scipy), so
fixing those might also be useful.

Chuck
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://mail.scipy.org/pipermail/numpy-discussion/attachments/20111118/e1815615/attachment-0001.html 


More information about the NumPy-Discussion mailing list