[SCILAB] Convolution et transformée de Fourier

[exes74]

Bonjour,

je cherche actuellement à comparer le temps de calcul d'un produit de convolution en utilisant la convolution directe et en utilisant la convolution par transformées de Fourier. Le code que j'ai développé est ci-dessous, or il ne me donne pas la même courbe selon le type de convolution utilisée, je n'arrive pas à trouver mon erreur, quelqu'un peut m'aider ?

Code: Tout sélectionner


//Initialisation
xdel(winsid());     // Fermeture des fenetres
clear               // Suppression des variables
clc                 // effacement des lignes de la console
lines(0);

// définition des constantes
N=300 ; //nombre de points,
Te=0.1e-3 ; // période d'échantillonange
F=100 ; // fréquence du signal

t=Te*(0:N-1);// description du vecteur temps

f1=cos(2*%pi*F*t);//description du vecteur signal f1
f2=sin(2*%pi*F*t);//description du vecteur signal f2


////////////////////////////////////////////////////////////////////////////
//Calcul par convolution directe
timer()

y = convol(f1, f2);

timer();

//Calcul par transformée de fourrier

timer()

z = fft( abs(fft(f1)) .* abs(fft(f2)) , 1);

timer();



// affichage sur la fenêtre 0
  xset('window',0); // sélection de la fenêtre 0
  xbasc(0); //effacer la fenetre 0
  xset("font size",4); //définition de la taille de la police
  xsetech([0, 0, 0.5,0.5]); // le haut gauche de la fenetre coupée en 4
 
  plot(f1) // affichage de s en fonction du rang
  xtitle('Sinusoide 1');

// affichage du module brut de forme
  xsetech([0,0.5,0.5,0.5]) ;// le bas gauche de la fenetre coupée en 4
  plot(f2) ;
  xtitle('Sinusoide 2');

// affichage du module corrigé de 0 à fe/2
  xsetech([0.5,0,0.5,0.5]) ; //le haut droit de la fenetre coupée en 4
  plot(y);
  xtitle('Convolution directe');

// zoom sur la raie, de la fréquence nulle à fe/5
  xsetech([0.5,0.5,0.5,0.5]) ; //le bas droit de la fenetre coupée en 4
  plot(z);
  xtitle('Convolution par transformées de Fourier');



Merci par avance !

[busard_des_roseaux]

bonsoir,

c'est normal les valeurs absolues abs dans la seconde expression ?

[exes74]

Bonsoir,

c'était un test, mais le résultat est le même avec ou sans l'utilisation des valeurs absolues.

[kazeriahm]

attends ce qu'il fait comparer c'est f*g (convolution) et F;)1(F(f)*F(g)) (*-> multiplication) non ?
toi tu compares le premier et F(f fois g)...

d'ailleurs (sous Matlab en tout cas) quand tu lui demandes plot(z) avec z un tableau de complexes, il te trace Im(z) en fonction de Re(z)

[exes74]

convol(f,g) => convolution directe



F;)1(F(f) .* F(g)) => c'est une convolution aussi, en passant par une double transformée de Fourier (sous Scilab .* => mutiplication membre à membre)

donc je compare bien deux convolutions. Non ?

C'est peut être plus clair comme ça:

Code: Tout sélectionner

//Initialisation
xdel(winsid());     // Fermeture des fenetres
clear               // Suppression des variables
clc                 // effacement des lignes de la console
lines(0);

// définition des constantes
N=300 ; //nombre de points,
Te= 1 ; //0.1e-3 ; // période d'échantillonange
F=100 ; // fréquence du signal

t=Te*(0:N-1);// description du vecteur temps

f1=2*cos(2*%pi*F*t);//description du vecteur signal f1
f2=sin(2*%pi*F*t);//description du vecteur signal f2

////////////////////////////////////////////////////////////////////////////
//Calcul par convolution directe
timer()

y = convol(f1,f2);

timer();

//Calcul par transformée de fourrier

timer()

F1=fft(f1,-1);
F2=fft(f2,-1);

z = fft( (F1 .* F2), 1);

timer();



// affichage sur la fenêtre 0
  xset('window',0); // sélection de la fenêtre 0
  xbasc(0); //effacer la fenetre 0
  xset("font size",4); //définition de la taille de la police
 
  xsetech([0, 0, 0.3,0.3]); // le haut gauche de la fenetre coupée en 4
  plot(f1)
  xtitle('Sinusoide 1');

  xsetech([0.3,0,0.3,0.3]) ;// le bas gauche de la fenetre coupée en 4
  plot(F1) ;
  xtitle('TF de la sinusoide 1');
 
  xsetech([0, 0.3, 0.3,0.3]); // le haut gauche de la fenetre coupée en 4
  plot(f2)
  xtitle('Sinusoide 2');

  xsetech([0.3,0.3,0.3,0.3]) ;// le bas gauche de la fenetre coupée en 4
  plot(F2) ;
  xtitle('TF de la sinusoide 2');

  xsetech([0.0,0.7,0.3,0.3]) ; //le haut droit de la fenetre coupée en 4
  plot(y);
  xtitle('Convolution directe');

  xsetech([0.3,0.7,0.3,0.3]) ; //le bas droit de la fenetre coupée en 4
  plot(z);
  xtitle('Convolution par transformées de Fourier');





[kazeriahm]

au temps pour moi, j'ai mal lu ton code

il faut rajouter des 0 artificiels aux deux vecteurs que tu veux convoler comme expliqué dans la page "convol" de l'aide Matlab :

Algorithm

The convolution theorem says, roughly, that convolving two sequences is the same as multiplying their Fourier transforms. In order to make this precise, it is necessary to pad the two vectors with zeros and ignore roundoff error. Thus, if

X = fft([x zeros(1,length(y)-1)])

and

Y = fft([y zeros(1,length(x)-1)])

then conv(x,y) = ifft(X.*Y)
See Also

[exes74]

Bonjour,

effectivement, en remplissant avec des zéros ça marche bcp mieux.

Une autre solution aurait été d'utiliser des signaux bornés.

Merci beaucoup !