1re Anne

Informatique
2011-2012
Bureaux dtudes
Traitement Numrique du Signal
Digital Signal Processing
M. Frikel - J. Fadili - Pierre-Emmanul Morant
GREYC, UMR 6072 CNRS, 6, Boulevard Marchal Juin, 14050 Caen Cedex
Table des matires
1 BE I : Echantillonnage de signaux et transforme de Fourier discrete
(TFD) 4
1.1 But du Bureau detude : . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
1.2 Travail prparatoire : . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
1.3 Travail raliser : . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
2 BE II : La Transforme de Fourier Discrte (DFT) 8
2.1 Partie I : . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
2.2 Partie II : . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
3 BE III : Filtres discrets et quation de rcurrence 11
3.1 But Du bureau dtudes . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
3.2 Travail prparatoire . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
3.2.1 Filtres tudis : . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
3.2.2 Travail prparatoire : . . . . . . . . . . . . . . . . . . . . . . . . . . 12
3.3 Etude avec Matlab . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
3.3.1 Vrication du travail prparatoire : . . . . . . . . . . . . . . . . . 13
3.3.2 Fonctionnement de lquation rcurrente : . . . . . . . . . . . . . . 13
3.4 Etude avec Simulink . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
3.4.1 Acquisition et visualisation des signaux avec Simulink : . . . . . . 13
3.4.2 Filtres lmentaires : . . . . . . . . . . . . . . . . . . . . . . . . . . 15
3.5 Annexe . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
4 BE IV : Mise en oeuvre dun ltrage numerique Pour supprimer le
parasite dun signal sonore 18
4.1 But de la manipulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
4.2 Analyse du signal . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
4.2.1 Matriel utilis : . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
4.2.2 Travail raliser : . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
4.3 Conception du ltre . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
4.3.1 Filtrage avec un ltre RII de Butterworth : . . . . . . . . . . . . . . 22
4.3.2 Filtrage avec un ltre RII de Chebyshev : . . . . . . . . . . . . . . 24
4.3.3 Filtrage avec un ltre RIF : . . . . . . . . . . . . . . . . . . . . . . 25
2
5 BE V : Gnration et dtection des signaux DTMF : Dual Tone Multi-
Frequency 26
5.1 Le problme . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
5.2 generation de signaux discrets : ltre generateur (ou formeur) . . . . . . . 27
5.2.1 Cas gnral : . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
5.2.2 Cas des signaux sinusodaux : . . . . . . . . . . . . . . . . . . . . . 29
5.2.3 Ralisation : Gnration de signaux sinusodaux : . . . . . . . . . . 29
5.2.4 Gnration dun signal DTMF : . . . . . . . . . . . . . . . . . . . . 31
5.3 Dtection - algorithme de goertzel . . . . . . . . . . . . . . . . . . . . . . . 32
5.3.1 Dtection dune sinusode : . . . . . . . . . . . . . . . . . . . . . . . 32
5.3.2 Rappel sur la TFD : . . . . . . . . . . . . . . . . . . . . . . . . . . 32
5.3.3 Algorithme de Goertzel : . . . . . . . . . . . . . . . . . . . . . . . . 33
5.3.4 Implmentation de lalgorithme de Goertzel : . . . . . . . . . . . . . 35
5.3.5 Ralisation : Dtection des frquences DTMF : . . . . . . . . . . . . 36
6 BE VI : Filtres rponse impulsionnelle innie RII : Transposition
Continu-Discret 38
6.1 Filtre rjecteur . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
6.2 Etude dun ltre du 2me ordre . . . . . . . . . . . . . . . . . . . . . . . . 40
6.2.1 Fonction de transfert dun ltre continu . . . . . . . . . . . . . . . . 40
6.2.2 Fonction de transfert dun ltre discret . . . . . . . . . . . . . . . . 41
6.2.3 Etude dun ltre du 3me ordre . . . . . . . . . . . . . . . . . . . . 43
7 BE VII : Dbruitage dun lectrocardiogramme (ECG) par ltrage 45
7.1 Analyse dun enregistrement dun ECG : . . . . . . . . . . . . . . . . . . . 45
7.2 Filtrage du bruit : . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
3
Chapitre 1
BE I : Echantillonnage de signaux et
transforme de Fourier discrete (TFD)
1.1 But du Bureau detude :
Mettre en pratique lutilisation de la TFD dans les deux cas suivants :
Estimer la transforme de Fourier dun signal
Raliser lanalyse spectrale (mise en vidence de raies) de signaux inconnus : analyseur
de spectre. indications matrielles
Ltude sera eectue avec le logiciel Matlab pour cela quelques normes sont respecter :
Utilisation dun chier "script" :
Les direntes tapes du travail seront ralises par lintermdiaire dun chier "script"
(.m file). Pour partir sur des bases saines, il est recommand de commencer ce chier
par les instructions clc (clear command qui rinitialise la fentre de commande), clear
(qui eace les variables stockes dans lespace de travail) et close all (qui eace toutes
les gures). Les gures ventuellement seront sauvegardes dans des chiers.
Chargement des donnes :
Les donnes sont fournies dans des chiers de donnes dextension .mat, ces chiers
contiennent la fois les mesures sous forme dun tableau de valeurs {x
k
} et la fr-
quence dchantillonnage utilise note F
s
. Par exemple, pour le chier de donnes
sig1_rect.mat :
Le chargement est ralis par la commande load sig1_rect.
Dans lespace de travail (workspace) sont gnres deux variables : F
s
la frquence dchan-
tillonnage du signal et x
k
le tableau de points rsultant de lchantillonnage.
Visualisation graphiques :
Comme nous travaillons en chantillonn, il sera prfrable dutiliser la commande stem
plutt que la commande plot. Comme sur toute commande de Matlab, une aide en ligne
est disponible avec la commande help.
4
Elle est calculable avec la commande fft. Lorsque le nombre de points est une puissance
de 2, alors lalgorithme de transforme de Fourier rapide est utilis, autrement il sagit
dun calcul classique des exponentielles avec sinus et cosinus et, pour le travail de ce BE,
cela na aucune importance.
La lecture de laide en ligne de Matlab est fortement conseille. Il est intressant de
complter ces travaux par la ralisation personnelle dun algorithme de FFT mais cela
relve dun enseignement dinformatique, algorithmique ou gnie informatique.
1.2 Travail prparatoire :
1. Rappeler la transforme de Fourier de rect(t). En dduire celle de rect(t/) puis
celle de x(t) = rect((t t
0
)/). On posera TF[x(t)] = X(f).
2. Le signal x(t) est chantillonn la priode Ts. Avec la courbe 1 ci-dessous calculer
la relation entre Ts et pour que la condition de Shannon soit remplie 10% puis
5%.
3. Lchantillonnage de x(t) dure un temps T > et fourni N chantillons : T = N.T
s
.
A partir des N chantillons {x
k
} nous calculons la transforme de Fourier discrte
5
{X
n
}. Combien dchantillons comporte {X
n
} ? Quelle est la relation entre X
n
et
X(f) ?
1.3 Travail raliser :
Etude du chier :
1. Charger les donnes du signal sig1_rect, le reprsenter graphiquement en talon-
nant correctement laxe des abscisses. En dduire T, T
s
et . La condition de Shannon
a-t-elle t respecte ?
2. calculer X
n
la TFD de x
k
.
3. Reprsenter graphiquement son module en talonnant correctement laxe des abs-
cisses en frquences. Nous pouvons nous contenter de la moiti des points : pourquoi ?
4. Comparer la valeur thorique.
Interpolation de la TFD :
1. Nous ne disposons pas de susamment de points dans lestimation de X(f). Pour
cela nous allons procder la manipulation suivante :
Constituer un nouveau signal x
k1
partir de x
k
auquel on ajoute M.N zros. (on
utilisera M=10)
2. Calculer la TFD X
n1
de ce nouveau signal.
3. La reprsenter sur la mme gure que X
n
(commande hold on).
4. Comparer les rsultats et interprter la manipulation eectue.
Analyse spectrale dun signal par la TFD :
Signal sinusodal : Le signal tudi est contenu dans le chier sig1_quidonc.
Analyse spectrale lmentaire :
1. Charger les donnes du chier . Reprsenter le signal en talonnant correctement
laxe des abscisses. Quelles indications pouvez vous en retirer ?
2. Calculer la TFD X
n
du signal x
k
et reprsenter son module. Quen dduisez vous ?
Amlioration de la prcision frquentielle :
1. Compte tenu de ltude prcdente, construire un signal x
k1
constitu du signal x
k
auquel on concatnera M.N zros (M=10).
2. Calculer X
n1
la TFD de x
k1
.
3. Reprsenter simultanment (commande hold on) les modules de X
n
et X
n1
en talon-
nant correctement laxe des abscisses (utiliser des couleurs direntes et commencer
par reprsenter X
n1
qui, ayant plus de points, risque de masquer X
n
).
6
4. Pouvez-vous en dduire une forme analytique du signal tudi ?
Intrt dune fentre de pondration :
1. Crer le signal x
kh
en pondrant x
k
par une fentre de hanning (commande hanning)
2. Reprendre lanalyse prcdente avec ce nouveau chier soit :
3. Avec la commande fft, calculer la TFD X
nh
du signal x
kh
.
4. Construire un signal xkh1 constitu du signal xkh auquel on concatnera M.N zros.
5. Calculer X
nh1
la TFD de x
kh1
.
6. Sur une troisime gure, reprsenter simultanment les modules de X
nh
et X
nh1
en
talonnant correctement laxe des abscisses
7. Pouvez-vous en dduire une forme analytique du signal tudi ?
8. Etudier la dirence en plaant une fentre de Hamming.
Signal inconnu :
1. Reprendre lanalyse prcdente avec le signal dont les donnes sont contenues dans
le chier sig2_quidonc. Pouvez-vous fournir une expression analytique du signal ?
2. Inuence de la priode dchantillonnage
Le mme signal est chantillonn la priode F
s1
=
F
s
4
et ces donnes sont contenues
dans le chier de donnes sig3_quidonc.
3. Quel est le rsultat de lanalyse spectrale ralise avec ces donnes ? Justier les
dirences observes par rapport lanalyse prcdente.
4. Ce problme naurait jamais d intervenir exprimentalement. Pourquoi ?
7
Chapitre 2
BE II : La Transforme de Fourier
Discrte (DFT)
Objectif : Comprendre le principe de la Transforme de Fourier Discrte (DFT), connatre
certaines de ses proprits et comparer sa complexit avec des algorithmes rapides (fft).
2.1 Partie I :
Montrez que la transforme de Fourier Discrte sur N points, dnie par :
X(k) =
N1

