加载中…
个人资料
微博
  • 博客等级：
  • 博客积分：
  • 博客访问：
  • 关注人气：
  • 获赠金笔：0支
  • 赠出金笔：0支
  • 荣誉徽章：
正文 字体大小：大 中 小

算法：FFT的物理意义及C代码

(2018-11-22 11:11:38)
分类： 信息技术研究

转：FFT结果的物理意义

    FFT是离散傅立叶变换的快速算法，可以将一个信号变换到频域。有些信号在时域上是很难看出什么特征的，但是如果变换到频域之后，就很容易看出特征了。这 就是很多信号分析采用FFT变换的原因。另外，FFT可以将一个信号的频谱提取出来，这在频谱分析方面也是经常用的。 
    虽然很多人都知道FFT是什么，可以用来做什么，怎么去做，但是却不知道FFT之后的结果是什意思、如何决定要使用多少点来做FFT。 
   一个模拟信号，经过ADC采样之后，就变成了数字信号。采样定理告诉我们，采样频率要大于信号频率的两倍，这些我就不在此罗嗦了。 
    采样得到的数字信号，就可以做FFT变换了。N个采样点，经过FFT之后，就可以得到N个点的FFT结果。为了方便进行FFT运算，通常N取2的整数次方。 其中，这N个采样点又称FFT size。
    假设采样频率为Fs，信号频率F，采样点数为N。那么FFT之后结果就是一个为N点的复数。每一个点就对应着一个频率点。这个点的模值，就是该频率值下的 幅度特性。具体跟原始信号的幅度有什么关系呢？假设原始信号的峰值为A，那么FFT的结果的每个点（除了第一个点直流分量之外）的模值就是A的N/2倍。 而第一个点就是直流分量，它的模值就是直流分量的N倍。而每个点的相位呢，就是在该频率下的信号的相位。第一个点表示直流分量（即0Hz），而最后一个点 N的再下一个点（实际上这个点是不存在的，这里是假设的第N+1个点，可以看做是将第一个点分做两半分，另一半移到最后）则表示采样频率Fs，这中间被 N-1个点平均分成N等份，每个点的频率依次增加。例如某点n所表示的频率为： 。由上面的公式可以看出，Fn所能分辨到频率为 Fs/N，如果采样频率Fs为1024Hz，采样点数为1024点，则可以分辨到1Hz。1024Hz的采样率采样1024点，刚好是1秒，也就是说，采样1秒时间的信号并做FFT，则结果可以分析到1Hz，如果采样2秒时间的信号并做FFT，则结果可以分析到0.5Hz。如果要提高频率分辨力（fft resolution)，则必须增加采样点数，也即采样时间。频率分辨率和采样时间是倒数关系。假设FFT之后某点n用复数a+bi表示，那么这个复数的模就是 ，相位就是 。根据以上的结果，就可以计算出n点（n≠1，且n<=N/2）对应的信号的表达式为： ，即 。对于n=1点的信号，是直流分量，幅度即为A1/N。由于FFT结果的对称性，通常我们只使用前半部分的结果，即小于采样频率一半的结果。 
    好了，说了半天，看着公式也晕，下面以一个实际的信号来做说明。 
    假设我们有一个信号，它含有2V的直流分量，频率为50Hz、相位为-30度、幅度为3V的交流信号，以及一个频率为75Hz、相位为90度、幅度为1.5V的交流信号。用数学表达式就是如下： 
S=2+3*cos(2*pi*50*t-pi*30/180)+1.5*cos(2*pi*75*t+pi*90/180) 
    式中cos参数为弧度，所以-30度和90度要分别换算成弧度。我们以256Hz的采样率对这个信号进行采样，总共采样256点。按照我们上面的分析，Fn=(n-1)*Fs/N，我们可以知道，每两个点之间的间距就是1Hz，第n个点的频率就是n-1。我们的信号有3个频率：0Hz、50Hz、75Hz，应该分别在第1个点、第51个点、第76个点上出现峰值，其它各点应该接近0。实际情况如何呢？
我们来看看FFT的结果的模值如图所示。 

图1 FFT结果
    从图中我们可以看到，在第1点、第51点、和第76点附近有比较大的值。我们分别将这三个点附近的数据拿上来细看：
1点： 512+0i
2点： -2.6195E-14 - 1.4162E-13i 
3点： -2.8586E-14 - 1.1898E-13i 
50点：-6.2076E-13 - 2.1713E-12i
51点：332.55 - 192i
52点：-1.6707E-12 - 1.5241E-12i 
75点：-2.2199E-13 -1.0076E-12i
76点：3.4315E-12 + 192i
77点：-3.0263E-14 +7.5609E-13i
   
    很明显，1点、51点、76点的值都比较大，它附近的点值都很小，可以认为是0，即在那些频率点上的信号幅度为0。接着，我们来计算各点的幅度值。分别计算这三个点的模值，
结果如下：
1点： 512
51点：384
76点：192
    按照公式，可以计算出直流分量为：512/N=512/256=2；50Hz信号的幅度为：384/(N/2)=384/(256/2)=3；75Hz信号的幅度为192/(N/2)=192/(256/2)=1.5。可见，从频谱分析出来的幅度是正确的。
