convolution, a partir de la FFT
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:
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. |
