试验matlab中基DIT-FFT的实现

1、1/ 13电子科技大、实验室名称：数字信号处理实验室二、实验项目名称：FFT的实现三、实验原理:j通常令eN2.直接计算DFT的问题及FFT的基本思想：由DFT的定义可以看出，在xn为复数序列的情况下，完全直接运算N点DFT需要（N-1）2次复数乘法和N（ N-1）次加法。因此，对于一些相当大的N值 （如1024）来说，直接计算它的DFT所作的计算量是很大的。FFT的基本思想在于，将原有的N点序列分成两个较短的序列，这些序列的DFT可以很简单的组合起来得到原序列的DFT例如，若N为偶数，将原有的N点序列分成两个（N/2）点序列，那么计算N点DFT将只需要约（N/2）2- 2=N2/2次复数乘法

2、。即比直接计算少作一半乘法。因子（N/2）2表示直接计算（N/2） 点DFT所需要的乘法次数，而乘数2代表必须完成两个DFT上述处理方法可以 反复使用，即（N/2）点的DFT计算也可以化成两个（N/4）点的DFT （假定N/2为偶数），从而又少作一半的乘法。这样一级一级的划分下去一直到最后就划分 成两点的FFT运算的情况。3.基2按时间抽取（DIT）的FFT算法思想：设序列长度为N2L，L为整数（如果序列长度不满足此条件，通过在后面 补零让其满足）。学生姓名:学 号：2010013080指导教师.FFT算法思想：1.DFT的定义：对于有限长离散数字信号x n,On N-1,DFT求得。DFT的

3、定义为：1j2-nk0eN,k=0,1,其离散谱xk可以由离散付氏变换（NXknN-1WN，称为旋转因子。2/ 13将长度为N才的序列xn(n 0,1,., N 1)，先按n的奇偶分成两组:X2r X1r，r=0,1，,N/2-1x2r 1 X2rDFT化为:N/2 1X2rWN寫，则上式可以写成:r 0Xk Xik WrkX2k(k=0,1,-,N/2-1)可以看出，X1k, X2k分别为从Xk中取出的N/2点偶数点和奇数点序列的N/2点DFT值，所以，一个N点序列的DFT可以用两个N/2点序列的DFT组合 而成。但是，从上式可以看出，这样的组合仅表示出了Xk前N/2点的DFT值,还需要继续

4、利用X1k, X2k表示Xk的后半段本算法推导才完整。利用旋转因NN/2 1r(_Nk)N/2 1NX1- kX1rWN/：X1rWNk/2Xjk，同样，k Xzik2r 0r 02(k=0,1,N/2-1)，所以后半段(k=N/2,N-1)的DFT值可以用前半段k值 表达式获得，中间还利用到wNk)wNwkwk，得到后半段的Xk值表达式为：Xk X1k wNkX2k(k=0,1,N/2-1)。这样，通过计算两个N/2点序列x,n, X2n的N/2点DFlX/k, X2k，可以组合得到N点序列的DFT值Xk，其组合过程如下图所示:N 1IXk DFTx nxnw,kn 0N/2 1wkX2rW

5、N2rkr 0N /2 1wX2rWN怯N/2 1x2r WN2rkr 0N/2 1x2r 1WN(2r1)kr 0N/2 1Xj rw(rkr 0N/2 1X1门 wNk/2r 0上式中利用了旋转因子的可约性，即：wN2rkwNk/2。又令N/2 1XdkX1rWNr；,X2kr 0子的周期性，有：wNk/2WN駡N/2)，则后半段的DFTfi表达式:3/ 13比如，一个N = 8点的FFT运算按照这种方法来计算FFT可以用下面的流程图来表示:4基2按频率抽取（DIF）的FFT算法思想:设序列长度为N 2L,L为整数（如果序列长度不满足此条件，通过在后面 补零让其满足） 。在把Xk按k的奇偶

