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

1第7章

MATLAB在數(shù)字信號處理中的應(yīng)用2書山有路勤為徑，學(xué)海無涯苦作舟。好好學(xué)習(xí)，天天向上3引言數(shù)字信號處理中包含離散時間信號與系統(tǒng)的一部分知識。本章主要研究內(nèi)容分為兩部分：離散時間信號與系統(tǒng)時域分析離散時間信號與系統(tǒng)頻域分析在數(shù)字信號處理中還有一個非常大的知識點：數(shù)字濾波器設(shè)計，留待下學(xué)期開設(shè)本課程時，可以學(xué)習(xí)。4一、離散時間信號與系統(tǒng)時域分析5包含內(nèi)容信號分析典型信號表示impseq(),stepseq()7.1序列相加，相乘7.2序列合成和截取7.3序列移位和周期延拓7.6離散序列卷積系統(tǒng)分析7.4系統(tǒng)響應(yīng)7.5系統(tǒng)線性判定67.1離散信號的產(chǎn)生及時域處理時域離散信號用x(n)表示（離散序列）。x和n同時使用才能完整地表示一個序列：時間變量n（表示采樣位置）只能取整數(shù)；向量x（表示采樣點的幅度）與n一一對應(yīng)。由于n序列是按整數(shù)遞增的,可簡單地用其初值ns決定,因為它的終值nf取決于ns和x的長度length(x),故可寫成：

n=[ns:nf]或 n=[ns:ns

length(x)

1]7典型信號表示impseq(),stepseq()8例7.1序列的相加和相乘給出兩個序列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。解：編程的思路是把序列長度延拓到覆蓋n1和n2的范圍，這樣才能把兩序列的時間變量對應(yīng)起來，然后進(jìn)行對應(yīng)元素的運算。9例程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)的時間變量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.*xa210結(jié)果顯示11例7.2序列的合成和截取用例6.13的結(jié)果編寫產(chǎn)生矩形序列RN(n)的程序。序列起點為n0，矩形序列起點為n1，長度為N（n0，n1，N由鍵盤輸入）。并用它截取一個復(fù)正弦序列exp(jπn/8).解：建模:矩形序列可看成兩個階躍序列之差。用MATLAB邏輯關(guān)系產(chǎn)生矩形序列x2(n)。而用它截取任何序列相當(dāng)于元素群相乘x2.*x，也稱為加窗運算。序列的合成和截取就是相加和相乘。12結(jié)果顯示13例7.2序列的合成和截取用例6.13的結(jié)果編寫產(chǎn)生矩形序列RN(n)的程序。序列起點為n0，矩形序列起點為n1，長度為N（n0，n1，N由鍵盤輸入）。并用它截取一個復(fù)正弦序列exp(jπn/8).解：建模:矩形序列可看成兩個階躍序列之差。用MATLAB邏輯關(guān)系產(chǎn)生矩形序列x2(n)。而用它截取任何序列相當(dāng)于元素群相乘x2.*x，也稱為加窗運算。序列的合成和截取就是相加和相乘。14程序要點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; %對復(fù)正弦序列加矩形窗（元素群乘）15程序顯示結(jié)果16例7.3序列的移位和周期延拓已知 ，利用MATLAB生成并圖示 表示x(n)以8為周期的延拓）和 解：方法1，利用矩陣乘法和冒號運算 x=[1234];y=x'*ones(1，3);方法2，采用求余函數(shù)mod，y=x(mod(n,M)+1)可實現(xiàn)對x(n)以M為周期的周期延拓。加1是因為MATLAB向量下標(biāo)只能從1開始，17程序?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);

end18周期化實現(xiàn)方法1：xc=xn(mod(n,8)+1); %產(chǎn)生x(n)的周期延拓

xcm=xn(mod(n-m,8)+1); %產(chǎn)生x(n)移位后的周期延拓19周期化實現(xiàn)方法2x=[1234];y=x'*ones(1,3);y=111222333444y1=y(:);%12341234123420程序結(jié)果21例7.6離散序列的卷積計算給出兩個序列 和 ,計算其卷積y(n)，并圖示各輸入輸出序列。解：在例6.4中，已經(jīng)給出了直接調(diào)用MATLAB的卷積函數(shù)conv的方法，也給出了自編卷積計算程序的方法，要注意的是本例時間變量的設(shè)定和移位方法。在本例中，設(shè)定n為從零開始，向量x和h的長度分別為Nx=20和Nh=10；結(jié)果向量y的長度為length(y)=Nx+Nh-1。22程序要點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)長度Ny2=length(y2);%y2(n)長度23程序顯示結(jié)果24例7.4離散系統(tǒng)對信號的響應(yīng)本題給定6階低通數(shù)字濾波器的系統(tǒng)函數(shù)，求它在下列輸入序列x(n)下的輸出序列y(n)。解：本題的計算原理見例6.14，在這里用工具箱函數(shù)filter來解。如果已知系統(tǒng)函數(shù)H(z)=B(z)/A(z)，則filter函數(shù)可求出系統(tǒng)對輸入信號x(n)的響應(yīng)y(n)。y

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

