FFT et FFT Inverse sous Matlab

[fabio123]

Bonsoir,

je ne sais pas si je suis sur le bon forum comme ça concerne les transformées de Fourier ...

je rencontre un problème sous Matlab à propos des FFT. 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 précise que pour avoir une meilleure précision, je fais du "zero padding" en rajoutant 100 points de plus au signal de départ avant d'effectuer la première 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');

%% No zero padding
%y_f=y;

% 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');

J'essaie de voir l'effet produit par la fonction "fftshift" (qui centre le spectre obtenu autour de la fréquence 0) avant d'effectuer la transformée inverse.

Je m'aperçois malheureusement qu'en utilisant cette fonction, je ne retrouve pas le signal original.

Voici la figure représentant le spectre sans utiliser "fftshift" :

*** Image sur un hébergeur externe : merci de respecter les règles de ce forum ***

On distingue bien les 2 diracs convolués en 0 et f_sampling.


ici la figure représentant le spectre avec "fftshif" :

*** Image sur un hébergeur externe : merci de respecter les règles de ce forum ***


Maintenant la FFT inverse du spectre obtenu sans "fftshift" :

*** Image sur un hébergeur externe : merci de respecter les règles de ce forum ***

On retrouve bien le signal de départ.

Maintenant, la FFT inverse du spectre obtenu avec "fftshift" :

*** Image sur un hébergeur externe : merci de respecter les règles de ce forum ***


Je ne comprends pas pourquoi cette FFT inverse ne donne pas le signal de départ ?

Le spectre obtenu avec fftshift représente bien pourtant la Transformée de Fourier d'un signal cosinusoïdal avec les 2 dirac situés en f = - f-signal et f = f_signal .

Toute aide est la bienvenue
-----

[acx01b]

% Frequency input signal and time vector
f_signal=10; % 10 hertz
duration = 1 % 1 seconde
f_sampling = 100; % 100 hertz

nb_points = duration * f_sampling; 

step_time=1/f_sampling;
t=linspace(0,duration - step_time,nb_points); %  [0,0.01,0.02...,0.99]

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

%% Adding 100 more points for y
n_zero_padding=100;
y_f=[y,zeros(1,n_zero_padding)];

%% No zero padding
%y_f=y;  

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

% Perform Forward FFT
z=fft(y_f);


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

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

% Perform Inverse FFT
z_inv = real(ifft(z));
figure(4);

% on pourrait aussi plot la phase
subplot(2,1,1);
plot(abs(z_inv)); 
subplot(2,1,2);
plot(angle(z_inv)); 




title('Inverse FFT with zero padding');

% Inverse FFT with centered spectrum
z_inv_centered=real(ifft(z_centered));  
figure(5);
%plot(abs(z_inv_centered));

subplot(2,1,1);
plot(abs(z_inv_centered)); 
subplot(2,1,2);
plot(angle(z_inv_centered)); 

title('Inverse FFT centered with zero padding');

[acx01b]

si y est un vecteur de réels

z = fft(y)

z est un vecteur de complexe

y2 = ifft(fft(y))

y2 est également un vecteur de complexe

et real(y2) est égal à y aux erreurs d'arrondi près



j'ai donc oublié de mettre des abs :

plot(abs(z))
...
plot(abs(z_centered))

[fabio123]

Bonsoir,

la FFT d'une fonction paire et périodique est purement réelle, c'est pour cela que j'arrondis le calcul de la première FFT en prenant :

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

Est-ce cette différence ( entre le calcul de

z=fft(y_f)

et de

z=real(fft(y_f))

) qui fait que je ne retrouve pas, après la fft inverse du spectre obtenu avec fftshif, le signal initial ?

c'est-à-dire avoir la figure suivante :

*** Image sur un hébergeur externe : merci de respecter les règles de ce forum ***

au lieu de ce signal d'entrée :

*** Image sur un hébergeur externe : merci de respecter les règles de ce forum ***

mais sur 200 points

Est-ce un problème de résolution spectrale ou de perte d'informations, ou à cause simplement du décalage autour de f=0 avec fftshift ?

Merci pour vos éclaircissements

[A voir en vidéo sur Futura]

[gg0]

A vérifier : échantillonner une fonction paire ne donne ps obligatoirement un résultat pair.

Cordialement.

[fabio123]

Bonjour,

je confirme que la FFT d'un signal réel, de symétrie paire et périodique est purement réelle.

Je sais que le problème vient de l'utilisation de "fftshift" mais impossible de savoir pourquoi ...

Merci pour vos éclaircissements

[gg0]

Je confirme qu'il ne faut pas confondre calcul exact avec calcul approché

[fabio123]

oui, finalement le "zero padding" fait apparaitre des termes imaginaires non négligeables ...

