现在的位置: 首页 > 综合 > 正文

用于ARM上的FFT与IFFT源代码（C语言，不依赖特定平台）

2018年06月07日 ⁄ 综合 ⁄ 共 38869字 ⁄ 字号 小 中 大 ⁄ 评论关闭

http://blog.csdn.net/piaojun_pj/article/details/5911914

代码在2011年全国电子大赛结束后（2011年9月3日）发布，多个版本，注释详细。

[cpp] view
plaincopy

  1. /******************************************************************************* 
  2. ** 程序名称：快速傅里叶变换(FFT)  
  3. ** 程序描述：本程序实现快速傅里叶变换  
  4. ** 程序作者：宋元瑞  
  5. ** 最后修改：2011年4月5日  
  6. *******************************************************************************/  
  7. #include <stdio.h>  
  8. #include <math.h>  
  9.   
  10. #define PI 3.141592653589   //圆周率，12位小数   
  11. #define N 8                 //傅里叶变换的点数   
  12. #define M 3                 //蝶形运算的级数，N = 2^M   
  13. typedef double ElemType;    //原始数据序列的数据类型,可以在这里设置  
  14.   
  15. typedef struct              //定义复数结构体   
  16. {  
  17.     ElemType real,imag;  
  18. }complex;  
  19.   
  20. complex data[N];            //定义存储单元，原始数据与负数结果均使用之   
  21. ElemType result[N];         //存储FFT后复数结果的模  
  22.   
  23. //变址   
  24. void ChangeSeat(complex *DataInput)  
  25. {  
  26.     int nextValue,nextM,i,k,j=0;  
  27.     complex temp;  
  28.       
  29.     nextValue=N/2;                  //变址运算，即把自然顺序变成倒位序，采用雷德算法  
  30.     nextM=N-1;  
  31.     for (i=0;i<nextM;i++)  
  32.     {  
  33.         if (i<j)                 //如果i<j,即进行变址  
  34.         {  
  35.             temp=DataInput[j];  
  36.             DataInput[j]=DataInput[i];  
  37.             DataInput[i]=temp;  
  38.         }  
  39.         k=nextValue;                //求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. //变址  
  50. void ChangeSeat(complex *DataInput) 
  51. { 
  52.     complex Temp[N]; 
  53.     int i,n,New_seat; 
  54.     for(i=0; i<N; i++) 
  55.     { 
  56.         Temp[i].real = DataInput[i].real; 
  57.         Temp[i].imag = DataInput[i].imag; 
  58.     } 
  59.     for(i=0; i<N; i++) 
  60.     { 
  61.         New_seat = 0; 
  62.         for(n=0;n<M;n++) 
  63.         { 
  64.             New_seat = New_seat | (((i>>n) & 0x01) << (M-n-1)); 
  65.         } 
  66.         DataInput[New_seat].real = Temp[i].real; 
  67.         DataInput[New_seat].imag = Temp[i].imag; 
  68.     } 
  69. } 
  70. */  
  71. //复数乘法   
  72. complex XX_complex(complex a, complex b)  
  73. {  
  74.     complex temp;  
  75.       
  76.     temp.real = a.real * b.real-a.imag*b.imag;  
  77.     temp.imag = b.imag*a.real + a.imag*b.real;  
  78.       
  79.     return temp;  
  80. }  
  81.   
  82. //FFT  
  83. void FFT(void)  
  84. {  
  85.     int L=0,B=0,J=0,K=0;  
  86.     int step=0;  
  87.     ElemType P=0,T=0;  
  88.     complex W,Temp_XX;  
  89.     //ElemType TempResult[N];  
  90.       
  91.     ChangeSeat(data);  
  92.     for(L=1; L<=M; L++)  
  93.     {  
  94.         B = 1<<(L-1);//B=2^(L-1)  
  95.         for(J=0; J<=B-1; J++)  
  96.         {  
  97.             P = (1<<(M-L))*J;//P=2^(M-L) *J   
  98.             step = 1<<L;//2^L  
  99.             for(K=J; K<=N-1; K=K+step)  
  100.             {  
  101.                 W.real =  cos(2*PI*P/N);  
  102.                 W.imag = -sin(2*PI*P/N);  
  103.                   
  104.                 Temp_XX = XX_complex(data[K+B],W);  
  105.                 data[K+B].real = data[K].real - Temp_XX.real;  
  106.                 data[K+B].imag = data[K].imag - Temp_XX.imag;  
  107.                   
  108.                 data[K].real = data[K].real + Temp_XX.real;  
  109.                 data[K].imag = data[K].imag + Temp_XX.imag;  
  110.             }  
  111.         }  
  112.     }  
  113. }  
  114. void IFFT(void)  
  115. {  
  116.     int L=0,B=0,J=0,K=0;  
  117.     int step=0;  
  118.     ElemType P=0,T=0;  
  119.     complex W,Temp_XX;  
  120.     //ElemType TempResult[N];  
  121.       
  122.     ChangeSeat(data);  
  123.     for(L=1; L<=M; L++)  
  124.     {  
  125.         B = 1<<(L-1);//B=2^(L-1)  
  126.         for(J=0; J<=B-1; J++)  
  127.         {  
  128.             P = (1<<(M-L))*J;//P=2^(M-L) *J   
  129.             step = 1<<L;//2^L  
  130.             for(K=J; K<=N-1; K=K+step)  
  131.             {  
  132.                 W.real =  cos(2*PI*P/N);  
  133.                 W.imag =  sin(2*PI*P/N);//逆运算，这里跟FFT符号相反   
  134.                   
  135.                 Temp_XX = XX_complex(data[K+B],W);  
  136.                 data[K+B].real = data[K].real - Temp_XX.real;  
  137.                 data[K+B].imag = data[K].imag - Temp_XX.imag;  
  138.                   
  139.                 data[K].real = data[K].real + Temp_XX.real;  
  140.                 data[K].imag = data[K].imag + Temp_XX.imag;  
  141.             }  
  142.         }  
  143.     }  
  144. }  
  145. int main(int argc, char *argv[])  
  146. {  
  147.     int i = 0;  
  148.     for(i=0; i<N; i++)//制造输入序列   
  149.     {  
  150.         data[i].real = sin(2*PI*i/N);  
  151.         printf("%lf ",data[i]);  
  152.     }  
  153.     printf("\n\n");  
  154.       
  155.       
  156.     FFT();//进行FFT计算   
  157.     printf("\n\n");  
  158.     for(i=0; i<N; i++)  
  159.         {printf("%lf ",sqrt(data[i].real*data[i].real+data[i].imag*data[i].imag));}  
  160.       
  161.     IFFT();//进行FFT计算   
  162.     printf("\n\n");  
  163.     for(i=0; i<N; i++)  
  164.         {printf("%lf ",data[i].real/N);}  
  165.     printf("\n");  
  166.     /*for(i=0; i<N; i++) 
  167.         {printf("%lf ",data[i].imag/N);} 
  168.     printf("\n");*/  
  169.     /*for(i=0; i<N; i++) 
  170.         {printf("%lf ",sqrt(data[i].real*data[i].real+data[i].imag*data[i].imag)/N);}*/  
  171.     return 0;  
  172. }  


 

[cpp] view
plaincopy

  1. /******************************************************************************* 
  2. ** 程序名称：快速傅里叶变换(FFT)  
  3. ** 程序描述：本程序实现快速傅里叶变换  
  4. ** 性能提升：修正了IFFT的bug,使用宏定义改变N大小  
  5. ** 程序版本：V6.5 
  6. ** 程序作者：宋元瑞  
  7. ** 最后修改：2011年5月16日  
  8. *******************************************************************************/  
  9. #include <stdio.h>  
  10. #include <math.h>  
  11.   
  12. #define PI  3.14159265358979323846264338327950288419716939937510    //圆周率，50位小数   
  13. #define PI2 6.28318530717958647692528676655900576839433879875021  
  14. #define N 1024              //傅里叶变换的点数   
  15. #define M 10                //蝶形运算的级数，N = 2^M   
  16. #define Npart2 512          //创建正弦函数表时取PI的1/2   
  17. #define Npart4 256          //创建正弦函数表时取PI的1/4   
  18. #define FFT_RESULT(x)   (sqrt(data[x].real*data[x].real+data[x].imag*data[x].imag))  
  19. #define IFFT_RESULT(x)  (data[x].real/N)  
  20. typedef float ElemType;     //原始数据序列的数据类型,可以在这里设置  
  21.   
  22. typedef struct              //定义复数结构体   
  23. {  
  24.     ElemType real,imag;  
  25. }complex;  
  26.   
  27. complex data[N];            //定义存储单元，原始数据与负数结果均使用之   
  28. ElemType SIN_TABLE[Npart4+1];  
  29.   
  30. //产生模拟原始数据输入   
  31. void InputData(void)  
  32. {  
  33.     int i;  
  34.     for(i=0; i<N; i++)//制造输入序列   
  35.     {  
  36.         data[i].real = sin(2*PI*i/N);   //正弦波   
  37.         data[i].imag = 0;  
  38.     }  
  39. }  
  40. //创建正弦函数表   
  41. void CREATE_SIN_TABLE(void)  
  42. {  
  43.     int i=0;   
  44.     for(i=0; i<=Npart4; i++)  
  45.     {  
  46.         SIN_TABLE[i] = sin(PI*i/Npart2);//SIN_TABLE[i] = sin(PI2*i/N);  
  47.     }  
  48. }  
  49.   
  50. ElemType Sin_find(ElemType x)  
  51. {  
  52.     int i = (int)(N*x);  
  53.     i = i>>1;  
  54.     if(i>Npart4)//注意：i已经转化为0~N之间的整数了！   
  55.     {//不会超过N/2   
  56.         i = Npart2 - i;//i = i - 2*(i-Npart4);  
  57.     }   
  58.     return SIN_TABLE[i];  
  59. }  
  60. ElemType Cos_find(ElemType x)  
  61. {  
  62.     int i = (int)(N*x);  
  63.     i = i>>1;  
  64.     if(i<Npart4)//注意：i已经转化为0~N之间的整数了！   
  65.     {//不会超过N/2   
  66.         //i = Npart4 - i;  
  67.         return SIN_TABLE[Npart4 - i];  
  68.     }   
  69.     else//i>Npart4 && i<N/2  
  70.     {  
  71.         //i = i - Npart4;  
  72.         return -SIN_TABLE[i - Npart4];  
  73.     }  
  74. }  
  75.   
  76. //变址   
  77. void ChangeSeat(complex *DataInput)  
  78. {  
  79.     int nextValue,nextM,i,k,j=0;  
  80.     complex temp;  
  81.       
  82.     nextValue=N/2;                  //变址运算，即把自然顺序变成倒位序，采用雷德算法  
  83.     nextM=N-1;  
  84.     for (i=0;i<nextM;i++)  
  85.     {  
  86.         if (i<j)                 //如果i<j,即进行变址  
  87.         {  
  88.             temp=DataInput[j];  
  89.             DataInput[j]=DataInput[i];  
  90.             DataInput[i]=temp;  
  91.         }  
  92.         k=nextValue;                //求j的下一个倒位序  
  93.         while (k<=j)             //如果k<=j,表示j的最高位为1  
  94.         {  
  95.             j=j-k;                  //把最高位变成0  
  96.             k=k/2;                  //k/2，比较次高位，依次类推，逐个比较，直到某个位为0  
  97.         }  
  98.         j=j+k;                      //把0改为1  
  99.     }                                         
  100. }                                             
  101.   
  102. //复数乘法   
  103. /*complex XX_complex(complex a, complex b) 
  104. { 
  105.     complex temp; 
  106.      
  107.     temp.real = a.real * b.real-a.imag*b.imag; 
  108.     temp.imag = b.imag*a.real + a.imag*b.real; 
  109.      
  110.     return temp; 
  111. }*/  
  112.   
  113. //FFT运算函数   
  114. void FFT(void)  
  115. {  
  116.     int L=0,B=0,J=0,K=0;  
  117.     int step=0, KB=0;  
  118.     //ElemType P=0;  
  119.     ElemType angle;  
  120.     complex W,Temp_XX;  
  121.       
  122.     ChangeSeat(data);//变址   
  123.     //CREATE_SIN_TABLE();  
  124.     for(L=1; L<=M; L++)  
  125.     {  
  126.         step = 1<<L;//2^L  
  127.         B = step>>1;//B=2^(L-1)  
  128.         for(J=0; J<B; J++)  
  129.         {  
  130.             //P = (1<<(M-L))*J;//P=2^(M-L) *J   
  131.             angle = (double)J/B;            //这里还可以优化   
  132.             W.imag =  -Sin_find(angle);     //用C++该函数课声明为inline   
  133.             W.real =   Cos_find(angle);     //用C++该函数课声明为inline   
  134.             //W.real =  cos(angle*PI);  
  135.             //W.imag = -sin(angle*PI);  
  136.             for(K=J; K<N; K=K+step)  
  137.             {  
  138.                 KB = K + B;  
  139.                 //Temp_XX = XX_complex(data[KB],W);  
  140.                 //用下面两行直接计算复数乘法，省去函数调用开销   
  141.                 Temp_XX.real = data[KB].real * W.real-data[KB].imag*W.imag;  
  142.                 Temp_XX.imag = W.imag*data[KB].real + data[KB].imag*W.real;  
  143.                   
  144.                 data[KB].real = data[K].real - Temp_XX.real;  
  145.                 data[KB].imag = data[K].imag - Temp_XX.imag;  
  146.                   
  147.                 data[K].real = data[K].real + Temp_XX.real;  
  148.                 data[K].imag = data[K].imag + Temp_XX.imag;  
  149.             }  
  150.         }  
  151.     }  
  152. }  
  153. //IFFT运算函数   
  154. void IFFT(void)  
  155. {  
  156.     int L=0,B=0,J=0,K=0;  
  157.     int step=0, KB=0;  
  158.     //ElemType P=0;  
  159.     ElemType angle;  
  160.     complex W,Temp_XX;  
  161.       
  162.     ChangeSeat(data);//变址   
  163.     //CREATE_SIN_TABLE();  
  164.     for(L=1; L<=M; L++)  
  165.     {  
  166.         step = 1<<L;//2^L  
  167.         B = step>>1;//B=2^(L-1)  
  168.         for(J=0; J<B; J++)  
  169.         {  
  170.             //P = (1<<(M-L))*J;//P=2^(M-L) *J   
  171.             angle = (double)J/B;            //这里还可以优化   
  172.               
  173.             W.imag =   Sin_find(angle);     //用C++该函数课声明为inline   
  174.             W.real =   Cos_find(angle);     //用C++该函数课声明为inline   
  175.             //W.real =  cos(angle*PI);  
  176.             //W.imag = -sin(angle*PI);  
  177.             for(K=J; K<N; K=K+step)  
  178.             {  
  179.                 KB = K + B;  
  180.                 //Temp_XX = XX_complex(data[KB],W);  
  181.                 //用下面两行直接计算复数乘法，省去函数调用开销   
  182.                 Temp_XX.real = data[KB].real * W.real-data[KB].imag*W.imag;  
  183.                 Temp_XX.imag = W.imag*data[KB].real + data[KB].imag*W.real;  
  184.                   
  185.                 data[KB].real = data[K].real - Temp_XX.real;  
  186.                 data[KB].imag = data[K].imag - Temp_XX.imag;  
  187.                   
  188.                 data[K].real = data[K].real + Temp_XX.real;  
  189.                 data[K].imag = data[K].imag + Temp_XX.imag;  
  190.             }  
  191.         }  
  192.     }  
  193. }  
  194. //主函数   
  195. int main(int argc, char *argv[])  
  196. {  
  197.     int i = 0;  
  198.       
  199.     CREATE_SIN_TABLE(); //创建正弦函数表 ,这句只需在程序开始时执行一次   
  200.       
  201.     InputData();        //输入原始数据 ，此处用公式模拟；实际应用时为AD采样数据   
  202.     FFT();              //进行 FFT计算   
  203.       
  204.     for(i=0; i<N; i++)  
  205.         {printf("%f ",FFT_RESULT(i));}/**/  
  206.       
  207.     IFFT();//进行 IFFT计算   
  208.     for(i=0; i<N; i++)  
  209.         {printf("%f ",IFFT_RESULT(i));}/**/  
  210.       
  211.     return 0;  
  212. }  


 

[cpp] view
plaincopy

  1. /******************************************************************************* 
  2. ** 程序名称：快速傅里叶变换(FFT) 
  3. ** 程序描述：本程序实现快速傅里叶变换及其逆变换 
  4. ** 性能提升：修正了FFT的bug,变量重新命名,并将 N_FFT改为动态值 
  5. ** 程序版本：V6.6 
  6. ** 程序作者：宋元瑞 
  7. ** 最后修改：2011年5月16日 
  8. *******************************************************************************/  
  9. #include <stdio.h>  
  10. #include <stdlib.h>  
  11. #include <math.h>  
  12.   
  13. //#define OUTPRINT printf  
  14. //#define DEL /##/  
  15.   
  16. #define RESULT(x) sqrt(data_of_N_FFT[x].real*data_of_N_FFT[x].real+data_of_N_FFT[x].imag*data_of_N_FFT[x].imag)  
  17. #define PI  3.14159265358979323846264338327950288419716939937510    //圆周率，50位小数  
  18. #define PI2 6.28318530717958647692528676655900576839433879875021  
  19. int N_FFT=0;                //傅里叶变换的点数  
  20. int M_of_N_FFT=0;           //蝶形运算的级数，N = 2^M  
  21. int Npart2_of_N_FFT=0;      //创建正弦函数表时取PI的1/2  
  22. int Npart4_of_N_FFT=0;      //创建正弦函数表时取PI的1/4  
  23.   
  24. typedef float ElemType;     //原始数据序列的数据类型,可以在这里设置  
  25.   
  26. typedef struct              //定义复数结构体  
  27. {  
  28.     ElemType real,imag;  
  29. }complex_of_N_FFT,*ptr_complex_of_N_FFT;  
  30.   
  31. ptr_complex_of_N_FFT data_of_N_FFT=NULL;//开辟存储单元，原始数据与负数结果均使用之  
  32. ElemType* SIN_TABLE_of_N_FFT=NULL;  
  33.   
  34. //产生模拟原始数据输入  
  35. void InputData(void)  
  36. {  
  37.     int i;  
  38.     for (i=0; i<N_FFT; i++)//制造输入序列  
  39.     {  
  40.         data_of_N_FFT[i].real = sin(2*PI*i/N_FFT);  //正弦波  
  41.         data_of_N_FFT[i].imag = 0;  
  42.         printf("%f ",data_of_N_FFT[i].real);  
  43.     }  
  44. }  
  45.   
  46. //创建正弦函数表  
  47. void CREATE_SIN_TABLE(void)  
  48. {  
  49.     int i=0;  
  50.     for (i=0; i<=Npart4_of_N_FFT; i++)  
  51.     {  
  52.         SIN_TABLE_of_N_FFT[i] = sin(PI*i/Npart2_of_N_FFT);//SIN_TABLE[i] = sin(PI2*i/N);  
  53.     }  
  54. }  
  55.   
  56.   
  57. ElemType Sin_find(ElemType x)  
  58. {  
  59.     int i = (int)(N_FFT*x);  
  60.     i = i>>1;  
  61.     if (i>Npart4_of_N_FFT)//注意：i已经转化为0~N之间的整数了！  
  62.     {  
  63.         //不会超过N/2  
  64.         i = Npart2_of_N_FFT - i;//i = i - 2*(i-Npart4);  
  65.     }  
  66.     return SIN_TABLE_of_N_FFT[i];  
  67. }  
  68.   
  69. ElemType Cos_find(ElemType x)  
  70. {  
  71.     int i = (int)(N_FFT*x);  
  72.     i = i>>1;  
  73.     if (i<Npart4_of_N_FFT)//注意：i已经转化为0~N之间的整数了！  
  74.     {  
  75.         //不会超过N/2  
  76.         //i = Npart4 - i;  
  77.         return SIN_TABLE_of_N_FFT[Npart4_of_N_FFT - i];  
  78.     }  
  79.     else//i>Npart4 && i<N/2  
  80.     {  
  81.         //i = i - Npart4;  
  82.         return -SIN_TABLE_of_N_FFT[i - Npart4_of_N_FFT];  
  83.     }  
  84. }  
  85.   
  86. //变址  
  87. void ChangeSeat(complex_of_N_FFT *DataInput)  
  88. {  
  89.     int nextValue,nextM,i,k,j=0;  
  90.     complex_of_N_FFT temp;  
  91.   
  92.     nextValue=N_FFT/2;                  //变址运算，即把自然顺序变成倒位序，采用雷德算法  
  93.     nextM=N_FFT-1;  
  94.     for (i=0;i<nextM;i++)  
  95.     {  
  96.         if (i<j)                 //如果i<j,即进行变址  
  97.         {  
  98.             temp=DataInput[j];  
  99.             DataInput[j]=DataInput[i];  
  100.             DataInput[i]=temp;  
  101.         }  
  102.         k=nextValue;                //求j的下一个倒位序  
  103.         while (k<=j)             //如果k<=j,表示j的最高位为1  
  104.         {  
  105.             j=j-k;                  //把最高位变成0  
  106.             k=k/2;                  //k/2，比较次高位，依次类推，逐个比较，直到某个位为0  
  107.         }  
  108.         j=j+k;                      //把0改为1  
  109.     }  
  110. }  
  111.   
  112. //复数乘法  
  113. /*complex XX_complex(complex a, complex b) 
  114. { 
  115.     complex temp; 
  116.  
  117.     temp.real = a.real * b.real-a.imag*b.imag; 
  118.     temp.imag = b.imag*a.real + a.imag*b.real; 
  119.  
  120.     return temp; 
  121. }*/  
  122.   
  123. //FFT运算函数  
  124. void FFT(void)  
  125. {  
  126.     int L=0,B=0,J=0,K=0;  
  127.     int step=0, KB=0;  
  128.     //ElemType P=0;  
  129.     ElemType angle;  
  130.     complex_of_N_FFT W,Temp_XX;  
  131.   
  132.     ChangeSeat(data_of_N_FFT);//变址  
  133.     //CREATE_SIN_TABLE();  
  134.     for (L=1; L<=M_of_N_FFT; L++)  
  135.     {  
  136.         step = 1<<L;//2^L  
  137.         B = step>>1;//B=2^(L-1)  
  138.         for (J=0; J<B; J++)  
  139.         {  
  140.             //P = (1<<(M-L))*J;//P=2^(M-L) *J  
  141.             angle = (double)J/B;            //这里还可以优化  
  142.             W.imag =  -Sin_find(angle);     //用C++该函数课声明为inline  
  143.             W.real =   Cos_find(angle);     //用C++该函数课声明为inline  
  144.             //W.real =  cos(angle*PI);  
  145.             //W.imag = -sin(angle*PI);  
  146.             for (K=J; K<N_FFT; K=K+step)  
  147.             {  
  148.                 KB = K + B;  
  149.                 //Temp_XX = XX_complex(data[KB],W);  
  150.                 //用下面两行直接计算复数乘法，省去函数调用开销  
  151.                 Temp_XX.real = data_of_N_FFT[KB].real * W.real-data_of_N_FFT[KB].imag*W.imag;  
  152.                 Temp_XX.imag = W.imag*data_of_N_FFT[KB].real + data_of_N_FFT[KB].imag*W.real;  
  153.   
  154.                 data_of_N_FFT[KB].real = data_of_N_FFT[K].real - Temp_XX.real;  
  155.                 data_of_N_FFT[KB].imag = data_of_N_FFT[K].imag - Temp_XX.imag;  
  156.   
  157.                 data_of_N_FFT[K].real = data_of_N_FFT[K].real + Temp_XX.real;  
  158.                 data_of_N_FFT[K].imag = data_of_N_FFT[K].imag + Temp_XX.imag;  
  159.             }  
  160.         }  
  161.     }  
  162. }  
  163.   
  164. //IFFT运算函数  
  165. void IFFT(void)  
  166. {  
  167.     int L=0,B=0,J=0,K=0;  
  168.     int step=0, KB=0;  
  169.     //ElemType P=0;  
  170.     ElemType angle;  
  171.     complex_of_N_FFT W,Temp_XX;  
  172.   
  173.     ChangeSeat(data_of_N_FFT);//变址  
  174.     //CREATE_SIN_TABLE();  
  175.     for (L=1; L<=M_of_N_FFT; L++)  
  176.     {  
  177.         step = 1<<L;//2^L  
  178.         B = step>>1;//B=2^(L-1)  
  179.         for (J=0; J<B; J++)  
  180.         {  
  181.             //P = (1<<(M-L))*J;//P=2^(M-L) *J  
  182.             angle = (double)J/B;            //这里还可以优化  
  183.   
  184.             W.imag =   Sin_find(angle);     //用C++该函数课声明为inline  
  185.             W.real =   Cos_find(angle);     //用C++该函数课声明为inline  
  186.             //W.real =  cos(angle*PI);  
  187.             //W.imag = -sin(angle*PI);  
  188.             for (K=J; K<N_FFT; K=K+step)  
  189.             {  
  190.                 KB = K + B;  
  191.                 //Temp_XX = XX_complex(data[KB],W);  
  192.                 //用下面两行直接计算复数乘法，省去函数调用开销  
  193.                 Temp_XX.real = data_of_N_FFT[KB].real * W.real-data_of_N_FFT[KB].imag*W.imag;  
  194.                 Temp_XX.imag = W.imag*data_of_N_FFT[KB].real + data_of_N_FFT[KB].imag*W.real;  
  195.   
  196.                 data_of_N_FFT[KB].real = data_of_N_FFT[K].real - Temp_XX.real;  
  197.                 data_of_N_FFT[KB].imag = data_of_N_FFT[K].imag - Temp_XX.imag;  
  198.   
  199.                 data_of_N_FFT[K].real = data_of_N_FFT[K].real + Temp_XX.real;  
  200.                 data_of_N_FFT[K].imag = data_of_N_FFT[K].imag + Temp_XX.imag;  
  201.             }  
  202.         }  
  203.     }  
  204. }  
  205.   
  206. //初始化FFT程序  
  207. //N_FFT是 FFT的点数，必须是2的次方  
  208. void Init_FFT(int N_of_FFT)  
  209. {  
  210.     int i=0;  
  211.     int temp_N_FFT=1;  
  212.     N_FFT = N_of_FFT;                   //傅里叶变换的点数 ，必须是 2的次方  
  213.     M_of_N_FFT = 0;                 //蝶形运算的级数，N = 2^M  
  214.     for (i=0; temp_N_FFT<N_FFT; i++)  
  215.     {  
  216.         temp_N_FFT = 2*temp_N_FFT;  
  217.         M_of_N_FFT++;  
  218.     }  
  219.     printf("\n%d\n",M_of_N_FFT);  
  220.     Npart2_of_N_FFT = N_FFT/2;      //创建正弦函数表时取PI的1/2  
  221.     Npart4_of_N_FFT = N_FFT/4;      //创建正弦函数表时取PI的1/4  
  222.   
  223.     //ptr_complex_of_N_FFT data_of_N_FFT=NULL;//开辟存储单元，原始数据与负数结果均使用之  
  224.     data_of_N_FFT = (ptr_complex_of_N_FFT)malloc(N_FFT * sizeof(complex_of_N_FFT));  
  225.     //ptr_complex_of_N_FFT SIN_TABLE_of_N_FFT=NULL;  
  226.     SIN_TABLE_of_N_FFT = (ElemType *)malloc((Npart4_of_N_FFT+1) * sizeof(ElemType));  
  227.   
  228.     CREATE_SIN_TABLE();             //创建正弦函数表  
  229. }  
  230.   
  231. //结束 FFT运算，释放相关内存  
  232. void Close_FFT(void)  
  233. {  
  234.     if (data_of_N_FFT != NULL)  
  235.     {  
  236.         free(data_of_N_FFT);          //释放内存  
  237.         data_of_N_FFT = NULL;  
  238.     }  
  239.     if (SIN_TABLE_of_N_FFT != NULL)  
  240.     {  
  241.         free(SIN_TABLE_of_N_FFT);     //释放内存  
  242.         SIN_TABLE_of_N_FFT = NULL;  
  243.     }  
  244. }  
  245.   
  246. //主函数  
  247. int main(int argc, char *argv[])  
  248. {  
  249.     int i = 0;  
  250.   
  251.     Init_FFT(1024);     //初始化各项参数，此函数只需执行一次  
  252.   
  253.     InputData();        //输入原始数据  
  254.     FFT();              //进行 FFT计算  
  255.   
  256.     printf("\n\n");  
  257.     for (i=0; i<N_FFT; i++)  
  258.     {  
  259.         printf("%f ",RESULT(i));  
  260.     }  
  261.   
  262.     IFFT();//进行 IFFT计算  
  263.     printf("\n\n");  
  264.     for (i=0; i<N_FFT; i++)  
  265.     {  
  266.         printf("%f ",data_of_N_FFT[i].real/N_FFT);  
  267.     }  
  268.   
  269.     Close_FFT();        //结束 FFT运算，释放相关内存  
  270.   
  271.     return 0;  
  272. }  


 

[cpp] view
plaincopy

  1. /******************************************************************************* 
  2. ** 模块名称：快速傅里叶变换(FFT)  
  3. ** 模块描述：本程序实现快速傅里叶变换  
  4. ** 性能提升：已达到网上同类程序最高速度  
  5. ** 模块版本：V6.0 
  6. ** 模块作者：宋元瑞  
  7. ** 最后修改：2011年5月6日 
  8. *******************************************************************************/  
  9. /******************************************************************************* 
  10. **  程序说明： 
  11.     FFT运算输入参数 source_Data(i) 是一个N大小的数组（注意是小括号） 
  12.     FFT运算输出结果 result_Data(i) 是一个N大小的数组（注意是小括号） 
  13. **  调用举例： 
  14. int main(int argc, char *argv[]) 
  15. { 
  16.     int i = 0; 
  17.     ///以下为FFT运算的调用方式  
  18.     FFT_prepare();      //为FFT做好准备，此函数只需程序开始时执行一次即可，切勿写在循环中  
  19.     while(1) 
  20.     { 
  21.         for(i=0;i<N_FFT;i++) //输入原始数据  
  22.             {source_Data(i) = sin(2*PI*i/N_FFT);}//注意inputData后面是小括号  
  23.         FFT();              //进行FFT计算  
  24.         //输出结果：XXX =  result_Data(i); 
  25.     } 
  26.     return 0; 
  27. } 
  28. *******************************************************************************/   
  29. #ifndef SYRFFT_H  
  30. //#pragma once  
  31.   
  32. //#include <stdio.h>  
  33. #include <math.h>  
  34.   
  35. #define FFT_prepare() CREATE_SIN_TABLE();   //创建正弦函数表   
  36. #define source_Data(i) data_FFT[i].imag     //接收输入数据的数组，大小为N   
  37. #define result_Data(i) sqrt(data_FFT[i].real*data_FFT[i].real+data_FFT[i].imag*data_FFT[i].imag)//FFT结果   
  38.   
  39. #define PI  3.14159265358979323846264338327950288419716939937510    //圆周率，50位小数   
  40. #define PI2 6.28318530717958647692528676655900576839433879875021  
  41. #define N_FFT 1024              //傅里叶变换的点数   
  42. #define M_of_N_FFT 10               //蝶形运算的级数，N = 2^M   
  43. #define Npart2_of_N_FFT 512         //创建正弦函数表时取PI的1/2   
  44. #define Npart4_of_N_FFT 256         //创建正弦函数表时取PI的1/4   
  45.   
  46. typedef float ElemType;     //原始数据序列的数据类型,可以在这里设置  
  47.   
  48. typedef struct              //定义复数结构体   
  49. {  
  50.     ElemType real,imag;  
  51. }complex_of_N_FFT;  
  52.   
  53. complex_of_N_FFT data_FFT[N_FFT];       //存放要进行FFT运输的数据，运算结果也存放这里   
  54. ElemType SIN_TABLE[Npart4_of_N_FFT+1];  
  55.   
  56. //产生模拟原始数据输入   
  57. /* 
  58. void InputData(void) 
  59. { 
  60.     int i; 
  61.     for(i=0; i<N_FFT; i++)       //制造输入序列  
  62.     { 
  63.         source_Data(i) = sin(2*PI*i/N_FFT); //模拟输入正弦波  
  64.         //data_FFT[i].imag = sin(2*PI*i/N); //正弦波  
  65.         //printf("%f ",data_FFT[i].imag); 
  66.     } 
  67. }*/  
  68. //创建正弦函数表   
  69. void CREATE_SIN_TABLE(void)  
  70. {  
  71.     int i=0;   
  72.     for(i=0; i<=Npart4_of_N_FFT; i++)  
  73.     {  
  74.         SIN_TABLE[i] = sin(PI*i/Npart2_of_N_FFT);//SIN_TABLE[i] = sin(PI2*i/N);  
  75.     }  
  76. }  
  77.   
  78. ElemType Sin_find(ElemType x)  
  79. {  
  80.     int i = (int)(N_FFT*x);  
  81.     i = i>>1;  
  82.     if(i>Npart4_of_N_FFT)//注意：i已经转化为0~N之间的整数了！   
  83.     {//不会超过N/2   
  84.         i = Npart2_of_N_FFT - i;//i = i - 2*(i-Npart4);  
  85.     }   
  86.     return SIN_TABLE[i];  
  87. }  
  88. ElemType Cos_find(ElemType x)  
  89. {  
  90.     int i = (int)(N_FFT*x);  
  91.     i = i>>1;  
  92.     if(i<Npart4_of_N_FFT)//注意：i已经转化为0~N之间的整数了！   
  93.     {//不会超过N/2   
  94.         //i = Npart4 - i;  
  95.         return SIN_TABLE[Npart4_of_N_FFT - i];  
  96.     }   
  97.     else//i>Npart4 && i<N/2  
  98.     {  
  99.         //i = i - Npart4;  
  100.         return -SIN_TABLE[i - Npart4_of_N_FFT];  
  101.     }  
  102. }  
  103.   
  104. //变址   
  105. void ChangeSeat(complex_of_N_FFT *DataInput)  
  106. {  
  107.     //变址前data_FFT[i].imag存储了输入数据，变址后data_FFT[i].real存储了输入数据  
  108.     int i,n,New_seat;  
  109.       
  110.     for(i=0; i<N_FFT; i++)  
  111.     {  
  112.         New_seat = 0;  
  113.         for(n=0;n<M_of_N_FFT;n++)  
  114.         {  
  115.             New_seat = New_seat | (((i>>n) & 0x01) << (M_of_N_FFT-n-1));  
  116.         }  
  117.         DataInput[New_seat].real = DataInput[i].imag;  
  118.     }  
  119.     for(i=0; i<N_FFT; i++)  
  120.     {  
  121.         DataInput[i].imag = 0;  
  122.     }  
  123. }                                                        
  124.   
  125. //复数乘法   
  126. /* 
  127. complex_of_N_FFT XX_complex(complex_of_N_FFT a, complex_of_N_FFT b) 
  128. { 
  129.     complex_of_N_FFT temp; 
  130.      
  131.     temp.real = a.real * b.real-a.imag*b.imag; 
  132.     temp.imag = b.imag*a.real + a.imag*b.real; 
  133.      
  134.     return temp; 
  135. }*/  
  136.   
  137. //FFT运算函数   
  138. void FFT(void)  
  139. {  
  140.     int L=0,B=0,J=0,K=0;  
  141.     int step=0, KB=0;  
  142.     //ElemType P=0;  
  143.     ElemType angle;  
  144.     complex_of_N_FFT W,Temp_XX;  
  145.       
  146.     ChangeSeat(data_FFT);   //变址   
  147.       
  148.     for(L=1; L<=M_of_N_FFT; L++)  
  149.     {  
  150.         step = 1<<L;//2^L  
  151.         B = step>>1;//B=2^(L-1)  
  152.         for(J=0; J<B; J++)  
  153.         {  
  154.             //P = (1<<(M-L))*J;//P=2^(M-L) *J   
  155.             angle = (double)J/B;            //这里还可以优化   
  156.             W.imag =  -Sin_find(angle);     //用C++该函数课声明为inline   
  157.             W.real =   Cos_find(angle);     //用C++该函数课声明为inline   
  158.             //W.real =  cos(angle*PI);  
  159.             //W.imag = -sin(angle*PI);  
  160.             for(K=J; K<N_FFT; K=K+step)  
  161.             {  
  162.                 KB = K + B;  
  163.                 //Temp_XX = XX_complex(data_FFT[KB],W);  
  164.                 //用下面两行直接计算复数乘法，省去函数调用开销   
  165.                 Temp_XX.real = data_FFT[KB].real * W.real-data_FFT[KB].imag*W.imag;  
  166.                 Temp_XX.imag = W.imag*data_FFT[KB].real + data_FFT[KB].imag*W.real;  
  167.                   
  168.                 data_FFT[KB].real = data_FFT[K].real - Temp_XX.real;  
  169.                 data_FFT[KB].imag = data_FFT[K].imag - Temp_XX.imag;  
  170.                   
  171.                 data_FFT[K].real = data_FFT[K].real + Temp_XX.real;  
  172.                 data_FFT[K].imag = data_FFT[K].imag + Temp_XX.imag;  
  173.             }  
  174.         }  
  175.     }  
  176. }  
  177.   
  178. #define SYRFFT_H  
  179. #endif  

 

[cpp] view
plaincopy

  1. /******************************************************************************* 
  2. ** 程序名称：快速傅里叶变换(FFT)  
  3. ** 程序描述：本程序实现快速傅里叶变换  
  4. ** 程序版本：V6.0 
  5. ** 程序作者：宋元瑞  
  6. ** 最后修改：2011年5月6日 
  7. *******************************************************************************/  
  8.   
  9. #include <stdio.h>  
  10. #include "syrFFT_H.h"   //包含FFT运算头文件   
  11.   
  12. //主函数   
  13. int main(int argc, char *argv[])  
  14. {  
  15.     int i = 0;  
  16.       
  17.     //以下3句为FFT运算的调用函数   
  18.     FFT_prepare();      //为FFT做好准备，此函数只需程序开始时执行一次即可，切勿写在循环中   
  19.     //while(1)  
  20.     {  
  21.         for(i=0;i<N_FFT;i++)//模拟输入   
  22.             {source_Data(i) = sin(2*PI*i/N_FFT);}//注意inputData后面是小括号   
  23.         FFT();  
  24.   
  25.         for(i=0; i<N_FFT; i++)//输出结果   
  26.             {printf("%lf ",result_Data(i));}  
  27.     }   
  28.       
  29.     return 0;  
  30. }  

 


 

[cpp] view
plaincopy

  1. #ifndef FFT_H  
  2.   
  3. #include <stdlib.h>  
  4. #include <math.h>  
  5.   
  6. #define RESULT(x) sqrt(data_of_N_FFT[x].real*data_of_N_FFT[x].real+data_of_N_FFT[x].imag*data_of_N_FFT[x].imag)  
  7. #define PI  3.14159265358979323846264338327950288419716939937510    //圆周率，50位小数   
  8. #define PI2 6.28318530717958647692528676655900576839433879875021  
  9. int N_FFT=0;                //傅里叶变换的点数   
  10. int M_of_N_FFT=0;           //蝶形运算的级数，N = 2^M   
  11. int Npart2_of_N_FFT=0;      //创建正弦函数表时取PI的1/2   
  12. int Npart4_of_N_FFT=0;      //创建正弦函数表时取PI的1/4   
  13.   
  14. typedef float ElemType;     //原始数据序列的数据类型,可以在这里设置  
  15.   
  16. typedef struct              //定义复数结构体   
  17. {  
  18.     ElemType real,imag;  
  19. }complex_of_N_FFT,*ptr_complex_of_N_FFT;  
  20.   
  21. ptr_complex_of_N_FFT data_of_N_FFT=NULL;//开辟存储单元，原始数据与负数结果均使用之   
  22. ElemType* SIN_TABLE_of_N_FFT=NULL;  
  23.   
  24.   
  25. //创建正弦函数表   
  26. void CREATE_SIN_TABLE(void)  
  27. {  
  28.     int i=0;   
  29.     for(i=0; i<=Npart4_of_N_FFT; i++)  
  30.     {  
  31.         SIN_TABLE_of_N_FFT[i] = sin(PI*i/Npart2_of_N_FFT);//SIN_TABLE[i] = sin(PI2*i/N);  
  32.     }  
  33. }  
  34.   
  35. ElemType Sin_find(ElemType x)  
  36. {  
  37.     int i = (int)(N_FFT*x);  
  38.     i = i>>1;  
  39.     if(i>Npart4_of_N_FFT)//注意：i已经转化为0~N之间的整数了！   
  40.     {//不会超过N/2   
  41.         i = Npart2_of_N_FFT - i;//i = i - 2*(i-Npart4);  
  42.     }   
  43.     return SIN_TABLE_of_N_FFT[i];  
  44. }  
  45. ElemType Cos_find(ElemType x)  
  46. {  
  47.     int i = (int)(N_FFT*x);  
  48.     i = i>>1;  
  49.     if(i<Npart4_of_N_FFT)//注意：i已经转化为0~N之间的整数了！   
  50.     {//不会超过N/2   
  51.         //i = Npart4 - i;  
  52.         return SIN_TABLE_of_N_FFT[Npart4_of_N_FFT - i];  
  53.     }   
  54.     else//i>Npart4 && i<N/2  
  55.     {  
  56.         //i = i - Npart4;  
  57.         return -SIN_TABLE_of_N_FFT[i - Npart4_of_N_FFT];  
  58.     }  
  59. }  
  60.   
  61. //变址   
  62. void ChangeSeat(complex_of_N_FFT *DataInput)  
  63. {  
  64.     int nextValue,nextM,i,k,j=0;  
  65.     complex_of_N_FFT temp;  
  66.       
  67.     nextValue=N_FFT/2;                  //变址运算，即把自然顺序变成倒位序，采用雷德算法  
  68.     nextM=N_FFT-1;  
  69.     for (i=0;i<nextM;i++)  
  70.     {  
  71.         if (i<j)                 //如果i<j,即进行变址  
  72.         {  
  73.             temp=DataInput[j];  
  74.             DataInput[j]=DataInput[i];  
  75.             DataInput[i]=temp;  
  76.         }  
  77.         k=nextValue;                //求j的下一个倒位序  
  78.         while (k<=j)             //如果k<=j,表示j的最高位为1  
  79.         {  
  80.             j=j-k;                  //把最高位变成0  
  81.             k=k/2;                  //k/2，比较次高位，依次类推，逐个比较，直到某个位为0  
  82.         }  
  83.         j=j+k;                      //把0改为1  
  84.     }                                         
  85. }                                             
  86.   
  87. //复数乘法   
  88. /*complex XX_complex(complex a, complex b) 
  89. { 
  90.     complex temp; 
  91.      
  92.     temp.real = a.real * b.real-a.imag*b.imag; 
  93.     temp.imag = b.imag*a.real + a.imag*b.real; 
  94.      
  95.     return temp; 
  96. }*/  
  97.   
  98. //FFT运算函数   
  99. void FFT(void)  
  100. {  
  101.     int L=0,B=0,J=0,K=0;  
  102.     int step=0, KB=0;  
  103.     //ElemType P=0;  
  104.     ElemType angle;  
  105.     complex_of_N_FFT W,Temp_XX;  
  106.       
  107.     ChangeSeat(data_of_N_FFT);//变址   
  108.     //CREATE_SIN_TABLE();  
  109.     for(L=1; L<=M_of_N_FFT; L++)  
  110.     {  
  111.         step = 1<<L;//2^L  
  112.         B = step>>1;//B=2^(L-1)  
  113.         for(J=0; J<B; J++)  
  114.         {  
  115.             //P = (1<<(M-L))*J;//P=2^(M-L) *J   
  116.             angle = (double)J/B;            //这里还可以优化   
  117.             W.imag =  -Sin_find(angle);     //用C++该函数课声明为inline   
  118.             W.real =   Cos_find(angle);     //用C++该函数课声明为inline   
  119.             //W.real =  cos(angle*PI);  
  120.             //W.imag = -sin(angle*PI);  
  121.             for(K=J; K<N_FFT; K=K+step)  
  122.             {  
  123.                 KB = K + B;  
  124.                 //Temp_XX = XX_complex(data[KB],W);  
  125.                 //用下面两行直接计算复数乘法，省去函数调用开销   
  126.                 Temp_XX.real = data_of_N_FFT[KB].real * W.real-data_of_N_FFT[KB].imag*W.imag;  
  127.                 Temp_XX.imag = W.imag*data_of_N_FFT[KB].real + data_of_N_FFT[KB].imag*W.real;  
  128.                   
  129.                 data_of_N_FFT[KB].real = data_of_N_FFT[K].real - Temp_XX.real;  
  130.                 data_of_N_FFT[KB].imag = data_of_N_FFT[K].imag - Temp_XX.imag;  
  131.                   
  132.                 data_of_N_FFT[K].real = data_of_N_FFT[K].real + Temp_XX.real;  
  133.                 data_of_N_FFT[K].imag = data_of_N_FFT[K].imag + Temp_XX.imag;  
  134.             }  
  135.         }  
  136.     }  
  137. }  
  138. //IFFT运算函数   
  139. void IFFT(void)  
  140. {  
  141.     int L=0,B=0,J=0,K=0;  
  142.     int step=0, KB=0;  
  143.     //ElemType P=0;  
  144.     ElemType angle;  
  145.     complex_of_N_FFT W,Temp_XX;  
  146.       
  147.     ChangeSeat(data_of_N_FFT);//变址   
  148.     //CREATE_SIN_TABLE();  
  149.     for(L=1; L<=M_of_N_FFT; L++)  
  150.     {  
  151.         step = 1<<L;//2^L  
  152.         B = step>>1;//B=2^(L-1)  
  153.         for(J=0; J<B; J++)  
  154.         {  
  155.             //P = (1<<(M-L))*J;//P=2^(M-L) *J   
  156.             angle = (double)J/B;            //这里还可以优化   
  157.               
  158.             W.imag =   Sin_find(angle);     //用C++该函数课声明为inline   
  159.             W.real =   Cos_find(angle);     //用C++该函数课声明为inline   
  160.             //W.real =  cos(angle*PI);  
  161.             //W.imag = -sin(angle*PI);  
  162.             for(K=J; K<N_FFT; K=K+step)  
  163.             {  
  164.                 KB = K + B;  
  165.                 //Temp_XX = XX_complex(data[KB],W);  
  166.                 //用下面两行直接计算复数乘法，省去函数调用开销   
  167.                 Temp_XX.real = data_of_N_FFT[KB].real * W.real-data_of_N_FFT[KB].imag*W.imag;  
  168.                 Temp_XX.imag = W.imag*data_of_N_FFT[KB].real + data_of_N_FFT[KB].imag*W.real;  
  169.                   
  170.                 data_of_N_FFT[KB].real = data_of_N_FFT[K].real - Temp_XX.real;  
  171.                 data_of_N_FFT[KB].imag = data_of_N_FFT[K].imag - Temp_XX.imag;  
  172.                   
  173.                 data_of_N_FFT[K].real = data_of_N_FFT[K].real + Temp_XX.real;  
  174.                 data_of_N_FFT[K].imag = data_of_N_FFT[K].imag + Temp_XX.imag;  
  175.             }  
  176.         }  
  177.     }  
  178. }  
  179.   
  180. //初始化FFT程序   
  181. //N_FFT是 FFT的点数，必须是2的次方   
  182. void Init_FFT(int N_of_FFT)  
  183. {  
  184.     int i=0;  
  185.     int temp_N_FFT=1;  
  186.     N_FFT = N_of_FFT;                   //傅里叶变换的点数 ，必须是 2的次方   
  187.     M_of_N_FFT = 0;                 //蝶形运算的级数，N = 2^M   
  188.     for(i=0; temp_N_FFT<N_FFT; i++)  
  189.     {  
  190.         temp_N_FFT = 2*temp_N_FFT;  
  191.         M_of_N_FFT++;  
  192.     }  
  193.     //printf("\n%d\n",M_of_N_FFT);  
  194.     Npart2_of_N_FFT = N_FFT/2;      //创建正弦函数表时取PI的1/2   
  195.     Npart4_of_N_FFT = N_FFT/4;      //创建正弦函数表时取PI的1/4   
  196.       
  197.     //ptr_complex_of_N_FFT data_of_N_FFT=NULL;//开辟存储单元，原始数据与负数结果均使用之   
  198.     data_of_N_FFT = (ptr_complex_of_N_FFT)malloc(N_FFT * sizeof(complex_of_N_FFT));  
  199.     //ptr_complex_of_N_FFT SIN_TABLE_of_N_FFT=NULL;  
  200.     SIN_TABLE_of_N_FFT = (ElemType *)malloc((Npart4_of_N_FFT+1) * sizeof(ElemType));  
  201.   
  202.     CREATE_SIN_TABLE();             //创建正弦函数表   
  203. }  
  204. //结束 FFT运算，释放相关内存   
  205. void Close_FFT(void)  
  206. {  
  207.     if(data_of_N_FFT != NULL)  
  208.     {  
  209.         free(data_of_N_FFT);          //释放内存  
  210.         data_of_N_FFT = NULL;  
  211.     }  
  212.     if(SIN_TABLE_of_N_FFT != NULL)  
  213.     {  
  214.         free(SIN_TABLE_of_N_FFT);     //释放内存  
  215.         SIN_TABLE_of_N_FFT = NULL;  
  216.     }  
  217. }  
  218.   
  219. #define FFT_H  
  220. #endif  
[cpp] view
plaincopy

  1. /******************************************************************************* 
  2. ** 程序名称：快速傅里叶变换(FFT)  
  3. ** 程序描述：本程序实现快速傅里叶变换及其逆变换  
  4. ** 性能提升：修正了FFT的bug  
  5. ** 程序版本：V6.6 
  6. ** 程序作者：宋元瑞  
  7. ** 最后修改：2011年5月16日  
  8. *******************************************************************************/  
  9.   
  10. #include "jxust_fft6_6.h"  
  11.   
  12. //产生模拟原始数据输入 ,在实际应用时替换为AD采样数据   
  13. void InputData(void)  
  14. {  
  15.     int i;  
  16.     for(i=0; i<N_FFT; i++)//制造输入序列   
  17.     {  
  18.         data_of_N_FFT[i].real = sin(2*PI*i/N_FFT);  //输入采样数据   
  19.         data_of_N_FFT[i].imag = 0;  
  20.     }  
  21. }  
  22.   
  23. //主函数 ，示例如何调用FFT运算   
  24. int main(int argc, char *argv[])  
  25. {  
  26.     int i = 0;  
  27.       
  28.     Init_FFT(1024);     //①初始化各项参数，此函数只需执行一次 ; 参数为FFT的点数，必须为2的次方   
  29.       
  30.     InputData();        //②输入原始数据 ,此处在实际应用时替换为AD采样数据   
  31.     FFT();              //③进行 FFT计算   
  32.     //for(i=0; i<N_FFT; i++)//④输出FFT频谱结果  sqrt(data_of_N_FFT[i].real*data_of_N_FFT[i].real+data_of_N_FFT[i].imag*data_of_N_FFT[i].imag)  
  33.         //{printf("%f ",RESULT(i));}  
  34.       
  35.     IFFT();//进行 IFFT计算   
  36.     //for(i=0; i<N_FFT; i++)//(可选步骤)⑤输出 IFFT结果 ，滤波时会用到；暂时不用   
  37.         //{printf("%f ",data_of_N_FFT[i].real/N_FFT);}  
  38.     Close_FFT();    //结束 FFT运算，释放相关内存 ;此函数在彻底结束FFT运算后调用，   
  39.       
  40.     return 0;  
  41. }  


 

[cpp] view
plaincopy

  1. #ifndef SYRFFT_6_55_H  
  2.   
  3. #include <math.h>  
  4.   
  5. #define FFT_RESULT(x)   (sqrt(data_of_N_FFT[x].real*data_of_N_FFT[x].real+data_of_N_FFT[x].imag*data_of_N_FFT[x].imag))  
  6. #define IFFT_RESULT(x)  (data_of_N_FFT[x].real/N_FFT)  
  7.   
  8. #define PI  3.14159265358979323846264338327950288419716939937510    //圆周率，50位小数   
  9. #define PI2 6.28318530717958647692528676655900576839433879875021  
  10. #define N_FFT           1024        //傅里叶变换的点数   
  11. #define M_of_N_FFT      10          //蝶形运算的级数，N = 2^M   
  12. #define Npart2_of_N_FFT 512         //创建正弦函数表时取PI的1/2   
  13. #define Npart4_of_N_FFT 256         //创建正弦函数表时取PI的1/4   
  14.   
  15. typedef float ElemType;     //原始数据序列的数据类型,可以在这里设置  
  16.   
  17. typedef struct              //定义复数结构体   
  18. {  
  19.     ElemType real,imag;  
  20. }complex_of_N_FFT,*ptr_complex_of_N_FFT;  
  21.   
  22. complex_of_N_FFT data_of_N_FFT[N_FFT];          //定义存储单元，原始数据与负数结果均使用之   
  23. ElemType SIN_TABLE_of_N_FFT[Npart4_of_N_FFT+1];  
  24.   
  25. //创建正弦函数表   
  26. void CREATE_SIN_TABLE(void)  
  27. {  
  28.     int i=0;   
  29.     for(i=0; i<=Npart4_of_N_FFT; i++)  
  30.     {  
  31.         SIN_TABLE_of_N_FFT[i] = sin(PI*i/Npart2_of_N_FFT);//SIN_TABLE[i] = sin(PI2*i/N);  
  32.     }  
  33. }  
  34.   
  35. ElemType Sin_find(ElemType x)  
  36. {  
  37.     int i = (int)(N_FFT*x);//注意：i已经转化为0~N之间的整数了！   
  38.     i = i>>1;// i = i / 2;  
  39.     if(i>Npart4_of_N_FFT)  
  40.     {//根据FFT相关公式，sin()参数不会超过PI， 即i不会超过N/2   
  41.         i = Npart2_of_N_FFT - i;//i = i - 2*(i-Npart4);  
  42.     }   
  43.     return SIN_TABLE_of_N_FFT[i];  
  44. }  
  45. ElemType Cos_find(ElemType x)  
  46. {  
  47.     int i = (int)(N_FFT*x);//注意：i已经转化为0~N之间的整数了！   
  48.     i = i>>1;  
  49.     if(i < Npart4_of_N_FFT )   
  50.     { //不会超过N/2   
  51.         //i = Npart4 - i;  
  52.         return SIN_TABLE_of_N_FFT[Npart4_of_N_FFT - i];  
  53.     }   
  54.     else //i > Npart4 && i < N/2   
  55.     {  
  56.         //i = i - Npart4;  
  57.         return -SIN_TABLE_of_N_FFT[i - Npart4_of_N_FFT];  
  58.     }  
  59. }  
  60.   
  61. //变址   
  62. void ChangeSeat(complex_of_N_FFT *DataInput)  
  63. {  
  64.     int nextValue,nextM,i,k,j=0;  
  65.     complex_of_N_FFT temp;  
  66.       
  67.     nextValue=N_FFT/2;                  //变址运算，即把自然顺序变成倒位序，采用雷德算法  
  68.     nextM=N_FFT-1;  
  69.     for (i=0;i<nextM;i++)  
  70.     {  
  71.         if (i<j)                 //如果i<j,即进行变址  
  72.         {  
  73.             temp=DataInput[j];  
  74.             DataInput[j]=DataInput[i];  
  75.             DataInput[i]=temp;  
  76.         }  
  77.         k=nextValue;                //求j的下一个倒位序  
  78.         while (k<=j)             //如果k<=j,表示j的最高位为1  
  79.         {  
  80.             j=j-k;                  //把最高位变成0  
  81.             k=k/2;                  //k/2，比较次高位，依次类推，逐个比较，直到某个位为0  
  82.         }  
  83.         j=j+k;                      //把0改为1  
  84.     }                                         
  85. }                                             
  86.   
  87. //复数乘法   
  88. /*complex XX_complex(complex a, complex b) 
  89. { 
  90.     complex temp; 
  91.      
  92.     temp.real = a.real * b.real-a.imag*b.imag; 
  93.     temp.imag = b.imag*a.real + a.imag*b.real; 
  94.      
  95.     return temp; 
  96. }*/  
  97.   
  98. //FFT运算函数   
  99. void FFT(void)  
  100. {  
  101.     int L=0,B=0,J=0,K=0;  
  102.     int step=0, KB=0;  
  103.     //ElemType P=0;  
  104.     ElemType angle;  
  105.     complex_of_N_FFT W,Temp_XX;  
  106.       
  107.     ChangeSeat(data_of_N_FFT);//变址   
  108.     //CREATE_SIN_TABLE();  
  109.     for(L=1; L<=M_of_N_FFT; L++)  
  110.     {  
  111.         step = 1<<L;//2^L  
  112.         B = step>>1;//B=2^(L-1)  
  113.         for(J=0; J<B; J++)  
  114.         {  
  115.             //P = (1<<(M-L))*J;//P=2^(M-L) *J   
  116.             angle = (double)J/B;            //这里还可以优化   
  117.             W.imag =  -Sin_find(angle);     //用C++该函数课声明为inline   
  118.             W.real =   Cos_find(angle);     //用C++该函数课声明为inline   
  119.             //W.real =  cos(angle*PI);  
  120.             //W.imag = -sin(angle*PI);  
  121.             for(K=J; K<N_FFT; K=K+step)  
  122.             {  
  123.                 KB = K + B;  
  124.                 //Temp_XX = XX_complex(data[KB],W);  
  125.                 //用下面两行直接计算复数乘法，省去函数调用开销   
  126.                 Temp_XX.real = data_of_N_FFT[KB].real * W.real-data_of_N_FFT[KB].imag*W.imag;  
  127.                 Temp_XX.imag = W.imag*data_of_N_FFT[KB].real + data_of_N_FFT[KB].imag*W.real;  
  128.                   
  129.                 data_of_N_FFT[KB].real = data_of_N_FFT[K].real - Temp_XX.real;  
  130.                 data_of_N_FFT[KB].imag = data_of_N_FFT[K].imag - Temp_XX.imag;  
  131.                   
  132.                 data_of_N_FFT[K].real = data_of_N_FFT[K].real + Temp_XX.real;  
  133.                 data_of_N_FFT[K].imag = data_of_N_FFT[K].imag + Temp_XX.imag;  
  134.             }  
  135.         }  
  136.     }  
  137. }  
  138.   
  139. //IFFT运算函数   
  140. void IFFT(void)  
  141. {  
  142.     int L=0,B=0,J=0,K=0;  
  143.     int step=0, KB=0;  
  144.     //ElemType P=0;  
  145.     ElemType angle;  
  146.     complex_of_N_FFT W,Temp_XX;  
  147.       
  148.     ChangeSeat(data_of_N_FFT);//变址   
  149.     //CREATE_SIN_TABLE();  
  150.     for(L=1; L<=M_of_N_FFT; L++)  
  151.     {  
  152.         step = 1<<L;//2^L  
  153.         B = step>>1;//B=2^(L-1)  
  154.         for(J=0; J<B; J++)  
  155.         {  
  156.             //P = (1<<(M-L))*J;//P=2^(M-L) *J   
  157.             angle = (double)J/B;            //这里还可以优化   
  158.               
  159.             W.imag =   Sin_find(angle);     //用C++该函数课声明为inline   
  160.             W.real =   Cos_find(angle);     //用C++该函数课声明为inline   
  161.             //W.real =  cos(angle*PI);  
  162.             //W.imag = -sin(angle*PI);  
  163.             for(K=J; K<N_FFT; K=K+step)  
  164.             {  
  165.                 KB = K + B;  
  166.                 //Temp_XX = XX_complex(data[KB],W);  
  167.                 //用下面两行直接计算复数乘法，省去函数调用开销   
  168.                 Temp_XX.real = data_of_N_FFT[KB].real * W.real-data_of_N_FFT[KB].imag*W.imag;  
  169.                 Temp_XX.imag = W.imag*data_of_N_FFT[KB].real + data_of_N_FFT[KB].imag*W.real;  
  170.                   
  171.                 data_of_N_FFT[KB].real = data_of_N_FFT[K].real - Temp_XX.real;  
  172.                 data_of_N_FFT[KB].imag = data_of_N_FFT[K].imag - Temp_XX.imag;  
  173.                   
  174.                 data_of_N_FFT[K].real = data_of_N_FFT[K].real + Temp_XX.real;  
  175.                 data_of_N_FFT[K].imag = data_of_N_FFT[K].imag + Temp_XX.imag;  
  176.             }  
  177.         }  
  178.     }  
  179. }  
  180.   
  181. //初始化FFT程序   
  182. void Init_FFT()  
  183. {  
  184.     CREATE_SIN_TABLE();             //创建正弦函数表   
  185. }  
  186.   
  187. //结束 FFT运算，释放相关内存   
  188. void Close_FFT(void)  
  189. {  
  190.       
  191. }  
  192.   
  193. #define SYRFFT_6_55_H  
  194. #endif  

 

[cpp] view
plaincopy

  1. /******************************************************************************* 
  2. ** 程序名称：快速傅里叶变换(FFT)  
  3. ** 程序描述：本程序实现快速傅里叶变换  
  4. ** 性能提升：修正了IFFT的bug,变量重新命名,使用宏定义改变N大小  
  5. ** 程序版本：V6.55 
  6. ** 程序作者：宋元瑞  
  7. ** 最后修改：2011年5月22日  
  8. *******************************************************************************/  
  9. #include "syrFFT_6_55.h"   
  10.   
  11. //产生模拟原始数据输入 ,在实际应用时替换为AD采样数据   
  12. void InputData(void)  
  13. {  
  14.     int i;  
  15.     for(i=0; i<N_FFT; i++)//制造输入序列   
  16.     {  
  17.         data_of_N_FFT[i].real = sin(2*PI*i/N_FFT);  //输入采样数据   
  18.         data_of_N_FFT[i].imag = 0;  
  19.     }  
  20. }  
  21.   
  22. //主函数 ，示例如何调用FFT运算   
  23. int main(int argc, char *argv[])  
  24. {  
  25.     int i = 0;  
  26.       
  27.     Init_FFT();         //①初始化各项参数，此函数只需执行一次 ; 修改FFT的点数去头文件的宏定义处修改   
  28.       
  29.     InputData();        //②输入原始数据 ,此处在实际应用时替换为AD采样数据   
  30.     FFT();              //③进行 FFT计算   
  31.     //for(i=0; i<N_FFT; i++)//④输出FFT频谱结果  sqrt(data_of_N_FFT[i].real*data_of_N_FFT[i].real+data_of_N_FFT[i].imag*data_of_N_FFT[i].imag)  
  32.         //{printf("%f ",FFT_RESULT(i));}  
  33.       
  34.     IFFT();//进行 IFFT计算   
  35.     //for(i=0; i<N_FFT; i++)//(可选步骤)⑤输出 IFFT结果 ，滤波时会用到；暂时不用   
  36.         //{printf("%f ",IFFT_RESULT(i));}  
  37.       
  38.     Close_FFT();    //结束 FFT运算，此版本此句无用，可不写   
  39.       
  40.     return 0;  
  41. }  
返回
【上篇】
【下篇】

作者:

抱歉!评论已关闭.

返回首页

Copyright © 2013-2018 学步园  保留所有权利.
软文销售 QQ客服：2265327166 点击这里给我发消息（其他合作也可洽谈）
必威体育
必威电竞