《MATLAB輔助現(xiàn)代工程數(shù)字信號處理》課件第2章

第2章離散時間信號與系統(tǒng)的頻域分析

2.1離散時間序列的傅里葉變換(DTFT)

2.2離散周期信號的傅里葉級數(shù)(DFS)

2.3離散傅里葉變換(DFT)

2.4快速傅里葉變換(FFT)

2.5Z變換

2.6線性時不變離散系統(tǒng)的頻域分析

2.7小結(jié)2.1離散時間序列的傅里葉變換(DTFT)

如果信號在頻域上是離散的，則信號在時域上就表現(xiàn)為周期性的時間函數(shù)。相反，如果信號在時域上是離散的，則其在頻域上就表現(xiàn)為周期性的頻率函數(shù)。因此，一個離散周期序列的頻譜一定既是周期的又是離散的。本章先介紹離散時間傅里葉變換(DiscreteTimeFourierTransform，DTFT)，再介紹離散周期序列及其傅里葉級數(shù)(DFS)。2.1.1DTFT的定義

序列x(n)的離散時間傅里葉變換(DTFT)定義為(2.1)由X(ej(ω+2π))=X(ejω)，故X(ejω)是ω的周期函數(shù)，周期為2π。DTFT成立的充分必要條件是序列x(n)絕對可和，即滿足(2.2)DTFT存在逆變換，稱為離散時間傅里葉逆變換(IDTFT):(2.3)式(2.1)和式(2.3)為序列x(n)的傅里葉變換對，即離散時間信號的傅里葉變換對。式(2.1)為正變換，式(2.3)為反變換或稱為逆變換。

