实验2_FFT的实现

电 子 科 技 大 学实 验 报 告学生姓名： 学 号： 指导教师：一、实验室名称：数字信号处理实验室 二、实验项目名称：FFT 的实现 三、实验原理：一．FFT 算法思想：1．DFT 的定义：对于有限长离散数字信号{x[n]}，0 ≤ n ≤ N-1，其离散谱{x[k]}可以由离散付氏变换（DFT ）求得。

DFT 的定义为：21[][]N jnk Nn X k x n eπ--==∑，k=0,1,…N-1通常令2jNN eW π-=，称为旋转因子。

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]=N 2/2次复数乘法。

即比直接计算少作一半乘法。

因子（N/2）2表示直接计算（N/2）点DFT 所需要的乘法次数，而乘数2代表必须完成两个DFT 。

上述处理方法可以反复使用，即（N/2）点的DFT 计算也可以化成两个（N/4）点的DFT （假定N/2为偶数），从而又少作一半的乘法。

这样一级一级的划分下去一直到最后就划分成两点的FFT 运算的情况。

3．基2按时间抽取（DIT ）的FFT 算法思想：设序列长度为2L N =，L 为整数（如果序列长度不满足此条件，通过在后面补零让其满足）。

将长度为2L N =的序列[](0,1,...,1)x n n N =-，先按n 的奇偶分成两组：12[2][][21][]x r x r x r x r =+=，r=0,1,…,N/2-1DFT 化为：1/21/212(21)0/21/21221200/21/211/22/2[]{[]}[][2][21][][][][]N N N nk rk r kNNNn r r N N rk k rkNNNr r N N rk k rkN NN r r X k DFT x n x n Wx r Wx r W x r W Wx r W x r WWx r W ---+===--==--=====++=+=+∑∑∑∑∑∑∑上式中利用了旋转因子的可约性，即：2/2rkrkNN W W =。

又令/21/2111/222/2[][],[][]N NrkrkNN r r X k x r WX k x r W--====∑∑，则上式可以写成： 12[][][]kN X k X k W X k =+(k=0,1,…,N/2-1)可以看出，12[],[]X k X k 分别为从[]X k 中取出的N/2点偶数点和奇数点序列的N/2点DFT 值，所以，一个N 点序列的DFT 可以用两个N/2点序列的DFT 组合而成。

但是，从上式可以看出，这样的组合仅表示出了[]X k 前N/2点的DFT 值，还需要继续利用12[],[]X k X k 表示[]X k 的后半段本算法推导才完整。

利用旋转因子的周期性，有：(/2)/2/2rkr k N N N W W +=，则后半段的DFT 值表达式：/21/21()211/21/2100[][][][]2N N N r k rk N N r r N X k x r W x r W X k --+==+===∑∑，同样，22[][]2N X k X k += （k=0,1,…,N/2-1），所以后半段（k=N/2,…,N-1）的DFT 值可以用前半段k 值表达式获得，中间还利用到()22NNk k k NNW W W W +==-，得到后半段的[]X k 值表达式为：12[][][]kN X k X k W X k =-（k=0,1,…,N/2-1）。

这样，通过计算两个N/2点序列12[],[]x n x n 的N/2点DFT 12[],[]X k X k ，可以组合得到N 点序列的DFT 值[]X k ，其组合过程如下图所示：1[]X k 12[][]kN X k W X k +2[]X k nkN W -1 12[][]k N X k W X k -比如，一个N = 8点的FFT 运算按照这种方法来计算FFT 可以用下面的流程图来表示：x(0)x(1)x(2)x(3)x(4)x(5)x(6)x(7)X(7)X(6)X(5)X(4)X(3)X(2)X(1)X(0)4．基2按频率抽取（DIF ）的FFT 算法思想：设序列长度为2L N =，L 为整数（如果序列长度不满足此条件，通过在后面补零让其满足）。

