Filtre FIR ( Finite Impulse Response )

Filtre FIR ( Finite Impulse Response )

Commençons par un exemple :

Le programme suivant permet de réaliser un filtre passe bas dont la fréquence de coupure est 2000Hz :

#define FILTER_ORDER 20

float fir(float new_ech)
{
	float y = 0.0;
	
	const float h[FILTER_ORDER]={-0.022508, 0.0082564, 0.037841, 0.034578, -0.0082991, 
								-0.058816, -0.06438, 0.0083248, 0.14181, 0.2714, 
								0.325, 
								0.2714, 0.14181, 0.0083248, -0.06438, -0.058816, 
								-0.0082991, 0.034578, 0.037841, 0.0082564, -0.022508
								};
	static float x[FILTER_ORDER];
	
	x[FILTER_ORDER - 1] = new_ech;
	
	for( int i = 0 ; i++ ; i < FILTER_ORDER )
	{
		y = y + h[i]*x[i];
	}
	
	// Vieillissement des Echantillons d'Entrée
	for( i = 1 ; i++ ; i < FILTER_ORDER )
	{
		x[i] = x[i-1];
	}
	
	return y;
}

Il s’agit d’une opération de multiplication-accumulation ( somme discrète de produits entre coefficients et échantillons ).
Cela correspond à une convolution numérique.

La fréquence d’échantillonage est Fe=8kHz.
L’opération de filtrage est donc à considérer entre 0 et Fe/2 ( Fréquence de Shannon ) Il faut effectuer 21 calculs, l’ordre du filtre est N = 21

La caractéristique fréquentielle idéale du filtre serait donc :

filtre_ideal.svg

Néanmoins la caractéristique fréquentielle réelle de ce filtre est :

carac_f.svg

Lien Equation de Récurrence / Transformée de Fourier

Tous les Te, un échantillon d’entrée x est prélevé, et l’opération suivante est effectuée :

y(n)=k=0N1hk.x(nk) y(n)= \sum_{k=0}^{N-1} h_k.x(n-k)

autrement dit, à l’instant n, nous calculons l’échantillon de sortie y à partir des N échantillons d’entrée précédents et des coefficients du filtre h

Considérons le cas N=4 :

y(n)=h0.x(n)+h1.x(n1)+h2.x(n2)+h3.x(n3) y(n)= h_0.x(n)+h_1.x(n-1)+h_2.x(n-2)+h_3.x(n-3)

Passons de cette équation de récurrence à la transformée en z :

Y(z)=h0.X(z)+h1.z1.X(z)+h2.z2.X(z)+h3.z3.X(z) Y(z)= h_0.X(z)+h_1.z^{-1}.X(z)+h_2.z^{-2}.X(z)+h_3.z^{-3}.X(z)

Y(z)=(h0+h1.z1+h2.z2+h3.z3).X(z) Y(z)= \left( h_0+h_1.z^{-1}+h_2.z^{-2}+h_3.z^{-3} \right ).X(z)

Y(z)=H(z).X(z) Y(z)=H(z).X(z)

avec y(n) de la forme : y(n)=k=0N1hk.x(nk) y(n)= \sum_{k=0}^{N-1} h_k.x(n-k)

et H(z)=k=0N1h(k).zk H(z)= \sum_{k=0}^{N-1} h(k).z^{-k}

Rappelons la définition de la transformée en z : z=esTe z= e^{sT_e}

s’agissant d’un filtre, nous nous intéressons au comportement de H en fonction de la fréquence.
Nous nous trouvons donc dans le cas particulier où s=jω s=j\omega
z=ejω.Te z= e^{j\omega.T_e}

Ainsi, la Transformée de Fourier Discrète ( TFD ) de H est :

H(f)=k=0N1h(k).ej2πfkTe \boxed{ H(f) = \sum_{k=0}^{N-1} h(k).e^{-j2\pi f k T_e } }


