Filtres RIF : filtrage par convolution et par transformée de Fourier rapide

2.a. Signal discret et transformée de Fourier à temps discret

Soit x(t) un signal continu (fonction à valeurs réelles) échantillonné avec une période T_e. Le signal discret (ou signal échantillonné) correspondant est la suite

$$x_n=x\left(n{T}_e\right)\text{(1)}$$

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

$$\tilde{X}\left(f\right)={\int }_{-\infty }^{+\infty }x\left(t\right)exp(-i2\pi ft)dt\text{(2)}$$

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

$$\tilde{X}\left(f\right)\simeq \sum _{n=-\infty }^{+\infty }{x}_{n}exp(-i2\pi fn{T}_e)\text{(3)}$$

Posons :

$$X\left(f\right)=\sum _{n=-\infty }^{+\infty }{x}_{n}exp(-i2\pi fn{T}_e)\text{(4)}$$

$X\left(f\right)$ est la transformée de Fourier à temps discret (TFTD) du signal discret. Cette fonction à valeurs complexes définit le spectre du signal échantillonné.

2.b. Transformée en Z

La transformée en Z bilatérale du signal discret est une fonction de la variable complexe Z définie par :

$$X\left(Z\right)=\sum _{n=-\infty }^{\infty }{x}_n{Z}^{-n}\text{(5)}$$

Si l'on pose :

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

on a :

$$X\left(f\right)=X\left(Z\left(f\right)\right)\text{(7)}$$

