《數字信號處理》課件-第4章

4.1引言

4.2提高DFT運算效率的基本途徑4.3基2時分FFT算法

4.4基2頻分FFT算法

4.5IDFT的快速算法

4.6實序列DFT的有效計算方法

4.7線性調頻Z(Chirp-Z)變換算法習題與上機題

離散傅里葉變換(DFT)是信號分析與處理中的一種重要變換。因為直接計算DFT的計算量與變換區(qū)間長度N的平方成正比，當N較大時，計算量太大，所以在快速傅里葉變換(FastFourierTransform，F(xiàn)FT)出現(xiàn)以前，直接用DFT進行譜分析和信號的實時處理是不切實際的。4.1引言自從1965年圖基(J.W.Tukey)和庫利(J.W.Cooley)在《計算數學》(MathematicsComputation)雜志上發(fā)表了著名的《機器計算傅里葉級數的一種算法》論文后，桑德-圖基等快速算法相繼出現(xiàn)，很快形成了一套高效運算方法，這種算法使DFT的運算效率提高了幾個數量級，為數字信號處理技術應用于各種信號的實時處理創(chuàng)造了良好的條件，大大推動了數字信號處理技術的發(fā)展。多年來，人們繼續(xù)尋求更快、更靈活的好算法，目前除了基2、基4算法外，還有分裂基、混合基等算法。對于N=1024，直接計算需要復數乘法1048576次，而幾種常用FFT算法使計算量大大降低，例如基2算法需要復數乘法5120次，基4算法需要復數乘法3840次，

分裂基算法需要復數乘法3413次。本章重點介紹基2快速傅里葉變換算法。設x(n)為N點有限長序列，其DFT為

(4.2.1)

其反傅里葉變換(IDFT)為

(4.2.2)4.2提高DFT運算效率的基本途徑正變換與反變換的差別只在于的指數符號不同，以及差一個乘數因子，因此我們以下只討論DFT的運算量，IDFT的運算量與DFT的基本相同。

一般來說，X(k)為復數(只有x(n)為實數且條件偶對稱時，X(k)才為實數)，x(n)多數情況下也為復數(例如：信號經A/D采樣，數字下變頻芯片之后輸出為I、Q正交兩路，詳見第9章)。因此，計算一點X(k)需要N次復數乘法和N-1次復數加法。由于有N點X(k)，因此完成N點序列的DFT共需要N2次復數乘法，即

CM=N2

(4.2.3)

以及N(N-1)次復數加法，即

CA=N(N-1)≈N2

(4.2.4)實際實現(xiàn)中復數乘法和加法都要依賴實數乘法和加法。由于一次復數乘法需要4次實數乘法和兩次實數加法，一次復數加法需要2次實數加法，因此直接計算DFT的實數乘法次數為

RM=4N2

(4.2.5)

實數加法次數為

RA=2N(N-1)+2N2≈4N2

(4.2.6)所以，直接計算DFT時的運算量與N2成正比。設N=1024，則直接計算DFT需要100多萬次(1048576次)復數乘法。在實時計算較寬帶寬的頻譜時，即使使用目前運算速度最快的DSP芯片也無法滿足計算要求。

觀察式(4.2.1)，要想提高DFT的運算效率，只有利用DFT定義式中旋轉因子的特點。根據第3章所講內容，我們可知有以下幾個特點:

(1)共軛對稱性：。

(2)周期性：。

(3)可約性：。利用共軛對稱性，可以合并式(4.2.1)中的乘法項，從而減少了近一半的乘法次數。

利用旋轉因子的周期性和可約性，可以將長序列的DFT分解成短序列的DFT。由于DFT的運算量與N2成正比，N越小，計算量越小。因此將長序列DFT變?yōu)槎绦蛄蠨FT是減少DFT運算量的根本途徑?？焖俑道锶~變換算法正是基于此思路而發(fā)展起來的。

將長序列的DFT劃分為短序列的DFT有兩種劃分方法：一種是將時域序列按奇偶劃分，稱之為時間抽取(DecimationinTime，DIT)算法或基2時分算法;一種是將頻域DFT按奇偶劃分，稱之為頻率抽取(DecimationinFrequency，DIF)算法或基2頻分算法。4.3.1基2時分蝶式運算定理

設序列x(n)的長度為N，且滿足N=2M，M為整數，其N點DFT為X(k)=DFT［x(n)］，0≤k，n≤N-1。按n的奇偶將x(n)分解為兩個點的子序列

x1(r)=x(2r)，r=0，1，…，(4.3.1)

