我应该如何解释 numpy.fft.rfft2 的输出？
How should I interpret the output of numpy.fft.rfft2?
显然 rfft2 函数只是计算输入矩阵的离散 fft。但是，如何解释给定的输出索引？给定输出索引，我在看哪个傅里叶系数？
我对输出的大小特别困惑。对于 n x n 矩阵，输出似乎是 n x (n/2)+1 矩阵（对于偶数 n）。为什么方阵最后会变成非方阵傅立叶变换？
来自文档：
https://docs.scipy.org/doc/numpy/reference/generated/numpy.fft.rfft2.html
This is really just rfftn with different default behavior. For more
details see rfftn.
numpy.fft.rfft2(a, s=None, axes=(-2, -1), norm=None)[source]
对比
numpy.fft.rfftn(a, s=None, axes=None, norm=None)[source]
https://docs.scipy.org/doc/numpy/reference/generated/numpy.fft.rfftn.html#numpy.fft.rfftn
Notes
The transform for real input is performed over the last transformation
axis, as by rfft, then the transform over the remaining axes is
performed as by fftn. The order of the output is as for rfft for the
final transformation axis, and as for fftn for the remaining
transformation axes.
See fft for details, definitions and conventions used.
https://docs.scipy.org/doc/numpy/reference/generated/numpy.fft.fft.html#numpy.fft.fft
我没有使用过 2d fftn，但我创建这个是为了解释 1d fft，这可能会让您对 2d 输出的解释有一些了解：
import math
import numpy as np
PERIOD = 30
SIFT = 2 # integer from 1 to PERIOD/2
def fourier_series(array, period, sift):
# Create an array of length data period; then reverse its order
signal = (array[-period:])[::-1]
# Extract amplitude and phase in (a + bi) complex form
complex_fft = np.fft.fft(signal)
''' Calculate amplitude, phase, frequency, and velocity '''
# Define empty lists for later use
amplitude = []
phase = []
frequency = []
velocity = []
# Extract real and imaginary coefficients from complex scipy output
for n in range(period, 0, -1):
amplitude.append(complex_fft.real[-n])
phase.append(complex_fft.imag[-n])
# The final equation will need to be divided by period
# I do it here so that it is calculated once saving cycles
amplitude = [(x/period) for x in amplitude]
# Extract the carrier
carrier = max(amplitude)
# The frequency is a helper function of fft
# It only has access to the length of the data set
frequency.append(np.fft.fftfreq(signal.size, 1))
# Convert frequency array to list
frequency = frequency[-1]
# Velocity is 2*pi*frequency; I do this here once to save cpu time
velocity = [x*2*math.pi for x in frequency]
''' Calculate the Full Spectrum Sinusoid '''
# Recombine ALL elements in the form An*sin(2*pi(Fn) + Pn) for full spectrum
full_spectrum = 0
for m in range(1, period+1):
full_spectrum += amplitude[-m]*(1+math.sin(velocity[-m] + phase[-m]))
''' Calculate the Filtered Sinusoid '''
# Normalize user sift input as an integer
sift = int(sift)
# If sift is more than half of the period, return full spectrum
if sift >= period/2:
filtered_transform = full_spectrum
# If sift is 0 or 1, return the carrier
else:
filtered_transform = carrier
# For every whole number of sift over 1, but less than 0.5*period:
# Add an 2 elements to the sinusoid respective of
# a negative and positive frequency pair
if sift > 1:
for m in range(1, sift):
p = period - m
filtered_transform += amplitude[-m]*(1+math.sin(velocity[-m] + phase[-m]))
filtered_transform += amplitude[-p]*(1+math.sin(velocity[-p] + phase[-p]))
''' Print Elements and Return FFT'''
if 1:
print('**********************************')
print('Carrier: %.3f' % amplitude[-period])
print(['%.2f' % x for x in amplitude])
print(['%.2f' % x for x in velocity])
print(['%.2f' % x for x in phase])
return filtered_transform, carrier, full_spectrum
stochastic = # Your 1d input array
y, y_carrier, y_full = fourier_series(stochastic, PERIOD, SIFT)
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 的输出和完整 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
我们也可以选择保留输出的上半部分而不是左半部分，方法是使用 np.fft.rfft2 的 axes 参数：
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 for np.fft.rfftn 所说，NumPy 在指定的最后一个轴上执行实数 FFT，在其他轴上执行复数 FFT。
当然，rfft2_result 中仍有一些冗余：我们可以丢弃第一列的下半部分和最后一列的下半部分，并且仍然能够使用相同的对称性重建它们像之前一样。位置 [0, 0]、[0, 3]、[3, 0] 和 [3, 3] 的条目都将是实数，因此我们可以丢弃它们的虚部。但这会给我们留下不太方便的数组表示形式。
还要注意 fft 输出中系数的顺序：
根据文档：默认情况下，第一个元素是 0 频率分量的系数（实际上是数组的总和或平均值），从第二个元素开始，我们有正频率的递增顺序的系数，并且从 n/2+1 开始，它们是按降序排列的负频率。要查看长度为 10 的数组的频率：
np.fft.fftfreq(10)
输出是：
array([ 0. , 0.1, 0.2, 0.3, 0.4, -0.5, -0.4, -0.3, -0.2, -0.1])
使用 np.fft.fftshift(cf)，其中 cf=np.fft.fft(array)，输出被移动，使其对应于此频率排序：
array([-0.5, -0.4, -0.3, -0.2, -0.1, 0. , 0.1, 0.2, 0.3, 0.4])
这对绘图有意义。
二维情况下也是一样。 fft2 和 rfft2 的区别正如其他人所解释的那样。
显然 rfft2 函数只是计算输入矩阵的离散 fft。但是，如何解释给定的输出索引？给定输出索引，我在看哪个傅里叶系数？
我对输出的大小特别困惑。对于 n x n 矩阵，输出似乎是 n x (n/2)+1 矩阵（对于偶数 n）。为什么方阵最后会变成非方阵傅立叶变换？
来自文档：
https://docs.scipy.org/doc/numpy/reference/generated/numpy.fft.rfft2.html
This is really just rfftn with different default behavior. For more details see rfftn.
numpy.fft.rfft2(a, s=None, axes=(-2, -1), norm=None)[source]
对比
numpy.fft.rfftn(a, s=None, axes=None, norm=None)[source]
https://docs.scipy.org/doc/numpy/reference/generated/numpy.fft.rfftn.html#numpy.fft.rfftn
Notes
The transform for real input is performed over the last transformation axis, as by rfft, then the transform over the remaining axes is performed as by fftn. The order of the output is as for rfft for the final transformation axis, and as for fftn for the remaining transformation axes.
See fft for details, definitions and conventions used.
https://docs.scipy.org/doc/numpy/reference/generated/numpy.fft.fft.html#numpy.fft.fft
我没有使用过 2d fftn，但我创建这个是为了解释 1d fft，这可能会让您对 2d 输出的解释有一些了解：
import math
import numpy as np
PERIOD = 30
SIFT = 2 # integer from 1 to PERIOD/2
def fourier_series(array, period, sift):
# Create an array of length data period; then reverse its order
signal = (array[-period:])[::-1]
# Extract amplitude and phase in (a + bi) complex form
complex_fft = np.fft.fft(signal)
''' Calculate amplitude, phase, frequency, and velocity '''
# Define empty lists for later use
amplitude = []
phase = []
frequency = []
velocity = []
# Extract real and imaginary coefficients from complex scipy output
for n in range(period, 0, -1):
amplitude.append(complex_fft.real[-n])
phase.append(complex_fft.imag[-n])
# The final equation will need to be divided by period
# I do it here so that it is calculated once saving cycles
amplitude = [(x/period) for x in amplitude]
# Extract the carrier
carrier = max(amplitude)
# The frequency is a helper function of fft
# It only has access to the length of the data set
frequency.append(np.fft.fftfreq(signal.size, 1))
# Convert frequency array to list
frequency = frequency[-1]
# Velocity is 2*pi*frequency; I do this here once to save cpu time
velocity = [x*2*math.pi for x in frequency]
''' Calculate the Full Spectrum Sinusoid '''
# Recombine ALL elements in the form An*sin(2*pi(Fn) + Pn) for full spectrum
full_spectrum = 0
for m in range(1, period+1):
full_spectrum += amplitude[-m]*(1+math.sin(velocity[-m] + phase[-m]))
''' Calculate the Filtered Sinusoid '''
# Normalize user sift input as an integer
sift = int(sift)
# If sift is more than half of the period, return full spectrum
if sift >= period/2:
filtered_transform = full_spectrum
# If sift is 0 or 1, return the carrier
else:
filtered_transform = carrier
# For every whole number of sift over 1, but less than 0.5*period:
# Add an 2 elements to the sinusoid respective of
# a negative and positive frequency pair
if sift > 1:
for m in range(1, sift):
p = period - m
filtered_transform += amplitude[-m]*(1+math.sin(velocity[-m] + phase[-m]))
filtered_transform += amplitude[-p]*(1+math.sin(velocity[-p] + phase[-p]))
''' Print Elements and Return FFT'''
if 1:
print('**********************************')
print('Carrier: %.3f' % amplitude[-period])
print(['%.2f' % x for x in amplitude])
print(['%.2f' % x for x in velocity])
print(['%.2f' % x for x in phase])
return filtered_transform, carrier, full_spectrum
stochastic = # Your 1d input array
y, y_carrier, y_full = fourier_series(stochastic, PERIOD, SIFT)
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 的输出和完整 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
我们也可以选择保留输出的上半部分而不是左半部分，方法是使用 np.fft.rfft2 的 axes 参数：
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 for np.fft.rfftn 所说，NumPy 在指定的最后一个轴上执行实数 FFT，在其他轴上执行复数 FFT。
当然，rfft2_result 中仍有一些冗余：我们可以丢弃第一列的下半部分和最后一列的下半部分，并且仍然能够使用相同的对称性重建它们像之前一样。位置 [0, 0]、[0, 3]、[3, 0] 和 [3, 3] 的条目都将是实数，因此我们可以丢弃它们的虚部。但这会给我们留下不太方便的数组表示形式。
还要注意 fft 输出中系数的顺序：
根据文档：默认情况下，第一个元素是 0 频率分量的系数（实际上是数组的总和或平均值），从第二个元素开始，我们有正频率的递增顺序的系数，并且从 n/2+1 开始，它们是按降序排列的负频率。要查看长度为 10 的数组的频率：
np.fft.fftfreq(10)
输出是：
array([ 0. , 0.1, 0.2, 0.3, 0.4, -0.5, -0.4, -0.3, -0.2, -0.1])
使用 np.fft.fftshift(cf)，其中 cf=np.fft.fft(array)，输出被移动，使其对应于此频率排序：
array([-0.5, -0.4, -0.3, -0.2, -0.1, 0. , 0.1, 0.2, 0.3, 0.4])
这对绘图有意义。
二维情况下也是一样。 fft2 和 rfft2 的区别正如其他人所解释的那样。