8点DIF-FFT程序编写课程设计

1数字信号处理简介
1.1数字信号处理定义
数字信号是用数字序列表示的信号，数字信号处理就是通过计算机或专用处理设备，用数值计算等数字方式对数字序列进行各种处理，将数字信号变换成符合要求的某种形式。

数字信号处理主要包括数字滤波和数字频谱分析两大部分。

例如，对数字信号进行滤波，限制其频带或滤除噪声和干扰，以提取和增强信号的有用分量；对信号进行频谱分析或功率分析，了解信号的频谱组成，以对信号进行识别。

当然，凡是用数字方式对信号进行滤波，变换，增强，压缩，估计和识别等都是数字信号处理的研究范畴。

数字信号处理在理论上所涉及的范围及其广泛。

数字领域中的微积分，概率统计，随机过程，高等代数，数值分析，复变函数和各种变换（如傅里叶变换，Z变换，离散傅里叶变换，小波变换等）都是它的基本工具，网络理论，信号与系统等则是它的理论基础。

在科学发展上，数字信号处理又和最优控制，通信理论等紧密相连，目前已成为人工智能，模式识别，神经网络等新兴学科的重要理论基础，其实现技术又和计算机科学和微电子技术密不可分。

因此，数字信号处理是把经典的理论基础体系作为自身的理论基础，同时又使自己成为一系列新兴学科的理论基础。

1.2数字信号处理的特点及实现方法
与模拟信号处理相比，数字信号处理具有高精度、高稳定性、灵活性好、易于大规模集成等显著的优点。

数字信号处理的主要研究对象是数字信号，且采用数值运算的方法达到处理的目的。

数字信号处理的实现方法基本上可以分为软件实现方法、硬件实现方法和软硬件想结合的实现方法。

数字信号处理的理论、算法和实现方法三者是密不可分的。

1 FFT 算法
1.1DFT 的定义
对于有限长离散数字信号{x[n]}，0 ≤ n ≤ N-1，其离散谱{X[k]}可以由离散付氏变换（DFT ）求得。

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

1.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 运算的情况。

1.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-1）
DFT 化为：
1/21/21
2(21)0
/21
/21
221200/21
/21
1/22/2
[]{[]}[][2][21][][][][]N N N nk rk r k
N
N
N
n r r N N rk k rk
N
N
N
r r N N rk k rk
N N
N r r X k D FT x n x n W
x r W
x r W x r W W x r W x r W
W
x r W ---+===--==--====
=
+
+=+=
+∑
∑
∑
∑
∑
∑
∑
上式中利用了旋转因子的可约性，即：2/2rk
rk
N
N W W =。

又令
/21
/21
11/2
22/2
[][],[][]N N rk rk
N N r r X k x r W
X k x r W --===
=
∑
∑
，则上式可以写成：
12[][][]k
N 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/2rk r k N N N W W +=，
则后半段的DFT 值表达式：
/21
/21
(
)
2
11/2
1/210
[
][][][]2
N
N N r k rk
N N r r N X k x r W x r W X k --+==+=
=
=∑
∑
22[
][]2
N X k X k += （k=0,1,…,N/2-1）
所以后半段（k=N/2,…,N-1）的DFT 值可以用前半段k 值表达式获得，中间还利用到
()
2
2N N
k k
k
N N W W W
W
+==-，得到后半段的[]X k 值表达式为：12[][][]
k N 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[][]k N X k W X k +
2[]X k nk N W -1 12[][]k N X k W X k -
1.4基2按频率抽取（DIF ）的FFT 算法
设序列长度为2L N =，L 为整数（如果序列长度不满足此条件，通过在后面补零让其满足）。

在把[]X k 按k 的奇偶分组之前，把输入按n 的顺序分成前后两半：
1
/21
1
0/2
/21
/21
()200
/21
2
[]{[]}[][][][][]2
[[][]],0,1,...,1
2
N N N nk
nk nk
N
N
N
n n n N N N N n k
nk N
N
n n N N k
nk
N
N n X k D FT x n x n W
x n W
x n W N x n W
x n W N x n x n W W k N ---===--+
==-===
=
+
=
+
+
=
++
=-∑∑
∑
∑
∑
∑
因为2
1N
N
W
=-，则有2(1)N k
k
N
W =-，所以：
/21
[][[](1)[]],0,1,...,12
N k
nk
N n N X k x n x n W k N -===
+-+
=-∑
按k 的奇偶来讨论，k 为偶数时：
/21
20
[2][[][]],0,1,...,12
N rn
N n N X r x n x n W k N -==
++
=-∑
k 为奇数时：/21
(21)0
[21][[][]],0,1,...,12
N r n
N
n N X r x n x n W k N -+=+=-+
=-∑
前面已经推导过2/2rk
rk
N
N W W =，所以上面的两个等式可以写为：
/21
/20
[2][[][]],0,1,...,/212N rn
N n N X r x n x n W r N -==
++
=-∑
/21
/20
[21]{[[][]]},0,1,...,/212
N 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]的前半段和后半段之差再与n
N W 相乘序列的N/2点DFT 值。

