Table des matières Python

Conception d'un filtre par placement des zéros et des pôles

1. Introduction

Ce document montre comment un filtre numérique récursif peut être conçu en plaçant ses zéros et ses pôles. Cette méthode est très efficace pour concevoir des filtres dont la fréquence de coupure est très basse (une fraction infime de la fréquence d'échantillonnage). On verra en particulier comment obtenir un filtre passe-bas ou passe-bande avec un domaine intégrateur très large, s'étendant d'une fréquence très faible jusqu'à la fréquence de Nyquist.

2. Principe

On considère un filtre numérique linéaire défini par la relation de récurrence suivante :

yn=k=0Nbkxn-k-k=1Makyn-k(1)

Les M+1 coefficients an et les N+1 bn sont réels.

Pour étudier la réponse fréquentielle du filtre, on introduit la variable complexe Z définie par :

z=exp(i2πfTe)=exp(iΩ)(2)

f est la fréquence, Te la période d'échantillonnage. Pour simplifier les notations, on utilise aussi la pulsation réduite Ω. Lorsque la fréquence varie entre 0 et la fréquence de Nyquist (moitié de la fréquence d'échantillonnage), la pulsation Ω varie entre 0 et π. Le nombre Z se déplace donc sur le demi-cercle unité.

La réponse fréquentielle du filtre est obtenue avec la fonction de transfert en Z :

H(z)=b0+b1z-1+b2z-2++bNz-N1+a1z-1+a2z-2++aMz-M(3)

Pour une démonstration de ce résultat, voir Conception et mise en œuvre des filtres numériques.

On peut écrire la fonction de transfert comme un rapport de deux polynômes :

H(z)=b0zN+b1zN-1+b2zN-1+bNzM+a1zM-1+a2zM-2++aM(4)

Les zéros sont les N racines du numérateur; on les notera qi. Les pôles sont les racines du dénominateur; on les notera pi. Pour que le filtre soit stable, il faut et il suffit que tous les pôles soient à l'intérieur (strictement) du cercle unité. Les coefficients du filtre étant réels, chaque pôle (ou zéro) non réel est associé à son complexe conjugué.

Voici par exemple une fonction de transfert bi-quadratique (M=N=2) écrite avec ses pôles et zéros :

H(z)=b0(z-q1)(z-q2)(z-p1)(z-p2)(5)

La méthode de placement des pôles ([1],[2]) consiste à définir le filtre en plaçant directement les pôles et les zéros.

Supposons que l'on souhaite annuler H pour une certaine pulsation Ωa non nulle. Il faut définir un zéro de module 1 et d'argument Ωa, et son conjugué :

q1=exp(iΩa)(6)q2=exp(-iΩa)(7)

Si la pulsation est nulle, un seul zéro réel suffit. Si l'on veut obtenir pour cette pulsation un gain faible mais non nul, il suffit de donner à ce zéro un module différent de 1 (inférieur ou supérieur à 1).

Pour définir une résonance, c'est-à-dire un maximum du gain, il faut agir sur les pôles. Supposons que l'on souhaite obtenir une résonance pour une pulsation Ωb. Si cette pulsation est non nulle, il faut définir deux pôles conjugués :

p1=r1exp(iΩb)(8)p2=r1exp(-iΩb)(9)