然后再来计算相位信息。直流信号没有相位可言，不用管它。先计算50Hz信号的相位，atan2(-192, 332.55)=-0.5236,结果是弧度，换算为角度就是180*(-0.5236)/pi=-30.0001。再计算75Hz信号的相位，atan2(192, 3.4315E-12)=1.5708弧度，换算成角度180*1.5708/pi=90.0002。可见，相位也是对的。
根据FFT结果以及上面的分析计算，我们就可以写出信号的表达式了，它就是我们开始提供的信号。 
    总结：假设采样频率为Fs，采样点数为N，做FFT之后，某一点n（n从1开始）表示的频率为：Fn=(n-1)*Fs/N；该点的模值除以N/2就是对应该频率下的信号的幅度（对于直流信号是除以N）；该点的相位即是对应该频率下的信号的相位。相位的计算可用函数atan2(b,a)计算。atan2(b,a)是求坐标为(a,b)点的角度值，范围从-pi到pi。要精确到xHz，则需要采样长度为1/x秒的信号，并做FFT。要提高频率分辨率，就需要增加采样点数，这在一些实际的应用中是不现实的，需要在较短的时间内完成分析。解决这个问题的方法有频率细分法，比较简单的方法是采样比较短时间的信号，然后在后面补充一定数量的0，使其长度达到需要的点数，再做FFT，这在一定程度上能够提高频率分辨力。

c代码：

文件fft.c

 

#include

#include "fft.h"


float *SIN_TAB;//定义正弦表的存放空间

int FFT_N = 1024;//定义采样点大小 

 

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

{

 struct compx c;

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

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

 return(c);

}


 

void create_sin_tab(float *sin_t,int PointNum)

{

  int i;

  SIN_TAB=sin_t;

  FFT_N=PointNum;

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

    SIN_TAB[i]=sin(2*PI*i/FFT_N);

}

 

float sin_tab(float pi)

{

  int n=0;

  float a=0;

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


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

    a=SIN_TAB[n];

  else if(n>FFT_N/4&&n

    {

     n-=FFT_N/4;

     a=SIN_TAB[FFT_N/4-n];

    }

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

    {

     n-=FFT_N/2;

     a=-SIN_TAB[n];

   }

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

    {

     n=FFT_N-n;

     a=-SIN_TAB[n];

   }


  return a;

}

 

float cos_tab(float pi)

{

   float a,pi2;

   pi2=pi+PI/2;

   if(pi2>2*PI)

     pi2-=2*PI;

   a=sin_tab(pi2);

   return a;

}

 

void FFT(struct compx *xin)

{

  int f,m,i,k,l,j=0;

  register int nv2,nm1;

  struct compx u,w,t;


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

   nm1=FFT_N-1;

   for(i=0;i

   {

    if(i

     {

      t=xin[j];

      xin[j]=xin[i];

      xin[i]=t;

     }

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

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

     {

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

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

     }

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

  }


  {

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

    f=FFT_N;

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

           ;

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

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

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

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

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

    u.imag=0.0;

    //w.real=cos(PI/lei);                  //不适用查表法计算sin值和cos值

    // w.imag=-sin(PI/lei);

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

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

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

     {

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

       {

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

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

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

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

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

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

       }

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

     }

   }

  }

 

}

 

fft.h

#ifndef FFT_H

#define FFT_H

                                                //定义福利叶变换的点数

#define PI 3.1415926535897932384626433832795028841971               //定义圆周率值

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

extern void create_sin_tab(float *sin_t,int PointNum);

extern void FFT(struct compx *xin);

#endif // FFT_H

程序主要分为两个部分，第一部分主要采用雷德算法来实现倒位序，第二部分主要是利用已经生成好的旋转矩阵做蝶形变换。

举例：工程应用 
为了找到基波位置，我们对第一次平均值做fft变换。第一次平均值即反映了ADC采样的值得整体趋势，只是在幅值上有所降低。ADC的定时器触发的采样频率为20KHz，40个采样点一次平均后的采样频率为500Hz，故fft的采样频率为500Hz（也即fft变换后最大频率不会达到500Hz）。 
fft的采样点为1024点。故fft的分辨率约为0.49Hz。下面对具体程序和实际情况分析。 
fft handle.c

void FFT_handle(u16 Input_voltage)

{


    int i = 0;


    s[t].real = Input_voltage;

    s[t].imag = 0;

    t++;

    if(t >= NPT)

    {

        t = 0;

        value_ACsignalMAX = 0;

        FFT(s);

        for(i=0;i

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

        //    s[i].real=(u16)(sqrt(s[i].real*s[i].real+s[i].imag*s[i].imag)/(i==0?NPT:NPT/2));

            table[i]=(u16)(sqrt(s[i].real*s[i].real+s[i].imag*s[i].imag)/(i==0?NPT:NPT/2));


            if(i == 0)

            {

                value_DCsignal = table[i];

            }

            else if(0 < i && i < NPT / 2)

            {

                if(table[i] > value_ACsignalMAX)

                {

                    if(i < 15) // 滤除高频干扰

                    {

                        value_ACsignalMAX = table[i];

                        MAX_palce = i;

                    }

                    else{}

                }

            }

            else{}    


                                                                jww=value_DCsignal;

            jwww=value_ACsignalMAX;

 

        }

首先要求出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。 

喜欢

0

阅读┊ 收藏 ┊ 喜欢▼ ┊打印┊举报/Report
  

新浪BLOG意见反馈留言板　欢迎批评指正

新浪简介 | About Sina | 广告服务 | 联系我们 | 招聘信息 | 网站律师 | SINA English | 产品答疑

Copyright © 1996 - 2022 SINA Corporation, All Rights Reserved

新浪公司 版权所有