I recoded parts of the Fourier analysis example that uses sunspots data found
here. It needed updating for most recent Python. All graphs match the original output.
import numpy as np
from scipy import fft, array
from matplotlib import pyplot as gplt
tempdata = np.loadtxt("sunspots.dat")
year=tempdata[:,0]
wolfer=tempdata[:,1]
gplt.plot(year,wolfer)
gplt.show()
gplt.plot(year[0:50],wolfer[0:50])
gplt.show()
Y=fft(wolfer)
gplt.plot(Y.real,Y.imag, 'r+')
gplt.xlim(-4000,2000)
gplt.ylim(-4000,4000)
gplt.show()
n=len(Y)
power = np.abs(Y[1:(n/2)])**2
nyquist=1./2
freq=array(range(n/2))/(n/2.0)*nyquist
freq += 0.00001
period=1./freq
gplt.plot(1/period[1:len(period)], power)
gplt.xlabel('cycles/year')
gplt.title('Periodogram')
gplt.show()
gplt.plot(1/period[1:40], power[1:40])
gplt.show()
gplt.plot(period[1:len(period)], power)
gplt.xlim(0,40)
gplt.show()
No comments:
Post a Comment