如何将两个 2D RFFT 数组（FFTPACK）相乘以与 NumPy 的 FFT 兼容？
How to multiply two 2D RFFT arrays (FFTPACK) to be compatible with NumPy's FFT?
我正在尝试将用 fftpack_rfft2d()（SciPy 的 FFTPACK RFFT）转换的两个二维数组相乘，结果与我从 scipy_rfft2d() 得到的结果不兼容（SciPy 的 FFT RFFT）。
下图分享了脚本的输出，其中显示：
- 两个输入数组的初始化值；
- 两个数组在使用 SciPy 对 RFFT 的 FFT 实现进行转换后使用
scipy_rfft2d()，然后是乘法输出在使用 scipy_irfft2d() 向后转换后输出；
- 使用 SciPy 的 FFTPACK 实现与
fftpack_rfft2d() 和 fftpack_irfft2d() 的 RFFT 相同；
np.allclose() 的测试结果，检查两个乘法的结果在用各自的 IRFFT 实现转换回来后是否相同。
为了清楚起见，红色矩形显示逆变换IRFFT后的乘法结果：左边的矩形使用SciPy的FFT IRFFT；右边的矩形，SciPy的FFTPACK IRFFT。当与 FFTPACK 版本的乘法固定时，它们应该呈现相同的数据。
我认为 FFTPACK 版本的乘法结果不正确，因为 scipy.fftpack returns 结果 RFFT 数组中的实部和虚部不同于来自 scipy.fft:
的 RFFT
- 我相信来自 scipy.fftpack returns 的 RFFT 是一个数组，其中一个元素包含实部，下一个元素包含其虚部；
- 在scipy.fft的RFFT中，每个元素都是一个复数，因此能够同时包含实部和虚部；
如有错误请指正！我还想指出，由于 scipy.fftpack 不提供转换二维数组的函数，如 rfft2() 和 irfft2()，我提供我自己在下面代码中的实现：
import numpy as np
from scipy import fftpack as scipy_fftpack
from scipy import fft as scipy_fft
# SCIPY RFFT 2D
def scipy_rfft2d(matrix):
fftRows = [scipy_fft.rfft(row) for row in matrix]
return np.transpose([scipy_fft.fft(row) for row in np.transpose(fftRows)])
# SCIPY IRFFT 2D
def scipy_irfft2d(matrix, s):
fftRows = [scipy_fft.irfft(row) for row in matrix]
return np.transpose([scipy_fft.ifft(row) for row in np.transpose(fftRows)])
# FFTPACK RFFT 2D
def fftpack_rfft2d(matrix):
fftRows = [scipy_fftpack.rfft(row) for row in matrix]
return np.transpose([scipy_fftpack.rfft(row) for row in np.transpose(fftRows)])
# FFTPACK IRFFT 2D
def fftpack_irfft2d(matrix):
fftRows = [scipy_fftpack.irfft(row) for row in matrix]
return np.transpose([scipy_fftpack.irfft(row) for row in np.transpose(fftRows)])
print('\n#################### INPUT DATA ###################\n')
# initialize two 2D arrays with random data for testing
in1 = np.array([[0, 0, 0, 0], \
[0, 255, 255, 0], \
[0, 0, 255, 255], \
[0, 0, 0, 0]])
print('\nin1 shape=', in1.shape, '\n', in1)
in2 = np.array([[0, 0, 0, 0], \
[0, 0, 255, 0], \
[0, 255, 255, 0], \
[0, 255, 0, 0]])
print('\nin2 shape=', in2.shape, '\n', in2)
print('\n############### SCIPY: 2D RFFT (MULT) ###############\n')
# transform both inputs with SciPy RFFT for 2D
scipy_rfft1 = scipy_fft.rfftn(in1)
scipy_rfft2 = scipy_fft.rfftn(in2)
print('* Output from scipy_fft.rfftn():')
print('scipy_fft1 shape=', scipy_rfft1.shape, '\n', scipy_rfft1.real)
print('\nscipy_fft2 shape=', scipy_rfft2.shape, '\n', scipy_rfft2.real)
# perform multiplication between two 2D arrays from SciPy RFFT
scipy_rfft_mult = scipy_rfft1 * scipy_rfft2
# perform inverse RFFT for 2D arrays using SciPy
scipy_data = scipy_fft.irfftn(scipy_rfft_mult, in1.shape) # passing shape guarantees the output will have the original data size
print('\n* Output from scipy_fft.irfftn():')
print('scipy_data shape=', scipy_data.shape, '\n', scipy_data)
print('\n############### FFTPACK: 2D RFFT (MULT) ###############\n')
# transform both inputs with FFTPACK RFFT for 2D
fftpack_rfft1 = fftpack_rfft2d(in1)
fftpack_rfft2 = fftpack_rfft2d(in2)
print('* Output from fftpack_rfft2d():')
print('fftpack_rfft1 shape=', fftpack_rfft1.shape, '\n', fftpack_rfft1)
print('\nfftpack_rfft2 shape=', fftpack_rfft2.shape, '\n', fftpack_rfft2)
# TODO: perform multiplication between two 2D arrays from FFTPACK RFFT
fftpack_rfft_mult = fftpack_rfft1 * fftpack_rfft2 # this doesn't work
# perform inverse RFFT for 2D arrays using FFTPACK
fftpack_data = fftpack_irfft2d(fftpack_rfft_mult)
print('\n* Output from fftpack_irfft2d():')
print('fftpack_data shape=', fftpack_data.shape, '\n', fftpack_data)
print('\n##################### RESULT #####################\n')
# compare FFTPACK result with SCIPY
print('\nIs fftpack_data equivalent to scipy_data?', np.allclose(fftpack_data, scipy_data), '\n')
假设我的猜测是正确的，对于乘以从 fftpack_rfft2d() 生成的两个二维数组的函数，正确的实现是什么？请记住：生成的数组必须能够用 fftpack_irfft2d().
转换回来
只邀请解决二维问题的答案。那些对如何乘以 1D FFTPACK 数组感兴趣的人可以 check this thread.
你的假设是正确的。 FFTPACK returns 格式为
的单个实向量中的所有系数
[y(0),Re(y(1)),Im(y(1)),...,Re(y(n/2))] if n is even
[y(0),Re(y(1)),Im(y(1)),...,Re(y(n/2)),Im(y(n/2))] if n is odd
其中 scipy.rfft returns 一个复数向量
[y(0),Re(y(1)) + 1.0j*Im(y(1)),...,Re(y(n/2) + 1.0j*Im(y(n/2)))]
所以你需要使用适当的步长来形成一个向量，如下：
y_fft = np.cat([y_fftpack[0], y_fftpack[1:2:] + 1.0j*y_fftpack[2:2:]])
要将 2 个复数数组相乘，您必须进行复数乘法。
参见 https://en.m.wikipedia.org/wiki/Complex_number
运算部分的乘法
您不能只乘以实部，然后分别乘以虚部或拆分元素，这可能是您的 fftpack 矩阵 mul 产生垃圾的原因。
正确的函数：
import numpy as np
from scipy import fftpack as scipy_fftpack
from scipy import fft as scipy
# FFTPACK RFFT 2D
def fftpack_rfft2d(matrix):
fftRows = scipy_fftpack.fft(matrix, axis=1)
fftCols = scipy_fftpack.fft(fftRows, axis=0)
return fftCols
# FFTPACK IRFFT 2D
def fftpack_irfft2d(matrix):
ifftRows = scipy_fftpack.ifft(matrix, axis=1)
ifftCols = scipy_fftpack.ifft(ifftRows, axis=0)
return ifftCols.real
您计算 2D FFT 的方法有误。是的，可以使用 rfft() 计算第一个 FFT（按您的情况的列），但第二个 FFT 计算 必须 复杂 第一个 FFT（按列） 的输出，因此 rfft() 的输出必须转换成真实的复谱。此外，这意味着您 必须 使用 fft() 而不是 rfft()按行进行 FFT。因此，在两个计算中使用fft()更方便。
而且，你的输入数据是numpy二维数组，为什么要使用列表理解？直接用fftpack.fft()，这样快很多。
- 如果您已经只有由错误函数计算出的二维数组并且需要将它们相乘： 那么，我认为，尝试使用相同的方法从错误的二维 FFT 重建输入数据'wrong' 方式然后计算 正确的 2D FFT
======================================== ========================
新功能完整测试代码版本：
import numpy as np
from scipy import fftpack as scipy_fftpack
from scipy import fft as scipy_fft
# FFTPACK RFFT 2D
def fftpack_rfft2d(matrix):
fftRows = scipy_fftpack.fft(matrix, axis=1)
fftCols = scipy_fftpack.fft(fftRows, axis=0)
return fftCols
# FFTPACK IRFFT 2D
def fftpack_irfft2d(matrix):
ifftRows = scipy_fftpack.ifft(matrix, axis=1)
ifftCols = scipy_fftpack.ifft(ifftRows, axis=0)
return ifftCols.real
print('\n#################### INPUT DATA ###################\n')
# initialize two 2D arrays with random data for testing
in1 = np.array([[0, 0, 0, 0], \
[0, 255, 255, 0], \
[0, 0, 255, 255], \
[0, 0, 0, 0]])
print('\nin1 shape=', in1.shape, '\n', in1)
in2 = np.array([[0, 0, 0, 0], \
[0, 0, 255, 0], \
[0, 255, 255, 0], \
[0, 255, 0, 0]])
print('\nin2 shape=', in2.shape, '\n', in2)
print('\n############### SCIPY: 2D RFFT (MULT) ###############\n')
# transform both inputs with SciPy RFFT for 2D
scipy_rfft1 = scipy_fft.fftn(in1)
scipy_rfft2 = scipy_fft.fftn(in2)
print('* Output from scipy_fft.rfftn():')
print('scipy_fft1 shape=', scipy_rfft1.shape, '\n', scipy_rfft1)
print('\nscipy_fft2 shape=', scipy_rfft2.shape, '\n', scipy_rfft2)
# perform multiplication between two 2D arrays from SciPy RFFT
scipy_rfft_mult = scipy_rfft1 * scipy_rfft2
# perform inverse RFFT for 2D arrays using SciPy
scipy_data = scipy_fft.irfftn(scipy_rfft_mult, in1.shape) # passing shape guarantees the output will
# have the original data size
print('\n* Output from scipy_fft.irfftn():')
print('scipy_data shape=', scipy_data.shape, '\n', scipy_data)
print('\n############### FFTPACK: 2D RFFT (MULT) ###############\n')
# transform both inputs with FFTPACK RFFT for 2D
fftpack_rfft1 = fftpack_rfft2d(in1)
fftpack_rfft2 = fftpack_rfft2d(in2)
print('* Output from fftpack_rfft2d():')
print('fftpack_rfft1 shape=', fftpack_rfft1.shape, '\n', fftpack_rfft1)
print('\nfftpack_rfft2 shape=', fftpack_rfft2.shape, '\n', fftpack_rfft2)
# TODO: perform multiplication between two 2D arrays from FFTPACK RFFT
fftpack_rfft_mult = fftpack_rfft1 * fftpack_rfft2 # this doesn't work
# perform inverse RFFT for 2D arrays using FFTPACK
fftpack_data = fftpack_irfft2d(fftpack_rfft_mult)
print('\n* Output from fftpack_irfft2d():')
print('fftpack_data shape=', fftpack_data.shape, '\n', fftpack_data)
print('\n##################### RESULT #####################\n')
# compare FFTPACK result with SCIPY
print('\nIs fftpack_data equivalent to scipy_data?', np.allclose(fftpack_data, scipy_data), '\n')
输出为：
#################### INPUT DATA ###################
in1 shape= (4, 4)
[[ 0 0 0 0]
[ 0 255 255 0]
[ 0 0 255 255]
[ 0 0 0 0]]
in2 shape= (4, 4)
[[ 0 0 0 0]
[ 0 0 255 0]
[ 0 255 255 0]
[ 0 255 0 0]]
############### SCIPY: 2D RFFT (MULT) ###############
* Output from scipy_fft.rfftn():
scipy_fft1 shape= (4, 4)
[[1020. -0.j -510. +0.j 0. -0.j -510. -0.j]
[-510.-510.j 0. +0.j 0. +0.j 510.+510.j]
[ 0. -0.j 0.+510.j 0. -0.j 0.-510.j]
[-510.+510.j 510.-510.j 0. -0.j 0. -0.j]]
scipy_fft2 shape= (4, 4)
[[1020. -0.j -510.-510.j 0. -0.j -510.+510.j]
[-510. +0.j 510.+510.j 0.-510.j 0. -0.j]
[ 0. -0.j 0. +0.j 0. -0.j 0. -0.j]
[-510. -0.j 0. +0.j 0.+510.j 510.-510.j]]
* Output from scipy_fft.irfftn():
scipy_data shape= (4, 4)
[[130050. 65025. 65025. 130050.]
[ 65025. 0. 0. 65025.]
[ 65025. 0. 0. 65025.]
[130050. 65025. 65025. 130050.]]
############### FFTPACK: 2D RFFT (MULT) ###############
* Output from fftpack_rfft2d():
fftpack_rfft1 shape= (4, 4)
[[1020. -0.j -510. +0.j 0. -0.j -510. +0.j]
[-510.-510.j 0. +0.j 0. +0.j 510.+510.j]
[ 0. +0.j 0.+510.j 0. +0.j 0.-510.j]
[-510.+510.j 510.-510.j 0. +0.j 0. +0.j]]
fftpack_rfft2 shape= (4, 4)
[[1020. -0.j -510.-510.j 0. -0.j -510.+510.j]
[-510. +0.j 510.+510.j 0.-510.j 0. +0.j]
[ 0. +0.j 0. +0.j 0. +0.j 0. +0.j]
[-510. +0.j 0. +0.j 0.+510.j 510.-510.j]]
* Output from fftpack_irfft2d():
fftpack_data shape= (4, 4)
[[130050.+0.j 65025.+0.j 65025.+0.j 130050.+0.j]
[ 65025.+0.j 0.+0.j 0.+0.j 65025.+0.j]
[ 65025.+0.j 0.+0.j 0.+0.j 65025.+0.j]
[130050.+0.j 65025.+0.j 65025.-0.j 130050.+0.j]]
##################### RESULT #####################
Is fftpack_data equivalent to scipy_data? True
是对的：只用复数FFT就简单多了（虽然他的实现复杂的没必要，直接用scipy.fftpack.fft2）。正如我在评论中所说，最好的选择是切换到 scipy.fft，这样使用起来更简单； fftpack 已弃用，取而代之的是它。
但是，如果您确实需要使用 fftpack，并且您确实希望通过使用 rfft 函数来节省一些计算时间，那么这是正确的方法。它需要在沿另一个维度计算 fft 之前，将 rfft 函数的实数值输出转换为复数值数组。使用此解决方案，下面的 fftpack_rfft2d 输出其输入的一半 2D FFT，另一半是冗余的。
import numpy as np
from scipy import fftpack
# FFTPACK RFFT 2D
def fftpack_rfft1d(matrix):
assert not (matrix.shape[1] & 0x1)
tmp = fftpack.rfft(matrix, axis=1)
assert tmp.dtype == np.dtype('float64')
return np.hstack((tmp[:, [0]], np.ascontiguousarray(tmp[:, 1:-1]).view(np.complex128), tmp[:, [-1]]))
def fftpack_rfft2d(matrix):
return fftpack.fft(fftpack_rfft1d(matrix), axis=0)
# FFTPACK IRFFT 2D
def fftpack_irfft1d(matrix):
assert matrix.dtype == np.dtype('complex128')
tmp = np.hstack((matrix[:, [0]].real, np.ascontiguousarray(matrix[:, 1:-1]).view(np.float64), matrix[:, [-1]].real))
return fftpack.irfft(tmp, axis=1)
def fftpack_irfft2d(matrix):
return fftpack_irfft1d(fftpack.ifft(matrix, axis=0))
######
# test data
in1 = np.random.randn(256,256)
in2 = np.random.randn(256,256)
# fftpack.fft2
gt_result = fftpack.ifft2(fftpack.fft2(in1) * fftpack.fft2(in2)).real
# fftpack_rfft2d
our_result = fftpack_irfft2d(fftpack_rfft2d(in1) * fftpack_rfft2d(in2) )
# compare
print('\nIs our result equivalent to the ground truth?', np.allclose(gt_result, our_result), '\n')
[此代码仅适用于偶数大小的图像，我没有费心让它通用，请参阅 here 了解如何操作）。
尽管如此，由于此解决方案需要数据副本，因此它实际上比仅使用普通的复值 FFT (fftpack.fft2) 更慢，尽管它进行的计算更少：
import time
tic = time.perf_counter()
for i in range(100):
fftpack.fft(in1)
toc = time.perf_counter()
print(f"fftpack.fft() takes {toc - tic:0.4f} seconds")
tic = time.perf_counter()
for i in range(100):
fftpack_rfft2d(in1)
toc = time.perf_counter()
print(f"fftpack_rfft2d() takes {toc - tic:0.4f} seconds")
输出：
fftpack.fft() takes 0.0442 seconds
fftpack_rfft2d() takes 0.0664 seconds
所以，确实，坚持fftpack.fft（或者如果可以的话，更确切地说scipy.fft.fft）。
除了@CrisLuengo 的回答（）。
性能测试
测试 fftpack.FFT 与 fftpack.RFFT - 1D
# test data
sz =50000
sz = fftpack.next_fast_len(sz)
in1 = np.random.randn(sz)
print(f"Input (len = {len(in1)}):", sep='\n')
rep = 1000
tic = time.perf_counter()
for i in range(rep):
spec1 = fftpack.fft(in1,axis=0)
toc = time.perf_counter()
print("", f"Spectrum FFT (len = {len(spec1)}):",
f"spec1 takes {10**6*((toc - tic)/rep):0.4f} us", sep="\n")
sz2 = sz//2 + 1
spec2 = np.empty(sz2, dtype=np.complex128)
tic = time.perf_counter()
for i in range(rep):
tmp = fftpack.rfft(in1)
assert tmp.dtype == np.dtype('float64')
if not sz & 0x1:
end = -1
spec2[end] = tmp[end]
else:
end = None
spec2[0] = tmp[0]
spec2[1:end] = tmp[1:end].view(np.complex128)
toc = time.perf_counter()
print("", f"Spectrum RFFT (len = {len(spec2)}):",
f"spec2 takes {10**6*((toc - tic)/rep):0.4f} us", sep="\n")
结果为
Input (len = 50000):
Spectrum FFT (len = 50000):
spec1 takes 583.5880 us
Spectrum RFFT (len = 25001):
spec2 takes 476.0843 us
- So，使用
fftpack.rfft() 并进一步将其输出转换为 complex 视图 比 [ 快 15-20% =18=] 用于大数组.
测试 fftpack.FFT 对比 fftpack.FFT2 - 2D
二维案例的类似测试：
# test data
sz = 5000
in1 = np.random.randn(sz, sz)
print(f"Input (len = {len(in1)}):", sep='\n')
rep = 1
tic = time.perf_counter()
for i in range(rep):
spec1 = np.apply_along_axis(fftpack.fft, 0, in1)
spec1 = np.apply_along_axis(fftpack.fft, 1, spec1)
toc = time.perf_counter()
print("", f"2D Spectrum FFT with np.apply_along_axis (len = {len(spec1)}):",
f"spec1 takes {10**0*((toc - tic)/rep):0.4f} s", sep="\n")
tic = time.perf_counter()
for i in range(rep):
spec2 = fftpack.fft(in1,axis=0)
spec2 = fftpack.fft(spec2,axis=1)
toc = time.perf_counter()
print("", f"2D Spectrum 2xFFT (len = {len(spec2)}):",
f"spec2 takes {10**0*((toc - tic)/rep):0.4f} s", sep="\n")
tic = time.perf_counter()
for i in range(rep):
spec3 = fftpack.fft2(in1)
toc = time.perf_counter()
print("", f"2D Spectrum FFT2 (len = {len(spec3)}):",
f"spec3 takes {10**0*((toc - tic)/rep):0.4f} s", sep="\n")
# compare
print('\nIs spec1 equivalent to the spec2?', np.allclose(spec1, spec2))
print('\nIs spec2 equivalent to the spec3?', np.allclose(spec2, spec3), '\n')
矩阵大小 = 5x5 的结果
Input (len = 5):
2D Spectrum FFT with np.apply_along_axis (len = 5):
spec1 takes 0.000183 s
2D Spectrum 2xFFT (len = 5):
spec2 takes 0.000010 s
2D Spectrum FFT2 (len = 5):
spec3 takes 0.000012 s
Is spec1 equivalent to the spec2? True
Is spec2 equivalent to the spec3? True
矩阵大小 = 500x500 的结果
Input (len = 500):
2D Spectrum FFT with np.apply_along_axis (len = 500):
spec1 takes 0.017626 s
2D Spectrum 2xFFT (len = 500):
spec2 takes 0.005324 s
2D Spectrum FFT2 (len = 500):
spec3 takes 0.003528 s
Is spec1 equivalent to the spec2? True
Is spec2 equivalent to the spec3? True
矩阵大小 = 5000x5000 的结果
Input (len = 5000):
2D Spectrum FFT with np.apply_along_axis (len = 5000):
spec1 takes 2.538471 s
2D Spectrum 2xFFT (len = 5000):
spec2 takes 0.846661 s
2D Spectrum FFT2 (len = 5000):
spec3 takes 0.574397 s
Is spec1 equivalent to the spec2? True
Is spec2 equivalent to the spec3? True
结论
从上面的测试看来，使用 fftpack.fft2() 对于更大的矩阵更有效。
使用np.apply_along_axis()是最慢的方法。
我正在尝试将用 fftpack_rfft2d()（SciPy 的 FFTPACK RFFT）转换的两个二维数组相乘，结果与我从 scipy_rfft2d() 得到的结果不兼容（SciPy 的 FFT RFFT）。
下图分享了脚本的输出，其中显示：
- 两个输入数组的初始化值；
- 两个数组在使用 SciPy 对 RFFT 的 FFT 实现进行转换后使用
scipy_rfft2d()，然后是乘法输出在使用scipy_irfft2d()向后转换后输出； - 使用 SciPy 的 FFTPACK 实现与
fftpack_rfft2d()和fftpack_irfft2d()的 RFFT 相同； np.allclose()的测试结果，检查两个乘法的结果在用各自的 IRFFT 实现转换回来后是否相同。
为了清楚起见，红色矩形显示逆变换IRFFT后的乘法结果：左边的矩形使用SciPy的FFT IRFFT；右边的矩形，SciPy的FFTPACK IRFFT。当与 FFTPACK 版本的乘法固定时，它们应该呈现相同的数据。
我认为 FFTPACK 版本的乘法结果不正确，因为 scipy.fftpack returns 结果 RFFT 数组中的实部和虚部不同于来自 scipy.fft:
的 RFFT- 我相信来自 scipy.fftpack returns 的 RFFT 是一个数组，其中一个元素包含实部，下一个元素包含其虚部；
- 在scipy.fft的RFFT中，每个元素都是一个复数，因此能够同时包含实部和虚部；
如有错误请指正！我还想指出，由于 scipy.fftpack 不提供转换二维数组的函数，如 rfft2() 和 irfft2()，我提供我自己在下面代码中的实现：
import numpy as np
from scipy import fftpack as scipy_fftpack
from scipy import fft as scipy_fft
# SCIPY RFFT 2D
def scipy_rfft2d(matrix):
fftRows = [scipy_fft.rfft(row) for row in matrix]
return np.transpose([scipy_fft.fft(row) for row in np.transpose(fftRows)])
# SCIPY IRFFT 2D
def scipy_irfft2d(matrix, s):
fftRows = [scipy_fft.irfft(row) for row in matrix]
return np.transpose([scipy_fft.ifft(row) for row in np.transpose(fftRows)])
# FFTPACK RFFT 2D
def fftpack_rfft2d(matrix):
fftRows = [scipy_fftpack.rfft(row) for row in matrix]
return np.transpose([scipy_fftpack.rfft(row) for row in np.transpose(fftRows)])
# FFTPACK IRFFT 2D
def fftpack_irfft2d(matrix):
fftRows = [scipy_fftpack.irfft(row) for row in matrix]
return np.transpose([scipy_fftpack.irfft(row) for row in np.transpose(fftRows)])
print('\n#################### INPUT DATA ###################\n')
# initialize two 2D arrays with random data for testing
in1 = np.array([[0, 0, 0, 0], \
[0, 255, 255, 0], \
[0, 0, 255, 255], \
[0, 0, 0, 0]])
print('\nin1 shape=', in1.shape, '\n', in1)
in2 = np.array([[0, 0, 0, 0], \
[0, 0, 255, 0], \
[0, 255, 255, 0], \
[0, 255, 0, 0]])
print('\nin2 shape=', in2.shape, '\n', in2)
print('\n############### SCIPY: 2D RFFT (MULT) ###############\n')
# transform both inputs with SciPy RFFT for 2D
scipy_rfft1 = scipy_fft.rfftn(in1)
scipy_rfft2 = scipy_fft.rfftn(in2)
print('* Output from scipy_fft.rfftn():')
print('scipy_fft1 shape=', scipy_rfft1.shape, '\n', scipy_rfft1.real)
print('\nscipy_fft2 shape=', scipy_rfft2.shape, '\n', scipy_rfft2.real)
# perform multiplication between two 2D arrays from SciPy RFFT
scipy_rfft_mult = scipy_rfft1 * scipy_rfft2
# perform inverse RFFT for 2D arrays using SciPy
scipy_data = scipy_fft.irfftn(scipy_rfft_mult, in1.shape) # passing shape guarantees the output will have the original data size
print('\n* Output from scipy_fft.irfftn():')
print('scipy_data shape=', scipy_data.shape, '\n', scipy_data)
print('\n############### FFTPACK: 2D RFFT (MULT) ###############\n')
# transform both inputs with FFTPACK RFFT for 2D
fftpack_rfft1 = fftpack_rfft2d(in1)
fftpack_rfft2 = fftpack_rfft2d(in2)
print('* Output from fftpack_rfft2d():')
print('fftpack_rfft1 shape=', fftpack_rfft1.shape, '\n', fftpack_rfft1)
print('\nfftpack_rfft2 shape=', fftpack_rfft2.shape, '\n', fftpack_rfft2)
# TODO: perform multiplication between two 2D arrays from FFTPACK RFFT
fftpack_rfft_mult = fftpack_rfft1 * fftpack_rfft2 # this doesn't work
# perform inverse RFFT for 2D arrays using FFTPACK
fftpack_data = fftpack_irfft2d(fftpack_rfft_mult)
print('\n* Output from fftpack_irfft2d():')
print('fftpack_data shape=', fftpack_data.shape, '\n', fftpack_data)
print('\n##################### RESULT #####################\n')
# compare FFTPACK result with SCIPY
print('\nIs fftpack_data equivalent to scipy_data?', np.allclose(fftpack_data, scipy_data), '\n')
假设我的猜测是正确的，对于乘以从 fftpack_rfft2d() 生成的两个二维数组的函数，正确的实现是什么？请记住：生成的数组必须能够用 fftpack_irfft2d().
只邀请解决二维问题的答案。那些对如何乘以 1D FFTPACK 数组感兴趣的人可以 check this thread.
你的假设是正确的。 FFTPACK returns 格式为
的单个实向量中的所有系数[y(0),Re(y(1)),Im(y(1)),...,Re(y(n/2))] if n is even
[y(0),Re(y(1)),Im(y(1)),...,Re(y(n/2)),Im(y(n/2))] if n is odd
其中 scipy.rfft returns 一个复数向量
[y(0),Re(y(1)) + 1.0j*Im(y(1)),...,Re(y(n/2) + 1.0j*Im(y(n/2)))]
所以你需要使用适当的步长来形成一个向量，如下：
y_fft = np.cat([y_fftpack[0], y_fftpack[1:2:] + 1.0j*y_fftpack[2:2:]])
要将 2 个复数数组相乘，您必须进行复数乘法。
参见 https://en.m.wikipedia.org/wiki/Complex_number
运算部分的乘法您不能只乘以实部，然后分别乘以虚部或拆分元素，这可能是您的 fftpack 矩阵 mul 产生垃圾的原因。
正确的函数：
import numpy as np
from scipy import fftpack as scipy_fftpack
from scipy import fft as scipy
# FFTPACK RFFT 2D
def fftpack_rfft2d(matrix):
fftRows = scipy_fftpack.fft(matrix, axis=1)
fftCols = scipy_fftpack.fft(fftRows, axis=0)
return fftCols
# FFTPACK IRFFT 2D
def fftpack_irfft2d(matrix):
ifftRows = scipy_fftpack.ifft(matrix, axis=1)
ifftCols = scipy_fftpack.ifft(ifftRows, axis=0)
return ifftCols.real
您计算 2D FFT 的方法有误。是的，可以使用 rfft() 计算第一个 FFT（按您的情况的列），但第二个 FFT 计算 必须 复杂 第一个 FFT（按列） 的输出，因此 rfft() 的输出必须转换成真实的复谱。此外，这意味着您 必须 使用 fft() 而不是 rfft()按行进行 FFT。因此，在两个计算中使用fft()更方便。
而且，你的输入数据是numpy二维数组，为什么要使用列表理解？直接用fftpack.fft()，这样快很多。
- 如果您已经只有由错误函数计算出的二维数组并且需要将它们相乘： 那么，我认为，尝试使用相同的方法从错误的二维 FFT 重建输入数据'wrong' 方式然后计算 正确的 2D FFT
======================================== ========================
新功能完整测试代码版本：
import numpy as np
from scipy import fftpack as scipy_fftpack
from scipy import fft as scipy_fft
# FFTPACK RFFT 2D
def fftpack_rfft2d(matrix):
fftRows = scipy_fftpack.fft(matrix, axis=1)
fftCols = scipy_fftpack.fft(fftRows, axis=0)
return fftCols
# FFTPACK IRFFT 2D
def fftpack_irfft2d(matrix):
ifftRows = scipy_fftpack.ifft(matrix, axis=1)
ifftCols = scipy_fftpack.ifft(ifftRows, axis=0)
return ifftCols.real
print('\n#################### INPUT DATA ###################\n')
# initialize two 2D arrays with random data for testing
in1 = np.array([[0, 0, 0, 0], \
[0, 255, 255, 0], \
[0, 0, 255, 255], \
[0, 0, 0, 0]])
print('\nin1 shape=', in1.shape, '\n', in1)
in2 = np.array([[0, 0, 0, 0], \
[0, 0, 255, 0], \
[0, 255, 255, 0], \
[0, 255, 0, 0]])
print('\nin2 shape=', in2.shape, '\n', in2)
print('\n############### SCIPY: 2D RFFT (MULT) ###############\n')
# transform both inputs with SciPy RFFT for 2D
scipy_rfft1 = scipy_fft.fftn(in1)
scipy_rfft2 = scipy_fft.fftn(in2)
print('* Output from scipy_fft.rfftn():')
print('scipy_fft1 shape=', scipy_rfft1.shape, '\n', scipy_rfft1)
print('\nscipy_fft2 shape=', scipy_rfft2.shape, '\n', scipy_rfft2)
# perform multiplication between two 2D arrays from SciPy RFFT
scipy_rfft_mult = scipy_rfft1 * scipy_rfft2
# perform inverse RFFT for 2D arrays using SciPy
scipy_data = scipy_fft.irfftn(scipy_rfft_mult, in1.shape) # passing shape guarantees the output will
# have the original data size
print('\n* Output from scipy_fft.irfftn():')
print('scipy_data shape=', scipy_data.shape, '\n', scipy_data)
print('\n############### FFTPACK: 2D RFFT (MULT) ###############\n')
# transform both inputs with FFTPACK RFFT for 2D
fftpack_rfft1 = fftpack_rfft2d(in1)
fftpack_rfft2 = fftpack_rfft2d(in2)
print('* Output from fftpack_rfft2d():')
print('fftpack_rfft1 shape=', fftpack_rfft1.shape, '\n', fftpack_rfft1)
print('\nfftpack_rfft2 shape=', fftpack_rfft2.shape, '\n', fftpack_rfft2)
# TODO: perform multiplication between two 2D arrays from FFTPACK RFFT
fftpack_rfft_mult = fftpack_rfft1 * fftpack_rfft2 # this doesn't work
# perform inverse RFFT for 2D arrays using FFTPACK
fftpack_data = fftpack_irfft2d(fftpack_rfft_mult)
print('\n* Output from fftpack_irfft2d():')
print('fftpack_data shape=', fftpack_data.shape, '\n', fftpack_data)
print('\n##################### RESULT #####################\n')
# compare FFTPACK result with SCIPY
print('\nIs fftpack_data equivalent to scipy_data?', np.allclose(fftpack_data, scipy_data), '\n')
输出为：
#################### INPUT DATA ###################
in1 shape= (4, 4)
[[ 0 0 0 0]
[ 0 255 255 0]
[ 0 0 255 255]
[ 0 0 0 0]]
in2 shape= (4, 4)
[[ 0 0 0 0]
[ 0 0 255 0]
[ 0 255 255 0]
[ 0 255 0 0]]
############### SCIPY: 2D RFFT (MULT) ###############
* Output from scipy_fft.rfftn():
scipy_fft1 shape= (4, 4)
[[1020. -0.j -510. +0.j 0. -0.j -510. -0.j]
[-510.-510.j 0. +0.j 0. +0.j 510.+510.j]
[ 0. -0.j 0.+510.j 0. -0.j 0.-510.j]
[-510.+510.j 510.-510.j 0. -0.j 0. -0.j]]
scipy_fft2 shape= (4, 4)
[[1020. -0.j -510.-510.j 0. -0.j -510.+510.j]
[-510. +0.j 510.+510.j 0.-510.j 0. -0.j]
[ 0. -0.j 0. +0.j 0. -0.j 0. -0.j]
[-510. -0.j 0. +0.j 0.+510.j 510.-510.j]]
* Output from scipy_fft.irfftn():
scipy_data shape= (4, 4)
[[130050. 65025. 65025. 130050.]
[ 65025. 0. 0. 65025.]
[ 65025. 0. 0. 65025.]
[130050. 65025. 65025. 130050.]]
############### FFTPACK: 2D RFFT (MULT) ###############
* Output from fftpack_rfft2d():
fftpack_rfft1 shape= (4, 4)
[[1020. -0.j -510. +0.j 0. -0.j -510. +0.j]
[-510.-510.j 0. +0.j 0. +0.j 510.+510.j]
[ 0. +0.j 0.+510.j 0. +0.j 0.-510.j]
[-510.+510.j 510.-510.j 0. +0.j 0. +0.j]]
fftpack_rfft2 shape= (4, 4)
[[1020. -0.j -510.-510.j 0. -0.j -510.+510.j]
[-510. +0.j 510.+510.j 0.-510.j 0. +0.j]
[ 0. +0.j 0. +0.j 0. +0.j 0. +0.j]
[-510. +0.j 0. +0.j 0.+510.j 510.-510.j]]
* Output from fftpack_irfft2d():
fftpack_data shape= (4, 4)
[[130050.+0.j 65025.+0.j 65025.+0.j 130050.+0.j]
[ 65025.+0.j 0.+0.j 0.+0.j 65025.+0.j]
[ 65025.+0.j 0.+0.j 0.+0.j 65025.+0.j]
[130050.+0.j 65025.+0.j 65025.-0.j 130050.+0.j]]
##################### RESULT #####################
Is fftpack_data equivalent to scipy_data? True
scipy.fftpack.fft2）。正如我在评论中所说，最好的选择是切换到 scipy.fft，这样使用起来更简单； fftpack 已弃用，取而代之的是它。
但是，如果您确实需要使用 fftpack，并且您确实希望通过使用 rfft 函数来节省一些计算时间，那么这是正确的方法。它需要在沿另一个维度计算 fft 之前，将 rfft 函数的实数值输出转换为复数值数组。使用此解决方案，下面的 fftpack_rfft2d 输出其输入的一半 2D FFT，另一半是冗余的。
import numpy as np
from scipy import fftpack
# FFTPACK RFFT 2D
def fftpack_rfft1d(matrix):
assert not (matrix.shape[1] & 0x1)
tmp = fftpack.rfft(matrix, axis=1)
assert tmp.dtype == np.dtype('float64')
return np.hstack((tmp[:, [0]], np.ascontiguousarray(tmp[:, 1:-1]).view(np.complex128), tmp[:, [-1]]))
def fftpack_rfft2d(matrix):
return fftpack.fft(fftpack_rfft1d(matrix), axis=0)
# FFTPACK IRFFT 2D
def fftpack_irfft1d(matrix):
assert matrix.dtype == np.dtype('complex128')
tmp = np.hstack((matrix[:, [0]].real, np.ascontiguousarray(matrix[:, 1:-1]).view(np.float64), matrix[:, [-1]].real))
return fftpack.irfft(tmp, axis=1)
def fftpack_irfft2d(matrix):
return fftpack_irfft1d(fftpack.ifft(matrix, axis=0))
######
# test data
in1 = np.random.randn(256,256)
in2 = np.random.randn(256,256)
# fftpack.fft2
gt_result = fftpack.ifft2(fftpack.fft2(in1) * fftpack.fft2(in2)).real
# fftpack_rfft2d
our_result = fftpack_irfft2d(fftpack_rfft2d(in1) * fftpack_rfft2d(in2) )
# compare
print('\nIs our result equivalent to the ground truth?', np.allclose(gt_result, our_result), '\n')
[此代码仅适用于偶数大小的图像，我没有费心让它通用，请参阅 here 了解如何操作）。
尽管如此，由于此解决方案需要数据副本，因此它实际上比仅使用普通的复值 FFT (fftpack.fft2) 更慢，尽管它进行的计算更少：
import time
tic = time.perf_counter()
for i in range(100):
fftpack.fft(in1)
toc = time.perf_counter()
print(f"fftpack.fft() takes {toc - tic:0.4f} seconds")
tic = time.perf_counter()
for i in range(100):
fftpack_rfft2d(in1)
toc = time.perf_counter()
print(f"fftpack_rfft2d() takes {toc - tic:0.4f} seconds")
输出：
fftpack.fft() takes 0.0442 seconds
fftpack_rfft2d() takes 0.0664 seconds
所以，确实，坚持fftpack.fft（或者如果可以的话，更确切地说scipy.fft.fft）。
除了@CrisLuengo 的回答（
性能测试
测试 fftpack.FFT 与 fftpack.RFFT - 1D
# test data
sz =50000
sz = fftpack.next_fast_len(sz)
in1 = np.random.randn(sz)
print(f"Input (len = {len(in1)}):", sep='\n')
rep = 1000
tic = time.perf_counter()
for i in range(rep):
spec1 = fftpack.fft(in1,axis=0)
toc = time.perf_counter()
print("", f"Spectrum FFT (len = {len(spec1)}):",
f"spec1 takes {10**6*((toc - tic)/rep):0.4f} us", sep="\n")
sz2 = sz//2 + 1
spec2 = np.empty(sz2, dtype=np.complex128)
tic = time.perf_counter()
for i in range(rep):
tmp = fftpack.rfft(in1)
assert tmp.dtype == np.dtype('float64')
if not sz & 0x1:
end = -1
spec2[end] = tmp[end]
else:
end = None
spec2[0] = tmp[0]
spec2[1:end] = tmp[1:end].view(np.complex128)
toc = time.perf_counter()
print("", f"Spectrum RFFT (len = {len(spec2)}):",
f"spec2 takes {10**6*((toc - tic)/rep):0.4f} us", sep="\n")
结果为
Input (len = 50000):
Spectrum FFT (len = 50000):
spec1 takes 583.5880 us
Spectrum RFFT (len = 25001):
spec2 takes 476.0843 us
- So，使用
fftpack.rfft()并进一步将其输出转换为complex视图 比 [ 快 15-20% =18=] 用于大数组.
测试 fftpack.FFT 对比 fftpack.FFT2 - 2D
二维案例的类似测试：
# test data
sz = 5000
in1 = np.random.randn(sz, sz)
print(f"Input (len = {len(in1)}):", sep='\n')
rep = 1
tic = time.perf_counter()
for i in range(rep):
spec1 = np.apply_along_axis(fftpack.fft, 0, in1)
spec1 = np.apply_along_axis(fftpack.fft, 1, spec1)
toc = time.perf_counter()
print("", f"2D Spectrum FFT with np.apply_along_axis (len = {len(spec1)}):",
f"spec1 takes {10**0*((toc - tic)/rep):0.4f} s", sep="\n")
tic = time.perf_counter()
for i in range(rep):
spec2 = fftpack.fft(in1,axis=0)
spec2 = fftpack.fft(spec2,axis=1)
toc = time.perf_counter()
print("", f"2D Spectrum 2xFFT (len = {len(spec2)}):",
f"spec2 takes {10**0*((toc - tic)/rep):0.4f} s", sep="\n")
tic = time.perf_counter()
for i in range(rep):
spec3 = fftpack.fft2(in1)
toc = time.perf_counter()
print("", f"2D Spectrum FFT2 (len = {len(spec3)}):",
f"spec3 takes {10**0*((toc - tic)/rep):0.4f} s", sep="\n")
# compare
print('\nIs spec1 equivalent to the spec2?', np.allclose(spec1, spec2))
print('\nIs spec2 equivalent to the spec3?', np.allclose(spec2, spec3), '\n')
矩阵大小 = 5x5 的结果
Input (len = 5):
2D Spectrum FFT with np.apply_along_axis (len = 5):
spec1 takes 0.000183 s
2D Spectrum 2xFFT (len = 5):
spec2 takes 0.000010 s
2D Spectrum FFT2 (len = 5):
spec3 takes 0.000012 s
Is spec1 equivalent to the spec2? True
Is spec2 equivalent to the spec3? True
矩阵大小 = 500x500 的结果
Input (len = 500):
2D Spectrum FFT with np.apply_along_axis (len = 500):
spec1 takes 0.017626 s
2D Spectrum 2xFFT (len = 500):
spec2 takes 0.005324 s
2D Spectrum FFT2 (len = 500):
spec3 takes 0.003528 s
Is spec1 equivalent to the spec2? True
Is spec2 equivalent to the spec3? True
矩阵大小 = 5000x5000 的结果
Input (len = 5000):
2D Spectrum FFT with np.apply_along_axis (len = 5000):
spec1 takes 2.538471 s
2D Spectrum 2xFFT (len = 5000):
spec2 takes 0.846661 s
2D Spectrum FFT2 (len = 5000):
spec3 takes 0.574397 s
Is spec1 equivalent to the spec2? True
Is spec2 equivalent to the spec3? True
结论
从上面的测试看来，使用 fftpack.fft2() 对于更大的矩阵更有效。
使用np.apply_along_axis()是最慢的方法。