如何将两个2D RFFT数组(FFTPACK)相乘以与NumPy的FFT兼容? [英] How to multiply two 2D RFFT arrays (FFTPACK) to be compatible with NumPy's FFT?

查看:97
本文介绍了如何将两个2D RFFT数组(FFTPACK)相乘以与NumPy的FFT兼容?的处理方法,对大家解决问题具有一定的参考价值,需要的朋友们下面随着小编来一起学习吧!

问题描述

我正试图将两个用fftpack_rfft2d()(SciPy的FFTPACK RFFT)转换的2D数组相乘,结果与我从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的scipy_rfft2d()用于RFFT的FFT实现进行了变换,然后是乘法的输出,然后通过scipy_irfft2d()向后变换了乘法的输出;
  • 使用SciPy的FFTPACK实现对具有fftpack_rfft2d()fftpack_irfft2d()的RFFT进行相同的操作;
  • 使用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 with scipy_irfft2d();
  • The same things using SciPy's FFTPACK implementation for RFFT with fftpack_rfft2d() and fftpack_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(),因此我在下面的代码中提供了自己的实现:/p>

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('\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()生成的两个2D数组相乘的函数的正确实现是什么?请记住:结果数组必须能够使用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.是的,可以使用 rfft()计算第一个FFT(按您所在的列),但是第二个FFT计算必须提供在复合物第一个FFT(按列)的strong>输出,因此必须将 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 2D数组,为什么要使用 list comprehension ?直接使用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.

  • 如果您仅具有由错误函数计算的2D数组并需要将它们相乘:然后,我认为,尝试使用相同的错误"方法从错误的2D 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('\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 

这篇关于如何将两个2D RFFT数组(FFTPACK)相乘以与NumPy的FFT兼容?的文章就介绍到这了,希望我们推荐的答案对大家有所帮助,也希望大家多多支持IT屋!

查看全文
登录 关闭
扫码关注1秒登录
发送“验证码”获取 | 15天全站免登陆