HG_REPOSITORIES/LPP/INSTRUMENTATION/SOLO_LFR/LFR_EMULATOR Commit - r0:90bacf12c042

Init

Alexis Jeandet -

r0:90bacf12c042 default

parent child

Context file:

r0:90bacf12c042

Collapse all files

SRC/basic_parameters.py created 644 +325 0

@@ -0,0 +1,325
	1	import numpy as np
	2	import math
	3	import cmath
	4
	5	import matplotlib.pyplot as plt
	6
	7	from bin16 import *
	8	from fft import *
	9
	10	import os
	11
	12	###########################################################
	13
	14	# (float)
	15	# it calculates the Spectral Matrice for the input signals, taking N points in consideration and an optional offset for the starting point
	16	def Spectral_Matrice(b1,b2,b3,e1,e2,fm,N,offset = 0):
	17
	18	#the input signals shall pass by a Hann window before being analysed with the Fourier Transform
	19	W = np.zeros(len(b1))
	20	W[offset:offset+N] = np.array([math.pow((math.sin(math.pi*i/(N-1))),2) for i in range(0,N)])
	21	b1W = b1 * W
	22	b2W = b2 * W
	23	b3W = b3 * W
	24	e1W = e1 * W
	25	e2W = e2 * W
	26
	27	B1 = np.fft.fft(b1W[offset:offset+N])/N
	28	B1 = B1[1:N/2+1]
	29	B2 = np.fft.fft(b2W[offset:offset+N])/N
	30	B2 = B2[1:N/2+1]
	31	B3 = np.fft.fft(b3W[offset:offset+N])/N
	32	B3 = B3[1:N/2+1]
	33	E1 = np.fft.fft(e1W[offset:offset+N])/N
	34	E1 = E1[1:N/2+1]
	35	E2 = np.fft.fft(e2W[offset:offset+N])/N
	36	E2 = E2[1:N/2+1]
	37
	38	#indices are SM[i][j][w] where i = line, j = colomn, w = frequency considered
	39	SM = [[B1B1.conjugate(), B1B2.conjugate(), B1B3.conjugate(), B1E1.conjugate(), B1*E2.conjugate()],
	40	[B2B1.conjugate(), B2B2.conjugate(), B2B3.conjugate(), B2E1.conjugate(), B2*E2.conjugate()],
	41	[B3B1.conjugate(), B3B2.conjugate(), B3B3.conjugate(), B3E1.conjugate(), B3*E2.conjugate()],
	42	[E1B1.conjugate(), E1B2.conjugate(), E1B3.conjugate(), E1E1.conjugate(), E1*E2.conjugate()],
	43	[E2B1.conjugate(), E2B2.conjugate(), E2B3.conjugate(), E2E1.conjugate(), E2*E2.conjugate()]]
	44
	45	w = 2math.pifm*np.arange(1,N/2+1)/N
	46
	47	return SM, w
	48
	49	# function that sums two SpectralMatrices
	50	def SumSpectralMatrices(M,N):
	51	Z = [[],[],[],[],[]]
	52	for i in range(5):
	53	for j in range(5):
	54	Z[i].append(M[i][j] + N[i][j])
	55	return Z
	56
	57	# (float)
	58	# function that takes the time average of several SpectralMatrices
	59	def Spectral_MatriceAvgTime(b1,b2,b3,e1,e2,fm,N):
	60	TBP = 4.0 # for the three cases of f in the LFR, TBP is equal to 4
	61	NSM = int(fm/N * TBP) # this gives the number of Spectral Matrices we should take into account for the average
	62	S,w = Spectral_Matrice(b1,b2,b3,e1,e2,fm,N,offset = 0)
	63	for k in range(1,NSM):
	64	Saux,w = Spectral_Matrice(b1,b2,b3,e1,e2,fm,N,offset = k*NSM)
	65	S = SumSpectralMatrices(S,Saux)
	66	for i in range(5):
	67	for j in range(5):
	68	S[i][j] = S[i][j]/NSM
	69	return S, w
	70
	71	# (float)
	72	# "Power spectrum density of the Magnetic Field"
	73	def PSD_B(S):
	74	PB = np.zeros(len(S[0][0]))
	75	PB = S[0][0] + S[1][1] + S[2][2]
	76	PB = abs(PB)
	77	return PB
	78
	79	# (float)
	80	# "Power spectrum density of the Electric Field"
	81	def PSD_E(S):
	82	PE = np.zeros(len(S[3][3]))
	83	PE = S[3][3] + S[4][4]
	84	PE = abs(PE)
	85	return PE
	86
	87	# (float)
	88	# "Ellipticity of the electromagnetic wave"
	89	def ellipticity(S):
	90	PB = PSD_B(S)
	91	e = np.zeros(len(PB))
	92	e = 2.0/PB * np.sqrt(S[0][1].imagS[0][1].imag + S[0][2].imagS[0][2].imag + S[1][2].imag*S[1][2].imag)
	93	return e
	94
	95	# (float)
	96	# "Degree of polarization of the electromagnetic wave"
	97	def degree_polarization(S):
	98	PB = PSD_B(S)
	99	d = np.zeros(len(PB))
	100	TrSCM2 = np.zeros(len(PB))
	101	TrSCM = np.zeros(len(PB))
	102	TrSCM2 = abs(S[0][0]S[0][0] + S[1][1]S[1][1] + S[2][2]S[2][2]) + 2 (abs(S[0][1]S[0][1]) + abs(S[0][2]S[0][2]) + abs(S[1][2]*S[1][2]))
	103	TrSCM = PB
	104	d = np.sqrt((3TrSCM2 - TrSCMTrSCM)/(2TrSCMTrSCM))
	105	return d
	106
	107	# (float)
	108	# "Normal wave vector"
	109	def normal_vector(S):
	110	n1 = np.zeros(len(S[0][0]))
	111	n2 = np.zeros(len(S[0][0]))
	112	n3 = np.zeros(len(S[0][0]))
	113	n1 = +1.0 * S[1][2].imag/np.sqrt(S[0][1].imagS[0][1].imag + S[0][2].imagS[0][2].imag + S[1][2].imag*S[1][2].imag)
	114	n2 = -1.0 * S[0][2].imag/np.sqrt(S[0][1].imagS[0][1].imag + S[0][2].imagS[0][2].imag + S[1][2].imag*S[1][2].imag)
	115	for i in range(len(n3)):
	116	n3[i] = math.copysign(1,S[0][1][i].imag)
	117	return [n1,n2,n3]
	118
	119	# (float)
	120	# this is a routine that takes a list (or vector) of 3 components for the n vector and gives back an example
	121	# of orthonormal matrix. This matrix may be used for the transformation matrix that takes every vector and puts it in the SCM referencial
	122	def Matrix_SCM(n):
	123
	124	n = np.array(n)
	125	n = n/np.linalg.norm(n)
	126
	127	v = np.array([np.random.random(),np.random.random(),np.random.random()])
	128
	129	v1 = v - np.dot(v,n)/np.dot(n,n) * n
	130	v1 = v1/np.linalg.norm(v1)
	131
	132	v2 = np.cross(n,v1)
	133	v2 = v2/np.linalg.norm(v2)
	134
	135	M = np.matrix([v1,v2,n]).transpose()
	136
	137	return M
	138
	139	# (float)
	140	# this is a routine that takes the values 'a' and 'b' and creates a B complex vector (phasor of magnetic field)
	141	# he then gives us the values of amplitude and phase for our cosines functions b1(t), b2(t), b3(t) for the SCM referencial (transformation using the R matrix)
	142	def MagneticComponents(a,b,w,R):
	143	B_PA = np.matrix([a,-b*1j,0]).transpose()
	144	B_SCM = R*B_PA
	145	return B_SCM
	146
	147	# (float)
	148	# this function takes some input values and gives back the corresponding Basic Parameters for the three possible fm's
	149	def SpectralMatrice_Monochromatic(a,b,n,E_para,f,fs,l):
	150	#'a' and 'b' are parameters related to the Magnetic Field
	151	#'n' is the normal vector to be considered
	152	#'E_para' is the degree of polarization of the propagating medium
	153	#'f' is the frequency of oscillation of the wave
	154	#'l' is the wavelength of the wave
	155	#'fm' is the original sampling frequency
	156
	157	time_step = 1.0/fs
	158	t_ini = 0
	159	t_fin = 8
	160	time_vec = np.arange(t_ini,t_fin,time_step)
	161
	162	# we must give the vector n, which has the information about the direction of k
	163	n = n/np.linalg.norm(n) #we ensure that n is an unitary vector
	164	k = 2.0 * math.pi / l
	165
	166	# the matrix R gives us an example of transformation to the SCM referencial
	167	# it was determined by the use of the normalized vector n and then two random vectors who form an orthogonal basis with it
	168	R = Matrix_SCM(n)
	169
	170	# now we transform the magnetic field to the SCM referencial
	171	B = MagneticComponents(a,b,2.0 * math.pi * f,R)
	172
	173	# now we have the caracteristics for our functions b1(t), b2(t) and b3(t) who shall enter in the LFR
	174	b1 = cmath.polar(B[0,0])[0] * np.cos(2.0 * math.pi * f * time_vec + cmath.polar(B[0,0])[1])
	175	b2 = cmath.polar(B[1,0])[0] * np.cos(2.0 * math.pi * f * time_vec + cmath.polar(B[1,0])[1])
	176	b3 = cmath.polar(B[2,0])[0] * np.cos(2.0 * math.pi * f * time_vec + cmath.polar(B[2,0])[1])
	177
	178	# the component of E who's parallel to n is of our choice. It represents the influence of the presence of charges where the electromagnetic
	179	# wave is travelling. However, the component of E who's perpendicular to n is determined by B and n as follows
	180
	181	E_orth = -1.0 * 2.0 * math.pi * f/k * np.cross(n,np.array(B.transpose()))
	182	E = E_orth + E_para
	183	E = E.transpose()
	184
	185	e1 = cmath.polar(E[0,0])[0] * np.cos(2.0 * math.pi * f * time_vec + cmath.polar(E[0,0])[1])
	186	e2 = cmath.polar(E[1,0])[0] * np.cos(2.0 * math.pi * f * time_vec + cmath.polar(E[1,0])[1])
	187
	188	N = 256
	189
	190	fm1 = 24576
	191	step1 = fs/fm1
	192	S1, w1 = Spectral_MatriceAvgTime(b1[::step1],b2[::step1],b3[::step1],e1[::step1],e2[::step1],fm1,N)
	193
	194	fm2 = 4096
	195	step2 = fs/fm2
	196	S2, w2 = Spectral_MatriceAvgTime(b1[::step2],b2[::step2],b3[::step2],e1[::step2],e2[::step2],fm2,N)
	197
	198	fm3 = 256
	199	step3 = fs/fm3
	200	S3, w3 = Spectral_MatriceAvgTime(b1[::step3],b2[::step3],b3[::step3],e1[::step3],e2[::step3],fm3,N)
	201
	202	question = raw_input("Do you wish to save the plots (Y/N): ")
	203	if question == 'Y':
	204	BasicParameters_plot(S1,w1,fm1,"fm_%s" % fm1)
	205	BasicParameters_plot(S2,w2,fm2,"fm_%s" % fm2)
	206	BasicParameters_plot(S3,w3,fm3,"fm_%s" % fm3)
	207
	208	return [S1,S2,S3], [w1,w2,w3]
	209
	210	# (float)
	211	# this function takes a Spectral Matrice in the input and gives a list with all the associated Basic Parameters
	212	def BasicParameters(S,w):
	213	PB = PSD_B(S)
	214	PE = PSD_E(S)
	215	d = degree_polarization(S)
	216	e = ellipticity(S)
	217	n = normal_vector(S)
	218	return [PB,PE,d,e,n]
	219
	220	# (float)
	221	# this function saves plots in .pdf for each of the Basic Parameters associated with the Spectral Matrice in the input
	222	def BasicParameters_plot(S,w,fm,mypath):
	223
	224	if not os.path.isdir(mypath):
	225	os.makedirs(mypath)
	226
	227	PB = PSD_B(S)
	228	plt.plot(w/(2*math.pi),PB)
	229	plt.title("Power spectrum density of the Magnetic Field (fm = %s Hz)" % fm)
	230	plt.xlabel("frequency [Hz]")
	231	plt.ylabel("PB(w)")
	232	plt.savefig('%s/PB.png' % mypath)
	233	plt.close()
	234
	235	PE = PSD_E(S)
	236	plt.plot(w/(2*math.pi),PE)
	237	plt.title("Power spectrum density of the Electric Field (fm = %s Hz)" % fm)
	238	plt.xlabel("frequency [Hz]")
	239	plt.ylabel("PE(w)")
	240	plt.savefig('%s/PE.png' % mypath)
	241	plt.close()
	242
	243	e = ellipticity(S)
	244	plt.plot(w/(2*math.pi),e)
	245	plt.title("Ellipticity of the electromagnetic wave (fm = %s Hz)" % fm)
	246	plt.xlabel("frequency [Hz]")
	247	plt.ylabel("e(w)")
	248	plt.savefig('%s/e.png' % mypath)
	249	plt.close()
	250
	251	d = degree_polarization(S)
	252	plt.plot(w/(2*math.pi),d)
	253	plt.title("Degree of polarization of the electromagnetic wave (fm = %s Hz)" % fm)
	254	plt.xlabel("frequency [Hz]")
	255	plt.ylabel("d(w)")
	256	plt.savefig('%s/d.png' % mypath)
	257	plt.close()
	258
	259	[n1, n2, n3] = normal_vector(S)
	260	print "n1: "
	261	print n1
	262	print "n2: "
	263	print n2
	264	print "n3: "
	265	print n3
	266
	267	# (float)
	268	# this function takes some input values and gives back the corresponding Basic Parameters for the three possible fm's
	269	def SignalsWaveMonochromatic(a,b,n,E_para,f,fs,l):
	270	#'a' and 'b' are parameters related to the Magnetic Field
	271	#'n' is the normal vector to be considered
	272	#'E_para' is the degree of polarization of the propagating medium
	273	#'f' is the frequency of oscillation of the wave
	274	#'l' is the wavelength of the wave
	275	#'fm' is the original sampling frequency
	276
	277	time_step = 1.0/fs
	278	t_ini = 0
	279	t_fin = 8
	280	time_vec = np.arange(t_ini,t_fin,time_step)
	281
	282	# we must give the vector n, which has the information about the direction of k
	283	n = n/np.linalg.norm(n) #we ensure that n is an unitary vector
	284	k = 2.0 * math.pi / l
	285
	286	print " 'a' = ", a
	287	print " 'b' = ", b
	288	print " n = ", n
	289	print " l = ", l
	290	print " f = ", f
	291	print "E_para = ", E_para
	292	print " "
	293
	294	# the matrix R gives us an example of transformation to the SCM referencial
	295	# it was determined by the use of the normalized vector n and then two random vectors who form an orthogonal basis with it
	296	R = Matrix_SCM(n)
	297
	298	# now we transform the magnetic field to the SCM referencial
	299	B = MagneticComponents(a,b,2.0 * math.pi * f,R)
	300
	301	# now we have the caracteristics for our functions b1(t), b2(t) and b3(t) who shall enter in the LFR
	302	b1 = cmath.polar(B[0,0])[0] * np.cos(2.0 * math.pi * f * time_vec + cmath.polar(B[0,0])[1])
	303	b2 = cmath.polar(B[1,0])[0] * np.cos(2.0 * math.pi * f * time_vec + cmath.polar(B[1,0])[1])
	304	b3 = cmath.polar(B[2,0])[0] * np.cos(2.0 * math.pi * f * time_vec + cmath.polar(B[2,0])[1])
	305
	306	print "b1(t) = ", cmath.polar(B[0,0])[0], "cos(w*t + ", cmath.polar(B[0,0])[1], ")"
	307	print "b2(t) = ", cmath.polar(B[1,0])[0], "cos(w*t + ", cmath.polar(B[1,0])[1], ")"
	308	print "b3(t) = ", cmath.polar(B[2,0])[0], "cos(w*t + ", cmath.polar(B[2,0])[1], ")"
	309
	310	# the component of E who's parallel to n is of our choice. It represents the influence of the presence of charges where the electromagnetic
	311	# wave is travelling. However, the component of E who's perpendicular to n is determined by B and n as follows
	312
	313	E_orth = -1.0 * 2.0 * math.pi * f/k * np.cross(n,np.array(B.transpose()))
	314	E = E_orth + E_para
	315	E = E.transpose()
	316
	317	e1 = cmath.polar(E[0,0])[0] * np.cos(2.0 * math.pi * f * time_vec + cmath.polar(E[0,0])[1])
	318	e2 = cmath.polar(E[1,0])[0] * np.cos(2.0 * math.pi * f * time_vec + cmath.polar(E[1,0])[1])
	319
	320	print "e1(t) = ", cmath.polar(E[0,0])[0], "cos(w*t + ", cmath.polar(E[0,0])[1], ")"
	321	print "e2(t) = ", cmath.polar(E[1,0])[0], "cos(w*t + ", cmath.polar(E[1,0])[1], ")"
	322	print "(the value of w is ", 2.0 * math.pi * f, ")"
	323	print " "
	324
	325	return [b1,b2,b3,e1,e2] No newline at end of file

