# How is the inverse Fourier transform of a Hermitian-symmetric input efficiently calculated?

The np.fft.rfft2 function exploits the Hermitian symmetry of the transform of a real input to more efficiently calculate the transform by omitting the negative frequency terms along the 1 axis (i.e., as opposed to the output of np.fft.fft2, which is ny$$\times$$nx for input ny$$\times$$nx, the output is ny$$\times$$(nx//2 + 1)). The output is then the left half of the Hermitian symmetric matrix representing the transform.

To efficiently transform back, the np.fft.irfft2 treats the input as the output of np.fft.rfft2. Indeed, if the input to np.fft.irfft2 is the left half of a Hermitian symmetric matrix, the output will be the same as if the full matrix were input to np.fft.ifft2:

import numpy as np
np.random.seed(1)
# randomly generate the left side of a Hermitian symmetric matrix
ny = 3
nx = 6 # must be even
h = np.fft.rfft2(np.random.standard_normal((ny, nx)))
h.round()
Out[2]:
array([[-3.+0.j,  0.-2.j,  3.-4.j, 16.+0.j],
[-1.-2.j,  3.+1.j, -1.-5.j,  1.-3.j],
[-1.+2.j,  0.+1.j,  4.+1.j,  1.+3.j]])
# get the right side and concatentate to get the full form
h_full = np.concatenate([h, np.flip(h[range(0, -ny, -1), 1:-1], axis=1).conj()], axis=1)
h_full.round()
Out[3]:
array([[-3.+0.j,  0.-2.j,  3.-4.j, 16.+0.j,  3.+4.j,  0.+2.j],
[-1.-2.j,  3.+1.j, -1.-5.j,  1.-3.j,  4.-1.j,  0.-1.j],
[-1.+2.j,  0.+1.j,  4.+1.j,  1.+3.j, -1.+5.j,  3.-1.j]])
# invert the partial form and the full form; show that the result is the same
np.allclose(np.fft.irfft2(h), np.fft.ifft2(h_full))
Out[4]: True


Unsurprisingly, if the input is not the left half of a Hermitian symmetric matrix, the same statement cannot be made, since the matrix cannot be completed into a Hermitian symmetric one:

# randomly generate the left side of an imaginary matrix
h = 1.j * np.random.standard_normal((ny, nx // 2 + 1)) + np.random.standard_normal((ny, nx // 2 + 1))
h.round()
Out[5]:
array([[-1.+0.j, -0.+1.j, -1.-1.j, -1.+1.j],
[-1.+1.j, -0.+1.j, -1.+1.j,  0.-1.j],
[ 2.-0.j,  1.-1.j, -0.-0.j, -1.+1.j]])
# get the right side the same way as before
h_full = np.concatenate([h, np.flip(h[range(0, -ny, -1), 1:-1], axis=1).conj()], axis=1)
h_full.round()
Out[6]:
array([[-1.+0.j, -0.+1.j, -1.-1.j, -1.+1.j, -1.+1.j, -0.-1.j],
[-1.+1.j, -0.+1.j, -1.+1.j,  0.-1.j, -0.+0.j,  1.+1.j],
[ 2.-0.j,  1.-1.j, -0.-0.j, -1.+1.j, -1.-1.j, -0.-1.j]])
# invert the partial form and the full form; show that the result is not the same
np.allclose(np.fft.irfft2(h), np.fft.ifft2(h_full))
Out[7]: False



How is the output calculated in this case?

• is that result defined meaningfully? If you have an input that does not belong to the domain of a mapping, what's the correct output? Commented Nov 12, 2022 at 9:51
• @MarcusMüller Non-Hermitian symmetric input is in the domain of np.fft.irfft2 - otherwise, the mapping wouldn't be defined and I would get an error. I would like to know how the inverse real transform operates on Hermitian symmetric input, which will let me know how numpy handles the non-Hermitian symmetric case. Commented Nov 12, 2022 at 17:26
• no, it isn't. Just because something doesn't throw an error doesn't mean it's sensible input, or in the domain of the mathematical function that is numerically implemented. Mathematically, your rfft2 is a mapping to Hermitian matrices. Hence, the inverse is only defined for Hermitian matrices. This really has nothing to do with Python exceptions being thrown – just logic. Commented Nov 12, 2022 at 17:28
• @MarcusMüller This question is not about the mathematical function that is being implemented, but the implementation of the algorithm itself. The algorithm accepts non-Hermitian symmetric matrices as input (although the mapping is not documented), and so non-Hermitian matrices are in the domain of the function in question Commented Nov 12, 2022 at 23:29
• Ah I might have simply misunderstood you then. So you're asking about the undocumented behavior? (I guess all we could do is look at the source code, and learn from that which part of the matrix is used. ) (I'd then also point out that depending on undocumented behavior is a recipe for pain once that behavior changes) Commented Nov 12, 2022 at 23:56

The 1D IRFFT (np.fft.irfft()) simply ignores the imaginary component of the first element, and also the imaginary component of the last element if the output size is even. This is the bit of source code where that happens. Note how the complex data is treated as real data with alternating real and imaginary components, then element 0 is copied, and n-1 elements starting at 2 (the real part of the second element) are copied. n here is the number of output values, which means that in the even case, the imaginary component of the last element is not included. [To appreciate this, you might want to notice that for the even case, the input has n+2 elements ((n+2)/2 complex values), whereas for the odd case, the input has only n+1 elements ((n+1)/2 complex values).]

The 2D IRFFT (np.fft.irfft2()) does a complex-valued IFFT over the first axis, and then a 1D IRFFT over the second axis. You can replicate its working in terms of complex IFFTs as follows:

import numpy as np

np.random.seed(1)

# randomly generate the left side of a Hermitian symmetric matrix
ny = 3
nx = 6 # must be even
h = 1.j * np.random.standard_normal((ny, nx // 2 + 1)) + np.random.standard_normal((ny, nx // 2 + 1))
h.round(2)

# 1D IFFT along first axis
p1 = np.fft.ifft(h, axis=0)
p1.round(2)

# Set complex elements of first and last column to 0
p1[:, 0] = p1[:, 0].real
p1[:, -1] = p1[:, -1].real
p1.round(2)

# Expand 2nd axes trough complex conjugate symmetry
p2 = np.concatenate([p1, np.flip(p1[:, 1:-1], axis=1).conj()], axis=1)
p2.round(2)

# 1D IFFT along second axis
p3 = np.fft.ifft(p2, axis=1)
p3.round(2)

# Take real component
p4 = p3.real
p4.round(2)

g = np.fft.irfft2(h, s=(ny, nx))
g.round(2)


If you run the code, you can see that the imaginary component of p3 is approximately 0, as expected. Then you see p4:

array([[ 0.13, -0.16,  0.36,  0.16, -0.76, -0.26],
[ 0.02,  0.57, -0.43,  0.05, -0.08, -0.19],
[-0.14, -0.2 , -0.27,  0.43,  0.46, -0.01]])


and g:

array([[ 0.13, -0.16,  0.36,  0.16, -0.76, -0.26],
[ 0.02,  0.57, -0.43,  0.05, -0.08, -0.19],
[-0.14, -0.2 , -0.27,  0.43,  0.46, -0.01]])


So how to construct h_full from h to match the results? This is a bit more complicated. One thing to note is that the first column of h needs to be complex conjugate symmetric, and the last column also if the rows are even in size. Because we accomplish this by setting the imaginary component to 0 in its IFFT, we are effectively removing the odd component of the signal. So, split the column in an even and an odd part, and discard the odd part:

h[:, 0] = (h[:, 0] + h[range(0, -ny, -1), 0].conj()) / 2


Let's try this:

q1 = h.copy()
q1[:, 0] = (q1[:, 0] + q1[range(0, -ny, -1), 0].conj()) / 2
q1[:, -1] = (q1[:, -1] + q1[range(0, -ny, -1), -1].conj()) / 2  # because nx is even
q1.round(2)

h_full = np.concatenate([q1, np.flip(q1[range(0, -ny, -1), 1:-1], axis=1).conj()], axis=1)
h_full.round(2)

np.allclose(np.fft.irfft2(h, s=(ny, nx)), np.fft.ifft2(h_full))


...which returns True.

PS: Note that irfft2() needs the s parameter specifying the output size. If not given, the assumption is that the output is even in size, but it is best to be explicit.

• Interesing! How would be the inverse process? That is - taking the fft2 of an image, and the trying to take half of the Fourier plane to get same result in irfft2? I'm currently asking this question here Commented Apr 16 at 16:34