式(2.1）又可視為X(ejω)的傅里葉級數(shù)展開式，式(2.3)確定的x(n)是傅里葉級數(shù)的系數(shù)。2.1.2DTFT的性質(zhì)

1.線性性

定義信號x1(n)和x2(n)的DTFT分別是X1(ejω)和X2(ejω)，并令x(n)=ax1(n)+bx2(n)，其中，a和b為常數(shù)，則

X(ejω)=aX1(ejω)+bX2(ejω)

(2.4)

2.時移與頻移性

令y(n)=x(n－n0)，則

Y(ejω)=e－jωn0X(ejω)

(2.5)

3.對稱性

將序列x(n)分成實部與虛部，即

x(n)=xR(n)+jxI(n)

(2.6)

將式(2.6)進行DTFT，得到

X(ejω)=XR(ejω)+jXI(ejω)

(2.7)

X(ejω)的實部XR(ejω)是偶函數(shù)，虛部XI(ejω)是奇函數(shù)，即

XR(ejω)=XR(e－jω)

(2.8)

XI(ejω)=－XI(e－jω)

(2.9)

4.時域卷積定理

設(shè)y(n)=x(n)*h(n)，則

Y(ejω)=X(ejω)·H(ejω)

(2.10)

該定理說明，兩序列卷積的DTFT服從相乘的關(guān)系。

線性時不變系統(tǒng)輸出y(n)的DTFT等于輸入信號x(n)的DTFT乘以系統(tǒng)單位脈沖響應(yīng)h(n)的DTFT。因此，求系統(tǒng)的輸出信號可以在時域用卷積公式計算，也可以在頻域按照式(2.10)，求出輸出的DTFT，再作逆DTFT求出輸出信號。

5.頻域卷積定理

若y(n)=x(n)h(n)，則

(2.11)

6.時域相關(guān)定理

若y(n)是x(n)和h(n)的相關(guān)函數(shù)，即則

7.帕斯維爾(Parseval)定理帕斯維爾定理如下：(2.13)帕斯維爾定理指出，信號時域的總能量等于頻域的總能量。需說明的是，這里頻域總能量是指|X(ejω)|2在一個周期中的積分再乘以1/2π。因此，|X(ejω)|2是信號的能量譜。|X(ejω)|2/2π為信號的能量譜密度。

8.維納-辛欽(Wiener-Khinchin)定理

若x(n)是功率信號，則其自相關(guān)函數(shù)的傅里葉變換為

(2.14)式中，X2N(ejω)為(2.15)若式(2.14)的右邊極限存在，則稱該極限為功率信號x(n)的功率譜Px(ejω)。式(2.14)稱為確定性信號的維納-辛欽定理，它說明功率信號x(n)的自相關(guān)函數(shù)和其功率譜是一對傅里葉變換。信號的總功率為(2.16)無論x(n)是實信號還是復(fù)信號，其功率譜Px(ejω)始終是ω的實函數(shù)，即功率譜失去了相位信息。相關(guān)函數(shù)和功率譜是描述隨機信號的重要統(tǒng)計量。

2.2離散周期信號的傅里葉級數(shù)(DFS)

2.2.1DFS的定義

一個周期為N的周期序列，即，k為任意整數(shù)，N為周期。該序列可用離散傅里葉級數(shù)來表示，即可用周期為N的正弦序列來表示

(2.17)將上式兩邊乘以，并對一個周期求和，有由于(2.18)故可得到

(2.19)是一個由N個獨立諧波分量組成的傅里葉級數(shù)，其周期為N。因此，時域上周期序列的離散傅里葉級數(shù)在頻域上仍是一個周期序列。式(2.17)和式(2.19)稱為周期序列的離散傅里葉級數(shù)(DFS)變換對，可表示為(2.20)(2.21)記，則DFS變換對可寫為

(2.22)(2.23)

DFS變換對公式表明，對于一個無窮長周期序列，只要知道一個周期的信號變化情況，就可以知道其他周期的情況。所以，這種無窮長序列實際上只有N個序列值的信息是有用的。因此，周期序列與有限長序列有著本質(zhì)的聯(lián)系。為了便于MATLAB實現(xiàn)計算，又因和均是周期為N的周期函數(shù)，故可設(shè)和代表序列和的主值區(qū)間序列，則式(2.22)與式(2.23)可寫為(2.24)(2.25)式中，(2.26)矩陣WN為正交酉矩陣，稱作DFS矩陣;W*N表示矩陣WN的復(fù)共軛。

【例2.1】

設(shè)，以N=4為周期進行周期延拓，求周期序列的離散傅里葉級數(shù)。

MATLAB程序如下:

%MATLABPROGRAM2-1

%實現(xiàn)離散傅里葉級數(shù)(DFS)的計算

%xn代表離散時間序列x(n)，N為離散時間序列x(n)的長度

%Xk為x(n)的傅里葉級數(shù)，且為其主值序列

xn=[0，1，2，3];N=4;%設(shè)定序列和周期

n=[0:1:N-1];k=[0:1:N-1];

%設(shè)定n和k

Wn=exp(-j*2*pi/N);

%設(shè)定Wn因子

nk=n′*k;

Wnnk=Wn.^nk;

%計算W矩陣

Xk=xn*Wnnk;

%計算DFS的系數(shù)Xk

disp(′xn=′);disp(xn);

disp(′Xk=′);disp(Xk);%顯示計算結(jié)果(系數(shù))

程序運行結(jié)果為

xn=

0

123

Xn=

6.0000－2.0000+2.0000i－2.0000－0.0000i－2.0000－2.0000i

【例2.2】

利用例2.1的計算結(jié)果，求離散傅里葉級數(shù)反變換。

MATLAB程序如下:

%MATLABPROGRAM2-2

%實現(xiàn)離散傅里葉級數(shù)反變換(IDFS)的計算

%Xk為x(n)的傅里葉級數(shù)，且為其主值序列，N為Xk的長度

%xn代表離散時間序列x(n)

Xk=[6.00－2.00+2.00i－2.00－0.00i－2.00－2.00i];%設(shè)定DFS主值序列

N=4;

n=[0:1:N－1];k=[0:1:N－1];%設(shè)定n和k

Wn=exp(－j*2*pi/N);

%設(shè)定Wn因子

nk=n′*k;

Wnnk=Wn.^(－nk);

%計算W矩陣

xn=Xk*Wnnk/N;

%計算xn，注意矩陣相乘順序

disp(′xn=′);disp(xn);

%顯示計算結(jié)果(系數(shù))程序運行結(jié)果為

xn=

0

1.0000－0.0000i

2.0000+0.0000i

3.0000+0.0000i

【例2.3】

設(shè)

其中，N為序列周期，L/N是占空比。繪出L=10，N=80的幅度和角度樣本。(m為整數(shù))

MATLAB程序如下:

%MATLABPROGRAM2-3

clc;

L=10;N=80;

n1=[0:L-1];

xn1=0.8*n1.*exp(-0.4*n1);

xn2=zeros(1，N-L);

xn=[xn1，xn2];

n=[0:1:N-1];k=[0:1:N-1];

Wn=exp(-j*2*pi/N);

nk=n′*k;

Wnnk=Wn.^nk;

Xk=xn*Wnnk;

magXk=fftshift(abs(Xk));

AngleXk=angle(Xk);

figure(1);

stem(k，magXk);

xlabel(′k′);ylabel(′|Xk|′);title(′AmplitudeofDFS′);

grid;

figure(2);

stem(k，AngleXk);

xlabel(′k′);ylabel(′Angle(Xk)′);title(′AngleofDFS′);

grid

程序運行結(jié)果如圖2.1所示。圖2.1離散傅里葉級數(shù)的幅值譜和相位譜2.2.2DFS的性質(zhì)

DFS在時域和頻域之間具有嚴格的對偶關(guān)系，下面列出幾個常用的性質(zhì)。

假設(shè)和都是周期為N的兩個周期序列，其離散傅里葉級數(shù)分別為DFS[

]和DFS[

]。

1.線性性

線性的定義如下：(2.27)式中，a，b為任意常數(shù)。

2.序列周期移位

由于及WNkn都是以N為周期的函數(shù)，故有

(2.28)

3.共軛對稱性對于復(fù)序列及其共軛復(fù)序列，有(2.29)復(fù)序列的離散傅里葉級數(shù)的共軛偶對稱分量為(2.30)復(fù)序列的離散傅里葉級數(shù)的共軛奇對稱分量為(2.31)

4.周期卷積

若，則(2.32)周期卷積中的和都是周期為N的周期序列，因而乘積也是周期為N的周期序列，且周期卷積中的求和只在一個周期上進行。

2.3離散傅里葉變換(DFT)

在計算機上實現(xiàn)信號頻譜分析時，要求信號在時域和頻域都是離散的，且都應(yīng)是有限長的。離散傅里葉級數(shù)DFS在時域和頻域都是離散的，但和都是以N為周期的無限長周期序列，因此，可取其中的一個周期作為主值序列樣本，形成離散傅里葉變換對，即DFT。DFS則是DFT進行周期延拓的結(jié)果。2.3.1DFT的定義

由離散傅里葉級數(shù)變換對(式(2.17)和式(2.19))可知，用離散周期序列中一個周期的N個樣本就能確定頻譜序列;同樣，只用一個周期中的N根譜線也可確定離散周期序列。在一個周期內(nèi)的N個樣本構(gòu)成有限長序列x(n)，在一個周期內(nèi)的N個樣本構(gòu)成有限長序列X(k)。x(n)與X(k)之間的變換關(guān)系定義為有限長序列的離散傅里葉變換(DiscreteFourierTransform，DFT)。

設(shè)x(n)是一個長度為N的有限長序列，則x(n)的離散傅里葉變換(DFT)為

(2.33)

X(k)的離散傅里葉逆變換(IDFT)為

(2.34)式中，

，通常稱式(2.33)和式(2.34)為離散傅里葉變換對。若把有限長序列x(n)看成周期序列的的主值序列，則相應(yīng)的離散傅里葉變換X(k)為相應(yīng)的周期序列的傅里葉級數(shù)的主值序列，即(2.35)(2.36)式中，RN(n)表示長度為N的矩形序列，見式(1.12)。離散傅里葉變換是一個長度為N的序列，對應(yīng)的離散頻率ω在0～2π之間，間隔相等，為2π/N。離散傅里葉變換具有唯一性，其物理意義表示序列x(n)的Z變換在單位圓上的等角距取樣。

在MATLAB中，用向量x和X分別表示序列x(n)和X(k)，類似于DFS的定義，可得矩陣表達式為

X=WNx

(2.37)

(2.38)式中，WN和W*N即為式(2.26)定義的變換矩陣，此處稱為DFT矩陣。

【例2.4】用MATLAB編寫函數(shù)dft和idft，分別實現(xiàn)DFT和IDFT的計算。MATLAB程序如下:

%MATLABPROGRAM2-4

functionXk=dft(xn，N)

%計算離散傅里葉變換

%Xk為離散序列x(n)的離散傅里葉變換

%xn為N點有限長序列x(n)

%N為離散序列x(n)的長度

n=[0:1:N－1];

k=[0:1:N－1];

Wn=exp(－j*2*pi/N);%Wn因子

nk=n′*k;

%產(chǎn)生一個含nk值的N維方陣

Wnnk=Wn.^nk;

%DFT矩陣

Xk=xn*Wnnk;

%DFT系數(shù)的行向量

%-----------------------------------------------------------------------------

functionxn=idft(Xk，N)

%計算逆離散傅里葉變換

%xn為N點有限長序列

%Xk為在區(qū)間[0，N-1]的DFT系數(shù)數(shù)組

%N為DFT的長度

n=[0:1:N－1];

%n的行向量

k=[0:1:N－1];

%k的行向量

Wn=exp(－j*2*pi/N);%Wn因子

nk=n′*k;

%產(chǎn)生一個含kp值的N乘N維矩陣

Wnnk=Wn.^(－nk);%IDFT矩陣

xn=Xk*Wnnk/N;%IDFT的行向量

【例2.5】

序列x(n)=sin(0.84πn)+sin(0.88πn)，0≤n≤90，試繪制x(n)及其離散傅里葉變換X(k)的幅值譜和相位譜。MATLAB程序如下:

%MATLABPROGRAM2-5

clc;

N=90;

n=[0:N－1];

xn=sin(0.84*pi*n)+sin(0.88*pi*n);

Xk=dft(xn，N);

magXk=abs(Xk);AngleXk=angle(Xk);

figure(1);

stem(n，xn);

xlabel(′n′);ylabel(′x(n)′);title(′序列x(n)′);

figure(2);

k=[0:length(magXk)－1];

stem(k，magXk);

xlabel(′k′);ylabel(′|X(k)|′);title(′DFT的幅值譜′);

figure(3);

stem(k，AngleXk);

xlabel(′k′);ylabel(′Angle(Xk)′);title(′DFT相位譜′);

程序運行結(jié)果如圖2.2～圖2.4所示。圖2.2離散周期序列x(n)圖2.3離散序列的離散傅里葉變換幅值譜圖2.4離散序列的離散傅里葉變換相位譜若序列長度N比較小，而頻率采樣間隔比較大，則不能直接觀測其頻率特性。為解決這一問題，通常采用原始序列xN(n)補零的方法，形成L(L≥N)點長的新離散序列xL(n)，即(2.39)其離散傅里葉變換為(2.40)

顯然，通過補零的方法，可滿足適合的頻率采樣間隔要求。

【例2.6】

對例2.5中的序列補零至長度為120，試繪制出x(n)及其離散傅里葉變換X(k)的幅值譜和相位譜。

MATLAB程序如下:

%MATLABPROGRAM2-6

clc;

N=90;

n=[0:N－1];

xn=sin(0.84*pi*n)+sin(0.88*pi*n);

%補零序列的產(chǎn)生及其傅里葉變換

L=120;

n1=[0:L－1];

xn=[xn，zeros(1，L－N)];

Xk=dft(xn，L);

magXk=abs(Xk);AngleXk=angle(Xk);

figure(1);

stem(n1，xn);

xlabel(′n′);ylabel(′x(n)′);title(′新序列x(n)′);

figure(2);

k=[0:length(magXk)－1];

stem(k，magXk);

xlabel(′k′);ylabel(′|X(k)|′);title(′DFT的幅值譜′);

figure(3);

stem(k，AngleXk);

xlabel(′k′);ylabel(′Angle(Xk)′);title(′DFT相位譜′);

程序運行結(jié)果如圖2.5～圖2.7所示。圖2.5補零新離散序列x(n)圖2.6新序列的離散傅里葉變換幅值譜圖2.7新序列的離散傅里葉變換相位譜補零可使數(shù)據(jù)N為2的整數(shù)倍，以便于使用快速傅里葉變換(FFT)，同時，還起到對原信號X(k)插值的作用。2.3.2DFT的性質(zhì)

1.線性性

如果x1(n)和x2(n)是兩個有限長序列，長度分別為N1和N2，假定

y(n)=ax1(n)+bx2(n)

式中，a和b為常數(shù)，y(n)的長度N=max[N1，N2]，則y(n)的N點DFT為

Y(k)=DFT[y(n)]=aX1(k)+bX2(k)，0≤k≤N－1

(2.41)

其中，X1(k)和X2(k)分別為x1(n)和x2(n)的N點DFT。

2.時域循環(huán)移位

設(shè)x(n)為長度為N的有限長序列，則x(n)的循環(huán)移位定義為

y(n)=x((n+m))NRN(N)

(2.42)

則

Y(k)=DFT[y(n)]=X(k)WN－mk

(2.43)

其中，X(k)=DFT[x(n)]，0≤k≤N－1。式(2.43)表明:序列在時域的循環(huán)移位將使頻域產(chǎn)生附加相移，幅頻特性保持不變。因此，循環(huán)移位又稱為圓周移位。序列向右移位也具有相同的性質(zhì)。

【例2.7】

用MATLAB編寫函數(shù)cirshift實現(xiàn)序列的循環(huán)移位的程序。

MATLAB程序如下:

%MATLABPROGRAM2-7

function y=cirshift(x，m，N);

%循環(huán)移位函數(shù)

%x:輸入序列

%m:移位位數(shù)

%N:序列長度

Iflength(x)>N

error(′Nmustbegreaterthanlength(x)′)

end

x=[x，zeros(1，length(x))];

n=[0:N－1];

n=sigmod((n－m)，N);

y=x(n+1);函數(shù)sigmod用來找出周期序列任意位置n所對應(yīng)的主值有限序列x(n)中的位置m。

function m=sigmod(n，N);

%用來找出周期序列任意位置n所對應(yīng)的主值有限序列x(n)中的位置

m=rem(n，N);

m=m+N;

m=rem(m，N);

3.頻域循環(huán)移位

若X(k)=DFT[x(n)]，0≤k≤N－1，且

Y(k)=X(k－l)NRN

則

y(n)=IDFT[Y(k)]=x(n)WN－ml

(2.44)

4.循環(huán)卷積

有限長序列x(n)和h(n)，長度分別為N1和N2，N=max[N1，N2]。x1(n)和x2(n)的N點DFT分別為X1(k)=DFT[x1(n)]，X2(k)=DFT[x2(n)]。

若

X(k)=X1(k)·X2(k)

則(2.45)稱式(2.45)表示的運算為x1(n)與x2(n)的循環(huán)卷積。循環(huán)卷積過程中，要求對x2(m)翻轉(zhuǎn)并循環(huán)移位。為區(qū)別于線性卷積，其運算符號用表示，即

(2.46)上式表明循環(huán)卷積滿足交換律。

【例2.8】

用MATLAB編寫函數(shù)circonvt實現(xiàn)序列的循環(huán)卷積。

MATLAB程序如下:

%MATLABPROGRAM2-8

function y=circonvt(x1，x2，N);

%循環(huán)卷積運算函數(shù)

%x1和x2為長度小于或等于N的有限長序列

%m:移位位數(shù)

%N:序列長度

iflength(x1)>N

error(′Nmustbegreaterthanlength(x1)′)

end

iflength(x2)>N

error(′Nmustbegreaterthanlength(x2)′)

end

x1=[x1，zeros(1，N－length(x1))];

x2=[x2，zeros(1，N－length(x2))];

xk1=dft(x1，N);

xk2=dft(x2，N);

Yk=xk1.*xk2;

y=idft(Yk，N);2.3.3DFT的應(yīng)用

DFT的快速算法FFT的出現(xiàn)，使DFT在數(shù)字通信、語音信號處理、圖像處理、功率譜估計、仿真、系統(tǒng)分析、雷達理論、光學(xué)以及數(shù)值分析等各個領(lǐng)域都得到了廣泛應(yīng)用。

1.用DFT計算線性卷積

在實際應(yīng)用中，為了分析時域離散線性系統(tǒng)或者對序列進行濾波處理等，需要計算兩個序列的線性卷積，為了提高運算速度，也希望用DFT。但是，DFT只能直接用來計算循環(huán)卷積。因此，必須導(dǎo)出線性卷積和循環(huán)卷積之間的關(guān)系以及循環(huán)卷積與線性卷積相等的條件。假設(shè)x1(n)和x2(n)都是有限長序列，長度分別是N和M。它們的線性卷積和循環(huán)卷積分別表示為(2.47)(2.48)其中，L≥max[N，M]。用DFT的循環(huán)卷積計算線性卷積的條件是將兩序列各自加零延長至L點，求其L點的循環(huán)卷積。當(dāng)L<(N+M－1)時，這種延拓將會產(chǎn)生混疊現(xiàn)象，其循環(huán)卷積與線性卷積不相等;當(dāng)L≥(N+M－1)時，其循環(huán)卷積與線性卷積相等，因而可用此種情況下的循環(huán)卷積來計算線性卷積。圖2.8為計算線性卷積框圖。圖2.8用DFT的循環(huán)卷積計算線性卷積

【例2.9】

已知序列x1(n)=2cos(2πn/5)，0≤n≤9;x2(n)=2e－0.6n，0≤n≤19。試計算x1(n)和x2(n)的線性卷積。

解

L=(N+M－1)=10+20－1=29

yl(n)=x1(n)*x2(n)=x1(n)x2(n)=IDFT(X1(k)X2(k))

MATLAB程序如下:

%MATLABPROGRAM2-9

n1=[0:9];

n2=[0:19];

x1=2*cos(2*pi*n1/5);

x2=2*exp(－0.6*n2);

L=length(x1)+length(x2)－1;

y1=circonvt(x1，x2，L);%利用循環(huán)卷積計算線性卷積

subplot(2，1，1);

n=[0:L－1];

stem(n，abs(y1));

xlabel(′時間序號n′);ylabel(′y(n)振幅′);title(′循環(huán)卷積計算y(n)′);

y2=conv(x1，x2);%利用線性卷積函數(shù)conv計算

subplot(2，1，2);

stem(n，abs(y2));

xlabel(′時間序號n′);ylabel(′y(n)振幅′);title(′線性卷積計算y(n)′);

程序運行結(jié)果如圖2.9所示。圖2.9不同算法的線性卷積計算結(jié)果

2.用DFT實現(xiàn)線性時不變系統(tǒng)

線性時不變系統(tǒng)的輸出響應(yīng)y(n)是輸入和單位采樣響應(yīng)的線性卷積，通過計算循環(huán)卷積的方法可以實現(xiàn)線性卷積，因此可用DFT實現(xiàn)線性時不變系統(tǒng)。

y(n)=x(n)*h(n)=x(n)h(n)=IDFT(X(k)H(k))

(2.49)

2.4快速傅里葉變換(FFT)

DFT是信號分析與處理中的一種重要變換。當(dāng)N較大時，直接計算DFT的計算量太大。在快速傅里葉變換(FFT)出現(xiàn)以前，直接用DFT算法進行譜分析和信號的實時處理是不切實際的。對實時性很強的信號處理來說，要求計算速度快，因此，需要改進DFT的計算方法，以大大減少運算次數(shù)。FFT作為DFT的一種快速算法，是數(shù)字信號處理領(lǐng)域中的一項重大突破。2.4.1FFT的基本思想

FFT不是一種與DFT不同的新變換，而是DFT的一種快速計算方法。N點序列x(n)的DFT變換對定義為(2.50)(2.51)式中，WN=稱為旋轉(zhuǎn)因子?？紤]x(n)為復(fù)數(shù)序列的一般情況，求N點X(k)需要N2次復(fù)數(shù)乘法，N(N－1)次復(fù)數(shù)加法，而1次復(fù)數(shù)乘法需要4次實數(shù)乘法和2次實數(shù)加法，當(dāng)N很大時，運算量很大。為提高計算機處理速度，利用旋轉(zhuǎn)因子WN的特性來減少運算量。WN具有明顯的周期性和對稱性。

