21ic电子技术开发论坛 › 单片机与嵌入式系统 › GigaDevice GD32 MCU › FFT快速傅里叶变换算法，51单片机都能用 ...
返回列表 发新帖我要提问

FFT快速傅里叶变换算法，51单片机都能用

[复制链接]
1962|10
 楼主| saservice 发表于 2023-9-28 06:28 | 显示全部楼层 |阅读模式
51单片机, fft, 傅里叶变换, 单片机, 快速傅里叶


复制
  1. /*********************************************************************
    

  2. 快速傅里叶变换C程序包
    

  3. 函数简介：此程序包是通用的快速傅里叶变换C语言函数，移植性强，以下部分不依
    

  4. 赖硬件。此程序包采用联合体的形式表示一个复数，输入为自然顺序的复
    

  5. 数（输入实数是可令复数虚部为0），输出为经过FFT变换的自然顺序的
    

  6. 复数.此程序包可在初始化时调用create_sin_tab()函数创建正弦函数表，
    

  7. 以后的可采用查表法计算耗时较多的sin和cos运算，加快可计算速度.与
    

  8. Ver1.1版相比较，Ver1.2版在创建正弦表时只建立了1/4个正弦波的采样值，
    

  9. 相比之下节省了FFT_N/4个存储空间
    

  10. 使用说明：使用此函数只需更改宏定义FFT_N的值即可实现点数的改变，FFT_N的
    

  11. 应该为2的N次方，不满足此条件时应在后面补0。若使用查表法计算sin值和
    

  12. cos值，应在调用FFT函数前调用create_sin_tab()函数创建正弦表
    

  13. 函数调用：FFT(Compx);
    

  14. 作    者：吉帅虎
    

  15. 时    间：2010-2-20
    

  16. 版    本：Ver1.2
    

  17. 参考文献：
    

  18. **********************************************************************/
    

  19. #include <math.h>
    

  20. #include "FFT.h"
    


  21. 

  22. struct compx Compx[FFT_N] = {0};                //FFT输入和输出：从Compx[0]开始存放，根据大小自己定义
    

  23. double SIN_TAB[FFT_N / 4 + 1];                        //定义正弦表的存放空间
    


  24. 

  25. /*******************************************************************
    

  26. 函数原型：struct compx EE(struct compx b1,struct compx b2)
    

  27. 函数功能：对两个复数进行乘法运算
    

  28. 输入参数：两个以联合体定义的复数a,b
    

  29. 输出参数：a和b的乘积，以联合体的形式输出
    

  30. *******************************************************************/
    

  31. struct compx EE(struct compx a, struct compx b)
    

  32. {
    

  33.         struct compx c;
    

  34.         c.real = a.real*b.real - a.imag*b.imag;
    

  35.         c.imag = a.real*b.imag + a.imag*b.real;
    

  36.         return(c);
    

  37. }
    


  38. 

  39. /******************************************************************
    

  40. 函数原型：void create_sin_tab(double *sin_t)
    

  41. 函数功能：创建一个正弦采样表，采样点数与福利叶变换点数相同
    

  42. 输入参数：*sin_t存放正弦表的数组指针
    

  43. 输出参数：无
    

  44. ******************************************************************/
    

  45. void create_sin_tab(double *sin_t)
    

  46. {
    

  47.         int i;
    

  48.         for (i = 0; i <= FFT_N / 4; i++)
    

  49.                 sin_t[i] = sin(2 * PI*i / FFT_N);
    

  50. }
    

  51. /******************************************************************
    

  52. 函数原型：void sin_tab(double pi)
    

  53. 函数功能：采用查表的方法计算一个数的正弦值
    

  54. 输入参数：pi 所要计算正弦值弧度值，范围0--2*PI，不满足时需要转换
    

  55. 输出参数：输入值pi的正弦值
    

  56. ******************************************************************/
    

  57. double sin_tab(double pi)
    

  58. {
    

  59.         int n;
    

  60.         double a = 0;
    

  61.         n = (int)(pi*FFT_N / 2 / PI);
    


  62. 

  63.         if (n >= 0 && n <= FFT_N / 4)
    

  64.                 a = SIN_TAB[n];
    

  65.         else if (n>FFT_N / 4 && n<FFT_N / 2)
    

  66.         {
    

  67.                 n -= FFT_N / 4;
    

  68.                 a = SIN_TAB[FFT_N / 4 - n];
    

  69.         }
    

  70.         else if (n >= FFT_N / 2 && n<3 * FFT_N / 4)
    

  71.         {
    

  72.                 n -= FFT_N / 2;
    

  73.                 a = -SIN_TAB[n];
    

  74.         }
    

  75.         else if (n >= 3 * FFT_N / 4 && n<3 * FFT_N)
    

  76.         {
    

  77.                 n = FFT_N - n;
    

  78.                 a = -SIN_TAB[n];
    

  79.         }
    


  80. 

  81.         return a;
    

  82. }
    

  83. /******************************************************************
    

  84. 函数原型：void cos_tab(double pi)
    

  85. 函数功能：采用查表的方法计算一个数的余弦值
    

  86. 输入参数：pi 所要计算余弦值弧度值，范围0--2*PI，不满足时需要转换
    

  87. 输出参数：输入值pi的余弦值
    

  88. ******************************************************************/
    

  89. double cos_tab(double pi)
    

  90. {
    

  91.         double a, pi2;
    

  92.         pi2 = pi + PI / 2;
    

  93.         if (pi2>2 * PI)
    

  94.                 pi2 -= 2 * PI;
    

  95.         a = sin_tab(pi2);
    

  96.         return a;
    

  97. }
    

  98. /*****************************************************************
    

  99. 函数原型：void FFT(struct compx *xin)
    

  100. 函数功能：对输入的复数组进行快速傅里叶变换（FFT）
    

  101. 输入参数：*xin复数结构体组的首地址指针，struct型
    

  102. 输出参数：无
    

  103. *****************************************************************/
    

  104. void FFT(struct compx *xin)
    

  105. {
    

  106.         register int f, m, nv2, nm1, i, k, l, j = 0;
    

  107.         struct compx u, w, t;
    


  108. 

  109.         nv2 = FFT_N / 2;                                        //变址运算，即把自然顺序变成倒位序，采用雷德算法
    

  110.         nm1 = FFT_N - 1;
    

  111.         for (i = 0; i < nm1; ++i)
    

  112.         {
    

  113.                 if (i < j)                                                //如果i<j,即进行变址
    

  114.                 {
    

  115.                         t = xin[j];
    

  116.                         xin[j] = xin[i];
    

  117.                         xin[i] = t;
    

  118.                 }
    

  119.                 k = nv2;                                                //求j的下一个倒位序
    

  120.                 while (k <= j)                                        //如果k<=j,表示j的最高位为1
    

  121.                 {
    

  122.                         j = j - k;                                        //把最高位变成0
    

  123.                         k = k / 2;                                        //k/2，比较次高位，依次类推，逐个比较，直到某个位为0
    

  124.                 }
    

  125.                 j = j + k;                                                //把0改为1
    

  126.         }
    


  127. 

  128.         {
    

  129.                 int le, lei, ip;                                //FFT运算核，使用蝶形运算完成FFT运算
    

  130.                 f = FFT_N;
    

  131.                 for (l = 1; (f = f / 2) != 1; ++l);                                //计算l的值，即计算蝶形级数
    

  132.                 for (m = 1; m <= l; m++)                                                // 控制蝶形结级数
    

  133.                 {   
    

  134.                         //m表示第m级蝶形，l为蝶形级总数l=log（2）N
    

  135.                         le = 2 << (m - 1);                                                        //le蝶形结距离，即第m级蝶形的蝶形结相距le点
    

  136.                         lei = le / 2;                               //同一蝶形结中参加运算的两点的距离
    

  137.                         u.real = 1.0;                                                                //u为蝶形结运算系数，初始值为1
    

  138.                         u.imag = 0.0;
    

  139.                         w.real = cos_tab(PI / lei);                                        //w为系数商，即当前系数与前一个系数的商
    

  140.                         w.imag = -sin_tab(PI / lei);
    

  141.                         for (j = 0; j <= lei - 1; j++)                                //控制计算不同种蝶形结，即计算系数不同的蝶形结
    

  142.                         {
    

  143.                                 for (i = j; i <= FFT_N - 1; i = i + le)        //控制同一蝶形结运算，即计算系数相同蝶形结
    

  144.                                 {
    

  145.                                         ip = i + lei;                                                //i，ip分别表示参加蝶形运算的两个节点
    

  146.                                         t = EE(xin[ip], u);                                        //蝶形运算，详见公式
    

  147.                                         xin[ip].real = xin[i].real - t.real;
    

  148.                                         xin[ip].imag = xin[i].imag - t.imag;
    

  149.                                         xin[i].real = xin[i].real + t.real;
    

  150.                                         xin[i].imag = xin[i].imag + t.imag;
    

  151.                                 }
    

  152.                                 u = EE(u, w);                                                        //改变系数，进行下一个蝶形运算
    

  153.                         }
    

  154.                 }
    

  155.         }
    

  156. }
    


  157. 


  158. 

  159. /*****************************************************************
    

  160. 函数原型：void Get_Result(struct compx *xin, double sample_frequency)
    

  161. 函数功能：求变换后结果的模值，存入复数的实部部分，频率存入复数的虚数部分，有效数据为前FFT_N/2个数
    

  162. 输入参数：*xin复数结构体组的首地址指针，struct型, sample_frequency: 采样频率
    

  163. 输出参数：无
    

  164. *****************************************************************/
    

  165. void Get_Result(struct compx *xin, double sample_frequency)
    

  166. {
    

  167.         int i = 0;
    

  168.         for (i = 0; i < FFT_N / 2; ++i)
    

  169.         {          
    

  170.                 //求变换后结果的模值，存入复数的实部部分
    

  171.                 xin[i].real = sqrt(xin[i].real*xin[i].real + xin[i].imag*xin[i].imag) / (FFT_N >> (i != 0));
    

  172.                 xin[i].imag = i * sample_frequency / FFT_N;
    

  173.         }
    

  174. }
    


  175. 

  176. /*****************************************************************
    

  177. 函数原型：void Refresh_Data(struct compx *xin, double wave_data)
    

  178. 函数功能：更新数据
    

  179. 输入参数：*xin复数结构体组的首地址指针, struct型, id: 标号, wave_data: 一个点的值
    

  180. 输出参数：无
    

  181. *****************************************************************/
    

  182. void Refresh_Data(struct compx *xin, int id, double wave_data)
    

  183. {
    

  184.         xin[id].real = wave_data;
    

  185.         xin[id].imag = 0;
    

  186. }
    


  187.

