Filtre IIR ( Infinite Impulse Response )

Un filtre numérique IIR ( Infinite Impulse Response ) résulte de la conversion d’un filtre analogique en filtre échantillonné.

A titre d’exemple, nous considérons le filtre analogique Passe Bas suivant :

Pente : -80 db/dec ( Filtre d’ordre 4 )

Nous allons dans un premier temps concevoir un filtre analogique de Butterworth pour répondre à ce cahier des charges.

1 - Filtres de Butterworth

L’objectif est d’obtenir un filtre de pente (-20*N dB/dec, N ordre du filtre), en ayant 3dB d’atténuation à la fréquence de coupure.

1.1 - Définition

$\boxed{ \left|H(\omega)\right|^2=\frac{1}{1+\left(\frac{\omega}{\omega_c} \right )^{2N}} }$

butterworth.py

from scipy import signal
from numpy import *
import matplotlib.pyplot as plt

N=1
cutoff_freq=5000

for N in range(1,6):
  b, a = signal.butter(N, cutoff_freq, 'low', analog=True)
  w, h = signal.freqs(b, a)
  plt.semilogx(w, 20 * np.log10(abs(h)),label=N)  

plt.title('Butterworth filter frequency response')
plt.xlabel('Frequency [radians / second]')
plt.ylabel('Amplitude [dB]')
plt.margins(0, 0.1)
plt.grid(which='both', axis='both')
plt.axvline(cutoff_freq, color='green') # cutoff frequency
plt.legend()
plt.show()

Cas N=1 :

Soit une cellule passe-bas du 1er ordre :

H(j\omega)=\frac{1}{1+j\frac{\omega}{\omega_c}}

\left|H(\omega)\right|=\frac{1}{\left|1+\frac{\omega}{\omega_c} \right|}

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

\left|H(\omega)\right|^2=\frac{1}{1+\left(\frac{\omega}{\omega_c}\right)^2}

On retrouve le cas particulier du filtre de Butterworth avec N=1.

Cas N=2 :

Soit une cellule passe-bas du 2nd ordre , avec $\xi =\frac{\sqrt{2}}{2}$ .
Ce cas particuliers permet de respecter l’atténuation de -3dB à $f_c$

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

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

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

…

\left|H(\omega)\right|^2=\frac{1}{1+\left(\frac{\omega}{\omega_c}\right)^4}

On retrouve le cas particulier du filtre de Butterworth avec N=2.

1.2 - Polynomes de Butterworth Normalisés

$s=j\frac{\omega}{\omega_c}$

$N \quad \text{Polynome}$
$1 \quad s+1$ $2 \quad s^2+1.4142s+1$
$3 \quad (s+1)(s^2+s+1)$
$4 \quad (s^2 + 0.765367s + 1)(s^2 + 1.847759s + 1)$
$5 \quad (s+1)(s^2+0.618s+1)(s^2+1.618s+1)$
$6 \quad (s^2 + 0.5176387s + 1)(s^2 + 1.414214 + 1)(s^2 +1.931852s + 1)$

etc…

2 - Conception d’un Filtre IIR

Nous allons un filtre numérique IIR, basé sur un filtre analogique de Butterworth, avec les caractéristiques suivantes :

Fréquence d’échantillonnage : $f_e=44100Hz$
Ordre N=4 ( atténuation -80 dB/décade )
Fréquence de coupure : $f_c=5000Hz$
Résonnance : Q=1

2.1 - Polynôme de Butterworth

$H(s)=\frac{Y(s)}{X(s)}=\frac{1}{s^2+\frac{1}{Q}.0.765367.s+1}.\frac{1}{s^2+\frac{1}{Q}.1.847759.s+1}$

$H(f)=\frac{Y(s)}{X(s)}=\frac{1}{\left (j\frac{f}{f_c}\right )^2+j\frac{1}{Q}.0.765367.\frac{f}{f_c}+1}. \frac{1}{\left (j\frac{f}{f_c}\right )^2+j\frac{1}{Q}.1.847759.\frac{f}{f_c}+1}$

2.2 - Transformation bilinéaire

Nous avons évoqué dans la partie echantillonnage une méthode d’intégration basée sur la méthode des trapèzes : la transformation bilinéaire.
En pratique il suffit de réaliser la substitution suivante afin de convertir un filtre analogique en un filtre numérique :

$s \rightarrow \frac{2}{T_e}\frac{1-z^{-1}}{1+z^{-1}}$

Conséquence :

un point appartenant à l’axe $s=j\omega$ dans le plan de Laplace est déplacé en un point appartenant au cercle de rayon unité $\exp{j\omega}$ dans le domaine Z.
Après transformation, la pulsation $\omega_a$ dans le domaine continu est déplacée à la pulsation $\omega_c$ dans le domaine numérique. Le lien entre ces deux pulsations est donné par:

