數(shù)字信號(hào)處理實(shí)驗(yàn)教材正文

1、數(shù)字信號(hào)處理實(shí)驗(yàn)指導(dǎo)書(shū)實(shí)驗(yàn)一 離散時(shí)間信號(hào)與系統(tǒng)的實(shí)現(xiàn)與表示由于MATLAB數(shù)值計(jì)算的特點(diǎn)，用它來(lái)分析離散時(shí)間信號(hào)與系統(tǒng)是很方便的。一、離散時(shí)間信號(hào)在MATLAB中，可以用一個(gè)列向量來(lái)表示一個(gè)有限長(zhǎng)度的序列，例如：在MATLAB中可以表示為：x = 2 1 0 2 3 -1 2 3;但是這種表示方法沒(méi)有包含抽樣的位置信息，要完全表示，要用和兩個(gè)向量表示，例如：（箭頭表示第0個(gè)抽樣點(diǎn)的位置）在MATLAB中可以表示為： n=-4 -3 -2 -1 0 1 2 3; x = 2 1 0 2 3 -1 2 3;由于內(nèi)存有限，MATLAB無(wú)法表示無(wú)限序列。下面介紹幾種典型離散信號(hào)的表示方法。1、單位脈

2、沖序列（1）在MATLAB中zeros(1,N)函數(shù)來(lái)實(shí)現(xiàn)有限區(qū)間的，例如：x= zeros(1,N);x(1)=1;（2）我們還可以通過(guò)編寫(xiě)impseq函數(shù)來(lái)實(shí)現(xiàn)，該函數(shù)代碼如下：function x,n= impseq(n0,n1,n2)%產(chǎn)生單位脈沖序列x(n)=(n-n0)%x,n= impseq(n0,n1,n2)% 1,n=n0% x =% 1,nn0%n=時(shí)間序列n=n1:n2;x=(n-n0)=0該函數(shù)產(chǎn)生如下信號(hào)，取值區(qū)間為，得到序列例1：輸入命令：x,n= impseq(0,-5,5)得到如下結(jié)果：x = 0 0 0 0 0 1 0 0 0 0 0n = -5 -4 -3

3、-2 -1 0 1 2 3 4 5 2、單位階躍序列（1）在MATLAB中ones(1,N)函數(shù)來(lái)實(shí)現(xiàn)有限區(qū)間的，例如：x= ones(1,N);（2）我們還可以通過(guò)編寫(xiě)impseq函數(shù)來(lái)實(shí)現(xiàn)，該函數(shù)代碼如下：function x,n= impseq(n0,n1,n2)%產(chǎn)生單位階躍序列x(n)=u(n-n0);n1n0n2%x,n= impseq(n0,n1,n2)% 1,nn0% x =% 1,n=0該函數(shù)產(chǎn)生如下信號(hào)，取值區(qū)間為，得到序列例2：輸入命令：x,n= impseq(0,-5,5)得到如下結(jié)果：x = 0 0 0 0 0 1 1 1 1 1 1n = -5 -4 -3 -2 -

4、1 0 1 2 3 4 53、正弦序列sin函數(shù)就可以產(chǎn)生正弦波，例如：n= 0:N-1;x=A*sin(n)4、實(shí)指數(shù)序列MATLAB實(shí)現(xiàn)：n= 0:N-1;x=a.n5、復(fù)指數(shù)序列MATLAB實(shí)現(xiàn)：n= 0:N-1;x=exp(lu+j*w)*n);6、隨機(jī)序列rand(1,N)：產(chǎn)生0,1上均勻分布的隨機(jī)序列；randn(1,N)：產(chǎn)生均值為0，方差為1的高斯隨機(jī)序列，即白噪聲序列。二、離散時(shí)間系統(tǒng)離散時(shí)間系統(tǒng)的階差分方程為：在Matlab中，提供函數(shù)filter，用來(lái)求解給定輸入和差分方程系數(shù)時(shí)，差分方程的數(shù)值解。調(diào)用格式為：Y=filter(B,A,X)其中：調(diào)用參數(shù)，分別為和各項(xiàng)差

5、分項(xiàng)系數(shù)；為輸入信號(hào)序列，注意必須保證系數(shù)不為零。輸出序列與輸入序列具有相同的長(zhǎng)度。返回參數(shù)為差分方程的數(shù)值解例3：已知某線性時(shí)不變系統(tǒng)的差分方程如下：當(dāng)輸入序列為時(shí)的輸出結(jié)果。 解：MATLAB程序如下： N=41;a=0.8 -0.44 0.36 0.22; %輸入序列X各項(xiàng)差分項(xiàng)系數(shù)b=1 0.7 -0.45 -0.6; %輸出序列Y各項(xiàng)差分項(xiàng)系數(shù)x=1 zeros(1,N-1); %單位脈沖序列n=0:1:N-1; %時(shí)間序列y=filter(a,b,x); %輸入單位脈沖序列得到單位脈沖響應(yīng)stem(n,y);xlabel(n);ylabel(幅度);運(yùn)行結(jié)果如下：三、實(shí)驗(yàn)內(nèi)容1、編

6、寫(xiě)程序用來(lái)產(chǎn)生下列基本序列，并繪制它們的波形圖。（1）單位脈沖序列，起點(diǎn)為0，終點(diǎn)為10，在處有一單位脈沖；（2）單位階躍序列，起點(diǎn)為0，終點(diǎn)為10，在前為0，后為1；（3）復(fù)指數(shù)序列，2、已知系統(tǒng)的差分方程為：求系統(tǒng)的單位脈沖響應(yīng)并顯示其波形圖（選?。?shí)驗(yàn)二 離散卷積和一、卷積和的定義信號(hào)的卷積是針對(duì)時(shí)域信號(hào)處理的一種分析方法。信號(hào)的卷積一般用于求取信號(hào)通過(guò)某系統(tǒng)后的響應(yīng)。在信號(hào)與系統(tǒng)中，我們通常求取某系統(tǒng)的單位沖激響應(yīng)，所求得的可作為系統(tǒng)的時(shí)域表征。任意系統(tǒng)的系統(tǒng)響應(yīng)可用卷積的方法求得：圓周卷積是數(shù)字信號(hào)處理的重要內(nèi)容，由于它也是卷積的一種，必須了解其定義，圓周卷積的條件為卷積長(zhǎng)度必須大于