复制
  1. #ifndef FFT_H
    

  2. #define FFT_H
    


  3. 

  4. #define FFT_N 16                                        //定义傅里叶变换的点数
    

  5. #define PI 3.14159265358979323846264338327950288419717  //定义圆周率值
    


  6. 

  7. struct compx { double real, imag; };                    //定义一个复数结构
    


  8. 

  9. extern struct compx Compx[];                                                        //FFT输入和输出：从Compx[0]开始存放，根据大小自己定义
    

  10. extern double SIN_TAB[];                                                                //正弦信号表
    


  11. 

  12. extern void Refresh_Data(struct compx *xin, int id, double wave_data);
    

  13. extern void create_sin_tab(double *sin_t);
    

  14. extern void FFT(struct compx *xin);
    

  15. extern void Get_Result(struct compx *xin, double sample_frequency);
    


  16. 

  17. #endif
    


  18.

复制
  1.  #define Sample_Frequency 800                        //采样频率 800Hz
    

  2.         #define Frequency 100                                        //测试信号 100Hz
    


  3. 

  4.         for (i = 0; i < FFT_N; ++i)                                //使用Refresh_Data(Compx, i, wave_data)函数填入数据，这里建立一个100Hz的方波测试信号
    

  5.         { 
    

  6.                 if (sin(2 * PI * Frequency * i / Sample_Frequency) > 0)
    

  7.                         Refresh_Data(Compx, i, 1);
    

  8.                 else if (sin(2 * PI * Frequency * i / Sample_Frequency) < 0)
    

  9.                         Refresh_Data(Compx, i, -1);
    

  10.                 else
    

  11.                         Refresh_Data(Compx, i, 0);
    

  12.         }                                                                                
    


  13. 

  14.         create_sin_tab(SIN_TAB);                                //建立正弦信号表, 之后将用查表法计算正弦值，加速计算
    

  15.         FFT(Compx);                                                                //快速傅里叶变换
    

  16.         Get_Result(Compx, Sample_Frequency);        //求变换后结果的模值，存入复数的实部部分，频率存入复数的虚数部分，有效数据为前FFT_N/2个数


