Filtrage numérique

2.a. Principe

La figure suivante montre le fonctionnement d'un filtrage numérique par convolution.

Figure pleine page

Seuls les 6 derniers échantillons de l'entrée x_n sont représentés. La convolution se fait avec un noyau, qui dans cet exemple comporte trois coefficients b[0],b[1],b[2]. Le signal filtré est noté y_n. Pour obtenir le dernier échantillon du signal filtré, on effectue une combinaison linéaire des trois derniers échantillons du signal d'entrée :

$$y_n=b_0{x}_n+b_1{x}_{n-1}+b_2{x}_{n-2}$$

Plus généralement, si N est le nombre d'élements du noyau, l'échantillon n de la sortie se calcule en faisant une combinaison linéaire des N derniers échantillons de l'entrée :

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

Cette relation définit la convolution discrète, et permet donc de réaliser un filtrage par convolution. L'effet du filtre dépend évidemment des coefficients choisis pour le noyau.

Dans le schéma ci-dessus, y_n et x_n sont considérés comme simultanés. Dans un filtre fonctionnant en temps réel, y_n est calculé dès que x_n est disponible et le temps de calcul doit être inférieur à la période d'échantillonnage. Si le temps de calcul est faible par rapport à T_e, y_n et x_n sont pratiquement simultanés.

Le filtre de convolution le plus simple est le filtre moyenneur, qui réalise la moyenne arithmétique des deux derniers échantillons :

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

Le filtre suivant définit y_n comme la moyenne arithmétique des N derniers échantillons :

$$y_n=\frac{1}{N}\sum _{k=0}^{N-1}{x}_{n-k}$$

2.b. Réponse impulsionnelle

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

$$\begin{array}{cc}\hfill & x_0=1\hfill \\ \hfill & x_n=0\left(n\ne 0\right)\hfill \end{array}$$

La sortie est alors :

$$\begin{array}{cc}\hfill & y_n=0(n<0)\hfill \\ \hfill & y_0=b_0{x}_0=b_0\hfill \\ \hfill & y_1=b_1{x}_0=b_1\hfill \\ \hfill & \cdots \hfill \\ \hfill & y_{N-1}=b_{N-1}{x}_0=b_{N-1}\hfill \\ \hfill & y_n=0\left(n\ge M\right)\hfill \end{array}$$

Les coefficients b₀,b₁... b_N-1 (le noyau du filtre) constituent donc la réponse impulsionnelle du filtre.

Un filtre de convolution est dit à réponse impulsionnelle finie (RIF) car le nombre de termes de la réponse impulsionnelle est fini.

Les filtres numériques RIF sont toujours stables : si l'entrée reste bornée, la sortie l'est aussi.

2.c. Réponse fréquentielle

Pour obtenir la réponse fréquentielle théorique du filtre, on considère un signal d'entrée sinusoïdal, de fréquence f et d'amplitude unité. L'échantillonnage à la période T_e donne le signal numérique suivant :

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

Le signal en sortie s'écrit :

$$\begin{array}{cc}\hfill y_n=& \sum _{k=0}^{N-1}{b}_k{x}_{n-k}\hfill \\ \hfill =& exp\left(i2\pi nf{T}_e\right)\sum _{k=0}^{N-1}{b}_{k}exp(-i2\pi fk{T}_e)\hfill \\ \hfill =& x_{n}H\left(Z\right)\hfill \end{array}$$

avec :

$$\begin{array}{ccc}\hfill & H\left(Z\right)=\sum _{k=0}^{N-1}{b}_k{Z}^{-k}\hfill & \hfill \text{(2)}\\ \hfill & Z=exp\left(i2\pi f{T}_e\right)\hfill & \hfill \text{(3)}\end{array}$$

On voit donc que la sortie (écrite sous forme complexe) est 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 finalement :

$$H_f\left(f\right)=\sum _{k=0}^{N-1}{b}_{k}exp(-i2\pi kf{T}_e)\text{(4)}$$

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

La représentation graphique de la réponse fréquentielle consiste à tracer le module de H_f et son argument en fonction de la fréquence sur l'intervalle [0,f_e/2].

