import numpy
from matplotlib.pyplot import *

[f,Ue,Us,phi]=numpy.loadtxt("reponseFreq-G30.txt",unpack=True,skiprows=1)
G=Us/Ue
GdB=20*numpy.log10(G)
figure(figsize=(8,8))
subplot(211)
plot(f,GdB,'-o')
xscale('log')
ylabel('$G_{dB}$')
grid()
subplot(212)
plot(f,phi,'-o')
xscale('log')
ylabel('$\phi$')
xlabel('f (Hz)')
grid()

                 

figure(figsize=(8,6))
plot(f,G)
gc=G.max()/numpy.sqrt(2)
xlim(0,2000)
plot([100,2000],[gc,gc],'k--')
xlabel("f (Hz)")
ylabel("G")
grid()
                  

[t,ue,us]=numpy.loadtxt("reponseEchelon-SP5-G30.txt",unpack=True,skiprows=1)
figure()
xlabel("t (s)")
ylabel("Tensions (V)")
plot(t,ue,label="e(t)")
plot(t,us,label="s(t)")
grid()
xlim(0,0.02)
ylim(-4,4)
legend(loc="upper right",fontsize=20)
                     

[t,ue,us]=numpy.loadtxt("reponseImp-G30.txt",unpack=True,skiprows=2)

figure()
xlabel("t (s)",fontsize=18)
ylabel("Tensions (V)",fontsize=18)
plot(t,ue,label="e(t)")
plot(t,us,label="s(t)")
grid()
xlim(0,0.02)
ylim(-11,11)
legend(loc="upper right",fontsize=20)
                

from numpy.fft import fft
N=len(ue)
T=t[N-1]
tfd_e = fft(ue)*2/N
spectre_e = numpy.absolute(tfd_e)
tfd_s = fft(us)*2/N
freq = numpy.arange(N)*1.0/T
spectre_s = numpy.absolute(tfd_s)

figure()
plot(freq,spectre_e,label="E")
plot(freq,spectre_s,label="S")
xlabel("f (Hz)",fontsize=18)
ylabel("Tensions (V)",fontsize=18)
grid()
legend(loc="upper right",fontsize=20)
axis([0,25000,0,spectre_s.max()])
                

gain=spectre_s/spectre_e
figure()
plot(freq,20*numpy.log10(gain))
plot(f,GdB,"k--")
xlabel("f (Hz)",fontsize=18)
ylabel("dB",fontsize=18)
xscale('log')
grid()
xlim(1e2,1e4)
ylim(-20,40)
                 

import scipy.signal
fe=20000.0
fo=1000.0
omega_a=numpy.pi*fo/(fe/2)
r=0.98
b0=(1-r*r)*0.5
b=[b0,0,-b0]
a=[1,-2*r*numpy.cos(omega_a),r*r]
w,h =scipy.signal.freqz(b,a)
g = numpy.absolute(h)
phase = numpy.unwrap(numpy.angle(h))
figure(figsize=(8,8))
subplot(211)
plot(w/(2*numpy.pi)*fe,g)
xlabel("f (Hz)")
ylabel("G")
grid()
subplot(212)
plot(w/(2*numpy.pi)*fe,phase)
xlabel("f (Hz)")
ylabel("phi")
grid()
                   

figure()
plot(w/(2*numpy.pi)*fe,20*numpy.log10(g))
xlabel('f (Hz)')
ylabel("GdB")
xscale('log')
axis([100,10000,-60,0])
grid()
                   

N=400
t=numpy.arange(N)*1.0/fe
y=numpy.zeros(N)
x=numpy.zeros(N)
x[0]=1
y[0]=b[0]*x[0]
y[1]=b[0]*x[1]+b[1]*x[0]-a[1]*y[0]
y[2]=b[0]*x[2]+b[1]*x[1]+b[2]*x[0]-a[1]*y[1]-a[2]*y[0]
for n in range(3,N):
    y[n] = b[0]*x[n]+b[1]*x[n-1]+b[0]*x[n-2]-a[1]*y[n-1]-a[2]*y[n-2]

figure()
plot(t,y)
xlabel("t (s)")
ylabel("h")
grid()