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

Question

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

>>> 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可以复制传递到的2d数组的输出：

>>> 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允许1D Fourier transformation使用matrix / 2d数组？

Answer 1

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

如果

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

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

可以通过以下方式验证

>>> 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。