令1[][][]2
N x n x n x n =++
，2[][[][]]2
n
N N x n x n x n W =-+
，则有：
/21
/21
1/2
2/2
[2][],[21][],0,1,...,
12
N N rn
rn
N N n n N X r x n W
X r x n W
r --===
+=
=-∑
∑
这样，也可以用两个N/2点DFT 来组合成一个N 点DFT ，组合过程如下图所示：
[]x n [][]2N x n x n ++
[]2
N x n +
-1 n
N W [[][]]2
n
N N x n x n W -+
1.5 按频率抽取的FFT 的特点
1．原位运算
在DIF-FFT 蝶形图中，取第m 级且两输入节点分别在第k ，j 行的蝶形为例，讨论DIF-FFT 的原位运算规律。

由图可得蝶形运算的关系式可表示为()m X k =11()()m m X k X j ---，
()m X j =[11()()m m X k X j ---]r
N W 。

有上式可得的m-1级的第k 行与第j 行的输出1()m X k -，
1()m X j -在运算流图中的作用是，
用来计算第m 级的第k 行和第j 行的输出()m X k ，()m X j ，这样当计算完()m X k ，()m X j 后，1()m X k -，1()m X j -在运算流图中就不在起作用，因此可以采用原位运算，把()m X k ，()m X j 直接存入原来存放1()m X k -，1()m X j -的存储单元。

同理可以把第m 级蝶形的N 个输出值直接存放在第m-1级蝶形输出的N 个存储单元中，这样从第一级的输入x （n ）开始到最后一级的输出X(k)，只需要N 个存储单元。

2．蝶形运算两节点之间的“距离”
第一级蝶形每个蝶形运算量节点的“距离”为4，第二级每个蝶形运算另节点的“距离”为2，第三级蝶形每个蝶形运算量节点的“距离”为1。

依次类推：对于N 等于2的L 次方的DIF-FFT ，可以得到第M 级蝶形每个蝶形运算量节点的“距离”为2的L-M 次方。

3．旋转因子的变化规律
以8点的FFT 为例，第一级蝶形，r=0,1，2；第二级蝶形，r=0,1；第三级的蝶形，r=0。

依次类推，对于M 级蝶形，旋转因子的指数为
r=12M J - ，J=0,1,2,3，……，121M --
这样就可以算出每一级的旋转因子。

对于M 级的任一蝶形运算所对应的旋转因子的指数，可以 如下方法得到：1将待求的蝶形输入节点中上面节点的行标号值k 写成L 位二进制数；2将此二进制数乘以2的M-1次方，即将L 位二进制数左移M-1位，右边的空位补零，然后从低位到高位截取L 位，即所得指数r 所对应的二进制数。

DIF程序编写
FFT程序包括变址（倒位序）和L级递推计算（N=L2，L为正整数）两部分。

1．变址
DIF-FFT是输出倒位序的变址处理，设x(i)表示存放自然顺序输入数据的内存单元，x （j）表示存放倒位序序数的内存单元，I、J=0,1,…,N-1,当I=J时，不用变址；当I ≠J时，需要变址；但是当i<j时，进行变址在先，故在I<J时，就不需要再变址了，否则变址两次等于不变。

其中本程序使用的“反向进位加法”。

也可以用bin2dec函数可以实现倒位序。

2．L级递推计算
整个L级递推过程由三个嵌套循环构成。

外层的一个循环控制L（L=N
2
log）级的顺序运算；内层的两个循环控制同一级（M相同）各蝶形结的运算，其中最内层循环控制同
一种（即
r
N
W中的r相同）蝶形结的运算。

其循环变量为I，I用来控制同一种蝶形结运算。

其步进值为蝶形结的间距值LE=M2，同一种蝶形结中参加运算的两节点的间距为LE1=2M L-点。

第二层循环，其循环变量J用来控制计算不同种（系数不同）的碟形结，J 的步进值为1。

也可以看出，最内层循环完成每级的蝶形结运算，第二层循环则完成系数k
N
W 的运算。

最外层循环，用循环变量M来控制运算的级数，M为1到L，步进值为1，当M 改变时，则LE1，LE和系数U都会改变。

