數(shù)字信號(hào)處理實(shí)驗(yàn)講義2012

1、數(shù)字信號(hào)處理講義 數(shù)字信號(hào)處理實(shí)驗(yàn)指導(dǎo)書   曲阜師范大學(xué)物理工程學(xué)院信號(hào)與信息處理教研室2012.0369 前 言數(shù)字信號(hào)處理是是電子信息工程專業(yè)的一門專業(yè)基礎(chǔ)課。本課程主要研究如何對信號(hào)進(jìn)行分析、變換、綜合、估計(jì)與識(shí)別等加工處理的基本理論和方法，數(shù)字信號(hào)處理實(shí)驗(yàn)是驗(yàn)證、鞏固和補(bǔ)充課堂講授的理論知識(shí)的必要環(huán)節(jié)。通過實(shí)驗(yàn)，使學(xué)生鞏固所學(xué)基本理論，掌握最基本的數(shù)字信號(hào)處理的理論和方法，提高綜合運(yùn)用所學(xué)知識(shí)，提高計(jì)算機(jī)編程的能力。進(jìn)一步加強(qiáng)學(xué)生獨(dú)立分析問題、解決問題的能力、綜合設(shè)計(jì)及創(chuàng)新能力的培養(yǎng)，同時(shí)注意培養(yǎng)學(xué)生實(shí)事求是、嚴(yán)肅認(rèn)真的科學(xué)作風(fēng)和良好的實(shí)驗(yàn)習(xí)慣

2、，為今后的工作打下良好的基礎(chǔ)。數(shù)字信號(hào)處理實(shí)驗(yàn)指導(dǎo)書針對每個(gè)實(shí)驗(yàn)介紹了MATLAB語言數(shù)字信號(hào)處理工具箱中的相應(yīng)函數(shù)，舉例并附有相應(yīng)的程序。為配合課堂理論教學(xué)，實(shí)驗(yàn)內(nèi)容安排仍從認(rèn)識(shí)性和驗(yàn)證性入手，逐步增加設(shè)計(jì)性和工程應(yīng)用性內(nèi)容，達(dá)到訓(xùn)練實(shí)驗(yàn)技能和積累工程實(shí)際應(yīng)用經(jīng)驗(yàn)之目的。數(shù)字信號(hào)處理實(shí)驗(yàn)成績占本門課程總成績的20%。實(shí) 驗(yàn) 要 求在實(shí)驗(yàn)過程中，要求學(xué)生做到：1、預(yù)習(xí)實(shí)驗(yàn)指導(dǎo)書有關(guān)部分，認(rèn)真做好實(shí)驗(yàn)內(nèi)容的準(zhǔn)備工作，就實(shí)驗(yàn)可能出現(xiàn)的情況提前作出思考和分析，需要計(jì)算的參數(shù)提前完成計(jì)算工作，并認(rèn)真寫出預(yù)習(xí)報(bào)告。2、仔細(xì)觀察實(shí)驗(yàn)過程中圖形隨參數(shù)的變化，記錄圖形變化的主要情況，作出必要說明和分析。3、認(rèn)

3、真書寫實(shí)驗(yàn)報(bào)告并在規(guī)定的時(shí)間內(nèi)把實(shí)驗(yàn)報(bào)告交給輔導(dǎo)教師。實(shí)驗(yàn)報(bào)告包括實(shí)驗(yàn)?zāi)康暮鸵?，?shí)驗(yàn)情況及其分析。對需要編程的實(shí)驗(yàn)，寫出程序設(shè)計(jì)說明，給出源程序框圖和清單以及實(shí)驗(yàn)結(jié)論、相關(guān)圖表，并做必要的分析總結(jié)。4、遵守機(jī)房紀(jì)律，服從輔導(dǎo)教師指揮，愛護(hù)實(shí)驗(yàn)設(shè)備，注意衛(wèi)生。5、實(shí)驗(yàn)課程不遲到。如有事不能出席，所缺實(shí)驗(yàn)一般不補(bǔ)。實(shí)驗(yàn)驗(yàn)收分為兩個(gè)部分。第一部分是上機(jī)操作，包括檢查程序運(yùn)行和即時(shí)提問。第二部分是提交書面的實(shí)驗(yàn)報(bào)告。每個(gè)實(shí)驗(yàn)都應(yīng)當(dāng)在規(guī)定的時(shí)間內(nèi)完成并檢查通過，過期視為未完成該實(shí)驗(yàn)，扣該實(shí)驗(yàn)操作成績。為避免期末集中檢查方式產(chǎn)生的諸多不良問題，希望同學(xué)們抓緊時(shí)間，合理安排，認(rèn)真完成。目 錄實(shí)驗(yàn)一 MAT

4、LAB基本命令1實(shí)驗(yàn)二 離散時(shí)間信號(hào)分析4實(shí)驗(yàn)三 離散系統(tǒng)分析17實(shí)驗(yàn)四 離散傅立葉變換30實(shí)驗(yàn)五 快速傅里葉變換（FFT）及其應(yīng)用37實(shí)驗(yàn)六 IIR濾波器設(shè)計(jì)42實(shí)驗(yàn)七 FIR濾波器設(shè)計(jì)55實(shí)驗(yàn)八 窗函數(shù)及其在譜分析中的作用59實(shí)驗(yàn)一 MATLAB基本命令一、實(shí)驗(yàn)?zāi)康模赫莆誐atlab基本命令的使用，以便在以后的實(shí)驗(yàn)和進(jìn)一步的學(xué)習(xí)中使用Matlab命令。實(shí)驗(yàn)任務(wù)：1、 熟悉Matlab 語言環(huán)境，學(xué)習(xí)基本命令的使用，并要求掌握。2、 熟悉Matlab中的信號(hào)處理Demo使用，并了解其中的程序結(jié)構(gòu)和命令。MATLAB代表Matrix Laboratory。 MATLAB 最早是為方便矩陣存取而

