数字信号处理课程设计-DFT与FFT计算速度比较分析.doc

燕山大学课 程 设 计 说 明 书题目： dft与fft计算速度比较分析 学院（系）： 电气工程学院 年级专业： 07级仪表3班 学 号： 学生姓名： 指导教师：教师职称： 完成时间： 2010年7月15日 燕山大学课程设计（论文）任务书院（系）： 基层教学单位： 学 号学生姓名专业（班级）07级仪表3班设计题目dft与fft计算速度比较分析设计技术参数用matlab实现dft与fftdft与fft的运算时间设计要求设计dft与fft程序，比较两种频谱分析方法与计算速度，并与理论值进行比较，分析误差产生原因。工作量先对两种算法进行介绍，包括推导过程及运算性质，然后用matlab实现两种算法，再分别对两种算法进行运算时间对比，并分析时间长短的原因。工作计划第一周周一至周三：去图书馆，上网查阅相关资料，了解相关知识。第一周周四至周五：利用matlab编写代码，实现所要求的程序。第二周周一至周二：调试代码，并写说明书第二周周三至周四：完善，修改说明书参考资料信号处理原理及应用，谢平等编著，机械工业出版社dft与fft在实用时的性能比较，邱宽明，赵胜凯matlab 6.x信号处理，邹鲲等编著，清华大学出版社数字信号处理基础及matlab实现，周辉, 董正宏编著，中国林业出版社指导教师签字基层教学单位主任签字说明：此表一式四份，学生、指导教师、基层教学单位、系部各一份。2010年 7 月 5日 燕山大学课程设计评审意见表指导教师评语：成绩： 指导教师： 2010年 7 月 16 日答辩小组评语： 成绩： 评阅人： 2010年 7 月 16 日课程设计总成绩：答辩小组成员签字：2010年 7 月 16 日 燕 山 大 学 课 程 设 计 说 明 书摘 要本说明书主要是在介绍两种用于信号处理的傅里叶变换算法dft（离散傅里叶变换）和fft（快速傅里叶变换），分别介绍了这两种运算的推导过程，并且对这两种变换作了简要的介绍，分析了各自的性质。然后通过matlab分别实现了这两种傅里叶变换，并对这两种变换进行了运算时间的比较分别对同一函数进行dft和fft计算两者的运行时间，并作图比较。本说明书的程序部分都是在matlab环境下进行的运算。matlab是矩阵实验室（matrix laboratory）的简称，是美国mathworks公司出品的商业数学软件，用于算法开发、数据可视化、数据分析以及数值计算的高级技术计算语言和交互式环境，主要包括matlab和simulink两大部分。matlab的基本数据单位是矩阵，它的指令表达式与数学、工程中常用的形式十分相似，故用matlab来解算问题要比用c，fortran等语言完成相同的事情简捷得多。在新的版本中也加入了对c，fortran，c+ ，java的支持，可以直接调用,用户也可以将自己编写的实用程序导入到matlab函数库中方便自己以后调用。通过本说明书，可以帮助读者进一步了解dft以及fft的一些相关知识，并能够学到一些matlab的知识。目 录第一章 离散傅里叶变换(dft)31.离散傅里叶变换32. 用dft对信号进行谱分析3第二章 离散傅里叶变换的快速算法(fft)31. 直接计算dft的问题和改善dft运算效率的基本途径32. 基2 fft算法33. 按时间抽取得fft算法的特点34. 按频率抽取fft算法(dif-fft)3第三章 用matlab实现dft和fft以及对两者运算时间的分析31.两种算法的时间计算32.分析两者运算时间的差异：3第四章 总结3第一章 离散傅里叶变换(dft)1.离散傅里叶变换所谓傅里叶变换就是在以时间为自变量的“信号”与以频率为自变量的“频谱”函数之间的某种变换关系。随时间自变量形式的不同，其傅里叶变换的形式也有不同，常用的两种傅里叶变换：周期序列的离散傅里叶级数(dfs)和非周期序列的傅里叶变换(dtft)，其表示式分别为： (1.1) (1.2)在实际工作中，当用数字计算机对信号进行频谱分析时，要求信号必须以有限长度的离散值作为输入，而计算所得的频谱值自然也是有限、离散的。上述两种形式的傅里叶变换中，dfs变换满足时、频域自变量的离散化，但其时间变量和频率变量又同时具有周期性；dtft变换满足时间自变量的有限长度(非周期能量有限信号)，但其频率变量为连续形式。可见，这两种变换都难以实际应用。考虑到dfs变换的时、频域形式虽是周期序列，但每个周期却只有n个独立的复值，知道其一个周期的内容即可得到其它的内容。因此，若从dfs变换的时、频域各取出一个周期，即可构造出时间和频率自变量皆为离散、有限长度的傅氏变换，这就是离散傅里叶变换(dft)的引出思想，下面进行具体推导。设是一个长度为m的有限长序列，由周期序列与有限长序列的本质联系，可以n()为周期将展开为无重叠的周期序列，即周期延拓为 (1.3)再利用式(4.3.1)对进行dfs变换，得到周期离散的频谱，取的主值序列，代入dfs反变换公式(4.1.3b)，即 (1.4)分析可见：在dfs正变换中，只要把一个周期内的乘以对应的，即可得任意k下的；同理，在dfs反变换中，仅用的一个周期的值，即可得到任意n下的。如果同时限制(4.3.1)式中的n和(4.3.4)式中的k，使其都只在区间内取值，就得到了一个周期的和一个周期的间的对应关系 (1.5) (1.6)式中，n为dft变换区间的长度，上两式即称为有限长序列的离散傅里叶变换对。(4.3.5)式称为离散傅里叶变换，简称dft；(4.3.6)式称为离散傅里叶逆变换(inverse discrete fourier transform)，简称idft。由上述dft的推导过程和定义可以看出，dft变换具有隐含周期性。一方面由于dft变换由dfs变换引出，另一方面由dft定义式中的满足下述周期性 m为整数 (1.7)使得有限长序列和均具有周期性，即满足因此，周期序列可以看作长度为n的有限长序列的周期延拓序列，如式 (4.3.3)，又可表示为 (1.8)式中是矩形序列，是余数运算表达式，即若 ()，则，表示以n为周期的周期序列。而可看作是从到的主值区间，表示为x(n)nn.图1.1有限长序列及其周期延拓 (4.3.9)上述关系如图1.1所示。同样可将表示为 上述分析同样说明：有限长序列的离散傅里叶变换，正好是的周期延拓序列的离散傅里叶级数系数的主值序列。也可以理解为，离散傅里叶变换(dft)是对序列的傅里叶变换(dtft)的频域采样，运用频域采样方法同样可以推导出dft的变换公式。2. 用dft对信号进行谱分析所谓信号的谱分析，就是计算信号的傅里叶变换。工程实际中，经常遇到的是连续信号，截取一段进行ft变换后其频谱函数也是连续函数，因此其计算过程不便于用计算机实现。而dft作为一种时域和频域均离散化的变换，适合数值运算，成为分析离散信号和系统的有力工具。对于连续信号和系统，可以通过对时域采样，得到离散时间序列，再对进行dft，因此是一种近似的谱分析。下面讨论用dft对连续信号进行谱分析的原理和误差问题。 设连续信号的持续时间为，最高频率为，其傅里叶变换为，如图1.2(a)所示。根据采样定理，对以采样间隔进行采样得，设共采样n点，由采样定理可知，当满足时，采样信号的傅里叶变换是的周期延拓，如图1.2 (b)所示。由于仍然是的连续周期函数，为实现dft，需要对其进行频域离散化，即在区间上等间隔采样n点，采样间隔为f，如图1.2 (c)所示。00 1 2 00图 1.2 dft信号频谱分析(a)(b)(c) 可见，对连续信号频谱进行数字处理时，既会遇到时域采样，也要处理频域采样。时域采样和频域采样有关的几个参数为：数据时间长度、时域采样频率、时域采样间隔t、频域采样点数n、频域采样间隔f，它们之间的关系为： (1.9)上述分析表明，在满足时域采样定理的条件下，可以通过对连续信号进行采样并进行dft，来近似地反映连续信号的频谱特性。由于dft变换需要进行时域采样和频域采样，因此这种近似必然带来频谱分析的一定误差，下面从三方面分析其误差原因和解决措施：(1) 混叠效应从上面的讨论可知，连续信号频谱进行数字处理时，均为有限长序列。而傅里叶变换理论指出，一般时宽有限的信号，其频谱无限长，例如单个矩形脉冲信号的频谱，反之亦然。即从理论上说，没有有限时宽的限带信号。而由处理技术的可实现性，实际上只能处理有限时宽信号。因此对频带无限的时域信号采样后，在频域中会出现混叠，形成频谱失真，不能反映原信号的全部信息，这就是混叠效应。所以用dft作频谱分析是近似分析，当然，对不同的场合要求，可以有不同的逼近程度，从工程角度讲这是允许的。为了减小频谱混叠效应，可以采用预滤波方法滤除一定的高频成分，使待处理信号的有效带宽小于折叠频率。为了进一步减小混叠效应，除了采用预滤波技术，通常采样频率 在。在已知信号的最高频率(即谱分析范围)时，为了避免在dft运算中发生频率混叠现象，要求采样速率满足下式： (2) 栅栏效应我们知道，n点dft是频率区间上对信号的频谱进行n点等间隔采样，而采样点之间的频谱函数值是不知道的。这就好像从(n+1)个栅栏缝隙中观看信号的频谱情况，仅得到n个缝隙中看到的频谱函数值。因此称这种现象为栅栏效应。由于栅栏效应，有可能漏掉(挡住)大的频谱分量。为了把原来被“栅栏”挡住的频谱分量检测出来，需要提高频率分辨率(用频率采样间隔f来描述，表示谱分析中能够分辨的两个频率分量的最小间隔)，也就是说：f越小，谱分析越接近原连续频谱，频率分辨率越高。按照(4.3.30)式，为使f减小，需要减小采样频率或增大采样点数n。由于采样频率的降低会引起谱分析范围的减少或者引起频谱混叠，所以在维持不变的条件下，频率分辨率可以增加采样点数n，因为，只有增加对信号的观察时间，才能增加n。和n可以按照下面两式进行选择： (1.10) (1.11)增加观察时间或者采样点数最简单的办法可以采用在原序列尾部补零的方法，改变序列长度n(即改变dft变换区间长度)，从而增加了频域采样点数和采样点位置，使原先漏掉的某些频谱分量被检测出来。对连续信号的谱分析，只要采样速率足够高，且采样点数满足频率分辨率要求，就可以认为dft变换后所得离散谱的包络近似代表原信号的频谱。 (3)截断效应 在谱分析时，对于持续时间很长的信号或者无限长信号，会由于采样点数太多而无法存储和计算，因此只能截取有限点进行dft。截断后的有限长序列表示为其中，为矩形窗函数。根据傅里叶变换的频域卷积定理，有 其中图 1.3 、及二者的卷积和分别是原序列和矩形窗序列的频谱，截断后序列的频谱为和的频域卷积，若令，采样频率为，则其、及其二者的卷积如图 1.3所示。由上述可见，截断后序列的频谱与原序列频谱必然有差别，这种差别对谱分析的影响主要表现在如下两个方面：(a)泄漏。由图1.3可知，原来序列的频谱是离散谱线时，经截断后，使原来的离散谱线向附近展宽，通常称这种展宽为泄漏。显然，泄露使频谱变模糊，使谱分辨率降低。(b)谱间干扰。在主谱线两边形成很多旁瓣，引起不同频率分量间的干扰(简称谱间干扰)，影响频谱分辨率，特别是强信号谱的旁瓣可能湮没弱信号的主谱线，或者把强信号谱的旁瓣误认为是另一信号的谱线，从而造成假信号，这样就会使谱分析产生较大偏差。由于上述两种影响是由对序列截断引起的，所以称之为截断效应。减小截断效应可以采用改善窗函数形状的方法，即用其它形状的窗函数代替矩形窗。 第二章 离散傅里叶变换的快速算法(fft)dft在数字信号处理中有很重要的作用，如频谱分析、线性卷积等，但由于直接计算dft的计算量与变换区间长度n的平方成正比，当n较大时，计算量太大，所以在快速傅里叶变换(fast fourier transform，简称fft)出现前，直接用dft进行谱分析和信号的实时处理是不切实际的。直到1965年库利(j.w.cooley)和图基(j.w.tukey)提出了dft运算的一种快速方法以后，情况才发生了根本的变化。多年来，人们不断改进和完善，形成了一系列新型fft算法，如基-2fft算法、分裂基fft算法、n为复合数的fft算法等。必须强调指出：fft并不是与dft不同的另外一种变换，而是为减少dft计算次数的一种快速算法。为了解fft高效算法的重要及实现思路，先介绍dft的运算特点，再具体讨论高效算法。1. 直接计算dft的问题和改善dft运算效率的基本途径有限长序列的n点dft为 (2.1)逆变换为 (2.2)二者的差别仅在于：的指数符号不同且差一个比例因子。因此，下面我们仅就(4.3.29)式算法进行讨论。一般由于、都是复数，也是复数，因此直接按(4.3.33)式计算某个值需要n次形式的复数乘法及次复数加法的运算。共有n个点(k从0取到)，所以完成全部dft的总计算量则为次复数相乘及次复数相加。可见，直接计算dft时，乘法次数与加法次数都是和成比例的。这样，当n很大时，所需的运算工作量非常可观。例如，n10点的dft，需要100次复数相乘，而n1024时，则需要1 048 576即一百多万次的复乘运算。这对于实时性很强的信号处理(如雷达信号处理)来说，必将对计算速度有十分苛刻的要求。为此，迫切需要改进dft的计算方法，以减少总的运算次数。从哪些方面能改进dft的运算以减少运算工作量呢? 如前所述，dft的运算量是与成正比的。显然，如果一个大点数n的dft能分解为若干小点数dft的组合，则可达到减少运算工作量的效果。充分利用系数的下列固有周期性和对称性，使dft运算中的有些项可以合并，从而使长序列的dft分解为更小点数的dft，提高运算效率。的对称性 的周期性 快速傅里叶变换算法正是基于上述基本思想而发展起来的。下面介绍最常用的基2 fft算法(，库利一图基算法)。2. 基2 fft算法基2 fft算法主要包括两类：按时间抽取(decimation-in-time，简称dit)法和按频率抽取(decmation-in-frequence，简称dif)法，首先介绍dit算法。设是列长为的输入序列，且，其中为整数。如果不满足这个条件，可以人为地加入若干零点来达到。将按n的奇偶分成两个子序列 (2.3)则(2.1)式可化为由于，故上式又可表示为 (2.4)其中和分别是及的点的dft (2.5) (2.6)(2.6)式表明了个n点的dft被分解为两个点的dft。但是这里有一个问题，即，的列长为，它们的dft，的点数也是，即，而却有n个点，所以按(2.4)式计算得到的只是 ()的前一半项数的结果，要用来表达全部的值还必须应用w系数的周期性，即 这样可得同理可得蝶形运算符号另外又考虑到的对称性因此，将上述公式代入(2.4)中，又可表达为 (2.7) (2.8)由上分析可见，只要求出区间内各个整数k值所对应的和值，即可求出区间内的全部值，这一点恰恰是fft能大量节省计算的关键所在。式(2.7)、(2.8)的运算可用信号流图符号表示，根据其形状称之为蝶形运算符号。由上图，要完成一个蝶形运算，需一次复数乘法和二次加(减)法。据此，一个n点的dft分解为两个点的dft(设此点的dft是由直接计算得到的，下面将可看到，它们还可被继续分解)，则各需(次复乘和次复加，两个dft则需要次复乘和次复加。将两个点的dft合成为n点的dft时，需要再进行个蝶形运算，即还需要次乘法和次加法运算。因此通过这样分解后，计算全部共需要次复乘和次复加。前已指出，直接计算n点dft需要次复乘和次复加。由此可见，仅仅作了一次分解，即可使计算量差不多节省了一半。既然这样分解是有效的，由于，仍然是偶数，所以可以进一步把每个点子序列再按上述dit方法进一步分解。与第一次分解相同，将按奇偶分解成两个n/4长的子序列和，即 那么，又可表示为 同理 上面式中 也可进行同样分解 而 上面式中 图（a） n点dit-dft运算流图（n=8） 这样，经过第二次分解，又将n/2点dft分解为两个n/4点dft蝶形与运算，如图4.3.14所示。依此类推，经过m-1次分解，最后将n点dft分解成n/2个2点dft。一个完整的8点dit-dft运算流图如上图所示。图中用到关系式。图中输入序列不是顺序排列，后面会看到，其排列是有规律的。图4.3.14中8点dft可分解为四个点的dft，其运算过程为：先作点的dft计算，再令相应的两个点dft的结果，合成点的dft，从而得到、。最后按(4.3.39)式，(4.3.40)式组合成n点的dft。这种方法，由于每一步分解都是按每级输入序列在时间上的次序是属于偶数还是奇数来分解为两个更短的子序列，所以称为“按时间抽取法”。l-1级l级l级基本运算碟形流图3. 按时间抽取得fft算法的特点为了最终任何点的dit-fft运算程序或者设计出硬件电路，下面介绍它的运算规律。(1) 原位运算 由所述算法原理及上图的信号流图可见，的计算中共有m级，每级有n/2个蝶形，每级计算是由n个复数节点值经个蝶形运算变成了另外n个复数据。第l级基本运算蝶形如图4.3.17所示，其运算关系为由信号流图4.3.15可见，为了计算第l级第p、q位置上的复数节点值，只需要第l-1级的p、q位置上的复数节点值，而和其它节点变量无关。因此，如果所有的的值己预先置存，那么除了运算的工作单元之外，只用n个寄存器就够了。因为每个蝶算是由两个寄存器中取出数据，而计算结果仍存放到此二寄存器中，该寄存单元中原存的内容，一经取用即可抹去，不影响以后的计算。每列的个蝶算全作完以后再开始下一列的蝶算。这种计算通常称为“原位”(inplace)运算，因为整个计算实际只需要一组(n个)复数存储单元，所以这种结构可以节省存储单元，降低设备成本。(2) 蝶形运算规律 如上所述，n点dit-fft运算流图中，每级都有个蝶形，每个蝶形都要乘以因子，称为旋转因子，r为旋转因子的指数。由图可见，运算级数不同，旋转因子和循环方式就不同。要编写计算程序，就要首先找出旋转因子与运算级数的关系。用l表示从左到右的运算级数，由图不难发现，第l级共有个不同的旋转因子。例如，当时各级旋转因子表示如下：l=1 ， ，一个旋转因子开始 dit-fft的运算程序流程图倒 序输 出结束l=2 ，两个旋转因子、l=3 ，四个旋转因子、对的一般情况，由于，所以第l级的旋转因子为，式中。由上式即可确定第l级的旋转因子。在此基础上可以得到第l级蝶形的一般运算规律。设序列经时间抽取后，存入数组x中。如果蝶形运算的两个输入数据相距b点，应用同址计算，则第l级蝶形运算可表示为基于上面的分析，可以归纳出dit-fft的运算规律：第l级每个蝶形的两个输入数据相距点；同一旋转因子对应着间隔为的个蝶形。根据上述规律，从输入端(第一级)开始，逐级进行，共作m级运算。在作第l级运算时，依次求出个不同的旋转因子，每求出一个旋转因子，就计算完它对应的所有个蝶形。这样可以用三重循环实现fft运算。程序流图如左图所示。(3)输入序列的序号及倒序规律 由上图可见，当按原位进行计算时，fft输出端的次序正好是顺序排列的，即。但这时输入却不能按自然顺序存入存储单元，而是按的顺序存入存储单元，这是由时间抽取法不断在时序中将长点的dft(按奇、偶)分解为短点dft引起的，可用下式表示：4. 按频率抽取fft算法(dif-fft)与按时间抽取的基2fft算法相对应，还存在另一种称为按频率抽取的算法dif-fft。这个方法不是将序列分解，而是将代表频域的输出序列 按其顺序是属于偶数还是奇数分解为越来越短的序列。对长度的序列，先将按n的顺序分成前后两半(注意，这里是纯前后而不是奇偶)，由后面分析可见则是按奇，偶分组。的dft可表示为由式中 可将x(k)分解为偶数组和奇数组，即分别令及， 则 则l级基本运算蝶形流图 上式表明，x(k)按照奇偶k值分为两组新序列的点dft，其偶数组是的点dft，奇数组是的点dft，其运算关系可用图4.3.21所示的蝶形运算流图符号来表示。与时间抽取法的推演过程一样，由于，仍是一个偶数，因此可以将点的dft输出再分解为偶数组与奇数组。这样就将点的dft进一步分解为2个点的dft。这两个点dft的输入分别是将点dft的输入和上下对半分开，通过蝶形运算而形成，情况与第一步分解相同。这样的分解可一直进行下去，经过m-1次分解后变成了求个两点的dft，而两点dft就是一个基本蝶形运算单元。由下图可见，这种算法是对x(k)按照奇偶抽取分解的结果，所以称之为频域抽取法fft。比较可知dit与dif两种算法既有相似又有差别，相似之处是：均可原位计算，共有m级运算，每级共有n/2个蝶形运算，所以两种算法的运算次数相同；不同的是dif-fft算法输入为自然顺序，输出为倒序排列，这与dit的情况正好相反。因此dif的输出序列要经过“整序”才能变为自然顺序输出，整序的规律和时间抽取法相同。另外，dif的蝶形运算与dit的蝶形运算略有不同，其差别在于dit蝶形先乘旋转因子后加减，dif中先加减后乘旋转因子。 （b） dif-fft的运算流图(n=8)实际上，比较图（a）与图（b）不难看出二者互为转置，即把图（a）反一个面并倒转信号流图的方向和交换输入与输出，即可得图（b）。同理，也可通过转置，从图（b）得到图（a）。概括地说：对于每一种按时间抽取的fft算法都存在一种按频率抽取的算法，二者互为转置。第三章 用matlab实现dft和fft以及对两者运算时间的分析1.两种算法的时间计算根据对dft以及对fft的定义，用matlab写出能够实现dft和fft的方程，分别保存为dftmatlab.m和fftmatlab.m文件，放到matlab的work文件夹内。打开matlab，输入以下代码，生成一个含有100元素的数组a：n=1:100;a=sin(n)+cos(n);分别用dft和fft对数组a进行傅里叶变换：a1=dftmatlab（a，100）；a2=fftmatlab（a）；原图fft变换后dft变换后比较观察两组数据,利用matlab里面的记录当前时间按函数clock和计算两过程时间差函数etime分别计算出对数组a进行dft变换和fft变换所用的时间：比较两图，发现当对函数进行fft时所用时间比dft所用的时间缩短很多。2.分析两者运算时间的差异：由前面介绍的按时间抽取的fft运算流图可见，每一级都由个蝶形运算构成。因此每一级运算都需要次复乘和n次复加(每个结加减各一次)。这样m级运算总共需要复乘数 (4.3.45)复加数 (4.3.46)实际计算量和这个数字稍有出入，因为由运算流图可见，这种情况共有次，这几个系数是都不用乘法运算的，但这种情况在直接计算dft中也是有的，且当n较大时，这些影响也较小。所以为了统一作比较我们不考虑以上特例。综上所述，可以得出如下结论；按时间抽取法所需的复乘数和复加数都是与成正比的，而直接计算dft时所需的复乘数为，复数加为n(n-1)次。当n1时，从而，dit-fft算法比直接计算dft的运算次数大大减小。例如时，0641282565121024641282565121024直接计算fft算法乘法次数(1000)n(取样点数)图 4.3.16 fft算法与直接计算dft所需乘法次数的比较曲线这样，就使运算效率提高200多倍。图4.3.16为fft算法和直接dft算法所需运算量与计算点数n的关系曲线。由此图更加直观地看出fft算法的优越性，显然，n越大时，优越性就越明显。随着运算次数的减少，fft的运算时间明显减少。第四章 总结过对这次为期两周的数字信号处理课程设计，使我对数字信号处理有了更深一步的了解，尤其是在fft和dft这方面的认识与了解。同时我也知道了，任何一个算法的出现或软件程序的出现，都饱含着相关人员的汗水。在老师的耐心指导以及通过去图书馆查阅资料，拓宽了我的知识面，学到了许多在书本上没有学到的东西。通过fft与dft两种算法的对比，我知道了对于任何事物，只要肯钻研，再掌握了基础知识以后就能够活学活用，就像fft算法一样，就是再充分掌握dft的基