Autrement dit, la transformée de Fourier à temps discret (c'est-à-dire le spectre du signal) s'obtient à partir de la transformée en Z en remplaçant Z par $exp(-i2\pi f{T}_e)$.

Les signaux réels sont nuls avant un certain instant (l'instant t=0 par convention) et on définit une transformée en Z unilatérale :

$$X\left(Z\right)=\sum _{n=0}^{\infty }{x}_n{Z}^{-n}\text{(8)}$$

2.c. Transformée de Fourier discrète

La somme $\mathrm{(4)}$ définissant la TFTD n'est pas calculable numériquement puisqu'elle comporte à nombre infini de termes. Tout au plus peut-on calculer cette somme pour N échantillons :

$$X\left(f\right)\approx \sum _{n=0}^{N-1}{x}_{n}exp(-i2\pi fn{T}_e)\text{(9)}$$

Dans cette écriture, on fait commencer la somme à n=0 par simplification, mais le bloc de N échantillons considéré peut en réalité commencer à un n quelconque. On s'intéresse à cette somme pour les fréquences suivantes :

$$f_k=k\frac{1}{N{T}_e}k=0,1,\dots N-1\text{(10)}$$

On posera aussi $T=N{T}_e$ et $f_k=\frac{k}{T}$. Cela nous amène à la définition de la transformée de Fourier discrète (TFD) :

$$X_k=\sum _{n=0}^{N-1}{x}_{n}exp(-i2\pi \frac{kn}{N})\text{(11)}$$

La TFD est donc une approximation de la TFTD calculée avec un nombre fini d'échantillons et pour des fréquences multiples de $\frac{1}{T}$.

La TFD constitue le spectre d'un signal échantillonné dont le nombre d'échantillons est fini, ce qui est toujours le cas pour les signaux réels.

La TFD peut être utilisée pour calculer les coefficients de Fourier d'une fonction périodique de période T dont on dispose de N échantillons sur cette durée. Cette application est exposée dans Transformée de Fourier discrète: série de Fourier. Dans cet usage, il peut arriver qu'on place un facteur $\frac{1}{N}$ dans la définition de la TFD. Dans le présent document, nous ne mettons pas ce facteur. La TFD ainsi définie correspond à ce que calcule la fonction numpy.fft.fft.

Les N échantillons peuvent être calculés à partir des N valeurs de la TFD par la transformée de Fourier discrète inverse [1] :

$$x_n=\frac{1}{N}\sum _{k=0}^{N-1}{X}_{k}exp\left(i2\pi \frac{kn}{N}\right)\text{(12)}$$

Par comparaison avec la TFD, il a un changement de signe de l'argument de l'exponentiel et le facteur $\frac{1}{N}$.

Pour un signal discret comportant N échantillons (qui est généralement une partie d'un signal plus grand), les N valeurs de la TFD $X_k$ constituent la représentation fréquentielle du signal (autrement dit son spectre) alors que les N valeur $x_n$ constituent sa représentation temporelle. La TFD permet de passer de la représentation temporelle à la représentation fréquentielle et la TFD inverse permet de faire le passage inverse.

L'utilisation de la TFD pour faire l'analyse spectrale d'un signal numérique est développée dans Analyse spectrale d'un signal numérique. Dans le présent document, nous allons utiliser la TFD pour faire une filtrage dans le domaine fréquentiel.

Le calcul de la TFD (et de la TFD inverse) de N nombres nécessite en principe le calcul de N termes, chacun étant une somme de N termes (avec un cosinus et un sinus à évaluer). Les propriétés de la TFD font qu'il existe un algorithme de calcul de complexité temporelle O(N ln(N)) pour la calculer, l'algorithme de transformée de Fourier rapide.

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 _{m=0}^{M-1}{h}_m{x}_{n-m}\text{(13)}$$

où les M coefficients h_m 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{(14)}\\ \hfill & x_n=0sin\ne 0\hfill & \hfill \text{(15)}\end{array}$$

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

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

C'est pourquoi la suite h_n est la réponse impulsionnelle du filtre. La réponse impulsionnelle est ici finie puisqu'elle comporte M valeurs (sous-entendu, elle est nulle à partir du rang M).

La relation $\mathrm{(13)}$ définit un fitrage par convolution. y est le produit de convolution (acyclique) entre x et h. Le nombre M est la longueur du signal h, c'est-à-dire la longueur de la réponse impulsionnelle.

3.b. Fonction de transfert en Z

Considérons alors la transformée en Z du signal en sortie du filtre :

$$\begin{array}{ccc}\hfill Y\left(Z\right)& =\sum _{n=-\infty }^{n=+\infty }\sum _{m=0}^{M-1}{h}_m{x}_{n-m}{Z}^{-n}\hfill & \hfill \text{(17)}\\ \hfill & =\sum _{n\text{'}=-\infty }^{n\text{'}=+\infty }\sum _{m=0}^{M-1}{h}_m{x}_{n\text{'}}{Z}^{-n\text{'}}{Z}^{-m}\hfill & \hfill \text{(18)}\\ \hfill & =\left(\sum _{m=0}^{M-1}{h}_m{Z}^{-m}\right)\sum _{n\text{'}=-\infty }^{n\text{'}=+\infty }{x}_{n\text{'}}{Z}^{-n\text{'}}\hfill & \hfill \text{(19)}\\ \hfill & =H\left(Z\right)X\left(Z\right)\hfill & \hfill \text{(20)}\end{array}$$

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

$$H\left(Z\right)=\sum _{m=0}^{M-1}{h}_m{Z}^{-m}\text{(21)}$$

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{(22)}$$

Le signal en sortie s'écrit :

$$\begin{array}{ccc}\hfill y_n=& \sum _{m=0}^{M-1}{h}_m{x}_{n-m}\hfill & \hfill \text{(23)}\\ \hfill =& exp\left(i2\pi nf{T}_e\right)\sum _{m=0}^{M-1}{h}_{m}exp(-i2\pi fm{T}_e)\hfill & \hfill \text{(24)}\\ \hfill =& exp\left(i2\pi nf{T}_e\right)H\left(Z\right)\hfill & \hfill \text{(25)}\\ \hfill =& x_{n}H\left(Z\right)\hfill & \hfill \text{(26)}\end{array}$$

avec

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

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

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

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_k=\frac{n}{M{T}_e}k=0,1,\dots M-1\text{(29)}$$

On a alors

$$H\left(f_k\right)=\sum _{m=0}^{M-1}{h}_{m}exp(-i2\pi m\frac{k}{M})\text{(30)}$$

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_m=\frac{1}{M}\sum _{k=0}^{M-1}H\left(f_k\right)exp\left(i2\pi k\frac{m}{M}\right)\text{(31)}$$

La conception d'un filtre RIF consiste à calculer M coefficients h_m pour obtenir une réponse fréquentielle du type voulu (par exemple un filtre passe-bas). La méthode la plus utilisée (et la plus simple) consiste à partir d'une réponse fréquentielle idéale, à calculer la réponse impulsionnelle correspondante (qui est infinie) et à la multiplier par une fonction de fenêtrage pour obtenir une réponse impulsionnelle finie. Voir Filtres à réponse impulsionnelle finie pour l'exposé de cette méthode et Conception et mise en œuvre des filtres numériques pour la conception de ce type de filtre avec la fonction scipy.signal.firwin.

3.d. Fitrage d'un signal en temps réel

Dans le cas d'un signal x_n délivré par un convertisseur analogique-numérique, le filtrage en temps réel consiste à calculer les valeurs de y_n en temps réel, c'est-à-dire que le calcul de y_n doit être terminé au plus tard lorsque la valeur x_n+1 est délivré par le convertisseur. La relation $\mathrm{(13)}$ permet de faire ce calcul sans difficulté, à condition que M ne soit pas trop grand pour que le calcul de la somme puisse se faire pendant la durée T_e (la période d'échantillonnage). Évidemment, le nombre M maximal dépend de la vitesse de calcul du processeur.

Il faut remarquer que le nombre d'échantillons total d'un signal permanent n'est pas connu a priori, ce qui ne pose pas de problème justement parce qu'on le filtre en temps réel.

4.a. Définition

Par opposition au filtrage en temps réel, le filtrage d'un signal de longueur déterminé consiste à filtrer un signal comportant un nombre N d'échantillons (supposé supérieur à M). Il peut s'agir de filtrer un signal mémorisé lors d'un enregistrement ou bien, comme nous le verrons plus loin, de filtrer un bloc d'un signal permament.

Supposons par convention que le signal à filtrer commence à l'indice n=0 et finisse à l'indice n=N-1. La relation $\mathrm{(13)}$ pose problème car x_n-m n'est pas défini lorsque n-m<0. Il existe trois manières de résoudre ce problème :

• 1. Commencer le calcul de y_n à partir du rang M-1, les M-1 premiers termes étant nuls.
• 2. Commencer le calcul de y_n mais n'utiliser qu'une partie de h_m pour les m premières valeurs.
• 3. Faire le calcul en utilisant une condition périodique : n-m signifie dans ce cas n-m modulo N.

La troisième possibilité définit le produit de convolution cyclique, dont l'intérêt apparaîtra plus loin.

La fonction scipy.signal.convolve permet de calculer la convolution de deux signaux de taille finie.

from scipy.signal import convolve
import numpy as np
from matplotlib.pyplot import *

h = [1,1,1]
x = np.arange(10)
y = convolve(x,h,mode='valid')
                

print(x)
--> array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])

print(y)
--> array([ 3,  6,  9, 12, 15, 18, 21, 24])

L'option mode='valid' correspond au cas 1 mais les premiers termes nuls ne sont pas renvoyés, donc le tableau final est plus court de M éléments que le tableau x.

y = convolve(x,h,mode='same')
                

print(y)
--> array([ 1,  3,  6,  9, 12, 15, 18, 21, 24, 17])

L'option mode='same' correspond au cas 2 mais la convolution est centrée. Lorsque M est impair (M=2P+1), la convolution centrée se définit par :

$$y_n=\sum _{p=-P}^P{h}_{p+P}{x}_{n+p}\text{(32)}$$

y = convolve(x,h,mode='full')
                

print(y)
--> array([ 0,  1,  3,  6,  9, 12, 15, 18, 21, 24, 17,  9])

4.b. Produit de convolution cyclique

Considérons deux signaux discrets comportant N échantillons. Le produit de convolution cyclique de ces deux signaux [1] est défini par

$${\left(x\otimes y\right)}_n=\sum _{m=0}^{M-1}{x}_m{y}_{n-m}\text{(33)}$$

où n-m désigne n-m modulo M. Avec le changement de variable $\ell =n-m$, on peut écrire ce produit de convolution :

$${\left(x\otimes y\right)}_n=\sum _{\ell =n}^{n-N-1}{x}_{n-\ell }{y}_{\ell }\text{(34)}$$

Compte tenu de la condition limite cyclique, la somme sur N termes peut aussi s'écrire :

$${\left(x\otimes y\right)}_n=\sum _{\ell =0}^{N-1}{x}_{n-\ell }{y}_{\ell }\text{(35)}$$

ce qui démontre que le produit de convolution cyclique est commutatif :

$$x\otimes y=y\otimes x\text{(36)}$$

4.c. Théorème de convolution

Considérons la transformée de Fourier discrète de la convolution cyclique à N éléments :

$$\begin{array}{cc}\hfill {TFD}_k\left(x\otimes y\right)& =\sum _{n=0}^{N-1}{\left(x\otimes y\right)}_{n}exp(-j2\pi \frac{nk}{N})\hfill \\ \hfill & =\sum _{n=0}^{N-1}\sum _{m=0}^{M-1}{x}_m{y}_{n-m}exp(-j2\pi \frac{nk}{N})\hfill \\ \hfill & =\sum _{m=0}^{N-1}{x}_m\sum _{n=0}^{N-1}{y}_{n-m}exp(-j2\pi \frac{nk}{N})\hfill \\ \hfill \end{array}$$

La TFD de y s'écrit :

$$Y_k=\sum _{n=0}^{N-1}{y}_{n}exp(-i2\pi \frac{kn}{N})\text{(37)}$$

d'où l'on déduit :

$$Y_{k}exp(-i2\pi \frac{mk}{N})=\sum _{n=0}^{N-1}{y}_{n}exp(-i2\pi \frac{k(n+m)}{N})=\sum _{n\text{'}=m}^{N-1+m}{y}_{n\text{'}-m}exp(-i2\pi \frac{kn\text{'}}{N})\text{(38)}$$

Si on ajoute N à n' dans l'exponentielle, celle-ci est inchangée. Il en est de même pour $y_{n\text{'}-m}$ car $n\text{'}-m$ est en fait égal à $n\text{'}-m$ module N. Il s'en suit que :

$$Y_{k}exp(-i2\pi \frac{mk}{N})=\sum _{n=0}^{N-1}{y}_{n-m}exp(-i2\pi \frac{kn}{N})\text{(39)}$$

On a finalement :

$$\begin{array}{cc}\hfill {TFD}_k\left(x\otimes y\right)& =\sum _{m=0}^{N-1}{x}_m{Y}_{k}exp(-i2\pi \frac{mk}{N})\hfill \\ \hfill & =\left(\sum _{m=0}^{N-1}{x}_{m}exp(-i2\pi \frac{mk}{N})\right)Y_k\hfill \\ \hfill & =X_k{Y}_k\hfill \end{array}$$

Ce résultat constitue le théorème de convolution : la TFD du produit de convolution cyclique de deux signaux de longueur N est égale au produit des TFD des deux signaux :

$${TFD}_k\left(x\otimes y\right)=X_k{Y}_k\text{(40)}$$

Formulé autrement, un produit de convolution dans le domaine temporel est équivalent à une multiplication dans le domaine fréquentiel.

Ce résultat suggère une réalisation possible du filtrage RIF dans le domaine fréquentiel, en multipliant la TFD du signal x par la TFD de la réponse impulsionnelle h, qui n'est rien d'autre que la réponse fréquentielle du filtre H. L'avantage de cette méthode de calcul réside dans la complexité $O\left(Nln\left(N\right)\right)$ de l'algorithme de transformée de Fourier rapide qui permet d'obtenir la TFD de x. Le calcul de cette TFD est suivie de N multiplications puis d'une TFD inverse. Le filtrage complet a donc une complexité $O\left(Nln\left(N\right)\right)$ , à comparer à la complexité O(N²) du filtrage par convolution.

Il y a cependant une difficulté : la convolution que l'on veut réaliser est celle donnée par la relation $\mathrm{(13)}$ et non pas la convolution cyclique.

4.d. Ajout de zéros

L'ajout de zéros aux signaux h et x permet de réaliser la convolution $\mathrm{(13)}$ au moyen d'une convolution cyclique. Initialement, il y a N valeurs x_n et M valeurs h_m. N est la longueur du signal à filtrer, M est la longueur de la réponse impulsionnelle. Tout d'abord, il faut compléter h par des zéros afin que les deux signaux aient la même longueur. De plus, il faut ajouter M-1 zéros au début de x afin que la convolution cylique appliquée sur les M-1 premiers éléments de x ne fasse pas intervenir la condition limite périodique. On obtient ainsi un filtrage par convolution qui correspond au cas 2 mentionné plus haut : pour les M premières échantillons de x (avant ajout de zéro), les échantillons manquants sont remplacés par des zéros.

Finalement, il faut ajouter des zéros aux deux signaux h et x afin que leur taille N' soit N+M-1 au minimum. On prend le nombre N' égal à une puissance de 2 le plus petit supérieur à N+M-1.

L'ajout de zéro dans le domaine temporel est une technique utilisée dans l'analyse spectrale d'un signal numérique afin d'améliorer la définition fréquentielle. En effet, ajouter des zéros dans le domaine temporel est équivalent à augmenter l'échantillonnage des fréquences dans le domaine fréquentiel, ce qui permet de s'approcher de la TFTD. Dans le cas du filtrage RIF, si on n'ajoute par le nombre minimal de zéro nécessaire, l'échantillonnage des fréquences de la TFD n'est pas assez fin et il se produit un phénomène de repliement temporel, analogue au phénomène de repliement fréquentiel qui se produit lorsque la fréquence d'échantillonnage d'un signal analogique n'est pas assez grande.

4.e. Filtrage dans le domaine fréquentiel

L'algorithme de filtrage dans le domaine fréquentiel, aussi nommé convolution par transformée de Fourier rapide, se déduit de ce qui précède. On dispose d'un signal discret x de taille N et d'un filtre de réponse impulsionnelle h de taille M.

• Déterminer la taille N' des signaux égale à une puissance de 2 et supérieure ou égale à N+M-1.
• Ajouter N'-M à la réponse impulsionnelle h (à la fin).
• Ajouter N'-N zéro au signal x. En principe, il faut ajouter M-1 zéros au début et le reste à la fin.
• Calculer les TFD X_k et H_k ($k=0,1,\dots N\text{'}-1$).
• Obtenir le signal y de longueur N' par transformée de Fourier inverse de X_kH_k.

Voyons un exemple. Supposons que l'on dispose d'un signal x de taille N=1000, par exemple une sinusoïde échantillonnée :

N = 1000
f = 83
fe = 10000 # fréquence d'échantillonnage
T = N/fe # durée du signal échantillonné
t = np.arange(N)/fe
x = np.sin(2*np.pi*f*t)
figure(figsize=(16,6))
plot(t,x)
grid()
xlabel('t')
ylabel('x')
                
                

fig1.pdf

Calculons les coefficients d'un filtre RIF passe-bas :

from scipy.signal import firwin,freqz
P = 100
M = 2*P+1
fc = 500 # fréquence de coupure
h = firwin(numtaps=M,cutoff=[fc/fe],window='hann',nyq=0.5)
                

et traçons sa réponse fréquentielle :

w,H=freqz(h,1000)
figure(figsize=(12,8))
subplot(211)
plot(w/(2*np.pi),20*np.log10(np.absolute(H)))
xlabel("f/fe")
ylabel("GdB")
grid()
subplot(212)      
plot(w/(2*np.pi),np.unwrap(np.angle(H)))
xlabel("f/fe")
ylabel("phase")
grid()
                

fig2.pdf

Le signal défini précédemment est dans la bande passante.

Voici tout d'abord le résultat du filtrage par convolution (simulation d'un fitrage en temps réel) :

y = convolve(x,h,mode='valid')
ty = (np.arange(len(y))+2*P)*1/fe
figure(figsize=(16,6))
plot(t,x)
plot(ty,y)
grid()
xlabel('t')
ylabel('x')
                

fig3.pdf

Pour faire le filtrage dans le domaine fréquentiel, commençons par ajouter des zéros à h afin d'atteindre la taille N :

from numpy.fft import fft,ifft
h1 = np.zeros(N)
h1[0:M] = h
X = fft(x)
H1 = fft(h1)
y1 = np.real(ifft(X*H1))
figure(figsize=(16,6))
plot(t,x)
plot(t,y1)
grid()
xlabel('t')
ylabel('x')
                

fig4.pdf

La transformée de Fourier discrète inverse renvoie en principe un tableau de nombres complexes dont la partie imaginaire est nulle mais des erreurs d'arrondis font qu'il existe tout de même une petite partie imaginaire; on doit donc prendre la partie réelle du tableau.

L'effet du repliement temporel est bien visible au début du signal. Il est donc impératif d'ajouter des zéros au signal x pour éliminer l'effet cyclique de la convolution. Voici la taille à atteindre par ajout de zéros :

Np = int(2**np.ceil(np.log2(N+M-1)))
                 

print(Np)
--> 2048

h2 = np.zeros(Np)
h2[0:M] = h
x2 = np.zeros(Np)
x2[M:M+N] = x
X2 = fft(x2)
H2 = fft(h2)
y2 = np.real(ifft(X2*H2))
figure(figsize=(16,6))
plot(x2,label='x2')
plot(y2,label='y2')
grid()
xlabel('n')
ylabel('x')
legend(loc='upper right')
                 

fig5.pdf

L'effet cyclique de la convolution a disparu mais le signal y présente des ondulations indésirables au début et surtout à la fin. Pour comprendre ce phénomène, voyons le spectre du signal x (après ajout de zéros), qui est le tracé du module de sa TFD en fonction de la fréquence, et la réponse fréquentielle du filtre :

f2 = np.arange(Np)*fe/Np
figure(figsize=(12,6))
subplot(211)
plot(f2,np.absolute(X2))
ylabel(r"$|X_2|$")
grid()
subplot(212)
plot(f2,np.absolute(H2))
xlabel('f')
ylabel(r"$|H_2|$")
grid()
                  

fig6.pdf

Voici un détail à basse fréquence :

figure(figsize=(12,6))
subplot(211)
plot(f2,np.absolute(X2))
ylabel(r"$|X_2|$")
grid()
xlim(0,1000)
subplot(212)
plot(f2,np.absolute(H2))
xlabel('f')
ylabel(r"$|H_2|$")
grid()
xlim(0,1000)

                  

fig7.pdf

La raie du spectre du signal sinusoïdal comporte des lobes secondaires dont l'intensité est grande. L'intensité de ces lobes peut être fortement réduite au moyen d'une fonction de fenêtrage. Voir Fenêtres pour l'analyse spectrale pour la définition de différentes fenêtres. Appliquons par exemple une fenêtre de Hamming (avant de calculer la TFD) :

from scipy.signal.windows import get_window
w = get_window('hamming',N)
a = 0.54 # facteur de normalisation
x2 = np.zeros(Np)
x2[M:M+N] = x*w
X2 = fft(x2)
figure(figsize=(12,6))
plot(f2,np.absolute(X2))
ylabel(r"$|X_2|$")
grid()
xlim(0,1000)
xlabel('f')
                    

fig8.pdf

Voyons le filtrage avec ce signal :

y2 = np.real(ifft(X2*H2))
figure(figsize=(16,6))
plot(x2,label='x2')
plot(y2,label='y2')
grid()
xlabel('n')
ylabel('x')
legend(loc='upper right')
                 

fig9.pdf

Les ondulations indésirables à la fin et au début ont disparu. Évidemment, l'inconvénient est la modulation de l'amplitude du signal due au fenêtrage. Nous allons voir dans la partie suivante comme se problème est résolu dans le cadre d'un filtrage avec une fenêtre mobile.

5.a. Fenêtre mobile

La technique de filtrage RIF par transformée de Fourier rapide a été développée dans la partie précédente mais dans le cas du filtrage d'un signal de taille N déterminée.

Il s'agit d'appliquer cette technique au moyen d'une fenêtre mobile pour filtrer un signal en temps réel. Comparé au filtrage dans le domaine temporel (filtrage par convolution), le filtrage dans le domaine fréquentiel (c'est-à-dire la multiplication des TFD) devrait fournir un gain notable de vitesse lorsque la longueur de la réponse impulsionnelle est grande (au moins une centaine). Cette méthode porte parfois le nom de convolution par transformée de Fourier, ce qui est une appelation incorrecte puisqu'il ne s'agit pas à proprement parler d'une convolution. Une dénomination plus convenable serait filtrage RIF par transformée de Fourier ou bien fitrage RIF dans le domaine fréquentiel.

Le signal est traité par bloc de N échantillons. L'extraction d'un bloc de taille N consiste à appliquer au signal une fenêtre de largeur N. Comme nous l'avons vu précédemment, cette fenêtre devra être pourvue d'une fonction de fenêtrage (ou fonction de pondération) avec une réduction de l'amplitude sur les bords, de manière à réduire les lobes secondaires dans le spectre du signal x. Il est courant de désigner par le mot fenêtre la combinaison de la fenêtre elle-même et de la fonction de fenêtrage. Deux positions successives de la fenêtre sont espacées d'un nombre Q d'échantillons et ces deux positions doivent se recouvrir, comme le montre le schéma suivant :

Figure pleine page

Ce fenêtrage doit avoir la propriété suivante : la somme des portions obtenues par le fenêtrage doit reconstituer le signal initial (à une constante multiplicative près). Si $w\left(n\right)$ désigne la fenêtre, en tant que fonction de l'indice étendue à l'ensemble des entiers relatifs, cette propriété s'écrit :

$$\sum _{m=-\infty }^{\infty }w(n-mQ)=1\forall n\text{(41)}$$

Une fenêtre vérifiant cette propriété pour un déplacement Q est une fenêtre COLA pour ce déplacement (constant-overlap-add) [2]. Par exemple, la fenêtre rectangulaire de longueur N est COLA pour Q=N (aussi pour Q=N/2). Cependant, nous avons vu plus que la fenêtre rectangulaire est à proscire en raison des lobes secondaires du spectre trop intenses.

D'après [2], la fenêtre de Hamming pour N impair est COLA pour Q=(N-1)/2. Vérifions numériquement cette propriété en appliquant la fenêtre mobile à un signal dont la valeur est constante :

Nt = 1000
x = np.ones(Nt)
N = 101
Q = int((N-1)/2)
w = get_window('hamming',N,fftbins=False)
x2 = np.zeros(Nt)
figure(figsize=(16,6))
for m in range(int(Nt/Q)-2):
    xf = np.zeros(Nt)
    xf[m*Q:m*Q+N] = x[m*Q:m*Q+N]*w
    plot(xf,'k--')
    x2 += xf
plot(x2,'r.')
xlabel('n')
grid()
                             

fig10.pdf

L'option fftbins = False permet de générer une fenêtre symétrique, alors que fftbins = True (valeur par défaut) génère une fenêtre dite périodique. Voici une comparaison de ces deux fenêtres pour une longueur impaire :

figure(figsize=(12,6))
plot(get_window('hamming',11,fftbins=True),'o-r',label='periodic')
plot(get_window('hamming',11,fftbins=False),'o-b',label='symmetric')
legend(loc='upper right')
                              

fig11.pdf

La première est dite périodique car elle s'applique aux échantillons d'un signal périodique échantillonné sur exactement une période (la valeur d'indice virtuel N doit être égale à celle d'indice 0). La fenêtre périodique ne vérifie pas du tout la propriété COLA. La figure ci-dessus montre que la fenêtre symétrique pour N impaire ne vérifie pas tout à fait la propriété COLA. Il faut la modifier de la manière suivante [2] :

w[0] /= 2
w[N-1] /=2
x2 = np.zeros(Nt)
figure(figsize=(16,6))
for m in range(int(Nt/Q)-2):
    xf = np.zeros(Nt)
    xf[m*Q:m*Q+N] = x[m*Q:m*Q+N]*w
    plot(xf,'k--')
    x2 += xf
plot(x2,'r-')
xlabel('n')
grid()
                               

fig12.pdf

Diviser le premier et le dernier élément de la fenêtre de Hamming par 2 (cas symétrique avec N impair) conduit bien à une fenêtre COLA. Cependant, la somme est supérieure à 1. Voici sa valeur :

print(x2[200:800].mean())
--> 1.0800000000000003

Dans la suite de ce document, nous utiliserons cette fenêtre (divisée par la valeur ci-dessus).

5.b. Algorithme de filtrage

Le filtrage en temps réel dans le domaine fréquentiel, ou filtrage RIF par transformée de Fourier, se fait de la manière suivante.

Configuration pour une fenêtre de Hamming :

• Définir un filtre par sa réponse impulsionnelle h de longueur M.
• Définir une taille de fenêtre N impaire.
• Définir le déplacement de la fenêtre Q=(N-1)/2.
• Calculer une fenêtre de hamming COLA pour ce déplacement, fenêtre w de taille N.
• Calculer N', la plus petite puissance de 2 supérieure ou égale à N+M-1.
• Ajouter N'-M zéro à la fin de h afin d'obtenir une réponse impulsionnelle de longueur N'.
• Calculer la TFD de h ainsi complétée par des zéros, notée H.

Dans la boucle de filtrage

• Sélectionner une portion du signal, notée x, de longueur N et lui appliquer la fonction w.
• Ajouter M zéros au début de x et N'-M-N à la fin, afin d'obtenir un signal de longueur N'.
• Calculer la TFD de x, notée X.
• Calculer $y=TF{D}^{-1}\left(HX\right)$.
• Recommencer après application d'un déplacement de la fenêtre de longueur Q.
• Additionner les blocs $y$ de taille N' afin d'obtenir le signal de sortie.

La figure suivante montre les blocs de taille N' qui sont obtenus par le filtrage et qui doivent être additionnés pour donner le signal de sortie du filtre.

Figure pleine page

Voici une implémentation de ce filtrage dans une classe.

class FiltreFFT:
    def __init__(self,h,N):
        self.M = len(h)
        if N%2==0: N+=1
        self.N = N
        self.Q = int((N-1)/2)
        self.w = get_window('hamming',N,fftbins=False)
        self.w[0] /= 2
        self.w[N-1] /= 2
        self.w /= 1.08
        self.Np = int(2**np.ceil(np.log2(self.N+self.M-1)))
        h2 = np.zeros(self.Np)
        h2[0:self.M] = h
        self.H = fft(h2)
        
        
    def filtrageFenetre(self,signal,indice):
        x2 = np.zeros(self.Np)
        x2[self.M:self.M+self.N] = signal[indice:indice+self.N]*self.w
        X = fft(x2)
        return np.real(ifft(X*self.H))
        
    def filtrageSignal(self,signal):
        Nt = len(signal)
        sortie = np.zeros(Nt)
        m = 0
        while m*self.Q+self.Np < Nt:
            y = self.filtrageFenetre(signal,m*self.Q)
            #print(len(y[np.where(y!=0)]))
            y2 = np.zeros(Nt)
            y2[m*self.Q:m*self.Q+self.Np] = y
            sortie += y2
            plot(y2,'k--')
            m += 1
        return sortie
            
        
                    

La fonction filtrageSignal effectue le filtrage d'un signal stocké dans un tableau mais elle simule le fonctionnement d'un filtre en temps réel. L'implémentation d'un véritable filtre temps réel sera cependant différente (et plus délicate) car on ne pourra pas dans ce cas créér un tableau y2.

Reprenons l'exemple précédent :

Nt = 10000
f = 83
fe = 10000 # fréquence d'échantillonnage
T = Nt/fe # durée du signal échantillonné
t = np.arange(Nt)/fe
signal = np.sin(2*np.pi*f*t)
N = 1025
P = 100
M = 2*P+1
fc = 500 # fréquence de coupure
h = firwin(numtaps=M,cutoff=[fc/fe],window='hann',nyq=0.5) 
filtre = FiltreFFT(h,N)
figure(figsize=(16,6))
sortie = filtre.filtrageSignal(signal)

plot(signal)
plot(sortie)
xlabel('n',fontsize=16)
grid()            
                    

fig13.pdf

La courbe en trait interrompu noir est le résultat du filtrage de chaque bloc. Voici un détail :

xlim(2000,4000)

                    

fig14.pdf

Pour comparaison, écrivons une fonction qui effectue le filtrage par convolution en simulant le fonctionnement d'un filtre temps réel :

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

Voici le résultat de ce filtrage par convolution :

sortie_conv = filtrage_convolution(signal,h)
figure(figsize=(16,6))
plot(signal)
plot(sortie_conv)
xlim(2000,4000)
xlabel('n',fontsize=16)
grid()  
                    

fig15.pdf

On remarque que le retard entre les deux signaux n'est pas le même. Pour le filtrage par convolution, ce retard est PT_e (où M=2P+1). Dans le cas du filtrage dans le domaine fréquentiel, le retard est beaucoup plus grand (attention, on ne voit sur ces courbes que le retard modulo la période).

Voici le filtrage d'un signal dont la fréquence est égale à la fréquence de coupure du filtre :

f = fc
signal = np.sin(2*np.pi*f*t)
figure(figsize=(16,6))
sortie = filtre.filtrageSignal(signal)
plot(signal)
plot(sortie)
xlabel('n',fontsize=16)
grid() 
                     

fig16.pdf

Comme prévu, le gain vaut 1/2. Voici un détail :

xlim(2000,2500)
                     

fig17.pdf

5.c. Filtrage en temps réel

Le filtrage en temps réel consiste à filtrer en continu un signal permanent au fur et à mesure que le convertisseur analogique-numérique fournit les échantillons. Temps réel signifie qu'il y a une contrainte de temps : le temps de calcul doit être assez faible pour que le retard entre le signal de sortie et le signal d'entrée reste constant. Cependant, ce retard peut être grand, beaucoup plus grand que le retard d'un filtrage par convolution.

L'implémentation du filtrage en temps réel nécessite un tampon servant à calculer la somme des blocs filtrés. Ce tampon peut être qualifié d'accumulateur. Lorsqu'un bloc de taille N' est obtenu (par la fonction filtrageFenetre), sont contenu doit être additionné dans l'accumulateur.

La longueur des blocs à sommer est N' et le déplacement entre deux blocs successifs est Q. La figure suivante montre ces blocs :

Figure pleine page

À chaque fois qu'un bloc de taille N' est filtré, il y a une nouvelle portion du signal de sortie, de longueur Q, dont le traitement est terminé. Or la durée de numérisation d'un bloc de longueur Q est QT_e. Il faut donc que le calcul complet du filtrage de ce bloc prenne un temps inférieur à QT_e. Cela signifie aussi que l'accumulateur est en réalité un tampon circulaire de taille N' dans lequel, à chaque itération, les Q premiers échantillons (ceux dont la somme est terminée) sont remplacés par les Q derniers (qui sont initialement nuls). Il s'agit donc d'une pile FIFO (First In = First Out). Avant d'être retiré, les Q premiers échantillons doivent être soit stockés en mémoire, soit envoyés vers un convertisseur numérique-analogique. Pour implémenter le tampon circulaire, on utilise un tableau de taille N' auquel on associe un indice first qui indique la position dans le tableau du premier élément. Les Q premiers éléments de l'accumulateur sont les Q éléments qui suive l'élément d'indice first. La figure suivante montre deux configurations successives de l'accumulateur :

Figure pleine page

À chaque itération, le tableau de taille N' obtenu par la TFD inverse est ajouté à l'accumulateur, le premier élément du tableau étant ajouté à l'élément d'indice first (second état représenté ci-dessus). Les entrées et sorties de blocs de taille Q sont simplifiées si Q est un sous multiple de N'. Il faut que Q soit une puissance de 2, auquel cas N-1 est également une puissance de 2. C'est bien le cas dans l'exemple précédent puisqu'on avait N=1025, N'=2048 et Q=512. Cela conduit à la règle suivantes : choisir N=2^a+1 et Q=2^a-1. Si la longueur de la réponse impulsionnelle est inférieure à N, cela conduit à N'=2^a+1.

Les échantillons qui proviennent du convertisseur analogique-numérique doivent être stockés dans un tampon circulaire de longueur N, qui contient en permanence les N derniers échantillons. On associe à ce tampon un indice debut qui marque le début du tampon (premier des N derniers échantillons). Tous les Q échantillons, on copie les N valeurs de ce tampon vers un tableau x de taille N, que l'on complète par des zéros pour obtenir un tableau de taille N', puis la TFD de ce tableau est calculée et enfin multipliée par celle de h. La TFD inverse fournit les N' valeurs qui sont additionnées dans l'accumulateur. La position de début de lecture des N valeurs dans le tampon est incrémenté de Q à chaque itération.

Nous réécrivons la classe précédente.

class FiltreFFT2:
    def __init__(self,h,a):
        self.M = len(h)
        N = 2**a+1
        self.N = N
        self.Q = 2**(a-1)
        self.w = get_window('hamming',N,fftbins=False)
        self.w[0] /= 2
        self.w[N-1] /= 2
        self.w /= 1.08
        self.Np = int(2**np.ceil(np.log2(self.N+self.M-1)))
        h2 = np.zeros(self.Np)
        h2[0:self.M] = h
        self.H = fft(h2)
        self.accum = np.zeros(self.Np)
        self.first = 0
        self.tamponSignal = np.zeros(self.N)
        self.debut = -1
        self.compt = 0
        self.sortie = np.zeros((0))
        
    def filtrageFenetre(self,x):
        x2 = np.zeros(self.Np)
        x2[self.M:self.M+self.N] = x*self.w
        X = fft(x2)
        return np.real(ifft(X*self.H))
        
    def entreeEchantillon(self,u):
        # u : échantillon
        self.debut = (self.debut+1)%self.N
        self.tamponSignal[self.debut] = u
        self.compt += 1
        if self.compt == self.Q:
            self.compt = 0
            x = np.zeros(self.N)
            i1 = self.N-self.debut
            x[0:i1] = self.tamponSignal[self.debut:self.N]
            x[i1:self.N] = self.tamponSignal[0:self.debut]
            y = self.filtrageFenetre(x) # tableau de longueur N'
            self.sortie = np.concatenate((self.sortie,self.accum[self.first:self.first+self.Q])) # sortie des Q premiers éléments de l'accumulateur 
            self.accum[self.first:self.first+self.Q] = np.zeros(self.Q) # et mise à zéro pour les Q derniers 
            self.first = (self.first+self.Q)%self.Np
            i = self.first
            n = 0
            while n<self.Np: # addition de y dans l'accumulateur (à optimiser)
                self.accum[i] += y[n]
                i = (i+1)%self.Np
                n += 1
            
    def sortieSignal(self):
        return self.sortie    
             
                    

La fonction entreeEchantillon est appelée pour entrer les échantillons un par un. Ces échantillons sont stockés dans tamponSignal. Tous les Q échantillons, on déclenche le filtrage des N derniers échantillons, c'est-à-dire du contenu de tamponSignal. Remarquons que cette implémentation n'est pas optimisée en ce qui concerne l'utilisation de la mémoire. On pourrait se passer du tableau x et recopier le contenu du tampon directement dans un tableau de taille N'. Pour plus de clarté, nous avons préférer garder la fonction filtrageFenetre, qui prend en argument un tableau de taille N.

On reprend l'exemple précédent :

Nt = 10000
f = 83
fe = 10000 # fréquence d'échantillonnage
T = Nt/fe # durée du signal échantillonné
t = np.arange(Nt)/fe
signal = np.sin(2*np.pi*f*t)
a = 10
P = 100
M = 2*P+1
fc = 500 # fréquence de coupure
h = firwin(numtaps=M,cutoff=[fc/fe],window='hann',nyq=0.5) 
filtre = FiltreFFT2(h,a)
for i in range(Nt):
    filtre.entreeEchantillon(signal[i])
sortie = filtre.sortieSignal()
figure(figsize=(16,6))
plot(signal)
plot(sortie)
xlabel('n',fontsize=16)
grid()
                    

fig18.pdf

xlim(2000,4000)
                    

fig19.pdf

Voici un autre exemple : le filtrage d'un signal périodique avec des harmoniques :

Nt = 10000
f = 300
fe = 10000 # fréquence d'échantillonnage
T = Nt/fe # durée du signal échantillonné
t = np.arange(Nt)/fe
signal = np.sin(2*np.pi*f*t) + 0.5*np.cos(2*2*np.pi*f*t) + 0.3*np.sin(3*2*np.pi*f*t+np.pi/3)
a = 10
P = 100
M = 2*P+1
fc = 500 # fréquence de coupure
h = firwin(numtaps=M,cutoff=[fc/fe],window='hann',nyq=0.5) 
filtre = FiltreFFT2(h,a)
for i in range(Nt):
    filtre.entreeEchantillon(signal[i])
sortie = filtre.sortieSignal()
figure(figsize=(16,6))
plot(signal)
plot(sortie)
xlabel('n',fontsize=16)
grid()
xlim(2000,2500)                  
                    

fig20.pdf

Cet exemple met en évidence la très grande sélectivité de ce filtre, dont la réponse impulsionnelle a 201 éléments. Pour une taille aussi grande, il est très probable qu'une implémentation du filtrage par transformée de Fourier rapide (FFT) sur un microcontroleur soit beaucoup plus rapide que le filtrage par convolution. Selon [2], le filtrage par FFT est plus rapide que la convolution à partir d'une taille de 128 environ et la convolution est plus rapide pour les réponses impulsionnelles très courtes (moins de 16). Cependant, seuls des tests de vitesse sur um matériel donné et une implémentation donnée permettra de conclure précisément sur la méthode optimale en fonction de la taille de la réponse impulsionnelle.

En plus de permettre un filtrage plus efficace avec un grand nombre de coefficients, le filtrage RIF par FFT a un autre avantage : le spectre du signal est calculé en temps réel (sur une fenêtre de taille N se déplaçant par saut de longueur Q). Il est donc possible de réaliser un filtre dont la réponse fréquentielle s'adapte au propriétés du signal. On pourrait par exemple concevoir un filtre qui détecte les signaux quasi périodiques et effectue une atténuation (ou un renforcement) de leurs harmoniques.