Nothing Special   »   [go: up one dir, main page]

Introfiltres PDF

Télécharger au format pdf ou txt
Télécharger au format pdf ou txt
Vous êtes sur la page 1sur 18

Frédéric Legrand Licence Creative Commons 1

Introduction aux ltres numériques

1. Introduction
Le ltrage du signal numérique intervient dans de nombreux domaines : transmission
du signal, traitement du son, des images, traitement des données scientiques, etc. Ce
document explique le principe de fonctionnement d'un ltre numérique et donne quelques
exemples de ltres.
On se limite au cas des signaux unidimensionnels. On considère principalement des
ltres à réponse impulsionnelle nie (ltres RIF). Une étude plus détaillée des ltres RIF
est présentée dans Filtres à réponse impulsionnelle nie. Voir aussi Exemples de ltres
RIF.
On verra aussi un exemple de ltre à réponse impulsionnelle innie, le ltre intégra-
teur.

2. Signal échantillonné
Lorsqu'un signal analogique x(t) est échantillonné avec une période d'échantillonnage
Te , on obtient la suite suivante :

xn = x(nTe ) (1)
où n est un entier. La suite xn constitue un signal numérique.
Lorsque le spectre de x(t) a une fréquence maximale fmax , la condition de Nyquist-
Shannon permet de s'assurer que toute l'information du signal analogique est présente
dans le signal échantillonné : la fréquence d'échantillonnage fe = 1/Te doit être supérieure
à 2fmax .

3. Filtres à réponse impulsionnelle nie


3.a. Dénition
Un ltre numérique à réponse impulsionnelle nie ([1]) réalise le calcul d'un signal de
sortie yn en fonction du signal d'entrée xn de la manière suivante :
N −1
(2)
X
yn = bk xn−k
k=0

La valeur de la sortie à l'instant n est donc une combinaison linéaire des valeurs de
l'entrée à l'instant n et aux N − 1 instants antérieurs. Si les valeurs de l'entrée ne sont
disponibles qu'à partir de l'indice 0, la première valeur de la sortie est yN −1 , puisqu'il
faut au moins N éléments en entrée pour faire le calcul. Pour voir la signication des
coecients hk , considérons en entrée une impulsion, dénie par :

x0 = 1 (3)
xn = 0 (n 6= 0) (4)
Frédéric Legrand Licence Creative Commons 2

La sortie est alors :

yn = 0 (n < 0) (5)
y 0 = b0 x 0 = b0 (6)
y 1 = b1 x 0 = b1 (7)
··· (8)
yN −1 = bN −1 x0 = bN −1 (9)
yn = 0 (n ≥ N ) (10)
Les coecients b0 , b1 ...bN −1 constituent donc la réponse impulsionnelle du ltre.
L'opération dénie par la relation (2) est un produit de convolution. La sortie est
donc le produit de convolution de l'entrée par la réponse impulsionnelle. Ce type de
ltrage est aussi appelé ltrage par convolution.
Le produit de convolution peut être calculé avec la fonction scipy.signal.convolve,
en suivant la syntaxe suivante :
y = scipy.signal.convolve(x,b,mode='valid')

L'option mode='valid' permet de limiter le calcul aux points valides : le premier


point calculé (yN −1 dans la relation (2)) est obtenu par combinaison des N premiers
points du tableau x, le dernier point par combinaison des N derniers points. Le tableau
y a ainsi N − 1 points de moins que le tableau x.
On peut aussi utiliser l'option mode='same' qui conserve le nombre de points. Pour
ce faire, les N/2 premiers points sont calculés en utilisant une partie de la réponse
impulsionnelle, de même pour les N/2 derniers points.

3.b. Réponse fréquentielle


Considérons un signal d'entrée sinusoïdal, de fréquence f et d'amplitude unité :

xn = exp(i2π f nTe ) (11)


Le signal en sortie s'écrit :

