import operator import math import time import matplotlib.mlab as mlab import matplotlib.pyplot as plt import numpy as np class LFSR: def __init__(self, n, seed): self.tapsDict = {2 : (1,2),3: (1,3),4: (1,4), 5: (2,5),6: (1,6), 7: (3,7), 8: (2,3,4,8),9: (4,9),10: (3,10),11: (2,11),12: (1,4,6,12), 13: (1,3,4,13),14: (1,6,10,14),15: (1,15),16: (1,3,12,16),17: (3,17), 18: (7,18),19: (1,2,5,19),20: (3,20),21: (2,21),22: (1,22),23: (5,23), 24: (1,2,7,24),25: (3,25),26: (1,2,6,26),27: (1,2,5,27),28: (3,28), 29: (2,29),30: (1,2,23,30),31: (3,31),32: (1,2,22,32),33: (13,33), 34: (1,2,27,34)} self.n = n self.taps = self.tapsDict[self.n] self.shift_register = self.get_bits(seed, n) def shift(self): new_bit = sum([self.shift_register[tap-1] for tap in self.taps])%2 self.shift_register = self.shift_register[:-1] self.shift_register.insert(0, new_bit) return new_bit def random(self, length): bitstring="" for i in range(length): bitstring += str(self.shift()) return float(int(bitstring, 2))/pow(2,self.n) def get_bits(self, k, bits_to_get): """returns the first bits_to_get bits of k starting with the lowest""" bits = [] for i in xrange(bits_to_get): bits.append(k%2) k = k >> 1 return bits if __name__== "__main__": num_samples = 10000 my_LFSR = LFSR(16, int(1000 * time.time())) samples1 = [my_LFSR.random(16) for i in range(num_samples)] samples2 = [my_LFSR.random(16) for i in range(num_samples)] g_pairs = [[math.sqrt(-2*math.log(x1))*math.sin(2*math.pi*x2), math.sqrt(-2*math.log(x1))*math.cos(2*math.pi*x2)] for x1, x2 in zip(samples1, samples2)] g_samples = reduce(operator.__add__, g_pairs) plt.hist(g_samples) plt.show() uniform_m1 = sum(samples1)/len(samples1) uniform_m2 = sum([x**2 for x in samples1])/len(samples1) uniform_m3 = sum([x**3 for x in samples1])/len(samples1) gauss_m1 = sum(g_samples)/len(g_samples) gauss_m2 = sum([x**2 for x in g_samples])/len(g_samples) gauss_m3 = sum([x**3 for x in g_samples])/len(g_samples) uniform_C1 = uniform_m1 uniform_C2 = uniform_m2 - uniform_m1**2 uniform_C3 = uniform_m3 - 3*uniform_m1*uniform_m2 + 2*uniform_m1**3 gauss_C1 = gauss_m1 gauss_C2 = gauss_m2 - gauss_m1**2 gauss_C3 = gauss_m3 - 3*gauss_m1*gauss_m2 + 2*gauss_m1**3 print "Uniform cumulants: " print "C1: ", uniform_C1 print "C2: ", uniform_C2 print "C3: ", uniform_C3 print "Gaussian cumulants: " print "C1: ", gauss_C1 print "C2: ", gauss_C2 print "C3: ", gauss_C3