螺紋FFT在numpy的/ SciPy的Enthought的Python

Question

快速傅立葉變換（FFT）不帶螺紋的。 Enthought Python與英特爾MKL數字庫一起提供，該庫可以線程化FFT。如何獲得這些例程？

Answer 1

新的和改進的版本，其在處理輸入和輸出數組任意進步。 默認情況下，這一個現在不在位並創建一個新的數組。 它模仿Numpy FFT例程，只是它具有不同的標準化。

''' Wrapper to MKL FFT routines ''' 

import numpy as _np 
import ctypes as _ctypes 

mkl = _ctypes.cdll.LoadLibrary("mk2_rt.dll") 
_DFTI_COMPLEX = _ctypes.c_int(32) 
_DFTI_DOUBLE = _ctypes.c_int(36) 
_DFTI_PLACEMENT = _ctypes.c_int(11) 
_DFTI_NOT_INPLACE = _ctypes.c_int(44) 
_DFTI_INPUT_STRIDES = _ctypes.c_int(12) 
_DFTI_OUTPUT_STRIDES = _ctypes.c_int(13) 

def fft2(a, out=None): 
    ''' 
Forward two-dimensional double-precision complex-complex FFT. 
Uses the Intel MKL libraries distributed with Enthought Python. 
Normalisation is different from Numpy! 
By default, allocates new memory like 'a' for output data. 
Returns the array containing output data. 
''' 

    assert a.dtype == _np.complex128 
    assert len(a.shape) == 2 

    inplace = False 

    if out is a: 
     inplace = True 

    elif out is not None: 
     assert out.dtype == _np.complex128 
     assert a.shape == out.shape 
     assert not _np.may_share_memory(a, out) 

    else: 
     out = _np.empty_like(a) 

    Desc_Handle = _ctypes.c_void_p(0) 
    dims = (_ctypes.c_int*2)(*a.shape) 

    mkl.DftiCreateDescriptor(_ctypes.byref(Desc_Handle), _DFTI_DOUBLE, _DFTI_COMPLEX, _ctypes.c_int(2), dims) 

    #Set input strides if necessary 
    if not a.flags['C_CONTIGUOUS']: 
     in_strides = (_ctypes.c_int*3)(0, a.strides[0]/16, a.strides[1]/16) 
     mkl.DftiSetValue(Desc_Handle, _DFTI_INPUT_STRIDES, _ctypes.byref(in_strides)) 

    if inplace: 
     #Inplace FFT 
     mkl.DftiCommitDescriptor(Desc_Handle) 
     mkl.DftiComputeForward(Desc_Handle, a.ctypes.data_as(_ctypes.c_void_p)) 

    else: 
     #Not-inplace FFT 
     mkl.DftiSetValue(Desc_Handle, _DFTI_PLACEMENT, _DFTI_NOT_INPLACE) 

     #Set output strides if necessary 
     if not out.flags['C_CONTIGUOUS']: 
      out_strides = (_ctypes.c_int*3)(0, out.strides[0]/16, out.strides[1]/16) 
      mkl.DftiSetValue(Desc_Handle, _DFTI_OUTPUT_STRIDES, _ctypes.byref(out_strides)) 

     mkl.DftiCommitDescriptor(Desc_Handle) 
     mkl.DftiComputeForward(Desc_Handle, a.ctypes.data_as(_ctypes.c_void_p), out.ctypes.data_as(_ctypes.c_void_p)) 

    mkl.DftiFreeDescriptor(_ctypes.byref(Desc_Handle)) 

    return out 

