Filtres à réponse impulsionnelle finie

3.a. Définition

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

$$y_n=\sum _{k=0}^{M-1}{h}_k{x}_{n-k}\text{(8)}$$

où les M coefficients h_k sont des constantes réelles. La valeur du signal discret y_n à l'instant n est donc obtenue par une combinaison linéaire des N valeurs précédentes du signal x_n. 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 :

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

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

$$y_n=h_n\text{(11)}$$

C'est pourquoi la suite h_n 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 :

$$\begin{array}{ccc}\hfill Y\left(Z\right)& =\sum _{n=-\infty }^{n=+\infty }\sum _{k=0}^{M-1}{h}_k{x}_{n-k}{Z}^{-n}\hfill & \hfill \text{(12)}\\ \hfill & =\sum _{n\text{'}=-\infty }^{n\text{'}=+\infty }\sum _{k=0}^{M-1}{h}_k{x}_{n\text{'}}{Z}^{-n\text{'}}{Z}^{-k}\hfill & \hfill \text{(13)}\\ \hfill & =\left(\sum _{k=0}^{M-1}{h}_k{Z}^{-k}\right)\sum _{n\text{'}=-\infty }^{n\text{'}=+\infty }{x}_{n\text{'}}{Z}^{-n\text{'}}\hfill & \hfill \text{(14)}\\ \hfill & =H\left(Z\right)X\left(Z\right)\hfill & \hfill \text{(15)}\end{array}$$

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

$$H\left(Z\right)=\sum _{k=0}^{M-1}{h}_k{Z}^{-k}\text{(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é :

$$x_n=exp\left(i2\pi fn{T}_e\right)\text{(17)}$$

Le signal en sortie s'écrit :

$$\begin{array}{ccc}\hfill y_n=& \sum _{k=0}^{M-1}{h}_k{x}_{n-k}\hfill & \hfill \text{(18)}\\ \hfill =& exp\left(i2\pi nf{T}_e\right)\sum _{k=0}^{M-1}{h}_{k}exp(-i2\pi fk{T}_e)\hfill & \hfill \text{(19)}\\ \hfill =& exp\left(i2\pi nf{T}_e\right)H\left(Z\right)\hfill & \hfill \text{(20)}\\ \hfill =& x_{n}H\left(Z\right)\hfill & \hfill \text{(21)}\end{array}$$

avec

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

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

$$H\left(f\right)=\sum _{k=0}^{M-1}{h}_{k}exp(-i2\pi kf{T}_e)\text{(23)}$$

On remarque que cette réponse fréquentielle est une fonction de fT_e, 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é f_e=1/T_e (la fréquence d'échantillonnage), considérons les M fréquences suivantes :

$$f_n=\frac{n}{M{T}_e}\text{(24)}$$

On a alors

$$H\left(f_n\right)=\sum _{k=0}^{M-1}{h}_{k}exp(-i2\pi k\frac{n}{M})\text{(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 :

$$h_k=\frac{1}{M}\sum _{n=0}^{M-1}H\left(f_n\right)exp\left(i2\pi n\frac{k}{M}\right)\text{(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 :

$$y_n=\frac{1}{2}(x_n+x_{n-1})\text{(27)}$$

Sa fonction de transfert en Z est :

$$H\left(Z\right)=\frac{1}{2}(1+Z^{-1})\text{(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()
                

figA.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()
                

figB.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\left(f\right)=cos\left(\pi f{T}_e\right)exp(-i\pi f{T}_e)\text{(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()
                

figC.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 T_e, 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()

                

figD.pdf

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

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 :

$$\begin{array}{ccc}\hfill & H\left(f\right)=G\left(f\right)exp(-i\phi \left(f\right))\hfill & \hfill \text{(30)}\\ \hfill & \phi \left(f\right)=2\pi f\tau \hfill & \hfill \text{(31)}\end{array}$$

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

$$2\pi ft-\phi \left(f\right)=2\pi f(t-\tau )\text{(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\left(t\right)=x(t-\tau )\text{(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=f_eτ (P entier). La fonction de transfert en Z d'un filtre à phase linéaire dans la bande passante est donc :

$$H_{bp}\left(Z\right)=Z^{-P}\text{(34)}$$

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

$$H\left(f\right)=exp(-i2\pi fP{T}_e)\sum _{k=0}^{M-1}{h}_{k}exp(-i2\pi (k-P)f{T}_e)\text{(35)}$$

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

$$G\left(f\right)=\sum _{k=-P}^{M-1-P}{h}_{k+P}exp(-i2\pi kf{T}_e)\text{(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 :

$$h_{k+P}=h_{-k+P}\text{(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-\frac{1}{2}=f_e\tau$.

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

Soit la suite définie par :

$$g_k=h_{k+P}(-P\le k\le P)\text{(38)}$$

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

$$H\left(f\right)=exp(-i2\pi f\tau )\sum _{k=-P}^P{g}_{k}exp(-i2\pi kf{T}_e)=exp(-i2\pi f\tau )G\left(f\right)\text{(39)}$$

En considérant la limite $P\to \infty$, on obtient alors :

$$G\left(f\right)=\sum _{k=-\infty }^{+\infty }{g}_{k}exp(-i2\pi kf{T}_e)\text{(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 g_k sont les coefficients de Fourier de la fonction G(f), fonction de période f_e=1/T_e. On a ainsi :

$$g_k=\frac{1}{f_e}{\int }_0^{f_e}G\left(f\right)exp\left(i2\pi k\frac{f}{f_e}\right)df\text{(41)}$$

Puisque g_-k=g_k, 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 g_k puis on tronque la série de Fourier à un rang P de manière à obtenir la réponse impulsionnelle finie définie par :

$$h_k=g_{k-P}\left(0\le k\le 2P\right)\text{(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 [-f_e/2,f_e/2] par :

$$\begin{array}{ccc}\hfill G\left(f\right)& =1si-f_c\le f\le f_c\hfill & \hfill \text{(43)}\\ \hfill & =0si{f}_c<\left|f\right|\le \frac{f_e}{2}\hfill & \hfill \text{(44)}\end{array}$$

Les coefficients de Fourier de cette fonction sont :

$$g_k=\frac{1}{k\pi }sin\left(2\pi k\frac{f_c}{f_e}\right)\text{(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=\frac{f_c}{f_e}$ :

$$g_k=2asinc\left(2\pi ka\right)\text{(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()
                

figE.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()
                

figF.pdf

phase = np.unwrap(np.angle(Hf(f)))
figure(figsize=(6,4))
plot(f,phase,'r')
xlabel('f/fe')
ylabel('Phase')
grid()
                

figG.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()
                

figH.pdf

phase = np.unwrap(np.angle(Hf(f)))
figure(figsize=(6,4))
plot(f,phase,'r')
xlabel('f/fe')
ylabel('Phase')
grid()
                

figI.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 τ=PT_e par rapport à l'entrée. Cela se comprend aisément en remarquant que le point prépondérant pour le calcul de y_n est x_n-P (terme g₀). 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 y_2P puisqu'il faut au moins 2P+1 valeurs antérieures de x_n pour appliquer la convolution. Ce premier point est en retard de τ=PT_e par rapport à l'entrée. Annuler le retard revient donc à attribuer l'instant (2P-P)T_e=PT_e à ce premier point. Ainsi, le dernier point calculé y_Q-2P correspond à l'instant (Q-P)T_e. 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()
                

figJ.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 τ=PT_e 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.