C语言写的FFT代码

#include <stdio.h>

#include <math.h>

#include <stdlib.h>



#define   N     8     //64  输入样本总数

#define    M     3   //DFT运算层数     //2^m=N  

#define    PI    3.1415926



float   twiddle[N/2] = {1.0, 0.707, 0.0, -0.707};

float   x_r[N] = {1, 1, 1, 1, 0, 0, 0, 0};  //输入数据，此处设为8个

float   x_i[N];                         //N=8





/**

* 初始化输出虚部

*/

static void fft_init( void )

{

    int i;

    for(i=0; i<N; i++)   x_i[i] = 0.0;

}



/**

* 反转算法.将时域信号重新排序.

* 这个算法有改进的空间

*/

static void bitrev( void )

{

    int    p=1, q, i;

    int    bit_rev[ N ];  

    float   xx_r[ N ];   

   

    bit_rev[ 0 ] = 0;

    while( p < N )

    {

       for(q=0; q<p; q++)  

       {

           bit_rev[ q ]     = bit_rev[ q ] * 2;

           bit_rev[ q + p ] = bit_rev[ q ] + 1;

       }

       p *= 2;

    }

   

    for(i=0; i<N; i++)   xx_r[ i ] = x_r[ i ];   

   

    for(i=0; i<N; i++)   x_r[i] = xx_r[ bit_rev[i] ];

}



void fft( void )

{   fp = fopen("log2.txt", "a+");//此处

    int     cur_layer, gr_num, i, k, p;        //cur_layer代表正要计算的当前层，gr_num代表当前层的颗粒数

    float   tmp_real, tmp_imag, temp;   // 临时变量, 记录实部

    float   tw1, tw2;// 旋转因子,tw1为旋转因子的实部cos部分, tw2为旋转因子的虚部sin部分.

      

    int    step;      // 步进

    int    sample_num;   // 颗粒的样本总数(各层不同, 因为各层颗粒的输入不同)

   

    /* 对层循环 */

    for(cur_layer=1; cur_layer<=M; cur_layer++)

    {      

       /* 求当前层拥有多少个颗粒(gr_num) */

       gr_num = 1;

       i = M - cur_layer;

       while(i > 0)

       {

           i--;

           gr_num *= 2;

       }

      

       /* 每个颗粒的输入样本数N' */

       sample_num    = (int)pow(2, cur_layer);

       /* 步进. 步进是N'/2 */

       step       = sample_num/2;

      

       /*  */

       k = 0;

      

       /* 对颗粒进行循环 */

       for(i=0; i<gr_num; i++)

       {

           /*

            * 对样本点进行循环, 注意上限和步进

            */

           for(p=0; p<sample_num/2; p++)

           {   

              // 旋转因子, 需要优化...   

              tw1 = cos(2*PI*p/pow(2, cur_layer));

              tw2 = -sin(2*PI*p/pow(2, cur_layer));

              

              tmp_real = x_r[k+p];

              tmp_imag = x_i[k+p];

              temp = x_r[k+p+step];

              

              /* 蝶形算法 */

              x_r[k+p]   = tmp_real + ( tw1*x_r[k+p+step] - tw2*x_i[k+p+step] );

              x_i[k+p]   = tmp_imag + ( tw2*x_r[k+p+step] + tw1*x_i[k+p+step] );

              /* X[k] = A(k)+WB(k)

               * X[k+N/2] = A(k)-WB(k) 的性质可以优化这里*/

              /*旋转因子, 需要优化...

              tw1 = cos(2*PI*(p+step)/pow(2, cur_layer));

              tw2 = -sin(2*PI*(p+step)/pow(2, cur_layer));

              x_r[k+p+step] = tmp_real + ( tw1*temp - tw2*x_i[k+p+step] );

              x_i[k+p+step] = tmp_imag + ( tw2*temp + tw1*x_i[k+p+step] );*/

        x_r[k+p+step]   = tmp_real - ( tw1* temp - tw2*x_i[k+p+step] );

              x_i[k+p+step]   = tmp_imag - ( tw2* temp + tw1*x_i[k+p+step] );

              

              printf("k=%d, x_r[k]=%f, x_i[k]=%f\n", k+p, x_r[k+p], x_i[k+p]);

              printf("k=%d, x_r[k]=%f, x_i[k]=%f\n", k+p+step, x_r[k+p+step], x_i[k+p+step]);

           }

           /* 开跳!:) */

           k += 2*step;

       }   

    }

}



void display( void )

{

    printf("\n\n");

    int   i;  

    for(i=0; i<N; i++)

       printf("%f\t%f\n", x_r[i], x_i[i]);

}



int main( void )

{

    fft_init( );                //初始化

    bitrev( );                //将输入直接按FFT计算要求排序，如8点FFT计算，排序为x[0]、x[4]、x[2]、x[6]、x[1]、x[5]、x[3]、x[7]

    fft( );                        //进行FFT计算

    display( );                //显示计算结果

   

    system( "pause" );

    return 1;

}