Conception et mise en œuvre des filtres numériques

2.a. Relation de récurrence

On considère un signal numérique x_n, obtenu avec une période d'échantillonnage T_e.

Le filtrage numérique linéaire consiste à obtenir un nouveau signal numérique y_n en utilisant une relation de récurrence de la forme suivante :

$$y_n=\sum _{k=0}^{N-1}{b}_k{x}_{n-k}-\sum _{k=1}^{M-1}{a}_k{y}_{n-k}\text{(1)}$$

Les coefficients a_n et b_n sont réels. D'une manière générale, un échantillon de la sortie est obtenu en faisant une combinaison linéaire des N échantillons précédents de l'entrée et des M-1 échantillons précédents de la sortie.

Considérons par exemple le filtre à 5 coefficients suivant :

$$y_n=b_0{x}_n+b_1{x}_{n-1}+b_2{x}_{n-2}-a_1{y}_{n-1}-a_2{y}_{n-2}\text{(2)}$$

Il peut être représenté sous la forme d'un schéma bloc :

Figure pleine page

2.b. Réponse fréquentielle

Lors de l'étude d'un filtre numérique, on est amené à rechercher comment le filtre répond à une entrée sinusoïdale, c'est-à-dire à un signal sinusoïdal échantillonné à la période T_e. Un tel signal a la forme suivante :

$$x_n=exp\left(i2\pi fn{T}_e\right)=z^n\text{(3)}$$

où l'on a introduit la variable complexe :

$$z=exp\left(i2\pi f{T}_e\right)\text{(4)}$$

Le décalage d'une unité dans le temps s'écrit :

$$x_{n-1}=z^{-1}{x}_n\text{(5)}$$

Sur le schéma bloc, ce décalage est donc représenté par l'opérateur z^-1.

Reprenons l'exemple précédent et écrivons la récurrence sous la forme :

$$y_n+a_1{y}_{n-1}+a_2{y}_{n-2}=b_0{x}_n+b_1{x}_{n-1}+b_2{x}_{n-2}\text{(6)}$$

Le terme de droite, combinaison linéaire des entrées, s'écrit pour une entrée sinusoïdale :

$$b_0{x}_n+b_1{x}_{n-1}+b_2{x}_{n-2}=z(b_0+b_1{z}^{-1}+b_2{z}^{-2})\text{(7)}$$

Si l'on admet que la sortie est sinusoïdale (en régime stationnaire), on peut l'écrire sous la forme :

$$y_n=H\left(z\right)x_n=H\left(z\right)z^n\text{(8)}$$

H(z) est par définition la fonction de transfert en z du filtre. En reportant l'expression de y_n dans la relation de récurrence, on obtient l'expression de la fonction de transfert en fonction des coefficients du filtre :

$$H\left(z\right)=\frac{b_0+b_1{z}^{-1}+b_2{z}^{-2}+\cdots }{1+a_1{z}^{-1}+a_2{z}^{-2}+\cdots }\text{(9)}$$

Pour obtenir la réponse fréquentielle du filtre, c'est-à-dire le gain et le déphasage en fonction de la fréquence, il suffit de remplacer z par l'expression $\mathrm{(4)}$.

On voit que cette réponse fréquentielle est une fonction du rapport de la fréquence par la fréquence d'échantillonnage.

2.c. Réponse impulsionnelle et stabilité

Une impulsion est un signal numérique défini par :

$$\begin{array}{ccc}\hfill & x_0=1\hfill & \hfill \text{(10)}\\ \hfill & x_n=0sin\ne 0\hfill & \hfill \text{(11)}\end{array}$$

La réponse impulsionnelle est la sortie y_n obtenue pour cette entrée.

La réponse impulsionnelle est dite finie si le nombre d'échantillons non nuls est fini. Un filtre est stable si la réponse impulsionnelle tend vers zéro lorsque n tend vers l'infini. Un filtre à réponse impulsionnele finie (RIF) est donc toujours stable.