(1)

序列的正反Z變換

其中，符號表示取z變換，z是復(fù)變量。相應(yīng)地，單邊z變換定義為：Z變換與其逆變換33MATLAB符號數(shù)學(xué)工具箱提供了計算離散時間信號單邊z變換的函數(shù)ztrans和z反變換函數(shù)iztrans，其語句格式分別為Z=ztrans(x)x=iztrans(z)上式中的x和Z分別為時域表達(dá)式和z域表達(dá)式的符號表示，可通過sym函數(shù)來定義。MATLAB實現(xiàn)-方法1（符號數(shù)學(xué)工具箱）34

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

>>x=sym('a^n*cos(pi*n)');>>Z=ztrans(x);>>simplify(Z)ans=z/(z+a)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符號工具箱中的表示，反變換后的函數(shù)形式為：36例7.14用符號運算工具箱解z變換解：無限長度時間序列的z變換和逆z變換都屬于符號運算的范圍。MATLAB的symbolic（符號運算）工具箱已提供了這種函數(shù)。如果讀者已在計算機上安裝了這個工具箱，可以鍵入以下程序。MATLAB程序q714.m其特點是程序的開始要指定符號自變量symsznaNw0 %規(guī)定z,n,a…為符號變量37程序分析：symszn a Nw0%規(guī)定z,n,a為符號變量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á)式的形式38如果信號的z域表示式是有理函數(shù)，進(jìn)行z反變換的另一個方法是對進(jìn)行部分分式展開，然后求各簡單分式的z反變換.如果X(z)的有理分式表示為：

MATLAB實現(xiàn)-方法2（residuez函數(shù)實現(xiàn)反變換）39

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

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

3，終點nf=nf1+nf2=2。 n=ns:nf。由x和n即可得出X(z)。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) %求出n4647例7.8求z多項式分式的逆變換設(shè)系統(tǒng)函數(shù)為 ，輸入例7.7中的x2信號，用z變換計算輸出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變換，得48例7.8z多項式分式逆變換(續(xù))由程序算出nsy=-1，留數(shù)、極點分別為r=-57.7581和204.7581p=0.7791和0.3209k=-150-30代入得49Z反變換法求解系統(tǒng)輸出(編程實現(xiàn))Y(z)=X(z)W(z)x=[2,4,3,5];nsx=-2;%輸入序列及初始時間nfx=nsx+length(x)-1;%計算序列終止時間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變換的多項式乘積A=Aw;

