基于MATLAB數(shù)字信號(hào)處理實(shí)驗(yàn)

1、第4章 MATLAB在信號(hào)處理(xn ho ch l)中的應(yīng)用 4.1 信號(hào)及其表示4.2 信號(hào)的基本(jbn)運(yùn)算 4.3 信號(hào)的能量和功率 4.4 線性時(shí)不變系統(tǒng)4.5 線性時(shí)不變系統(tǒng)的響應(yīng)4.6 線性時(shí)不變系統(tǒng)的頻率響應(yīng) 4.7 傅里葉(Fourier)變換4.8 IIR數(shù)字濾波器的設(shè)計(jì)方法4.9 FIR數(shù)字濾波器設(shè)計(jì)共四十一頁4.1 信號(hào)(xnho)及其表示4.1.1連續(xù)時(shí)間(shjin)信號(hào)的表示 連續(xù)時(shí)間信號(hào)：時(shí)間變化連續(xù)。如y=x(t) 離散時(shí)間信號(hào)(序列)：時(shí)間離散，如x(nT)=x(t)|t=nT.4.1.2工具箱中的信號(hào)產(chǎn)生函數(shù)函數(shù)名功能函數(shù)名功能sawtooth產(chǎn)生鋸齒波

2、或三角波信號(hào)pulstran產(chǎn)生沖激串square產(chǎn)生方波信號(hào)rectpule產(chǎn)生非周期的方波信號(hào) sinc產(chǎn)生sinc函數(shù)波形tripuls產(chǎn)生非周期的三角波信號(hào)chirp產(chǎn)生調(diào)頻余弦信號(hào)diric產(chǎn)生Dirichlet或周期sinc函數(shù)gauspuls產(chǎn)生高斯正弦脈沖信號(hào)gmonopuls產(chǎn)生高斯單脈沖信號(hào)vco電壓控制振蕩器共四十一頁4.1.3離散(lsn)時(shí)間信號(hào)的表示 在MATLAB中，離散時(shí)間信號(hào)x(n)的表示(biosh)：需用一個(gè)向量x表示(biosh)序列幅值，用另一個(gè)等長(zhǎng)的定位時(shí)間變量n，才能完整地表示(biosh)一個(gè)序列。 例4-10 繪制離散時(shí)間信號(hào)的棒狀圖。其中x(

3、-1)=-1, x(0)=1, x(1)=2, x(2)=1, x(3)=0, x(4)=-1。MATLAB源程序?yàn)椋簄=-3:5; %定位時(shí)間變量x=0,0,-1,1,2,1,-1,0,0;stem(n,x); grid; % 繪制棒狀圖line(-3,5,0,0); %畫x軸線xlabel(n); ylabel(xn)運(yùn)行結(jié)果如圖4.11所示。圖 4.11 離散時(shí)間信號(hào)圖形共四十一頁4.1.4幾種常用(chn yn)離散時(shí)間信號(hào)的表示1單位(dnwi)脈沖序列直接實(shí)現(xiàn)：x=zeros(1,N); x(1,n0)=1;2單位階躍序列 直接實(shí)現(xiàn)：n=ns:nf; x=(n-n0)=0;共四十一

4、頁3實(shí)指數(shù)(zhsh)序列直接(zhji)實(shí)現(xiàn)：n=ns:nf; x=a.n;4復(fù)指數(shù)序列直接實(shí)現(xiàn)：n=ns:nf; x=exp(sigema+jw)*n);5正(余)弦序列直接實(shí)現(xiàn)：n=ns:nf; x=cos(w*n+sita);共四十一頁4.2 信號(hào)的基本(jbn)運(yùn)算4.2.1信號(hào)(xnho)的相加與相乘 y(n)=x1(n)+x2(n) y(n)=x1(n)x2(n) MATLAB實(shí)現(xiàn)：y=x1+x2; y=x1.*x24.2.2序列移位與周期延拓運(yùn)算序列移位：y(n)=x(n-m)。MATLAB實(shí)現(xiàn)：y=x; ny=nx-m序列周期延拓：y(n)=x(n)M，MATLAB實(shí)現(xiàn)：ny

