How to multiply two 2D RFFT arrays (FFTPACK) to be compatible with NumPy's FFT?
Your hypothesis is correct. FFTPACK returns all coefficients in a single real vector in the format
[y(0),Re(y(1)),Im(y(1)),...,Re(y(n/2))] if n is even
[y(0),Re(y(1)),Im(y(1)),...,Re(y(n/2)),Im(y(n/2))] if n is odd
where scipy.rfft returns a complex vector
[y(0),Re(y(1)) + 1.0j*Im(y(1)),...,Re(y(n/2) + 1.0j*Im(y(n/2)))]
so you need to form a vector using the proper stride, as follows:
y_fft = np.cat([y_fftpack[0], y_fftpack[1:2:] + 1.0j*y_fftpack[2:2:]])
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('\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')
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
@Andrei is right: it's much simpler to just use the complex-valued FFT (though his implementation is unnecessarily complicated, just use scipy.fftpack.fft2). As I said in a comment, the best option is to switch over to scipy.fft, which is simpler to use; fftpack is deprecated in favour of it.
However, if you do need to use fftpack, and you do want to save some computational time by using the rfft function, then this is the right way to do so. It requires converting the real-valued output of the rfft function to a complex-valued array before computing the fft along the other dimension. With this solution, fftpack_rfft2d below outputs half the 2D FFT of its input, with the other half being redundant.
import numpy as np
from scipy import fftpack
# FFTPACK RFFT 2D
def fftpack_rfft1d(matrix):
assert not (matrix.shape[1] & 0x1)
tmp = fftpack.rfft(matrix, axis=1)
assert tmp.dtype == np.dtype('float64')
return np.hstack((tmp[:, [0]], np.ascontiguousarray(tmp[:, 1:-1]).view(np.complex128), tmp[:, [-1]]))
def fftpack_rfft2d(matrix):
return fftpack.fft(fftpack_rfft1d(matrix), axis=0)
# FFTPACK IRFFT 2D
def fftpack_irfft1d(matrix):
assert matrix.dtype == np.dtype('complex128')
tmp = np.hstack((matrix[:, [0]].real, np.ascontiguousarray(matrix[:, 1:-1]).view(np.float64), matrix[:, [-1]].real))
return fftpack.irfft(tmp, axis=1)
def fftpack_irfft2d(matrix):
return fftpack_irfft1d(fftpack.ifft(matrix, axis=0))
######
# test data
in1 = np.random.randn(256,256)
in2 = np.random.randn(256,256)
# fftpack.fft2
gt_result = fftpack.ifft2(fftpack.fft2(in1) * fftpack.fft2(in2)).real
# fftpack_rfft2d
our_result = fftpack_irfft2d(fftpack_rfft2d(in1) * fftpack_rfft2d(in2) )
# compare
print('\nIs our result equivalent to the ground truth?', np.allclose(gt_result, our_result), '\n')
[This code only works for even-sized images, I didn't bother making it generic, see here for how to do so).
Nonetheless, as this solution requires copies of the data, it's actually slower than just using a normal, complex-valued FFT (fftpack.fft2), even though it does fewer computations:
import time
tic = time.perf_counter()
for i in range(100):
fftpack.fft(in1)
toc = time.perf_counter()
print(f"fftpack.fft() takes {toc - tic:0.4f} seconds")
tic = time.perf_counter()
for i in range(100):
fftpack_rfft2d(in1)
toc = time.perf_counter()
print(f"fftpack_rfft2d() takes {toc - tic:0.4f} seconds")
outputs:
fftpack.fft() takes 0.0442 seconds
fftpack_rfft2d() takes 0.0664 seconds
So, indeed, stick to fftpack.fft (or rather scipy.fft.fft if you can).