畢業(yè)論文-快速傅里葉變換算法及其在信號處理中的應(yīng)用

1、武漢工程大學(xué)畢業(yè)設(shè)計(論文)說明書 2015 屆畢業(yè)設(shè)計(論文) 題 目快速傅里葉變換算法及其在信號處理中的應(yīng)用專 業(yè) 班 級2011電子信息工程02 學(xué) 號1104030231 姓 名周汝耀 指 導(dǎo) 教 師華夏講師 學(xué) 院 名 稱電氣信息學(xué)院 2011 年 6 月 9 日摘 要快速傅氏變換（FFT），是離散傅氏變換的快速算法，它是根據(jù)離散傅氏變換的奇、偶、虛、實等特性，對離散傅立葉變換的算法進(jìn)行改進(jìn)獲得的。它對傅氏變換的理論并沒有新的發(fā)現(xiàn)，但是對于在計算機(jī)系統(tǒng)或者說數(shù)字系統(tǒng)中應(yīng)用離散傅立葉變換，可以說是進(jìn)了一大步。傅里葉變換的理論與方法在“數(shù)理方程”、“線性系統(tǒng)分析”、“信號處理、仿真”等很

2、多學(xué)科領(lǐng)域都有著廣泛應(yīng)用,由于計算機(jī)只能處理有限長度的離散的序列,所以真正在計算機(jī)上運(yùn)算的是一種離散傅里葉變換.雖然傅里葉運(yùn)算在各方面計算中有著重要的作用，但是它的計算過于復(fù)雜，大量的計算對于系統(tǒng)的運(yùn)算負(fù)擔(dān)過于龐大，使得一些對于耗電量少，運(yùn)算速度慢的系統(tǒng)對其敬而遠(yuǎn)之，然而，快速傅里葉變換的產(chǎn)生，使得傅里葉變換大為簡化，在不犧牲耗電量的條件下提高了系統(tǒng)的運(yùn)算速度，增強(qiáng)了系統(tǒng)的綜合能力，提高了運(yùn)算速度，因此快速傅里葉變換在生產(chǎn)和生活中都有著非常重要的作用，對于學(xué)習(xí)掌握都有著非常大的意義。 關(guān)鍵詞：快速傅氏變換；快速算法；簡化；廣泛應(yīng)用AbstractFast Fourier Transform (

3、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 has not found a new, but in the computer sy

4、stem 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 have a wide range of applications, as the com

5、puter 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 complicated, a lot of computing system for calculati

6、ng 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 increase the speed of computing systems,

7、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 algorithm; simplified; widely used 目 錄摘要.IAbstr

8、act . II1.緒論1.1 選題背景11.2 課題研究的意義22.快速傅里葉變換原理及性質(zhì)2.1快速傅里葉變換原理32.2快速傅里葉變換的優(yōu)越性42.3快速傅里葉變換的意義43.快速傅里葉變換的算法3.1快速傅里葉變換算法63.2 Cooley=Tukey FFT算法83.3 Rader-Brennr FFT算法93.4 Goertsel 算法104.快速傅里葉變換在信號處理中的理論應(yīng)用4.1利用FFT計算連續(xù)時間信號的傅里葉變換134.2利用FFT計算離散信號的線性卷積174.3利用FFT進(jìn)行離散信號壓縮194.4利用FFT對離散信號進(jìn)行濾波224.5利用FFT提取離散信號中的最強(qiáng)正弦分

9、量245.快速傅里葉變換在數(shù)字信號分析與處理的實際應(yīng)用5.1快速傅里葉變換在喇曼光譜信號噪聲平滑中的應(yīng)用295.2采用異步實現(xiàn)的快速傅里葉變換處理器315.3快速傅里葉算法在哈特曼夏克傳感器波前重構(gòu)算法中的應(yīng)用33致謝36參考文獻(xiàn)37 1 緒論傅立葉變換在生產(chǎn)生活中的重要性非常突出，它將原本難以處理的時域信號相比比較容易地轉(zhuǎn)換成易于分析的頻域信號，我們可以利用一些專業(yè)工具對這些頻域信號進(jìn)行處理、加工，使信號轉(zhuǎn)化為可以對其進(jìn)行各種數(shù)學(xué)變換的數(shù)學(xué)公式，對其進(jìn)行處理。最后還可以根據(jù)傅立葉反變換將這些頻域信號轉(zhuǎn)換成原來的時域信號，這是一種特殊的積分變換。它能夠?qū)M足一定條件的某個函數(shù)表示成為正弦基函數(shù)