(1)WN的周期性:

(2)WN的對稱性:(2.52)(2.53)經(jīng)過周期性與對稱性簡化之后，容易發(fā)現(xiàn)DFT運算中存在著不必要的重復(fù)計算，F(xiàn)FT避免了這種重復(fù)，大大簡化了運算過程，這就是FFT的基本思想。2.4.2FFT算法及其實現(xiàn)

FFT算法有多種形式，但基本上分為兩大類:時間抽取法FFT(DecimationInTimeFFT，DITFFT)和頻域抽取法FFT(DecimationInFrequencyFFT，DIFFFT)。

FFT算法有兩個發(fā)展方向:一是針對N等于2的整數(shù)次冪的算法，如基2算法、基4算法、實因子算法和分裂基算法等;二是N不等于2的整數(shù)次冪的算法，如Winograd算法及素因子算法等。本節(jié)重點討論FFT的基2算法和分裂基算法。

1.時間抽取基2FFT算法

對于長度為N＝2M(M為正整數(shù))的DFT運算，可按奇偶分解為

(2.54)由此，得到X(k)的前N/2點的值(2.55)對于X(k)的后N/2點的值(2.56)式(2.55)和式(2.56)又稱為蝶形運算關(guān)系式(由于運算流程看起來像是一只蝴蝶，故稱這種運算為蝶形運算)。由于N/2是偶數(shù)，可以類似地將N/2點序列再分解為兩個N/4點序列，直至最后是2點序列的DFT。由于每次分解均將序列從時域中按奇偶提取，且每次都是一分為二的，故稱為時間抽取基2-FFT算法。