6、分组之前，把输入按n的顺序分成前后两半:N 1N/2 1N 1Xk DFTx nxnWNnkxnWr?kxnWr?kn 0n 0n N/2N/2 1kN/2 1N(nN)kxnWNnkxn -WN2n 0n 02N/2 1N-kkxn xn -WN2gWNk,k 0,1,.,N 1kWNX2kx(0)x(7)X(0)WNX2kx(1) X(1)x(2)X(2)x(3)X(3)XXx(5) X(5)x(6)X(6)4/ 13n 02NN因为WN21，则有WN( 1)k，所以：5/ 13前面已经推导过 Wjk，所以上面的两个等式可以写为:N/2 1NX2rxn xnNWNn/2,r0,1,.,N/

7、2 1n 02N/2 1zX2r 1n 0 xn xny曲WNnr/2，r0,1，.，N/2 1通过上面的推导，Xk的偶数点值 X2r和奇数点值 X2r 1分别可以由 组合而成的N/2点的序列来求得，其中偶数点值 X2r为输入xn的前半段和 后半段之和序列的N/2点DFT值，奇数点值 X2r 1为输入xn的前半段和后 半段之差再与WNJ相乘序列的N/2点DFT值。令张如0号，x2nxn号Wn，则有:Xk按k的奇偶来讨论,X2rN/2 1xn(n 0k为偶数时：N/2 1x nn 01)kxn却曲,k 0,1，.，N 1xny gwN2rn, k0,1,.,N 1k为奇数时：X2rN/2 11n