Si Ωb=0, un seul pôle réel suffit. Le module r1 doit être strictement inférieur à 1 pour que le filtre soit stable. Plus ce module est proche de 1, plus la résonance est forte (plus le maximum est haut). On peut dans certains cas choisir un module égal à 1, ce qui donne une instabilité pour la pulsation Ωb (le gain tend vers l'infini pour cette pulsation).

Les pôles et les zéros sont représentés graphiquement dans le plan complexe. Les pôles sont réprésentés par des croix, les zéros par des cercles.

planZ-1.svgFigure pleine page

Lorsque les pôles et les zéros sont définis, il faut déterminer la constante b0 en fonction du gain souhaité, par exemple le gain maximal pour un filtre passe-bas.

La méthode de placement des pôles et zéros est une méthode qualitative. Il faut tracer la réponse fréquentielle pour obtenir le comportement quantitatif du filtre.

3. Filtres du premier ordre

3.a. Intégrateur

Un intégrateur parfait est défini avec un pôle à fréquence nulle p1=1 et un zéro à la fréquence de Nyquist q1=-1 :

planZ-2.svgFigure pleine pageH(z)=z+1z-1=1+z-11-z-1(10)

La seconde écriture permet d'obtenir la relation de récurrence, sachant qu'au numérateur z-1 correspond à un retard d'une unité pour l'entrée, au dénominateur un retard d'une unité pour la sortie :

yn=yn-1+xn+xn-1(11)

La réponse fréquentielle est :

H(Ω)=eiΩ+1eiΩ-1(12)

Le gain est infini à pulsation nulle, ce qui signifie que le filtre est instable si le signal d'entrée possède une fréquence nulle dans son spectre. L'instabilité vient de la présence d'un pôle sur le cercle unité. Voici le tracé du gain et du déphasage :

import numpy
import scipy.signal
from matplotlib.pyplot import *
b=[1,1]
a=[1,-1]
w,h=scipy.signal.freqz(b,a)

figure()
subplot(211)
plot(w/(2*numpy.pi),20*numpy.log10(numpy.absolute(h)))
xlabel("f/fe")
ylabel("GdB")
grid()
subplot(212)      
plot(w/(2*numpy.pi),numpy.angle(h)/numpy.pi)
xlabel("f/fe")
ylabel("phase/pi")
axis([0,0.5,-1,1])
grid()
                   
figA.svgFigure pleine page

Pour ce filtre, il n'y a pas de gain maximal pour fixer la constante b0. On fixera cette constante en fonction de l'objectif.

Supposons que l'on cherche à réaliser une intégration vraie, qui correspond à la relation en temps continu suivante :

y(t)=0tx(t)dt(13)

Dans ce cas, la relation de récurrence est

yn=yn-1+Te2(xn+xn-1)(14)

qui correspond à la méthode des trapèzes pour le calcul numérique d'une intégrale. Sa fonction de transfert est :

H(z)=Te2z+1z-1=Te21+z-11-z-1(15)

Le filtre intégrateur joue un rôle important dans la conception des filtres, car il sert de référence dans la méthode de transformation bilinéaire, expliquée dans Filtres à réponse impulsionnelle infinie.

3.b. Filtre passe-bas

En partant du filtre intégrateur précédent, on peut déplacer le pôle sur l'axe des réels pour lui donner un module r inférieur à 1, ce qui a pour effet de stabiliser le filtre.

planZ-3.svgFigure pleine pageH(z)=b0z+1z-r(16)

Il s'agit d'un filtre passe-bas dont le gain dans la bande passante (z=1) est :

G=b021-r(17)

On peut par exemple choisir b0 pour que G=1.

La relation de récurrence est :

yn=ryn-1+b0(xn+xn-1)(18)

Voici un exemple :

r=0.8
G=1
b0 = (1-r)*G/2
b=[b0,b0]
a=[1,-r]
w,h=scipy.signal.freqz(b,a)

figure()
subplot(211)
plot(w/(2*numpy.pi),20*numpy.log10(numpy.absolute(h)))
xlabel("f/fe")
ylabel("GdB")
grid()
subplot(212)      
plot(w/(2*numpy.pi),numpy.angle(h)/numpy.pi)
xlabel("f/fe")
ylabel("phase/pi")
axis([0,0.5,-1,1])
grid()
                 
figB.svgFigure pleine page

Le filtre obtenu est similaire à une filtre du premier ordre calculé par transformation bilinéaire d'une fonction de transfert analogique. La méthode de placement des pôles et zéros permet d'obtenir facilement un filtre passe-bas ayant une très basse fréquence de coupure, par exemple :

r=0.99
G=1
b0 = (1-r)*G/2
b=[b0,b0]
a=[1,-r]
w,h=scipy.signal.freqz(b,a)

figure()
subplot(211)
plot(w/(2*numpy.pi),20*numpy.log10(numpy.absolute(h)))
xlabel("f/fe")
ylabel("GdB")
grid()
subplot(212)      
plot(w/(2*numpy.pi),numpy.angle(h)/numpy.pi)
xlabel("f/fe")
ylabel("phase/pi")
axis([0,0.5,-1,1])
grid()
                 
figC.svgFigure pleine page

On obtient ainsi un intégrateur stable, dont le gain à fréquence nulle est fini. Pour cette utilisation, il peut être utile d'ajuster la constante b0 pour obtenir une intégration vraie du signal d'entrée. On doit alors avoir dans la plage de fréquence d'intégration :

b0z+1z-r=Te2z+1z-1(19)

ce qui donne :

b0=Te2z-rz-1(20)

Si r est très proche de 1, et si z est différent de 1 (fréquence non nulle), la fraction est pratiquement égale à 1. Cela démontre le comportement intégrateur et indique le choix de la constante :

b0=Te2(21)

On obtient finalement la relation de récurrence suivante pour une utilisation de ce filtre en tant qu'intégrateur :

yn=ryn-1+Te2(xn+xn-1)(22)

La constante r est inférieure à 1 mais proche de 1, d'autant plus que l'on veut une fréquence de coupure basse. Pour r=1, on retrouve l'intégrateur parfait, instable à fréquence nulle.

Pour une utilisation en tant que filtre passe-bas avec un gain unité dans la bande passante, la relation de récurrence est :

yn=ryn-1+1-r2(xn+xn-1)(23)

On peut chercher à obtenir ce filtre passe-bas directement à partir de l'équation différentielle du système à temps continu. Pour une pulsation analogique de coupure ωc et un gain dans la bande passante g, cette équation est :

dy(t)dt+ωcy(t)=ωcgx(t)(24)

La méthode d'intégration d'Euler s'écrit :

yn=yn-1+ωcTe2(gxn-1-yn-1)(25)=(1-ωcTe)yn-1+ωcTegxn-1(26)=ryn-1+(1-r)gxn-1(27)

Cette relation est différente de la relation (23) puisque la méthode d'Euler est une méthode à un pas. La fonction de transfert correspondante a l'inconvénient d'être démunie de zéro. Voyons sa réponse fréquentielle :

b=[0,1-r]
a=[1,-r]
w,h=scipy.signal.freqz(b,a)

figure()
subplot(211)
plot(w/(2*numpy.pi),20*numpy.log10(numpy.absolute(h)))
xlabel("f/fe")
ylabel("GdB")
grid()
subplot(212)      
plot(w/(2*numpy.pi),numpy.angle(h)/numpy.pi)
xlabel("f/fe")
ylabel("phase/pi")
axis([0,0.5,-1,1])
grid()
                 
figC1.svgFigure pleine page

Cette réponse est très différente de la réponse fréquentielle du filtre analogique. En particulier, on ne retrouve pas du tout le comportement intégrateur dans la bande atténuante. La solution consiste à appliquer une méthode numérique à deux pas, du type Adams-Moulton ([3]), comme on le fait pour l'intégrateur (méthode des trapèzes) :

yn=(1-ωcTe)yn-1+ωcTe12g(xn+xn-1)(28)

qui est bien de la forme (23) si l'on pose r=1-ωcTe et g=1. La relation réalisant une intégration vraie dans la bande atténuante s'obtient en posant :

ωcg=1(29)

3.c. Filtre passe-haut

Un filtre passe-haut est obtenu avec un pôle p1=-r pour la fréquence de Nyquist et un zéro q1=1 pour la fréquence nulle.

planZ-4.svgFigure pleine pageH(z)=b0z-1z+r(30)

La relation de récurrence est :

yn=-ryn-1+b0(xn-xn-1)(31)

Le gain dans la bande passante (z=-1) est :

G=b021-r(32)

qui est la même relation que pour le filtre passe-bas. Voici un exemple :

r=0.2
G=1
b0 = (1-r)*G/2
b=[b0,-b0]
a=[1,r]
w,h=scipy.signal.freqz(b,a)

figure()
subplot(211)
plot(w/(2*numpy.pi),20*numpy.log10(numpy.absolute(h)))
xlabel("f/fe")
ylabel("GdB")
grid()
subplot(212)      
plot(w/(2*numpy.pi),numpy.angle(h)/numpy.pi)
xlabel("f/fe")
ylabel("phase/pi")
axis([0,0.5,-1,1])
grid()
                 
figD.svgFigure pleine page

Pour obtenir un filtre bloquant la composante continue, il faut abaisser la fréquence de coupure. Cela se fait en plaçant le pôle en p1=r, avec r très proche de 1 :

planZ-5.svgFigure pleine pageH(z)=b0z-1z-r(33)

Le gain dans la bande passante est :

G=b021+r(34)

La relation de récurrence est :

yn=ryn-1+b0(xn-xn-1)(35)

Voici un exemple. Si est r est très proche de 1, on a pratiquement b0=1.

r=0.999
G=1
b0 = (1+r)*G/2
b=[b0,-b0]
a=[1,-r]
w,h=scipy.signal.freqz(b,a)

figure()
subplot(211)
plot(w/(2*numpy.pi),20*numpy.log10(numpy.absolute(h)))
xlabel("f/fe")
ylabel("GdB")
grid()
subplot(212)      
plot(w/(2*numpy.pi),numpy.angle(h)/numpy.pi)
xlabel("f/fe")
ylabel("phase/pi")
axis([0,0.5,-1,1])
grid()
                 
figE.svgFigure pleine page

3.d. Filtre passe-bande

Un filtre passe-bande peut être obtenu en associant en série un filtre passe-bas et un filtre passe-haut :

planZ-6.svgFigure pleine pageH(z)=b0z-1z-r1z+1z-r2=b01-z-21-(r1+r2)z-1+r1r2z-2(36)
r1=0.9
r2=0.95

b=[1,0,-1]
a=[1,-(r1+r2),r1*r2]
w,h=scipy.signal.freqz(b,a)

figure()
subplot(211)
plot(w/(2*numpy.pi),20*numpy.log10(numpy.absolute(h)))
xlabel("f/fe")
ylabel("GdB")
grid()
subplot(212)      
plot(w/(2*numpy.pi),numpy.angle(h)/numpy.pi)
xlabel("f/fe")
ylabel("phase/pi")
axis([0,0.5,-1,1])
grid()
                 
figF.svgFigure pleine page

Une application de ce type de filtre est un intégrateur avec un gain nul à fréquence nulle, obtenu avec r1 et r2 très proches de 1.

r1=0.98
r2=0.99

b=[1,0,-1]
a=[1,-(r1+r2),r1*r2]
w,h=scipy.signal.freqz(b,a)

figure()
subplot(211)
plot(w/(2*numpy.pi),20*numpy.log10(numpy.absolute(h)))
xlabel("f/fe")
ylabel("GdB")
grid()
subplot(212)      
plot(w/(2*numpy.pi),numpy.angle(h)/numpy.pi)
xlabel("f/fe")
ylabel("phase/pi")
axis([0,0.5,-1,1])
grid()
                 
figG.svgFigure pleine page

Voyons la réponse impulsionnelle de ce filtre :

x = numpy.zeros(1000)
x[0] = 1.0
zi = scipy.signal.lfiltic(b,a,y=[0,0],x=[0,0])
[y,zf] = scipy.signal.lfilter(b,a,x,zi=zi)

figure()
plot(y)
xlabel("n")
ylabel("y")
grid()
                   
figG1.svgFigure pleine page

Il faut environ 600 échantillons pour atteindre le régime stationnaire. Pour obtenir une réponse plus rapide, il faut augmenter la fréquence de la bande passante, en réduisant r1 et r2.

Pour obtenir une intégration vraie dans la bande atténuante, il faut choisir la constante b0 pour que l'égalité suivante soit à peu près vérifiée dans la bande atténuant :

b0z-1z-r1z+1z-r2=Te2z+1z-1(37)

Lorsque r1 et r2 sont très proches de 1 et lorsque z est différent de 1 (fréquence non nulle), cette égalité est approximativement vérifiée si :

b0=Te2(38)

Voici la relation de récurrence de cet intégrateur :

yn=(r1+r2)yn-1-r1r2yn-2+Te2(xn-xn-2)(39)

Si le signal d'entrée comporte une composante de fréquence nulle (composante continue), l'intégration se fait sans décalage en sortie. L'inconvénient est la lenteur de la réponse transitoire. En pratique, il faut ajuster la fréquence de la bande passante juste en dessous de la plus basse fréquence du signal.

4. Filtres du second ordre

4.a. Filtre passe-bande

Le filtre passe-bande précédent était défini avec deux pôles réels. Pour obtenir un résonateur, on définit deux pôles complexes conjugués à la pulsation de résonance, avec deux zéros q1=1 et q2=-1 qui permettent d'annuler le gain à fréquence nulle et à la fréquence de Nyquist :

planZ-7.svgFigure pleine pageH(z)=b0(z-1)(z+1)(z-reiΩa)(z-re-iΩa)=b01-z-21-2rcos(Ωa)z-1+r2z-2(40)

Voici un exemple, avec une fréquence de résonance égale à 1/5 ième de la fréquence de Nyquist :

omega_a=numpy.pi/5
r=0.95
b=[1,0,-1]
a=[1,-2*r*numpy.cos(omega_a),r*r]
w,h=scipy.signal.freqz(b,a)

figure()
subplot(211)
plot(w/(2*numpy.pi),20*numpy.log10(numpy.absolute(h)))
xlabel("f/fe")
ylabel("GdB")
grid()
subplot(212)      
plot(w/(2*numpy.pi),numpy.angle(h)/numpy.pi)
xlabel("f/fe")
ylabel("phase/pi")
axis([0,0.5,-1,1])
grid()

                  
figH.svgFigure pleine page

Il faudra bien sûr choisir une constante b0 adéquate pour obtenir un gain unité dans la bande passante.

Voyons la réponse impulsionnelle de ce filtre :

x = numpy.zeros(100)
x[0] = 1.0
zi = scipy.signal.lfiltic(b,a,y=[0,0],x=[0,0])
[y,zf] = scipy.signal.lfilter(b,a,x,zi=zi)

figure()
stem(y)
xlabel("n")
ylabel("y")
grid()
                   
figI.svgFigure pleine page

Elle présente des oscillations à la fréquence de résonance, égale à un dixième de la fréquence d'échantillonnage. Elle tend bien vers zéro puisque ce filtre est stable (ses pôles ont un module strictement inférieur à 1).

Références
[1]  J.O. Smith,  Introduction to digital filters, with audio applications,  (BookSurge, 2007)
[2]  Tan Li, Jiang Jean,  Digital signal processing : fundamentals and applications,  (Elsevier, 2013)
[3]  J.P. Demailly,  Analyse numérique et équations différentielles,  (P.U.G., 1991)
Creative Commons LicenseTextes et figures sont mis à disposition sous contrat Creative Commons.