How to multiply two 2D RFFT arrays (FFTPACK) to be compatible with NumPy's FFT?
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.
Solution 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('\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
Solution 2:
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:]])
Solution 3:
@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).