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

by

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

More Related Contents:

  1. Peak-finding algorithm for Python/SciPy
  2. Plotting a fast Fourier transform in Python
  3. Python Scipy FFT wav files
  4. Convolving image with kernel in Fourier domain
  5. Improving FFT performance in Python
  6. Get frequency with highest amplitude from FFT
  7. What is the difference between flatten and ravel functions in numpy?
  8. Computing the correlation coefficient between two multi-dimensional arrays
  9. How does condensed distance matrix work? (pdist)
  10. ImportError in importing from sklearn: cannot import name check_build
  11. binning data in python with scipy/numpy
  12. Calculating Pearson correlation and significance in Python
  13. Calculating the area under a curve given a set of coordinates, without knowing the function
  14. Sparse matrix slicing using list of int
  15. What do all the distributions available in scipy.stats look like?
  16. Python 3D polynomial surface fit, order dependent
  17. Efficiently create sparse pivot tables in pandas?
  18. Map each list value to its corresponding percentile
  19. Creating a Confidence Ellipses in a sccatterplot using matplotlib
  20. scipy, lognormal distribution – parameters
  21. Most efficient way to calculate radial profile
  22. Efficient distance calculation between N points and a reference in numpy/scipy
  23. Distance between point and a line (from two points)
  24. Python/Scipy 2D Interpolation (Non-uniform Data)
  25. Zero pad numpy array
  26. Fitting a Weibull distribution using Scipy
  27. Numpy transpose of 1D array not giving expected result
  28. Fit mixture of two gaussian/normal distributions to a histogram from one set of data, python
  29. How to remove repititve pattern from an image using FFT
  30. Speedup scipy griddata for multiple interpolations between two irregular grids

Leave a Comment