10、的線性組合或者積分。然而，它在運(yùn)算上過于復(fù)雜，過于宏大的運(yùn)算過程，對于一些相對簡單的低功耗處理器來說，難以自如應(yīng)對，因此，快速傅里葉變換則顯出了它的優(yōu)越性?？焖俑凳献儞Q（FFT），即離散傅氏變換的快速算法，它是根據(jù)離散傅氏變換的奇、偶、虛、實等特性，對離散傅立葉變換的算法進(jìn)行改進(jìn)獲得的1。對于計算機(jī)處理信號方面上是一大進(jìn)步。系統(tǒng)的速度不但取決于其本身的速度，而且還在相當(dāng)大的程度上取決于運(yùn)用的算法，算法運(yùn)算量的大小直接影響到對設(shè)備的控制質(zhì)量。通過傅立葉變換（DFT），運(yùn)用測試軟件進(jìn)行檢測，我們可以看出，快速傅里葉變換大大的提高了運(yùn)算速度，它為各系統(tǒng)的設(shè)計方案提供了簡單算法，有著非常重要的意義。1

11、. 1 選題背景近十多年來，數(shù)字信號處理技術(shù)同大規(guī)模集成電路、數(shù)字計算機(jī)等，都有了突飛猛進(jìn)的發(fā)展，日新月異，早已成為了一門具有強(qiáng)大生命力的技術(shù)科學(xué)。因為它本身就具有一系列的優(yōu)點(diǎn)，所以能夠有效地促進(jìn)工程技術(shù)領(lǐng)域的技術(shù)改造和科學(xué)發(fā)展，應(yīng)用領(lǐng)域也更加地廣泛、深入，越來越受到人們的重視。在信號處理中，離散傅里葉變換（Discrete Fourier Transform，DFT）是比較常用的變換方法之一，它在各種數(shù)字信號處理系統(tǒng)中扮演著及其重要的角色。由于離散傅里葉變換（DFT）而發(fā)現(xiàn)了頻率離散化，可以直接用它來分析信號的頻譜、計算濾波器的頻率響應(yīng)、以及實現(xiàn)信號通過線系統(tǒng)的卷積運(yùn)算等，因而在信號的譜分析

12、等方面有著非常大的作用。傅里葉變換已經(jīng)有一百多年的歷史了，我們熟知頻域分析往往比時域分析更優(yōu)越，不僅簡單明了，而且易于分析較為復(fù)雜的信號。但需要用較為精準(zhǔn)的數(shù)字方法，即DFT，進(jìn)行譜分析，在快速傅氏變換（FFT）出現(xiàn)以前是不切實際的。由于DFT的計算量太大，即使運(yùn)用計算機(jī)也很難對問題進(jìn)行實時的有效處理，所以DFT并沒有得到真正的應(yīng)用。直到1965年庫利（J.W.Cooly）和圖基（J.W.Tukey）首次發(fā)現(xiàn)DFT的一種快速算法，局面才發(fā)生根本性的變化。繼庫利和圖基算法出來之后，桑德（G.Sander）等快速算法也相繼出現(xiàn)，又經(jīng)過其他學(xué)者一步步改進(jìn)，很快就出現(xiàn)了通用型的快速傅里葉變換，簡稱FF

13、T?？焖俑道锶~變換（Fast Fourier Transform，F(xiàn)FT）并非與離散傅里葉變換完全不同的另一種變換，而是為了減少DFT計算次數(shù)而誕生的一種快速、有效的算法。應(yīng)當(dāng)指出的是，也是因為當(dāng)時電子數(shù)字計算機(jī)的“落后”條件也促成了這個算法的提出。它使得DFT的運(yùn)算量大大的縮小簡化，它推動了近30年來信號處理技術(shù)止步不前的前進(jìn)發(fā)展，成為了數(shù)字信號處理應(yīng)用領(lǐng)域里強(qiáng)有力的工具，為DFT乃至數(shù)字信號處理技術(shù)的實際應(yīng)用創(chuàng)造了良好的條件，從而使DFT在實際使用中得以廣泛的應(yīng)用2。數(shù)字信號處理器（DSP），是一種可編程的高性能處理器。近年來發(fā)展尤為迅速，它不僅應(yīng)用于數(shù)字信號處理方面，而且在圖像處理、語音

