理論課第7章yyxMATLAB在數(shù)字信號(hào)處理中的應(yīng)用精講課件

1第7章

MATLAB在數(shù)字信號(hào)處理中的應(yīng)用1第7章

MATLAB在數(shù)字信號(hào)處理中的應(yīng)用2書(shū)山有路勤為徑，學(xué)海無(wú)涯苦作舟。好好學(xué)習(xí)，天天向上2書(shū)山有路勤為徑，好好學(xué)習(xí)，天天向上3引言數(shù)字信號(hào)處理中包含離散時(shí)間信號(hào)與系統(tǒng)的一部分知識(shí)。本章主要研究?jī)?nèi)容分為兩部分：離散時(shí)間信號(hào)與系統(tǒng)時(shí)域分析離散時(shí)間信號(hào)與系統(tǒng)頻域分析在數(shù)字信號(hào)處理中還有一個(gè)非常大的知識(shí)點(diǎn)：數(shù)字濾波器設(shè)計(jì)，留待下學(xué)期開(kāi)設(shè)本課程時(shí)，可以學(xué)習(xí)。3引言數(shù)字信號(hào)處理中包含離散時(shí)間信號(hào)與系統(tǒng)的一部分知識(shí)。4一、離散時(shí)間信號(hào)與系統(tǒng)時(shí)域分析4一、離散時(shí)間信號(hào)與系統(tǒng)時(shí)域分析5包含內(nèi)容信號(hào)分析典型信號(hào)表示impseq(),stepseq()7.1序列相加，相乘7.2序列合成和截取7.3序列移位和周期延拓7.6離散序列卷積系統(tǒng)分析7.4系統(tǒng)響應(yīng)7.5系統(tǒng)線性判定5包含內(nèi)容信號(hào)分析67.1離散信號(hào)的產(chǎn)生及時(shí)域處理時(shí)域離散信號(hào)用x(n)表示（離散序列）。x和n同時(shí)使用才能完整地表示一個(gè)序列：時(shí)間變量n（表示采樣位置）只能取整數(shù)；向量x（表示采樣點(diǎn)的幅度）與n一一對(duì)應(yīng)。由于n序列是按整數(shù)遞增的,可簡(jiǎn)單地用其初值ns決定,因?yàn)樗慕K值nf取決于ns和x的長(zhǎng)度length(x),故可寫(xiě)成：n=[ns:nf]或 n=[ns:nslength(x)1]67.1離散信號(hào)的產(chǎn)生及時(shí)域處理時(shí)域離散信號(hào)用x(n)表7典型信號(hào)表示impseq(),stepseq()7典型信號(hào)表示impseq(),stepseq()8例7.1序列的相加和相乘給出兩個(gè)序列x1(n)和x2(n)。x1=[0,1,2,3,4,3,2,1,0];n1=[-2:6];x2=[2,2,0,0,0,-2,-2];n2=[2:8];要求它們的和ya及乘積yp。解：編程的思路是把序列長(zhǎng)度延拓到覆蓋n1和n2的范圍，這樣才能把兩序列的時(shí)間變量對(duì)應(yīng)起來(lái)，然后進(jìn)行對(duì)應(yīng)元素的運(yùn)算。8例7.1序列的相加和相乘給出兩個(gè)序列x1(n)和x9例程x1=[0,1,2,3,4,3,2,1,0];ns1=-2; %給定x1及ns1x2=[2,2,0,0,0,-2,-2];ns2=2; %給定x2及ns2nf1=ns1+length(x1)-1;nf2=ns2+length(x2)-1;ny=min(ns1,ns2):max(nf1,nf2);%y(n)的時(shí)間變量xa1=zeros(1,length(ny));xa2=xa1;%延拓序列初始化xa1(find((ny>=ns1)&(ny<=nf1)==1))=x1;%給xa1賦值x1xa2(find((ny>=ns2)&(ny<=nf2)==1))=x2;%給xa2賦值x2ya=xa1+xa2%序列相加yp=xa1.*xa29例程x1=[0,1,2,3,4,3,2,1,0];ns1=10結(jié)果顯示10結(jié)果顯示11例7.2序列的合成和截取用例6.13的結(jié)果編寫(xiě)產(chǎn)生矩形序列RN(n)的程序。序列起點(diǎn)為n0，矩形序列起點(diǎn)為n1，長(zhǎng)度為N（n0，n1，N由鍵盤(pán)輸入）。并用它截取一個(gè)復(fù)正弦序列exp(jπn/8).解：建模:矩形序列可看成兩個(gè)階躍序列之差。用MATLAB邏輯關(guān)系產(chǎn)生矩形序列x2(n)。而用它截取任何序列相當(dāng)于元素群相乘x2.*x，也稱為加窗運(yùn)算。序列的合成和截取就是相加和相乘。11例7.2序列的合成和截取用例6.13的結(jié)果編寫(xiě)12結(jié)果顯示12結(jié)果顯示13例7.2序列的合成和截取用例6.13的結(jié)果編寫(xiě)產(chǎn)生矩形序列RN(n)的程序。序列起點(diǎn)為n0，矩形序列起點(diǎn)為n1，長(zhǎng)度為N（n0，n1，N由鍵盤(pán)輸入）。并用它截取一個(gè)復(fù)正弦序列exp(jπn/8).解：建模:矩形序列可看成兩個(gè)階躍序列之差。用MATLAB邏輯關(guān)系產(chǎn)生矩形序列x2(n)。而用它截取任何序列相當(dāng)于元素群相乘x2.*x，也稱為加窗運(yùn)算。序列的合成和截取就是相加和相乘。13例7.2序列的合成和截取用例6.13的結(jié)果編寫(xiě)14程序要點(diǎn)n=n0:n1+N+5; %生成自變量數(shù)組u=[(n-n1)>=0];%產(chǎn)生單位階躍序列（u(n-n1)）x1=[(n-n1)>=0]-[(n-n1-N)>=0] %用階躍序列之差產(chǎn)生矩形序列x2=[(n>=n1)&(n<(N+n1))]; %用邏輯式產(chǎn)生矩形序列x3=exp(j*n*pi/8).*x2; %對(duì)復(fù)正弦序列加矩形窗（元素群乘）14程序要點(diǎn)n=n0:n1+N+5; %生成自變量數(shù)組15程序顯示結(jié)果15程序顯示結(jié)果16例7.3序列的移位和周期延拓已知 ，利用MATLAB生成并圖示 表示x(n)以8為周期的延拓）和 解：方法1，利用矩陣乘法和冒號(hào)運(yùn)算 x=[1234];y=x'*ones(1，3);方法2，采用求余函數(shù)mod，y=x(mod(n,M)+1)可實(shí)現(xiàn)對(duì)x(n)以M為周期的周期延拓。加1是因?yàn)镸ATLAB向量下標(biāo)只能從1開(kāi)始，16例7.3序列的移位和周期延拓已知 ，利用M17程序?qū)崿F(xiàn)分析構(gòu)造x(n)x1=(0.8).^n;x2=[(n>=0)&(n<M)];xn=x1.*x2;移位：fork=m+1:m+M

xm(k)=xn(k-m);

end17程序?qū)崿F(xiàn)分析構(gòu)造x(n)18周期化實(shí)現(xiàn)方法1：xc=xn(mod(n,8)+1); %產(chǎn)生x(n)的周期延拓

xcm=xn(mod(n-m,8)+1); %產(chǎn)生x(n)移位后的周期延拓18周期化實(shí)現(xiàn)方法1：19周期化實(shí)現(xiàn)方法2x=[1234];y=x'*ones(1,3);y=111222333444y1=y(:);%12341234123419周期化實(shí)現(xiàn)方法220程序結(jié)果20程序結(jié)果21例7.6離散序列的卷積計(jì)算給出兩個(gè)序列 和 ,計(jì)算其卷積y(n)，并圖示各輸入輸出序列。解：在例6.4中，已經(jīng)給出了直接調(diào)用MATLAB的卷積函數(shù)conv的方法，也給出了自編卷積計(jì)算程序的方法，要注意的是本例時(shí)間變量的設(shè)定和移位方法。在本例中，設(shè)定n為從零開(kāi)始，向量x和h的長(zhǎng)度分別為Nx=20和Nh=10；結(jié)果向量y的長(zhǎng)度為length(y)=Nx+Nh-1。21例7.6離散序列的卷積計(jì)算給出兩個(gè)序列 和 22程序要點(diǎn)x1=(0.9).^n;%x1(n)%產(chǎn)生x2(n)=x1(n-m)x2=zeros(1,Nx+m);fork=m+1:m+Nx

x2(k)=x1(k-m);endnh=0:Nh-1;h1=ones(1,Nh);%h1(n)h2=h1;%h2(n)y1=conv(x1,h1);%y1(n)y2=conv(x2,h2);%y2(n)Ny1=length(y1);%y1(n)長(zhǎng)度Ny2=length(y2);%y2(n)長(zhǎng)度22程序要點(diǎn)x1=(0.9).^n;%x1(n)23程序顯示結(jié)果23程序顯示結(jié)果24例7.4離散系統(tǒng)對(duì)信號(hào)的響應(yīng)本題給定6階低通數(shù)字濾波器的系統(tǒng)函數(shù)，求它在下列輸入序列x(n)下的輸出序列y(n)。解：本題的計(jì)算原理見(jiàn)例6.14，在這里用工具箱函數(shù)filter來(lái)解。如果已知系統(tǒng)函數(shù)H(z)=B(z)/A(z)，則filter函數(shù)可求出系統(tǒng)對(duì)輸入信號(hào)x(n)的響應(yīng)y(n)。y

