手把手教你理解(FFT)

手把手教你理解(FFT)实验原理傅里叶变换是一种将信号从时域变换到频域的变换形式，是信号处理的重要分析工具。离散傅里叶变换(DFT)是傅里叶变换在离散系统中的表示形式。但是DFT的计算量非常大，FFT就是DFT的一种快速算法,FFT将DFT的N?步运算减少至(N/2)log?N步。离散信号x(n)的傅里叶变换可以表示为：N-1X(R)=。［观jV=O式中的Wn称为蝶形因子，利用它的对称性和周期性可以减少运算量。一般而言,FFT算法分为时间抽取(DIT)和频率抽取(DIF)两大类。两者的区别是蝶形因子出现的位置不同，前者中蝶形因子出现在输入端，后者中出现在输出端。本实验以时间抽取方法为例。前者和后者如图所示：时间抽取FFT是将N点输入序列x(n)按照偶数项和奇数项分解为偶序列和奇序列。偶序列为：x(0),x(2),x(4),...,x(N-2)；奇序列为：x(l),x(3),x(5),...,x(N-l)o这样x(n)的N点DFT可写成：N/2-1 jV/2-1X（k）=£.伽）叱严+£x（2〃+l）W严於TOC\o"1-5"\h\z”=0 "=0考虑到Wn的性质，即叱=［严”竹2=严/(N⑵二％因此有：N/2-1 N/2-1X伙)=£x(2〃)w為+叱£x(2〃+l)叱2H=0 /1=0或者写成：X伙)二瞅)+W刚由于Y(k)与Z(k)的周期为N/2,并且利用Wn的对称性和周期性，即:W严4=_wk可得：X(k+N/2)二除)—W；Z(k)对Y(k)与Z(k)继续以同样的方式分解下去，就可以使一个N点的DFT最终用一组2点的DFT來计算。在基数为2的FFT中，总共有log2(N)级运算，每级中有N/2个2点FFT蝶形运算。单个蝶形运算示意图如下:Wn-1

Wn-1以N=8为例，时间抽取FFT的信号流图如下:x(0)x(4)xx(0)x(4)x⑵x(6)x(l)x(5)x⑶x(7)X(0)X(l)X⑵X⑶X⑷X⑸X⑹X⑺从上图可以看出，输出序列是按自然顺■列的,而输入序列的顺序则是“比特反转”方式排列的。也就是说，将字号用二进制表示，然后将二进制数以相反方向排列，再以这个数作为序号。如011变成110,那么第3个输入值和第六个输入值就要交换位置了。本实验中釆用了一种比较常用有效的方法完成这一步工作雷徳算法。

/*PiogiamforFFT*/#include<math.h>#include<stdio.h>#deflneN64constfloatpi=3.1415926;/*FFT点数*//*输入信号序列/*输入信号序列*//*输出频谱序列*/floatv_ie[64],y_im[64];floatvout[64];/*蝶形因子/*蝶形因子*//*蝶形运算的级数，即Log2(N)*//*临时变量*/iiitm;floatt_re,t_im,v_re?v_im;iiita,b,c;voidpaixu(floaty_re[N].floaty_im[N])/*用雷徳算法对输入信号序列进行倒序重排*/尸0；foi(i=0;i<N;i++)t_re=y_ie|j];y-ieU]=y_re[i];y_imU]=y_mi[i];y_re[i]=t_ie;}k=N/2;while((k<=^j)&(k>0)){j=j-k；k=k/2;

