《數(shù)字信號處理》課程設計報告-數(shù)字濾波器的設計與信號的濾波處理

《數(shù)字信號處理》課程設計報告目錄TOC\o"1-4"\h\z\u一、設計題目 1二、設計目的 1三、課程設計的內容 1四、課程設計需要完成的工作 1五、設計過程與結果 15.1三組原始數(shù)據(jù)的時域圖像 15.2三組原始數(shù)據(jù)的頻譜圖 25.3數(shù)字濾波器的設計與信號濾波 25.3.1IIR數(shù)字濾波器設計 35.3.1.1.Butterworth低通濾波器 35.3.1.2．切比雪夫I型低通濾波器 35.3.1.3．切比雪夫II型低通濾波器 45.3.1.4．橢圓低通濾波器 65.3.2FIR（窗函數(shù)法）濾波器設計 65.3.2.1．Hamming窗函數(shù)濾波器 65.3.2.2．頻率采樣設計法、最優(yōu)等波紋設計法 75.4結果分析 7六、濾波器分析 8七、個人思考 107.1FIR濾波器的窗口設計法 107.2IIR數(shù)字濾波器設計 12八、拓展 148.1FDATool界面說明 158.2FDATool性能分析介紹 158.3FDATool濾波器系數(shù)導出 15九、收獲與感想 16十、MATLAB程序代碼 16《數(shù)字信號處理》課程設計報告——數(shù)字濾波器的設計與信號的濾波處理一、設計題目應用Matlab軟件對在實驗室采集的信號進行頻譜分析并設計數(shù)字濾波器對信號進行濾波處理。二、設計目的為了鞏固所學的數(shù)字信號處理理論知識，使學生對信號的采集、處理、讀取存儲等有一個系統(tǒng)的掌握和理解，學會數(shù)字濾波器的設計與應用。三、課程設計的內容采集個人的呼吸流量電壓信號、左右肺電壓信號，共三組信號，時長5分鐘，采樣率為7.812Hz，以獲得處理所需信號。畫出采樣信號的時域波形和頻譜分析圖，并設計合適的數(shù)字濾波器對采集的信號數(shù)據(jù)進行濾波，畫出濾波后信號的時域波形，并與濾波前信號圖像作對比，分析信號的變化。并通過對比幾種數(shù)字濾波器的濾波結果，分析它們各自的特點。四、課程設計需要完成的工作電壓信號的采集；信號的頻譜分析；數(shù)字濾波器的設計；對信號進行濾波處理；對濾波前后信號進行對比分析；對比幾種數(shù)字濾波器的濾波結果，分析它們各自的特點。五、設計過程與結果5.1三組原始數(shù)據(jù)的時域圖像由圖像可以看出，三組數(shù)據(jù)原始數(shù)據(jù)具有較多的噪聲（噪聲出現(xiàn)的原因將在后面的信號分析結果中進行說明），從圖像上難以對數(shù)據(jù)特點進行分析，因此需要進行合適的濾波后才能進行數(shù)據(jù)的分析。為了得出數(shù)//據(jù)的主頻率一遍進行濾波設計，下面將數(shù)據(jù)進行FFT，畫出頻譜圖進行分析。5.2三組原始數(shù)據(jù)的頻譜圖從三組數(shù)據(jù)進行FFT后的圖像中可以看出，三組數(shù)據(jù)的有效頻率都在0到0.5Hz之間，0.5Hz以上的可視為噪聲，應選擇低通濾波器對噪聲加以濾除。這個結果從生理學上來說是正確的，因為根據(jù)資料顯示，人呼吸的頻率約在0.2Hz到0.6Hz之間，因此對信號進行頻譜分析后會得出這樣的圖像結果。下面根據(jù)這個結果進行低通濾波器的設計。5.3數(shù)字濾波器的設計與信號濾波數(shù)字濾波器主要分為兩種：FIR(FiniteImpulseResponse)濾波器（又名：“有限長單位沖激響應濾波器”）以及IIR(InfiniteImpulseResponse)數(shù)字濾波器（又名“無限脈沖響應數(shù)字濾波器”，或“遞歸濾波器”）。IIR數(shù)字濾波器與FIR數(shù)字濾波器的區(qū)別主要有：1、單位響應IIR數(shù)字濾波器單位響應為無限脈沖序列，而FIR數(shù)字濾波器單位響應為有限的；FIR濾波器，也就是“非遞歸濾波器”，沒有引入反饋。這種濾波器的脈沖響應是有限的。2、幅頻特性IIR數(shù)字濾波器幅頻特性精度很高，不是線性相位的，可以應用于對相位信息不敏感的音頻信號上；FIR數(shù)字濾波器的幅頻特性精度較之于IIR數(shù)字濾波器低，但是線性相位，就是不同頻率分量的信號經(jīng)過FIR濾波器后他們的時間差不變，這是很好的性質。3、實時信號處理FIR數(shù)字濾波器是有限的單位響應也有利于對數(shù)字信號的處理，便于編程，用于計算的時延也小，這對實時的信號處理很重要。5.3.1IIR數(shù)字濾波器設計 5.3.1.1.Butterworth低通濾波器5.3.1.2．切比雪夫I型低通濾波器%通帶波紋Rp=0.025的圖像：%通帶波紋Rp較小效果較好，以下給出Rp=1的圖像：%從結果圖可以看出，Rp越大，濾波器在通帶處的波動越大。5.3.1.3．切比雪夫II型低通濾波器%阻帶衰減As=60的圖像： %阻帶衰減As較大效果較好，以下給出As=15的圖像：%從結果圖可以看出，As越小，濾波器在阻帶處的波動越大。5.3.1.4．橢圓低通濾波器5.3.2FIR（窗函數(shù)法）濾波器設計5.3.2.1．Hamming窗函數(shù)濾波器利用Matlab可以進行的窗函數(shù)濾波器類型還有很多，例如：（1）w=boxcar(M)在數(shù)組w中產(chǎn)生M點的矩形窗函數(shù)。（2）w=triang(M)在數(shù)組w中產(chǎn)生M點的Bartlett（三角形）窗函數(shù)。（3）w=hann(M)在數(shù)組w中產(chǎn)生M點的Hanning窗函數(shù)。（4）w=hamming(M)在數(shù)組w中產(chǎn)生M點的Hamming窗函數(shù)。（5）w=blackman(M)在數(shù)組w中產(chǎn)生M點的Blackman窗函數(shù)。（6）w=kaiser(M，beta)在數(shù)組w中產(chǎn)生M點的kaiser窗函數(shù)。（7）……這些窗函數(shù)需要根據(jù)不同的實際條件進行選擇，這里就不再一一贅述了。5.3.2.2．頻率采樣設計法、最優(yōu)等波紋設計法FIR濾波器的設計方法還有很多，例如頻率采樣設計法、最優(yōu)等波紋設計法等，這些方法在這里就不去設計了，查閱相關書籍資料便可得知詳情。5.4結果分析從信號原始圖上看，在310-820信號點之間出現(xiàn)了比較強烈噪聲干擾，特別是左肺電極電壓信號，結合下面的FFT分析，噪聲信號集中在2Hz左右，具體原因難以說明，據(jù)了解，這次試驗中采集的信號大多都存在著這樣一種現(xiàn)象，估計原因為以下幾點：心電信號的干擾：左肺信號測試電極位于心臟處，因此更容易受到心電信號的影響，產(chǎn)生噪聲干擾。呼吸速率的影響：由于第一次進行這個實驗的測試，全程用嘴巴呼吸，有很大程度上的緊張、呼吸不習慣等原因存在，特別卡礙事一段時間后，呼吸頻率會變得急促，后來經(jīng)過調整便能恢復，且呼吸頻率也會直接影響心跳速率，因此，信號在中段處的噪聲會特別明顯。測試過程中身體的移動可能導致噪聲的產(chǎn)生。從對信號進行FFT處理后得出的圖像上看，有效信號集中在0.1Hz-0.6Hz之間，這個結果從生理學上來說是正確的，因為根據(jù)資料顯示，人呼吸的頻率約在0.2Hz到0.6Hz之間，因此對信號進行頻譜分析后會得出這樣的圖像結果。在濾波器的設計過程中，濾波器設計的參數(shù)選擇很重要，不同的參數(shù)對濾波的效果產(chǎn)生很大的影響，因此在進行濾波器設計的過程中，必須花心思去考慮好設計參數(shù)，以便得到更好的濾波器。方才采用了4種不同的濾波器對信號進行了濾波處理，以下選取Butterworth濾波前后圖像進行對比放大：左肺信號：右肺信號：呼吸信號：從濾波前后的圖像對比可以看出，信號的噪聲已經(jīng)得到較好的濾除，圖像清晰，便于進行以后所需要的信號數(shù)據(jù)分析的工作。六、濾波器分析方才使用了IIR數(shù)字濾波器中三種最常用的濾波器類型，共4種濾波器，它們在濾波階數(shù)設置一樣的情況下對采集的信號做出濾波，現(xiàn)對比4種濾波器的濾波效果，分析它們各自的特點：四種濾波器濾波后的圖像：四種濾波器的幅頻特性曲線：由以上結果圖分析可知，濾波器取相同階數(shù)時：Butterworth低通濾波器具有單調下降的幅頻特性，但過渡帶最寬。兩種類型的Chebyshev濾波器的過渡帶寬度相等，比Butterworth濾波器的過渡帶窄，但比橢圓濾波器的過渡帶寬。ChebyshevI型濾波器在通帶具有等波紋幅頻特性，阻帶內單調；ChebyshevII型濾波器通帶內單調，二阻帶內等波紋。橢圓濾波器阻帶和通帶都是等波紋的，但過渡帶最窄。七、個人思考Matlab軟件是一款功能很強大的軟件，利用Matlab軟件可以方便地進行濾波器的設計。課程設計的主要目的在于鞏固所學的數(shù)字信號處理理論知識，因此如果僅僅是利用Matlab簡單地進行濾波器設計并處理數(shù)據(jù)，這難以達到課程設計的主要目的，因此，為了更好地去理解各種濾波器，也鍛煉自己的編寫程序的能力，下面通過書本上的理論知識加上Matlab函數(shù)對FIR、IIR中的集中常見濾波器進行設計，以加深對濾波器的理解。7.1FIR濾波器的窗口設計法窗函數(shù)名過渡帶寬最小阻帶衰減近似值準確值矩形4π1.8π21dBBartlett/巴特利特86.1π25dBHanning/漢寧86.2π44dBHamming/哈明86.6π53dBBlackman/布萊克曼1211π74dBMatlab提供了命令用于產(chǎn)生這些窗函數(shù)：（1）w=boxcar(M)在數(shù)組w中產(chǎn)生M點的矩形窗函數(shù)。（2）w=triang(M)在數(shù)組w中產(chǎn)生M點的Bartlett（三角形）窗函數(shù)。（3）w=hann(M)在數(shù)組w中產(chǎn)生M點的Hanning窗函數(shù)。（4）w=hamming(M)在數(shù)組w中產(chǎn)生M點的Hamming窗函數(shù)。（5）w=blackman(M)在數(shù)組w中產(chǎn)生M點的Blackman窗函數(shù)。（6）w=kaiser(M，beta)在數(shù)組w中產(chǎn)生M點的kaiser窗函數(shù)。利用這些指令便可以用Matlab設計基于窗函數(shù)法的FIR濾波器。為了方便使用，在這里給出一個理想低通脈沖響應hd(n)，編寫程序如下：functionhd=ideal_lp(wc,M);%wc為截止頻率，M為窗函數(shù)長度alpha=(M-1)/2;n=[0:1:(M-1)];m=n-alpha;fc=wc/pi;hd=fc*sinc(fc*m);end例子：設計一個數(shù)字FIR低通濾波器，器技術指標如下：wp=0.2π，ws=0.3π，Rp=0.25dB，As=50dB題解: 從表格可知，用Hamming和Blackman窗函數(shù)都能提供50dB的阻帶衰減，現(xiàn)選用Hamming窗函數(shù)，因為它能給出比較小的過渡帶，因此有較低的階。盡管在Hamming窗函數(shù)濾波器的設計中沒有用到通帶波紋值Rp，但是還是必須從設計中校核真正的波紋值，以證實波紋確實符合題目要求。設計結果如下：Matlab代碼運行結果：>>M=67>>Rp=0.0394>>As=52結果顯示真實的阻帶衰減是52dB，通帶波紋是0.0394dB，顯然滿足題目設計要求。我們知道，Kaiser窗函數(shù)對同一個M值可以給出不同的過渡帶寬，它的這種性能使其可以和Hamming窗相匹敵，下面根據(jù)上述指標，選擇Kaiser窗設計這個低通濾波器，設計結果如下：運行結果：>>M=61>>beta=4.5513>>Rp=0.0442>>As=527.2IIR數(shù)字濾波器設計IIR數(shù)字濾波器具有無限長脈沖響應，因此能偶與模擬濾波器相匹敵；一般來說，所有的模擬濾波器都有無限長脈沖響應。因此，IIR數(shù)字濾波器設計的基本方法是利用復值映射將模擬濾波器變換為數(shù)字濾波器。復值映射的變換方法有很多種，如脈沖響應不變變換、雙線性變換、匹配Z變換等，在實際中采用的最為普遍的方法是雙線性變換法。在Matlab中，有buttord、cheb2ord、cheb2ord和ellipord函數(shù)可以在已知設計指標時給出最小濾波階數(shù)N和截止頻率wn；以下程序嘗試采用定義中給出的公式編程確定這些參數(shù)，以加深對雙線性變換的理解。下面利用雙線性變換法對三種最廣泛使用的濾波器進行設計。例子：設計一個數(shù)字FIR低通濾波器，器技術指標如下：wp=0.2π，ws=0.3π，Rp=1dB，As=15dB（1）巴特沃斯低通濾波器：設計結果：>>ButterworthFilterOrder=6>>直接型>>b=0.00060.00350.00870.01160.00870.00350.0006>>a=1.0000-3.31434.9501-4.14332.0275-0.54580.0628>>級聯(lián)型>>BA= 1.00002.00001.00001.0000-0.94590.2342 1.00002.00481.00481.0000-1.05410.3753 1.00001.99520.99521.0000-1.31430.7149>>b0= 5.7969e-04則這個濾波器系統(tǒng)函數(shù)是：H（2）數(shù)字切比雪夫I型、II型低通濾波器：運行結果： %切比雪夫I型>>ChebyshevIFilterOrder=4>>直接型>>b=0.00180.00730.01100.00730.0018>>a=1.0000-3.05433.8290-2.29250.5507>>級聯(lián)型>>BA= 1.00002.00001.00001.0000-1.55480.6493 1.00002.00001.00001.0000-1.49960.8482>>b0=0.0018則這個濾波器系統(tǒng)函數(shù)是：H%切比雪夫II型>>ChebyshevIIFilterOrder=4>>直接型>>b=0.1797-0.09160.2525-0.09160.1797>>a=1.0000-1.55081.3423-0.47070.1079>>級聯(lián)型>>BA= 1.00000.55741.00001.0000-0.41830.1503 1.0000-1.06711.00001.0000-1.13250.7183>>b0=0.1797則這個濾波器系統(tǒng)函數(shù)是：H（3）數(shù)字橢圓低通濾波器：運行結果：>>EllipticFilterOrder=3>>直接型>>b=0.1214-0.0511-0.05110.1214>>a=1.0000-2.11121.7843-0.5325>>級聯(lián)型>>BA= 1.00001.000001.0000-0.618301.0000-1.42111.00001.0000-1.49280.8612>>b0= 0.1214則這個濾波器系統(tǒng)函數(shù)是：H三類濾波器的比較：原型階次N最小阻帶衰減/dB巴特沃斯615切比雪夫425橢圓327在給定了的相同的濾波器設計參數(shù)后，從階數(shù)上看，橢圓濾波器給出了最低階數(shù)的設計而巴特沃斯濾波器給出了三者中最高的階數(shù)設計，但是我們也知道，橢圓濾波器在通帶內非線性相位最為嚴重，因此需要根據(jù)實際情況取舍，選擇合適的濾波器設計。八、拓展除了采用編程的方法實現(xiàn)濾波器之外，Matlab中自帶工具箱FDATool也能很方便快速的實現(xiàn)濾波器的設計。另外FPGA、DSP等實現(xiàn)數(shù)字濾波器算法時，經(jīng)常要用到濾波器系數(shù)，采用FDATool工具箱也能快速的得到濾波器系數(shù)。方法1：在Matlab主界面輸入FDATool，回車即可調用FDATool工具箱。方法2：單擊Matlab主窗口下方的“Start”按鈕，如圖B.1所示，選擇菜單“ToolBox”

