武汉理工大学8点基于DIFFFT.docx

武汉理工大学8点基于DIFFFT
课程设计任务书

学生姓名：李嘉辛专业班级：电信1206
指导教师：黄朝兵工作单位： 信息工程学院

题 目：8点基于DIF的FFT的实现
初始条件：
具备Matlab编程能力；
熟悉基于DIF的FFT的实现原理；
提供编程所需要的计算机一台
要求完成的主要任务：（包括课程设计工作量及其技术要求，以及说明书撰写等具体要求）
1、独立编写一个8点的基于DIF的FFT实现程序，不能使用matlab自带的FFT实现函数
2、并调用该函数实现16点的FFT运算，用matlab自带函数对运行结果结果进行验证
3、完成符合学校要求的设计说明书
时间安排：
一周，其中3天程序设计，2天程序调试
指导教师签名： 年 月 日

系主任（或责任教师）签名： 年 月 日

15 / 1715 / 17
目录
摘要 1
1概述 1
1.1数字信号处理定义 1
1.2数字信号处理的特点及实现方法 1
2 理论分析 2
2.1 DFT的定义 2
2.2 直接计算DFT的问题及FFT思想 2
2.3基2按时间抽取（DIT）的FFT算法 2
2.4基2按频率抽取（DIF）的FFT算法 4
2.5 按频率抽取的FFT的特点 6
2.5.1 原位运算 6
2.5.2 蝶形运算两节点之间的“距离” 6
2.5.3 旋转因子的变化规律 6
3 程序设计 7
3.1变址 7
3.2 L级递推计算 7
4结果及分析 9
5心得体会 12
参考文献 13

摘要
快速傅里叶变换（FFT）是离散傅里叶变换（DFT）的快速算法，FFT算法通过利用旋转因子的性质，将一个大点数DFT化成几个小点数DFT，就可以大大减少运算量。DIF-FFT是利用频率抽选的FFT算法，在Matlab中可以通过三重循环语句实现。
关键词：FFT，蝶形运算，倒序排列

1概述
1.1数字信号处理定义
数字信号是用数字序列表示的信号，数字信号处理就是通过计算机或专用处理设备，用数值计算等数字方式对数字序列进行各种处理，将数字信号变换成符合要求的某种形式。数字信号处理主要包括数字滤波和数字频谱分析两大部分。例如，对数字信号进行滤波，限制其频带或滤除噪声和干扰，以提取和增强信号的有用分量；对信号进行频谱分析或功率分析，了解信号的频谱组成，以对信号进行识别。当然，凡是用数字方式对信号进行滤波，变换，增强，压缩，估计和识别等都是数字信号处理的研究范畴。
数字信号处理在理论上所涉及的范围及其广泛。数字领域中的微积分，概率统计，随机过程，高等代数，数值分析，复变函数和各种变换（如傅里叶变换，Z变换，离散傅里叶变换，小波变换等）都是它的基本工具，网络理论，信号及系统等则是它的理论基础。在科学发展上，数字信号处理又和最优控制，通信理论等紧密相连，目前已成为人工智能，模式识别，神经网络等新兴学科的重要理论基础，其实现技术又和计算机科学和微电子技术密不可分。因此，数字信号处理是把经典的理论基础体系作为自身的理论基础，同时又使自己成为一系列新兴学科的理论基础。
1.2数字信号处理的特点及实现方法
及模拟信号处理相比，数字信号处理具有高精度、高稳定性、灵活性好、易于大规模集成等显著的优点。
数字信号处理的主要研究对象是数字信号，且采用数值运算的方法达到处理的目的。数字信号处理的实现方法基本上可以分为软件实现方法、硬件实现方法和软硬件想结合的实现方法。数字信号处理的理论、算法和实现方法三者是密不可分的。