%分母不變y(n)=IZT[Y(Z)][r,p,k]=residuez(B,A) %求留數(shù)r,極點p及直接項k nf=input('終點時間nf='); %要求鍵入終點時間n=min(nsx,nsy):nf; %生成總時間數(shù)組%求無限序列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;50另外一種方法(filter函數(shù))：filter(B,A,x)51例7.9離散時間傅里葉變換取周期的正弦信號，作8點采樣，求它的連續(xù)頻譜。然后對該信號進(jìn)行N個周期延拓，再求它的連續(xù)頻譜。把N無限增大，比較分析其結(jié)果。解：先求離散傅立葉變換的MATLAB子程序最后得到X=x*exp(-j*w*n‘)。有了子程序，本例就沒有什么難度了。52DTFT程序：function[X]=dtft(x,w)%計算離散時間傅立葉變換%[X]=dtft(x,n,w),%X=在w頻率點上的DTFT數(shù)組,%x=沿n的有限長度序列,%n=樣本位置向量%w=頻率點位置向量n=1:length(x);ewn=exp(-n'*w*i);X=x*ewn;53正弦信號x0的DTFT變換x0=sin(2*pi*[1:8]/8)*5; %x0是8點行向量dt=2*pi/8;w=linspace(-2*pi,2*pi,1000)/dt; %w是1000點行向量X0=dtft(x0,w)*dt; %求得頻率響應(yīng)X054正弦信號延拓N個周期后，DTFTN=4; %延拓周期數(shù)x1=reshape(x0'*ones(1,N),1,N*length(x0)); %延拓后的時域信號x1X1=dtft(x1,w)*dt; %求x1的頻率響應(yīng)X1%延拓100次后%x1=reshape(x0'*ones(1,100),1,100*length(x0)); %延拓后的時域信號x1%X1=dtft(x1,w)/100*dt; %求x1的頻率響應(yīng)X155重復(fù)無窮次disp('重復(fù)無窮次的八點信號的離散時間傅立葉變換-傅立葉級數(shù)')pause,X2=fft(x0*dt); %離散傅立葉變換w1=2*pi*[0:length(x0)-1]/length(x0); %離散頻點向量subplot(3,1,3),stem([-w1,w1],[abs(X2),abs(X2)]),grid,axis([min(w),max(w),0,max(abs(X2))]),grid56例7.9離散時間傅里葉變換2

程序運行結(jié)果執(zhí)行程序q709并按提示鍵入N=4，所得圖形如圖7.10所示。N取得愈大，其峰值愈大，寬度愈窄。當(dāng)N取得很大時，會出現(xiàn)內(nèi)存不足的問題，這是用矩陣乘法做傅里葉變換的缺點。另外，因為那時峰值點處的寬度很窄，也會出現(xiàn)所選頻點對不上峰值點的問題。所以對于N無限增大的情況，必須用fft函數(shù)來求。這時用連續(xù)頻譜也沒有意義了。這里用同樣的橫坐標(biāo)把幾種頻譜進(jìn)行對比，使讀者更好地理解其關(guān)系。57程序顯示結(jié)果：58例7.10時域采樣頻率與頻譜混疊分別以采樣頻率fs=1000Hz，400Hz和200Hz對xa(t)進(jìn)行等間隔采樣，計算并圖示三種采樣頻率下的采樣信號及其幅頻特性解：程序分別設(shè)定4種采樣頻率fs=10kHz，1kHz，400Hz和200Hz，對xa(t)進(jìn)行采樣，得到采樣序列xa(t)，xa1(n)，xa2(n)，xa3(n)，畫出其幅度頻譜。采樣時間區(qū)間均為0.1秒。為了便于比較，畫出了幅度歸一化的幅頻曲線，如圖7.11所示。59例7.10采樣頻率與頻譜混疊(續(xù))由于由以上關(guān)系式可見，采樣信號的頻譜函數(shù)是原模擬信號頻譜函數(shù)的周期延拓，延拓周期為2

/T。如果以頻率f為自變量（

=2

f），則以采樣頻率fs=1/T為延拓周期。對頻帶限于fc的模擬信號xa(t)，只有當(dāng)fs≥2fc時，采樣后 才不會發(fā)生頻譜混疊失真。這就是著名的采樣定理程序要點：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);%近似模擬信號頻譜60結(jié)果161結(jié)果26263例7.12梳狀濾波器零極點和幅特性梳狀濾波器系統(tǒng)函數(shù)有如下兩種類型。FIR型：IIR型：freqz數(shù)字濾波器頻率特性計算和繪制函數(shù)zplaneH(z)的零-極點圖繪制。解：調(diào)用函數(shù)freqz和zplane很容易寫出程序q712.m。64程序分析freqzfreqz(b,a,w)zplane65頻率響應(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);66零極點及繪圖subplot(2,2,1);zplane(b,a0);title('FIR梳狀濾波器零點圖');subplot(2,2,2);zplane(b,a1);title('IIR梳狀濾波器零、極點圖,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梳狀濾波器零、極點圖,a=0.9');subplot(2,2,2);zplane(b,a3);title('IIR梳狀濾波器零、極點圖,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('ω/π')67結(jié)果顯示68結(jié)果顯示69例7.13低通濾波及時域卷積定理輸入信號x(n)=cos(0.04

n)+cos(0.08

n)+cos(0.4

n)+0.3

(n)，0≤n≤63通過低通濾波器，計算濾波器對x(n)的響應(yīng)輸出y(n)，并圖示x(n)和y(n)，觀察濾波效果。解：如前所述，只要求出H(z)=B(z)/A(z)的分子和分母多項式系數(shù)向量B和A，則可調(diào)用濾波器直接Ⅱ型實現(xiàn)函數(shù)filter對輸入信號x(n)進(jìn)行濾波。