2.d. Filtrage d'un signal en mémoire

Nous n'allons pas effectuer un filtrage en temps réel, mais le filtrage d'un signal numérique stocké en mémoire après la conversion analogique/numérique. On remarque que le filtrage par convolution nécessite au moins N termes en entrée pour s'appliquer.

La fonction suivante effectue le fitrage par convolution :

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
             

Les N-1 premiers échantillons de la sortie sont nuls car le filtrage ne commence qu'à partir de l'échantillon d'indice N-1. Afin de simuler le fonctionnement d'un filtre temps réel, l'échantillon y_n de la sortie est considéré comme simultané avec l'échantillon x_n de l'entrée.

2.e. Filtre moyenneur

Considérons le filtre suivant, qui définit y_n comme la moyenne des 10 derniers échantillons de l'entrée :

$$y_n=\frac{1}{10}\sum _{k=0}^9{x}_{n-k}$$

La suite de ces coefficients (qui est aussi sa réponse impulsionnelle) est placée dans un tableau :

import numpy
b = numpy.ones(10)*1/10
                

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

import scipy.signal
w,H=scipy.signal.freqz(b)
                

Le tableau w contient les pulsations (500 points par défaut) et le tableau H contient les valeurs complexes de la réponse fréquentielle. Voici comment tracer le gain et le déphasage en fonction de la fréquence :

from matplotlib.pyplot import *               
figure()
subplot(211)
plot(w/(2*numpy.pi),numpy.absolute(H))
ylabel("G")
grid()
subplot(212)
plot(w/(2*numpy.pi),numpy.unwrap(numpy.angle(H)))
xlabel("f/fe")
ylabel("déphasage")
grid()
                
                

figA0.pdf

Le filtre moyenneur est donc un filtre passe-bas. Pour les filtres numériques, il est d'usage de définir la fréquence de coupure comme la fréquence pour laquelle le gain vaut 1/2. La fréquence de coupure de ce filtre est donc 0,06f_e. Le fait d'avoir multiplié la moyenne par 1/10 permet d'obtenir un gain de 1 à très basse fréquence. Le déphasage dans la bande passante varie parfaitement linéairement avec la fréquence, ce qui est la propriété recherchée pour le filtrage des signaux périodiques, car elle implique que le retard est constant.

Le filtre moyenneur est rarement utilisé car les caractéristiques de sa bande atténuante sont mauvaises : le gain comporte des rebonds dans la bande atténuante, ce qui fait que la décroissance du gain est très lente. Par exemple, à une fréquence 0,25f_e, soit 4 fois la fréquence de coupure, le gain est de 0,14.

Nous allons voir qu'il est possible de définir un filtre passe-bas avec une décroissance du gain dans la bande atténuante beaucoup plus rapide.

2.f. Filtre passe-bas

Pour définir un filtre passe-bas, il faut tout d'abord calculer le rapport de la fréquence de coupure sur la fréquence d'échantillonnage (qui doit être inférieur à 1/2) :

$$a=\frac{f_c}{f_e}\text{(5)}$$

On appelle filtre idéal un filtre qui laisse passer sans atténuation toutes les sinusoïdes de fréquences inférieures à la fréquence de coupure, et qui annule toutes les sinusoïdes de fréquences supérieures à la coupure. Son déphasage doit varier linéairement avec la fréquence dans la bande passante. Un filtre passe-bas idéal est irréalisable car sa réponse impulsionnelle est infinie. Elle est donnée par :

$$g\left(k\right)=2a\frac{sin\left(2\pi ka\right)}{2\pi ka}=2asinc\left(2ka\right)\text{(6)}$$

La fonction sinc est la fonction sinus cardinal définie par $sinc\left(u\right)=sin\left(\pi u\right)/\left(\pi u\right)$. Voici par exemple la réponse impulsionnelle pour a=0,1 (en fait une partie de cette réponse) :