2理论分析
2.1DFT的定义
对于有限长离散数字信号{x[n]}，0 £n£N-1，其离散谱{X[k]}可以由离散付氏变换（DFT）求得。DFT的定义为：
，k=0,1,…N-1 （2.1）
通常令，称为旋转因子。
2.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运算的情况。
2.3基2按时间抽取（DIT）的FFT算法
设序列长度为，L为整数（如果序列长度不满足此条件，通过在后面补零使其满足）。
将长度为的序列，先按n的奇偶分成两组：
（r=0,1,…,N/2-1） （2.2）
DFT化为：
(2.3)
上式中利用了旋转因子的可约性，即：。又令，则上式可以写成：
(k=0,1,…,N/2-1) (2.4)
可以看出，分别为从中取出的N/2点偶数点和奇数点序列的N/2点DFT值，所以，一个N点序列的DFT可以用两个N/2点序列的DFT组合而成。但是，从上式可以看出，这样的组合仅表示出了前N/2点的DFT值，还需要继续利用表示的后半段本算法推导才完整。利用旋转因子的周期性，有：，则后半段的DFT值表达式：
(2.5)
（k=0,1,…,N/2-1）(2.6)
所以后半段（k=N/2,…,N-1）的DFT值可以用前半段k值表达式获得，中间还利用到，得到后半段的值表达式为：（k=0,1,…,N/2-1）。
这样，通过计算两个N/2点序列的N/2点DFT，可以组合得到N点序列的DFT值，其组合过程如下图所示：

-1
图2.1 FFT形成过程
2.4基2按频率抽取（DIF）的FFT算法
设序列长度为，L为整数（如果序列长度不满足此条件，通过在后面补零使其满足）。
在把按k的奇偶分组之前，把输入按n的顺序分成前后两半：
(2.7)
因为，则有，所以：
X[k]=n=0N/2-1[xn+-1kx[n+N2]]WNnk,k=0,1,…,N-1(2.8)
按k的奇偶来讨论，k为偶数时：
X[2r]=n=0N/2-1[xn+x[n+N2]]WN2rn,k=0,1,…,N-1(2.9)
k为奇数时：
X[2r+1]=n=0N/2-1[xn-x[n+N2]]WN(2r+1)n,k=0,1,…,N-1 (2.10)
前面已经推导过，所以上面的两个等式可以写为：
X[2r]=n=0N/2-1[xn+x[n+N2]]WN/2rn,k=0,1,…,N/2-1 (2.11)
X[2r+1]=WN/2rnn=0N/2-1[xn-x[n+N2]]WNrn,k=0,1,…,N/2-1(2.12)
通过上面的推导，X[k]的偶数点值X[2r]和奇数点值X[2r+1]分别可以由组合而成的N/2点的序列来求得，其中偶数点值X[2r]为输入x[n]的前半段和后半段之和序列的N/2点DFT值，奇数点值X[2r+1]为输入x[n]的前半段和后半段之差再及WNnk相乘序列的N/2点DFT值。
令， (2.13)
(2.14)
则有：
(2.15)
这样，也可以用两个N/2点DFT来组合成一个N点DFT，组合过程如下图所示：

-1
图2.2 DIF-FFT算法
这样可以把一个N点DFT分解为两个N/2点DFT的组合，两个N/2点DFT还可以继续分解，设N=2M，则经过M-1次分解，最后可以分解成为N/2个两点DFT，可以由一个蝶形运算来求解。例如8点DIF-FFT蝶形运算图如图2.2

图2.3 8点DIF-FFT蝶形图
输出序列的排列规律不是从小到大按顺序的，而是按照倒叙规则排序的，即先将0-7转换为二进制数，然后将二进制数左右倒序，再转为十进制就可以得到新的数列，即：0，4，2，6，1，5，3，7。

2.5 按频率抽取的FFT的特点
2.5.1 原位运算
在DIF-FFT蝶形图中，取第m级且两输入节点分别在第k，j行的蝶形为例，讨论DIF-FFT的原位运算规律。由图可得蝶形运算的关系式可表示为=，=[]。有上式可得的m-1级的第k行及第j行的输出，在运算流图中的作用是，用来计算第m级的第k行和第j行的输出，，这样当计算完，后，，在运算流图中就不在起作用，因此可以采用原位运算，把，直接存入原来存放，的存储单元。同理可以把第m级蝶形的N个输出值直接存放在第m-1级蝶形输出的N个存储单元中，这样从第一级的输入x（n）开始到最后一级的输出X(k)，只需要N个存储单元。
2.5.2蝶形运算两节点之间的“距离”
第一级蝶形每个蝶形运算量节点的“距离”为4，第二级每个蝶形运算另节点的“距离”为2，第三级蝶形每个蝶形运算量节点的“距离”为1。依次类推：对于N等于2的L次方的DIF-FFT，可以得到第M级蝶形每个蝶形运算量节点的“距离”为2的L-M次方。
2.5.3旋转因子的变化规律
以8点的FFT为例，第一级蝶形，r=0,1，2；第二级蝶形，r=0,1；第三级的蝶形，r=0。依次类推，对于M级蝶形，旋转因子的指数为
r=J∙2M-1，J=0,1,2,3，……， (2.16)
这样就可以算出每一级的旋转因子。对于M级的任一蝶形运算所对应的旋转因子的指数，可以 如下方法得到：1将待求的蝶形输入节点中上面节点的行标号值k写成L位二进制数；2将此二进制数乘以2的M-1次方，即将L位二进制数左移M-1位，右边的空位补零，然后从低位到高位截取L位，即所得指数r所对应的二进制数。

