8点基于DIF-FFT课程设计(共14页)

1、武汉理工大学信号分析与处理课程设计1数字信号处理简介1.1 数字信号处理概念相对于模拟信号，数字信号是时间、幅值都离散的信号，数字信号处理就是通过计算机或专用处理设备，用数值计算等数字方式对数字信号进行各种处理，将数字信号变换成符合要求的某种形式。数字信号处理主要包括数字滤波和数字频谱分析两大部分。例如，对数字信号进行滤波，限制其频带或滤除噪声和干扰，以提取和增强信号的有用分量；对信号进行频谱分析或功率分析，了解信号的频谱组成，以对信号进行识别。当然，凡是用数字方式对信号进行滤波，变换，增强，压缩，估计和识别等都是数字信号处理的研究范畴。数字信号处理在理论上所涉及的范围极其广泛。数字领域中的微

2、积分，概率统计，随机过程，高等代数，数值分析，复变函数和各种变换（如傅里叶变换，Z变换，离散傅里叶变换，小波变换等）都是它的基本工具，网络理论，信号与系统等则是它的理论基础。在科学发展上，数字信号处理又和最优控制，通信理论等紧密相连，目前已成为人工智能，模式识别，神经网络等新兴学科的重要理论基础，其实现技术又和计算机科学和微电子技术密不可分。因此，数字信号处理是把经典的理论基础体系作为自身的理论基础，同时又使自己成为一系列新兴学科的理论基础。1.2数字信号处理的特点及实现方法与模拟信号处理相比，数字信号处理具有高精度、高稳定性、灵活性好、易于大规模集成等显着的优点。数字信号处理的主要研究对象是

3、数字信号，且采用数值运算的方法达到处理的目的。数字信号处理的实现方法基本上可以分为软件实现方法、硬件实现方法和软硬件相结合的实现方法。数字信号处理的理论、算法和实现方法三者是密不可分的。 2 FFT算法介绍2.1 DFT的定义对于有限长离散数字信号xn，0 n N-1，其离散谱Xk可以由离散傅氏变换（DFT）求得。DFT的定义为：，k=0,1,N-1通常令，称为旋转因子。2.2 FFT基本思想由DFT的定义可以看出，在xn为复数序列的情况下，直接运算N点DFT需要N的二次方次复数乘法和N（N-1）次加法。因此，对于一些相当大的N值（如1024）来说，直接计算它的DFT所作的计算量是很大的，对实

4、现快速计算是很不利的。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运算的情况。虽然这样做，要求序列的

5、长度为（n=0,1,）,但在原序列后面补零使其满足条件长度，并不影响对原序列的频域分析，仅仅只是抽样点更密了一些。DIT-FFT和DIF-FFT正是基于这种思想。2 DIF-FFT2.1 DIF-FFT过程推导设序列长度为，L为整数（如果序列长度不满足此条件，通过在后面补零让其满足）。在把按k的奇偶分组之前，把输入按n的顺序分成前后两半：因为，则有，所以：按k的奇偶来讨论，k为偶数时：k为奇数时：前面已经推导过，所以上面的两个等式可以写为：通过上面的推导，的偶数点值和奇数点值分别可以由组合而成的N/2点的序列来求得，其中偶数点值为输入xn的前半段和后半段之和序列的N/2点DFT值，奇数点值为输

6、入xn的前半段和后半段之差再与相乘序列的N/2点DFT值。令，则有：这样，也可以用两个N/2点DFT来组合成一个N点DFT，组合过程如图2-1-1所示： -1 图2-1-1由于，N/2还是偶数，所以可以把N/2点的DFT奇、偶分组，这样就把一个N/2点的DFT分解成两个N/4点的DFT。当N=8时，采用DIF-FFT运算二次分解的运算流图如图2-1-2所示：图2-1-2 8点DIF-FFT流程图2.5 DIF-FFT的特点2.5.1 原位运算在DIF-FFT蝶形图中，取第m级且两输入节点分别在第k，j行的蝶形为例，讨论DIF-FFT的原位运算规律。由图可得蝶形运算的关系式可表示为=，= 。有上

