Download Analyse du signal: première approche

Analyse du signal
Thierry Audibert
[email protected]
20 décembre 2013
Disponible sur http://www.univenligne.fr sous le nom AnalyseSignal.pdf
Table des matières
1
Objectifs et mode d’emploi
2
2
Transformée de Fourier
2.1 Transformée de Fourier : point de vue élémentaire . . . . . . .
2.2 Produit de convolution (fonctions intégrables) . . . . . . . . .
2.3 Transformée de Fourier et convolution : formulaire . . . . . .
2.4 Les corrigés et démonstrations des sous-sections (2.1) et (2.2)
.
.
.
.
3
3
4
7
8
3
Étude du signal
3.1 Première approche du spectre . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.2 Échantillonnage, théorème de Shannon et Nyquist . . . . . . . . . . . . . . . . .
3.3 Corrigés des exercices des sections (3.1) et (3.2) . . . . . . . . . . . . . . . . .
12
12
16
22
4
Transformations de Fourier discrètes
4.1 Transformation de Fourier discrète . . . . . . .
4.2 Les transformations de Fourier rapides . . . . .
4.3 Transformées en cosinus discrètes . . . . . . .
4.4 Transformée en cosinus discrète en dimension 2
4.5 Corrigés des exercices de la section (4) . . . .
.
.
.
.
.
30
30
31
35
37
38
L’analyse du son
5.1 Enregistrement aux formats .wav et .text avec goldwave . . . . . . . . . . . . . .
42
42
5
1
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
1
Objectifs et mode d’emploi
Le cours sur les séries de Fourier nous apprend que l’on peut représenter une fonction ou un signal
T −périodique (et continu par morceaux) comme somme d’une série trigonométrique
∞
X
2iπnt
cn e T
−∞
dans laquelle les fonctions trigonométriques ont des périodes qui divisent T.
La plupart des signaux n’admettent pourtant pas une représentation aussi simple. Ainsi, les sons
les plus courants, qu’ils soient produits par un instrument, par la voix, par le fonctionnement d’un
moteur, ne sont pas des signaux périodiques. Nous verrons pourtant que l’on peut, en faisant varier
la durée d’observation, faire apparaı̂tre la plupart des sons comme une superposition de signaux
sinusoı̈daux. Les mêmes outils mathématiques et numériques nous permettront d’étudier des signaux qui nécessitent une observation sur de grandes durées comme les variations des hauteurs de
marée. On conçoit avec ces deux exemples que la définition du signal va devenir assez générale :
pour nous ce sera une fonction du temps, suffisamment régulière pour se prêter à nos calculs.
On propose ici une approche quasi-expérimentale qui va nous permettre de débroussailler le
problème. C’est pourquoi il n’est pas indispensable de connaı̂tre tous les résultats de la section
2 sur la transformation de Fourier avant d’aborder la section 3 et le cœur de l’analyse du signal.
Chaque section est suivie d’un résumé (sous forme de formulaire ou de bilan provisoire) ce qui
devrait faciliter une lecture avec aller-retours.
• Plan de lecture possible
1. Commencer par la définition de la transformée de Fourier page 3 et l’exercice 1.
2. On pourra découvrir ensuite la notion de spectre d’un signal avec l’exercice 9 page 13,
3. et poursuivre par l’étude de la problématique de d’échantillonnage avec le théorème de
Nyquist et Shannon (forme simplifiée) en section 3.2.
4. La transformation de Fourier discrète et l’algorithme de transformation rapide sont présentées
en section 4.
5. Des applications concrètes sont proposées dans la section ??. C’est à mon avis là que les
choses s’éclairent.
2
2
Transformée de Fourier
2.1
Transformée de Fourier : point de vue élémentaire
Définition 1 On définit la transformée de Fourier d’une fonction intégrable sur R en posant
Z
ˆ
e−iωx f (x) dx
(2.1)
F(f )(ω) = f (ω) =
R
La définition varie dans la littérature d’un facteur multiplicatif autour de la pulsation et on peut
aussi poser
Z
ˆ
e−2πiωx f (x) dx.
F(f )(ω) = f (ω) =
R
Exercice 1 existence, premières propriétés
Soit f une fonction intégrable sur R.
1. Montrer que sa transformée de Fourier est définie pour tout réel ω.
2. Montrer que fˆ est continue et bornée sur R, donner une majoration de ||fˆ||∞ .
R
3. Calculer la transformée de Fourier de la fonction caractéristique de l’intervalle [a, b] :
(
χ(t) = 1 si a ≤ t ≤ b
χ(t) = 0 sinon
Est-elle intégrable ?
corrigé en 2.4.1
Exercice 2 lemme de Riemann Lebesgue
Soit f une fonction continue par morceaux sur R. On veut montrer que
lim fˆ(ω) = 0;
ω→+∞
1. On suppose que f est une fonction en escalier sur [a, b]. Montrer que
Z b
lim
eiαt f (t) dt
α→+∞ a
a pour limite 0 (α désigne un réel).
2. Montrer que le même résultat est vrai pour toute fonction f, continue par morceaux sur
l’intervalle [a, b].
3. On considère maintenant f intégrable sur R. Montrer que
(a) pour tout ε > 0, il existe un intervalle compact [a, b] tel que
Z
Z
|f | ≤ ε,
|f | ≤ ε,
]−∞,a]
[b,+∞[
(b) En déduire que
Z
∞
lim
α→∞ −∞
3
eiαt f (t) dt = 0.
corrigé en 2.4.2
Exercice 3 dérivation
1. Donner une CS pour que la transformée de Fourier de f soit dérivable. Donner alors une
expression de sa dérivée.
2. Donner une CS pour que la transformée de Fourier de f 0 soit définie. La calculer dans ce
cas.
corrigé en 2.4.3
Exercice 4 transformée de Fourier d’une Gaussienne
2
Soit f définie par f (t) = e−at , avec a > 0.
1. Montrer que sa transformée de Fourier est de classe C ∞.
2. Montrer que fˆ vérifie une équation différentielle linéaire et la calculer.
corrigé en 2.4.4
2.2
Produit de convolution (fonctions intégrables)
1
Définition 2 produit de convolution
Soient f et g deux fonctions continue par morceaux sur R. Lorsque la fonction t → f (x − t)g(t)
est intégrable, on pose
Z
f ∗ g(x) =
f (x − t)g(t) dt
(2.2)
R
On appelle produit de convolution de f et g cette fonction f ∗ g
Exercice 5 existence et propriétés du produit de convolution
Soient f et g continues par morceaux sur R. On suppose f intégrable et g bornée.
1. Montrer que f ∗ g et g ∗ f sont définies sur R et que f ∗ g = g ∗ f
2. Montrer que si f ou g est continue, le produit de convolution f ∗ g est aussi continu.
3. En admettant au besoin la validité d’une formule du type
Z Z
Z Z
f (x, t)dt dx =
f (x, t)dx dt
R
R
R
exprimer une relation entre f[
∗ g, fˆ et ĝ.
corrigé en 2.4.5
1
4
R
Exercice 6 convolution par une gaussienne
2
On note g la fonction définie sur R par g(x) = e−x et, pour s > 0, on pose
1 x
θs (x) = √ g
.
s
s π
1. Tracer les fonctions θ1/n , pour n = 1, 2, ..5, et étudier la limite de la suite de fonctions
(θ1/n )n .
2. Calculer l’intégrale de θs sur R. On rappelle que
Z
√
2
e−t dt = π.
R
3. Pour f continue par morceaux et intégrable sur R, on pose
Z
Z
f ? θs (x) =
f (t)θs (x − t) dt =
f (x − t)θs (t) dt.
R
R
(a) Justifier l’existence et l’égalité de ces expressions.
(b) On suppose que f est, de plus, continue en x0 et bornée sur R. Montrer que
lim f ? θs (x0 ) − f (x0 ) = 0.
s→0
On pourra observer que
Z
f ? θs (x0 ) − f (x0 ) =
Z
f (x0 − t)θs (t) dt −
R
f (x0 )θs (t) dt.
R
(c) Étudier la convergence uniforme.
corrigé en 2.4.6
Théorème 1 formule de Fubini pour les fonctions intégrables 1
On suppose f continue et intégrable sur I × J. Si
– pour toutRx, la fonction y → f (x, y) est intégrable sur J,
– g : x → J f (x, y) dy est continue par morceaux et intégrable sur I,
alors,
ZZ
Z
f = g.
I×J
I
Remarque : avec les hypothèses symétriques :
– pour toutRy, la fonction x → f (x, y) est intégrable sur I,
– h : y → I f (x, y) dx est continue par morceaux et intégrable sur J,
on obtient
ZZ
Z
Z
f= g=
h.
I×J
I
J
1. On comparera au théorème de Fubini pour les séries numériques...
5
Exercice 7 transformation de Fourier, formule d’échange
A toute fonction numérique f, continue par morceaux et intégrable sur R, on associe sa transformée de Fourier définie par
Z
ˆ
f (t)e−iωt dt.
f (ω) =
R
1. Justifier que fˆ est continue et bornée.
2. Montrer que si f et g sont intégrables sur R, il en va de même pour les fonctions f × ĝ et
g × fˆ.
3. On suppose f et g intégrables et continues sur R. Montrer que la fonction de deux variables
φ(x, y) = f (x)g(y)e−ixy
est continue et intégrable sur R2 .
4. Sous ces hypothèses, montrer que
Z
Z
f × ĝ =
g × fˆ.
R
R
Exercice 8 transformée de Fourier et convolution
On considère ici deux fonctions numériques f et g, continues et intégrables sur R.
1. On suppose que, de plus, l’une des deux fonctions f ou g est bornée sur R. Montrer que,
pour tout x ∈ R, la fonction t → f (x − t)g(t) est intégrable sur R et que la fonction
Z
h:x→
f (x − t)g(t) dt
R
est continue sur R. On commencera par le cas le plus simple.
2. Montrer que la fonction de deux variables
(x, t) → f (x − t)g(t)
est intégrable sur R2 .
3. On définit l’intégrale de Fourier d’une fonction intégrable θ, en posant
Z
θ̂ω) =
e−iωs θ(s) ds.
R
Montrer que
fˆ(ω) × ĝ(ω) =
ZZ
φ avec φ(x, t) = f (x − t)e−iω(x−t) g(t)e−iωt .
R2
Peut on établir la formule ĥ(ω) = fˆ(ω) × ĝ(ω), lorsque h est la fonction définie en 1, en
apportant éventuellement des hypothèses supplémentaires ?
6
2.3
Transformée de Fourier et convolution : formulaire
• Si f est une fonction intégrable sur R, sa transformée de Fourier est définie sur R, continue et
bornée sur R sa limite est nulle en ±∞ : limω→+∞ fˆ(ω) = 0 et on a les formules suivantes :
• lorsque la transformée est définie par
F(f )(ω) = fˆ(ω) =
R
e−iωx f (x) dx,
R
Si f et x → xf (x) sont intégrables,
dfˆ(ω)
dω
=
Si f et f 0 sont intégrables,
fb0 (ω)
=
Si ga (t) = e
−at2
−i
Si f, g, f ∗ g sont (définies) et intégrables,
• lorsque la transformée est définie par
gba (ω)
=
f[
∗g
=
F(f )(ω) = fˆ(ω) =
e−iωx x f (x) dx
R
i ω fˆ(ω)
r
,
R
−ω 2
r
π
π
4a
=
e
g
(ω)
a
a 1/4a
fˆ × ĝ.
R
e−2πiωx f (x) dx,
R
Si f et x → xf (x) sont intégrables,
dfˆ(ω)
dω
=
Si f et f 0 sont intégrables,
fb0 (ω)
=
−2 i π
R
e−2 i π ωx x f (x) dx
R
2, i π ω fˆ(ω)
2
2
r
Si ga (t) = e−at ,
gba (ω)
=
Si f, g, f ∗ g sont (définies) et intégrables,
f[
∗g
=
1
7
−π 2 ω
r
π
π
e a
=
g 2 (ω)
a
a π /a
fˆ × ĝ.
2.4
Les corrigés et démonstrations des sous-sections (2.1) et (2.2)
Corrigé n˚ 2.4.1 de l’exercice 1 (intégrale de Fourier)
1. Si f est intégrable sur R, t → f (t)e−iωt l’est aussi : elles ont le même module.
2. Pour tout ω ∈ R,
Z
Z
−iωt
ˆ
dt ≤
|f (t)| dt
|f (ω)| = f (t)e
R
R
on a donc ||fˆ||∞ ≤ ||f ||1 .
3.
Z
e
χ̂(ω) =
−iωt
Z
b
χ(t) dt =
e−iωt dt
a
R
(b − a)
si ω = 0
i −iωt b
 e
si ω 6= 0
a
ω
(b − a)
si ω = 0
= 2
b−a
 sin b−a ω e−i 2 ω si ω 6= 0
2
ω
=
sin b−a w 2
Cette fonction a pour module 2 dont on sait qu’elle n’est pas intégrable sur R
w
(mais elle admet une intégrale impropre).
Corrigé n˚ 2.4.2 corrigé 2 (lemme de Riemann Lebesgue)
1. Les fonctions en escalier : si φ est une fonction en escalier attachée à la subdivision (tj )j du
segment [a, b], on a
Z b
n−1 Z t
X k+1
int
int
=
φ(t)e
dt
α
e
dt
k
a
≤
n−1
X Z tk+1
k=0
tk
k=0
tk
int
e k+1 − eintk dt ≤ 2 sup |αk ||b − a|.
|αk | in
n
La limite est bien 0 lorsque n → +∞.
2. Considérons maintenant f continue par morceaux sur [a, b] et ε > 0. Il existe une fonction
en escalier qui vérifie ||f − φ|| ∞ ≤ ε/|b − a|. Nous avons par ailleurs
[a,b]
Z b
Z b
Z b
int
int
int
φ(t)e dt .
f (t)e dt ≤ (f (t) − φ(t))e dt + a
a
a
– Le premier terme du membre de droite est majoré par ε;
– Pour le second, on sait qu’il existe un rang Nε à partir duquel
Z b
int φ(t)e dt ≤ ε.
a
8
En conséquence, pour tout ε > 0 il existe Nε tel que
Z b
int f (t)e dt ≤ 2ε.
n≥ε⇒
a
3. Considérons maintenant une fonction f intégrable sur R ou sur un intervalle quelconque
I ⊂ R, et ε > 0. Il existe un segment [a, b] ⊂ I tel que
Z
Z
|f | ≤ ε.
|f | −
I
[a,b]
On a :
Z
Z
f (t)eint dt = Z
f (t)eint dt + f (t)eint dt
[a,b]
I\[a,b]
Z
≤
f (t)eint dt + ε.
[a,b]
I
Comme dans la démonstration précédente, il existe un rang Nε à partir duquel
Z b
int f
(t)e
dt
≤ ε,
a
et on conclut de la même façon.
Corrigé n˚ 2.4.3 de l’exercice 3 (dérivation et transformation de Fourier)
1. On
une fonction f intégrable sur R de transformée de Fourier définie par fˆ(ω) =
R considère
−iωt
dt. La fonction (ω, t) → f (t)e−iωt satisfait clairement au théorème de contiR f (t)e
nuité d’une intégrale à paramètre
– g := t → f (t)e−iωt est intégrable sur R;
– ω → f (t)e−iωt est continue ;
– il existe une majoration uniforme |f (t)e−iωt | ≤ |f (t)| par une fonction intégrable sur R.
Par contre, pour la dérivation sous le signe somme, rien n’assure que
∂
g(ω) = −iωf (t)e−iωt
∂ω
vérifie les mêmes hypothèses. On supposera donc pour aller plus loin dans les calculs que
de plus, t → tf (t) est intégrable sur R. Alors
∂
– g := t →
g(ω) = −i t f (t)e−iωt est intégrable sur R;
∂ω
∂
– ω→
g(ω)f (t)e−iωt est continue ;
∂ω
∂
– il existe une majoration uniforme g(ω) ≤ |t f (t)| par une fonction intégrable sur R.
∂ω
Donc, si f (t) et tf (t) sont intégrables sur R, fˆ est de classe C 1 et
dfˆ(ω)
= −i
dω
Z
e−iωx x f (x) dx
R
9
2. On suppose maintenant que f est intégrable et dérivable sur R et que f 0 est intégrable elle
aussi. On peut donc définir la transformée de Fourier de f 0 qui est
Z
0
0
ˆ
e−iωx f 0 (x) dx.
F(f )(ω) = f (ω) =
R
On est alors tenté d’intégrer par parties. Cela donnerait
Z
−iωx
+∞
0
F(f )(ω) = e
f (x) −∞ + iω e−iωx f (x) dx.
R
Dans cette formule les deux intégrales sont définies. En conséquence la partie toute intégrée
a elle même un sens lorsque f et f 0 sont toutes deux intégrables. Pour se convaincre que
cette limite est nulle on intègre les mêmes fonctions sur l’intervalle [0, +∞[ par exemple.
Cela montre que
Z +∞
Z ∞
+∞
e−iωx f 0 (x) dx = e−iωx f (x) 0 + iω
e−iωx f (x) dx
0
0
et là on voit clairement que si la partie toute intégrée a une limite, lim+∞ f = 0.
3. Le formulaire est en page 7
Corrigé n˚ 2.4.4 de l’exercice 4 (transformée d’une gaussienne)
R
R
2
1. fˆ(ω) = R e−iωt e−at dt = R h(ω, t) dt.
La fonction h admet des dérivées partielles par rapport à ω à tous les ordres qui sont
∂k
2
h(ω, t) = (−i t)k e−iωt e−at = hk (ω, t).
k
∂ω
Pour tout k ≥ 0, on a
– t → hk (ω, t) intégrable sur R;
– t → hk (ω, t) continue sur [−A, A];
2
– il existe φk , intégrable sur R telle que pour (ω, t) ∈ [−A, A]×R, |hk (ω, t)| = |t|k e−at ≤
2
Ak e−at = φk (t).
R
Ainsi, fˆ est de classe C ∞ sur tout [−A, A] et donc sur R et fˆ0 (ω) = R h1 (ω, t) dt.
2. En dérivant, nous avons
Z
Z 2
fˆ0 (ω) =
h1 (ω, t) dt = −i
te−at e−iωt dt
R
"R 2
#∞
Z −at2
e−at −iω t
e
= −i
e
+i
(−iω)e−iωt dt
−2a
R −2a
−∞
ω
ce qui nous donne l’équation différentielle y 0 (ω) = −
y(ω) assortie de la condition
−2a
R −at2
initiale y(0) = fˆ(0) = R e
dt, d’où, après calcul de l’intégrale de Gauss
r
π − ω2
ˆ
f (ω) =
e 4a
(2.3)
a
10
Corrigé n˚ 2.4.5 de l’exercice 5 (convolution)
1.
2.
Corrigé n˚ 2.4.6 de l’exercice 6 (convolution par une gaussienne)
1.
2.
1
11
3 Étude du signal
1
Pour l’étude du signal on préférera la définition
F(f )(ν) = fˆ(ν) =
Z
e−2πiνx f (x) dx
R
et on se référera donc au deuxième formulaire de la page 7.
3.1
Première approche du spectre
On propose ici une première exploration de la notion de spectre d’un signal. On aborde la question avec le calcul et l’observation graphique des transformées de Fourier de fonctions à support
compact et sommes de sinusoı̈des. La notion de spectre est alors simple et intuitive.
12
Exercice 9
2
ce que la transformée de Fourier nous apprend des fréquences...
1. On considère λ > 0, a > 0 et la fonction à support compact
cos(2πλt), si t ∈ [−a, a]
fa,λ (t) =
0
sinon .
(a) Calculer sa transformée de Fourier ; on choisira la définition
Z
ˆ
e−2πiωx f (x) dx.
F(f )(ω) = f (ω) =
R
(b) Calculer puis représenter graphiquement dans une feuille Maple ou Scilab les transformées de Fourier de fa,λ pour λ = 0.7 et a = 10, 20, 50, par exemple. Qu’observet-on ? Justifier les observations.
(c) La fonction fd
a,λ est elle intégrable sur R? La fonction
2iπωx
ω → fd
a,λ (ω) e
admet elle une intégrale impropre sur R? Le cas échéant, calculez la transformée de
Fourier inverse de fd
a,λ .
2.
(a) On définit
g := (a, t) → (2 cos (2 π λ1 t) + 5 cos (2 π λ2 t) + cos (2 π λ3 t)) H (a − t) H (a + t)
où H désigne la fonction échelon (H(a − t)H(a + 1) est égale à 1 sur [−a, a] nulle
en dehors). Représenter graphiquement la transformée de Fourier de t → g(a, t) sous
MAPLE ou Scilab. Qu’observe-t-on ? Que se passerait il avec l’autre définition usuelle
de la transformée de Fourier ?
(b) On pose φN (t) = fN,λ (t).
Soit ε > 0. Montrer que la suite des transformées de Fourier
φc
N
N
, est uni-
formément bornée sur le complémentaire de ]λ − ε, λ + ε[∪] − λ − ε, −λ + ε[.
c
Évaluez φc
N (±λ). Que dire du comportement asymptotique de (φN (±λ))N ?
Même question pour la transformée d’une combinaison linéaire
X
Ai cos (2 π λi t) .
i
corrigé en 3.3.1
Exercice 10 Le spectre permet il de reconstituer le signal ?
On considère une fonction g continue par morceaux définie sur l’intervalle [0, T ], T > 0.
1. Montrer que les coefficients de Fourier cn (g̃) n ∈ Z, de son prolongement T −périodique
s’expriment simplement en fonction de la transformée de Fourier de g.
2. fichier : AnTransFourierAnSignal.mws
13
2. Montrer que l’on peut exprimer g sous la forme d’une série de fonctions dont les termes
n
dépendent des seules ( !) valeurs de son spectre, échantillonnées en les points , n ∈ Z.
T
Préciser le type de convergence de cette série.
3. Exemple calculatoire et graphique : choisir un signal triangulaire sur [0, T ] arbitraire, calculer formellement sa transformée de Fourier, écrire un bref programme donnant les sommes
partielles de la série de la question précédente et comparer au signal d’origine.
corrigé en 3.3.2
14
Bilan provisoire
1
• La transformée de Fourier d’un signal sinusoı̈dal de durée finie 2a et de fréquence est un sinus
λ
cardinal dont les maxima (en valeur absolue) sont atteints en ±λ et vérifient :
fd
λ,a (±λ)
∼
a→+∞
a.
40
30
20
10
–3
–2
–1
0
1
2
3
w
Transformée de Fourier d’un signal sinusoı̈dal de durée de vie finie
• La transformée de Fourier d’une combinaison linéaire de signaux sinusoı̈daux de durée finie
permet donc de repérer les différentes fréquences qui interviennent dans ce signal.
15
10
5
–8
–6
–4
–2
0
2
4
6
8
w
–5
–10
–15
Transformée de Fourier d’un signal composite de durée de vie finie
• La transformée de Fourier d’un signal composite, continu par morceaux , de durée finie, permet
de reconstituer le signal d’origine comme somme d’une série de Fourier dont les coefficients sont
obtenus par échantillonnage du spectre et qui converge en moyenne quadratique (et converge
normalement si g est de plus de classe C 1 par morceaux et continue) :
2iπnt
∞
1 X n
g̃(t) =
ĝ
e T
T n=−∞
T
(3.1)
La série converge simplement si g est de classe C 1 par morceaux et alors :
2iπnt
∞
g̃(t+) + g̃(t−)
1 X n
=
ĝ
e T
2
T n=−∞
T
15
(3.2)
3.2 Échantillonnage, théorème de Shannon et Nyquist
• Vocabulaire
Nous dirons qu’un signal est analogique lorsqu’il est modélisé par une grandeur susceptible de
prendre des valeurs appartenant à un intervalle de R (on parle d’ensemble continu). C’est le cas
lorsque le signal est produit par des composants électroniques, des instruments de mesure...
Nous dirons qu’un signal est numérique lorsque l’ensemble des valeurs qu’il est susceptible de
prendre est discret (chaque élément de l’ensemble est isolé des autres 3 ) . C’est le cas d’un signal
obtenu par échantillonnage ou produit par un algorithme...
• En pratique, lorsqu’on veut traiter un signal, par exemple déterminer son spectre, on est le plus
souvent amené à le discrétiser, à l’échantillonner. C’est à dire que l’on mesure des valeurs en
des instants le plus souvent régulièrement espacés t0 + kT. On parle alors d’échantillonnage de
1
fréquence fe = . Ce que nous proposons ici, c’est d’étudier les problèmes que pose l’échantillonT
nage et de présenter quelques solutions.
Exercice 11 repliement de spectre
1. On considère ici deux signaux sinusoı̈daux f (t) = cos(ω1 t) et g(t) = cos(ω2 t) avec
|ω1 | =
6 |ω1 |. On mesure ces signaux à des instants k T, k ∈ N.
1
les valeurs f (k T ) et g (k T ) seront elles confonPour quelles valeurs de T ou de fe =
T
dues pour tout k ∈ Z?
ωi
Vérifier que si fe > |f1 | + |f2 |, avec fi =
, ce phénomène ne se produit pas.
2π
Construire une figure comme celle-ci où deux signaux de fréquences f1 et f2 sont indiscernables par un échantillonnage de fréquence fe .
2. De telles valeurs existent elles lorsque f (t) = cos(ω1 t + φ1 ) et g(t) = cos(ω2 t + φ2 )?
corrigé en 3.3.3
3. a est isolé dans D ssi ∃r > 0, ]a − r, a + r[∩D = {a}.
16
Exercice 12 un premier pas vers le théorème de Shannon (échantillonnage)
P
iωk t , un échantillonnage
On veut montrer que pour tout signal signal g de la forme g(t) = N
k=0 ak e
1
ωk
de fréquence fe = > 2 sup |fk | (avec fk =
) permet de reconstituer le signal g comme
T
2π
somme d’une série de fonctions.
sin x
.
x
1. Soient λ ∈ R,, T > 0 et fλ la fonction 2π−périodique qui coı̈ncide avec x → eiλx sur
[−π, π[.
On notera s le prolongement à R de la fonction x →
(a) Donner une représentation graphique de Im(fλ ) sur [−5π, 5π] pour λ = 1/2, λ =
9/10...
(b) Calculer la série de Fourier de fλ et étudier avec soin sa convergence. On précisera en
particulier la somme en les points x = 0, π.
1
En déduire une expression de
puis de de cotan x comme somme d’une série.
sin x
2. Soient (a0 , a1 , ..., aN ) une suite de complexes et (ω0 , ω1 , ..., ωN ) une suite de réels.
(a) Montrer que la fonction g définie par
g(t) =
N
X
ak eiωk t
k=0
vérifie la formule d’interpolation :
X πt
− nπ g(nT )
g(t) =
s
T
n∈Z
π
1
=
.
|ωk |
2|fk |
On commencera par établir le résultat pour un terme eiωk t ...
dès que T est inférieur à la fréquence de Nyquist : T ≤
corrigé en 3.3.4
17
Bilan provisoire
P
iωk t ,
• Fréquence de Nyquist : on appelle fréquence de Nyquist d’un signal g(t) = N
k=0 ak e
1
ωk
la fréquence f0 =
= 2 sup |fk | (avec fk =
). C’est la fréquence d’échantillonnage au
T0
2π
dessus de laquelle il est théoriquement possible de reconstituer un signal comme somme d’une
série (exercice 12), en dessous de laquelle ce n’est en général pas possible (exercice 11).
Nous généraliserons cette définition en section ??.
• Théorème de Shannon, version élémentaire
Théorème 2
P
iωk t . Pour toute fréquence supérieure à la fréquence
Soit g un signal de la forme g(t) = N
k=0 ak e
1
de Nyquist, fe = > 2 sup |fk |, g vérifie la formule d’interpolation de Shannon :
T
X πt
g(t) =
s
− nπ g(nT )
(3.3)
T
n∈Z
Démonstration avec l’exercice 12.
Nous généraliserons ce théorème en section ??.
• Pourquoi parle-t-on de temps et de fréquence ?
Rapprochons la formule 3.2 de la formule de Shannon 3.3
2iπnt
∞
1 X n
g̃(t) =
ĝ
e T
T n=−∞
T
X πt
g(t) =
s
− nπ g(nT )
T
n∈Z
Celle de gauche reconstitue le signal à partir d’un échantillonnage des fréquences, celle de droite,
à partir d’un échantillonnage des temps.
• Que faire en pratique ?
Bien évidemment dans des conditions expérimentales ou numériques, il n’est ni possible d’échantillonner
sur une infinité de termes, ni de sommer sur une infinité de termes dans la somme.
En pratique, on filtre le signal pour éliminer les fréquences dépassant un certain seuil. Connaissant
ce seuil (par exemple 22000Hz) on dispose d’une fréquence de Nyquist (44000 Hz, dans l’exemple
précédent).
On choisit alors une fréquence d’échantillonnage supérieure à la fréquence de Nyquist et on
prélève un nombre d’échantillons aux instants 0,T,2T,...,nT, suffisant pour couvrir l’intervalle. On
calcule alors la somme partielle construite sur ces échantillons. La feuille de travail Maple qui suit
illustre cela.
18
La formule de Shannon
>
(2.1)
On se propose de reconstituer un signal à partir d'un échantillonnage. On construit une
fonction g pour tester la convergence de la formule de Shannon.
>
(2.2)
fny
est la fréquence de Nyquist de la fonction g.
ech(g,f,n) permet d'échantillonner 2n+1 valeurs de la fonction g aux instants 0, T, 2T, ..., 2nT, à
la fréquence f=1/T.
sha(e,f,t) calcule la somme partielle de la série de Shannon construite sur l'échantillon e
(comportant un nombre impair de termes), de fréquence f en t;
>
(2.3)
On teste avec la fréquence de Nyquist et en faisant varier la taille de l'échantillon
>
19
>
On teste à une fréquence supérieure à la fréquence de Nyquist et en faisant varier la taille de
l'échantillon. Plus la fréquence est grande, plus l'échantillon doit être de grande taille pour
couvrir un intervalle de temps donné
>
>
20
>
>
21
3.3
Corrigés des exercices des sections (3.1) et (3.2)
Corrigé n˚ 3.3.1 de l’exercice 9
1.
(a) Le calcul donne (penser à écrire la fonction cos comme somme d’exponentielles) :
pour ω 6= ±λ :
"
#a
Z a
e2iπ(λ−ω)x
e−2iπ(λ+ω)x
−2 iπ ω x
cos (2 π λ x) e
dx =
−
4iπ(λ − ω)
4iπ(ω + λ)
−a
a
sin (2π(λ − ω)a) sin (2π(λ + ω)a)
=
+
2π(λ − ω)
2π(λ + ω)
et pour ω = ±λ, par parité et passage à la limite, une transformée de Fourier étant
toujours continue,
Z a
sin (4πλa)
cos (2 π λ x) e−2 iπ λ x dx =a +
4πλ
−a
(b) Comme la fonction signal est paire, sa transformée de Fourier est réelle. La question
avait donc un sens et voila ce que l’on trouve en superposant les transformées de
Fourier pour λ = 0.7 et a = 10, 20, 40 :
Il apparait que
– les oscillations sont très amorties en dehors d’un voisinage de ±λ = ±0.7 qui est la
fréquence du signal ;
22
– plus la plage d’observation du signal est large (2a), plus l’amplitude maximale atteinte en ±λ est grande avec comme on peut l’observer graphiquement ou le déduire
du calcul précédent,
fˆ(±λ) ∼ a
a→+∞
(c)
2.
(a) Le calcul est sans mystère, simple superposition qui donne avec a = 100 et
g : x 7→ 1 + 2 cos (1.0 π x) − 3 cos (3.0 π x) + 4 cos (6 π x) ,
On retrouve ici les fréquences et des amplitudes maximales aux points ω = ±λi qui
sont proportionnelles aux coefficients des signaux de période λi .
sin x
. Nous avons donc en dehors de ±λ,
(b) Notons s la fonction x →
x
sin (2π(λ − ω)a) sin (2π(λ + ω)a)
fd
+
= a (s(2π(λ − ω)a) + s(2π(λ + ω)a))
λ,a (ω) =
2π(λ − ω)
2π(λ + ω)
23
• Comme pour |x| =
6 0, |s(x)| ≤
1
, nous avons
|x|
d
fλ,a (ω) ≤
1
1
+
4π|λ − ω| 4π|λ + ω|
et si ω ∈]λ
/ − ε, λ + ε[∪] − λ − ε, −λ + ε[, ie : |λ − ω| ≥ ε et |λ + ω| ≥ ε cela donne
1
d
.
fλ,a (ω) ≤
2πε
• Nous avons déjà noté que pour un signal sinusoı̈dal
fˆ(±λ)
∼
a→+∞
a
il vient donc
fd
λ,N (±λ)
∼
N →+∞
N.
• Cela explique l’observation faite avec le signal composite : les valeurs en les points ±λi
sont équivalentes à Ai × a où Ai est l’amplitude de la composante Ai cos(λi t)
Corrigé n˚ 3.3.2 de l’exercice 10
1. Les coefficients de Fourier de g̃ sont donnés par la formule :
Z
Z
1
1
1
1 T −2inπ/T
−2inπ/T
e
g(t) dt =
e
g(t) dt = ĝ
.
cn (g̃) =
T 0
T R
T
T
2. Comme g̃ est continue par morceaux , sa série de Fourier converge en moyenne quadratique
vers g̃ sur R. On écrit donc
2iπnt
2iπnt
∞
1 X n
g̃(t) =
cn (g̃)e T =
e T
ĝ
T n=−∞
T
n=−∞
∞
X
– Lorsque g est continue par morceaux , la série converge en moyenne quadratique vers g
sur [0, T ];
– lorsque g est de classe C 1 par morceaux, elle converge simplement vers la régularisée de
g(0+) + g(T −)
g (en 0 et T la limite est donc
;
2
1
– lorsque g est continue et de classe C par morceaux, elle converge normalement vers g.
3. Exploration numérique et graphique :
(a) On considère par exemple un signal triangulaire sur l’intervalle de temps [0, T ]. Connaissant son spectre (ou transformée de Fourier) le signal est reconstitué avec la formule
g̃(t) =
∞
1 X n 2iπnt
ĝ
e T
T n=−∞
T
T2 ˆ
1 −1 + 2 e−iπ wT − e−2 iT π w
1
ˆ
g̃(0)
=
, g̃(ω) =
=
(1 − cos πωT )e−iπωT
2
2
4
4
π w
2 πω 2
24
restart;
with(plots);
g :=(T, x)-> Heaviside(x)*Heaviside((1/2)*T-x)*x
+ Heaviside(T-x)*Heaviside(x-(1/2)*T)*(T-x);
G := plot(g(2, x), x = -2 .. 4, color = red, thickness = 3);
Fg := (T, w)-> Int(exp(-(2*I)*Pi*w*t)*t, t = 0 .. (1/2)*T)
+ Int(exp(-(2*I)*Pi*w*t)*(T-t), t = (1/2)*T .. T);
c0 := simplify(value(Fg(T, 0));
simplify(value(Fg(T, w)));
1/4 T 2
−1 + 2 e−iπ wT − e−2 iT π w
π 2 w2
Tg :=(N, T, t) -> simplify(
((1/4)*Tˆ2
+sum(value(Fg(T, k/T))*exp((2*I)*Pi*k*t/T), k = 1.. N)
+sum(value(Fg(T, -k/T))*exp(-(2*I)*Pi*k*t/T), k = 1.. N))/T);
Tg(5, T, t);
1
πt
πt
πt
2
T −1800 cos 2
− 200 cos 6
− 72 cos 10
+ 225 π π −2
900
T
T
T
1/4
25
plot({Im(simplify(value(Fg(2, w)))), Re(simplify(value(Fg(2, w))))},
w = -3 .. 3, thickness = 2)
plot(evalc(Tg(7, 2, t)), t = -2 .. 4, color = black, thickness = 2);
display({%, G});
26
Corrigé n˚ 3.3.3 de l’exercice 11
1. Dire que pour tout k ∈ Z, f (kT ) = g(kT ) c’est dire que


∃nk ∈ Z, ω1 k T = ω2 k T + 2nk π
∀k ∈ Z, ou


∃nk ∈ Z, ω1 k T = −ω2 k T + 2nk π



∃nk ∈ Z, k T



∀k ∈ Z, ou




∃nk ∈ Z, k T
=
nk
2 nk π
=
ω1 − ω2
f1 − f2
=
2 nk π
nk
=
ω1 + ω2
f1 + f2
1
n1
=
et on vérifie sans peine que de telles valeurs
fe
f1 ± f2
de T ou de la fréquence d’échantillonnage conduiront à des valeurs identiques de f (kT ) =
g(kT ) (on aura alors des solutions avec nk = kn1 ).
Pour être certain que ce phénomène ne se produira pas on choisira
En faisant k = 1, il vient T =
T <
1
1
≤
soit fe > |f1 | + |f2 |
|f1 | + |f2 |
|f1 ± f2 |
Pour la figure on a choisi f 1 := 1.1; f 2 := 3.06; T =
1
2
=
;
fe
f1 + f2
2. Dans ce cas , dire que pour tout k ∈ Z, f (kT ) = g(kT ) c’est dire que


∃nk ∈ Z, (ω1 − ω2 )k T + (φ1 − φ2 ) = 2nk π
∀k ∈ Z, ou


∃nk ∈ Z, (ω1 + ω2 )k T + (φ1 + φ2 ) = 2nk π
En faisant k = 0 on voit qu’il est nécessaire que φ1 − φ2 = 2n0 π ou φ1 + φ2 = 2n0 π
Chacune de ces conditions est suffisante et conduit à une condition sur T qui s’exprime
n1 − n0
.
T =
f1 ± f2
27
Corrigé n˚ 3.3.4 correction de l’exercice 12
1.
(a) Figure avec
Rf :=(L,x) − > sin(L*x) :
plot(Rf(0.9,x-2*Pi*floor(x/(2*Pi)+1/2) ),x=-5*Pi..5*Pi, discont=true) ;
F IGURE 1 – Im(f0.9 )
(b) si λ ∈ Z, c’est clair : tous les coefficients cn (f ) sont nuls sauf cλ et la fonction est sa
propre somme de Fourier dès que n ≥ |λ|.
Sinon
"
#π
Z 2π
1
1
ei(λ−n)
iλx −inx
cn (f ) =
e e
dx =
= s((λ − n)π).
2π 0
2π i(λ − n)
−π
La fonction fλ est de période 2π et de classe C 1 par morceaux. Sa série de Fourier
converge simplement vers sa régularisée. Cela donne
fλ (x) =
∞
X
s((λ − n)π)einx , x ∈ R/π(2Z + 1),
n=−∞
∞
X
sin(λπ)
fλ (x+) + fλ (x−)
= cos(λπ) =
, si x = (2k + 1)π.
2
(λ − n)π
n=−∞
De cela on déduit avec y = λπ :
cot(y) =
X
1
y
+2
.
y
y 2 − n2 π 2
n≥1
2. Considérons la fonction de période 2π dont la restriction à ] − π, π], est f (t) = eiλx . Sa
série de Fourier en un point de continuité x ∈] − π, π[ est :
X
eiλx =
s((λ − n)π)einx
(3.4)
n∈Z
On en déduit la formule de Shannon en posant x = ωk T et λ = t/T :
28
eiλx = eiωk t =
X
s((t/T − n)π)einωk T .
n∈Z
Ou encore
g(t) =
X
s((t/T − n)π)g(nT ).
n∈Z
On prendra garde au fait que |x| = |ωk |T impose T ≤
π
1
=
.
|ωk |
2|fk |
Enfin, lorsque g est comme dans l’énoncé combinaison linéaire de fonctions de ce type, il
π
vient encore, pour tout T tel que 0 < T <
, et pour t ∈ R, la formule de Shannon
sup |ωk |
(convergence simple)
X
g(t) =
s((t/T − n)π)g(nT ).
n∈Z
29
4
Transformations de Fourier discrètes
Se pose maintenant la question du calcul effectif de la transformée de Fourier d’un signal temporel
(ou fonction) f à partir d’un échantillonnage. Nous allons voir comment une approximation par
la méthode des rectangles fait apparaı̂tre naturellement la transformée de Fourier discrète. Cette
relation entre intégrale de Fourier et transformée de Fourier serait de peu d’intérêt si nous ne
disposions pas d’algorithmes de calculs rapides pour cette dernière. Nous donnons un aperçu de
ces méthodes en indiquant leur utilisation sous Maple et Scilab.
4.1
Transformation de Fourier discrète
Se pose maintenant la question du calcul effectif de la transformée de Fourier d’un signal temporel
représenté par une fonction f. Quelques observations préalables s’imposent.
– Le signal f est observé sur une période restreinte, [T0 , T1 ].
– On ne dispose en pratique que d’un échantillonnage de la fonction f entre des instants T0 et T1 .
C’est à dire que l’on connaı̂t seulement des valeurs f (tk ) pour t0 = T0 , ..., tN = T1 .
– On veut calculer
Z T1
Z ∞
f (t) e−2 i π νt dt ≈
f (t) e−2 i π ν t dt.
fˆ(ν) =
−∞
T0
– On obtient une valeur approchée de cette dernière intégrale avec une somme de Riemann :
Z
T1
f (t)e
−2 i π ν t
N
−1
X
(tk+1 − tk )f (tk )e−2 i π tk ν
dt ≈
T0
k=0
– Comme on aura choisi une subdivision régulière pour échantillonner, tk+1 − tk =
vient :
Z T1
N −1
k ∆T
∆ T −2 i π ν T0 X
f (t) e−2 i π ν t dt ≈
e
f (tk ) e−2 i π ν N
N
T0
T1 − T0
, il
N
(4.1)
k=0
– Il va de soi que l’on ne peut pas calculer toutes les valeurs fˆ(ν). Il va falloir choisir. On en
calculera N. On va vite voir lesquelles, pourquoi et comment après avoir défini la notion de
transformation de Fourier discrète.
Définition 3 Soit (αk )0≤k≤N −1 une suite de N nombres. On appelle transformée de Fourier
discrète de cette suite, la suite finie de N nombres également, définie par
T F D(α)(n) =
N
−1
X
−2iπ
αk e
kn
N , n = 0, 1, ..., N − 1.
(4.2)
k=0
Exercice 13
1. Expliciter la matrice ΩN de la transformation de Fourier discrète à l’aide de ωN
1
et montrer que l’inverse de ΩN est ΩN .
N
30
−2iπ
=e N .
2. Calculer la transformée de Fourier discrète de α = (f (tk ))0≤k≤N −1 , la subdivision étant
définie comme ci-dessus et établir une relation asymptotique entre T F D(α)(n) et fˆ(ν)
pour ν bien choisi.
voir corrigé en 4.5.1
4.2
Les transformations de Fourier rapides
Le calcul de la TFD pourrait s’avérer laborieux pour de grandes valeurs de N : le produit matricevecteur brutalement programmé demande O(N 2 ) opérations. En 1965, redécouvrant une méthode
déjà mise en œuvre par Gauss, Cooley et Tuckey, ont remarqué que l’on pouvait réduire ce calcul
à O(N ln(N )) opérations... Depuis, de nombreuses variantes de cet algorithme, connues sous le
nom de transformées de Fourier rapides (TFR ou FFT), ont vu le jour, les plus récentes visant à
s’adapter à l’architecture des processeurs pour réduire encore les temps de calculs 4 et limiter les
erreurs lors des calculs en flottants.
• Sans algorithme rapide, le calcul des N sommes
Y [n] =
N
−1
X
nk
X[k]ωN
k=0
demande O(N 2 ) opérations comme nous l’avons déjà remarqué ; il y a bien sûr plusieurs façons
de l’envisager et nous n’insisterons pas sur cette question autrement que par l’exercice 14 qui suit :
Exercice 14 TFD sans algorithme rapide
1. Dans l’exemple qui suit, où le même algorithme est proposé en Maple et en Scilab, dénombrer
les additions et les multiplications (on prendra garde au décalage des indices, ici Y =
[Y [1], ..., Y [N ]] est indexé de 1 à N et non pas de 0 à N-1 et non pas [Y [0], ..., Y [N − 1]]
comme dans nos formules).
2. Justifier que l’on calcule bien la transformée de Fourier discrète du vecteur X.
3. Y aurait il inconvénient à construire la matrice ΩN pour ensuite effectuer le calcul Y =
ΩN X? Comparer les opérations.
4. à distinguer du nombre d’opérations
31
Avec Maple
N
X
w
#
:=
:=
:=
au
8;
Array(1 .. N, i-->x[i-1]);
z;
lieu de exp(-(2*I)*Pi/N) pour vérification rapide
wn := 1;
Y := Array(1 .. N);
for n to N do
Y[n] := 0;
wnk := 1;
for k to N do
Y[n] := Y[n]+wnk*X[k];
wnk := wnk*wn
end do;
wn := wn*w
end do;
Y(7);
x0 + z 6 x1 + z 12 x2 + z 18 x3 + z 24 x4 + z 30 x5 + z 36 x6 + z 42 x7
Avec N :=5 et w :=exp(-(2*I)*Pi/N) en place de w :=z, nous obtenons pour Y(4),
3
6
9
12
x0 + e−2/5 iπ x1 + e−2/5 iπ x2 + e−2/5 iπ x3 + e−2/5 iπ
x4
Comme toujours, la vérification est immédiate avec un logiciel de calcul formel ; on voit ce que l’on
a programmé : c’est l’avenir !
voir corrigé en (4.5.2)
32
• Diviser pour régner
Supposons que N soit un entier de la forme N = N1 N2 et reprenons le calcul de
Y [n] =
N
−1
X
nk
X[k]ωN
k=0
en posant pour n, k tels que 0 ≤ n, k ≤ N − 1 :
(
n = n1 + n2 N2 , 0 ≤ n1 ≤ N2 − 1
k = k1 + k2 N1 , 0 ≤ k1 ≤ N1 − 1
ce qui nous donne
k1 n1 k1 n2 k2 n1
kn
+
+ k2 n2
=
+
N
N
N1
N2
k1 n1
k1 n2
k2 n1
nk
wN
= wN
× wN
× wN
×1
1
2
Y [n1 + n2 N2 ] =
NX
1 −1

NX
2 −1

k1 =0


k2 n1  k1 n1  n2 k1
X[k1 + k2 N1 ]wN
wN
ωN1
2
(4.3)
k2 =0
Ainsi, pour chaque valeur de n (ou chaque couple (n1 , n2 ), nous avons à calculer N1 TFD de taille
N1 :


NX
2 −1
k2 n1 

X[k1 + k2 N1 ]wN
2
k2 =0
k1 n1
à multiplier ce résultat par wN
et à calculer la TFD de taille N2 du vecteur obtenu. Cela conduit
à une procédure récursive qui s’écrit plus facilement lorsque N1 est invariablement égal à 2, ce qui
suppose que N est de la forme N = 2p .
Exercice 15 l’algorithme de Cooley et Tuckey pour N = 2p .
N
.
2
1. Réécrire dans un tel cas la formule (4.3) en observant que k1 et n1 ne prennent que deux
valeurs.
On se propose ici d’expliciter l’algorithme de calcul rapide lorsque N1 = 2 et N2 =
2. Écrire un algorithme récursif du calcul de Y lorsque Y est une puissance de 2. Préciser la
condition d’arrêt ou l’initialisation.
3. Déterminez la complexité de votre algorithme en nombre de multiplications par exemples.
S’il n’est pas en O(N ln N ) réécrivez le !
4. Réécrire cet algorithme dans une procédure itérative avec une boucle while.
voir corrigé en (4.5.3)
33
Bilan
• Soit N un entier naturel non nul. La transformée de Fourier discrète (TFD ou DFT) d’ordre N
définie par
!
N
−1
X
nk
(Xk )k ∈ CN → (Yn )n =
Xk ωN
∈ CN ,
k=0
n
est bijective et son inverse (transformée de Fourier discrète inverse) est :
!
N
−1
X
1
−nk
(Yn )n ∈ CN → (Xk )k =
Yn ωN
∈ CN ,
N
n=0
La matrice de la TFD (ou DFT) :

1
1
1
1

 1
ωN
ωN 2
ωN 3


 1
ωN 2
ωN 4
ωN 6

ΩN = 
 1
ωN 3
ωN 6
ωN 9


 ..
..
..
..
 .
.
.
.

1 ωN (N −1) ωN 2(N −1) ωN 3(N −1)
k
...
...
...
...
..
.
...

1

ωN (N −1) 


2(N
−1)

ωN

,
3(N
−1)

ωN



..

.

ωN (N −1)
2
1
ΩN .
d’inverse Ω−1
N =
N
Attention cette représentation matricielle, utile pour les démonstrations, n’intervient pas dans les
calculs effectifs (comme dans la plupart des cas).
n
• Lorsqu’on dispose d’un échantillonnage de taille N de f sur [T0 , T1 ], pour toute valeur ν =
,
∆T
Z
T1
f (t) e
T0
−2 i π ν t
∆ T −2 i π ν T0
e
T F D(α)(ν ∆T ) + O
dt =
N
1
N
(4.4)
• L’intérêt de cette formule c’est qu’on dispose pour le calcul d’algorithmes de complexité O(N ln N )
lorsque N est une puissance de 2. On se ramène à ce cas en complétant (par interpolation le plus
souvent) les données dont on dispose pour se ramener à 2N points. On dispose par ailleurs d’algorithmes rapides plus sophistiqués qui permettent de traiter les relevés en N points sans que N ne
soit une puissance de 2 (voir l’article [2], disponible sur la toile, qui décrit ce type d’algorithme et
en détaille un qui est par ailleurs implémenté sous Scilab - voir l’aide de Scilab version 5.4.0).
34
4.3
Transformées en cosinus discrètes
A voir ; utile pour les images...
La ressemblance avec la section précédente, où de la transformée de Fourier discrète il s’agissait,
n’a rien de fortuit.
On rappelle que l’espace C2π muni du produit scalaire
Z 2π
1
¯ dt
f (t)g(t)
ψ(f, g) =
2π 0
est un espace préhilbertien.
Exercice 16 Vérifier que les trois familles de fonctions
{t → eint /n ∈ Z}
{t → cos(nt)/n ∈ N} ∪ {t → sin(nt)/n ∈ N∗ },
nt
/n ∈ N},
2
sont orthogonales. Préciser les normes de chacune de ces foncions.
{t → cos(
Exercice 17 On considère une fonction de classe C 1 sur l’intervalle [0, T ].
1. coefficients de Fourier
– On considère la fonction g, T −périodique égale à f sur [0, T [. Donner une expression
des coefficients de Fourier de g en fonction de ceux de g 0 .
– Donner une condition suffisante sur f pour que la convergence de cette série de Fourier
soit une convergence normale.
2. coefficients en cosinus...
– Montrer que la fonction h de période 2T, paire, égale à f sur [0, T [, est limite uniforme
de polynômes de la forme
X
a0 /2 +
ak cos(kπt/T ).
k
– Donner une expression des coefficients de Fourier de h en fonction de ceux de h0 . Comparer avec les coefficients de Fourier de g.
3. Peut on en dire autant pour des polynômes 5
X
bk sin(kπt/T )?
k
• Transformée en cosinus discrète en dimension 1
On définit la transformée en cosinus discrète (TCD, en français et DCT, en anglais) comme étant
l’application linéaire de CN dans lui-même définie par :
T CD = p ∈ CN → G ∈ CN ,
5. penser à la façon dont nous avons étudié l’équation de la chaleur
35
Gj
= αcj
n−1
X
k=0
πj
pk cos (2k + 1)
2N
(4.5)
où les coefficients cj sont

c = √1
j
2

cj = 1
sij = 0,
sinon.
Le choix de α sera fixé par la suite. L’exercice qui suit a pour but de déterminer la bijection
réciproque et son expression :
T CDI(G)k0 = pk0
N
−1
X
= β
j=0
πj
cj Gj cos (2k0 + 1)
2N
(4.6)
Exercice 18 Calcul de la transformation inverse
On se propose montrer que T CD est un automorphisme de CN et d’en calculer la bijection réciproque.
La démarche est celle mise en œuvre pour le calcul de l’inverse de la T F D et les résultats forts
ressemblants.
1. Donner une expression de la matrice DN associé à l’automorphisme T CD. 6
2. Soit k0 , un entier fixé, tel que 0 ≤ k0 ≤ N − 1. On se propose de calculer
N
−1
X
j=0
πj
cj Gj cos (2k0 + 1)
.
2N
(a) Exprimer l’expression ci-dessus en fonctions des deux sommes
σ1 =
N
−1
X
c2j cos
j=0
σ2 =
N
−1
X
j=0
πj
(k + k0 + 1)
N
πj
c2j cos (2k0 + 1)
.
N
(b) Calculer chacune d’elles avec soin (en distinguant suivant que k = k0 ou pas, que
k + k0 est pair ou impair.
(c) En déduire que la T CD et sa matrice DN sont inversibles, donner une expression de
chacune d’elles. Que dire des colonnes de DN ?
3. Choisir α pour que les coefficients α et β de la formule (4.6), soient égaux.
4. Calculer
2
2
2 (cos(1/16 π)) + 2 (cos(3/16 π)) + 2
5
cos( π)
16
par exemple et montrer que la matrice DN est orthogonale.
6. fichier MAPLE : TCD.mws
36
2
2
7
+ 2 cos( π) ,
16
4.4
Transformée en cosinus discrète en dimension 2
On introduit ici la TCD et son inverse en dimension 2 et on essaiera de montrer comment elle
intervient pour la compression du signal audio ou vidéo.
On définit de façon analogue la transformée en cosinus discrète pour un signal p(x, y) dépendant
de deux paramètres discrets. Penser pour fixer les idées que la matrice carrée
P = [p(x − 1, y − 1)]1≤x≤N,1≤y≤N ,
formée des coefficients des (p(x, y))0≤x<N,1≤y<N représente une image de N 2 points (ou pixels),
x et y représentant les abscisse et ordonnée de ces points, p(x, y) un décrivant une intensité de
couleur par exemple...
T CD2 (i, j)
=
T CDI2 (x, y)
=
N
−1 N
−1
X
X
πi
πj
1
p(x, y) cos (2x + 1)
ci cj
cos (2y + 1)
(4.7)
2N
2N
2N
x=0 y=0
N −1 N −1
πi
πj
1 X X
√
ci cj T CD2 (i, j) cos (2x + 1)
cos (2y + 1)
(4.8)
2N
2N
2N i=0 j=0
√
avec toujours,

c = √1
j
2

cj = 1
sij = 0,
sinon.
Exercice 19 On note toujours DN la matrice de la TCD en dimension 1, qui est donc une application de Cn dans lui-même.
1. Soit P la matrice carrée de taille N définie par
P = [p(x − 1, y − 1)]1≤x≤N,1≤y≤N .
Exprimer T CD2 (P ) en fonction de P, de DN et de sa transposée.
2. En déduire que T CD2 est inversible et vérifier que son inverse est bien donnée par la formule (4.8).
3. Calculer, lorsque, P est une matrice 1024 × 512, le nombre d’additions, de multiplications
nécessaires pour calculer sa T CD.
4. On partage la matrice en blocs 8 × 8. Combien faut il d’additions, de multiplications pour
calculer le T CD de chaque bloc ? Quel est le coût total ?
37
4.5
Corrigés des exercices de la section (4)
Corrigé n˚ 4.5.1 de l’exercice 13
1. Pour n = 0, 1, ..., N − 1 :
T F D(α)(n) =
N
−1
X
αk e
−2iπ
kn N −1
X
kn
N =
αk ωN
,
k=0
k=0
(n−1)(k−1)
la matrice de cette application linéaire est donc définie par ΩN, n, k = ωN
ΩN

1
1
1
1
...






=






1
ωN
ωN 2
ωN 3
...
1
ωN 2
ωN 4
ωN 6
...
1
ωN 3
ωN 6
ωN 9
...
..
.
..
.
..
.
..
.
..
.
1 ωN (N −1) ωN 2(N −1) ωN 3(N −1) . . .
, soit

1

ωN (N −1) 


2(N
−1)

ωN


3(N
−1)

ωN



..

.

2
ωN (N −1)
C’est la matrice de Vandermonde attachée aux racines N ièmes de l’unité. Elle est donc
inversible. Pour déterminer son inverse, un calcul explicite pour N = 2, 3, 4 permet de
conjecturer. Calculons ΩN ΩN pour vérifier
N
N
X
X
(n−1)(j−1) (j−1)(k−1)
ωN
ωN
ΩN ΩN n,k =
ΩN n j ΩN j k =
j=1
j=1
Comme ω N est le conjugué de ωN , on a
N
X
(j−1)(n−k)
ΩN ΩN n,k =
ωN
=
j=1
(
N
0
n−k
=1
si ωN
sinon (somme des termes d0 une s.g.).
La racine ωN étant primitive (ses puissances engendrent le groupe des racines de l’unité),
seuls les termes diagonaux sont non nuls. On a donc
Ω−1
N =
2. Rappelons que tk = t0 + (T1 − T0 )
0, 1, ..., N − 1 :
T F D(α)(n) =
N
−1
X
1
ΩN .
N
k
k
tk − t0
soit
=
. Il vient donc, pour n =
N
N
T1 − T0
kn
f (tk )ωN
k=0
=
N
−1
X
k=0
38
f (tk )e
−2iπ
kn
N
Si nous rapprochons cela de la formule (4.1)
Z
T1
f (t) e−2 i π ν t dt =
T0
lim
N →+∞
nous avons pour chaque valeur ν =
N
−1
X
T F D(α)(n) =
N −1
k ∆T
∆ T −2 i π ν T0 X
e
f (tk ) e−2 i π ν N
N
!
,
k=0
n
,
∆T
kn
f (tk )ωN
=
k=0
N
−1
X
k=0
f (tk )e
−2iπ
N
kn
=
N
−1
X
f (tk )e−2 i π k ν
∆T
N
k=0
en d’autres termes et en se souvenant que la méthode des rectangles converge en O
Z
T1
−2 i π ν t
f (t) e
T0
∆ T −2 i π ν T0
dt =
e
T F D(α)(ν ∆T ) + O
N
1
N
1
N
:
Corrigé n˚ 4.5.2 de l’exercice 14
1. Le programme MAPLE avec X1 , X2 , ..., Xn et w = ωN = e−2iπ/N donnés :
Le programme
le total des opérations
w := exp(-(2*I)*Pi/N);
wn := 1;
Y := Array(1 .. N);
for n to N do
Y[n] := 0;
wnk := 1;
for k to N do
Y[n] := Y[n]+wnk*X[k];
wnk := wnk*wn
end do;
wn := wn*w
end do
N 2 additions et 2N 2 multiplications dans la boucle
interne ;
N multiplications dans la boucle externe ;
2. Ce que l’on calcule : la preuve formelle est dans le cours d’informatique [1].
3. Coût de la construction de la matrice ΩN :
39
Corrigé n˚ 4.5.3 de l’exercice 15
1. On suppose que N1 = 2. Comme
(
n = n1 + n2 N2 ,
k = k1 + k2 N1 = k1 + 2 k2 ,
cela donne
0 ≤ n1 ≤ N 2 − 1
0 ≤ k1 ≤ 1.
kn
k1 n1 k1 n2 k2 n1
+ k2 n2
=
+
+
N
N
2
N2
k2 n1
k1 n1
k1 n2
nk
×1
× wN
wN
= wN
× wN
2
1
La formule (4.3) se réécrit donc, avec ωN1 = ω2 = −1,
Y [n1 + n2 N2 ] =
NX
1 −1


k1 =0
=
NX
2 −1
1
X

k1 =0
=

k2 n1  k1 n1  n2 k1
X[k1 + k2 N1 ]wN
wN
ωN1
2
(4.9)
k2 =0
NX
2 −1



NX
2 −1


k2 n1  k1 n1 
X[k1 + 2 k2 ]wN
wN
(−1)n2 k1
2
(4.10)
k2 =0


k 2 n1 
X[2 k2 ]wN
+ (−1)n2 
2
k2 =0
NX
2 −1

k2 n1  n1
X[2 k2 + 1]wN
wN
2
k2 =0
(4.11)
On calcule donc pour (n1 , n2 ) fixé, (avec 0 ≤ n1 ≤ N2 − 1, et 0 ≤ n2 ≤ 1), deux TFD de
N
taille N2 = , celles des composantes paires et impaires de X.
2
2. Voir l’algorithme pages suivantes. Le programme Scilab est détaillé dans [1]
3. On notera Cp le nombre de multiplications lors d’un appel avec N = 2p du programme
récursif en MAPLE qui est présenté dans le tableau.
(
C0 =
0
N
On a, puisque N2 =
= 2p−1 ,
2
Cp+1 = 2 Cp + 3 × 2p
Ce qui donne Cp = 3 × p × 2p−1 que l’on exprime encore en fonction de N :
40
3
× ln N × N
2 ln 2
Le programme avec MAPLE
TFDrapide := proc (X)
local N, N2, r, Y, Y0, Y1, X0, X1, w, W, n1, n2, k2;
N := Dimension(X);
if N = 1
then
Y := X;
else
N2 := iquo(N, 2, r);
if r <> 0 then error "2 doit diviser N" end if;
X0 := Array([seq(X[2*k2-1], k2 = 1 .. N2)]);
Y0 := TFDrapide(X0);
X1 := Array([seq(X[2*k2], k2 = 1 .. N2)]);
Y1 := TFDrapide(X1);
Y := Array(1 .. N);
w := exp(-(2*I)*Pi/N);
W := 1;
for n1 from 0 to N2-1 do
Y[n1+1]
:= Y0[n1+1]+Y1[n1+1]*W;
Y[n1+1+N2] := Y0[n1+1]-Y1[n1+1]*W;
W
:= W*w
end do
end if;
Y
end proc
Il serait judicieux d’introduire un compteur pour vérifier tout ce qui précède (à placer dans la
boucle for et à déclarer en variable globale sans oublier de l’initialiser avant chaque appel).
41
5
L’analyse du son
On se propose d’aborder l’étude des sons du point de vue de l’analyse du signal. Pour cela nous
allons enregistrer quelques sons et montrer comment fonctionnent les outils que nous avons mis
en place.
• Quelques remarques préliminaires et désordonnées
– Un son est détecté par l’intermédiaire d’une membrane (située dans un micro ou dans votre
oreille) qui enregistre des variations de pression de l’air ou du fluide dans lequel elle est immergée. Sur le plan physique enregistrer un son consiste donc à fournir des relevés de pression
sous une forme analogique t → P (t) ou numérique (P (k T ))k .
– Seules les variations rapides de pression semblent détectées par nos instruments (un plongeur
n’entend pas de bruit lié directement à la variation de pression accompagnant sa descente, un
micro placé dans un caisson hyperbare ne nous permettra pas non plus de repérer la pression
ambiante contrairement à un manomètre) ;
– un son est produit et a fortiori enregistré pendant un intervalle de temps fini [τ0 , τ1 ]; on lui
associera donc au mieux des restrictions ou troncatures de signaux virtuels périodiques (cas du
diapason) ;
– en dehors de cas très particuliers, mais dont la reconnaissance va se révéler fondamentale, un
signal sonore n’est en général pas périodique.
– On distingue pourtant des sons graves et aigus ce qui conduit aussi à penser aussi le son en
termes de fréquences.
– Nous nous proposons d’étudier un enregistrement sonore. Cela nous conduira à numériser (ou
échantillonner) notre enregistrement ; cela signifie que l’intervalle d’étude [τ0 , τ1 ] est discrétisé
et que la pression est relevée avec une fréquence proportionnelle au nombre de termes de la subdivision, ce qui induit une perte d’information. Ce phénomène a déjà été étudié dans l’abstrait
(fréquence de Nyquist, théorème de Shannon page18).
• Les outils
On pourra procéder à un enregistrement sonore avec le logiciel livré avec la carte son de l’ordinateur, utiliser un logiciel comme goldwave, télécharger un fichier son quelconque... Pour l’analyse
du son nous allons proposer de visualiser et analyser avec goldwave, Maple et Scilab.
5.1
Enregistrement aux formats .wav et .text avec goldwave
Commençons par enregistrer un son quelconque sous goldwave. On peut enregistrer le son sous le
format .wav dans un fichier .wav ou .txt. La première ligne du fichier .txt est
[48000Hz, Channels: 2, Samples: 343680, Flags: 0]
ce qui signifie que la fréquence d’échantillonnage est de 48000Hz, que le son a deux canaux
(stéréo) et que l’on a prélevé 343680 échantillons, ce qui couvre environ 7 secondes. On présente
ici, la première figure, un relevé des variations de pression enregistrées sur une durée d’environ
4.5 secondes (le texte, remarquable, est pouet, pouet, pouet, pouet, I stop it). On reconnaitra la
prononciation du I : a-i. Comme on l’a déjà observé, à part la répétition de pouet pouet, ce n’est
pas vraiment périodique...
42
On reprend donc le même fichier et on observe le deuxième pouet pendant 0.4 s puis pendant
0.045s et là, on commence à y voir plus clair.
43
44
Références
[1] T HIERRY AUDIBERT, A MAR O USSALAH
Informatique, Programmation et Calcul Scientifique
Ellipses, 2013
[2] M ATTEO F RIGO AND S TEVEN G. J OHNSON
The Design and Implementation of FFTW3 rapide
Proc. IEEE, vol. 93, no. 2, pp. 216-231 (2005).
disponible en pdf sur le site de la revue IEEE :
http ://ieeexplore.ieee.org/search/
http ://www.fftw.org/
[3] C. G ASQUET, P. W ITOMSKI .
Analyse de Fourier et Applications
Masson, 1990
[4] H. R EINHARD .
Elements de mathématiques du signal
tome 1 - Signaux déterministes
Dunod, 1995
[5] W. RUDIN .
Analyse réelle et complexe
Masson, 1980
[6] B ERNARD S IMON .
La marée 1. La Marée océanique cotière
Institut océanographique Collection : SYNTHESES (2007)
[7] L. S CHWARTZ .
Méthodes mathématiques pour les sciences Physiques
Hermann, 1993
[8] L. S CHWARTZ .
Analyse III - Calcul Intégral
Hermann, 1993
[9] L. S CHWARTZ .
Analyse IV - Applications de la Théorie de la Mesure
Hermann, 1993
[10] W IKIPEDIA .
Transformée de Fourier rapide
http ://fr.wikipedia.org/wiki/Transformée de Fourier rapide
consultée novembre 2012
45
Index
algorithme
Cooley et Tuckey, 34
convolution, 4, 6
dérivation
transformée de Fourier, 3
Fourier
transformée de, 3
transformée discrète, 30
fréquence
de Nyquist, 18
Fubini
la formule de, 5
Lebesgue
lemme de Riemann-Lebesgue, 3
lemme
de Riemann-Lebesgue, 3
Nyquist
fréquence de, 18
produit
de convolution, 4
Riemann
lemme de Riemann-Lebesgue, 3
Shannon
formule de, 17, 18
théorème, 18
TFD, 30
théorème
de Fubini, 5
de Shannon, 17, 18
transformée
de Fourier, 3
de Fourier Discrète, 30
46