a=0.1
k=numpy.arange(-50,50)
g=2*a*numpy.sinc(2*k*a)
figure(figsize=(12,4))
stem(k,g)
xlabel("k")
ylabel("g")
grid()

                        

figA1.pdf

Cette réponse impulsionnelle est infinie et on ne peut donc pas construire un filtre de convolution idéal. La réponse impulsionnelle est rendue finie par troncature. Soit P l'indice de troncature. La réponse impulsionnelle finie comporte alors N=2P+1 termes, correspondant aux indices k allant de -P à P. Par ailleurs, la réponse impulsionnelle permettant d'obtenir le noyau du filtre doit être définie pour des indices positifs ou nuls, ce qui s'obtient en décalant les indices de +P. Finalement, la réponse impulsionnelle est définie par :

$$b_k=g(k-P)\text{(7)}$$

l'indice k variant de 0 à 2P.

L'indice de troncature devra être d'autant plus grand que a est faible. Il faut donc éviter une fréquence d'échantillonnage trop grande par rapport à la fréquence de coupure (le sur-échantillonnage n'est pas toujours souhaitable en filtrage numérique). Un critère simple est Pa>1. Pour un coefficient a donné, la sélectivité du filtre augmente avec P. À la limite $P\to \infty$, on obtient un filtre idéal.

Voici un exemple de filtre passe-bas et sa réponse impulsionnelle :

a=0.1
P = 10
b = numpy.zeros(2*P+1)
n = range(2*P+1)

        
for k in range(2*P+1):
    b[k] = 2*a*numpy.sinc(2*(k-P)*a)

figure(figsize=(8,4)) 
stem(n,b)
xlabel("n")
ylabel("b")
grid()                               

                         

figA2.pdf

On voit que le maximum de la réponse impulsionnelle (le maximum du sinus cardinal) est décalé de P échantillons par rapport à l'échantillon de sortie calculé. Cela signifie que le signal de sortie sera retardé de P échantillons.

Voici la réponse fréquentielle de ce filtre :

w,H=scipy.signal.freqz(b)

figure()
subplot(211)
plot(w/(2*numpy.pi),numpy.absolute(H))
ylabel("G")
grid()
subplot(212)
plot(w/(2*numpy.pi),numpy.unwrap(numpy.angle(H)))
xlabel("f/fe")
ylabel("déphasage")
grid()
                         

figA3.pdf

Le gain à la fréquence de coupure est égal à 1/2, et non pas à $1/\sqrt{2}$ comme il est d'usage pour un filtre analogique.

Pour tracer le déphasage, nous avons utilisé la fonction numpy.unwrap, qui permet d'éviter que l'angle soit ramené dans l'intervalle [-π,π]. En effet, le retard de la sortie par rapport à l'entrée peut dépasser la moitié de la période. On remarque que le déphasage dans la bande passante est parfaitement linéaire par rapport à la fréquence. En revanche, le gain dans la bande passante présente un maximum alors qu'il est préférable d'avoir un gain décroissant. Pour éliminer cet effet, on doit multiplier la réponse impulsionnelle par une fenêtre, c'est-à-dire une fonction qui réduit progressivement les valeurs de la réponse lorsqu'on s'approche du bord de l'intervalle. Voici par exemple la fenêtre de Hamming :

import scipy.signal
fen = scipy.signal.get_window("hamming",b.size)
figure(figsize=(8,4))
stem(n,fen)
xlabel("n")
ylabel("hamming")
grid()
              

figB.pdf

Voici la nouvelle réponse impulsionnelle :

b = b*fen
figure(figsize=(8,4))
stem(n,b)
xlabel("n")
ylabel("b")
grid()
              

figC.pdf

et la réponse fréquentielle :

w,H=scipy.signal.freqz(b)

figure()
subplot(211)
plot(w/(2*numpy.pi),numpy.absolute(H))
ylabel("G")
grid()
subplot(212)
plot(w/(2*numpy.pi),numpy.unwrap(numpy.angle(H)))
xlabel("f/fe")
ylabel("phase")
grid()
                         

figC2.pdf

Avec la fenêtre de Hamming, le gain est bien décroissant.