[应用相关] C语言实现fft理论基础与工程应用的实例分析

[复制链接]
2291|11
 楼主| xuanhuanzi 发表于 2018-1-31 19:02 | 显示全部楼层 |阅读模式
一、理论分析

快速傅里叶变换(Fast Fourier Transform)是离散傅里叶变换的一种快速算法,简称FFT,通过FFT可以将一个信号从时域变换到频域。

模拟信号经过A/D转换变为数字信号的过程称为采样。为保证采样后信号的频谱形状不失真,采样频率必须大于信号中最高频率成分的2倍,这称之为采样定理。

假设采样频率为fs,采样点数为N,那么FFT结果就是一个N点的复数,每一个点就对应着一个频率点,某一点n(n从1开始)表示的频率为:fn=(n-1)*fs/N。

举例说明:用1kHz的采样频率采样128点,则FFT结果的128个数据即对应的频率点分别是0,1k/128,2k/128,3k/128,…,127k/128 Hz。

这个频率点的幅值为:该点复数的模值除以N/2(n=1时是直流分量,其幅值是该点的模值除以N)。

      下面先来简要分析下封装好的FFT的C程序包
  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(s);
  14. 作    者:吉帅虎
  15. 时    间:2010-2-20
  16. 版    本:Ver1.2
  17. 参考文献:
  18. **********************************************************************/
  19. #include <math.h>
  20. #include "fft.h"

  21. float *SIN_TAB;//定义正弦表的存放空间
  22. int FFT_N = 1024;//定义采样点大小
  23. /*******************************************************************
  24. 函数原型:struct compx EE(struct compx b1,struct compx b2)
  25. 函数功能:对两个复数进行乘法运算
  26. 输入参数:两个以联合体定义的复数a,b
  27. 输出参数:a和b的乘积,以联合体的形式输出
  28. *******************************************************************/
  29. struct compx EE(struct compx a,struct compx b)
  30. {
  31. struct compx c;
  32. c.real=a.real*b.real-a.imag*b.imag;
  33. c.imag=a.real*b.imag+a.imag*b.real;
  34. return(c);
  35. }

  36. /******************************************************************
  37. 函数原型:void create_sin_tab(float *sin_t,int PointNum)
  38. 函数功能:创建一个正弦采样表,采样点数与福利叶变换点数相同
  39. 输入参数:*sin_t存放正弦表的数组指针,PointNum采样点数
  40. 输出参数:无
  41. ******************************************************************/
  42. void create_sin_tab(float *sin_t,int PointNum)
  43. {
  44.   int i;
  45.   SIN_TAB=sin_t;
  46.   FFT_N=PointNum;
  47.   for(i=0;i<=FFT_N/4;i++)
  48.     SIN_TAB[i]=sin(2*PI*i/FFT_N);
  49. }
  50. /******************************************************************
  51. 函数原型:void sin_tab(float pi)
  52. 函数功能:采用查表的方法计算一个数的正弦值
  53. 输入参数:pi 所要计算正弦值弧度值,范围0--2*PI,不满足时需要转换
  54. 输出参数:输入值pi的正弦值
  55. ******************************************************************/
  56. float sin_tab(float pi)
  57. {
  58.   int n=0;
  59.   float a=0;
  60.    n=(int)(pi*FFT_N/2/PI);

  61.   if(n>=0&&n<=FFT_N/4)
  62.     a=SIN_TAB[n];
  63.   else if(n>FFT_N/4&&n<FFT_N/2)
  64.     {
  65.      n-=FFT_N/4;
  66.      a=SIN_TAB[FFT_N/4-n];
  67.     }
  68.   else if(n>=FFT_N/2&&n<3*FFT_N/4)
  69.     {
  70.      n-=FFT_N/2;
  71.      a=-SIN_TAB[n];
  72.    }
  73.   else if(n>=3*FFT_N/4&&n<3*FFT_N)
  74.     {
  75.      n=FFT_N-n;
  76.      a=-SIN_TAB[n];
  77.    }

  78.   return a;
  79. }
  80. /******************************************************************
  81. 函数原型:void cos_tab(float pi)
  82. 函数功能:采用查表的方法计算一个数的余弦值
  83. 输入参数:pi 所要计算余弦值弧度值,范围0--2*PI,不满足时需要转换
  84. 输出参数:输入值pi的余弦值
  85. ******************************************************************/
  86. float cos_tab(float pi)
  87. {
  88.    float a,pi2;
  89.    pi2=pi+PI/2;
  90.    if(pi2>2*PI)
  91.      pi2-=2*PI;
  92.    a=sin_tab(pi2);
  93.    return a;
  94. }
  95. /*****************************************************************
  96. 函数原型:void FFT(struct compx *xin)
  97. 函数功能:对输入的复数组进行快速傅里叶变换(FFT)
  98. 输入参数:*xin复数结构体组的首地址指针,struct型
  99. 输出参数:无
  100. *****************************************************************/
  101. void FFT(struct compx *xin)
  102. {
  103.   int f,m,i,k,l,j=0;
  104.   register int nv2,nm1;
  105.   struct compx u,w,t;

  106.    nv2=FFT_N/2;                  //变址运算,即把自然顺序变成倒位序,采用雷德算法
  107.    nm1=FFT_N-1;
  108.    for(i=0;i<nm1;i++)
  109.    {
  110.     if(i<j)                    //如果i<j,即进行变址
  111.      {
  112.       t=xin[j];
  113.       xin[j]=xin[i];
  114.       xin[i]=t;
  115.      }
  116.     k=nv2;                    //求j的下一个倒位序
  117.     while(k<=j)               //如果k<=j,表示j的最高位为1
  118.      {
  119.       j=j-k;                 //把最高位变成0
  120.       k=k/2;                 //k/2,比较次高位,依次类推,逐个比较,直到某个位为0
  121.      }
  122.    j=j+k;                   //把0改为1
  123.   }

  124.   {
  125.    int le,lei,ip;                            //FFT运算核,使用蝶形运算完成FFT运算
  126.     f=FFT_N;
  127.    for(l=1;(f=f/2)!=1;l++)                  //计算l的值,即计算蝶形级数
  128.            ;
  129.   for(m=1;m<=l;m++)                         // 控制蝶形结级数
  130.    {                                        //m表示第m级蝶形,l为蝶形级总数l=log(2)N
  131.     le=2<<(m-1);                            //le蝶形结距离,即第m级蝶形的蝶形结相距le点
  132.     lei=le/2;                               //同一蝶形结中参加运算的两点的距离
  133.     u.real=1.0;                             //u为蝶形结运算系数,初始值为1
  134.     u.imag=0.0;
  135.     //w.real=cos(PI/lei);                  //不适用查表法计算sin值和cos值
  136.     // w.imag=-sin(PI/lei);
  137.     w.real=cos_tab(PI/lei);                //w为系数商,即当前系数与前一个系数的商
  138.     w.imag=-sin_tab(PI/lei);
  139.     for(j=0;j<=lei-1;j++)                  //控制计算不同种蝶形结,即计算系数不同的蝶形结
  140.      {
  141.       for(i=j;i<=FFT_N-1;i=i+le)           //控制同一蝶形结运算,即计算系数相同蝶形结
  142.        {
  143.         ip=i+lei;                          //i,ip分别表示参加蝶形运算的两个节点
  144.         t=EE(xin[ip],u);                   //蝶形运算,详见公式
  145.         xin[ip].real=xin[i].real-t.real;
  146.         xin[ip].imag=xin[i].imag-t.imag;
  147.         xin[i].real=xin[i].real+t.real;
  148.         xin[i].imag=xin[i].imag+t.imag;
  149.        }
  150.       u=EE(u,w);                          //改变系数,进行下一个蝶形运算
  151.      }
  152.    }
  153.   }
  154. }

  155. fft.c
  1. #ifndef FFT_H
  2. #define FFT_H
  3.                                                  //定义福利叶变换的点数
  4. #define PI 3.1415926535897932384626433832795028841971               //定义圆周率值

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

  6. extern void create_sin_tab(float *sin_t,int PointNum);
  7. extern void FFT(struct compx *xin);

  8. #endif // FFT_H

  9. fft.h


 楼主| xuanhuanzi 发表于 2018-1-31 19:03 | 显示全部楼层
   程序主要分为两个部分,第一部分主要采用雷德算法来实现倒位序,第二部分主要是利用已经生成好的旋转矩阵做蝶形变换。(程序来源于《数字信号处理》清华大学出版社,程佩青,按时间抽选的基-2 FFT蝶形图)

