Fast fourier transform

Soyez le premier à donner votre avis sur cette source.

Vue 21 656 fois - Téléchargée 1 518 fois

Télécharger le projet
Commenter

Description

Un petit exemple de la FAST FOURIER TRANSFORM sur la fonction f(x) = x*(1-x).

N'hésitez pas à mettre des commentaires. Si vous avez des idées pour améliorer mon code n'hésitez pas non plus.

Source / Exemple : 


// FAST FOURIER TRANSFORM
// Exemple de la FFT sur la fonction f(x) = x*(1-x)

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#define  swap(a,b)  norm=(a); (a)=(b); (b)=norm

//reel[] et imag[i] sont la liste des réelles et des imaginaires
// sign = 1 donne la transformée de Fourier
// sign = -1 donne la transformée de Fourier inverse

void fft(double *reel, double *imag, int log2n, int sign) {
  
  int n, m, m2, i, j, k, l;
  double  c1, c2, norm, norm2, cphi, sphi;

  n = 1<<log2n;

  /* Inversement des bits */
  for(i=0; i<n; i++) {
    
    for(j=log2n-1, m=0, k=i; j>=0; j--, k>>=1) m += (k&1)<<j;
    
    if(m>i) {
      swap(reel[i],reel[m]);
      swap(imag[i],imag[m]);
    }
  }
  
  /* normalisation de la transformée de Fourier */
  norm = 1.0/sqrt((double)n);
  for(i=0; i<n ;i++) {
    reel[i] *= norm;
    imag[i] *= norm;
  }
  
  /* calcul de la FFT */
  for(j=0; j < log2n; j++) {
    m = 1<<j;  m2 = 2*m;
    c1 = 1.0; 
    c2 = 0.0;
    cphi = cos(sign*2.0*M_PI/((double)m2));
    sphi = sin(sign*2.0*M_PI/((double)m2));
    for(k=0; k<m; k++) {
      for(i=k; i<n; i+=m2) {
	l = i + m;
	norm  = c1*reel[l] - c2*imag[l];
	norm2 = c1*imag[l] + c2*reel[l];
	reel[l] = reel[i] - norm;
	imag[l] = imag[i] - norm2;
	reel[i] += norm;
	imag[i] += norm2;
      }
      norm  = c1*cphi - c2*sphi; // Calcul de exp(2 pi k/m) avec
      norm2 = c1*sphi + c2*cphi; // le théorème d'addition
      c1 = norm;  c2 = norm2;
    }
  }
  
}
int main(int argc, char *argv[])
{
  int n=16, k=4, i;
  double re[16], im[16];
  double h,x;
  FILE *fichier;
  fichier = fopen("FFT.dat","w");
  
  h = 1.0/(n-1);
  
  printf("Calcul des Points:\n");
  for(i=0; i<n; i++) {
    x = h*i;
    re[i] = x*(1-x);
    im[i] = 0;
    
    printf(" % lf %+lf i\n",re[i],im[i]);
    
  }
  
  fft(re,im,k,+1);
  printf("Transformation:\n");
  for(i=0; i<n; i++) {
    printf(" % lf %+lf i\n",re[i],im[i]);
    fprintf(fichier,"%d %lf %lf\n",i,re[i],im[i]);//on enregistre dans FFT.dat
  }
  fclose(fichier);
  
  fft(re,im,k,-1); 
  printf("FFT inverse:\n");
  for(i=0; i<n; i++) {
  printf(" % lf %+lf i\n",re[i],im[i]);
  }
  system("PAUSE");	
  return 0;
}
Télécharger le projet

Codes Sources

A voir également

Ajouter un commentaire

Commentaires

Signaler
luhtor
Messages postés
2023
Date d'inscription
mardi 24 septembre 2002
Statut
Membre
Dernière intervention
28 juillet 2008
 4 -
La FFT, elle a besoin de calculer les coefficients de fourrier ?
Pour les calculer, c'est un calcul d'intégrale, et j'arrive pas à comprendre ou est ce que tu fais ce calcul.
J'ai du mal à comprendre la fonction, tu pourrais détailler un peu ? car j'aimerais bien l'utiliser, mais je voudrais d'abord savoir comment elle fonctionne pour pouvoir l'adapter.

Merci
Signaler
Arnaud16022
Messages postés
1329
Date d'inscription
vendredi 15 août 2003
Statut
Membre
Dernière intervention
16 juin 2010
 2 -
idem...
le systeme est utilisable sur n'importe quelle fonction ou crée/optimisé spécialement pour celle la?
Signaler
Jarod1980
Messages postés
273
Date d'inscription
samedi 5 juillet 2003
Statut
Membre
Dernière intervention
31 mars 2015
 1 -
Pour ce qui est de l'explication de la FFT je vais mettre dès que j'ai un peu de temps un PDF explicatif avec toutes les formules. Cette routine FFT est utilisable avec n'importe quelle fonction. Ici j'ai choisis comme exemple f(x) = x*(1-x) mais vous pouvez utiliser une autre fonction, par exemple :

sin(2*M_PI*3*i/((double)n)) + 1.0;
Signaler
Arnaud16022
Messages postés
1329
Date d'inscription
vendredi 15 août 2003
Statut
Membre
Dernière intervention
16 juin 2010
 2 -
yeah
vive adobe
je l'attends avec impatience
:)
Signaler
Arnaud16022
Messages postés
1329
Date d'inscription
vendredi 15 août 2003
Statut
Membre
Dernière intervention
16 juin 2010
 2 -
je viens de dl le pdf (c pas trop tot..)
c'est bien galère
je vais regarder ca mais les intégrales et moi ca fait au mois 10
merci qd meme
ad
Afficher les 10 commentaires
S'inscrire maintenant

Vous n'êtes pas encore membre ?

inscrivez-vous, c'est gratuit et ça prend moins d'une minute !

Les membres obtiennent plus de réponses que les utilisateurs anonymes.

Le fait d'être membre vous permet d'avoir un suivi détaillé de vos demandes et codes sources.

Le fait d'être membre vous permet d'avoir des options supplémentaires.

Du même auteur (Jarod1980)