#**************************************************************
#* This script generates a PWL file for simulation.	      *	
#* The file is a pseudo-Gaussian Noise piecewise-linear	      *	
#* function for approximating a noisy resistor in simulation  *
#**************************************************************
import os
import sys
import math 
import random
import numpy as np
import matplotlib.pyplot as plt
import scipy
import scipy.fftpack
from pylab import *
from scipy import pi
graph_ON=False
default_update=False
lowpass_ON=False
highpass_ON=False
settings_Text_List=['-1','-1','-1','-1','-1','-1','-1','-1','-1']
default_List={'TEMP':-1,'RES':-1,'FREQ_LOW':-1,'FREQ_HIGH':-1,'SimSTART':-1,'SimEND':-1,'LowPASS':-1,'HighPASS':-1,'TimeMULT':-1}
update_List={'TEMP':-1,'RES':-1,'FREQ_LOW':-1,'FREQ_HIGH':-1,'SimSTART':-1,'SimEND':-1,'LowPASS':-1,'HighPASS':-1,'TimeMULT':-1}
#Initialize variables
KB=1.3806488e-23 #Boltzman Constant
H=6.62e-34
#Modifiable variables
TEMP=300.3 #Temp in Kelvin
RES=1e3 #Resistor Value
FREQ_LOW=1e3 #Lower Bound of Bandwidth
FREQ_HIGH=1000e9 #Upper Bound of Bandwidth
SimSTART=0 #Time start of simulation
SimEND=1 #Time end of simulation
TimeMULT=1e-15 #Magnitude of time step.
Sum_Of_Squares=0 #This variable will be used to check that random data is correct. Sanity Check.
LowPASS = 0


#################################
# Check for user command inputs #
#################################
if len(sys.argv)>1:
	print "-------------------"
	print "User Defined Inputs"
	print "-------------------\n"
	for i,cmd in enumerate(sys.argv):
		if cmd=='-g':
			print 'Graphing Enabled'
			graph_ON=True
		elif cmd=='-lpe':
			print 'Low-Pass Filter Enabled'
			lowpass_ON=True
		elif cmd=='-hpe':
			print 'High-Pass Filter Enabled'
			highpass_ON=True
		elif cmd=='-s':
			default_update=True
		
		elif cmd=='-t':
			try:
				TEMP = float(sys.argv[i+1])
				update_List['TEMP'] = TEMP
				print 'Temp=' + str(TEMP)
				settings_Text_List[0]=str(TEMP)+'\n'
				#break
			except IndexError:
				print 'ERROR: When using t command, a temp (K) must follow!...'
				exit()
		elif cmd=='-r':
			try:
				RES = float(sys.argv[i+1])
				print 'Res=' + str(RES)
				update_List['RES'] = RES
				settings_Text_List[1]=str(RES)+'\n'
			except IndexError:
				print 'ERROR: When using r command, a resistance value must follow!...'
				exit()
		elif cmd=='-fl':
			try:
				FREQ_LOW = float(sys.argv[i+1])
				print 'Low Frequency=' + str(FREQ_LOW)
				update_List['FREQ_LOW'] = FREQ_LOW
				settings_Text_List[2]=str(FREQ_LOW)+'\n'
			except IndexError:
				print 'ERROR: When using fl command, a frequency value in Hz must follow!...'
				exit()
		elif cmd=='-fh':
			try:
				FREQ_HIGH = float(sys.argv[i+1])
				print 'High Frequency=' + str(FREQ_HIGH)
				update_List['FREQ_HIGH'] = FREQ_HIGH
				settings_Text_List[3]=str(FREQ_HIGH)+'\n'
			except IndexError:
				print 'ERROR: When using fh command, a frequency value in Hz must follow!...'
				exit()
		elif cmd=='-ss':
			try:
				SimSTART = float(sys.argv[i+1])
				print 'Sim Start Time=' + str(SimSTART)
				update_List['SimSTART'] = SimSTART
				settings_Text_List[4]=str(SimSTART)+'\n'
			except IndexError:
				print 'ERROR: When using ss command, a time value in femto-seconds must follow!...'
				exit()
		elif cmd=='-se':
			try:
				SimEND = float(sys.argv[i+1])
				print 'Sim End Time=' + str(SimEND)
				update_List['SimEND'] = SimEND
				settings_Text_List[5]=str(SimEND)+'\n'
			except IndexError:
				print 'ERROR: When using se command, a time value in femto-seconds must follow!...'
				exit()
		elif cmd=='-lp':
			try:
				LowPASS = float(sys.argv[i+1])
				print 'Low Pass Filter -3dB frequency=' + str(LowPASS) + 'Hz'
				update_List['LowPASS'] = LowPASS
				settings_Text_List[6]=str(LowPASS)+'\n'
			except IndexError:
				print 'ERROR: When using lp command, a -3dB frequency in Hz must be specified!...'
				exit()
		elif cmd=='-hp':
			try:
				HighPASS = float(sys.argv[i+1])
				print 'High Pass Filter -3dB frequency=' + str(HighPASS) + 'Hz'
				update_List['HighPASS'] = HighPASS
				settings_Text_List[7]=str(HighPASS)+'\n'
			except IndexError:
				print 'ERROR: When using hp command, a -3dB frequency in Hz must be specified!...'
				exit()
		elif cmd=='-ts':
			try:
				TimeMULT = float(sys.argv[i+1])
				update_List['TimeMULT'] = TimeMULT
				print 'Time Step=' + str(TimeMULT)
				settings_Text_List[8]=str(TimeMULT)+'\n'
				#break
			except IndexError:
				print 'ERROR: When using ts command, a time step (s) must follow!...'
				exit()	

		elif cmd=='-h':
			print "Usage: python gauss.py [-g] [-lpe] [-hpe] [-t temperature(K)] [-r resistance(ohms)] [-fl Min Band(Hz)] [-fh Max Band(Hz)] [-ss Sim Start(fS)] [-se Sim End(fS)] [-lp Low-Pass fc(Hz)] [-hp Low-Pass fc(Hz)] [-ts Time Step (s)]" 			
			print "----------------------------------"
			print "| -g   | Enable Graphing         |"
			print "| -lpe | Enable Low Pass Filter  |"
			print "| -hpe | Enable High Pass Filter |"
			print "----------------------------------"
			exit()
