numpy.fft numpy.fft2 和二维数组的 FFTW

Question

我正在尝试使用 C FFTW 库重现 numpy.fft.fft 和 numpy.fft.fft2 的输出。

>>> b
array([1, 2, 3, 4, 5, 6])
>>> type(b)
<class 'numpy.ndarray'>
>>> b.shape
(6,)
>>> np.fft.fft(b)
array([21.+0.j        , -3.+5.19615242j, -3.+1.73205081j, -3.+0.j        ,
       -3.-1.73205081j, -3.-5.19615242j])

可以通过以下方式获得此输出：

int N = 10;
double in[] = {1,2,3,4,5,6,0,0,0,0};
fftw_complex *out;
fftw_plan p;
out = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * (N/2 +1));
p = fftw_plan_dft_r2c_1d(6, in, out, FFTW_ESTIMATE);
fftw_execute(p);
fftw_destroy_plan(p);
fftw_free(out);

同样，可以重现传递给 numpy.fft.fft2 的二维数组的输出：

>>> a
array([[1, 2],
       [3, 4],
       [5, 6]])
>>> a.shape
(3, 2)
>>> np.fft.fft2(a)
array([[21.+0.j        , -3.+0.j        ],
       [-6.+3.46410162j,  0.+0.j        ],
       [-6.-3.46410162j,  0.+0.j        ]])

对应的C++代码为（只改一行）

p = fftw_plan_dft_r2c_2d(3, 2, in, out, FFTW_ESTIMATE);

我遇到了一个 Python 代码，它将 2d array 传递给 numpy.fft.fft

>>> a
array([[1, 2],
       [3, 4],
       [5, 6]])
>>> a.shape
(3, 2)
>>> np.fft.fft(a)
array([[ 3.+0.j, -1.+0.j],
       [ 7.+0.j, -1.+0.j],
       [11.+0.j, -1.+0.j]])

我正在尝试找出如何使用 FFTW APIs 来实现这一点。关于如何重现这个的任何线索？或者为什么 numpy 允许 matrix/2d 数组的 1D Fourier transformation？

Answer 1

为什么 NumPy 允许将二维数组传递给一维 FFT？目标是能够同时计算多个单独的一维信号的 FFT。

如果

>>> a = np.array([[1, 2], [3, 4], [5, 6]])
>>> A = np.fft.fft(a)

那么 A 的第一行将是 a 第一行的一维 FFT。 A的第二行是a的第二行的一维FFT，依此类推。

这可以通过

验证

>>> np.fft.fft(a[0, :])
array([ 3.+0.j, -1.+0.j])

>>> A[0, :]
array([ 3.+0.j, -1.+0.j])

或

>>> np.fft.fft(a[1, :])
array([ 7.+0.j, -1.+0.j])

>>> A[1, :]
array([ 7.+0.j, -1.+0.j])

要在 FFTW 中执行相同的操作，您可以针对不同的行多次执行计划，或者使用 fftw_plan fftw_plan_many_dft。

numpy.fft numpy.fft2 和二维数组的 FFTW

numpy.fft numpy.fft2 and FFTW for 2D arrays

python

numpy

fft

fftw