毕业论文-快速傅里叶变换算法及其在信号处理中的应用(共45页)

1、武汉工程大学毕业设计(论文)说明书PAGE 2015 届毕业设计(b y sh j)(论文) 题 目快速傅里叶变换算法及其在信号处理中的应用专 业 班 级2011电子信息工程02 学 号1104030231 姓 名周汝耀 指 导 教 师华夏讲师 学 院 名 称电气信息学院 2011 年 6 月 9 日摘 要快速傅氏变换（FFT），是离散傅氏变换的快速算法，它是根据离散傅氏变换的奇、偶、虚、实等特性，对离散傅立叶变换的算法进行(jnxng)改进获得的。它对傅氏变换的理论并没有新的发现，但是对于在计算机系统(xtng)或者说数字系统中应用离散傅立叶变换，可以说是进了一大步。傅里叶变换的理论(lln

2、)与方法在“数理方程(fngchng)”、“线性系统分析”、“信号处理、仿真”等很多学科领域都有着广泛应用,由于计算机只能处理有限长度的离散的序列,所以真正在计算机上运算的是一种离散傅里叶变换.虽然傅里叶运算在各方面计算中有着重要的作用，但是它的计算过于复杂，大量的计算对于系统的运算负担过于庞大，使得一些对于耗电量少，运算速度慢的系统对其敬而远之，然而，快速傅里叶变换的产生，使得傅里叶变换大为简化，在不牺牲耗电量的条件下提高了系统的运算速度，增强了系统的综合能力，提高了运算速度，因此快速傅里叶变换在生产和生活中都有着非常重要的作用，对于学习掌握都有着非常大的意义。 关键词：快速傅氏变换；快速算

3、法；简化；广泛应用 AbstractFast Fourier Transform (FFT), is a discrete fast Fourier transform algorithm, which is based on the Discrete Fourier Transform of odd and even, false, false, and other characteristics of the Discrete Fourier Transform algorithms improvements obtained. Its Fourier transform theory h

4、as not found a new, but in the computer system or the application of digital systems Discrete Fourier Transform can be said to be a big step into. Fourier transform theory and methods in the mathematical equation and linear systems analysis and signal processing, simulation, and many other areas hav

5、e a wide range of applications, as the computer can only handle a limited length of the sequence of discrete, so true On the computers operation is a discrete Fourier transform. Fourier Although all aspects of computing in the calculation has an important role, but its calculation was too complicate

6、d, a lot of computing system for calculating the burden is too large for some Less power consumption, the slow speed of operation of its system at arms length, however, have the fast Fourier transform, Fourier transform greatly simplifying the making, not in power at the expense of the conditions to

7、 increase the speed of computing systems, and enhance the system The comprehensive ability to improve the speed of operation, the Fast Fourier Transform in the production and life have a very important role in learning to master all have great significance.KeyWords：Fast Fourier Transform; fast algor

8、ithm; simplified; widely used 目 录摘要(zhiyo).IAbstract . II TOC o 1-3 h z u HYPERLINK l _Toc385634561 1.绪论(xln) HYPERLINK l _Toc385634562 1.1 选题(xun t)背景1 HYPERLINK l _Toc385634563 1.2 课题研究的意义2 HYPERLINK l _Toc385634567 HYPERLINK l _Toc31137 2.快速傅里叶变换原理及性质 HYPERLINK l _Toc385634568 2.1快速傅里叶变换原理3 HYPER

9、LINK l _Toc385634569 2.2快速傅里叶变换的优越性4 HYPERLINK l _Toc385634570 HYPERLINK l _Toc9754 2.3快速傅里叶变换的意义4 HYPERLINK l _Toc385634573 HYPERLINK l _Toc27753 3.快速傅里叶变换的算法 HYPERLINK l _Toc385634574 HYPERLINK l _Toc8124 3.1快速傅里叶变换算法6 HYPERLINK l _Toc385634575 3.2 Cooley=Tukey FFT算法8 HYPERLINK l _Toc385634576 3.3

10、 Rader-Brennr FFT算法9 HYPERLINK l _Toc385634576 3.4 Goertsel 算法10 HYPERLINK l _Toc385634577 HYPERLINK l _Toc27753 4.快速傅里叶变换在信号处理中的理论应用 HYPERLINK l _Toc385634578 HYPERLINK l _Toc8124 4.1利用FFT计算连续时间信号的傅里叶变换13 HYPERLINK l _Toc385634579 4.2利用FFT计算离散信号的线性卷积17 HYPERLINK l _Toc385634580 4.3利用FFT进行离散信号压缩19 H

11、YPERLINK l _Toc385634579 4.4利用FFT对离散信号进行滤波22 HYPERLINK l _Toc385634579 4.5利用FFT提取离散信号中的最强正弦分量24 HYPERLINK l _Toc385634581 HYPERLINK l _Toc31137 5.快速傅里叶变换在数字信号分析与处理的实际应用 HYPERLINK l _Toc385634582 5.1快速傅里叶变换在喇曼光谱信号噪声平滑中的应用29 HYPERLINK l _Toc385634583 5.2采用异步实现的快速傅里叶变换处理器31 HYPERLINK l _Toc385634584 HY