3 程序设计
FFT程序包括变址（倒位序）和L级递推计算（N=，L为正整数）两部分。
3.1变址
DIF-FFT是输出倒位序的变址处理，设x(i)表示存放自然顺序输入数据的内存单元，x（j）表示存放倒位序序数的内存单元，I、J=0,1,…,N-1,当I=J时，不用变址；当IJ时，需要变址；但是当I<J时，进行变址在先，故在I>J时，就不需要再变址了，否则变址两次等于不变。其中本程序使用的“反向进位加法”。也可用bin2dec函数可以实现倒以位序。
3.2L级递推计算
整个L级递推过程由三个嵌套循环构成。外层的一个循环控制L（L=）级的顺序运算；内层的两个循环控制同一级（M相同）各蝶形结的运算，其中最内层循环控制同一种（即中的r相同）蝶形结的运算。其循环变量为I，I用来控制同一种蝶形结运算。其步进值为蝶形结的间距值LE=，同一种蝶形结中参加运算的两节点的间距为LE1=点。第二层循环，其循环变量J用来控制计算不同种（系数不同）的碟形结，J的步进值为1。也可以看出，最内层循环完成每级的蝶形结运算，第二层循环则完成系数的运算。最外层循环，用循环变量M来控制运算的级数，M为1到L，步进值为1，当M改变时，则LE1，LE和系数U都会改变。
MATLAB实现的代码：
function [Xk]=DIF_FFT(xn,N);
%实现对输入序列的DIF-FFT的基2算法，点数小于等于2的次幂
n=nextpow2(length(xn)); %取最小的满足2次幂的n
N=2^n;
if length(xn)<N
xn=[xn,zeros(1,N-length(xn))];
end
%蝶形运算算法
M=log2(N); %M=3 或 4
%第一层循环 分解步骤
for m=0:M-1
Num_Group=2^m; %每一次分解中组的个数
IntV_Group=N/2^m; %每一次分解中组之间的间距
IntV_Unit=N/2^(m+1); %每一组运算单元的间距
Count=IntV_Unit-1; %运算单元循环次数第三次循环次数
Wn=exp(-j*2*pi/IntV_Group);%旋转因子
%第二层循环组的循环
for g=1:Num_Group
dx1=(g-1)*IntV_Group; %第g组中运算量1的偏移量
dx2=(g-1)*IntV_Group+IntV_Unit; %第g组中运算量2的偏移量
%第三层循环组内运算的循环
for r=0:Count
k=r+1; %“组内”序列的下标
b=xn(k+dx1);
xn(k+dx1)=b+xn(k+dx2); %第m级，第g组的蝶形运算式1
xn(k+dx2)=[b-xn(k+dx2)]*Wn^r; %第m级，第g组的蝶形运算式2
end
end
end
%序列排序开始
n1=fliplr(dec2bin([0:N-1])); %码位倒置十进制转二进制并倒序
n2=bin2dec(n1); %倒序后二进制转十进制
%倒序循环输出
for i=1:N
Xk(i)=xn(n2(i)+1);
end
4结果及分析
编写主函数，在主函数中输入一个16点的序列分别调用自己编写的FFT函数，和MATLAB本身系统的FFT函数并比较两个结果是否相等，以判断自己编写的FFT程序是否正确
xn=[0:15];m=1:16；N=16
x1=DIF_FFT(xn,N)
x2=fft(xn)
x3=abs(x1);x4=abs(x2);
x5=angle(x1);x6=angle(x2);
figure('NumberTitle', 'off', 'Name', '电信1206 李嘉辛');
subplot(3,1,1)
stem(m,xn);title('输入的离散序列')
subplot(3,1,2)
stem(m,x3);title('经过DIF_FFT后得到的频谱的幅度')
subplot(3,1,3)
stem(m,x5);title('经过DIF_FFT后得到的频谱的相位')
figure('NumberTitle', 'off', 'Name', '电信1206 李嘉辛');
subplot(3,1,1)
stem(m,xn);title('输入的离散序列')
subplot(3,1,2)
stem(m,x4);title('经过库函数fft后得到的频谱的幅度')
subplot(3,1,3)
stem(m,x6);title('经过库函数fft后得到的频谱的相位')