→“Filter

Design”

→“Filter

Design

&

Analysis

Tool（FDATool）”命令，打開FDATool。其界面如下所示：8.1FDATool界面說明FDATool（FilterDesign&AnalysisTool）是Matlab信號處理工具箱里專用的濾波器設計分析工具，Matlab6.0以上的版本還專門增加了濾波器設計工具箱（FilterDesignToolbox）。FDATool可以設計幾乎所有的基本的常規(guī)濾波器，包括fir和iir的各種設計方法。它操作簡單，方便靈活。FDATool界面總共分兩大部分，一部分是DesignFilter，在界面的下半部，用來設置濾波器的設計參數(shù)，另一部分則是特性區(qū)，在界面的上半部分，用來顯示濾波器的各種特性。DesignFilter部分主要分為：FilterType（濾波器類型）選項，包括lowpass（低通）、highpass（高通）、bandpass（帶通）、bandstop（帶阻）和特殊的fir濾波器。DesignMethod（設計方法）選項，包括IIR數(shù)字濾波器的butterworth（巴特沃思）法、chebyshevtypei（切比雪夫i型）法、chebyshevtypeii（切比雪夫ii型）法、elliptic（橢圓濾波器）法和fir濾波器的equiripple法、least-squares（最小乘方）法、window（窗函數(shù)）法。FilterOrder（濾波器階數(shù)）選項，定義濾波器的階數(shù)，包括specifyorder（指定階數(shù)）和minimumorder（最小階數(shù)）。在specifyorder中填入所要設計的濾波器的階數(shù)（n階濾波器，specifyorder＝n-1），如果選擇minimumorder則Matlab根據(jù)所選擇的濾波器類型自動使用最小階數(shù)。FrenquencySpecifications選項，可以詳細定義頻帶的各參數(shù)，包括采樣頻率fs和頻帶的截止頻率。它的具體選項由filtertype選項和designmethod選項決定，例如bandpass（帶通）濾波器需要定義fstop1（下阻帶截止頻率）、fpass1（通帶下限截止頻率）、fpass2（通帶上限截止頻率）、fstop2（上阻帶截止頻率），而lowpass（低通）濾波器只需要定義fstop1、fpass1。采用窗函數(shù)設計濾波器時，由于過渡帶是由窗函數(shù)的類型和階數(shù)所決定的，所以只需要定義通帶截止頻率，而不必定義阻帶參數(shù)。MagnitudeSpecifications選項，可以定義幅值衰減的情況。例如設計帶通濾波器時，可以定義wstop1（頻率fstop1處的幅值衰減）、wpass（通帶范圍內的幅值衰減）、wstop2（頻率fstop2處的幅值衰減）。當采用窗函數(shù)設計時，通帶截止頻率處的幅值衰減固定為6db，所以不必定義。WindowSpecifications選項，當選取采用窗函數(shù)設計時，該選項可定義，它包含了各種窗函數(shù)。8.2FDATool性能分析介紹在設計完所需的濾波器后，往往需要對設計好的FIR濾波器進行相關的性能分析，以便了解該濾波器是否滿足設計要求。分析操作步驟如下：（1）選擇FDATool的菜單“Analysis”→“Magnitude

