S.S.I.I., 2013-14, n°6, Créer des filtres sur mesure pour compresser S.S.I.I., 2013-14, n°6, : Créer des filtres sur mesure pour compresser 1 Créer un

11S.S.I.I., 2013-14, n°6, S.S.I.I., 2013-14, n°6, : Créer des filtres sur mesure pour compresserCréer des filtres sur mesure pour compresser

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

0 fef

0 fef

0 fef

4

4

H1(f)4

0

0 fef

spectre (R échantillons)

R/4

R/8

3*R/16

En appliquant le principe du cours n° 5 pour compresser et décompresser le signal suivant dans un facteur 4 :


Quel est le taux de compression envisageable ici pour le signal dont le spectre est donné ci-dessous?

0 fef

spectre

0 fef

0 fef

0 fef

2

2

2


En fait, si on sait créer le filtre H2 ci-dessous, on atteint un taux de compression de 4 au lieu de 2 !

0 fef

spectre

0 fef

0 fef

0 fef

H24

4

4

Car la condition de Shannon généralisée est respectée : la largeur du spectre du signal étant inférieure à fe/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)

1

011110

R

k knkRnRnnn ehehehehs

1

0*

R

k knkk knk ehehehs


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 :

1

0

/2)(]1..0,[R

k

fkfik

TFDk

eehfHRkhh

ZffHfH

ZHH

ehRmfHH

RmHhfftH

e

Rmm

R

k

Rmkikem

m

),()(

,

)/(

]1..0,[)(1

0

/2

1. H = fft(h) est périodique de période R

2. h= ifft(H) est périodique, de période R

ZRTthth

Zhh

eHkThh

RkhHiffth

e

Rkk

R

m

Rmkimek

k

),()(

,

)(

]1..0,[)(1

0

/2

1

0

/2)(]1..0,[1 R

m

Rtmfim

TFDm

eeHthRmHH

3. En effet, les transformations fft et ifft sont presque identiques, au signe des exponentielles près. On vérifie que:

)()(]1..0,[ hiffthfftRkhh k


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-1car

Et 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

)./2cos(2

/2/2

/22/2

/)(2/2

RmkH

eHeH

eeHeH

eHeH

m

RmkimR

Rmkim

RmkikimR

Rmkim

RkmRimR

Rmkim

.)/2cos(212/

12/0

R

m mRk RkmHHHh

]1..0,[)(]1..0,[ RmEefftEetRkee mk

]1..0,[)(]1..0,[ RmSsfftSetRkss mk


La TFD du produit de convolution s= h*e est le produit cartésien des TFD de e et de h : 1..0, RkEHS kkk

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 à filtrer

C.Q.F.D.

]1..0,[ Rnss n

ehehsR

mmnmn *

1

0

]1..0,[ Rnee n

kk

R

m

R

v

Rkviv

Rkmim

R

m

R

n

Rmnkimn

Rkmim

R

n

R

m

Rmnkimn

Rkmim

R

n

RnkiR

mmnmk

k

EH

eeeh

eeeh

eeeh

eehS

RkSsfftS

1

0

1

0

/2/2

1

0

1

0

/)(2/2

1

0

1

0

/)(2/2

1

0

/21

0

]1..0,[)(

avec v= n-m quand n-m >0 et v=n-m+R quand n-m<0, puisque en-m=en-m+R.


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

1..0,

1..0,

*1

0

RkR

kfE

R

kfH

R

kfS

RkEHS

ehehs

eee

kkk

R

mmnmn


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)'])


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))")


Réponse fréquentielle du filtre de coefficients réels h=fftshift(real(ifft(H))) tracée sur M=1024 valeurs

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")

h = fftshift(real(ifft(H))) revient à permuter les deux moitiés du vecteur h


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)];


Exemple : calcul et d’utilisation du filtre H2 avec Scilab

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

// filtre passe bande 1000Hz-2000Hz// gain 4, R=64, fe=8000HzR=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

Documents

S.S.I.I., 2013-14, n°6, Créer des filtres sur mesure pour compresser S.S.I.I., 2013-14, n°6, : Créer des filtres sur mesure pour compresser 1 Créer un