电 子 科 技 大 学
实 验 报 告
学生姓名: 学 号:2010013080 指导教师
一、实验室名称:数字信号处理实验室 二、实验项目名称:FFT的实现 三、实验原理:
一.FFT算法思想:
1.DFT的定义:
对于有限长离散数字信号{x[n]},0 ? n ? N-1,其离散谱{x[k]}可以由离散付氏变换(DFT)求得。DFT的定义为:
X[k]??x[n]en?0N?1?j2?nkN,k=0,1,…N-1
通常令e?j2?N?WN,称为旋转因子。
2.直接计算DFT的问题及FFT的基本思想:
由DFT的定义可以看出,在x[n]为复数序列的情况下,完全直接运算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次复数乘法。即比直接计算少作一半乘法。因子(N/2)2表示直接计算(N/2)点DFT所需要的乘法次数,而乘数2代表必须完成两个DFT。上述处理方法可以反复使用,即(N/2)点的DFT计算也可以化成两个(N/4)点的DFT(假定N/2为偶数),从而又少作一半的乘法。这样一级一级的划分下去一直到最后就划分成两点的FFT运算的情况。
3.基2按时间抽取(DIT)的FFT算法思想:
设序列长度为N?2L,L为整数(如果序列长度不满足此条件,通过在后面补零让其满足)。
将长度为N?2L的序列x[n](n?0,1,...,N?1),先按n的奇偶分成两组:
x[2r]?x1[r]x[2r?1]?x2[r]DFT化为:
N?1n?0,r=0,1,…,N/2-1
X[k]?DFT{x[n]}??x[n]WN/2?1nkNN/2?1??r?0x[2r]W2rkNN/2?1??r?0(2r?1)kx[2r?1]WN???r?0x1[r]Wx1[r]W2rkN?WkNN/2?1?r?0r?02rkx2[r]WN
N/2?1?r?0rkN/2?WkNN/2?1?rkx2[r]WN/2上式中利用了旋转因子的可约性,即:WNN/?22rkrk?WN/2。又令
X1[k]??x[r]W1r?01rkNN,/X2[k]?2?x[r]Wr?0?/21rk2N,则上式可以写成: /2kX[k]?X1[k]?WNX2[k](k=0,1,…,N/2-1)
可以看出,X1[k],X2[k]分别为从X[k]中取出的N/2点偶数点和奇数点序列的N/2点DFT值,所以,一个N点序列的DFT可以用两个N/2点序列的DFT组合而成。但是,从上式可以看出,这样的组合仅表示出了X[k]前N/2点的DFT值,还需要继续利用X1[k],X2[k]表示X[k]的后半段本算法推导才完整。利用旋转因子的周期性,有:WN/2?WN/2rkr(k?N/2),则后半段的DFT值表达式:
NN/2?1N/2?1r(?k)NNrk2X[?k]?X2[k] ,同样,X1[?k]??x1[r]WN/2??x1[r]WN?X[k]2/2122r?0r?0(k=0,1,…,N/2-1),所以后半段(k=N/2,…,N-1)的DFT值可以用前半段k值表达式获得,中间还利用到WN(N?k)2得到后半段的X[k]值表达式?WWk??Wk,
N2Nk为:X[k]?X1[k]?WN。 X2[k](k=0,1,…,N/2-1)
这样,通过计算两个N/2点序列x1[n],x2[n]的N/2点DFTX1[k],X2[k],可以组合得到N点序列的DFT值X[k],其组合过程如下图所示:
k X1[k] X1[k]?WNX2[k]
knk X2[k] WN -1 X1[k]?WNX2[k]
比如,一个N = 8点的FFT运算按照这种方法来计算FFT可以用下面的流程图来表示:
x(0)W0x(1)W0x(2)W0x(3)x(4)W0x(5)W0x(6)W0x(7)W2X(0)X(1)X(2)W2W0W1W2W3X(3)X(4)X(5)X(6)X(7)
4.基2按频率抽取(DIF)的FFT算法思想:
设序列长度为N?2L,L为整数(如果序列长度不满足此条件,通过在后面补零让其满足)。
在把X[k]按k的奇偶分组之前,把输入按n的顺序分成前后两半:
X[k]?DFT{x[n]}??x[n]Wn?0N/2?1N?1nkNN/2?1??n?0x[n]WnkN?n?N/2?N?1x[n]WNnk??因为WN2N?n?0n?0x[n]WNnk?N/2?1?n?0)kN(n?N2x[n?]WN2
N/2?1?kNN2[x[n]?x[n?]WN]?WNnk,k?0,1,...,N?12Nk2N??1,则有W?(?1)k,所以:
N/2?1X[k]???n?0Nnk[x[n]?(?1)x[n?]]?WN,k?0,1,...,N?1
2k按k的奇偶来讨论,k为偶数时:
N/2?1X[2r]??n?0[x[n]?x[n?N2rn]]?WN,k?0,1,...,N?1 2N/2?1k为奇数时:X[2r?1]?前面已经推导过WNN/2?1?n?0N(2r?1)n[x[n]?x[n?]]?WN,k?0,1,...,N?1
22rkrk?WN/2,所以上面的两个等式可以写为:
X[2r]??n?0[x[n]?x[n?Nrn]]?WN/2,r?0,1,...,N/2?1 2Nnnr]]?WN}WN/2,r?0,1,...,N/2?1 2N/2?1X[2r?1]??{[x[n]?x[n?n?0通过上面的推导,X[k]的偶数点值X[2r]和奇数点值X[2r?1]分别可以由组合而成的N/2点的序列来求得,其中偶数点值X[2r]为输入x[n]的前半段和后半段之和序列的N/2点DFT值,奇数点值X[2r?1]为输入x[n]的前半段和后半段之差再与WN相乘序列的N/2点DFT值。
令x1[n]?x[n]?x[n?N/2?1nNNWNn,则有: ],x2[n]?[x[n]?x[n?]]?22N/2?1X[2r]??n?0x1[n]?WrnN/2,X[2r?1]??n?0rnx2[n]?WN/2,r?0,1,...,N?1 2这样,也可以用两个N/2点DFT来组合成一个N点DFT,组合过程如下图所示:
Nx[n] x[n]?x[n?]
2
NNn x[n?] -1 WNn [x[n]?x[n?]]WN
22
二.在FFT计算中使用到的MATLAB命令:
函数fft(x)可以计算R点序列的R点DFT值;而fft(x,N)则计算R点序列的N点DFT,若R>N,则直接截取R点DFT的前N点,若R 四、实验目的: 离散傅氏变换(DFT)的目的是把信号由时域变换到频域,从而可以在频域分析处理信息,得到的结果再由逆DFT变换到时域。FFT是DFT的一种快速算法。在数字信号处理系统中,FFT作为一个非常重要的工具经常使用,甚至成为DSP运算能力的一个考核因素。 本实验通过直接计算DFT,利用FFT算法思想计算DFT,以及使用MATLAB函数中的FFT命令计算离散时间信号的频谱,以加深对离散信号的DFT变换及FFT算法的理解。 五、实验内容: a) 计算实数序列x(n)?cos5?n,0?n?256的256点DFT。 16b) 计算周期为1kHz的方波序列(占空比为50%,幅度取为+/-512,采样 频率为25kHz,取256点长度) 256点DFT。 六、实验器材(设备、元器件): 安装MATLAB软件的PC机一台,DSP实验演示系统一套。 七、实验步骤: (1) 先利用DFT定义式,编程直接计算2个要求序列的DFT值。 (2) 利用MATLAB中提供的FFT函数,计算2个要求序列的DFT值。 (3) (拓展要求)不改变序列的点数,仅改变DFT计算点数(如变为计 算1024点DFT值),观察画出来的频谱与前面频谱的差别,并解释这种差别。通过这一步骤的分析,理解频谱分辨力的概念,解释如何提高频谱分辨力。 (4) 利用FFT的基本思想(基2-DIT或基2-DIF),自己编写FFT计 算函数,并用该函数计算要求序列的DFT值。并对前面3个结果进行对比。 (5) (拓展要求)尝试对其他快速傅立叶变换算法(如Goertzel算法) 进行MATLAB编程实现,并用它来计算要求的序列的DFT值。并与前面的结果进行对比。 (6) (拓展要求)在提供的DSP实验板上演示要求的2种序列的FFT 算法(基2-DIT),用示波器观察实际计算出来的频谱结果,并与理论结果对比。