綜上所述，只要求出兩個N/2點的DFT，即X1(k)和X2(k)，再經(jīng)過蝶形運算就可求出全部X(k)的值，運算量大大減少了。圖2.10和圖2.11給出了N=8點DFT的時間抽取算法的流程分解圖。圖2.10時間抽取基2-FFT算法的第一次分解圖(N=8)圖2.11時間抽取基2-FFT算法的第二次分解圖(N=8)

2.頻域抽取法FFT

在基2快速算法中，頻域抽取法FFT(DIF-FFT)也是一種常用的快速算法。

設(shè)序列x(n)長度為N＝2M(M為正整數(shù))，首先將x(n)按n的順序分成前后兩半，再把輸出X(k)按序號k的奇偶分組，其DFT可表示為式中，WNNk/2=(－1)k。分別令k=2r和k=2r+1，而r=0,1,…,N/2，于是得

(2.57)(2.58)令(2.59)(2.60)則(2.61)(2.62)這樣，就將N點DFT分解成了兩個N/2點的DFT。同

理，按此方法繼續(xù)將每個N/2點DFT的輸出再分解為偶數(shù)組與奇數(shù)組，直到得到兩點的DFT。圖2.12和圖2.13以N=8點DFT為例，給出了頻域抽取算法的分解流程。

