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

最新推荐文章于 2021-11-11 04:45:39 发布
最新推荐文章于 2021-11-11 04:45:39 发布
阅读量4.3k 收藏 9
点赞数

http://blog.csdn.net/syrchina/article/details/6670517

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

[cpp]  view plain copy
  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 plain copy
  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 plain copy
  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 plain copy
  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 plain copy
  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 plain copy
  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 plain copy
  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 plain copy
  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 plain copy
  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. }  



 

anakin4
关注
MATLAB、FPGA、STM32中调用FFT计算频率、幅值及相位差
贾saisai的博客
09-13 1748
最近在学习如何在STM32中调用FFT。
Qt和MDK-ARM平台下的fft
linux驱动开发大纲
03-15 600
QT平台 使用的库：fftw-3.3.5-dll32 .pro文件加入如下内容 LIBS += -L$$PWD/./ -llibfftw3-3 INCLUDEPATH += $$PWD/. DEPENDPATH += $$PWD/. 使用 #include "fftw3.h" #define N 256 fftw_complex *in, *out;//声明复数类型的输入数据in1_c和FFT变换的结果out1_c resolver::resolver(QObject *parent) :
ARM的DSP增强型扩展
Java-Jinguo
02-12 696
　　现在，很多新兴的应用领域涌现了许多新的算法标准，这些算法对于处理器提出了更高的性能和控制要求。信号处理需要处理器提供高峰值性能，但这部分在整个算法中的比例有减少的趋势。对于开发包含高性能算法的大规模应用的设计团队来说，有很多种可供选择的方案。市场竞争的压力使得选择一个对高性能有充分保障的处理器平台来实现高性能的功能变得非常重要。但是选择超出需求性能很多的处理器平台对于系统成本和电源消耗也是有着...
基于ARM的fft算法实现
07-13
该程序是在arm板子上对512点的fft算法实现
用于ARM上的FFT与IFFT源代码-C语言
weixin_34162401的博客
07-17 380
/********************************************************************************* 程序名称：快速傅里叶变换(FFT) ** 程序描述：本程序实现快速傅里叶变换 ** 程序作者：宋元瑞 ** 最后修改：2011年4月5日 ************************************************...
C语言实现的FFT与IFFT源代码，不依赖特定平台
乙酸氧铍的博客
11-11 7240
目录源码FFT.cFFT.h使用方法初始化输入数据FFT 快速傅里叶变换解算FFT结果使用python绘制FFT波形IFFT 快速傅里叶逆变换解算IFFT结果 Windows 10 20H2 Visual Studio 2015 Python 3.8.12 (default, Oct 12 2021, 03:01:40) [MSC v.1916 64 bit (AMD64)] :: Anaconda, Inc. on win32 小改自用于ARM上的FFT与IFFT源代码（C语言，不依赖特定平台）——
FFT.rar_C FFT IFFT_c语言实现FFT_fft_fft ifft c_fft ifft c语言
09-24
标题中的"FFT.rar_C FFT IFFT_c语言实现FFT_fft_fft ifft c_fft ifft c语言"指的是一个关于快速傅里叶变换（Fast Fourier Transform, FFT）和逆快速傅里叶变换（Inverse Fast Fourier Transform, IFFT）的C语言实现...
fft_ifft.rar_FFT反变换_fft ifft c++_fft的反变换_ifft源代码_傅里叶
09-23
在实际应用中，需要注意的是，FFT和IFFT对于输入序列的长度有特定的要求，即必须是2的幂次。如果原始序列长度不是2的幂次，通常需要通过填充零或者截断来调整长度。此外，对于实数序列，可以使用更高效的实数FFT算法...
fft.rar_fft 定点_fft ifft c语言_fft 频率_ifft
09-19
**标题与描述解析** ...综上所述，这个压缩包文件很可能包含了用C语言编写的定点FFT和IFFT算法的源代码，以及可能的频率抽选功能。这些工具对于进行数字信号处理、滤波、频谱分析等任务非常有用。
fft.rar_C语言 FFT_c语言 fft_fft c++_fft c语言_fft_ifft
09-14
标题中的"fft.rar"指的是一个使用C语言实现快速傅里叶变换(FFT)的代码库，其中包含了C语言版本的FFT以及逆快速傅里叶变换(IFFT)。这个压缩包可能是一个开发项目或者学习资源，提供了详细的注释，便于理解和应用。 ...
快速傅里叶变换（fft\ifft）C语言实现
06-16
根据网上已有的资源，经过自己的整理的快速傅里叶变换，希望能够为信号处理，信号的时域、频域转换需求的开发，研究人员带来方便。 敬请批评指正！
可以用于arm的fft算法
05-19
可以用于arm的fft算法，对arm做了优化处理，速度很快、
嵌入式系统/ARM技术中的嵌入式系统中FFT算法研究
12-08
摘　要：首先分析实数FFT算法的推导过程，然后给出一种具体实现FFT算法的C语言程序，可以直接应用于需要FFT运算的单片机或DSP等嵌入式系统中。 关键词：嵌入式系统 FFT算法 单片机 DSP 目前国内有关数字信号处理的教材在讲解快速傅里叶变换（FFT）时，都是以复数FFT为重点，实数FFT算法都是一笔带过，书中给出的具体实现程序多为BASIC或FORTRAN程序并且多数不能真正运行。鉴于目前在许多嵌入式系统中要用到FFT运算，如以DSP为核心的交流采样系统、频谱分析、相关分析等。本人结合自己的实际开发经验，研究了实数的FFT算法并给出具体的C语言函数，读者可以直接应用于自己的系统中
FFTW3 windows和arm 版本库
最新发布
09-05
使用方式和编译方法可参考博客：https://blog.csdn.net/lmgandxka/article/details/136811326 也可直接下载该资源包使用. 描述：FFTW是一个 C 子程序库，用于计算一维或多维、任意输入大小、实数和复数数据（以及偶数/奇数数据，即离散余弦/正弦变换或 DCT/DST）的离散傅里叶变换 (DFT)。 也就是说它可移植到任何具有 C 编译器的平台。 特点： Speed. (Supports SSE/SSE2/Altivec, since version 3.0. Version 3.3.1 supports AVX and ARM Neon.) Both one-dimensional and multi-dimensional transforms. Arbitrary-size transforms. (Sizes with small prime factors are best, but FFTW uses O(N log N) algorithms even for prime sizes.) ………………
arm 查表法FFT 256点和512点
04-01
单片机 arm中用到的快速傅立叶函数，查表法提高效率，512点可以满足绝大部分的要求。
fft_ifft.rar_C语言 DSP IFFT_FFT_IFFT.rar_ifft dsp_混合基fft_混合基·fft i
09-24
标题中的“fft_ifft.rar”和“C语言 DSP IFFT_FFT_IFFT.rar_ifft dsp_混合基fft_混合基·fft i”暗示了这个压缩包包含的是与快速傅里叶变换（FFT）和逆快速傅里叶变换（IFFT）相关的C语言源代码，特别是针对数字信号...
用于ARM上的FFT与IFFT源代码（C语言，不依赖特定平台）（转）
weixin_33724570的博客
08-06 208
源：用于ARM上的FFT与IFFT源代码（C语言，不依赖特定平台） 代码在2011年全国电子大赛结束后（2011年9月3日）发布，多个版本，注释详细。 /******************************************************************************* ** 程序名称：快速傅里叶变换(FFT) ** 程序描述：本程序实现快速...
评论 3
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值