在把[]X k 按k 的奇偶分组之前，把输入按n 的顺序分成前后两半：1/211/2/21/21()200/212[]{[]}[][][][][]2[[][]],0,1,...,12N N N nknk nkNNN n n n N N N N n k nk N Nn n N N k nk N N n X k DFT x n x n Wx n Wx n W N x n W x n W N x n x n W W k N ---===--+==-====+=++=++=-∑∑∑∑∑∑因为21N NW=-，则有2(1)N k k NW=-，所以：/21[][[](1)[]],0,1,...,12N knkN n N X k x n x n W k N -===+-+=-∑按k 的奇偶来讨论，k 为偶数时：/2120[2][[][]],0,1,...,12N rn N n N X r x n x n W k N -==++=-∑k 为奇数时：/21(21)0[21][[][]],0,1,...,12N r nN n N X r x n x n W k N -+=+=-+=-∑前面已经推导过2/2rkrk NN W W =，所以上面的两个等式可以写为：/21/20[2][[][]],0,1,...,/212N rnN n N X r x n x n W r N -==++=-∑/21/20[21]{[[][]]},0,1,...,/212N n nr N N n N X r x n x n W W r N -=+=-+=-∑ 通过上面的推导，[]X k 的偶数点值[2]X r 和奇数点值[21]X r +分别可以由组合而成的N/2点的序列来求得，其中偶数点值[2]X r 为输入x[n]的前半段和后半段之和序列的N/2点DFT 值，奇数点值[21]X r +为输入x[n]的前半段和后半段之差再与nN W 相乘序列的N/2点DFT 值。

令1[][][]2N x n x n x n =++，2[][[][]]2nN N x n x n x n W =-+，则有：/21/211/22/20[2][],[21][],0,1,...,12N N rnrnN N n n NX r x n WX r x n W r --===+==-∑∑这样，也可以用两个N/2点DFT 来组合成一个N 点DFT ，组合过程如下图所示：[]x n [][]2Nx n x n ++[]2N x n + -1 n N W [[][]]2nN N x n x n W -+二．在FFT 计算中使用到的MATLAB 命令：函数fft(x)可以计算R 点序列的R 点DFT 值；而fft(x,N)则计算R 点序列的N 点DFT ，若R>N ，则直接截取R 点DFT 的前N 点，若R<N ，则x 先进行补零扩展为N 点序列再求N 点DFT 。

函数ifft(X)可以计算R 点的谱序列的R 点IDFT 值；而ifft(X,N)同fft(x,N)的情况。

四、实验目的：离散傅氏变换（DFT ）的目的是把信号由时域变换到频域，从而可以在频域分析处理信息，得到的结果再由逆DFT 变换到时域。

FFT 是DFT 的一种快速算法。

在数字信号处理系统中，FFT 作为一个非常重要的工具经常使用，甚至成为DSP 运算能力的一个考核因素。

本实验通过直接计算DFT ，利用FFT 算法思想计算DFT ，以及使用MATLAB 函数中的FFT 命令计算离散时间信号的频谱，以加深对离散信号的DFT 变换及FFT 算法的理解。

五、实验内容：a) 计算实数序列5()cos,025616x n n n π=≤≤的256点DFT 。

b) 计算周期为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 ），用示波器观察实际计算出来的频谱结果，并与理论结果对比。

八、实验数据及结果分析：注：本次实验在寝室电脑上完成，所用MATLAB版本为MATLAB R2010b 程序：（2）对要求的2种序列进行基2－DIT和基2－DIF FFT算法程序序列生成的代码同（1）。

这里给出FFT算法的程序：（3）对要求的2种序列用MATLAB中提供的FFT函数进行计算的程结果：（1）对2种要求的序列直接进行DFT计算的频域波形50100150200250300-1-0.500.51Time index nA m p l i t u d e Sinusoidal, time-domain sequenceFrequency index kM a g n i t u d eMagnitude of DFT samples(direct method)050100150200250300-1000-50005001000Time index nA m p l i t u d ePeriodic square wave, time-domain sequence4Frequency index kM a g n i t u d eMagnitude of DFT samples(direct method)（2） 对2种要求的序列进行基2－DIT 和基2－DIF FFT 算法频域波形（3） 对2种要求的序列用MATLAB 中提供的FFT 函数计算的频域波形。