5、編寫的一種軟件，現(xiàn)已發(fā)展成為一個(gè)具有高性能數(shù)值計(jì)算和可視化功能的科學(xué)計(jì)算環(huán)境。MATLAB集成了數(shù)值分析、矩陣計(jì)算、信號(hào)處理和圖形等眾多功能，具有簡單易用的環(huán)境，問題的提出和解答只須以數(shù)學(xué)方式表達(dá)和描述，不需要大量原始而傳統(tǒng)的編程過程。 MATLAB 一些簡單命令：1. Demo命令：MATLAB例程演示初次學(xué)習(xí)MATLAB可以使用demo命令來熟悉MATLAB語言環(huán)境。你可以在MATLAB的窗口中直接敲demo，它將帶你學(xué)習(xí)使用MATLAB。2. Help命令：在線幫助。語法格式為Help topic topic 為要取得的幫助主題。例如：help fft 將顯示fft( )函數(shù)的說明。3.

6、 編寫MATLAB程序：MATLAB 的程序以 .M 結(jié)尾，你可以使用MATLAB窗口中File-àNewàM-file 的菜單命令來建立一個(gè).m空白文件，編寫完成后，你保存一下，然后可以在MATLAB窗口中敲此文件名執(zhí)行該文件。 函數(shù)編寫 script編寫4. 常用的MATLAB命令歸類：通用命令：Help 在線幫助Demo 例程演示Path 路徑設(shè)置Who 列出當(dāng)前變量What 列出當(dāng)前路徑下的.m文件Clear 清除所有變量Clc 清除當(dāng)前命令窗口命令Clf 清除當(dāng)前圖形Quit 退出MATLAB運(yùn)算符：加法： +減法： 乘法： *點(diǎn)乘： .*指數(shù)： . 除法： 點(diǎn)除

7、： .關(guān)系運(yùn)算符：小于 < 大于 > 小于等于 <= 大于等于 >= 等于 = 不等于 =基本矩陣zeros 零陣 ones 1陣 信號(hào)運(yùn)算abs取模值angle相角real實(shí)部imag虛部conv卷積fft快速傅立葉變換ifft 快速傅立葉反變換圖形函數(shù)figure生成一圖框axis設(shè)置坐標(biāo)軸hold保持原圖text在圖上標(biāo)記文字gtext通過鼠標(biāo)點(diǎn)擊方式在圖上標(biāo)記文字plot畫圖subplot 子圖title添加圖名xlabel給X軸增加文本標(biāo)注ylabel給Y軸增加文本標(biāo)注grid設(shè)置網(wǎng)格線stem離散序列數(shù)據(jù)桿狀圖bar條形圖polar極坐標(biāo)實(shí)驗(yàn)二 離散時(shí)間信號(hào)

8、分析實(shí)驗(yàn)?zāi)康模和ㄟ^Matlab實(shí)現(xiàn)基本信號(hào)的表示和運(yùn)算。理解信號(hào)的構(gòu)成和性質(zhì)。實(shí)驗(yàn)任務(wù)：1 學(xué)習(xí)使用基本信號(hào)在Matlab 中的基本表示方法。2 了解本實(shí)驗(yàn)中的擴(kuò)展函數(shù)的結(jié)構(gòu)和組成，實(shí)現(xiàn)自己設(shè)計(jì)的擴(kuò)展函數(shù)的編程： （1）寫一新的卷積函數(shù)conv_m，它可求出帶下標(biāo)的序列卷積。（2）修改給出的evenodd函數(shù)，使他能接受任意序列并把它分解成共軛對稱分量和共軛反對稱分量。3 掌握復(fù)雜信號(hào)的表示方法；并對復(fù)雜信號(hào)進(jìn)行編程實(shí)現(xiàn)。 （1）產(chǎn)生下列序列并繪出離散圖：（a） x(n)=2(n+2)-(n-4) -5n5（b） x(n)=nu(n)-u(n-10)+10e-0.3(n-10)u(n-10)-

9、u(n-2) 0n20（c） x(n)=cos (0.04n)+0.2w(n) 0n50 其中w(n)是均值為零、方差為1的白噪聲序列（d） （n）=x1 x1 x1 x1,其中x1=5 4 3 2 1 (e) 已知x(n)=u(n)-u(n-10),要求將它分解成奇偶序列。(2) 用Matlab產(chǎn)生一矩形脈沖串，連續(xù)時(shí)間為2s,脈沖個(gè)數(shù)為10，并繪制其波形圖。4 掌握復(fù)序列x(n)的實(shí)部、虛部、幅值和相位離散圖，并編程實(shí)現(xiàn)。 產(chǎn)生復(fù)信號(hào)： x(n)=e(-0.1+j0.3)n -10n10并畫出復(fù)序列x(n)的實(shí)部、虛部、幅值和相位離散圖。一離散信號(hào)在數(shù)字信號(hào)處理（DSP）中，所有的信號(hào)都是

10、離散（時(shí)間）信號(hào)，因此首先應(yīng)解決在MATLAB中如何表示離散信號(hào)。設(shè)一模擬信號(hào)經(jīng)A/D變換后，得到序列信號(hào)x(n)=x(n)= , x (-1), x ( 0 ) , x (1 ),由于MATLAB對下標(biāo)的約定為從1開始遞增，因此要表示x(n),一般應(yīng)采用兩個(gè)矢量，如:n=-3,-2,-1,0,1,2,3,4,5y=1,-1,3,2,0,4,5,2,1這表示了一個(gè)含9個(gè)采樣點(diǎn)的矢量：y(n)=x(-3), x (-2) , x (-1) , x (0) ,x (1 ) , x (5) 。通常情況下，序列值從x(0)開始，因此一個(gè)N點(diǎn)序列x(n)=x (0) ,x (1) , ,x (N-1)可

