Traitement Numérique Du Signal - 2020
Traitement Numérique Du Signal - 2020
Traitement Numérique Du Signal - 2020
net/publication/280805015
CITATIONS READS
0 9,641
1 author:
Assia Kourgli
University of Science and Technology Houari Boumediene
43 PUBLICATIONS 153 CITATIONS
SEE PROFILE
Some of the authors of this publication are also working on these related projects:
All content following this page was uploaded by Assia Kourgli on 15 June 2020.
- Fenêtres de pondération
- Exercices
- Transformées en Z
- Propriétés de la TZ
- TZ rationnelles
- Caractéristiques des FN
- Exercices
- Exercices
- Synthèse des filtres RII par la méthode des pôles et des zéros
- Exercices
- Exercices
- Analyse multi-résolution
- Exemples d’applications
- Exercices
Travaux Pratiques
5. Filtrage Multi-Cadences
Un signal est la représentation physique de l’information qu’il transporte de sa source à son destinataire. Il
sert de vecteur à une information. Il constitue la manifestation physique d’une grandeur mesurable (courant,
tension, force, température, pression, etc.). Les signaux sont des grandeurs électriques variant en fonction du
temps x(t) obtenues à l’aide de capteurs. Sur le plan analytique : Un signal sera une fonction d'une variable réelle,
en général le temps.
Amplitude
Onde acoustique : délivré par un microphone (parole, musique, …) 0
-0.5
Signaux biologiques : EEG, ECG
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
Temps (s)
Tension aux bornes composant électronique
Voyelle sur 25 ms
0.5
Signaux géophysiques : vibrations sismiques
Amplitude
0
Images, Vidéos 0.954 0.956 0.958 0.96 0.962 0.964 0.966 0.968 0.97 0.972 0.974
Temps (s)
Remarque: Tout signal physique comporte une composante aléatoire (perturbation externe, bruit, erreur de mesure,
etc …). Le bruit est défini comme tout phénomène perturbateur gênant la perception ou l’interprétation d’un
signal, par analogie avec les nuisances acoustiques (interférence, bruit de fond, etc.). La différentiation entre le
signal et le bruit est artificielle et dépend de l’intérêt de l’utilisateur : les ondes électromagnétiques d’origine
galactique sont du bruit pour un ingénieur des télécommunications par satellites et un signal utile pour les
radioastronomes.
Les fonctions du traitement du signal peuvent se diviser en deux catégories : l’élaboration des signaux
(incorporation des informations)et l’interprétation des signaux (extraction des informations). Les principales
fonctions intégrées dans ces deux parties sont les suivantes [1]:
Élaboration des signaux : synthèse, modulation, codage/compression, etc.
Interprétation des signaux : filtrage, détection, identification, analyse, mesure, etc.
La numérisation d’un signal est l’opération qui consiste à faire passer un signal de la représentation dans le
domaine des temps et des amplitudes continus au domaine des temps et des amplitudes discrets. Cette opération
de numérisation d’un signal peut être décomposée en deux étapes principales : échantillonnage et quantification.
La restitution (ou l’interpolation) constitue le processus inverse qui intervient lors du passage du signal
numérique au signal analogique : commande d’un actionneur. Ces trois étapes sont indissociables. En effet, le
signal, étant le support physique d’une information, doit conserver au cours de ces modifications tout le contenu
informatif initial. Cette condition, ajoutée à la notion de coût limite d’un système, va être à la base de la
numérisation des signaux et de l’étude du traitement numérique [1].
Les traitements numériques sont aisément réalisés grâce à des additionneurs, des multiplieurs numériques, et
des mémoires. En outre, les systèmes numériques possèdent de nombreux avantages comparées à ceux
analogiques, entre autres [2]:
- Simplicité: Les systèmes numériques sont intrinsèquement plus simples à analyser (et donc à synthétiser) que
les systèmes analogiques
- Possibilités de traitement accrues: Il est possible de réaliser, en numérique, des opérations beaucoup plus
complexes qu’en analogique, notamment des opérations non-linéaires.
- Robustesse aux bruits. Les systèmes numériques sont par essence insensibles aux bruits parasites
électromagnétiques. Le transcodage de l’information sous forme numérique joue un peu le rôle de « firewall ».
- Précision et stabilité. Puisque les seuls « bruits » sont liés à la précision des calculs, cette dernière dépend
uniquement du calculateur utilisé ; elle est insensible à la température et ne varie pas avec l’âge du système.
- Flexibilité. Dans un grand nombre de systèmes numériques, le traitement est défini par un logiciel chargé en
mémoire. Il est dès lors très facile de modifier ce traitement, sans devoir modifier la machine qui le réalise.
Rappelons que les signaux déterministes renferment une information dont l'évolution en fonction du temps
peut être parfaitement prédite par un modèle mathématique (au contraire des signaux aléatoires/stochastiques).
Nous présentons dans cette section quelques fonctions mathématiques supports de signaux élémentaires et
utilisées tout au long du cours de traitement du signal numérique. Nous y introduisons, également, quelques
outils permettant de caractériser et d'établir une analyse temporelle des signaux numériques.
Rappelons qu'un signal à temps discret provient souvent de l’échantillonnage à la cadence fe = 1/Te, d’un
signal x(t) déterministe à temps continu qui est supposé à bande limitée (−fe/2, fe/2). Nous noterons x(n) = x(nTe).
Fonction Signe
1
0.8
0.6
0.4
0.2
- Fonction signe 0
-0.2
1 n≥0
-0.4
sgn(n) = -0.6
− 1 n<0 -0.8
-1
-10 -8 -6 -4 -2 0 2 4 6 8 10
Fonction Echellon
1
0.9
- Fonction échelon (unité) 0.8
0.7
1 n≥0 0.6
0 n<0 0.4
0.3
0.2
0.1
0
-10 -8 -6 -4 -2 0 2 4 6 8 10
Fonction Porte
1
0.9
0.7
0.6
1 − N /2 ≤ n ≤ N /2 0.5
Π ( n) = 0.4
N +1
0 ailleurs 0.3
0.2
0.1
0
-10 -8 -6 -4 -2 0 2 4 6 8 10
Fonction Rectangle
1
0.9
0.8
0.6
1 0 ≤ n ≤ N −1
0.5
rect(n / N ) = 0.4
0 ailleurs 0.3
0.2
0.1
0
-10 -8 -6 -4 -2 0 2 4 6 8 10
0.9
1 n=0 0.8
δ ( n) = = U (n) − U (n − 1) 0.7
0 n≠0 0.6
0.5
0.4
0.2
0.1
0
-10 -8 -6 -4 -2 0 2 4 6 8 10
Fonction Sinc
1
- Fonction sinus cardinal
0.8
sin(πθ n)
sin c(θ n) =
0.6
πθ n 0.4
0.2
-0.2
-0.4
-50 -40 -30 -20 -10 0 10 20 30 40 50
Energie et puissance
Toute transmission d’information s’accompagne de transferts d’énergie. En effet, les signaux continus ou
discrets sont essentiellement caractérisés par l’énergie ou la puissance qu’ils véhiculent. Ce sont les seules
grandeurs physiques auxquelles sont sensibles les détecteurs. Beaucoup de capteurs physiques mesurent une
énergie ou une quantité quadratique. Par exemple, les capteurs optiques mesurent une intensité, les compteurs
d’électricité mesurent une énergie, etc. Compte tenu de la définition fondamentale, l’énergie du signal entre les
instants t et t+dt est : |x(t)|2 dt (puissance instantanée multipliée par le temps).
+∞
2
Soit un signal x(n) à temps discret, tel que x(n)existe et converge. Alors le signal est dit à énergie finie
+∞
−∞
Ex =
2
et la valeur de cette somme est appelée énergie du signal : x(n)
−∞
Exemples:
x(n) = Rect(n/N) énergie finie. x(n) = a (constante) et x(n) = A sin(2πf0n) ne sont pas à énergie finie
Pour un signal périodique, cette somme ne converge pas. On peut néanmoins définir la puissance d'un
signal x(n) périodique de période N par :
1 N / 2−1 N −1
Px = 1
2
Px = x ( n)
2
x(n) ou
N −N / 2 2. N −N
Dans le cas général, on parle de signaux à puissance moyenne finie définie par:
1 N / 2−1 N −1
Px = limN →∞ 1
2
Px = lim N →∞ x( n)
2
x(n) ou
N −N / 2 2. N −N
Exemples:
Il existe des signaux ni périodiques, ni d'énergie finie, pour lesquels la puissance ne peut être définie, comme par
exemple la rampe x(n)=n. Il s'agit là de définitions mathématiques, en pratique, un signal mesuré ne l'est jamais
sur un intervalle de temps infini. On peut commencer à visualiser un signal à un instant qu'on prendra comme
origine des temps, et dans ce cas on arrêtera son examen au bout d'un temps Tobs :
Nobs
Ex = x (n )
2
n =0
Remarques
Signal à énergie finie puissance nulle Signal à puissance finie énergie infinie
Le calcul de l'énergie ou la puissance permet d'obtenir une première caractérisation du signal. Par ailleurs,
la théorie du signal a largement développé des méthodes d’étude basées sur la corrélation pour caractériser le
comportement temporel du signal.
Corrélation et auto-corrélation
La fonction de corrélation permet de mesurer le degré de ressemblance entre deux signaux en fonction d’un
décalage. Considérons x(n) et y (n) deux signaux d'énergie finie, la fonction d'intercorrélation Rx,y(k) est définie
∞
par: Rxy (k ) = x(n) y (n − k )
n = −∞
*
1 N / 2 −1
Pour des signaux à puissance moyenne finie, elle vaut : Rxy (k ) = lim x(n) y* (n − k )
N →∞ N n = − N / 2
Exemples:Soient un signal aléatoire et sa version décalée de 50s. On remarque que les signaux se ressemblent le
plus quand y(n) est décalé de 50 secondes.
Intercorrélation maximum à 50 s
3 3 200
x(n) y(n)
2
2
150
1
1
100
0
0
-1
50
-1
-2
0
-3 -2
-4
0 50 100 150 200 250 -3 -50
0 50 100 150 200 250 -250 -200 -150 -100 -50 0 50 100 150 200 250
FEI, USTHB [assiakourgli@gmail.com / http://perso.usthb.dz/~akourgli/
8
Traitement Numérique du Signal Master ST M1 2019/2020
Pour l'auto-corrélation, on remplace y(n) par x(n) on obtient l'expression de l'auto-corrélation pour les
∞
signaux à énergie finie: Rx (k ) = x(n)x (n − k)
n=−∞
*
L'auto-corrélation permet de détecter des régularités, des profils répétés dans un signal comme un signal
périodique perturbé par beaucoup de bruit (Voir TP n°0)
Propriétés :
Auto-corrélation des signaux périodiques : Le calcul sur une seule période suffit. L’auto-corrélation d’un
signal périodique est-elle même périodique. Par définition, le signal périodique ressemble parfaitement à lui-
même, décalé d’une ou plusieurs périodes.
x(n) périodique Autocorrélation de x(n)
1 1
- signaux périodiques
1 N −n
Rx (k ) = x(n) x* (n − k )
0.8 0.8
0.4 0.4
0.2 0.2
0 0
0 0.05 0.1 -0.1 -0.05 0 0.05 0.1
Exemples d'application : Extraction d'un signal noyé dans du bruit, mesure d'un temps ou retard, détection d'un
signal périodique(Voir TP n° 0).L'exemple ci-dessous illustre l'auto-corrélation d'un signal sinusoïdal d'amplitude
1 noyé dans du bruit Gaussien de variance 1.
4 200
sinusoide bruitée Autocorrélation de x(n)
3
150
2
100
1
50
0
-1
-50
-2
-3 -100
0 0.002 0.004 0.006 0.008 0.01 0.012 0.014 -0.015 -0.01 -0.005 0 0.005 0.01 0.015
La corrélation est largement utilisée dans les systèmes radar. Ainsi, pour détecter un avion, on envoie une
impulsion, puis on reçoit une version retardée, atténuée et bruitée de cette impulsion . L'intercorrélation su signal
reçu et émis présentera un pic à l'instant correspondant au retard.
1.5 10
signal émis signal reçu: bruité et retardé de 50S Intercorrélation du signal reçu et émis avec pic en 50s
1 8
1
6
0.8 0.5
4
0.6 0
0.4 -0.5
0
0.2 -1
-2
0 -1.5
0 50 100 150 200 250 300 350 400 450 500 0 50 100 150 200 250 300 350 400 450 500 -4
0 50 100 150 200 250 300 350 400 450 500
Remarques : La notion de bruit est relative, elle dépend du contexte. Le rapport signal/bruit désigne la qualité de
la transmission d'une information par rapport aux parasites. Il est défini par:
3
1.5
SNR=0 db SNR=20 db
2 1
SNRdb=20 Log(PS/PB)
1 0.5
0 0
-1 -0.5
-2 -1
-3 -1.5
0 0.002 0.004 0.006 0.008 0.01 0.012 0.014 0 0.002 0.004 0.006 0.008 0.01 0.012 0.014
Un système linéaire est un modèle de système qui applique un opérateur linéaire à un signal d'entrée.
C'est une abstraction mathématique très utile en automatique, traitement du signal, mécanique et
télécommunications. Les systèmes linéaires sont ainsi fréquemment utilisés pour décrire un système non linéaire
en ignorant les petites non-linéarités. Un système est discret, si à la suite d'entrée discrète x(n) correspond une
suite de sortie discrète y(n).
x(n) y(n)
Système discret
- Si l'entrée x(n) produit une sortie y(n), quand on applique une entrée k.x(n) , la sortie sera k.y(n). Si deux entrées
x1(n) et x2(n) engendrent deux sorties y1(n) et y2(n) alors x1(n) + x2(n) engendrera y1(n) + y2(n)
Exemple
y1 réponse à x1 y2 réponse à x2 y1+y2 réponse à x1+x2
20 20 20
18 18 18
16 16 16
14 14 14
12 12 12
10 10 10
8 8 8
6 6 6
4 4 4
2 2 2
0 0 0
0 10 20 30 40 50 0 10 20 30 40 50 0 10 20 30 40 50
- S’il y a invariance dans le temps, une translation de l'entrée (x(n)x(n-m)) se traduira par une même translation
dans le temps de la sortie (y(n)y(n-m)).
x(n-m) y(n-m)
Système discret
Si le système est invariant, cela implique que le système réagit de la même façon quel que soit l’instant auquel
nous appliquons ses excitations. Cette propriété exprime que la caractéristique du système ne dépend pas de
l’origine du temps, on parle encore de stationnarité.
Convolution
Si les hypothèses de linéarité et d'invariance temporelle sont vérifiées, on peut caractériser le système par sa
réponse impulsionnelle h(n).
On peut en déduire l'effet d'une entrée quelconque sous la forme d'une convolution. Cettedernière est
l’opération de traitement de signal la plus fondamentale. Elle indique que la valeur du signal de sortie à l’instant
n est obtenue par la sommation (intégrale) pondérée des valeurs passées du signal d'excitation x(n). La fonction
de pondération est précisément la réponse impulsionnelle h(n):
∞ ∞
y ( n) = h( n) ∗ x( n) = h( m) x ( n − m) =
m = −∞
x (m) h(n − m)
m = −∞
La réponse impulsionnelle h(n) est le signal qu'on obtient en sortie y(n)=h(n) si on applique en entrée une
impulsion "de Dirac'' x(n)=δ(n). Le Dirac est l'élément neutre de l'opération de convolution:
Exemple 1:
On distingue 3 cas :
1
x(n)
0.8
0.6
0.4
Exemple 2 : 0.2
0
N /2 0 50 100 150 200 250 300
1 1
h(n) = Π (n) y (n) =
N + 1 N +1
x ( n + m)
N + 1 m=− N / 2 3
y(n)
2
(Voir TP n°0)
1
0
0 50 100 150 200 250 300
Calculer y(n)=x(n)*h(n)
Y(n)={2, 3, 5, 10, 3, 9}
Remarques:
- Les opérations de corrélation et convolution sont liées. Mathématiquement, on peut écrire une relation qui
permet d’exprimer la fonction de corrélation comme un produit de convolution (et réciproquement).
On peut donc considérer l'opération d'un SLI comme une mesure de la corrélation entre deux signaux (x*(-n) et
h(n)). En fait, le signal de sortie est "construit" à partir des composantes fréquentielles communes au signal
d'entrée et à la réponse impulsionnelle.
- Si on applique à un SLIT une entrée sinusoïdale réelle ou complexe de fréquence f0, alors, la sortie sera une
sinusoïde dont l'amplitude et la phase pourront être modifiées mais qui conservera la même forme (une
sinusoïde) et la même fréquence f0. On dit que les sinusoïdes sont les fonctions propres des SLIT.
- Un système linéaire invariant est un système dont le comportement dans le temps, peut-être décrit par une
M N
équation aux différences : ai y(n − i) = bi x(n − i) ,
i =0 i =0
Une contrainte importante pour la formalisation de nombreux problèmes est de respecter la notion de
causalité (les effets ne peuvent pas précéder la cause). Dans le cas des SLIT, cette causalité se traduit par le fait que
pour: h(n) = 0 pour n<0.
n +∞
x(n) = 0, n < n0 alors y(n) = 0, n < n0 h( n) = 0, n < 0 , y (n) = x (m) h( n − m) = x ( n − m) h( m) ,
m = −∞ m =0
n
- si h et x sont causaux y ( n ) = h(n − m) x(m)
m =0
Remarque:
Nous pouvons envisager mémoriser les signaux d’entrée et faire un traitement de ceux-ci en temps différé, les
systèmes utilisés ne sont plus alors nécessairement causaux car pour élaborer la sortie à l’instant ni, nous
disposons en mémoire des entrées aux instants suivants. C’est souvent le cas en traitement d’image, en traitement
de parole effectué après mémorisation du signal à traiter.
Une autre notion fondamentale est la stabilité des systèmes. La propriété de stabilité des systèmes bouclés
est non seulement une performance mais une exigence pour le bon fonctionnement d’une boucle d’asservissement
FEI, USTHB [assiakourgli@gmail.com / http://perso.usthb.dz/~akourgli/
13
Traitement Numérique du Signal Master ST M1 2019/2020
ou de régulation. Une boucle instable est une boucle inutilisable. La définition la plus courante de cette stabilité
est la suivante : 7
Système stable
10
Système instable
9
6.5
8
6
7
5.5
6
h(n) 4.5
4
3
4
2
3.5
1
3 0
0 5 10 15 20 25 30 35 40 45 50 0 5 10 15 20 25 30 35 40 45 50
On dit qu'un système est stable si, en lui appliquant une entrée bornée quelconque, la sortie reste bornée,
ce qui implique dans le cas des SLIT:
h( n) < ∞
n
- Si les ai sont ≠ de 0, le système est dit récursif (RII), il est non récursif s'il ne dépend que des x(n-i) (RIF)
K
- Si le système est à réponse impulsionnelle de durée finie (RIF), alors : y ( n) = h ( m) x ( n − m)
m=0
Dans ce cas, le système numérique est une fenêtre centrée sur les K plus récents échantillons.
+∞
- Si le système est à réponse impulsionnelle de durée infinie (RII) : y ( n ) = h ( m) x ( n − m)
m=0
Dans ce cas, il est nécessaire de connaître tous les échantillons présents et passés, le système à une mémoire de
longueur infinie.
avec comme réponse impulsionnelle h(n)=δ(n)+a1δ (n-1)+a2δ (n-2)+.......+akδ (n-k) qui, on peut le constater,est bien
finie.
Exemple 2y(n)=x(n)+a1y(n-1) est l'équation aux différence finies d'un filtre RII
y(n)=x(n)+a1x(n-1)+a12x(n-2))+a13x(n-3)+a14x(n-3)+.....+ a1my(n-m)
En poursuivant le procédé à l'infini y(n) dépend d'une infinité de x(n-k) ce qui en fait un filtre RII.
Exemple d'application
Les séquences x(n) (réel) et y(n) représentent respectivement l’entrée et la sortie d’un système discret. Pour
chaque cas, identifiez celles représentant
a) des systèmes linéaires, b) des systèmes causals, c) des systèmes invariants aux translations de n,
d) des systèmes assurément ou possiblement stables (en fonctions des constantes)
1. y(n) = x(n) + bx(n-1) 2. y(n) = x(n) + bx(n+1) 6. y(n) = bx(n) b : constante réelle
3. y(n) = nx(n) 5. y(n) = x(n)e n 7. y(n) = |x(n)|
4. y(n) = x(n) sin(2πf0n) N : constante entière
a) Tous sauf 6 et 7 b) Tous sauf 2
c) Les systèmes 1, 2, 6 et 7 d) 1 (b finie), 2 (b finie), 4, 6 (b finie) et 7.
Si X(f) = fonction réel ⇔x(t) est paire Si X(f) = fonction imaginaire pure ⇔x(t) est impaire
Remarques :
X ( f ) = A ( f ) + B ( f ) et ϕ = arg( X ( f )) = arctg
2 2 B( f )
A( f )
La TF d'un signal périodique est divergente, mais on peut définir une TF au sens des distributions en utilisant
la décomposition en Série de Fourier. Le résultat correspond à un spectre de raies (non continu):
T /2 . +∞
X ( f ) = x (t )e − 2π j f t dt
1
sachant que: Cn =
T x(t ) exp( −2πjnf 0 t )dt
−T / 2
et
−∞
alors X ( f ) = lim( T .C n )
T → +∞
+∞s +∞s
Pour les signaux à énergie finie, la TF conserve l'énergie (relation de Parseval) : E =
x(t ) dt = X( f )
2 2
x df
−∞ −∞
On peut donc définir une notion d'énergie par unité de fréquence, la densité spectrale d'énergie (DSE). La DES
est la TF de l'autocorrélation (Thèorème de Wiener-Kintchine)
+∞s
R (τ )e − 2 π j f τ dτ
2
Sx ( f ) = X ( f ) = x
−∞
Pour les signaux à puissance moyenne finie, on définit alors une densité spectrale de puissance (DSP):
2
X(f )
Px ( f ) = lim
T →∞ T
- La propriété de changement d'échelle indique que plus le support temporel d'une fonction est étroit plus le
support de sa TF est large.
FEI, USTHB [assiakourgli@gmail.com / http://perso.usthb.dz/~akourgli/
15
Traitement Numérique du Signal Master ST M1 2019/2020
- La translation d'un signal temporel se traduit par un déphasage en fréquence. Une translation en fréquence
équivaut à une modulation temporelle.
- La propriété de dualité permet d’obtenir facilement de nouvelles paires de transformées de FOURIER à partir
des paires déjà connues.
Exemple d'Application
Soient x(t)= Α Π (t ) et y(t) = Α Λ (t ) signal porte et signal triangulaire d’amplitude A>0 et de largeur θ.
θ θ /2
• Tracer x(t) et y(t) sur le même graphe
• Utiliser les dérivées pour trouver X(f) et Y(f) qui seront représentés sur le même graphe puis commenter et
interpréter les deux graphes pour A=1 et θ=20.
• Soit z(t)=cos (2πf0t), tracer Z(f), Z1(f) la TF de x(t).z(t) et Z2(f) la TF de y(t).z(t) sur le même graphe et comparer
Z(f) et Z1(f) puis Z(f) et Z2(f).
20
Rectangle
Solution 18 Triangle
X ( f ) = Aθ sin c ( fθ ) θ
Y(f ) = A sin c 2 ( f θ / 2 ) 16
2 14
1 1
Z ( f ) = δ ( f − f0 ) + δ ( f + f0 ) 12
2 2
A A
Z 1 ( f ) = θ sin c (θ ( f − f 0 )) + θ sin c (θ ( f + f 0 ))
10
2 2 8
6
A A
Z2( f ) = θ sin c 2 (θ ( f − f 0 )) + θ sin c 2 (θ ( f + f 0 ))
4 4 4
0
Exercices d'application -0.5 0 0.5
+∞ +∞ +∞
j 2πn t
1. Soit le signal x(t ) = δ (t − nT) ,
−∞
montrer que : e
−∞
T
= T δ (t − nT ) et
−∞
déterminer la TF de x(t)
+∞ +∞
1
Solution δ (t − nT ) → X ( f ) = T δ ( f − nf
n= −∞
TF
n =−∞
0 )
Principales propriétés de la TF
1
o Changement d'échelle : x(at ) →
TF
X ( f / a)
a
d n x(t ) TF
o Dérivation : n
→(2π j f ) n X ( f )
dt
x ( −t ) →
TF
X (− f )
o Inversion et conjugaison :
x * (t ) →
TF
X * (− f )
o Dirac : δ (t − t0 ) →
TF
e−2π j f t 0
δ (t)
TF
→1
1 1 1
o Echelon et signe: U (t) →
TF
+ δ( f ) Sgn( f ) =
2πjf 2 πjf
+∞ +∞
o Périodiques : C
n=−∞
n exp(2πjnf0t ) → Cnδ ( f − nf0 )
TF
n=−∞
+∞
1 +∞
o Peigne de Dirac : δ (t − nT ) → X ( f ) = δ ( f − nf0 )
TF
n =−∞ T n=−∞
1 1
o cos( 2πf 0 t ) →
TF
δ ( f − f0 ) + δ ( f + f0 )
2 2
1 1
o sin(2πf 0t ) → δ ( f − f0 ) − δ ( f + f0 )
TF
2j 2j
Série 0
3. Les signaux suivants sont-ils à énergie finie, à puissance moyenne finie, ou ni l’un, ni l’autre ? Calculer, dans
chaque cas, l’énergie totale et la puissance moyenne totale (a>0).
• Arect(n/N+1) Asin(2̟f0n) Asin(2̟f0n).U(n) U(n)
• n.U(n) Ae u(n)
-an Ae -an
4. Calculer l'autocorrélation des signaux suivants pour les cas continus et discrets
• Arect(n/N) Asin(2̟f0n) 5 δN(n) Bcos(2̟f0n)
7. Soit le signal x(n)=e-an U(n) transmis à travers le système h(n) dont la réponse impulsionnelle est h(n)
est donnée par h(0)=h(1)=h(2)=1/3.
• Calculer l’énergie et la puissance du signal d'entrée
• Calculer et tracer l’autocorrélation de h(n)
• Déterminer le signal y(n) résultant de la convolution numérique x(n)*h(n).
8.Etudier la linéarité, la causalité, l'invariance et la stabilité des systèmes définis par les équations aux
différences finies:
Linéarité Invariance Causalité Stabilité
y(n)= 2 .n. x(n-1)
y(n+1)=2 x(n+1)2
y(n)=|x(n-1)-0.5x(n+1)|
y(n)=1/en. x(n)
Solutions
2. U(n)=1/2(sgn(n)+1)
3.E=A2.(N+1) Pm=0, E=∞ Pm=A2/2, E=∞Pm=A2/4, E=∞ Pm=1/2
E=∞ Pm=∞, E=A /(1-e- ) Pm=0 , E=∞ Pm=∞,
2 2a E=2A N(1+N) Pm=0
2
4. A N.ΛN(k)
2 A /2.cos(2πf0k)
2 25δN(k) B2/2.cos(2πf0k)
5. x ( n ) = e − a ( n − n 0 ) + e − a ( n − n1) 6. f(n)*f(n)= E02.(n+1) pour n≥0 et 0 ailleurs
7. E=1/(1-e-2a) P=0 Rh(0)=1/3, Rh(±1)=2/9, Rh(±2)=1/9 y(n)=1/3 (x(n)+x(n-1)+x(n-2)
8. O-N-O-N N-O-O-N N-O-N-O O-N-O-N
Exercices supplémentaires
1. Donner l’expression du signal x(n) = Arect[(n-n0)/N+1]=A ∏ (n − n0 ) à l’aide du signal signe seulement.
N +1
3. Calculer et esquisser graphiquement pour les cas n0 < n1 et n0> n1 le produit de convolution z(n) =
x(n)*y(n) pour les cas suivants :
X(n) = A[δ(n+n0) + δ(n-n0)] et Y(n) = B δ(n) + ½B[δ(n+n1) + δ(n-n1)]
E 2 (n + 1 + N ) − N ≤ n ≤ 0
E N .Λ N (n ) = − E 2 (n − 1 − N ) 0 ≤ n ≤ N
4. La fonction triangulaire est définie de la manière suivante: 2
0
ailleurs
- Vérifier analytiquement et graphiquement la relation E2.N. ΛN(n) = E.ΠN+1(n) * E.ΠN+1 (n),
- En déduire l'auto-corrélation du signal et son énergie (devoir à rendre)
5. Soient x(n) et h(n) deux signaux numériques provenant respectivement de l’échantillonnage d’un signal
x et de la réponse impulsionnelle h d’un système :
x(n)={0,0,0,0.5,1.5,0.5,1.5,0,0,0} et h(n)={0.5,0.5}
- Calculer l’énergie de chaque signal
- Calculer et tracer l’autocorrélation de h(n)
- Calculer la séquence y résultant de la convolution numérique x(n)*h(n).
- Interpréter ces résultats du point de vue des plages de fréquences éliminées et conservées.
- Quel est le signal d’entrée qui permettrait de connaître le signal h(n) ?
- Proposer un signal h(n) permettant de réaliser un filtrage passe haut du signal x(n)
6. On désire déterminer la vitesse de propagation d’une impulsion dans un câble coaxial sans perte.
L’intensité est alors une fonction de la longueur de câble traversé ‘l’ et du temps ‘t’. Les mesures relevées donnent
l’évolution de cette intensité pour l=0 correspondant au signal en entrée du câble et pour une longueur l =100m
en fonction du temps.
A A
Câble coaxial
0.1 t en µs 0.5 t en µs
I. Rappels
x(0)
x(1) x(4)
Un signal discret s(n) est une suite de N échantillons régulièrement
x(-3) x(-2) x(-1)
espacés de Te secondes: x(0),x(Te),x(2Te),…........, x((N-1)Te) où x(2) x(3)
x(5)
Fe=1/Te est la fréquence d’échantillonnage. Le tracé graphique
d'un signal discrétisé en temps peut s'effectuer simplement à l'aide
de la fonction stem sous matlab. n
x(6) x(7)
- L’énergie d’un signal x(n) est fournie sous matlab par sum(x.^2). Concernant la puissance moyenne, il faut
diviser l’énergie par le nombre d’éléments de x(n).
- Pour la corrélation et la convolution, on utilisera, respectivement, les fonctions xcorr et conv. A noter que la
convolution ou la corrélation de x et h de durée respective N et M est un signal y(n) de durée (N+M-1)
- La fonction b=m+ s*randn(N,1) permet de générer un vecteur bruit b de distribution pseudo normale
(Gaussienne) de taille N de moyenne m et de variance s2 dont la puissance est Ps = m2 +s2 .
1. Le programme suivant permet de générer un Dirac en 0 : δ(n) =1 pour n=0 et vaut 0 ailleurs
clc ; clear all ; close all ;
t=-10:10;
x=[zeros(1,10),1,zeros(1,10)];
stem(t,x);
axis([-10 10 -0.5 1.5]);
title('Dirac');
xlabel('n');ylabel('Amplitude');
2. Le programme suivant permet de générer un échelon U(n)=1 pour n≥0 et 0 pour n<0
clc ; clear all ; close all ;
t=-20:20;
x=[zeros(1,20),ones(1,21)];
stem(t,x);
title('Echelon unite');
xlabel('n');ylabel('Amplitude');
3. Pour générer N=128 échantillons d'une sinusoïde de fréquence f0=1000, on peut procéder de la façon suivante,
choisir une fréquence d'échantillonnage : Fe = 8000 (le pas de temps Te=1/Fe) Créer le vecteur des temps : t =
(0:N-1)Te. Te. Calculer les échantillons: x = cos(2*pi*t*f0) ; Puis, regarder le résultat : plot(x) ou plot(t,x).
Ce qui nous donne :
clc ; clear all ; close all ;
N=128; fo=1000; Fe=8000; Te=1/Fe;
t=(0:N-1)*Te; x=cos(2*pi*fo*t);
plot(t,x) ;
figure; plot(x);
figure; stem(t,x);
4.Générer et visualiser le signal x composé de 50 échantillons d’une sinusoïde de fréquence f0 =0.1 avec fe=10.f0
- Calculer et afficher son auto-corrélation et comparer avec l'auto-corrélation théorique.
- Retrouver les caractéristiques du signal (puissance et fréquence) à partir de l'auto-corrélation.
- Générer le signal z=x+b ou b est un bruit avec (m=0 , s=0.5) . Calculer l'auto-corrélation du bruit et commenter.
- Calculer et visualiser l'auto-corrélation de z pour s=0.5, s=1 puis s=2 en commentant.
- Retrouver dans chacun des cas précédents la fréquence du signal à partir de l'auto-corrélation de z.
Dans la réalité, les signaux n'ont pas toujours une forme simple en raison de la nature de l'information qu'ils
portent. Dans de tels cas, la représentation du signal en fonction de la fréquence est peut s'avérer très utile. Pour
cela, on fait appel à la transformée de Fourier. Elle a pour but de mettre en évidence des caractéristiques du signal
non évidentes dans la représentation temporelle : les propriétés fréquentielles (spectrales). L’utilisation de cette
description fréquentielle permet de caractériser simplement les filtres linéaires, et faciliter leur étude.
Dans le but de calculer la transformée de Fourier X(f) d’un signal x(t) à l’aide d’un ordinateur, celui-ci
n’ayant qu’un nombre limité de mots de taille finie, on est amené à discrétiser le signal (échantillonnage), à
tronquer temporellement ce signal et à discrétiser sa transformée de Fourier[1].
Te
Soit x(t) un signal analogique de transformée de Fourier X(f). Echantillonner le signal x(t) consiste à choisir
une fréquence Fe et de construire un nouveau signal avec les x(nTe) avec n un entier et Te=1/fe.
Cette expression montre que le spectre Xe(f) est périodique de période feet
-fmax fmax f
qu’il est la somme des répliques (copies) du spectre original X(f) décalées
de nfe . L'échantillonnage dans le domaine temporel se traduit par une
Xe(f)
"périodisation" de période fe dans le domaine fréquentiel.
Théorème de Shannon
On considère que x(t) est un signal réel dont le spectre est borné en fréquence, de fréquence maximale fmax
soit ∀ f > f max X(f) = 0
Xe(f)
Deux cas de configurations :
- fe> 2 fmax
- fe< 2 fmax
Remarques :Si le support du spectre X(f) n’est pas borné (s’étale sur l’axe réel) il y a un repliement du spectre des
échantillons (aliasing), on ne peut pas isoler le spectre original à partir de celui des échantillons.Dans la pratique,
on ne peut pas se contenter de prendre une fréquence d´échantillonnage égale à la fréquence de Nyquist (fe/2>fmax),
mais en prendre une supérieure car on ne peut réaliser un filtre passe-bas idéal avec une fréquence de coupure
très nette. Par exemple, pour numériser la parole dans le réseau téléphonique, on utilise une fréquence
d'échantillonnage 8kHz alors que le spectre de la voix est compris entre 300Hz et 3400Hz.
Filtre anti-repliementLes signaux étudiés en réalité sont rarement à support fréquentiel borné, c’est-à-dire que
fmax = infinie. C’est par exemple le cas d’un signal rectangulaire périodique dont les raies fréquentielles s’étendent
à l’infini ou encore un signal bruité. Ceci implique que quelle que soit la fréquence d’échantillonnage il y aura
repliement de spectre puisque fe> 2 fmax = ∞ est une condition impossible à réaliser. Pour remédier à ce problème,
on utilise à l’entrée d’un système numérique un filtre passe-bas appelé filtre anti-repliement ou anti-aliasing. Ce
filtre est analogique, idéalisé il doit avoir un gain de 1 sur une bande de fréquence Fe, centrée en zéro. Son rôle
va être de limiter le contenu spectral du signal à la partie utile. Il va participer aussi à limiter l’influence du bruit
[3].
Exemple d'application
On échantillonne un signal sinusoïdal de fréquence 200Hz avec une fréquence d’échantillonnage fe = 500Hz
puis avec fe = 300Hz. Quel signal obtient-on lors d’une reconstruction parfaite dans les deux cas ?
= cos 2
Solution
= cos 2
fe=500Hz t)avec f0=200Hz
fe=300Hz t)avec f0=100Hz
Lorsque le signal à traiter n’est plus analogique mais numérique, la relation de la TF devient :
. +∞
TF {x (t )} = X ( f ) = x(t )e − 2π j f t dt
−∞
+∞ +∞ +∞ +∞
L'échantillonnage périodise le spectre du signal avec une période de répétition feainsi Xe(f)=Xe(f+fe), par
ailleurs, l'amplitude est multiplié par un facteur fe. Sachant que tout signal périodique peut être décomposé en
séries de Fourier, on a :
fe / 2
+∞ 1
x(nTe ) = C n = X ( f ) e 2π j n Te f df
Xe( f ) = C e
n = −∞
n
− 2π j n Te f Avec
fe − fe / 2
e
Cette transformée de Fourier appliquée aux signaux discrets est donc une fonction à fréquence continue,
périodique de période fe. Il est d’usage de la représenter sur un intervalle de longueur fe, de -fe/2 à +fe/2.
Cependant, si on veut calculer la TF d'un signal discret à l'aide d'un calculateur, on se retrouve confronté
aux problèmes suivants : Le calcul de la TF nécessite une infinité de points de mesures x(n) (pas toujours possible
dans la pratique : contraintes temps réel, etc.). En outre, le calculateur ne peut calculer une TFTD car sa réponse
fréquentielle est forcément discrète = un nombre fini de points fréquentiel alors que f varie continûment [5]. La
solution est de limiter la durée de x(n) i.e. considérer un nombre fini N de points temporels et de discrétiser la
fréquence (considérer un nombre fini L de points fréquentiels) d'où la TFD.
Cette transformée, popularisée par son calcul rapide (TFR ou FFT : Fast Fourier Transform), fait
correspondre une suite de N valeurs à une autre de suite de N valeurs numériques également.
On considère un signal numérique s(n) défini par N échantillons temporels, obtenus par échantillonnage
avec la période Te. La numérisation du signal X(f) passe par l'échantillonnage de X(f). On divise l'intervalle fe par
N, ainsi X(f) est échantillonné à la cadence ∆f=fe/N=1/NTe. Ce dernier résultat entraîne une périodicité du signal
temporel de T0=1/∆f = NTe.
x(t) X(f)
t f
x(nTe) Xfe(f)
fe
Te t f
XT0(nTe)
Xfe(k∆f)
t ∆f=fe/N f
T0=NTe
1 N / 2−1
1 N / 2−1
x ( n) =
X k e 2π j n k / N
N k =− N / 2
Et x ( n) =
fe
X
k =− N / 2
e (k ) e 2π j n k / N .∆f =
N
X
k =− N / 2
e ( k ) e 2π j n k / N
N −1
Ce qui nous permet d'obtenir la TFD et la TFD inverse :
X ( k ) =
n =0
x ( n ) e − 2π j n k / N
1 N / 2 −1
x(n) =
N k =− N / 2
X ( k ) e 2π j n k / N
On remarquer aisément que l’on perd toute notion de temps relatif aux échantillons. Nous obtenons au final
une relation entre une suite indexée par une variable entière n et une suite indexée par k. Les N termes X(k)
correspondent à une approximation (à un facteur multiplicatif Te près) de latransformée de Fourier de ce signal
aux N points de fréquence k ∆f, avec k entre 0 et N −1, c'est-à-dire f entre 0 et fe.
Exemple
+∞
- Calculons d'abord la TFTD, ce qui nous donne X ( f ) = x(nT )e
n=−∞
e
−2π j f nTe
= 1 + e −6πjfTe = 2 cos(3πfTe )e −3πjfTe
- Calcul maintenant la TFD sur N=4 échantillons (4 échantillons de la TFD à partir de 4 échantillons du signal)
3
x(n)e = x(0) + x(1) + x(2) + x(3) = 2
0
n3=0
x(n)e −πjk / 2 = x(0) + x(1)e −πj / 2 + x(2)e −πj + x(3)e −3πjk / 2 = 2 cos(3π / 4)e −3πj / 4
3
X (k ) = x(n)e −2πjnk / 4 = n =0 3
n =0
n =0
x(n)e −πjnk = 2 cos(3π / 2)e −3πj / 2
3
n =0
x(n)e −3πjnk / 2 = 2 cos(9π / 4)e −9πj / 4
2 2
TFTD TFTD
1.8
1.5
1.6
1
1.4
0.5
1.2
Module
Phase
1 0
0.8
-0.5
0.6
-1
0.4
-1.5
0.2
0 -2
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
On peut observer que les quatre échantillons de la TFD (en rouge) se superposent à la courbe de la TFTD (en
bleu). On confirme que la TFD n'est que l'échantillonnage de la TFTD limitée à N. On note en outre, que la
précision fréquentielle est de ∆f=fe/N. Pour améliorer cette précision, il faudrait diminuer le pas en fréquence.
Remarque x(n) est périodique de période N et X(k) est aussi périodique de période N :
N −1 2π ( n + N ) k N −1 2πnk 2πNk
1 j 1 j j
x(n + N ) =
N
X ( k )e
k =0
N
=
N
X ( k )e
k =0
N
e N
= sn
Sachant que x(n) et X(k) sont calculés sur le même nombre de points N, on peut augmenter la précision,
par la technique du zéro-padding : on calcule la TFD sur un nombre NF pouvant être bien plus grand que le
nombre de points N disponible du signal (NF>>N). La figure suivante en donne un exemple (Voir TP n°1).
4 4
TFD d une porte (largeur 4) calculee sur 4 points TFD d une porte (largeur 4) calculee sur 16 points
3.5 3.5
3 3
2.5 2.5
2 2
1.5 1.5
1 1
0.5 0.5
0 0
-0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 -0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5
4 4
TFD d une porte (largeur 4) calculee sur 32 points TFD d une porte (largeur 4) calculee sur 64 points
3.5 3.5
3 3
2.5 2.5
2 2
1.5 1.5
1 1
0.5 0.5
0 0
-0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5 -0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5
Autres propriétés :Toutes les propriétés se déduisent des propriétés de la transformée continue en se rappelant
que chaque signal manipulé, de durée finie, doit être considéré comme une période d'un signal périodique, et
cela en temps et en fréquence. La conséquence en est que la translation d'un signal ( lors d'une convolution ou de
corrélation) se traduit par un décalage circulaire [6]
−2π j k m
o Décalage temporel : x ( n − m ) TFD
→ X ( k ) e N
2π j k0 n
o Décalage fréquentiel : x ( n ) e N
TFD
→ X ( k − k 0 )
1
o Changement d'échelle : x(an)
→ TFD
X (k / a)
a
o Convolution périodique : x (n) ∗ h(n) TFD
→ X (k ).H (k )
o δ (n) TFD
→ 1 δ (n − m) TFD
→ e −2πjmk / N
N −1 N −1
1
o Par ailleurs, l'énergie se conserve: x(n) = X (k )
2 2
0 N 0
La TFD est restée un outil peu utilisé jusqu’à l’apparition d’algorithmes « rapides » permettant son calcul.
Le plus connu est du à Cooley et Tuckey et date de 1965. Le calcul direct de la TFD sur N points nécessite N2
opérations complexes. L’algorithme proposé réduit à Nlog2(N) le nombre d’opérations. Ainsi, pour N = 1024, le
temps de calcul de l'algorithme rapide peut être 100 fois plus court que le calcul utilisant la formule de définition
de la TFD.Pour en donner un exemple, prenons le cas de N=8, il faut calculer :
N −1 2π n k 7 2π n k
−j −j
S k = sn e N
= sn e N
n=0 n=0
2π 7
Soit : W N = exp − j alors S k = s nW Nnk
N n =0
Pour N=8, explicitons la relation précédente :
Les facteurs Wn présentent un certain nombre de propriétés dont certaines sont mises à profit dans
l’algorithme : n+ N
nN
WN = 1, WN
N /2 n
= −1, WN = WN
S 0 W80 W80 W80 W80 W80 W80 W80 W80 s0
0 Imaginaire
S1 W8 W81 W82 W83 W84 W85 W86 W87 s1 W86
S 2 W80 W82 W84 W86 W80 W82 W84 W86 s2 W85 X W87
0 X X
S 3 = W8 W83 W86 W81 W84 W87 W82 W85 s3 W84 W80=W88
X X
S W 0 W 4
W 0
W 4
W 0
W 4
W 0
W84 s4 Réel
4 80 8
5
8
2
8
7
8
4
8
1
8
6
X
S 5 W8 W 8 W 8 W 8 W 8 W8 W8 W83 s5 W83
X
W81
S W 0 W 6
W 4
W 2
W 0
W 6
W 4
W82 s6 W82
6 8 8 8 8 8 8 8
S 7 W80 W 8
7
W 8
6
W 8
5
W 8
4
W8
3
W8
2
W81 s7
L’algorithme suppose que N est pair : posons N=2 P. Introduisons les 2 sous-suites de sn en fonction de la
parité de n.
u n = {s2 n }
n = 0 ,..., P −1 N −1 P −1 P −1
S k = snWNnk = uiW22Pik + viW2(P2i +1) k
vn = {s2 n +1}n =0,..., P −1 n=0 i =0 i =0
P −1 P −1
S k = uiWPik + WNk viWPik
On obtient ainsi : i =0 i =0
S k = U k + W Vk k
N
Par ailleurs, N −1 P −1 P −1
S ( k + P ) = s nWN( k + P ) n = uiW22Pi ( k + P ) + viW2(P2i +1)( k + P )
n =0 i =0 i =0
P −1 P −1
S ( k + P ) = uiWPikWPiP + WNkWNP viWPikWPiP
i =0 i =0
P −1 P −1
S ( k + P ) = uiWPik .1 + WNk (−1) viWPik .1
i =0 i =0
S ( k + P ) = U k − W Vk
k
N
Le calcul de la FFT revient donc à calculer Uk et Vk qui sont les TFD sur P points des suites de termes de
rang pair et impair. On s’aperçoit qu’il ne reste qu’à exprimer les Uk et Vk.
s0 U0 S0
+ + +
s2
W2
- + U2 + S1
0
+ W4 - U4 +
s4 W2 S2
- W4
2 - U6 +
s6 S3
+ + W -
s1 V1 S4
- + 8 -
s3 W2 V3 0
S5
+ 0
- -
W4 V5
s5 W2
- - W8
1
- S6
2
s7 W4 V7 2 S7
W8
3
Or, ce sont des TFD sur P points, qui peuvent subir le même schéma que précédemment à condition que P
soit pair. On peut réitérer le processus à chaque sous-niveausi N est une puissance de 2 (Dans le cas contraire, on
rajoute autant de 0 que nécessaire pour obtenir une puissance de 2).
L’algorithme papillon de FFT dit aussi à entrelacement temporel peut s’écrire sous forme matricielle. On
obtient :
S 0 1 W8 0 0 0 0 0 0 0 1 0 W40 0 0 0 0 0 1 0 0 0 W2
0
0 0 0 s 0
S
1 1 − W8
0 0 0
0 0 0 0 0 0 0 1 0 W4 0 0 0 0 0 1 0 0 0 W2 0 0 s1
S 2 0 0 1 W8
2
0 0 0 0 1 0 − W40 0 0 0 0 0 0 0 1 0 0 0 W2
0
0 s 2
0
S 3 = 0 0 1 − W8
2
0 0 0 0 0 1 0 − W40 0 0 0 0 0 0 0 1 0 0 0 W2 s 3
× ×
S 4 0 0 0 0 1 W8
1
0 0 0 0 0 0 1 0 W42 0 1 0 0 0 − W2
0
0 0 0 s 4
1 − W8 − W2
1 0
S 5 0 0 0 0 0 0 0 0 0 0 0 1 0 W42 0 1 0 0 0 0 0 s5
S 0 0 0 0 0 0 1
3
W8 0 0 0 0 1 0 − W42 0 0 0 1 0 0 0 − W2
0
0 s 6
6
7 0
S 0 0 0 0 0 1 − W8 3 0 0 0 0 0 1 0 − W42 0 0 0 1 0 0 0 − W2 s 7
0
Le nombre d’éléments d’une séquence transformée par la TFD est implicitement limité, la fenêtre
intrinsèque à la transformée discrète de Fourier est donc la fenêtre rectangulaire de durée T0=NTe.
TFTD x(n) . Re ct (n) ≈ X ( f ) * N sin c( fT0 )
T0 = NTe
La troncation du signal échantillonné par une fenêtre de largeur T0 a pour effet de convoluer le spectre avec
un sinus cardinal qui s’annule tous les 1/T0 avec T0 = NTe soit tous les fe/N.
Exemple :
En rouge est illustré X(k) le module de la TFD de {x(n) = e2πjf0n}, pour n= {0, . . . ,N − 1}, avec N = 16 et f0 = 0,2. En
bleu X(f), le module de la TFTD de x(n). L’allure de X(f) fait apparaître un lobe principal de largeur 2/N (ou
2fe/N) autour de la fréquence f0 et des lobes secondaires de largeur 1/N (ou fe/N).
16
TFD
14 TFTD
12
10
Phase
0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Dans le cas général, le spectre, obtenu par transformée de Fourier discrète, est donc un ensemble de fonctions
sinc(T0f ) centrées sur les fréquences qui composent le signal théorique initial.
20 16
TFD
18
14 TFTD
16
12
14
12 10
3 fréquences
3 sinc
10 8
8
6
6
4
4
2
2
0 0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Exemple:
Le tracé du spectre par TFD montre que si l'écart en valeur absolue entre f0 et f1 est supérieur à fe/N, il sera
possible de distinguer les deux fréquences sur le tracé. Cette résolution en fréquence est liée au nombre de points
N du signal (voir les 2 figure ci-dessous)
16 16
fe=1 N=16 f0=0.16 f1=0.2 A0=A1 fe=1 N=16 f0=0.18 f1=0.3 A0=A1
14 14
12 12
10 10
8
8
6
6
4
4
2
2
0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Un masquage du lobe principal de la composante en f1 peut survenir en raison des ondulations présentes
dans le spectre de A0exp(2j̟f0n). Une « fréquence » d’amplitude faible au voisinage d’une d’amplitude plus élevée
sera masquée par le premier lobe secondaire.
16
fe=1 N=16 f0=0.18 f1=0.3 A0=10*A1
14
12
10
0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Ce phénomène se révèlera gênant lorsque le spectre est composé de nombreuses raies, proches les unes des
autres. L'importance des lobes secondaires peut être réduite par l'emploi d'autres fenêtres.
Mais cela se fera au détriment de la séparation de « fréquences » très voisines mais d’amplitude semblables
car les 2 raies seront confondues dans un lobe principal plus large (la fenêtre rectangulaire possède le lobe
principal le plus étroit).
6. Fenêtres de pondération
Lors de l’analyse spectrale d’un signal de longue durée, nous n’avons accès, en pratique, qu’à une portion
limitée de ce signal. Le spectre obtenu correspond donc au spectre du signal à analyser auquel une fenêtre a été
préalablement multipliée [7]. Pour ne pas altérer le spectre original, il faudrait que W(f) (spectre de la fenêtre) se
rapproche le plus possible d’une distribution de Dirac. La distribution de Dirac étant l’élément neutre du produit
de convolution. Il y a deux éléments importants pour se rapprocher de la distribution de Dirac. La finesse du lobe
10
principale et la hauteur des lobes secondaires.
9
7
est grande, c’est-à-dire que l’on peut séparer des raies proches. Et plus les
6
lobes secondaires sont élevés plus on dégrade la forme du spectre et la
5
détection d’un signal d’amplitude faible en présence d’un signal
4
d’amplitude élevée sera ardue [8]. 3
0
- Fenêtre Rectangulaire -0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5
Pour la fenêtre rectangulaire, on voit que la finesse du lobe principale (2 fe/N).peut être réglée par le nombre
d’échantillons N. Ainsi, plus on observe le signal longtemps, plus la résolution du spectre augmente ce qui parait
logique. Par contre, λR varie très peu en fonction de N (-13dB) donc toujours une distorsion de spectre.
- Fenêtre Triangulaire 10
TFD Fen réctangulaire
9 TFD Fen Triangulaire
Pour obtenir la transformée de Fourier de la fenêtre triangulaire de
8
largeur N, rappelons que la convolution de deux signaux
7
rectangulaires donne un signal triangulaire. Ainsi, on peut
6
exprimer cette fenêtre sous la forme d’une convolution de deux
5
rectangles de largeur N.Te/2. On observe une atténuation des lobes
4
secondaires (-24dB) par rapport à la fenêtre rectangulaire[8].
3
Malheureusement, ceci se fait au prix de l’élargissement du pic
2
central (4fe/N).
1
0
-0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5
- Autres fenêtres
10
Rectangle TFD Rectangle
Hanning 9 TFD Hanning
1 Hamming TFD Hamming
8
Blackman TFD Blackman
7
0.8
6
0.6 5
0.4
3
2
0.2
1
0 0
5 10 15 20 25 30 35 40 45 50 -0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5
Dans un problème d’analyse spectrale, on utilise généralement plusieurs fenêtres l’une après l’autre afin d’obtenir
un bon compromis résolution/déformation.
FEI, USTHB [assiakourgli@gmail.com / http://perso.usthb.dz/~akourgli/
33
Traitement Numérique du Signal Master ST M1 2019/2020
Série n°1
1. On considère un signal de parole de durée 1mn et ayant une bande passante de 10 kHz échantillonne comme
suit :
xe(t)
xr(t)
x(0)
2.Calculer la transformée de Fourier à temps discret (TFTD) de x(n)=δ(n) + 6δ(n –1) + 3δ(n –2)
0.7 0.7
- Calculer l'énergie du signal.
0.6 0.6
- Déterminer et tracer l'auto-corrélation du signal S1.
0.5 0.5
- On décale le signal S1 de 2 vers la droite, sans
0.4 0.4
faire de calcul, tracer son auto-corrélation.
- Déterminer la TFTD de S1 et tracer son module. 0.3 0.3
- Calculer la TFD pour tout N. Tracer la TFD pour N=6. 0.2 0.2
7. Etant donné les signaux s(n)={1,-2,3,2} et v(n) ={-2,1,2,3}, déterminer w(n)=s(n)*v(n) par
- la méthode directe
- la méthode de la TFD
1.5 1.5 1
8
7
1 0.8
6
1
0.5 0.6
5
4
0 0.4
0.5 3
-0.5 2 0.2
1
-1 0
0 0 5 10 15 20 25 30 0 5 10 15 20 25 30 0 5 10 15 20 25 30
0 5 10 15 20 25 30
1.5 1.5
30
14
1 1
25 12
0.5 0.5 20 10
8
15
0 0 6
10
4
-0.5 -0.5
5
2
-1 -1 0 0
0 5 10 15 20 25 30 0 5 10 15 20 25 30 0 5 10 15 20 25 30 0 5 10 15 20 25 30
1.5 1.5 12 16
10 14
1
12
1 8
0.5 10
6 8
0
4 6
0.5
4
-0.5 2
2
0
0 -1 0
0 5 10 15 20 25 30 0 5 10 15 20 25 30 0 5 10 15 20 25 30 0 5 10 15 20 25 30
1.5 1.5
15
7
1
6
1
0.5 10 5
4
0 3
0.5 5
2
-0.5
1
-1 0 0
0 5 10 15 20 25 30 0 5 10 15 20 25 30 0 5 10 15 20 25 30 0 5 10 15 20 25 30
9. Un signal analogique est échantillonné à une fréquence fe=7500 Hz et N échantillons sont collectés.
- Quelle est la résolution fréquentielle pour N=1250?
- Déterminer N permettant d'atteindre une résolution fréquentielle de 4.5 Hz.
10. Un signal analogique est échantillonné à une fréquence fe= 500Hz et N=980 échantillons sont collectés.
On veut connaître la valeur du spectre à 120 Hz.
- Quel indice k de la TFD est le plus proche de 120Hz?
- Quel est le nombre minimum de 0 à rajouter obtenir une valeur de la TFD à 120 Hz exactement?
- Donner alors la valeur de l'indice k correspondant.
11. Soit x(n) = sin(2̟ f1 n) + 0.5 sin(2̟ f2 n) avec 0 ≤ n ≤ 127,f1 = 0.223, f2 = 0.240
- On veut employer un fenêtrage, calculer la largeur du lobe principale de chaque fenêtre
- Déterminer alors le fenêtrage permettant de visualiser les 2 sinusoïdes.
12. Soit x(n)=A0 e2πjf0n + A0 e2πjf1n pour n (0:N-1) où f0=0.25, et f1=0.70 et fe=2 N=20.
- Donner l’allure de la TFTD et de la TFD.
- Quel est le but du fenêtrage ?
- Quelle est la meilleure fenêtre à utiliser pour ce signal ? Justifier.
Solutions
50
2. X(f)=1+6.e-2πjfTe+3.e-4πjfTe
sin(pi*N*f)/sin(pi*f)
N*sinc(N*f))
40
30
-10
-20
0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5
1 1 − e − aN
4. TFTD X(f ) = TFD X ( k ) =
1 − e −a −2π j fTe 1 − e − a − 2π j k / N
5. E=3, Rs1(0)=2, Rs1(1)= Rs1(-1)=1, identique, SS1(f)=2cos(πfTe)e-πjfTe, SS1(k)=2cos(πk/N)e-πjk/N
SS1=[2, 1.73 1,0, 1, 1.73] , zéro padding, Rajouter des points entre la TFD de N=6.
6. N=2, X(k)=[1, -1] N=3, X(k)= 3, 1 (−1 − 3 j ), 1 (−1 + 3 j ) N=4, X(k)= [6 , 2 j − 2 , − 2 , − 2 − 2 j]
2
2
Réel
W4
0
W 4
3
W4
6
4
9
W W 4
0
W4
3
W 4
2
W 1 j − 1 − j
4
1
-1
W41= -j X
Exercices supplémentaires
1. Calculer la TFD du signal x(n)=1 pour 0≤n≤6 et x(n)=0 pour 7≤n≤29 puis représenter la.
3. Calculer la transformée de Fourier discrète (TFD) de la suite x(n) formée de N = 8 points (n ∈ [0,7]),obtenue en
échantillonnant à la fréquence fe = 16 Hz le signal s(t) = 2sin(8̟t)+8cos(4̟t)
x(n) y(n)
4.On considère le système suivant : h(n)
On suppose que H(0)=H(1)=1 et H(2)=H(3)=0.
- Déterminer et tracer h(n) pour 0 ≤ n≤ 3
- On pose H(4)=H(5)=……= H(N-1)=0. Donner le nom de cette technique.
- Quelle est le but de cette opération ? Donner, alors, l’allure approximative de h(n) pour N assez grand.
0
0 5 10 15 20
6. On souhaite calculer la TFD du signal suivant : e2πjf0n +0.1 e2πjf11n avec f0=0.26, f1=0.3 , N=100 (fe=1) et NF=1024.
A cette fin, on teste différentes fenêtres (Rectangulaire, Hamming et Hanning). Les TFD obtenues sont illustrées
ci-dessous.
- Identifier les fenêtres utilisées (en justifiant vos réponses)
- Quelle est la meilleure fenêtre pour ce signal ? Justifier
0.7 0.7 1
0.9
0.6 0.6
0.8
0.6
0.4 0.4
0.5
0.3 0.3
0.4
0.2
0.1 0.1
0.1
0 0 0
-0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5 -0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5 -0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5
7. On donne les signaux et les TFD suivantes, associer chaque signal à sa TFD
9.Pour chacun des cas suivants, on désire choisir le fenêtrage adéquat permettant d’analyser le contenu
fréquentiel (fe=1,N=20):
- x(n)=2 exp(2πjf0n) +exp(2πjf1n) f0=0.3f1= 0.36
- x(n)=exp(2πjf0n) +exp(2πjf1n) f0=0.3f1= 0.46 - x(n)=2 exp(2πjf0n) +0.2 exp(2πjf1n) f0=0.3f1= 0.42
8
Solutions: 7
1. X (0) = 7 X (k ≠ 0) = e 5
sin(πk / 30) 4
0
0 5 10 15 20 25 30
S 0 W80 W80 W80 W80 W80 W80 W80 W80 s0 W80 W80 W80 W80 W80 W80 W80 W80 s0
0 Imaginaire
1 2 3 4 5 6
W87 s1 W80 1 2 3 4 5 6
W87 s1 W86
S1 W8 W8 W8 W8 W 8 W 8 W 8 W8 W8 W8 W8 W
8 W8
3. S 2 W80 W8
2
W8
4
W8
6
W 8
8
W8
10
W8
12 14
W8 s2 W80 W8
2
W8
4
W8
6
W8
0
W
8
2
W8
4
W8 s 2
6 W85 X W87
0 3 6 9 12 15 18 21 X X
S3 = W8 W W W W W W W8 s3 W80 W 3
W 6
W 1
W 4
W7
W 2
W85 s3
8 8 8 8 8 8
= 8 8 8 8 8 8
S W 0 W 4
W 8
W 12
W 16
W 20
W 24
W828
s W 0 W 4
W 0
W 4
W 0
W4
W 0
W8 s 4
4 W84 W80
4 80 8
5
8
10
8
15
8
20
8
25
8
30 35
4 80 8
5
8
2
8
7
8
4
8
1
8
6
X X
S5 W8 W8 W8 W8 W8 W8 W8 W8 s5 W8 W8 W8 W8 W8 W
8 W8 W83 s5 Réel
S W 0 W 6
W 12
W 18
W 24
W 30
W 36
W842 s6 W80 W 6
W 4
W 2
W 0
W6
W 4
W8 s6
2
6 8 8 8 8 8 8 8 8 8 8 8 8 8 X X
W83 X W81
S 7 W80 W8
7
W8
14
W8
21
W8
28
W8
35
W8
42
W849 s7 W80 W8
7
W8
6
W8
5
W8
4
W
8
3
W8
2
W81 s7
W82
[
sn = 8 2 + 4 2 0 − 2 − 4 2 − 8 2 − 4 2 0 − 2 + 4 2 W8k = 1
] 2
2
(1− j) − j
2
2
(−1− j) −1
2
2
(−1+ j) j
2
2
(1+ j)
4. Voir examen 14/15
5.Voir interro 15/16
6. Hamming, Hanning, Rectangulaire, Meilleure : Hamming
7. (1,3), (2,8), (3,4), (4,2), (5,7), (6,5), (7,1) et (8,6)
8. Voir examen 15/16
9. Voir examen 16/17
Rappels
3. Programmes à réaliser
I. Reprendre l'exemple vu en cours x(n)= [1 0 0 1], calculer, afficher et visualiser sa TFD pour NF=4 puis NF=32
et commenter.
- Pourquoi parle-t-on de fenêtre réctangulaire ? A-t-on vraiment multiplié par une fenêtre ?
- Que se passe-t-il si l'on remplace f0 par 1600?
- On remet f0 à 1680 et on rajoute au programme les lignes suivantes
fen=hanning(N);
xx=x.*fen'; y=fft(xx,NF); yy=fftshift(y/N);
plot(axe_f,abs(yy),'r'); grid; title('TFD Fenêtre Hanning');
Remarque Deux raies d’un spectre sont considérées comme séparables, si le maximum de l’une correspond au
premier minimum nul de l’autre soit |f1-f2| >Largeur du lobe principal/2=L∆f/2 (Voir tableau page 31) .
Rappels
1. Transformée de Fourier Discrète
La TFD d’ordre N d’un signal numérique s(kTe), k=0…N-1 est définie par :
k N −1
X (k ) = X f e = x(nTe )e−2π j k n / N , k = −N / 2.......N / 2
N n=0
Quelques astuces
1. Pour des figures séparées, aller Outils ->Préférences->Console IPython->Graphise->Sortie: Automatique.
Fermer spyder et ré-ouvrir.
2. Pour effacer toutes la variables avant exécution Préférences->Exécuter->Supprimer toutes les variables avant
exécution
3. Ctl+L efface la console
E=0.0
for i in x:
E+=pow(i,2)
t = np.linspace(0.0, N*Te, N)
z = np.exp(2*np.pi*1j*f0*t)
Remarque Deux raies d’un spectre sont considérées comme séparables, si le maximum de l’une correspond au
premier minimum nul de l’autre soit |f1-f2| >Largeur du lobe principal/2=L∆f/2 (Voir tableau page 31) .
La transformée de Fourier est un outil précieux d'analyse et de traitement des signaux. Cependant, dans
certains problèmes (comme le filtrage numérique), les limites de la TF sont vite atteintes. La transformée en Z,
qui s'applique aux signaux discrets, généralise la TF et permet de dépasser ces limites [10]. Elle est tout-à-fait
analogue à la transformée de Laplace, mais plus facile à utiliser. Ce type de transformée permet de décrire
aisément les signaux à temps discret et la réponse des systèmes linéaires invariants soumis à des entrées diverses.
C'est un outil qui permet de calculer la réponse impulsionnelle d’un système linéaire invariant décrit par une
équation aux différences finies. Elle permet l'interprétation directe des caractéristiques des signaux et des filtres
dans le domaine des fréquences [9].
1. Transformée en Z
+∞
• Définition : La TZ est la généralisation de la TFTD ( X ( f ) = x(nT )e
n=−∞
e
−2π j f nTe
).
En effet, comme cette transformation est une série infinie, elle n’existera que pour lesvaleurs de z pour lesquelles
cette série converge.
Exemples :
- x(n) = δ(n)X(z) = 1,
X ( z ) = 1 + 2 z −1 + 3z −2 + 5z −3 + 2 z −5 RDC=¢-{0}
La série des puissances introduite dans l’équation de définition de la TZ ne converge que pour un sous-ensemble
du plan complexe. Ce sous-ensemble est appelé région de convergence (RDC) ou domaine de convergence. Une
région de convergence correspond à l’ensemble des valeurs de z telles que X(z) soit définie et à valeurs finies.
Spécifier le domaine de convergence de la transformée est tout aussi important que la transformée elle même [9].
• Condition d'existence : La transformée existe si la série converge. Pour cela, on utilise le critère de Cauchy
+∞
lim U n
n →∞
1/ n
< 1 sur la convergence des séries géométriques S = U
n =0
n
L'ensemble des valeurs de la variable complexe z pour lesquelles la série converge est appelée Région De
Convergence (RDC): +∞
RDC = z ∈ C / x(n).z − n < +∞
n = −∞
Exemple : TZ{u(n)}
+∞ N −1
1 − z −N 1
U ( z) = z −n
= lim z −n
= lim −1
, la limite est finie si z < 1 ⇔ z > 1 U ( z ) = pour z > 1
N →∞ N →∞ 1 − z −1 1 − z −1
n =0 n =0
De façon générale, on montre que la RDC est un anneau de convergence centré sur l'origine défini par :
−1 / n
Im(z)
avec r1 = lim x( n) et r2 = lim x ( − n )
1/ n
r1 < z < r2
n → +∞ n → +∞
r2 Re(z)
RDC
Exemples
+∞
a n si n ≥ 0 z
- Soit a > 0, x ( n ) = X(z)= a n
z −n =
z−a
, convergente pour z > a .
0 si non n =0
0 si n ≥ 0 z
- Soitb > 0, y ( n ) = Y(z)= − b n
z −n =
z −b
, convergente pour z < b .
− b si n < 0
n
n<0
a n si n ≥ 0 z z
- Soient a >0 ,b> 0, w ( z ) = W(z)= − , convergente pour b > z > a .
b si n < 0
n z −a z −b
2. Propriétés de la TZ
Les propriétés qui sont les plus utilisées sont résumées comme suit :
k −1
- Théorème de l'avance : x(n + k ) → z . X ( z) −
TZ k
x(n) z
n =0
k −n
RDC : identique
z
- Multiplication par an : a x(n) → X
n TZ
RDC : a .r1 < z < a .r2
a
dX ( z )
- Théorème de dérivation : n. x ( n ) →
TZ
− z. RDC : identique
dz
Exemples
1 jw 0 n
1. x ( n ) = cos( wo n)U ( n) = (e + e − jw 0 n )U ( n)
2
1 1 1 1 − cos( w0 ) z −1
X (z) = + = avec z >1
2 1 − e jw0 z −1 1 − e − jw 0 z −1 1 − 2 cos(w0 ) z −1 + z −2
2. Calculer la transformée en z des fonctions discrètes suivantes. Vérifier que les théorèmes de la valeur initiale
et finale s'appliquent : x(n) = 0,8n u(n) et y(n) = n0,8nu(n).
z
d
z z − 0 .8
Y (z) = − z
0 .8 z
X (z) = =
z − 0.8 dz (z − 0 .8 )2
z 0 .8 z
x(0) = lim =1 y (0) = lim =0
z →∞ z − 0.8 z →∞ ( z − 0 . 8 )2
z( z − 1) 0 .8 z ( z − 1)
x(∞) = lim =0 y ( ∞ ) = lim =0
z →∞ z − 0.8 z→∞ ( z − 0 . 8 )2
Quelques TZ
− a U ( − n − 1)
n 1 z < a
1 − az −1
cos(ω 0 nTe )U ( n ) 1 − z −1 cos(ω 0Te ) z >1
1 − 2 z −1 cos(ω 0Te ) + z − 2
sin(ω0 nTe )U (n) z −1 sin(ω0Te ) z >1
1 − 2 z −1 cos(ω0Te ) + z −2
a n cos( ω 0 nT e )U ( n ) 1 − az −1 cos(ω 0Te ) z > a
1 − 2az −1 cos(ω 0Te ) + a 2 z − 2
a n sin(ω 0 nTe )U ( n ) az −1 sin(ω 0Te ) z > a
1 − 2az −1 cos(ω 0Te ) + a 2 z − 2
Les systèmes linéaires invariants décrits par une équation aux différences finies possèdent une
transformée en Z rationnelle c'est ainsi que celles-ci vont s’écrire comme le rapport de deux polynômes en z-1.
M N M N
On peut caractériser un système LI par h(n) ou par la transformée en Z (H(z)) de sa réponse impulsionnelle
h(n), encore appelée fonction de transfert du système.
b .z i
−i
b N ( z ) b0 M − N ∏ (z − z ) = K z
N
∏ (z − z )
N
H (z) = i =0
= 0 z M −N = z i =1 M −N i =1
i i
∏ (z − p ) ∏ (z − p )
M M M
a0 D ( z ) a0
a .z
i =0
i
−i
i =1 i i =1 i
On appelle zéros, les valeurs de z pour lesquelles H(z)=0 et on appelle pôles, les valeurs de z pour
lesquelles H(z) est infini (annule le dénominateur). C'est ainsi que H(z) possède N zéros (zi), M pôles (pi). Si M>N,
elle possède (M-N) zéros en 0, sinon (N-M) pôles en 0.
Ainsi, la position de ses pôles et de ses zéros ( +le facteur d’amplitude K=b0/a0) va nous fournir une
description complète de H(z) (par conséquent de h(n) et H(f)) donc du comportement du système. H(z) peut donc
être représentée sous la forme d’un cercle modélisant la position des pôle set des zéros dans le plan complexe.
Exemple
1 Im(z)
3z − 2
H ( z) =
( z − 1)(z + 0.5) -1
X X
Un zéro en 2/3 et deux pôles p1= -0.5 et p2= 1 1 Re(z)
Remarques
-1
- Dans la plupart des systèmes, les ai et le bi sont réels les pôles et les zéros sont soient réels soient des paires
de complexes conjuguées.
- Rappelons que le rayon d'un système causal se trouve à l'extérieur d'un cercle. Par ailleurs, s'il est stable :
+∞
n
h ( n ) < ∞ , puisque H ( z ) =
h(n).z − n , il suffit donc que z=1 fasse partie de la RDC.
n = −∞
- Pour un système causal et stable, tous les pôles sont à l’intérieur du cercle unité (|pi|<1, ∀i). Le domaine de
convergence ne peut contenir de pôles puisque la TZ ne converge pas aux pôles. S’il est anti-causal, il sera stable
si les pôles sont à l’extérieur du cercle unité.
N
- Si le filtre est non-récursif H ( z ) = b .z
i =0
i
−i
. Un filtre RIF a tous ses pôles à l’origine et sera donc toujours stable.
15 0.8
0.6
10
0.4
200 0.2 1
5
0 0.9
150
0 -0.2 0.8
-5
-0.6 0.6
50
-0.8 0.5
-10
0 2 4 6 8 10 12 14 16 18
0
-1 0.4
0 2 4 6 8 10 12 14 16 18
0.3
-50
0.2
-100
0.1
-150
0 5 10 15 20 25 1 Im(z) 0
0 5 10 15 20 25
X X 300
X 250
200
1 X XX X X X 150
-1 1 100
0.8
0.6
Re(z) 50
0.4 X 0
0.2
X
0
X -100
-50
-0.2 0 2 4 6 8 10 12 14 16 18
-0.4
-0.6
-1
1
1
-0.8
1 0.9
-1 0.8
0 5 10 15 20 25 0.8
0.8
0.6
0.7
0.6
0.4 0.6
0.2 0.4
0.5
0 0.2 0.4
0.3
-0.2 0
0.2
-0.4
-0.2
0.1
-0.6
-0.4
0
0 5 10 15 20 25
-0.8
0 5 10 15 20 25
-0.6
0 2 4 6 8 10 12 14 16 18
- A un pôle pi simple ou multiple va correspondre une réponse impulsionnelle qui converge si |pi|<1. Elle
divergera dans le cas contraire, soit si |pi|>1 .
- Sachant qu'à chaque pôle complexe est associéun pôle conjugué cela donnera une réponse impulsionnelle h(n)
oscillante(cosinus ou sinus) amortiesi|pi=1,2|< 1 ou divergentesi |pi=1,2|> 1.
- Dans un système à phase minimale, tous les zéros sont à l’intérieur du cercle unité (|zi|<1, ∀i).
On suppose que le cercle unité (|z|=1) ∈ RDC de X (z). On restreint le calcul de X(z) au cercle unité en
posant z = e 2πj f Te.
Lorsqu'un zéro est placé sur un point donné du plan en z, la réponse fréquentielle sera de 0 au point considéré.
Un pôle quant à lui produira un pic au point correspondant. Plus les pôles ou les zéros sont proches du cercle
unité, plus ils influencent la réponse en fréquence [11].
– un zéro sur le cercle unité introduit une annulation du module pour la fréquence correspondant
– Un zéro au voisinage du cercle unité introduit une atténuation dans le module de la réponse en fréquence.
Atténuation d’autant plus importante que le zéro est proche du cercle unité.
– Un pôle sur le cercle unité introduit une résonance infinie dans le module de lar éponse en fréquence pour la
fréquence correspondante.
– Un pôle au voisinage du cercle unité introduit une résonance d’autant plus importante dans le module de la
réponse en fréquence que le pôle est proche du cercle unité.
Exemples: 1 pôle en 0 et un zéro en 0.7. Pour obtenir l'allure de H(f), on divise le vecteur du numérateur (en rouge)
sur celui du dénominateur (en mauve). Et pour le suivant un pôle 0.8 et un zéro en -0.5
f=fe /4
1.8
1
0.8 1.6
0.6
1.4
0.4
1.2
0.2
Imaginary Part
-0.2
x 1
0.8
-0.4
f=0 H(f)
0.6
-0.6
-0.8 0.4
-1
0.2
-1 -0.5 0 0.5 1
Real Part
0
0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5
8
7
0.8
0.6
6
0.4
5
0.2
Imaginary Part
f=fe 0
-0.2
x f=0 H(f) 4
3
-0.4
-0.6
2
-0.8
1
-1
-1 -0.5 0 0.5 1
Real Part 0
0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5
Autres Exemples
1)Sur la figure ci-dessous le cercle complet correspond à une fréquence d'échantillonnage fe. Des pôles proches
du cercle unité sont à l'origine de larges pics tandis que des zéros proches ou sur le cercle unité produisent des
minima. Ce tracé nous permettra d'identifier la nature du filtre.
f=fe /4
12
1
0.8 X
10
0.6
0.4
8
Partie Imaginaire
0.2
f=fe /2 0 6
f=0
-0.2
-0.4 4
-0.6
-0.8 2
X
-1
-1 -0.5 0 0.5 1 0
0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5
Partie Réelle
f=fe/4 f=fe/2
FEI, USTHB [assiakourgli@gmail.com / http://perso.usthb.dz/~akourgli/
51
Traitement Numérique du Signal Master ST M1 2019/2020
On peut aussi avoir une idée sur son comportement général : passe-bas, passe-hautou passe-bande, connaître sa
ou ses fréquences de coupure.
f=fe/4
2) 1 10
0.8 f=fe/6 9
0.6 X 8
0.4 7
Partie Imaginaire
0.2
2 60° 6
f=fe/2 0
2 f=0
5
-0.2
4
-0.4
3
-0.6
X 2
-0.8
1
-1
0
-1 -0.5 0 0.5 1 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5
z 2 − 2z + 1
3) H ( z) = Zéros double en z=-1, pôles p12 = ±0.6e j 72°
z − 0.371z + 0.36
2
- Un zéro double en z = 1 ⇒ |H(f)| = 0 pour f = 0 -Des pôles proches du cercle unité ⇒ maxima.
3
1
0.8
2.5
0.6
X
0.4
2
Partie Imaginaire
0.2
2
0
1.5
-0.2
-0.4 1
-0.6 X
-0.8 0.5
-1
-1 -0.5 0 0.5 1 0
Partie Réelle 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5
Remarque : Puisque les coefficients du filtre sont réels, les pôles et zéros sont réels (sur l'axe des réels) ou paires de
complexes conjugués.
La transformée en Z présente l’avantage d’être plus facilement inversible que la transformée de Fourier.
Le passage de la TZ vers h(n) peut se faire par le biais de la transformées en Z de signaux élémentaires connus à
condition qu'il soit possible d’écrire H(z) comme la combinaison de transformées élémentaires. Dans le cas
contraire, on peut employer l’intégration sur un contour fermé en utilisant le calcul des résidus, ou le
développement en puissance de z et de z-1, ou encore le développement en fractions élémentaires [12].
1. La relation générale de la transformée en z inverse est donnée par l’équation donnée par l'intégrale de Cauchy :
1
2.π . j
n −1
x(n) = X ( z ).z .dz , où C est un contour fermé parcouru dans le sens inverse des aiguilles d’une montre
C
contenant l’origine.
Re s[ z n−1 X ( z )] z = pi =
1 d m−1
(m − 1)! dz m−1
[
(z − pi )m z n−1 X ( z ) z = pi ]
pi=ea Re s[ z n −1 X ( z )] z = e = [z n ]z = e = e an .u ( n )
z
Exemple : X ( z ) = a
z − ea
a
2. Transformée inverse par division polynômiale :Il est possible de calculer la transformée en Z inverse selon les
puissances croissantes de z-1(système causal) ou selon les puissances décroissantes de z (système anti-causal).
X ( z ) = Cn z − n TZ
−1
Exemples
→x ( n) = Cn
n
( )
∞ ∞
1
= 1 + z −3 + z −6 + ............. h(n) = δ (n − 3k )
3
- y(n)=y(n-3)+x(n) H ( z ) = −3
= z −k
1− z k =0 k =0
1 1-a.z-1
-1+a.z-1 1 + a.z-1 + a2.z-2+ ......
0 +a.z-1
-a.z-1 + a2.z-2 x[0] x[1] x[2]
0 + a2.z-2 ......
1
On obtient : −1
= 1 + a.z −1 + a 2 .z − 2 x(n) = a n .u (n)
1 − a.z
z -a+z
-z+a-1.z2 -a-1.z –a-2.z2 –a-3.z3 -…………
-1 2
0 +a .z
-a-1.z2+a-2.z3 x[-1] x[-∞]
0 + a-2.z3
……
Notons que la division peut se réaliser sans faire apparaître une expression analytique générale.
3. L’idée générale de cette approche consiste à trouver pour une fonction X(z) complexe un développement en
fonctions en Z plus simples et pour lesquelles une transformée inverse est connue:
X ( z) = X i ( z) TZ
→x(n) = xi (n)
−1
i i
où les Xi(z) sont des fonctions dont les TZ-1 sont connues (Voir page 38).
Un filtre numérique est constitué d'un groupement de circuits logiques astreints à un processus de calcul
(ou algorithme) qui confère à ce filtre une fonction déterminée (passe-bas, passe-haut, passe-bande, réjecteur de
bande, intégrateur[y(n)=(x(n)+x(n-1))/2], différentiateur[y(n)=(x(n)-x(n-1))/2], ...). Il faut souligner que certains
filtres ne sont pas conçus pour arrêter une fréquence, mais pour modifier légèrement le gain à différentes
fréquences, comme les égaliseurs. Ce sont tous des systèmes linéaires, discrets, invariants dans le temps et
unidimensionnels. De plus, pour qu’ils soient physiquement réalisables, il faut qu’ils soient nécessairement
causaux.
1.4
Les filtres représentés ci-dessus sont idéaux. Dans un cas réel il n'est
|H(f)|
pas possible d'obtenir une fréquence de coupure aussi raide. Le 1+δ1 1.2
Remarque : Idéalement, il est souhaitable qu'un filtre possède une phase linéaire dans la bande passante. Une
phase linéaire assurera un même déphasage pour toutes les fréquences (pas de distorsion). Les filtres FIR peuvent
générer des filtres à phase linéaire à la condition que la réponse impulsionnelle soit symétrique.
H ( f ) = R( f )e − jφ ( f ) avec φ ( f ) = φ 0 + 2π fτ
Et la dérivée de cette dernière par rapport à f fournit le ‘retard de groupe ‘, défini donc par :
dφ ( f )
β =−
df
et qui correspond au retard subi par le signal après être passé par un filtre. Si la phase est linéaire (filtres RIF
symétrique), le retard est constant et le signal à la sortie aura donc une distorsion minimale puisque l’effet de la
phase sur le signal sera un simple décalage temporel (primordial dans un système audio)
Exemple y (n) =
1
(x(n) + 2 x(n − 1) + x(n − 2)) h(n) = 1 (δ (n) + 2δ (n − 1) + δ (n − 2) ) et H ( f ) = e −2 π j f cos 2 (πf )
4 4
1 0 2
0.9 1.8
-0.5
0.8 1.6
0.7 -1 1.4
Retard de Phase de H(f)
Module de H(f)
0.6 1.2
Phase de H(f)
-1.5
0.5 1
-2
0.4 0.8
0.2 0.4
-3
0.1 0.2
0 -3.5 0
0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5
fréquence (Hz) fréquence (Hz) fréquence (Hz)
RII ou RIF ?
Les filtres RII, on l'avantage qu'ils sont efficaces. Avec très peu de pôles et zéros on peut assurer la plupart
des réponses fréquentielles dont on peut avoir besoin dans les applications audio. Cependant, le filtre étant
récursif, les erreurs de précision numérique deviennent une question d'importance, car ils peuvent s'amplifier et
devenir hors de contrôle, d'abord dans la forme de bruit, mais éventuellement dans la forme d'instabilité. La
forme de la réponse impulsionnelle n'est pas facile à déterminer, non plus, car elle est définie indirectement par
les pôles et zéros de H(f).
Par contre, les filtres RIF n'ont jamais des problèmes d'instabilité, car la sortie n'est qu'une somme finie
d’échantillons de l'entrée. Cependant, quand la réponse impulsionnelle est longue, le numéro d'opérations peut
devenir un facteur décisif quand il faut choisir entre RIF ou RII. Un autre avantage des RIF est le retard de groupe
constant, qui permet d'avoir une distorsion de phase minimale sur le signal traité [13].
L’application d’un filtre numérique implique le calcul de la sortie y(n) à l’instant t=nTe à partir des sorties et
entrées précédentes plus la valeur courante de l’entrée.
M N
…
b0 b1 b2 .. bN
….. y(n)
Un filtre numérique est généralement constitué des éléments suivants :un ou plusieurs organes de retard
(ce sont des registres à décalage jouant le rôle de mémoires retardées), pilotés par une horloge de période; des
opérateurs arithmétiques (additionneurs et multiplieurs); des registres fournissant les coefficients de pondération
du filtre [18]. Ainsi, à chaque top d’horloge Te, les valeurs des registres subissent un décalage permettant de
calculer la nouvelle sortie. Cette structure peut être aussi réalisée par logiciel.
En employant la fonction de transfert H(z), diverses structures peuvent être utilisées : structure directe
(implémentation de l’équation aux différences), structure canonique (structure directe avec minimisation des
composants) et structure en éléments simples.
Série n°2
≠ 0 N1 ≤ n ≤ N 2
1.Une séquence finie x(n) est définie par : x ( n) = 1 0.5 0.25
= 0 ailleurs
3. Calculer la transformée en z, X(z), et esquisser la carte des pôles et zéros ainsi que la ROC pour chacune des
séquences suivantes :
x ( n ) = (0.5 ) u ( n ) + (0 .25 ) u ( n ) , y ( n ) = (0.25 ) u ( n ) + (0.5 ) u ( − n − 1) , z ( n ) = (0 .5 ) u ( n ) + (0 .25 ) u ( − n − 1)
n n n n n n
az − 1
5. Soit H(z) la fonction de transfert d’un SLIT causal avec : H ( z ) = avec a réel
z−a
Déterminer les valeurs de a pour lesquelles H(z) correspond à un système stable. Prendre une valeur de a=0.5.
Représenter alors les pôles et zéros de la fonction, la région de convergence. Donner et tracer|H(f)|.
7. Etablir les correspondances entre les diagrammes pôles zéros et les réponses enfréquence pour une fréquence
d’échantillonnage fe=1 en justifiant vos choix :
1 1 1 1
0.4
X
0.4 0.4 0.4
Partie Imaginaire
Partie Imaginaire
Partie Imaginaire
Partie Imaginaire
-1 -1 -1
X
-1
-1 -0.5 0 0.5 1 -1 -0.5 0 0.5 1 -1 -0.5 0 0.5 1 -1 -0.5 0 0.5 1
Partie Réelle Partie Réelle Partie Réelle Partie Réelle
7 2 2 8
1.8 1.8
6 7
1.6 1.6
6
5 1.4 1.4
1.2 5
4 1.2
1 1 4
3
0.8
0.8
3
2 0.6
0.6
2
0.4
0.4
1
0.2 1
0.2
0 0
0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0 0
0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5
8. On suppose que le tracé des pôles et des zéros de ce système est lesuivant :
- Est-ce un filtre RIF ou RII ? (Justifier votre réponse)
- Donner l’allure approximative de H(f)
- Déterminer H(z) puis déterminer et tracer h(n) 3
X
- A partir de h(n), étudier la stabilité, la causalité et l’invariance de ce filtre.
- Calculer et tracer H(f) pour au moins 3 valeurs
- Ce filtre possède-t-il un retard de groupe constant (justifier)
- Déterminer sa réponse pour une entrée échelon x(n)=U(n).
1
9. Trouver la séquence y(n) qui a comme transformée en z : Y ( z ) = .
6 − 5 z −1 + z − 2
2 2z2
10. De quelle fonction x(n), la fonction: X ( z ) = est la transformée en Z
z 2 − 0.8 2 z + 0.64
z −2 − 2z −1 + 1
11. Trouver h(n) correspondant à la transformée en z suivante : H ( z ) = 2 −2
a z − 2az −1 + 1
Solutions
4. (1-z-N)/(1-z-1)5. Pôle z=a, stable si |a|<1. Le module vaut toujours 1 (cellule passe-tout).
11
n
11
n
π
U (n) − U (n) 10. x(n) = 0.8 cos (n − 1)U (n)
n
9. y (n) =
22 33 4
n −1 n−2
11. x (n) = (n + 1)a U (n + 1) − 2na U (n) + (n − 1)a U (n − 1)
n
Exercices supplémentaires
1. Soit un système linéaire invariant dans le temps (SLID) dont la réponse impulsionnelle h(n) est telle que :
h(n)=1 pour 0 ≤ n ≤ 3 et 0 ailleurs
Calculer la réponse y(n) à la suite x(n) définie par :
- x(n)=anpour 0 ≤ n ≤ 5, avec a=0.7 et x(n)=0 ailleurs.
- x(n)=cos(2πn/8) pour 0≤ n ≤ 7 et x(n)=0 ailleurs.
Solution y(0) = 1 y(1) = 1+a = 1.7 y(2) = 1+a+a2 = 2.19 y(3) = 1+a+a2+a3 = 2.533 y(4) = a+a2+a3+a4 = 1.7731 y(5) = a2+a3+a4+a5 = 1.24117
y(6) = a3+a4+a5 = 0.75117 y(7) = a4+a5 = 0.40817 y(8) = a5 = 0.16807
y(0) = 1 y(1) = 1+cos(p/4) = 1.707 y(2) = 1.707 y(3) = 1 y(4) =-1 y(5) = -2.414 y(6) = -2.414 y(7) = -1 y(8) = 0 y(9) = 0.707 y(10) = 0.707
2. Soit y(n) = x(n) + ax(n-1)+ by(n-1), l’équation aux différences d’un système discret causal.
a) Trouvez h(n), la réponse impulsionnelle de ce système ; pour quelles valeurs de a et b le système est- il stable
b) Trouvez la réponse impulsionnelle du système formé par la mise en série de deux systèmes h(n).
c) Même question pour la mise en parallèle de deux systèmes h(n).
n n-1
Réponses :a) h(n) = b u(n) + ab u(n-1); stable pour a finie et |b|<1.
n n-1 2 n-2 n n-1
b) h(n)*h(n) = (n+1)b u(n) + (2nab + (n-1)a b )u(n-1) c) 2h(n) = 2b u(n) + 2ab u(n-1)
3.On considère un système linéaire régi par l’équation aux différences suivante :
y(n) =(x(n+m)+ x(n+m-1)+ x(n+m-2))/3 où m est un paramètre entier.
-Montrer que ce système est linéaire invariant dans le temps. Etudier la causalité et la stabilité
selon les valeurs de m.
-Calculer la fonction de transfert H(z) de ce système pour m=0 et m=1.
-En déduire la réponse fréquentielle.
Réponses Le système est causal si m<1 et toujours stable. Pour m=0,H(z)=1/3(1+z-1+z-2),
H(f) =0.33(1+ 2 cos(2πfTe)e-2 jπfTe (filtre RIF passe-bas à phase linéaire). Pour m=1, de même filtre moyenneur
5. Déterminer en utilisant la décomposition en éléments simples, la forme du signal x(n) dont la TZ est donnée
par :
z2 z2
X ( z) = 2 avec z >2 et X (z) = avec a <1
z − 3z + 2 z 2 − (a + 1) z + a
Réponses : x(n) = (2n+1 -1)u(n) x(n)= 1/(1-a)+ ana/(a-1)=(1-an+1)/(1-a)
z z
6. On considère la transformée : X ( z) = + on pose z0 = e (r + jθ ) , déterminer x(n).
z − z0 z − z * 0
Réponse x(n) = 2e rn cos(nθ ) U (n)
7. En utilisant la méthode des résidus puis celle de décompositions en éléments dans la TZ est connue, déterminer
z − z0
h(n) dont la transformée est : H ( z ) = z 0 : réel , p0 = ρe jθ
( z − p0 )( z − p0 )
*
ρ n −1 z0 z0
Réponse h ( n ) = 1 − cos θ sin( nθ ) + sin θ cos( n θ ) U ( n − 1)
sin θ ρ ρ
y ( n ) = x ( n − 1) − z 0 x ( n − 2 ) + 2 ρ cos θ y ( n − 1) − ρ y ( n − 2)
2
h(n) 0.3
1/3
On suppose que h(n) est donné comme ci-contre. 0.25
0.15
- Tracer son auto corrélation si l'on suppose qu'il est périodique de période 9.
0.1
- Etudier la causalité
0.05
- Est ce un filtre RIF ou RII
0
-A partir de l'expression de h(n), déduire le rôle de ce filtre : 1 2 3 4 5 6 7 8 9 10
9.On considère que l'équation aux récurrences du système suivant est donnéecomme suit :
y(n)=0.9 y(n-1) - 0.81 y(n-2) + x(n) + 2 x(n-1) + x(n-2)
- Etudier la causalité et l'invariance de ce système
- Est-ce un filtre RIF ou RII ?
- Déterminer H(z) et donner le tracé des pôles et zéros, en déduire le rôle de ce filtre :
- Donner les allures approximatives de h(n) et |H(f)|
- Déterminer h(n) et tracer la pour les 3 premières valeurs
- On suppose que fe= 6 kHz, quelle sera la sortie du filtre si l'on donne en entrée : un signal bruité par une
sinusoïde de 500 Hz puis un signal composé de 2 sinusoïdes l'une de 500 Hz et l'autre de 3000Hz
10. On considère le système LIT décrit par l'équation aux récurrences suivantes :
x(n)=δ(n)+2δ(n-1)+3δ(n-2)+4δ(n-3) et h(n)=9δ(n)+7δ(n-1)+4δ(n-2)+δ(n-3)
- Ce système est-il invariant ? Justifier.
- Etudier la stabilité et la causalité de x(n), h(n) et y(n).
- Les signaux x(n), h(n) et y(n) sont-ils à énergie finie (ou infinie) ? Justifier
- Déterminer y(n) directement et tracer le.
- Déterminer y(n) en passant par la TZ de x(n) et h(n)
- On considère que x(n) et h(n) sont périodiques de période 4, déterminer y(n) de 2 façons.
- Quelle est lien entre la convolution et l’autocorrélation ?
0.8
11. On suppose que le tracé des pôles et des zéros d'un système est le suivant :
0.6
- Donner les allures approximatives de h(n) et |H(f)| 0.4
0.9
puis en déduire le rôle du filtre
Imaginary Part
0.2 60°
2
0
- Déterminer H(z) (On supposera un gain de 1 en z= -1)
-0.2
- Déterminer les coefficients du filtre -0.4
- On suppose que fe= 3 kHz, quelle sera la sortie du filtre si l'on donne en entrée
-0.6 :
un signal bruité par une sinusoïde de 500 Hz : -0.8
-1
FEI, USTHB [assiakourgli@gmail.com / http://perso.usthb.dz/~akourgli/ -1 -0.5 0 0.5 1
Real Part
60
Traitement Numérique du Signal Master ST M1 2019/2020
Imaginaire
- Donner les allures approximatives de h(n) et |H(f)| 0
- Déterminer H(z) (gain de 1 à f=0)
- Déterminer les coefficients du filtre -0.5
- Esquisser la TF de la sortie Y(f) si la TF de l'entrée X(f) est la suivante:
1.4 -1
1.2 -1 -0.5 0 0.5 1
Réel
1
0.8
X(f)
0.6
0.4
0.2
0
0 0.1 0.2 0.3 0.4 0.5
fréquence
0.6
- Déterminer les pôles et les zéros 0.4
Partie Imaginaire
-0.2
x0.8
- En déduire le rôle du filtre : -0.4
0.2
-0.6
-1
-1 -0.5 0 0.5 1
Im(z) Im(z)
dont le RDC est donné comme suit :
RDC
- Pour R= 2
- Pour R=0.8 R
R
RDC Re(z) Re(z)
Solutions
7 à 15 voir interros et examens années précédentes
FEI, USTHB [assiakourgli@gmail.com / http://perso.usthb.dz/~akourgli/
61
Traitement Numérique du Signal Master ST M1 2019/2020
Rappel : Soit H(z) la transformée en z d’un filtre numérique donné dont la décomposition sous forme fraction
rationnelle est donnée par :
b0 + b1 z −1 + .............. +b N z − N
H ( z) =
1 + a1 z −1 + ................ + aM z −M
1. Déterminer h(n), les pôles et zéros et esquisser H(f). En déduire le rôle de H(f), puis calculer le retard de
groupe (en préparation).
2. Vérifier ces réponses par matlab.
3. Que peut-on dire sur la stabilité, la nature et le retard de groupe de ce filtre?
4. Si l'on remplace l'un des coefficients 0.5 par 1, que devient le retard de groupe?
5. Quelle serait la sortie d'un tel filtre si l'entrée était constante?
6. Rétablir les valeurs par défauts et rajouter les lignes concernant le fichier audio puis commenter.
7. Prendre la même portion du signal et observer sa TF avant et après filtrage en commentant.
8. Créer un signal composé de la somme de 2 sinusoïdes de fréquences 0.1 et 0.4 puis observer le signal
avant et après filtrage.
1 4
9. Refaire le même travail pour y(n)= x( n − i )
5 i =0
Rappel : Soit H(z) la transformée en z d’un filtre numérique donné dont la décomposition sous forme fraction
rationnelle est donnée par :
b0 + b1 z −1 + .............. +b N z − N
H ( z) =
1 + a1 z −1 + ................ + aM z −M
Grâce à la seule connaissance du vecteur b et du vecteur a, on peut analyser tout filtre et :
- Déterminer les pôles et les zéros du filtre (et étudier sa stabilité)
- Déterminer la réponde impulsionnelle ou indicielle
- Déterminer la réponse fréquentielle et le retard de groupe (dérivée de la phase), etc.
1. Déterminer h(n), les pôles et zéros et esquisser H(f). En déduire le rôle de H(f), puis calculer le retard de
groupe (en préparation).
2. Vérifier ces réponses par matlab.
3. Que peut-on dire sur la stabilité, la nature et le retard de groupe de ce filtre?
4. Si l'on remplace l'un des coefficients 0.5 par 1, que devient le retard de groupe?
5. Quelle serait la sortie d'un tel filtre si l'entrée était constante?
6. Rétablir les valeurs par défauts et rajouter les lignes concernant le fichier audio puis commenter.
7. Prendre la même portion du signal et observer sa TF avant et après filtrage en commentant.
8. Créer un signal composé de la somme de 2 sinusoïdes de fréquences 0.1 et 0.4 puis observer le signal
avant et après filtrage.
1 4
9. Refaire le même travail pour y(n)= x( n − i )
5 i =0
Le traitement numérique a conduit à une amélioration importante des dispositifs de filtrage linéaire
notamment en termes de fiabilité, de reproductibilité, de souplesse et de complexité des fonctions réalisables. En
outre, les filtres numériques ont aussi d'autres propriétés difficiles qu'il n'est pas aisé de mettre en œuvre dans le
cas des filtres analogiques, entre autres : le filtrage numérique en temps réel (transmissions numériques, codage
des sons MP3, synthèse de parole, télévision numérique par exemple).La plupart des modèles de filtres
analogiques peuvent ainsi être reproduits sous forme numérique. Les éléments physiques (résistance, capacité,
inductance, amplificateurs opérationnels) sont en quelque sorte transposés en éléments logiques [17].
La synthèse d’un filtre est un ensemble de processus qui débute par la définition des caractéristiques du
filtre, jusqu’à sa réalisation informatique et/ou électronique, en passant par la détermination de ses coefficients.
Pour synthétiser un filtre numérique, on considère connu le gabarit du filtre analogique et on cherche un système
numérique caractérisée par une fonction de transfert H(z) à insérer dans le circuit ci-dessus permettant de
satisfaire le gabarit analogique.
x(t) y(t)
Système analogique
Ha(p)
La détermination de la fonction de transfert d’un filtre numérique, par une méthode directe, n’est pas
toujours très simple. Par contre, le problème qui consiste à transformer un filtre analogique en un filtre numérique
est relativement simple. De ce fait, de nombreuses méthodes sont proposées pour concevoir un filtre numérique
à partir du filtre analogique équivalent. Dans tous les cas, la synthèse d’un filtre numérique est une approximation
d’un filtre analogique idéal équivalent. Il est nécessaire de contraindre un certain nombre de paramètres.
Les fonctions modèles utilisées pour la synthèse des filtres sont soit la réponse impulsionnelle soit la
réponse en fréquence (celle-ci est préférée) de filtres analogiques connus. Si l'on emploie la réponse
impulsionnelle, les éléments h(n) de la réponse impulsionnelle numérique sont obtenus en calculant h(t), la
réponse impulsionnelle du filtre analogique, aux instants t=nTe.
0.015
fait, il n'est plus possible d’obtenir un filtre passe-bas idéal (droit et 0.005
-0.005
filtres qui vont pouvoir être réellement synthétisés n'ont pas de
-0.01
réponse fréquentielle correspondant à la fonction porte, mais 0 5 10 15
t(s)
pourront s'en rapprocher.
2 10
1.8
8
1.6
6
Réponse impulsionnelle h(n)
1.4
1.2
4
H(f) idéal
0.8
TF 2
0.6
0
0.4
Troncature
-2
0.2
0 -4
-0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5 -10 -8 -6 -4 -2 0 2 4 6 8 10
frequence (Hz) t(s)
2 10
1.8
8
1.6
1.4 6
Réponse impulsionnelle h(n)
1.2
4
H(f) réel
0.8 2
0.6
0
0.4 TF -1
-2
0.2
0 -4
-0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5 -10 -8 -6 -4 -2 0 2 4 6 8 10
frequence (Hz) t(s)
Comme l'illustre la figure suivante, la troncature (multiplication par une porte de largeur NTe) du sinc dans
le domaine temporel se traduira par une convolution dans le domaine fréquentiel du filtre idéal avec un sinc
s'annulant tous les NTe. Pour de grandes valeurs de N, les sinc dans la bande passante se compenseront les uns
les autres mais autour des points de discontinuité (fréquence de coupure), les ondulations restent apparentes.
0.06
0.04
0.02
-0.02
-5 -4 -3 -2 -1 0 1 2 3 4 5
0.06
0.04
0.02
-0.02
-5 -4 -3 -2 -1 0 1 2 3 4 5
On peut observer que les différences vis-à-vis du filtre idéal (soit la fonction porte) sont principalement les
ondulations dans la bande passante et dans la bande atténuée ainsi que la largeur de la transition.
C'est ainsi que les spécifications du filtre vont être définies par un gabarit fréquentiel linéaire ou en dB
(décibels). Ce gabarit indique la ou les fréquences de coupure, la largeur de la bande de transition minimale
souhaitée, le maximum d’ondulation de la bande passante et de la bande atténuée, la fréquence d'échantillonnage
et éventuellement l'ordre maximal permis.
En pratique, plus les fréquences fa et fp sont proches, plus l'ordre du filtre devra être élevé. Pour un filtre idéal, ces
valeurs seraient confondues
L'emploi des filtres RIF peut se révéler attrayant eu égard à ses nombreux avantages : stabilité inconditionnelle
(Tous les pôles sont en 0), phase linéaire possible. Néanmoins, ils présentent l'inconvénient de nécessiter un plus
grand nombre de coefficients que les filtres RII pour obtenir les mêmes caractéristiques fréquentielles à cause de
l'absence de pôles hors 0. Ainsi, toute fonction de filtrage numérique stable et causale peut être approchée par la
fonction de transfert d'un filtre RIF.
Rappelons que la sortie d'un filtre RIF va s'exprimer comme une combinaison linéaire d'un ensemble fini
d'éléments d'entrée :
N −1 N −1 N −1
y(n) = bi x(n − i) d’où H ( z) = bn .z −n H ( f ) = bn .e −2π j f nTe
i =0 n=0 n =0
Ainsi, les coefficients de pondération ne sont rien d'autre que les valeurs de la réponse impulsionnelle du filtre.
Ces coefficients constituent les coefficients du développement en série de Fourier de la fonction de transfert H(f)
(voir TFTD chapitre 1)
Du fait qu'un filtre RIF possède une fonction de transfert polynomiale (non rationnelle), il ne peut être obtenu par
transposition d'un filtre continu. Les deux méthodes les plus utilisées pour l’approximation des filtres
numériques RIF sont alors:
o Développement par série de Fourier : cette série est ensuite tronquée par des fonctions fenêtres pour
limiter la réponse impulsionnelle. Les coefficients de Fourier coïncident avec les échantillons de la réponse
impulsionnelle du Filtre.
o Echantillonnage de la réponse fréquentielle : Cette méthode fait appel à la TFD. Celle-ci est appliquée aux
coefficients recherchés bipour obtenir une suite fréquentielle qui corresponde à la réponse fréquentielle
du filtre.
Il existe d'autres méthodes telles les méthodes d’optimisation qui sont basées sur la minimisation d’un critère
d’erreur entre la courbe réelle et le filtre idéal [9].
3. Méthode de la fenêtre
Cette technique consiste, connaissant l'expression analytique H(f) de la réponse fréquentielle continue (dont la
formulation mathématique connue)à approcher, à déterminer par utilisation de la transformée de Fourier à temps
discret inverse, la réponse impulsionnelle. Cette réponse temporelle non causale obtenue sera retardée pour la
rendre causale [18]. Ainsi :
1. A partir du gabarit idéal du filtre, on détermine les coefficients du filtre en limitant le calcul à N valeurs
réparties symétriquement autour de n=0. Puis, on calcule de la TFTD inverse du filtre idéal qui nous permettra
de retrouver les échantillons de la réponse impulsionnelle soient les coefficients du filtre :
1 fe / 2
H ( f )e 2 π j f n Te df N impair ( Filtre de type I )
f e − fe / 2
h( n ) = fe / 2
1 H ( f )e − j π f Te e 2 π j f n Te df N pair ( Filtre de type II )
f e − f/ 2
e
Remarques :
• Le filtre de type I est le plus employé (C’est celui que nous développons dans ce qui suit) hormis pour les
différenciateurs.
• Le filtre de type II possède un zéro en -1 (Nombre de zéros est impair donc répartition de zéros en paires
conjuguées + un zéro en fe/2 puisque une somme de cosinus) , il ne peut donc être employé pour un passe-
haut ou un coupe-bande.
• Il existe des filtres de type III (N impair) et IV (N pair) permettant d’obtenir une réponse impulsionnelle anti-
symétrique avec déphasage linéaire également (somme de sinus au lieu de cosinus). Ils possèdent tous deux
des zéros en 1 (sin en f=0 est nulle), leur utilisation n’est donc pas adaptée aux filtres passe-bas. En outre le
filtre de type III possède également un 0 en -1, on ne peut l’employer que pour un passe-bande.
2.Cette méthode produit une série infinie de coefficients, on limite, alors la réponse impulsionnelle à N
échantillons (troncature). Sachant que la troncature induit des ondulations, on peut faire appel aux fenêtres de
pondération pour les atténuer. Ainsi, la réponse impulsionnelle idéale h(n) sera multipliée par la fenêtre discrète
wN(n) de longueur N:
3. Il ne reste plus qu'à décaler la réponse impulsionnelle h'(n) pour avoir une solution causale.
Remarque : Pour le choix de fc il faudra faire attention aux fréquences de coupure à prendre en compte. Afin d'avoir
de bons résultats lors de la synthèse, ce ne sont pas les fréquences de coupure du filtre idéal qu'il faut utiliser mais
il faut déplacer celles-ci afin de les centrer dans la zone de transition. Pour un passe-bas, l'augmenter de la demi
zone de transition (∆f), soit fp+∆f et pour un passe-haut, la diminuer de la demi zone de transition. Pour un passe-
bande, diminuer la première fréquence de coupure de la demi zone de transition et augmenter la seconde de la
demi zone de transition, pour un rejecteur de bande, on fera l'inverse.
|H(f)|
Exemple :
fe / 2
1 2 π j f nT e
h(n) =
fe H ( f )e
− fe / 2
df
]
fc
1 1
e 2π j f nTe df = e2π j f nTe
fc
On pose fc=(fp+fa)/2 h(n) =
f e − fc 2πjnfeTe
− fc
On normalise les fréquences par rapport à fe/2 (ou fe) c.a.d que l’on remplace partout fc par fc/(fe/2), on obtient
alors :
Tout aussi facilement, on peut déterminer les réponses impulsionnelles d’un passe-haut (1-H(f)), d’un passe-
bande (différence de 2 passe-bas) et d’un coupe-bande donnés.
Les valeurs de la réponse impulsionnelle idéale h(n) sont données au tableau suivant. Les fréquences fc
indiquées dans ce tableau (fréquences de coupure désirées) s'expriment également en fréquences normalisées
(divisées parfe/2).
h'N(n)=h(n).w(n)
N −1
1 si n ≤ sin( Nπf )
si w(n) = 2 W( f ) = H 'N ( f ) = H ( f ) *W ( f )
0 ailleurs sin(πf )
La troncature temporelle introduit des ondulations et induit une zone de transition moins rapide déterminée
par la largeur du lobe principal. Un compromis est à faire entre la raideur et l'amplitude des ondulations. Notons
que cette méthode donne des ondulations de même amplitude dans la bande passante et dans la bande atténuée.
1.4
Passe-bande avec Fen Rectangulaire (N=51)
Passe-bande avec Fen Hanning (N=51)
1.2
Passe-bande avec Fen Rectangulaire (N=91)
Remarque: Si N augmente, l'étendue des oscillations diminue
(la réponse est plus plate) et la largeur de la bande de 1
grande.
Pour le choix de la fenêtre de pondération, on procédera comme suit : En fonction de l'atténuation δ2 requise
dans la spécification du filtre, on choisira le type de fenêtre w(n) à utiliser. Puis en fonction de la largeur de la
zone de transition ∆f (spécifiée aussi au départ) et du type de la fenêtre w(n),on déterminera la longueur de la
réponse impulsionnelle N.
Exemple : On veut synthétiser un filtre passe-bas de fréquence de coupure fc = fe/10 avec ∆f =fe/5 et une
1.5
ondulation en bande atténuée > 50 db (voir TP n°3)
sin( π n / 5 )
b- h ( n ) = 0.2
0.5
πn / 5
Imaginary Part
16
16
0 X
c- On choisit w(n) comme étant la fenêtre de Hamming
-0.5
-1.5
-1.5 -1 -0.5 0 0.5 1 1.5
Real Part
Traitement Numérique du Signal Master ST M1 2019/2020
sin(π n / 5 )
e- on calcule les valeurs h(n)= 0.2 pour -8≤n≤8
πn / 5
N -8 -7 -6 -5 -4 -3 -2 -1 0
h(n) -0,0399 -0,0456 -0,0329 0 0,0493 0,1064 0,1597 0,1974 0,2110
N 1 2 3 4 5 6 7 8
h(n) 0,1974 0,1597 0,1064 0,0493 0 -0,0329 -0,0456 -0,0399 0.3
Réponse impusionnelle h(n)
0.25 Réponse impusionnelle h'(n)
0.2
f- On multiplie h(n) par w(n) pour trouver h’N(n)=h(n).w(n)
0.15
-0.05
0 5 10 15 20
N 0 1 2 3 4 5 6 7 8
h(n) -0,0399 -0,0456 -0,0329 0 0,0493 0,1064 0,1597 0,1974 0,2110
h'N(n) -0,0031 -0,0050 -0,0067 0 0,0255 0,0730 0,1325 0,1826 0,2023
N 9 10 11 12 13 14 15 16
h(n) 0,1974 0,1597 0,1064 0,0493 0 -0,0329 -0,0456 -0,0399
h'N(n) 0,1826 0,1325 0,0730 0,0255 0 -0,006 -0,005 -0,0030
1.4
1.5
Passe-bas avec Fen Rectangulaire (N=17)
Passe-bas avec Fen Hamming (N=17)
1.2
Passe-bas idéal
1
1
0.5
0.8
Imaginary Part
16
16
0
X
0.6
-0.5
0.4
-1
0.2
-1.5
0
0 200 400 600 800 1000 1200 1400 1600 1800 2000 -1 -0.5 0 0.5 1 1.5 2 2.5
Real Part
fc1=400 fe=4000Hz
On peut remarquer l'emplacement de zéros autour du premier 0 en 1 permet d'atténuer les lobes secondaire et de
maintenir une réponse cste autour du zéro.
FEI, USTHB [assiakourgli@gmail.com / http://perso.usthb.dz/~akourgli/
73
Traitement Numérique du Signal Master ST M1 2019/2020
La méthode de synthèse par échantillonnage en fréquence est appliquée depuis la réponse fréquentielle d’un
filtre continu idéal H(f) dont on ne connait pas la formule mathématique (on ne peut alors calculer h(n) par TF
inverse de H(f)). On utilise alors la transformation de Fourier Discrète inverse. C'est-à-dire que l'on
"échantillonne" la réponse désirée dans le domaine fréquentiel, on obtient N points de cette réponse fréquentielle
auxquels on fait correspondre N points de la réponse temporelle équivalente obtenus par TFD inverse [18] comme
suit :
H (k ) = H ( f ) f =k / N k = −( N − 1) / 2 à ( N − 1) / 2
1 ( N −1) / 2
puis on applique la TFD inverse h(n) = H (k )e 2π j k n / N
N k =−( N −1) / 2
Cette méthode de synthèse est très simple et permet de réaliser toute forme de filtre (chose qu'on ne peut
réaliser avec la méthode précédente). Cependant, cette méthode de synthèse ne garantit que les points
fréquentiels H(k). Entre ces points, la valeur de H(f) n'est pas maitrisée, il peut y avoir des oscillations qui ne sont
pas également réparties avec un maximum d’erreur entre la réponse idéale et la réponse obtenue se situant autour
de la bande de transition.
Pour obtenir la réponse en fréquence du filtre finalement obtenu, on peut par exemple appliquer une TFD à la
réponse impulsionnelle h(n) de taille N obtenue, après avoir ajouté un grand nombre de zéros. Par ailleurs, du
fait de l'emploi d'une TFD inverse sur N points, la réponse impulsionnelle h(n) obtenue est périodique de période
N bien que la réponse impulsionnelle idéale souhaitée ne soit pas de durée limitée.
Exemple : On cherche à réaliser sous forme numérique un filtre passe-bas idéal de fréquence de coupure fc=fe/10
avec ∆f<fe/16. On prend donc N=17, ce qui nous donne ∆f =0.0588 (voir TP n°3)
On a N = 17, H(0)=H(−1)=H(1)=1 et H(2)= H(−2)=...= H(8)= H(−8)=0. On peut en déduire, par transformation de
0.2
Fourier discrète (TFD), les valeurs de la réponse impulsionnelle Réponse impusionnelle h(n)
0.15
1 ( N −1) / 2
H (k )e 2π j k n / N h( n) = (1 + e − 2πjn / 17 + e 2πjn / 17 )
1
h(n) =
N k =−( N −1) / 2 17
0.1
0.05
0.8 1
0.6
0.4 0.8
H(0) H(1)
0.2 16
Imaginary Part
16 0.6
0
X
-0.2
-0.4
0.4
-0.6
-0.8
0.2
-1
Sachant que H(k)*=H(-k) pour un signal h(n) réel, on peut, de manière générale, démontrer que :
2πn
( N −1) / 2
1 pour –(N-1)/2≤ n ≤(N-1)/2
h( n) = H (0) + 2 H (k ) cos
N k =1 N
La réduction des oscillations peut aussi s'obtenir à travers le fenêtrage. Ci-dessous les coefficients du filtre h(n)
suivis de ceux obtenus après fenêtrage de Hamming h'N(n).
n 0 1 2 3 4 5 6 7 8
h(n) -0,0257 -0,0269 -0,0153 0,0114 0,0514 0,0985 0,1430 0,1749 0,1865
h'N(n) -0,0020 -0,0031 -0,0033 0,0041 0,0279 0,0705 0,1237 0,1688 0,1865
n 9 10 11 12 13 14 15 16
h(n) 0,1749 0,1430 0,0985 0,0514 0,01139 -0,0153 -0,0269 -0,0257
h'N(n) 0,1688 0,1237 0,0705 0,0279 0,0041 -0,0032 -0,0031 -0,0020
1.5 1.4
Passe-bas avec Méthode Echantillonnage (Rect)
Passe-bas avec Méthode Echantillonnage (Hamm)
1.2 Passe-bas idéal
1
1
0.5
Imaginary Part
0.8
16
0 X
0.6
-0.5
0.4
-1
0.2
-1.5
Remarques: Pour le choix de N, il faut veiller à ce que fe/N soit inférieur à ∆f.
Par ailleurs, on peut quelques peu atténuerles ondulations en adoucissant les transitions dans la transmittance
du filtre. Pour cela, on introduira 0.5 entre 1 et 0, Mais de ce fait ∆faugmentera, on ajustera, alors, la valeur de N
en conséquence (2∆f<fe/16). Ce qui nous fournira une valeur de N de 33 et ∆f =0.0303. On prendra H(3)=
H(−3)=0.5. Ainsi H(0)=H(−1)=H(1)=H(2)=H(−2)=1 et H(3)=H(−3)=0.5et H(4)=H(-4)=.......= H(32)= H(−32)=0.
h ( n) =
1
(1 + 2 cos(2πn / 33) + 2 cos(4πn / 33) + cos(6πn / 33) pour -16 ≤ n ≤ 16
33
1.4
Passe-bas avec Méthode Echantillonnage F=17
Passe-bas avec Méthode Echantillonnage F=33
0.3 1.2 Passe-bas avec Méthode Echantillonnage F=33 + H(4)=0.5
Réponse impusionnelle h(n) pour N=17
0.25 Réponse impusionnelle h(n)pour N=33 1
0.2
0.8 H(0) H(1 ) H(2)
0.15
0.6
0.1
H(3)
0.4
0.05
0.2
0
-0.05 0
-20 -10 0 10 20 0 500 1000 1500 2000
Le gabarit du filtre est supposé connu. On cherche les N coefficients de la réponse impulsionnelle du
filtre RIFdont la réponse en fréquence H(f) entre dans le Gabarit :h(k) pour k = 0,..., N -1
Le principe consiste à mettre en œuvre une méthode itérative qui cherche à minimiser à chaque étape la «
distance » du filtre au gabarit. Puis, à chaque étape on détermine et on calcule la réponse en fréquence par
Transformée de Fourier Discrète. Puis, on utilise la technique du Zéro-Padding pour sur-échantillonner la
réponse en fréquence et vérifier précisément la distance au gabarit.
1.4 1.4
|H(f)|
1.2 1.2
1+δ1
1 1
0.8 0.8
Amplitude
Amplitude
0.6 0.6
0.4 0.4
0.2 0.2
0 0
0 0.2 0.4 0.6 0.8 1 1.2 1.4 0 0.2 0.4 0.6 0.8 1 1.2 1.4
Bande passante Bande de transition Bande atténuée Bande passante Bande de transition Bande atténuée
Les points qui n’entrent pas dans le Gabarit sont déplacés pour obtenir une nouvelle fonction de transfert
qui vérifie le Gabarit. Par transformée de Fourier inverse on obtient la réponse impulsionnelle de durée k*N
points d’un filtre qui entre dans le gabarit. On la tronque à N points pour obtenir un nouveau filtre et on vérifie
qu’il entre dans le Gabarit en calculant la distance au Gabarit qui est la moyenne des écarts, on la désire la plus
petite possible [9] ;
La réalisation concrète d'un filtre numérique consistera en fait à matérialiser l'algorithme de calcul pour la
structure retenue. On aura la possibilité de travailler : Soit en logique câblée (assemblage d'organes logiques, tels
que portes, mémoires, etc ...),soit en logique programmée (organisation autour d'un processeur de traitement du
signal (DSP) ou, même, utilisation d'un microprocesseur(micro-ordinateur) standard).
Série n° 3
1. Obtenir les coefficients d'un filtre à RIF passe-bas par la méthode de fenêtrage pour obtenir les spécifications
suivantes :
- Fréquence de coupure idéale : fc=1.75 kHz
- Largeur de transition : ∆f=0.5 kHz
- Atténuation en bande atténuée : A=-20 log10(δ)> 51 dB avec δ =min(δ1,δ2)
- Fréquence d'échantillonnage : fe=8 kHz
2. On souhaite approcher un filtre idéal passe-haut par un filtre à réponse impulsionnelle finie, synthétisé par la
méthode du fenêtrage. Ce filtre doit répondre aux spécifications suivantes :
Fréquence de coupure fc=2 kHz
Largeur de transition : ∆f=0.5 kHz
Atténuation en bande atténuée : A=-20 log10(δ)> 40 dB avec δ =min(δ1,δ2)
Fréquence d'échantillonnage : fe=8 kHz
- Déterminer l’expression mathématique exacte de h’(n).
- Calculer h’(0) et tracer approximativement h’(n).
- Quel est l’intérêt du fenêtrage et quel est son inconvénient ?
- Quel est l’inconvénient de cette technique de synthèse des filtres?
- Citer un avantage et un inconvénient de la synthèse par des filtres RIF de même que par les RII.
3. Soit un filtre passe-bas decoupure 0.25 et ∆f<0.08 (fe=1 ). On souhaite synthétiser un filtre d’ordre par la
méthode de l’échantillonnage fréquentielle.
4. Un filtre numérique est à implanter. La fréquence d'échantillonnage est fixée à 8kHz. Idéalement, on veut
supprimer la bande fréquentielle de 1kHz à 2 kHz.
- Tracez, pour ce filtre idéal, le module de H(f)
-Si on utilise la méthode d'échantillonnage en fréquence en prenant N points sur cette réponse idéale,
quel serait h(n) pour N=21.
- Quelle est l'inconvénient de cette méthode? comment y remédier
5. On suppose le filtre h(n) dont le module du spectre H(f) est donné ci-dessous :
On veut déterminer le filtre numérique h(n) équivalent par la méthode de l'échantillonnage fréquentiel. Le filtre
doit répondre aux spécifications suivantesfc1=fe / 5 ,fc2=2 fc1et une largeur de transition ∆f<fe / 20.Prendre
fe=4,2Khz
- Déterminer N, ∆fpuis H(k) et donner son tracé
- Déterminer h(n)
-Si on utilise la méthode des fenêtres avec une atténuation en bande atténuée Aa>40 db et une zone de transition
de ∆f=fe/7 donner l’expression de h(n) exacte (avec application numérique)
FEI, USTHB [assiakourgli@gmail.com / http://perso.usthb.dz/~akourgli/
78
Traitement Numérique du Signal Master ST M1 2019/2020
Solutions
Les coefficients du filtre sont symétriques, il suffira par conséquent de calculer les valeur de h(0)à h(26).
sin(π n f c )
2. ∆f=0.5 fc= 0.125 h(n) = − f c pour k≠0 et h(0)=1-fc A>40 Hanning N=51
πnfc
sin(π n / 2 ) 2πn
h' (n) = −0.5 (0.5 + 0.5 cos( )) h' (0) = 0.5 -25≤n≤25
πn / 2 50
3. fc =0.25, fe=1 et N=13 ainsi : (fe/N=0.07, 2fe/N =0.15, 3fe=0.23, 4fe/N=0.30) H(0)=H(1)=H(-1)=H(2)=H(-2)=H(3)=H(-3)=1
et les autres à 0
1 ( N −1) / 2
2πkn 1 6
2πkn
N Π
h( n ) = H (0) + 2 H (k ) cos( ) = H (0) + 2 H (k ) cos( ) ( n)
N k =1 N 13 k =1 N
4. fc1 =1, fc2 =2 fe=8 et N=21 ainsi : (fe/N=381, 2fe/N =762, 3fe=1143, 4fe/N=1523, 5fe/N=1905, 6fe/N=2286) H(3)=H(-
3)=H(4)=H(-4)=H(5)=H(-5)=0 et les autres à 1.
1 ( N −1) / 2
2πkn 1 2
2πkn 10
2πkn
N Π
h(n) = H (0) + 2 H (k ) cos( ) = H (0) + 2 H (k ) cos( ) + 2 H (k ) cos( ) ( n)
N k =1 N 21 k =1 N k =6 N
La TF de ce H(n), va effectivement passer par ces points mais en dehors de ces points le comportement du filtre peut dévier
du filtre idéal. Pour améliorer, il faudra prendre plus de points, d’où un filtre plus long.
Hanning, N=23, h(n)=0.8 sinc(0.8n)- 0.4 sinc(0.4n) et ℎ = 0.8 0.8 0.4 0.4 0.5 0.5cos 2 /22 Π n
h’(0)=0.4
6. y(n)=b3(x(n-3)+x(n-4))+b2(x(n-2)+x(n-5))+b1(x(n-1)+x(n-6))+b0(x(n)+x(n-7))
Type II , pas adéquat pour un passe-haut (un zéro en -1)
Exercices supplémentaires
1. Au cours de la transmission d'un signal numérique (échantillonné à une fréquence de 2,5 kHz) , il a été affecté
d'un bruit localisé entre les bandes de fréquence 350 Hz et 550 Hz. On veut éliminer le bruit par l'emploi d'un
filtre RIF possédant une possédant une bande de transition ∆f=100 Hz. Concevoir ce filtre :
A] Par la méthode du fenêtrage. On souhaite une atténuation en bande atténuée A=-20 log10(δ)> 20 dB.
- Tracer H(f) idéal
- Déterminer h(n) et l'ordre N du filtre
- Calculer h(0), h(1)=h(-1)
- Tracer alors approximativement H(f)
2. Au cours de la transmission d'un signal numérique (échantillonné à une fréquence de 5 kHz) , il a été affecté
par un bruit sinusoïdal de fréquence f0=250 Hz. On veut éliminer le bruit par l'emploi d'un filtre possédant une
bande de transition ∆f=±50 Hz à -3db. Concevoir un filtre RIF par la méthode du fenêtrage. On souhaite une
atténuation en bande atténuée A=-20 log10(δ)> 40 dB avec fc1,2=f0± ∆f
-fc fc
1. Les tracés du module de la réponse en fréquence de 0 à fe/2 et des pôles et zéros d’un filtre numérique
synthétisé par la méthode de l’échantillonnage fréquentiel sont donnés comme suit :
1.4
6
1.2
4
1
Imaginary Part 2
0.8
14
0
0.6
-2
0.4
-4
0.2
-6
0
0 500 1000 1500 2000 0 2 4
Real Part
But du TP : Dans ce TP, on teste deux méthode différents pour synthétiser un filtre RIF: méthode de séries de
Fourier (ou méthode des fenêtres) et méthode de la TFD dite par échantillonnage fréquentielle.
I. Rappels
Un filtre de réponse impulsionnelle finie (RIF) possède une fonction de transfert polynomiale. Il ne peut
pas être obtenu par transposition d’un filtre continu, comme cela est fait pour les filtres RII. Les filtres RIF
présentent l’inconvénient de nécessiter un grand nombre de coefficients pour obtenir les mêmes caractéristiques
fréquentielles. Mais par contre, ils sont inconditionnellement stables. On peut synthétiser des filtres RIF à phase
linéaire, c’est-à-dire à temps de propagation de groupe constant.
1.4
Le gabarit d’un filtre n’est autre que l’ensemble des caractéristiques 1+δ1
du filtre, à savoir : 1
1-δ1
- Le gain du filtre dans la bande passante. 0.8
Amplitude
- L’atténuation du filtre en bande coupée fa.
- La fréquence de coupure fc, on l’exprime souvent sous forme 0.6
normalisée par rapport à la fréquence d’échantillonnage. 0.4
- La largeur de bande de transition ∆fsouhaitée qui doit être la plus petite fpfa
possible. δ2
0.2
- La détermination des coefficients d’un filtre RIF par la méthode de la fenêtre est réalisée par la fonction Matlab
FIR1. Pour utiliser la technique d'échantillonnage de la réponse fréquentielle, on emploiera la fonction Matlab
FIR2.
Le principe de définition du cahier des charges d'un filtre récursif (RII) se déroule de la même manière que
pour un filtre non récursif (RIF) mais la synthèse se fera de différente manière. L'intérêt d'employer un filtre RII
réside principalement dans la possibilité d'obtenir une bande de transition étroite pour un ordre raisonnable bien
qu'il présente d'une part, un risque d'instabilité du à une grande sensibilité numérique des coefficients mais que
l'on peut toutefois contrôler en déterminant une structure mieux adaptée, et d'autre part, une variation de phase
fortement non linéaire.
Rappelons que pour un filtre RII, la sortie s'exprime comme une combinaison linéaire d'un ensemble fini
d'éléments d'entrées et de sortie: N
N M
−i
b .z i
y ( n ) = bi x ( n − i ) − ai y ( n − i ) d’où H ( z) = i =0
M
i =0 i =1 1 + ai . z − i
i =1
La méthode directe consiste à placer des pôles et des zéros aux fréquences utiles mais la plus courante (indirecte
dite aussi de transposition) est l’utilisation des méthodes de synthèse des filtres analogiques aboutissant à une
fonction H(p) correspondant aux spécifications. Une fonction permettant le passage du plan p dans le plan z (p =
fct(z)) est ensuite utilisée pour obtenir H(z). Cette fonction doit maintenir la stabilité du filtre analogique et
maintenir, au mieux, les caractéristiques de la réponse fréquentielle H(f) du filtre numérique.
De ce fait, on peut, par un placement adéquat des pôles et des zéros, obtenir un filtre sélectif simple.
o A partir des spécifications du filtre, placer des pôles et des zéros sur le cercle (ou à partir des zéros et des
pôles de H(p) calculer les zéros et les pôles de H(z)).
FEI, USTHB [assiakourgli@gmail.com / http://perso.usthb.dz/~akourgli/
84
Traitement Numérique du Signal Master ST M1 2019/2020
o Si c’est un passe-bas, le caractère passe-bas est renforcé en complétant l’ensemble des zéros par autant de
zéros que possible en z=-1 tout en conservant à H(z) son caractère propre (degré du numérateur égal au degré
du dénominateur). Si c’est un passe-haut, placer les zéros en z=1.
o On ajuste le coefficient K pour obtenir un gain unitaire à z= 1 pour un passe-bas, à z=-1 pour un passe- haut.
X
•
0.8
0.4 r
Imaginary Part
0.2
2
K ( z + 1)( z + 1) K ( z + 2 z + 1)2 θ 2
H ( z) = = 2
0
( z − re )( z − re ) z − 2 zr cos θ + r 2
− jθ jθ -0.2
-0.4
1 + r − 2r cos θ
-0.6
∆f
2
θ = ( f c / f e )× 360 °
X
K= r = 1 − 3 db π
-0.8
avec et -1
4 fe -1 -0.5 0
Real Part
0.5 1
0.8 X
0.6
K ( z − 1)( z − 1) K ( z 2 − 2 z + 1) 0.4 r
H ( z) = = 2
Imaginary Part
0.2
( z − re − jθ )( z − re jθ ) z 2 − 2zr cosθ + r 2 0
-0.2
θ 2
1 + r 2 + 2r cosθ -0.4
avec K= -0.6
X
4 -0.8
-1
-1 -0.5 0 0.5 1
Real Part
K ( z − 1)( z + 1) K ( z 2 − 1) 0.8 X
H ( z) = = 2 0.6
( z − re − jθ )( z − re jθ ) z − 2 zr cosθ + r
2 0.4
r
Imaginary Part
0.2
θ 2
2 jθ jθ
− 2r cosθ e + r
0
2
e -0.2
avec K= -0.4
e 2 jθ − 1 -0.6
-0.8 X
-1
-1 -0.5 0 0.5 1
K ( z − e )( z − e ) K ( z − 2 cos θ z + 1)
− jθ jθ 2 0.8
X
H ( z) = = 2 0.6
0.2
θ
1 + r 2 − 2r cos θ '
2
0
avec K= -0.2
2 − 2 cos θ -0.4
-0.6
-0.8 X
• Pour un filtre passe-bas de premier ordre, nous avons : -1
-1 -0.5 0 0.5 1
Real Part
K ( z / α + 1/ α ) 1
z − (1 − K ) / α 0.6
0.4
Partie Imaginaire
0.2
0 α X
-0.2
K (z / α −1/ α )
-0.6
H ( z) = -0.8
z − (1 − K ) / α
-1
-1 -0.5 0 0.5 1
Partie Réelle
Exemple: Trouver, grâce à la position des pôles et des zéros, la fonction de 1.4
transfert et l’équation de récurrence d’un filtre numérique simple qui a les 1.2
caractéristiques suivantes : 1
0.2
Pour rejeter les composants à 50Hz, on place une paire de zéros aux 0
0 50 100 150 200 250
fc=50 fe=500Hz
points du cercle unité qui correspondent à 50Hz. Ces points se situent à 1
0.4
de coupure, deux pôles complexes conjugués sont placés sur le cercle de
Partie Imaginaire
0.2
-0.2
La relation entre la largeur de bande et le rayon est donc applicable, par -0.6
-0.8
conséquent, la norme des pôles est 0,937.
-1
-1 -0.5 0 0.5 1
Depuis la figure, la fonction de transfert du filtre est donnée par : Partie Réelle
y(n)≈x(n)-1.618x(n-1)+x(n-2)+1.5161y(n-1)-0.878y(n-2)
Pour concevoir un filtre analogique, on peut employer des filtres passifs obtenus par combinaison de
résistances, de condensateurs et/ou de bobines ou encore utiliser des filtres actifs comportant un élément
amplificateur (transistor, AO, etc.) qui permet donc de modifier les amplitudes des signaux. Ainsi, un filtre passif
de base (de premier ordre) peut être composé d'une cellule RC ou d'une cellule RL.
FEI, USTHB [assiakourgli@gmail.com / http://perso.usthb.dz/~akourgli/
86
Traitement Numérique du Signal Master ST M1 2019/2020
- les filtres passifs (L,C, quartz, etc ;) sont utilisés pour les hautes fréquences
- les filtres actifs (R,C, ALI) utilisent des amplificateurs linéaires integrés (ALI) limités à 1Mhz
- les filtres à capacité commutés (R et C intégrés, ALI, Interrupteur commandé MOS) ont des fréquences
(<10 Mhz) programmables
Ci-dessous, sont donnés quelques exemples de structure pour des filtres passifs passe-bas de premier et second
ordre suivis de filtres actifs passe-bas de premier et second ordre aussi.
Pour obtenir des filtres d'ordre N plus élevé, on emploiera des filtres du premier et du deuxième ordremis en
cascade. Les pentes asymptotiques seront proportionnelles au nombre de cellules (+N ↑ et +∆f↓). A noter
également que les filtres passe-bas ou passe-haut peuvent avoir un nombre entier d'ordre (1, 2, 3...) tandis que les
filtres passe-bande ou coupe-bande ne peuvent qu'avoir un ordre pair (2, 4, 6, ...) car ils sont formés de paires de
cellules : 2 cellules RC ou une cellule RC et une cellule RL.
Pour synthétiser des filtres analogiques répondant à un gabarit donné, on choisira parmi un ensemble de filtres
analogiques testés et éprouvés donc connus pour leurs propriétés en termes de pente d’atténuation et
d’ondulation dans la bande passante et atténuée telles que :
Remarque : Le filtre de Bessel est intéressant pour filtrer les signaux large bande (modulations HF haut débit :
PSK, 8-PSK, OFDM…) en préservant les phases. Il n’a pas intérêt dans le domaine du filtrage numérique.
1. Un filtre de Butterworth est caractérisé par le fait que la réponse d’amplitude est maximalement plate dans
la bande passante et monotonement décroissante à partir d’une certaine fréquence (fréquence de coupure).
L’amplitude d’un filtre analogique Passe-bas de Butterworth d’ordre N est définie par l’expression :
H (ω ) = H (ω ) =
1 ou 1
2N 2N
ω ω
1 + 1 + ε .
2
ωc ωc
ωc est la pulsations de coupure limiteset ε estun paramètre de conception qui fixe la région de tolérance dans
la bande passanteδ1=1/(1+ε2)1/2.Le filtre analogique aura N pôles placés sur lecercle unité (à partie réelle
négative).
L’ordre N est déterminé par la zone de transition (ωa-ωp)et les ondulations permises en bande passante AP=20
log(1+δ1) et bande atténuée (Aa=-logδ2)
10 Aa Ap
1 1
Aa
log (10 − 1) /(10 − 1) log 2 − 1 / 2 − 1
10
log10 − 1
10
= δ2 δ1
N ≥=
N≥ 2 log(ω a / ω p ) 2 log(ω a / ω p )
2 log(ωa / ω p )
Remarque :
Plus N est élevé, plus grande est la sélectivité (zone de transition ∆f étroite) mais un filtre RII de fait de sa
récursivité a une période d’initialisation avant le régime permanent qui augmente avec l’ordre du filtre. En outre,
son déphasage augmente avec l’ordre du filtre.
Déterminer le filtre et N.
Filtre :Butterworth N ≥=
(
log (10 4 ) / 0.2589 )
= 4.80 ≈ 5
2 log(3 / 1)
2. Un filtre de Tchebychev est caractérisé généralement par une ondulation équilibrée dans la bande coupée.
La forme analytique du module de réponse fréquentielle d’un filtre analogique de Tchebychev (ou
Chebyshev) d’ordre N est donnée par :
H (ω ) =
1
ω
1 + ε 2 .T N2
ωc
L’ordre N est déterminé par la zone de transition (ωa-ωp) et les ondulations permises en bande passante AP=20
log(1+δ1) et bande atténuée (Aa=-logδ2) *
Aa
10
Ap
10 10 Aa 10
log 2.10 20 / 10 − 1
Aa Ap Ap
log 10 − 1
10 − 1 + 10 10 − 1
10 − 1 − 1
N≥ N≥
log ω a / ω p + (ω / ω p ) − 1
2
log ω a / ω p +
(ω a / ω p ) − 1
2
a
Exemple
Aa
10
Ap
log 2.10 20 . / 10 − 1
2.100 / 0.5089
N≥ = = 3.39 ≈ 4
log ω a / ω p +
(ω a /ωp )
2
− 1
3+ 8
Dans la pratique, on utilise trois valeurs d'ondulation, ε = 0.1 dB, 0.5 dB et 1 dB. Pour ces 3 valeurs HN(p) est
donnée :
On rappelle que les fonctions de transfert HN(p) des filtres analogiques polynômiaux (Butterworth,
Tchebychev,Bessel, Cauer, etc.) sont données pour des fréquences de coupure normalisées et uniquement pour
desfiltres passe-bas [2].
1 p ω
Pulsation centrale
Passe-bas
p = p /ωA Passe-bande p = + A ω A = ω A1 ω A 2
B ωA p
−1 Largeur de Bande
1 p ω A
p = ωA / p p = + B = (ω A2 − ω A1 )
Passe-haut Coupe-bande
B ωA p
A partir de l'expression de HN(p) normalisée, on dénormalise en remplaçant p par les valeurs données au
tableau précédent permettant d'aboutir à une fonction de transfert dénormalisée H(p).
La méthode de l’invariance impulsionnelle consiste à effectuer la synthèse d’un filtre numérique dont la
réponse impulsionnelle h(n) est l’échantillonnage de la réponse impulsionnelle h(t) du filtre analogique
équivalent [18] :
- on recherche la fonction de transfert H(z) du filtre numérique qui a pour réponse impulsionnelle la suite h(n).
h(n)
ha(t)
Rappelons toutefois que l'échantillonnage temporel se traduit par une périodisation en fréquence telle que :
−1
Ha ( p)
L =nTe
→ha (t) t→h(n) = Teha (nTe )
TZ
→H(z) ou H ( z ) = Te H a ( p)
pôles pi de Ha ( p )
Résidus −1 pT
1− z e e p = pi
Exemple 1
On veut synthétiser un filtre numérique qui possède la même réponse impulsionnelle qu'un filtre passe-bas
analogique du 1er ordre de transmittance normalisée: H N ( p ) = 1
(1 + p )
1 ωc
On commence par dénormaliser H a ( p ) = H N ( p ) p = p / ωc = =
1 + ( p / ωc ) ωc + p
1 z
h ( n ) = Teω c e − nTe ω c pour n≥0 H ( z ) = Teω c −Te ωc
= Teω c
1− e z −1
z − e −Teωc
1
Il s'avère qu'il y a égalité entre les 2 réponses à un coefficient K près, on posera donc : H ( z ) = KTeω c −Te ω c
1− e z −1
FEI, USTHB [assiakourgli@gmail.com / http://perso.usthb.dz/~akourgli/
91
Traitement Numérique du Signal Master ST M1 2019/2020
1 − e − Te / τ
Sachant quep=2πjf, pour f=0, Ha(0)=1 et z=e2πjfTe, pour f=0, H(0)=KTeωc/(1-e-Teωc ) H ( z) =
1 − e −Te / τ z −1
Remarque : La réponse du filtre numérique sera proche de celle du filtre analogique dans la bande [-fe/2, fe/2] si
le filtre analogique a une réponse fréquentielle nulle en dehors de cette bande. Cette méthode est utile seulement
dans le cas de filtres analogiques à bande limitée.
En effet, cette méthode n'est applicable que si l'on est en mesure de définir une fréquence d'échantillonnage
adéquate [19]. Pour cela, il faut que la réponse en fréquence du filtre continu soit nulle (ou presque) au-delà d'une
certaine fréquence de valeur finie. Cette méthode ne peut donc s'appliquer aux filtres de types passe haut ou
coupe-bande (fmaxest ∞ donc on ne peut définir fe>2fmax).
Exemple 2
On considère le filtre analogique d’ordre 3 associé à la fonction d’approximation de Butterworth. Réaliser le filtre
numérique passe-bas correspondant en utilisant la méthode de l’invariance impulsionnelle. La fréquence de
coupure est fc=1 kHz et la fréquence d’échantillonnage fe=10 kHz.
Rappelons que le filtre de Butterworth d’ordre 3 est formé de 2 cellules d’ordre 1 et 2, respectivement :
1 1
H N ( p) = H ( p) =
(
(1 + p) 1 + p + p )
2 alors
(
(1 + p / ωc ) 1 + p / ωc + ( p / ωc )2 )
Pour éviter de décomposer en éléments simple, on peut calculer H(z) en utilisant les relations entre
transformées de Laplace et en z pour les systèmes du premier et du second ordres.
H(p) Te.H(z)
1 1
Te _ aTe −1
p+a 1− e z
p+a 1 − e − aTe z −1 cos(ω0Te )
Te
( p + a )2 + ω0 2 1 − 2e −aTe z −1 cos(ω 0Te ) + e − 2aTe z − 2
ω0 e − aTe z −1 sin(ω 0Te )
Te
( p + a )2 + ω0 2 1 − 2e − aTe z −1 cos(ω 0Te ) + e − 2 aTe z − 2
ωc ωc
p+ −
ωc ωc p ωc
H ( p) = − = − ωc 2 2
(ω c + p ) ω ω
2 2
(ω c + p ) ω
2
ω2
p+ c +3 c
p+ c +3 c
2 4 2 4
3ω cωc
p+
ωc ωc
H ( p) = −ω 2 + 2
(ωc + p ) c ω c 2 ωc2 3 ωc
2
ω2
p+ +3 p+ +3 c
2 4 2 4
Il reste à déterminer le gain K en adaptant les gains des deux filtres analogique et numérique pour une
fréquence donnée. Comme il s'agit d'un filtre passe-bas, on supposera l'égalité pour f=0 (p=0 et z=1)
Cette méthode a pour objectif de faire coïncider au mieux les domaines analogique et numérique. Les filtres
qui en sont dérivés sont plus stables que sont obtenus à travers l'emploi de la méthode de la variance
impulsionnelle. Mais, en contrepartie, elle introduit une distorsion sur l'axe des fréquences.
On sait qu’un filtre analogique est caractérisé par sa fonction de transfert H(p) et un filtre numérique est
défini par sa fonction de transfert H(z) avec z= e2̟jfTe= epTe. Ce qui implique que p=Ln(z)/Te. Dans les deux cas, cette
transformation entraîne une relation non linéaire entre les fréquences fA du domaine analogique et les fréquences
fN du domaine numérique.
Pour conserver le caractère polynomialdes fonctions de transfert, une approximation de ln(z) par les séries
de Laurent est employée :
1 − z −1 1 1 − z −1 3 1 1 − z −1 5
ln( z ) ≈ 2 + + + .........
1 + z −1 3 1 + z −1 5 1 + z −1
Ainsi, on peut montrer qu’on obtient la correspondance entre fréquence analogique fA(ou pulsation
analogique ωA) et la fréquence numérique fN (ou ωN) par :
2 ω N Te 2 π f N
ωA = tg = tg
Te 2 Te f e
Cette équation montre qu'il n’y a pas égalité entre pulsation analogique
et pulsation discrète et que la relation les liant n’est pas non plus linéaire fe π fN
fA = tg
puisqu'il y a distorsion des fréquences, y compris du retard de groupe. π fe
Pour faire la synthèse d’un filtre numérique par la transformation bilinéaire, on procède comme suit :
- On définit les caractéristiques souhaitées du filtre numérique (fréquence d’échantillonnage, de coupure, etc.)
2 π fN
- On calcule les pulsations analogiques ωA correspondant aux pulsations numériques ω A = tg
Te f e
- On détermine le gabarit du filtre analogique HN(p) normalisée d'ordre n (Chebyshev, Butterworth, etc.) qui
servira de modèle au filtre numérique et on écrit la fonction de transfert dénormalisée H(p) de ce filtre analogique
(qu’il faut recalculer en fonction de ωA).
2 1 − z −1
- On applique la transformation bilinéaireà H(p) en remplaçant p = , ce qui donne la fonction H(z).
Te 1 + z −1
Exemple 1
On désire concevoir un filtre passe-bas numérique de premier ordre à partir d’une fonction de transfert d'un filtre
RC dans le domaine continu (HN(p)=1/(1+p)). La fréquence de coupure désirée est fN=30 Hz et la fréquence
d’échantillonnage est fe=150 Hz.
2
0.7265
2 π 2 1 Te
On calcule d'abord ω A = tg = 0.7265 puis H ( p ) = =
Te 5 Te (1 + p ) p= p / ω A 2
p + 0.7265
Te
2
0.7265
Te 0.7265( z + 1) 0.7265( z + 1)
H ( z ) = H ( p) p= 2 z −1 = = =
2 z −1 2
Te z +1
+ 0.7265 ( z − 1) + 0.7265( z + 1) 1.7265 z − 0.2735
Te z + 1 Te
2 π f
On peut remarquer que f A = tg = e 0.7265 = 34.69 ≠ f N = 30
Te (2π ) 5 π
Exemple 2 : On désire réaliser un filtre numérique passe-bas du second ordre avecles caractéristiques suivantes
[18]: fréquence de coupure fN= 500 Hz (gain de 1), fréquence d'échantillonnage fe=5 kHz,
2 ω N Te 2 π 2
ωA = tg = tg = 0.325
Te 2 Te 10 Te
2
0.035 p
0.1075 p 0.1075 pω A Te
H N ( p) = H ( p) = 2 =
p + 0.1075 p + 1
2
p + 0.1075 pω A + ω A
2
2 4
p 2 + 0.035 p + 2 0.1056
Te Te
2 2 z −1 z −1
0.035 0.035
Te Te z + 1 z +1
H ( z ) = H ( p) p= 2 z −1 = =
z −1 z −1
2 2
Te z +1 2 z −1 2 2 z −1 4
+ 0.035 + 2 0.1056 + 0.035 + 0.1056
Te z + 1 Te Te z + 1 Te z +1 z +1
Pour un filtre numérique passe-bas ou passe-haut, on peut déterminer H(z) directement à partir de HN(p) en
π fN
utilisant la pulsation normalisée Ω c = tg et remplacer p dans HN(p) comme suit:
fe
1 z −1
- pour un passe-bas par p = soit H ( z ) = H N ( p ) p = 1 z −1
Ωc z +1 Ω c z +1
z +1
- pour un passe-haut par p = Ω c soit H ( z ) = H N ( p ) p = Ω z +1
z −1 c
z −1
- pour un passe bande ou un rejecteur de bande, la détermination de ωA s'obtient par la pulsation centrale
ω A = ω A1ω A 2
Remarque: Puisqu'il n'est pas possible de faire coïncider précisément la réponse fréquentielle du filtre analogique
(sur[0, +∞[ et numérique (sur [0, fe/2] la transformation bilinéaire aura pour objet de faire coïncider au mieuxles
deux réponses. Cette méthode donne d'assez bons résultats à condition que les fréquences caractéristiques ne
soient pas trop proches de la demi-fréquence d'échantillonnage.
Série 4
Exercice 1 : Déterminer, grâce à la position des pôles et des zéros, la fonction de transfert, les coefficients du filtre
et l’équation de récurrence d’un filtre numérique passe-bande qui a les caractéristiques suivantes :
1
0.4
- Fréquence d’échantillonnage : 500Hz
Partie Imaginaire
0.2
Solution 0
-0.2
- Passe bande Placer un zéro en 1 et un zéro en -1 -0.4
K ( z − 1)( z + 1) K ( z 2 − 1)
Remplacer θ et r dans l’équation H ( z ) =
-1
- =
z 2 − 2 zr cos θ + r 2
-1 -0.5 0 0.5 1
( z − re − jθ )( z − re jθ ) Partie Réelle
K ( z 2 − 1) K (1 − z −2 )
H ( z) = = e 2 jθ − 2 r cos θ e jθ + r 2
- z 2 + 0 .878 1 + 0 .878 z − 2 avec K =
2 jθ
e −1
- b0=K b1=0 b2=-K a0=1 a1=0 a2=0.878
- L’équation de récurrence est : y(n)= K.x(n)-K.x(n-2)-0.878y(n-2)
Exercice 2 : On souhaite approcher un filtre idéal passe-haut de second ordre par un filtre à réponse
impulsionnelle infinie, synthétisé par la méthode des pôles et des zéros. Ce filtre doit répondre aux spécifications
suivantes
- Fréquence de coupurefc=2 kHz - Largeur de transition à 3db: ∆f=240 Hz
- Fréquence d'échantillonnage : fe=8 kHz
- Tracer H(f) idéal
- Donner l’expression théorique de H(z)
- Placer les pôles et les zéros pour obtenir une approximation de H(f)
- Déterminer l’expression mathématique exacte de H(z) en déterminant K
- Déterminer l’équation aux différences.
K ( z −1)(z −1) K ( z 2 − 2z + 1) K ( z 2 − 2z + 1)
H ( z) = = =
θ = ( f c / f e )× 360 ° = 2 / 8 * 360 = 90° ( z − re− jθ )(z − re jθ ) z 2 − 2zr cosθ + r 2 z 2 + 0.82
∆f 3 db 1 + r 2 + 2r cosθ
r = 1− π = 0.906 avec K= = 0.455
fe
4
Exercice 3 : Considérer la fonction de transfert passe bas dénormalisée H(p) suivant : 2.25/(p2+0.3p+2.25).
Utiliser la méthode de l'invariance impulsionnelle puis celle des pôles et zéros pour transformer ce filtre en
numérique.
Solution
2.25 1.4925 1.292z
H ( p) = = 1.5075 in → H ( z) = 2
variance impulsionnelle
p + 0.3 p + 2.25
2
( p + 0.15) + 1.4925
2 2
z − 0.1347z + 0.7408
ω0 e − aTe z −1 sin(ω 0Te )
Te
( p + a ) + ω0
2 2
1 − 2e −aTe z −1 cos(ω 0Te ) + e −2 aTe z −2
1.292 z ( z + 1)
Remarque On peut renforcer le caractère passe-bas en rajoutant un zéro en -1 soit H ( z ) =
z − 0.1347 z + 0.7408
2
Exercice 4 : Considérer la fonction de transfert passe bas H(p) dénormalisée suivante : 6/(p2+5p+6). Utiliser la
méthode de l’invariance impulsionnelle pour transformer ce filtre en numérique.
6 1 6 −6
Solution H A ( p ) = = + h a ( t ) = 6 e − 2 t U ( t ) − 6 e − 3 t U (t )
p+2 p+3 p+2 p+3
6Te − 6Te
h ( n ) = 6Te e −2 nTe U ( n ) − 6Te e −3 nTe U ( n ) H ( z ) = 2Te −1
+
1− e z 1 − e 3Te z −1
Exercice 5 : On considère le filtre passe-bas analogique normalisé à un pôle unique HN(p)=1/(p+1). On considère
que la fréquence de coupure normalisée à -3db se produit pour fc=fe /10 . Par transformation bilinéaire, trouver le
filtre numérique passe-bas équivalent H(z).
2 π 0.65 ωA 0.65 / Te
ωA = tg = H ( p) = H N ( p) p / ω = =
Te 10 Te A
p + ω A p + 0.65 / Te
2 1 − z −1 0.65 / Te 0.245(1 + z −1 )
p= H ( z ) = H ( p ) 2 1− z −1 = =
Te 1 + z −1 p=
Te 1+ z −1 2 1 − z −1 1 − 0.509z −1
+ 0.65 / T
Te 1 + z −1
e
Exercice 6: Soit la fonction de transfert passe bas dénormalisée H(p) suivante : (p+0.1)/((p+0.1)2+ωa2). Sachant
que ωa=4, utiliser la méthode de la transformée bilinéaire pour transformer ce filtre en numérique. Le filtre
numérique aura sa résonnance pour ωN = π/(2.Te).
1 − z −1
4 + 0.1
2 π 1 − z −1 1 + −1 0.125 + 0.006 z −1 − 0.118z −2
ωa = 4 = tg Te = 0.5 p = 4 H ( z) = z =
Te 4 1 + z −1 1 − z −1
2
1 + 0.0006z −1 − 0.950z −2
4 −1
+ 0.1 + 16
1+ z
Exercices supplémentaires
1. Au cours de la transmission d'un signal numérique (échantillonné à une fréquence de 5 kHz) , il a été affecté
d'un bruit sinusoïdal de fréquence f0=250 Hz. On veut éliminer le bruit par l'emploi d'un filtre possédant une
bande de transition ∆f=±50 Hz à -3db.
Concevoir un filtre RII de second ordre par placement des pôles et zéros
- Tracer H(f) idéal
- Déterminer θ et R puis donner le racé des pôles et zéros
- Déterminer H(z) et le gain K.
- Donner l'équation aux récurrences
- Citer un inconvénient des filtres RII de 1er ordre et de 2ème ordre
H a ( p)
Solution3 pôles p1=-0.1 , p2,3=-0.1± 3j H ( z ) = Te
pôles pi de Ha ( p )
Résidus −1 pT
1− z e e p = pi
0 .5 0 .5
H A ( p) = +
p + 0 .1 − 3 j p + 0 .1 + 3 j
0.5Te 0.5Te 1 − e −0.1Te cos(3Te ) z −1
H (z) = + = T
1 − e 0.1Te e3 jTe z −1 1 − e −0.1Te e3 jTe z −1 1 − 2e −0.1Te cos(3Te ) z −1 + e −0.2Te z −2
e
1 −1 −1
H ( p) = + + ha (t ) = e − tU (t ) − te − 2 tU (t ) − e − 2 tU (t )
p + 1 ( p + 2) 2
( p + 2)
4. Soit à construire un filtre numérique passe-bas échantillonné à la fréquence fe=4kHz par la méthode de
l'invariance impulsionnelle. La fonctionmodèle est la réponse en fréquence d'un filtre passe-bas de type
Butterworth du 2ème ordredont la fréquence de coupure à -3dB est égale à 500Hz.
1 ω c2
H N ( p) = H ( p) = H N ( p) p = p / ωc
= 2
p2 + 2 p +1 p + 2ω c p + ω c2
ω c2 ω c2 ω c2 ωc / 2
H ( p) = 2 = = = 2ω c
p + 2ω c p + ω c2
2 2 2
ω ω
2 2
2 ω c2 1 1
p+
ωc +
p+ ω c + c p+ ω c + c
2 2 2 2 2 2
2ωc z −1 sin(ωcTe / 2 )e −ωcTe / 2
H ( z) = Te ω c T e = π / 4 ω c T e / 2 = 0 .555
1 − 2z −1 cos(ωcTe / 2 )e −ωcTe / 2
+ z −2 e −2ωcTe / 2
5.On considère le filtre analogique d’ordre 3 associé à la fonction d’approximation de Butterworth. Réaliser le
filtre numérique passe-bas correspondant en utilisant la transformation bilinéaire. La fréquence de coupure est
fc=1 kHz et la fréquence d’échantillonnage fe=10 kHz.
1 1
H N ( p) = H ( p) = H N ( p) p= p / ωc =
(
(1 + p) 1 + p + p 2
) (1 + p / ωc ) 1 + p / ωc + ( p / ωc )2 ( )
2 π 2 * 0.325
ωc = tg = H ( z ) = H ( p) p= 2 z −1
Te 10 Te Te z +1
0.034( z + 1)3
H ( z) =
(1,325z − 0.675)(1.43z 2 − 1.79z + 0.78) Remarque :
H ( z) = H N ( p) p=
1 z −1
Ω c z +1
= H N ( p) p=
1 z −1
0 .325 z +1
6. En employant un filtre de Butterworth passe-bas de second ordre, concevoirun filtre numérique correspondant
(passe-bas) en utilisant la transformation bilinéaire. La fréquence de coupure est fc=100 kHz et la fréquence
d’échantillonnage fe=1 kHz.
ωa = 2T 0.325rad / s 1 ω c2
H N ( p) = H ( z ) = H N ( p) p = p / ωc
=
e p2 + 2 p +1 p 2 + 2ω c p + ω c2
0.1056( z + 1) 2
H ( z) = H ( p) p= 2 z −1 =
Te z +1 1.65z 2 − 1.8z + 0.65
7.On cherche à réaliser un filtre numérique équivalent au filtre analogique dénormalisé de transmittance :
1
H A ( p) =
- 1 + 0. 2 p
- Déterminer la réponse impulsionnelle ha(t)
- Calculer la réponse en z de ce filtre obtenue par transformation bilinéaire pour une fréquence
d’échantillonnage Te = 0.2.
- Quelle est sa fréquence de coupure ?
- Calculer la réponse en z d’un filtre similaire de même fréquence de coupure que le filtre analogique.
8.En employant un filtre de Butterworth passe-bas de second ordre, concevoir un filtre numérique correspondant
(passe-bas) en utilisant la transformation bilinéaire. La fréquence de coupure est fc=200 Hz et la fréquence
d’échantillonnage fe=1 kHz.
1
H N ( p) =
p + 2 p +1
2
- Déterminer H(z)
- Etablir le tracé des pôles et zéros
- Comparer les fréquences analogiques et numériques pour fc=200 Hz et fc=400 Hz
- Si fc=400 Hz, quel problème cela posera-t-il?
- Quelle méthode alors utiliser? Quel est son inconvénient?
ωa = 2T 0.7265rad / s 1 ω c2
H N ( p) = H ( z ) = H N ( p) p = p / ωc
=
e p2 + 2 p +1 p 2 + 2ω c p + ω c2
9. On désire choisir, dans chacun des cas suivants, la meilleure approche de synthèse de filtres numériques RII
En employant un filtre de Butterworth passe-bas de second ordre, concevoir un filtre numérique correspondant passe-
hauten utilisant la transformation bilinéaire. La fréquence de coupure désirée estfc=250 Hz et la fréquence
d’échantillonnage fevaut 1kHz.
1
H N ( p) =
p + 2 p +1
2
1. Déterminer H(z)
2. Calculer la fréquence de coupure analogique et justifier la différence avec celle numérique.
But du TP : Dans ce TP, on synthétisera un filtre RII par la méthode des pôles et zéros puis par la méthode de
l'invariance impulsionnelle et enfin par transformation bilinéaire en utilisant des filtres analogiques (Chebychev
et Butterworth).
I.On commence par choisir un filtre analogique normalisé Hn(p) (en continu) d'ordre N puis on le dénormalise
pour créer H(p) continu. Ensuite, on utilise une des méthodes de transformation permettant trouver les
paramètres du filtre numérique
clc;clear all; close all;
Dans les systèmes considérés précédemment, n’était considéré qu’une seule fréquence ou cadence
d’échantillonnage fe =1/Te. Mais il arrive que la fréquence maximale fmax soit diminuée ( filtrage passe-bas ) ou
augmentée (le cas d’une modulation). Dans les deux cas, il est plus judicieux d’adapter la fréquence
d’échantillonnage afin de minimiser le temps de calcul. Il s’agira d’une décimation ou sous-échantillonnage
(augmentation de fmax) et d’une interpolation ou sur-échantillonnage (diminution de fmax). A noter, également que
le sur-échantillonnage est souvent utilisé pour augmenter le SNR dans les opérations de CNA. Tous deux
interviennent dans le changement de fréquence et dans les techniques des bancs de filtres
Exemple : Supposons que l’on veuille synthétiser un filtre passe-bas avec les spécifications
suivantes :fp=100Hz fa=300Hz Aa>50 db fe=20kHz.
L’emploi de la fenêtre de Hamming nous donne N=331, soit un nombre de coefficients important. Or le
filtre occupe une petite portion de la bande passante (passe-bas avec fc=1% fe). On peut donc sous-échantillonner
et puis sur-échantillonner sans perte d’informations.
Ainsi en prenant fe=fe/25=800Hz cela nous donnerait N=15 << 331 et on obtient, ainsi, un filtre plus rapide.
C’est le principe des bancs de filtres qui revient à diviser la bande passante en plusieurs filtres à bande étroite
donc des filtres que l’on pourra facilement décimer ce qui représente donc un gain de temps.
Rappelons que lorsqu’on souhaite échantillonner un signal, nous devons respecter le théorème de Shannon à
savoir prendre une fréquence d’échantillonnage fe/ 2 >fmax. Si ce signal a subi des traitements qui ont eu pour
résultat, entres autres, la modification de la fréquence maximale. On parlera alors de sous-échantillonnage ou de
sur-échantillonnage.
x(n) D xD (m )
1 0 2
m
échantillonner directement à fe/D, mais rien ne garantit que fe /D >2 fmax. Il y a donc un risque de repliement
spectral.
Exemple : Le signal suivant est le résultat de la somme de 2 sinusoïdes de fréquences 100 et 300Hz. On le
décime de 2 et de 6. On remarque que lorsque fe=fe/6 il y a repliement puisque la condition de Shanon n’est plus
respectée
2 0.5
0.4
1
0.3
0
0.2
-1
0.1
-2 0
0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1 -1000 -800 -600 -400 -200 0 200 400 600 800 1000
2 0.5
0.4
1
0.3
0
0.2
-1
0.1
-2 0
0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1 -500 -400 -300 -200 -100 0 100 200 300 400 500
2 0.8
1 0.6
0 0.4
-1 0.2
-2 0
0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1 -200 -150 -100 -50 0 50 100 150 200
Pour sous-échantillonner correctement, il va falloir supprimer le contenu fréquentiel de x(n) entre fe /2D et
fe /D puisque le spectre devient périodique de période fe/D et ce avant la décimation.
2 0.5
0.4
1
0.3
0
0.2
-1
0.1
-2 0
0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1 -1000 -800 -600 -400 -200 0 200 400 600 800 1000
2 0.5
0.4
1
0.3
0
0.2
-1
0.1
-2 0
0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1 -500 -400 -300 -200 -100 0 100 200 300 400 500
1 0.5
0.4
0.5
0.3
0
0.2
-0.5
0.1
-1 0
FEI, USTHB
0 0.01 [assiakourgli@gmail.com / http://perso.usthb.dz/~akourgli/
0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1 -200 -150 -100 -50 0 50 100 150 200
103
Traitement Numérique du Signal Master ST M1 2019/2020
Rappels :
Peigne de Dirac : ∑+
)/&+ , -. périodique de période -. décomposable en série de Fourier telle que
∑+ -. = ∑+ En discret ∑ , %9 = ∑:&1
7 ;
0) = )/&+ , )/&+ 4 '/ 4
1 1 56) 1 56'
83 <
23 23 :
&)
Fonction de transfert @ # = ∑+ ∑
CDEF;
Réponse en fréquence
= % = .∑, %9 @ =! ∗ ∑ ,> ?
K3 :&1 'K3
: '/ :
@ = ∑ !> ? Périodisation LM
K3 :&1 'K3 NM
: '/ : O
x(n) D xD (m)
n 2 . D DD …2D 2
Théoriquement, si le signal d’origine a été bien échantillonné, la reconstruction parfaite est possible, donc
le sur-échantillonnage aussi. Cependant, la formule de reconstruction (avec les sinus cardinaux) est compliquée
à mettre en pratique. On fait appel à des idées plus simples [25]:
Le signal d’origine n’a pas de contenu fréquentiel au-delà de fe/2. Le processus d’interpolation est
susceptible d’en ajouter. On le fait donc suivre d’un filtrage passe-bas à la coupure fe/2 pour supprimer le contenu
ajouté. C’est la raison pour laquelle la méthode d’interpolation n’a pas beaucoup d’importance.
Ce qui en fréquence nous donne : !" = !Q# = 4 56K"2. = 4 56K"/K. R ce qui signifie que le signal original
et sur-échantillonné ont le même spectre, on ne fait que parcourir le cercle D fois.
0.4
1
0.3
0
0.2
-1
0.1
-2 0
0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1 -1000 -800 -600 -400 -200 0 200 400 600 800 1000
2 0.25
0.2
1
0.15
0
0.1
-1
0.05
-2 0
0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1 -2000 -1500 -1000 -500 0 500 1000 1500 2000
2 0.08
1 0.06
0 0.04
-1 0.02
-2 0
0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1 -4000 -3000 -2000 -1000 0 1000 2000 3000 4000
Après sur-échantillonnage, la forme du spectre reste la même aucune énergie n’a été ajoutée au signal),
mais la fréquence d’échantillonnage est maintenant de Dfe. Le sur-échantillonnage fait apparaître ce que l’on
appelle des spectres miroirs et pour obtenir une signal dont la forme temporelle correspondrait au signal x(n) que
l’on aurait échantillonné à une fréquence Dfe, il faut supprimer ces spectres miroirs et donc filtrer passe-bas le
signal sur-échantillonné y(m) avec un filtre ayant une bande atténuée à partir de fe/2 .
2 0.5
0.4
1
0.3
0
0.2
-1
0.1
-2 0
0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1 -1000 -800 -600 -400 -200 0 200 400 600 800 1000
2 0.5
0.4
1
0.3
0
0.2
-1
0.1
-2 0
0 0.02 0.04 0.06 0.08 0.1 0.12 -2000 -1500 -1000 -500 0 500 1000 1500 2000
2 0.4
1 0.3
0 0.2
-1 0.1
-2 0
0 0.02 0.04 0.06 0.08 0.1 0.12 -4000 -3000 -2000 -1000 0 1000 2000 3000 4000
2. Décomposition polyphases
La décomposition polyphases intervient dans les bancs de filtres, elle permet une représentation simple
pour le sur et sous échantillonnage [21]. Quand un signal subit une décimation d’un facteur 2, cela revient à
séparer les échantillons avec un indice pair de ceux avec un indice impair, ou encore, à séparer deux phases
distinctes ayant des délais (aussi appelé écarts de phase) différents dans le vecteur [22]. Si on considère une
décimation d’un facteur M, le principe est exactement le même mais avec M phases distinctes.
S # = ∑+
)/&+ ℎ 2 # & ) +# &1 ∑∞= ∞ℎ 2 1 # 2
En posant : V # = ∑+
)/&+ ℎ 2 #
&)
et V1 # = ∑+
)/&+ ℎ 2 1 # &) alors S # = V0 Q#2 R # 1 V1 Q#2 R
P0(z2)
H(z) x(n) y(n)
z-1
x(n) y(n)
P1(z2)
x(n) P0(z2)
y(n)
-1
z
2
P1(z )
FEI, USTHB [assiakourgli@gmail.com / http://perso.usthb.dz/~akourgli/
106
Traitement Numérique du Signal Master ST M1 2019/2020
P0(z4) P0(z4)
x(n) y(n) x(n) y(n)
-1
z z -1
4
P1(z ) P1(z4)
z-1 z-1
P2(z4)
P2(z4)
z-1
P3(z4) z-1
P3(z4)
D’une manière plus générale, et pour un entier M donné, H(z) peut être décomposé en :
B/ #
S # = ∑:&1 W V Q#9 R oùVB # = ∑+ W #
W )/&+ ℎ 9
Cette manière d’écrire la fonction de transfert H(z) est appelée la représentation polyphase du filtre H(z) et
les filtres Ek(z) sont dit ses composantes polyphases. Elle permet de réaliser une fonction de transfert équivalente
avec un système dont la plupart des éléments fonctionnent à une cadence réduite.
Décimation: Prenons le cas du décimateur [20], si on considère un filtre FIR de longueur N multiple de M
(Rajouter des zéros si nécessaire pour cela) ayant pour réponse impulsionnelle h(n), dont la transformée en Z est
Y
&1
donnée par S: # = ∑[&1
)/ ℎ # &) = ∑:&1
B/ ∑Z/ ℎ X9
<
W # & Z:*B !
Y
&1
S: # = ∑:&1
B/ #
&B ∑<
/ ℎ X9 W #: &Z
= ∑:&1
B/ # VB #
&B :
Où VB # = ∑Z/ ℎ X9 W # &Z
[/:&1
En considérant le schéma bloc d’un banc de filtres, on voit que le filtrage est classiquement appliqué avant
la décimation. Cela signifie donc que le signal est filtré par les M filtres et que l’on ne garde ensuite qu’une phase
parmi les M de chaque résultat de filtrage. Cette manière de faire la transformation n’est pas efficace car elle
entraîne beaucoup de calculs qui ne seront pas exploités par la suite. La notion de représentation polyphase
revient à inverser ce schéma en faisant en premier lieu l’opération de sous-échantillonnage, puis une opération
de filtrage.
↑( H(# " )
⇔
H(#) ↑(
x(n) y(n) x(n) y(n)
Soit le signal x(n) à décimer de 3 avec un filtre passe-bas h(n) dont la fonction de transfert est
La première étape est de décomposer le filtre h(n) en trois filtres P0 (z3), P1 (z3) et P2 (z3). Considérons la
première composante du filtre polyphases P0 (z3)= h(0)+h(3)z-3+ h(6) z-6+h(9)z-9
Ainsi
y0(n) = x(0)h(0)
+ x(1)h(0)z-1 + x(2)h(0)z-2 + [x(3)h(0)+ x(0)h(3)]z-3
+ [x(4)h(0)+x(1)h(3)]z-4 + [x(5)h(0)+x(2)h(3)]z-5 + [x(6)h(0)+x(3)h(3)+ x(0)h(6)]z-6
+[x(7)h(0)+x(4)h(3)+ x(1)h(6)]z-7+[x(8)h(0)+x(5)h(3)+x(2)h(6)]z-8+[x(9)h(0)+x(6)h(3)+x(3)h(6)+x(0)h(9)]z-9
+[x(10)h(0)+x(7)h(3)+x(4)h(6)+x(1)h(9)] z-10 + [x(11)h(0)+x(8)h(3)+x(5)h(6)+x(2)h(9)] z-11
+ [x(12)h(0)+x(9)h(3)+x(6)h(6)+x(3)h(9)] z-12 + ……;
Puisque, le filtre est suivi d’une décimation par 3, on ne retiendra que les éléments de puissance multiple
de 3 (en rouge) ce qui au final revient à
3 P0(z)
x(n) y0(n)
FEI, USTHB [assiakourgli@gmail.com / http://perso.usthb.dz/~akourgli/
108
Traitement Numérique du Signal Master ST M1 2019/2020
Soit le x(n) à interpoler de 3 avec un filtre passe-bas h(n) ayant pour fonction de transfert:
La première étape est de décomposer le filtre h(n) en trois filtres P0 (z3), P1 (z3) et P2 (z3). Considérons la
première composante du filtre polyphase P0 (z3)= h(0)+h(3)z-3+ h(6) z-6+h(9)z-9
Après interpolation x0(n) = {x(0),0, 0, x(1), 0, 0, x(2), 0, 0, x(3)} X0(z)=x(0)+x(1)z-3+ x(2) z-6+x(3)z-9 +x(4)z-12 +…
Sachant que P0 (z3)= h(0)+h(3)z-3+ h(6) z-6+h(9)z-9 P0(z) = h(0)+h(3)z-1+ h(6) z-2+h(9)z-3
Après Interpolation par 3 y0(n) = x(0)h(0)+ 0 z-1 + 0 z-2 + [x(1)h(0)+ x(0)h(3)]z-3 + 0 z-4 + 0 z-5 +
[x(2)h(0)+x(1)h(3)+ x(0)h(6)]z-6 + 0 z-7 + 0 z-8 + [x(3)h(0)+x(2)h(3)+x(1)h(6)+x(0)h(9)] z-9 + 0 z-10 + 0 z-11 + +
[x(4)h(0)+x(3)h(3)+x(2)h(6)+x(1)h(9)] z-12 +..
Ainsi, considérons un décimateur de 4 équipé d’un filtre anti-recouvrement H(z) auquel on applique la
décomposition polyphases :
En appliquant la première identité noble, on obtient une forme permutée permettant d’effectuer le filtrage
à des taux plus faibles
Ainsi, chaque ligne de la décomposition polyphases constitue une réponse à un taux M fois plus faible,
déphasée d’un échantillon par rapport à la ligne précédente.
Interpolation : En suivant une démarche similaire on peut construire la matrice polyphases d’un banc de
filtres de synthèse. Dans ce cas, l’opération de sur échantillonnage qui précède le filtrage est remplacée par une
étape de filtrage suivie d’un sur-échantillonnage. On opère de la même façon pour l’interpolateur équipé d’un
filtre d’interpolation H(z).
Exemple : Décomposition polyphases pour une interpolation D=4. On utilisera la structure transposée plus
pratique dans ce cas.
En appliquant la deuxième identité noble, on gagne en termes de temps de traitement puisque tout le
filtrage s’effectue au taux le plus faible.
Ce type de matrice polyphases, correspondant à un banc de filtres de synthèse, est parfois appelé matrice
FEI, USTHB [assiakourgli@gmail.com / http://perso.usthb.dz/~akourgli/
110
Traitement Numérique du Signal Master ST M1 2019/2020
polyphases de type 2 par opposition avec les matrices polyphases de type 1 qui correspondent à un banc de filtres
d’analyse [22] .
Remarque : A chaque instant d’échantillonnage en sortie du système une seule composante polyphase est
non nulle. Cela conduit au modèle de commutateur où on ne considère à chaque instant d’échantillonnage que
la composante non nulle [24]. Un commutateur en sortie peut jouer le rôle des délais et additionneurs. En effet,
l'interpolation par D = 4, dans cet exemple, insère 3 zéros sur 4. A un instant donné, une seule branche contribue
à la sommation et les autres sont à zéro.
3. Systèmes multi-cadences
Le traitement multi-cadences, traitement qui implique plusieurs fréquences d’échantillonnage, peut avoir
plusieurs motivations. Il peut être le fruit de contraintes liées à une application dans laquelle les flux de signaux
numériques à traiter et à générer doivent l’être à des fréquences d’échantillonnage différentes [20]. Cette situation
peut se rencontrer par exemple dans les applications audio où co-existent plusieurs fréquences
d’échantillonnage : 32kHz (diffusion); 44,1 kHZ (disque compact numérique) ; 48kHz (bande sonore numérique) ;
96kHZ ; …
Cette approche peut aussi permettre d’améliorer les caractéristiques d’implantation d’un algorithme en
adoptant pour chacune de ses étapes un échantillonnage que l’on appelle critique dans le sens où la fréquence de
traitement est choisie égale à la fréquence de Nyquist. Ainsi la chaîne de traitement correspondant à l’émetteur
bande de base d’un modulateur DPQSK où l’on a trois fréquences de traitement : la fréquence bit, la fréquence
symbole et la fréquence d’échantillonnage [20]
Supposons que l’on veuille changer numériquement la fréquence d’échantillonnage d’un signal x(n) et
passer d’une fréquence d’échantillonnagef1 à une fréquence d’échantillonnage f2 où f2/f1=L/M. Cela conduit donc
à sur-échantillonner le signal x(n) d’un facteur L puis à le décimer d’un facteur M. Le sur-échantillonnage par L
doit être suivie d’un filtre passe-bas coupant à fe/2L pour supprimer les spectres dus au sur-échantillonnage où
fe est la fréquence d’échantillonnage du signal de sortie. La décimation par M doit être précédée d’un filtrage anti-
repliement ayant une fréquence de coupure de fe/2M, si fe est la fréquence d’échantillonnage du signal d’entrée.
On ne conservera alors qu’un seul des deux filtres passe-bas : celui dont la fréquence de coupure est la plus basse
↑c HL(#) HM(#) ↓9
x(n) y(n)=xL/M(n)
Exemple : On souhaite sur une réduction de fréquence d’échantillonnage d’un signal audio d’un facteur
2/3 (passage de 48kHz à 32kHz). Voici la structure polyphase conduisant évaluer le filtre à la fréquence la plus
basse possible. En appliquant la décomposition polyphase pour le sur-échantillonnage, on obtient la structure ci-
contre). A noter que décimation et interpolation ne sont pas commutatives sauf si L et M sont premiers entre eux
[20].
On commence par garder le filtre passe-bas dont la fréquence de coupure est la plus petite. Puis on
décompose le filtre suivant l'interpolation.
3. Bancs de filtres
Un banc de filtre est un ensemble de filtres numériques travaillant en parallèle et découpant la bande de
fréquence en K sous bandes. Le traitement du signal numérique traditionnel utilise des blocs comme :
l’additionneur, multiplicateurs et des retardateurs. Dans les systèmes banc de filtres, nous avons deux nouveaux
blocs que sont le décimateur et l’interpolateur [23].
- Un banc de filtres est un ensemble de filtres, avec une entrée ou une sortie commune
- Un banc de filtres est un ensemble de filtres (H0(z);….;HM-1(z)) tel que tel que les bandes passantes des filtres
forment (+/-) une partition de l’ensemble des fréquences et leurs réponses en fréquence se déduisent les unes
FEI, USTHB [assiakourgli@gmail.com / http://perso.usthb.dz/~akourgli/
113
Traitement Numérique du Signal Master ST M1 2019/2020
des autres par certaines relations telles que SB # = S A#4 & G = S Q#I B R
CDEF
<
- Chaque sortie du filtre d’un banc est à bande étroite, on peut donc la décimer pour la traiter plus rapidement.
- Un banc de filtres est dit “à reconstruction parfaite” s’il existe un ensemble de filtres (dits “de synthèse”) qui,
lorsqu’on les applique aux sorties du banc initial et qu’on somme leurs sorties, reconstruit parfaitement le
signal d’entrée (à un retard près).
Ces deux cas sont représentés sur la figure ci-dessus. Ce système est appelé banc de filtres d’analyse, et les
filtres Hk(z) sont les filtres d’analyse. Ce banc décompose le signal x(n) en M signaux vi(n) appelés signaux de
sous-bandes. Le système de gauche est appelé banc de filtres de synthèse et les filtres Fk(z) sont les filtres de
synthèse. Il combine les M signaux wk(n) en un seul signal y(n)[21]
Les bancs de filtres d’analyse et de synthèse sont généralement associés, le premier décompose un signal
pour appliquer un traitement à chaque signal de sous-bande et le second re-combine les signaux de sous-bandes
traités pour construire le signal modifié. Un tel système d’analyse/synthèse est appelé banc de filtres à
reconstruction parfaite quand en l’absence de tout traitement dans les sous-bandes, le signal de sortie y(n)=x(n-
k).Ces filtres sont alors dits bi-orthogonaux, ils sont étroitement liés à la théorie des ondelettes objet du prochain
chapitre.
Dans un banc de filtres à 2 canaux on sépare le signal d’entrée en un signal passe-haut et un signal passe-
bas, respectivement d(n) et a(n). On appelle H0(z) et H1(z) les filtres passe bas et passe haut du banc d’analyse et
F0(z) et F1(z) les filtres passe bas et passe haut du banc de synthèse. On appelle X(z) la transformée en Z du signal
x(n). L’inconvénient de cette structure est que le nombre d’échantillon est multiplié par 2 en sortie du banc
d’analyse. Si les filtres du banc d’analyse H0(z) et H1(z) ont des gabarits correspondant à ceux représentés ci-
contre, il est alors possible d’introduire un échantillonnage critique en décimant par deux en sortie des filtres du
banc du d’analyse.
0.5
0
0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5
On remarquera qu’il n’y a pas de filtre anti-repliement lors de la décimation, en fait on se sert du fait que le
spectre soit à bande limitée (ici la moitié de la bande de fréquence). C’est un cas particulier où le signal est à bande
limité égale à fe/M. Ainsi on peut alors récupérer le signal original à partir du signal sous échantillonné. De m^me
pour le sur échantillonnage puisque on se sert de filtres de reconstruction qui eux aussi occupent la moitié de la
bande de fréquence.
!g # = d # h # !g # = h # QS # ! # S # ! # R
1
!g1 # = h1 # QS1 # ! # S1 # ! # R
1
!g # = Qh # S # h1 # S1 # R! # Qh # S # h1 # S1 # R! #
1 1
Rappelons que l’expression de la sortie du banc de filtre est donnée par [20,21]:
1 1
!g # = Qh # S # h1 # S1 # R! # Qh # S # h1 # S1 # R! #
2 2
h # S # h1 # S1 # = 2 # &B
h # S # h1 # S1 # =0
ou h # = S1 # &1 4 h1 # = S # &1
• Une première famille de filtres dite Q.M.F. (Quadrature Mirror Filter) permet d'éliminer le terme en X(-z)
• La deuxième famille de filtres dits C.Q.F. (Conjugate Quadrature Filter) ont les propriétés suivantes:
Série n°5
Solution
0,2,4,6,8,…. 0,3,6,9,12,…. 0,4,8,12,16,….
0,0,1,0,2,0,3,0,4,0,…. ou 0,0,1,1,2,2,3,3,4,4,…. ou 0,0.5,1,1.5,2,2.5,3,3.5,4,…
0,2,4,6,8,…. puis 0,0,2,0,4,0,6,0,8,8,…. En inversant 0,0,1,0,2,0,3,0,4,0,…. puis 0,1,2,3,4,….
0,3,6,9,12,…. puis 0,0,3,0,6,0,9,0,12,0,…. En inversant
0,0,1,0,2,0,3,0,4,0,5,0,6,0,7,0,8,0,9,0,10,0,11,0 12,0…. puis 0,0, 3,0,6,0,9,0,12,..
Non Invariance de la décimation 0,2,4,6,8,…. Ou 1,3,5,7,9,….
2. Soit le signal suivant, on veut le décimer par 2 puis l'interpoler par 2 également
0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17
- Donner le tracé du signal décimé puis interpolé
- On considère que la fréquence du signal original est 10kHz, donner la fréquence à la sortie du décimateur
puis celle de l'interpolateur.
Solution
1 2 3 4 5 6 7 8 9
0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17
Après décimation . = 5WS# Après interpolation . = 10WS#
↓2 ↓2
- Préciser le rôle de H(z) dans chaque cas
-150 -50 50 150 Hz
H(z)
x(n) y(n) x(n) y(n)
↓4 H(z) ↓4
x(n) y(n) x(n) y(n)
↑2 ↑2 H(#)
x(n) y(n) x(n) y(n)
↑4 ↑4
-
H(#)
x(n) y(n) x(n) y(n)
Solution
-600 -450 -300 -150 -75 -50 50 75 150 300 450 600 Hz
C'est pour cette raison qu'on place un filtre passe-bas à . /2 = ./ 29 Période de -75 à 75 Hz
H(z) ↓2
x(n) y(n)
-75 -50 50 75 Hz
↓4
Période de -37,5 à 37,5 Hz
H(z)
x(n) y(n)
↑2 H(#)
x(n) y(n)
-300 -50 50 300H
↑4 H(#)
x(n) y(n)
j = .
= 75, j = .l
= 37,5 j = = 150, j = = 150
4. Donner l'expression des filtres polyphases pour les 4 cas de l'exercice précédent ainsi que la décomposition
polyphases initiale et finale pour chaque cas.
Solution
- M= 2 S # = ∑+
)/&+ ℎ # &) =∑+
)/&+ ℎ 2 # & ) +# &1 ∑+
)/&+ ℎ 2 1 #& )
On pose V # = ∑+
)/&+ ℎ 2 # &) V1 # = ∑+
)/&+ ℎ 2 1 # &) S # = V # # &1 V1 #
- M=4
S # = ∑+
)/&+ ℎ 4 # &l) +∑+
)/&+ ℎ 4 1 # &l)&1 +∑+
)/&+ ℎQ4 2 R# &l)& ∑+
)/&+ ℎ 4 3 # &l)&
On pose V # = ∑+
)/&+ ℎ 4 # &) V1 # = ∑+
)/&+ ℎ 4 1 # &) nV # = ∑+
)/&+ ℎ 4 2 # &)
V # = ∑+
)/&+ ℎ 4 3 # &) S # = V #l # &1 V1 # l + # & V # l #& V #l
5. Soit le filtre S # = 1 0.5 # &1 0.5 # & 0.9 # & 0.8 # &l 1.2 # &n
- S # =V # # &1 V1 #
Solution
V # = ∑+ )/&+ ℎ 2 #
&)
= 1 0.5 # &1 0.8 # &
V1 # = ∑)/&+ ℎ 2
+
1 # &) = 0.5 0.9 # &1 1.2 # &
- S # = V # # &1 V1 # #& V #
V # = ∑)/&+ ℎ 3 # = 1 0.9 # &1
+ &)
V1 # = ∑+
)/&+ ℎ 3 1 # &) = 0.5 0.8 # &1
V # = ∑+ )/&+ ℎ 3 2 # &) = 0.5 1.2 # &1
FEI, USTHB [assiakourgli@gmail.com / http://perso.usthb.dz/~akourgli/
121
Traitement Numérique du Signal Master ST M1 2019/2020
0.45
1.5
0.4
1
0.35
0.5
0.3
0 0.25
0.2
-0.5
0.15
-1
0.1
-1.5
0.05
-2 0
0 0.002 0.004 0.006 0.008 0.01 -2000 -1500 -1000 -500 0 500 1000 1500 2000
Solution
7. Un signal audio est enregistré a été transmis à une fréquence de 30 kHz, à la réception avant de le convertir
en analogique, on opère un changement de fréquence telle que la nouvelle fréquence soit 45Hz.
’ ’
fc= f1/2= fe /6=15 fc= fe /4=22.5
FEI, USTHB [assiakourgli@gmail.com / http://perso.usthb.dz/~akourgli/
122
Traitement Numérique du Signal Master ST M1 2019/2020
8. Considérer une porte de largeur 5, vérifier les propriétés liant les foncions de transfert pour la décimation
et l'interpolation par 2
Interpolation par M @ # = ! # :
Solution
1
@ # = qQ # 1 # 1/ 1 # &1/ # &1 R Q #1 # 1/ 1 # &1/ # &1 Rr = #1 1 # &1
2
Solution
S =s s1 4 & 56K23
s 4 &l56K23 s 4 &u56K23 ⋯
&
S1 =s s1 4 56K23
s 4 l56K23 s 4u56K23 ⋯
S =s s1 4 & 56K23
s 4 &l56K23 s 4 &u56K23 ⋯
S1 =s s1 4 & 56K23
s 4 &l56K23 s 4 &u56K23 ⋯
Solution
11. On souhaite réaliser un banc de filtres pour M=3, donner la structure d'analyse et reconstruction
correspondante.
Exercices supplémentaires
1. On dispose d’un signal audio échantillonné à f1=48 kHz que l’on souhaite convertir au standard CD de
= yC x
K1 w xn
K
f2=44,1 kHz par décimation et interpolation. Solution
2. Soit le signal suivant :h(n) = sin(2*pi*150*n) + sin(2*pi*350*n) et le module de la TFD dont les tracés
respectifs sont donnés ci-dessus :
- Déterminer Te.
- Quel est l’intérêt de l’interpolation ?
- Tracer le signal (les 20 premières valeurs) et la TF obtenue pour une interpolation de 2 puis de 4(sans
filtrage passe-bas)
- Pourquoi emploie-t-on un filtre passe-bas lors l’interpolation, où le place-t-on ?
- Quel est le type du filtre de Tchebeychev employé dans ce cas et pourquoi ?
- Donner la décomposition polyphases finale pour une interpolation de 2 en donnant l’expression des
filtres polyphases
x(n) 3 y(n)
-fmax fmax=fe/5
- Tracer la TFTD obtenue après décimation (fe=30kHz)
- On veut placer un filtre anti-repliement, ou doit-on le placer ? (justifier), tracer la TFTD à nouveau
- Donner la décomposition polyphases initiale et finale en donnant l’expression de chaque filtre.
Solution : Examen 18/19
4. Soit le signal suivant :h(n) = sin(2*pi*150*n) + sin(2*pi*350*n) et le module de la TFD dont les tracés
respectifs sont donnés ci-après :
2 0.7
1.5 0.6
1
0.5
0.5
0.4
0
0.3
-0.5
0.2
-1
-1.5 0.1
-2 0
0 0.005 0.01 0.015 0.02 -1000 -500 0 500 1000
- Tracer le signal et la TF obtenue pour une décimation de 2 puis de 4(avec filtrage passe-bas)
- Donner la décomposition polyphases finale pour une décimation de 4 en donnant l’expression des
filtres polyphases
- Quel est l’intérêt de la décomposition polyphases
- Pourquoi emploie-t-on un filtre passe-bas lors de la décimation
- Solution : Examen 17/18
TP n°5
1. Interpolation et Décimation
1. Etablir le lien entre période du signal original et celui sous-échantillonné, faire de même pour les
fréquences.
2. Expliquer la différence entre decimate(x,6); entre x(1:6:N);
L=2;Y1 = fft(x1,NF);Y1=fftshift(Y1);axef=L*Fe*(-1/2:1/NF:1/2-1/NF);
subplot(3,2,4); plot(axef,abs(Y1/N1)); grid
L=4 ;Y2 = fft(x2,NF); Y2=fftshift(Y2);axef=L*Fe*(-1/2:1/NF:1/2-1/NF);
subplot(3,2,6); plot(axef,abs(Y2/N2)); grid
5. Pour expérimenter avec le changement de fe et ses conséquences : Téléchargez le son Vous avez du
courrier en attente.wav.
- Ecoutez et Observez son spectre.
- Décimez-le (sans filtre passe-bas), écoutez et observez le spectre du résultat.
- Interpolez-le de deux manières différentes (zéro et échelon). Écoutez et observez son spectre.
2. Décomposition polyphases
3. Bancs de filtres
La plupart des signaux du monde réel ne sont pas stationnaires, et c’est justement dans l’évolution de leurs
caractéristiques (statistiques, fréquentielles, temporelles, spatiales) que réside l’essentiel de l’information qu’ils
contiennent. Les signaux vocaux et les images en sont des exemples courants [ondelettes]. Rappelons que
l’analyse de Fourier permet une caractérisation globale du signal (on intègre de -∞ à +∞), on perd toute
localisation temporelle ou spatiale, l’idéal est de faire appel à une transformation qui nous apporte l’information
sur le contenu fréquentiel tout en préservant la localisation (temporelle ou spatiale) afin d’obtenir une
représentation temps/fréquence ou espace/échelle du signal. C’est ainsi que deux théories ont été élaborées, la
transformée de Fourier à fenêtre puis la transformée continue par ondelettes.
Les ondelettes, famille de fonctions déduites d’une même fonction, appelée ondelette mère, par opérations
de translations et dilatations, ont trouvé, de par la puissance de leur théorie, des applications dans de nombreux
domaines aussi variés que les mathématiques (analyse, probabilités, fractales), le traitement du signal
(compression, astronomie, sismique), la physique (mécanique quantique, turbulence).
Où h(n) est une fenêtre de pondération (Rectangulaire, Bartlett, Hanning, Hamming, Gaussienne etc.)
+∞ 2
− j 2πfτ
S x (t , f ) = x(τ )h (τ − t )e dτ
*
−∞
Il faut relever n’existe pas qu’une seule TFCT puisqu’elle dépend de : La durée de la fenêtre (choisie pour
que le signal soit supposé stationnaire sur cette durée), la forme de la fenêtre (compromis largeur-hauteur des
lobes), le taux de recouvrement entre les fenêtres.
Le rôle de la fenêtre h(t)(dont l’énergie doit valoir 1) est de découper un voisinage de longueur L du point t,
dans lequel le contenu fréquentiel est analysé. On conçoit qu’il y a un compromis entre la longueur L de h(t) qui
représente la résolution temporelle, qui induit une résolution fréquentielle en fe/L, et la capacité de la TFCT à
suivre des modulations plus ou moins rapides. Ces 2 résolutions évoluent en inverse l’une de l’autre. Il a été
montré (principe d’incertitude d’Heisenberg) que la fenêtre qui réalise le meilleur compromis temps-fréquence
est la fenêtre gaussienne.
Exemple : Analyse d’un signal chirp allant de f1 puis f2 et inversement (voir TP n°6)
0.5 0.5
0 0
-0.5 -0.5
-1 -1
0 0.5 1 1.5 2 0 0.5 1 1.5 2
FFT de y1 FFT de y2
30 30
20 20
10 10
0 0
-100 -50 0 50 100 -100 -50 0 50 100
Comme illustré ci-haut, on perd toute localisation temporelle puisqu’on obtient la même TF pour les deux
signaux. Pour y remédier on calcule la TFCT en effectuant la TFD pour différents intervalles et recouvrements
TFCT avec Tfen=100 Recouv=20%, NF=1024 TFCT avec Tfen=100 Recouv=20%, NF=1024
1.4 1.4
1.2 1.2
1 1
Time
Time
0.8 0.8
0.6 0.6
0.4 0.4
0 20 40 60 80 100 0 20 40 60 80 100
Frequency (Hz) Frequency (Hz)
TFCT avec Tfen=10 Recouv=20%, NF=1024 TFCT avec Tfen=10 Recouv=20%, NF=1024
1.5 1.5
Time
Time
1 1
0.5 0.5
0 20 40 60 80 100 0 20 40 60 80 100
Frequency (Hz) Frequency (Hz)
Les TFDs sur chaque fenêtre glissante fournissent le spectrogramme qui permet d’adapter la TF à la
FEI, USTHB [assiakourgli@gmail.com / http://perso.usthb.dz/~akourgli/
129
Traitement Numérique du Signal Master ST M1 2019/2020
caractérisation des signaux non stationnaires. On obtient, alors, une représentation temps-fréquence permettant
de localiser la distribution de l’énergie simultanément en temps et fréquence. Rappelons que la longueur de la
fenêtre choisie va conditionner le nombre de points fréquentiels (résolution fréquentielle) et la résolution
temporelle (spectre moyen sur la fenêtre).
Il est clair que pour ce cas, la TFCT opérant avec une taille de fenêtre unique ne permet pas de localiser
précisément chaque fréquence. Celle-ci devrait s’adapter en fonction de l’évolution du signal. L’idéal serait de
pouvoir choisir une fenêtre et une forme d’onde (signal oscillant dans une fenêtre temporelle donnée) que l’on
pourrait dilater (pour les basses fréquences) et contracter (fréquences élevées) à volonté.
2. Ondelettes continues
Elle a été introduite par Jean Morlet en 1981 pour résoudre des problèmes de signaux sismiques en recherche
pétrolière.
Partant d’une fenêtre h (dite fonction mère) ayant pour symbole ψ dépendant de t, on peut générer un
ensemble de fonctions de base similaire par dilatation (indice a) et translation (indice b) d'un seul
prototypez{,|
1 s
:
z{,| = zA G ~>0
√~ ~
Où a > 0 est un paramètre d’échelle de contraction (a < 1) ou de dilatation (a > 1) de la fenêtre et b une
translation de la fenêtre. On notera que la norme de ψa,b est conservée lors du changement de facteur d’échelle
(voir figure ci-dessous):
*+
1 s
€z{,| € = • ‚z A G‚ ƒ = ‖z‖
&+ ~ ~
Exemple :
1 0 ≤ ≤
1
La fonction de Haar z =… 1 1
≤ ≤1
0 ~ XX4‡ˆ
s ≤ ≤s
1 {
√{
= z> ? =…
1 ‰&|
s ≤ ≤s ~
ainsiz{,| 1 {
√{ {
√{
0 ~ XX4‡ˆ
1.5
Haar a=0.5, b=0.5
Haar a=1, b=2
1 Haar a=3, b=4
Amplitude de l ondelette de Haar
0.5
-0.5
-1
-1.5
0 1 2 3 4 5 6 7 8
temps
Exemples d’ondelettes
C
z =Š >1 ? 4
7
1 ‰C & C
C‹
√ 5 Š C
- L’ondelette chapeau mexicain (Paul à l’ordre 2)
7C
z = 4 4&
1 & 56K‰
C‹C
Š√ 5
- L’ondelette de Morlet
L’ondelette de Morlet est une sinusoïde complexe modulée par une gaussienne. L’ondelette de Paul décroît
plus vite que celle de Morlet et autorise des localisations en temps plus précises. L’ondelette obtenue par dérivée
d’une gaussienne permet des localisations temporelles d’une qualité légèrement inférieure à celle de Paul [25]
En décomposant le signal x(t) sur cette famille, on obtient ainsi les coefficients d’ondelettesI-Œ,• ~, s qui
caractérisent le coefficient de la décomposition du signal x(t) dans cette base, soit :
+
1 *+
s ∗
I-Œ,• ~, s = • z{,|
∗
ƒ = • zA G ƒ
√~ &+ ~
&+
Pour un facteur d’échelle assez grand, la représentation des coefficients d’ondelettes en fonction de b, la
position, donne une représentation de “la forme générale de la fonction”. Par contre un facteur d’échelle faible
correspond à une représentation des singularités.
La transformée en ondelettes est un opérateur linéaire, invariant par translation, et par dilatation. Quelle que
soit l’échelle et quel que soit l’endroit, l’analyse du signal se fait avec la même fonction. La transformée en
ondelettes d’un signal n’est pas unique, elle dépend de l’ondelette mère utilisée. En effet, l’ondelette mère ψ (t)
devra avoir une bonne localisation (nulle en dehors d’un certain intervalle), et devra être oscillante le nombre de
moments nuls correspond au nombre d’oscillations).
On peut montrer que si la fonction analysante (l’ondelette) est correctement choisie, la transformation en
ondelettes est inversible. Le signal x(t) peut être reconstruit après double intégration suivant le facteur d’échelle
a et le paramètre de translation b :
*+ *+
1
=• • I-Œ,• ~, s z{,|
∗
ƒ~ ƒs
&+ &+ ~
Exemple d’application : Analyse du signal chirp avec l’ondelette de Haar (Voir TP n°6)
Absolute Values of Ca,b Coefficients for a = 1 2 3 4 5 ...
127
120
Les coefficients de l’ondelette sont 113
106
représentés en fonction du temps où les 99
valeurs les plus élevées sont de couleur claire. 92
85
On note qu’un coefficient à une amplitude 78
d’autant plus grande que l’ondelette
scales a
71
64
ressemble au signal sur la portion analysée. 57
Lorsque la fenêtre est étroite (onde étroite), on 50
43
observe les hautes fréquences et lorsqu’elle 36
est large, ce sont les coefficients des basses 29
22
fréquences qui sont élevées. 15
8
1
50 100 150 200 250 300 350 400
By scale Values of Ca,b Coefficients for a = 1 2 3 4 5 ...
time (or space) b
Toutefois bien que la transformation
continue en ondelettes permette d’obtenir
une bonne localisation temps-fréquence,
toutefois, l’information engendrée par la 2
COEFS
0
transformation est infiniment redondante -2
~ = 2&6 4 s = W2&6
Si la fonction x(t) est discrétisée, en supposant une période d’échantillonnage égale à 1, pour des raisons de
simplicité, l’équation s’écrit alors : 0Œ Q2&6 , W2&6 R = 2C ∑) zQ26 WR
E
La construction de telles ondelettes peut s'aborder comme un problème de choix de base de décomposition
de signal, mais aussi comme un problème de décomposition en sous-bandes (voir chapitre précédent).
L’idée de l’analyse multirésolution développée par Meyer et Mallat permet d'analyser un signal à différentes
échelles à travers des opérateurs linéaires à des niveaux de résolutions correspondant à différentes bandes de
fréquences spatiales.
Le signal est décomposé de manière itérative par un
banc de filtres. A chaque niveau, le spectre du signal est
partage en deux bandes par un filtre passe-bas et un par
filtre passe-haut. Le filtre passe-bas (fonction d’échelle)
donne les informations grossières tandis que le filtre passe-
haut (ondelettes) encode les détails. Le processus est répété
jusqu’à obtenir la résolution souhaitée. L’avantage de cette
méthode est de n’utiliser que deux filtres. Quand le
processus est arrêté, le signal peut être défini par le jeu de
coefficients correspondant [25]: {PB3, PH3,PH2,PH1}
B
ƒ6*1 • • = P ℎ1 •2 W•~6 •W•
B
L’ondelette doit être très compacte ainsi, ses coefficients doivent être, pour la plupart, proche de zéro. Cette
condition dépend de quelques paramètres comme la régularité de la fonction, le nombre de moments nuls et la
taille de son support. La compacité du support impose que les filtres (h0) qui vont engendrer la fonction d’échelle
(donnant une approximation grossière du signal) et l’ondelette (h1 fournissant les détails) soient à réponse
impulsionnelle finie.
Pour connaître les filtres ℎ , ℎ1 , , 4 1 , dans le cas de bases orthogonales, il suffit de connaître
seulement le filtre ℎ‘[k] puisque les trois autres filtres se calculent à partir de ce dernier (filtres de reconstruction
identiques aux filtres d’analyse mais retournés dans le temps):
Remarque : Dans ce qui suit nous adopterons la convention suivante le premier élément correspond à
l'indice 1 (et non pas 0) pour être en conformité avec les TPs fait sous matlab
Ainsi, considérons par exemple l’ondelette de Daubechies à 6 coefficients telle que la fonction d’échelle h soit :
ℎ = [0.2352 0.5706 0.3252 -0.0955 -0.0604 0.0249],
ℎ1 est obtenu en retournant ℎ dans le temps et en inversant le signe des coefficients impairs, ce qui nous donne :
ℎ1 = [-0.0249 -0.0604 0.0955 0.3252 -0.5706 0.2352]
Filtres de reconstruction identiques aux filtres d’analyse mais retournés dans le temps):
f0= [0.0249 -0.0604 -0.0955 0.3252 0.5706 0.2352], f1= [0.2352 -0.5706 0.3252 0.0955 -0.0604 -0.0249 ]
Ci – dessous sont donnés les coeffcients d’ondelettes de Daubechies 3 avec le tracé des pôles et zéros ainsi que
des réponses en fréquence (Voir TP n°6)
Filtre passe-bas de décomposition Filtre passe-haut de décomposition Filtre passe-bas de décomposition Filtre passe-haut de décomposition
1 1
1
1
Imaginary Part
Imaginary Part
0.5 0.5 0.5
3 5 5 3
0 0
0 0
-0.5
-1
-0.5 -0.5 -1
0 2 4 6 0 2 4 6 -1 1 02 -1
-0.5 0 0.5 1
Real Part Real Part
Filtre passe-bas de Reconstruction Filtre passe-haut de Reconstruction
1 1 Filtre passe-bas de Reconstruction Filtre passe-haut de Reconstruction
1
1
Imaginary Part
Imaginary Part
0.5 0.5 0.5
3 5 5 3
0 0
0 0
-0.5
-1
-0.5 -0.5
0 2 4 6 0 2 4 6 -1
-1 -0.5 0 0.5 1 -2 -1 0 1
Real Part Real Part
Filtre passe-bas et passe-haut de décomposition Filtre passe-bas et passe-haut de Reconstruction
1.5 1.5
1 1
0.5 0.5
0 0
0 0.1 0.2 0.3 0.4 0.5 0 0.1 0.2 0.3 0.4 0.5
Exemple d’application 1 : On suppose que le signal à analser est une rampe et que l’ondelette est celle de Haar.
Décomposition
~6*1 • • = P ℎ •2 W•~6 •W• = P •W 2 •~6 •W• ƒ6*1 • • = P ℎ1 •2 W•~6 •W• = P 1 •W 2 •~6 •W•
B B B B
Approximation Détail
Niveau 0 0,1,2,3,4,5,6,7
,, , ,
1 n ˜ 1 1 1 1 1
√ √ √ √ √ √ √ √
Niveau 1
, , ,
Reconstruction:
,, , ,
1 n ˜ 1 1 1 1 1
√ √ √ √ √ √ √ √
Niveau 1
, , ,
, 0, 0, ,0, , 0, , 0,
1 n ˜ 1 1 1 1 1
√ √ √ √ √ √ √ √
0, , 0, , 0, ,0 ,0
1 1 5 5 9 9 13 13 1 1 1 1 1 1 1 1
, , , , , , ,
2 2 2 2 2 2 2 2
, , , , , , ,
Niveau 0 0,1,2,3,4,5,6,7
Exemple d’application 2 :
On suppose que le filtrage passe-bas ℎ est un opérateur de moyennage entre 2 coefficients et que le passe-haut
ℎ1 est un opérateur qui encode la demi-différence entre 2 coefficients [cours 2013] tels que ℎ ={0.5, 0.5} , ℎ1 ={-
0.5, 0.5} , =2*{0.5, 0.5} , 1 =2*{0.5, -0.5}
Signal original : 11, 9, 5, 7
Décomposition Reconstruction
Approximation Détail Approximation Détail
L’opérateur ℎ donne une représentation approximative alors que ℎ1 donne une représentation plus précise
Il est possible de reconstruire le signal à partir des éléments de la décomposition en sommant le résultat de
l’application des filtres et 1.
Remarque : Le lifting scheme permet une implantation très simple des décompositions en ondelettes et de
leurs opérations inverses en employant la factorisation polyphases [21]
Exemple 3 Voici le résultat de la décomposition du signal ecg, ne sont gardés à la fin que la dernière
approximation et tous les niveaux de détails grâce auxquels on pourra reconstruire le signal original (TP n°6)
0 0
-1000 -1000
0 500 1000 1500 2000 0 500 1000 1500 2000
Coefficients Détail N 1 Coefficients Détail N 2
20 100
0 0
-20 -100
0 500 1000 1500 0 200 400 600
Coefficients Détail N 3 Coefficients d Approximation N 3
500 5000
0 0
-500 -5000
0 100 200 300 0 100 200 300
Les ondelettes sont bien adaptées au débruitage des signaux. Le principe est simple, on décompose le signal
bruité et l’on force à zéro les coefficients qui ne passent pas un seuil. Puis l’on reconstruit le signal (Voir TP n°6).
De la même façon on peut compresser le signal.
Pour une image, on appliquera deux filtres le long des lignes l’un passe-haut (H), l’autre passe-bas (L),
ensuite, on ne considère qu’une colonne sur deux et on applique à nouveau les 2 filtres. En ne considérant qu’une
ligne sur deux, on obtient 4 images de taille réduite de moitié. En réitérant ce processus sur l-image
d’approximation (LL), on obtient une pyramide d’images.
Cet exemple donne le résultat d’une compression jpeg et jpeg2000 (Ondelettes) avec un bpp de 0.2
4 56KH ‰ œ•‡ˆ ž
1. Déterminer la TFCT du signal suivant (h est une fenêtre porte de largeur T):
= › 56K ‰
4 C œ•‡ˆ >
=
⎧ - Q -R4 & œ•‡ˆ ž
2
Solution :
56 K&KH ‰
1
⎪ - Q -R4 & 56 K&KC ‰
œ•‡ˆ
2
>
⎪
-
- -
> ? A> ? -G 4
& 56 K&KH > 0 ?
⎨ 0
2 4
2 0 2 1
⎪ - - -
0.6
1
1 0 ≤ ≤
1
2. Soient les ondelettes suivantes : 0.5
z =… 1
0.5
≤ ≤1
1
TF de ψ(t)
0.4
- La fonction de Haar
0 ~ XX4‡ˆ
ψ(t)
0
0.3
-0.5
0.2
7C
z =£ 1 4& C
-1
0.1
- L’ondelette chapeau mexicain
-1.5 0
0 0.5 1 -1000 -500 0 500 1000
Calculer leur TF et tracer la t f
Ondelette Chapeau Mexicain TF de l Ondelette Chapeau Mexicain
1 0.1
0.09
0.8
Solution : 0.08
Haar zg
Dv
¥¦§C > ?
=¤ 4 &56K
0.6 0.07
C
Dv
TF de ψ(t)
0.06
0.4
C
ψ(t)
0.05
Chapeau mexicain zg
vC
= ¨4 4
0.2
&
0.04
C 0 0.03
0.02
-0.2
0.01
-0.4 0
-10 -5 0 5 10 -2 -1 0 1 2
t f
3. Soit x(t) = A U(t-t0), montrer que les coefficients d’ondelettes sont nuls en dehors d’un domaine
triangulaire du plan. On suppose que l’on utilise une ondelette de Haar de largeur T centrée en 0 ? Calculer
alors les coefficients.
Solution
5. Soit le filtre passe bas de décomposition suivant : ℎ ={-0.1294, 0.2241, 0.8365, 0.4830}, déterminer le filtre
passe-haut de décomposition et les filtres de reconstruction
7. Soit le signal suivant, donner et tracer sa décomposition en ondelettes si la fonction d'échelle est ℎ =[0.5,
0.5] alors la fonction d'ondelettes est ℎ1 =[-0.5, 0.5]
Décomposition
√– √– √– √– √– √– √– √–
Décomposition
Reconstruction
Energie
Exercices Supplémentaires
1. Soit le signal suivant : x(n)={0, 1, -2, 3, -4, 5, -6, 7 }
√ √
- Donner le schéma de reconstruction permettant de passer d’un niveau à un autre.
- Reconstruire le signal à partir du niveau 1.
- Donner le signal décomposé aux niveaux 2 puis 3.
- Vérifier que l’énergie se conserve à tous les niveaux (1, 2 et 3)
- Proposer un moyen de compresser ce signal.
Solution Voir examen 17/18
- L’énergie se conserve-t-elle ?
- Reconstruire le signal à partir du niveau 2.
Solution Voir Ratt 17/18
TP n°6
1. De la TFCT aux ondelettes
clc; clear all; close all;
Te=0.005; Fe=1/Te; Tfen=50; Trec=2; NF=1024 ; t=0:Te:2; f0=1; f1=50 ; y1=chirp(t,f0,2,f1);
figure;subplot(2,2,1); plot(t,y1);title('Chirp de f0=1Hz à f1=50Hz');
y2=chirp(t,f1,2,f0);subplot(2,2,2); plot(t,y2);title('Chirp de f1=50Hz à f0=1Hz');
yy1=fft(y1,NF);yy1=fftshift(yy1);
axef=Fe*(-1/2:1/NF:1/2-1/NF);subplot(2,2,3);plot(axef,abs(yy1));title ('FFT de y1');
yy2=fft(y2,NF);yy2=fftshift(yy2);
axef=Fe*(-1/2:1/NF:1/2-1/NF);subplot(2,2,4);plot(axef,abs(yy2));title ('FFT de y2');
figure;subplot(2,2,1); spectrogram(y1,100,20,NF,Fe);
title('TFCT avec Tfen=100 Recouv=20%, NF=1024');
subplot(2,2,2); spectrogram(y2,100,20,NF,Fe);
title('TFCT avec Tfen=100 Recouv=20%, NF=1024');
subplot(2,2,3); spectrogram(y1,10,2,NF,Fe);
title('TFCT avec Tfen=10 Recouv=20%, NF=1024');
subplot(2,2,4); spectrogram(y2,10,2,NF,Fe);
title('TFCT avec Tfen=10 Recouv=20%, NF=1024');
figure;subplot(2,1,1);COEFS = cwt(y1,1:32,'haar','plot');
subplot(2,1,2);COEFS = cwt(y1,1:129,'haar','plot');
figure; COEFS = cwt(y1,1:32,'haar','3Dlvl');
1. La TF permet-elle de localiser temporellement les fréquences
2. L’emploi de la TFCT permet-il de régler le problème ?
3. Quelle est la fenêtre employée par la TFCT (aider vous du help sur spectogram) et pourquoi ?
4. Faire varier Tfen, quel est l’inconvénient alors ?
5. L’emploi de l’ondelette continue (CWT) permet-il de pallier ce problème ?
6. Quel est l’inconvénient principal ? comment y remédie-t-on ?
Références