processing_chain.py
178 lines
| 5.3 KiB
| text/x-python
|
PythonLexer
/ SRC / processing_chain.py
Alexis Jeandet
|
r0 | #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] | ||||