11、簡單地表示成y=x (0), x (1) , , x (N-1)這時(shí)y的下標(biāo)為1N。二基本信號(hào)表示 1單位取樣序列這一函數(shù)可利用MATLAB的zeros函數(shù)實(shí)現(xiàn)：x=zeros(1,N);x(1)=1 還可以借助于關(guān)系操作符實(shí)現(xiàn)：n=1：N;x=n=1;如要產(chǎn)生 則可采用MATLAB實(shí)現(xiàn)：n=n1：n2;x=(n-n0)=0;2單位階躍序列這一函數(shù)可利用MATLAB的ones函數(shù)實(shí)現(xiàn)：x=ones(1,N);還借助于關(guān)系操作符“>=”來實(shí)現(xiàn)。如要產(chǎn)生在n1n0n2上的單位階躍序列則可采用MATLAB實(shí)現(xiàn)：n=n1:n2;x=(n-n0)>=0;3實(shí)指數(shù)序列x(n)=an 采用MAT

12、LAB實(shí)現(xiàn)：n=0:N-1;x=a.n;4.復(fù)指數(shù)序列采用MATBLAB實(shí)現(xiàn)：n=0:N-1;x=exp(lu+j*w0)*n);5正（余）弦序列采用MATLAB實(shí)現(xiàn)：n=0:N-1;x=cos(w0*n+Q);6隨機(jī)序列MATLAB中提供了兩類（偽）隨機(jī)信號(hào)：rand(1,N)產(chǎn)生0，1上均勻分布的隨機(jī)矢量；randn(1,N)產(chǎn)生均值為0，方差為1的高斯隨機(jī)序列，也就是白噪聲序列。其它分布的隨機(jī)數(shù)可通過上述隨機(jī)數(shù)的變換而產(chǎn)生。7周期序列x(n)=x(n+N) 例如，設(shè)x1表示x序列中一個(gè)周期的序列，要產(chǎn)生P個(gè)周期的x序列，可以把它重復(fù)4次，用MALAB實(shí)現(xiàn)：x=x1,x1,x1,.,x1;

13、但是高明的方法是利用MATLAB的強(qiáng)有力的下標(biāo)能力。先產(chǎn)生一個(gè)包含P行x(n)值的矩陣，然后用結(jié)構(gòu)（：）把它的P行串接起來成為一個(gè)長行，不過這種結(jié)構(gòu)只能用于列向，所以我們往往還要用矩陣轉(zhuǎn)置把它擴(kuò)展成行向。>>xtilde=x1*ones(1,P); %P列x1>>xtilde=xtilde（：）； %長的列向量>>xtilde=xtilde; 長的行向量三 序列操作1信號(hào)加x(n)=x1(n)+x2(n)采用MATLAB實(shí)現(xiàn)：x=x1+x2;注意：當(dāng)x1和x2序列的長度不等或其位置不對應(yīng)時(shí)，信號(hào)相加就不是這么簡單。首先應(yīng)使x1、x2具有相等的長度，然后兩者對

