1、實驗三 用雙線性變換法設計IIR數字濾波器實驗項目名稱：用雙線性變換法設計IIR數字濾波器實驗項目性質：驗證性實驗所屬課程名稱：數字信號處理實驗計劃學時：2一. 實驗目的(1) 熟悉用雙線性變換法設計IIR數字濾波器的原理與方法。(2) 掌握數字濾波器的計算機仿真方法。(3) 通過觀察對實際心電圖信號的濾波作用，獲得數字濾波的感性知識。二. 實驗內容和要求(1) 用雙線性變換法設計一個巴特沃斯低通IIR數字濾波器。設計指標參數為：在通帶內頻率低于0.2時，最大衰減小于1dB；在阻帶內0.3，頻率區間上，最小衰減大與15dB。(2) 以0.02為采樣間隔，打印出數字濾波器在頻率區間0，/2的幅頻

2、響應特性曲線。(3) 用所設計的濾波器對實際心電圖信號采樣序列（在本實驗后面給出）進行仿真濾波處理，并打印出濾波前后的心電圖信號波形圖，觀察總結濾波作用與效果。三. 實驗主要儀器設備和材料計算機，MATLAB6.5或以上版本四. 實驗方法、步驟及結果測試(1) 復習有關巴特沃斯模擬濾波器設計和用雙線性變換法設計IIR數字濾波器的內容，用雙線性變換法設計數字濾波器系統函數。其中滿足本實驗要求的數字濾波器系統函數為： （3.1）式中： （3.2）根據設計指標，調用MATLAB信號處理工具箱buttord和butter，也可以得到。圖3-1 濾波器的組成由公式（3.1）和（3.2）可見，濾波器由三個

3、二階濾波器、和級聯而成，如圖3-1所示。(2) 編寫濾波器仿真程序，計算對心電圖采樣序列x(n)的響應序列y(n)。設yk(n)為第k級二階濾波器Hk(z)的輸出序列，yk-1(n)為輸入序列，如圖3-1所示。由（3.2）式可得到差分方程為：（3.3）當k=1時，yk-1(n)=x(n)。所以H(z)對x(n)的總響應序列y(n)可以用順序迭代算法得到。即依次對k=1，2，3，求解差分方程（3.3），最后得到y3(n)=y(n)。仿真程序就是實現上述求解差分方程和順序迭代算法的通用程序。也可以直接調用MATLAB filter函數實現仿真。(3) 在通用計算機上運行仿真濾波程序，并調用通用繪圖

4、子程序，完成實驗內容（2）和（3）。(4) 本實驗中用到的心電圖信號采用序列x(n)人體心電圖信號在測量過程往往受到工業高頻干擾，所以必須經過低通濾波處理后，才能作為判斷心臟功能的有用信息。下面給出一實際心電圖信號采樣序列樣本x(n)，其中存在高頻干擾。在實驗中，以x(n)作為輸入序列，濾除其中的干擾成分。x(n)=-4， -2， 0， -4， -6， -4， -2， -4， -6， -6， -4， -4， -6， -6， -2， 6， 12， 8， 0，-16， -38，-60，-84，-90，-66，-32， -4， -2， -4， 8， 12， 12， 10， 6， 6， 6， 4， 0

5、， 0， 0， 0， 0， -2， -4， 0， 0， 0， -2， -2， 0，0， -2， -2， -2， -2， 0MATLAB程序清單：%實驗三，用雙線性變換法設計IIR數字濾波器x=-4,-2,0,-4,-6,-4,-2,-4,-6,-6,-4,-4,-6,-6,-2,6,12,8,0,-16,-38,-60,-84,-90,-66,-32,-4,-2,-4,8,12,12,10,6,6,6,4,0,0,0,0,0,-2,-4,0,0,0,-2,-2,0,0,-2,-2,-2,-2,0;k=1;close all;figure(1);subplot(2,2,1);n=0:55; %更

6、正stem(n,x,'.');axis(0,56,-100,50); %更正hold on;n=0:60;m=zeros(61);plot(n,m);xlabel('n');ylabel('x(n)');title('心電圖信號采集序列x(n)');B=0.09036,2*0.09036,0.09036;A=1.2686,-0.7051;A1=1.0106,-0.3583;A2=0.9044,-0.2155;while(k<=3) y=filter(B,A,x); %The function is to filte the