7、式可得的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级蝶形每个蝶形运

8、算量节点的“距离”为2的L-M次方。 2.5.3 旋转因子的变化规律以8点的FFT为例，第一级蝶形，r=0,1，2；第二级蝶形，r=0,1；第三级的蝶形，r=0。依次类推，对于M级蝶形，旋转因子的指数为r=，J=0,1,2,3，这样就可以算出每一级的旋转因子。对于M级的任一蝶形运算所对应的旋转因子的指数，可以 如下方法得到：1将待求的蝶形输入节点中上面节点的行标号值k写成L位二进制数；2将此二进制数乘以2的M-1次方，即将L位二进制数左移M-1位，右边的空位补零，然后从低位到高位截取L位，即所得指数r所对应的二进制数。3 DIF程序编写与验证3.1 DIF程序编写FFT程序包括变址（倒位序）和

9、L级递推计算（N=，L为正整数）两部分。3.1.1 变址DIF-FFT是输出倒位序的变址处理，设x(i)表示存放自然顺序输入数据的内存单元，x（j）表示存放倒位序序数的内存单元，I、J=0,1,N-1,当I=J时，不用变址；当I J时，需要变址；但是当ij时，进行变址在先，故在IJ时，就不需要再变址了，否则变址两次等于不变。其中本程序使用的“反向进位加法”。也可以用bin2dec函数可以实现倒位序。3.1.2 L级递推计算整个L级递推过程由三个嵌套循环构成。外层的一个循环控制L（L=）级的顺序运算；内层的两个循环控制同一级（M相同）各蝶形结的运算，其中最内层循环控制同一种（即中的r相同）蝶形结

10、的运算。其循环变量为I，I用来控制同一种蝶形结运算。其步进值为蝶形结的间距值LE=，同一种蝶形结中参加运算的两节点的间距为LE1=点。第二层循环，其循环变量J用来控制计算不同种（系数不同）的碟形结，J的步进值为1。也可以看出，最内层循环完成每级的蝶形结运算，第二层循环则完成系数的运算。最外层循环，用循环变量M来控制运算的级数，M为1到L，步进值为1，当M改变时，则LE1，LE和系数U都会改变。3.2 程序编写流图 可按如下程序流图编写程序： 开始 原序列补零至 计算序列级数M= m=0 否mM-1 是本级中蝶形图个数group= 每个蝶形图中元素个数uint=8/group 旋转因子Wn=ex

