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

《Matlab在數(shù)字信號(hào)處理中的應(yīng)用.ppt》由會(huì)員分享，可在線閱讀，更多相關(guān)《Matlab在數(shù)字信號(hào)處理中的應(yīng)用.ppt（73頁珍藏版）》請(qǐng)?jiān)谘b配圖網(wǎng)上搜索。

第 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è)等長的定位時(shí)間變量n。x和n同時(shí)使用才能完整地表示一個(gè)序列,由于n序列是按整數(shù)遞增的,可簡(jiǎn)單地用其初值ns決定,因?yàn)樗慕K值nf取決于ns 和x的長度length(x),故可寫成： ? n = [ns:nf] 或 n = [ns: ns?length(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。 解：編程的思路是把序列長度延拓到覆蓋n1和n2的范圍，這樣才能把兩序列的時(shí)間變量對(duì)應(yīng)起來，然后進(jìn)行對(duì)應(yīng)元素的運(yùn)算。,例7.2 序列的合成和截取,用例6.13的結(jié)果編寫產(chǎn)生矩形序列RN(n)的程序。序列起點(diǎn)為n0，矩形序列起點(diǎn)為n1，長度為N（n0，n1，N由鍵盤輸入）。并用它截取一個(gè)復(fù)正弦序列exp(jπn/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ì)算原理見例6.14，在這里用工具箱函數(shù)filter來解。如果已知系統(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的長度分別為Nx=20和Nh=10；結(jié)果向量y的長度為length(y)=Nx+Nh-1。,,,求z的逆變換的方法,對(duì)于z變換分式 可以用部分分式法或長除法求其反變換。 用函數(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ù)來求得。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無限增大，比較分析其結(jié)果。 解： 先求離散傅立葉變換的MATLAB子程序 最后得到X = x*exp(-j*w*n‘)。 有了子程序，本例就沒有什么難度了。,,例7.9 離散時(shí)間傅里葉變換2,? 程序運(yùn)行結(jié)果 執(zhí)行程序q709并按提示鍵入N = 4，所得圖形如圖7.10所示。N取得愈大，其峰值愈大，寬度愈窄。當(dāng)N取得很大時(shí)，會(huì)出現(xiàn)內(nèi)存不足的問題，這是用矩陣乘法做傅里葉變換的缺點(diǎn)。另外，因?yàn)槟菚r(shí)峰值點(diǎn)處的寬度很窄，也會(huì)出現(xiàn)所選頻點(diǎn)對(duì)不上峰值點(diǎn)的問題。所以對(duì)于N無限增大的情況，必須用fft函數(shù)來求。這時(shí)用連續(xù)頻譜也沒有意義了。這里用同樣的橫坐標(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)系式可見，采樣信號(hào)的頻譜函數(shù)是原模擬信號(hào)頻譜函數(shù)的周期延拓，延拓周期為2?/T。如果以頻率f為自變量（? = 2?f），則以采樣頻率fs = 1/T為延拓周期。對(duì)頻帶限于fc的模擬信號(hào)xa(t)，只有當(dāng)fs≥2fc時(shí)，采樣后 才不會(huì)發(fā)生頻譜混疊失真。這就是著名的采樣定理,,,例7.11 由離散序列恢復(fù)模擬信號(hào),用時(shí)域內(nèi)插公式 其中 模擬用理想低通濾波器恢復(fù)的過程，觀察恢復(fù)波形，計(jì)算出最大恢復(fù)誤差。 解：這個(gè)公式與卷積公式相像，可以用向量和矩陣乘法來解決。,,,例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.04?n) + cos(0.08?n) + cos(0.4?n) + 0.3?(n)，0≤n≤63 通過低通濾波器，計(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變換,解：無限長度時(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來驗(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) = 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)系。,,,,,,例7.17 頻域與時(shí)域采樣對(duì)偶性,由于頻域在[0，2?]上的采樣點(diǎn)數(shù)N（N = 32）小于x(n)的長度M（M = 40），所以，產(chǎn)生時(shí)域混疊現(xiàn)象，不能由X1(k)恢復(fù)原序列x(n)。只有滿足N≥M時(shí)，可由頻域采樣X1(k)得到原序列x(n)。 這就是頻域采樣定理。對(duì)N≥M的情況，請(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)，所以使卷積速度大大提高。框圖如下:,,例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)的近似頻譜。其誤差主要來自以下因素： ① 截?cái)嘈?yīng)（頻譜泄露和譜間干擾） ② 頻譜混疊失真 因素①使譜分辨率（能分辨開的兩根譜線間的最小間距）降低，并產(chǎn)生譜間干擾； 因素②使折疊頻率（fs / 2）附近的頻譜產(chǎn)生較大失真。,例7.19 用DFT求連續(xù)信號(hào)頻譜,加大截取長度Tp可提高頻率分辨率；選擇合適的窗函數(shù)可降低譜間干擾；而頻譜混疊失真要通過提高采樣頻率fs和（或）預(yù)濾波（在采樣之前濾除折疊頻率以外的頻率成分）來改善。 編寫程序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é)果。由圖可見，由于截?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為L6階矩陣，表示式中的系數(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】中介紹的方法和公式同樣可以用來求離散域的數(shù)字濾波器的系統(tǒng)函數(shù)。不同的地方僅僅在于節(jié)點(diǎn)方程中出現(xiàn)了作為系數(shù)的符號(hào)變量z?1，它將出現(xiàn)在系數(shù)矩陣中。MATLAB是不能處理上標(biāo)變量的，因此在程序中設(shè)q=z ?1，在計(jì)算完成后再人工地把結(jié)果中的q恢復(fù)為z?1。,例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)來截取hd(n)的適當(dā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)的長度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，其窗口長度M由過渡帶寬度B = 0.15? 決定，Blackman窗設(shè)計(jì)的濾波器過渡帶寬度為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)稍有富裕，過渡帶寬度和阻帶最小衰減剛好滿足指標(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)可以互換，過渡帶寬度與通帶波紋和阻帶衰減指標(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é)果如圖 由圖可見，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ì)的全過程的各個(gè)步驟，MATLAB都提供了相應(yīng)的工具箱函數(shù)，使IIR數(shù)字濾波器設(shè)計(jì)變得非常簡(jiǎn)單。本節(jié)主要結(jié)合例題介紹這些IIR濾波器設(shè)計(jì)的工具箱函數(shù)。 IIR數(shù)字濾波器的設(shè)計(jì)步驟由以下的流程圖來表示。下面以巴特沃斯濾波器設(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ù)以上公式編寫的。因此就無需再記憶這些公式了。,,模擬轉(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ù)用來實(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)曲線。 程序中的核心語句是以下兩條： [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）可見，對(duì)脈沖響應(yīng)不變法，采樣頻率Fs越高（T越?。殳B越?。挥蓤D7.30（c）可見，對(duì)雙線性變換法，無頻率混疊，但存在非線性失真。,例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)要說明如下： ? [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。,,,,,