N −1
(12)
X
yn = bk xn−k
k=0
N −1
(13)
X
= exp(i2π n f Te ) bk exp(−i2π f k Te )
k=0
=xn H(Z) (14)
avec :

N −1
(15)
X
H(Z) = bk Z −k
k=0
Z = exp(i2π f Te ) (16)
Frédéric Legrand Licence Creative Commons 3

On voit donc que la sortie est dans ce cas proportionnelle à l'entrée. H(z) est la fonction
de transfert en Z , analogue à la fonction de transfert des signaux continus. La réponse
fréquentielle est nalement :
N −1
(17)
X
Hf (f ) = bk exp(−i2π k f Te )
k=0

On remarque que cette réponse fréquentielle est en fait une fonction de f Te , qui est
le rapport de la fréquence du signal sinusoïdal sur la fréquence d'échantillonnage fe =
1/Te . D'après le théorème de Shannon, ce rapport doit être inférieur à 1/2. La réponse
fréquentielle est de période fe mais en pratique on la trace sur l'intervalle [0, fe /2].
La fonction suivante renvoit le module et l'argument de la réponse fréquentielle. La
réponse impulsionnelle et le nombre de points de la courbe sont fournis en argument.

import math
import cmath
import numpy

def reponseFreq(b,nf):
f = numpy.arange(0.0,0.5,0.5/nf)
g = numpy.zeros(f.size)
phi = numpy.zeros(f.size)
for m in range(f.size):
H = 0.0
for k in range(b.size):
H += b[k]*cmath.exp(-1j*2*math.pi*k*f[m])
g[m] = abs(H)
phi[m] = cmath.phase(H)
return (f,g,numpy.unwrap(phi))

On peut aussi utiliser la fonction scipy.signal.freqz.

4. Exemples
4.a. Synthèse d'un signal bruité
Pour tester les diérents ltres, nous allons utiliser un signal périodique auquel on
ajoute un bruit gaussien :
On commence par dénir le signal par ses harmoniques :

def signal(t):
return 1.0*math.cos(2*math.pi*t) \
+0.5*math.cos(3*2*math.pi*t+math.pi/3) \
+0.2*math.cos(5*2*math.pi*t+math.pi/5) \
+0.2*math.cos(10*2*math.pi*t)
Frédéric Legrand Licence Creative Commons 4

Au moment d'échantillonner ce signal, on ajoute un bruit gaussien :

import numpy
import random
from matplotlib.pyplot import *

fe=500
T=2.0
N=int(fe*T)
te = T/N
x = numpy.zeros(N)
t = numpy.zeros(N)
sigma = 0.05
for k in range(N):
t[k] = k*te
x[k] = signal(t[k])+random.gauss(0,sigma)
figure(figsize=(10,4))
plot(t,x)

On calcule aussi le spectre du signal discret au moyen de sa transformée de Fourier


discrète :

from numpy.fft import fft


tfd = fft(x)
freq = numpy.arange(x.size)*1.0/T
spectre = 10*numpy.log10(numpy.absolute(tfd))
spectre = spectre-spectre.max()
figure(figsize=(10,4))
plot(freq,spectre,'r')
axis([0,fe/2,spectre.min(),spectre.max()])
xlabel("f")
ylabel("AdB")
grid()
Frédéric Legrand Licence Creative Commons 5