7、singal x x=y; if k=2 A=A1; end if k=3 A=A2; end k=k+1;endsubplot(2,2,3)n=0:55; %更正stem(n,y,'.');axis(0,56,-15,5);hold on;n=0:60;m=zeros(61);plot(n,m);xlabel('n');ylabel('y(n)');title('三級濾波后的心電圖信號');%求數字濾波器的幅頻特性A=0.09036,0.1872,0.09036;B1=1,-1.2686,0.7051;B2=1,-1.0106,

8、0.3583;B3=1,-0.9044,0.2155;H1,w=freqz(A,B1,100);H2,w=freqz(A,B2,100);H3,w=freqz(A,B3,100);H4=H1.*(H2);H=H4.*(H3);mag=abs(H);db=20*log10(mag+eps)/max(mag);subplot(2,2,2)plot(w/pi,db);axis(0,0.5,-50,10);grid on; %更正title('濾波器的幅頻響應曲線');程序運行結果：直接運行程序，結果輸出濾波器幅頻特性曲線圖，有噪聲的心電圖采集信號波形圖和經過三級二階濾波器濾波后的心電

9、圖信號波形圖，見圖3-2，對比分析就可以看出低通濾波器濾除信號中高頻噪聲的濾波效果。圖3-2 心電圖信號波形圖五. 實驗報告要求(1) 簡述實驗原理及目的。(2) 由所打印的| H(ejw) |特性曲線及設計過程簡述雙線性變換法的特點。(3) 對比濾波前后的心電圖信號波形，說明數字濾波器的濾波過程與濾波作用。(4) 簡要回答思考題。六. 思考題(1) 用雙線性變換法設計數字濾波器過程中，變換公式中T的取值對設計結果有無影響？為什么？實驗四 用窗函數法設計FIR數字濾波器實驗項目名稱：用窗函數法設計FIR數字濾波器實驗項目性質：驗證性實驗所屬課程名稱：數字信號處理實驗計劃學時：2一. 實驗目的(

10、1) 掌握用窗函數法設計FIR數字濾波器的原理與方法。(2) 熟悉線性相位FIR數字濾波器的特性。(3) 了解各種窗函數對濾波特性的影響。二. 實驗內容和要求(1) 復習用窗函數法設計FIR數字濾波器一節內容，閱讀本實驗原理，掌握設計步驟。(2) 用升余弦窗設計一線性相位低通FIR數字濾波器，截止頻率。窗口長度N =15，33。要求在兩種窗口長度情況下，分別求出，打印出相應的幅頻特性和相頻特性曲線，觀察3dB帶寬和20dB帶寬?？偨Y窗口長度N 對濾波器特性的影響。設計低通FIR數字濾波器時，一般以理想低通濾波特性為逼近函數，即其中(3) ，用四種窗函數設計線性相位低通濾波器，繪制相應的幅頻特性

11、曲線，觀察3dB帶寬和20dB帶寬以及阻帶最小衰減，比較四種窗函數對濾波器特性的影響。三. 實驗主要儀器設備和材料計算機，MATLAB6.5或以上版本四. 實驗方法、步驟及結果測試如果所希望的濾波器的理想的頻率響應函數為，則其對應的單位脈沖響應為 （4.1）窗函數設計法的基本原理是用有限長單位脈沖響應序列逼近。由于往往是無限長序列，而且是非因果的，所以用窗函數將截斷，并進行加權處理，得到：（4.2）就作為實際設計的FIR數字濾波器的單位脈沖響應序列，其頻率響應函數為（4.3）式中，N為所選窗函數的長度。我們知道，用窗函數法設計的濾波器性能取決于窗函數的類型及窗口長度N的取值。設計過程中，要根據