14、處理、通信等領(lǐng)域得到廣泛的應(yīng)用。之前通用的微處理器在運(yùn)算速度上已經(jīng)很難適應(yīng)信號實時處理的高要求。DSP處理器中集成了高速的乘法硬件，能快速、準(zhǔn)確地進(jìn)行大量數(shù)據(jù)的乘法以及加法的運(yùn)算。數(shù)字信號處理區(qū)別于普通的科學(xué)計算與分析，它強(qiáng)調(diào)運(yùn)算的實時性。除了需要普通微處理器所強(qiáng)調(diào)的高速運(yùn)算和控制能力之外，鑒于實時數(shù)字信號處理的特點(diǎn)，在處理器結(jié)構(gòu)、指令系統(tǒng)、指令流程上做了很大程度上的改進(jìn)。1. 2 課題研究的意義如上所述，基于對DSP的快速傅里葉變換算法的研究，從而使FFT算法能夠有效地在DSP芯片上實現(xiàn)。DSP芯片的出現(xiàn)，使FFT的實現(xiàn)更加方便。多數(shù)的DSP芯片都能夠在一個指令周期內(nèi)完成一次乘法和加法，并且

15、提供了專門的FFT指令，完成一次指令的周期只需10ns，使得FFT算法在DSP芯片上實現(xiàn)的速度更加快速?？焖俑道锶~變換為頻譜分析、卷積、相關(guān)數(shù)字濾波器設(shè)計與實現(xiàn)與功率譜計算、傳遞函數(shù)建模、圖像處理等，提供了快速有效的運(yùn)算方法。FFT技術(shù)應(yīng)用DSP芯片，從而可以提供使調(diào)制、解調(diào)、壓縮、解壓縮的數(shù)據(jù)傳輸更為高效的信號處理解決方案，因而廣泛應(yīng)用于雷達(dá)、圖像處理、通信、生物醫(yī)學(xué)和聲納領(lǐng)域。2.快速傅里葉變換原理及性質(zhì) 數(shù)字信號中的傅里葉變換,通常是采用離散型傅里葉變換(DFT)。DFT 存在的缺點(diǎn)就是計算量太大,不易進(jìn)行實時處理。比如，計算一個N 點(diǎn)的DFT ,一般需要次復(fù)數(shù)乘法和N(N-1)次復(fù)數(shù)加

16、法運(yùn)算.因此,當(dāng)N較大或要求對信號進(jìn)行實時處理時,往往很難實現(xiàn)達(dá)到所需的運(yùn)算速度。1965年,J.W.Cooly和J.W.Tukey發(fā)現(xiàn)了DFT的一種快速算法,經(jīng)過后來學(xué)者的進(jìn)一步改進(jìn), 很快便形成了一套高效的運(yùn)算方法,即現(xiàn)在通用的快速傅里葉變換, 簡稱FFT( The Fast Fourier Transform)?？焖俑道锶~變換的實質(zhì)是利用式(3-1)中的權(quán)函數(shù)的對稱性和周期性,把N點(diǎn)DFT進(jìn)行一系列分解和組合,使整個DFT的計算過程變成一系列疊代運(yùn)算過程,使DFT的運(yùn)算量大大簡化,為DFT及數(shù)字信號的實時處理和應(yīng)用創(chuàng)造了非常良好的條件3。2. 1 快速傅里葉變換原理快速傅里葉變換原理：1

17、. 將長序列DFT分解為短序列的DFT2. 利用旋轉(zhuǎn)因子的對稱性、周期性、可約性。將時域序列逐次分解為一組子序列，依據(jù)旋轉(zhuǎn)因子的特性，由子序列的DFT來實現(xiàn)整個序列的DFT4。其中：快速傅里葉變換分為兩種，分為基2時間抽取算法和基2頻率抽取算法基2-時間抽取(Decimation in time)FFT算法其中:r=0,1,2 （2-1）基2-頻率抽取(Decimation in frequency)FFT算法 （2-2）2. 2 快速傅里葉變換的優(yōu)越性設(shè)為項的復(fù)數(shù)序列，依據(jù)DFT變換，任一的計算都需要有次復(fù)數(shù)乘法和（）次復(fù)數(shù)加法，而且一次復(fù)數(shù)乘法等同于四次實數(shù)乘法和兩次實數(shù)加法，同樣的，一次

18、復(fù)數(shù)加法等同于兩次實數(shù)加法，即使我們把一次復(fù)數(shù)乘法和一次復(fù)數(shù)加法定義成一次“運(yùn)算”（四次實數(shù)乘法和四次實數(shù)加法），那么求出項復(fù)數(shù)序列的, 即N點(diǎn)DFT變換大約就需要次運(yùn)算。當(dāng)點(diǎn)甚至更多的時候，需要N2=1048576次運(yùn)算。在FFT中，利用WN的對稱性和周期性，把一個N項序列（設(shè),為正整數(shù)），分為兩個項的子序列，而且每個點(diǎn)的DFT變換需要次運(yùn)算，再運(yùn)用N次運(yùn)算把兩個點(diǎn)的DFT 變換重新組合成一個N點(diǎn)的DFT變換。如此變換以后，總的運(yùn)算次數(shù)就變成了。承接上面的例子，當(dāng)時，總的運(yùn)算次數(shù)就變成了525312次，這樣看來，節(jié)省了大約50%的運(yùn)算量。而如果我們將這種“一分為二”的思想不斷進(jìn)行下去，直到分