MATLAB实现的代码：
function [Xk]=DIF_FFT_2(xn,N);
%本程序对输入序列实现DIF-FFT基2算法，点数取小于等于长度的2的幂次
N=8;
n=log2(2^nextpow2(length(xn))); %求的x长度对应的2的最低幂次n
N=2^n;
if length(xn)<N
xn=[xn,zeros(1,N-length(xn))]; %若长度不是2的幂，补0到2的整数幂end
%蝶形运算开始
M=log2(N); %“级”的数量
for m=0:M-1 %“级”循环开始
Num_of_Group=2^m; %每一级中组的个数
Interval_of_Group=N/2^m; %每一级中组与组之间的间距
Interval_of_Unit=N/2^(m+1); %每一组中相关运算单元之间的间距
Cycle_Count= Interval_of_Unit -1; %每一组内的循环次数
Wn=exp(-j*2*pi/Interval_of_Group);%旋转因子
for g=1:Num_of_Group %“组”循环开始
Interval_1=(g-1)*Interval_of_Group;
%第g组中蝶形运算变量1的偏移量
Interval_2=(g-1)*Interval_of_Group+Interval_of_Unit;
%第g组中蝶形运算变量2的偏移量
for r=0:Cycle_Count; %“组内”循环开始
k=r+1; %“组内”序列的下标
xn(k+Interval_1)=xn(k+Interval_1)+xn(k+Interval_2);
%第m级，第g组的蝶形运算式1
xn(k+Interval_2)=[xn(k+Interval_1)-xn(k+Interval_2)-xn(k+Interval_2)]*Wn^r;
%第m级，第g组的蝶形运算式2 end
end
end
%序列排序开始
n1=fliplr(dec2bin([0:N-1]));
%码位倒置步骤1：将码位转换为二进制，再进行倒序
n2=[bin2dec(n1)];
%码位倒置步骤2：将码位转换为十进制后翻转
for i=1:N
Xk(i)=xn(n2(i)+1); %根据码位倒置的结果，将xn重新排序,存入Xk中
End
程序正确认证：
编写主函数，在主函数中输入一个序列分别调用自己编写的FFT函数，和MATLAB本身系统的FFT函数并比较两个结果是否相等，以判断自己编写的FFT程序是否正确
xn=[0:7];m=1:8；N=8
x1=DIF_FFT(xn,N)
x2=fft(xn)
x3=abs(x1);x4=abs(x2);
x5=angle(x1);x6=angle(x2);
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 (2)
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后得到的频谱的相位')
通过观察比较，得到的序列各点的值以及直观的通过图形，可以得到自己编写的DIF_FFT函数实现对序列进行FFT变换得到的结果与库函数FFT得到的结果是一样的。

说明DIF_FFT子程序是正确的。

从图中也可以看出有限长序列通过FFT后得到的频域为离散的。

从理论讲，有限长序列经过离散傅里叶变换后，得到的频谱为离散的，从而也说明了FFT是DFT的优化方法，也属于DFT。

这个程序可以实现基2FFT,但是如果想在运行时直接输入要变换的点数就不行，必须在调用FFT函数前现将要算的序列定义好，这是这个程序的不足之处。

但是该程序可以计算不是2的整数次幂的序列。

所以在主程序中，输入序列必须给出才能进行FFT变换。

当使用编写的程序进行8点的DIF-FFT计算时结果如下：
》xn=[1 2 3 4 5 6 7 8];N=8;DIF_FFT(xn,N)
Ans=
Columns 1 through 6
36.0-4.0000+9.6569i -4.0000+4.0000i -4.0000-1.6569i
Columns 7 through 8
-4.0000-4.0000i -4.0000-9.6569i
当调用matlab自带的FFT程序进行相同的8点的FFT计算时结果如下：
》xn=[1 2 3 4 5 6 7 8];fft[xn]
Ans=
Columns 1 through 6
37.0-4.0000+9.6569i -4.0000+4.0000i -4.0000-1.6569i
Columns 7 through 8
-4.0000-4.0000i -4.0000-9.6569i
两者结果相同，故编写的程序正确。

体会
通过做这次程序设计，我对MATLAB编程有了进一步的掌握，对数字信号处理在MATLAB中的实现有了更深的体会，对数字信号处理的的理论知识有了更深刻的认识，在学习基2FFT算法时有许多地方不理解比如如何进行倒位序，蝶形运算是怎么回事等，通过重新复习相关知识，编写程序对以前一些迷惑的地方理解的更透彻，学会了理论与实践相结合的方法。

理论是实践的基础，只有掌握了相关的理论知识才能更好更轻松的实践。

当理论某些细节不是很理解时，可以通过编程仿真来实现，将仿真结果与理论结合起来进行对比理解这样会容易点。

不过这限于一些小的地方，当很多都不懂时，就不行了因为你很多地方不懂时就不能实现编程了。

如果自己根据FFT的计算方法，计算一个序列经FFT 变换后的结果，会花很长的实践但是直接调用已编好的程序，会轻松的得到结果。

可见，使用程序的好处。

在以后的学习中，可以多根据理论，在编程上去实现，这样也有助于学习。

武汉理工大学《数字信号处理》课程设计报告书
参考文献
[1] 程佩青. 数字信号处理教程.清华大学出版社. 2008. 5
[2] 刘卫国. MA TLAB程序设计教程.中国水利水电出版社. 2008. 6
[3] 余成波.数字信号处理及其MA TLAB实现.清华大学出版社. 2008. 2
[4] 陈怀琛.数字信号处理教程：MATLAB释义及实现.电子工业出版社 . 2006.6
[5] 王洪元.MATLAB语言及其在电子信息工程中的应用.清华大学出版社. 2005.9
[6] 王嘉梅.基于MA TLAB的数字信号处理及实践开发.西安电子科技大学出版社2007.12。