Pour étudier la stabilité d'un filtre, il faut déterminer les pôles de sa fonction de transfert en Z, c'est-à-dire les racines de son dénominateur. Un filtre à réponse impulsionnelle infinie (RII) est stable si et seulement si tous les pôles de sa fonction de transfert ont un module strictement inférieur à 1. Lorsqu'un filtre RII est stable, sa réponse impulsionnelle décroît de manière exponentielle avec n, si bien qu'elle peut être en pratique considérée comme finie.

Un filtre RII instable a une réponse impulsionnelle qui présente soit un comportement oscillatoire, soit une divergence vers l'infini.

3.a. Filtre à réponse impulsionnelle finie

Lorsque les coefficients a_n sont tous nuls, la sortie s'exprime comme une combinaison linéaire des échantillons de l'entrée :

$$y_n=\sum _{k=0}^{N-1}{b}_k{x}_{n-k}\text{(12)}$$

L'opération ainsi réalisée est un produit de convolution, c'est pourquoi on parle dans ce cas de filtrage par convolution. Ce type de filtrage est aussi employé en traitement d'image.

La réponse impulsionnelle est :

$$y_n=b_n\text{(13)}$$

Autrement dit, les coefficients b_n constituent la réponse impulsionnelle du filtre. Cette réponse impulsionnelle est dite finie car le nombre d'échantillons non nuls en sortie est fini, égal au nombre de coefficients b_n non nuls.

Le type le plus répandu de filtre à réponse impulsionnelle finie (filtre RIF) est le filtre RIF à phase linéaire, dont les coefficients sont obtenus par série de Fourier. Voir le document Filtres à réponse impulsionnelle finie pour une explication détaillée de cette méthode. On se contente ici de montrer comment calculer un tel filtre avec la fonction scipy.signal.firwin.

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

Pour définir un filtre passe-bas, on procède de la manière suivante :

P=20
b1 = scipy.signal.firwin(numtaps=2*P+1,cutoff=[0.1],window='hann',nyq=0.5)
                  

numptaps est le nombre de coefficients du filtre. On définit un nombre impair, pour faciliter la suite de l'étude. cutoff est une liste de fréquences de coupure (une seule pour un passe-bas). Ces fréquences sont relatives à la fréquence d'échantillonnage. La fréquence de Nyquist est donc 0.5, ce qui est précisé par l'argument nyq (la valeur par défaut est 1). Enfin l'argument window permet de choisir une fenêtre multiplicative pour les coefficients. Par exemple, une fenêtre de Hanning permet de réduire les ondulations dans la bande passante, qui seraient trop importantes avec une fenêtre rectangulaire.

Il est intéressant de tracer la réponse impulsionnelle :

figure(figsize=(10,5))
stem(b1)
xlabel("n")
ylabel("b")
grid()
                    

Figure pleine page

Il s'agit d'une fonction sinus cardinal dont le maximum est à l'indice P. Lors du calcul de y_n, ce maximum est centré sur x_n-P. On peut donc s'attendre à un décalage PT_e entre l'entrée et la sortie.

La réponse fréquentielle peut être obtenue avec la fonction freqz :

w,h=scipy.signal.freqz(b1)

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.unwrap(numpy.angle(h)))
xlabel("f/fe")
ylabel("phase")
grid()
                   

Figure pleine page

Les filtres RIF peuvent être très sélectifs lorsque le nombre de coefficients est élevé. Par ailleurs la phase est parfaitement linéaire par rapport à la fréquence dans la bande passante. Les signaux dont le spectre est dans la bande passante sont donc transmis sans distorsion, avec un décalage constant, égal ici à PT_e.

Pour définir un filtre passe-bande ou passe-bas, il faut donner la valeur False à l'argument pass_zero. Voici par exemple un filtre RIF passe-bande :

P=20
b2 = scipy.signal.firwin(numtaps=2*P+1,cutoff=[0.1,0.3],pass_zero=False,window='hann',nyq=0.5)
w,h=scipy.signal.freqz(b2)
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.unwrap(numpy.angle(h)))
xlabel("f/fe")
ylabel("phase")
grid()
                   

Figure pleine page

3.b. Filtre à réponse impulsionnelle infinie

Lorsque les coefficients a_n sont non nuls, on parle de filtre récursif :

$$y_n=\sum _{k=0}^{N-1}{b}_k{x}_{n-k}-\sum _{k=1}^{M-1}{a}_k{y}_{n-k}\text{(14)}$$

