Analysis of Discrete Fourier Transform and its Applications

subplot(3,1,1)

plot(real(FS));

subplot(3,1,2);

plot(imag(FS));

subplot(3,1,3);

plot(abs(FS));

1

Les coefficients de la fft peuvent être vu comme les coefficients de la décomposition d’un vecteur

dans la base formée par les exponentielles complexes discrètes. En effet pour tout entier k entre

1 et N on peut construire le vecteur ek de longueur N défini par : ek(n) = e

(ek)1(cid:54)k(cid:54)N forme une base orthogonale pour le produit scalaire défini par

2iπ(k−1)(n−1)

N

. La famille

(cid:104) f , g(cid:105) =

N

∑

k=1

f [k]g[k].

(3)

1√

N

N. Si on

Cette base n’est pas orthonormée dans la mesure où les vecteurs ont tous une norme

dans la transformée directe on aurait une Base Orthonormée (BO). Comme

avait un facteur

tous les vecteurs ont la même norme, on peut faire l’abus de langage consistant à dire que les

exponentielles complexes discrètes forment une BO pour le produit scalaire usuelle. Une autre ruse

consiste à modifier le produit scalaire par un facteur 1

N pour e,n faire une BO.

√

2. Construire pour k = 0, 1, 2, 5, 10, le vecteur ek et visualiser sa transformée de Fourier à l’aide du

programme précédent.

3. Afficher la fft du vecteur suivant :

>>k=7;

>>N=1024;t=[0:1/N:1-1/N];S3=sin(2pit);

et interpréter le résultats.

4. Tester d’autres valeurs de k, remplacer le sinus par un cosinus, introduisez une phase φ et essayer

d’anticiper le résultat. Tester en particulier les valeurs k = 1020, 1021, 1022. Qu’en pensez vous ?

5. Justifier que si ( f [k])k est un vecteur dont les composantes sont réelles, le module de la fft admet

une symétrie. Plus précisément montrer que pour tout k > 0, ˆf [k + 1] = ˆf [N − k + 1].

6. Est ce encore vrai si le vecteur f n’est pas réel ?

2iπk·

N

et eN−k+1 = e− 2iπk·

N