n=0
x[n]e
(j
2
N
)nk
(2.1)
peut scrire comme le produit dune matrice W, de taille NxN, avec le vecteur colonne
x[n], o llment W(n, k) de la matrice W est :
W(k, N) = e
(j
2
N
)nk
(2.2)
(k est le numro de la ligne, qui vaut entre 0 et N 1, et n est le numro de la colonne,
qui vaut aussi entre 0 et N 1). Faites en particulier le cas N = 4, et montrez que la
matrice W est donne par :
_
_
_
_
_
1 1 1 1
1 j 1 j
1 1 1 1
1 j 1 j
_
_
_
_
_
Dans MATLAB, calculez la transforme de Fourier de
x = [1 3 1 2]
8
(Notation MATLAB) en faisant le produit matriciel Wx
T
(T dnote la transposition) puis
en calculant la fft de x avec fft(x).
Vous devriez obtenir le mme rsultat.
Pour obtenir la Transforme inverse (i.e. retrouver x[n] partir de X(k)), on applique
lquation de la transforme inverse :
x[n] =
1
N
N1

k=0
X(k)e
(j
2
N
)nk
(2.3)
qui peut aussi scrire sous forme matricielle. crivez cette matrice inverse pour la cas N
= 4, et appliquez-la la transforme de Fourier du signal x[n] plus haut - vous devriez
retrouver le signal x[n].
Comparez avec la fonction ifft.
2.2 Partie II :
Il sagit ici de gnrer une sinusode pure, de frquence donne, partir de la Transforme
de Fourier inverse. Notez dabord quil faut 2 exponentielles complexes pour gnrer une
sinusode (ici, un cosinus) :
cos[n] =
e
jn
+ e
jn
2
Notez aussi que la frquence des 2 exponentielles complexes requises est la mme, mais
de signe inverse et .
Notez que lexponentielle complexe : e
(j
2
N
)n(Nk)
.
est quivalente lexponentielle complexe : e
(j
2
N
)nk
.
On tire de tout ceci quen mettant deux termes non-nuls dans la transforme de Fourier
discrte, le premier k et le second (son complexe conjugu) N k, on obtient une
sinusode discrte de frquence normalise :
=
2
N
k
Dans MATLAB, en prenant par exemple N = 64 et k = 3 :
N = 64 ;
k = 3 ;
X = [zeros(1,k) N/2 zeros(1,N-(2*k-1)-2) N/2 zeros(1,k-1)] ;
x = real(ifft(X)) ;
plot(x) ;
Le rsultat est un cosinus pur qui fait exactement k = 3 cycles en N = 64 chantillons.
9
tudiez bien le code MATLAB ci-dessus, notamment la ligne qui gnre les coecients X(k)
(troisime ligne), puis ralisez ce qui suit.
En dnissant dabord les coecients de Fourier X(k) appropris tel que montr dans le
code MATLAB ci-dessus, gnrez par transforme inverse (ifft) les signaux x[n] suivants :
1. un cosinus pur de dure 1 seconde, de frquence 500 Hz et damplitude 1, pour
un systme chantillonn 8000 chantillons/seconde (vous devez donc convertir
la dure en nombre dchantillons, et la frquence (500 Hz) en une valeur k qui
dtermine la position des raies non-nulles dans le spectre).
2. un sinus pur de dure
1
2
seconde, de frquence 800 Hz et damplitude 1, pour un
systme chantillonn 8000 chantillons/seconde.
10
Chapitre 3
BE III : Filtres discrets et quation de
rcurrence
3.1 But Du bureau dtudes
Au niveau de la programmation, un ltre numrique est une simple quation rcurrente.
En utilisant lenvironnement de Matlab pour la ralisation de programmes dtude et
celui de Simulink avec la carte son E/S pour les manipulations, nous allons vrier les
proprits et le fonctionnement de quelques ltres lmentaires.
3.2 Travail prparatoire
3.2.1 Filtres tudis :
Nous allons nous intresser 4 types de ltres :
1. Le drivateur numrique :
H(z) =
1 z
1
T
s
(3.1)
2. Filtre intgrateur :
p
T
s
1 z
1
(3.2)
Avec la transposition prcdente, la transposition numrique de lintgrateur est
ralise par :
Remarque : pour interprter les rsultats, bien rchir laction dun tel ltre sur
un signal ayant une composante continue.
11
3. Moyenne mobile :
Cest un ltre numrique qui ralise la moyenne pondre des k dernires mesures
eectues (ltre MA pour Mobile Average). Sa fonction de transfert chantillonne
est :
H(z) = 1 + z
1
+ z
2
+ ... + z
(k1)
(3.3)
La moyenne mobile tudie sera prise avec 5 points ( k = 5 )
4. Filtre notch lmentaire :
Il a pour fonction de supprimer une pulsation particulire
0
. Il comprend deux
ples complexes conjugus et deux zros associs :
Zeros e
j
0
T
s
Poles re
j
0
T
s
avec r < 1 (3.4)
ce qui donne la fonction de transfert :
H(z) =
1 2cos(
0
T
s
).z
1
+ z
2
1 2.r.cos(
0
T
s
).z
1
+ r
2
.z
2
(3.5)
Dans la manipulation on choisira f
0
600Hz et la manipulation sera faite pour deux
valeurs de r : r = 0.9 et r = 0.7.
3.2.2 Travail prparatoire :
Pour tous les ltres tudis faire :
1. La reprsentation des ples et zros
2. Ltude frquentielle sommaire
3. Lcriture de lquation rcurrente qui permettra dimplanter le ltre
3.3 Etude avec Matlab
Une premire tude de ces ltres sera ralise avec un chier "script" de Matlab (.m
file) dj utilis lors du BE n1 ( Voir le texte correspondant pour les prcautions
lmentaires ). Ce chier sera eectuer selon lossature suivante :
12
3.3.1 Vrication du travail prparatoire :
1. Dnir une fonction de transfert associe au ltre (tf)
2. Eectuer la reprsentation des ples et zros (pzmap) et du module de la rponse
frquentielle (bode).
3.3.2 Fonctionnement de lquation rcurrente :
Complter le chier par :
1. Charger un chier de donnes : voir le tableau ci-dessous.
2. Filtrer ces donnes grce lquation rcurrente (boucle for) et comparer le rsultat
celui fourni par linstruction filter de Matlab
Drivateur Intgrateur la moyenne mobile Filtre notch
Sinus600.mat Sinus600.mat Sinus600.mat Sinus300.mat
Sinus1800.mat Sinus1800.mat Sinus1800.mat Sinus600.mat
Carre600.mat Carre600.mat Carre600.mat Sinus1800.mat
Scie600.mat Scie600.mat Scie600.mat
3.4 Etude avec Simulink
Bien imprgns du secret des fonctions lmentaires de Matlab, le travail peut tre repris
avec une interface utilisateur graphique : Simulink. La partie matrielle dacquisition des
signaux peut tre faite avec le DSP de la carte son de lordinateur ce qui vite, dans cette
phase de dcouverte, la mise en oeuvre de moyens techniques plus ddis mais plus lourds
comme les cartes E/S, DSP, processeurs,... utiliss usuellement pour ce type dapplications.
3.4.1 Acquisition et visualisation des signaux avec Simulink :
Nous utiliserons de prfrence les blocs prsents dans la bibliothque Signal Processing
Blockset. Quelques lments de base sur lutilisation de Simulink sont rappels en an-
nexe.
Achage dun signal :
Ralise avec le bloc Signal Processing Blocksets > Signal processing sinks >
Time Scope. Sa prise en main est aise, il est possible dy dnir plusieurs voies. Matlab
fournit ventuellement une aide en ligne
Lecture dun signal avec la carte son :
Le signal mesurer doit tre plac sur lentre line du PC (che arrire bleue) et peut
provenir dun gnrateur de fonctions en utilisant un cble BNC-jack3.5. Il est ncessaire
13
de sassurer, dans la conguration du PC, que cette entre est bien active et non lentre
du microphone comme cest souvent le cas par dfaut.
Dans Simulink, lacquisition sera ralise par le bloc Signal Processing Blocksets >
Signal Processing Sources > From Wave Device. Ce bloc permet dactiver le conver-
tisseur analogique-numrique de la carte son. Ses rglages sont les suivants :
Priode dchantillonnage : nous choisissons "User-defined" = dnie par lutilisateur.
Sa valeur sera F
s
quil faudra donc dnir au pralable dans lenvironnement Matlab.
Un codage 16 bits
Une utilisation dun buer tampon intermdiaire : ici une "frame" de 64 chantillons
est susante.
Un retard de pas trop lev, 0.2s est une valeur raisonnable.
Le type des donnes est double
Raliser le schma ci-dessous. Connecter un gnrateur fournissant un signal sinusodal
denviron 600 Hz, choisir une frquence dchantillonnage Fs=8000Hz et visualiser le signal
au Scope. La dbuerisation est utile pour robtenir le signal chantillonn : bloc Signal
14
Processing Blocksets > Signal Management > Buffers.
Pour raliser une manipulation correcte il faut :
Utiliser des signaux dentre respectant la norme wav cest dire que leur valeur doit
tre comprise entre -0,5V et +0,5V
Paramtrer correctement la simulation soit :
Un temps de simulation pas trop long.
Nous sommes en chantillonn pur, il ny a pas de systme continu et il est inutile
dutiliser un solveur dquations do : discrete (non continuous states).
Utiliser donc un pas xe : Fixed-step
Ce pas xe sera la priode dchantillonnage :
1
F
s
.
Signal provenant dun chier enregistr :
Ces chiers sont la norme wav cest dire quils respectent lamplitude pour les donnes
et contiennent la frquence dchantillonnage laquelle ils ont t crs.
Raliser le schma ci-dessous et visualiser le signal qui doit tre priori une sinusode de
frquence 1000Hz.
3.4.2 Filtres lmentaires :
Lquation rcurrente simplante sous Simulink avec le bloc :
Signal Processing Blocksets > Filtering > Filter Designs > Digital
Filter.
En utilisant ce bloc et une entre avec un gnrateur basse frquence, vrier les rsul-
tats dj observs avec Matlab. On pourra en particulier dans le cas du ltre drivateur
regarder son action sur des signaux de forme carre et triangulaire.
3.5 Annexe
Dmarrage de Simulink
15
Cration dun chier de travail
Ouverture dun chier de travail dj existant
Bibliothque
Blocs de la bibliothque rangs par fonctions
16
17
Chapitre 4
BE IV : Mise en oeuvre dun ltrage
numerique Pour supprimer le parasite
dun signal sonore
4.1 But de la manipulation
Nous disposons dun signal musical sur lequel est prsent un parasite (siement). Le but
est de dbarrasser ce signal du parasite en procdant de la manire suivante :
1. Identier la ou les frquences du parasite en utilisant une FFT fentre glissante.
2. Concevoir un ltre dont la rponse en frquence supprime le parasite en altrant le
moins possible le signal musical.
3. Vrier son fonctionnement avec Simulink.
4.2 Analyse du signal
La premire tape consiste analyser le signal pour dterminer la ou les frquences para-
sites. Lanalyse frquentielle se fait lanalyseur de FFT de Simulink qui permet une analyse
dynamique par fentre glissante puisque le signal varie dans le temps.
4.2.1 Matriel utilis :
Nous utiliserons Matlab mais pour intgrer ltude temporelle il nous est ncessaire dutili-
ser sa "toolbox" Simulink. Ce logiciel comporte une bote outils DSP Blockset (Digi-
tal Signal Processing Traitement numrique du signal) qui fournit les blocs de traitement
de signal pour des applications orientes vers lutilisation de processeurs spcialiss.
Dans cette bibliothque deux lments sont particulirement intressants pour notre ap-
plication :
18
1. Spectrum scope : analyseur de spectre en temps rel.
2. From Wave File : lecture en temps rel dun chier son au format .wav
3. Dautres lments peuvent tre intressants :
4. Vector scope : reprsentation des signaux temporels.
5. To Wave Device : envoi du signal sur les hauts parleurs
4.2.2 Travail raliser :
Avec Simulink, raliser lapplication Analyse.mdl qui en utilise les trois blocs : From
Wave File (source), Spectrum scope (analyse frquentielle) et To Wave Device
(audition du fichier).
Paramtrage du bloc From Wave File :
Nom du chier .wav utilis
Les donnes sont fournies par "frame" ou par paquets. Ici une "frame" de 256 chan-
tillons.
19
La lecture du chier se fait par paquets, il est lgant de placer ici le mme nombre
dchantillons que dans la "frame"
Type de donnes manipules, ici obligatoirement type double
Paramtrage du bloc Spectrum Scope :
Il y a quatre onglets :
Proprits du Scope, essentiel au bon fonctionnement de lanalyse
Proprits de lachage : la prsentation du rsultat
Proprits des axes : choix des chelles
Proprits des lignes : style (pointill, continu ...), couleur ....
Paramtrage du bloc To Wave Device :
Ayant notre disposition des ordinateurs PC standards, le paramtrage sera proche de
celui indiqu sur la gure ci-contre.
Pour ce TP :
On utilisera le chier le chier Lockwood2.wav qui contient de la musique bruite chan-
tillonne 8 kHz. Choisir et justier la longueur de FFT ainsi que le nombre de moyennes
spectrales (Spectrum Scope). Choisir les paramtres dachage et des axes. En ralit,
ceux-ci seront choisis en faisant quelques essais de simulation pour obtenir un achage
confortable. Pour les lignes on peut laisser lachage par dfaut (continu, noir).
Les paramtres de la simulation sont choisis de manire standard avec un temps de simu-
lation denviron 60s (ce qui correspond peu prs la longueur de lenregistrement) ou
linni mais dans ce cas il faudra penser larrter manuellement. Pour viter les pertes
de temps, puisque nous ne traitons que des signaux discrets, le solveur dquations peut
20
tre paramtr pas xe et avec uniquement des tats discrets.
Dans le mme souci, il est prfrable de supprimer toutes les sauvegardes de donnes dans
lespace de travail Data Import/Export.
Eectuer des essais de simulation et aner les rglages : achage, axes, choix de la
longueur de FFT,...
Les rsultats :
En choisissant un achage non persistant (temps rel) et coutant le signal, reprer le
spectre de la perturbation.
Refaire la manipulation avec un achage persistant (les courbes successives sont mmo-
rises sur lcran). Vrier que la perturbation est distincte du signal et dterminer la
bande de frquence qui la concerne.
21
4.3 Conception du ltre
Le ltre concevoir est un ltre passe-bas qui doit supprimer les frquences de la perturba-
tion. Pour cela, Matlab fournit un grand nombre doutils performants mais, pour illustrer
la dmarche, nous allons concevoir " la main" notre propre ltre tout en utilisant les
outils danalyse de Matlab.
Les dirents calculs de mise au point seront mens en utilisant un cher script de Matlab :
un chier Passe_Bas.m crer de toutes pices.
4.3.1 Filtrage avec un ltre RII de Butterworth :
Ordre du ltre :
La premire question qui se pose est : quel est lordre du ltre qui correspond notre
problme ? La rponse peut tre trouve par essais successifs de ltres de plus en plus
complexes cependant, une rponse peut tre trouve de manire systmatique programme
dans la commande Matlab buttord.
[N, Wn] = BUTTORD(Wp, Ws, Rp, Rs)
Laide en ligne nous indique que cette commande retourne :
N : le nombre minimal dchantillons du ltre
22
W
n
: la frquence de coupure 3dB du ltre
Pour les charges suivantes :
Une attnuation de R
p
(en dB) pour la frquence W
p
en bande passante
Une attnuation de R
s
(en dB) pour la frquence Ws en bande coupe
Toutes les frquences sont normalises par la demi frquence dchantillonnage (F
s
/2).
Sur la manipulation prcdente, dterminer la frquence Wp pour laquelle nous souhai-
tons une attnuation infrieure 10% et la frquence Ws pour laquelle lattnuation sera
suprieure 90%. En dduire lordre et la frquence de coupure 3dB du ltre.
Calcul des coecients du ltre :
Ceux-ci peuvent tre directement calculs par la commande :
[B,A] = BUTTER(N,Wn)
o B et A contiennent respectivement les coecients du numrateur et du dnomina-
teur du ltre de Butterworth dordre N (N+1 coecients) et de frquence de coupure
normalise W
n
.
Vrication :
A laide de la commande Bode, calculer la rponse en frquence du ltre et reprsenter
23
son module sur une gure.
Modier le chier Simulink analyse.mdl en y incluant le bloc de ltrage : Signal
Processings>Filter Designs>Digital Filter pour obtenir un chier Filtrage.mdl
et vrier que le ltrage se ralise correctement.
4.3.2 Filtrage avec un ltre RII de Chebyshev :
La mme manipulation est reprise avec une ralisation de Chebyshev
Ordre du ltre :
Une rponse peut tre trouve de manire systmatique programme dans la commande
Matlab cheb1ord :
[N, Wp] = CHEB1ORD(Wp, Ws, Rp, Rs)
Laide en ligne nous indique que cette commande retourne :
N : le nombre minimal dchantillons du ltre
W
p
: la frquence de coupure 3dB du ltre
Pour les charges suivantes :
24
Une attnuation de R
p
(en dB) pour la frquence W
p
en bande passante
Une attnuation de R
s
(en dB) pour la frquence W
s
en bande coupe
Toutes les frquences sont normalises par la demi frquence dchantillonnage (Fs/2).
Sur la manipulation prcdente, dterminer la frquence Wp pour laquelle nous souhai-
tons une attnuation infrieure 10% et la frquence Ws pour laquelle lattnuation sera
suprieure 90%. En dduire lordre et la frquence de coupure 3dB du ltre.
Calcul des coecients du ltre :
Calculs par la commande :
[B,A] = CHEBY1(N,R,Wn)
o B et A contiennent respectivement les coecients du numrateur et du dnominateur
du ltre de Butterworth dordre N (N+1 coecients) et de frquence de coupure norma-
lise Wn. Par rapport au ltre de Butterworth apparat le paramtre R qui est le taux
dondulation maximal de la rponse en frquence en bade passante. Une valeur de lordre
de 2% est juge acceptable.
Vrication :
A laide de la commande Bode, calculer la rponse en frquence du ltre et superposer
son module sur la gure prcdente.
Vrier que le ltrage ralise correctement avec le chier Filtrage.mdl prcdemment
cr. Comparer lordre des ltres de Butterworth et de Chebyshev permettant dobtenir
des actions similaires..
4.3.3 Filtrage avec un ltre RIF :
Le mme problme est repris en utilisant un ltre rponse impulsionnelle nie
Nous pouvons concevoir un ltre RIF par la commande :
B = FIR1(N,Wn)
O B contient les coecients de la squence de pondration du ltre dordre N et de
frquence de coupure 3dB normalise W
n
.
En superposant sur la mme gure la rponse de ce ltre celle des ltres RII prcdents,
dterminer par essais successifs les paramtres N et W
n
qui donnent une rponse peu
prs quivalente.
Implanter le ltre et vrier son fonctionnement. Comparer son ordre ceux des ltre
prcdents.
25
Chapitre 5
BE V : Gnration et dtection des
signaux DTMF : Dual Tone
Multi-Frequency
5.1 Le problme
Le DTMF (Dual Tone Multi-Frequency) est utilis en tlphonie pour le codage du clavier
du tlphone.
Lors de ltablissement dune communication tlphonique, il est ncessaire de transmettre
au standard le ndemand. Sous quelle forme cette indication est-elle transmise ? Dans
les tlphones touche modernes le clavier est vu comme une matrice et lappui sur
une touche envoie un signal qui est le mlange de deux sinusodes. Avec deux groupes de
frquences respectivement hautes et basses nous pouvons coder 16 touches selon le schma
suivant :
Ceci correspond aux signaux temporels et aux spectres associs ports sur la gure sui-
vante.
Dun point de vue pratique deux problmes se posent :
A lmission : comment raliser par programmation les direntes sinusodes ?
A la rception : comment les dtecter ?
Pour cela nous utiliserons la technique du numrique, les signaux utiliss sont des signaux
chantillonns et la norme internationale xe la frquence dchantillonnage F
s
= 8000
Hz.
26
5.2 generation de signaux discrets : ltre generateur
(ou formeur)
5.2.1 Cas gnral :
Tout signal discret peut tre modlis par la rponse impulsionnelle dun ltre (ltre
formeur ou ltre gnrateur) selon le schma : Pour un signal ayant le ltre gnrateur de
fonction de transfert lquation rcurrente entre-sortie scrit :
H(z) =
b
0
+ b
1
z
1
1 + a
1
z
1
+ a
2
z
2
(5.1)
x
k
= b
0
u
k
+ b
1
u
k1
a
1
x
k1
a
2
x
k2
o u(t) est lentre du systme. Ici, u(t) = (t) cest dire que u
k
= 0 pour k = 0.
Pour gnrer le signal x(t) prsum causal, il sut de faire fonctionner cette quation
rcurrente avec des conditions initiales nulles et seul u
0
= 0. Dun point de vue implanta-
tion, cela oblige faire systmatiquement un test sur k pour attribuer la valeur u
k
et u
k1
alors quelle est pratiquement toujours nulle. Pour viter ce test il est possible dcrire le
problme autrement
27
Pour k 0 x
k+2
= a
1
x
k+1
a
2
x
k
Pour dbuter cet algorithme, il faut les deux conditions initiales qui se calculent avec
lquation rcurrente :
28
x
0
= b
0
x
1
= b
1
b
0
.a
1
Tout signal peut tre gnr par une quation rcurrente issue de son ltre formeur
et un choix judicieux des conditions initiales. La dicult est de trouver le ltre formeur
associ un signal donn.
5.2.2 Cas des signaux sinusodaux :
Pour un signal de forme sinusodale le cas est trs simple. Nous savons que :
Quelque soit la phase, sans importance dans le problme du DTMF, un signal sinusodal
est ainsi trs simple gnrer par son ltre formeur.
Algorithme :
Calcul de a
1
= 2cos(2f
0
T
s
)
Calcul des conditions initiales x
0
= b
0
; x
1
= b
1
b
0
a
1
En boucle excuter lquation rcurrente :
x(t) = a
1
x(t 1) x(t 2)
5.2.3 Ralisation : Gnration de signaux sinusodaux :
Programmation de lquation de rcurrence :
Nous utilisons, comme dans les BE prcdents, un chier script de Matlab. On pourra se
rfrer ces BE pour la mise en oeuvre dun tel chier.
En choisissant une frquence dchantillonnage de 8000 Hz pour les frquences f
0
= 440Hz,
523Hz et 1046Hz calculer et tester des signaux de dure 0.25s. Pour cela on eectuera les
oprations suivantes :
29
Calcul des coecients b
0
, b
1
et a
1
pour gnrer un cosinus damplitude 0.4
Calcul des chantillons du cosinus
Calcul des coecients b
0
, b
1
et a
1
pour gnrer un sinus damplitude 0.4
Calcul des chantillons du sinus
Reprsenter les deux signaux sur le mme graphe, les comparer
Avec la fonction sound les couter grce la carte son
Programmation graphique :
En utilisant Simulink implanter lalgorithme de gnration dun sinus. Les blocs "re-
tard" ncessaires sont Signal Processing Blockset > Signal Operations > Delay. Les pa-
ramtres des dirents blocs seront calculs par le chier script prcdent. Visualiser le
signal : (Signal Processing Blockset > Signal Processing Sinks > Time Scope)
et lcouter (Signal Processing Blockset > Signal Processing Sinks > To Wave Device).
En ralit ce travail est dj programm sous forme complte dans le bloc Simulink
: Signal Processing Blockset > Filtering > Filter Designs > Digital Filter.
Ajouter ce bloc et comparer son fonctionnement avec le bloc prcdent.
Signaux du DTMF :
Calcul des paramtres :
Avec un chier script, calculer pour les 8 frquences du DTMF les conditions initiales x
0
;
x
1
et le coecient a
1
. On utilisera pour cela des vecteurs de dimension 8, les frquences
du DTMF tant traditionnellement classes de la plus basse (697Hz index 1) vers la
plus haute (1633Hz index 8). Lamplitude des oscillations obtenues sera xe 0.2.
Ralisation dun bloc lmentaire :
Cet exercice sera prtexte utiliser une fonctionnalit de Simulink : crer des sous-
systmes et leur attribuer un masque de dialogue :
Par copier/coller rcuprer lun des schmas de gnration dune sinusode de ltape
prcdente.
Lui attribuer le coecient a
1
(k) et les conditions initiales x
0
(k) et x
1
(k).
Avec la souris, slectionner lensemble des lments du systme puis crer un sous sys-
tme : menu Edit > Create Subsystem.
Sur le systme nouvellement cr avec la commande Edit > Mask Subsystem nous
pouvons crer un interface graphique avec le sous-programme comme lindique la fentre
suivante :
Dnition des paramtres
Nom de la variable : texte destin rappeler quelle doit tre lentre.
Variable : porter ici le nom utilis dans le sous-systme. Pour nous ce sera k, lindex
de la frquence.
30
Cela sut pour crer un masque simple associ un sous-systme. Une fois ceci
sauvegard, un clic sur le bloc ouvrira une fentre de dialogue invitant prciser
le nde la frquence utilise. Ensuite, le bloc ira rechercher dans lespace de travail
les variables a
1
(k), x
0
(k) et x
1
(k) dsires ce qui implique dexcuter au pralable le
chier script qui les calcule.
5.2.4 Gnration dun signal DTMF :
Dupliquer le bloc ralis, puis sommer la sortie des deux blocs. Initialiser chaque bloc
pour produire une frquence dirente et observer les signaux obtenus dun point de
vue temporel : (Signal Processing Blockset > Signal Processing Sinks > Time
Scope) et frquentiel (Signal Processing Blockset > Signal Processing Sinks >
Spectrum Scope).
PARTIE II
31
5.3 Dtection - algorithme de goertzel
5.3.1 Dtection dune sinusode :
Le problme peu sembler trivial : il sut dun ltre slectif la frquence donne. Cepen-
dant, pour obtenir une slectivit susante, ce ltre va requrir un nombre apprciable
de coecients (par exemple, pour sparer correctement les raies 1336Hz et 1477Hz, il faut
un ltre de Butterworth dordre au moins gal 25, ceci faire fonctionner pour chaque
raie soit 8 fois).
Lautre alternative est de calculer une FFT cependant lalgorithme est assez lourd pro-
grammer et il nous fournit beaucoup trop de renseignements puisquil calcule les N raies
de la bande de frquence de Shannon. Pour dtecter une sinusode, nous navons besoin
que dune raie du spectre pour savoir si oui ou non elle est prsente, cest lide de Goertzel
qui le conduit la lumire de cette rexion une simplication extrme du calcul de la
TFD.
5.3.2 Rappel sur la TFD :
Avec N points, la TFD dun signal est dnie par :
X
n
=
N1

