import numpy
from matplotlib.pyplot import *

data = numpy.loadtxt("lampeNa-7-10kHz.txt")
t = data[0]
u = data[1]
te = t[1]-t[0]
fe = 1.0/te
N = t.size

figure()
plot(t,u)
xlabel("t (s)")
ylabel("u (V)")
grid()
axis([0,0.1,0,1])
                 

import numpy.fft


spectre = numpy.absolute(numpy.fft.fft(u))/N
f = numpy.arange(N)*fe/N
figure()
plot(f,spectre,'r')
xlabel("f (Hz)")
ylabel("S")
grid()
axis([-100,2000,0,spectre.max()])
                  

import scipy.signal

fe = 1000.0
fc = 20.0
P = 200
h = scipy.signal.firwin(numtaps=P,cutoff=[fc/fe],nyq=0.5,window='hann')
w,H = scipy.signal.freqz(h)
figure()
plot(w/(2*numpy.pi)*fe,numpy.absolute(H))
xlabel("f (Hz)")
ylabel("G")
grid()
axis([0,100,0,1.2])
                     

figure()
plot(w/(2*numpy.pi)*fe,numpy.unwrap(numpy.angle(H)))
xlabel("f (Hz)")
ylabel("phi")
grid()
axis([0,100,-20,0])
                      

fc = 5.0
def H(f):
    return 1.0/(1.0+numpy.sqrt(2)*1j*f/fc-(f/fc)**2)
f = numpy.linspace(0,100,1000)
Ht = H(f)
figure()
subplot(211)
plot(f,numpy.absolute(Ht))
xlabel("f (Hz)")
ylabel("G")
grid()
subplot(212)
plot(f,numpy.angle(Ht))
xlabel("f (Hz)")
ylabel("phi")
grid()
                     

f = numpy.linspace(0,5,200)
Ht = H(f)
figure()
subplot(211)
plot(f,numpy.absolute(Ht))
xlabel("f (Hz)")
ylabel("G")
grid()
subplot(212)
plot(f,numpy.angle(Ht))
xlabel("f (Hz)")
ylabel("phi")
grid()
                     

data = numpy.loadtxt("Na-8.txt")
t = data[0]
u = data[1]
N = t.size
delta = t*10.0/9
figure()
plot(delta,u)
xlabel("delta (microns)")
ylabel("u (V)")
axis([0,delta.max(),0,0.6])
grid()
                 

axis([320,360,0,0.6])
                 

v = scipy.signal.hamming(N)*(u-u.mean())
spectre = numpy.absolute(numpy.fft.fft(v))/N
d = delta[1]-delta[0]
sigma = numpy.arange(N)*1.0/(N*d)
figure()
plot(sigma,spectre)
xlabel("sigma (1/microns)")
ylabel("S")
Smax = spectre.max()
axis([0,3,0,Smax])
grid()
                 

figure()
specgram(v,NFFT=5000,Fs=1.0/d,noverlap=0)
ylabel("sigma (1/microns)")
axis([0,1000,0,3])
                      

data = numpy.loadtxt("Na-9.txt")
t = data[0]
u = data[1]
N = t.size
delta = t*10.0/9
figure()
plot(delta,u)
xlabel("delta (microns)")
ylabel("u (V)")
axis([0,delta.max(),0,1.0])
grid()
                 

u=u-u.mean()
v = scipy.signal.hamming(N)*u
spectre = numpy.absolute(numpy.fft.fft(v))/N
d = delta[1]-delta[0]
sigma = numpy.arange(N)*1.0/(N*d)
figure()
plot(sigma,spectre)
xlabel("sigma (1/microns)")
ylabel("S")
Smax=spectre.max()
axis([0,3,0,Smax])
grid()
                 

N1 = int((2500.0-100.0)/4000*N)
N2 = int((2500.0+100.0)/4000*N)
v = scipy.signal.hamming(N2-N1)*u[N1:N2]
spectre = numpy.absolute(numpy.fft.fft(v))/(N2-N1)
sigma = numpy.arange(N2-N1)*1.0/((N2-N1)*d)
figure()
plot(sigma,spectre)
xlabel("sigma (1/microns)")
ylabel("S")
title("Sodium")
Smax = spectre.max()
axis([0,3,0,Smax]) 
grid()
                  

data = numpy.loadtxt("Na-fast-7-filtre.txt")
t = data[0]
u = data[1]
N = t.size
delta = t*10.0/9
figure()
plot(delta,u)
xlabel("delta (microns)")
ylabel("u (V)")
axis([400,420,0,1.0])
grid()
                        

delta = t*10.0/9
v = scipy.signal.hamming(N)*(u-u.mean())
spectre = numpy.absolute(numpy.fft.fft(v))/N
d = delta[1]-delta[0]
sigma = numpy.arange(N)*1.0/(N*d)
figure()
plot(sigma,spectre)
xlabel("sigma (1/microns)")
ylabel("S")
axis([0,3,0,0.005])
grid()
                 