1.第一个要解决的问题就是码位倒序,以雷德算法为例。

下面假如使用A[I]存的是顺序位序,而B[J]存的是倒位序。I<J的时候需要变序,I>J的时候就不用,不然就白忙活了。

例如   N = 8 的时候,
倒位序 顺序          二进制表示      倒位序 顺序
0 0                                       000          000
4 1                                       100          001
2 2                                       010          010           
6 3                                       110          011
1 4                                        001         100
5 5                                        101         101
3 6                                        011         110
7 7                                        111          111

     由上面的表可以看出,按自然顺序排列的二进制数,其下面一个数总是比其上面一个数大1,即下面一个数是上面一个数在最低位加1并向高位进位而得到的。而倒位序二进制数的下面一个数是上面一个数在最高位加1并由高位向低位进位而得到。
I、J都是从0开始,若已知某个倒位序J,要求下一个倒位序数,则应先判断J的最高位是否为0,这可与k=N/2相比较,因为N/2总是等于100..的。如果k>J,则J的最高位为0,只要把该位变为1(J与k=N/2相加即可),就得到下一个倒位序数;如果K<=J,则J的最高位为1,可将最高位变为0(J与k=N/2相减即可)。然后还需判断次高位,这可与k=N\4相比较,若次高位为0,则需将它变为1(加N\4即可)其他位不变,既得到下一个倒位序数;若次高位是1,则需将它也变为0。然后再判断下一位。具体代码实现参见代码fft.c。
2.第二个要解决的问题就是蝶形运算
241127438878453.jpg
①第1级(第1列)每个蝶形的两节点“距离”为1,第2级每个蝶形的两节点“距离”为2,第3级每个蝶形的两节点“距离”为4,第4级每个蝶形的两节点“距离”为8。由此推得,第m级蝶形运算,每个蝶形的两节点“距离”L=2m-1。