14、齊，最后進(jìn)行相加。2信號(hào)乘x(n)=x1(n)*x2(n)這是一種樣本對樣本的相乘，也即點(diǎn)乘運(yùn)算，在MATLAB中可采用.*（數(shù)組乘法）來實(shí)現(xiàn)，但兩序列x1、x2也應(yīng)經(jīng)過處理。3改變比例y(n)=ax(n)=ax(n)采用MATLAB實(shí)現(xiàn)：y=alpha*x;4移位y(n)=x(n-k)5折疊y(n)=x(-n)它將序列x(n)在n=0處倒轉(zhuǎn)，在MATLAB中可直接用fliplr函數(shù)實(shí)現(xiàn)。6取樣和它不同于信號(hào)加。采用MATLAB實(shí)現(xiàn)：y=sum(x(n1:n2);7取樣積采用MATLAB實(shí)現(xiàn)：y=prod(x(n1:n2);8卷積和 y(n)=x(n)*h(n)=可直接采用MATLAB中的函數(shù)

15、conv，即y=conv(x,h)；它默認(rèn)序列從n=0開始。但如果序列是從一負(fù)值開始，即如x(n):nxbnnxeh(n):nhbnnhe其中nxb<0或nhb<0，或兩者同時(shí)為負(fù)時(shí)，就不能直接采用conv函數(shù)。通過卷積結(jié)果序列 y(n):nybnnye且 nyb=nxb+nhb , nye=nxe+nhe這樣我們就可以構(gòu)成一新的卷積函數(shù)conv_m，它可求出帶下標(biāo)的序列卷積。9信號(hào)能量采用MATLAB實(shí)現(xiàn)：Px=sum(abs(x).2); Px=sum(x.*conj(x);10奇偶綜合任何一個(gè)序列x(n)都可分解成偶對稱部分xe(n)和奇對稱部分x0(n)，即x(n)=xe(

16、n)+x0(n)其中xe(n)=x(n)+x(-n)x0(n)=x(n)-x(-n)這樣我們可以設(shè)計(jì)一函數(shù)evenodd，完成將任一給定序列x(n) 分解成xe(n)和 x0(n)。四MATLAB信號(hào)生成函數(shù)1square（）方波發(fā)生器。調(diào)用格式為：square(T) 返回一周期為2的方波，采樣時(shí)刻由向量T指定。square(T,DUTY) 產(chǎn)生一給定占空比的方波，占空比DUTY是信號(hào)為正值的比例。2 sawtooth（）鋸齒波和方波發(fā)生器。調(diào)用格式為：sawtooth(T) 返回一周期為2的鋸齒波，采樣時(shí)刻由向量T指定。sawtooth(T,WIDTH) 產(chǎn)生一三角波，WIDTH指定最大值出

17、現(xiàn)的地方，其取值在0到1之間，1對應(yīng)2。3 diric（）周期sinc函數(shù)（Dirichlet函數(shù)）。調(diào)用格式：Y=diric(X,N) 返回一大小與X相同的矩陣，其元素為Dirichlet函數(shù)。N必須為正整數(shù)，該函數(shù)將0到2等間隔分成N等分。4 pulstran （）通過對連續(xù)函數(shù)或脈沖原型進(jìn)行采樣而得到脈沖序列的發(fā)生器5 rectplus（）非周期矩形波發(fā)生器6 tripuls（）非周期三角脈沖發(fā)生器五MATLAB七個(gè)基本的擴(kuò)展函數(shù)1單位取樣序列impseq.mfunction x,n=impseq(n0,n1,n2)%Generates x(n)=delta(n-n0);n1<=n

18、,n0<-n2%-%x,n=impseq(n0,n1,n2)%if(n0<n1)|(n0>n2)|(n1>n2) error(arguments must satisfy n1<=n0<=n2)end n=n1:n2;x=(n-n0)=0;2單位階躍序列stepseq.mfunctionx,n=stepseq(n0,n1,n2)%Generates x(n)=u(n0);n1<=n0<=n2%-%x,n=stepseq(n0,n1,n2)%if(n0<n1)|(n0>n2)|(n1>n2) error(arguments mus

19、t satisfy n1<=n0<=n2)end n=n1:n2;x=(n-n0)>=0;3信號(hào)加sigadd.mfunctiony,n=sigadd(x1,n1,x2,n2)%implements y(n)=x1(n)+x2(n)%-%y,n=sigadd(x1,n1,x2,n2)%y=sum sequence over n, which includes n1 and n2%x1=first sequence over n1%x2=second sequence over n2(n2 can be different from n1)%n=min(min(n1),min(

20、n2):max(max(n1),max(n2);y1=zeros(1,length(n);y2=y1;y1(find(n>=min(n1)&(n<max(n1)= =1)=x1;y2(find(n>=min(n2)&(n<=max(n2)= =1)=x2;y=y1+y24信號(hào)乘sigmult.mfunctiony,n=sigmult(x1,n1,x2,n2)%implements y(n)=x1(n)*x2(n)%-%y,n=sigmult(x1,n1,x2,n2)%y=product sequence over n ,which includes n1

21、 and n2%x1=first sequence over n1%x2=second sequence over n2(n2 can be different from n1)%n=min(min(n1),min(n2):max(max(n1),max(n2);y1=zeros(1,length(n);y2=y1;y1(find (n>=min(n1)&(n<=max(n1)= =1)=x1;y2(find(n>=min(n2)&(n<=max(n2)= =1)=x2;y=y1.*y25移位sigshift.m functiony,n=sigshift

22、(x,m,n0)%implements y(n)=x(n-n0)%-%y,n=sigshift(x,m,n0)%n=m+n0;y=x;6折疊sigfold.m functiony,n=sigfold(x,n)%implements y(n)=x(-n)%-%y,n=sigsfold(x,n)%y=fliplr(x);n=-fliplr(n);7 奇偶綜合evenodd.mfunction xe,xo,m=evenodd(x,n)%Real signal decomposition into even and odd parts%-%xe,xo,m=evenodd(x,n)%if any (im

23、ag(x)=0) error(x is not a real sequence)end m=-fliplr(n);m1=min(m,n);m2=max(m,n);m=m1:m2;nm=n(1)-m(1);n1=1:length(n);x1=zeros(1,length(m);x1(n1+nm)=x;x=x1;xe=0.5*(x+fliplr(x);xo=0.5*(x-fliplr(x);【例2-1】實(shí)驗(yàn)任務(wù)3.（1）產(chǎn)生下列序列并繪出離散圖：x(n)=2(n+2)-(n-4) -5n5解： 利用MATLAB及信號(hào)處理工具箱函數(shù)，再加上前面構(gòu)造的幾個(gè)函數(shù)如sigadd,sigmult等，可很容易

24、編寫出可直接執(zhí)行的MATLAB程序ex31.m:%Example 3.1%a) x(n)=2*delta(n+2)-delta(n-4), -5<=n<=5%figure(1);clfn=-5:5;x=2* impseq (-2,-5,5)-impseq(4,-5,5);subplot(2,2,1);stem(n,x);title (Sequence in Example 3.1a);【例2-2】 設(shè)線性時(shí)不變（LTI）系統(tǒng)的單位沖激響應(yīng)為 h(n)=(0.9)nu(n)輸入序列為 x(n)=u(n)-u(n-10) 求系統(tǒng)輸出y(n)。解： 系統(tǒng)輸出y(n)為輸入x(n)與系統(tǒng)單

25、位沖激響應(yīng)h(n)的卷積，可直接采用conv-m函數(shù)求出輸出序列。MATLAB程序?yàn)閑x22.m:%Example2-2%x(n)=u(n)-u(n-10);h(n)=(0.9)n*u(n) n=0:04%y(n)=conv(x,h)%n=-5:50;u1=stepseq(0,-5,50);u2=stepseq(10,-5,50);%input x(n)x=u1-u2;%impulse response h(n)h=(0.9).n).*u1; figure(1)subplot (3,1,1);stem(n,x);axis(-5,50,0,2)title(Input Sequence)ylabe

26、l(x(n)subplot(3,1,2);stem(n,h);axis(-5,50,0,2)ylabel(h(n)%output responsey,ny=conv_m(x,n,h,n);subplot(3,1,3);stem(ny,y);axis(-5,50,0,8)gtext(Impulse Response)xlabel(n),ylabel(y(n)執(zhí)行結(jié)果分析所得的序列的離散圖.實(shí)驗(yàn)報(bào)告要求1. 實(shí)驗(yàn)報(bào)告包括目的、要求、內(nèi)容、步驟、結(jié)果、總結(jié)，形成完整實(shí)驗(yàn)報(bào)告。實(shí)驗(yàn)步驟不是書本內(nèi)容的復(fù)制，而是自己結(jié)合實(shí)驗(yàn)內(nèi)容進(jìn)行探索的過程。2. 提交打印版實(shí)驗(yàn)報(bào)告，A4紙張打印，并附源程序清單。3.

27、提交的報(bào)告和源程序中標(biāo)注清楚姓名學(xué)號(hào)專業(yè)等基本信息，禁止抄襲。實(shí)驗(yàn)三 離散系統(tǒng)分析實(shí)驗(yàn)?zāi)康模和ㄟ^用Matlab實(shí)現(xiàn)信號(hào)的分析和表示方法，掌握離散系統(tǒng)的分析方法。實(shí)驗(yàn)任務(wù)：1. 了解系統(tǒng)的Matlab描述和轉(zhuǎn)換2. 掌握分析線性時(shí)不變系統(tǒng)的方法，并編程實(shí)現(xiàn)。給定系統(tǒng)H(z)=-0.2z/(z2+0.8),(1) 求出H(z)的幅頻響應(yīng)和相頻響應(yīng)；(2) 繪制極零點(diǎn)圖；(3) 求出并繪出該系統(tǒng)的單位抽樣響應(yīng)；(4) 令x(n)=u(n)，求出并繪出系統(tǒng)的單位階躍響應(yīng)y(n)。MATLAB系統(tǒng)的描述和轉(zhuǎn)換MATLAB的信號(hào)處理工具提供了幾種線性時(shí)不變系統(tǒng)的模型，我們可以靈活的選擇合適的模型以滿足不同

28、應(yīng)用場合的需要。信號(hào)處理工具箱提供的離散時(shí)間系統(tǒng)模型常用于數(shù)字濾波器，這些模型包括傳遞函數(shù)、零極點(diǎn)增益形式、狀態(tài)空間形式、部分分式、二階級(jí)聯(lián)形式、格狀結(jié)構(gòu)和卷積矩陣等。1 傳遞函數(shù)形式在MATLAB信號(hào)處理工具箱中，數(shù)字濾波器的傳遞函數(shù)具有以下形式：這種形式稱為DSP（Digital Signal Processing）,它和控制工具箱的z傳遞函數(shù)模型形式不同。這一模型用z傳遞函數(shù)分子和分母多項(xiàng)式的系數(shù)構(gòu)成兩個(gè)向量來確定，即 num=b(1) b(2)b(nb+1)den=a(1) a(2)a(na+1)式中nb和na分別為分子num和den的階數(shù)在MATLAB中，要定義一個(gè)DSP形式的傳遞函

29、數(shù)可利用FILT。調(diào)用格式： sys=filt(num,den) sys=filt(num,den,Ts)2 零極點(diǎn)增益形式 在MATLAB信號(hào)處理工具箱中，離散時(shí)間系統(tǒng)（數(shù)字濾波器）的零極點(diǎn)形式具有以下形式：這一模型用零點(diǎn)向量和極點(diǎn)向量及增益來表示，即 q(z)=q(1) ,q(2),q(n)p(z)=p(1) ,p(2),p(n)MATLAB函數(shù)poly和roots可實(shí)現(xiàn)傳遞函數(shù)形式和零極點(diǎn)增益形式的轉(zhuǎn)換。3 狀態(tài)空間形式離散系統(tǒng)的狀態(tài)空間表達(dá)式為： x(n+1)=Ax(n)+Bu(n) y(n)=Cx(n)+Du(n)式中，u為輸入，y為輸出,x為狀態(tài)向量。在MATLAB信號(hào)處理工具箱中

30、，用矩陣A、B、C、D表示數(shù)字濾波器模型（和控制工具箱一樣）MATLAB信號(hào)處理工具箱還提供傳遞函數(shù)、零極點(diǎn)和狀態(tài)空間三種模型之間的相互轉(zhuǎn)換函數(shù)，它們是zp2ft、zp2ss 、tf2ss、 tf2zp、ss2tf和 ss2zp4 部分分式離散時(shí)間系統(tǒng)（數(shù)字濾波器）的傳遞函數(shù)H(z)也可以用下面部分分式展開來描述：（1）若H(z)部包含重極點(diǎn)（2）若H(z)包含多重極點(diǎn)，如p(j)有Sr個(gè)重極點(diǎn)，H(z)部分分式展開式中包含相應(yīng)的形式為：MATLAB信號(hào)處理工具箱中，函數(shù)RESIDUEZ()用來把傳遞函數(shù)形式轉(zhuǎn)變?yōu)椴糠址质秸归_式。調(diào)用格式為：r,p,k=residuez(b,a) 其中，b,a

31、分別為z傳遞函數(shù)的分子、分母；p為傳遞函數(shù)極點(diǎn)向量；r為相應(yīng)極點(diǎn)的留數(shù)向量；k為傳遞函數(shù)商向量。若原傳遞函數(shù)為真分?jǐn)?shù)，此項(xiàng)為空向量。Residuez也可用于將部分分式展開式轉(zhuǎn)換為傳遞函數(shù)形式，調(diào)用格式為：b,a=residuez(r,p,k) 可用函數(shù)residuez將H(z)展成部分分式的形式，然后利用公式直接求H(z)的z反變換。5 二階級(jí)聯(lián)形式任何傳遞函數(shù)可寫為二階因子級(jí)聯(lián)形式，即即H(z)可以化成多個(gè)二階因式相乘，L為二階因式的個(gè)數(shù)。這樣一個(gè)離散系統(tǒng)可用一個(gè)L×6的數(shù)組SOS表示，每一行表示一個(gè)二階因子，分子三個(gè)系數(shù)和分母三個(gè)系數(shù)。 這種形式用于濾波器的二階節(jié)級(jí)聯(lián)結(jié)構(gòu)中。MA

32、TLAB中提供了一般模型和SOS形式的相互轉(zhuǎn)換函數(shù)zp2sos、sos2zp、sos2os、sos2tf、6 格狀結(jié)構(gòu)濾波器的實(shí)現(xiàn)通常采用格狀結(jié)構(gòu)。對于N階全零點(diǎn)和全極點(diǎn)濾波器采用格型結(jié)構(gòu)，可以用多項(xiàng)式a(n),n=1,2,N來描述，也可以通過格型結(jié)構(gòu)反射系數(shù)k(n)，n=1,2,.,N來描述。 MATLAB信號(hào)處理工具箱提供了濾波器格狀結(jié)構(gòu)和轉(zhuǎn)換函數(shù)latcfilt、poly2rc、rc2poly、latc2tf和tf2latc。7 卷積矩陣若兩個(gè)向量x和h，滿足下面等式的矩陣C稱為卷積矩陣。C·x=h*x （式36）當(dāng)濾波器的單位采樣序列為h(n)，濾波器輸入為x(n)，濾波器輸

33、出信號(hào)為： y(n)=h*x=C·x可見，卷積矩陣是數(shù)字濾波器的另一種描述形式。MATLAB信號(hào)處理工具箱提供了計(jì)算卷積矩陣函數(shù)convmtx。調(diào)用格式為： C=convmtx(h,n)式中，h為m×1列向量；n為x列向量維數(shù)；c為(m+n-1)×n階卷積矩陣。【例31】已知濾波器的傳遞函數(shù) 試求該濾波器的零極點(diǎn)和狀態(tài)空間表達(dá)式和部分分式展開式。解：用MATLAB編程如下%MATLAB Program 3-1%input parametersnum=2 3;den=1 0.4 1;f=filt(num,den)disp(Convert to zero-pole-g

34、ain)z,p,k=tf2zp(num,den)disp(Convert to state space)A,B,C,D=tf2ss(num,den)transform partial-fraction expansionr,p1,k1=residuez(num,den)MATLAB系統(tǒng)分析MATLAB信號(hào)處理工具箱提供了許多工具函數(shù)來分析濾波器（系統(tǒng)）的特性，包括時(shí)域分析，對任意輸入響應(yīng)、脈沖響應(yīng)等；頻域分析，幅值響應(yīng)、相位響應(yīng)、零極點(diǎn)位置；其他特性，群延遲等。濾波器時(shí)域和頻域分析是設(shè)計(jì)各類濾波器、評價(jià)濾波器性能和濾波器應(yīng)用基礎(chǔ)。這里簡單介紹一下有關(guān)函數(shù)。一時(shí)間響應(yīng)分析工具filter,fil

35、tic,fftfilt和impulse1 Filter（）函數(shù)filter用于實(shí)現(xiàn)IIR和FIR濾波器對數(shù)據(jù)濾波，常用來計(jì)算濾波器對輸入的響應(yīng)，函數(shù)調(diào)用格式為：y=filter(b,a,x)其中，b，a分別是濾波器系統(tǒng)函數(shù)H(z)的分子和分母多項(xiàng)式系數(shù)向量；x為濾波器輸入向量或矩陣； 函數(shù)filter可用于全極點(diǎn)、全零點(diǎn)、零極點(diǎn)系統(tǒng)。2 Filtic（）函數(shù)filtic用于從系統(tǒng)過去值y和x求系統(tǒng)狀態(tài)初始值，調(diào)用格式：z=filtic(b,a,y,x)z=filtic(b,a,y)其中，y=y(n-1),y(n-2),.,y(-na)為輸出y的過去向量；x=x(n-1),x(n-2),.,x(

36、-nb)為輸出x的過去向量;b,a分別為分子、分母多項(xiàng)式系數(shù)向量；z為系統(tǒng)初始狀態(tài)。如果，n<=-1時(shí)，x(n)=0,則filtic函數(shù)中不用指定x?！纠?2】求解 其中 ，初始條件為y(-1)=4和y(-2)=10。解：用MATLAB計(jì)算全響應(yīng)：n=0:7; x=(1/4).n; xic=1,2;b=1;a=1,1.5,0.5;format longy1=filter(b,a,x,xic);Matlab提供了filtic函數(shù)計(jì)算對于上例：Y=4,10;xic=filter(b,a,Y)3Fftfilt（）函數(shù)fftfilt利用效率高的基于FFT重疊相加法實(shí)現(xiàn)對數(shù)據(jù)濾波，該函數(shù)只適用于F

37、IR濾波器，函數(shù)調(diào)用格式為：y=fftfilt(b,x)y=fftfilt(b,x,n)式中，b為FIR濾波器系數(shù)向量；x為輸入數(shù)據(jù)；n為FFT長度，省時(shí)，函數(shù)選擇最佳的FFT長度；y為濾波器輸出。Y=fftfilt(b,x)等價(jià)于y=filter(b,1,x)。【例33】：FIR低通濾波器截止頻率為200Hz，采樣頻率Fs=1000Hz，對信號(hào)x(t)=sin(2f1t)+sin(2f2t)濾波，f1=50Hz,f2=250Hz,求濾波器輸出。%Matlab Program 3.2%Create Orignal signal dataN=1000;Fs=1000;n=0：N-1;t=n/.F

38、s;x=sin(2*pi*50*t)+sin(2*pi*240*t);%Create a FIR filter with lawpassb=fir1(40,200/500);%Filtering the datayfft=fftfilt(b,x,256);Outputn1=81:241;t1=t(n1);x1=x(n1);y1=yfft(n1);subplot(221);plot(t1,x1);tilte(Original Signal);subplot(222);plot(t1,y1);title(Signal after the filter);grid;執(zhí)行上面這個(gè)程序，分析程序結(jié)果。4

39、Impz（）impz用于產(chǎn)生數(shù)字系統(tǒng)的脈沖響應(yīng)?；菊{(diào)用格式為：h,t=impz(b,a,n,Fs)b,a分別為系統(tǒng)傳遞函數(shù)分子和分母多項(xiàng)式系數(shù)向量；n為采樣點(diǎn)數(shù)；Fs為采樣周期，缺省值1；h為單位脈沖響應(yīng)向量；t為和h對應(yīng)的時(shí)間向量。當(dāng)函數(shù)輸出缺省時(shí)，繪制系統(tǒng)脈沖響應(yīng)圖；當(dāng)n缺省時(shí)，函數(shù)自動(dòng)選擇n值。在Matlab中，h=impz(b,a,n)等價(jià)于下面語句：x=1 zeros(1,n);h=filter(b,a,x);因此，不少場合下，用后一種方法求系統(tǒng)的脈沖響應(yīng)更為靈活。5 Conv（）線性時(shí)不變系統(tǒng)的輸入輸出關(guān)系還可通過單位脈沖響應(yīng)h(n)表示：y(n)=x(n)*h(n)=所以可直接

40、采用MATLAB中的函數(shù)conv或者用可求出帶下標(biāo)的序列卷積的函數(shù)conv_m來求系統(tǒng)對輸入的響應(yīng)。二 頻率響應(yīng)1 Freqs()函數(shù)freqs用于求模擬濾波器的頻率響應(yīng)。調(diào)用格式為h=freqs(b,a,)； h,=freqs(b,a,n); freqs(b,a)其中，b,a分別為模擬濾波器傳遞函數(shù)分子和分母多項(xiàng)式向量；n為頻率點(diǎn)數(shù)(整數(shù)),缺省200；h為頻率響應(yīng),復(fù)數(shù)；頻率向量，實(shí)數(shù)。函數(shù)缺省時(shí)，繪制模擬濾波器的幅頻響應(yīng)圖和相頻響應(yīng)圖。模擬濾波器的系統(tǒng)函數(shù)形式為輸出h為模擬濾波器的頻率響應(yīng)H(ej)。2Freqz（）函數(shù)freqz用于求數(shù)字濾波器的頻率響應(yīng)。調(diào)用格式為h,=freqz(b

41、,a,n); h,f=freqz(b,a,n,Fs);h,=freqz(b,a,n,whole); h,f=freqz(b,a,n,whole,Fs);h=freqz(b,a,f,Fs); freqz(b,a);說明：（1）函數(shù)輸入：b,a分別為數(shù)字濾波器傳遞函數(shù)分子和分母多項(xiàng)式向量；n為頻率點(diǎn)數(shù)(整數(shù)),最好為2的冪，缺省512；Fs為采樣頻率，Hz；f為給定的頻率矢量；whole 表示返回的值為包含z整個(gè)單位圓頻率矢量，即02；缺省時(shí)，僅為包含z平面上半單位圓（0）之間等間距N個(gè)點(diǎn)頻率矢量。函數(shù)返回：h為頻率響應(yīng),復(fù)數(shù)；為n點(diǎn)頻率向量（弧度），返回范圍與輸入?yún)⒘縲hole有關(guān)；f為n點(diǎn)頻率

42、向量(Hz)。函數(shù)缺省時(shí)，繪制模擬濾波器的幅頻響應(yīng)圖和相頻響應(yīng)圖。（2）該函數(shù)適用于和式31形式的數(shù)字濾波器。（3）該函數(shù)執(zhí)行補(bǔ)零的FFT算法，即h=fft(b,n)./fft(a,n)（4）Matlab提供兩個(gè)函數(shù)ABS和ANGLE由頻率響應(yīng)H(ej)求幅頻響應(yīng)| H(ej)|和相頻響應(yīng)argH(ej).（5）函數(shù)freqz輸出的頻率向量在02之間，為了獲得一個(gè)濾波器真正的相頻特性圖，要對相位角進(jìn)行修正。為此Matlab信號(hào)處理工具箱提供一個(gè)函數(shù)UNWRAP.3 UNWRAP.函數(shù)UNWARP展開由函數(shù)freqz產(chǎn)生的頻率.調(diào)用格式為：P=unwrap()三 零極點(diǎn)圖系統(tǒng)零極點(diǎn)位置決定了系統(tǒng)

43、穩(wěn)定性和性能，因此，考察系統(tǒng)零極點(diǎn)位置是分析系統(tǒng)特性重要方面之一。Matlab信號(hào)處理工具箱提供繪制線性離散時(shí)間系統(tǒng)零極點(diǎn)圖的工具函數(shù)zplane.調(diào)用格式為： zplane(z,p) zplane(b,a)其中，z,p分別為系統(tǒng)零、極點(diǎn)向量；b,a為系統(tǒng)z傳遞函數(shù)分子、分母多項(xiàng)式系數(shù)向量。四 群延遲 由信號(hào)傳輸不失真條件，濾波器相頻特性應(yīng)該是一條經(jīng)過原點(diǎn)直線，即 ()= td td為常數(shù) （式38）但一般實(shí)際系統(tǒng)不滿足上述條件，衡量實(shí)際濾波器相位平均延遲用群延遲。群延遲定義為信號(hào)通過系統(tǒng)的平均延遲對頻率的函數(shù)，即濾波器相頻特性圖上切線的負(fù)斜率 （式39）Matlab信號(hào)處理工具箱提供計(jì)算群延

44、遲函數(shù)GRPDELAY，調(diào)用格式為：gd,=grpdelay(b,a,n); gd,f=grpdelay(b,a,n,Fs);gd,=grpdelay(b,a,n,whole);gd,f=grpdelay(b,a,n,wholeFs);其中，gd為群延遲，其他各項(xiàng)意義同函數(shù)freqz。函數(shù)輸出項(xiàng)缺省時(shí)繪制群延遲圖?！纠?4】：編一函數(shù)實(shí)現(xiàn)：求系統(tǒng)的絕對的和相對的幅度響應(yīng)，相位響應(yīng)及群延遲響應(yīng)。解：函數(shù)freqz_m是對freqz的修正，它給出了系統(tǒng)的絕對的和相對的幅度響應(yīng)，相位響應(yīng)及群延遲響應(yīng)。程序如下：function db,mag,pha,grd,w=freqz_m(b,a)% db,ma

45、g,pha,grd,w=freqz_m(b,a)%db=0pi區(qū)間內(nèi)的相對振幅mag=0pi區(qū)間內(nèi)的絕對振幅pha=0pi區(qū)間內(nèi)的相位響應(yīng)grd=0pi區(qū)間內(nèi)的群延遲w=0pi區(qū)間內(nèi)的501個(gè)頻率樣本向量b=H(z)的分子多項(xiàng)式的系數(shù)（對FIR：b=h）a=H(z)的分母多項(xiàng)式的系數(shù)(對FIR:a=1)H,w=freqz(b,a,1000,whole);H=(H(1:1:501);w=(w(1:1:501);mag=abs(H);db=20*log10(mag+eps)/max(mag);pha=angle(H); grd=grpdelay(b,a,w);【例35】 已知濾波器的傳遞函數(shù) H(

46、z)=0.15/(1-0.8z-1)輸入信號(hào)為x(t)=2sin(0.05t)+w(t)，w(t)為隨機(jī)信號(hào)，幅值為0.2。試?yán)L出濾波器的輸出信號(hào)波形。用兩種方法求濾波器的輸出，MATLAB編寫程序如下：%MATLAB Program 3-3b=0.15;a=1 0.8;%create input of filter n=0:100;x=2*sin(0.05*pi*n)+0.2*randn(1,101);%Calculate response according to convolutionimp=1,zeros(100,1);h=impz(b,a,101);yc=conv(h,x);yc=y

47、c(1:101);%Calculate response according to matlab functiony1=filter(b,a,x);plot(n,x,r,n,y1,b,n,y,m);xlabel(n);ylabel(x y yc);title(Time response);grid;分析程序執(zhí)行結(jié)果。實(shí)驗(yàn)報(bào)告要求1. 實(shí)驗(yàn)報(bào)告包括目的、要求、內(nèi)容、步驟、結(jié)果、總結(jié)，形成完整實(shí)驗(yàn)報(bào)告。實(shí)驗(yàn)步驟不是書本內(nèi)容的復(fù)制，而是自己結(jié)合實(shí)驗(yàn)內(nèi)容進(jìn)行探索的過程。2. 提交打印版實(shí)驗(yàn)報(bào)告，A4紙張打印，并附源程序清單。3. 提交的報(bào)告和源程序中標(biāo)注清楚姓名學(xué)號(hào)專業(yè)等基本信息，禁止抄襲。實(shí)驗(yàn)四 離

48、散傅立葉變換實(shí)驗(yàn)?zāi)康模豪斫怆x散傅立葉變換及其性質(zhì)，與線性卷積和線性移位相比較掌握循環(huán)卷積和循環(huán)移位的實(shí)現(xiàn)方法。實(shí)驗(yàn)任務(wù)：1. 根據(jù)給出的函數(shù)，編程實(shí)現(xiàn)某一序列DFT, IDFT。已知x1為5點(diǎn)序列，x2為14點(diǎn)序列，各作15點(diǎn)的DFT，然后兩個(gè)DFT相乘，再求乘積的IDFT，設(shè)結(jié)果為y(n)，分析y(n)哪些點(diǎn)對應(yīng)于x1*x2應(yīng)該得到的點(diǎn)？2. 編程實(shí)現(xiàn)某一序列的循環(huán)移位和循環(huán)卷積。已知x1=(n+1)R4(n)，h(n)=R4(n-2)，試求x1和h(n)的循環(huán)卷積，并作圖。3. 線性卷積與循環(huán)卷積的關(guān)系已知x1=1 2 2，x2=1 2 3 4，試計(jì)算并分析結(jié)果（1）x1x2；(2)x1x

49、2；(3) x1x2；(4)x1*x2.4DFT的奇偶虛實(shí)性已知序列x(n)=10(0.8)n，序列長度N=21,繪出x(n)傅立葉變換的實(shí)部和虛部并驗(yàn)證其奇偶性。4. 執(zhí)行exam41例程序，分析其結(jié)果?！纠?41】為說明高密度頻譜和高分辯頻譜之間的差異，考慮 x(n)=cos(0.48n)+cos(0.52n) （1）取x(n)(0n10)時(shí)，求x(n)的DFT X(k)； （2）將（1）中的x(n)以補(bǔ)零方式使x(n)加長到0n100,求X(k) （3）取x(n)(0n100),求X(k).解 MATLAB程序?yàn)?exam41.m:%MATLAB Program 4-1%part 1%s

50、pecttum based on the first 10 amples of x(n)n=0:1:10;x=cos(0.48*pi*n)+cos(0.52*pi*n);figure(1)n1=0:1:9;y1=x(1:1:10);subplot(2,1,1);stem(n1,y1);title('signal x(n),0<=n<=9');xlabel('n')%axis(0,10,-4,4);Y1=fft(y1);magY1=abs(Y1(1:1:6);k1=0:1:5;w1=2*pi/10*k1;Subplot(2,1,2);stem(w1/p

51、i,magY1);title('Samples of DTET Magnitude');Xlabel('frequency in pi units')axis(0,1,0,10)%part 2%High density spectrum (100 samples )based on the first 10 samples of x(n)figure(2)n3=0:1:99;y3=x(1:1:10)zeros(1,90);subplot(2,1,1);stem(n3,y3);title('signal x(n),0<=n<=9+90 zeros');xlabel('n')axis(0,100,-2.5,2.5);Y3=fft(y3);magY3=abs(Y3(1:1:51);k3=0:1:50;w3=2*pi/100*k3;subplot(2,1,2);stem(w3/pi,magY3