Algorithme de fft en c

Soyez le premier à donner votre avis sur cette source.

Vue 31 975 fois - Téléchargée 1 417 fois

Description

Les fonctionnalites:
DIF Radix2, DIF Radix4
DIT Radix2, DIT Radix4
FFT complex, IFFT complex
FFT real, IFFT real

Le mode d'emploi:
lguo@lguo-enst:~$ unzip -r Algo.zip
lguo@lguo-enst:~$ cd Algo/DIF
lguo@lguo-enst:~/Algo/DIF$ gcc -o FFT_DIF_R2 FFT_DIF_R2.c -lm
lguo@lguo-enst:~/Algo/DIF$ ./FFT_DIF_R2

lguo@lguo-enst:~$ cd Algo/FFT_Test_Linux
lguo@lguo-enst:~/Algo/FFT_Test_Linux$ make
lguo@lguo-enst:~/Algo/FFT_Test_Linux$ ./FFT

Source / Exemple :


/***********************************
// DIT radix-4  FFT  complex 
//
// 1. Maximium points are 1024.
//
// 2. The last stage is 2 DFT ( for 8, 32, 128, 512...)
// or all stages are 4 DFT ( for 4, 16, 64, 256, 1024 ...).
//
// 3. The functions special for FFT real .
//
// 24 juillet 2007
// purcharse*gmail.com

                                                                        • /
#include "FFT.h" static complex multicomplex( complex b1, complex b2) /* multiplication of complex */ { complex b3; b3.real=b1.real*b2.real-b1.imag*b2.imag; b3.imag=b1.real*b2.imag+b1.imag*b2.real; return(b3); } static int mylog2(int N) /* Max(N) = 4098 */ { int k=0; if (N>>12) { k+=12; N>>=12; } if (N>>8) { k+=8; N>>=8; } if (N>>4) { k+=4; N>>=4; } if (N>>2) { k+=2; N>>=2; } if (N>>1) { k+=1; N>>=1; } return k ; } static void BitReverse(complex *xin, int N) { int LH, i, j, k; struct complex tmp; LH=N/2; j = N/2; for( i = 1; i <= (N -2); i++) { if(i < j) { tmp = xin[j]; xin[j] = xin[i]; xin[i] = tmp; } k = LH; while(j >= k) { j = j-k; k = k/2; } j = j + k; } } static void DFT_2(complex *b1, complex *b2) { struct complex tmp; tmp = *b1; (*b1).real = (*b1).real + (*b2).real; (*b1).imag = (*b1).imag + (*b2).imag; (*b2).real = tmp.real - (*b2).real; (*b2).imag = tmp.imag - (*b2).imag; } static void DFT_4(complex* b0, complex* b1, complex* b2, complex* b3) { /*variables locales*/ struct complex temp[4]; /*calcul x1*/ temp[0].real=(*b0).real+(*b1).real; temp[0].imag=(*b0).imag+(*b1).imag; /*calcul x2*/ temp[1].real=(*b0).real-(*b1).real; temp[1].imag=(*b0).imag-(*b1).imag; /*calcul x3*/ temp[2].real=(*b2).real+(*b3).real; temp[2].imag=(*b2).imag+(*b3).imag; /*calcul x4 + multiplication with -j*/ temp[3].imag=(*b3).real-(*b2).real; temp[3].real=(*b2).imag-(*b3).imag; /*the last stage*/ (*b0).real=temp[0].real+temp[2].real; (*b0).imag=temp[0].imag+temp[2].imag; (*b1).real=temp[1].real+temp[3].real; (*b1).imag=temp[1].imag+temp[3].imag; (*b2).real=temp[0].real-temp[2].real; (*b2).imag=temp[0].imag-temp[2].imag; (*b3).real=temp[1].real-temp[3].real; (*b3).imag=temp[1].imag-temp[3].imag; } static void FFT_R4(complex *xin, int N, int m) { int i, L, j; double ps1, ps2, ps3; int le,B; struct complex w[4]; for( L = 1; L <= m; L++) { le = pow(4 ,L); B = le/4; /*the distance of buttefly*/ for(j = 0; j <= B-1 ; j++) { // ps0 = (TWICEPI/N) * 0 * j; // w[0].real = cos(ps0); // w[0].imag = -sin(ps0); ps1 = ((TWICEPI)/le)*2*j; w[1].real = cos(ps1); w[1].imag = -sin(ps1); ps2 = (TWICEPI/le)*j; w[2].real = cos(ps2); w[2].imag = -sin(ps2); ps3 = (TWICEPI/le)*3*j; w[3].real = cos(ps3); w[3].imag = -sin(ps3); for(i = j; i <= N-1; i = i + le) /* controle those same butteflies*/ { /* multiple with W */ // xin[i] = multicomplex(xin[i], w[0]); xin[i + B] = multicomplex(xin[i + B], w[1]); xin[i + 2*B] = multicomplex(xin[i + 2*B], w[2]); xin[i + 3*B] = multicomplex(xin[i + 3*B], w[3]); /* DFT-4 */ DFT_4(xin + i, xin + i + B, xin + i + 2*B, xin + i + 3*B); } } /* printf("*****N°%d **********\n", L); for(i=0;i<N;i++) { printf("%.8f\t\t",xin[i].real); printf("%.8f\n",xin[i].imag); }
  • /
} }//end of FFT_R4 static void FFT_L2(complex *xin, int N) { /* For the last stage 2 DFT*/ int j, B; double p, ps ; struct complex w; B = N/2; for(j = 0; j <= B - 1; j++) { ps = (TWICEPI/N)*j; w.real = cos(ps); w.imag = -sin(ps); /* multiple avec W */ xin[j+ B] = multicomplex(xin[j + B], w); DFT_2(xin + j ,xin + j + B); } }//end of FFT_L2 /**************** FFT complex ***************/ void RunFFT(complex *xin, int N) { int m, i; BitReverse(xin, N); m = mylog2(N); if( (m%2) == 0 ) { /*All the stages are radix 4*/ FFT_R4(xin, N, m/2); } else { /*the last stage is radix 2*/ FFT_R4(xin, N, m/2); FFT_L2(xin, N); } } void RunIFFT(complex *xin,int N) { /* inverse FFT */ int i; for(i=0; i < N + 1 ; i++) { xin[i].imag = -xin[i].imag; } RunFFT(xin,N); for(i = 0; i < N + 1 ; i++) { xin[i].real = xin[i].real/N; xin[i].imag = -xin[i].imag/N; } } /************** FFT real ****************/ void RunFFTR(complex *xin, int N) { int i; double ps; complex *Realin; complex Realtmp1; complex Realtmp2; complex w; Realin = (complex *)malloc((N/2)*sizeof(complex)); /*** X(n)= A(2n) +j*A(2n+1) ***/ for(i = 0 ; i < N/2 ; i++) { Realin[i].real = xin[2*i].real; Realin[i].imag = xin[2*i + 1].real; } RunFFT(Realin, N/2); for(i = 0; i < N/2 ; i++) { /***** factor w *****/ ps = (TWICEPI/N)*i; w.real = cos(ps); w.imag = -sin(ps); /***** conjugue *****/ Realtmp1.real = Realin[(N/2) - i].real; Realtmp1.imag = -Realin[(N/2) - i].imag; if(i == 0) { Realtmp1.real = Realin[0].real; Realtmp1.imag = -Realin[0].imag; } /***** part 2 *****/ Realtmp2.real = Realin[i].imag - Realtmp1.imag; Realtmp2.imag = Realtmp1.real - Realin[i].real; if(i > 0) Realtmp2 = multicomplex(Realtmp2, w); /**** part 1 ****/ Realtmp1.real = Realin[i].real + Realtmp1.real; Realtmp1.imag = Realin[i].imag + Realtmp1.imag; xin[i].real = (Realtmp1.real + Realtmp2.real)/2; xin[i].imag = (Realtmp1.imag + Realtmp2.imag)/2; xin[N/2 + i].real = (Realtmp1.real - Realtmp2.real)/2; xin[N/2 + i].imag = (Realtmp1.imag - Realtmp2.imag)/2; } } void RunIFFTR(complex *xin, int N) { int i; double ps; complex *Realin; complex Realtmp1; complex Realtmp2; complex w; Realin = (complex *)malloc((N/2)*sizeof(complex)); /*** X(n)= A(2n) +j*A(2n+1) ***/ for(i = 0 ; i < N/2 ; i++) { Realin[i].real = xin[2*i].real; Realin[i].imag = xin[2*i + 1].real; } RunIFFT(Realin, N/2); for(i = 0; i < N/2 ; i++) { /***** factor w *****/ ps = (TWICEPI/N)*i; w.real = cos(ps); w.imag = -sin(ps); /***** conjugue *****/ Realtmp1.real = Realin[(N/2) - i].real; Realtmp1.imag = -Realin[(N/2) - i].imag; if(i == 0) { Realtmp1.real = Realin[0].real; Realtmp1.imag = -Realin[0].imag; } /***** part 2 *****/ Realtmp2.real = Realin[i].imag - Realtmp1.imag; Realtmp2.imag = Realtmp1.real - Realin[i].real; if(i > 0) Realtmp2 = multicomplex(Realtmp2, w); /**** part 1 ****/ Realtmp1.real = Realin[i].real + Realtmp1.real; Realtmp1.imag = Realin[i].imag + Realtmp1.imag; xin[i].real = (Realtmp1.real + Realtmp2.real)/2; xin[i].imag = (Realtmp1.imag + Realtmp2.imag)/2; xin[N/2 + i].real = (Realtmp1.real - Realtmp2.real)/2; xin[N/2 + i].imag = (Realtmp1.imag - Realtmp2.imag)/2; } }

Codes Sources

A voir également

Ajouter un commentaire

Commentaires

cs_JCDjcd
Messages postés
1138
Date d'inscription
mardi 10 juin 2003
Statut
Membre
Dernière intervention
25 janvier 2009
2
cs_JCDjcd
Messages postés
1138
Date d'inscription
mardi 10 juin 2003
Statut
Membre
Dernière intervention
25 janvier 2009
2
Bon alors en fait du fait un FFT avec le bit reverse.
Cet algorithme est plus difficile a implementer.
Ton code n'est pas lisible, on sait pas ou l'on va.
Meme si je connais le principe, je ne comprends pas.
On peut faire une FFT directement.
A toute fin utile, voici ma version :
# //-------------------------------------------------
# // on decompose le polynome P = A0 + B0 avec les coefficients
# // de degres pairs (en A0) et impairs (en B0)
# // alors on a :
# // P(w) = A0(w) + B0(w)
# // = A(w^2) + w.B(w^2)
# // or w^2 est une racine nieme, donc on se ramene a un cas
# // de taille deux fois moins grand
# static void polyFFT(
# P_COMPLEX res,
# P_COMPLEX tabCoeff,
# int stepi,
# int n,
# P_COMPLEX w
# )
# {
#
# if(1 == n)
# {
# // P(z)=constante
# res->re = tabCoeff->re;
# res->im = tabCoeff->im;
# }
# else
# {
# int k;
# COMPLEX w2,wPowk;
#
#
# // strategie diviser pour regner
#
# // calcul du carre de <w>
# MulComplex(&w2,w,w);
# // appel recursif
# polyFFT(res ,tabCoeff ,2*stepi,n/2,&w2); // coefficients pairs
# polyFFT(res + n/2,tabCoeff + stepi ,2*stepi,n/2,&w2); // coefficients impairs
#
# // on calcule la FFT a partir de ce que l'on vient de calculer
# wPowk.re = 1.;
# wPowk.im = 0.;
# for(k=0;k<n/2;k++)
# {
# COMPLEX A,B;
#
# A = res[k];
# B = res[k + n/2];
# MulComplex(&B,&B,&wPowk);
#
# AddComplex(res + k ,&A,&B);
# SubComplex(res + k + n/2,&A,&B);
#
# MulComplex(&wPowk,&wPowk,w);
# }
# }
# } // polyFFT()

bien plus courte (normalement bien indentee)
cs_Patrice99
Messages postés
1222
Date d'inscription
jeudi 23 août 2001
Statut
Membre
Dernière intervention
9 septembre 2018

Voir aussi en VB 2005 Express :
VBWaveComp : Le comparateur de spectre audio en VB .Net
Vers un "benchmark" de la compression audio
www.vbfrance.com/code.aspx?ID=5319
ghuysmans99
Messages postés
2501
Date d'inscription
jeudi 14 juillet 2005
Statut
Contributeur
Dernière intervention
5 juin 2016
1
pourquoi tu ne fais pas toute la compil dans un makefile ?

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.