x2(r)=x(2r+1)，r=0，1，…，(4.3.2)4.3基2時分FFT算法若X1(k)=DFT［x1(r)］，X2(k)=DFT［x2(r)］，0≤k≤

-1，則有

(4.3.3a)

(4.3.3b)

證明x(n)的N點DFT為

利用旋轉因子的可約性，即，

X(k)可以寫成

由于X1(k)和X2(k)都隱含周期性，周期為，因此上式可以寫為

當0≤k≤

-1時

當≤k≤N-1時

令，則有

由于，因此

最后

由證明過程可得，N點DFT可以分解為兩個點的DFT，它們按式(4.3.3)又組合成一個N點的DFT。式(4.3.3)也稱為蝶式運算公式。由蝶式運算定理可以看出，要完成一個蝶形運算(即計算一次蝶式運算定理)，需要一次復數乘法和兩次復數加法，計算N點DFT共需計算兩個點的DFT和個蝶形運算。計算點DFT需要次復數乘法和次復數加法。當N1時，計算N點DFT共需由此可見，僅經過一次分解，就使運算量減少近一半。如果N=2M，則我們可以繼續(xù)分解，直至分解到1點DFT，運算量即可大大減少。4.3.2基2時分的蝶形流圖與計算量分析

1.蝶形流圖

基2時分FFT蝶式運算定理的運算公式可以用圖4.3.1(a)所示的信號流圖表示，由于這個圖形呈蝶形，故稱為蝶形流圖，圖(b)是其簡化形式。一個蝶形流圖是基2時分算法的一個基本單元。圖4.3.1DIT-FFT蝶式運算信號流圖(a)蝶形流圖;(b)蝶形流圖簡化形式圖中左面兩支為輸入，中間以一個點表示加、減運算，右上支為相加輸出，右下支為相減輸出。如果在某一支路上信號需要進行乘法運算，則在該支路上標以箭頭，并將

相乘的系數標在箭頭邊。這樣，蝶式運算定理的運算公式(4.3.3)所表示的運算，可用圖4.3.1(b)所表示的“蝶形結”來表示。采用這種表示法，可將基2時分的分解過程用計算流圖來表示。

【例4.3.1】將8點序列的DFT用基2時分FFT算法進行分解。

解

(1)第一次分解。

將8點序列x(n)分解為兩個4點序列x1(n)和x2(n)，x1(n)為偶數序列，x2(n)為奇數序列，即

x1(n)={x(0)，x(2)，x(4)，x(6)}

x2(n)={x(1)，x(3)，x(5)，x(7)}

4點序列x1(n)和x2(n)分別做4點DFT得到X1(k)和X2(k)，由X1(k)和X2(k)通過蝶形構造獲得8點DFTX(k)。第一次分解的旋轉因子為(k=0，1，2，3)，運算流圖如圖4.3.2所示。

(2)第二次分解。將兩個4點序列x1(n)和x2(n)分解為四個2點序列x3(n)、x4(n)、和。

圖4.3.2

N點DFT基2時分第一次分解運算流圖(N=8)

2點序列分別做2點DFT得到X3(k)、X4(k)、、

，由四個2點DFT通過蝶形構造獲得兩個4點DFT，再通過蝶形構造獲得8點DFTX(k)。第二次分解的旋轉因

子為，運算流圖如圖4.3.3所示。圖4.3.3

N點DFT基2時分第二次分解運算流圖(N=8)

(3)第三次分解。

將四個2點序列分解為8個1點序列，其中

由于1點序列的DFT值為其序列本身，因此在最后一次分解后，流圖中已經沒有直接計算DFT的環(huán)節(jié)。第三次分解旋轉因子為，運算流圖如圖4.3.4所示。

由以上例子我們可以看出，由于每一次分解都是按輸入序列在時域上的次序是偶數還是奇數來抽取的，最終分解成N個1點DFT，因此稱為基2時分FFT。圖4.3.4

N點DFT基2時分FFT運算流圖(N=8)

2.計算量分析

對于任何一個2的整數冪N=2M，總是可以通過M次分解直到1點的DFT運算。這樣的M次分解，就構成了從x(n)到X(k)的M級運算過程。從圖4.3.4可以看出，每一級運算都由個蝶形運算構成。由于每一個蝶形運算需要一次復數乘法和兩次復數加法，因此每一級分解所需的復數乘法和加法次數分別為

N點DFT(M級分解)所需總的復數乘法和加法次數分別為

(4.3.4)

(4.3.5)我們知道直接計算N點DFT需要的復數乘法為N2次，復數加法為N(N-1)次，當N1時，基2時分FFT算法將比直接計算DFT的運算量大大減少?；?時分算法與直接計算N點DFT的復數乘法次數比為