=filter(B,A,x)由差分方程可得到H(z)的分子和分母多項(xiàng)式系數(shù)向量A和B,再給出輸入向量x即可。系統(tǒng)分析24例7.4離散系統(tǒng)對(duì)信號(hào)的響應(yīng)本題給定6階低通數(shù)字濾25程序要點(diǎn)%設(shè)定系統(tǒng)參數(shù)A,BB=0.0003738*conv([1,1],conv([1,1],conv([1,1],conv([1,1],conv([1,1],[1,1])))))A=conv([1,-1.2686,0.7051],conv([1,-1.0106,0.3583],[1,-0.9044,0.2155]))x1=[n==0]; %產(chǎn)生輸入信號(hào)x1(n)y1=filter(B,A,x1); %對(duì)x1(n)的響應(yīng)x2=[(n-m)==0]; %產(chǎn)生輸入信號(hào)x2(n)y2=filter(B,A,x2); %對(duì)x2(n)的響應(yīng)x3=[n>=0]; %產(chǎn)生輸入信號(hào)x3(n)y3=filter(B,A,x3); %對(duì)x3(n)的響應(yīng)x4=[(n>=0)&(n<32)]; %產(chǎn)生輸入信號(hào)x4(n)y4=filter(B,A,x4); %對(duì)x4(n)的響應(yīng)x5=exp(j*pi*n/8); %產(chǎn)生輸入信號(hào)x5(n)y5=filter(B,A,x5); %對(duì)x5(n)的響應(yīng)25程序要點(diǎn)%設(shè)定系統(tǒng)參數(shù)A,B26顯示結(jié)果26顯示結(jié)果27例7.5系統(tǒng)線性性質(zhì)驗(yàn)證設(shè)系統(tǒng)差分方程為y(n)=x(n)+0.8y(n-1)要求用程序驗(yàn)證系統(tǒng)的線性性質(zhì)。解：產(chǎn)生兩種輸入序列，分別乘以常數(shù)后1。分別激勵(lì)系統(tǒng)，再求輸出之和；2。先相加，再激勵(lì)系統(tǒng)求輸出；對(duì)兩個(gè)結(jié)果進(jìn)行比較，方法是求它們之差，按誤差的絕對(duì)值是否極小進(jìn)行判斷。27例7.5系統(tǒng)線性性質(zhì)驗(yàn)證設(shè)系統(tǒng)差分方程為28程序要點(diǎn)B=1;A=[1,-0.8];x1=0.8.^n; %產(chǎn)生輸入信號(hào)x1(n)x=[(n>=0)&(n<32)];x1=x1.*x;y1=filter(B,A,x1); %對(duì)x1(n)的響應(yīng)y1(n)x2=[(n-m)==0];y2=filter(B,A,x2); %對(duì)x2(n)的響應(yīng)y2(n)x3=5*x1+3*x2;y3=filter(B,A,x3); %對(duì)5x1(n)+3x2(n)的響應(yīng)y3(n)y=5*y1+3*y2; %y(n)=5y1(n)+3y2(n)28程序要點(diǎn)B=1;29程序顯示結(jié)果29程序顯示結(jié)果30二、離散時(shí)間信號(hào)與系統(tǒng)頻域分析30二、離散時(shí)間信號(hào)與系統(tǒng)頻域分析31包含內(nèi)容Z變換序列z變換與逆變換7.7Z反變換法求解系統(tǒng)響應(yīng)7.8離散時(shí)間傅里葉變換7.9離散傅里葉變換7.15系統(tǒng)頻域響應(yīng)7.12,7.1331包含內(nèi)容Z變換32

(1)

序列的正反Z變換

其中，符號(hào)表示取z變換，z是復(fù)變量。相應(yīng)地，單邊z變換定義為：Z變換與其逆變換32(1)

序列的正反Z變換

其中，符號(hào)表示取z33MATLAB符號(hào)數(shù)學(xué)工具箱提供了計(jì)算離散時(shí)間信號(hào)單邊z變換的函數(shù)ztrans和z反變換函數(shù)iztrans，其語(yǔ)句格式分別為Z=ztrans(x)x=iztrans(z)上式中的x和Z分別為時(shí)域表達(dá)式和z域表達(dá)式的符號(hào)表示，可通過(guò)sym函數(shù)來(lái)定義。MATLAB實(shí)現(xiàn)-方法1（符號(hào)數(shù)學(xué)工具箱）33MATLAB符號(hào)數(shù)學(xué)工具箱提供了計(jì)算離散時(shí)間信號(hào)單邊z變34

【例1】試用ztrans函數(shù)求下列函數(shù)的z變換。

>>x=sym('a^n*cos(pi*n)');>>Z=ztrans(x);>>simplify(Z)ans=z/(z+a)34【例1】試用ztrans函數(shù)求下列函數(shù)的z變換。>35

【例2】試用iztrans函數(shù)求下列函數(shù)的z反變換。

>>Z=sym('(8*z-19)/(z^2-5*z+6)');>>x=iztrans(Z);>>simplify(x)ans=-19/6*charfcn[0](n)+5*3^(n-1)+3*2^(n-1)-19/6*charfcn[0](n)+5*3^(n-1)+3*2^(n-1)charfcn[0](n)是函數(shù)在MATLAB符號(hào)工具箱中的表示，反變換后的函數(shù)形式為：35【例2】試用iztrans函數(shù)求下列函數(shù)的z反變換。36例7.14用符號(hào)運(yùn)算工具箱解z變換解：無(wú)限長(zhǎng)度時(shí)間序列的z變換和逆z變換都屬于符號(hào)運(yùn)算的范圍。MATLAB的symbolic（符號(hào)運(yùn)算）工具箱已提供了這種函數(shù)。如果讀者已在計(jì)算機(jī)上安裝了這個(gè)工具箱，可以鍵入以下程序。MATLAB程序q714.m其特點(diǎn)是程序的開(kāi)始要指定符號(hào)自變量symsznaNw0 %規(guī)定z,n,a…為符號(hào)變量36例7.14用符號(hào)運(yùn)算工具箱解z變換解：無(wú)限長(zhǎng)度時(shí)間序37程序分析：symszn a Nw0%規(guī)定z,n,a為符號(hào)變量y1=a^n; %給出y的表示式Y(jié)1=ztrans(y1) %求y的z變換X1=-3*z^-1/(2-5*z^-1+2*z^-2); %給出X的表示式x1=iztrans(X1) %求X的逆Z變換simplify(x);%可顯示變換后表達(dá)式的形式37程序分析：symszn a Nw0%規(guī)定z38如果信號(hào)的z域表示式是有理函數(shù)，進(jìn)行z反變換的另一個(gè)方法是對(duì)進(jìn)行部分分式展開(kāi)，然后求各簡(jiǎn)單分式的z反變換.如果X(z)的有理分式表示為：

MATLAB實(shí)現(xiàn)-方法2（residuez函數(shù)實(shí)現(xiàn)反變換）38如果信號(hào)的z域表示式是有理函數(shù)，進(jìn)行z反變換的另一個(gè)方法39

MATLAB信號(hào)處理工具箱提供了一個(gè)對(duì)進(jìn)行部分分式展開(kāi)的函數(shù)residuez，其語(yǔ)句格式為：[R,P,K]=residuez(B,A)其中，B，A分別表示X(z)的分子與分母多項(xiàng)式的系數(shù)向量，分子與分母多項(xiàng)式按照z^-1升冪排列。R為部分分式的系數(shù)向量；P為極點(diǎn)向量；K為多項(xiàng)式的系數(shù)。若X(z)為有理真分式，則K為零。39MATLAB信號(hào)處理工具箱提供了一個(gè)對(duì)進(jìn)行部分分式展40【例1】試用MATLAB命令進(jìn)行部分分式展開(kāi)，并求出其z反變換。解：MATLAB源程序?yàn)?gt;>B=[18];>>A=[18,3,-4,-1];>>[R,P,K]=residuez(B,A)R=0.36000.24000.400040【例1】試用MATLAB命令進(jìn)行部分分式展開(kāi)，并求出其z41

