Inverse Fourier transform - porting from MATLAB to Python

Question

I am porting a script from MATLAB to Python, but I am failing when it comes to the inverse Fourier transform.

The function in MATLAB (ifft) includes a 'symflag', which treats the data as conjugate symmetric and ensures that the output is real. However, the equivalent function in numpy does not have a similar option, and therefore my outputs are different.

The 'symflag' clearly covers up for some of my mathematical inadequacies, and I would be grateful for assistance in resolving this issue. My aim would be to have a solution in Python that closely matches the solution in MATLAB.

In Matlab

Y = [complex(1), complex(2), complex(3), complex(3), complex(2), complex(1)]

y = ifft(Y, 'symmetric')

y = 2.33    -0.50   -0.16   0   -0.16   -0.50

In Python

Y = np.array([1, 2, 3, 3, 2, 1], dtype='complex')

y = ifft(Y)

y = ?

Many thanks in advance.

Answer 1

For a vector Y of length $N$, where $N$ is even, conjugate symmetry requires the following property to hold: $$ Y[k] = Y[N-k]^* $$ for $k=1,...,N-1$ and $Y_0$ is real. (See this.) Here $^*$ denotes the complex conjugate.

The array [1,2,3,3,2,1] this conjugate symmetry property doesn't hold. What did Matlab do when you used 'symmetric' with an input vector that wasn't conjugate symmetric? It pretends [1,2,3,3,2,1] is nearly conjugate symmetric$^\dagger$ and just replaces it with the closest conjugate symmetric vector, which is, [1,2,3,3,3,2].

If you try this in Python:

Y = np.array([1,2,3,3,3,2])
y = np.fft.ifft(Y)

your answer will match Matlab's answer with the 'symmetric' option for [1,2,3,3,2,1]. Note there there will be some rounding issues and you can verify that the imaginary part of y in Python is practically zero. Here's my output:

In [8]: np.fft.ifft(np.array([1,2,3,3,3,2]))
Out[8]: 
array([  2.33333333e+00 +0.00000000e+00j,
        -5.00000000e-01 -1.48029737e-16j,
        -1.66666667e-01 +0.00000000e+00j,
         1.11022302e-16 +2.90348171e-16j,
        -1.66666667e-01 +0.00000000e+00j,  -5.00000000e-01 -1.42318434e-16j])

In [9]: np.real(np.fft.ifft(np.array([1,2,3,3,3,2])))
Out[9]: 
array([  2.33333333e+00,  -5.00000000e-01,  -1.66666667e-01,
         1.11022302e-16,  -1.66666667e-01,  -5.00000000e-01])

In [10]: np.imag(np.fft.ifft(np.array([1,2,3,3,3,2])))
Out[10]: 
array([  0.00000000e+00,  -1.48029737e-16,   0.00000000e+00,
         2.90348171e-16,   0.00000000e+00,  -1.42318434e-16])

And now to your main question: How do you port this to Python?

Step 1: Take the input vector Y=[1,2,3,3,2,1]

Step 2: Force it to be conjugate symmetric by looping through and changing the right half of the vector to be equal to the complex conjugate of the left half. Call this modified vector Y_new = [1,2,3,3,3,2].

Step 3: Compute np.fft.ifft(Y_new) and throw away its imaginary part. (You can do an assert and ensure the imaginary part has a very tiny norm.)

Alternatively, when $N$ is even, Steps 2 and 3 can be replaced by a quicker np.fft.irfft(Y[0:N/2+1]) as endolith suggests in the comment below.

See also: https://blogs.mathworks.com/steve/2010/07/16/complex-surprises-from-fft/ for an example when $N$ is odd.

$^\dagger$ AFAIK, mathematically, there is no such thing is "nearly conjugate symmetric." A vector is either conjugate symmetric or not! I wish Matlab's documentation was more explicit about what it is doing under to hood.