$\omega_a = \frac{2}{T_e}\tan\left( \omega_c\frac{T_e}{2}\right)$

2.3 - Méthode de conception du filtre

1- A partir de l’ordre du filtre indiqué par le cahier des charges, on en déduit les polynômes de Butterworth correspondants.
2- Comme la transformée bilinéaire induit une distorsion non linéaire de l’axe fréquentiel, on corrige la fréquence de coupure souhaitée ( prewarping)

$\omega_a = \frac{2}{T_e}\tan\left( \omega_c\frac{T_e}{2}\right)$

3- On applique la transformée bilinéaire

$s \rightarrow \frac{2}{T_e}\frac{1-z^{-1}}{1+z^{-1}}$

Filtre IIR d’Ordre 2

synthese_iir_ordre2.py

synthese_iir_ordre2.py

import numpy as np
from numpy import pi
from scipy.signal import butter, lti, dlti
from scipy import signal
import matplotlib.pyplot as plt

fc = 440.0
fs = 44100.0
Ts = 1/fs
Q = 1.0
#==============================================
# Forme standard : 
#        a2 * s^2 + a1 * s + a0
# H(s) = ----------------------
#        b2 * s^2 + b1 * s + b0  
#==============================================
# Polynome butterworth normalise
# (1 / (s^2 + (1/Q) * 1.4142s + 1)) 
a0 = 1.0
a1 = 0
a2 = 0
b0 = 1.0
b1 = 1.4142/Q
b2 = 1.0
#==============================================
# prewarp(a0, a1, a2)
wa = 2.0*fs*np.tan(pi*fc/fs);
a2 = a2/(wa**2)
a1 = a1/wa 
# prewarp(b0, b1, b2)
b2 = b2/(wa**2)
b1 = b1/wa
#==============================================
# Transformée Bilinéaire
#        1 + alpha1 * z^(-1) + alpha2 * z^(-2)
# H(z) = -------------------------------------
#        1 + beta1 * z^(-1) + beta2 * z^(-2)

#        z^2 + alpha1 * z + alpha2 
# H(z) = --------------------------
#        z^2 + beta1 * z + beta2 

ad = 4.0*a2*fs**2+ 2.0*a1*fs + a0;
bd = 4.0*b2*fs**2+ 2.0*b1*fs + b0;
k  = ad/bd
S0_beta1  = (2.0*b0 - 8.0*b2*fs**2)/bd
S0_beta2  = (4.0*b2 * fs**2 - 2.0*b1*fs + b0)/bd
S0_alpha1 = (2.0*a0 - 8.0*a2*fs**2)/ad
S0_alpha2 = (4.0*a2*fs**2 - 2.0*a1*fs + a0)/ad
S0_k = k
#========================================
num0 = [1,S0_alpha1,S0_alpha2]  
den0 = [1,S0_beta1,S0_beta2]

Hd0 = dlti(S0_k*np.array(num0),np.array(den0))

print("Hd=",Hd)
print("k :",S0_k)
print("---------STAGE 0---------")
print("S0_beta1 :",S0_beta1)
print("S0_beta2 :",S0_beta2)
print("S0_alpha1 :",S0_alpha1)
print("S0_alpha2 :",S0_alpha2)

w, H_1jw = Hd0.freqresp()
plt.semilogx(fs*w/(2*np.pi),20*np.log10(abs(H_1jw)), label="discrete bilinear manual")

b, a = signal.butter(2, fc, 'low', analog=True)
w, h = signal.freqs(b, a)
plt.semilogx(w, 20*np.log10(abs(h)),label="analog") 

plt.ylabel("Module $H(e^{j\omega})\ en\ dB$")
plt.xlabel("Frequence [Hz]")
plt.ylim([-60, 4])
plt.legend()
plt.grid()
plt.show()

Filtre IIR d’Ordre 4

synthese_iir_ordre4.py

synthese_iir_ordre4.py

import numpy as np
from numpy import pi
from scipy.signal import butter, lti, dlti
from scipy import signal
import matplotlib.pyplot as plt

fc = 440.0
fs = 44100.0
Ts = 1/fs
Q = 1
#==============================================
# Forme standard : 
#        a2 * s^2 + a1 * s + a0
# H(s) = ----------------------
#        b2 * s^2 + b1 * s + b0  
#==============================================
# Polynome butterworth normalise
# (1 / (s^2 + (1/Q) * 0.765367s + 1)) * (1 / (s^2 + (1/Q) * 1.847759s +1))
# Nous avons 2 polynomes de degré 2, S0 = Stage 0, S1 = Stage 1
S0_a0 = 1.0
S0_a1 = 0
S0_a2 = 0
S0_b0 = 1.0
S0_b1 = 0.765367/Q
S0_b2 = 1.0

