Synthèse de Filtres Analogiques

Ce tutorial montre comment synthétiser et analyser des filtres analogiques avec python.

Home Electronic .ipynb
In [561]:
%matplotlib inline
import numpy as np
import scipy.signal as sig
import matplotlib.pyplot as plt

En électronique, nous souhaitons réaliser des filtres avec un gabarit particulier. Ce gabarit est spécifié via plusieurs paramètres caractérisant le module du filtre $|H(j\omega)|$:

  • Le type de filtre: passe-bas, passe-haut, passe-bande, etc.
  • $f_c$ (Hz): la fréquence de coupure (cutoff frequency) du filtre définissant la bande passante (pass band).
  • $T_c$ (dB): le gain minimum dans la bande passante
  • $f_s$ (Hz): la première fréquence rejetée définissant le début de la bande atténuée (stop band)
  • $T_s$ (dB): le gain maximum dans la bande atténuée.

Concernant la phase, nous souhaitons dans l'idéal que cette dernière soit linéaire dans la bande passante (retard de groupe constant). La problématique de synthèse d'un filtre consiste alors à déterminer l'ordre $N$ ainsi que les coefficients $a_k$ et $b_k$ de la fonction de transfert de manière à respecter le gabarit du filtre. Dans les parties suivantes, nous allons présenter plusieurs structures de filtres très utilisées.

Pour illustrer les caractéristiques des ces structures, nous allons les utiliser pour synthétiser un filtre passe-bas d'ordre $N=3$ ayant les paramètres suivants:

In [562]:
N = 3          # ordre
fc = 100       # fréquence de coupure
wc = 2*np.pi*fc   #pulsation de courpure
Tc = 3         # Gain minimum dans la bande passante en dB
Ts = 40        # Gain maximum dans la bande attenuée en dB

# Base des pulsations et des fréquences pour l'analyse frequentielle
w = np.logspace(1,4,num=100)
f = w/(2*np.pi)
In [563]:
[b,a]=sig.butter(N,wc,btype="low",analog=True)
H_butter=sig.lti(b,a)

Module du filtre

Le module d'un filtre de Butterworth s'exprime sous la forme:

$$|H(j\omega)|=\frac{1}{\sqrt{1+(\frac{\omega}{\omega_c})^{2N}}}$$

où $\omega_c=2\pi f_c$ désigne la pulsation de coupure à $-3dB$ c-a-d $|H(j\omega_c)|=\frac{1}{\sqrt{2}}$. Le filtre de Butterworth possède la propriété d'avoir un module le plus plat que possible dans la bande passante (maximally flat). De plus, le module possède la propriété d'être monotone (sans ondulation) à la fois dans la bande passante et dans la bande atténuée.

In [564]:
[w, r_w] = H_butter.freqresp(w)
module = np.abs(r_w)

plt.loglog(f,module)
plt.axhline(1/np.sqrt(2), color='g')
plt.axvline(fc, color='g')
plt.ylabel("Module")
plt.xlabel("Fréquence (Hz)");

Pôles du filtre

Les pôles d'un filtre de Butterworth s'obtiennent en ne conservant que les solutions de l'équation $1+(s/j\omega_c)^{2N}=0$ situées sur la gauche du plan complexe. Mathématiquement, ces pôles s'expriment sous la forme ($k=0,\cdots,N-1$):

$$p_k=\omega_c e^{j\varphi_k}$$

avec $\varphi_k=\frac{\pi}{2N}(2k+N+1)$. Dans le cas présent, nous avons donc 3 pôles en $\omega_c e^{2j\pi/3}$, $-\omega_c$ et $\omega_c e^{4j\pi/3}$.

In [565]:
poles = H_butter.poles
zeros = H_butter.zeros