12、PERLINK l _Toc9754 5.3快速傅里叶算法在哈特曼夏克传感器波前重构算法中的应用33 HYPERLINK l _Toc385634587 致谢36 HYPERLINK l _Toc385634586 参考文献37 1 绪论(xln)傅立叶变换在生产生活中的重要性非常(fichng)突出，它将原本难以(nny)处理的时域信号相比比较容易地转换成易于分析的频域信号，我们可以利用一些专业工具对这些频域信号进行处理、加工，使信号转化为可以对其进行各种数学变换的数学公式，对其进行处理。最后还可以根据傅立叶反变换将这些频域信号转换成原来的时域信号，这是一种特殊的积分变换。它能够将满足一定条

13、件的某个函数表示成为正弦基函数的线性组合或者积分。然而，它在运算上过于复杂，过于宏大的运算过程，对于一些相对简单的低功耗处理器来说，难以自如应对，因此，快速傅里叶变换则显出了它的优越性。快速傅氏变换（FFT），即离散傅氏变换的快速算法，它是根据离散傅氏变换的奇、偶、虚、实等特性，对离散傅立叶变换的算法进行改进获得的1。对于计算机处理信号方面上是一大进步。系统的速度不但取决于其本身的速度，而且还在相当大的程度上取决于运用的算法，算法运算量的大小直接影响到对设备的控制质量。通过傅立叶变换（DFT），运用测试软件进行检测，我们可以看出，快速傅里叶变换大大的提高了运算速度，它为各系统的设计方案提供了简

14、单算法，有着非常重要的意义。1. 1 选题背景近十多年来，数字信号处理技术同大规模集成电路、数字计算机等，都有了突飞猛进的发展，日新月异，早已成为了一门具有强大生命力的技术科学。因为它本身就具有一系列的优点，所以能够有效地促进工程技术领域的技术改造和科学发展，应用领域也更加地广泛、深入，越来越受到人们的重视。在信号处理中，离散傅里叶变换（Discrete Fourier Transform，DFT）是比较常用的变换方法之一，它在各种数字信号处理系统中扮演着及其重要的角色。由于离散傅里叶变换（DFT）而发现了频率离散化，可以直接用它来分析信号的频谱、计算滤波器的频率响应、以及实现信号通过线系统的

15、卷积运算等，因而在信号的谱分析等方面有着非常大的作用。傅里叶变换已经有一百多年的历史了，我们熟知频域分析往往比时域分析更优越，不仅简单明了，而且易于分析较为复杂的信号。但需要用较为精准的数字方法，即DFT，进行谱分析，在快速傅氏变换（FFT）出现以前是不切实际的。由于DFT的计算量太大，即使运用计算机也很难对问题进行实时的有效处理，所以DFT并没有得到真正的应用。直到1965年库利（J.W.Cooly）和图基（J.W.Tukey）首次发现DFT的一种快速算法，局面才发生根本性的变化。继库利和图基算法出来之后，桑德（G.Sander）等快速算法也相继出现，又经过其他学者一步步改进，很快就出现了通

16、用型的快速傅里叶变换，简称FFT。快速傅里叶变换（Fast Fourier Transform，FFT）并非与离散傅里叶变换完全不同的另一种变换，而是为了减少DFT计算次数而诞生的一种快速、有效的算法。应当指出的是，也是因为当时电子数字计算机的“落后”条件也促成了这个算法的提出。它使得DFT的运算量大大的缩小简化，它推动了近30年来信号处理技术止步不前的前进发展，成为了数字信号处理应用领域里强有力的工具，为DFT乃至数字信号处理技术的实际应用创造了良好的条件，从而使DFT在实际使用中得以广泛的应用2。数字(shz)信号处理器（DSP），是一种(y zhn)可编程的高性能处理器。近年来发展(fz

17、hn)尤为迅速，它不仅应用于数字信号处理方面，而且在图像处理、语音处理、通信等领域得到广泛的应用。之前通用的微处理器在运算速度上已经很难适应信号实时处理的高要求。DSP处理器中集成了高速的乘法硬件，能快速、准确地进行大量数据的乘法以及加法的运算。数字信号处理区别于普通的科学计算与分析，它强调运算的实时性。除了需要普通微处理器所强调的高速运算和控制能力之外，鉴于实时数字信号处理的特点，在处理器结构、指令系统、指令流程上做了很大程度上的改进。1. 2 课题研究的意义如上所述，基于对DSP的快速傅里叶变换算法的研究，从而使FFT算法能够有效地在DSP芯片上实现。DSP芯片的出现，使FFT的实现更加方