將N點DFT分解為兩個N/2點DFT的組合(N=8)。圖2.12N點DFT分解為兩個N/2點DFT的組合圖2.13N/2點DFT分解為兩個N/4點DFT的組合頻域抽取法與時間抽取法的異同:

(1)頻域抽取法輸入是自然順序的，輸出是倒位序的;時間抽取法正好相反。

(2)頻域抽取法的基本蝶形與時間抽取法的基本蝶形有所不同。

(3)頻域抽取法運算量與時間抽取法相同。

(4)頻域抽取法與時間抽取法的基本蝶形是互為轉(zhuǎn)置的。

3.分裂基FFT算法

分裂基(split-radix)算法又稱基2/4算法或混合基算法，它既和基2算法有關(guān)，也和基4算法有關(guān)。當(dāng)序列長度N=4M(M為正整數(shù))時，基4算法比基2算法更有效，它將N點序列抽取成4個子序列(x(4n)，x(4n+1)，x(4n+2)，x(4n+3)，n=0，1，…，N/4－1)，再對4個序列進行DFT的蝶形運算。

分裂基算法的基本思路是對偶序號輸出使用基2算法，對奇序號輸出使用基4算法。此算法在針對N=2M的算法中具有最少的乘法次數(shù)和加法次數(shù)，并且和Cooley-Tukey算法有同樣好的結(jié)構(gòu)，因此被認為是最好的FFT算法。對于N=2M點DFT，其偶序號輸出項按基2算法，即(2.63)

對奇序號輸出項采用基4算法，即(2.64)(2.65)式中，r=0，1，…，N/4－1。上面三式構(gòu)成了分裂基算法的L型算法結(jié)構(gòu)。在MATLAB的信號處理工具箱中，內(nèi)部函數(shù)fft和ifft用于實現(xiàn)快速傅里葉變換和其逆變換。函數(shù)fft用于快速傅里葉變換，調(diào)用方式為

y=fft(x);

y=fft(x，N);

y=fft(x)利用FFT算法計算序列x的離散傅里葉變換。當(dāng)x為矩陣時，y為矩陣x每一列的FFT。當(dāng)x長度為2的整數(shù)次冪時，函數(shù)fft采用高速的基2FFT算法，否則采用混合基算法。

y=fft(x，N)采用N點FFT。當(dāng)序列x長度小于N時，函數(shù)fft自動對序列尾部補零，構(gòu)成N點數(shù)據(jù);當(dāng)x長度大于N時，函數(shù)fft自動截取序列前面N點數(shù)據(jù)進行FFT。

函數(shù)ifft用于快速傅里葉反變換，調(diào)用方式為

y=ifft(x);

y=ifft(x，N);式中，x為離散序列的傅里葉變換X(k)，y是IDFT[X(k)]，N為正整數(shù)。若x為向量且長度小于N，則函數(shù)將x補零至N;若向量x的長度大于N，則函數(shù)截斷x，使之長度等于N。若x為矩陣，則按同樣的方法對x進行處理。

【例2.10】

若已知序列x(n)=[3,6,8,4,5,2,1,9]，求DFT[x(n)]。

MATLAB程序如下:

%MATLABPROGRAM2-10

N=8;

n=[0:N-1];

xn=[3，6，8，4，5，2，1，9];

Xk=fft(xn);

程序運行結(jié)果如下:

Xk=

38.00004.3640－6.2929i－1.0000+5.0000i－8.3640+7.7071i

－4.0000－8.3640－7.7071i－1.0000－5.0000i

4.3640+6.2929i2.4.3FFT的應(yīng)用

FFT算法在眾多領(lǐng)域中得到了廣泛的應(yīng)用，包括譜分析、線性濾波、相關(guān)計算等。

1.用FFT進行譜分析