SRC/basic_parameters_Int.py created 644 +296 0

@@ -0,0 +1,296
	1	import numpy as np
	2	import math
	3	import cmath
	4
	5	import matplotlib.pyplot as plt
	6
	7	from bin16 import *
	8	from fft import *
	9
	10	import os
	11
	12	# (integers)
	13	# it calculates the Spectral Matrice for the input signals, taking N points in consideration and an optional offset for the starting point
	14	def Spectral_Matrice_Int(b1,b2,b3,e1,e2,fm,N,offset = 0):
	15
	16	#the input signals shall pass by a Hann window before being analysed with the Fourier Transform
	17	W = np.zeros(len(b1))
	18	W[offset:offset+N] = np.array([math.pow((math.sin(math.pi*i/(N-1))),2) for i in range(0,N)])
	19	b1W = b1 * W
	20	b2W = b2 * W
	21	b3W = b3 * W
	22	e1W = e1 * W
	23	e2W = e2 * W
	24
	25	#this is for keeping stuff as integers
	26	for i in range(len(b1)):
	27	b1W[i] = int(b1W[i])
	28	b2W[i] = int(b2W[i])
	29	b3W[i] = int(b3W[i])
	30	e1W[i] = int(e1W[i])
	31	e2W[i] = int(e2W[i])
	32
	33	#remembering that fft_CT already divides the FFT by N
	34	B1 = fft_CT(b1W[offset:offset+N])
	35	B1 = B1[1:N/2+1]
	36	B2 = fft_CT(b2W[offset:offset+N])
	37	B2 = B2[1:N/2+1]
	38	B3 = fft_CT(b3W[offset:offset+N])
	39	B3 = B3[1:N/2+1]
	40	E1 = fft_CT(e1W[offset:offset+N])
	41	E1 = E1[1:N/2+1]
	42	E2 = fft_CT(e2W[offset:offset+N])
	43	E2 = E2[1:N/2+1]
	44
	45	#indices are SM[i][j][w] where i = line, j = colomn, w = frequency considered
	46	SM = [[B1B1.conjugate(), B1B2.conjugate(), B1B3.conjugate(), B1E1.conjugate(), B1*E2.conjugate()],
	47	[B2B1.conjugate(), B2B2.conjugate(), B2B3.conjugate(), B2E1.conjugate(), B2*E2.conjugate()],
	48	[B3B1.conjugate(), B3B2.conjugate(), B3B3.conjugate(), B3E1.conjugate(), B3*E2.conjugate()],
	49	[E1B1.conjugate(), E1B2.conjugate(), E1B3.conjugate(), E1E1.conjugate(), E1*E2.conjugate()],
	50	[E2B1.conjugate(), E2B2.conjugate(), E2B3.conjugate(), E2E1.conjugate(), E2*E2.conjugate()]]
	51
	52	#this we can keep as float
	53	w = 2math.pifm*np.arange(1,N/2+1)/N
	54
	55	return SM, w
	56
	57	# function that sums two SpectralMatrices
	58	def SumSpectralMatrices(M,N):
	59	Z = [[],[],[],[],[]]
	60	for i in range(5):
	61	for j in range(5):
	62	Z[i].append(M[i][j] + N[i][j])
	63	return Z
	64
	65	# (integers)
	66	# function that takes the time average of several SpectralMatrices
	67	def Spectral_MatriceAvgTime_Int(b1,b2,b3,e1,e2,fm,N):
	68	TBP = 4.0 # for the three cases of f in the LFR, TBP is equal to 4
	69	NSM = int(fm/N * TBP) # this gives the number of Spectral Matrices we should take into account for the average
	70	S,w = Spectral_Matrice_Int(b1,b2,b3,e1,e2,fm,N,offset = 0)
	71	for k in range(1,NSM):
	72	Saux,w = Spectral_Matrice_Int(b1,b2,b3,e1,e2,fm,N,offset = k*NSM)
	73	S = SumSpectralMatrices(S,Saux)
	74	for i in range(5):
	75	for j in range(5):
	76	for k in range(len(S[i][j])):
	77	S[i][j][k] = complex(int(S[i][j][k].real/NSM),int(S[i][j][k].imag/NSM))
	78	return S, w
	79
	80	# (integers)
	81	# "Power spectrum density of the Magnetic Field"
	82	# it's being coded over 32 bits for now
	83	def PSD_B_Int(S):
	84	PB = np.zeros(len(S[0][0]))
	85	PB = S[0][0] + S[1][1] + S[2][2]
	86	PB = abs(PB)
	87	return PB
	88
	89	# (integers)
	90	# "Power spectrum density of the Electric Field"
	91	# it's being coded over 32 bits for now
	92	def PSD_E_Int(S):
	93	PE = np.zeros(len(S[3][3]))
	94	PE = S[3][3] + S[4][4]
	95	PE = abs(PE)
	96	return PE
	97
	98	# (integers)
	99	# "Ellipticity of the electromagnetic wave"
	100	# Coding over 4 bits, from 0 to 1
	101	def ellipticity_Int(S):
	102	PB = PSD_B_Int(S)
	103	e = np.zeros(len(PB))
	104	e = 2.0/PB * np.sqrt(1.0(S[0][1].imagS[0][1].imag + S[0][2].imagS[0][2].imag + S[1][2].imagS[1][2].imag))
	105	for i in range(len(e)):
	106	if math.isnan(e[i]):
	107	continue
	108	e[i] = quantN(e[i],0,1,4)
	109	e[i] = truncN(e[i],4)
	110	return e
	111
	112	# (integers)
	113	# "Degree of polarization of the electromagnetic wave"
	114	# Coding over 3 bits, from 0 to 1
	115	def degree_polarization_Int(S):
	116	PB = PSD_B_Int(S)
	117	d = np.zeros(len(PB))
	118	TrSCM2 = np.zeros(len(PB))
	119	TrSCM = np.zeros(len(PB))
	120	TrSCM2 = abs(S[0][0]S[0][0] + S[1][1]S[1][1] + S[2][2]S[2][2]) + 2 (abs(S[0][1]S[0][1]) + abs(S[0][2]S[0][2]) + abs(S[1][2]*S[1][2]))
	121	TrSCM = PB
	122	d = np.sqrt((3.0TrSCM2 - 1.0TrSCMTrSCM)/(2.0TrSCM*TrSCM))
	123	for i in range(len(d)):
	124	if not(math.isnan(d[i])):
	125	d[i] = quantN(d[i],0,1,3)
	126	d[i] = truncN(d[i],3)
	127	return d
	128
	129	# (integers)
	130	# "Normal wave vector"
	131	# Coding over 8 bits for n1 and n2, 1 bit for the n3 (0 if n3 = +1, 1 if n3 = -1)
	132	def normal_vector_Int(S):
	133	n1 = np.zeros(len(S[0][0]))
	134	n2 = np.zeros(len(S[0][0]))
	135	n3 = np.zeros(len(S[0][0]))
	136	n1 = +1.0 * S[1][2].imag/np.sqrt(S[0][1].imagS[0][1].imag + S[0][2].imagS[0][2].imag + S[1][2].imag*S[1][2].imag)
	137	n2 = -1.0 * S[0][2].imag/np.sqrt(S[0][1].imagS[0][1].imag + S[0][2].imagS[0][2].imag + S[1][2].imag*S[1][2].imag)
	138	for i in range(len(n3)):
	139	n3[i] = math.copysign(1,S[0][1][i].imag)
	140	for i in range(len(n1)):
	141	if not(math.isnan(n1[i])):
	142	n1[i] = quantN(n1[i],0,1,8)
	143	n1[i] = truncN(n1[i],8)
	144	if not(math.isnan(n2[i])):
	145	n2[i] = quantN(n2[i],0,1,8)
	146	n2[i] = truncN(n2[i],8)
	147	if not(math.isnan(n3[i])):
	148	if n3[i] == -1:
	149	n3[i] = 1
	150	if n3[i] == +1:
	151	n3[i] = 0
	152	return [n1,n2,n3]
	153
	154	# (integers)
	155	# "Z-component of the normalized Poynting flux"
	156	# Coding over 8 bits
	157	def poynting_vector_Int(S):
	158	return poynting
	159
	160	# (integers)
	161	# "Phase velocity estimator"
	162	# Coding over 8 bits
	163	def phase_velocity_Int(S):
	164	return vp
	165
	166	# (integers)
	167	# "Autocorrelation"
	168	# it's being coded over 32 bits for now
	169	def autocorrelation(S):
	170	return [S[0][0],S[1][1],S[2][2],S[3][3],S[4][4]]
	171
	172	# (integers)
	173	# Normalized cross correlations
	174	# Coding over 8 bits, for values between -1 and +1
	175	def cross_correlations_Int(S):
	176
	177	R12 = 1.0 * S[0][1].real/np.sqrt(1.0 * S[0][0] * S[1][1])
	178	R12 = quantN(R12,-1,+1,8)
	179	R12 = truncN(R12,8)
	180	R13 = 1.0 * S[0][2].real/np.sqrt(1.0 * S[0][0] * S[2][2])
	181	R13 = quantN(R13,-1,+1,8)
	182	R13 = truncN(R13,8)
	183	R23 = 1.0 * S[1][2].real/np.sqrt(1.0 * S[1][1] * S[2][2])
	184	R23 = quantN(R23,-1,+1,8)
	185	R23 = truncN(R23,8)
	186	R45 = 1.0 * S[3][4].real/np.sqrt(1.0 * S[3][3] * S[4][4])
	187	R45 = quantN(R45,-1,+1,8)
	188	R45 = truncN(R45,8)
	189	R14 = 1.0 * S[0][3].real/np.sqrt(1.0 * S[0][0] * S[3][3])
	190	R14 = quantN(R14,-1,+1,8)
	191	R14 = truncN(R14,8)
	192	R15 = 1.0 * S[0][4].real/np.sqrt(1.0 * S[0][0] * S[4][4])
	193	R15 = quantN(R15,-1,+1,8)
	194	R15 = truncN(R15,8)
	195	R24 = 1.0 * S[1][3].real/np.sqrt(1.0 * S[1][1] * S[3][3])
	196	R24 = quantN(R24,-1,+1,8)
	197	R24 = truncN(R24,8)
	198	R25 = 1.0 * S[1][4].real/np.sqrt(1.0 * S[1][1] * S[4][4])
	199	R25 = quantN(R25,-1,+1,8)
	200	R25 = truncN(R25,8)
	201	R34 = 1.0 * S[2][3].real/np.sqrt(1.0 * S[2][2] * S[3][3])
	202	R34 = quantN(R34,-1,+1,8)
	203	R34 = truncN(R34,8)
	204	R35 = 1.0 * S[2][4].real/np.sqrt(1.0 * S[2][2] * S[4][4])
	205	R35 = quantN(R35,-1,+1,8)
	206	R35 = truncN(R35,8)
	207
	208	I12 = 1.0 * S[0][1].imag/np.sqrt(1.0 * S[0][0] * S[1][1])
	209	I12 = quantN(I12,-1,+1,8)
	210	I12 = truncN(I12,8)
	211	I13 = 1.0 * S[0][2].imag/np.sqrt(1.0 * S[0][0] * S[2][2])
	212	I13 = quantN(I13,-1,+1,8)
	213	I13 = truncN(I13,8)
	214	I23 = 1.0 * S[1][2].imag/np.sqrt(1.0 * S[1][1] * S[2][2])
	215	I23 = quantN(I23,-1,+1,8)
	216	I23 = truncN(I23,8)
	217	I45 = 1.0 * S[3][4].imag/np.sqrt(1.0 * S[3][3] * S[4][4])
	218	I45 = quantN(I45,-1,+1,8)
	219	I45 = truncN(I45,8)
	220	I14 = 1.0 * S[0][3].imag/np.sqrt(1.0 * S[0][0] * S[3][3])
	221	I14 = quantN(I14,-1,+1,8)
	222	I14 = truncN(I14,8)
	223	I15 = 1.0 * S[0][4].imag/np.sqrt(1.0 * S[0][0] * S[4][4])
	224	I15 = quantN(I15,-1,+1,8)
	225	I15 = truncN(I15,8)
	226	I24 = 1.0 * S[1][3].imag/np.sqrt(1.0 * S[1][1] * S[3][3])
	227	I24 = quantN(I24,-1,+1,8)
	228	I24 = truncN(I24,8)
	229	I25 = 1.0 * S[1][4].imag/np.sqrt(1.0 * S[1][1] * S[4][4])
	230	I25 = quantN(I25,-1,+1,8)
	231	I25 = truncN(I25,8)
	232	I34 = 1.0 * S[2][3].imag/np.sqrt(1.0 * S[2][2] * S[3][3])
	233	I34 = quantN(I34,-1,+1,8)
	234	I34 = truncN(I34,8)
	235	I35 = 1.0 * S[2][4].imag/np.sqrt(1.0 * S[2][2] * S[4][4])
	236	I35 = quantN(I35,-1,+1,8)
	237	I35 = truncN(I35,8)
	238
	239	return [R12,R13,R23,R45,R14,R15,R24,R25,R34,R35],[I12,I13,I23,I45,I14,I15,I24,I25,I34,I35]
	240
	241	# (integers)
	242	# this function takes a Spectral Matrice in the input and gives a list with all the associated Basic Parameters
	243	def BasicParameters_Int(S,w):
	244	PB = PSD_B_Int(S)
	245	PE = PSD_E_Int(S)
	246	d = degree_polarization_Int(S)
	247	e = ellipticity_Int(S)
	248	n = normal_vector_Int(S)
	249	return [PB,PE,d,e,n]
	250
	251	# (integers)
	252	# this function saves plots in .pdf for each of the Basic Parameters associated with the Spectral Matrice in the input
	253	def BasicParameters_plot_Int(S,w,fm,mypath):
	254
	255	if not os.path.isdir(mypath):
	256	os.makedirs(mypath)
	257
	258	PB = PSD_B_Int(S)
	259	plt.plot(w/(2*math.pi),PB)
	260	plt.title("Power spectrum density of the Magnetic Field (fm = %s Hz)" % fm)
	261	plt.xlabel("frequency [Hz]")
	262	plt.ylabel("PB(w)")
	263	plt.savefig('%s/PB.png' % mypath)
	264	plt.close()
	265
	266	PE = PSD_E_Int(S)
	267	plt.plot(w/(2*math.pi),PE)
	268	plt.title("Power spectrum density of the Electric Field (fm = %s Hz)" % fm)
	269	plt.xlabel("frequency [Hz]")
	270	plt.ylabel("PE(w)")
	271	plt.savefig('%s/PE.png' % mypath)
	272	plt.close()
	273
	274	e = ellipticity_Int(S)
	275	plt.plot(w/(2*math.pi),e)
	276	plt.title("Ellipticity of the electromagnetic wave (fm = %s Hz)" % fm)
	277	plt.xlabel("frequency [Hz]")
	278	plt.ylabel("e(w)")
	279	plt.savefig('%s/e.png' % mypath)
	280	plt.close()
	281
	282	d = degree_polarization_Int(S)
	283	plt.plot(w/(2*math.pi),d)
	284	plt.title("Degree of polarization of the electromagnetic wave (fm = %s Hz)" % fm)
	285	plt.xlabel("frequency [Hz]")
	286	plt.ylabel("d(w)")
	287	plt.savefig('%s/d.png' % mypath)
	288	plt.close()
	289
	290	[n1,n2,n3] = normal_vector_Int(S)
	291	print "n1"
	292	print n1
	293	print "n2"
	294	print n2
	295	print "n3"
	296	print n3 No newline at end of file

