Table des matières Python

Filtres à réponse impulsionnelle finie

1. introduction

Ce document constitue une introduction aux filtres à réponse impulsionnelle finie (filtres RIF). La méthode de conception d'un filtre RIF (à phase linéaire) par série de Fourier est introduite à travers un exemple.

Les filtres RIF sont aussi des filtres de convolution car leur réalisation dans le domaine temporel se fait par un produit de convolution entre le signal échantillonné et la réponse impulsionnelle du filtre.

2. Signal discret et transformée en Z

Soit x(t) un signal continu échantillonné avec une période Te. Le signal discret correspondant est la suite

xn=x(nTe)(1)

Considérons la transformée de Fourier du signal continu :

X˜(f)=-+x(t)exp(-i2πft)dt(2)

Une approximation de la transformée de Fourier est obtenue à partir du signal échantillonné par la méthode des rectangles :

X˜(f)n=-n=+xnexp(-i2πfnTe)(3)

Posons :

X(f)=n=-n=+xnexp(-i2πfnTe)(4)

X(f) est la transformée de Fourier à temps discret (TFTD) du signal discret. Lorsque la somme se fait sur un nombre N d'échantillons correspondant à une durée T, on obtient la transformée de Fourier discrète (TFD). Voici la définition de la TFD :

Xk=n=0Nxnexp(-i2πknN)(5)

Alors que la TFTD est définie pour toute fréquence f (nombre réel), le terme Xk de la TFD correspond à la fréquence fk=kT . Remarquons que la TFTD n'est pas calculable numériquement puisqu'elle se définit par une somme infinie. En d'autres termes, la TFD est la réalisation de la TFTD.

On pose :

Z=exp(i2πfTe)(6)

De manière équivalente à la TFTD, il est intéressant de définir la transformée en Z ([1],[2]) du signal discret :

X(Z)=n=-n=+xnZ-n(7)

La transformée en Z est ainsi l'équivalent pour les signaux discrets de la transformée de Fourier pour les signaux continus. En reportant l'expression (6) de Z dans la transformée en Z, on retrouve en effet la transformée de Fourier à temps discret (4). Cependant, la transformée en Z ne se limite pas à un moyen d'obtenir la TFTD car sa définition est étendue pour tout nombre complexe Z (sous réserve de convergence de la série) et non pas seulement pour les nombres complexes du cercle unité.

3. Système linéaire causal à réponse impulsionnelle finie

3.a. Définition

Un système linéaire causal à réponse impulsionnelle finie transforme un signal discret xn en un signal yn par la relation :

yn=k=0M-1 hkxn-k(8)

où les M coefficients hk sont des constantes réelles. La valeur du signal discret yn à l'instant n est donc obtenue par une combinaison linéaire des N valeurs précédentes du signal xn. Ce système est dit causal car l'état de la sortie ne dépend que des états antérieurs de l'entrée.

L'impulsion unité est le signal discret défini par :

x0=1(9)xn=0sin0(10)

Pour une impulsion en entrée, le système linéaire fournit en sortie le signal :

yn=hn(11)

C'est pourquoi la suite hn est appelée la réponse impulsionnelle (discrète) du système. La réponse impulsionnelle est ici finie puisqu'elle comporte M valeurs.

On remarque que le signal de sortie est le produit de convolution du signal d'entrée par la réponse impulsionnelle.

3.b. Fonction de transfert en Z

Ce système constitue un filtre à réponse impulsionnelle finie (filtre RIF). Considérons alors la transformée en Z du signal de sortie :

Y(Z)=n=-n=+k=0M-1 hkxn-kZ-n(12)=n'=-n'=+k=0M-1 hkxn'Z-n'Z-k(13)=(k=0M-1hkZ-k)n'=-n'=+xn'Z-n'(14)=H(Z)X(Z)(15)

H(Z) est la fonction de transfert en Z définie par :

H(Z)=k=0M-1hkZ-k(16)

La transformée en Z du signal de sortie est donc le produit de la transformée en Z du signal d'entrée par la transformée en Z de la réponse impulsionnelle (les transformées de Fourier des signaux continus vérifient un théorème analogue).

3.c. Réponse fréquentielle

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