(4.3.6)復數加法次數比為

(4.3.7)若N=210，則直接計算需要的復數乘法和加法次數分別為CM=10242=1048576，CA=1024(1024-1)=1047552，基2時分需要的復數乘法和加法次數分別為CM=5120，CA=10240,復數乘法次數比為αM≈，即基2時分FFT的復數乘法運算效率提高了約200倍。圖4.3.5顯示了直接計算DFT與基2時分FFT算法所需乘法次數的比較曲線，顯然，N越大，F(xiàn)FT算法的效率越高。圖4.3.5直接計算DFT與基2時分FFT算法所需乘法次數的比較曲線

4.2.4基2時分FFT算法的運算規(guī)律及編程思想由基2時分FFT運算流圖，可以看出這種算法的運算規(guī)律。

1.原位計算

由圖4.3.4的基2時分FFT運算流圖可以看出，同一級運算中，每個蝶形的兩個輸入數據只對計算本蝶形有用，而且每個蝶形的輸入、輸出數據結點又同在一條水平線上，這就意味著計算完一個蝶形后，所得輸出數據可立即存入原輸入數據所占用的存儲單元。這樣，經過M級運算后，原來存放輸入序列數據的N個存儲單元中便依次存放X(k)的N個值。

這種利用同一存儲單元存儲蝶形計算輸入、輸出數據的方法稱為原位運算。

2.旋轉因子的變化規(guī)律

在基2時分FFT運算流圖中，旋轉因子的變化和兩個節(jié)點之間的距離都呈現(xiàn)一定的規(guī)律，p稱為旋轉因子的指數。設L表示從左到右的運算級數，即L=1，2，…，M。設N=23=8，L=1，2，3，第L級共有2L-1個不同的旋轉因子，即

L=1時，;

L=2時，;

L=3時，。

對N=2M的一般情況，L=1，2，…，M，第L級的旋轉因子為

(4.3.8)在實際編程實現(xiàn)中，旋轉因子有三種計算方法。

(1)直接計算。在用高級編程語言實現(xiàn)FFT，且對運算時間沒有要求的理論分析場合下，

一般采用直接計算求旋轉因子。根據歐拉公式有

如果采用相同的旋轉因子計算完后存儲起來的辦法，則共需計算個，即N個正弦函數。

(2)遞推計算。根據旋轉因子的特點，即

可以利用遞推的辦法求解旋轉因子。根據上式，在每一級只需計算一個旋轉因子，其余的旋轉因子均可通過上式求出。因此N=2M點FFT共需計算2M個正弦函數。

(3)查表。在要求實時計算FFT的應用場合，用DSP匯編或FPGA編程實現(xiàn)FFT算法時，一般用查表方法。由于N點基2時分FFT算法中共有個不同的旋轉因子，因此事先計算好個旋轉因子表以供查詢，可以省去計算正弦函數的時間。

3.輸入序列的逆序

基2時分FFT算法輸入序列的排序看起來似乎很亂，但仔細分析就會發(fā)現(xiàn)這種排序是

有規(guī)律的，我們稱之為比特逆序。

對于N=2M，N個順序數可以用M位二進制數(nM-1nM-2

…n1n0)2來表示，即(n)10=(nM-1nM-2…n0)2，定義

ρM(n)=(n0n1…nM-1)2

(4.3.9)

為n的M位比特逆序數。

【例4.3.2】計算n=2，12的二進制比特逆序數ρ4(n)。

解

造成基2時分FFT算法的輸入序列順序逆序的原因是將輸入序列x(n)按下標n的奇偶不斷分解造成的。表4.3.1列出了N=8的自然順序數和比特逆序數。由表顯而易見，只要

將自然順序的二進制數按位倒置，就能得到比特逆序數。

在實際應用中，一般先按自然順序數將輸入序列存入存儲單元A(m)，m=0，1，…，N-1。

為了得到比特逆序的排列，可以通過變址運算來完成，變址過程如圖4.3.6所示。表4.3.1N=8的自然順序數和比特逆序數圖4.3.6比特逆序的變址過程(N=8)由表4.3.1可以看出，自然順序的次序增加是在二進制數低位加1，二進制進位由低向高;比特逆序的次序增加在二進制數高位加1，進位由高向低。因此，實際求數

n(n=0，1，2，…，N-1)的比特逆序ρM(n)步驟如下:

(1)0的比特逆序為0，因此ρM(0)=0。

(2)求n+1的逆序ρM(n+1)。如ρM(n)的二進制數最高位為0，則將該位賦1，循環(huán)結束。

