这个你参考下：

1. #include <reg52.h>
   
2. #include <intrins.h>
   
3. 
   
4. 
   
5. #include<math.h>
   
6. 
   
7. #define PI 3.14159               //定义圆周率值
   
8. #define FFT_N 32                                                  //定义福利叶变换的点数
   
9. 
   
10. struct compx {float real,imag;};                                    //定义一个复数结构
    
11. struct compx s[FFT_N];                                              //FFT输入和输出：从S[1]开始存放，根据大小自己定义
    
12. 
    
13. 
    
14. 
    
15. struct compx EE(struct compx a,struct compx b)     
    
16. {
    
17. struct compx c;
    
18. c.real=a.real*b.real-a.imag*b.imag;
    
19. c.imag=a.real*b.imag+a.imag*b.real;
    
20. return(c);
    
21. }
    
22. 
    
23. 
    
24. void FFT(struct compx *xin)
    
25. {
    
26.   int f,m,nv2,nm1,i,k,l,j=0;
    
27.   struct compx u,w,t;
    
28.   
    
29.    nv2=FFT_N/2;                  //变址运算，即把自然顺序变成倒位序，采用雷德算法
    
30.    nm1=FFT_N-1;
    
31.    for(i=0;i<nm1;i++)      
    
32.    {
    
33.     if(i<j)                    //如果i<j,即进行变址
    
34.      {
    
35.       t=xin[j];         
    
36.       xin[j]=xin[i];
    
37.       xin[i]=t;
    
38.      }
    
39.     k=nv2;                    //求j的下一个倒位序
    
40.     while(k<=j)               //如果k<=j,表示j的最高位为1  
    
41.      {         
    
42.       j=j-k;                 //把最高位变成0
    
43.       k=k/2;                 //k/2，比较次高位，依次类推，逐个比较，直到某个位为0
    
44.      }
    
45.    j=j+k;                   //把0改为1
    
46.   }
    
47.                         
    
48.   {
    
49.    int le,lei,ip;                            //FFT运算核，使用蝶形运算完成FFT运算
    
50.     f=FFT_N;
    
51.    for(l=1;(f=f/2)!=1;l++)                  //计算l的值，即计算蝶形级数
    
52.            ;
    
53.   for(m=1;m<=l;m++)                         // 控制蝶形结级数
    
54.    {                                        //m表示第m级蝶形，l为蝶形级总数l=log（2）N
    
55.     le=2<<(m-1);                            //le蝶形结距离，即第m级蝶形的蝶形结相距le点
    
56.     lei=le/2;                               //同一蝶形结中参加运算的两点的距离
    
57.     u.real=1.0;                             //u为蝶形结运算系数，初始值为1
    
58.     u.imag=0.0;
    
59.     w.real=cos(PI/lei);                     //w为系数商，即当前系数与前一个系数的商
    
60.     w.imag=-sin(PI/lei);
    
61.     for(j=0;j<=lei-1;j++)                   //控制计算不同种蝶形结，即计算系数不同的蝶形结
    
62.      {
    
63.       for(i=j;i<=FFT_N-1;i=i+le)            //控制同一蝶形结运算，即计算系数相同蝶形结
    
64.        {
    
65.         ip=i+lei;                           //i，ip分别表示参加蝶形运算的两个节点
    
66.         t=EE(xin[ip],u);                    //蝶形运算，详见公式
    
67.         xin[ip].real=xin[i].real-t.real;
    
68.         xin[ip].imag=xin[i].imag-t.imag;
    
69.         xin[i].real=xin[i].real+t.real;
    
70.         xin[i].imag=xin[i].imag+t.imag;
    
71.        }
    
72.       u=EE(u,w);                           //改变系数，进行下一个蝶形运算
    
73.      }
    
74.    }
    
75.   }
    
76. 
    
77. }
    
78. 
    
79. 
    
80. void main()  
    
81. {
    
82.   unsigned char i;
    
83.   for(i=0;i<FFT_N;i++)                           //给结构体赋值
    
84.   {
    
85.      s[i].real=sin(2*PI*i/FFT_N); //实部为正弦波FFT_N点采样，赋值为1
    
86.      s[i].imag=0;                                //虚部为0
    
87.   }
    
88. 
    
89.   FFT(s);                                        //进行快速福利叶变换
    
90. 
    
91.   for(i=0;i<FFT_N;i++)                           //求变换后结果的模值，存入复数的实部部分
    
92.   s[i].real=sqrt(s[i].real*s[i].real+s[i].imag*s[i].imag);
    
93. 
    
94.    while(1);
    
95. }
    

复制代码