python - 我应该如何解释 numpy.fft.rfft2 的输出？

Question

显然，rfft2 函数只是计算输入矩阵的离散 fft。但是，如何解释给定的输出索引？给定输出索引，我在看哪个傅里叶系数？
我对输出的大小特别困惑。对于 n x n 矩阵，输出似乎是 n x (n/2)+1 矩阵(对于偶数 n)。为什么方阵最终会进行非方阵傅立叶变换？

Answer 1

numpy.fft.rfft2 的输出只是标准二维 FFT 的左半部分(加一列)，由 numpy.fft.fft2 计算得出.不需要 rfft2提供结果的右半部分，因为实数数组的 FFT 具有 natural and simple symmetry ，因此可以使用该对称性从左半部分导出完整 FFT 的右半部分。

这里举个例子来说明。首先，为了便于重现和查看，我将设置 NumPy 的随机状态和打印选项:

In [1]: import numpy as np

In [2]: np.set_printoptions(precision=3, suppress=True, linewidth=128)

In [3]: random = np.random.RandomState(seed=15206)

让我们创建一个真正的输入数组，有 6 行和 6 列:

In [4]: x = random.randn(6, 6)

In [5]: x
Out[5]: 
array([[ 1.577,  0.426,  0.322, -0.891, -0.793,  0.017],
       [ 0.238,  0.603, -0.094, -0.087, -0.936, -1.139],
       [-0.583,  0.394,  0.323, -1.384,  1.255,  0.457],
       [-0.186,  0.687, -0.815, -0.54 ,  0.762, -0.674],
       [-1.604, -0.557,  1.933, -1.122, -0.516, -1.51 ],
       [-1.683, -0.006, -1.648, -0.016,  1.145,  0.809]])

现在看一下完整的 FFT(使用 fft2，而不是 rfft2):

In [6]: fft2_result = np.fft.fft2(x)

In [7]: fft2_result
Out[7]: 
array([[ -5.834+0.j   ,   1.084-2.33j ,  -6.504-3.884j,   3.228-0.j   ,  -6.504+3.884j,   1.084+2.33j ],
       [  1.475-3.311j,   1.865-3.699j,   2.777-0.095j,  -2.570-1.152j,   4.705-3.373j,   4.555-3.657j],
       [  2.758+3.339j,  -3.512+0.398j,   5.824-4.045j,   1.149-3.705j,   0.661-2.127j,  12.368+1.464j],
       [  1.326-0.j   ,   1.191-4.479j,  -3.263+6.19j ,   8.939-0.j   ,  -3.263-6.19j ,   1.191+4.479j],
       [  2.758-3.339j,  12.368-1.464j,   0.661+2.127j,   1.149+3.705j,   5.824+4.045j,  -3.512-0.398j],
       [  1.475+3.311j,   4.555+3.657j,   4.705+3.373j,  -2.570+1.152j,   2.777+0.095j,   1.865+3.699j]])

注意这里有一个对称性:对于任何索引 i和 j与 0 <= i < 6和 0 <= j < 6 , fft2_result[i, j]是 fft_result[-i, -j] 的复共轭.例如:

In [8]: fft2_result[2, 4]
Out[8]: (0.66075993512998199-2.127249005984857j)

In [9]: fft2_result[-2, -4].conj()
Out[9]: (0.66075993512998199-2.127249005984857j)

这意味着我们不需要包含输出的右半部分，因为它可以从左半部分导出。我们可以通过只计算完整 FFT 的左半部分来节省内存，也许还可以节省一点时间。这正是rfft2做:

In [10]: rfft2_result = np.fft.rfft2(x)

In [11]: rfft2_result
Out[11]: 
array([[ -5.834+0.j   ,   1.084-2.33j ,  -6.504-3.884j,   3.228+0.j   ],
       [  1.475-3.311j,   1.865-3.699j,   2.777-0.095j,  -2.570-1.152j],
       [  2.758+3.339j,  -3.512+0.398j,   5.824-4.045j,   1.149-3.705j],
       [  1.326-0.j   ,   1.191-4.479j,  -3.263+6.19j ,   8.939-0.j   ],
       [  2.758-3.339j,  12.368-1.464j,   0.661+2.127j,   1.149+3.705j],
       [  1.475+3.311j,   4.555+3.657j,   4.705+3.373j,  -2.570+1.152j]])

注意 rfft2_result火柴fft2_result[:, :4] ，至少达到数值错误:

In [12]: np.allclose(rfft2_result, fft2_result[:, :4])
Out[12]: True

我们还可以选择保留输出的上半部分而不是左半部分，方法是使用 axes np.fft.rfft2 的参数:

In [13]: np.fft.rfft2(x, axes=[1, 0])
Out[13]: 
array([[ -5.834+0.j   ,   1.084-2.33j ,  -6.504-3.884j,   3.228-0.j   ,  -6.504+3.884j,   1.084+2.33j ],
       [  1.475-3.311j,   1.865-3.699j,   2.777-0.095j,  -2.570-1.152j,   4.705-3.373j,   4.555-3.657j],
       [  2.758+3.339j,  -3.512+0.398j,   5.824-4.045j,   1.149-3.705j,   0.661-2.127j,  12.368+1.464j],
       [  1.326+0.j   ,   1.191-4.479j,  -3.263+6.19j ,   8.939-0.j   ,  -3.263-6.19j ,   1.191+4.479j]])

作为documentation对于 np.fft.rfftn表示，NumPy 在指定的最后一个轴上执行实数 FFT，在其他轴上执行复数 FFT。

当然，rfft2_result 还是有些冗余的:我们可以丢弃第一列的下半部分和最后一列的下半部分，并且仍然能够使用与以前相同的对称性来重建它们。以及位置 [0, 0] 的条目, [0, 3] , [3, 0]和 [3, 3]都是真实的，所以我们可以丢弃它们的虚部。但这会给我们留下不太方便的数组表示。