1. #include "uart.h"
    2. #include "fft.h" 
    3. void main()
    4. {
    5.     u8 Freq=0;                //峰值对应频率
    7.     //从y=200*sin(2*pi*40*t)中取得数据点                                         
    8.     s16 code DATA[128]={0,118 ,    -190 ,    190 ,    -118 ,    0 ,    118 ,    -190 ,    
    9.     190 ,-118 ,    0 ,    118 ,    -190 ,    190 ,    -118,     0 ,    
    10.     118 ,-190 ,    190 ,    -118 ,    0 ,    118 ,    -190 ,    190 ,    
    11.     -118,0 ,118 ,    -190 ,    190 ,    -118,     0 ,    118 ,    
    12.     -190 ,    190 ,    -118 ,    0 ,    118 ,    -190 ,    190 ,    -118 ,    
    13.     0 ,    118 ,    -190 ,    190 ,    -118,     0 ,    118 ,    -190 ,    
    14.     190 ,    -118 ,    0 ,    118 ,    -190 ,    190 ,    -118 ,    0 ,    
    15.     118 ,    -190 ,    190 ,    -118 ,    0 ,    118 ,    -190 ,    190 ,
    16.     -118 ,    0 ,    118 ,    -190 ,    190 ,    -118 ,    0 ,    118 ,
    17.     -190 ,    190 ,    -118 ,    0 ,    118 ,    -190 ,    190 ,    -118 ,    
    18.     0 ,    118 ,    -190 ,    190 ,    -118 ,    0 ,    118 ,    -190 ,    
    19.     190 ,    -118 ,    0 ,    118 ,    -190 ,    190 ,    -118 ,    0 ,    
    20.     118 ,    -190 ,    190 ,    -118 ,    0 ,    118 ,    -190 ,    190 ,    
    21.     -118 ,    0 ,    118 ,    -190 ,    190 ,    -118 ,    0 ,    118 ,    
    22.     -190 ,    190 ,    -118 ,    0 ,    118 ,    -190 ,    190 ,    -118,    
    23.     0 ,    118 ,    -190,     190 ,    -118,     0 ,    118 ,    -190 };
    25.     u8 ii,Max_Place=0;               //用Max_Place保存频谱中最大峰值在数组中的位置
    26.     float Out_abs[NN/2];               //用于保存模值，因为频谱图关于Fs/2对称，所以只求一半
    27.     complex y[NN];
    28.     for (ii = 0; ii < NN; ++ii)
    29.     {
    30.         y[ii].real = DATA[ii];
    31.         y[ii].imag = 0;
    32.     }
    33.     fft(NN, y);                      //FFT运算
    34.     c_abs(y, Out_abs, NN/2);      //求模
    35.     Max_Place=Mag(Out_abs);          //求最大模所在位置
    36.     Freq=(Max_Place-1)*Fs/(NN-1);          //得到最大模对应频率
    37.     UartInit();
    38.     printf("Freq: %bd,max_place: %bd\r\n",Freq,Max_Place);
    39.     while(1)
    40.     {
    42.     }
    43. }
    1. #ifndef __FFT_H__  
    2. #define __FFT_H__  
    4. typedef struct complex //复数类型  
    5. {  
    6.   float real;       //实部  
    7.   float imag;       //虚部  
    8. }complex;  
    10. #define PI 3.1415926535  
    11. #define NN 16         //采样点数
    12. #define Fs 100        //采样频率
    14. #define u8 unsigned char
    15. #define s8 signed char
    16. #define u16 unsigned int
    17. #define s16 signed int
    20. ///////////////////////////////////////////In file fft.c  
    22. void c_plus(complex a,complex b,complex *c);//复数加  
    23. void c_mul(complex a,complex b,complex *c) ;//复数乘  
    24. void c_sub(complex a,complex b,complex *c); //复数减法  
    26. void fft(int N,complex f[]);//傅立叶变换 输出也存在数组f中  
    27. void c_abs(complex f[],float out[],int n);//复数数组取模
    28. u16 Mag(float outRes[]);  //求最大值所在位置
    29. ////////////////////////////////////////////  
    32. #endif
    1. #include "math.h" 
    2. #include "fft.h" 
    6. //精度0.0001弧度  
    7. void c_abs(complex f[],float out[],int n)  
    8. {  
    9.   int i = 0;  
    10.   float t;  
    11.   for(i=0;i<n;i++)  
    12.   {  
    13.     t = f[i].real * f[i].real + f[i].imag * f[i].imag;  
    14.     out[i] = sqrt(t);  
    15.   }   
    16. }  
    19. void c_plus(complex a,complex b,complex *c)  
    20. {  
    21.   c->real = a.real + b.real;  
    22.   c->imag = a.imag + b.imag;  
    23. }  
    25. void c_sub(complex a,complex b,complex *c)  
    26. {  
    27.   c->real = a.real - b.real;  
    28.   c->imag = a.imag - b.imag;   
    29. }  
    31. void c_mul(complex a,complex b,complex *c)  
    32. {  
    33.   c->real = a.real * b.real - a.imag * b.imag;  
    34.   c->imag = a.real * b.imag + a.imag * b.real;     
    35. }  
    39. void Wn_i(int n,int i,complex *Wn)  
    40. {  
    41.   Wn->real = cos(2*PI*i/n);   
    42.   Wn->imag = -sin(2*PI*i/n);  
    43. }  
    45. //傅里叶变化  
    46. void fft(int N,complex f[])  
    47. {  
    48.   complex t,wn;//中间变量  
    49.   int i,j,k,m,n,l,r,M;  
    50.   int la,lb,lc;  
    51.   /*----计算分解的级数M=log2(N)----*/  
    52.   for(i=N,M=1;(i=i/2)!=1;M++);   
    53.   /*----按照倒位序重新排列原信号----*/  
    54.   for(i=1,j=N/2;i<=N-2;i++)  
    55.   {  
    56.     if(i<j)  
    57.     {  
    58.       t=f[j];  
    59.       f[j]=f[i];  
    60.       f[i]=t;  
    61.     }  
    62.     k=N/2;  
    63.     while(k<=j)  
    64.     {  
    65.       j=j-k;  
    66.       k=k/2;  
    67.     }  
    68.     j=j+k;  
    69.   }  
    71.   /*----FFT算法----*/  
    72.   for(m=1;m<=M;m++)  
    73.   {  
    74.     la=pow(2,m); //la=2^m代表第m级每个分组所含节点数       
    75.     lb=la/2;    //lb代表第m级每个分组所含碟形单元数  
    76.                  //同时它也表示每个碟形单元上下节点之间的距离  
    77.     /*----碟形运算----*/  
    78.     for(l=1;l<=lb;l++)  
    79.     {  
    80.       r=(l-1)*pow(2,M-m);     
    81.       for(n=l-1;n<N-1;n=n+la) //遍历每个分组，分组总数为N/la  
    82.       {  
    83.         lc=n+lb;  //n,lc分别代表一个碟形单元的上、下节点编号       
    84.         Wn_i(N,r,&wn);//wn=Wnr  
    85.         c_mul(f[lc],wn,&t);//t = f[lc] * wn复数运算  
    86.         c_sub(f[n],t,&(f[lc]));//f[lc] = f[n] - f[lc] * Wnr  
    87.         c_plus(f[n],t,&(f[n]));//f[n] = f[n] + f[lc] * Wnr  
    88.       }  
    89.     }  
    90.   }  
    91. } 
    93. u16 Mag(float outRes[])
    94. {
    95.     u8 i=1;
    96.     u8 max,place;
    97.     max=outRes[i];
    98.     place=i;
    99.     for(i=2;i<NN/2;i++)
    100.     {
    101.         if(outRes[i]>max)
    102.         {
    103.             max=outRes[i];
    104.             place=i;
    105.         }        
    106.     }
    107.     return place;    
    108. }