18、便。多数的DSP芯片都能够在一个指令周期内完成一次乘法和加法，并且提供了专门的FFT指令，完成一次指令的周期只需10ns，使得FFT算法在DSP芯片上实现的速度更加快速。快速傅里叶变换为频谱分析、卷积、相关数字滤波器设计与实现与功率谱计算、传递函数建模、图像处理等，提供了快速有效的运算方法。FFT技术应用DSP芯片，从而可以提供使调制、解调、压缩、解压缩的数据传输更为高效的信号处理解决方案，因而广泛应用于雷达、图像处理、通信、生物医学和声纳领域。2.快速傅里叶变换(binhun)原理及性质(xngzh) 数字信号中的傅里叶变换,通常是采用离散型傅里叶变换(DFT)。DFT 存在(cnzi)的缺

19、点就是计算量太大,不易进行实时处理。比如，计算一个N 点的DFT ,一般需要次复数乘法和N(N-1)次复数加法运算.因此,当N较大或要求对信号进行实时处理时,往往很难实现达到所需的运算速度。1965年,J.W.Cooly和J.W.Tukey发现了DFT的一种快速算法,经过后来学者的进一步改进, 很快便形成了一套高效的运算方法,即现在通用的快速傅里叶变换, 简称FFT( The Fast Fourier Transform)。快速傅里叶变换的实质是利用式(3-1)中的权函数的对称性和周期性,把N点DFT进行一系列分解和组合,使整个DFT的计算过程变成一系列叠代运算过程,使DFT的运算量大大简化,

20、为DFT及数字信号的实时处理和应用创造了非常良好的条件3。2. 1 快速傅里叶变换原理快速傅里叶变换原理：1. 将长序列DFT分解为短序列的DFT2. 利用旋转因子的对称性、周期性、可约性。将时域序列逐次分解为一组子序列，依据旋转因子的特性，由子序列的DFT来实现整个序列的DFT4。其中：快速傅里叶变换分为两种，分为基2时间抽取算法和基2频率抽取算法基2-时间抽取(Decimation in time)FFT算法其中:r=0,1,2 （2-1）基2-频率抽取(Decimation in frequency)FFT算法 （2-2）2. 2 快速(kui s)傅里叶变换的优越性设为项的复数(fsh

21、)序列，依据(yj)DFT变换，任一的计算都需要有次复数乘法和（）次复数加法，而且一次复数乘法等同于四次实数乘法和两次实数加法，同样的，一次复数加法等同于两次实数加法，即使我们把一次复数乘法和一次复数加法定义成一次“运算”（四次实数乘法和四次实数加法），那么求出项复数序列的, 即N点DFT变换大约就需要次运算。当点甚至更多的时候，需要N2=1048576次运算。在FFT中，利用WN的对称性和周期性，把一个N项序列（设,为正整数），分为两个项的子序列，而且每个点的DFT变换需要次运算，再运用N次运算把两个点的DFT 变换重新组合成一个N点的DFT变换。如此变换以后，总的运算次数就变成了。承接上面

22、的例子，当时，总的运算次数就变成了525312次，这样看来，节省了大约50%的运算量。而如果我们将这种“一分为二”的思想不断进行下去，直到分成两两一组的DFT运算单元，那么N点的 DFT变换就只需要次的运算，在1024点时，运算量仅有10240次，是先前的直接算法的1%，点数越多，运算量的节约就越大，这就是 FFT的优越性.当然，FFT提高了运算速度,但是,也对参与运算的样本序列作出了限制,即要求样本数为2N点.离散傅里叶变换DFT则无上述限制5。2. 3 快速傅里叶变换的意义 傅立叶变换的物理意义：傅立叶变换是数字信号处理技术领域一项很重要的算法。要知道傅立叶变换算法的意义，首先要了解傅立叶

23、变换原理的意义。傅立叶变换原理表明：任何连续测量的时序或信号，都能够表示成为不同频率的正弦波信号的无限叠加。而利用该原理而创立的傅立叶变换算法则利用直接能测量到的原始信号，并以累加方式来计算该信号中不同正弦波信号的频率、相位和振幅。与傅立叶变换算法(sun f)对应的是反傅立叶变换算法。该反变换从本质上说也就是一种累加处理(chl)，这样便可以将单独改变的正弦波信号(xnho)转换成一个信号。因此，也可以说，傅立叶变换将原来难以处理的时域信号转换成了易于分析的频域信号（信号的频谱），我们可以利用一些专业工具对这些频域信号进行加工、处理。最后还可以利用傅立叶反变换将这些频域信号转换成原来的时域信

24、号。从现代数学的眼光来看，傅里叶变换其实就是一种特殊的积分变换。它能够将满足一定条件下的某个函数表示成为正弦基函数的线性组合或者积分形式。在不同的研究领域里，傅里叶变换具有多种形式各异的变体形式，如连续傅里叶变换和离散傅里叶变换。在数学领域，尽管最初傅立叶分析是作为热过程的解析分析的工具，但是其思想方法仍然具有典型的还原论和分析主义的特征。任意的函数通过一定的分解，都能够表示为正弦函数的线性组合的形式，而正弦函数在物理上是被充分研究而相对简单的函数类：1. 傅立叶变换是线性算子,若赋予适当的范数,它还是酉算子;2. 傅立叶变换的逆变换容易求出,而且形式与正变换非常类似;3. 正弦基函数是微分运