data = numpy.loadtxt("Hg-5.txt")
t = data[0]
u = data[1]
N = t.size
delta = t*10.0/9
figure()
plot(delta,u)
xlabel("delta (microns)")
ylabel("u (V)")
axis([0,delta.max(),0,2])
grid()
                 

u = u-u.mean()
v = scipy.signal.hamming(N)*u
spectre = numpy.absolute(numpy.fft.fft(v))/N
d = delta[1]-delta[0]
sigma = numpy.arange(N)*1.0/(N*d)
figure()
plot(sigma,spectre)
xlabel("sigma (1/microns)")
ylabel("S")
Smax = spectre.max()
axis([0,3,0,Smax])
grid()
                 

d_max = delta[N-1]
N1 = int((1780.0-100.0)/d_max*N)
N2 = int((1780.0+100)/d_max*N)
v = scipy.signal.hamming(N2-N1)*u[N1:N2]
spectre = numpy.absolute(numpy.fft.fft(v))/(N2-N1)
sigma = numpy.arange(N2-N1)*1.0/((N2-N1)*d)
figure()
plot(sigma,spectre)
xlabel("sigma (1/microns)")
ylabel("S")
title("Mercure")
Smax = spectre.max()
axis([0,3,0,Smax]) 
grid()
                  

figure()
plot(1/sigma,spectre)
xlabel("lambda (microns)")
ylabel("S")
grid()
axis([0.2,1.5,0,Smax])
                  

L = [0.365,0.405,0.436,0.546,0.577,0.579]
for l in L:
    plot([l,l],[0,Smax],'r')
axis([0.2,0.8,0,Smax])
                   

data = numpy.loadtxt("Hg-Cd-3.txt")
t = data[0]
u = data[1]
N = t.size
delta = t*10.0/9
figure()
plot(delta,u)
xlabel("delta (microns)")
ylabel("u (V)")
axis([0,delta.max(),0,2])
grid()
                 

u=u-u.mean()
d_max = delta[N-1]
N1 = int((325.0-100.0)/d_max*N)
N2 = int((325.0+100)/d_max*N)
v = scipy.signal.hamming(N2-N1)*u[N1:N2]
spectre = numpy.absolute(numpy.fft.fft(v))/(N2-N1)
sigma = numpy.arange(N2-N1)*1.0/((N2-N1)*d)
figure()
plot(1.0/sigma,spectre)
xlabel("lambda (microns)")
ylabel("S")
title("Mercure-Cadmium")
Smax = spectre.max()
axis([0.2,0.8,0,Smax]) 
grid()
                  

data = numpy.loadtxt("incandescence-QI-2.txt")
t = data[0]
u = data[1]
N = t.size
delta = t*10.0/9
figure()
plot(delta,u)
xlabel("delta (microns)")
ylabel("u (V)")
axis([40,65,0,5.0])
grid()
                 

u = u-u.mean()
spectre = numpy.absolute(numpy.fft.fft(u))/N
d = delta[1]-delta[0]
sigma = numpy.arange(N)*1.0/(N*d)
figure()
plot(1/sigma,spectre)
xlabel("lambda (microns)")
ylabel("S")
title("Incandescence Quartz-Iode")
axis([0.2,1.5,0,2e-2])
grid()
                 

def analyse(fichier):
    data = numpy.loadtxt(fichier)
    t = data[0]
    u = data[1]
    u = u-u.mean()
    N = t.size
    print(N)
    delta = t*10.0/9
    figure()
    subplot(211)
    plot(delta,u)
    xlabel("delta (microns)")
    ylabel("u (V)")
    grid()
    v = scipy.signal.hamming(N)*u
    spectre = numpy.absolute(numpy.fft.fft(v))/N
    d = delta[1]-delta[0]
    sigma = numpy.arange(N)*1.0/(N*d)
    subplot(212)
    plot(1/sigma,spectre)
    xlabel("lambda (microns)")
    ylabel("S")
    Smax = spectre.max()
    axis([0.4,1.2,0,Smax])
    grid()
    return (sigma,spectre)
                  

(sigma,S) = analyse("incandescence-QI-4.txt")
                  

(sigma_rouge,S_rouge) = analyse("incandescence-QI-filtreRouge-2.txt")
                  

(sigma_jaune,S_jaune) = analyse("incandescence-QI-filtreJaune-1.txt")
                  

(sigma_cyan,S_cyan) = analyse("incandescence-QI-filtreCyan-1.txt")
                  

A_rouge = 1.0-S_rouge/S
A_jaune = 1.0-S_jaune/S
A_cyan = 1.0-S_cyan/S
figure()
plot(1/sigma_rouge,A_rouge,'r')
plot(1/sigma_jaune,A_jaune,'g')
plot(1/sigma_cyan,A_cyan,'b')
xlabel("lambda (microns)")
ylabel("A")
grid()
axis([0.4,1.2,0,1])