angle=np.arange(-np.pi,np.pi,0.01)
plt.plot(poles.real, poles.imag, 'x', markersize=5)
plt.plot(zeros.real, zeros.imag,  'o', markersize=5) 
plt.plot(wc*np.sin(angle), wc*np.cos(angle), 'g--') 
plt.axis('scaled')
plt.xlabel("Partie Réelle");
plt.ylabel("Partie Imaginaire");
plt.axis([-1.5*wc,1.5*wc,-1.5*wc,1.5*wc]);
In [566]:
[b,a]=sig.cheby1(N,Tc,wc,btype="low",analog=True)
H_cheby=sig.lti(b,a)

Module du filtre

Le module d'un filtre de Chebyshev s'exprime sous la forme:

$$|H(j\omega)|=\frac{1}{\sqrt{1+\epsilon^2 T_N^2\left(\frac{\omega}{\omega_c}\right)}}$$

où $\omega_c=2\pi f_c$ désigne la pulsation de coupure, $\epsilon$ est un paramètre du filtre lié au gain minimum dans la bande passante $T_c$ ($T_c=-10\log_{10} (1+\epsilon^2)$), et $T_N(x)$ correspond au polynôme de Chebyshev d'ordre $N$ défini par:

$$V_N(x)=\cos(N\cos^{-1}(x))$$

Le filtre de Chebyshev possède la propriété d'avoir un module avec ondulation (ripple) dans la bande passante et d'être monotone dans la bande atténuée. En laissant passer une ondulation dans la bande passante, le filtre de Chebyshev permet d'obtenir de meilleurs performances dans la bande attenuée que le filtre de Butterworth. En particulier, pour une fréquence de coupure et un gain minimum tolérée dans la bande passante identiques, le filtre de Chebyshev atténue plus fortement le signal dans la bande rejetée.

In [567]:
[w, r_w] = H_cheby.freqresp(w)
module = np.abs(r_w)

Tc_val_nat=10**(-Tc/20)
plt.loglog(f,module)
plt.axhline(Tc_val_nat, color='g')
plt.axvline(fc, color='g')
plt.ylabel("Module")
plt.xlabel("Fréquence (Hz)");

Pôles du filtre

Les pôles s'obtiennent en conservant les solutions de l'équation $1+\epsilon^2 T_N^2\left(\frac{\omega}{\omega_c}\right)=0$ situées sur la gauche du plan complexe. Il est possible de montrer que ces pôles appartiennent à une ellipse centrée à l'origine et ayant des axes mineur et majeur de longueur respective

$$a=\omega_c (\alpha^{1/N}-\alpha^{-1/N})$$$$b=\omega_c (\alpha^{1/N}+\alpha^{-1/N})$$

où $\alpha=\epsilon^{-1}+\sqrt{1+\epsilon^{-2}}$.

In [568]:
poles = H_cheby.poles
zeros = H_cheby.zeros

epsilon = np.sqrt(10**(0.1*Tc)-1)

alpha = epsilon**(-1)+np.sqrt(1+ (epsilon**-2))
a = wc*(alpha**(1/N)-alpha**(-1/N))
b = wc*(alpha**(1/N)+alpha**(-1/N))
angle=np.arange(-np.pi,np.pi,0.01)

plt.plot(poles.real, poles.imag, 'x', markersize=5)
plt.plot(zeros.real, zeros.imag,  'o', markersize=5) 
plt.plot(0.5*a*np.sin(angle),0.5*b*np.cos(angle),'g--')
plt.axis('scaled')
plt.xlabel("Partie Réelle");
plt.ylabel("Partie Imaginaire");
plt.axis([-1.5*wc,1.5*wc,-1.5*wc,1.5*wc]);
In [569]:
[b,a]=sig.bessel(3, wc, btype='low', norm='delay', analog=True)
H_bessel=sig.lti(b,a)

Module du filtre

Le module d'un filtre de Bessel s'exprime sous la forme:

$$|H(j\omega)|=\frac{K_0}{B_N\left(\frac{\omega}{\omega_c}\right)}$$