On repère sur ce spectre les raies correspondant aux harmoniques dénies dans le signal.
Le bruit se manifeste à toute fréquence (c'est un bruit blanc). La fréquence d'échan-
tillonnage a été choisie très largement supérieure aux fréquences utiles du signal. On
parle dans ce cas de sur-échantillonnage. Cette condition est nécessaire pour eectuer un
ltrage ecace du bruit.

4.b. Filtre moyenneur


Ce ltre calcule la moyenne arithmétique de deux valeurs consécutives :
1
yn = (xn + xn−1 ) (18)
2
On trace le module et l'argument de sa réponse fréquentielle :

b = numpy.array([0.5,0.5])
(f,g,phi)=reponseFreq(b,100)
figure(figsize=(6,4))
plot(f,g)
xlabel('f/fe')
ylabel('|H|')
axis([0,0.5,0,1])
grid()
Frédéric Legrand Licence Creative Commons 6

figure(figsize=(6,4))
plot(f,phi)
xlabel('f/fe')
ylabel('Phase')
axis([0,0.5,-3,3])
grid()

Le ltre moyenneur est un ltre passe-bas (peu sélectif). Le déphasage varie linéairement
avec la fréquence.
Frédéric Legrand Licence Creative Commons 7

Pour obtenir le signal discret ltré, il sut d'eectuer la convolution avec la réponse
impulsionnelle. Le premier point calculé correspond à l'instant Te , c'est pourquoi l'échelle
de temps doit être modiée.

import scipy.signal
b = np.array([0.5,0.5])
y = scipy.signal.convolve(x,b,mode='valid')
ny = y.size
ty = te+np.arange(ny)

figure(figsize=(10,4))
plot(ty,y)
xlabel('t')
ylabel('y')
grid()

On voit que le bruit est légèrement réduit par ce ltre passe-bas. Pour augmenter l'ef-
cacité du ltrage, il faut augmenter l'ordre N du ltre (la longueur de la réponse
impulsionnelle).
Une première idée consiste à calculer une moyenne sur un plus grand nombre de
points. Voici la réponse fréquentielle pour une moyenne sur 10 points :

b = numpy.ones(10)/10.0
(f,g,phi)=reponseFreq(b,100)
figure(figsize=(6,4))
plot(f,g)
xlabel('f/fe')
ylabel('|H|')
axis([0,0.5,0,1])
grid()
Frédéric Legrand Licence Creative Commons 8

Il y a des rebonds dans la bande atténuée qui en font un très mauvais ltre. Nous allons
voir comment modier les coecients de la réponse impulsionnelle pour obtenir une
meilleure réponse.

4.c. Filtre passe-bas gaussien


Pour réduire le bruit d'un signal, on utilise le plus souvent un ltre gaussien, dont la
réponse impulsionnelle est une gaussienne.
La gaussienne est la fonction dénie par :
t2
 
1
f (t) = √ exp − 2 (19)
σ 2π 2σ
On choisit une longueur de réponse impulsionnelle impaire N = 2P + 1, ce qui permet
de centrer la gaussienne sur l'indice k = P . La réponse impulsionnelle est :

bk = f (k − P ) (20)
Si P est donné, on choisit σ pour avoir une valeur  très faible sur les bords :
P
σ=√ (21)
−2 ln 
La réponse impulsionnelle est normalisée pour que la somme de ses coecients soit égale
à 1. Voici un exemple :

P_gauss=10
b_gauss = numpy.zeros(2*P_gauss+1)
epsilon=0.01
sigma=P_gauss/math.sqrt(-2.0*math.log(epsilon))
som = 0.0
for k in range(2*P_gauss+1):
Frédéric Legrand Licence Creative Commons 9

b_gauss[k] = math.exp(-(k-P_gauss)**2/(2*sigma**2))
som += b_gauss[k]
b_gauss = b_gauss/som
(f,g,phi)=reponseFreq(b_gauss,500)
figure(figsize=(6,4))
plot(f,g)
xlabel('f/fe')
ylabel('|H|')
axis([0,0.5,0,1])
grid()

La réponse fréquentielle d'un ltre gaussien est aussi une gausienne. Voyons la phase :

figure(figsize=(6,4))
plot(f,phi)
xlabel('f/fe')
ylabel('Phase')
grid()
Frédéric Legrand Licence Creative Commons 10

