import numpy
from matplotlib.pyplot import *
import numpy.fft
import scipy.signal

[t,u] = numpy.loadtxt("sinus100Hz-fe10kHz-8bits.txt2")
Ne=len(t)
te=t[1]-t[0]
fe = 1.0/te
tfd = numpy.fft.fft(u*u)
freq = numpy.arange(Ne)*1.0/(Ne*te)
figure()
plot(freq,10*numpy.log10(numpy.absolute(tfd)))
xlabel("f (Hz)")
ylabel("dB")
title("8 bits")
grid()
             

[t,u] = numpy.loadtxt("sinus100Hz-fe10kHz-10bits.txt2")
tfd = numpy.fft.fft(u*u)
freq = numpy.arange(Ne)*1.0/(Ne*te)
figure()
plot(freq,10*numpy.log10(numpy.absolute(tfd)))
xlabel("f (Hz)")
ylabel("dB")
title("10 bits")
grid()
             

[t,u] = numpy.loadtxt("sinus100Hz-fe10kHz-12bits.txt2")
tfd = numpy.fft.fft(u*u)
freq = numpy.arange(Ne)*1.0/(Ne*te)
figure()
plot(freq,10*numpy.log10(numpy.absolute(tfd)))
xlabel("f (Hz)")
ylabel("dB")
title("12 bits")
grid()
             


[t,u] = numpy.loadtxt("sinus100Hz-fe100kHz-8bits.txt2")
Ne=len(t)
te=t[1]-t[0]
fe = 1.0/te
tfd = numpy.fft.fft(u*u)
freq = numpy.arange(Ne)*1.0/(Ne*te)
figure()
plot(freq,10*numpy.log10(numpy.absolute(tfd)))
xlabel("f (Hz)")
ylabel("dB")
title("8 bits")
grid()
             

fc = 4.0e3
b = scipy.signal.firwin(numtaps=100,cutoff=[fc/fe],window='hann',nyq=0.5)
[w,h] = scipy.signal.freqz(b,[1])
figure()
plot(w/(2*numpy.pi),20*numpy.log10(numpy.abs(h)),'r')
xlabel('f/fe')
ylabel('GdB')
grid()
             

u1 = scipy.signal.lfilter(b,[1],u)
tfd = numpy.fft.fft(u1*u1)
freq = numpy.arange(Ne)*1.0/(Ne*te)
figure()
plot(freq,10*numpy.log10(numpy.absolute(tfd)))
xlabel("f (Hz)")
ylabel("dB")
title("8 bits")
grid()
             

u2 = u1[0::10]
fe2 = fe/10
Ne2 = len(u2)
tfd = numpy.fft.fft(u2*u2)
freq = numpy.arange(Ne2)*1.0/(Ne2/fe2)
figure()
plot(freq,10*numpy.log10(numpy.absolute(tfd)))
xlabel("f (Hz)")
ylabel("dB")
title("8 bits 100kHz+filtrage+reduction")
grid()
             

def derive(x,te):
    a=[te]
    b=[1,-1]
    return scipy.signal.lfilter(b,a,x)
            

[t,x] =  numpy.loadtxt("sinus100Hz-fe10kHz-10bits.txt2")
y = derive(x,t[1]-t[0])
figure()
plot(t,x,label="x")
plot(t,y/1000,label="dx/dt")
xlabel("t (s)")
grid()
axis([0.1,0.15,-5,5])
            

[t,x] =  numpy.loadtxt("sinus100Hz-fe100kHz-10bits.txt2")
y = derive(x,t[1]-t[0])
figure()
plot(t,x,label="x")
plot(t,y/1000,label="dx/dt")
xlabel("t (s)")
grid()
axis([0.1,0.15,-5,5])
            

fc = 4.0e3
fe = 1e5
b = scipy.signal.firwin(numtaps=100,cutoff=[fc/fe],window='hann',nyq=0.5)
y = scipy.signal.lfilter(b,[1],x)
y1 = derive(y,t[1]-t[0])
figure()
plot(t,x,'b',label="x")
plot(t,y1/1000,'r',label="dx/dt")
xlabel("t (s)")
grid()
axis([0.1,0.15,-5,5])
             

y2 = y[0::10]
t2 = t[0::10]
y3 = derive(y2,t2[1]-t2[0])
figure()
plot(t,x,'b',label="x")
plot(t2,y3/1000,'g',label="dx/dt")
xlabel("t (s)")
grid()
axis([0.1,0.15,-5,5])
             

figure()
plot(t,y1/1000,'r',label="sans reduction")
plot(t2,y3/1000,'g',label="avec reduction")
xlabel("t (s)")
legend(loc="upper right")
grid()
axis([0.1,0.103,1,1.6])
              

[t,x] =  numpy.loadtxt("sinus100Hz-fe10kHz-10bits.txt2")
fc = 4.0e3
fe = 1e4
b = scipy.signal.firwin(numtaps=100,cutoff=[fc/fe],window='hann',nyq=0.5)
y = scipy.signal.lfilter(b,[1],x)
y1 = derive(y,t[1]-t[0])
figure()
plot(t,y1/1000,'r',label="dx/dt")
xlabel("t (s)")
grid()
axis([0.106,0.11,1,1.6])