深入淺出數(shù)字信號處理下部（6-9章）目錄\h第6章快速傅里葉變換\h6.1DFT的困境與機遇\h6.1.1DFT面臨的主要問題\h6.1.2DFT高效運算的可能性\h6.2基-2DITFFT算法\h6.2.1算法原理\h6.2.2算法特點\h6.2.3運算量分析\h6.3其他常用的FFT算法\h6.3.1基-2DIFFFT\h6.3.2基-4FFT\h6.3.3分裂基FFT\h6.3.4算法比較\h6.4現(xiàn)實應用中的考慮\h6.4.1旋轉因子的生成\h6.4.2IDFT的快速計算\h6.4.3實信號的FFT\h6.5應用實例及其MATLAB實現(xiàn)\h6.6兩類特殊的頻率分析算法\h6.6.1Goertzel算法\h6.6.2ChirpZ算法\h6.7信號處理雜談：高斯的遺憾\h6.8本章小結\h第7章數(shù)字濾波器概述\h7.1數(shù)字濾波器的基本概念\h7.1.1數(shù)字濾波器與信號處理過程\h7.1.2數(shù)字濾波器的描述\h7.1.3數(shù)字濾波器的分類\h7.2實際數(shù)字濾波器的性能參數(shù)\h7.2.1頻域性能參數(shù)\h7.2.2時域性能參數(shù)\h7.3數(shù)字濾波問題的一般解決方案\h7.3.1濾波器性能參數(shù)的確定\h7.3.2濾波器類型的選擇\h7.3.3濾波器系數(shù)的計算\h7.3.4濾波器的實現(xiàn)結構\h7.3.5有限字長對濾波器的影響\h7.3.6濾波器的軟硬件實現(xiàn)\h7.4數(shù)字濾波器的評估\h7.5從低通濾波器出發(fā)\h7.6信號處理雜談：數(shù)字濾波器的通俗化解釋\h7.7本章小結\h第8章有限沖激響應濾波器\h8.1FIR濾波器的主要特征\h8.1.1基本特征\h8.1.2線性相位特性\h8.2FIR濾波器系數(shù)的計算方法\h8.2.1窗函數(shù)法\h8.2.2頻率采樣法\h8.2.3最優(yōu)化方法\h8.3FIR濾波器的實現(xiàn)結構\h8.3.1橫向結構\h8.3.2線性相位結構\h8.3.3頻率采樣結構\h8.3.4快速卷積結構\h8.3.5實現(xiàn)結構的選擇\h8.4有限字長的影響\h8.4.1系數(shù)的量化誤差\h8.4.2運算的舍入誤差\h8.5兩類特殊的FIR濾波器\h8.5.1滑動平均濾波器\h8.5.2半帶濾波器\h8.6應用實例及其MATLAB實現(xiàn)\h8.7信號處理雜談：最優(yōu)化方法的發(fā)展\h8.8本章小結\h第9章無限沖激響應濾波器\h9.1IIR濾波器的主要特征\h9.2IIR濾波器系數(shù)計算方法\h9.2.1基礎知識\h9.2.2沖激不變法\h9.2.3雙線性變換法\h9.2.4零極點放置法\h9.2.5計算方法的選擇\h9.3IIR濾波器的實現(xiàn)結構\h9.3.1直接型\h9.3.2串聯(lián)型\h9.3.3并聯(lián)型\h9.4有限字長的影響\h9.4.1系數(shù)量化對極點位置的擾動\h9.4.2系數(shù)量化導致的頻率響應誤差\h9.4.3運算舍入誤差的影響\h9.5幾類特殊的IIR濾波器\h9.5.1諧振器\h9.5.2陷波器\h9.5.3梳狀濾波器\h9.6應用實例及其MATLAB實現(xiàn)\h9.7信號處理雜談：IIR濾波器沉浮錄\h9.8本章小結\h附錄AMATLAB語言基礎\hA.1變量\hA.2數(shù)組、向量與矩陣\hA.3部分特殊變量和常數(shù)\hA.4部分常用運算符\hA.5程序結構\hA.6部分基本數(shù)學函數(shù)和作圖函數(shù)\h附錄BMATLAB信號處理工具箱常用函數(shù)\h附錄C分貝（dB）\hC.1分貝的定義\hC.2分貝的優(yōu)點\hC.3一些有用的數(shù)值\h附錄D時域采樣與頻域延拓第6章快速傅里葉變換在前面的討論中已經(jīng)知道，相比傅里葉家族中的其他變換，DFT是實踐性最強的算法，因為它在時間和頻率上都是離散的，非常便于計算機等數(shù)字系統(tǒng)的實現(xiàn)。但DFT在具體的工程實現(xiàn)中是否就一帆風順呢？如果有問題的話，癥結主要在哪里呢？科學家們又是如何解決的呢？這些是本章討論的主要問題。本章在6.1節(jié)簡要討論了DFT在工程實現(xiàn)中面臨的主要問題，并指出了總體的解決思路。6.2節(jié)詳細介紹了按時間抽取的DFT快速算法，即基-2DITFFT算法，不僅討論了算法的原理，而且討論了算法的運算規(guī)律及運算量的減少等問題。6.3節(jié)簡要介紹了其他幾種常用的FFT算法，包括基-2DIFFFT、基-4FFT、分裂基FFT等。6.4節(jié)討論了FFT算法在現(xiàn)實應用中的一些考慮，給出了FFT算法的兩種實用技巧。6.5節(jié)給出了一個FFT算法的應用實例，并給出了完整的MATLAB代碼。6.6節(jié)討論了兩種特殊的頻率分析算法，分別是Goertzel算法和ChirpZ算法。6.7節(jié)以雜談的形式，介紹了高斯與FFT算法的一些歷史典故。最后在6.8節(jié)小結了全章內(nèi)容。6.1DFT的困境與機遇在工程實踐中，DFT面臨的最主要的問題是運算效率很低，特別是當點數(shù)比較大的時候，很難用于實時處理。為此，多年來許多科學家致力于對DFT的快速算法的研究，提出了多種DFT的快速運算方案。本節(jié)主要分析DFT運算所需的運算量及總體的解決思路。6.1.1DFT面臨的主要問題根據(jù)定義，信號x（n）的N點DFT為：相應地，X（k）的N點IDFT為：僅從運算量的角度看，式（6.1）和式（6.2）除了相差一個常數(shù)因子1/N外，DFT正變換和逆變換的運算量完全一樣，而且計算中只包含乘法和加法運算。如果x（n）為復數(shù)信號，則根據(jù)式（6.1）做N點DFT需要N次復數(shù)乘法和N（N－1）次復數(shù)加法。另外我們知道，復數(shù)運算也同樣要通過實數(shù)運算來實現(xiàn)，并且每一個復數(shù)乘法包括4次實數(shù)乘法和2次實數(shù)加法，每一個復數(shù)加法包括2次實數(shù)加法。于是，對于某一個k值，計算X（k）需要4N次實數(shù)乘法和2N＋2（N－1）＝2（2N－1）次實數(shù)加法。一個完整的N點DFT共需要4N次實數(shù)乘法和2N（2N－1）次實數(shù)加法。N點DFT的運算次數(shù)如表6.1所列。表6.1N點DFT的運算量分析當N>>1時，復數(shù)加法次數(shù)（N－1）N≈N，實數(shù)加法次數(shù)2N（2N－1）≈4N。所以無論是從復數(shù)的角度，還是從實數(shù)的角度，都可以說N點DFT的乘法和加法次數(shù)均與N成正比。當N比較小時，運算的問題可能還不突出；但當N較大時，運算量非常大。比如，當N＝4時，僅需要16次的復數(shù)乘法和12次的復數(shù)加法，或者說64次的實數(shù)乘法和56次的實數(shù)加法。而當N＝1024時，則需要1048576次的復數(shù)乘法和加法，或者說4194304次的實數(shù)乘法和加法，所需的運算量非常巨大，實時實現(xiàn)非常困難。在工程實際中，N取1024是一個非常普遍的情況，很多情況下，還要比這個取值大得多。因此，DFT應用于實際最大的問題在于運算量太大，很難滿足實時性的要求。6.1.2DFT高效運算的可能性為了提高DFT的運算效率，使其在工程實踐中發(fā)揮作用，科學家們進行了大量的研究。在長期研究的基礎上，發(fā)現(xiàn)在式（6.1）的計算中實際上存在大量的冗余計算，并且冗余的計算主要來自旋轉因子如下的一些特點：①對稱性：②周期性：③可約性：④特殊值：先通過一個例子來認識如何通過旋轉因子的這些特點來提高DFT的運算效率。假定x（n）為4點長的信號，對其做4點的DFT，則由式（6.1）可寫為：很明顯，如果按照原始的定義式，計算4點的DFT需要16次的復數(shù)乘法和12次的復數(shù)加法。但利用旋轉因子式（6.3）～式（6.6）所示的4個特點，式（6.7）可簡化為：由式（6.8）可以看出，此時計算4點的DFT僅需要1次復數(shù)乘法和8次復數(shù)加法，比按原始定義計算的運算量要明顯減少。從原始定義出發(fā)，DFT的運算量與N成正比，顯然N越小越有利，因此也就很自然地希望能將大點數(shù)的DFT分解為多個小點數(shù)的DFT來計算。而在這種分解過程中，充分利用旋轉因子的對稱性、周期性、可約性及特殊值等特點，為降低DFT運算量提供了可能性。這是各種提高DFT運算效率算法的基本出發(fā)點。由于DFT快速算法的重要性，在數(shù)字信號處理中，專門用快速傅里葉變換（FastFourierTransform，F(xiàn)FT）來表示DFT的快速計算。由此也可知，F(xiàn)FT并不是一種新的變換，而只是DFT的快速算法。FFT算法主要分為兩大類：按時間抽取方法（Decimation-In-Time，DIT）和按頻率抽取方法（Decimation-In-Frequency，DIF）。根據(jù)抽取長度，F(xiàn)FT又分為基-2、基-4等算法，在一個DFT運算中還可同時采用不同基的混合基算法，也稱為分裂基算法。6.2基-2DITFFT算法按時間抽取算法是指在時間上將信號長度逐步減小的算法，基-2算法則是指以2為基底降低信號長度的過程。基-2DITFFT算法通常也稱為庫利—圖基算法，是最早提出的DFT快速算法，是數(shù)字信號處理的發(fā)展歷史上具有里程碑意義的算法。我們先介紹算法的原理，然后再進一步討論算法的運算規(guī)律。6.2.1算法原理基-2DITFFT算法要求信號長度N＝2，其中M為整數(shù)，若信號長度不滿足這個條件，可以通過補零來滿足要求。為了計算長度為N的信號x（n）的DFT，可以在時間上以2為基底，先將x（n）分解為兩個長度均為N/2點的信號x0（n）和x1（n），其中x0（n）為x（n）的偶數(shù)點部分，x1（n）為x（n）的奇數(shù)點部分，用數(shù)學公式表示如下：然后再分別計算x0（n）的N/2點DFTX0（k）和x1（n）的N/2點DFTX1（k），用數(shù)學公式表示如下：最后對X0（x）和X1（k）進行某種組合，可以得到x（n）的DFTX（k）。對照X（k）的定義式：利用旋轉因子的如下特性：式（6.11）可變形為：由式（6.10）和式（6.13），可以得到X（k）與X0（k）及X1（k）有如下關系：對照式（6.10）和式（6.14），細心的讀者可能會問，X0（k）和X1（k）中k的取值范圍是0～（N/2－1），但X（k）中k的取值范圍卻是0～（N－1），在0～（N/2－1）范圍內(nèi)，根據(jù)X0（k）和X1（k），利用式（6.14）計算得到X（k）很好理解，但在N/2～（N－1）范圍時，式（6.14）還能適用嗎？回答同樣是肯定的，因為由DFT的周期性可知，X0（k）和X1（k）的周期均為N/2，所以只要知道0～（N/2－1）范圍內(nèi)X0（k）和X1（k）的取值，就意味著也知道了N/2～（N－1）范圍內(nèi)X0（k）和X1（k）的取值，因而同樣可以用式（6.14）計算X（k）。式（6.14）的運算關系可以用圖6.1所示的流程圖表示，因為其形狀如蝴蝶一般，故也稱為蝶形運算。當支路上沒有標出系數(shù)時，表示該支路的系數(shù)為1。由圖可以看出，完成一個基本的蝶形運算，需要兩次復數(shù)乘法。但利用旋轉因子的如下特性：式（6.14）可改寫為：圖6.1所示的流程圖可以簡化為圖6.2所示的形狀，此時完成一個基本蝶形運算只需要一次復數(shù)乘法，從而進一步提高了運算效率。圖6.1蝶形運算符號圖6.2只需要一次復數(shù)乘法的蝶形運算采用圖6.2所示的形式，可以將前面討論的分解過程表示于圖6.3中。圖中的取值為N＝2＝8。圖6.3將N點DFT分解為兩個N/2點DFT的按時間抽取信號流程圖從直觀上可以理解，采用這種將大點數(shù)DFT分解為小點數(shù)DFT的方法，然后再充分利用旋轉因子的特點，可以提高DFT的運算效率。實際的情況又是怎樣呢？由圖6.3可知，通過一次分解之后，N點的DFT包括如下的運算：兩個N/2點的DFT和N/2個蝶形運算。每個蝶形運算包括一次復數(shù)乘法和兩次復數(shù)加法。分解后的運算量如表6.2所列。一次分解后，總的復數(shù)乘法次數(shù)為N（N＋1）/2≈N/2，總的加法次數(shù)為N/2。而直接進行N點DFT的話，由表6.1可知，需要N次復數(shù)乘法和近似N次復數(shù)加法。這表明，一次分解后，運算量近似為直接計算的1/2。表6.2一次分解后N點DFT的運算量由此可見，上述這種先在時間上按奇偶數(shù)將大點數(shù)DFT分解為小點數(shù)DFT，然后再對小點數(shù)DFT結果進行組合的方法，的確可以提高DFT的運算效率。那么很自然地，如果N/2仍然是偶數(shù)，則對圖6.3中的兩個N/2點的DFT，也可以按照式（6.9）～式（6.14）重復同樣的過程，將每個N/2點的DFT都分解為兩個N/4點的DFT。若用x2（n）和x3（n）分別表示x0（n）的偶數(shù)點部分和奇數(shù)點部分，X2（k）和X3（k）分別表示其N/4點的DFT，則與式（6.14）類似，X0（k）可寫為：用x4（n）和x5（n）分別表示x1（n）的偶數(shù)點部分和奇數(shù)點部分，X4（k）和X5（k）分別表示其N/4點的DFT，則與式（6.17）類似，X1（k）可寫為：考慮到，經(jīng)過第2次分解后，圖6.3所示的流程圖變?yōu)閳D6.4所示的流程圖。圖6.4將N點DFT分解為4個N/4點DFT的按時間抽取信號流程圖若N/4仍然為偶數(shù)，則還可以按照同樣的步驟，分解為兩個N/8點的DFT。這樣的過程一直持續(xù)下去，直到DFT的點數(shù)變?yōu)?點。實際上，2點的DFT仍然可以分解為兩個1點的DFT，因為1點的DFT就是信號本身，兩個1點的DFT，或者說一個2點的DFT，正好構成一個基本的蝶形運算。因此，通常將大點數(shù)DFT按級分解，直到DFT的點數(shù)變?yōu)?。在圖6.4所示的例子中，經(jīng)過兩次分解之后，N/4＝8/4＝2，這2點的DFT構成一個最基本的蝶形運算。圖6.5給出了8點DFT全部分解過程的運算流程圖。圖6.58點基-2DITFFT算法完整流程圖6.2.2算法特點從6.2.1小節(jié)的討論過程可以看出，基-2DITFFT算法的推導雖然看起來復雜，但卻非常有規(guī)律，下面就進一步探討FFT算法的特點。1．級的概念在前面的推導過程中，先將N點的DFT分解為兩個N/2點的DFT，再將兩個N/2點的DFT分解為4個N/4點的DFT，再將4個N/4點的DFT分解為8個N/8點的DFT，直至分解為N/2個2點的DFT。2點的DFT構成一個基本的蝶形運算。每分解一次，稱為一級運算，由于N＝2，所以N點的DFT可以分解為M級蝶形運算，如圖6.5所示。圖中N＝2＝8，總共有3級蝶形運算，依次為m＝1級、m＝2級、m＝3級。每一級蝶形運算都包括N/2個基本的蝶形運算。2．同址運算同址運算也稱為原位運算，指的是在運算的過程中，輸入和輸出占用相同的存儲地址，從而可以節(jié)約存儲空間，進而降低硬件成本。FFT運算是典型的同址運算。觀察圖6.5所示的完整的FFT運算流圖，可以看出每個蝶形運算的輸入和輸出并不交叉。任何一個蝶形運算的兩個輸入經(jīng)過蝶形運算后就失去了利用價值，不再需要保存，因此可以實現(xiàn)原位運算。FFT是由多級蝶形運算構成，若輸入數(shù)據(jù)最初保存在A（0）～A（N－1）的存儲單元中，第一級蝶形運算完成之后，運算結果保存在A（0）～A（N－1）中，原來保存在A（0）～A（N－1）中的原始輸入數(shù)據(jù)被覆蓋掉了。上一級的輸入為下一級的輸出，每一級蝶形運算均為同址運算，最終的DFT結果也保存在A（0）～A（N－1）中，中間各級蝶形運算的結果都沒有被保存。這即是說，存放數(shù)據(jù)的存儲單元始終為N個。3．旋轉因子在圖6.5所示的基-2DITFFT算法流程圖中，每一級都有N/2個蝶形運算，每個蝶形運算都要乘以旋轉因子。每一級旋轉因子都不相同，但排列卻很有規(guī)律，如表6.3所列。表6.3旋轉因子排列規(guī)律表6.3中給出了每級不同旋轉因子的個數(shù)及分組。這里所謂的組，指的是具有相同結構及旋轉因子的一類蝶形運算，每級蝶形運算都包括N/2個蝶形運算，但在不同的級中，分組是不一樣的，如在m＝1級分成了4組，在m＝2級分成了2組，在m＝3級分成了1組。更一般意義上說，第m＝L級，可分為N/2組。4．倒位序規(guī)律由圖6.5可以看出，對于同址運算結構，F(xiàn)FT的輸出X（k）是按正常順序排列在存儲單元中，即按照X（0）、X（1）、X（2）、…、X（N－1）的順序排列。由于要對輸入x（n）逐級進行奇、偶分解，結果導致了輸入信號不是按自然順序存儲的，看起來好像“雜亂無章”，實際上卻是有規(guī)律的，我們稱之為倒位序。所謂倒位序，就是將二進制數(shù)的最高有效位到最低有效位的位序進行顛倒排列而得到的二進制數(shù)。例如，十進制數(shù)6，用二進制可以表示為（110）2，其倒位序為（011）2，即二進制數(shù)3。下面來看看，每一級的奇偶分解是如何導致了輸入的倒位序排列的。如果十進制數(shù)n，用二進制數(shù)表示為（n2n1n0）2。第1次分解時，偶數(shù)在上，奇數(shù)在下，如圖6.3所示。此時，若最低位n0為0，則對應為偶數(shù)；n0為1，則對應為奇數(shù)。對n0為0的部分再進行分解，同樣是偶數(shù)在上，奇數(shù)在下。很明顯，此時若次低位n1為0，則對應為偶數(shù)；n1為1，則對應為奇數(shù)。對n0為1的部分也進行同樣的分解。這一級的分解即為第2次分解，如圖6.4所示。依次類推，繼續(xù)分解相當于逐位判決。最終按位歸類排列就可得到倒位序的排列順序。這種在時間上按照奇數(shù)、偶數(shù)不斷分解的過程，可用圖6.6所示的樹狀圖來描述，從中可以更清楚地理解輸入信號倒位序排列形成的原因。圖6.6描述倒位序的樹狀圖需要說明的是，這里的二進制數(shù)用3位來描述，隱含的前提是N＝2＝8。對于更一般的情況，N＝2，二進制數(shù)要用M位來描述，過程是完全一樣的。5．倒位序的實現(xiàn)一般實際運算中，總是先按自然順序將輸入信號存入存儲單元，然后通過整序或者說變址，得到輸入信號倒位序的排列。如果輸入信號的序號n用二進制數(shù)（n2n1n0）表示，其倒位序二進制數(shù)為（n0n1n2），對應的十進制數(shù)用表示。這樣，在原來自然順序時應該存儲x（n）的單元，倒位序后要存儲。表6.4列出了N＝8時的自然順序二進制數(shù)及相應的倒位序二進制數(shù)。表6.4N＝8時自然順序與倒位序的對應關系把按自然順序排列的輸入信號，轉換成FFT所需的倒位序排列，其過程如圖6.7所示。當時，不必調(diào)換。如n＝0時，對應的也為0，此時不必調(diào)換。當時，必須將原來存放x（n）的單元與原來存放的單元互換。如n＝1時，對應的，原來存放x（1）的單元為A（1），原來存放x（4）的單元為A（4），此時要將A（1）與A（4）兩個單元中的存儲內(nèi)容互換，調(diào)換后A（1）中存放x（4），A（4）中存放x（1）。為了避免把已調(diào)換過的數(shù)據(jù)再次調(diào)換，保證只調(diào)換一次，只需要看是否比n小，若比n小，則意味著此時x（n）在前邊已經(jīng)和互相調(diào)換過了。只有當時，才要調(diào)換。這樣就得到輸入所需的倒位序順序。需要說明的是，圖6.7所示的整序運算也是同址運算。在整序的過程中，數(shù)據(jù)是成對的互換，因此不需要額外的存儲器。圖6.7從正常順序到倒位序的整序過程6.2.3運算量分析從前面算法的推導過程可以看出，F(xiàn)FT算法是一個非常純粹的數(shù)學技巧，幾乎不涉及物理上的意義。因此對FFT算法，運算量的減少才是最為關鍵的。在表6.2中已經(jīng)證實，通過將大點數(shù)DFT分解為小點數(shù)DFT來計算，可以節(jié)省運算量。那么，對N＝2點的FFT來說，經(jīng)過M級分解之后，運算量的情況又如何呢？對于M級蝶形運算，每一級都由N/2個蝶形構成，因此全部N點的FFT共有M·N/2個蝶形運算。每個蝶形運算需要一次復數(shù)乘法和兩次復數(shù)加法。所以，N點的FFT所需的復數(shù)乘法次數(shù)mF和復數(shù)加法次數(shù)aF分別為：而由表6.1可知，N點DFT所需的復數(shù)乘法次數(shù)和復數(shù)加法次數(shù)分別為N和N（N－1）。直接計算DFT與FFT算法復數(shù)乘法的運算量之比為：表6.5給出了N為不同取值情況下的復數(shù)乘法運算次數(shù)及對應的比值。由表可以看出，點數(shù)越大，F(xiàn)FT算法的優(yōu)勢越明顯。直接計算DFT與FFT算法復數(shù)加法的運算量比較與表6.5所列的復數(shù)乘法的運算量比較結果完全類似，此處不再另外給出。表6.5直接計算DFT與FFT算法的復數(shù)乘法運算量比較6.3其他常用的FFT算法在6.2節(jié)中比較詳細地介紹了基-2DITFFT算法。除此之外，其他常用的FFT算法還有基-2DIFFFT、基-4FFT及分裂基FFT等。下面對這3種FFT算法的原理簡要介紹一下。6.3.1基-2DIFFFT如果將FFT算法也看成是一個系統(tǒng)的話，其輸入為x（n），輸出為X（k）。在基-2DIT算法中是將輸入x（n）按偶數(shù)點和奇數(shù)點來分解，類似地，也可以將輸出X（k）按偶數(shù)點和奇數(shù)點來分解。這也即是說，先計算偶數(shù)點的X（k），再計算奇數(shù)點的X（k）。遵循的同樣是將大點數(shù)DFT分解為小點數(shù)DFT這個直觀的原則?；?2DIF同樣要求N＝2。第1步是將X（k）按偶數(shù)點和奇數(shù)點劃分為X0（k）和X1（k），用數(shù)學公式表示如下：第2步是分別計算X0（k）和X1（k），用數(shù)學公式表示如下：式（6.23）和式（6.24）中，n和k的取值范圍均為0，1，2，…，N/2－1。式（6.23）和式（6.24）所示的運算可以用如圖6.8所示的蝶形圖表示，圖中N＝8。與圖6.3所示的按時間抽取相比，此時在計算N/2點DFT之前要先對信號x（n）的前半部分和后半部分進行組合。圖6.8將N點DFT分解為兩個N/2點DFT的按頻率抽取信號流程圖按照式（6.22）～式（6.24）所示的分解方法，可以再將2個N/2點的DFT分解為4個N/4點的DFT，再分解為8個N/8點的DFT，一直分解下去，直到分解為2點的DFT。圖6.9給出了8點DIFFFT算法的完整流程圖。由圖可以看出，輸入x（n）是正常順序的，輸出X（k）是倒位序的。這與圖6.5所示的按時間抽取方法正好相反。圖6.98點基-2DIFFFT算法完整流程圖6.3.2基-4FFT在前面基-2FFT算法的討論過程中，先是將N點DFT分解為2個N/2點DFT，再分解為4個N/4點DFT，再分解為8個N/8點DFT，一直分解下去。每一次的分解后DFT的點數(shù)變?yōu)樵瓉睃c數(shù)的1/2，因而稱為基-2FFT。類似的，基-4FFT算法指的是每一次的分解后點數(shù)變?yōu)樵瓉睃c數(shù)的1/4。最開始是N點的DFT，第1次分解后是4個N/4點的DFT，第2次分解后是16個N/16點的DFT，按照這個規(guī)律一直分解下去。在基-2FFT算法中，要求N＝2，最后一直分解到2點的DFT為最小單位。類似，基-4FFT算法中，要求N＝4，最后一直分解到4點的DFT為最小單位。基-2FFT算法有DIT和DIF兩種形式，基-4FFT算法也同樣有DIT和DIF兩種形式。這里僅以DIT為例來介紹基-4FFT算法。第1步是按如下方式將x（n）分為4段，每段長度為N/4：第2步是分別計算X0（k）、X1（k）、X2（k）、X3（k）：第3步是根據(jù)X0（k）、X1（k）、X2（k）、X3（k）組合得到X（k），推導過程與基-2算法類似，下面僅給出結論：式（6.27）中k的取值范圍為0，1，2，…，N/4－1?；?4FFT算法的基本蝶形運算如圖6.10所示。圖6.10基-4FFT算法基本蝶形運算流程圖經(jīng)過一次分解后的蝶形圖如圖6.11所示。圖中N＝16。圖6.11將N點DFT分解為4個N/4點DFT的按時間抽取信號流程圖按照式（6.25）～式（6.27）所示的分解方法，可以再將4個N/4點的DFT分解為16個N/16點的DFT，以此類推，一直分解下去，直到分解為4點的DFT。圖6.12給出了16點DITFFT算法的完整流程圖。圖中輸入為倒位序，輸出為正常順序，其原因與基-2DITFFT算法的完全一致。圖6.1216點基-4DITFFT算法完整流程圖6.3.3分裂基FFT所謂的分裂基FFT算法實際上是一種綜合使用基-2和基-4的FFT算法。它要求N＝2，基本原理仍然是將大點數(shù)的DFT分解為多個小點數(shù)的DFT，但在具體做法上有所不同。先是如基-2算法一般將x（n）分解為偶數(shù)點部分和奇數(shù)點部分，然后再將奇數(shù)點部分分解為兩部分。也即是說，分裂基算法對x（n）的偶數(shù)點部分是按基-2算法來分解，對奇數(shù)點部分是按基-4算法來分解。分裂基FFT算法也同樣可以按時間抽取或者按頻率抽取。下面以按時間抽取為例來介紹。第1步，對x（n）的偶數(shù)點部分進行基-2分解，奇數(shù)點部分進行基-4分解，用數(shù)學公式表示如下：第2步，分別計算X0（k）、X1（k）、X2（k）：第3步，根據(jù)X0（k）、X1（k）、X2（k）組合得到X（k），推導過程與基-2及基-4算法類似，下面僅給出結論：式（6.30）中k的取值范圍為0，1，2，…，N/4－1。分裂基FFT算法的基本蝶形運算如圖6.13所示。圖6.13分裂基FFT算法基本蝶形運算流程圖經(jīng)過一次分解后的蝶形圖如圖6.14所示，圖中N＝16。由圖可見，一次分解使整個運算由一個N/2點DFT、2個N/4點DFT及N/4個基本的分裂基蝶形運算組成。圖6.14分裂基算法按時間抽取一次分解之后的信號流程圖像基-2和基-4算法一樣，還可以對x0（n）、x1（n）和x2（n）做類似的分解，即將x0（n）分解為一個N/4點DFT和兩個N/8點DFT，將x1（n）分解為一個N/8點DFT和兩個N/16點DFT，將x2（n）也同樣分解為一個N/8點DFT和兩個N/16點DFT。按照這種方法一直分解下去，直到DFT點數(shù)如果是2或4，這樣就可以直接表示成基-2或基-4的基本蝶形運算，此時分解結束。例如對N＝2＝16點的DFT，第1級分解中，16點的x（n）分解為一個8點的x0（n）和4點的x1（n）及x2（n）。第2級分解中，8點的x0（n）可分解為4點的x3（n）和2點的x4（n）及x5（n），而x1（n）及x2（n）均為4點，可直接用基-4的基本蝶形運算實現(xiàn)，因而無需再分解。經(jīng)過第2級分解后，x3（n）為4點，可直接用基-4的基本蝶形運算實現(xiàn)，也無需再分解，x4（n）及x5（n）均為2點，可直接用基-2的基本蝶形運算實現(xiàn)，因而也無需再分解。也即是說，16點的DFT用分裂基FFT算法，只需要兩級分解。在第1級分解中，還包括了4個基本的分裂基蝶形運算；在第2級分解中，包括了2個基本的分裂基蝶形運算。圖6.15給出了N＝16時分裂基FFT算法的結構示意圖。圖6.1516點分裂基FFT算法結構示意圖根據(jù)圖6.15可以很方便地畫出N＝16時的分裂基FFT算法完整流程圖，如圖6.16所示。此時的輸入為倒位序，輸出為正常順序。原因與基-2DITFFT算法的完全一致。圖6.1616點分裂基FFT算法按時間抽取完整流程圖6.3.4算法比較為了更方便比較基-2、基-4及分裂基FFT算法的特點，先畫出N＝16時的基-2DITFFT算法流程圖，如圖6.17所示。圖6.1716點基-2DITFFT算法完整流程圖對照圖6.17、圖6.16和圖6.12，可以發(fā)現(xiàn)，無論是基-2算法、基-4算法還是分裂基算法，蝶形運算圖看起來都差不多，最主要的區(qū)別在于旋轉因子的不同?；?2算法規(guī)律性最強，運算效率也已經(jīng)比較高，那么為什么還要引入基-4算法及分裂基算法呢？在本章的討論中一直強調(diào)，F(xiàn)FT算法更多的是數(shù)學上的運算技巧，而少有物理上的意義。為此，我們主要還是要從運算量的角度來比較基-2、基-4及分裂基算法的優(yōu)缺點。基-2算法的運算量在前面已經(jīng)分析過，下面再分析一下基-4算法及分裂基算法的運算量。1．基-4FFT算法運算量分析由圖6.10所示的基-4算法基本蝶形圖可知，每一個基本的蝶形運算包括3次復數(shù)乘法及8次復數(shù)加法。由圖6.12所示基-4算法的完整流程圖可知，N＝4點的DFT，總共包括M級蝶形運算，每級蝶形運算包括N/4個基本的蝶形運算，并且第1級的蝶形運算沒有乘法，因此所需的復數(shù)乘法次數(shù)mF和復數(shù)加法次數(shù)aF分別為：表6.6中第3列給出了N為不同取值時基-4算法所需的復數(shù)乘法的次數(shù)，符號“-”表示不能進行基-4分解。對比式（6.19）與式（6.31）可以發(fā)現(xiàn)，基-4算法所需的復數(shù)乘法次數(shù)比基-2算法減少約25％?；?2算法和基-4算法所需的復數(shù)加法次數(shù)完全相同。表6.6基-2、基-4及分裂基算法復數(shù)乘法次數(shù)2．分裂基算法運算量分析由圖6.13所示的基本分裂基蝶形圖可知，每一個基本的蝶形運算包括2次復數(shù)乘法和6次復數(shù)加法?？梢宰C明，對于N＝2點的DFT來說，其所需的基本蝶形運算由下式確定：式中BM表示N＝2時所需的基本蝶形運算數(shù)目，于是，分裂基算法所需的復數(shù)乘法次數(shù)mF為：所需的復數(shù)加法次數(shù)aF仍然是Nlog2N。詳細的推導過程不再給出，感興趣的讀者可參考相關文獻。表6.6中第4列給出了N為不同取值時分裂算法所需的運算量。對比式（6.19）、式（6.31）與式（6.34）可以發(fā)現(xiàn)，分裂基算法所需的復數(shù)乘法次數(shù)比基-2算法減少約33％，比基-4算法減少約11％。分裂接算法與基-2算法和基-4算法所需的復數(shù)加法次數(shù)完全相同。對比圖6.2所示的基-2基本蝶形運算和圖6.10所示的基-4基本蝶形運算，每個基-4的基本蝶形運算由兩級基-2基本蝶形運算構成，每級都包括2個基本的基-2蝶形運算，不同的是，基-4算法更合理地利用了旋轉因子的特性，從而減少1次復數(shù)乘法。因此，在滿足基-4算法運算條件的前提下，基-4算法效率比基-2算法更高。分裂基算法充分利用了基-2算法中每一級的每一組蝶形運算中上半部分均沒有旋轉因子，而只在下半部分有旋轉因子的特點，重新優(yōu)化了旋轉因子的排列結構，同時還用上了基-4基本蝶形運算效率高于基-2基本蝶形運算的特點，從而進一步提高了運算效率。在已知的N＝2的各種算法中，分裂基算法的運算效率最高?；?2、基-4和分裂基算法都是常用的FFT算法?；?4算法雖然運算效率高于基-2算法，但其所要求的N＝4的條件比基-2算法要苛刻，在不滿足條件的情況下，補零的個數(shù)過多，反過來降低了運算效率。分裂基算法的運算效率也高于基-2算法，但編程的復雜度也比較高。因此，在實際應用中最常用的還是基-2算法。最后還要說明的是，基-2算法、基-4算法和分裂基算法都是同址運算。無論是基-2、基-4還是分裂基算法，按時間抽取與按頻率抽取所需的運算量完全相同。6.4現(xiàn)實應用中的考慮前面主要以數(shù)學公式和流程圖的方式介紹了基-2、基-4及分裂基FFT算法的基本原理。流程圖的方法能非常清晰地描述信號的變化情況，是理解FFT算法最有效的方式，但在具體實現(xiàn)時還有些問題需要考慮，比如倒位序的整序和旋轉因子的生成等。倒位序的問題已經(jīng)在6.2.2小節(jié)討論過，這里再介紹一下旋轉因子的生成。另外，在實際的應用過程中，充分利用FFT算法的特點，也可進一步提高整個的處理效率。這里主要介紹了IDFT的快速計算和實信號的FFT兩個技巧。6.4.1旋轉因子的生成在實際的FFT運算中，通過計算＝cos（2πn/N）－jsin（2πn/N）來生成旋轉因子。由于計算正弦和余弦函數(shù)值的運算量很大，程序運行時，旋轉因子的產(chǎn)生將對運算速度產(chǎn)生較大的影響。在對時間要求比較苛刻時，應該對此進行改進。一種有效的方法是查表法。查表法的思想是在進行FFT運算前，將所有用到的旋轉因子，n＝0，1，2，…，N/2－1計算出來，并存儲在計算機中，作為旋轉因子表。在程序執(zhí)行時不必再進行計算，而是直接查表得到所需的旋轉因子值。這樣可使運算速度大大提高，當然查表法的缺點是需要占用較多的計算機內(nèi)存。6.4.2IDFT的快速計算前面討論的FFT算法可以快速地根據(jù)x（n）計算得到X（k）。在實際應用中，往往還需要根據(jù)X（k）計算得到x（n），這時同樣可以用FFT算法快速得到。IDFT的數(shù)學公式為：觀察式（6.35），為了要利用FFT算法，首先是式中的旋轉因子應變?yōu)?。因為為復?shù)，這步的變化可以通過求共軛來實現(xiàn)。因此，先對式（6.35）兩邊求共軛，得到：式（6.36）可以用FFT算法實現(xiàn)。但此時得到的并不是x（n），而是x（n）的共軛。因此，還需要對得到的x（n）再取一次共軛即可得到x（n）。用數(shù)學公式表示如下：式（6.37）表明，為了利用FFT算法由X（k）計算得到x（n），第1步是對X（k）求共軛；第2步是對第1步得到的結果進行FFT運算；第3步是對第2步得到的結果求共軛；第4步是對第3步得到的結果除以點數(shù)N。這個過程的框圖如圖6.18所示。圖6.18用FFT算法快速計算IDFT6.4.3實信號的FFT前面討論的FFT算法，由于旋轉因子為復數(shù)，在第1級蝶形運算之后，無論輸入x（n）是復信號還是實信號，所有的運算均為復數(shù)運算。因此FFT算法都默認輸入x（n）是一個復數(shù)。但現(xiàn)實世界中的信號通常都為實信號，按照前面的方法計算實信號的FFT，最直觀的辦法是把x（n）看作是一個虛部為零的復信號。顯然，這樣x（n）虛部的信息沒有得到充分的利用，從而降低了運算效率。因此，對于實信號的FFT，還可以考慮用另外的方法來進一步提高效率。下面分兩種情況來討論。1．兩個實信號的FFT若輸入x1（n）與x2（n）均為N點的實信號，其DFT分別為X1（k）和X2（k）。將x1（n）看作是復數(shù)x（n）的實部，x2（n）看作是復數(shù)x（n）的虛部，于是復數(shù)x（n）可寫為：對于N點的復信號x（n），其DFTX（k）可用FFT快速計算，此時沒有冗余計算。再根據(jù)X（k），可以得到所要計算的X1（k）和X2（k）。這就是計算兩個實信號FFT的基本思路。要注意的是，由于X1（k）和X2（k）也均為復數(shù)，因而不能由X（k）的實部和虛部簡單求得。為了根據(jù)X（k）求得X1（k）和X2（k），最好是先用x（n）來完全表示x1（n）和x2（n）。由式（6.38），x1（n）和x2（n）可用x（n）完全表示如下：于是X1（k）和X2（k）可分別寫為：由DFT的復共軛特性可知，DFT［（x（n）］＝X（N－k），于是可將上式改寫為：利用式（6.41），可根據(jù)X（k）得到X1（k）和X2（k）。簡單再回顧一下用FFT快速計算兩個實信號DFT的過程。第1步是將兩個實信號x1（n）和x2（n）組合成一個復數(shù)x（n），其中x1（n）為實部，x2（n）為虛部；第2步是用FFT算法快速計算x（n）的DFTX（k）；第3步是根據(jù)X（k）利用式（6.41）分別計算X1（k）和X2（k）。從直觀上就能理解用上述方法能提高運算效率，但具體數(shù)值是多少呢？直接將虛部補零的方法計算兩個實信號的FFT所需的復數(shù)乘法mF1和復數(shù)加法aF1分別為：用組合成復數(shù)的方法計算所需的復數(shù)乘法mF2和復數(shù)加法aF2分別為：以乘法為例，提高的運算效率為：圖6.19給出了N為不同取值時運算效率的提高曲線。由圖可以看出，隨著N的增加，用組合成復數(shù)的方法可以提高約45％的運算效率。圖6.19用組合成復數(shù)的方法快速計算兩個實信號DFT的運算量提高曲線2．一個2N點的實信號的FFT前面討論了兩個N點實信號DFT的快速計算方法，其中最基本的思想是將兩個實信號組合成一個復信號，這樣就充分利用了虛部的信息，從而降低了約一半的運算量。若要求一個2N點實信號的FFT，情況又如何呢？借鑒前面的方法，一個2N的實信號同樣也可以組合成一個復信號，然后充分利用虛部的信息，從而降低運算量。具體做法是，第1步，將2N點的實信號x（n）分解為兩個N點的實信號x0（n）和x1（n），其中x0（n）為x（n）的偶數(shù)點部分，x1（n）為x（n）的奇數(shù)點部分，用數(shù)學公式表示如下：第2步，將兩個N點的實信號x0（n）和x1（n）組合成一個N點的復信號y（n），其中x0（n）為實部，x1（n）為虛部，用數(shù)學公式表示如下：第3步，用FFT算法計算y（n）的DFTY（k），并根據(jù)式（6.41）的結果，有：第4步，根據(jù)X0（k）和X1（k）的結果，組合計算可得到X（k）。具體的方法與基-2DITFFT算法第1級的分解完全一致，其運算流程如圖6.20所示。用數(shù)學公式，X（k）可寫為：圖6.20用FFT算法快速計算一個2N點實信號的DFT因為充分利用了虛部的信息，從直觀上能理解用上述方法能提高運算效率，但具體數(shù)值是多少呢？直接將虛部補零的方法計算實信號x（n）的FFT所需的復數(shù)乘法mF1和復數(shù)加法aF1分別為：用組合成復數(shù)的方法計算所需的復數(shù)乘法mF2和復數(shù)加法aF2分別為：以乘法為例，提高的運算效率為：式（6.55）的結果與式（6.46）的結果非常類似，這是因為兩者都是將實信號的FFT轉換為復信號的FFT，內(nèi)核完全一致。6.5應用實例及其MATLAB實現(xiàn)FFT算法是最經(jīng)典的數(shù)字信號處理算法之一，其應用非常廣泛。在第5章的應用實例中，使用的實際上就是FFT算法。從概念上來說，DFT算法更為重要，在具體實現(xiàn)時，往往都是用FFT算法。從這個意義上說，F(xiàn)FT算法的應用即是DFT算法的應用。這里我們要舉例的是FFT算法的另外一種應用。學完FFT算法之后，估計很多朋友會有困惑：FFT算法非?；A，幾乎各種主流的計算機編程語言都有現(xiàn)成的FFT算法軟件包，在工程實現(xiàn)時，只需調(diào)用軟件包即可，為什么還要對FFT算法的原理啰嗦半天呢？我們知道，在數(shù)字信號處理的工程實現(xiàn)中，往往用DSP芯片來實現(xiàn)具體的算法，DSP芯片的價格和片內(nèi)存儲器關系很大。在通信等消費電子產(chǎn)品中，由于價格因素非常敏感，在選用DSP的時候，往往不能選擇片內(nèi)存儲器很大的DSP。對每一種型號的DSP，硬件開發(fā)商都會給出優(yōu)化后的FFT算法源代碼，一般情況下直接調(diào)用就可以。但是，在某個具體的應用中，已經(jīng)根據(jù)主要的任務需求，確定了某款DSP芯片，但這款DSP的存儲容量最大只能做8192點的FFT，而在應用中卻要求能做16384點的FFT，這時怎么辦呢？如果對FFT算法原理不了解的話，這時肯定是無法下手的。如果對FFT算法原理非常熟悉的話，很自然地會考慮將16384點的數(shù)據(jù)先分為2個8192點，先對每個8192點的數(shù)據(jù)直接進行FFT，再分別對這兩次的結果進行組合，即可得到16384點的FFT結果，這是按時間抽取的思路。或者先對輸入數(shù)據(jù)進行某種組合，再將組合后的16384點數(shù)據(jù)分為兩個8192點，分別對其進行FFT，這是按頻率抽取的思路。在本例中，用按頻率抽取的思路更為方便，其流程圖如圖6.21所示。組合后的信號分別為：圖6.21大點數(shù)FFT流程圖完成兩個8192點FFT之后，X0（k）對應著X（k）的偶數(shù)部分，X1（k）對應著X（k）的奇數(shù)部分。這樣就得到了16384點的FFT結果。實際上，上述這種因為芯片內(nèi)存不足而導致FFT需要分解的問題，可統(tǒng)稱為大點數(shù)FFT的問題。雖然半導體器件發(fā)展很快，內(nèi)存不斷增大，但人們的欲望也同樣無止境，因而在工程中總有可能面臨大點數(shù)FFT的問題。按照圖6.21所示的流程圖進行仿真。仿真中，輸入信號為16384點的單頻復正弦信號，頻率為ω0＝0.25π，噪聲為高斯白噪聲。圖6.22給出了組合后的輸入信號x0（n）和x1（n），圖6.22（a）為x0（n）的實部，圖6.22（b）為x1（n）的實部。由于受到噪聲污染，x0（n）和x1（n）看起來均雜亂無章。相對而言，x1（n）還較有規(guī)律，雖然受到噪聲污染，但還是能明顯看出關于4096點的對稱關系。圖6.22組合后的輸入信號分別對x0（n）和x1（n）做8192點FFT，并以X0（k）為X（k）的偶數(shù)部分，X1（k）為X（k）的奇數(shù)部分，這樣的結果稱為分解得到的16384點FFT結果。仿真中，可直接對輸入信號做16384點的FFT，得到的結果稱為直接結果。將直接結果與分解結果相減，得到計算誤差，其幅度如圖6.23所示。由圖可以看出，誤差的幅度非常之小，在10量級，可忽略不計。這也即是說，直接調(diào)用16384點的FFT與按照圖6.21分解計算16384點的FFT結果完全一致。但圖6.21所示的分解的方法，雖然在效率上比直接調(diào)用的方法要低一些，但對DSP芯片的內(nèi)存要求大為降低。圖6.23直接計算與分解計算的誤差完成上述大點數(shù)FFT功能的完整MATLAB代碼如下：6.6兩類特殊的頻率分析算法前面介紹的FFT算法，比直接計算DFT可以節(jié)省大量的運算量。但在有些應用場合中，只需要計算某個或者某幾個頻率點的X（k）值，而不需要計算所有的X（k）值，此時用FFT算法可能并不比直接計算DFT運算量更小。另外，DFT算法是傅里葉變換在單位圓上等間隔采樣的結果，但在某些應用場合，只需要分析信號在某一個很小頻段的頻譜值，而且要對這很小頻段范圍的頻譜值進行比較精細的分析，此時也無法直接利用FFT算法。對于這些特殊的情況，需要有專門的快速算法。6.6.1Goertzel算法對于長度為N＝1024的信號x（n），若只需要計算X（128）這一個頻率點的頻譜值，用FFT算法，要先計算全部的1024點X（k），然后只取X（128）這一個值。此時所需的復數(shù)乘法次數(shù)為5120，所需的復數(shù)加法次數(shù)為10240。也就是說，計算全部的1024點X（k）和只計算X（128）這一個頻點的值，所需的運算量是一樣的。很明顯，這時用FFT的方法浪費的運算量太多。如果只需要一個頻點的頻譜值，直接計算的話也僅需要1024次的復數(shù)乘法和1023次的復數(shù)加法。因此，這時用FFT的方法并不合適。此時，用Goertzel算法可以有更高的運算效率。1．算法原理在第5章中已經(jīng)講過，一個N點的DFT可以等效為N個窄帶的濾波器組。X（k）為x（n）通過第k個窄帶濾波器hk（n）在n＝N時刻的輸出，即：式中，hk（n）為：對hk（n）進行Z變換：于是，為了計算X（k）＝y(tǒng)k（N），還可以根據(jù)式（6.58）所定義的差分方程來進行遞推運算，用數(shù)學公式表示如下：遞推的初始條件為yk（-1）＝0。先根據(jù)yk（-1）得到y(tǒng)k（0），再遞推得到y(tǒng)k（1），再得到y(tǒng)k（2），最后一直遞推得到y(tǒng)k（N）。此即為所要計算的X（k）。由式（6.59）可知，每一次遞推需要一次復數(shù)乘法和一次復數(shù)加法。計算一個X（k）需要遞推N次，也即總共需要N次復數(shù)乘法和N次復數(shù)加法。這比直接按公式計算X（k）多一次復數(shù)加法，但此時只需要計算或者存儲，相比計算或者存儲要簡單得多。對式（6.58）所示的系統(tǒng)函數(shù)進行簡單的變形，還可以進一步提高運算效率。式（6.58）所示的濾波器在處有一個極點，因而可以看作是一個諧振器。對其增加一個共軛極點，式（6.58）可變?yōu)椋菏剑?.60）所示系統(tǒng)的流程圖如圖6.24所示，該系統(tǒng)寫成遞推的差分方程式如下：圖6.24Goertzel算法流程圖遞推的初始條件為vk（-2）＝vk（-1）＝0。為了計算yk（N），只需要通過式（6.61）遞推得到vk（N）和vk（N－1）即可。在式（6.61）中，系數(shù)均為實數(shù)，因而可以節(jié)約運算量。若輸入x（n）為復數(shù)，則由初始條件根據(jù)式（6.61）得到vk（N）和vk（N－1）需要N次遞推，每次遞推需要2次的實數(shù)乘法和2次復數(shù)加法，也即是總共需要2N次實數(shù)乘法和2N次復數(shù)加法。由vk（N）和vk（N－1），根據(jù)式（6.62）只需要遞推一次即可得到y(tǒng)k（N），此時需要1次復數(shù)乘法和1次復數(shù)加法。也即是說，為了計算某個X（k），式（6.61）需要遞推N次，式（6.62）只需要遞推1次。所需的總運算量為2N次實數(shù)乘法，1次復數(shù)乘法，2N＋1次復數(shù)加法。都換算為實數(shù)運算的情況，總共需要2N＋4次實數(shù)乘法，4N＋4次實數(shù)加法。換算成實數(shù)運算，式（6.59）的遞推需要4N次的實數(shù)乘法和4N次的實數(shù)加法。而直接計算DFT需要4N次的實數(shù)乘法和4N－2次的實數(shù)加法。由此可看出，式（6.60）所示系統(tǒng)的運算量要優(yōu)于式（6.58）所示的運算量，也優(yōu)于直接按公式計算DFT。而且，式（6.60）所示的系統(tǒng)同樣保留了旋轉因子運算或者存儲簡單的優(yōu)點。利用式（6.60）所示的系統(tǒng)計算X（k）的方法稱為Goertzel算法。如果計算很少的幾個頻率點的頻譜值，Goertzel算法的效率高于FFT算法，也高于直接計算DFT的算法。2．應用實例電話機是日常生活中的常用物品，其按鍵的排列如圖6.25所示。打電話之前要先進行撥號，其工作的實質(zhì)是每個按鍵產(chǎn)生兩個單頻信號，一個稱為行頻，一個稱為列頻。每個按鍵對應的行頻和列頻也都在圖6.25中標出。一串電話號碼，在頻譜上就對應了一系列的行頻和列頻。在交換機這一端，通過檢測相應的頻率，就可以恢復出所撥打的電話號碼，從而實現(xiàn)通話。在頻率檢測的過程中就要用到Goertzel算法。圖6.25電話機按鍵與頻率的對應關系由于每個按鍵對應了兩個頻率，每個電話包含多個數(shù)字，因此圖6.25所示電話按鍵所產(chǎn)生的信號也稱為雙音多頻（DualToneMultipleFrequency，DTMF）信號。國際電信聯(lián)盟對DTMF信號的要求是，每個按鍵的持續(xù)時間是100ms，其中真正代表按鍵的音頻信號持續(xù)時間至少為45ms，最長為55ms，100ms中的其他時間是靜音，以便將連續(xù)的按鍵區(qū)別開來。由圖6.25可以看出，典型的DTMF信號頻率范圍為700～1700Hz，采樣頻率fs選擇為8kHz，滿足采樣頻率要求。若采樣點數(shù)為N＝205，則根據(jù)k＝N×f/fs，可計算出行頻及列頻所應的k值，如表6.7所列。表6.7DTMF頻率及其對應的k值因此，對于每個DTMF信號，只需要計算出上述8個頻點的頻譜值，即可判斷到底是哪個按鍵被按下，其流程如圖6.26所示。比如說，電話機上按鍵“7”被按下之后，DTMF信號包含兩個頻率的信號，分別為行頻852Hz和列頻1209Hz，此時，僅有X（22）和X（31）的幅度值比較大，X（18）、X（20）、X（24）、X（34）、X（38）及X（42）的幅度值都很小。在實際中可通過設定門限來進行判決。對于不同的DTMF信號，重復相同的過程即可。圖6.26DTMF信號檢測流程圖6.6.2ChirpZ算法假定信號x（n）為窄帶信號，采樣頻率為10kHz，信號的頻率范圍為2～2.5kHz。若采樣點數(shù)N為1024，則用FFT算法計算X（k），有用的頻譜為X（205）～X（256）。也即是說1024點的DFT，僅有52點是有用的，其余的972點的結果都沒有用。對于帶寬范圍內(nèi)的頻譜，若想要進行更精細的分析，最直接的方法是增加頻域的采樣點數(shù)，比如說做2048點的FFT，此時因為帶寬很小，這樣帶來的冗余計算更多。此時，用ChirpZ算法可以更好地分析感興趣帶寬內(nèi)的頻譜。更一般地說，DFT可以看作是對Z變換結果X（z）在單位圓上的等間隔采樣。在信號處理中，除了計算整個單位圓上的X（z）之外，也常常需要計算X（z）在z平面上其他曲線上的取值，這時的快速計算都要用到ChirpZ算法。1．算法原理假定我們希望在一組點zk上計算x（n）的Z變換的值，用數(shù)學公式可表示為：zk為z平面上一條弧線上的采樣點。該弧線從某一點開始，然后向原點盤旋或者離開原點向外盤旋，用數(shù)學公式表示為：圖6.27列出了z0、R0及φ0為不同取值情況下所對應的z平面上的弧線。若R0＞1，則弧線離開原點向外盤旋；若R0＜1，則弧線由外向原點盤旋；若R0＝1，則弧線為一半徑為r0的圓弧。更進一步，若r0＝1、θ0＝0、φ0為2π/N，并且L＝N，則對應的弧線為單位圓。若r0＝1、θ0≠0，則r0應的弧線為單位圓上的一段圓弧，這即是前面所討論的窄帶信號的情況，此時可以通過減小采樣間隔φ0來對信號進行更精細的分析。圖6.27z平面上不同的弧線將式（6.64）中的點zk代入式（6.63）所示的Z變換表達式，可以得到：式中，。如何快速計算式（6.65），正是ChirpZ算法所要解決的問題。如果將式（6.65）中的kn寫為：則式（6.65）可變形為：式中：式（6.67）中的求和可以看作是信號g（n）通過單位沖激響應為h（n）的濾波器的輸出，其中，h（n）為：于是式（6.67）可寫為：式中，y（k）表示信號g（k）通過系統(tǒng)h（k）的輸出，用數(shù)學公式表示為：這樣，只要計算出y（k），就很容易得到X（zk）。而用卷積表示的信號與系統(tǒng)的相互作用，可以利用FFT來快速實現(xiàn)。圖6.28給出了用卷積計算X（zk）的實現(xiàn)框圖。圖6.28ChirpZ算法的運算流程當R0＝1時，為復指數(shù)信號，其頻率為ω＝φ0n/2。很明顯，此信號的頻率隨時間的變化而線性變化，這樣的信號在信號處理中稱為線性調(diào)頻信號，被廣泛應用于雷達、聲納等系統(tǒng)中。因此，利用式（6.67）計算的Z變換稱為ChirpZ算法。由上述的算法原理的分析可知，ChirpZ算法最關鍵的是將Z變換分解為一個信號通過一個系統(tǒng)，然后再利用FFT來快速計算卷積，從而實現(xiàn)了Z變換的快速計算。2．實現(xiàn)步驟由式（6.71）可知，信號g（n）是有限時寬的，其取值范圍為［0，N－1］，而線性系統(tǒng)h（n）則是無限時寬的。幸運的是，按式（6.71）計算y（k）在［0，L－1］范圍內(nèi)的取值只需要h（n）的一部分，即只要用到［-（N－1），（L－1）］范圍內(nèi)的h（n）。此時h（n）也可以看作是有限時寬的，長度為N＋L－1，如圖6.29（b）所示。從卷積的觀點看，g（n）的點數(shù)為N，h（n）的點數(shù)為N＋L－1，則g（n）*h（n）的點數(shù)為2N＋L－2，因而用圓周卷積代替線性卷積且不產(chǎn)生混疊失真的條件是圓周卷積的點數(shù)應≥2N＋L－2，但是，由于只需要前L個值正確，對以后的其他值是否有混疊失真并不感興趣，這樣可將圓周卷積的點數(shù)縮減到最小為N＋L－1。當然，為了進行基-2FFT運算，圓周卷積的點數(shù)應取為M≥N＋L－1，同時又滿足M＝2的最小M。這樣可先將h（n）補零，補到點數(shù)等于M，也就是從n＝L開始，補M－（N＋L－1）個零點，然后將補零后的h（n）以M為周期進行周期延拓，再取主值區(qū)間的值，從而得到進行圓周卷積的一個信號，如圖6.29（c）所示。進行圓周卷積的另外一個信號只需將g（n）補零，使之點數(shù)為M即可，如圖6.29（a）所示。圖6.29ChirpZ算法實現(xiàn)時各階段的波形示意圖由此，可列出ChirpZ算法實現(xiàn)的具體步驟如下：①選擇一個最小的整數(shù)M，使其滿足M≥N＋L－1，同時滿足M＝2，以便采用基-2FFT算法。②將式（6.68）所示的N點信號g（k）補零，變?yōu)镸點的信號，用數(shù)學公式可表示為：其波形示意圖如圖6.29（a）所示。對式（6.72）的信號做M點的DFT，得到③形成M點的信號h（k）。如前所述，在k＝0到L－1一段取h（k）＝V，在k＝L到M－N段取h（k）為0，在k＝M－N＋1到M－1段取h（k）為V的周期延拓V，用數(shù)學公式表示如下：其波形示意圖如圖6.29（c）所示。對式（6.74）的信號做M點的DFT，得到④將式（6.73）所示的G（m）和式（6.75）所示的H（m）相乘，得到Y（m）＝G（m）·H（m）。顯然，Y（m）為M點的離散信號。⑤用FFT法求Y（m）的M點IDFT，得到h（m）與g（n）的圓周卷積，用數(shù)學公式表示如下：其中y（k）的前L個值等于h（k）與g（k）的線性卷積，其波形示意圖如圖6.29（d）所示。y（k）中k≥L的值沒有意義，直接舍棄即可。⑥最后求X（zk），有：3．運算量分析ChirpZ算法所需的復數(shù)乘法如下：①形成M點信號需要N次復數(shù)乘法。系數(shù)可事先計算好并預存在存儲器內(nèi)。②計算g（k）的M點FFT，需要（M/2）log2M次復數(shù)乘法。③計算h（k）的M點FFT，需要（M/2）log2M次復數(shù)乘法。④G（m）與H（m）相乘需要M次復數(shù)乘法。⑤計算Y（m）的M點IFFT，需要（M/2）log2M次復數(shù)乘法。⑥最后形成L點的輸出X（zk），需要L次復數(shù)乘法。綜上所述，ChirpZ算法所需的總的復數(shù)乘法次數(shù)為：而根據(jù)式（6.63）直接計算X（zk）則需要N·L次復數(shù)乘法。當N及L都比較大時（例如，N、L均大于50），ChirpZ算法的運算效率要高得多。由以上討論可以看出，ChirpZ算法非常靈活，它的輸入信號點數(shù)N和輸出點數(shù)L可以不相等，且N和L均可為任意數(shù)；各zk點間的角度間隔φ0可以是任意的，因而頻率分辨率可以調(diào)整；計算Z變換的弧線可以不是圓而是螺旋線；起始點z0可任意選定，也就是說可從任意頻率開始對輸入數(shù)據(jù)進行分析，非常便于做窄帶信號的高分辨率分析。4．應用實例在5.8.1小節(jié)用DFT對雷達回波頻譜進行分析的實例中看到，加海明窗后，v0＝300m/s和v1＝306m/s這兩個目標很難分辨出來。這時自然希望對目標出現(xiàn)的地方再進行更精細的分析。增加FFT點數(shù)是一種解決方法，但這種方法帶來的冗余數(shù)據(jù)太多。另外一種效率更高的方法是ChirpZ算法。在第5章的分析中可知，v0＝300m/s對應的譜線出現(xiàn)在k0＝80，即模擬頻率f0＝20000Hz，v2＝360m/s對應的譜線出現(xiàn)在k2＝96，即模擬頻率f2＝24000Hz。也即是說，為了更詳細地分析目標速度信息，只需要對f0～f2頻率范圍內(nèi)的頻譜進行更詳細地分析即可。留一點富余量，選擇起始的k0為70，結束的k2為99，總共是30根譜線范圍做ChirpZ變換，如圖6.30所示。圖6.30雷達測速中ChirpZ算法在z平面上的弧線ChirpZ算法中，這段頻率范圍的采樣點數(shù)變?yōu)?00點，其結果如圖6.31所示。由圖可以看出，經(jīng)過ChirpZ變換之后，v0＝300m/s和v1＝306m/s這兩個目標也能非常清楚地分辨出來。由圖6.31可知，v0出現(xiàn)在k0＝100，對應的測量速度為：與仿真設定的值完全一致。v1出現(xiàn)在k1＝116，對應的測量速度為306m/s；v2出現(xiàn)在k2＝260，對應的測量速度為360m/s。v1和v2的測量值與仿真設定的數(shù)值完全一致。當然，這個結果是非常理想化的，因為經(jīng)過ChirpZ算法增加f0～f2頻率范圍內(nèi)的采樣點數(shù)之后，三個目標都正好在整數(shù)值的譜線上，測量結果才非常好。一般的情況下，會有一定的誤差。圖6.31雷達測速中ChirpZ算法的結果6.7信號處理雜談：高斯的遺憾快速傅里葉變換（FFT）是一種實現(xiàn)高效DFT運算的快速算法。以1965年庫利（J.W.Cooley）和圖基（J.W.Tukey）發(fā)表在《MathematicsofComputation》雜志上的題為《AnAlgorithmfortheMachineCalculationofComplexFourierSeries》的論文為標志，F(xiàn)FT算法成為數(shù)字信號處理這門學科興起的一座極為重