[PHYSIQUE] TP Filtrage numérique

SPE/PHYSIQUE/TP_FILTRAGE_NUMERIQUE/main.py

@@ -0,0 +1,197 @@
+# -*- coding: utf-8 -*-
+"""
+action d'un filtre passe-bas d'ordre 1 sur un créneau
+avec un contrôle externe sur la pulsation de coupure.
+"""
+
+import numpy as np
+import scipy.integrate as integr
+from math import *
+import cmath
+import matplotlib.pyplot as plt
+from matplotlib.widgets import Slider,Button
+###################
+# ZONE A MODIFIER #
+###################
+
+Nb = 73
+T = 1
+w = 2*pi/T
+w0 = 1.0*w
+Q = .5
+Te = 0.01
+D = 4.56 * T
+Nb = int(D//Te)
+
+def carre(t):
+    if t-t//T < T/2:
+        return 1.
+    else:
+        return 0.
+
+def triangle(t):
+    if t%T < T/2:
+        return 2*(t%(T/2))
+    elif t%T > T/2:
+        return 1 - 2*(t%(T/2))
+    else: return 1
+p=.3
+def PWM(t):
+    if t%T < T*p:
+        return 1.
+    else:
+        return 0.
+
+sinabs = lambda t : np.abs(np.sin(w*t))
+sin = lambda t : np.sin(w*t)
+fonc = sinabs
+##########################
+# FIN DE ZONE A MODIFIER #
+##########################   
+   
+    
+# échantillonage signal d'entrée
+xn = np.zeros(Nb)
+for k in range(0,Nb):
+    xn[k] = fonc(k*Te)
+
+
+###################
+# ZONE A MODIFIER #
+###################
+    
+# calcul du signal filtré avec linéarisation d'Euler p=1/Te.(1-z^-1)
+def passe_(w0_, Q_):
+    yn = np.zeros(len(xn))
+    for i in range(1, len(xn)):
+        yn[i] = 0
+    return yn
+
+def passe_bas(w0_, Q_):
+    yn = np.zeros(len(xn))
+    for i in range(1, len(xn)):
+        yn[i] = (yn[i-1]+Te*w0_*xn[i])/(1+Te*w0_)
+    return yn
+
+def passe_bas_bil(w0_, Q_):
+    yn = np.zeros(len(xn))
+    for i in range(1, len(xn)):
+        yn[i] = ((2/(Te*w0_) - 1)*yn[i-1] + xn[i] + xn[i-1])/(1+2/(Te*w0_))
+    return yn
+
+def passe_haut(w0_, Q_):
+    yn = np.zeros(len(xn))
+    for i in range(1, len(xn)):
+        yn[i] = (yn[i-1]+xn[i]-xn[i-1])/(1+Te*w0_)
+    return yn
+
+def passe_haut_bil(w0_, Q_):
+    yn = np.zeros(len(xn))
+    for i in range(1, len(xn)):
+        yn[i] = ((1-Te*w0_/2)*yn[i-1]+xn[i]-xn[i-1])/(1+Te*w0_/2)
+    return yn
+
+def passe_bande(w0_, Q_):
+    yn = np.zeros(len(xn))
+    for i in range(2, len(xn)):
+        yn[i] = ((1+2*Q_/(Te*w0_))*yn[i-1] - (Q_/(Te*w0_))*yn[i-2] + xn[i] - xn[i-1])/(1 + Q_*Te*w0_ + Q_/(Te*w0_))
+    return yn
+
+def passe_bande_bil(w0_, Q_):
+    yn = np.zeros(len(xn))
+    for i in range(2, len(xn)):
+        yn[i] = ((4*Q_/(Te*w0_) -Q_*Te*w0_)*yn[i-1] + (1 - 2*Q_/(Te*w0_) - Q_*Te*w0_*.5)*yn[i-2] + xn[i] - xn[i-2])/(1 + Q_*Te*w0_*.5 + 2*Q_/(Te*w0_))
+    return yn
+
+def passe_passe_bande_bil(w0_, Q_):
+    yn = np.zeros(len(passe_bande_bil(w0_, Q_)))
+    for i in range(2, len(xn)):
+        yn[i] = ((4*Q_/(Te*w0_) -Q_*Te*w0_)*yn[i-1] + (1 - 2*Q_/(Te*w0_) - Q_*Te*w0_*.5)*yn[i-2] + xn[i] - xn[i-2])/(1 + Q_*Te*w0_*.5 + 2*Q_/(Te*w0_))
+    return yn
+Fonction_passe=passe_bande_bil
+    
+##########################
+# FIN DE ZONE A MODIFIER #
+##########################  
+
+
+#calcul de la FFT et troncage à la moitié pour ne pas visualiser le spectre miroir
+def fast_fourier(z):
+    zf = abs(np.fft.fft(z))
+    return(zf[0:Nb/2])
+
+#spectre du signal d'entrée
+xft = fast_fourier(xn)
+    
+# échelle des temps
+tt = np.linspace(0,D,Nb)
+
+#filtrage passe-bas
+yn =Fonction_passe(w0, Q)
+
+#spectre du signal de sortie
+yft = fast_fourier(yn)
+
+# échelle des fréquences
+f = np.linspace(0,1/(2*Te),Nb//2)
+
+fig,ax=plt.subplots()
+#plt.subplots()
+plt.subplots_adjust(bottom=0.25,wspace=0.3)
+
+
+plt.subplot(2,2,1)
+plt.title(u'''signal d'entrée échantillonné''')
+#dessin du signal échantillonné
+plt.plot(tt,xn,'bo')
+#plt.ylim(-0.2,1.2)
+
+plt.subplot(2,2,2)
+plt.title(u'réponse')
+#dessin du la réponse
+plt.grid()
+lrep,=plt.plot(tt,yn,'k--')
+
+
+plt.subplot(2,2,3)
+plt.title(u'''spectre du signal d'entrée''')
+#dessin du spectre du signal échantillonné
+plt.plot(f,xft)
+
+plt.subplot(2,2,4)
+plt.title(u'''spectre du signal de sortie''')
+#dessin du spectre du signal échantillonné
+lf,=plt.plot(f,yft,'r')
+
+plt.show()       
+
+#boîte contenant le contrôle dynamique
+axcolor='yellow'
+axw0=plt.axes([0.15,0.1,0.6,0.03],axisbg=axcolor)
+
+sw0=Slider(axw0,'$\omega_0/\omega$',0.01,20,valinit=w0/w)
+#action lorsque le contrôle dynamique est modifié
+#attention à changer les appels à la fonction si on change le nom passe_bas
+
+def update(val):
+    global w0
+    global Q
+    w0=w*sw0.val
+    lrep.set_ydata(Fonction_passe(w0, Q))
+    lf.set_ydata(fast_fourier(Fonction_passe(w0, Q)))
+    fig.canvas.draw_idle()
+
+sw0.on_changed(update)
+
+axQ = plt.axes([0.15,0.1-.05,0.6,0.03],axisbg=axcolor)
+sQ = Slider(axQ,'$Q$',0.001,5,valinit=.5)
+
+def update_Q(val):
+    global w0
+    global Q
+    Q = val
+    lrep.set_ydata(Fonction_passe(w0, Q))
+    lf.set_ydata(fast_fourier(Fonction_passe(w0, Q)))
+    fig.canvas.draw_idle()
+
+sQ.on_changed(update_Q)