SRC/bin16.py created 644 +59 0

@@ -0,0 +1,59
	1	import math
	2
	3	def quant16(x,liminf,limsup):
	4	p = 1.0/pow(2,16) * (limsup - liminf)
	5	x_quant = int(round(x/p))
	6	return x_quant
	7
	8	def quantN(x,liminf,limsup,N):
	9	p = 1.0/pow(2,N) * (limsup - liminf)
	10	x_quant = int(round(x/p))
	11	return x_quant
	12
	13	def trunc16(x):
	14	if x > (pow(2,15) - 1):
	15	z = pow(2,15) - 1
	16	z = int(z)
	17	return z
	18	if x < (-1 * pow(2,15)):
	19	z = -1 * pow(2,15)
	20	z = int(z)
	21	return z
	22	return x
	23
	24	def truncN(x,N,signed = 0):
	25	if signed == 0:
	26	if x > (pow(2,N)):
	27	z = pow(2,N) - 1
	28	z = int(z)
	29	return z
	30	if x < 0:
	31	z = 0
	32	return z
	33	return x
	34	if signed == 1:
	35	if x > (pow(2,N-1) - 1):
	36	z = pow(2,N-1) - 1
	37	z = int(z)
	38	return z
	39	if x < (-1 * pow(2,N-1)):
	40	z = -1 * pow(2,N-1)
	41	z = int(z)
	42	return z
	43	return x
	44
	45	def add16(x,y):
	46	z = x + y
	47	z = trunc16(z)
	48	return z
	49
	50	def prod16(x,y):
	51	z = x * y
	52	z = trunc16(z)
	53	return z
	54
	55	def div16(x,y):
	56	z = x/y
	57	z = trunc16(z)
	58	return z
	59

