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