import numpy as np
from numpy import pi 
from scipy.integrate import *
import matplotlib.pyplot as plt

#paramètres des filtres
H0=1
E0=1
fc=1000
fr=1000
omegac=2*pi*fc
Ne=500 #nombre d'échantillons

def passeBas(f,tf,choixAffichage):
    plt.figure()
    plt.title("Passe-bas, f="+str(f))
    T=np.linspace(0,tf,2000)
    #entrée
    def entree(t):
        return E0*np.cos(2*pi*f*t)
    if choixAffichage[3]:
        plt.plot(T,entree(T),"-",label="e(t)")
    #réponse du filtre analogique (équation différentielle)
    def ED(s,t):
        return omegac*(H0*entree(t)-s)        
    if choixAffichage[1] :
        repC=odeint(ED,E0,T)
        plt.plot(T,repC,"-",label="s(t) [ana]")  
    # réponse du filtre numérique (équation discrète)
    if choixAffichage[0]:
        #échantillonnage des temps
        Tech=np.linspace(0,tf,Ne)
        fe=Ne/tf 
        repD=np.zeros(Ne)
        repD[0]=E0
        alpha=fe/(fe+omegac)
        for i in range(1,Ne):
            repD[i]=alpha*repD[i-1]+(1-alpha)*H0*entree(Tech[i])
        plt.plot(Tech,repD,"o",label="s(t) [num]")
    #régime établi : fonction de transfert
    if choixAffichage[2]:
        etab=np.real(H0*E0*np.exp(2*pi*f*T*1j)/(1+1j*f/fc))
        plt.plot(T,etab,"-",label="établi")
    plt.legend()
    plt.show()

def passeHaut(f,tf,choixAffichage,tronque=False):
    plt.figure()
    plt.title("Passe-haut, f="+str(f))
    T=np.linspace(0,tf,2000)
    #entrée
    def entree(t):
        return E0*np.cos(2*pi*f*t)
    if choixAffichage[3]:
        plt.plot(T,entree(T),"-",label="e(t)")
    #réponse du filtre analogique (équation différentielle)
    def ED(s,t):
        return -H0*E0*2*pi*f*np.sin(2*pi*f*t)-s*omegac            
    if choixAffichage[1]:
        repC=odeint(ED,E0,T)
        plt.plot(T,repC,"-",label="s(t) [ana]")  
    # réponse du filtre numérique (équation discrète)
    if choixAffichage[0]:
        #échantillonnage des temps
        Tech=np.linspace(0,tf,Ne)
        fe=Ne/tf
        repD=np.zeros(Ne)
        repD[0]=E0
        alpha=fe/(fe+omegac)
        for i in range(1,Ne):
            repD[i]=alpha*(repD[i-1]+H0*(entree(Tech[i])-entree(Tech[i-1])))
        plt.plot(Tech,repD,"o",label="s(t) [num]")
    #régime établi : fonction de transfert
    etab=np.real(H0*E0*1j*f/fc*np.exp(2*pi*f*T*1j)/(1+1j*f/fc))
    if tronque:
        Y=5*np.max(etab)
        plt.ylim(-Y,Y)
    if choixAffichage[2]:
        plt.plot(T,etab,"-",label="établi")
    plt.legend()
    plt.show()

def passeBande(Q,f,tf,choixAffichage,tronque=False):
    titre = "Passe-bande, f="+str(f)+", Q="+str(Q)
    plt.figure()
    plt.title(titre)
    T=np.linspace(0,tf,2000)
    #entrée
    def entree(t):
        return E0*np.cos(2*pi*f*t)
    if choixAffichage[3]:
        plt.plot(T,entree(T),"-",label="e(t)")
    #réponse du filtre analogique (équation différentielle)
    def ED(S,t):
        return np.array([S[1],-2*pi*fr/Q*(H0*E0*2*pi*f*np.sin(2*pi*f*t)+S[1]+Q*2*pi*fr*S[0])])        
    # réponse du filtre numérique (équation discrète)
    if choixAffichage[1] :
        Rep=odeint(ED,[E0,0], T)
        repC=Rep[:,0]
        plt.plot(T,repC,"-",label="s(t) [ana]") 
    if choixAffichage[0]:
        #échantillonnage des temps
        Tech=np.linspace(0,tf,Ne)
        fe=Ne/tf 
        repD=np.zeros(Ne)
        repD[0]=E0
        repD[1]=E0
        alpha=2*pi*fr/fe
        for i in range(2,Ne):
            w=H0*alpha/Q*(entree(Tech[i-1])-entree(Tech[i-2]))
            w+=(2-alpha/Q-alpha**2)*repD[i-1]
            w-=(1-alpha/Q)*repD[i-2]
            repD[i]=w
        plt.plot(Tech,repD,"o",label="s(t) [num]")
    #régime établi : fonction de transfert
    etab=np.real(H0*E0*np.exp(2*pi*f*T*1j)/(1+1j*Q*(f/fr-fr/f)))
    if tronque:
        Y=5*np.max(etab)
        plt.ylim(-Y,Y)
    if choixAffichage[2]:
        plt.plot(T,etab,"-",label="établi")
    plt.legend()
    plt.show()   

#choisir le type d'affichage
afficheEntree=True
solAna=True
solNum=True
regimeEtabli=True
choixAffichage=[solNum,solAna,regimeEtabli,afficheEntree]

#Choisir une ligne à décommenter
#NB : pour la fréquence 10 il y a souvent à faire attention #au choix de la fenêtre d'affichage et à l'échantillonnage
#passeBas(10,0.1,choixAffichage)
#zoomons car peu de différence transitoire établi
#passeBas(10,0.001,choixAffichage)
#passeBas(1000,0.005,choixAffichage)
#passeBas(10000,0.001,choixAffichage)

#passeHaut(10,0.1,choixAffichage)
#zoomons, peu de différence transitoire établi
#passeHaut(10,0.1,choixAffichage,tronque=True)
#passeHaut(1000,0.001,choixAffichage)
#passeHaut(10000,0.00015,choixAffichage)


#pour les deux premiers (fréquence 10), 
#la solution discrète ne marche que si on augmente Ne
passeBande(10,10,0.1,choixAffichage)  
passeBande(10,10,0.1,choixAffichage,tronque=True)
passeBande(10,1000,0.001,choixAffichage)
passeBande(10,10000,0.01,choixAffichage)
passeBande(1,10000,0.002,choixAffichage)