convolution, NRiC

convolution, NRiC - C - Programmation

Marsh Posté le 10-12-2009 à 18:35:47    

Bonjour,
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])  
 
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, ca m'avancerait bien!
merci
 
dans le code, ans et la transformée de fourrier de Data, temp et la transformée de fourrier de réponse,
ans (et temp) est stockée 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 "
 

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


 

Reply

Marsh Posté le 10-12-2009 à 18:35:47   

Reply

Marsh Posté le 10-12-2009 à 21:06:06    

Tes deux vecteurs temp et ans sont à valeurs complexes, stockés partie réelle, partie imaginaire, partie réelle, partie imaginaire etc...
Ca fait (a+ib)*(c+id), redispatché en partie réelle, partie imaginaire.
 
Whé je sais encore le faire  [:zcoold]

Reply

Marsh Posté le 11-12-2009 à 12:39:13    

autant pour moi je n'avais pas fait attention au fait que ca pouvait etre une multiplication de complexe.
 
En effet j'étais persuadé que l'algo reelfft renvoyer des reels, et que les seules complexes étaient stockes dans data[0] data[1],
j'ai lu trop vite la description de l'algorithme et l'anglais m'a joué un vilain tour.
 
Merci pour ton aide, c'était pas grand chose mais j'étais tellement persuadé de travaillé en réelle que je ne comprenais pas cette opération  :sweat:

Reply

Sujets relatifs:

Leave a Replay

Make sure you enter the(*)required information where indicate.HTML code is not allowed