比较好的C语言实现FFT程序，加入了时间测试。

#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#define DDC_PI 3.1415926
/////////////////////////////////////////////////////////////////////////

int IsPowerOfTwo ( unsigned x )//=======判断是否为2的整数倍
{
    if ( x < 2 )
        return false;
    if ( x & (x-1) )        // Thanks to 'byang' for this cute trick!
        return false;
    return true;
}

unsigned NumberOfBitsNeeded ( unsigned PowerOfTwo ) 
//===========计算需要运算的点数是2的几次方，返回指数=============
{
    unsigned i;
    if ( PowerOfTwo < 2 )
    {
        fprintf (
            stderr,
            ">>> Error in fftmisc.c: argument %d to NumberOfBitsNeeded is too
small./n",
            PowerOfTwo );
        exit(1);
    }
    for ( i=0; ; i++ )
    {
        if ( PowerOfTwo & (1 << i) )
            return i;
    }
}

 

inline unsigned ReverseBits ( unsigned index, unsigned NumBits )
///===========?????????????
{
    unsigned i, rev;
    for ( i=rev=0; i < NumBits; i++ )
    {
        rev = (rev << 1) | (index & 1);
        index >>= 1;
    }
    return rev;
}

double Index_to_frequency ( unsigned NumSamples, unsigned Index )
{
    if ( Index >= NumSamples )
        return 0.0;
    else if ( Index <= NumSamples/2 )
        return (double)Index / (double)NumSamples;
    return -(double)(NumSamples-Index) / (double)NumSamples;
}

void fft_float (
    unsigned  NumSamples,
    int       InverseTransform,
    double    *RealIn,
    double    *ImagIn,
    double    *RealOut,
    double    *ImagOut )
{
    unsigned NumBits;    /* Number of bits needed to store indices */
    unsigned i, j, k, n;
    unsigned BlockSize, BlockEnd;
    double angle_numerator = 2.0 * DDC_PI;
    double tr, ti;     /* temp real, temp imaginary */
    /*if ( !IsPowerOfTwo(NumSamples) )
    {
        fprintf (
            stderr,
            "Error in fft():  NumSamples=%u is not power of two/n",
            NumSamples );
        exit(1);
    }*/
    //if ( InverseTransform )
        //angle_numerator = -angle_numerator;
    //CHECKPOINTER ( RealIn );//=============???????
    //CHECKPOINTER ( RealOut );
    //CHECKPOINTER ( ImagOut );
    NumBits = NumberOfBitsNeeded ( NumSamples );
    /*
    **   Do simultaneous data copy and bit-reversal ordering into outputs...
    */
    for ( i=0; i < NumSamples; i++ )
    {
        j = ReverseBits ( i, NumBits );
        RealOut[j] = RealIn[i];
        //ImagOut[j] = (ImagIn == NULL) ? 0.0 : ImagIn[i];
  ImagOut[j] = ImagIn[i];
    }
    /*
    **   Do the FFT itself...
    */
    BlockEnd = 1;
    for ( BlockSize = 2; BlockSize <= NumSamples; BlockSize <<= 1 )
    {
        double delta_angle = angle_numerator / (double)BlockSize;
        double sm2 = sin ( -2 * delta_angle );
        double sm1 = sin ( -delta_angle );
        double cm2 = cos ( -2 * delta_angle );
        double cm1 = cos ( -delta_angle );
        double w = 2 * cm1;
        double ar[3], ai[3];
        double temp;
        for ( i=0; i < NumSamples; i += BlockSize )
        {
            ar[2] = cm2;
            ar[1] = cm1;
            ai[2] = sm2;
            ai[1] = sm1;
            for ( j=i, n=0; n < BlockEnd; j++, n++ )
            {
                ar[0] = w*ar[1] - ar[2];
                ar[2] = ar[1];
                ar[1] = ar[0];
                ai[0] = w*ai[1] - ai[2];
                ai[2] = ai[1];
                ai[1] = ai[0];
                k = j + BlockEnd;
                tr = ar[0]*RealOut[k] - ai[0]*ImagOut[k];
                ti = ar[0]*ImagOut[k] + ai[0]*RealOut[k];
                RealOut[k] = RealOut[j] - tr;
                ImagOut[k] = ImagOut[j] - ti;
                RealOut[j] += tr;
                ImagOut[j] += ti;
            }
        }
        BlockEnd = BlockSize;
    }
    /*
    **   Need to normalize if inverse transform...
    */
   /*if ( InverseTransform )
    {
        double denom = (double)NumSamples;
        for ( i=0; i < NumSamples; i++ )
        {
            RealOut[i] /= denom;
            ImagOut[i] /= denom;
        }
    }*/
}

 

 

 

 

/*****************main programe********************/

 

unsigned int   Num=65535;
double RealIn[65535],RealOut[65535];
double ImagIn[65535],ImagOut[65535];

void main()
{
clock_t start,finish;

int i;
start=clock();
for(int j = 0;j<10;j++)
{
for(i=0;i<Num;i++)
{
//s[i].real=sin(pp*i/32);
RealIn[i]= sin(DDC_PI*i/32);
ImagIn[i]=0.0;
//printf("%f ",RealIn[i]);
}

fft_float (
    Num,
    0,
    RealIn,
    ImagIn,
    RealOut,
    ImagOut );
}
finish=clock();
 printf("运行时间是%f/n",difftime(finish,start));//表示的是毫秒即千分之一秒

/*for(i=0;i<16;i++)
{
printf("%.4f",RealOut[i]);
printf("+%.4fj/n",ImagOut[i]);

}*/

}