xn=exp(i2πfnTe)(17)

Le signal en sortie s'écrit :

yn=k=0M-1 hkxn-k(18)=exp(i2πnfTe)k=0M-1hkexp(-i2πfkTe)(19)=exp(i2πnfTe)H(Z)(20)=xnH(Z)(21)

avec

Z=exp(i2πfTe)(22)

La réponse fréquentielle du filtre se déduit donc de la fonction de transfert en Z :

H(f)=k=0M-1hkexp(-i2πkfTe)(23)

On remarque que cette réponse fréquentielle est une fonction de fTe, qui est le rapport de la fréquence du signal sinusoïdal par la fréquence d'échantillonnage. D'après le théorème de Shannon, ce rapport doit être inférieur à 1/2.

En remarquant que la réponse fréquentielle est de périodicité fe=1/Te (la fréquence d'échantillonnage), considérons les M fréquences suivantes :

fn=nMTe(24)

On a alors

H(fn)=k=0M-1hkexp(-i2πknM)(25)

La réponse fréquentielle pour ces fréquences est donc obtenue par transformée de Fourier discrète (TFD) de la réponse impulsionnelle. Il s'en suit que la réponse impulsionnelle s'obtient par transformée de Fourier discrète inverse de la réponse fréquentielle :

hk = 1Mn=0M-1H(fn)exp(i2πnkM)(26)

3.d. Exemple

Considérons comme exemple simple un filtre qui réalise la moyenne arithmétique de deux valeurs consécutives de l'entrée :

yn=12(xn+xn-1)(27)

Sa fonction de transfert en Z est :

H(Z)=12(1+Z-1)(28)

Voyons le tracé sa réponse fréquentielle avec Python :

import numpy as np
from matplotlib.pyplot import *
def H(Z):
    return 0.5*(1.0+Z**(-1))
def Hf(f):
    return H(np.exp(1j*2*np.pi*f))
f = np.arange(start=0.0,stop=0.5,step=0.01)
figure(figsize=(6,4))
plot(f,np.absolute(Hf(f)),'r')
xlabel('f/fe')
ylabel('|H|')
axis([0,0.5,0,1])
grid()
                
figAfigA.pdf
phase = np.unwrap(np.angle(Hf(f)))
figure(figsize=(6,4))
plot(f,phase,'r')
xlabel('f/fe')
ylabel('Phase')
axis([0,0.5,-3,3])
grid()
                
figBfigB.pdf

Le filtre moyenneur est un filtre passe-bas (assez peu sélectif). Le déphasage varie linéairement avec la fréquence. Cela se vérifie sur l'expression suivante de la réponse fréquentielle :

H(f)=cos(πfTe)exp(-iπfTe)(29)

Comme exemple d'action de ce filtre, considérons un signal sinusoïdal auquel on ajoute du bruit :

import numpy.random as random
def signal(t):
    return np.cos(2*np.pi*t)
fe = 100
te=1.0/fe
t = np.arange(start=0.0,stop=5.0,step=te)
x = signal(t) +0.05*random.uniform(-1.0,1.0,len(t))
figure(figsize=(10,4))
plot(t,x)
xlabel('t')
ylabel('x')
axis([0,5,-1.5,1.5])
grid()
                
figCfigC.pdf

Pour obtenir le signal discret filtré, il faut effectuer la convolution avec la réponse impulsionnelle. La fonction scipy.signal.convolve permet de le faire. On utilise l'option mode='valid' , qui restreint le calcul aux points valides. Ainsi le premier point calculé est la moyenne des deux premiers et le dernier point est la moyenne des deux derniers. Le premier point calculé correspond à l'instant Te, c'est pourquoi l'échelle de temps doit être modifiée.

import scipy.signal
h = np.array([0.5,0.5])
y = scipy.signal.convolve(x,h,mode='valid')
ty = np.arange(len(y))*te
figure(figsize=(10,4))
plot(ty,y)
xlabel('t')
ylabel('y')
axis([0,5,-1.5,1.5])
grid()

                
figDfigD.pdf

L'effet de ce filtre passe-bas et une petite réduction du bruit.

4. Filtres RIF à phase linéaire

4.a. Définition et propriétés

Pour un filtre à phase linéaire, le déphasage est une fonction linéaire de la fréquence. La réponse fréquentielle a donc la forme suivante :

H(f)=G(f)exp(-iφ(f))(30)φ(f)=2πfτ(31)

La phase d'une composante de fréquence f devient en sortie du filtre :

2πft-φ(f)=2πf(t-τ)(32)

Toutes les fréquences du signal subissent le même décalage τ en traversant le filtre. τ est le temps de propagation. Si toutes les composantes spectrales sont dans la bande passante, pour laquelle G(f)=1, on obtient :

y(t)=x(t-τ)(33)

La forme du signal n'est pas modifiée par le filtrage en bande passante. On définit le retard en nombre d'échantillons par P=feτ (P entier). La fonction de transfert en Z d'un filtre à phase linéaire dans la bande passante est donc :

Hbp(Z)=Z-P(34)

En isolant le terme contenant la phase, la réponse fréquentielle s'écrit d'après l'expression(15) :

H(f)=exp(-i2πfPTe)k=0M-1hkexp(-i2π(k-P)fTe)(35)

Après un changement de variable dans la somme, on en déduit l'expression du gain :

G(f)=k=-PM-1-Phk+Pexp(-i2πkfTe)(36)

On veut que G(f) soit réel. Il faut donc que la réponse impulsionnelle soit symétrique par rapport à l'indice P, c'est-à-dire :

hk+P=h-k+P(37)

On en déduit que M=2P+1. M est donc impair. Il est toutefois possible d'obtenir un filtre RIF à phase linéaire avec M pair en posant P-12=feτ .

4.b. Synthèse par série de Fourier

Soit la suite définie par :

gk=hk+P(-PkP)(38)

Elle vérifie la propriété g-k=gk. La réponse fréquentielle s'écrit

H(f)=exp(-i2πfτ)k=-PPgkexp(-i2πkfTe)=exp(-i2πfτ)G(f)(39)

En considérant la limite P , on obtient alors :

G(f)=k=-+gkexp(-i2πkfTe)(40)

On obtient un filtre à phase linéaire mais à réponse impulsionnelle infini (impossible à réaliser). Cette expression a toutefois l'avantage de montrer que les coefficients gk sont les coefficients de Fourier de la fonction G(f), fonction de période fe=1/Te. On a ainsi :

gk=1fe0feG(f)exp(i2πkffe)df(41)

Puisque g-k=gk, ces coefficients sont réels et la fonction G(f) est paire.

Pour une fonction de gain G(f) souhaitée, on calcule les coefficients de Fourier gk puis on tronque la série de Fourier à un rang P de manière à obtenir la réponse impulsionnelle finie définie par :

hk=gk-P(0k2P)(42)

Cette méthode revient à appliquer un fenêtrage rectangulaire aux coefficients de Fourier. La fonction G(f) effectivement obtenue après troncature est donc le produit de convolution de la fonction souhaitée par la transformée de Fourier de la fenêtre. En pratique, il est préférable d'appliquer d'autres types de fenêtres, comme les fenêtres de Hann ou de Hamming.

4.c. Exemple : filtre RIF passe-bas

La fonction de gain souhaitée est définie sur l'intervalle [-fe/2,fe/2] par :

G(f)=1si-fcffc(43)=0sifc<|f|fe2(44)

Les coefficients de Fourier de cette fonction sont :

gk=1kπsin(2πkfcfe)(45)

Le résultat peut s'exprimer avec la fonction sinus cardinal et ne dépend que du rapport de la fréquence de coupure sur la fréquence d'échantillonnage a=fcfe :

gk=2asinc(2πka)(46)

Considérons le cas d'une troncature par une fenêtre rectangulaire au rang P. On fixe la fréquence de coupure, le rang P de troncature et on calcule la réponse impulsionnelle :

a=0.1
P=20
def g(k):
    return 2*a*np.sinc(2*(k-P)*a)
N=2*P+1
liste_k = np.arange(start=0,stop=N)
h = g(liste_k)
figure(figsize=(10,4))
vlines(liste_k,[0],h,'r')
xlabel('k')
ylabel('hk')
grid()
                
figEfigE.pdf

La réponse fréquentielle est donnée par la fonction suivante :

def Hf(f):
    s = 0.0
    for k in range(N):
       s += h[k]*np.exp(-1j*2*np.pi*k*f)
    return s
                

On trace le gain et le déphasage :

f = np.arange(start=0.0,stop=0.5,step=0.001)

figure(figsize=(6,4))
plot(f,np.absolute(Hf(f)),'r')
xlabel('f/fe')
ylabel('|H|')
axis([0,0.5,0,2])
grid()
                
figFfigF.pdf
phase = np.unwrap(np.angle(Hf(f)))
figure(figsize=(6,4))
plot(f,phase,'r')
xlabel('f/fe')
ylabel('Phase')
grid()
                
figGfigG.pdf

On constate que la phase est bien linéaire dans la bande passante, mais le gain présente des ondulations très fortes. Dans la bande atténuée, il y a des discontinuités de π de la phase. Bien sûr, les différences par rapport à la fonction de transfert souhaitée sont dûes à la troncature de la réponse impulsionnelle.

Essayons une troncature par une fenêtre de Hann :

h = h*scipy.signal.hann(N)
figure(figsize=(6,4))
plot(f,np.absolute(Hf(f)),'r')
xlabel('f/fe')
ylabel('|H|')
axis([0,0.5,0,2])
grid()
                
figHfigH.pdf
phase = np.unwrap(np.angle(Hf(f)))
figure(figsize=(6,4))
plot(f,phase,'r')
xlabel('f/fe')
ylabel('Phase')
grid()
                
figIfigI.pdf

Les ondulations dans la bande passante et dans la bande atténuée sont considérablement réduites. La linéarité de la phase dans la bande passante est toujours assurée.

Lorsque ce filtre est utilisé dans un système de traitement du signal en temps réel, pour un signal situé dans la bande passante, la sortie présente un retard de τ=PTe par rapport à l'entrée. Cela se comprend aisément en remarquant que le point prépondérant pour le calcul de yn est xn-P (terme g0). Pour obtenir un filtre très sélectif, il faut augmenter P. Si le retard τ doit rester fixe, il faut parallèlement augmenter la fréquence d'échantillonnage.

Pour l'analyse de données expérimentales (nombre échantillons fini), le retard τ n'a aucune signification physique et peut être annulé de la manière suivante. Soit Q le nombre de points de l'échantillon. Le premier point de la sortie calculé est y2P puisqu'il faut au moins 2P+1 valeurs antérieures de xn pour appliquer la convolution. Ce premier point est en retard de τ=PTe par rapport à l'entrée. Annuler le retard revient donc à attribuer l'instant (2P-P)Te=PTe à ce premier point. Ainsi, le dernier point calculé yQ-2P correspond à l'instant (Q-P)Te. L'application du filtrage revient donc à enlever les P premiers instants et les P derniers instants de l'échantillon. Voyons cela sur le signal déjà défini plus haut.

y = scipy.signal.convolve(x,h,mode='valid')
N = len(y)
ty = P*te+te*np.arange(N)
figure(figsize=(10,4))
plot(t,x,'b')
plot(ty,y,'r')
xlabel('t')
ylabel('y')
axis([0,5,-1.5,1.5])
grid()
                
figJfigJ.pdf

Le signal traité est bien synchrone avec le signal d'origine; les P premiers et les P derniers points sont perdus. Pour un filtrage numérique en temps réel, la sortie serait retardée, c'est-à-dire décalée de τ=PTe vers la droite.

Remarquer l'utilisation de l'option mode='valid' qui restreint le calcul de la convolution aux points valides. Ainsi le premier point de la convolution est calculé à partir des 2P+1 premiers points de x et le dernier point est calculé à partir des 2P1+1 derniers points de x.

Références
[1]  M. Bellanger,  Traitement numérique du signal,  (Dunod, 1998)
[2]  E. Tisserand, J.F. Pautex, P. Schweitzer,  Analyse et traitement des signaux,  (Dunod, 2008)
Creative Commons LicenseTextes et figures sont mis à disposition sous contrat Creative Commons.