如何将两个 2D RFFT 数组 (FFTPACK) 相乘以与 NumPy 的 FFT 兼容？

Question

如何将两个 2D RFFT 数组 (FFTPACK) 相乘以与 NumPy 的 FFT 兼容？

kar*_*lip 10 python fft scipy fftpack numpy-ndarray

我正在尝试将用fftpack_rfft2d()（SciPy 的 FFTPACK RFFT）转换的两个二维数组相乘，结果与我从scipy_rfft2d()（SciPy 的 FFT RFFT）得到的结果不兼容。

下图共享了脚本的输出，其中显示：

两个输入数组的初始化值；
使用 SciPy 对 RFFT 的 FFT 实现进行变换后的两个数组scipy_rfft2d()，然后是乘法的输出，然后使用向后变换scipy_irfft2d()；
使用 SciPy 的 FFTPACK 实现 RFFT 与fftpack_rfft2d()和相同的事情fftpack_irfft2d()；
使用该测试np.allclose()的结果检查两个乘法的结果在使用其各自的 IRFFT 实现转换回来后是否相同。

为了清楚起见，红色矩形显示逆变换 IRFFT 后的乘法结果：左侧的矩形使用 SciPy 的 FFT IRFFT；右边的矩形，SciPy 的 FFTPACK IRFFT。当与 FFTPACK 版本的乘法固定时，它们应该显示相同的数据。

我认为与 FFTPACK 版本的乘法结果不正确，因为scipy.fftpack返回结果 RFFT 数组中的实部和虚部与来自scipy.fft的 RFFT 不同：

我相信来自scipy.fftpack 的 RFFT返回一个数组，其中一个元素包含实部，下一个元素包含其虚部；
在 scipy.fft 的RFFT 中，每个元素都是一个复数，因此能够同时保存实部和虚部；

如果我错了，请纠正我！我还想指出，由于scipy.fftpack不提供用于转换 2D 数组的函数，例如rfft2()and 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')

Run Code Online (Sandbox Code Playgroud)

假设我的猜测是正确的，那么将两个由生成的二维数组相乘的函数的正确实现是fftpack_rfft2d()什么？请记住：生成的数组必须能够使用fftpack_irfft2d().

仅邀请解决二维问题的答案。那些对如何乘以 1D FFTPACK 数组感兴趣的人可以查看此线程。

Answer 1

And*_*hei 4

正确的功能：

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

Run Code Online (Sandbox Code Playgroud)

您以错误的方式计算了 2D FFT。是的，可以使用rfft()计算第一个 FFT （在您的情况下按列），但第二个 FFT 计算必须在第一个 FFT （按列）的复数输出上提供，因此rfft()的输出必须被转换成真正的复谱。此外，这意味着您必须使用fft()而不是rfft()进行第二个按行 FFT。因此，在这两种计算中都使用fft()更方便。

此外，您的输入数据为numpy 2D 数组，为什么要使用列表理解？直接使用fftpack.fft()，这样快很多。

如果您已经只有由错误函数计算的 2D 数组，并且需要将它们相乘：那么，我认为，尝试使用相同的“错误”方式从错误的 2D FFT 重建输入数据，然后计算正确的 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')
Run Code Online (Sandbox Code Playgroud)
输出是：

#################### 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
Run Code Online (Sandbox Code Playgroud)

归档时间：	5 年，6 月前
查看次数：	722 次
最近记录：	5 年，6 月前