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

1数字信号处理简介

1.1数字信号处理定义

数字信号是用数字序列表示的信号，数字信号处理就是通过计算机或专用处理设备，用数值计算等数字方式对数字序列进行各种处理，将数字信号变换成符合要求的某种形式。数字信号处理主要包括数字滤波和数字频谱分析两大部分。例如，对数字信号进行滤波，其频带或滤除噪声和干扰，以提取和增强信号的有用分量；对信号进行频谱分析或功率分析，了解信号的频谱组成，以对信号进行识别。当然，凡是用数字方式对信号进行滤波，变换，增强，压缩，估计和识别等都是数字信号处理的研究范畴。

数字信号处理在理论上所涉及的范围及其广泛。数字领域中的微积分，概率统计，随机过程，高等代数，数值分析，复变函数和各种变换（如傅里叶变换，Z变换，离散傅里叶变换，小波变换等）都是它的基本工具，网络理论，信号与系统等则是它的理论基础。在科学发展上，数字信号处理又和最优控制，通信理论等紧密相连，目前已成为人工智能，模式识别，神经网络等新兴学科的重要理论基础，其实现技术又和计算机科学和微电子技术密不可分。因此，数字信号处理是把经典的理论基础体系作为自身的理论基础，同时又使自己成为一系列新兴学科的理论基础。

1.2数字信号处理的特点及实现方法

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

数字信号处理的主要研究对象是数字信号，且采用数值运算的方法达到处理的目的。数字信号处理的实现方法基本上可以分为软件实现方法、硬件实现方法和软硬件想结合的实现方法。数字信号处理的理论、算法和实现方法三者是密不可分的。

1

武汉理工大学《数字信号处理》课程设计报告书

1 FFT算法

1.1DFT的定义

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

N1X[k]通常令ej2Nx[n]en0j2Nnk，k=0,1,„N-1

WN，称为旋转因子。

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]=N2/2次复数乘法。即比直接计算少作一半乘法。因子（N/2）2表示直接计算（N/2）点DFT所需要的乘法次数，而乘数2代表必须完成两个DFT。上述处理方法可以反复使用，即（N/2）点的DFT计算也可以化成两个（N/4）点的DFT（假定N/2为偶数），从而又少作一半的乘法。这样一级一级的划分下去一直到最后就划分成两点的FFT运算的情况。

1.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] （r=0,1,„,N/2-1）

DFT化为：

2

武汉理工大学《数字信号处理》课程设计报告书

N1N/21N/21X[k]DFT{x[n]}N/21n0x[n]WnkN2rkr0x[2r]W2rkNr0x[2r1]WN(2r1)kN/21r0N/21x1[r]Wx1[r]W2rkNWWkNr0N/21x2[r]WN

r0rkN/2kNr0x2[r]WN/22rkrk上式中利用了旋转因子的可约性，即：WNN/21N/21WN/2。又令

rkX1[k]r0x1[r]WrkN/2,X2[k]r0x2[r]WN/2rk，则上式可以写成：

X[k]X1[k]WNX2[k] (k=0,1,„,N/2-1)

k可以看出，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/2则后半段的DFT值表达式：

X1[N2N/21rkr(kN/2)，

k]r0x1[r]W2N/2r(Nk)N/21r0x1[r]WN/2X1[k]

rkX2[N2k]X2[k] （k=0,1,„,N/2-1）

所以后半段（k=N/2,„,N-1）的DFT值可以用前半段k值表达式获得，中间还利用到

WN2(Nk)NWN2WkWk，得到后半段的X[k]值表达式为：X[k]X1[k]WNkX2[k]（k=0,1,„,N/2-1）。

这样，通过计算两个N/2点序列x1[n],x2[n]的N/2点DFTX1[k],X2[k]，可以组合得到N点序列的DFT值X[k]，其组合过程如下图所示：

X1[k] X1[k]WNkX2[k]

X2[k] WNnk -1 X1[k]WNkX2[k]

1.4基2按频率抽取（DIF）的FFT算法

3

武汉理工大学《数字信号处理》课程设计报告书

设序列长度为N2L，L为整数（如果序列长度不满足此条件，通过在后面补零让其满足）。

在把X[k]按k的奇偶分组之前，把输入按n的顺序分成前后两半：

N1N/21nkNN1X[k]DFT{x[n]}N/21N/21x[n]Wn0(nn0N2)kx[n]WnkNnN/2x[n]WNnkn0N/21x[n]WnkNn0x[nNkN2]WNnk