SRC/fft.py created 644 +52 0

@@ -0,0 +1,52
	1	import math
	2	import numpy as np
	3	from bin16 import *
	4
	5	def power_of_2(number):
	6	if number == 0:
	7	return 0
	8	while(number % 2 == 0):
	9	number = number / 2
	10	if number > 1:
	11	return 0
	12	return 1
	13
	14	#here we simply do the FFT algorithm, using the Cooley-Tukey idea of divide and conquer and quantizing the values (working with integers)
	15	#note that the output result will be already divided by N, which is not what happens for the np.fft routine
	16	def fft_CT(x):
	17
	18	N = len(x)
	19
	20	if N == 1:
	21	return x
	22
	23	x_e = x[::2]
	24	x_o = x[1::2]
	25	X_e = fft_CT(x_e)
	26	X_o = fft_CT(x_o)
	27
	28	X = np.zeros(N,'complex')
	29
	30	M = N/2
	31	for k in range(M):
	32	X[k] = X_e[k] + X_o[k] * pow(math.e,-2jmath.pik/N)
	33	a = X[k].real/2
	34	b = X[k].imag/2
	35	a = int(a)
	36	a = trunc16(a)
	37	b = int(b)
	38	b = trunc16(b)
	39	c = complex(a,b)
	40	X[k] = c
	41	for k in range(M,N):
	42	X[k] = X_e[k-M] - X_o[k-M] * pow(math.e,-2jmath.pi(k-M)/N)
	43	a = X[k].real/2
	44	b = X[k].imag/2
	45	a = int(a)
	46	a = trunc16(a)
	47	b = int(b)
	48	b = trunc16(b)
	49	c = complex(a,b)
	50	X[k] = c
	51	return X
	52	No newline at end of file

SRC/fft_integerVsfloat.py created 644 +50 0

@@ -0,0 +1,50
	1	import math
	2	import numpy as np
	3	from bin16 import *
	4	from fft import *
	5	import matplotlib.pyplot as plt
	6
	7	alfa = 8192.0/8000.0
	8
	9	t_i = 0.0*alfa
	10	t_f = 8.0*alfa
	11	t_a = 3.5*alfa
	12	t_b = 4.5*alfa
	13	time_step = 1.0/1000
	14
	15	time_vec = np.arange(t_i,t_f,time_step)
	16
	17	Afloat = np.zeros(len(time_vec))
	18	for i in range(int(t_a/time_step),int(t_b/time_step)):
	19	Afloat[i] = 1.0
	20
	21	# plt.figure(1)
	22	# plt.plot(time_vec,Afloat,'r')
	23	# plt.ylim(0,2)
	24	# plt.show()
	25
	26	Afloat_FFT = np.fft.fft(Afloat)/len(time_vec)
	27
	28	# plt.figure(2)
	29	# plt.plot(abs(Afloat_FFT[0:100]),'g')
	30	# plt.show()
	31
	32	Aint = np.zeros(len(time_vec))
	33	for i in range(len(time_vec)):
	34	Aint[i] = quant16(Afloat[i],-2,+2)
	35
	36	# plt.figure(3)
	37	# plt.plot(time_vec,Aint,'r')
	38	# plt.show()
	39
	40	Aint_FFT = fft_CT(Aint)
	41
	42	plt.figure(4)
	43	plt.plot(abs(Afloat_FFT[0:100]),'g')
	44	plt.scatter(range(100),1.0 * abs(Aint_FFT[0:100])/16384)
	45	plt.title("Comparing the FFT's")
	46	plt.show()
	47
	48
	49
	50

SRC/filters.py created 644 +91 0

@@ -0,0 +1,91
	1	# -- coding: utf-8 --
	2	"""
	3	Created on Tue Nov 13 15:53:23 2012
	4
	5	@author: pedroluizcoelhorodrigues
	6	"""
	7	import numpy as np
	8	from bin16 import *
	9
	10	class Filter_Coefficients():
	11	"""List of coefficients of a IIR filter of order 2"""
	12	def __init__(self,a0=0,a1=0,a2=0,b0=0,b1=0,b2=0):
	13	self.b0 = b0
	14	self.b1 = b1
	15	self.b2 = b2
	16	self.a0 = a0
	17	self.a1 = a1
	18	self.a2 = a2
	19
	20	class Filter(object):
	21	"""Basic Filter of Order 2"""
	22
	23	def __init__(self,Coefficients,x1=0,x2=0,y1=0,y2=0):
	24	self.Coefficients = Coefficients
	25	self.x1 = x1
	26	self.x2 = x2
	27	self.y1 = y1
	28	self.y2 = y2
	29
	30	def print_coefficients(self):
	31	print 'a0 = ', self.Coefficients.a0
	32	print 'a1 = ', self.Coefficients.a1
	33	print 'a2 = ', self.Coefficients.a2
	34	print 'b0 = ', self.Coefficients.b0
	35	print 'b1 = ', self.Coefficients.b1
	36	print 'b2 = ', self.Coefficients.b2, "\n"
	37
	38	def convolve(self,input):
	39
	40	B = np.array([self.Coefficients.b0,self.Coefficients.b1,self.Coefficients.b2])
	41	A = np.array([self.Coefficients.a1,self.Coefficients.a2])
	42
	43	y = np.zeros(len(input))
	44	y[0] = round(1.0/(self.Coefficients.a0) * (np.dot(B,np.array([input[0],self.x1,self.x2])) - np.dot(A,np.array([self.y1,self.y2]))))
	45	y[0] = trunc16(y[0])
	46	y[1] = round(1.0/(self.Coefficients.a0) * (np.dot(B,np.array([input[1],input[0],self.x1])) - np.dot(A,np.array([y[0],self.y1]))))
	47	y[1] = trunc16(y[1])
	48
	49	for j in range(2,len(input)):
	50	y[j] = trunc16(round(1.0/(self.Coefficients.a0) * (np.dot(B,np.array([input[j],input[j-1],input[j-2]])) - np.dot(A,np.array([y[j-1],y[j-2]])))))
	51
	52	self.x1 = input[j]
	53	self.x2 = input[j-1]
	54	self.y1 = y[j]
	55	self.y2 = y[j-1]
	56
	57	return y
	58
	59	def reset(self):
	60	self.x1 = 0
	61	self.x2 = 0
	62	self.y1 = 0
	63	self.y2 = 0
	64
	65	def impulsive_response(self, N):
	66	delta = np.zeros(N)
	67	delta[0] = 1
	68	return self.convolve(delta)
	69
	70	class FilterGlobal(object):
	71	"""Global filter with other basic filters (of order 2) in cascade"""
	72	def __init__(self,Filters):
	73	self.Filter = Filters
	74
	75	def convolve(self,x):
	76	y = np.zeros(len(x))
	77	y = self.Filter[0].convolve(x)
	78	if len(self.Filter) > 1:
	79	for j in range(1,len(self.Filter)):
	80	y = self.Filter[j].convolve(y)
	81	return y
	82
	83	def impulsive_response(self, N):
	84	delta = np.zeros(N)
	85	delta[0] = 1
	86	return self.convolve(delta)
	87
	88	def reset(self):
	89	for j in range(len(self.Filter)):
	90	self.Filter[j].reset()
	91	No newline at end of file

SRC/processing_chain.py created 644 +178 0