如ρM(n)的二進制數從高位開始數前幾位均是1，則把為1的最低位的右邊為0的位賦1，同時將為1的位賦0，循環(huán)結束。

綜合以上求比特逆序步驟和圖4.3.6，將自然順序的輸入序列通過比特逆序變址變?yōu)楸忍啬嫘蜉斎氪涡虻某绦蛄鲌D如圖4.3.7所示。圖4.3.7比特逆序程序流圖圖4.3.7所示的程序流圖是求j的逆序，同時將存儲單元內自然順序排列的輸入序列變?yōu)槟嫘蚺判?，i存儲的是j-1的逆序數。由于j=0，N-1的逆序既是數據本身，因此程序循

環(huán)從j=1到j=N-2。s存放的是二進制數每一位的權，每次j循環(huán)開始，s為二進制最高位的權。程序流圖中虛線框部分是將為1的位賦0的過程。

4.編程思路與程序流圖

觀察基2時分FFT蝶形流圖，可歸納出一些對編程有用的運算規(guī)律：第L級運算中，每個蝶形的兩個輸入數據相距2L-1個點，同一個旋轉因子對應著間隔為2L點的2M-L個蝶形?？偨Y上述運算規(guī)律，可以采用以下方法。首先從第一級(輸入端)開始，逐級進行，共進行M級運算。在進行第L級運算時，依次求出2L-1個不同的旋轉因子，每求出一個旋轉因子，就計算完它對應的所有2M-L個運算蝶。因此基2時分FFT程序由逆序和三重循環(huán)組成。三重循環(huán)與流圖的對應關系為：最外層循環(huán)(循環(huán)變量為L)對應蝶形流圖的級數M，中間循環(huán)(循環(huán)變量為k)對應流圖中每一級不同旋轉因子對應的運算蝶數目2L-1，最內層循環(huán)(循環(huán)變量為J)對應每一級相同旋轉因子運算蝶的順序，程序流圖如圖4.3.8所示，C語言程序見附錄。圖4.3.8基2時分FFT算法程序流圖

【例4.3.3】設x(n)為N點有限序列，N為偶數，其N點DFT為X(k)，試用X(k)表示如下序列的N點DFT。

解

(1)解法1。

(2)解法2。用基2時分蝶式運算定理將x2(n)分解為兩個

點的序列x20(n)和x21(n)。

根據蝶式運算定理有

因此有

同理，將x(n)分解為兩個點的序列x0(n)和x1(n)。

因為

所以

令，有

因為

所以

4.4.1基2頻分蝶式運算定理

設序列x(n)的長度為N=2M，其N點DFT為X(k)=DFT［x(n)］，0≤k≤N-1，按k的奇偶將X(k)分解為兩個點的子序列X1(k)和X2(k)。

4.4基2頻分FFT算法若x1(n)=IDFT［X1(r)］，x2(n)=IDFT［X2(r)］，0≤n≤

-1，則

(4.4.1)

(4.4.2)證明由證明過程可得，N點序列x(n)按式(4.4.1)和式(4.4.2)分解為兩個點的子序列，這兩個點子序列的DFT正好是N點DFTX(k)的偶數組和奇數組。由蝶式運算定理可以看出，要完成一個蝶形運算，需要一次復數乘法和兩次復數加法，計算N點DFT共需計算兩個點DFT和次蝶形運算。計算點DFT需要次復數乘法和次復數加法。當N>>1時，計算N點DFT共需

次復數加法。由此可見，僅經過一次分解，就使運算量減少近一半。如果N=2M，則我們可以繼續(xù)分解，直至分解

到1點DFT，運算量即可大大減少。4.4.2基2頻分的蝶形流圖與計算量分析

1.蝶形流圖

基2頻分FFT蝶式運算定理的運算公式可以用圖4.4.1所示的信號流圖表示，一個蝶形流圖是基2頻分算法的一個基本單元。圖4.4.1DIF-FFT蝶式運算信號流圖

【例4.4.1】將8點序列的DFT用基2頻分FFT算法進行分解。

解

(1)第一次分解。

將8點序列x(n)通過蝶形構造分解為兩個4點序列x1(n)和x2(n),分別做4點DFT得到X1(k)和X2(k)，即

X1(k)={X(0)，X(2)，X(4)，X(6)}

X2(k)={X(1)，X(3)，X(5)，X(7)}

由X1(k)和X2(k)獲得8點DFTX(k)，旋轉因子為(n=0，1，2，3)，運算流圖如圖4.4.2所示。圖4.4.2

N點DFT基2頻分第一次分解運算流圖(N=8)

(2)第二次分解。

