Discussion :

[diogene]

Si on part du principe qu'il faut d'abord remplir le tableau dataout vide avec les valeurs datain avant d'utiliser SWAP, ce n'est vraimment pas interessant en nombres de calcul.

Non : on fait le swap dans la foulée.
Je me suis sans doute mal exprimé : on utilise la même méthode que pour le swap (c.a.d. passer par un stockage intermédiaire ).
ça donne

Code :

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

1
2
3
4
5
6
7
8
9
10
   if (j >= i) 				
   {  
      double prov ;		
       prov = datain[j-1];
       dataout[j-1]=datain[i-1];
       dataout[i-1]=prov;
       prov = datain[j];
       dataout[j]=datain[i];
       dataout[i]=prov ;
   }

Et ça marche que les tableaux soient les mêmes ou non. Le surcoût de calcul sur la fft est absolument négligeable.

[nakor]

Bonjour,
une petite question au niveau du remplissage du tableau DATA[2*N]

Il est conseille de remplir le tableau d'une maniere particuliere pour les fonctions paires : la partie positive dans les premieres valeurs de data (data[0], data[1]...) et la partie negative a la fin du tableau (data[n],data[n-1]...)

Quel est l'interet d'une telle repartition ? est ce que l'on va gagner en temps de computation ?

Si on suit cette repartition des points, une fonction Lorentziennehttp://fr.wikipedia.org/wiki/Fonction_lorentzienne
serait ecrite de cette maniere dans mon tableau ?

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
 
 
 
double *func1 (unsigned long nn, double r) //r=largeur a mi-hauteur
{                                                        //nn : nombre de point complexe
 
   double j0=0 ;         //j0 : abscisse du sommet
   double *ptr,*ub,*fp;
   unsigned long nncp = nn<<1;
   unsigned long nah =(unsigned long)((nn)/2);
 
   fp=(double*)(malloc(nncp*sizeof(double))); 
      ub=fp+nncp;
 
   for (ptr=fp;ptr<ub;ptr++)*ptr=0.0
 
   for (ptr=fp;ptr<ub-nah;ptr+=2) 
   {
      *ptr=0,5*r/(PI*(i-jo)*(i-jo)+(0,5*r)*(0,5*r));
      i+=2; 
   } 
 
   for (ptr=ub;ptr>=ub-nah;ptr-=2) 
   {
      *ptr=0,5*r/(PI*(-i-jo)*(-i-jo)+(0,5*r)*(0,5*r));
      i+=2; 
   } 
 
   return fp;      
}

[diogene]

On ne gagne rien en temps de calcul.
La question est de savoir quoi faire des données pour t<0 si il y en a (comme dans une fonction paire).
Il faut savoir que la transformée de Fourier discrète (TFD) calculée par la FFT, est en fait la transformée de Fourier de la fonction rendue périodique (parce que la TFD est elle-même échantillonnée) par décalage et superposition.

Un petit dessin :

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
Autrement dit, si j'ai cette fonction échantillonnée sur N points :
         ****
        *    *
       *      *
...****   0   ****...

Sa TFD est la TF de 
        ****            ****            ****
       *    *          *    *          *    *
      *      *        *      *        *      *
...***        ********   0    ********    N   ****....

Une période de cette fonction est donc de 0 à N :
          **            **
            *          *       
             *        *         
          0   ********    N-1
La partie négative est reproduite à la fin de la partie positive.

Si j'utilise comme donnée : 
                ****
               *    *
              *      *
          ****        ****
          0              N-1

Les données sont décalées de N/2 points.
La conséquence est un gradient de phase sur la TFD.
Comme un décalage de 1 point entraine un gradient de phase de 2pi/N
par point sur la TFD, cela entraine un gradient de phase de pi par point
sur la TFD.

[nakor]

Merci, c'est beaucoup plus facile de comprendre a l'aide du schema.

Sinon, un autre point a eclaircir : il semblerait que lorsque cet algorythme a ete ecrit, il etait tres longt de calculer cos(n*x) ou sin(n*x), on utilisait donc des relations trigonometrique pour arriver au meme resultat.

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
 
while (n>mmax) 
{
   istep=mmax <<1;
   theta=isign*(6.28318530717959/mmax);
   wtemp=sin(0.5*theta);
   wpr=-2.0*wtemp*wtemp;
   wpi=sin(theta);
   wr=1.0;
   wi=0.0;
   for (m=1;m<mmax;m+=2) 
   {
      for (i=m; i<=n;i+=istep) 
      {  [...]
      }
      wr=(wtemp=wr)*wpr-wi*wpi+wr;
      wi=wi*wpr+wtemp*wpi+wi;
   }
   mmax=istep;
}

Mais maintenant, il semble que ce n'est plus tres utile, ni interessant d'utiliser ces relations.

On peut donc modifier le code comme suit ?

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
 
while (n>mmax) 
{
   istep=mmax <<1;
   theta=isign*(6.28318530717959/mmax);
   //wtemp=sin(0.5*theta);
   //wpr=-2.0*wtemp*wtemp;
   //wpi=sin(theta);
   wr=1.0;
   wi=0.0;
   for (m=1;m<mmax;m+=2) 
   {
      for (i=m; i<=n;i+=istep) 
      {  [...]
      }
      //wr=(wtemp=wr)*wpr-wi*wpi+wr;
      wr=cos(m*theta);
     //wi=wi*wpr+wtemp*wpi+wi;
      wi=sin(m*theta);
   }
   mmax=istep;
}

[diogene]

Il ne me semble pas que cela soit souhaitable :
Ici, les sin/cos sont calculés avec au total 8 opérations float (* et +-). Il est certain pour moi que les fonctions sin/cos de la bibliothèque (sans compter le temps d'appel et de retour des fonctions) doivent demander plus de 4 opérations float chacune.
Je pense donc qu'on perdra du temps et qu'on n'en tirera aucun avantage.