图4.1 DIF-FFT结果图4.2库函数FFT结果
通过观察比较，得到的序列各点的值以及直观的通过图形，可以得到自己编写的DIF_FFT函数实现对序列进行FFT变换得到的结果及库函数FFT得到的结果是一样的。说明DIF_FFT子程序是正确的。从图中也可以看出有限长序列通过FFT后得到的频域为离散的。从理论讲，有限长序列经过离散傅里叶变换后，得到的频谱为离散的，从而也说明了FFT是DFT的优化方法，也属于DFT。
这个程序可以实现基2-FFT,但是如果想在运行时直接输入要变换的点数就不行，必须在调用FFT函数前现将要算的序列定义好，这是这个程序的不足之处。但是该程序可以计算不是2的整数次幂的序列。所以在主程序中，输入序列必须给出才能进行FFT变换。
当使用编写的程序进行16点的DIF-FFT计算时结果如下：
》xn=[1 2 3 4 5 6 7 8];N=8;DIF_FFT(xn,N)
>> xn=[0:15];N=16;DIF_FFT(xn,N)
ans =
1.0e+02 *
Columns 1 through 4
1.2000 + 0.0000i -0.0800 + 0.4022i -0.0800 + 0.1931i -0.0800 + 0.1197i
Columns 5 through 8
-0.0800 + 0.0800i -0.0800 + 0.0535i -0.0800 + 0.0331i -0.0800 + 0.0159i
Columns 9 through 12
-0.0800 + 0.0000i -0.0800 - 0.0159i -0.0800 - 0.0331i -0.0800 - 0.0535i
Columns 13 through 16
-0.0800 - 0.0800i -0.0800 - 0.1197i -0.0800 - 0.1931i -0.0800 - 0.4022i
当调用matlab自带的FFT程序进行相同的16点的FFT计算时结果如下：
>> xn=[0:15];fft(xn)
ans =
1.0e+02 *
Columns 1 through 4
1.2000 + 0.0000i -0.0800 + 0.4022i -0.0800 + 0.1931i -0.0800 + 0.1197i
Columns 5 through 8
-0.0800 + 0.0800i -0.0800 + 0.0535i -0.0800 + 0.0331i -0.0800 + 0.0159i
Columns 9 through 12
-0.0800 + 0.0000i -0.0800 - 0.0159i -0.0800 - 0.0331i -0.0800 - 0.0535i
Columns 13 through 16
-0.0800 - 0.0800i -0.0800 - 0.1197i -0.0800 - 0.1931i -0.0800 - 0.4022i
两者结果相同，故编写的程序正确。

5心得体会
通过做这次程序设计，我对MATLAB编程有了进一步的掌握，对数字信号处理在MATLAB中的实现有了更深的体会，对数字信号处理的的理论知识有了更深刻的认识，在学习基2FFT算法时有许多地方不理解比如如何进行倒位序，蝶形运算是怎么回事等，通过重新复习相关知识，编写程序对以前一些迷惑的地方理解的更透彻，学会了理论及实践相结合的方法。理论是实践的基础，只有掌握了相关的理论知识才能更好更轻松的实践。当理论某些细节不是很理解时，可以通过编程仿真来实现，将仿真结果及理论结合起来进行对比理解这样会容易点。不过这限于一些小的地方，当很多都不懂时，就不行了因为你很多地方不懂时就不能实现编程了。如果自己根据FFT的计算方法，计算一个序列经FFT变换后的结果，会花很长的实践但是直接调用已编好的程序，会轻松的得到结果。可见，使用程序的好处。在以后的学习中，可以多根据理论，在编程上去实现，这样也有助于学习。

参考文献
[1] 刘泉，阙大顺．数字信号处理原理及实现．北京：电子工业出版社，2005．
[2] 程佩青. 数字信号处理教程. 清华大学出版社. 2008. 5
[3] 余成波. 数字信号处理及其MATLAB实现. 清华大学出版社. 2008. 2
[4] 王洪元. MATLAB语言及其在电子信息工程中的应用. 清华大学出版社 . 2005.9
[5] 王嘉梅. 基于MATLAB的数字信号处理及实践开发. 西安电子科技大学出版社 2007.12

本科生课程设计成绩评定表

姓名
李嘉辛
性别
男
专业班级
电信1206
课程设计题目：8点基于DIF的FFT的实现
课程设计答辩或质疑记录：
成绩评定依据：

指导教师签字：_____________
2015年1 月20日