我正在尝试在 Python
中实现基于 FFT 的子像素移位(转换)算法。傅立叶移位定理允许通过以下方式将阵列平移一个子像素数量:
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 直接在 scipy.weave.inline
中使用复杂的 numpy 数组?
谢谢
最佳答案
根据fftw manual ,你可以在fftw.h
之前导入complex.h
,这样可以保证fftw_complex
对应原生的C数据类型。我很确定 numpy 数据类型也保证(或实际上可能)与 native C 数据类型兼容。
在这种情况下,您可以访问指向数组数据的指针作为 a.data_as(ctypes.c_void_p)
。不幸的是,ctypes 无法识别复杂类型,但希望转换为 void 指针可以解决问题。
执行此操作时,必须注意数组 a
以 C 连续方式存储,在创建时由参数 order='C'
指定数组。
关于python - 直接在 scipy.weave.inline 中对复杂的 numpy 数组进行 FFTW3,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/13055090/