Revenons au cas idéal, pour lequel le nombre d’échantillons serait infini :

Hideal(f)=k=+h(k).ej2πfkTe H_{ideal}(f) = \sum_{k=-\infty}^{+\infty} h(k).e^{-j2\pi f k T_e }

Soit la fréquence normalisée fn=fFe f_n = \frac{f}{F_e}

Hideal(fn)=k=+h(k).ej2πfnk H_{ideal}(f_n) = \sum_{k=-\infty}^{+\infty} h(k).e^{-j2\pi f_n k}

Calcul de la Transformée de Fourier Discrète Inverse :

hideal(k)=1+1Hideal(fn)ej2πfnkdfn h_{ideal}(k) = \int_{-1}^{+1} H_{ideal}(f_n) e^{j2\pi f_n k }df_n

Considérons le filtre passe bas ideal :

  • Hideal(fn)=1 H_{ideal}(f_n) = 1 pour 0<fn<fcn 0 < f_n < fc_n avec fcn=FcFe fc_n=\frac{Fc}{Fe}
  • Hideal(fn)=0 H_{ideal}(f_n) = 0 pour fcn<fn<1 fc_n < f_n < 1

La réponse impulsionnelle est donc

hideal(k)=fcn+fcn1.ej2πkfndfn h_{ideal}(k) = \int_{-fc_n}^{+fc_n} 1. e^{j2\pi k f_n} df_n

hideal(k)=2.FcFesin(2πkFcFe)2πkFcFe=2.FcFesinc(2πkFcFe) \boxed{ h_{ideal}(k) = 2.\frac{F_c}{F_e} \frac{sin\left(2\pi k \frac{F_c}{F_e}\right)}{2\pi k \frac{F_c}{F_e}} = 2.\frac{F_c}{F_e} sinc\left(2\pi k \frac{F_c}{F_e}\right) }

h_ideale_1.svg

Réalisation du Filtre

Problème : cela correspond à une séquence infinie ( donc impossible à réaliser avec une boucle de calculs) .
Il faut donc tronquer n pour avoir une réponse finie ( on ne peut pas réaliser une multiplication-accumulation d’une infinité d’échantillons à chaque période ) .

Pour réaliser le filtre, nous allons multiplier hideal(n) h_{ideal}(n) par une fenêtre rectangle.
Cela aura bien sûr comme effet de dégrader l’effet du filtrage, mais c’est une condition nécéssaire pour le réaliser.

En effet La Transformée de fourier d’un produit de signaux temporels correspond à la convolution des tranformées de fourier de chacun des signaux.

h(n)=hideal(n).rect(n) h(n)=h_{ideal}(n).rect(n)

avec :

  • rect(n) = 1 si 0<nN1 0 < n \leqslant N-1
  • rect(n) = 0 si n>N1 n > N-1

Ex pour N=4:

h(n)=hd(0).1+hd(1).1+hd(2).1+hd(3).1+hd(4).0+hd(5).0+... h(n)=h_d(0).1 + h_d(1).1 + h_d(2).1 + h_d(3).1 + h_d(4).0 + h_d(5).0 + ...

ATTENTION : il s’agit bien ici d’une multiplication en temporel

h_ideale_2.svg

Ainsi, dans le cas où N=21, nous retenons 21 coefficients de hideal(n)=2.FcFesinc(2πnFcFe) h_{ideal}(n) = 2.\frac{F_c}{F_e} sinc\left(2\pi n \frac{F_c}{F_e}\right)

h_reelle_2.svg

Notons qu’en numérique la notion de causalité est relative, étant donné que l’on traite “N échantillons précédents mémorisés”

Influence de l’ordre du filtre

Bien évidemment, en augmentant l’ordre du filtre, les caractéristiques sont meilleures ( pente plus prononcée notamment ).
Attention toutefois, il faut pouvoir réaliser le calcul entre 2 acquisitions d’un échantillon d’entrée.

Influence du type de fenêtre