19、成兩兩一組的DFT運(yùn)算單元，那么N點(diǎn)的 DFT變換就只需要次的運(yùn)算，在1024點(diǎn)時，運(yùn)算量僅有10240次，是先前的直接算法的1%，點(diǎn)數(shù)越多，運(yùn)算量的節(jié)約就越大，這就是 FFT的優(yōu)越性.當(dāng)然，F(xiàn)FT提高了運(yùn)算速度,但是,也對參與運(yùn)算的樣本序列作出了限制,即要求樣本數(shù)為2N點(diǎn).離散傅里葉變換DFT則無上述限制5。2. 3 快速傅里葉變換的意義 傅立葉變換的物理意義：傅立葉變換是數(shù)字信號處理技術(shù)領(lǐng)域一項很重要的算法。要知道傅立葉變換算法的意義，首先要了解傅立葉變換原理的意義。傅立葉變換原理表明：任何連續(xù)測量的時序或信號，都能夠表示成為不同頻率的正弦波信號的無限疊加。而利用該原理而創(chuàng)立的傅立葉變換算

20、法則利用直接能測量到的原始信號，并以累加方式來計算該信號中不同正弦波信號的頻率、相位和振幅。與傅立葉變換算法對應(yīng)的是反傅立葉變換算法。該反變換從本質(zhì)上說也就是一種累加處理，這樣便可以將單獨(dú)改變的正弦波信號轉(zhuǎn)換成一個信號。因此，也可以說，傅立葉變換將原來難以處理的時域信號轉(zhuǎn)換成了易于分析的頻域信號（信號的頻譜），我們可以利用一些專業(yè)工具對這些頻域信號進(jìn)行加工、處理。最后還可以利用傅立葉反變換將這些頻域信號轉(zhuǎn)換成原來的時域信號。從現(xiàn)代數(shù)學(xué)的眼光來看，傅里葉變換其實就是一種特殊的積分變換。它能夠?qū)M足一定條件下的某個函數(shù)表示成為正弦基函數(shù)的線性組合或者積分形式。在不同的研究領(lǐng)域里，傅里葉變換具有多種

21、形式各異的變體形式，如連續(xù)傅里葉變換和離散傅里葉變換。在數(shù)學(xué)領(lǐng)域，盡管最初傅立葉分析是作為熱過程的解析分析的工具，但是其思想方法仍然具有典型的還原論和分析主義的特征。任意的函數(shù)通過一定的分解，都能夠表示為正弦函數(shù)的線性組合的形式，而正弦函數(shù)在物理上是被充分研究而相對簡單的函數(shù)類：1. 傅立葉變換是線性算子,若賦予適當(dāng)?shù)姆稊?shù),它還是酉算子;2. 傅立葉變換的逆變換容易求出,而且形式與正變換非常類似;3. 正弦基函數(shù)是微分運(yùn)算的本征函數(shù),從而使得線性微分方程的求解可以轉(zhuǎn)化為常系數(shù)的代數(shù)方程的求解.在線性時不變雜的卷積運(yùn)算為簡單的乘積運(yùn)算,從而提供了計算卷積的一種簡單手段;4. 離散形式的傅立葉的物