k=0
x
k
e
(j
2
N
)nk
(5.2)
les {xk} tant les N chantillons du signal (k [0; N 1]) et les {X
n
} les N chantillons
de la TFD (n [0; N 1]) qui estiment la TF du signal chantillonn. Lcart frquentiel
entre deux chantillons successifs de cette TFD (prcision frquentielle) est gal :
32

TFD
=
F
s
N
=
1
NT
s
5.3.3 Algorithme de Goertzel :
Nous recherchons une sinusode, nous avons besoin de calculer uniquement la composante
X
n
telle que :
n =
F
0

TFD
=
f
0
f
s
N
Formulation partir de la TFD :
En reprenant les notations de la TFD nous pouvons rcrire la formule de la TFD :
X
n
=
N1

k=0
x
k
e
(j
2
N
)nk
W
N
= e
j
2
N
(5.3)
X
n
=
N1

k=0
x
k
W
nk
N
X
n
=
N1

k=0
x
k
e
(j
2
N
)nk
W
nN
N
= 1 (5.4)
X
n
=
N1

k=0
x
k
W
n(Nk)
N
Solution rcurrente
Posons : Q
m
=
m

k=0
x
k
W
n(m+1k)
N
la solution recherche est : X
n
= Q
N1
.
Nous pouvons crire :
Q
m
=
_
m1