y=filter(B,A,x)70程序要點輸入信號構(gòu)造%產(chǎn)生輸入信號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;71程序要點%求H(z)分子分母多項式系數(shù)向量B和Ab=[1,2,1];%(1+z-1)2的展開系數(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);72程序要點濾波及圖形%對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('輸入信號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('輸出信號y(n)');xlabel('n');ylabel('y(n)')k=0:N-1;f=2*k/N;subplot(3,2,2);plot(f,abs(X))title('輸入信號x(n)的幅頻曲線');xlabel('ω/π');ylabel('|FT[x(n)]|')axis([0,0.5,0,max(abs(X))]);subplot(3,2,6);plot(f,abs(Y))title('輸出信號y(n)的幅頻曲線');xlabel('ω/π');ylabel('|FT[y(n)]|')axis([0,0.5,0,max(abs(Y))]);73結(jié)果顯示：74程序要點[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))]);75顯示結(jié)果76此處講一下-離散系統(tǒng)時頻域綜合分析771、系統(tǒng)函數(shù)的零極點分析離散時間系統(tǒng)的系統(tǒng)函數(shù)定義為系統(tǒng)零狀態(tài)響應(yīng)的z變換與激勵的z變換之比:

如果系統(tǒng)函數(shù)的有理函數(shù)表示式為離散系統(tǒng)時頻域綜合分析78在MATLAB中系統(tǒng)函數(shù)的零極點就可通過函數(shù)roots得到，也可借助DSP工具箱中的函數(shù)tf2zp得到，tf2zp的語句格式為：[R,P,K]=tf2zp(B,A)其中，B與A分別表示分子與分母多項式的系數(shù)向量。它的作用是將H(z)的有理分式表示式轉(zhuǎn)換為零極點增益形式：MATLAB實現(xiàn)離散系統(tǒng)時頻域綜合分析79離散系統(tǒng)時頻域綜合分析【例1】已知一離散因果LTI系統(tǒng)的系統(tǒng)函數(shù)為：試用MATLAB命令求該系統(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因此，零點為，極點為：。求解：81離散系統(tǒng)時頻域綜合分析若要獲得系統(tǒng)函數(shù)的零極點分布圖，可直接應(yīng)用zplane函數(shù)，其語句格式為：zplane(B,A)其中，B與A分別表示的分子和分母多項式的系數(shù)向量。它的作用是在Z平面上畫出單位圓、零點與極點。。82離散系統(tǒng)時頻域綜合分析【例2】已知一離散因果LTI系統(tǒng)的系統(tǒng)函數(shù)為：試用MATLAB命令繪出該系統(tǒng)的零極點分布圖。。83求解：。MATLAB源程序為：>>B=[1,0,-0.36];>>A=[1,-1.52,0.68];>>zplane(B,A),gridon>>legend('零點','極點')>>title('零極點分布圖')程序運行結(jié)果如圖所示。84離散系統(tǒng)時頻域綜合分析系統(tǒng)函數(shù)的零極點分布與其時域特性的關(guān)系：在離散系統(tǒng)中，z變換建立了時域函數(shù)與z域函數(shù)之間的對應(yīng)關(guān)系。因此，z變換的函數(shù)從形式可以反映的部分內(nèi)在性質(zhì)。我們通過討論H(z)的一階極點情況，來說明系統(tǒng)函數(shù)的零極點分布與系統(tǒng)時域特性的關(guān)系。85離散系統(tǒng)時頻域綜合分析MATLAB求解單位抽樣響應(yīng)可利用函數(shù)filter，

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

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

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

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

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

返回值w包含范圍內(nèi)的N個頻率等分點；返回值H則是離散時間系統(tǒng)頻率響應(yīng)。91離散系統(tǒng)時頻域綜合分析另一種形式為：[H,w]=freqz(B,A,N,’whole’)與第一種方式不同之處在于角頻率的范圍擴展到92離散系統(tǒng)時頻域綜合分析【例1】試用MATLAB命令繪制以下系統(tǒng)的頻率響應(yīng)曲線。解：利用函數(shù)freqz計算出然后利用函數(shù)abs和angle分別求出幅頻特性與相頻特性，最后利用plot命令繪出曲線。93離散系統(tǒng)時頻域綜合分析MATLAB源程序為：b=[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)相頻特性曲線')94結(jié)果顯示：95最后一個內(nèi)容DFT變換967.3離散傅里葉變換（DFT）定義DFT：用類似于例7.9中的方法，可把（7.3）式寫成矩陣乘法運算。其中，xn為序列行向量，Wnk是一N×N階方陣，而 稱為旋轉(zhuǎn)因子。977.3離散傅里葉變換（DFT）用矩陣乘法計算N點DFT的程序如下。

MATLAB程序q73a.m%用矩陣乘法計算N點DFTclear;closeallxn=input('請輸入序列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; %計算N點DFT98例7.15序列的離散傅立葉變換求復(fù)正弦序列

余弦序列

正弦序列的離散傅立葉變換，分別按N=16和N=8進(jìn)行計算。繪出幅頻特性曲線，進(jìn)行比較討論。99程序分析clear;closeallN=16;N1=8;%%產(chǎn)生序列x1(n),計算DFT[x1(n)]n=0:N-1;x1n=exp(j*pi*n/8);%產(chǎn)生x1(n)X1k=fft(x1n,N);%計算N點DFT[x1(n)]Xk1=fft(x1n,N1);%計算N1點DFT[x1(n)]%%產(chǎn)生序列x2(n),計算DFT[x2(n)]x2n=cos(pi*n/8);X2k=fft(x2n,N);%計算N點DFT[x2(n)]Xk2=fft(x2n,N1);%計算N1點DFT[x1(n)]%%產(chǎn)生序列x3(n),計算DFT[x3(n)]x3n=sin(pi*n/8);X3k=fft(x3n,N);%計算N點DFT[x3(n)]Xk3=fft(x3n,N1);%計算N1點DFT[x1(n)]%100例7.15序列的離散傅立葉變換在截取16點時，得到的是完整的余弦波形；而截取8點時，得到的是半截的余弦波形，當(dāng)然有大量的諧波成分。101本題分析DFT求解結(jié)果會分析102例7.16驗證N點DFT的物理意義（1） ，繪出幅頻曲線和相頻曲線。（2）計算并圖示x(n)的8點DFT。（3）計算并圖示x(n)的16點DFT。解:序列x(n)的Ｎ點DFT的物理意義是 在[0，2

]上進(jìn)行Ｎ點等間隔采樣。程序先密集采樣，繪制出幅頻曲線圖。然后再分別做8點和16點DFT來驗證這個采樣關(guān)系。結(jié)果1：103結(jié)果2104結(jié)果3105106對本章《MATLAB在數(shù)字信號處理中的應(yīng)用》的內(nèi)容講解，告一段落。107例7.10時域采樣頻率與頻譜混疊分別以采樣頻率fs=1000Hz，400Hz和200Hz對xa(t)進(jìn)行等間隔采樣，計算并圖示三種采樣頻率下的采樣信號及其幅頻特性解：程序分別設(shè)定4種采樣頻率fs=10kHz，1kHz，400Hz和200Hz，對xa(t)進(jìn)行采樣，得到采樣序列xa(t)，xa1(n)，xa2(n)，xa3(n)，畫出其幅度頻譜。采樣時間區(qū)間均為0.1秒。為了便于比較，畫出了幅度歸一化的幅頻曲線，如圖7.11所示。108例7.10采樣頻率與頻譜混疊(續(xù))由于由以上關(guān)系式可見，采樣信號的頻譜函數(shù)是原模擬信號頻譜函數(shù)的周期延拓，延拓周期為2

/T。如果以頻率f為自變量（

=2

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

其中 模擬用理想低通濾波器恢復(fù)的過程，觀察恢復(fù)波形，計算出最大恢復(fù)誤差。解：這個公式與卷積公式相像，可以用向量和矩陣乘法來解決。110例7.11由離散序列恢復(fù)模擬信號xa=x*g(TNM)=x*G其中G=sinc(Fs*TNM)M表示在兩個采樣點之間增加的間隔數(shù)，使輸出更密，更接近模擬信號。111例7.16驗證N點DFT的物理意義（1） ，繪出幅頻曲線和相頻曲線。（2）計算并圖示x(n)的8點DFT。（3）計算并圖示x(n)的16點DFT。解:序列x(n)的Ｎ點DFT的物理意義是 在[0，2

]上進(jìn)行Ｎ點等間隔采樣。程序先密集采樣，繪制出幅頻曲線圖。然后再分別做8點和16點DFT來驗證這個采樣關(guān)系。112例7.17頻域與時域采樣對偶性（1）產(chǎn)生三角波序列（2）對M=40，計算x(n)的64點DFT，并圖示x(n)和X(k)=DFT[x(n)]，k=0,1,…,63。（3）對（2）中所得X(k)在[0，2

]上進(jìn)行32點抽樣得（4）求 的32點IDFT，即

（5）繪出 的波形圖，評述它與x(n)的關(guān)系。113例7.17頻域與時域采樣對偶性由于頻域在[0，2

]上的采樣點數(shù)N（N=32）小于x(n)的長度M（M=40），所以，產(chǎn)生時域混疊現(xiàn)象，不能由X1(k)恢復(fù)原序列x(n)。只有滿足N≥M時，可由頻域采樣X1(k)得到原序列x(n)。這就是頻域采樣定理。對N≥M的情況，請讀者自己編程驗證。114例7.18快速卷積快速卷積就是根據(jù)DFT的循環(huán)卷積性質(zhì)，將時域卷積轉(zhuǎn)換為頻域相乘，最后再進(jìn)行IDFT得到時域卷積序列y(n)。其中時域和頻域之間的變換均用FFT實現(xiàn)，所以使卷積速度大大提高。框圖如下:

115例7.19用DFT求連續(xù)信號頻譜在計算機上用DFT對模擬信號進(jìn)行譜分析時，只能以有限大的采樣頻率fs對模擬信號采樣有限點樣本序列（等價于截取模擬信號一段進(jìn)行采樣）作DFT變換，得到模擬信號的近似頻譜。其誤差主要來自以下因素：①截斷效應(yīng)（頻譜泄露和譜間干擾）②頻譜混疊失真因素①使譜分辨率（能分辨開的兩根譜線間的最小間距）降低，并產(chǎn)生譜間干擾；因素②使折疊頻率（fs/2）附近的頻譜產(chǎn)生較大失真。116例7.19用DFT求連續(xù)信號頻譜加大截取長度Tp可提高頻率分辨率；選擇合適的窗函數(shù)可降低譜間干擾；而頻譜混疊失真要通過提高采樣頻率fs和（或）預(yù)濾波（在采樣之前濾除折疊頻率以外的頻率成分）來改善。編寫程序q719.m驗證截斷效應(yīng)及加窗的改善作用，先選取以下參數(shù)：

采樣頻率fs=400Hz，T=1/fs

采樣信號序列

對x(n)作4096點DFT作為的近似頻譜Xa(jf)。117例7.19用DFT求連續(xù)信號頻譜如圖7.19所示。圖中X1(jf)，X4(jf)和X8(jf)分別表示Tp=0.04s，0.04*4s和0.04*8s時的譜分析結(jié)果。由圖可見，由于截斷使原頻譜中的單頻譜線展寬（又稱之為泄漏），Tp越大泄漏越小，頻率分辨率越高。Tp=0.04s時，25Hz與50Hz兩根譜線已分辨不清了。所以實際譜分析的截取時間Tp是由頻率分辨率決定的。另外，在本應(yīng)為零的頻段上出現(xiàn)了一些參差不齊的小譜包（稱為譜間干擾）。譜間干擾的大小取決于加窗的類型。用矩形窗比用Hamming窗的頻率分辨率高（泄漏小），但譜間干擾剛好相反。118例7.20IIR濾波器直接型的轉(zhuǎn)換程序調(diào)用了信號處理工具箱函數(shù)tf2sos和擴展函數(shù)dir2par，

[sos,g]=tf2sos(B,A)實現(xiàn)從直接型到級聯(lián)型（二階分割形式）的轉(zhuǎn)換。g為式中的增益，sos為L×6階矩陣，表示式中的系數(shù)。119例7.20IIR濾波器直接型的轉(zhuǎn)換

[Cp,Bp,Ap]=dir2par(B,A)實現(xiàn)從直接型到并聯(lián)型的轉(zhuǎn)換。B為直接型H(z)的分子多項式系數(shù)向量，A為直接型H(z)的分母多項式系數(shù)向量；Cp，Bp，Ap的含義與擴展函數(shù)dir2par中的C，B，A相同。dir2par中又調(diào)用了復(fù)共軛對比較函數(shù)cplxcomp。由于dir2par和cplxcomp是文獻(xiàn)[7]中開發(fā)的，不屬于MATLAB工具箱函數(shù)，所以將其M文件清單附在程序q720.m之后。120例7.20IIR濾波器直接型的轉(zhuǎn)換根據(jù)計算結(jié)果,級聯(lián)型H(z)表達(dá)式：

級聯(lián)型結(jié)構(gòu)圖。121例7.20IIR濾波器直接型的轉(zhuǎn)換并聯(lián)型結(jié)構(gòu)H(z)表達(dá)式并聯(lián)型結(jié)構(gòu)圖122例7.21直接型結(jié)構(gòu)到格型梯形結(jié)構(gòu)

tf2latc函數(shù)實現(xiàn)直接型到格型轉(zhuǎn)換

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

K=tf2latc(1,A)求出全極點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ù)項B(1)歸一化）。123例7.21直接型結(jié)構(gòu)到格型梯形結(jié)構(gòu)直接型系統(tǒng)函數(shù)轉(zhuǎn)換為格型梯形結(jié)構(gòu)124例7.22FIR濾波器直接型到其他型系統(tǒng)函數(shù)為調(diào)用信號處理工具箱函數(shù)tf2sos和tf2latc,給變元A賦值1,B=[2,13/12,5/4,2/3]即可.級聯(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.3333125例7.22FIR濾波器直接型到其他型得出級聯(lián)型為格型結(jié)構(gòu)為126例7.23FIR格型到直接型轉(zhuǎn)換給定K=[2,1/4,1/2,1/3],用函數(shù)latc2tf即可由B=latc2tf(K)得到B,寫出直接型結(jié)構(gòu)127例7.24系統(tǒng)函數(shù)的計算機推導(dǎo)數(shù)字濾波器的網(wǎng)絡(luò)結(jié)構(gòu)圖實際上也是一種信號流圖。因此【例6.20】中介紹的方法和公式同樣可以用來求離散域的數(shù)字濾波器的系統(tǒng)函數(shù)。不同的地方僅僅在于節(jié)點方程中出現(xiàn)了作為系數(shù)的符號變量z

1，它將出現(xiàn)在系數(shù)矩陣中。MATLAB是不能處理上標(biāo)變量的，因此在程序中設(shè)q=z

1，在計算完成后再人工地把結(jié)果中的q恢復(fù)為z

1。128例7.24的結(jié)構(gòu)圖與方程129例7.24的方程的矩陣形式由此可以求出系統(tǒng)函數(shù)130例7.24解出的系統(tǒng)函數(shù)程序運行的結(jié)果如果加入一個激勵x(n)

A

(n)，則得出1317.5FIR數(shù)字濾波器設(shè)計濾波器的特性指標(biāo)用絕對值δ1,δ2表示;用分貝Rp,Rs表示132（1）窗函數(shù)法設(shè)計FIR濾波器先根據(jù)

c和N求出相應(yīng)的理想濾波器單位脈沖響應(yīng)hd(n)。第二步要選擇合適的窗函數(shù)w(n)來截取hd(n)的適當(dāng)長度（即階數(shù)），以保證實現(xiàn)要求的阻帶衰減；最后得到FIR濾波器單位脈沖響應(yīng)h(n)=hd(n).*w(n),即其系數(shù)。133（2）等波紋最佳一致逼近法（2）等波紋最佳一致逼近法:信號處理工具箱采用remez算法實現(xiàn)線性相位FIR數(shù)字濾波器的等波紋最佳一致逼近設(shè)計。其優(yōu)點是，設(shè)計指標(biāo)相同時，使濾波器階數(shù)最低；或階數(shù)相同時，使通帶最平坦，阻帶最小衰減最大；通帶和阻帶均為等波紋形式，最適合設(shè)計片段常數(shù)特性的濾波器。其調(diào)用格式如下：

b=remez(N,f,m,w,'ftype')其中N由remezord函數(shù)求出:

[N,fo,mo,w]=remezord(f,m,dev,Fs)輸入變元dev為各逼近頻段允許的波紋振幅。remez函數(shù)可直接調(diào)用remezord返回的參數(shù)如下：

b=remez(N,fo,mo,w)134例7.25窗函數(shù)法設(shè)計數(shù)字濾波器分別用矩形窗和Hamming窗設(shè)計線性相位FIR低通濾波器。要求通帶截止頻率

c=

/4，單位脈沖響應(yīng)h(n)的長度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)。135例7.25窗函數(shù)法設(shè)計數(shù)字濾波器本題中，

c=

/4，N=21，所以線性相位理想低通濾波器的單位脈沖響應(yīng)為:為了滿足線性相位FIR濾波器條件h(n)=h(N-1-n)，要求

=(N-1)/2=10。信號處理工具箱中有窗生成函數(shù)boxcar，hamming，hanning和blackman等。136例7.25窗函數(shù)法設(shè)計數(shù)字濾波器對兩種窗函數(shù)的設(shè)計結(jié)果分別如右圖7.25-1和圖7.25-2所示。137工具箱設(shè)計函數(shù)fir1和fir2MATLAB提供了基于窗函數(shù)法的FIR濾波器設(shè)計函數(shù)fir1和fir2,其功能及用法如下。

fir1功能：標(biāo)準(zhǔn)頻率響應(yīng)形狀。格式：b=fir1(N,wc,‘ftype’,window)。當(dāng)wc=[wc1，wc2]時，是的帶通濾波器。當(dāng)ftype=high時，設(shè)計高通FIR濾波器；當(dāng)ftype=stop時，設(shè)計帶阻FIR濾波器。

fir2功能：任意頻率響應(yīng)形狀。格式：b=fir2(N,f,m,window)138例7.26窗函數(shù)法設(shè)計帶通濾波器使用fir1函數(shù)b=fir1(N,wc,window)編程參數(shù)

c為行向量

c

=[

lp/

，

hp/

]根據(jù)阻帶最小衰減Rs

=60dB選擇窗函數(shù)類型和階次。可以查上面列出的“窗函數(shù)設(shè)計濾波器時的階數(shù)選擇表”。選blackman窗，其濾波器阻帶最小衰減可達(dá)到74dB，其窗口長度M由過渡帶寬度B=0.15

決定，Blackman窗設(shè)計的濾波器過渡帶寬度為12

/M，故M取80。因M=N+1，所以濾波器階數(shù)N=79。139例7.27用remez函數(shù)低通濾波器解：先由題意計算設(shè)計參數(shù)f=[1/4,5/16]，m=[1,0];dev的計算稍復(fù)雜一些，由于所以有了這幾個參數(shù)就可以調(diào)用remezord和remez函數(shù)了.140例7.27用remez函數(shù)低通濾波器橫線為-3dB，兩條豎線分別位于頻率

/4和5

/16。顯然，通帶指標(biāo)稍有富裕，過渡帶寬度和阻帶最小衰減剛好滿足指標(biāo)要求。程序輸出的幅頻特性141例7.28remez函數(shù)設(shè)計高通濾波器觀察等波紋逼近法中加權(quán)系數(shù)w(

)及濾波器階數(shù)N的作用和影響。期望逼近的濾波器通帶為[3

/4，

]，阻帶為[0，23

/32]。解:在濾波器設(shè)計中，技術(shù)指標(biāo)越高，實現(xiàn)濾波器的階數(shù)也就越高。另外，對固定的階數(shù)，通帶與阻帶指標(biāo)可以互換，過渡帶寬度與通帶波紋和阻帶衰減指標(biāo)可以互換。取f=[0,3/4,23/32,1]，m=[0,0,1,1]。其余參數(shù)分三種情況進(jìn)行設(shè)計，①N=30，w=[1,1]；②N=30，w=[1,5];③N=60，w=[1,1]。142例7.28remez函數(shù)設(shè)計高通濾波器程序運行結(jié)果如圖由圖可見，w較大的頻段逼近精度較高；w較小的頻段逼近精度較低。N較大時逼近精度較高,N較小時逼近精度較低。1437.6IIR數(shù)字濾波器設(shè)計IIR數(shù)字濾波器設(shè)計的主要方法是先設(shè)計低通模擬濾波器，進(jìn)行頻率變換，將其轉(zhuǎn)換為相應(yīng)的（高通、帶通等）模擬濾波器，再轉(zhuǎn)換為高通、帶通或帶阻數(shù)字濾波器。對設(shè)計的全過程的各個步驟，MATLAB都提供了相應(yīng)的工具箱函數(shù)，使IIR數(shù)字濾波器設(shè)計變得非常簡單。本節(jié)主要結(jié)合例題介紹這些IIR濾波器設(shè)計的工具箱函數(shù)。IIR數(shù)字濾波器的設(shè)計步驟由以下的流程圖來表示。下面以巴特沃斯濾波器設(shè)計函數(shù)為典型，介紹此流程圖中函數(shù)的功能和用法。144IIR數(shù)字濾波器設(shè)計流程圖模擬低通濾波器原型設(shè)計Buttap,cheb1ap,cheb2apbesselap，ellipap函數(shù)頻率變換(變?yōu)楦咄ǎ瑤?，帶阻?lp2lp，lp2hp，lp2bp，lp2bs模擬數(shù)字變換bilinearimpinvar合為一步的設(shè)計函數(shù)butter，cheb1，cheb2，ellip，besself求最小階數(shù)NButtord,

cheb1ordCheb2ord,ellipord145巴特沃斯濾波器設(shè)計流程（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è)計函數(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)146集成的數(shù)字濾波器設(shè)計函數(shù)把（2）、（3）、（4）合為一步的數(shù)字濾波器設(shè)計函數(shù)bu