#! /usr/bin/python
##! /Users/dhw/anaconda/bin/python
# fftnoise.py -- generate Gaussian noise and measure/plot its FFT
# fftnoise.py n amp
#   n = length of array
#   amp = rms amplitude of the noise

import sys
import numpy as np
import matplotlib.pyplot as plt
from math import *

np.random.seed(12)

n=int(sys.argv[1])
amp=float(sys.argv[2])

x=np.arange(n)/float(n)
y=np.random.normal(scale=amp,size=n)

z=np.fft.fft(y)/n
#z=np.fft.rfft(y)/n		# uncomment for real-to-complex fft
f=np.arange(z.size)
zr=np.real(z)
zi=np.imag(z)

modz=np.absolute(z)
meansqr1=np.sum(y**2)/n
# for complex fft
meansqr2=np.sum(modz**2)
# for real fft
# meansqr2=2*np.sum(modz**2)-(modz[modz.size-1])**2-(modz[0])**2
rms1=sqrt(meansqr1)
rms2=sqrt(meansqr2)
print rms1, rms2

plt.subplot(2,1,1)
plt.plot(x,y)
plt.plot(x,y,'bo')

plt.subplot(2,1,2)
plt.plot(f,zr,'ro',f,zi,'b^')

plt.show()
plt.savefig('fftnoise2.eps')