k=0
x
k
W
n(mk)
N
_
W
n
N
+ x
m
W
n
N
Q
m
= Q
m1
W
n
N
+ x
m
W
n
N
Soit la relation de rcurrence :
(1 q
1
W
n
N
)Q
m
= x
m
W
n
N
33
Les Q
m
se dduisent ainsi facilement des x
m
par un ltrage grce une cellule du premier
ordre. Linconvnient est que le ple associ est W
n
N
cest dire quil est complexe. Cela
ncessite dans les calculs des coecients complexes et peu de processeurs permettent
aisment cette manipulation. Une petite manipulation va remdier ce souci.
Un signal intermdiaire - le rsonateur :
En remarquant que : W
n
N
= W
n
N
nous crivons :
(1 q
1
W
n
N
)(1 q
1
W
n
N
)Q
m
= (1 q
1
W
n
N
)x
m
W
n
N
(1 (W
n
N
+ W
n
N
)q
1
+ q
2
)Q
m
= (1 2cos(
2n
N
)q
1
+ q
2
)Q
m
= (1 q
1
W
n
N
)x
m
W
n
N
(5.5)
En faisant intervenir un signal intermdiaire V (t) :
V
m
=
1
1 2cos(
2n
N
)q
1
+ q
2
x
m
Q
m
= (W
n
N
q
1
)V
m
(5.6)
Le signal ltr tant rel, les coecients du ltre tant eux aussi rels le signal V
m
est
un signal rel. Ce signal V (t) est un signal chantillonn obtenu par ltrage de x(t) par
le ltre de fonction de transfert H(z) :
V (z) = H(z)X(z)
H(z) =
z
2
z
2
2cos(
2n
N
)z
1
+ 1
(5.7)
Ce ltre possde deux ples complexes conjugus e
j
2n
N
)
, deux ples sur le cercle unit
do la rsonance pour la frquence f =
n
NT
s
n
N
F
s
= F
0
. Cest la frquence de la sinusode
recherche.
Le rsultat nal :
Le rsultat qui nous intresse est Q
N1
. Le rsonateur fonctionne bien la frquence
F
s
= 8000Hz (T
s
= 125s) alors que le rsultat qui nous intresse ne se calcule que tous
les N.T
s
:
Q
N1
= (W
n
N
q
1
)V
N1
= W
n
N
V
N1
V
N2
Nous nous intressons au module ou module au carr :
34
|Q
N1
|
2
= (W
n
N
V
N1
V
N2
)(W
n
N
V
N1
V
N2
)
= V
2
N1
+ V
2
N2
2cos(
2n
N
)V
N1
V
N2
|X
n
|
2
= V
2
N1
+ V
2
N2
+ a
1
V
N1
V
N2
5.3.4 Implmentation de lalgorithme de Goertzel :
En permanence tourne lalgorithme qui reconstruit le signal V (t). Cest un ltrage du
second ordre des donnes reues avec une quation de rcurrence telle que :
V
m
= x
m
a
1
.V
m1
V
m2
a
1
= 2cos(
2n
N
) (5.8)
Au Nime chantillonnage on peut dterminer lamplitude de la raie recherche soit :
|X
n
|
2
= V
2
N1
+ V
2
N2
a
1
V
N1
V
N2
(5.9)
si la raie nest pas prsente, cette valeur sera nulle (en pratique faible). Pour les chantillons
suivants, on rinitialise lalgorithme pour reprendre le calcul tant quil y a un signal :
lalgorithme est rinitialis tous les N.T
s
.
En pratique, on peut aussi faire fonctionner lensemble la priode dchantillonnage T
s
,
le rsultat du calcul ntant pris en compte que tous les instants N.T
s
.
35
Norme internationale :
Dans la pratique la norme utilise est celle du Bell Laboratory soit : F
s
= 8000 Hz ;
N = 205 (
TFD
=
F
s
N
39Hz ce qui est susant pour sparer les raies du codage
DTMF. Le tableau des valeurs de n utiliser est le suivant (la valeur pratique est lentier
immdiatement le plus proche) :
f
0
(en Hz) 697 770 852 941 1209 1336 1477 1633
n 18 20 22 24 31 34 38 42
Le processeur ralisant lalgorithme dispose dau plus 125 s pour cela, pour dcider de la
prsence ou non dune raie il faudra au moins 26ms. Les signaux DTMF durant environ
100 ms, le calcul est ralis 4 5 fois lors de la rception dun signal.
5.3.5 Ralisation : Dtection des frquences DTMF :
Programmation de lquation de rcurrence :
Avec un nouveau chier script on implmentera lalgorithme de Goertzel selon les grandes
lignes suivantes :
F
s
= 8000 Hz, N = 205
Calculer les valeurs de n du tableau et en dduire le vecteur a
1
(k) des valeurs utiliser
dans le rsonateur (k nde la frquence DTMF).
Crer un signal superposition de deux sinusodes du DTMF
Implmenter lalgorithme de calcul de Goertzel (Remarque : pour une programmation
simplie, on pourra utiliser un tableau V
m
(3, 8) 1 ligne par rsonateur, la colonne
correspondant aux 3 valeurs ncessaires la rcursivit. Il peut tre intressant dutiliser
les produits points).
Acher graphiquement le rsultat de la simulation. Changer les frquences du signal
pour vrier le fonctionnement de lalgorithme.
Ralisation avec Simulink :
On utilisera les a
1
(k) calculs par le chier prcdent.
Un rsonateur :
Avec Simulink, raliser lalgorithme de calcul. Pour se simplier lexistence, on fera tra-
vailler la mme frquence le rsonateur et le calcul de |X(n)|
2
. Il est important de
rinitialiser lalgorithme tous les N coups et pour cela, il faut utiliser les retards du rso-
nateur avec loption Reset, cette entre sera rinitialise par un gnrateur dimpulsions
de cadence N.T
s
(Simulink > Sources > Pulse Generator).
Pour tester son fonctionnement, on utilisera un gnrateur BF appliqu sur lentre ligne
de lordinateur. Faire varier la frquence et observer la sortie.
36
Dtecteur complet :
Dupliquer 7 fois le bloc, attribuer chacun une frquence du DTMF (il peut tre
astucieux de les mettre dans lordre).
Utiliser comme signal dentre celui fournit par le chier NumeroBE.wav
Visualiser sur un mme scope le signal dentre et les 8 sorties des rsonateurs
En dduire le signal le nde tlphone enregistr dans le chier .wav
37
Chapitre 6
BE VI : Filtres rponse
impulsionnelle innie RII :
Transposition Continu-Discret
Les mthodes de transposition sont des mthode de synthse de ltres RII : elles partent du
principe que le problme de la synthse des ltres a dj t largement dvelopp dans le
domaine du signal continu. Il a t tabli plusieurs mthodes dapproximation polynmiale
dun gabarit aboutissant aux ltres de Butterworth, Tchebychev, Cauer, Legendre.... Le
ltrage des signaux discrets reprend ces mthodes en cherchant simplement une technique
de passage du continu au discret qui prserverai au mieux les caractristiques du ltre.
Principe de la transposition :
Nous supposons possder un ltre continu dont la fonction de transfert isochrone H(jw)
constitue une approximation satisfaisante du gabarit du ltre recherch. Il nous reste
chercher le ltre discret de transfert isochrone G(e
jwT
s
) quivalent. H(p) G(z).
Il nous faut une mthode laquelle nous imposons les contraintes suivantes :
- Etre simple.
- Faire en sorte que G(e
jwT
s
) soit la plus proche possible de H(jw).
- Obtenir une fonction de transfert discrte G(z) qui puisse simplanter sous forme dqua-
tion rcurrente ce qui ncessite une forme de rapport de polynmes B(z)/A(z).
6.1 Filtre rjecteur
Un ltre continu rjecteur de bande a comme fonction de transfert :
G(p) =
p
2
+ 1
(p + 1)
2
(6.1)
38
Fig. 6.1 rejet
Le module de sa rponse en frquence est indiqu sur la gure ??. Ce ltre travaille en
frquences basses et sa ralisation est dlicate. Nous allons rechercher un ltre discret
ralisant la mme fonction grce une mthode de transposition.
1. La priode dchantillonnage du ltre discret est : T
s
= 1s.
Pour synthtiser le ltre numrique correspondant, nous avons utilis la transposi-
tion bilinaire et nous avons obtenu la fonction suivante :
H
1
(z) =
5z
2
6z + 5
9z
2
6z + 1
(6.2)
- Quels sont les points singuliers responsables de leet recherch (eet rjecteur de
bande). Les reprsenter dans le plan z. On calculera le module de la rponse en
frquence pour = 0, =
max
(maximum autoris par le thorme de Shannon)
et =
R
(pulsation rejete).
- Comparer aux caractristiques du ltre continu. Comment sinterprtent ces rsul-
tats ?
2. Transformation avec pr-dcalage (Transformation de Tustin) : Le rsultat nest pas
compltement satisfaisant. Eectuer les oprations ncessaires pour quil le devienne.
Tracer le module de la FFT de la nouvelle fonction de transfert.
39
10
2
10
1
10
0
10
1
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
frquence
M
o
d
u
l
e

d
e

l
a

r

p
o
n
s
e

e
n

f
r

q
u
e
n
c
e
Fig. 6.2 Gain de la transmittance en frquence
6.2 Etude dun ltre du 2me ordre
6.2.1 Fonction de transfert dun ltre continu
La fonction de transfert en p dun ltre analogique de Butterworth de second dordre pour
une pulsation de coupure
c
= 1, est :
H
2
(p) =
1
1 +

2p + p
2
(6.3)
On veut la fontion de transfert en frquence H
2
(f). On a la relation p = j2f, do :
H
2
(f) =
1
1 + j2

2f (2f)
2
(6.4)
On value le module | H
2
(f) | :
| H
2
(f) |=
1
_
1 + (2f)
4
(6.5)
La frquence de coupure f
c
est obtenue 3dB quand on est au maximum du gain divis
par

2. Le gain maximum est de 1, donc :
| H
2
(f
c
) |=
1
_
1 + (2f
c
)
4
=
1

2
(6.6)
On en dduit f
c
=
1
2
.
40
6.2.2 Fonction de transfert dun ltre discret
Transposition bilinaire :
On veut maintenant dterminer G
2
(z) la fonction de transfert en z dun ltre numrique
de pas temporel T
s
ralisant le mme ltrage que prcdemment.
On utilise la transformation bilinaire pour passer de la variable p la variable z :
p =
2
T
s
.
1 z
1
1 + z
1
(6.7)
Lquation de rcurrence temporelle reliant le signal de sortie y
n
au signal numrique
dentre x
n
scrit sous la forme :
y
n
=
1
1 +

2 +
2
[x
n
+ 2x
n1
+ x
n2

2(1
2
)y
n1
(1

2 +
2
)y
n2
]
avec =
2
T
s
.
y
n
=
1
1 +
2
T
s

2 +
_
2
T
s
_
2
[x
n
+ 2x
n1
+ x
n2

2(1 +
_
2
T
s
_
2
)y
n1
(1
2
T
s

2 +
_
2
T
s
_
2
)y
n2
]
G
2
(z) =
1 + 2z
1
+ z
2
1 +

2 +
2
+ 2(1
2
)z
n1
+ (1

2 +
2
)z
n2
On veut la rponse impulsionnelle h(t). En signal dentre on doit donc avoir une impulsion
de Dirac. Pour reprsenter une impulsion de Dirac sous matlab, on cre le vecteur suivant :
> x=zeros(1,10000) ;
> x(1)=1 ;
Et pour dterminer la rponse impulsionnelle, daprs la dfinition de la
fonction "filter" :
>B=1/(1+alpha*sqrt(2)+alpha
2
)*[1 2 1] ;
>A=1/(1+alpha*sqrt(2)+alpha
2
)*[ (1+alpha*sqrt(2)+alpha
2
)
2*(1-alpha
2
) (1-alpha*sqrt(2)+alpha
2
] ;
41
>y=filter(B,A,x)
- Pour passer de la fonction de transfert en z du ltre numrique une transmittance en
frquence, on pose : z = e
j2T
s
f
.
En dduire la rponse en frquence du ltre numrique.
- Comparer les modules des ltres analogiques et numriques | H
2
(f) | et | G
2
(f) |. On
prend T
s
= 0.01s.
- Calculer lerreur quadratique entre le modle analogique et le modle numrique.
- Conclusion ?
Invariance impulsionnelle
Une autre ide pour raliser la transposition continu-discret est de se rappeler quun
systme est caractris par sa rponse impulsionnelle. Nous pouvons donc prendre comme
systme discret celui dont la rponse impulsionnelle est la version chantillonne de celle du
systme continu dorigine. Dun point de vue mathmatique lopration est donc simple :
pour un systme continu de rponse impulsionnelle h(t) et de fonction de transfert H(p)
on obtiendra un systme discret de fonction de transfert G(z) correspondant une rponse
impulsionnelle g(t).
1. Dterminer et tracer la rponse impulsionelle du ltre numrique. Ce ltre est t-il
stable ? Justier votre rponse.
2. Retrouver la rponse en frquence du ltre numrique partir de sa rponse impul-
sionnelle. Comparer le rsultat obtenu avec celui obtenu prcdement (frquence de
coupure, pente).
3. On souhaite dsormais avoir une frquence de coupure f
c
= f
0
quelconque. On doit
alors transformer lquation de rcurrence temporelle en changeant par :

0
=

2f
0
4. Crer un programme Matlab qui pour direntes valeurs de f
0
calcule la rponse
impulsionnelle et en dduire par Transforme de Fourier la courbe du gain en fonction
de la frquence.
Tracer les courbes de gains obtenues our direntes f
0
, on choisit :
f
0
= [0.01Hz 0.1Hz 1Hz 10Hz].
Que constatez vous ?
5. Etudier les sorties de ce ltre numrique pour dirents signaux basses et hautes
frquences.
soit le signal dentre x(t) gnrant des hautes et basses frquences :
42
x(t) =
_
sin
_
2
T
1
t
_
+ sin
_
2
T
2
t
__
e
t
T
s
(6.8)
o T
1
et T
2
sont des contantes de temps et T
s
la priode dchantillonnage. 1/T
1
et 1/T
2
reprsentent des frquences. On veut des basses et hautes frquences. On
choisit par exemple :
T
1
= 0.1s : pour les hautes frquences.
T
2
= 0.0001s : pour les basses frquences.
6. Tracer les rponses y(t) au signal x(t) pour les valeurs suivantes :
1. f
0
= 0.1Hz, T
1
= 0.1s, T
2
= 0.01s
2. f
0
= 0.1Hz, T
1
= 10s, T
2
= 0.01s
3. f
0
= 0.01Hz, T
1
= 0.1s, T
2
= 0.01s
7. Tracer les gains pour direntes valeurs de T
1
et T
2
.
T
1
T
2
1 s 0.00001 s
10 s 0.00001 s
0.1 s 0.0001 s
0.05 s 0.0001 s
0.03 s 0.0001 s
0.025 s 0.0001 s
Approximation de la drive
Faire la transposition continu-dscret en utilisant cette fois-ci la mthode de lapproxima-
tion de la drive :
p
z 1
T
s
Etudier la qualit de cette transposition.
6.2.3 Etude dun ltre du 3me ordre
Soit le ltre numrique dquation de rcurrence temporelle suivante :
y
n
=
1
3
y
n2
+
1
6
x
n

1
2
x
n1
+
1
2
x
n2

1
6
x
n3
(6.9)
1. Tracer la rponse impulsionnelle et la rponse en frquence du gain de ce ltre.
conclusions. De quel type de ltre sagit t-il ?
43
2. On tudie maintenant leet du ltre dun signal sur un signal contenant les basses
et les hautes frquences. On reprend le mme signal x(t) que prcdement.
On trace direntes courbes pour direntes valeurs de T
1
et T
2
. Conclusions.
44
Chapitre 7
BE VII : Dbruitage dun
lectrocardiogramme (ECG) par ltrage
7.1 Analyse dun enregistrement dun ECG :
A laide de la commande Matlab load(nom du fichier) chargez le chier "data.mat"
qui se trouve dans lespace de travail (work). Ce chier contient la variable "datas" qui
correspond aux chantillons de 20 secondes denregistrement dun ECG (lectrocardio-
gramme). Malheureusement, ce signal a t perturb par un bruit haute frquence.
On vous demande de ltrer ce signal an disoler le signal utile. Notons que laspect visuel
tant important dans ce type de signal biomdical, on veillera choisir le ltre de faon
ce que la forme du signal soit prserve.
Dans un premier temps, an danlayser ce signal ECG :
1. Spcier la frquence laquelle le signal a t chantillonn.
2. Acher le spectre de ce signal et dindiquer les commandes utilises.
3. Noter les caractristiques de votre ltre ainsi que le raisonnement qui a t men et
les commandes Matlab utilises.
4. Donner la valeur de lamplitude A et la frquence de la composante fondamentale
(A cos(wt)) du signal ltr partir de lachage du contenu frquentiel et mettez
en relation la frquence fondamentale trouve avec lvolution temporelle du signal
ltr.
7.2 Filtrage du bruit :
La gure ci-dessous reprsente lachage de la rponse impulsionnelle dun ltre.
1. Donnez lexpression de la fonction de transfert oprationnelle du ltre ?
45
2. A laide de la fonction Matlab permettant dacher les ples et les zros dune
fonction de transfert, prdisez lallure de la fonction de transfert du ltre. Expliquez
votre raisonnement et indiquez les lignes de commande Matlab utilises.
3. Vriez que vos prdictions sont correctes (indiquez les lignes de commande Matlab
utilises).
4. Ajustez le gain du ltre an quil ne modie pas le continu.
5. Donnez les caractristiques de ce ltre (type, bande passante 3dB, bande blo-
quante).
6. Est-ce que ce ltre conviendrait pour analyser un signal ECG pollu par des hautes
frquences qui a t chantillonn 1000Hz et dont la gamme de frquences utiles
stendrait jusque 50Hz ?
46