P=0.5000-0.3333-0.3333K=[]從運(yùn)行結(jié)果可知，表示系統(tǒng)有一個(gè)二重極點(diǎn)。所以，X(z)的部分分式展開(kāi)為41P=從運(yùn)行結(jié)果可知，表示系統(tǒng)有一個(gè)二重極點(diǎn)。42教材中求z的逆變換的方法對(duì)于z變換分式可以用部分分式法或長(zhǎng)除法求其反變換。用函數(shù)residuez可以求出它的極點(diǎn)留數(shù)分解其中 [r,p,k]=residuez(B,A)其反變換為：42教材中求z的逆變換的方法對(duì)于z變換分式43注意幾點(diǎn)：[r,p,k]=residuez(B,A),A首項(xiàng)不能為零，即分母z最高次項(xiàng)不能為0k只有當(dāng)分式為假分式時(shí)會(huì)有值，即對(duì)應(yīng)沖激序列部分（FIR）；真分式，值為零還有一些細(xì)節(jié)需要處理：43注意幾點(diǎn)：[r,p,k]=residuez(B,A),A44例7.7有限序列的z和逆z變換兩序列x1=[1,2,3]，n1=[-1:1]及x2=[2,4,3,5]，n2=[-2:1]，求出x1與x2及其卷積x的z變換。解：其z變換可寫(xiě)成兩個(gè)多項(xiàng)式乘積 可用conv函數(shù)來(lái)求得。n數(shù)組要自己判別。n的起點(diǎn)ns=ns1+ns2=3，終點(diǎn)nf=nf1+nf2=2。 n=ns:nf。由x和n即可得出X(z)。44例7.7有限序列的z和逆z變換兩序列x1=[1,45程序分析x1=[1,2,3];ns1=-1;

%設(shè)定x1和ns1nf1=ns1+length(x1)-1; %nf1可以算出x2=[2,4,3,5];ns2=-2; %設(shè)定x2和ns2nf2=ns1+length(x1)-1; %nf2可以算出x=conv(x1,x2) %求出xn=(ns1+ns2):(nf1+nf2) %求出n45程序分析x1=[1,2,3];ns1=-1; %464647例7.8求z多項(xiàng)式分式的逆變換設(shè)系統(tǒng)函數(shù)為 ，輸入例7.7中的x2信號(hào)，用z變換計(jì)算輸出y(n)解：由例7.7可知 ，故 Y(z)=X(z)W(z)=其中nsy=分母分子中z的最高冪次之差。調(diào)用[r,p,k]=residuez(B,A)，可由B，A求出r，p，k，進(jìn)而求逆z變換，得47例7.8求z多項(xiàng)式分式的逆變換設(shè)系統(tǒng)函數(shù)為 48例7.8z多項(xiàng)式分式逆變換(續(xù))由程序算出nsy=-1，留數(shù)、極點(diǎn)分別為r=-57.7581和204.7581p=0.7791和0.3209k=-150-30代入得48例7.8z多項(xiàng)式分式逆變換(續(xù))由程序算出nsy49Z反變換法求解系統(tǒng)輸出(編程實(shí)現(xiàn))Y(z)=X(z)W(z)x=[2,4,3,5];nsx=-2;%輸入序列及初始時(shí)間nfx=nsx+length(x)-1;%計(jì)算序列終止時(shí)間Bw=-3;nsbw=-1;%系統(tǒng)函數(shù)的分子系數(shù),及z的最高次數(shù)Aw=[2,-2.2,0.5];nsaw=0; %系統(tǒng)函數(shù)的分母系數(shù),及z的最高次數(shù)B=conv(-3,x); %輸入與分子z變換的多項(xiàng)式乘積A=Aw;

%分母不變y(n)=IZT[Y(Z)][r,p,k]=residuez(B,A) %求留數(shù)r,極點(diǎn)p及直接項(xiàng)k nf=input('終點(diǎn)時(shí)間nf='); %要求鍵入終點(diǎn)時(shí)間n=min(nsx,nsy):nf; %生成總時(shí)間數(shù)組%求無(wú)限序列yi和直接序列ydyi=(r(1)*p(1).^(n-nsy)+r(2)*p(2).^(n-nsy)).*stepseq(nsy,n(1),nf);yd=k(1)*impseq(nsy,n(1),nf)+k(2)*impseq(-1-nsy,n(1),nf);y=yi+yd;49Z反變換法求解系統(tǒng)輸出(編程實(shí)現(xiàn))Y(z)=X(z)W(50另外一種方法(filter函數(shù))：filter(B,A,x)50另外一種方法(filter函數(shù))：filter(B,A,51例7.9離散時(shí)間傅里葉變換取周期的正弦信號(hào)，作8點(diǎn)采樣，求它的連續(xù)頻譜。然后對(duì)該信號(hào)進(jìn)行N個(gè)周期延拓，再求它的連續(xù)頻譜。把N無(wú)限增大，比較分析其結(jié)果。解：先求離散傅立葉變換的MATLAB子程序最后得到X=x*exp(-j*w*n‘)。有了子程序，本例就沒(méi)有什么難度了。51例7.9離散時(shí)間傅里葉變換取周期的正弦信號(hào)，作8點(diǎn)52DTFT程序：function[X]=dtft(x,w)%計(jì)算離散時(shí)間傅立葉變換%[X]=dtft(x,n,w),%X=在w頻率點(diǎn)上的DTFT數(shù)組,%x=沿n的有限長(zhǎng)度序列,%n=樣本位置向量%w=頻率點(diǎn)位置向量n=1:length(x);ewn=exp(-n'*w*i);X=x*ewn;52DTFT程序：function[X]=dtft(x,w53正弦信號(hào)x0的DTFT變換x0=sin(2*pi*[1:8]/8)*5; %x0是8點(diǎn)行向量dt=2*pi/8;w=linspace(-2*pi,2*pi,1000)/dt; %w是1000點(diǎn)行向量X0=dtft(x0,w)*dt; %求得頻率響應(yīng)X053正弦信號(hào)x0的DTFT變換x0=sin(2*pi*[154正弦信號(hào)延拓N個(gè)周期后，DTFTN=4; %延拓周期數(shù)x1=reshape(x0'*ones(1,N),1,N*length(x0)); %延拓后的時(shí)域信號(hào)x1X1=dtft(x1,w)*dt; %求x1的頻率響應(yīng)X1%延拓100次后%x1=reshape(x0'*ones(1,100),1,100*length(x0)); %延拓后的時(shí)域信號(hào)x1%X1=dtft(x1,w)/100*dt; %求x1的頻率響應(yīng)X154正弦信號(hào)延拓N個(gè)周期后，DTFTN=4; %延拓周55重復(fù)無(wú)窮次disp('重復(fù)無(wú)窮次的八點(diǎn)信號(hào)的離散時(shí)間傅立葉變換-傅立葉級(jí)數(shù)')pause,X2=fft(x0*dt); %離散傅立葉變換w1=2*pi*[0:length(x0)-1]/length(x0); %離散頻點(diǎn)向量subplot(3,1,3),stem([-w1,w1],[abs(X2),abs(X2)]),grid,axis([min(w),max(w),0,max(abs(X2))]),grid55重復(fù)無(wú)窮次disp('重復(fù)無(wú)窮次的八點(diǎn)信號(hào)的離散時(shí)間傅立56例7.9離散時(shí)間傅里葉變換2程序運(yùn)行結(jié)果執(zhí)行程序q709并按提示鍵入N=4，所得圖形如圖7.10所示。N取得愈大，其峰值愈大，寬度愈窄。當(dāng)N取得很大時(shí)，會(huì)出現(xiàn)內(nèi)存不足的問(wèn)題，這是用矩陣乘法做傅里葉變換的缺點(diǎn)。另外，因?yàn)槟菚r(shí)峰值點(diǎn)處的寬度很窄，也會(huì)出現(xiàn)所選頻點(diǎn)對(duì)不上峰值點(diǎn)的問(wèn)題。所以對(duì)于N無(wú)限增大的情況，必須用fft函數(shù)來(lái)求。這時(shí)用連續(xù)頻譜也沒(méi)有意義了。這里用同樣的橫坐標(biāo)把幾種頻譜進(jìn)行對(duì)比，使讀者更好地理解其關(guān)系。56例7.9離散時(shí)間傅里葉變換2程序運(yùn)行結(jié)果57程序顯示結(jié)果：57程序顯示結(jié)果：58例7.10時(shí)域采樣頻率與頻譜混疊分別以采樣頻率fs=1000Hz，400Hz和200Hz對(duì)xa(t)進(jìn)行等間隔采樣，計(jì)算并圖示三種采樣頻率下的采樣信號(hào)及其幅頻特性解：程序分別設(shè)定4種采樣頻率fs=10kHz，1kHz，400Hz和200Hz，對(duì)xa(t)進(jìn)行采樣，得到采樣序列xa(t)，xa1(n)，xa2(n)，xa3(n)，畫(huà)出其幅度頻譜。采樣時(shí)間區(qū)間均為0.1秒。為了便于比較，畫(huà)出了幅度歸一化的幅頻曲線，如圖7.11所示。58例7.10時(shí)域采樣頻率與頻譜混疊分別以采樣頻率fs=59例7.10采樣頻率與頻譜混疊(續(xù))由于由以上關(guān)系式可見(jiàn)，采樣信號(hào)的頻譜函數(shù)是原模擬信號(hào)頻譜函數(shù)的周期延拓，延拓周期為2/T。如果以頻率f為自變量（

=2f），則以采樣頻率fs=1/T為延拓周期。對(duì)頻帶限于fc的模擬信號(hào)xa(t)，只有當(dāng)fs≥2fc時(shí)，采樣后 才不會(huì)發(fā)生頻譜混疊失真。這就是著名的采樣定理59例7.10采樣頻率與頻譜混疊(續(xù))由于程序要點(diǎn)：t=0:1/fs:0.1;%fs代不同的取樣頻率

xa=exp(-a*t).*sin(b*t);k=0:511;f=fs*k/512; %由wk=2πk/512=2πfT求得模擬頻率fXa=dtft(xa,2*pi*k/512);%近似模擬信號(hào)頻譜60程序要點(diǎn)：t=0:1/fs:0.1;%fs代不同的取樣結(jié)果161結(jié)果161結(jié)果262結(jié)果26263例7.12梳狀濾波器零極點(diǎn)和幅特性梳狀濾波器系統(tǒng)函數(shù)有如下兩種類型。FIR型：IIR型：freqz數(shù)字濾波器頻率特性計(jì)算和繪制函數(shù)zplaneH(z)的零-極點(diǎn)圖繪制。解：調(diào)用函數(shù)freqz和zplane很容易寫(xiě)出程序q712.m。63例7.12梳狀濾波器零極點(diǎn)和幅特性梳狀濾波器系統(tǒng)函數(shù)64程序分析freqzfreqz(b,a,w)zplane64程序分析freqz65頻率響應(yīng)b=[1,0,0,0,0,0,0,0,-1];a0=1;a1=[1,0,0,0,0,0,0,0,-(0.8)^8];a2=[1,0,0,0,0,0,0,0,-(0.9)^8];a3=[1,0,0,0,0,0,0,0,-(0.98)^8];[H,w]=freqz(b,a0);[H1,w1]=freqz(b,a1);[H2,w2]=freqz(b,a2);[H3,w3]=freqz(b,a3);65頻率響應(yīng)b=[1,0,0,0,0,0,0,0,-1];66零極點(diǎn)及繪圖subplot(2,2,1);zplane(b,a0);title('FIR梳狀濾波器零點(diǎn)圖');subplot(2,2,2);zplane(b,a1);title('IIR梳狀濾波器零、極點(diǎn)圖,a=0.8');subplot(2,2,3);plot(w/pi,abs(H));title('FIR梳狀濾波器幅頻響應(yīng)曲線');ylabel('幅度');xlabel('ω/π');subplot(2,2,4);plot(w1/pi,abs(H1));title('IIR梳狀濾波器幅頻響應(yīng)曲線,a=0.8');ylabel('幅度');xlabel('ω/π');figure(2);subplot(2,2,1);zplane(b,a2);title('IIR梳狀濾波器零、極點(diǎn)圖,a=0.9');subplot(2,2,2);zplane(b,a3);title('IIR梳狀濾波器零、極點(diǎn)圖,a=0.98)');subplot(2,2,3);plot(w2/pi,abs(H2));title('IIR梳狀濾波器幅頻響應(yīng)曲線,a=0.9');ylabel('幅度');xlabel('ω/π');subplot(2,2,4);plot(w3/pi,abs(H3));title('IIR梳狀濾波器幅頻響應(yīng)曲線,a=0.98');ylabel('幅度');xlabel('ω/π')66零極點(diǎn)及繪圖subplot(2,2,1);zplane(67結(jié)果顯示67結(jié)果顯示68結(jié)果顯示68結(jié)果顯示69例7.13低通濾波及時(shí)域卷積定理輸入信號(hào)x(n)=cos(0.04n)+cos(0.08n)+cos(0.4n)+0.3(n)，0≤n≤63通過(guò)低通濾波器，計(jì)算濾波器對(duì)x(n)的響應(yīng)輸出y(n)，并圖示x(n)和y(n)，觀察濾波效果。解：如前所述，只要求出H(z)=B(z)/A(z)的分子和分母多項(xiàng)式系數(shù)向量B和A，則可調(diào)用濾波器直接Ⅱ型實(shí)現(xiàn)函數(shù)filter對(duì)輸入信號(hào)x(n)進(jìn)行濾波。

y=filter(B,A,x)69例7.13低通濾波及時(shí)域卷積定理輸入信號(hào)x(n)70程序要點(diǎn)輸入信號(hào)構(gòu)造%產(chǎn)生輸入信號(hào)x(n)n=0:255;N=4096;x=cos(0.04*pi*n)+cos(0.08*pi*n)+cos(0.4*pi*n);w=randn(size(x));%產(chǎn)生正態(tài)零均值噪聲x=x+0.3*w;70程序要點(diǎn)輸入信號(hào)構(gòu)造71程序要點(diǎn)%求H(z)分子分母多項(xiàng)式系數(shù)向量B和Ab=[1,2,1];%(1+z-1)2的展開(kāi)系數(shù)B=0.0003738*conv(conv(b,b),b);%嵌套調(diào)用卷積函數(shù)conva1=[1,-1.2686,0.7051];a2=[1,-1.0106,0.3583];a3=[1,-0.9044,0.2155];A=conv(conv(a1,a2),a3);71程序要點(diǎn)%求H(z)分子分母多項(xiàng)式系數(shù)向量B和A72程序要點(diǎn)濾波及圖形%對(duì)x(n)濾波y=filter(B,A,x);%繪圖X=fft(x,N);Y=fft(y,N);subplot(3,2,1);stem(x,'.')axis([0,max(n)/4,min(x),max(x)]);line([0,max(n)],[0,0])title('輸入信號(hào)x(n)');xlabel('n');ylabel('x(n)')subplot(3,2,5);stem(y,'.')axis([0,max(n)/4,min(y),max(y)]);line([0,max(n)],[0,0])title('輸出信號(hào)y(n)');xlabel('n');ylabel('y(n)')k=0:N-1;f=2*k/N;subplot(3,2,2);plot(f,abs(X))title('輸入信號(hào)x(n)的幅頻曲線');xlabel('ω/π');ylabel('|FT[x(n)]|')axis([0,0.5,0,max(abs(X))]);subplot(3,2,6);plot(f,abs(Y))title('輸出信號(hào)y(n)的幅頻曲線');xlabel('ω/π');ylabel('|FT[y(n)]|')axis([0,0.5,0,max(abs(Y))]);72程序要點(diǎn)濾波及圖形73結(jié)果顯示：73結(jié)果顯示：74程序要點(diǎn)[h,f]=freqz(B,A,N,'whole');figure(2)subplot(3,2,1);plot(f/pi,abs(h))title('濾波器幅頻響應(yīng)曲線');xlabel('ω/π');ylabel('H幅度')axis([0,0.5,0,max(abs(h))]);Ym=h'.*X;subplot(3,2,2);plot(f/pi,abs(Ym))title('|FT[x(n)]FT[h(n)]|');xlabel('ω/π');ylabel('Ym幅度')axis([0,0.5,0,max(abs(Ym))]);74程序要點(diǎn)[h,f]=freqz(B,A,N,'whole75顯示結(jié)果75顯示結(jié)果76此處講一下-離散系統(tǒng)時(shí)頻域綜合分析76此處講一下-離散系統(tǒng)時(shí)頻域綜合分析771、系統(tǒng)函數(shù)的零極點(diǎn)分析離散時(shí)間系統(tǒng)的系統(tǒng)函數(shù)定義為系統(tǒng)零狀態(tài)響應(yīng)的z變換與激勵(lì)的z變換之比:

如果系統(tǒng)函數(shù)的有理函數(shù)表示式為離散系統(tǒng)時(shí)頻域綜合分析771、系統(tǒng)函數(shù)的零極點(diǎn)分析離散時(shí)間系統(tǒng)的系統(tǒng)函數(shù)定義為系統(tǒng)78在MATLAB中系統(tǒng)函數(shù)的零極點(diǎn)就可通過(guò)函數(shù)roots得到，也可借助DSP工具箱中的函數(shù)tf2zp得到，tf2zp的語(yǔ)句格式為：[R,P,K]=tf2zp(B,A)其中，B與A分別表示分子與分母多項(xiàng)式的系數(shù)向量。它的作用是將H(z)的有理分式表示式轉(zhuǎn)換為零極點(diǎn)增益形式：MATLAB實(shí)現(xiàn)離散系統(tǒng)時(shí)頻域綜合分析78在MATLAB中系統(tǒng)函數(shù)的零極點(diǎn)就可通過(guò)函數(shù)roots得79離散系統(tǒng)時(shí)頻域綜合分析【例1】已知一離散因果LTI系統(tǒng)的系統(tǒng)函數(shù)為：試用MATLAB命令求該系統(tǒng)的零極點(diǎn)。

。79離散系統(tǒng)時(shí)頻域綜合分析【例1】已知一離散因果LTI系統(tǒng)80>>B=[1,0.32];>>A=[1,1,0.16];>>[R,P,K]=tf2zp(B,A)R=-0.3200P=-0.8000-0.2000K=1因此，零點(diǎn)為，極點(diǎn)為：。求解：80>>B=[1,0.32];，極點(diǎn)為：。求解：81離散系統(tǒng)時(shí)頻域綜合分析若要獲得系統(tǒng)函數(shù)的零極點(diǎn)分布圖，可直接應(yīng)用zplane函數(shù)，其語(yǔ)句格式為：zplane(B,A)其中，B與A分別表示的分子和分母多項(xiàng)式的系數(shù)向量。它的作用是在Z平面上畫(huà)出單位圓、零點(diǎn)與極點(diǎn)。。81離散系統(tǒng)時(shí)頻域綜合分析若要獲得系統(tǒng)函數(shù)的零極點(diǎn)分布圖，可82離散系統(tǒng)時(shí)頻域綜合分析【例2】已知一離散因果LTI系統(tǒng)的系統(tǒng)函數(shù)為：試用MATLAB命令繪出該系統(tǒng)的零極點(diǎn)分布圖。。82離散系統(tǒng)時(shí)頻域綜合分析【例2】已知一離散因果LTI系統(tǒng)83求解：。MATLAB源程序?yàn)椋?gt;>B=[1,0,-0.36];>>A=[1,-1.52,0.68];>>zplane(B,A),gridon>>legend('零點(diǎn)','極點(diǎn)')>>title('零極點(diǎn)分布圖')程序運(yùn)行結(jié)果如圖所示。83求解：。MATLAB源程序?yàn)椋?4離散系統(tǒng)時(shí)頻域綜合分析系統(tǒng)函數(shù)的零極點(diǎn)分布與其時(shí)域特性的關(guān)系：在離散系統(tǒng)中，z變換建立了時(shí)域函數(shù)與z域函數(shù)之間的對(duì)應(yīng)關(guān)系。因此，z變換的函數(shù)從形式可以反映的部分內(nèi)在性質(zhì)。我們通過(guò)討論H(z)的一階極點(diǎn)情況，來(lái)說(shuō)明系統(tǒng)函數(shù)的零極點(diǎn)分布與系統(tǒng)時(shí)域特性的關(guān)系。84離散系統(tǒng)時(shí)頻域綜合分析系統(tǒng)函數(shù)的零極點(diǎn)分布與其時(shí)域特性的85離散系統(tǒng)時(shí)頻域綜合分析MATLAB求解單位抽樣響應(yīng)可利用函數(shù)filter，

filter函數(shù)的常用語(yǔ)句格式為：y=filter(b,a,x)表示由向量b和a組成的系統(tǒng)對(duì)輸入x進(jìn)行濾波，系統(tǒng)的輸出為y;

2、系統(tǒng)時(shí)域響應(yīng)分析85離散系統(tǒng)時(shí)頻域綜合分析MATLAB求解單位抽樣響應(yīng)86離散系統(tǒng)時(shí)頻域綜合分析MATLAB另一種求單位抽樣響應(yīng)的方法是利用控制系統(tǒng)工具箱提供的函數(shù)impz來(lái)實(shí)現(xiàn)。impz函數(shù)的常用語(yǔ)句格式為impz(b,a,N)其中，參數(shù)N通常為正整數(shù)，代表計(jì)算單位抽樣響應(yīng)的樣值個(gè)數(shù)。86離散系統(tǒng)時(shí)頻域綜合分析MATLAB另一種求單位抽樣響應(yīng)87離散系統(tǒng)時(shí)頻域綜合分析【例1】試用MATLAB命令畫(huà)出系統(tǒng)函數(shù)的零極點(diǎn)分布圖、以及對(duì)應(yīng)的時(shí)域單位抽樣響應(yīng)的波形。

b1=[1,0];a1=[1,-0.8];subplot(121)zplane(b1,a1)title('極點(diǎn)在單位圓內(nèi)的正實(shí)數(shù)')subplot(122)impz(b1,a1,30);gridon;87離散系統(tǒng)時(shí)頻域綜合分析【例1】試用MATLAB命令畫(huà)出88結(jié)果顯示：88結(jié)果顯示：89離散系統(tǒng)時(shí)頻域綜合分析3、離散時(shí)間LTI系統(tǒng)的頻率特性分析

離散時(shí)間系統(tǒng)的頻率響應(yīng)定義為：其中，稱為離散時(shí)間系統(tǒng)的幅頻特性；稱為離散時(shí)間系統(tǒng)的相頻特性。

89離散系統(tǒng)時(shí)頻域綜合分析3、離散時(shí)間LTI系統(tǒng)的頻率特性分90離散系統(tǒng)時(shí)頻域綜合分析MATLAB提供了求離散時(shí)間系統(tǒng)頻響特性的函數(shù)freqz，調(diào)用freqz的格式主要有兩種。一種形式為[H,w]=freqz(B,A,N)其中，B與A分別表示的分子和分母多項(xiàng)式的系數(shù)向量；N為正整數(shù)，默認(rèn)值為512；

返回值w包含范圍內(nèi)的N個(gè)頻率等分點(diǎn)；返回值H則是離散時(shí)間系統(tǒng)頻率響應(yīng)。90離散系統(tǒng)時(shí)頻域綜合分析MATLAB提供了求離散時(shí)間系統(tǒng)頻91離散系統(tǒng)時(shí)頻域綜合分析另一種形式為：[H,w]=freqz(B,A,N,’whole’)與第一種方式不同之處在于角頻率的范圍擴(kuò)展到91離散系統(tǒng)時(shí)頻域綜合分析另一種形式為：92離散系統(tǒng)時(shí)頻域綜合分析【例1】試用MATLAB命令繪制以下系統(tǒng)的頻率響應(yīng)曲線。解：利用函數(shù)freqz計(jì)算出然后利用函數(shù)abs和angle分別求出幅頻特性與相頻特性，最后利用plot命令繪出曲線。92離散系統(tǒng)時(shí)頻域綜合分析【例1】試用MATLAB命令繪制93離散系統(tǒng)時(shí)頻域綜合分析MATLAB源程序?yàn)椋篵=[1-0.960.9028];a=[1-1.560.8109];[H,w]=freqz(b,a,400,'whole');Hm=abs(H);Hp=angle(H);subplot(211)plot(w,Hm),gridonxlabel('\omega(rad/s)'),ylabel('Magnitude')title('離散系統(tǒng)幅頻特性曲線')subplot(212)plot(w,Hp),gridonxlabel('\omega(rad/s)'),ylabel('Phase')title('離散系統(tǒng)相頻特性曲線')93離散系統(tǒng)時(shí)頻域綜合分析MATLAB源程序?yàn)椋?4結(jié)果顯示：94結(jié)果顯示：95最后一個(gè)內(nèi)容DFT變換95最后一個(gè)內(nèi)容DFT變換967.3離散傅里葉變換（DFT）定義DFT：用類似于例7.9中的方法，可把（7.3）式寫(xiě)成矩陣乘法運(yùn)算。其中，xn為序列行向量，Wnk是一N×N階方陣，而 稱為旋轉(zhuǎn)因子。967.3離散傅里葉變換（DFT）定義DFT：977.3離散傅里葉變換（DFT）用矩陣乘法計(jì)算N點(diǎn)DFT的程序如下。

