Discussion :

[nakor]

Bonjour,
(mon précédent message s'étant noyer dans les abimes du forum algorythme j'essaye ici, si quelqu'un s'est déjà intéresser a la convolution de NRiC, il pourra sans doute m'aider)

je cherche a réaliser la convolution d'une série de données (data) avec une fonction réponse.

Mais j'ai un peu de mal a comprendre l'algorithme de NRiC. Pourtant sur le principe, ça devrait être plutôt simple, en utilisant la FFT, la convolution devrait être le produit en chaque point de la transformée de fourier de data et de réponse non ? en tenant compte du facteur d'echelle, 1/N.

Pourtant dans NRiC, ce n'est pas ce qui est effectue, on multiplie
FFT(data[i])*FFT(reponse[i]) puis on y soustrait FFT(data[i+1])*FFT(reponse[i+1])
(FFT(X) etant la FFT de data, on ne fait pas appelle a la FFT a chaque opération)

Et pour la deconvolution (-1 dans le code) c'est encore plus étrange, la ou devrait trouver une division pour revenir a data, c'est toujours une multiplication.

Si vous pouvez m'éclairer un peu, ça m'avancerait bien!
merci

dans le code (NRiC, $12.3), ans et la transformée de fourrier de Data, temp et la transformée de fourrier de réponse,

ans (et temp) sont stockées tel que
" in array ans[0..n-1]) by the positive frequency half of their complex Fourier transform. The real-valued first and last components of the complex transform are returned as elements ans[0] and ans[1], respectively "

no2 vaut n/2

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
 
if (isign == 1) 
{
	for (i=2;i<n;i+=2) 
	{ //Multiply FFTs to convolve.
		tmp=ans[i];
		ans[i]=(ans[i]*temp[i]-ans[i+1]*temp[i+1])/no2;
		ans[i+1]=(ans[i+1]*temp[i]+tmp*temp[i+1])/no2;
//j'aurai plutôt écrit directement 
//ans[i]=ans[i]*temp[i]/n;
	}
	ans[0]=ans[0]*temp[0]/no2;
	ans[1]=ans[1]*temp[1]/no2;
} 
else if (isign == -1) 
{
	for (i=2;i<n;i+=2) 
	{// Divide FFTs to deconvolve.
		if ((mag2=SQR(temp[i])+SQR(temp[i+1])) == 0.0)
		throw("Deconvolving at response zero in convlv");
		tmp=ans[i];
		ans[i]=(ans[i]*temp[i]+ans[i+1]*temp[i+1])/mag2/no2;
		ans[i+1]=(ans[i+1]*temp[i]-tmp*temp[i+1])/mag2/no2;
//j'aurai plutôt écrit directement 
//ans[i]=ans[i]/tmp[i]/n;
 
	}
	if (temp[0] == 0.0 || temp[1] == 0.0)
	throw("Deconvolving at response zero in convlv");
	ans[0]=ans[0]/temp[0]/no2;
	ans[1]=ans[1]/temp[1]/no2;
} 
else throw("No meaning for isign in convlv");
 
realft(ans,-1); Inverse transform back to time domain.

[diogene]

Pourtant dans NRiC, ce n'est pas ce qui est effectue, on multiplie
FFT(data[i])*FFT(reponse[i]) puis on y soustrait FFT(data[i+1])*FFT(reponse[i+1])

Je ne comprend pas ce que cette expression veut dire. Est-ce

FFT(data)[i]*FFT(reponse)[i] puis on y soustrait FFT(data)[i+1]*FFT(reponse)[i+1] ?

Fais-tu allusion à ce code ?

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

ans[i]=(ans[i]*temp[i]-ans[i+1]*temp[i+1])/no2

[nakor]

C'était pour dire que l'on prend la transformée de Fourier de data (dans le code ans) en un point, et qu'on la multiplie a la transformée de Fourier de réponse au même point.
Oui c'est la même chose que ce morceau de code mais je pensai que ans et temp ne parlaient pas trop si on ne sait pas qu'ils sont la transformée de Fourier de data et réponse.
Bien sur ça ne veut pas dire qu'on appelle la FFT à chaque points, c'est plus une notation (maladroite !) pour expliquer qu'on utilise les transformées de Fourier de data et de réponse.

[diogene]

devrait être le produit en chaque point de la transformée de fourier de data et de réponse non

C'est bien ce qui est écrit dans le code : attention, les FFT donnent comme résultat des complexes et il s'agit là de la multiplication complexe des points:

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
tmp=ans[i];
ans[i]=(ans[i]*temp[i]-ans[i+1]*temp[i+1])/no2;
ans[i+1]=(ans[i+1]*temp[i]+tmp*temp[i+1])/no2;

De même,

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
tmp=ans[i];
ans[i]=(ans[i]*temp[i]+ans[i+1]*temp[i+1])/mag2/no2;
ans[i+1]=(ans[i+1]*temp[i]-tmp*temp[i+1])/mag2/no2;

représente bien la division : (a+ib)/(c+id) = (a+ib)(c-id)/(c*c+d*d)

[nakor]

J'avais compris qu'avec l'algorithme reelfft, la partie imaginaire était uniquement stockés dans data[0] et data[1]
"The real-valued first and last components of the complex transform are returned as elements ans[0] and ans[1], respectively"

euh j'ai vraiment pas prêter attention a ce qui était écrit. désolé :S

Sinon en admettant que les donnes sont stockes reel,im,reel,im, c'est vrai que c'est logique de faire cette opération.

Mais en acceptant ça, je ne trouve pas les résultats que je pense devoir obtenir.
la convolution de deux fonctions réelles devrait donner une fonction réelle non ? et ce n'est pas ce que j'obtiens, avec cette algorithme je trouve une partie imaginaire

[diogene]

la convolution de deux fonctions réelles devrait donner une fonction réelle non ?

Oui

et ce n'est pas ce que j'obtiens, avec cette algorithme je trouve une partie imaginaire

Où trouves-tu les valeurs que tu qualifies d'imaginaires ?
Apparemment, la FFT utilisée est dédiée aux suites réelles :
La FFT directe ne calcule alors que la moitié des points (les autres étant complexes conjugués des premiers) La FFT inverse part de cette supposition sur les données et ne donne que des valeurs réelles. Les valeurs que tu penses imaginaires sont des valeurs réelles.

[Complément]

Cette méthode, passant par les transformées de Fourier discrète, calcule la convolution circulaire de deux suites de P points et donne une suite de P points.
Pour le genre d'application que je pense que tu as, il te faut une convolution linéaire.
La convolution linéaire d'une suite de N points et d'une suite de M points donne une suite de N+M-1 points. Si tu veux utiliser un algo de convolution circulaire pour calculer la convolution linéaire, il faut étendre tes suites à P points, P>= N+M-1 compatible avec ton algo, en les complétant par des zéros. Tu obtiendras la convolution linéaire sur P points.