où $\omega_c=2\pi f_c$ désigne la fréquence de coupure du filtre, $K_0$ est un facteur de normalisation permettant d'obtenir un gain unitaire ($|H(0)|=1$), et $B_N(x)$ correspond à un polynôme de Bessel d'ordre $N$. Un polynôme de Bessel d'ordre $N$ peut être défini par la relation de récurrence suivante:

$$ \begin{align} B_0(x)&=1\\\ B_1(x)&=x+1\\\ B_N(x)&=(2N-1)B_{N-1}(x)+x^2B_{N-2}(x) \end{align} $$
In [570]:
[w, r_w] = H_bessel.freqresp(w)
module = np.abs(r_w)

plt.loglog(f,module)
plt.axvline(fc, color='g')
plt.ylabel("Module")
plt.xlabel("Fréquence (Hz)");

Retard de groupe

Le filtre de Bessel possède la propriété de maximiser la linéarité de la phase dans la bande passante. Nous obtenons alors un retard de groupe le plus plat que possible dans la même bande.

In [571]:
retard_groupe = -np.diff(np.unwrap(np.angle(r_w)))/np.diff(w)

plt.figure()
plt.semilogx(f[1:],retard_groupe)
plt.axvline(fc, color='g')
plt.ylabel("Retard groupe (s)")
plt.xlabel("Fréquence (Hz)");
In [572]:
[b,a]=sig.ellip(N, Tc, Ts, wc, btype='low', analog=True)
H_ellip=sig.lti(b,a)

Module du filtre

Le module d'un filtre Elliptique s'exprime sous la forme:

$$|H(j\omega)|=\frac{1}{\sqrt{1+\epsilon^2 U_N^2\left(\mu,\frac{\omega}{\omega_c}\right)}}$$

où $\omega_c=2\pi f_c$ désigne la pulsation de coupure, $\epsilon$ est un paramètre du filtre lié au gain minimum dans la bande passante, $U_N(x,y)$ est la fonction elliptique rationnelle d'ordre N et $\mu$ est un paramètre lié au gain maximum dans la bande atténuée.

Le filtre Elliptique possède la propriété d'avoir un module avec ondulation (ripple) à la fois dans la bande passante et dans la bande atténuée. De plus ce filtre possède la propriété de respecter le gabarit du filtre tout en minimisant la bande de transition c-a-d la distance entre $f_c$ et $f_s$.

In [573]:
[w, r_w] = H_ellip.freqresp(w)
module = np.abs(r_w)

Tc_val_nat=10**(-Tc/20)
Ts_val_nat=10**(-Ts/20)

plt.loglog(f,module)
plt.axhline(Tc_val_nat, color='g')
plt.axhline(Ts_val_nat, color='g')
plt.axvline(fc, color='g')
plt.ylabel("Module")
plt.xlabel("Fréquence (Hz)");

Comparaison des différents filtres

Les figures suivantes présentent le module et l'argument des 4 différents filtres.

In [574]:
H_list = [H_butter, H_cheby, H_bessel, H_ellip]
nom_list = ["Butter","Chebychev","Bessel", "Elliptique"]


for indice in range(len(H_list)):
    H = H_list[indice]
    nom = nom_list[indice]
    [w,r_w] =H.freqresp(w)
    
    module = np.abs(r_w)
    argument = 180*np.unwrap(np.angle(r_w))/np.pi
    retard_groupe = -np.diff(np.unwrap(np.angle(r_w)))/np.diff(w)
    
    plt.figure(1)
    plt.loglog(f,module,label=nom)
    
    plt.figure(2)
    plt.semilogx(f,argument,label=nom)
    
    plt.figure(3)
    plt.semilogx(f[1:],retard_groupe,label=nom)
    
plt.figure(1)
plt.ylabel("Module")
plt.xlabel("Fréquence (Hz)");
plt.legend()

plt.figure(2)
plt.ylabel("Argument (deg)")
plt.xlabel("Fréquence (Hz)");
plt.legend()

plt.figure(3)
plt.ylabel("Retard de groupe (sec)")
plt.xlabel("Fréquence (Hz)");
plt.ylim([0,0.01])
plt.legend();