from scipy import *
import Gnuplot
from scipy import fftpack
import math
import copy

y = [[], [], [], []]
power = copy.deepcopy(y)
period = copy.deepcopy(y)
p_data = copy.deepcopy(y)
nyquist=1./2

for i in range(32):
    y[0].append(math.sin(i/100.)
for i in range(1,4):
    y[i] = copy.deepcopy(y[0])

# randomised y
for i in xrange(32):
    y[1][i] += random.random()

# random inserts and doubled range
for i in range(1, 64, 2):
    if ( random.random() >= 0.85 ) and ( i <= len(y[3]) ):
        y[2].insert(i, y[2][i-1])
    y[3].insert(i, y[3][i-1])

for i in range(0, 4):
    Y = fft(y[i])
    ni = len(Y)/2
    nf = len(Y)/2.0
    freq = array(range(ni))/(nf)*nyquist

    power[i] = abs( Y[i][1:ni] )**2
    period[i] = 1./freq
    p_data[i] = Gnuplot.PlotItems.Data(period[i][1:len(period[i])], power[i],y[i])

Gnuplot.plot(p_data[0], p_data[1], p_data[2], p_data[3])
Gnuplot.xaxis((0,64))
Gnuplot.ytitle('|FFT|**2')
Gnuplot.grid("off")
Gnuplot.output('sunspot_period.png','png medium transparent picsize 600 400')