8、 0 xn却曲k 0,1,.,N 1N/2 1X2rx1n gC2,X2r 1N/2 1X2n(Wn2, rc . N ,0,1,., 12示：xn Nn 0n 0这样，也可以用两个N/2点DFT来组合成一个N点DFT6/ 13二.在FFT计算中使用到的MATLA命令：函数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的目的是把信号

9、由时域变换到频域，从而可以在频域分析处理信息，得到的结果再由逆DFT变换到时域。FFT是DFT的一种快速算法。在数字信号处理系统中，FFT作为一个非常重要的工具经常使用，甚至成为DSP运算能力的一个考核因素。本实验通过直接计算DFT利用FFT算法思想计算DFT以及使用MATLAB函数中的FFT命令计算离散时间信号的频谱，以加深对离散信号的DFT变换及FFT算法的理解。五、实验内容:cosP0 n 256的256点DF聽b）计算周期为1kHz的方波序列（占空比为50%，幅度取为+/-512，采样频率为25kHz,取256点长度）256点DFT六、实验器材（设备、元器件）：安装MATLAB软件的P

10、C机一台，DSP实验演示系统一套。七、实验步骤:先利用DFT定义式，编程直接计算2个要求序列的DFT值。 利用MATLAB提供的FFT函数，计算2个要求序列的DFT值。（拓展要求）不改变序列的点数，仅改变DFT计算点数（如变为计 算1024点DFT值），观察画出来的频谱与前面频谱的差别，并解 释这种差别。通过这一步骤的分析，理解频谱分辨力的概念，解 释如何提高频谱分辨力。利用FFT的基本思想（基2-DIT或基2-DIF），自己编写FFT计 算函数，并用该函数计算要求序列的DFT值。并对前面3个结果 进行对比。（拓展要求）尝试对其他快速傅立叶变换算法（如Goertzel算法） 进行MATLAB程

11、实现，并用它来计算要求的序列的DFT值。并与 前面的结果进行对比。（拓展要求）在提供的DSP实验板上演示要求的2种序列的FFT算法（基2-DIT），用示波器观察实际计算出来的频谱结果，并 与理论结果对比。八、实验数据及结果分析：程序：a）计算实数序列x（n）(1)(2(5)(6)7/ 13(1)对要求的2种序列直接进行DFT计算的程序%第一种序列的计算N=0:255; X=cos(5*pi*N/16);for a=1:256Y(a)=0;for b=1:256Y(a)=Y(a)+X(b)*exp(-j*2*pi*(b-1)*(a-1)/256);end end subplot(2,1,1)st

12、em(N,abs(Y)title(DFT卩？a 1?)subplot(2,1,2)Y2=fft(X);stem(N,Y2)title(FFT卩?a 1?)%第二种序列的计算N=0:1/(1000*25):255/(1000*25);X=512*square(2*pi*N*1000);for a=1:256Y(a)=0;for b=1:256Y(a)=Y(a)+X(b)*exp(-j*2*pi*(b-1)*(a-1)/256);end end Y%?-ffto ie ?卩?a 1?卩?士e f=0:255;Y1=fft(X); subplot(2,1,1) stem(f,Y)title(DFT卩

13、?a 1?)subplot(2,1,2)stem(f,Y1)title(FFT卩?a 1?)2) 对要求的2种序列进行基2DIT和基2DIF FFT算法程序8/ 13%基-2DIT-FFT的算法%?u -2-DIT-FFT clear clc tic %?e ? e ?D6a Dx21 a ? a 2M?N=input(? e ?e ?卩?e yN=：x=i np ut(? e ?e ?D6a Dx(?D 6 ?-? 6 ?a ?n=0:N-1)=);M=nextpow2(length(x);N=2M;n=0:N-1;x=x,zeros(1,N-length(x);%?e ?e ?卩?D6a

14、Dx?DD?D?D)LH=N/2;j1=LH;N1=N-2;for i=1:N1if (i=k)j1=j1-k;k=k/2;endj1=j1+k;endB=2A(L-1);%i ? 6 ?Dyx a 6 6X6卩?-?-p=i*2A(M-L);fork=i：2AL：(2A(M-L)-1)*(2L)+i%i ?6 ? i ?6 ?DyX a 6 6X6卩?-? . ?u L?D i ?6?DyX a6 6 X 6 3?2A(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

15、*pi*p/N);endend end stem(n,x)for L=1:M%?e y卩?-?-for i=0:B-19/ 13title(? u -2-DIT-FFT卩? a 1?)time=toc%基-2-DIF-FFT的算法%?u -2-DIT-FFT clear clc tic %?e ? e ?D6a Dx21 a ? a 2M?N=input( ? e ?e ?! ?e yN=);x=i np ut(? e ?e ?D6a Dx(?D 6 ?-? r ?a ?n=0:N-1)=);M=nextpow2(length(x);N=2M;n=0:N-1;x=x,zeros(1,N-len

16、gth(x);%?e ?e ?! ?D6a Dx?DD?D?D6LH=N/2;j1=LH;N1=N-2; for i=1:N1if (i=k)j1=j1-k;k=k/2;endj1=j1+k;endp=i*2A(M-L);fork=i：2AL：(2A(M-L)-1)*(2AL)+i%i ?6 ? i ?6 ?DyX a66X6! ?-? . ?u L?D i ?6?DyX a6 6 X 6 3?2A(M-L) ?temp=x(k+1);x(k+1)=temp+x(k+1+B)*exp(-j*2*pi*p/N);for L=1:M%?e y卩?-?-B=2A(L-1);for i=0:B-1%1

17、 ? 6 ?Dyx a 6 6X6(! ?-?-10/ 13x(k+1+B)=te mp-x(k+1+B)*ex p(-j*2* pi*p/N); endend(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);p lot(N*2/256,abs(y1),o , w*2/4096,abs(y2)结果：(1)对2种要求的序列直接进行DFT计算的频域波形(2) 对2种要求的序列进行基2-DIT和基2-DIF FFT算法频域波形(3)

18、 对2种要求的序列用MATLAB中提供的FFT函数计算的频域波形。endstem( n,x) title(?time=tocu -2-DIT-FFT卩？a 1?)1501BLIL一一一一-：+ ；FI I；.-.-77 T： T T r r ! ! -J*-J*-DFT的结果1005050100150200250300FFT的结果11/ 13140-200基-2-DIT-FFT的结果12010080604020050100150200250300DFT的结果4x 10 x 10-201 1 1 1L L1JJIraj一秤F二t单1丄E忙讣吧1 1rt1rFFT的结果42050100150200250430012/ 1313/ 134（4） （拓展要求）分析利用上面的方法画出的信号频谱与理论计算