S1_a0 = 1.0
S1_a1 = 0
S1_a2 = 0
S1_b0 = 1.0
S1_b1 = 1.847759/Q
S1_b2 = 1.0
#==============================================
# TRAITEMENT STAGE 0
#==============================================
a0 = S0_a0
a1 = S0_a1
a2 = S0_a2
b0 = S0_b0
b1 = S0_b1
b2 = S0_b2
#==============================================
# prewarp(a0, a1, a2)
wa = 2.0*fs*np.tan(pi*fc/fs);
a2 = a2/(wa**2)
a1 = a1/wa 
# prewarp(b0, b1, b2)
b2 = b2/(wa**2)
b1 = b1/wa
#==============================================
# Transformée Bilinéaire
#        1 + alpha1 * z^(-1) + alpha2 * z^(-2)
# H(z) = -------------------------------------
#        1 + beta1 * z^(-1) + beta2 * z^(-2)

#        z^2 + alpha1 * z + alpha2 
# H(z) = --------------------------
#        z^2 + beta1 * z + beta2 

ad = 4.0*a2*fs**2+ 2.0*a1*fs + a0;
bd = 4.0*b2*fs**2+ 2.0*b1*fs + b0;
k  = ad/bd
S0_beta1  = (2.0*b0 - 8.0*b2*fs**2)/bd
S0_beta2  = (4.0*b2 * fs**2 - 2.0*b1*fs + b0)/bd
S0_alpha1 = (2.0*a0 - 8.0*a2*fs**2)/ad
S0_alpha2 = (4.0*a2*fs**2 - 2.0*a1*fs + a0)/ad
S0_k = k
#==============================================
# TRAITEMENT STAGE 1
#========================================
a0 = S1_a0
a1 = S1_a1
a2 = S1_a2
b0 = S1_b0
b1 = S1_b1
b2 = S1_b2

# prewarp(a0, a1, a2, fc, fs)
wa = 2.0*fs*np.tan(pi*fc/fs);
a2 = a2/(wa**2)
a1 = a1/wa 
# prewarp(b0, b1, b2, fc, fs)
b2 = b2/wa**2
b1 = b1/wa
#==============================================
# Transformée Bilinéaire
#        1 + alpha1 * z^(-1) + alpha2 * z^(-2)
# H(z) = -------------------------------------
#        1 + beta1 * z^(-1) + beta2 * z^(-2)

#        z^2 + alpha1 * z + alpha2 
# H(z) = --------------------------
#        z^2 + beta1 * z + beta2 

ad = 4.0*a2*fs**2+ 2.0*a1*fs + a0;
bd = 4.0*b2*fs**2+ 2.0*b1*fs + b0;
k  = ad/bd
S1_beta1  = (2.0*b0 - 8.0*b2*fs**2)/bd
S1_beta2  = (4.0*b2 * fs**2 - 2.0*b1*fs + b0)/bd
S1_alpha1 = (2.0*a0 - 8.0*a2*fs**2)/ad
S1_alpha2 = (4.0*a2*fs**2 - 2.0*a1*fs + a0)/ad
S1_k = k
#==============================================
num0 = [1,S0_alpha1,S0_alpha2]  
den0 = [1,S0_beta1,S0_beta2]

num1 = [1,S1_alpha1,S1_alpha2]  
den1 = [1,S1_beta1,S1_beta2]

num=np.convolve(np.array(num0),np.array(num1))
den=np.convolve(np.array(den0),np.array(den1))
K=S0_k*S1_k

Hd=dlti(K*num,den)

print("Hd=",Hd)
print("k :",K)
print("---------STAGE 0---------")
print("S0_beta1 :",S0_beta1)
print("S0_beta2 :",S0_beta2)
print("S0_alpha1 :",S0_alpha1)
print("S0_alpha2 :",S0_alpha2)
print("---------STAGE 1---------")
print("S1_beta1 :",S1_beta1)
print("S1_beta2 :",S1_beta2)
print("S1_alpha1 :",S1_alpha1)
print("S1_alpha2 :",S1_alpha2)

w, H_1jw = Hd.freqresp()
plt.semilogx(fs*w/(2*np.pi), 20*np.log10(abs(H_1jw)), label="discrete bilinear manual")


b, a = signal.butter(4, fc, 'low', analog=True)
w, h = signal.freqs(b, a)
plt.semilogx(w, 20 * np.log10(abs(h)),label="analog") 

plt.ylabel("Module $H(e^{j\omega})\ en\ dB$")
plt.xlabel("Frequence [Hz]")
plt.ylim([-120, 4])
plt.legend()
plt.grid()
plt.show()
#==============================================

Filtre FIR ( Finite Impulse Response ) Comparaison Filtres FIR / Filtres IIR