############################
# Read in default settings #
#			   #	
# Open gauss.njn and read/ #
# modify settings          #
############################

#File structure line by line:  'TEMP','RES','FREQ_LOW','FREQ_HIGH','SimSTART','SimEND','LowPASS','HighPASS','TimeMULT'
f_settings_file = open('gauss.njn','r')
f_settings=f_settings_file.readlines()
f_settings_file.close()
for key, value in update_List.items():
	if value==-1: #Value not defined. Pull from file
		if key=='TEMP':
			update_List['TEMP']=float(f_settings[0])
			TEMP=float(f_settings[0])
			print "Temp set default:"+str(update_List['TEMP']) + "K"
			settings_Text_List[0]=str(update_List['TEMP'])+'\n'
		elif key=='RES':
			update_List['RES']=float(f_settings[1])
			RES=float(f_settings[1])
			print "RES set default:"+str(update_List['RES']) + " ohms"
			settings_Text_List[1]=str(update_List['RES'])+'\n'
		elif key=='FREQ_LOW':
			update_List['FREQ_LOW']=float(f_settings[2])
			FREQ_LOW=float(f_settings[2])
			print "Lower Bandwidth set default:"+str(update_List['FREQ_LOW']) + " Hz"
			settings_Text_List[2]=str(update_List['FREQ_LOW'])+'\n'			
		elif key=='FREQ_HIGH':
			update_List['FREQ_HIGH']=float(f_settings[3])
			FREQ_HIGH=float(f_settings[3])
			print "Upper Bandwidth set default:"+str(update_List['FREQ_HIGH']) + " Hz"
			settings_Text_List[3]=str(update_List['FREQ_HIGH'])+'\n'
		elif key=='SimSTART':
			update_List['SimSTART']=float(f_settings[4])
			SimSTART=float(f_settings[4])
			print "Simulation Start set default:"+str(update_List['SimSTART'])+" fS"
			settings_Text_List[4]=str(update_List['SimSTART'])+'\n'
		elif key=='SimEND':
			update_List['SimEND']=float(f_settings[5])
			SimEND=float(f_settings[5])	
			print "Simulation End set default:"+str(update_List['SimEND'])+" fS"
			settings_Text_List[5]=str(update_List['SimEND'])+'\n'
		elif key=='LowPASS':
			update_List['LowPASS']=float(f_settings[6])
			LowPASS=float(f_settings[6])	
			print "Low Pass Filter Set -3dB frequency set default:"+str(update_List['LowPASS'])+" Hz"
			settings_Text_List[6]=str(update_List['LowPASS'])+'\n'
		elif key=='HighPASS':
			update_List['HighPASS']=float(f_settings[7])
			HighPASS=float(f_settings[7])	
			print "High Pass Filter Set -3dB frequency set default:"+str(update_List['HighPASS'])+" Hz"
			settings_Text_List[7]=str(update_List['HighPASS'])+'\n'
		elif key=='TimeMULT':
			update_List['TimeMULT']=float(f_settings[8])
			TimeMULT=float(f_settings[8])	
			print "Time Step set default:"+str(update_List['TimeMULT'])+" s"
			settings_Text_List[8]=str(update_List['TimeMULT'])+'\n'
if default_update:
	print "Updating defaults"
	f_settings = open('gauss.njn','w')
	f_settings.writelines(settings_Text_List)
	f_settings.close()

################################
# Generate Pseudo Random Noise #
################################
#Standard Deviation for Thermal Voltage
#Technically Mean is 0 for thermal noise because it fluctuates around 0.
MEAN=0
SSPSD=4*KB*TEMP*RES
SSPSD_sqrt=math.sqrt(SSPSD)
STD_DEV=math.sqrt(4*KB*TEMP*RES*(FREQ_HIGH-FREQ_LOW))