def ifft2(a, out=None): 
    ''' 
Backward two-dimensional double-precision complex-complex FFT. 
Uses the Intel MKL libraries distributed with Enthought Python. 
Normalisation is different from Numpy! 
By default, allocates new memory like 'a' for output data. 
Returns the array containing output data. 
''' 

    assert a.dtype == _np.complex128 
    assert len(a.shape) == 2 

    inplace = False 

    if out is a: 
     inplace = True 

    elif out is not None: 
     assert out.dtype == _np.complex128 
     assert a.shape == out.shape 
     assert not _np.may_share_memory(a, out) 

    else: 
     out = _np.empty_like(a) 

    Desc_Handle = _ctypes.c_void_p(0) 
    dims = (_ctypes.c_int*2)(*a.shape) 

    mkl.DftiCreateDescriptor(_ctypes.byref(Desc_Handle), _DFTI_DOUBLE, _DFTI_COMPLEX, _ctypes.c_int(2), dims) 

    #Set input strides if necessary 
    if not a.flags['C_CONTIGUOUS']: 
     in_strides = (_ctypes.c_int*3)(0, a.strides[0]/16, a.strides[1]/16) 
     mkl.DftiSetValue(Desc_Handle, _DFTI_INPUT_STRIDES, _ctypes.byref(in_strides)) 

    if inplace: 
     #Inplace FFT 
     mkl.DftiCommitDescriptor(Desc_Handle) 
     mkl.DftiComputeBackward(Desc_Handle, a.ctypes.data_as(_ctypes.c_void_p)) 

    else: 
     #Not-inplace FFT 
     mkl.DftiSetValue(Desc_Handle, _DFTI_PLACEMENT, _DFTI_NOT_INPLACE) 

     #Set output strides if necessary 
     if not out.flags['C_CONTIGUOUS']: 
      out_strides = (_ctypes.c_int*3)(0, out.strides[0]/16, out.strides[1]/16) 
      mkl.DftiSetValue(Desc_Handle, _DFTI_OUTPUT_STRIDES, _ctypes.byref(out_strides)) 

     mkl.DftiCommitDescriptor(Desc_Handle) 
     mkl.DftiComputeBackward(Desc_Handle, a.ctypes.data_as(_ctypes.c_void_p), out.ctypes.data_as(_ctypes.c_void_p)) 

    mkl.DftiFreeDescriptor(_ctypes.byref(Desc_Handle)) 

    return out 

Answer 2

下面的代碼在Windows 7旗艦版64位對我的作品與Enthought 7.3-1（64位）。我沒有對它進行基準測試，但它肯定會同時使用所有內核，而不是僅使用一個。

from ctypes import * 

class Mkl_Fft: 
    c_double_p = POINTER(c_double) 

    def __init__(self,num_threads=8): 
     self.dfti = cdll.LoadLibrary("mk2_rt.dll") 
     self.dfti.MKL_Set_Num_Threads(num_threads) 
     self.Create = self.dfti.DftiCreateDescriptor_d_md 
     self.Commit = self.dfti.DftiCommitDescriptor 
     self.ComputeForward = self.dfti.DftiComputeForward 

    def fft(self,a): 
     Desc_Handle = c_void_p(0) 
     dims = (c_int*2)(*a.shape) 
     DFTI_COMPLEX = c_int(32) 
     rank = 2 

     self.Create(byref(Desc_Handle), DFTI_COMPLEX, rank, dims) 
     self.Commit(Desc_Handle) 
     self.ComputeForward(Desc_Handle, a.ctypes.data_as(self.c_double_p)) 

用法：

import numpy as np 
a = np.ones((32,32), dtype = complex128) 
fft = Mkl_Fft() 
fft.fft(a) 

Answer 3

我原來的答覆的清潔器版本如下：

from ctypes import * 

mkl = cdll.LoadLibrary("mk2_rt.dll") 
c_double_p = POINTER(c_double) 
DFTI_COMPLEX = c_int(32) 
DFTI_DOUBLE = c_int(36) 

def fft2(a): 
    Desc_Handle = c_void_p(0) 
    dims = (c_int*2)(*a.shape) 

    mkl.DftiCreateDescriptor(byref(Desc_Handle), DFTI_DOUBLE, DFTI_COMPLEX, 2, dims) 
    mkl.DftiCommitDescriptor(Desc_Handle) 
    mkl.DftiComputeForward(Desc_Handle, a.ctypes.data_as(c_void_p)) 
    mkl.DftiFreeDescriptor(byref(Desc_Handle)) 

    return a 

def ifft2(a): 
    Desc_Handle = c_void_p(0) 
    dims = (c_int*2)(*a.shape) 

    mkl.DftiCreateDescriptor(byref(Desc_Handle), DFTI_DOUBLE, DFTI_COMPLEX, 2, dims) 
    mkl.DftiCommitDescriptor(Desc_Handle) 
    mkl.DftiComputeBackward(Desc_Handle, a.ctypes.data_as(c_void_p)) 
    mkl.DftiFreeDescriptor(byref(Desc_Handle)) 

    return a