手把手教你理解(FFT)

手把手教你理解(FFT)实验原理傅里叶变换是一种将信号从时域变换到频域的变换形式，是信号处理的重要分析工具。离散傅里叶变换（DFT）是傅里叶变换在离散系统中的表示形式。但是DFT的计算量非常大，FFT就是DFT的一种快速算法,FFT将DFT的N2步运算减少至(N/2)log2N步。离散信号x(n)的傅里叶变换可以表示为:10][)(NNnkNWnxkXNjNeW/2式中的WN称为蝶形因子，利用它的对称性和周期性可以减少运算量。一般而言，FFT算法分为时间抽取（DIT）和频率抽取（DIF）两大类。两者的区别是蝶形因子出现的位置不同，前者中蝶形因子出现在输入端，后者中出现在输出端。本实验以时间抽取方法为例。前者和后者如图所示：时间抽取FFT是将N点输入序列x(n)按照偶数项和奇数项分解为偶序列和奇序列。偶序列为：x(0),x(2),x(4),…,x(N-2)；奇序列为：x(1),x(3),x(5),…,x(N-1)。这样x(n)的N点DFT可写成：12/0)12(12/02122)(NnknNNnnkNWnxWnxkX考虑到WN的性质，即2/)2//(22/)2(2][NNjNjNWeeW因此有：12/02/12/02/122)(NnnkNkNNnnkNWnxWWnxkX或者写成：kZWkYkXkN)(由于Y(k)与Z(k)的周期为N/2，并且利用WN的对称性和周期性，即：kNNkNWW2/可得：kZWkYNkXkN)2/(对Y(k)与Z(k)继续以同样的方式分解下去，就可以使一个N点的DFT最终用一组2点的DFT来计算。在基数为2的FFT中，总共有log2(N)级运算，每级中有N/2个2点FFT蝶形运算。单个蝶形运算示意图如下：以N＝8为例，时间抽取FFT的信号流图如下：W0W0W2W0W2W0W1W2W3x(0)x(4)x(2)x(6)x(1)x(5)x(3)x(7)X(7)X(6)X(5)X(4)X(3)X(2)X(1)X(0)W0W0W0从上图可以看出，输出序列是按自然顺序排列的，而输入序列的顺序则是“比特反转”方式排列的。也就是说，将序号用二进制表示，然后将二进制数以相反方向排列，再以这个数作为序号。如011变成110，那么第3个输入值和第六个输入值就要交换位置了。本实验中采用了一种比较常用有效的方法完成这一步工作＿＿雷德算法。编程思路流程图：/*ProgramforFFT*/#includemath.h#includestdio.h#defineN64constfloatpi=3.1415926;/*FFT点数*/floatx_re[64],x_im[64];/*输入信号序列*/floaty_re[64],y_im[64];/*输出频谱序列*/floatyout[64];floatw_re,w_im;/*蝶形因子*/intm;/*蝶形运算的级数,即Log2(N)*/floatt_re,t_im,v_re,v_im;/*临时变量*/intj,i,k,f,n;inta,b,c;voidpaixu(floaty_re[N],floaty_im[N]){/*用雷德算法对输入信号序列进行倒序重排*/j=0;for(i=0;iN;i++){if(ij){t_re=y_re[j];t_im=y_im[j];y_re[j]=y_re[i];y_im[j]=y_im[i];y_re[i]=t_re;y_im[i]=t_im;}k=N/2;while((k=j)&(k0)){j=j-k;k=k/2;}j=j+k;}}voidfft(floaty_re[N],floaty_im[N]){/*计算蝶形运算的级数log2(N)*/f=N;for(m=1;(f=f/2)!=1;m++);/***FFT***/for(n=1;n=m;n++){a=1;/*a=2的n次方*/for(i=0;in;i++)a=a*2;b=a/2;v_re=1.0;/*蝶形因子*/v_im=0.0;w_re=cos(pi/b);w_im=-sin(pi/b);for(j=0;jb;j++)/*蝶形运算*/{for(i=j;iN;i=i+a){c=i+b;t_re=y_re[c]*v_re-y_im[c]*v_im;t_im=y_re[c]*v_im+y_im[c]*v_re;y_re[c]=y_re[i]-t_re;y_im[c]=y_im[i]-t_im;y_re[i]=y_re[i]+t_re;y_im[i]=y_im[i]+t_im;}t_re=v_re*w_re-v_im*w_im;t_im=v_re*w_im+v_im*w_re;v_re=t_re;v_im=t_im;}}}voidmain(){/*初始化数据空间*/for(i=0;iN;i++){x_re[i]=0;x_im[i]=0;}//输入的模拟信号for(i=0;i=N;i++){//x_re[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_im[i]=0;}/*复制到输出数组*/for(i=0;iN;i++){y_re[i]=x_re[i];y_im[i]=x_im[i];}paixu(y_re,y_im);//对输入数据进行倒序fft(y_re,y_im);//对输入数据的点数进行FFT运算for(i=0;iN/2;i++){yout[i]=sqrt(y_re[i]*y_re[i]+y_im[i]*y_im[i])*2/N;//计算幅值}while(1);}实验步骤1.以64点FFT的信号流图为例，理解FFT算法的过程；2.在CCS3.3环境中打开本实验的工程（Example_fft.pjt），编译并重建.out输出文件，然后通过仿真器把执行代码下载到DSP芯片中；3.运行程序；4.选择view-graph-time/frequency…。设置对话框中的参数:其中“StartAddress”设为“x_re”，“Acquisitionbuffersize”和“DisplayDatasize”都设为“64”，并且把“DSPDataType”设为“32-bitfloatingpoint”（如图），同样方法观察经DFT变换后的输出序列“y_re”的波形，“StartAddress”改为“y_re”，其余参数不变（如图）；采样定理告诉我们，采样频率要大于信号频率的两倍，假设采样频率为Fs，信号频率F，采样点数为N。模糊知识点的解释：1.采样出的点经过FFT运算后，各自点对应的频率是多少？答：也可以看做是将第一个点分做两半分，另一半移到最后）则表示采样频率Fs，这中间被N-1个点平均分成N等份，每个点的频率依次增加。例如某点n所表示的频率为：Fn=(n-1)*Fs/N。举例：Fn所能分辨到频率为为Fs/N，如果采样频率Fs为1024Hz，采样点数为1024点，则可以分辨到1Hz。024Hz的采样率采样1024点，刚好是1秒，也就是说，采样1秒时间的信号并做FFT，则结果可以分析到1Hz，如果采样2秒时间的信号并做FFT，则结果可以分析到0.5Hz。如果要提高频率分辨力，则必须增加采样点数，也即采样时间。频率分辨率和采样时间是倒数关系。2.采样出的点经过FFT运算后，各自点对应的频率是多少？答：每一个点就对应着一个频率点。这个点的模值，就是该频率值下的幅度特性。具体跟原始信号的幅度有什么关系呢？假设原始信号的峰值为A，那么FFT的结果的每个点（除了第一个点直流分量之外）的模值就是A的N/2倍。而第一个点就是直流分量，它的模值就是直流分量的N倍。3.信号的相位怎么通过FFT计算出的点推算出来？答：假设FFT之后某点n用复数a+bi表示，那么这个复数的模就是An=根号a*a+b*b，相位就是Pn=atan2(b,a)。（atan2）总结：根据以上的结果，就可以计算出n点（n≠1，且n=N/2）对应的信号的表达式为：An/(N/2)*cos(2*pi*Fn*t+Pn)，即2*An/N*cos(2*pi*Fn*t+Pn)。对于n=1点的信号，是直流分量，幅度即为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-1)*Fs/N，我们可以知道，每两个点之间的间距就是1Hz，第n个点的频率就是n-1。我们的信号有3个频率：0Hz、50Hz、75Hz，应该分别在第1个点、第51个点、第76个点上出现峰值，其它各点应该接近0。实际情况如何呢？我们来看看FFT的结果的模值如图所示。很明显，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.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。总结：假设采样频率为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语言的版本用VC6.0运行结果：（具体程序给我留言）模拟产生信号采样点值：test[i]=(float)10+1024*sin(2.0*i*PI/N-PI*30/180)+512*sin(4.0*i*PI/N+PI*90/180)+256*sin(6.0*i*PI/N+PI*40/180);运行结果：