voici la figure obtenue après une FFT et une FFT inverse sans passer par fftshift :

*** Image sur un hébergeur externe : merci de respecter les règles de ce forum ***

Je retrouve le signal de départ avec les zéros (100 points en +)

reste à résoudre le problème du décalage avec fftshift pour la transformée inverse

merci pour votre aide

[fabio123]

désolé, voici l'image en question :

[gg0]

Je ne sais pas trop ce que tu as fait, mais pour avoir eu des ennuis et avoir été obligé de rajouter fftshift, je sais que cette fonction sert à traiter une FFT de signal pair (construit sur un intervalle centré en 0 - dans mon cas, il s'agissait d'une gaussienne). Si tu es parti d'un échantillon pris seulement pour les valeurs positives, il n'y a pas besoin du fftshift. Pour voir son intérêt, il va falloir partir de ton cos échantillonné de -a à a, pas de 0 à b.

Cordialement.

NB : je ne rentre pas dans les détail, ça fait des années que je n'utilise plus scilab.

[acx01b]

il n'est pas bien mon code ? tu l'as essayé au moins ?

les deux points importants :
- les constantes que je définis sont à mon avis beaucoup plus claires
- le abs et angle de fft(y) si c'est réel alors angle vaut 0 ou pi
- le real(ifft(..)) qui ne sont pas présents dans ton code
- je ne sais pas si tu sais mais ifft(z) c'est pareil que conj(fft(conj(z)))/length(z)

[fabio123]

Je te remercie pour ton code, j'y vois un peu + clair.

je ne comprenais pas la figure que j'obtenais avec la fft inverse de la fft que j'ai décalée avec fftshift.

En fait, fftshift décale les 2 diracs et les centre autour du point situé à "n_points_fft/2" (et non autour de la fréquence 0), et donc en faisant la fft inverse,
j'obtiens un signal modulé (un produit de 2 cosinus) avec ces 2 fréquences très proches l'une de l'autre.

Merci pour votre aide

[invite3abecd25]

bonjour mes chères amis ,j'espère que tout le monde va bien .
Alors mes amis au aujourd’hui j'ai trouver une autre problème le voila je cherche a tracer la fonction xp en fonction de x
xp = (1/Z)*A
avec:
A=<gq,E0> avec p varie entre 1 et 100
Z=<gq,fm0> zmn <fm0,gq> avec p varie entre 1 et 50 et m varie entre 1 et 300

les expressions:
E0=-sqrt(2/a*b)*sin(pi*x/a)
gp=-sin(p*pi*x/a1) sur intervalle [0,a1] et gp=-sin(p*pi(x-a2)/(a-a2)) sur intervalle [a2,a]

zmn=1i*nu0*omega/(sqrt((m*pi/a)^2-k0^2))

les constants de mon problème
a=22.9e-3; b=10.16e-3 ;w=a/9; a1=6.8e-3; a2=a1+w; a3=a-(a1+a2) nu0=4*pi*1e-7 f=1e9; e0=1/(36*pi)*1e-9; c=1/sqrt(e0*nu0); omega=2*pi*f;
e k0=omega/c

j'ai crée cette code mais ne donne pas de résultats , s'ils vos plait les amis donne moi de l'aide, soit me corrige le code ou bien me une autre code simple

je compte sur vos mes chères amis
voila le code

clear all
clc;
a=22.9e-3;
b=10.16e-3 ;
w=a/9;
a1=6.8e-3;
a2=a1+w;
a3=a-(a1+a2)
nu0=4*pi*1e-7;
f=1e9;
e0=1/(36*pi)*1e-9;
c=1/sqrt(e0*nu0);
omega=2*pi*f;
k0=omega/c
P=10;
M=2;

xx=[0:0.1e-3:22.9e-3];
for s=1:length(xx)
x=xx(s);

for m=2:M

for p=1
A(p,1)=sqrt(2/a*b)*quad(@(x)sin(pi*x/a).*sin(p*pi*x/a1),0,a1)+sqrt(2/a*b)*quad(@(x)sin(pi*x/a).*sin(p*pi*(x-a2)/(a-a2)),a2,a);
Fm0_gp(p,1)= sqrt(2/a*b)*quad(@(x)sin(m*pi*x/a).*sin(p*pi*x/a1),0,a1)+sqrt(2/a*b)*quad(@(x)sin(m*pi*x/a).*sin(p*pi*(x-a2)/(a-a2)),a2,a);


end

ZTE=1i*nu0*omega/(sqrt((m*pi/a)^2-k0^2));
B= Fm0_gp*ZTE* Fm0_gp';
xp(1,s)=inv(B)*A;
X(1,s)=x;
s
end

end
plot(X,abs(xp),'.');
grid on
ylabel('xp');
xlabel('X');


et merci d'avance