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