回复

收藏 举报

相关帖子
AloneKaven 发表于 2023-9-29 22:19 来自手机 | 显示全部楼层
宏定义改点数，太好了！
回复

收藏 举报

AloneKaven 发表于 2023-9-29 22:19 来自手机 | 显示全部楼层
上次有个64点改128直接给我愁死
回复

收藏 举报

tpgf 发表于 2023-10-13 10:09 | 显示全部楼层
快速傅里叶变换即利用计算机计算离散傅里叶变换(DFT)的高效、快速计算方法的统称，简称FFT
回复

收藏 举报

guanjiaer 发表于 2023-10-13 13:10 | 显示全部楼层
采用这种算法能使计算机计算离散傅里叶变换所需要的乘法次数大为减少
回复

收藏 举报

keaibukelian 发表于 2023-10-13 13:26 | 显示全部楼层
快速傅氏变换(FFT)，是离散傅氏变换的快速算法，它是根据离散傅氏变换的奇、偶、虚、实等特性，对离散傅立叶变换的算法进行改进获得的
回复

收藏 举报

heimaojingzhang 发表于 2023-10-13 14:25 | 显示全部楼层
被变换的抽样点数N越多，FFT算法计算量的节省就越显著
回复

收藏 举报

paotangsan 发表于 2023-10-14 09:03 | 显示全部楼层
它是根据离散傅氏变换的奇、偶、虚、实等特性,对离散傅立叶变换的算法进行改进获得的
回复

收藏 举报

renzheshengui 发表于 2023-10-14 09:48 | 显示全部楼层
FFT的基本思想是把原始的N点序列，依次分解成一系列的短序列。充分利用DFT计算式中指数因子 所具有的对称性质和周期性质，进而求出这些短序列相应的DFT并进行适当组合，达到删除重复计算，减少乘法运算和简化结构的目的
回复

收藏 举报

ql2000 发表于 2023-10-14 14:19 | 显示全部楼层
谢谢楼主分享
回复

收藏 举报

chenjun89 发表于 2023-10-17 08:09 来自手机 | 显示全部楼层
用肯定能用，只是效率有点低而已。
回复

收藏 举报

返回列表 发新帖
高级模式
B Color Image Link Quote Code Smilies
您需要登录后才可以回帖 登录 | 注册

本版积分规则

saservice

43

主题

1565

帖子

2

粉丝
快速回复 在线客服 返回列表 返回顶部