La réponse impulsionnelle d'un tel filtre est infinie, c'est-à-dire qu'elle comporte un nombre infini d'échantillons non nuls (du moins en théorie). Il peut arriver que le filtre soit instable, c'est-à-dire que sa réponse impulsionnelle ne tende pas vers zéro.

Une méthode de calcul de ce type de filtre consiste à partir de la fonction de transfert d'un filtre analogique et à effectuer une transformation bilinéaire pour obtenir la fonction de transfert en z, d'où l'on déduit les coefficients a_n et b_n. Cette méthode est expliquée en détail dans le document Filtre à réponse impulsionnelle infinie.

La fonction scipy.signal.iitfilter permet d'obtenir un filtre RII à partir de fonctions de transfert analogiques standard (Butterworth, Chebychev, etc). Voici par exemple comment obtenir un filtre passe-bas d'ordre 2 avec une fréquence de coupure f_c/f_e=0.1 :

b3,a3 = scipy.signal.iirfilter(N=2,Wn=[0.1*2],btype="lowpass",ftype="butter")
                  

print(a3)
--> array([ 1.       , -1.1429805,  0.4128016])

print(b3)
--> array([ 0.06745527,  0.13491055,  0.06745527])

La fréquence de coupure relative a été multipliée par deux, car la fréquence de Nyquist est par convention égale à 1.

Ce filtre comporte 3 coefficients b_n et 3 coefficients a_n (en comptant a₀=1). Voici sa réponse fréquentielle :

w,h=scipy.signal.freqz(b3,a3)

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.unwrap(numpy.angle(h)))
xlabel("f/fe")
ylabel("phase")
grid()
                   

Figure pleine page

Les filtres RII ont l'avantage de comporter peu de coefficients, d'introduire un faible déphasage (à peu près linéaire dans la bande passante), et d'être très proches des filtres analogiques. En revanche, ils ne permettent pas d'obtenir la très forte sélectivité des filtres RIF, car les filtres RII d'ordre élevé sont instables. À nombre de coefficients identique, ils sont toutefois plus efficaces que les filtres RIF. Les filtres RII sont très utilisés en traitement du son.

L'inconvénient des filtres RII est la possibilité d'obtenir un filtre instable lorsque l'ordre est élevé (à partir de l'ordre 6 environ). Pour savoir si un filtre est stable, il faut calculer ses pôles :

(zeros,poles,gain) = scipy.signal.tf2zpk(b3,a3) 
                     

et vérifier que leur module est strictement inférieur à 1 :

print(numpy.absolute(poles))
--> array([ 0.64249638,  0.64249638])

4.a. Introduction

Le filtrage consiste à appliquer la relation de récurrence $\mathrm{(1)}$ pour calculer chaque échantillon du signal de sortie.

L'implémentation du filtrage en temps réel, sur un processeur de signal (DSP) ou sur microcontrôleur, pose des problèmes spécifiques liés aux ressources disponibles sur ce type d'unité (quantité de mémoire, vitesse de calcul, calculs en virgule flottante, etc).

On s'intéresse ici au filtrage réalisé en Python, sur un signal entièrement disponible en mémoire, sous forme d'un tableau d'échantillons. On montrera l'effet de différents filtres sur un signal obtenu par numérisation à la fréquence d'échantillonnage f_e=10 kHz d'un signal périodique comportant 3 harmoniques. 5000 échantillons ont été acquis. La fréquence du fondamental est f=400,3000 Hz. Le signal a été stocké dans un fichier texte. On commence donc par récupérer les échantillons, avec les instants correspondant :

[t,x] = numpy.loadtxt("signal-400.3.txt")

figure(figsize=(8,5))
plot(t,x)
axis([0,0.02,-2,2])
xlabel("t (s)")
ylabel("x (V)")
grid()

                 

figF.pdf

4.b. Réalisation d'un filtrage RIF par convolution

L'opération à réaliser sur le signal est un produit de convolution :

$$y_n=\sum _{k=0}^{N-1}{b}_k{x}_{n-k}\text{(15)}$$