7、每個(gè)序列長(zhǎng)度和。二、卷積工具箱函數(shù)序列的卷積和運(yùn)算， MTLAB提供的計(jì)算兩個(gè)離散序列卷積的函數(shù)conv，其調(diào)用方式為：y=conv(f1,f2)其中：調(diào)用參數(shù)為卷積運(yùn)算所需的兩個(gè)序列，返回值是卷積結(jié)果。MATLAB conv( )函數(shù)不需要給定的非零樣值的時(shí)間序號(hào)，也不返回卷積和序列的時(shí)間序號(hào)。即的返回參數(shù)中只有卷積的結(jié)果，沒(méi)有的取值范圍。此外，conv( )假定都是從開(kāi)始，這就限制了它的應(yīng)用范圍（由離散序列卷積的性質(zhì)可知，當(dāng)序列和的起始點(diǎn)都為時(shí)，的取值范圍為至length(x)+length(h)-2）。因此，要對(duì)從任意值開(kāi)始的序列進(jìn)行卷積和運(yùn)算，同時(shí)正確標(biāo)識(shí)出函數(shù)conv( )的計(jì)算結(jié)果

8、各量，還須構(gòu)造序列和的對(duì)應(yīng)序號(hào)向量。例1：已知序列：，求其卷積和。解：MATLAB程序如下：n1=-2:2;f1=ones(1,length(k1);n2=0:3;f2=2.n2;f=conv(f1,f2)運(yùn)行結(jié)果為：f = 1 3 7 15 15 14 12 8由于用conv函數(shù)求出卷積后沒(méi)有給出新序列所對(duì)應(yīng)的時(shí)間變量。為此，我們可以構(gòu)造一個(gè)實(shí)用點(diǎn)函數(shù)dconv( )，該函數(shù)實(shí)現(xiàn)的功能為：它在完成conv函數(shù)功能的同時(shí)，還產(chǎn)生一個(gè)對(duì)應(yīng)新序列的時(shí)間變量，實(shí)現(xiàn)了序號(hào)向量的返回，并且可以繪制出序列和時(shí)域波形圖。dconv函數(shù)MATLAB程序如下所示：function f,n=dconv(f1,f2

9、,n1,n2)%the function of compute f=f1*f2%f: 卷積序列f(n)對(duì)應(yīng)的非零樣值向量%n: 序列f(n)的對(duì)因序號(hào)向量%f1: 序列f1(n)非零樣值向量%f2: 序列f2(n)非零樣值向量%n1: 序列f1(n)的對(duì)應(yīng)序號(hào)向量%n2: 序列f2(n)的對(duì)應(yīng)序號(hào)向量f=conv(f1,f2) %計(jì)算序列f1與f2的卷積和fn0=n1(1)+n2(1); %計(jì)算序列f非零樣值的起點(diǎn)位置n3=length(f1)+length(f2)-2; %計(jì)算卷積和f的非零樣值的寬度n=n0:n0+n3 %確定卷積和f非零樣值的序號(hào)向量 subplot(2,2,1);ste

10、m(n1,f1);title(f1(n);xlabel(n) ;ylabel(f1(n) %在子圖1繪序列f1(k)時(shí)域波形圖subplot(2,2,2); stem(n2,f2); title(f2(n); xlabel(n); ylabel(f2(n); %在子圖2繪序列f2(k)時(shí)波形圖subplot(2,2,3); stem(n,f); title(f1(n)與f2(n的卷積和f(n); xlabel(n); ylabel(f(n); %在子圖3繪序列f(k)的波形圖例2：對(duì)于下面兩個(gè)離散序列，我們可以調(diào)用dconv( )函數(shù)來(lái)求離散序列的卷積。解：MATLAB命令如下：f1=1 2

11、-1 -3 -2 4 ;n1=1:6;f2=ones(1,5);n2=-2:2; %f1，f2兩個(gè)信號(hào)的向量表示f,n=dconv(f1,f2,n1,n2) %子函數(shù)調(diào)用運(yùn)行程序如下所示：f = 1 3 2 -1 -3 0 -2 -1 2 4n = -1 0 1 2 3 4 5 6 7 8f = 1 3 2 -1 -3 0 -2 -1 2 4n = -1 0 1 2 3 4 5 6 7 8需要注意的是，調(diào)用conv( )函數(shù)計(jì)算序列卷積時(shí)，該函數(shù)將向量和以外的序列樣值均視為零，因此，當(dāng)序列或?yàn)闊o(wú)限長(zhǎng)序列時(shí)調(diào)用conv( )函數(shù)就可能出現(xiàn)誤差。如果碰到無(wú)限長(zhǎng)序列時(shí)候，我們必須將其截?cái)嗖拍芮髱氲?/p>

12、conv( )函數(shù)中。此時(shí)，函數(shù)將把截?cái)鄥^(qū)域外的區(qū)間視為零，故conv( )計(jì)算出的卷積只有部分是真實(shí)的。三、實(shí)驗(yàn)內(nèi)容1、對(duì)下面三個(gè)序列，用conv函數(shù)驗(yàn)證卷積特性（交換律、結(jié)合律、分配律） 交換律 結(jié)合律 分配律 其中： 2、用dconv( )函數(shù)求例1的離散卷積和。實(shí)驗(yàn)三 離散系統(tǒng)的z域分析一、z域模型對(duì)于一個(gè)線性的離散系統(tǒng)，MALAB提供了6種表示方法：傳遞函數(shù)法，零極點(diǎn)增益法，狀態(tài)空間法，部分分式展開(kāi)法，二階部分積和Lattice結(jié)構(gòu)。每種方法都適用于不同的場(chǎng)合，考慮到理論知識(shí)的內(nèi)容，在此盡介紹傳遞函數(shù)法，零極點(diǎn)增益法和帶余數(shù)的部分分式展開(kāi)法。1、傳遞函數(shù)法已知，則系統(tǒng)傳遞函數(shù)為：在M

13、ATLAB中用兩個(gè)向量表示該系統(tǒng)。調(diào)用格式為：num=b1,b2,bm,bm+1den=1,a1,an-1,an其中： num，den分別為分子多項(xiàng)式與分母多項(xiàng)式的系數(shù)，都是按的降冪進(jìn)行排列的。2、零極點(diǎn)增益法MATLAB信號(hào)處理工具箱提供的tf2zp和zp2tf函數(shù)可進(jìn)行系統(tǒng)函數(shù)的不同表示形式的轉(zhuǎn)換。正冪有理多項(xiàng)式表示的系統(tǒng)函數(shù)為：用零點(diǎn)、極點(diǎn)和常數(shù)表示系統(tǒng)函數(shù)的零極增益模型為：（1）由傳函模型得到零極增益模型MATLAB工具箱提供的tf2zp函數(shù)可以實(shí)現(xiàn)由系統(tǒng)的傳函模型到零極增益模型的轉(zhuǎn)換。調(diào)用格式為：z,p,k=tf2zp(b,a)其中：調(diào)用參數(shù)為，分別表示的分子和分母多項(xiàng)式的系數(shù)（從高

14、到底排列）。返回參數(shù)為，分別表示零點(diǎn)、極點(diǎn)和增益通過(guò)tf2zp求出系統(tǒng)的零極點(diǎn)后可根據(jù)極點(diǎn)的位置判斷系統(tǒng)的穩(wěn)定性。（2）由零極增益模型得到傳函模型MATLAB工具箱提供的zp2tf 函數(shù)可以實(shí)現(xiàn)由系統(tǒng)的零極增益模型到傳函模型的轉(zhuǎn)換。調(diào)用格式為：b,a=zp2tf(z,p,k)其中：調(diào)用參數(shù)為，分別表示零點(diǎn)、極點(diǎn)和增益返回參數(shù)為，分別表示的分子和分母多項(xiàng)式的系數(shù)（從高到底排列）。例1：試將系統(tǒng)函數(shù)，表示為零極增益形式。解：用正冪有理多項(xiàng)式表示的系統(tǒng)函數(shù)為：MATLAB程序如下：%確定一階因子形式的系統(tǒng)函數(shù)b=1 0 0.04 0;a=1 -0.8 0.16 -0.128;z,p,k=tf2zp(

15、b,a);disp(零點(diǎn));disp(z);disp(極點(diǎn))；disp(p);disp(常數(shù))；disp(k);程序運(yùn)行結(jié)果為：零點(diǎn) 0 0+0.2000i 0-0.2000i極點(diǎn) 0.8000 0.0000+0.4000i 0.0000-0.4000i常數(shù) 1由此可以寫(xiě)出系統(tǒng)函數(shù)的零極增益形式：（3）零極點(diǎn)圖MATLAB工具箱提供了zplane函數(shù)用于繪制系統(tǒng)的零極點(diǎn)圖。調(diào)用格式為：zplane(b,a)其中：調(diào)用參數(shù)為，分別表示的分子和分母多項(xiàng)式的系數(shù)（從高到底排列）。由zplane函數(shù)繪制得到的零極點(diǎn)圖中，用符號(hào)表示零點(diǎn)；表示極點(diǎn)；符號(hào)旁的數(shù)字表示零極點(diǎn)的階數(shù)；虛線畫(huà)的圓表示單位圓。我們

16、可以通過(guò)平面的零極點(diǎn)分布來(lái)判斷系統(tǒng)的穩(wěn)定性。例2：畫(huà)出例1系統(tǒng)的零極點(diǎn)圖，判別其穩(wěn)定性解：MATLAB程序如下：b=1 0 0.04 0;a=1 -0.8 0.16 -0.128;zplane(b,a)程序運(yùn)行結(jié)果為：上由圖可知，該系統(tǒng)極點(diǎn)全部位于單位圓內(nèi)，故系統(tǒng)穩(wěn)定。3、部分分式展開(kāi)法（1）由一般式得到部分分式展開(kāi)式若給出的一般式：可以由函數(shù)residuez得到的部分分式展開(kāi)式。調(diào)用形式為：r,p,k=resifuez(b,a)其中：調(diào)用參數(shù)為，分別表示的分子和分母多項(xiàng)式的系數(shù)（從高到底排列）。返回參數(shù)為，其中為極點(diǎn)向量，為極點(diǎn)對(duì)應(yīng)的系數(shù)向量，為余數(shù)多項(xiàng)式向量。（2）由部分分式展開(kāi)式得到一般

17、式若給出的部分分式展開(kāi)式：可以由函數(shù)residuez得到的一般式。調(diào)用格式為：b,a=residuez(r,p,k)其中：調(diào)用參數(shù)為，其中為極點(diǎn)向量，為極點(diǎn)對(duì)應(yīng)的系數(shù)向量，為余數(shù)多項(xiàng)式向量。返回參數(shù)為，分別表示的分子和分母多項(xiàng)式的系數(shù)（從高到底排列）。所以上述部分分式展開(kāi)式的residuez調(diào)用參數(shù)分別為：r=r1 r2 r3 r4p=p1 p2 p3 p3k=k1 k2同一極點(diǎn)在向量中出現(xiàn)了兩次，表示是一個(gè)二階的重極點(diǎn)。例3：試用MATLAB計(jì)算的部分分式展開(kāi)式。解： MATLAB程序如下： b=1.5,0.98,-2.608,1.2,-0.144; a=1,-1.4,0.6,-0.072;

18、 r,p,k=residuez(b,a) %部分分式展開(kāi) disp(留書(shū))；disp(r) disp(極點(diǎn))；disp(p) disp(常數(shù))；disp(k);程序運(yùn)行結(jié)果為：留數(shù)0.7000+0.0000i 0.5000-0.0000i 0.3000極點(diǎn)0.6000+0.0000i 0.6000-0.0000i 0.2000常數(shù)0 2由此可以寫(xiě)出部分分式展開(kāi)的結(jié)果： 二、z變換和逆z變換變換是時(shí)域離散信號(hào)和系統(tǒng)分析及設(shè)計(jì)的重要的數(shù)學(xué)工具。序列的變換為：對(duì)于給定的，使收斂的所有值的集合稱為變換的收斂域，即滿足：的反變換定義為：MATLAB的符號(hào)數(shù)學(xué)工具箱提供了計(jì)算正變換的函數(shù)ztrans和計(jì)算

19、逆變換的函數(shù)iztrans。調(diào)用形式為：F=ztrans(f)該調(diào)用方式用于求符號(hào)函數(shù)的變換，返回函數(shù)的自變量為；F=ztrans(f,w)該調(diào)用方式用于求符號(hào)函數(shù)的變換，返回函數(shù)的自變量為；F=ztrans(f,n,w)該調(diào)用方式用于對(duì)自變量為的符號(hào)函數(shù)求變換，返回函數(shù)的自變量為。f=iztrans(F)該調(diào)用方式用于對(duì)自變量為的符號(hào)函數(shù)求逆變換，返回函數(shù)的自變量為；f=iztrans(F,n)該調(diào)用方式用于對(duì)自變量為的符號(hào)函數(shù)求逆變換，返回函數(shù)的自變量為；f=iztrans(F,w,n)該調(diào)用方式用于對(duì)自變量為的符號(hào)函數(shù)求逆變換，返回函數(shù)的自變量為。例4：已知序列，求其變換。解：MATLA

20、B程序如下：syms nf=sym(2(-n); %定義序列F=ztrans(f) %求z變換程序運(yùn)行結(jié)果為：F =2*z/(2*z-1)即：例5：已知一離散系統(tǒng)的系統(tǒng)函數(shù)，求其沖激響應(yīng)。解：MATLAB程序如下：syms n zH=sym(z/(z2+3*z+2);h=iztrans(H,n) %求逆z變換程序運(yùn)行結(jié)果為：h =(-1)k-(-2)k，即：對(duì)象函數(shù)求逆變換，還可以利用residuez( )對(duì)象函數(shù)作部分分式展開(kāi)，然后按部分分式展開(kāi)法求得原函數(shù)。在此不做介紹。三、離散系統(tǒng)的頻率響應(yīng)分析若離散系統(tǒng)是穩(wěn)定的，其系統(tǒng)函數(shù)的收斂域應(yīng)包含單位圓，離散系統(tǒng)的頻率響應(yīng)即為單位圓上（）的系統(tǒng)函

21、數(shù)，即：其中：為系統(tǒng)的幅頻特性，為系統(tǒng)的相頻特性。在MATLAB中，利用freqz( )函數(shù)可方便的求出系統(tǒng)的頻率響應(yīng)。調(diào)用格式為：freqz(b,a)該調(diào)用方式將繪制系統(tǒng)在范圍內(nèi)的幅頻特性和相頻特性圖；其中：分別為系統(tǒng)函數(shù)的分子、分母多項(xiàng)式的系數(shù)向量。freqz(b,a,whole)該調(diào)用方式將繪制系統(tǒng)在范圍內(nèi)的幅頻特性和相頻特性圖。freqz(b,a,N)該調(diào)用方式將繪制系統(tǒng)在范圍內(nèi)個(gè)頻率等分點(diǎn)的幅頻特性和相頻特性圖，的缺省值為512。freqz(b,a,N,whole)該調(diào)用方式將繪制系統(tǒng)在范圍內(nèi)個(gè)頻率等分點(diǎn)的幅頻特性和相頻特性圖。此外，還有如下相類(lèi)似的四種調(diào)用形式。H,w= freqz

22、(b,a)H,w=freqz(b,a,whole) H,w=freqz(b,a,N) H,w=freqz(b,a,N,whole) 其中，返回向量包含了離散系統(tǒng)頻率響應(yīng)在（或）范圍內(nèi)各頻率點(diǎn)處的值，返回向量則包含了在（或）范圍內(nèi)個(gè)（的缺省值為512）頻率等分點(diǎn)。利用這些調(diào)用方式MATLAB并不直接繪制系統(tǒng)的頻率特性圖，但可由向量用abs、angle、plot等函數(shù)來(lái)繪制幅頻特性和相頻特性圖。例6：已知一離散系統(tǒng)的系統(tǒng)函數(shù)，試?yán)L制頻率特性圖。解：MATLAB程序如下：b=0 1 0;a=1 0.3 0.2;freqz(b,a,whole)程序運(yùn)行結(jié)果為：四、實(shí)驗(yàn)內(nèi)容1、已知一因果的LTI系統(tǒng)的函

23、數(shù)為，求系統(tǒng)的零極增益模型，部分分式模型，并用函數(shù)zplane畫(huà)出系統(tǒng)的零極點(diǎn)分布，判斷系統(tǒng)的穩(wěn)定性。2、已知因果系統(tǒng)的系統(tǒng)函數(shù)為：求系統(tǒng)的單位脈沖響應(yīng)并顯示其波形圖；若輸入為，求系統(tǒng)的零狀態(tài)響應(yīng)并顯示其波形圖。實(shí)驗(yàn)四 離散傅立葉變換（DFT）DFT是數(shù)字信號(hào)處理中最重要的數(shù)學(xué)工具之一。其實(shí)質(zhì)是對(duì)有限長(zhǎng)序列頻譜的離散化，即通過(guò)DFT使時(shí)域有限長(zhǎng)序列與頻域有限長(zhǎng)序列相對(duì)應(yīng)，從而可在頻率域用計(jì)算機(jī)進(jìn)行信號(hào)處理。更重要的是DFT有多種快速算法（FFT），可使信號(hào)處理速度提高好幾倍，使數(shù)字信號(hào)的實(shí)時(shí)處理得以實(shí)現(xiàn)。因此，DFT即具有重要的理論意義又有廣泛的實(shí)際應(yīng)用價(jià)值。熟悉DFT的定義、物理意義和重要性

24、質(zhì)，有助于正確使用DFT解決數(shù)字信號(hào)處理的實(shí)際問(wèn)題。本次試驗(yàn)主要學(xué)習(xí)用MATLAB工具箱函數(shù)，以序列的時(shí)域和頻域波形直觀的驗(yàn)證DFT的物理意義及頻域采樣理論。一、DFT1、DFT定義設(shè)序列的長(zhǎng)度為，則的點(diǎn)離散傅立葉變換對(duì)定義為：正變換： （1）逆變換： （2）其中，成為變換區(qū)間的長(zhǎng)度。把（1）式改寫(xiě)為矩陣乘法運(yùn)算：其中為序列行向量，是階矩陣，通常稱為旋轉(zhuǎn)因子矩陣。上式用MATLAB的矩陣運(yùn)算表示為：2、MATLAB程序點(diǎn)DFT的MATLAB程序如下：用矩陣乘法計(jì)算N點(diǎn)DFTclear; close allxn=input(請(qǐng)輸入序列x=)N=length(xn);n=0:N-1; k=n; n

25、k=n*k; %生成0:N-1*0:N-1方陣WN=exp(-j*2*pi/N);Wnk=WN.nk; %產(chǎn)生旋轉(zhuǎn)因子X(jué)k=xn*Wnk; %計(jì)算N點(diǎn)DFT只要輸入序列，運(yùn)行該程序，即可實(shí)現(xiàn)的點(diǎn)DFT優(yōu)點(diǎn)：概念清楚，編程簡(jiǎn)單缺點(diǎn)：占用內(nèi)存大，運(yùn)行速度低，不太實(shí)用二、FFT快速傅里葉變換簡(jiǎn)稱FFT，它并不是與DFT不同的另一種變換，而是為了減少DFT計(jì)算次數(shù)的一種快速有效的算法。FFT算法基本可以分為兩類(lèi)：按時(shí)間抽取的FFT（DIT-FFT）算法和按頻率抽取的FFT（DIF-FFT）。若，其中為整數(shù)，這時(shí)候的FFT稱為基2FFT。關(guān)于FFT算法的具體原理在這里不作介紹，請(qǐng)大家自己查閱相關(guān)資料。M

26、ATLAB基礎(chǔ)部分提供了fft， ifft， fft2和ifft2等快速計(jì)算傅立葉變換的函數(shù)，它使DFT的運(yùn)算速度量提高了若干的數(shù)量級(jí)。信號(hào)處理工具箱提供了等有關(guān)的變換函數(shù)。為了簡(jiǎn)化程序，后面的例題均直接調(diào)用這些函數(shù)。和用于計(jì)算一維快速正逆傅立葉變換，其調(diào)用格式如下：X=fft(x,N)用于計(jì)算序列向量的點(diǎn)DFT。默認(rèn)時(shí)fft函數(shù)自動(dòng)按的長(zhǎng)度計(jì)算DFT。當(dāng)為2的整數(shù)次冪時(shí)，fft按基2算法計(jì)算，否則用混合基算法。 調(diào)用格式相仿。例1：基本序列的離散傅立葉變換計(jì)算已知一下序列：復(fù)正弦序列：余弦序列：正弦序列：分別對(duì)和計(jì)算以上序列的點(diǎn)DFT，并繪出幅頻特性曲線。解：直接產(chǎn)生序列調(diào)用函數(shù)求解本題的程

27、序如下：clear;close allN=16;N1=8;n=0:N-1;k=0:N1-1; %產(chǎn)生序列x1(n),計(jì)算DFTx1(n)x1n=exp(j*pi*n/8); %產(chǎn)生序列x1(n)X1k=fft(x1n,N); %計(jì)算N點(diǎn)DFTx1(n)Xk1=fft(x1n,N1); %計(jì)算N1點(diǎn)DFTx1(n)%產(chǎn)生序列x2(n),計(jì)算DFTx2(n)x2n=cos(pi*n/8); X2k=fft(x2n,N); %計(jì)算N點(diǎn)DFTx2(n)Xk2=fft(x2n,N1); %計(jì)算N1點(diǎn)DFTx2(n)%產(chǎn)生序列x3(n),計(jì)算DFTx3(n)x3n=sin(pi*n/8); X3k=fft

28、(x3n,N); %計(jì)算N點(diǎn)DFTx3(n)Xk3=fft(x3n,N1); %計(jì)算N1點(diǎn)DFTx3(n)用，語(yǔ)句畫(huà)出這三種信號(hào)、兩種采樣序列長(zhǎng)度、共六種情況的曲線。圖形劃分及標(biāo)注語(yǔ)句從略。運(yùn)行結(jié)果如下：三、驗(yàn)證N點(diǎn)DFT物理意義長(zhǎng)度為的序列的傅立葉變換和點(diǎn)離散傅立葉變換的定義式如下：式中，比較上面兩式可以得出的傅立葉變換與的離散傅立葉變換的關(guān)系式為：簡(jiǎn)言之，序列點(diǎn)DFT的物理意義就是對(duì)在上進(jìn)行點(diǎn)等間隔采樣。為了驗(yàn)證這種關(guān)系，程序先對(duì)密集采樣，繪制出幅頻曲線圖。然后再分別做8點(diǎn)和16點(diǎn)DFT來(lái)驗(yàn)證如上采樣關(guān)系。為了能直觀的看出和之間的采樣關(guān)系，再的幅度和相位圖中同時(shí)畫(huà)出了的幅頻和相頻曲線。例2

29、：（1），繪出幅頻曲線和相頻曲線（2）計(jì)算并圖示的8點(diǎn)DFT（3）計(jì)算并圖示的16點(diǎn)DFT解：MATLAB程序如下：%計(jì)算并圖示DFTx(n)clear;close allN1=8;N2=16; %兩種FFT變換長(zhǎng)度n=0:N1-1;w=2*pi*(0:2047)/2048Xw=(1-exp(-j*4*w)./(1-exp(-j*w) %對(duì)x(n)頻譜函數(shù)采樣2048點(diǎn)subplot(3,2,1);plot(w/pi,abs(Xw)title(x(n)的幅頻曲線);xlabel(w/pi);ylabel(幅度);subplot(3,2,2);plot(w/pi,angle(Xw)title(x

30、(n)的相頻曲線);axis(0,2,-pi,pi);line(0,2,0,0)xlabel(w/pi);ylabel(相位(rad);xn=(n=0)&(n4) %產(chǎn)生序列x(n)=xnX1k=fft(xn,N1); %計(jì)算N1點(diǎn)的DFTx(n)X2k=fft(xn,N2); %計(jì)算N2點(diǎn)的DFTx(n)figure(2)k1=0:N1-1subplot(3,2,1);stem(k1,abs(X1k),.) %畫(huà)N1點(diǎn)離散頻譜幅度title(N1點(diǎn)DFTx(n)X1k)xlabel(k(w=2pik/N1);ylabel(X1k幅度);hold onplot(N1/2*w/pi,abs(Ww

31、) %疊加上連續(xù)頻譜幅度subplot(3,2,2);stem(k1,angle(X1k),.) %畫(huà)N1點(diǎn)離散頻譜相位title(X1k的相頻曲線);axis(0,N1,-pi,pi);line(0,N1,0,0)xlabel(k(w=2pik/N1);ylabel(相位(rad);hold on plot(N1/2*w/pi,angle(Xw) %疊加上連續(xù)頻譜相位subplot(3,2,3);stem(k1,angle(X2k),.) %畫(huà)N2點(diǎn)離散頻譜相位運(yùn)行結(jié)果為：四、實(shí)驗(yàn)內(nèi)容1、設(shè)是長(zhǎng)度為的矩形序列，用matlab程序分析FFT取不同長(zhǎng)度時(shí)的頻譜變化2、設(shè)完成如下要求：（1）計(jì)算其

32、DTFT，并畫(huà)出區(qū)間的波形（2）計(jì)算4點(diǎn)DFT，并把結(jié)果顯示在（1）所畫(huà)的圖形中（3）對(duì)補(bǔ)零，計(jì)算64點(diǎn)DFT，并顯示結(jié)果。3、思考：用DFT理論解釋例1中為何兩種N值下的DFT結(jié)果差別如此之大？實(shí)驗(yàn)五 IIR數(shù)字濾波器的設(shè)計(jì)IIR數(shù)字濾波器設(shè)計(jì)的主要方法是先設(shè)計(jì)低通模擬濾波器，進(jìn)行頻率變換，將其轉(zhuǎn)換成相應(yīng)的（高通的、帶通的）模擬濾波器，再轉(zhuǎn)換成為高通、帶通或帶阻數(shù)字濾波器。對(duì)設(shè)計(jì)的全過(guò)程的各個(gè)步驟，MATLAB都提供了相應(yīng)的工具箱函數(shù)，使IIR數(shù)字濾波器設(shè)計(jì)變得非常簡(jiǎn)單。一、工具箱函數(shù)介紹 IIR數(shù)字濾波器的設(shè)計(jì)步驟可由下圖的流程圖來(lái)表示。濾波器系數(shù)模擬低通濾波器原型設(shè)計(jì)buttap，che

33、blap，cheb2ap，besselap，ellipap頻率變換（變?yōu)楦咄?、帶通、帶阻等）lp2lp，lp2hp，lp2bp，lp2bs求最小階數(shù)buttord，cheblord，cheb2ord，ellipord合為一步的設(shè)計(jì)函數(shù)butter，cheb1，cheb2，ellip，besself模擬數(shù)字變換bilinearimpinvar設(shè)計(jì)指標(biāo)濾波器系數(shù)圖1 IIR數(shù)字濾波器設(shè)計(jì)步驟流程圖上圖清楚的表示了5類(lèi)20個(gè)信號(hào)處理工具箱函數(shù)的作用。下面以巴特沃斯濾波器設(shè)計(jì)函數(shù)為典型，介紹此流圖中函數(shù)的功能和用法。其它類(lèi)型的濾波器設(shè)計(jì)函數(shù)用法可以類(lèi)推。1、求最小階數(shù)N的函數(shù)buttord調(diào)用格式為：

34、N,wc= buttord(wp,ws,Ap,As,s)其中：調(diào)用參數(shù)為濾波器指標(biāo)；返回參數(shù)為巴特沃斯模擬濾波器的階數(shù)及頻率參數(shù)【注】及均在區(qū)間歸一化，單位為弧度/秒；去掉最后的變?cè)?，得到的就是?shù)字濾波器，調(diào)用格式為： N,wc= buttord(wp,ws,Ap,As)其中：調(diào)用參數(shù)為濾波器指標(biāo)；返回參數(shù)為巴特沃斯數(shù)字濾波器的階數(shù)及頻率參數(shù)2、模擬低通濾波器原型設(shè)計(jì)函數(shù)buttap調(diào)用格式為：z,p,k= buttap(N)其中：調(diào)用參數(shù)為butterworth濾波器階數(shù)；返回參數(shù)分別為濾波器的零點(diǎn)、極點(diǎn)和增益。得到后，很容易求出濾波器系數(shù)。3、模擬頻率變換函數(shù)lp2lp調(diào)用格式為： Bt

35、,At= lp2lp(B,A,wo)其中：調(diào)用參數(shù)為單位截止頻率的模擬低通濾波器分母和分子多項(xiàng)式的系數(shù)；（單位：弧度/秒）為低通濾波器所期望的截止頻率； 返回參數(shù)為得到的頻率為的低通濾波器的分母和分子多項(xiàng)式的系數(shù)。4、模擬數(shù)字變換函數(shù)模擬數(shù)字變換函數(shù)包括雙線性變換函數(shù)bilinear和脈沖響應(yīng)不變法函數(shù)impinvar兩種。（1）雙線性變換函數(shù)bilinear MATLAB工具箱提供了函數(shù)bilinear，使用雙線性變換法實(shí)現(xiàn)模擬濾波器的數(shù)字化。調(diào)用格式為：Bz,Az= bilinear (B,A,Fs)其中：調(diào)用參數(shù)為模擬濾波器傳遞函數(shù)分母和分子多項(xiàng)式的系數(shù)；為抽樣頻率； 返回參數(shù)為用雙線性

36、變換法設(shè)計(jì)得到的數(shù)字濾波器分母和分子多項(xiàng)式的系數(shù)。（2）脈沖響應(yīng)不變法函數(shù)impinvarMATLAB工具箱提供了函數(shù)impinvar，使用脈沖響應(yīng)不變法實(shí)現(xiàn)模擬濾波器的數(shù)字化。調(diào)用格式為： Bz,Az=impinvar(B,A,Fs)其中：調(diào)用參數(shù)為模擬濾波器傳遞函數(shù)分母和分子多項(xiàng)式的系數(shù)；為抽樣頻率；其默認(rèn)值為。 返回參數(shù)為用脈沖響應(yīng)不變法設(shè)計(jì)得到的數(shù)字濾波器分母和分子多項(xiàng)式系數(shù)。5、把（2），（3），（4）合為一步的數(shù)字濾波器設(shè)計(jì)函數(shù)butter（N,wc,ftype）B,A=butter(N,Wc, ftype)其中：調(diào)用參數(shù)為階數(shù)和截止頻率； 返回參數(shù)設(shè)計(jì)得到的數(shù)字濾波器分母和分子多

37、項(xiàng)式的系數(shù)；ftype用來(lái)控制設(shè)計(jì)得到的濾波器是低通、高通、帶通還是帶阻濾波器。B,A=butter（N,wc）設(shè)計(jì)低通或帶通數(shù)字濾波器系數(shù)。 B,A=butter（N,wc,high）設(shè)計(jì)高通數(shù)字濾波器系數(shù)。 B,A=butter（N,wc,stop）設(shè)計(jì)帶阻數(shù)字濾波器系數(shù)?！咀ⅰ浚?）采用的是雙線性變換函數(shù)bilinear，如果用脈沖響應(yīng)不變法就得分開(kāi)來(lái)作。由此可見(jiàn)：通常情況之下，有了（1）和（5）兩類(lèi)函數(shù)，數(shù)字濾波器設(shè)計(jì)問(wèn)題就解決了。二、IIR濾波器設(shè)計(jì)1、低通巴特沃斯模擬濾波器設(shè)計(jì)低通巴特沃斯模擬濾波器的系統(tǒng)函數(shù)一般形式如下：極點(diǎn)：由此可見(jiàn)，低通巴特沃斯濾波器的系統(tǒng)函數(shù)完全由階數(shù)和截止

38、頻率決定。而和是由濾波器設(shè)計(jì)指標(biāo)決定的，其計(jì)算公式如下：，?。核O(shè)計(jì)的濾波器的通帶指標(biāo)剛好滿足，阻帶指標(biāo)富裕。?。核O(shè)計(jì)的濾波器的阻帶指標(biāo)剛好滿足，通帶指標(biāo)富裕。MATLAB工具箱buttord，butter就是根據(jù)以上關(guān)系式編寫(xiě)的，因此設(shè)計(jì)的時(shí)候不需要再記憶和使用這些公式了。例1：設(shè)計(jì)一個(gè)低通巴特沃斯模擬濾波器指標(biāo)如下：通帶截止頻率：，通帶最大衰減：通帶截止頻率：，阻帶最小衰減：解：MATLAB程序如下：%低通巴特沃斯低通濾波器設(shè)計(jì)clear;close allfp=3400;fs=4000;Ap=3;As=40; %輸入濾波器指標(biāo)N,fc=buttord(fp,fs,Ap,As,s) %計(jì)

39、算階數(shù)N和3dB截止頻率fcB,A=butter(N,fc,s); %設(shè)計(jì)低通巴特沃斯模擬濾波器hf,f=freqs(B,A,1024); %計(jì)算模擬濾波器頻率響應(yīng)，freqs為工具函數(shù)plot(f,20*log10(abs(hf)/abs(hf(1)grid;xlabel(f/Hz);ylabel(幅度(dB);axis(0,4000,-40,5)line(0,4000,-3,-3);line(0,4000,-3,-3)line(3400,3400,-90,5)2、模擬低通濾波器轉(zhuǎn)換成數(shù)字低通濾波器模擬濾波器離散化的基本方法有脈沖響應(yīng)不變法和雙線性變換法。1、脈沖響應(yīng)不變法函數(shù)impinva

40、r 對(duì)只有單階極點(diǎn)的階模擬濾波器，用部分分式展開(kāi)為：則脈沖響應(yīng)不變法?。篗ATLAB工具箱函數(shù)impinvar用來(lái)實(shí)現(xiàn)這個(gè)轉(zhuǎn)換。2、雙線性變換函數(shù)bilinear雙線性變換原理是用替換中的值，得到。MATLAB工具箱函數(shù)bilinear用來(lái)實(shí)現(xiàn)這個(gè)轉(zhuǎn)換。例2：已知一模擬濾波器的系統(tǒng)函數(shù)為，分別用脈沖響應(yīng)不變法和雙線性變換法將轉(zhuǎn)換成數(shù)字濾波器系統(tǒng)函數(shù)，并圖示和的幅頻響應(yīng)曲線。分別取采樣頻率和觀察脈沖響應(yīng)不變法中存在的頻率混疊失真和雙線性變換法中存在的非線性頻率失真。解：MATLAB程序如下：%用脈沖響應(yīng)不變法和雙線性變換法將模擬濾波器離散化clear;close allb=1000;a=1,10

41、00;w=0:1000*2*pi; %設(shè)定模擬頻率hf,w=freqs(b,a,w); %計(jì)算模擬濾波器的頻響函數(shù)subplot(2,3,1) %畫(huà)出模擬濾波器的幅頻特性plot(w/2/pi,abs(hf);grid;title(a)xlabel(f(Hz);ylabel(幅度);Fs0=1000,500;for m=1:2 Fs=Fs0(m) %T=0.001s及T=0.002s d,c=impinvar(b,a,Fs) %用impinvar函數(shù)實(shí)現(xiàn)離散化 wd=0:512*pi/512 %設(shè)定數(shù)字歸一化頻率 hw1=freqz(d,c,wd) % 計(jì)算數(shù)字濾波器頻響函數(shù)subplot(2

42、,3,2);plot(wd/pi,abs(hw1)/abs(hw1(1);hold on; %畫(huà)出數(shù)字濾波器幅頻特性endgrid;xlabel(f/Hz);text(0.25,0.88,T=0.002s);text(0.12,0.54,T=0.001s);for m=1:2 Fs=Fs0(m) %T=0.001s及T=0.002s f,e=bilinear(b,a,Fs) %用bilinear函數(shù)實(shí)現(xiàn)離散化 wd=0:512*pi/512 %設(shè)定數(shù)字歸一化頻率 hw2=freqz(f,e,wd) % 計(jì)算數(shù)字濾波器頻響函數(shù)subplot(2,3,3);plot(wd/pi,abs(hw2)/

43、abs(hw2(1);hold on; %畫(huà)出數(shù)字濾波器幅頻特性endgrid;xlabel(f/Hz);text(0.5,0.74,T=0.002s);text(0.12,0.34,T=0.001s);三、實(shí)驗(yàn)內(nèi)容1、運(yùn)行例題一及例題二程序；顯示運(yùn)行結(jié)果；根據(jù)運(yùn)行結(jié)果說(shuō)明脈沖響應(yīng)不變法和雙線性變換法各自的特點(diǎn)。2、用雙線性法設(shè)計(jì)一個(gè)Butterworth低通數(shù)字濾波器，頻率在處衰減，在處衰減，采樣頻率為，確定系統(tǒng)函數(shù)。并圖示和的幅頻響應(yīng)曲線?！咎崾尽繑?shù)字低通雙線性變換與脈沖不變法程序：wp=500*2*pi;ws=150*2*pi;Ap=3;As=15;Fs=2000;n,Wn=buttor

44、d(wp,ws,Ap,As,s);b,a=butter(n,Wn,s);bn1,an1=bilinear(b,a,2000); %雙線性變換H1,w=freqz(bn1,an1);%采用脈沖響應(yīng)不變法進(jìn)行離散化處理bn2,an2=impinvar(b,a,2000); %脈沖響應(yīng)不變法H2,w=freqz(bn2,an2)實(shí)驗(yàn)六 FIR數(shù)字濾波器的設(shè)計(jì) IIR數(shù)字濾波器的設(shè)計(jì)簡(jiǎn)單方便，特別是采用雙線性變換法設(shè)計(jì)的數(shù)字濾波器不存在頻譜混疊問(wèn)題，效果較好。但是IIR數(shù)字濾波器有一個(gè)明顯的缺點(diǎn)，就是其相位特性一般是非線性的。如果濾波器在有效傳輸頻帶內(nèi)的相位特性不是線性的，將造成有用信號(hào)的傳輸失真。在

45、有些實(shí)際應(yīng)用場(chǎng)合，例如數(shù)據(jù)傳輸，圖像處理等，對(duì)濾波器的線性相位特性要求頗為嚴(yán)格，所以在這些場(chǎng)合中IIR濾波器一般是不能勝任的。FIR濾波器最大的優(yōu)點(diǎn)是容易設(shè)計(jì)成線性相位特性，其幅度特性可以隨意設(shè)計(jì)，而且不存在穩(wěn)定性問(wèn)題。設(shè)FIR濾波器的單位脈沖響應(yīng)為，則其系統(tǒng)函數(shù)為：上式表明，F(xiàn)IR濾波器的系統(tǒng)函數(shù)為的多項(xiàng)式，而IIR濾波器的系統(tǒng)函數(shù)為的有理分式形式。因此，F(xiàn)IR濾波器在S平面上找不到與之響應(yīng)的模擬系統(tǒng)函數(shù)，也就是說(shuō)，F(xiàn)IR濾波器的設(shè)計(jì)不能借用模擬濾波器設(shè)計(jì)的一套成熟理論。MATLAB信號(hào)處理工具箱提供的FIR數(shù)字濾波器的設(shè)計(jì)方法有以下兩種：窗函數(shù)法，等波紋最佳一致逼近法。本次實(shí)驗(yàn)主要介紹窗函

46、數(shù)法。一、MATLAB提供的窗函數(shù)窗函數(shù)法設(shè)計(jì)濾波器時(shí)，階數(shù)選擇如下表所示：表1 六種窗函數(shù)特性的表窗函數(shù)旁瓣峰值近似過(guò)渡帶精確過(guò)渡帶最小阻帶衰減矩形窗巴特利特漢寧海明布拉克曼凱澤窗Matlab提供了相應(yīng)的子程序來(lái)實(shí)現(xiàn)本節(jié)中的窗函數(shù)，如下：wd=boxcar(N) %數(shù)組中返回點(diǎn)矩形窗函數(shù)wd=triang(N) %數(shù)組中返回點(diǎn)Bartlett（三角）窗函數(shù)wd=hanning(N) %數(shù)組中返回點(diǎn)漢寧窗函數(shù)wd=hamming(N) %數(shù)組中返回點(diǎn)海明窗函數(shù)wd=blackman(N) %數(shù)組中返回點(diǎn)布拉克曼窗函數(shù)wd=kaiser(N,beta) %數(shù)組中返回給定beta值時(shí)點(diǎn)凱澤窗函數(shù)【

47、注】以上函數(shù)的輸入變?cè)话阒灰昂瘮?shù)的長(zhǎng)度就夠了，只有凱澤窗還需要規(guī)定值。輸出變?cè)褪侵行闹禋w一化為1的窗函數(shù)序列，它是列向量。利用語(yǔ)句plot就可以畫(huà)出窗函數(shù)序列的形狀。二、窗函數(shù)法設(shè)計(jì)FIR DF的步驟1、設(shè)計(jì)步驟（1）根據(jù)對(duì)過(guò)渡帶寬及阻帶衰減要求，選擇窗函數(shù)的類(lèi)型并估計(jì)窗口長(zhǎng)度（或階數(shù)）窗函數(shù)類(lèi)型可根據(jù)其阻帶最小衰減的條件獨(dú)立選擇，因?yàn)槠溟L(zhǎng)度對(duì)阻帶最小衰減沒(méi)有影響。在確定窗函數(shù)之后，就可根據(jù)過(guò)渡帶寬小于給定指標(biāo)的條件，確定所擬用的窗函數(shù)長(zhǎng)度。（2）根據(jù)和求出相應(yīng)的理想濾波器單位脈沖響應(yīng)。在一般情況下，是不能用封閉公式表示的，需要采用數(shù)值方法。從到采樣點(diǎn)，用IDFT來(lái)求出。（3）計(jì)算FIR

48、數(shù)字濾波器單位脈沖響應(yīng)。由MATLAB得出的是一個(gè)無(wú)限長(zhǎng)序列，因而不可用FIR濾波器實(shí)現(xiàn)。所以要選擇合適的窗函數(shù)來(lái)截取的適當(dāng)長(zhǎng)度（即階數(shù)），以保證實(shí)現(xiàn)要求的阻帶衰減；在matlab中用點(diǎn)乘命令表示為：（4）驗(yàn)證技術(shù)指標(biāo)是否滿足要求。為了計(jì)算數(shù)字濾波器在頻域中的特性，可以用freqz子程序。例1：分別用矩形窗和Hamming窗設(shè)計(jì)線性相位的FIR數(shù)字低通濾波器。要求通帶截止頻率，單位脈沖響應(yīng)的長(zhǎng)度。繪出及其幅頻響應(yīng)特性曲線。解：【建?！?用窗函數(shù)法設(shè)計(jì)FIR濾波器時(shí)，先求出相應(yīng)的理想低通濾波器單位脈沖響應(yīng)，在根據(jù)阻帶最小衰減選擇合適的窗函數(shù)，最后得到FIR數(shù)字濾波器單位脈沖響應(yīng)。本題中，所以線性相位理想低通濾波器的單位脈沖響應(yīng)為：為了滿足線性相位FIR濾波器條件，要求【程序】%窗函數(shù)法設(shè)計(jì)FIR低