5、=nxs:nxf；y=x(mod(ny,M)+1)4.2.3 序列翻褶與序列累加運(yùn)算序列翻褶：y(n)=x(-n)。MATLAB可實(shí)現(xiàn): y=fliplr(x)序列累加的數(shù)學(xué)描述為： MATLAB實(shí)現(xiàn)：y=cumsum(x)共四十一頁4.2.4 兩序列(xli)的卷積運(yùn)算兩序列(xli)卷積運(yùn)算： MATLAB實(shí)現(xiàn)：y=conv(x1,x2)。序列x1(n)和x2(n)必須長(zhǎng)度有限。 4.2.5 兩序列的相關(guān)運(yùn)算兩序列相關(guān)運(yùn)算： 。MATLAB實(shí)現(xiàn)：y=xcorr(x1,x2)。共四十一頁4.3 信號(hào)的能量(nngling)和功率1.信號(hào)(xnho)能量數(shù)字定義：MATLAB實(shí)現(xiàn): E=sum

6、(x.*conj(x); 或 E=sum(abs(x).2);數(shù)字定義：2. 信號(hào)功率MATLAB實(shí)現(xiàn): P=sum(x.*conj(x)/N; 或 E=sum(abs(x).2)/N;共四十一頁4.4 線性時(shí)不變系統(tǒng)(xtng)4.4.1 系統(tǒng)(xtng)的描述1常系數(shù)線性微分/差分方程2系統(tǒng)傳遞函數(shù)3零極點(diǎn)增益模型連續(xù)系統(tǒng)： 連續(xù)系統(tǒng)： 離散系統(tǒng)： 離散系統(tǒng)： 共四十一頁4極點(diǎn)(jdin)留數(shù)模型離散系統(tǒng)： 連續(xù)(linx)系統(tǒng)： 5二次分式模型連續(xù)系統(tǒng)： 離散系統(tǒng)： 6狀態(tài)空間模型連續(xù)系統(tǒng)： 離散系統(tǒng)： 共四十一頁4.4.2 系統(tǒng)模型的轉(zhuǎn)換(zhunhun)函數(shù) 在MATLAB中，用so

7、s、ss、tf、zp分別表示二次分式模型、狀態(tài)空間模型、傳遞函數(shù)模型和零極點(diǎn)(jdin)增益模型。其中sos表示二次分式，g為比例系數(shù)，sos為L(zhǎng)6的矩陣，即 （415） 1ss2tf函數(shù)格式：num, den=ss2tf(A,B,C,D,iu)功能：將指定輸入量iu的線性系統(tǒng)(A,B,C,D)轉(zhuǎn)換為傳遞函數(shù)模型num,den。2zp2tf函數(shù)格式：num,den=zp2tf(z,p,k)功能：將給定系統(tǒng)的零極點(diǎn)增益模型轉(zhuǎn)換為傳遞函數(shù)模型，z、p、k分別為零點(diǎn)列向量、極點(diǎn)列向量和增益系數(shù)。 共四十一頁線性系統(tǒng)模型(mxng)的變換函數(shù)函數(shù)名功能說明函數(shù)名功能說明ss2tf 狀態(tài)空間模型轉(zhuǎn)換為傳

8、遞函數(shù)模型zp2tf零極點(diǎn)增益模型轉(zhuǎn)換為傳遞函數(shù)模型ss2zp狀態(tài)空間模型轉(zhuǎn)換為零極點(diǎn)增益模型zp2ss零極點(diǎn)增益模型轉(zhuǎn)換為狀態(tài)空間模型ss2sos狀態(tài)空間模型轉(zhuǎn)換為二次分式模型zp2sos零極點(diǎn)增益模型轉(zhuǎn)換為二次分式模型tf2ss傳遞函數(shù)模型轉(zhuǎn)換為狀態(tài)空間模型sos2tf二次分式模型轉(zhuǎn)換為傳遞函數(shù)模型tf2zp傳遞函數(shù)模型轉(zhuǎn)換為零極點(diǎn)增益模型sos2zp二次分式模型轉(zhuǎn)換為零極點(diǎn)增益模型tf2sos傳遞函數(shù)模型轉(zhuǎn)換為二次分式模型sos2ss二次分式模型轉(zhuǎn)換為狀態(tài)空間模型共四十一頁例4-18 求離散(lsn)時(shí)間系統(tǒng)的零、極點(diǎn)向量和增益系數(shù)。在命令(mng lng)窗口輸入： num=2,3;

9、den=1,0.4,1; num,den=eqtflength(num,den);%使長(zhǎng)度相等 z,p,k=tf2zp(num,den)屏幕顯示為z = 0 -1.5000p = -0.2000 + 0.9798i -0.2000 - 0.9798ik = 2共四十一頁4.4.3 系統(tǒng)(xtng)互聯(lián)與系統(tǒng)(xtng)結(jié)構(gòu)MATLAB實(shí)現(xiàn)(shxin)函數(shù)series( ) 格式：A,B,C,D=series(A1,B1,C1,D1,A2,B2,C2,D2) 或 num,den=series(num1,den1,num2,den2)將系統(tǒng)1、系統(tǒng)2級(jí)聯(lián)，可得到級(jí)聯(lián)連接的傳遞函數(shù)形式為：1. 系

10、統(tǒng)的級(jí)聯(lián)共四十一頁MATLAB實(shí)現(xiàn)(shxin)函數(shù)parallel( )格式：A,B,C,D=parallel(A1,B1,C1,D1,A2,B2,C2,D2) 或 num,den=parallel(num1,den1,num2,den2)2. 系統(tǒng)(xtng)的并聯(lián)將系統(tǒng)1、系統(tǒng)2并聯(lián)，可得到并聯(lián)連接的傳遞函數(shù)形式為：3. 兩個(gè)系統(tǒng)的反饋連接函數(shù)feedback格式：A,B,C,D=feedback(A1,B1,C1,D1,A2,B2,C2,D2,sign) 或 num,den=feedback(num1,den1,num2,den2,sign)將系統(tǒng)1和系統(tǒng)2進(jìn)行反饋連接，sign表示反

11、饋方式(默認(rèn)值為-1)； 當(dāng)sig=+1時(shí)表示正反饋；當(dāng)sig=-1時(shí)表示負(fù)反饋。共四十一頁例4-19 求兩個(gè)單輸入(shr)單輸出子系統(tǒng)的級(jí)聯(lián)、并聯(lián)和反饋(fnku)后系統(tǒng)的傳遞函數(shù)。MATLAB源程序?yàn)椋簄um1=1; den1=1,1; %系統(tǒng)1num2=2; den2=1,2; %系統(tǒng)2nums,dens=series(num1,den1,num2,den2) %實(shí)現(xiàn)兩個(gè)系統(tǒng)級(jí)聯(lián)nump,denp=parallel(num1,den1,num2,den2) %實(shí)現(xiàn)兩個(gè)系統(tǒng)并聯(lián) numf,denf=feedback(num1,den1,num2,den2) %實(shí)現(xiàn)兩個(gè)系統(tǒng)反饋程序運(yùn)行結(jié)果

12、為：nums = 0 0 2 ； dens = 1 3 2nump = 0 3 4 ； denp = 1 3 2numf = 0 1 2 ； denf = 1 3 4因此，各系統(tǒng)的傳遞函數(shù)分別為：共四十一頁4.5 線性時(shí)不變系統(tǒng)(xtng)的響應(yīng)4.5.1線性時(shí)不變系統(tǒng)(xtng)的時(shí)域響應(yīng)1連續(xù)LTI系統(tǒng)的響應(yīng)2離散LTI系統(tǒng)的響應(yīng)用MATLAB中的卷積函數(shù)conv( )來實(shí)現(xiàn)。 用MATLAB中的卷積函數(shù)conv( )來實(shí)現(xiàn)。 共四十一頁格式：y,x=lsim(a,b,c,d,u,t)功能(gngnng)：返回連續(xù)LTI系統(tǒng) (2)對(duì)任意輸入的離散LTI系統(tǒng)響應(yīng)函數(shù)dlsim( )格式(g

13、 shi)：y,x=dlsim(a,b,c,d,u)功能：返回離散LTI系統(tǒng) 對(duì)任意輸入時(shí)系統(tǒng)的輸出響應(yīng)y和狀態(tài)記錄x，其中u給出每個(gè)輸入的時(shí)序列，一般情況下u為一個(gè)矩陣；t用于指定仿真的時(shí)間軸，它應(yīng)為等間隔。對(duì)輸入序列u的響應(yīng)y和狀態(tài)記錄x。3時(shí)域響應(yīng)函數(shù)(1)對(duì)任意輸入的連續(xù)LTI系統(tǒng)響應(yīng)函數(shù)lsim( )共四十一頁4.5.2 LTI系統(tǒng)(xtng)的單位沖激響應(yīng)1. 求連續(xù)LTI系統(tǒng)的單位沖激響應(yīng)函數(shù)impulse( )格式：Y,T = impulse(sys)或impulse(sys)功能：返回系統(tǒng)的響應(yīng)Y和時(shí)間向量T，自動(dòng)選擇仿真的時(shí)間范圍。其中(qzhng)sys可為系統(tǒng)傳遞函數(shù)、

14、零極增益模型或狀態(tài)空間模型。2. 求離散系統(tǒng)的單位沖激響應(yīng)函數(shù)dimpulse( )格式：y,x=dimpulse(num,den)功能：返回項(xiàng)式傳遞函數(shù)的單位沖激響應(yīng)y向量和時(shí)間狀態(tài)歷史記錄x向量。共四十一頁4.5.3 時(shí)域響應(yīng)的其它函數(shù)1. 求連續(xù)LTI系統(tǒng)的零輸入響應(yīng)函數(shù)initial( )格式：y,t,x=initial(a,b,c,d,x0)功能：計(jì)算出連續(xù)時(shí)間LTI系統(tǒng)由于初始狀態(tài)x0所引起的零輸入響應(yīng)y。其中x為狀態(tài)記錄，t為仿真所用的采樣(ci yn)時(shí)間向量。2. 求離散系統(tǒng)的零輸入響應(yīng)函數(shù)dinitial( )格式：y,x,n=dinitial(a,b,c,d,x0)功能：

15、計(jì)算離散(lsn)時(shí)間LTI系統(tǒng)由初始狀態(tài)x0所引起的零輸入響應(yīng)y和狀態(tài)響應(yīng)響應(yīng)x，取樣點(diǎn)數(shù)由函數(shù)自動(dòng)選取。n為仿真所用的點(diǎn)數(shù)。3. 求連續(xù)系統(tǒng)的單位階躍響應(yīng)函數(shù)step( )格式：Y,T = step(sys)功能：返回系統(tǒng)的單位階躍響應(yīng)Y和仿真所用的時(shí)間向量T，自動(dòng)選擇仿真的時(shí)間范圍。其中sys可為系統(tǒng)傳遞函數(shù)(TF)、零極增益模型(ZPK)或狀態(tài)空間模型(SS)。4. 求離散系統(tǒng)的單位階躍響應(yīng)函數(shù)dstep( )格式：y,x= dstep (num,den)功能：返回多項(xiàng)式傳遞函數(shù)G(z)=num(z)/den(z)表示的系統(tǒng)單位階躍響應(yīng)。共四十一頁4.6線性時(shí)不變系統(tǒng)(xtng)的頻率

16、響應(yīng)1求模擬濾波器Ha(s)的頻率響應(yīng)函數(shù)freqs( )格式：Hfreqs(B,A,W) 功能：計(jì)算(j sun)由向量W(rad/s)指定的頻率點(diǎn)上模擬濾器系統(tǒng)函數(shù)Ha(s)的頻率響應(yīng)Ha(j)，結(jié)果存于H向量中。例4-31 已知某模擬濾波器的系統(tǒng)函數(shù)求該模擬濾波器的頻率響應(yīng)。MATLAB源程序如下。B=1;A=1 2.6131 3.4142 2.6131 1;W=0:0.1:2*pi*5;freqs(B,A,W)圖4.30 模擬濾波器的頻率響應(yīng)共四十一頁 例4-32 已知某濾波器的系統(tǒng)(xtng)函數(shù)為求該濾波器的頻率響應(yīng)(pn l xin yn)。MATLAB源程序?yàn)椋築=1 0 0