@@ -0,0 +1,178
	1	#the signal processing chain in the LFR analyser.
	2
	3	import numpy as np
	4	from filters import *
	5	from bin16 import *
	6	from fft import *
	7	import math
	8	import cmath
	9	import matplotlib.pyplot as plt
	10	from basic_parameters import *
	11	import os
	12
	13	# The parameters who describe the magnetic field B are 'a' and 'b',
	14	# who shall form a complex phasor given by: [a -b1j 0] exp(1j * (w*t + phi))
	15	a = np.random.random() # we take a value at random between 0 and 1
	16	b = np.sqrt(1 - a*a) # we want to force a^2 + b^2 = 1
	17
	18	# 'f' is the frequency of oscillation of our electromagnetic wave, who's monochromatic for now
	19	f = 2500
	20	# 'l' is the wavelength of the electromagnetic wave
	21	l = 3000
	22	k = 2.0 * math.pi / l
	23	# is the propagating vector, who's related to the k vector
	24	n = [1,1,1]
	25	# is a vector who tells us the degree of polarization of the medium where the electromagnetic wave is being propagated
	26	E_para = [0,0,0]
	27	# 'fs' is the original sampling frequency (the one who enters the LFR)
	28	fs = 98304
	29	t_ini = 0
	30	t_fin = 10
	31	time_step = 1.0/fs
	32	time_vec = np.arange(t_ini, t_fin, time_step)
	33
	34	# the matrix R gives us an example of transformation to the SCM referencial
	35	# it was determined by the use of the normalized vector n and then two random vectors who form an orthogonal basis with it
	36	R = Matrix_SCM(n)
	37
	38	# now we transform the magnetic field to the SCM referencial
	39	B = MagneticComponents(a,b,2.0 * math.pi * f,R)
	40
	41	# now we have the caracteristics for our functions b1(t), b2(t) and b3(t) who shall enter in the LFR
	42	b1 = cmath.polar(B[0,0])[0] * np.cos(2.0 * math.pi * f * time_vec + cmath.polar(B[0,0])[1])
	43	b2 = cmath.polar(B[1,0])[0] * np.cos(2.0 * math.pi * f * time_vec + cmath.polar(B[1,0])[1])
	44	b3 = cmath.polar(B[2,0])[0] * np.cos(2.0 * math.pi * f * time_vec + cmath.polar(B[2,0])[1])
	45
	46	# the component of E who's parallel to n is of our choice. It represents the influence of the presence of charges where the electromagnetic
	47	# wave is travelling. However, the component of E who's perpendicular to n is determined by B and n as follows
	48
	49	E_orth = -1.0 * 2.0 * math.pi * f/k * np.cross(n,np.array(B.transpose()))
	50	E = E_orth + E_para
	51	E = E.transpose()
	52
	53	e1 = cmath.polar(E[0,0])[0] * np.cos(2.0 * math.pi * f * time_vec + cmath.polar(E[0,0])[1])
	54	e2 = cmath.polar(E[1,0])[0] * np.cos(2.0 * math.pi * f * time_vec + cmath.polar(E[1,0])[1])
	55
	56	liminfB = -1
	57	limsupB = +1
	58	liminfE = -10000000
	59	limsupE = +10000000
	60
	61	b1 = b1 + np.random.uniform(liminfB/100,limsupB/100,len(time_vec))
	62	b2 = b2 + np.random.uniform(liminfB/100,limsupB/100,len(time_vec))
	63	b3 = b3 + np.random.uniform(liminfB/100,limsupB/100,len(time_vec))
	64	e1 = e1 + np.random.uniform(liminfE/100,limsupE/100,len(time_vec))
	65	e2 = e2 + np.random.uniform(liminfE/100,limsupE/100,len(time_vec))
	66
	67	plt.figure(1)
	68	plt.plot(b1[0:len(time_vec)/5000],'r')
	69	plt.plot(b2[0:len(time_vec)/5000],'g')
	70	plt.plot(b3[0:len(time_vec)/5000],'b')
	71	plt.show()
	72
	73	#we want the input signal to be quantized (integer values)
	74	for i in range(len(time_vec)):
	75	b1[i] = quant16(b1[i],liminfB,limsupB)
	76	b2[i] = quant16(b2[i],liminfB,limsupB)
	77	b3[i] = quant16(b3[i],liminfB,limsupB)
	78	e1[i] = quant16(e1[i],liminfE,limsupE)
	79	e2[i] = quant16(e2[i],liminfE,limsupE)
	80
	81	#we filter the signal with a digital filter
	82
	83	#Filter 0
	84	a0_0 = -128
	85	a0_1 = 189
	86	a0_2 = -111
	87	b0_0 = 58
	88	b0_1 = -66
	89	b0_2 = 58
	90
	91	CF0 = Filter_Coefficients(a0_0, a0_1, a0_2, b0_0, b0_1, b0_2)
	92	F0 = Filter(CF0)
	93
	94	#Filter 1
	95	a1_0 = -128
	96	a1_1 = 162
	97	a1_2 = -81
	98	b1_0 = 58
	99	b1_1 = -57
	100	b1_2 = 58
	101
	102	CF1 = Filter_Coefficients(a1_0, a1_1, a1_2, b1_0, b1_1, b1_2)
	103	F1 = Filter(CF1)
	104
	105	#Filter 2
	106	a2_0 = -128
	107	a2_1 = 136
	108	a2_2 = -55
	109	b2_0 = 29
	110	b2_1 = -17
	111	b2_2 = 29
	112
	113	CF2 = Filter_Coefficients(a2_0, a2_1, a2_2, b2_0, b2_1, b2_2)
	114	F2 = Filter(CF2)
	115
	116	#Filter 3
	117	a3_0 = -128
	118	a3_1 = 114
	119	a3_2 = -33
	120	b3_0 = 15
	121	b3_1 = 4
	122	b3_2 = 15
	123
	124	CF3 = Filter_Coefficients(a3_0, a3_1, a3_2, b3_0, b3_1, b3_2)
	125	F3 = Filter(CF3)
	126
	127	#Filter 4
	128	a4_0 = -128
	129	a4_1 = 100
	130	a4_2 = -20
	131	b4_0 = 15
	132	b4_1 = 24
	133	b4_2 = 15
	134
	135	CF4 = Filter_Coefficients(a4_0, a4_1, a4_2, b4_0, b4_1, b4_2)
	136	F4 = Filter(CF4)
	137
	138	Filters = [F0,F1,F2,F3,F4]
	139	F_10k = FilterGlobal(Filters)
	140
	141	#this filtering assure us of not having aliasing problems we downsampling to fm = 24576 Hz
	142	b1_filtered = F_10k.convolve(b1)
	143	b2_filtered = F_10k.convolve(b2)
	144	b3_filtered = F_10k.convolve(b3)
	145	e1_filtered = F_10k.convolve(e1)
	146	e2_filtered = F_10k.convolve(e2)
	147
	148	#we downsample the signal to get to a sampling frequency fm = 24576Hz
	149	b1_24576 = b1_filtered[::4]
	150	b2_24576 = b2_filtered[::4]
	151	b3_24576 = b3_filtered[::4]
	152	e1_24576 = e1_filtered[::4]
	153	e2_24576 = e2_filtered[::4]
	154
	155	S_24576, w_24576 = Spectral_MatriceAvgTime_Int(b1_24576,b2_24576,b3_24576,e1_24576,e2_24576,24576,256)
	156	PB_24576 = PSD_B_Int(S_24576)
	157	PE_24576 = PSD_E_Int(S_24576)
	158
	159	# #Here we're doing simply a downsample, without filtering before. We should pass by a CIC filter instead.
	160	# b1_4096 = b1_24576[::6]
	161	# b2_4096 = b2_24576[::6]
	162	# b3_4096 = b3_24576[::6]
	163	# e1_4096 = e1_24576[::6]
	164	# e2_4096 = e2_24576[::6]
	165	#
	166	# #Here we're doing simply a downsample, without filtering before. We should pass by a CIC filter instead.
	167	# b1_256 = b1_24576[::96]
	168	# b2_256 = b2_24576[::96]
	169	# b3_256 = b3_24576[::96]
	170	# e1_256 = e1_24576[::96]
	171	# e2_256 = e2_24576[::96]
	172	#
	173	# #Here we're doing simply a downsample, without filtering before. We should pass by a CIC filter instead.
	174	# b1_16 = b1_4096[::256]
	175	# b2_16 = b2_4096[::256]
	176	# b3_16 = b3_4096[::256]
	177	# e1_16 = e1_4096[::256]
	178	# e2_16 = e2_4096[::256]

SRC/test_cases1.py created 644 +222 0