尸j+k;voidfft(floaty_ie[N].floatv_im[N]){/*计算蝶形运算的级数log2(N)*/f=N;foi(m=l;(仁f/2)!=l;m十十);t_ie=v_ie*w_ie-v_mi*w_ini;t_ini=v_re*w_im+v_mi*w_ie;v_ie=t_ie;***fft***/***fft***/foi(n=l;n<=m;n++){a=l; /*a=2foi(n=l;n<=m;n++){a=l; /*a=2的n次方*/foi(i=0;i<n;i++)a=a*2;b=a/2;v_ie=1.0; /*蝶形因子*/v_un=0.0;w_re=cos(pi/b);w_im=-sin(pi/b);foi(j=0;j<b;j卄)/*蝶形运算*/c=i十b;t_ie=v_ie[c]*v_ie-y_nn[c]t_im=y_ie[c]*v_im+y_im[c]*v_ie;y_re[c]=y_ie[i]-t_re;v_re[i]=y_ie[i]+t_re;v_un=t_im;}}}voidmaui(){/*初始化数据空间*/for(i=0;i<N;i十十){x_re[i]=0;x_im[i]=0;}〃输入的模拟信号for(i=0;i<=N;i++){//x_ie[i]=(2*cos(-16*pi*i/N)+cos(-10*pi*i/N)+5*cos(-24*pi*i/N));x_re[i]=4*cos(-10*pi*i/N)+2*cos(-16*pi*i/N);x_un[i]=0;}/*复制到输出数组*/for(i=0;i<N;i十十){y_ie[i]=x_ie[i];y_mi[i]=x_mi[i];}paixu(v_re,y_im); 〃对输入数据进行倒序fft(y_re,y_im); 〃对输入数据的点数进行FFT运算for(i=0;i<N/2;i-H-){yout[i]=sqrt(y_ie[i]*y_re[i]+y_im[i]*v_mi[i])*2/N; 〃计算幅值while(l);}实验步骤以64点FFT的信号流图为例，理解FFT算法的过程；在CCS3.3环境中打开本实验的工程(Example_fflpjt),编译并重建.out输出文件，然后通过仿真器把执行代码下载到DSP芯片中；运行程序；选择view->grapli->time/fiequency...。设置对话框中的参数：其中“StartAddress设为"x_re”，"Acquisitionbuffersize"和“DisplayDatasize”都设为“64"，并且把“DSPDataType”设为“32-bit floatuigpoint”(如图)，SSGraphPropertyDialogDisplayTypeSingleTimeAGraphTitleGraphicalDisplayStartAddressx_rePageDataAcq,uisitionBufferSize64言*Irkereirient1DisplayD：at色Size64DSPDataType32"bitfloatingpointSajnpliiigRateCHz)1PlotDataPromLefttoRightLefJisplayYesAutoscaleOnDCValue0AxesDisplOnTimeDisplayUnit5StatusBarDisplayOnV|QK|QancJ|Help|同样方法观察经DFT变换后的输出序列“y」e”的波形，“StartAddress"改为“y」e”，其余参数不变（如图）;釆样定理告诉我们，釆样频率要大于信号频率的两借，假设采样频率为Fs,信号频率F,釆样点数为N。模糊知识点的解释：1•采样出的点经过FFT运算后，各自点对应的频率是多少？答：也可以看做是将第一个点分做两半分，另一半移到最后）则表示釆样频率Fs,这中间被N-1个点平均分成N等份，每个点的频率依次增加。例如某点n所表示的频率为：Fn=（n-l）*Fs/No举例：Fn所能分辨到频率为为Fs/N,如果釆样频率Fs为1024Hz,釆样点数为1024点，则可以分辨到1Hz。024Hz的釆样率釆样1024点，刚好是1秒，也就是说，釆样1秒时间的信号并做FFT,则结果可以分析到1Hz,如果釆样2秒时间的信号并做FFT,则结果可以分析到0.5Hzo如果要提高频率分辨力，则必须增加采样点数，也即采样时间。频率分辨率和采样时间是倒数关系。2•采样出的点经过FFT运算后，各自点对应的频率是多少？答：每一个点就对应着一个频率点。这个点的模值，就是该频率值下的幅度特性。具体跟原始信号的幅度有什么关系呢？假设原始信号的峰值为A,那么FFT的结果的每个点（除了第一个点直流分量之外）的模值就是A的N/2倍。而第一个点就是直流分量，它的模值就是直流分量的N倍。

3•信号的相位怎么通过FFT计算出的点推算出来？答：假设FFT之后某点n用复数a+bi表示，那么这个复数的模就是An二根号a*a+b*b,相位就是Pn=atan2(b,a)。(atan2)0幵七求相位为什么差90度menggxingg|分类；工程技术科学|浏览429次 2012-06-060分李到；S分李到；S举报举报▼ 2012-06-061L丄提问者采纳正弦和余弦之间的相位正好差90S，可能是这个原因；如果你的信号模型是正弦，你所用的方法求得的是余弦型的，那么就正好差90度了，没问题的,据我说知，FFT求得的相位应该是余弦型的atan2返回给定的X及Y坐标值的反正切值。反正切的角度值等于X轴正方向与通过原点和给定坐标点(¥坐标,邂标)的射线之间笊夹角。结果以弧度表示芥介于叩i到pi之间(不包括叩i)。总结：根据以上的结果，就可以计算出n点(nHl,且nUN/2)对应的信号的表达式为：An/(N/2)*cos(2*pi*Fn*t+Pn),即2*An/N*cos(2*pi*Fn*t+Pn)。对于nh点的信号，是直流分量，幅度即为A1/N。由于FFT结果的对称性，通常我们只使用前半部分的结果，即小于采样频率一半的结果。借用高人用matlab做的一个来理解：假设我们有一个信号，它含有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-l)*Fs/N,我们可以知道，每两个点之间的间距就是1Hz,第n个点的频率就是n-l。我们的信号有3个频率：0Hz、50Hz、75Hz,应该分别在第1个点、第51个点、第76个点上出现峰值，其它各点应该接近0。实际情况如何呢？很明显，1点、51点、76点的值都比较大，它附近的点值都很小，可以认为是0,即在那些频率点上的信号幅度为0。接着，我们來计算各点的幅度值。分别计算这三个点的模值，结果如下：1点：51251点：38476点：192按照公式，可以计算出直流分量为：512/N二512/256二2；50Hz信号的幅度为:384/(N/2)=384/(256/2)=3：75Hz信号的幅度为192/(N/2)=192/(256/2)=1.5O可见,从频谱分析出来的幅度是正确的。然后再來计算相位信息。直流信号没有相位可言，不用管它。先计算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。总结：假设釆样频率为Fs,采样点数为N,做FFT之后，某一点n(n从1开始)表示的频率为：Fn=(n-l)*Fs/N；该点的模值除以N/2就是对应该频率下的信号的幅度(对于直流信号是除以N)；该点的相位即是对应该频率下的信号的相位。相位的计算可用函数atan2(b,a)计算。atan2(b,a)是求坐标为(a,b)点的角度值，范围从-pi到pi。要精确到xHz,则需要采样长度为1/x秒的信号，并做FFT。要提高频率分辨率，就需要增加釆样点数，这在一些实际的应用中是不现实的，需要在较短的时间内完成分析。解决这个问题的方法有频率细分法，比较简单的方法是采样比较短时间的信号，然后在后面补充一定数量的0,使其长度达到需要的点数，再做FFT,这在一定程度上能够提高频率分辨力。我写的另外一种C语言的版本用VC6.0运行结果：(具体程序给我留言)模拟