17、0 0 0 0 0 1;A=1;freqz(B,A)該程序運(yùn)行所繪出的幅頻與相頻性曲線如圖4.31所示。圖4.31濾波器幅度和相位曲線 2求數(shù)字濾波器H(z)的頻率響應(yīng)函數(shù)freqz( )格式：H=freqz(B,A,W)功能：計(jì)算由向量W(rad)指定的數(shù)字頻率點(diǎn)上(通常指在H(z)的頻率響應(yīng)H(ejw )。 共四十一頁3濾波函數(shù)filter格式：y=filter(B,A,x)功能：對(duì)向量x中的數(shù)據(jù)進(jìn)行濾波處理，即差分方程求解，產(chǎn)生輸出序列向量y。B和A分別為數(shù)字濾波器系統(tǒng)函數(shù)H(z)的分子(fnz)和分母多項(xiàng)式系數(shù)向量。例4-33 設(shè)系統(tǒng)(xtng)差分方程為MATLAB源程序?yàn)椋築=1;

18、 A=1,-0.8;N=0:31; x=0.8.n;y=filter(B,A,x);subplot(2,1,1);stem(x)subplot(2,1,2);stem(y)該程序運(yùn)行所得結(jié)果如圖4.32所示。，求該系統(tǒng)對(duì)信號(hào)的響應(yīng)。圖4.32系統(tǒng)對(duì)信號(hào)的響應(yīng) 共四十一頁4.7傅里葉(Fourier)變換(binhun)4.7.1連續(xù)時(shí)間(shjin)、連續(xù)頻率傅里葉變換4.7.2 連續(xù)時(shí)間、離散頻率傅里葉級(jí)數(shù)正變換： 逆變換： 正變換： 逆變換： 共四十一頁4.7.3 時(shí)間離散、連續(xù)頻率(pnl)序列傅里葉變換4.7.4 離散時(shí)間(shjin)、離散頻率離散傅里葉級(jí)數(shù)4.7.5離散時(shí)間、離散頻