@@ -0,0 +1,222
	1	#Operations with float numbers#
	2
	3	import numpy as np
	4	from filters import *
	5	from bin16 import *
	6	from fft import *
	7	import math
	8	import cmath
	9	import matplotlib.pyplot as plt
	10	from basic_parameters import *
	11	import os
	12
	13	################################################# Simulations ##########################################################
	14
	15	#This creates a new folder to save the plots
	16	mypath = 'plots1'
	17	if not os.path.isdir(mypath):
	18	os.makedirs(mypath)
	19
	20	fs = 98304 # is the sampling frequency of our signals
	21	time_step = 1.0/fs
	22	t_ini = 0
	23	t_fin = 10
	24	time_vec = np.arange(t_ini,t_fin,time_step)
	25	fig = 0 #index for the plots1
	26
	27	f = 2500 # is the frequency of our signal of test (which gives the w = 2 * pi * f)
	28
	29	print "These are the results obtained for the basic parameters on a monochromatic electromagnetic wave"
	30	print "The rounding errors of the LFR are still not taken into account - we're operating with float numbers \n"
	31
	32	# we must give the vector n, which has the information about the direction of k
	33	n = np.array([1,1,1]) # for now we take this example
	34	n = n/np.linalg.norm(n)
	35	l = 3000 # this is the wavelength of our wave (which is not necessarily determined by c = l * f)
	36	k = 2.0 * math.pi / l
	37
	38	print "The n vector is:", n
	39	print "The wavelength l is:", l, "[m]"
	40	print "The wave number k is then:", k, "[m]^-1 \n"
	41
	42	# the values of 'a' and 'b' tell us the intensity of the magnetic field
	43	a = np.random.random() # we take a value at random between 0 and 1
	44	b = np.sqrt(1 - a*a) # we want to force a^2 + b^2 = 1
	45	phi = 0
	46
	47	print "The parameters who describe the magnetic field B are 'a' and 'b', who shall form a complex phasor given by: [a -b1j 0] exp(1j * (w*t + phi))"
	48	print "We're considering phi equal to ", phi
	49	print "The vector for the present case is: [", a, ", -j", b, ", 0]"
	50	print "('a' and 'b' were chosen randomly between 0 and 1) \n"
	51
	52	# the matrix R gives us an example of transformation to the SCM referencial
	53	# it was determined by the use of the normalized vector n and then two random vectors who form an orthogonal basis with it
	54	R = Matrix_SCM(n)
	55
	56	# now we transform the magnetic field to the SCM referencial
	57	B = MagneticComponents(a,b,2.0 * math.pi * f,R)
	58
	59	print "From the n vector we can construct an Orthonormal Matrix R capable of taking the components of B and taking them into the SCM reference"
	60	print "R = ", R
	61	print "The values of the components of B in the SCM reference are then: "
	62	print "B = ", B
	63	print " "
	64
	65	# now we have the caracteristics for our functions b1(t), b2(t) and b3(t) who shall enter in the LFR
	66	# b1(t) = cmath.polar(B[0])[0] * cos(w*t + cmath.polar(B[0])[1] + phi)
	67	# b2(t) = cmath.polar(B[1])[0] * cos(w*t + cmath.polar(B[1])[1] + phi)
	68	# b3(t) = cmath.polar(B[2])[0] * cos(w*t + cmath.polar(B[2])[1] + phi)
	69
	70	b1 = cmath.polar(B[0,0])[0] * np.cos(2.0 * math.pi * f * time_vec + cmath.polar(B[0,0])[1] + phi)
	71	b2 = cmath.polar(B[1,0])[0] * np.cos(2.0 * math.pi * f * time_vec + cmath.polar(B[1,0])[1] + phi)
	72	b3 = cmath.polar(B[2,0])[0] * np.cos(2.0 * math.pi * f * time_vec + cmath.polar(B[2,0])[1] + phi)
	73
	74	print "Now we finally have the functions b1(t), b2(t) and b3(t) who are going to serve as examples for inputs in the LFR analyzer"
	75	print "b1(t) = ", cmath.polar(B[0,0])[0], "cos(w*t + ", cmath.polar(B[0,0])[1], " + ", phi, ")"
	76	print "b2(t) = ", cmath.polar(B[1,0])[0], "cos(w*t + ", cmath.polar(B[1,0])[1], " + ", phi, ")"
	77	print "b3(t) = ", cmath.polar(B[2,0])[0], "cos(w*t + ", cmath.polar(B[2,0])[1], " + ", phi, ")"
	78	print "(the value of w is,", 2.0 * math.pi * f, ")"
	79	print " "
	80
	81	graph = raw_input("Do you wish to plot the graphs for the magnetic field signals? (Y/N): ")
	82	if graph == 'Y':
	83	print "Next we plot the graph for the three magnetic field signals: b1(t) [red], b2(t) [blue], b3(t) [green]"
	84	fig = fig + 1
	85	plt.figure(fig)
	86	p1, = plt.plot(1000*time_vec[0:len(time_vec)/2000],b1[0:len(time_vec)/2000],'r')
	87	p2, = plt.plot(1000*time_vec[0:len(time_vec)/2000],b2[0:len(time_vec)/2000],'b')
	88	p3, = plt.plot(1000*time_vec[0:len(time_vec)/2000],b3[0:len(time_vec)/2000],'g')
	89	plt.title("The three input signals for the Magnetic Field")
	90	plt.xlabel("time [ms]")
	91	plt.ylabel("b(t)")
	92	plt.legend([p1, p2, p3], ["b1(t)", "b2(t)", "b3(t)"])
	93	plt.show()
	94	print " "
	95
	96	# the component of E who's parallel to n is of our choice. It represents the influence of the presence of charges where the electromagnetic
	97	# wave is travelling. However, the component of E who's perpendicular to n is determined by B and n as follows
	98
	99	E_orth = -1.0 * 2.0 * math.pi * f/k * np.cross(n,np.array(B.transpose()))
	100	E_para = np.zeros(3) # for now, we'll consider that E_para is zero
	101	E = E_orth + E_para
	102	E = E.transpose()
	103
	104	e1 = cmath.polar(E[0,0])[0] * np.cos(2.0 * math.pi * f * time_vec + cmath.polar(E[0,0])[1])
	105	e2 = cmath.polar(E[1,0])[0] * np.cos(2.0 * math.pi * f * time_vec + cmath.polar(E[1,0])[1])
	106
	107	print "The electric field has two components: one parallel to the n vector (E_para) and one orthogonal (E_orth)"
	108	print "E_para represents the degree of polarization of the medium where the wave is propagating"
	109	print "E_para = ", E_para
	110	print "E_orth is fully determined by the B vector representing the magnetic field"
	111	print "E_orth = -w/k * (n x B)"
	112	print "E_orth = ", E_orth
	113	print " "
	114
	115	print "We have then the functions e1(t), e2(t)who are going to serve as examples for inputs in the LFR analyzer"
	116	print "e1(t) = ", cmath.polar(E[0,0])[0], "cos(w*t + ", cmath.polar(E[0,0])[1], " + ", phi, ")"
	117	print "e2(t) = ", cmath.polar(E[1,0])[0], "cos(w*t + ", cmath.polar(E[1,0])[1], " + ", phi, ")"
	118	print "(the value of w is,", 2.0 * math.pi * f, ")"
	119	print " "
	120
	121	graph = raw_input("Do you wish to plot the graphs for the electric field signals? (Y/N): ")
	122	if graph == 'Y':
	123	print "Next we plot the graph for two of the electric field signals: e1(t) [red], e2(t) [blue]"
	124	print "(The electric field is being considered in the SCM referencial, which is not totally right, but should work for now)"
	125	fig = fig + 1
	126	plt.figure(fig)
	127	p1, = plt.plot(1000*time_vec[0:len(time_vec)/2000],e1[0:len(time_vec)/2000],'r')
	128	p2, = plt.plot(1000*time_vec[0:len(time_vec)/2000],e2[0:len(time_vec)/2000],'b')
	129	plt.title("The two input signals for the Electric Field")
	130	plt.xlabel("time [ms]")
	131	plt.ylabel("e(t)")
	132	plt.legend([p1, p2], ["e1(t)", "e2(t)"])
	133	plt.show()
	134	print " "
	135
	136	#This is the number of points used in the FFT
	137	N = 256
	138
	139	#We may add some noise to the input signals
	140	b1 = b1 + 10.0 * np.random.random(len(b1))/100
	141	b2 = b2 + 10.0 * np.random.random(len(b2))/100
	142	b3 = b3 + 10.0 * np.random.random(len(b3))/100
	143	e1 = e1 + 10.0 * 10000000.0 * np.random.random(len(e1))/100
	144	e2 = e2 + 10.0 * 10000000.0 * np.random.random(len(e2))/100
	145
	146	graph = raw_input("Do you wish to plot the graphs for the noised magnetic and electric field signals? (Y/N): ")
	147	if graph == 'Y':
	148	print "Next we plot the graph for the three magnetic field signals: b1(t) [red], b2(t) [blue], b3(t) [green]"
	149	print "As well as for the two electric field signals: e1(t) [red], e2(t) [blue]"
	150	fig = fig + 1
	151	plt.figure(fig)
	152	p1, = plt.plot(1000*time_vec[0:len(time_vec)/2000],b1[0:len(time_vec)/2000],'r')
	153	p2, = plt.plot(1000*time_vec[0:len(time_vec)/2000],b2[0:len(time_vec)/2000],'b')
	154	p3, = plt.plot(1000*time_vec[0:len(time_vec)/2000],b3[0:len(time_vec)/2000],'g')
	155	plt.title("The three input signals for the Magnetic Field")
	156	plt.xlabel("time [ms]")
	157	plt.ylabel("b(t)")
	158	plt.legend([p1, p2, p3], ["b1(t)", "b2(t)", "b3(t)"])
	159	plt.show()
	160
	161	fig = fig + 1
	162	plt.figure(fig)
	163	p1, = plt.plot(1000*time_vec[0:len(time_vec)/2000],e1[0:len(time_vec)/2000],'r')
	164	p2, = plt.plot(1000*time_vec[0:len(time_vec)/2000],e2[0:len(time_vec)/2000],'b')
	165	plt.title("The two input signals for the Electric Field")
	166	plt.xlabel("time [ms]")
	167	plt.ylabel("e(t)")
	168	plt.legend([p1, p2], ["e1(t)", "e2(t)"])
	169	plt.show()
	170	print " "
	171
	172	S, w = Spectral_Matrice(b1,b2,b3,e1,e2,fs,N,offset = 0)
	173
	174	print "Now we calculate the Spectral Matrice for the ensemble of signals entering the LFR."
	175	print "(For now, we're just taking one spectral matrice, for all the values of frequencies. The next step will be to average them in time and frequence)"
	176	print "Then, from this Spectral Matrice we calculate the basic parameters:"
	177
	178	PB = PSD_B(S)
	179	print "Power spectrum density of the Magnetic Field"
	180	fig = fig + 1
	181	plt.figure(fig)
	182	plt.plot(w/(2*math.pi),PB)
	183	plt.title("Power spectrum density of the Magnetic Field")
	184	plt.xlabel("frequency [Hz]")
	185	plt.ylabel("PB(w)")
	186	plt.axvline(x = f, color = 'r', linestyle='--')
	187	plt.savefig('plots1/PB.png')
	188
	189	PE = PSD_E(S)
	190	print "Power spectrum density of the Electric Field"
	191	fig = fig + 1
	192	plt.figure(fig)
	193	plt.plot(w/(2*math.pi),PE)
	194	plt.title("Power spectrum density of the Electric Field")
	195	plt.xlabel("frequency [Hz]")
	196	plt.ylabel("PE(w)")
	197	plt.axvline(x = f, color = 'r', linestyle='--')
	198	plt.savefig('plots1/PE.png')
	199
	200	e = ellipticity(S)
	201	print "Ellipticity of the electromagnetic wave"
	202	fig = fig + 1
	203	plt.figure(fig)
	204	plt.plot(w/(2*math.pi),e)
	205	plt.title("Ellipticity of the electromagnetic wave")
	206	plt.xlabel("frequency [Hz]")
	207	plt.ylabel("e(w)")
	208	plt.axvline(x = f, color = 'r', linestyle = '--')
	209	plt.axhline(y = 2.0/(a/b + b/a), color = 'r', linestyle='--')
	210	plt.savefig('plots1/e.png')
	211
	212	d = degree_polarization(S)
	213	print "Degree of polarization of the electromagnetic wave"
	214	fig = fig + 1
	215	plt.figure(fig)
	216	plt.plot(w/(2*math.pi),d)
	217	plt.title("Degree of polarization of the electromagnetic wave")
	218	plt.xlabel("frequency [Hz]")
	219	plt.ylabel("d(w)")
	220	plt.savefig('plots1/d.png')
	221
	222	#[n1,n2,n3] = normal_vector(S) No newline at end of file

SRC/test_cases2.py created 644 +238 0