Dans la bande passante (et même au delà), la phase varie linéairement avec la fréquence.
Il y a donc un décalage temporel constant entre l'entrée et la sortie, égal à P Te .
Voyons l'eet sur le signal bruité déni plus haut. L'échelle de temps est calculée de
manière à annuler le décalage entre l'entrée et la sortie (cela n'est pas faisable dans un
ltre fonctionnant en temps réel).

y = scipy.signal.convolve(x,b_gauss,mode='valid')
ny = y.size
ty = np.arange(ny)+P_gauss*te
figure(figsize=(10,4))
plot(t,x,'b')
plot(ty,y,'r')
xlabel('t')
ylabel('y')
grid()
Frédéric Legrand Licence Creative Commons 11

Le bruit est éliminé sans que la forme du signal ne soit modiée. On remarque que P
points au début et à la n du signal ne sont pas traités. Dans un ltre temps-réel, cela
ne pose pas de problème car les 2P points non ltrés sont au début du signal.

4.d. Filtre passe-bas en sinus cardinal


Le ltre passe-bas gaussien est ecace pour éliminer le bruit. Sur sa réponse fréquen-
tielle, on voit en eet qu'il élimine totalement le bruit sur la bande de fréquence entre
0.15 et 0.5. En revanche, il n'est pas très sélectif. Une meilleure sélectivité signie une
décroissance plus rapide du gain après la coupure. Cela peut être obtenu par un ltre
dont la réponse impulsionnelle est un sinus cardinal. On dénit une fréquence de coupure
fc et on pose :
fc
a= (22)
fe
qui doit être inférieur à 1/2. La réponse impulsionnelle est dénie par :

bk = 2a sinc(2π(k − P )a) (23)

P=10
b = numpy.zeros(2*P+1)
def sinc(u):
if u==0:
return 1.0
else:
return math.sin(u)/u
a=0.15
for k in range(2*P+1):
b[k] = 2*a*sinc(2*math.pi*(k-P)*a)
(f,g,phi)=reponseFreq(b,500)
figure(figsize=(6,4))
plot(f,g)
xlabel('f/fe')
ylabel('|H|')
axis([0,0.5,0,1.2])
grid()
Frédéric Legrand Licence Creative Commons 12

Ce ltre est beaucoup plus sélectif que le ltre gaussien, mais a l'inconvénient de com-
porter des ondulations dans la bande passante et dans la bande atténuée.

figure(figsize=(6,4))
plot(f,phi)
xlabel('f/fe')
ylabel('Phase')
grid()
Frédéric Legrand Licence Creative Commons 13

La phase est bien linéaire dans la bande passante. Voyons l'eet sur le signal :

y = scipy.signal.convolve(x,b,mode='valid')
ny = y.size
ty = np.arange(ny)+P*te
figure(figsize=(10,4))
plot(t,x,'b')
plot(ty,y,'r')
xlabel('t')
ylabel('y')
grid()

Ce ltre est moins ecace que le ltre gaussien pour réduire le bruit, en raison de la
présence de rebonds dans la bande atténuée.

4.e. Filtre dérivateur


La dérivation est une opération courante en traitement du signal. Elle est utilisée pour
calculer la vitesse d'un phénomène, ou détecter des variations rapides dans un signal ou
dans une image.
La manière la plus simple de dénir une dérivation est :
1
yn = (xn − xn−1 ) (24)
Te
Sa réponse impulsionnelle est (1, −1). Voyons sa réponse fréquentielle :

b = numpy.array([1.0,-1.0])
(f,g,phi)=reponseFreq(b,500)
figure(figsize=(6,4))
plot(f,g)
xlabel('f/fe')
ylabel('|H|')
grid()
Frédéric Legrand Licence Creative Commons 14

figure(figsize=(6,4))
plot(f,phi)
xlabel('f/fe')
ylabel('Phase')
grid()

Le ltre dérivateur est passe-haut. Voyons son eet sur le signal bruité :
y_deriv = scipy.signal.convolve(x,b,mode='valid')
Frédéric Legrand Licence Creative Commons 15

