from numpy import *
from pylab import *

M = 24							# number of bits in shift register
bits = [0]*M					# initialize register with all zeros
bits[-1] = 1					# populate the last element of the register

# load lags data from file
lags = loadtxt("lag_values.txt", dtype = {'names': ('M', 'i'), 'formats': ('i2', 'S9')})

trow = lags[M-2]				# extract relevant row
data = trow[1]					# extract string of tap indices
taps = [0]*M					# initialize taps with all zeros

# set taps
start = 0
while True:
	end = data.find(",", start)						# find the next comma
	if end == -1:
		taps[int16(data[start:len(data)]) - 1] = 1	# extract the last tap index
		break
	else:
		taps[int16(data[start:end]) - 1] = 1		# extract tap index
		start = end + 1								# set the new start search value

def f(x): return 2**x								# initialize powers of 2
twos = map(f, range(M))

def y1(x1, x2): return sqrt(-2*log(x1))*sin(x2)		# define transformations
def y2(x1, x2): return sqrt(-2*log(x1))*cos(x2)		

all_x1 = []											# initialize output
all_x2 = []
all_y1 = []
all_y2 = []

N = 2000											# number of randoms generated
leap = M											# to reduce correlation of LFSR, leap forward M iterations

for i in range(N):
	for j in range(leap):
		bits.insert(0, mod(dot(bits, taps), 2))		# increment the register
		del bits[-1]								# push all bits to the right
	x1 = float(dot(bits, twos))/(2**M - 1)			# convert from binary to decimal, uniform distribution from 0 to 1
	all_x1.append(x1)								# save x1
	
	for j in range(leap):
		bits.insert(0, mod(dot(bits, taps), 2))		# increment the register
		del bits[-1]								# push all bits to the right
	x2 = float(dot(bits, twos))*2*pi/(2**M - 1)		# convert from binary to decimal, uniform distribution from 0 to 2*pi
	all_x2.append(x2)								# save x2
	
	all_y1.append(y1(x1, x2))						# compute and save y1, y2
	all_y2.append(y2(x1, x2))

# compute first, second, and third cumulants
C1_1 = float(sum(all_y1))/N
print 'The first cumulant of y1 is ' + str(C1_1)
C1_2 = float(sum(all_y2))/N
print 'The first cumulant of y2 is ' + str(C1_2)
C2_1 = float(sum([x**2 - C1_1**2 for x in all_y1]))/N
print 'The second cumulant of y1 is ' + str(C2_1)
C2_2 = float(sum([x**2 - C1_1**2 for x in all_y2]))/N
print 'The second cumulant of y2 is ' + str(C2_2)
C3_1 = float(sum([x**3 - 3*C2_1*C1_1 + 2*C1_1**3 for x in all_y1]))/N
print 'The third cumulant of y1 is ' + str(C3_1)
C3_2 = float(sum([x**3 - 3*C2_1*C1_1 + 2*C1_1**3 for x in all_y2]))/N
print 'The third cumulant of y2 is ' + str(C3_2)

# plot and save distribution
fig = figure()
ax = fig.add_subplot(111)
ax.scatter(all_y1, all_y2)
ax.set_aspect('equal')
title('LFSR with $M = 24$ mapped to Gaussian')
xlim([-3, 3])
ylim([-3, 3])
savefig('Gaussian_LFSR_M24_ff24.pdf')
show()