②对于16点的FFT,第1级有16组蝶形,每组有1个蝶形;第2级有4组蝶形,每组有2个蝶形;第3级有2组蝶形,每组有4个蝶形;第4级有1组蝶形,每组有8个蝶形。由此可推出,对于N点的FFT,第m级有N/2L组蝶形,每组有L=2m-1个蝶形。

③旋转因子的确定

为提高FFT的运算速度,我们可以事先建立一个旋转因子数组,然后通过查表法来实现

complex code WN[N](旋转因子数组)
为节省CPU计算时间,旋转因子采用查表处理
根据实际FFT的点数N,该表数据需自行修改
以下结果通过Excel自动生成
WN[k].real=cos(2*PI/N*k)
WN[k].img=-sin(2*PI/N*k)


 楼主| xuanhuanzi 发表于 2018-1-31 19:04 | 显示全部楼层
二、工程应用

     为了找到基波位置,我们对第一次平均值做fft变换。第一次平均值即反映了ADC采样的值得整体趋势,只是在幅值上有所降低。ADC的定时器触发的采样频率为20KHz,40个采样点一次平均后的采样频率为500Hz,故fft的采样频率为500Hz(也即fft变换后最大频率不会达到500Hz)。

fft的采样点为1024点。故fft的分辨率约为0.49Hz。下面对具体程序和实际情况分析。
  1. void FFT_handle(u16 Input_voltage)
  2. {
  3.    
  4.     int i = 0;
  5.    
  6.     s[t].real = Input_voltage;
  7.     s[t].imag = 0;
  8.     t++;
  9.     if(t >= NPT)
  10.     {
  11.         t = 0;
  12.         value_ACsignalMAX = 0;
  13.         FFT(s);
  14.         for(i=0;i<NPT / 2;i++)
  15.         {                               //求变换后结果的模值,存入复数的实部部分
  16.         //    s[i].real=(u16)(sqrt(s[i].real*s[i].real+s[i].imag*s[i].imag)/(i==0?NPT:NPT/2));
  17.             table[i]=(u16)(sqrt(s[i].real*s[i].real+s[i].imag*s[i].imag)/(i==0?NPT:NPT/2));

  18.             if(i == 0)
  19.             {
  20.                 value_DCsignal = table[i];
  21.             }
  22.             else if(0 < i && i < NPT / 2)
  23.             {
  24.                 if(table[i] > value_ACsignalMAX)
  25.                 {
  26.                     if(i < 15) // 滤除高频干扰
  27.                     {
  28.                         value_ACsignalMAX = table[i];
  29.                         MAX_palce = i;
  30.                     }
  31.                     else{}
  32.                 }
  33.             }
  34.             else{}   
  35.    
  36.                                                                 jww=value_DCsignal;
  37.             jwww=value_ACsignalMAX;
  38.         }

  39. fft handle.c
首先要求出fft变化后的各频率点的模值,i=0对应直流分量,其次找到最大频率点(即基波),但程序界定i<15,故滤除高于15*1.49=7.4Hz的分量。
在debug全速运行时,程序稳定状态下测得Max i=9,故对应4.4Hz。

用信号发生器测得被ADC采样波形和将其做fft变换。可以看出基波频率为4.148Hz,而示波器fft变换后最大模值点为4.19Hz。

由此得出结论:模拟测试4.14Hz与程序结果4.4Hz(i=9)基本一致,误差来源于fft采样分辨率0.49Hz。 1.gif 2.gif

ylslib 发表于 2018-1-31 21:43 | 显示全部楼层
FFT一般用在哪些场合呢?声音合成吗?
wahahaheihei 发表于 2018-1-31 23:28 来自手机 | 显示全部楼层
这种消耗资源多吗
 楼主| xuanhuanzi 发表于 2018-2-6 18:00 | 显示全部楼层
 楼主| xuanhuanzi 发表于 2018-2-6 18:00 | 显示全部楼层
ylslib 发表于 2018-1-31 21:43
FFT一般用在哪些场合呢?声音合成吗?

滤波,声音可以变声,或者调音。
wahahaheihei 发表于 2018-2-7 18:51 | 显示全部楼层
看着其实很深奥的。
plsbackup 发表于 2018-2-7 22:45 | 显示全部楼层
能够分析1024吗?
plsbackup 发表于 2018-2-7 22:47 | 显示全部楼层
正常都是做到128的FFT。
tkyl01 发表于 2018-2-10 13:55 | 显示全部楼层
本帖最后由 tkyl01 于 2018-2-10 13:56 编辑

点数太少对低频效果不好,变换出频率误差较大,点数太多滞后比较明显,所以对于低频应用比如几HZ还是"寻峰+IIR"既简单又快速,比如计步应用,如果频率较高FFT还是很准的,STM32F4 的DSP库加上硬浮点快速计算等,1024点的FFT不到1ms 。
skawu 发表于 2018-2-10 14:46 | 显示全部楼层
学习了!!!
您需要登录后才可以回帖 登录 | 注册

本版积分规则

183

主题

2331

帖子

3

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