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()

No comments:

Post a Comment