probleme implémentation algorithme FFT

**philo69** · 26/03/2007, 15h18

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

**souviron34** · 26/03/2007, 16h10

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.

**diogene** · 26/03/2007, 19h20

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

**philo69** · 26/03/2007, 21h33

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!!

**souviron34** · 27/03/2007, 01h00

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.

**diogene** · 27/03/2007, 13h32

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.

**jako2** · 30/04/2007, 16h53

[-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

**jako2** · 30/04/2007, 16h58

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

**souviron34** · 30/04/2007, 17h33

D'abord pense à la balise CODE (#) et pas de langages sms ici s'il te plaît..

Merci

**diogene** · 02/05/2007, 13h54

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

**sofiya** · 07/05/2007, 22h01

[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]);
}

**jako2** · 08/05/2007, 17h33

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]);
}
]

probleme implémentation algorithme FFT

C

Vue hybride

Discussions similaires

Partager

Partager