經(jīng)函數(shù)fft求得的序列一般是復(fù)序列，通常會對此進行頻譜分析，而頻譜分析涉及到幅值和相位。MATLAB提供了求復(fù)數(shù)幅值和相位的函數(shù)，即abs(x)和angle(x)。函數(shù)abs(x)用于計算復(fù)向量x的幅值;angle(x)用于計算復(fù)向量x的相角(介于－π和π之間，以弧度表示)。

【例2.11】

設(shè)模擬信號x(t)=2sin4πt+5cos8πt，以t=0.01n(n=0，1，…，N－1)進行采樣，試用fft函數(shù)對其進行頻譜分析。N分別為(1)N=45;(2)N=50;(3)N=55;(4)N=60。MATLAB程序如下:

%MATLABPROGRAM2-11

%計算N=45的FFT并繪出其幅頻曲線

N=45;n=[0:N－1];t=0.01*n;

q=n*2*pi/N;

x=2*sin(4*pi*t)+5*cos(8*pi*t);

y=fft(x，N);

subplot(2，2，1);

stem(q，abs(y));

title(′FFTN=45幅頻曲線′);

%計算N=50的FFT并繪出其幅頻曲線

N=50;n=[0:N－1];t=0.01*n;

q=n*2*pi/N;

x=2*sin(4*pi*t)+5*cos(8*pi*t);

y=fft(x，N);

subplot(2，2，2);

stem(q，abs(y));

title(′FFTN=50幅頻曲線′);

%計算N=55的FFT并繪出其幅頻曲線

N=55;n=[0:N-1];t=0.01*n;

q=n*2*pi/N;

x=2*sin(4*pi*t)+5*cos(8*pi*t);

y=fft(x，N);

subplot(2，2，3);

stem(q，abs(y));

title(′FFTN=55幅頻曲線′);

%計算N=60的FFT并繪出其幅頻曲線

N=60;n=[0:N-1];t=0.01*n;

q=n*2*pi/N;

x=2*sin(4*pi*t)+5*cos(8*pi*t);

y=fft(x，N);

subplot(2，2，4);

stem(q，abs(y));

title(′FFTN=60幅頻曲線′);

程序運行結(jié)果如圖2.14所示。圖2.14不同N值下的采樣序列的DFT的幅頻曲線由t=0.01n進行取樣可得，采樣頻率fs=100Hz。連續(xù)信號的最高模擬角頻率Ω=8π，由Ω=2πf可得，最高頻率為8π/2π=4Hz，滿足采樣定理的要求。采樣序列為

2.用FFT實現(xiàn)線性卷積

在MATLAB中，當(dāng)N<50時，可以用conv函數(shù)計算卷積;當(dāng)N較大時，用FFT法計算更有效。設(shè)離散序列x1(n)和x2(n)的長度分別為N1和N2，計算二者的卷積為

x1(n)*x2(n)=IFFT[FFT(x1(n))·FFT(x2(n))]

【例2.12】

已知序列x(n)=R4(n)，

(1)用conv函數(shù)求x(n)與x(n)的線性卷積y(n),并繪出圖形;

(2)用FFT求x(n)與x(n)的4點循環(huán)卷積y1(n),并繪出圖形;

(3)用FFT求x(n)與x(n)的8點循環(huán)卷積y2(n),并繪出圖形。

MATLAB程序如下:

%MATLABPROGRAM2-12

N1=4;N2=8;

n1=[0:1:N1－1]; n2=[0:1:N2－1];

x=[1，1，1，1];

%構(gòu)造序列x(n)

x1=[1，1，1，1，0，0，0，0];%在序列x(n)后補4個零

subplot(2，2，1);

stem(n1，x);

title(′序列x(n)′);

grid;

y1=conv(x，x);

%y1為x(n)與x(n)的線性卷積

subplot(2，2，2);

stem(0:1:length(y1)－1，y1);

title(′x(n)與x(n)線性卷積′);

grid;

X2=fft(x);

%計算x(n)與x(n)的4點循環(huán)卷積

Y2=X2.*X2;

y2=ifft(Y2);

subplot(2，2，3);

stem(n1，y2);

title(′x(n)與x(n)的4點循環(huán)卷積′);

grid;

X3=fft(x1);

%計算x(n)與x(n)的8點循環(huán)卷積

Y3=X3.*X3;

y3=ifft(Y3);

subplot(2，2，4);

stem(n2，y3);

title(′x(n)與x(n)的8點循環(huán)卷積′);

grid;

程序運行結(jié)果如圖2.15所示。圖2.15用FFT計算線性卷積

3.利用FFT計算線性相關(guān)

由第1章介紹的相關(guān)函數(shù)和卷積的關(guān)系可知，利用FFT可以快速計算線性相關(guān)。設(shè)離散序列x1(n)和x2(n)的長度分別為N1和N2，計算二者相關(guān)函數(shù)的步驟如下:

(1)將離散序列x1(n)和x2(n)通過增加零值延長到N=N1+N2－1;

(2)利用FFT計算x1(n)和x2(n)各自的N點DFT;

(3)計算R(k)=X1(k)X*2(k)，其中，X*2(k)為X2(k)的共軛;

(4)利用IFFT計算R(k)的傅里葉反變換，即可得到二者的相關(guān)函數(shù)r(n)。2.5Z變換

利用差分方程可求離散系統(tǒng)的結(jié)構(gòu)及瞬態(tài)解，為了分析系統(tǒng)的另外一些重要特性，如穩(wěn)定性和頻率響應(yīng)等，需要研究離散時間系統(tǒng)的Z變換，它是分析離散系統(tǒng)和離散信號的重要工具。

2.5.1Z變換的定義

一個離散序列x(n)的Z變換定義為(2.66)式中，z為復(fù)變量，以其實部為橫坐標、虛部為縱坐標構(gòu)成的平面為z平面。這種變換也稱為雙邊Z變換，與此相應(yīng)的還有單邊Z變換。單邊Z變換只是對單邊序列(n≥0)進行變換的Z變換，定義為(2.67)單邊Z變換只在少數(shù)情況下與雙邊Z變換有所區(qū)別，大多數(shù)情況下，可以把單邊Z變換看成是雙邊Z變換的一種特例，即單邊Z變換是因果序列情況下的雙邊Z變換。這種單邊Z變換的求和是從零到無限大，因此，對于因果序列，用兩種Z變換定義計算出的結(jié)果是一樣的。本書中如不特別說明，均是用雙邊Z變換對信號進行的分析和變換。

