1
$\begingroup$

I'm trying to understand how Python's fft2 is different from irfft2. Basically, I have some continuous Fourier 2D samples that I need to put in a (n,n) grid (this can be done by a simply count-in-cell gridding or convolutional gridding). This results in the following uv-plane (named in the code uv_grid, see npy file):

Example using simulated gridded data

2D Fourier plane

If I do an IFFT2:

import numpy as np
image = np.fft.ifftshift(np.fft.ifft2(np.fft.fftshift(uv_grid)))

I get this image:

IFFT2

Now the question is - if I get rid from the redundant values of the gridded Fourier plane, can I get the same image?

First, I took the Fourier plane and I took half of the pixels in the x-axis and add a column such that shape is (n, n//2 + 1).

hermitian_symmetry_grid = np.hstack(
    [
        uv_grid[:, 0:n // 2],
        np.zeros((n, 1), dtype=np.complex64)
    ]
)
hermitian_symmetry_grid = np.flip(hermitian_symmetry_grid, axis=-1).conj()

Resulting in:

Hermitian

Now if I do

image_from_symmetry = np.fft.ifftshift(np.fft.irfft2(
    np.fft.fftshift(hermitian_symmetry_grid, axes=-2), s=(n,n)
))

I get this image:

IRFFT2 Image

which has clearly visible artifacts and is upside down and is not the same as the IFFT2 of the original grid. How can I get the same exact image?

Example using random data

import numpy as np

n = 256
image = np.random.random((n,n))

uv_grid = np.fft.fftshift(np.fft.fft2(np.fft.ifftshift(image)))

hermitian_symmetry_grid = np.hstack(
    [
        uv_grid[:, n // 2::],
        np.zeros((n, 1), dtype=np.complex64)
    ]
)

image_from_symmetry = np.fft.fftshift(np.fft.irfft2(
    np.fft.ifftshift(hermitian_symmetry_grid, axes=-2), s=(n,n))
)

np.allclose(image, image_from_symmetry)

(The last line returns False)

There was a bug in the gridding procedure that made the data not being properly hermitian-symmetric. It's all consistent after fixing that bug and using the @Cris Luengo code (see answer).

$\endgroup$
12
  • $\begingroup$ rfft is for real data since the spectra of real data is mirrored. 2D data has positive and negative frequencies in the vertical and horizontal direction, so my guess is that you would need to take a quarter plane of the 2D spectrum and use irfft2, probably starting at 0 in both directions. $\endgroup$
    – Baddioes
    Apr 15 at 17:11
  • $\begingroup$ You mean taking the fourth quadrant and putting it on the upper part of the (n,n//2 +1) matrix? What about the bottom part of the matrix? $\endgroup$ Apr 15 at 17:48
  • $\begingroup$ If you have an $M$-by-$N$ 2D FFT of a real 2D signal, I would assume irfft2 is expecting an $\frac{M}{2}$-by-$\frac{N}{2}$ input, as the other parts of the spectrum can be computed via combinations of conjugations in the vertical and/or horizontal dimension. $\endgroup$
    – Baddioes
    Apr 15 at 18:03
  • $\begingroup$ Nope, irfft2 expects a (m, n//2-1) input. $\endgroup$ Apr 15 at 18:06
  • $\begingroup$ Interesting, not super familiar with that function, so I'll have to look into how that works. If your spectrum is centered at DC, you need to ifftshift the spectrum, as a spectrum centered at DC is one that has been shifted via fftshift, and ifftshift inverts this process. You would then irfft2, and there's no need for another shift. Try that and if that doesn't work let me know. $\endgroup$
    – Baddioes
    Apr 15 at 18:39

2 Answers 2

1
$\begingroup$

Things are a lot simpler if we don't use so many fftshift calls.

The RFFT for real-valued inputs produces the left half of the general FFT. This includes the zero frequency at index 0, and the highest frequency at or below Nyquist at index n//2 (for odd-length signals, the frequency for index n//2 is below the Nyquist frequency, for even-length signals it is at the Nyquist frequency).

So to reproduce that from the output of fft2 is quite simple:

import numpy as np

n = 256
image = np.random.random((n,n))

uv_grid = np.fft.fft2(image)
assert np.allclose(image, np.fft.ifft2(uv_grid))

half_uv_grid = uv_grid[:, :n // 2 + 1]  # +1 because end index is excluded

image_from_half_uv_grid = np.fft.irfft2(half_uv_grid, s=(n,n))
assert np.allclose(image, image_from_half_uv_grid)
$\endgroup$
2
  • $\begingroup$ I have been able to work this out even with the fft and ifft shifts. However, if you try the same with the real data (which is the aim of the question) does not work. I would appreciate if you can give a try to that data. $\endgroup$ Apr 16 at 20:36
  • $\begingroup$ I will mark this as correct, since all the problem was in the gridding procedure. Cheers. $\endgroup$ Apr 16 at 21:43
0
$\begingroup$

Running the following code gets me quite close.

image = np.fft.ifftshift(np.fft.ifft2(np.fft.fftshift(uv_grid)))
n = image.shape[0]
RFFT = np.fft.fftshift(np.fft.rfft2(image))
aRFFT = np.abs(RFFT)

hermitian_symmetry_grid = np.hstack(
    [
        uv_grid[:,0:n//4],uv_grid[:, n // 2:3*n//4+1]
    ]
)

image_from_symmetry = np.fft.fftshift(np.fft.irfft2(np.fft.ifftshift(
    hermitian_symmetry_grid), s=image.shape
))

plt.figure(1)
imgplot = plt.imshow(np.real(image),cmap=None,vmin=-0.0025,vmax=0.0175)  # Display the magnitude of the complex image
plt.colorbar()  # Add a colorbar to the current figure
plt.title('Original Image')
plt.show()

plt.figure(2)
imgplot2 = plt.imshow(image_from_symmetry,cmap=None,vmin=-0.0025,vmax=0.0175)
plt.colorbar()
plt.title('Image from Hermitian Symmetry')
plt.show()

plt.figure(3)
imgplot3 = plt.imshow(np.abs(RFFT))
plt.colorbar()
plt.title('RFFT Spectrum')
plt.show()

plt.figure(4)
imgplot3 = plt.imshow(np.abs(hermitian_symmetry_grid))
plt.colorbar()
plt.title('Hermitian Symmetry Grid')
plt.show()

There is still some minor artifacting that occurs, which is hard to say what is causing that. It could be the irfft2 function, it could be finite precision effects, etc. However, this should get you much closer.

EDIT: See updated code

$\endgroup$
4
  • $\begingroup$ Thanks for your answer. However, this is not correct. The result from irfft2 should be real. And there should be no need to plot the absolute values. If you don't plot the abs the result is just a flat image. Additionally, I'm not sure if taking the left part of the complete plane is correct either. For example, try to do an rfft2 to an image and check where the data is (it will be on the left corners of the Fourier plane). In your example, the ifftshift will put the data on the right corners, which seems incorrect to me. $\endgroup$ Apr 16 at 15:35
  • $\begingroup$ @MiguelCárcamo see the updated code. This doesn't use abs() on the irfft2() image. This gets much closer, but there is still some minor artifacting as a result of irfft2(), which I haven't figured out. $\endgroup$
    – Baddioes
    Apr 16 at 18:36
  • $\begingroup$ I have checked this and the result is the same, with unfortunately same artifacts. $\endgroup$ Apr 16 at 19:19
  • $\begingroup$ @MiguelCárcamo I’m not really sure what to say. It shouldn’t be the same because what you had before was not in the output format of rfft2(), but I edited it so it now is. When I run the provided code, I get the images to be scaled much more closely, with a minor worsening in the artifacts using irfft2(). $\endgroup$
    – Baddioes
    Apr 16 at 19:56

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.