22、理系統(tǒng)內(nèi),頻率是個不變的性質(zhì),從而系統(tǒng)對于復(fù)雜激勵的響應(yīng)可以通過組合其對不同頻率正弦信號的響應(yīng)來獲取;5. 著名的卷積定理指出:傅立葉變換可以化復(fù)變換可以利用數(shù)字計算機(jī)快速的算出(其算法稱為快速傅立葉變換算法(FFT)。正是由于上述的良好性質(zhì),傅里葉變換在物理學(xué)、數(shù)論、組合數(shù)學(xué)、信號處理、概率、統(tǒng)計、密碼學(xué)、聲學(xué)、光學(xué)等領(lǐng)域都有著廣泛的應(yīng)用6。3. 快速傅里葉變換的算法3.1 快速傅里葉變換算法 快速傅里葉變換FFT 是離散傅里葉變換DFT 的一種快速算法，只有FFT 才能在現(xiàn)實中有實際應(yīng)用的意義。 (3-1)由(3-1)式可知，對每一個n，計算X()須作N次復(fù)數(shù)乘法及N-1次復(fù)數(shù)加法，要完成

23、這組變換共需次乘法及N(N-1)次復(fù)數(shù)加法。但以下介紹的快速傅里葉變換的算法，可大大減少運(yùn)算次數(shù)，提高工作效率。當(dāng)時，n和k可用二進(jìn)制數(shù)表示： (3-2) （3-3）又記,則（3-1）式可改寫為 （3-4）式中： (3-5)因為所以（3-2）可改成 （3-6） (3-7)則式（3-7）即為式（3-6）的分解形式。將初始數(shù)據(jù)代入式（3-7）的第一個等式，可得每一組計算數(shù)據(jù)，一般將痗L-1組計算數(shù)據(jù)代入式（3-7）的第L個等式，計算后可得第L組計算數(shù)據(jù)（L，），計算公式也可表示為 = （3-8）式中 （3-9） 根據(jù)式（3-8），第L個數(shù)組中每個 的計算只依賴于上一個數(shù)組的兩個數(shù)據(jù)這兩個數(shù)據(jù)的標(biāo)號

24、相差，即，而且這兩個數(shù)據(jù)只用于計算第L個數(shù)組中標(biāo)號的數(shù)據(jù)（等號右端為二進(jìn)制數(shù)）。當(dāng)分別取和時，分別有。因此，用上一組的兩個數(shù)據(jù)計算所得的兩個新數(shù)據(jù)仍可儲存在原來位置，計算過程中只需要N個存儲器。將與稱為第L個數(shù)組中的對偶結(jié)點(diǎn)對。計算每個對偶結(jié)點(diǎn)對只需一次乘法，事實上由式（3-8）可得 （3-10）式中： ；別為式（3-9）中取，時對應(yīng)的P值。因，于是對偶結(jié)點(diǎn)的有如下關(guān)系： ,因此式（3-8）可表示為 （3-11） P的求法：在中，i寫成二進(jìn)制數(shù)右移位，就成為顛倒位序得式（3-7）中，前面的r個等式，每個等式均對應(yīng)一組數(shù)據(jù)進(jìn)行計算，每組數(shù)據(jù)都有N/對結(jié)點(diǎn)，根據(jù)式（3-11），每對結(jié)點(diǎn)只需作次乘法

25、和次加法，因此，每組數(shù)據(jù)只需N/2次乘法和N次加法，因而完成r組數(shù)據(jù)的計算共需Nr/2次乘法和Nr次加法。3. 2 Cooley-Tukey FFT算法的核心是將一層運(yùn)算映射為二層運(yùn)算。作點(diǎn)時，若不是素數(shù)，則可分解為，那么由的 （3-12）通過映射： （3-13）可得到 而，可化簡為 （3-14）從而式（3-12）轉(zhuǎn)化為 （3-15）其中。以點(diǎn)為例，映射方式為：，則計算流圖如圖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

26、 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 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 點(diǎn)FFT算法3. 3 Rader-Brenner FFT算法Rader-Brenner 算法是類似于 Cooley-Tukey 算法，但是采用的旋轉(zhuǎn)因子都是純

27、虛數(shù)，以增加加法運(yùn)算和降低了數(shù)值穩(wěn)定性為代價減少了乘法運(yùn)算。這方法之后被split-radix variant of Cooley-Tukey所取代，與Rader-Brenner算法相比，有一樣多的乘法量，卻有較少的加法量，且不犧牲數(shù)值的準(zhǔn)確性7。 Bruun以及QFT算法是不斷的把DFT分成許多較小的DFT運(yùn)算。(Rader-Brenner以及QFT算法是為了2的指數(shù)所設(shè)計的算法，但依然可以適用在可分解的整數(shù)上。Bruun算法則可以運(yùn)用在可被分成偶數(shù)個運(yùn)算的數(shù)字)。尤其是Bruun算法，把FFT看成是 ，并把它分解成 與 的形式。 另一個從多項式觀點(diǎn)的快速傅立葉變換法是Winograd 算法

28、。此算法把 分 解成cyclotomic多項式，而這些多項式的系數(shù)通常為1，0，-1。 這樣只需要很少的乘法量(如果有需要的話)，所以winograd是可以得到最少乘法量的快速傅立葉算法，對于較小的數(shù)字，可以找出有效率的算方式。更精確地說，winograd算法讓DFT可以用2K點(diǎn)的DFT來簡化，但減少乘法量的同時，也增加了非常多的加法量。Winograd也可以利用剩余值定理來簡化DFT。 Rader算法提出了利用點(diǎn)數(shù)為N(N為質(zhì)數(shù))的DFT進(jìn)行長度為N-1的回旋折積來表示原本的DFT，如此就可利用折積用一對基本的FFT來計算 DFT。另一個prime-size的FFT算法為chirp-Z算法。

29、此法也是將DFT用折積來表示，此法與Rader算法相比，能運(yùn)用在更一般的轉(zhuǎn)換 上，其轉(zhuǎn)換的基礎(chǔ)為Z轉(zhuǎn)換8。3.4 Goertsel 算法如前所述, 點(diǎn)時域序列的離散付里葉變換式為， （3-16）這點(diǎn)頻域序列是同時被算出的，不可能只計算其中某一個或幾個指定點(diǎn)。Goetzel 算法是為了解決這個問題而提出的。這個算法把離散付里葉變換看作一組濾波器，將輸入端的時域序列與其中一個濾波器的沖激響應(yīng)序列進(jìn)行卷積運(yùn)算，求濾波器的輸出序列，即得序列的一點(diǎn)。這種算法利用旋轉(zhuǎn)因子的周期性，使DFT運(yùn)算化為線性濾波運(yùn)算9。 由于 故式（3-16）可化為 （3-17）定義序列為 （3-18）可見是由兩個序列卷積而得到

30、的序列。 （3-19）其中，是輸入的點(diǎn)序列，另一個序列被看作濾波器的沖激響應(yīng)序列。 （3-20）對比式(3-17）和式(3-18），可知：按式（3-19）進(jìn)行卷積運(yùn)算，當(dāng)時，濾波器的輸出就是： （3-21）對式(3-20)進(jìn)行Z變換，可得濾波器的系統(tǒng)函數(shù) （3-22）這是一個一階系統(tǒng)。圖3-2示出這個系統(tǒng)的信號流圖，相應(yīng)的差分方程為 ， (3-23) 按照此式進(jìn)行遞推運(yùn)算，到了 時刻，即可依據(jù)式（3-21）得到。 按照式(3-22)進(jìn)行運(yùn)算時，可先算好旋轉(zhuǎn)因子，儲存起來。每次遞推包含一次復(fù)數(shù)乘法。按式(3-16)直接計算點(diǎn)離散付里葉變換,需要次實數(shù)乘法和次實數(shù)加法。按照上述Goertzel算法

31、，所需的實數(shù)乘法和實數(shù)加法都是次。所以當(dāng)不大時，上述算法的效率稍差10。下面介紹改進(jìn)的Goertzel 算法，這種算法所需的實數(shù)乘法次數(shù)約為直接方法的一半。圖3-2 用一階系統(tǒng)實現(xiàn)Goertzel 算法 圖 3-3 用二階系統(tǒng)實現(xiàn)Goert算法把式（3-22）的分子和分母都乘上因子，就得到第個濾波器的系統(tǒng)函數(shù)為 （3-24）與此相應(yīng)的信號流圖示于圖3-3。由式（3-24)可見，濾波器是一個二階系統(tǒng)，有一對復(fù)數(shù)共軛極點(diǎn)和一個復(fù)數(shù)零點(diǎn)。為了便于運(yùn)算，在圖3-3所示的流圖中，設(shè)立狀態(tài)變量 和。按照圖3-3計算時，步驟有二，即： 1.實現(xiàn)一對復(fù)數(shù)極點(diǎn)輸入點(diǎn)依次取 ，進(jìn)行遞推運(yùn)算。每次運(yùn)算中，更新狀態(tài)變

32、量和。作次迭代所需的計算量是次實數(shù)乘法和次實數(shù)加法。 2.實現(xiàn)復(fù)數(shù)零點(diǎn)。是一個點(diǎn)序列， 。在點(diǎn)上。計算狀態(tài)變量和。這時，按照圖3-2算出濾波器的輸出，此即。所需的計算量是4次實數(shù)乘法和次實數(shù)加法11。 綜上所述，計算一點(diǎn)需要進(jìn)行次實數(shù)乘法和次實數(shù)加法。這種算法要求的乘法次數(shù)約為直接算法的一半。在這種較為有效的方案中，仍具有這樣的優(yōu)點(diǎn)，即必須計算和存儲的系數(shù)只有和。還要說明圖3-3所示的算法的另一個優(yōu)點(diǎn)。當(dāng)輸入序列為實序列時，離散付里葉變換序列是對稱的，即。容易證明，圖3-3的網(wǎng)絡(luò)形式在計算時和計算時具有完全相同的極點(diǎn)，但前者的零點(diǎn)系數(shù)與后者的零點(diǎn)系數(shù)成復(fù)共軛關(guān)系。由于零點(diǎn)僅在最后的迭代中實現(xiàn)，

33、所以諸極點(diǎn)要求的次乘法和 次加法可以用來計算離散付里葉變換的兩個值。因此，若用Goertzel 算法計算離散付里葉變換的所有個點(diǎn)的值，需要的乘法次數(shù)近似為，加法次數(shù)近似為。然而，它同直接計算離散付里葉變換一樣，計算量仍然正比于。4.快速傅里葉變換在信號處理中的理論應(yīng)用4.1 連續(xù)時間信號的快速傅里葉變換設(shè)是連續(xù)時間信號，并假設(shè)時，則其傅里葉變換由下式給出： （4-1） 令是一固定的正實數(shù)，是一個固定的正整數(shù)。當(dāng)時，利用FFT算法可計算。已知一個固定的時間間隔，選擇讓T足夠小，使得每一個T秒的間隔內(nèi)，的變化很小，則式中積分可近似為 （4-2）假設(shè)足夠大，對于所有的整數(shù)，幅值很小，則式（4-2）變

34、為 （4-3）當(dāng)時，式（4-2）兩邊的值為 （4-4）其中代表抽樣信號的點(diǎn)。最后令，則上式變?yōu)?（4-5）首先用FFT算法求出，然后可用上式求出時的。應(yīng)該強(qiáng)調(diào)的是，式（4-3）只是一個近似表示，計算得到的只是一個近似值。通過取更小的抽樣間隔，或者增加點(diǎn)數(shù)，可以得到更精確的值。如果時，幅度譜很小，對應(yīng)于奈奎斯特抽樣頻率，抽樣間隔選擇比較合適。如果已知信號只在時間區(qū)間內(nèi)存在，可以通過對時的抽樣信號補(bǔ)零，使足夠大12。將連續(xù)時間傅立葉變換進(jìn)行數(shù)字近似，用函數(shù)fft（快速傅立葉算法）高效地計算這個近似值。很多信號都能用（4-1）式連續(xù)時間傅立葉變換（CTFT）來表示。利用MATLAB可以計算（CTFT

35、）積分的數(shù)值近似。利用在密集的等間隔t的樣本上的求和來近似這個積分，就可以用函數(shù)fft高效地計算這個近似值。所用的近似式是根據(jù)積分的定義得到的，即 （4-6）對于一般信號，在足夠小的下，上式右邊的和式是對于CTFT積分的一個好的近似。若信號對于和為零，那么這個近似式就能寫成 （4-7）式中，N為一整數(shù)?？梢岳煤瘮?shù)fft對一組離散的頻率計算上式中的和式。如果N個樣本是存在向量x內(nèi)的話，那么調(diào)用函數(shù)X=tau*fft(x)就可以計算出 （4-8）式中 以及N假設(shè)為偶數(shù)。為了計算高效，fft在負(fù)的頻率樣本之前先產(chǎn)生正頻率樣本。為了將頻率樣本置于上升的順序，能用函數(shù)fftshift。為了將存入X中的

