FFT avec Matlab - fonction fftshift

[fabio123]

Bonjour,

j'essaie de comprendre un résultat que j'obtiens avec Matlab sur des FFT et plus particulièrement avec la fonction fftshift

Je calcule d'abord la FFT d'un signal simple ( y=cos(2*pi*f_signal*t) puis la transformée inverse pour essayer de retrouver le signal initial.

Je compare les résultats sans et avec du "zero padding" en rajoutant 100 points de plus au signal de départ avant d'effectuer la FFT.

voici le code :

% Frequency input signal and time vector
f_signal=10;
nb_points = 100;
f_sampling=nb_points*f_signal;
step_time=1/f_sampling;
t=linspace(0,1/f_signal,nb_points);


%% Original signal with 100 points
y=cos(2*pi*f_signal*t);

%% Adding 100 more points for y
n_zero_padding=100;
y_f=padarray(y,[0 n_zero_padding],'post');

% Plot Input Signal
figure(1);
plot(y);
title('Input Signal');

% Perform Forward FFT
z=real(fft(y_f));

% Center FFT
z_centered=fftshift(z);
figure(2);
plot(z_centered);
title('FFT centrered with zero padding');

% Without centering
figure(3);
plot(z);
title('FFT with zero padding');

% Perform Inverse FFT
z_inv=ifft(z);  
figure(4);
plot(z_inv);
title('Inverse FFT with zero padding');

% Inverse FFT with centered spectrum
z_inv_centered=ifft(z_centered);  
figure(5);
plot(z_inv_centered);
title('Inverse FFT centered with zero padding');

Je ne comprends pas le résultat obtenu sur la figure 5 (transformée inverse d'un spectre décalé avec fftshift);

voici le premier spectre obtenu (non décalé) (figure(3)) :

Spectre_not_shifted.png

voici le premier spectre décalé avec fftshift (figure(2)) :

Spectre_shifted.png

et la TF inverse du spectre décalé ci-dessus (figure(5)) :

Inverse_FFT_of_Shifted_Spectre.png

Je ne comprends pas cette dernière figure. Si je reprends la théorie, le fait de décaler le spectre avec fftshift conduit à l'obtention de 2 Dirac centré autour de f_sampling/2, plus précisément en f_sampling/2-f_signal et f_sampling/2+f_signal.

J'ai donc la somme de ces 2 Dirac et je fais la TF Inverse de cette somme : je devrais donc obtenir pour la TF inverse la somme de 2 cosinus avec des fréquences très proches : on dirait que j'obtiens sur la figure(5) le produit de 2 signaux cosinus mais je ne vois pas où est mon erreur d'interprétation ...

Si quelqu'un pouvait m'expliquer ce résultat ?

Merci par avance
-----

[Neluge]

Bonjour,

Je ne vois pas d'erreur d'interprétation. En vertu des formules de trigo, une somme de cosinus, c'est pareil qu'un produit de cosinus. Ce que vous observez correspond a :

[fabio123]

merci, j'avais oublié cette formule.

Juste une dernière chose, J'ai fait un plot de la somme de 2 cosinus pour essayer de retrouver la dernière figure ci-dessus :

t=linspace(0,2/f_signal-2/(f_signal*2*nb_points),2*nb_po ints);
y1=cos(2*pi*(f_sampling/2-f_signal)*t);
y2=cos(2*pi*(f_sampling/2+f_signal)*t);
sum_y=y1+y1;
figure(1);
plot(sum_y);

J'obtiens la figure suivante :



Comme vous pouvez voir, il y a un facteur 4 en amplitude par rapport à la figure ci-dessus.

Pourriez-vous me dire d'où vient ce facteur puisque :

TF(cos(2*pi*f_s)) = 1/2 delta(f-f_s) + 1/2 delta(f+f_s)

Merci par avance pour votre aide

[Neluge]

Bonjour,

D'abord je reviens sur mon post précédent, votre interprétation ne me semble, finalement, plus très logique... En effet, pour effectuer la transformée de Fourier inverse, Matlab considère que la première moitié du tableau correspond aux fréquences positives (de 0 à f_sampling/2) et la seconde moitié aux fréquences négatives (à partir de -f_sampling/2). fftshift ne fait qu'intervertir les deux moitiés du tableau, pour pouvoir visualiser un spectre centré en 0, mais quand vous donnez à la fonction ifft le tableau représenté sur la figure 2 (spectre décalé), le premier pic sera associé à la fréquence positive , et le second à une fréquence négative, soit . La transformée de Fourier inverse correspondante devrait être... une simple sinusoïde de fréquence ... non ?

Sinon, une idée peut-être pour cette diminution d'amplitude : quand vous faites ifft(fft(y_f)) (càd transformée de Fourier puis transformée inverse), et qu'on trace le résultat (plot), on voit qu'on retombe bien sur le signal y_f de départ. Logique. Mais quand on fait ifft(real(fft(y_f))), le signal sur lequel on retombe admet deux différences par rapport au signal de départ :
-> la sinusoïde s'étale sur les 200 points (et non plus juste sur les 100 premiers) ;
-> l'amplitude a diminué d'un facteur 2.
Il n'est pas anormal qu'on ne retombe pas sur le signal de départ : la transformée de Fourier d'un cos fenêtré est rigoureusement constituée de deux "diracs" (en fait des sinus cardinaux) s'il est centré en 0, mais ici, il est défini entre 0 et T, donc le résultat est un nombre complexe (sinus cardinaux multipliés par un facteur de phase). En ne retenant que la partie réelle, on perd des informations de phase, donc il est compréhensible qu'en faisant la TF inverse, on obtienne un signal différent du signal du départ. Mais cela pourrait peut-être aider si l'on pouvait s'expliquer pourquoi la sinusoïde n'est plus tronquée et son amplitude est divisée par 2 ?

Bref, je ne vous apporte peut-être pas beaucoup de réponses, mais ce sont peut-être des pistes de réflexion...

[A voir en vidéo sur Futura]

[Neluge]

Bonjour,

En réfléchissant j'ai trouvé une explication aux problèmes soulevés dans le second paragraphe : si l'on considère un signal sinusoidal de pulsation omega0 défini entre -T et T (donc un cos multiplié par une fonction "porte" entre -T et T), la théorie dit que sa transformée de Fourier sera la somme de deux demi-diracs (TF du cos) convolée a un sinus cardinal (TF du carré). Mais ici, le signal est défini a partir de t=0 : il est décalé de T. La théorie dit donc que sa TF est multipliée par un facteur exp(-i omega T) :

(l'étoile représente le produit de convolution). La, on peut vérifier qu'en faisant la TF inverse, on retrouve toujours le signal de départ.
Mais que se passe-t-il quand vous prenez la partie réelle ? ("real(fft(y_f))")
devient . En faisant la TF inverse (la TF inverse d'un cos étant constituée de deux diracs en temps... c'est symétrique), on obtient donc :

ou j'ai noté H la fonction "porte". La convolution par un dirac a pour effet de translater le signal (de -T ou de +T ici) : ce qu'on observe donc en faisant ifft(y_f), c'est un demi-cosinus fenêtré centré en -T (défini entre -2T et 0) et un demi-cosinus centré en T (défini entre 0 et 2T).
Ce qui explique déja pourquoi l'amplitude est divisée d'un facteur 2. Enfin, quand vous tracez le résultat ("plot"), Matlab le représente entre 0 et 4T (on a considéré, au départ, un signal sinusoidal fenêtré entre 0 et 2T, et on l'a complété par des zéros entre 2T et 4T - "zero-padding"). Entre 0 et 2T, on observe donc le premier demi-cosinus ; entre 2T et 4T, on observe l'autre demi-cosinus (de la même façon qu'en faisant la FFT, la seconde moitié du tableau correspond en fait aux fréquences négatives).

Voila, désolé de répondre a des questions que vous n'avez pas posées, mais ceci explique déja pourquoi l'amplitude est au moins diminuée d'un facteur 2 (a défaut de comprendre ce facteur 4... et pourquoi on n'a pas juste un cosinus) : c'est dû au fait que vous prenez la partie réelle de la FFT avant de faire la transformée inverse.

[Krati]

Bonjour FABIO123 !

J'aurais besoin de ton aide. En effet, en 2012, tu a posté un message sur l'écoulement de l'air et tu as modélisé cet écoulement à l'aide d'un logiciel sympa dont je cherche le nom pour pouvoir l'utiliser!

Merci de bien vouloir me répondre!

PS : voici le lien de la discussion sur l'écoulement de l'air : http://forums.futura-sciences.com/ph...ml#post5193912