电子科技大学-DSP-实验二-FFT的实现.docx

《电子科技大学-DSP-实验二-FFT的实现.docx》由会员分享，可在线阅读，更多相关《电子科技大学-DSP-实验二-FFT的实现.docx（7页珍藏版）》请在优知文库上搜索。

1、孑科技大学实验才一、实验室名称：数字信号处理实验室二、实验工程名称：FFT的实现三、实验原理：.FFT算法思想：1 .DFT的定义：对于有限长离散数字信号xn,0nNT,其离散谱xk可以由离散付氏变换（DFT）求得。DFT的定义为：v-_匹成X伙=Z4e,k=o,1,-N-I=0通常令J月=%,称为旋转因子。2 .直接计算DFT的问题及FFT的根本思想：由DFT的定义可以看出，在xn为复数序列的情况下，完全直接运算N点DFT需要（N-I）2次复数乘法和N（N-I）次加法。因此，对于一些相当大的N值（如1024）来说，直接计算它的DFT所作的计算量是很大的。FFT的根本思想在于，将原有的N点序列

2、分成两个较短的序列，这些序列的DFT可以很简单的组合起来得到原序列的DFT。例如，假设N为偶数，将原有的N点序列分成两个（N/2）点序列，那么计算N点DFT将只需要约（N2）2-2=M2次复数乘法。即比直接计算少作一半乘法。因子（N/2）2表示直接计算（N/2）点DFT所需要的乘法次数，而乘数2代表必须完成两个DFT。上述处理方法可以反复使用，即（N/2）点的DFT计算也可以化成两个（N/4）点的DFT（假定N/2为偶数），从而又少作一半的乘法。这样一级一级的划分下去一直到最后就划分成两点的FFT运算的情况。3 .基2按时间抽取（Dm的FFT算法思想：设序列长度为N=21.为整数（如果序列长度

3、不满足此条件，通过在后面补零让其满足）。将长度为N=2的序列伍=先按n的奇偶分成两组：42r=Xjr,r=0,l,-,N2-lx2r+l=x2rDFT化为：上式中利用了旋转因子的可约性，即：Wrk=Wf2。又令N/2-12-lXW=xirW2tX2k=X/团闱2,那么上式可以写成：r=0r=0Xk=Xyk+W,X2k(k=0,1,N2-1)可以看出，X伙,X2伙分别为从X灯中取出的N/2点偶数点和奇数点序列的N/2点DFT值，所以，一个N点序列的DFT可以用两个N/2点序列的DFT组合而成。但是，从上式可以看出，这样的组合仅表示出了X网前N/2点的DFT值，还需要继续利用XiX2k表示X眉的后