Le nombre de coefficients N=2P+1 a été défini par numtaps dans la fonction firwin. On remarque que les 2P premières valeurs de y_n ne peuvent pas être calculées avec cette somme. Il y deux manières de procéder pour ces 2P premiers échantillons : soit on les calcule en utilisant une partie des coefficients, soit on ne les filtre pas et on annule la sortie. Nous allons adopter la seconde solution, plus facile à implémenter. Voici une fonction réalisant le filtrage :

def filtrage_convolution(x,b):
    N = len(b)
    ne = len(x)
    y = numpy.zeros(ne)
    for n in range(N-1,ne):
        accum = 0.0
        for k in range(N):
            accum += b[k]*x[n-k]
        y[n] = accum
    return y
                   

On utilise le filtre RIF passe-bas défini plus haut, puis on trace le signal de départ avec le signal filtré :

y = filtrage_convolution(x,b1)
                
figure(figsize=(8,6))
subplot(211)
plot(t,x)
axis([0,0.02,-2,2])
xlabel("t (s)")
ylabel("x(V)")
grid()
subplot(212)
plot(t,y)
axis([0,0.02,-2,2])
xlabel("t (s)")
ylabel("y(V)")
grid()
                

figG1.pdf

Le filtre passe-bas a une fréquence de coupure de 1 kHz. L'harmonique de rang 3 (fréquence 1200 Hz) est donc légèrement réduite, mais l'essentiel du signal se trouve dans la bande passante. On constate qu'il y a un décalage entre la sortie et l'entrée. C'est une conséquence de la variation linéaire du déphasage avec la fréquence dans la bande passante. Pour ce type de filtre RIF, le décalage est τ=PT_e, ici égal à 20 ms.

La fonction utilisée ci-dessus est lente car la boucle est écrite en python. Le temps de calcul est proportionnel au nombre d'échantillons multiplié par le nombre de coefficients du filtre. Pour traiter rapidement des grandes listes d'échantillons, il est préférable d'utiliser la fonction scipy.signal.convolve, de la manière suivante :

y = scipy.signal.convolve(x,b1,mode='same')
                

Le mode same permet d'obtenir une sortie de même taille que l'entrée.

figure(figsize=(8,6))
subplot(211)
plot(t,x)
axis([0,0.02,-2,2])
xlabel("t (s)")
ylabel("x(V)")
grid()
subplot(212)
plot(t,y)
axis([0,0.02,-2,2])
xlabel("t (s)")
ylabel("y(V)")
grid()
                

figG2.pdf

Le signal de sortie n'est pas décalé par rapport au signal d'entrée. Cela est du au fait que la convolution appliquée est en fait centrée. Si N=2P+1 est le nombre de coefficients, la relation utilisée est :

$$y_n=\sum _{k=-P}^P{b}_{P-k}{x}_{n+k}\text{(16)}$$

Pour les P premiers et les P derniers points, une partie de la somme est calculée.

Ce mode de calcul de la convolution, appelé convolution centrée, convient bien au filtrage des images, où l'on veut que chaque pixel filtré coïncide avec le pixel d'origine, mais il ne correspond pas à la réalité du filtrage numérique en temps-réel, où les échantillons du futur ne sont évidemment pas disponibles lorsqu'on calcule y_n.

Si l'on souhaite utiliser la fonction convolve et reproduire le décalage, on peut utiliser le mode valid, qui permet de limiter le calcul de la convolution aux points valides, c'est-à-dire à tous les points sauf les P premiers et les P derniers. Il faut aussi construire une nouvelle échelle de temps pour tracer la sortie, en remarquant que le premier point calculé est en fait celui d'indice 2P.

y = scipy.signal.convolve(x,b1,mode='valid')
Te = t[1]-t[0]
ty = (numpy.arange(y.size)+2*P)*Te
figure(figsize=(8,4))
plot(t,x,'b')
plot(ty,y,'r')
axis([0,0.02,-2,2])
xlabel("t (s)")
ylabel("y(V)")
grid()

                 

figG3.pdf

On obtient ainsi la sortie avec le décalage qu'elle aurait dans un filtre temps-réel. Le retard de la sortie est ici de 2 ms.

4.c. Réalisations d'un filtrage récursif