基于Cooley-Tukey的FFT算法的递归C实现

信号处理课程需要用到FFT变换，于是找资料写了一个C程序，网上的资料不少，但是感觉没什么大用，基本上大部分是MATLAB写的，在Introduction to Algorithm的第三十章有一个算法伪代码图，小心wn这个是每一次都变的而不是一直是2pi/N,我犯了这个错误，一直不大对，后来在cprogramming上的以为朋友指出了这个错误，之后犯的错误是结果一直都是不对称，找来找去不知道错在哪里，第二天才发现关键错在

re_w = re_w*W_RE - im_w*W_IM;
 im_w = im_w*W_RE + re_w*W_IM;

应该替换为

save_re = re_w;save_im = im_w;
re_w = save_re*W_RE(len_x) - save_im*W_IM(len_x);
im_w = save_im*W_RE(len_x) + save_re*W_IM(len_x);

好了，贴上我的代码，大家有什么问题欢迎讨论

#include<stdio.h>
#include<math.h>

#define PI 3.1415926
#define N 16 //需要处理的离散信号的总个数

#define W_RE(n) (cos(-2.0*PI/n))
#define W_IM(n) (sin(-2.0*PI/n))

void recursive_fft(double [],int,double [],double []);

int main(){

    double x[N] = {
    0.5751,
    0.4514,
    0.0439,
    0.0272,
    0.3127,
    0.0129,
    0.3840,
    0.6831,
    0.0928,
    0.0353,
    0.6124,
    0.6085,
    0.0158,
    0.0164,
    0.1901,
    0.5869
    };
    double y_re[N],y_im[N];

    recursive_fft(x,N,y_re,y_im);

    /*输出结果*/
    for(int i = 0;i<N;i++){
        printf("%lf",y_re[i]);
        if(y_im[i])
            printf(" + %lfi",y_im[i]);
        printf("\n");
    }
    return 0;
}

/*x是离散信号值，len_x为其长度
 *y_re和y_im分别是结果的实部和虚部
 * */
void recursive_fft(double x[],int len_x,\
                        double y_re[],double y_im[]){

    double re_w = 1.0,im_w = 0,save_re,save_im;
    double x_even[len_x/2],x_odd[len_x/2];
    double y1_re[len_x/2] ,y1_im[len_x/2];
    double y2_re[len_x/2] ,y2_im[len_x/2];
    int index_even = 0,index_odd = 0;

    for(int i=0;i<len_x/2;i++){
            y1_re[i] = y1_im[i] = y2_re[i] = y2_im[i] = 0;
    }

    if(1 == len_x){
        y_re[0] = x[0];y_im[0] = 0;
        return ;
    }

    /*按奇偶分组*/
    for(int i = 0;i<len_x;i++){
        if(i % 2 == 0)/*偶数*/
            x_even[index_even++] = x[i];
        else
            x_odd[index_odd++] = x[i];
    }
    /*处理奇数*/
    recursive_fft(x_odd,index_odd,y2_re,y2_im);
    /*处理偶数*/
    recursive_fft(x_even,index_even,y1_re,y1_im);

    for(int k = 0;k<len_x/2;k++){
        y_re[k] = y2_re[k]*re_w - im_w*y2_im[k] + y1_re[k];
        y_im[k] = y2_im[k]*re_w + y2_re[k]*im_w + y1_im[k];
        y_re[k+len_x/2] = im_w*y2_im[k] - y2_re[k]*re_w + y1_re[k];
        y_im[k+len_x/2] = -1*y2_im[k]*re_w - y2_re[k]*im_w + y1_im[k];
        save_re = re_w;save_im = im_w;
        re_w = save_re*W_RE(len_x) - save_im*W_IM(len_x);
        im_w = save_im*W_RE(len_x) + save_re*W_IM(len_x);

    }
}