Fast fourier transform

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

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)