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()
Saturday, December 18, 2010
Fourier Transform and Sunspots Data
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.
Subscribe to:
Post Comments (Atom)
No comments:
Post a Comment