1

我正在尝试在Python. 傅里叶位移定理允许通过以下方式将阵列平移一个子像素量： 1. 正向 FFT 阵列 2. 将阵列乘以傅里叶空间中的线性相位斜坡 3. 逆 FFT 阵列

该算法很容易在 python 中使用 numpy/scipy 实现，但对于 256**2 数组，每班次的速度非常慢（~10 毫秒）。我试图通过使用 scipy.weave.inline 直接从 python 调用 c 代码来加快速度。

但是，我在将复杂的 numpy 数组传递给 FFTW 时遇到了麻烦。c 代码如下所示：

    #include <fftw3.h>
    #include <stdlib.h>

    #define INVERSE +1
    #define FORWARD -1


    fftw_complex *i, *o;
    int n, m;
    fftw_plan pf, pi;
    #line 22 "test_scipy_weave.py" 

    i = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * xdim*ydim);
    o = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * xdim*ydim);

    pf = fftw_plan_dft_2d(xdim, ydim, i, o, -1, FFTW_PATIENT);
    pi = fftw_plan_dft_2d(xdim, ydim, o, i,  1, FFTW_PATIENT);

    # Copy data to fftw_complex array. How to use python arrays directly
    for (n=0; n<xdim;n++){
        for (m=0; m<ydim; m++){
            i[n*xdim+m][0]=a[n*xdim+m].real();
            i[n*xdim+m][1]=a[n*xdim+m].imag();
        }
    }

    fftw_execute(pf);

    /* Mult by linear phase ramp here */

    fftw_execute(pi);

    for (n=0; n<xdim;n++){
        for (m=0; m<ydim; m++){
            b[n*xdim+m] = std::complex<double>([in*xdim+m][0], i[n*xdim+m][1]);
        }
    }

    fftw_destroy_plan(p);

所以你可以看到我必须将存储在 numpy 数组“a”中的数据复制到 fftw_complex 数组“i”中。最后，我必须将结果“i”复制到输出 numpy 数组“b”中。直接在 fftw 中使用 numpy 数组“a”和“b”会更有效，但我无法让它工作。

有没有人知道如何让 fftw 直接使用复杂的 numpy 数组scipy.weave.inline？

谢谢

4

1 回答 1

0

根据fftw 手册，可以 import complex.hbefore fftw.h，这样可以保证fftw_complex对应原生 C 数据类型。我很确定 numpy 数据类型也保证（或在实践中可能）与本机 C 数据类型兼容。

在这种情况下，您可以访问指向数组数据的指针a.data_as(ctypes.c_void_p)。不幸的是 ctypes 不能识别复杂类型，但希望转换为 void 指针可以解决问题。

执行此操作时，您必须注意您的数组a以 C 连续方式存储，由order='C'创建数组时的参数指定。

于 2012-10-25T07:11:49.230 回答