4、半段本算法推导才完整。利用旋转因子的周期性，有：W/2=,那么后半段的DFT值表达式：NN/2-Ir(Nk)N12-1NX%+即=W1%4=1r=X1Z:，同样，X2-+kX2k2r=0r=02(k=0,1,N/2-1),所以后半段(k=N2,N-D的DFT值可以用前半段k值表达式获得,中间还利用到w：A=W4wJ-8，得到后半段的X伙值表达式为：Xk=Xxk-WX2k(k=0,l,N2-l)o这样，通过计算两个N/2点序列卬叫马川的N/2点DFTXJ灯,X2k，可以组合得到N点序列的DFT值X灯，其组合过程如下列图所示：X2Wwk-1Xik-WX2k比方,一个N=8点的FFT运算按照这种方法

5、来计算FFT可以用下面的流程图来表示：4 .基2按频率抽取(DlF)的FFT算法思想：设序列长度为N=21.为整数(如果序列长度不满足此条件，通过在后面补零让其满足)。在把Xk按k的奇偶分组之前，把输入按n的顺序分成前后两半：NN_k因为WM=-1,那么有卬工=(T),所以：按k的奇偶来讨论，k为偶数时：N12-1Nk为奇数时：X2r+1=Z3川一和+力叱尸叫左=(U,.,NTM=o2前面已经推导过W；=Wf%,所以上面的两个等式可以写为：通过上面的推导，Xk的偶数点值X2r和奇数点值X2r+1分别可以由组合而成的N/2点的序列来求得，其中偶数点值X2为输入xn的前半段和后半段之和序列的N/2

6、点DFT值，奇数点值X2r+1为输入xn的前半段和后半段之差再与M相乘序列的N/2点DFT值。令XJm=X+#+y,x2n=xn-xn-W,那么有：乙乙这样，也可以用两个N/2点DFT来组合霹一个N卓DFT,组合过程如下列图所示：NN4+-1也3川一加+万吗二.在FFT计算中使用到的MAT1.AB命令：函数fft（x）可以计算R点序列的R点DFT值；而fft（x,N）那么计算R点序列的N点DFT,假设RN,那么直接截取R点DFT的前N点，假设RN,那么X先进行补零扩展为N点序列再求N点DFT。函数ifft（X）可以计算R点的谱序列的R点IDFT值；而ifft（X,N）同fft（x,N）的情况。

7、四、实验目的：离散傅氏变换（DFT）的目的是把信号由时域变换到频域，从而可以在频域分析处理信息，得到的结果再由逆DFT变换到时域。FFT是DFT的一种快速算法。在数字信号处理系统中，FFT作为一个非常重要的工具经常使用，甚至成为DSP运算能力的一个考核因素。本实验通过直接计算DFT,利用FFT算法思想计算DFT,以及使用MAT1.AB函数中的FFT命令计算离散时间信号的频谱，以加深对离散信号的DFT变换及FFT算法的理解。五、实验内容：a）计算实数序歹Ux（）=COS葛几0m256的256点DFTob）计算周期为IkHz的方波序列（占空比为50%,幅度取为+/-512,采样频率为25kHz,取

8、256点长度）256点DFT。六、实验器材（设备、元器件）：安装MAT1.AB软件的PC机一台，DSP实验演示系统一套。七、实验步骤：（1）先利用DFT定义式，编程直接计算2个要求序列的DFT值。（2）利用MAT1.AB中提供的FFT函数，计算2个要求序列的DFT值。（3）（拓展要求）不改变序列的点数，仅改变DFT计算点数（如变为计算1024点DFT值），观察画出来的频谱与前面频谱的差异，并解释这种差异。通过这一步骤的分析，理解频谱分辨力的概念，解释如何提高频谱分辨力。（4）利用FFT的根本思想基2DIT或基2DIF）,自己编写FFT计算函数，并用该函数计算要求序列的DFT值。并对前面3个结果

9、进行比照。（5）（拓展要求）尝试对其他快速傅立叶变换算法（如GOCrtZeI算法）进行MAT1.AB编程实现，并用它来计算要求的序列的DFT值。并与前面的结果进行比照。（6）（拓展要求）在提供的DSP实验板上演示要求的2种序列的FFT算法（基2-DIT）,用示波器观察实际计算出来的频谱结果，并与理论结果比照。八、实验数据及结果分析：程序：(1) 对要求的2种序列直接进行DFT计算的程序functionlXk=DFT（xn,n）Xk=zeros(l,n);if(length(xn)n)xn_tmp=xn,zeros(1,n-length(xn);elsexn_tmp=xn;endfori=0:n

10、-lforj=0m-lXk(i+l)=Xk(i+l)+xn_tmp(j+1)*exp(-1i*2*pi*j*in);endend(2) 对要求的2种序列进行基2-DIT和基2-DIFFFT算法程序functionXk=DIT(xn,n)if(length(xn)n)xn_tmp=xn,zeros(1,n-length(xn);elsexn_tmp=xn;endx1=xn_tmp(1,2:2:n);x2=xn_tmp(1,1:2:n-1);WNK=exp(-1i*2*pi*(0:n/2-1)n);Xlk=DFT(xl,n2);X2k=DFT(x2,n2);Xk_l=Xlk+X2k.*WNK;Xk

11、_2=Xlk-X2k.*WNK;Xk=Xk_l,Xk_2J;functionlXk=D!F(xn,n)if(length(xn)n)xn_tmp=xn,zeros(1,n-length(xn)J;elsexn_tmp=xn;endxn1=xn_tmp(l:n/2);xn2=xn-tmp(n2+l:n);xkl=DFT(xnl,n2);xk2=DFT(xn2,n2);WNK=exp(-1i*2*pi*(0:n/2-1)n);Xk=zeros(l,n);Xk(2:2:n)=xkl+xk2;Xk(l:2:n-l)=(xkl-xk2).*WNK;(3) 对要求的2种序列用MAT1.AB中提供的FFT函

12、数进行计算的程序clear;clc;n=0:256;kl=0:255;k2=0J023;N1=256;N2=1024;%cos256点DFTxn_cos=cos(5*pi.*n/16);Xk_cos=DFT(xn_cos,256);Xk_cos_FFT=fft(xn_cos,256);Xk_cos_1024=DFT(xn_cos,1024);figure(1),set(gcf,name,cos序列变换)；subplot(4,1,1),stem(n,xn_cos),axis(0,256,-1.2,1.2),title(xn,),xlabel(,n);subplot(4,1,2),stem(k1,

13、abs(Xk_cos),axis(0,255,l.2*min(abs(Xk-cos),1.2*max(abs(Xk-cos)J),title(,256点DFT求出的IXkl)Xlabel(k);subplot(4,1,3),stem(k1,abs(Xk_cos_FFT),axis(0,255,1.2*min(abs(Xk_cos_FFT),1.2*max(abs(Xk.cos.FFT),title(,256点FFT求出的IXkl)XIabel(k);subplot(4,1,4),stem(k2,abs(Xk_cos_1024),axis(0,l023,1.2*min(abs(Xk_cos_l0

14、24),1.2*max(abs(Xk_cos_1024)J),title(11024点DFT求出的IXkI)Xlabel(k);%IKhz方蕨序列256点DFTFs=25OOO;n_squ=0:l/Fs:256/Fs;xn_squ=512*square(2*pi*1000*n_squ);Xk_squ=DFT(xn_squ,N1);Xk_squ_1024=DFT(xn_squ,N2);Xk_squ_FFT=fft(xn_squ,N1);figure,set(gcf,name；方波序列变换);subplot(4,1,1),stem(n,xn-squ),axis(0,255,-550,550),title(,1kHz方波序列,),XlabeI(,n);subplot(4,l,2),stem(kl,abs(Xk_squ),axis(0,255,1.2*min(abs(Xk_squ),1.2*max(abs(Xk_squ)J),title(,256点DFT求出的IXkI)XlabeI(k);subplot(4,1,3),stem(k1,abs(Xk_squ_FFT),axis(0,255,1.2*min(abs(Xk_squ_FFT),1.2*max(abs(Xk_squ_FFT),titleC256点FFT求出的IXkl),xla