N = 1000
om = 2*pi*arange(1,N/2+1)/N
t = arange(0,N)
a=1./arange(1,N+1)
phi = rand(N)*2.*np.pi
x=0.
for i in range(1,N/2):
    kern=np.real(np.exp(1j*(t*om[i]+phi[i])))

    x = x+(a[i]**2)*kern
xf=x
figure(figsize=(6,8))
x=0.
for i in range(1,N/2):
    kern=np.real(np.exp(1j*(t*om[i]+phi[i])))

    x = x+(a[i]**2)*kern
    if i<=6:
        subplot(3,2,i)
        plot(x,label='Fit so far')
        plot((a[i]**2)*kern,label='Fit %d'%i)
        plot(xf,label='Signal')
        legend(loc=4,fontsize='small')
        ylim([-0.9,0.5])
savefig('images/FourierFit.png')

plot(x)

loglog(a)