36、的樣本排列成使就是對于，在上求得的CTFT,可用X=fftshift(tau*fft(x)。例4-1 利用FFT計算傅里葉變換如圖4-1所示的信號其傅里葉變換為：利用下面的命令，可得到的近似值和準(zhǔn)確值。 圖4-1 連續(xù)時間信號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*k*gamma)*X;%計算真實值X(w)w=0.05:0.05:10

37、;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|)運(yùn)行程序后輸入N=128，T=0.1，此時，得到實際的和近似的傅里葉變換的幅度譜如圖4-2所示，此時近似值已經(jīng)相當(dāng)準(zhǔn)確。通過增加NT可以增加更多的細(xì)節(jié)，減少T使得到的值更精確。再次運(yùn)行程序后輸入N=512，T=0.05，此時，得到實際的和近似的傅里葉變換的幅度譜如圖4-3所示。圖4-2 N=128，T=0.1時的幅度譜

38、圖4-3 N=512，T=0.05時的幅度譜4.2 FFT計算離散信號的線性卷積已知兩個離散時間信號與，取，對和右端補(bǔ)零，使得 （4-9）利用FFT算法可以求得和的L點(diǎn)DFT，分別是和，利用DTFT卷積性質(zhì)，卷積等于乘積的L點(diǎn)DFT反變換，這也可以通過FFT算法得到。例4-2 利用FFT計算線性卷積已知，其中為單位階躍序列，信號如圖4-4所示。由于當(dāng)時，很小，故可以取為17；N取10，。利用下面的Matlab命令，可得到、的卷積圖形如圖4-4所示。subplot(3,1,1);n=0:16;x=0.8.n;stem(n,x);xlabel(n);ylabel(xn);subplot(3,1,2

39、);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 信號xn、yn及其卷積zn=xn*yn利用下面的Matlab命令，可得到信號xn、yn的幅度譜與相位譜如圖4-5所示。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(

40、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)(弧度)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,angle(Y);axis(0 25 -3 3);xlabel(k);ylabel(Angle(Yk)(弧度)圖4-5 信號xn、yn的幅度譜與相位譜4.3

41、 FFT進(jìn)行離散信號壓縮利用FFT算法對離散信號進(jìn)行壓縮的步驟如下：1）通過采樣將信號離散化；2）對離散化信號進(jìn)行傅里葉變換；3）對變換后的系數(shù)進(jìn)行處理，將絕對值小于某一閾值的系數(shù)置為0，保留剩余的系數(shù)；4）利用IFFT算法對處理后的信號進(jìn)行逆傅里葉變換13。 例4-3 對單位區(qū)間上的下列連續(xù)信號以采樣頻率進(jìn)行采樣，將其離散化為個采樣值用FFT分解信號，對信號進(jìn)行小波壓縮，然后重構(gòu)信號。令絕對值最小的80%系數(shù)為0，得到重構(gòu)信號圖形如圖4-6a)所示，均方差為0.0429，相對誤差為0.0449；令絕對值最小的90%系數(shù)為0，得到重構(gòu)信號圖形如圖4-6b)所示，均方差為0.0610，相對誤差為