Response”，啟動幅頻響應，x軸為頻率，y軸為幅度值（單位為dB）。（2）選擇菜單“Analysis”→“Phase

Response”，啟動相頻響應分析；（3）選擇菜單“Analysis”→“Magnitude

and

Phase

Response”來啟動濾波器幅頻特性與相頻特性的比較分析；（4）選擇菜單“Analysis”→“Group

Delay

Response”，啟動群時延分析。FDATool還提供了以下幾種分析工具：群時延響應分析（GroupDelayResponse）。沖激響應分析（Impulse

Response）；階躍響應分析（Step

Response）；零極點圖分析（Pole/Zero

Plot）。8.3FDATool濾波器系數(shù)導出為導出設計好的濾波器系數(shù)，選擇FDATool菜單的“File”→“Export”命令，打開Export（導出）對話框，如下圖所示。在該窗口中，選擇導出到工作區(qū)（Workplace）。這時濾波器系數(shù)就存入到一個一維變量Num中了。這時Num中的元素是以小數(shù)形式出現(xiàn)的。九、收獲與感想濾波器的設計在實際生活中必不可少，當然也比課程設計的題目要復雜得多，但是它們設計的成功與否，歸結到底還是參數(shù)的設置以及濾波器類型的選擇。通過前面的濾波器設計，我們已經(jīng)知道每一種濾波器都有各自的優(yōu)缺點，在實際情況下到底該選取哪一種濾波器，主要取決于每種濾波器的優(yōu)點在設計中的重要性。因此，濾波器的設計必須做到以下幾點：確定所需的濾波器的用途，并對其進行性能需求分析，以明確該系統(tǒng)所需的各種性能指標；必須有系統(tǒng)的設計思維，把每一個細節(jié)都放到整個系統(tǒng)中考慮，考慮整個系統(tǒng)的可行性、完整性、穩(wěn)定性和功能的實現(xiàn)；經(jīng)過這次的課程設計，我對《數(shù)字信號處理》這門課程有了更進一步的認識，對所學的知識也有了更深一層的理解和掌握，不再僅僅是停留在“學過”這個層面上，還學會了如何在實際應用中運用自己所學到的知識?；蛟S在平時的課堂上，學習理論知識的時候，我們會覺得“這種東西”有什么用啊，因而沒有過多地去在意它，當進行這個課程設計的時候，次貿然發(fā)現(xiàn)其實那些理論是多么的重要，而且在實際運用中并不像理論課上那樣的枯燥乏味，更多的是需要靈活地去運用所學的理論解決實際問題，不再像寫作業(yè)那樣是幾乎的過程那樣幾乎是定勢的一種計算。在設計的過程中遇到過不少的困難，然后需要查閱相關的資料并學習相關的知識一解決問題，這是個考驗耐心的過程，在這個過程中必須保持良好的心態(tài)，切記心浮氣躁，不然問題只會無法有效地解決，造成惡性循環(huán)。面對一個問題，或許思索良久都沒有好的解決辦法，這是不妨放開思維，大膽地嘗試，說不定通過另一種思路問題就迎刃而解了，有時候只要稍微轉一下彎，不要鉆牛角尖，結果就會不一樣了?？偠灾?，通過這次的課程設計，不僅讓我鞏固了課堂上學習到的專業(yè)知識，也鍛煉了我各方面的能力，讓我收獲了很多。十、MATLAB程序代碼5.1畫出三組原始數(shù)據(jù)的時域圖像：Matlab代碼：clearA=xlsread('data.xls');a=A(:,1)';b=A(:,2)';c=A(:,3)';subplot(3,1,1);plot(a);xlim([100,1000]);title('OriginalData_a');subplot(3,1,2);plot(b);xlim([100,1000]);title('OriginalData_b');subplot(3,1,3);plot(c);xlim([100,1000]);title('OriginalData_c');5.2畫出三組原始數(shù)據(jù)的頻譜圖：Matlab代碼：fs=7.812;N=1000;%采樣頻率和數(shù)據(jù)點數(shù)n=0:N-1;t=n/fs;%時間序列ya=fft(a,N);%對信號進行快速Fourier變換mag1=abs(ya);%求得快速Fourier變換后的振幅yb=fft(b,N);mag2=abs(yb);yc=fft(c,N);mag3=abs(yc);f=n*fs/N;%頻率序列subplot(3,1,1);plot(f,mag1);%繪出隨頻率變化的振幅xlim([0,7.812/2]);ylim([0,10]);title('FFTofdata_a');subplot(3,1,2);plot(f,mag2);xlim([0,7.812/2]);ylim([0,10]);title('FFTofdata_b');subplot(3,1,3);plot(f,mag3);xlim([0,7.812/2]);ylim([0,10]);title('FFTofdata_c');xlabel('frequency/Hz');濾波器的設計與信號濾波：5.3.1IIR數(shù)字濾波器設計： 5.3.1.1．巴特沃斯低通濾波器：Matlab代碼：fs=7.812;Wn=2*0.5/fs;%歸一化N=20;%設置階數(shù)[B,A]=butter(N,Wn,'low');X=filter(B,A,a);%a可以換成b、c以對b、c組數(shù)據(jù)濾波[H,w]=freqz(B,A,512,'whole');mag=abs(H);db=20*log10((mag+eps)/max(mag));pha=angle(H);subplot(2,2,1);plot(X);xlim([100,1000]);title('用butterworth濾波后圖像');subplot(2,2,2);plot(w/pi,mag);xlim([0,1]);title('幅頻特性曲線');xlabel('單位為π');subplot(2,2,3);plot(w/pi,pha);xlim([0,1]);title('相頻特性曲線');xlabel('單位為π');subplot(2,2,4);plot(w/pi,db);xlim([0,1]);title('幅頻特性曲線(dB)');xlabel('單位為π');5.3.1.2．切比雪夫I型低通濾波器：Matlab代碼：fs=7.812;wn=2*0.5/fs;N=10;%歸一化，并設定階數(shù)Rp=0.025;%通帶波紋[B,A]=cheby1(N,Rp,wn);X=filter(B,A,a);[H,w]=freqz(B,A,512,'whole');mag=abs(H);db=20*log10((mag+eps)/max(mag));pha=angle(H);subplot(2,2,1);plot(X);xlim([100,1000]);title('用切比雪夫I型濾波后圖像');subplot(2,2,2);plot(w/pi,mag);xlim([0,1]);title('幅頻特性曲線');xlabel('單位為π');subplot(2,2,3);plot(w/pi,pha);xlim([0,1]);title('相頻特性曲線');xlabel('單位為π');subplot(2,2,4);plot(w/pi,db);xlim([0,1]);title('幅頻特性曲線(dB)');xlabel('單位為π');5.3.1.3．切比雪夫II型低通濾波器：Matlab代碼：fs=7.812;wn=2*0.5/fs;N=10;As=60;%阻帶衰減[B,A]=cheby2(N,As,wn);X=filter(B,A,a);[H,w]=freqz(B,A,512,'whole');mag=abs(H);db=20*log10((mag+eps)/max(mag));pha=angle(H);subplot(2,2,1);plot(X);xlim([100,1000]);title('用切比雪夫II型濾波后圖像');subplot(2,2,2);plot(w/pi,mag);xlim([0,1]);title('幅頻特性曲線');xlabel('單位為π');subplot(2,2,3);plot(w/pi,pha);xlim([0,1]);title('相頻特性曲線');xlabel('單位為π');subplot(2,2,4);plot(w/pi,db);xlim([0,1]);title('幅頻特性曲線(dB)');xlabel('單位為π');5.3.1.4．橢圓低通濾波器：Matlab代碼：fs=7.812;wn=2*0.5/fs;Rp=0.025;As=60;N=5;%設定設計參數(shù)[B,A]=ellip(N,Rp,As,wn);X=filter(B,A,a);[H,w]=freqz(B,A,512,'whole');mag=abs(H);db=20*log10((mag+eps)/max(mag));pha=angle(H);subplot(2,2,1);plot(X);xlim([100,1000]);title('用ellip濾波后圖像');subplot(2,2,2);plot(w/pi,mag);xlim([0,1]);title('幅頻特性曲線');xlabel('單位為π');subplot(2,2,3);plot(w/pi,pha);xlim([0,1]);title('相頻特性曲線');xlabel('單位為π');subplot(2,2,4);plot(w/pi,db);xlim([0,1]);title('幅頻特性曲線(dB)');xlabel('單位為π');5.3.2．FIR（窗函數(shù)法）濾波器設計：5.3.2.1．Hamming窗函數(shù)濾波器：Matlab代碼：wp=0.4;ws=0.6;Rp=0.2;As=60;fs=7.812;delta=ws-wp;N=ceil(6.6*pi/delta);wn=(wp+ws)/2;[B,A]=fir1(N,wn/pi,hamming(N+1));%選擇窗函數(shù)X=filter(B,A,a);[H,w]=freqz(B,A,512,'whole');mag=abs(H);db=20*log10((mag+eps)/max(mag));pha=angle(H);subplot(2,2,1);plot(X);xlim([100,1000]);title('用hamming窗函數(shù)濾波后圖像');subplot(2,2,2);plot(w/pi,mag);xlim([0,1]);title('幅頻特性曲線');xlabel('單位為π');subplot(2,2,3);plot(w/pi,pha);xlim([0,1]);title('相頻特性曲線');xlabel('單位為π');subplot(2,2,4);plot(w/pi,db);xlim([0,1]);title('幅頻特性曲線(dB)');xlabel('單位為π');7.1Hamming窗函數(shù)濾波器設計：Matlab代碼：clearwp=0.2*pi;ws=0.3*pi;tr_width=ws-wp;M=ceil(6.6*pi/tr_width)+1n=[0:1:(M-1)];wc=(ws+wp)/2;hd=ideal_lp(wc,M);w_ham=(hamming(M))';h=hd.*w_ham;[H,w]=freqz(h,[1],1000,'whole');mag=abs(H);db=20*log10((mag+eps)/max(mag));delta_w=2*pi/1000;Rp=-(min(db(1:1:(wp/delta_w+1))))As=-round(max(db(ws/delta_w+1:1:501)))%%%%%%%%%%plots%%%%%%%%%%subplot(2,2,1);stem(n,hd);title('IdealImpulseResponse');axis([0M-1-0.10.3]);xlabel('n');ylabel('h_d(n)');subplot(2,2,2);stem(n,w_ham);title('HammingWindow');axis([0M-101.1]);xlabel('n');ylabel('w(n)');subplot(2,2,3);stem(n,h);title('ActualimpulseResponse');axis([0M-1-0.10.3]);xlabel('n');ylabel('h(n)');subplot(2,2,4);plot(w/pi,db,w/pi,-50,'-k');title('MagnitudeResponseindB');axis([01-12010]);xlabel('frequencyinpiunits');ylabel('dB');7.1Kaiser窗函數(shù)濾波器設計：Matlab代碼：clearwp=0.2*pi;ws=0.3*pi;tr_width=ws-wp;As=50;M=ceil((As-7.95)/(2.285*tr_width)+1)+1n=[0:1:(M-1)];beta=0.1102*(As-8.7)wc=(ws+wp)/2;hd=ideal_lp(wc,M);w_kai=(kaiser(M,beta))';h=hd.*w_kai;[H,w]=freqz(h,[1],1000,'whole');mag=abs(H);db=20*log10((mag+eps)/max(mag));delta_w=2*pi/1000;Rp=-(min(db(1:1:(wp/delta_w+1))))As=-round(max(db(ws/delta_w+1:1:501)))%%%%%%%%%%plots%%%%%%%%%%subplot(2,2,1);stem(n,hd);title('IdealImpulseResponse');axis([0M-1-0.10.3]);xlabel('n');ylabel('h_d(n)');subplot(2,2,2);stem(n,w_kai);title('kaiserWindow');axis([0M-101.1]);xlabel('n');ylabel('w(n)');subplot(2,2,3);stem(n,h);title('ActualimpulseResponse');axis([0M-1-0.10.3]);xlabel('n');ylabel('h(n)');subplot(2,2,4);plot(w/pi,db,w/pi,-50,'-k');title('MagnitudeResponseindB');axis([01-12010]);xlabel('frequency