MATLAB程序q73a.m%用矩陣乘法計(jì)算N點(diǎn)DFTclear;closeallxn=input('請(qǐng)輸入序列x=');N=length(xn); %n=0:N-1;k=n;nk=n'*k; %生成N×N方陣WN=exp(-j*2*pi/N);Wnk=WN.^nk; %產(chǎn)生旋轉(zhuǎn)因子矩陣Xk=xn*Wnk; %計(jì)算N點(diǎn)DFT977.3離散傅里葉變換（DFT）用矩陣乘法計(jì)算N點(diǎn)DF98例7.15序列的離散傅立葉變換求復(fù)正弦序列

余弦序列

正弦序列的離散傅立葉變換，分別按N=16和N=8進(jìn)行計(jì)算。繪出幅頻特性曲線，進(jìn)行比較討論。98例7.15序列的離散傅立葉變換求復(fù)正弦序列99程序分析clear;closeallN=16;N1=8;%%產(chǎn)生序列x1(n),計(jì)算DFT[x1(n)]n=0:N-1;x1n=exp(j*pi*n/8);%產(chǎn)生x1(n)X1k=fft(x1n,N);%計(jì)算N點(diǎn)DFT[x1(n)]Xk1=fft(x1n,N1);%計(jì)算N1點(diǎn)DFT[x1(n)]%%產(chǎn)生序列x2(n),計(jì)算DFT[x2(n)]x2n=cos(pi*n/8);X2k=fft(x2n,N);%計(jì)算N點(diǎn)DFT[x2(n)]Xk2=fft(x2n,N1);%計(jì)算N1點(diǎn)DFT[x1(n)]%%產(chǎn)生序列x3(n),計(jì)算DFT[x3(n)]x3n=sin(pi*n/8);X3k=fft(x3n,N);%計(jì)算N點(diǎn)DFT[x3(n)]Xk3=fft(x3n,N1);%計(jì)算N1點(diǎn)DFT[x1(n)]%99程序分析clear;closeall100例7.15序列的離散傅立葉變換在截取16點(diǎn)時(shí)，得到的是完整的余弦波形；而截取8點(diǎn)時(shí)，得到的是半截的余弦波形，當(dāng)然有大量的諧波成分。100例7.15序列的離散傅立葉變換在截取16點(diǎn)時(shí)，得到101本題分析DFT求解結(jié)果會(huì)分析101本題分析DFT求解102例7.16驗(yàn)證N點(diǎn)DFT的物理意義（1） ，繪出幅頻曲線和相頻曲線。（2）計(jì)算并圖示x(n)的8點(diǎn)DFT。（3）計(jì)算并圖示x(n)的16點(diǎn)DFT。解:序列x(n)的Ｎ點(diǎn)DFT的物理意義是 在[0，2]上進(jìn)行Ｎ點(diǎn)等間隔采樣。程序先密集采樣，繪制出幅頻曲線圖。然后再分別做8點(diǎn)和16點(diǎn)DFT來(lái)驗(yàn)證這個(gè)采樣關(guān)系。102例7.16驗(yàn)證N點(diǎn)DFT的物理意義（1） 結(jié)果1：103結(jié)果1：103結(jié)果2104結(jié)果2104結(jié)果3105結(jié)果3105106對(duì)本章《MATLAB在數(shù)字信號(hào)處理中的應(yīng)用》的內(nèi)容講解，告一段落。106對(duì)本章《MATLAB在數(shù)字信號(hào)處理中的應(yīng)用》的內(nèi)容講解107例7.10時(shí)域采樣頻率與頻譜混疊分別以采樣頻率fs=1000Hz，400Hz和200Hz對(duì)xa(t)進(jìn)行等間隔采樣，計(jì)算并圖示三種采樣頻率下的采樣信號(hào)及其幅頻特性解：程序分別設(shè)定4種采樣頻率fs=10kHz，1kHz，400Hz和200Hz，對(duì)xa(t)進(jìn)行采樣，得到采樣序列xa(t)，xa1(n)，xa2(n)，xa3(n)，畫(huà)出其幅度頻譜。采樣時(shí)間區(qū)間均為0.1秒。為了便于比較，畫(huà)出了幅度歸一化的幅頻曲線，如圖7.11所示。107例7.10時(shí)域采樣頻率與頻譜混疊分別以采樣頻率fs108例7.10采樣頻率與頻譜混疊(續(xù))由于由以上關(guān)系式可見(jiàn)，采樣信號(hào)的頻譜函數(shù)是原模擬信號(hào)頻譜函數(shù)的周期延拓，延拓周期為2/T。如果以頻率f為自變量（

