Créer un filtre de réponse fréquentielle donnée

1. spectre (R échantillons) R/4 3*R/16 R/8 f 0 fe f 0 fe f 0 fe H1(f) 4 0 f 0 fe Créer un filtre de réponse fréquentielle donnée Jean-Paul Stromboni, Polytech'Nice Sophia, S.I. 3ème année Cours n° 6, novembre 2013, durée : 50 mn, avec vidéoprojecteur En appliquant le principe du cours n° 5 pour compresser et décompresser le signal suivant dans un facteur 4 :

2. spectre f 0 fe f 0 fe f 0 fe 2 f 0 fe Quel est le taux de compression envisageable ici pour le signal dont le spectre est donné ci-dessous?

3. spectre f 0 fe f 0 fe f 0 fe 4 H2 f 0 fe En fait, si on sait créer le filtre H2 ci-dessous, on atteint un taux de compression de 4 au lieu de 2 ! Car la condition de Shannon généralisée est respectée : la largeur du spectre du signal étant inférieure à fe/4 !!

4. Pour réaliser H2 à l’aide d’un filtre programmé, on peut utiliser un filtre non récursif dont voici l’équation : • Le vecteur e=(en, n=0..N-1) contient les N échantillons du signal à filtrer, ou entrée du filtre • Le vecteur h=(hm, m= 0..R-1) contient les R coefficients du filtre (coefficients réels) • Le vecteur s=(sn, n=0..N-1) contient la sortie du filtre ou signal filtré, chaque valeur sn est calculée par une itération de l’équation ci-dessus • R est la taille du filtre • L’équation est un produit de convolution (symbole ‘*’): • h contient la réponse impulsionnelle du filtre, c’est-à-dire que s= h pour une entrée impulsion (e0=1, en=0 si n!=0). • Pour en savoir plus: • Il s’agit d’un filtre linéaire et stationnaire, en anglais Linear Time Invariant. Pour mieux comprendre : • sn=en+ en-1 est linéaire et stationnaire • sn= sin (en-1) est non linéaire • sn=en+n en-1 est non stationnaire • Il s’agit d’un filtre non récursif, ou à réponse impulsionnelle de durée finie (Finite Impulse Response FIR en anglais) : • Par contre, sn=sn-1+en-1, est un filtre récursif (ou IIR)

5. Pour calculer la réponse fréquentielle d’un filtre récursif, on note les propriétés suivantes de la TFD directe et inverse : • H = fft(h) est périodique de période R • h= ifft(H) est périodique, de période R • En effet, les transformations fft et ifft sont presque identiques, au signe des exponentielles près. On vérifie que:

6. Pour calculer la réponse fréquentielle d’un filtre récursif, on note les propriétés suivantes de la TFD directe et inverse : • Si le vecteur H est réel, soit Hm réel pour m=0..R-1, à quelle condition le vecteur h=ifft(H) est il réel ?Réponse : il suffit que Hm=HR-m, pour m=0..R-1carEt par conséquent hk est réel, pour k=0..R-1: • noter que ce que l’on vient d’établir pour h et H, est vrai également pour en particulier ek est périodique de période R, • et pour • On utilisera dans la suite e et E, h et H et s et S ainsi définis, et de tailles R

7. La TFD du produit de convolution s= h*e est le produit cartésien des TFD de e et de h : Voici la démonstration, qui utilise la périodicité de la TFD inverse.Soit l ’équation du filtre Soit le signal filtréet le signal à filtrerC.Q.F.D. avec v= n-m quand n-m >0 et v=n-m+R quand n-m<0, puisque en-m=en-m+R.

8. Voici donc l’effet sur le spectre d’un filtre non récursif • E contient le spectre du signal à filtrer • S contient le spectre du signal filtré • H= fft(h) est la réponse fréquentielle du filtre dont les coefficients réels sont dans le vecteur h • La réponse fréquentielle H est la transformée de Fourier de la réponse impulsionnelle h • Si on fixe H à la valeur désirée H2, h= ifft(H) calcule les coefficients du filtre de réponse fréquentielle H=H2 • Et ces coefficients seront bien réels si on a pris la précaution de choisir Hm= HR-m, m=0..R-1 • Attention ! on impose uniquement R valeurs sur la réponse fréquentielle, aux fréquence kfe/R, k=0..R-1, il faudra vérifier H(f) entre ces fréquences

