我应该如何解释numpy.fft.rfft2的输出? [英] How should I interpret the output of numpy.fft.rfft2?

问题描述

显然，rfft2函数只是计算输入矩阵的离散fft.但是，如何解释输出的给定索引?给定输出的索引，我在看哪个傅里叶系数?
输出的大小让我特别困惑.对于n×n矩阵，输出似乎是n×(n/2)+1矩阵(对于偶数n).为什么方形矩阵以非方形傅立叶变换结束?

Obviously the rfft2 function simply computes the discrete fft of the input matrix. However how do I interpret a given index of the output? Given an index of the output, which Fourier coefficient am I looking at?
I am especially confused by the sizes of the output. For an n by n matrix, the output seems to be an n by (n/2)+1 matrix (for even n). Why does a square matrix ends up with a non-square fourier transform?

推荐答案

The output of numpy.fft.rfft2 is simply the left half (plus one column) of a standard two-dimensional FFT, as computed by numpy.fft.fft2. There's no need for rfft2 to supply the right half of the result, because the FFT of a real array has a natural and simple symmetry, and the right half of the full FFT can therefore be derived from the left half using that symmetry.

这里是一个例子，以进行说明.首先，为了易于复制和查看，我将设置NumPy的随机状态和打印选项:

Here's an example, to illustrate. First, to make it easy to reproduce and easy to see, I'll set NumPy's random state and printing options:

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列的实际输入数组:

Let's create a real input array, with 6 rows and 6 columns:

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):

Now take a look at the full FFT (using fft2, not 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]])

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

Notice that there's a symmetry here: for any indices i and j with 0 <= i < 6 and 0 <= j < 6, fft2_result[i, j] is the complex conjugate of fft_result[-i, -j]. For example:

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的作用:

That means that we don't need to include the right-hand half of the output, since it can be derived from the left half. We can save memory and perhaps also a tiny bit of time by only computing the left half of the full FFT. And that's exactly what rfft2 does:

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]相匹配，至少不超过数字错误:

Notice that rfft2_result matches fft2_result[:, :4], at least up to numerical error:

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

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

We could also choose to keep the top half of the output rather than the left half, by using the axes argument to 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]])

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

As the documentation for np.fft.rfftn says, NumPy performs a real FFT over the last axis specified, and complex FFTs over the other axes.

当然，在rfft2_result中仍然存在一些冗余:我们可以丢弃第一列的下半部分和最后一列的下半部分，并且仍然能够使用与以前相同的对称性来重构它们.并且位置[0, 0]，[0, 3]，[3, 0]和[3, 3]的条目都是真实的，因此我们可以丢弃它们的虚部.但这会使我们不方便使用数组表示形式.

Of course, there's still some redundancy in rfft2_result: we could discard the bottom half of the first column, and the bottom half of the last column, and still be able to reconstruct them using the same symmetry as before. And the entries at positions [0, 0], [0, 3], [3, 0] and [3, 3] are all going to be real, so we could discard their imaginary parts. But that would leave us with a less convenient array representation.

这篇关于我应该如何解释numpy.fft.rfft2的输出?的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！