## Paramètres importants
import numpy as np
import matplotlib.pyplot as plt 
import scipy.integrate as si 

f=12 #fréquence du signal
T=1/f #période du signal 
w=2*np.pi/T #pulsation du signal

fe=15 #fréquence d'échantillonnage
Te=1/fe

print(f"La fréquence du signal est f={f:.2f} Hz")   # écriture formatée avec le print(f"...")
print(f"La fréquence d'échantillonnage est fe={fe:.2f} Hz")  # ici variable {var:.2f} avec 2 décimales


## Echantillonnage d'un signal
 
# On cherche à modéliser l'échantillonnage d'un signal sinusoïdal connu de fréquence $f$ à une fréquence d'échantillonnage $f_e$ et à comparer le signal échantillonné au signal analogique initial.
# 
# On va chercher à observer $N$ périodes et on définit alors deux numpy_array temps qui contiennent respectivement 
# -les points d'échantillonnage pour t
# -un "continuum" (donc un grand nombre de points) pour t2

N=10 #nombre de période
t=np.linspace(0,((N*T)//Te)*Te,int((N*T)//Te)+1) #on génère les instants "échantillonnés"
t2=np.linspace(0,((N*T)//Te)*Te,N*200) #ainsi qu'une échelle de temps plus "continu" pour visualiser le signal analogique


# On rappelle que du fait du critère de Shannon, on ne peut observer correctement que les fréquences $f\in[0,f_e/2]$.
# Si cette condition n'est pas vérifiée, il y aura repliement de spectre et on aura une fréquence $fvue$ perçue dans l'intervalle qui sera erronnée. 

# On trace :
# - la fonction d'entrée analogique d'origine
# - les points relevés par échantillonnage
# - la sinusoïde perçue dans l'intervalle $[0,f_e/2]$ (il faut alors dans un premier temps calculer $\texttt{fvue}$
plt.close('all')

plt.plot(t2,np.cos(w*t2))
plt.plot(t,np.cos(w*t),linestyle='none',marker='o',markersize=12,color='g')
if (f%(fe)<(fe/2)):
    fvue=f%(fe)
else:
    fvue=fe-f%(fe)
print(f"La fréquence détectée (dans l'intervalle [0,fe/2] après d'échantillonnage est f={fvue:.2f} Hz")
plt.plot(t2,np.cos(2*np.pi*t2*fvue),color='g')
plt.show()

## Filtrage numérique
# 
# On cherche dans un second temps à modéliser l'effet d'un filtrage numérique sur un signal analogique de fréquence $f$ échantillonné.
# 
# Soit un passe-bas d'ordre 1 de fréquence de coupure $f_c$ et de gain $H_0=1$.
# La relation discrètre qui existe les valeurs échantillonnées s_n et celle d'entrée e_n est 
# s_{n+1}=s_n (1-T_e/{\tau})+T_e/{\tau} e_n

f=1 #fréquence du signal
T=1/f #période du signal 
w=2*np.pi/T #pulsation du signal

fe=4 #fréquence d'échantillonnage (on s'assure ici de respecter Shannon)
Te=1/fe #fe=4 =>problème  ; fe=100 OK

fc=1 #fréquence de coupure du filtre 
tau=1/(2*np.pi*fc) #temps caractéristique


N=10 #nombre de périodes qu'on souhaite visualiser 
tn=np.linspace(0,((N*T)//Te)*Te,int((N*T)//Te)+1) #on génère les instants "échantillonné"
tcontinu=np.linspace(0,((N*T)//Te)*Te,N*200) #ainsi qu'une échelle de temps plus "continu" pour visualiser le signal analogique

def e(t): #on définit la fonction d'entrée analogique qu'on va échantillonner
    return np.cos(w*t)


##Réponse à un signal sinusoïdal
# On suppose qu'on part de s(0)=0.

def s(t):
    phi=-np.arctan(f/fc)
    H=1/np.sqrt(1+(f/fc)**2)
    return H*np.cos(w*t+phi)

en=e(tn) #(on utilise ici la vectorisation de numpy)
sn=np.zeros(len(tn))
for i in range(1,len(tn)):
    sn[i]=Te/tau*(en[i-1]-sn[i-1])+sn[i-1]

plt.figure()
plt.plot(tn,sn,linestyle='none',marker='o',markersize=12,color='g')
plt.plot(tcontinu,s(tcontinu))
plt.show()

## Réponse à un créneau
def e2(t): #on teste cette fois avec un créneau
    if (t%T)<T/2:
        return 1
    else:
        return -1
e2=np.vectorize(e2)

def eqdiff(X,t):
    return (e2(t)-X)/tau

cini=[0]
s2_sol=si.odeint(eqdiff,cini,tcontinu)
e2n=e2(tn)
s2n=np.zeros(len(tn))

for i in range(1,len(tn)):
    s2n[i]=Te/tau*(e2n[i-1]-s2n[i-1])+s2n[i-1]

plt.figure()
plt.plot(tn,s2n,linestyle='none',marker='o',markersize=12,color='g')
plt.plot(tcontinu,s2_sol)
plt.show()