12、對阻帶最小衰減和過渡帶寬度的要求選擇合適的窗函數類型和窗口長度N 。各種類型的窗函數可達到的阻帶最小衰減和過渡帶寬度見表4.1。表4.1 各種窗函數的基本參數窗函數旁瓣峰值幅度/dB過渡帶寬阻帶最小衰減/dB矩形窗-134/N-12三角形窗-258/N-25漢寧窗-318/N-44哈明窗-418/N-53不萊克曼窗-5712/N-74凱塞窗(=7.865)-5710/N-80這樣選定窗函數類型和長度N之后，求出單位脈沖響應，并按照式（4.3）求出。是否滿足要求，要進行演算。一般在尾部加零使長度滿足2的整數次冪，以便用FFT計算。如果要觀察細節，補零點數增多即可。如果不滿足要求，則要重新選擇窗函

13、數類型和長度N ，再次驗算，直至滿足要求。如果要求線性相位特性，則還必須滿足根據上式中的正、負號和長度N的奇偶性又將線性相位FIR濾波器分成四類。要根據所設計的濾波特性正確選擇其中一類，例如，要設計線性相位低通特性，可以選擇這一類，而不能選擇這一類。主程序框圖如圖4.1所示。其中幅度特性要求用dB表示。開始讀入窗口長度N計算hd(n)調用窗函數子程序求w(n)調用子程序（函數）計算H(k)=DFTh(n)調用繪圖子程序（函數）繪制H(k)幅度相位曲線結束圖4-1 主程序框圖計算h(n)= hd(n) w(n)設畫圖時，用打印幅度特性。第k點對應的頻率。為使曲線包絡更接近的幅度特性曲線，DFT變

14、換區間要選大些。例如窗口長度N=33時，可通過在末尾補零的方法，使長度變為64，再進行64點DFT，則可以得到更精確的幅度衰減特性曲線。下面給出MATLAB主程序：%實驗四，用窗函數法設計FIR數字濾波器b=1;close all;i=0;while(b); temp=menu('選擇窗函數長度N','N=10','N=15','N=20','N=25','N=30','N=33','N=35','N=40','N=45','N

15、=50','N=55','N=60','N=64'); menu1=10,15,20,25,30,33,35,40,45,50,55,60,64; N=menu1(temp); temp=menu('選擇逼近理想低通濾波器截止頻率Wc','Wc=pi/4','Wc=pi/2','Wc=3*pi/4','Wc=pi','Wc=0.5','Wc=1.0','Wc=1.5','Wc=2.0','

16、Wc=2.5','Wc=3.0'); menu2=pi/4,pi/2,3*pi/4,pi,0.5,1,1.5,2,2.5,3; w=menu2(temp); n=0:(N-1); hd=ideal(w,N); %得到理想低通濾波器 k=menu('請選擇窗口類型:','boxcar','hamming','hanning','blackman'); if k=1 B=boxcar(N); string='Boxcar','N=',num2str(N); els

17、e if k=2 B=hamming(N); string='Hamming','N=',num2str(N); else if k=3 B=hanning(N); string='Hanning','N=',num2str(N); else if k=4 B=blackman(N); string='Blackman','N=',num2str(N); end end end end h=hd.*(B)' %得到FIR數字濾波器 H,m=freqz(h,1,1024,'whole&

18、#39;); %求其頻率響應 mag=abs(H); %得到幅值 db=20*log10(mag+eps)/max(mag); pha=angle(H); %得到相位 i=i+1; figure(i) subplot(2,2,1); n=0:N-1; stem(n,h,'.'); axis(0,N-1,-0.1,0.3); hold on; n=0:N-1; x=zeros(N); plot(n,x,'-'); xlabel('n'); ylabel('h(n)'); title('實際低通濾波器的h(n)');

19、text(0.3*N),0.27,string); hold off; subplot(2,2,2); plot(m/pi,db); axis(0,1,-100,0); xlabel('w/pi'); ylabel('dB'); title('衰減特性(dB)'); grid; subplot(2,2,3); plot(m,pha); hold on; n=0:7; x=zeros(8); plot(n,x,'-'); title('相頻特性'); xlabel('頻率(rad)'); ylabel

20、('相位(rad)'); axis(0,3.15,-4,4); subplot(2,2,4); plot(m,mag); title('頻率特性'); xlabel('頻率W(rad)'); ylabel('幅值'); axis(0,3.15,0,1.5); text(0.9,1.2,string); b=menu('Do You want To Continue ?','Yes','No'); if b=2 b=0; endendtemp=menu('Close All F

21、igure ?','Yes','No');if temp=1 close allend%實驗四，用窗函數法設計FIR數字濾波器%實驗四中的子函數：產生理想低通濾波器單位脈沖響應hd(n)新建.m文件，把下列綠色區域語句復制，保存為ideal.mfunction hd=ideal(w,N);alpha=(N-1)/2;n=0:(N-1);m=n-alpha+eps;hd=sin(w*m)./(pi*m);新建.m文件，把下列紅色區域語句復制，保存為main.m注意：ideal.m 與main.m 一定要在同一文件夾內。close all;i=0;N=10;

22、 w=pi/4; n=0:(N-1); hd=ideal(w,N); %得到理想低通濾波器 B=boxcar(N); h=hd.*(B)' %得到FIR數字濾波器 H,m=freqz(h,1,1024,'whole'); %求其頻率響應 mag=abs(H); %得到幅值 db=20*log10(mag+eps)/max(mag); pha=angle(H); %得到相位 i=i+1; figure(i) subplot(2,2,1); n=0:N-1; stem(n,h,'.'); axis(0,N-1,-0.1,0.3); hold on; n=0:

23、N-1; x=zeros(N); plot(n,x,'-'); xlabel('n'); ylabel('h(n)'); title('實際低通濾波器的h(n)'); text(0.3*N),0.27,'boxcar,N=10'); hold off; subplot(2,2,2); plot(m/pi,db); axis(0,1,-100,0); xlabel('w/pi'); ylabel('dB'); title('衰減特性(dB)'); grid; subpl

24、ot(2,2,3); plot(m,pha); hold on; n=0:7; x=zeros(8); plot(n,x,'-'); title('相頻特性'); xlabel('頻率(rad)'); ylabel('相位(rad)'); axis(0,3.15,-4,4); subplot(2,2,4); plot(m,mag); title('頻率特性'); xlabel('頻率W(rad)'); ylabel('幅值'); axis(0,3.15,0,1.5); text(0.

25、9,1.2,'boxcar,N=10');程序運行結果：運行程序，根據實驗內容要求和程序提示選擇你要進行的實驗參數。三個實驗參數選定后，程序運行輸出用所選窗函數設計的實際FIR低通數字濾波器的單位脈沖響應h(n)、幅頻衰減特性(20lg|H(ejw)|)、相頻特性及幅頻特性|H(ejw)|的波形，h(n)和|H(ejw)|圖中標出了所選窗函數類型及其長度N值。對四種窗函數（N=15和N=33）的程序運行結果如圖4-2到圖4-9所示，由圖可以看出用各種窗函數設計的FIR濾波器的阻帶最小衰減及過渡帶均與教材中一致。在通帶內均為嚴格相位特性。五. 實驗報告要求(1) 簡述實驗原理及目

26、的。(2) 按照實驗步驟以及要求，比較各種情況下的濾波性能，說明窗口長度N和窗函數類型對濾波特性的影響。(3) 總結用窗函數法設計FIR濾波器的主要特點。(4) 簡要回答思考題。六. 思考題(1) 如果給定通帶截止頻率和阻帶截止頻率以及阻帶最小衰減，如何用窗函數法設計線性相位低通濾波器，寫出設計步驟。(2) 如果要求用窗函數法設計帶通濾波器，而且給定上、下邊帶截止頻率為和，試求理想帶通的單位脈沖響應。附錄一 MATLAB信號處理工具箱函數MATLAB信號處理工具箱包含了許多信號處理函數。這些函數的格式和使用說明很容易通過HELP命令查到。調用這些函數很容易驗證數字信號處理習題答案，使數字濾波器

27、分析與設計的程序非常簡單。例如，已知系統函數H(z)，手工畫其幅頻響應曲線和相頻曲線是相當困難的。但只要調用函數freqz(b，a)就可以畫出來(b和a分別為H(z)的分子和分母多項式系數)。為了便于讀者分類快速查詢，下面按功能列出MATLAB信號處理工具箱函數。一. 表附1-1 波形產生函數名功能sawtooth產生鋸齒波或三角波square產生方波sinc產生sinc或者函數diric產生Dirichlet或者周期sinc函數二. 表附1-2 濾波器分析和實現函數名功能abs求絕對值（幅值）angle求相角conv求卷積fftfilt重疊相加法FFT濾波器實現filter直接濾波器實現fi

28、ltfilt零相位數字濾波filticfilter函數初始條件選擇freqs模擬濾波器頻率響應freqspace頻率響應中的頻率間隔freqz數字濾波器頻率響應grpdelay平均濾波器延遲（群延遲）impz數字濾波器的沖激響應zplane離散系統零極點圖三. 表附1-3 線性系統變換函數名功能convmtx卷積矩陣poly2rc從多項式系數中計算反射系數rc2poly從反射系數中計算多項式系數residuezZ變換部分分式展開或留數計算sos2ss變系統二階分割形式為狀態空間形式sos2tf變系統二階分割形式為傳遞函數形式sos2zp變系統二階分割形式為零極點增益形式ss2sos變系統狀態空

29、間形式為二階分割形式ss2tf變系統狀態空間形式為傳遞函數形式ss2zp變系統狀態空間形式為零極點增益形式tf2ss變系統傳遞函數形式為狀態空間形式tf2zp變系統傳遞函數形式為零極點增益形式zp2sos變系統零極點增益形式為二階分割形式zp2ss變系統零極點增益形式為狀態空間形式zp2tf變系統零極點增益形式為傳遞函數形式四. 表附1-4 IIR濾波器設計函數名功能besselfBessel（貝塞爾）模擬濾波器設計butterButterworth（比特沃思）濾波器設計cheby1Chebyshev （切比雪夫）I型濾波器設計cheby2Chebyshev （切比雪夫）II型濾波器設計ell

30、ip橢圓濾波器設計yulewalk遞歸數字濾波器設計五. 表附1-5 IIR濾波器階的選擇函數名功能buttordButterworth（比特沃思）濾波器階的選擇cheb1ordChebyshev （切比雪夫）I型濾波器階的選擇cheb2ordChebyshev （切比雪夫）II型濾波器階的選擇ellipord橢圓濾波器階的選擇六. 表附1-6 FIR濾波器設計函數名功能fir1基于窗函數的FIR濾波器設計標準響應fir2基于窗函數的FIR濾波器設計任意響應firls最小二乘FIR濾波器設計intfilt內插FIR濾波器設計remezParks-McCellan最優FIR濾波器設計remezo

31、rdParks-McCellan最優FIR濾波器階估計七. 表附1-7 變換函數名功能czt線性調頻Z變換dct離散余弦變換（DCT）idct逆離散余弦變換dftmtx離散傅立葉變換矩陣fft一維快速傅立葉變換ifft一維逆快速傅立葉變換fftshift重新排列FFT的輸出hilbertHilbert（希爾伯特）變換八. 表附1-8 統計信號處理函數名功能cov協方差矩陣xcov互協方差函數估計corrcoef相關系數矩陣xcorr互相關函數估計cohere相關函數平方幅值估計csd互譜密度（CSD）估計psd信號功率譜密度（PSD）估計tfe從輸入輸出中估計傳遞函數九. 表附1-9 窗函數函

32、數名功能boxcar矩形窗triang三角窗bartlettBartlett（巴特利特）窗hammingHamming（哈明）窗hanningHanning（漢寧）窗blackmanBlackman（布萊克曼）窗chebwinChebyshev （切比雪夫）窗kaiserKaiser（凱澤）窗十. 表附1-10 參數化建模函數名功能invfreqs模擬濾波器擬合頻率響應invfreqz離散濾波器擬合頻率響應prony利用Prony法的離散濾波器擬合時間響應stmcb利用Steiglitz-Mcbride迭代方法求線性模型levinsonLevinson-Durbin遞歸算法lpc線性預測系數十一. 表附1-11 特殊操作函數名功能rceps實倒譜和最小相位重構cceps倒譜分析和最小相位重構decimate降低序列的取樣速率interp提高取樣速率（內插）resample改變取樣速率medifiltl一維中值濾波deconv反卷積和多項式除法modulate通訊仿真中的調制domod通訊仿真中的解調vco電壓控制振蕩器specgram頻譜分析十二. 表附1-12 模擬原形濾波器設計函數名功能besselapBessel（貝塞爾）模