Nn0[x[n]x[nN2NkN2]WNk2]WN,k0,1,...,N1因为W2N1，则有W(1)，所以：

[x[n](1)x[nkN/21X[k]n0N2]]WN,k0,1,...,N1

nk按k的奇偶来讨论，k为偶数时：

N/21X[2r]n0[x[n]x[nN2N2]]WN,k0,1,...,N1

2rnN/21k为奇数时：X[2r1]前面已经推导过WN2rkn0[x[n]x[n]]WN(2r1)n,k0,1,...,N1

WN/2，所以上面的两个等式可以写为：

N/21rkX[2r]n0[x[n]x[nN2N2]]WN/2,r0,1,...,N/21

rnN/21X[2r1]n0{[x[n]x[n]]WN}WN/2,r0,1,...,N/21

nnr通过上面的推导，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1nN2]，x2[n][x[n]x[nN/21N2]]WN，则有：

nX[2r]

n0x1[n]WrnN/2,X[2r1]n0x2[n]W4

rnN/2,r0,1,...,N21

武汉理工大学《数字信号处理》课程设计报告书

这样，也可以用两个N/2点DFT来组合成一个N点DFT，组合过程如下图所示：

x[n] x[n]x[nN2]

x[nN2] -1 WN [x[n]x[nn

N2]]WN

n1.5 按频率抽取的FFT的特点

1．原位运算

在DIF-FFT蝶形图中，取第m级且两输入节点分别在第k，j行的蝶形为例，讨论DIF-FFT的原位运算规律。由图可得蝶形运算的关系式可表示为Xm(k)=Xm1(k)Xm1(j)，

Xm(j)=[Xm1(k)Xm1(j)]WN 。有上式可得的m-1级的第k行与第j行的输出Xm1(k)，

r用来计算第m级的第k行和第j行的输出Xm(k)，Xm(j)，Xm1(j)在运算流图中的作用是，

这样当计算完Xm(k)，Xm(j)后，Xm1(k)，Xm1(j)在运算流图中就不在起作用，因此可以采用原位运算，把Xm(k)，Xm(j)直接存入原来存放Xm1(k)，Xm1(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=J2M1，J=0,1,2,3，„„，2M11

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

5

武汉理工大学《数字信号处理》课程设计报告书

DIF程序编写

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

DIF-FFT是输出倒位序的变址处理，设x(i)表示存放自然顺序输入数据的内存单元，x（j）表示存放倒位序序数的内存单元，I、J=0,1,…,N-1,当I=J时，不用变址；当I J时，需要变址；但是当i整个L级递推过程由三个嵌套循环构成。外层的一个循环控制L（L=log2N）级的顺序运算；内层的两个循环控制同一级（M相同）各蝶形结的运算，其中最内层循环控制同一种（即

WNr中的r相同）蝶形结的运算。其循环变量为I，I用来控制同一种蝶形结运算。

其步进值为蝶形结的间距值LE=2M，同一种蝶形结中参加运算的两节点的间距为LE1=2ML点。第二层循环，其循环变量J用来控制计算不同种（系数不同）的碟形结，J

k的步进值为1。也可以看出，最内层循环完成每级的蝶形结运算，第二层循环则完成系数WN的运算。最外层循环，用循环变量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)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; %每一级中组与组之间的间距

6

武汉理工大学《数字信号处理》课程设计报告书

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

7

武汉理工大学《数字信号处理》课程设计报告书

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

8

武汉理工大学《数字信号处理》课程设计报告书

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变换后的结果，会花很长的实践但是直接调用已编好的程序，会轻松的得到结果。可见，使用程序的好处。在以后的学习中，可以多根据理论，在编程上去实现，这样也有助于学习。

9

[1] 程佩青. [2] 刘卫国. [3] 余成波.[4] 陈怀琛.[5] 王洪元.[6] 王嘉梅.

武汉理工大学《数字信号处理》课程设计报告书

参考文献

数字信号处理教程. 清华大学出版社. 2008. 5

MATLAB程序设计教程 . 中国水利水电出版社 . 2008. 6 数字信号处理及其MATLAB实现. 清华大学出版社. 2008. 2

数字信号处理教程：MATLAB释义及实现. 电子工业出版社 . 2006.6 MATLAB语言及其在电子信息工程中的应用. 清华大学出版社 . 2005.9 基于MATLAB的数字信号处理及实践开发. 西安电子科技大学出版社 2007.1210