將兩個4點序列分別通過蝶形構造分解為四個2點序列，分別做2點DFT得到

由2點DFT即可獲得8點DFTX(k)，第二次分解的旋轉因子為(n=0，1)，運算流圖如圖4.4.3所示。圖4.4.3

N點DFT基2頻分第二次分解運算流圖(N=8)

(3)第三次分解。

將四個2點序列通過蝶形構造分解為8個1點序列，其中

由于1點序列的DFT值為其序列本身，因此，在最后一次分解后，流圖中已經沒有直接計算DFT的環(huán)節(jié)。第三次分解旋轉因子為，運算流圖如圖4.4.4所示。圖4.4.4

N點DFT基2頻分FFT運算流圖(N=8)由以上例子我們可以看出，由于每一次分解都是按輸入序列在頻域上的次序是偶數還是奇數來抽取的，最終分解成N個1點DFT，因此稱為基2頻分FFT。

2.計算量分析

對于任何一個長度為2的整數冪(N=2M)序列，其DFT總是可以通過M次分解到1點的DFT運算。這樣的M次分解，就構成了從x(n)到X(k)的M級運算過程。從圖4.4.4可以看出，每一級運算都由個蝶形運算構成。由于每一個蝶形運算需要一次復數乘法和兩次復數加法，因此N點DFT(M級分解)所需總的復數乘法和加法次數分別為

(4.4.3)

CA=NM=Nlog2N

(4.4.4)由于基2時分FFT算法和基2頻分FFT算法抽取的域不同，因此它們的蝶形構造發(fā)生的地方、基本運算和蝶形流圖都不相同?；?時分是時域抽取，頻域蝶形構造，輸入逆序，輸出順序;基2頻分是頻域抽取，時域蝶形構造，輸入順序，輸出逆序。但兩者的運算量是相同的。具體的編程思想不再贅述。

在前面章節(jié)的學習中我們或者使用基2時分FFT算法，或者使用基2頻分FFT算法將某一個序列的DFT計算出來。下面我們舉個例子，在同一個序列中某些步驟用基2時分，某些步驟用基2頻分。

【例4.4.2】在N=8的DFT中，分兩步做快速算法，第一步用基2頻分將8點DFT化為兩個4點DFT，第二步用基2時分將這兩個4點DFT化為四個2點DFT，并畫出流圖。

解根據題目要求，第一步使用基2頻分，構造4點序列x1(n)和x2(n)，其DFT分別是

X1(k)={X(0)，X(2)，X(4)，X(6)}，

X2(k)={X(1)，X(3)，X(5)，X(7)}

運算流圖如圖4.4.5所示。圖4.4.5例4.4.2第一步基2頻分FFT第一次分解運算流圖第二步使用基2時分，以x1(n)為例，將其按照n的奇偶分成2點序列x3(n)和x4(n),由2點DFTX3(k)和X4(k)通過蝶形運算即可構造出X1(k)，如圖4.4.6所示。同理可得X2(k)。

最后，將前兩步的結果合并起來，運算流圖如圖4.4.7所示。圖4.4.6例4.4.2第二步基2時分FFT分解運算流圖圖4.4.7例4.4.2運算流圖長度為N的有限長序列x(n)的離散傅里葉變換為

(4.5.1)

其離散傅里葉逆變換為

(4.5.2)4.5IDFT的快速算法比較以上兩式，可見離散傅里葉變換對具有對稱性。因此可以利用DFT的快速算法(FFT)來實現(xiàn)IDFT。具體來說，就是將FFT運算流圖中的旋轉因子取共軛(虛部取反)，最后輸出乘以。此時，輸入為X(k)，輸出即為x(n)。前面討論的序列x(n)的DFT快速算法都是復數運算，包括序列x(n)也認為是復數。在實際應用中，信號有可能是實序列，任何實數都可看成虛部為零的復數。例如，求某實信號x(n)的頻譜，可認為是將實信號加上數值為零的虛部變成復信號(x(n)+j0)，再用FFT求其離散傅里葉變換。這種做法很不經濟，因為把實序列變成復序列，存儲器要增加一倍，且計算機運行時，即使虛部為零，也要進行涉及虛部的運算，浪費了運算量。

合理的解決方法是利用復數FFT對實數進行有效計算。4.6實序列DFT的有效計算方法

1.一個N點FFT同時計算兩個N點實序列的DFT

【例4.6.1】設x1(n)和x2(n)是彼此獨立的兩個N點實序列，且X1(k)=DFT［x1(n)］，X2(k)=DFT［x2(n)］。

試通過一次FFT運算同時求得X1(k)和X2(k)。

