code numerical recipe

Version imprimable

09/02/2007, 21h59
deubelte

code numerical recipe

Si vous connaissez le site Numerical recipe, que pensez vous de leurs codes C++?
en terme de rigeur, rapidité...?
merci
09/02/2007, 22h10
roulious

La qualité (en terme code pur, pas de performance) est à mon avis déplorable. C'est principalement du C sur lequel ils ont jeté deux ou trois fonctionalités du C++. Les variables sont toutes déclarées en début de fonction, ils mettent du using namespace std partout et ne font aucune utilisation de la bibliothèque standard.

Le code est par contre conforme à l'algorithme décrit dans le texte, que ce soit au niveau de la correction ou de la complexité.

La base du code est là, mais si tu veux le modifier un tant soit peu, mieux vaut commencer par le ré-écrire un peu plus proprement.

en fait, voila pourquoi je pose cette question. LE code ci dessous est tiré de NR.
Il permet de calculer la FFT d'un 2^n points .

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
36
37
38
39
40
41
42
43
44
45
46
47
 
void four1(float data[],unsigned long nn,int isign)
{
    unsigned long n,mmax,m,j,istep,i;
    double wtemp,wr,wpr,wpi,wi,theta;
    float tempr,tempi;
 
  n=nn<<1;
   j=1;
   for (i=1;i<n;i+=2){
         if(j>i){
             SWAP(data[j],data[i]);
             SWAP(data[j+1],data[i+1]);
         }
         m=nn;
         while (m>=2 && j>m){
             j-=m;
             m>>=1;
         }
         j+=m;
   }
   mmax=2;
   while(n>mmax){
       istep=mmax<<1;cout<<istep;
       theta=isign*(6.283185/mmax);
       wtemp=sin(0.5*theta);
       wpr=-2*wtemp*wtemp;
       wpi=sin(theta);
       wr=1;
       wi=0;
       for(m=1;m<mmax;m+=2){
            for(i=m;i<=n;i+=istep){
            j=i+mmax;
            tempr=wr*data[j]-wi*data[j+1];
            tempi=wr*data[j+1]+wi*data[j];
            data[j]=data[i]-tempr;
            data[j+1]=data[i+1]-tempi;
            data[i]+=tempr;
            data[i+1]+=tempi;
        }
      wr=(wtemp=wr)*wpr-wi*wpi+wr;
      wi=wi*wpr+wtemp*wpi+wi;
}
mmax=istep;
   }
 cout<<data[0]<<";"<<data[1]<<";"<<data[2]<<";" <<data[3]; <<";"<<data[4] ;
  }

Le problème est que j'arrive pas a l'utiliser. Alors si vous pouviez m'aider, ca serait bien.
Je ne sais pas comment appeller la fonction four1 et comment obtenir les outputs.

Normalement, on devrait faire comme ca:

Code:

1
2
3
4
5
 
int main(){
float A[4]={3,2,4,5}
four1(A,4,1) 
}

La réponse serait normalement : ([14, -1., 0, -1.]) mais c'est pas ce que j'obtiens.
Je sais que c'est difficile de comprendre un code aussi long sans avoir lu auparavant un pseudo-algo, mais si vous pouvez jetter un coup d'oeil ca serait bien.

09/02/2007, 23h46
thewho

Bonjour,

Il faut aussi lire ce qui est précisé dans NC.

Charge ce fichier pdf:
http://www.ma.utexas.edu/documentati...cpdf/c12-2.pdf

four1 utilise des nombres complexes arrangés dans un tableau de float.
09/02/2007, 23h49
deubelte

j'ai vu merci
mais je crois que leur code ne fonctionne pas.
10/02/2007, 00h08
Matthieu Brucher

La réponse n'est pas forcément celle que tu donnes... Avant de dire que ça ne marche pas, relis un cours sur la FFT. - Je n'ai pas fini mon tuto sur la FFT, mais regarde la partie 1-C http://miles.developpez.com/tutoriels/algo/fft/ -
10/02/2007, 00h43
loufoque

Pourquoi ne pas utiliser FFTW comme tout le monde ?
10/02/2007, 13h24
thewho

Bonjour,

Citation:

Envoyé par deubelte

j'ai vu merci
mais je crois que leur code ne fonctionne pas.

Peut-être que leur code ne fonctionne pas, je ne l'ai pas épluché.

MAIS, au moins, utilise des données au format attendu par la fonction pour tes tests, ce qui n'est pas du tout le cas avec l'exemple que tu nous as donné, ce qui, bien entendu, ne saurait donner un résultat cohérent.

Puisque tu avais vu ce qu'il fallait, pourquoi ne pas avoir utilisé des données correctes ?
10/02/2007, 16h52
deubelte

[CITATION]
Pourquoi ne pas utiliser FFTW comme tout le monde ?
[/CITATION]

parce que je ne sais pas ou ca se trouve.
Ou ca se trouve?

[CITATION]
ce qui n'est pas du tout le cas avec l'exemple que tu nous as donné, ce qui, bien entendu, ne saurait donner un résultat cohérent.
[/CITATION]

ah bonn?? je croyais que j'avais donné le bon format. Quel est le format qu'il faut donner?
float data[]
c'est bien un truc du style:
data[4]={3,2,4,5};
Non?
10/02/2007, 16h52
deubelte

enfin,
quelles sont les balises pour les citations?
10/02/2007, 17h15
loufoque

Tu tapes FFTW dans un moteur de recherche et magie, il te dit où ça se trouve.
Par contre je croyais que c'était LGPL, pas GPL. Ça peut donc être problématique.
10/02/2007, 18h56
thewho

Citation:

Envoyé par deubelte

ah bonn?? je croyais que j'avais donné le bon format. Quel est le format qu'il faut donner?
float data[]
c'est bien un truc du style:
data[4]={3,2,4,5};
Non?

Oui,la déclaration correspond bien à ce que tu fais,

MAIS tu as écrit plus haut que tu avais vu le papier que je t'ai signalé, et le format, pour 4 valeurs, ce n'est pas ce que tu utilises, j'ai donc un sérieux doute sur ton affirmation.

D'ailleurs, je suis même étonné que tu aies pu écrire

Citation:

Envoyé par deubelte

La réponse serait normalement : ([14, -1., 0, -1.]) mais c'est pas ce que j'obtiens.

Car, avec l'appel comme tu dis l'avoir fait,

Citation:

Envoyé par deubelte

int main(){
float A[4]={3,2,4,5}
four1(A,4,1)
}

tu devrais avoir un joli "access violation", car avec le 2ème paramètre = 4, il faut un tableau d'au moins 2*4 = 8 éléments.

Relis le papier, et pas seulement le code.

Et j'ajoute : certes, il est presque certain qu'il y a des erreurs dans les codes de Numerical Recipes, mais avant d'accuser ce code, il est nécessaire de l'utiliser correctement.
10/02/2007, 20h32
olivier1978

Salut,

Il me semble qu'il y a aussi de la FFT dans la GSL. Qui elle, est GPL ...
http://www.gnu.org/software/gsl/manu...ier-Transforms
11/02/2007, 09h33
grob1212

J'ai utilisé la fft de NR il y a 7 ans, et je confirme qu'elle fonctionne parfaitement bien (comparaison effectuée des résultats avec ceux fournis par matlab).