Slightly different FFT results from Matlab fft and Scipy fft

Question

I've been making a routine which measures the phase difference between two spectra using NumPy/Scipy.

I already had the routine written in Matlab, so I basically re-implemented the function and the corresponding unit test using NumPy. However, I found that the unit test fails because scipy.fftpack.fft is introducing some small numerical errors:

import numpy as np
import scipy.fftpack.fft
x = np.array([0.0, 1.0, 2.0, 3.0, 4.0, 3.0, 2.0, 1.0])
X = scipy.fftpack.fft(x)

In this case, since the time-domain signal is symmetric, the expected output is

[16.0000   -6.8284         0   -1.1716         0   -1.1716         0   -6.8284]

as shown in the following Matlab code:

>> x = [0.0, 1.0, 2.0, 3.0, 4.0, 3.0, 2.0, 1.0];
>> X = fft(x)

X =

   16.0000   -6.8284         0   -1.1716         0   -1.1716         0   -6.8284

The result should not contain any imaginary components, based on DSP theory. However, the scipy result is as follows:

array([ 16.00000000 +0.00000000e+00j,  -6.82842712 -2.22044605e-16j,
         0.00000000 -0.00000000e+00j,  -1.17157288 -2.22044605e-16j,
         0.00000000 +0.00000000e+00j,  -1.17157288 +2.22044605e-16j,
         0.00000000 +0.00000000e+00j,  -6.82842712 +2.22044605e-16j])

Why does scipy.fftpack.fft introduce small imaginary components? I really want to avoid this issue. Could anyone give me a suggestion?

Answer 1

For one thing, scipy.fftpack.fft is guaranteed to always return a complex result, whereas the result of MATLAB's fft function is sometimes real and sometimes complex, depending on whether there is a non-zero imaginary component. However, that doesn't explain why the result of scipy.fftpack.fft actually contains non-zero imaginary components, whereas the result of MATLAB's fft function does not.

I suspect that the underlying reason for the difference has to do with the fact that MATLAB's fft function is apparently based on FFTW, whereas scipy and numpy use FFTPACK due to licensing restrictions.

pyfftw, however, does provide Python bindings to FFTW. If we compare the imaginary components of the results for FFTPACK and FFTW:

from pyfftw.interfaces import scipy_fftpack as fftw

Fx1 = fftpack.fft(x)
print(Fx1.imag)
# [  0.00000000e+00  -2.22044605e-16  -0.00000000e+00  -2.22044605e-16
#    0.00000000e+00   2.22044605e-16   0.00000000e+00   2.22044605e-16]
print(Fx1.imag == 0)
# [ True False  True False  True False  True False]

Fx2 = fftw.fft(x)
print(Fx2.imag)
# [ 0.  0.  0.  0.  0.  0.  0.  0.]
print(Fx2.imag == 0)
# [ True  True  True  True  True  True  True  True]

we see that the imaginary component of the FFTW result compares exactly equal to zero, whereas FFTPACK has a tiny amount of floating-point rounding error.

Beyond that, I have no idea why FFTW's implementation suffers less from rounding error than FFTPACK, but in any case it's important to note that these rounding errors are small enough that they normally don't cause problems (you know you shouldn't be testing for exact equality between float values, right?).

Usually you would simply take the real component of the result, e.g.:

scipy.fftpack.fft(x).real

If these errors are a problem then you could switch to using pyfftw instead of numpy/scipy, but if your code is that sensitive to rounding error then it probably means you're doing something wrong anyway.