python - How to multiply two 2D RFFT arrays (FFTPACK) to be compatible with NumPy's FFT?

Question

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:

• 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.

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.

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:

• 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;

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), '
')

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().

Only answers that address the problem in 2-dimensions are invited. Those interested in how to multiply 1D FFTPACK arrays can check this thread.

Answer 1

Correct functions:

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

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.

Moreover, you have input data as a numpy 2D arrays, why do you use list comprehension? Use fftpack.fft() directly, this is much faster.

• 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

================================================================

The full testing code with new functions version:

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), '
')

The output is:

####################     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