@@ -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 = [[B1*B1.conjugate(), B1*B2.conjugate(), B1*B3.conjugate(), B1*E1.conjugate(), B1*E2.conjugate()], | |||
|
40 | [B2*B1.conjugate(), B2*B2.conjugate(), B2*B3.conjugate(), B2*E1.conjugate(), B2*E2.conjugate()], | |||
|
41 | [B3*B1.conjugate(), B3*B2.conjugate(), B3*B3.conjugate(), B3*E1.conjugate(), B3*E2.conjugate()], | |||
|
42 | [E1*B1.conjugate(), E1*B2.conjugate(), E1*B3.conjugate(), E1*E1.conjugate(), E1*E2.conjugate()], | |||
|
43 | [E2*B1.conjugate(), E2*B2.conjugate(), E2*B3.conjugate(), E2*E1.conjugate(), E2*E2.conjugate()]] | |||
|
44 | ||||
|
45 | w = 2*math.pi*fm*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].imag*S[0][1].imag + S[0][2].imag*S[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((3*TrSCM2 - TrSCM*TrSCM)/(2*TrSCM*TrSCM)) | |||
|
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].imag*S[0][1].imag + S[0][2].imag*S[0][2].imag + S[1][2].imag*S[1][2].imag) | |||
|
114 | n2 = -1.0 * S[0][2].imag/np.sqrt(S[0][1].imag*S[0][1].imag + S[0][2].imag*S[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 |
@@ -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 = [[B1*B1.conjugate(), B1*B2.conjugate(), B1*B3.conjugate(), B1*E1.conjugate(), B1*E2.conjugate()], | |||
|
47 | [B2*B1.conjugate(), B2*B2.conjugate(), B2*B3.conjugate(), B2*E1.conjugate(), B2*E2.conjugate()], | |||
|
48 | [B3*B1.conjugate(), B3*B2.conjugate(), B3*B3.conjugate(), B3*E1.conjugate(), B3*E2.conjugate()], | |||
|
49 | [E1*B1.conjugate(), E1*B2.conjugate(), E1*B3.conjugate(), E1*E1.conjugate(), E1*E2.conjugate()], | |||
|
50 | [E2*B1.conjugate(), E2*B2.conjugate(), E2*B3.conjugate(), E2*E1.conjugate(), E2*E2.conjugate()]] | |||
|
51 | ||||
|
52 | #this we can keep as float | |||
|
53 | w = 2*math.pi*fm*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].imag*S[0][1].imag + S[0][2].imag*S[0][2].imag + S[1][2].imag*S[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.0*TrSCM2 - 1.0*TrSCM*TrSCM)/(2.0*TrSCM*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].imag*S[0][1].imag + S[0][2].imag*S[0][2].imag + S[1][2].imag*S[1][2].imag) | |||
|
137 | n2 = -1.0 * S[0][2].imag/np.sqrt(S[0][1].imag*S[0][1].imag + S[0][2].imag*S[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 |
@@ -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 |
@@ -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,-2j*math.pi*k/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,-2j*math.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 |
@@ -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 |
@@ -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 |
@@ -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 -b*1j 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] |
@@ -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 -b*1j 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 |
@@ -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 |
@@ -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 -b*1j 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 |
@@ -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.0*liminfB/100,1.0*limsupB/100,len(b1)) | |||
|
38 | b2 = b2 + np.random.uniform(1.0*liminfB/100,1.0*limsupB/100,len(b2)) | |||
|
39 | b3 = b3 + np.random.uniform(1.0*liminfB/100,1.0*limsupB/100,len(b3)) | |||
|
40 | e1 = e1 + np.random.uniform(1.0*liminfE/100,1.0*limsupE/100,len(e1)) | |||
|
41 | e2 = e2 + np.random.uniform(1.0*liminfE/100,1.0*limsupE/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/(2*math.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/(2*math.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/(2*math.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/(2*math.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 |
@@ -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
You need to be logged in to leave comments.
Login now