|
|
#the signal processing chain in the LFR analyser.
|
|
|
|
|
|
import numpy as np
|
|
|
from filters import *
|
|
|
from bin16 import *
|
|
|
from fft import *
|
|
|
import math
|
|
|
import cmath
|
|
|
import matplotlib.pyplot as plt
|
|
|
from basic_parameters import *
|
|
|
import os
|
|
|
|
|
|
# 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))
|
|
|
a = np.random.random() # we take a value at random between 0 and 1
|
|
|
b = np.sqrt(1 - a*a) # we want to force a^2 + b^2 = 1
|
|
|
|
|
|
# 'f' is the frequency of oscillation of our electromagnetic wave, who's monochromatic for now
|
|
|
f = 2500
|
|
|
# 'l' is the wavelength of the electromagnetic wave
|
|
|
l = 3000
|
|
|
k = 2.0 * math.pi / l
|
|
|
# is the propagating vector, who's related to the k vector
|
|
|
n = [1,1,1]
|
|
|
# is a vector who tells us the degree of polarization of the medium where the electromagnetic wave is being propagated
|
|
|
E_para = [0,0,0]
|
|
|
# 'fs' is the original sampling frequency (the one who enters the LFR)
|
|
|
fs = 98304
|
|
|
t_ini = 0
|
|
|
t_fin = 10
|
|
|
time_step = 1.0/fs
|
|
|
time_vec = np.arange(t_ini, t_fin, time_step)
|
|
|
|
|
|
# the matrix R gives us an example of transformation to the SCM referencial
|
|
|
# it was determined by the use of the normalized vector n and then two random vectors who form an orthogonal basis with it
|
|
|
R = Matrix_SCM(n)
|
|
|
|
|
|
# now we transform the magnetic field to the SCM referencial
|
|
|
B = MagneticComponents(a,b,2.0 * math.pi * f,R)
|
|
|
|
|
|
# now we have the caracteristics for our functions b1(t), b2(t) and b3(t) who shall enter in the LFR
|
|
|
b1 = cmath.polar(B[0,0])[0] * np.cos(2.0 * math.pi * f * time_vec + cmath.polar(B[0,0])[1])
|
|
|
b2 = cmath.polar(B[1,0])[0] * np.cos(2.0 * math.pi * f * time_vec + cmath.polar(B[1,0])[1])
|
|
|
b3 = cmath.polar(B[2,0])[0] * np.cos(2.0 * math.pi * f * time_vec + cmath.polar(B[2,0])[1])
|
|
|
|
|
|
# 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
|
|
|
# wave is travelling. However, the component of E who's perpendicular to n is determined by B and n as follows
|
|
|
|
|
|
E_orth = -1.0 * 2.0 * math.pi * f/k * np.cross(n,np.array(B.transpose()))
|
|
|
E = E_orth + E_para
|
|
|
E = E.transpose()
|
|
|
|
|
|
e1 = cmath.polar(E[0,0])[0] * np.cos(2.0 * math.pi * f * time_vec + cmath.polar(E[0,0])[1])
|
|
|
e2 = cmath.polar(E[1,0])[0] * np.cos(2.0 * math.pi * f * time_vec + cmath.polar(E[1,0])[1])
|
|
|
|
|
|
liminfB = -1
|
|
|
limsupB = +1
|
|
|
liminfE = -10000000
|
|
|
limsupE = +10000000
|
|
|
|
|
|
b1 = b1 + np.random.uniform(liminfB/100,limsupB/100,len(time_vec))
|
|
|
b2 = b2 + np.random.uniform(liminfB/100,limsupB/100,len(time_vec))
|
|
|
b3 = b3 + np.random.uniform(liminfB/100,limsupB/100,len(time_vec))
|
|
|
e1 = e1 + np.random.uniform(liminfE/100,limsupE/100,len(time_vec))
|
|
|
e2 = e2 + np.random.uniform(liminfE/100,limsupE/100,len(time_vec))
|
|
|
|
|
|
plt.figure(1)
|
|
|
plt.plot(b1[0:len(time_vec)/5000],'r')
|
|
|
plt.plot(b2[0:len(time_vec)/5000],'g')
|
|
|
plt.plot(b3[0:len(time_vec)/5000],'b')
|
|
|
plt.show()
|
|
|
|
|
|
#we want the input signal to be quantized (integer values)
|
|
|
for i in range(len(time_vec)):
|
|
|
b1[i] = quant16(b1[i],liminfB,limsupB)
|
|
|
b2[i] = quant16(b2[i],liminfB,limsupB)
|
|
|
b3[i] = quant16(b3[i],liminfB,limsupB)
|
|
|
e1[i] = quant16(e1[i],liminfE,limsupE)
|
|
|
e2[i] = quant16(e2[i],liminfE,limsupE)
|
|
|
|
|
|
#we filter the signal with a digital filter
|
|
|
|
|
|
#Filter 0
|
|
|
a0_0 = -128
|
|
|
a0_1 = 189
|
|
|
a0_2 = -111
|
|
|
b0_0 = 58
|
|
|
b0_1 = -66
|
|
|
b0_2 = 58
|
|
|
|
|
|
CF0 = Filter_Coefficients(a0_0, a0_1, a0_2, b0_0, b0_1, b0_2)
|
|
|
F0 = Filter(CF0)
|
|
|
|
|
|
#Filter 1
|
|
|
a1_0 = -128
|
|
|
a1_1 = 162
|
|
|
a1_2 = -81
|
|
|
b1_0 = 58
|
|
|
b1_1 = -57
|
|
|
b1_2 = 58
|
|
|
|
|
|
CF1 = Filter_Coefficients(a1_0, a1_1, a1_2, b1_0, b1_1, b1_2)
|
|
|
F1 = Filter(CF1)
|
|
|
|
|
|
#Filter 2
|
|
|
a2_0 = -128
|
|
|
a2_1 = 136
|
|
|
a2_2 = -55
|
|
|
b2_0 = 29
|
|
|
b2_1 = -17
|
|
|
b2_2 = 29
|
|
|
|
|
|
CF2 = Filter_Coefficients(a2_0, a2_1, a2_2, b2_0, b2_1, b2_2)
|
|
|
F2 = Filter(CF2)
|
|
|
|
|
|
#Filter 3
|
|
|
a3_0 = -128
|
|
|
a3_1 = 114
|
|
|
a3_2 = -33
|
|
|
b3_0 = 15
|
|
|
b3_1 = 4
|
|
|
b3_2 = 15
|
|
|
|
|
|
CF3 = Filter_Coefficients(a3_0, a3_1, a3_2, b3_0, b3_1, b3_2)
|
|
|
F3 = Filter(CF3)
|
|
|
|
|
|
#Filter 4
|
|
|
a4_0 = -128
|
|
|
a4_1 = 100
|
|
|
a4_2 = -20
|
|
|
b4_0 = 15
|
|
|
b4_1 = 24
|
|
|
b4_2 = 15
|
|
|
|
|
|
CF4 = Filter_Coefficients(a4_0, a4_1, a4_2, b4_0, b4_1, b4_2)
|
|
|
F4 = Filter(CF4)
|
|
|
|
|
|
Filters = [F0,F1,F2,F3,F4]
|
|
|
F_10k = FilterGlobal(Filters)
|
|
|
|
|
|
#this filtering assure us of not having aliasing problems we downsampling to fm = 24576 Hz
|
|
|
b1_filtered = F_10k.convolve(b1)
|
|
|
b2_filtered = F_10k.convolve(b2)
|
|
|
b3_filtered = F_10k.convolve(b3)
|
|
|
e1_filtered = F_10k.convolve(e1)
|
|
|
e2_filtered = F_10k.convolve(e2)
|
|
|
|
|
|
#we downsample the signal to get to a sampling frequency fm = 24576Hz
|
|
|
b1_24576 = b1_filtered[::4]
|
|
|
b2_24576 = b2_filtered[::4]
|
|
|
b3_24576 = b3_filtered[::4]
|
|
|
e1_24576 = e1_filtered[::4]
|
|
|
e2_24576 = e2_filtered[::4]
|
|
|
|
|
|
S_24576, w_24576 = Spectral_MatriceAvgTime_Int(b1_24576,b2_24576,b3_24576,e1_24576,e2_24576,24576,256)
|
|
|
PB_24576 = PSD_B_Int(S_24576)
|
|
|
PE_24576 = PSD_E_Int(S_24576)
|
|
|
|
|
|
# #Here we're doing simply a downsample, without filtering before. We should pass by a CIC filter instead.
|
|
|
# b1_4096 = b1_24576[::6]
|
|
|
# b2_4096 = b2_24576[::6]
|
|
|
# b3_4096 = b3_24576[::6]
|
|
|
# e1_4096 = e1_24576[::6]
|
|
|
# e2_4096 = e2_24576[::6]
|
|
|
#
|
|
|
# #Here we're doing simply a downsample, without filtering before. We should pass by a CIC filter instead.
|
|
|
# b1_256 = b1_24576[::96]
|
|
|
# b2_256 = b2_24576[::96]
|
|
|
# b3_256 = b3_24576[::96]
|
|
|
# e1_256 = e1_24576[::96]
|
|
|
# e2_256 = e2_24576[::96]
|
|
|
#
|
|
|
# #Here we're doing simply a downsample, without filtering before. We should pass by a CIC filter instead.
|
|
|
# b1_16 = b1_4096[::256]
|
|
|
# b2_16 = b2_4096[::256]
|
|
|
# b3_16 = b3_4096[::256]
|
|
|
# e1_16 = e1_4096[::256]
|
|
|
# e2_16 = e2_4096[::256]
|
|
|
|