解首先將x1(n)和x2(n)分別作為一個N點復序列x(n)的實部及虛部，即

x(n)=x1(n)+jx2(n)通過調用一個N點的FFT算法即可獲得x(n)的DFT序列X(k)。

利用離散傅里葉變換的共軛對稱性可獲得X1(k)和X2(k)，即

下面分析計算量。通過調用N點FFT算法求X(k)需要的復數乘法次數和復數加法次數分別為

根據X(k)求X1(k)和X2(k)需要的復數乘法次數和復數加法次數分別為

CM2=N-1，CA2=2(N-1)

因此，該方法總共需要的復數乘法次數和復數加法次數分別為

如果調用兩次N點FFT算法直接計算X1(k)和X2(k)，所需的運算量為

CM=Nlog2N，CA=2Nlog2N

設N=1024，則直接調用兩次FFT計算需要的CM=10240，CA=20480。利用上述算法則有CM=6143，CA=12286，計算量減少近半。

2.一個N點FFT計算一個2N點實序列的DFT

【例4.6.2】設x(n)是2N點的實序列，試設計用一次N點的FFT算法完成計算X(k)的高效算法。

解將x(n)分解為偶數序列x1(n)和奇數序列x2(n)，即

x1(n)=x(2n)，n=0，1，…，N-1

x2(n)=x(2n+1)，n=0，1，…，N-1將x1(n)和x2(n)分別作為一個N點復序列y(n)的實部及虛部，即

y(n)=x1(n)+jx2(n)

通過N點FFT運算可得到y(tǒng)(n)的N點DFTY(k)，根據前面的討論可得

根據基2時分蝶式運算定理得

下面分析計算量。通過調用N點FFT算法求Y(k)需要的復數乘法次數和復數加法次數分別為

根據Y(k)求X1(k)和X2(k)需要的復數乘法次數和復數加法次數分別為

CM2=N-1，

CA2=2(N-1)

利用基2時分蝶式運算定理求X(k)需要的復數乘法次數和復數加法次數分別為

CM3=N，

CA3=2N因此，該方法總共需要的復數乘法次數和復數加法次數分別為

直接計算2N點FFT計算量為

CM=Nlog22N，

CA=2Nlog22N

設N=1024，則直接計算CM=11264，CA=22528。

利用上述算法則有CM=7167，CA=14334。前面已經介紹過，采用FFT算法可以很快計算出全部N點DFT值，即Z變換X(z)在z平面單位圓上的全部等間隔取樣值。而實際中也許不需要計算整個單位圓上Z變換的取樣值。例如對于窄帶信號，只需要對信號所在的一段頻帶進行分析，這時希望頻譜的取樣集中在這一頻帶內，以獲得較高的分辨率，而頻帶以外的部分可不考慮。另外，實際中有時也對非

單位圓上的取樣值感興趣。4.7線性調頻Z(Chirp-Z)變換算法例如需要知道Z變換的極點所在處的頻率，如極點位置離單位圓較遠，則其單位圓上的頻譜就很平滑，這時很難從中識別出極點所在處的頻率，此時，如果取樣不是沿單位圓而是沿一條接近這些極點的弧線進行，則極點所在頻率上的頻譜將出現(xiàn)明顯的尖峰，由此可較準確地測定極點頻率。

Z變換采用螺旋線取樣是一種適合于以上需要的變換，且可以采用FFT來快速計算，這種變換也稱做線性調頻Z變換(也稱Chirp-Z變換或CZT)，它是適用于這種更為一般情況下由x(n)求X(zk)的快速算法。4.7.1Chirp-Z變換的基本原理

已知x(n)(0≤n≤N-1)為有限長序列，其Z變換為

(4.7.1)

為適應z可以沿z平面上的一段螺旋線做等分角的取樣，令z的取樣點為

zk=AW-k，k=0，1，…，M-1

(4.7.2)其中，M為所要分析的復頻譜的點數(不一定等于N);A和W是任意復數，可表示為

(4.7.3)

(4.7.4)圖4.7.1CZT在z平面的螺旋線取樣將式(4.7.3)和式(4.7.4)帶入式(4.7.2)中，可得

(4.7.5)

有

取樣點在z平面上所沿的周線如圖4.7.1所示。由以上討論及圖可見:

(1)A0表示起始取樣點的矢量半徑長度，通常A0≤1，否則z0將處于單位圓外。

(2)ω0表示起始取樣點z0的相角，可以是正值也可以是負值。

(3)φ0表示兩相鄰取樣點之間的等分角度，φ0為正，表示zk的路徑是逆時針旋轉的;φ0為負，表示zk的路徑是順時針旋轉的。

