Matlab在數(shù)字信號(hào)處理中的應(yīng)用.ppt

第 7章 在數(shù)字信號(hào)處理中的應(yīng)用,7.1 離散信號(hào)的產(chǎn)生及時(shí)域處理,時(shí)域離散信號(hào)用x(n)表示,時(shí)間變量n（表示采樣位置）只能取整數(shù)。因此,x(n)是一個(gè)離散序列,以后簡(jiǎn)稱序列。用一個(gè)向量x不足以表示序列值x(n)。必須再用另一個(gè)等長(zhǎng)的定位時(shí)間變量n。x和n同時(shí)使用才能完整地表示一個(gè)序列,由于n序列是按整數(shù)遞增的,可簡(jiǎn)單地用其初值ns決定,因?yàn)樗慕K值nf取決于ns 和x的長(zhǎng)度length(x),故可寫成： n = ns:nf 或 n = ns: nslength(x)1,例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)算。,例7.2 序列的合成和截取,用例6.13的結(jié)果編寫產(chǎn)生矩形序列RN(n)的程序。序列起點(diǎn)為n0，矩形序列起點(diǎn)為n1，長(zhǎng)度為N（n0，n1，N由鍵盤輸入）。并用它截取一個(gè)復(fù)正弦序列exp(jn/8) . 解：建模:矩形序列可看成兩個(gè)階躍序列之差。 用MATLAB邏輯關(guān)系產(chǎn)生矩形序列x2(n)。而用它截取任何序列相當(dāng)于元素群相乘x2.*x，也稱為加窗運(yùn)算。序列的合成和截取就是相加和相乘。,例7.3 序列的移位和周期延拓,已知 ，利用MATLAB生成并圖示 表示x(n)以8為周期的延拓）和 解：方法1，利用矩陣乘法和冒號(hào)運(yùn)算 x=1 2 3 4;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開始，,例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即可。,例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)行判斷。,例7.6 離散序列的卷積計(jì)算,給出兩個(gè)序列 和 ,計(jì)算其卷積y(n)，并圖示各輸入輸出序列。 解：在例6.4中，已經(jīng)給出了直接調(diào)用MATLAB的卷積函數(shù)conv的方法，也給出了自編卷積計(jì)算程序的方法，要注意的是本例時(shí)間變量的設(shè)定和移位方法。在本例中，設(shè)定n為從零開始，向量x和h的長(zhǎng)度分別為Nx=20和Nh=10；結(jié)果向量y的長(zhǎng)度為length(y)=Nx+Nh-1。,求z的逆變換的方法,對(duì)于z變換分式 可以用部分分式法或長(zhǎng)除法求其反變換。 用函數(shù)residuez可以求出它的極點(diǎn)留數(shù)分解 其中 r, p, k = residuez (B, A) 其反變換為：,例7.7 有限序列的z和逆z變換,兩序列x1 = 1,2,3， n1 = -1:1 及x2 = 2,4,3,5，n2 = -2:1，求出x1與x2及其卷積x的z變換。 解：其z變換可寫成 兩個(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)。,例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變換，得,例7.8 z多項(xiàng)式分式逆變換(續(xù)),由程序算出nsy =-1，留數(shù)、極點(diǎn)分別為 r = -57.7581 和 204.7581 p = 0.7791 和 0.3209 k = -150 -30 代入 得,例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)有什么難度了。,例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)系。,例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)，畫出其幅度頻譜。采樣時(shí)間區(qū)間均為0.1秒。為了便于比較，畫出了幅度歸一化的幅頻曲線，如圖7.11所示。,例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)fs2fc時(shí)，采樣后 才不會(huì)發(fā)生頻譜混疊失真。這就是著名的采樣定理,例7.11 由離散序列恢復(fù)模擬信號(hào),用時(shí)域內(nèi)插公式 其中 模擬用理想低通濾波器恢復(fù)的過(guò)程，觀察恢復(fù)波形，計(jì)算出最大恢復(fù)誤差。 解：這個(gè)公式與卷積公式相像，可以用向量和矩陣乘法來(lái)解決。,例7.11 由離散序列恢復(fù)模擬信號(hào),xa = x * g (TNM) = x * G 其中 G = sinc(Fs * TNM) M表示在兩個(gè)采樣點(diǎn)之間增加的間隔數(shù)，使輸出更密，更接近模擬信號(hào)。,例7.12 梳狀濾波器零極點(diǎn)和幅特性,梳狀濾波器系統(tǒng)函數(shù)有如下兩種類型。 FIR型： IIR型： freqz 數(shù)字濾波器頻率特性計(jì)算和繪制函數(shù) zplane H(z)的零-極點(diǎn)圖繪制。 解：調(diào)用函數(shù)freqz和zplane 很容易寫出程序q712.m。,例7.13 低通濾波及時(shí)域卷積定理,輸入信號(hào) x(n) = cos(0.04n) + cos(0.08n) + cos(0.4n) + 0.3(n)，0n63 通過(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),例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)是程序的開始要指定符號(hào)自變量 syms z n a N w0 % 規(guī)定z,n,a為符號(hào)變量,7.3 離散傅里葉變換（DFT）,定義DFT： 用類似于例7.9中的方法，可把（7.3）式寫成矩陣乘法運(yùn)算。 其中，xn為序列行向量，Wnk是一NN階方陣， 而 稱為旋轉(zhuǎn)因子。,7.3 離散傅里葉變換（DFT）,用矩陣乘法計(jì)算N點(diǎn)DFT的程序如下。 MATLAB程序q73a.m %用矩陣乘法計(jì)算N點(diǎn)DFT clear;close all xn=input(請(qǐng)輸入序列x= ); N = length(xn); % n=0:N-1; k=n; nk=n*k; % 生成NN方陣 WN=exp(-j*2*pi/N); Wnk=WN.nk; % 產(chǎn)生旋轉(zhuǎn)因子矩陣 Xk=xn*Wnk; % 計(jì)算N點(diǎn)DFT,例7.15 序列的離散傅立葉變換,求復(fù)正弦序列 余弦序列 正弦序列 的離散傅立葉變換，分別按N =16和N =8進(jìn)行計(jì)算。繪出幅頻特性曲線， 進(jìn)行比較討論。,例7.15 序列的離散傅立葉變換,在截取16點(diǎn)時(shí)，得到的是完整的余弦波形；而截取8點(diǎn)時(shí)，得到的是半截的余弦波形，當(dāng)然有大量的諧波成分。,例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)系。,例7.17 頻域與時(shí)域采樣對(duì)偶性,（1）產(chǎn)生三角波序列 （2）對(duì)M = 40，計(jì)算x(n)的64點(diǎn)DFT，并圖示x(n)和X(k) = DFTx(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)系。,例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)。只有滿足NM時(shí)，可由頻域采樣X(jué)1(k)得到原序列x(n)。 這就是頻域采樣定理。對(duì)NM的情況，請(qǐng)讀者自己編程驗(yàn)證。,例7.18 快速卷積,快速卷積就是根據(jù)DFT的循環(huán)卷積性質(zhì)，將時(shí)域卷積轉(zhuǎn)換為頻域相乘，最后再進(jìn)行IDFT得到時(shí)域卷積序列y(n)。其中時(shí)域和頻域之間的變換均用FFT實(shí)現(xiàn)，所以使卷積速度大大提高?？驁D如下:,例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)（頻譜泄露和譜間干擾） 頻譜混疊失真 因素使譜分辨率（能分辨開的兩根譜線間的最小間距）降低，并產(chǎn)生譜間干擾； 因素使折疊頻率（fs / 2）附近的頻譜產(chǎn)生較大失真。,例7.19 用DFT求連續(xù)信號(hào)頻譜,加大截取長(zhǎng)度Tp可提高頻率分辨率；選擇合適的窗函數(shù)可降低譜間干擾；而頻譜混疊失真要通過(guò)提高采樣頻率fs和（或）預(yù)濾波（在采樣之前濾除折疊頻率以外的頻率成分）來(lái)改善。 編寫程序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)。,例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間干擾剛好相反。,例7.20 IIR濾波器直接型的轉(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ù)。,例7.20 IIR濾波器直接型的轉(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中開發(fā)的，不屬于MATLAB工具箱函數(shù)，所以將其M文件清單附在程序q720.m之后。,例7.20 IIR濾波器直接型的轉(zhuǎn)換,根據(jù)計(jì)算結(jié)果, 級(jí)聯(lián)型H(z)表達(dá)式： 級(jí)聯(lián)型結(jié)構(gòu)圖。,例7.20 IIR濾波器直接型的轉(zhuǎn)換,并聯(lián)型結(jié)構(gòu)H(z)表達(dá)式 并聯(lián)型結(jié)構(gòu)圖,例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)歸一化）。,例7.21 直接型結(jié)構(gòu)到格型梯形結(jié)構(gòu),直接型系統(tǒng)函數(shù) 轉(zhuǎn)換為格型梯形結(jié)構(gòu),例7.22 FIR濾波器直接型到其他型,系統(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.0000 0.5360 0 1.0000 0 0 1.0000 0.0057 0.6219 1.0000 0 0 g = 2 格型結(jié)構(gòu)系數(shù)（反射系數(shù)）： K = 0.2500 0.5000 0.3333,例7.22 FIR濾波器直接型到其他型,得出級(jí)聯(lián)型為 格型結(jié)構(gòu)為,例7.23 FIR格型到直接型轉(zhuǎn)換,給定K = 2,1/4,1/2,1/3 ,用函數(shù) latc2tf即可 由B=latc2tf(K)得到B,寫出直接型結(jié)構(gòu),例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=z 1，在計(jì)算完成后再人工地把結(jié)果中的q恢復(fù)為z1。,例7.24的結(jié)構(gòu)圖與方程,例7.24的方程的矩陣形式,由此可以求出系統(tǒng)函數(shù),例7.24解出的系統(tǒng)函數(shù),程序運(yùn)行的結(jié)果 如果加入一個(gè)激勵(lì)x(n)A(n)，則得出,7.5 FIR數(shù)字濾波器設(shè)計(jì),濾波器的特性指標(biāo) 用絕對(duì)值1,2表示; 用分貝Rp,Rs表示,（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ù)。,（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),例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)。,例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等。,例7.25 窗函數(shù)法設(shè)計(jì)數(shù)字濾波器,對(duì)兩種窗函數(shù)的設(shè)計(jì)結(jié)果分別如右圖7.25-1和圖7.25-2所示。,工具箱設(shè)計(jì)函數(shù)fir1和fir2,MATLAB提供了基于窗函數(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),例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。,例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ù)了.,例7.27 用remez函數(shù)低通濾波器,橫線為-3dB，兩條豎線分別位于頻率/4和5/16。顯然，通帶指標(biāo)稍有富裕，過(guò)渡帶寬度和阻帶最小衰減剛好滿足指標(biāo)要求。,程序輸出的幅頻特性,例7.28 remez函數(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。,例7.28 remez函數(shù)設(shè)計(jì)高通濾波器,程序運(yùn)行結(jié)果如圖 由圖可見(jiàn)，w較大的頻段逼近精度較高；w較小的頻段逼近精度較低。N較大時(shí)逼近精度較高,N較小時(shí)逼近精度較低 。,7.6 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)單。本節(jié)主要結(jié)合例題介紹這些IIR濾波器設(shè)計(jì)的工具箱函數(shù)。 IIR數(shù)字濾波器的設(shè)計(jì)步驟由以下的流程圖來(lái)表示。下面以巴特沃斯濾波器設(shè)計(jì)函數(shù)為典型，介紹此流程圖中函數(shù)的功能和用法。,IIR數(shù)字濾波器設(shè)計(jì)流程圖,模擬低通濾波器原型設(shè)計(jì) Buttap,cheb1ap,cheb2ap besselap，ellipap函數(shù),頻率變換(變?yōu)楦咄?，帶通，帶阻? lp2lp，lp2hp，lp2bp，lp2bs,模擬數(shù)字變換 bilinear impinvar,合為一步的設(shè)計(jì)函數(shù) butter，cheb1，cheb2，ellip， besself,求最小階數(shù)N Buttord, cheb1ord Cheb2ord,ellipord,巴特沃斯濾波器設(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),集成的數(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命令查閱。,例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ù)以上公式編寫的。因此就無(wú)需再記憶這些公式了。,模擬轉(zhuǎn)換為數(shù)字：脈沖響應(yīng)不變法,模擬濾波器離散化的基本方法有脈沖響應(yīng)不變法和雙線性變換法。 脈沖響應(yīng)不變法及impinvar函數(shù) 單極點(diǎn)的N階模擬濾波器Ha(s)，用部分分式展開為 脈沖響應(yīng)不變法的數(shù)字化結(jié)果為 工具箱函數(shù)impinvar可實(shí)現(xiàn)以上計(jì)算，格式為 Bz, Az = impinvar(B, A, Fs),模擬轉(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)是存在非線性頻率失真。,例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ù)離散化,例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ú)頻率混疊，但存在非線性失真。,例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ù)cheb1ord和cheby1。,例7.32 橢圓帶通數(shù)字濾波器設(shè)計(jì),設(shè)計(jì)一個(gè)橢圓帶通數(shù)字濾波器，要求如下 Rp = 1 dB， Rs = 60 dB 解：調(diào)用ellipord函數(shù)和ellip函數(shù)，代入?yún)?shù) wp = 0.25，0.45；ws = 0.15，0.55； Rp = 1；Rs = 60; 即可得到本題的程序 。,例7.33 高通和帶通數(shù)字濾波器設(shè)計(jì),用雙線性變換法從Ha(s)設(shè)計(jì)四階帶通巴特沃斯數(shù)字濾波器，并圖示 （設(shè)采樣周期T=1s）。指標(biāo)如下： 解：本題主要涉及三個(gè)步驟： （1）由數(shù)字濾波器指標(biāo)求相應(yīng)的模擬濾波器指標(biāo)； （2）模擬濾波器頻率變換（因?yàn)橐呀o定階數(shù)和模擬濾波器歸一化低通原型）； （3）由相應(yīng)的模擬濾波器到數(shù)字濾波器轉(zhuǎn)換（雙線性變換法）。,例7.33 高通和帶通數(shù)字濾波器設(shè)計(jì),首先要設(shè)計(jì)出與該指標(biāo)相應(yīng)的四階Butterworth模擬濾波器。然后，調(diào)用bilinear函數(shù)將其轉(zhuǎn)換成數(shù)字濾波器。對(duì)雙線性變換法，由數(shù)字邊界頻率求相應(yīng)的模擬邊界頻率時(shí)，一定要考慮預(yù)畸變矯正。這樣，最終設(shè)計(jì)結(jié)果才能滿足所給指標(biāo)。 步驟（1）： 設(shè)計(jì)高通濾波器時(shí)，模擬高通3dB截止頻率為 設(shè)計(jì)帶通濾波器時(shí)，模擬帶通3dB截止頻率為,例7.33 高通和帶通數(shù)字濾波器設(shè)計(jì),步驟（2）： 可調(diào)用MATLAB頻率變換函數(shù)lp2lp，lp2hp，lp2bp，分別實(shí)現(xiàn)從模擬低通到模擬低通、高通、帶通、帶阻的頻率變換。本題用到的lp2hp和lp2bp的格式簡(jiǎn)要說(shuō)明如下： Bt，At = lp2hp(B, A, wc) 將系數(shù)向量為B和A的模擬濾波器歸一化低通原型（3dB截止頻率為1rad/s）變換成3dB截止頻率為wc的高通模擬濾波器，返回高通模擬濾波器系數(shù)向量Bt和At。 步驟（3）： 調(diào)用bilinear函數(shù)將濾波器系數(shù)向量Bt和At轉(zhuǎn)換到數(shù)字濾波器Bz和Az。,