N ) et cos( 2iπk·

L’espace vectoriel engendré par ek+1 = e

est le même que celui engendré

par sin( 2iπk·

N ). Ansi pour tout k > 1, les composantes d’indice k + 1 et N − k + 1 de

la fft peuvent être interprétées comme la composante d’un vecteur sur les sinusoides de fréquence

2πk

N .

On appelle basse fréquences, les valeurs de ˆf [k] pour les petites valeurs de k. Elles correspondent

aux coeffcients de décomposition d’un vecteur f sur des exponentielles discrète complexes qui

oscillent peu. On appelle hautes fréquences, les coeffcients associés aux exponentielles qui oscille

beaucoup. Les exemples précédents indiquent que la fft de matlab dispose les hautes fréquence au

centre et les basses fréquences à gauche et à droite.

Ces concepts de basses fréquences et de hautes fréquences discrètes sont discutables dans la mesure

où chaque coefficients correspond à une infinité de fréquences différentes.

7. A quoi correspond la composante ˆf [1] ?

8. Expliquer ce que renvoie le programme suivant :

function Exo1(S,k)

plot(S);hold on;

FC=0*S;

FS=fft(S);

FC(k+1)=FS(k+1);

FC(N-k+1)=FS(k+1);

C=real(ifft(FC));

plot(C,’r’);hold off;

9. Après l’avoir testé sur différents vecteurs, expliquer ce que renvoie le programme suivant :

function Exo2(S,k)

close all;

FS=fft(S);

for j=1:k

plot(S);hold on;

FApp=FS;

FApp(j+1:N-j+1)=0;

App=real(ifft(Fapp));

plot(App,’r’);hold off;

pause(0.1);

end

10. Que se passe t-il si on utilise cette fonction si S est un vecteur ek ou une sinusoide de la forme

S[n] = sin( 2iπ(k−1)n

+ φ ) ?

N

11. La commande fftshift est une commande qui a pour but de dsiposer les coefficients de la FFT

différamment en mettant au centre les coefficients dits basse fréquence au centre et haute fréquence

aux extrémités. Attention cette command n’effectue pas de FFT elle échange juste des coefficients.

Tester la commande ffshift sur plusieurs exemples.

2

4 FFT et note de musique

La transformée de Fourier est un des principaux outils en traitement du son. Elle permet d’analyser

et de mesurer les différentes fréquences présentes dans un son. Si le son est musical on peut associer

à chaque fréquence une note.

Télécharger sur ma page web un son de saxo, de violon et un mélange de différents instruments. Ces

fichiers sont au format .wav. Pour transformer ces fichiers en vecteur matlab utiliser la commande

suivante :

>>[Y,Fs,B]=wavread(’Mix11.wav’);

La variable Y est un vecteur (parfois une matrice avec deux lignes si le son est enregistré en stéréo),

contenant la position de la membrane d’un haut parleur par rapport à sa position de repos à chaque

échantillon de temps. Y prend des valeurs réelles entre −1 et 1.

La variable Fs contient la fréquence d’échantillonnage du son exprimé en Hertz. C’est cette fré-

quence qui permet de passer d’un son à un vecteur et d’un vecteur à un son. Si Fs = 44000 cela

signifie que le son est échantilloné à 44000 Hz, autrement dit deux valeurs consécutives de Y repré-

44000 secondes.

sentent les positions de la membrane d’un haut parleur à deux instants séparés de

La variable B ne sera pas utilisé et contient le nombre de bits sur lesquels sont codés les compo-

santes de Y .

1

12. Visualiser le vecteur Y dans sa totalité pour chacun des 3 sons.

Calculer la DFT d’un vecteur représentant un son musical sur un temps long n’a pas d’intérêt ma-

jeur. En effet, la transformée de Fourier représente un vecteur comme une somme de sinusoides

discrètes. Quand plusieurs notes sont jouées successivement il est rarement intéressant de repré-

senter le son comme une somme de sinusoides. Ainsi dans la pratique, on découpe de le son en

sous-vecteurs, appelés fenêtres sur lesquels on calcule des DFT.

13. Dans les différents vecteurs issus de son musicaux, extraire des sous vecteurs de taille 1024 ou

2048 et les visualiser. A quelle durée de son ces extraits correspondent-ils ?

14. Visualiser les FFT de ces extraits.

Le premier coefficient de la FFT correspond à la somme des coefficients de f , il mesure donc la

moyenne du vecteur f . Le dernier coefficient correspond à la fréquence d’échantillonnage : par

exemple 44000 Hz. Les autres coefficients correspondent à des fréquences régulièrement espacées

entre 0 et FS. Si on calcule une FFT sur 44001 points, chaque coefficient correspond à un intervalle

fréquenciel de 1 Hz. Ainsi le 101 ième coefficient correspond à une fréquence de 100 Hz.

Si on fait une FFT sur 441 points chaque coefficient correspond à un intervalle de 10 Hz. Ainsi plus

la fenêtre est grande plus on a de précision fréquentielle mais moins on a de précision spatiale.

15. Dans les notes jouées par le saxophone et le violon on remarque plusieurs coefficients de Fourier

importants. A quoi correspondent ils ? Déterminer les notes jouées.

16. A l’aide la commande wavwrite qui s’utilise de manière suivante :

>>wavwrite(Y,FS,’toto.wav’);

générer un petit fichier wav d’une durée de 2 secondes avec une note pure, par exemple le la 440.

C’est à dire à 440 Hz.

5 Filtrage et convolution

Si f et g sont deux vecteurs de longueur N, la convolution circulaire de f et g est le vecteur h = f (cid:63)g

de longuer N défini par :

h[n]

N

∑

k=1

f [k]g[n − k] =

N

∑

k=1

f [n − k]g[k]

(4)

où f et g sont considérés comme des vecteurs périodiques de période N, c’est à dire tels que

f [l] = f [l + N] et g[l] = g[l + N], ∀l ∈ Z. Un calcul élémentaire montre que ∀k ∈ Z, ˆh[k] = ˆf [k] ˆg[k].

17. Ecrire une fonction matlab

function h=ConvFourier(f,g)

qui calcule la convolution de deux vecteurs en utilisant la FFT.

18. Faites plusieurs tests et afficher les résultats.

19. Quel est le résultat de la convolution entre deux vecteurs ek et el ?

20. Déterminer l’élément neutre pour la convolution discrète circulaire.

21. Quel est le résultat de la convolution entre un vecteur S de taille N et un vecteur de la forme

sin(2kπ/N) ?

On construit une gaussienne discrète à l’aide la fonction suivante :

3

function G=Gaussienne(N,s)

t=[0:1/(N-1):1]*2-1;

G=exp(-t.^2/s);

G=G/sum(G);

22. Afficher plusieurs gaussiennes en testant différentes valeurs de s.

23. Pour différentes valeurs de s afficher également en parallèle le module de la FFT de la gaussienne

en centrant les basses fréquences (en utilisant fftshif.)

24. A quoi ressemble le module de la FFT d’une gaussienne discrète ?

Dans la plupart des convolution, on réalise la convolution entre un objet d’intérêt (Signal ou Image)

avec un vecteur lié à un effet (Flou, dérivation, moyennage local). Ce vecteur lié à la transformation

est appelé filtre.

25. Appliquer la fonction suivante à différents vecteurs

function Chaleur(S)

N=length(S);

for k=1:50

G=Gaussienne(N,k/5);

C=ConvFourier(S,G);

plot(S);hold on;plot(C,’r’);hold off;pause(0.1);

end

26. A quelle opération physique correspond la convolution par une gaussienne ?

27. A quelle opération correspond la convolution par le vecteur dec suivant :

>>dec=0*S;

>>dec(30)=1;

28. A quelle opération correspond la convolution par le vecteur d suivant :

>>d=0*S;

>>d(1)=1;d(2)=-1;

6 Transformée de Fourier 2D

Pour Matlab, une image de taille (m, n) est une matrice I de taille (m, n), et la valeur de I(i, j)

correspond à la valeur du niveau de gris de l’image au pixel (i, j).

Matlab est capable de lire à peu près tous les formats standards d’images.

6.1 Exemples : visualisation d’une image

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Chargement d’une image en Matlab:

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

load gatlin2;

% -> L’image est chargee dans la variable X

%Autres images:

%load clown; load gatlin; load mandrill;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Visualisation:

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

imagesc(X);

colormap gray;

%pour voir l’image en niveaux de gris

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Pour ouvrir une deuxieme figure:

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

figure(2);

colormap gray;

XX=imread(’cameraman.tif’);

imshow(XX);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

On peut aussi utiliser la commande imwrite, puis visualiser l’image sauvegardée avec son éditeur

d’image préféré (souvent plus pratique que ceux de Matlab).

Sources d’erreur classique : Il faut faire attention que les commandes imwrite ou imread suppose

que les images sont codées comme des entiers non signés (uint8 permet de faire la conversion), alors

4

que imagesc lit les images au format double, et surtout que vous ferez tous vos calculs au format

double. Si vous ne faîtes pas attention au format de vos données, vous obtiendrez de mangifiques

images noires en sortie de vos codes . . .

Choix d’une image : Vous pouvez continuer le TP avec votre image favorite (à condition qu’elle

soit de taille raisonnable). Sachant que le TP va être basé sur la transformée de Fourier, et qu’avec

matlab vous disposez de la FFT, donner une condition nécessaire pour que la taille soit raison-

nable ?

Si vous choisissez une image couleur, n’oubliez pas de la convertir en niveau de gris (e.g. a=mean(a,3)).

6.2 Un premier exemple de fonction Matlab : bruitage d’une image

function out = bruitage_gaussien(I,s)

%bruitage gaussien (d’ecart-type s) d’une image I

% bruitage_gaussien(I,s)

[m,n]=size(I);

J=zeros(m,n);

%creation du bruit gaussien

J=s*randn(m,n);

%bruitage

out=I+J;

7 Transformée de Fourier d’une image

En dimension 2, La TFD de l’image { f [k, l]} (k, l = 0, . . . , N) est définie par

ˆf [k, l] =

N

∑

n=1

N

∑

n=1

f[m, n]ω −(m−1)(k−1)

N

ω −(l−1)(n−1)

N

La transformée inverse est donnée par

ˆf [k, l] =

1

N2

N

∑

n=1

N

∑

n=1

ˆf [m, n]ω (m−1)(k−1)

N

ω (n−1)(l−1)

N

7.1 Visualisation du spectre d’une image synthétique

On va commencer par observer le spectre d’images synthétiques.

Sinusoïde

x=[0:255]’;

y=2pi(8/256)*ones(1,256);

a=sin(x*y);

29. Visualiser a dans une première figure. Pour visualiser son spectre, on peut faire :

A=fft2(a);

A=fftshift(A);

figure

imagesc(1+log(abs(A)));

colomap gray;

30. Expliquer le spectre observé (nombre de points lumineux, distance entre ces points,. . .).

Carré

b=zeros(256,256);

b(63:191,63:191)=1;

31. Visaliser b et son spectre comme précédemment. Commentaires ?

Sinusoïde fenêtrée Avec les notations précédentes :

c=b.sin(xy);

32. Visualiser c et son spectre comme précédemment. Commentaires ?

5

7.2 Visualisation du spectre d’une image réelle

33. Charger l’image clown. La mettre dans la variable a.

Calculer la TFD de l’image a (fonction f f t2). La stocker dans la variable A

Visualiser log(1 + abs(A)).

34. A quoi correspondent les structures horizontales, verticales, obliques, observées dans la TFD ?

35. Refaire l’opération avec d’autres images (gatlin, mandrill, ...). Commentaires ?

8 Phase et module de la transformée de Fourier

36. Nous allons illustrer l’importance respective de la phase et du module des coefficients de la Trans-

formée de Fourier. Nous commençons par échanger les phases des TFD de deux images.

– Choisisser une image b différente de a. Calculer les TFD A et B respectivement de a et b. Rem-

placer les coefficients de B par leurs modules, puis multiplier les par la phase des coefficients

correspondants de a (on pourra utiliser les fonctions abs et angle). Utiliser la fonction i f f t2, et

afficher le résultat.

– Commentaires ?

9 Filtrage passe-haut et images

Les contours dans une image constituent une information essentielle utilisée par notre système

visuel pour comprendre le contenu d’une image. L’extraction des contours est un point clé du

traitement d’images. Comme ces contours correspondent à des discontinuités, une méthode simple

pour les mettre en évidence consiste à appliquer les filtres suivants, passe-haut : Si I est une image,

>>f=0*I;f(1,1)=1;f(1,2)=-1;

>>g=f’;

Pour appliquer les filtres on effectue une convolution en passant par la transformée de Fourier 2D :

function C=ConvFourier2D(I,h)

C=real(ifft2(fft2(I).*fft2(h)));

Dans le cas présent le passage par Fourier n’est pas nécessaire car les filtres que nous utilisons sont

de tout petit support mais c’est toujours un moyen possible.

37. Construire et visualiser

>>g1=ConvFourier2D(I,f);

>>g2=ConvFourier2D(I,g);

puis construire une image par sommation des carrés des images g1 et g2 :

>>q=sqrt{g1^2+g2^2};

Il s’agit d’une approximation discrète du module du gradient des images. Appliquer cette procédure

à différentes images (éventuellement brutiée).

L’inconvénient d’une telle méthode est sa forte snesibilité au bruit et aux textures des images. On

procède donc en général à uen régularisation préalable de l’image. Effectuez l’opérération précé-

dente après avoir convolué l’image par un filtre gaussien 2D. Un tel filtre se construit comme le

produit de deux gaussiennes 1D.

6