CODEWITHSUNDEEP
pythonnumpyfft
J Webster
J Webster

Reputation: 219

Unexpected behavior of numpy.fft.fft with high precision numbers

I have the following code...Note the two lines under # generate sine curve. One uses a higher precision value for 2pi than the other, they should still give near identical results though.

import numpy as np
import matplotlib.pyplot as plt


t1 = np.arange(0., 1., .01)

# generate sine curve
y1 = np.sin(6.28318*5.*t1)  
#y1 = np.sin(6.283185307179586*5.*t1) # equivalent to np.sin(2*np.pi*t1)

# calculate the fft (no averaging!) of the time series
ffty = np.fft.fft(y1)

fig, ax_list = plt.subplots(3,1)
ax_list[0].plot(t1,y1, '.-')

ax_list[1].plot(ffty.real, '.-', label='Real Part')
ax_list[1].legend()

ax_list[2].plot(ffty.imag, '.-', label='Imag Part')
ax_list[2].legend()


plt.show()

If you run the code with the lower precision 6.28318 you get the expected result for the fft... enter image description here

However, if you run the code with the higher precision 6.283185307179586 which is equal to 2.*numpy.pi, you get the unexpected result below... the real part is drastically wrong...The amplitudes are way off, it's not symmetric, it doesn't make any sense. enter image description here

I'm at a loss as to what is causing this. Anyone have any ideas?

Upvotes: 4

Views: 2001

Answers (2)

user2699
user2699

Reputation: 3147

As @Cris Luengo said you need to look at the scale of the y-axis to accurately compare two plots. Another way to do this is to plot both of the things you're trying to compare on the same figure, as I've done below.

The magnitude of the FFT is displayed, using a log scale, and it's quite evident that using fewer significant figures of pi does indeed result in a lower accuracy result. Most of the values aren't exactly zero, as is to be expected when using floating point numbers, but using more significant figures gives many orders of magnitude improvement, which is not immediately apparent when the FFTs are plotted separately. overlayed fft for various precisions

code used:

import numpy as np
import matplotlib.pyplot as plt


t1 = np.arange(0., 1., .01)

values = {
'low':6.28318,
'higher':6.283185307179586,
'highest':2*numpy.pi,
}
styles = {
'low':'-',
'higher':'-',
'highest':'.-'
}

fig, ax_list = plt.subplots(3,1)
for name, tau in values.items():
    y1 = np.sin(tau*5.*t1)
    ffty = np.fft.fft(y1)

    ax_list[0].plot(t1,y1, styles[name], label=name)
    ax_list[1].plot(abs(ffty.real), styles[name],label=name)
    ax_list[2].plot(abs(ffty.imag), styles[name], label=name)

[ax.legend() for ax in ax_list]
ax_list[0].set_title('time domain')
ax_list[1].set_title('real part')
ax_list[2].set_title('imaginary part')
ax_list[1].set_yscale('log')
ax_list[2].set_yscale('log')
plt.draw()

Upvotes: 6

Cris Luengo
Cris Luengo

Reputation: 60615

This is totally expected behavior. Computers use floating-point computations, which are inherently imprecise.

Note the y-axis for your real result. If no numerical inaccuracy existed, the real component would be identically 0. With your "higher precision" result, the real part is almost identical to 0 (1e-14 is very close to the precision of double-precision floats). With a lower precision, the real part becomes much larger (though still much, much smaller than the imaginary part). Because of the larger numbers, there is more structure as well (i.e. the error is not given by rounding errors, but by an actual feature of your input data, a period that is slightly shorter than ideal).

Upvotes: 7

Related Questions