Z變換存在的條件是式(2.66)或式(2.67)等號右邊級數(shù)收斂，即要求級數(shù)絕對可和:(2.68)

Z變換與DTFT的關(guān)系可表示為(2.69)

LTI系統(tǒng)傳遞函數(shù)所涉及的Z變換都為z的有理函數(shù)，因此可表示為(2.70)或(2.71)式中，z=ξl稱為零點，z=λl稱為極點。當(dāng)N>M時，上式在z=0處有額外N－M個零點，當(dāng)N<M時，有額外M－N個極點。2.5.2Z變換的收斂域

一般序列的Z變換并不一定對任何z值都存在，z平面上使式(2.68)左邊級數(shù)收斂的區(qū)域稱為Z變換的收斂域(ROC)。級數(shù)一致收斂的條件是絕對值可和。但|z|值在一定范圍內(nèi)才能滿足絕對可和的條件，這個范圍一般表示為

Rx－≤|z|≤Rx+

(2.72)

這就是收斂域，即一個以Rx－和Rx+為半徑的兩個圓所圍成的環(huán)形區(qū)域。Rx－和Rx+稱為收斂半徑，Rx－和Rx+的大小，即收斂域的位置與具體序列有關(guān)，特殊情況為Rx－等于0，Rx+為無窮大，這時圓環(huán)變成圓或空心圓。

關(guān)于Z變換的收斂域，主要介紹以下四種序列。

1.有限長序列序列，其Z變換為

X(z)是有限項的級數(shù)和，若級數(shù)每一項有界，則有限項和也有界，所以有限長序列Z變換的收斂域為0<|z|<∞的z平面。

2.右邊序列

右邊序列只在n≥n1時有值，在n<n1時，x(n)=0，其Z變換為(2.74)收斂域為|z|>Rx－，即收斂半徑Rx－以外的z平面。右邊序列中最重要的一種序列是因果序列。因果序列只在n≥0時有值，在n<0時，x(n)=0，即n1≥0的右邊序列。

3.左邊序列

序列x(n)只在n≤n2時有值，在n>n2時，x(n)=0，其Z變換為

(2.75)收斂域為|z|<Rx+，即在收斂半徑為Rx+的圓內(nèi)。

【例2.13】

求x(n)=anu(－n－1)的Z變換及其收斂域。

解X(z)存在要求|a－1z|<1，即收斂域為|z|<|a|。因此，x(n)的Z變換為

4.雙邊序列

雙邊序列可看做一個左邊序列和一個右邊序列之和。因此，雙邊序列Z變換的收斂域是這兩個序列Z變換收斂域的公共部分，其Z變換表示為(2.76)

如果Rx+>Rx－，則存在公共的收斂區(qū)間，X(z)的收斂域為Rx－<|z|<Rx+。如果Rx+<Rx－，則不存在公共的收斂區(qū)間，X(z)無收斂域，不存在Z變換。由此可知，Z變換的收斂域具有如下特點:

(1)收斂域是一個圓環(huán)，有時可向內(nèi)收縮到原點，有時可向外擴展到∞，只有x(n)=δ(n)的收斂域是整個z平面。

(2)在收斂域內(nèi)沒有極點，X(z)在收斂域內(nèi)的每一點上都是解析函數(shù)。

【例2.14】

求x(n)=c|n|的Z變換。

解

若|c|<1，則存在公共的收斂區(qū)域若|c|>1，則無公共的收斂區(qū)域，此時，x(n)向兩邊發(fā)散。

下面列出幾種常見的序列及其Z變換。

(1)單位沖激序列:

(2)單位階躍序列:

(3)單位斜變序列:

(4)單位指數(shù)序列:(5)單邊正余弦序列:2.5.3Z變換的性質(zhì)

Z變換的許多重要性質(zhì)在數(shù)字信號處理中常常要用到。

1.線性性

設(shè)X(z)=Z[x(n)]，收斂域為Rx－<|z|<Rx+;Y(z)=Z[y(n)]，

收斂域為Ry－<|z|<Ry+。則

Z(ax(n)+by(n))=aX(z)+bY(z)

收斂域為

max(Rx－，Ry－)<|z|<min(Rx+，Ry+)

(2.77)

2.序列的時移

設(shè)X(z)=Z[x(n)]，收斂域為Rx－<|z|<Rx+。則

Z[x(n－n0)]=z－n0X(z)，收斂域為Rx－<|z|<Rx+

(2.78)

3.序列指數(shù)加權(quán)(z域尺度變換)

設(shè)X(z)=Z[x(n)]，收斂域為Rx－<|z|<Rx+。則

Z[anx(n)]=X(z/a)，收斂域為|a|Rx－<|z|<|a|Rx+

(2.79)

4.序列的線性加權(quán)(z域微分)

設(shè)X(z)=Z[x(n)]，收斂域為Rx－<|z|<Rx+。則

(2.80)

(2.81)

收斂域唯一可能的變化是加上或去掉零或無窮。

5.初值定理

X(z)是因果序列x(n)的Z變換，則。，收斂域為Rx－<|z|<Rx+，收斂域為Rx－<|z|<Rx+

6.終值定理

X(z)是因果序列x(n)的Z變換，則

終值定理只有在極限存在時才能用，此時，X(z)的極點必須在單位圓內(nèi)(若位于單位圓上，則只能位于z=1，且是一階極點)。

7.時域卷積定理

(2.82)

卷積的Z變換的收斂域至少是原序列Z變換的收斂域的交集。當(dāng)出現(xiàn)零、極點相抵時，收斂域可能會擴大。卷積性質(zhì)是Z變換的重要性質(zhì)，使用Z變換可計算兩信號的卷積。

8.z域卷積定理

設(shè)X(z)=Z[x(n)]，收斂域為Rx－<|z|<Rx+;Y(z)=Z[y(n)]，收斂域為Ry－<|z|<Ry+。則(2.83)其中，C是圍繞原點的閉合曲線，位于收斂域重疊部分內(nèi)。收斂域為Rx－Ry－<|z|<Rx+Ry+。2.5.4逆Z變換

已知函數(shù)X(z)及其收斂域，反過來求序列x(n)的變換稱為逆Z變換，常用Z－1[X(z)]表示。若則逆Z變換為(2.84)