=2f），則以采樣頻率fs=1/T為延拓周期。對(duì)頻帶限于fc的模擬信號(hào)xa(t)，只有當(dāng)fs≥2fc時(shí)，采樣后 才不會(huì)發(fā)生頻譜混疊失真。這就是著名的采樣定理108例7.10采樣頻率與頻譜混疊(續(xù))由于109例7.11由離散序列恢復(fù)模擬信號(hào)用時(shí)域內(nèi)插公式

其中 模擬用理想低通濾波器恢復(fù)的過(guò)程，觀察恢復(fù)波形，計(jì)算出最大恢復(fù)誤差。解：這個(gè)公式與卷積公式相像，可以用向量和矩陣乘法來(lái)解決。109例7.11由離散序列恢復(fù)模擬信號(hào)用時(shí)域內(nèi)插公式110例7.11由離散序列恢復(fù)模擬信號(hào)xa=x*g(TNM)=x*G其中G=sinc(Fs*TNM)M表示在兩個(gè)采樣點(diǎn)之間增加的間隔數(shù)，使輸出更密，更接近模擬信號(hào)。110例7.11由離散序列恢復(fù)模擬信號(hào)xa=x*111例7.16驗(yàn)證N點(diǎn)DFT的物理意義（1） ，繪出幅頻曲線和相頻曲線。（2）計(jì)算并圖示x(n)的8點(diǎn)DFT。（3）計(jì)算并圖示x(n)的16點(diǎn)DFT。解:序列x(n)的Ｎ點(diǎn)DFT的物理意義是 在[0，2]上進(jìn)行Ｎ點(diǎn)等間隔采樣。程序先密集采樣，繪制出幅頻曲線圖。然后再分別做8點(diǎn)和16點(diǎn)DFT來(lái)驗(yàn)證這個(gè)采樣關(guān)系。111例7.16驗(yàn)證N點(diǎn)DFT的物理意義（1） 112例7.17頻域與時(shí)域采樣對(duì)偶性（1）產(chǎn)生三角波序列（2）對(duì)M=40，計(jì)算x(n)的64點(diǎn)DFT，并圖示x(n)和X(k)=DFT[x(n)]，k=0,1,…,63。（3）對(duì)（2）中所得X(k)在[0，2]上進(jìn)行32點(diǎn)抽樣得（4）求 的32點(diǎn)IDFT，即