(4)W0的大小表示螺旋線的伸展率，當W0<1時，隨著k的增加螺旋線外伸;當W0>1時,隨著k的增加螺旋線內縮；W0=1表示半徑為A0的一段圓弧，若A0=1，這段圓弧則是單位圓的一部分。當

時，各個zk均勻等間隔地分布在單位圓上，此時CZT即為DFT。

將式(4.7.2)中的zk代入Z變換表達式(4.7.1)中，可得

(4.7.6)

顯然，同直接計算DFT情況相仿，按照式(4.7.6)計算出全部M點取樣值需要NM次復數乘法和(N-1)M次復數加法，當N及M較大時，計算量迅速增加。如果通過一定的變換，如采用布魯斯坦(Bluestein)等式，將以上運算轉換為卷積形式，則可采用FFT進行，這樣即可大大提高計算速度。布魯斯坦等式為

(4.7.7)將上式代入式(4.7.6)，可得

令

(4.7.8)

(4.7.9)

則有

(4.7.10)式(4.7.10)說明，如對信號x(n)先進行一次加權處理，加權系數為，然后通過一個單位脈沖響應為h(n)的線性時不變系統(tǒng)，最后對該系統(tǒng)的前M點輸出再作一次的加權，就可得到全部M點螺旋線取樣值。運算流程如圖4.7.2所示。圖4.7.2線性調頻Z變換運算流程由于系統(tǒng)的單位脈沖響應與頻率隨時間成線性增加的線性調頻信號(即ChirpSignal)相似，因此這種算法稱為Chirp-Z變換。

4.7.2Chirp-Z變換的實現(xiàn)步驟

由式(4.7.10)可以看出，由于輸入信號g(n)是長度為N的有限長序列，序列是無限長的，而計算0~M-1點卷積g(k)*h(k)所需要的h(n)是取值在n=-(N-1)~(M-1)那一部分的值，如圖4.7.3(a)所示。因此，可認為h(n)是一個有限長序列，長度為L=N+M-1。所以，Chirp-Z變換為兩個有限長序列的線性卷積

g(k)*h(k),可用循環(huán)卷積通過FFT來實現(xiàn)，即

g(k)*h(k)=IDFT［G(r)H(r)］根據第3章所學知識，g(k)*h(k)的長度為2N+M-2，因而用循環(huán)卷積求線性卷積時，循環(huán)卷積的長度最短為2N+M-2。由于我們只需要求0~M-1處的線性卷積結果，因此選循環(huán)卷積的長度為兩個序列中最長的長度，即L=N+M-1點。另外，h(n)的主值序列(0≤n≤L-1)可由h(n)作周期延拓后取0≤n≤L-1部分值獲得，如圖4.7.3(b)所示。

將h(n)與g(n)作循環(huán)卷積后，其輸出的前M個值就是CZT的M個值。這個循環(huán)卷積的過程可在頻域上通過FFT實現(xiàn)。圖4.7.3

h(n)的主值序列求解過程總之，CZT運算步驟如下：

(1)選擇一個最小整數L，使其滿足L≥N+M-1，同時滿足L=2m，以便采用基2FFT算法。

(2)對x(n)加權并補零得到g(n)：

(4.7.11)并利用FFT算法求此序列的L點DFTG(r)。

(4.7.12)

(3)求h(n)的主值序列:

(4.7.13)并利用FFT算法求此序列的L點DFTHm(r):

(4.7.14)

(4)將Hm(r)和G(r)相乘，得到Q(r)=G(r)Hm(r)。

(5)用FFT算法求L點Q(r)的IDFT，得到hm(n)與g(n)的循環(huán)卷積。

(4.7.15)其中，q(n)的前M-1點等于h(n)與g(n)的線性卷積結果(與第3章所講的線性卷積與循環(huán)卷積的關系不同，這里的hm(n)是h(n)以L為周期延拓后的主值序列)。

(6)求X(zk):

(4.7.16)4.7.3Chirp-Z變換運算量的估算

采用上小節(jié)的算法后，CZT所需的運算量如下:

(1)形成L點序列，由于補零關系，

g(n)只有N點有序列值，需要N次復數乘法。系數可以通過遞推求得。令

則

初始條件。

因此通過遞推求Cn需要2N次復數乘法。

(2)形成L點序列hm(n)，由于它是由在

-N+1≤n≤M-1段內的序列值構成的，是偶對稱序列，如果設N≥M，則只需求0≤n≤N-1段內的N點序列值。求

也可以用遞推方式求解，因此需要復數乘法2N次。

(3)計算Hm(r)、G(r)和q(n)共需要3次L點FFT算法，因此需要復數乘法次數為。

