IdentifiantMot de passe
Loading...
Mot de passe oublié ?Je m'inscris ! (gratuit)
Navigation

Inscrivez-vous gratuitement
pour pouvoir participer, suivre les réponses en temps réel, voter pour les messages, poser vos propres questions et recevoir la newsletter

C Discussion :

probleme implémentation algorithme FFT


Sujet :

C

  1. #1
    Membre habitué
    Profil pro
    Inscrit en
    Mars 2007
    Messages
    10
    Détails du profil
    Informations personnelles :
    Localisation : France

    Informations forums :
    Inscription : Mars 2007
    Messages : 10
    Par défaut probleme implémentation algorithme FFT
    Bonjour,

    Je suis en train d'étudier la FFT. J'ai crée un fichier "wav" pour un signal sinusoidal de 1KHz (16 bits signé, échantiolloné à 44100). J'ai supprimé ensuite le "header" 44octets. Le résultat un fichier PCM que j'ai appellé "1khz.data". Les données dans ce fichier sont organisée comme ça (par exemple):
    [C011 0000 2324 0000 ......] le codage est bien sûr little-indian (lsb msb).

    J'ai utilisé l'algorithme de "Numerical recipes" que j'ai adapté un peu (changement très mineur quoi).
    Voici mon code:

    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
    36
    37
    38
    39
    40
    41
    42
    43
    44
    45
    46
    47
    48
    49
    50
    51
    52
    53
    54
    55
    56
    57
    58
    59
    60
    61
    62
    63
    64
    65
    66
    67
    68
    69
    70
    71
    72
    73
    74
    75
    76
    77
    78
    79
    80
    81
    82
    83
    84
    85
    86
    87
    88
    89
    90
    91
    92
    93
    94
    95
    96
    97
    98
    99
    100
    101
    102
    103
    104
    105
    106
    107
    108
    109
    110
    111
    112
    113
    114
    115
    116
    117
    #include <stdio.h>
    #include <stdlib.h>
    #include <math.h>
     
    void FFT(int number);
    #define SWAP(a,b) tempr=(a);(a)=(b);(b)=tempr
     
     
    double tableau[2048]; // tableau pour 1024 points (les imagineurs =0)
     
     
     
     
    void main()
    {
      FILE *fichier=NULL;
      long lsb, msb, i;
      long  data, amplitude;
     
      fichier=fopen("1Khz.data", "rb"); 
     
     
      for (i=0; i<2048; i++)
      {data=fgetc(fichier);
       lsb=data; //codage little indian
       data=fgetc(fichier);
       msb=data;
       data=msb*256+lsb; // la donnee 
     
       //conversion signé nonsigné
        data=data-32768;
       if (data<0) data =data+65536;
     
     
       tableau[i]=data; // ranger dans le tableau
       }
     
      FFT(1024); //on fait la FFT!
     
      i=0;
      while (i<1024)  // la moitié de la tableau (fréquences positives)
      { 
    	  amplitude =sqrt(tableau[i]*tableau[i]+tableau[i+1]*tableau[i+1]);
    	printf("amplitude de fréquence No %ld  est: %ld\n", i, amplitude);
    	i=i+2;
      }
     
     
      fclose(fichier);
      system("PAUSE");    
     
    }
     
     
     
     
    void FFT( int number)
    {   
        long n,mmax,m,j,istep,i;
        long wtemp,wr,wpr,wpi,wi,theta,tempr,tempi;
     
        n=number* 2; 
     
     
    j=0;
        for (i=0;i<n/2;i+=2) {
            if (j > i) {
     
                SWAP(tableau[j],tableau[i]);
     
                SWAP(tableau[j+1],tableau[i+1]);
     
                if((j/2)<(n/4)){
     
                    SWAP(tableau[(n-(i+2))],tableau[(n-(j+2))]);
     
                    SWAP(tableau[(n-(i+2))+1],tableau[(n-(j+2))+1]);
                }
            }
            m=n/2;
            while (m >= 2 && j >= m) {
                j -= m;
                m = m/2;
            }
            j += m;
        }
     
     
        mmax=2;
     
        while (n > mmax)
        {
            istep = mmax<<  1;
            theta=(2*3.141592653589793/mmax);
            wtemp=sin(0.5*theta);
            wpr = -2.0*wtemp*wtemp;
            wpi=sin(theta);
            wr=1.0;
            wi=0.0;
            //internal loops
            for (m=1;m<mmax;m+=2) {
                for (i= m;i<=n;i+=istep) {
                    j=i+mmax;
                    tempr=wr*tableau[j-1]-wi*tableau[j];
                    tempi=wr*tableau[j]+wi*tableau[j-1];
                    tableau[j-1]=tableau[i-1]-tempr;
                    tableau[j]=tableau[i]-tempi;
                    tableau[i-1] += tempr;
                    tableau[i] += tempi;
                }
                wr=(wtemp=wr)*wpr-wi*wpi+wr;
                wi=wi*wpr+wtemp*wpi+wi;
            }
            mmax=istep;
        }
     
    }
    en principe j'obtiens l'amplitude maximal vers les indices No 22, 23, 24 (44100/1024=43Hz, 1000/43=23). l'amplitude sera nul sur les autres indices. Hélas, ça marchait pas!! J'ai obtenu même des amplitudes négatifs!

    merci d'avance pour votre aide

  2. #2
    Expert confirmé

    Profil pro
    Inscrit en
    Janvier 2007
    Messages
    10 610
    Détails du profil
    Informations personnelles :
    Âge : 67
    Localisation : France

    Informations forums :
    Inscription : Janvier 2007
    Messages : 10 610
    Billets dans le blog
    2
    Par défaut
    deux choses de prime abord...

    1) la lisbilité et l'assurance dans les fonctions mathématiques :

    ne jamais présumer de l'ordre d'exécution des opérations, et donc mettre les parenthèses partout où c'est nécessaire.

    Exemple :

    Code : Sélectionner tout - Visualiser dans une fenêtre à part
    1
    2
        amplitude =sqrt((tableau[i]*tableau[i])+(tableau[i+1]*tableau[i+1]));
    Ensuite, sqrt renvoie un double, que tu mets allégrement sans cast dans un long.. Pas trop étonnant que ce soit négatif....

    tableau[i] = (double)data

    amplitude = (long)(srqrt(..)+0.5) si entier le plus proche ou (long)sqrt(..) si troncature.

  3. #3
    Expert confirmé
    Avatar de diogene
    Homme Profil pro
    Enseignant Chercheur
    Inscrit en
    Juin 2005
    Messages
    5 761
    Détails du profil
    Informations personnelles :
    Sexe : Homme
    Localisation : France, Essonne (Île de France)

    Informations professionnelles :
    Activité : Enseignant Chercheur
    Secteur : Enseignement

    Informations forums :
    Inscription : Juin 2005
    Messages : 5 761
    Par défaut
    Code : Sélectionner tout - Visualiser dans une fenêtre à part
     long wtemp,wr,wpr,wpi,wi,theta,tempr,tempi;
    toutes ces variables doivent être des double
    Code : Sélectionner tout - Visualiser dans une fenêtre à part
    1
    2
    3
       //conversion signé nonsigné
        data=data-32768;
       if (data<0) data =data+65536;
    Pourquoi les données doivent elles être non signées ? En fait, il faut les rendre signées

  4. #4
    Membre habitué
    Profil pro
    Inscrit en
    Mars 2007
    Messages
    10
    Détails du profil
    Informations personnelles :
    Localisation : France

    Informations forums :
    Inscription : Mars 2007
    Messages : 10
    Par défaut
    j'ai refut le code avec l'algo original, mais ça marche toujrs pas!!

    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
    36
    37
    38
    39
    40
    41
    42
    43
    44
    45
    46
    47
    48
    49
    50
    51
    52
    53
    54
    55
    56
    57
    58
    59
    60
    61
    62
    63
    64
    65
    66
    67
    68
    69
    70
    71
    72
    73
    74
    75
    76
    77
    78
    79
    80
    81
    82
    83
    84
    85
    86
    87
    88
    89
    90
    91
    92
    93
    94
    95
    96
    97
    #include <stdio.h>
    #include <stdlib.h>
    #include <math.h>
     
     
     
    #define SWAP(a,b) tempr=(a);(a)=(b);(b)=tempr
    void FFT(float data[],unsigned long nn,int isign);
     
    float tableau[2050];
     
     
     
     
     
    int main(void)
    {
      FILE *fichier=NULL;
      long i, c;
      double   amplitude;
      long data1;
     
      fichier=fopen("1Khz.data", "rb"); 
     
      for (i=1; i<2049; i++)
      {  c = fgetc(fichier);
         data1 |= ((c & 0xFF) << (8 * 0));
         c = fgetc(fichier);
         data1 |= ((c & 0xFF) << (8 * 1));
         tableau[i] = data1;
      } 
     
      FFT(tableau,1024,1);
     
      i=1;
      while (i<1024)
      {   amplitude=sqrt((tableau[i]*tableau[i])+(tableau[i+1]*tableau[i+1]));
      printf("l'amplitude de la fréquence %ld est: %ld\n", i, amplitude);
      i=i+2;
      }
     
     
      fclose(fichier);
      system("PAUSE");    
     
    }
     
     
     
     
     
    void FFT(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;
           theta=isign*(6.28318530717959/mmax);
           wtemp=sin(0.5*theta);
           wpr=-2*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){
                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;
       }
     
      }
    A L'AIIIIDE!!

  5. #5
    Expert confirmé

    Profil pro
    Inscrit en
    Janvier 2007
    Messages
    10 610
    Détails du profil
    Informations personnelles :
    Âge : 67
    Localisation : France

    Informations forums :
    Inscription : Janvier 2007
    Messages : 10 610
    Billets dans le blog
    2
    Par défaut
    Primo pourquoi mélanger des float et des doubles ??

    Tu prends l'un, ou tu prends l'autre.. En général, en C, pour les maths, c'est mieux de prendre des doubles, pour plusieurs raisons, dont le fait que c'est le type "naturel" de C pour les nombres flottants (si tu passes un float à une fonction, le compilateur le promeut automatiquement en double à l'entrée et le déprécie en float à la sortie), et d'autre part toutes les fonctions de maths utilisent des double par défaut. Il existe les versions en float, mais comme de toutes façons c'est le type par défaut, c'est mieux d'utiliser des doubles, sauf si tu as des problèmes spécifiques de place mémoire ET que tu es sûr que tes nombres ne dépassent pas FLOAT_MAX.

    Secondo, dans ton printf après amplitude, tu imprimes un double (tu as transformé amplitude) avec un format entier (tu as oublié de changer le format). Donc je te recommande %g ou %lf ..


    Tertio soit cohérent aussi dans tes notations dans la FFT. Des fois tu mets 1.0, d'autres fois -2.. Mets partout des .0. Et de plus, quand tu fais une opération mixant int et double, fait le cast du int vers double, surtout si il y a une division en jeu (voir theta et wpr).

    Enfin dans l'algo de FFT, et bien que ça ne change rien, il y a encore 2 lignes où les parenthèses ne sont pas complètes (tempr et tempi).

    Et dans la ligne sur wr, essayes de ne pas utiliser d'écriture cryptique comme :

    Code : Sélectionner tout - Visualiser dans une fenêtre à part
    1
    2
     
          wr=(wtemp=wr)*wpr-wi*wpi+wr;
    dont je ne sais même pas si ça marche, et déjà contente-toi de faire du code qui marche.

    Une fois qu'il marche, tu pourras éventuellement penser à l'optimiser.

  6. #6
    Expert confirmé
    Avatar de diogene
    Homme Profil pro
    Enseignant Chercheur
    Inscrit en
    Juin 2005
    Messages
    5 761
    Détails du profil
    Informations personnelles :
    Sexe : Homme
    Localisation : France, Essonne (Île de France)

    Informations professionnelles :
    Activité : Enseignant Chercheur
    Secteur : Enseignement

    Informations forums :
    Inscription : Juin 2005
    Messages : 5 761
    Par défaut
    Le tableau est mal construit:
    Code : Sélectionner tout - Visualiser dans une fenêtre à part
    1
    2
    3
    4
    5
    6
    7
    8
    9
      for (i=1; i<2049; i++)
      {  c = fgetc(fichier);
         data1 |= ((c & 0xFF) << (8 * 0));
    //  data1 = (c & 0xFF)
         c = fgetc(fichier);
         data1 |= ((c & 0xFF) << (8 * 1));
         tableau[i] = data1;
    // Il faut rétablir le signe dans data1 avant le rangement dans le tableau
      }
    Une fonction doit être correctement testée avant d'être utilisée : Pour tester la FFT, il est préférable de remplir le tableau par une sinusoïde et de vérifier la sortie de la FFT. Les sinusoïdes de test doivent faire exactement une, deux,... périodes dans le tableau de façon à avoir un spectre d'une raie dont la position et l'amplitude attendues sont parfaitement connues. Sinon, on a un mélange des problèmes de lecture du fichier et de vérification de la FFT qui sont imbriqués.

  7. #7
    Membre habitué
    Profil pro
    Inscrit en
    Mars 2007
    Messages
    10
    Détails du profil
    Informations personnelles :
    Localisation : France

    Informations forums :
    Inscription : Mars 2007
    Messages : 10
    Par défaut
    Merci diogene, souviron34, Emmanuel

    oui mon tableau contient une seule sinusoide de 1khz échantillonée à 44100. Pour 1024 échantillons, en principe je dois avoir l'amplitude maximal à l'indice No 23 du tableau ( 44100/1024 =43, 1000/43=23).

    Pourquoi, stp, le tableau est mal construit?
    Les données PCM dans mon fichier sont organisées de cette façon:

    [C011 0000 2324 0000 ..........] le codage est little-indian donc lsb=C0, msb=11 ...etc, chaque échantillon contient deux parties, chacune sur 16 bits signés. Les 0000 sont la partie imaginaire de l'échantillon.

    Est-ce que il faut envoyer les données à la FFT tels quels ou il faut faire une quelconque conversion avant?
    Les résultats des calculs de la FFT seront elles dans le même tableau de données de départ (entrées) ou dans un autre tableau?

    J'ai rectifié les types des variables selon les idications de souviron34. Maintenant le programme affiche des FLOTS positifs, donc c'est déjà pas mal!! mais je nai pas encore trouvé un amplitude maximal à l'indice 23!!

    Merci

  8. #8
    Expert confirmé
    Avatar de diogene
    Homme Profil pro
    Enseignant Chercheur
    Inscrit en
    Juin 2005
    Messages
    5 761
    Détails du profil
    Informations personnelles :
    Sexe : Homme
    Localisation : France, Essonne (Île de France)

    Informations professionnelles :
    Activité : Enseignant Chercheur
    Secteur : Enseignement

    Informations forums :
    Inscription : Juin 2005
    Messages : 5 761
    Par défaut
    Pourquoi, stp, le tableau est mal construit?
    Data1 est un long. Tu construit sa valeur par des "OU" sur 16 bits. Comme tes long font plus de 16 bits, les poids forts sont à 0, et les valeurs obtenus sont toutes positives : il faut étendre le bit de signe sur les poids forts.
    Le premier |= sur data1 est de plus faux
    Est-ce que il faut envoyer les données à la FFT tels quels ou il faut faire une quelconque conversion avant?
    Il faut tenir compte des particularités de TA FFT : Est ce que celle-ci traite des signaux réels ou des signaux complexes (dans ce cas, où doivent être placés les parties imaginaires)
    Les résultats des calculs de la FFT seront elles dans le même tableau de données de départ (entrées) ou dans un autre tableau?
    La plupart écrase le tableau d'entrée avec le résultat du calcul. Ceci est à l'évidence ton cas puisque le prototype de ta FFT ne comporte qu'un tableau. L'ordre des fréquences dans le tableau peut par contre varier ( de -f/2 à +f/2 ou de 0 à f/2 puis -f/2 à 0).

    Teste la FFT sur un signal de synthèse. Cela te permettra d'en déduire facilement le mode d'emploi de TA fft. Celle-ci mise au point, tu pourras passer au problème du remplissage correct du tableau.

    Bon courage

  9. #9
    Membre habitué
    Profil pro
    Inscrit en
    Mars 2007
    Messages
    10
    Détails du profil
    Informations personnelles :
    Localisation : France

    Informations forums :
    Inscription : Mars 2007
    Messages : 10
    Par défaut
    Merci diogene,

    ça marche impec avec un signal de synthese!
    l'algo de la FFT est celui de Numerical recipes un peu modifié
    http://www.codeproject.com/cpp/howtofft.asp
    voici le code:

    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
    36
    37
    38
    39
    40
    41
    42
    43
    44
    45
    46
    47
    48
    49
    50
    51
    52
    53
    54
    55
    56
    57
    58
    59
    60
    61
    62
    63
    64
    65
    66
    67
    68
    69
    70
    71
    72
    73
    74
    75
    76
    77
    78
    79
    80
    81
    82
    83
    84
    85
    86
    87
    88
    89
    90
    91
    92
    93
    94
    95
    96
    97
    98
    99
    100
    101
    102
    103
    104
    105
    106
    107
    108
    109
    110
    111
    112
    113
    114
    115
    116
    117
    118
    119
    120
    121
    122
    123
    124
    125
    #include <stdio.h>
    #include <stdlib.h>
    #include <math.h>
    #define PI 3.1415926535897932385
    #define SWAP(a,b) tempr=(a);(a)=(b);(b)=tempr
     
    void FFT(float data[],unsigned long N,int isign);
     
     
    int main()
    { 
    	int i, j=5; // j pour le signal carre
    	float sample, amplitude, A[2048];
     
     
     
      /* géneration d'un signal sinusoidal 1000hz echantilloné à 44100*/
     
     
    	for (i=0; i<1024; i++)
    	{ sample= 32000*sin((i/44.1)*2*PI); // 44100:1000=44.1 (44.1 echantillons par periode)
    	  A[2*i]=sample;
    	  A[2*i+1]=0;
    	}
     
     
    	 /* *******  ou un signal  rectangle  (plusieurs harmoniques) 
     
     
     
      for( i = 0 ; i < 1024 ; i++)
        {
          if( i % 44 == 0 ) // 44 echantillons par période
    	{
    	  if ( j == 5 ) // amplitude = 5 volt
    	    j = -5;
    	  else
    	    j = 5;
    	}
          A[2*i] = j;
    	  A[2*i+1] =0;
     
        }     *******  */
     
     
     
    FFT(A,1024,1);  
     
     
    /* Affiche la TF du signal d'entrée */
     
      i = 0 ; 
      while ( i < 256)    //les 256 premières fréquences
        {   
            amplitude = sqrt(A[i]*A[i]+A[i+1]*A[i+1]);
    		printf(" l'amplitude à l'indice %ld est: %f \n",i, amplitude);
    		i=i+2;
        }
     
    system("PAUSE");
    return 0;
    }
     
     
     
    void FFT(float data[], unsigned long N, int isign)
    {
     
        unsigned long n,mmax,m,j,istep,i;
        double wtemp,wr,wpr,wpi,wi,theta,tempr,tempi;
     
        n=N * 2; 
     
     
        j=0;
        for (i=0;i<n/2;i+=2) {
            if (j > i) {
     
                SWAP(data[j],data[i]);
     
                SWAP(data[j+1],data[i+1]);
     
                if((j/2)<(n/4)){
     
                    SWAP(data[(n-(i+2))],data[(n-(j+2))]);
     
                    SWAP(data[(n-(i+2))+1],data[(n-(j+2))+1]);
                }
            }
            m=n/2;
            while (m >= 2 && j >= m) {
                j -= m;
                m = m/2;
            }
            j += m;
        }
     
        mmax=2;
        //external loop
        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;
            //internal loops
            for (m=1;m<mmax;m+=2) {
                for (i= m;i<=n;i+=istep) {
                    j=i+mmax;
                    tempr=wr*data[j-1]-wi*data[j];
                    tempi=wr*data[j]+wi*data[j-1];
                    data[j-1]=data[i-1]-tempr;
                    data[j]=data[i]-tempi;
                    data[i-1] += tempr;
                    data[i] += tempi;
                }
                wr=(wtemp=wr)*wpr-wi*wpi+wr;
                wi=wi*wpr+wtemp*wpi+wi;
            }
            mmax=istep;
        }
    }
    (J'ai essayé aussi l'algo Cooley-Tuckey ça marche bien).

    par contre j'arrive toujours pas à faire la FFT de mon fichier dont j'ai parlé au début de la discussion!!!!

  10. #10
    Expert confirmé
    Avatar de diogene
    Homme Profil pro
    Enseignant Chercheur
    Inscrit en
    Juin 2005
    Messages
    5 761
    Détails du profil
    Informations personnelles :
    Sexe : Homme
    Localisation : France, Essonne (Île de France)

    Informations professionnelles :
    Activité : Enseignant Chercheur
    Secteur : Enseignement

    Informations forums :
    Inscription : Juin 2005
    Messages : 5 761
    Par défaut
    ça marche impec avec un signal de synthese!
    dans ce cas, les données d'entrée sont mal construites. Pour un 16 bits en PCM , sauf erreur de ma part, les données sont en complément à 2. Si tu les lis dans un long, initialisé à zéro, il faut restituer le signe sur le long avant de le ranger dans le tableau de float. La FFT que tu utilises semble classer dans le tableau la partie réelle puis la partie imaginaire. L'élément suivant du tableau doit être mis à 0.
    Code : Sélectionner tout - Visualiser dans une fenêtre à part
    1
    2
    3
    4
    5
    6
    7
    8
      for (i=0; i<2048; i +=2)
      { 
        data1 = fgetc(fichier);
        data1 = data1 +256*fgetc(fichier);
        if(data1>32767) data1 -= 65536;
        tableau[i] = data1;
        tableau[i+1]= 0;
      }

  11. #11
    Nouveau candidat au Club
    Inscrit en
    Avril 2007
    Messages
    3
    Détails du profil
    Informations forums :
    Inscription : Avril 2007
    Messages : 3
    Par défaut besoin d'aide
    [-mod- Merci de relire les règles du forum. Il est considéré comme tout à fait grossier de squatter ainsi un fil de discussion de quelqu'un d'autre, même si il est en relation avec ta question. Le principe, c'est un fil par question. Recréé un nouveau sujet, sinon, je serais obligé de supprimer tes posts de ce fil.]
    bonjour!
    je suis entrain d'essayer de programmer la fft ki est sur le livre numerica recepes.moi mon but c juste donner un petit signal et voir si ca marche .donc j'ai ajouter un programme pricipal(main).mais j'ai tjr un probleme ,d'apres ce que g compris ,il suffit de donner des valeur pr isign,exempl 1 por fft directe et -1 inverse???mais apres la simulation -1 me donne NAN(not a number).quelq'un peux me dir ou est le probleme dans mon programme pricipal ou bien m'aider!!ici g choisis une fonction sinus.g vraiment besoin d'aide.c urgent

    [-mod- Merci de réécrire ceci dans un français lisible. L'usage du SMS et autre abréviations abusives est interdit]

    merci d'avance
    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
    36
    37
    38
    39
    40
    41
    42
    43
    44
    45
    46
    47
    48
    49
    50
    51
    52
    53
    54
    55
    56
    57
    58
    59
    60
    61
    62
    63
    64
    65
    66
    67
    68
    #include<math.h>
    #include<stdio.h>
    #include<stdlib.h>
     
    #define swap(a,b)tempr=(a);(a)=(b);(b)=tempr
     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=n>>1;
           while(m>=2 && j>m)
           {
                      j-=m;
                      m>>=1;
                      }
                      j+=m;
                      }
                      mmax=2;
                      while (n>mmax)
                      {
                            istep=mmax<<1;
                            theta=isign*(6.2831/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)
                            {
                                                  j=i+mmax;
                                                  tempr=wr*data[j]-wi*data[j+1];
                                                  tempr=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;
                           }
                           }
     
      int main  ()
    {   
         int i,nn,isign;
    float x[3];
    nn=4;
         for(i=0;i<3;i++)
         x[i]=sin(2*3.14*i);
         four1(x,nn,isign);
         for(i=0;i<3;i++)
          printf("data[%d]=%f\n",i,x[i]);
     
         }
    [-mod- Je rappelle que le code doit être posté entre des balises de code : ]

    le nombre de point =10,et le nombre d'echantillon=nn=8

  12. #12
    Nouveau candidat au Club
    Inscrit en
    Avril 2007
    Messages
    3
    Détails du profil
    Informations forums :
    Inscription : Avril 2007
    Messages : 3
    Par défaut besoin d'aide
    bonjour!
    je suis entrain d'essayer de programmer la fft ki est sur le livre numerica recepes.moi mon but c juste donner un petit signal et voir si ca marche .donc j'ai ajouter un programme pricipal(main).mais j'ai tjr un probleme ,d'apres ce que g compris ,il suffit de donner des valeur pr isign,exempl 1 por fft directe et -1 inverse???mais apres la simulation -1 me donne NAN(not a number).quelq'un peux me dir ou est le probleme dans mon programme pricipal ou bien m'aider!!ici g choisis une fonction sinus.g vraiment besoin d'aide.c urgent
    merci d'avance
    #include<math.h>
    #include<stdio.h>
    #include<stdlib.h>

    #define swap(a,b)tempr=(a);(a)=(b);(b)=tempr
    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=n>>1;
    while(m>=2 && j>m)
    {
    j-=m;
    m>>=1;
    }
    j+=m;
    }
    mmax=2;
    while (n>mmax)
    {
    istep=mmax<<1;
    theta=isign*(6.2831/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)
    {
    j=i+mmax;
    tempr=wr*data[j]-wi*data[j+1];
    tempr=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;
    }
    }

    int main ()
    {
    int i,nn,isign;
    float x[3];
    nn=4;
    for(i=0;i<3;i++)
    x[i]=sin(2*3.14*i);
    four1(x,nn,isign);
    for(i=0;i<3;i++)
    printf("data[%d]=%f\n",i,x[i]);

    }


    le nombre de point =3,et le nombre d'echantillon=nn=8

  13. #13
    Expert confirmé

    Profil pro
    Inscrit en
    Janvier 2007
    Messages
    10 610
    Détails du profil
    Informations personnelles :
    Âge : 67
    Localisation : France

    Informations forums :
    Inscription : Janvier 2007
    Messages : 10 610
    Billets dans le blog
    2
    Par défaut
    D'abord pense à la balise CODE (#) et pas de langages sms ici s'il te plaît..

    Merci




  14. #14
    Expert confirmé
    Avatar de diogene
    Homme Profil pro
    Enseignant Chercheur
    Inscrit en
    Juin 2005
    Messages
    5 761
    Détails du profil
    Informations personnelles :
    Sexe : Homme
    Localisation : France, Essonne (Île de France)

    Informations professionnelles :
    Activité : Enseignant Chercheur
    Secteur : Enseignement

    Informations forums :
    Inscription : Juin 2005
    Messages : 5 761
    Par défaut
    Il y a des erreurs dans le code de la FFT :
    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
    36
    37
    38
    39
    40
    41
    42
    43
    44
    45
    46
    47
    48
    49
    50
    51
    #define swap(a,b)tempr=(a);(a)=(b);(b)=tempr
    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=0;                                 //j=1
      for (i=0;i<n;i+=2)                   //i=1
      {
        if (j>i)
        {
          swap(data[j],data[i]);
          swap(data[j+1],data[i+1]);
        }
        m=n>>1;
        while(m>=2 && j>=m)                 //j>m
        {
          j-=m;
          m>>=1;
        }
        j+=m;
      }
      mmax=2;
      while (n>mmax)
      {
        istep=mmax<<1;
        theta=isign*(6.2831853/mmax);       // 2*pi
        wtemp=sin(0.5*theta);
        wpr=-2.0*wtemp*wtemp;
        wpi=sin(theta);
        wr=1.0;
        wi=0.0;
        for(m=0;m<mmax;m+=2)                // m=1
        {
          for(i=m;i<n;i+=istep)             //i<=n
          {
            j=i+mmax;
            tempr=wr*data[j]-wi*data[j+1];
            tempi=wr*data[j+1]+wi*data[j];  //tempr
            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;
      }
    }
    Les lignes modifiées sont signalées par un commentaire
    Pour tester :
    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
    #define N 32  // Nombre de points de la FFT.
                  // DOIT être une puissance de 2
    #define F 11  // Fréquence du signal de test.
                  //Positive,  négative ou nulle
    int main (void)
    {
      int i;
      float x[2*N];
      float phi = 3.1415927*F/N;
      for(i=0;i<2*N;i+= 2)
      {
        x[i]=cos(phi*i);    // ou cos(phi*i)  // ou sin(phi*i)
        x[i+1]=sin(phi*i);  // et 0.0;        // et 0.0
      }
      four1(x,N,1);
     
    return 0;
    }
    - si on a x[i]=cos(phi*i) et x[i+1]=sin(phi*i) comme signal de test, une raie d'amplitude N doit apparaitre en position 2*N+2*F si F<0 et en 2*F si F>=0. Le reste du tableau est pratiquement à 0
    - Avec x[i]=cos(phi*i) et x[i+1]=0, deux raies d'amplitude N/2 sont en position 2*|F| et 2*N+2*|F|
    - Avec x[i]=sin(phi*i) et x[i+1]=0, deux raies d'amplitude N/2 et de signe opposé sont en position 2*|F|+1 et 2*N+2*|F|+1

  15. #15
    Futur Membre du Club
    Inscrit en
    Mai 2007
    Messages
    5
    Détails du profil
    Informations forums :
    Inscription : Mai 2007
    Messages : 5
    Par défaut fft inverse
    [CODE][CODE]
    bonjour!!
    cette fft marche tres bien.mais j'arrive toujour pas a la modifier en fft inverse.quelqu'un peux m'aider??me dir ou dois-je changer??comment fair??
    mecri d'avance






    #include <math.h>
    #define PI 3.1415926535897932385
    #define NMAX 8
    #include <stdio.h>
    #include <stdlib.h>
    /* Cette fonction effectue le "bit reversal" sur un indice
    donné
    @param: l'indice sur lequel s'effectue le "bit reversal"
    @param: le nombre de bit dont on a besoin pour écrire chaque
    indice
    @return: l'indice obtenu après le "bit reversal"
    */
    unsigned int tf_reversion(unsigned int indice,int nb)
    {
    unsigned int p, k,i,size;
    size=32;
    p=indice >> nb-1;
    for(i=1;i<nb;i++)
    {
    k=indice << i+size-nb;
    k=k >> size-1;
    k=k << i;
    p+=k;
    }
    return(p);
    }
    /* Cette fonction calcule la fft
    @param: logarithme en base 2 de N(nombres de points)
    @param: valeur utilisé dans l'algorithme
    @param: tableau des échantillons (partie réelle + partie
    imaginaire)
    @param: tableau des calculs (partie réelle + partie
    imaginaire)
    */
    void tf_cooley_tukey(int r,unsigned int n,double
    f[NMAX][2],double t[NMAX][2])
    {
    int k;
    unsigned int ne,ns,n1,n2,i,j,tj;
    double arg,tr1,ti1,tr2,ti2,cr,ci;
    double tcos[NMAX],tsin[NMAX];
    // préparation des calculs (il n'existe pas de complexe en c
    ,on doit traiter les imaginaires et les réels séparément.
    n2=n/2;
    arg=PI/n2;
    for(i=0;i<n2;i++)
    {
    tcos[i]=cos(-arg*i);
    tcos[i+n2]=-tcos[i];
    tsin[i]=sin(-arg*i);
    tsin[i+n2]=-tsin[i];
    }
    for (i=0;i<n;i++)
    {// "bit reversal"
    j=tf_reversion(i,r);
    t[j][0]=f[i][0];
    t[j][1]=f[i][1];
    }
    ne=1;n1=2;ns=n;
    for (k=0;k<r;k++)
    {
    printf("Transformation T%d :\n",k);
    for (i=0;i<n;i++)
    printf("( %12.5e ) + ( %12.5e ) i\n",t[i][0],t[i][1]);
    for (i=0;i<=n-n1;i+=n1)
    for(j=i;j<=i+ne-1;j++)
    {
    tj=j*n2;
    while (tj>=n) tj-=n;
    cr=tcos[tj];
    ci=tsin[tj];
    // opérateur "butterfly" 1ere partie
    tr1=t[j][0]+cr*t[j+ne][0]-ci*t[j+ne][1];
    ti1=t[j][1]+cr*t[j+ne][1]+ci*t[j+ne][0];
    tj=(j+ne)*n2;
    while (tj>=n) tj-=n;
    cr=tcos[tj];
    ci=tsin[tj];
    // opérateur "butterfly" 2eme partie
    tr2=t[j][0]+cr*t[j+ne][0]-ci*t[j+ne][1];
    ti2=t[j][1]+cr*t[j+ne][1]+ci*t[j+ne][0];
    t[j][0]=tr1;t[j][1]=ti1;
    t[j+ne][0]=tr2;t[j+ne][1]=ti2;
    }
    ne=n1;n1=2*ne;ns=n2;n2=ns/2;
    }
    }
    main()
    {
    int r;
    unsigned int n,i;
    double f[NMAX][2],t[NMAX][2];
    printf("Méthode de Cooley-Tukey\n");
    f[0][0]=1;f[0][1]=0;
    f[1][0]=2;f[1][1]=0;
    f[2][0]=3;f[2][1]=0;
    f[3][0]=4;f[3][1]=0;
    f[4][0]=5;f[4][1]=0;
    f[5][0]=6;f[5][1]=0;
    f[6][0]=7;f[6][1]=0;
    f[7][0]=8;f[7][1]=0;
    r=3;
    n=1;
    21
    for(i=1;i<=r;i++)
    n*=2;
    tf_cooley_tukey(r,n,f,t);
    printf("Transformation T%d :\n",r);
    for (i=0;i<n;i++)
    printf("F%u = ( %12.5e ) + ( %12.5e )
    i\n",i,t[i][0],t[i][1]);
    }

  16. #16
    Nouveau candidat au Club
    Inscrit en
    Avril 2007
    Messages
    3
    Détails du profil
    Informations forums :
    Inscription : Avril 2007
    Messages : 3
    Par défaut fft
    bonjour!!
    voici le programme de la fft qui donne des résultats,mais j'arrive toujours pas à faire sa fft inverse, je ne sais pas comment proceder ,quelqu'un pourait me dire où faire le changement.
    Merci de votre aide
    [le programme:
    unsigned int tf_reversion (unsigned int indice,int nb)
    {
    unsigned int p,k,i,size;
    size=32;
    p=indice>>nb-1;
    for(i=1;i<nb;i++)
    {
    k=indice<<i+size-nb;
    k=k>>size-1;
    k=k<<i;
    p+=k;
    }
    return(p);
    }
    void tf_cooley_tukey(int r,unsigned int n,double f[Nmax][2],double t[Nmax][2])
    {
    int k;
    unsigned int ne,ns,n1,n2,i,j,tj;
    double arg,tr1,ti1,tr2,ti2,cr,ci;
    double tcos[Nmax],tsin[Nmax];
    n2=n/2;
    arg=pi/n2;
    for (i=0;i<n2;i++)
    {
    tcos[i]=cos(-arg*i);
    tcos[i+n2]=-tcos[i];
    tsin[i]=sin(-arg*i);
    tsin[i+n2]=-tsin[i];
    }
    for (i=0;i<n;i++)
    {
    j=tf_reversion(i,r);
    t[j][0]=f[i][0];
    t[j][1]=f[i][1];
    }
    ne=1;n2=2;ns=n;
    for (k=0;k<r;k++)
    {
    printf("transformation T%d :\n",k);
    for(i=0;i<n;i++)
    printf("%12.5e)+(%12.5e) i\n",t[i][0],t[i][1]);
    for (i=0;i<=n-n1;i+=n1)
    for(j=i;j<=i+ne-1;j++)
    {
    tj=j*n2 ;
    while (tj>=n) tj-=n;
    cr=tcos[tj];
    ci=tsin[tj];
    tr1=t[j][0]+cr*t[j+ne][0]-ci*t[j+ne][1];
    ti1=t[j][1]+cr*t[j+ne][1]+ci*t[j+ne][0];
    tj=(j+ne)*n2;
    while (tj>=n) tj-=n;
    cr=tcos[tj];
    ci=tsin[tj];
    tr2=t[j][0]+cr*t[j+ne][0]-ci*t[j+ne][1];
    ti2=t[j][1]+cr*t[j+ne][1]+ci*t[j+ne][0];
    t[j][0]=tr1;t[j][1]=ti1;
    t[j+ne][0]=tr2;t[j+ne][1]=ti2;
    }
    ne=n1;n1=2*ne;ns=n2;n2=ns/2;
    }
    }
    main( )
    {
    int r;
    unsigned int n,i;
    double f[Nmax][2],t[Nmax][2];
    printf("methode de cooley-tukey \n");
    f[0][0]=1;f[0][1]=0;
    f[1][0]=2;f[1][1]=0;
    f[2][0]=3;f[2][1]=0;
    f[3][0]=4;f[3][1]=0;
    f[4][0]=5;f[4][1]=0;
    f[5][0]=6;f[5][1]=0;
    f[6][0]=7;f[6][1]=0;
    f[7][0]=8;f[7][1]=0;
    r=3;
    n=1;
    for(i=1;i<=r;i++)
    n*=2;
    tf_cooley_tukey(r,n,f,t);
    printf("transformation T%d:\n",r);
    for(i=0;i<n;i++)
    printf("F%u=(%12.5e)+(%12.5e)i\n",i,t[i][0],t[i][1]);
    }
    ]

Discussions similaires

  1. Algorithme FFT - Un petit coup de pouce?
    Par curlcie dans le forum Débuter
    Réponses: 0
    Dernier message: 26/10/2010, 20h04
  2. probleme avec la FFT
    Par Iven1 dans le forum C
    Réponses: 9
    Dernier message: 15/06/2009, 20h09
  3. vitesse d'un algorithme FFT
    Par ABN84 dans le forum Algorithmes et structures de données
    Réponses: 1
    Dernier message: 03/05/2009, 13h12
  4. Probleme implémentation bd
    Par Bardack dans le forum Hibernate
    Réponses: 1
    Dernier message: 21/03/2007, 15h20
  5. Implémenter un fft 2d
    Par jeje99 dans le forum Traitement du signal
    Réponses: 7
    Dernier message: 28/01/2005, 13h36

Partager

Partager
  • Envoyer la discussion sur Viadeo
  • Envoyer la discussion sur Twitter
  • Envoyer la discussion sur Google
  • Envoyer la discussion sur Facebook
  • Envoyer la discussion sur Digg
  • Envoyer la discussion sur Delicious
  • Envoyer la discussion sur MySpace
  • Envoyer la discussion sur Yahoo