（5）繪出 的波形圖，評(píng)述它與x(n)的關(guān)系。112例7.17頻域與時(shí)域采樣對(duì)偶性（1）產(chǎn)生三角波序113例7.17頻域與時(shí)域采樣對(duì)偶性由于頻域在[0，2]上的采樣點(diǎn)數(shù)N（N=32）小于x(n)的長(zhǎng)度M（M=40），所以，產(chǎn)生時(shí)域混疊現(xiàn)象，不能由X1(k)恢復(fù)原序列x(n)。只有滿足N≥M時(shí)，可由頻域采樣X(jué)1(k)得到原序列x(n)。這就是頻域采樣定理。對(duì)N≥M的情況，請(qǐng)讀者自己編程驗(yàn)證。113例7.17頻域與時(shí)域采樣對(duì)偶性由于頻域在[0，2114例7.18快速卷積快速卷積就是根據(jù)DFT的循環(huán)卷積性質(zhì)，將時(shí)域卷積轉(zhuǎn)換為頻域相乘，最后再進(jìn)行IDFT得到時(shí)域卷積序列y(n)。其中時(shí)域和頻域之間的變換均用FFT實(shí)現(xiàn)，所以使卷積速度大大提高?？驁D如下:

114例7.18快速卷積快速卷積就是根據(jù)DFT的循環(huán)115例7.19用DFT求連續(xù)信號(hào)頻譜在計(jì)算機(jī)上用DFT對(duì)模擬信號(hào)進(jìn)行譜分析時(shí)，只能以有限大的采樣頻率fs對(duì)模擬信號(hào)采樣有限點(diǎn)樣本序列（等價(jià)于截取模擬信號(hào)一段進(jìn)行采樣）作DFT變換，得到模擬信號(hào)的近似頻譜。其誤差主要來(lái)自以下因素：①截?cái)嘈?yīng)（頻譜泄露和譜間干擾）②頻譜混疊失真因素①使譜分辨率（能分辨開(kāi)的兩根譜線間的最小間距）降低，并產(chǎn)生譜間干擾；因素②使折疊頻率（fs/2）附近的頻譜產(chǎn)生較大失真。115例7.19用DFT求連續(xù)信號(hào)頻譜在計(jì)算機(jī)上用DFT116例7.19用DFT求連續(xù)信號(hào)頻譜加大截取長(zhǎng)度Tp可提高頻率分辨率；選擇合適的窗函數(shù)可降低譜間干擾；而頻譜混疊失真要通過(guò)提高采樣頻率fs和（或）預(yù)濾波（在采樣之前濾除折疊頻率以外的頻率成分）來(lái)改善。編寫(xiě)程序q719.m驗(yàn)證截?cái)嘈?yīng)及加窗的改善作用，先選取以下參數(shù)：采樣頻率fs=400Hz，T=1/fs采樣信號(hào)序列對(duì)x(n)作4096點(diǎn)DFT作為的近似頻譜Xa(jf)。116例7.19用DFT求連續(xù)信號(hào)頻譜加大截取長(zhǎng)度Tp可117例7.19用DFT求連續(xù)信號(hào)頻譜如圖7.19所示。圖中X1(jf)，X4(jf)和X8(jf)分別表示Tp=0.04s，0.04*4s和0.04*8s時(shí)的譜分析結(jié)果。由圖可見(jiàn)，由于截?cái)嗍乖l譜中的單頻譜線展寬（又稱之為泄漏），Tp越大泄漏越小，頻率分辨率越高。Tp=0.04s時(shí)，25Hz與50Hz兩根譜線已分辨不清了。所以實(shí)際譜分析的截取時(shí)間Tp是由頻率分辨率決定的。另外，在本應(yīng)為零的頻段上出現(xiàn)了一些參差不齊的小譜包（稱為譜間干擾）。譜間干擾的大小取決于加窗的類型。用矩形窗比用Hamming窗的頻率分辨率高（泄漏?。V間干擾剛好相反。117例7.19用DFT求連續(xù)信號(hào)頻譜如圖7.19所示。118例7.20IIR濾波器直接型的轉(zhuǎn)換程序調(diào)用了信號(hào)處理工具箱函數(shù)tf2sos和擴(kuò)展函數(shù)dir2par，

