Fast fourier transform

Soyez le premier à donner votre avis sur cette source.

Vue 23 025 fois - Téléchargée 1 635 fois

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

Codes Sources

A voir également

Ajouter un commentaire

Commentaires

Messages postés
475
Date d'inscription
vendredi 27 juin 2003
Statut
Membre
Dernière intervention
20 septembre 2010

C'est quoi k=4 ?
Est-ce qu'on peut trouver la fondamentale dun signal avec cette source et comment?
Messages postés
1
Date d'inscription
mardi 7 mars 2006
Statut
Membre
Dernière intervention
8 mars 2006

Est ce que quelqu'un a essayé la librairie FFTW (fftw.org)...
j'essai actuellement mais je comprend pas trop trop... si quelqu'un la connait et pouvait en faire un exemple comme celui-ci dessus ca serait genial...

:)
merci
Messages postés
1
Date d'inscription
dimanche 9 janvier 2005
Statut
Membre
Dernière intervention
7 mars 2006

bonjour, je trouve ce code tres interessant, je voudrais savoir si je peux y apporter des modification afin de traiter que des signaux à valeur reel et si oui quelles modification dois-je apporter.
merci
cordialement
Messages postés
1
Date d'inscription
vendredi 18 mars 2005
Statut
Membre
Dernière intervention
25 février 2006

Le pdf que tu a mis concerne la DFT(discrete Fourier Transform) et non la FFT(Fast Fourier Transform)
Voila les liens ou sont expliqués :
la DFT:http://www.spectrumsdi.com/ch8.pdf
la FFT:http://www.spectrumsdi.com/ch12.pdf
Ces pdf expliquent relativement bien comment tout cela marche.
Messages postés
2
Date d'inscription
dimanche 6 février 2005
Statut
Membre
Dernière intervention
6 juin 2005

Pourrais tu mettre ton pdf sur un serveur libre ... Ca serait genial ( mon sujet de stage est proche de la FFT et etant a l'etranger il est difficile pour moi de devenir membre) D'ailleur ils en profitent bien pour se faire de la thune sur notre dos ...C'est + proche de win que de linux...
Afficher les 10 commentaires

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.