9. Vérifier avec Scilab que H est symétrique par rapport à R/2, et que imag(ifft(H)) = 0 et donc h=real(ifft(H)) fe=8000; R=64; // H symétrique par rapport à R/2 H=4*[zeros(1,R/8),ones(1,1+R/8),zeros(1,-1+R/2),ones(1,1+R/8),zeros(1,-1+R/8)]; //étude de h=ifft(H) h=ifft(H); t=[0:R-1]/fe; plot2d(t',[real(h'),imag(h')]) e=gce(); e.children(1).thickness=3; xgrid(); xtitle("vérification: imag(ifft(H))=0",... "temps (s)","donc h=real(ifft(H))") h1=legend(['real(h)';'imag(h)'])

10. Réponse fréquentielle du filtre de coefficients réels h=real(ifft(H)) tracée sur M=1024 valeurs au lieu de 64 • Pour tracer la réponse fréquentielle du filtre de coefficients h=real(ifft(H)), il suffit de tracer abs(fft(h)) • Pour tracer M=16*R valeurs au lieu de R, il suffit d’aug-menter le vecteur h de 15*R coefficients nuls : M=16*R; fe=8000; fM=[0:M-1]*fe/M; h=real(ifft(H)); hM=zeros(1,M); hM(1:R)=h; plot2d(fM,abs(fft(hM))) xgrid(); xtitle(["tracé de … h=real(ifft(H))) sur",string(M),"points"] ... ,"fréquence (Hz)","H=abs(fft(h))")

11. Réponse fréquentielle du filtre de coefficients réels h=fftshift(real(ifft(H))) tracée sur M=1024 valeurs h = fftshift(real(ifft(H))) revient à permuter les deux moitiés du vecteur h clf(); M=16*R; h=fftshift(ifft(H)); fM=[0:M-1]*fe/M; hM=zeros(1,M);hM(1:R)=real(h); plot2d(fM,abs(fft(hM))) xgrid(); xtitle(["tracé de abs(fft(fftshift(real(ifft(H))))),", string(M),…" points"],"fréquence (Hz)","H")

12. Il suffit d’arrondir la forme de la réponse fréquentielle spécifiée dans le vecteur H pour atténuer les oscillations résiduelles. H=4*[zeros(1,R/8-1),0.1,0.5,0.9,ones(1,R/8-3),0.9,0.5,0.1, ... zeros(1,-3+R/2),0.1,0.5,0.9,ones(1,R/8-3),0.9,0.5,0.1, ... zeros(1,-2+R/8)];

13. La réponse fréquentielle du filtre est définie dans le vecteur H, les coefficients du filtre sont calculés dans le vecteur h, on filtre ‘piano.wav’, on compare spectrogrammes et énergies avant et après filtrage Exemple : calcul et d’utilisation du filtre H2 avec Scilab // filtre passe bande 1000Hz-2000Hz // gain 4, R=64, fe=8000Hz R=64; fe=8000; n=0:R-1; fr=n*fe/R; H=4*[zeros(1,R/8),ones(1,1+R/8), ... zeros(1,-1+R/2),ones(1,1+R/8), ... zeros(1,-1+R/8)]; plot2d3(fr,H) xgrid xtitle(['H2,avec R=',string(R)], ... 'fréquence (Hz)’, ‘H’) //calcul des coefficients du filtre h=fftshift(real(ifft(H))); plot2d3(n/fe,h) xtitle('coefficients du filtre',... 'temps (s)',... 'h=fftshift(real(ifft(H)))') xgrid(); // filtrage [y,fe]=wavread('piano.wav'); disp(fe) // fe=8000 sound(y,fe) yf= convol(h,y); wavwrite(yf,fe,'pianofilt.wav') sound(yf,fe) //Spectrogrammes (Goldwave) // énergie Ey=(y*y')/2 // énergie y = 163.96 Eyf=(yf*yf')/2 // énergie yf =89.62