[Numpy-discussion] FFTS for numpy's FFTs (was: Re: Choosing between NumPy and SciPy functions)

[Sturla Molden]

On 29/10/14 10:48, Eelco Hoogendoorn wrote:

 > Id rather have us discuss how to facilitate the integration of
 > as many possible fft libraries with numpy behind a maximally uniform
 > interface, rather than having us debate which fft library is 'best'.

I am happy with the NumPy interface. There are minor differences between 
the SciPy and NumPy FFT interfaces (e.g. for rfft, see below). 
Personally I prefer the NumPy interface because it makes it easier to 
map Fourier coeffs to frequencies.

One thing we could do, without too much hassle, is to use FFTs from MKL 
or Accelerate (vDSP) if we link with these libararies for BLAS/LAPACK.

MKL has an API compatible with FFTW, so FFTW and MKL can be supported 
with the same C code. FFTW and MKL also have a Fortran 77 API which we 
could wrap with f2py (no Fortran compiler are needed). It is actually 
possible to use the FFTs in FFTW and MKL from Python without any C 
coding at all. We just need to add a Python interface on top of the f2py 
wrappers, which is similar to what we do for scipy.linalg.

The FFTs in Accelerate have a very simple C interface, but only support 
power-of-two array sizes, so we would need to use them with Bluestein's 
algorithm. Again, because of their simplicity, it is possible to wrap 
these FFT functions with f2py.

We cannot bundle NumPy or SciPy binaries with FFTW due to GPL [*], but 
as I understand it we already have permission from Intel to bundle 
binary wheels linked with MKL. Accelerate is a system library, so that 
does not pose a license problem.

[*] Actually, we could, but the binaries would be tainted with a viral 
license.



 >>> a = np.random.rand(8)

 >>> scipy.fftpack.rfft(a)[:,None]

array([[ 3.47756851],
        [-0.45869926],
        [-0.21730867],
        [-0.43763425],
        [-0.67338213],
        [-0.28799   ],
        [ 0.17321793],
        [-0.31514119]])

 >>> np.fft.rfft(a)[:,None]

array([[ 3.47756851+0.j        ],
        [-0.45869926-0.21730867j],
        [-0.43763425-0.67338213j],
        [-0.28799000+0.17321793j],
        [-0.31514119+0.j        ]])





Sturla