42、0.0638。 a) 絕對值最小的80%系數(shù)為0的重構(gòu)信號（FFT） b) 絕對值最小的90%系數(shù)為0的重構(gòu)信號（FFT）圖4-6 用FFT壓縮后的重構(gòu)信號相關(guān)Matlab程序如下:function wc=compress(w,r)%壓縮函數(shù)compress.m%輸入信號數(shù)據(jù)w,壓縮率r%輸出壓縮后的信號數(shù)據(jù)if(r1) error(r 應(yīng)該介于0和1之間!);end;N=length(w);Nr=floor(N*r);ww=sort(abs(w);tol=abs(ww(Nr+1);wc=(abs(w)=tol).*w; function unbiased_variance,error=fft

43、comp(t,y,r)%利用FFT做離散信號壓縮%輸入時間t,原信號y,以及壓縮率rif(r1) error(r 應(yīng)該介于0和1之間!);end;fy=fft(y);fyc=compress(fy,r); %調(diào)用壓縮函數(shù)compress.myc=ifft(fyc);plot(t,y,r,t,yc,b);legend(原信號,重構(gòu)信號);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);unb

44、iased_variance,error=fftcomp(t,f,0.8) unbiased_variance = 0.0429error =0.0449如果用Harr尺度函數(shù)和Harr小波分解信號，對信號進(jìn)行小波壓縮，然后重構(gòu)信號。令絕對值最小的80%系數(shù)為0，得到重構(gòu)信號圖形如圖4-7a)所示，均方差為0.0584，相對誤差為0.0611；令絕對值最小的90%系數(shù)為0，得到重構(gòu)信號圖形如圖4-7b)所示，均方差為0.1136，相對誤差為0.1190。 a) 絕對值最小的80%系數(shù)為0的重構(gòu)信號（Harr） b) 絕對值最小的90%系數(shù)為0的重構(gòu)信號（Harr）圖4-7 用Harr小波壓縮后

