麻烦大家看一下这个FFT程序
看了算法导论第30章后写的，按理说对于数据
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
求FFT应该得到
4.6485
0.0176 + 0.3122i
1.1116 + 0.0427i
1.6777 - 0.1353i
-0.2339 + 1.3898i
0.3651 - 1.2589i
-0.4325 + 0.2072i
-0.1312 + 0.3763i
-0.1950
-0.1312 - 0.3763i
-0.4325 - 0.2072i
0.3651 + 1.2589i
-0.2339 - 1.3898i
1.6777 + 0.1353i
1.1116 - 0.0427i
0.0176 - 0.3122i
可是我的结果不大对，大家帮忙看看，谢谢了
#include<stdio.h>#include<stdlib.h>#include<math.h>#define PI 3.1415926#define N 16 //需要处理的离散信号的总个数#define W_RE (cos(-2.0*PI/N))#define W_IM (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; 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,y1_re,y1_im); /*处理偶数*/ recursive_fft(x_even,index_even,y2_re,y2_im); for(int k = 0;k<len_x/2;k++){ y_re[k] = y1_re[k]*re_w - im_w*y1_im[k] + y2_re[k]; y_im[k] = y1_im[k]*re_w + y1_re[k]*im_w + y2_im[k]; y_re[k+len_x/2] = im_w*y1_im[k] - y1_re[k]*re_w + y2_re[k]; y_im[k+len_x/2] = -1*y1_im[k]*re_w - y1_re[k]*im_w + y2_im[k]; re_w = re_w*W_RE - im_w*W_IM; im_w = im_w*W_RE + re_w*W_IM; }}