11、p（-j*2*pi/group）利用组内计算 计算完毕 MATLAB实现的代码：1、在M序列中编写DIF-FFT函数：function Xk=DIF_FFT(xn,N);%本程序对输入序列实现DIF-FFT基2算法，%点数取小于等于长度的2的幂次 n=log2(2nextpow2(length(xn); %求的x长度对应的2的最低幂次nN=2n; 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=

12、2m; %每一级中组的个数 Interval_of_Group=N/2m; %每一级中组与组之间的间距 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+

13、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)*Wnr; %第m级，第g组的蝶形运算式2 end endend%序列排序开始n1=fliplr(dec2bin(0:N-1);%码位倒置步骤1：将码位转换为二进制

14、，再进行倒序n2=bin2dec(n1); %码位倒置步骤2：将码位转换为十进制后翻转for i=1:N Xk(i)=xn(n2(i)+1); %根据码位倒置的结果，将xn重新排序,存入Xk中end3.2 程序正确认证编写主函数，在主函数中输入一个8点序列分别调用自己编写的DIF-FFT函数，和MATLAB本身系统的FFT函数,并比较两个结果是否相等，以判断自己编写的FFT程序是否正确xn=1:8;m=1:8；N=8x1=DIF_FFT(xn,N)x2=fft(xn)x3=abs(x1);x4=abs(x2);x5=angle(x1);x6=angle(x2);subplot(3,1,1)st

15、em(m,xn),grid;title(输入的离散序列)subplot(3,1,2)stem(m,x3),grid;title(经过DIF_FFT后得到的幅度谱)subplot(3,1,3)stem(m,x5),grid;title(经过DIF-FFT后得到的相位谱)figure (2)subplot(3,1,1)stem(m,xn),grid;title(输入的离散序列)subplot(3,1,2)stem(m,x4),grid;title(调用库函数fft后得到的幅度谱)subplot(3,1,3)stem(m,x6),grid;title(调用库函数fft后得到的相位谱)通过观察比较，得

16、到的序列各点的值以及直观的通过图形，可以得到自己编写的DIF_FFT函数实现对序列进行FFT变换得到的结果与库函数FFT得到的结果是一样的。说明DIF_FFT子程序是正确的。从图中也可以看出有限长序列通过FFT后得到的频域为离散的。从理论讲，有限长序列经过离散傅里叶变换后，得到的频谱为离散的，从而也说明了FFT是DFT的优化方法，也属于DFT。这个程序可以实现基2FFT,但是如果想在运行时直接输入要变换的点数就不行，必须在调用FFT函数前现将要算的序列定义好，这是这个程序的不足之处。但是该程序可以计算不是2的整数次幂的序列。所以在主程序中，输入序列必须给出才能进行FFT变换。当使用编写的程序进

17、行8点的DIF-FFT计算时结果如下：xn=1 2 3 4 5 6 7 8;N=8;DIF_FFT(xn,N)x1 = Columns 1 through 7 36.0000 -4.0000 + 9.6569i -4.0000 + 4.0000i -4.0000 + 1.6569i -4.0000 -4.0000 - 1.6569i -4.0000 - 4.0000i Column 8 -4.0000 - 9.6569i当调用matlab自带的FFT程序进行相同的8点的FFT计算时结果如下：xn=1 2 3 4 5 6 7 8;fft(xn)ans = Columns 1 through 7

18、36.0000 -4.0000 + 9.6569i -4.0000 + 4.0000i -4.0000 + 1.6569i -4.0000 -4.0000 - 1.6569i -4.0000 - 4.0000i Column 8 -4.0000 - 9.6569i两者结果相同，故编写的程序正确。只需将上述程序中“N=8” 改为“N=16”即可求得16点的FFT，程序设计如下：xn=0:15;m=1:16;N=16x1=DIF_FFT(xn,N)x1 =Columns 1 through 7 1.2000 -0.0800 + 0.4022i -0.0800 + 0.1931i -0.0800 +

19、 0.1197i -0.0800 + 0.0800i -0.0800 + 0.0535i -0.0800 + 0.0331iColumns 8 through 14 -0.0800 + 0.0159i -0.0800 -0.0800 - 0.0159i -0.0800 - 0.0331i -0.0800 - 0.0535i -0.0800 - 0.0800i -0.0800 - 0.1197iColumns 15 through 16 -0.0800 - 0.1931i -0.0800 - 0.4022i当调用matlab自带的FFT程序进行相同的8点的FFT计算时结果如下：ans =Colu

20、mns 1 through 7 1.2000 -0.0800 + 0.4022i -0.0800 + 0.1931i -0.0800 + 0.1197i -0.0800 + 0.0800i -0.0800 + 0.0535i -0.0800 + 0.0331iColumns 8 through 14 -0.0800 + 0.0159i -0.0800 -0.0800 - 0.0159i -0.0800 - 0.0331i -0.0800 - 0.0535i -0.0800 - 0.0800i -0.0800 - 0.1197iColumns 15 through 16 -0.0800 - 0.1931i -0.0800 - 0.4022i显然程序执行结果一致。将两者图形比较如下所示：图3-2-1由fft所得的函数谱图形 图3-2-1由DIF