L’utilisation d’une fenêtre rectangulaire induit le phénomène de gibbs en fréquentiel, ce qui altère la caractéristique.

Il existe d’autres fenêtres pour des caractèristiques meilleures.

  • Fenêtre Rectangulaire ( scipy : boxcar ) : r(n)=1 r(n)= 1
  • Fenêtre de Hamming ( scipy : hamming ) : h(n)=0.540.46.cos(2πnN1) h(n)= 0.54-0.46.\cos(\frac{2\pi n}{N-1})
  • Fenêtre de Blackman ( scipy : blackman ) : b(n)=0.420.5.cos(2πnN)+0.08.cos(4πnN) b(n)= 0.42-0.5.\cos(\frac{2\pi n}{N})+0.08.\cos(\frac{4\pi n}{N})
from scipy.signal import get_window
import matplotlib.pyplot as plt

N = 1000

r = get_window("boxcar", N+1, fftbins=False)
plt.plot(r, label="boxcar")

h = get_window("hamming", N+1, fftbins=False)
plt.plot(h, label="hamming")

b = get_window("blackman", N+1, fftbins=False)
plt.plot(b, label="blackman")

plt.grid()
plt.xlabel("$n$")
plt.xlim([0, N])
plt.legend()
plt.show()

windows_fir.png

Considérons un filtre FIR Passe Bas d’ ordre 128, de fréquence de coupure fc=5000Hz, à une fréquence d' échantillonnage fs=44100Hz.

from scipy import signal
import matplotlib.pyplot as plt
import numpy as np
from scipy.signal import freqz, get_window

N = 128
fc = 5000
fs = 44100

n_vect = np.arange(N+1) #0 1 2 3 ...100

# Fenêtre Idéale
# !! dans numpy, sinc(x)=sin(pi.x)(pi.x)
h_d =2*(fc/fs)*np.sinc((2*(fc/fs))*(n_vect-int((N+1)/2)))

# Fenêtre Rectangulaire  
f, H = freqz(h_d, fs=fs)
h_dB = 20*np.log10(abs(H))

plt.subplot(211)
plt.plot(f,h_dB, label="rectangle")

h_Phase = np.unwrap(np.arctan2(np.imag(H),np.real(H)))
plt.subplot(212)
plt.plot(f,h_Phase, label="rectangle")

# Fenêtre de Hamming  
w = get_window("hamming", N+1, fftbins=False)
h=w*h_d
f, H = freqz(h, fs=fs)
h_dB = 20*np.log10(abs(H))
plt.subplot(211)
plt.plot(f,h_dB, label="Hamming")

h_Phase = np.unwrap(np.arctan2(np.imag(H),np.real(H)))
plt.subplot(212)
plt.plot(f,h_Phase,label="Hamming" )

# Fenêtre de Blackman  
w = get_window("blackman", N+1, fftbins=False)
h=w*h_d
f, H = freqz(h, fs=fs)
h_dB = 20*np.log10(abs(H))
plt.subplot(211)
plt.plot(f,h_dB, label="Blackman")

h_Phase = np.unwrap(np.arctan2(np.imag(H),np.real(H)))
plt.subplot(212)
plt.plot(f,h_Phase,label="Blackman" )

# MISE EN FORME DES TRACES

plt.subplot(211)
plt.ylabel("Module $H(e^{j\omega})\ en\ dB$")
plt.xlabel("Frequence [Hz]")
plt.xlim([0, fs/2])
plt.grid()
plt.legend()

plt.subplot(212)
plt.ylabel("Phase $H(e^{j\omega})\ en\ rad$")
plt.xlabel("Frequence [Hz]")
plt.xlim([0, fs/2])
plt.grid()
plt.legend()

plt.show()
#=========================================================

filtre_fir.png

REMARQUE : Le Filtre FIR a une Phase Linéaire, ainsi le filtre fir ne provoque pas de distorsion de phase.