实验十：FFT的实现与应用 信息工程学院 网络工程系 强文萍

一、实验目的 了解计算DFT算法存在的问题及改进途径。 掌握几种DFT算法（时间抽取算法DIT算法，频率抽取算法DIF算法）。 利用FFT计算线性卷积与圆周卷积；

二、实验仪器 微型计算机 MATLAB软件

三、实验原理 1．FFT算法 有限长序列通过离散傅里叶变换(DFT)将其频域离散化成有限长序列，但其计算量太大(与N的平方成正比）, 很难实时地处理问题, 因此引出了快速傅里叶变换(FFT)。 FFT并不是一种新的变换形式,它只是DFT的一种快速算法.并且根据对序列分解与选取方法的不同而产生了FFT的多种算法。 DFT的快速算法—FFT是数字信号处理的基本方法和基本技术，是必须牢牢掌握的。

三、实验原理 时间抽选FFT算法 前半部分 称为蝶形运算，WNk称为旋转因子 该算法遵循的准则： 对时间奇偶分；对频率前后分。 若取序列长度为2的幂次，则可以对序列逐级奇偶分解，最终分解为2点的DFT,再反向逐级利用蝶形算法最终得到原序列的DFT。 后半部分 前半部分

三、实验原理 群 群 群 蝶 这种算法的流图特点是： （1）基本运算单元都是蝶形 x(0) x(4) x(2) x(6) x(1) x(5)

三、实验原理 （2）同址（原位）计算 这是由蝶形运算带来的好处，每一级蝶形运算的结果 Xm+1(p)无须另外存储，只要再存入Xm(p)中即可，Xm+1(q) 亦然。这样将大大节省存储单元。 （3）倒位序 输入为“倒位序”（码位倒置）排列，输出按自然序排 列，因而对输入要进行“变址”计算（即码位倒置计算）。 “变址”实际上是一种“整序”的行为，目的是保证“同址”。

三、实验原理 FFT的应用 凡是利用付里叶变换来进行分析、综合、变换的地方，都可以利用FFT算法来减少其计算量。 FFT主要应用在 1、快速卷积 2、快速相关 3、频谱分析

三、实验原理 2．计算圆周卷积和线性卷积 1) 线性卷积 的长度为 两个有限长序列的线性卷积 线性卷积的长度

三、实验原理  2) 圆周卷积 两个有限长序列的N点圆周卷积 变换域计算方法 N 原理： x1(n) x2(n) IFFT FFT 的长度为 两个有限长序列的N点圆周卷积 变换域计算方法 N  原理： x1(n) x2(n) IFFT FFT N

三、实验原理 3) 用圆周卷积代替线性卷积 圆周卷积和线性卷积之间的关系： N

三、实验原理 4) 用FFT计算圆周卷积、线性卷积的步骤 圆周卷积 线性卷积

三、实验原理 分析线性卷积的运算量： 总计算量 次复数乘法

三、实验原理 线性卷积的运算量： 设x(n)的长度N1>N2 y(n)的长度为N=N1+N2-1 乘法运算量为N1·N

三、实验原理 比较两者乘法运算量 当两序列长度相当时 若长度较短时，如为8，16，32时,圆周卷积的时间大于直接线性卷积的结果； 当M=64时，两者的运算速度相当， 当M超过64以后，M越长圆周卷积的速度越快。 当输入序列x(n)较长时，可用重叠相加法和重叠保留法，进行分段卷积，保障圆周卷积的优势。

四、实验内容 1. 下面是实现基2的DIT-FFT算法的一段程序，分析代码，填写出来缺少的语句。 利用程序求序列x=[2,0,1,2]的DFT，和fft(x)的结果比较，验证程序的正确性