[sos,g]=tf2sos(B,A)實(shí)現(xiàn)從直接型到級(jí)聯(lián)型（二階分割形式）的轉(zhuǎn)換。g為式中的增益，sos為L(zhǎng)×6階矩陣，表示式中的系數(shù)。118例7.20IIR濾波器直接型的轉(zhuǎn)換程序調(diào)用了信號(hào)處119例7.20IIR濾波器直接型的轉(zhuǎn)換[Cp,Bp,Ap]=dir2par(B,A)實(shí)現(xiàn)從直接型到并聯(lián)型的轉(zhuǎn)換。B為直接型H(z)的分子多項(xiàng)式系數(shù)向量，A為直接型H(z)的分母多項(xiàng)式系數(shù)向量；Cp，Bp，Ap的含義與擴(kuò)展函數(shù)dir2par中的C，B，A相同。dir2par中又調(diào)用了復(fù)共軛對(duì)比較函數(shù)cplxcomp。由于dir2par和cplxcomp是文獻(xiàn)[7]中開(kāi)發(fā)的，不屬于MATLAB工具箱函數(shù)，所以將其M文件清單附在程序q720.m之后。119例7.20IIR濾波器直接型的轉(zhuǎn)換[Cp,B120例7.20IIR濾波器直接型的轉(zhuǎn)換根據(jù)計(jì)算結(jié)果,級(jí)聯(lián)型H(z)表達(dá)式：

級(jí)聯(lián)型結(jié)構(gòu)圖。120例7.20IIR濾波器直接型的轉(zhuǎn)換根據(jù)計(jì)算結(jié)果,121例7.20IIR濾波器直接型的轉(zhuǎn)換并聯(lián)型結(jié)構(gòu)H(z)表達(dá)式并聯(lián)型結(jié)構(gòu)圖121例7.20IIR濾波器直接型的轉(zhuǎn)換并聯(lián)型結(jié)構(gòu)H(z122例7.21直接型結(jié)構(gòu)到格型梯形結(jié)構(gòu)tf2latc函數(shù)實(shí)現(xiàn)直接型到格型轉(zhuǎn)換

[K,C]=tf2latc(B,A)求出零-極點(diǎn)IIR系統(tǒng)格型梯形結(jié)構(gòu)的格型參數(shù)向量K和梯形參數(shù)向量C（用A(1)歸一化）。注意，當(dāng)系統(tǒng)函數(shù)在單位圓上有極點(diǎn)時(shí)發(fā)生錯(cuò)誤。

K=tf2latc(1,A)求出全極點(diǎn)IIR系統(tǒng)的格型結(jié)構(gòu)參數(shù)向量K。如果使用格式[K,C]=tf2latc(1,A)，返回的系數(shù)C為標(biāo)量。

K=tf2latc(B)求出FIR系統(tǒng)的格型梯形結(jié)構(gòu)參數(shù)（反射系數(shù)）向量K（用H(z)的常數(shù)項(xiàng)B(1)歸一化）。122例7.21直接型結(jié)構(gòu)到格型梯形結(jié)構(gòu)tf2lat123例7.21直接型結(jié)構(gòu)到格型梯形結(jié)構(gòu)直接型系統(tǒng)函數(shù)轉(zhuǎn)換為格型梯形結(jié)構(gòu)123例7.21直接型結(jié)構(gòu)到格型梯形結(jié)構(gòu)直接型系統(tǒng)函數(shù)124例7.22FIR濾波器直接型到其他型系統(tǒng)函數(shù)為調(diào)用信號(hào)處理工具箱函數(shù)tf2sos和tf2latc,給變?cè)狝賦值1,B=[2,13/12,5/4,2/3]即可.級(jí)聯(lián)型結(jié)構(gòu)系數(shù)sos=1.00000.536001.0000001.00000.00570.62191.000000g=2格型結(jié)構(gòu)系數(shù)（反射系數(shù)）：K=0.25000.50000.3333124例7.22FIR濾波器直接型到其他型系統(tǒng)函數(shù)為125例7.22FIR濾波器直接型到其他型得出級(jí)聯(lián)型為格型結(jié)構(gòu)為125例7.22FIR濾波器直接型到其他型得出級(jí)聯(lián)型為126例7.23FIR格型到直接型轉(zhuǎn)換給定K=[2,1/4,1/2,1/3],用函數(shù)latc2tf即可由B=latc2tf(K)得到B,寫(xiě)出直接型結(jié)構(gòu)126例7.23FIR格型到直接型轉(zhuǎn)換給定K=[2127例7.24系統(tǒng)函數(shù)的計(jì)算機(jī)推導(dǎo)數(shù)字濾波器的網(wǎng)絡(luò)結(jié)構(gòu)圖實(shí)際上也是一種信號(hào)流圖。因此【例6.20】中介紹的方法和公式同樣可以用來(lái)求離散域的數(shù)字濾波器的系統(tǒng)函數(shù)。不同的地方僅僅在于節(jié)點(diǎn)方程中出現(xiàn)了作為系數(shù)的符號(hào)變量z1，它將出現(xiàn)在系數(shù)矩陣中。MATLAB是不能處理上標(biāo)變量的，因此在程序中設(shè)q=z1，在計(jì)算完成后再人工地把結(jié)果中的q恢復(fù)為z1。127例7.24系統(tǒng)函數(shù)的計(jì)算機(jī)推導(dǎo)數(shù)字濾波器的網(wǎng)絡(luò)128例7.24的結(jié)構(gòu)圖與方程128例7.24的結(jié)構(gòu)圖與方程129例7.24的方程的矩陣形式由此可以求出系統(tǒng)函數(shù)129例7.24的方程的矩陣形式由此可以求出系統(tǒng)函數(shù)130例7.24解出的系統(tǒng)函數(shù)程序運(yùn)行的結(jié)果如果加入一個(gè)激勵(lì)x(n)A(n)，則得出130例7.24解出的系統(tǒng)函數(shù)程序運(yùn)行的結(jié)果1317.5FIR數(shù)字濾波器設(shè)計(jì)濾波器的特性指標(biāo)用絕對(duì)值δ1,δ2表示;用分貝Rp,Rs表示1317.5FIR數(shù)字濾波器設(shè)計(jì)濾波器的特性指標(biāo)132（1）窗函數(shù)法設(shè)計(jì)FIR濾波器先根據(jù)c和N求出相應(yīng)的理想濾波器單位脈沖響應(yīng)hd(n)。第二步要選擇合適的窗函數(shù)w(n)來(lái)截取hd(n)的適當(dāng)長(zhǎng)度（即階數(shù)），以保證實(shí)現(xiàn)要求的阻帶衰減；最后得到FIR濾波器單位脈沖響應(yīng)h(n)=hd(n).*w(n),即其系數(shù)。132（1）窗函數(shù)法設(shè)計(jì)FIR濾波器先根據(jù)c和N求出相應(yīng)的133（2）等波紋最佳一致逼近法（2）等波紋最佳一致逼近法:信號(hào)處理工具箱采用remez算法實(shí)現(xiàn)線性相位FIR數(shù)字濾波器的等波紋最佳一致逼近設(shè)計(jì)。其優(yōu)點(diǎn)是，設(shè)計(jì)指標(biāo)相同時(shí)，使濾波器階數(shù)最低；或階數(shù)相同時(shí)，使通帶最平坦，阻帶最小衰減最大；通帶和阻帶均為等波紋形式，最適合設(shè)計(jì)片段常數(shù)特性的濾波器。其調(diào)用格式如下：

b=remez(N,f,m,w,'ftype')其中N由remezord函數(shù)求出:[N,fo,mo,w]=remezord(f,m,dev,Fs)輸入變?cè)猟ev為各逼近頻段允許的波紋振幅。remez函數(shù)可直接調(diào)用remezord返回的參數(shù)如下：

b=remez(N,fo,mo,w)133（2）等波紋最佳一致逼近法（2）等波紋最佳一致逼近法:134例7.25窗函數(shù)法設(shè)計(jì)數(shù)字濾波器分別用矩形窗和Hamming窗設(shè)計(jì)線性相位FIR低通濾波器。要求通帶截止頻率c=/4，單位脈沖響應(yīng)h(n)的長(zhǎng)度N=21。繪出h(n)及其幅頻響應(yīng)特性曲線。先求出相應(yīng)的理想濾波器（本例應(yīng)為理想低通）單位脈沖響應(yīng)hd(n)，再根據(jù)阻帶最小衰減選擇合適的窗函數(shù)w(n)，最后得到FIR濾波器單位脈沖響應(yīng)h(n)=hd(n).*w(n)。134例7.25窗函數(shù)法設(shè)計(jì)數(shù)字濾波器分別用矩形窗和H135例7.25窗函數(shù)法設(shè)計(jì)數(shù)字濾波器本題中，c=/4，N=21，所以線性相位理想低通濾波器的單位脈沖響應(yīng)為:為了滿足線性相位FIR濾波器條件h(n)=h(N-1-n)，要求

=(N-1)/2=10。信號(hào)處理工具箱中有窗生成函數(shù)boxcar，hamming，hanning和blackman等。135例7.25窗函數(shù)法設(shè)計(jì)數(shù)字濾波器本題中，c=136例7.25窗函數(shù)法設(shè)計(jì)數(shù)字濾波器對(duì)兩種窗函數(shù)的設(shè)計(jì)結(jié)果分別如右圖7.25-1和圖7.25-2所示。136例7.25窗函數(shù)法設(shè)計(jì)數(shù)字濾波器對(duì)兩種窗函數(shù)的設(shè)137工具箱設(shè)計(jì)函數(shù)fir1和fir2MATLAB提供了基于窗函數(shù)法的FIR濾波器設(shè)計(jì)函數(shù)fir1和fir2,其功能及用法如下。fir1功能：標(biāo)準(zhǔn)頻率響應(yīng)形狀。格式：b=fir1(N,wc,‘ftype’,window)。當(dāng)wc=[wc1，wc2]時(shí)，是的帶通濾波器。當(dāng)ftype=high時(shí)，設(shè)計(jì)高通FIR濾波器；當(dāng)ftype=stop時(shí)，設(shè)計(jì)帶阻FIR濾波器。fir2功能：任意頻率響應(yīng)形狀。格式：b=fir2(N,f,m,window)137工具箱設(shè)計(jì)函數(shù)fir1和fir2MATLAB提供了基于138例7.26窗函數(shù)法設(shè)計(jì)帶通濾波器使用fir1函數(shù)b=fir1(N,wc,window)編程參數(shù)c為行向量c

=[lp/，hp/]根據(jù)阻帶最小衰減Rs

=60dB選擇窗函數(shù)類型和階次?？梢圆樯厦媪谐龅摹按昂瘮?shù)設(shè)計(jì)濾波器時(shí)的階數(shù)選擇表”。選blackman窗，其濾波器阻帶最小衰減可達(dá)到74dB，其窗口長(zhǎng)度M由過(guò)渡帶寬度B=0.15決定，Blackman窗設(shè)計(jì)的濾波器過(guò)渡帶寬度為12/M，故M取80。因M=N+1，所以濾波器階數(shù)N=79。138例7.26窗函數(shù)法設(shè)計(jì)帶通濾波器使用fir1函數(shù)139例7.27用remez函數(shù)低通濾波器解：先由題意計(jì)算設(shè)計(jì)參數(shù)f=[1/4,5/16]，m=[1,0];dev的計(jì)算稍復(fù)雜一些，由于所以有了這幾個(gè)參數(shù)就可以調(diào)用remezord和remez函數(shù)了.139例7.27用remez函數(shù)低通濾波器解：先由題意計(jì)140例7.27用remez函數(shù)低通濾波器橫線為-3dB，兩條豎線分別位于頻率/4和5/16。顯然，通帶指標(biāo)稍有富裕，過(guò)渡帶寬度和阻帶最小衰減剛好滿足指標(biāo)要求。程序輸出的幅頻特性140例7.27用remez函數(shù)低通濾波器橫線為-3dB141例7.28remez函數(shù)設(shè)計(jì)高通濾波器觀察等波紋逼近法中加權(quán)系數(shù)w()及濾波器階數(shù)N的作用和影響。期望逼近的濾波器通帶為[3/4，]，阻帶為[0，23/32]。解:在濾波器設(shè)計(jì)中，技術(shù)指標(biāo)越高，實(shí)現(xiàn)濾波器的階數(shù)也就越高。另外，對(duì)固定的階數(shù)，通帶與阻帶指標(biāo)可以互換，過(guò)渡帶寬度與通帶波紋和阻帶衰減指標(biāo)可以互換。取f=[0,3/4,23/32,1]，m=[0,0,1,1]。其余參數(shù)分三種情況進(jìn)行設(shè)計(jì)，①N=30，w=[1,1]；②N=30，w=[1,5];③N=60，w=[1,1]。141例7.28remez函數(shù)設(shè)計(jì)高通濾波器觀察等波紋逼142例7.28remez函數(shù)設(shè)計(jì)高通濾波器程序運(yùn)行結(jié)果如圖由圖可見(jiàn)，w較大的頻段逼近精度較高；w較小的頻段逼近精度較低。N較大時(shí)逼近精度較高,N較小時(shí)逼近精度較低。142例7.28remez函數(shù)設(shè)計(jì)高通濾波器程序運(yùn)行結(jié)1437.6IIR數(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)單。本節(jié)主要結(jié)合例題介紹這些IIR濾波器設(shè)計(jì)的工具箱函數(shù)。IIR數(shù)字濾波器的設(shè)計(jì)步驟由以下的流程圖來(lái)表示。下面以巴特沃斯濾波器設(shè)計(jì)函數(shù)為典型，介紹此流程圖中函數(shù)的功能和用法。1437.6IIR數(shù)字濾波器設(shè)計(jì)IIR數(shù)字濾波器設(shè)計(jì)的144IIR數(shù)字濾波器設(shè)計(jì)流程圖模擬低通濾波器原型設(shè)計(jì)Buttap,cheb1ap,cheb2apbesselap，ellipap函數(shù)頻率變換(變?yōu)楦咄?，帶通，帶阻?lp2lp，lp2hp，lp2bp，lp2bs模擬數(shù)字變換bilinearimpinvar合為一步的設(shè)計(jì)函數(shù)butter，cheb1，cheb2，ellip，besself求最小階數(shù)NButtord,

