//////////////////////////////////////// //////////////////
// Internal definitions
# Define Pi (double) 3.14159265359
/* Complex Number */
Typedef struct
{
Double re;
Double im;
} Complex;
/* Complex add */
Complex add (complex C1, complex C2)
{
Complex C;
C. Re = c1.re + c2.re;
C. Im = c1.im + c2.im;
Return C;
}
/* Complex substract */
Complex sub (complex C1, complex C2)
{
Complex C;
C. Re = c1.re-c2.re;
C. Im = c1.im-c2.im;
Return C;
}
/* Complex multiple */
Complex MUL (complex C1, complex C2)
{
Complex C;
C. Re = c1.re * c2.re-c1.im * c2.im;
C. Im = c1.re * c2.im + c2.re * c1.im;
Return C;
}
//////////////////////////////////////// //////////////////
This is the FFT /************************************* ***************
FFT ()
Parameters:
TD is the time domain Value
FD is a frequency range value.
Power is a power of 2
Return Value:
None
Note:
This function implements fast Fourier transformation.
**************************************** ************/
Void FFT (complex * TD, complex * FD, int power)
{
Int count;
Int I, J, K, bfsize, P;
Double angle;
Complex * w, * X1, * X2, * X;
/* Calculate the points of Fourier transformation */
Count = 1 <power;
/* Allocate memory required for calculation */
W = (complex *) malloc (sizeof (complex) * count/2 );
X1 = (complex *) malloc (sizeof (complex) * count );
X2 = (complex *) malloc (sizeof (complex) * count );
/* Calculate the weighting coefficient */
For (I = 0; I <count/2; I ++)
{
Angle =-I * pI * 2/count;
W [I]. Re = cos (angle );
W [I]. Im = sin (angle );
}
/* Write the time domain point into the memory */
Memcpy (x1, TD, sizeof (complex) * count );
/* Butterfly operation */
For (k = 0; k <power; k ++)
{
For (j = 0; j <1 <K; j ++)
{
Bfsize = 1 <(power-k );
For (I = 0; I <bfsize/2; I ++)
{
P = J * bfsize;
X2 [I + P] = add (x1 [I + P], X1 [I + P + bfsize/2]);
X2 [I + P + bfsize/2] = MUL (sub (x1 [I + P], X1 [I + P + bfsize/2]), W [I * (1 <k)]);
}
}
X = x1;
X1 = x2;
X2 = X;
}
/* Re-Sort */
For (j = 0; j <count; j ++)
{
P = 0;
For (I = 0; I <power; I ++)
{
If (J & (1 <I) P + = 1 <(power-i-1 );
}
FD [J] = x1 [p];
}
/* Release the memory */
Free (w );
Free (X1 );
Free (X2 );
}
/*************************************** *************
IFFT ()
Parameters:
FD is a frequency range value.
TD is the time domain Value
Power is a power of 2
Return Value:
None
Note:
This function uses fast Fourier transformation to implement Fourier inverse transformation.
**************************************** ************/
Void IFFT (complex * FD, complex * TD, int power)
{
Int I, count;
Complex * X;
/* Calculate the number of points of Fourier inverse transformation */
Count = 1 <power;
/* Allocate memory required for calculation */
X = (complex *) malloc (sizeof (complex) * count );
/* Write the frequency point into the memory */
Memcpy (x, FD, sizeof (complex) * count );
/* Calculate the frequency-domain point concatenation */
For (I = 0; I <count; I ++)
X [I]. Im =-X [I]. Im;
/* Call FFT */
FFT (x, TD, power );
/* Calculate the time-domain point's concatenation */
For (I = 0; I <count; I ++)
{
TD [I]. Re/= count;
TD [I]. Im =-Td [I]. Im/count;
}
/* Release the memory */
Free (X );
}