@@ -0,0 +1,238
	1	#Operations with float numbers#
	2
	3	import numpy as np
	4	from filters import *
	5	from bin16 import *
	6	from fft import *
	7	import math
	8	import cmath
	9	import matplotlib.pyplot as plt
	10	from basic_parameters import *
	11	import os
	12
	13	################################################# Simulations ##########################################################
	14
	15	#This creates a new folder to save the plots
	16	mypath = 'plots2'
	17	if not os.path.isdir(mypath):
	18	os.makedirs(mypath)
	19
	20	fs = 98304 # is the sampling frequency of our signals
	21	time_step = 1.0/fs
	22	t_ini = 0
	23	t_fin = 10
	24	time_vec = np.arange(t_ini,t_fin,time_step)
	25	fig = 0 #index for the plots
	26
	27	f = 500 # is the frequency of our signal of test (which gives the w = 2 * pi * f)
	28
	29	# we must give the vector n, which has the information about the direction of k
	30	n = np.array([1,1,1]) # for now we take this example
	31	n = n/np.linalg.norm(n)
	32	l = 3000 # this is the wavelength of our wave (which is not necessarily determined by c = l * f)
	33	k = 2.0 * math.pi / l
	34
	35	# the values of 'a' and 'b' tell us the intensity of the magnetic field
	36	a = np.random.random() # we take a value at random between 0 and 1
	37	b = np.sqrt(1 - a*a) # we want to force a^2 + b^2 = 1
	38	phi = 0
	39
	40	# the matrix R gives us an example of transformation to the SCM referencial
	41	# it was determined by the use of the normalized vector n and then two random vectors who form an orthogonal basis with it
	42	R = Matrix_SCM(n)
	43
	44	# now we transform the magnetic field to the SCM referencial
	45	B = MagneticComponents(a,b,2.0 * math.pi * f,R)
	46
	47	# now we have the caracteristics for our functions b1(t), b2(t) and b3(t) who shall enter in the LFR
	48	# b1(t) = cmath.polar(B[0])[0] * cos(w*t + cmath.polar(B[0])[1] + phi)
	49	# b2(t) = cmath.polar(B[1])[0] * cos(w*t + cmath.polar(B[1])[1] + phi)
	50	# b3(t) = cmath.polar(B[2])[0] * cos(w*t + cmath.polar(B[2])[1] + phi)
	51
	52	b1 = cmath.polar(B[0,0])[0] * np.cos(2.0 * math.pi * f * time_vec + cmath.polar(B[0,0])[1] + phi)
	53	b2 = cmath.polar(B[1,0])[0] * np.cos(2.0 * math.pi * f * time_vec + cmath.polar(B[1,0])[1] + phi)
	54	b3 = cmath.polar(B[2,0])[0] * np.cos(2.0 * math.pi * f * time_vec + cmath.polar(B[2,0])[1] + phi)
	55
	56	# the component of E who's parallel to n is of our choice. It represents the influence of the presence of charges where the electromagnetic
	57	# wave is travelling. However, the component of E who's perpendicular to n is determined by B and n as follows
	58
	59	E_orth = -1.0 * 2.0 * math.pi * f/k * np.cross(n,np.array(B.transpose()))
	60	E_para = np.zeros(3) # for now, we'll consider that E_para is zero
	61	E = E_orth + E_para
	62	E = E.transpose()
	63
	64	e1 = cmath.polar(E[0,0])[0] * np.cos(2.0 * math.pi * f * time_vec + cmath.polar(E[0,0])[1])
	65	e2 = cmath.polar(E[1,0])[0] * np.cos(2.0 * math.pi * f * time_vec + cmath.polar(E[1,0])[1])
	66
	67	#This is the number of points used in the FFT
	68	N = 256
	69
	70	#We may add some noise to the input signals
	71	# b1 = b1 + np.random.random(len(b1))/10
	72	# b2 = b2 + np.random.random(len(b2))/10
	73	# b3 = b3 + np.random.random(len(b3))/10
	74	# e1 = e1 + np.random.random(len(e1))/10
	75	# e2 = e2 + np.random.random(len(e2))/10
	76
	77	fm1 = 24576
	78	fm2 = 4096
	79	fm3 = 256
	80
	81	S1, w1 = Spectral_MatriceAvgTime(b1[::4],b2[::4],b3[::4],e1[::4],e2[::4],fm1,N)
	82	S2, w2 = Spectral_MatriceAvgTime(b1[::24],b2[::24],b3[::24],e1[::24],e2[::24],fm2,N)
	83	S3, w3 = Spectral_MatriceAvgTime(b1[::384],b2[::384],b3[::384],e1[::384],e2[::384],fm3,N)
	84
	85	print "S1"
	86
	87	PB1 = PSD_B(S1)
	88	print "Power spectrum density of the Magnetic Field"
	89	fig = fig + 1
	90	plt.figure(fig)
	91	plt.plot(w1/(2*math.pi),PB1)
	92	plt.title("Power spectrum density of the Magnetic Field (fm = %s Hz)" % fm1)
	93	plt.xlabel("frequency [Hz]")
	94	plt.ylabel("PB(w)")
	95	plt.axvline(x = f, color = 'r', linestyle='--')
	96	plt.savefig('plots2/PB1.pdf')
	97
	98	PE1 = PSD_E(S1)
	99	print "Power spectrum density of the Electric Field"
	100	fig = fig + 1
	101	plt.figure(fig)
	102	plt.plot(w1/(2*math.pi),PE1)
	103	plt.title("Power spectrum density of the Electric Field (fm = %s Hz)" % fm1)
	104	plt.xlabel("frequency [Hz]")
	105	plt.ylabel("PE(w)")
	106	plt.axvline(x = f, color = 'r', linestyle='--')
	107	plt.savefig('plots2/PE1.pdf')
	108
	109	e1 = ellipticity(S1)
	110	print "Ellipticity of the electromagnetic wave"
	111	fig = fig + 1
	112	plt.figure(fig)
	113	plt.plot(w1/(2*math.pi),e1)
	114	plt.title("Ellipticity of the electromagnetic wave (fm = %s Hz)" % fm1)
	115	plt.xlabel("frequency [Hz]")
	116	plt.ylabel("e(w)")
	117	plt.axvline(x = f, color = 'r', linestyle = '--')
	118	plt.axhline(y = 2.0/(a/b + b/a), color = 'r', linestyle='--')
	119	plt.savefig('plots2/e1.pdf')
	120
	121	d1 = degree_polarization(S1)
	122	print "Degree of polarization of the electromagnetic wave"
	123	fig = fig + 1
	124	plt.figure(fig)
	125	plt.plot(w1/(2*math.pi),d1)
	126	plt.title("Degree of polarization of the electromagnetic wave (fm = %s Hz)" % fm1)
	127	plt.xlabel("frequency [Hz]")
	128	plt.ylabel("d(w)")
	129	plt.axvline(x = f, color = 'r', linestyle='--')
	130	plt.savefig('plots2/d1.pdf')
	131
	132	#[n11,n21,n31] = normal_vector(S1)
	133
	134	print "S2"
	135
	136	PB2 = PSD_B(S2)
	137	print "Power spectrum density of the Magnetic Field"
	138	fig = fig + 1
	139	plt.figure(fig)
	140	plt.plot(w2/(2*math.pi),PB2)
	141	plt.title("Power spectrum density of the Magnetic Field (fm = %s Hz)" % fm2)
	142	plt.xlabel("frequency [Hz]")
	143	plt.ylabel("PB(w)")
	144	plt.axvline(x = f, color = 'r', linestyle='--')
	145	plt.savefig('plots2/PB2.pdf')
	146	plt.close()
	147
	148	PE2 = PSD_E(S2)
	149	print "Power spectrum density of the Electric Field"
	150	fig = fig + 1
	151	plt.figure(fig)
	152	plt.plot(w2/(2*math.pi),PE2)
	153	plt.title("Power spectrum density of the Electric Field (fm = %s Hz)" % fm2)
	154	plt.xlabel("frequency [Hz]")
	155	plt.ylabel("PE(w)")
	156	plt.axvline(x = f, color = 'r', linestyle='--')
	157	plt.savefig('plots2/PE2.pdf')
	158	plt.close()
	159
	160	e2 = ellipticity(S2)
	161	print "Ellipticity of the electromagnetic wave"
	162	fig = fig + 1
	163	plt.figure(fig)
	164	plt.plot(w2/(2*math.pi),e2)
	165	plt.title("Ellipticity of the electromagnetic wave (fm = %s Hz)" % fm2)
	166	plt.xlabel("frequency [Hz]")
	167	plt.ylabel("e(w)")
	168	plt.axvline(x = f, color = 'r', linestyle = '--')
	169	plt.axhline(y = 2.0/(a/b + b/a), color = 'r', linestyle='--')
	170	plt.savefig('plots2/e2.pdf')
	171	plt.close()
	172
	173	d2 = degree_polarization(S2)
	174	print "Degree of polarization of the electromagnetic wave"
	175	fig = fig + 1
	176	plt.figure(fig)
	177	plt.plot(w2/(2*math.pi),d2)
	178	plt.title("Degree of polarization of the electromagnetic wave (fm = %s Hz)" % fm2)
	179	plt.xlabel("frequency [Hz]")
	180	plt.ylabel("d(w)")
	181	plt.axvline(x = f, color = 'r', linestyle='--')
	182	plt.savefig('plots2/d2.pdf')
	183	plt.close()
	184
	185	#[n12,n22,n32] = normal_vector(S2)
	186
	187	print "S3"
	188
	189	PB3 = PSD_B(S3)
	190	print "Power spectrum density of the Magnetic Field"
	191	fig = fig + 1
	192	plt.figure(fig)
	193	plt.plot(w3/(2*math.pi),PB3)
	194	plt.title("Power spectrum density of the Magnetic Field (fm = %s Hz)" % fm3)
	195	plt.xlabel("frequency [Hz]")
	196	plt.ylabel("PB(w)")
	197	plt.axvline(x = f, color = 'r', linestyle='--')
	198	plt.savefig('plots2/PB3.pdf')
	199	plt.close()
	200
	201	PE3 = PSD_E(S3)
	202	print "Power spectrum density of the Electric Field "
	203	fig = fig + 1
	204	plt.figure(fig)
	205	plt.plot(w3/(2*math.pi),PE3)
	206	plt.title("Power spectrum density of the Electric Field (fm = %s Hz)" % fm3)
	207	plt.xlabel("frequency [Hz]")
	208	plt.ylabel("PE(w)")
	209	plt.axvline(x = f, color = 'r', linestyle='--')
	210	plt.savefig('plots2/PE3.pdf')
	211	plt.close()
	212
	213	e3 = ellipticity(S3)
	214	print "Ellipticity of the electromagnetic wave"
	215	fig = fig + 1
	216	plt.figure(fig)
	217	plt.plot(w3/(2*math.pi),e3)
	218	plt.title("Ellipticity of the electromagnetic wave (fm = %s Hz)" % fm3)
	219	plt.xlabel("frequency [Hz]")
	220	plt.ylabel("e(w)")
	221	plt.axvline(x = f, color = 'r', linestyle = '--')
	222	plt.axhline(y = 2.0/(a/b + b/a), color = 'r', linestyle='--')
	223	plt.savefig('plots2/e3.pdf')
	224	plt.close()
	225
	226	d3 = degree_polarization(S3)
	227	print "Degree of polarization of the electromagnetic wave"
	228	fig = fig + 1
	229	plt.figure(fig)
	230	plt.plot(w3/(2*math.pi),d3)
	231	plt.title("Degree of polarization of the electromagnetic wave (fm = %s Hz)" % fm3)
	232	plt.xlabel("frequency [Hz]")
	233	plt.ylabel("d(w)")
	234	plt.axvline(x = f, color = 'r', linestyle='--')
	235	plt.savefig('plots2/d3.pdf')
	236	plt.close()
	237
	238	#[n13,n23,n33] = normal_vector(S3) No newline at end of file