45、的重構(gòu)信號相關(guān)Matlab程序如下function unbiased_variance,error=daubcomp(t,y,n,r)%利用Daubechies系列小波做離散信號壓縮%輸入時間t,原信號y,分解層數(shù)n,以及壓縮率r%輸出原信號和壓縮后重構(gòu)信號的圖像,以及重構(gòu)均方差和相對l2誤差if(r1) error(r應(yīng)該介于0和1之間!);end;c,l=wavedec(y,n,db1);cc=compress(c,r); %調(diào)用壓縮函數(shù)compress.myc=waverec(cc,l,db1);plot(t,y,r,t,yc,b);legend(原信號,重構(gòu)信號);unbiased_v

46、ariance=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(t,f,8,0.8) unbiased_variance = 0.0584error = 0.0611結(jié)論：在信號沒有突變、快變化或者大致上具有周期性的信號，用FFT可以處理得很好（甚至比小波還要好）。4. 4 利用FFT對離散信號進(jìn)行濾波利用FFT算法對信號進(jìn)行濾波的步驟如下：1）通過采樣將信號離

47、散化；2）對離散化信號進(jìn)行傅里葉變換；3）對變換后的系數(shù)進(jìn)行處理，將絕對值大于某一閾值的系數(shù)置為0，保留剩余的系數(shù)；4）利用IFFT算法對處理后的信號進(jìn)行逆傅里葉變換14。例4-4 對被白噪聲污染的信進(jìn)行頻譜分析，從中鑒別出有用的信號。要求：將信號的幅度換算成實際的幅度，信號的頻率換算成實際的頻率。相關(guān)Matlab程序如下fs=1000;%采樣頻率N=512;%數(shù)據(jù)點(diǎn)數(shù)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);%對信號進(jìn)行FFT變換magY=abs(Y);%求得FFT