cheb1ordCheb2ord,ellipord144IIR數(shù)字濾波器設(shè)計(jì)流程圖145巴特沃斯濾波器設(shè)計(jì)流程（1）求最小階數(shù)N的函數(shù)buttord[N,wc]=buttord(wp,ws,Rp,Rs,‘s’)根據(jù)濾波器指標(biāo)wp，ws，Rp，Rs，求出巴特沃斯模擬濾波器的階數(shù)N及頻率參數(shù)wc，此處wp，ws及wc均以弧度/秒為單位。（2）得到N后,調(diào)用設(shè)計(jì)函數(shù)buttap[z，p，k]=buttap(N)得到[z,p,k]后，很容易求出濾波器系數(shù)B，A。（3）調(diào)用模擬頻率變換函數(shù)lp2lp[Bt,At]=lp2lp(B,A,wo)（4）調(diào)用模擬數(shù)字變換函數(shù)[Bd,Ad]=bilinear(B,A,Fs)145巴特沃斯濾波器設(shè)計(jì)流程（1）求最小階數(shù)N的函數(shù)butt146集成的數(shù)字濾波器設(shè)計(jì)函數(shù)把（2）、（3）、（4）合為一步的數(shù)字濾波器設(shè)計(jì)函數(shù)butter(N,wc,'ftype')[B,A]=butter(N,wc)設(shè)計(jì)低通或帶通數(shù)字濾波器系數(shù)B，A（當(dāng)為帶通濾波器時(shí)，第(1)類函數(shù)由wp=[wp1,wp2]會(huì)自動(dòng)生成wc=[w1,w2]）。

[B,A]=butter(N,wc,'high')設(shè)計(jì)高通數(shù)字濾波器系數(shù)B，A。

[B,A]=butter(N,wc,'stop')設(shè)計(jì)帶阻數(shù)字濾波器系數(shù)B，A。butter(N,wc,'ftype')還有零極增益和狀態(tài)空間形式，讀者可用help命令查閱。146集成的數(shù)字濾波器設(shè)計(jì)函數(shù)把（2）、（3）、（4）合為一147例7.29巴特沃斯模擬濾波器設(shè)計(jì)設(shè)計(jì)一個(gè)低通巴特沃斯模擬濾波器，指標(biāo)如下。通帶頻率：fp=3400Hz，最大衰減：Rp

=3dB阻帶頻率：fs

=4000Hz，最小衰減：As=40dB解:它的系統(tǒng)函數(shù)完全由階數(shù)N和3dB截止頻率c決定。而N和c是由濾波器設(shè)計(jì)指標(biāo)決定的。取c=c1，通帶指標(biāo)剛好，阻帶指標(biāo)富裕；取c=c2，則阻帶指標(biāo)剛好，通帶指標(biāo)富裕。MATLAB工具箱函數(shù)buttord，butter就是根據(jù)以上公式編寫(xiě)的。因此就無(wú)需再記憶這些公式了。

147例7.29巴特沃斯模擬濾波器設(shè)計(jì)設(shè)計(jì)一個(gè)低通巴特148模擬轉(zhuǎn)換為數(shù)字：脈沖響應(yīng)不變法模擬濾波器離散化的基本方法有脈沖響應(yīng)不變法和雙線性變換法。脈沖響應(yīng)不變法及impinvar函數(shù)單極點(diǎn)的N階模擬濾波器Ha(s)，用部分分式展開(kāi)為脈沖響應(yīng)不變法的數(shù)字化結(jié)果為工具箱函數(shù)impinvar可實(shí)現(xiàn)以上計(jì)算，格式為

[Bz,Az]=impinvar(B,A,Fs)148模擬轉(zhuǎn)換為數(shù)字：脈沖響應(yīng)不變法模擬濾波器離散化的基本方149模擬轉(zhuǎn)換為數(shù)字：雙線性變換法雙線性變換法函數(shù)bilinear雙線性變換法的原理是用 代換Ha(s)中的s值，得到H(z)。bilinear函數(shù)用來(lái)實(shí)現(xiàn)這個(gè)轉(zhuǎn)換。其使用格式為

[Bz,Az]=bilinear(B,A,Fs)脈沖響應(yīng)不變法的缺點(diǎn)是存在頻率混疊失真。雙線性變換法可完全消除頻率混疊失真，缺點(diǎn)是存在非線性頻率失真。149模擬轉(zhuǎn)換為數(shù)字：雙線性變換法雙線性變換法函數(shù)bil150例7.30模擬低通轉(zhuǎn)換為數(shù)字低通已知一模擬濾波器的系統(tǒng)函數(shù)為分別用脈沖響應(yīng)不變法和雙線性變換法將Ha(s)轉(zhuǎn)換成數(shù)字濾波器系統(tǒng)函數(shù)H(z)，并圖示Ha(s)和H(z)的幅頻響應(yīng)曲線。程序中的核心語(yǔ)句是以下兩條：[d,c]=impinvar(b,a,Fs) %用impinvar函數(shù)離散化[f,e]=bilinear(b,a,Fs) %用bilinear函數(shù)離散化150例7.30模擬低通轉(zhuǎn)換為數(shù)字低通已知一模擬濾波器的151例7.30模擬低通轉(zhuǎn)換為數(shù)字低通圖形結(jié)果如圖7.30所示。由圖7.30（b）可見(jiàn)，對(duì)脈沖響應(yīng)不變法，采樣頻率Fs越高（T越?。?，混疊越??；由圖7.30（c）可見(jiàn)，對(duì)雙線性變換法，無(wú)頻率混疊，但存在非線性失真。151例7.30模擬低通轉(zhuǎn)換為數(shù)字低通圖形結(jié)果如圖7.3152例7.31切比雪夫Ⅱ數(shù)字濾波器設(shè)計(jì)解：切比雪夫Ⅰ型濾波器通帶內(nèi)為等波紋，阻帶內(nèi)單調(diào)下降；切比雪夫Ⅱ型濾波器通帶內(nèi)為單調(diào)下降，阻帶內(nèi)等波紋。調(diào)用cheb2ord函數(shù)和cheby2函數(shù)使切比雪夫Ⅱ型設(shè)計(jì)變得非常簡(jiǎn)單。先用[N,wc]=Cheb2ord(wp,ws,Rp,Rs)求出N和wc，提供函數(shù)cheby2的輸入變?cè)?，再由[B,A]=cheby2(N,Rp,wc)設(shè)計(jì)切比雪夫Ⅱ型數(shù)字濾波器。B和A分別為H(z)的分子和分母多項(xiàng)式系數(shù)。對(duì)切比雪夫Ⅰ型濾波器，同樣有相應(yīng)的工具箱函數(shù)c