实验2matlab中基2DITFFT的实现1

《实验2matlab中基2DITFFT的实现1》由会员分享，可在线阅读，更多相关《实验2matlab中基2DITFFT的实现1（17页珍藏版）》请在装配图网上搜索。

1、电 子 科 技 大 学 实 验 报 告 学生姓名: 学 号:2010013080 指导教师 一、实验室名称:数字信号处理实验室 二、实验项目名称:FFT的实现 三、实验原理: 一(FFT算法思想: 1(DFT的定义: 对于有限长离散数字信号xn，0 , n , N-1，其离散谱xk可以由离散付氏变换(DFT)求得。DFT的定义为: 2,N,1,jnkNXkxne,，k=0,1,N-1 n,02,jN通常令，称为旋转因子。 eW,N2(直接计算DFT的问题及FFT的基本思想: 由DFT的定义可以看出，在xn为复数序列的情况下，完全直接运算N点2DFT需要(N-1)次复数乘法和N(N-1)次加法。

2、因此，对于一些相当大的N值(如1024)来说，直接计算它的DFT所作的计算量是很大的。 FFT的基本思想在于，将原有的N点序列分成两个较短的序列，这些序列的DFT可以很简单的组合起来得到原序列的DFT。例如，若N为偶数，将原有的N22点序列分成两个(N/2)点序列，那么计算N点DFT将只需要约(N/2) ?2=N/22次复数乘法。即比直接计算少作一半乘法。因子(N/2)表示直接计算(N/2)点DFT所需要的乘法次数，而乘数2代表必须完成两个DFT。上述处理方法可以反复使用，即(N/2)点的DFT计算也可以化成两个(N/4)点的DFT(假定N/2为偶数)，从而又少作一半的乘法。这样一级一级的划分

3、下去一直到最后就划分成两点的FFT运算的情况。 3(基2按时间抽取(DIT)的FFT算法思想: LN,2设序列长度为，L为整数(如果序列长度不满足此条件，通过在后面补零让其满足)。 L将长度为的序列，先按n的奇偶分成两组: N,2xnnN(0,1,.,1),xrxr2,1，r=0,1,N/2-1 xrxr21，,2DFT化为: NNN,1/21/21nkrkrk2(21)，XkDFTxnxnWxrWxrW221,，,NNNnrr,000NN/21/21,22rkkrk,，xrWWxrW,12NNN rr,00NN/21/21,rkkrk,，xrWWxrW,1/22/2NNNrr,002rkrk

4、WW,上式中利用了旋转因子的可约性，即:。又令NN/2NN/21/21,rkrk，则上式可以写成: XkxrWXkxrW,11/222/2NNrr,00k(k=0,1,N/2-1) XkXkWXk,，12N可以看出，分别为从中取出的N/2点偶数点和奇数点序列XkXk,Xk12的N/2点DFT值，所以，一个N点序列的DFT可以用两个N/2点序列的DFT组合而成。但是，从上式可以看出，这样的组合仅表示出了前N/2点的DFT值，Xk还需要继续利用XkXk,表示的后半段本算法推导才完整。利用旋转因Xk12rkrkN(/2)，WW,子的周期性，有:，则后半段的DFT值表达式: NN/2/2N/21/21

5、,NNrk()，NNrk2XkXk，,，同样， ，,XkxrWxrWXk22,11/21/21NN22,00rr(k=0,1,N/2-1)，所以后半段(k=N/2,N-1)的DFT值可以用前半段k值NN()，kkk22表达式获得，中间还利用到WWWW,，得到后半段的值表达式XkNNk为:(k=0,1,N/2-1)。 XkXkWXk,12Nxnxn,XkXk,这样，通过计算两个N/2点序列的N/2点DFT，可以1212组合得到N点序列的DFT值Xk，其组合过程如下图所示: k XkXkWXk，112Nnkk -1 XkWXkWXk,212NN比如，一个N = 8点的FFT运算按照这种方法来计算F

6、FT可以用下面的流程图来表示: x(0)X(0)0Wx(1)X(1)0Wx(2)X(2)2W0W0x(3)X(3)W1x(4)X(4)W0W2x(5)X(5)W0W3x(6)X(6)W2W0Wx(7)X(7) 4(基2按频率抽取(DIF)的FFT算法思想: LN,2设序列长度为，L为整数(如果序列长度不满足此条件，通过在后面补零让其满足)。 在把按k的奇偶分组之前，把输入按n的顺序分成前后两半: XkNNN,1/211nknknkXkDFTxnxnWxnWxnW,，,NNNnnnN,00/2NNN/21/21,()nk，Nnk2,，xnWxnW,NN 2nn,00NN/21,kNnk2,，,0

7、,1,.,1xnxnWWkN ,NN2n,0NNkk22W,1W,(1)因为，则有，所以: NNN/21,Nknk(1),0,1,.,1,，,，, XkxnxnWkN ,N2n,0按k的奇偶来讨论，k为偶数时: N/21,N2rn2,0,1,.,1,，, XrxnxnWkN ,N2n,0N/21,N(21)rn，21,0,1,.,1，,，, XrxnxnWkNk为奇数时: ,N2n,02rkrkWW,前面已经推导过，所以上面的两个等式可以写为: NN/2N/21,Nrn2,0,1,.,/21,，, XrxnxnWrN ,N/22n,0N/21,Nnnr21,0,1,.,/21，,，, Xrxn

8、xnWWrN ,NN/22n,0Xr2Xr21，通过上面的推导，的偶数点值和奇数点值分别可以由XkXr2组合而成的N/2点的序列来求得，其中偶数点值为输入xn的前半段和Xr21，后半段之和序列的N/2点DFT值，奇数点值为输入xn的前半段和后nW半段之差再与相乘序列的N/2点DFT值。 NNNnxnxnxnW,， xnxnxn,，令，则有: 2N122NN/21/21,Nrnrn2,21,0,1,.,1,，, XrxnWXrxnWr ,1/22/2NN2nn,00这样，也可以用两个N/2点DFT来组合成一个N点DFT，组合过程如下图所示: Nxnxn， xn2NNnnxn，xnxnW,， -1

9、 WNN22二(在FFT计算中使用到的MATLAB命令: 函数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)的情况。 四、实验目的: 离散傅氏变换(DFT)的目的是把信号由时域变换到频域，从而可以在频域分析处理信息，得到的结果再由逆DFT变换到时域。FFT是DFT的一种快速算法。在数字信号处理系统中，FFT作为一个非常重要的工具经常使用，甚至成为DSP运算能力的一个考核

10、因素。 本实验通过直接计算DFT，利用FFT算法思想计算DFT，以及使用MATLAB函数中的FFT命令计算离散时间信号的频谱，以加深对离散信号的DFT变换及FFT算法的理解。 五、实验内容: 5,a) 计算实数序列的256点DFT。 ,xnnn()cos,025616b) 计算周期为1kHz的方波序列(占空比为50,，幅度取为，/-512，采样频率为25kHz，取256点长度) 256点DFT。 六、实验器材(设备、元器件): 安装MATLAB软件的PC机一台，DSP实验演示系统一套。 七、实验步骤: (1) 先利用DFT定义式，编程直接计算2个要求序列的DFT值。 (2) 利用MATLAB中

11、提供的FFT函数，计算2个要求序列的DFT值。 (3) (拓展要求)不改变序列的点数，仅改变DFT计算点数(如变为计算1024点DFT值)，观察画出来的频谱与前面频谱的差别，并解释这种差别。通过这一步骤的分析，理解频谱分辨力的概念，解释如何提高频谱分辨力。 (4) 利用FFT的基本思想(基2,DIT或基2,DIF)，自己编写FFT计算函数，并用该函数计算要求序列的DFT值。并对前面3个结果进行对比。 (5) (拓展要求)尝试对其他快速傅立叶变换算法(如Goertzel算法)进行MATLAB编程实现，并用它来计算要求的序列的DFT值。并与前面的结果进行对比。 (6) (拓展要求)在提供的DSP实

12、验板上演示要求的2种序列的FFT算法(基2,DIT)，用示波器观察实际计算出来的频谱结果，并与理论结果对比。 八、实验数据及结果分析: 程序: (1) 对要求的2种序列直接进行DFT计算的程序 %第一种序列的计算 N=0:255; X=cos(5*pi*N/16); for a=1:256 Y(a)=0; for b=1:256 Y(a)=Y(a)+X(b)*exp(-j*2*pi*(b-1)*(a-1)/256); end end subplot(2,1,1) stem(N,abs(Y) title(DFT?) subplot(2,1,2) Y2=fft(X); stem(N,Y2) tit

13、le(FFT?) %第二种序列的计算 N=0:1/(1000*25):255/(1000*25); X=512*square(2*pi*N*1000); for a=1:256 Y(a)=0; for b=1:256 Y(a)=Y(a)+X(b)*exp(-j*2*pi*(b-1)*(a-1)/256); end end Y %?,fft? f=0:255; Y1=fft(X); subplot(2,1,1) stem(f,Y) title(DFT?) subplot(2,1,2) stem(f,Y1) title(FFT?) (2) 对要求的2种序列进行基2,DIT和基2,DIF FFT算法

14、程序 %基-2DIT-FFT的算法 %?-2-DIT-FFT clear clc tic %?x?2M? N=input(N=); x=input(x(?,?n=0:N-1)=); M=nextpow2(length(x); N=2M; n=0:N-1; x=x,zeros(1,N-length(x); %?x? LH=N/2; j1=LH; N1=N-2; for i=1:N1 if(i=k) j1=j1-k; k=k/2; end j1=j1+k; end %? for L=1:M %?,? B=2(L-1); for i=0:B-1 %?,? p=i*2(M-L); for k=i:2L

15、:(2(M-L)-1)*(2L)+i %?,?.?L?2(M-L)? temp=x(k+1); x(k+1)=temp+x(k+1+B)*exp(-j*2*pi*p/N); x(k+1+B)=temp-x(k+1+B)*exp(-j*2*pi*p/N); end end end stem(n,x) title(?-2-DIT-FFT?) time=toc %基-2-DIF-FFT的算法 %?-2-DIT-FFT clear clc tic %?x?2M? N=input(N=); x=input(x(?,?n=0:N-1)=); M=nextpow2(length(x); N=2M; n=0:

16、N-1; x=x,zeros(1,N-length(x); %?x? LH=N/2; j1=LH; N1=N-2; for i=1:N1 if(i=k) j1=j1-k; k=k/2; end j1=j1+k; end %? for L=1:M %?,? B=2(L-1); for i=0:B-1 %?,? p=i*2(M-L); for k=i:2L:(2(M-L)-1)*(2L)+i %?,?.?L?2(M-L)? temp=x(k+1); x(k+1)=temp+x(k+1+B)*exp(-j*2*pi*p/N); x(k+1+B)=temp-x(k+1+B)*exp(-j*2*pi*p

17、/N); end end end stem(n,x) title(?-2-DIT-FFT?) time=toc (3) 对要求的2种序列用MATLAB中提供的FFT函数进行计算的程序 N=0:255; w=0:4095; X=512*square(2*pi*1000*N/25000); y1=fft(X); y2=fft(X,4096); plot(N*2/256,abs(y1),o, w*2/4096,abs(y2) 结果:(1)对2种要求的序列直接进行DFT计算的频域波形 (2)对2种要求的序列进行基2,DIT和基2,DIF FFT算法频域波形 (3)对2种要求的序列用MATLAB中提供的

18、FFT函数计算的频域波形。 DFT的结果150100500050100150200250300FFT的结果150100500-50050100150200250300基-2-DIT-FFT的结果140120100806040200-200501001502002503004DFT的结果x 106420-20501001502002503004FFT的结果x 106420-20501001502002503004x 106543210-1-2050100150200250300(4)(拓展要求)分析利用上面的方法画出的信号频谱与理论计算出来的频谱之间的差异，并解释这种差异。 答:对于第一个函数，

19、理论上的频谱应该有很多的褶皱的。而用DFT画出，这就是因为采样没有采到其它的来的频谱就只有两点才有值，其他的值都为0点。 (5)(拓展要求)保持序列点数不变，改变DFT计算点数(变为1024点)，观察频谱的变化，并分析这种变化，由此讨论如何提高频谱分辨力的问题。 答:如下图所示，我们将第一个函数的采样点提高了。这个时候，我们可以看到其实提高点数，就是增加在频域的采样点数，在频域0到2*pi范围内的采样点数。我们增加FFT变换的点数，实际上就是增加了采样率。 14012010080604020000.511.524x 10987654321000.511.52九、实验结论: 答:DFT就是将FT采样后的的序列点数。我们可以通过增加采样的点数来提高采样的精度。同时，通过本次试验，我自己写了基2-DIT-FFT和基2-DIF-FFT的代码。对于基2的FFT有了更加深刻的理解。通过对已知函数的DFT不同点数的计算，我们 十、总结及心得体会: 这次试验，我自己写了基2-DIT-FFT和基2-DIF-FFT的代码。对于基2的FFT有了更加深刻的理解。通过对已知函数的DFT不同点数的计算，我对频率分辨与DFT的关系了有了更加深刻的认识 十一、对本实验过程及方法、手段的改进建议: 报告评分: 指导教师签字: