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

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

第 7章 在數(shù)字信號處理中的應(yīng)用,7.1 離散信號的產(chǎn)生及時域處理,時域離散信號用x(n)表示,時間變量n（表示采樣位置）只能取整數(shù)。因此,x(n)是一個離散序列,以后簡稱序列。用一個向量x不足以表示序列值x(n)。必須再用另一個等長的定位時間變量n。x和n同時使用才能完整地表示一個序列,由于n序列是按整數(shù)遞增的,可簡單地用其初值ns決定,因為它的終值nf取決于ns 和x的長度length(x),故可寫成： ? n = [ns:nf] 或 n = [ns: ns?length(x)?1],例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)元素的運算。,例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，也稱為加窗運算。序列的合成和截取就是相加和相乘。,,,,例7.3 序列的移位和周期延拓,已知 ，利用MATLAB生成并圖示 表示x(n)以8為周期的延拓）和 解：方法1，利用矩陣乘法和冒號運算 x=[1 2 3 4];y=x*ones(1，3); 方法2，采用求余函數(shù)mod， y = x(mod(n, M)+1)可實現(xiàn)對x(n) 以M為周期的周期延拓。加1是因為MATLAB向量下標(biāo)只能從1開始，,,,,例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即可。,例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)行判斷。,例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。,,,求z的逆變換的方法,對于z變換分式 可以用部分分式法或長除法求其反變換。 用函數(shù)residuez可以求出它的極點留數(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變換可寫成 兩個多項式乘積 可用conv函數(shù)來求得。n數(shù)組要自己判別。n的起點ns = ns1 + ns2 = ?3，終點nf = nf1 + nf2 = 2。 n=ns:nf。由x和n即可得出X(z)。,,,,例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變換，得,,,,,例7.8 z多項式分式逆變換(續(xù)),由程序算出nsy =-1，留數(shù)、極點分別為 r = -57.7581 和 204.7581 p = 0.7791 和 0.3209 k = -150 -30 代入 得,,,例7.9 離散時間傅里葉變換,取周期的正弦信號，作8點采樣，求它的連續(xù)頻譜。然后對該信號進(jìn)行N個周期延拓，再求它的連續(xù)頻譜。把N無限增大，比較分析其結(jié)果。 解： 先求離散傅立葉變換的MATLAB子程序 最后得到X = x*exp(-j*w*n‘)。 有了子程序，本例就沒有什么難度了。,,例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)系。,例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所示。,例7.10 采樣頻率與頻譜混疊(續(xù)),由于 由以上關(guān)系式可見，采樣信號的頻譜函數(shù)是原模擬信號頻譜函數(shù)的周期延拓，延拓周期為2?/T。如果以頻率f為自變量（? = 2?f），則以采樣頻率fs = 1/T為延拓周期。對頻帶限于fc的模擬信號xa(t)，只有當(dāng)fs≥2fc時，采樣后 才不會發(fā)生頻譜混疊失真。這就是著名的采樣定理,,,例7.11 由離散序列恢復(fù)模擬信號,用時域內(nèi)插公式 其中 模擬用理想低通濾波器恢復(fù)的過程，觀察恢復(fù)波形，計算出最大恢復(fù)誤差。 解：這個公式與卷積公式相像，可以用向量和矩陣乘法來解決。,,,例7.11 由離散序列恢復(fù)模擬信號,xa = x * g (TNM) = x * G 其中 G = sinc(Fs * TNM) M表示在兩個采樣點之間增加的間隔數(shù)，使輸出更密，更接近模擬信號。,,例7.12 梳狀濾波器零極點和幅特性,梳狀濾波器系統(tǒng)函數(shù)有如下兩種類型。 FIR型： IIR型： freqz 數(shù)字濾波器頻率特性計算和繪制函數(shù) zplane H(z)的零-極點圖繪制。 解：調(diào)用函數(shù)freqz和zplane 很容易寫出程序q712.m。,,,例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),例7.14 用符號運算工具箱解z變換,解：無限長度時間序列的z變換和逆z變換都屬于符號運算的范圍。MATLAB的symbolic（符號運算）工具箱已提供了這種函數(shù)。如果讀者已在計算機上安裝了這個工具箱，可以鍵入以下程序。 MATLAB程序q714.m 其特點是程序的開始要指定符號自變量 syms z n a N w0 % 規(guī)定z,n,a…為符號變量,7.3 離散傅里葉變換（DFT）,定義DFT： 用類似于例7.9中的方法，可把（7.3）式寫成矩陣乘法運算。 其中，xn為序列行向量，Wnk是一NN階方陣， 而 稱為旋轉(zhuǎn)因子。,,,,,7.3 離散傅里葉變換（DFT）,用矩陣乘法計算N點DFT的程序如下。 ? MATLAB程序q73a.m %用矩陣乘法計算N點DFT clear;close all xn=input(請輸入序列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; % 計算N點DFT,例7.15 序列的離散傅立葉變換,求復(fù)正弦序列 余弦序列 正弦序列 的離散傅立葉變換，分別按N =16和N =8進(jìn)行計算。繪出幅頻特性曲線， 進(jìn)行比較討論。,,,,例7.15 序列的離散傅立葉變換,在截取16點時，得到的是完整的余弦波形；而截取8點時，得到的是半截的余弦波形，當(dāng)然有大量的諧波成分。,例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)系。,,,例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)系。,,,,,,例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的情況，請讀者自己編程驗證。,,例7.18 快速卷積,快速卷積就是根據(jù)DFT的循環(huán)卷積性質(zhì)，將時域卷積轉(zhuǎn)換為頻域相乘，最后再進(jìn)行IDFT得到時域卷積序列y(n)。其中時域和頻域之間的變換均用FFT實現(xiàn)，所以使卷積速度大大提高?？驁D如下:,,例7.19 用DFT求連續(xù)信號頻譜,在計算機上用DFT對模擬信號進(jìn)行譜分析時，只能以有限大的采樣頻率fs對模擬信號采樣有限點樣本序列（等價于截取模擬信號一段進(jìn)行采樣）作DFT變換，得到模擬信號的近似頻譜。其誤差主要來自以下因素： ① 截斷效應(yīng)（頻譜泄露和譜間干擾） ② 頻譜混疊失真 因素①使譜分辨率（能分辨開的兩根譜線間的最小間距）降低，并產(chǎn)生譜間干擾； 因素②使折疊頻率（fs / 2）附近的頻譜產(chǎn)生較大失真。,例7.19 用DFT求連續(xù)信號頻譜,加大截取長度Tp可提高頻率分辨率；選擇合適的窗函數(shù)可降低譜間干擾；而頻譜混疊失真要通過提高采樣頻率fs和（或）預(yù)濾波（在采樣之前濾除折疊頻率以外的頻率成分）來改善。 編寫程序q719.m驗證截斷效應(yīng)及加窗的改善作用，先選取以下參數(shù)： ? 采樣頻率fs = 400Hz，T = 1/fs ? 采樣信號序列 ? 對x(n)作4096點DFT作為的近似頻譜Xa(jf)。,,例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窗的頻率分辨率高（泄漏?。?，但譜間干擾剛好相反。,例7.20 IIR濾波器直接型的轉(zhuǎn)換,程序調(diào)用了信號處理工具箱函數(shù)tf2sos和擴展函數(shù)dir2par， ? [sos, g] = tf2sos(B, A) 實現(xiàn)從直接型到級聯(lián)型（二階分割形式）的轉(zhuǎn)換。g為式中的增益，sos為L6階矩陣，表示式中的系數(shù)。,,例7.20 IIR濾波器直接型的轉(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之后。,例7.20 IIR濾波器直接型的轉(zhuǎn)換,根據(jù)計算結(jié)果, 級聯(lián)型H(z)表達(dá)式： 級聯(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ù)實現(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)歸一化）。,例7.21 直接型結(jié)構(gòu)到格型梯形結(jié)構(gòu),直接型系統(tǒng)函數(shù) 轉(zhuǎn)換為格型梯形結(jié)構(gòu),,,例7.22 FIR濾波器直接型到其他型,系統(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.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濾波器直接型到其他型,得出級聯(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ù)的計算機推導(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。,例7.24的結(jié)構(gòu)圖與方程,例7.24的方程的矩陣形式,由此可以求出系統(tǒng)函數(shù),,,例7.24解出的系統(tǒng)函數(shù),程序運行的結(jié)果 如果加入一個激勵x(n)?A?(n)，則得出,,,,7.5 FIR數(shù)字濾波器設(shè)計,濾波器的特性指標(biāo) 用絕對值δ1,δ2表示; 用分貝Rp,Rs表示,,,（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ù)。,,（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),例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)。,例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等。,,例7.25 窗函數(shù)法設(shè)計數(shù)字濾波器,對兩種窗函數(shù)的設(shè)計結(jié)果分別如右圖7.25-1和圖7.25-2所示。,,,工具箱設(shè)計函數(shù)fir1和fir2,MATLAB提供了基于窗函數(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),例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。,例7.27 用remez函數(shù)低通濾波器,解：先由題意計算設(shè)計參數(shù) f = [1/4,5/16]，m = [1,0]; dev的計算稍復(fù)雜一些，由于 所以 有了這幾個參數(shù)就可以調(diào)用remezord和remez函數(shù)了.,,,例7.27 用remez函數(shù)低通濾波器,橫線為-3dB，兩條豎線分別位于頻率?/4和5?/16。顯然，通帶指標(biāo)稍有富裕，過渡帶寬度和阻帶最小衰減剛好滿足指標(biāo)要求。,程序輸出的幅頻特性,例7.28 remez函數(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]。,例7.28 remez函數(shù)設(shè)計高通濾波器,程序運行結(jié)果如圖 由圖可見，w較大的頻段逼近精度較高；w較小的頻段逼近精度較低。N較大時逼近精度較高,N較小時逼近精度較低 。,,7.6 IIR數(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ù)的功能和用法。,IIR數(shù)字濾波器設(shè)計流程圖,模擬低通濾波器原型設(shè)計 Buttap,cheb1ap,cheb2ap besselap，ellipap函數(shù),頻率變換(變?yōu)楦咄?，帶通，帶阻? lp2lp，lp2hp，lp2bp，lp2bs,模擬數(shù)字變換 bilinear impinvar,,,,合為一步的設(shè)計函數(shù) butter，cheb1，cheb2，ellip， besself,求最小階數(shù)N Buttord, cheb1ord Cheb2ord,ellipord,,,,,,巴特沃斯濾波器設(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),集成的數(shù)字濾波器設(shè)計函數(shù),把（2）、（3）、（4）合為一步的數(shù)字濾波器設(shè)計函數(shù)butter(N, wc, ftype) ? [B, A] = butter (N, wc) 設(shè)計低通或帶通數(shù)字濾波器系數(shù)B，A（當(dāng)為帶通濾波器時，第(1)類函數(shù)由wp = [wp1, wp2]會自動生成wc = [w1, w2]）。 ? [B, A] = butter (N, wc, high) 設(shè)計高通數(shù)字濾波器系數(shù)B，A。 ? [B, A] = butter (N, wc, stop) 設(shè)計帶阻數(shù)字濾波器系數(shù)B，A。 butter(N, wc, ftype)還有零極增益和狀態(tài)空間形式，讀者可用help命令查閱。,例7.29 巴特沃斯模擬濾波器設(shè)計,設(shè)計一個低通巴特沃斯模擬濾波器，指標(biāo)如下。 通帶頻率：fp = 3400Hz，最大衰減：Rp = 3dB 阻帶頻率：fs = 4000Hz，最小衰減：As = 40dB 解:它的系統(tǒng)函數(shù)完全由階數(shù)N和3dB截止頻率?c決定。而N和?c是由濾波器設(shè)計指標(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ù) 單極點的N階模擬濾波器Ha(s)，用部分分式展開為 脈沖響應(yīng)不變法的數(shù)字化結(jié)果為 工具箱函數(shù)impinvar可實現(xiàn)以上計算，格式為 ? [Bz, Az] = impinvar(B, A, Fs),,,模擬轉(zhuǎn)換為數(shù)字：雙線性變換法,? 雙線性變換法函數(shù)bilinear 雙線性變換法的原理是用 代換Ha(s)中的s值，得到H(z)。bilinear函數(shù)用來實現(xiàn)這個轉(zhuǎn)換。其使用格式為 ? [Bz, Az] = bilinear(B, A, Fs) 脈沖響應(yīng)不變法的缺點是存在頻率混疊失真。雙線性變換法可完全消除頻率混疊失真，缺點是存在非線性頻率失真。,,例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）可見，對脈沖響應(yīng)不變法，采樣頻率Fs越高（T越?。殳B越??；由圖7.30（c）可見，對雙線性變換法，無頻率混疊，但存在非線性失真。,例7.31 切比雪夫Ⅱ數(shù)字濾波器設(shè)計,解：切比雪夫Ⅰ型濾波器通帶內(nèi)為等波紋，阻帶內(nèi)單調(diào)下降；切比雪夫Ⅱ型濾波器通帶內(nèi)為單調(diào)下降，阻帶內(nèi)等波紋。調(diào)用cheb2ord函數(shù)和cheby2函數(shù)使切比雪夫Ⅱ型設(shè)計變得非常簡單。 先用[N, wc] = Cheb2ord(wp, ws, Rp, Rs)求出N和 wc，提供函數(shù)cheby2的輸入變元，再由[B, A] = cheby2(N, Rp, wc)設(shè)計切比雪夫Ⅱ型數(shù)字濾波器。B和A分別為H(z)的分子和分母多項式系數(shù)。 對切比雪夫Ⅰ型濾波器，同樣有相應(yīng)的工具箱函數(shù)cheb1ord和cheby1。,例7.32 橢圓帶通數(shù)字濾波器設(shè)計,設(shè)計一個橢圓帶通數(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è)計,用雙線性變換法從Ha(s)設(shè)計四階帶通巴特沃斯數(shù)字濾波器，并圖示 （設(shè)采樣周期T=1s）。指標(biāo)如下： 解：本題主要涉及三個步驟： （1）由數(shù)字濾波器指標(biāo)求相應(yīng)的模擬濾波器指標(biāo)； （2）模擬濾波器頻率變換（因為已給定階數(shù)和模擬濾波器歸一化低通原型）； （3）由相應(yīng)的模擬濾波器到數(shù)字濾波器轉(zhuǎn)換（雙線性變換法）。,,例7.33 高通和帶通數(shù)字濾波器設(shè)計,首先要設(shè)計出與該指標(biāo)相應(yīng)的四階Butterworth模擬濾波器。然后，調(diào)用bilinear函數(shù)將其轉(zhuǎn)換成數(shù)字濾波器。對雙線性變換法，由數(shù)字邊界頻率求相應(yīng)的模擬邊界頻率時，一定要考慮預(yù)畸變矯正。這樣，最終設(shè)計結(jié)果才能滿足所給指標(biāo)。 步驟（1）： ? 設(shè)計高通濾波器時，模擬高通3dB截止頻率為 ? 設(shè)計帶通濾波器時，模擬帶通3dB截止頻率為,,例7.33 高通和帶通數(shù)字濾波器設(shè)計,步驟（2）： 可調(diào)用MATLAB頻率變換函數(shù)lp2lp，lp2hp，lp2bp，分別實現(xiàn)從模擬低通到模擬低通、高通、帶通、帶阻的頻率變換。本題用到的lp2hp和lp2bp的格式簡要說明如下： ? [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。,,,,,