@@ -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