SRC/test_cases3.py created 644 +36 0

@@ -0,0 +1,36
	1	#Operations with float numbers#
	2
	3	import numpy as np
	4	from filters import *
	5	from bin16 import *
	6	from fft import *
	7	import math
	8	import cmath
	9	import matplotlib.pyplot as plt
	10	from basic_parameters import *
	11
	12	# The parameters who describe the magnetic field B are 'a' and 'b',
	13	# who shall form a complex phasor given by: [a -b1j 0] exp(1j * (w*t + phi))
	14	a = np.random.random() # we take a value at random between 0 and 1
	15	b = np.sqrt(1 - a*a) # we want to force a^2 + b^2 = 1
	16
	17	# 'f' is the frequency of oscillation of our electromagnetic wave, who's monochromatic for now
	18	f = 2500
	19
	20	# 'l' is the wavelength of the electromagnetic wave
	21	l = 3000
	22
	23	# is the propagating vector, who's related to the k vector
	24	n = [1,1,1]
	25
	26	# is a vector who tells us the degree of polarization of the medium where the electromagnetic wave is being propagated
	27	E_para = [0,0,0]
	28
	29	# 'fm' is the original sampling frequency (the one who enters the LFRß)
	30	fm = 98304
	31
	32	# this function will give us:
	33	# - A list of Spectral Matrices: S = [S1, S2, S3]
	34	# - A list of Angular Frequencies: w = [w1, w2, w3]
	35	# - The possibility of plotting and saving the graphs for each value of fm (24576 Hz, 4096 Hz, 256 Hz)
	36	S,w = SpectralMatrice_Monochromatic(a,b,n,E_para,f,fm,l) No newline at end of file

SRC/test_cases4.py created 644 +110 0

@@ -0,0 +1,110
	1	#Operations with integers#
	2
	3	import numpy as np
	4	from bin16 import *
	5	from fft import *
	6	import math
	7	import cmath
	8	import matplotlib.pyplot as plt
	9	from basic_parameters_Int import *
	10	from basic_parameters import *
	11
	12	#This is the number of points used for the FFT
	13	N = 256
	14
	15	#These are the limits of our quantization
	16	liminfB = -1
	17	limsupB = +1
	18	liminfE = -10000000
	19	limsupE = +10000000
	20
	21	#############################
	22
	23	#Example of signals taken from a simulation
	24	a = np.random.rand()
	25	b = np.sqrt(1 - a*a)
	26	n = [1,1,1]
	27	E_para = [0,0,0]
	28	f = 1500
	29	fs = 98304
	30	l = 30000
	31
	32	[b1,b2,b3,e1,e2] = SignalsWaveMonochromatic(a,b,n,E_para,f,fs,l)
	33
	34	#############################
	35
	36	#We add noise
	37	b1 = b1 + np.random.uniform(1.0liminfB/100,1.0limsupB/100,len(b1))
	38	b2 = b2 + np.random.uniform(1.0liminfB/100,1.0limsupB/100,len(b2))
	39	b3 = b3 + np.random.uniform(1.0liminfB/100,1.0limsupB/100,len(b3))
	40	e1 = e1 + np.random.uniform(1.0liminfE/100,1.0limsupE/100,len(e1))
	41	e2 = e2 + np.random.uniform(1.0liminfE/100,1.0limsupE/100,len(e2))
	42
	43	#This will take our float-number signals and quantize them into integers going from -32768 to 32767
	44	#The liminf and limsup is what decides the maximum and minimum values for our quantizations
	45	#quant16 is a function from the bin16 library
	46
	47	b1Q = np.zeros(len(b1))
	48	b2Q = np.zeros(len(b2))
	49	b3Q = np.zeros(len(b3))
	50	e1Q = np.zeros(len(e1))
	51	e2Q = np.zeros(len(e2))
	52
	53	for i in range(len(b1)):
	54	b1Q[i] = quant16(b1[i],liminfB,limsupB)
	55	b2Q[i] = quant16(b2[i],liminfB,limsupB)
	56	b3Q[i] = quant16(b3[i],liminfB,limsupB)
	57	e1Q[i] = quant16(e1[i],liminfE,limsupE)
	58	e2Q[i] = quant16(e2[i],liminfE,limsupE)
	59
	60	#We calculate the Spectral Matrices (averaged in time) for the case of float-numbers operations and also for the integer-numbers operations
	61	Si,wi = Spectral_MatriceAvgTime_Int(b1Q[::24],b2Q[::24],b3Q[::24],e1Q[::24],e2Q[::24],4096,N)
	62	Sf,wf = Spectral_MatriceAvgTime(b1[::24],b2[::24],b3[::24],e1[::24],e2[::24],4096,N)
	63
	64	[PBf,PEf,df,ef,nf] = BasicParameters(Sf,wf)
	65	[PBi,PEi,di,ei,ni] = BasicParameters_Int(Si,wi)
	66
	67	mypath = "plotsFloatVsInt"
	68	if not os.path.isdir(mypath):
	69	os.makedirs(mypath)
	70
	71	fm = 4096
	72
	73	alfa = 1.0 * math.pow(2,30)
	74	plt.plot(wf/(2*math.pi),PBf,'g')
	75	plt.scatter(wf/(2math.pi),1.0 PBi/alfa)
	76	plt.title("Power spectrum density of the Magnetic Field (fm = %s Hz)" % fm)
	77	plt.xlabel("frequency [Hz]")
	78	plt.ylabel("PB(w)")
	79	plt.savefig('%s/PB.png' % mypath)
	80	plt.close()
	81
	82	alfa = 1.0 * math.pow(2,30)/math.pow(10,14)
	83	plt.plot(wf/(2*math.pi),PEf,'g')
	84	plt.scatter(wf/(2math.pi),1.0 PEi/alfa)
	85	plt.title("Power spectrum density of the Electric Field (fm = %s Hz)" % fm)
	86	plt.xlabel("frequency [Hz]")
	87	plt.ylabel("PE(w)")
	88	plt.savefig('%s/PE.png' % mypath)
	89	plt.close()
	90
	91	alfa = 16
	92	plt.plot(wf/(2*math.pi),ef,'g')
	93	plt.scatter(wf/(2math.pi),1.0 ei/alfa)
	94	plt.title("Ellipticity of the electromagnetic wave (fm = %s Hz)" % fm)
	95	plt.xlabel("frequency [Hz]")
	96	plt.ylabel("e(w)")
	97	plt.savefig('%s/e.png' % mypath)
	98	plt.close()
	99
	100	alfa = 8
	101	plt.plot(wf/(2*math.pi),df, 'g')
	102	plt.scatter(wf/(2math.pi),1.0 di/alfa)
	103	plt.title("Degree of polarization of the electromagnetic wave (fm = %s Hz)" % fm)
	104	plt.xlabel("frequency [Hz]")
	105	plt.ylabel("d(w)")
	106	plt.savefig('%s/d.png' % mypath)
	107	plt.close()
	108
	109
	110

SRC/testing_FFT.py created 644 +40 0

@@ -0,0 +1,40
	1	import matplotlib.pyplot as plt
	2	from bin16 import *
	3	from fft import *
	4	import math
	5
	6	#These are the first parameters for choosing the time step of our vector
	7	fs = 98304
	8	time_step = 1.0/fs
	9	t_ini = 0
	10	t_fin = 0.25
	11	time_vec = np.arange(t_ini, t_fin, time_step)
	12
	13	#The two frequences we're going to use for our different signals
	14	f1 = 4000
	15	f2 = 800
	16
	17	#This ensures me of getting the higher value of N who's a power of 2 and smaller than the length of the time_vec
	18	N = int(math.pow(2,int(math.log(len(time_vec),2))))
	19
	20	x1 = 0.5 * np.cos(2 * np.pi * f1 * time_vec)
	21	x2 = 0.5 * np.cos(2 * np.pi * f2 * time_vec)
	22
	23	liminf = -1
	24	limsup = +1
	25
	26	#This will take our float-number signals, who were going from -0.5 to 0.5 (they were all sines), and quantize them into integers going from -32768 to 32767
	27	#The liminf and limsup is what decides the maximum and minimum values for our quantizations
	28	#quant16 is a function from the bin16 library
	29	for i in range(len(time_vec)):
	30	x1[i] = quant16(x1[i],liminf,limsup)
	31	x2[i] = quant16(x2[i],liminf,limsup)
	32
	33	#Now we obtain the values for the Fourier Transforms of both signals, with integer values as well (rounded in the function fft_CT)
	34	X1 = fft_CT(x1[0:N])
	35	X2 = fft_CT(x2[0:N])
	36
	37	plt.figure(1)
	38	plt.plot(abs(X1),'r')
	39	plt.plot(abs(X2),'g')
	40	plt.show() No newline at end of file

General Comments 0

Write
Preview

You need to be logged in to leave comments. Login now

No TODOs yet

	Site-wide shortcuts
/	Use quick search box
g h	Goto home page
g g	Goto my private gists page
g G	Goto my public gists page
g 0-9	Goto bookmarked items from 0-9
n r	New repository page
n g	New gist page

	Repositories
g s	Goto summary page
g c	Goto changelog page
g f	Goto files page
g F	Goto files page with file search activated
g p	Goto pull requests page
g o	Goto repository settings
g O	Goto repository access permissions settings
t s	Toggle sidebar on some pages