四、实验内容 function y=myditfft(x) %本程序对输入序列x实现DIT-FFT基2算法，点数取大于等于x长度的2的幂次。 m=nextpow2(length(x)); N=2^m; %求x的长度对应的2的最低幂次m if length(x)<N x=[x,zeros(1,N-length(x))]; %若x长度不是2的幂，补零到2的整数幂 end nxd=bin2dec(fliplr(dec2bin([1:N]-1,m)))+1; %求1:2^m数列的倒序 y= ; C=N/2; B=1; ; for D=0:m-1 % 逐级计算，C:级内群数 for c=0:C-1 % 逐群计算，B:群内蝶数 for b=0:B-1 % 逐蝶计算 POS = b + c * (B*2)+1; % 获得第一个数据的位置 TMP = y(POS + B) * WN^(b*C) ; y(POS) = xiang1 + TMP ; y(POS + B) = xiang1 - TMP ; C=C/2; B=B*2; % 级数总加一级，级内群数减少一倍

四、实验内容 2. 读入实验9中的音频文件w1.wav，分别截取n=64,1024,8096点，形成序列x1,x2,x3，分别利用 题1中的算法myditfft(x) 下面给出的依据定义求DFT的函数mydft(x) 求序列的DFT，并且比较所花销的时间 function y=mydft(x) % x为给定时间序列 % y为x的离散傅立叶变换 N=length(x); % 输入序列的长度 n=0:N-1; k=n; % 确定时域位置序列n和频域位置序列k WN=exp(-j*2*pi/N); % 计算DFT所需的旋转因子 nk=n'*k; WNnk=WN.^nk; % 构成旋转因子矩阵 Xk=x*WNnk; y=Xk; % 按DFT定义计算x的傅立叶变换

四、实验内容 比较算法时间 读入wav文件 clear [x, fs, nbits]=wavread('w1.wav',65536); tic,X=myditfft(x);toc %测试myditfft子程序所需运行时间 tic,X=mydft(x);toc %测试mydft子程序所需运行时间 读入wav文件 clear [x, fs, nbits]=wavread('w1.wav',65536); x1=x(1:64); %截短 将算法花销的时间填入下表 时间 Myditfft() Mydft() 64点 1024点 8096点

四、实验内容 3. 已知两序列 (1) 利用FFT计算 ，N=10,11,12,13 ，在figure(1)中画出4个图 3. 已知两序列 (1) 利用FFT计算 ，N=10,11,12,13 ，在figure(1)中画出4个图 (2) 直接计算两序列的线性卷积，在figure(2)中画出结果，比较看，和(1)中的哪些圆周卷积结果是一致的？ N

下面是对输入序列x1和x2实现N点圆周卷积的函数，请同学们补充完整。 function y=NCircleConv(x1,x2,N) %本程序对输入序列x1和x2实现N点圆周卷积，利用FFT算法 l1=length(x1); l2=length(x2); if N>l1 x1=[x1,zeros(1,N-l1)] end if N>l2 X1=fft(x1); X2=fft(x2); ; 调用这个函数完成本题。

直接计算两序列线性卷积的结果

五、实验报告要求 1.整理实验数据。 2.画出实验波形，并与各对应的频域的图形相比较。 3.完成题目后面的问题。

D表示当前级，c表示当前群，b表示当前蝶形单元。 C表示当前级内的群数，B表示当前群内的蝶形单元数。 FOR(D=0，C=N/2， B=1; D<R; D++){ // 逐级计算，C:级内群数 FOR(c=0; c<C; C++){ // 逐群计算，B:群内蝶数 FOR(b=0; b<B; b++){ // 逐蝶计算 POS = b + c * (B*2) //获得第一个数据的位置 TMP = DATA[POS + B] * W[b*C] ; DATA[POS] = DATA[POS] + TMP ; DATA[POS + B] = DATA[POS] - TMP ; } //END FOR{} } //END FOR{} C /= 2; //级数总加一级，级内群数减少一倍 B *= 2; //级数总加一级，群内蝶数增加一倍 } //END FOR{} N=2R是FFT的点数。 D表示当前级，c表示当前群，b表示当前蝶形单元。 C表示当前级内的群数，B表示当前群内的蝶形单元数。