25、算的本征函数,从而使得线性微分方程的求解可以转化为常系数的代数方程的求解.在线性时不变杂的卷积运算为简单的乘积运算,从而提供了计算卷积的一种简单手段;4. 离散形式的傅立叶的物理系统内,频率是个不变的性质,从而系统对于复杂激励的响应可以通过组合其对不同频率正弦信号的响应来获取;5. 著名的卷积定理指出:傅立叶变换可以化复变换可以利用数字计算机快速的算出(其算法称为快速傅立叶变换算法(FFT)。正是由于上述的良好性质,傅里叶变换在物理学、数论、组合数学、信号处理、概率、统计、密码学、声学、光学等领域都有着广泛的应用6。快速(kui s)傅里叶变换的算法3.1 快速(kui s)傅里叶变换算法 快

26、速(kui s)傅里叶变换FFT 是离散傅里叶变换DFT 的一种快速算法，只有FFT 才能在现实中有实际应用的意义。 QUOTE (3-1)由(3-1)式可知，对每一个n，计算X()须作N次复数乘法及N-1次复数加法，要完成这组变换共需 QUOTE 次乘法及N(N-1)次复数加法。但以下介绍的快速傅里叶变换的算法，可大大减少运算次数，提高工作效率。当时，n和k可用二进制数表示： QUOTE (3-2) （3-3）又记,则（3-1）式可改写为 （3-4）式中： (3-5)因为(yn wi)所以(suy)（3-2）可改成(i chn) （3-6） (3-7)则式（3-7）即为式（3-6）的分解形式

27、。将初始数据代入式（3-7）的第一个等式，可得每一组计算数据，一般将痗L-1组计算数据代入式（3-7）的第L个等式，计算后可得第L组计算数据（L，），计算公式也可表示为 = （3-8）式中 （3-9） 根据式（3-8），第L个数组中每个 的计算只依赖于上一个数组的两个数据这两个数据的标号相差，即，而且这两个数据只用于计算第L个数组中标号的数据（等号右端为二进制数）。当分别取和时，分别有。因此，用上一组的两个数据计算所得的两个新数据仍可储存在原来位置，计算过程中只需要N个存储器。将与称为第L个数组中的对偶结点对。计算每个对偶结点对只需一次乘法，事实上由式（3-8）可得 （3-10）式中： ；别为

28、式（3-9）中取，时对应(duyng)的P值。因，于是对偶(du u)结点的有如下(rxi)关系： ,因此式（3-8）可表示为 （3-11） P的求法：在中，i写成二进制数右移位，就成为颠倒位序得式（3-7）中，前面的r个等式，每个等式均对应一组数据进行计算，每组数据都有N/对结点，根据式（3-11），每对结点只需作次乘法和次加法，因此，每组数据只需N/2次乘法和N次加法，因而完成r组数据的计算共需Nr/2次乘法和Nr次加法。3. 2 Cooley-Tukey FFT算法的核心是将一层运算映射为二层运算。作点时，若不是素数，则可分解为，那么由的 （3-12）通过映射： （3-13）可得到(d

29、do) 而，可化简为 （3-14）从而(cng r)式（3-12）转化(zhunhu)为 （3-15）其中。以点为例，映射方式为：，则计算流图如图3-1所示。 n1 k2f0 0 W0 0 F0f4 1 W0 1 F5f8 2 W0 2 F10f12 3 W0 3 F15f16 4 W0 0 F1f1 0 W1 1 F6f5 1 W2 2 F11f9 2 W3 3 F16f13 3 f17 4 W0 0 F2 W2 1 F7f2 0 W4 2 F12f6 1 W6 3 F17f10 2f14 3 W0 0 F3f18 4 W3 1 F8 W6 2 F13f3 0 W9 3 F18f7 1f11

30、 2 W0 0 F4f15 3 W4 1 F9f19 4 W8 2 F14 W12 3 F19k1=0n2=0n2=1n2=2n2=3k1=1k1=2k1=3k1=4图3-1 Cooley-Tukey 20 点FFT算法(sun f)3. 3 Rader-Brenner FFT算法(sun f)Rader-Brenner 算法是类似于 Cooley-Tukey 算法，但是采用的旋转因子都是纯虚数，以增加加法运算和降低了数值稳定性为代价减少了乘法运算。这方法之后被split-radix variant of Cooley-Tukey所取代，与Rader-Brenner算法相比，有一样(yyng)

31、多的乘法量，却有较少的加法量，且不牺牲数值的准确性7。 Bruun以及QFT算法是不断的把DFT分成许多较小的DFT运算。(Rader-Brenner以及QFT算法是为了2的指数所设计的算法，但依然可以适用在可分解的整数上。Bruun算法则可以运用在可被分成偶数个运算的数字)。尤其是Bruun算法，把FFT看成是 ，并把它分解成 与 的形式。 另一个从多项式观点的快速傅立叶变换(binhun)法是Winograd 算法。此算法把 分 解成cyclotomic多项式，而这些多项式的系数通常为1，0，-1。 这样只需要很少的乘法量(如果有需要的话)，所以(suy)winograd是可以得到最少乘法

32、量的快速傅立叶算法，对于较小的数字，可以找出有效率的算方式。更精确地说，winograd算法让DFT可以用2K点的DFT来简化，但减少乘法量的同时(tngsh)，也增加了非常多的加法量。Winograd也可以利用剩余值定理来简化DFT。 Rader算法提出了利用点数为N(N为质数)的DFT进行长度为N-1的回旋折积来表示原本的DFT，如此就可利用折积用一对基本的FFT来计算 DFT。另一个prime-size的FFT算法为chirp-Z算法。此法也是将DFT用折积来表示，此法与Rader算法相比，能运用在更一般的转换 上，其转换的基础为Z转换8。3.4 Goertsel 算法如前所述, 点时域

33、序列的离散付里叶变换式为， （3-16）这点频域序列是同时被算出的，不可能只计算其中某一个或几个指定点。Goetzel 算法是为了解决这个问题而提出的。这个算法把离散付里叶变换看作一组滤波器，将输入端的时域序列与其中一个滤波器的冲激响应序列进行卷积运算，求滤波器的输出序列，即得序列的一点。这种算法利用旋转因子的周期性，使DFT运算化为线性滤波运算9。 由于 故式（3-16）可化为 （3-17）定义序列为 （3-18）可见(kjin)是由两个(lin )序列卷积而得到的序列。 （3-19）其中(qzhng)，是输入的点序列，另一个序列被看作滤波器的冲激响应序列。 （3-20）对比式(3-17）和

34、式(3-18），可知：按式（3-19）进行卷积运算，当时，滤波器的输出就是： （3-21）对式(3-20)进行Z变换，可得滤波器的系统函数 （3-22）这是一个一阶系统。图3-2示出这个系统的信号流图，相应的差分方程为 ， (3-23) 按照此式进行递推运算，到了 时刻，即可依据式（3-21）得到。 按照式(3-22)进行运算时，可先算好旋转因子，储存起来。每次递推包含一次复数乘法。按式(3-16)直接计算点离散付里叶变换,需要次实数乘法和次实数加法。按照上述Goertzel算法，所需的实数乘法和实数加法都是次。所以当不大时，上述算法的效率稍差10。下面介绍改进的Goertzel 算法，这种算

35、法所需的实数乘法次数约为直接方法的一半。图3-2 用一阶系统实现(shxin)Goertzel 算法 图 3-3 用二阶系统(xtng)实现(shxin)Goert算法把式（3-22）的分子和分母都乘上因子，就得到第个滤波器的系统函数为 （3-24）与此相应的信号流图示于图3-3。由式（3-24)可见，滤波器是一个二阶系统，有一对复数共轭极点和一个复数零点。为了便于运算，在图3-3所示的流图中，设立状态变量 和。按照图3-3计算时，步骤有二，即： 1.实现一对复数极点输入点依次取 ，进行递推运算。每次运算中，更新状态变量和。作次迭代所需的计算量是次实数乘法和次实数加法。 2.实现复数零点。是一

36、个点序列， 。在点上。计算状态变量和。这时，按照图3-2算出滤波器的输出，此即。所需的计算量是4次实数乘法和次实数加法11。 综上所述，计算一点需要进行次实数乘法和次实数加法。这种算法要求的乘法次数约为直接算法的一半。在这种较为有效的方案中，仍具有这样的优点，即必须计算和存储的系数只有和。还要说明(shumng)图3-3所示的算法的另一个优点。当输入(shr)序列为实序列时，离散付里叶变换序列是对称(duchn)的，即。容易证明，图3-3的网络形式在计算时和计算时具有完全相同的极点，但前者的零点系数与后者的零点系数成复共轭关系。由于零点仅在最后的迭代中实现，所以诸极点要求的次乘法和 次加法可以

37、用来计算离散付里叶变换的两个值。因此，若用Goertzel 算法计算离散付里叶变换的所有个点的值，需要的乘法次数近似为，加法次数近似为。然而，它同直接计算离散付里叶变换一样，计算量仍然正比于。4.快速傅里叶变换在信号处理中的理论应用4.1 连续时间信号的快速傅里叶变换设是连续时间信号，并假设时，则其傅里叶变换由下式给出： （4-1） 令是一固定的正实数，是一个固定的正整数。当时，利用FFT算法可计算。已知一个固定的时间间隔，选择让T足够小，使得每一个T秒的间隔内，的变化很小，则式中积分可近似为 （4-2）假设(jish)足够(zgu)大，对于所有的整数(zhngsh)，幅值很小，则式（4-2）

38、变为 （4-3）当时，式（4-2）两边的值为 （4-4）其中代表抽样信号的点。最后令，则上式变为 （4-5）首先用FFT算法求出，然后可用上式求出时的。应该强调的是，式（4-3）只是一个近似表示，计算得到的只是一个近似值。通过取更小的抽样间隔，或者增加点数，可以得到更精确的值。如果时，幅度谱很小，对应于奈奎斯特抽样频率，抽样间隔选择比较合适。如果已知信号只在时间区间内存在，可以通过对时的抽样信号补零，使足够大12。将连续时间傅立叶变换进行数字(shz)近似，用函数fft（快速傅立叶算法）高效地计算这个近似值。很多信号都能用（4-1）式连续时间傅立叶变换（CTFT）来表示。利用MATLAB可以计

39、算（CTFT）积分的数值近似。利用在密集(mj)的等间隔t的样本上的求和来近似这个积分，就可以用函数fft高效地计算这个近似值。所用的近似式是根据积分的定义得到的，即 （4-6）对于一般(ybn)信号，在足够小的下，上式右边的和式是对于CTFT积分的一个好的近似。若信号对于和为零，那么这个近似式就能写成 （4-7）式中，N为一整数。可以利用函数fft对一组离散的频率计算上式中的和式。如果N个样本是存在向量x内的话，那么调用函数X=tau*fft(x)就可以计算出 （4-8）式中 以及N假设为偶数。为了计算高效，fft在负的频率样本之前先产生正频率样本。为了将频率样本置于上升的顺序，能用函数ff

40、tshift。为了将存入X中的的样本排列成使就是对于，在上求得的CTFT,可用X=fftshift(tau*fft(x)。例4-1 利用(lyng)FFT计算傅里叶变换如图4-1所示的信号(xnho)其傅里叶变换(binhun)为：利用下面的命令，可得到的近似值和准确值。 图4-1 连续时间信号x(t)N=input(Input N:);T=input(Input T:);%计算X(w)近似值t=0:T:2;x=t-1 zeros(1,N-length(t);X=fft(x);gamma=2*pi/(N*T);k=0:10/gamma;Xapp=(1-exp(-i*k*gamma*T)/(i*

41、k*gamma)*X;%计算真实值X(w)w=0.05:0.05:10;Xact=exp(-i*w)*2*i.*(w.*cos(w)-sin(w)./(w.*w);plot(k*gamma,abs(Xapp(1:length(k),o,w,abs(Xact);legend(近似值,真实值);xlabel(频率(rad/s);ylabel(|X|)运行程序后输入N=128，T=0.1，此时，得到实际的和近似的傅里叶变换的幅度谱如图4-2所示，此时近似值已经相当准确。通过增加NT可以增加更多的细节，减少T使得到的值更精确。再次运行程序后输入N=512，T=0.05，此时，得到实际的和近似的傅里叶变

42、换的幅度谱如图4-3所示。图4-2 N=128，T=0.1时的幅度(fd)谱 图4-3 N=512，T=0.05时的幅度(fd)谱4.2 FFT计算(j sun)离散信号的线性卷积已知两个离散时间(shjin)信号与，取，对和右端补零，使得(sh de) （4-9）利用FFT算法可以求得和的L点DFT，分别是和，利用DTFT卷积性质，卷积等于乘积的L点DFT反变换，这也可以通过FFT算法得到。例4-2 利用FFT计算线性卷积已知，其中为单位阶跃序列，信号如图4-4所示。由于当时，很小，故可以取为17；N取10，。利用下面的Matlab命令，可得到、的卷积图形如图4-4所示。subplot(3,

43、1,1);n=0:16;x=0.8.n;stem(n,x);xlabel(n);ylabel(xn);subplot(3,1,2);n=0:15;y=ones(1,10) zeros(1,6); stem(n,y);xlabel(n);ylabel(yn)subplot(3,1,3);L=26;n=0:L-1;X=fft(x,L);Y=fft(y,L);Z=X.*Y;z=ifft(Z,L);stem(n,z);xlabel(n);ylabel(zn)图4-4 信号(xnho)xn、yn及其卷积zn=xn*yn利用(lyng)下面的Matlab命令，可得到信号xn、yn的幅度谱与相位谱如图4-5

44、所示。subplot(2,2,1);L=26;k=0:L-1;n=0:16;x=0.8.n;X=fft(x,L);stem(k,abs(X);axis(0 25 0 5);xlabel(k);ylabel(|Xk|)subplot(2,2,2);stem(k,angle(X);axis(0 25 -1 1);xlabel(k);ylabel(Angle(Xk)(弧度(hd)subplot(2,2,3);y=ones(1,10);Y=fft(y,L);stem(k,abs(Y);axis(0 25 0 10);xlabel(k);ylabel(|Yk|)subplot(2,2,4);stem(k

45、,angle(Y);axis(0 25 -3 3);xlabel(k);ylabel(Angle(Yk)(弧度(hd)图4-5 信号(xnho)xn、yn的幅度谱与相位谱4.3 FFT进行离散信号(xnho)压缩利用FFT算法对离散信号进行压缩的步骤如下：1）通过采样将信号离散化；2）对离散化信号进行傅里叶变换；3）对变换后的系数进行处理，将绝对值小于某一阈值的系数置为0，保留剩余的系数；4）利用IFFT算法对处理后的信号进行逆傅里叶变换13。 例4-3 对单位区间上的下列连续信号以采样频率进行采样，将其离散化为个采样值用FFT分解信号，对信号进行小波压缩，然后重构信号。令绝对值最小的80%系

46、数为0，得到(d do)重构信号图形如图4-6a)所示，均方差为0.0429，相对误差为0.0449；令绝对值最小的90%系数(xsh)为0，得到重构信号图形如图4-6b)所示，均方差(fn ch)为0.0610，相对误差为0.0638。 a) 绝对值最小的80%系数为0的重构信号（FFT） b) 绝对值最小的90%系数为0的重构信号（FFT）图4-6 用FFT压缩后的重构信号相关Matlab程序如下:function wc=compress(w,r)%压缩函数compress.m%输入信号数据w,压缩率r%输出压缩后的信号数据if(r1) error(r 应该介于0和1之间!);end;N=

47、length(w);Nr=floor(N*r);ww=sort(abs(w);tol=abs(ww(Nr+1);wc=(abs(w)=tol).*w; function unbiased_variance,error=fftcomp(t,y,r)%利用FFT做离散信号压缩%输入时间t,原信号y,以及压缩率rif(r1) error(r 应该(ynggi)介于0和1之间!);end;fy=fft(y);fyc=compress(fy,r); %调用压缩(y su)函数compress.myc=ifft(fyc);plot(t,y,r,t,yc,b);legend(原信号(xnho),重构信号);

48、unbiased_variance=norm(y-yc)/sqrt(length(t);error=norm(y-yc)/norm(y);输入以下Matlab命令：t=(0:255)/256;f=t+cos(4*pi*t)+1/2*sin(8*pi*t);unbiased_variance,error=fftcomp(t,f,0.8) unbiased_variance = 0.0429error =0.0449如果用Harr尺度函数和Harr小波分解信号，对信号进行小波压缩，然后重构信号。令绝对值最小的80%系数为0，得到重构信号图形如图4-7a)所示，均方差为0.0584，相对误差为0.0

49、611；令绝对值最小的90%系数为0，得到重构信号图形如图4-7b)所示，均方差为0.1136，相对误差为0.1190。 a) 绝对值最小的80%系数为0的重构信号（Harr） b) 绝对值最小的90%系数为0的重构信号（Harr）图4-7 用Harr小波压缩(y su)后的重构信号相关Matlab程序(chngx)如下function unbiased_variance,error=daubcomp(t,y,n,r)%利用Daubechies系列小波做离散信号(xnho)压缩%输入时间t,原信号y,分解层数n,以及压缩率r%输出原信号和压缩后重构信号的图像,以及重构均方差和相对l2误差if(

50、r1) error(r应该介于0和1之间!);end;c,l=wavedec(y,n,db1);cc=compress(c,r); %调用压缩函数compress.myc=waverec(cc,l,db1);plot(t,y,r,t,yc,b);legend(原信号,重构信号);unbiased_variance=norm(y-yc)/sqrt(length(t);error=norm(y-yc)/norm(y);输入以下Matlab命令：=(0:255)/256;f=t+cos(4*pi*t)+1/2*sin(8*pi*t);unbiased_variance,error=daubcomp(

51、t,f,8,0.8) unbiased_variance = 0.0584error = 0.0611结论：在信号没有突变、快变化或者大致上具有周期性的信号，用FFT可以处理得很好（甚至比小波还要好）。4. 4 利用FFT对离散信号进行滤波利用FFT算法对信号进行滤波的步骤如下：1）通过采样将信号离散化；2）对离散化信号进行傅里叶变换；3）对变换后的系数进行处理，将绝对值大于某一阈值的系数置为0，保留剩余的系数；4）利用IFFT算法对处理后的信号进行逆傅里叶变换14。例4-4 对被白噪声污染(zo shn w rn)的信进行频谱分析(fnx)，从中鉴别出有用的信号。要求：将信号的幅度换算成实际

52、的幅度，信号的频率换算成实际的频率。相关(xinggun)Matlab程序如下fs=1000;%采样频率N=512;%数据点数n=0:N-1;randn(1,N);t=0:1/fs:(N-1)/fs;%采样时间序列f0=100;%信号频率x=2+3*cos(100*pi*t-pi/6)+(3/2)*cos(150*pi*t+pi/2);subplot(3,1,1);plot(t,x);xlabel(t);ylabel(sin(2+3*cos(100*pi*t-pi/6)+(3/2)*cos(150*pi*t+pi/2);title(时域信号);Y=fft(x,N);%对信号进行FFT变换mag

53、Y=abs(Y);%求得FFT变换后的幅度angY=angle(Y)*180/pi;%求得FFT变换后的相位f=n*fs/N;%频率序列subplot(3,1,2);plot(f,magY);%画出幅频响应xlabel(f);ylabel(幅度);title(N=512);grid;subplot(3,1,3);plot(f,angY);%画出相频响应(xingyng)xlabel(f);ylabel(相位(xingwi)title(N=512);grid;图4-8 滤波(lb)信号的相位，幅度4.5 利用FFT提取离散信号中的最强正弦分量这里最强是指在信号中某个正弦分量的振幅远大于其它正弦分

54、量的振幅。可以对求点来确定信号中是否有最强的正弦分量。如果信号是连续时间形式的，首先还要对其进行抽样，得到离散时间形式的信号，根据Nyquist定理，抽样间隔应满，其中是中的最大频率分量。要判断信号中是否包含最强正弦分量，采样数据至少要包含该分量一个完整周期的数据15。例4-5 太阳耀斑数据的分析太阳耀斑活动的周期是11年，这个事实可以通过提取耀斑数据的最强正弦分量加以证实。耀斑数据可以从比利时皇家(hungji)天文台(Royal Observatory of Belgium)太阳耀斑数据索引中心（Sunspot Index Data Center, SIDC）网站下载。网址是：http:/

55、sidc.oma.be/index.php3.下载(xi zi)后的数据存放在文件“monthssn.dat”中，里面有四列数据，第一年是日期，第三列是太阳耀斑的平均数，第四列平滑后太阳耀斑的平均数，可以得到从1749年到当前年月（2006年4月）的耀斑数据。本次分析选择(xunz)第三列1981年1月作为开始日期，2005年12月作为结束日期，共25年300个月份的数据。为此先把相关数据复制到Excel表格的第一列中，然后保存到Matlab所在目录下，并命名为“sunspotdata.csv”。然后输入以下命令，得到耀斑曲线如图9所示16。spd=csvread(sunspotdata.cs

56、v,0,0,0 0 299 0);plot(spd);grid;xlabel(月数);ylabel(耀斑平均数);axis(0 300 0 200)图4-9 1981年1月至2005年12月太阳耀斑的平均数由上图可见，太阳耀斑的活动确实具有周期性，但周期的准确值不明显。可以通过数第一个峰值和第二个峰值之间的月份来估计周期的值。查验表中的数据得，第一个峰值为200.3，出现在第116个月（1990年8月），第二个峰值为170.1，出现在第235个月（2000年7月），所以周期是235-116=119月，和实际值132月比较接近。下面(xi mian)利用来分析(fnx)。首先(shuxin)，从

57、图4-9中可以看出从到整个区间的平均值不为0，为了更方便地分析，需要减去该均值，得到结果如下： （4-12）然后对进行傅里叶变换，得到。在Matlab中输入以下命令，得到的幅度谱的图形如图10所示。x=spd-sum(spd)/300;X=fft(x);k=0:299;plot(k,abs(X);grid;xlabel(k);ylabel(|Xk|)图4-10 的幅度(fd)谱由图4-10可见(kjin)，的幅度谱中有一个清晰(qngx)的尖峰，这就证实中确实包含一个最强正弦分量。为了确定尖峰所对应的频率，用火柴棒图重画当时的图形，在Matlab中输入以下命令，得到图形如图4-11所示。k=0

58、:10;X=X(1:11);stem(k,abs(X);xlabel(k);ylabel(|Xk|)图4-11 当时的火柴棒图可以看出上图中有两个强谱分量，频率分别为弧度/月和，周期分别为150月和100月。由于的数据长度不是其中强正弦分量的整数倍，谱分量在出现了泄露。要消除泄露，需要使数据的长度正好是其中强正弦分量的周期（太阳耀斑活动的周期，也即）的整数倍。为此，重新取中从1至264之间的数据进行分析，令 （4-13）然后(rnhu)对进行(jnxng)傅里叶变换，得到。输入以下命令(mng lng)，得到时的火柴棒图如图4-12所示。y=spd(1:264);z=y-sum(y)/264;

59、Z=fft(z);k=0:10;Z=Z(1:11);stem(k,abs(Z);xlabel(k);ylabel(|Zk|)图4-12 当时的火柴棒图注意到的峰值仍然在，但此时比其附近的值大得多，这说明没有出现泄露，频率为弧度/月的正弦分量就是或者中唯一的最强正弦分量。该频率对应周期月，正好等于11年。时的值主要是由于太阳耀斑数据中的噪声引起的。然而和其它值相比，、和的值相对较大，这说明除了周期为11年的最强分量，还包括其它正弦分量，此时太阳耀斑的活动并不是一个只包含单一频率的纯正弦运动。5. 快速(kui s)傅里叶变换在数字信号分析(fnx)与处理的实际(shj)应用5.1 快速傅里叶变换

60、在喇曼光谱信号噪声平滑中的应用电探测系统是光信号的转换、传输及处理的系统. 系统的各个部分在工作时总会受到一些无用信号的干扰,给光谱峰的检测判别及进一步的数据处理带来了不利因素.对光谱信号进行数字滤波,以获得更真实的光谱信息，显得格外重要. 目前最为通用和有效的信号滤波处理方法是快速傅立叶变换方法.纯水是一种较弱的喇曼散射介质,需要专用的喇曼散射光谱仪器才能获得高信噪比的喇曼光谱.我们以增强型的CCD 探头为探测器,结合普通的分光单色仪,在YAG 激光器532 nm 激光线的激励下获得低信噪比的纯水的喇曼光谱. 信噪比较差的喇曼光谱经过FFT 变换后,用FFT 的逆变换将滤除噪声后的频谱信号转