ny = y_deriv.size
ty = te+np.arange(ny)
figure(figsize=(10,4))
plot(ty,y_deriv,'r')
xlabel('t')
ylabel('y')
grid()

Le bruit est tellement amplié par le ltre dérivateur qu'il couvre complètement le signal
utile. Cela est dû au caractère passe-haut du dérivateur. Pour dériver un signal bruité, il
faut impérativement eectuer un ltrage passe-bas au préalable, par exemple un ltrage
gaussien :

y = scipy.signal.convolve(x,b_gauss,mode='valid')
y_deriv = scipy.signal.convolve(y,b,mode='valid')
ny = y_deriv.size
ty = np.arange(ny)+(P_gauss+1)*te
figure(figsize=(10,4))
plot(ty,y_deriv,'r')
xlabel('t')
ylabel('y')
grid()
Frédéric Legrand Licence Creative Commons 16

5. Filtre intégrateur
Un ltre intégrateur peut être obtenu par la relation suivante :
xn + xn−1
yn = yn−1 + Te (25)
2
Voir Filtres intégrateur et dérivateur pour la justication de cette relation.
Cette relation est diérente de la relation (2) d'un ltre à réponse impulsionnelle
nie : la sortie à l'instant n dépend de l'état de la sortie à l'instant n − 1. Ce type de
ltre est appelé ltre récursif. Pour l'appliquer à un signal, il faut xer une condition
initiale, ici la valeur y0 .
Appliquons ce ltre à l'impulsion unité, avec la condition initiale y0 = 0. On obtient :

x0 1
y1 = y0 + = (26)
Te Te
1
y2 = y1 = (27)
Te
1
y3 = y2 = (28)
Te
··· (29)
La réponse impulsionnelle est un échelon. C'est un exemple de réponse impulsionnelle
innie. Un ltre est stable si sa réponse impulsionnelle tend vers zéro. Le ltre intégrateur
est donc instable, mais sa réponse impulsionnelle est tout de même bornée (stabilité
marginale).
La réponse fréquentielle de ce ltre peut être obtenue en suivant la même méthode
qu'en 3.b. On obtient ainsi la fonction de transfert en Z :
Te 1 + Z −1
H(Z) = (30)
2 1 − Z −1
et la réponse fréquentielle se déduit en posant :

Z = exp(i2π f Te ) (31)
D'une manière générale, la fonction de transfert en Z d'un ltre linéaire est une frac-
tion rationnelle en Z −1 . La fonction scipy.signal.freqz permet d'obtenir la réponse
fréquentielle d'un ltre déni par sa fonction de transfert en Z :

b=[0.5,0.5]
a=[1,-1]
[w,h] = scipy.signal.freqz(b,a)
figure(figsize=(6,4))
plot(w/(2*math.pi),numpy.abs(h))
xlabel('f/fe')
ylabel('|H|')
grid()
Frédéric Legrand Licence Creative Commons 17

L'échelle logarithmique (diagramme de Bode) est préférable :


figure(figsize=(6,4))
plot(numpy.log10(w/(2*math.pi)),20*numpy.log10(numpy.abs(h)))
xlabel('log(f/fe)')
ylabel('GdB')
grid()

Voyons l'application de ce ltre au signal bruité. La fonction scipy.signal.lfilter


permet de le faire. La condition initiale par défaut est nulle.
Frédéric Legrand Licence Creative Commons 18

y = scipy.signal.lfilter(b,a,x)
figure(figsize=(10,4))
plot(t,y)
xlabel('t')
ylabel('y')
grid()

Au contraire du ltre dérivateur, le ltre intégrateur réduit le bruit, en raison de son


caractère passe-bas.

Références
[1] M. Bellanger, Traitement numérique du signal, (Dunod, 1998)

Vous aimerez peut-être aussi