print('************************************')
print('* Thermal Noise PWL File Generator *')
print('************************************')
print('Mean='+str(MEAN)+'; Standard Deviation='+str(STD_DEV)+'\n\n')
print('vrms per root Hz='+str(SSPSD_sqrt)+' V/sqrt(Hz)')
print('vrms per Hz='+str(SSPSD)+' V/Hz')
print('vrms='+str(STD_DEV)+ 'V')
raw_input("Press Enter to continue...")
print('Generating PWL file...')

#Generate Noise Values using Gaussian Distribution
noise_LIST=np.random.normal(0, STD_DEV, (SimEND+1-SimSTART))

######################
# Optional Filtering #
#		     #	
# IIR Single Pole    #
######################

if lowpass_ON:
#Low pass filter random signals
	sig_LP=[]
	ts=TimeMULT
	fs=1/ts
	fcHz=LowPASS
	fc=fcHz/fs
	#fc is between 0 and 0.5
	tau=math.exp(-2*math.pi*fc)
	a0=1-tau
	b1=tau
	print ('ts='+str(ts))
	print ('fs='+str(fs))
	print ('fc='+str(fc))
	print ('y[n]=a0*x[n]+b1*y[n-1]')
	print ('y[n]='+str(a0)+'*x[n]+'+str(b1)+'*y[n-1]')
	ynminus1 = 0
	for i,xn in enumerate(noise_LIST):
		#print str(xnminus1)
		yn=a0*xn+b1*ynminus1
		#print ('y['+str(i)+']='+str(yn))
		sig_LP.append(yn)
		ynminus1=yn
	noise_LIST=sig_LP

if highpass_ON:
#High pass filter random signals
# 'Single Pole IIR Filter'
	sig_HP=[]
	ts=TimeMULT
	fs=1/ts
	fcHz=HighPASS
	fc=fcHz/fs
	#fc is between 0 and 0.5
	tau=math.exp(-2*math.pi*fc)
	a0=(1+tau)/2
	a1=-a0	
	b1=tau
	
	print ('ts='+str(ts))
	print ('fs='+str(fs))
	print ('fc='+str(fc))
	print ('y[n]=a0*x[n]+a1*x[n-1]+b1*y[n-1]')
	print ('y[n]='+str(a0)+'*x[n]+'+str(a1)+'x[n-1]*+'+str(b1)+'*y[n-1]')
	xnminus1 = 0
	ynminus1 = 0
	for i,xn in enumerate(noise_LIST):
		#print str(xnminus1)
		yn=a0*xn+a1*xnminus1+b1*ynminus1
		#print ('y['+str(i)+']='+str(yn))
		sig_HP.append(yn)
		xnminus1=xn
		ynminus1=yn
	noise_LIST=sig_HP

#################
# Write to File #
#################
time_step=0
Sum_Squares=0
f = open('gaussian.pwl','w')
for t in noise_LIST:
	f.write(str(time_step*TimeMULT)+'\t'+str(t)+'\n')
	Sum_Squares += np.square(t)
	time_step+=1
f.close()

	
rms_volt=math.sqrt(Sum_Squares/(SimEND+1-SimSTART))
volt_per_root_Hz=rms_volt/math.sqrt(FREQ_HIGH-FREQ_LOW)
print "\n***************"
print "PWL Value Check"
print "***************"
print "\nvrms="+str(rms_volt)+" V"
print "\nvrms per root Hz="+str(volt_per_root_Hz)+" V/sqrt(Hz)"


##################
# Generate Plots #
##################
if graph_ON:
	print '\n%%%%%%%%%%%%'
	print   '% Plotting %'
	print   '%%%%%%%%%%%%'		
	#Plot Data
	plt.subplot(2,1,1)
	plt.title('Thermal Noise: T='+str(TEMP)+'K, R='+str(RES)+' ohms')
	#plt.text(2, STD_DEV*5, "numpy Plot", family="serif")
	#plt.text(2, STD_DEV*4, "T="+str(TEMP), family="serif")
	plt.xlabel('time['+str(TimeMULT)+'*s]')
	plt.ylabel('vrms')
	plt.ylim(-6*STD_DEV,6*STD_DEV)
	plt.plot(noise_LIST)
	
	#PSD Check
	plt.subplot(2,1,2)
	n = len(noise_LIST)
	timestep = TimeMULT#Inverse of Sampling Rate
	psd_VALS=plt.psd(noise_LIST, NFFT=len(noise_LIST), Fs=1/timestep, scale_by_freq=False)
	plt.xlim([0,1/(2*TimeMULT)])
	plt.show()
	psd_SUM_SQRS=0	
	j=0
	for i,ps in enumerate(psd_VALS[0]):
	        
		psd_SUM_SQRS = psd_SUM_SQRS + ps
		j= j+1
		
	print 'Average Power via Spectral Analysis:'
	print 'vrms per root Hz='+str(math.sqrt(psd_SUM_SQRS/(j)*1/((FREQ_HIGH-FREQ_LOW)*2)))+' V/sqrt(Hz)' ###????????
exit()