逆Z變換是一個對X(z)zn－1進行的圍線積分，積分路徑C是一條在X(z)收斂環(huán)域(Rx－，Rx+)以內(nèi)反時針方向繞原點一周的單圍線。直接計算圍線積分(留數(shù)法)比較麻煩，一般不采用此法求逆Z變換。求解逆Z變換的常用方法有:冪級數(shù)法(長除法)和部分分式法。冪級數(shù)法是按Z變換的定義，將X(z)寫成冪級數(shù)形式，級數(shù)的系數(shù)就是序列x(n)。如果x(n)是右邊序列，則級數(shù)應(yīng)是負冪級數(shù);如果x(n)是左邊序列，則級數(shù)是正冪級數(shù)。

對于大多數(shù)單階極點的序列，常用部分分式展開法求逆Z變換。設(shè)x(n)的Z變換X(z)是有理函數(shù)，分母多項式是N階，分子多項式是M階，將X(z)展開成一些簡單常用的部分分式之和，通過查表求得各部分的逆變換，再相加即得到原序列x(n)。

設(shè)X(z)只有N個一階極點，則可展開成

(2.85)或?qū)懗捎^察以上兩式，X(z)/z在z=0的極點留數(shù)就是系數(shù)A0，在z=zm的極點留數(shù)就是系數(shù)Am。即(2.86)

求出Am系數(shù)(m=0，1，2，…，N)后，很容易求得x(n)序列。在MATLAB中，提供了函數(shù)impz和函數(shù)residuez來求得逆Z變換。

1.函數(shù)impz

函數(shù)impz提供了時域序列的樣本，該序列假定為因果，其三種調(diào)用格式為

[h，t]=impz(num，den);

[h，t]=impz(num，den，L);

[h，t]=impz(num，den，L，F(xiàn)T);

其中，num和den是按z－1的升冪排列的分子和分母多項式系數(shù)的行向量;L是所求逆變換的樣本數(shù);h是包含從樣本n=0開始的逆變換的樣本向量;t是h的長度;FT是單位為Hz的給定抽樣頻率，默認值為1。

【例2.15】

計算，求|z|>1的逆Z變換。

解有理分式X(z)的分子和分母多項式都按z的降冪排列為MATLAB程序如下:

%MATLABPROGRAM2-15

num=[0，1，1];

den=[3，－4，1];

[x，t]=impz(num，den);

disp(′x(n)樣本序號′);disp(t′);%顯示輸出參數(shù)

disp(′x(n)樣本向量′);disp(x′);

程序運行結(jié)果為

x(n)樣本長度:

012345678910

x(n)樣本向量:

00.33330.77780.92590.97530.99180.99730.99910.99970.99991.0000

2.函數(shù)residuez

函數(shù)residuez適合計算離散系統(tǒng)有理函數(shù)的留數(shù)和極點，可以用于求解序列的逆Z變換。

函數(shù)residuez的基本調(diào)用格式為

[r，p，c]=residuez(b，a);

其中，輸入?yún)?shù)b=[b0，b1，…，bM]為分子多項式的系數(shù);a=[a0，a1，…，aN]為分母多項式的系數(shù)，這些多項式都按z的降冪排列。

輸出參數(shù)r是極點的留數(shù)，p是極點，c是無窮項多項式的系數(shù)項，僅當(dāng)M≥N時存在。用MATLAB的函數(shù)residuez將X(z)分解為式(2.70)或式(2.71)的形式，它們包含Z變換的基本形式，由Z變換的基本形式求出x(n)，即(2.87)

【例2.16】

計算

的逆Z變換。

解有理分式X(z)的分子和分母多項式都按z的降冪排列為MATLAB程序如下:

%MATLABPROGRAM2-16

b=[0，1];a=[2，－3，1];%多項式的系數(shù)

[r，p，c]=residuez(b，a);

%求留數(shù)、極點和系數(shù)項

disp(′留數(shù)r:′);disp(r′);%顯示輸出參數(shù)

disp(′極點p:′);disp(p′);

disp(′系數(shù)項c:′);disp(c′);

程序運行結(jié)果為

留數(shù):1－1

極點:1.00000.5000

系數(shù)項:

根據(jù)程序運行結(jié)果，得

(1)當(dāng)1<|z|<∞時，Rx－=1，兩個極點z1=1，，|z1|≤Rx－,|z2|<Rx－。(2)當(dāng)0<|z|<1/2時，Rx+=1/2，|z1|>Rx+，|z2|≥Rx+(3)當(dāng)1/2<|z|<1時，Rx－=1/2，Rx+=1，|z1|≥Rx+，|z2|≤Rx－2.5.5利用Z變換求解差分方程

求解差分方程就是求系統(tǒng)響應(yīng)y(n)。在第1章中介紹了差分方程的遞推解法，下面介紹Z變換解法。這種方法將差分方程變成了代數(shù)方程，使求解過程變得簡單了。

設(shè)N階線性常系數(shù)差分方程為

(1)求穩(wěn)態(tài)解:若輸入序列x(n)的初始狀態(tài)為零，對式(2.88)求Z變換，再求逆Z變換，則求得穩(wěn)態(tài)解y1(n)。式中，H(z)稱為系統(tǒng)函數(shù)，且具有式(2.70)的形式，即

y1(n)=Z－1[Y(z)]

(2)求暫態(tài)解:對于N階差分方程，求暫態(tài)解必須已知N個初始條件。若系統(tǒng)的輸入為零，則計算由初始狀態(tài)引起的響應(yīng)y2(n)。

(3)求全解:

y(n)=y1(n)+y2(n)在MATLAB中，可以利用函數(shù)filter和filtic求解差分方程的全解，即離散時間系統(tǒng)在輸入和初始狀態(tài)作用下的響應(yīng)。關(guān)于函數(shù)filter和filtic具體見第3章。

【例2.17】

求解系統(tǒng)差分方程y(n)=2x(n)－3x(n－1)+7x(n－4)，x(n)=0.8nu(n)。

解方程兩邊取Z變換:

H(z)=2－3z－1+7z－4

MATLAB程序如下:

%MATLABPROGRAM2-17

num=[2，－3，0，0，7];

N=30;

n=[0:N－1];

x=0.8.^n;

y=filter(num，1，x);

stem(n，y);

grid;

程序運行結(jié)果如圖2.16所示。圖2.16差分方程的輸出序列

【例2.18】