19、率離散傅里葉變換(DFT)正變換： 逆變換： 正變換： 逆變換： 正變換： 逆變換： 共四十一頁1一維快速正傅里葉變換函數(shù)fft格式：X=fft(x, N)功能：采用(ciyng)FFT算法計(jì)算序列向量x的N點(diǎn)DFT變換， 當(dāng)N缺省時(shí)，fft函數(shù)自動(dòng)按x的長(zhǎng)度計(jì)算DFT。當(dāng)N為2整數(shù)次冪時(shí)，fft按基-2算法計(jì)算，否則用混合算法。2一維快速逆傅里葉變換函數(shù)ifft格式：x=ifft(X, N)功能：采用FFT算法計(jì)算序列向量X的N點(diǎn)IDFT變換。例4-36 用快速傅里葉變換FFT計(jì)算下面(xi mian)兩個(gè)序列的卷積。, 并測(cè)試直接卷積和快速卷積的時(shí)間。圖4.35 快速卷積框圖共四十一頁MA

20、TLAB程序(部分)：%線性卷積xn= sin(0.4*1:15);%對(duì)序列(xli)x(n)賦值, M=15hn= 0.9.(1:20); %對(duì)序列h(n)賦值, N=20yn=conv(xn,hn); % 直接調(diào)用函數(shù)conv計(jì)算卷積%園周卷積L=pow2(nextpow2(M+N-1); Xk=fft(xn,L);Hk=fft(hn,L);Yk=Xk.*Hk; yn=ifft(Yk,L); 圖4.36 x(n),h(n)及其線性卷積波形(b xn)共四十一頁4.8 IIR數(shù)字(shz)濾波器的設(shè)計(jì)方法1. 數(shù)字(shz)濾波器的頻率響應(yīng)函數(shù)幅度響應(yīng)：相位響應(yīng)：圖4.37 理想低通、高通、

21、帶通、帶阻數(shù)字濾波器幅度特性 共四十一頁2. 濾波器的技術(shù)指標(biāo) 幅度響應(yīng)(xingyng)指標(biāo)、相位響應(yīng)(xingyng)指標(biāo) 圖4.38 數(shù)字(shz)低通濾波器的幅度特性 通帶要求： 阻帶要求： 通帶最大衰減：阻帶最小衰減：共四十一頁4.8.1沖激響應(yīng)不變法(bin f)2.MATLAB信號(hào)處理工箱中的專用函數(shù)impinvar( )：格式：BZ,AZ =impinvar(B,A,Fs) 功能：把具有B,A模擬(mn)濾波器傳遞函數(shù)模型轉(zhuǎn)換成采樣頻率為Fs(Hz)的數(shù)字濾波器的傳遞函數(shù)模型BZ,AZ。采樣頻率Fs的默認(rèn)值為Fs=1。1. 沖激響應(yīng)不變法設(shè)計(jì)IIR數(shù)字濾波器的基本原理： 例4-

22、37 MATLAB源程序如下：num=1; %模擬濾波器系統(tǒng)函數(shù)的分子den=1,sqrt(5),2,sqrt(2),1; %模擬濾波器系統(tǒng)函數(shù)的分母num1,den1=impinvar(num,den) %求數(shù)字低通濾波器的系統(tǒng)函數(shù)程序的執(zhí)行結(jié)果如下：num1 = -0.0000 0.0942 0.2158 0.0311den1 = 1.0000 -2.0032 1.9982 -0.7612 0.1069共四十一頁MATLAB信號(hào)處理(xn ho ch l)工具箱中的專用雙線性變換函數(shù)bilinear( )格式：numd,dendbilinear(num,den,Fs)功能：把模擬濾波器的傳

23、遞函數(shù)模型轉(zhuǎn)換成數(shù)字濾波器的傳遞函數(shù)模型。4.8.2雙線性變換法雙線性變換(binhun)利用頻率變換(binhun)關(guān)系： 例4-38 MATLAB源程序如下： num=1; %模擬濾波器系統(tǒng)函數(shù)的分子 den=1,sqrt(3),sqrt(2),1; %模擬濾波器系統(tǒng)函數(shù)的分母 num1,den1=bilinear(num,den,1) %求數(shù)字濾波器的傳遞函數(shù)運(yùn)算的結(jié)果如下：num1 = 0.0533 0.1599 0.1599 0.0533den1 = 1.0000 -1.3382 0.9193 -0.1546共四十一頁4.8.3 IIR數(shù)字(shz)濾波器的頻率變換設(shè)計(jì)法1. IIR

24、數(shù)字濾波器的頻率變換設(shè)計(jì)法的基本原理 根據(jù)濾波器設(shè)計(jì)要求，設(shè)計(jì)模擬原型低通濾波器，然后(rnhu)進(jìn)行頻率變換，將其轉(zhuǎn)換為相應(yīng)的模擬濾波器(高通、帶通等)，最后利用沖激響應(yīng)不變法或雙線性變換法，將模擬濾波器數(shù)字化成相應(yīng)的數(shù)字濾波器。 圖4.39 IIR數(shù)字濾波器MATLAB設(shè)計(jì)步驟流程圖 共四十一頁1MATLAB的典型(dinxng)設(shè)計(jì)利用在MATLAB設(shè)計(jì)IIR數(shù)字(shz)濾波器可分以下幾步來實(shí)現(xiàn) (1)按一定規(guī)則將數(shù)字濾波器的技術(shù)指標(biāo)轉(zhuǎn)換為模擬低通濾波器的技術(shù)指標(biāo)；(2)根據(jù)轉(zhuǎn)換后的技術(shù)指標(biāo)使用濾波器階數(shù)函數(shù)，確定濾波器的最小階數(shù)N和截止頻率Wc；(3)利用最小階數(shù)N產(chǎn)生模擬低通濾波原

25、型；(4)利用截止頻率Wc把模擬低通濾波器原型轉(zhuǎn)換成模擬低通、高通、帶通或帶阻濾波器；(5)利用沖激響應(yīng)不變法或雙線性不變法把模擬濾波器轉(zhuǎn)換成數(shù)字濾波器。例4-39 設(shè)計(jì)一個(gè)數(shù)字信號(hào)處理系統(tǒng)，它的采樣率為Fs100Hz，希望在該系統(tǒng)中設(shè)計(jì)一個(gè)Butterworth型高通數(shù)字濾波器，使其通帶中允許的最小衰減為0.5dB，阻帶內(nèi)的最小衰減為40dB，通帶上限臨界頻率為30Hz，阻帶下限臨界頻率為40Hz。 共四十一頁MATLAB源程序設(shè)計(jì)如下： %把數(shù)字濾波器的頻率特征轉(zhuǎn)換成模擬濾波器的頻率特征 wp=30*2*pi;ws=40*2*pi;rp=0.5;rs=40;Fs=100; N,Wc=but

26、tord(wp,ws,rp,rs,s); %選擇濾波器的最小階數(shù) Z,P,K=buttap(N); %創(chuàng)建Butterworth低通濾波器原型 A,B,C,D=zp2ss(Z,P,K); %零極點(diǎn)增益模型轉(zhuǎn)換為狀態(tài)(zhungti)空間模型 AT,BT,CT,DT=lp2hp(A,B,C,D,Wc); %實(shí)現(xiàn)低通向高通的轉(zhuǎn)變 num1,den1=ss2tf(AT,BT,CT,DT); %狀態(tài)空間模型轉(zhuǎn)換為傳遞函數(shù)模型 %運(yùn)用雙線性變換法把模擬濾波器轉(zhuǎn)換成數(shù)字濾波器 num2,den2=bilinear(num1,den1,100); H,W=freqz(num2,den2); %求頻率響應(yīng) p

27、lot(W*Fs/(2*pi),abs(H);grid; %繪出頻率響應(yīng)曲線 xlabel(頻率/Hz);ylabel(幅值)程序運(yùn)行結(jié)果如圖4.40所示。 共四十一頁2MATLAB的直接(zhji)設(shè)計(jì)圖4.39 IIR數(shù)字(shz)濾波器MATLAB設(shè)計(jì)步驟流程圖 例4-41 試設(shè)計(jì)一個(gè)帶阻IIR數(shù)字濾波器，其具體的要求是：通帶的截止頻率：wp1650Hz、wp2850Hz；阻帶的截止頻率：ws1700Hz、ws2800Hz；通帶內(nèi)的最大衰減為rp0.1dB；阻帶內(nèi)的最小衰減為rs50dB；采樣頻率為Fs2000Hz。MATLAB源程序設(shè)計(jì)如下： wp1=650;wp2=850;ws1=7

28、00;ws2=800;rp=0.1;rs=50;Fs=2000; wp=wp1,wp2/(Fs/2);ws=ws1,ws2/(Fs/2); %利用Nyquist頻率頻率歸一化 N,wc=ellipord(wp,ws,rp,rs,z); %求濾波器階數(shù) num,den=ellip(N,rp,rs,wc,stop);%求濾波器傳遞函數(shù) H,W=freqz(num,den); %繪出頻率響應(yīng)曲線 plot(W*Fs/(2*pi),abs(H);grid; xlabel(頻率/Hz);ylabel(幅值)該程序運(yùn)行后的幅頻響應(yīng)曲線如圖4.42所示。 共四十一頁4.9 FIR數(shù)字(shz)濾波器設(shè)計(jì)格式

29、：w = boxcar(M) 功能：返回M點(diǎn)矩形窗序列。MATLAB信號(hào)處理(xn ho ch l)工具箱中的窗函數(shù)法設(shè)計(jì)FIR數(shù)字濾波器的專用命令fir1( )。格式：Bfir1(N,wc)功能：設(shè)計(jì)一個(gè)具有線性相位的N階(N點(diǎn))的低通FIR數(shù)字濾波器，返回的向量B為濾波器的系數(shù)(單位沖激響應(yīng)序列)，其長(zhǎng)度為N+1。4.9.1窗函數(shù)設(shè)計(jì)法窗函數(shù)設(shè)計(jì)的基本原理： h(n)=w(n)hd(n) w(n)為窗函數(shù)， hd(n)理想數(shù)字濾波器的單位沖激響應(yīng)。 在MATLAB信號(hào)處理工具箱中為用戶提供了Boxcar (矩形)、Bartlet(巴特利特)、Hanning(漢寧)等窗函數(shù)，這些窗函數(shù)的調(diào)用

30、格式相同。 FIR數(shù)字濾波器的單位沖激響應(yīng)h(n)滿足偶(奇)對(duì)稱 h(n)=h(N-n-1) 或 h(n)=-h(N-n-1)FIR數(shù)字濾波器具有線性相位：或共四十一頁例4-43 用矩形窗設(shè)計(jì)線性相位FIR低通濾波器。該濾波器的通帶截止頻率wc=pi/4，單位脈沖響h(n)的長(zhǎng)度M=21。并繪出h(n)及其幅度響應(yīng)特性(txng)曲線。MATLAB源程序?yàn)椋篗=21; wc=pi/4; % 理想低通濾波器參數(shù)n=0:M-1; r=(M-1)/2;nr=n-r+eps*(n-r)=0);hdn=sin(wc*nr)/pi./nr; % 計(jì)算理想低通單位脈沖響應(yīng)hd(n)if rem(M,2)=

31、0, hdn(r+1)=wc/pi; end% M為奇數(shù)時(shí)，處理n=r點(diǎn)的0/0型wn1=boxcar(M); % 矩形窗hn1=hdn.*wn1; % 加窗subplot(2,1,1);stem(n,hn1,.); line(0,20,0,0);xlabel(n),ylabel(h(n),title(矩形窗設(shè)計(jì)的h(n);hw1=fft(hn1,512); w1=2*0:511/512;%求頻譜subplot(2,1,2), plot(w1,20*log10(abs(hw1)xlabel(w/pi), ylabel(幅度(dB); title(幅度特性(dB);程序運(yùn)行結(jié)果如圖4.44所示。共四十一頁4.9.2頻率(pnl)抽樣法1. 頻率抽樣法的基本原理 對(duì)所期望的濾波器的