如何将两个 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).
I'm trying to multiply two 2D arrays that were transformed with fftpack_rfft2d() (SciPy's FFTPACK RFFT) and the result is not compatible with what I get from scipy_rfft2d() (SciPy's FFT RFFT).
下图共享脚本的输出，显示:
The image below shares the output of the script, which displays:
- 两个输入数组的初始化值；
- 使用
scipy_rfft2d()使用 SciPy 的 FFT 实现对 RFFT 进行变换后的两个数组，然后是使用scipy_irfft2d()向后变换后的乘法输出； - 使用 SciPy 的 FFTPACK 实现 RFFT 与
fftpack_rfft2d()和fftpack_irfft2d()相同的事情； - 使用
np.allclose()的测试结果，用于检查两个乘法的结果在使用各自的 IRFFT 实现转换回来后是否相同.
- The initialization values of both input arrays;
- Both arrays after they were transform with SciPy's FFT implementation for RFFT using
scipy_rfft2d(), followed by the output of the multiplication after its transformed backwards withscipy_irfft2d(); - The same things using SciPy's FFTPACK implementation for RFFT with
fftpack_rfft2d()andfftpack_irfft2d(); - The result of a test with
np.allclose()that checks if the result of both multiplications are the same after they were transformed back with their respective implementations for IRFFT.
为了清楚起见，红色矩形显示逆变换IRFFT后的乘法结果:左边的矩形使用SciPy的FFT IRFFT；右边的矩形，SciPy 的 FFTPACK IRFFT.当与 FFTPACK 版本的乘法固定时，它们应该呈现相同的数据.
Just to be clear, the red rectangles display the multiplication result after the inverse transform IRFFT: the rectangle on the left uses SciPy's FFT IRFFT; the rectangle on the right, SciPy's FFTPACK IRFFT. They should present the same data when the multiplication with the FFTPACK version is fixed.
我认为与 FFTPACK 版本的乘法结果不正确，因为 scipy.fftpack 返回的结果 RFFT 数组中的实部和虚部与 scipy.fft 中的 RFFT 不同强>:
I think the multiplication result with the FFTPACK version is not correct because scipy.fftpack returns the real and imaginary parts in the resulting RFFT array differently than the RFFT from scipy.fft:
- 我相信 scipy.fftpack 中的 RFFT 返回一个数组，其中一个元素包含实部，下一个元素包含其虚部；
- 在 scipy.fft 的 RFFT 中，每个元素都是一个复数，因此能够同时保存实部和虚部；
- I believe that RFFT from scipy.fftpack returns an array where one element contains the real part and the next element holds its imaginary counterpart;
- In RFFT from scipy.fft, each element is a complex number and therefore is able to hold the real and imaginary parts simultaneously;
如果我错了，请纠正我！我还想指出，由于 scipy.fftpack 不提供用于转换 2D 数组的函数，例如 rfft2() 和 irfft2()，我在下面的代码中提供了我自己的实现:
Please correct me if I'm wrong! I would also like to point out that since scipy.fftpack doesn't provide functions for transforming 2D arrays like rfft2() and irfft2(), I'm providing my own implementations in the code below:
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('
#################### INPUT DATA ###################
')
# 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('
in1 shape=', in1.shape, '
', in1)
in2 = np.array([[0, 0, 0, 0],
[0, 0, 255, 0],
[0, 255, 255, 0],
[0, 255, 0, 0]])
print('
in2 shape=', in2.shape, '
', in2)
print('
############### SCIPY: 2D RFFT (MULT) ###############
')
# 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, '
', scipy_rfft1.real)
print('
scipy_fft2 shape=', scipy_rfft2.shape, '
', 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('
* Output from scipy_fft.irfftn():')
print('scipy_data shape=', scipy_data.shape, '
', scipy_data)
print('
############### FFTPACK: 2D RFFT (MULT) ###############
')
# 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, '
', fftpack_rfft1)
print('
fftpack_rfft2 shape=', fftpack_rfft2.shape, '
', 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('
* Output from fftpack_irfft2d():')
print('fftpack_data shape=', fftpack_data.shape, '
', fftpack_data)
print('
##################### RESULT #####################
')
# compare FFTPACK result with SCIPY
print('
Is fftpack_data equivalent to scipy_data?', np.allclose(fftpack_data, scipy_data), '
')
假设我的猜测是正确的，那么对于一个将 fftpack_rfft2d() 生成的两个二维数组相乘的函数来说，正确的实现是什么?请记住:结果数组必须能够使用 fftpack_irfft2d() 转换回来.
Assuming my guess is correct, what would be the correct implementation for a function that multiplies two 2D arrays that were generated from fftpack_rfft2d()? Remember: the resulting array must be able to be transformed back with fftpack_irfft2d().
仅邀请解决二维问题的答案.那些对如何乘以一维 FFTPACK 数组感兴趣的人可以查看此线程.
Only answers that address the problem in 2-dimensions are invited. Those interested in how to multiply 1D FFTPACK arrays can check this thread.
推荐答案
正确的功能:
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.是的，第一个 FFT(在您的情况下按列)可以使用 rfft() 计算，但是第二个 FFT 计算必须在 complex 上提供strong> 第一个 FFT(按列) 的输出，因此 rfft() 的输出必须转换为真正的复频谱.此外，这意味着，您必须使用 fft() 而不是 rfft() 用于按行进行的第二个 FFT.因此，在两种计算中都使用 fft() 更方便.
You calculated the 2D FFT in wrong way. Yes, the first FFT (by columns in your case) can be calculated using rfft(), but the second FFT calculation must be provided on the complex output of the first FFT (by columns), so the output of the rfft() must be converted into true complex spectrum. Moreover, this mean, that you must use fft() instead of rfft() for the second FFT by rows. Consiquently, it is more convenient to use fft() in both calculations.
此外，您输入的数据是 numpy 二维数组，为什么要使用 列表推导?直接使用fftpack.fft()，这样会快很多.
Moreover, you have input data as a numpy 2D arrays, why do you use list comprehension? Use fftpack.fft() directly, this is much faster.
- 如果您已经只有由错误函数计算出的二维数组并且需要将它们相乘:那么，我认为，尝试使用相同的错误"方式从错误的二维 FFT 中重建输入数据，并且然后计算正确的2D FFT
- If you already have only 2D arrays calculated by wrong functions and need multiply them: then, my opinion, to try reconstruct the input data from the wrong 2D FFT using the same 'wrong' way and then calculate correct 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('
#################### INPUT DATA ###################
')
# 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('
in1 shape=', in1.shape, '
', in1)
in2 = np.array([[0, 0, 0, 0],
[0, 0, 255, 0],
[0, 255, 255, 0],
[0, 255, 0, 0]])
print('
in2 shape=', in2.shape, '
', in2)
print('
############### SCIPY: 2D RFFT (MULT) ###############
')
# 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, '
', scipy_rfft1)
print('
scipy_fft2 shape=', scipy_rfft2.shape, '
', 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('
* Output from scipy_fft.irfftn():')
print('scipy_data shape=', scipy_data.shape, '
', scipy_data)
print('
############### FFTPACK: 2D RFFT (MULT) ###############
')
# 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, '
', fftpack_rfft1)
print('
fftpack_rfft2 shape=', fftpack_rfft2.shape, '
', 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('
* Output from fftpack_irfft2d():')
print('fftpack_data shape=', fftpack_data.shape, '
', fftpack_data)
print('
##################### RESULT #####################
')
# compare FFTPACK result with SCIPY
print('
Is fftpack_data equivalent to scipy_data?', np.allclose(fftpack_data, scipy_data), '
')
输出为:
#################### 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
这篇关于如何将两个 2D RFFT 数组 (FFTPACK) 相乘以与 NumPy 的 FFT 兼容?的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
