Slightly different FFT results from Matlab fft and Scipy fft

Question

I am doing a routine that measures the phase difference between two spectra using NumPy / Scipy.

I already had a program written in Matlab, so I basically reimplemented the function and the corresponding unit test with NumPy. However, I found that the unit test fails because scipy.fftpack.fft introduces 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 symmetrical, the expected output

 [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 

As a result, there should be no imaginary components based on the 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 scipy.fftpack.fft introduce small imaginary components? I really want to avoid this problem. Can anyone give me a suggestion?

Answer 1

On the one hand, scipy.fftpack.fft guaranteed to always return a complex result, while the result of the MATLAB fft function is sometimes large and sometimes complex, depending on the presence of a nonzero imaginary component. However, this does not explain why the result of scipy.fftpack.fft really contains non-zero imaginary components, while the result of the MATLAB fft function fft not work.

I suspect that the main reason for the difference is due to the fact that the MATLAB fft function is apparently based on FFTW , while scipy and numpy use FFTPACK due to licensing restrictions.

pyfftw , however, provides 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 is compared exactly equal to zero, while FFTPACK has a tiny number of round-off errors with floating point.

In addition, I have no idea why the FFTW implementation suffers less from rounding errors than FFTPACK, but in any case it is important to note that these rounding errors are small enough so that they usually do not cause problems (you know, t check for exact equality between float values, right?).

Usually you just accept the real component of the result, for example:

 scipy.fftpack.fft(x).real 

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