#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]