(4)計算Q(r)=G(r)Hm(r)需要L次復數乘法。

(5)計算X(zk)=，0≤k≤M-1需要M次復數乘法。綜上所述，Chirp-Z變換總共需要的復數乘法次數為

(4.7.17)

前面討論過直接計算X(zk)需要復數乘法為NM次，當N=925，M=100時，直接計算需要CM=92500次，采用快速算法后所需的復數乘法次數為CM=21109，計算量減少較多?？傊?，與DFT相比，CZT有以下特點：

(1)輸入序列的長度N與輸出序列的長度M不需要相等；(2)N及M不必是合成數，二者均可為素數；

(3)zk點的角度間隔φ0是任意的，因此頻率分辨率也是任意的；

(4)周線不必是z平面上的圓；

(5)起始點z0可任意選定，因此可以從任意頻率上開始對輸入數據進行窄帶高分辨率分析；

(6)若A=1，M=N，即使N為素數，也可用Chirp-Z變換計算DFT。4.7.4用Matlab計算Chirp-Z變換

Matlab信號處理工具箱提供了內部函數czt用于實現(xiàn)Chirp-Z變換，調用格式如下：

y=czt(x，M，W，A);

y=czt(x);

y=czt(x，M，W，A)用于計算指定參數M、W、A下的Chirp-Z變換，而y=czt(x)則使用默認參數進行Chirp-Z變換，默認參數為M=length(x)，W=exp(j*2*pi/m)，A=1，此時即計算DFT。

【例4.7.1】某一低通濾波器的零、極點分布如圖4.7.4所示，利用Chirp-Z變換觀察濾波器零點特性。

clear；fs=2000; ％采樣頻率為2000Hz

fp=150；

p1=0.9*exp(-j*2*pi*fp/fs);

p2=p1′;p=［p1，p2］′; ％極點在150Hz處

fz1=500;

fz2=300;

z1=1.2*exp(-j*2*pi*fz1/fs)；z2=z1′;

z3=1.2*exp(-j*2*pi*fz2/fs);z4=z3′;圖4.7.4濾波器零、極點分布圖

z=［z1，z2z3z4］′;

％零點在300Hz和500Hz處的單位圓外，半徑為1.2

［b，a］=zp2tf(z，p，1);

ww=0：0.005*pi：2*pi;

h1=freqz(b，a，ww);

hn1=ifft(h1);hn=real(hn1)；

N=length(ww)； ％濾波器單位脈沖響應長度(401點)

f2=700;f1=50；M=201；

W=exp(-j*2*pi*(f2-f1)/(M*fs));

％在z平面ω=2πf1~2πf2之間取M點

A=1.19*exp(j*2*pi*f1/fs); ％取樣圓弧半徑為1.19

yzz=czt(hn，M，W，A); ％求CZT由于零點在單位圓外，半徑為1.2處，因此CZT的圓弧線半徑為1.19，起始點ω=2π×50,終止點為ω=2π×700，取樣點數為201點，如圖4.7.5所示。濾波器單位圓上的幅頻響應如圖4.7.6所示，由于單位圓距離零點較遠，兩個零點處的頻率在幅頻響應里不明顯。由于CZT取樣曲線接近零點，因此在300Hz和500Hz處的頻率在CZT中非常明顯，CZT如圖4.7.7所示。由于CZT的取樣圓弧離極點位置較遠，因此極點處的頻率在CZT曲線中已不明顯。圖4.7.5CZT取樣螺旋線圖圖4.7.6濾波器幅頻響應圖圖4.7.7濾波器CZT圖用Matlab按快速算法的實現(xiàn)步驟編寫的CZT程序如下(參數設置與上相同):

L=M+N-1;

gn=zeros(1，L);

gn(1：N)=A.^(-(0：N-1)).*W.^((0：N-1).^2./2).*hn;

gf=fft(gn，L);

hhn(1：M)=W.^(-(0：M-1).^2./2);

hhn(M+1：L)=W.^(-(L-(M：L-1)).^2./2);

hhnf=fft(hhn);

Qr=hhnf.*gf;

qn=ifft(Qr);

yzz=qn(1：M).*W.^((0：M-1).^2./2)

【例4.7.2】設信號為x(n)，它由4個頻率分別為60Hz、64Hz、95Hz和100Hz的正弦序列組合而成。采樣頻率為600Hz，樣點數為200點。試用CZT觀察信號頻譜。

解為比較起見，分別用200點的DFT和CZT觀察信號頻譜。

圖4.7.8所示為x(n)的DFT，由于點數太少，60Hz與64Hz、95H