Friday, December 24, 2010
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.
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()
Subscribe to:
Posts (Atom)